w3ft204.f

Go to the documentation of this file.

C> @file
C> @brief Convert (361,181) grid to (93,68) mercator grid.
C> @author Ralph Jones @date 1994-05-18
 
C> Convert a n. s. hemisphere 1.0 degree lat.,lon. 361 by
C> 181 grid to a national - hawaii (mercator) 93*68 awips 204
C> grid.
C>
C> ### Program History Log:
C> Date | Programmer | Comment
C> -----|------------|--------
C> 1994-05-18 | Ralph Jones | Initial.
C>
C> @param[in] ALOLA 361*181 grid 1.0 deg. lat,lon grid n. hemi.
C> 65341 point grid. 360 * 181 one degree grib grid 3 was flipped, greenwish added
C> to right side.
C> @param[in] INTERP 1 linear interpolation , ne.1 biquadratic
C> @param[out] AMERC 93*68 grid national - hawaii (mercator) 6324 point grid
C> is awips grid type 204
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 1994-05-18
      SUBROUTINE w3ft204(ALOLA,AMERC,INTERP)
C
       parameter(npts=6324,ii=93,jj=68)
       parameter(alatin=20.000)
       parameter(pi=3.1416)
       parameter(dx=160000.0)
       parameter(alat1=-25.000)
       parameter(alon1=110.000)
C
       REAL        WLON(NPTS),    XLAT(NPTS)
       REAL        XI(II,JJ),     XJ(II,JJ)
       REAL        XII(NPTS),     XJJ(NPTS)
       REAL        ALOLA(361,181), 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
      radpd  = pi / 180.0
      degpr  = 180.0 / pi
      clain  = cos(radpd * alatin)
      dellon = dx / (rerth * clain)
      djeo   = (alog(tan(0.5*((alat1+90.0)*radpd))))/dellon
C
      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) + 91.0
 400  CONTINUE
C
      iswt   = 1
      intrpo = interp
      GO TO 1000
C
C     AFTER THE 1ST CALL TO W3FT204 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
        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