w3ft210.f

Go to the documentation of this file.

C> @file
C> @brief Convert (361,91) grid to (25,25) mercator grid.
C> @author Ralph Jones @date 1993-10-19
 
C> Convert a northern hemisphere 1.0 degree lat.,lon. 361 by
C> 91 grid to a regional - puerto rico (mercator) 25*25 awips 210
C> grid.
C>
C> ### Program History Log:
C> Date | Programmer | Comment
C> -----|------------|--------
C> 1993-10-19 | Ralph Jones | Initial.
 
C> @param[in] ALOLA 361*91 grid 1.0 deg. lat,lon grid n. hemi.
C> 32851 point grid. 360 * 181 one degree grib grid 3 was flipped, greenwish added
C> to right side and cut to 361 * 91.
C> @param[in] INTERP 1 linear interpolation , ne.1 biquadratic
C> @param[out] AMERC 25*25 grid of northern mercator 625 point grid is awips grid type 210
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. 20 other array
C> are saved and reused on the next call.
C>
C> @author Ralph Jones @date 1993-10-19
      SUBROUTINE w3ft210(ALOLA,AMERC,INTERP)
C
       parameter(npts=625,ii=25,jj=25)
       parameter(alatin=20.000)
       parameter(pi=3.1416)
       parameter(dx=80000.0)
       parameter(alat1=9.000)
       parameter(alon1=283.000)
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), AMERC(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
C      DATA  DEGPR /57.2957795/
       DATA  rerth /6.3712e+6/
       DATA  intrpo/99/
       DATA  iswt  /0/
C
      degpr  = 180.0 / pi
      radpd  = pi / 180.0
      clain  = cos(radpd * alatin)
      dellon = dx / (rerth * clain)
      djeo   = (alog(tan(0.5*((alat1+90.0)*radpd))))/dellon
      lin = .false.
      IF (interp.EQ.1) lin = .true.
C
      IF  (iswt.EQ.1)  GO TO  900
C
        deg    = 1.0
C
C     NEXT 32 LINES OF CODE PUTS SUBROUTINE W3FB09 IN LINE
C
      DO 100 j = 1,jj
         DO 100 i = 1,ii
             xi(i,j) = i
             xj(i,j) = j
 100     CONTINUE
C
      DO 200 kk = 1,npts
        xlat(kk) = 2.0*atan(exp(dellon*(djeo + xjj(kk)-1.)))
     &              * degpr - 90.0
 200  CONTINUE
C
      DO 300 kk = 1,npts
        wlon(kk) = (xii(kk) -1.0) * dellon * degpr + alon1
 300  CONTINUE
C
      DO 400 kk = 1,npts
        w1(kk)  = wlon(kk) + 1.0
        w2(kk)  = xlat(kk) + 1.0
 400  CONTINUE
C
      iswt   = 1
      intrpo = interp
      GO TO 1000
C
C     AFTER THE 1ST CALL TO W3FT210 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 (.NOT.lin) THEN
        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
      END IF
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
        amerc(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
         amerc(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