w3fb07.f

Go to the documentation of this file.

C> @file
C> @brief Grid coords to lat/lon for grib.
C> @author John Stackpole @date 1988-01-01
 
C> Converts the coordinates of a location on earth given in a
C> grid coordinate system overlaid on a polar stereographic map pro-
C> jection true at 60 degrees n or s latitude to the
C> natural coordinate system of latitude/longitude
C> w3fb07() is the reverse of w3fb06().
C> uses grib specification of the location of the grid
C>
C> Program history log:
C> - John Stackpole 1988-01-01
C> - Ralph Jones 1990-04-12 Convert to cray cft77 fortran.
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] ALON1 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 60 deg lat
C> must be set negative if using
C> southern hemisphere projection; real*4
C> 190500.0 lfm grid,
C> 381000.0 nh pe grid, -381000.0 sh pe grid, etc.
C> @param[in] ALONV 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> for example:
C> 255.0 for lfm grid,
C> 280.0 nh pe grid, 100.0 sh pe grid, etc.
C> @param[out] ALAT Latitude in degrees (negative in southern hemi.).
C> @param[out] ALON East longitude in degrees, 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-01-01
      SUBROUTINE w3fb07(XI,XJ,ALAT1,ALON1,DX,ALONV,ALAT,ALON)
C
         DATA  rerth /6.3712e+6/,pi/3.1416/
         DATA  ss60  /1.86603/
C
C        PRELIMINARY VARIABLES AND REDIFINITIONS
C
C        H = 1 FOR NORTHERN HEMISPHERE; = -1 FOR SOUTHERN
C
C        REFLON IS LONGITUDE UPON WHICH THE POSITIVE X-COORDINATE
C        DRAWN THROUGH THE POLE AND TO THE RIGHT LIES
C        ROTATED AROUND FROM ORIENTATION (Y-COORDINATE) LONGITUDE
C        DIFFERENTLY IN EACH HEMISPHERE
C
         IF (dx.LT.0) THEN
           h      = -1.0
           dxl    = -dx
           reflon = alonv - 90.0
         ELSE
           h      = 1.0
           dxl    = dx
           reflon = alonv - 270.0
         ENDIF
C
         radpd  = pi    / 180.0
         degprd = 180.0 / pi
         rebydx = rerth / dxl
C
C        RADIUS TO LOWER LEFT HAND (LL) CORNER
C
         ala1 =  alat1 * radpd
         rmll = rebydx * cos(ala1) * ss60/(1. + h * sin(ala1))
C
C        USE LL POINT INFO TO LOCATE POLE POINT
C
         alo1 = (alon1 - reflon) * radpd
         polei = 1. - rmll * cos(alo1)
         polej = 1. - h * rmll * sin(alo1)
C
C        RADIUS TO THE I,J POINT (IN GRID UNITS)
C
         xx =  xi - polei
         yy = (xj - polej) * h
         r2 =  xx**2 + yy**2
C
C        NOW THE MAGIC FORMULAE
C
         IF (r2.EQ.0) THEN
           alat = h * 90.
           alon = reflon
         ELSE
           gi2    = (rebydx * ss60)**2
           alat   = degprd * h * asin((gi2 - r2)/(gi2 + r2))
           arccos = acos(xx/sqrt(r2))
           IF (yy.GT.0) THEN
             alon = reflon + degprd * arccos
           ELSE
             alon = reflon - degprd * arccos
           ENDIF
         ENDIF
         IF (alon.LT.0) alon = alon + 360.
C
      RETURN
      END