subroutine varan(gk,wk,xk,yk,fk,kk,axy,bxy,fxy,xi,yj,
     +                 ixx,ix,jy,ibc,jbc,iax,iby,lulog,ierr)
c
c     This routine performs a variational analysis of the
c     observations gk to give fxy on an evenly spaced grid.
c     Observations are from 1 to k, and analysis grid points
c     are from 0 to ix and 0 to jy.  If gk is outside the
c     analysis grid, it is set to eval, and ignored in the
c     analysis.
c
c     Input:
c         kk    = No. of observations
c         gk    = Value of observation k
c         wk    = Analysis weight for observation k
c         xk    = x-coordinate of observation k
c         yk    = y-coordinate of observation k
c         ixx   = Maximum x-dimension of analysis grid arrays
c         ix    = No. of x analysis grid intervals
c         jy    = No. of y analysis grid intervals
c         axy   = Filter parameter in x-direction
c         bxy   = Filter parameter in y-direction
c         iax   = flag for treatment of x-filter parameters (see below**)
c         iby   = flag for treatment of y-filter parameters (see below**)
c         xi    = x-coordinates of analysis grid (must be evenly spaced)
c         yj    = y-coordinates of analysis grid (must be evenly spaced)
c         ibc,jbc = boundary condition type in x,y directions
c                   (0= zero 2nd derivative condition)
c                   (1= periodic condition, f(0)=f(ix) or f(jy) )
c         lulog   = unit number for writing log messages
c
c     Output:
c        fk    = Analysis grid values of fxy interpolated to
c                the observation points
c        fxy   = Analyzed function
c        ierr  = Error flag
c                0= Analysis completed normally
c                1= not enough observations inside analysis domain
c                2= Improper grid spacing
c                3= Iteration did not converge
c                4= Local variables not dimensioned large enough
c
c     **Further explanation of filter flag:
c     If iax=0 then it is assumed that axy are already calculated.
c
c     If iax .gt. 0 then, axy are calculated as follows. The
c     minimum allowable filter wavelength is for a half-power wavelength
c     of 2dx. If the data are sparse, axy are increased
c     using the average (iax=1) or minimum (iax=2) distance to the nearest data
c     at each grid point. 
c     The averaging or minimum computation is in y-direction.
c
c     If iax .lt. 0, then axy are calculated so that the half-power wavelength
c     of the filter is dx*abs(iax)/10. 
c
      dimension gk(kk),fk(kk),wk(kk),xk(kk),yk(kk)
      dimension axy(0:ixx,0:jy),bxy(0:ixx,0:jy),fxy(0:ixx,0:jy)
      dimension xi(0:ix),yj(0:jy)
c
c     Local variables
      parameter (lx=200,ly=200)
      dimension fxyt(-1:lx,-1:ly),dcdf(-1:lx,-1:ly)
      dimension idcell(0:lx,0:ly)
c
      parameter (lk=100000)
      dimension indxk(lk),indyk(lk),dk(lk)
      dimension wtint(lk,2,2)
c
      character *30 label
c
      common /cons/ eval,dx,dy
c
c     Set value for flagging observations outside of analysis domain
      eval = -999.9
c
c     Set percentage cost function decrease required for convergence
      pcmax = 0.000000001
c
c     Specify initial value of cost function parameter
      alpha = 0.1
c
c     Specify adjustment fraction for alpha
      aadj = 0.1
c
c     Specify maximum number of times to move down the gradient
c     without updating the gradient
      mmax=1
c
c     Specify maximum number of iterations
      mit=30000
c
c     Set itestp=0 for no array print outs on log file
c               =1 for data only print outs on log file
c               =2 for data and analyses print outs on log file
c     and set interval for printing information during iterations
      itestp=0
      modint=100
      modit = 25
c
c     Check dimensions of local work arrays
      if (kk .gt. lk) then
         write(lulog,969) kk,lk
 969     format('kk,lk: ',2i10)
         print*,kk,lk
	 write(lulog,970) 
  970    format(/,' Too many data points, increase lk in sub. varan')
	 ierr = 4
	 return
      endif
c
      if (ixx+1 .gt. lx .or. jy+1 .gt. ly) then
         print*,'ixx,lx,jy,ly : ',ixx,lx,jy,ly
	 write(lulog,975) 
  975    format(/,' Analysis grid too large,',
     +            ' increase lx and/or ly in sub. varan')
	 ierr = 4
	 return
      endif
c
c     Calculate grid spacing
      dx = xi(1)-xi(0)
      dy = yj(1)-yj(0)
      if (dx .le. 0.0 .or. dy .le. 0.0) then
         print*,dx,dy
         ierr=2
         return
      endif
c
c     Determine how many observation points are in the analysis grid,
c     and calculate grid parameters
      ko = 0
      do 10 k=1,kk
         if (gk(k) .ne. eval) then
            if (xk(k) .ge. xi(0) .and. xk(k) .le. xi(ix) .and.
     +          yk(k) .ge. yj(0) .and. yk(k) .le. yj(jy)      ) then
               ko = ko + 1
c
c              Calculate x,y indicies of closest analysis grid point
c              to the lower left of the current observation point
               indxk(k) = ifix( (xk(k)/dx) )
               indyk(k) = ifix( (yk(k)/dy) )
c
               if (indxk(k) .gt. ix-1) indxk(k) = ix-1
               if (indyk(k) .gt. jy-1) indyk(k) = jy-1
c
               indx = indxk(k)
               indy = indyk(k)
c
c              Calculate weights for interpolating grid values
c              to observation points.
               x1 = xi(indx  )
               x2 = xi(indx+1)
               y1 = yj(indy  )
               y2 = yj(indy+1)
c
               wtint(k,1,1) = (x2-xk(k))*(y2-yk(k))/(dx*dy)
               wtint(k,2,1) = (xk(k)-x1)*(y2-yk(k))/(dx*dy)
               wtint(k,1,2) = (x2-xk(k))*(yk(k)-y1)/(dx*dy)
               wtint(k,2,2) = (xk(k)-x1)*(yk(k)-y1)/(dx*dy)
            else
               gk(k) = eval
               indxk(k) = -1
               indyk(k) = -1
               wtint(k,1,1) = 0.0
               wtint(k,2,1) = 0.0
               wtint(k,1,2) = 0.0
               wtint(k,2,2) = 0.0
            endif
         endif
   10 continue
c
c     Identify grid cells with at least one data point in them
      do 11 j=0,jy
      do 11 i=0,ix
         idcell(i,j) = 0
   11 continue
c
      do 12 k=1,kk
         if (indxk(k) .lt. 0 .or. indyk(k) .lt. 0) go to 12
         ii = indxk(k)+1
         jj = indyk(k)+1
         idcell(ii,jj) = idcell(ii,jj) + 1
   12 continue
c
      if (itestp .ge. 1) then
         label='Grid Cells with Data'
         izero=1
         call iprint(idcell,lx,ix,jy,izero,lulog,label)
      endif
c
      if (ko .gt. 1) then
         write(lulog,200) ko,kk
  200    format(/,1x,i4,' out of ',i4,
     +          ' observations accepted in analysis.')
      else
         write(lulog,205)
  205    format(/,' Not enough observations in analysis domain.')
         ierr=1
         return
      endif
c
      if (itestp .ge. 1) then
c        Print data points
         label=' Input Data'
         call gprint(xk,yk,gk,kk,lulog,label)
      endif
c
c     Calculate filter parameters
      call abcal(axy,bxy,xi,yj,ixx,ix,jy,iax,iby)
c
      write(lulog,206) iax,iby
  206 format(/,' Filter parameters calculated with iax=',i4,' iby=',i4)
c  
c     Copy fxy to a temporary arrray
      do 15 j=0,jy
      do 15 i=0,ix
         fxyt(i,j) = fxy(i,j)
   15 continue
c
c     Apply boundary conditions to fxy
      call bcapp(fxyt,lx,ly,ix,jy,ibc,jbc)
c
      if (itestp .ge. 2) then
         iscale=0
         label=' fxy after bc applied'
         call eprint(fxyt,lx,ly,ix,jy,lulog,label,iscale)
      endif
c
c     Interpolate fxyt to observation points
      call xytoob(fxyt,lx,ly,fk,lk,kk,
     +            indxk,indyk,wtint)
c
      if (itestp .ge. 2) then
         label = ' First guess at ob points'
         call gprint(xk,yk,fk,kk,lulog,label)
      endif
c
c     Evaluate the initial value of the cost function
      call cfun(fk,gk,wk,kk,fxyt,axy,bxy,lx,ly,ixx,ix,jy,
     +          cdat,cpenx,cpeny,ctot)
c
      idum=0
      write(lulog,210) cdat,cpenx,cpeny,ctot,idum
  210 format(/,' cd=',e11.4,' cpx=',e11.4,' cpy=',e11.4,
     +         ' ctot=',e11.4,'  it=',i4)
c
c     Calculate the derivative of the cost function with
c     respect to fij.
      call gradc(fk,gk,wk,lk,kk,fxyt,axy,bxy,lx,ly,ixx,ix,jy,
     +           indxk,indyk,wtint,dcdf)
c
      if (itestp .ge. 2) then
         iscale=1
         label = 'Cost function gradient'
         call eprint(dcdf,lx,ly,ix,jy,lulog,label,iscale)
      endif
c
c     Start iterative solution for fxy
      cgnorm = dx*dy
c
c     Initialize old value of cost function for convergence test
      cfold = ctot
c
      do 98 m=1,mit
         if (mod(m,modit) .eq. 0) then
            write(lulog,215) m,alpha,ctot
  215       format(' Start iteration ',i4,' alpha=',f7.4,' ctot=',e12.5)
         endif
c
         ifinish=1
         gsign = -1.0
         do 18 mm=1,mmax
 1000       continue
c
c           Save previous value of cost function
            ctotp = ctot
c
c           Adjust fxyt using cost function gradient
            do 20 j=0,jy
            do 20 i=0,ix
               fxyt(i,j) = fxyt(i,j) + gsign*alpha*dcdf(i,j)/cgnorm
   20       continue
c 
c           Apply boundary conditions to fxy
            call bcapp(fxyt,lx,ly,ix,jy,ibc,jbc)
c
c           Interpolate fxyt to observation points
            call xytoob(fxyt,lx,ly,fk,lk,kk,
     +                  indxk,indyk,wtint)
c
c           Evaluate the cost function
            call cfun(fk,gk,wk,kk,fxyt,axy,bxy,lx,ly,ixx,ix,jy,
     +                cdat,cpenx,cpeny,ctot)
c
c            write(lulog,220) cdat,cpenx,cpeny,ctot,m,mm
c  220       format(' cd=',e11.4,' cpx=',e11.4,' cpy=',e11.4,
c     +             ' ctot=',e11.4,' it=',i4,' gs=',i2)
c
            if (ifinish .eq. 0) go to 1500
c
            if (ctot .gt. ctotp) then
c              Went too far down the gradient, undo the last adjustment
c              and move on to the next iteration. 
               ifinish = 0
               gsign = 1.0
               go to 1000
            endif
   18    continue
c
 1500    continue
         if (ifinish .eq. 1) then
c           The minimum was not reached during the last iteration.
c           Make alpha bigger for the next iteration.
            alpha = alpha*(1.0+aadj)
         else
c           The minimum was reached too soon during the last iteration.
c           Make alpha smaller for the next iteration.
            alpha = alpha*(1.0-aadj)
         endif
c
         if (mod(m,modint) .eq. 0 .and. itestp .ge. 2) then
            iscale=0
            label=' fxy after bc applied'
            call eprint(fxyt,lx,ly,ix,jy,lulog,label,iscale)
         endif
c
c        Calculate data deviations at observation points
         do 22 k=1,kk
            dk(k) = fk(k)-gk(k)
   22    continue
c
         if (mod(m,modint) .eq. 0 .and. itestp .ge. 1) then
            label = ' f-g at ob points'
            call gprint(xk,yk,dk,kk,lulog,label)
         endif
c
c        Calculate the derivative of the cost function with
c        respect to fij.
         call gradc(fk,gk,wk,lk,kk,fxyt,axy,bxy,lx,ly,ixx,ix,jy,
     +              indxk,indyk,wtint,dcdf)
C
         if (mod(m,modint) .eq. 0 .and. itestp .ge. 2) then
            iscale=1
            label = 'Cost function gradient'
            call eprint(dcdf,lx,ly,ix,jy,lulog,label,iscale)
         endif
c
c        Check for convergence
         if (ctot .eq. 0.0) go to 2000
c
         delc = (cfold-ctot)/cfold
c         write(6,410) cfold,ctot,delc,pcmax
c  410    format(' cfold,ctot,delc,pcmax: ',4(e11.4,1x))
         cfold=ctot
         if (delc .lt. pcmax .and. mm .gt. 1) go to 2000
   98 continue
c
      ierr=3
      write(lulog,225) m
  225 format(/,' Analysis did not converge after ',i6,' iterations.')
      go to 3000
c
 2000 continue
      write(lulog,230) m
  230 format(/' Analysis converged after ',i4,' iterations.')
c
c     Calculate the average difference between the data and the analysis
      dkrms = 0.0
      do k=1,kk
         dkrms = dkrms + dk(k)**2
      enddo
      dkrms = sqrt(dkrms/float(kk))
      write(lulog,236) dkrms,kk
  236 format('RMS of data fit: ',f8.2,' no. data points=',i6)
c
 3000 continue
c     Move fxyt back to fxy
      do 99 j=0,jy
      do 99 i=0,ix
         fxy(i,j) = fxyt(i,j)
   99 continue
c
      return
      end
      subroutine gradc(fk,gk,wk,lk,kk,fxyt,axy,bxy,lx,ly,ixx,ix,jy,
     +                 indxk,indyk,wtint,dcdf)
c     This routine calculates the gradient of the cost function
c     with respect to fxyt. 
c     
      dimension fk(kk),gk(kk),wk(kk)
      dimension fxyt(-1:lx,-1:ly),dcdf(-1:lx,-1:ly)
      dimension axy( 0:ixx,0:jy),bxy(0:ixx,0:jy)
      common /cons/ eval,dx,dy
c
      dimension indxk(lk),indyk(lk)
      dimension wtint(lk,2,2)
c
c     Initialize the gradient to zero
      do 10 j=-1,jy+1
      do 10 i=-1,ix+1
         dcdf(i,j) = 0.0
   10 continue
c
c     Calculate the contribution from the data misfit term
      do 20 k=1,kk
         if (indxk(k) .lt. 0) go to 20
         kx = indxk(k)
         ky = indyk(k)
         dd = wk(k)*(fk(k)-gk(k))*dx*dy
         dcdf(kx  ,ky  ) = dcdf(kx  ,ky  ) + dd*wtint(k,1,1)
         dcdf(kx+1,ky  ) = dcdf(kx+1,ky  ) + dd*wtint(k,2,1)
         dcdf(kx  ,ky+1) = dcdf(kx  ,ky+1) + dd*wtint(k,1,2)
         dcdf(kx+1,ky+1) = dcdf(kx+1,ky+1) + dd*wtint(k,2,2)
   20 continue
c
c     Calculate the contribution from the x-penalty term
      dxm4 = dx*dy*(dx**(-4))
      do 30 j=0,jy
         do 35 i=1,ix-1
            dcdf(i,j) = dcdf(i,j) + dxm4*(
     +      -2.0*axy(i  ,j)*(fxyt(i+1,j)+fxyt(i-1,j)-2.*fxyt(i  ,j))+
     +           axy(i-1,j)*(fxyt(i  ,j)+fxyt(i-2,j)-2.*fxyt(i-1,j))+
     +           axy(i+1,j)*(fxyt(i+2,j)+fxyt(i  ,j)-2.*fxyt(i+1,j)) )
   35    continue
c
         i=0
         dcdf(i,j) = dcdf(i,j) + dxm4*(
     +   -2.0*axy(i  ,j)*(fxyt(i+1,j)+fxyt(i-1,j)-2.*fxyt(i  ,j))+
     +        axy(i+1,j)*(fxyt(i+2,j)+fxyt(i  ,j)-2.*fxyt(i+1,j)) )
c
         i=ix
         dcdf(i,j) = dcdf(i,j) + dxm4*(
     +   -2.0*axy(i  ,j)*(fxyt(i+1,j)+fxyt(i-1,j)-2.*fxyt(i  ,j))+
     +        axy(i-1,j)*(fxyt(i  ,j)+fxyt(i-2,j)-2.*fxyt(i-1,j)) )
   30 continue
c
c     Calculate the contribution from the y-penalty term
      dym4 = dx*dy*(dy**(-4))
      do 40 i=0,ix
         do 45 j=1,jy-1
            dcdf(i,j) = dcdf(i,j) + dym4*(
     +      -2.0*bxy(i,  j)*(fxyt(i,j+1)+fxyt(i,j-1)-2.*fxyt(i,  j))+
     +           bxy(i,j-1)*(fxyt(i,j  )+fxyt(i,j-2)-2.*fxyt(i,j-1))+
     +           bxy(i,j+1)*(fxyt(i,j+2)+fxyt(i,j  )-2.*fxyt(i,j+1)) )
   45    continue
c
         j=0
         dcdf(i,j) = dcdf(i,j) + dym4*(
     +   -2.0*bxy(i,j  )*(fxyt(i,j+1)+fxyt(i,j-1)-2.*fxyt(i,  j))+
     +        bxy(i,j+1)*(fxyt(i,j+2)+fxyt(i,j  )-2.*fxyt(i,j+1)) )
c
         j=jy
         dcdf(i,j) = dcdf(i,j) + dym4*(
     +   -2.0*bxy(i,j  )*(fxyt(i,j+1)+fxyt(i,j-1)-2.*fxyt(i,  j))+
     +        bxy(i,j-1)*(fxyt(i,j  )+fxyt(i,j-2)-2.*fxyt(i,j-1)) )
   40 continue
c         
      return
      end
      subroutine cfun(fk,gk,wk,kk,fxy,axy,bxy,lx,ly,ixx,ix,jy,
     +                cdat,cpenx,cpeny,ctot)
c     This routine calculates the total cost function (ctot)
c     and the contributions from the data misfit term (cdat)
c     and the x and y contributions to the penalty term
c     (cpenx,cpeny) which imposes smoothness.
c
      dimension fk(kk),gk(kk),wk(kk)
      dimension fxy(-1:lx,-1:ly)
      dimension axy( 0:ixx,0:jy),bxy(0:ixx,0:jy)
      common /cons/ eval,dx,dy
c
c     Calculate the data misfit term
      cdat = 0.0
      do 10 k=1,kk
         if (gk(k) .ne. eval) then
            cdat = cdat + wk(k)*(fk(k)-gk(k))**2
         endif
   10 continue
      cdat = 0.5*cdat*dx*dy
c
c     Calculate the x,y penalty terms
      cpenx = 0.0
      cpeny = 0.0
      do 20 j=0,jy
      do 20 i=0,ix
         cpenx = cpenx + axy(i,j)*
     +                   ( fxy(i+1,j)+fxy(i-1,j)-2.0*fxy(i,j) )**2
         cpeny = cpeny + bxy(i,j)*
     +                   ( fxy(i,j+1)+fxy(i,j-1)-2.0*fxy(i,j) )**2
   20 continue
c
      cpenx = cpenx*dx*dy/(dx**4)
      cpeny = cpeny*dx*dy/(dy**4)
c
      ctot = cdat + cpenx + cpeny
c
      return
      end
      subroutine xytoob(fxy,lx,ly,fk,lk,kk,
     +                  indxk,indyk,wtint)
c
c     This routine bi-linearly interpolates fxy on the
c     evenly spaced grid to the observation points to give fk.
c
      dimension fxy(-1:lx,-1:ly)
      dimension fk(lk)
c
      dimension indxk(lk),indyk(lk)
      dimension wtint(lk,2,2)
c
      common /cons/ eval,dx,dy
c
      do 10 k=1,kk
         if (indxk(k) .lt. 0) then
            fk(k) = eval
         else
           indx = indxk(k)
           indy = indyk(k)
c
           fk(k) = wtint(k,1,1)*fxy(indx  ,indy  ) +
     +             wtint(k,2,1)*fxy(indx+1,indy  ) +
     +             wtint(k,1,2)*fxy(indx  ,indy+1) +
     +             wtint(k,2,2)*fxy(indx+1,indy+1)
         endif
   10 continue
c
      return
      end
      subroutine bcapp(fxy,lx,ly,ix,jy,ibc,jbc)
c     This routine applies the boundary conditions to the
c     array fxy. The boundary condition types are determined by
c     ibc as follows:
c          ibc = 0:  Second derivative condition
c              = 1:  Periodic condition
c              =-1:  fxy set to zero outside of domain
c
      dimension fxy(-1:lx,-1:ly)
      real avg
c
      if (ibc .eq. 0) then
c        Apply 2nd derivative condition in x-direction.
c
c        Fill in column of points outside domain
c        assuming d2f/dx2=0 at boundary
         avg=0.0
         do  9 j=0,jy
            avg=avg+fxy(0,j)
 9       continue
         avg=avg/(jy+1)
c         print*,'avg =',avg

         do 10 j=0,jy
            j8=j+jy/2
            if (j8.gt.jy)j8=j8-jy
ck            fxy(  -1,j) = avg
ck            fxy(   0,j) = avg
ck            fxy(ix+1,j) = fxy(ix-1,j)
ck            fxy(ix  ,j) = fxy(ix-1,j)
            fxy(  -1,j) = fxy(1,j8)
c            fxy(  -1,j) = 2.0*fxy( 0,j) - fxy(   1,j)
            fxy(ix+1,j) = 2.0*fxy(ix,j) - fxy(ix-1,j)
   10    continue
      endif
c
      if (jbc .eq. 0) then
c        Apply 2nd derivative condition in y-direction.
c        Fill in row of points outside domain
c        assuming d2f/dy2=0 at boundary.
         do 20 i=0,ix
            fxy(i,  -1) = 2.0*fxy(i, 0) - fxy(i,   1)
            fxy(i,jy+1) = 2.0*fxy(i,jy) - fxy(i,jy-1)
   20    continue
      endif
c
      if (ibc .eq. 1) then
c        Impose x-periodic conditions at right boundary and
c        fill in column of points outside the domain
         do 15 j=0,jy
            fxy(ix  ,j) = fxy(   0,j)
            fxy(ix+1,j) = fxy(   1,j)
            fxy(  -1,j) = fxy(ix-1,j)
   15    continue
c
c        Correction for mixed BC type (ibc=1,jbc=0)
         if (jbc .eq. 0) then
            fxy(ix  ,  -1) = fxy(0   ,  -1)
            fxy(ix  ,jy+1) = fxy(0   ,jy+1)
            fxy(  -1,  -1) = fxy(ix-1,  -1)
            fxy(  -1,jy+1) = fxy(ix-1,jy+1)
            fxy(ix+1,  -1) = fxy(   1,  -1)
            fxy(ix+1,jy+1) = fxy(   1,jy+1)
         endif
      endif
c
      if (jbc .eq. 1) then
c        Impose y-periodic conditions at top boundary and
c        fill in row of points outside the domain
         do 25 i=-1,ix+1
            fxy(i,jy  ) = fxy(i,   0)
            fxy(i,jy+1) = fxy(i,   1)
            fxy(i,  -1) = fxy(i,jy-1)
   25    continue
      endif
c
      return
      end
      subroutine abcal(axy,bxy,xi,yj,ixx,ix,jy,iax,iby)
c     This routine calculates the filter parameters axy,bxy.
c     The parameters iax,iby determine the method used to determine
c     axy,bxy.
c
      dimension axy(0:ixx,0:jy),bxy(0:ixx,0:jy)
      dimension xi(0:ix),yj(0:jy)
c
      dx = xi(1)-xi(0)
      dy = yj(1)-yj(0)
c
      pi = 3.14159
c
      if (iax .lt. 0) then
	 hppt = abs(0.1*float(iax))
c 
         do 10 j=0,jy
         do 10 i=0,ix
            axy(i,j) = ( (hppt*dx)/(2.0*pi) )**4
   10    continue
      endif
c
      if (iby .lt. 0) then
	 hppt = abs(0.1*float(iby))
c 
         do 20 j=0,jy
         do 20 i=0,ix
            bxy(i,j) = ( (hppt*dy)/(2.0*pi) )**4
   20    continue
      endif
c
      return
      end
      subroutine gprint(xk,yk,gk,kk,lulog,label)
c     This routine prints the value of gk at the observation
c     locations xk,yk.
c
      dimension xk(kk),yk(kk),gk(kk)
c
      character *30 label
c
      write(lulog,200) label
  200 format(/,a30,/,'  n     x    ', '   y   ','  g(x,y)  ')
c
      do 10 k=1,kk
         write(lulog,210) k,xk(k),yk(k),gk(k)
  210    format(1x,i4,1x,3(f7.2,1x))
   10 continue
c
      return
      end
      subroutine aprint(axy,ixx,ix,jy,lulog,luwin,label,iscale )
c     This routine prints the array axy. If iscale=1,
c     the array values scaled by the maximum array value
c     are printed. 
c
      dimension axy(0:ixx,0:jy)
      character *30 label
c
      if (iscale .eq. 1) then
c        Scale the data by the maximum value
         scale = 0.0
         do 10 j=0,jy
         do 10 i=0,ix
            atemp = axy(i,j)
            atemp = abs(atemp)
            if (atemp .gt. scale) scale = atemp
   10    continue
         if (scale .eq. 0.0) scale = 1.0
      else
         scale = 1.0
      endif
c
      write(lulog,100) label,scale
      write(luwin,100) label,scale
  100 format(/,1x,a30,'  Scale factor=',e10.3,/)
c
      do 20 i=0,ix
         write(lulog,110) i,(axy(i,j)/scale,j=0,jy)
	 write(luwin,110) i,(axy(i,j)/scale,j=0,jy)
  110    format(1x,i2,1x,21(f6.2))
   20 continue
c
      write(lulog,120) (j,j=0,jy)
      write(luwin,120) (j,j=0,jy)
  120 format(/,'  y:  ',21(i2,4x))
c
      return
      end
      subroutine iprint(iaxy,ixx,ix,jy,izero,lulog,label)
c     This routine prints the array integer axy.
c     If izero=0, then the printing starts at element
c     zero. Otherwise it starts at element 1. 
c
      dimension iaxy(0:ixx,0:jy)
      character *30 label
c
      if (izero .eq. 0) then
         istart=0
         jstart=0
      else
         istart=1
         jstart=1
      endif
c
      write(lulog,100) label
  100 format(/,1x,a30)
c
      do 20 j=jy,jstart,-1
         write(lulog,110) j,(iaxy(i,j),i=istart,ix)
  110    format(1x,i2,1x,20(i4))
   20 continue
c
      write(lulog,120) (i,i=istart,ix)
  120 format(/,6x,21(i2,2x))
c
      return
      end
      subroutine eprint(axy,lx,ly,ix,jy,lulog,label,iscale)
c     This routine prints the array axy (extended to include
c     one point beyond the analysis domain). If iscale=1,
c     the array values scaled by the maximum array value
c     are printed. 
c
      dimension axy(-1:lx,-1:ly)
      character *30 label
c
      if (iscale .eq. 1) then
c        Scale the data by the maximum value
         scale = 0.0
         do 10 j=-1,jy+1
         do 10 i=-1,ix+1
            atemp = axy(i,j)
            atemp = abs(atemp)
            if (atemp .gt. scale) scale = atemp
   10    continue
         if (scale .eq. 0.0) scale = 1.0
      else
         scale = 1.0
      endif
c
      write(lulog,100) label,scale
  100 format(/,1x,a30,'  Scale factor=',e10.3,/)
c
      do 20 j=jy+1,-1,-1
         write(lulog,110) j,(axy(i,j)/scale,i=-1,ix+1)
  110    format(1x,i2,1x,21(f6.2))
   20 continue
c
      write(lulog,120) (i,i=-1,ix+1)
  120 format(/,6x,21(i2,4x))
c
      return
      end
      subroutine xytora(x,y,r,a)
c     This routine calculates the radius and azimuth in degrees
c     (r,a) from x and y in a coordinate system where azimuth
c     is measured counterclockwise from the +x axis.
c
      rtd = 180.0/3.14159
c
      r = sqrt(x*x + y*y)
      if (r .eq. 0.0) then
         a = 0.0
         return
      endif
c
      a = rtd*acos(x/r)
      if (y .lt. 0.0) a = 360.0-a
c
      return
      end
      subroutine ratoxy(r,a,x,y)
c     This routine calculates x,y from the radius and azimuth 
c     in degrees (r,a) in a coordinate system where azimuth
c     is measured counterclockwise from the +x axis.
c
      dtr = 3.14159/180.0
c
      x = r*cos(dtr*a)
      y = r*sin(dtr*a)
c
      return
      end
      subroutine lltoxy(tlat,tlon,slat,slon,x,y)
c     This routine calculates the x,y coordinates in km of the point
c     tlat,tlon, assuming the origin is at slat,slon. All lat/lon
c     are input in degrees N and E.
c
      pi = 3.14159
      cy = 111.0
      cx = 111.0*cos(pi*slat/180.0)
c
      y = cy*(tlat-slat)
      x = cx*(tlon-slon)
c
      return
      end
      subroutine xytoll(x,y,slat,slon,tlat,tlon)
c     This routine calculates the lat/lon of the point xc,yc
c     assuming the origin is at slat,slon. The x,y input are in km 
c     and the tlat,tlon output are input in degrees N and E.
c
      pi = 3.14159
      cy = 111.0
      cx = 111.0*cos(pi*slat/180.0)
c
      tlat = slat + y/cy
      tlon = slon + x/cx
c
      return
      end