WPS/geogrid/src/module_map_utils.f90

   1 MODULE map_utils
   2
   3 ! Module that defines constants, data structures, and
   4 ! subroutines used to convert grid indices to lat/lon
   5 ! and vice versa.
   6 !
   7 ! SUPPORTED PROJECTIONS
   8 ! ---------------------
   9 ! Cylindrical Lat/Lon (code = PROJ_LATLON)
  10 ! Mercator (code = PROJ_MERC)
  11 ! Lambert Conformal (code = PROJ_LC)
  12 ! Gaussian (code = PROJ_GAUSS)
  13 ! Polar Stereographic (code = PROJ_PS)
  14 ! Rotated Lat/Lon (code = PROJ_ROTLL)
  15 !
  16 ! REMARKS
  17 ! -------
  18 ! The routines contained within were adapted from routines
  19 ! obtained from NCEP's w3 library.  The original NCEP routines were less
  20 ! flexible (e.g., polar-stereo routines only supported truelat of 60N/60S)
  21 ! than what we needed, so modifications based on equations in Hoke, Hayes, and
  22 ! Renninger (AFGWC/TN/79-003) were added to improve the flexibility.
  23 ! Additionally, coding was improved to F90 standards and the routines were
  24 ! combined into this module.
  25 !
  26 ! ASSUMPTIONS
  27 ! -----------
  28 !  Grid Definition:
  29 !    For mercator, lambert conformal, and polar-stereographic projections,
  30 !    the routines within assume the following:
  31 !
  32 !       1.  Grid is dimensioned (i,j) where i is the East-West direction,
  33 !           positive toward the east, and j is the north-south direction,
  34 !           positive toward the north.
  35 !       2.  Origin is at (1,1) and is located at the southwest corner,
  36 !           regardless of hemispere.
  37 !       3.  Grid spacing (dx) is always positive.
  38 !       4.  Values of true latitudes must be positive for NH domains
  39 !           and negative for SH domains.
  40 !
  41 !     For the latlon and Gaussian projection, the grid origin may be at any
  42 !     of the corners, and the deltalat and deltalon values can be signed to
  43 !     account for this using the following convention:
  44 !       Origin Location        Deltalat Sign      Deltalon Sign
  45 !       ---------------        -------------      -------------
  46 !        SW Corner                  +                   +
  47 !        NE Corner                  -                   -
  48 !        NW Corner                  -                   +
  49 !        SE Corner                  +                   -
  50 !
  51 !  Data Definitions:
  52 !       1. Any arguments that are a latitude value are expressed in
  53 !          degrees north with a valid range of -90 -> 90
  54 !       2. Any arguments that are a longitude value are expressed in
  55 !          degrees east with a valid range of -180 -> 180.
  56 !       3. Distances are in meters and are always positive.
  57 !       4. The standard longitude (stdlon) is defined as the longitude
  58 !          line which is parallel to the grid's y-axis (j-direction), along
  59 !          which latitude increases (NOT the absolute value of latitude, but
  60 !          the actual latitude, such that latitude increases continuously
  61 !          from the south pole to the north pole) as j increases.
  62 !       5. One true latitude value is required for polar-stereographic and
  63 !          mercator projections, and defines at which latitude the
  64 !          grid spacing is true.  For lambert conformal, two true latitude
  65 !          values must be specified, but may be set equal to each other to
  66 !          specify a tangent projection instead of a secant projection.
  67 !
  68 ! USAGE
  69 ! -----
  70 ! To use the routines in this module, the calling routines must have the
  71 ! following statement at the beginning of its declaration block:
  72 !   USE map_utils
  73 !
  74 ! The use of the module not only provides access to the necessary routines,
  75 ! but also defines a structure of TYPE (proj_info) that can be used
  76 ! to declare a variable of the same type to hold your map projection
  77 ! information.  It also defines some integer parameters that contain
  78 ! the projection codes so one only has to use those variable names rather
  79 ! than remembering the acutal code when using them.  The basic steps are
  80 ! as follows:
  81 !
  82 !   1.  Ensure the "USE map_utils" is in your declarations.
  83 !   2.  Declare the projection information structure as type(proj_info):
  84 !         TYPE(proj_info) :: proj
  85 !   3.  Populate your structure by calling the map_set routine:
  86 !         CALL map_set(code,lat1,lon1,knowni,knownj,dx,stdlon,truelat1,truelat2,proj)
  87 !       where:
  88 !         code (input) = one of PROJ_LATLON, PROJ_MERC, PROJ_LC, PROJ_PS,
  89 !                        PROJ_GAUSS, or PROJ_ROTLL
  90 !         lat1 (input) = Latitude of grid origin point (i,j)=(1,1)
  91 !                         (see assumptions!)
  92 !         lon1 (input) = Longitude of grid origin
  93 !         knowni (input) = origin point, x-location
  94 !         knownj (input) = origin point, y-location
  95 !         dx (input) = grid spacing in meters (ignored for LATLON projections)
  96 !         stdlon (input) = Standard longitude for PROJ_PS and PROJ_LC,
  97 !               deltalon (see assumptions) for PROJ_LATLON,
  98 !               ignored for PROJ_MERC
  99 !         truelat1 (input) = 1st true latitude for PROJ_PS, PROJ_LC, and
 100 !                PROJ_MERC, deltalat (see assumptions) for PROJ_LATLON
 101 !         truelat2 (input) = 2nd true latitude for PROJ_LC,
 102 !                ignored for all others.
 103 !         proj (output) = The structure of type (proj_info) that will be fully
 104 !                populated after this call
 105 !
 106 !   4.  Now that the proj structure is populated, you may call either
 107 !       of the following routines:
 108 !
 109 !       latlon_to_ij(proj, lat, lon, i, j)
 110 !       ij_to_latlon(proj, i, j, lat, lon)
 111 !
 112 !       It is incumbent upon the calling routine to determine whether or
 113 !       not the values returned are within your domain's bounds.  All values
 114 !       of i, j, lat, and lon are REAL values.
 115 !
 116 !
 117 ! REFERENCES
 118 ! ----------
 119 !  Hoke, Hayes, and Renninger, "Map Preojections and Grid Systems for
 120 !       Meteorological Applications." AFGWC/TN-79/003(Rev), Air Weather
 121 !       Service, 1985.
 122 !
 123 !  NCAR MM5v3 Modeling System, REGRIDDER program, module_first_guess_map.F
 124 !  NCEP routines w3fb06, w3fb07, w3fb08, w3fb09, w3fb11, w3fb12
 125 !
 126 ! HISTORY
 127 ! -------
 128 ! 27 Mar 2001 - Original Version
 129 !               Brent L. Shaw, NOAA/FSL (CSU/CIRA)
 130 !
 131 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 132
 133    use constants_module
 134    use misc_definitions_module
 135    use module_debug
 136
 137    ! Define some private constants
 138    INTEGER, PRIVATE, PARAMETER :: HIGH = 8
 139
 140    TYPE proj_info
 141
 142       INTEGER          :: code     ! Integer code for projection TYPE
 143       INTEGER          :: nlat     ! For Gaussian -- number of latitude points
 144                                    !  north of the equator
 145       INTEGER          :: nlon     !
 146                                    !
 147       INTEGER          :: ixdim    ! For Rotated Lat/Lon -- number of mass points
 148                                    !  in an odd row
 149       INTEGER          :: jydim    ! For Rotated Lat/Lon -- number of rows
 150       INTEGER          :: stagger  ! For Rotated Lat/Lon -- mass or velocity grid
 151       REAL             :: phi      ! For Rotated Lat/Lon -- domain half-extent in
 152                                    !  degrees latitude
 153       REAL             :: lambda   ! For Rotated Lat/Lon -- domain half-extend in
 154                                    !  degrees longitude
 155       REAL             :: lat1     ! SW latitude (1,1) in degrees (-90->90N)
 156       REAL             :: lon1     ! SW longitude (1,1) in degrees (-180->180E)
 157       REAL             :: lat0     ! For Cassini, latitude of projection pole
 158       REAL             :: lon0     ! For Cassini, longitude of projection pole
 159       REAL             :: dx       ! Grid spacing in meters at truelats, used
 160                                    !  only for ps, lc, and merc projections
 161       REAL             :: dy       ! Grid spacing in meters at truelats, used
 162                                    !  only for ps, lc, and merc projections
 163       REAL             :: latinc   ! Latitude increment for cylindrical lat/lon
 164       REAL             :: loninc   ! Longitude increment for cylindrical lat/lon
 165                                    !  also the lon increment for Gaussian grid
 166       REAL             :: dlat     ! Lat increment for lat/lon grids
 167       REAL             :: dlon     ! Lon increment for lat/lon grids
 168       REAL             :: stdlon   ! Longitude parallel to y-axis (-180->180E)
 169       REAL             :: truelat1 ! First true latitude (all projections)
 170       REAL             :: truelat2 ! Second true lat (LC only)
 171       REAL             :: hemi     ! 1 for NH, -1 for SH
 172       REAL             :: cone     ! Cone factor for LC projections
 173       REAL             :: polei    ! Computed i-location of pole point
 174       REAL             :: polej    ! Computed j-location of pole point
 175       REAL             :: rsw      ! Computed radius to SW corner
 176       REAL             :: rebydx   ! Earth radius divided by dx
 177       REAL             :: knowni   ! X-location of known lat/lon
 178       REAL             :: knownj   ! Y-location of known lat/lon
 179       REAL             :: re_m     ! Radius of spherical earth, meters
 180       REAL             :: rho0     ! For Albers equal area
 181       REAL             :: nc       ! For Albers equal area
 182       REAL             :: bigc     ! For Albers equal area
 183       LOGICAL          :: init     ! Flag to indicate if this struct is
 184                                    !  ready for use
 185       LOGICAL          :: wrap     ! For Gaussian -- flag to indicate wrapping
 186                                    !  around globe?
 187       REAL, POINTER, DIMENSION(:) :: gauss_lat  ! Latitude array for Gaussian grid
 188
 189    END TYPE proj_info
 190
 191  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 192  CONTAINS
 193  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 194
 195    SUBROUTINE map_init(proj)
 196       ! Initializes the map projection structure to missing values
 197
 198       IMPLICIT NONE
 199       TYPE(proj_info), INTENT(INOUT)  :: proj
 200
 201       proj%lat1     = -999.9
 202       proj%lon1     = -999.9
 203       proj%lat0     = -999.9
 204       proj%lon0     = -999.9
 205       proj%dx       = -999.9
 206       proj%dy       = -999.9
 207       proj%latinc   = -999.9
 208       proj%loninc   = -999.9
 209       proj%stdlon   = -999.9
 210       proj%truelat1 = -999.9
 211       proj%truelat2 = -999.9
 212       proj%phi      = -999.9
 213       proj%lambda   = -999.9
 214       proj%ixdim    = -999
 215       proj%jydim    = -999
 216       proj%stagger  = HH
 217       proj%nlat     = 0
 218       proj%nlon     = 0
 219       proj%hemi     = 0.0
 220       proj%cone     = -999.9
 221       proj%polei    = -999.9
 222       proj%polej    = -999.9
 223       proj%rsw      = -999.9
 224       proj%knowni   = -999.9
 225       proj%knownj   = -999.9
 226       proj%re_m     = EARTH_RADIUS_M
 227       proj%init     = .FALSE.
 228       proj%wrap     = .FALSE.
 229       proj%rho0     = 0.
 230       proj%nc       = 0.
 231       proj%bigc     = 0.
 232       nullify(proj%gauss_lat)
 233
 234    END SUBROUTINE map_init
 235
 236
 237    SUBROUTINE map_set(proj_code, proj, lat1, lon1, lat0, lon0, knowni, knownj, dx, latinc, &
 238                       loninc, stdlon, truelat1, truelat2, nlat, nlon, ixdim, jydim, &
 239                       stagger, phi, lambda, r_earth)
 240       ! Given a partially filled proj_info structure, this routine computes
 241       ! polei, polej, rsw, and cone (if LC projection) to complete the
 242       ! structure.  This allows us to eliminate redundant calculations when
 243       ! calling the coordinate conversion routines multiple times for the
 244       ! same map.
 245       ! This will generally be the first routine called when a user wants
 246       ! to be able to use the coordinate conversion routines, and it
 247       ! will call the appropriate subroutines based on the
 248       ! proj%code which indicates which projection type this is.
 249
 250       IMPLICIT NONE
 251
 252       ! Declare arguments
 253       INTEGER, INTENT(IN)               :: proj_code
 254       INTEGER, INTENT(IN), OPTIONAL     :: nlat
 255       INTEGER, INTENT(IN), OPTIONAL     :: nlon
 256       INTEGER, INTENT(IN), OPTIONAL     :: ixdim
 257       INTEGER, INTENT(IN), OPTIONAL     :: jydim
 258       INTEGER, INTENT(IN), OPTIONAL     :: stagger
 259       REAL, INTENT(IN), OPTIONAL        :: latinc
 260       REAL, INTENT(IN), OPTIONAL        :: loninc
 261       REAL, INTENT(IN), OPTIONAL        :: lat1
 262       REAL, INTENT(IN), OPTIONAL        :: lon1
 263       REAL, INTENT(IN), OPTIONAL        :: lat0
 264       REAL, INTENT(IN), OPTIONAL        :: lon0
 265       REAL, INTENT(IN), OPTIONAL        :: dx
 266       REAL, INTENT(IN), OPTIONAL        :: stdlon
 267       REAL, INTENT(IN), OPTIONAL        :: truelat1
 268       REAL, INTENT(IN), OPTIONAL        :: truelat2
 269       REAL, INTENT(IN), OPTIONAL        :: knowni
 270       REAL, INTENT(IN), OPTIONAL        :: knownj
 271       REAL, INTENT(IN), OPTIONAL        :: phi
 272       REAL, INTENT(IN), OPTIONAL        :: lambda
 273       REAL, INTENT(IN), OPTIONAL        :: r_earth
 274       TYPE(proj_info), INTENT(OUT)      :: proj
 275
 276       INTEGER :: iter
 277       REAL :: dummy_lon1
 278       REAL :: dummy_lon0
 279       REAL :: dummy_stdlon
 280
 281       ! First, verify that mandatory parameters are present for the specified proj_code
 282       IF ( proj_code == PROJ_LC ) THEN
 283          IF ( .NOT.PRESENT(truelat1) .OR. &
 284               .NOT.PRESENT(truelat2) .OR. &
 285               .NOT.PRESENT(lat1) .OR. &
 286               .NOT.PRESENT(lon1) .OR. &
 287               .NOT.PRESENT(knowni) .OR. &
 288               .NOT.PRESENT(knownj) .OR. &
 289               .NOT.PRESENT(stdlon) .OR. &
 290               .NOT.PRESENT(dx) ) THEN
 291             PRINT '(A,I2)', 'The following are mandatory parameters for projection code : ', proj_code
 292             PRINT '(A)', ' truelat1, truelat2, lat1, lon1, knowni, knownj, stdlon, dx'
 293             call mprintf(.true.,ERROR,'MAP_INIT')
 294          END IF
 295       ELSE IF ( proj_code == PROJ_PS ) THEN
 296          IF ( .NOT.PRESENT(truelat1) .OR. &
 297               .NOT.PRESENT(lat1) .OR. &
 298               .NOT.PRESENT(lon1) .OR. &
 299               .NOT.PRESENT(knowni) .OR. &
 300               .NOT.PRESENT(knownj) .OR. &
 301               .NOT.PRESENT(stdlon) .OR. &
 302               .NOT.PRESENT(dx) ) THEN
 303             PRINT '(A,I2)', 'The following are mandatory parameters for projection code : ', proj_code
 304             PRINT '(A)', ' truelat1, lat1, lon1, knonwi, knownj, stdlon, dx'
 305             call mprintf(.true.,ERROR,'MAP_INIT')
 306          END IF
 307       ELSE IF ( proj_code == PROJ_PS_WGS84 ) THEN
 308          IF ( .NOT.PRESENT(truelat1) .OR. &
 309               .NOT.PRESENT(lat1) .OR. &
 310               .NOT.PRESENT(lon1) .OR. &
 311               .NOT.PRESENT(knowni) .OR. &
 312               .NOT.PRESENT(knownj) .OR. &
 313               .NOT.PRESENT(stdlon) .OR. &
 314               .NOT.PRESENT(dx) ) THEN
 315             PRINT '(A,I2)', 'The following are mandatory parameters for projection code : ', proj_code
 316             PRINT '(A)', ' truelat1, lat1, lon1, knonwi, knownj, stdlon, dx'
 317             call mprintf(.true.,ERROR,'MAP_INIT')
 318          END IF
 319       ELSE IF ( proj_code == PROJ_ALBERS_NAD83 ) THEN
 320          IF ( .NOT.PRESENT(truelat1) .OR. &
 321               .NOT.PRESENT(truelat2) .OR. &
 322               .NOT.PRESENT(lat1) .OR. &
 323               .NOT.PRESENT(lon1) .OR. &
 324               .NOT.PRESENT(knowni) .OR. &
 325               .NOT.PRESENT(knownj) .OR. &
 326               .NOT.PRESENT(stdlon) .OR. &
 327               .NOT.PRESENT(dx) ) THEN
 328             PRINT '(A,I2)', 'The following are mandatory parameters for projection code : ', proj_code
 329             PRINT '(A)', ' truelat1, truelat2, lat1, lon1, knonwi, knownj, stdlon, dx'
 330             call mprintf(.true.,ERROR,'MAP_INIT')
 331          END IF
 332       ELSE IF ( proj_code == PROJ_MERC ) THEN
 333          IF ( .NOT.PRESENT(truelat1) .OR. &
 334               .NOT.PRESENT(lat1) .OR. &
 335               .NOT.PRESENT(lon1) .OR. &
 336               .NOT.PRESENT(knowni) .OR. &
 337               .NOT.PRESENT(knownj) .OR. &
 338               .NOT.PRESENT(dx) ) THEN
 339             PRINT '(A,I2)', 'The following are mandatory parameters for projection code : ', proj_code
 340             PRINT '(A)', ' truelat1, lat1, lon1, knowni, knownj, dx'
 341             call mprintf(.true.,ERROR,'MAP_INIT')
 342          END IF
 343       ELSE IF ( proj_code == PROJ_LATLON ) THEN
 344          IF ( .NOT.PRESENT(latinc) .OR. &
 345               .NOT.PRESENT(loninc) .OR. &
 346               .NOT.PRESENT(knowni) .OR. &
 347               .NOT.PRESENT(knownj) .OR. &
 348               .NOT.PRESENT(lat1) .OR. &
 349               .NOT.PRESENT(lon1) ) THEN
 350             PRINT '(A,I2)', 'The following are mandatory parameters for projection code : ', proj_code
 351             PRINT '(A)', ' latinc, loninc, knowni, knownj, lat1, lon1'
 352             call mprintf(.true.,ERROR,'MAP_INIT')
 353          END IF
 354       ELSE IF ( proj_code == PROJ_CYL ) THEN
 355          IF ( .NOT.PRESENT(latinc) .OR. &
 356               .NOT.PRESENT(loninc) .OR. &
 357               .NOT.PRESENT(stdlon) ) THEN
 358             PRINT '(A,I2)', 'The following are mandatory parameters for projection code : ', proj_code
 359             PRINT '(A)', ' latinc, loninc, stdlon'
 360             call mprintf(.true.,ERROR,'MAP_INIT')
 361          END IF
 362       ELSE IF ( proj_code == PROJ_CASSINI ) THEN
 363          IF ( .NOT.PRESENT(latinc) .OR. &
 364               .NOT.PRESENT(loninc) .OR. &
 365               .NOT.PRESENT(lat1) .OR. &
 366               .NOT.PRESENT(lon1) .OR. &
 367               .NOT.PRESENT(lat0) .OR. &
 368               .NOT.PRESENT(lon0) .OR. &
 369               .NOT.PRESENT(knowni) .OR. &
 370               .NOT.PRESENT(knownj) .OR. &
 371               .NOT.PRESENT(stdlon) ) THEN
 372             PRINT '(A,I2)', 'The following are mandatory parameters for projection code : ', proj_code
 373             PRINT '(A)', ' latinc, loninc, lat1, lon1, knowni, knownj, lat0, lon0, stdlon'
 374             call mprintf(.true.,ERROR,'MAP_INIT')
 375          END IF
 376       ELSE IF ( proj_code == PROJ_GAUSS ) THEN
 377          IF ( .NOT.PRESENT(nlat) .OR. &
 378               .NOT.PRESENT(lat1) .OR. &
 379               .NOT.PRESENT(lon1) .OR. &
 380               .NOT.PRESENT(loninc) ) THEN
 381             PRINT '(A,I2)', 'The following are mandatory parameters for projection code : ', proj_code
 382             PRINT '(A)', ' nlat, lat1, lon1, loninc'
 383             call mprintf(.true.,ERROR,'MAP_INIT')
 384          END IF
 385       ELSE IF ( proj_code == PROJ_ROTLL ) THEN
 386          IF ( .NOT.PRESENT(ixdim) .OR. &
 387               .NOT.PRESENT(jydim) .OR. &
 388               .NOT.PRESENT(phi) .OR. &
 389               .NOT.PRESENT(lambda) .OR. &
 390               .NOT.PRESENT(lat1) .OR. &
 391               .NOT.PRESENT(lon1) .OR. &
 392               .NOT.PRESENT(stagger) ) THEN
 393             PRINT '(A,I2)', 'The following are mandatory parameters for projection code : ', proj_code
 394             PRINT '(A)', ' ixdim, jydim, phi, lambda, lat1, lon1, stagger'
 395             call mprintf(.true.,ERROR,'MAP_INIT')
 396          END IF
 397       ELSE
 398          PRINT '(A,I2)', 'Unknown projection code: ', proj_code
 399          call mprintf(.true.,ERROR,'MAP_INIT')
 400       END IF
 401
 402       ! Check for validity of mandatory variables in proj
 403       IF ( PRESENT(lat1) ) THEN
 404          IF ( ABS(lat1) .GT. 90. ) THEN
 405             PRINT '(A)', 'Latitude of origin corner required as follows:'
 406             PRINT '(A)', '    -90N <= lat1 < = 90.N'
 407             call mprintf(.true.,ERROR,'MAP_INIT')
 408          ENDIF
 409       ENDIF
 410
 411       IF ( PRESENT(lon1) ) THEN
 412          dummy_lon1 = lon1
 413          IF ( ABS(dummy_lon1) .GT. 180.) THEN
 414             iter = 0
 415             DO WHILE (ABS(dummy_lon1) > 180. .AND. iter < 10)
 416                IF (dummy_lon1 < -180.) dummy_lon1 = dummy_lon1 + 360.
 417                IF (dummy_lon1 > 180.) dummy_lon1 = dummy_lon1 - 360.
 418                iter = iter + 1
 419             END DO
 420             IF (abs(dummy_lon1) > 180.) THEN
 421                PRINT '(A)', 'Longitude of origin required as follows:'
 422                PRINT '(A)', '   -180E <= lon1 <= 180W'
 423                call mprintf(.true.,ERROR,'MAP_INIT')
 424             ENDIF
 425          ENDIF
 426       ENDIF
 427
 428       IF ( PRESENT(lon0) ) THEN
 429          dummy_lon0 = lon0
 430          IF ( ABS(dummy_lon0) .GT. 180.) THEN
 431             iter = 0
 432             DO WHILE (ABS(dummy_lon0) > 180. .AND. iter < 10)
 433                IF (dummy_lon0 < -180.) dummy_lon0 = dummy_lon0 + 360.
 434                IF (dummy_lon0 > 180.) dummy_lon0 = dummy_lon0 - 360.
 435                iter = iter + 1
 436             END DO
 437             IF (abs(dummy_lon0) > 180.) THEN
 438                PRINT '(A)', 'Longitude of pole required as follows:'
 439                PRINT '(A)', '   -180E <= lon0 <= 180W'
 440                call mprintf(.true.,ERROR,'MAP_INIT')
 441             ENDIF
 442          ENDIF
 443       ENDIF
 444
 445       IF ( PRESENT(dx) ) THEN
 446          IF ((dx .LE. 0.).AND.(proj_code .NE. PROJ_LATLON)) THEN
 447             PRINT '(A)', 'Require grid spacing (dx) in meters be positive!'
 448             call mprintf(.true.,ERROR,'MAP_INIT')
 449          ENDIF
 450       ENDIF
 451
 452       IF ( PRESENT(stdlon) ) THEN
 453          dummy_stdlon = stdlon
 454          IF ((ABS(dummy_stdlon) > 180.).AND.(proj_code /= PROJ_MERC)) THEN
 455             iter = 0
 456             DO WHILE (ABS(dummy_stdlon) > 180. .AND. iter < 10)
 457                IF (dummy_stdlon < -180.) dummy_stdlon = dummy_stdlon + 360.
 458                IF (dummy_stdlon > 180.) dummy_stdlon = dummy_stdlon - 360.
 459                iter = iter + 1
 460             END DO
 461             IF (abs(dummy_stdlon) > 180.) THEN
 462                PRINT '(A)', 'Need orientation longitude (stdlon) as: '
 463                PRINT '(A)', '   -180E <= stdlon <= 180W'
 464                call mprintf(.true.,ERROR,'MAP_INIT')
 465             ENDIF
 466          ENDIF
 467       ENDIF
 468
 469       IF ( PRESENT(truelat1) ) THEN
 470          IF (ABS(truelat1).GT.90.) THEN
 471             PRINT '(A)', 'Set true latitude 1 for all projections!'
 472             call mprintf(.true.,ERROR,'MAP_INIT')
 473          ENDIF
 474       ENDIF
 475
 476       CALL map_init(proj)
 477       proj%code  = proj_code
 478       IF ( PRESENT(lat1) )     proj%lat1     = lat1
 479       IF ( PRESENT(lon1) )     proj%lon1     = dummy_lon1
 480       IF ( PRESENT(lat0) )     proj%lat0     = lat0
 481       IF ( PRESENT(lon0) )     proj%lon0     = dummy_lon0
 482       IF ( PRESENT(latinc) )   proj%latinc   = latinc
 483       IF ( PRESENT(loninc) )   proj%loninc   = loninc
 484       IF ( PRESENT(knowni) )   proj%knowni   = knowni
 485       IF ( PRESENT(knownj) )   proj%knownj   = knownj
 486       IF ( PRESENT(dx) )       proj%dx       = dx
 487       IF ( PRESENT(stdlon) )   proj%stdlon   = dummy_stdlon
 488       IF ( PRESENT(truelat1) ) proj%truelat1 = truelat1
 489       IF ( PRESENT(truelat2) ) proj%truelat2 = truelat2
 490       IF ( PRESENT(nlat) )     proj%nlat     = nlat
 491       IF ( PRESENT(nlon) )     proj%nlon     = nlon
 492       IF ( PRESENT(ixdim) )    proj%ixdim    = ixdim
 493       IF ( PRESENT(jydim) )    proj%jydim    = jydim
 494       IF ( PRESENT(stagger) )  proj%stagger  = stagger
 495       IF ( PRESENT(phi) )      proj%phi      = phi
 496       IF ( PRESENT(lambda) )   proj%lambda   = lambda
 497       IF ( PRESENT(r_earth) )  proj%re_m     = r_earth
 498
 499       IF ( PRESENT(dx) ) THEN
 500          IF ( (proj_code == PROJ_LC) .OR. (proj_code == PROJ_PS) .OR. &
 501               (proj_code == PROJ_PS_WGS84) .OR. (proj_code == PROJ_ALBERS_NAD83) .OR. &
 502               (proj_code == PROJ_MERC) ) THEN
 503             proj%dx = dx
 504             IF (truelat1 .LT. 0.) THEN
 505                proj%hemi = -1.0
 506             ELSE
 507                proj%hemi = 1.0
 508             ENDIF
 509             proj%rebydx = proj%re_m / dx
 510          ENDIF
 511       ENDIF
 512
 513       pick_proj: SELECT CASE(proj%code)
 514
 515          CASE(PROJ_PS)
 516             CALL set_ps(proj)
 517
 518          CASE(PROJ_PS_WGS84)
 519             CALL set_ps_wgs84(proj)
 520
 521          CASE(PROJ_ALBERS_NAD83)
 522             CALL set_albers_nad83(proj)
 523
 524          CASE(PROJ_LC)
 525             IF (ABS(proj%truelat2) .GT. 90.) THEN
 526                proj%truelat2=proj%truelat1
 527             ENDIF
 528             CALL set_lc(proj)
 529
 530          CASE (PROJ_MERC)
 531             CALL set_merc(proj)
 532
 533          CASE (PROJ_LATLON)
 534
 535          CASE (PROJ_GAUSS)
 536             CALL set_gauss(proj)
 537
 538          CASE (PROJ_CYL)
 539             CALL set_cyl(proj)
 540
 541          CASE (PROJ_CASSINI)
 542             CALL set_cassini(proj)
 543
 544          CASE (PROJ_ROTLL)
 545
 546       END SELECT pick_proj
 547       proj%init = .TRUE.
 548
 549       RETURN
 550
 551    END SUBROUTINE map_set
 552
 553
 554    SUBROUTINE latlon_to_ij(proj, lat, lon, i, j)
 555       ! Converts input lat/lon values to the cartesian (i,j) value
 556       ! for the given projection.
 557
 558       IMPLICIT NONE
 559       TYPE(proj_info), INTENT(IN)          :: proj
 560       REAL, INTENT(IN)                     :: lat
 561       REAL, INTENT(IN)                     :: lon
 562       REAL, INTENT(OUT)                    :: i
 563       REAL, INTENT(OUT)                    :: j
 564
 565       IF (.NOT.proj%init) THEN
 566          PRINT '(A)', 'You have not called map_set for this projection!'
 567          call mprintf(.true.,ERROR,'LATLON_TO_IJ')
 568       ENDIF
 569
 570       SELECT CASE(proj%code)
 571
 572          CASE(PROJ_LATLON)
 573             CALL llij_latlon(lat,lon,proj,i,j)
 574
 575          CASE(PROJ_MERC)
 576             CALL llij_merc(lat,lon,proj,i,j)
 577
 578          CASE(PROJ_PS)
 579             CALL llij_ps(lat,lon,proj,i,j)
 580
 581          CASE(PROJ_PS_WGS84)
 582             CALL llij_ps_wgs84(lat,lon,proj,i,j)
 583
 584          CASE(PROJ_ALBERS_NAD83)
 585             CALL llij_albers_nad83(lat,lon,proj,i,j)
 586
 587          CASE(PROJ_LC)
 588             CALL llij_lc(lat,lon,proj,i,j)
 589
 590          CASE(PROJ_GAUSS)
 591             CALL llij_gauss(lat,lon,proj,i,j)
 592
 593          CASE(PROJ_CYL)
 594             CALL llij_cyl(lat,lon,proj,i,j)
 595
 596          CASE(PROJ_CASSINI)
 597             CALL llij_cassini(lat,lon,proj,i,j)
 598
 599          CASE(PROJ_ROTLL)
 600             CALL llij_rotlatlon(lat,lon,proj,i,j)
 601
 602          CASE DEFAULT
 603             PRINT '(A,I2)', 'Unrecognized map projection code: ', proj%code
 604             call mprintf(.true.,ERROR,'LATLON_TO_IJ')
 605
 606       END SELECT
 607
 608       RETURN
 609
 610    END SUBROUTINE latlon_to_ij
 611
 612
 613    SUBROUTINE ij_to_latlon(proj, i, j, lat, lon)
 614       ! Computes geographical latitude and longitude for a given (i,j) point
 615       ! in a grid with a projection of proj
 616
 617       IMPLICIT NONE
 618       TYPE(proj_info),INTENT(IN)          :: proj
 619       REAL, INTENT(IN)                    :: i
 620       REAL, INTENT(IN)                    :: j
 621       REAL, INTENT(OUT)                   :: lat
 622       REAL, INTENT(OUT)                   :: lon
 623
 624       IF (.NOT.proj%init) THEN
 625          PRINT '(A)', 'You have not called map_set for this projection!'
 626          call mprintf(.true.,ERROR,'IJ_TO_LATLON')
 627       ENDIF
 628       SELECT CASE (proj%code)
 629
 630          CASE (PROJ_LATLON)
 631             CALL ijll_latlon(i, j, proj, lat, lon)
 632
 633          CASE (PROJ_MERC)
 634             CALL ijll_merc(i, j, proj, lat, lon)
 635
 636          CASE (PROJ_PS)
 637             CALL ijll_ps(i, j, proj, lat, lon)
 638
 639          CASE (PROJ_PS_WGS84)
 640             CALL ijll_ps_wgs84(i, j, proj, lat, lon)
 641
 642          CASE (PROJ_ALBERS_NAD83)
 643             CALL ijll_albers_nad83(i, j, proj, lat, lon)
 644
 645          CASE (PROJ_LC)
 646             CALL ijll_lc(i, j, proj, lat, lon)
 647
 648          CASE (PROJ_CYL)
 649             CALL ijll_cyl(i, j, proj, lat, lon)
 650
 651          CASE (PROJ_CASSINI)
 652             CALL ijll_cassini(i, j, proj, lat, lon)
 653
 654          CASE (PROJ_ROTLL)
 655             CALL ijll_rotlatlon(i, j, proj, lat, lon)
 656
 657          CASE DEFAULT
 658             PRINT '(A,I2)', 'Unrecognized map projection code: ', proj%code
 659             call mprintf(.true.,ERROR,'IJ_TO_LATLON')
 660
 661       END SELECT
 662       RETURN
 663    END SUBROUTINE ij_to_latlon
 664
 665
 666    SUBROUTINE set_ps(proj)
 667       ! Initializes a polar-stereographic map projection from the partially
 668       ! filled proj structure. This routine computes the radius to the
 669       ! southwest corner and computes the i/j location of the pole for use
 670       ! in llij_ps and ijll_ps.
 671       IMPLICIT NONE
 672
 673       ! Declare args
 674       TYPE(proj_info), INTENT(INOUT)    :: proj
 675
 676       ! Local vars
 677       REAL                              :: ala1
 678       REAL                              :: alo1
 679       REAL                              :: reflon
 680       REAL                              :: scale_top
 681
 682       ! Executable code
 683       reflon = proj%stdlon + 90.
 684
 685       ! Compute numerator term of map scale factor
 686       scale_top = 1. + proj%hemi * SIN(proj%truelat1 * rad_per_deg)
 687
 688       ! Compute radius to lower-left (SW) corner
 689       ala1 = proj%lat1 * rad_per_deg
 690       proj%rsw = proj%rebydx*COS(ala1)*scale_top/(1.+proj%hemi*SIN(ala1))
 691
 692       ! Find the pole point
 693       alo1 = (proj%lon1 - reflon) * rad_per_deg
 694       proj%polei = proj%knowni - proj%rsw * COS(alo1)
 695       proj%polej = proj%knownj - proj%hemi * proj%rsw * SIN(alo1)
 696
 697       RETURN
 698
 699    END SUBROUTINE set_ps
 700
 701
 702    SUBROUTINE llij_ps(lat,lon,proj,i,j)
 703       ! Given latitude (-90 to 90), longitude (-180 to 180), and the
 704       ! standard polar-stereographic projection information via the
 705       ! public proj structure, this routine returns the i/j indices which
 706       ! if within the domain range from 1->nx and 1->ny, respectively.
 707
 708       IMPLICIT NONE
 709
 710       ! Delcare input arguments
 711       REAL, INTENT(IN)               :: lat
 712       REAL, INTENT(IN)               :: lon
 713       TYPE(proj_info),INTENT(IN)     :: proj
 714
 715       ! Declare output arguments
 716       REAL, INTENT(OUT)              :: i !(x-index)
 717       REAL, INTENT(OUT)              :: j !(y-index)
 718
 719       ! Declare local variables
 720
 721       REAL                           :: reflon
 722       REAL                           :: scale_top
 723       REAL                           :: ala
 724       REAL                           :: alo
 725       REAL                           :: rm
 726
 727       ! BEGIN CODE
 728
 729       reflon = proj%stdlon + 90.
 730
 731       ! Compute numerator term of map scale factor
 732
 733       scale_top = 1. + proj%hemi * SIN(proj%truelat1 * rad_per_deg)
 734
 735       ! Find radius to desired point
 736       ala = lat * rad_per_deg
 737       rm = proj%rebydx * COS(ala) * scale_top/(1. + proj%hemi *SIN(ala))
 738       alo = (lon - reflon) * rad_per_deg
 739       i = proj%polei + rm * COS(alo)
 740       j = proj%polej + proj%hemi * rm * SIN(alo)
 741
 742       RETURN
 743
 744    END SUBROUTINE llij_ps
 745
 746
 747    SUBROUTINE ijll_ps(i, j, proj, lat, lon)
 748
 749       ! This is the inverse subroutine of llij_ps.  It returns the
 750       ! latitude and longitude of an i/j point given the projection info
 751       ! structure.
 752
 753       IMPLICIT NONE
 754
 755       ! Declare input arguments
 756       REAL, INTENT(IN)                    :: i    ! Column
 757       REAL, INTENT(IN)                    :: j    ! Row
 758       TYPE (proj_info), INTENT(IN)        :: proj
 759
 760       ! Declare output arguments
 761       REAL, INTENT(OUT)                   :: lat     ! -90 -> 90 north
 762       REAL, INTENT(OUT)                   :: lon     ! -180 -> 180 East
 763
 764       ! Local variables
 765       REAL                                :: reflon
 766       REAL                                :: scale_top
 767       REAL                                :: xx,yy
 768       REAL                                :: gi2, r2
 769       REAL                                :: arccos
 770
 771       ! Begin Code
 772
 773       ! Compute the reference longitude by rotating 90 degrees to the east
 774       ! to find the longitude line parallel to the positive x-axis.
 775       reflon = proj%stdlon + 90.
 776
 777       ! Compute numerator term of map scale factor
 778       scale_top = 1. + proj%hemi * SIN(proj%truelat1 * rad_per_deg)
 779
 780       ! Compute radius to point of interest
 781       xx = i - proj%polei
 782       yy = (j - proj%polej) * proj%hemi
 783       r2 = xx**2 + yy**2
 784
 785       ! Now the magic code
 786       IF (r2 .EQ. 0.) THEN
 787          lat = proj%hemi * 90.
 788          lon = reflon
 789       ELSE
 790          gi2 = (proj%rebydx * scale_top)**2.
 791          lat = deg_per_rad * proj%hemi * ASIN((gi2-r2)/(gi2+r2))
 792          arccos = ACOS(xx/SQRT(r2))
 793          IF (yy .GT. 0) THEN
 794             lon = reflon + deg_per_rad * arccos
 795          ELSE
 796             lon = reflon - deg_per_rad * arccos
 797          ENDIF
 798       ENDIF
 799
 800       ! Convert to a -180 -> 180 East convention
 801       IF (lon .GT. 180.) lon = lon - 360.
 802       IF (lon .LT. -180.) lon = lon + 360.
 803
 804       RETURN
 805
 806    END SUBROUTINE ijll_ps
 807
 808
 809    SUBROUTINE set_ps_wgs84(proj)
 810       ! Initializes a polar-stereographic map projection (WGS84 ellipsoid)
 811       ! from the partially filled proj structure. This routine computes the
 812       ! radius to the southwest corner and computes the i/j location of the
 813       ! pole for use in llij_ps and ijll_ps.
 814
 815       IMPLICIT NONE
 816
 817       ! Arguments
 818       TYPE(proj_info), INTENT(INOUT)    :: proj
 819
 820       ! Local variables
 821       real :: h, mc, tc, t, rho
 822
 823       h = proj%hemi
 824
 825       mc = cos(h*proj%truelat1*rad_per_deg)/sqrt(1.0-(E_WGS84*sin(h*proj%truelat1*rad_per_deg))**2.0)
 826       tc = sqrt(((1.0-sin(h*proj%truelat1*rad_per_deg))/(1.0+sin(h*proj%truelat1*rad_per_deg)))* &
 827                 (((1.0+E_WGS84*sin(h*proj%truelat1*rad_per_deg))/(1.0-E_WGS84*sin(h*proj%truelat1*rad_per_deg)))**E_WGS84 ))
 828
 829       ! Find the i/j location of reference lat/lon with respect to the pole of the projection
 830       t = sqrt(((1.0-sin(h*proj%lat1*rad_per_deg))/(1.0+sin(h*proj%lat1*rad_per_deg)))* &
 831                (((1.0+E_WGS84*sin(h*proj%lat1*rad_per_deg))/(1.0-E_WGS84*sin(h*proj%lat1*rad_per_deg)) )**E_WGS84 ) )
 832       rho = h * (A_WGS84 / proj%dx) * mc * t / tc
 833       proj%polei = rho * sin((h*proj%lon1 - h*proj%stdlon)*rad_per_deg)
 834       proj%polej = -rho * cos((h*proj%lon1 - h*proj%stdlon)*rad_per_deg)
 835
 836       RETURN
 837
 838    END SUBROUTINE set_ps_wgs84
 839
 840
 841    SUBROUTINE llij_ps_wgs84(lat,lon,proj,i,j)
 842       ! Given latitude (-90 to 90), longitude (-180 to 180), and the
 843       ! standard polar-stereographic projection information via the
 844       ! public proj structure, this routine returns the i/j indices which
 845       ! if within the domain range from 1->nx and 1->ny, respectively.
 846
 847       IMPLICIT NONE
 848
 849       ! Arguments
 850       REAL, INTENT(IN)               :: lat
 851       REAL, INTENT(IN)               :: lon
 852       REAL, INTENT(OUT)              :: i !(x-index)
 853       REAL, INTENT(OUT)              :: j !(y-index)
 854       TYPE(proj_info),INTENT(IN)     :: proj
 855
 856       ! Local variables
 857       real :: h, mc, tc, t, rho
 858
 859       h = proj%hemi
 860
 861       mc = cos(h*proj%truelat1*rad_per_deg)/sqrt(1.0-(E_WGS84*sin(h*proj%truelat1*rad_per_deg))**2.0)
 862       tc = sqrt(((1.0-sin(h*proj%truelat1*rad_per_deg))/(1.0+sin(h*proj%truelat1*rad_per_deg)))* &
 863                 (((1.0+E_WGS84*sin(h*proj%truelat1*rad_per_deg))/(1.0-E_WGS84*sin(h*proj%truelat1*rad_per_deg)))**E_WGS84 ))
 864
 865       t = sqrt(((1.0-sin(h*lat*rad_per_deg))/(1.0+sin(h*lat*rad_per_deg))) * &
 866                (((1.0+E_WGS84*sin(h*lat*rad_per_deg))/(1.0-E_WGS84*sin(h*lat*rad_per_deg)))**E_WGS84))
 867
 868       ! Find the x/y location of the requested lat/lon with respect to the pole of the projection
 869       rho = (A_WGS84 / proj%dx) * mc * t / tc
 870       i = h *  rho * sin((h*lon - h*proj%stdlon)*rad_per_deg)
 871       j = h *(-rho)* cos((h*lon - h*proj%stdlon)*rad_per_deg)
 872
 873       ! Get i/j relative to reference i/j
 874       i = proj%knowni + (i - proj%polei)
 875       j = proj%knownj + (j - proj%polej)
 876
 877       RETURN
 878
 879    END SUBROUTINE llij_ps_wgs84
 880
 881
 882    SUBROUTINE ijll_ps_wgs84(i, j, proj, lat, lon)
 883
 884       ! This is the inverse subroutine of llij_ps.  It returns the
 885       ! latitude and longitude of an i/j point given the projection info
 886       ! structure.
 887
 888       implicit none
 889
 890       ! Arguments
 891       REAL, INTENT(IN)                    :: i    ! Column
 892       REAL, INTENT(IN)                    :: j    ! Row
 893       REAL, INTENT(OUT)                   :: lat     ! -90 -> 90 north
 894       REAL, INTENT(OUT)                   :: lon     ! -180 -> 180 East
 895       TYPE (proj_info), INTENT(IN)        :: proj
 896
 897       ! Local variables
 898       real :: h, mc, tc, t, rho, x, y
 899       real :: chi, a, b, c, d
 900
 901       h = proj%hemi
 902       x = (i - proj%knowni + proj%polei)
 903       y = (j - proj%knownj + proj%polej)
 904
 905       mc = cos(h*proj%truelat1*rad_per_deg)/sqrt(1.0-(E_WGS84*sin(h*proj%truelat1*rad_per_deg))**2.0)
 906       tc = sqrt(((1.0-sin(h*proj%truelat1*rad_per_deg))/(1.0+sin(h*proj%truelat1*rad_per_deg))) * &
 907                 (((1.0+E_WGS84*sin(h*proj%truelat1*rad_per_deg))/(1.0-E_WGS84*sin(h*proj%truelat1*rad_per_deg)))**E_WGS84 ))
 908
 909       rho = sqrt((x*proj%dx)**2.0 + (y*proj%dx)**2.0)
 910       t = rho * tc / (A_WGS84 * mc)
 911
 912       lon = h*proj%stdlon + h*atan2(h*x,h*(-y))
 913
 914       chi = PI/2.0-2.0*atan(t)
 915       a = 1./2.*E_WGS84**2. + 5./24.*E_WGS84**4. +  1./40.*E_WGS84**6.  +    73./2016.*E_WGS84**8.
 916       b =                     7./24.*E_WGS84**4. + 29./120.*E_WGS84**6. + 54113./40320.*E_WGS84**8.
 917       c =                                           7./30.*E_WGS84**6.  +    81./280.*E_WGS84**8.
 918       d =                                                                  4279./20160.*E_WGS84**8.
 919
 920       lat = chi + sin(2.*chi)*(a + cos(2.*chi)*(b + cos(2.*chi)*(c + d*cos(2.*chi))))
 921       lat = h * lat
 922
 923       lat = lat*deg_per_rad
 924       lon = lon*deg_per_rad
 925
 926       RETURN
 927
 928    END SUBROUTINE ijll_ps_wgs84
 929
 930
 931    SUBROUTINE set_albers_nad83(proj)
 932       ! Initializes an Albers equal area map projection (NAD83 ellipsoid)
 933       ! from the partially filled proj structure. This routine computes the
 934       ! radius to the southwest corner and computes the i/j location of the
 935       ! pole for use in llij_albers_nad83 and ijll_albers_nad83.
 936
 937       IMPLICIT NONE
 938
 939       ! Arguments
 940       TYPE(proj_info), INTENT(INOUT)    :: proj
 941
 942       ! Local variables
 943       real :: h, m1, m2, q1, q2, theta, q, sinphi
 944
 945       h = proj%hemi
 946
 947       m1 = cos(h*proj%truelat1*rad_per_deg)/sqrt(1.0-(E_NAD83*sin(h*proj%truelat1*rad_per_deg))**2.0)
 948       m2 = cos(h*proj%truelat2*rad_per_deg)/sqrt(1.0-(E_NAD83*sin(h*proj%truelat2*rad_per_deg))**2.0)
 949
 950       sinphi = sin(proj%truelat1*rad_per_deg)
 951       q1 = (1.0-E_NAD83**2.0) * &
 952            ((sinphi/(1.0-(E_NAD83*sinphi)**2.0)) - 1.0/(2.0*E_NAD83) * log((1.0-E_NAD83*sinphi)/(1.0+E_NAD83*sinphi)))
 953
 954       sinphi = sin(proj%truelat2*rad_per_deg)
 955       q2 = (1.0-E_NAD83**2.0) * &
 956            ((sinphi/(1.0-(E_NAD83*sinphi)**2.0)) - 1.0/(2.0*E_NAD83) * log((1.0-E_NAD83*sinphi)/(1.0+E_NAD83*sinphi)))
 957
 958       if (proj%truelat1 == proj%truelat2) then
 959          proj%nc = sin(proj%truelat1*rad_per_deg)
 960       else
 961          proj%nc = (m1**2.0 - m2**2.0) / (q2 - q1)
 962       end if
 963
 964       proj%bigc = m1**2.0 + proj%nc*q1
 965
 966       ! Find the i/j location of reference lat/lon with respect to the pole of the projection
 967       sinphi = sin(proj%lat1*rad_per_deg)
 968       q = (1.0-E_NAD83**2.0) * &
 969            ((sinphi/(1.0-(E_NAD83*sinphi)**2.0)) - 1.0/(2.0*E_NAD83) * log((1.0-E_NAD83*sinphi)/(1.0+E_NAD83*sinphi)))
 970
 971       proj%rho0 = h * (A_NAD83 / proj%dx) * sqrt(proj%bigc - proj%nc * q) / proj%nc
 972       theta = proj%nc*(proj%lon1 - proj%stdlon)*rad_per_deg
 973
 974       proj%polei = proj%rho0 * sin(h*theta)
 975       proj%polej = proj%rho0 - proj%rho0 * cos(h*theta)
 976
 977       RETURN
 978
 979    END SUBROUTINE set_albers_nad83
 980
 981
 982    SUBROUTINE llij_albers_nad83(lat,lon,proj,i,j)
 983       ! Given latitude (-90 to 90), longitude (-180 to 180), and the
 984       ! standard projection information via the
 985       ! public proj structure, this routine returns the i/j indices which
 986       ! if within the domain range from 1->nx and 1->ny, respectively.
 987
 988       IMPLICIT NONE
 989
 990       ! Arguments
 991       REAL, INTENT(IN)               :: lat
 992       REAL, INTENT(IN)               :: lon
 993       REAL, INTENT(OUT)              :: i !(x-index)
 994       REAL, INTENT(OUT)              :: j !(y-index)
 995       TYPE(proj_info),INTENT(IN)     :: proj
 996
 997       ! Local variables
 998       real :: h, q, rho, theta, sinphi
 999
1000       h = proj%hemi
1001
1002       sinphi = sin(h*lat*rad_per_deg)
1003
1004       ! Find the x/y location of the requested lat/lon with respect to the pole of the projection
1005       q = (1.0-E_NAD83**2.0) * &
1006            ((sinphi/(1.0-(E_NAD83*sinphi)**2.0)) - 1.0/(2.0*E_NAD83) * log((1.0-E_NAD83*sinphi)/(1.0+E_NAD83*sinphi)))
1007
1008       rho = h * (A_NAD83 / proj%dx) * sqrt(proj%bigc - proj%nc * q) / proj%nc
1009       theta = proj%nc * (h*lon - h*proj%stdlon)*rad_per_deg
1010
1011       i = h*rho*sin(theta)
1012       j = h*proj%rho0 - h*rho*cos(theta)
1013
1014       ! Get i/j relative to reference i/j
1015       i = proj%knowni + (i - proj%polei)
1016       j = proj%knownj + (j - proj%polej)
1017
1018       RETURN
1019
1020    END SUBROUTINE llij_albers_nad83
1021
1022
1023    SUBROUTINE ijll_albers_nad83(i, j, proj, lat, lon)
1024
1025       ! This is the inverse subroutine of llij_albers_nad83.  It returns the
1026       ! latitude and longitude of an i/j point given the projection info
1027       ! structure.
1028
1029       implicit none
1030
1031       ! Arguments
1032       REAL, INTENT(IN)                    :: i    ! Column
1033       REAL, INTENT(IN)                    :: j    ! Row
1034       REAL, INTENT(OUT)                   :: lat     ! -90 -> 90 north
1035       REAL, INTENT(OUT)                   :: lon     ! -180 -> 180 East
1036       TYPE (proj_info), INTENT(IN)        :: proj
1037
1038       ! Local variables
1039       real :: h, q, rho, theta, beta, x, y
1040       real :: a, b, c
1041
1042       h = proj%hemi
1043
1044       x = (i - proj%knowni + proj%polei)
1045       y = (j - proj%knownj + proj%polej)
1046
1047       rho = sqrt(x**2.0 + (proj%rho0 - y)**2.0)
1048       theta = atan2(x, proj%rho0-y)
1049
1050       q = (proj%bigc - (rho*proj%nc*proj%dx/A_NAD83)**2.0) / proj%nc
1051
1052       beta = asin(q/(1.0 - log((1.0-E_NAD83)/(1.0+E_NAD83))*(1.0-E_NAD83**2.0)/(2.0*E_NAD83)))
1053       a = 1./3.*E_NAD83**2. + 31./180.*E_NAD83**4. + 517./5040.*E_NAD83**6.
1054       b =                     23./360.*E_NAD83**4. + 251./3780.*E_NAD83**6.
1055       c =                                            761./45360.*E_NAD83**6.
1056
1057       lat = beta + a*sin(2.*beta) + b*sin(4.*beta) + c*sin(6.*beta)
1058
1059       lat = h*lat*deg_per_rad
1060       lon = proj%stdlon + theta*deg_per_rad/proj%nc
1061
1062       RETURN
1063
1064    END SUBROUTINE ijll_albers_nad83
1065
1066
1067    SUBROUTINE set_lc(proj)
1068       ! Initialize the remaining items in the proj structure for a
1069       ! lambert conformal grid.
1070
1071       IMPLICIT NONE
1072
1073       TYPE(proj_info), INTENT(INOUT)     :: proj
1074
1075       REAL                               :: arg
1076       REAL                               :: deltalon1
1077       REAL                               :: tl1r
1078       REAL                               :: ctl1r
1079
1080       ! Compute cone factor
1081       CALL lc_cone(proj%truelat1, proj%truelat2, proj%cone)
1082
1083       ! Compute longitude differences and ensure we stay out of the
1084       ! forbidden "cut zone"
1085       deltalon1 = proj%lon1 - proj%stdlon
1086       IF (deltalon1 .GT. +180.) deltalon1 = deltalon1 - 360.
1087       IF (deltalon1 .LT. -180.) deltalon1 = deltalon1 + 360.
1088
1089       ! Convert truelat1 to radian and compute COS for later use
1090       tl1r = proj%truelat1 * rad_per_deg
1091       ctl1r = COS(tl1r)
1092
1093       ! Compute the radius to our known lower-left (SW) corner
1094       proj%rsw = proj%rebydx * ctl1r/proj%cone * &
1095              (TAN((90.*proj%hemi-proj%lat1)*rad_per_deg/2.) / &
1096               TAN((90.*proj%hemi-proj%truelat1)*rad_per_deg/2.))**proj%cone
1097
1098       ! Find pole point
1099       arg = proj%cone*(deltalon1*rad_per_deg)
1100       proj%polei = proj%hemi*proj%knowni - proj%hemi * proj%rsw * SIN(arg)
1101       proj%polej = proj%hemi*proj%knownj + proj%rsw * COS(arg)
1102
1103       RETURN
1104
1105    END SUBROUTINE set_lc
1106
1107
1108    SUBROUTINE lc_cone(truelat1, truelat2, cone)
1109
1110    ! Subroutine to compute the cone factor of a Lambert Conformal projection
1111
1112       IMPLICIT NONE
1113
1114       ! Input Args
1115       REAL, INTENT(IN)             :: truelat1  ! (-90 -> 90 degrees N)
1116       REAL, INTENT(IN)             :: truelat2  !   "   "  "   "     "
1117
1118       ! Output Args
1119       REAL, INTENT(OUT)            :: cone
1120
1121       ! Locals
1122
1123       ! BEGIN CODE
1124
1125       ! First, see if this is a secant or tangent projection.  For tangent
1126       ! projections, truelat1 = truelat2 and the cone is tangent to the
1127       ! Earth's surface at this latitude.  For secant projections, the cone
1128       ! intersects the Earth's surface at each of the distinctly different
1129       ! latitudes
1130       IF (ABS(truelat1-truelat2) .GT. 0.1) THEN
1131          cone = ALOG10(COS(truelat1*rad_per_deg)) - &
1132                 ALOG10(COS(truelat2*rad_per_deg))
1133          cone = cone /(ALOG10(TAN((45.0 - ABS(truelat1)/2.0) * rad_per_deg)) - &
1134                 ALOG10(TAN((45.0 - ABS(truelat2)/2.0) * rad_per_deg)))
1135       ELSE
1136          cone = SIN(ABS(truelat1)*rad_per_deg )
1137       ENDIF
1138
1139       RETURN
1140
1141    END SUBROUTINE lc_cone
1142
1143
1144    SUBROUTINE ijll_lc( i, j, proj, lat, lon)
1145
1146    ! Subroutine to convert from the (i,j) cartesian coordinate to the
1147    ! geographical latitude and longitude for a Lambert Conformal projection.
1148
1149    ! History:
1150    ! 25 Jul 01: Corrected by B. Shaw, NOAA/FSL
1151    !
1152       IMPLICIT NONE
1153
1154       ! Input Args
1155       REAL, INTENT(IN)              :: i        ! Cartesian X coordinate
1156       REAL, INTENT(IN)              :: j        ! Cartesian Y coordinate
1157       TYPE(proj_info),INTENT(IN)    :: proj     ! Projection info structure
1158
1159       ! Output Args
1160       REAL, INTENT(OUT)             :: lat      ! Latitude (-90->90 deg N)
1161       REAL, INTENT(OUT)             :: lon      ! Longitude (-180->180 E)
1162
1163       ! Locals
1164       REAL                          :: inew
1165       REAL                          :: jnew
1166       REAL                          :: r
1167       REAL                          :: chi,chi1,chi2
1168       REAL                          :: r2
1169       REAL                          :: xx
1170       REAL                          :: yy
1171
1172       ! BEGIN CODE
1173
1174       chi1 = (90. - proj%hemi*proj%truelat1)*rad_per_deg
1175       chi2 = (90. - proj%hemi*proj%truelat2)*rad_per_deg
1176
1177       ! See if we are in the southern hemispere and flip the indices
1178       ! if we are.
1179       inew = proj%hemi * i
1180       jnew = proj%hemi * j
1181
1182       ! Compute radius**2 to i/j location
1183       xx = inew - proj%polei
1184       yy = proj%polej - jnew
1185       r2 = (xx*xx + yy*yy)
1186       r = SQRT(r2)/proj%rebydx
1187
1188       ! Convert to lat/lon
1189       IF (r2 .EQ. 0.) THEN
1190          lat = proj%hemi * 90.
1191          lon = proj%stdlon
1192       ELSE
1193
1194          ! Longitude
1195          lon = proj%stdlon + deg_per_rad * ATAN2(proj%hemi*xx,yy)/proj%cone
1196          lon = AMOD(lon+360., 360.)
1197
1198          ! Latitude.  Latitude determined by solving an equation adapted
1199          ! from:
1200          !  Maling, D.H., 1973: Coordinate Systems and Map Projections
1201          ! Equations #20 in Appendix I.
1202
1203          IF (chi1 .EQ. chi2) THEN
1204             chi = 2.0*ATAN( ( r/TAN(chi1) )**(1./proj%cone) * TAN(chi1*0.5) )
1205          ELSE
1206             chi = 2.0*ATAN( (r*proj%cone/SIN(chi1))**(1./proj%cone) * TAN(chi1*0.5))
1207          ENDIF
1208          lat = (90.0-chi*deg_per_rad)*proj%hemi
1209
1210       ENDIF
1211
1212       IF (lon .GT. +180.) lon = lon - 360.
1213       IF (lon .LT. -180.) lon = lon + 360.
1214
1215       RETURN
1216
1217    END SUBROUTINE ijll_lc
1218
1219
1220    SUBROUTINE llij_lc( lat, lon, proj, i, j)
1221
1222    ! Subroutine to compute the geographical latitude and longitude values
1223    ! to the cartesian x/y on a Lambert Conformal projection.
1224
1225       IMPLICIT NONE
1226
1227       ! Input Args
1228       REAL, INTENT(IN)              :: lat      ! Latitude (-90->90 deg N)
1229       REAL, INTENT(IN)              :: lon      ! Longitude (-180->180 E)
1230       TYPE(proj_info),INTENT(IN)      :: proj     ! Projection info structure
1231
1232       ! Output Args
1233       REAL, INTENT(OUT)             :: i        ! Cartesian X coordinate
1234       REAL, INTENT(OUT)             :: j        ! Cartesian Y coordinate
1235
1236       ! Locals
1237       REAL                          :: arg
1238       REAL                          :: deltalon
1239       REAL                          :: tl1r
1240       REAL                          :: rm
1241       REAL                          :: ctl1r
1242
1243
1244       ! BEGIN CODE
1245
1246       ! Compute deltalon between known longitude and standard lon and ensure
1247       ! it is not in the cut zone
1248       deltalon = lon - proj%stdlon
1249       IF (deltalon .GT. +180.) deltalon = deltalon - 360.
1250       IF (deltalon .LT. -180.) deltalon = deltalon + 360.
1251
1252       ! Convert truelat1 to radian and compute COS for later use
1253       tl1r = proj%truelat1 * rad_per_deg
1254       ctl1r = COS(tl1r)
1255
1256       ! Radius to desired point
1257       rm = proj%rebydx * ctl1r/proj%cone * &
1258            (TAN((90.*proj%hemi-lat)*rad_per_deg/2.) / &
1259             TAN((90.*proj%hemi-proj%truelat1)*rad_per_deg/2.))**proj%cone
1260
1261       arg = proj%cone*(deltalon*rad_per_deg)
1262       i = proj%polei + proj%hemi * rm * SIN(arg)
1263       j = proj%polej - rm * COS(arg)
1264
1265       ! Finally, if we are in the southern hemisphere, flip the i/j
1266       ! values to a coordinate system where (1,1) is the SW corner
1267       ! (what we assume) which is different than the original NCEP
1268       ! algorithms which used the NE corner as the origin in the
1269       ! southern hemisphere (left-hand vs. right-hand coordinate?)
1270       i = proj%hemi * i
1271       j = proj%hemi * j
1272
1273       RETURN
1274    END SUBROUTINE llij_lc
1275
1276
1277    SUBROUTINE set_merc(proj)
1278
1279       ! Sets up the remaining basic elements for the mercator projection
1280
1281       IMPLICIT NONE
1282       TYPE(proj_info), INTENT(INOUT)       :: proj
1283       REAL                                 :: clain
1284
1285
1286       !  Preliminary variables
1287
1288       clain = COS(rad_per_deg*proj%truelat1)
1289       proj%dlon = proj%dx / (proj%re_m * clain)
1290
1291       ! Compute distance from equator to origin, and store in the
1292       ! proj%rsw tag.
1293
1294       proj%rsw = 0.
1295       IF (proj%lat1 .NE. 0.) THEN
1296          proj%rsw = (ALOG(TAN(0.5*((proj%lat1+90.)*rad_per_deg))))/proj%dlon
1297       ENDIF
1298
1299       RETURN
1300
1301    END SUBROUTINE set_merc
1302
1303
1304    SUBROUTINE llij_merc(lat, lon, proj, i, j)
1305
1306       ! Compute i/j coordinate from lat lon for mercator projection
1307
1308       IMPLICIT NONE
1309       REAL, INTENT(IN)              :: lat
1310       REAL, INTENT(IN)              :: lon
1311       TYPE(proj_info),INTENT(IN)    :: proj
1312       REAL,INTENT(OUT)              :: i
1313       REAL,INTENT(OUT)              :: j
1314       REAL                          :: deltalon
1315
1316       deltalon = lon - proj%lon1
1317       IF (deltalon .LT. -180.) deltalon = deltalon + 360.
1318       IF (deltalon .GT. 180.) deltalon = deltalon - 360.
1319       i = proj%knowni + (deltalon/(proj%dlon*deg_per_rad))
1320       j = proj%knownj + (ALOG(TAN(0.5*((lat + 90.) * rad_per_deg)))) / &
1321              proj%dlon - proj%rsw
1322
1323       RETURN
1324
1325    END SUBROUTINE llij_merc
1326
1327
1328    SUBROUTINE ijll_merc(i, j, proj, lat, lon)
1329
1330       ! Compute the lat/lon from i/j for mercator projection
1331
1332       IMPLICIT NONE
1333       REAL,INTENT(IN)               :: i
1334       REAL,INTENT(IN)               :: j
1335       TYPE(proj_info),INTENT(IN)    :: proj
1336       REAL, INTENT(OUT)             :: lat
1337       REAL, INTENT(OUT)             :: lon
1338
1339
1340       lat = 2.0*ATAN(EXP(proj%dlon*(proj%rsw + j-proj%knownj)))*deg_per_rad - 90.
1341       lon = (i-proj%knowni)*proj%dlon*deg_per_rad + proj%lon1
1342       IF (lon.GT.180.) lon = lon - 360.
1343       IF (lon.LT.-180.) lon = lon + 360.
1344       RETURN
1345
1346    END SUBROUTINE ijll_merc
1347
1348
1349    SUBROUTINE llij_latlon(lat, lon, proj, i, j)
1350
1351       ! Compute the i/j location of a lat/lon on a LATLON grid.
1352       IMPLICIT NONE
1353       REAL, INTENT(IN)             :: lat
1354       REAL, INTENT(IN)             :: lon
1355       TYPE(proj_info), INTENT(IN)  :: proj
1356       REAL, INTENT(OUT)            :: i
1357       REAL, INTENT(OUT)            :: j
1358
1359       REAL                         :: deltalat
1360       REAL                         :: deltalon
1361
1362       ! Compute deltalat and deltalon as the difference between the input
1363       ! lat/lon and the origin lat/lon
1364       deltalat = lat - proj%lat1
1365       deltalon = lon - proj%lon1
1366
1367       ! Compute i/j
1368       i = deltalon/proj%loninc
1369       j = deltalat/proj%latinc
1370
1371       i = i + proj%knowni
1372       j = j + proj%knownj
1373
1374       RETURN
1375
1376    END SUBROUTINE llij_latlon
1377
1378
1379    SUBROUTINE ijll_latlon(i, j, proj, lat, lon)
1380
1381       ! Compute the lat/lon location of an i/j on a LATLON grid.
1382       IMPLICIT NONE
1383       REAL, INTENT(IN)             :: i
1384       REAL, INTENT(IN)             :: j
1385       TYPE(proj_info), INTENT(IN)  :: proj
1386       REAL, INTENT(OUT)            :: lat
1387       REAL, INTENT(OUT)            :: lon
1388
1389       REAL                         :: i_work, j_work
1390       REAL                         :: deltalat
1391       REAL                         :: deltalon
1392
1393       i_work = i - proj%knowni
1394       j_work = j - proj%knownj
1395
1396       ! Compute deltalat and deltalon
1397       deltalat = j_work*proj%latinc
1398       deltalon = i_work*proj%loninc
1399
1400       lat = proj%lat1 + deltalat
1401       lon = proj%lon1 + deltalon
1402
1403       RETURN
1404
1405    END SUBROUTINE ijll_latlon
1406
1407
1408    SUBROUTINE set_cyl(proj)
1409
1410       implicit none
1411
1412       ! Arguments
1413       type(proj_info), intent(inout) :: proj
1414
1415       proj%hemi = 1.0
1416
1417    END SUBROUTINE set_cyl
1418
1419
1420    SUBROUTINE llij_cyl(lat, lon, proj, i, j)
1421
1422       implicit none
1423
1424       ! Arguments
1425       real, intent(in) :: lat, lon
1426       real, intent(out) :: i, j
1427       type(proj_info), intent(in) :: proj
1428
1429       ! Local variables
1430       real :: deltalat
1431       real :: deltalon
1432
1433       ! Compute deltalat and deltalon as the difference between the input
1434       ! lat/lon and the origin lat/lon
1435       deltalat = lat - proj%lat1
1436 !      deltalon = lon - proj%stdlon
1437       deltalon = lon - proj%lon1
1438
1439       if (deltalon <   0.) deltalon = deltalon + 360.
1440       if (deltalon > 360.) deltalon = deltalon - 360.
1441
1442       ! Compute i/j
1443       i = deltalon/proj%loninc
1444       j = deltalat/proj%latinc
1445
1446       if (i <= 0.)              i = i + 360./proj%loninc
1447       if (i > 360./proj%loninc) i = i - 360./proj%loninc
1448
1449       i = i + proj%knowni
1450       j = j + proj%knowni
1451
1452    END SUBROUTINE llij_cyl
1453
1454
1455    SUBROUTINE ijll_cyl(i, j, proj, lat, lon)
1456
1457       implicit none
1458
1459       ! Arguments
1460       real, intent(in) :: i, j
1461       real, intent(out) :: lat, lon
1462       type(proj_info), intent(in) :: proj
1463
1464       ! Local variables
1465       real :: deltalat
1466       real :: deltalon
1467       real :: i_work, j_work
1468
1469       i_work = i - proj%knowni
1470       j_work = j - proj%knownj
1471
1472       if (i_work < 0.)              i_work = i_work + 360./proj%loninc
1473       if (i_work >= 360./proj%loninc) i_work = i_work - 360./proj%loninc
1474
1475       ! Compute deltalat and deltalon
1476       deltalat = j_work*proj%latinc
1477       deltalon = i_work*proj%loninc
1478
1479       lat = deltalat + proj%lat1
1480 !      lon = deltalon + proj%stdlon
1481       lon = deltalon + proj%lon1
1482
1483       if (lon < -180.) lon = lon + 360.
1484       if (lon >  180.) lon = lon - 360.
1485
1486    END SUBROUTINE ijll_cyl
1487
1488
1489    SUBROUTINE set_cassini(proj)
1490
1491       implicit none
1492
1493       ! Arguments
1494       type(proj_info), intent(inout) :: proj
1495
1496       ! Local variables
1497       real :: comp_lat, comp_lon
1498       logical :: global_domain
1499
1500       proj%hemi = 1.0
1501
1502       ! Try to determine whether this domain has global coverage
1503       if (abs(proj%lat1 - proj%latinc/2. + 90.) < 0.001 .and. &
1504           abs(mod(proj%lon1 - proj%loninc/2. - proj%stdlon,360.)) < 0.001) then
1505          global_domain = .true.
1506       else
1507          global_domain = .false.
1508       end if
1509
1510       if (abs(proj%lat0) /= 90. .and. .not.global_domain) then
1511          call rotate_coords(proj%lat1,proj%lon1,comp_lat,comp_lon,proj%lat0,proj%lon0,proj%stdlon,-1)
1512          proj%lat1 = comp_lat
1513          proj%lon1 = comp_lon
1514       end if
1515
1516    END SUBROUTINE set_cassini
1517
1518
1519    SUBROUTINE llij_cassini(lat, lon, proj, i, j)
1520
1521       implicit none
1522
1523       ! Arguments
1524       real, intent(in) :: lat, lon
1525       real, intent(out) :: i, j
1526       type(proj_info), intent(in) :: proj
1527
1528       ! Local variables
1529       real :: comp_lat, comp_lon
1530
1531       ! Convert geographic to computational lat/lon
1532       if (abs(proj%lat0) /= 90.) then
1533          call rotate_coords(lat,lon,comp_lat,comp_lon,proj%lat0,proj%lon0,proj%stdlon,-1)
1534       else
1535          comp_lat = lat
1536          comp_lon = lon
1537       end if
1538
1539       ! Convert computational lat/lon to i/j
1540       call llij_cyl(comp_lat, comp_lon, proj, i, j)
1541
1542    END SUBROUTINE llij_cassini
1543
1544
1545    SUBROUTINE ijll_cassini(i, j, proj, lat, lon)
1546
1547       implicit none
1548
1549       ! Arguments
1550       real, intent(in) :: i, j
1551       real, intent(out) :: lat, lon
1552       type(proj_info), intent(in) :: proj
1553
1554       ! Local variables
1555       real :: comp_lat, comp_lon
1556
1557       ! Convert i/j to computational lat/lon
1558       call ijll_cyl(i, j, proj, comp_lat, comp_lon)
1559
1560       ! Convert computational to geographic lat/lon
1561       if (abs(proj%lat0) /= 90.) then
1562          call rotate_coords(comp_lat,comp_lon,lat,lon,proj%lat0,proj%lon0,proj%stdlon,1)
1563       else
1564          lat = comp_lat
1565          lon = comp_lon
1566       end if
1567
1568    END SUBROUTINE ijll_cassini
1569
1570
1571    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1572    ! Purpose: Converts between computational and geographic lat/lon for Cassini
1573    !
1574    ! Notes: This routine was provided by Bill Skamarock, 2007-03-27
1575    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1576    SUBROUTINE rotate_coords(ilat,ilon,olat,olon,lat_np,lon_np,lon_0,direction)
1577
1578       IMPLICIT NONE
1579
1580       REAL, INTENT(IN   ) :: ilat, ilon
1581       REAL, INTENT(  OUT) :: olat, olon
1582       REAL, INTENT(IN   ) :: lat_np, lon_np, lon_0
1583       INTEGER, INTENT(IN  ), OPTIONAL :: direction
1584       ! >=0, default : computational -> geographical
1585       ! < 0          : geographical  -> computational
1586
1587       REAL :: rlat, rlon
1588       REAL :: phi_np, lam_np, lam_0, dlam
1589       REAL :: sinphi, cosphi, coslam, sinlam
1590
1591       ! Convert all angles to radians
1592       phi_np = lat_np * rad_per_deg
1593       lam_np = lon_np * rad_per_deg
1594       lam_0  = lon_0  * rad_per_deg
1595       rlat = ilat * rad_per_deg
1596       rlon = ilon * rad_per_deg
1597
1598       IF (PRESENT(direction) .AND. (direction < 0)) THEN
1599          ! The equations are exactly the same except for one small difference
1600          ! with respect to longitude ...
1601          dlam = PI - lam_0
1602       ELSE
1603          dlam = lam_np
1604       END IF
1605       sinphi = COS(phi_np)*COS(rlat)*COS(rlon-dlam) + SIN(phi_np)*SIN(rlat)
1606       cosphi = SQRT(1.-sinphi*sinphi)
1607       coslam = SIN(phi_np)*COS(rlat)*COS(rlon-dlam) - COS(phi_np)*SIN(rlat)
1608       sinlam = COS(rlat)*SIN(rlon-dlam)
1609       IF ( cosphi /= 0. ) THEN
1610          coslam = coslam/cosphi
1611          sinlam = sinlam/cosphi
1612       END IF
1613       olat = deg_per_rad*ASIN(sinphi)
1614       olon = deg_per_rad*(ATAN2(sinlam,coslam)-dlam-lam_0+lam_np)
1615       ! Both of my F90 text books prefer the DO-EXIT form, and claim it is faster
1616       ! when optimization is turned on (as we will always do...)
1617       DO
1618          IF (olon >= -180.) EXIT
1619          olon = olon + 360.
1620       END DO
1621       DO
1622          IF (olon <=  180.) EXIT
1623          olon = olon - 360.
1624       END DO
1625
1626    END SUBROUTINE rotate_coords
1627
1628
1629    SUBROUTINE llij_rotlatlon(lat, lon, proj, i, j)
1630
1631       IMPLICIT NONE
1632
1633       ! Arguments
1634       REAL, INTENT(IN) :: lat, lon
1635       REAL, INTENT(OUT) :: i, j
1636       TYPE (proj_info), INTENT(IN) :: proj
1637
1638       ! Local variables
1639       REAL(KIND=HIGH) :: dphd,dlmd !Grid increments, degrees
1640       INTEGER :: ii,jj,jmt,ncol,nrow
1641       REAL(KIND=HIGH) :: glatd  !Geographic latitude, positive north
1642       REAL(KIND=HIGH) :: glond  !Geographic longitude, positive west
1643       REAL(KIND=HIGH) :: col,d1,d2,d2r,dlm,dlm1,dlm2,dph,glat,glon,    &
1644               pi,r2d,row,tlat,tlat1,tlat2,              &
1645               tlon,tlon1,tlon2,tph0,tlm0,x,y,z
1646
1647       glatd = lat
1648       glond = -lon
1649
1650       dphd = proj%phi/REAL((proj%jydim-1)/2)
1651       dlmd = proj%lambda/REAL(proj%ixdim-1)
1652
1653       pi = ACOS(-1.0)
1654       d2r = pi/180.
1655       r2d = 1./d2r
1656
1657       jmt = proj%jydim/2+1
1658
1659       glat = glatd*d2r
1660       glon = glond*d2r
1661       dph = dphd*d2r
1662       dlm = dlmd*d2r
1663       tph0 = proj%lat1*d2r
1664       tlm0 = -proj%lon1*d2r
1665
1666       x = COS(tph0)*COS(glat)*COS(glon-tlm0)+SIN(tph0)*SIN(glat)
1667       y = -COS(glat)*SIN(glon-tlm0)
1668       z = COS(tph0)*SIN(glat)-SIN(tph0)*COS(glat)*COS(glon-tlm0)
1669       tlat = r2d*ATAN(z/SQRT(x*x+y*y))
1670       tlon = r2d*ATAN(y/x)
1671
1672       row = tlat/dphd+jmt
1673       col = tlon/dlmd+proj%ixdim
1674
1675       nrow = NINT(row)
1676       ncol = NINT(col)
1677
1678       tlat = tlat*d2r
1679       tlon = tlon*d2r
1680
1681       IF (proj%stagger == HH) THEN
1682
1683          IF ((abs(MOD(nrow,2)) == 1 .AND. abs(MOD(ncol,2)) == 1) .OR. &
1684              (MOD(nrow,2) == 0 .AND. MOD(ncol,2) == 0)) THEN
1685
1686             tlat1 = (nrow-jmt)*dph
1687             tlat2 = tlat1+dph
1688             tlon1 = (ncol-proj%ixdim)*dlm
1689             tlon2 = tlon1+dlm
1690
1691             dlm1 = tlon-tlon1
1692             dlm2 = tlon-tlon2
1693             d1 = ACOS(COS(tlat)*COS(tlat1)*COS(dlm1)+SIN(tlat)*SIN(tlat1))
1694             d2 = ACOS(COS(tlat)*COS(tlat2)*COS(dlm2)+SIN(tlat)*SIN(tlat2))
1695
1696             IF (d1 > d2) THEN
1697                nrow = nrow+1
1698                ncol = ncol+1
1699             END IF
1700
1701          ELSE
1702             tlat1 = (nrow+1-jmt)*dph
1703             tlat2 = tlat1-dph
1704             tlon1 = (ncol-proj%ixdim)*dlm
1705             tlon2 = tlon1+dlm
1706             dlm1 = tlon-tlon1
1707             dlm2 = tlon-tlon2
1708             d1 = ACOS(COS(tlat)*COS(tlat1)*COS(dlm1)+SIN(tlat)*SIN(tlat1))
1709             d2 = ACOS(COS(tlat)*COS(tlat2)*COS(dlm2)+SIN(tlat)*SIN(tlat2))
1710
1711             IF (d1 < d2) THEN
1712                nrow = nrow+1
1713             ELSE
1714                ncol = ncol+1
1715             END IF
1716          END IF
1717
1718       ELSE IF (proj%stagger == VV) THEN
1719
1720          IF ((MOD(nrow,2) == 0 .AND. abs(MOD(ncol,2)) == 1) .OR. &
1721              (abs(MOD(nrow,2)) == 1 .AND. MOD(ncol,2) == 0)) THEN
1722             tlat1 = (nrow-jmt)*dph
1723             tlat2 = tlat1+dph
1724             tlon1 = (ncol-proj%ixdim)*dlm
1725             tlon2 = tlon1+dlm
1726             dlm1 = tlon-tlon1
1727             dlm2 = tlon-tlon2
1728             d1 = ACOS(COS(tlat)*COS(tlat1)*COS(dlm1)+SIN(tlat)*SIN(tlat1))
1729             d2 = ACOS(COS(tlat)*COS(tlat2)*COS(dlm2)+SIN(tlat)*SIN(tlat2))
1730
1731             IF (d1 > d2) THEN
1732                nrow = nrow+1
1733                ncol = ncol+1
1734             END IF
1735
1736          ELSE
1737             tlat1 = (nrow+1-jmt)*dph
1738             tlat2 = tlat1-dph
1739             tlon1 = (ncol-proj%ixdim)*dlm
1740             tlon2 = tlon1+dlm
1741             dlm1 = tlon-tlon1
1742             dlm2 = tlon-tlon2
1743             d1 = ACOS(COS(tlat)*COS(tlat1)*COS(dlm1)+SIN(tlat)*SIN(tlat1))
1744             d2 = ACOS(COS(tlat)*COS(tlat2)*COS(dlm2)+SIN(tlat)*SIN(tlat2))
1745
1746             IF (d1 < d2) THEN
1747                nrow = nrow+1
1748             ELSE
1749                ncol = ncol+1
1750             END IF
1751          END IF
1752       END IF
1753
1754
1755 !!! Added next line as a Kludge - not yet understood why needed
1756         if (ncol .le. 0) ncol=ncol-1
1757
1758       jj = nrow
1759       ii = ncol/2
1760
1761       IF (proj%stagger == HH) THEN
1762          IF (abs(MOD(jj,2)) == 1) ii = ii+1
1763       ELSE IF (proj%stagger == VV) THEN
1764          IF (MOD(jj,2) == 0) ii=ii+1
1765       END IF
1766
1767       i = REAL(ii)
1768       j = REAL(jj)
1769
1770    END SUBROUTINE llij_rotlatlon
1771
1772
1773    SUBROUTINE ijll_rotlatlon(i, j, proj, lat,lon)
1774
1775       IMPLICIT NONE
1776
1777       ! Arguments
1778       REAL, INTENT(IN) :: i, j
1779       REAL, INTENT(OUT) :: lat, lon
1780       TYPE (proj_info), INTENT(IN) :: proj
1781
1782       ! Local variables
1783       INTEGER :: ih,jh
1784       INTEGER :: midcol,midrow,ncol
1785       REAL :: dphd,dlmd !Grid increments, degrees
1786       REAL(KIND=HIGH) :: arg1,arg2,d2r,fctr,glatr,glatd,glond,pi, &
1787               r2d,tlatd,tlond,tlatr,tlonr,tlm0,tph0
1788
1789       ih = NINT(i)
1790       jh = NINT(j)
1791
1792       dphd = proj%phi/REAL((proj%jydim-1)/2)
1793       dlmd = proj%lambda/REAL(proj%ixdim-1)
1794
1795       pi = ACOS(-1.0)
1796       d2r = pi/180.
1797       r2d = 1./d2r
1798       tph0 = proj%lat1*d2r
1799       tlm0 = -proj%lon1*d2r
1800
1801       midrow = (proj%jydim+1)/2
1802       midcol = proj%ixdim
1803
1804 !      IF (proj%stagger == HH) THEN
1805
1806 !!!         ncol = 2*ih-1+MOD(jh+1,2)
1807          ncol = 2*ih-1+abs(MOD(jh+1,2))
1808          tlatd = (jh-midrow)*dphd
1809          tlond = (ncol-midcol)*dlmd
1810        IF (proj%stagger == VV) THEN
1811           if (mod(jh,2) .eq. 0) then
1812              tlond = tlond - DLMD
1813           else
1814              tlond = tlond + DLMD
1815           end if
1816        END IF
1817
1818       tlatr = tlatd*d2r
1819       tlonr = tlond*d2r
1820       arg1 = SIN(tlatr)*COS(tph0)+COS(tlatr)*SIN(tph0)*COS(tlonr)
1821       glatr = ASIN(arg1)
1822
1823       glatd = glatr*r2d
1824
1825       arg2 = COS(tlatr)*COS(tlonr)/(COS(glatr)*COS(tph0))-TAN(glatr)*TAN(tph0)
1826       IF (ABS(arg2) > 1.) arg2 = ABS(arg2)/arg2
1827       fctr = 1.
1828       IF (tlond > 0.) fctr = -1.
1829
1830       glond = tlm0*r2d+fctr*ACOS(arg2)*r2d
1831
1832       lat = glatd
1833       lon = -glond
1834
1835       IF (lon .GT. +180.) lon = lon - 360.
1836       IF (lon .LT. -180.) lon = lon + 360.
1837
1838    END SUBROUTINE ijll_rotlatlon
1839
1840
1841    SUBROUTINE set_gauss(proj)
1842
1843       IMPLICIT NONE
1844
1845       ! Argument
1846       type (proj_info), intent(inout) :: proj
1847
1848       !  Initialize the array that will hold the Gaussian latitudes.
1849
1850       IF ( ASSOCIATED( proj%gauss_lat ) ) THEN
1851          DEALLOCATE ( proj%gauss_lat )
1852       END IF
1853
1854       !  Get the needed space for our array.
1855
1856       ALLOCATE ( proj%gauss_lat(proj%nlat*2) )
1857
1858       !  Compute the Gaussian latitudes.
1859
1860       CALL gausll( proj%nlat*2 , proj%gauss_lat )
1861
1862       !  Now, these could be upside down from what we want, so let's check.
1863       !  We take advantage of the equatorial symmetry to remove any sort of
1864       !  array re-ordering.
1865
1866       IF ( ABS(proj%gauss_lat(1) - proj%lat1) .GT. 0.01 ) THEN
1867          proj%gauss_lat = -1. * proj%gauss_lat
1868       END IF
1869
1870       !  Just a sanity check.
1871
1872       IF ( ABS(proj%gauss_lat(1) - proj%lat1) .GT. 0.01 ) THEN
1873          PRINT '(A)','Oops, something is not right with the Gaussian latitude computation.'
1874          PRINT '(A,F8.3,A)','The input data gave the starting latitude as ',proj%lat1,'.'
1875          PRINT '(A,F8.3,A)','This routine computed the starting latitude as +-',ABS(proj%gauss_lat(1)),'.'
1876          PRINT '(A,F8.3,A)','The difference is larger than 0.01 degrees, which is not expected.'
1877          call mprintf(.true.,ERROR,'Gaussian_latitude_computation')
1878       END IF
1879
1880    END SUBROUTINE set_gauss
1881
1882
1883    SUBROUTINE gausll ( nlat , lat_sp )
1884
1885       IMPLICIT NONE
1886
1887       INTEGER                            :: nlat , i
1888       REAL (KIND=HIGH) , PARAMETER       :: pi = 3.141592653589793
1889       REAL (KIND=HIGH) , DIMENSION(nlat) :: cosc , gwt , sinc , colat , wos2 , lat
1890       REAL             , DIMENSION(nlat) :: lat_sp
1891
1892       CALL lggaus(nlat, cosc, gwt, sinc, colat, wos2)
1893
1894       DO i = 1, nlat
1895          lat(i) = ACOS(sinc(i)) * 180._HIGH / pi
1896          IF (i.gt.nlat/2) lat(i) = -lat(i)
1897       END DO
1898
1899       lat_sp = REAL(lat)
1900
1901    END SUBROUTINE gausll
1902
1903
1904    SUBROUTINE lggaus( nlat, cosc, gwt, sinc, colat, wos2 )
1905
1906       IMPLICIT NONE
1907
1908       !  LGGAUS finds the Gaussian latitudes by finding the roots of the
1909       !  ordinary Legendre polynomial of degree NLAT using Newton's
1910       !  iteration method.
1911
1912       !  On entry:
1913             integer NLAT ! the number of latitudes (degree of the polynomial)
1914
1915       !  On exit: for each Gaussian latitude
1916       !     COSC   - cos(colatitude) or sin(latitude)
1917       !     GWT    - the Gaussian weights
1918       !     SINC   - sin(colatitude) or cos(latitude)
1919       !     COLAT  - the colatitudes in radians
1920       !     WOS2   - Gaussian weight over sin**2(colatitude)
1921
1922       REAL (KIND=HIGH) , DIMENSION(nlat) :: cosc , gwt , sinc , colat  , wos2
1923       REAL (KIND=HIGH) , PARAMETER       :: pi = 3.141592653589793
1924
1925       !  Convergence criterion for iteration of cos latitude
1926
1927       REAL , PARAMETER :: xlim  = 1.0E-14
1928
1929       INTEGER :: nzero, i, j
1930       REAL (KIND=HIGH) :: fi, fi1, a, b, g, gm, gp, gt, delta, c, d
1931
1932       !  The number of zeros between pole and equator
1933
1934       nzero = nlat/2
1935
1936       !  Set first guess for cos(colat)
1937
1938       DO i=1,nzero
1939          cosc(i) = SIN( (i-0.5)*pi/nlat + pi*0.5 )
1940       END DO
1941
1942       !  Constants for determining the derivative of the polynomial
1943       fi  = nlat
1944       fi1 = fi+1.0
1945       a   = fi*fi1 / SQRT(4.0*fi1*fi1-1.0)
1946       b   = fi1*fi / SQRT(4.0*fi*fi-1.0)
1947
1948       !  Loop over latitudes, iterating the search for each root
1949
1950       DO i=1,nzero
1951          j=0
1952
1953          !  Determine the value of the ordinary Legendre polynomial for
1954          !  the current guess root
1955
1956          DO
1957             CALL lgord( g, cosc(i), nlat )
1958
1959             !  Determine the derivative of the polynomial at this point
1960
1961             CALL lgord( gm, cosc(i), nlat-1 )
1962             CALL lgord( gp, cosc(i), nlat+1 )
1963             gt = (cosc(i)*cosc(i)-1.0) / (a*gp-b*gm)
1964
1965             !  Update the estimate of the root
1966
1967             delta   = g*gt
1968             cosc(i) = cosc(i) - delta
1969
1970             !  If convergence criterion has not been met, keep trying
1971
1972             j = j+1
1973             IF( ABS(delta).GT.xlim ) CYCLE
1974
1975             !  Determine the Gaussian weights
1976
1977             c      = 2.0 *( 1.0-cosc(i)*cosc(i) )
1978             CALL lgord( d, cosc(i), nlat-1 )
1979             d      = d*d*fi*fi
1980             gwt(i) = c *( fi-0.5 ) / d
1981             EXIT
1982
1983          END DO
1984
1985       END DO
1986
1987       !  Determine the colatitudes and sin(colat) and weights over sin**2
1988
1989       DO i=1,nzero
1990          colat(i)= ACOS(cosc(i))
1991          sinc(i) = SIN(colat(i))
1992          wos2(i) = gwt(i) /( sinc(i)*sinc(i) )
1993       END DO
1994
1995       !  If NLAT is odd, set values at the equator
1996
1997       IF( MOD(nlat,2) .NE. 0 ) THEN
1998          i       = nzero+1
1999          cosc(i) = 0.0
2000          c       = 2.0
2001          CALL lgord( d, cosc(i), nlat-1 )
2002          d       = d*d*fi*fi
2003          gwt(i)  = c *( fi-0.5 ) / d
2004          colat(i)= pi*0.5
2005          sinc(i) = 1.0
2006          wos2(i) = gwt(i)
2007       END IF
2008
2009       !  Determine the southern hemisphere values by symmetry
2010
2011       DO i=nlat-nzero+1,nlat
2012          cosc(i) =-cosc(nlat+1-i)
2013          gwt(i)  = gwt(nlat+1-i)
2014          colat(i)= pi-colat(nlat+1-i)
2015          sinc(i) = sinc(nlat+1-i)
2016          wos2(i) = wos2(nlat+1-i)
2017       END DO
2018
2019    END SUBROUTINE lggaus
2020
2021
2022    SUBROUTINE lgord( f, cosc, n )
2023
2024       IMPLICIT NONE
2025
2026       !  LGORD calculates the value of an ordinary Legendre polynomial at a
2027       !  specific latitude.
2028
2029       !  On entry:
2030       !     cosc - COS(colatitude)
2031       !     n      - the degree of the polynomial
2032
2033       !  On exit:
2034       !     f      - the value of the Legendre polynomial of degree N at
2035       !              latitude ASIN(cosc)
2036
2037       REAL (KIND=HIGH) :: s1, c4, a, b, fk, f, cosc, colat, c1, fn, ang
2038       INTEGER :: n, k
2039
2040       !  Determine the colatitude
2041
2042       colat = ACOS(cosc)
2043
2044       c1 = SQRT(2.0_HIGH)
2045       DO k=1,n
2046          c1 = c1 * SQRT( 1.0 - 1.0/(4*k*k) )
2047       END DO
2048
2049       fn = n
2050       ang= fn * colat
2051       s1 = 0.0
2052       c4 = 1.0
2053       a  =-1.0
2054       b  = 0.0
2055       DO k=0,n,2
2056          IF (k.eq.n) c4 = 0.5 * c4
2057          s1 = s1 + c4 * COS(ang)
2058          a  = a + 2.0
2059          b  = b + 1.0
2060          fk = k
2061          ang= colat * (fn-fk-2.0)
2062          c4 = ( a * (fn-b+1.0) / ( b * (fn+fn-a) ) ) * c4
2063       END DO
2064
2065       f = s1 * c1
2066
2067    END SUBROUTINE lgord
2068
2069
2070    SUBROUTINE llij_gauss (lat, lon, proj, i, j)
2071
2072       IMPLICIT NONE
2073
2074       REAL    , INTENT(IN)  :: lat, lon
2075       REAL    , INTENT(OUT) :: i, j
2076       TYPE (proj_info), INTENT(IN) :: proj
2077
2078       INTEGER :: n , n_low
2079       LOGICAL :: found = .FALSE.
2080       REAL    :: diff_1 , diff_nlat
2081
2082       !  The easy one first, get the x location.  The calling routine has already made
2083       !  sure that the necessary assumptions concerning the sign of the deltalon and the
2084       !  relative east/west'ness of the longitude and the starting longitude are consistent
2085       !  to allow this easy computation.
2086
2087       i = ( lon - proj%lon1 ) / proj%loninc + 1.
2088
2089       !  Since this is a global data set, we need to be concerned about wrapping the
2090       !  fields around the globe.
2091
2092 !      IF      ( ( proj%loninc .GT. 0 ) .AND. &
2093 !                ( FLOOR((lon-proj%lon1)/proj%loninc) + 1 .GE. proj%ixdim ) .AND. &
2094 !                ( lon + proj%loninc .GE. proj%lon1 + 360 ) ) THEN
2095 !! BUG: We may need to set proj%wrap, but proj is intent(in)
2096 !! WHAT IS THIS USED FOR?
2097 !!        proj%wrap = .TRUE.
2098 !         i = proj%ixdim
2099 !      ELSE IF ( ( proj%loninc .LT. 0 ) .AND. &
2100 !                ( FLOOR((lon-proj%lon1)/proj%loninc) + 1 .GE. proj%ixdim ) .AND. &
2101 !                ( lon + proj%loninc .LE. proj%lon1 - 360 ) ) THEN
2102 ! ! BUG: We may need to set proj%wrap, but proj is intent(in)
2103 ! ! WHAT IS THIS USED FOR?
2104 ! !        proj%wrap = .TRUE.
2105 !         i = proj%ixdim
2106 !      END IF
2107
2108       !  Yet another quicky test, can we find bounding values?  If not, then we may be
2109       !  dealing with putting data to a polar projection, so just give them them maximal
2110       !  value for the location.  This is an OK assumption for the interpolation across the
2111       !  top of the pole, given how close the longitude lines are.
2112
2113       IF ( ABS(lat) .GT. ABS(proj%gauss_lat(1)) ) THEN
2114
2115          diff_1    = lat - proj%gauss_lat(1)
2116          diff_nlat = lat - proj%gauss_lat(proj%nlat*2)
2117
2118          IF ( ABS(diff_1) .LT. ABS(diff_nlat) ) THEN
2119             j = 1
2120          ELSE
2121             j = proj%nlat*2
2122          END IF
2123
2124       !  If the latitude is between the two bounding values, we have to search and interpolate.
2125
2126       ELSE
2127
2128          DO n = 1 , proj%nlat*2 -1
2129             IF ( ( proj%gauss_lat(n) - lat ) * ( proj%gauss_lat(n+1) - lat ) .LE. 0 ) THEN
2130                found = .TRUE.
2131                n_low = n
2132                EXIT
2133             END IF
2134          END DO
2135
2136          !  Everything still OK?
2137
2138          IF ( .NOT. found ) THEN
2139             PRINT '(A)','Troubles in river city.  No bounding values of latitude found in the Gaussian routines.'
2140             call mprintf(.true.,ERROR,'Gee_no_bounding_lats_Gaussian')
2141          END IF
2142
2143          j = ( ( proj%gauss_lat(n_low) - lat                     ) * ( n_low + 1 ) + &
2144                ( lat                   - proj%gauss_lat(n_low+1) ) * ( n_low     ) ) / &
2145                ( proj%gauss_lat(n_low) - proj%gauss_lat(n_low+1) )
2146
2147       END IF
2148
2149    END SUBROUTINE llij_gauss
2150
2151 END MODULE map_utils