geogrid/src/module_map_utils.F

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