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