| blob | 113208a55721faf683fbdf2ef8097caa575eaaac |
1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 ! Module: process_tile
3 !
4 ! Description:
5 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
6 module process_tile_module
8 use module_debug
11 contains
14 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
15 ! Name: process_tile
16 !
17 ! Purpose: To process a tile, whose lower-left corner is at
18 ! (tile_i_min, tile_j_min) and whose upper-right corner is at
19 ! (tile_i_max, tile_j_max), of the model grid given by which_domain
20 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
21 subroutine process_tile(which_domain, grid_type, dynopt, &
22 dummy_start_dom_i, dummy_end_dom_i, &
23 dummy_start_dom_j, dummy_end_dom_j, &
24 dummy_start_patch_i, dummy_end_patch_i, &
25 dummy_start_patch_j, dummy_end_patch_j, &
26 extra_col, extra_row)
28 use bitarray_module
29 use hash_module
30 use llxy_module
31 use misc_definitions_module
32 use output_module
33 use smooth_module
34 use source_data_module
36 implicit none
38 ! Arguments
39 integer, intent(in) :: which_domain, dynopt, &
40 dummy_start_dom_i, dummy_end_dom_i, dummy_start_dom_j, dummy_end_dom_j, &
41 dummy_start_patch_i, dummy_end_patch_i, dummy_start_patch_j, dummy_end_patch_j
42 logical, intent(in) :: extra_col, extra_row
43 character (len=1), intent(in) :: grid_type
45 ! Local variables
46 integer :: i, j, k, kk, istatus, ifieldstatus, idomcatstatus, field_count
47 integer :: min_category, max_category, min_level, max_level, &
48 smth_opt, smth_passes, num_landmask_categories
49 integer :: start_dom_i, end_dom_i, start_dom_j, end_dom_j, end_dom_stag_i, end_dom_stag_j
50 integer :: start_patch_i, end_patch_i, start_patch_j, end_patch_j, end_patch_stag_i, end_patch_stag_j
51 integer :: start_mem_i, end_mem_i, start_mem_j, end_mem_j, end_mem_stag_i, end_mem_stag_j
52 integer :: sm1, em1, sm2, em2
53 integer :: istagger
54 integer, dimension(MAX_LANDMASK_CATEGORIES) :: landmask_value
55 real :: sum, dominant, msg_fill_val, topo_flag_val, mass_flag, land_total, water_total
56 real, dimension(16) :: corner_lats, corner_lons
57 real, pointer, dimension(:,:) :: xlat_array, xlon_array, &
58 xlat_array_u, xlon_array_u, &
59 xlat_array_v, xlon_array_v, &
60 xlat_array_corner, xlon_array_corner, &
61 clat_array, clon_array, &
62 xlat_array_subgrid, xlon_array_subgrid, &
63 f_array, e_array, &
64 mapfac_array_m_x, mapfac_array_u_x, mapfac_array_v_x, &
65 mapfac_array_m_y, mapfac_array_u_y, mapfac_array_v_y, &
66 mapfac_array_x_subgrid, mapfac_array_y_subgrid, &
67 sina_array, cosa_array
68 real, pointer, dimension(:,:) :: xlat_ptr, xlon_ptr, mapfac_ptr_x, mapfac_ptr_y, landmask, dominant_field
69 real, pointer, dimension(:,:,:) :: field, slp_field
70 logical :: is_water_mask, only_save_dominant, halt_on_missing
71 character (len=19) :: datestr
72 character (len=128) :: fieldname, gradname, domname, landmask_name
73 character (len=256) :: temp_string
74 type (bitarray) :: processed_domain
75 type (hashtable) :: processed_fieldnames
76 character (len=128), dimension(2) :: dimnames
77 integer :: sub_x, sub_y
78 integer :: opt_status
80 ! Probably not all of these nullify statements are needed...
81 nullify(xlat_array)
82 nullify(xlon_array)
83 nullify(xlat_array_u)
84 nullify(xlon_array_u)
85 nullify(xlat_array_v)
86 nullify(xlon_array_v)
87 nullify(xlat_array_corner)
88 nullify(xlon_array_corner)
89 nullify(clat_array)
90 nullify(clon_array)
91 nullify(xlat_array_subgrid)
92 nullify(xlon_array_subgrid)
93 nullify(f_array)
94 nullify(e_array)
95 nullify(mapfac_array_m_x)
96 nullify(mapfac_array_u_x)
97 nullify(mapfac_array_v_x)
98 nullify(mapfac_array_m_y)
99 nullify(mapfac_array_u_y)
100 nullify(mapfac_array_v_y)
101 nullify(mapfac_array_x_subgrid)
102 nullify(mapfac_array_y_subgrid)
103 nullify(sina_array)
104 nullify(cosa_array)
105 nullify(xlat_ptr)
106 nullify(xlon_ptr)
107 nullify(mapfac_ptr_x)
108 nullify(mapfac_ptr_y)
109 nullify(landmask)
110 nullify(dominant_field)
111 nullify(field)
112 nullify(slp_field)
114 datestr = '0000-00-00_00:00:00'
115 field_count = 0
116 mass_flag=1.0
118 ! The following pertains primarily to the C grid
119 ! Determine whether only (n-1)th rows/columns should be computed for variables
120 ! on staggered grid. In a distributed memory situation, not every tile should
121 ! have only (n-1)th rows/columns computed, or we end up with (n-k)
122 ! rows/columns when there are k patches in the y/x direction
123 if (extra_col) then
124 start_patch_i = dummy_start_patch_i ! The seemingly pointless renaming of start
125 end_patch_i = dummy_end_patch_i - 1 ! naming convention with modified end_patch variables,
126 end_patch_stag_i = dummy_end_patch_i ! variables is so that we can maintain consistent
127 ! which are marked as intent(in)
128 start_mem_i = start_patch_i - HALO_WIDTH
129 end_mem_i = end_patch_i + HALO_WIDTH
130 end_mem_stag_i = end_patch_stag_i + HALO_WIDTH
131 else
132 start_patch_i = dummy_start_patch_i
133 end_patch_i = dummy_end_patch_i
134 end_patch_stag_i = dummy_end_patch_i
136 start_mem_i = start_patch_i - HALO_WIDTH
137 end_mem_i = end_patch_i + HALO_WIDTH
138 end_mem_stag_i = end_patch_stag_i + HALO_WIDTH
139 end if
141 if (extra_row) then
142 start_patch_j = dummy_start_patch_j
143 end_patch_j = dummy_end_patch_j - 1
144 end_patch_stag_j = dummy_end_patch_j
146 start_mem_j = start_patch_j - HALO_WIDTH
147 end_mem_j = end_patch_j + HALO_WIDTH
148 end_mem_stag_j = end_patch_stag_j + HALO_WIDTH
149 else
150 start_patch_j = dummy_start_patch_j
151 end_patch_j = dummy_end_patch_j
152 end_patch_stag_j = dummy_end_patch_j
154 start_mem_j = start_patch_j - HALO_WIDTH
155 end_mem_j = end_patch_j + HALO_WIDTH
156 end_mem_stag_j = end_patch_stag_j + HALO_WIDTH
157 end if
159 start_dom_i = dummy_start_dom_i
160 if (grid_type == 'C') then
161 end_dom_i = dummy_end_dom_i - 1
162 end_dom_stag_i = dummy_end_dom_i
163 else if (grid_type == 'E') then
164 end_dom_i = dummy_end_dom_i
165 end_dom_stag_i = dummy_end_dom_i
166 end if
168 start_dom_j = dummy_start_dom_j
169 if (grid_type == 'C') then
170 end_dom_j = dummy_end_dom_j - 1
171 end_dom_stag_j = dummy_end_dom_j
172 else if (grid_type == 'E') then
173 end_dom_j = dummy_end_dom_j
174 end_dom_stag_j = dummy_end_dom_j
175 end if
177 ! Allocate arrays to hold all lat/lon fields; these will persist for the duration of
178 ! the process_tile routine
179 ! For C grid, we have M, U, and V points
180 ! For E grid, we have only M and V points
181 allocate(xlat_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
182 allocate(xlon_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
183 allocate(xlat_array_v(start_mem_i:end_mem_i, start_mem_j:end_mem_stag_j))
184 allocate(xlon_array_v(start_mem_i:end_mem_i, start_mem_j:end_mem_stag_j))
185 if (grid_type == 'C') then
186 allocate(xlat_array_u(start_mem_i:end_mem_stag_i, start_mem_j:end_mem_j))
187 allocate(xlon_array_u(start_mem_i:end_mem_stag_i, start_mem_j:end_mem_j))
188 allocate(clat_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
189 allocate(clon_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
190 allocate(xlat_array_corner(start_mem_i:end_mem_stag_i, start_mem_j:end_mem_stag_j))
191 allocate(xlon_array_corner(start_mem_i:end_mem_stag_i, start_mem_j:end_mem_stag_j))
192 end if
193 nullify(xlat_array_subgrid)
194 nullify(xlon_array_subgrid)
195 nullify(mapfac_array_x_subgrid)
196 nullify(mapfac_array_y_subgrid)
198 ! Initialize hash table to track which fields have been processed
199 call hash_init(processed_fieldnames)
201 !
202 ! Calculate lat/lon for every point in the tile (XLAT and XLON)
203 ! The xlat_array and xlon_array arrays will be used in processing other fields
204 !
205 call mprintf(.true.,STDOUT,' Processing XLAT and XLONG')
207 if (grid_type == 'C') then
208 call get_lat_lon_fields(xlat_array, xlon_array, start_mem_i, &
209 start_mem_j, end_mem_i, end_mem_j, M)
210 call get_lat_lon_fields(xlat_array_v, xlon_array_v, start_mem_i, &
211 start_mem_j, end_mem_i, end_mem_stag_j, V)
212 call get_lat_lon_fields(xlat_array_u, xlon_array_u, start_mem_i, &
213 start_mem_j, end_mem_stag_i, end_mem_j, U)
214 call get_lat_lon_fields(xlat_array_corner, xlon_array_corner, start_mem_i, &
215 start_mem_j, end_mem_stag_i, end_mem_stag_j, CORNER)
216 call get_lat_lon_fields(clat_array, clon_array, start_mem_i, &
217 start_mem_j, end_mem_i, end_mem_j, M, comp_ll=.true.)
219 corner_lats(1) = xlat_array(start_patch_i,start_patch_j)
220 corner_lats(2) = xlat_array(start_patch_i,end_patch_j)
221 corner_lats(3) = xlat_array(end_patch_i,end_patch_j)
222 corner_lats(4) = xlat_array(end_patch_i,start_patch_j)
224 corner_lats(5) = xlat_array_u(start_patch_i,start_patch_j)
225 corner_lats(6) = xlat_array_u(start_patch_i,end_patch_j)
226 corner_lats(7) = xlat_array_u(end_patch_stag_i,end_patch_j)
227 corner_lats(8) = xlat_array_u(end_patch_stag_i,start_patch_j)
229 corner_lats(9) = xlat_array_v(start_patch_i,start_patch_j)
230 corner_lats(10) = xlat_array_v(start_patch_i,end_patch_stag_j)
231 corner_lats(11) = xlat_array_v(end_patch_i,end_patch_stag_j)
232 corner_lats(12) = xlat_array_v(end_patch_i,start_patch_j)
234 call xytoll(real(start_patch_i)-0.5, real(start_patch_j)-0.5, corner_lats(13), corner_lons(13), M)
235 call xytoll(real(start_patch_i)-0.5, real(end_patch_j)+0.5, corner_lats(14), corner_lons(14), M)
236 call xytoll(real(end_patch_i)+0.5, real(end_patch_j)+0.5, corner_lats(15), corner_lons(15), M)
237 call xytoll(real(end_patch_i)+0.5, real(start_patch_j)-0.5, corner_lats(16), corner_lons(16), M)
239 corner_lons(1) = xlon_array(start_patch_i,start_patch_j)
240 corner_lons(2) = xlon_array(start_patch_i,end_patch_j)
241 corner_lons(3) = xlon_array(end_patch_i,end_patch_j)
242 corner_lons(4) = xlon_array(end_patch_i,start_patch_j)
244 corner_lons(5) = xlon_array_u(start_patch_i,start_patch_j)
245 corner_lons(6) = xlon_array_u(start_patch_i,end_patch_j)
246 corner_lons(7) = xlon_array_u(end_patch_stag_i,end_patch_j)
247 corner_lons(8) = xlon_array_u(end_patch_stag_i,start_patch_j)
249 corner_lons(9) = xlon_array_v(start_patch_i,start_patch_j)
250 corner_lons(10) = xlon_array_v(start_patch_i,end_patch_stag_j)
251 corner_lons(11) = xlon_array_v(end_patch_i,end_patch_stag_j)
252 corner_lons(12) = xlon_array_v(end_patch_i,start_patch_j)
254 else if (grid_type == 'E') then
255 call get_lat_lon_fields(xlat_array, xlon_array, start_mem_i, &
256 start_mem_j, end_mem_i, end_mem_j, HH)
257 call get_lat_lon_fields(xlat_array_v, xlon_array_v, start_mem_i, &
258 start_mem_j, end_mem_i, end_mem_stag_j, VV)
260 corner_lats(1) = xlat_array(start_patch_i,start_patch_j)
261 corner_lats(2) = xlat_array(start_patch_i,end_patch_j)
262 corner_lats(3) = xlat_array(end_patch_i,end_patch_j)
263 corner_lats(4) = xlat_array(end_patch_i,start_patch_j)
265 corner_lats(5) = xlat_array_v(start_patch_i,start_patch_j)
266 corner_lats(6) = xlat_array_v(start_patch_i,end_patch_stag_j)
267 corner_lats(7) = xlat_array_v(end_patch_i,end_patch_stag_j)
268 corner_lats(8) = xlat_array_v(end_patch_i,start_patch_j)
270 corner_lats(9) = 0.0
271 corner_lats(10) = 0.0
272 corner_lats(11) = 0.0
273 corner_lats(12) = 0.0
275 corner_lats(13) = 0.0
276 corner_lats(14) = 0.0
277 corner_lats(15) = 0.0
278 corner_lats(16) = 0.0
280 corner_lons(1) = xlon_array(start_patch_i,start_patch_j)
281 corner_lons(2) = xlon_array(start_patch_i,end_patch_j)
282 corner_lons(3) = xlon_array(end_patch_i,end_patch_j)
283 corner_lons(4) = xlon_array(end_patch_i,start_patch_j)
285 corner_lons(5) = xlon_array_v(start_patch_i,start_patch_j)
286 corner_lons(6) = xlon_array_v(start_patch_i,end_patch_stag_j)
287 corner_lons(7) = xlon_array_v(end_patch_i,end_patch_stag_j)
288 corner_lons(8) = xlon_array_v(end_patch_i,start_patch_j)
290 corner_lons(9) = 0.0
291 corner_lons(10) = 0.0
292 corner_lons(11) = 0.0
293 corner_lons(12) = 0.0
295 corner_lons(13) = 0.0
296 corner_lons(14) = 0.0
297 corner_lons(15) = 0.0
298 corner_lons(16) = 0.0
300 end if
302 ! Initialize the output module now that we have the corner point lats/lons
303 call output_init(which_domain, 'OUTPUT FROM GEOGRID V3.7.1', '0000-00-00_00:00:00', grid_type, dynopt, &
304 corner_lats, corner_lons, &
305 start_dom_i, end_dom_i, start_dom_j, end_dom_j, &
306 start_patch_i, end_patch_i, start_patch_j, end_patch_j, &
307 start_mem_i, end_mem_i, start_mem_j, end_mem_j, &
308 extra_col, extra_row)
310 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, &
311 'XLAT_M', datestr, real_array = xlat_array)
312 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, &
313 'XLONG_M', datestr, real_array = xlon_array)
314 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_stag_j, 1, 1, &
315 'XLAT_V', datestr, real_array = xlat_array_v)
316 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_stag_j, 1, 1, &
317 'XLONG_V', datestr, real_array = xlon_array_v)
318 if (grid_type == 'C') then
319 call write_field(start_mem_i, end_mem_stag_i, start_mem_j, end_mem_j, 1, 1, &
320 'XLAT_U', datestr, real_array = xlat_array_u)
321 call write_field(start_mem_i, end_mem_stag_i, start_mem_j, end_mem_j, 1, 1, &
322 'XLONG_U', datestr, real_array = xlon_array_u)
323 call write_field(start_mem_i, end_mem_stag_i, start_mem_j, end_mem_stag_j, 1, 1, &
324 'XLAT_C', datestr, real_array = xlat_array_corner)
325 call write_field(start_mem_i, end_mem_stag_i, start_mem_j, end_mem_stag_j, 1, 1, &
326 'XLONG_C', datestr, real_array = xlon_array_corner)
327 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, &
328 'CLAT', datestr, real_array = clat_array)
329 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, &
330 'CLONG', datestr, real_array = clon_array)
332 if (associated(clat_array)) deallocate(clat_array)
333 if (associated(clon_array)) deallocate(clon_array)
335 end if
338 !
339 ! Calculate map factor for current domain
340 !
341 if (grid_type == 'C') then
342 call mprintf(.true.,STDOUT,' Processing MAPFAC')
344 allocate(mapfac_array_m_x(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
345 allocate(mapfac_array_m_y(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
346 call get_map_factor(xlat_array, xlon_array, mapfac_array_m_x, mapfac_array_m_y, start_mem_i, &
347 start_mem_j, end_mem_i, end_mem_j)
348 ! Global WRF uses map scale factors in X and Y directions, but "regular" WRF uses a single MSF
349 ! on each staggering. In the case of regular WRF, we can assume that MAPFAC_MX = MAPFAC_MY = MAPFAC_M,
350 ! and so we can simply write MAPFAC_MX as the MAPFAC_M field. Ultimately, when global WRF is
351 ! merged into the WRF trunk, we will need only two map scale factor fields for each staggering,
352 ! in the x and y directions, and these will be the same in the case of non-Cassini projections
353 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'MAPFAC_M', &
354 datestr, real_array = mapfac_array_m_x)
355 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'MAPFAC_MX', &
356 datestr, real_array = mapfac_array_m_x)
357 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'MAPFAC_MY', &
358 datestr, real_array = mapfac_array_m_y)
360 allocate(mapfac_array_v_x(start_mem_i:end_mem_i, start_mem_j:end_mem_stag_j))
361 allocate(mapfac_array_v_y(start_mem_i:end_mem_i, start_mem_j:end_mem_stag_j))
362 call get_map_factor(xlat_array_v, xlon_array_v, mapfac_array_v_x, mapfac_array_v_y, start_mem_i, &
363 start_mem_j, end_mem_i, end_mem_stag_j)
364 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_stag_j, 1, 1, 'MAPFAC_V', &
365 datestr, real_array = mapfac_array_v_x)
366 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_stag_j, 1, 1, 'MAPFAC_VX', &
367 datestr, real_array = mapfac_array_v_x)
368 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_stag_j, 1, 1, 'MAPFAC_VY', &
369 datestr, real_array = mapfac_array_v_y)
371 allocate(mapfac_array_u_x(start_mem_i:end_mem_stag_i, start_mem_j:end_mem_j))
372 allocate(mapfac_array_u_y(start_mem_i:end_mem_stag_i, start_mem_j:end_mem_j))
373 call get_map_factor(xlat_array_u, xlon_array_u, mapfac_array_u_x, mapfac_array_u_y, start_mem_i, &
374 start_mem_j, end_mem_stag_i, end_mem_j)
375 call write_field(start_mem_i, end_mem_stag_i, start_mem_j, end_mem_j, 1, 1, 'MAPFAC_U', &
376 datestr, real_array = mapfac_array_u_x)
377 call write_field(start_mem_i, end_mem_stag_i, start_mem_j, end_mem_j, 1, 1, 'MAPFAC_UX', &
378 datestr, real_array = mapfac_array_u_x)
379 call write_field(start_mem_i, end_mem_stag_i, start_mem_j, end_mem_j, 1, 1, 'MAPFAC_UY', &
380 datestr, real_array = mapfac_array_u_y)
382 end if
385 !
386 ! Coriolis parameters (E and F)
387 !
388 call mprintf(.true.,STDOUT,' Processing F and E')
390 allocate(f_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
391 allocate(e_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
393 call get_coriolis_parameters(xlat_array, f_array, e_array, &
394 start_mem_i, start_mem_j, end_mem_i, end_mem_j)
396 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'E', &
397 datestr, real_array = e_array)
398 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'F', &
399 datestr, real_array = f_array)
401 if (associated(f_array)) deallocate(f_array)
402 if (associated(e_array)) deallocate(e_array)
405 !
406 ! Rotation angle (SINALPHA and COSALPHA)
407 !
408 if (grid_type == 'C') then
409 call mprintf(.true.,STDOUT,' Processing ROTANG')
411 allocate(sina_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
412 allocate(cosa_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
414 call get_rotang(xlat_array, xlon_array, cosa_array, sina_array, &
415 start_mem_i, start_mem_j, end_mem_i, end_mem_j)
417 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'SINALPHA', &
418 datestr, real_array = sina_array)
419 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'COSALPHA', &
420 datestr, real_array = cosa_array)
422 if (associated(sina_array)) deallocate(sina_array)
423 if (associated(cosa_array)) deallocate(cosa_array)
424 end if
426 ! Every field up until now should probably just be processed regardless of what the user
427 ! has specified for fields to be processed.
428 ! Hereafter, we process user-specified fields
430 !
431 ! First process the field that we will derive a landmask from
432 !
433 call get_landmask_field(geog_data_res(which_domain), landmask_name, is_water_mask, landmask_value, istatus)
435 do kk=1,MAX_LANDMASK_CATEGORIES
436 if (landmask_value(kk) == INVALID) then
437 num_landmask_categories = kk-1
438 exit
439 end if
440 end do
441 if (kk > MAX_LANDMASK_CATEGORIES) num_landmask_categories = MAX_LANDMASK_CATEGORIES
443 if (istatus /= 0) then
444 call mprintf(.true.,WARN,'No field specified for landmask calculation. Will set landmask=1 at every grid point.')
446 allocate(landmask(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
447 landmask = 1.
448 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'LANDMASK', &
449 datestr, landmask)
451 else
453 allocate(landmask(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
454 landmask = 1.
456 call mprintf(.true.,STDOUT,' Processing %s', s1=trim(landmask_name))
458 call get_missing_fill_value(landmask_name, msg_fill_val, istatus)
459 if (istatus /= 0) msg_fill_val = NAN
461 call get_halt_on_missing(landmask_name, halt_on_missing, istatus)
462 if (istatus /= 0) halt_on_missing = .false.
464 ! Do we calculate a dominant category for this field?
465 call get_domcategory_name(landmask_name, domname, only_save_dominant, idomcatstatus)
467 temp_string = ' '
468 temp_string(1:128) = landmask_name
469 call hash_insert(processed_fieldnames, temp_string)
471 call get_max_categories(landmask_name, min_category, max_category, istatus)
472 allocate(field(start_mem_i:end_mem_i, start_mem_j:end_mem_j, min_category:max_category))
474 if (.not. only_save_dominant) then
475 field_count = field_count + 1
476 call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
477 i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=landmask_name)
478 else
479 field_count = field_count + 1
480 call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
481 i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=domname)
482 end if
484 if (grid_type == 'C') then
485 call calc_field(landmask_name, field, xlat_array, xlon_array, M, &
486 start_mem_i, end_mem_i, start_mem_j, end_mem_j, &
487 min_category, max_category, processed_domain, 1, landmask=landmask, sr_x=1, sr_y=1)
488 else if (grid_type == 'E') then
489 call calc_field(landmask_name, field, xlat_array, xlon_array, HH, &
490 start_mem_i, end_mem_i, start_mem_j, end_mem_j, &
491 min_category, max_category, processed_domain, 1, landmask=landmask, sr_x=1, sr_y=1)
492 end if
494 ! Find fractions
495 do i=start_mem_i, end_mem_i
496 do j=start_mem_j, end_mem_j
497 sum = 0.
498 do k=min_category,max_category
499 sum = sum + field(i,j,k)
500 end do
501 if (sum > 0.0) then
502 do k=min_category,max_category
503 field(i,j,k) = field(i,j,k) / sum
504 end do
505 else
506 do k=min_category,max_category
507 field(i,j,k) = msg_fill_val
508 end do
509 end if
510 end do
511 end do
513 ! If user wants to halt when a missing value is found in output field, check now
514 if (halt_on_missing) then
515 do i=start_mem_i, end_mem_i
516 do j=start_mem_j, end_mem_j
517 ! Only need to examine k=1
518 if (field(i,j,1) == msg_fill_val) then
519 call mprintf(.true.,ERROR,' Missing value encountered in output field. Quitting.')
520 end if
521 end do
522 end do
523 end if
525 if (is_water_mask) then
526 call mprintf(.true.,STDOUT,' Calculating landmask from %s ( WATER =', &
527 newline=.false.,s1=trim(landmask_name))
528 else
529 call mprintf(.true.,STDOUT,' Calculating landmask from %s ( LAND =', &
530 newline=.false.,s1=trim(landmask_name))
531 end if
532 do k = 1, num_landmask_categories
533 call mprintf(.true.,STDOUT,' %i',newline=.false.,i1=landmask_value(k))
534 if (k == num_landmask_categories) call mprintf(.true.,STDOUT,')')
535 end do
537 ! Calculate landmask
538 if (is_water_mask) then
539 do i=start_mem_i, end_mem_i
540 do j=start_mem_j, end_mem_j
541 water_total = -1.
542 do k=1,num_landmask_categories
543 if (landmask_value(k) >= min_category .and. landmask_value(k) <= max_category) then
544 if (field(i,j,landmask_value(k)) /= msg_fill_val) then
545 if (water_total < 0.) water_total = 0.
546 water_total = water_total + field(i,j,landmask_value(k))
547 else
548 water_total = -1.
549 exit
550 end if
551 end if
552 end do
553 if (water_total >= 0.0) then
554 if (water_total < 0.50) then
555 landmask(i,j) = 1.
556 else
557 landmask(i,j) = 0.
558 end if
559 else
560 landmask(i,j) = -1.
561 end if
562 end do
563 end do
564 else
565 do i=start_mem_i, end_mem_i
566 do j=start_mem_j, end_mem_j
567 land_total = -1.
568 do k=1,num_landmask_categories
569 if (landmask_value(k) >= min_category .and. landmask_value(k) <= max_category) then
570 if (field(i,j,landmask_value(k)) /= msg_fill_val) then
571 if (land_total < 0.) land_total = 0.
572 land_total = land_total + field(i,j,landmask_value(k))
573 else
574 land_total = -1.
575 exit
576 end if
577 end if
578 end do
579 if (land_total >= 0.0) then
580 if (land_total > 0.50) then
581 landmask(i,j) = 1.
582 else
583 landmask(i,j) = 0.
584 end if
585 else
586 landmask(i,j) = -1.
587 end if
588 end do
589 end do
590 end if
592 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'LANDMASK', &
593 datestr, landmask)
595 ! If we should only save the dominant category, then no need to write out fractional field
596 if (.not.only_save_dominant .or. (idomcatstatus /= 0)) then
598 ! Finally, we may be asked to smooth the fractional field
599 call get_smooth_option(landmask_name, smth_opt, smth_passes, istatus)
600 if (istatus == 0) then
602 if (grid_type == 'C') then
603 if (smth_opt == ONETWOONE) then
604 call one_two_one(field, &
605 start_patch_i, end_patch_i, &
606 start_patch_j, end_patch_j, &
607 start_mem_i, end_mem_i, &
608 start_mem_j, end_mem_j, &
609 min_category, max_category, &
610 smth_passes, msg_fill_val)
611 else if (smth_opt == SMTHDESMTH) then
612 call smth_desmth(field, &
613 start_patch_i, end_patch_i, &
614 start_patch_j, end_patch_j, &
615 start_mem_i, end_mem_i, &
616 start_mem_j, end_mem_j, &
617 min_category, max_category, &
618 smth_passes, msg_fill_val)
619 else if (smth_opt == SMTHDESMTH_SPECIAL) then
620 call smth_desmth_special(field, &
621 start_patch_i, end_patch_i, &
622 start_patch_j, end_patch_j, &
623 start_mem_i, end_mem_i, &
624 start_mem_j, end_mem_j, &
625 min_category, max_category, &
626 smth_passes, msg_fill_val)
627 end if
628 else if (grid_type == 'E') then
629 if (smth_opt == ONETWOONE) then
630 call one_two_one_egrid(field, &
631 start_patch_i, end_patch_i, &
632 start_patch_j, end_patch_j, &
633 start_mem_i, end_mem_i, &
634 start_mem_j, end_mem_j, &
635 min_category, max_category, &
636 smth_passes, msg_fill_val, 1.0)
637 else if (smth_opt == SMTHDESMTH) then
638 call smth_desmth_egrid(field, &
639 start_patch_i, end_patch_i, &
640 start_patch_j, end_patch_j, &
641 start_mem_i, end_mem_i, &
642 start_mem_j, end_mem_j, &
643 min_category, max_category, &
644 smth_passes, msg_fill_val, 1.0)
645 else if (smth_opt == SMTHDESMTH_SPECIAL) then
646 call mprintf(.true.,WARN,'smth-desmth_special is not currently implemented for NMM. '// &
647 'No smoothing will be done.')
648 end if
649 end if
651 end if
653 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, &
654 min_category, max_category, trim(landmask_name), &
655 datestr, real_array=field)
656 end if
658 if (idomcatstatus == 0) then
659 allocate(dominant_field(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
661 if (.not. only_save_dominant) then
662 field_count = field_count + 1
663 call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
664 i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=domname)
665 end if
667 do i=start_mem_i, end_mem_i
668 do j=start_mem_j, end_mem_j
669 if ((landmask(i,j) == 1. .and. is_water_mask) .or. &
670 (landmask(i,j) == 0. .and. .not.is_water_mask)) then
671 dominant = 0.
672 dominant_field(i,j) = real(min_category-1)
673 do k=min_category,max_category
674 do kk=1,num_landmask_categories
675 if (k == landmask_value(kk)) exit
676 end do
677 if (field(i,j,k) > dominant .and. kk > num_landmask_categories) then
678 dominant_field(i,j) = real(k)
679 dominant = field(i,j,k)
680 end if
681 end do
682 else
683 dominant = 0.
684 dominant_field(i,j) = real(min_category-1)
685 do k=min_category,max_category
686 do kk=1,num_landmask_categories
687 if (field(i,j,k) > dominant .and. k == landmask_value(kk)) then
688 dominant_field(i,j) = real(k)
689 dominant = field(i,j,k)
690 end if
691 end do
692 end do
693 end if
694 end do
695 end do
697 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, trim(domname), &
698 datestr, dominant_field)
700 deallocate(dominant_field)
701 end if
703 deallocate(field)
704 end if
706 !
707 ! Now process all other fields specified by the user
708 !
709 call reset_next_field()
710 ifieldstatus = 0
711 do while (ifieldstatus == 0)
712 call get_next_fieldname(fieldname, ifieldstatus)
714 ! There is another field in the GEOGRID.TBL file
715 if (ifieldstatus == 0) then
716 temp_string(1:128) = fieldname
718 call get_source_opt_status(fieldname, 0, opt_status)
720 ! If this field is still to be processed
721 if (.not. hash_search(processed_fieldnames, temp_string) .and. opt_status == 0) then
723 call hash_insert(processed_fieldnames, temp_string)
724 call mprintf(.true.,STDOUT,' Processing %s', s1=trim(fieldname))
726 call get_output_stagger(fieldname, istagger, istatus)
727 dimnames(:) = 'null'
728 call get_subgrid_dim_name(which_domain, fieldname, dimnames, &
729 sub_x, sub_y, istatus)
731 if (istagger == M .or. (sub_x > 1) .or. (sub_y > 1)) then
732 sm1 = start_mem_i
733 em1 = end_mem_i
734 sm2 = start_mem_j
735 em2 = end_mem_j
736 xlat_ptr => xlat_array
737 xlon_ptr => xlon_array
738 mapfac_ptr_x => mapfac_array_m_x
739 mapfac_ptr_y => mapfac_array_m_y
740 else if (istagger == U) then ! In the case that extra_cols = .false.
741 sm1 = start_mem_i ! we should have that end_mem_stag is
742 em1 = end_mem_stag_i ! the same as end_mem, so we do not need
743 sm2 = start_mem_j ! to check extra_cols or extra rows here
744 em2 = end_mem_j
745 xlat_ptr => xlat_array_u
746 xlon_ptr => xlon_array_u
747 mapfac_ptr_x => mapfac_array_u_x
748 mapfac_ptr_y => mapfac_array_u_y
749 else if (istagger == V) then
750 sm1 = start_mem_i
751 em1 = end_mem_i
752 sm2 = start_mem_j
753 em2 = end_mem_stag_j
754 xlat_ptr => xlat_array_v
755 xlon_ptr => xlon_array_v
756 mapfac_ptr_x => mapfac_array_v_x
757 mapfac_ptr_y => mapfac_array_v_y
758 else if (istagger == HH) then ! E grid
759 sm1 = start_mem_i
760 em1 = end_mem_i
761 sm2 = start_mem_j
762 em2 = end_mem_j
763 xlat_ptr => xlat_array
764 xlon_ptr => xlon_array
765 mapfac_ptr_x => mapfac_array_m_x
766 mapfac_ptr_y => mapfac_array_m_y
767 else if (istagger == VV) then ! E grid
768 sm1 = start_mem_i
769 em1 = end_mem_i
770 sm2 = start_mem_j
771 em2 = end_mem_stag_j
772 xlat_ptr => xlat_array_v
773 xlon_ptr => xlon_array_v
774 mapfac_ptr_x => mapfac_array_v_x
775 mapfac_ptr_y => mapfac_array_v_y
776 end if
778 if (sub_x > 1) then
779 sm1 = (start_mem_i - 1)*sub_x + 1
780 if (extra_col) then
781 em1 = (end_mem_i + 1)*sub_x
782 else
783 em1 = (end_mem_i )*sub_x
784 end if
785 end if
786 if (sub_y > 1)then
787 sm2 = (start_mem_j - 1)*sub_y + 1
788 if (extra_row) then
789 em2 = (end_mem_j + 1)*sub_y
790 else
791 em2 = (end_mem_j )*sub_y
792 end if
793 end if
795 !BUG: This should probably be moved up to where other lat/lon fields are calculated, and we should
796 ! just determine whether we will have any subgrids or not at that point
797 if ((sub_x > 1) .or. (sub_y > 1)) then
798 ! if (associated(xlat_array_subgrid)) deallocate(xlat_array_subgrid)
799 ! if (associated(xlon_array_subgrid)) deallocate(xlon_array_subgrid)
800 ! if (associated(mapfac_array_x_subgrid)) deallocate(mapfac_array_x_subgrid)
801 ! if (associated(mapfac_array_y_subgrid)) deallocate(mapfac_array_y_subgrid)
802 allocate(xlat_array_subgrid(sm1:em1,sm2:em2))
803 allocate(xlon_array_subgrid(sm1:em1,sm2:em2))
804 allocate(mapfac_array_x_subgrid(sm1:em1,sm2:em2))
805 allocate(mapfac_array_y_subgrid(sm1:em1,sm2:em2))
806 call get_lat_lon_fields(xlat_array_subgrid, xlon_array_subgrid, &
807 sm1, sm2, em1, em2, M, sub_x=sub_x, sub_y=sub_y)
808 xlat_ptr => xlat_array_subgrid
809 xlon_ptr => xlon_array_subgrid
810 call get_map_factor(xlat_ptr, xlon_ptr, mapfac_array_x_subgrid, &
811 mapfac_array_y_subgrid, sm1, sm2, em1, em2)
812 mapfac_ptr_x => mapfac_array_x_subgrid
813 mapfac_ptr_y => mapfac_array_y_subgrid
814 end if
816 call get_missing_fill_value(fieldname, msg_fill_val, istatus)
817 if (istatus /= 0) msg_fill_val = NAN
819 call get_halt_on_missing(fieldname, halt_on_missing, istatus)
820 if (istatus /= 0) halt_on_missing = .false.
822 ! Destination field type is CONTINUOUS
823 if (iget_fieldtype(fieldname,istatus) == CONTINUOUS) then
824 call get_max_levels(fieldname, min_level, max_level, istatus)
825 allocate(field(sm1:em1, sm2:em2, min_level:max_level))
827 field_count = field_count + 1
828 call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
829 i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=fieldname)
831 if ((sub_x > 1) .or. (sub_y > 1)) then
832 call calc_field(fieldname, field, xlat_ptr, xlon_ptr, istagger, &
833 sm1, em1, sm2, em2, min_level, max_level, &
834 processed_domain, 1, sr_x=sub_x, sr_y=sub_y)
835 else
836 call calc_field(fieldname, field, xlat_ptr, xlon_ptr, istagger, &
837 sm1, em1, sm2, em2, min_level, max_level, &
838 processed_domain, 1, landmask=landmask, sr_x=sub_x, sr_y=sub_y)
839 end if
841 ! If user wants to halt when a missing value is found in output field, check now
842 if (halt_on_missing) then
843 do i=sm1, em1
844 do j=sm2, em2
845 ! Only need to examine k=1
846 if (field(i,j,1) == msg_fill_val) then
847 call mprintf(.true.,ERROR,' Missing value encountered in output field. Quitting.')
848 end if
849 end do
850 end do
851 end if
853 ! We may be asked to smooth the fractional field
854 call get_smooth_option(fieldname, smth_opt, smth_passes, istatus)
855 if (istatus == 0) then
857 if (grid_type == 'C') then
858 if (smth_opt == ONETWOONE) then
859 call one_two_one(field, &
860 start_patch_i, end_patch_i, &
861 start_patch_j, end_patch_j, &
862 sm1, em1, &
863 sm2, em2, &
864 min_level, max_level, &
865 smth_passes, msg_fill_val)
866 else if (smth_opt == SMTHDESMTH) then
867 call smth_desmth(field, &
868 start_patch_i, end_patch_i, &
869 start_patch_j, end_patch_j, &
870 sm1, em1, &
871 sm2, em2, &
872 min_level, max_level, &
873 smth_passes, msg_fill_val)
874 else if (smth_opt == SMTHDESMTH_SPECIAL) then
875 call smth_desmth_special(field, &
876 start_patch_i, end_patch_i, &
877 start_patch_j, end_patch_j, &
878 sm1, em1, &
879 sm2, em2, &
880 min_level, max_level, &
881 smth_passes, msg_fill_val)
882 end if
884 else if (grid_type == 'E') then
886 if (trim(fieldname) == 'HGT_M' ) then
887 topo_flag_val=1.0
888 mass_flag=1.0
889 else if (trim(fieldname) == 'HGT_V') then
890 topo_flag_val=1.0
891 mass_flag=0.0
892 else
893 topo_flag_val=0.0
894 end if
896 if (smth_opt == ONETWOONE) then
897 call one_two_one_egrid(field, &
898 start_patch_i, end_patch_i, &
899 start_patch_j, end_patch_j, &
900 sm1, em1, &
901 sm2, em2, &
902 min_level, max_level, &
903 smth_passes, topo_flag_val, mass_flag)
904 else if (smth_opt == SMTHDESMTH) then
905 call smth_desmth_egrid(field, &
906 start_patch_i, end_patch_i, &
907 start_patch_j, end_patch_j, &
908 sm1, em1, &
909 sm2, em2, &
910 min_level, max_level, &
911 smth_passes, topo_flag_val, mass_flag)
912 else if (smth_opt == SMTHDESMTH_SPECIAL) then
913 call mprintf(.true.,WARN,'smth-desmth_special is not currently implemented for NMM. '// &
914 'No smoothing will be done.')
915 end if
917 end if
919 end if
921 call write_field(sm1, em1, sm2, em2, &
922 min_level, max_level, trim(fieldname), datestr, real_array=field)
924 ! Do we calculate directional derivatives from this field?
925 call get_dfdx_name(fieldname, gradname, istatus)
926 if (istatus == 0) then
927 allocate(slp_field(sm1:em1,sm2:em2,min_level:max_level))
929 field_count = field_count + 1
930 call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
931 i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=gradname)
933 if (grid_type == 'C') then
934 call calc_dfdx(field, slp_field, sm1, sm2, min_level, em1, em2, max_level, mapfac_ptr_x)
935 else if (grid_type == 'E') then
936 call calc_dfdx(field, slp_field, sm1, sm2, min_level, em1, em2, max_level)
937 end if
938 call write_field(sm1, em1, sm2, em2, &
939 min_level, max_level, trim(gradname), datestr, real_array=slp_field)
940 deallocate(slp_field)
941 end if
942 call get_dfdy_name(fieldname, gradname, istatus)
943 if (istatus == 0) then
944 allocate(slp_field(sm1:em1,sm2:em2,min_level:max_level))
946 field_count = field_count + 1
947 call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
948 i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=gradname)
950 if (grid_type == 'C') then
951 call calc_dfdy(field, slp_field, sm1, sm2, min_level, em1, em2, max_level, mapfac_ptr_y)
952 else if (grid_type == 'E') then
953 call calc_dfdy(field, slp_field, sm1, sm2, min_level, em1, em2, max_level)
954 end if
955 call write_field(sm1, em1, sm2, em2, &
956 min_level, max_level, trim(gradname), datestr, real_array=slp_field)
957 deallocate(slp_field)
958 end if
960 deallocate(field)
962 ! Destination field type is CATEGORICAL
963 else
964 call get_max_categories(fieldname, min_category, max_category, istatus)
965 allocate(field(sm1:em1, sm2:em2, min_category:max_category))
967 ! Do we calculate a dominant category for this field?
968 call get_domcategory_name(fieldname, domname, only_save_dominant, idomcatstatus)
970 if (.not. only_save_dominant) then
971 field_count = field_count + 1
972 call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
973 i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=fieldname)
974 else
975 field_count = field_count + 1
976 call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
977 i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=domname)
978 end if
980 if ((sub_x > 1) .or. (sub_y > 1)) then
981 call calc_field(fieldname, field, xlat_ptr, xlon_ptr, istagger, &
982 sm1, em1, sm2, em2, min_category, max_category, &
983 processed_domain, 1, sr_x=sub_x, sr_y=sub_y)
984 else
985 call calc_field(fieldname, field, xlat_ptr, xlon_ptr, istagger, &
986 sm1, em1, sm2, em2, min_category, max_category, &
987 processed_domain, 1, landmask=landmask, sr_x=sub_x, sr_y=sub_y)
988 end if
990 ! Find fractions
991 do i=sm1, em1
992 do j=sm2, em2
993 sum = 0.
994 do k=min_category,max_category
995 sum = sum + field(i,j,k)
996 end do
997 if (sum > 0.0) then
998 do k=min_category,max_category
999 field(i,j,k) = field(i,j,k) / sum
1000 end do
1001 else
1002 do k=min_category,max_category
1003 field(i,j,k) = msg_fill_val
1004 end do
1005 end if
1006 end do
1007 end do
1009 ! If user wants to halt when a missing value is found in output field, check now
1010 if (halt_on_missing) then
1011 do i=sm1, em1
1012 do j=sm2, em2
1013 ! Only need to examine k=1
1014 if (field(i,j,1) == msg_fill_val) then
1015 call mprintf(.true.,ERROR,' Missing value encountered in output field. Quitting.')
1016 end if
1017 end do
1018 end do
1019 end if
1021 ! If we should only save the dominant category, then no need to write out fractional field
1022 if (.not.only_save_dominant .or. (idomcatstatus /= 0)) then
1024 ! Finally, we may be asked to smooth the fractional field
1025 call get_smooth_option(fieldname, smth_opt, smth_passes, istatus)
1026 if (istatus == 0) then
1027 if (grid_type == 'C') then
1028 if (smth_opt == ONETWOONE) then
1029 call one_two_one(field, &
1030 start_patch_i, end_patch_i, &
1031 start_patch_j, end_patch_j, &
1032 sm1, em1, &
1033 sm2, em2, &
1034 min_category, max_category, &
1035 smth_passes, msg_fill_val)
1036 else if (smth_opt == SMTHDESMTH) then
1037 call smth_desmth(field, &
1038 start_patch_i, end_patch_i, &
1039 start_patch_j, end_patch_j, &
1040 sm1, em1, &
1041 sm2, em2, &
1042 min_category, max_category, &
1043 smth_passes, msg_fill_val)
1044 else if (smth_opt == SMTHDESMTH_SPECIAL) then
1045 call smth_desmth_special(field, &
1046 start_patch_i, end_patch_i, &
1047 start_patch_j, end_patch_j, &
1048 sm1, em1, &
1049 sm2, em2, &
1050 min_category, max_category, &
1051 smth_passes, msg_fill_val)
1052 end if
1053 else if (grid_type == 'E') then
1054 if (smth_opt == ONETWOONE) then
1055 call one_two_one_egrid(field, &
1056 start_patch_i, end_patch_i, &
1057 start_patch_j, end_patch_j, &
1058 sm1, em1, &
1059 sm2, em2, &
1060 min_category, max_category, &
1061 smth_passes, msg_fill_val, 1.0)
1062 else if (smth_opt == SMTHDESMTH) then
1063 call smth_desmth_egrid(field, &
1064 start_patch_i, end_patch_i, &
1065 start_patch_j, end_patch_j, &
1066 sm1, em1, &
1067 sm2, em2, &
1068 min_category, max_category, &
1069 smth_passes, msg_fill_val, 1.0)
1070 else if (smth_opt == SMTHDESMTH_SPECIAL) then
1071 call mprintf(.true.,WARN,'smth-desmth_special is not currently implemented for NMM. '// &
1072 'No smoothing will be done.')
1073 end if
1074 end if
1075 end if
1077 call write_field(sm1, em1, sm2, em2, &
1078 min_category, max_category, trim(fieldname), datestr, real_array=field)
1079 end if
1081 if (idomcatstatus == 0) then
1082 call mprintf(.true.,STDOUT,' Processing %s', s1=trim(domname))
1083 allocate(dominant_field(sm1:em1, sm2:em2))
1085 if (.not. only_save_dominant) then
1086 field_count = field_count + 1
1087 call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
1088 i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=domname)
1089 end if
1091 do i=sm1, em1
1092 do j=sm2, em2
1093 dominant = 0.
1094 dominant_field(i,j) = real(min_category-1)
1095 do k=min_category,max_category
1096 if (field(i,j,k) > dominant .and. field(i,j,k) /= msg_fill_val) then
1097 dominant_field(i,j) = real(k)
1098 dominant = field(i,j,k)
1099 ! else
1100 ! dominant_field(i,j) = nint(msg_fill_val)
1101 ! Maybe we should put an else clause here to set the category equal to the missing fill value?
1102 ! BUG: The problem here seems to be that, when we set a fraction equal to the missing fill value
1103 ! above, if the last fractional index we process here has been filled, we think that the dominant
1104 ! category should be set to the missing fill value. Perhaps we could do some check to only
1105 ! assign the msg_fill_val if no other valid category has been assigned? But this may still not
1106 ! work if the missing fill value is something like 0.5. Somehow use bitarrays, perhaps, to remember
1107 ! which points are missing and which just happen to have the missing fill value?
1108 end if
1109 end do
1110 if (dominant_field(i,j) == real(min_category-1)) dominant_field(i,j) = msg_fill_val
1111 end do
1112 end do
1113 call write_field(sm1, em1, sm2, em2, 1, 1, &
1114 trim(domname), datestr, dominant_field)
1115 deallocate(dominant_field)
1116 end if
1118 deallocate(field)
1120 if ((sub_x > 1) .or. (sub_y > 1)) then
1121 if (associated(xlat_array_subgrid)) deallocate(xlat_array_subgrid)
1122 if (associated(xlon_array_subgrid)) deallocate(xlon_array_subgrid)
1123 if (associated(mapfac_array_x_subgrid)) deallocate(mapfac_array_x_subgrid)
1124 if (associated(mapfac_array_y_subgrid)) deallocate(mapfac_array_y_subgrid)
1125 end if
1127 end if
1129 end if
1130 end if
1132 end do
1134 ! Close output
1135 call output_close()
1137 call hash_destroy(processed_fieldnames)
1139 ! Free up memory
1140 if (associated(xlat_array)) deallocate(xlat_array)
1141 if (associated(xlon_array)) deallocate(xlon_array)
1142 if (grid_type == 'C') then
1143 if (associated(xlat_array_u)) deallocate(xlat_array_u)
1144 if (associated(xlon_array_u)) deallocate(xlon_array_u)
1145 if (associated(xlat_array_corner)) deallocate(xlat_array_corner)
1146 if (associated(xlon_array_corner)) deallocate(xlon_array_corner)
1147 if (associated(mapfac_array_u_x)) deallocate(mapfac_array_u_x)
1148 if (associated(mapfac_array_u_y)) deallocate(mapfac_array_u_y)
1149 end if
1150 if (associated(xlat_array_v)) deallocate(xlat_array_v)
1151 if (associated(xlon_array_v)) deallocate(xlon_array_v)
1152 if (associated(mapfac_array_m_x)) deallocate(mapfac_array_m_x)
1153 if (associated(mapfac_array_m_y)) deallocate(mapfac_array_m_y)
1154 if (associated(mapfac_array_v_x)) deallocate(mapfac_array_v_x)
1155 if (associated(mapfac_array_v_y)) deallocate(mapfac_array_v_y)
1156 if (associated(landmask)) deallocate(landmask)
1157 if (associated(xlat_array_subgrid)) deallocate(xlat_array_subgrid)
1158 if (associated(xlon_array_subgrid)) deallocate(xlon_array_subgrid)
1159 if (associated(mapfac_array_x_subgrid)) deallocate(mapfac_array_x_subgrid)
1160 if (associated(mapfac_array_y_subgrid)) deallocate(mapfac_array_y_subgrid)
1162 nullify(xlat_ptr)
1163 nullify(xlon_ptr)
1165 end subroutine process_tile
1168 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1169 ! Name: calc_field
1170 !
1171 ! Purpose: This routine fills in the "field" array with interpolated source
1172 ! data. When multiple resolutions of source data are available, an appropriate
1173 ! resolution is chosen automatically. The specified field may either be a
1174 ! continuous field or a categorical field.
1175 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1176 recursive subroutine calc_field(fieldname, field, xlat_array, xlon_array, istagger, &
1177 start_i, end_i, start_j, end_j, start_k, end_k, &
1178 processed_domain, ilevel, landmask, sr_x, sr_y)
1180 use bitarray_module
1181 use interp_module
1182 use llxy_module
1183 use misc_definitions_module
1184 use proc_point_module
1185 use queue_module
1186 use source_data_module
1188 implicit none
1190 ! Arguments
1191 integer, intent(in) :: start_i, end_i, start_j, end_j, start_k, end_k, ilevel, istagger
1192 real, dimension(start_i:end_i, start_j:end_j), intent(in) :: xlat_array, xlon_array
1193 real, dimension(start_i:end_i, start_j:end_j, start_k:end_k), intent(inout) :: field
1194 real, dimension(start_i:end_i, start_j:end_j), intent(in), optional :: landmask
1195 integer, intent(in), optional :: sr_x, sr_y
1196 character (len=128), intent(in) :: fieldname
1197 type (bitarray), intent(inout) :: processed_domain
1199 ! Local variables
1200 integer :: start_src_k, end_src_k
1201 integer :: i, j, k, ix, iy, itype
1202 integer :: user_iproj, istatus
1203 integer :: opt_status
1204 real :: mask_val
1205 real :: temp
1206 real :: scale_factor
1207 real :: msg_val, msg_fill_val, threshold, src_dx, src_dy, dom_dx, dom_dy
1208 real :: user_stand_lon, user_truelat1, user_truelat2, user_dxkm, user_dykm, &
1209 user_known_x, user_known_y, user_known_lat, user_known_lon
1210 real, pointer, dimension(:,:,:) :: data_count
1211 integer, pointer, dimension(:) :: interp_type
1212 integer, pointer, dimension(:) :: interp_opts
1213 character (len=128) :: interp_string
1214 type (bitarray) :: bit_domain, level_domain
1215 type (queue) :: point_queue, tile_queue
1216 type (q_data) :: current_pt
1218 nullify(data_count)
1219 nullify(interp_type)
1220 nullify(interp_opts)
1222 ! If this is the first trip through this routine, we need to allocate the bit array that
1223 ! will persist through all recursive calls, tracking which grid points have been assigned
1224 ! a value.
1225 if (ilevel == 1) call bitarray_create(processed_domain, end_i-start_i+1, end_j-start_j+1)
1227 ! Find out if this "priority level" (given by ilevel) exists
1228 call check_priority_level(fieldname, ilevel, istatus)
1230 ! A bad status indicates that that no data for priority level ilevel is available, and thus, that
1231 ! no further levels will be specified. We are done processing for this level.
1232 if (istatus /= 0) then
1233 if (ilevel == 1) call bitarray_destroy(processed_domain)
1234 return
1235 end if
1237 ! Before proceeding with processing for this level, though, process for the next highest priority level
1238 ! of source data
1239 call calc_field(fieldname, field, xlat_array, xlon_array, istagger, start_i, end_i, &
1240 start_j, end_j, start_k, end_k, processed_domain, ilevel+1, landmask, sr_x, sr_y)
1242 ! At this point, all levels of source data with higher priority have been processed, and we can assign
1243 ! values to all grid points that have not already been given values from higher-priority data
1245 call get_source_opt_status(fieldname, ilevel, opt_status)
1246 if (opt_status == 0) then
1248 ! Find out the projection of the data for this "priority level" (given by ilevel)
1249 call get_data_projection(fieldname, user_iproj, user_stand_lon, user_truelat1, user_truelat2, &
1250 user_dxkm, user_dykm, user_known_x, user_known_y, user_known_lat, &
1251 user_known_lon, ilevel, istatus)
1253 ! A good status indicates that there is data for this priority level, so we store the projection
1254 ! of that data on a stack. The projection will be on the top of the stack (and hence will be
1255 ! the "active" projection) once all higher-priority levels have been processed
1256 call push_source_projection(user_iproj, user_stand_lon, user_truelat1, user_truelat2, &
1257 user_dxkm, user_dykm, user_dykm, user_dxkm, user_known_x, user_known_y, &
1258 user_known_lat, user_known_lon)
1260 ! Initialize point processing module
1261 call proc_point_init()
1263 ! Initialize queues
1264 call q_init(point_queue)
1265 call q_init(tile_queue)
1267 ! Determine whether we will be processing categorical data or continuous data
1268 itype = iget_source_fieldtype(fieldname, ilevel, istatus)
1269 call get_interp_option(fieldname, ilevel, interp_string, istatus)
1270 interp_type => interp_array_from_string(interp_string)
1271 interp_opts => interp_options_from_string(interp_string)
1273 ! Also, check whether we will be using the cell averaging interpolator for continuous fields
1274 if (index(interp_string,'average_gcell') /= 0 .and. itype == CONTINUOUS) then
1275 call get_gcell_threshold(interp_string, threshold, istatus)
1276 if (istatus == 0) then
1277 call get_source_resolution(fieldname, ilevel, src_dx, src_dy, istatus)
1278 if (istatus == 0) then
1279 call get_domain_resolution(dom_dx, dom_dy)
1280 if (gridtype == 'C') then
1281 if (threshold*max(src_dx,src_dy)*111. <= max(dom_dx,dom_dy)/1000.) then
1282 itype = SP_CONTINUOUS
1283 allocate(data_count(start_i:end_i,start_j:end_j,start_k:end_k))
1284 data_count = 0.
1285 end if
1286 else if (gridtype == 'E') then
1287 if (max(src_dx,src_dy) >= threshold*max(dom_dx,dom_dy)) then
1288 itype = SP_CONTINUOUS
1289 allocate(data_count(start_i:end_i,start_j:end_j,start_k:end_k))
1290 data_count = 0.
1291 end if
1292 end if
1293 end if
1294 end if
1295 end if
1297 call get_missing_value(fieldname, ilevel, msg_val, istatus)
1298 if (istatus /= 0) msg_val = NAN
1299 call get_missing_fill_value(fieldname, msg_fill_val, istatus)
1300 if (istatus /= 0) msg_fill_val = NAN
1302 call get_masked_value(fieldname, ilevel, mask_val, istatus)
1303 if (istatus /= 0) mask_val = -1.
1305 if (itype == CONTINUOUS .or. itype == SP_CONTINUOUS) then
1306 call get_source_levels(fieldname, ilevel, start_src_k, end_src_k, istatus)
1307 if (istatus /= 0) return
1308 end if
1310 ! Initialize bitarray used to track which points have been visited and assigned values while
1311 ! processing *this* priority level of data
1312 call bitarray_create(bit_domain, end_i-start_i+1, end_j-start_j+1)
1313 call bitarray_create(level_domain, end_i-start_i+1, end_j-start_j+1)
1315 ! Begin by placing a point in the tile_queue
1316 current_pt%lat = xlat_array(start_i,start_j)
1317 current_pt%lon = xlon_array(start_i,start_j)
1318 current_pt%x = start_i
1319 current_pt%y = start_j
1320 call q_insert(tile_queue, current_pt)
1322 ! While there are still grid points in tiles that have not yet been processed
1323 do while (q_isdata(tile_queue))
1325 ! Take a point from the outer queue and place it in the point_queue for processing
1326 current_pt = q_remove(tile_queue)
1328 ! If this level of data is categorical (i.e., is given as an array of category indices),
1329 ! then first try to process the entire tile in one call to accum_categorical. Any grid
1330 ! points that are not given values by accum_categorical and that lie within the current
1331 ! tile of source data are individually assigned values in the inner loop
1332 if (itype == CATEGORICAL) then
1334 ! Have we already visited this point? If so, this tile has already been processed by
1335 ! accum_categorical.
1336 if (.not. bitarray_test(bit_domain, current_pt%x-start_i+1, current_pt%y-start_j+1)) then
1337 call q_insert(point_queue, current_pt)
1338 if (.not. have_processed_tile(current_pt%lat, current_pt%lon, fieldname, ilevel)) then
1339 call accum_categorical(current_pt%lat, current_pt%lon, istagger, field, &
1340 start_i, end_i, start_j, end_j, start_k, end_k, &
1341 fieldname, processed_domain, level_domain, &
1342 ilevel, msg_val, mask_val, sr_x, sr_y)
1343 ! BUG: Where do we mask out those points that are on land/water when masked=land/water is set?
1344 end if
1345 call bitarray_set(bit_domain, current_pt%x-start_i+1, current_pt%y-start_j+1)
1346 end if
1348 else if (itype == SP_CONTINUOUS) then
1350 ! Have we already visited this point? If so, this tile has already been processed by
1351 ! accum_continuous.
1352 if (.not. bitarray_test(bit_domain, current_pt%x-start_i+1, current_pt%y-start_j+1)) then
1353 call q_insert(point_queue, current_pt)
1354 if (.not. have_processed_tile(current_pt%lat, current_pt%lon, fieldname, ilevel)) then
1355 call accum_continuous(current_pt%lat, current_pt%lon, istagger, field, data_count, &
1356 start_i, end_i, start_j, end_j, start_k, end_k, &
1357 fieldname, processed_domain, level_domain, &
1358 ilevel, msg_val, mask_val, sr_x, sr_y)
1359 ! BUG: Where do we mask out those points that are on land/water when masked=land/water is set?
1360 end if
1361 call bitarray_set(bit_domain, current_pt%x-start_i+1, current_pt%y-start_j+1)
1362 end if
1364 else if (itype == CONTINUOUS) then
1366 ! Have we already visited this point? If so, the tile containing this point has already been
1367 ! processed in the inner loop.
1368 if (.not. bitarray_test(bit_domain, current_pt%x-start_i+1, current_pt%y-start_j+1)) then
1369 call q_insert(point_queue, current_pt)
1370 call bitarray_set(bit_domain, current_pt%x-start_i+1, current_pt%y-start_j+1)
1371 end if
1373 end if
1375 ! This inner loop, where all grid points contained in the current source tile are processed
1376 do while (q_isdata(point_queue))
1377 current_pt = q_remove(point_queue)
1378 ix = current_pt%x
1379 iy = current_pt%y
1381 ! Process the current point
1382 if (itype == CONTINUOUS .or. itype == SP_CONTINUOUS) then
1384 ! Have we already assigned this point a value from this priority level?
1385 if (.not. bitarray_test(level_domain, ix-start_i+1, iy-start_j+1)) then
1387 ! If the point was already assigned a value from a higher-priority level, no
1388 ! need to assign a new value
1389 if (bitarray_test(processed_domain, ix-start_i+1, iy-start_j+1)) then
1390 call bitarray_set(level_domain, ix-start_i+1, iy-start_j+1)
1392 ! Otherwise, need to assign values from this level of source data if we can
1393 else
1394 if (present(landmask) .and. (istagger == M .or. istagger == HH)) then
1395 if (landmask(ix,iy) /= mask_val) then
1396 do k=start_src_k,end_src_k
1397 temp = get_point(current_pt%lat, current_pt%lon, k, &
1398 fieldname, ilevel, interp_type, interp_opts, msg_val)
1399 if (temp /= msg_val) then
1400 field(ix, iy, k) = temp
1401 call bitarray_set(level_domain, ix-start_i+1, iy-start_j+1)
1402 if (itype == SP_CONTINUOUS) data_count(ix, iy, k) = 1.0
1403 else
1404 field(ix, iy, k) = msg_fill_val
1405 end if
1406 end do
1407 else
1408 do k=start_k,end_k
1409 field(ix,iy,k) = msg_fill_val
1410 end do
1411 end if
1412 else
1413 do k=start_src_k,end_src_k
1414 temp = get_point(current_pt%lat, current_pt%lon, k, &
1415 fieldname, ilevel, interp_type, interp_opts, msg_val)
1416 if (temp /= msg_val) then
1417 field(ix, iy, k) = temp
1418 call bitarray_set(level_domain, ix-start_i+1, iy-start_j+1)
1419 if (itype == SP_CONTINUOUS) data_count(ix, iy, k) = 1.0
1420 else
1421 field(ix, iy, k) = msg_fill_val
1422 end if
1423 end do
1424 end if
1425 end if
1426 end if
1428 else if (itype == CATEGORICAL) then
1430 ! Have we already assigned this point a value from this priority level?
1431 if (.not.bitarray_test(level_domain, ix-start_i+1, iy-start_j+1)) then
1433 ! If the point was already assigned a value from a higher-priority level, no
1434 ! need to assign a new value
1435 if (bitarray_test(processed_domain, ix-start_i+1, iy-start_j+1)) then
1436 call bitarray_set(level_domain, ix-start_i+1, iy-start_j+1)
1438 ! Otherwise, the point was apparently not given a value when accum_categorical
1439 ! was called for the current tile, and we need to assign values from this
1440 ! level of source data if we can
1441 else
1442 if (present(landmask) .and. (istagger == M .or. istagger == HH)) then
1443 if (landmask(ix,iy) /= mask_val) then
1444 temp = get_point(current_pt%lat, current_pt%lon, 1, &
1445 fieldname, ilevel, interp_type, interp_opts, msg_val)
1447 do k=start_k,end_k
1448 field(ix,iy,k) = 0.
1449 end do
1451 if (temp /= msg_val) then
1452 if (int(temp) >= start_k .and. int(temp) <= end_k) then
1453 field(ix, iy, int(temp)) = field(ix, iy, int(temp)) + 1.
1454 else
1455 call mprintf(.true.,WARN,' Attempted to assign an invalid category '// &
1456 '%i to grid point (%i, %i)', i1=int(temp), i2=ix, i3=iy)
1457 end if
1458 call bitarray_set(level_domain, ix-start_i+1, iy-start_j+1)
1459 end if
1461 else
1462 do k=start_k,end_k
1463 field(ix,iy,k) = 0.
1464 end do
1465 end if
1466 else
1467 temp = get_point(current_pt%lat, current_pt%lon, 1, &
1468 fieldname, ilevel, interp_type, interp_opts, msg_val)
1470 do k=start_k,end_k
1471 field(ix,iy,k) = 0.
1472 end do
1474 if (temp /= msg_val) then
1475 if (int(temp) >= start_k .and. int(temp) <= end_k) then
1476 field(ix, iy, int(temp)) = field(ix, iy, int(temp)) + 1.
1477 else
1478 call mprintf(.true.,WARN,' Attempted to assign an invalid category '// &
1479 '%i to grid point (%i, %i)', i1=int(temp), i2=ix, i3=iy)
1480 end if
1481 call bitarray_set(level_domain, ix-start_i+1, iy-start_j+1)
1482 end if
1483 end if
1484 end if
1485 end if
1487 end if
1489 ! Scan neighboring points, adding them to the appropriate queue based on whether they
1490 ! are in the current tile or not
1491 if (iy > start_j) then
1492 if (ix > start_i) then
1494 ! Neighbor with relative position (-1,-1)
1495 call process_neighbor(ix-1, iy-1, bit_domain, point_queue, tile_queue, &
1496 xlat_array, xlon_array, start_i, end_i, start_j, end_j, ilevel)
1497 end if
1499 ! Neighbor with relative position (0,-1)
1500 call process_neighbor(ix, iy-1, bit_domain, point_queue, tile_queue, &
1501 xlat_array, xlon_array, start_i, end_i, start_j, end_j, ilevel)
1503 if (ix < end_i) then
1505 ! Neighbor with relative position (+1,-1)
1506 call process_neighbor(ix+1, iy-1, bit_domain, point_queue, tile_queue, &
1507 xlat_array, xlon_array, start_i, end_i, start_j, end_j, ilevel)
1508 end if
1509 end if
1511 if (ix > start_i) then
1513 ! Neighbor with relative position (-1,0)
1514 call process_neighbor(ix-1, iy, bit_domain, point_queue, tile_queue, &
1515 xlat_array, xlon_array, start_i, end_i, start_j, end_j, ilevel)
1516 end if
1518 if (ix < end_i) then
1520 ! Neighbor with relative position (+1,0)
1521 call process_neighbor(ix+1, iy, bit_domain, point_queue, tile_queue, &
1522 xlat_array, xlon_array, start_i, end_i, start_j, end_j, ilevel)
1523 end if
1525 if (iy < end_j) then
1526 if (ix > start_i) then
1528 ! Neighbor with relative position (-1,+1)
1529 call process_neighbor(ix-1, iy+1, bit_domain, point_queue, tile_queue, &
1530 xlat_array, xlon_array, start_i, end_i, start_j, end_j, ilevel)
1531 end if
1533 ! Neighbor with relative position (0,+1)
1534 call process_neighbor(ix, iy+1, bit_domain, point_queue, tile_queue, &
1535 xlat_array, xlon_array, start_i, end_i, start_j, end_j, ilevel)
1536 if (ix < end_i) then
1538 ! Neighbor with relative position (+1,+1)
1539 call process_neighbor(ix+1, iy+1, bit_domain, point_queue, tile_queue, &
1540 xlat_array, xlon_array, start_i, end_i, start_j, end_j, ilevel)
1541 end if
1542 end if
1544 end do
1546 end do
1548 if (itype == SP_CONTINUOUS) then
1549 itype = CONTINUOUS
1550 if (present(landmask) .and. (istagger == M .or. istagger == HH)) then
1551 do j=start_j,end_j
1552 do i=start_i,end_i
1553 if (landmask(i,j) /= mask_val) then
1554 do k=start_k,end_k
1555 if (data_count(i,j,k) > 0.) then
1556 field(i,j,k) = field(i,j,k) / data_count(i,j,k)
1557 else
1558 if (.not.bitarray_test(processed_domain, i-start_i+1, j-start_j+1)) then
1559 field(i,j,k) = msg_fill_val
1560 end if
1561 end if
1562 end do
1563 else
1564 if (.not.bitarray_test(processed_domain, i-start_i+1, j-start_j+1)) then
1565 do k=start_k,end_k
1566 field(i,j,k) = msg_fill_val
1567 end do
1568 end if
1569 end if
1570 end do
1571 end do
1572 else
1573 do k=start_k,end_k
1574 do j=start_j,end_j
1575 do i=start_i,end_i
1576 if (data_count(i,j,k) > 0.) then
1577 field(i,j,k) = field(i,j,k) / data_count(i,j,k)
1578 else
1579 if (.not.bitarray_test(processed_domain, i-start_i+1, j-start_j+1)) then
1580 field(i,j,k) = msg_fill_val
1581 end if
1582 end if
1583 end do
1584 end do
1585 end do
1586 end if
1587 deallocate(data_count)
1589 else if (itype == CATEGORICAL) then
1590 if (present(landmask) .and. (istagger == M .or. istagger == HH)) then
1591 do j=start_j,end_j
1592 do i=start_i,end_i
1593 if (landmask(i,j) == mask_val) then
1594 do k=start_k,end_k
1595 field(i,j,k) = 0.
1596 end do
1597 end if
1598 end do
1599 end do
1600 end if
1601 end if
1603 deallocate(interp_type)
1604 deallocate(interp_opts)
1607 ! We may need to scale this field by a constant
1608 call get_field_scale_factor(fieldname, ilevel, scale_factor, istatus)
1609 if (istatus == 0) then
1610 do i=start_i, end_i
1611 do j=start_j, end_j
1612 if (bitarray_test(level_domain,i-start_i+1,j-start_j+1) .and. &
1613 .not. bitarray_test(processed_domain,i-start_i+1,j-start_j+1)) then
1614 do k=start_k,end_k
1615 if (field(i,j,k) /= msg_fill_val) then
1616 field(i,j,k) = field(i,j,k) * scale_factor
1617 end if
1618 end do
1619 end if
1620 end do
1621 end do
1622 end if
1625 ! Now add the points that were assigned values at this priority level to the complete array
1626 ! of points that have been assigned values
1627 call bitarray_merge(processed_domain, level_domain)
1629 call bitarray_destroy(bit_domain)
1630 call bitarray_destroy(level_domain)
1631 call q_destroy(point_queue)
1632 call q_destroy(tile_queue)
1633 call proc_point_shutdown()
1635 ! Remove the projection of the current level of source data from the stack, thus "activating"
1636 ! the projection of the next highest level
1637 call pop_source_projection()
1639 else
1640 call mprintf(.true.,STDOUT,' Important note: could not open input dataset for priority level %i, '// &
1641 'but this level is optional.', i1=ilevel)
1642 call mprintf(.true.,LOGFILE,' Important note: could not open input dataset for priority level %i, '// &
1643 'but this level is optional.', i1=ilevel)
1644 end if
1646 ! If this is the last level of the recursion, we can also deallocate processed_domain
1647 if (ilevel == 1) call bitarray_destroy(processed_domain)
1649 end subroutine calc_field
1652 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1653 ! Name: get_lat_lon_fields
1654 !
1655 ! Purpose: To calculate the latitude and longitude for every gridpoint in the
1656 ! tile of the model domain. The caller may specify that the grid for which
1657 ! values are computed is staggered or unstaggered using the "stagger"
1658 ! argument.
1659 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1660 subroutine get_lat_lon_fields(xlat_arr, xlon_arr, start_mem_i, &
1661 start_mem_j, end_mem_i, end_mem_j, stagger, comp_ll, &
1662 sub_x, sub_y)
1664 use llxy_module
1665 use misc_definitions_module
1667 implicit none
1669 ! Arguments
1670 integer, intent(in) :: start_mem_i, start_mem_j, end_mem_i, &
1671 end_mem_j, stagger
1672 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(out) :: xlat_arr, xlon_arr
1673 logical, optional, intent(in) :: comp_ll
1674 integer, optional, intent(in) :: sub_x, sub_y
1676 ! Local variables
1677 integer :: i, j
1678 real :: rx, ry
1680 rx = 1.0
1681 ry = 1.0
1682 if (present(sub_x)) rx = real(sub_x)
1683 if (present(sub_y)) ry = real(sub_y)
1685 do i=start_mem_i, end_mem_i
1686 do j=start_mem_j, end_mem_j
1687 call xytoll(real(i-1)/rx+1.0, real(j-1)/ry+1.0, &
1688 xlat_arr(i,j), xlon_arr(i,j), stagger, comp_ll=comp_ll)
1689 end do
1690 end do
1692 end subroutine get_lat_lon_fields
1695 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1696 ! Name: get_map_factor
1697 !
1698 ! Purpose: Given the latitude field, this routine calculates map factors for
1699 ! the grid points of the specified domain. For different grids (e.g., C grid,
1700 ! E grid), the latitude array should provide the latitudes of the points for
1701 ! which map factors are to be calculated.
1702 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1703 subroutine get_map_factor(xlat_arr, xlon_arr, mapfac_arr_x, mapfac_arr_y, &
1704 start_mem_i, start_mem_j, end_mem_i, end_mem_j)
1706 use constants_module
1707 use gridinfo_module
1708 use misc_definitions_module
1709 use map_utils
1711 implicit none
1713 ! Arguments
1714 integer, intent(in) :: start_mem_i, start_mem_j, end_mem_i, end_mem_j
1715 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(in) :: xlat_arr, xlon_arr
1716 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(out) :: mapfac_arr_x
1717 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(out) :: mapfac_arr_y
1719 ! Local variables
1720 integer :: i, j
1721 real :: n, colat, colat0, colat1, colat2, comp_lat, comp_lon
1723 !
1724 ! Equations for map factor given in Principles of Meteorological Analysis,
1725 ! Walter J. Saucier, pp. 32-33
1726 !
1728 ! Lambert conformal projection
1729 if (iproj_type == PROJ_LC) then
1730 if (truelat1 /= truelat2) then
1731 colat1 = rad_per_deg*(90.0 - truelat1)
1732 colat2 = rad_per_deg*(90.0 - truelat2)
1733 n = (log(sin(colat1)) - log(sin(colat2))) &
1734 / (log(tan(colat1/2.0)) - log(tan(colat2/2.0)))
1736 do i=start_mem_i, end_mem_i
1737 do j=start_mem_j, end_mem_j
1738 colat = rad_per_deg*(90.0 - xlat_arr(i,j))
1739 mapfac_arr_x(i,j) = sin(colat2)/sin(colat)*(tan(colat/2.0)/tan(colat2/2.0))**n
1740 mapfac_arr_y(i,j) = mapfac_arr_x(i,j)
1741 end do
1742 end do
1744 else
1745 colat0 = rad_per_deg*(90.0 - truelat1)
1747 do i=start_mem_i, end_mem_i
1748 do j=start_mem_j, end_mem_j
1749 colat = rad_per_deg*(90.0 - xlat_arr(i,j))
1750 mapfac_arr_x(i,j) = sin(colat0)/sin(colat)*(tan(colat/2.0)/tan(colat0/2.0))**cos(colat0)
1751 mapfac_arr_y(i,j) = mapfac_arr_x(i,j)
1752 end do
1753 end do
1755 end if
1757 ! Polar stereographic projection
1758 else if (iproj_type == PROJ_PS) then
1760 do i=start_mem_i, end_mem_i
1761 do j=start_mem_j, end_mem_j
1762 mapfac_arr_x(i,j) = (1.0 + sin(rad_per_deg*abs(truelat1)))/(1.0 + sin(rad_per_deg*sign(1.,truelat1)*xlat_arr(i,j)))
1763 mapfac_arr_y(i,j) = mapfac_arr_x(i,j)
1764 end do
1765 end do
1767 ! Mercator projection
1768 else if (iproj_type == PROJ_MERC) then
1769 colat0 = rad_per_deg*(90.0 - truelat1)
1771 do i=start_mem_i, end_mem_i
1772 do j=start_mem_j, end_mem_j
1773 colat = rad_per_deg*(90.0 - xlat_arr(i,j))
1774 mapfac_arr_x(i,j) = sin(colat0) / sin(colat)
1775 mapfac_arr_y(i,j) = mapfac_arr_x(i,j)
1776 end do
1777 end do
1779 ! Global cylindrical projection
1780 else if (iproj_type == PROJ_CYL) then
1782 do i=start_mem_i, end_mem_i
1783 do j=start_mem_j, end_mem_j
1784 if (abs(xlat_arr(i,j)) == 90.0) then
1785 mapfac_arr_x(i,j) = 0. ! MSF actually becomes infinite at poles, but
1786 ! the values should never be used there; by
1787 ! setting to 0, we hope to induce a "divide
1788 ! by zero" error if they are
1789 else
1790 mapfac_arr_x(i,j) = 1.0 / cos(xlat_arr(i,j)*rad_per_deg)
1791 end if
1792 mapfac_arr_y(i,j) = 1.0
1793 end do
1794 end do
1796 ! Rotated global cylindrical projection
1797 else if (iproj_type == PROJ_CASSINI) then
1799 if (abs(pole_lat) == 90.) then
1800 do i=start_mem_i, end_mem_i
1801 do j=start_mem_j, end_mem_j
1802 if (abs(xlat_arr(i,j)) >= 90.0) then
1803 mapfac_arr_x(i,j) = 0. ! MSF actually becomes infinite at poles, but
1804 ! the values should never be used there; by
1805 ! setting to 0, we hope to induce a "divide
1806 ! by zero" error if they are
1807 else
1808 mapfac_arr_x(i,j) = 1.0 / cos(xlat_arr(i,j)*rad_per_deg)
1809 end if
1810 mapfac_arr_y(i,j) = 1.0
1811 end do
1812 end do
1813 else
1814 do i=start_mem_i, end_mem_i
1815 do j=start_mem_j, end_mem_j
1816 call rotate_coords(xlat_arr(i,j),xlon_arr(i,j), &
1817 comp_lat, comp_lon, &
1818 pole_lat, pole_lon, stand_lon, &
1819 -1)
1820 if (abs(comp_lat) >= 90.0) then
1821 mapfac_arr_x(i,j) = 0. ! MSF actually becomes infinite at poles, but
1822 ! the values should never be used there; by
1823 ! setting to 0, we hope to induce a "divide
1824 ! by zero" error if they are
1825 else
1826 mapfac_arr_x(i,j) = 1.0 / cos(comp_lat*rad_per_deg)
1827 end if
1828 mapfac_arr_y(i,j) = 1.0
1829 end do
1830 end do
1831 end if
1833 else if (iproj_type == PROJ_ROTLL) then
1835 do i=start_mem_i, end_mem_i
1836 do j=start_mem_j, end_mem_j
1837 mapfac_arr_x(i,j) = 1.0
1838 mapfac_arr_y(i,j) = 1.0
1839 end do
1840 end do
1842 end if
1844 end subroutine get_map_factor
1847 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1848 ! Name: get_coriolis_parameters
1849 !
1850 ! Purpose: To calculate the Coriolis parameters f and e for every gridpoint in
1851 ! the tile of the model domain
1852 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1853 subroutine get_coriolis_parameters(xlat_arr, f, e, &
1854 start_mem_i, start_mem_j, end_mem_i, end_mem_j)
1856 use constants_module
1858 implicit none
1860 ! Arguments
1861 integer, intent(in) :: start_mem_i, start_mem_j, end_mem_i, end_mem_j
1862 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(in) :: xlat_arr
1863 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(out) :: f, e
1865 ! Local variables
1866 integer :: i, j
1868 do i=start_mem_i, end_mem_i
1869 do j=start_mem_j, end_mem_j
1871 f(i,j) = 2.0*OMEGA_E*sin(rad_per_deg*xlat_arr(i,j))
1872 e(i,j) = 2.0*OMEGA_E*cos(rad_per_deg*xlat_arr(i,j))
1874 end do
1875 end do
1877 end subroutine get_coriolis_parameters
1880 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1881 ! Name: get_rotang
1882 !
1883 ! Purpose: To calculate the sine and cosine of rotation angle.
1884 !
1885 ! NOTES: The formulas used in this routine come from those in the
1886 ! vecrot_rotlat() routine of the original WRF SI.
1887 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1888 subroutine get_rotang(xlat_arr, xlon_arr, cosa, sina, &
1889 start_mem_i, start_mem_j, end_mem_i, end_mem_j)
1891 use constants_module
1892 use gridinfo_module
1894 implicit none
1896 ! Arguments
1897 integer, intent(in) :: start_mem_i, start_mem_j, end_mem_i, end_mem_j
1898 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(in) :: xlat_arr, xlon_arr
1899 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(out) :: cosa, sina
1901 ! Local variables
1902 integer :: i, j
1903 real :: alpha, d_lon
1905 do i=start_mem_i, end_mem_i
1906 do j=start_mem_j+1, end_mem_j-1
1907 d_lon = xlon_arr(i,j+1)-xlon_arr(i,j-1)
1908 if (d_lon > 180.) then
1909 d_lon = d_lon - 360.
1910 else if (d_lon < -180.) then
1911 d_lon = d_lon + 360.
1912 end if
1914 alpha = atan2(-cos(xlat_arr(i,j)*RAD_PER_DEG) * (d_lon*RAD_PER_DEG), &
1915 ((xlat_arr(i,j+1)-xlat_arr(i,j-1))*RAD_PER_DEG))
1916 sina(i,j) = sin(alpha)
1917 cosa(i,j) = cos(alpha)
1918 end do
1919 end do
1921 do i=start_mem_i, end_mem_i
1922 d_lon = xlon_arr(i,start_mem_j+1)-xlon_arr(i,start_mem_j)
1923 if (d_lon > 180.) then
1924 d_lon = d_lon - 360.
1925 else if (d_lon < -180.) then
1926 d_lon = d_lon + 360.
1927 end if
1929 alpha = atan2(-cos(xlat_arr(i,start_mem_j)*RAD_PER_DEG) * (d_lon*RAD_PER_DEG), &
1930 ((xlat_arr(i,start_mem_j+1)-xlat_arr(i,start_mem_j))*RAD_PER_DEG))
1931 sina(i,start_mem_j) = sin(alpha)
1932 cosa(i,start_mem_j) = cos(alpha)
1933 end do
1935 do i=start_mem_i, end_mem_i
1936 d_lon = xlon_arr(i,end_mem_j)-xlon_arr(i,end_mem_j-1)
1937 if (d_lon > 180.) then
1938 d_lon = d_lon - 360.
1939 else if (d_lon < -180.) then
1940 d_lon = d_lon + 360.
1941 end if
1943 alpha = atan2(-cos(xlat_arr(i,end_mem_j)*RAD_PER_DEG) * (d_lon*RAD_PER_DEG), &
1944 ((xlat_arr(i,end_mem_j)-xlat_arr(i,end_mem_j-1))*RAD_PER_DEG))
1945 sina(i,end_mem_j) = sin(alpha)
1946 cosa(i,end_mem_j) = cos(alpha)
1947 end do
1949 end subroutine get_rotang
1952 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1953 ! Name: process_neighbor
1954 !
1955 ! Purpose: This routine, give the x/y location of a point, determines whether
1956 ! the point has already been processed, and if not, which processing queue
1957 ! the point should be placed in.
1958 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1959 subroutine process_neighbor(ix, iy, bit_domain, point_queue, tile_queue, &
1960 xlat_array, xlon_array, &
1961 start_i, end_i, start_j, end_j, ilevel)
1963 use bitarray_module
1964 use misc_definitions_module
1965 use proc_point_module
1966 use queue_module
1968 implicit none
1970 ! Arguments
1971 integer, intent(in) :: ix, iy, start_i, end_i, start_j, end_j, ilevel
1972 real, dimension(start_i:end_i, start_j:end_j), intent(in) :: xlat_array, xlon_array
1973 type (bitarray), intent(inout) :: bit_domain
1974 type (queue), intent(inout) :: point_queue, tile_queue
1976 ! Local variables
1977 type (q_data) :: process_pt
1978 logical :: is_in_tile
1980 ! If the point has already been visited, no need to do anything more.
1981 if (.not. bitarray_test(bit_domain, ix-start_i+1, iy-start_j+1)) then
1983 ! Create a queue item for the current point
1984 process_pt%lat = xlat_array(ix,iy)
1985 process_pt%lon = xlon_array(ix,iy)
1986 process_pt%x = ix
1987 process_pt%y = iy
1989 is_in_tile = is_point_in_tile(process_pt%lat, process_pt%lon, ilevel)
1991 ! If the point is in the current tile, add it to the list of points
1992 ! to be processed in the inner loop
1993 if (is_in_tile) then
1994 call q_insert(point_queue, process_pt)
1995 call bitarray_set(bit_domain, ix-start_i+1, iy-start_j+1)
1997 ! Otherwise, we will process this point later. Add it to the list for
1998 ! the outer loop.
1999 else
2000 call q_insert(tile_queue, process_pt)
2001 end if
2003 end if
2005 end subroutine process_neighbor
2008 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2009 ! Name: calc_dfdy
2010 !
2011 ! Purpose: This routine calculates df/dy for the field in src_arr, and places
2012 ! the result in dst_array.
2013 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2014 subroutine calc_dfdy(src_arr, dst_arr, start_mem_i, start_mem_j, start_mem_k, &
2015 end_mem_i, end_mem_j, end_mem_k, mapfac)
2017 ! Modules
2018 use gridinfo_module
2020 implicit none
2022 ! Arguments
2023 integer, intent(in) :: start_mem_i, start_mem_j, start_mem_k, end_mem_i, end_mem_j, end_mem_k
2024 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j,start_mem_k:end_mem_k), intent(in) :: src_arr
2025 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j,start_mem_k:end_mem_k), intent(out) :: dst_arr
2026 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(in), optional :: mapfac
2028 ! Local variables
2029 integer :: i, j, k
2031 if (present(mapfac)) then
2032 do k=start_mem_k,end_mem_k
2033 do i=start_mem_i, end_mem_i
2034 do j=start_mem_j+1, end_mem_j-1
2035 dst_arr(i,j,k) = (src_arr(i,j+1,k) - src_arr(i,j-1,k))/(2.*dykm*mapfac(i,j))
2036 end do
2037 end do
2039 do i=start_mem_i, end_mem_i
2040 dst_arr(i,start_mem_j,k) = (src_arr(i,start_mem_j+1,k) - src_arr(i,start_mem_j,k))/(dykm*mapfac(i,j))
2041 end do
2043 do i=start_mem_i, end_mem_i
2044 dst_arr(i,end_mem_j,k) = (src_arr(i,end_mem_j,k) - src_arr(i,end_mem_j-1,k))/(dykm*mapfac(i,j))
2045 end do
2046 end do
2047 else
2048 do k=start_mem_k,end_mem_k
2049 do i=start_mem_i, end_mem_i
2050 do j=start_mem_j+1, end_mem_j-1
2051 dst_arr(i,j,k) = (src_arr(i,j+1,k) - src_arr(i,j-1,k))/(2.*dykm)
2052 end do
2053 end do
2055 do i=start_mem_i, end_mem_i
2056 dst_arr(i,start_mem_j,k) = (src_arr(i,start_mem_j+1,k) - src_arr(i,start_mem_j,k))/(dykm)
2057 end do
2059 do i=start_mem_i, end_mem_i
2060 dst_arr(i,end_mem_j,k) = (src_arr(i,end_mem_j,k) - src_arr(i,end_mem_j-1,k))/(dykm)
2061 end do
2062 end do
2063 end if
2065 end subroutine calc_dfdy
2068 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2069 ! Name: calc_dfdx
2070 !
2071 ! Purpose: This routine calculates df/dx for the field in src_arr, and places
2072 ! the result in dst_array.
2073 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2074 subroutine calc_dfdx(src_arr, dst_arr, start_mem_i, start_mem_j, &
2075 start_mem_k, end_mem_i, end_mem_j, end_mem_k, mapfac)
2077 ! Modules
2078 use gridinfo_module
2080 implicit none
2082 ! Arguments
2083 integer, intent(in) :: start_mem_i, start_mem_j, start_mem_k, end_mem_i, end_mem_j, end_mem_k
2084 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j, start_mem_k:end_mem_k), intent(in) :: src_arr
2085 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j, start_mem_k:end_mem_k), intent(out) :: dst_arr
2086 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(in), optional :: mapfac
2088 ! Local variables
2089 integer :: i, j, k
2091 if (present(mapfac)) then
2092 do k=start_mem_k, end_mem_k
2093 do i=start_mem_i+1, end_mem_i-1
2094 do j=start_mem_j, end_mem_j
2095 dst_arr(i,j,k) = (src_arr(i+1,j,k) - src_arr(i-1,j,k))/(2.*dxkm*mapfac(i,j))
2096 end do
2097 end do
2099 do j=start_mem_j, end_mem_j
2100 dst_arr(start_mem_i,j,k) = (src_arr(start_mem_i+1,j,k) - src_arr(start_mem_i,j,k))/(dxkm*mapfac(i,j))
2101 end do
2103 do j=start_mem_j, end_mem_j
2104 dst_arr(end_mem_i,j,k) = (src_arr(end_mem_i,j,k) - src_arr(end_mem_i-1,j,k))/(dxkm*mapfac(i,j))
2105 end do
2106 end do
2107 else
2108 do k=start_mem_k, end_mem_k
2109 do i=start_mem_i+1, end_mem_i-1
2110 do j=start_mem_j, end_mem_j
2111 dst_arr(i,j,k) = (src_arr(i+1,j,k) - src_arr(i-1,j,k))/(2.*dxkm)
2112 end do
2113 end do
2115 do j=start_mem_j, end_mem_j
2116 dst_arr(start_mem_i,j,k) = (src_arr(start_mem_i+1,j,k) - src_arr(start_mem_i,j,k))/(dxkm)
2117 end do
2119 do j=start_mem_j, end_mem_j
2120 dst_arr(end_mem_i,j,k) = (src_arr(end_mem_i,j,k) - src_arr(end_mem_i-1,j,k))/(dxkm)
2121 end do
2122 end do
2123 end if
2125 end subroutine calc_dfdx
2127 end module process_tile_module