share/module_llxy.F

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