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