Context Navigation

module_map_utils.f90 @ 897

Last change on this file since 897 was 207, checked in by aslmd, 14 years ago
MESOSCALE: A GENERAL CLEAN-UP FOLLOWING UPDATING THE USER MANUAL. EVERYTHING ESSENTIAL IS IN MESOSCALE (much lighter than before). EVERYTHING FOR DEVELOPPERS OR EXPERTS IS IN MESOSCALE_DEV.
File size: 40.3 KB

Line
1	!dis
2	!dis Open Source License/Disclaimer, Forecast Systems Laboratory
3	!dis NOAA/OAR/FSL, 325 Broadway Boulder, CO 80305
4	!dis
5	!dis This software is distributed under the Open Source Definition,
6	!dis which may be found at http://www.opensource.org/osd.html.
7	!dis
8	!dis In particular, redistribution and use in source and binary forms,
9	!dis with or without modification, are permitted provided that the
10	!dis following conditions are met:
11	!dis
12	!dis - Redistributions of source code must retain this notice, this
13	!dis list of conditions and the following disclaimer.
14	!dis
15	!dis - Redistributions in binary form must provide access to this
16	!dis notice, this list of conditions and the following disclaimer, and
17	!dis the underlying source code.
18	!dis
19	!dis - All modifications to this software must be clearly documented,
20	!dis and are solely the responsibility of the agent making the
21	!dis modifications.
22	!dis
23	!dis - If significant modifications or enhancements are made to this
24	!dis software, the FSL Software Policy Manager
25	!dis (softwaremgr@fsl.noaa.gov) should be notified.
26	!dis
27	!dis THIS SOFTWARE AND ITS DOCUMENTATION ARE IN THE PUBLIC DOMAIN
28	!dis AND ARE FURNISHED "AS IS." THE AUTHORS, THE UNITED STATES
29	!dis GOVERNMENT, ITS INSTRUMENTALITIES, OFFICERS, EMPLOYEES, AND
30	!dis AGENTS MAKE NO WARRANTY, EXPRESS OR IMPLIED, AS TO THE USEFULNESS
31	!dis OF THE SOFTWARE AND DOCUMENTATION FOR ANY PURPOSE. THEY ASSUME
32	!dis NO RESPONSIBILITY (1) FOR THE USE OF THE SOFTWARE AND
33	!dis DOCUMENTATION; OR (2) TO PROVIDE TECHNICAL SUPPORT TO USERS.
34	!dis
35	!dis
36
37	MODULE map_utils
38
39	! Module that defines constants, data structures, and
40	! subroutines used to convert grid indices to lat/lon
41	! and vice versa.
42	!
43	! SUPPORTED PROJECTIONS
44	! ---------------------
45	! Cylindrical Lat/Lon (code = PROJ_LATLON)
46	! Mercator (code = PROJ_MERC)
47	! Lambert Conformal (code = PROJ_LC)
48	! Polar Stereographic (code = PROJ_PS)
49	!
50	! REMARKS
51	! -------
52	! The routines contained within were adapted from routines
53	! obtained from the NCEP w3 library. The original NCEP routines were less
54	! flexible (e.g., polar-stereo routines only supported truelat of 60N/60S)
55	! than what we needed, so modifications based on equations in Hoke, Hayes, and
56	! Renninger (AFGWC/TN/79-003) were added to improve the flexibility.
57	! Additionally, coding was improved to F90 standards and the routines were
58	! combined into this module.
59	!
60	! ASSUMPTIONS
61	! -----------
62	! Grid Definition:
63	! For mercator, lambert conformal, and polar-stereographic projections,
64	! the routines within assume the following:
65	!
66	! 1. Grid is dimensioned (i,j) where i is the East-West direction,
67	! positive toward the east, and j is the north-south direction,
68	! positive toward the north.
69	! 2. Origin is at (1,1) and is located at the southwest corner,
70	! regardless of hemispere.
71	! 3. Grid spacing (dx) is always positive.
72	! 4. Values of true latitudes must be positive for NH domains
73	! and negative for SH domains.
74	!
75	! For the latlon projection, the grid origin may be at any of the
76	! corners, and the deltalat and deltalon values can be signed to
77	! account for this using the following convention:
78	! Origin Location Deltalat Sign Deltalon Sign
79	! --------------- ------------- -------------
80	! SW Corner + +
81	! NE Corner - -
82	! NW Corner - +
83	! SE Corner + -
84	!
85	! Data Definitions:
86	! 1. Any arguments that are a latitude value are expressed in
87	! degrees north with a valid range of -90 -> 90
88	! 2. Any arguments that are a longitude value are expressed in
89	! degrees east with a valid range of -180 -> 180.
90	! 3. Distances are in meters and are always positive.
91	! 4. The standard longitude (stdlon) is defined as the longitude
92	! line which is parallel to the y-axis (j-direction), along
93	! which latitude increases (NOT the absolute value of latitude, but
94	! the actual latitude, such that latitude increases continuously
95	! from the south pole to the north pole) as j increases.
96	! 5. One true latitude value is required for polar-stereographic and
97	! mercator projections, and defines at which latitude the
98	! grid spacing is true. For lambert conformal, two true latitude
99	! values must be specified, but may be set equal to each other to
100	! specify a tangent projection instead of a secant projection.
101	!
102	! USAGE
103	! -----
104	! To use the routines in this module, the calling routines must have the
105	! following statement at the beginning of its declaration block:
106	! USE map_utils
107	!
108	! The use of the module not only provides access to the necessary routines,
109	! but also defines a structure of TYPE (proj_info) that can be used
110	! to declare a variable of the same type to hold your map projection
111	! information. It also defines some integer parameters that contain
112	! the projection codes so one only has to use those variable names rather
113	! than remembering the acutal code when using them. The basic steps are
114	! as follows:
115	!
116	! 1. Ensure the "USE map_utils" is in your declarations.
117	! 2. Declare the projection information structure as type(proj_info):
118	! TYPE(proj_info) :: proj
119	! 3. Populate your structure by calling the map_set routine:
120	! CALL map_set(code,lat1,lon1,dx,stdlon,truelat1,truelat2,nx,ny,proj)
121	! where:
122	! code (input) = one of PROJ_LATLON, PROJ_MERC, PROJ_LC, or PROJ_PS
123	! lat1 (input) = Latitude of grid origin point (i,j)=(1,1)
124	! (see assumptions!)
125	! lon1 (input) = Longitude of grid origin
126	! dx (input) = grid spacing in meters (ignored for LATLON projections)
127	! stdlon (input) = Standard longitude for PROJ_PS and PROJ_LC,
128	! deltalon (see assumptions) for PROJ_LATLON,
129	! ignored for PROJ_MERC
130	! truelat1 (input) = 1st true latitude for PROJ_PS, PROJ_LC, and
131	! PROJ_MERC, deltalat (see assumptions) for PROJ_LATLON
132	! truelat2 (input) = 2nd true latitude for PROJ_LC,
133	! ignored for all others.
134	! nx = number of points in east-west direction
135	! ny = number of points in north-south direction
136	! proj (output) = The structure of type (proj_info) that will be fully
137	! populated after this call
138	!
139	! 4. Now that the proj structure is populated, you may call any
140	! of the following routines:
141	!
142	! latlon_to_ij(proj, lat, lon, i, j)
143	! ij_to_latlon(proj, i, j, lat, lon)
144	! truewind_to_gridwind(lon, proj, ugrid, vgrid, utrue, vtrue)
145	! gridwind_to_truewind(lon, proj, utrue, vtrue, ugrid, vgrid)
146	! compare_projections(proj1, proj2, same_proj)
147	!
148	! It is incumbent upon the calling routine to determine whether or
149	! not the values returned are within your domain bounds. All values
150	! of i, j, lat, and lon are REAL values.
151	!
152	!
153	! REFERENCES
154	! ----------
155	! Hoke, Hayes, and Renninger, "Map Preojections and Grid Systems for
156	! Meteorological Applications." AFGWC/TN-79/003(Rev), Air Weather
157	! Service, 1985.
158	!
159	! NCAR MM5v3 Modeling System, REGRIDDER program, module_first_guess_map.F
160	! NCEP routines w3fb06, w3fb07, w3fb08, w3fb09, w3fb11, w3fb12
161	!
162	! HISTORY
163	! -------
164	! 27 Mar 2001 - Original Version
165	! Brent L. Shaw, NOAA/FSL (CSU/CIRA)
166	! 02 Apr 2001 - Added routines to rotate winds from true to grid
167	! and vice versa.
168	! Brent L. Shaw, NOAA/FSL (CSU/CIRA)
169	! 09 Apr 2001 - Added compare_projections routine to compare two
170	! sets of projection parameters.
171	!
172	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
173
174	IMPLICIT NONE
175
176	! Define some private constants
177	REAL, PRIVATE, PARAMETER :: pi = 3.1415927
178	REAL, PRIVATE, PARAMETER :: deg_per_rad = 180./pi
179	REAL, PRIVATE, PARAMETER :: rad_per_deg = pi / 180.
180
181	!! Mean Earth Radius in m. The value below is consistent
182	!! with NCEPs routines and grids.
183	!REAL, PUBLIC, PARAMETER :: earth_radius_m = 6371200.
184	REAL, PUBLIC, PARAMETER :: earth_radius_m = 3397200.
185
186	! Define public parameters
187
188	! Projection codes for proj_info structure:
189	INTEGER, PUBLIC, PARAMETER :: PROJ_LATLON = 0
190	INTEGER, PUBLIC, PARAMETER :: PROJ_MERC = 3
191	INTEGER, PUBLIC, PARAMETER :: PROJ_LC = 1
192	INTEGER, PUBLIC, PARAMETER :: PROJ_PS = 2
193
194
195	! Define data structures to define various projections
196
197	TYPE proj_info
198
199	INTEGER :: code ! Integer code for projection type
200	REAL :: lat1 ! SW latitude (1,1) in degrees (-90->90N)
201	REAL :: lon1 ! SW longitude (1,1) in degrees (-180->180E)
202	REAL :: dx ! Grid spacing in meters at truelats, used
203	! only for ps, lc, and merc projections
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	LOGICAL :: init ! Flag to indicate if this struct is
216	! ready for use
217	INTEGER :: nx
218	INTEGER :: ny
219	END TYPE proj_info
220
221	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
222	CONTAINS
223	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
224	SUBROUTINE map_init(proj)
225	! Initializes the map projection structure to missing values
226
227	IMPLICIT NONE
228	TYPE(proj_info), INTENT(INOUT) :: proj
229
230	proj%lat1 = -999.9
231	proj%lon1 = -999.9
232	proj%dx = -999.9
233	proj%stdlon = -999.9
234	proj%truelat1 = -999.9
235	proj%truelat2 = -999.9
236	proj%hemi = 0.0
237	proj%cone = -999.9
238	proj%polei = -999.9
239	proj%polej = -999.9
240	proj%rsw = -999.9
241	proj%init = .FALSE.
242	proj%nx = -99
243	proj%ny = -99
244	END SUBROUTINE map_init
245	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
246	SUBROUTINE map_set(proj_code,lat1,lon1,dx,stdlon,truelat1,truelat2, &
247	idim,jdim,proj)
248	! Given a partially filled proj_info structure, this routine computes
249	! polei, polej, rsw, and cone (if LC projection) to complete the
250	! structure. This allows us to eliminate redundant calculations when
251	! calling the coordinate conversion routines multiple times for the
252	! same map.
253	! This will generally be the first routine called when a user wants
254	! to be able to use the coordinate conversion routines, and it
255	! will call the appropriate subroutines based on the
256	! proj%code which indicates which projection type this is.
257
258	IMPLICIT NONE
259
260	! Declare arguments
261	INTEGER, INTENT(IN) :: proj_code
262	REAL, INTENT(IN) :: lat1
263	REAL, INTENT(IN) :: lon1
264	REAL, INTENT(IN) :: dx
265	REAL, INTENT(IN) :: stdlon
266	REAL, INTENT(IN) :: truelat1
267	REAL, INTENT(IN) :: truelat2
268	INTEGER, INTENT(IN) :: idim
269	INTEGER, INTENT(IN) :: jdim
270	TYPE(proj_info), INTENT(OUT) :: proj
271
272	! Local variables
273
274
275	! Executable code
276
277	! First, check for validity of mandatory variables in proj
278	IF ( ABS(lat1) .GT. 90.001 ) THEN
279	PRINT '(A)', 'Latitude of origin corner required as follows:'
280	PRINT '(A)', ' -90N <= lat1 < = 90.N'
281	STOP 'MAP_INIT'
282	ENDIF
283	IF ( ABS(lon1) .GT. 180.) THEN
284	PRINT '(A)', 'Longitude of origin required as follows:'
285	PRINT '(A)', ' -180E <= lon1 <= 180W'
286	STOP 'MAP_INIT'
287	ENDIF
288	IF ((dx .LE. 0.).AND.(proj_code .NE. PROJ_LATLON)) THEN
289	PRINT '(A)', 'Require grid spacing (dx) in meters be positive!'
290	STOP 'MAP_INIT'
291	ENDIF
292	IF ((ABS(stdlon) .GT. 180.).AND.(proj_code .NE. PROJ_MERC)) THEN
293	PRINT '(A)', 'Need orientation longitude (stdlon) as: '
294	PRINT '(A)', ' -180E <= lon1 <= 180W'
295	STOP 'MAP_INIT'
296	ENDIF
297	IF (ABS(truelat1).GT.90.) THEN
298	PRINT '(A)', 'Set true latitude 1 for all projections!'
299	STOP 'MAP_INIT'
300	ENDIF
301
302	CALL map_init(proj)
303	proj%code = proj_code
304	proj%lat1 = lat1
305	proj%lon1 = lon1
306	proj%dx = dx
307	proj%stdlon = stdlon
308	proj%truelat1 = truelat1
309	proj%truelat2 = truelat2
310	proj%nx = idim
311	proj%ny = jdim
312	IF (proj%code .NE. PROJ_LATLON) THEN
313	proj%dx = dx
314	IF (truelat1 .LT. 0.) THEN
315	proj%hemi = -1.0
316	ELSE
317	proj%hemi = 1.0
318	ENDIF
319	proj%rebydx = earth_radius_m / dx
320	ENDIF
321	pick_proj: SELECT CASE(proj%code)
322
323	CASE(PROJ_PS)
324	!PRINT '(A)', 'Setting up POLAR STEREOGRAPHIC map...'
325	CALL set_ps(proj)
326
327	CASE(PROJ_LC)
328	!PRINT '(A)', 'Setting up LAMBERT CONFORMAL map...'
329	IF (ABS(proj%truelat2) .GT. 90.) THEN
330	PRINT '(A)', 'Second true latitude not set, assuming a tangent'
331	PRINT '(A,F10.3)', 'projection at truelat1: ', proj%truelat1
332	proj%truelat2=proj%truelat1
333	ELSE
334	! Ensure truelat1 < truelat2
335	proj%truelat1 = MIN(truelat1,truelat2)
336	proj%truelat2 = MAX(truelat1,truelat2)
337	ENDIF
338	CALL set_lc(proj)
339
340	CASE (PROJ_MERC)
341	!PRINT '(A)', 'Setting up MERCATOR map...'
342	CALL set_merc(proj)
343
344	CASE (PROJ_LATLON)
345	!PRINT '(A)', 'Setting up CYLINDRICAL EQUIDISTANT LATLON map...'
346	! Convert lon1 to 0->360 notation
347	IF (proj%lon1 .LT. 0.) proj%lon1 = proj%lon1 + 360.
348
349	CASE DEFAULT
350	PRINT '(A,I2)', 'Unknown projection code: ', proj%code
351	STOP 'MAP_INIT'
352
353	END SELECT pick_proj
354	proj%init = .TRUE.
355	RETURN
356	END SUBROUTINE map_set
357	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
358	SUBROUTINE latlon_to_ij(proj, lat, lon, i, j)
359	! Converts input lat/lon values to the cartesian (i,j) value
360	! for the given projection.
361
362	IMPLICIT NONE
363	TYPE(proj_info), INTENT(IN) :: proj
364	REAL, INTENT(IN) :: lat
365	REAL, INTENT(IN) :: lon
366	REAL, INTENT(OUT) :: i
367	REAL, INTENT(OUT) :: j
368
369	IF (.NOT.proj%init) THEN
370	PRINT '(A)', 'You have not called map_set for this projection!'
371	STOP 'LATLON_TO_IJ'
372	ENDIF
373
374	SELECT CASE(proj%code)
375
376	CASE(PROJ_LATLON)
377	CALL llij_latlon(lat,lon,proj,i,j)
378
379	CASE(PROJ_MERC)
380	CALL llij_merc(lat,lon,proj,i,j)
381
382	CASE(PROJ_PS)
383	CALL llij_ps(lat,lon,proj,i,j)
384
385	CASE(PROJ_LC)
386	CALL llij_lc(lat,lon,proj,i,j)
387
388	CASE DEFAULT
389	PRINT '(A,I2)', 'Unrecognized map projection code: ', proj%code
390	STOP 'LATLON_TO_IJ'
391
392	END SELECT
393	RETURN
394	END SUBROUTINE latlon_to_ij
395	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
396	SUBROUTINE ij_to_latlon(proj, i, j, lat, lon)
397	! Computes geographical latitude and longitude for a given (i,j) point
398	! in a grid with a projection of proj
399
400	IMPLICIT NONE
401	TYPE(proj_info),INTENT(IN) :: proj
402	REAL, INTENT(IN) :: i
403	REAL, INTENT(IN) :: j
404	REAL, INTENT(OUT) :: lat
405	REAL, INTENT(OUT) :: lon
406
407	IF (.NOT.proj%init) THEN
408	PRINT '(A)', 'You have not called map_set for this projection!'
409	STOP 'IJ_TO_LATLON'
410	ENDIF
411	SELECT CASE (proj%code)
412
413	CASE (PROJ_LATLON)
414	CALL ijll_latlon(i, j, proj, lat, lon)
415
416	CASE (PROJ_MERC)
417	CALL ijll_merc(i, j, proj, lat, lon)
418
419	CASE (PROJ_PS)
420	CALL ijll_ps(i, j, proj, lat, lon)
421
422	CASE (PROJ_LC)
423	CALL ijll_lc(i, j, proj, lat, lon)
424
425	CASE DEFAULT
426	PRINT '(A,I2)', 'Unrecognized map projection code: ', proj%code
427	STOP 'IJ_TO_LATLON'
428
429	END SELECT
430	RETURN
431	END SUBROUTINE ij_to_latlon
432	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
433	SUBROUTINE set_ps(proj)
434	! Initializes a polar-stereographic map projection from the partially
435	! filled proj structure. This routine computes the radius to the
436	! southwest corner and computes the i/j location of the pole for use
437	! in llij_ps and ijll_ps.
438	IMPLICIT NONE
439
440	! Declare args
441	TYPE(proj_info), INTENT(INOUT) :: proj
442
443	! Local vars
444	REAL :: ala1
445	REAL :: alo1
446	REAL :: reflon
447	REAL :: scale_top
448
449	! Executable code
450	reflon = proj%stdlon + 90.
451
452	! Cone factor
453	proj%cone = 1.0
454
455	! Compute numerator term of map scale factor
456	scale_top = 1. + proj%hemi * SIN(proj%truelat1 * rad_per_deg)
457
458	! Compute radius to lower-left (SW) corner
459	ala1 = proj%lat1 * rad_per_deg
460	proj%rsw = proj%rebydxCOS(ala1)scale_top/(1.+proj%hemi*SIN(ala1))
461
462	! Find the pole point
463	alo1 = (proj%lon1 - reflon) * rad_per_deg
464	proj%polei = 1. - proj%rsw * COS(alo1)
465	proj%polej = 1. - proj%hemi * proj%rsw * SIN(alo1)
466	PRINT '(A,2F10.1)', 'Computed (I,J) of pole point: ',proj%polei,proj%polej
467	RETURN
468	END SUBROUTINE set_ps
469	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
470	SUBROUTINE llij_ps(lat,lon,proj,i,j)
471	! Given latitude (-90 to 90), longitude (-180 to 180), and the
472	! standard polar-stereographic projection information via the
473	! public proj structure, this routine returns the i/j indices which
474	! if within the domain range from 1->nx and 1->ny, respectively.
475
476	IMPLICIT NONE
477
478	! Delcare input arguments
479	REAL, INTENT(IN) :: lat
480	REAL, INTENT(IN) :: lon
481	TYPE(proj_info),INTENT(IN) :: proj
482
483	! Declare output arguments
484	REAL, INTENT(OUT) :: i !(x-index)
485	REAL, INTENT(OUT) :: j !(y-index)
486
487	! Declare local variables
488
489	REAL :: reflon
490	REAL :: scale_top
491	REAL :: ala
492	REAL :: alo
493	REAL :: rm
494
495	! BEGIN CODE
496
497
498	reflon = proj%stdlon + 90.
499
500	! Compute numerator term of map scale factor
501
502	scale_top = 1. + proj%hemi * SIN(proj%truelat1 * rad_per_deg)
503
504	! Find radius to desired point
505	ala = lat * rad_per_deg
506	rm = proj%rebydx * COS(ala) * scale_top/(1. + proj%hemi *SIN(ala))
507	alo = (lon - reflon) * rad_per_deg
508	i = proj%polei + rm * COS(alo)
509	j = proj%polej + proj%hemi * rm * SIN(alo)
510
511	RETURN
512	END SUBROUTINE llij_ps
513	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
514	SUBROUTINE ijll_ps(i, j, proj, lat, lon)
515
516	! This is the inverse subroutine of llij_ps. It returns the
517	! latitude and longitude of an i/j point given the projection info
518	! structure.
519
520	IMPLICIT NONE
521
522	! Declare input arguments
523	REAL, INTENT(IN) :: i ! Column
524	REAL, INTENT(IN) :: j ! Row
525	TYPE (proj_info), INTENT(IN) :: proj
526
527	! Declare output arguments
528	REAL, INTENT(OUT) :: lat ! -90 -> 90 North
529	REAL, INTENT(OUT) :: lon ! -180 -> 180 East
530
531	! Local variables
532	REAL :: reflon
533	REAL :: scale_top
534	REAL :: xx,yy
535	REAL :: gi2, r2
536	REAL :: arccos
537
538	! Begin Code
539
540	! Compute the reference longitude by rotating 90 degrees to the east
541	! to find the longitude line parallel to the positive x-axis.
542	reflon = proj%stdlon + 90.
543
544	! Compute numerator term of map scale factor
545	scale_top = 1. + proj%hemi * SIN(proj%truelat1 * rad_per_deg)
546
547	! Compute radius to point of interest
548	xx = i - proj%polei
549	yy = (j - proj%polej) * proj%hemi
550	r2 = xx2 + yy2
551
552	! Now the magic code
553	IF (r2 .EQ. 0.) THEN
554	lat = proj%hemi * 90.
555	lon = reflon
556	ELSE
557	gi2 = (proj%rebydx * scale_top)**2.
558	lat = deg_per_rad * proj%hemi * ASIN((gi2-r2)/(gi2+r2))
559	arccos = ACOS(xx/SQRT(r2))
560	IF (yy .GT. 0) THEN
561	lon = reflon + deg_per_rad * arccos
562	ELSE
563	lon = reflon - deg_per_rad * arccos
564	ENDIF
565	ENDIF
566
567	! Convert to a -180 -> 180 East convention
568	IF (lon .GT. 180.) lon = lon - 360.
569	IF (lon .LT. -180.) lon = lon + 360.
570	RETURN
571
572	END SUBROUTINE ijll_ps
573	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
574	SUBROUTINE set_lc(proj)
575	! Initialize the remaining items in the proj structure for a
576	! lambert conformal grid.
577
578	IMPLICIT NONE
579
580	TYPE(proj_info), INTENT(INOUT) :: proj
581
582	REAL :: arg
583	REAL :: deltalon1
584	REAL :: tl1r
585	REAL :: ctl1r
586
587	! Compute cone factor
588	CALL lc_cone(proj%truelat1, proj%truelat2, proj%cone)
589	! PRINT '(A,F8.6)', 'Computed cone factor: ', proj%cone
590	! Compute longitude differences and ensure we stay out of the
591	! forbidden "cut zone"
592	deltalon1 = proj%lon1 - proj%stdlon
593	IF (deltalon1 .GT. +180.) deltalon1 = deltalon1 - 360.
594	IF (deltalon1 .LT. -180.) deltalon1 = deltalon1 + 360.
595
596	! Convert truelat1 to radian and compute COS for later use
597	tl1r = proj%truelat1 * rad_per_deg
598	ctl1r = COS(tl1r)
599
600	! Compute the radius to our known lower-left (SW) corner
601	proj%rsw = proj%rebydx * ctl1r/proj%cone * &
602	(TAN((90.proj%hemi-proj%lat1)rad_per_deg/2.) / &
603	TAN((90.proj%hemi-proj%truelat1)rad_per_deg/2.))**proj%cone
604
605	! Find pole point
606	arg = proj%cone(deltalon1rad_per_deg)
607	proj%polei = 1. - proj%hemi * proj%rsw * SIN(arg)
608	proj%polej = 1. + proj%rsw * COS(arg)
609	!PRINT '(A,2F10.3)', 'Computed pole i/j = ', proj%polei, proj%polej
610	RETURN
611	END SUBROUTINE set_lc
612	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
613	SUBROUTINE lc_cone(truelat1, truelat2, cone)
614
615	! Subroutine to compute the cone factor of a Lambert Conformal projection
616
617	IMPLICIT NONE
618
619	! Input Args
620	REAL, INTENT(IN) :: truelat1 ! (-90 -> 90 degrees N)
621	REAL, INTENT(IN) :: truelat2 ! " " " " "
622
623	! Output Args
624	REAL, INTENT(OUT) :: cone
625
626	! Locals
627
628	! BEGIN CODE
629
630	! First, see if this is a secant or tangent projection. For tangent
631	! projections, truelat1 = truelat2 and the cone is tangent to the
632	! Earth surface at this latitude. For secant projections, the cone
633	! intersects the Earth surface at each of the distinctly different
634	! latitudes
635	IF (ABS(truelat1-truelat2) .GT. 0.1) THEN
636
637	! Compute cone factor following:
638	cone=(ALOG(COS(truelat1rad_per_deg))-ALOG(COS(truelat2rad_per_deg))) / &
639	(ALOG(TAN((90.-ABS(truelat1))rad_per_deg0.5 ))- &
640	ALOG(TAN((90.-ABS(truelat2))rad_per_deg0.5 )) )
641	ELSE
642	cone = SIN(ABS(truelat1)*rad_per_deg )
643	ENDIF
644	RETURN
645	END SUBROUTINE lc_cone
646	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
647	SUBROUTINE ijll_lc( i, j, proj, lat, lon)
648
649	! Subroutine to convert from the (i,j) cartesian coordinate to the
650	! geographical latitude and longitude for a Lambert Conformal projection.
651
652	! History:
653	! 25 Jul 01: Corrected by B. Shaw, NOAA/FSL
654	!
655	IMPLICIT NONE
656
657	! Input Args
658	REAL, INTENT(IN) :: i ! Cartesian X coordinate
659	REAL, INTENT(IN) :: j ! Cartesian Y coordinate
660	TYPE(proj_info),INTENT(IN) :: proj ! Projection info structure
661
662	! Output Args
663	REAL, INTENT(OUT) :: lat ! Latitude (-90->90 deg N)
664	REAL, INTENT(OUT) :: lon ! Longitude (-180->180 E)
665
666	! Locals
667	REAL :: inew
668	REAL :: jnew
669	REAL :: r
670	REAL :: chi,chi1,chi2
671	REAL :: r2
672	REAL :: xx
673	REAL :: yy
674
675	! BEGIN CODE
676
677	chi1 = (90. - proj%hemiproj%truelat1)rad_per_deg
678	chi2 = (90. - proj%hemiproj%truelat2)rad_per_deg
679
680	! See if we are in the southern hemispere and flip the indices
681	! if we are.
682	IF (proj%hemi .EQ. -1.) THEN
683	inew = -i + 2.
684	jnew = -j + 2.
685	ELSE
686	inew = i
687	jnew = j
688	ENDIF
689
690	! Compute radius**2 to i/j location
691	xx = inew - proj%polei
692	yy = proj%polej - jnew
693	r2 = (xxxx + yyyy)
694	r = SQRT(r2)/proj%rebydx
695
696	! Convert to lat/lon
697	IF (r2 .EQ. 0.) THEN
698	lat = proj%hemi * 90.
699	lon = proj%stdlon
700	ELSE
701
702	! Longitude
703	lon = proj%stdlon + deg_per_rad * ATAN2(proj%hemi*xx,yy)/proj%cone
704	lon = AMOD(lon+360., 360.)
705
706	! Latitude. Latitude determined by solving an equation adapted
707	! from:
708	! Maling, D.H., 1973: Coordinate Systems and Map Projections
709	! Equations #20 in Appendix I.
710
711	IF (chi1 .EQ. chi2) THEN
712	chi = 2.0ATAN( ( r/TAN(chi1) )(1./proj%cone) TAN(chi1*0.5) )
713	ELSE
714	chi = 2.0ATAN( (rproj%cone/SIN(chi1))*(1./proj%cone) TAN(chi1*0.5))
715	ENDIF
716	lat = (90.0-chideg_per_rad)proj%hemi
717
718	ENDIF
719
720	IF (lon .GT. +180.) lon = lon - 360.
721	IF (lon .LT. -180.) lon = lon + 360.
722	RETURN
723	END SUBROUTINE ijll_lc
724	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
725	SUBROUTINE llij_lc( lat, lon, proj, i, j)
726
727	! Subroutine to compute the geographical latitude and longitude values
728	! to the cartesian x/y on a Lambert Conformal projection.
729
730	IMPLICIT NONE
731
732	! Input Args
733	REAL, INTENT(IN) :: lat ! Latitude (-90->90 deg N)
734	REAL, INTENT(IN) :: lon ! Longitude (-180->180 E)
735	TYPE(proj_info),INTENT(IN) :: proj ! Projection info structure
736
737	! Output Args
738	REAL, INTENT(OUT) :: i ! Cartesian X coordinate
739	REAL, INTENT(OUT) :: j ! Cartesian Y coordinate
740
741	! Locals
742	REAL :: arg
743	REAL :: deltalon
744	REAL :: tl1r
745	REAL :: rm
746	REAL :: ctl1r
747
748
749	! BEGIN CODE
750
751	! Compute deltalon between known longitude and standard lon and ensure
752	! it is not in the cut zone
753	deltalon = lon - proj%stdlon
754	IF (deltalon .GT. +180.) deltalon = deltalon - 360.
755	IF (deltalon .LT. -180.) deltalon = deltalon + 360.
756
757	! Convert truelat1 to radian and compute COS for later use
758	tl1r = proj%truelat1 * rad_per_deg
759	ctl1r = COS(tl1r)
760
761	! Radius to desired point
762	rm = proj%rebydx * ctl1r/proj%cone * &
763	(TAN((90.proj%hemi-lat)rad_per_deg/2.) / &
764	TAN((90.proj%hemi-proj%truelat1)rad_per_deg/2.))**proj%cone
765
766	arg = proj%cone(deltalonrad_per_deg)
767	i = proj%polei + proj%hemi * rm * SIN(arg)
768	j = proj%polej - rm * COS(arg)
769
770	! Finally, if we are in the southern hemisphere, flip the i/j
771	! values to a coordinate system where (1,1) is the SW corner
772	! (what we assume) which is different than the original NCEP
773	! algorithms which used the NE corner as the origin in the
774	! southern hemisphere (left-hand vs. right-hand coordinate?)
775	IF (proj%hemi .EQ. -1.) THEN
776	i = 2. - i
777	j = 2. - j
778	ENDIF
779	RETURN
780	END SUBROUTINE llij_lc
781	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
782	SUBROUTINE set_merc(proj)
783
784	! Sets up the remaining basic elements for the mercator projection
785
786	IMPLICIT NONE
787	TYPE(proj_info), INTENT(INOUT) :: proj
788	REAL :: clain
789
790
791	! Preliminary variables
792
793	clain = COS(rad_per_deg*proj%truelat1)
794	proj%dlon = proj%dx / (earth_radius_m * clain)
795
796	! Compute distance from equator to origin, and store in the
797	! proj%rsw tag.
798
799	proj%rsw = 0.
800	IF (proj%lat1 .NE. 0.) THEN
801	proj%rsw = (ALOG(TAN(0.5((proj%lat1+90.)rad_per_deg))))/proj%dlon
802	ENDIF
803	RETURN
804	END SUBROUTINE set_merc
805	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
806	SUBROUTINE llij_merc(lat, lon, proj, i, j)
807
808	! Compute i/j coordinate from lat lon for mercator projection
809
810	IMPLICIT NONE
811	REAL, INTENT(IN) :: lat
812	REAL, INTENT(IN) :: lon
813	TYPE(proj_info),INTENT(IN) :: proj
814	REAL,INTENT(OUT) :: i
815	REAL,INTENT(OUT) :: j
816	REAL :: deltalon
817
818	deltalon = lon - proj%lon1
819	IF (deltalon .LT. -180.) deltalon = deltalon + 360.
820	IF (deltalon .GT. 180.) deltalon = deltalon - 360.
821	i = 1. + (deltalon/(proj%dlon*deg_per_rad))
822	j = 1. + (ALOG(TAN(0.5((lat + 90.) rad_per_deg)))) / &
823	proj%dlon - proj%rsw
824
825	RETURN
826	END SUBROUTINE llij_merc
827	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
828	SUBROUTINE ijll_merc(i, j, proj, lat, lon)
829
830	! Compute the lat/lon from i/j for mercator projection
831
832	IMPLICIT NONE
833	REAL,INTENT(IN) :: i
834	REAL,INTENT(IN) :: j
835	TYPE(proj_info),INTENT(IN) :: proj
836	REAL, INTENT(OUT) :: lat
837	REAL, INTENT(OUT) :: lon
838
839
840	lat = 2.0ATAN(EXP(proj%dlon(proj%rsw + j-1.)))*deg_per_rad - 90.
841	lon = (i-1.)proj%dlondeg_per_rad + proj%lon1
842	IF (lon.GT.180.) lon = lon - 360.
843	IF (lon.LT.-180.) lon = lon + 360.
844	RETURN
845	END SUBROUTINE ijll_merc
846	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
847	SUBROUTINE llij_latlon(lat, lon, proj, i, j)
848
849	! Compute the i/j location of a lat/lon on a LATLON grid.
850	IMPLICIT NONE
851	REAL, INTENT(IN) :: lat
852	REAL, INTENT(IN) :: lon
853	TYPE(proj_info), INTENT(IN) :: proj
854	REAL, INTENT(OUT) :: i
855	REAL, INTENT(OUT) :: j
856
857	REAL :: deltalat
858	REAL :: deltalon
859	REAL :: lon360
860	REAL :: latinc
861	REAL :: loninc
862
863	! Extract the latitude and longitude increments for this grid
864	! (e.g., 2.5 deg for NCEP reanalysis) from the proj structure, where
865	! loninc is saved in the stdlon tag and latinc is saved in truelat1
866
867	latinc = proj%truelat1
868	loninc = proj%stdlon
869
870	! Compute deltalat and deltalon as the difference between the input
871	! lat/lon and the origin lat/lon
872
873	deltalat = lat - proj%lat1
874
875	! To account for issues around the dateline, convert the incoming
876	! longitudes to be 0->360.
877	IF (lon .LT. 0) THEN
878	lon360 = lon + 360.
879	ELSE
880	lon360 = lon
881	ENDIF
882	deltalon = lon360 - proj%lon1
883	IF (deltalon .LT. 0) deltalon = deltalon + 360.
884
885	! Compute i/j
886	i = deltalon/loninc + 1.
887	j = deltalat/latinc + 1.
888	RETURN
889	END SUBROUTINE llij_latlon
890	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
891	SUBROUTINE ijll_latlon(i, j, proj, lat, lon)
892
893	! Compute the lat/lon location of an i/j on a LATLON grid.
894	IMPLICIT NONE
895	REAL, INTENT(IN) :: i
896	REAL, INTENT(IN) :: j
897	TYPE(proj_info), INTENT(IN) :: proj
898	REAL, INTENT(OUT) :: lat
899	REAL, INTENT(OUT) :: lon
900
901	REAL :: deltalat
902	REAL :: deltalon
903	REAL :: lon360
904	REAL :: latinc
905	REAL :: loninc
906
907	! Extract the latitude and longitude increments for this grid
908	! (e.g., 2.5 deg for NCEP reanalysis) from the proj structure, where
909	! loninc is saved in the stdlon tag and latinc is saved in truelat1
910
911	latinc = proj%truelat1
912	loninc = proj%stdlon
913
914	! Compute deltalat and deltalon
915
916	deltalat = (j-1.)*latinc
917	deltalon = (i-1.)*loninc
918	lat = proj%lat1 + deltalat
919	lon = proj%lon1 + deltalon
920
921	IF ((ABS(lat) .GT. 90.).OR.(ABS(deltalon) .GT.360.)) THEN
922	! Off the earth for this grid
923	lat = -999.
924	lon = -999.
925	ELSE
926	lon = lon + 360.
927	lon = AMOD(lon,360.)
928	IF (lon .GT. 180.) lon = lon -360.
929	ENDIF
930
931	RETURN
932	END SUBROUTINE ijll_latlon
933	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
934	SUBROUTINE gridwind_to_truewind(lon,proj,ugrid,vgrid,utrue,vtrue)
935
936	! Subroutine to convert a wind from grid north to true north.
937
938	IMPLICIT NONE
939
940	! Arguments
941	REAL, INTENT(IN) :: lon ! Longitude of point in degrees
942	TYPE(proj_info),INTENT(IN):: proj ! Projection info structure
943	REAL, INTENT(IN) :: ugrid ! U-component, grid-relative
944	REAL, INTENT(IN) :: vgrid ! V-component, grid-relative
945	REAL, INTENT(OUT) :: utrue ! U-component, earth-relative
946	REAL, INTENT(OUT) :: vtrue ! V-component, earth-relative
947
948	! Locals
949	REAL :: alpha
950	REAL :: diff
951
952	IF ((proj%code .EQ. PROJ_PS).OR.(proj%code .EQ. PROJ_LC))THEN
953	diff = lon - proj%stdlon
954	IF (diff .GT. 180.) diff = diff - 360.
955	IF (diff .LT.-180.) diff = diff + 360.
956	alpha = diff * proj%cone * rad_per_deg * SIGN(1.,proj%truelat1)
957	utrue = vgrid * SIN(alpha) + ugrid * COS(alpha)
958	vtrue = vgrid * COS(alpha) - ugrid * SIN(alpha)
959	ELSEIF ((proj%code .EQ. PROJ_MERC).OR.(proj%code .EQ. PROJ_LATLON))THEN
960	utrue = ugrid
961	vtrue = vgrid
962	ELSE
963	PRINT '(A)', 'Unrecognized projection.'
964	STOP 'GRIDWIND_TO_TRUEWIND'
965	ENDIF
966
967	RETURN
968	END SUBROUTINE gridwind_to_truewind
969	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
970	SUBROUTINE truewind_to_gridwind(lon, proj, utrue, vtrue, ugrid, vgrid)
971
972	! Subroutine to compute grid-relative u/v wind components from the earth-
973	! relative values for a given projection.
974
975	IMPLICIT NONE
976
977	! Arguments
978	REAL, INTENT(IN) :: lon
979	TYPE(proj_info),INTENT(IN) :: proj
980	REAL, INTENT(IN) :: utrue
981	REAL, INTENT(IN) :: vtrue
982	REAL, INTENT(OUT) :: ugrid
983	REAL, INTENT(OUT) :: vgrid
984
985	! Locals
986	REAL :: alpha
987	REAL :: diff
988
989	IF ((proj%code .EQ. PROJ_PS).OR.(proj%code .EQ. PROJ_LC))THEN
990
991	diff = proj%stdlon - lon
992	IF (diff .GT. 180.) diff = diff - 360.
993	IF (diff .LT.-180.) diff = diff + 360.
994	alpha = diff * proj%cone * rad_per_deg * SIGN(1.,proj%truelat1)
995	ugrid = vtrue * SIN(alpha) + utrue * COS(alpha)
996	vgrid = vtrue * COS(alpha) - utrue * SIN(alpha)
997
998	ELSEIF((proj%code .EQ. PROJ_MERC).OR.(proj%code .EQ. PROJ_LATLON)) THEN
999	ugrid = utrue
1000	vgrid = vtrue
1001	ELSE
1002	PRINT '(A)', 'Unrecognized map projection.'
1003	STOP 'TRUEWIND_TO_GRIDWIND'
1004	ENDIF
1005	RETURN
1006	END SUBROUTINE truewind_to_gridwind
1007	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1008	SUBROUTINE compare_projections(proj1, proj2, same_proj)
1009
1010	! Subroutine to compare two proj_info structures to determine if the
1011	! maps are the same.
1012
1013	IMPLICIT NONE
1014
1015	! Arguments
1016	TYPE(proj_info), INTENT(IN) :: proj1
1017	TYPE(proj_info), INTENT(IN) :: proj2
1018	LOGICAL, INTENT(OUT) :: same_proj
1019
1020	! Locals
1021
1022	same_proj = .false.
1023
1024	! Make sure both structures are initialized
1025
1026	IF ((.NOT. proj1%init).OR.(.NOT.proj2%init)) THEN
1027	PRINT '(A)', 'COMPARE_PROJECTIONS: Map_Set not called yet!'
1028	same_proj = .false.
1029	ELSE
1030	same_proj = .true.
1031	ENDIF
1032
1033	! Check projection type
1034	IF (same_proj) THEN
1035	IF (proj1%code .NE. proj2%code) THEN
1036	PRINT '(A)', 'COMPARE_PROJECTIONS: Different projection type.'
1037	ELSE
1038	same_proj = .true.
1039	ENDIF
1040	ENDIF
1041
1042	! Check corner lat/lon
1043	IF (same_proj) THEN
1044	IF ( (ABS(proj1%lat1-proj2%lat1) .GT. 0.001) .OR. &
1045	(ABS(proj1%lon1-proj2%lon1) .GT. 0.001) ) THEN
1046	PRINT '(A)', 'COMPARE_PROJECTIONS: Different lat1/lon1'
1047	same_proj = .false.
1048	ENDIF
1049	ENDIF
1050
1051	! Compare dx
1052	IF ((same_proj).AND.(proj1%code .NE. PROJ_LATLON)) THEN
1053	IF ( ABS(proj1%dx-proj2%dx).GT.1. ) THEN
1054	PRINT '(A)', 'COMPARE_PROJECTIONS: Different dx'
1055	same_proj = .false.
1056	ENDIF
1057	ENDIF
1058	! Compare dimensions
1059	IF ((same_proj).AND.( (proj1%nx .NE. proj2%nx).OR. &
1060	(proj1%ny .NE. proj2%ny) ) ) THEN
1061	PRINT '(A)', 'COMPARE_PROJECTIONS: Different dimensions'
1062	same_proj = .false.
1063	ENDIF
1064
1065	IF ((same_proj).AND.(proj1%code .EQ. PROJ_LATLON)) THEN
1066	IF ( (proj1%dlat .NE. proj2%dlat).OR. &
1067	(proj1%dlon .NE. proj2%dlon)) THEN
1068	PRINT '(A)', 'COMPARE_PROJECTIONS: Different dlat/dlon'
1069	same_proj = .false.
1070	ENDIF
1071	ENDIF
1072
1073	! Compare stdlon for Polar and LC projections
1074	IF ( (same_proj).AND. ( (proj1%code .EQ. PROJ_PS).OR. &
1075	(proj1%code .EQ. PROJ_LC) ) ) THEN
1076	IF (proj1%stdlon .NE. proj2%stdlon) THEN
1077	PRINT '(A)', 'COMPARE_PROJECTIONS: Different stdlon.'
1078	same_proj = .false.
1079	ENDIF
1080	ENDIF
1081	! Compare true latitude1 if polar stereo, mercator, or lambert
1082	IF ( (same_proj).AND. ( (proj1%code .EQ. PROJ_PS) .OR. &
1083	(proj1%code .EQ. PROJ_LC) .OR. &
1084	(proj1%code .EQ. PROJ_MERC) ) ) THEN
1085	IF (proj1%truelat1 .NE. proj2%truelat1) THEN
1086	PRINT '(A)', 'COMPARE_PROJECTIONS: Different truelat1'
1087	same_proj = .false.
1088	ENDIF
1089	ENDIF
1090
1091	! Compare true latitude2 if LC
1092	IF ( (same_proj).AND.(proj1%code .EQ. PROJ_LC)) THEN
1093	IF (proj1%truelat2 .NE. proj2%truelat2) THEN
1094	PRINT '(A)', 'COMPARE_PROJECTIONS: Different truelat2'
1095	same_proj = .false.
1096	ENDIF
1097	ENDIF
1098
1099	RETURN
1100	END SUBROUTINE compare_projections
1101	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1102	SUBROUTINE compute_msf_lc(lat,truelat1,truelat2,msf)
1103
1104	! Computes the map scale factor for a Lambert Conformal grid at a given
1105	! latitude.
1106
1107	IMPLICIT NONE
1108	REAL, INTENT(IN) :: lat ! latitude where msf is requested
1109	REAL, INTENT(IN) :: truelat1
1110	REAL, INTENT(IN) :: truelat2
1111	REAL, INTENT(OUT) :: msf
1112
1113	REAL :: cone
1114	REAL :: psi1, psix, pole
1115
1116	CALL lc_cone(truelat1,truelat2, cone)
1117	IF (truelat1 .GE. 0.) THEN
1118	psi1 = (90. - truelat1) * rad_per_deg
1119	pole =90.
1120	ELSE
1121	psi1 = (90. + truelat1) * rad_per_deg
1122	pole = -90.
1123	ENDIF
1124	psix = (pole - lat)*rad_per_deg
1125	msf = (SIN(psi1)/SIN(psix)) &
1126	* ((TAN(psix0.5)/TAN(psi10.5))**cone)
1127	RETURN
1128	END SUBROUTINE compute_msf_lc
1129	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1130	SUBROUTINE compute_msf_ps(lat,truelat1,msf)
1131
1132	! Computes the map scale factor for a Polar Stereographic grid at a given
1133	! latitude.
1134
1135	IMPLICIT NONE
1136	REAL, INTENT(IN) :: lat ! latitude where msf is requested
1137	REAL, INTENT(IN) :: truelat1
1138	REAL, INTENT(OUT) :: msf
1139
1140	REAL :: psi1, psix, pole
1141
1142	IF (truelat1 .GE. 0.) THEN
1143	psi1 = (90. - truelat1) * rad_per_deg
1144	pole =90.
1145	ELSE
1146	psi1 = (90. + truelat1) * rad_per_deg
1147	pole = -90.
1148	ENDIF
1149	psix = (pole - lat)*rad_per_deg
1150	msf = ((1.+COS(psi1))/(1.0 + COS(psix)))
1151	RETURN
1152	END SUBROUTINE compute_msf_ps
1153	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1154	END MODULE map_utils

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

Download in other formats:

Original Format