w3fb08.f

Go to the documentation of this file.

C> @file
C> @brief Lat/lon to merc (i,j) for grib.
C> @author John Stackpole @date 1988-03-01
 
C> Converts a location on earth given in
C> the coordinate system of latitude/longitude to an (i,j)
C> coordinate system overlaid on a mercator map projection
C> w3fb08() is the reverse of w3fb09()
C> uses grib specification of the location of the grid.
C>
C> Program history log:
C> - John Stackpole 1988-03-01
C> - Ralph Jones 1990-04-12 Convert to cray cft77 fortran.
C>
C> @param[in] ALAT Latitude in degrees (negative in southern hemis).
C> @param[in] ALON East longitude in degrees, real*4.
C> @param[in] ALAT1 Latitude of lower left corner of grid (point (1,1)).
C> @param[in] ALON1 Longitude of lower left corner of grid (point (1,1))
C> all real*4.
C> @param[in] ALATIN The latitude at which the mercator cylinder
C> intersects the earth.
C> @param[in] DX Mesh length of grid in meters at alatin.
C> @param[out] XI I coordinate of the point specified by alat, alon.
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-03-01
      SUBROUTINE w3fb08(ALAT,ALON,ALAT1,ALON1,ALATIN,DX,XI,XJ)
C
         DATA  rerth /6.3712e+6/, pi/3.1416/
C
C        PRELIMINARY VARIABLES AND REDIFINITIONS
C
         radpd  = pi    / 180.0
         degpr  = 180.0 / pi
         clain  = cos(radpd*alatin)
         dellon = dx   / (rerth*clain)
C
C        GET DISTANCE FROM EQUATOR TO ORIGIN ALAT1
C
         djeo = 0.
         IF (alat1.NE.0.)
     &     djeo = (alog(tan(0.5*((alat1+90.0)*radpd))))/dellon
C
C        NOW THE I AND J COORDINATES
C
         xi = 1. + ((alon - alon1)/(dellon*degpr))
         xj = 1. + (alog(tan(0.5*((alat + 90.) * radpd))))/
     &               dellon
     &             - djeo
C
      RETURN
      END