| blob | 78aac5735c8d6a3afd9ef54226913b02a2b9af02 |
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
78 datestr = '0000-00-00_00:00:00'
79 field_count = 0
80 mass_flag=1.0
82 ! The following pertains primarily to the C grid
83 ! Determine whether only (n-1)th rows/columns should be computed for variables
84 ! on staggered grid. In a distributed memory situation, not every tile should
85 ! have only (n-1)th rows/columns computed, or we end up with (n-k)
86 ! rows/columns when there are k patches in the y/x direction
87 if (extra_col) then
88 start_patch_i = dummy_start_patch_i ! The seemingly pointless renaming of start
89 end_patch_i = dummy_end_patch_i - 1 ! naming convention with modified end_patch variables,
90 end_patch_stag_i = dummy_end_patch_i ! variables is so that we can maintain consistent
91 ! which are marked as intent(in)
92 start_mem_i = start_patch_i - HALO_WIDTH
93 end_mem_i = end_patch_i + HALO_WIDTH
94 end_mem_stag_i = end_patch_stag_i + HALO_WIDTH
95 else
96 start_patch_i = dummy_start_patch_i
97 end_patch_i = dummy_end_patch_i
98 end_patch_stag_i = dummy_end_patch_i
100 start_mem_i = start_patch_i - HALO_WIDTH
101 end_mem_i = end_patch_i + HALO_WIDTH
102 end_mem_stag_i = end_patch_stag_i + HALO_WIDTH
103 end if
105 if (extra_row) then
106 start_patch_j = dummy_start_patch_j
107 end_patch_j = dummy_end_patch_j - 1
108 end_patch_stag_j = dummy_end_patch_j
110 start_mem_j = start_patch_j - HALO_WIDTH
111 end_mem_j = end_patch_j + HALO_WIDTH
112 end_mem_stag_j = end_patch_stag_j + HALO_WIDTH
113 else
114 start_patch_j = dummy_start_patch_j
115 end_patch_j = dummy_end_patch_j
116 end_patch_stag_j = dummy_end_patch_j
118 start_mem_j = start_patch_j - HALO_WIDTH
119 end_mem_j = end_patch_j + HALO_WIDTH
120 end_mem_stag_j = end_patch_stag_j + HALO_WIDTH
121 end if
123 start_dom_i = dummy_start_dom_i
124 if (grid_type == 'C') then
125 end_dom_i = dummy_end_dom_i - 1
126 end_dom_stag_i = dummy_end_dom_i
127 else if (grid_type == 'E') then
128 end_dom_i = dummy_end_dom_i
129 end_dom_stag_i = dummy_end_dom_i
130 end if
132 start_dom_j = dummy_start_dom_j
133 if (grid_type == 'C') then
134 end_dom_j = dummy_end_dom_j - 1
135 end_dom_stag_j = dummy_end_dom_j
136 else if (grid_type == 'E') then
137 end_dom_j = dummy_end_dom_j
138 end_dom_stag_j = dummy_end_dom_j
139 end if
141 ! Allocate arrays to hold all lat/lon fields; these will persist for the duration of
142 ! the process_tile routine
143 ! For C grid, we have M, U, and V points
144 ! For E grid, we have only M and V points
145 allocate(xlat_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
146 allocate(xlon_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
147 allocate(xlat_array_v(start_mem_i:end_mem_i, start_mem_j:end_mem_stag_j))
148 allocate(xlon_array_v(start_mem_i:end_mem_i, start_mem_j:end_mem_stag_j))
149 if (grid_type == 'C') then
150 allocate(xlat_array_u(start_mem_i:end_mem_stag_i, start_mem_j:end_mem_j))
151 allocate(xlon_array_u(start_mem_i:end_mem_stag_i, start_mem_j:end_mem_j))
152 allocate(clat_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
153 allocate(clon_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
154 end if
155 nullify(xlat_array_subgrid)
156 nullify(xlon_array_subgrid)
157 nullify(mapfac_array_x_subgrid)
158 nullify(mapfac_array_y_subgrid)
160 ! Initialize hash table to track which fields have been processed
161 call hash_init(processed_fieldnames)
163 !
164 ! Calculate lat/lon for every point in the tile (XLAT and XLON)
165 ! The xlat_array and xlon_array arrays will be used in processing other fields
166 !
167 call mprintf(.true.,STDOUT,' Processing XLAT and XLONG')
169 if (grid_type == 'C') then
170 call get_lat_lon_fields(xlat_array, xlon_array, start_mem_i, &
171 start_mem_j, end_mem_i, end_mem_j, M)
172 call get_lat_lon_fields(xlat_array_v, xlon_array_v, start_mem_i, &
173 start_mem_j, end_mem_i, end_mem_stag_j, V)
174 call get_lat_lon_fields(xlat_array_u, xlon_array_u, start_mem_i, &
175 start_mem_j, end_mem_stag_i, end_mem_j, U)
176 call get_lat_lon_fields(clat_array, clon_array, start_mem_i, &
177 start_mem_j, end_mem_i, end_mem_j, M, comp_ll=.true.)
179 corner_lats(1) = xlat_array(start_patch_i,start_patch_j)
180 corner_lats(2) = xlat_array(start_patch_i,end_patch_j)
181 corner_lats(3) = xlat_array(end_patch_i,end_patch_j)
182 corner_lats(4) = xlat_array(end_patch_i,start_patch_j)
184 corner_lats(5) = xlat_array_u(start_patch_i,start_patch_j)
185 corner_lats(6) = xlat_array_u(start_patch_i,end_patch_j)
186 corner_lats(7) = xlat_array_u(end_patch_stag_i,end_patch_j)
187 corner_lats(8) = xlat_array_u(end_patch_stag_i,start_patch_j)
189 corner_lats(9) = xlat_array_v(start_patch_i,start_patch_j)
190 corner_lats(10) = xlat_array_v(start_patch_i,end_patch_stag_j)
191 corner_lats(11) = xlat_array_v(end_patch_i,end_patch_stag_j)
192 corner_lats(12) = xlat_array_v(end_patch_i,start_patch_j)
194 call xytoll(real(start_patch_i)-0.5, real(start_patch_j)-0.5, corner_lats(13), corner_lons(13), M)
195 call xytoll(real(start_patch_i)-0.5, real(end_patch_j)+0.5, corner_lats(14), corner_lons(14), M)
196 call xytoll(real(end_patch_i)+0.5, real(end_patch_j)+0.5, corner_lats(15), corner_lons(15), M)
197 call xytoll(real(end_patch_i)+0.5, real(start_patch_j)-0.5, corner_lats(16), corner_lons(16), M)
199 corner_lons(1) = xlon_array(start_patch_i,start_patch_j)
200 corner_lons(2) = xlon_array(start_patch_i,end_patch_j)
201 corner_lons(3) = xlon_array(end_patch_i,end_patch_j)
202 corner_lons(4) = xlon_array(end_patch_i,start_patch_j)
204 corner_lons(5) = xlon_array_u(start_patch_i,start_patch_j)
205 corner_lons(6) = xlon_array_u(start_patch_i,end_patch_j)
206 corner_lons(7) = xlon_array_u(end_patch_stag_i,end_patch_j)
207 corner_lons(8) = xlon_array_u(end_patch_stag_i,start_patch_j)
209 corner_lons(9) = xlon_array_v(start_patch_i,start_patch_j)
210 corner_lons(10) = xlon_array_v(start_patch_i,end_patch_stag_j)
211 corner_lons(11) = xlon_array_v(end_patch_i,end_patch_stag_j)
212 corner_lons(12) = xlon_array_v(end_patch_i,start_patch_j)
214 else if (grid_type == 'E') then
215 call get_lat_lon_fields(xlat_array, xlon_array, start_mem_i, &
216 start_mem_j, end_mem_i, end_mem_j, HH)
217 call get_lat_lon_fields(xlat_array_v, xlon_array_v, start_mem_i, &
218 start_mem_j, end_mem_i, end_mem_stag_j, VV)
220 corner_lats(1) = xlat_array(start_patch_i,start_patch_j)
221 corner_lats(2) = xlat_array(start_patch_i,end_patch_j)
222 corner_lats(3) = xlat_array(end_patch_i,end_patch_j)
223 corner_lats(4) = xlat_array(end_patch_i,start_patch_j)
225 corner_lats(5) = xlat_array_v(start_patch_i,start_patch_j)
226 corner_lats(6) = xlat_array_v(start_patch_i,end_patch_stag_j)
227 corner_lats(7) = xlat_array_v(end_patch_i,end_patch_stag_j)
228 corner_lats(8) = xlat_array_v(end_patch_i,start_patch_j)
230 corner_lats(9) = 0.0
231 corner_lats(10) = 0.0
232 corner_lats(11) = 0.0
233 corner_lats(12) = 0.0
235 corner_lats(13) = 0.0
236 corner_lats(14) = 0.0
237 corner_lats(15) = 0.0
238 corner_lats(16) = 0.0
240 corner_lons(1) = xlon_array(start_patch_i,start_patch_j)
241 corner_lons(2) = xlon_array(start_patch_i,end_patch_j)
242 corner_lons(3) = xlon_array(end_patch_i,end_patch_j)
243 corner_lons(4) = xlon_array(end_patch_i,start_patch_j)
245 corner_lons(5) = xlon_array_v(start_patch_i,start_patch_j)
246 corner_lons(6) = xlon_array_v(start_patch_i,end_patch_stag_j)
247 corner_lons(7) = xlon_array_v(end_patch_i,end_patch_stag_j)
248 corner_lons(8) = xlon_array_v(end_patch_i,start_patch_j)
250 corner_lons(9) = 0.0
251 corner_lons(10) = 0.0
252 corner_lons(11) = 0.0
253 corner_lons(12) = 0.0
255 corner_lons(13) = 0.0
256 corner_lons(14) = 0.0
257 corner_lons(15) = 0.0
258 corner_lons(16) = 0.0
260 end if
262 ! Initialize the output module now that we have the corner point lats/lons
263 call output_init(which_domain, 'OUTPUT FROM GEOGRID V3.3.1', '0000-00-00_00:00:00', grid_type, dynopt, &
264 corner_lats, corner_lons, &
265 start_dom_i, end_dom_i, start_dom_j, end_dom_j, &
266 start_patch_i, end_patch_i, start_patch_j, end_patch_j, &
267 start_mem_i, end_mem_i, start_mem_j, end_mem_j, &
268 extra_col, extra_row)
270 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, &
271 'XLAT_M', datestr, real_array = xlat_array)
272 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, &
273 'XLONG_M', datestr, real_array = xlon_array)
274 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_stag_j, 1, 1, &
275 'XLAT_V', datestr, real_array = xlat_array_v)
276 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_stag_j, 1, 1, &
277 'XLONG_V', datestr, real_array = xlon_array_v)
278 if (grid_type == 'C') then
279 call write_field(start_mem_i, end_mem_stag_i, start_mem_j, end_mem_j, 1, 1, &
280 'XLAT_U', datestr, real_array = xlat_array_u)
281 call write_field(start_mem_i, end_mem_stag_i, start_mem_j, end_mem_j, 1, 1, &
282 'XLONG_U', datestr, real_array = xlon_array_u)
283 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, &
284 'CLAT', datestr, real_array = clat_array)
285 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, &
286 'CLONG', datestr, real_array = clon_array)
288 if (associated(clat_array)) deallocate(clat_array)
289 if (associated(clon_array)) deallocate(clon_array)
291 end if
294 !
295 ! Calculate map factor for current domain
296 !
297 if (grid_type == 'C') then
298 call mprintf(.true.,STDOUT,' Processing MAPFAC')
300 allocate(mapfac_array_m_x(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
301 allocate(mapfac_array_m_y(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
302 call get_map_factor(xlat_array, xlon_array, mapfac_array_m_x, mapfac_array_m_y, start_mem_i, &
303 start_mem_j, end_mem_i, end_mem_j)
304 ! Global WRF uses map scale factors in X and Y directions, but "regular" WRF uses a single MSF
305 ! on each staggering. In the case of regular WRF, we can assume that MAPFAC_MX = MAPFAC_MY = MAPFAC_M,
306 ! and so we can simply write MAPFAC_MX as the MAPFAC_M field. Ultimately, when global WRF is
307 ! merged into the WRF trunk, we will need only two map scale factor fields for each staggering,
308 ! in the x and y directions, and these will be the same in the case of non-Cassini projections
309 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'MAPFAC_M', &
310 datestr, real_array = mapfac_array_m_x)
311 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'MAPFAC_MX', &
312 datestr, real_array = mapfac_array_m_x)
313 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'MAPFAC_MY', &
314 datestr, real_array = mapfac_array_m_y)
316 allocate(mapfac_array_v_x(start_mem_i:end_mem_i, start_mem_j:end_mem_stag_j))
317 allocate(mapfac_array_v_y(start_mem_i:end_mem_i, start_mem_j:end_mem_stag_j))
318 call get_map_factor(xlat_array_v, xlon_array_v, mapfac_array_v_x, mapfac_array_v_y, start_mem_i, &
319 start_mem_j, end_mem_i, end_mem_stag_j)
320 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_stag_j, 1, 1, 'MAPFAC_V', &
321 datestr, real_array = mapfac_array_v_x)
322 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_stag_j, 1, 1, 'MAPFAC_VX', &
323 datestr, real_array = mapfac_array_v_x)
324 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_stag_j, 1, 1, 'MAPFAC_VY', &
325 datestr, real_array = mapfac_array_v_y)
327 allocate(mapfac_array_u_x(start_mem_i:end_mem_stag_i, start_mem_j:end_mem_j))
328 allocate(mapfac_array_u_y(start_mem_i:end_mem_stag_i, start_mem_j:end_mem_j))
329 call get_map_factor(xlat_array_u, xlon_array_u, mapfac_array_u_x, mapfac_array_u_y, start_mem_i, &
330 start_mem_j, end_mem_stag_i, end_mem_j)
331 call write_field(start_mem_i, end_mem_stag_i, start_mem_j, end_mem_j, 1, 1, 'MAPFAC_U', &
332 datestr, real_array = mapfac_array_u_x)
333 call write_field(start_mem_i, end_mem_stag_i, start_mem_j, end_mem_j, 1, 1, 'MAPFAC_UX', &
334 datestr, real_array = mapfac_array_u_x)
335 call write_field(start_mem_i, end_mem_stag_i, start_mem_j, end_mem_j, 1, 1, 'MAPFAC_UY', &
336 datestr, real_array = mapfac_array_u_y)
338 end if
341 !
342 ! Coriolis parameters (E and F)
343 !
344 call mprintf(.true.,STDOUT,' Processing F and E')
346 allocate(f_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
347 allocate(e_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
349 call get_coriolis_parameters(xlat_array, f_array, e_array, &
350 start_mem_i, start_mem_j, end_mem_i, end_mem_j)
352 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'E', &
353 datestr, real_array = e_array)
354 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'F', &
355 datestr, real_array = f_array)
357 if (associated(f_array)) deallocate(f_array)
358 if (associated(e_array)) deallocate(e_array)
361 !
362 ! Rotation angle (SINALPHA and COSALPHA)
363 !
364 if (grid_type == 'C') then
365 call mprintf(.true.,STDOUT,' Processing ROTANG')
367 allocate(sina_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
368 allocate(cosa_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
370 call get_rotang(xlat_array, xlon_array, cosa_array, sina_array, &
371 start_mem_i, start_mem_j, end_mem_i, end_mem_j)
373 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'SINALPHA', &
374 datestr, real_array = sina_array)
375 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'COSALPHA', &
376 datestr, real_array = cosa_array)
378 if (associated(sina_array)) deallocate(sina_array)
379 if (associated(cosa_array)) deallocate(cosa_array)
380 end if
382 ! Every field up until now should probably just be processed regardless of what the user
383 ! has specified for fields to be processed.
384 ! Hereafter, we process user-specified fields
386 !
387 ! First process the field that we will derive a landmask from
388 !
389 call get_landmask_field(geog_data_res(which_domain), landmask_name, is_water_mask, landmask_value, istatus)
391 do kk=1,MAX_LANDMASK_CATEGORIES
392 if (landmask_value(kk) == INVALID) then
393 num_landmask_categories = kk-1
394 exit
395 end if
396 end do
397 if (kk > MAX_LANDMASK_CATEGORIES) num_landmask_categories = MAX_LANDMASK_CATEGORIES
399 if (istatus /= 0) then
400 call mprintf(.true.,WARN,'No field specified for landmask calculation. Will set landmask=1 at every grid point.')
402 allocate(landmask(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
403 landmask = 1.
404 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'LANDMASK', &
405 datestr, landmask)
407 else
409 allocate(landmask(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
410 landmask = 1.
412 call mprintf(.true.,STDOUT,' Processing %s', s1=trim(landmask_name))
414 call get_missing_fill_value(landmask_name, msg_fill_val, istatus)
415 if (istatus /= 0) msg_fill_val = NAN
417 call get_halt_on_missing(landmask_name, halt_on_missing, istatus)
418 if (istatus /= 0) halt_on_missing = .false.
420 ! Do we calculate a dominant category for this field?
421 call get_domcategory_name(landmask_name, domname, only_save_dominant, idomcatstatus)
423 temp_string = ' '
424 temp_string(1:128) = landmask_name
425 call hash_insert(processed_fieldnames, temp_string)
427 call get_max_categories(landmask_name, min_category, max_category, istatus)
428 allocate(field(start_mem_i:end_mem_i, start_mem_j:end_mem_j, min_category:max_category))
430 if (.not. only_save_dominant) then
431 field_count = field_count + 1
432 call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
433 i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=landmask_name)
434 else
435 field_count = field_count + 1
436 call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
437 i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=domname)
438 end if
440 if (grid_type == 'C') then
441 call calc_field(landmask_name, field, xlat_array, xlon_array, M, &
442 start_mem_i, end_mem_i, start_mem_j, end_mem_j, &
443 min_category, max_category, processed_domain, 1, landmask=landmask, sr_x=1, sr_y=1)
444 else if (grid_type == 'E') then
445 call calc_field(landmask_name, field, xlat_array, xlon_array, HH, &
446 start_mem_i, end_mem_i, start_mem_j, end_mem_j, &
447 min_category, max_category, processed_domain, 1, landmask=landmask, sr_x=1, sr_y=1)
448 end if
450 ! Find fractions
451 do i=start_mem_i, end_mem_i
452 do j=start_mem_j, end_mem_j
453 sum = 0.
454 do k=min_category,max_category
455 sum = sum + field(i,j,k)
456 end do
457 if (sum > 0.0) then
458 do k=min_category,max_category
459 field(i,j,k) = field(i,j,k) / sum
460 end do
461 else
462 do k=min_category,max_category
463 field(i,j,k) = msg_fill_val
464 end do
465 end if
466 end do
467 end do
469 ! If user wants to halt when a missing value is found in output field, check now
470 if (halt_on_missing) then
471 do i=start_mem_i, end_mem_i
472 do j=start_mem_j, end_mem_j
473 ! Only need to examine k=1
474 if (field(i,j,1) == msg_fill_val) then
475 call mprintf(.true.,ERROR,' Missing value encountered in output field. Quitting.')
476 end if
477 end do
478 end do
479 end if
481 if (is_water_mask) then
482 call mprintf(.true.,STDOUT,' Calculating landmask from %s ( WATER =', &
483 newline=.false.,s1=trim(landmask_name))
484 else
485 call mprintf(.true.,STDOUT,' Calculating landmask from %s ( LAND =', &
486 newline=.false.,s1=trim(landmask_name))
487 end if
488 do k = 1, num_landmask_categories
489 call mprintf(.true.,STDOUT,' %i',newline=.false.,i1=landmask_value(k))
490 if (k == num_landmask_categories) call mprintf(.true.,STDOUT,')')
491 end do
493 ! Calculate landmask
494 if (is_water_mask) then
495 do i=start_mem_i, end_mem_i
496 do j=start_mem_j, end_mem_j
497 water_total = -1.
498 do k=1,num_landmask_categories
499 if (landmask_value(k) >= min_category .and. landmask_value(k) <= max_category) then
500 if (field(i,j,landmask_value(k)) /= msg_fill_val) then
501 if (water_total < 0.) water_total = 0.
502 water_total = water_total + field(i,j,landmask_value(k))
503 else
504 water_total = -1.
505 exit
506 end if
507 end if
508 end do
509 if (water_total >= 0.0) then
510 if (water_total < 0.50) then
511 landmask(i,j) = 1.
512 else
513 landmask(i,j) = 0.
514 end if
515 else
516 landmask(i,j) = -1.
517 end if
518 end do
519 end do
520 else
521 do i=start_mem_i, end_mem_i
522 do j=start_mem_j, end_mem_j
523 land_total = -1.
524 do k=1,num_landmask_categories
525 if (landmask_value(k) >= min_category .and. landmask_value(k) <= max_category) then
526 if (field(i,j,landmask_value(k)) /= msg_fill_val) then
527 if (land_total < 0.) land_total = 0.
528 land_total = land_total + field(i,j,landmask_value(k))
529 else
530 land_total = -1.
531 exit
532 end if
533 end if
534 end do
535 if (land_total >= 0.0) then
536 if (land_total > 0.50) then
537 landmask(i,j) = 1.
538 else
539 landmask(i,j) = 0.
540 end if
541 else
542 landmask(i,j) = -1.
543 end if
544 end do
545 end do
546 end if
548 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'LANDMASK', &
549 datestr, landmask)
551 ! If we should only save the dominant category, then no need to write out fractional field
552 if (.not.only_save_dominant .or. (idomcatstatus /= 0)) then
554 ! Finally, we may be asked to smooth the fractional field
555 call get_smooth_option(landmask_name, smth_opt, smth_passes, istatus)
556 if (istatus == 0) then
558 if (grid_type == 'C') then
559 if (smth_opt == ONETWOONE) then
560 call one_two_one(field, &
561 start_patch_i, end_patch_i, &
562 start_patch_j, end_patch_j, &
563 start_mem_i, end_mem_i, &
564 start_mem_j, end_mem_j, &
565 min_category, max_category, &
566 smth_passes, msg_fill_val)
567 else if (smth_opt == SMTHDESMTH) then
568 call smth_desmth(field, &
569 start_patch_i, end_patch_i, &
570 start_patch_j, end_patch_j, &
571 start_mem_i, end_mem_i, &
572 start_mem_j, end_mem_j, &
573 min_category, max_category, &
574 smth_passes, msg_fill_val)
575 else if (smth_opt == SMTHDESMTH_SPECIAL) then
576 call smth_desmth_special(field, &
577 start_patch_i, end_patch_i, &
578 start_patch_j, end_patch_j, &
579 start_mem_i, end_mem_i, &
580 start_mem_j, end_mem_j, &
581 min_category, max_category, &
582 smth_passes, msg_fill_val)
583 end if
584 else if (grid_type == 'E') then
585 if (smth_opt == ONETWOONE) then
586 call one_two_one_egrid(field, &
587 start_patch_i, end_patch_i, &
588 start_patch_j, end_patch_j, &
589 start_mem_i, end_mem_i, &
590 start_mem_j, end_mem_j, &
591 min_category, max_category, &
592 smth_passes, msg_fill_val, 1.0)
593 else if (smth_opt == SMTHDESMTH) then
594 call smth_desmth_egrid(field, &
595 start_patch_i, end_patch_i, &
596 start_patch_j, end_patch_j, &
597 start_mem_i, end_mem_i, &
598 start_mem_j, end_mem_j, &
599 min_category, max_category, &
600 smth_passes, msg_fill_val, 1.0)
601 else if (smth_opt == SMTHDESMTH_SPECIAL) then
602 call mprintf(.true.,WARN,'smth-desmth_special is not currently implemented for NMM. '// &
603 'No smoothing will be done.')
604 end if
605 end if
607 end if
609 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, &
610 min_category, max_category, trim(landmask_name), &
611 datestr, real_array=field)
612 end if
614 if (idomcatstatus == 0) then
615 allocate(dominant_field(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
617 if (.not. only_save_dominant) then
618 field_count = field_count + 1
619 call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
620 i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=domname)
621 end if
623 do i=start_mem_i, end_mem_i
624 do j=start_mem_j, end_mem_j
625 if ((landmask(i,j) == 1. .and. is_water_mask) .or. &
626 (landmask(i,j) == 0. .and. .not.is_water_mask)) then
627 dominant = 0.
628 dominant_field(i,j) = real(min_category-1)
629 do k=min_category,max_category
630 do kk=1,num_landmask_categories
631 if (k == landmask_value(kk)) exit
632 end do
633 if (field(i,j,k) > dominant .and. kk > num_landmask_categories) then
634 dominant_field(i,j) = real(k)
635 dominant = field(i,j,k)
636 end if
637 end do
638 else
639 dominant = 0.
640 dominant_field(i,j) = real(min_category-1)
641 do k=min_category,max_category
642 do kk=1,num_landmask_categories
643 if (field(i,j,k) > dominant .and. k == landmask_value(kk)) then
644 dominant_field(i,j) = real(k)
645 dominant = field(i,j,k)
646 end if
647 end do
648 end do
649 end if
650 end do
651 end do
653 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, trim(domname), &
654 datestr, dominant_field)
656 deallocate(dominant_field)
657 end if
659 deallocate(field)
660 end if
662 !
663 ! Now process all other fields specified by the user
664 !
665 call reset_next_field()
666 ifieldstatus = 0
667 do while (ifieldstatus == 0)
668 call get_next_fieldname(fieldname, ifieldstatus)
670 ! There is another field in the GEOGRID.TBL file
671 if (ifieldstatus == 0) then
672 temp_string(1:128) = fieldname
674 ! If this field is still to be processed
675 if (.not. hash_search(processed_fieldnames, temp_string)) then
677 call hash_insert(processed_fieldnames, temp_string)
678 call mprintf(.true.,STDOUT,' Processing %s', s1=trim(fieldname))
680 call get_output_stagger(fieldname, istagger, istatus)
681 dimnames(:) = 'null'
682 call get_subgrid_dim_name(which_domain, fieldname, dimnames, &
683 sub_x, sub_y, istatus)
685 if (istagger == M .or. (sub_x > 1) .or. (sub_y > 1)) then
686 sm1 = start_mem_i
687 em1 = end_mem_i
688 sm2 = start_mem_j
689 em2 = end_mem_j
690 xlat_ptr => xlat_array
691 xlon_ptr => xlon_array
692 mapfac_ptr_x => mapfac_array_m_x
693 mapfac_ptr_y => mapfac_array_m_y
694 else if (istagger == U) then ! In the case that extra_cols = .false.
695 sm1 = start_mem_i ! we should have that end_mem_stag is
696 em1 = end_mem_stag_i ! the same as end_mem, so we do not need
697 sm2 = start_mem_j ! to check extra_cols or extra rows here
698 em2 = end_mem_j
699 xlat_ptr => xlat_array_u
700 xlon_ptr => xlon_array_u
701 mapfac_ptr_x => mapfac_array_u_x
702 mapfac_ptr_y => mapfac_array_u_y
703 else if (istagger == V) then
704 sm1 = start_mem_i
705 em1 = end_mem_i
706 sm2 = start_mem_j
707 em2 = end_mem_stag_j
708 xlat_ptr => xlat_array_v
709 xlon_ptr => xlon_array_v
710 mapfac_ptr_x => mapfac_array_v_x
711 mapfac_ptr_y => mapfac_array_v_y
712 else if (istagger == HH) then ! E grid
713 sm1 = start_mem_i
714 em1 = end_mem_i
715 sm2 = start_mem_j
716 em2 = end_mem_j
717 xlat_ptr => xlat_array
718 xlon_ptr => xlon_array
719 mapfac_ptr_x => mapfac_array_m_x
720 mapfac_ptr_y => mapfac_array_m_y
721 else if (istagger == VV) then ! E grid
722 sm1 = start_mem_i
723 em1 = end_mem_i
724 sm2 = start_mem_j
725 em2 = end_mem_stag_j
726 xlat_ptr => xlat_array_v
727 xlon_ptr => xlon_array_v
728 mapfac_ptr_x => mapfac_array_v_x
729 mapfac_ptr_y => mapfac_array_v_y
730 end if
732 if (sub_x > 1) then
733 sm1 = (start_mem_i - 1)*sub_x + 1
734 if (extra_col) then
735 em1 = (end_mem_i + 1)*sub_x
736 else
737 em1 = (end_mem_i )*sub_x
738 end if
739 end if
740 if (sub_y > 1)then
741 sm2 = (start_mem_j - 1)*sub_y + 1
742 if (extra_row) then
743 em2 = (end_mem_j + 1)*sub_y
744 else
745 em2 = (end_mem_j )*sub_y
746 end if
747 end if
749 !BUG: This should probably be moved up to where other lat/lon fields are calculated, and we should
750 ! just determine whether we will have any subgrids or not at that point
751 if ((sub_x > 1) .or. (sub_y > 1)) then
752 ! if (associated(xlat_array_subgrid)) deallocate(xlat_array_subgrid)
753 ! if (associated(xlon_array_subgrid)) deallocate(xlon_array_subgrid)
754 ! if (associated(mapfac_array_x_subgrid)) deallocate(mapfac_array_x_subgrid)
755 ! if (associated(mapfac_array_y_subgrid)) deallocate(mapfac_array_y_subgrid)
756 allocate(xlat_array_subgrid(sm1:em1,sm2:em2))
757 allocate(xlon_array_subgrid(sm1:em1,sm2:em2))
758 allocate(mapfac_array_x_subgrid(sm1:em1,sm2:em2))
759 allocate(mapfac_array_y_subgrid(sm1:em1,sm2:em2))
760 call get_lat_lon_fields(xlat_array_subgrid, xlon_array_subgrid, &
761 sm1, sm2, em1, em2, M, sub_x=sub_x, sub_y=sub_y)
762 xlat_ptr => xlat_array_subgrid
763 xlon_ptr => xlon_array_subgrid
764 call get_map_factor(xlat_ptr, xlon_ptr, mapfac_array_x_subgrid, &
765 mapfac_array_y_subgrid, sm1, sm2, em1, em2)
766 mapfac_ptr_x => mapfac_array_x_subgrid
767 mapfac_ptr_y => mapfac_array_y_subgrid
768 end if
770 call get_missing_fill_value(fieldname, msg_fill_val, istatus)
771 if (istatus /= 0) msg_fill_val = NAN
773 call get_halt_on_missing(fieldname, halt_on_missing, istatus)
774 if (istatus /= 0) halt_on_missing = .false.
776 ! Destination field type is CONTINUOUS
777 if (iget_fieldtype(fieldname,istatus) == CONTINUOUS) then
778 call get_max_levels(fieldname, min_level, max_level, istatus)
779 allocate(field(sm1:em1, sm2:em2, min_level:max_level))
781 field_count = field_count + 1
782 call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
783 i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=fieldname)
785 if ((sub_x > 1) .or. (sub_y > 1)) then
786 call calc_field(fieldname, field, xlat_ptr, xlon_ptr, istagger, &
787 sm1, em1, sm2, em2, min_level, max_level, &
788 processed_domain, 1, sr_x=sub_x, sr_y=sub_y)
789 else
790 call calc_field(fieldname, field, xlat_ptr, xlon_ptr, istagger, &
791 sm1, em1, sm2, em2, min_level, max_level, &
792 processed_domain, 1, landmask=landmask, sr_x=sub_x, sr_y=sub_y)
793 end if
795 ! If user wants to halt when a missing value is found in output field, check now
796 if (halt_on_missing) then
797 do i=sm1, em1
798 do j=sm2, em2
799 ! Only need to examine k=1
800 if (field(i,j,1) == msg_fill_val) then
801 call mprintf(.true.,ERROR,' Missing value encountered in output field. Quitting.')
802 end if
803 end do
804 end do
805 end if
807 ! We may be asked to smooth the fractional field
808 call get_smooth_option(fieldname, smth_opt, smth_passes, istatus)
809 if (istatus == 0) then
811 if (grid_type == 'C') then
812 if (smth_opt == ONETWOONE) then
813 call one_two_one(field, &
814 start_patch_i, end_patch_i, &
815 start_patch_j, end_patch_j, &
816 sm1, em1, &
817 sm2, em2, &
818 min_level, max_level, &
819 smth_passes, msg_fill_val)
820 else if (smth_opt == SMTHDESMTH) then
821 call smth_desmth(field, &
822 start_patch_i, end_patch_i, &
823 start_patch_j, end_patch_j, &
824 sm1, em1, &
825 sm2, em2, &
826 min_level, max_level, &
827 smth_passes, msg_fill_val)
828 else if (smth_opt == SMTHDESMTH_SPECIAL) then
829 call smth_desmth_special(field, &
830 start_patch_i, end_patch_i, &
831 start_patch_j, end_patch_j, &
832 sm1, em1, &
833 sm2, em2, &
834 min_level, max_level, &
835 smth_passes, msg_fill_val)
836 end if
838 else if (grid_type == 'E') then
840 if (trim(fieldname) == 'HGT_M' ) then
841 topo_flag_val=1.0
842 mass_flag=1.0
843 else if (trim(fieldname) == 'HGT_V') then
844 topo_flag_val=1.0
845 mass_flag=0.0
846 else
847 topo_flag_val=0.0
848 end if
850 if (smth_opt == ONETWOONE) then
851 call one_two_one_egrid(field, &
852 start_patch_i, end_patch_i, &
853 start_patch_j, end_patch_j, &
854 sm1, em1, &
855 sm2, em2, &
856 min_level, max_level, &
857 smth_passes, topo_flag_val, mass_flag)
858 else if (smth_opt == SMTHDESMTH) then
859 call smth_desmth_egrid(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, topo_flag_val, mass_flag)
866 else if (smth_opt == SMTHDESMTH_SPECIAL) then
867 call mprintf(.true.,WARN,'smth-desmth_special is not currently implemented for NMM. '// &
868 'No smoothing will be done.')
869 end if
871 end if
873 end if
875 call write_field(sm1, em1, sm2, em2, &
876 min_level, max_level, trim(fieldname), datestr, real_array=field)
878 ! Do we calculate directional derivatives from this field?
879 call get_dfdx_name(fieldname, gradname, istatus)
880 if (istatus == 0) then
881 allocate(slp_field(sm1:em1,sm2:em2,min_level:max_level))
883 field_count = field_count + 1
884 call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
885 i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=gradname)
887 if (grid_type == 'C') then
888 call calc_dfdx(field, slp_field, sm1, sm2, min_level, em1, em2, max_level, mapfac_ptr_x)
889 else if (grid_type == 'E') then
890 call calc_dfdx(field, slp_field, sm1, sm2, min_level, em1, em2, max_level)
891 end if
892 call write_field(sm1, em1, sm2, em2, &
893 min_level, max_level, trim(gradname), datestr, real_array=slp_field)
894 deallocate(slp_field)
895 end if
896 call get_dfdy_name(fieldname, gradname, istatus)
897 if (istatus == 0) then
898 allocate(slp_field(sm1:em1,sm2:em2,min_level:max_level))
900 field_count = field_count + 1
901 call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
902 i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=gradname)
904 if (grid_type == 'C') then
905 call calc_dfdy(field, slp_field, sm1, sm2, min_level, em1, em2, max_level, mapfac_ptr_y)
906 else if (grid_type == 'E') then
907 call calc_dfdy(field, slp_field, sm1, sm2, min_level, em1, em2, max_level)
908 end if
909 call write_field(sm1, em1, sm2, em2, &
910 min_level, max_level, trim(gradname), datestr, real_array=slp_field)
911 deallocate(slp_field)
912 end if
914 deallocate(field)
916 ! Destination field type is CATEGORICAL
917 else
918 call get_max_categories(fieldname, min_category, max_category, istatus)
919 allocate(field(sm1:em1, sm2:em2, min_category:max_category))
921 ! Do we calculate a dominant category for this field?
922 call get_domcategory_name(fieldname, domname, only_save_dominant, idomcatstatus)
924 if (.not. only_save_dominant) then
925 field_count = field_count + 1
926 call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
927 i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=fieldname)
928 else
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=domname)
932 end if
934 if ((sub_x > 1) .or. (sub_y > 1)) then
935 call calc_field(fieldname, field, xlat_ptr, xlon_ptr, istagger, &
936 sm1, em1, sm2, em2, min_category, max_category, &
937 processed_domain, 1, sr_x=sub_x, sr_y=sub_y)
938 else
939 call calc_field(fieldname, field, xlat_ptr, xlon_ptr, istagger, &
940 sm1, em1, sm2, em2, min_category, max_category, &
941 processed_domain, 1, landmask=landmask, sr_x=sub_x, sr_y=sub_y)
942 end if
944 ! Find fractions
945 do i=sm1, em1
946 do j=sm2, em2
947 sum = 0.
948 do k=min_category,max_category
949 sum = sum + field(i,j,k)
950 end do
951 if (sum > 0.0) then
952 do k=min_category,max_category
953 field(i,j,k) = field(i,j,k) / sum
954 end do
955 else
956 do k=min_category,max_category
957 field(i,j,k) = msg_fill_val
958 end do
959 end if
960 end do
961 end do
963 ! If user wants to halt when a missing value is found in output field, check now
964 if (halt_on_missing) then
965 do i=sm1, em1
966 do j=sm2, em2
967 ! Only need to examine k=1
968 if (field(i,j,1) == msg_fill_val) then
969 call mprintf(.true.,ERROR,' Missing value encountered in output field. Quitting.')
970 end if
971 end do
972 end do
973 end if
975 ! If we should only save the dominant category, then no need to write out fractional field
976 if (.not.only_save_dominant .or. (idomcatstatus /= 0)) then
978 ! Finally, we may be asked to smooth the fractional field
979 call get_smooth_option(fieldname, smth_opt, smth_passes, istatus)
980 if (istatus == 0) then
981 if (grid_type == 'C') then
982 if (smth_opt == ONETWOONE) then
983 call one_two_one(field, &
984 start_patch_i, end_patch_i, &
985 start_patch_j, end_patch_j, &
986 sm1, em1, &
987 sm2, em2, &
988 min_category, max_category, &
989 smth_passes, msg_fill_val)
990 else if (smth_opt == SMTHDESMTH) then
991 call smth_desmth(field, &
992 start_patch_i, end_patch_i, &
993 start_patch_j, end_patch_j, &
994 sm1, em1, &
995 sm2, em2, &
996 min_category, max_category, &
997 smth_passes, msg_fill_val)
998 else if (smth_opt == SMTHDESMTH_SPECIAL) then
999 call smth_desmth_special(field, &
1000 start_patch_i, end_patch_i, &
1001 start_patch_j, end_patch_j, &
1002 sm1, em1, &
1003 sm2, em2, &
1004 min_category, max_category, &
1005 smth_passes, msg_fill_val)
1006 end if
1007 else if (grid_type == 'E') then
1008 if (smth_opt == ONETWOONE) then
1009 call one_two_one_egrid(field, &
1010 start_patch_i, end_patch_i, &
1011 start_patch_j, end_patch_j, &
1012 sm1, em1, &
1013 sm2, em2, &
1014 min_category, max_category, &
1015 smth_passes, msg_fill_val, 1.0)
1016 else if (smth_opt == SMTHDESMTH) then
1017 call smth_desmth_egrid(field, &
1018 start_patch_i, end_patch_i, &
1019 start_patch_j, end_patch_j, &
1020 sm1, em1, &
1021 sm2, em2, &
1022 min_category, max_category, &
1023 smth_passes, msg_fill_val, 1.0)
1024 else if (smth_opt == SMTHDESMTH_SPECIAL) then
1025 call mprintf(.true.,WARN,'smth-desmth_special is not currently implemented for NMM. '// &
1026 'No smoothing will be done.')
1027 end if
1028 end if
1029 end if
1031 call write_field(sm1, em1, sm2, em2, &
1032 min_category, max_category, trim(fieldname), datestr, real_array=field)
1033 end if
1035 if (idomcatstatus == 0) then
1036 call mprintf(.true.,STDOUT,' Processing %s', s1=trim(domname))
1037 allocate(dominant_field(sm1:em1, sm2:em2))
1039 if (.not. only_save_dominant) then
1040 field_count = field_count + 1
1041 call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
1042 i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=domname)
1043 end if
1045 do i=sm1, em1
1046 do j=sm2, em2
1047 dominant = 0.
1048 dominant_field(i,j) = real(min_category-1)
1049 do k=min_category,max_category
1050 if (field(i,j,k) > dominant .and. field(i,j,k) /= msg_fill_val) then
1051 dominant_field(i,j) = real(k)
1052 dominant = field(i,j,k)
1053 ! else
1054 ! dominant_field(i,j) = nint(msg_fill_val)
1055 ! Maybe we should put an else clause here to set the category equal to the missing fill value?
1056 ! BUG: The problem here seems to be that, when we set a fraction equal to the missing fill value
1057 ! above, if the last fractional index we process here has been filled, we think that the dominant
1058 ! category should be set to the missing fill value. Perhaps we could do some check to only
1059 ! assign the msg_fill_val if no other valid category has been assigned? But this may still not
1060 ! work if the missing fill value is something like 0.5. Somehow use bitarrays, perhaps, to remember
1061 ! which points are missing and which just happen to have the missing fill value?
1062 end if
1063 end do
1064 if (dominant_field(i,j) == real(min_category-1)) dominant_field(i,j) = msg_fill_val
1065 end do
1066 end do
1067 call write_field(sm1, em1, sm2, em2, 1, 1, &
1068 trim(domname), datestr, dominant_field)
1069 deallocate(dominant_field)
1070 end if
1072 deallocate(field)
1074 if ((sub_x > 1) .or. (sub_y > 1)) then
1075 if (associated(xlat_array_subgrid)) deallocate(xlat_array_subgrid)
1076 if (associated(xlon_array_subgrid)) deallocate(xlon_array_subgrid)
1077 if (associated(mapfac_array_x_subgrid)) deallocate(mapfac_array_x_subgrid)
1078 if (associated(mapfac_array_y_subgrid)) deallocate(mapfac_array_y_subgrid)
1079 end if
1081 end if
1083 end if
1084 end if
1086 end do
1088 ! Close output
1089 call output_close()
1091 call hash_destroy(processed_fieldnames)
1093 ! Free up memory
1094 if (associated(xlat_array)) deallocate(xlat_array)
1095 if (associated(xlon_array)) deallocate(xlon_array)
1096 if (grid_type == 'C') then
1097 if (associated(xlat_array_u)) deallocate(xlat_array_u)
1098 if (associated(xlon_array_u)) deallocate(xlon_array_u)
1099 if (associated(mapfac_array_u_x)) deallocate(mapfac_array_u_x)
1100 if (associated(mapfac_array_u_y)) deallocate(mapfac_array_u_y)
1101 end if
1102 if (associated(xlat_array_v)) deallocate(xlat_array_v)
1103 if (associated(xlon_array_v)) deallocate(xlon_array_v)
1104 if (associated(mapfac_array_m_x)) deallocate(mapfac_array_m_x)
1105 if (associated(mapfac_array_m_y)) deallocate(mapfac_array_m_y)
1106 if (associated(mapfac_array_v_x)) deallocate(mapfac_array_v_x)
1107 if (associated(mapfac_array_v_y)) deallocate(mapfac_array_v_y)
1108 if (associated(landmask)) deallocate(landmask)
1109 if (associated(xlat_array_subgrid)) deallocate(xlat_array_subgrid)
1110 if (associated(xlon_array_subgrid)) deallocate(xlon_array_subgrid)
1111 if (associated(mapfac_array_x_subgrid)) deallocate(mapfac_array_x_subgrid)
1112 if (associated(mapfac_array_y_subgrid)) deallocate(mapfac_array_y_subgrid)
1114 nullify(xlat_ptr)
1115 nullify(xlon_ptr)
1117 end subroutine process_tile
1120 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1121 ! Name: calc_field
1122 !
1123 ! Purpose: This routine fills in the "field" array with interpolated source
1124 ! data. When multiple resolutions of source data are available, an appropriate
1125 ! resolution is chosen automatically. The specified field may either be a
1126 ! continuous field or a categorical field.
1127 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1128 recursive subroutine calc_field(fieldname, field, xlat_array, xlon_array, istagger, &
1129 start_i, end_i, start_j, end_j, start_k, end_k, &
1130 processed_domain, ilevel, landmask, sr_x, sr_y)
1132 use bitarray_module
1133 use interp_module
1134 use llxy_module
1135 use misc_definitions_module
1136 use proc_point_module
1137 use queue_module
1138 use source_data_module
1140 implicit none
1142 ! Arguments
1143 integer, intent(in) :: start_i, end_i, start_j, end_j, start_k, end_k, ilevel, istagger
1144 real, dimension(start_i:end_i, start_j:end_j), intent(in) :: xlat_array, xlon_array
1145 real, dimension(start_i:end_i, start_j:end_j, start_k:end_k), intent(inout) :: field
1146 real, dimension(start_i:end_i, start_j:end_j), intent(in), optional :: landmask
1147 integer, intent(in), optional :: sr_x, sr_y
1148 character (len=128), intent(in) :: fieldname
1149 type (bitarray), intent(inout) :: processed_domain
1151 ! Local variables
1152 integer :: start_src_k, end_src_k
1153 integer :: i, j, k, ix, iy, itype
1154 integer :: user_iproj, istatus
1155 real :: mask_val
1156 real :: temp
1157 real :: scale_factor
1158 real :: msg_val, msg_fill_val, threshold, src_dx, src_dy, dom_dx, dom_dy
1159 real :: user_stand_lon, user_truelat1, user_truelat2, user_dxkm, user_dykm, &
1160 user_known_x, user_known_y, user_known_lat, user_known_lon
1161 real, pointer, dimension(:,:,:) :: data_count
1162 integer, pointer, dimension(:) :: interp_type
1163 character (len=128) :: interp_string
1164 type (bitarray) :: bit_domain, level_domain
1165 type (queue) :: point_queue, tile_queue
1166 type (q_data) :: current_pt
1168 ! If this is the first trip through this routine, we need to allocate the bit array that
1169 ! will persist through all recursive calls, tracking which grid points have been assigned
1170 ! a value.
1171 if (ilevel == 1) call bitarray_create(processed_domain, end_i-start_i+1, end_j-start_j+1)
1173 ! Find out the projection of the data for this "priority level" (given by ilevel)
1174 call get_data_projection(fieldname, user_iproj, user_stand_lon, user_truelat1, user_truelat2, &
1175 user_dxkm, user_dykm, user_known_x, user_known_y, user_known_lat, &
1176 user_known_lon, ilevel, istatus)
1178 ! A bad status indicates that that no data for priority level ilevel is available, and thus, that
1179 ! no further levels will be specified. We are done processing for this level.
1180 if (istatus /= 0) then
1181 if (ilevel == 1) call bitarray_destroy(processed_domain)
1182 return
1183 end if
1185 ! A good status indicates that there is data for this priority level, so we store the projection
1186 ! of that data on a stack. The projection will be on the top of the stack (and hence will be
1187 ! the "active" projection) once all higher-priority levels have been processed
1188 call push_source_projection(user_iproj, user_stand_lon, user_truelat1, user_truelat2, &
1189 user_dxkm, user_dykm, user_dykm, user_dxkm, user_known_x, user_known_y, &
1190 user_known_lat, user_known_lon)
1192 ! Before proceeding with processing for this level, though, process for the next highest priority level
1193 ! of source data
1194 call calc_field(fieldname, field, xlat_array, xlon_array, istagger, start_i, end_i, &
1195 start_j, end_j, start_k, end_k, processed_domain, ilevel+1, landmask, sr_x, sr_y)
1197 ! At this point, all levels of source data with higher priority have been processed, and we can assign
1198 ! values to all grid points that have not already been given values from higher-priority data
1200 ! Initialize point processing module
1201 call proc_point_init()
1203 ! Initialize queues
1204 call q_init(point_queue)
1205 call q_init(tile_queue)
1207 ! Determine whether we will be processing categorical data or continuous data
1208 itype = iget_source_fieldtype(fieldname, ilevel, istatus)
1209 call get_interp_option(fieldname, ilevel, interp_string, istatus)
1210 interp_type => interp_array_from_string(interp_string)
1212 ! Also, check whether we will be using the cell averaging interpolator for continuous fields
1213 if (index(interp_string,'average_gcell') /= 0 .and. itype == CONTINUOUS) then
1214 call get_gcell_threshold(interp_string, threshold, istatus)
1215 if (istatus == 0) then
1216 call get_source_resolution(fieldname, ilevel, src_dx, src_dy, istatus)
1217 if (istatus == 0) then
1218 call get_domain_resolution(dom_dx, dom_dy)
1219 if (gridtype == 'C') then
1220 if (threshold*max(src_dx,src_dy)*111. <= max(dom_dx,dom_dy)/1000.) then
1221 itype = SP_CONTINUOUS
1222 allocate(data_count(start_i:end_i,start_j:end_j,start_k:end_k))
1223 data_count = 0.
1224 end if
1225 else if (gridtype == 'E') then
1226 if (max(src_dx,src_dy) >= threshold*max(dom_dx,dom_dy)) then
1227 itype = SP_CONTINUOUS
1228 allocate(data_count(start_i:end_i,start_j:end_j,start_k:end_k))
1229 data_count = 0.
1230 end if
1231 end if
1232 end if
1233 end if
1234 end if
1236 call get_missing_value(fieldname, ilevel, msg_val, istatus)
1237 if (istatus /= 0) msg_val = NAN
1238 call get_missing_fill_value(fieldname, msg_fill_val, istatus)
1239 if (istatus /= 0) msg_fill_val = NAN
1241 call get_masked_value(fieldname, ilevel, mask_val, istatus)
1242 if (istatus /= 0) mask_val = -1.
1244 if (itype == CONTINUOUS .or. itype == SP_CONTINUOUS) then
1245 call get_source_levels(fieldname, ilevel, start_src_k, end_src_k, istatus)
1246 if (istatus /= 0) return
1247 end if
1249 ! Initialize bitarray used to track which points have been visited and assigned values while
1250 ! processing *this* priority level of data
1251 call bitarray_create(bit_domain, end_i-start_i+1, end_j-start_j+1)
1252 call bitarray_create(level_domain, end_i-start_i+1, end_j-start_j+1)
1254 ! Begin by placing a point in the tile_queue
1255 current_pt%lat = xlat_array(start_i,start_j)
1256 current_pt%lon = xlon_array(start_i,start_j)
1257 current_pt%x = start_i
1258 current_pt%y = start_j
1259 call q_insert(tile_queue, current_pt)
1261 ! While there are still grid points in tiles that have not yet been processed
1262 do while (q_isdata(tile_queue))
1264 ! Take a point from the outer queue and place it in the point_queue for processing
1265 current_pt = q_remove(tile_queue)
1267 ! If this level of data is categorical (i.e., is given as an array of category indices),
1268 ! then first try to process the entire tile in one call to accum_categorical. Any grid
1269 ! points that are not given values by accum_categorical and that lie within the current
1270 ! tile of source data are individually assigned values in the inner loop
1271 if (itype == CATEGORICAL) then
1273 ! Have we already visited this point? If so, this tile has already been processed by
1274 ! accum_categorical.
1275 if (.not. bitarray_test(bit_domain, current_pt%x-start_i+1, current_pt%y-start_j+1)) then
1276 call q_insert(point_queue, current_pt)
1277 if (.not. have_processed_tile(current_pt%lat, current_pt%lon, fieldname, ilevel)) then
1278 call accum_categorical(current_pt%lat, current_pt%lon, istagger, field, &
1279 start_i, end_i, start_j, end_j, start_k, end_k, &
1280 fieldname, processed_domain, level_domain, &
1281 ilevel, msg_val, mask_val, sr_x, sr_y)
1282 ! BUG: Where do we mask out those points that are on land/water when masked=land/water is set?
1283 end if
1284 call bitarray_set(bit_domain, current_pt%x-start_i+1, current_pt%y-start_j+1)
1285 end if
1287 else if (itype == SP_CONTINUOUS) then
1289 ! Have we already visited this point? If so, this tile has already been processed by
1290 ! accum_continuous.
1291 if (.not. bitarray_test(bit_domain, current_pt%x-start_i+1, current_pt%y-start_j+1)) then
1292 call q_insert(point_queue, current_pt)
1293 if (.not. have_processed_tile(current_pt%lat, current_pt%lon, fieldname, ilevel)) then
1294 call accum_continuous(current_pt%lat, current_pt%lon, istagger, field, data_count, &
1295 start_i, end_i, start_j, end_j, start_k, end_k, &
1296 fieldname, processed_domain, level_domain, &
1297 ilevel, msg_val, mask_val, sr_x, sr_y)
1298 ! BUG: Where do we mask out those points that are on land/water when masked=land/water is set?
1299 end if
1300 call bitarray_set(bit_domain, current_pt%x-start_i+1, current_pt%y-start_j+1)
1301 end if
1303 else if (itype == CONTINUOUS) then
1305 ! Have we already visited this point? If so, the tile containing this point has already been
1306 ! processed in the inner loop.
1307 if (.not. bitarray_test(bit_domain, current_pt%x-start_i+1, current_pt%y-start_j+1)) then
1308 call q_insert(point_queue, current_pt)
1309 call bitarray_set(bit_domain, current_pt%x-start_i+1, current_pt%y-start_j+1)
1310 end if
1312 end if
1314 ! This inner loop, where all grid points contained in the current source tile are processed
1315 do while (q_isdata(point_queue))
1316 current_pt = q_remove(point_queue)
1317 ix = current_pt%x
1318 iy = current_pt%y
1320 ! Process the current point
1321 if (itype == CONTINUOUS .or. itype == SP_CONTINUOUS) then
1323 ! Have we already assigned this point a value from this priority level?
1324 if (.not. bitarray_test(level_domain, ix-start_i+1, iy-start_j+1)) then
1326 ! If the point was already assigned a value from a higher-priority level, no
1327 ! need to assign a new value
1328 if (bitarray_test(processed_domain, ix-start_i+1, iy-start_j+1)) then
1329 call bitarray_set(level_domain, ix-start_i+1, iy-start_j+1)
1331 ! Otherwise, need to assign values from this level of source data if we can
1332 else
1333 if (present(landmask) .and. (istagger == M .or. istagger == HH)) then
1334 if (landmask(ix,iy) /= mask_val) then
1335 do k=start_src_k,end_src_k
1336 temp = get_point(current_pt%lat, current_pt%lon, k, &
1337 fieldname, ilevel, interp_type, msg_val)
1338 if (temp /= msg_val) then
1339 field(ix, iy, k) = temp
1340 call bitarray_set(level_domain, ix-start_i+1, iy-start_j+1)
1341 if (itype == SP_CONTINUOUS) data_count(ix, iy, k) = 1.0
1342 else
1343 field(ix, iy, k) = msg_fill_val
1344 end if
1345 end do
1346 else
1347 do k=start_k,end_k
1348 field(ix,iy,k) = msg_fill_val
1349 end do
1350 end if
1351 else
1352 do k=start_src_k,end_src_k
1353 temp = get_point(current_pt%lat, current_pt%lon, k, &
1354 fieldname, ilevel, interp_type, msg_val)
1355 if (temp /= msg_val) then
1356 field(ix, iy, k) = temp
1357 call bitarray_set(level_domain, ix-start_i+1, iy-start_j+1)
1358 if (itype == SP_CONTINUOUS) data_count(ix, iy, k) = 1.0
1359 else
1360 field(ix, iy, k) = msg_fill_val
1361 end if
1362 end do
1363 end if
1364 end if
1365 end if
1367 else if (itype == CATEGORICAL) then
1369 ! Have we already assigned this point a value from this priority level?
1370 if (.not.bitarray_test(level_domain, ix-start_i+1, iy-start_j+1)) then
1372 ! If the point was already assigned a value from a higher-priority level, no
1373 ! need to assign a new value
1374 if (bitarray_test(processed_domain, ix-start_i+1, iy-start_j+1)) then
1375 call bitarray_set(level_domain, ix-start_i+1, iy-start_j+1)
1377 ! Otherwise, the point was apparently not given a value when accum_categorical
1378 ! was called for the current tile, and we need to assign values from this
1379 ! 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 temp = get_point(current_pt%lat, current_pt%lon, 1, &
1384 fieldname, ilevel, interp_type, msg_val)
1386 do k=start_k,end_k
1387 field(ix,iy,k) = 0.
1388 end do
1390 if (temp /= msg_val) then
1391 if (int(temp) >= start_k .and. int(temp) <= end_k) then
1392 field(ix, iy, int(temp)) = field(ix, iy, int(temp)) + 1.
1393 else
1394 call mprintf(.true.,WARN,' Attempted to assign an invalid category '// &
1395 '%i to grid point (%i, %i)', i1=int(temp), i2=ix, i3=iy)
1396 end if
1397 call bitarray_set(level_domain, ix-start_i+1, iy-start_j+1)
1398 end if
1400 else
1401 do k=start_k,end_k
1402 field(ix,iy,k) = 0.
1403 end do
1404 end if
1405 else
1406 temp = get_point(current_pt%lat, current_pt%lon, 1, &
1407 fieldname, ilevel, interp_type, msg_val)
1409 do k=start_k,end_k
1410 field(ix,iy,k) = 0.
1411 end do
1413 if (temp /= msg_val) then
1414 if (int(temp) >= start_k .and. int(temp) <= end_k) then
1415 field(ix, iy, int(temp)) = field(ix, iy, int(temp)) + 1.
1416 else
1417 call mprintf(.true.,WARN,' Attempted to assign an invalid category '// &
1418 '%i to grid point (%i, %i)', i1=int(temp), i2=ix, i3=iy)
1419 end if
1420 call bitarray_set(level_domain, ix-start_i+1, iy-start_j+1)
1421 end if
1422 end if
1423 end if
1424 end if
1426 end if
1428 ! Scan neighboring points, adding them to the appropriate queue based on whether they
1429 ! are in the current tile or not
1430 if (iy > start_j) then
1431 if (ix > start_i) then
1433 ! Neighbor with relative position (-1,-1)
1434 call process_neighbor(ix-1, iy-1, bit_domain, point_queue, tile_queue, &
1435 xlat_array, xlon_array, start_i, end_i, start_j, end_j, ilevel)
1436 end if
1438 ! Neighbor with relative position (0,-1)
1439 call process_neighbor(ix, iy-1, bit_domain, point_queue, tile_queue, &
1440 xlat_array, xlon_array, start_i, end_i, start_j, end_j, ilevel)
1442 if (ix < end_i) then
1444 ! Neighbor with relative position (+1,-1)
1445 call process_neighbor(ix+1, iy-1, bit_domain, point_queue, tile_queue, &
1446 xlat_array, xlon_array, start_i, end_i, start_j, end_j, ilevel)
1447 end if
1448 end if
1450 if (ix > start_i) then
1452 ! Neighbor with relative position (-1,0)
1453 call process_neighbor(ix-1, iy, bit_domain, point_queue, tile_queue, &
1454 xlat_array, xlon_array, start_i, end_i, start_j, end_j, ilevel)
1455 end if
1457 if (ix < end_i) then
1459 ! Neighbor with relative position (+1,0)
1460 call process_neighbor(ix+1, iy, bit_domain, point_queue, tile_queue, &
1461 xlat_array, xlon_array, start_i, end_i, start_j, end_j, ilevel)
1462 end if
1464 if (iy < end_j) then
1465 if (ix > start_i) then
1467 ! Neighbor with relative position (-1,+1)
1468 call process_neighbor(ix-1, iy+1, bit_domain, point_queue, tile_queue, &
1469 xlat_array, xlon_array, start_i, end_i, start_j, end_j, ilevel)
1470 end if
1472 ! Neighbor with relative position (0,+1)
1473 call process_neighbor(ix, iy+1, bit_domain, point_queue, tile_queue, &
1474 xlat_array, xlon_array, start_i, end_i, start_j, end_j, ilevel)
1475 if (ix < end_i) then
1477 ! Neighbor with relative position (+1,+1)
1478 call process_neighbor(ix+1, iy+1, bit_domain, point_queue, tile_queue, &
1479 xlat_array, xlon_array, start_i, end_i, start_j, end_j, ilevel)
1480 end if
1481 end if
1483 end do
1485 end do
1487 if (itype == SP_CONTINUOUS) then
1488 itype = CONTINUOUS
1489 if (present(landmask) .and. (istagger == M .or. istagger == HH)) then
1490 do j=start_j,end_j
1491 do i=start_i,end_i
1492 if (landmask(i,j) /= mask_val) then
1493 do k=start_k,end_k
1494 if (data_count(i,j,k) > 0.) then
1495 field(i,j,k) = field(i,j,k) / data_count(i,j,k)
1496 else
1497 if (.not.bitarray_test(processed_domain, i-start_i+1, j-start_j+1)) then
1498 field(i,j,k) = msg_fill_val
1499 end if
1500 end if
1501 end do
1502 else
1503 if (.not.bitarray_test(processed_domain, i-start_i+1, j-start_j+1)) then
1504 do k=start_k,end_k
1505 field(i,j,k) = msg_fill_val
1506 end do
1507 end if
1508 end if
1509 end do
1510 end do
1511 else
1512 do k=start_k,end_k
1513 do j=start_j,end_j
1514 do i=start_i,end_i
1515 if (data_count(i,j,k) > 0.) then
1516 field(i,j,k) = field(i,j,k) / data_count(i,j,k)
1517 else
1518 if (.not.bitarray_test(processed_domain, i-start_i+1, j-start_j+1)) then
1519 field(i,j,k) = msg_fill_val
1520 end if
1521 end if
1522 end do
1523 end do
1524 end do
1525 end if
1526 deallocate(data_count)
1528 else if (itype == CATEGORICAL) then
1529 if (present(landmask) .and. (istagger == M .or. istagger == HH)) then
1530 do j=start_j,end_j
1531 do i=start_i,end_i
1532 if (landmask(i,j) == mask_val) then
1533 do k=start_k,end_k
1534 field(i,j,k) = 0.
1535 end do
1536 end if
1537 end do
1538 end do
1539 end if
1540 end if
1542 deallocate(interp_type)
1545 ! We may need to scale this field by a constant
1546 call get_field_scale_factor(fieldname, ilevel, scale_factor, istatus)
1547 if (istatus == 0) then
1548 do i=start_i, end_i
1549 do j=start_j, end_j
1550 if (bitarray_test(level_domain,i-start_i+1,j-start_j+1) .and. &
1551 .not. bitarray_test(processed_domain,i-start_i+1,j-start_j+1)) then
1552 do k=start_k,end_k
1553 if (field(i,j,k) /= msg_fill_val) then
1554 field(i,j,k) = field(i,j,k) * scale_factor
1555 end if
1556 end do
1557 end if
1558 end do
1559 end do
1560 end if
1563 ! Now add the points that were assigned values at this priority level to the complete array
1564 ! of points that have been assigned values
1565 call bitarray_merge(processed_domain, level_domain)
1567 call bitarray_destroy(bit_domain)
1568 call bitarray_destroy(level_domain)
1569 call q_destroy(point_queue)
1570 call q_destroy(tile_queue)
1571 call proc_point_shutdown()
1572 ! If this is the last level of the recursion, we can also deallocate processed_domain
1573 if (ilevel == 1) call bitarray_destroy(processed_domain)
1575 ! Remove the projection of the current level of source data from the stack, thus "activating"
1576 ! the projection of the next highest level
1577 call pop_source_projection()
1579 end subroutine calc_field
1582 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1583 ! Name: get_lat_lon_fields
1584 !
1585 ! Purpose: To calculate the latitude and longitude for every gridpoint in the
1586 ! tile of the model domain. The caller may specify that the grid for which
1587 ! values are computed is staggered or unstaggered using the "stagger"
1588 ! argument.
1589 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1590 subroutine get_lat_lon_fields(xlat_arr, xlon_arr, start_mem_i, &
1591 start_mem_j, end_mem_i, end_mem_j, stagger, comp_ll, &
1592 sub_x, sub_y)
1594 use llxy_module
1595 use misc_definitions_module
1597 implicit none
1599 ! Arguments
1600 integer, intent(in) :: start_mem_i, start_mem_j, end_mem_i, &
1601 end_mem_j, stagger
1602 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(out) :: xlat_arr, xlon_arr
1603 logical, optional, intent(in) :: comp_ll
1604 integer, optional, intent(in) :: sub_x, sub_y
1606 ! Local variables
1607 integer :: i, j
1608 real :: rx, ry
1610 rx = 1.0
1611 ry = 1.0
1612 if (present(sub_x)) rx = real(sub_x)
1613 if (present(sub_y)) ry = real(sub_y)
1615 do i=start_mem_i, end_mem_i
1616 do j=start_mem_j, end_mem_j
1617 call xytoll(real(i-1)/rx+1.0, real(j-1)/ry+1.0, &
1618 xlat_arr(i,j), xlon_arr(i,j), stagger, comp_ll=comp_ll)
1619 end do
1620 end do
1622 end subroutine get_lat_lon_fields
1625 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1626 ! Name: get_map_factor
1627 !
1628 ! Purpose: Given the latitude field, this routine calculates map factors for
1629 ! the grid points of the specified domain. For different grids (e.g., C grid,
1630 ! E grid), the latitude array should provide the latitudes of the points for
1631 ! which map factors are to be calculated.
1632 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1633 subroutine get_map_factor(xlat_arr, xlon_arr, mapfac_arr_x, mapfac_arr_y, &
1634 start_mem_i, start_mem_j, end_mem_i, end_mem_j)
1636 use constants_module
1637 use gridinfo_module
1638 use misc_definitions_module
1639 use map_utils
1641 implicit none
1643 ! Arguments
1644 integer, intent(in) :: start_mem_i, start_mem_j, end_mem_i, end_mem_j
1645 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(in) :: xlat_arr, xlon_arr
1646 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(out) :: mapfac_arr_x
1647 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(out) :: mapfac_arr_y
1649 ! Local variables
1650 integer :: i, j
1651 real :: n, colat, colat0, colat1, colat2, comp_lat, comp_lon
1653 !
1654 ! Equations for map factor given in Principles of Meteorological Analysis,
1655 ! Walter J. Saucier, pp. 32-33
1656 !
1658 ! Lambert conformal projection
1659 if (iproj_type == PROJ_LC) then
1660 if (truelat1 /= truelat2) then
1661 colat1 = rad_per_deg*(90.0 - truelat1)
1662 colat2 = rad_per_deg*(90.0 - truelat2)
1663 n = (log(sin(colat1)) - log(sin(colat2))) &
1664 / (log(tan(colat1/2.0)) - log(tan(colat2/2.0)))
1666 do i=start_mem_i, end_mem_i
1667 do j=start_mem_j, end_mem_j
1668 colat = rad_per_deg*(90.0 - xlat_arr(i,j))
1669 mapfac_arr_x(i,j) = sin(colat2)/sin(colat)*(tan(colat/2.0)/tan(colat2/2.0))**n
1670 mapfac_arr_y(i,j) = mapfac_arr_x(i,j)
1671 end do
1672 end do
1674 else
1675 colat0 = rad_per_deg*(90.0 - truelat1)
1677 do i=start_mem_i, end_mem_i
1678 do j=start_mem_j, end_mem_j
1679 colat = rad_per_deg*(90.0 - xlat_arr(i,j))
1680 mapfac_arr_x(i,j) = sin(colat0)/sin(colat)*(tan(colat/2.0)/tan(colat0/2.0))**cos(colat0)
1681 mapfac_arr_y(i,j) = mapfac_arr_x(i,j)
1682 end do
1683 end do
1685 end if
1687 ! Polar stereographic projection
1688 else if (iproj_type == PROJ_PS) then
1690 do i=start_mem_i, end_mem_i
1691 do j=start_mem_j, end_mem_j
1692 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)))
1693 mapfac_arr_y(i,j) = mapfac_arr_x(i,j)
1694 end do
1695 end do
1697 ! Mercator projection
1698 else if (iproj_type == PROJ_MERC) then
1699 colat0 = rad_per_deg*(90.0 - truelat1)
1701 do i=start_mem_i, end_mem_i
1702 do j=start_mem_j, end_mem_j
1703 colat = rad_per_deg*(90.0 - xlat_arr(i,j))
1704 mapfac_arr_x(i,j) = sin(colat0) / sin(colat)
1705 mapfac_arr_y(i,j) = mapfac_arr_x(i,j)
1706 end do
1707 end do
1709 ! Global cylindrical projection
1710 else if (iproj_type == PROJ_CYL) then
1712 do i=start_mem_i, end_mem_i
1713 do j=start_mem_j, end_mem_j
1714 if (abs(xlat_arr(i,j)) == 90.0) then
1715 mapfac_arr_x(i,j) = 0. ! MSF actually becomes infinite at poles, but
1716 ! the values should never be used there; by
1717 ! setting to 0, we hope to induce a "divide
1718 ! by zero" error if they are
1719 else
1720 mapfac_arr_x(i,j) = 1.0 / cos(xlat_arr(i,j)*rad_per_deg)
1721 end if
1722 mapfac_arr_y(i,j) = 1.0
1723 end do
1724 end do
1726 ! Rotated global cylindrical projection
1727 else if (iproj_type == PROJ_CASSINI) then
1729 if (abs(pole_lat) == 90.) then
1730 do i=start_mem_i, end_mem_i
1731 do j=start_mem_j, end_mem_j
1732 if (abs(xlat_arr(i,j)) >= 90.0) then
1733 mapfac_arr_x(i,j) = 0. ! MSF actually becomes infinite at poles, but
1734 ! the values should never be used there; by
1735 ! setting to 0, we hope to induce a "divide
1736 ! by zero" error if they are
1737 else
1738 mapfac_arr_x(i,j) = 1.0 / cos(xlat_arr(i,j)*rad_per_deg)
1739 end if
1740 mapfac_arr_y(i,j) = 1.0
1741 end do
1742 end do
1743 else
1744 do i=start_mem_i, end_mem_i
1745 do j=start_mem_j, end_mem_j
1746 call rotate_coords(xlat_arr(i,j),xlon_arr(i,j), &
1747 comp_lat, comp_lon, &
1748 pole_lat, pole_lon, stand_lon, &
1749 -1)
1750 if (abs(comp_lat) >= 90.0) then
1751 mapfac_arr_x(i,j) = 0. ! MSF actually becomes infinite at poles, but
1752 ! the values should never be used there; by
1753 ! setting to 0, we hope to induce a "divide
1754 ! by zero" error if they are
1755 else
1756 mapfac_arr_x(i,j) = 1.0 / cos(comp_lat*rad_per_deg)
1757 end if
1758 mapfac_arr_y(i,j) = 1.0
1759 end do
1760 end do
1761 end if
1763 else if (iproj_type == PROJ_ROTLL) then
1765 do i=start_mem_i, end_mem_i
1766 do j=start_mem_j, end_mem_j
1767 mapfac_arr_x(i,j) = 1.0
1768 mapfac_arr_y(i,j) = 1.0
1769 end do
1770 end do
1772 end if
1774 end subroutine get_map_factor
1777 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1778 ! Name: get_coriolis_parameters
1779 !
1780 ! Purpose: To calculate the Coriolis parameters f and e for every gridpoint in
1781 ! the tile of the model domain
1782 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1783 subroutine get_coriolis_parameters(xlat_arr, f, e, &
1784 start_mem_i, start_mem_j, end_mem_i, end_mem_j)
1786 use constants_module
1788 implicit none
1790 ! Arguments
1791 integer, intent(in) :: start_mem_i, start_mem_j, end_mem_i, end_mem_j
1792 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(in) :: xlat_arr
1793 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(out) :: f, e
1795 ! Local variables
1796 integer :: i, j
1798 do i=start_mem_i, end_mem_i
1799 do j=start_mem_j, end_mem_j
1801 f(i,j) = 2.0*OMEGA_E*sin(rad_per_deg*xlat_arr(i,j))
1802 e(i,j) = 2.0*OMEGA_E*cos(rad_per_deg*xlat_arr(i,j))
1804 end do
1805 end do
1807 end subroutine get_coriolis_parameters
1810 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1811 ! Name: get_rotang
1812 !
1813 ! Purpose: To calculate the sine and cosine of rotation angle.
1814 !
1815 ! NOTES: The formulas used in this routine come from those in the
1816 ! vecrot_rotlat() routine of the original WRF SI.
1817 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1818 subroutine get_rotang(xlat_arr, xlon_arr, cosa, sina, &
1819 start_mem_i, start_mem_j, end_mem_i, end_mem_j)
1821 use constants_module
1822 use gridinfo_module
1824 implicit none
1826 ! Arguments
1827 integer, intent(in) :: start_mem_i, start_mem_j, end_mem_i, end_mem_j
1828 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(in) :: xlat_arr, xlon_arr
1829 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(out) :: cosa, sina
1831 ! Local variables
1832 integer :: i, j
1833 real :: alpha, d_lon
1835 do i=start_mem_i, end_mem_i
1836 do j=start_mem_j+1, end_mem_j-1
1837 d_lon = xlon_arr(i,j+1)-xlon_arr(i,j-1)
1838 if (d_lon > 180.) then
1839 d_lon = d_lon - 360.
1840 else if (d_lon < -180.) then
1841 d_lon = d_lon + 360.
1842 end if
1844 alpha = atan2(-cos(xlat_arr(i,j)*RAD_PER_DEG) * (d_lon*RAD_PER_DEG), &
1845 ((xlat_arr(i,j+1)-xlat_arr(i,j-1))*RAD_PER_DEG))
1846 sina(i,j) = sin(alpha)
1847 cosa(i,j) = cos(alpha)
1848 end do
1849 end do
1851 do i=start_mem_i, end_mem_i
1852 d_lon = xlon_arr(i,start_mem_j+1)-xlon_arr(i,start_mem_j)
1853 if (d_lon > 180.) then
1854 d_lon = d_lon - 360.
1855 else if (d_lon < -180.) then
1856 d_lon = d_lon + 360.
1857 end if
1859 alpha = atan2(-cos(xlat_arr(i,start_mem_j)*RAD_PER_DEG) * (d_lon*RAD_PER_DEG), &
1860 ((xlat_arr(i,start_mem_j+1)-xlat_arr(i,start_mem_j))*RAD_PER_DEG))
1861 sina(i,start_mem_j) = sin(alpha)
1862 cosa(i,start_mem_j) = cos(alpha)
1863 end do
1865 do i=start_mem_i, end_mem_i
1866 d_lon = xlon_arr(i,end_mem_j)-xlon_arr(i,end_mem_j-1)
1867 if (d_lon > 180.) then
1868 d_lon = d_lon - 360.
1869 else if (d_lon < -180.) then
1870 d_lon = d_lon + 360.
1871 end if
1873 alpha = atan2(-cos(xlat_arr(i,end_mem_j)*RAD_PER_DEG) * (d_lon*RAD_PER_DEG), &
1874 ((xlat_arr(i,end_mem_j)-xlat_arr(i,end_mem_j-1))*RAD_PER_DEG))
1875 sina(i,end_mem_j) = sin(alpha)
1876 cosa(i,end_mem_j) = cos(alpha)
1877 end do
1879 end subroutine get_rotang
1882 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1883 ! Name: process_neighbor
1884 !
1885 ! Purpose: This routine, give the x/y location of a point, determines whether
1886 ! the point has already been processed, and if not, which processing queue
1887 ! the point should be placed in.
1888 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1889 subroutine process_neighbor(ix, iy, bit_domain, point_queue, tile_queue, &
1890 xlat_array, xlon_array, &
1891 start_i, end_i, start_j, end_j, ilevel)
1893 use bitarray_module
1894 use misc_definitions_module
1895 use proc_point_module
1896 use queue_module
1898 implicit none
1900 ! Arguments
1901 integer, intent(in) :: ix, iy, start_i, end_i, start_j, end_j, ilevel
1902 real, dimension(start_i:end_i, start_j:end_j), intent(in) :: xlat_array, xlon_array
1903 type (bitarray), intent(inout) :: bit_domain
1904 type (queue), intent(inout) :: point_queue, tile_queue
1906 ! Local variables
1907 type (q_data) :: process_pt
1908 logical :: is_in_tile
1910 ! If the point has already been visited, no need to do anything more.
1911 if (.not. bitarray_test(bit_domain, ix-start_i+1, iy-start_j+1)) then
1913 ! Create a queue item for the current point
1914 process_pt%lat = xlat_array(ix,iy)
1915 process_pt%lon = xlon_array(ix,iy)
1916 process_pt%x = ix
1917 process_pt%y = iy
1919 is_in_tile = is_point_in_tile(process_pt%lat, process_pt%lon, ilevel)
1921 ! If the point is in the current tile, add it to the list of points
1922 ! to be processed in the inner loop
1923 if (is_in_tile) then
1924 call q_insert(point_queue, process_pt)
1925 call bitarray_set(bit_domain, ix-start_i+1, iy-start_j+1)
1927 ! Otherwise, we will process this point later. Add it to the list for
1928 ! the outer loop.
1929 else
1930 call q_insert(tile_queue, process_pt)
1931 end if
1933 end if
1935 end subroutine process_neighbor
1938 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1939 ! Name: calc_dfdy
1940 !
1941 ! Purpose: This routine calculates df/dy for the field in src_arr, and places
1942 ! the result in dst_array.
1943 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1944 subroutine calc_dfdy(src_arr, dst_arr, start_mem_i, start_mem_j, start_mem_k, &
1945 end_mem_i, end_mem_j, end_mem_k, mapfac)
1947 ! Modules
1948 use gridinfo_module
1950 implicit none
1952 ! Arguments
1953 integer, intent(in) :: start_mem_i, start_mem_j, start_mem_k, end_mem_i, end_mem_j, end_mem_k
1954 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j,start_mem_k:end_mem_k), intent(in) :: src_arr
1955 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j,start_mem_k:end_mem_k), intent(out) :: dst_arr
1956 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(in), optional :: mapfac
1958 ! Local variables
1959 integer :: i, j, k
1961 if (present(mapfac)) then
1962 do k=start_mem_k,end_mem_k
1963 do i=start_mem_i, end_mem_i
1964 do j=start_mem_j+1, end_mem_j-1
1965 dst_arr(i,j,k) = (src_arr(i,j+1,k) - src_arr(i,j-1,k))/(2.*dykm*mapfac(i,j))
1966 end do
1967 end do
1969 do i=start_mem_i, end_mem_i
1970 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))
1971 end do
1973 do i=start_mem_i, end_mem_i
1974 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))
1975 end do
1976 end do
1977 else
1978 do k=start_mem_k,end_mem_k
1979 do i=start_mem_i, end_mem_i
1980 do j=start_mem_j+1, end_mem_j-1
1981 dst_arr(i,j,k) = (src_arr(i,j+1,k) - src_arr(i,j-1,k))/(2.*dykm)
1982 end do
1983 end do
1985 do i=start_mem_i, end_mem_i
1986 dst_arr(i,start_mem_j,k) = (src_arr(i,start_mem_j+1,k) - src_arr(i,start_mem_j,k))/(dykm)
1987 end do
1989 do i=start_mem_i, end_mem_i
1990 dst_arr(i,end_mem_j,k) = (src_arr(i,end_mem_j,k) - src_arr(i,end_mem_j-1,k))/(dykm)
1991 end do
1992 end do
1993 end if
1995 end subroutine calc_dfdy
1998 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1999 ! Name: calc_dfdx
2000 !
2001 ! Purpose: This routine calculates df/dx for the field in src_arr, and places
2002 ! the result in dst_array.
2003 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2004 subroutine calc_dfdx(src_arr, dst_arr, start_mem_i, start_mem_j, &
2005 start_mem_k, end_mem_i, end_mem_j, end_mem_k, mapfac)
2007 ! Modules
2008 use gridinfo_module
2010 implicit none
2012 ! Arguments
2013 integer, intent(in) :: start_mem_i, start_mem_j, start_mem_k, end_mem_i, end_mem_j, end_mem_k
2014 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j, start_mem_k:end_mem_k), intent(in) :: src_arr
2015 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j, start_mem_k:end_mem_k), intent(out) :: dst_arr
2016 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(in), optional :: mapfac
2018 ! Local variables
2019 integer :: i, j, k
2021 if (present(mapfac)) then
2022 do k=start_mem_k, end_mem_k
2023 do i=start_mem_i+1, end_mem_i-1
2024 do j=start_mem_j, end_mem_j
2025 dst_arr(i,j,k) = (src_arr(i+1,j,k) - src_arr(i-1,j,k))/(2.*dxkm*mapfac(i,j))
2026 end do
2027 end do
2029 do j=start_mem_j, end_mem_j
2030 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))
2031 end do
2033 do j=start_mem_j, end_mem_j
2034 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))
2035 end do
2036 end do
2037 else
2038 do k=start_mem_k, end_mem_k
2039 do i=start_mem_i+1, end_mem_i-1
2040 do j=start_mem_j, end_mem_j
2041 dst_arr(i,j,k) = (src_arr(i+1,j,k) - src_arr(i-1,j,k))/(2.*dxkm)
2042 end do
2043 end do
2045 do j=start_mem_j, end_mem_j
2046 dst_arr(start_mem_i,j,k) = (src_arr(start_mem_i+1,j,k) - src_arr(start_mem_i,j,k))/(dxkm)
2047 end do
2049 do j=start_mem_j, end_mem_j
2050 dst_arr(end_mem_i,j,k) = (src_arr(end_mem_i,j,k) - src_arr(end_mem_i-1,j,k))/(dxkm)
2051 end do
2052 end do
2053 end if
2055 end subroutine calc_dfdx
2057 end module process_tile_module