w3ft06.f

Go to the documentation of this file.

C> @file
C> @brief Convert (145,37) to (65,65) s. hemi. grid.
C> @author Ralph Jones @date 1984-06-18
 
C> Convert a southern 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. s.; The mesh
C> length is 381 km. and the oriention is 260 deg. w.(100e).
C>
C> ### Program History Log:
C> Date | Programmer | Comment
C> -----|------------|--------
C> 1984-06-18 | 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 deg 2.5 lat,lon grid s. hemi.
C> 5365 point grid is type 30 or 1e hex o.n. 84.
C> @param[in] LINEAR 1 linear interpolation , ne.1 biquadratic.
C> @param[out] APOLA 65*65 grid of southern hemi.
C> 4225 point grid is type 28 or 1c hex o.n. 84.
C> @param[out] W1 65*65 scratch field.
C> @param[out] W2 65*65 scratch field. FT06F001 Error message
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 w3fc10() 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 verion w3ft06v() on the cray
C> it has 3 parameters in the call, runs about times faster, uses
C> more memory.
C>
C> @author Ralph Jones @date 1984-06-18
      SUBROUTINE w3ft06(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 w3ft06 ***)
C
         lin = .false.
         IF (linear.EQ.1) lin = .true.
         IF  (iswt.EQ.0)  GO TO  300
C
C        TEST TO SEE IF W1 OR W2 WAS 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
         GO TO 800
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 600 PUTS SUBROUTINE W3FB03 IN LINE
C
         DO  600 j=1,65
           xj = j - 33
           xj2 = xj * xj
           DO  600 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  500
 400         CONTINUE
             xlong = degprd * atan2(xj,xi)
             wlon = xlong -10.0
             IF (wlon.LT.0.0)  wlon = wlon + 360.0
             xlat = asin((gi2-r2)/(gi2+r2))*degprd
             xlat = -xlat
 500         CONTINUE
             xlat = xlat + 90.0
             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
  600      CONTINUE
C
         DO 700  kk=1,10
           savew1(kk)=w1(kk)
           savew2(kk)=w2(kk)
 700     CONTINUE
C
           iswt = 1
C
 800     CONTINUE
C
         DO  1900  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  900
           ip2 = i + 2
           im1 = i - 1
           jy(4) = j + 2
           jy(1) = j - 1
           xi2tm = xdeli*(xdeli-1.)*.25
           xj2tm = xdelj*(xdelj-1.)*.25
 900     CONTINUE
         IF ((i.LT.2).OR.(j.LT.2)) GO TO 1000
         IF ((i.GT.142).OR.(j.GT.34)) GO TO  1000
C     QUADRATIC (LINEAR TOO) OK W/O FURTHER ADO SO GO TO 1500
         GO TO  1500
C
 1000    CONTINUE
         IF (i.EQ.1) GO TO 1100
         IF (i.EQ.144) GO TO 1200
         ip2 = i+2
         im1 = i-1
         GO TO  1300
C
 1100    CONTINUE
         ip2 = 3
         im1 = 144
         GO TO  1300
C
 1200    CONTINUE
         ip2 = 2
         im1 = 143
C
 1300    CONTINUE
         ip1 = i + 1
         IF (lin)  GO TO  1400
         IF ((j.LT.2).OR.(j.GE.36))  xj2tm=0.
C.....DO NOT ALLOW POINT OFF GRID,CYCLIC OR NOT
         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
 1400    CONTINUE
C.....DO NOT ALLOW POINT OFF GRID,CYCLIC OR NOT
         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
 1500    CONTINUE
C
         IF (.NOT.lin)  GO TO  1700
C
C        LINEAR INTERPOLATION
C
          DO 1600 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)
 1600    CONTINUE
C
         apola(kk) = eras(2) + (eras(3) - eras(2)) * xdelj
         GO TO  1900
C
 1700    CONTINUE
C
C        QUADRATIC INTERPOLATION
C
         DO 1800  k = 1,4
           j1 = jy(k)
C.....DO NOT ALLOW POINT OFF GRID,CYCLIC OR NOT
           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
 1800    CONTINUE
C
         apola(kk) = eras(2)+(eras(3)-eras(2))*xdelj+(eras(1)-
     &               eras(2)-eras(3)+eras(4))*xj2tm
C
 1900    CONTINUE
C
C        SET POLE POINT, WMO STANDARD FOR U OR V
C
         apola(2113) = alola(73,1)
C
         RETURN
       END