da_map_set.inc
References to this file elsewhere.
1 subroutine da_map_set(proj_code,lat1,lon1,knowni,knownj,dx,stdlon,truelat1,truelat2,proj)
2 ! Given a partially filled proj_info structure, this routine computes
3 ! polei, polej, rsw, and cone (if LC projection) to complete the
4 ! structure. This allows us to eliminate redundant calculations when
5 ! calling the coordinate conversion routines multiple times for the
6 ! same map.
7 ! This will generally be the first routine called when a user wants
8 ! to be able to use the coordinate conversion routines, and it
9 ! will call the appropriate subroutines based on the
10 ! proj%code which indicates which projection type this is.
11
12 implicit none
13
14 integer, intent(in) :: proj_code
15 real, intent(in) :: lat1
16 real, intent(in) :: lon1
17 real, intent(in) :: dx
18 real, intent(in) :: stdlon
19 real, intent(in) :: truelat1
20 real, intent(in) :: truelat2
21 real, intent(in) :: knowni , knownj
22 type(proj_info), intent(out) :: proj
23
24 if (trace_use_dull) call da_trace_entry("da_map_set")
25
26 ! First, check for validity of mandatory variables in proj
27 if (ABS(lat1) > 90.0) then
28 call da_error(__FILE__,__LINE__, &
29 (/"Latitude of origin corner required as follows: -90N <= lat1 < = 90N"/))
30 end if
31 if (ABS(lon1) > 180.0) then
32 call da_error(__FILE__,__LINE__, &
33 (/"Longitude of origin required as follows: -180E <= lon1 <= 180W"/))
34 end if
35 if ((dx .LE. 0.0).AND.(proj_code .NE. PROJ_LATLON)) then
36 call da_error(__FILE__,__LINE__, &
37 (/"Require grid spacing (dx) in meters be positive!"/))
38 end if
39 if ((ABS(stdlon) > 180.0).AND.(proj_code .NE. PROJ_MERC)) then
40 call da_error(__FILE__,__LINE__, &
41 (/"Need orientation longitude (stdlon) as: -180E <= lon1 <= 180W"/))
42 end if
43 if (ABS(truelat1)>90.0) then
44 call da_error(__FILE__,__LINE__, &
45 (/"Set true latitude 1 for all projections!"/))
46 end if
47
48 call da_map_init(proj)
49 proj%code = proj_code
50 proj%lat1 = lat1
51 proj%lon1 = lon1
52 proj%knowni = knowni
53 proj%knownj = knownj
54 proj%dx = dx
55 proj%stdlon = stdlon
56 proj%truelat1 = truelat1
57 proj%truelat2 = truelat2
58 if (proj%code .NE. PROJ_LATLON) then
59 proj%dx = dx
60 if (truelat1 < 0.0) then
61 proj%hemi = -1.0
62 else
63 proj%hemi = 1.0
64 end if
65 proj%rebydx = earth_radius_m / dx
66 end if
67 pick_proj: select case(proj%code)
68
69 case(PROJ_PS)
70 if (print_detail_map) then
71 write(unit=stdout,fmt='(A)') 'Setting up POLAR STEREOGRAPHIC map...'
72 end if
73 call da_set_ps(proj)
74
75 case(PROJ_LC)
76 if (print_detail_map) then
77 write(unit=stdout,fmt='(A)') 'Setting up LAMBERT CONFORMAL map...'
78 end if
79 if (ABS(proj%truelat2) > 90.0) then
80 if (print_detail_map) then
81 write(unit=stdout,fmt='(A)') 'Second true latitude not set, assuming a tangent'
82 write(unit=stdout,fmt='(A,F10.3)') 'projection at truelat1: ', proj%truelat1
83 end if
84 proj%truelat2=proj%truelat1
85 end if
86 call da_set_lc(proj)
87
88 case (PROJ_MERC)
89 if (print_detail_map) then
90 write(unit=stdout,fmt='(A)') 'Setting up MERCATOR map...'
91 end if
92 call da_set_merc(proj)
93
94 case (PROJ_LATLON)
95 if (print_detail_map) then
96 write(unit=stdout,fmt='(A)') 'Setting up CYLinDRICAL EQUIDISTANT LATLON map...'
97 end if
98 ! Convert lon1 to 0->360 notation
99 if (proj%lon1 < 0.0) proj%lon1 = proj%lon1 + 360.0
100
101 case default
102 write(unit=message(1),fmt='(A,I2)') 'Unknown projection code: ', proj%code
103 call da_error(__FILE__,__LINE__,message(1:1))
104 end select pick_proj
105 proj%init = .true.
106
107 if (trace_use_dull) call da_trace_exit("da_map_set")
108
109 end subroutine da_map_set
110
111