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