w3fb11.f

Go to the documentation of this file.

C> @file
C> @brief Lat/lon to lambert(i,j) for grib.
C> @author John Stackpole @date 1988-11-25
 
C> Converts the coordinates of a location on Earth given in
C> the natural coordinate system of latitude/longitude to a grid
C> coordinate system overlaid on a lambert conformal tangent cone
C> projection true at a given n or s latitude. w3fb11() is the reverse
C> of w3fb12(). uses grib specification of the location of the grid
C>
C> Program history log:
C> - John Stackpole 1988-11-25
C> - Ralph Jones 1990-04-12 Convert to cft77 fortran.
C> - Ralph Jones 1994-04-28 Add save statement.
C>
C> @param[in] ALAT Latitude in degrees (negative in southern hemis).
C> @param[in] ELON East longitude in degrees, real*4.
C> @param[in] ALAT1 Latitude of lower left point of grid (point (1,1)).
C> @param[in] ELON1 Longitude of lower left point of grid (point (1,1))
C> all real*4.
C> @param[in] DX Mesh length of grid in meters at tangent latitude.
C> @param[in] ELONV The orientation of the grid. i.e.,
C> the east longitude value of the vertical meridian
C> which is parallel to the y-axis (or columns of
C> of the grid) along which latitude increases as
C> the y-coordinate increases. real*4
C> this is also the meridian (on the back side of the
C> tangent cone) along which the cut is made to lay
C> the cone flat.
C> @param[in] ALATAN The latitude at which the lambert cone is tangent to
C> (touching) the spherical Earth. Set negative to indicate a
C> southern hemisphere projection.
C> @param[out] XI I coordinate of the point specified by alat, elon
C> @param[out] XJ J coordinate of the point; both real*4
C>
C> @note Formulae and notation loosely based on hoke, hayes,
C> and renninger's "map projections and grid systems...", march 1981
C> afgwc/tn-79/003.
C>
C> @author John Stackpole @date 1988-11-25
      SUBROUTINE w3fb11(ALAT,ELON,ALAT1,ELON1,DX,ELONV,ALATAN,XI,XJ)
C
         SAVE
C
         DATA  rerth /6.3712e+6/, pi/3.14159/
C
C        PRELIMINARY VARIABLES AND REDIFINITIONS
C
C        H = 1 FOR NORTHERN HEMISPHERE; = -1 FOR SOUTHERN
C
         IF (alatan.GT.0) THEN
           h = 1.
         ELSE
           h = -1.
         ENDIF
C
         radpd  = pi    / 180.0
         rebydx = rerth / dx
         alatn1 = alatan * radpd
         an     = h * sin(alatn1)
         cosltn = cos(alatn1)
C
C        MAKE SURE THAT INPUT LONGITUDES DO NOT PASS THROUGH
C        THE CUT ZONE (FORBIDDEN TERRITORY) OF THE FLAT MAP
C        AS MEASURED FROM THE VERTICAL (REFERENCE) LONGITUDE.
C
         elon1l = elon1
         IF ((elon1 - elonv).GT.180.)
     &     elon1l = elon1 - 360.
         IF ((elon1 - elonv).LT.(-180.))
     &     elon1l = elon1 + 360.
C
         elonl = elon
         IF ((elon  - elonv).GT.180.)
     &     elonl  = elon  - 360.
         IF ((elon - elonv).LT.(-180.))
     &     elonl = elon + 360.
C
         elonvr = elonv * radpd
C
C        RADIUS TO LOWER LEFT HAND (LL) CORNER
C
         ala1 =  alat1 * radpd
         rmll = rebydx * (((cosltn)**(1.-an))*(1.+an)**an) *
     &           (((cos(ala1))/(1.+h*sin(ala1)))**an)/an
C
C        USE LL POINT INFO TO LOCATE POLE POINT
C
         elo1 = elon1l * radpd
         arg = an * (elo1-elonvr)
         polei = 1. - h * rmll * sin(arg)
         polej = 1. + rmll * cos(arg)
C
C        RADIUS TO DESIRED POINT AND THE I J TOO
C
         ala =  alat * radpd
         rm = rebydx * ((cosltn**(1.-an))*(1.+an)**an) *
     &           (((cos(ala))/(1.+h*sin(ala)))**an)/an
C
         elo = elonl * radpd
         arg = an*(elo-elonvr)
         xi = polei + h * rm * sin(arg)
         xj = polej - rm * cos(arg)
C
C        IF COORDINATE LESS THAN 1
C        COMPENSATE FOR ORIGIN AT (1,1)
C
         IF (xi.LT.1.)  xi = xi - 1.
         IF (xj.LT.1.)  xj = xj - 1.
C
      RETURN
      END