w3fb12.f

Go to the documentation of this file.

C> @file
C> @brief Lambert(i,j) to lat/lon for grib.
C> @author John Stackpole @date 1988-11-25
 
C> Converts the coordinates of a location on Earth given in a
C> grid coordinate system overlaid on a lambert conformal tangent
C> cone projection true at a given N or S latitude to the
C> natural coordinate system of latitude/longitude
C> w3fb12() is the reverse of w3fb11().
C> 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] XI I coordinate of the point real*4
C> @param[in] XJ J coordinate of the point real*4
C> @param[in] ALAT1 Latitude of lower left point of grid (point 1,1)
C> latitude <0 for southern hemisphere; real*4
C> @param[in] ELON1 Longitude of lower left point of grid (point 1,1)
C> east longitude used throughout; 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> the grid) along which latitude increases as
C> the y-coordinate increases. real*4
C> this is also the meridian (on the other 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> (touches or osculates) the spherical Earth.
C> set negative to indicate a
C> southern hemisphere projection; real*4
C>
C> @param[out] ALAT Latitude in degrees (negative in southern hemi.)
C> @param[out] ELON East longitude in degrees, real*4
C> @param[out] IERR
C> - .eq. 0 if no problem
C> - .ge. 1 if the requested xi,xj point is in the
C> forbidden zone, i.e. off the lambert map
C> in the open space where the cone is cut.
C> - if ierr.ge.1 then alat=999. and elon=999.
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 w3fb12(XI,XJ,ALAT1,ELON1,DX,ELONV,ALATAN,ALAT,ELON,
     &                               IERR)
C
         LOGICAL NEWMAP
C
         SAVE
C
         DATA  rerth /6.3712e+6/, pi/3.14159/, oldrml/99999./
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
         piby2  = pi     / 2.0
         radpd  = pi     / 180.0
         degprd = 1.0    / radpd
         rebydx = rerth  / dx
         alatn1 = alatan * radpd
         an     = h * sin(alatn1)
         cosltn = cos(alatn1)
C
C        MAKE SURE THAT INPUT LONGITUDE DOES 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
         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 RMLL TO TEST IF MAP AND GRID UNCHANGED FROM PREVIOUS
C        CALL TO THIS CODE.  THUS AVOID UNNEEDED RECOMPUTATIONS.
C
         IF (rmll.EQ.oldrml) THEN
           newmap = .false.
         ELSE
           newmap = .true.
           oldrml = rmll
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)
         ENDIF
C
C        RADIUS TO THE I,J POINT (IN GRID UNITS)
C              YY REVERSED SO POSITIVE IS DOWN
C
         xx = xi - polei
         yy = polej - xj
         r2 = xx**2 + yy**2
C
C        CHECK THAT THE REQUESTED I,J IS NOT IN THE FORBIDDEN ZONE
C           YY MUST BE POSITIVE UP FOR THIS TEST
C
         theta = pi*(1.-an)
         beta = abs(atan2(xx,-yy))
         ierr = 0
         IF (beta.LE.theta) THEN
           ierr = 1
           alat = 999.
           elon = 999.
           IF (.NOT.newmap)  RETURN
         ENDIF
C
C        NOW THE MAGIC FORMULAE
C
         IF (r2.EQ.0) THEN
           alat = h * 90.0
           elon = elonv
         ELSE
C
C          FIRST THE LONGITUDE
C
           elon = elonv + degprd * atan2(h*xx,yy)/an
           elon = amod(elon+360., 360.)
C
C          NOW THE LATITUDE
C          RECALCULATE THE THING ONLY IF MAP IS NEW SINCE LAST TIME
C
           IF (newmap) THEN
             aninv = 1./an
             aninv2 = aninv/2.
             thing = ((an/rebydx) ** aninv)/
     &         ((cosltn**((1.-an)*aninv))*(1.+ an))
           ENDIF
           alat = h*(piby2 - 2.*atan(thing*(r2**aninv2)))*degprd
         ENDIF
C
C        FOLLOWING TO ASSURE ERROR VALUES IF FIRST TIME THRU
C         IS OFF THE MAP
C
         IF (ierr.NE.0) THEN
           alat = 999.
           elon = 999.
           ierr = 2
         ENDIF
         RETURN
         END