c  Program absolute velocity grid.  Microsoft Fortran 4.1, July 9, 1990
c  Ported to Sun Fortran, July 12, 1993.
c  Produces a grid of velocities.  Input is from standard binary files
c  via DOPIN.  Depth bins are specified ADCP bins, and variable depth
c  penetration is allowed for.  Only data above a specified acceptance
c  level in pct is used; 35% is recommended.
c  Oct 17, 1994:  Horizontal grid coordinate can be latitude, time, or
c  a special coordinate system.  These are found in subroutine special
c  at end of this file.  Routines written so far are:
c      a linear function of latitude,longitude, and time 
c      a clockwise angle from north relative to specified center point
c  Four variable options are possible:
c      vel - zonal and meridional velocities and number of profiles averaged
c            are output in files u..., v..., n...
c      qua - percent good and Automatic Gain Control values are output to
c            p..., a...
c      tem - temperature at instrument depth is output to t... 
c      ver - vertical velocity and number of profiles averaged are output to 
c            w..., n...  [Note that vertical velocity is generally noise]
      parameter (maxlat=2000,maxdep=50,max=maxlat*maxdep)
      real v(maxdep,maxlat),u(maxdep,maxlat),ave(maxdep,maxlat)
      real*8 time(maxlat),pos1(maxlat),pos2(maxlat)
      real*8 ttime,tpos1,tpos2,tzero
      integer npos(maxlat)
      integer ngrid(maxdep,maxlat)
      real r(5)
      integer inp(5)
      real*8 dinp(5),dbl,time1,time2,cengrd,delta,delta2
      character name*50,ext*3,var*3,type*3,formt*50
      common /arrayu/u
      common /arrayv/v
      common /arrayn/n
      common /fun$e/iin,iout
      common /fun$std/iinstd,ioutstd

c  Special coordinate variables shared between SPECIAL and ABSVGRD
      real spec1,spec2
      common /speclim/cengrd,delta2,spec1,spec2


c  Header parameters from DOPIN
c  ...transducer characteristics
      character*2 trans
c  ...profiler setup commands from AD menu
      integer*2 epings
      byte nbins,binwid,pulwid,blank,pcntmin
      real*4 avgint
c  ...date and cruise ID
      byte iyr,imo,idy,ihr,imn,isc,ihd
      character*60 cruise
      character*20 notes
      common /tranc2m/trans
      common /setupcom/avgint,epings,nbins,binwid,pulwid,blank,pcntmin
      common /datecom/iyr,imo,idy,ihr,imn,isc,ihd,cruise,notes

c  Data parameters from DOPIN
c  ...Raw date, corrected date, # of pings, recent headingX100, ave tempX1000
      real*8 rdate,cdate
      integer*2 npings,irechead,iavetemp
c  ...Velocity data (mm/s)
      parameter (vmax=64)
      integer*2 vx(vmax),vy(vmax),vz(vmax),ve(vmax)
      byte agc(vmax),pct(vmax)
c  ...Bottom tracking and navigation data
      integer*2 vxbot,vybot,depth,vxship,vyship
      real*8 lat,long
c  ...Calibration factorX10000 and heading biasX1000
      integer*2 calib,headbias
      common /cdatecom/rdate,cdate,npings,irechead,iavetemp
      common /velcom/vx,vy,vz,ve,agc,pct
      common /trackcom/lat,long,vxbot,vybot,depth,vxship,vyship
      common /calibcom/calib,headbias
 

      idsk=10
      undef=999.9
      call ioask

c  Get variable and grid type
    5 call aenter(var,3,' Variable desired (vel,qua,tem,ver): ')
      if(var.ne.'vel'.and.var.ne.'qua'.and.
     +   var.ne.'tem'.and.var.ne.'ver')go to 500
      call aenter(type,3,' Enter grid desired (lat,tim,spc): ')
      if(type.ne.'lat'.and.type.ne.'tim'.and.type.ne.'spc')go to 500

c  Get grid parameters
c  ...In latitude
      call renter(r,2,' Enter S and N latitude limits (deg): ')
      if(r(1).ge.r(2))then
         alat1=-90.0
         alat2=+90.0
        else
         alat1=r(1)
         alat2=r(2)
        endif
      call denter(dinp,2,
     +   ' Enter time bracket (decimal days, 0s for none): ')
      if(dinp(1).ge.dinp(2))then
         time1=-1.0e35
         time2=+1.0e35
        else
         time1=dinp(1)
         time2=dinp(2)
        endif
      call denter(tzero,1,' Zero time: ')
      if(type.eq.'spc')then
c    ...get special coordinates set up
        iflag=1
        call special(iflag,lat,long,cdate)
        endif
      call denter(delta,1,' Enter bin width: ')
      if(type.eq.'lat') a=1.0+(alat2-alat1)*1.0001/delta
      if(type.eq.'tim') a=1.0+(time2-time1)*1.0001/delta
      if(type.eq.'spc') a=1.0+(spec2-spec1)*1.0001/delta
      ngrd=nint(a)
c  ...subdivide grid by ndiv to give better averaging
      call ienter(ndiv,1,' Subgrid division factor: ')
      ngrd2=ngrd*ndiv
      delta2=delta/float(ndiv)
c     ...cengrd is center of first sub-bin
      if(type.eq.'lat')then
        cengrd=alat1-delta2*0.5*(ndiv-1)
        a=alat1
        b=alat2
        endif
      if(type.eq.'tim')then
        cengrd=time1-delta2*0.5*(ndiv-1)
        a=time1
        b=time2
        endif
      if(type.eq.'spc')then
        cengrd=spec1-delta2*0.5*(ndiv-1)
        a=spec1
        b=spec2
        endif
      write(iout,1020)a,b,delta,ngrd,ndiv,delta2,ngrd2,cengrd
 1020 format(' Grid spans',2f9.4,/,
     &       6x,' bin width',f9.4,' gives',i4,' grid bins',/
     &       6x,i4,' subdiv gives',f9.4,' sub-width and',i4,
     &       ' sub-bins, first centered at',f9.4)
c  ...In depth
      call ienter(inp,3,' Enter top, bottom depth bin and increment: ')
      minbin=inp(1)
      maxbin=inp(2)
      ndelbin=inp(3)
      ndep=1+(maxbin-minbin)/ndelbin
      if(var.eq.'tem')ndep=1
      maxbin=minbin+(ndep-1)*ndelbin
      write(iout,1010)minbin,maxbin,ndelbin
 1010 format(' Using bin range',2i3,' by',i2)
      i=ngrd2*ndep
      if(i.gt.max)stop ' ***** Grid too big *****'

c  Get acceptance threshholds; note that ranges straddling dateline
c  are included in possiblities.
      call renter(r,2,' Enter W and E longitude limits: ')
      along=r(1)
      blong=r(2)
      call ienter(minpct,1,' Enter minimum ADCP percent in %: ')
      call ienter(minave,1,' Enter minimum valid subgrid entries: ')

c  Zero arrays
   10 do 30 j=1,maxlat
        do 20 i=1,maxdep
          v(i,j)=0.0
          u(i,j)=0.0
          ngrid(i,j)=0
          ave(i,j)=0.0
   20     continue
        time(j)=0.0
        pos1(j)=0.0
        pos2(j)=0.0
        npos(j)=0
   30   continue
      icnt=0
      jcnt=0

c  Get ADCP input file template
      call aenter(name,50,' Enter input filename (with extension): ')
      nc=lenstr(name)
      if(nc.eq.0)go to 5
      ext=name(nc-2:nc)
      read(ext,1000)jext

C  Prepare to open new ADCP file: signal initial read with irec=0
  100 write(ext,1000)jext
 1000 format(i3)
      call zerblk(ext)
      name(nc-2:nc)=ext
      write(iout,1100)name
      if(iout.ne.ioutstd)write(*,1100)name
 1100 format(' opening input ADCP file ',a50)
      jext=jext+1
      irec=0

c  Read data
c  ...Get new record
  200 call dopin(irec,name,idsk)
      if(irec.eq.-2)go to 300
      if(irec.le.0)then
        write(iout,2000)
 2000   format(' End of file')
        close(idsk)
        go to 100
        endif
      icnt=icnt+1
c  ...Check acceptance criterion
      if(type.eq.'lat')then
         if(cdate.lt.time1)go to 200
         if(cdate.gt.time2)go to 300
        elseif(type.eq.'tim')then
         if(lat.lt.alat1)go to 200
         if(lat.gt.alat2)go to 200
         if(cdate.gt.time2)go to 300
        elseif(type.eq.'spc')then
         if(cdate.lt.time1)go to 200
         if(cdate.gt.time2)go to 300
         if(lat.lt.alat1)go to 200
         if(lat.gt.alat2)go to 200
        else
         stop
         endif
      if((var.eq.'vel').and.(vxship.gt.29000))go to 200
c     ...Longitude test also valid across dateline
	      a=(blong-long)*(long-along)/(blong-along)
      if(a.lt.0.0)go to 200
      jcnt=jcnt+1

c  Process record
c  ...Find grid bin
      if(type.eq.'lat') igrd=1.5+(lat-cengrd)/delta2
      if(type.eq.'tim') igrd=1.5+(cdate-cengrd)/delta2
      if(type.eq.'spc')then
        iflag=2
        call special(iflag,lat,long,cdate)
        igrd=iflag
        endif
      if(igrd.lt.1.or.igrd.gt.ngrd2)go to 200

ctest
ctest      write(iout,*)' averaging',cdate
ctest

c  ...Load horizontal velocities into grid in cm/s
      if(var.eq.'vel')then
        do 210 i=1,ndep
          ibin=minbin+ndelbin*(i-1)
c      ...Test quality
          n=pct(ibin)
          if(n.lt.minpct)go to 210
c      ...Load data
          u(i,igrd)=u(i,igrd)+0.1*vx(ibin)
          v(i,igrd)=v(i,igrd)+0.1*vy(ibin)
          ngrid(i,igrd)=ngrid(i,igrd)+1
  210     continue
        go to 250
        endif

c  ...Or load qualities into grid in origional percent and gain
      if(var.eq.'qua')then
        do 220 i=1,ndep
          ibin=minbin+ndelbin*(i-1)
c      ...Load data
          u(i,igrd)=u(i,igrd)+float(pct(ibin))
          a=float(agc(ibin))
          if(a.lt.0.0)a=a+256.0
          v(i,igrd)=v(i,igrd)+a
          ngrid(i,igrd)=ngrid(i,igrd)+1
  220     continue
        go to 250
        endif

c  ...Or load temperatures into top of grid in degrees C
      if(var.eq.'tem')then
c    ...Load data
        u(1,igrd)=u(1,igrd)+0.001*float(iavetemp)
        ngrid(1,igrd)=ngrid(1,igrd)+1
        go to 250
        endif

c  ...Or load vertical velocity into grid in cm/s
      if(var.eq.'ver')then
        do 230 i=1,ndep
          ibin=minbin+ndelbin*(i-1)
c      ...Test quality
          n=pct(ibin)
          if(n.lt.minpct)go to 230
c      ...Load data
          u(i,igrd)=u(i,igrd)+0.1*vz(ibin)
          ngrid(i,igrd)=ngrid(i,igrd)+1
  230     continue
        go to 250
        endif

c  Load time and position into labling arrays
  250 time(igrd)=time(igrd)+cdate-tzero
      pos1(igrd)=pos1(igrd)+lat
c  ...Longitude wrap-around put at 0 degrees for ranges
c  ...stradling dateline
      if((along.gt.blong).and.(long.lt.0.0))long=long+360.0
      pos2(igrd)=pos2(igrd)+long
      npos(igrd)=npos(igrd)+1

c  Go for new record
      go to 200


c  All data read; finish up
  300 continue

c  Average subgrid into grid.  i is grid latbin, j is depth, k is
c  subgrid latbin.  Divide each subgrid by number of data.
c  ...For every grid bin...
      do 360 i=1,ngrd
        k2=i*ndiv
        k1=k2-ndiv+1
c    ...And every depth...
        do 340 j=1,ndep
          n=0
          utemp=0.0
          vtemp=0.0
c      ...Average all the subgrids/ndata.
          do 330 k=k1,k2
            a=ngrid(j,k)
            if(a.le.0.0)go to 330
            utemp=utemp+u(j,k)/a
            vtemp=vtemp+v(j,k)/a
            n=n+1
  330       continue
c      ...Set grid to undefined if not enough valid subgrids
          ave(j,i)=n
          if(n.ge.minave)then
             u(j,i)=utemp/n
             v(j,i)=vtemp/n
            else
             u(j,i)=undef
             v(j,i)=undef
             endif
  340     continue
c    ...Also average extra labling arrays over subgrids
        n=0
        ttime=0.0
        tpos1=0.0
        tpos2=0.0
        do 350 k=k1,k2
          a=npos(k)
          if(a.le.0.0)go to 350
          ttime=ttime+time(k)/a
          tpos1=tpos1+pos1(k)/a
          tpos2=tpos2+pos2(k)/a
          n=n+1
  350     continue
c    ...Set lable arrays to undefined if not enough valid subgrids
        if(n.ge.minave)then
           a=n
           time(i)=ttime/a
           pos1(i)=tpos1/a
           pos2(i)=tpos2/a
c       ...Put longitude wrap-around back to 180 degrees
           if(pos2(i).gt.180.0)pos2(i)=pos2(i)-360.0
          else
           time(i)=undef
           pos1(i)=undef
           pos2(i)=undef
           endif
  360   continue

c  Output final files
      write(iout,*)' Output file name, first character to be'
      call aenter(name,50,' replaced by u,v,n or p,a or t: ')
      call aenter(formt,50,' Output format: ')
      if(type.eq.'lat') dbl=alat1
      if(type.eq.'tim') dbl=time1-tzero
      if(type.eq.'spc') dbl=spec1
      if(var.eq.'vel')then
        name(1:1)='u'
        call output(name,maxdep,ndep,ngrd,u,ave,dbl,delta,
     +              time,pos1,pos2,idsk,formt)
        name(1:1)='v'
        call output(name,maxdep,ndep,ngrd,v,ave,dbl,delta,
     +              time,pos1,pos2,idsk,formt)
        name(1:1)='n'
        call output(name,maxdep,ndep,ngrd,ave,ave,dbl,delta,
     +              time,pos1,pos2,idsk,formt)
        endif
      if(var.eq.'qua')then
        name(1:1)='p'
        call output(name,maxdep,ndep,ngrd,u,ave,dbl,delta,
     +              time,pos1,pos2,idsk,formt)
        name(1:1)='a'
        call output(name,maxdep,ndep,ngrd,v,ave,dbl,delta,
     +              time,pos1,pos2,idsk,formt)
        name(1:1)='n'
        call output(name,maxdep,ndep,ngrd,ave,ave,dbl,delta,
     +              time,pos1,pos2,idsk,formt)
        endif
      if(var.eq.'tem')then
        name(1:1)='t'
        call output(name,maxdep,ndep,ngrd,u,ave,dbl,delta,
     +              time,pos1,pos2,idsk,formt)
        endif
      if(var.eq.'ver')then
        name(1:1)='w'
        call output(name,maxdep,ndep,ngrd,u,ave,dbl,delta,
     +              time,pos1,pos2,idsk,formt)
        name(1:1)='n'
        call output(name,maxdep,ndep,ngrd,ave,ave,dbl,delta,
     +              time,pos1,pos2,idsk,formt)
        endif
      write(iout,3000)icnt,jcnt
 3000 format(i8,' Records read,',i8,' records averaged')
      go to 10

c  Finish
  500 end

c  ***************************************************************************

      subroutine output(name,maxdep,ndep,ngrd,data,ave,
     +                  dbl,delta,time,pos1,pos2,idsk,formt)
c  Write out data in ADCP format.
      real data(maxdep,ngrd),ave(maxdep,ngrd)
      real*8 time(ngrd),pos1(ngrd),pos2(ngrd)
      real*8 dbl,delta
      character name*50,formt*50
      call file2(i,2,name,idsk,' ',0)
      if(i.ne.0)go to 100
      write(idsk,1050)ndep,ngrd
 1050 format('ndep=',i2,' ngrd=',i3)
      do 10 j=1,ngrd
        a=dbl+float(j-1)*delta
        write(idsk,formt)a,(data(i,j),i=1,ndep),ave(1,j),
     +                    time(j),pos1(j),pos2(j)
c 1060   format(f9.4,10f9.3,6(/,9x,10f9.3))
   10   continue
      close (idsk)
  100 return
      end


c  *******************************************************************

      subroutine xspecial(iflag,lat,long,cdate)
c  Gets special coordinates for program ABSVGRD.
c  On entry:
c    iflag = 1 sets up special coordinates
c    iflag = 2 iflag set to coordinate bin of passed lat,long,cdate
c  This version produces a combined coordinate which is a linear
c  combination of time, latitude, and longitude.

      real*8 lat,long,cdate
      real r(2)
      real*8 dinp(4),dbl
      common /fun$e/iin,iout
      common /fun$std/iinstd,ioutstd

c  Special coordinate variables shared between SPECIAL and ABSVGRD
      real spec1,spec2
      real*8 cengrd,delta2
      common /speclim/cengrd,delta2,spec1,spec2
c  Special coordinate variables saved internally
      real*8 cmean,clat,clong,ctime
      common /specvar/cmean,clat,clong,ctime
 
c  iflag=1: set up coordinate coefficients and limits
      if(iflag.eq.1)then 
        write(iout,*)' Enter combined coordinate coefficients'
        call denter(dinp,4,' mean,lat(degs),long(degs),time(days): ')
        cmean=dinp(1)
        clat=dinp(2)
        clong=dinp(3)
        ctime=dinp(4)
        write(iout,9010)cmean,clat,clong,ctime
 9010   format(4f15.6)
        call renter(r,2,' Enter combined coordinate limits: ')
        spec1=r(1)
        spec2=r(2)
        go to 400
        endif

c  iflag=2: calculate bin from passed lat,long,cdate
      if(iflag.eq.2)then
        dbl= cmean + clat*lat + clong*long + ctime*cdate
        igrd=1.5+(dbl-cengrd)/delta2
        write(iout,9000)lat,long,cdate,dbl,igrd
 9000   format(4f13.6,i12)
        iflag=igrd
        go to 400
        endif
 
c  Bad iflag value
      write(*,2000)iflag
 2000 format(' Unrecognised flag value: ',i8)
      stop

c  Return
  400 return
      end

c  *******************************************************************

      subroutine special(iflag,lat,long,cdate)
c  Gets special coordinates for program ABSVGRD.
c  On entry:
c    iflag = 1 sets up special coordinates
c    iflag = 2 iflag set to coordinate bin of passed lat,long,cdate
c  This version produces a location on the Moving Ship Tomography
c  cruise circle: angle clockwise in degrees from due north, with
c  successive passes giving angles greater than 360 degrees.

      real*8 lat,long,cdate
      real r(2)
      real*8 dinp(4)
      common /fun$e/iin,iout
      common /fun$std/iinstd,ioutstd

c  Special coordinate variables shared between SPECIAL and ABSVGRD
      real spec1,spec2
      real*8 cengrd,delta2
      common /speclim/cengrd,delta2,spec1,spec2
c  Special coordinate variables saved internally
      real*8 zlat,zlong,wrap1,wrap2
      common /specvar/zlat,zlong,wrap1,wrap2
 
c  iflag=1: set up coordinate coefficients and limits
      if(iflag.eq.1)then 
        call denter(dinp,2,' Enter center lat,long in degrees: ')
        zlat=dinp(1)
        zlong=dinp(2)
        zfact=cos(0.0174533*zlat)
        call denter(dinp,2,' Enter two wrap-around times: ')
        wrap1=dinp(1)
        wrap2=dinp(2)
        write(iout,9010)zlat,zlong,wrap1,wrap2
 9010   format(4f15.6)
        call renter(r,2,' Enter angle limits: ')
        spec1=r(1)
        spec2=r(2)
        wrap=0.0
        go to 400
        endif

c  iflag=2: calculate bin from passed lat,long,cdate
      if(iflag.eq.2)then
c    ...Find actual angle, compensating for narrowed longitudes
        sine=zfact*(long-zlong)
        cosine=lat-zlat
        a=57.29578*atan2(sine,cosine)
c    ...Check for update of wrap factor: function wrap is at 180 deg
        if(wrap.eq.0.0.and.cdate.gt.wrap1.and.a.lt.140.0)wrap=360.0
        if(wrap.eq.360.0.and.cdate.gt.wrap2.and.a.lt.140.0)wrap=720.0
        a=a+wrap
  100   igrd=1.5+(a-cengrd)/delta2
ctest        write(iout,9000)lat,long,cdate,a,igrd
ctest 9000   format(4f13.6,i12)
        iflag=igrd
        go to 400
        endif
 
c  Bad iflag value
      write(*,2000)iflag
 2000 format(' Unrecognised flag value: ',i8)
      stop

c  Return
  400 return
      end