| blob | 89be1a1236ecab5d4fbac7d671d355027693288e |
1 module process_domain_module
3 contains
5 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
6 ! Name: process_domain
7 !
8 ! Purpose:
9 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
10 subroutine process_domain(n, extra_row, extra_col)
12 use date_pack
13 use gridinfo_module
14 use interp_option_module
15 use misc_definitions_module
16 use module_debug
17 use storage_module
19 implicit none
21 ! Arguments
22 integer, intent(in) :: n
23 logical, intent(in) :: extra_row, extra_col
25 ! Local variables
26 integer :: i, t, dyn_opt, &
27 we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
28 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
29 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
30 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
31 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
32 idiff, n_times, &
33 west_east_dim, south_north_dim, bottom_top_dim, map_proj, &
34 is_water, is_lake, is_ice, is_urban, i_soilwater, &
35 grid_id, parent_id, i_parent_start, j_parent_start, &
36 i_parent_end, j_parent_end, parent_grid_ratio, sub_x, sub_y, num_land_cat, process_bdy_width
37 real :: cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
38 dom_dx, dom_dy, pole_lat, pole_lon
39 real, dimension(16) :: corner_lats, corner_lons
40 real, pointer, dimension(:,:) :: landmask
41 real, pointer, dimension(:,:) :: xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v
42 logical, allocatable, dimension(:) :: got_this_field, got_const_field
43 character (len=19) :: valid_date, temp_date
44 character (len=128) :: title, mminlu
45 character (len=128), allocatable, dimension(:) :: output_flags, td_output_flags
47 ! Compute number of times that we will process
48 call geth_idts(end_date(n), start_date(n), idiff)
49 call mprintf((idiff < 0),ERROR,'Ending date is earlier than starting date in namelist for domain %i.', i1=n)
51 n_times = idiff / interval_seconds
53 ! Check that the interval evenly divides the range of times to process
54 call mprintf((mod(idiff, interval_seconds) /= 0),WARN, &
55 'In namelist, interval_seconds does not evenly divide '// &
56 '(end_date - start_date) for domain %i. Only %i time periods '// &
57 'will be processed.', i1=n, i2=n_times)
59 ! Initialize the storage module
60 call mprintf(.true.,LOGFILE,'Initializing storage module')
61 call storage_init()
63 !
64 ! Do time-independent processing
65 !
66 call get_static_fields(n, dyn_opt, west_east_dim, south_north_dim, bottom_top_dim, map_proj, &
67 we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
68 we_patch_s, we_patch_e, &
69 we_patch_stag_s, we_patch_stag_e, &
70 sn_patch_s, sn_patch_e, &
71 sn_patch_stag_s, sn_patch_stag_e, &
72 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
73 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
74 mminlu, num_land_cat, is_water, is_lake, is_ice, is_urban, i_soilwater, &
75 grid_id, parent_id, &
76 i_parent_start, j_parent_start, i_parent_end, j_parent_end, &
77 parent_grid_ratio, sub_x, sub_y, &
78 cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
79 pole_lat, pole_lon, dom_dx, dom_dy, landmask, xlat, xlon, xlat_u, xlon_u, &
80 xlat_v, xlon_v, corner_lats, corner_lons, title)
83 allocate(output_flags(num_entries))
84 allocate(got_const_field(num_entries))
86 do i=1,num_entries
87 output_flags(i) = ' '
88 got_const_field(i) = .false.
89 end do
91 ! This call is to process the constant met fields (SST or SEAICE, for example)
92 ! That we process constant fields is indicated by the first argument
93 call process_single_met_time(.true., temp_date, n, extra_row, extra_col, xlat, xlon, &
94 xlat_u, xlon_u, xlat_v, xlon_v, landmask, &
95 title, dyn_opt, &
96 west_east_dim, south_north_dim, &
97 we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
98 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
99 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
100 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
101 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
102 got_const_field, &
103 map_proj, mminlu, num_land_cat, is_water, is_lake, is_ice, is_urban, i_soilwater, &
104 grid_id, parent_id, i_parent_start, &
105 j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
106 cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
107 pole_lat, pole_lon, parent_grid_ratio, sub_x, sub_y, &
108 corner_lats, corner_lons, output_flags, 0)
110 !
111 ! Begin time-dependent processing
112 !
114 allocate(td_output_flags(num_entries))
115 allocate(got_this_field (num_entries))
117 ! Loop over all times to be processed for this domain
118 do t=0,n_times
120 call geth_newdate(valid_date, trim(start_date(n)), t*interval_seconds)
121 temp_date = ' '
123 if (mod(interval_seconds,3600) == 0) then
124 write(temp_date,'(a13)') valid_date(1:10)//'_'//valid_date(12:13)
125 else if (mod(interval_seconds,60) == 0) then
126 write(temp_date,'(a16)') valid_date(1:10)//'_'//valid_date(12:16)
127 else
128 write(temp_date,'(a19)') valid_date(1:10)//'_'//valid_date(12:19)
129 end if
131 call mprintf(.true.,STDOUT, ' Processing %s', s1=trim(temp_date))
132 call mprintf(.true.,LOGFILE, 'Preparing to process output time %s', s1=temp_date)
134 do i=1,num_entries
135 td_output_flags(i) = output_flags(i)
136 got_this_field(i) = got_const_field(i)
137 end do
139 if (t > 0) then
140 process_bdy_width = process_only_bdy
141 else
142 process_bdy_width = 0
143 end if
145 call process_single_met_time(.false., temp_date, n, extra_row, extra_col, xlat, xlon, &
146 xlat_u, xlon_u, xlat_v, xlon_v, landmask, &
147 title, dyn_opt, &
148 west_east_dim, south_north_dim, &
149 we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
150 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
151 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
152 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
153 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
154 got_this_field, &
155 map_proj, mminlu, num_land_cat, is_water, is_lake, is_ice, is_urban, i_soilwater, &
156 grid_id, parent_id, i_parent_start, &
157 j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
158 cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
159 pole_lat, pole_lon, parent_grid_ratio, sub_x, sub_y, &
160 corner_lats, corner_lons, td_output_flags, process_bdy_width)
162 end do ! Loop over n_times
165 deallocate(td_output_flags)
166 deallocate(got_this_field)
168 deallocate(output_flags)
169 deallocate(got_const_field)
171 call storage_delete_all()
173 end subroutine process_domain
176 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
177 ! Name: get_static_fields
178 !
179 ! Purpose:
180 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
181 subroutine get_static_fields(n, dyn_opt, west_east_dim, south_north_dim, bottom_top_dim, &
182 map_proj, &
183 we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
184 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
185 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
186 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
187 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
188 mminlu, num_land_cat, is_water, is_lake, is_ice, is_urban, i_soilwater, &
189 grid_id, parent_id, &
190 i_parent_start, j_parent_start, i_parent_end, j_parent_end, &
191 parent_grid_ratio, sub_x, sub_y, &
192 cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
193 pole_lat, pole_lon, dom_dx, dom_dy, landmask, xlat, xlon, xlat_u, xlon_u, &
194 xlat_v, xlon_v, corner_lats, corner_lons, title)
196 use gridinfo_module
197 use input_module
198 use llxy_module
199 use parallel_module
200 use storage_module
201 use module_debug
203 implicit none
205 ! Arguments
206 integer, intent(in) :: n
207 integer, intent(inout) :: dyn_opt, west_east_dim, south_north_dim, bottom_top_dim, &
208 map_proj, &
209 we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
210 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
211 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
212 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
213 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
214 is_water, is_lake, is_ice, is_urban, i_soilwater, grid_id, parent_id, &
215 i_parent_start, j_parent_start, i_parent_end, j_parent_end, &
216 parent_grid_ratio, sub_x, sub_y, num_land_cat
217 real, pointer, dimension(:,:) :: landmask
218 real, intent(inout) :: cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
219 dom_dx, dom_dy, pole_lat, pole_lon
220 real, pointer, dimension(:,:) :: xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v
221 real, dimension(16), intent(out) :: corner_lats, corner_lons
222 character (len=128), intent(inout) :: title, mminlu
224 ! Local variables
225 integer :: istatus, i, j, k, sp1, ep1, sp2, ep2, sp3, ep3, &
226 lh_mult, rh_mult, bh_mult, th_mult, subx, suby
227 integer :: we_mem_subgrid_s, we_mem_subgrid_e, &
228 sn_mem_subgrid_s, sn_mem_subgrid_e
229 integer :: we_patch_subgrid_s, we_patch_subgrid_e, &
230 sn_patch_subgrid_s, sn_patch_subgrid_e
231 real, pointer, dimension(:,:,:) :: real_array
232 character (len=3) :: memorder
233 character (len=128) :: grid_type, datestr, cname, stagger, cunits, cdesc
234 character (len=128), dimension(3) :: dimnames
235 type (fg_input) :: field
237 ! Initialize the input module to read static input data for this domain
238 call mprintf(.true.,LOGFILE,'Opening static input file.')
239 call input_init(n, istatus)
240 call mprintf((istatus /= 0),ERROR, 'input_init(): Error opening input for domain %i.', i1=n)
242 ! Read global attributes from the static data input file
243 call mprintf(.true.,LOGFILE,'Reading static global attributes.')
244 call read_global_attrs(title, datestr, grid_type, dyn_opt, west_east_dim, &
245 south_north_dim, bottom_top_dim, &
246 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
247 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
248 map_proj, mminlu, num_land_cat, &
249 is_water, is_lake, is_ice, is_urban, i_soilwater, &
250 grid_id, parent_id, i_parent_start, &
251 j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
252 cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, &
253 truelat2, pole_lat, pole_lon, parent_grid_ratio, &
254 corner_lats, corner_lons, sub_x, sub_y)
256 we_dom_s = 1
257 sn_dom_s = 1
258 if (grid_type(1:1) == 'C') then
259 we_dom_e = west_east_dim - 1
260 sn_dom_e = south_north_dim - 1
261 else if (grid_type(1:1) == 'E') then
262 we_dom_e = west_east_dim
263 sn_dom_e = south_north_dim
264 end if
266 ! Given the full dimensions of this domain, find out the range of indices
267 ! that will be worked on by this processor. This information is given by
268 ! my_minx, my_miny, my_maxx, my_maxy
269 call parallel_get_tile_dims(west_east_dim, south_north_dim)
271 ! Must figure out patch dimensions from info in parallel module
272 if (nprocs > 1 .and. .not. do_tiled_input) then
274 we_patch_s = my_minx
275 we_patch_stag_s = my_minx
276 we_patch_e = my_maxx - 1
277 sn_patch_s = my_miny
278 sn_patch_stag_s = my_miny
279 sn_patch_e = my_maxy - 1
281 if (gridtype == 'C') then
282 if (my_x /= nproc_x - 1) then
283 we_patch_e = we_patch_e + 1
284 we_patch_stag_e = we_patch_e
285 else
286 we_patch_stag_e = we_patch_e + 1
287 end if
288 if (my_y /= nproc_y - 1) then
289 sn_patch_e = sn_patch_e + 1
290 sn_patch_stag_e = sn_patch_e
291 else
292 sn_patch_stag_e = sn_patch_e + 1
293 end if
294 else if (gridtype == 'E') then
295 we_patch_e = we_patch_e + 1
296 sn_patch_e = sn_patch_e + 1
297 we_patch_stag_e = we_patch_e
298 sn_patch_stag_e = sn_patch_e
299 end if
301 end if
303 ! Compute multipliers for halo width; these must be 0/1
304 if (my_x /= 0) then
305 lh_mult = 1
306 else
307 lh_mult = 0
308 end if
309 if (my_x /= (nproc_x-1)) then
310 rh_mult = 1
311 else
312 rh_mult = 0
313 end if
314 if (my_y /= 0) then
315 bh_mult = 1
316 else
317 bh_mult = 0
318 end if
319 if (my_y /= (nproc_y-1)) then
320 th_mult = 1
321 else
322 th_mult = 0
323 end if
325 we_mem_s = we_patch_s - HALO_WIDTH*lh_mult
326 we_mem_e = we_patch_e + HALO_WIDTH*rh_mult
327 sn_mem_s = sn_patch_s - HALO_WIDTH*bh_mult
328 sn_mem_e = sn_patch_e + HALO_WIDTH*th_mult
329 we_mem_stag_s = we_patch_stag_s - HALO_WIDTH*lh_mult
330 we_mem_stag_e = we_patch_stag_e + HALO_WIDTH*rh_mult
331 sn_mem_stag_s = sn_patch_stag_s - HALO_WIDTH*bh_mult
332 sn_mem_stag_e = sn_patch_stag_e + HALO_WIDTH*th_mult
334 ! Initialize a proj_info type for the destination grid projection. This will
335 ! primarily be used for rotating Earth-relative winds to grid-relative winds
336 call set_domain_projection(map_proj, stand_lon, truelat1, truelat2, &
337 dom_dx, dom_dy, dom_dx, dom_dy, west_east_dim, &
338 south_north_dim, real(west_east_dim)/2., &
339 real(south_north_dim)/2.,cen_lat, cen_lon, &
340 cen_lat, cen_lon)
342 ! Read static fields using the input module; we know that there are no more
343 ! fields to be read when read_next_field() returns a non-zero status.
344 istatus = 0
345 do while (istatus == 0)
346 call read_next_field(sp1, ep1, sp2, ep2, sp3, ep3, cname, cunits, cdesc, &
347 memorder, stagger, dimnames, subx, suby, &
348 real_array, istatus)
349 if (istatus == 0) then
351 call mprintf(.true.,LOGFILE, 'Read in static field %s.',s1=cname)
353 ! We will also keep copies in core of the lat/lon arrays, for use in
354 ! interpolation of the met fields to the model grid.
355 ! For now, we assume that the lat/lon arrays will have known field names
356 if (index(cname, 'XLAT_M') /= 0 .and. &
357 len_trim(cname) == len_trim('XLAT_M')) then
358 allocate(xlat(we_mem_s:we_mem_e,sn_mem_s:sn_mem_e))
359 xlat(we_patch_s:we_patch_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
360 call exchange_halo_r(xlat, &
361 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, 1, 1, &
362 we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, 1, 1)
364 else if (index(cname, 'XLONG_M') /= 0 .and. &
365 len_trim(cname) == len_trim('XLONG_M')) then
366 allocate(xlon(we_mem_s:we_mem_e,sn_mem_s:sn_mem_e))
367 xlon(we_patch_s:we_patch_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
368 call exchange_halo_r(xlon, &
369 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, 1, 1, &
370 we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, 1, 1)
372 else if (index(cname, 'XLAT_U') /= 0 .and. &
373 len_trim(cname) == len_trim('XLAT_U')) then
374 allocate(xlat_u(we_mem_stag_s:we_mem_stag_e,sn_mem_s:sn_mem_e))
375 xlat_u(we_patch_stag_s:we_patch_stag_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
376 call exchange_halo_r(xlat_u, &
377 we_mem_stag_s, we_mem_stag_e, sn_mem_s, sn_mem_e, 1, 1, &
378 we_patch_stag_s, we_patch_stag_e, sn_patch_s, sn_patch_e, 1, 1)
380 else if (index(cname, 'XLONG_U') /= 0 .and. &
381 len_trim(cname) == len_trim('XLONG_U')) then
382 allocate(xlon_u(we_mem_stag_s:we_mem_stag_e,sn_mem_s:sn_mem_e))
383 xlon_u(we_patch_stag_s:we_patch_stag_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
384 call exchange_halo_r(xlon_u, &
385 we_mem_stag_s, we_mem_stag_e, sn_mem_s, sn_mem_e, 1, 1, &
386 we_patch_stag_s, we_patch_stag_e, sn_patch_s, sn_patch_e, 1, 1)
388 else if (index(cname, 'XLAT_V') /= 0 .and. &
389 len_trim(cname) == len_trim('XLAT_V')) then
390 allocate(xlat_v(we_mem_s:we_mem_e,sn_mem_stag_s:sn_mem_stag_e))
391 xlat_v(we_patch_s:we_patch_e,sn_patch_stag_s:sn_patch_stag_e) = real_array(:,:,1)
392 call exchange_halo_r(xlat_v, &
393 we_mem_s, we_mem_e, sn_mem_stag_s, sn_mem_stag_e, 1, 1, &
394 we_patch_s, we_patch_e, sn_patch_stag_s, sn_patch_stag_e, 1, 1)
396 else if (index(cname, 'XLONG_V') /= 0 .and. &
397 len_trim(cname) == len_trim('XLONG_V')) then
398 allocate(xlon_v(we_mem_s:we_mem_e,sn_mem_stag_s:sn_mem_stag_e))
399 xlon_v(we_patch_s:we_patch_e,sn_patch_stag_s:sn_patch_stag_e) = real_array(:,:,1)
400 call exchange_halo_r(xlon_v, &
401 we_mem_s, we_mem_e, sn_mem_stag_s, sn_mem_stag_e, 1, 1, &
402 we_patch_s, we_patch_e, sn_patch_stag_s, sn_patch_stag_e, 1, 1)
404 else if (index(cname, 'LANDMASK') /= 0 .and. &
405 len_trim(cname) == len_trim('LANDMASK')) then
406 allocate(landmask(we_mem_s:we_mem_e,sn_mem_s:sn_mem_e))
407 landmask(we_patch_s:we_patch_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
408 call exchange_halo_r(landmask, &
409 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, 1, 1, &
410 we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, 1, 1)
412 end if
414 if (subx > 1) then
415 we_mem_subgrid_s = (we_mem_s + HALO_WIDTH*lh_mult - 1)*subx - HALO_WIDTH*lh_mult + 1
416 we_mem_subgrid_e = (we_mem_e + (1-rh_mult) - HALO_WIDTH*rh_mult )*subx + HALO_WIDTH*rh_mult
417 we_patch_subgrid_s = (we_patch_s - 1)*subx + 1
418 we_patch_subgrid_e = (we_patch_e + (1-rh_mult) )*subx
419 end if
420 if (suby > 1) then
421 sn_mem_subgrid_s = (sn_mem_s + HALO_WIDTH*bh_mult - 1)*suby - HALO_WIDTH*bh_mult + 1
422 sn_mem_subgrid_e = (sn_mem_e + (1-th_mult) - HALO_WIDTH*th_mult )*suby + HALO_WIDTH*th_mult
423 sn_patch_subgrid_s = (sn_patch_s - 1)*suby + 1
424 sn_patch_subgrid_e = (sn_patch_e + (1-th_mult) )*suby
425 end if
427 ! Having read in a field, we write each level individually to the
428 ! storage module; levels will be reassembled later on when they
429 ! are written.
430 do k=sp3,ep3
431 field%header%sr_x=subx
432 field%header%sr_y=suby
433 field%header%version = 1
434 field%header%date = start_date(n)
435 field%header%time_dependent = .false.
436 field%header%mask_field = .false.
437 field%header%forecast_hour = 0.0
438 field%header%fg_source = 'geogrid_model'
439 field%header%field = cname
440 field%header%units = cunits
441 field%header%description = cdesc
442 field%header%vertical_coord = dimnames(3)
443 field%header%vertical_level = k
444 field%header%array_order = memorder
445 field%header%is_wind_grid_rel = .true.
446 field%header%array_has_missing_values = .false.
447 if (gridtype == 'C') then
448 if (subx > 1 .or. suby > 1) then
449 field%map%stagger = M
450 field%header%dim1(1) = we_mem_subgrid_s
451 field%header%dim1(2) = we_mem_subgrid_e
452 field%header%dim2(1) = sn_mem_subgrid_s
453 field%header%dim2(2) = sn_mem_subgrid_e
454 else if (trim(stagger) == 'M') then
455 field%map%stagger = M
456 field%header%dim1(1) = we_mem_s
457 field%header%dim1(2) = we_mem_e
458 field%header%dim2(1) = sn_mem_s
459 field%header%dim2(2) = sn_mem_e
460 else if (trim(stagger) == 'U') then
461 field%map%stagger = U
462 field%header%dim1(1) = we_mem_stag_s
463 field%header%dim1(2) = we_mem_stag_e
464 field%header%dim2(1) = sn_mem_s
465 field%header%dim2(2) = sn_mem_e
466 else if (trim(stagger) == 'V') then
467 field%map%stagger = V
468 field%header%dim1(1) = we_mem_s
469 field%header%dim1(2) = we_mem_e
470 field%header%dim2(1) = sn_mem_stag_s
471 field%header%dim2(2) = sn_mem_stag_e
472 end if
473 else if (gridtype == 'E') then
474 if (trim(stagger) == 'M') then
475 field%map%stagger = HH
476 else if (trim(stagger) == 'V') then
477 field%map%stagger = VV
478 end if
479 field%header%dim1(1) = we_mem_s
480 field%header%dim1(2) = we_mem_e
481 field%header%dim2(1) = sn_mem_s
482 field%header%dim2(2) = sn_mem_e
483 end if
485 allocate(field%valid_mask)
487 if (subx > 1 .or. suby > 1) then
488 allocate(field%r_arr(we_mem_subgrid_s:we_mem_subgrid_e,&
489 sn_mem_subgrid_s:sn_mem_subgrid_e))
490 field%r_arr(we_patch_subgrid_s:we_patch_subgrid_e,sn_patch_subgrid_s:sn_patch_subgrid_e) = &
491 real_array(sp1:ep1,sp2:ep2,k)
492 call exchange_halo_r(field%r_arr, &
493 we_mem_subgrid_s, we_mem_subgrid_e, sn_mem_subgrid_s, sn_mem_subgrid_e, 1, 1, &
494 we_patch_subgrid_s, we_patch_subgrid_e, sn_patch_subgrid_s, sn_patch_subgrid_e, 1, 1)
495 call bitarray_create(field%valid_mask, &
496 (we_mem_subgrid_e-we_mem_subgrid_s)+1, &
497 (sn_mem_subgrid_e-sn_mem_subgrid_s)+1)
498 do j=1,(sn_mem_subgrid_e-sn_mem_subgrid_s)+1
499 do i=1,(we_mem_subgrid_e-we_mem_subgrid_s)+1
500 call bitarray_set(field%valid_mask, i, j)
501 end do
502 end do
504 else if (field%map%stagger == M .or. &
505 field%map%stagger == HH .or. &
506 field%map%stagger == VV) then
507 allocate(field%r_arr(we_mem_s:we_mem_e,&
508 sn_mem_s:sn_mem_e))
509 field%r_arr(we_patch_s:we_patch_e,sn_patch_s:sn_patch_e) = real_array(sp1:ep1,sp2:ep2,k)
510 call exchange_halo_r(field%r_arr, &
511 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, 1, 1, &
512 we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, 1, 1)
513 call bitarray_create(field%valid_mask, &
514 (we_mem_e-we_mem_s)+1, &
515 (sn_mem_e-sn_mem_s)+1)
516 do j=1,(sn_mem_e-sn_mem_s)+1
517 do i=1,(we_mem_e-we_mem_s)+1
518 call bitarray_set(field%valid_mask, i, j)
519 end do
520 end do
521 else if (field%map%stagger == U) then
522 allocate(field%r_arr(we_mem_stag_s:we_mem_stag_e,&
523 sn_mem_s:sn_mem_e))
524 field%r_arr(we_patch_stag_s:we_patch_stag_e,sn_patch_s:sn_patch_e) = real_array(sp1:ep1,sp2:ep2,k)
525 call exchange_halo_r(field%r_arr, &
526 we_mem_stag_s, we_mem_stag_e, sn_mem_s, sn_mem_e, 1, 1, &
527 we_patch_stag_s, we_patch_stag_e, sn_patch_s, sn_patch_e, 1, 1)
528 call bitarray_create(field%valid_mask, &
529 (we_mem_stag_e-we_mem_stag_s)+1, &
530 (sn_mem_e-sn_mem_s)+1)
531 do j=1,(sn_mem_e-sn_mem_s)+1
532 do i=1,(we_mem_stag_e-we_mem_stag_s)+1
533 call bitarray_set(field%valid_mask, i, j)
534 end do
535 end do
536 else if (field%map%stagger == V) then
537 allocate(field%r_arr(we_mem_s:we_mem_e,&
538 sn_mem_stag_s:sn_mem_stag_e))
539 field%r_arr(we_patch_s:we_patch_e,sn_patch_stag_s:sn_patch_stag_e) = real_array(sp1:ep1,sp2:ep2,k)
540 call exchange_halo_r(field%r_arr, &
541 we_mem_s, we_mem_e, sn_mem_stag_s, sn_mem_stag_e, 1, 1, &
542 we_patch_s, we_patch_e, sn_patch_stag_s, sn_patch_stag_e, 1, 1)
543 call bitarray_create(field%valid_mask, &
544 (we_mem_e-we_mem_s)+1, &
545 (sn_mem_stag_e-sn_mem_stag_s)+1)
546 do j=1,(sn_mem_stag_e-sn_mem_stag_s)+1
547 do i=1,(we_mem_e-we_mem_s)+1
548 call bitarray_set(field%valid_mask, i, j)
549 end do
550 end do
551 end if
553 nullify(field%modified_mask)
555 call storage_put_field(field)
557 end do
559 end if
560 end do
562 ! Done reading all static fields for this domain
563 call input_close()
565 end subroutine get_static_fields
568 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
569 ! Name: process_single_met_time
570 !
571 ! Purpose:
572 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
573 subroutine process_single_met_time(do_const_processing, &
574 temp_date, n, extra_row, extra_col, xlat, xlon, &
575 xlat_u, xlon_u, xlat_v, xlon_v, landmask, &
576 title, dyn_opt, &
577 west_east_dim, south_north_dim, &
578 we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
579 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
580 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
581 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
582 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
583 got_this_field, &
584 map_proj, mminlu, num_land_cat, is_water, is_lake, is_ice, is_urban, i_soilwater, &
585 grid_id, parent_id, i_parent_start, &
586 j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
587 cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
588 pole_lat, pole_lon, parent_grid_ratio, sub_x, sub_y, &
589 corner_lats, corner_lons, output_flags, process_bdy_width)
591 use bitarray_module
592 use gridinfo_module
593 use interp_module
594 use interp_option_module
595 use llxy_module
596 use misc_definitions_module
597 use module_debug
598 use output_module
599 use parallel_module
600 use read_met_module
601 use rotate_winds_module
602 use storage_module
604 implicit none
606 ! Arguments
607 logical, intent(in) :: do_const_processing
608 integer, intent(in) :: n, dyn_opt, west_east_dim, south_north_dim, map_proj, &
609 we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
610 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
611 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
612 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
613 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
614 is_water, is_lake, is_ice, is_urban, i_soilwater, &
615 grid_id, parent_id, i_parent_start, j_parent_start, &
616 i_parent_end, j_parent_end, parent_grid_ratio, sub_x, sub_y, num_land_cat, &
617 process_bdy_width
618 ! BUG: Should we be passing these around as pointers, or just declare them as arrays?
619 real, pointer, dimension(:,:) :: landmask
620 real, intent(in) :: dom_dx, dom_dy, cen_lat, moad_cen_lat, cen_lon, stand_lon, &
621 truelat1, truelat2, pole_lat, pole_lon
622 real, dimension(16), intent(in) :: corner_lats, corner_lons
623 real, pointer, dimension(:,:) :: xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v
624 logical, intent(in) :: extra_row, extra_col
625 logical, dimension(:), intent(inout) :: got_this_field
626 character (len=19), intent(in) :: temp_date
627 character (len=128), intent(in) :: mminlu
628 character (len=128), dimension(:), intent(inout) :: output_flags
630 ! BUG: Move this constant to misc_definitions_module?
631 integer, parameter :: BDR_WIDTH = 3
633 ! Local variables
634 integer :: istatus, iqstatus, fg_idx, idx, idxt, i, j, bottom_top_dim, &
635 sm1, em1, sm2, em2, sm3, em3, &
636 sp1, ep1, sp2, ep2, sp3, ep3, &
637 sd1, ed1, sd2, ed2, sd3, ed3, &
638 u_idx, bdr_wdth
639 integer :: nmet_flags
640 integer :: num_metgrid_soil_levs
641 integer, pointer, dimension(:) :: soil_levels
642 real :: rx, ry
643 real :: threshold
644 logical :: do_gcell_interp
645 integer, pointer, dimension(:) :: u_levels, v_levels
646 real, pointer, dimension(:,:) :: halo_slab
647 real, pointer, dimension(:,:,:) :: real_array
648 character (len=19) :: output_date
649 character (len=128) :: cname, title
650 character (len=MAX_FILENAME_LEN) :: input_name
651 character (len=128), allocatable, dimension(:) :: met_flags
652 type (fg_input) :: field, u_field, v_field
653 type (met_data) :: fg_data
656 ! For this time, we need to process all first-guess filename roots. When we
657 ! hit a root containing a '*', we assume we have hit the end of the list
658 fg_idx = 1
659 if (do_const_processing) then
660 input_name = constants_name(fg_idx)
661 else
662 input_name = fg_name(fg_idx)
663 end if
664 do while (input_name /= '*')
666 call mprintf(.true.,STDOUT, ' %s', s1=input_name)
667 call mprintf(.true.,LOGFILE, 'Getting input fields from %s', s1=input_name)
669 ! Do a first pass through this fg source to get all mask fields used
670 ! during interpolation
671 call get_interp_masks(trim(input_name), do_const_processing, temp_date)
673 istatus = 0
675 ! Initialize the module for reading in the met fields
676 call read_met_init(trim(input_name), do_const_processing, temp_date, istatus)
678 if (istatus == 0) then
680 ! Process all fields and levels from the current file; read_next_met_field()
681 ! will return a non-zero status when there are no more fields to be read.
682 do while (istatus == 0)
684 call read_next_met_field(fg_data, istatus)
686 if (istatus == 0) then
688 ! Find index into fieldname, interp_method, masked, and fill_missing
689 ! of the current field
690 idxt = num_entries + 1
691 do idx=1,num_entries
692 if ((index(fieldname(idx), trim(fg_data%field)) /= 0) .and. &
693 (len_trim(fieldname(idx)) == len_trim(fg_data%field))) then
695 got_this_field(idx) = .true.
697 if (index(input_name,trim(from_input(idx))) /= 0 .or. &
698 (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
699 idxt = idx
700 end if
702 end if
703 end do
704 idx = idxt
705 if (idx > num_entries) idx = num_entries ! The last entry is a default
707 ! Do we need to rename this field?
708 if (output_name(idx) /= ' ') then
709 fg_data%field = output_name(idx)(1:9)
711 idxt = num_entries + 1
712 do idx=1,num_entries
713 if ((index(fieldname(idx), trim(fg_data%field)) /= 0) .and. &
714 (len_trim(fieldname(idx)) == len_trim(fg_data%field))) then
716 got_this_field(idx) = .true.
718 if (index(input_name,trim(from_input(idx))) /= 0 .or. &
719 (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
720 idxt = idx
721 end if
723 end if
724 end do
725 idx = idxt
726 if (idx > num_entries) idx = num_entries ! The last entry is a default
727 end if
729 ! Do a simple check to see whether this is a global dataset
730 ! Note that we do not currently support regional Gaussian grids
731 if ((fg_data%iproj == PROJ_LATLON .and. abs(fg_data%nx * fg_data%deltalon - 360.) < 0.0001) &
732 .or. (fg_data%iproj == PROJ_GAUSS)) then
733 bdr_wdth = BDR_WIDTH
734 allocate(halo_slab(1-BDR_WIDTH:fg_data%nx+BDR_WIDTH,1:fg_data%ny))
736 halo_slab(1:fg_data%nx, 1:fg_data%ny) = &
737 fg_data%slab(1:fg_data%nx, 1:fg_data%ny)
739 halo_slab(1-BDR_WIDTH:0, 1:fg_data%ny) = &
740 fg_data%slab(fg_data%nx-BDR_WIDTH+1:fg_data%nx, 1:fg_data%ny)
742 halo_slab(fg_data%nx+1:fg_data%nx+BDR_WIDTH, 1:fg_data%ny) = &
743 fg_data%slab(1:BDR_WIDTH, 1:fg_data%ny)
745 deallocate(fg_data%slab)
746 else
747 bdr_wdth = 0
748 halo_slab => fg_data%slab
749 nullify(fg_data%slab)
750 end if
752 call mprintf(.true.,LOGFILE,'Processing %s at level %f.',s1=fg_data%field,f1=fg_data%xlvl)
754 call push_source_projection(fg_data%iproj, fg_data%xlonc, fg_data%truelat1, &
755 fg_data%truelat2, fg_data%dx, fg_data%dy, fg_data%deltalat, &
756 fg_data%deltalon, fg_data%starti, fg_data%startj, &
757 fg_data%startlat, fg_data%startlon, earth_radius=fg_data%earth_radius*1000.)
759 ! Initialize fg_input structure to store the field
760 field%header%version = 1
761 field%header%date = fg_data%hdate//' '
762 if (do_const_processing) then
763 field%header%time_dependent = .false.
764 else
765 field%header%time_dependent = .true.
766 end if
767 field%header%forecast_hour = fg_data%xfcst
768 field%header%fg_source = 'FG'
769 field%header%field = ' '
770 field%header%field(1:9) = fg_data%field
771 field%header%units = ' '
772 field%header%units(1:25) = fg_data%units
773 field%header%description = ' '
774 field%header%description(1:46) = fg_data%desc
775 call get_z_dim_name(fg_data%field,field%header%vertical_coord)
776 field%header%vertical_level = nint(fg_data%xlvl)
777 field%header%sr_x = 1
778 field%header%sr_y = 1
779 field%header%array_order = 'XY '
780 field%header%is_wind_grid_rel = fg_data%is_wind_grid_rel
781 field%header%array_has_missing_values = .false.
782 nullify(field%r_arr)
783 nullify(field%valid_mask)
784 nullify(field%modified_mask)
786 if (output_this_field(idx) .and. flag_in_output(idx) /= ' ') then
787 output_flags(idx) = flag_in_output(idx)
788 end if
790 ! If we should not output this field, just list it as a mask field
791 if (output_this_field(idx)) then
792 field%header%mask_field = .false.
793 else
794 field%header%mask_field = .true.
795 end if
797 !
798 ! Before actually doing any interpolation to the model grid, we must check
799 ! whether we will be using the average_gcell interpolator that averages all
800 ! source points in each model grid cell
801 !
802 do_gcell_interp = .false.
803 if (index(interp_method(idx),'average_gcell') /= 0) then
805 call get_gcell_threshold(interp_method(idx), threshold, istatus)
806 if (istatus == 0) then
807 if (fg_data%dx == 0. .and. fg_data%dy == 0. .and. &
808 fg_data%deltalat /= 0. .and. fg_data%deltalon /= 0.) then
809 fg_data%dx = abs(fg_data%deltalon)
810 fg_data%dy = abs(fg_data%deltalat)
811 else
812 ! BUG: Need to more correctly handle dx/dy in meters.
813 fg_data%dx = fg_data%dx / 111000. ! Convert meters to approximate degrees
814 fg_data%dy = fg_data%dy / 111000.
815 end if
816 if (gridtype == 'C') then
817 if (threshold*max(fg_data%dx,fg_data%dy)*111. <= max(dom_dx,dom_dy)/1000.) &
818 do_gcell_interp = .true.
819 else if (gridtype == 'E') then
820 if (threshold*max(fg_data%dx,fg_data%dy) <= max(dom_dx,dom_dy)) &
821 do_gcell_interp = .true.
822 end if
823 end if
824 end if
826 ! Interpolate to U staggering
827 if (output_stagger(idx) == U) then
829 call storage_query_field(field, iqstatus)
830 if (iqstatus == 0) then
831 call storage_get_field(field, iqstatus)
832 call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// &
833 ' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl))
834 if (associated(field%modified_mask)) then
835 call bitarray_destroy(field%modified_mask)
836 nullify(field%modified_mask)
837 end if
838 else
839 allocate(field%valid_mask)
840 call bitarray_create(field%valid_mask, we_mem_stag_e-we_mem_stag_s+1, sn_mem_e-sn_mem_s+1)
841 end if
843 ! Save a copy of the fg_input structure for the U field so that we can find it later
844 if (is_u_field(idx)) call dup(field, u_field)
846 allocate(field%modified_mask)
847 call bitarray_create(field%modified_mask, we_mem_stag_e-we_mem_stag_s+1, sn_mem_e-sn_mem_s+1)
849 if (do_const_processing .or. field%header%time_dependent) then
850 call interp_met_field(input_name, fg_data%field, U, M, &
851 field, xlat_u, xlon_u, we_mem_stag_s, we_mem_stag_e, sn_mem_s, sn_mem_e, &
852 we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
853 halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
854 field%modified_mask, process_bdy_width)
855 else
856 call mprintf(.true.,INFORM,' - already processed this field from constant file.')
857 end if
859 ! Interpolate to V staggering
860 else if (output_stagger(idx) == V) then
862 call storage_query_field(field, iqstatus)
863 if (iqstatus == 0) then
864 call storage_get_field(field, iqstatus)
865 call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// &
866 ' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl))
867 if (associated(field%modified_mask)) then
868 call bitarray_destroy(field%modified_mask)
869 nullify(field%modified_mask)
870 end if
871 else
872 allocate(field%valid_mask)
873 call bitarray_create(field%valid_mask, we_mem_e-we_mem_s+1, sn_mem_stag_e-sn_mem_stag_s+1)
874 end if
876 ! Save a copy of the fg_input structure for the V field so that we can find it later
877 if (is_v_field(idx)) call dup(field, v_field)
879 allocate(field%modified_mask)
880 call bitarray_create(field%modified_mask, we_mem_e-we_mem_s+1, sn_mem_stag_e-sn_mem_stag_s+1)
882 if (do_const_processing .or. field%header%time_dependent) then
883 call interp_met_field(input_name, fg_data%field, V, M, &
884 field, xlat_v, xlon_v, we_mem_s, we_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
885 we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
886 halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
887 field%modified_mask, process_bdy_width)
888 else
889 call mprintf(.true.,INFORM,' - already processed this field from constant file.')
890 end if
892 ! Interpolate to VV staggering
893 else if (output_stagger(idx) == VV) then
895 call storage_query_field(field, iqstatus)
896 if (iqstatus == 0) then
897 call storage_get_field(field, iqstatus)
898 call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// &
899 ' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl))
900 if (associated(field%modified_mask)) then
901 call bitarray_destroy(field%modified_mask)
902 nullify(field%modified_mask)
903 end if
904 else
905 allocate(field%valid_mask)
906 call bitarray_create(field%valid_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1)
907 end if
909 ! Save a copy of the fg_input structure for the U field so that we can find it later
910 if (is_u_field(idx)) call dup(field, u_field)
912 ! Save a copy of the fg_input structure for the V field so that we can find it later
913 if (is_v_field(idx)) call dup(field, v_field)
915 allocate(field%modified_mask)
916 call bitarray_create(field%modified_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1)
918 if (do_const_processing .or. field%header%time_dependent) then
919 call interp_met_field(input_name, fg_data%field, VV, M, &
920 field, xlat_v, xlon_v, we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
921 we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
922 halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
923 field%modified_mask, process_bdy_width)
924 else
925 call mprintf(.true.,INFORM,' - already processed this field from constant file.')
926 end if
928 ! All other fields interpolated to M staggering for C grid, H staggering for E grid
929 else
931 call storage_query_field(field, iqstatus)
932 if (iqstatus == 0) then
933 call storage_get_field(field, iqstatus)
934 call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// &
935 ' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl))
936 if (associated(field%modified_mask)) then
937 call bitarray_destroy(field%modified_mask)
938 nullify(field%modified_mask)
939 end if
940 else
941 allocate(field%valid_mask)
942 call bitarray_create(field%valid_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1)
943 end if
945 allocate(field%modified_mask)
946 call bitarray_create(field%modified_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1)
948 if (do_const_processing .or. field%header%time_dependent) then
949 if (gridtype == 'C') then
950 call interp_met_field(input_name, fg_data%field, M, M, &
951 field, xlat, xlon, we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
952 we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
953 halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
954 field%modified_mask, process_bdy_width, landmask)
956 else if (gridtype == 'E') then
957 call interp_met_field(input_name, fg_data%field, HH, M, &
958 field, xlat, xlon, we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
959 we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
960 halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
961 field%modified_mask, process_bdy_width, landmask)
962 end if
963 else
964 call mprintf(.true.,INFORM,' - already processed this field from constant file.')
965 end if
967 end if
969 call bitarray_merge(field%valid_mask, field%modified_mask)
971 deallocate(halo_slab)
973 ! Store the interpolated field
974 call storage_put_field(field)
976 call pop_source_projection()
978 end if
979 end do
981 call read_met_close()
983 call push_source_projection(fg_data%iproj, fg_data%xlonc, fg_data%truelat1, &
984 fg_data%truelat2, fg_data%dx, fg_data%dy, fg_data%deltalat, &
985 fg_data%deltalon, fg_data%starti, fg_data%startj, &
986 fg_data%startlat, fg_data%startlon, earth_radius=fg_data%earth_radius*1000.)
988 !
989 ! If necessary, rotate winds to earth-relative for this fg source
990 !
992 call storage_get_levels(u_field, u_levels)
993 call storage_get_levels(v_field, v_levels)
995 if (associated(u_levels) .and. associated(v_levels)) then
996 u_idx = 1
997 do u_idx = 1, size(u_levels)
998 u_field%header%vertical_level = u_levels(u_idx)
999 call storage_get_field(u_field, istatus)
1000 v_field%header%vertical_level = v_levels(u_idx)
1001 call storage_get_field(v_field, istatus)
1003 if (associated(u_field%modified_mask) .and. &
1004 associated(v_field%modified_mask)) then
1006 if (u_field%header%is_wind_grid_rel) then
1007 if (gridtype == 'C') then
1008 call map_to_met(u_field%r_arr, u_field%modified_mask, &
1009 v_field%r_arr, v_field%modified_mask, &
1010 we_mem_stag_s, sn_mem_s, &
1011 we_mem_stag_e, sn_mem_e, &
1012 we_mem_s, sn_mem_stag_s, &
1013 we_mem_e, sn_mem_stag_e, &
1014 xlon_u, xlon_v, xlat_u, xlat_v)
1015 else if (gridtype == 'E') then
1016 call map_to_met_nmm(u_field%r_arr, u_field%modified_mask, &
1017 v_field%r_arr, v_field%modified_mask, &
1018 we_mem_s, sn_mem_s, &
1019 we_mem_e, sn_mem_e, &
1020 xlat_v, xlon_v)
1021 end if
1022 end if
1024 call bitarray_destroy(u_field%modified_mask)
1025 call bitarray_destroy(v_field%modified_mask)
1026 nullify(u_field%modified_mask)
1027 nullify(v_field%modified_mask)
1028 call storage_put_field(u_field)
1029 call storage_put_field(v_field)
1030 end if
1032 end do
1034 deallocate(u_levels)
1035 deallocate(v_levels)
1037 end if
1039 call pop_source_projection()
1041 else
1042 if (do_const_processing) then
1043 call mprintf(.true.,WARN,'Couldn''t open file %s for input.',s1=input_name)
1044 else
1045 call mprintf(.true.,WARN,'Couldn''t open file %s for input.',s1=trim(input_name)//':'//trim(temp_date))
1046 end if
1047 end if
1049 fg_idx = fg_idx + 1
1050 if (do_const_processing) then
1051 input_name = constants_name(fg_idx)
1052 else
1053 input_name = fg_name(fg_idx)
1054 end if
1055 end do ! while (input_name /= '*')
1057 !
1058 ! Rotate winds from earth-relative to grid-relative
1059 !
1061 call storage_get_levels(u_field, u_levels)
1062 call storage_get_levels(v_field, v_levels)
1064 if (associated(u_levels) .and. associated(v_levels)) then
1065 u_idx = 1
1066 do u_idx = 1, size(u_levels)
1067 u_field%header%vertical_level = u_levels(u_idx)
1068 call storage_get_field(u_field, istatus)
1069 v_field%header%vertical_level = v_levels(u_idx)
1070 call storage_get_field(v_field, istatus)
1072 if (gridtype == 'C') then
1073 call met_to_map(u_field%r_arr, u_field%valid_mask, &
1074 v_field%r_arr, v_field%valid_mask, &
1075 we_mem_stag_s, sn_mem_s, &
1076 we_mem_stag_e, sn_mem_e, &
1077 we_mem_s, sn_mem_stag_s, &
1078 we_mem_e, sn_mem_stag_e, &
1079 xlon_u, xlon_v, xlat_u, xlat_v)
1080 else if (gridtype == 'E') then
1081 call met_to_map_nmm(u_field%r_arr, u_field%valid_mask, &
1082 v_field%r_arr, v_field%valid_mask, &
1083 we_mem_s, sn_mem_s, &
1084 we_mem_e, sn_mem_e, &
1085 xlat_v, xlon_v)
1086 end if
1088 end do
1090 deallocate(u_levels)
1091 deallocate(v_levels)
1093 end if
1095 if (do_const_processing) return
1097 !
1098 ! Now that we have all degribbed fields, we build a 3-d pressure field, and fill in any
1099 ! missing levels in the other 3-d fields
1100 !
1101 call mprintf(.true.,LOGFILE,'Filling missing levels.')
1102 call fill_missing_levels(output_flags)
1104 call mprintf(.true.,LOGFILE,'Creating derived fields.')
1105 call create_derived_fields(gridtype, fg_data%hdate, fg_data%xfcst, &
1106 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
1107 we_mem_stag_s, we_mem_stag_e, sn_mem_stag_s, sn_mem_stag_e, &
1108 got_this_field, output_flags)
1110 !
1111 ! Check that every mandatory field was found in input data
1112 !
1113 do i=1,num_entries
1114 if (is_mandatory(i) .and. .not. got_this_field(i)) then
1115 call mprintf(.true.,ERROR,'The mandatory field %s was not found in any input data.',s1=fieldname(i))
1116 end if
1117 end do
1119 !
1120 ! Before we begin to write fields, if debug_level is set high enough, we
1121 ! write a table of which fields are available at which levels to the
1122 ! metgrid.log file, and then we check to see if any fields are not
1123 ! completely covered with data.
1124 !
1125 call storage_print_fields()
1126 call find_missing_values()
1128 !
1129 ! All of the processing is now done for this time period for this domain;
1130 ! now we simply output every field from the storage module.
1131 !
1133 title = 'OUTPUT FROM METGRID V3.3.1'
1135 ! Initialize the output module for this domain and time
1136 call mprintf(.true.,LOGFILE,'Initializing output module.')
1137 output_date = temp_date
1138 if ( .not. nocolons ) then
1139 if (len_trim(temp_date) == 13) then
1140 output_date(14:19) = ':00:00'
1141 else if (len_trim(temp_date) == 16) then
1142 output_date(17:19) = ':00'
1143 end if
1144 else
1145 if (len_trim(temp_date) == 13) then
1146 output_date(14:19) = '_00_00'
1147 else if (len_trim(temp_date) == 16) then
1148 output_date(17:19) = '_00'
1149 end if
1150 endif
1152 call output_init(n, title, output_date, gridtype, dyn_opt, &
1153 corner_lats, corner_lons, &
1154 we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
1155 we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, &
1156 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
1157 extra_col, extra_row)
1159 call get_bottom_top_dim(bottom_top_dim)
1161 ! Add in a flag to tell real that we have seven new msf fields
1162 nmet_flags = num_entries + 1
1163 allocate(met_flags(nmet_flags))
1164 do i=1,num_entries
1165 met_flags(i) = output_flags(i)
1166 end do
1167 if (gridtype == 'C') then
1168 met_flags(num_entries+1) = 'FLAG_MF_XY'
1169 else
1170 met_flags(num_entries+1) = ' '
1171 end if
1173 ! Find out how many soil levels or layers we have; this assumes a field named ST
1174 field % header % field = 'ST'
1175 nullify(soil_levels)
1176 call storage_get_levels(field, soil_levels)
1178 if (.not. associated(soil_levels)) then
1179 field % header % field = 'SOILT'
1180 nullify(soil_levels)
1181 call storage_get_levels(field, soil_levels)
1182 end if
1184 if (.not. associated(soil_levels)) then
1185 field % header % field = 'STC_WPS'
1186 nullify(soil_levels)
1187 call storage_get_levels(field, soil_levels)
1188 end if
1190 if (associated(soil_levels)) then
1191 num_metgrid_soil_levs = size(soil_levels)
1192 deallocate(soil_levels)
1193 else
1194 num_metgrid_soil_levs = 0
1195 end if
1197 ! First write out global attributes
1198 call mprintf(.true.,LOGFILE,'Writing global attributes to output.')
1199 call write_global_attrs(title, output_date, gridtype, dyn_opt, west_east_dim, &
1200 south_north_dim, bottom_top_dim, &
1201 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
1202 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
1203 map_proj, mminlu, num_land_cat, &
1204 is_water, is_lake, is_ice, is_urban, i_soilwater, &
1205 grid_id, parent_id, i_parent_start, &
1206 j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
1207 cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
1208 pole_lat, pole_lon, parent_grid_ratio, sub_x, sub_y, &
1209 corner_lats, corner_lons, num_metgrid_soil_levs, &
1210 met_flags, nmet_flags, process_bdy_width)
1212 deallocate(met_flags)
1214 call reset_next_field()
1216 istatus = 0
1218 ! Now loop over all output fields, writing each to the output module
1219 do while (istatus == 0)
1220 call get_next_output_field(cname, real_array, &
1221 sm1, em1, sm2, em2, sm3, em3, istatus)
1222 if (istatus == 0) then
1224 call mprintf(.true.,LOGFILE,'Writing field %s to output.',s1=cname)
1225 call write_field(sm1, em1, sm2, em2, sm3, em3, &
1226 cname, output_date, real_array)
1227 deallocate(real_array)
1229 end if
1231 end do
1233 call mprintf(.true.,LOGFILE,'Closing output file.')
1234 call output_close()
1236 ! Free up memory used by met fields for this valid time
1237 call storage_delete_all_td()
1239 end subroutine process_single_met_time
1242 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1243 ! Name: get_interp_masks
1244 !
1245 ! Purpose:
1246 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1247 subroutine get_interp_masks(fg_prefix, is_constants, fg_date)
1249 use interp_option_module
1250 use read_met_module
1251 use storage_module
1253 implicit none
1255 ! Arguments
1256 logical, intent(in) :: is_constants
1257 character (len=*), intent(in) :: fg_prefix, fg_date
1259 ! BUG: Move this constant to misc_definitions_module?
1260 integer, parameter :: BDR_WIDTH = 3
1262 ! Local variables
1263 integer :: i, istatus, idx, idxt
1264 type (fg_input) :: mask_field
1265 type (met_data) :: fg_data
1267 istatus = 0
1269 call read_met_init(fg_prefix, is_constants, fg_date, istatus)
1271 do while (istatus == 0)
1273 call read_next_met_field(fg_data, istatus)
1275 if (istatus == 0) then
1277 ! Find out which METGRID.TBL entry goes with this field
1278 idxt = num_entries + 1
1279 do idx=1,num_entries
1280 if ((index(fieldname(idx), trim(fg_data%field)) /= 0) .and. &
1281 (len_trim(fieldname(idx)) == len_trim(fg_data%field))) then
1283 if (index(fg_prefix,trim(from_input(idx))) /= 0 .or. &
1284 (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
1285 idxt = idx
1286 end if
1288 end if
1289 end do
1290 idx = idxt
1291 if (idx > num_entries) idx = num_entries ! The last entry is a default
1293 ! Do we need to rename this field?
1294 if (output_name(idx) /= ' ') then
1295 fg_data%field = output_name(idx)(1:9)
1297 idxt = num_entries + 1
1298 do idx=1,num_entries
1299 if ((index(fieldname(idx), trim(fg_data%field)) /= 0) .and. &
1300 (len_trim(fieldname(idx)) == len_trim(fg_data%field))) then
1302 if (index(fg_prefix,trim(from_input(idx))) /= 0 .or. &
1303 (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
1304 idxt = idx
1305 end if
1307 end if
1308 end do
1309 idx = idxt
1310 if (idx > num_entries) idx = num_entries ! The last entry is a default
1311 end if
1313 do i=1,num_entries
1314 if (interp_mask(i) /= ' ' .and. (trim(interp_mask(i)) == trim(fg_data%field))) then
1316 mask_field%header%version = 1
1317 mask_field%header%date = ' '
1318 mask_field%header%date = fg_date
1319 if (is_constants) then
1320 mask_field%header%time_dependent = .false.
1321 else
1322 mask_field%header%time_dependent = .true.
1323 end if
1324 mask_field%header%mask_field = .true.
1325 mask_field%header%forecast_hour = 0.
1326 mask_field%header%fg_source = 'degribbed met data'
1327 mask_field%header%field = trim(fg_data%field)//'.mask'
1328 mask_field%header%units = '-'
1329 mask_field%header%description = '-'
1330 mask_field%header%vertical_coord = 'none'
1331 mask_field%header%vertical_level = 1
1332 mask_field%header%sr_x = 1
1333 mask_field%header%sr_y = 1
1334 mask_field%header%array_order = 'XY'
1335 mask_field%header%dim1(1) = 1
1336 mask_field%header%dim1(2) = fg_data%nx
1337 mask_field%header%dim2(1) = 1
1338 mask_field%header%dim2(2) = fg_data%ny
1339 mask_field%header%is_wind_grid_rel = .true.
1340 mask_field%header%array_has_missing_values = .false.
1341 mask_field%map%stagger = M
1343 ! Do a simple check to see whether this is a global lat/lon dataset
1344 ! Note that we do not currently support regional Gaussian grids
1345 if ((fg_data%iproj == PROJ_LATLON .and. abs(fg_data%nx * fg_data%deltalon - 360.) < 0.0001) &
1346 .or. (fg_data%iproj == PROJ_GAUSS)) then
1347 allocate(mask_field%r_arr(1-BDR_WIDTH:fg_data%nx+BDR_WIDTH,1:fg_data%ny))
1349 mask_field%r_arr(1:fg_data%nx, 1:fg_data%ny) = &
1350 fg_data%slab(1:fg_data%nx, 1:fg_data%ny)
1352 mask_field%r_arr(1-BDR_WIDTH:0, 1:fg_data%ny) = &
1353 fg_data%slab(fg_data%nx-BDR_WIDTH+1:fg_data%nx, 1:fg_data%ny)
1355 mask_field%r_arr(fg_data%nx+1:fg_data%nx+BDR_WIDTH, 1:fg_data%ny) = &
1356 fg_data%slab(1:BDR_WIDTH, 1:fg_data%ny)
1358 else
1359 allocate(mask_field%r_arr(1:fg_data%nx,1:fg_data%ny))
1360 mask_field%r_arr = fg_data%slab
1361 end if
1363 nullify(mask_field%valid_mask)
1364 nullify(mask_field%modified_mask)
1366 call storage_put_field(mask_field)
1368 exit
1370 end if
1371 end do
1373 if (associated(fg_data%slab)) deallocate(fg_data%slab)
1375 end if
1377 end do
1379 call read_met_close()
1381 end subroutine get_interp_masks
1384 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1385 ! Name: interp_met_field
1386 !
1387 ! Purpose:
1388 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1389 subroutine interp_met_field(input_name, short_fieldnm, ifieldstagger, istagger, &
1390 field, xlat, xlon, sm1, em1, sm2, em2, &
1391 sd1, ed1, sd2, ed2, &
1392 slab, minx, maxx, miny, maxy, bdr, do_gcell_interp, &
1393 new_pts, process_bdy_width, landmask)
1395 use bitarray_module
1396 use interp_module
1397 use interp_option_module
1398 use llxy_module
1399 use misc_definitions_module
1400 use storage_module
1402 implicit none
1404 ! Arguments
1405 integer, intent(in) :: ifieldstagger, istagger, &
1406 sm1, em1, sm2, em2, &
1407 sd1, ed1, sd2, ed2, &
1408 minx, maxx, miny, maxy, bdr, &
1409 process_bdy_width
1410 real, dimension(minx:maxx,miny:maxy), intent(in) :: slab
1411 real, dimension(sm1:em1,sm2:em2), intent(in) :: xlat, xlon
1412 real, dimension(sm1:em1,sm2:em2), intent(in), optional :: landmask
1413 logical, intent(in) :: do_gcell_interp
1414 character (len=9), intent(in) :: short_fieldnm
1415 character (len=MAX_FILENAME_LEN), intent(in) :: input_name
1416 type (fg_input), intent(inout) :: field
1417 type (bitarray), intent(inout) :: new_pts
1419 ! Local variables
1420 integer :: i, j, idx, idxt, orig_selected_proj, interp_mask_status, &
1421 interp_land_mask_status, interp_water_mask_status, process_width
1422 integer, pointer, dimension(:) :: interp_array
1423 real :: rx, ry, temp
1424 real, pointer, dimension(:,:) :: data_count
1425 type (fg_input) :: mask_field, mask_water_field, mask_land_field
1427 ! Find index into fieldname, interp_method, masked, and fill_missing
1428 ! of the current field
1429 idxt = num_entries + 1
1430 do idx=1,num_entries
1431 if ((index(fieldname(idx), trim(short_fieldnm)) /= 0) .and. &
1432 (len_trim(fieldname(idx)) == len_trim(short_fieldnm))) then
1433 if (index(input_name,trim(from_input(idx))) /= 0 .or. &
1434 (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
1435 idxt = idx
1436 end if
1437 end if
1438 end do
1439 idx = idxt
1440 if (idx > num_entries) then
1441 call mprintf(.true.,WARN,'Entry in METGRID.TBL not found for field %s. '// &
1442 'Default options will be used for this field!', s1=short_fieldnm)
1443 idx = num_entries ! The last entry is a default
1444 end if
1446 if (process_bdy_width == 0) then
1447 process_width = max(ed1-sd1+1, ed2-sd2+1)
1448 else
1449 process_width = process_bdy_width
1450 ! Add two extra rows/cols to accommodate staggered fields: one extra row/col for
1451 ! averaging to mass points in real, and one beyond that for averaging during
1452 ! wind rotation
1453 if (ifieldstagger /= M) process_width = process_width + 2
1454 end if
1456 field%header%dim1(1) = sm1
1457 field%header%dim1(2) = em1
1458 field%header%dim2(1) = sm2
1459 field%header%dim2(2) = em2
1460 field%map%stagger = ifieldstagger
1461 if (.not. associated(field%r_arr)) then
1462 allocate(field%r_arr(sm1:em1,sm2:em2))
1463 end if
1465 interp_mask_status = 1
1466 interp_land_mask_status = 1
1467 interp_water_mask_status = 1
1469 if (interp_mask(idx) /= ' ') then
1470 mask_field%header%version = 1
1471 mask_field%header%forecast_hour = 0.
1472 mask_field%header%field = trim(interp_mask(idx))//'.mask'
1473 mask_field%header%vertical_coord = 'none'
1474 mask_field%header%vertical_level = 1
1476 call storage_get_field(mask_field, interp_mask_status)
1478 end if
1479 if (interp_land_mask(idx) /= ' ') then
1480 mask_land_field%header%version = 1
1481 mask_land_field%header%forecast_hour = 0.
1482 mask_land_field%header%field = trim(interp_land_mask(idx))//'.mask'
1483 mask_land_field%header%vertical_coord = 'none'
1484 mask_land_field%header%vertical_level = 1
1486 call storage_get_field(mask_land_field, interp_land_mask_status)
1488 end if
1489 if (interp_water_mask(idx) /= ' ') then
1490 mask_water_field%header%version = 1
1491 mask_water_field%header%forecast_hour = 0.
1492 mask_water_field%header%field = trim(interp_water_mask(idx))//'.mask'
1493 mask_water_field%header%vertical_coord = 'none'
1494 mask_water_field%header%vertical_level = 1
1496 call storage_get_field(mask_water_field, interp_water_mask_status)
1498 end if
1500 interp_array => interp_array_from_string(interp_method(idx))
1503 !
1504 ! Interpolate using average_gcell interpolation method
1505 !
1506 if (do_gcell_interp) then
1507 allocate(data_count(sm1:em1,sm2:em2))
1508 data_count = 0.
1510 if (interp_mask_status == 0) then
1511 call accum_continuous(slab, &
1512 minx, maxx, miny, maxy, 1, 1, bdr, &
1513 field%r_arr, data_count, &
1514 sm1, em1, sm2, em2, 1, 1, &
1515 istagger, &
1516 new_pts, missing_value(idx), interp_mask_val(idx), interp_mask_relational(idx), mask_field%r_arr)
1517 else
1518 call accum_continuous(slab, &
1519 minx, maxx, miny, maxy, 1, 1, bdr, &
1520 field%r_arr, data_count, &
1521 sm1, em1, sm2, em2, 1, 1, &
1522 istagger, &
1523 new_pts, missing_value(idx), -1.) ! The -1 is the maskval, but since we
1524 ! we do not give an optional mask, no
1525 ! no need to worry about -1s in data
1526 end if
1528 orig_selected_proj = iget_selected_domain()
1529 call select_domain(SOURCE_PROJ)
1530 do j=sm2,em2
1531 do i=sm1,em1
1533 if (abs(i - sd1) >= process_width .and. (abs(i - ed1) >= process_width) .and. &
1534 abs(j - sd2) >= process_width .and. (abs(j - ed2) >= process_width)) then
1535 field%r_arr(i,j) = fill_missing(idx)
1536 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1537 cycle
1538 end if
1540 if (present(landmask)) then
1542 if (landmask(i,j) /= masked(idx)) then
1543 if (data_count(i,j) > 0.) then
1544 field%r_arr(i,j) = field%r_arr(i,j) / data_count(i,j)
1545 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1546 else
1548 if (interp_mask_status == 0) then
1549 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1550 minx, maxx, miny, maxy, bdr, missing_value(idx), &
1551 mask_relational=interp_mask_relational(idx), &
1552 mask_val=interp_mask_val(idx), mask_field=mask_field%r_arr)
1553 else
1554 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1555 minx, maxx, miny, maxy, bdr, missing_value(idx))
1556 end if
1558 if (temp /= missing_value(idx)) then
1559 field%r_arr(i,j) = temp
1560 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1561 end if
1563 end if
1564 else
1565 field%r_arr(i,j) = fill_missing(idx)
1566 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1567 end if
1569 if (.not. bitarray_test(new_pts, i-sm1+1, j-sm2+1) .and. &
1570 .not. bitarray_test(field%valid_mask, i-sm1+1, j-sm2+1)) then
1571 field%r_arr(i,j) = fill_missing(idx)
1573 ! Assume that if missing fill value is other than default, then user has asked
1574 ! to fill in any missing values, and we can consider this point to have
1575 ! received a valid value
1576 if (fill_missing(idx) /= NAN) call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1577 end if
1579 else
1581 if (data_count(i,j) > 0.) then
1582 field%r_arr(i,j) = field%r_arr(i,j) / data_count(i,j)
1583 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1584 else
1586 if (interp_mask_status == 0) then
1587 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1588 minx, maxx, miny, maxy, bdr, missing_value(idx), &
1589 mask_relational=interp_mask_relational(idx), &
1590 mask_val=interp_mask_val(idx), mask_field=mask_field%r_arr)
1591 else
1592 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1593 minx, maxx, miny, maxy, bdr, missing_value(idx))
1594 end if
1596 if (temp /= missing_value(idx)) then
1597 field%r_arr(i,j) = temp
1598 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1599 end if
1601 end if
1603 if (.not. bitarray_test(new_pts, i-sm1+1, j-sm2+1) .and. &
1604 .not. bitarray_test(field%valid_mask, i-sm1+1, j-sm2+1)) then
1605 field%r_arr(i,j) = fill_missing(idx)
1607 ! Assume that if missing fill value is other than default, then user has asked
1608 ! to fill in any missing values, and we can consider this point to have
1609 ! received a valid value
1610 if (fill_missing(idx) /= NAN) call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1611 end if
1613 end if
1615 end do
1616 end do
1617 call select_domain(orig_selected_proj)
1618 deallocate(data_count)
1620 !
1621 ! No average_gcell interpolation method
1622 !
1623 else
1625 orig_selected_proj = iget_selected_domain()
1626 call select_domain(SOURCE_PROJ)
1627 do j=sm2,em2
1628 do i=sm1,em1
1630 if (abs(i - sd1) >= process_width .and. (abs(i - ed1) >= process_width) .and. &
1631 abs(j - sd2) >= process_width .and. (abs(j - ed2) >= process_width)) then
1632 field%r_arr(i,j) = fill_missing(idx)
1633 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1634 cycle
1635 end if
1637 if (present(landmask)) then
1639 if (masked(idx) == MASKED_BOTH) then
1641 if (landmask(i,j) == 0) then ! WATER POINT
1643 if (interp_land_mask_status == 0) then
1644 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1645 minx, maxx, miny, maxy, bdr, missing_value(idx), &
1646 mask_relational=interp_land_mask_relational(idx), &
1647 mask_val=interp_land_mask_val(idx), mask_field=mask_land_field%r_arr)
1648 else
1649 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1650 minx, maxx, miny, maxy, bdr, missing_value(idx))
1651 end if
1653 else if (landmask(i,j) == 1) then ! LAND POINT
1655 if (interp_water_mask_status == 0) then
1656 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1657 minx, maxx, miny, maxy, bdr, missing_value(idx), &
1658 mask_relational=interp_water_mask_relational(idx), &
1659 mask_val=interp_water_mask_val(idx), mask_field=mask_water_field%r_arr)
1660 else
1661 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1662 minx, maxx, miny, maxy, bdr, missing_value(idx))
1663 end if
1665 end if
1667 else if (landmask(i,j) /= masked(idx)) then
1669 if (interp_mask_status == 0) then
1670 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1671 minx, maxx, miny, maxy, bdr, missing_value(idx), &
1672 mask_relational=interp_mask_relational(idx), &
1673 mask_val=interp_mask_val(idx), mask_field=mask_field%r_arr)
1674 else
1675 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1676 minx, maxx, miny, maxy, bdr, missing_value(idx))
1677 end if
1679 else
1680 temp = missing_value(idx)
1681 end if
1683 ! No landmask for this field
1684 else
1686 if (interp_mask_status == 0) then
1687 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1688 minx, maxx, miny, maxy, bdr, missing_value(idx), &
1689 mask_relational=interp_mask_relational(idx), &
1690 mask_val=interp_mask_val(idx), mask_field=mask_field%r_arr)
1691 else
1692 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1693 minx, maxx, miny, maxy, bdr, missing_value(idx))
1694 end if
1696 end if
1698 if (temp /= missing_value(idx)) then
1699 field%r_arr(i,j) = temp
1700 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1701 else if (present(landmask)) then
1702 if (landmask(i,j) == masked(idx)) then
1703 field%r_arr(i,j) = fill_missing(idx)
1704 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1705 end if
1706 end if
1708 if (.not. bitarray_test(new_pts, i-sm1+1, j-sm2+1) .and. &
1709 .not. bitarray_test(field%valid_mask, i-sm1+1, j-sm2+1)) then
1710 field%r_arr(i,j) = fill_missing(idx)
1712 ! Assume that if missing fill value is other than default, then user has asked
1713 ! to fill in any missing values, and we can consider this point to have
1714 ! received a valid value
1715 if (fill_missing(idx) /= NAN) call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1716 end if
1718 end do
1719 end do
1720 call select_domain(orig_selected_proj)
1721 end if
1723 deallocate(interp_array)
1725 end subroutine interp_met_field
1728 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1729 ! Name: interp_to_latlon
1730 !
1731 ! Purpose:
1732 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1733 function interp_to_latlon(rlat, rlon, istagger, interp_method_list, slab, &
1734 minx, maxx, miny, maxy, bdr, source_missing_value, &
1735 mask_field, mask_relational, mask_val)
1737 use interp_module
1738 use llxy_module
1740 implicit none
1742 ! Arguments
1743 integer, intent(in) :: minx, maxx, miny, maxy, bdr, istagger
1744 integer, dimension(:), intent(in) :: interp_method_list
1745 real, intent(in) :: rlat, rlon, source_missing_value
1746 real, dimension(minx:maxx,miny:maxy), intent(in) :: slab
1747 real, intent(in), optional :: mask_val
1748 real, dimension(minx:maxx,miny:maxy), intent(in), optional :: mask_field
1749 character(len=1), intent(in), optional :: mask_relational
1751 ! Return value
1752 real :: interp_to_latlon
1754 ! Local variables
1755 real :: rx, ry
1757 interp_to_latlon = source_missing_value
1759 call lltoxy(rlat, rlon, rx, ry, istagger)
1760 if (rx >= minx+bdr-0.5 .and. rx <= maxx-bdr+0.5) then
1761 if (present(mask_field) .and. present(mask_val) .and. present (mask_relational)) then
1762 interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, 1, 1, source_missing_value, &
1763 interp_method_list, 1, mask_relational, mask_val, mask_field)
1764 else if (present(mask_field) .and. present(mask_val)) then
1765 interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, 1, 1, source_missing_value, &
1766 interp_method_list, 1, maskval=mask_val, mask_array=mask_field)
1767 else
1768 interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, 1, 1, source_missing_value, &
1769 interp_method_list, 1)
1770 end if
1771 else
1772 interp_to_latlon = source_missing_value
1773 end if
1775 if (interp_to_latlon == source_missing_value) then
1777 ! Try a lon in the range 0. to 360.; all lons in the xlon
1778 ! array should be in the range -180. to 180.
1779 if (rlon < 0.) then
1780 call lltoxy(rlat, rlon+360., rx, ry, istagger)
1781 if (rx >= minx+bdr-0.5 .and. rx <= maxx-bdr+0.5) then
1782 if (present(mask_field) .and. present(mask_val) .and. present(mask_relational)) then
1783 interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, &
1784 1, 1, source_missing_value, &
1785 interp_method_list, 1, mask_relational, mask_val, mask_field)
1786 else if (present(mask_field) .and. present(mask_val)) then
1787 interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, &
1788 1, 1, source_missing_value, &
1789 interp_method_list, 1, maskval=mask_val, mask_array=mask_field)
1790 else
1791 interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, &
1792 1, 1, source_missing_value, &
1793 interp_method_list, 1)
1794 end if
1795 else
1796 interp_to_latlon = source_missing_value
1797 end if
1799 end if
1801 end if
1803 return
1805 end function interp_to_latlon
1808 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1809 ! Name: get_bottom_top_dim
1810 !
1811 ! Purpose:
1812 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1813 subroutine get_bottom_top_dim(bottom_top_dim)
1815 use interp_option_module
1816 use list_module
1817 use storage_module
1819 implicit none
1821 ! Arguments
1822 integer, intent(out) :: bottom_top_dim
1824 ! Local variables
1825 integer :: i, j
1826 integer, pointer, dimension(:) :: field_levels
1827 character (len=32) :: z_dim
1828 type (fg_input), pointer, dimension(:) :: headers
1829 type (list) :: temp_levels
1831 ! Initialize a list to store levels that are found for 3-d fields
1832 call list_init(temp_levels)
1834 ! Get a list of all time-dependent fields (given by their headers) from
1835 ! the storage module
1836 call storage_get_td_headers(headers)
1838 !
1839 ! Given headers of all fields, we first build a list of all possible levels
1840 ! for 3-d met fields (excluding sea-level, though).
1841 !
1842 do i=1,size(headers)
1843 call get_z_dim_name(headers(i)%header%field, z_dim)
1845 ! We only want to consider 3-d met fields
1846 if (z_dim(1:18) == 'num_metgrid_levels') then
1848 ! Find out what levels the current field has
1849 call storage_get_levels(headers(i), field_levels)
1850 do j=1,size(field_levels)
1852 ! If this level has not yet been encountered, add it to our list
1853 if (.not. list_search(temp_levels, ikey=field_levels(j), ivalue=field_levels(j))) then
1854 if (field_levels(j) /= 201300) then
1855 call list_insert(temp_levels, ikey=field_levels(j), ivalue=field_levels(j))
1856 end if
1857 end if
1859 end do
1861 deallocate(field_levels)
1863 end if
1865 end do
1867 bottom_top_dim = list_length(temp_levels)
1869 call list_destroy(temp_levels)
1870 deallocate(headers)
1872 end subroutine get_bottom_top_dim
1875 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1876 ! Name: fill_missing_levels
1877 !
1878 ! Purpose:
1879 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1880 subroutine fill_missing_levels(output_flags)
1882 use interp_option_module
1883 use list_module
1884 use module_debug
1885 use module_mergesort
1886 use storage_module
1888 implicit none
1890 ! Arguments
1891 character (len=128), dimension(:), intent(inout) :: output_flags
1893 ! Local variables
1894 integer :: i, ii, j, ix, jx, k, lower, upper, temp, istatus
1895 integer, pointer, dimension(:) :: union_levels, field_levels
1896 real, pointer, dimension(:) :: r_union_levels
1897 character (len=128) :: clevel
1898 type (fg_input) :: lower_field, upper_field, new_field, search_field
1899 type (fg_input), pointer, dimension(:) :: headers, all_headers
1900 type (list) :: temp_levels
1901 type (list_item), pointer, dimension(:) :: keys
1903 ! Initialize a list to store levels that are found for 3-d fields
1904 call list_init(temp_levels)
1906 ! Get a list of all fields (given by their headers) from the storage module
1907 call storage_get_td_headers(headers)
1908 call storage_get_all_headers(all_headers)
1910 !
1911 ! Given headers of all fields, we first build a list of all possible levels
1912 ! for 3-d met fields (excluding sea-level, though).
1913 !
1914 do i=1,size(headers)
1916 ! Find out what levels the current field has
1917 call storage_get_levels(headers(i), field_levels)
1918 do j=1,size(field_levels)
1920 ! If this level has not yet been encountered, add it to our list
1921 if (.not. list_search(temp_levels, ikey=field_levels(j), ivalue=field_levels(j))) then
1922 if (field_levels(j) /= 201300) then
1923 call list_insert(temp_levels, ikey=field_levels(j), ivalue=field_levels(j))
1924 end if
1925 end if
1927 end do
1929 deallocate(field_levels)
1931 end do
1933 if (list_length(temp_levels) > 0) then
1935 !
1936 ! With all possible levels stored in a list, get an array of levels, sorted
1937 ! in decreasing order
1938 !
1939 i = 0
1940 allocate(union_levels(list_length(temp_levels)))
1941 do while (list_length(temp_levels) > 0)
1942 i = i + 1
1943 call list_get_first_item(temp_levels, ikey=union_levels(i), ivalue=temp)
1944 end do
1945 call mergesort(union_levels, 1, size(union_levels))
1947 allocate(r_union_levels(size(union_levels)))
1948 do i=1,size(union_levels)
1949 r_union_levels(i) = real(union_levels(i))
1950 end do
1952 !
1953 ! With a sorted, complete list of levels, we need
1954 ! to go back and fill in missing levels for each 3-d field
1955 !
1956 do i=1,size(headers)
1958 !
1959 ! Find entry in METGRID.TBL for this field, if one exists; if it does, then the
1960 ! entry may tell us how to get values for the current field at the missing level
1961 !
1962 do ii=1,num_entries
1963 if (fieldname(ii) == headers(i)%header%field) exit
1964 end do
1965 if (ii <= num_entries) then
1966 call dup(headers(i),new_field)
1967 nullify(new_field%valid_mask)
1968 nullify(new_field%modified_mask)
1969 call fill_field(new_field, ii, output_flags, r_union_levels)
1970 end if
1972 end do
1974 deallocate(union_levels)
1975 deallocate(r_union_levels)
1976 deallocate(headers)
1978 call storage_get_td_headers(headers)
1980 !
1981 ! Now we may need to vertically interpolate to missing values in 3-d fields
1982 !
1983 do i=1,size(headers)
1985 call storage_get_levels(headers(i), field_levels)
1987 ! If this isn't a 3-d array, nothing to do
1988 if (size(field_levels) > 1) then
1990 do k=1,size(field_levels)
1991 call dup(headers(i),search_field)
1992 search_field%header%vertical_level = field_levels(k)
1993 call storage_get_field(search_field,istatus)
1994 if (istatus == 0) then
1995 JLOOP: do jx=search_field%header%dim2(1),search_field%header%dim2(2)
1996 ILOOP: do ix=search_field%header%dim1(1),search_field%header%dim1(2)
1997 if (.not. bitarray_test(search_field%valid_mask, &
1998 ix-search_field%header%dim1(1)+1, &
1999 jx-search_field%header%dim2(1)+1)) then
2001 call dup(search_field, lower_field)
2002 do lower=k-1,1,-1
2003 lower_field%header%vertical_level = field_levels(lower)
2004 call storage_get_field(lower_field,istatus)
2005 if (bitarray_test(lower_field%valid_mask, &
2006 ix-search_field%header%dim1(1)+1, &
2007 jx-search_field%header%dim2(1)+1)) &
2008 exit
2010 end do
2012 call dup(search_field, upper_field)
2013 do upper=k+1,size(field_levels)
2014 upper_field%header%vertical_level = field_levels(upper)
2015 call storage_get_field(upper_field,istatus)
2016 if (bitarray_test(upper_field%valid_mask, &
2017 ix-search_field%header%dim1(1)+1, &
2018 jx-search_field%header%dim2(1)+1)) &
2019 exit
2021 end do
2022 if (upper <= size(field_levels) .and. lower >= 1) then
2023 search_field%r_arr(ix,jx) = real(abs(field_levels(upper)-field_levels(k))) &
2024 / real(abs(field_levels(upper)-field_levels(lower))) &
2025 * lower_field%r_arr(ix,jx) &
2026 + real(abs(field_levels(k)-field_levels(lower))) &
2027 / real(abs(field_levels(upper)-field_levels(lower))) &
2028 * upper_field%r_arr(ix,jx)
2029 call bitarray_set(search_field%valid_mask, &
2030 ix-search_field%header%dim1(1)+1, &
2031 jx-search_field%header%dim2(1)+1)
2032 end if
2033 end if
2034 end do ILOOP
2035 end do JLOOP
2036 else
2037 call mprintf(.true.,ERROR, &
2038 'This is bad, could not get %s at level %i.', &
2039 s1=trim(search_field%header%field), i1=field_levels(k))
2040 end if
2041 end do
2043 end if
2045 deallocate(field_levels)
2047 end do
2049 end if
2051 call list_destroy(temp_levels)
2052 deallocate(all_headers)
2053 deallocate(headers)
2055 end subroutine fill_missing_levels
2058 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2059 ! Name: create_derived_fields
2060 !
2061 ! Purpose:
2062 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2063 subroutine create_derived_fields(arg_gridtype, hdate, xfcst, &
2064 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
2065 we_mem_stag_s, we_mem_stag_e, sn_mem_stag_s, sn_mem_stag_e, &
2066 created_this_field, output_flags)
2068 use interp_option_module
2069 use list_module
2070 use module_mergesort
2071 use storage_module
2073 implicit none
2075 ! Arguments
2076 integer, intent(in) :: we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
2077 we_mem_stag_s, we_mem_stag_e, sn_mem_stag_s, sn_mem_stag_e
2078 real, intent(in) :: xfcst
2079 logical, dimension(:), intent(inout) :: created_this_field
2080 character (len=1), intent(in) :: arg_gridtype
2081 character (len=24), intent(in) :: hdate
2082 character (len=128), dimension(:), intent(inout) :: output_flags
2084 ! Local variables
2085 integer :: idx, i, j, istatus
2086 type (fg_input) :: field
2088 ! Initialize fg_input structure to store the field
2089 field%header%version = 1
2090 field%header%date = hdate//' '
2091 field%header%time_dependent = .true.
2092 field%header%mask_field = .false.
2093 field%header%forecast_hour = xfcst
2094 field%header%fg_source = 'Derived from FG'
2095 field%header%field = ' '
2096 field%header%units = ' '
2097 field%header%description = ' '
2098 field%header%vertical_level = 0
2099 field%header%sr_x = 1
2100 field%header%sr_y = 1
2101 field%header%array_order = 'XY '
2102 field%header%is_wind_grid_rel = .true.
2103 field%header%array_has_missing_values = .false.
2104 nullify(field%r_arr)
2105 nullify(field%valid_mask)
2106 nullify(field%modified_mask)
2108 !
2109 ! Check each entry in METGRID.TBL to see whether it is a derive field
2110 !
2111 do idx=1,num_entries
2112 if (is_derived_field(idx)) then
2114 created_this_field(idx) = .true.
2116 call mprintf(.true.,INFORM,'Going to create the field %s',s1=fieldname(idx))
2118 ! Intialize more fields in storage structure
2119 field%header%field = fieldname(idx)
2120 call get_z_dim_name(fieldname(idx),field%header%vertical_coord)
2121 field%map%stagger = output_stagger(idx)
2122 if (arg_gridtype == 'E') then
2123 field%header%dim1(1) = we_mem_s
2124 field%header%dim1(2) = we_mem_e
2125 field%header%dim2(1) = sn_mem_s
2126 field%header%dim2(2) = sn_mem_e
2127 else if (arg_gridtype == 'C') then
2128 if (output_stagger(idx) == M) then
2129 field%header%dim1(1) = we_mem_s
2130 field%header%dim1(2) = we_mem_e
2131 field%header%dim2(1) = sn_mem_s
2132 field%header%dim2(2) = sn_mem_e
2133 else if (output_stagger(idx) == U) then
2134 field%header%dim1(1) = we_mem_stag_s
2135 field%header%dim1(2) = we_mem_stag_e
2136 field%header%dim2(1) = sn_mem_s
2137 field%header%dim2(2) = sn_mem_e
2138 else if (output_stagger(idx) == V) then
2139 field%header%dim1(1) = we_mem_s
2140 field%header%dim1(2) = we_mem_e
2141 field%header%dim2(1) = sn_mem_stag_s
2142 field%header%dim2(2) = sn_mem_stag_e
2143 end if
2144 end if
2146 call fill_field(field, idx, output_flags)
2148 end if
2149 end do
2152 end subroutine create_derived_fields
2155 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2156 ! Name: fill_field
2157 !
2158 ! Purpose:
2159 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2160 subroutine fill_field(field, idx, output_flags, all_level_list)
2162 use interp_option_module
2163 use list_module
2164 use module_mergesort
2165 use storage_module
2167 implicit none
2169 ! Arguments
2170 integer, intent(in) :: idx
2171 type (fg_input), intent(inout) :: field
2172 character (len=128), dimension(:), intent(inout) :: output_flags
2173 real, dimension(:), intent(in), optional :: all_level_list
2175 ! Local variables
2176 integer :: i, j, istatus, isrclevel
2177 integer, pointer, dimension(:) :: all_list
2178 real :: rfillconst, rlevel, rsrclevel
2179 type (fg_input) :: query_field
2180 type (list_item), pointer, dimension(:) :: keys
2181 character (len=128) :: asrcname
2182 logical :: filled_all_lev
2184 filled_all_lev = .false.
2186 !
2187 ! Get a list of all levels to be filled for this field
2188 !
2189 keys => list_get_keys(fill_lev_list(idx))
2191 do i=1,list_length(fill_lev_list(idx))
2193 !
2194 ! First handle a specification for levels "all"
2195 !
2196 if (trim(keys(i)%ckey) == 'all') then
2198 ! We only want to fill all levels if we haven't already filled "all" of them
2199 if (.not. filled_all_lev) then
2201 filled_all_lev = .true.
2203 query_field%header%time_dependent = .true.
2204 query_field%header%field = ' '
2205 nullify(query_field%r_arr)
2206 nullify(query_field%valid_mask)
2207 nullify(query_field%modified_mask)
2209 ! See if we are filling this level with a constant
2210 call get_constant_fill_lev(keys(i)%cvalue, rfillconst, istatus)
2211 if (istatus == 0) then
2212 if (present(all_level_list)) then
2213 do j=1,size(all_level_list)
2214 call create_level(field, real(all_level_list(j)), idx, output_flags, rfillconst=rfillconst)
2215 end do
2216 else
2217 query_field%header%field = level_template(idx)
2218 nullify(all_list)
2219 call storage_get_levels(query_field, all_list)
2220 if (associated(all_list)) then
2221 do j=1,size(all_list)
2222 call create_level(field, real(all_list(j)), idx, output_flags, rfillconst=rfillconst)
2223 end do
2224 deallocate(all_list)
2225 end if
2226 end if
2228 ! Else see if we are filling this level with a constant equal
2229 ! to the value of the level
2230 else if (trim(keys(i)%cvalue) == 'vertical_index') then
2231 if (present(all_level_list)) then
2232 do j=1,size(all_level_list)
2233 call create_level(field, real(all_level_list(j)), idx, output_flags, &
2234 rfillconst=real(all_level_list(j)))
2235 end do
2236 else
2237 query_field%header%field = level_template(idx)
2238 nullify(all_list)
2239 call storage_get_levels(query_field, all_list)
2240 if (associated(all_list)) then
2241 do j=1,size(all_list)
2242 call create_level(field, real(all_list(j)), idx, output_flags, rfillconst=real(all_list(j)))
2243 end do
2244 deallocate(all_list)
2245 end if
2246 end if
2248 ! Else, we assume that it is a field from which we are copying levels
2249 else
2250 if (present(all_level_list)) then
2251 do j=1,size(all_level_list)
2252 call create_level(field, real(all_level_list(j)), idx, output_flags, &
2253 asrcname=keys(i)%cvalue, rsrclevel=real(all_level_list(j)))
2254 end do
2255 else
2256 query_field%header%field = keys(i)%cvalue ! Use same levels as source field, not level_template
2257 nullify(all_list)
2258 call storage_get_levels(query_field, all_list)
2259 if (associated(all_list)) then
2260 do j=1,size(all_list)
2261 call create_level(field, real(all_list(j)), idx, output_flags, &
2262 asrcname=keys(i)%cvalue, rsrclevel=real(all_list(j)))
2263 end do
2264 deallocate(all_list)
2266 else
2268 ! If the field doesn't have any levels (or does not exist) then we have not
2269 ! really filled all levels at this point.
2270 filled_all_lev = .false.
2271 end if
2272 end if
2274 end if
2275 end if
2277 !
2278 ! Handle individually specified levels
2279 !
2280 else
2282 read(keys(i)%ckey,*) rlevel
2284 ! See if we are filling this level with a constant
2285 call get_constant_fill_lev(keys(i)%cvalue, rfillconst, istatus)
2286 if (istatus == 0) then
2287 call create_level(field, rlevel, idx, output_flags, rfillconst=rfillconst)
2289 ! Otherwise, we are filling from another level
2290 else
2291 call get_fill_src_level(keys(i)%cvalue, asrcname, isrclevel)
2292 rsrclevel = real(isrclevel)
2293 call create_level(field, rlevel, idx, output_flags, &
2294 asrcname=asrcname, rsrclevel=rsrclevel)
2296 end if
2297 end if
2298 end do
2300 if (associated(keys)) deallocate(keys)
2302 end subroutine fill_field
2305 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2306 ! Name: create_level
2307 !
2308 ! Purpose:
2309 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2310 subroutine create_level(field_template, rlevel, idx, output_flags, &
2311 rfillconst, asrcname, rsrclevel)
2313 use storage_module
2314 use interp_option_module
2316 implicit none
2318 ! Arguments
2319 type (fg_input), intent(inout) :: field_template
2320 real, intent(in) :: rlevel
2321 integer, intent(in) :: idx
2322 character (len=128), dimension(:), intent(inout) :: output_flags
2323 real, intent(in), optional :: rfillconst, rsrclevel
2324 character (len=128), intent(in), optional :: asrcname
2326 ! Local variables
2327 integer :: i, j, istatus
2328 integer :: sm1,em1,sm2,em2
2329 type (fg_input) :: query_field
2331 !
2332 ! Check to make sure optional arguments are sane
2333 !
2334 if (present(rfillconst) .and. (present(asrcname) .or. present(rsrclevel))) then
2335 call mprintf(.true.,ERROR,'A call to create_level() cannot be given specifications '// &
2336 'for both a constant fill value and a source level.')
2338 else if ((present(asrcname) .and. .not. present(rsrclevel)) .or. &
2339 (.not. present(asrcname) .and. present(rsrclevel))) then
2340 call mprintf(.true.,ERROR,'Neither or both of optional arguments asrcname and '// &
2341 'rsrclevel must be specified to subroutine create_level().')
2343 else if (.not. present(rfillconst) .and. &
2344 .not. present(asrcname) .and. &
2345 .not. present(rsrclevel)) then
2346 call mprintf(.true.,ERROR,'A call to create_level() must be given either a specification '// &
2347 'for a constant fill value or a source level.')
2348 end if
2350 query_field%header%time_dependent = .true.
2351 query_field%header%field = field_template%header%field
2352 query_field%header%vertical_level = rlevel
2353 nullify(query_field%r_arr)
2354 nullify(query_field%valid_mask)
2355 nullify(query_field%modified_mask)
2357 call storage_query_field(query_field, istatus)
2358 if (istatus == 0) then
2359 call mprintf(.true.,INFORM,'%s at level %f already exists; leaving it alone.', &
2360 s1=field_template%header%field, f1=rlevel)
2361 return
2362 end if
2364 sm1 = field_template%header%dim1(1)
2365 em1 = field_template%header%dim1(2)
2366 sm2 = field_template%header%dim2(1)
2367 em2 = field_template%header%dim2(2)
2369 !
2370 ! Handle constant fill value case
2371 !
2372 if (present(rfillconst)) then
2374 field_template%header%vertical_level = rlevel
2375 allocate(field_template%r_arr(sm1:em1,sm2:em2))
2376 allocate(field_template%valid_mask)
2377 allocate(field_template%modified_mask)
2378 call bitarray_create(field_template%valid_mask, em1-sm1+1, em2-sm2+1)
2379 call bitarray_create(field_template%modified_mask, em1-sm1+1, em2-sm2+1)
2381 field_template%r_arr = rfillconst
2383 do j=sm2,em2
2384 do i=sm1,em1
2385 call bitarray_set(field_template%valid_mask, i-sm1+1, j-sm2+1)
2386 end do
2387 end do
2389 call storage_put_field(field_template)
2391 if (output_this_field(idx) .and. flag_in_output(idx) /= ' ') then
2392 output_flags(idx) = flag_in_output(idx)
2393 end if
2395 !
2396 ! Handle source field and source level case
2397 !
2398 else if (present(asrcname) .and. present(rsrclevel)) then
2400 query_field%header%field = ' '
2401 query_field%header%field = asrcname
2402 query_field%header%vertical_level = rsrclevel
2404 ! Check to see whether the requested source field exists at the requested level
2405 call storage_query_field(query_field, istatus)
2407 if (istatus == 0) then
2409 ! Read in requested field at requested level
2410 call storage_get_field(query_field, istatus)
2411 if ((query_field%header%dim1(1) /= field_template%header%dim1(1)) .or. &
2412 (query_field%header%dim1(2) /= field_template%header%dim1(2)) .or. &
2413 (query_field%header%dim2(1) /= field_template%header%dim2(1)) .or. &
2414 (query_field%header%dim2(2) /= field_template%header%dim2(2))) then
2415 call mprintf(.true.,ERROR,'Dimensions for %s do not match those of %s. This is '// &
2416 'probably because the staggerings of the fields do not match.', &
2417 s1=query_field%header%field, s2=field_template%header%field)
2418 end if
2420 field_template%header%vertical_level = rlevel
2421 allocate(field_template%r_arr(sm1:em1,sm2:em2))
2422 allocate(field_template%valid_mask)
2423 allocate(field_template%modified_mask)
2424 call bitarray_create(field_template%valid_mask, em1-sm1+1, em2-sm2+1)
2425 call bitarray_create(field_template%modified_mask, em1-sm1+1, em2-sm2+1)
2427 field_template%r_arr = query_field%r_arr
2429 ! We should retain information about which points in the field are valid
2430 do j=sm2,em2
2431 do i=sm1,em1
2432 if (bitarray_test(query_field%valid_mask, i-sm1+1, j-sm2+1)) then
2433 call bitarray_set(field_template%valid_mask, i-sm1+1, j-sm2+1)
2434 end if
2435 end do
2436 end do
2438 call storage_put_field(field_template)
2440 if (output_this_field(idx) .and. flag_in_output(idx) /= ' ') then
2441 output_flags(idx) = flag_in_output(idx)
2442 end if
2444 else
2445 call mprintf(.true.,INFORM,'Couldn''t find %s at level %f to fill level %f of %s.', &
2446 s1=asrcname,f1=rsrclevel,f2=rlevel,s2=field_template%header%field)
2447 end if
2449 end if
2451 end subroutine create_level
2454 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2455 ! Name: accum_continuous
2456 !
2457 ! Purpose: Sum up all of the source data points whose nearest neighbor in the
2458 ! model grid is the specified model grid point.
2459 !
2460 ! NOTE: When processing the source tile, those source points that are
2461 ! closest to a different model grid point will be added to the totals for
2462 ! such grid points; thus, an entire source tile will be processed at a time.
2463 ! This routine really processes for all model grid points that are
2464 ! within a source tile, and not just for a single grid point.
2465 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2466 subroutine accum_continuous(src_array, &
2467 src_min_x, src_max_x, src_min_y, src_max_y, src_min_z, src_max_z, bdr_width, &
2468 dst_array, n, &
2469 start_i, end_i, start_j, end_j, start_k, end_k, &
2470 istagger, &
2471 new_pts, msgval, maskval, mask_relational, mask_array, sr_x, sr_y)
2473 use bitarray_module
2474 use misc_definitions_module
2476 implicit none
2478 ! Arguments
2479 integer, intent(in) :: start_i, end_i, start_j, end_j, start_k, end_k, istagger, &
2480 src_min_x, src_max_x, src_min_y, src_max_y, src_min_z, src_max_z, bdr_width
2481 real, intent(in) :: maskval, msgval
2482 real, dimension(src_min_x:src_max_x, src_min_y:src_max_y, src_min_z:src_max_z), intent(in) :: src_array
2483 real, dimension(start_i:end_i, start_j:end_j, start_k:end_k), intent(inout) :: dst_array, n
2484 real, dimension(src_min_x:src_max_x, src_min_y:src_max_y), intent(in), optional :: mask_array
2485 integer, intent(in), optional :: sr_x, sr_y
2486 type (bitarray), intent(inout) :: new_pts
2487 character(len=1), intent(in), optional :: mask_relational
2489 ! Local variables
2490 integer :: i, j
2491 integer, pointer, dimension(:,:,:) :: where_maps_to
2492 real :: rsr_x, rsr_y
2494 rsr_x = 1.0
2495 rsr_y = 1.0
2496 if (present(sr_x)) rsr_x = real(sr_x)
2497 if (present(sr_y)) rsr_y = real(sr_y)
2499 allocate(where_maps_to(src_min_x:src_max_x,src_min_y:src_max_y,2))
2500 do i=src_min_x,src_max_x
2501 do j=src_min_y,src_max_y
2502 where_maps_to(i,j,1) = NOT_PROCESSED
2503 end do
2504 end do
2506 call process_continuous_block(src_array, where_maps_to, &
2507 src_min_x, src_min_y, src_min_z, src_max_x, src_max_y, src_max_z, &
2508 src_min_x+bdr_width, src_min_y, src_min_z, &
2509 src_max_x-bdr_width, src_max_y, src_max_z, &
2510 dst_array, n, start_i, end_i, start_j, end_j, start_k, end_k, &
2511 istagger, &
2512 new_pts, rsr_x, rsr_y, msgval, maskval, mask_relational, mask_array)
2514 deallocate(where_maps_to)
2516 end subroutine accum_continuous
2519 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2520 ! Name: process_continuous_block
2521 !
2522 ! Purpose: To recursively process a subarray of continuous data, adding the
2523 ! points in a block to the sum for their nearest grid point. The nearest
2524 ! neighbor may be estimated in some cases; for example, if the four corners
2525 ! of a subarray all have the same nearest grid point, all elements in the
2526 ! subarray are added to that grid point.
2527 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2528 recursive subroutine process_continuous_block(tile_array, where_maps_to, &
2529 src_min_x, src_min_y, src_min_z, src_max_x, src_max_y, src_max_z, &
2530 min_i, min_j, min_k, max_i, max_j, max_k, &
2531 dst_array, n, &
2532 start_x, end_x, start_y, end_y, start_z, end_z, &
2533 istagger, &
2534 new_pts, sr_x, sr_y, msgval, maskval, mask_relational, mask_array)
2536 use bitarray_module
2537 use llxy_module
2538 use misc_definitions_module
2540 implicit none
2542 ! Arguments
2543 integer, intent(in) :: min_i, min_j, min_k, max_i, max_j, max_k, &
2544 src_min_x, src_min_y, src_min_z, src_max_x, src_max_y, src_max_z, &
2545 start_x, end_x, start_y, end_y, start_z, end_z, istagger
2546 integer, dimension(src_min_x:src_max_x,src_min_y:src_max_y,2), intent(inout) :: where_maps_to
2547 real, intent(in) :: sr_x, sr_y, maskval, msgval
2548 real, dimension(src_min_x:src_max_x,src_min_y:src_max_y,src_min_z:src_max_z), intent(in) :: tile_array
2549 real, dimension(src_min_x:src_max_x,src_min_y:src_max_y), intent(in), optional :: mask_array
2550 real, dimension(start_x:end_x,start_y:end_y,start_z:end_z), intent(inout) :: dst_array, n
2551 type (bitarray), intent(inout) :: new_pts
2552 character(len=1), intent(in), optional :: mask_relational
2554 ! Local variables
2555 integer :: orig_selected_domain, x_dest, y_dest, i, j, k, center_i, center_j
2556 real :: lat_corner, lon_corner, rx, ry
2558 ! Compute the model grid point that the corners of the rectangle to be
2559 ! processed map to
2560 ! Lower-left corner
2561 if (where_maps_to(min_i,min_j,1) == NOT_PROCESSED) then
2562 orig_selected_domain = iget_selected_domain()
2563 call select_domain(SOURCE_PROJ)
2564 call xytoll(real(min_i), real(min_j), lat_corner, lon_corner, istagger)
2565 call select_domain(1)
2566 call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
2567 rx = (rx - 1.0)*sr_x + 1.0
2568 ry = (ry - 1.0)*sr_y + 1.0
2569 call select_domain(orig_selected_domain)
2570 if (real(start_x) <= rx .and. rx <= real(end_x) .and. &
2571 real(start_y) <= ry .and. ry <= real(end_y)) then
2572 where_maps_to(min_i,min_j,1) = nint(rx)
2573 where_maps_to(min_i,min_j,2) = nint(ry)
2574 else
2575 where_maps_to(min_i,min_j,1) = OUTSIDE_DOMAIN
2576 end if
2577 end if
2579 ! Upper-left corner
2580 if (where_maps_to(min_i,max_j,1) == NOT_PROCESSED) then
2581 orig_selected_domain = iget_selected_domain()
2582 call select_domain(SOURCE_PROJ)
2583 call xytoll(real(min_i), real(max_j), lat_corner, lon_corner, istagger)
2584 call select_domain(1)
2585 call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
2586 rx = (rx - 1.0)*sr_x + 1.0
2587 ry = (ry - 1.0)*sr_y + 1.0
2588 call select_domain(orig_selected_domain)
2589 if (real(start_x) <= rx .and. rx <= real(end_x) .and. &
2590 real(start_y) <= ry .and. ry <= real(end_y)) then
2591 where_maps_to(min_i,max_j,1) = nint(rx)
2592 where_maps_to(min_i,max_j,2) = nint(ry)
2593 else
2594 where_maps_to(min_i,max_j,1) = OUTSIDE_DOMAIN
2595 end if
2596 end if
2598 ! Upper-right corner
2599 if (where_maps_to(max_i,max_j,1) == NOT_PROCESSED) then
2600 orig_selected_domain = iget_selected_domain()
2601 call select_domain(SOURCE_PROJ)
2602 call xytoll(real(max_i), real(max_j), lat_corner, lon_corner, istagger)
2603 call select_domain(1)
2604 call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
2605 rx = (rx - 1.0)*sr_x + 1.0
2606 ry = (ry - 1.0)*sr_y + 1.0
2607 call select_domain(orig_selected_domain)
2608 if (real(start_x) <= rx .and. rx <= real(end_x) .and. &
2609 real(start_y) <= ry .and. ry <= real(end_y)) then
2610 where_maps_to(max_i,max_j,1) = nint(rx)
2611 where_maps_to(max_i,max_j,2) = nint(ry)
2612 else
2613 where_maps_to(max_i,max_j,1) = OUTSIDE_DOMAIN
2614 end if
2615 end if
2617 ! Lower-right corner
2618 if (where_maps_to(max_i,min_j,1) == NOT_PROCESSED) then
2619 orig_selected_domain = iget_selected_domain()
2620 call select_domain(SOURCE_PROJ)
2621 call xytoll(real(max_i), real(min_j), lat_corner, lon_corner, istagger)
2622 call select_domain(1)
2623 call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
2624 rx = (rx - 1.0)*sr_x + 1.0
2625 ry = (ry - 1.0)*sr_y + 1.0
2626 call select_domain(orig_selected_domain)
2627 if (real(start_x) <= rx .and. rx <= real(end_x) .and. &
2628 real(start_y) <= ry .and. ry <= real(end_y)) then
2629 where_maps_to(max_i,min_j,1) = nint(rx)
2630 where_maps_to(max_i,min_j,2) = nint(ry)
2631 else
2632 where_maps_to(max_i,min_j,1) = OUTSIDE_DOMAIN
2633 end if
2634 end if
2636 ! If all four corners map to same model grid point, accumulate the
2637 ! entire rectangle
2638 if (where_maps_to(min_i,min_j,1) == where_maps_to(min_i,max_j,1) .and. &
2639 where_maps_to(min_i,min_j,1) == where_maps_to(max_i,max_j,1) .and. &
2640 where_maps_to(min_i,min_j,1) == where_maps_to(max_i,min_j,1) .and. &
2641 where_maps_to(min_i,min_j,2) == where_maps_to(min_i,max_j,2) .and. &
2642 where_maps_to(min_i,min_j,2) == where_maps_to(max_i,max_j,2) .and. &
2643 where_maps_to(min_i,min_j,2) == where_maps_to(max_i,min_j,2) .and. &
2644 where_maps_to(min_i,min_j,1) /= OUTSIDE_DOMAIN) then
2645 x_dest = where_maps_to(min_i,min_j,1)
2646 y_dest = where_maps_to(min_i,min_j,2)
2648 ! If this grid point was already given a value from higher-priority source data,
2649 ! there is nothing to do.
2650 ! if (.not. bitarray_test(processed_pts, x_dest-start_x+1, y_dest-start_y+1)) then
2652 ! If this grid point has never been given a value by this level of source data,
2653 ! initialize the point
2654 if (.not. bitarray_test(new_pts, x_dest-start_x+1, y_dest-start_y+1)) then
2655 do k=min_k,max_k
2656 dst_array(x_dest,y_dest,k) = 0.
2657 end do
2658 end if
2660 ! Sum all the points whose nearest neighbor is this grid point
2661 if (present(mask_array) .and. present(mask_relational)) then
2662 do i=min_i,max_i
2663 do j=min_j,max_j
2664 do k=min_k,max_k
2665 ! Ignore masked/missing values in the source data
2666 if (tile_array(i,j,k) /= msgval) then
2667 if (mask_relational == ' ' .and. mask_array(i,j) /= maskval) then
2668 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2669 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2670 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2671 else if (mask_relational == '<' .and. mask_array(i,j) >= maskval) then
2672 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2673 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2674 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2675 else if (mask_relational == '>' .and. mask_array(i,j) <= maskval) then
2676 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2677 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2678 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2679 end if
2680 end if
2681 end do
2682 end do
2683 end do
2684 else if (present(mask_array)) then
2685 do i=min_i,max_i
2686 do j=min_j,max_j
2687 do k=min_k,max_k
2688 ! Ignore masked/missing values in the source data
2689 if ((tile_array(i,j,k) /= msgval) .and. &
2690 (mask_array(i,j) /= maskval)) then
2691 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2692 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2693 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2694 end if
2695 end do
2696 end do
2697 end do
2698 else
2699 do i=min_i,max_i
2700 do j=min_j,max_j
2701 do k=min_k,max_k
2702 ! Ignore masked/missing values in the source data
2703 if ((tile_array(i,j,k) /= msgval)) then
2704 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2705 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2706 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2707 end if
2708 end do
2709 end do
2710 end do
2711 end if
2713 ! end if
2715 ! Rectangle is a square of four points, and we can simply deal with each of the points
2716 else if (((max_i - min_i + 1) <= 2) .and. ((max_j - min_j + 1) <= 2)) then
2717 do i=min_i,max_i
2718 do j=min_j,max_j
2719 x_dest = where_maps_to(i,j,1)
2720 y_dest = where_maps_to(i,j,2)
2722 if (x_dest /= OUTSIDE_DOMAIN) then
2724 ! if (.not. bitarray_test(processed_pts, x_dest-start_x+1, y_dest-start_y+1)) then
2725 if (.not. bitarray_test(new_pts, x_dest-start_x+1, y_dest-start_y+1)) then
2726 do k=min_k,max_k
2727 dst_array(x_dest,y_dest,k) = 0.
2728 end do
2729 end if
2731 if (present(mask_array) .and. present(mask_relational)) then
2732 do k=min_k,max_k
2733 ! Ignore masked/missing values
2734 if (tile_array(i,j,k) /= msgval) then
2735 if (mask_relational == ' ' .and. mask_array(i,j) /= maskval) then
2736 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2737 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2738 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2739 else if (mask_relational == '<' .and. mask_array(i,j) >= maskval) then
2740 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2741 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2742 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2743 else if (mask_relational == '>' .and. mask_array(i,j) <= maskval) then
2744 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2745 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2746 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2747 end if
2748 end if
2749 end do
2750 else if (present(mask_array)) then
2751 do k=min_k,max_k
2752 ! Ignore masked/missing values
2753 if ((tile_array(i,j,k) /= msgval) .and. &
2754 (mask_array(i,j) /= maskval)) then
2755 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2756 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2757 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2758 end if
2759 end do
2760 else
2761 do k=min_k,max_k
2762 ! Ignore masked/missing values
2763 if ((tile_array(i,j,k) /= msgval)) then
2764 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2765 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2766 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2767 end if
2768 end do
2769 end if
2770 ! end if
2772 end if
2773 end do
2774 end do
2776 ! Not all corners map to the same grid point, and the rectangle contains more than
2777 ! four points
2778 else
2779 center_i = (max_i + min_i)/2
2780 center_j = (max_j + min_j)/2
2782 ! Recursively process lower-left rectangle
2783 call process_continuous_block(tile_array, where_maps_to, &
2784 src_min_x, src_min_y, src_min_z, &
2785 src_max_x, src_max_y, src_max_z, &
2786 min_i, min_j, min_k, &
2787 center_i, center_j, max_k, &
2788 dst_array, n, &
2789 start_x, end_x, start_y, end_y, start_z, end_z, &
2790 istagger, &
2791 new_pts, sr_x, sr_y, msgval, maskval, mask_relational, mask_array)
2793 if (center_i < max_i) then
2794 ! Recursively process lower-right rectangle
2795 call process_continuous_block(tile_array, where_maps_to, &
2796 src_min_x, src_min_y, src_min_z, &
2797 src_max_x, src_max_y, src_max_z, &
2798 center_i+1, min_j, min_k, max_i, &
2799 center_j, max_k, &
2800 dst_array, n, &
2801 start_x, end_x, start_y, &
2802 end_y, start_z, end_z, &
2803 istagger, &
2804 new_pts, sr_x, sr_y, msgval, maskval, mask_relational, mask_array)
2805 end if
2807 if (center_j < max_j) then
2808 ! Recursively process upper-left rectangle
2809 call process_continuous_block(tile_array, where_maps_to, &
2810 src_min_x, src_min_y, src_min_z, &
2811 src_max_x, src_max_y, src_max_z, &
2812 min_i, center_j+1, min_k, center_i, &
2813 max_j, max_k, &
2814 dst_array, n, &
2815 start_x, end_x, start_y, &
2816 end_y, start_z, end_z, &
2817 istagger, &
2818 new_pts, sr_x, sr_y, msgval, maskval, mask_relational, mask_array)
2819 end if
2821 if (center_i < max_i .and. center_j < max_j) then
2822 ! Recursively process upper-right rectangle
2823 call process_continuous_block(tile_array, where_maps_to, &
2824 src_min_x, src_min_y, src_min_z, &
2825 src_max_x, src_max_y, src_max_z, &
2826 center_i+1, center_j+1, min_k, max_i, &
2827 max_j, max_k, &
2828 dst_array, n, &
2829 start_x, end_x, start_y, &
2830 end_y, start_z, end_z, &
2831 istagger, &
2832 new_pts, sr_x, sr_y, msgval, maskval, mask_relational, mask_array)
2833 end if
2834 end if
2836 end subroutine process_continuous_block
2838 end module process_domain_module