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