Context Navigation

source: branches/publications/ORCHIDEE_IPSLCM5A2.1.r5307/ORCHIDEE/src_global/module_llxy.f90 @ 7683

Last change on this file since 7683 was 3447, checked in by josefine.ghattas, 8 years ago
Integrated branch ORCHIDEE-DRIVER in the trunk. See #244
File size: 73.8 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: