!dis
!dis Open Source License/Disclaimer, Forecast Systems Laboratory
!dis NOAA/OAR/FSL, 325 Broadway Boulder, CO 80305
!dis
!dis This software is distributed under the Open Source Definition,
!dis which may be found at http://www.opensource.org/osd.html.
!dis
!dis In particular, redistribution and use in source and binary forms,
!dis with or without modification, are permitted provided that the
!dis following conditions are met:
!dis
!dis - Redistributions of source code must retain this notice, this
!dis list of conditions and the following disclaimer.
!dis
!dis - Redistributions in binary form must provide access to this
!dis notice, this list of conditions and the following disclaimer, and
!dis the underlying source code.
!dis
!dis - All modifications to this software must be clearly documented,
!dis and are solely the responsibility of the agent making the
!dis modifications.
!dis
!dis - If significant modifications or enhancements are made to this
!dis software, the FSL Software Policy Manager
!dis (softwaremgr@fsl.noaa.gov) should be notified.
!dis
!dis THIS SOFTWARE AND ITS DOCUMENTATION ARE IN THE PUBLIC DOMAIN
!dis AND ARE FURNISHED "AS IS." THE AUTHORS, THE UNITED STATES
!dis GOVERNMENT, ITS INSTRUMENTALITIES, OFFICERS, EMPLOYEES, AND
!dis AGENTS MAKE NO WARRANTY, EXPRESS OR IMPLIED, AS TO THE USEFULNESS
!dis OF THE SOFTWARE AND DOCUMENTATION FOR ANY PURPOSE. THEY ASSUME
!dis NO RESPONSIBILITY (1) FOR THE USE OF THE SOFTWARE AND
!dis DOCUMENTATION; OR (2) TO PROVIDE TECHNICAL SUPPORT TO USERS.
!dis
!dis
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)
! Polar Stereographic (code = PROJ_PS)
!
! 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 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, or PROJ_PS
! 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)
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
IMPLICIT NONE
! Define some private constants
REAL, PRIVATE, PARAMETER :: pi = 3.1415927
REAL, PRIVATE, PARAMETER :: deg_per_rad = 180./pi
REAL, PRIVATE, PARAMETER :: rad_per_deg = pi / 180.
! Mean Earth Radius in m. The value below is consistent
! with NCEP's routines and grids.
! REAL, PUBLIC , PARAMETER :: earth_radius_m = 6371200. ! from Brent
REAL, PUBLIC , PARAMETER :: earth_radius_m = 6370000. ! same as MM5 system
REAL, PUBLIC , PARAMETER :: radians_per_degree = pi / 180.
! Define public parameters
! Projection codes for proj_info structure:
INTEGER, PUBLIC, PARAMETER :: PROJ_LATLON = 0
INTEGER, PUBLIC, PARAMETER :: PROJ_MERC = 1
INTEGER, PUBLIC, PARAMETER :: PROJ_LC = 3
INTEGER, PUBLIC, PARAMETER :: PROJ_PS = 5
! Define data structures to define various projections
TYPE proj_info
INTEGER :: code ! Integer code for projection type
REAL :: lat1 ! SW latitude (1,1) in degrees (-90->90N)
REAL :: lon1 ! SW longitude (1,1) in degrees (-180->180E)
REAL :: dx ! Grid spacing in meters at truelats, used
! only for ps, lc, and merc projections
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
LOGICAL :: init ! Flag to indicate if this struct is
! ready for use
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%dx = -999.9
proj%stdlon = -999.9
proj%truelat1 = -999.9
proj%truelat2 = -999.9
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%init = .FALSE.
END SUBROUTINE map_init
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE map_set(proj_code,lat1,lon1,knowni,knownj,dx,stdlon,truelat1,truelat2,proj)
! 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
REAL, INTENT(IN) :: lat1
REAL, INTENT(IN) :: lon1
REAL, INTENT(IN) :: dx
REAL, INTENT(IN) :: stdlon
REAL, INTENT(IN) :: truelat1
REAL, INTENT(IN) :: truelat2
REAL, INTENT(IN) :: knowni , knownj
TYPE(proj_info), INTENT(OUT) :: proj
! Local variables
! Executable code
! First, check for validity of mandatory variables in proj
IF ( ABS(lat1) .GT. 90. ) THEN
WRITE(0,'(A)') 'Latitude of origin corner required as follows:'
WRITE(0,'(A)') ' -90N <= lat1 < = 90.N'
STOP 'MAP_INIT'
ENDIF
IF ( ABS(lon1) .GT. 180.) THEN
WRITE(0,'(A)') 'Longitude of origin required as follows:'
WRITE(0,'(A)') ' -180E <= lon1 <= 180W'
STOP 'MAP_INIT'
ENDIF
IF ((dx .LE. 0.).AND.(proj_code .NE. PROJ_LATLON)) THEN
WRITE(0,'(A)') 'Require grid spacing (dx) in meters be positive!'
STOP 'MAP_INIT'
ENDIF
IF ((ABS(stdlon) .GT. 180.).AND.(proj_code .NE. PROJ_MERC)) THEN
WRITE(0,'(A)') 'Need orientation longitude (stdlon) as: '
WRITE(0,'(A)') ' -180E <= lon1 <= 180W'
STOP 'MAP_INIT'
ENDIF
IF (ABS(truelat1).GT.90.) THEN
WRITE(0,'(A)') 'Set true latitude 1 for all projections!'
STOP 'MAP_INIT'
ENDIF
CALL map_init(proj)
proj%code = proj_code
proj%lat1 = lat1
proj%lon1 = lon1
proj%knowni = knowni
proj%knownj = knownj
proj%dx = dx
proj%stdlon = stdlon
proj%truelat1 = truelat1
proj%truelat2 = truelat2
IF (proj%code .NE. PROJ_LATLON) THEN
proj%dx = dx
IF (truelat1 .LT. 0.) THEN
proj%hemi = -1.0
ELSE
proj%hemi = 1.0
ENDIF
proj%rebydx = earth_radius_m / dx
ENDIF
pick_proj: SELECT CASE(proj%code)
CASE(PROJ_PS)
WRITE(0,'(A)') 'Setting up POLAR STEREOGRAPHIC map...'
CALL set_ps(proj)
CASE(PROJ_LC)
WRITE(0,'(A)') 'Setting up LAMBERT CONFORMAL map...'
IF (ABS(proj%truelat2) .GT. 90.) THEN
WRITE(0,'(A)') 'Second true latitude not set, assuming a tangent'
WRITE(0,'(A,F10.3)') 'projection at truelat1: ', proj%truelat1
proj%truelat2=proj%truelat1
ENDIF
CALL set_lc(proj)
CASE (PROJ_MERC)
WRITE(0,'(A)') 'Setting up MERCATOR map...'
CALL set_merc(proj)
CASE (PROJ_LATLON)
WRITE(0,'(A)') 'Setting up CYLINDRICAL EQUIDISTANT LATLON map...'
! Convert lon1 to 0->360 notation
IF (proj%lon1 .LT. 0.) proj%lon1 = proj%lon1 + 360.
CASE DEFAULT
WRITE(0,'(A,I2)') 'Unknown projection code: ', proj%code
STOP 'MAP_INIT'
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
WRITE(0,'(A)') 'You have not called map_set for this projection!'
STOP '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)
i = i + proj%knowni - 1.0
j = j + proj%knownj - 1.0
CASE(PROJ_PS)
CALL llij_ps(lat,lon,proj,i,j)
CASE(PROJ_LC)
CALL llij_lc(lat,lon,proj,i,j)
i = i + proj%knowni - 1.0
j = j + proj%knownj - 1.0
CASE DEFAULT
WRITE(0,'(A,I2)') 'Unrecognized map projection code: ', proj%code
STOP 'LATLON_TO_IJ'
END SELECT
RETURN
END SUBROUTINE latlon_to_ij
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE ij_to_latlon(proj, ii, jj, 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) :: ii
REAL, INTENT(IN) :: jj
REAL, INTENT(OUT) :: lat
REAL, INTENT(OUT) :: lon
REAL :: i, j
IF (.NOT.proj%init) THEN
WRITE(0,'(A)') 'You have not called map_set for this projection!'
STOP 'IJ_TO_LATLON'
ENDIF
i = ii
j = jj
SELECT CASE (proj%code)
CASE (PROJ_LATLON)
CALL ijll_latlon(i, j, proj, lat, lon)
CASE (PROJ_MERC)
i = ii - proj%knowni + 1.0
j = jj - proj%knownj + 1.0
CALL ijll_merc(i, j, proj, lat, lon)
CASE (PROJ_PS)
CALL ijll_ps(i, j, proj, lat, lon)
CASE (PROJ_LC)
i = ii - proj%knowni + 1.0
j = jj - proj%knownj + 1.0
CALL ijll_lc(i, j, proj, lat, lon)
CASE DEFAULT
WRITE(0,'(A,I2)') 'Unrecognized map projection code: ', proj%code
STOP '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)
WRITE(0,'(A,2F10.1)')'Computed (I,J) of pole point: ',proj%polei,proj%polej
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_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)
WRITE(0,'(A,F8.6)') 'Computed cone factor: ', 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 = 1. - proj%hemi * proj%rsw * SIN(arg)
proj%polej = 1. + proj%rsw * COS(arg)
WRITE(0,'(A,2F10.3)') 'Computed pole i/j = ', proj%polei, proj%polej
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.
IF (proj%hemi .EQ. -1.) THEN
inew = -i + 2.
jnew = -j + 2.
ELSE
inew = i
jnew = j
ENDIF
! 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?)
IF (proj%hemi .EQ. -1.) THEN
i = 2. - i
j = 2. - j
ENDIF
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 / (earth_radius_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) :: lat
REAL, INTENT(IN) :: lon
TYPE(proj_info),INTENT(IN) :: proj
REAL,INTENT(OUT) :: i
REAL,INTENT(OUT) :: j
REAL :: deltalon
deltalon = lon - proj%lon1
IF (deltalon .LT. -180.) deltalon = deltalon + 360.
IF (deltalon .GT. 180.) deltalon = deltalon - 360.
i = 1. + (deltalon/(proj%dlon*deg_per_rad))
j = 1. + (ALOG(TAN(0.5*((lat + 90.) * rad_per_deg)))) / &
proj%dlon - proj%rsw
RETURN
END SUBROUTINE llij_merc
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE ijll_merc(i, j, proj, lat, lon)
! Compute the lat/lon from i/j for mercator projection
IMPLICIT NONE
REAL,INTENT(IN) :: i
REAL,INTENT(IN) :: j
TYPE(proj_info),INTENT(IN) :: proj
REAL, INTENT(OUT) :: lat
REAL, INTENT(OUT) :: lon
lat = 2.0*ATAN(EXP(proj%dlon*(proj%rsw + j-1.)))*deg_per_rad - 90.
lon = (i-1.)*proj%dlon*deg_per_rad + proj%lon1
IF (lon.GT.180.) lon = lon - 360.
IF (lon.LT.-180.) lon = lon + 360.
RETURN
END SUBROUTINE ijll_merc
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE llij_latlon(lat, lon, proj, i, j)
! Compute the i/j location of a lat/lon on a LATLON grid.
IMPLICIT NONE
REAL, INTENT(IN) :: lat
REAL, INTENT(IN) :: lon
TYPE(proj_info), INTENT(IN) :: proj
REAL, INTENT(OUT) :: i
REAL, INTENT(OUT) :: j
REAL :: deltalat
REAL :: deltalon
REAL :: lon360
REAL :: latinc
REAL :: loninc
! Extract the latitude and longitude increments for this grid
! (e.g., 2.5 deg for NCEP reanalysis) from the proj structure, where
! loninc is saved in the stdlon tag and latinc is saved in truelat1
latinc = proj%truelat1
loninc = proj%stdlon
! Compute deltalat and deltalon as the difference between the input
! lat/lon and the origin lat/lon
deltalat = lat - proj%lat1
! To account for issues around the dateline, convert the incoming
! longitudes to be 0->360.
IF (lon .LT. 0) THEN
lon360 = lon + 360.
ELSE
lon360 = lon
ENDIF
deltalon = lon360 - proj%lon1
! Compute i/j
i = deltalon/loninc + 1.
j = deltalat/latinc + 1.
RETURN
END SUBROUTINE llij_latlon
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE ijll_latlon(i, j, proj, lat, lon)
! Compute the lat/lon location of an i/j on a LATLON grid.
IMPLICIT NONE
REAL, INTENT(IN) :: i
REAL, INTENT(IN) :: j
TYPE(proj_info), INTENT(IN) :: proj
REAL, INTENT(OUT) :: lat
REAL, INTENT(OUT) :: lon
REAL :: deltalat
REAL :: deltalon
REAL :: lon360
REAL :: latinc
REAL :: loninc
! Extract the latitude and longitude increments for this grid
! (e.g., 2.5 deg for NCEP reanalysis) from the proj structure, where
! loninc is saved in the stdlon tag and latinc is saved in truelat1
latinc = proj%truelat1
loninc = proj%stdlon
! Compute deltalat and deltalon
deltalat = (j-1.)*latinc
deltalon = (i-1.)*loninc
lat = proj%lat1 + deltalat
lon = proj%lon1 + deltalon
IF ((ABS(lat) .GT. 90.).OR.(ABS(deltalon) .GT.360.)) THEN
! Off the earth for this grid
lat = -999.
lon = -999.
ELSE
lon = lon + 360.
lon = AMOD(lon,360.)
IF (lon .GT. 180.) lon = lon -360.
ENDIF
RETURN
END SUBROUTINE ijll_latlon
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
END MODULE map_utils