w3ft05.f

Go to the documentation of this file.

C> @file
C> @brief Convert (145,37) to (65,65) n. hemi. grid
C> @author Ralph Jones @date 1985-04-08
 
C> Convert a northern hemisphere 2.5 degree lat.,lon. 145 by
C> 37 grid to a polar stereographic 65 by 65 grid. The polar
C> stereographic map projection is true at 60 deg. n. , The mesh
C> length is 381 km. and the oriention is 80 deg. w.
C>
C> ### Program History Log:
C> Date | Programmer | Comment
C> -----|------------|--------
C> 1985-04-08 | Ralph Jones | Initial.
C> 1991-07-30 | Ralph Jones | convert to cray cft77 fortran.
C> 1992-05-02 | Ralph Jones | add save.
C>
C> @param[in] ALOLA 145*37 grid 2.5 lat,lon grid n. hemi.
C> 5365 point grid is type 29 or 1d hex o.n. 84
C> @param[in] LINEAR 1 linear interpolation , ne.1 biquadratic
C> @param[out] APOLA 65*65 grid of northern hemi.
C> 4225 point grid is type 27 or 1b hex o.n. 84
C> @param[out] W1 65*65 scratch field
C> @param[out] W2 65*65 scratch field
C>
C> @remark
C> - 1. W1 and w2 are used to store sets of constants which are
C> reusable for repeated calls to the subroutine. If they are
C> over written by the user, a warning message will be printed
C> and w1 and w2 will be recomputed.
C> - 2. Wind components are not rotated to the 65*65 grid orientation
C> after interpolation. You may use w3fc08() to do this.
C> - 3. The grid points values on the equator have been extrapolated
C> outward to all the grid points outside the equator on the 65*65
C> grid (about 1100 points).
C> - 4. You should use the cray vectorized version w3ft05v on the cray
C> it has 3 parameters in the call, runs about 10 times faster. Uses
C> more memory.
C>
C> @author Ralph Jones @date 1985-04-08
      SUBROUTINE w3ft05(ALOLA,APOLA,W1,W2,LINEAR)
C
       REAL            ALOLA(145,37)
       REAL            APOLA(4225)
       REAL            ERAS(4)
       REAL            SAVEW1(10)
       REAL            SAVEW2(10)
       REAL            W1(4225)
       REAL            W2(4225)
C
       INTEGER         JY(4)
       INTEGER         OUT
C
       LOGICAL         LIN
C
       SAVE
C
       DATA  degprd/57.2957795/
       DATA  earthr/6371.2/
       DATA  iswt  /0/
       DATA  out   /6/
C
 4000  FORMAT ( 52h *** warning , w1 or w2 scratch files over written ,,
     &          43h i will restore them , burning up cpu time,,
     &          14h in w3ft05 ***)
C
         lin = .false.
         IF (linear.EQ.1) lin = .true.
C
         IF  (iswt.EQ.0)  GO TO  300
C
C        TEST W1 AND W2 TO SEE IF THEY WERE WRITTEN OVER
C
         DO 100  kk=1,10
           IF (savew1(kk).NE.w1(kk)) GO TO  200
           IF (savew2(kk).NE.w2(kk)) GO TO  200
  100    CONTINUE
         GOTO  1000
C
  200    CONTINUE
         WRITE (out,4000)
C
  300    CONTINUE
         deg   = 2.5
         nn    = 0
         xmesh = 381.0
         gi2   = (1.86603*earthr) / xmesh
         gi2   = gi2 * gi2
C
C        DO LOOP 800 PUTS SUBROUTINE W3FB01 IN LINE
C
         DO  800  j = 1,65
           xj  = j - 33
           xj2 = xj * xj
           DO  800  i=1,65
             xi = i - 33
             r2 = xi*xi + xj2
             IF  (r2.NE.0.0)  GO TO  400
             wlon = 0.0
             xlat = 90.0
             GO TO  700
 400         CONTINUE
             xlong = degprd * atan2(xj,xi)
             IF (xlong.GE.0.0)  GO TO 500
             wlon  = -10.0 - xlong
             IF (wlon.LT.0.0)  wlon = wlon + 360.0
             GO TO  600
C
 500         CONTINUE
             wlon = 350.0 - xlong
 600         CONTINUE
             xlat = asin((gi2-r2)/(gi2+r2))*degprd
 700         CONTINUE
             IF  (wlon.GT.360.0)  wlon = wlon - 360.0
             IF  (wlon.LT.0.0)   wlon = wlon + 360.0
             nn     = nn + 1
             w1(nn) = ( 360.0 - wlon ) / deg + 1.0
             w2(nn) = xlat / deg + 1.0
 800       CONTINUE
C
         DO 900  kk = 1,10
           savew1(kk) = w1(kk)
           savew2(kk) = w2(kk)
 900     CONTINUE
C
         iswt = 1
C
 1000    CONTINUE
C
         DO  2100  kk = 1,4225
           i     = w1(kk)
           j     = w2(kk)
           fi    = i
           fj    = j
           xdeli = w1(kk) - fi
           xdelj = w2(kk) - fj
           ip1   = i + 1
           jy(3) = j + 1
           jy(2) = j
           IF (lin)  GO TO 1100
           ip2   = i + 2
           im1   = i - 1
           jy(4) = j + 2
           jy(1) = j - 1
           xi2tm = xdeli * (xdeli-1.) * 0.25
           xj2tm = xdelj * (xdelj-1.) * 0.25
C
 1100    CONTINUE
           IF ((i.LT.2).OR.(j.LT.2))    GO TO  1200
           IF ((i.GT.142).OR.(j.GT.34)) GO TO  1200
C
C     QUADRATIC (LINEAR TOO) OK W/O FURTHER ADO SO GO TO 1700
C
           GO TO  1700
C
 1200    CONTINUE
           IF (i.EQ.1)   GO TO  1300
           IF (i.EQ.144) GO TO 1400
           ip2 = i + 2
           im1 = i - 1
           GO TO  1500
C
 1300    CONTINUE
           ip2 = 3
           im1 = 144
           GO TO   1500
C
 1400    CONTINUE
           ip2 = 2
           im1 = 143
C
 1500   CONTINUE
          ip1 = i + 1
          IF (lin)  GO TO 1600
          IF ((j.LT.2).OR.(j.GE.36))  xj2tm=0.
C.....DO NOT ALLOW POINT OFF GRID
          IF (ip2.LT.1)   ip2 = 1
          IF (im1.LT.1)   im1 = 1
          IF (ip2.GT.145) ip2 = 145
          IF (im1.GT.145) im1 = 145
C
 1600   CONTINUE
C.....DO NOT ALLOW POINT OFF GRID
          IF (i.LT.1)     i   = 1
          IF (ip1.LT.1)   ip1 = 1
          IF (i.GT.145)   i   = 145
          IF (ip1.GT.145) ip1 = 145
C
 1700    CONTINUE
         IF (.NOT.lin)  GO TO  1900
C
C        LINEAR INTERPLOATION
C
         DO 1800 k = 2,3
           j1 = jy(k)
           IF  (j1.LT.1)  j1 = 1
           IF  (j1.GT.37) j1 = 37
           eras(k) = (alola(ip1,j1) - alola(i,j1)) * xdeli + alola(i,j1)
 1800    CONTINUE
C
           apola(kk) = eras(2) + (eras(3) - eras(2)) * xdelj
           GO TO  2100
C
 1900    CONTINUE
C
C        QUADRATIC INTERPOLATION
C
         DO 2000 k = 1,4
           j1 = jy(k)
C.....DO NOT ALLOW POINT OFF GRID
           IF (j1.LT.1) j1 = 1
           IF (j1.GT.37) j1 = 37
           eras(k) = (alola(ip1,j1)-alola(i,j1))*xdeli+alola(i,j1)+
     &               (alola(im1,j1)-alola(i,j1)-alola(ip1,j1)+
     &                alola(ip2,j1))*xi2tm
 2000    CONTINUE
C
         apola(kk) = eras(2)+(eras(3)-eras(2))*xdelj+(eras(1)-
     &               eras(2)-eras(3)+eras(4)) * xj2tm
C
 2100    CONTINUE
C
C        SET POLE POINT , WMO STANDARD FOR U OR V
C
         apola(2113) = alola(73,37)
C
         RETURN
       END