w3ft202.f

Go to the documentation of this file.

C> @file
C> @brief Convert (361,91) grid to (65,43) n. hemi. grid
C> @author Ralph Jones @date 1994-05-18
 
C> Convert a northern hemisphere 1.0 degree lat.,lon. 361 by
C> 91 grid to a polar stereographic 65 by 43 grid. The polar
C> stereographic map projection is true at 60 deg. n. , The mesh
C> length is 190.5 km. and the oriention is 105 deg. w.
C>
C> ### Program History Log:
C> Date | Programmer | Comment
C> -----|------------|--------
C> 1994-05-18 | Ralph Jones | Initial.
C>
C> @param[in] ALOLA 361*91 grid 1.0 lat,lon grid n. hemisphere 32851 point
C> grid is o.n. 84 type ?? or ?? hex
C> @param[in] INTERP 1 linear interpolation , ne.1 biquadratic
C> @param[out] APOLA 65*43 grid of northern hemisphere. 2795 point grid is
C> awips grid type 202
C>
C> @note
C> - 1. W1 and w2 are used to store sets of constants which are
C> reusable for repeated calls to the subroutine.
C> - 2. Wind components are not rotated to the 65*43 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*43
C> grid (about 1100 points).
C>
C> @author Ralph Jones @date 1994-05-18
      SUBROUTINE w3ft202(ALOLA,APOLA,INTERP)
C
       parameter(npts=2795,ii=65,jj=43)
       parameter(orient=105.0,ipole=33,jpole=45)
       parameter(xmesh=190.5)
C
       REAL        R2(NPTS),      WLON(NPTS)
       REAL        XLAT(NPTS),    XI(II,JJ),   XJ(II,JJ)
       REAL        XII(NPTS),     XJJ(NPTS),   ANGLE(NPTS)
       REAL        ALOLA(361,91), APOLA(NPTS), ERAS(NPTS,4)
       REAL        W1(NPTS),      W2(NPTS)
       REAL        XDELI(NPTS),   XDELJ(NPTS)
       REAL        XI2TM(NPTS),   XJ2TM(NPTS)
C
       INTEGER     IV(NPTS),      JV(NPTS),    JY(NPTS,4)
       INTEGER     IM1(NPTS),     IP1(NPTS),   IP2(NPTS)
C
       LOGICAL     LIN
C
       SAVE
C
       equivalence(xi(1,1),xii(1)),(xj(1,1),xjj(1))
C
       DATA  degprd/57.2957795/
       DATA  earthr/6371.2/
       DATA  intrpo/99/
       DATA  iswt  /0/
C
      lin = .false.
      IF (interp.EQ.1) lin = .true.
C
      IF  (iswt.EQ.1)  GO TO  900
C
        deg    = 1.0
        gi2    = (1.86603 * earthr) / xmesh
        gi2    = gi2 * gi2
C
C     NEXT 32 LINES OF CODE PUTS SUBROUTINE W3FB01 IN LINE
C
      DO 100 j = 1,jj
         xj1 = j - jpole
         DO 100 i = 1,ii
             xi(i,j) = i - ipole
             xj(i,j) = xj1
 100     CONTINUE
C
      DO 200 kk = 1,npts
        r2(kk)   = xjj(kk) * xjj(kk) + xii(kk) * xii(kk)
        xlat(kk) = degprd *
     &      asin((gi2 - r2(kk)) / (gi2 + r2(kk)))
 200  CONTINUE
C
      DO 300 kk = 1,npts
        angle(kk) = degprd * atan2(xjj(kk),xii(kk))
 300  CONTINUE
C
      DO 400 kk = 1,npts
        IF (angle(kk).LT.0.0) angle(kk) = angle(kk) + 360.0
 400  CONTINUE
C
      DO 500 kk = 1,npts
        wlon(kk) = 270.0 + orient - angle(kk)
 500  CONTINUE
C
      DO 600 kk = 1,npts
        IF (wlon(kk).LT.0.0)   wlon(kk) = wlon(kk) + 360.0
 600  CONTINUE
C
      DO 700 kk = 1,npts
        IF (wlon(kk).GE.360.0) wlon(kk) = wlon(kk) - 360.0
 700  CONTINUE
C
      DO 800 kk = 1,npts
        w1(kk)  = (360.0 - wlon(kk)) / deg + 1.0
        w2(kk)  = xlat(kk) / deg + 1.0
 800  CONTINUE
C
      iswt   = 1
      intrpo = interp
      GO TO 1000
C
C     AFTER THE 1ST CALL TO W3FT202 TEST INTERP, IF IT HAS
C     CHANGED RECOMPUTE SOME CONSTANTS
C
  900 CONTINUE
        IF (interp.EQ.intrpo) GO TO 2100
        intrpo = interp
C
 1000 CONTINUE
        DO 1100 k = 1,npts
          iv(k)    = w1(k)
          jv(k)    = w2(k)
          xdeli(k) = w1(k) - iv(k)
          xdelj(k) = w2(k) - jv(k)
          ip1(k)   = iv(k) + 1
          jy(k,3)  = jv(k) + 1
          jy(k,2)  = jv(k)
 1100   CONTINUE
C
      IF (lin) GO TO 2100
C
      DO 1200 k = 1,npts
        ip2(k)   = iv(k) + 2
        im1(k)   = iv(k) - 1
        jy(k,1)  = jv(k) - 1
        jy(k,4)  = jv(k) + 2
        xi2tm(k) = xdeli(k) * (xdeli(k) - 1.0) * .25
        xj2tm(k) = xdelj(k) * (xdelj(k) - 1.0) * .25
 1200 CONTINUE
C
 2100 CONTINUE
      IF (lin) THEN
C
C     LINEAR INTERPOLATION
C
      DO 2200 kk = 1,npts
        eras(kk,2) = (alola(ip1(kk),jy(kk,2))-alola(iv(kk),jy(kk,2)))
     &             * xdeli(kk) + alola(iv(kk),jy(kk,2))
        eras(kk,3) = (alola(ip1(kk),jy(kk,3))-alola(iv(kk),jy(kk,3)))
     &             * xdeli(kk) + alola(iv(kk),jy(kk,3))
 2200 CONTINUE
C
      DO 2300 kk = 1,npts
        apola(kk) = eras(kk,2) + (eras(kk,3) - eras(kk,2))
     &            * xdelj(kk)
 2300 CONTINUE
C
      ELSE
C
C     QUADRATIC INTERPOLATION
C
      DO 2400 kk = 1,npts
        eras(kk,1)=(alola(ip1(kk),jy(kk,1))-alola(iv(kk),jy(kk,1)))
     &            * xdeli(kk) + alola(iv(kk),jy(kk,1)) +
     &            ( alola(im1(kk),jy(kk,1)) - alola(iv(kk),jy(kk,1))
     &            - alola(ip1(kk),jy(kk,1))+alola(ip2(kk),jy(kk,1)))
     &            * xi2tm(kk)
        eras(kk,2)=(alola(ip1(kk),jy(kk,2))-alola(iv(kk),jy(kk,2)))
     &            * xdeli(kk) + alola(iv(kk),jy(kk,2)) +
     &            ( alola(im1(kk),jy(kk,2)) - alola(iv(kk),jy(kk,2))
     &            - alola(ip1(kk),jy(kk,2))+alola(ip2(kk),jy(kk,2)))
     &            * xi2tm(kk)
        eras(kk,3)=(alola(ip1(kk),jy(kk,3))-alola(iv(kk),jy(kk,3)))
     &            * xdeli(kk) + alola(iv(kk),jy(kk,3)) +
     &            ( alola(im1(kk),jy(kk,3)) - alola(iv(kk),jy(kk,3))
     &            - alola(ip1(kk),jy(kk,3))+alola(ip2(kk),jy(kk,3)))
     &            * xi2tm(kk)
        eras(kk,4)=(alola(ip1(kk),jy(kk,4))-alola(iv(kk),jy(kk,4)))
     &            * xdeli(kk) + alola(iv(kk),jy(kk,4)) +
     &            ( alola(im1(kk),jy(kk,4)) - alola(iv(kk),jy(kk,4))
     &            - alola(ip1(kk),jy(kk,4))+alola(ip2(kk),jy(kk,4)))
     &            * xi2tm(kk)
 2400      CONTINUE
C
       DO 2500 kk = 1,npts
         apola(kk) = eras(kk,2) + (eras(kk,3) - eras(kk,2))
     &             * xdelj(kk)  + (eras(kk,1) - eras(kk,2)
     &             - eras(kk,3) + eras(kk,4)) * xj2tm(kk)
 2500  CONTINUE
C
      ENDIF
C
      RETURN
      END