var/obsproc/src/module_map_utils.F90

   1 !dis
   2 !dis    Open Source License/Disclaimer, Forecast Systems Laboratory
   3 !dis    NOAA/OAR/FSL, 325 Broadway Boulder, CO 80305
   4 !dis
   5 !dis    This software is distributed under the Open Source Definition,
   6 !dis    which may be found at http://www.opensource.org/osd.html.
   7 !dis
   8 !dis    In particular, redistribution and use in source and binary forms,
   9 !dis    with or without modification, are permitted provided that the
  10 !dis    following conditions are met:
  11 !dis
  12 !dis    - Redistributions of source code must retain this notice, this
  13 !dis    list of conditions and the following disclaimer.
  14 !dis
  15 !dis    - Redistributions in binary form must provide access to this
  16 !dis    notice, this list of conditions and the following disclaimer, and
  17 !dis    the underlying source code.
  18 !dis
  19 !dis    - All modifications to this software must be clearly documented,
  20 !dis    and are solely the responsibility of the agent making the
  21 !dis    modifications.
  22 !dis
  23 !dis    - If significant modifications or enhancements are made to this
  24 !dis    software, the FSL Software Policy Manager
  25 !dis    (softwaremgr@fsl.noaa.gov) should be notified.
  26 !dis
  27 !dis    THIS SOFTWARE AND ITS DOCUMENTATION ARE IN THE PUBLIC DOMAIN
  28 !dis    AND ARE FURNISHED "AS IS."  THE AUTHORS, THE UNITED STATES
  29 !dis    GOVERNMENT, ITS INSTRUMENTALITIES, OFFICERS, EMPLOYEES, AND
  30 !dis    AGENTS MAKE NO WARRANTY, EXPRESS OR IMPLIED, AS TO THE USEFULNESS
  31 !dis    OF THE SOFTWARE AND DOCUMENTATION FOR ANY PURPOSE.  THEY ASSUME
  32 !dis    NO RESPONSIBILITY (1) FOR THE USE OF THE SOFTWARE AND
  33 !dis    DOCUMENTATION; OR (2) TO PROVIDE TECHNICAL SUPPORT TO USERS.
  34 !dis
  35 !dis
  36
  37 MODULE map_utils
  38
  39 ! Module that defines constants, data structures, and
  40 ! subroutines used to convert grid indices to lat/lon
  41 ! and vice versa.
  42 !
  43 ! SUPPORTED PROJECTIONS
  44 ! ---------------------
  45 ! Cylindrical Lat/Lon (code = PROJ_LATLON)
  46 ! Mercator (code = PROJ_MERC)
  47 ! Lambert Conformal (code = PROJ_LC)
  48 ! Polar Stereographic (code = PROJ_PS)
  49 !
  50 ! REMARKS
  51 ! -------
  52 ! The routines contained within were adapted from routines
  53 ! obtained from NCEP's w3 library.  The original NCEP routines were less
  54 ! flexible (e.g., polar-stereo routines only supported truelat of 60N/60S)
  55 ! than what we needed, so modifications based on equations in Hoke, Hayes, and
  56 ! Renninger (AFGWC/TN/79-003) were added to improve the flexibility.
  57 ! Additionally, coding was improved to F90 standards and the routines were
  58 ! combined into this module.
  59 !
  60 ! ASSUMPTIONS
  61 ! -----------
  62 !  Grid Definition:
  63 !    For mercator, lambert conformal, and polar-stereographic projections,
  64 !    the routines within assume the following:
  65 !
  66 !       1.  Grid is dimensioned (i,j) where i is the East-West direction,
  67 !           positive toward the east, and j is the north-south direction,
  68 !           positive toward the north.
  69 !       2.  Origin is at (1,1) and is located at the southwest corner,
  70 !           regardless of hemispere.
  71 !       3.  Grid spacing (dx) is always positive.
  72 !       4.  Values of true latitudes must be positive for NH domains
  73 !           and negative for SH domains.
  74 !
  75 !     For the latlon projection, the grid origin may be at any of the
  76 !     corners, and the deltalat and deltalon values can be signed to
  77 !     account for this using the following convention:
  78 !       Origin Location        Deltalat Sign      Deltalon Sign
  79 !       ---------------        -------------      -------------
  80 !        SW Corner                  +                   +
  81 !        NE Corner                  -                   -
  82 !        NW Corner                  -                   +
  83 !        SE Corner                  +                   -
  84 !
  85 !  Data Definitions:
  86 !       1. Any arguments that are a latitude value are expressed in
  87 !          degrees north with a valid range of -90 -> 90
  88 !       2. Any arguments that are a longitude value are expressed in
  89 !          degrees east with a valid range of -180 -> 180.
  90 !       3. Distances are in meters and are always positive.
  91 !       4. The standard longitude (stdlon) is defined as the longitude
  92 !          line which is parallel to the grid's y-axis (j-direction), along
  93 !          which latitude increases (NOT the absolute value of latitude, but
  94 !          the actual latitude, such that latitude increases continuously
  95 !          from the south pole to the north pole) as j increases.
  96 !       5. One true latitude value is required for polar-stereographic and
  97 !          mercator projections, and defines at which latitude the
  98 !          grid spacing is true.  For lambert conformal, two true latitude
  99 !          values must be specified, but may be set equal to each other to
 100 !          specify a tangent projection instead of a secant projection.
 101 !
 102 ! USAGE
 103 ! -----
 104 ! To use the routines in this module, the calling routines must have the
 105 ! following statement at the beginning of its declaration block:
 106 !   USE map_utils
 107 !
 108 ! The use of the module not only provides access to the necessary routines,
 109 ! but also defines a structure of TYPE (proj_info) that can be used
 110 ! to declare a variable of the same type to hold your map projection
 111 ! information.  It also defines some integer parameters that contain
 112 ! the projection codes so one only has to use those variable names rather
 113 ! than remembering the acutal code when using them.  The basic steps are
 114 ! as follows:
 115 !
 116 !   1.  Ensure the "USE map_utils" is in your declarations.
 117 !   2.  Declare the projection information structure as type(proj_info):
 118 !         TYPE(proj_info) :: proj
 119 !   3.  Populate your structure by calling the map_set routine:
 120 !         CALL map_set(code,lat1,lon1,knowni,knownj,dx,stdlon,truelat1,truelat2,proj)
 121 !       where:
 122 !         code (input) = one of PROJ_LATLON, PROJ_MERC, PROJ_LC, or PROJ_PS
 123 !         lat1 (input) = Latitude of grid origin point (i,j)=(1,1)
 124 !                         (see assumptions!)
 125 !         lon1 (input) = Longitude of grid origin
 126 !         knowni (input) = origin point, x-location
 127 !         knownj (input) = origin point, y-location
 128 !         dx (input) = grid spacing in meters (ignored for LATLON projections)
 129 !         stdlon (input) = Standard longitude for PROJ_PS and PROJ_LC,
 130 !               deltalon (see assumptions) for PROJ_LATLON,
 131 !               ignored for PROJ_MERC
 132 !         truelat1 (input) = 1st true latitude for PROJ_PS, PROJ_LC, and
 133 !                PROJ_MERC, deltalat (see assumptions) for PROJ_LATLON
 134 !         truelat2 (input) = 2nd true latitude for PROJ_LC,
 135 !                ignored for all others.
 136 !         proj (output) = The structure of type (proj_info) that will be fully
 137 !                populated after this call
 138 !
 139 !   4.  Now that the proj structure is populated, you may call either
 140 !       of the following routines:
 141 !
 142 !       latlon_to_ij(proj, lat, lon, i, j)
 143 !       ij_to_latlon(proj, i, j, lat, lon)
 144 !
 145 !       It is incumbent upon the calling routine to determine whether or
 146 !       not the values returned are within your domain's bounds.  All values
 147 !       of i, j, lat, and lon are REAL values.
 148 !
 149 !
 150 ! REFERENCES
 151 ! ----------
 152 !  Hoke, Hayes, and Renninger, "Map Preojections and Grid Systems for
 153 !       Meteorological Applications." AFGWC/TN-79/003(Rev), Air Weather
 154 !       Service, 1985.
 155 !
 156 !  NCAR MM5v3 Modeling System, REGRIDDER program, module_first_guess_map.F
 157 !  NCEP routines w3fb06, w3fb07, w3fb08, w3fb09, w3fb11, w3fb12
 158 !
 159 ! HISTORY
 160 ! -------
 161 ! 27 Mar 2001 - Original Version
 162 !               Brent L. Shaw, NOAA/FSL (CSU/CIRA)
 163 !
 164 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 165
 166 IMPLICIT NONE
 167
 168   ! Define some private constants
 169   REAL, PRIVATE, PARAMETER    :: pi = 3.1415927
 170   REAL, PRIVATE, PARAMETER    :: deg_per_rad = 180./pi
 171   REAL, PRIVATE, PARAMETER    :: rad_per_deg = pi / 180.
 172
 173   ! Mean Earth Radius in m.  The value below is consistent
 174   ! with NCEP's routines and grids.
 175 ! REAL, PUBLIC , PARAMETER    :: earth_radius_m = 6371200. ! from Brent
 176   REAL, PUBLIC , PARAMETER    :: earth_radius_m = 6370000. ! same as MM5 system
 177   REAL, PUBLIC , PARAMETER    :: radians_per_degree = pi / 180.
 178
 179   ! Define public parameters
 180
 181   ! Projection codes for proj_info structure:
 182   INTEGER, PUBLIC, PARAMETER  :: PROJ_LATLON = 0
 183   INTEGER, PUBLIC, PARAMETER  :: PROJ_MERC = 1
 184   INTEGER, PUBLIC, PARAMETER  :: PROJ_LC = 3
 185   INTEGER, PUBLIC, PARAMETER  :: PROJ_PS = 5
 186
 187
 188   ! Define data structures to define various projections
 189
 190   TYPE proj_info
 191
 192     INTEGER          :: code     ! Integer code for projection type
 193     REAL             :: lat1    ! SW latitude (1,1) in degrees (-90->90N)
 194     REAL             :: lon1    ! SW longitude (1,1) in degrees (-180->180E)
 195     REAL             :: dx       ! Grid spacing in meters at truelats, used
 196                                  ! only for ps, lc, and merc projections
 197     REAL             :: dlat     ! Lat increment for lat/lon grids
 198     REAL             :: dlon     ! Lon increment for lat/lon grids
 199     REAL             :: stdlon   ! Longitude parallel to y-axis (-180->180E)
 200     REAL             :: truelat1 ! First true latitude (all projections)
 201     REAL             :: truelat2 ! Second true lat (LC only)
 202     REAL             :: hemi     ! 1 for NH, -1 for SH
 203     REAL             :: cone     ! Cone factor for LC projections
 204     REAL             :: polei    ! Computed i-location of pole point
 205     REAL             :: polej    ! Computed j-location of pole point
 206     REAL             :: rsw      ! Computed radius to SW corner
 207     REAL             :: rebydx   ! Earth radius divided by dx
 208     REAL             :: knowni   ! X-location of known lat/lon
 209     REAL             :: knownj   ! Y-location of known lat/lon
 210     LOGICAL          :: init     ! Flag to indicate if this struct is
 211                                  ! ready for use
 212   END TYPE proj_info
 213
 214 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 215 CONTAINS
 216 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 217   SUBROUTINE map_init(proj)
 218     ! Initializes the map projection structure to missing values
 219
 220     IMPLICIT NONE
 221     TYPE(proj_info), INTENT(INOUT)  :: proj
 222
 223     proj%lat1 =    -999.9
 224     proj%lon1 =    -999.9
 225     proj%dx    =    -999.9
 226     proj%stdlon =   -999.9
 227     proj%truelat1 = -999.9
 228     proj%truelat2 = -999.9
 229     proj%hemi     = 0.0
 230     proj%cone     = -999.9
 231     proj%polei    = -999.9
 232     proj%polej    = -999.9
 233     proj%rsw      = -999.9
 234     proj%knowni   = -999.9
 235     proj%knownj   = -999.9
 236     proj%init     = .FALSE.
 237
 238   END SUBROUTINE map_init
 239 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 240   SUBROUTINE map_set(proj_code,lat1,lon1,knowni,knownj,dx,stdlon,truelat1,truelat2,proj)
 241     ! Given a partially filled proj_info structure, this routine computes
 242     ! polei, polej, rsw, and cone (if LC projection) to complete the
 243     ! structure.  This allows us to eliminate redundant calculations when
 244     ! calling the coordinate conversion routines multiple times for the
 245     ! same map.
 246     ! This will generally be the first routine called when a user wants
 247     ! to be able to use the coordinate conversion routines, and it
 248     ! will call the appropriate subroutines based on the
 249     ! proj%code which indicates which projection type  this is.
 250
 251     IMPLICIT NONE
 252
 253     ! Declare arguments
 254     INTEGER, INTENT(IN)               :: proj_code
 255     REAL, INTENT(IN)                  :: lat1
 256     REAL, INTENT(IN)                  :: lon1
 257     REAL, INTENT(IN)                  :: dx
 258     REAL, INTENT(IN)                  :: stdlon
 259     REAL, INTENT(IN)                  :: truelat1
 260     REAL, INTENT(IN)                  :: truelat2
 261     REAL, INTENT(IN)                  :: knowni , knownj
 262     TYPE(proj_info), INTENT(OUT)      :: proj
 263
 264     ! Local variables
 265
 266
 267     ! Executable code
 268
 269     ! First, check for validity of mandatory variables in proj
 270     IF ( ABS(lat1) .GT. 90. ) THEN
 271       WRITE(0,'(A)') 'Latitude of origin corner required as follows:'
 272       WRITE(0,'(A)') '    -90N <= lat1 < = 90.N'
 273       STOP 'MAP_INIT'
 274     ENDIF
 275     IF ( ABS(lon1) .GT. 180.) THEN
 276       WRITE(0,'(A)') 'Longitude of origin required as follows:'
 277       WRITE(0,'(A)') '   -180E <= lon1 <= 180W'
 278       STOP 'MAP_INIT'
 279     ENDIF
 280     IF ((dx .LE. 0.).AND.(proj_code .NE. PROJ_LATLON)) THEN
 281       WRITE(0,'(A)') 'Require grid spacing (dx) in meters be positive!'
 282       STOP 'MAP_INIT'
 283     ENDIF
 284     IF ((ABS(stdlon) .GT. 180.).AND.(proj_code .NE. PROJ_MERC)) THEN
 285       WRITE(0,'(A)') 'Need orientation longitude (stdlon) as: '
 286       WRITE(0,'(A)') '   -180E <= lon1 <= 180W'
 287       STOP 'MAP_INIT'
 288     ENDIF
 289     IF (ABS(truelat1).GT.90.) THEN
 290       WRITE(0,'(A)') 'Set true latitude 1 for all projections!'
 291       STOP 'MAP_INIT'
 292     ENDIF
 293
 294     CALL map_init(proj)
 295     proj%code  = proj_code
 296     proj%lat1 = lat1
 297     proj%lon1 = lon1
 298     proj%knowni = knowni
 299     proj%knownj = knownj
 300     proj%dx    = dx
 301     proj%stdlon = stdlon
 302     proj%truelat1 = truelat1
 303     proj%truelat2 = truelat2
 304     IF (proj%code .NE. PROJ_LATLON) THEN
 305       proj%dx = dx
 306       IF (truelat1 .LT. 0.) THEN
 307         proj%hemi = -1.0
 308       ELSE
 309         proj%hemi = 1.0
 310       ENDIF
 311       proj%rebydx = earth_radius_m / dx
 312     ENDIF
 313     pick_proj: SELECT CASE(proj%code)
 314
 315       CASE(PROJ_PS)
 316         WRITE(0,'(A)') 'Setting up POLAR STEREOGRAPHIC map...'
 317         CALL set_ps(proj)
 318
 319       CASE(PROJ_LC)
 320         WRITE(0,'(A)') 'Setting up LAMBERT CONFORMAL map...'
 321         IF (ABS(proj%truelat2) .GT. 90.) THEN
 322           WRITE(0,'(A)') 'Second true latitude not set, assuming a tangent'
 323           WRITE(0,'(A,F10.3)') 'projection at truelat1: ', proj%truelat1
 324           proj%truelat2=proj%truelat1
 325         ENDIF
 326         CALL set_lc(proj)
 327
 328       CASE (PROJ_MERC)
 329         WRITE(0,'(A)') 'Setting up MERCATOR map...'
 330         CALL set_merc(proj)
 331
 332       CASE (PROJ_LATLON)
 333         WRITE(0,'(A)') 'Setting up CYLINDRICAL EQUIDISTANT LATLON map...'
 334         ! Convert lon1 to 0->360 notation
 335         IF (proj%lon1 .LT. 0.) proj%lon1 = proj%lon1 + 360.
 336
 337       CASE DEFAULT
 338         WRITE(0,'(A,I2)') 'Unknown projection code: ', proj%code
 339         STOP 'MAP_INIT'
 340
 341     END SELECT pick_proj
 342     proj%init = .TRUE.
 343     RETURN
 344   END SUBROUTINE map_set
 345 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 346   SUBROUTINE latlon_to_ij(proj, lat, lon, i, j)
 347     ! Converts input lat/lon values to the cartesian (i,j) value
 348     ! for the given projection.
 349
 350     IMPLICIT NONE
 351     TYPE(proj_info), INTENT(IN)          :: proj
 352     REAL, INTENT(IN)                     :: lat
 353     REAL, INTENT(IN)                     :: lon
 354     REAL, INTENT(OUT)                    :: i
 355     REAL, INTENT(OUT)                    :: j
 356
 357     IF (.NOT.proj%init) THEN
 358       WRITE(0,'(A)') 'You have not called map_set for this projection!'
 359       STOP 'LATLON_TO_IJ'
 360     ENDIF
 361
 362     SELECT CASE(proj%code)
 363
 364       CASE(PROJ_LATLON)
 365         CALL llij_latlon(lat,lon,proj,i,j)
 366
 367       CASE(PROJ_MERC)
 368         CALL llij_merc(lat,lon,proj,i,j)
 369         i = i + proj%knowni - 1.0
 370         j = j + proj%knownj - 1.0
 371
 372       CASE(PROJ_PS)
 373         CALL llij_ps(lat,lon,proj,i,j)
 374
 375       CASE(PROJ_LC)
 376         CALL llij_lc(lat,lon,proj,i,j)
 377         i = i + proj%knowni - 1.0
 378         j = j + proj%knownj - 1.0
 379
 380       CASE DEFAULT
 381         WRITE(0,'(A,I2)') 'Unrecognized map projection code: ', proj%code
 382         STOP 'LATLON_TO_IJ'
 383
 384     END SELECT
 385
 386     RETURN
 387   END SUBROUTINE latlon_to_ij
 388 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 389   SUBROUTINE ij_to_latlon(proj, ii, jj, lat, lon)
 390     ! Computes geographical latitude and longitude for a given (i,j) point
 391     ! in a grid with a projection of proj
 392
 393     IMPLICIT NONE
 394     TYPE(proj_info),INTENT(IN)          :: proj
 395     REAL, INTENT(IN)                    :: ii
 396     REAL, INTENT(IN)                    :: jj
 397     REAL, INTENT(OUT)                   :: lat
 398     REAL, INTENT(OUT)                   :: lon
 399     REAL         :: i, j
 400
 401     IF (.NOT.proj%init) THEN
 402       WRITE(0,'(A)') 'You have not called map_set for this projection!'
 403       STOP 'IJ_TO_LATLON'
 404     ENDIF
 405
 406     i = ii
 407     j = jj
 408
 409     SELECT CASE (proj%code)
 410
 411       CASE (PROJ_LATLON)
 412         CALL ijll_latlon(i, j, proj, lat, lon)
 413
 414       CASE (PROJ_MERC)
 415         i = ii - proj%knowni + 1.0
 416         j = jj - proj%knownj + 1.0
 417         CALL ijll_merc(i, j, proj, lat, lon)
 418
 419       CASE (PROJ_PS)
 420         CALL ijll_ps(i, j, proj, lat, lon)
 421
 422       CASE (PROJ_LC)
 423
 424         i = ii - proj%knowni + 1.0
 425         j = jj - proj%knownj + 1.0
 426         CALL ijll_lc(i, j, proj, lat, lon)
 427
 428       CASE DEFAULT
 429         WRITE(0,'(A,I2)') 'Unrecognized map projection code: ', proj%code
 430         STOP 'IJ_TO_LATLON'
 431
 432     END SELECT
 433     RETURN
 434   END SUBROUTINE ij_to_latlon
 435 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 436   SUBROUTINE set_ps(proj)
 437     ! Initializes a polar-stereographic map projection from the partially
 438     ! filled proj structure. This routine computes the radius to the
 439     ! southwest corner and computes the i/j location of the pole for use
 440     ! in llij_ps and ijll_ps.
 441     IMPLICIT NONE
 442
 443     ! Declare args
 444     TYPE(proj_info), INTENT(INOUT)    :: proj
 445
 446     ! Local vars
 447     REAL                              :: ala1
 448     REAL                              :: alo1
 449     REAL                              :: reflon
 450     REAL                              :: scale_top
 451
 452     ! Executable code
 453     reflon = proj%stdlon + 90.
 454
 455     ! Compute numerator term of map scale factor
 456     scale_top = 1. + proj%hemi * SIN(proj%truelat1 * rad_per_deg)
 457
 458     ! Compute radius to lower-left (SW) corner
 459     ala1 = proj%lat1 * rad_per_deg
 460     proj%rsw = proj%rebydx*COS(ala1)*scale_top/(1.+proj%hemi*SIN(ala1))
 461
 462     ! Find the pole point
 463     alo1 = (proj%lon1 - reflon) * rad_per_deg
 464     proj%polei = proj%knowni - proj%rsw * COS(alo1)
 465     proj%polej = proj%knownj - proj%hemi * proj%rsw * SIN(alo1)
 466     WRITE(0,'(A,2F10.1)')'Computed (I,J) of pole point: ',proj%polei,proj%polej
 467     RETURN
 468   END SUBROUTINE set_ps
 469 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 470   SUBROUTINE llij_ps(lat,lon,proj,i,j)
 471     ! Given latitude (-90 to 90), longitude (-180 to 180), and the
 472     ! standard polar-stereographic projection information via the
 473     ! public proj structure, this routine returns the i/j indices which
 474     ! if within the domain range from 1->nx and 1->ny, respectively.
 475
 476     IMPLICIT NONE
 477
 478     ! Delcare input arguments
 479     REAL, INTENT(IN)               :: lat
 480     REAL, INTENT(IN)               :: lon
 481     TYPE(proj_info),INTENT(IN)     :: proj
 482
 483     ! Declare output arguments
 484     REAL, INTENT(OUT)              :: i !(x-index)
 485     REAL, INTENT(OUT)              :: j !(y-index)
 486
 487     ! Declare local variables
 488
 489     REAL                           :: reflon
 490     REAL                           :: scale_top
 491     REAL                           :: ala
 492     REAL                           :: alo
 493     REAL                           :: rm
 494
 495     ! BEGIN CODE
 496
 497
 498     reflon = proj%stdlon + 90.
 499
 500     ! Compute numerator term of map scale factor
 501
 502     scale_top = 1. + proj%hemi * SIN(proj%truelat1 * rad_per_deg)
 503
 504     ! Find radius to desired point
 505     ala = lat * rad_per_deg
 506     rm = proj%rebydx * COS(ala) * scale_top/(1. + proj%hemi *SIN(ala))
 507     alo = (lon - reflon) * rad_per_deg
 508     i = proj%polei + rm * COS(alo)
 509     j = proj%polej + proj%hemi * rm * SIN(alo)
 510
 511     RETURN
 512   END SUBROUTINE llij_ps
 513 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 514   SUBROUTINE ijll_ps(i, j, proj, lat, lon)
 515
 516     ! This is the inverse subroutine of llij_ps.  It returns the
 517     ! latitude and longitude of an i/j point given the projection info
 518     ! structure.
 519
 520     IMPLICIT NONE
 521
 522     ! Declare input arguments
 523     REAL, INTENT(IN)                    :: i    ! Column
 524     REAL, INTENT(IN)                    :: j    ! Row
 525     TYPE (proj_info), INTENT(IN)        :: proj
 526
 527     ! Declare output arguments
 528     REAL, INTENT(OUT)                   :: lat     ! -90 -> 90 North
 529     REAL, INTENT(OUT)                   :: lon     ! -180 -> 180 East
 530
 531     ! Local variables
 532     REAL                                :: reflon
 533     REAL                                :: scale_top
 534     REAL                                :: xx,yy
 535     REAL                                :: gi2, r2
 536     REAL                                :: arccos
 537
 538     ! Begin Code
 539
 540     ! Compute the reference longitude by rotating 90 degrees to the east
 541     ! to find the longitude line parallel to the positive x-axis.
 542     reflon = proj%stdlon + 90.
 543
 544     ! Compute numerator term of map scale factor
 545     scale_top = 1. + proj%hemi * SIN(proj%truelat1 * rad_per_deg)
 546
 547     ! Compute radius to point of interest
 548     xx = i - proj%polei
 549     yy = (j - proj%polej) * proj%hemi
 550     r2 = xx**2 + yy**2
 551
 552     ! Now the magic code
 553     IF (r2 .EQ. 0.) THEN
 554       lat = proj%hemi * 90.
 555       lon = reflon
 556     ELSE
 557       gi2 = (proj%rebydx * scale_top)**2.
 558       lat = deg_per_rad * proj%hemi * ASIN((gi2-r2)/(gi2+r2))
 559       arccos = ACOS(xx/SQRT(r2))
 560       IF (yy .GT. 0) THEN
 561         lon = reflon + deg_per_rad * arccos
 562       ELSE
 563         lon = reflon - deg_per_rad * arccos
 564       ENDIF
 565     ENDIF
 566
 567     ! Convert to a -180 -> 180 East convention
 568     IF (lon .GT. 180.) lon = lon - 360.
 569     IF (lon .LT. -180.) lon = lon + 360.
 570     RETURN
 571
 572   END SUBROUTINE ijll_ps
 573 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 574   SUBROUTINE set_lc(proj)
 575     ! Initialize the remaining items in the proj structure for a
 576     ! lambert conformal grid.
 577
 578     IMPLICIT NONE
 579
 580     TYPE(proj_info), INTENT(INOUT)     :: proj
 581
 582     REAL                               :: arg
 583     REAL                               :: deltalon1
 584     REAL                               :: tl1r
 585     REAL                               :: ctl1r
 586
 587     ! Compute cone factor
 588     CALL lc_cone(proj%truelat1, proj%truelat2, proj%cone)
 589     WRITE(0,'(A,F8.6)') 'Computed cone factor: ', proj%cone
 590     ! Compute longitude differences and ensure we stay out of the
 591     ! forbidden "cut zone"
 592     deltalon1 = proj%lon1 - proj%stdlon
 593     IF (deltalon1 .GT. +180.) deltalon1 = deltalon1 - 360.
 594     IF (deltalon1 .LT. -180.) deltalon1 = deltalon1 + 360.
 595
 596     ! Convert truelat1 to radian and compute COS for later use
 597     tl1r = proj%truelat1 * rad_per_deg
 598     ctl1r = COS(tl1r)
 599
 600     ! Compute the radius to our known lower-left (SW) corner
 601     proj%rsw = proj%rebydx * ctl1r/proj%cone * &
 602            (TAN((90.*proj%hemi-proj%lat1)*rad_per_deg/2.) / &
 603             TAN((90.*proj%hemi-proj%truelat1)*rad_per_deg/2.))**proj%cone
 604
 605     ! Find pole point
 606     arg = proj%cone*(deltalon1*rad_per_deg)
 607     proj%polei = 1. - proj%hemi * proj%rsw * SIN(arg)
 608     proj%polej = 1. + proj%rsw * COS(arg)
 609     WRITE(0,'(A,2F10.3)') 'Computed pole i/j = ', proj%polei, proj%polej
 610     RETURN
 611   END SUBROUTINE set_lc
 612 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 613   SUBROUTINE lc_cone(truelat1, truelat2, cone)
 614
 615   ! Subroutine to compute the cone factor of a Lambert Conformal projection
 616
 617     IMPLICIT NONE
 618
 619     ! Input Args
 620     REAL, INTENT(IN)             :: truelat1  ! (-90 -> 90 degrees N)
 621     REAL, INTENT(IN)             :: truelat2  !   "   "  "   "     "
 622
 623     ! Output Args
 624     REAL, INTENT(OUT)            :: cone
 625
 626     ! Locals
 627
 628     ! BEGIN CODE
 629
 630     ! First, see if this is a secant or tangent projection.  For tangent
 631     ! projections, truelat1 = truelat2 and the cone is tangent to the
 632     ! Earth's surface at this latitude.  For secant projections, the cone
 633     ! intersects the Earth's surface at each of the distinctly different
 634     ! latitudes
 635     IF (ABS(truelat1-truelat2) .GT. 0.1) THEN
 636       cone = ALOG10(COS(truelat1*rad_per_deg)) - &
 637              ALOG10(COS(truelat2*rad_per_deg))
 638       cone = cone /(ALOG10(TAN((45.0 - ABS(truelat1)/2.0) * rad_per_deg)) - &
 639              ALOG10(TAN((45.0 - ABS(truelat2)/2.0) * rad_per_deg)))
 640     ELSE
 641        cone = SIN(ABS(truelat1)*rad_per_deg )
 642     ENDIF
 643   RETURN
 644   END SUBROUTINE lc_cone
 645 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 646   SUBROUTINE ijll_lc( i, j, proj, lat, lon)
 647
 648   ! Subroutine to convert from the (i,j) cartesian coordinate to the
 649   ! geographical latitude and longitude for a Lambert Conformal projection.
 650
 651   ! History:
 652   ! 25 Jul 01: Corrected by B. Shaw, NOAA/FSL
 653   !
 654     IMPLICIT NONE
 655
 656     ! Input Args
 657     REAL, INTENT(IN)              :: i        ! Cartesian X coordinate
 658     REAL, INTENT(IN)              :: j        ! Cartesian Y coordinate
 659     TYPE(proj_info),INTENT(IN)    :: proj     ! Projection info structure
 660
 661     ! Output Args
 662     REAL, INTENT(OUT)             :: lat      ! Latitude (-90->90 deg N)
 663     REAL, INTENT(OUT)             :: lon      ! Longitude (-180->180 E)
 664
 665     ! Locals
 666     REAL                          :: inew
 667     REAL                          :: jnew
 668     REAL                          :: r
 669     REAL                          :: chi,chi1,chi2
 670     REAL                          :: r2
 671     REAL                          :: xx
 672     REAL                          :: yy
 673
 674     ! BEGIN CODE
 675
 676     chi1 = (90. - proj%hemi*proj%truelat1)*rad_per_deg
 677     chi2 = (90. - proj%hemi*proj%truelat2)*rad_per_deg
 678
 679     ! See if we are in the southern hemispere and flip the indices
 680     ! if we are.
 681     IF (proj%hemi .EQ. -1.) THEN
 682       inew = -i + 2.
 683       jnew = -j + 2.
 684     ELSE
 685       inew = i
 686       jnew = j
 687     ENDIF
 688
 689     ! Compute radius**2 to i/j location
 690     xx = inew - proj%polei
 691     yy = proj%polej - jnew
 692     r2 = (xx*xx + yy*yy)
 693     r = SQRT(r2)/proj%rebydx
 694
 695     ! Convert to lat/lon
 696     IF (r2 .EQ. 0.) THEN
 697       lat = proj%hemi * 90.
 698       lon = proj%stdlon
 699     ELSE
 700
 701       ! Longitude
 702       lon = proj%stdlon + deg_per_rad * ATAN2(proj%hemi*xx,yy)/proj%cone
 703       lon = AMOD(lon+360., 360.)
 704
 705       ! Latitude.  Latitude determined by solving an equation adapted
 706       ! from:
 707       !  Maling, D.H., 1973: Coordinate Systems and Map Projections
 708       ! Equations #20 in Appendix I.
 709
 710       IF (chi1 .EQ. chi2) THEN
 711         chi = 2.0*ATAN( ( r/TAN(chi1) )**(1./proj%cone) * TAN(chi1*0.5) )
 712       ELSE
 713         chi = 2.0*ATAN( (r*proj%cone/SIN(chi1))**(1./proj%cone) * TAN(chi1*0.5))
 714       ENDIF
 715       lat = (90.0-chi*deg_per_rad)*proj%hemi
 716
 717     ENDIF
 718
 719     IF (lon .GT. +180.) lon = lon - 360.
 720     IF (lon .LT. -180.) lon = lon + 360.
 721     RETURN
 722     END SUBROUTINE ijll_lc
 723 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 724   SUBROUTINE llij_lc( lat, lon, proj, i, j)
 725
 726   ! Subroutine to compute the geographical latitude and longitude values
 727   ! to the cartesian x/y on a Lambert Conformal projection.
 728
 729     IMPLICIT NONE
 730
 731     ! Input Args
 732     REAL, INTENT(IN)              :: lat      ! Latitude (-90->90 deg N)
 733     REAL, INTENT(IN)              :: lon      ! Longitude (-180->180 E)
 734     TYPE(proj_info),INTENT(IN)      :: proj     ! Projection info structure
 735
 736     ! Output Args
 737     REAL, INTENT(OUT)             :: i        ! Cartesian X coordinate
 738     REAL, INTENT(OUT)             :: j        ! Cartesian Y coordinate
 739
 740     ! Locals
 741     REAL                          :: arg
 742     REAL                          :: deltalon
 743     REAL                          :: tl1r
 744     REAL                          :: rm
 745     REAL                          :: ctl1r
 746
 747
 748     ! BEGIN CODE
 749
 750     ! Compute deltalon between known longitude and standard lon and ensure
 751     ! it is not in the cut zone
 752     deltalon = lon - proj%stdlon
 753     IF (deltalon .GT. +180.) deltalon = deltalon - 360.
 754     IF (deltalon .LT. -180.) deltalon = deltalon + 360.
 755
 756     ! Convert truelat1 to radian and compute COS for later use
 757     tl1r = proj%truelat1 * rad_per_deg
 758     ctl1r = COS(tl1r)
 759
 760     ! Radius to desired point
 761     rm = proj%rebydx * ctl1r/proj%cone * &
 762          (TAN((90.*proj%hemi-lat)*rad_per_deg/2.) / &
 763           TAN((90.*proj%hemi-proj%truelat1)*rad_per_deg/2.))**proj%cone
 764
 765     arg = proj%cone*(deltalon*rad_per_deg)
 766     i = proj%polei + proj%hemi * rm * SIN(arg)
 767     j = proj%polej - rm * COS(arg)
 768
 769     ! Finally, if we are in the southern hemisphere, flip the i/j
 770     ! values to a coordinate system where (1,1) is the SW corner
 771     ! (what we assume) which is different than the original NCEP
 772     ! algorithms which used the NE corner as the origin in the
 773     ! southern hemisphere (left-hand vs. right-hand coordinate?)
 774     IF (proj%hemi .EQ. -1.) THEN
 775       i = 2. - i
 776       j = 2. - j
 777     ENDIF
 778     RETURN
 779   END SUBROUTINE llij_lc
 780 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 781   SUBROUTINE set_merc(proj)
 782
 783     ! Sets up the remaining basic elements for the mercator projection
 784
 785     IMPLICIT NONE
 786     TYPE(proj_info), INTENT(INOUT)       :: proj
 787     REAL                                 :: clain
 788
 789
 790     !  Preliminary variables
 791
 792     clain = COS(rad_per_deg*proj%truelat1)
 793     proj%dlon = proj%dx / (earth_radius_m * clain)
 794
 795     ! Compute distance from equator to origin, and store in the
 796     ! proj%rsw tag.
 797
 798     proj%rsw = 0.
 799     IF (proj%lat1 .NE. 0.) THEN
 800       proj%rsw = (ALOG(TAN(0.5*((proj%lat1+90.)*rad_per_deg))))/proj%dlon
 801     ENDIF
 802     RETURN
 803   END SUBROUTINE set_merc
 804 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 805   SUBROUTINE llij_merc(lat, lon, proj, i, j)
 806
 807     ! Compute i/j coordinate from lat lon for mercator projection
 808
 809     IMPLICIT NONE
 810     REAL, INTENT(IN)              :: lat
 811     REAL, INTENT(IN)              :: lon
 812     TYPE(proj_info),INTENT(IN)    :: proj
 813     REAL,INTENT(OUT)              :: i
 814     REAL,INTENT(OUT)              :: j
 815     REAL                          :: deltalon
 816
 817     deltalon = lon - proj%lon1
 818     IF (deltalon .LT. -180.) deltalon = deltalon + 360.
 819     IF (deltalon .GT. 180.) deltalon = deltalon - 360.
 820     i = 1. + (deltalon/(proj%dlon*deg_per_rad))
 821     j = 1. + (ALOG(TAN(0.5*((lat + 90.) * rad_per_deg)))) / &
 822            proj%dlon - proj%rsw
 823
 824     RETURN
 825   END SUBROUTINE llij_merc
 826 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 827   SUBROUTINE ijll_merc(i, j, proj, lat, lon)
 828
 829     ! Compute the lat/lon from i/j for mercator projection
 830
 831     IMPLICIT NONE
 832     REAL,INTENT(IN)               :: i
 833     REAL,INTENT(IN)               :: j
 834     TYPE(proj_info),INTENT(IN)    :: proj
 835     REAL, INTENT(OUT)             :: lat
 836     REAL, INTENT(OUT)             :: lon
 837
 838
 839     lat = 2.0*ATAN(EXP(proj%dlon*(proj%rsw + j-1.)))*deg_per_rad - 90.
 840     lon = (i-1.)*proj%dlon*deg_per_rad + proj%lon1
 841     IF (lon.GT.180.) lon = lon - 360.
 842     IF (lon.LT.-180.) lon = lon + 360.
 843     RETURN
 844   END SUBROUTINE ijll_merc
 845 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 846   SUBROUTINE llij_latlon(lat, lon, proj, i, j)
 847
 848     ! Compute the i/j location of a lat/lon on a LATLON grid.
 849     IMPLICIT NONE
 850     REAL, INTENT(IN)             :: lat
 851     REAL, INTENT(IN)             :: lon
 852     TYPE(proj_info), INTENT(IN)  :: proj
 853     REAL, INTENT(OUT)            :: i
 854     REAL, INTENT(OUT)            :: j
 855
 856     REAL                         :: deltalat
 857     REAL                         :: deltalon
 858     REAL                         :: lon360
 859     REAL                         :: latinc
 860     REAL                         :: loninc
 861
 862     ! Extract the latitude and longitude increments for this grid
 863     ! (e.g., 2.5 deg for NCEP reanalysis) from the proj structure, where
 864     ! loninc is saved in the stdlon tag and latinc is saved in truelat1
 865
 866     latinc = proj%truelat1
 867     loninc = proj%stdlon
 868
 869     ! Compute deltalat and deltalon as the difference between the input
 870     ! lat/lon and the origin lat/lon
 871
 872     deltalat = lat - proj%lat1
 873
 874     ! To account for issues around the dateline, convert the incoming
 875     ! longitudes to be 0->360.
 876     IF (lon .LT. 0) THEN
 877       lon360 = lon + 360.
 878     ELSE
 879       lon360 = lon
 880     ENDIF
 881     deltalon = lon360 - proj%lon1
 882
 883     ! Compute i/j
 884     i = deltalon/loninc + 1.
 885     j = deltalat/latinc + 1.
 886     RETURN
 887     END SUBROUTINE llij_latlon
 888 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 889   SUBROUTINE ijll_latlon(i, j, proj, lat, lon)
 890
 891     ! Compute the lat/lon location of an i/j on a LATLON grid.
 892     IMPLICIT NONE
 893     REAL, INTENT(IN)             :: i
 894     REAL, INTENT(IN)             :: j
 895     TYPE(proj_info), INTENT(IN)  :: proj
 896     REAL, INTENT(OUT)            :: lat
 897     REAL, INTENT(OUT)            :: lon
 898
 899     REAL                         :: deltalat
 900     REAL                         :: deltalon
 901     REAL                         :: lon360
 902     REAL                         :: latinc
 903     REAL                         :: loninc
 904
 905     ! Extract the latitude and longitude increments for this grid
 906     ! (e.g., 2.5 deg for NCEP reanalysis) from the proj structure, where
 907     ! loninc is saved in the stdlon tag and latinc is saved in truelat1
 908
 909     latinc = proj%truelat1
 910     loninc = proj%stdlon
 911
 912     ! Compute deltalat and deltalon
 913
 914     deltalat = (j-1.)*latinc
 915     deltalon = (i-1.)*loninc
 916     lat = proj%lat1 + deltalat
 917     lon = proj%lon1 + deltalon
 918
 919     IF ((ABS(lat) .GT. 90.).OR.(ABS(deltalon) .GT.360.)) THEN
 920       ! Off the earth for this grid
 921       lat = -999.
 922       lon = -999.
 923     ELSE
 924       lon = lon + 360.
 925       lon = AMOD(lon,360.)
 926       IF (lon .GT. 180.) lon = lon -360.
 927     ENDIF
 928
 929     RETURN
 930   END SUBROUTINE ijll_latlon
 931
 932 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 933 END MODULE map_utils