#include <stdio.h>
#include <string.h>
#include <math.h>

#define  BAD  -1.0e30
#define  PI    3.1415
#define  SW9   5.4
#define  SW10  6.0
#define  SW12  7.2
#define  SW14  8.3


main()
 {
 FILE *fp,*fplog,*fpp,*fpin,*fpE,*fpout,*fprough,*fppl;
 int i,is,iti,i0,i1,i2,isuf,iav,iave,ipl,Dib,j,j0,j1,j2,ii,jj,k,ki,kf,kk,m,o,ch,chi,chf,len,exptno;
 int ictd1,ictd2;
 int Ndata[340],Tdatalog[340],Tdata[340],Mdata,Nxdata,Nydata,Nexpt,Npl,Nrain,Nrainsub,Ncntr,Nlz,lzpl[10];
 int drp,drpctd,drptot,drop[340],cast[340],ctddrop[340],ctdcast,dropi,dropf,Ndrp[50];
 int idate[50],iTdata[50],iend[50],badctd[340][50];
 int itopo,cntr,clr;
 int ilonmin,ilonmax,ilon,Nlat;
 int smooth,Di,choice,choici,choicf,plot,output,fit,wt;
 int dic,sgnx,sgny,lab,logaxis,test;
 int Sdata[3],SdataH[3];
 int badv,badT[50];
 int iE[340],istand[340],node[340],iCW[340],iCCW[340],mode[340];
 int surfrefl[340],bottrefl[340],surfinten[340],bottintens[340];
 char str[72],str1[72],str2[72],str3[72],str4[72],expt[30],font[20],flag[6],Ndat[9],dir[50][50],ocean[50][50];
 char ifile[120],ofile[120],pfile[120],sfile[120],topodir[75];
 char ctddir[50][50],ctdpre[50][50],ctdsuf[50][50];
 
 float yr,mo,da,hr,mn,lat0,lon0,latd,latm,lond,lonm;
 float lat[340],lon[340],zb[340],startda[340],bottomda[340],endda[340];
 float dOcart,dO1,Orad,sinO,cosO,sinlat,coslat,latkm,lonkm,f;
 float zmn[340],zmx,zwkbmn[340],zwkbmx;
 float HKEwkb[340],Emn,Emx,roughmn,roughmx;

 float latb[505],lonb[1205],zbb[1205][505],zbmn,zbmx;

 float Dzin[50],Dz[3],delz[3],kz[1025],lz[1025];
 float dEstat,Estat[5],dtopostat,topostat[8];
 float nE[5],nstand[5],nt[8],ndown[5],nup[5],nmode[5],nnode[5];
 float nstandE[5][5],ndownE[5][5],nupE[5][5],nmodeE[5][5],nnodeE[5][5];
 float nEm[5][5],nstandm[5][5],ndownm[5][5],nupm[5][5],nnodem[5][5];
 float nstandt[8][8],ndownt[8][8],nupt[8][8],nmodet[8][8],nnodet[8][8];
 float dlatbin,latbin[20],istandlat[20],idownlat[20],iuplat[20],itotlat[20];
 float imode1lat[20],imode2lat[20],imode3lat[20],imode4pluslat[20],imodetotlat[20];

 float lz0,lz1,kz0,kz1,Rkz,lhc,ufit,vfit;
 float zi,zf,zf1,hi,hf,zmid,zctdi,zctdf,sigi,sigf,signi,signf;
 float uav,vav,Zav,Zzav,ZzTav;
 float varu,varv,varZ,varZz,varZzT;
 float varSu,varSv,varSZ,varSZz,varSZzT;
 float varSSu,varSSv,varSSZ;
 float varSSuz,varSSvz,varSSZz,varSSZzT;
 float wndo,frac,Rz,junk,sini,sinf;
 float zbb1[340],varzbb5[340],varzbb15[340];
 
 float latmin,latmax,dlat,dsublat;
 float lonmin,lonmax,dlon,dsublon,dellon;
 float xmin,xmax,dx,dsubx,xmid,delx;
 float ymin,ymax,dy,dsuby,ymid,dely;
 float xplmin,xplmax,delxpl,xplmid;
 float yplmin,yplmax,delypl,yplmid;
 float cplmin,cplmax,delcpl,cplmid;
 float dcntr,zcntr;
 float dsubcntr,zsubcntr;
 float zcntrmin,zcntrmax,cntrmin,cntrmax;
 double xcntr[6],ycntr[6],fcntr[6];
 float r1,r2,r3,r4,xx,XX,XXo,XXc,yy,cc,rpl,xpl1,xpl2,ypl1,ypl2;
 double root,xax;
 double xpl[1600],ypl[1600],Wt[1600];
 float gray,red,green,blue,rainbo[350][4],lw,brown[4];

 badv=0;
 roughmn=0.;
 roughmx=3.;
 /* bad = -1e30; */
 brown[1]=0.4;   brown[2]=0.259;  brown[3]=0.040;
 /**********************************************parameter file ***/
 fpp=fopen("rotlat.p","r");
 fscanf(fpp,"%s",expt);
 printf("%s \n",expt);
 fscanf(fpp,"%s %s",str,font);
 fscanf(fpp,"%s %d %d",str,&smooth,&Di);
 fscanf(fpp,"%s %d",str,&plot);
 fscanf(fpp,"%s %f %f",str,&lat0,&lon0);
 latkm=PI*6370/180.;
 sinlat=sin(PI*lat0/180.);
 coslat=cos(Orad);
 f=2*PI*sinlat/12./3600.;
 lonkm=latkm*coslat;
 fscanf(fpp,"%s %f %f %f %f",str,&latmin,&latmax,&dlat,&dsublat);
 fscanf(fpp,"%s %f %f %f %f",str,&lonmin,&lonmax,&dlon,&dsublon);
 fscanf(fpp,"%s %d %d",str,&ki,&kf);
 fscanf(fpp,"%s",str);
 for (k=1;k<=kf;k=k+1) 
    fscanf(fpp,"%d %s %d %s %s %s %d %d %d %f %s",
           &j,dir[k],&Ndrp[k],ctddir[k],ctdpre[k],ctdsuf[k],&idate[k],&iTdata[k],&iend[k],&Dzin[k],ocean[k]);
 fclose(fpp);
 /*****************************************************************setup stat***/
 dEstat=dtopostat=1.;
 for (i=0;i<=4;i=i+1)
  {
  Estat[i]=dEstat*i;
  nE[i]=nstand[i]=nt[i]=ndown[i]=nup[i]=nmode[i]=nnode[i]=0.;
  for (j=0;j<=4;j=j+1)
     nstandE[i][j]=ndownE[i][j]=nupE[i][j]=nmodeE[i][j]=nnodeE[i][j]=0.;
  } /* next i */
 for (i=1;i<=3;i=i+1)
  {
  for (j=0;j<=4;j=j+1)
     nEm[i][j]=nstandm[i][j]=ndownm[i][j]=nupm[i][j]=nnodem[i][j]=0.;
  } /* next i */
 for (i=0;i<=7;i=i+1)
  {
  topostat[i]=dtopostat*i;
  nt[i]=0.;
  for (j=0;j<=7;j=j+1)
     nstandt[i][j]=ndownt[i][j]=nupt[i][j]=nmodet[i][j]=nnodet[i][j]=0.;
  } /* next i */
 /***************************************************setup lat bins***/
 Nlat=15;
 dlatbin=10.;
 for (o=1;o<=Nlat;o=o+1)
  {
  latbin[o]=-70.+(o-1)*dlatbin;
  istandlat[o]=idownlat[o]=iuplat[o]=itotlat[o]=0.;
  imode1lat[o]=imode2lat[o]=imode3lat[o]=imode4pluslat[o]=0.;
  } /* next o */
 /**********************************************INPUT LOOP***/
  /***********************************************************color***/
    sprintf(ifile,"%s","/home/kunze/plotlib/color/blackblueyellowred.dat");
    fp=fopen(ifile,"r");
    fscanf(fp,"%d",&Nrain);
    Nrainsub=Nrain;
    for (o=1;o<=Nrain;o=o+1)
       fscanf(fp,"%f %f %f",&rainbo[o][1],&rainbo[o][2],&rainbo[o][3]);
    fclose(fp);
    /************************************************************************section loop***/
        Emn=roughmn=100000000000000.;
	Emx=roughmx=-Emn;
        for (k=ki;k<=kf;k=k+1)
         {
         /************************************************************************input***/
         sprintf(ifile,"%s%s%s","/data/kunze/oldarchive/archive/ladcp/data/ladcp/",dir[k],"/roughlog");
         fprough=fopen(ifile,"r");
         fscanf(fprough,"%d",&Ndrp[0]);

         sprintf(ifile,"%s%s%s","/data/kunze/oldarchive/archive/ladcp/data/ladcp/",dir[k],"/rotlog");
         fpin=fopen(ifile,"r");
         sprintf(ifile,"%s%s%s","/data/kunze/oldarchive/archive/ladcp/data/ladcp/",dir[k],"/Elog");
         fpE=fopen(ifile,"r");

         fscanf(fpin,"%d",&Ndrp[k]);
	 /*
	 printf("%s %d \n",ifile,Ndrp[k]);
	 */
         for (drp=1;drp<=Ndrp[k];drp=drp+1)
          {
          fscanf(fpin,"%d %d %f %f %f %f %f %d %d %d %d %d %d %d %d %d %d",
                       &drop[drp],&cast[drp],&lat[drp],&lon[drp],&zmn[drp],&zb[drp],&zwkbmn[drp],
                       &iE[drp],&istand[drp],&node[drp],&iCW[drp],&iCCW[drp],&mode[drp],&surfrefl[drp],
	               &bottrefl[drp],&surfinten[drp],&bottintens[drp]);
	  /*
	  if (ch == 8 && k == 5) 
	    printf("%d%s%d %d %d %d %d \n",drp,"/",Ndrp[k],surfrefl[drp],bottrefl[drp],surfinten[drp],bottintens[drp]);
	  */
	  nE[0]=nE[0]+1.;
          /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * input E * * */
          fscanf(fpE,"%d %d %f %f %f %f %f %f",
                       &drop[0],&cast[0],&lat[0],&lon[0],&zmn[0],&zb[0],&zwkbmn[0],
                       &HKEwkb[drp]); 
	  if (HKEwkb[drp] < Emn) Emn=HKEwkb[drp];
	  if (HKEwkb[drp] > Emx) Emx=HKEwkb[drp];
          /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * input topo roughness * * */
          fscanf(fprough,"%d %f %f %f %f %f %f %f",
			  &j,&lat[0],&lon[0],&zb[0],&zbb1[drp],&xx,&varzbb5[drp],&varzbb15[drp]);
	  if (varzbb5[drp] < roughmn) roughmn=varzbb5[drp];
	  if (varzbb5[drp] > roughmx) roughmx=varzbb5[drp];
          while (lon[drp] < 0.) lon[drp]=lon[drp]+360.;
          while (lon[drp] > 360.) lon[drp]=lon[drp]-360.;
	  j=lat[drp]/sqrt(lat[drp]*lat[drp]);
	  if (j < 0)
	   {
	   o=iCW[drp];
	   iCW[drp]=iCCW[drp];
	   iCCW[drp]=o;
	   } /* end if (j */
	  /************************************************lat binning***/
	  for (o=1;o<=Nlat;o=o+1)
	   {
	   if (lat[drp] >= latbin[o] && lat[drp] <= latbin[o]+dlatbin)
	    {
	    if (istand[drp] == 1) istandlat[o]=istandlat[o]+1.;
	    if (iCW[drp] == 1) idownlat[o]=idownlat[o]+1.;
	    if (iCCW[drp] == 1) iuplat[o]=iuplat[o]+1.;

	    if (mode[drp] == 1) imode1lat[o]=imode1lat[o]+1.;
	    if (mode[drp] == 2) imode2lat[o]=imode2lat[o]+1.;
	    if (mode[drp] == 3) imode3lat[o]=imode3lat[o]+1.;
	    if (mode[drp] >= 4) imode4pluslat[o]=imode4pluslat[o]+1.;
	    } /* end if (lat */
	   } /* next o */
	  /************************************************pdf binning***/
	  /*
	if (lat[drp] > 10. || lat[drp] < -10.)
	 {	  
	 */
	  if (istand[drp] == 0) nstand[0]=nstand[0]+1.;
	   if (istand[drp] == 1) nstand[1]=nstand[1]+1.;
	  if (iCW[drp] == 0) ndown[0]=ndown[0]+1.;
	   if (iCW[drp] == 1) ndown[1]=ndown[1]+1.;
	  if (iCCW[drp] == 0) nup[0]=nup[0]+1.;
	   if (iCCW[drp] == 1) nup[1]=nup[1]+1.;
	  if (mode[drp] == 1) nmode[1]=nmode[1]+1.;
	   if (mode[drp] == 2) nmode[2]=nmode[2]+1.;
	    if (mode[drp] >= 3) nmode[3]=nmode[3]+1.;
	  if (node[drp] == -1) nnode[1]=nnode[1]+1.;
	   if (node[drp] == 0) nnode[2]=nnode[2]+1.;
	    if (node[drp] >= 1) nnode[3]=nnode[3]+1.;
	  /************************************************E binning***/
	  for (i=0;i<=3;i=i+1)
	   {
	   cc=log10(HKEwkb[drp]);
	   if (cc >= Estat[i] && cc <= Estat[i]+dEstat) nE[i]=nE[i]+1.;

	   if (istand[drp] == 0 && cc >= Estat[i] && cc <= Estat[i]+dEstat) nstandE[i][0]=nstandE[i][0]+1.;
 	    if (istand[drp] == 1 && cc >= Estat[i] && cc <= Estat[i]+dEstat) nstandE[i][1]=nstandE[i][1]+1.;
	   if (iCW[drp] == 0 && cc >= Estat[i] && cc <= Estat[i]+dEstat) ndownE[i][0]=ndownE[i][0]+1.;
	    if (iCW[drp] == 1 && cc >= Estat[i] && cc <= Estat[i]+dEstat) ndownE[i][1]=ndownE[i][1]+1.;
	   if (iCCW[drp] == 0 && cc >= Estat[i] && cc <= Estat[i]+dEstat) nupE[i][0]=nupE[i][0]+1.;
	    if (iCCW[drp] == 1 && cc >= Estat[i] && cc <= Estat[i]+dEstat) nupE[i][1]=nupE[i][1]+1.;
	   if (mode[drp] == 1 && cc >= Estat[i] && cc <= Estat[i]+dEstat) nmodeE[i][1]=nmodeE[i][1]+1.;
	    if (mode[drp] == 2 && cc >= Estat[i] && cc <= Estat[i]+dEstat) nmodeE[i][2]=nmodeE[i][2]+1.;
	     if (mode[drp] >= 3 && cc >= Estat[i] && cc <= Estat[i]+dEstat) nmodeE[i][3]=nmodeE[i][3]+1.;
	   if (node[drp] == -1 && cc >= Estat[i] && cc <= Estat[i]+dEstat) nnodeE[i][1]=nnodeE[i][1]+1.;
	    if (node[drp] == 0 && cc >= Estat[i] && cc <= Estat[i]+dEstat) nnodeE[i][2]=nnodeE[i][2]+1.;
	     if (node[drp] >= 1 && cc >= Estat[i] && cc <= Estat[i]+dEstat) nnodeE[i][3]=nnodeE[i][3]+1.;
	   } /* next i */
	  /************************************************mode binning***/
	  for (i=1;i<=3;i=i+1)
	   {
	   m=mode[drp];
	   cc=log10(HKEwkb[drp]);
	   if (cc >= Estat[0] && cc <= Estat[0]+dEstat && m == i) nEm[i][0]=nEm[i][0]+1.;
	    if (cc >= Estat[1] && cc <= Estat[1]+dEstat && m == i) nEm[i][1]=nEm[i][1]+1.;
	     if (cc >= Estat[2] && cc <= Estat[2]+dEstat && m == i) nEm[i][2]=nEm[i][2]+1.;
	      if (cc >= Estat[3] && cc <= Estat[3]+dEstat && m == i) nEm[i][3]=nEm[i][3]+1.;
	   if (istand[drp] == 0 && m == i) nstandm[i][0]=nstandm[i][0]+1.;
 	    if (istand[drp] == 1 && m == i) nstandm[i][1]=nstandm[i][1]+1.;
	   if (iCW[drp] == 0 && m == i) ndownm[i][0]=ndownm[i][0]+1.;
 	    if (iCW[drp] == 1 && m == i) ndownm[i][1]=ndownm[i][1]+1.;
	   if (iCCW[drp] == 0 && m == i) nupm[i][0]=nupm[i][0]+1.;
 	    if (iCCW[drp] == 1 && m == i) nupm[i][1]=nupm[i][1]+1.;
	   if (node[drp] == -1 && m == i) nnodem[i][1]=nnodem[i][1]+1.;
	    if (node[drp] == 0 && m == i) nnodem[i][2]=nnodem[i][2]+1.;
	     if (node[drp] >= 1 && m == i) nnodem[i][3]=nnodem[i][3]+1.;
	   } /* next i */
	  /************************************************topo binning***/
	  for (i=1;i<=6;i=i+1)
	   {
	   cc=log10(varzbb5[drp]);   
	   if (cc >= topostat[i] && cc <= topostat[i]+dtopostat) nt[i]=nt[i]+1.;

	   if (istand[drp] == 0 && cc >= topostat[i] && cc <= topostat[i]+dtopostat) nstandt[i][0]=nstandt[i][0]+1.;
	    if (istand[drp] == 1 && cc >= topostat[i] && cc <= topostat[i]+dtopostat) nstandt[i][1]=nstandt[i][1]+1.;
	   if (iCW[drp] == 0 && cc >= topostat[i] && cc <= topostat[i]+dtopostat) ndownt[i][0]=ndownt[i][0]+1.;
	    if (iCW[drp] == 1 && cc >= topostat[i] && cc <= topostat[i]+dtopostat) ndownt[i][1]=ndownt[i][1]+1.;
	   if (iCCW[drp] == 0 && cc >= topostat[i] && cc <= topostat[i]+dtopostat) nupt[i][0]=nupt[i][0]+1.;
	    if (iCCW[drp] == 1 && cc >= topostat[i] && cc <= topostat[i]+dtopostat) nupt[i][1]=nupt[i][1]+1.;
	   if (mode[drp] == 1 && cc >= topostat[i] && cc <= topostat[i]+dtopostat) nmodet[i][1]=nmodet[i][1]+1.;
	    if (mode[drp] == 2 && cc >= topostat[i] && cc <= topostat[i]+dtopostat) nmodet[i][2]=nmodet[i][2]+1.;
	     if (mode[drp] >= 3 && cc >= topostat[i] && cc <= topostat[i]+dtopostat) nmodet[i][3]=nmodet[i][3]+1.;
	   if (node[drp] == -1 && cc >= topostat[i] && cc <= topostat[i]+dtopostat) nnodet[i][1]=nnodet[i][1]+1.;
	    if (node[drp] == 0 && cc >= topostat[i] && cc <= topostat[i]+dtopostat) nnodet[i][2]=nnodet[i][2]+1.;
	     if (node[drp] == 1 && cc >= topostat[i] && cc <= topostat[i]+dtopostat) nnodet[i][3]=nnodet[i][3]+1.;
	   } /* next i */
	  /*
        } /* end if (lat *	  
	  */
          } /* next drp */
         fclose(fpin);
         fclose(fpE);
         fclose(fprough);
         } /* next k */
        /*****************************************************************stats***/
	for (o=1;o<=Nlat-1;o=o+1)
	 {
	 itotlat[o]=istandlat[o]+idownlat[o]+iuplat[o];
	 imodetotlat[o]=imode1lat[o]+imode2lat[o]+imode3lat[o]+imode4pluslat[o];
         printf("%2.0f%s%2.0f %3.2f %3.2f %3.2f  %3.2f %3.2f %3.2f %3.2f \n",
		 latbin[o],"-",latbin[o]+dlatbin,istandlat[o]/itotlat[o],idownlat[o]/itotlat[o],iuplat[o]/itotlat[o],
		 imode1lat[o]/imodetotlat[o],imode2lat[o]/imodetotlat[o],imode3lat[o]/imodetotlat[o],imode4pluslat[o]/imodetotlat[o]);
	 } /* next o */
        /********************************************************************/
        /********************************************************************/
       /********************************************************************************/
       /**********************************************************PLOT rot***/
       /***************************************set scales****/
       xmin=latmin;
       xmax=latmax;
       dx=dlat;
       dsubx=dsublat;
       delx=xmax-xmin;
       xmid=(xmin+xmax)/2;
 
       ymin=0.;
       ymax=1.2;
       dy=0.2;
       dsuby=0.1;
       dely=ymax-ymin;
       ymid=(ymin+ymax)/2;
 
       xplmin=100.;
       xplmax=500.;
       delxpl=xplmax-xplmin;
       xplmid=(xplmin+xplmax)/2;

       yplmax=700.;
       yplmin=yplmax-200.;
       delypl=yplmax-yplmin;
       yplmid=(yplmin+yplmax)/2;
        /***********************************************make postscript file***/
        dic=0;
	gray=0.0;
        sprintf(sfile,"%s","rotlat.ps");
        fppl=fopen(sfile,"w");
        dic=1;
        printf("%s \n",sfile);
       frame(fppl,&dic,&gray,&xplmin,&xplmax,&yplmin,&yplmax);
        fprintf(fppl,"%s%s%s \n","/",font," findfont 12 scalefont setfont");
        fprintf(fppl,"%5.1f %5.1f %s \n",xplmin+SW12,yplmax+5.,"m");
        /************************************************************axes***/
       frame(fp,&dic,&gray,&xplmin,&xplmax,&yplmin,&yplmax);
        fprintf(fppl,"%s%s%s \n","/",font," findfont 9 scalefont setfont");
        sgnx=sgny=1;
	gray=0.;
        logaxis=0;
        lab=-1;
       upperaxis(fppl,&gray,&red,&green,&blue,&xplmin,&xplmax,&yplmax,&xmin,&xmax,&dx,&dsubx,&sgnx,font,&lab,&logaxis);

       loweraxis(fppl,&gray,&red,&green,&blue,&xplmin,&xplmax,&yplmin,&xmin,&xmax,&dx,&dsubx,&sgnx,font,&lab,&logaxis);
        fprintf(fppl,"%5.1f %5.1f %s \n",xplmid-30.,yplmin-25.,"m");
        fprintf(fppl,"%s \n","(  LATITUDE) show");
	xpl[0]=xmin-dx;
	while (xpl[0] < xmax)
	 {
         xpl[0]=xpl[0]+dx;
	 xx=xpl[0];
        xxpl(&xx,&xmin,&xplmin,&xmax,&xplmax);
	 xpl[1]=xx;
         fprintf(fppl,"%5.1f %5.1f %s \n",xpl[1]-3.*SW12/2.,yplmin-12.,"m");
         if (xpl[0] > 0.) fprintf(fppl,"%s%3.0f%s \n","(",xpl[0],"N) show");
         if (xpl[0] == 0.) fprintf(fppl,"%s%3.0f%s \n","(",xpl[0],") show");
         if (xpl[0] < 0.) fprintf(fppl,"%s%3.0f%s \n","(",-xpl[0],"S) show");
	 } /* end while (xpl */

	lab=0.;
       leftaxis(fppl,&yplmin,&yplmax,&xplmin,&ymin,&ymax,&dy,&dsuby,&sgny,font,&lab,&logaxis);
        fprintf(fppl,"%5.1f %5.1f %s \n",xplmin-30.,yplmid-3.*SW12/2.,"m");
        fprintf(fppl,"%5.2f %s \n",90.,"rotate");
        fprintf(fppl,"%s \n","(ROT) show");
        fprintf(fppl,"%5.2f %s \n",-90.,"rotate");

	lab=-1;
       rightaxis(fppl,&yplmin,&yplmax,&xplmax,&ymin,&ymax,&dy,&dsuby,&sgny,font,&lab,&logaxis);
        /********************************************************************zero line***/
        Npl=2;
        lw=0.5;
        ypl[1]=yplmin;
        ypl[2]=yplmax;
        xx=0.0;
       xxpl(&xx,&xmin,&xplmin,&xmax,&xplmax);
        xpl[1]=xpl[2]=xx;
       dot(fppl,&Npl,&lw,&gray,&red,&green,&blue,&xplmin,&xplmax,&yplmin,&yplmax,xpl,ypl);
        /************************************************************************end PLOT setup***/
        /*****************************************************************************/
        /*****************************************************************************/
        /*****************************************************************************/
        /**********************************************************plotting***/
        Npl=5;
        lw=0.01;
	for (o=1;o<=Nlat-1;o=o+1)
         {
         xx=latbin[o];
         if (xx < xmin) xx=xmin;
         if (xx > xmax) xx=xmax;
        xxpl(&xx,&xmin,&xplmin,&xmax,&xplmax);
	 xpl[1]=xpl[2]=xpl[5]=xx; 
         xx=latbin[o]+dlatbin;
         if (xx < xmin) xx=xmin;
         if (xx > xmax) xx=xmax;
        xxpl(&xx,&xmin,&xplmin,&xmax,&xplmax);
	 xpl[3]=xpl[4]=xx; 

	 ypl[1]=ypl[4]=ypl[5]=yplmin;
         yy=istandlat[o]/itotlat[o];
         if (yy < ymin) yy=ymin;
         if (yy > ymax) yy=ymax;
        yypl(&yy,&ymin,&yplmin,&ymax,&yplmax);
	 ypl[2]=ypl[3]=yy;
	 gray=0.;
         red=green=blue=0.;
	 if (itotlat[o] > 0.)
           shape(fppl,&Npl,&lw,&gray,&red,&green,&blue,&xplmin,&xplmax,&yplmin,&yplmax,xpl,ypl);
         fprintf(fppl,"%5.1f %5.1f %s \n",(xpl[2]+xpl[3])/2.-3.*SW12/2.,yplmax+5.,"m");
         fprintf(fppl,"%s%3.0f%s \n","(",itotlat[o],") show");

	 ypl[1]=ypl[4]=ypl[5]=ypl[2];
         yy=(istandlat[o]+idownlat[o])/itotlat[o];
         if (yy < ymin) yy=ymin;
         if (yy > ymax) yy=ymax;
        yypl(&yy,&ymin,&yplmin,&ymax,&yplmax);
	 ypl[2]=ypl[3]=yy;
	 gray=-1.;
         red=green=blue=0.;
	 blue=1.0;
	 if (itotlat[o] > 0.)
           shape(fppl,&Npl,&lw,&gray,&red,&green,&blue,&xplmin,&xplmax,&yplmin,&yplmax,xpl,ypl);

	 ypl[1]=ypl[4]=ypl[5]=ypl[2];
         yy=1.;
         if (yy < ymin) yy=ymin;
         if (yy > ymax) yy=ymax;
        yypl(&yy,&ymin,&yplmin,&ymax,&yplmax);
	 ypl[2]=ypl[3]=yy;
	 gray=-1.;
         red=green=blue=0.;
	 red=1.0;
	 if (itotlat[o] > 0.)
           shape(fppl,&Npl,&lw,&gray,&red,&green,&blue,&xplmin,&xplmax,&yplmin,&yplmax,xpl,ypl);
         } /* next o */
	/*
	printf("%s %3.2f %3.2f \n","Emn/mx:",Emn,Emx);
	printf("%s %3.2f %3.2f \n","roughmn/mx:",roughmn,roughmx);
	*/
	/*************************************************re-frame***/
	gray=0.;
       frame(fppl,&dic,&gray,&xplmin,&xplmax,&yplmin,&yplmax);
        lab=-1;
       upperaxis(fppl,&gray,&red,&green,&blue,&xplmin,&xplmax,&yplmax,&xmin,&xmax,&dx,&dsubx,&sgnx,font,&lab,&logaxis);
       loweraxis(fppl,&gray,&red,&green,&blue,&xplmin,&xplmax,&yplmin,&xmin,&xmax,&dx,&dsubx,&sgnx,font,&lab,&logaxis);
       leftaxis(fppl,&yplmin,&yplmax,&xplmin,&ymin,&ymax,&dy,&dsuby,&sgny,font,&lab,&logaxis);
       rightaxis(fppl,&yplmin,&yplmax,&xplmax,&ymin,&ymax,&dy,&dsuby,&sgny,font,&lab,&logaxis);
        /********************************************************************zero line***/
        Npl=2;
        lw=0.5;
        ypl[1]=yplmin;
        ypl[2]=yplmax;
        xx=180.0;
       xxpl(&xx,&xmin,&xplmin,&xmax,&xplmax);
        xpl[1]=xpl[2]=xx;
       dot(fppl,&Npl,&lw,&gray,&red,&green,&blue,&xplmin,&xplmax,&yplmin,&yplmax,xpl,ypl);

        xpl[1]=xplmin;
        xpl[2]=xplmax;
        yy=0.0;
       yypl(&yy,&ymin,&yplmin,&ymax,&yplmax);
        ypl[1]=ypl[2]=yy;
       dot(fppl,&Npl,&lw,&gray,&red,&green,&blue,&xplmin,&xplmax,&yplmin,&yplmax,xpl,ypl);
        /**************************************************wrapup***/
        fprintf(fppl,"\n %s \n","showpage");
        fclose(fppl);
       /********************************************************************************/
       /**********************************************************PLOT mode***/
       /***************************************set scales****/
       xmin=latmin;
       xmax=latmax;
       dx=dlat;
       dsubx=dsublat;
       delx=xmax-xmin;
       xmid=(xmin+xmax)/2;
 
       ymin=0.;
       ymax=1.2;
       dy=0.2;
       dsuby=0.1;
       dely=ymax-ymin;
       ymid=(ymin+ymax)/2;
 
       xplmin=100.;
       xplmax=500.;
       delxpl=xplmax-xplmin;
       xplmid=(xplmin+xplmax)/2;

       yplmax=700.;
       yplmin=yplmax-200.;
       delypl=yplmax-yplmin;
       yplmid=(yplmin+yplmax)/2;
        /***********************************************make postscript file***/
        dic=0;
	gray=0.0;
        sprintf(sfile,"%s","modelat.ps");
        fppl=fopen(sfile,"w");
        dic=1;
        printf("%s \n",sfile);
       frame(fppl,&dic,&gray,&xplmin,&xplmax,&yplmin,&yplmax);
        fprintf(fppl,"%s%s%s \n","/",font," findfont 12 scalefont setfont");
        fprintf(fppl,"%5.1f %5.1f %s \n",xplmin+SW12,yplmax+5.,"m");
        /************************************************************axes***/
       frame(fp,&dic,&gray,&xplmin,&xplmax,&yplmin,&yplmax);
        fprintf(fppl,"%s%s%s \n","/",font," findfont 9 scalefont setfont");
        sgnx=sgny=1;
	gray=0.;
        logaxis=0;
        lab=-1;
       upperaxis(fppl,&gray,&red,&green,&blue,&xplmin,&xplmax,&yplmax,&xmin,&xmax,&dx,&dsubx,&sgnx,font,&lab,&logaxis);

       loweraxis(fppl,&gray,&red,&green,&blue,&xplmin,&xplmax,&yplmin,&xmin,&xmax,&dx,&dsubx,&sgnx,font,&lab,&logaxis);
        fprintf(fppl,"%5.1f %5.1f %s \n",xplmid-30.,yplmin-25.,"m");
        fprintf(fppl,"%s \n","(  LATITUDE) show");
	xpl[0]=xmin-dx;
	while (xpl[0] < xmax)
	 {
         xpl[0]=xpl[0]+dx;
	 xx=xpl[0];
        xxpl(&xx,&xmin,&xplmin,&xmax,&xplmax);
	 xpl[1]=xx;
         fprintf(fppl,"%5.1f %5.1f %s \n",xpl[1]-3.*SW12/2.,yplmin-12.,"m");
         if (xpl[0] > 0.) fprintf(fppl,"%s%3.0f%s \n","(",xpl[0],"N) show");
         if (xpl[0] == 0.) fprintf(fppl,"%s%3.0f%s \n","(",xpl[0],") show");
         if (xpl[0] < 0.) fprintf(fppl,"%s%3.0f%s \n","(",-xpl[0],"S) show");
	 } /* end while (xpl */

	lab=0.;
       leftaxis(fppl,&yplmin,&yplmax,&xplmin,&ymin,&ymax,&dy,&dsuby,&sgny,font,&lab,&logaxis);
        fprintf(fppl,"%5.1f %5.1f %s \n",xplmin-30.,yplmid-3.*SW12/2.,"m");
        fprintf(fppl,"%5.2f %s \n",90.,"rotate");
        fprintf(fppl,"%s \n","(MODE) show");
        fprintf(fppl,"%5.2f %s \n",-90.,"rotate");

	lab=-1;
       rightaxis(fppl,&yplmin,&yplmax,&xplmax,&ymin,&ymax,&dy,&dsuby,&sgny,font,&lab,&logaxis);
        /********************************************************************zero line***/
        Npl=2;
        lw=0.5;
        ypl[1]=yplmin;
        ypl[2]=yplmax;
        xx=0.0;
       xxpl(&xx,&xmin,&xplmin,&xmax,&xplmax);
        xpl[1]=xpl[2]=xx;
       dot(fppl,&Npl,&lw,&gray,&red,&green,&blue,&xplmin,&xplmax,&yplmin,&yplmax,xpl,ypl);
        /************************************************************************end PLOT setup***/
        /*****************************************************************************/
        /*****************************************************************************/
        /*****************************************************************************/
        /**********************************************************plotting***/
        Npl=5;
        lw=0.01;
	for (o=1;o<=Nlat-1;o=o+1)
         {
         xx=latbin[o];
         if (xx < xmin) xx=xmin;
         if (xx > xmax) xx=xmax;
        xxpl(&xx,&xmin,&xplmin,&xmax,&xplmax);
	 xpl[1]=xpl[2]=xpl[5]=xx; 
         xx=latbin[o]+dlatbin;
         if (xx < xmin) xx=xmin;
         if (xx > xmax) xx=xmax;
        xxpl(&xx,&xmin,&xplmin,&xmax,&xplmax);
	 xpl[3]=xpl[4]=xx; 

	 ypl[1]=ypl[4]=ypl[5]=yplmin;
         yy=imode1lat[o]/imodetotlat[o];
         if (yy < ymin) yy=ymin;
         if (yy > ymax) yy=ymax;
        yypl(&yy,&ymin,&yplmin,&ymax,&yplmax);
	 ypl[2]=ypl[3]=yy;
	 gray=0.;
         red=green=blue=0.;
	 if (imodetotlat[o] > 0.)
           shape(fppl,&Npl,&lw,&gray,&red,&green,&blue,&xplmin,&xplmax,&yplmin,&yplmax,xpl,ypl);
         fprintf(fppl,"%5.1f %5.1f %s \n",(xpl[2]+xpl[3])/2.-3.*SW12/2.,yplmax+5.,"m");
         fprintf(fppl,"%s%3.0f%s \n","(",imodetotlat[o],") show");

	 ypl[1]=ypl[4]=ypl[5]=ypl[2];
         yy=(imode1lat[o]+imode2lat[o])/imodetotlat[o];
         if (yy < ymin) yy=ymin;
         if (yy > ymax) yy=ymax;
        yypl(&yy,&ymin,&yplmin,&ymax,&yplmax);
	 ypl[2]=ypl[3]=yy;
	 gray=-1.;
         red=green=blue=0.;
	 blue=1.0;
	 if (imodetotlat[o] > 0.)
           shape(fppl,&Npl,&lw,&gray,&red,&green,&blue,&xplmin,&xplmax,&yplmin,&yplmax,xpl,ypl);

	 ypl[1]=ypl[4]=ypl[5]=ypl[2];
         yy=(imode1lat[o]+imode2lat[o]+imode3lat[o])/imodetotlat[o];
         if (yy < ymin) yy=ymin;
         if (yy > ymax) yy=ymax;
        yypl(&yy,&ymin,&yplmin,&ymax,&yplmax);
	 ypl[2]=ypl[3]=yy;
	 gray=-1.;
         red=green=blue=0.;
	 green=1.0;
	 if (imodetotlat[o] > 0.)
           shape(fppl,&Npl,&lw,&gray,&red,&green,&blue,&xplmin,&xplmax,&yplmin,&yplmax,xpl,ypl);

	 ypl[1]=ypl[4]=ypl[5]=ypl[2];
         yy=1.;
         if (yy < ymin) yy=ymin;
         if (yy > ymax) yy=ymax;
        yypl(&yy,&ymin,&yplmin,&ymax,&yplmax);
	 ypl[2]=ypl[3]=yy;
	 gray=-1.;
         red=green=blue=0.;
	 red=1.0;
	 if (itotlat[o] > 0.)
           shape(fppl,&Npl,&lw,&gray,&red,&green,&blue,&xplmin,&xplmax,&yplmin,&yplmax,xpl,ypl);
         } /* next o */
	/*************************************************re-frame***/
	gray=0.;
       frame(fppl,&dic,&gray,&xplmin,&xplmax,&yplmin,&yplmax);
        lab=-1;
       upperaxis(fppl,&gray,&red,&green,&blue,&xplmin,&xplmax,&yplmax,&xmin,&xmax,&dx,&dsubx,&sgnx,font,&lab,&logaxis);
       loweraxis(fppl,&gray,&red,&green,&blue,&xplmin,&xplmax,&yplmin,&xmin,&xmax,&dx,&dsubx,&sgnx,font,&lab,&logaxis);
       leftaxis(fppl,&yplmin,&yplmax,&xplmin,&ymin,&ymax,&dy,&dsuby,&sgny,font,&lab,&logaxis);
       rightaxis(fppl,&yplmin,&yplmax,&xplmax,&ymin,&ymax,&dy,&dsuby,&sgny,font,&lab,&logaxis);
        /********************************************************************zero line***/
        Npl=2;
        lw=0.5;
        ypl[1]=yplmin;
        ypl[2]=yplmax;
        xx=180.0;
       xxpl(&xx,&xmin,&xplmin,&xmax,&xplmax);
        xpl[1]=xpl[2]=xx;
       dot(fppl,&Npl,&lw,&gray,&red,&green,&blue,&xplmin,&xplmax,&yplmin,&yplmax,xpl,ypl);

        xpl[1]=xplmin;
        xpl[2]=xplmax;
        yy=0.0;
       yypl(&yy,&ymin,&yplmin,&ymax,&yplmax);
        ypl[1]=ypl[2]=yy;
       dot(fppl,&Npl,&lw,&gray,&red,&green,&blue,&xplmin,&xplmax,&yplmin,&yplmax,xpl,ypl);
        /**************************************************wrapup***/
        fprintf(fppl,"\n %s \n","showpage");
        } /*END*/