#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;
 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],nsurf[5],nbott[5];
 float nstandE[5][5],ndownE[5][5],nupE[5][5],nmodeE[5][5],nnodeE[5][5],nsurfE[5][5],nbottE[5][5];
 float nEm[5][5],nstandm[5][5],ndownm[5][5],nupm[5][5],nnodem[5][5],nsurfm[5][5],nbottm[5][5];
 float nstandt[8][8],ndownt[8][8],nupt[8][8],nmodet[8][8],nnodet[8][8],nsurft[8][8],nbottt[8][8];

 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("rotmap.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]=nsurf[i]=nbott[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]=nsurfE[i][j]=nbottE[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]=nsurfm[i][j]=nbottm[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]=nsurft[i][j]=nbottt[i][j]=0.;
  } /* next i */
 /**********************************************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);
    /********************************************************************************/
    /********************************************************************************/
    /**********************************************************PLOT map***/
       /***************************************set scales****/
       xmin=lonmin;
       xmax=lonmax;
       dx=dlon;
       dsubx=dsublon;
       delx=xmax-xmin;
       xmid=(xmin+xmax)/2;
 
       ymin=latmin;
       ymax=latmax;
       dy=dlat;
       dsuby=dsublat;
       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;
       chi=1;
       chf=14;
       for (ch=chi;ch<=chf;ch=ch+1)
        {
        /*************************************************make postscript file***/
        dic=0;
	gray=0.0;
        if (ch == 1) sprintf(sfile,"%s","rotmapE.ps");
        if (ch == 2) sprintf(sfile,"%s","rotmapstand.ps");
        if (ch == 3) sprintf(sfile,"%s","rotmapnode.ps");
        if (ch == 4) sprintf(sfile,"%s","rotmaprot.ps");
        if (ch == 5) sprintf(sfile,"%s","rotmapmode.ps");
        if (ch == 6) sprintf(sfile,"%s","rotmaprough.ps");
        if (ch == 7) sprintf(sfile,"%s","rotmapsurfrefl.ps");
        if (ch == 8) sprintf(sfile,"%s","rotmapbottrefl.ps");
        if (ch == 9) sprintf(sfile,"%s","rotmapsurfintens.ps");
        if (ch == 10) sprintf(sfile,"%s","rotmapbottintens.ps");
        if (ch == 11) sprintf(sfile,"%s","rotmapmode1.ps");
        if (ch == 12) sprintf(sfile,"%s","rotmapmode2.ps");
        if (ch == 13) sprintf(sfile,"%s","rotmapmode3.ps");
        if (ch == 14) sprintf(sfile,"%s","rotmapmode4plus.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");
	if (ch == 1) fprintf(fppl,"%s \n","(E^1/2) show");
	if (ch == 2) fprintf(fppl,"%s \n","(E^1/2 standing) show");
	if (ch == 3) fprintf(fppl,"%s \n","(E^1/2 node) show");
	if (ch == 4) fprintf(fppl,"%s \n","(E^1/2 down/up) show");
	if (ch == 5) fprintf(fppl,"%s \n","(E^1/2 mode) show");
	if (ch == 6) fprintf(fppl,"%s \n","(roughness) show");
	if (ch == 7) fprintf(fppl,"%s \n","(E^1/2 surface reflection) show");
	if (ch == 8) fprintf(fppl,"%s \n","(roughness bottom reflection) show");
	if (ch == 9) fprintf(fppl,"%s \n","(E^1/2 surface intensification) show");
	if (ch == 10) fprintf(fppl,"%s \n","(roughness bottom intensification) show");
	if (ch == 11) fprintf(fppl,"%s \n","(E^1/2 mode 1) show");
	if (ch == 12) fprintf(fppl,"%s \n","(E^1/2 mode 2) show");
	if (ch == 13) fprintf(fppl,"%s \n","(E^1/2 mode 3) show");
	if (ch == 14) fprintf(fppl,"%s \n","(E^1/2 mode >= 4) show");
        /***************************************************bathymetry gray bar***/
        itopo=1;
        cplmin=xplmax+5;
        cplmax=cplmin+10;
        delcpl=cplmax-cplmin;
        cplmid=(cplmin+cplmax)/2;
        dic=1;
        gray=0.0;
       frame(fp,&dic,&gray,&cplmin,&cplmax,&yplmin,&yplmax);
        Npl=5;
        lw=0.5;
        gray=0.0;
        red=green=blue=0.0;
        if (clr == 1) gray=-1.;
        xpl[1]=xpl[4]=xpl[5]=cplmin;
        xpl[2]=xpl[3]=cplmax;
        Nrainsub=Nrain;
        for (o=1;o<=Nrainsub;o=o+1)
         {
         ypl[1]=ypl[2]=ypl[5]=yplmin+(o-1)*delypl/Nrainsub;
         ypl[3]=ypl[4]=yplmin+o*delypl/Nrainsub;
         gray=1.0-0.8*o/Nrain;
         if (clr == 1) red=rainbo[o][1];
         if (clr == 1) green=rainbo[o][2];
         if (clr == 1) blue=rainbo[o][3];
        shape(fp,&Npl,&lw,&gray,&red,&green,&blue,&cplmin,&cplmax,&yplmin,&yplmax,xpl,ypl);
         } /* next o */
        fprintf(fp,"%s%s%s \n","/",font," findfont 9 scalefont setfont");
        Npl=2;
        red=green=blue=0.0;
        if (itopo == 0)
         {
         zbmn=-6000.;
         zbmx=0;
         dsubcntr=1000;
         for (o=-6000;o<=0;o=o+dsubcntr)
          {
          gray=1.0;
          xpl[1]=cplmax;
          xpl[2]=cplmin;
          ypl[1]=ypl[2]=yplmin+delypl*(o-zbmn)/(zbmx-zbmn);
          if (cntr == 1)
            solid(fp,&Npl,&lw,&gray,&red,&green,&blue,&cplmin,&cplmax,&yplmin,&yplmax,xpl,ypl);
          gray=0.0;
          xpl[1]=xpl[2]=cplmax;
         solid(fp,&Npl,&lw,&gray,&red,&green,&blue,&cplmin,&cplmax,&yplmin,&yplmax,xpl,ypl);
          fprintf(fp,"%5.1f %5.1f %s \n",xpl[1]+3,ypl[1]-5,"m");
          if (o != 0) fprintf(fp,"%s%d%s \n","(",-o,") show");
          if (o == 0) fprintf(fp,"%s%d%s \n","(",-o," m) show");
          } /* next o */
         } /* end if (itopo */
        if (itopo >= 1)
         {
         zbmn=-0.04;
         zbmx=0.04;
         for (o=-4;o<=4;o=o+2)
          {
          gray=1.0;
          xpl[1]=cplmax;
          xpl[2]=cplmin;
          ypl[1]=ypl[2]=yplmin+delypl*(o/100.-zbmn)/(zbmx-zbmn);
          if (cntr == 1)
            solid(fp,&Npl,&lw,&gray,&red,&green,&blue,&cplmin,&cplmax,&yplmin,&yplmax,xpl,ypl);
          gray=0.0;
          xpl[1]=xpl[2]=cplmax;
         solid(fp,&Npl,&lw,&gray,&red,&green,&blue,&cplmin,&cplmax,&yplmin,&yplmax,xpl,ypl);
          fprintf(fp,"%5.1f %5.1f %s \n",xpl[1]+3,ypl[1]-5,"m");
          xx=o/100.;
          if (o != 10) fprintf(fp,"%s%3.2f%s \n","(",xx,") show");
          if (o == 10) fprintf(fp,"%s%3.2f%s \n","(",xx," m/m) show");
          } /* next o */
         } /* end if (itopo */
        dic=0;
        gray=0.0;
       frame(fp,&dic,&gray,&cplmin,&cplmax,&yplmin,&yplmax);
        fprintf(fp,"%5.1f %5.1f %s \n",cplmax+1.5*SW12,(yplmin+yplmax)/2.-5.,"m");
        fprintf(fp,"%s \n","(      h) show");
        fprintf(fp,"%5.1f %5.1f %s \n",cplmax+1.5*SW12,(yplmin+yplmax)/2.-9.,"m");
        fprintf(fp,"%s \n","(        x) show");
        /************************************************************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","(LONGITUDE) 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. && xpl[0] < 180.) fprintf(fppl,"%s%3.0f%s \n","(",xpl[0],"E) show");
         if (xpl[0] == 0.) fprintf(fppl,"%s%3.0f%s \n","(",xpl[0],") show");
         if (xpl[0] == 180.) fprintf(fppl,"%s%3.0f%s \n","(",xpl[0],") show");
         if (xpl[0] > 180. && xpl[0] < 360.) fprintf(fppl,"%s%3.0f%s \n","(",360.-xpl[0],"W) show");
         if (xpl[0] == 360.) fprintf(fppl,"%s%3.0f%s \n","(",360.-xpl[0],") show");
	 } /* end while (xpl */

       leftaxis(fppl,&yplmin,&yplmax,&xplmin,&ymin,&ymax,&dy,&dsuby,&sgny,font,&lab,&logaxis);
        fprintf(fppl,"%5.1f %5.1f %s \n",xplmin-30.,yplmid-9.*SW12/2.,"m");
        fprintf(fppl,"%5.2f %s \n",90.,"rotate");
        fprintf(fppl,"%s \n","(LATITUDE) show");
        fprintf(fppl,"%5.2f %s \n",-90.,"rotate");
	ypl[0]=ymin-dy;
	while (ypl[0] < ymax)
	 {
         ypl[0]=ypl[0]+dy;
	 yy=ypl[0];
        yypl(&yy,&ymin,&yplmin,&ymax,&yplmax);
	 ypl[1]=yy;
         fprintf(fppl,"%5.1f %5.1f %s \n",xplmin-4.*SW12,ypl[1]-5.,"m");
         if (ypl[0] < 0.) fprintf(fppl,"%s%3.0f%s \n","(",-ypl[0],"S) show");
         if (ypl[0] == 0.) fprintf(fppl,"%s%3.0f%s \n","(",ypl[0],") show");
         if (ypl[0] > 0.) fprintf(fppl,"%s%3.0f%s \n","(",ypl[0],"N) show");
	 } /* end while (xpl */
       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);
       /*********************************************************************bathymetry***/
       /**************************************************input Smith and Sandwell bathy***/  
       test=0;
       if (test == 0)
        {
        i=400;
        while (i > lonmin) i=i-40;
        ilonmin=i;
        lonmin=ilonmin;
        i=-400;
        while (i < lonmax) i=i+40;
        ilonmax=i;
        lonmax=ilonmax;
        for (ilon=ilonmin;ilon<=ilonmax-40;ilon=ilon+40)
         {
         i=ilon;
         while (i < 0) i=i+360;
         while (i > 360) i=i-360;

         sprintf(topodir,"%s%d%s%d%s","/data/kunze/oldarchive/archive/ssbathy/dir",i,"E",i+40,"E");
         for (j=-75;j<=70;j=j+5)
          {
          if ((j >= latmin && j <= latmax) || (j+5 >= latmin && j+5 <= latmax))
           {
           if (j < 0) sprintf(ifile,"%s%s%d%s%d%s%d%s%d%s",topodir,"/ss",-j-5,"S",-j,"S",i,"E",i+40,"E.topo");
           if (j >= 0) sprintf(ifile,"%s%s%d%s%d%s%d%s%d%s",topodir,"/ss",j+5,"N",j,"N",i,"E",i+40,"E.topo");
           fpin=fopen(ifile,"r");
           fscanf(fpin,"%d %d",&Nydata,&Nxdata);
           /*
           printf("  %s %d %d \n",ifile,Nxdata,Nydata);
           */
           for (jj=1;jj<=Nydata;jj=jj+1)
            {
            for (ii=1;ii<=Nxdata;ii=ii+1)
             {
             fscanf(fpin,"%f %f %f",&lonb[ii],&latb[jj],&zbb[ii][jj]);
             } /* next ii */
            } /* next jj */
           fclose(fpin);
           /**************************************************color bathymetry depth***/
           if (itopo == 0)
            {
            zbmn=-5500.;
            zbmx=0.0;
            Dib=10;
            for (ii=1+Dib;ii<=Nxdata;ii=ii+Dib)
             {
             if ((lonb[ii] >= lonmin && lonb[ii] <= lonmax) || (lonb[ii]-360. >= lonmin && lonb[ii]-360 <= lonmax))
              {
              lw=0.5;
              gray=0.0;
              red=green=blue=0.0;
              if (clr == 1) gray=-1.0;
              Npl=5;
              for (jj=1+Dib;jj<=Nydata;jj=jj+Dib)
               {
               if (latb[jj] > latmin && latb[jj] < latmax)
                {
                if (lonb[ii] >= lonmin && lonb[ii] <= lonmax) xx=lonb[ii-Dib];
                if (lonb[ii]-360. >= lonmin && lonb[ii]-360. <= lonmax) xx=lonb[ii-Dib]-360.;
                if (xx < xmin) xx=xmin;
                if (xx > xmax) xx=xmax;
               xxpl(&xx,&xmin,&xplmin,&xmax,&xplmax);
                xpl[1]=xpl[2]=xpl[5]=xx;

                if (lonb[ii] >= lonmin && lonb[ii] <= lonmax) xx=lonb[ii];
                if (lonb[ii]-360. >= lonmin && lonb[ii]-360. <= lonmax) xx=lonb[ii]-360.;
                if (xx < xmin) xx=xmin;
                if (xx > xmax) xx=xmax;
               xxpl(&xx,&xmin,&xplmin,&xmax,&xplmax);
                xpl[3]=xpl[4]=xx;
  
                yy=latb[jj-Dib];
                if (yy < ymin) yy=ymin;
                if (yy > ymax) yy=ymax;
               yypl(&yy,&ymin,&yplmin,&ymax,&yplmax);
                ypl[1]=ypl[4]=ypl[5]=yy;
                yy=latb[jj];
                if (yy < ymin) yy=ymin;
                if (yy > ymax) yy=ymax;
               yypl(&yy,&ymin,&yplmin,&ymax,&yplmax);
                ypl[2]=ypl[3]=yy;
 
                zbb[0][0]=(zbb[ii][jj]+zbb[ii-Dib][jj]+zbb[ii][jj-Dib]+zbb[ii-Dib][jj-Dib])/4.;
                zbb[0][0]=zbb[ii-Dib/2][jj-Dib/2];
                xx=1.+(Nrainsub-1)*(zbb[0][0]-zbmn)/(zbmx-zbmn);
                o=1;
                while (o < xx && o < Nrainsub) o=o+1;
                gray=0.5-0.5*o/Nrain;
                if (clr == 1) red=rainbo[o][1];
                if (clr == 1) green=rainbo[o][2];
                if (clr == 1) blue=rainbo[o][3];
                if (zbb[ii-Dib/2][jj-Dib/2] != -9999)
                  shape(fp,&Npl,&lw,&gray,&red,&green,&blue,&xplmin,&xplmax,&yplmin,&yplmax,xpl,ypl);
                } /* end if (latb */
               } /* next jj */
              } /* end if (lonb */
             } /* next ii */
            } /* end if (itopo */
          /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * zonal slope * * */
          if (itopo >= 1)
           {
           zbmn=-0.04;
           zbmx=0.04;
      	   Dib=10;
          for (ii=1+Dib;ii<=Nxdata;ii=ii+Dib)
           {
           if ((lonb[ii] >= lonmin && lonb[ii] <= lonmax) || (lonb[ii]-360. >= lonmin && lonb[ii]-360 <= lonmax))
            {
            lw=0.5;
            gray=0.0;
            red=green=blue=0.0;
            if (clr == 1) gray=-1.0;
            Npl=5;
            for (jj=1+Dib;jj<=Nydata;jj=jj+Dib)
             {
             if (latb[jj] > latmin && latb[jj] < latmax)
              {
              if (lonb[ii] >= lonmin && lonb[ii] <= lonmax) xx=lonb[ii-Dib];
              if (lonb[ii]-360. >= lonmin && lonb[ii]-360. <= lonmax) xx=lonb[ii-Dib]-360.;
              if (xx < xmin) xx=xmin;
              if (xx > xmax) xx=xmax;
             xxpl(&xx,&xmin,&xplmin,&xmax,&xplmax);
              xpl[1]=xpl[2]=xpl[5]=xx;
              if (lonb[ii] >= lonmin && lonb[ii] <= lonmax) xx=lonb[ii];
              if (lonb[ii]-360. >= lonmin && lonb[ii]-360. <= lonmax) xx=lonb[ii]-360.;
              if (xx < xmin) xx=xmin;
              if (xx > xmax) xx=xmax;
             xxpl(&xx,&xmin,&xplmin,&xmax,&xplmax);
              xpl[3]=xpl[4]=xx;
 
              yy=latb[jj-Dib];
              if (yy < ymin) yy=ymin;
              if (yy > ymax) yy=ymax;
             yypl(&yy,&ymin,&yplmin,&ymax,&yplmax);
              ypl[1]=ypl[4]=ypl[5]=yy;
              yy=latb[jj];
              if (yy < ymin) yy=ymin;
              if (yy > ymax) yy=ymax;
             yypl(&yy,&ymin,&yplmin,&ymax,&yplmax);
              ypl[2]=ypl[3]=yy;
 
   	      coslat=cos(PI*latb[jj]/180.);
	      lonkm=latkm*coslat;
              if (itopo == 1) zbb[0][0]=(zbb[ii][jj]-zbb[ii-Dib][jj])/(lonb[ii]-lonb[ii-Dib])/lonkm/1000.;
              if (itopo == 2) zbb[0][0]=(zbb[ii][jj]-zbb[ii][jj-Dib])/(latb[jj]-latb[jj-Dib])/latkm/1000.;
              xx=1.+(Nrainsub-1)*(zbb[0][0]-zbmn)/(zbmx-zbmn);
              o=1;
              while (o < xx && o < Nrainsub) o=o+1;
              gray=1.0-0.8*o/Nrain;
              if (clr == 1) red=rainbo[o][1];
              if (clr == 1) green=rainbo[o][2];
              if (clr == 1) blue=rainbo[o][3];
              if (zbb[ii-Dib][jj] != -9999 && zbb[ii][jj] != -9999 && zbb[ii][jj-Dib] != -9999)
                shape(fp,&Npl,&lw,&gray,&red,&green,&blue,&xplmin,&xplmax,&yplmin,&yplmax,xpl,ypl);
               } /* end if (latb */
              } /* next jj */
             } /* end if (lonb */
            } /* next ii */
           } /* end if (itopo */
          /**************************************************brown land***/
          if (clr == 0)
           {
            lw=0.5;
            gray=-1.0;
            red=green=blue=0.0;
            red=brown[1];
            green=brown[2];
            blue=brown[3];
            Npl=5;
            for (ii=1+Dib;ii<=Nxdata;ii=ii+Dib)
             {
             if ((lonb[ii] >= lonmin && lonb[ii] <= lonmax) || (lonb[ii]-360. >= lonmin && lonb[ii]-360 <= lonmax))
              {
              for (jj=1+Dib;jj<=Nydata;jj=jj+Dib)
               {
               if (latb[jj] > latmin && latb[jj] < latmax)
                {
                if (lonb[ii] >= lonmin && lonb[ii] <= lonmax) xx=lonb[ii-Dib];
                if (lonb[ii]-360. >= lonmin && lonb[ii]-360. <= lonmax) xx=lonb[ii-Dib]-360.;
                if (xx < xmin) xx=xmin;
                if (xx > xmax) xx=xmax;
               xxpl(&xx,&xmin,&xplmin,&xmax,&xplmax);
                xpl[1]=xpl[2]=xpl[5]=xx;
  
                if (lonb[ii] >= lonmin && lonb[ii] <= lonmax) xx=lonb[ii];
                if (lonb[ii]-360. >= lonmin && lonb[ii]-360. <= lonmax) xx=lonb[ii]-360.;
                if (xx < xmin) xx=xmin;
                if (xx > xmax) xx=xmax;
               xxpl(&xx,&xmin,&xplmin,&xmax,&xplmax);
                xpl[3]=xpl[4]=xx;
  
                yy=latb[jj-Dib];
                if (yy < ymin) yy=ymin;
                if (yy > ymax) yy=ymax;
               yypl(&yy,&ymin,&yplmin,&ymax,&yplmax);
                ypl[1]=ypl[4]=ypl[5]=yy;
  
                yy=latb[jj];
                if (yy < ymin) yy=ymin;
                if (yy > ymax) yy=ymax;
               yypl(&yy,&ymin,&yplmin,&ymax,&yplmax);
                ypl[2]=ypl[3]=yy;
 
                if (zbb[ii][jj]+zbb[ii][jj-Dib]+zbb[ii-Dib][jj]+zbb[ii-Dib][jj-Dib] > 0.)
                  shape(fp,&Npl,&lw,&gray,&red,&green,&blue,&xplmin,&xplmax,&yplmin,&yplmax,xpl,ypl);
                } /* end if (latb */
               } /* next jj */
              } /* end if (lonb */
             } /* next ii */
            } /* end if (clr */
           /***************************************************************contour***/
           if (cntr == 1)
            {
            lw=0.5;
            gray=1.0;
            red=blue=green=0.0;
            Ncntr=5;
            for (ii=1+Dib;ii<=Nxdata;ii=ii+Dib)
             {
             if (lonb[ii] > lonmin && lonb[ii] < lonmax)
              {
              for (jj=1+Dib;jj<=Nydata;jj=jj+Dib)
               {
               if (latb[jj] > latmin && latb[jj] < latmax)
                {
                xx=lonb[ii-Dib];
                if (lonb[ii-Dib] > 180. && xmax <= 0.) xx=lonb[ii-Dib]-360.;
                if (xx < xmin) xx=xmin;
                if (xx > xmax) xx=xmax;
               xxpl(&xx,&xmin,&xplmin,&xmax,&xplmax);
                xcntr[1]=xcntr[2]=xcntr[5]=xx;
  
                xx=lonb[ii];
                if (lonb[ii] > 180. && xmax <= 0.) xx=lonb[ii]-360.;
                if (xx < xmin) xx=xmin;
                if (xx > xmax) xx=xmax;
               xxpl(&xx,&xmin,&xplmin,&xmax,&xplmax);
                xcntr[3]=xcntr[4]=xx;
  
                yy=latb[jj-Dib];
                if (yy < ymin) yy=ymin;
                if (yy > ymax) yy=ymax;
               yypl(&yy,&ymin,&yplmin,&ymax,&yplmax);
                ycntr[1]=ycntr[4]=ycntr[5]=yy;
  
                yy=latb[jj];
                if (yy < ymin) yy=ymin;
                if (yy > ymax) yy=ymax;
               yypl(&yy,&ymin,&yplmin,&ymax,&yplmax);
                ycntr[2]=ycntr[3]=yy;
  
                fcntr[1]=fcntr[5]=-zbb[ii-Dib][jj-Dib];
                fcntr[2]=-zbb[ii-Dib][jj];
                fcntr[3]=-zbb[ii][jj];
                fcntr[4]=-zbb[ii][jj-Dib];
               contour(fp,&Ncntr,&gray,&red,&blue,&green,&xplmin,&xplmax,&yplmin,&yplmax,xcntr,ycntr,fcntr,&dcntr,&dsubcntr);
                } /* end if (latb */
               } /* next jj */
              } /* end if (lonb */
             } /* next ii */
            } /* end if (cntr */
           } /* end if (lonmin */
          } /* next j */
         } /* next ilon */
        } /* end if (test */
       /**********************************************************************************end bathy***/
       /************************************************************************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 */
	  /************************************************pdf binning***/
	  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.;
	  if (surfrefl[drp] == 0 && surfinten[drp] == 0) nsurf[0]=nsurf[0]+1.;
	   if (surfrefl[drp] == 1 || surfinten[drp] == 1) nsurf[1]=nsurf[1]+1.;
	  if (bottrefl[drp] == 0 && bottintens[drp] == 0) nbott[0]=nbott[0]+1.;
	   if (bottrefl[drp] == 1 || bottintens[drp] == 1) nbott[1]=nbott[1]+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.;
	   if (surfrefl[drp] == 0 && surfinten[drp] == 0 && cc >= Estat[i] && cc <= Estat[i]+dEstat) nsurfE[i][0]=nsurfE[i][0]+1.;
	   if ((surfrefl[drp] == 1 || surfinten[drp] == 1) && cc >= Estat[i] && cc <= Estat[i]+dEstat) nsurfE[i][1]=nsurfE[i][1]+1.;
	   if (bottrefl[drp] == 0 && bottintens[drp] == 0 && cc >= Estat[i] && cc <= Estat[i]+dEstat) nbottE[i][0]=nbottE[i][0]+1.;
	   if ((bottrefl[drp] == 1 || bottintens[drp] == 1) && cc >= Estat[i] && cc <= Estat[i]+dEstat) nbottE[i][1]=nbottE[i][1]+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.;
	   if (surfrefl[drp] == 0 && surfinten[drp] == 0 && m == i) nsurfm[i][0]=nsurfm[i][0]+1.;
	   if ((surfrefl[drp] == 1 || surfinten[drp] == 0) && m == i) nsurfm[i][1]=nsurfm[i][1]+1.;
	   if (bottrefl[drp] == 0 && bottintens[drp] == 1 && m == i) nbottm[i][0]=nbottm[i][0]+1.;
	   if ((bottrefl[drp] == 1 || bottintens[drp] == 1) && m == i) nbottm[i][1]=nbottm[i][1]+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.;
	   if (surfrefl[drp] == 0 && surfinten[drp] == 0 && cc >= topostat[i] && cc <= topostat[i]+dtopostat) nsurft[i][0]=nsurft[i][0]+1.;
	   if ((surfrefl[drp] == 1 || surfinten[drp] == 1) && cc >= topostat[i] && cc <= topostat[i]+dtopostat) nsurft[i][1]=nsurft[i][1]+1.;
	   if (bottrefl[drp] == 0 && bottintens[drp] == 0 && cc >= topostat[i] && cc <= topostat[i]+dtopostat) nbottt[i][0]=nbottt[i][0]+1.;
	   if ((bottrefl[drp] == 1 || bottintens[drp] == 1) && cc >= topostat[i] && cc <= topostat[i]+dtopostat) nbottt[i][1]=nbottt[i][1]+1.;
	   } /* next i */
          } /* next drp */
         fclose(fpin);
         fclose(fpE);
         fclose(fprough);
         /********************************************************************/
         /********************************************************************/
         /********************************************************************/
         /**********************************************************plotting***/
         Npl=7;
	 lw=0.01;
	 gray=-1.;
         for (drp=1;drp<=Ndrp[k];drp=drp+1)
          {
          yy=lat[drp];
          if (yy < ymin) yy=ymin;
          if (yy > ymax) yy=ymax;
         yypl(&yy,&ymin,&yplmin,&ymax,&yplmax);
 
          xx=lon[drp];
          if (xx < xmin) xx=xmin;
          if (xx > xmax) xx=xmax;
         xxpl(&xx,&xmin,&xplmin,&xmax,&xplmax);
	  rpl=0.0;
	  gray=-1.;
          red=green=blue=0.;
	  /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * sqrt(E) * * */
	  if (ch == 1 && HKEwkb[drp] > 0.) rpl=sqrt(HKEwkb[drp]/100.);
	  if (ch == 1 && iE[drp] == 1) blue=0.8;
	  if (ch == 1 && iE[drp] == 2) blue=0.6;
	  if (ch == 1 && iE[drp] == 2) red=0.7;
	  if (ch == 1 && iE[drp] >= 3) red=1.0;
	  /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * standing * * */
	  /*
	  if (ch == 2 && istand[drp] == 1) rpl=iE[drp];
	  */
	  if (ch == 2 && HKEwkb[drp] > 0. && istand[drp] == 1) rpl=sqrt(HKEwkb[drp]/100.);
	  if (ch == 2 && iCW[drp] != 0) rpl=0.;
	  if (ch == 2 && iCCW[drp] != 0) rpl=0.;
	  cc=(log10(HKEwkb[drp]))/4.;
	  o=1;
	  while (o < Nrain && o < Nrain*cc) o=o+1;
	  if (ch == 2 && istand[drp] == 1) red=rainbo[o][1];
	  if (ch == 2 && istand[drp] == 1) green=rainbo[o][2];
	  if (ch == 2 && istand[drp] == 1) blue=rainbo[o][3];
	  if (ch == 2 && istand[drp] == 0) rpl=red=green=blue=0.;
	  /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * node * * */
	  /*
	  if (ch == 3 && node[drp] != 0) rpl=1.;
	  */
	  if (ch == 3 && HKEwkb[drp] > 0. && node[drp] != 0) rpl=sqrt(HKEwkb[drp]/100.);
	  if (ch == 3 && node[drp] == 1) red=1.;
	  if (ch == 3 && node[drp] == -1) blue=1.;
	  if (ch == 3 && node[drp] == -1) green=red=0.4;
	  /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * rotary * * */
	  /*
	  if (ch == 4 && istand[drp] != 0) rpl=iE[drp];
	  if (ch == 4 && iCW[drp] != 0) rpl=iE[drp];
	  if (ch == 4 && iCCW[drp] != 0) rpl=iE[drp];
	  */
	  if (ch == 4 && istand[drp] != 0) rpl=0.;
	  if (ch == 4 && HKEwkb[drp] > 0. && iCW[drp] != 0) rpl=sqrt(HKEwkb[drp]/100.);
	  if (ch == 4 && HKEwkb[drp] > 0. && iCCW[drp] != 0) rpl=sqrt(HKEwkb[drp]/100.);
	  if (ch == 4 && iCCW[drp] == 1) red=1.;
	  if (ch == 4 && iCW[drp] == 1) blue=1.;
	  if (ch == 3 && iCW[drp] == 1) green=red=0.4;
	  if (ch == 4 && istand[drp] == 1) green=0.6;
	  /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * mode * * */
	  /*
	  if (ch == 5 && mode[drp] != 0) rpl=iE[drp];
	  */
	  if (ch == 5 && HKEwkb[drp] > 0. && mode[drp] != 0) rpl=sqrt(HKEwkb[drp]/100.);
	  if (ch == 5 && mode[drp] == 1) gray=0.;
	  if (ch == 5 && mode[drp] == 2) blue=1.;
	  if (ch == 5 && mode[drp] == 2) green=red=0.4;
	  if (ch == 5 && mode[drp] == 3) green=0.6;
	  if (ch == 5 && mode[drp] > 3) red=1.0;
	  /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * topo roughness * * */
	  if (ch == 6) 
	   {
	   if (varzbb5[drp] > 0 && HKEwkb[drp] > 0.) rpl=sqrt(HKEwkb[drp]/100.);
	   cc=(log10(sqrt(varzbb5[drp]))-roughmn)/(roughmx-roughmn);
	   o=1;
	   while (o < cc*Nrain && o < Nrain) o=o+1;
	   red=rainbo[o][1];
	   green=rainbo[o][2];
	   blue=rainbo[o][3];
	   } /* end if (ch */
	  /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * surf refl * * */
	  if (ch == 7 && HKEwkb[drp] > 0. && surfrefl[drp] == 1) rpl=sqrt(HKEwkb[drp]/100.);
	  if (ch == 7 && surfrefl[drp] == 1) red=1.0;
	  /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * bottom refl * * */
	  cc=(log10(sqrt(varzbb5[drp]))-1.)/3.;
	  if (ch == 8 && cc > 0. && bottrefl[drp] == 1) rpl=4.*cc;
	  if (ch == 8 && bottrefl[drp] == 1) red=1.0;
	  /* * * * * * * * * * * * * * * * * * * * * * * * * * * surf intensification * * */
	  if (ch == 9 && HKEwkb[drp] > 0. && surfinten[drp] == 1) rpl=sqrt(HKEwkb[drp]/100.);
	  if (ch == 9 && surfinten[drp] == 1) red=1.0;
	  /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * bottom intensification * * */
	  cc=(log10(sqrt(varzbb5[drp]))-1.)/3.;
	  if (ch == 10 && cc > 0. && bottintens[drp] == 1) rpl=4.*cc;
	  if (ch == 10 && bottintens[drp] == 1) red=1.0;
	  /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * mode 1 * * */
	  if (ch == 11 && HKEwkb[drp] > 0. && mode[drp] == 1) rpl=sqrt(HKEwkb[drp]/100.);
	  if (ch == 11 && mode[drp] == 1) red=1.0;
	  /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * mode 2 * * */
	  if (ch == 12 && HKEwkb[drp] > 0. && mode[drp] == 2) rpl=sqrt(HKEwkb[drp]/100.);
	  if (ch == 12 && mode[drp] == 2) red=1.0;
	  /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * mode == 3 * * */
	  if (ch == 13 && HKEwkb[drp] > 0. && mode[drp] == 3) rpl=sqrt(HKEwkb[drp]/100.);
	  if (ch == 13 && mode[drp] == 3) red=1.0;
	  /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * mode > 3 * * */
	  if (ch == 14 && HKEwkb[drp] > 0. && mode[drp] > 3) rpl=sqrt(HKEwkb[drp]/100.);
	  if (ch == 14 && mode[drp] > 3) red=1.0;
	  /************************************************************shape***/
          for (o=1;o<=Npl;o=o+1)
	   { 
	   Orad=2*PI*(o-1)/(Npl-1);
	   cosO=cos(Orad);
	   sinO=sin(Orad);
	   xpl[o]=xx+rpl*cosO;
	   ypl[o]=yy+rpl*sinO;
           } /* next o */
        shape(fppl,&Npl,&lw,&gray,&red,&green,&blue,&xplmin,&xplmax,&yplmin,&yplmax,xpl,ypl);
          } /* next drp */
         } /* next k */
	/*
	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 lines***/
        Npl=2;
        lw=0.8;
        ypl[1]=yplmin;
        ypl[2]=yplmax;
	for (o=1;o<=17;o=o+1)
         {
         xx=20.0*o;
        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);
	 } /* next o */

        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);
	/**********************************************************************************print stats***/
	if (ch == 1)
	 {
	 /**********************************************************************************print Ebin stats***/
         printf("%s \n"," ");
         printf("%s \n","ENERGY");
         rpl=nE[0]+nE[1]+nE[2]+nE[3]+nE[4];
         printf("%s \n","E TOT 0-1 1-2 2-3 3-4");
         printf("%2.1f %2.1f %2.1f %2.1f %2.1f \n",Estat[0],Estat[1],Estat[2],Estat[3],Estat[4]);
         printf("%s %5.0f  %4.3f %4.3f %4.3f %4.3f %4.3f \n","nE TOT :",rpl,nE[0]/rpl,nE[1]/rpl,nE[2]/rpl,nE[3]/rpl,nE[4]/rpl);
         printf("%s \n"," ");

         nstand[4]=nstand[0]+nstand[1];
         printf("%s %5.0f  %4.3f %4.3f \n","nstand TOT no yes:",nstand[4],nstand[0]/nstand[4],nstand[1]/nstand[4]);
         nstandE[0][4]=nstandE[0][0]+nstandE[0][1];
         printf("  %s %5.0f  %4.3f %4.3f \n","nstand-low-E TOT no yes:",nstandE[0][4],nstandE[0][0]/nstandE[0][4],nstandE[0][1]/nstandE[0][4]);
         nstandE[0][4]=nstandE[1][0]+nstandE[1][1];
         printf("  %s %5.0f  %4.3f %4.3f \n","nstand-E1 TOT no yes:",nstandE[0][4],nstandE[1][0]/nstandE[0][4],nstandE[1][1]/nstandE[0][4]);
         nstandE[0][4]=nstandE[2][0]+nstandE[2][1];
         printf("  %s %5.0f  %4.3f %4.3f \n","nstand-E2 TOT no yes:",nstandE[0][4],nstandE[2][0]/nstandE[0][4],nstandE[2][1]/nstandE[0][4]);
         nstandE[0][4]=nstandE[3][0]+nstandE[3][1];
         printf("  %s %5.0f  %4.3f %4.3f \n","nstand-E3 TOT no yes:",nstandE[0][4],nstandE[3][0]/nstandE[0][4],nstandE[3][1]/nstandE[0][4]);
         printf("%s \n"," ");

         ndown[4]=ndown[0]+ndown[1];
         printf("%s %5.0f  %4.3f %4.3f \n","ndown TOT no yes:",ndown[4],ndown[0]/ndown[4],ndown[1]/ndown[4]);
         ndownE[0][4]=ndownE[0][0]+ndownE[0][1];
         printf("  %s %5.0f  %4.3f %4.3f \n","ndown-low-E TOT no yes:",ndownE[0][4],ndownE[0][0]/ndownE[0][4],ndownE[0][1]/ndownE[0][4]);
         ndownE[1][4]=ndownE[1][0]+ndownE[1][1];
         printf("  %s %5.0f  %4.3f %4.3f \n","ndown-E1 TOT no yes:",ndownE[1][4],ndownE[1][0]/ndownE[1][4],ndownE[1][1]/ndownE[1][4]);
         ndownE[2][4]=ndownE[2][0]+ndownE[2][1];
         printf("  %s %5.0f  %4.3f %4.3f \n","ndown-E2 TOT no yes:",ndownE[2][4],ndownE[2][0]/ndownE[2][4],ndownE[2][1]/ndownE[2][4]);
         ndownE[3][4]=ndownE[3][0]+ndownE[3][1];
         printf("  %s %5.0f  %4.3f %4.3f \n","ndown-E3 TOT no yes:",ndownE[3][4],ndownE[3][0]/ndownE[3][4],ndownE[3][1]/ndownE[3][4]);
         printf("%s \n"," ");

         nup[4]=nup[0]+nup[1];
         printf("%s %5.0f  %4.3f %4.3f \n","nup TOT no yes:",nup[4],nup[0]/nup[4],nup[1]/nup[4]);
         nupE[0][4]=nupE[0][0]+nupE[0][1];
         printf("  %s %5.0f  %4.3f %4.3f \n","nup-low-E TOT no yes:",nupE[0][4],nupE[0][0]/nupE[0][4],nupE[0][1]/nupE[0][4]);
         nupE[1][4]=nupE[1][0]+nupE[1][1];
         printf("  %s %5.0f  %4.3f %4.3f \n","nup-E1 TOT no yes:",nupE[1][4],nupE[1][0]/nupE[1][4],nupE[1][1]/nupE[1][4]);
         nupE[2][4]=nupE[2][0]+nupE[2][1];
         printf("  %s %5.0f  %4.3f %4.3f \n","nup-E2 TOT no yes:",nupE[2][4],nupE[2][0]/nupE[2][4],nupE[2][1]/nupE[2][4]);
         nupE[3][4]=nupE[3][0]+nupE[3][1];
         printf("  %s %5.0f  %4.3f %4.3f \n","nup-E3 TOT no yes:",nupE[3][4],nupE[3][0]/nupE[3][4],nupE[3][1]/nupE[3][4]);
         printf("%s \n"," ");

         nmode[0]=nmode[1]+nmode[2]+nmode[3];
         printf("%s %5.0f  %4.3f %4.3f %4.3f \n","nmode TOT 1 2 >3:",nmode[0],nmode[1]/nmode[0],nmode[2]/nmode[0],nmode[3]/nmode[0]);
         nmode[0]=nmodeE[0][1]+nmodeE[0][2]+nmodeE[0][3];
         printf("   %s %5.0f  %4.3f %4.3f %4.3f \n","nmode-low-E TOT 1 2 >3:",nmode[0],nmodeE[0][1]/nmode[0],nmodeE[0][2]/nmode[0],nmodeE[0][3]/nmode[0]);
         nmode[0]=nmodeE[1][1]+nmodeE[1][2]+nmodeE[1][3];
         printf("   %s %5.0f  %4.3f %4.3f %4.3f \n","nmode-E1 TOT 1 2 >3:",nmode[0],nmodeE[1][1]/nmode[0],nmodeE[1][2]/nmode[0],nmodeE[1][3]/nmode[0]);
         nmode[0]=nmodeE[2][1]+nmodeE[2][2]+nmodeE[2][3];
         printf("   %s %5.0f  %4.3f %4.3f %4.3f \n","nmode-E2 TOT 1 2 >3:",nmode[0],nmodeE[2][1]/nmode[0],nmodeE[2][2]/nmode[0],nmodeE[2][3]/nmode[0]);
         nmode[0]=nmodeE[3][1]+nmodeE[3][2]+nmodeE[3][3];
         printf("   %s %5.0f  %4.3f %4.3f %4.3f \n","nmode-E3 TOT 1 2 >3:",nmode[0],nmodeE[3][1]/nmode[0],nmodeE[3][2]/nmode[0],nmodeE[3][3]/nmode[0]);
         printf("%s \n"," ");

         nnode[0]=nnode[1]+nnode[2]+nnode[3];
         printf("%s %5.0f  %4.3f %4.3f %4.3f \n","nnode TOT -1 0 +1:",nnode[0],nnode[1]/nnode[0],nnode[2]/nnode[0],nnode[3]/nnode[0]);
         nnode[0]=nnodeE[0][1]+nnodeE[0][2]+nnodeE[0][3];
         printf("   %s %5.0f  %4.3f %4.3f %4.3f \n","nnode TOT-low-E -1 0 +1:",nnode[0],nnodeE[0][1]/nnode[0],nnodeE[0][2]/nnode[0],nnodeE[0][3]/nnode[0]);
         nnode[0]=nnodeE[1][1]+nnodeE[1][2]+nnodeE[1][3];
         printf("   %s %5.0f  %4.3f %4.3f %4.3f \n","nnode TOT-E1 -1 0 +1:",nnode[0],nnodeE[1][1]/nnode[0],nnodeE[1][2]/nnode[0],nnodeE[1][3]/nnode[0]);
         nnode[0]=nnodeE[2][1]+nnodeE[2][2]+nnodeE[2][3];
         printf("   %s %5.0f  %4.3f %4.3f %4.3f \n","nnode TOT-E2 -1 0 +1:",nnode[0],nnodeE[2][1]/nnode[0],nnodeE[2][2]/nnode[0],nnodeE[2][3]/nnode[0]);
         nnode[0]=nnodeE[3][1]+nnodeE[3][2]+nnodeE[3][3];
         printf("   %s %5.0f  %4.3f %4.3f %4.3f \n","nnode TOT-E3 -1 0 +1:",nnode[0],nnodeE[3][1]/nnode[0],nnodeE[3][2]/nnode[0],nnodeE[3][3]/nnode[0]);
         printf("%s \n"," ");

         rpl=nsurf[0]+nsurf[1];
         printf("%s %5.0f  %4.3f %4.3f \n","nsurf TOT 0 +1:",rpl,nsurf[0]/rpl,nsurf[1]/rpl);
         rpl=nsurfE[0][0]+nsurfE[0][1];
         printf("   %s %5.0f  %4.3f %4.3f \n","nsurf TOT E0 no yes:",
		    rpl,nsurfE[0][0]/rpl,nsurfE[0][1]/rpl);
         rpl=nsurfE[1][0]+nsurfE[1][1];
         printf("   %s %5.0f  %4.3f %4.3f \n","nsurf TOT E1 no yes:",
		    rpl,nsurfE[1][0]/rpl,nsurfE[1][1]/rpl);
         rpl=nsurfE[2][0]+nsurfE[2][1];
         printf("   %s %5.0f  %4.3f %4.3f \n","nsurf TOT E2 no yes:",
		    rpl,nsurfE[2][0]/rpl,nsurfE[2][1]/rpl);
         rpl=nsurfE[3][0]+nsurfE[3][1];
         printf("   %s %5.0f  %4.3f %4.3f \n","nsurf TOT E3 no yes:",
		    rpl,nsurfE[3][0]/rpl,nsurfE[3][1]/rpl);
         printf("%s \n"," ");

         nbott[3]=nbott[0]+nbott[1];
         printf("%s %5.0f  %4.3f %4.3f \n","nbott TOT 0 +1:",nbott[3],nbott[0]/nbott[3],nbott[1]/nbott[3]);
         rpl=nbottE[0][0]+nbottE[0][1];
         printf("   %s %5.0f  %4.3f %4.3f \n","nbott TOT E0 no yes:",
		    rpl,nbottE[0][0]/rpl,nbottE[0][1]/rpl);
         rpl=nbottE[1][0]+nbottE[1][1];
         printf("   %s %5.0f  %4.3f %4.3f \n","nbott TOT E1 no yes:",
		    rpl,nbottE[1][0]/rpl,nbottE[1][1]/rpl);
         rpl=nbottE[2][0]+nbottE[2][1];
         printf("   %s %5.0f  %4.3f %4.3f \n","nbott TOT E2 no yes:",
		    rpl,nbottE[2][0]/rpl,nbottE[2][1]/rpl);
         rpl=nbottE[3][0]+nbottE[3][1];
         printf("   %s %5.0f  %4.3f %4.3f \n","nbott TOT E3 no yes:",
		    rpl,nbottE[3][0]/rpl,nbottE[3][1]/rpl);
         printf("%s \n"," ");
	 /**********************************************************************************print modebin stats***/
         printf("%s \n","MODES");
         rpl=nEm[1][0]+nEm[1][1]+nEm[1][2]+nEm[1][3];
         printf("  %s %5.0f  %4.3f %4.3f %4.3f %4.3f \n","nstand-mode 1 TOT E0 1 2 3:",rpl,nEm[1][0]/rpl,nEm[1][1]/rpl,nEm[1][2]/rpl,nEm[1][3]/rpl);
         rpl=nEm[2][0]+nEm[2][1]+nEm[2][2]+nEm[2][3];
         printf("  %s %5.0f  %4.3f %4.3f %4.3f %4.3f \n","nstand-mode 2 TOT E0 1 2 3:",rpl,nEm[2][0]/rpl,nEm[2][1]/rpl,nEm[2][2]/rpl,nEm[2][3]/rpl);
         rpl=nEm[3][0]+nEm[3][1]+nEm[3][2]+nEm[3][3];
         printf("  %s %5.0f  %4.3f %4.3f %4.3f %4.3f \n","nstand-mode 3 TOT E0 1 2 3:",rpl,nEm[3][0]/rpl,nEm[3][1]/rpl,nEm[3][2]/rpl,nEm[3][3]/rpl);
         printf("%s \n"," ");

         rpl=nstandm[1][0]+nstandm[1][1];
         printf("  %s %5.0f  %4.3f %4.3f \n","nstand-mode 1 TOT no yes:",rpl,nstandm[1][0]/rpl,nstandm[1][1]/rpl);
         rpl=nstandm[2][0]+nstandm[2][1];
         printf("  %s %5.0f  %4.3f %4.3f \n","nstand-mode 2 TOT no yes:",rpl,nstandm[2][0]/rpl,nstandm[2][1]/rpl);
         rpl=nstandm[3][0]+nstandm[3][1];
         printf("  %s %5.0f  %4.3f %4.3f \n","nstand-mode 3 TOT no yes:",rpl,nstandm[3][0]/rpl,nstandm[3][1]/rpl);
         printf("%s \n"," ");

         rpl=ndownm[1][0]+ndownm[1][1];
         printf("  %s %5.0f  %4.3f %4.3f \n","ndown-mode 1 TOT no yes:",rpl,ndownm[1][0]/rpl,ndownm[1][1]/rpl);
         rpl=ndownm[2][0]+ndownm[2][1];
         printf("  %s %5.0f  %4.3f %4.3f \n","ndown-mode 2 TOT no yes:",rpl,ndownm[2][0]/rpl,ndownm[2][1]/rpl);
         rpl=ndownm[3][0]+ndownm[3][1];
         printf("  %s %5.0f  %4.3f %4.3f \n","ndown-mode 3 TOT no yes:",rpl,ndownm[3][0]/rpl,ndownm[3][1]/rpl);
         printf("%s \n"," ");

         rpl=nupm[1][0]+nupm[1][1];
         printf("  %s %5.0f  %4.3f %4.3f \n","nup-mode 1 TOT no yes:",rpl,nupm[1][0]/rpl,nupm[1][1]/rpl);
         rpl=nupm[2][0]+nupm[2][1];
         printf("  %s %5.0f  %4.3f %4.3f \n","nup-mode 2 TOT no yes:",rpl,nupm[2][0]/rpl,nupm[2][1]/rpl);
         rpl=nupm[3][0]+nupm[3][1];
         printf("  %s %5.0f  %4.3f %4.3f \n","nup-mode 3 TOT no yes:",rpl,nupm[3][0]/rpl,nupm[3][1]/rpl);
         printf("%s \n"," ");

         rpl=nnodem[1][1]+nnodem[1][2]+nnodem[1][3];
         printf("  %s %5.0f  %4.3f %4.3f %4.3f \n","nnode-mode 1 TOT -1 0 1:",rpl,nnodem[1][1]/rpl,nnodem[1][2]/rpl,nnodem[1][3]/rpl);
         rpl=nnodem[2][1]+nnodem[2][2]+nnodem[2][3];
         printf("  %s %5.0f  %4.3f %4.3f %4.3f \n","nnode-mode 2 TOT -1 0 1:",rpl,nnodem[2][1]/rpl,nnodem[2][2]/rpl,nnodem[2][3]/rpl);
         rpl=nnodem[3][1]+nnodem[3][2]+nnodem[3][3];
         printf("  %s %5.0f  %4.3f %4.3f %4.3f \n","nnode-mode 3 TOT -1 0 1:",rpl,nnodem[3][1]/rpl,nnodem[3][2]/rpl,nnodem[3][3]/rpl);
         printf("%s \n"," ");

         rpl=nsurfm[1][0]+nsurfm[1][1];
         printf("   %s %5.0f  %4.3f %4.3f \n","nsurf TOT mode 1 no yes:",
		    rpl,nsurfm[1][0]/rpl,nsurfm[1][1]/rpl);
         rpl=nsurfm[2][0]+nsurfm[2][1];
         printf("   %s %5.0f  %4.3f %4.3f \n","nsurf TOT mode 2 no yes:",
		    rpl,nsurfm[2][0]/rpl,nsurfm[2][1]/rpl);
         rpl=nsurfm[3][0]+nsurfm[3][1];
         printf("   %s %5.0f  %4.3f %4.3f \n","nsurf TOT mode 3+ no yes:",
		    rpl,nsurfm[3][0]/rpl,nsurfm[3][1]/rpl);
         printf("%s \n"," ");

         rpl=nbottm[1][0]+nbottm[1][1];
         printf("   %s %5.0f  %4.3f %4.3f \n","nbott TOT mode 1 no yes:",
		    rpl,nbottm[1][0]/rpl,nbottm[1][1]/rpl);
         rpl=nbottm[2][0]+nbottm[2][1];
         printf("   %s %5.0f  %4.3f %4.3f \n","nbott TOT mode 2 no yes:",
		    rpl,nbottm[2][0]/rpl,nbottm[2][1]/rpl);
         rpl=nsurfm[3][0]+nbottm[3][1];
         printf("   %s %5.0f  %4.3f %4.3f \n","nbott TOT mode 3+ no yes:",
		    rpl,nbottm[3][0]/rpl,nbottm[3][1]/rpl);
         printf("%s \n"," ");
	 /**********************************************************************************print topo stats***/
         printf("%s \n","TOPO");
         rpl=nt[1]+nt[2]+nt[3]+nt[4]+nt[5]+nt[6];
         printf("%s \n","topo TOT 1-2 2-3 3-4 4-5 5-6");
         printf("%2.1f %2.1f %2.1f %2.1f %2.1f %2.1f \n",topostat[1],topostat[2],topostat[3],topostat[4],topostat[5],topostat[6]);
         printf("%s %5.0f  %4.3f %4.3f %4.3f %4.3f %4.3f %4.3f \n","ntopo TOT :",rpl,nt[1]/rpl,nt[2]/rpl,nt[3]/rpl,nt[4]/rpl,nt[5]/rpl,nt[6]/rpl);
         printf("%s \n"," ");

         nstandE[0][4]=nstandt[1][0]+nstandt[1][1];
         printf("  %s %5.0f  %4.3f %4.3f \n","nstand-t1 TOT no yes:",nstandE[0][4],nstandt[1][0]/nstandE[0][4],nstandt[1][1]/nstandE[0][4]);
         nstandE[0][4]=nstandt[2][0]+nstandt[2][1];
         printf("  %s %5.0f  %4.3f %4.3f \n","nstand-t2 TOT no yes:",nstandE[0][4],nstandt[2][0]/nstandE[0][4],nstandt[2][1]/nstandE[0][4]);
         nstandE[0][4]=nstandt[3][0]+nstandt[3][1];
         printf("  %s %5.0f  %4.3f %4.3f \n","nstand-t3 TOT no yes:",nstandE[0][4],nstandt[3][0]/nstandE[0][4],nstandt[3][1]/nstandE[0][4]);
         nstandE[0][4]=nstandt[4][0]+nstandt[4][1];
         printf("  %s %5.0f  %4.3f %4.3f \n","nstand-t4 TOT no yes:",nstandE[0][4],nstandt[4][0]/nstandE[0][4],nstandt[4][1]/nstandE[0][4]);
         nstandE[0][4]=nstandt[5][0]+nstandt[5][1];
         printf("  %s %5.0f  %4.3f %4.3f \n","nstand-t5 TOT no yes:",nstandE[0][4],nstandt[5][0]/nstandE[0][4],nstandt[5][1]/nstandE[0][4]);
         nstandE[0][4]=nstandt[6][0]+nstandt[6][1];
         printf("  %s %5.0f  %4.3f %4.3f \n","nstand-t6 TOT no yes:",nstandE[0][4],nstandt[6][0]/nstandE[0][4],nstandt[6][1]/nstandE[0][4]);
         printf("%s \n"," ");

         ndownE[0][4]=ndownt[1][0]+ndownt[1][1];
         printf("  %s %5.0f  %4.3f %4.3f \n","ndown-t1 TOT no yes:",ndownE[0][4],ndownt[1][0]/ndownE[0][4],ndownt[1][1]/ndownE[0][4]);
         ndownE[0][4]=ndownt[2][0]+ndownt[2][1];
         printf("  %s %5.0f  %4.3f %4.3f \n","ndown-t2 TOT no yes:",ndownE[0][4],ndownt[2][0]/ndownE[0][4],ndownt[2][1]/ndownE[0][4]);
         ndownE[0][4]=ndownt[3][0]+ndownt[3][1];
         printf("  %s %5.0f  %4.3f %4.3f \n","ndown-t3 TOT no yes:",ndownE[0][4],ndownt[3][0]/ndownE[0][4],ndownt[3][1]/ndownE[0][4]);
         ndownE[0][4]=ndownt[4][0]+ndownt[4][1];
         printf("  %s %5.0f  %4.3f %4.3f \n","ndown-t4 TOT no yes:",ndownE[0][4],ndownt[4][0]/ndownE[0][4],ndownt[4][1]/ndownE[0][4]);
         ndownE[0][4]=ndownt[5][0]+ndownt[5][1];
         printf("  %s %5.0f  %4.3f %4.3f \n","ndown-t5 TOT no yes:",ndownE[0][4],ndownt[5][0]/ndownE[0][4],ndownt[5][1]/ndownE[0][4]);
         ndownE[0][4]=ndownt[6][0]+ndownt[6][1];
         printf("  %s %5.0f  %4.3f %4.3f \n","ndown-t6 TOT no yes:",ndownE[0][4],ndownt[6][0]/ndownE[0][4],ndownt[6][1]/ndownE[0][4]);
         printf("%s \n"," ");

         nupE[0][4]=nupt[1][0]+nupt[1][1];
         printf("  %s %5.0f  %4.3f %4.3f \n","nup-t1 TOT no yes:",nupE[0][4],nupt[1][0]/nupE[0][4],nupt[1][1]/nupE[0][4]);
         nupE[0][4]=nupt[2][0]+nupt[2][1];
         printf("  %s %5.0f  %4.3f %4.3f \n","nup-t2 TOT no yes:",nupE[0][4],nupt[2][0]/nupE[0][4],nupt[2][1]/nupE[0][4]);
         nupE[0][4]=nupt[3][0]+nupt[3][1];
         printf("  %s %5.0f  %4.3f %4.3f \n","nup-t3 TOT no yes:",nupE[0][4],nupt[3][0]/nupE[0][4],nupt[3][1]/nupE[0][4]);
         nupE[0][4]=nupt[4][0]+nupt[4][1];
         printf("  %s %5.0f  %4.3f %4.3f \n","nup-t4 TOT no yes:",nupE[0][4],nupt[4][0]/nupE[0][4],nupt[4][1]/nupE[0][4]);
         nupE[0][4]=nupt[5][0]+nupt[5][1];
         printf("  %s %5.0f  %4.3f %4.3f \n","nup-t5 TOT no yes:",nupE[0][4],nupt[5][0]/nupE[0][4],nupt[5][1]/nupE[0][4]);
         nupE[0][4]=nupt[6][0]+nupt[6][1];
         printf("  %s %5.0f  %4.3f %4.3f \n","nup-t6 TOT no yes:",nupE[0][4],nupt[6][0]/nupE[0][4],nupt[6][1]/nupE[0][4]);
         printf("%s \n"," ");

         nmode[0]=nmodet[1][1]+nmodet[1][2]+nmodet[1][3];
         printf("   %s %5.0f  %4.3f %4.3f %4.3f \n","nmode-t1 TOT 1 2 >3:",nmode[0],nmodet[1][1]/nmode[0],nmodet[1][2]/nmode[0],nmodet[1][3]/nmode[0]);
         nmode[0]=nmodet[2][1]+nmodet[2][2]+nmodet[2][3];
         printf("   %s %5.0f  %4.3f %4.3f %4.3f \n","nmode-t2 TOT 1 2 >3:",nmode[0],nmodet[2][1]/nmode[0],nmodet[2][2]/nmode[0],nmodet[2][3]/nmode[0]);
         nmode[0]=nmodet[3][1]+nmodet[3][2]+nmodet[3][3];
         printf("   %s %5.0f  %4.3f %4.3f %4.3f \n","nmode-t3 TOT 1 2 >3:",nmode[0],nmodet[3][1]/nmode[0],nmodet[3][2]/nmode[0],nmodet[3][3]/nmode[0]);
         nmode[0]=nmodet[4][1]+nmodet[4][2]+nmodet[4][3];
         printf("   %s %5.0f  %4.3f %4.3f %4.3f \n","nmode-t4 TOT 1 2 >3:",nmode[0],nmodet[4][1]/nmode[0],nmodet[4][2]/nmode[0],nmodet[4][3]/nmode[0]);
         nmode[0]=nmodet[5][1]+nmodet[5][2]+nmodet[5][3];
         printf("   %s %5.0f  %4.3f %4.3f %4.3f \n","nmode-t5 TOT 1 2 >3:",nmode[0],nmodet[5][1]/nmode[0],nmodet[5][2]/nmode[0],nmodet[5][3]/nmode[0]);
         nmode[0]=nmodet[6][1]+nmodet[6][2]+nmodet[6][3];
         printf("   %s %5.0f  %4.3f %4.3f %4.3f \n","nmode-t6 TOT 1 2 >3:",nmode[0],nmodet[6][1]/nmode[0],nmodet[6][2]/nmode[0],nmodet[6][3]/nmode[0]);
         printf("%s \n"," ");

         nnode[0]=nnodet[1][1]+nnodet[1][2]+nnodet[1][3];
         printf("   %s %5.0f  %4.3f %4.3f %4.3f \n","nnode TOT-t1 -1 0 +1:",nnode[0],nnodet[1][1]/nnode[0],nnodet[1][2]/nnode[0],nnodet[1][3]/nnode[0]);
         nnode[0]=nnodet[2][1]+nnodet[2][2]+nnodet[2][3];
         printf("   %s %5.0f  %4.3f %4.3f %4.3f \n","nnode TOT-t2 -1 0 +1:",nnode[0],nnodet[2][1]/nnode[0],nnodet[2][2]/nnode[0],nnodet[2][3]/nnode[0]);
         nnode[0]=nnodet[3][1]+nnodet[3][2]+nnodet[3][3];
         printf("   %s %5.0f  %4.3f %4.3f %4.3f \n","nnode TOT-t3 -1 0 +1:",nnode[0],nnodet[3][1]/nnode[0],nnodet[3][2]/nnode[0],nnodet[3][3]/nnode[0]);
         nnode[0]=nnodet[4][1]+nnodet[4][2]+nnodet[4][3];
         printf("   %s %5.0f  %4.3f %4.3f %4.3f \n","nnode TOT-t4 -1 0 +1:",nnode[0],nnodet[4][1]/nnode[0],nnodet[4][2]/nnode[0],nnodet[4][3]/nnode[0]);
         nnode[0]=nnodet[5][1]+nnodet[5][2]+nnodet[5][3];
         printf("   %s %5.0f  %4.3f %4.3f %4.3f \n","nnode TOT-t5 -1 0 +1:",nnode[0],nnodet[5][1]/nnode[0],nnodet[5][2]/nnode[0],nnodet[5][3]/nnode[0]);
         nnode[0]=nnodet[6][1]+nnodet[6][2]+nnodet[6][3];
         printf("   %s %5.0f  %4.3f %4.3f %4.3f \n","nnode TOT-t6 -1 0 +1:",nnode[0],nnodet[6][1]/nnode[0],nnodet[6][2]/nnode[0],nnodet[6][3]/nnode[0]);
         printf("%s \n"," ");

         rpl=nsurft[1][0]+nsurft[1][1];
         printf("   %s %5.0f  %4.3f %4.3f \n","nsurf TOT topo 1 no yes:",
		    rpl,nsurft[1][0]/rpl,nsurft[1][1]/rpl);
         rpl=nsurft[2][0]+nsurft[2][1];
         printf("   %s %5.0f  %4.3f %4.3f \n","nsurf TOT topo 2 no yes:",
		    rpl,nsurft[2][0]/rpl,nsurft[2][1]/rpl);
         rpl=nsurft[3][0]+nsurft[3][1];
         printf("   %s %5.0f  %4.3f %4.3f \n","nsurf TOT topo 3 no yes:",
		    rpl,nsurft[3][0]/rpl,nsurft[3][1]/rpl);
         rpl=nsurft[4][0]+nsurft[4][1];
         printf("   %s %5.0f  %4.3f %4.3f \n","nsurf TOT topo 4 no yes:",
		    rpl,nsurft[4][0]/rpl,nsurft[4][1]/rpl);
         rpl=nsurft[5][0]+nsurft[5][1];
         printf("   %s %5.0f  %4.3f %4.3f \n","nsurf TOT topo 5 no yes:",
		    rpl,nsurft[5][0]/rpl,nsurft[5][1]/rpl);
         rpl=nsurft[6][0]+nsurft[6][1];
         printf("   %s %5.0f  %4.3f %4.3f \n","nsurf TOT topo 6 no yes:",
		    rpl,nsurft[6][0]/rpl,nsurft[6][1]/rpl);
         printf("%s \n"," ");

         rpl=nbottt[1][0]+nbottt[1][1];
         printf("   %s %5.0f  %4.3f %4.3f \n","nbott TOT topo 1 no yes:",
		    rpl,nbottt[1][0]/rpl,nbottt[1][1]/rpl);
         rpl=nbottt[2][0]+nbottt[2][1];
         printf("   %s %5.0f  %4.3f %4.3f \n","nbott TOT topo 2 no yes:",
		    rpl,nbottt[2][0]/rpl,nbottt[2][1]/rpl);
         rpl=nbottt[3][0]+nbottt[3][1];
         printf("   %s %5.0f  %4.3f %4.3f \n","nbott TOT topo 3 no yes:",
		    rpl,nbottt[3][0]/rpl,nbottt[3][1]/rpl);
         rpl=nbottt[4][0]+nbottt[4][1];
         printf("   %s %5.0f  %4.3f %4.3f \n","nbott TOT topo 4 no yes:",
		    rpl,nbottt[4][0]/rpl,nbottt[4][1]/rpl);
         rpl=nbottt[5][0]+nbottt[5][1];
         printf("   %s %5.0f  %4.3f %4.3f \n","nbott TOT topo 5 no yes:",
		    rpl,nbottt[5][0]/rpl,nbottt[5][1]/rpl);
         rpl=nbottt[6][0]+nbottt[6][1];
         printf("   %s %5.0f  %4.3f %4.3f \n","nbott TOT topo 6 no yes:",
		    rpl,nbottt[6][0]/rpl,nbottt[6][1]/rpl);
	 } /* end if (ch */
	} /* next ch */
       } /*END*/