w3fb05.f

Go to the documentation of this file.

C> @file
C> @brief Grid coordinates to latitude, longitude.
C> @author Ralph Jones @date 1986-07-17
C>
C> Converts the coordinates of a location from the grid(i,j)
C> coordinate system overlaid on the polar stereographic map projec-
C> tion true at 60 degrees n or s latitude to the natural coordinate
C> system of latitude/longitude on the earth. w3fb05() is the reverse
C> of w3fb04().
C>
C> Program history log:
C> - Ralph Jones 1986-07-17
C> - Ralph Jones 1989-11-01 Change to cray cft77 fortran.
C>
C> @param[in] XI I of the point relative to the north or s. pole
C> @param[in] XJ J of the point relative to the north or s. pole
C> @param[in] XMESHL Mesh length of grid in km at 60 degrees(<0 if sh)
C> (190.5 lfm grid, 381.0 nh pe grid,-381.0 sh pe grid)
C> @param[in] ORIENT Orientation west longitude of the grid
C> (105.0 lfm grid, 80.0 nh pe grid, 260.0 sh pe grid)
C> @param[out] ALAT Latitude in degrees  (<0 if sh)
C> @param[out] ALONG West longitude in degrees
C>
C> @note All parameters in the calling statement must be
C> real. the range of allowable latitudes is from a pole to
C> 30 degrees into the opposite hemisphere.
C> the grid used in this subroutine has its origin (i=0,j=0)
C> at the pole, so if the user's grid has its origin at a point
C> other than a pole, a translation is required to get i and j for
C> input into w3fb05(). the subroutine grid is oriented so that
C> gridlines of i=constant are parallel to a west longitude sup-
C> plied by the user. the earth's radius is taken to be 6371.2 km.
C>
C> @note This code will not vectorize, it is normaly used in a
C> double do loop with w3ft01(), w3ft00(), etc. to vectorize it,
C> put it in line, put w3ft01(), w3ft00(), etc. in line.
C>
C> @author Ralph Jones @date 1986-07-17
      SUBROUTINE w3fb05(XI,XJ,XMESHL,ORIENT,ALAT,ALONG)
C
      DATA  degprd/57.2957795/
      DATA  earthr/6371.2/
C
      gi2   = ((1.86603 * earthr) / (xmeshl))**2
      r2    = xi * xi + xj * xj
C
      IF (r2.EQ.0.0) THEN
        along = 0.0
        alat  = 90.0
        IF (xmeshl.LT.0.0) alat = -alat
        RETURN
      ELSE
        alat  = asin((gi2 - r2) / (gi2 + r2)) * degprd
        angle = degprd * atan2(xj,xi)
        IF (angle.LT.0.0) angle = angle + 360.0
      ENDIF
C
      IF (xmeshl.GE.0.0) THEN
        along = 270.0 + orient - angle
C
      ELSE
C
        along = angle + orient - 270.0
        alat  = -(alat)
      ENDIF
C
      IF (along.LT.0.0)   along = along + 360.0
      IF (along.GE.360.0) along = along - 360.0
C
      RETURN
C
      END