| blob | 5609909296e4f3af6ecb57c6f8f018270a60c663 |
1 !---------------------------------------------------------------------------
2 ! Code copied from SI
3 !---------------------------------------------------------------------------
5 !dis
6 !dis Open Source License/Disclaimer, Forecast Systems Laboratory
7 !dis NOAA/OAR/FSL, 325 Broadway Boulder, CO 80305
8 !dis
9 !dis This software is distributed under the Open Source Definition,
10 !dis which may be found at http://www.opensource.org/osd.html.
11 !dis
12 !dis In particular, redistribution and use in source and binary forms,
13 !dis with or without modification, are permitted provided that the
14 !dis following conditions are met:
15 !dis
16 !dis - Redistributions of source code must retain this notice, this
17 !dis list of conditions and the following disclaimer.
18 !dis
19 !dis - Redistributions in binary form must provide access to this
20 !dis notice, this list of conditions and the following disclaimer, and
21 !dis the underlying source code.
22 !dis
23 !dis - All modifications to this software must be clearly documented,
24 !dis and are solely the responsibility of the agent making the
25 !dis modifications.
26 !dis
27 !dis - If significant modifications or enhancements are made to this
28 !dis software, the FSL Software Policy Manager
29 !dis (softwaremgr@fsl.noaa.gov) should be notified.
30 !dis
31 !dis THIS SOFTWARE AND ITS doCUMENTATION ARE in THE PUBLIC doMAin
32 !dis AND ARE FURNISHED "AS IS." THE AUTHORS, THE unitED STATES
33 !dis GOVERNMENT, ITS inSTRUMENTALITIES, OFFICERS, EMPLOYEES, AND
34 !dis AGENTS MAKE NO WARRANTY, EXPRESS OR IMPLIED, AS TO THE useFULNESS
35 !dis OF THE SOFTWARE AND doCUMENTATION FOR ANY Purpose. THEY ASsumE
36 !dis NO RESPONSIBILITY (1) FOR THE use OF THE SOFTWARE AND
37 !dis doCUMENTATION; OR (2) TO PROVIDE TECHNICAL SUPPORT TO useRS.
38 !dis
39 !dis
41 !module map_utils
43 ! Module that defines constants, data structures, and
44 ! subroutines used to convert grid indices to lat/lon
45 ! and vice versa.
46 !
47 ! SUPPORTED PROJECTIONS
48 ! ---------------------
49 ! Cylindrical Lat/Lon (code = PROJ_LATLON)
50 ! Mercator (code = PROJ_MERC)
51 ! Lambert Conformal (code = PROJ_LC)
52 ! Polar Stereographic (code = PROJ_PS)
53 !
54 ! REMARKS
55 ! -------
56 ! The routines contained within were adapted from routines
57 ! obtained from NCEP's w3 library. The original NCEP routines were less
58 ! flexible (e.g., polar-stereo routines only supported truelat of 60N/60S)
59 ! than what we needed, so modifications based on equations in Hoke, Hayes, and
60 ! Renninger (AFGWC/TN/79-003) were added to improve the flexibility.
61 ! Additionally, coding was improved to F90 standards and the routines were
62 ! combined into this module.
63 !
64 ! Assumptions
65 ! -----------
66 ! Grid Definition:
67 ! For mercator, lambert conformal, and polar-stereographic projections,
68 ! the routines within assume the following:
69 !
70 ! 1. Grid is dimensioned (i,j) where i is the East-West direction,
71 ! positive toward the east, and j is the north-south direction,
72 ! positive toward the north.
73 ! 2. Origin is at (1,1) and is located at the southwest corner,
74 ! regardless of hemispere.
75 ! 3. Grid spacing (dx) is always positive.
76 ! 4. Values of true latitudes must be positive for NH domains
77 ! and negative for SH domains.
78 !
79 ! For the latlon projection, the grid origin may be at any of the
80 ! corners, and the deltalat and deltalon values can be signed to
81 ! account for this using the following convention:
82 ! Origin Location Deltalat Sign Deltalon Sign
83 ! --------------- ------------- -------------
84 ! SW Corner + +
85 ! NE Corner - -
86 ! NW Corner - +
87 ! SE Corner + -
88 !
89 ! Data Definitions:
90 ! 1. Any arguments that are a latitude value are expressed in
91 ! degrees north with a valid range of -90 -> 90
92 ! 2. Any arguments that are a longitude value are expressed in
93 ! degrees east with a valid range of -180 -> 180
94 ! 3. Distances are in meters and are always positive.
95 ! 4. The standard longitude (stdlon) is defined as the longitude
96 ! line which is parallel to the grid's y-axis (j-direction), along
97 ! which latitude increases (NOT the absolute value of latitude, but
98 ! the actual latitude, such that latitude increases continuously
99 ! from the south pole to the north pole) as j increases.
100 ! 5. One true latitude value is required for polar-stereographic and
101 ! mercator projections, and defines at which latitude the
102 ! grid spacing is true. For lambert conformal, two true latitude
103 ! values must be specified, but may be set equal to each other to
104 ! specify a tangent projection instead of a secant projection.
105 !
106 ! USAGE
107 ! -----
108 ! To use the routines in this module, the calling routines must have the
109 ! following statement at the beginning of its declaration block:
110 ! use map_utils
111 !
112 ! The use of the module not only provides access to the necessary routines,
113 ! but also defines a structure of type (proj_info) that can be used
114 ! to declare a variable of the same type to hold your map projection
115 ! information. It also defines some integer parameters that contain
116 ! the projection codes so one only has to use those variable names rather
117 ! than remembering the acutal code when using them. The basic steps are
118 ! as follows:
119 !
120 ! 1. Ensure the "use map_utils" is in your declarations.
121 ! 2. Declare the projection information structure as type(proj_info):
122 ! type(proj_info) :: proj
123 ! 3. Populate your structure by calling the map_set routine:
124 ! call map_set(code,lat1,lon1,knowni,knownj,dx,stdlon,truelat1,truelat2,proj)
125 ! where:
126 ! code (input) = one of PROJ_LATLON, PROJ_MERC, PROJ_LC, or PROJ_PS
127 ! lat1 (input) = Latitude of grid origin point (i,j)=(1,1)
128 ! (see assumptions!)
129 ! lon1 (input) = Longitude of grid origin
130 ! knowni (input) = origin point, x-location
131 ! knownj (input) = origin point, y-location
132 ! dx (input) = grid spacing in meters (ignored for LATLON projections)
133 ! stdlon (input) = Standard longitude for PROJ_PS and PROJ_LC,
134 ! deltalon (see assumptions) for PROJ_LATLON,
135 ! ignored for PROJ_MERC
136 ! truelat1 (input) = 1st true latitude for PROJ_PS, PROJ_LC, and
137 ! PROJ_MERC, deltalat (see assumptions) for PROJ_LATLON
138 ! truelat2 (input) = 2nd true latitude for PROJ_LC,
139 ! ignored for all others.
140 ! proj (output) = The structure of type (proj_info) that will be fully
141 ! populated after this call
142 !
143 ! 4. Now that the proj structure is populated, you may call either
144 ! of the following routines:
145 !
146 ! da_latlon_to_ij(proj, lat, lon, i, j)
147 ! da_xyll(proj, i, j, lat, lon)
148 !
149 ! It is incumbent upon the calling routine to determine whether or
150 ! not the values returned are within your domain's bounds. All values
151 ! of i, j, lat, and lon are real values.
152 !
153 !
154 ! REFERENCES
155 ! ----------
156 ! Hoke, Hayes, and Renninger, "Map Preojections and Grid Systems for
157 ! Meteorological Applications." AFGWC/TN-79/003(Rev), Air Weather
158 ! Service, 1985.
159 !
160 ! NCAR MM5v3 Modeling System, REGRIDDER program, module_first_guess_map.F
161 ! NCEP routines w3fb06, w3fb07, w3fb08, w3fb09, w3fb11, w3fb12
162 !
163 ! HISTORY
164 ! -------
165 ! 27 Mar 2001 - Original Version
166 ! Brent L. Shaw, NOAA/FSL (CSU/CIRA)
167 !
168 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
170 ! Define some private constants
172 real, private, parameter :: deg_per_rad = 180.0/pi
173 real, private, parameter :: rad_per_deg = pi / 180.0
175 ! Mean Earth Radius in m. The value below is consistent
176 ! with NCEP's routines and grids.
177 ! real, public , parameter :: earth_radius_m = 6371200.0 ! from Brent
178 real, public , parameter :: earth_radius_m = 6370000.0
179 real, public , parameter :: radians_per_degree = pi / 180.0
181 ! Define public parameters
183 ! Projection codes for proj_info structure:
184 integer, public, parameter :: PROJ_LATLON = 0
185 integer, public, parameter :: PROJ_MERC = 1
186 integer, public, parameter :: PROJ_LC = 3
187 integer, public, parameter :: PROJ_PS = 5
190 ! Define data structures to define various projections
192 type proj_info
193 integer :: code ! integer code for projection type
194 real :: lat1 ! SW latitude (1,1) in degrees (-90->90N)
195 real :: lon1 ! SW longitude (1,1) in degrees (-180->180E)
196 real :: dx ! Grid spacing in meters at truelats, used
197 ! only for ps, lc, and merc projections
198 real :: dlat ! Lat increment for lat/lon grids
199 real :: dlon ! Lon increment for lat/lon grids
200 real :: stdlon ! Longitude parallel to y-axis (-180->180E)
201 real :: truelat1 ! First true latitude (all projections)
202 real :: truelat2 ! Second true lat (LC only)
203 real :: hemi ! 1 for NH, -1 for SH
204 real :: cone ! Cone factor for LC projections
205 real :: polei ! Computed i-location of pole point
206 real :: polej ! Computed j-location of pole point
207 real :: rsw ! Computed radius to SW corner
208 real :: rebydx ! Earth radius divided by dx
209 real :: knowni ! X-location of known lat/lon
210 real :: knownj ! Y-location of known lat/lon
211 real :: latinc ! Latitude increments in degrees
212 real :: loninc ! Longitude increments in degrees
213 logical :: init ! Flag to indicate if this struct is
214 ! ready for use
215 end type proj_info
217 type(proj_info) :: map_info, map_info_ens