/WPS/geogrid/src/module_map_utils.F
FORTRAN Legacy | 2210 lines | 1165 code | 507 blank | 538 comment | 37 complexity | 9089163e9f96298a6c09e995f476fd9e MD5 | raw file
Possible License(s): AGPL-1.0
Large files files are truncated, but you can click here to view the full file
- MODULE map_utils
- ! Module that defines constants, data structures, and
- ! subroutines used to convert grid indices to lat/lon
- ! and vice versa.
- !
- ! SUPPORTED PROJECTIONS
- ! ---------------------
- ! Cylindrical Lat/Lon (code = PROJ_LATLON)
- ! Mercator (code = PROJ_MERC)
- ! Lambert Conformal (code = PROJ_LC)
- ! Gaussian (code = PROJ_GAUSS)
- ! Polar Stereographic (code = PROJ_PS)
- ! Rotated Lat/Lon (code = PROJ_ROTLL)
- !
- ! REMARKS
- ! -------
- ! The routines contained within were adapted from routines
- ! obtained from NCEP's w3 library. The original NCEP routines were less
- ! flexible (e.g., polar-stereo routines only supported truelat of 60N/60S)
- ! than what we needed, so modifications based on equations in Hoke, Hayes, and
- ! Renninger (AFGWC/TN/79-003) were added to improve the flexibility.
- ! Additionally, coding was improved to F90 standards and the routines were
- ! combined into this module.
- !
- ! ASSUMPTIONS
- ! -----------
- ! Grid Definition:
- ! For mercator, lambert conformal, and polar-stereographic projections,
- ! the routines within assume the following:
- !
- ! 1. Grid is dimensioned (i,j) where i is the East-West direction,
- ! positive toward the east, and j is the north-south direction,
- ! positive toward the north.
- ! 2. Origin is at (1,1) and is located at the southwest corner,
- ! regardless of hemispere.
- ! 3. Grid spacing (dx) is always positive.
- ! 4. Values of true latitudes must be positive for NH domains
- ! and negative for SH domains.
- !
- ! For the latlon and Gaussian projection, the grid origin may be at any
- ! of the corners, and the deltalat and deltalon values can be signed to
- ! account for this using the following convention:
- ! Origin Location Deltalat Sign Deltalon Sign
- ! --------------- ------------- -------------
- ! SW Corner + +
- ! NE Corner - -
- ! NW Corner - +
- ! SE Corner + -
- !
- ! Data Definitions:
- ! 1. Any arguments that are a latitude value are expressed in
- ! degrees north with a valid range of -90 -> 90
- ! 2. Any arguments that are a longitude value are expressed in
- ! degrees east with a valid range of -180 -> 180.
- ! 3. Distances are in meters and are always positive.
- ! 4. The standard longitude (stdlon) is defined as the longitude
- ! line which is parallel to the grid's y-axis (j-direction), along
- ! which latitude increases (NOT the absolute value of latitude, but
- ! the actual latitude, such that latitude increases continuously
- ! from the south pole to the north pole) as j increases.
- ! 5. One true latitude value is required for polar-stereographic and
- ! mercator projections, and defines at which latitude the
- ! grid spacing is true. For lambert conformal, two true latitude
- ! values must be specified, but may be set equal to each other to
- ! specify a tangent projection instead of a secant projection.
- !
- ! USAGE
- ! -----
- ! To use the routines in this module, the calling routines must have the
- ! following statement at the beginning of its declaration block:
- ! USE map_utils
- !
- ! The use of the module not only provides access to the necessary routines,
- ! but also defines a structure of TYPE (proj_info) that can be used
- ! to declare a variable of the same type to hold your map projection
- ! information. It also defines some integer parameters that contain
- ! the projection codes so one only has to use those variable names rather
- ! than remembering the acutal code when using them. The basic steps are
- ! as follows:
- !
- ! 1. Ensure the "USE map_utils" is in your declarations.
- ! 2. Declare the projection information structure as type(proj_info):
- ! TYPE(proj_info) :: proj
- ! 3. Populate your structure by calling the map_set routine:
- ! CALL map_set(code,lat1,lon1,knowni,knownj,dx,stdlon,truelat1,truelat2,proj)
- ! where:
- ! code (input) = one of PROJ_LATLON, PROJ_MERC, PROJ_LC, PROJ_PS,
- ! PROJ_GAUSS, or PROJ_ROTLL
- ! lat1 (input) = Latitude of grid origin point (i,j)=(1,1)
- ! (see assumptions!)
- ! lon1 (input) = Longitude of grid origin
- ! knowni (input) = origin point, x-location
- ! knownj (input) = origin point, y-location
- ! dx (input) = grid spacing in meters (ignored for LATLON projections)
- ! stdlon (input) = Standard longitude for PROJ_PS and PROJ_LC,
- ! deltalon (see assumptions) for PROJ_LATLON,
- ! ignored for PROJ_MERC
- ! truelat1 (input) = 1st true latitude for PROJ_PS, PROJ_LC, and
- ! PROJ_MERC, deltalat (see assumptions) for PROJ_LATLON
- ! truelat2 (input) = 2nd true latitude for PROJ_LC,
- ! ignored for all others.
- ! proj (output) = The structure of type (proj_info) that will be fully
- ! populated after this call
- !
- ! 4. Now that the proj structure is populated, you may call either
- ! of the following routines:
- !
- ! latlon_to_ij(proj, lat, lon, i, j)
- ! ij_to_latlon(proj, i, j, lat, lon)
- !
- ! It is incumbent upon the calling routine to determine whether or
- ! not the values returned are within your domain's bounds. All values
- ! of i, j, lat, and lon are REAL values.
- !
- !
- ! REFERENCES
- ! ----------
- ! Hoke, Hayes, and Renninger, "Map Preojections and Grid Systems for
- ! Meteorological Applications." AFGWC/TN-79/003(Rev), Air Weather
- ! Service, 1985.
- !
- ! NCAR MM5v3 Modeling System, REGRIDDER program, module_first_guess_map.F
- ! NCEP routines w3fb06, w3fb07, w3fb08, w3fb09, w3fb11, w3fb12
- !
- ! HISTORY
- ! -------
- ! 27 Mar 2001 - Original Version
- ! Brent L. Shaw, NOAA/FSL (CSU/CIRA)
- !
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- use constants_module
- use misc_definitions_module
- use module_debug
- ! Define some private constants
- INTEGER, PRIVATE, PARAMETER :: HIGH = 8
-
- TYPE proj_info
-
- INTEGER :: code ! Integer code for projection TYPE
- INTEGER :: nlat ! For Gaussian -- number of latitude points
- ! north of the equator
- INTEGER :: nlon !
- !
- INTEGER :: nxmin ! Starting x-coordinate of periodic, regular lat/lon dataset
- INTEGER :: nxmax ! Ending x-coordinate of periodic, regular lat/lon dataset
- INTEGER :: ixdim ! For Rotated Lat/Lon -- number of mass points
- ! in an odd row
- INTEGER :: jydim ! For Rotated Lat/Lon -- number of rows
- INTEGER :: stagger ! For Rotated Lat/Lon -- mass or velocity grid
- REAL :: phi ! For Rotated Lat/Lon -- domain half-extent in
- ! degrees latitude
- REAL :: lambda ! For Rotated Lat/Lon -- domain half-extend in
- ! degrees longitude
- REAL :: lat1 ! SW latitude (1,1) in degrees (-90->90N)
- REAL :: lon1 ! SW longitude (1,1) in degrees (-180->180E)
- REAL :: lat0 ! For Cassini, latitude of projection pole
- REAL :: lon0 ! For Cassini, longitude of projection pole
- REAL :: dx ! Grid spacing in meters at truelats, used
- ! only for ps, lc, and merc projections
- REAL :: dy ! Grid spacing in meters at truelats, used
- ! only for ps, lc, and merc projections
- REAL :: latinc ! Latitude increment for cylindrical lat/lon
- REAL :: loninc ! Longitude increment for cylindrical lat/lon
- ! also the lon increment for Gaussian grid
- REAL :: dlat ! Lat increment for lat/lon grids
- REAL :: dlon ! Lon increment for lat/lon grids
- REAL :: stdlon ! Longitude parallel to y-axis (-180->180E)
- REAL :: truelat1 ! First true latitude (all projections)
- REAL :: truelat2 ! Second true lat (LC only)
- REAL :: hemi ! 1 for NH, -1 for SH
- REAL :: cone ! Cone factor for LC projections
- REAL :: polei ! Computed i-location of pole point
- REAL :: polej ! Computed j-location of pole point
- REAL :: rsw ! Computed radius to SW corner
- REAL :: rebydx ! Earth radius divided by dx
- REAL :: knowni ! X-location of known lat/lon
- REAL :: knownj ! Y-location of known lat/lon
- REAL :: re_m ! Radius of spherical earth, meters
- REAL :: rho0 ! For Albers equal area
- REAL :: nc ! For Albers equal area
- REAL :: bigc ! For Albers equal area
- LOGICAL :: init ! Flag to indicate if this struct is
- ! ready for use
- LOGICAL :: wrap ! For Gaussian -- flag to indicate wrapping
- ! around globe?
- LOGICAL :: comp_ll ! Work in computational lat/lon space for Cassini
- REAL, POINTER, DIMENSION(:) :: gauss_lat ! Latitude array for Gaussian grid
-
- END TYPE proj_info
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- CONTAINS
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- SUBROUTINE map_init(proj)
- ! Initializes the map projection structure to missing values
-
- IMPLICIT NONE
- TYPE(proj_info), INTENT(INOUT) :: proj
-
- proj%lat1 = -999.9
- proj%lon1 = -999.9
- proj%lat0 = -999.9
- proj%lon0 = -999.9
- proj%dx = -999.9
- proj%dy = -999.9
- proj%latinc = -999.9
- proj%loninc = -999.9
- proj%stdlon = -999.9
- proj%truelat1 = -999.9
- proj%truelat2 = -999.9
- proj%phi = -999.9
- proj%lambda = -999.9
- proj%ixdim = -999
- proj%jydim = -999
- proj%stagger = HH
- proj%nlat = 0
- proj%nlon = 0
- proj%nxmin = 1
- proj%nxmax = 43200
- proj%hemi = 0.0
- proj%cone = -999.9
- proj%polei = -999.9
- proj%polej = -999.9
- proj%rsw = -999.9
- proj%knowni = -999.9
- proj%knownj = -999.9
- proj%re_m = EARTH_RADIUS_M
- proj%init = .FALSE.
- proj%wrap = .FALSE.
- proj%rho0 = 0.
- proj%nc = 0.
- proj%bigc = 0.
- proj%comp_ll = .FALSE.
- nullify(proj%gauss_lat)
-
- END SUBROUTINE map_init
- SUBROUTINE map_set(proj_code, proj, lat1, lon1, lat0, lon0, knowni, knownj, dx, dy, latinc, &
- loninc, stdlon, truelat1, truelat2, nlat, nlon, ixdim, jydim, nxmin, nxmax, &
- stagger, phi, lambda, r_earth)
- ! Given a partially filled proj_info structure, this routine computes
- ! polei, polej, rsw, and cone (if LC projection) to complete the
- ! structure. This allows us to eliminate redundant calculations when
- ! calling the coordinate conversion routines multiple times for the
- ! same map.
- ! This will generally be the first routine called when a user wants
- ! to be able to use the coordinate conversion routines, and it
- ! will call the appropriate subroutines based on the
- ! proj%code which indicates which projection type this is.
-
- IMPLICIT NONE
-
- ! Declare arguments
- INTEGER, INTENT(IN) :: proj_code
- INTEGER, INTENT(IN), OPTIONAL :: nlat
- INTEGER, INTENT(IN), OPTIONAL :: nlon
- INTEGER, INTENT(IN), OPTIONAL :: ixdim
- INTEGER, INTENT(IN), OPTIONAL :: jydim
- INTEGER, INTENT(IN), OPTIONAL :: nxmin
- INTEGER, INTENT(IN), OPTIONAL :: nxmax
- INTEGER, INTENT(IN), OPTIONAL :: stagger
- REAL, INTENT(IN), OPTIONAL :: latinc
- REAL, INTENT(IN), OPTIONAL :: loninc
- REAL, INTENT(IN), OPTIONAL :: lat1
- REAL, INTENT(IN), OPTIONAL :: lon1
- REAL, INTENT(IN), OPTIONAL :: lat0
- REAL, INTENT(IN), OPTIONAL :: lon0
- REAL, INTENT(IN), OPTIONAL :: dx
- REAL, INTENT(IN), OPTIONAL :: dy
- REAL, INTENT(IN), OPTIONAL :: stdlon
- REAL, INTENT(IN), OPTIONAL :: truelat1
- REAL, INTENT(IN), OPTIONAL :: truelat2
- REAL, INTENT(IN), OPTIONAL :: knowni
- REAL, INTENT(IN), OPTIONAL :: knownj
- REAL, INTENT(IN), OPTIONAL :: phi
- REAL, INTENT(IN), OPTIONAL :: lambda
- REAL, INTENT(IN), OPTIONAL :: r_earth
- TYPE(proj_info), INTENT(OUT) :: proj
- INTEGER :: iter
- REAL :: dummy_lon1
- REAL :: dummy_lon0
- REAL :: dummy_stdlon
-
- ! First, verify that mandatory parameters are present for the specified proj_code
- IF ( proj_code == PROJ_LC ) THEN
- IF ( .NOT.PRESENT(truelat1) .OR. &
- .NOT.PRESENT(truelat2) .OR. &
- .NOT.PRESENT(lat1) .OR. &
- .NOT.PRESENT(lon1) .OR. &
- .NOT.PRESENT(knowni) .OR. &
- .NOT.PRESENT(knownj) .OR. &
- .NOT.PRESENT(stdlon) .OR. &
- .NOT.PRESENT(dx) ) THEN
- PRINT '(A,I2)', 'The following are mandatory parameters for projection code : ', proj_code
- PRINT '(A)', ' truelat1, truelat2, lat1, lon1, knowni, knownj, stdlon, dx'
- call mprintf(.true.,ERROR,'MAP_INIT')
- END IF
- ELSE IF ( proj_code == PROJ_PS ) THEN
- IF ( .NOT.PRESENT(truelat1) .OR. &
- .NOT.PRESENT(lat1) .OR. &
- .NOT.PRESENT(lon1) .OR. &
- .NOT.PRESENT(knowni) .OR. &
- .NOT.PRESENT(knownj) .OR. &
- .NOT.PRESENT(stdlon) .OR. &
- .NOT.PRESENT(dx) ) THEN
- PRINT '(A,I2)', 'The following are mandatory parameters for projection code : ', proj_code
- PRINT '(A)', ' truelat1, lat1, lon1, knonwi, knownj, stdlon, dx'
- call mprintf(.true.,ERROR,'MAP_INIT')
- END IF
- ELSE IF ( proj_code == PROJ_PS_WGS84 ) THEN
- IF ( .NOT.PRESENT(truelat1) .OR. &
- .NOT.PRESENT(lat1) .OR. &
- .NOT.PRESENT(lon1) .OR. &
- .NOT.PRESENT(knowni) .OR. &
- .NOT.PRESENT(knownj) .OR. &
- .NOT.PRESENT(stdlon) .OR. &
- .NOT.PRESENT(dx) ) THEN
- PRINT '(A,I2)', 'The following are mandatory parameters for projection code : ', proj_code
- PRINT '(A)', ' truelat1, lat1, lon1, knonwi, knownj, stdlon, dx'
- call mprintf(.true.,ERROR,'MAP_INIT')
- END IF
- ELSE IF ( proj_code == PROJ_ALBERS_NAD83 ) THEN
- IF ( .NOT.PRESENT(truelat1) .OR. &
- .NOT.PRESENT(truelat2) .OR. &
- .NOT.PRESENT(lat1) .OR. &
- .NOT.PRESENT(lon1) .OR. &
- .NOT.PRESENT(knowni) .OR. &
- .NOT.PRESENT(knownj) .OR. &
- .NOT.PRESENT(stdlon) .OR. &
- .NOT.PRESENT(dx) ) THEN
- PRINT '(A,I2)', 'The following are mandatory parameters for projection code : ', proj_code
- PRINT '(A)', ' truelat1, truelat2, lat1, lon1, knonwi, knownj, stdlon, dx'
- call mprintf(.true.,ERROR,'MAP_INIT')
- END IF
- ELSE IF ( proj_code == PROJ_MERC ) THEN
- IF ( .NOT.PRESENT(truelat1) .OR. &
- .NOT.PRESENT(lat1) .OR. &
- .NOT.PRESENT(lon1) .OR. &
- .NOT.PRESENT(knowni) .OR. &
- .NOT.PRESENT(knownj) .OR. &
- .NOT.PRESENT(dx) ) THEN
- PRINT '(A,I2)', 'The following are mandatory parameters for projection code : ', proj_code
- PRINT '(A)', ' truelat1, lat1, lon1, knowni, knownj, dx'
- call mprintf(.true.,ERROR,'MAP_INIT')
- END IF
- ELSE IF ( proj_code == PROJ_LATLON ) THEN
- IF ( .NOT.PRESENT(latinc) .OR. &
- .NOT.PRESENT(loninc) .OR. &
- .NOT.PRESENT(knowni) .OR. &
- .NOT.PRESENT(knownj) .OR. &
- .NOT.PRESENT(lat1) .OR. &
- .NOT.PRESENT(lon1) ) THEN
- PRINT '(A,I2)', 'The following are mandatory parameters for projection code : ', proj_code
- PRINT '(A)', ' latinc, loninc, knowni, knownj, lat1, lon1'
- call mprintf(.true.,ERROR,'MAP_INIT')
- END IF
- ELSE IF ( proj_code == PROJ_CYL ) THEN
- IF ( .NOT.PRESENT(latinc) .OR. &
- .NOT.PRESENT(loninc) .OR. &
- .NOT.PRESENT(stdlon) ) THEN
- PRINT '(A,I2)', 'The following are mandatory parameters for projection code : ', proj_code
- PRINT '(A)', ' latinc, loninc, stdlon'
- call mprintf(.true.,ERROR,'MAP_INIT')
- END IF
- ELSE IF ( proj_code == PROJ_CASSINI ) THEN
- IF ( .NOT.PRESENT(latinc) .OR. &
- .NOT.PRESENT(loninc) .OR. &
- .NOT.PRESENT(lat1) .OR. &
- .NOT.PRESENT(lon1) .OR. &
- .NOT.PRESENT(lat0) .OR. &
- .NOT.PRESENT(lon0) .OR. &
- .NOT.PRESENT(knowni) .OR. &
- .NOT.PRESENT(knownj) .OR. &
- .NOT.PRESENT(stdlon) ) THEN
- PRINT '(A,I2)', 'The following are mandatory parameters for projection code : ', proj_code
- PRINT '(A)', ' latinc, loninc, lat1, lon1, knowni, knownj, lat0, lon0, stdlon'
- call mprintf(.true.,ERROR,'MAP_INIT')
- END IF
- ELSE IF ( proj_code == PROJ_GAUSS ) THEN
- IF ( .NOT.PRESENT(nlat) .OR. &
- .NOT.PRESENT(lat1) .OR. &
- .NOT.PRESENT(lon1) .OR. &
- .NOT.PRESENT(loninc) ) THEN
- PRINT '(A,I2)', 'The following are mandatory parameters for projection code : ', proj_code
- PRINT '(A)', ' nlat, lat1, lon1, loninc'
- call mprintf(.true.,ERROR,'MAP_INIT')
- END IF
- ELSE IF ( proj_code == PROJ_ROTLL ) THEN
- IF ( .NOT.PRESENT(ixdim) .OR. &
- .NOT.PRESENT(jydim) .OR. &
- .NOT.PRESENT(phi) .OR. &
- .NOT.PRESENT(lambda) .OR. &
- .NOT.PRESENT(lat1) .OR. &
- .NOT.PRESENT(lon1) .OR. &
- .NOT.PRESENT(stagger) ) THEN
- PRINT '(A,I2)', 'The following are mandatory parameters for projection code : ', proj_code
- PRINT '(A)', ' ixdim, jydim, phi, lambda, lat1, lon1, stagger'
- call mprintf(.true.,ERROR,'MAP_INIT')
- END IF
- ELSE
- PRINT '(A,I2)', 'Unknown projection code: ', proj_code
- call mprintf(.true.,ERROR,'MAP_INIT')
- END IF
-
- ! Check for validity of mandatory variables in proj
- IF ( PRESENT(lat1) ) THEN
- IF ( ABS(lat1) .GT. 90. ) THEN
- PRINT '(A)', 'Latitude of origin corner required as follows:'
- PRINT '(A)', ' -90N <= lat1 < = 90.N'
- call mprintf(.true.,ERROR,'MAP_INIT')
- ENDIF
- ENDIF
-
- IF ( PRESENT(lon1) ) THEN
- dummy_lon1 = lon1
- IF ( ABS(dummy_lon1) .GT. 180.) THEN
- iter = 0
- DO WHILE (ABS(dummy_lon1) > 180. .AND. iter < 10)
- IF (dummy_lon1 < -180.) dummy_lon1 = dummy_lon1 + 360.
- IF (dummy_lon1 > 180.) dummy_lon1 = dummy_lon1 - 360.
- iter = iter + 1
- END DO
- IF (abs(dummy_lon1) > 180.) THEN
- PRINT '(A)', 'Longitude of origin required as follows:'
- PRINT '(A)', ' -180E <= lon1 <= 180W'
- call mprintf(.true.,ERROR,'MAP_INIT')
- ENDIF
- ENDIF
- ENDIF
-
- IF ( PRESENT(lon0) ) THEN
- dummy_lon0 = lon0
- IF ( ABS(dummy_lon0) .GT. 180.) THEN
- iter = 0
- DO WHILE (ABS(dummy_lon0) > 180. .AND. iter < 10)
- IF (dummy_lon0 < -180.) dummy_lon0 = dummy_lon0 + 360.
- IF (dummy_lon0 > 180.) dummy_lon0 = dummy_lon0 - 360.
- iter = iter + 1
- END DO
- IF (abs(dummy_lon0) > 180.) THEN
- PRINT '(A)', 'Longitude of pole required as follows:'
- PRINT '(A)', ' -180E <= lon0 <= 180W'
- call mprintf(.true.,ERROR,'MAP_INIT')
- ENDIF
- ENDIF
- ENDIF
-
- IF ( PRESENT(dx) ) THEN
- IF ((dx .LE. 0.).AND.(proj_code .NE. PROJ_LATLON)) THEN
- PRINT '(A)', 'Require grid spacing (dx) in meters be positive!'
- call mprintf(.true.,ERROR,'MAP_INIT')
- ENDIF
- ENDIF
-
- IF ( PRESENT(stdlon) ) THEN
- dummy_stdlon = stdlon
- IF ((ABS(dummy_stdlon) > 180.).AND.(proj_code /= PROJ_MERC)) THEN
- iter = 0
- DO WHILE (ABS(dummy_stdlon) > 180. .AND. iter < 10)
- IF (dummy_stdlon < -180.) dummy_stdlon = dummy_stdlon + 360.
- IF (dummy_stdlon > 180.) dummy_stdlon = dummy_stdlon - 360.
- iter = iter + 1
- END DO
- IF (abs(dummy_stdlon) > 180.) THEN
- PRINT '(A)', 'Need orientation longitude (stdlon) as: '
- PRINT '(A)', ' -180E <= stdlon <= 180W'
- call mprintf(.true.,ERROR,'MAP_INIT')
- ENDIF
- ENDIF
- ENDIF
-
- IF ( PRESENT(truelat1) ) THEN
- IF (ABS(truelat1).GT.90.) THEN
- PRINT '(A)', 'Set true latitude 1 for all projections!'
- call mprintf(.true.,ERROR,'MAP_INIT')
- ENDIF
- ENDIF
-
- CALL map_init(proj)
- proj%code = proj_code
- IF ( PRESENT(lat1) ) proj%lat1 = lat1
- IF ( PRESENT(lon1) ) proj%lon1 = dummy_lon1
- IF ( PRESENT(lat0) ) proj%lat0 = lat0
- IF ( PRESENT(lon0) ) proj%lon0 = dummy_lon0
- IF ( PRESENT(latinc) ) proj%latinc = latinc
- IF ( PRESENT(loninc) ) proj%loninc = loninc
- IF ( PRESENT(knowni) ) proj%knowni = knowni
- IF ( PRESENT(knownj) ) proj%knownj = knownj
- IF ( PRESENT(nxmin) ) proj%nxmin = nxmin
- IF ( PRESENT(nxmax) ) proj%nxmax = nxmax
- IF ( PRESENT(dx) ) proj%dx = dx
- IF ( PRESENT(dy) ) THEN
- proj%dy = dy
- ELSE IF ( PRESENT(dx) ) THEN
- proj%dy = dx
- END IF
- IF ( PRESENT(stdlon) ) proj%stdlon = dummy_stdlon
- IF ( PRESENT(truelat1) ) proj%truelat1 = truelat1
- IF ( PRESENT(truelat2) ) proj%truelat2 = truelat2
- IF ( PRESENT(nlat) ) proj%nlat = nlat
- IF ( PRESENT(nlon) ) proj%nlon = nlon
- IF ( PRESENT(ixdim) ) proj%ixdim = ixdim
- IF ( PRESENT(jydim) ) proj%jydim = jydim
- IF ( PRESENT(stagger) ) proj%stagger = stagger
- IF ( PRESENT(phi) ) proj%phi = phi
- IF ( PRESENT(lambda) ) proj%lambda = lambda
- IF ( PRESENT(r_earth) ) proj%re_m = r_earth
-
- IF ( PRESENT(dx) ) THEN
- IF ( (proj_code == PROJ_LC) .OR. (proj_code == PROJ_PS) .OR. &
- (proj_code == PROJ_PS_WGS84) .OR. (proj_code == PROJ_ALBERS_NAD83) .OR. &
- (proj_code == PROJ_MERC) ) THEN
- proj%dx = dx
- IF (truelat1 .LT. 0.) THEN
- proj%hemi = -1.0
- ELSE
- proj%hemi = 1.0
- ENDIF
- proj%rebydx = proj%re_m / dx
- ENDIF
- ENDIF
- pick_proj: SELECT CASE(proj%code)
-
- CASE(PROJ_PS)
- CALL set_ps(proj)
- CASE(PROJ_PS_WGS84)
- CALL set_ps_wgs84(proj)
- CASE(PROJ_ALBERS_NAD83)
- CALL set_albers_nad83(proj)
-
- CASE(PROJ_LC)
- IF (ABS(proj%truelat2) .GT. 90.) THEN
- proj%truelat2=proj%truelat1
- ENDIF
- CALL set_lc(proj)
-
- CASE (PROJ_MERC)
- CALL set_merc(proj)
-
- CASE (PROJ_LATLON)
-
- CASE (PROJ_GAUSS)
- CALL set_gauss(proj)
-
- CASE (PROJ_CYL)
- CALL set_cyl(proj)
-
- CASE (PROJ_CASSINI)
- CALL set_cassini(proj)
-
- CASE (PROJ_ROTLL)
-
- END SELECT pick_proj
- proj%init = .TRUE.
- RETURN
- END SUBROUTINE map_set
- SUBROUTINE latlon_to_ij(proj, lat, lon, i, j)
- ! Converts input lat/lon values to the cartesian (i,j) value
- ! for the given projection.
-
- IMPLICIT NONE
- TYPE(proj_info), INTENT(IN) :: proj
- REAL, INTENT(IN) :: lat
- REAL, INTENT(IN) :: lon
- REAL, INTENT(OUT) :: i
- REAL, INTENT(OUT) :: j
-
- IF (.NOT.proj%init) THEN
- PRINT '(A)', 'You have not called map_set for this projection!'
- call mprintf(.true.,ERROR,'LATLON_TO_IJ')
- ENDIF
-
- SELECT CASE(proj%code)
-
- CASE(PROJ_LATLON)
- CALL llij_latlon(lat,lon,proj,i,j)
-
- CASE(PROJ_MERC)
- CALL llij_merc(lat,lon,proj,i,j)
-
- CASE(PROJ_PS)
- CALL llij_ps(lat,lon,proj,i,j)
- CASE(PROJ_PS_WGS84)
- CALL llij_ps_wgs84(lat,lon,proj,i,j)
-
- CASE(PROJ_ALBERS_NAD83)
- CALL llij_albers_nad83(lat,lon,proj,i,j)
-
- CASE(PROJ_LC)
- CALL llij_lc(lat,lon,proj,i,j)
-
- CASE(PROJ_GAUSS)
- CALL llij_gauss(lat,lon,proj,i,j)
-
- CASE(PROJ_CYL)
- CALL llij_cyl(lat,lon,proj,i,j)
- CASE(PROJ_CASSINI)
- CALL llij_cassini(lat,lon,proj,i,j)
- CASE(PROJ_ROTLL)
- CALL llij_rotlatlon(lat,lon,proj,i,j)
-
- CASE DEFAULT
- PRINT '(A,I2)', 'Unrecognized map projection code: ', proj%code
- call mprintf(.true.,ERROR,'LATLON_TO_IJ')
-
- END SELECT
- RETURN
- END SUBROUTINE latlon_to_ij
- SUBROUTINE ij_to_latlon(proj, i, j, lat, lon)
- ! Computes geographical latitude and longitude for a given (i,j) point
- ! in a grid with a projection of proj
-
- IMPLICIT NONE
- TYPE(proj_info),INTENT(IN) :: proj
- REAL, INTENT(IN) :: i
- REAL, INTENT(IN) :: j
- REAL, INTENT(OUT) :: lat
- REAL, INTENT(OUT) :: lon
-
- IF (.NOT.proj%init) THEN
- PRINT '(A)', 'You have not called map_set for this projection!'
- call mprintf(.true.,ERROR,'IJ_TO_LATLON')
- ENDIF
- SELECT CASE (proj%code)
-
- CASE (PROJ_LATLON)
- CALL ijll_latlon(i, j, proj, lat, lon)
-
- CASE (PROJ_MERC)
- CALL ijll_merc(i, j, proj, lat, lon)
-
- CASE (PROJ_PS)
- CALL ijll_ps(i, j, proj, lat, lon)
- CASE (PROJ_PS_WGS84)
- CALL ijll_ps_wgs84(i, j, proj, lat, lon)
-
- CASE (PROJ_ALBERS_NAD83)
- CALL ijll_albers_nad83(i, j, proj, lat, lon)
-
- CASE (PROJ_LC)
- CALL ijll_lc(i, j, proj, lat, lon)
-
- CASE (PROJ_CYL)
- CALL ijll_cyl(i, j, proj, lat, lon)
-
- CASE (PROJ_CASSINI)
- CALL ijll_cassini(i, j, proj, lat, lon)
-
- CASE (PROJ_ROTLL)
- CALL ijll_rotlatlon(i, j, proj, lat, lon)
-
- CASE DEFAULT
- PRINT '(A,I2)', 'Unrecognized map projection code: ', proj%code
- call mprintf(.true.,ERROR,'IJ_TO_LATLON')
-
- END SELECT
- RETURN
- END SUBROUTINE ij_to_latlon
- SUBROUTINE set_ps(proj)
- ! Initializes a polar-stereographic map projection from the partially
- ! filled proj structure. This routine computes the radius to the
- ! southwest corner and computes the i/j location of the pole for use
- ! in llij_ps and ijll_ps.
- IMPLICIT NONE
-
- ! Declare args
- TYPE(proj_info), INTENT(INOUT) :: proj
-
- ! Local vars
- REAL :: ala1
- REAL :: alo1
- REAL :: reflon
- REAL :: scale_top
-
- ! Executable code
- reflon = proj%stdlon + 90.
-
- ! Compute numerator term of map scale factor
- scale_top = 1. + proj%hemi * SIN(proj%truelat1 * rad_per_deg)
-
- ! Compute radius to lower-left (SW) corner
- ala1 = proj%lat1 * rad_per_deg
- proj%rsw = proj%rebydx*COS(ala1)*scale_top/(1.+proj%hemi*SIN(ala1))
-
- ! Find the pole point
- alo1 = (proj%lon1 - reflon) * rad_per_deg
- proj%polei = proj%knowni - proj%rsw * COS(alo1)
- proj%polej = proj%knownj - proj%hemi * proj%rsw * SIN(alo1)
- RETURN
- END SUBROUTINE set_ps
- SUBROUTINE llij_ps(lat,lon,proj,i,j)
- ! Given latitude (-90 to 90), longitude (-180 to 180), and the
- ! standard polar-stereographic projection information via the
- ! public proj structure, this routine returns the i/j indices which
- ! if within the domain range from 1->nx and 1->ny, respectively.
-
- IMPLICIT NONE
-
- ! Delcare input arguments
- REAL, INTENT(IN) :: lat
- REAL, INTENT(IN) :: lon
- TYPE(proj_info),INTENT(IN) :: proj
-
- ! Declare output arguments
- REAL, INTENT(OUT) :: i !(x-index)
- REAL, INTENT(OUT) :: j !(y-index)
-
- ! Declare local variables
-
- REAL :: reflon
- REAL :: scale_top
- REAL :: ala
- REAL :: alo
- REAL :: rm
-
- ! BEGIN CODE
-
- reflon = proj%stdlon + 90.
-
- ! Compute numerator term of map scale factor
-
- scale_top = 1. + proj%hemi * SIN(proj%truelat1 * rad_per_deg)
-
- ! Find radius to desired point
- ala = lat * rad_per_deg
- rm = proj%rebydx * COS(ala) * scale_top/(1. + proj%hemi *SIN(ala))
- alo = (lon - reflon) * rad_per_deg
- i = proj%polei + rm * COS(alo)
- j = proj%polej + proj%hemi * rm * SIN(alo)
-
- RETURN
- END SUBROUTINE llij_ps
- SUBROUTINE ijll_ps(i, j, proj, lat, lon)
-
- ! This is the inverse subroutine of llij_ps. It returns the
- ! latitude and longitude of an i/j point given the projection info
- ! structure.
-
- IMPLICIT NONE
-
- ! Declare input arguments
- REAL, INTENT(IN) :: i ! Column
- REAL, INTENT(IN) :: j ! Row
- TYPE (proj_info), INTENT(IN) :: proj
-
- ! Declare output arguments
- REAL, INTENT(OUT) :: lat ! -90 -> 90 north
- REAL, INTENT(OUT) :: lon ! -180 -> 180 East
-
- ! Local variables
- REAL :: reflon
- REAL :: scale_top
- REAL :: xx,yy
- REAL :: gi2, r2
- REAL :: arccos
-
- ! Begin Code
-
- ! Compute the reference longitude by rotating 90 degrees to the east
- ! to find the longitude line parallel to the positive x-axis.
- reflon = proj%stdlon + 90.
-
- ! Compute numerator term of map scale factor
- scale_top = 1. + proj%hemi * SIN(proj%truelat1 * rad_per_deg)
-
- ! Compute radius to point of interest
- xx = i - proj%polei
- yy = (j - proj%polej) * proj%hemi
- r2 = xx**2 + yy**2
-
- ! Now the magic code
- IF (r2 .EQ. 0.) THEN
- lat = proj%hemi * 90.
- lon = reflon
- ELSE
- gi2 = (proj%rebydx * scale_top)**2.
- lat = deg_per_rad * proj%hemi * ASIN((gi2-r2)/(gi2+r2))
- arccos = ACOS(xx/SQRT(r2))
- IF (yy .GT. 0) THEN
- lon = reflon + deg_per_rad * arccos
- ELSE
- lon = reflon - deg_per_rad * arccos
- ENDIF
- ENDIF
-
- ! Convert to a -180 -> 180 East convention
- IF (lon .GT. 180.) lon = lon - 360.
- IF (lon .LT. -180.) lon = lon + 360.
- RETURN
-
- END SUBROUTINE ijll_ps
- SUBROUTINE set_ps_wgs84(proj)
- ! Initializes a polar-stereographic map projection (WGS84 ellipsoid)
- ! from the partially filled proj structure. This routine computes the
- ! radius to the southwest corner and computes the i/j location of the
- ! pole for use in llij_ps and ijll_ps.
- IMPLICIT NONE
-
- ! Arguments
- TYPE(proj_info), INTENT(INOUT) :: proj
-
- ! Local variables
- real :: h, mc, tc, t, rho
- h = proj%hemi
- mc = cos(h*proj%truelat1*rad_per_deg)/sqrt(1.0-(E_WGS84*sin(h*proj%truelat1*rad_per_deg))**2.0)
- tc = sqrt(((1.0-sin(h*proj%truelat1*rad_per_deg))/(1.0+sin(h*proj%truelat1*rad_per_deg)))* &
- (((1.0+E_WGS84*sin(h*proj%truelat1*rad_per_deg))/(1.0-E_WGS84*sin(h*proj%truelat1*rad_per_deg)))**E_WGS84 ))
- ! Find the i/j location of reference lat/lon with respect to the pole of the projection
- t = sqrt(((1.0-sin(h*proj%lat1*rad_per_deg))/(1.0+sin(h*proj%lat1*rad_per_deg)))* &
- (((1.0+E_WGS84*sin(h*proj%lat1*rad_per_deg))/(1.0-E_WGS84*sin(h*proj%lat1*rad_per_deg)) )**E_WGS84 ) )
- rho = h * (A_WGS84 / proj%dx) * mc * t / tc
- proj%polei = rho * sin((h*proj%lon1 - h*proj%stdlon)*rad_per_deg)
- proj%polej = -rho * cos((h*proj%lon1 - h*proj%stdlon)*rad_per_deg)
- RETURN
- END SUBROUTINE set_ps_wgs84
- SUBROUTINE llij_ps_wgs84(lat,lon,proj,i,j)
- ! Given latitude (-90 to 90), longitude (-180 to 180), and the
- ! standard polar-stereographic projection information via the
- ! public proj structure, this routine returns the i/j indices which
- ! if within the domain range from 1->nx and 1->ny, respectively.
-
- IMPLICIT NONE
-
- ! Arguments
- REAL, INTENT(IN) :: lat
- REAL, INTENT(IN) :: lon
- REAL, INTENT(OUT) :: i !(x-index)
- REAL, INTENT(OUT) :: j !(y-index)
- TYPE(proj_info),INTENT(IN) :: proj
-
- ! Local variables
- real :: h, mc, tc, t, rho
- h = proj%hemi
- mc = cos(h*proj%truelat1*rad_per_deg)/sqrt(1.0-(E_WGS84*sin(h*proj%truelat1*rad_per_deg))**2.0)
- tc = sqrt(((1.0-sin(h*proj%truelat1*rad_per_deg))/(1.0+sin(h*proj%truelat1*rad_per_deg)))* &
- (((1.0+E_WGS84*sin(h*proj%truelat1*rad_per_deg))/(1.0-E_WGS84*sin(h*proj%truelat1*rad_per_deg)))**E_WGS84 ))
- t = sqrt(((1.0-sin(h*lat*rad_per_deg))/(1.0+sin(h*lat*rad_per_deg))) * &
- (((1.0+E_WGS84*sin(h*lat*rad_per_deg))/(1.0-E_WGS84*sin(h*lat*rad_per_deg)))**E_WGS84))
- ! Find the x/y location of the requested lat/lon with respect to the pole of the projection
- rho = (A_WGS84 / proj%dx) * mc * t / tc
- i = h * rho * sin((h*lon - h*proj%stdlon)*rad_per_deg)
- j = h *(-rho)* cos((h*lon - h*proj%stdlon)*rad_per_deg)
- ! Get i/j relative to reference i/j
- i = proj%knowni + (i - proj%polei)
- j = proj%knownj + (j - proj%polej)
-
- RETURN
- END SUBROUTINE llij_ps_wgs84
- SUBROUTINE ijll_ps_wgs84(i, j, proj, lat, lon)
-
- ! This is the inverse subroutine of llij_ps. It returns the
- ! latitude and longitude of an i/j point given the projection info
- ! structure.
-
- implicit none
-
- ! Arguments
- REAL, INTENT(IN) :: i ! Column
- REAL, INTENT(IN) :: j ! Row
- REAL, INTENT(OUT) :: lat ! -90 -> 90 north
- REAL, INTENT(OUT) :: lon ! -180 -> 180 East
- TYPE (proj_info), INTENT(IN) :: proj
- ! Local variables
- real :: h, mc, tc, t, rho, x, y
- real :: chi, a, b, c, d
- h = proj%hemi
- x = (i - proj%knowni + proj%polei)
- y = (j - proj%knownj + proj%polej)
- mc = cos(h*proj%truelat1*rad_per_deg)/sqrt(1.0-(E_WGS84*sin(h*proj%truelat1*rad_per_deg))**2.0)
- tc = sqrt(((1.0-sin(h*proj%truelat1*rad_per_deg))/(1.0+sin(h*proj%truelat1*rad_per_deg))) * &
- (((1.0+E_WGS84*sin(h*proj%truelat1*rad_per_deg))/(1.0-E_WGS84*sin(h*proj%truelat1*rad_per_deg)))**E_WGS84 ))
- rho = sqrt((x*proj%dx)**2.0 + (y*proj%dx)**2.0)
- t = rho * tc / (A_WGS84 * mc)
- lon = h*proj%stdlon + h*atan2(h*x,h*(-y))
- chi = PI/2.0-2.0*atan(t)
- a = 1./2.*E_WGS84**2. + 5./24.*E_WGS84**4. + 1./40.*E_WGS84**6. + 73./2016.*E_WGS84**8.
- b = 7./24.*E_WGS84**4. + 29./120.*E_WGS84**6. + 54113./40320.*E_WGS84**8.
- c = 7./30.*E_WGS84**6. + 81./280.*E_WGS84**8.
- d = 4279./20160.*E_WGS84**8.
- lat = chi + sin(2.*chi)*(a + cos(2.*chi)*(b + cos(2.*chi)*(c + d*cos(2.*chi))))
- lat = h * lat
- lat = lat*deg_per_rad
- lon = lon*deg_per_rad
- RETURN
-
- END SUBROUTINE ijll_ps_wgs84
- SUBROUTINE set_albers_nad83(proj)
- ! Initializes an Albers equal area map projection (NAD83 ellipsoid)
- ! from the partially filled proj structure. This routine computes the
- ! radius to the southwest corner and computes the i/j location of the
- ! pole for use in llij_albers_nad83 and ijll_albers_nad83.
- IMPLICIT NONE
-
- ! Arguments
- TYPE(proj_info), INTENT(INOUT) :: proj
-
- ! Local variables
- real :: h, m1, m2, q1, q2, theta, q, sinphi
- h = proj%hemi
- m1 = cos(h*proj%truelat1*rad_per_deg)/sqrt(1.0-(E_NAD83*sin(h*proj%truelat1*rad_per_deg))**2.0)
- m2 = cos(h*proj%truelat2*rad_per_deg)/sqrt(1.0-(E_NAD83*sin(h*proj%truelat2*rad_per_deg))**2.0)
- sinphi = sin(proj%truelat1*rad_per_deg)
- q1 = (1.0-E_NAD83**2.0) * &
- ((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)))
- sinphi = sin(proj%truelat2*rad_per_deg)
- q2 = (1.0-E_NAD83**2.0) * &
- ((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)))
- if (proj%truelat1 == proj%truelat2) then
- proj%nc = sin(proj%truelat1*rad_per_deg)
- else
- proj%nc = (m1**2.0 - m2**2.0) / (q2 - q1)
- end if
- proj%bigc = m1**2.0 + proj%nc*q1
- ! Find the i/j location of reference lat/lon with respect to the pole of the projection
- sinphi = sin(proj%lat1*rad_per_deg)
- q = (1.0-E_NAD83**2.0) * &
- ((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)))
- proj%rho0 = h * (A_NAD83 / proj%dx) * sqrt(proj%bigc - proj%nc * q) / proj%nc
- theta = proj%nc*(proj%lon1 - proj%stdlon)*rad_per_deg
- proj%polei = proj%rho0 * sin(h*theta)
- proj%polej = proj%rho0 - proj%rho0 * cos(h*theta)
- RETURN
- END SUBROUTINE set_albers_nad83
- SUBROUTINE llij_albers_nad83(lat,lon,proj,i,j)
- ! Given latitude (-90 to 90), longitude (-180 to 180), and the
- ! standard projection information via the
- ! public proj structure, this routine returns the i/j indices which
- ! if within the domain range from 1->nx and 1->ny, respectively.
-
- IMPLICIT NONE
-
- ! Arguments
- REAL, INTENT(IN) :: lat
- REAL, INTENT(IN) :: lon
- REAL, INTENT(OUT) :: i !(x-index)
- REAL, INTENT(OUT) :: j !(y-index)
- TYPE(proj_info),INTENT(IN) :: proj
-
- ! Local variables
- real :: h, q, rho, theta, sinphi
- h = proj%hemi
- sinphi = sin(h*lat*rad_per_deg)
- ! Find the x/y location of the requested lat/lon with respect to the pole of the projection
- q = (1.0-E_NAD83**2.0) * &
- ((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)))
- rho = h * (A_NAD83 / proj%dx) * sqrt(proj%bigc - proj%nc * q) / proj%nc
- theta = proj%nc * (h*lon - h*proj%stdlon)*rad_per_deg
- i = h*rho*sin(theta)
- j = h*proj%rho0 - h*rho*cos(theta)
- ! Get i/j relative to reference i/j
- i = proj%knowni + (i - proj%polei)
- j = proj%knownj + (j - proj%polej)
- RETURN
- END SUBROUTINE llij_albers_nad83
- SUBROUTINE ijll_albers_nad83(i, j, proj, lat, lon)
-
- ! This is the inverse subroutine of llij_albers_nad83. It returns the
- ! latitude and longitude of an i/j point given the projection info
- ! structure.
-
- implicit none
-
- ! Arguments
- REAL, INTENT(IN) :: i ! Column
- REAL, INTENT(IN) :: j ! Row
- REAL, INTENT(OUT) :: lat ! -90 -> 90 north
- REAL, INTENT(OUT) :: lon ! -180 -> 180 East
- TYPE (proj_info), INTENT(IN) :: proj
- ! Local variables
- real :: h, q, rho, theta, beta, x, y
- real :: a, b, c
- h = proj%hemi
- x = (i - proj%knowni + proj%polei)
- y = (j - proj%knownj + proj%polej)
- rho = sqrt(x**2.0 + (proj%rho0 - y)**2.0)
- theta = atan2(x, proj%rho0-y)
- q = (proj%bigc - (rho*proj%nc*proj%dx/A_NAD83)**2.0) / proj%nc
- beta = asin(q/(1.0 - log((1.0-E_NAD83)/(1.0+E_NAD83))*(1.0-E_NAD83**2.0)/(2.0*E_NAD83)))
- a = 1./3.*E_NAD83**2. + 31./180.*E_NAD83**4. + 517./5040.*E_NAD83**6.
- b = 23./360.*E_NAD83**4. + 251./3780.*E_NAD83**6.
- c = 761./45360.*E_NAD83**6.
- lat = beta + a*sin(2.*beta) + b*sin(4.*beta) + c*sin(6.*beta)
- lat = h*lat*deg_per_rad
- lon = proj%stdlon + theta*deg_per_rad/proj%nc
- RETURN
-
- END SUBROUTINE ijll_albers_nad83
- SUBROUTINE set_lc(proj)
- ! Initialize the remaining items in the proj structure for a
- ! lambert conformal grid.
-
- IMPLICIT NONE
-
- TYPE(proj_info), INTENT(INOUT) :: proj
-
- REAL :: arg
- REAL :: deltalon1
- REAL :: tl1r
- REAL :: ctl1r
-
- ! Compute cone factor
- CALL lc_cone(proj%truelat1, proj%truelat2, proj%cone)
-
- ! Compute longitude differences and ensure we stay out of the
- ! forbidden "cut zone"
- deltalon1 = proj%lon1 - proj%stdlon
- IF (deltalon1 .GT. +180.) deltalon1 = deltalon1 - 360.
- IF (deltalon1 .LT. -180.) deltalon1 = deltalon1 + 360.
-
- ! Convert truelat1 to radian and compute COS for later use
- tl1r = proj%truelat1 * rad_per_deg
- ctl1r = COS(tl1r)
-
- ! Compute the radius to our known lower-left (SW) corner
- proj%rsw = proj%rebydx * ctl1r/proj%cone * &
- (TAN((90.*proj%hemi-proj%lat1)*rad_per_deg/2.) / &
- TAN((90.*proj%hemi-proj%truelat1)*rad_per_deg/2.))**proj%cone
-
- ! Find pole point
- arg = proj%cone*(deltalon1*rad_per_deg)
- proj%polei = proj%hemi*proj%knowni - proj%hemi * proj%rsw * SIN(arg)
- proj%polej = proj%hemi*proj%knownj + proj%rsw * COS(arg)
-
- RETURN
- END SUBROUTINE set_lc
- SUBROUTINE lc_cone(truelat1, truelat2, cone)
-
- ! Subroutine to compute the cone factor of a Lambert Conformal projection
-
- IMPLICIT NONE
-
- ! Input Args
- REAL, INTENT(IN) :: truelat1 ! (-90 -> 90 degrees N)
- REAL, INTENT(IN) :: truelat2 ! " " " " "
-
- ! Output Args
- REAL, INTENT(OUT) :: cone
-
- ! Locals
-
- ! BEGIN CODE
-
- ! First, see if this is a secant or tangent projection. For tangent
- ! projections, truelat1 = truelat2 and the cone is tangent to the
- ! Earth's surface at this latitude. For secant projections, the cone
- ! intersects the Earth's surface at each of the distinctly different
- ! latitudes
- IF (ABS(truelat1-truelat2) .GT. 0.1) THEN
- cone = ALOG10(COS(truelat1*rad_per_deg)) - &
- ALOG10(COS(truelat2*rad_per_deg))
- cone = cone /(ALOG10(TAN((45.0 - ABS(truelat1)/2.0) * rad_per_deg)) - &
- ALOG10(TAN((45.0 - ABS(truelat2)/2.0) * rad_per_deg)))
- ELSE
- cone = SIN(ABS(truelat1)*rad_per_deg )
- ENDIF
- RETURN
- END SUBROUTINE lc_cone
- SUBROUTINE ijll_lc( i, j, proj, lat, lon)
-
- ! Subroutine to convert from the (i,j) cartesian coordinate to the
- ! geographical latitude and longitude for a Lambert Conformal projection.
-
- ! History:
- ! 25 Jul 01: Corrected by B. Shaw, NOAA/FSL
- !
- IMPLICIT NONE
-
- ! Input Args
- REAL, INTENT(IN) :: i ! Cartesian X coordinate
- REAL, INTENT(IN) :: j ! Cartesian Y coordinate
- TYPE(proj_info),INTENT(IN) :: proj ! Projection info structure
-
- ! Output Args
- REAL, INTENT(OUT) :: lat ! Latitude (-90->90 deg N)
- REAL, INTENT(OUT) :: lon ! Longitude (-180->180 E)
-
- ! Locals
- REAL :: inew
- REAL :: jnew
- REAL :: r
- REAL :: chi,chi1,chi2
- REAL :: r2
- REAL :: xx
- REAL :: yy
-
- ! BEGIN CODE
-
- chi1 = (90. - proj%hemi*proj%truelat1)*rad_per_deg
- chi2 = (90. - proj%hemi*proj%truelat2)*rad_per_deg
-
- ! See if we are in the southern hemispere and flip the indices
- ! if we are.
- inew = proj%hemi * i
- jnew = proj%hemi * j
-
- ! Compute radius**2 to i/j location
- xx = inew - proj%polei
- yy = proj%polej - jnew
- r2 = (xx*xx + yy*yy)
- r = SQRT(r2)/proj%rebydx
-
- ! Convert to lat/lon
- IF (r2 .EQ. 0.) THEN
- lat = proj%hemi * 90.
- lon = proj%stdlon
- ELSE
-
- ! Longitude
- lon = proj%stdlon + deg_per_rad * ATAN2(proj%hemi*xx,yy)/proj%cone
- lon = AMOD(lon+360., 360.)
-
- ! Latitude. Latitude determined by solving an equation adapted
- ! from:
- ! Maling, D.H., 1973: Coordinate Systems and Map Projections
- ! Equations #20 in Appendix I.
-
- IF (chi1 .EQ. chi2) THEN
- chi = 2.0*ATAN( ( r/TAN(chi1) )**(1./proj%cone) * TAN(chi1*0.5) )
- ELSE
- chi = 2.0*ATAN( (r*proj%cone/SIN(chi1))**(1./proj%cone) * TAN(chi1*0.5))
- ENDIF
- lat = (90.0-chi*deg_per_rad)*proj%hemi
-
- ENDIF
-
- IF (lon .GT. +180.) lon = lon - 360.
- IF (lon .LT. -180.) lon = lon + 360.
-
- RETURN
- END SUBROUTINE ijll_lc
- SUBROUTINE llij_lc( lat, lon, proj, i, j)
-
- ! Subroutine to compute the geographical latitude and longitude values
- ! to the cartesian x/y on a Lambert Conformal projection.
-
- IMPLICIT NONE
-
- ! Input Args
- REAL, INTENT(IN) :: lat ! Latitude (-90->90 deg N)
- REAL, INTENT(IN) :: lon ! Longitude (-180->180 E)
- TYPE(proj_info),INTENT(IN) :: proj ! Projection info structure
-
- ! Output Args
- REAL, INTENT(OUT) :: i ! Cartesian X coordinate
- REAL, INTENT(OUT) :: j ! Cartesian Y coordinate
-
- ! Locals
- REAL :: arg
- REAL :: deltalon
- REAL :: tl1r
- REAL :: rm
- REAL :: ctl1r
-
-
- ! BEGIN CODE
-
- ! Compute deltalon between known longitude and standard lon and ensure
- ! it is not in the cut zone
- deltalon = lon - proj%stdlon
- IF (deltalon .GT. +180.) deltalon = deltalon - 360.
- IF (deltalon .LT. -180.) deltalon = deltalon + 360.
-
- ! Convert truelat1 to radian and compute COS for later use
- tl1r = proj%truelat1 * rad_per_deg
- ctl1r = COS(tl1r)
-
- ! Radius to desired point
- rm = proj%rebydx * ctl1r/proj%cone * &
- (TAN((90.*proj%hemi-lat)*rad_per_deg/2.) / &
- TAN((90.*proj%hemi-proj%truelat1)*rad_per_deg/2.))**proj%cone
-
- arg = proj%cone*(deltalon*rad_per_deg)
- i = proj%polei + proj%hemi * rm * SIN(arg)
- j = proj%polej - rm * COS(arg)
-
- ! Finally, if we are in the southern hemisphere, flip the i/j
- ! values to a coordinate system where (1,1) is the SW corner
- ! (what we assume) which is different than the original NCEP
- ! algorithms which used the NE corner as the origin in the
- ! southern hemisphere (left-hand vs. right-hand coordinate?)
- i = proj%hemi * i
- j = proj%hemi * j
- RETURN
- END SUBROUTINE llij_lc
- SUBROUTINE set_merc(proj)
-
- ! Sets up the remaining basic elements for the mercator projection
-
- IMPLICIT NONE
- TYPE(proj_info), INTENT(INOUT) :: proj
- REAL :: clain
-
-
- ! Preliminary variables
-
- clain = COS(rad_per_deg*proj%truelat1)
- proj%dlon = proj%dx / (proj%re_m * clain)
-
- ! Compute distance from equator to origin, and store in the
- ! proj%rsw tag.
-
- proj%rsw = 0.
- IF (proj%lat1 .NE. 0.) THEN
- proj%rsw = (ALOG(TAN(0.5*((proj%lat1+90.)*rad_per_deg))))/proj%dlon
- ENDIF
- RETURN
- END SUBROUTINE set_merc
- SUBROUTINE llij_merc(lat, lon, proj, i, j)
-
- ! Compute i/j coordinate from lat lon for mercator projection
-
- IMPLICIT NONE
- REAL, INTENT(IN) …
Large files files are truncated, but you can click here to view the full file