| blob | 1fe930380407e5a35145eaffe84078f680023bf4 |
1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 ! Module: process_tile
3 !
4 ! Description:
5 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
6 module process_tile_module
8 use module_debug
11 contains
14 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
15 ! Name: process_tile
16 !
17 ! Purpose: To process a tile, whose lower-left corner is at
18 ! (tile_i_min, tile_j_min) and whose upper-right corner is at
19 ! (tile_i_max, tile_j_max), of the model grid given by which_domain
20 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
21 subroutine process_tile(which_domain, grid_type, dynopt, &
22 dummy_start_dom_i, dummy_end_dom_i, &
23 dummy_start_dom_j, dummy_end_dom_j, &
24 dummy_start_patch_i, dummy_end_patch_i, &
25 dummy_start_patch_j, dummy_end_patch_j, &
26 extra_col, extra_row)
28 use bitarray_module
29 use hash_module
30 use llxy_module
31 use misc_definitions_module
32 use output_module
33 use smooth_module
34 use source_data_module
36 implicit none
38 ! Arguments
39 integer, intent(in) :: which_domain, dynopt, &
40 dummy_start_dom_i, dummy_end_dom_i, dummy_start_dom_j, dummy_end_dom_j, &
41 dummy_start_patch_i, dummy_end_patch_i, dummy_start_patch_j, dummy_end_patch_j
42 logical, intent(in) :: extra_col, extra_row
43 character (len=1), intent(in) :: grid_type
45 ! Local variables
46 integer :: i, j, k, kk, istatus, ifieldstatus, idomcatstatus, field_count
47 integer :: min_category, max_category, min_level, max_level, &
48 smth_opt, smth_passes, num_landmask_categories
49 integer :: start_dom_i, end_dom_i, start_dom_j, end_dom_j, end_dom_stag_i, end_dom_stag_j
50 integer :: start_patch_i, end_patch_i, start_patch_j, end_patch_j, end_patch_stag_i, end_patch_stag_j
51 integer :: start_mem_i, end_mem_i, start_mem_j, end_mem_j, end_mem_stag_i, end_mem_stag_j
52 integer :: sm1, em1, sm2, em2
53 integer :: istagger
54 integer, dimension(MAX_LANDMASK_CATEGORIES) :: landmask_value
55 real :: sum, dominant, msg_fill_val, topo_flag_val, mass_flag, land_total, water_total
56 real, dimension(16) :: corner_lats, corner_lons
57 real, pointer, dimension(:,:) :: xlat_array, xlon_array, &
58 xlat_array_u, xlon_array_u, &
59 xlat_array_v, xlon_array_v, &
60 xlat_array_corner, xlon_array_corner, &
61 clat_array, clon_array, &
62 xlat_array_subgrid, xlon_array_subgrid, &
63 f_array, e_array, &
64 mapfac_array_m_x, mapfac_array_u_x, mapfac_array_v_x, &
65 mapfac_array_m_y, mapfac_array_u_y, mapfac_array_v_y, &
66 mapfac_array_x_subgrid, mapfac_array_y_subgrid, &
67 sina_array, cosa_array
68 real, pointer, dimension(:,:) :: xlat_ptr, xlon_ptr, mapfac_ptr_x, mapfac_ptr_y, landmask, dominant_field
69 real, pointer, dimension(:,:,:) :: field, slp_field
70 logical :: is_water_mask, only_save_dominant, halt_on_missing, &
71 is_subgrid_var, sub_status
72 character (len=19) :: datestr
73 character (len=128) :: fieldname, gradname, domname, landmask_name
74 character (len=256) :: temp_string
75 type (bitarray) :: processed_domain
76 type (hashtable) :: processed_fieldnames
77 character (len=128), dimension(2) :: dimnames
78 integer :: sub_x, sub_y
79 integer :: opt_status
81 ! Probably not all of these nullify statements are needed...
82 nullify(xlat_array)
83 nullify(xlon_array)
84 nullify(xlat_array_u)
85 nullify(xlon_array_u)
86 nullify(xlat_array_v)
87 nullify(xlon_array_v)
88 nullify(xlat_array_corner)
89 nullify(xlon_array_corner)
90 nullify(clat_array)
91 nullify(clon_array)
92 nullify(xlat_array_subgrid)
93 nullify(xlon_array_subgrid)
94 nullify(f_array)
95 nullify(e_array)
96 nullify(mapfac_array_m_x)
97 nullify(mapfac_array_u_x)
98 nullify(mapfac_array_v_x)
99 nullify(mapfac_array_m_y)
100 nullify(mapfac_array_u_y)
101 nullify(mapfac_array_v_y)
102 nullify(mapfac_array_x_subgrid)
103 nullify(mapfac_array_y_subgrid)
104 nullify(sina_array)
105 nullify(cosa_array)
106 nullify(xlat_ptr)
107 nullify(xlon_ptr)
108 nullify(mapfac_ptr_x)
109 nullify(mapfac_ptr_y)
110 nullify(landmask)
111 nullify(dominant_field)
112 nullify(field)
113 nullify(slp_field)
115 datestr = '0000-00-00_00:00:00'
116 field_count = 0
117 mass_flag=1.0
119 ! The following pertains primarily to the C grid
120 ! Determine whether only (n-1)th rows/columns should be computed for variables
121 ! on staggered grid. In a distributed memory situation, not every tile should
122 ! have only (n-1)th rows/columns computed, or we end up with (n-k)
123 ! rows/columns when there are k patches in the y/x direction
124 if (extra_col) then
125 start_patch_i = dummy_start_patch_i ! The seemingly pointless renaming of start
126 end_patch_i = dummy_end_patch_i - 1 ! naming convention with modified end_patch variables,
127 end_patch_stag_i = dummy_end_patch_i ! variables is so that we can maintain consistent
128 ! which are marked as intent(in)
129 start_mem_i = start_patch_i - HALO_WIDTH
130 end_mem_i = end_patch_i + HALO_WIDTH
131 end_mem_stag_i = end_patch_stag_i + HALO_WIDTH
132 else
133 start_patch_i = dummy_start_patch_i
134 end_patch_i = dummy_end_patch_i
135 end_patch_stag_i = dummy_end_patch_i
137 start_mem_i = start_patch_i - HALO_WIDTH
138 end_mem_i = end_patch_i + HALO_WIDTH
139 end_mem_stag_i = end_patch_stag_i + HALO_WIDTH
140 end if
142 if (extra_row) then
143 start_patch_j = dummy_start_patch_j
144 end_patch_j = dummy_end_patch_j - 1
145 end_patch_stag_j = dummy_end_patch_j
147 start_mem_j = start_patch_j - HALO_WIDTH
148 end_mem_j = end_patch_j + HALO_WIDTH
149 end_mem_stag_j = end_patch_stag_j + HALO_WIDTH
150 else
151 start_patch_j = dummy_start_patch_j
152 end_patch_j = dummy_end_patch_j
153 end_patch_stag_j = dummy_end_patch_j
155 start_mem_j = start_patch_j - HALO_WIDTH
156 end_mem_j = end_patch_j + HALO_WIDTH
157 end_mem_stag_j = end_patch_stag_j + HALO_WIDTH
158 end if
160 start_dom_i = dummy_start_dom_i
161 if (grid_type == 'C') then
162 end_dom_i = dummy_end_dom_i - 1
163 end_dom_stag_i = dummy_end_dom_i
164 else if (grid_type == 'E') then
165 end_dom_i = dummy_end_dom_i
166 end_dom_stag_i = dummy_end_dom_i
167 end if
169 start_dom_j = dummy_start_dom_j
170 if (grid_type == 'C') then
171 end_dom_j = dummy_end_dom_j - 1
172 end_dom_stag_j = dummy_end_dom_j
173 else if (grid_type == 'E') then
174 end_dom_j = dummy_end_dom_j
175 end_dom_stag_j = dummy_end_dom_j
176 end if
178 ! Allocate arrays to hold all lat/lon fields; these will persist for the duration of
179 ! the process_tile routine
180 ! For C grid, we have M, U, and V points
181 ! For E grid, we have only M and V points
182 allocate(xlat_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
183 allocate(xlon_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
184 allocate(xlat_array_v(start_mem_i:end_mem_i, start_mem_j:end_mem_stag_j))
185 allocate(xlon_array_v(start_mem_i:end_mem_i, start_mem_j:end_mem_stag_j))
186 if (grid_type == 'C') then
187 allocate(xlat_array_u(start_mem_i:end_mem_stag_i, start_mem_j:end_mem_j))
188 allocate(xlon_array_u(start_mem_i:end_mem_stag_i, start_mem_j:end_mem_j))
189 allocate(clat_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
190 allocate(clon_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
191 allocate(xlat_array_corner(start_mem_i:end_mem_stag_i, start_mem_j:end_mem_stag_j))
192 allocate(xlon_array_corner(start_mem_i:end_mem_stag_i, start_mem_j:end_mem_stag_j))
193 end if
194 nullify(xlat_array_subgrid)
195 nullify(xlon_array_subgrid)
196 nullify(mapfac_array_x_subgrid)
197 nullify(mapfac_array_y_subgrid)
199 ! Initialize hash table to track which fields have been processed
200 call hash_init(processed_fieldnames)
202 !
203 ! Calculate lat/lon for every point in the tile (XLAT and XLON)
204 ! The xlat_array and xlon_array arrays will be used in processing other fields
205 !
206 call mprintf(.true.,STDOUT,' Processing XLAT and XLONG')
208 if (grid_type == 'C') then
209 call get_lat_lon_fields(xlat_array, xlon_array, start_mem_i, &
210 start_mem_j, end_mem_i, end_mem_j, M)
211 call get_lat_lon_fields(xlat_array_v, xlon_array_v, start_mem_i, &
212 start_mem_j, end_mem_i, end_mem_stag_j, V)
213 call get_lat_lon_fields(xlat_array_u, xlon_array_u, start_mem_i, &
214 start_mem_j, end_mem_stag_i, end_mem_j, U)
215 call get_lat_lon_fields(xlat_array_corner, xlon_array_corner, start_mem_i, &
216 start_mem_j, end_mem_stag_i, end_mem_stag_j, CORNER)
217 call get_lat_lon_fields(clat_array, clon_array, start_mem_i, &
218 start_mem_j, end_mem_i, end_mem_j, M, comp_ll=.true.)
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_u(start_patch_i,start_patch_j)
226 corner_lats(6) = xlat_array_u(start_patch_i,end_patch_j)
227 corner_lats(7) = xlat_array_u(end_patch_stag_i,end_patch_j)
228 corner_lats(8) = xlat_array_u(end_patch_stag_i,start_patch_j)
230 corner_lats(9) = xlat_array_v(start_patch_i,start_patch_j)
231 corner_lats(10) = xlat_array_v(start_patch_i,end_patch_stag_j)
232 corner_lats(11) = xlat_array_v(end_patch_i,end_patch_stag_j)
233 corner_lats(12) = xlat_array_v(end_patch_i,start_patch_j)
235 call xytoll(real(start_patch_i)-0.5, real(start_patch_j)-0.5, corner_lats(13), corner_lons(13), M)
236 call xytoll(real(start_patch_i)-0.5, real(end_patch_j)+0.5, corner_lats(14), corner_lons(14), M)
237 call xytoll(real(end_patch_i)+0.5, real(end_patch_j)+0.5, corner_lats(15), corner_lons(15), M)
238 call xytoll(real(end_patch_i)+0.5, real(start_patch_j)-0.5, corner_lats(16), corner_lons(16), M)
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_u(start_patch_i,start_patch_j)
246 corner_lons(6) = xlon_array_u(start_patch_i,end_patch_j)
247 corner_lons(7) = xlon_array_u(end_patch_stag_i,end_patch_j)
248 corner_lons(8) = xlon_array_u(end_patch_stag_i,start_patch_j)
250 corner_lons(9) = xlon_array_v(start_patch_i,start_patch_j)
251 corner_lons(10) = xlon_array_v(start_patch_i,end_patch_stag_j)
252 corner_lons(11) = xlon_array_v(end_patch_i,end_patch_stag_j)
253 corner_lons(12) = xlon_array_v(end_patch_i,start_patch_j)
255 else if (grid_type == 'E') then
256 call get_lat_lon_fields(xlat_array, xlon_array, start_mem_i, &
257 start_mem_j, end_mem_i, end_mem_j, HH)
258 call get_lat_lon_fields(xlat_array_v, xlon_array_v, start_mem_i, &
259 start_mem_j, end_mem_i, end_mem_stag_j, VV)
261 corner_lats(1) = xlat_array(start_patch_i,start_patch_j)
262 corner_lats(2) = xlat_array(start_patch_i,end_patch_j)
263 corner_lats(3) = xlat_array(end_patch_i,end_patch_j)
264 corner_lats(4) = xlat_array(end_patch_i,start_patch_j)
266 corner_lats(5) = xlat_array_v(start_patch_i,start_patch_j)
267 corner_lats(6) = xlat_array_v(start_patch_i,end_patch_stag_j)
268 corner_lats(7) = xlat_array_v(end_patch_i,end_patch_stag_j)
269 corner_lats(8) = xlat_array_v(end_patch_i,start_patch_j)
271 corner_lats(9) = 0.0
272 corner_lats(10) = 0.0
273 corner_lats(11) = 0.0
274 corner_lats(12) = 0.0
276 corner_lats(13) = 0.0
277 corner_lats(14) = 0.0
278 corner_lats(15) = 0.0
279 corner_lats(16) = 0.0
281 corner_lons(1) = xlon_array(start_patch_i,start_patch_j)
282 corner_lons(2) = xlon_array(start_patch_i,end_patch_j)
283 corner_lons(3) = xlon_array(end_patch_i,end_patch_j)
284 corner_lons(4) = xlon_array(end_patch_i,start_patch_j)
286 corner_lons(5) = xlon_array_v(start_patch_i,start_patch_j)
287 corner_lons(6) = xlon_array_v(start_patch_i,end_patch_stag_j)
288 corner_lons(7) = xlon_array_v(end_patch_i,end_patch_stag_j)
289 corner_lons(8) = xlon_array_v(end_patch_i,start_patch_j)
291 corner_lons(9) = 0.0
292 corner_lons(10) = 0.0
293 corner_lons(11) = 0.0
294 corner_lons(12) = 0.0
296 corner_lons(13) = 0.0
297 corner_lons(14) = 0.0
298 corner_lons(15) = 0.0
299 corner_lons(16) = 0.0
301 end if
303 ! Initialize the output module now that we have the corner point lats/lons
304 call output_init(which_domain, 'OUTPUT FROM GEOGRID V4.5', '0000-00-00_00:00:00', grid_type, dynopt, &
305 corner_lats, corner_lons, &
306 start_dom_i, end_dom_i, start_dom_j, end_dom_j, &
307 start_patch_i, end_patch_i, start_patch_j, end_patch_j, &
308 start_mem_i, end_mem_i, start_mem_j, end_mem_j, &
309 extra_col, extra_row)
311 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, &
312 'XLAT_M', datestr, real_array = xlat_array)
313 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, &
314 'XLONG_M', datestr, real_array = xlon_array)
315 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_stag_j, 1, 1, &
316 'XLAT_V', datestr, real_array = xlat_array_v)
317 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_stag_j, 1, 1, &
318 'XLONG_V', datestr, real_array = xlon_array_v)
319 if (grid_type == 'C') then
320 call write_field(start_mem_i, end_mem_stag_i, start_mem_j, end_mem_j, 1, 1, &
321 'XLAT_U', datestr, real_array = xlat_array_u)
322 call write_field(start_mem_i, end_mem_stag_i, start_mem_j, end_mem_j, 1, 1, &
323 'XLONG_U', datestr, real_array = xlon_array_u)
324 call write_field(start_mem_i, end_mem_stag_i, start_mem_j, end_mem_stag_j, 1, 1, &
325 'XLAT_C', datestr, real_array = xlat_array_corner)
326 call write_field(start_mem_i, end_mem_stag_i, start_mem_j, end_mem_stag_j, 1, 1, &
327 'XLONG_C', datestr, real_array = xlon_array_corner)
328 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, &
329 'CLAT', datestr, real_array = clat_array)
330 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, &
331 'CLONG', datestr, real_array = clon_array)
333 if (associated(clat_array)) deallocate(clat_array)
334 if (associated(clon_array)) deallocate(clon_array)
336 end if
339 !
340 ! Calculate map factor for current domain
341 !
342 if (grid_type == 'C') then
343 call mprintf(.true.,STDOUT,' Processing MAPFAC')
345 allocate(mapfac_array_m_x(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
346 allocate(mapfac_array_m_y(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
347 call get_map_factor(xlat_array, xlon_array, mapfac_array_m_x, mapfac_array_m_y, start_mem_i, &
348 start_mem_j, end_mem_i, end_mem_j)
349 ! Global WRF uses map scale factors in X and Y directions, but "regular" WRF uses a single MSF
350 ! on each staggering. In the case of regular WRF, we can assume that MAPFAC_MX = MAPFAC_MY = MAPFAC_M,
351 ! and so we can simply write MAPFAC_MX as the MAPFAC_M field. Ultimately, when global WRF is
352 ! merged into the WRF trunk, we will need only two map scale factor fields for each staggering,
353 ! in the x and y directions, and these will be the same in the case of non-Cassini projections
354 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'MAPFAC_M', &
355 datestr, real_array = mapfac_array_m_x)
356 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'MAPFAC_MX', &
357 datestr, real_array = mapfac_array_m_x)
358 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'MAPFAC_MY', &
359 datestr, real_array = mapfac_array_m_y)
361 allocate(mapfac_array_v_x(start_mem_i:end_mem_i, start_mem_j:end_mem_stag_j))
362 allocate(mapfac_array_v_y(start_mem_i:end_mem_i, start_mem_j:end_mem_stag_j))
363 call get_map_factor(xlat_array_v, xlon_array_v, mapfac_array_v_x, mapfac_array_v_y, start_mem_i, &
364 start_mem_j, end_mem_i, end_mem_stag_j)
365 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_stag_j, 1, 1, 'MAPFAC_V', &
366 datestr, real_array = mapfac_array_v_x)
367 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_stag_j, 1, 1, 'MAPFAC_VX', &
368 datestr, real_array = mapfac_array_v_x)
369 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_stag_j, 1, 1, 'MAPFAC_VY', &
370 datestr, real_array = mapfac_array_v_y)
372 allocate(mapfac_array_u_x(start_mem_i:end_mem_stag_i, start_mem_j:end_mem_j))
373 allocate(mapfac_array_u_y(start_mem_i:end_mem_stag_i, start_mem_j:end_mem_j))
374 call get_map_factor(xlat_array_u, xlon_array_u, mapfac_array_u_x, mapfac_array_u_y, start_mem_i, &
375 start_mem_j, end_mem_stag_i, end_mem_j)
376 call write_field(start_mem_i, end_mem_stag_i, start_mem_j, end_mem_j, 1, 1, 'MAPFAC_U', &
377 datestr, real_array = mapfac_array_u_x)
378 call write_field(start_mem_i, end_mem_stag_i, start_mem_j, end_mem_j, 1, 1, 'MAPFAC_UX', &
379 datestr, real_array = mapfac_array_u_x)
380 call write_field(start_mem_i, end_mem_stag_i, start_mem_j, end_mem_j, 1, 1, 'MAPFAC_UY', &
381 datestr, real_array = mapfac_array_u_y)
383 end if
386 !
387 ! Coriolis parameters (E and F)
388 !
389 call mprintf(.true.,STDOUT,' Processing F and E')
391 allocate(f_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
392 allocate(e_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
394 call get_coriolis_parameters(xlat_array, f_array, e_array, &
395 start_mem_i, start_mem_j, end_mem_i, end_mem_j)
397 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'E', &
398 datestr, real_array = e_array)
399 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'F', &
400 datestr, real_array = f_array)
402 if (associated(f_array)) deallocate(f_array)
403 if (associated(e_array)) deallocate(e_array)
406 !
407 ! Rotation angle (SINALPHA and COSALPHA)
408 !
409 if (grid_type == 'C') then
410 call mprintf(.true.,STDOUT,' Processing ROTANG')
412 ! Mass-staggered points
413 allocate(sina_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
414 allocate(cosa_array(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
416 call get_rotang(xlat_array, xlon_array, cosa_array, sina_array, &
417 start_mem_i, start_mem_j, end_mem_i, end_mem_j)
419 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'SINALPHA', &
420 datestr, real_array = sina_array)
421 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'COSALPHA', &
422 datestr, real_array = cosa_array)
424 if (associated(sina_array)) deallocate(sina_array)
425 if (associated(cosa_array)) deallocate(cosa_array)
427 ! U-staggered points
428 allocate(sina_array(start_mem_i:end_mem_stag_i, start_mem_j:end_mem_j))
429 allocate(cosa_array(start_mem_i:end_mem_stag_i, start_mem_j:end_mem_j))
431 call get_rotang(xlat_array_u, xlon_array_u, cosa_array, sina_array, &
432 start_mem_i, start_mem_j, end_mem_stag_i, end_mem_j)
434 call write_field(start_mem_i, end_mem_stag_i, start_mem_j, end_mem_j, 1, 1, 'SINALPHA_U', &
435 datestr, real_array = sina_array)
436 call write_field(start_mem_i, end_mem_stag_i, start_mem_j, end_mem_j, 1, 1, 'COSALPHA_U', &
437 datestr, real_array = cosa_array)
439 if (associated(sina_array)) deallocate(sina_array)
440 if (associated(cosa_array)) deallocate(cosa_array)
442 ! V-staggered points
443 allocate(sina_array(start_mem_i:end_mem_i, start_mem_j:end_mem_stag_j))
444 allocate(cosa_array(start_mem_i:end_mem_i, start_mem_j:end_mem_stag_j))
446 call get_rotang(xlat_array_v, xlon_array_v, cosa_array, sina_array, &
447 start_mem_i, start_mem_j, end_mem_i, end_mem_stag_j)
449 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_stag_j, 1, 1, 'SINALPHA_V', &
450 datestr, real_array = sina_array)
451 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_stag_j, 1, 1, 'COSALPHA_V', &
452 datestr, real_array = cosa_array)
454 if (associated(sina_array)) deallocate(sina_array)
455 if (associated(cosa_array)) deallocate(cosa_array)
456 end if
458 ! Every field up until now should probably just be processed regardless of what the user
459 ! has specified for fields to be processed.
460 ! Hereafter, we process user-specified fields
462 !
463 ! First process the field that we will derive a landmask from
464 !
465 call get_landmask_field(geog_data_res(which_domain), landmask_name, is_water_mask, landmask_value, istatus)
467 do kk=1,MAX_LANDMASK_CATEGORIES
468 if (landmask_value(kk) == INVALID) then
469 num_landmask_categories = kk-1
470 exit
471 end if
472 end do
473 if (kk > MAX_LANDMASK_CATEGORIES) num_landmask_categories = MAX_LANDMASK_CATEGORIES
475 if (istatus /= 0) then
476 call mprintf(.true.,WARN,'No field specified for landmask calculation. Will set landmask=1 at every grid point.')
478 allocate(landmask(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
479 landmask = 1.
480 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'LANDMASK', &
481 datestr, landmask)
483 else
485 allocate(landmask(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
486 landmask = 1.
488 call mprintf(.true.,STDOUT,' Processing %s', s1=trim(landmask_name))
490 call get_missing_fill_value(landmask_name, msg_fill_val, istatus)
491 if (istatus /= 0) msg_fill_val = NAN
493 call get_halt_on_missing(landmask_name, halt_on_missing, istatus)
494 if (istatus /= 0) halt_on_missing = .false.
496 ! Do we calculate a dominant category for this field?
497 call get_domcategory_name(landmask_name, domname, only_save_dominant, idomcatstatus)
499 temp_string = ' '
500 temp_string(1:128) = landmask_name
501 call hash_insert(processed_fieldnames, temp_string)
503 call get_max_categories(landmask_name, min_category, max_category, istatus)
504 allocate(field(start_mem_i:end_mem_i, start_mem_j:end_mem_j, min_category:max_category))
506 if (.not. only_save_dominant) then
507 field_count = field_count + 1
508 call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
509 i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=landmask_name)
510 else
511 field_count = field_count + 1
512 call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
513 i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=domname)
514 end if
516 if (grid_type == 'C') then
517 call calc_field(landmask_name, field, xlat_array, xlon_array, M, &
518 start_mem_i, end_mem_i, start_mem_j, end_mem_j, &
519 min_category, max_category, processed_domain, 1, landmask=landmask, sr_x=1, sr_y=1)
520 else if (grid_type == 'E') then
521 call calc_field(landmask_name, field, xlat_array, xlon_array, HH, &
522 start_mem_i, end_mem_i, start_mem_j, end_mem_j, &
523 min_category, max_category, processed_domain, 1, landmask=landmask, sr_x=1, sr_y=1)
524 end if
526 ! If user wants to halt when a missing value is found in output field, check now
527 if (halt_on_missing) then
528 do i=start_mem_i, end_mem_i
529 do j=start_mem_j, end_mem_j
530 ! Only need to examine k=1
531 if (field(i,j,1) == msg_fill_val) then
532 call mprintf(.true.,ERROR,' Missing value encountered in output field. Quitting.')
533 end if
534 end do
535 end do
536 end if
538 ! Find fractions
539 do i=start_mem_i, end_mem_i
540 do j=start_mem_j, end_mem_j
541 sum = 0.
542 do k=min_category,max_category
543 sum = sum + field(i,j,k)
544 end do
545 if (sum > 0.0) then
546 do k=min_category,max_category
547 field(i,j,k) = field(i,j,k) / sum
548 end do
549 else
550 do k=min_category,max_category
551 field(i,j,k) = msg_fill_val
552 end do
553 end if
554 end do
555 end do
557 if (is_water_mask) then
558 call mprintf(.true.,STDOUT,' Calculating landmask from %s ( WATER =', &
559 newline=.false.,s1=trim(landmask_name))
560 else
561 call mprintf(.true.,STDOUT,' Calculating landmask from %s ( LAND =', &
562 newline=.false.,s1=trim(landmask_name))
563 end if
564 do k = 1, num_landmask_categories
565 call mprintf(.true.,STDOUT,' %i',newline=.false.,i1=landmask_value(k))
566 if (k == num_landmask_categories) call mprintf(.true.,STDOUT,')')
567 end do
569 ! Calculate landmask
570 if (is_water_mask) then
571 do i=start_mem_i, end_mem_i
572 do j=start_mem_j, end_mem_j
573 water_total = -1.
574 do k=1,num_landmask_categories
575 if (landmask_value(k) >= min_category .and. landmask_value(k) <= max_category) then
576 if (field(i,j,landmask_value(k)) /= msg_fill_val) then
577 if (water_total < 0.) water_total = 0.
578 water_total = water_total + field(i,j,landmask_value(k))
579 else
580 water_total = -1.
581 exit
582 end if
583 end if
584 end do
585 if (water_total >= 0.0) then
586 if (water_total < 0.50) then
587 landmask(i,j) = 1.
588 else
589 landmask(i,j) = 0.
590 end if
591 else
592 landmask(i,j) = -1.
593 end if
594 end do
595 end do
596 else
597 do i=start_mem_i, end_mem_i
598 do j=start_mem_j, end_mem_j
599 land_total = -1.
600 do k=1,num_landmask_categories
601 if (landmask_value(k) >= min_category .and. landmask_value(k) <= max_category) then
602 if (field(i,j,landmask_value(k)) /= msg_fill_val) then
603 if (land_total < 0.) land_total = 0.
604 land_total = land_total + field(i,j,landmask_value(k))
605 else
606 land_total = -1.
607 exit
608 end if
609 end if
610 end do
611 if (land_total >= 0.0) then
612 if (land_total > 0.50) then
613 landmask(i,j) = 1.
614 else
615 landmask(i,j) = 0.
616 end if
617 else
618 landmask(i,j) = -1.
619 end if
620 end do
621 end do
622 end if
624 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, 'LANDMASK', &
625 datestr, landmask)
627 ! If we should only save the dominant category, then no need to write out fractional field
628 if (.not.only_save_dominant .or. (idomcatstatus /= 0)) then
630 ! Finally, we may be asked to smooth the fractional field
631 call get_smooth_option(landmask_name, smth_opt, smth_passes, istatus)
632 if (istatus == 0) then
634 if (grid_type == 'C') then
635 if (smth_opt == ONETWOONE) then
636 call one_two_one(field, &
637 start_patch_i, end_patch_i, &
638 start_patch_j, end_patch_j, &
639 start_mem_i, end_mem_i, &
640 start_mem_j, end_mem_j, &
641 min_category, max_category, &
642 smth_passes, msg_fill_val)
643 else if (smth_opt == SMTHDESMTH) then
644 call smth_desmth(field, &
645 start_patch_i, end_patch_i, &
646 start_patch_j, end_patch_j, &
647 start_mem_i, end_mem_i, &
648 start_mem_j, end_mem_j, &
649 min_category, max_category, &
650 smth_passes, msg_fill_val)
651 else if (smth_opt == SMTHDESMTH_SPECIAL) then
652 call smth_desmth_special(field, &
653 start_patch_i, end_patch_i, &
654 start_patch_j, end_patch_j, &
655 start_mem_i, end_mem_i, &
656 start_mem_j, end_mem_j, &
657 min_category, max_category, &
658 smth_passes, msg_fill_val)
659 end if
660 else if (grid_type == 'E') then
661 if (smth_opt == ONETWOONE) then
662 call one_two_one_egrid(field, &
663 start_patch_i, end_patch_i, &
664 start_patch_j, end_patch_j, &
665 start_mem_i, end_mem_i, &
666 start_mem_j, end_mem_j, &
667 min_category, max_category, &
668 smth_passes, msg_fill_val, 1.0)
669 else if (smth_opt == SMTHDESMTH) then
670 call smth_desmth_egrid(field, &
671 start_patch_i, end_patch_i, &
672 start_patch_j, end_patch_j, &
673 start_mem_i, end_mem_i, &
674 start_mem_j, end_mem_j, &
675 min_category, max_category, &
676 smth_passes, msg_fill_val, 1.0)
677 else if (smth_opt == SMTHDESMTH_SPECIAL) then
678 call mprintf(.true.,WARN,'smth-desmth_special is not currently implemented for NMM. '// &
679 'No smoothing will be done.')
680 end if
681 end if
683 end if
685 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, &
686 min_category, max_category, trim(landmask_name), &
687 datestr, real_array=field)
688 end if
690 if (idomcatstatus == 0) then
691 allocate(dominant_field(start_mem_i:end_mem_i, start_mem_j:end_mem_j))
693 if (.not. only_save_dominant) then
694 field_count = field_count + 1
695 call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
696 i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=domname)
697 end if
699 do i=start_mem_i, end_mem_i
700 do j=start_mem_j, end_mem_j
701 if ((landmask(i,j) == 1. .and. is_water_mask) .or. &
702 (landmask(i,j) == 0. .and. .not.is_water_mask)) then
703 dominant = 0.
704 dominant_field(i,j) = real(min_category-1)
705 do k=min_category,max_category
706 do kk=1,num_landmask_categories
707 if (k == landmask_value(kk)) exit
708 end do
709 if (field(i,j,k) > dominant .and. kk > num_landmask_categories) then
710 dominant_field(i,j) = real(k)
711 dominant = field(i,j,k)
712 end if
713 end do
714 else
715 dominant = 0.
716 dominant_field(i,j) = real(min_category-1)
717 do k=min_category,max_category
718 do kk=1,num_landmask_categories
719 if (field(i,j,k) > dominant .and. k == landmask_value(kk)) then
720 dominant_field(i,j) = real(k)
721 dominant = field(i,j,k)
722 end if
723 end do
724 end do
725 end if
726 end do
727 end do
729 call write_field(start_mem_i, end_mem_i, start_mem_j, end_mem_j, 1, 1, trim(domname), &
730 datestr, dominant_field)
732 deallocate(dominant_field)
733 end if
735 deallocate(field)
736 end if
738 !
739 ! Now process all other fields specified by the user
740 !
741 call reset_next_field()
742 ifieldstatus = 0
743 do while (ifieldstatus == 0)
744 call get_next_fieldname(fieldname, ifieldstatus)
746 ! There is another field in the GEOGRID.TBL file
747 if (ifieldstatus == 0) then
748 temp_string(1:128) = fieldname
750 call get_source_opt_status(fieldname, 0, opt_status)
751 call get_output_stagger(fieldname, istagger, istatus)
752 dimnames(:) = 'null'
753 call get_subgrid_dim_name(which_domain, fieldname, dimnames, &
754 sub_x, sub_y, is_subgrid_var, istatus)
755 sub_status = (.not. is_subgrid_var .or. (sub_x > 0 .or. sub_y > 0))
757 ! If this field is still to be processed
758 if (.not. hash_search(processed_fieldnames, temp_string) .and. opt_status == 0 .and. sub_status) then
760 call hash_insert(processed_fieldnames, temp_string)
761 call mprintf(.true.,STDOUT,' Processing %s', s1=trim(fieldname))
763 if (istagger == M .or. (sub_x > 1) .or. (sub_y > 1)) then
764 sm1 = start_mem_i
765 em1 = end_mem_i
766 sm2 = start_mem_j
767 em2 = end_mem_j
768 xlat_ptr => xlat_array
769 xlon_ptr => xlon_array
770 mapfac_ptr_x => mapfac_array_m_x
771 mapfac_ptr_y => mapfac_array_m_y
772 else if (istagger == U) then ! In the case that extra_cols = .false.
773 sm1 = start_mem_i ! we should have that end_mem_stag is
774 em1 = end_mem_stag_i ! the same as end_mem, so we do not need
775 sm2 = start_mem_j ! to check extra_cols or extra rows here
776 em2 = end_mem_j
777 xlat_ptr => xlat_array_u
778 xlon_ptr => xlon_array_u
779 mapfac_ptr_x => mapfac_array_u_x
780 mapfac_ptr_y => mapfac_array_u_y
781 else if (istagger == V) then
782 sm1 = start_mem_i
783 em1 = end_mem_i
784 sm2 = start_mem_j
785 em2 = end_mem_stag_j
786 xlat_ptr => xlat_array_v
787 xlon_ptr => xlon_array_v
788 mapfac_ptr_x => mapfac_array_v_x
789 mapfac_ptr_y => mapfac_array_v_y
790 else if (istagger == HH) then ! E grid
791 sm1 = start_mem_i
792 em1 = end_mem_i
793 sm2 = start_mem_j
794 em2 = end_mem_j
795 xlat_ptr => xlat_array
796 xlon_ptr => xlon_array
797 mapfac_ptr_x => mapfac_array_m_x
798 mapfac_ptr_y => mapfac_array_m_y
799 else if (istagger == VV) then ! E grid
800 sm1 = start_mem_i
801 em1 = end_mem_i
802 sm2 = start_mem_j
803 em2 = end_mem_stag_j
804 xlat_ptr => xlat_array_v
805 xlon_ptr => xlon_array_v
806 mapfac_ptr_x => mapfac_array_v_x
807 mapfac_ptr_y => mapfac_array_v_y
808 end if
810 if (sub_x > 1) then
811 sm1 = (start_mem_i - 1)*sub_x + 1
812 if (extra_col) then
813 em1 = (end_mem_i + 1)*sub_x
814 else
815 em1 = (end_mem_i )*sub_x
816 end if
817 end if
818 if (sub_y > 1)then
819 sm2 = (start_mem_j - 1)*sub_y + 1
820 if (extra_row) then
821 em2 = (end_mem_j + 1)*sub_y
822 else
823 em2 = (end_mem_j )*sub_y
824 end if
825 end if
827 !BUG: This should probably be moved up to where other lat/lon fields are calculated, and we should
828 ! just determine whether we will have any subgrids or not at that point
829 if ((sub_x > 1) .or. (sub_y > 1)) then
830 ! if (associated(xlat_array_subgrid)) deallocate(xlat_array_subgrid)
831 ! if (associated(xlon_array_subgrid)) deallocate(xlon_array_subgrid)
832 ! if (associated(mapfac_array_x_subgrid)) deallocate(mapfac_array_x_subgrid)
833 ! if (associated(mapfac_array_y_subgrid)) deallocate(mapfac_array_y_subgrid)
834 allocate(xlat_array_subgrid(sm1:em1,sm2:em2))
835 allocate(xlon_array_subgrid(sm1:em1,sm2:em2))
836 allocate(mapfac_array_x_subgrid(sm1:em1,sm2:em2))
837 allocate(mapfac_array_y_subgrid(sm1:em1,sm2:em2))
838 call get_lat_lon_fields(xlat_array_subgrid, xlon_array_subgrid, &
839 sm1, sm2, em1, em2, M, sub_x=sub_x, sub_y=sub_y)
840 xlat_ptr => xlat_array_subgrid
841 xlon_ptr => xlon_array_subgrid
842 call get_map_factor(xlat_ptr, xlon_ptr, mapfac_array_x_subgrid, &
843 mapfac_array_y_subgrid, sm1, sm2, em1, em2)
844 mapfac_ptr_x => mapfac_array_x_subgrid
845 mapfac_ptr_y => mapfac_array_y_subgrid
846 end if
848 call get_missing_fill_value(fieldname, msg_fill_val, istatus)
849 if (istatus /= 0) msg_fill_val = NAN
851 call get_halt_on_missing(fieldname, halt_on_missing, istatus)
852 if (istatus /= 0) halt_on_missing = .false.
854 ! Destination field type is CONTINUOUS
855 if (iget_fieldtype(fieldname,istatus) == CONTINUOUS) then
856 call get_max_levels(fieldname, min_level, max_level, istatus)
857 allocate(field(sm1:em1, sm2:em2, min_level:max_level))
859 field_count = field_count + 1
860 call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
861 i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=fieldname)
863 if ((sub_x > 1) .or. (sub_y > 1)) then
864 call calc_field(fieldname, field, xlat_ptr, xlon_ptr, istagger, &
865 sm1, em1, sm2, em2, min_level, max_level, &
866 processed_domain, 1, sr_x=sub_x, sr_y=sub_y)
867 else
868 call calc_field(fieldname, field, xlat_ptr, xlon_ptr, istagger, &
869 sm1, em1, sm2, em2, min_level, max_level, &
870 processed_domain, 1, landmask=landmask, sr_x=sub_x, sr_y=sub_y)
871 end if
873 ! If user wants to halt when a missing value is found in output field, check now
874 if (halt_on_missing) then
875 do i=sm1, em1
876 do j=sm2, em2
877 ! Only need to examine k=1
878 if (field(i,j,1) == msg_fill_val) then
879 call mprintf(.true.,ERROR,' Missing value encountered in output field. Quitting.')
880 end if
881 end do
882 end do
883 end if
885 ! We may be asked to smooth the fractional field
886 call get_smooth_option(fieldname, smth_opt, smth_passes, istatus)
887 if (istatus == 0) then
889 if (grid_type == 'C') then
890 if (smth_opt == ONETWOONE) then
891 call one_two_one(field, &
892 start_patch_i, end_patch_i, &
893 start_patch_j, end_patch_j, &
894 sm1, em1, &
895 sm2, em2, &
896 min_level, max_level, &
897 smth_passes, msg_fill_val)
898 else if (smth_opt == SMTHDESMTH) then
899 call smth_desmth(field, &
900 start_patch_i, end_patch_i, &
901 start_patch_j, end_patch_j, &
902 sm1, em1, &
903 sm2, em2, &
904 min_level, max_level, &
905 smth_passes, msg_fill_val)
906 else if (smth_opt == SMTHDESMTH_SPECIAL) then
907 call smth_desmth_special(field, &
908 start_patch_i, end_patch_i, &
909 start_patch_j, end_patch_j, &
910 sm1, em1, &
911 sm2, em2, &
912 min_level, max_level, &
913 smth_passes, msg_fill_val)
914 end if
916 else if (grid_type == 'E') then
918 if (trim(fieldname) == 'HGT_M' ) then
919 topo_flag_val=1.0
920 mass_flag=1.0
921 else if (trim(fieldname) == 'HGT_V') then
922 topo_flag_val=1.0
923 mass_flag=0.0
924 else
925 topo_flag_val=0.0
926 end if
928 if (smth_opt == ONETWOONE) then
929 call one_two_one_egrid(field, &
930 start_patch_i, end_patch_i, &
931 start_patch_j, end_patch_j, &
932 sm1, em1, &
933 sm2, em2, &
934 min_level, max_level, &
935 smth_passes, topo_flag_val, mass_flag)
936 else if (smth_opt == SMTHDESMTH) then
937 call smth_desmth_egrid(field, &
938 start_patch_i, end_patch_i, &
939 start_patch_j, end_patch_j, &
940 sm1, em1, &
941 sm2, em2, &
942 min_level, max_level, &
943 smth_passes, topo_flag_val, mass_flag)
944 else if (smth_opt == SMTHDESMTH_SPECIAL) then
945 call mprintf(.true.,WARN,'smth-desmth_special is not currently implemented for NMM. '// &
946 'No smoothing will be done.')
947 end if
949 end if
951 end if
953 call write_field(sm1, em1, sm2, em2, &
954 min_level, max_level, trim(fieldname), datestr, real_array=field)
956 ! Do we calculate directional derivatives from this field?
957 call get_dfdx_name(fieldname, gradname, istatus)
958 if (istatus == 0) then
959 allocate(slp_field(sm1:em1,sm2:em2,min_level:max_level))
961 field_count = field_count + 1
962 call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
963 i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=gradname)
965 if (grid_type == 'C') then
966 call calc_dfdx(field, slp_field, sm1, sm2, min_level, em1, em2, max_level, sub_x, mapfac_ptr_x)
967 else if (grid_type == 'E') then
968 call calc_dfdx(field, slp_field, sm1, sm2, min_level, em1, em2, max_level, sub_x)
969 end if
970 call write_field(sm1, em1, sm2, em2, &
971 min_level, max_level, trim(gradname), datestr, real_array=slp_field)
972 deallocate(slp_field)
973 end if
974 call get_dfdy_name(fieldname, gradname, istatus)
975 if (istatus == 0) then
976 allocate(slp_field(sm1:em1,sm2:em2,min_level:max_level))
978 field_count = field_count + 1
979 call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
980 i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=gradname)
982 if (grid_type == 'C') then
983 call calc_dfdy(field, slp_field, sm1, sm2, min_level, em1, em2, max_level, sub_y, mapfac_ptr_y)
984 else if (grid_type == 'E') then
985 call calc_dfdy(field, slp_field, sm1, sm2, min_level, em1, em2, max_level, sub_y)
986 end if
987 call write_field(sm1, em1, sm2, em2, &
988 min_level, max_level, trim(gradname), datestr, real_array=slp_field)
989 deallocate(slp_field)
990 end if
992 deallocate(field)
994 ! Destination field type is CATEGORICAL
995 else
996 call get_max_categories(fieldname, min_category, max_category, istatus)
997 allocate(field(sm1:em1, sm2:em2, min_category:max_category))
999 ! Do we calculate a dominant category for this field?
1000 call get_domcategory_name(fieldname, domname, only_save_dominant, idomcatstatus)
1002 if (.not. only_save_dominant) then
1003 field_count = field_count + 1
1004 call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
1005 i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=fieldname)
1006 else
1007 field_count = field_count + 1
1008 call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
1009 i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=domname)
1010 end if
1012 if ((sub_x > 1) .or. (sub_y > 1)) then
1013 call calc_field(fieldname, field, xlat_ptr, xlon_ptr, istagger, &
1014 sm1, em1, sm2, em2, min_category, max_category, &
1015 processed_domain, 1, sr_x=sub_x, sr_y=sub_y)
1016 else
1017 call calc_field(fieldname, field, xlat_ptr, xlon_ptr, istagger, &
1018 sm1, em1, sm2, em2, min_category, max_category, &
1019 processed_domain, 1, landmask=landmask, sr_x=sub_x, sr_y=sub_y)
1020 end if
1022 ! If user wants to halt when a missing value is found in output field, check now
1023 if (halt_on_missing) then
1024 do i=sm1, em1
1025 do j=sm2, em2
1026 ! Only need to examine k=1
1027 if (field(i,j,1) == msg_fill_val) then
1028 call mprintf(.true.,ERROR,' Missing value encountered in output field. Quitting.')
1029 end if
1030 end do
1031 end do
1032 end if
1034 ! Find fractions
1035 do i=sm1, em1
1036 do j=sm2, em2
1037 sum = 0.
1038 do k=min_category,max_category
1039 sum = sum + field(i,j,k)
1040 end do
1041 if (sum > 0.0) then
1042 do k=min_category,max_category
1043 field(i,j,k) = field(i,j,k) / sum
1044 end do
1045 else
1046 do k=min_category,max_category
1047 field(i,j,k) = msg_fill_val
1048 end do
1049 end if
1050 end do
1051 end do
1053 ! If we should only save the dominant category, then no need to write out fractional field
1054 if (.not.only_save_dominant .or. (idomcatstatus /= 0)) then
1056 ! Finally, we may be asked to smooth the fractional field
1057 call get_smooth_option(fieldname, smth_opt, smth_passes, istatus)
1058 if (istatus == 0) then
1059 if (grid_type == 'C') then
1060 if (smth_opt == ONETWOONE) then
1061 call one_two_one(field, &
1062 start_patch_i, end_patch_i, &
1063 start_patch_j, end_patch_j, &
1064 sm1, em1, &
1065 sm2, em2, &
1066 min_category, max_category, &
1067 smth_passes, msg_fill_val)
1068 else if (smth_opt == SMTHDESMTH) then
1069 call smth_desmth(field, &
1070 start_patch_i, end_patch_i, &
1071 start_patch_j, end_patch_j, &
1072 sm1, em1, &
1073 sm2, em2, &
1074 min_category, max_category, &
1075 smth_passes, msg_fill_val)
1076 else if (smth_opt == SMTHDESMTH_SPECIAL) then
1077 call smth_desmth_special(field, &
1078 start_patch_i, end_patch_i, &
1079 start_patch_j, end_patch_j, &
1080 sm1, em1, &
1081 sm2, em2, &
1082 min_category, max_category, &
1083 smth_passes, msg_fill_val)
1084 end if
1085 else if (grid_type == 'E') then
1086 if (smth_opt == ONETWOONE) then
1087 call one_two_one_egrid(field, &
1088 start_patch_i, end_patch_i, &
1089 start_patch_j, end_patch_j, &
1090 sm1, em1, &
1091 sm2, em2, &
1092 min_category, max_category, &
1093 smth_passes, msg_fill_val, 1.0)
1094 else if (smth_opt == SMTHDESMTH) then
1095 call smth_desmth_egrid(field, &
1096 start_patch_i, end_patch_i, &
1097 start_patch_j, end_patch_j, &
1098 sm1, em1, &
1099 sm2, em2, &
1100 min_category, max_category, &
1101 smth_passes, msg_fill_val, 1.0)
1102 else if (smth_opt == SMTHDESMTH_SPECIAL) then
1103 call mprintf(.true.,WARN,'smth-desmth_special is not currently implemented for NMM. '// &
1104 'No smoothing will be done.')
1105 end if
1106 end if
1107 end if
1109 call write_field(sm1, em1, sm2, em2, &
1110 min_category, max_category, trim(fieldname), datestr, real_array=field)
1111 end if
1113 if (idomcatstatus == 0) then
1114 call mprintf(.true.,STDOUT,' Processing %s', s1=trim(domname))
1115 allocate(dominant_field(sm1:em1, sm2:em2))
1117 if (.not. only_save_dominant) then
1118 field_count = field_count + 1
1119 call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', &
1120 i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=domname)
1121 end if
1123 do i=sm1, em1
1124 do j=sm2, em2
1125 dominant = 0.
1126 dominant_field(i,j) = real(min_category-1)
1127 do k=min_category,max_category
1128 if (field(i,j,k) > dominant .and. field(i,j,k) /= msg_fill_val) then
1129 dominant_field(i,j) = real(k)
1130 dominant = field(i,j,k)
1131 ! else
1132 ! dominant_field(i,j) = nint(msg_fill_val)
1133 ! Maybe we should put an else clause here to set the category equal to the missing fill value?
1134 ! BUG: The problem here seems to be that, when we set a fraction equal to the missing fill value
1135 ! above, if the last fractional index we process here has been filled, we think that the dominant
1136 ! category should be set to the missing fill value. Perhaps we could do some check to only
1137 ! assign the msg_fill_val if no other valid category has been assigned? But this may still not
1138 ! work if the missing fill value is something like 0.5. Somehow use bitarrays, perhaps, to remember
1139 ! which points are missing and which just happen to have the missing fill value?
1140 end if
1141 end do
1142 if (dominant_field(i,j) == real(min_category-1)) dominant_field(i,j) = msg_fill_val
1143 end do
1144 end do
1145 call write_field(sm1, em1, sm2, em2, 1, 1, &
1146 trim(domname), datestr, dominant_field)
1147 deallocate(dominant_field)
1148 end if
1150 deallocate(field)
1152 if ((sub_x > 1) .or. (sub_y > 1)) then
1153 if (associated(xlat_array_subgrid)) deallocate(xlat_array_subgrid)
1154 if (associated(xlon_array_subgrid)) deallocate(xlon_array_subgrid)
1155 if (associated(mapfac_array_x_subgrid)) deallocate(mapfac_array_x_subgrid)
1156 if (associated(mapfac_array_y_subgrid)) deallocate(mapfac_array_y_subgrid)
1157 end if
1159 end if
1161 end if
1162 end if
1164 end do
1166 ! Close output
1167 call output_close()
1169 call hash_destroy(processed_fieldnames)
1171 ! Free up memory
1172 if (associated(xlat_array)) deallocate(xlat_array)
1173 if (associated(xlon_array)) deallocate(xlon_array)
1174 if (grid_type == 'C') then
1175 if (associated(xlat_array_u)) deallocate(xlat_array_u)
1176 if (associated(xlon_array_u)) deallocate(xlon_array_u)
1177 if (associated(xlat_array_corner)) deallocate(xlat_array_corner)
1178 if (associated(xlon_array_corner)) deallocate(xlon_array_corner)
1179 if (associated(mapfac_array_u_x)) deallocate(mapfac_array_u_x)
1180 if (associated(mapfac_array_u_y)) deallocate(mapfac_array_u_y)
1181 end if
1182 if (associated(xlat_array_v)) deallocate(xlat_array_v)
1183 if (associated(xlon_array_v)) deallocate(xlon_array_v)
1184 if (associated(mapfac_array_m_x)) deallocate(mapfac_array_m_x)
1185 if (associated(mapfac_array_m_y)) deallocate(mapfac_array_m_y)
1186 if (associated(mapfac_array_v_x)) deallocate(mapfac_array_v_x)
1187 if (associated(mapfac_array_v_y)) deallocate(mapfac_array_v_y)
1188 if (associated(landmask)) deallocate(landmask)
1189 if (associated(xlat_array_subgrid)) deallocate(xlat_array_subgrid)
1190 if (associated(xlon_array_subgrid)) deallocate(xlon_array_subgrid)
1191 if (associated(mapfac_array_x_subgrid)) deallocate(mapfac_array_x_subgrid)
1192 if (associated(mapfac_array_y_subgrid)) deallocate(mapfac_array_y_subgrid)
1194 nullify(xlat_ptr)
1195 nullify(xlon_ptr)
1197 end subroutine process_tile
1200 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1201 ! Name: calc_field
1202 !
1203 ! Purpose: This routine fills in the "field" array with interpolated source
1204 ! data. When multiple resolutions of source data are available, an appropriate
1205 ! resolution is chosen automatically. The specified field may either be a
1206 ! continuous field or a categorical field.
1207 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1208 recursive subroutine calc_field(fieldname, field, xlat_array, xlon_array, istagger, &
1209 start_i, end_i, start_j, end_j, start_k, end_k, &
1210 processed_domain, ilevel, landmask, sr_x, sr_y)
1212 use bitarray_module
1213 use interp_module
1214 use llxy_module
1215 use misc_definitions_module
1216 use proc_point_module
1217 use queue_module
1218 use source_data_module
1220 implicit none
1222 ! Arguments
1223 integer, intent(in) :: start_i, end_i, start_j, end_j, start_k, end_k, ilevel, istagger
1224 real, dimension(start_i:end_i, start_j:end_j), intent(in) :: xlat_array, xlon_array
1225 real, dimension(start_i:end_i, start_j:end_j, start_k:end_k), intent(inout) :: field
1226 real, dimension(start_i:end_i, start_j:end_j), intent(in), optional :: landmask
1227 integer, intent(in), optional :: sr_x, sr_y
1228 character (len=128), intent(in) :: fieldname
1229 type (bitarray), intent(inout) :: processed_domain
1231 ! Local variables
1232 integer :: start_src_k, end_src_k
1233 integer :: i, j, k, ix, iy, itype
1234 integer :: user_iproj, istatus
1235 integer :: opt_status
1236 real :: mask_val
1237 real :: temp
1238 real :: scale_factor
1239 real :: msg_val, msg_fill_val, threshold, src_dx, src_dy, dom_dx, dom_dy
1240 real :: user_stand_lon, user_truelat1, user_truelat2, user_dxkm, user_dykm, &
1241 user_known_x, user_known_y, user_known_lat, user_known_lon
1242 real, pointer, dimension(:,:,:) :: data_count
1243 integer, pointer, dimension(:) :: interp_type
1244 integer, pointer, dimension(:) :: interp_opts
1245 character (len=128) :: interp_string
1246 type (bitarray) :: bit_domain, level_domain
1247 type (queue) :: point_queue, tile_queue
1248 type (q_data) :: current_pt
1250 nullify(data_count)
1251 nullify(interp_type)
1252 nullify(interp_opts)
1254 ! If this is the first trip through this routine, we need to allocate the bit array that
1255 ! will persist through all recursive calls, tracking which grid points have been assigned
1256 ! a value.
1257 if (ilevel == 1) call bitarray_create(processed_domain, end_i-start_i+1, end_j-start_j+1)
1259 ! Find out if this "priority level" (given by ilevel) exists
1260 call check_priority_level(fieldname, ilevel, istatus)
1262 ! A bad status indicates that that no data for priority level ilevel is available, and thus, that
1263 ! no further levels will be specified. We are done processing for this level.
1264 if (istatus /= 0) then
1265 if (ilevel == 1) call bitarray_destroy(processed_domain)
1266 return
1267 end if
1269 ! Before proceeding with processing for this level, though, process for the next highest priority level
1270 ! of source data
1271 call calc_field(fieldname, field, xlat_array, xlon_array, istagger, start_i, end_i, &
1272 start_j, end_j, start_k, end_k, processed_domain, ilevel+1, landmask, sr_x, sr_y)
1274 ! At this point, all levels of source data with higher priority have been processed, and we can assign
1275 ! values to all grid points that have not already been given values from higher-priority data
1277 call get_source_opt_status(fieldname, ilevel, opt_status)
1278 if (opt_status == 0) then
1280 ! Find out the projection of the data for this "priority level" (given by ilevel)
1281 call get_data_projection(fieldname, user_iproj, user_stand_lon, user_truelat1, user_truelat2, &
1282 user_dxkm, user_dykm, user_known_x, user_known_y, user_known_lat, &
1283 user_known_lon, ilevel, istatus)
1285 ! A good status indicates that there is data for this priority level, so we store the projection
1286 ! of that data on a stack. The projection will be on the top of the stack (and hence will be
1287 ! the "active" projection) once all higher-priority levels have been processed
1288 call push_source_projection(user_iproj, user_stand_lon, user_truelat1, user_truelat2, &
1289 user_dxkm, user_dykm, user_dykm, user_dxkm, user_known_x, user_known_y, &
1290 user_known_lat, user_known_lon)
1292 ! Initialize point processing module
1293 call proc_point_init()
1295 ! Initialize queues
1296 call q_init(point_queue)
1297 call q_init(tile_queue)
1299 ! Determine whether we will be processing categorical data or continuous data
1300 itype = iget_source_fieldtype(fieldname, ilevel, istatus)
1301 call get_interp_option(fieldname, ilevel, interp_string, istatus)
1302 interp_type => interp_array_from_string(interp_string)
1303 interp_opts => interp_options_from_string(interp_string)
1305 ! Also, check whether we will be using the cell averaging interpolator for continuous fields
1306 if (index(interp_string,'average_gcell') /= 0 .and. itype == CONTINUOUS) then
1307 call get_gcell_threshold(interp_string, threshold, istatus)
1308 if (istatus == 0) then
1309 call get_source_resolution(fieldname, ilevel, src_dx, src_dy, istatus)
1310 if (istatus == 0) then
1311 call get_domain_resolution(dom_dx, dom_dy)
1312 if (gridtype == 'C') then
1313 if (threshold*max(src_dx,src_dy)*111. <= max(dom_dx,dom_dy)/1000.) then
1314 itype = SP_CONTINUOUS
1315 allocate(data_count(start_i:end_i,start_j:end_j,start_k:end_k))
1316 data_count = 0.
1317 end if
1318 else if (gridtype == 'E') then
1319 if (max(src_dx,src_dy) >= threshold*max(dom_dx,dom_dy)) then
1320 itype = SP_CONTINUOUS
1321 allocate(data_count(start_i:end_i,start_j:end_j,start_k:end_k))
1322 data_count = 0.
1323 end if
1324 end if
1325 end if
1326 end if
1327 end if
1329 call get_missing_value(fieldname, ilevel, msg_val, istatus)
1330 if (istatus /= 0) msg_val = NAN
1331 call get_missing_fill_value(fieldname, msg_fill_val, istatus)
1332 if (istatus /= 0) msg_fill_val = NAN
1334 call get_masked_value(fieldname, ilevel, mask_val, istatus)
1335 if (istatus /= 0) mask_val = -1.
1337 if (itype == CONTINUOUS .or. itype == SP_CONTINUOUS) then
1338 call get_source_levels(fieldname, ilevel, start_src_k, end_src_k, istatus)
1339 if (istatus /= 0) return
1340 end if
1342 ! Initialize bitarray used to track which points have been visited and assigned values while
1343 ! processing *this* priority level of data
1344 call bitarray_create(bit_domain, end_i-start_i+1, end_j-start_j+1)
1345 call bitarray_create(level_domain, end_i-start_i+1, end_j-start_j+1)
1347 ! Begin by placing a point in the tile_queue
1348 current_pt%lat = xlat_array(start_i,start_j)
1349 current_pt%lon = xlon_array(start_i,start_j)
1350 current_pt%x = start_i
1351 current_pt%y = start_j
1352 call q_insert(tile_queue, current_pt)
1354 ! While there are still grid points in tiles that have not yet been processed
1355 do while (q_isdata(tile_queue))
1357 ! Take a point from the outer queue and place it in the point_queue for processing
1358 current_pt = q_remove(tile_queue)
1360 ! If this level of data is categorical (i.e., is given as an array of category indices),
1361 ! then first try to process the entire tile in one call to accum_categorical. Any grid
1362 ! points that are not given values by accum_categorical and that lie within the current
1363 ! tile of source data are individually assigned values in the inner loop
1364 if (itype == CATEGORICAL) then
1366 ! Have we already visited this point? If so, this tile has already been processed by
1367 ! accum_categorical.
1368 if (.not. bitarray_test(bit_domain, current_pt%x-start_i+1, current_pt%y-start_j+1)) then
1369 call q_insert(point_queue, current_pt)
1370 if (.not. have_processed_tile(current_pt%lat, current_pt%lon, fieldname, ilevel)) then
1371 call accum_categorical(current_pt%lat, current_pt%lon, istagger, field, &
1372 start_i, end_i, start_j, end_j, start_k, end_k, &
1373 fieldname, processed_domain, level_domain, &
1374 ilevel, msg_val, mask_val, sr_x, sr_y)
1375 ! BUG: Where do we mask out those points that are on land/water when masked=land/water is set?
1376 end if
1377 call bitarray_set(bit_domain, current_pt%x-start_i+1, current_pt%y-start_j+1)
1378 end if
1380 else if (itype == SP_CONTINUOUS) then
1382 ! Have we already visited this point? If so, this tile has already been processed by
1383 ! accum_continuous.
1384 if (.not. bitarray_test(bit_domain, current_pt%x-start_i+1, current_pt%y-start_j+1)) then
1385 call q_insert(point_queue, current_pt)
1386 if (.not. have_processed_tile(current_pt%lat, current_pt%lon, fieldname, ilevel)) then
1387 call accum_continuous(current_pt%lat, current_pt%lon, istagger, field, data_count, &
1388 start_i, end_i, start_j, end_j, start_k, end_k, &
1389 fieldname, processed_domain, level_domain, &
1390 ilevel, msg_val, mask_val, sr_x, sr_y)
1391 ! BUG: Where do we mask out those points that are on land/water when masked=land/water is set?
1392 end if
1393 call bitarray_set(bit_domain, current_pt%x-start_i+1, current_pt%y-start_j+1)
1394 end if
1396 else if (itype == CONTINUOUS) then
1398 ! Have we already visited this point? If so, the tile containing this point has already been
1399 ! processed in the inner loop.
1400 if (.not. bitarray_test(bit_domain, current_pt%x-start_i+1, current_pt%y-start_j+1)) then
1401 call q_insert(point_queue, current_pt)
1402 call bitarray_set(bit_domain, current_pt%x-start_i+1, current_pt%y-start_j+1)
1403 end if
1405 end if
1407 ! This inner loop, where all grid points contained in the current source tile are processed
1408 do while (q_isdata(point_queue))
1409 current_pt = q_remove(point_queue)
1410 ix = current_pt%x
1411 iy = current_pt%y
1413 ! Process the current point
1414 if (itype == CONTINUOUS .or. itype == SP_CONTINUOUS) then
1416 ! Have we already assigned this point a value from this priority level?
1417 if (.not. bitarray_test(level_domain, ix-start_i+1, iy-start_j+1)) then
1419 ! If the point was already assigned a value from a higher-priority level, no
1420 ! need to assign a new value
1421 if (bitarray_test(processed_domain, ix-start_i+1, iy-start_j+1)) then
1422 call bitarray_set(level_domain, ix-start_i+1, iy-start_j+1)
1424 ! Otherwise, need to assign values from this level of source data if we can
1425 else
1426 if (present(landmask) .and. (istagger == M .or. istagger == HH)) then
1427 if (landmask(ix,iy) /= mask_val) then
1428 do k=start_src_k,end_src_k
1429 temp = get_point(current_pt%lat, current_pt%lon, k, &
1430 fieldname, ilevel, interp_type, interp_opts, msg_val)
1431 if (temp /= msg_val) then
1432 field(ix, iy, k) = temp
1433 call bitarray_set(level_domain, ix-start_i+1, iy-start_j+1)
1434 if (itype == SP_CONTINUOUS) data_count(ix, iy, k) = 1.0
1435 else
1436 field(ix, iy, k) = msg_fill_val
1437 end if
1438 end do
1439 else
1440 do k=start_k,end_k
1441 field(ix,iy,k) = msg_fill_val
1442 end do
1443 end if
1444 else
1445 do k=start_src_k,end_src_k
1446 temp = get_point(current_pt%lat, current_pt%lon, k, &
1447 fieldname, ilevel, interp_type, interp_opts, msg_val)
1448 if (temp /= msg_val) then
1449 field(ix, iy, k) = temp
1450 call bitarray_set(level_domain, ix-start_i+1, iy-start_j+1)
1451 if (itype == SP_CONTINUOUS) data_count(ix, iy, k) = 1.0
1452 else
1453 field(ix, iy, k) = msg_fill_val
1454 end if
1455 end do
1456 end if
1457 end if
1458 end if
1460 else if (itype == CATEGORICAL) then
1462 ! Have we already assigned this point a value from this priority level?
1463 if (.not.bitarray_test(level_domain, ix-start_i+1, iy-start_j+1)) then
1465 ! If the point was already assigned a value from a higher-priority level, no
1466 ! need to assign a new value
1467 if (bitarray_test(processed_domain, ix-start_i+1, iy-start_j+1)) then
1468 call bitarray_set(level_domain, ix-start_i+1, iy-start_j+1)
1470 ! Otherwise, the point was apparently not given a value when accum_categorical
1471 ! was called for the current tile, and we need to assign values from this
1472 ! level of source data if we can
1473 else
1474 if (present(landmask) .and. (istagger == M .or. istagger == HH)) then
1475 if (landmask(ix,iy) /= mask_val) then
1476 temp = get_point(current_pt%lat, current_pt%lon, 1, &
1477 fieldname, ilevel, interp_type, interp_opts, msg_val)
1479 do k=start_k,end_k
1480 field(ix,iy,k) = 0.
1481 end do
1483 if (temp /= msg_val) then
1484 if (int(temp) >= start_k .and. int(temp) <= end_k) then
1485 field(ix, iy, int(temp)) = field(ix, iy, int(temp)) + 1.
1486 else
1487 call mprintf(.true.,WARN,' Attempted to assign an invalid category '// &
1488 '%i to grid point (%i, %i)', i1=int(temp), i2=ix, i3=iy)
1489 end if
1490 call bitarray_set(level_domain, ix-start_i+1, iy-start_j+1)
1491 end if
1493 else
1494 do k=start_k,end_k
1495 field(ix,iy,k) = 0.
1496 end do
1497 end if
1498 else
1499 temp = get_point(current_pt%lat, current_pt%lon, 1, &
1500 fieldname, ilevel, interp_type, interp_opts, msg_val)
1502 do k=start_k,end_k
1503 field(ix,iy,k) = 0.
1504 end do
1506 if (temp /= msg_val) then
1507 if (int(temp) >= start_k .and. int(temp) <= end_k) then
1508 field(ix, iy, int(temp)) = field(ix, iy, int(temp)) + 1.
1509 else
1510 call mprintf(.true.,WARN,' Attempted to assign an invalid category '// &
1511 '%i to grid point (%i, %i)', i1=int(temp), i2=ix, i3=iy)
1512 end if
1513 call bitarray_set(level_domain, ix-start_i+1, iy-start_j+1)
1514 end if
1515 end if
1516 end if
1517 end if
1519 end if
1521 ! Scan neighboring points, adding them to the appropriate queue based on whether they
1522 ! are in the current tile or not
1523 if (iy > start_j) then
1524 if (ix > start_i) then
1526 ! Neighbor with relative position (-1,-1)
1527 call process_neighbor(ix-1, iy-1, bit_domain, point_queue, tile_queue, &
1528 xlat_array, xlon_array, start_i, end_i, start_j, end_j, ilevel)
1529 end if
1531 ! Neighbor with relative position (0,-1)
1532 call process_neighbor(ix, iy-1, bit_domain, point_queue, tile_queue, &
1533 xlat_array, xlon_array, start_i, end_i, start_j, end_j, ilevel)
1535 if (ix < end_i) then
1537 ! Neighbor with relative position (+1,-1)
1538 call process_neighbor(ix+1, iy-1, bit_domain, point_queue, tile_queue, &
1539 xlat_array, xlon_array, start_i, end_i, start_j, end_j, ilevel)
1540 end if
1541 end if
1543 if (ix > start_i) then
1545 ! Neighbor with relative position (-1,0)
1546 call process_neighbor(ix-1, iy, bit_domain, point_queue, tile_queue, &
1547 xlat_array, xlon_array, start_i, end_i, start_j, end_j, ilevel)
1548 end if
1550 if (ix < end_i) then
1552 ! Neighbor with relative position (+1,0)
1553 call process_neighbor(ix+1, iy, bit_domain, point_queue, tile_queue, &
1554 xlat_array, xlon_array, start_i, end_i, start_j, end_j, ilevel)
1555 end if
1557 if (iy < end_j) then
1558 if (ix > start_i) then
1560 ! Neighbor with relative position (-1,+1)
1561 call process_neighbor(ix-1, iy+1, bit_domain, point_queue, tile_queue, &
1562 xlat_array, xlon_array, start_i, end_i, start_j, end_j, ilevel)
1563 end if
1565 ! Neighbor with relative position (0,+1)
1566 call process_neighbor(ix, iy+1, bit_domain, point_queue, tile_queue, &
1567 xlat_array, xlon_array, start_i, end_i, start_j, end_j, ilevel)
1568 if (ix < end_i) then
1570 ! Neighbor with relative position (+1,+1)
1571 call process_neighbor(ix+1, iy+1, bit_domain, point_queue, tile_queue, &
1572 xlat_array, xlon_array, start_i, end_i, start_j, end_j, ilevel)
1573 end if
1574 end if
1576 end do
1578 end do
1580 if (itype == SP_CONTINUOUS) then
1581 itype = CONTINUOUS
1582 if (present(landmask) .and. (istagger == M .or. istagger == HH)) then
1583 do j=start_j,end_j
1584 do i=start_i,end_i
1585 if (landmask(i,j) /= mask_val) then
1586 do k=start_k,end_k
1587 if (data_count(i,j,k) > 0.) then
1588 field(i,j,k) = field(i,j,k) / data_count(i,j,k)
1589 else
1590 if (.not.bitarray_test(processed_domain, i-start_i+1, j-start_j+1)) then
1591 field(i,j,k) = msg_fill_val
1592 end if
1593 end if
1594 end do
1595 else
1596 if (.not.bitarray_test(processed_domain, i-start_i+1, j-start_j+1)) then
1597 do k=start_k,end_k
1598 field(i,j,k) = msg_fill_val
1599 end do
1600 end if
1601 end if
1602 end do
1603 end do
1604 else
1605 do k=start_k,end_k
1606 do j=start_j,end_j
1607 do i=start_i,end_i
1608 if (data_count(i,j,k) > 0.) then
1609 field(i,j,k) = field(i,j,k) / data_count(i,j,k)
1610 else
1611 if (.not.bitarray_test(processed_domain, i-start_i+1, j-start_j+1)) then
1612 field(i,j,k) = msg_fill_val
1613 end if
1614 end if
1615 end do
1616 end do
1617 end do
1618 end if
1619 deallocate(data_count)
1621 else if (itype == CATEGORICAL) then
1622 if (present(landmask) .and. (istagger == M .or. istagger == HH)) then
1623 do j=start_j,end_j
1624 do i=start_i,end_i
1625 if (landmask(i,j) == mask_val) then
1626 do k=start_k,end_k
1627 field(i,j,k) = 0.
1628 end do
1629 end if
1630 end do
1631 end do
1632 end if
1633 end if
1635 deallocate(interp_type)
1636 deallocate(interp_opts)
1639 ! We may need to scale this field by a constant
1640 call get_field_scale_factor(fieldname, ilevel, scale_factor, istatus)
1641 if (istatus == 0) then
1642 do i=start_i, end_i
1643 do j=start_j, end_j
1644 if (bitarray_test(level_domain,i-start_i+1,j-start_j+1) .and. &
1645 .not. bitarray_test(processed_domain,i-start_i+1,j-start_j+1)) then
1646 do k=start_k,end_k
1647 if (field(i,j,k) /= msg_fill_val) then
1648 field(i,j,k) = field(i,j,k) * scale_factor
1649 end if
1650 end do
1651 end if
1652 end do
1653 end do
1654 end if
1657 ! Now add the points that were assigned values at this priority level to the complete array
1658 ! of points that have been assigned values
1659 call bitarray_merge(processed_domain, level_domain)
1661 call bitarray_destroy(bit_domain)
1662 call bitarray_destroy(level_domain)
1663 call q_destroy(point_queue)
1664 call q_destroy(tile_queue)
1665 call proc_point_shutdown()
1667 ! Remove the projection of the current level of source data from the stack, thus "activating"
1668 ! the projection of the next highest level
1669 call pop_source_projection()
1671 else
1672 call mprintf(.true.,STDOUT,' Important note: could not open input dataset for priority level %i, '// &
1673 'but this level is optional.', i1=ilevel)
1674 call mprintf(.true.,LOGFILE,' Important note: could not open input dataset for priority level %i, '// &
1675 'but this level is optional.', i1=ilevel)
1676 end if
1678 ! If this is the last level of the recursion, we can also deallocate processed_domain
1679 if (ilevel == 1) call bitarray_destroy(processed_domain)
1681 end subroutine calc_field
1684 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1685 ! Name: get_lat_lon_fields
1686 !
1687 ! Purpose: To calculate the latitude and longitude for every gridpoint in the
1688 ! tile of the model domain. The caller may specify that the grid for which
1689 ! values are computed is staggered or unstaggered using the "stagger"
1690 ! argument.
1691 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1692 subroutine get_lat_lon_fields(xlat_arr, xlon_arr, start_mem_i, &
1693 start_mem_j, end_mem_i, end_mem_j, stagger, comp_ll, &
1694 sub_x, sub_y)
1696 use llxy_module
1697 use misc_definitions_module
1699 implicit none
1701 ! Arguments
1702 integer, intent(in) :: start_mem_i, start_mem_j, end_mem_i, &
1703 end_mem_j, stagger
1704 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(out) :: xlat_arr, xlon_arr
1705 logical, optional, intent(in) :: comp_ll
1706 integer, optional, intent(in) :: sub_x, sub_y
1708 ! Local variables
1709 integer :: i, j
1710 real :: rx, ry
1712 rx = 1.0
1713 ry = 1.0
1714 if (present(sub_x)) rx = real(sub_x)
1715 if (present(sub_y)) ry = real(sub_y)
1717 do i=start_mem_i, end_mem_i
1718 do j=start_mem_j, end_mem_j
1719 call xytoll(real(i-0.5)/rx+0.5, real(j-0.5)/ry+0.5, &
1720 xlat_arr(i,j), xlon_arr(i,j), stagger, comp_ll=comp_ll)
1721 end do
1722 end do
1724 end subroutine get_lat_lon_fields
1727 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1728 ! Name: get_map_factor
1729 !
1730 ! Purpose: Given the latitude field, this routine calculates map factors for
1731 ! the grid points of the specified domain. For different grids (e.g., C grid,
1732 ! E grid), the latitude array should provide the latitudes of the points for
1733 ! which map factors are to be calculated.
1734 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1735 subroutine get_map_factor(xlat_arr, xlon_arr, mapfac_arr_x, mapfac_arr_y, &
1736 start_mem_i, start_mem_j, end_mem_i, end_mem_j)
1738 use constants_module
1739 use gridinfo_module
1740 use misc_definitions_module
1741 use map_utils
1743 implicit none
1745 ! Arguments
1746 integer, intent(in) :: start_mem_i, start_mem_j, end_mem_i, end_mem_j
1747 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(in) :: xlat_arr, xlon_arr
1748 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(out) :: mapfac_arr_x
1749 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(out) :: mapfac_arr_y
1751 ! Local variables
1752 integer :: i, j
1753 real :: n, colat, colat0, colat1, colat2, comp_lat, comp_lon
1755 !
1756 ! Equations for map factor given in Principles of Meteorological Analysis,
1757 ! Walter J. Saucier, pp. 32-33
1758 !
1760 ! Lambert conformal projection
1761 if (iproj_type == PROJ_LC) then
1762 if (truelat1 /= truelat2) then
1763 colat1 = rad_per_deg*(90.0 - truelat1)
1764 colat2 = rad_per_deg*(90.0 - truelat2)
1765 n = (log(sin(colat1)) - log(sin(colat2))) &
1766 / (log(tan(colat1/2.0)) - log(tan(colat2/2.0)))
1768 do i=start_mem_i, end_mem_i
1769 do j=start_mem_j, end_mem_j
1770 colat = rad_per_deg*(90.0 - xlat_arr(i,j))
1771 mapfac_arr_x(i,j) = sin(colat2)/sin(colat)*(tan(colat/2.0)/tan(colat2/2.0))**n
1772 mapfac_arr_y(i,j) = mapfac_arr_x(i,j)
1773 end do
1774 end do
1776 else
1777 colat0 = rad_per_deg*(90.0 - truelat1)
1779 do i=start_mem_i, end_mem_i
1780 do j=start_mem_j, end_mem_j
1781 colat = rad_per_deg*(90.0 - xlat_arr(i,j))
1782 mapfac_arr_x(i,j) = sin(colat0)/sin(colat)*(tan(colat/2.0)/tan(colat0/2.0))**cos(colat0)
1783 mapfac_arr_y(i,j) = mapfac_arr_x(i,j)
1784 end do
1785 end do
1787 end if
1789 ! Polar stereographic projection
1790 else if (iproj_type == PROJ_PS) then
1792 do i=start_mem_i, end_mem_i
1793 do j=start_mem_j, end_mem_j
1794 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)))
1795 mapfac_arr_y(i,j) = mapfac_arr_x(i,j)
1796 end do
1797 end do
1799 ! Mercator projection
1800 else if (iproj_type == PROJ_MERC) then
1801 colat0 = rad_per_deg*(90.0 - truelat1)
1803 do i=start_mem_i, end_mem_i
1804 do j=start_mem_j, end_mem_j
1805 colat = rad_per_deg*(90.0 - xlat_arr(i,j))
1806 mapfac_arr_x(i,j) = sin(colat0) / sin(colat)
1807 mapfac_arr_y(i,j) = mapfac_arr_x(i,j)
1808 end do
1809 end do
1811 ! Global cylindrical projection
1812 else if (iproj_type == PROJ_CYL) then
1814 do i=start_mem_i, end_mem_i
1815 do j=start_mem_j, end_mem_j
1816 if (abs(xlat_arr(i,j)) == 90.0) then
1817 mapfac_arr_x(i,j) = 0. ! MSF actually becomes infinite at poles, but
1818 ! the values should never be used there; by
1819 ! setting to 0, we hope to induce a "divide
1820 ! by zero" error if they are
1821 else
1822 mapfac_arr_x(i,j) = 1.0 / cos(xlat_arr(i,j)*rad_per_deg)
1823 end if
1824 mapfac_arr_y(i,j) = 1.0
1825 end do
1826 end do
1828 ! Rotated global cylindrical projection
1829 else if (iproj_type == PROJ_CASSINI) then
1831 if (abs(pole_lat) == 90.) then
1832 do i=start_mem_i, end_mem_i
1833 do j=start_mem_j, end_mem_j
1834 if (abs(xlat_arr(i,j)) >= 90.0) then
1835 mapfac_arr_x(i,j) = 0. ! MSF actually becomes infinite at poles, but
1836 ! the values should never be used there; by
1837 ! setting to 0, we hope to induce a "divide
1838 ! by zero" error if they are
1839 else
1840 mapfac_arr_x(i,j) = 1.0 / cos(xlat_arr(i,j)*rad_per_deg)
1841 end if
1842 mapfac_arr_y(i,j) = 1.0
1843 end do
1844 end do
1845 else
1846 do i=start_mem_i, end_mem_i
1847 do j=start_mem_j, end_mem_j
1848 call rotate_coords(xlat_arr(i,j),xlon_arr(i,j), &
1849 comp_lat, comp_lon, &
1850 pole_lat, pole_lon, stand_lon, &
1851 -1)
1852 if (abs(comp_lat) >= 90.0) then
1853 mapfac_arr_x(i,j) = 0. ! MSF actually becomes infinite at poles, but
1854 ! the values should never be used there; by
1855 ! setting to 0, we hope to induce a "divide
1856 ! by zero" error if they are
1857 else
1858 mapfac_arr_x(i,j) = 1.0 / cos(comp_lat*rad_per_deg)
1859 end if
1860 mapfac_arr_y(i,j) = 1.0
1861 end do
1862 end do
1863 end if
1865 else if (iproj_type == PROJ_ROTLL) then
1867 do i=start_mem_i, end_mem_i
1868 do j=start_mem_j, end_mem_j
1869 mapfac_arr_x(i,j) = 1.0
1870 mapfac_arr_y(i,j) = 1.0
1871 end do
1872 end do
1874 end if
1876 end subroutine get_map_factor
1879 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1880 ! Name: get_coriolis_parameters
1881 !
1882 ! Purpose: To calculate the Coriolis parameters f and e for every gridpoint in
1883 ! the tile of the model domain
1884 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1885 subroutine get_coriolis_parameters(xlat_arr, f, e, &
1886 start_mem_i, start_mem_j, end_mem_i, end_mem_j)
1888 use constants_module
1890 implicit none
1892 ! Arguments
1893 integer, intent(in) :: start_mem_i, start_mem_j, end_mem_i, end_mem_j
1894 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(in) :: xlat_arr
1895 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(out) :: f, e
1897 ! Local variables
1898 integer :: i, j
1900 do i=start_mem_i, end_mem_i
1901 do j=start_mem_j, end_mem_j
1903 f(i,j) = 2.0*OMEGA_E*sin(rad_per_deg*xlat_arr(i,j))
1904 e(i,j) = 2.0*OMEGA_E*cos(rad_per_deg*xlat_arr(i,j))
1906 end do
1907 end do
1909 end subroutine get_coriolis_parameters
1912 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1913 ! Name: get_rotang
1914 !
1915 ! Purpose: To calculate the sine and cosine of rotation angle.
1916 !
1917 ! NOTES: The formulas used in this routine come from those in the
1918 ! vecrot_rotlat() routine of the original WRF SI.
1919 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1920 subroutine get_rotang(xlat_arr, xlon_arr, cosa, sina, &
1921 start_mem_i, start_mem_j, end_mem_i, end_mem_j)
1923 use constants_module
1924 use gridinfo_module
1925 use map_utils, only : lc_cone
1927 implicit none
1929 ! Arguments
1930 integer, intent(in) :: start_mem_i, start_mem_j, end_mem_i, end_mem_j
1931 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(in) :: xlat_arr, xlon_arr
1932 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(out) :: cosa, sina
1934 ! Local variables
1935 integer :: i, j
1936 real :: cone
1937 real :: alpha, d_lon
1939 !
1940 ! Lambert conformal projection
1941 !
1942 if (iproj_type == PROJ_LC) then
1944 ! Rotation angle is given by the difference in longitude from the standard longitude
1945 ! scaled by the cone factor of the projection
1947 call lc_cone(truelat1, truelat2, cone)
1949 do i=start_mem_i, end_mem_i
1950 do j=start_mem_j, end_mem_j
1951 d_lon = stand_lon - xlon_arr(i,j)
1952 if (d_lon > 180.) then
1953 d_lon = d_lon - 360.
1954 else if (d_lon < -180.) then
1955 d_lon = d_lon + 360.
1956 end if
1958 alpha = d_lon * cone * RAD_PER_DEG
1960 sina(i,j) = sin(alpha)
1961 cosa(i,j) = cos(alpha)
1962 end do
1963 end do
1965 !
1966 ! Polar stereographic projection
1967 !
1968 else if (iproj_type == PROJ_PS) then
1970 ! Rotation angle is given by the difference in longitude from the standard longitude
1972 do i=start_mem_i, end_mem_i
1973 do j=start_mem_j, end_mem_j
1974 d_lon = stand_lon - xlon_arr(i,j)
1975 if (d_lon > 180.) then
1976 d_lon = d_lon - 360.
1977 else if (d_lon < -180.) then
1978 d_lon = d_lon + 360.
1979 end if
1981 alpha = d_lon * RAD_PER_DEG
1983 sina(i,j) = sin(alpha)
1984 cosa(i,j) = cos(alpha)
1985 end do
1986 end do
1988 !
1989 ! Mercator projection
1990 !
1991 else if (iproj_type == PROJ_MERC) then
1993 ! Rotation angle is always zero
1995 sina(:,:) = 0.0
1996 cosa(:,:) = 1.0
1998 !
1999 ! Global cylindrical projection
2000 !
2001 else if (iproj_type == PROJ_CYL) then
2003 ! Rotation angle is always zero
2005 sina(:,:) = 0.0
2006 cosa(:,:) = 1.0
2008 !
2009 ! Rotated global cylindrical projection
2010 !
2011 else if (iproj_type == PROJ_CASSINI) then
2013 ! Compute angles using finite differences
2015 do i=start_mem_i, end_mem_i
2016 do j=start_mem_j+1, end_mem_j-1
2017 d_lon = xlon_arr(i,j+1)-xlon_arr(i,j-1)
2018 if (d_lon > 180.) then
2019 d_lon = d_lon - 360.
2020 else if (d_lon < -180.) then
2021 d_lon = d_lon + 360.
2022 end if
2024 alpha = atan2(-cos(xlat_arr(i,j)*RAD_PER_DEG) * (d_lon*RAD_PER_DEG), &
2025 ((xlat_arr(i,j+1)-xlat_arr(i,j-1))*RAD_PER_DEG))
2026 sina(i,j) = sin(alpha)
2027 cosa(i,j) = cos(alpha)
2028 end do
2029 end do
2031 do i=start_mem_i, end_mem_i
2032 d_lon = xlon_arr(i,start_mem_j+1)-xlon_arr(i,start_mem_j)
2033 if (d_lon > 180.) then
2034 d_lon = d_lon - 360.
2035 else if (d_lon < -180.) then
2036 d_lon = d_lon + 360.
2037 end if
2039 alpha = atan2(-cos(xlat_arr(i,start_mem_j)*RAD_PER_DEG) * (d_lon*RAD_PER_DEG), &
2040 ((xlat_arr(i,start_mem_j+1)-xlat_arr(i,start_mem_j))*RAD_PER_DEG))
2041 sina(i,start_mem_j) = sin(alpha)
2042 cosa(i,start_mem_j) = cos(alpha)
2043 end do
2045 do i=start_mem_i, end_mem_i
2046 d_lon = xlon_arr(i,end_mem_j)-xlon_arr(i,end_mem_j-1)
2047 if (d_lon > 180.) then
2048 d_lon = d_lon - 360.
2049 else if (d_lon < -180.) then
2050 d_lon = d_lon + 360.
2051 end if
2053 alpha = atan2(-cos(xlat_arr(i,end_mem_j)*RAD_PER_DEG) * (d_lon*RAD_PER_DEG), &
2054 ((xlat_arr(i,end_mem_j)-xlat_arr(i,end_mem_j-1))*RAD_PER_DEG))
2055 sina(i,end_mem_j) = sin(alpha)
2056 cosa(i,end_mem_j) = cos(alpha)
2057 end do
2059 ! NMM Rotated latitude-longitude projection
2060 else if (iproj_type == PROJ_ROTLL) then
2062 ! No longer supported...
2063 call mprintf(.true.,ERROR,' The NMM rotated lat-lon grid is no longer supported. Rotation angles cannot be computed.')
2065 end if
2067 end subroutine get_rotang
2070 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2071 ! Name: process_neighbor
2072 !
2073 ! Purpose: This routine, give the x/y location of a point, determines whether
2074 ! the point has already been processed, and if not, which processing queue
2075 ! the point should be placed in.
2076 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2077 subroutine process_neighbor(ix, iy, bit_domain, point_queue, tile_queue, &
2078 xlat_array, xlon_array, &
2079 start_i, end_i, start_j, end_j, ilevel)
2081 use bitarray_module
2082 use misc_definitions_module
2083 use proc_point_module
2084 use queue_module
2086 implicit none
2088 ! Arguments
2089 integer, intent(in) :: ix, iy, start_i, end_i, start_j, end_j, ilevel
2090 real, dimension(start_i:end_i, start_j:end_j), intent(in) :: xlat_array, xlon_array
2091 type (bitarray), intent(inout) :: bit_domain
2092 type (queue), intent(inout) :: point_queue, tile_queue
2094 ! Local variables
2095 type (q_data) :: process_pt
2096 logical :: is_in_tile
2098 ! If the point has already been visited, no need to do anything more.
2099 if (.not. bitarray_test(bit_domain, ix-start_i+1, iy-start_j+1)) then
2101 ! Create a queue item for the current point
2102 process_pt%lat = xlat_array(ix,iy)
2103 process_pt%lon = xlon_array(ix,iy)
2104 process_pt%x = ix
2105 process_pt%y = iy
2107 is_in_tile = is_point_in_tile(process_pt%lat, process_pt%lon, ilevel)
2109 ! If the point is in the current tile, add it to the list of points
2110 ! to be processed in the inner loop
2111 if (is_in_tile) then
2112 call q_insert(point_queue, process_pt)
2113 call bitarray_set(bit_domain, ix-start_i+1, iy-start_j+1)
2115 ! Otherwise, we will process this point later. Add it to the list for
2116 ! the outer loop.
2117 else
2118 call q_insert(tile_queue, process_pt)
2119 end if
2121 end if
2123 end subroutine process_neighbor
2126 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2127 ! Name: calc_dfdy
2128 !
2129 ! Purpose: This routine calculates df/dy for the field in src_arr, and places
2130 ! the result in dst_array.
2131 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2132 subroutine calc_dfdy(src_arr, dst_arr, start_mem_i, start_mem_j, start_mem_k, &
2133 end_mem_i, end_mem_j, end_mem_k, sr_y, mapfac)
2135 ! Modules
2136 use llxy_module
2138 implicit none
2140 ! Arguments
2141 integer, intent(in) :: start_mem_i, start_mem_j, start_mem_k, end_mem_i, end_mem_j, end_mem_k
2142 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j,start_mem_k:end_mem_k), intent(in) :: src_arr
2143 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j,start_mem_k:end_mem_k), intent(out) :: dst_arr
2144 integer, intent(in) :: sr_y
2145 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(in), optional :: mapfac
2147 ! Local variables
2148 integer :: i, j, k
2149 real :: dom_dx, dom_dy
2151 ! Get domain resolution
2152 call get_domain_resolution(dom_dx, dom_dy)
2154 if (present(mapfac)) then
2155 do k=start_mem_k,end_mem_k
2156 do i=start_mem_i, end_mem_i
2157 do j=start_mem_j+1, end_mem_j-1
2158 dst_arr(i,j,k) = (src_arr(i,j+1,k) - src_arr(i,j-1,k))/(2.*dom_dy*mapfac(i,j)/sr_y)
2159 end do
2160 end do
2162 do i=start_mem_i, end_mem_i
2163 dst_arr(i,start_mem_j,k) = (src_arr(i,start_mem_j+1,k) - src_arr(i,start_mem_j,k))/(dom_dy*mapfac(i,j)/sr_y)
2164 end do
2166 do i=start_mem_i, end_mem_i
2167 dst_arr(i,end_mem_j,k) = (src_arr(i,end_mem_j,k) - src_arr(i,end_mem_j-1,k))/(dom_dy*mapfac(i,j)/sr_y)
2168 end do
2169 end do
2170 else
2171 do k=start_mem_k,end_mem_k
2172 do i=start_mem_i, end_mem_i
2173 do j=start_mem_j+1, end_mem_j-1
2174 dst_arr(i,j,k) = (src_arr(i,j+1,k) - src_arr(i,j-1,k))/(2.*dom_dy/sr_y)
2175 end do
2176 end do
2178 do i=start_mem_i, end_mem_i
2179 dst_arr(i,start_mem_j,k) = (src_arr(i,start_mem_j+1,k) - src_arr(i,start_mem_j,k))/(dom_dy/sr_y)
2180 end do
2182 do i=start_mem_i, end_mem_i
2183 dst_arr(i,end_mem_j,k) = (src_arr(i,end_mem_j,k) - src_arr(i,end_mem_j-1,k))/(dom_dy/sr_y)
2184 end do
2185 end do
2186 end if
2188 end subroutine calc_dfdy
2191 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2192 ! Name: calc_dfdx
2193 !
2194 ! Purpose: This routine calculates df/dx for the field in src_arr, and places
2195 ! the result in dst_array.
2196 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2197 subroutine calc_dfdx(src_arr, dst_arr, start_mem_i, start_mem_j, &
2198 start_mem_k, end_mem_i, end_mem_j, end_mem_k, sr_x, mapfac)
2200 ! Modules
2201 use llxy_module
2203 implicit none
2205 ! Arguments
2206 integer, intent(in) :: start_mem_i, start_mem_j, start_mem_k, end_mem_i, end_mem_j, end_mem_k
2207 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j, start_mem_k:end_mem_k), intent(in) :: src_arr
2208 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j, start_mem_k:end_mem_k), intent(out) :: dst_arr
2209 integer, intent(in) :: sr_x
2210 real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(in), optional :: mapfac
2212 ! Local variables
2213 integer :: i, j, k
2214 real :: dom_dx, dom_dy
2216 ! Get domain resolution
2217 call get_domain_resolution(dom_dx, dom_dy)
2219 if (present(mapfac)) then
2220 do k=start_mem_k, end_mem_k
2221 do i=start_mem_i+1, end_mem_i-1
2222 do j=start_mem_j, end_mem_j
2223 dst_arr(i,j,k) = (src_arr(i+1,j,k) - src_arr(i-1,j,k))/(2.*dom_dx*mapfac(i,j)/sr_x)
2224 end do
2225 end do
2227 do j=start_mem_j, end_mem_j
2228 dst_arr(start_mem_i,j,k) = (src_arr(start_mem_i+1,j,k) - src_arr(start_mem_i,j,k))/(dom_dx*mapfac(i,j)/sr_x)
2229 end do
2231 do j=start_mem_j, end_mem_j
2232 dst_arr(end_mem_i,j,k) = (src_arr(end_mem_i,j,k) - src_arr(end_mem_i-1,j,k))/(dom_dx*mapfac(i,j)/sr_x)
2233 end do
2234 end do
2235 else
2236 do k=start_mem_k, end_mem_k
2237 do i=start_mem_i+1, end_mem_i-1
2238 do j=start_mem_j, end_mem_j
2239 dst_arr(i,j,k) = (src_arr(i+1,j,k) - src_arr(i-1,j,k))/(2.*dom_dx/sr_x)
2240 end do
2241 end do
2243 do j=start_mem_j, end_mem_j
2244 dst_arr(start_mem_i,j,k) = (src_arr(start_mem_i+1,j,k) - src_arr(start_mem_i,j,k))/(dom_dx/sr_x)
2245 end do
2247 do j=start_mem_j, end_mem_j
2248 dst_arr(end_mem_i,j,k) = (src_arr(end_mem_i,j,k) - src_arr(end_mem_i-1,j,k))/(dom_dx/sr_x)
2249 end do
2250 end do
2251 end if
2253 end subroutine calc_dfdx
2255 end module process_tile_module