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 ! First, check for validity of mandatory variables in proj
25 if (ABS(lat1) > 90.) then
26 call da_error(__FILE__,__LINE__, &
27 (/"Latitude of origin corner required as follows: -90N <= lat1 < = 90.N"/))
28 end if
29 if (ABS(lon1) > 180.) then
30 call da_error(__FILE__,__LINE__, &
31 (/"Longitude of origin required as follows: -180E <= lon1 <= 180W"/))
32 end if
33 if ((dx .LE. 0.).AND.(proj_code .NE. PROJ_LATLON)) then
34 call da_error(__FILE__,__LINE__, &
35 (/"Require grid spacing (dx) in meters be positive!"/))
36 end if
37 if ((ABS(stdlon) > 180.).AND.(proj_code .NE. PROJ_MERC)) then
38 call da_error(__FILE__,__LINE__, &
39 (/"Need orientation longitude (stdlon) as: -180E <= lon1 <= 180W"/))
40 end if
41 if (ABS(truelat1)>90.) then
42 call da_error(__FILE__,__LINE__, &
43 (/"Set true latitude 1 for all projections!"/))
44 end if
45
46 call da_map_init(proj)
47 proj%code = proj_code
48 proj%lat1 = lat1
49 proj%lon1 = lon1
50 proj%knowni = knowni
51 proj%knownj = knownj
52 proj%dx = dx
53 proj%stdlon = stdlon
54 proj%truelat1 = truelat1
55 proj%truelat2 = truelat2
56 if (proj%code .NE. PROJ_LATLON) then
57 proj%dx = dx
58 if (truelat1 < 0.) then
59 proj%hemi = -1.0
60 else
61 proj%hemi = 1.0
62 end if
63 proj%rebydx = earth_radius_m / dx
64 end if
65 pick_proj: select case(proj%code)
66
67 case(PROJ_PS)
68 if (print_detail_map) then
69 write(unit=stdout,fmt='(A)') 'Setting up POLAR STEREOGRAPHIC map...'
70 end if
71 call da_set_ps(proj)
72
73 case(PROJ_LC)
74 if (print_detail_map) then
75 write(unit=stdout,fmt='(A)') 'Setting up LAMBERT CONFORMAL map...'
76 end if
77 if (ABS(proj%truelat2) > 90.) then
78 if (print_detail_map) then
79 write(unit=stdout,fmt='(A)') 'Second true latitude not set, assuming a tangent'
80 write(unit=stdout,fmt='(A,F10.3)') 'projection at truelat1: ', proj%truelat1
81 end if
82 proj%truelat2=proj%truelat1
83 end if
84 call da_set_lc(proj)
85
86 case (PROJ_MERC)
87 if (print_detail_map) then
88 write(unit=stdout,fmt='(A)') 'Setting up MERCATOR map...'
89 end if
90 call da_set_merc(proj)
91
92 case (PROJ_LATLON)
93 if (print_detail_map) then
94 write(unit=stdout,fmt='(A)') 'Setting up CYLinDRICAL EQUIDISTANT LATLON map...'
95 end if
96 ! Convert lon1 to 0->360 notation
97 if (proj%lon1 < 0.) proj%lon1 = proj%lon1 + 360.
98
99 case default
100 write(unit=message(1),fmt='(A,I2)') 'Unknown projection code: ', proj%code
101 call da_error(__FILE__,__LINE__,message(1:1))
102 end select pick_proj
103 proj%init = .true.
104
105 end subroutine da_map_set
106
107