subroutine tpwpcal(tpw,mx,my,nx,ny,
     +                   tlon1,tlon2,dlon,tlat1,tlat2,dlat,
     +                   slat,slon,cdir00,sdir00,itpw,itmax)
c
c     This routine calculates TPW variables centered at slat,slon 
c
c     This version of the routine calculates area averaged and ring average
c     TPW predictors for the SHIPS model. The TPW fields and related grid parameters
c     are provided through the calling parameters. 
c
c     Modified Feb 2023 (KM) for NCO specifications
c
c     Input variables:
c       tpw(mx,my) - TPW values (mm) on an evenly spaced lat/lon grid
c       mx,my       - Max dimensions of tpwxy
c       nx,ny       - No. of lon,lat points with TPW data
c       tlat1,tlat2 - Bottom and top latitudes with TPW data (deg N)
c       dlat        - lat increment (deg)
c       tlon1,tlon2 - Left and right longitudes (0 to 360 deg convention)
c       dlon        - Lon increment (deg)
c       slat        - Storm lat (deg N)
c       slon        - Storm lon (0-360 deg)    
c       cdir00      - Direction of storm motion in deg CCW from +x axis (east)
c       sdir00      - Direction of shear vector in deg CCW from +x axis (east)
c       itmax       - Dimension of vector of output variables (itpw)
c
c     Output:
c       itpw(0:itmax) - Vector of 21 TPW variables as described below.
c                       These are stormed in itpw(i) for i=0 to 20
c       
c     1-2 0-200 km average & standard deviation
c     3-4 200-400 km average & standard deviation
c     5-6 400-600 km average & standard deviation
c     7-8 600-800 km average & standard deviation
c     9-10 800-1000 km average & standard deviation
c     11-12 0-400 km average & standard deviation
c     13-14 0-600 km average & standard deviation
c     15-16 0-800 km average & standard deviation
c     17-18 0-1000 km average & standard deviation
c     19. %TPW less than 45 mm, r=0 to r=500 km, 90 deg quadrant 
c     centered upshear
c     20. TPW averaged r=0 to r=500 km, 90 deg quadrant centered upshear
c     21. TPW averaged r=0 to r=500 km
c
c****************************************************************************
c     
c     Passed variables
      dimension tpw(mx,my)
      dimension itpw(0:itmax)
c
c     Local variables
c
c     Cylindrical grid variables
      parameter (mr=101,mt=16)
      dimension r(0:mr),theta(0:mt)
      dimension xrt(0:mr,0:mt),yrt(0:mr,0:mt)
      dimension tpwrt(0:mr,0:mt)
c
c     write(6,*) 'mx,my:',mx,my,' nx,ny:',nx,ny
c     write(6,*) 'tlat1,tlat2,dlat: ',tlat1,tlat2,dlat
c     write(6,*) 'tlon1,tlon2,dlon: ',tlon1,tlon2,dlat
c     write(6,*) 'slon,slat: ',slon,slat
c     Calculate cylindrical grid info
      nr=60
      nt=16
      r0=0.0
      dr=20.0
      dt=360.0/float(nt)
      dtr=3.14159/180.0
c
      do ii=0,nr
         r(ii) = r0 + dr*float(ii)
      enddo
c
      do jj=0,nt
         theta(jj) = dt*float(jj)
      enddo
c
      do ii=0,nr
      do jj=0,nt
         xrt(ii,jj) = r(ii)*cos(dtr*theta(jj))
         yrt(ii,jj) = r(ii)*sin(dtr*theta(jj))
      enddo
      enddo
c
c     Initialize tpw variables to missing
      do k=0,itmax
         itpw(k) = 9999
      enddo
c
c     Calcualte dx,dy of lat/lon grid at storm center location
      dy = 111.1*dlat
      dx = 111.1*dlon*cos(dtr*slat)
c
c     Assuming storm center is orgin, calculate x,y of lower left
c     TPW grid point
      xg1 = dx*(tlon1-slon)/dlon
      yg1 = dy*(tlat1-slat)/dlat
c
c     Interpolate TPW field to cylindrical grid
      tmiss = 0.1
      do ii=0,nr
      do jj=0,nt
         call lintcft(tpw,tmiss,xg1,dx,yg1,dy,mx,my,nx,ny,
     +               xrt(ii,jj),yrt(ii,jj),fxyii)
         tpwrt(ii,jj) = fxyii
      enddo
      enddo
c     
c     Calculate the symmetric TPW averages and standard deviations
c
c     ++itpw(0),tpw(1)
      r1 = 0.0
      r2 = 200.0
      th1 = 0.0
      th2 = 360.0

      call caavgsd(tpwrt,tmiss,r1,r2,th1,th2,mr,mt,nr,nt,r,theta,tbar
     .     ,tsd)
c
      if (tbar .gt. 0.0) then
         itpw(0) = nint(10.0*tbar)
         itpw(1) = nint(10.0*tsd)
      endif
c
c     ++itpw(2),tpw(3)
      r1  = 200.0
      r2  = 400.0
      th1 = 0.0
      th2 = 360.0

      call caavgsd(tpwrt,tmiss,r1,r2,th1,th2,mr,mt,nr,nt,r,theta,tbar
     .     ,tsd)

c
      if (tbar .gt. 0.0) then
         itpw(2) = ifix(10.0*tbar + 0.5)
         itpw(3) = ifix(10.0*tsd  + 0.5)
      endif

c     ++itpw(4),tpw(5)
      r1  = 400.0
      r2  = 600.0
      th1 = 0.0
      th2 = 360.0

      call caavgsd(tpwrt,tmiss,r1,r2,th1,th2,mr,mt,nr,nt,r,theta,tbar
     .     ,tsd)

c
      if (tbar .gt. 0.0) then
         itpw(4) = ifix(10.0*tbar + 0.5)
         itpw(5) = ifix(10.0*tsd  + 0.5)
      endif

c     ++itpw(6),tpw(7)
      r1  = 600.0
      r2  = 800.0
      th1 = 0.0
      th2 = 360.0

      call caavgsd(tpwrt,tmiss,r1,r2,th1,th2,mr,mt,nr,nt,r,theta,tbar
     .     ,tsd)

c
      if (tbar .gt. 0.0) then
         itpw(6) = ifix(10.0*tbar + 0.5)
         itpw(7) = ifix(10.0*tsd  + 0.5)
      endif
c     ++itpw(8),tpw(9)
      r1  = 800.0
      r2  = 1000.0
      th1 = 0.0
      th2 = 360.0

      call caavgsd(tpwrt,tmiss,r1,r2,th1,th2,mr,mt,nr,nt,r,theta,tbar
     .     ,tsd)
c
      if (tbar .gt. 0.0) then
         itpw(8) = ifix(10.0*tbar + 0.5)
         itpw(9) = ifix(10.0*tsd  + 0.5)
      endif

c     ++itpw(10),tpw(11)
      r1  = 0.0
      r2  = 400.0
      th1 = 0.0
      th2 = 360.0

      call caavgsd(tpwrt,tmiss,r1,r2,th1,th2,mr,mt,nr,nt,r,theta,tbar
     .     ,tsd)
c
      if (tbar .gt. 0.0) then
         itpw(10) = ifix(10.0*tbar + 0.5)
         itpw(11) = ifix(10.0*tsd  + 0.5)
      endif
c     ++itpw(12),tpw(13)
      r1  = 0.0
      r2  = 600.0
      th1 = 0.0
      th2 = 360.0

      call caavgsd(tpwrt,tmiss,r1,r2,th1,th2,mr,mt,nr,nt,r,theta,tbar
     .     ,tsd)
c
      if (tbar .gt. 0.0) then
         itpw(12) = ifix(10.0*tbar + 0.5)
         itpw(13) = ifix(10.0*tsd  + 0.5)
      endif

c     ++itpw(14),tpw(15)
      r1  = 0.0
      r2  = 800.0
      th1 = 0.0
      th2 = 360.0

      call caavgsd(tpwrt,tmiss,r1,r2,th1,th2,mr,mt,nr,nt,r,theta,tbar
     .     ,tsd)
c
      if (tbar .gt. 0.0) then
         itpw(14) = ifix(10.0*tbar + 0.5)
         itpw(15) = ifix(10.0*tsd  + 0.5)
      endif

c     ++itpw(16),tpw(17)
      r1  = 0.0
      r2  = 1000.0
      th1 = 0.0
      th2 = 360.0

      call caavgsd(tpwrt,tmiss,r1,r2,th1,th2,mr,mt,nr,nt,r,theta,tbar
     .     ,tsd)
c
      if (tbar .gt. 0.0) then
         itpw(16) = ifix(10.0*tbar + 0.5)
         itpw(17) = ifix(10.0*tsd  + 0.5)
      endif
c
c     Calculate upshear-relative padry in 90 deg arc
c     ++itpw(18)
      r1 = 000.0
      r2 = 500.0
      tmin = 45.0
      th1 = (180.0+sdir00-45.0)
      th2 = th1 + 90.0
      if (th1 .gt. 360.0) th1=th1-360.0
      if (th2 .gt. 360.0) th2=th2-360.0
      if (th1 .lt.   0.0) th1=th1+360.0
      if (th2 .lt.   0.0) th2=th2+360.0
c        write(6,*) k,th1,th2
      call padry(tpwrt,tmin,tmiss,r1,r2,th1,th2,
     +           mr,mt,nr,nt,r,theta,tbar)
      if (tbar .ge. 0.0) then
         itpw(18) = ifix(10.0*tbar + 0.5)
      endif
c
c     Calculate upshear-relative avg TPW in 90 deg arc
c     ++itpw(19)
      r1 = 000.0
      r2 = 500.0
      th1 = (180.0+sdir00-45.0)
      th2 = th1 + 90.0
      if (th1 .gt. 360.0) th1=th1-360.0
      if (th2 .gt. 360.0) th2=th2-360.0
      if (th1 .lt.   0.0) th1=th1+360.0
      if (th2 .lt.   0.0) th2=th2+360.0
c        write(6,*) k,th1,th2
      call caavg(tpwrt,tmiss,r1,r2,th1,th2,mr,mt,nr,nt,r,theta,tbar)
      if (tbar .gt. 0.0) then
         itpw(19) = ifix(10.0*tbar + 0.5)
      endif
c
c     Calculate large area symmetric avg TPW
c     ++itpw(20)
      r1 = 000.0
      r2 = 500.0
      th1 = 0.0
      th2 = 360.0
      call caavg(tpwrt,tmiss,r1,r2,th1,th2,mr,mt,nr,nt,r,theta,tbar)
      if (tbar .gt. 0.0) then
         itpw(20) = ifix(10.0*tbar + 0.5)
      endif
c
      return
      end
      subroutine lintcft(fxy,fmin,x1,dx,y1,dy,mx,my,nx,ny,xi,yi,fxyii)
c     This routine bi-linearly interpolates fxy to the point             
c     (xi,yi) to give fxyii.                                             
c                                                                        
c     Input:   fxy(mx,my) - array to be interpolated                     
c              fmin       - minimum value of f to include in the interpolation
c              dx,dy      - spacing of evenly spaced x,y points          
c              x1,y1      - x,y coorindates of lower-left point          
c              mx,my      - max dimensions of fxy                        
c              nx,ny      - working dimensions of fxy                    
c              xi,yi      - coordinates of point to be interpolated      
c                                                                        
c     Output:  fxyii      - value of interpolated function at (xi,yi)    
c                                                                        
      dimension fxy(mx,my)                                               
c                                                                        
c     Specify value to set fxyii if all surrounding points 
c     are missing
      fmiss = -99.
c     
c     Find the indices of the lower-left point of the               
c     grid box containing the point (xi,yi)                              
      i0 = 1 + ifix( (xi-x1)/dx )                                        
      j0 = 1 + ifix( (yi-y1)/dy )                                        
c                                                                        
c     Check index bounds                                                 
      if (i0 .lt.    1) i0=   1                                          
      if (i0 .gt. nx-1) i0=nx-1                                          
      if (j0 .lt.    1) j0=   1                                          
      if (j0 .gt. ny-1) j0=ny-1                                          
c                                                                        
      i1 = i0+1                                                          
      j1 = j0+1                                                          
c                                                                        
c     Calculate normalized x,y distances                                 
      xn = ( (xi-x1) - dx*float(i0-1) )/dx                               
      yn = ( (yi-y1) - dy*float(j0-1) )/dy                               
c                                                                        
      if (xn .lt. 0.0) xn = 0.0                                          
      if (xn .gt. 1.0) xn = 1.0                                          
      if (yn .lt. 0.0) yn = 0.0                                          
      if (yn .gt. 1.0) yn = 1.0                                          
c                                                                        
c     Calculate coefficients for interpolation function                  
      f00 = fxy(i0,j0)                                                   
      f10 = fxy(i1,j0)                                                   
      f01 = fxy(i0,j1)                                                   
      f11 = fxy(i1,j1)                                                   
c                                                                        
c     a = f00                                                            
c     b = f10-f00                                                        
c     c = f01-f00                                                        
c     d = f00+f11-f10-f01                                                
c                                                                        
c     fxyii = a + b*xn + c*yn + d*xn*yn                                  
c                                                                        
      w00 =  1 + xn*yn - (xn + yn)
      w10 = xn - xn*yn
      w01 = yn - xn*yn
      w11 =      xn*yn
c
c     Check for missing values and exclude them from the interpolation
      if (f00 .le. fmin) w00=0.0
      if (f10 .le. fmin) w10=0.0
      if (f01 .le. fmin) w01=0.0
      if (f11 .le. fmin) w11=0.0
c   
      fxyii = w00*f00 + w10*f10 + w01*f01 + w11*f11
c   
      wsum = w00 + w10 + w01 + w11
      if (wsum .le. 0.0) then
         fxyii = fmiss
      else
         fxyii = fxyii/wsum
      endif
c
c     xyn = xn*yn
c     a0  = 1.0
c     write(6,222) a,b,c,d,a0,xn,yn,xyn,fxyii,
c    +             w00,w10,w01,w11,f00,f10,f01,f11,fxyjj
c 222 format(/,'  a, b, c, d,  1, x, y,xy,fii: ',9(f8.4,1x),/,
c    +         'w00,10,01,11,f00,10,01,11,fjj: ',9(f8.4,1x))
c
      return                                                             
      end
      subroutine caavg(frt,fmin,r1,r2,th1,th2,mr,mt,nr,nt,r,theta,fbar)
c     This routine calculates the area average of f(r,theta)
c     from r=r1 to r=r2. Values less than fmin are not included
c     in the average. 
c
      dimension frt(0:mr,0:mt)
      dimension r(0:mr),theta(0:mt)
c
c     Local variables
      parameter (mtt=1000)
      dimension tmask(0:mtt)
c
      if (mtt .lt. nt) then
         write(6,*) 'mtt too small in routine caavg'
         stop
      endif
c
c     Determine which theta values to include in the azimuthal average
      dth1 = 10000.0
      dth2 = 10000.0
      nt1  = -99
      nt2  = -99
      do n=0,nt
         dist1 = abs(th1-theta(n))
         if (dist1 .lt. dth1) then
            dth1 = dist1
            nt1  = n
         endif
c
         dist2 = abs(th2-theta(n))
         if (dist2 .lt. dth2) then
            dth2 = dist2
            nt2 = n
         endif
      enddo
c
      if (nt2 .gt. nt1) then
         do n=0,nt
            if (n .ge. nt1 .and. n .le. nt2) then
               tmask(n) = 1.0
            else
               tmask(n) = 0.0
            endif
         enddo
      elseif (nt1 .gt. nt2) then
         do n=0,nt
            if (n .gt. nt2 .and. n .lt. nt1) then
               tmask(n) = 0.0
            else
               tmask(n) = 1.0
            endif
         enddo
      else
         do n=0,nt
            if (n .eq. nt1 .and. n .eq. nt2) then
               tmask(n) = 1.0
            else
               tmask(n) = 0.0
            endif
         enddo
      endif
c
c     write(6,*) 'caavg: nt1,nt2',nt1,nt2
c     write(6,777) (theta(kk),kk=0,nt)
c     write(6,777) (tmask(kk),kk=0,nt)
c 777 format(20(f6.1))
c
c
c     Find indices corresponding to r1 and r2
      ir1 = 0
      ir2 = nr
c
      do i=nr,0,-1
         if (r(i) .le. r1) then
            ir1 = i
            exit
         endif
      enddo
c
      do i=1,nr
         if (r(i) .ge. r2) then
            ir2 = i
            exit
         endif
      enddo
c
      if (ir1 .ge. ir2) then
         fbar = 0.0
         return
      endif
c
      fbar = 0.0
      abar = 0.0
      do i=ir1,ir2
         do j=0,nt-1
            if (frt(i,j) .ge. fmin) then
               abar = abar + r(i)*tmask(j)
               fbar = fbar + r(i)*frt(i,j)*tmask(j)
            endif
         enddo
      enddo
c           
      if (abar .gt. 0.0) then
         fbar = fbar/abar
      else
         fbar = 0.0
      endif
c
      return
      end
      subroutine padry(frt,fmin,fmiss,r1,r2,th1,th2,
     +                  mr,mt,nr,nt,r,theta,fper)
c     This routine calculates the percent of area of f(r,theta)
c     from r=r1 to r=r2  that is less than fmin (fper). Values less than 
c     fmiss are not included in the calculation.
c
      dimension frt(0:mr,0:mt)
      dimension r(0:mr),theta(0:mt)
c
c     Local variables
      parameter (mtt=1000)
      dimension tmask(0:mtt)
c
      if (mtt .lt. nt) then
         write(6,*) 'mtt too small in routine padry'
         stop
      endif
c
c     Determine which theta values to include in the azimuthal average
      dth1 = 10000.0
      dth2 = 10000.0
      nt1  = -99
      nt2  = -99
      do n=0,nt
         dist1 = abs(th1-theta(n))
         if (dist1 .lt. dth1) then
            dth1 = dist1
            nt1  = n
         endif
c
         dist2 = abs(th2-theta(n))
         if (dist2 .lt. dth2) then
            dth2 = dist2
            nt2 = n
         endif
      enddo
c
      if (nt2 .gt. nt1) then
         do n=0,nt
            if (n .ge. nt1 .and. n .le. nt2) then
               tmask(n) = 1.0
            else
               tmask(n) = 0.0
            endif
         enddo
      elseif (nt1 .gt. nt2) then
         do n=0,nt
            if (n .gt. nt2 .and. n .lt. nt1) then
               tmask(n) = 0.0
            else
               tmask(n) = 1.0
            endif
         enddo
      else
         do n=0,nt
            if (n .eq. nt1 .and. n .eq. nt2) then
               tmask(n) = 1.0
            else
               tmask(n) = 0.0
            endif
         enddo
      endif
c
c     write(6,*) 'padry: nt1,nt2',nt1,nt2
c     write(6,777) (theta(kk),kk=0,nt)
c     write(6,777) (tmask(kk),kk=0,nt)
c 777 format(20(f6.1))
c
c
c     Find indices corresponding to r1 and r2
      ir1 = 0
      ir2 = nr
c
      do i=nr,0,-1
         if (r(i) .le. r1) then
            ir1 = i
            exit
         endif
      enddo
c
      do i=1,nr
         if (r(i) .ge. r2) then
            ir2 = i
            exit
         endif
      enddo
c
      if (ir1 .ge. ir2) then
         fper = 0.0
         return
      endif
c
      fper  = 0.0
      abar  = 0.0
      rcount = 0.0
      do i=ir1,ir2
         do j=0,nt-1
            if (frt(i,j) .ge. fmiss) then
               abar = abar + r(i)*tmask(j)
               rcount = rcount + 1.0*tmask(j)
               if (frt(i,j) .le. fmin) then
                  fper = fper + r(i)*tmask(j)
               endif
            endif
         enddo
      enddo
c           
      if (abar .gt. 0.0) then
         fper = 100.0*fper/abar
      else
         fper = -1.0
      endif
c
c     write(6,*) 'ir1,ir2,rcount=',ir1,ir2,rcount
      return
      end
      subroutine caavgsd(frt,fmin,r1,r2,th1,th2,mr,mt,nr,nt,r,theta,fbar
     .     ,fsigma)
c     This routine calculates the area average of f(r,theta)
c     from r=r1 to r=r2. Values less than fmin are not included
c     in the average. 
c
      dimension frt(0:mr,0:mt)
      dimension r(0:mr),theta(0:mt)
c
c     Local variables
      parameter (mtt=1000)
      dimension tmask(0:mtt)
c
      if (mtt .lt. nt) then
         write(6,*) 'mtt too small in routine caavg'
         stop
      endif
c
c     Determine which theta values to include in the azimuthal average
      dth1 = 10000.0
      dth2 = 10000.0
      nt1  = -99
      nt2  = -99
      do n=0,nt
         dist1 = abs(th1-theta(n))
         if (dist1 .lt. dth1) then
            dth1 = dist1
            nt1  = n
         endif
c
         dist2 = abs(th2-theta(n))
         if (dist2 .lt. dth2) then
            dth2 = dist2
            nt2 = n
         endif
      enddo
c
      if (nt2 .gt. nt1) then
         do n=0,nt
            if (n .ge. nt1 .and. n .le. nt2) then
               tmask(n) = 1.0
            else
               tmask(n) = 0.0
            endif
         enddo
      elseif (nt1 .gt. nt2) then
         do n=0,nt
            if (n .gt. nt2 .and. n .lt. nt1) then
               tmask(n) = 0.0
            else
               tmask(n) = 1.0
            endif
         enddo
      else
         do n=0,nt
            if (n .eq. nt1 .and. n .eq. nt2) then
               tmask(n) = 1.0
            else
               tmask(n) = 0.0
            endif
         enddo
      endif
c
c     write(6,*) 'caavg: nt1,nt2',nt1,nt2
c     write(6,777) (theta(kk),kk=0,nt)
c     write(6,777) (tmask(kk),kk=0,nt)
c 777 format(20(f6.1))
c
c
c     Find indices corresponding to r1 and r2
      ir1 = 0
      ir2 = nr
c
      do i=nr,0,-1
         if (r(i) .le. r1) then
            ir1 = i
            exit
         endif
      enddo
c
      do i=1,nr
         if (r(i) .ge. r2) then
            ir2 = i
            exit
         endif
      enddo
c
      if (ir1 .ge. ir2) then
         fbar = 0.0
         return
      endif
c
      fbar = 0.0
      abar = 0.0
      fsigma = 0.0
      do i=ir1,ir2
         do j=0,nt-1
            if (frt(i,j) .ge. fmin) then
               abar = abar + r(i)*tmask(j)
               fbar = fbar + r(i)*frt(i,j)*tmask(j)
               fsigma = fsigma + r(i)*frt(i,j)*frt(i,j)*tmask(j)
            endif
         enddo
      enddo
c           
      if (abar .gt. 0.0) then
         fbar = fbar/abar
      else
         fbar = 0.0
      endif
c     print*,'variance, mean**2',fsigma/abar,fbar**2
      if (fsigma/abar.gt.fbar**2)then
         fsigma=sqrt((fsigma)/abar -fbar**2)
      else
         fsigma=1.0
      end if

c
      return
      end