1module regions_mod
2 ! MODULE regions_mod (prefix='reg' category='7. Low-level data objects')
3 !
4 !:Purpose: This is a subset of the EnKF module that deals with
5 ! operations involving the generation of a set of regions.
6 ! Each region consists of a small number of observations
7 ! which is to be assimilated simultaneously in the
8 ! sequential algorithm.
9 !
10 use mathPhysConstants_mod
11 use midasMpi_mod
12 implicit none
13
14 private
15 public :: struct_reg
16 public :: reg_getlatitude, reg_getblock, reg_locatestn, reg_init_struct
17
18 ! REG_MXOBS : maximum number of observations in a region
19 ! REG_MXOBS must be at least as large as the number of observations
20 ! in an individual radiosonde (121 currently)
21 integer, parameter :: REG_MXOBS = 900
22
23 type :: struct_reg
24 ! horizontal radius in km of the area in which observations are selected
25 real(8) :: r0_km
26 ! horizontal radius in km of the area in which an observation can have an impact
27 real(8) :: r1_km
28 ! the vertical extend in log(P) of the area in which an observation can have an
29 ! impact
30 real(8) :: rz_logp
31 ! as r0km but in radians and adjusted considering that we want an integer
32 ! number of latitude bands
33 real(8) :: r0_rad
34 real(8) :: r1_rad ! equivalent of r1 in radians
35 ! the horizontal length scale for the hadamard product (0.5 r1)
36 real(8) :: scale_hor
37 ! the vertical length scale for the hadamard product (in log (p))
38 real(8) :: scale_ver
39 ! number of latitude bands
40 integer :: nlatband
41 ! estimated maximum number of regions that could possibly be formed
42 integer :: mxreg
43 end type struct_reg
44
45contains
46
47 subroutine reg_getblock(nlatband, r0, latmin, latmax, nlonblock)
48 !
49 !:Purpose: We have nlatband latitude bands (bounded by latitudes
50 ! latmin and latmax). For each band determine how many
51 ! longitidunal blocks of at most r0*root(2) radians it
52 ! contains. Output is in nlonblock.
53 !
54 implicit none
55
56 ! Arguments:
57 integer, intent(in) :: nlatband
58 real(8), intent(in) :: latmin(nlatband)
59 real(8), intent(in) :: latmax(nlatband)
60 real(8), intent(in) :: r0
61 integer, intent(out) :: nlonblock(nlatband)
62
63 ! Locals:
64 integer :: i
65 real(8) :: bar, dpie, lat
66
67 dpie=MPC_PI_R8
68
69 nlonblock(1)=1
70 nlonblock(nlatband)=1
71
72 do i = 2, nlatband-1
73 if (latmin(i)*latmax(i).gt.0) then
74 lat = min(abs(latmin(i)),abs(latmax(i)))
75 else
76 ! at the equator
77 lat = 0
78 end if
79 bar = (2**0.5)*r0/cos(lat)
80 nlonblock(i) = int((2.0D0*dpie)/bar)+1
81 end do
82
83 end subroutine reg_getblock
84
85
86 subroutine reg_getlatitude(r0, nlatband, latmin, latcenter, latmax)
87 !
88 !:Purpose: the circle is covered with nlatband latitude bands.
89 ! the polar caps (with radius r0 radians) form the first and
90 ! last band. Intermediate bands are of width r0*(2**0.5)
91 ! For reach band we have the southern most latitude latmin,
92 ! the central latitude latcenter (at the poles for the extreme
93 ! bands) and the northern most latitude latcenter
94 !
95 implicit none
96
97 ! Arguments:
98 integer, intent(in) :: nlatband
99 real(8), intent(in) :: r0
100 real(8), intent(out) :: latmin(nlatband)
101 real(8), intent(out) :: latcenter(nlatband)
102 real(8), intent(out) :: latmax(nlatband)
103
104 ! Locals:
105 integer :: i
106 real(8) :: dpie
107
108 dpie = MPC_PI_R8
109
110 latmin(1) = -0.5D0*dpie
111 latcenter(1) = -0.5D0*dpie
112 latmax(1) = -0.5D0*dpie+r0
113
114 if (mmpi_myid == 0) write(*,*) 'get latitude r0 and nlatband: ', r0, nlatband
115 if (mmpi_myid == 0) write(*,*) 'at ', 1, latmin(1), latcenter(1), latmax(1)
116
117 latmin(nlatband) = 0.5D0*dpie-r0
118 latcenter(nlatband) = 0.5D0*dpie
119 latmax(nlatband) = 0.5D0*dpie
120
121 if (mmpi_myid == 0) write(*,*) 'at ', nlatband, latmin(nlatband), latcenter(nlatband), latmax(nlatband)
122 do i = 2, nlatband-1
123 latmin(i) = latmax(i-1)
124 latmax(i) = latmin(i)+r0*(2**0.5)
125 latcenter(i) = latmin(i)+0.5D0*r0*(2**0.5)
126 if (mmpi_myid == 0) write(*,*) 'at ', i, latmin(i), latcenter(i), latmax(i)
127 end do
128
129 end subroutine reg_getlatitude
130
131
132 subroutine getr0(r0km, r0, nlatband)
133 !
134 !:Purpose: The user specifies a radius of r0km (in km) within which
135 ! stations can be taken together for simultaneous analysis
136 ! in a single batch. For simplicity of the search procedures
137 ! we - instead - use somewhat smaller squares with sides of
138 ! at most r0km*root(2). We do have circles at the two poles.
139 ! Here we compute the radius r0 in radians such that we arrive
140 ! exactly at nlatband latitude bands. The equator separates
141 ! two latitude bands.
142 !
143 implicit none
144
145 ! Arguments:
146 real(8), intent(in) :: r0km
147 real(8), intent(out) :: r0
148 integer, intent(out) :: nlatband
149
150 ! Locals:
151 real(8) :: dpie, eps, difference, r0new, dminpole, nbandr
152 integer :: nbandi
153
154 dpie = MPC_PI_R8
155
156 ! figure out if r0km*root(2) is a divisor of 10000-r0km
157
158 eps = 0.01
159 ! we put a circle with radius r0km around each pole.
160 dminpole = 10000.-r0km
161 ! potential fractional number of latitude bands.
162 nbandr = dminpole/(r0km*(2**0.5))
163 nbandi = nint(nbandr)
164 difference = abs(dble(nbandi)-nbandr)
165 if (difference.gt.eps) then
166 ! increasing number of bands to nbandi
167 nbandi = nbandi+1
168 end if
169 ! decrease r0km to r0new
170 r0new = 10000./(dble(nbandi)*2**0.5 + 1.)
171 ! compute the corresponding r0 in radians
172 r0 = (r0new/20000.)*dpie
173 ! two hemispheres with a pole
174 nlatband = 2*(nbandi+1)
175
176 end subroutine getr0
177
178 subroutine reg_init_struct(lsc, r0_km, r1_km, rz_logp)
179 !
180 !:Purpose: store and derive parameters related to the localization
181 ! in a structure.
182 !
183 implicit none
184
185 ! Arguments:
186 type(struct_reg), intent(out) :: lsc
187 real(8), intent(in) :: r0_km
188 real(8), intent(in) :: r1_km
189 real(8), intent(in) :: rz_logp
190
191 ! Locals:
192 real(8) :: dpie
193 integer :: ione
194
195 lsc%r0_km = r0_km
196 lsc%r1_km = r1_km
197 lsc%rz_logp = rz_logp
198
199 call getr0(lsc%r0_km, lsc%r0_rad, lsc%nlatband)
200 lsc%r1_rad = lsc%r1_km/6370.
201 lsc%scale_hor = 0.5*lsc%r1_rad
202 lsc%scale_ver = 0.5*lsc%rz_logp
203
204 dpie = MPC_PI_R8
205 lsc%mxreg = floor(dpie/(lsc%r1_rad+(2**0.5)*lsc%r0_rad)**2)
206
207 ione = 1
208 lsc%mxreg = max(lsc%mxreg, ione)
209
210 end subroutine reg_init_struct
211
212
213 subroutine reg_locatestn(r0, lat, lon, nlatband, nlonblock, &
214 nblockoffset, iblock)
215 !
216 !:Purpose: locate the lat-lon block in which position (lat,lon) falls
217 !
218 !:Arguments:
219 ! r0: radius of a region (in radians)
220 ! (lat,lon): latitude and longitude in radians
221 ! nlatband: number of latitude bands
222 ! nlonblock: number of longitudinal blocks at each latitude band
223 ! nblockoffset: offset (the number of blocks south of each
224 ! latitude band)
225 ! iblock: block number
226 !
227 implicit none
228
229 ! Arguments:
230 integer, intent(out) :: iblock
231 integer, intent(in) :: nlatband
232 integer, intent(in) :: nlonblock(nlatband)
233 integer, intent(in) :: nblockoffset(nlatband)
234 real(4), intent(inout) :: lat
235 real(4), intent(in) :: lon
236 real(8), intent(in) :: r0
237
238 ! Locals:
239 integer :: ilat, ilon
240 real(8) :: dpie
241
242 dpie = MPC_PI_R8
243
244 if (lat.le.(-0.5D0*dpie+r0)) then
245 ilat = 1
246 else if (lat.gt.(0.5D0*dpie-r0)) then
247 ilat = nlatband
248 else
249 lat = lat+0.5D0*dpie-r0
250 ilat = 1+ceiling(lat/(r0*(2**0.5)))
251 end if
252 ilat = min(max(ilat,1),nlatband)
253 ilon = ceiling((lon/(2.0D0*dpie))*nlonblock(ilat))
254 ! added for points at the Greenwich meridian
255 ilon = min(max(ilon,1),nlonblock(ilat))
256 iblock = nblockoffset(ilat)+ilon
257
258 end subroutine reg_locatestn
259
260end module regions_mod
