| blob | c4da0a891b483af186526ad4882785ffc8c46bc1 |
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
46 character (len=128), dimension(:), pointer :: geogrid_flags
48 ! CWH Initialize local pointer variables
49 nullify(landmask)
50 nullify(xlat)
51 nullify(xlon)
52 nullify(xlat_u)
53 nullify(xlon_u)
54 nullify(xlat_v)
55 nullify(xlon_v)
56 nullify(geogrid_flags)
58 ! Compute number of times that we will process
59 call geth_idts(end_date(n), start_date(n), idiff)
60 call mprintf((idiff < 0),ERROR,'Ending date is earlier than starting date in namelist for domain %i.', i1=n)
62 n_times = idiff / interval_seconds
64 ! Check that the interval evenly divides the range of times to process
65 call mprintf((mod(idiff, interval_seconds) /= 0),WARN, &
66 'In namelist, interval_seconds does not evenly divide '// &
67 '(end_date - start_date) for domain %i. Only %i time periods '// &
68 'will be processed.', i1=n, i2=n_times)
70 ! Initialize the storage module
71 call mprintf(.true.,LOGFILE,'Initializing storage module')
72 call storage_init()
74 !
75 ! Do time-independent processing
76 !
77 call get_static_fields(n, dyn_opt, west_east_dim, south_north_dim, bottom_top_dim, map_proj, &
78 we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
79 we_patch_s, we_patch_e, &
80 we_patch_stag_s, we_patch_stag_e, &
81 sn_patch_s, sn_patch_e, &
82 sn_patch_stag_s, sn_patch_stag_e, &
83 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
84 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
85 mminlu, num_land_cat, is_water, is_lake, is_ice, is_urban, i_soilwater, &
86 grid_id, parent_id, &
87 i_parent_start, j_parent_start, i_parent_end, j_parent_end, &
88 parent_grid_ratio, sub_x, sub_y, &
89 cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
90 pole_lat, pole_lon, dom_dx, dom_dy, landmask, xlat, xlon, xlat_u, xlon_u, &
91 xlat_v, xlon_v, corner_lats, corner_lons, title, geogrid_flags)
94 allocate(output_flags(num_entries))
95 allocate(got_const_field(num_entries))
97 do i=1,num_entries
98 output_flags(i) = ' '
99 got_const_field(i) = .false.
100 end do
102 ! This call is to process the constant met fields (SST or SEAICE, for example)
103 ! That we process constant fields is indicated by the first argument
104 call process_single_met_time(.true., temp_date, n, extra_row, extra_col, xlat, xlon, &
105 xlat_u, xlon_u, xlat_v, xlon_v, landmask, &
106 title, dyn_opt, &
107 west_east_dim, south_north_dim, &
108 we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
109 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
110 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
111 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
112 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
113 got_const_field, &
114 map_proj, mminlu, num_land_cat, is_water, is_lake, is_ice, is_urban, i_soilwater, &
115 grid_id, parent_id, i_parent_start, &
116 j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
117 cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
118 pole_lat, pole_lon, parent_grid_ratio, sub_x, sub_y, &
119 corner_lats, corner_lons, output_flags, geogrid_flags, 0)
121 !
122 ! Begin time-dependent processing
123 !
125 allocate(td_output_flags(num_entries))
126 allocate(got_this_field (num_entries))
128 ! Loop over all times to be processed for this domain
129 do t=0,n_times
131 call geth_newdate(valid_date, trim(start_date(n)), t*interval_seconds)
132 temp_date = ' '
134 if (mod(interval_seconds,3600) == 0) then
135 write(temp_date,'(a13)') valid_date(1:10)//'_'//valid_date(12:13)
136 else if (mod(interval_seconds,60) == 0) then
137 write(temp_date,'(a16)') valid_date(1:10)//'_'//valid_date(12:16)
138 else
139 write(temp_date,'(a19)') valid_date(1:10)//'_'//valid_date(12:19)
140 end if
142 call mprintf(.true.,STDOUT, ' Processing %s', s1=trim(temp_date))
143 call mprintf(.true.,LOGFILE, 'Preparing to process output time %s', s1=temp_date)
145 do i=1,num_entries
146 td_output_flags(i) = output_flags(i)
147 got_this_field(i) = got_const_field(i)
148 end do
150 if (t > 0) then
151 process_bdy_width = process_only_bdy
152 else
153 process_bdy_width = 0
154 end if
156 call process_single_met_time(.false., temp_date, n, extra_row, extra_col, xlat, xlon, &
157 xlat_u, xlon_u, xlat_v, xlon_v, landmask, &
158 title, dyn_opt, &
159 west_east_dim, south_north_dim, &
160 we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
161 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
162 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
163 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
164 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
165 got_this_field, &
166 map_proj, mminlu, num_land_cat, is_water, is_lake, is_ice, is_urban, i_soilwater, &
167 grid_id, parent_id, i_parent_start, &
168 j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
169 cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
170 pole_lat, pole_lon, parent_grid_ratio, sub_x, sub_y, &
171 corner_lats, corner_lons, td_output_flags, geogrid_flags, process_bdy_width)
173 end do ! Loop over n_times
176 deallocate(td_output_flags)
177 deallocate(got_this_field)
179 deallocate(output_flags)
180 deallocate(got_const_field)
182 if (associated(geogrid_flags)) deallocate(geogrid_flags)
184 call storage_delete_all()
186 end subroutine process_domain
189 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
190 ! Name: get_static_fields
191 !
192 ! Purpose:
193 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
194 subroutine get_static_fields(n, dyn_opt, west_east_dim, south_north_dim, bottom_top_dim, &
195 map_proj, &
196 we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
197 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
198 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
199 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
200 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
201 mminlu, num_land_cat, is_water, is_lake, is_ice, is_urban, i_soilwater, &
202 grid_id, parent_id, &
203 i_parent_start, j_parent_start, i_parent_end, j_parent_end, &
204 parent_grid_ratio, sub_x, sub_y, &
205 cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
206 pole_lat, pole_lon, dom_dx, dom_dy, landmask, xlat, xlon, xlat_u, xlon_u, &
207 xlat_v, xlon_v, corner_lats, corner_lons, title, geogrid_flags)
209 use gridinfo_module
210 use input_module
211 use llxy_module
212 use parallel_module
213 use storage_module
214 use module_debug
215 use list_module
217 implicit none
219 ! Arguments
220 integer, intent(in) :: n
221 integer, intent(inout) :: dyn_opt, west_east_dim, south_north_dim, bottom_top_dim, &
222 map_proj, &
223 we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
224 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
225 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
226 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
227 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
228 is_water, is_lake, is_ice, is_urban, i_soilwater, grid_id, parent_id, &
229 i_parent_start, j_parent_start, i_parent_end, j_parent_end, &
230 parent_grid_ratio, sub_x, sub_y, num_land_cat
231 real, pointer, dimension(:,:) :: landmask
232 real, intent(inout) :: cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
233 dom_dx, dom_dy, pole_lat, pole_lon
234 real, pointer, dimension(:,:) :: xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v
235 real, dimension(16), intent(out) :: corner_lats, corner_lons
236 character (len=128), intent(inout) :: title, mminlu
237 character (len=128), dimension(:), pointer :: geogrid_flags
239 ! Local variables
240 integer :: istatus, i, j, k, sp1, ep1, sp2, ep2, sp3, ep3, &
241 lh_mult, rh_mult, bh_mult, th_mult, subx, suby
242 integer :: we_mem_subgrid_s, we_mem_subgrid_e, &
243 sn_mem_subgrid_s, sn_mem_subgrid_e
244 integer :: we_patch_subgrid_s, we_patch_subgrid_e, &
245 sn_patch_subgrid_s, sn_patch_subgrid_e
246 real, pointer, dimension(:,:,:) :: real_array
247 character (len=3) :: memorder
248 character (len=128) :: grid_type, datestr, cname, stagger, cunits, cdesc, ctemp
249 character (len=128), dimension(3) :: dimnames
250 type (fg_input) :: field
251 type (list) :: static_list, flag_list
252 type (list_item), dimension(:), pointer :: static_list_array
254 call list_init(static_list)
256 ! CWH Initialize local pointer variables
257 nullify(real_array)
259 ! Initialize the input module to read static input data for this domain
260 call mprintf(.true.,LOGFILE,'Opening static input file.')
261 call input_init(n, istatus)
262 call mprintf((istatus /= 0),ERROR, 'input_init(): Error opening input for domain %i.', i1=n)
264 ! Read global attributes from the static data input file
265 call mprintf(.true.,LOGFILE,'Reading static global attributes.')
266 call read_global_attrs(title, datestr, grid_type, dyn_opt, west_east_dim, &
267 south_north_dim, bottom_top_dim, &
268 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
269 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
270 map_proj, mminlu, num_land_cat, &
271 is_water, is_lake, is_ice, is_urban, i_soilwater, &
272 grid_id, parent_id, i_parent_start, &
273 j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
274 cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, &
275 truelat2, pole_lat, pole_lon, parent_grid_ratio, &
276 corner_lats, corner_lons, sub_x, sub_y)
278 we_dom_s = 1
279 sn_dom_s = 1
280 if (grid_type(1:1) == 'C') then
281 we_dom_e = west_east_dim - 1
282 sn_dom_e = south_north_dim - 1
283 else if (grid_type(1:1) == 'E') then
284 we_dom_e = west_east_dim
285 sn_dom_e = south_north_dim
286 end if
288 ! Given the full dimensions of this domain, find out the range of indices
289 ! that will be worked on by this processor. This information is given by
290 ! my_minx, my_miny, my_maxx, my_maxy
291 call parallel_get_tile_dims(west_east_dim, south_north_dim)
293 ! Must figure out patch dimensions from info in parallel module
294 if (nprocs > 1 .and. .not. do_tiled_input) then
296 we_patch_s = my_minx
297 we_patch_stag_s = my_minx
298 we_patch_e = my_maxx - 1
299 sn_patch_s = my_miny
300 sn_patch_stag_s = my_miny
301 sn_patch_e = my_maxy - 1
303 if (gridtype == 'C') then
304 if (my_x /= nproc_x - 1) then
305 we_patch_e = we_patch_e + 1
306 we_patch_stag_e = we_patch_e
307 else
308 we_patch_stag_e = we_patch_e + 1
309 end if
310 if (my_y /= nproc_y - 1) then
311 sn_patch_e = sn_patch_e + 1
312 sn_patch_stag_e = sn_patch_e
313 else
314 sn_patch_stag_e = sn_patch_e + 1
315 end if
316 else if (gridtype == 'E') then
317 we_patch_e = we_patch_e + 1
318 sn_patch_e = sn_patch_e + 1
319 we_patch_stag_e = we_patch_e
320 sn_patch_stag_e = sn_patch_e
321 end if
323 end if
325 ! Compute multipliers for halo width; these must be 0/1
326 if (my_x /= 0) then
327 lh_mult = 1
328 else
329 lh_mult = 0
330 end if
331 if (my_x /= (nproc_x-1)) then
332 rh_mult = 1
333 else
334 rh_mult = 0
335 end if
336 if (my_y /= 0) then
337 bh_mult = 1
338 else
339 bh_mult = 0
340 end if
341 if (my_y /= (nproc_y-1)) then
342 th_mult = 1
343 else
344 th_mult = 0
345 end if
347 we_mem_s = we_patch_s - HALO_WIDTH*lh_mult
348 we_mem_e = we_patch_e + HALO_WIDTH*rh_mult
349 sn_mem_s = sn_patch_s - HALO_WIDTH*bh_mult
350 sn_mem_e = sn_patch_e + HALO_WIDTH*th_mult
351 we_mem_stag_s = we_patch_stag_s - HALO_WIDTH*lh_mult
352 we_mem_stag_e = we_patch_stag_e + HALO_WIDTH*rh_mult
353 sn_mem_stag_s = sn_patch_stag_s - HALO_WIDTH*bh_mult
354 sn_mem_stag_e = sn_patch_stag_e + HALO_WIDTH*th_mult
356 ! Initialize a proj_info type for the destination grid projection. This will
357 ! primarily be used for rotating Earth-relative winds to grid-relative winds
358 call set_domain_projection(map_proj, stand_lon, truelat1, truelat2, &
359 dom_dx, dom_dy, dom_dx, dom_dy, west_east_dim, &
360 south_north_dim, real(west_east_dim)/2., &
361 real(south_north_dim)/2.,cen_lat, cen_lon, &
362 cen_lat, cen_lon)
364 ! Read static fields using the input module; we know that there are no more
365 ! fields to be read when read_next_field() returns a non-zero status.
366 istatus = 0
367 do while (istatus == 0)
368 call read_next_field(sp1, ep1, sp2, ep2, sp3, ep3, cname, cunits, cdesc, &
369 memorder, stagger, dimnames, subx, suby, &
370 real_array, istatus)
371 if (istatus == 0) then
373 call mprintf(.true.,LOGFILE, 'Read in static field %s.',s1=cname)
374 call list_insert(static_list, ckey=cname, cvalue=cname)
376 ! We will also keep copies in core of the lat/lon arrays, for use in
377 ! interpolation of the met fields to the model grid.
378 ! For now, we assume that the lat/lon arrays will have known field names
379 if (index(cname, 'XLAT_M') /= 0 .and. &
380 len_trim(cname) == len_trim('XLAT_M')) then
381 allocate(xlat(we_mem_s:we_mem_e,sn_mem_s:sn_mem_e))
382 xlat(we_patch_s:we_patch_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
383 call exchange_halo_r(xlat, &
384 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, 1, 1, &
385 we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, 1, 1)
387 else if (index(cname, 'XLONG_M') /= 0 .and. &
388 len_trim(cname) == len_trim('XLONG_M')) then
389 allocate(xlon(we_mem_s:we_mem_e,sn_mem_s:sn_mem_e))
390 xlon(we_patch_s:we_patch_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
391 call exchange_halo_r(xlon, &
392 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, 1, 1, &
393 we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, 1, 1)
395 else if (index(cname, 'XLAT_U') /= 0 .and. &
396 len_trim(cname) == len_trim('XLAT_U')) then
397 allocate(xlat_u(we_mem_stag_s:we_mem_stag_e,sn_mem_s:sn_mem_e))
398 xlat_u(we_patch_stag_s:we_patch_stag_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
399 call exchange_halo_r(xlat_u, &
400 we_mem_stag_s, we_mem_stag_e, sn_mem_s, sn_mem_e, 1, 1, &
401 we_patch_stag_s, we_patch_stag_e, sn_patch_s, sn_patch_e, 1, 1)
403 else if (index(cname, 'XLONG_U') /= 0 .and. &
404 len_trim(cname) == len_trim('XLONG_U')) then
405 allocate(xlon_u(we_mem_stag_s:we_mem_stag_e,sn_mem_s:sn_mem_e))
406 xlon_u(we_patch_stag_s:we_patch_stag_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
407 call exchange_halo_r(xlon_u, &
408 we_mem_stag_s, we_mem_stag_e, sn_mem_s, sn_mem_e, 1, 1, &
409 we_patch_stag_s, we_patch_stag_e, sn_patch_s, sn_patch_e, 1, 1)
411 else if (index(cname, 'XLAT_V') /= 0 .and. &
412 len_trim(cname) == len_trim('XLAT_V')) then
413 allocate(xlat_v(we_mem_s:we_mem_e,sn_mem_stag_s:sn_mem_stag_e))
414 xlat_v(we_patch_s:we_patch_e,sn_patch_stag_s:sn_patch_stag_e) = real_array(:,:,1)
415 call exchange_halo_r(xlat_v, &
416 we_mem_s, we_mem_e, sn_mem_stag_s, sn_mem_stag_e, 1, 1, &
417 we_patch_s, we_patch_e, sn_patch_stag_s, sn_patch_stag_e, 1, 1)
419 else if (index(cname, 'XLONG_V') /= 0 .and. &
420 len_trim(cname) == len_trim('XLONG_V')) then
421 allocate(xlon_v(we_mem_s:we_mem_e,sn_mem_stag_s:sn_mem_stag_e))
422 xlon_v(we_patch_s:we_patch_e,sn_patch_stag_s:sn_patch_stag_e) = real_array(:,:,1)
423 call exchange_halo_r(xlon_v, &
424 we_mem_s, we_mem_e, sn_mem_stag_s, sn_mem_stag_e, 1, 1, &
425 we_patch_s, we_patch_e, sn_patch_stag_s, sn_patch_stag_e, 1, 1)
427 else if (index(cname, 'LANDMASK') /= 0 .and. &
428 len_trim(cname) == len_trim('LANDMASK')) then
429 allocate(landmask(we_mem_s:we_mem_e,sn_mem_s:sn_mem_e))
430 landmask(we_patch_s:we_patch_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
431 call exchange_halo_r(landmask, &
432 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, 1, 1, &
433 we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, 1, 1)
435 end if
437 if (subx > 1) then
438 we_mem_subgrid_s = (we_mem_s + HALO_WIDTH*lh_mult - 1)*subx - HALO_WIDTH*lh_mult + 1
439 we_mem_subgrid_e = (we_mem_e + (1-rh_mult) - HALO_WIDTH*rh_mult )*subx + HALO_WIDTH*rh_mult
440 we_patch_subgrid_s = (we_patch_s - 1)*subx + 1
441 we_patch_subgrid_e = (we_patch_e + (1-rh_mult) )*subx
442 end if
443 if (suby > 1) then
444 sn_mem_subgrid_s = (sn_mem_s + HALO_WIDTH*bh_mult - 1)*suby - HALO_WIDTH*bh_mult + 1
445 sn_mem_subgrid_e = (sn_mem_e + (1-th_mult) - HALO_WIDTH*th_mult )*suby + HALO_WIDTH*th_mult
446 sn_patch_subgrid_s = (sn_patch_s - 1)*suby + 1
447 sn_patch_subgrid_e = (sn_patch_e + (1-th_mult) )*suby
448 end if
450 ! Having read in a field, we write each level individually to the
451 ! storage module; levels will be reassembled later on when they
452 ! are written.
453 do k=sp3,ep3
454 field%header%sr_x=subx
455 field%header%sr_y=suby
456 field%header%version = 1
457 field%header%date = start_date(n)
458 field%header%time_dependent = .false.
459 field%header%mask_field = .false.
460 field%header%constant_field = .false.
461 field%header%forecast_hour = 0.0
462 field%header%fg_source = 'geogrid_model'
463 field%header%field = cname
464 field%header%units = cunits
465 field%header%description = cdesc
466 field%header%vertical_coord = dimnames(3)
467 field%header%vertical_level = k
468 field%header%array_order = memorder
469 field%header%is_wind_grid_rel = .true.
470 field%header%array_has_missing_values = .false.
471 if (gridtype == 'C') then
472 if (subx > 1 .or. suby > 1) then
473 field%map%stagger = M
474 field%header%dim1(1) = we_mem_subgrid_s
475 field%header%dim1(2) = we_mem_subgrid_e
476 field%header%dim2(1) = sn_mem_subgrid_s
477 field%header%dim2(2) = sn_mem_subgrid_e
478 else if (trim(stagger) == 'M') then
479 field%map%stagger = M
480 field%header%dim1(1) = we_mem_s
481 field%header%dim1(2) = we_mem_e
482 field%header%dim2(1) = sn_mem_s
483 field%header%dim2(2) = sn_mem_e
484 else if (trim(stagger) == 'U') then
485 field%map%stagger = U
486 field%header%dim1(1) = we_mem_stag_s
487 field%header%dim1(2) = we_mem_stag_e
488 field%header%dim2(1) = sn_mem_s
489 field%header%dim2(2) = sn_mem_e
490 else if (trim(stagger) == 'V') then
491 field%map%stagger = V
492 field%header%dim1(1) = we_mem_s
493 field%header%dim1(2) = we_mem_e
494 field%header%dim2(1) = sn_mem_stag_s
495 field%header%dim2(2) = sn_mem_stag_e
496 else if (trim(stagger) == 'CORNER') then
497 field%map%stagger = CORNER
498 field%header%dim1(1) = we_mem_stag_s
499 field%header%dim1(2) = we_mem_stag_e
500 field%header%dim2(1) = sn_mem_stag_s
501 field%header%dim2(2) = sn_mem_stag_e
502 end if
503 else if (gridtype == 'E') then
504 if (trim(stagger) == 'M') then
505 field%map%stagger = HH
506 else if (trim(stagger) == 'V') then
507 field%map%stagger = VV
508 end if
509 field%header%dim1(1) = we_mem_s
510 field%header%dim1(2) = we_mem_e
511 field%header%dim2(1) = sn_mem_s
512 field%header%dim2(2) = sn_mem_e
513 end if
515 allocate(field%valid_mask)
517 if (subx > 1 .or. suby > 1) then
518 allocate(field%r_arr(we_mem_subgrid_s:we_mem_subgrid_e,&
519 sn_mem_subgrid_s:sn_mem_subgrid_e))
520 field%r_arr(we_patch_subgrid_s:we_patch_subgrid_e,sn_patch_subgrid_s:sn_patch_subgrid_e) = &
521 real_array(sp1:ep1,sp2:ep2,k)
522 call exchange_halo_r(field%r_arr, &
523 we_mem_subgrid_s, we_mem_subgrid_e, sn_mem_subgrid_s, sn_mem_subgrid_e, 1, 1, &
524 we_patch_subgrid_s, we_patch_subgrid_e, sn_patch_subgrid_s, sn_patch_subgrid_e, 1, 1)
525 call bitarray_create(field%valid_mask, &
526 (we_mem_subgrid_e-we_mem_subgrid_s)+1, &
527 (sn_mem_subgrid_e-sn_mem_subgrid_s)+1)
528 do j=1,(sn_mem_subgrid_e-sn_mem_subgrid_s)+1
529 do i=1,(we_mem_subgrid_e-we_mem_subgrid_s)+1
530 call bitarray_set(field%valid_mask, i, j)
531 end do
532 end do
534 else if (field%map%stagger == M .or. &
535 field%map%stagger == HH .or. &
536 field%map%stagger == VV) then
537 allocate(field%r_arr(we_mem_s:we_mem_e,&
538 sn_mem_s:sn_mem_e))
539 field%r_arr(we_patch_s:we_patch_e,sn_patch_s:sn_patch_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_s, sn_mem_e, 1, 1, &
542 we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, 1, 1)
543 call bitarray_create(field%valid_mask, &
544 (we_mem_e-we_mem_s)+1, &
545 (sn_mem_e-sn_mem_s)+1)
546 do j=1,(sn_mem_e-sn_mem_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 else if (field%map%stagger == U) then
552 allocate(field%r_arr(we_mem_stag_s:we_mem_stag_e,&
553 sn_mem_s:sn_mem_e))
554 field%r_arr(we_patch_stag_s:we_patch_stag_e,sn_patch_s:sn_patch_e) = real_array(sp1:ep1,sp2:ep2,k)
555 call exchange_halo_r(field%r_arr, &
556 we_mem_stag_s, we_mem_stag_e, sn_mem_s, sn_mem_e, 1, 1, &
557 we_patch_stag_s, we_patch_stag_e, sn_patch_s, sn_patch_e, 1, 1)
558 call bitarray_create(field%valid_mask, &
559 (we_mem_stag_e-we_mem_stag_s)+1, &
560 (sn_mem_e-sn_mem_s)+1)
561 do j=1,(sn_mem_e-sn_mem_s)+1
562 do i=1,(we_mem_stag_e-we_mem_stag_s)+1
563 call bitarray_set(field%valid_mask, i, j)
564 end do
565 end do
566 else if (field%map%stagger == V) then
567 allocate(field%r_arr(we_mem_s:we_mem_e,&
568 sn_mem_stag_s:sn_mem_stag_e))
569 field%r_arr(we_patch_s:we_patch_e,sn_patch_stag_s:sn_patch_stag_e) = real_array(sp1:ep1,sp2:ep2,k)
570 call exchange_halo_r(field%r_arr, &
571 we_mem_s, we_mem_e, sn_mem_stag_s, sn_mem_stag_e, 1, 1, &
572 we_patch_s, we_patch_e, sn_patch_stag_s, sn_patch_stag_e, 1, 1)
573 call bitarray_create(field%valid_mask, &
574 (we_mem_e-we_mem_s)+1, &
575 (sn_mem_stag_e-sn_mem_stag_s)+1)
576 do j=1,(sn_mem_stag_e-sn_mem_stag_s)+1
577 do i=1,(we_mem_e-we_mem_s)+1
578 call bitarray_set(field%valid_mask, i, j)
579 end do
580 end do
581 else if (field%map%stagger == CORNER) then
582 allocate(field%r_arr(we_mem_stag_s:we_mem_stag_e,&
583 sn_mem_stag_s:sn_mem_stag_e))
584 field%r_arr(we_patch_stag_s:we_patch_stag_e,sn_patch_stag_s:sn_patch_stag_e) = real_array(sp1:ep1,sp2:ep2,k)
585 call exchange_halo_r(field%r_arr, &
586 we_mem_stag_s, we_mem_stag_e, sn_mem_stag_s, sn_mem_stag_e, 1, 1, &
587 we_patch_stag_s, we_patch_stag_e, sn_patch_stag_s, sn_patch_stag_e, 1, 1)
588 call bitarray_create(field%valid_mask, &
589 (we_mem_stag_e-we_mem_stag_s)+1, &
590 (sn_mem_stag_e-sn_mem_stag_s)+1)
591 do j=1,(sn_mem_stag_e-sn_mem_stag_s)+1
592 do i=1,(we_mem_stag_e-we_mem_stag_s)+1
593 call bitarray_set(field%valid_mask, i, j)
594 end do
595 end do
596 end if
598 nullify(field%modified_mask)
600 call storage_put_field(field)
602 end do
604 end if
605 end do
607 static_list_array => list_get_keys(static_list)
608 call list_init(flag_list)
609 do i=1,size(static_list_array)
610 istatus = 0
611 ctemp = 'FLAG_'//trim(static_list_array(i)%ckey)
612 call ext_get_dom_ti_integer_scalar(ctemp, istatus, suppress_errors=.true.)
613 if (istatus == 1) call list_insert(flag_list, ckey=ctemp, cvalue=ctemp)
614 end do
615 deallocate(static_list_array)
616 call list_destroy(static_list)
618 ! Done reading all static fields for this domain
619 call input_close()
621 static_list_array => list_get_keys(flag_list)
622 allocate(geogrid_flags(size(static_list_array)))
623 do i=1,size(static_list_array)
624 geogrid_flags(i) = static_list_array(i)%ckey
625 end do
626 deallocate(static_list_array)
628 call list_destroy(flag_list)
630 end subroutine get_static_fields
633 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
634 ! Name: process_single_met_time
635 !
636 ! Purpose:
637 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
638 subroutine process_single_met_time(do_const_processing, &
639 temp_date, n, extra_row, extra_col, xlat, xlon, &
640 xlat_u, xlon_u, xlat_v, xlon_v, landmask, &
641 title, dyn_opt, &
642 west_east_dim, south_north_dim, &
643 we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
644 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
645 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
646 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
647 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
648 got_this_field, &
649 map_proj, mminlu, num_land_cat, is_water, is_lake, is_ice, is_urban, i_soilwater, &
650 grid_id, parent_id, i_parent_start, &
651 j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
652 cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
653 pole_lat, pole_lon, parent_grid_ratio, sub_x, sub_y, &
654 corner_lats, corner_lons, output_flags, geogrid_flags, process_bdy_width)
656 use bitarray_module
657 use gridinfo_module
658 use interp_module
659 use interp_option_module
660 use llxy_module
661 use misc_definitions_module
662 use module_debug
663 use output_module
664 use parallel_module
665 use read_met_module
666 use rotate_winds_module
667 use storage_module
669 implicit none
671 ! Arguments
672 logical, intent(in) :: do_const_processing
673 integer, intent(in) :: n, dyn_opt, west_east_dim, south_north_dim, map_proj, &
674 we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
675 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
676 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
677 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
678 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
679 is_water, is_lake, is_ice, is_urban, i_soilwater, &
680 grid_id, parent_id, i_parent_start, j_parent_start, &
681 i_parent_end, j_parent_end, parent_grid_ratio, sub_x, sub_y, num_land_cat, &
682 process_bdy_width
683 ! BUG: Should we be passing these around as pointers, or just declare them as arrays?
684 real, pointer, dimension(:,:) :: landmask
685 real, intent(in) :: dom_dx, dom_dy, cen_lat, moad_cen_lat, cen_lon, stand_lon, &
686 truelat1, truelat2, pole_lat, pole_lon
687 real, dimension(16), intent(in) :: corner_lats, corner_lons
688 real, pointer, dimension(:,:) :: xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v
689 logical, intent(in) :: extra_row, extra_col
690 logical, dimension(:), intent(inout) :: got_this_field
691 character (len=19), intent(in) :: temp_date
692 character (len=128), intent(in) :: mminlu
693 character (len=128), dimension(:), intent(inout) :: output_flags
694 character (len=128), dimension(:), pointer :: geogrid_flags
696 ! BUG: Move this constant to misc_definitions_module?
697 integer, parameter :: BDR_WIDTH = 3
699 ! Local variables
700 integer :: istatus, iqstatus, fg_idx, idx, idxt, i, j, bottom_top_dim, &
701 sm1, em1, sm2, em2, sm3, em3, &
702 sp1, ep1, sp2, ep2, sp3, ep3, &
703 sd1, ed1, sd2, ed2, sd3, ed3, &
704 u_idx, bdr_wdth
705 integer :: nmet_flags
706 integer :: num_metgrid_soil_levs
707 integer, pointer, dimension(:) :: soil_levels
708 real :: rx, ry
709 real :: threshold
710 logical :: do_gcell_interp
711 integer, pointer, dimension(:) :: u_levels, v_levels
712 real, pointer, dimension(:,:) :: halo_slab
713 real, pointer, dimension(:,:,:) :: real_array
714 character (len=19) :: output_date
715 character (len=128) :: cname, title
716 character (len=MAX_FILENAME_LEN) :: input_name
717 character (len=128), allocatable, dimension(:) :: met_flags
718 type (fg_input) :: field, u_field, v_field
719 type (met_data) :: fg_data
721 ! CWH Initialize local pointer variables
722 nullify(soil_levels)
723 nullify(u_levels)
724 nullify(v_levels)
725 nullify(halo_slab)
726 nullify(real_array)
729 ! For this time, we need to process all first-guess filename roots. When we
730 ! hit a root containing a '*', we assume we have hit the end of the list
731 fg_idx = 1
732 if (do_const_processing) then
733 input_name = constants_name(fg_idx)
734 else
735 input_name = fg_name(fg_idx)
736 end if
737 do while (input_name /= '*')
739 call mprintf(.true.,STDOUT, ' %s', s1=input_name)
740 call mprintf(.true.,LOGFILE, 'Getting input fields from %s', s1=input_name)
742 ! Do a first pass through this fg source to get all mask fields used
743 ! during interpolation
744 call get_interp_masks(trim(input_name), do_const_processing, temp_date)
746 istatus = 0
748 ! Initialize the module for reading in the met fields
749 call read_met_init(trim(input_name), do_const_processing, temp_date, istatus)
751 if (istatus == 0) then
753 ! Process all fields and levels from the current file; read_next_met_field()
754 ! will return a non-zero status when there are no more fields to be read.
755 do while (istatus == 0)
757 call read_next_met_field(fg_data, istatus)
759 if (istatus == 0) then
761 ! Find index into fieldname, interp_method, masked, and fill_missing
762 ! of the current field
763 idxt = num_entries + 1
764 do idx=1,num_entries
765 if ((index(fieldname(idx), trim(fg_data%field)) /= 0) .and. &
766 (len_trim(fieldname(idx)) == len_trim(fg_data%field))) then
768 got_this_field(idx) = .true.
770 if (index(input_name,trim(from_input(idx))) /= 0 .or. &
771 (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
772 idxt = idx
773 end if
775 end if
776 end do
777 idx = idxt
778 if (idx > num_entries) idx = num_entries ! The last entry is a default
780 ! Do we need to rename this field?
781 if (output_name(idx) /= ' ') then
782 fg_data%field = output_name(idx)(1:9)
784 idxt = num_entries + 1
785 do idx=1,num_entries
786 if ((index(fieldname(idx), trim(fg_data%field)) /= 0) .and. &
787 (len_trim(fieldname(idx)) == len_trim(fg_data%field))) then
789 got_this_field(idx) = .true.
791 if (index(input_name,trim(from_input(idx))) /= 0 .or. &
792 (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
793 idxt = idx
794 end if
796 end if
797 end do
798 idx = idxt
799 if (idx > num_entries) idx = num_entries ! The last entry is a default
800 end if
802 ! Do a simple check to see whether this is a global dataset
803 ! Note that we do not currently support regional Gaussian grids
804 if ((fg_data%iproj == PROJ_LATLON .and. abs(fg_data%nx * fg_data%deltalon - 360.) < 0.0001) &
805 .or. (fg_data%iproj == PROJ_GAUSS)) then
806 bdr_wdth = BDR_WIDTH
807 allocate(halo_slab(1-BDR_WIDTH:fg_data%nx+BDR_WIDTH,1:fg_data%ny))
809 halo_slab(1:fg_data%nx, 1:fg_data%ny) = &
810 fg_data%slab(1:fg_data%nx, 1:fg_data%ny)
812 halo_slab(1-BDR_WIDTH:0, 1:fg_data%ny) = &
813 fg_data%slab(fg_data%nx-BDR_WIDTH+1:fg_data%nx, 1:fg_data%ny)
815 halo_slab(fg_data%nx+1:fg_data%nx+BDR_WIDTH, 1:fg_data%ny) = &
816 fg_data%slab(1:BDR_WIDTH, 1:fg_data%ny)
818 deallocate(fg_data%slab)
819 else
820 bdr_wdth = 0
821 halo_slab => fg_data%slab
822 nullify(fg_data%slab)
823 end if
825 call mprintf(.true.,LOGFILE,'Processing %s at level %f.',s1=fg_data%field,f1=fg_data%xlvl)
827 call push_source_projection(fg_data%iproj, fg_data%xlonc, fg_data%truelat1, &
828 fg_data%truelat2, fg_data%dx, fg_data%dy, fg_data%deltalat, &
829 fg_data%deltalon, fg_data%starti, fg_data%startj, &
830 fg_data%startlat, fg_data%startlon, earth_radius=fg_data%earth_radius*1000.)
832 ! Initialize fg_input structure to store the field
833 field%header%version = 1
834 field%header%date = fg_data%hdate//' '
835 if (do_const_processing) then
836 field%header%time_dependent = .false.
837 field%header%constant_field = .true.
838 else
839 field%header%time_dependent = .true.
840 field%header%constant_field = .false.
841 end if
842 field%header%forecast_hour = fg_data%xfcst
843 field%header%fg_source = 'FG'
844 field%header%field = ' '
845 field%header%field(1:9) = fg_data%field
846 field%header%units = ' '
847 field%header%units(1:25) = fg_data%units
848 field%header%description = ' '
849 field%header%description(1:46) = fg_data%desc
850 call get_z_dim_name(fg_data%field,field%header%vertical_coord)
851 field%header%vertical_level = nint(fg_data%xlvl)
852 field%header%sr_x = 1
853 field%header%sr_y = 1
854 field%header%array_order = 'XY '
855 field%header%is_wind_grid_rel = fg_data%is_wind_grid_rel
856 field%header%array_has_missing_values = .false.
857 nullify(field%r_arr)
858 nullify(field%valid_mask)
859 nullify(field%modified_mask)
861 if (output_this_field(idx) .and. flag_in_output(idx) /= ' ') then
862 output_flags(idx) = flag_in_output(idx)
863 end if
865 ! If we should not output this field, just list it as a mask field
866 if (output_this_field(idx)) then
867 field%header%mask_field = .false.
868 else
869 field%header%mask_field = .true.
870 end if
872 !
873 ! Before actually doing any interpolation to the model grid, we must check
874 ! whether we will be using the average_gcell interpolator that averages all
875 ! source points in each model grid cell
876 !
877 do_gcell_interp = .false.
878 if (index(interp_method(idx),'average_gcell') /= 0) then
880 call get_gcell_threshold(interp_method(idx), threshold, istatus)
881 if (istatus == 0) then
882 if (fg_data%dx == 0. .and. fg_data%dy == 0. .and. &
883 fg_data%deltalat /= 0. .and. fg_data%deltalon /= 0.) then
884 fg_data%dx = abs(fg_data%deltalon)
885 fg_data%dy = abs(fg_data%deltalat)
886 else
887 ! BUG: Need to more correctly handle dx/dy in meters.
888 fg_data%dx = fg_data%dx / 111000. ! Convert meters to approximate degrees
889 fg_data%dy = fg_data%dy / 111000.
890 end if
891 if (gridtype == 'C') then
892 if (threshold*max(fg_data%dx,fg_data%dy)*111. <= max(dom_dx,dom_dy)/1000.) &
893 do_gcell_interp = .true.
894 else if (gridtype == 'E') then
895 if (threshold*max(fg_data%dx,fg_data%dy) <= max(dom_dx,dom_dy)) &
896 do_gcell_interp = .true.
897 end if
898 end if
899 end if
901 ! Interpolate to U staggering
902 if (output_stagger(idx) == U) then
904 call storage_query_field(field, iqstatus)
905 if (iqstatus == 0) then
906 call storage_get_field(field, iqstatus)
907 call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// &
908 ' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl))
909 if (associated(field%modified_mask)) then
910 call bitarray_destroy(field%modified_mask)
911 nullify(field%modified_mask)
912 end if
913 else
914 allocate(field%valid_mask)
915 call bitarray_create(field%valid_mask, we_mem_stag_e-we_mem_stag_s+1, sn_mem_e-sn_mem_s+1)
916 end if
918 ! Save a copy of the fg_input structure for the U field so that we can find it later
919 if (is_u_field(idx)) call dup(field, u_field)
921 allocate(field%modified_mask)
922 call bitarray_create(field%modified_mask, we_mem_stag_e-we_mem_stag_s+1, sn_mem_e-sn_mem_s+1)
924 if (do_const_processing .or. field%header%time_dependent) then
925 call interp_met_field(input_name, fg_data%field, U, M, &
926 field, xlat_u, xlon_u, we_mem_stag_s, we_mem_stag_e, sn_mem_s, sn_mem_e, &
927 we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
928 halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
929 field%modified_mask, process_bdy_width)
930 else
931 call mprintf(.true.,INFORM,' - already processed this field from constant file.')
932 end if
934 ! Interpolate to V staggering
935 else if (output_stagger(idx) == V) then
937 call storage_query_field(field, iqstatus)
938 if (iqstatus == 0) then
939 call storage_get_field(field, iqstatus)
940 call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// &
941 ' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl))
942 if (associated(field%modified_mask)) then
943 call bitarray_destroy(field%modified_mask)
944 nullify(field%modified_mask)
945 end if
946 else
947 allocate(field%valid_mask)
948 call bitarray_create(field%valid_mask, we_mem_e-we_mem_s+1, sn_mem_stag_e-sn_mem_stag_s+1)
949 end if
951 ! Save a copy of the fg_input structure for the V field so that we can find it later
952 if (is_v_field(idx)) call dup(field, v_field)
954 allocate(field%modified_mask)
955 call bitarray_create(field%modified_mask, we_mem_e-we_mem_s+1, sn_mem_stag_e-sn_mem_stag_s+1)
957 if (do_const_processing .or. field%header%time_dependent) then
958 call interp_met_field(input_name, fg_data%field, V, M, &
959 field, xlat_v, xlon_v, we_mem_s, we_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
960 we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
961 halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
962 field%modified_mask, process_bdy_width)
963 else
964 call mprintf(.true.,INFORM,' - already processed this field from constant file.')
965 end if
967 ! Interpolate to VV staggering
968 else if (output_stagger(idx) == VV) then
970 call storage_query_field(field, iqstatus)
971 if (iqstatus == 0) then
972 call storage_get_field(field, iqstatus)
973 call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// &
974 ' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl))
975 if (associated(field%modified_mask)) then
976 call bitarray_destroy(field%modified_mask)
977 nullify(field%modified_mask)
978 end if
979 else
980 allocate(field%valid_mask)
981 call bitarray_create(field%valid_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1)
982 end if
984 ! Save a copy of the fg_input structure for the U field so that we can find it later
985 if (is_u_field(idx)) call dup(field, u_field)
987 ! Save a copy of the fg_input structure for the V field so that we can find it later
988 if (is_v_field(idx)) call dup(field, v_field)
990 allocate(field%modified_mask)
991 call bitarray_create(field%modified_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1)
993 if (do_const_processing .or. field%header%time_dependent) then
994 call interp_met_field(input_name, fg_data%field, VV, M, &
995 field, xlat_v, xlon_v, we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
996 we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
997 halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
998 field%modified_mask, process_bdy_width)
999 else
1000 call mprintf(.true.,INFORM,' - already processed this field from constant file.')
1001 end if
1003 ! All other fields interpolated to M staggering for C grid, H staggering for E grid
1004 else
1006 call storage_query_field(field, iqstatus)
1007 if (iqstatus == 0) then
1008 call storage_get_field(field, iqstatus)
1009 call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// &
1010 ' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl))
1011 if (associated(field%modified_mask)) then
1012 call bitarray_destroy(field%modified_mask)
1013 nullify(field%modified_mask)
1014 end if
1015 else
1016 allocate(field%valid_mask)
1017 call bitarray_create(field%valid_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1)
1018 end if
1020 allocate(field%modified_mask)
1021 call bitarray_create(field%modified_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1)
1023 if (do_const_processing .or. field%header%time_dependent) then
1024 if (gridtype == 'C') then
1025 call interp_met_field(input_name, fg_data%field, M, M, &
1026 field, xlat, xlon, we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
1027 we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
1028 halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
1029 field%modified_mask, process_bdy_width, landmask)
1031 else if (gridtype == 'E') then
1032 call interp_met_field(input_name, fg_data%field, HH, M, &
1033 field, xlat, xlon, we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
1034 we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
1035 halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
1036 field%modified_mask, process_bdy_width, landmask)
1037 end if
1038 else
1039 call mprintf(.true.,INFORM,' - already processed this field from constant file.')
1040 end if
1042 end if
1044 call bitarray_merge(field%valid_mask, field%modified_mask)
1046 deallocate(halo_slab)
1048 ! Store the interpolated field
1049 call storage_put_field(field)
1051 call pop_source_projection()
1053 end if
1054 end do
1056 call read_met_close()
1058 call push_source_projection(fg_data%iproj, fg_data%xlonc, fg_data%truelat1, &
1059 fg_data%truelat2, fg_data%dx, fg_data%dy, fg_data%deltalat, &
1060 fg_data%deltalon, fg_data%starti, fg_data%startj, &
1061 fg_data%startlat, fg_data%startlon, earth_radius=fg_data%earth_radius*1000.)
1063 !
1064 ! If necessary, rotate winds to earth-relative for this fg source
1065 !
1067 call storage_get_levels(u_field, u_levels)
1068 call storage_get_levels(v_field, v_levels)
1070 if (associated(u_levels) .and. associated(v_levels)) then
1071 u_idx = 1
1072 do u_idx = 1, size(u_levels)
1073 u_field%header%vertical_level = u_levels(u_idx)
1074 call storage_get_field(u_field, istatus)
1075 v_field%header%vertical_level = v_levels(u_idx)
1076 call storage_get_field(v_field, istatus)
1078 if (associated(u_field%modified_mask) .and. &
1079 associated(v_field%modified_mask)) then
1081 if (u_field%header%is_wind_grid_rel) then
1082 if (gridtype == 'C') then
1083 call map_to_met(u_field%r_arr, u_field%modified_mask, &
1084 v_field%r_arr, v_field%modified_mask, &
1085 we_mem_stag_s, sn_mem_s, &
1086 we_mem_stag_e, sn_mem_e, &
1087 we_mem_s, sn_mem_stag_s, &
1088 we_mem_e, sn_mem_stag_e, &
1089 xlon_u, xlon_v, xlat_u, xlat_v)
1090 else if (gridtype == 'E') then
1091 call map_to_met_nmm(u_field%r_arr, u_field%modified_mask, &
1092 v_field%r_arr, v_field%modified_mask, &
1093 we_mem_s, sn_mem_s, &
1094 we_mem_e, sn_mem_e, &
1095 xlat_v, xlon_v)
1096 end if
1097 end if
1099 call bitarray_destroy(u_field%modified_mask)
1100 call bitarray_destroy(v_field%modified_mask)
1101 nullify(u_field%modified_mask)
1102 nullify(v_field%modified_mask)
1103 call storage_put_field(u_field)
1104 call storage_put_field(v_field)
1105 end if
1107 end do
1109 deallocate(u_levels)
1110 deallocate(v_levels)
1112 end if
1114 call pop_source_projection()
1116 else
1117 if (do_const_processing) then
1118 call mprintf(.true.,WARN,'Couldn''t open file %s for input.',s1=input_name)
1119 else
1120 call mprintf(.true.,WARN,'Couldn''t open file %s for input.',s1=trim(input_name)//':'//trim(temp_date))
1121 end if
1122 end if
1124 fg_idx = fg_idx + 1
1125 if (do_const_processing) then
1126 input_name = constants_name(fg_idx)
1127 else
1128 input_name = fg_name(fg_idx)
1129 end if
1130 end do ! while (input_name /= '*')
1132 !
1133 ! Rotate winds from earth-relative to grid-relative
1134 !
1136 call storage_get_levels(u_field, u_levels)
1137 call storage_get_levels(v_field, v_levels)
1139 if (associated(u_levels) .and. associated(v_levels)) then
1140 u_idx = 1
1141 do u_idx = 1, size(u_levels)
1142 u_field%header%vertical_level = u_levels(u_idx)
1143 call storage_get_field(u_field, istatus)
1144 v_field%header%vertical_level = v_levels(u_idx)
1145 call storage_get_field(v_field, istatus)
1147 if (gridtype == 'C') then
1148 call met_to_map(u_field%r_arr, u_field%valid_mask, &
1149 v_field%r_arr, v_field%valid_mask, &
1150 we_mem_stag_s, sn_mem_s, &
1151 we_mem_stag_e, sn_mem_e, &
1152 we_mem_s, sn_mem_stag_s, &
1153 we_mem_e, sn_mem_stag_e, &
1154 xlon_u, xlon_v, xlat_u, xlat_v)
1155 else if (gridtype == 'E') then
1156 call met_to_map_nmm(u_field%r_arr, u_field%valid_mask, &
1157 v_field%r_arr, v_field%valid_mask, &
1158 we_mem_s, sn_mem_s, &
1159 we_mem_e, sn_mem_e, &
1160 xlat_v, xlon_v)
1161 end if
1163 end do
1165 deallocate(u_levels)
1166 deallocate(v_levels)
1168 end if
1170 if (do_const_processing) return
1172 !
1173 ! Now that we have all degribbed fields, we build a 3-d pressure field, and fill in any
1174 ! missing levels in the other 3-d fields
1175 !
1176 call mprintf(.true.,LOGFILE,'Filling missing levels.')
1177 call fill_missing_levels(output_flags)
1179 call mprintf(.true.,LOGFILE,'Creating derived fields.')
1180 call create_derived_fields(gridtype, fg_data%hdate, fg_data%xfcst, &
1181 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
1182 we_mem_stag_s, we_mem_stag_e, sn_mem_stag_s, sn_mem_stag_e, &
1183 got_this_field, output_flags)
1185 !
1186 ! Check that every mandatory field was found in input data
1187 !
1188 do i=1,num_entries
1189 if (is_mandatory(i) .and. .not. got_this_field(i)) then
1190 call mprintf(.true.,ERROR,'The mandatory field %s was not found in any input data.',s1=fieldname(i))
1191 end if
1192 end do
1194 !
1195 ! Before we begin to write fields, if debug_level is set high enough, we
1196 ! write a table of which fields are available at which levels to the
1197 ! metgrid.log file, and then we check to see if any fields are not
1198 ! completely covered with data.
1199 !
1200 call storage_print_fields()
1201 call find_missing_values()
1203 !
1204 ! All of the processing is now done for this time period for this domain;
1205 ! now we simply output every field from the storage module.
1206 !
1208 title = 'OUTPUT FROM METGRID V3.7.1'
1210 ! Initialize the output module for this domain and time
1211 call mprintf(.true.,LOGFILE,'Initializing output module.')
1212 output_date = temp_date
1213 if ( .not. nocolons ) then
1214 if (len_trim(temp_date) == 13) then
1215 output_date(14:19) = ':00:00'
1216 else if (len_trim(temp_date) == 16) then
1217 output_date(17:19) = ':00'
1218 end if
1219 else
1220 if (len_trim(temp_date) == 13) then
1221 output_date(14:19) = '_00_00'
1222 else if (len_trim(temp_date) == 16) then
1223 output_date(17:19) = '_00'
1224 end if
1225 endif
1227 call output_init(n, title, output_date, gridtype, dyn_opt, &
1228 corner_lats, corner_lons, &
1229 we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
1230 we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, &
1231 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
1232 extra_col, extra_row)
1234 call get_bottom_top_dim(bottom_top_dim)
1236 ! Add in a flag to tell real that we have seven new msf fields
1237 nmet_flags = num_entries + 1
1238 if (associated(geogrid_flags)) nmet_flags = nmet_flags + size(geogrid_flags)
1239 allocate(met_flags(nmet_flags))
1240 do i=1,num_entries
1241 met_flags(i) = output_flags(i)
1242 end do
1243 if (gridtype == 'C') then
1244 met_flags(num_entries+1) = 'FLAG_MF_XY'
1245 else
1246 met_flags(num_entries+1) = ' '
1247 end if
1248 if (associated(geogrid_flags)) then
1249 do i=1,size(geogrid_flags)
1250 met_flags(num_entries+1+i) = geogrid_flags(i)
1251 end do
1252 end if
1254 ! Find out how many soil levels or layers we have; this assumes a field named ST
1255 field % header % field = 'ST'
1256 nullify(soil_levels)
1257 call storage_get_levels(field, soil_levels)
1259 if (.not. associated(soil_levels)) then
1260 field % header % field = 'SOILT'
1261 nullify(soil_levels)
1262 call storage_get_levels(field, soil_levels)
1263 end if
1265 if (.not. associated(soil_levels)) then
1266 field % header % field = 'STC_WPS'
1267 nullify(soil_levels)
1268 call storage_get_levels(field, soil_levels)
1269 end if
1271 if (associated(soil_levels)) then
1272 num_metgrid_soil_levs = size(soil_levels)
1273 deallocate(soil_levels)
1274 else
1275 num_metgrid_soil_levs = 0
1276 end if
1278 ! First write out global attributes
1279 call mprintf(.true.,LOGFILE,'Writing global attributes to output.')
1280 call write_global_attrs(title, output_date, gridtype, dyn_opt, west_east_dim, &
1281 south_north_dim, bottom_top_dim, &
1282 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
1283 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
1284 map_proj, mminlu, num_land_cat, &
1285 is_water, is_lake, is_ice, is_urban, i_soilwater, &
1286 grid_id, parent_id, i_parent_start, &
1287 j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
1288 cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
1289 pole_lat, pole_lon, parent_grid_ratio, sub_x, sub_y, &
1290 corner_lats, corner_lons, num_metgrid_soil_levs, &
1291 met_flags, nmet_flags, process_bdy_width)
1293 deallocate(met_flags)
1295 call reset_next_field()
1297 istatus = 0
1299 ! Now loop over all output fields, writing each to the output module
1300 do while (istatus == 0)
1301 call get_next_output_field(cname, real_array, &
1302 sm1, em1, sm2, em2, sm3, em3, istatus)
1303 if (istatus == 0) then
1305 call mprintf(.true.,LOGFILE,'Writing field %s to output.',s1=cname)
1306 call write_field(sm1, em1, sm2, em2, sm3, em3, &
1307 cname, output_date, real_array)
1308 deallocate(real_array)
1310 end if
1312 end do
1314 call mprintf(.true.,LOGFILE,'Closing output file.')
1315 call output_close()
1317 ! Free up memory used by met fields for this valid time
1318 call storage_delete_all_td()
1320 end subroutine process_single_met_time
1323 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1324 ! Name: get_interp_masks
1325 !
1326 ! Purpose:
1327 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1328 subroutine get_interp_masks(fg_prefix, is_constants, fg_date)
1330 use interp_option_module
1331 use read_met_module
1332 use storage_module
1334 implicit none
1336 ! Arguments
1337 logical, intent(in) :: is_constants
1338 character (len=*), intent(in) :: fg_prefix, fg_date
1340 ! BUG: Move this constant to misc_definitions_module?
1341 integer, parameter :: BDR_WIDTH = 3
1343 ! Local variables
1344 integer :: i, istatus, idx, idxt
1345 type (fg_input) :: mask_field
1346 type (met_data) :: fg_data
1348 istatus = 0
1350 call read_met_init(fg_prefix, is_constants, fg_date, istatus)
1352 do while (istatus == 0)
1354 call read_next_met_field(fg_data, istatus)
1356 if (istatus == 0) then
1358 ! Find out which METGRID.TBL entry goes with this field
1359 idxt = num_entries + 1
1360 do idx=1,num_entries
1361 if ((index(fieldname(idx), trim(fg_data%field)) /= 0) .and. &
1362 (len_trim(fieldname(idx)) == len_trim(fg_data%field))) then
1364 if (index(fg_prefix,trim(from_input(idx))) /= 0 .or. &
1365 (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
1366 idxt = idx
1367 end if
1369 end if
1370 end do
1371 idx = idxt
1372 if (idx > num_entries) idx = num_entries ! The last entry is a default
1374 ! Do we need to rename this field?
1375 if (output_name(idx) /= ' ') then
1376 fg_data%field = output_name(idx)(1:9)
1378 idxt = num_entries + 1
1379 do idx=1,num_entries
1380 if ((index(fieldname(idx), trim(fg_data%field)) /= 0) .and. &
1381 (len_trim(fieldname(idx)) == len_trim(fg_data%field))) then
1383 if (index(fg_prefix,trim(from_input(idx))) /= 0 .or. &
1384 (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
1385 idxt = idx
1386 end if
1388 end if
1389 end do
1390 idx = idxt
1391 if (idx > num_entries) idx = num_entries ! The last entry is a default
1392 end if
1394 do i=1,num_entries
1395 if (interp_mask(i) /= ' ' .and. (trim(interp_mask(i)) == trim(fg_data%field))) then
1397 mask_field%header%version = 1
1398 mask_field%header%date = ' '
1399 mask_field%header%date = fg_date
1400 if (is_constants) then
1401 mask_field%header%time_dependent = .false.
1402 mask_field%header%constant_field = .true.
1403 else
1404 mask_field%header%time_dependent = .true.
1405 mask_field%header%constant_field = .false.
1406 end if
1407 mask_field%header%mask_field = .true.
1408 mask_field%header%forecast_hour = 0.
1409 mask_field%header%fg_source = 'degribbed met data'
1410 mask_field%header%field = trim(fg_data%field)//'.mask'
1411 mask_field%header%units = '-'
1412 mask_field%header%description = '-'
1413 mask_field%header%vertical_coord = 'none'
1414 mask_field%header%vertical_level = 1
1415 mask_field%header%sr_x = 1
1416 mask_field%header%sr_y = 1
1417 mask_field%header%array_order = 'XY'
1418 mask_field%header%dim1(1) = 1
1419 mask_field%header%dim1(2) = fg_data%nx
1420 mask_field%header%dim2(1) = 1
1421 mask_field%header%dim2(2) = fg_data%ny
1422 mask_field%header%is_wind_grid_rel = .true.
1423 mask_field%header%array_has_missing_values = .false.
1424 mask_field%map%stagger = M
1426 ! Do a simple check to see whether this is a global lat/lon dataset
1427 ! Note that we do not currently support regional Gaussian grids
1428 if ((fg_data%iproj == PROJ_LATLON .and. abs(fg_data%nx * fg_data%deltalon - 360.) < 0.0001) &
1429 .or. (fg_data%iproj == PROJ_GAUSS)) then
1430 allocate(mask_field%r_arr(1-BDR_WIDTH:fg_data%nx+BDR_WIDTH,1:fg_data%ny))
1432 mask_field%r_arr(1:fg_data%nx, 1:fg_data%ny) = &
1433 fg_data%slab(1:fg_data%nx, 1:fg_data%ny)
1435 mask_field%r_arr(1-BDR_WIDTH:0, 1:fg_data%ny) = &
1436 fg_data%slab(fg_data%nx-BDR_WIDTH+1:fg_data%nx, 1:fg_data%ny)
1438 mask_field%r_arr(fg_data%nx+1:fg_data%nx+BDR_WIDTH, 1:fg_data%ny) = &
1439 fg_data%slab(1:BDR_WIDTH, 1:fg_data%ny)
1441 else
1442 allocate(mask_field%r_arr(1:fg_data%nx,1:fg_data%ny))
1443 mask_field%r_arr = fg_data%slab
1444 end if
1446 nullify(mask_field%valid_mask)
1447 nullify(mask_field%modified_mask)
1449 call storage_put_field(mask_field)
1451 exit
1453 end if
1454 end do
1456 if (associated(fg_data%slab)) deallocate(fg_data%slab)
1458 end if
1460 end do
1462 call read_met_close()
1464 end subroutine get_interp_masks
1467 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1468 ! Name: interp_met_field
1469 !
1470 ! Purpose:
1471 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1472 subroutine interp_met_field(input_name, short_fieldnm, ifieldstagger, istagger, &
1473 field, xlat, xlon, sm1, em1, sm2, em2, &
1474 sd1, ed1, sd2, ed2, &
1475 slab, minx, maxx, miny, maxy, bdr, do_gcell_interp, &
1476 new_pts, process_bdy_width, landmask)
1478 use bitarray_module
1479 use interp_module
1480 use interp_option_module
1481 use llxy_module
1482 use misc_definitions_module
1483 use storage_module
1485 implicit none
1487 ! Arguments
1488 integer, intent(in) :: ifieldstagger, istagger, &
1489 sm1, em1, sm2, em2, &
1490 sd1, ed1, sd2, ed2, &
1491 minx, maxx, miny, maxy, bdr, &
1492 process_bdy_width
1493 real, dimension(minx:maxx,miny:maxy), intent(in) :: slab
1494 real, dimension(sm1:em1,sm2:em2), intent(in) :: xlat, xlon
1495 real, dimension(sm1:em1,sm2:em2), intent(in), optional :: landmask
1496 logical, intent(in) :: do_gcell_interp
1497 character (len=9), intent(in) :: short_fieldnm
1498 character (len=MAX_FILENAME_LEN), intent(in) :: input_name
1499 type (fg_input), intent(inout) :: field
1500 type (bitarray), intent(inout) :: new_pts
1502 ! Local variables
1503 integer :: i, j, idx, idxt, orig_selected_proj, interp_mask_status, &
1504 interp_land_mask_status, interp_water_mask_status, process_width
1505 integer, pointer, dimension(:) :: interp_array, interp_opts
1506 real :: rx, ry, temp
1507 real, pointer, dimension(:,:) :: data_count
1508 type (fg_input) :: mask_field, mask_water_field, mask_land_field
1509 !BPR BEGIN
1510 real, dimension(sm1:em1,sm2:em2) :: r_arr_cur_source
1511 !BPR END
1513 ! CWH Initialize local pointer variables
1514 nullify(interp_array)
1515 nullify(interp_opts)
1516 nullify(data_count)
1518 ! Find index into fieldname, interp_method, masked, and fill_missing
1519 ! of the current field
1520 idxt = num_entries + 1
1521 do idx=1,num_entries
1522 if ((index(fieldname(idx), trim(short_fieldnm)) /= 0) .and. &
1523 (len_trim(fieldname(idx)) == len_trim(short_fieldnm))) then
1524 if (index(input_name,trim(from_input(idx))) /= 0 .or. &
1525 (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
1526 idxt = idx
1527 end if
1528 end if
1529 end do
1530 idx = idxt
1531 if (idx > num_entries) then
1532 call mprintf(.true.,WARN,'Entry in METGRID.TBL not found for field %s. '// &
1533 'Default options will be used for this field!', s1=short_fieldnm)
1534 idx = num_entries ! The last entry is a default
1535 end if
1537 if (process_bdy_width == 0) then
1538 process_width = max(ed1-sd1+1, ed2-sd2+1)
1539 else
1540 process_width = process_bdy_width
1541 ! Add two extra rows/cols to accommodate staggered fields: one extra row/col for
1542 ! averaging to mass points in real, and one beyond that for averaging during
1543 ! wind rotation
1544 if (ifieldstagger /= M) process_width = process_width + 2
1545 end if
1547 field%header%dim1(1) = sm1
1548 field%header%dim1(2) = em1
1549 field%header%dim2(1) = sm2
1550 field%header%dim2(2) = em2
1551 field%map%stagger = ifieldstagger
1552 if (.not. associated(field%r_arr)) then
1553 allocate(field%r_arr(sm1:em1,sm2:em2))
1554 end if
1556 interp_mask_status = 1
1557 interp_land_mask_status = 1
1558 interp_water_mask_status = 1
1560 if (interp_mask(idx) /= ' ') then
1561 mask_field%header%version = 1
1562 mask_field%header%forecast_hour = 0.
1563 mask_field%header%field = trim(interp_mask(idx))//'.mask'
1564 mask_field%header%vertical_coord = 'none'
1565 mask_field%header%vertical_level = 1
1567 call storage_get_field(mask_field, interp_mask_status)
1569 end if
1570 if (interp_land_mask(idx) /= ' ') then
1571 mask_land_field%header%version = 1
1572 mask_land_field%header%forecast_hour = 0.
1573 mask_land_field%header%field = trim(interp_land_mask(idx))//'.mask'
1574 mask_land_field%header%vertical_coord = 'none'
1575 mask_land_field%header%vertical_level = 1
1577 call storage_get_field(mask_land_field, interp_land_mask_status)
1579 end if
1580 if (interp_water_mask(idx) /= ' ') then
1581 mask_water_field%header%version = 1
1582 mask_water_field%header%forecast_hour = 0.
1583 mask_water_field%header%field = trim(interp_water_mask(idx))//'.mask'
1584 mask_water_field%header%vertical_coord = 'none'
1585 mask_water_field%header%vertical_level = 1
1587 call storage_get_field(mask_water_field, interp_water_mask_status)
1589 end if
1591 interp_array => interp_array_from_string(interp_method(idx))
1592 interp_opts => interp_options_from_string(interp_method(idx))
1595 !
1596 ! Interpolate using average_gcell interpolation method
1597 !
1598 if (do_gcell_interp) then
1599 !BPR BEGIN
1600 !If a lower priority source of the current field has already been read
1601 !in, the results of that input are currently in field%r_arr
1602 !Pass COPY of field%r_arr into accum_continous because in accum_continuous
1603 !it will set the input variable to zero over points covered by the
1604 !current source and then determine the appropriate value to place at
1605 !that point. This will overwrite data already put in field%r_arr by
1606 !lower priority sources.
1607 !This is problematic if the current source results in missing values
1608 !over parts of the area covered by the current source where a lower
1609 !priority source has already provided a value. In this case, if one
1610 !passes in field%r_arr, it will overwrite the value provided by the
1611 !lower priority source with zero.
1612 r_arr_cur_source = field%r_arr
1613 !BPR END
1614 allocate(data_count(sm1:em1,sm2:em2))
1615 data_count = 0.
1617 if (interp_mask_status == 0) then
1618 call accum_continuous(slab, &
1619 minx, maxx, miny, maxy, 1, 1, bdr, &
1620 !BPR BEGIN
1621 !Pass COPY of field%r_arr instead of field%r_arr itself
1622 !field%r_arr, data_count, &
1623 r_arr_cur_source, data_count, &
1624 !BPR END
1625 sm1, em1, sm2, em2, 1, 1, &
1626 istagger, &
1627 new_pts, missing_value(idx), interp_mask_val(idx), interp_mask_relational(idx), mask_field%r_arr)
1628 else
1629 call accum_continuous(slab, &
1630 minx, maxx, miny, maxy, 1, 1, bdr, &
1631 !BPR BEGIN
1632 !Pass COPY of field%r_arr instead of field%r_arr itself
1633 !field%r_arr, data_count, &
1634 r_arr_cur_source, data_count, &
1635 !BPR END
1636 sm1, em1, sm2, em2, 1, 1, &
1637 istagger, &
1638 new_pts, missing_value(idx), -1.) ! The -1 is the maskval, but since we
1639 ! we do not give an optional mask, no
1640 ! no need to worry about -1s in data
1641 end if
1643 orig_selected_proj = iget_selected_domain()
1644 call select_domain(SOURCE_PROJ)
1645 do j=sm2,em2
1646 do i=sm1,em1
1648 if (abs(i - sd1) >= process_width .and. (abs(i - ed1) >= process_width) .and. &
1649 abs(j - sd2) >= process_width .and. (abs(j - ed2) >= process_width)) then
1650 field%r_arr(i,j) = fill_missing(idx)
1651 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1652 cycle
1653 end if
1655 if (present(landmask)) then
1657 if (landmask(i,j) /= masked(idx)) then
1658 if (data_count(i,j) > 0.) then
1659 !BPR BEGIN
1660 !accum_continuous is now passed a copy of field%r_arr (r_arr_cur_source)
1661 !instead of field%r_arr itself so that it does not set
1662 !field%r_arr to zero where the input source is marked as missing
1663 !field%r_arr(i,j) = field%r_arr(i,j) / data_count(i,j)
1664 field%r_arr(i,j) = r_arr_cur_source(i,j) / data_count(i,j)
1665 !BPR END
1666 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1667 else
1669 if (interp_mask_status == 0) then
1670 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, 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, interp_opts, slab, &
1676 minx, maxx, miny, maxy, bdr, missing_value(idx))
1677 end if
1679 if (temp /= missing_value(idx)) then
1680 field%r_arr(i,j) = temp
1681 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1682 end if
1684 end if
1685 else
1686 field%r_arr(i,j) = fill_missing(idx)
1687 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1688 end if
1690 if (.not. bitarray_test(new_pts, i-sm1+1, j-sm2+1) .and. &
1691 .not. bitarray_test(field%valid_mask, i-sm1+1, j-sm2+1)) then
1692 field%r_arr(i,j) = fill_missing(idx)
1694 ! Assume that if missing fill value is other than default, then user has asked
1695 ! to fill in any missing values, and we can consider this point to have
1696 ! received a valid value
1697 if (fill_missing(idx) /= NAN) call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1698 end if
1700 else
1702 if (data_count(i,j) > 0.) then
1703 !BPR BEGIN
1704 !accum_continuous is now passed a copy of field%r_arr (r_arr_cur_source)
1705 !instead of field%r_arr itself so that it does not set
1706 !field%r_arr to zero where the input source is marked as missing
1707 !field%r_arr(i,j) = field%r_arr(i,j) / data_count(i,j)
1708 field%r_arr(i,j) = r_arr_cur_source(i,j) / data_count(i,j)
1709 !BPR END
1710 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1711 else
1713 if (interp_mask_status == 0) then
1714 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
1715 minx, maxx, miny, maxy, bdr, missing_value(idx), &
1716 mask_relational=interp_mask_relational(idx), &
1717 mask_val=interp_mask_val(idx), mask_field=mask_field%r_arr)
1718 else
1719 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
1720 minx, maxx, miny, maxy, bdr, missing_value(idx))
1721 end if
1723 if (temp /= missing_value(idx)) then
1724 field%r_arr(i,j) = temp
1725 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1726 end if
1728 end if
1730 if (.not. bitarray_test(new_pts, i-sm1+1, j-sm2+1) .and. &
1731 .not. bitarray_test(field%valid_mask, i-sm1+1, j-sm2+1)) then
1732 field%r_arr(i,j) = fill_missing(idx)
1734 ! Assume that if missing fill value is other than default, then user has asked
1735 ! to fill in any missing values, and we can consider this point to have
1736 ! received a valid value
1737 if (fill_missing(idx) /= NAN) call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1738 end if
1740 end if
1742 end do
1743 end do
1744 call select_domain(orig_selected_proj)
1745 deallocate(data_count)
1747 !
1748 ! No average_gcell interpolation method
1749 !
1750 else
1752 orig_selected_proj = iget_selected_domain()
1753 call select_domain(SOURCE_PROJ)
1754 do j=sm2,em2
1755 do i=sm1,em1
1757 if (abs(i - sd1) >= process_width .and. (abs(i - ed1) >= process_width) .and. &
1758 abs(j - sd2) >= process_width .and. (abs(j - ed2) >= process_width)) then
1759 field%r_arr(i,j) = fill_missing(idx)
1760 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1761 cycle
1762 end if
1764 if (present(landmask)) then
1766 if (masked(idx) == MASKED_BOTH) then
1768 if (landmask(i,j) == 0) then ! WATER POINT
1770 if (interp_land_mask_status == 0) then
1771 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
1772 minx, maxx, miny, maxy, bdr, missing_value(idx), &
1773 mask_relational=interp_land_mask_relational(idx), &
1774 mask_val=interp_land_mask_val(idx), mask_field=mask_land_field%r_arr)
1775 else
1776 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
1777 minx, maxx, miny, maxy, bdr, missing_value(idx))
1778 end if
1780 else if (landmask(i,j) == 1) then ! LAND POINT
1782 if (interp_water_mask_status == 0) then
1783 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
1784 minx, maxx, miny, maxy, bdr, missing_value(idx), &
1785 mask_relational=interp_water_mask_relational(idx), &
1786 mask_val=interp_water_mask_val(idx), mask_field=mask_water_field%r_arr)
1787 else
1788 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
1789 minx, maxx, miny, maxy, bdr, missing_value(idx))
1790 end if
1792 end if
1794 else if (landmask(i,j) /= masked(idx)) then
1796 if (interp_mask_status == 0) then
1797 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
1798 minx, maxx, miny, maxy, bdr, missing_value(idx), &
1799 mask_relational=interp_mask_relational(idx), &
1800 mask_val=interp_mask_val(idx), mask_field=mask_field%r_arr)
1801 else
1802 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
1803 minx, maxx, miny, maxy, bdr, missing_value(idx))
1804 end if
1806 else
1807 temp = missing_value(idx)
1808 end if
1810 ! No landmask for this field
1811 else
1813 if (interp_mask_status == 0) then
1814 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
1815 minx, maxx, miny, maxy, bdr, missing_value(idx), &
1816 mask_relational=interp_mask_relational(idx), &
1817 mask_val=interp_mask_val(idx), mask_field=mask_field%r_arr)
1818 else
1819 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, interp_opts, slab, &
1820 minx, maxx, miny, maxy, bdr, missing_value(idx))
1821 end if
1823 end if
1825 if (temp /= missing_value(idx)) then
1826 field%r_arr(i,j) = temp
1827 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1828 else if (present(landmask)) then
1829 if (landmask(i,j) == masked(idx)) then
1830 field%r_arr(i,j) = fill_missing(idx)
1831 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1832 end if
1833 end if
1835 if (.not. bitarray_test(new_pts, i-sm1+1, j-sm2+1) .and. &
1836 .not. bitarray_test(field%valid_mask, i-sm1+1, j-sm2+1)) then
1837 field%r_arr(i,j) = fill_missing(idx)
1839 ! Assume that if missing fill value is other than default, then user has asked
1840 ! to fill in any missing values, and we can consider this point to have
1841 ! received a valid value
1842 if (fill_missing(idx) /= NAN) call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1843 end if
1845 end do
1846 end do
1847 call select_domain(orig_selected_proj)
1848 end if
1850 deallocate(interp_array)
1851 deallocate(interp_opts)
1853 end subroutine interp_met_field
1856 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1857 ! Name: interp_to_latlon
1858 !
1859 ! Purpose:
1860 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1861 function interp_to_latlon(rlat, rlon, istagger, interp_method_list, interp_opt_list, slab, &
1862 minx, maxx, miny, maxy, bdr, source_missing_value, &
1863 mask_field, mask_relational, mask_val)
1865 use interp_module
1866 use llxy_module
1868 implicit none
1870 ! Arguments
1871 integer, intent(in) :: minx, maxx, miny, maxy, bdr, istagger
1872 integer, dimension(:), intent(in) :: interp_method_list
1873 integer, dimension(:), intent(in) :: interp_opt_list
1874 real, intent(in) :: rlat, rlon, source_missing_value
1875 real, dimension(minx:maxx,miny:maxy), intent(in) :: slab
1876 real, intent(in), optional :: mask_val
1877 real, dimension(minx:maxx,miny:maxy), intent(in), optional :: mask_field
1878 character(len=1), intent(in), optional :: mask_relational
1880 ! Return value
1881 real :: interp_to_latlon
1883 ! Local variables
1884 real :: rx, ry
1886 interp_to_latlon = source_missing_value
1888 call lltoxy(rlat, rlon, rx, ry, istagger)
1889 if (rx >= minx+bdr-0.5 .and. rx <= maxx-bdr+0.5) then
1890 if (present(mask_field) .and. present(mask_val) .and. present (mask_relational)) then
1891 interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, 1, 1, source_missing_value, &
1892 interp_method_list, interp_opt_list, 1, mask_relational, mask_val, mask_field)
1893 else if (present(mask_field) .and. present(mask_val)) then
1894 interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, 1, 1, source_missing_value, &
1895 interp_method_list, interp_opt_list, 1, maskval=mask_val, mask_array=mask_field)
1896 else
1897 interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, 1, 1, source_missing_value, &
1898 interp_method_list, interp_opt_list, 1)
1899 end if
1900 else
1901 interp_to_latlon = source_missing_value
1902 end if
1904 if (interp_to_latlon == source_missing_value) then
1906 ! Try a lon in the range 0. to 360.; all lons in the xlon
1907 ! array should be in the range -180. to 180.
1908 if (rlon < 0.) then
1909 call lltoxy(rlat, rlon+360., rx, ry, istagger)
1910 if (rx >= minx+bdr-0.5 .and. rx <= maxx-bdr+0.5) then
1911 if (present(mask_field) .and. present(mask_val) .and. present(mask_relational)) then
1912 interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, &
1913 1, 1, source_missing_value, &
1914 interp_method_list, interp_opt_list, 1, &
1915 mask_relational, mask_val, mask_field)
1916 else if (present(mask_field) .and. present(mask_val)) then
1917 interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, &
1918 1, 1, source_missing_value, &
1919 interp_method_list, interp_opt_list, 1, &
1920 maskval=mask_val, mask_array=mask_field)
1921 else
1922 interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, &
1923 1, 1, source_missing_value, &
1924 interp_method_list, interp_opt_list, 1)
1925 end if
1926 else
1927 interp_to_latlon = source_missing_value
1928 end if
1930 end if
1932 end if
1934 return
1936 end function interp_to_latlon
1939 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1940 ! Name: get_bottom_top_dim
1941 !
1942 ! Purpose:
1943 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1944 subroutine get_bottom_top_dim(bottom_top_dim)
1946 use interp_option_module
1947 use list_module
1948 use storage_module
1950 implicit none
1952 ! Arguments
1953 integer, intent(out) :: bottom_top_dim
1955 ! Local variables
1956 integer :: i, j
1957 integer, pointer, dimension(:) :: field_levels
1958 character (len=32) :: z_dim
1959 type (fg_input), pointer, dimension(:) :: headers
1960 type (list) :: temp_levels
1962 !CWH Initialize local pointer variables
1963 nullify(field_levels)
1964 nullify(headers)
1966 ! Initialize a list to store levels that are found for 3-d fields
1967 call list_init(temp_levels)
1969 ! Get a list of all time-dependent fields (given by their headers) from
1970 ! the storage module
1971 call storage_get_td_headers(headers)
1973 !
1974 ! Given headers of all fields, we first build a list of all possible levels
1975 ! for 3-d met fields (excluding sea-level, though).
1976 !
1977 do i=1,size(headers)
1978 call get_z_dim_name(headers(i)%header%field, z_dim)
1980 ! We only want to consider 3-d met fields
1981 if (z_dim(1:18) == 'num_metgrid_levels') then
1983 ! Find out what levels the current field has
1984 call storage_get_levels(headers(i), field_levels)
1985 do j=1,size(field_levels)
1987 ! If this level has not yet been encountered, add it to our list
1988 if (.not. list_search(temp_levels, ikey=field_levels(j), ivalue=field_levels(j))) then
1989 if (field_levels(j) /= 201300) then
1990 call list_insert(temp_levels, ikey=field_levels(j), ivalue=field_levels(j))
1991 end if
1992 end if
1994 end do
1996 deallocate(field_levels)
1998 end if
2000 end do
2002 bottom_top_dim = list_length(temp_levels)
2004 call list_destroy(temp_levels)
2005 deallocate(headers)
2007 end subroutine get_bottom_top_dim
2010 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2011 ! Name: fill_missing_levels
2012 !
2013 ! Purpose:
2014 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2015 subroutine fill_missing_levels(output_flags)
2017 use interp_option_module
2018 use list_module
2019 use module_debug
2020 use module_mergesort
2021 use storage_module
2023 implicit none
2025 ! Arguments
2026 character (len=128), dimension(:), intent(inout) :: output_flags
2028 ! Local variables
2029 integer :: i, ii, j, ix, jx, k, lower, upper, temp, istatus
2030 integer, pointer, dimension(:) :: union_levels, field_levels
2031 real, pointer, dimension(:) :: r_union_levels
2032 character (len=128) :: clevel
2033 type (fg_input) :: lower_field, upper_field, new_field, search_field
2034 type (fg_input), pointer, dimension(:) :: headers, all_headers
2035 type (list) :: temp_levels
2036 type (list_item), pointer, dimension(:) :: keys
2038 ! CWH Initialize local pointer variables
2039 nullify(union_levels)
2040 nullify(field_levels)
2041 nullify(r_union_levels)
2042 nullify(headers)
2043 nullify(all_headers)
2044 nullify(keys)
2046 ! Initialize a list to store levels that are found for 3-d fields
2047 call list_init(temp_levels)
2049 ! Get a list of all fields (given by their headers) from the storage module
2050 call storage_get_td_headers(headers)
2051 call storage_get_all_headers(all_headers)
2053 !
2054 ! Given headers of all fields, we first build a list of all possible levels
2055 ! for 3-d met fields (excluding sea-level, though).
2056 !
2057 do i=1,size(headers)
2059 ! Find out what levels the current field has
2060 call storage_get_levels(headers(i), field_levels)
2061 do j=1,size(field_levels)
2063 ! If this level has not yet been encountered, add it to our list
2064 if (.not. list_search(temp_levels, ikey=field_levels(j), ivalue=field_levels(j))) then
2065 if (field_levels(j) /= 201300) then
2066 call list_insert(temp_levels, ikey=field_levels(j), ivalue=field_levels(j))
2067 end if
2068 end if
2070 end do
2072 deallocate(field_levels)
2074 end do
2076 if (list_length(temp_levels) > 0) then
2078 !
2079 ! With all possible levels stored in a list, get an array of levels, sorted
2080 ! in decreasing order
2081 !
2082 i = 0
2083 allocate(union_levels(list_length(temp_levels)))
2084 do while (list_length(temp_levels) > 0)
2085 i = i + 1
2086 call list_get_first_item(temp_levels, ikey=union_levels(i), ivalue=temp)
2087 end do
2088 call mergesort(union_levels, 1, size(union_levels))
2090 allocate(r_union_levels(size(union_levels)))
2091 do i=1,size(union_levels)
2092 r_union_levels(i) = real(union_levels(i))
2093 end do
2095 !
2096 ! With a sorted, complete list of levels, we need
2097 ! to go back and fill in missing levels for each 3-d field
2098 !
2099 do i=1,size(headers)
2101 !
2102 ! Find entry in METGRID.TBL for this field, if one exists; if it does, then the
2103 ! entry may tell us how to get values for the current field at the missing level
2104 !
2105 do ii=1,num_entries
2106 if (fieldname(ii) == headers(i)%header%field) exit
2107 end do
2108 if (ii <= num_entries) then
2109 call dup(headers(i),new_field)
2110 nullify(new_field%valid_mask)
2111 nullify(new_field%modified_mask)
2112 call fill_field(new_field, ii, output_flags, r_union_levels)
2113 end if
2115 end do
2117 deallocate(union_levels)
2118 deallocate(r_union_levels)
2119 deallocate(headers)
2121 call storage_get_td_headers(headers)
2123 !
2124 ! Now we may need to vertically interpolate to missing values in 3-d fields
2125 !
2126 do i=1,size(headers)
2128 call storage_get_levels(headers(i), field_levels)
2130 ! If this isn't a 3-d array, nothing to do
2131 if (size(field_levels) > 1) then
2133 do k=1,size(field_levels)
2134 call dup(headers(i),search_field)
2135 search_field%header%vertical_level = field_levels(k)
2136 call storage_get_field(search_field,istatus)
2137 if (istatus == 0) then
2138 JLOOP: do jx=search_field%header%dim2(1),search_field%header%dim2(2)
2139 ILOOP: do ix=search_field%header%dim1(1),search_field%header%dim1(2)
2140 if (.not. bitarray_test(search_field%valid_mask, &
2141 ix-search_field%header%dim1(1)+1, &
2142 jx-search_field%header%dim2(1)+1)) then
2144 call dup(search_field, lower_field)
2145 do lower=k-1,1,-1
2146 lower_field%header%vertical_level = field_levels(lower)
2147 call storage_get_field(lower_field,istatus)
2148 if (bitarray_test(lower_field%valid_mask, &
2149 ix-search_field%header%dim1(1)+1, &
2150 jx-search_field%header%dim2(1)+1)) &
2151 exit
2153 end do
2155 call dup(search_field, upper_field)
2156 do upper=k+1,size(field_levels)
2157 upper_field%header%vertical_level = field_levels(upper)
2158 call storage_get_field(upper_field,istatus)
2159 if (bitarray_test(upper_field%valid_mask, &
2160 ix-search_field%header%dim1(1)+1, &
2161 jx-search_field%header%dim2(1)+1)) &
2162 exit
2164 end do
2165 if (upper <= size(field_levels) .and. lower >= 1) then
2166 search_field%r_arr(ix,jx) = real(abs(field_levels(upper)-field_levels(k))) &
2167 / real(abs(field_levels(upper)-field_levels(lower))) &
2168 * lower_field%r_arr(ix,jx) &
2169 + real(abs(field_levels(k)-field_levels(lower))) &
2170 / real(abs(field_levels(upper)-field_levels(lower))) &
2171 * upper_field%r_arr(ix,jx)
2172 call bitarray_set(search_field%valid_mask, &
2173 ix-search_field%header%dim1(1)+1, &
2174 jx-search_field%header%dim2(1)+1)
2175 end if
2176 end if
2177 end do ILOOP
2178 end do JLOOP
2179 else
2180 call mprintf(.true.,ERROR, &
2181 'This is bad, could not get %s at level %i.', &
2182 s1=trim(search_field%header%field), i1=field_levels(k))
2183 end if
2184 end do
2186 end if
2188 deallocate(field_levels)
2190 end do
2192 end if
2194 call list_destroy(temp_levels)
2195 deallocate(all_headers)
2196 deallocate(headers)
2198 end subroutine fill_missing_levels
2201 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2202 ! Name: create_derived_fields
2203 !
2204 ! Purpose:
2205 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2206 subroutine create_derived_fields(arg_gridtype, hdate, xfcst, &
2207 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
2208 we_mem_stag_s, we_mem_stag_e, sn_mem_stag_s, sn_mem_stag_e, &
2209 created_this_field, output_flags)
2211 use interp_option_module
2212 use list_module
2213 use module_mergesort
2214 use storage_module
2216 implicit none
2218 ! Arguments
2219 integer, intent(in) :: we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
2220 we_mem_stag_s, we_mem_stag_e, sn_mem_stag_s, sn_mem_stag_e
2221 real, intent(in) :: xfcst
2222 logical, dimension(:), intent(inout) :: created_this_field
2223 character (len=1), intent(in) :: arg_gridtype
2224 character (len=24), intent(in) :: hdate
2225 character (len=128), dimension(:), intent(inout) :: output_flags
2227 ! Local variables
2228 integer :: idx, i, j, istatus
2229 type (fg_input) :: field
2231 ! Initialize fg_input structure to store the field
2232 field%header%version = 1
2233 field%header%date = hdate//' '
2234 field%header%time_dependent = .true.
2235 field%header%mask_field = .false.
2236 field%header%constant_field = .false.
2237 field%header%forecast_hour = xfcst
2238 field%header%fg_source = 'Derived from FG'
2239 field%header%field = ' '
2240 field%header%units = ' '
2241 field%header%description = ' '
2242 field%header%vertical_level = 0
2243 field%header%sr_x = 1
2244 field%header%sr_y = 1
2245 field%header%array_order = 'XY '
2246 field%header%is_wind_grid_rel = .true.
2247 field%header%array_has_missing_values = .false.
2248 nullify(field%r_arr)
2249 nullify(field%valid_mask)
2250 nullify(field%modified_mask)
2252 !
2253 ! Check each entry in METGRID.TBL to see whether it is a derive field
2254 !
2255 do idx=1,num_entries
2256 if (is_derived_field(idx)) then
2258 created_this_field(idx) = .true.
2260 call mprintf(.true.,INFORM,'Going to create the field %s',s1=fieldname(idx))
2262 ! Intialize more fields in storage structure
2263 field%header%field = fieldname(idx)
2264 call get_z_dim_name(fieldname(idx),field%header%vertical_coord)
2265 field%map%stagger = output_stagger(idx)
2266 if (arg_gridtype == 'E') then
2267 field%header%dim1(1) = we_mem_s
2268 field%header%dim1(2) = we_mem_e
2269 field%header%dim2(1) = sn_mem_s
2270 field%header%dim2(2) = sn_mem_e
2271 else if (arg_gridtype == 'C') then
2272 if (output_stagger(idx) == M) then
2273 field%header%dim1(1) = we_mem_s
2274 field%header%dim1(2) = we_mem_e
2275 field%header%dim2(1) = sn_mem_s
2276 field%header%dim2(2) = sn_mem_e
2277 else if (output_stagger(idx) == U) then
2278 field%header%dim1(1) = we_mem_stag_s
2279 field%header%dim1(2) = we_mem_stag_e
2280 field%header%dim2(1) = sn_mem_s
2281 field%header%dim2(2) = sn_mem_e
2282 else if (output_stagger(idx) == V) then
2283 field%header%dim1(1) = we_mem_s
2284 field%header%dim1(2) = we_mem_e
2285 field%header%dim2(1) = sn_mem_stag_s
2286 field%header%dim2(2) = sn_mem_stag_e
2287 end if
2288 end if
2290 call fill_field(field, idx, output_flags)
2292 end if
2293 end do
2296 end subroutine create_derived_fields
2299 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2300 ! Name: fill_field
2301 !
2302 ! Purpose:
2303 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2304 subroutine fill_field(field, idx, output_flags, all_level_list)
2306 use interp_option_module
2307 use list_module
2308 use module_mergesort
2309 use storage_module
2311 implicit none
2313 ! Arguments
2314 integer, intent(in) :: idx
2315 type (fg_input), intent(inout) :: field
2316 character (len=128), dimension(:), intent(inout) :: output_flags
2317 real, dimension(:), intent(in), optional :: all_level_list
2319 ! Local variables
2320 integer :: i, j, istatus, isrclevel
2321 integer, pointer, dimension(:) :: all_list
2322 real :: rfillconst, rlevel, rsrclevel
2323 type (fg_input) :: query_field
2324 type (list_item), pointer, dimension(:) :: keys
2325 character (len=128) :: asrcname
2326 logical :: filled_all_lev
2328 !CWH Initialize local pointer variables
2329 nullify(all_list)
2330 nullify(keys)
2332 filled_all_lev = .false.
2334 !
2335 ! Get a list of all levels to be filled for this field
2336 !
2337 keys => list_get_keys(fill_lev_list(idx))
2339 do i=1,list_length(fill_lev_list(idx))
2341 !
2342 ! First handle a specification for levels "all"
2343 !
2344 if (trim(keys(i)%ckey) == 'all') then
2346 ! We only want to fill all levels if we haven't already filled "all" of them
2347 if (.not. filled_all_lev) then
2349 filled_all_lev = .true.
2351 query_field%header%time_dependent = .true.
2352 query_field%header%field = ' '
2353 nullify(query_field%r_arr)
2354 nullify(query_field%valid_mask)
2355 nullify(query_field%modified_mask)
2357 ! See if we are filling this level with a constant
2358 call get_constant_fill_lev(keys(i)%cvalue, rfillconst, istatus)
2359 if (istatus == 0) then
2360 if (present(all_level_list)) then
2361 do j=1,size(all_level_list)
2362 call create_level(field, real(all_level_list(j)), idx, output_flags, rfillconst=rfillconst)
2363 end do
2364 else
2365 query_field%header%field = level_template(idx)
2366 nullify(all_list)
2367 call storage_get_levels(query_field, all_list)
2368 if (associated(all_list)) then
2369 do j=1,size(all_list)
2370 call create_level(field, real(all_list(j)), idx, output_flags, rfillconst=rfillconst)
2371 end do
2372 deallocate(all_list)
2373 end if
2374 end if
2376 ! Else see if we are filling this level with a constant equal
2377 ! to the value of the level
2378 else if (trim(keys(i)%cvalue) == 'vertical_index') then
2379 if (present(all_level_list)) then
2380 do j=1,size(all_level_list)
2381 call create_level(field, real(all_level_list(j)), idx, output_flags, &
2382 rfillconst=real(all_level_list(j)))
2383 end do
2384 else
2385 query_field%header%field = level_template(idx)
2386 nullify(all_list)
2387 call storage_get_levels(query_field, all_list)
2388 if (associated(all_list)) then
2389 do j=1,size(all_list)
2390 call create_level(field, real(all_list(j)), idx, output_flags, rfillconst=real(all_list(j)))
2391 end do
2392 deallocate(all_list)
2393 end if
2394 end if
2396 ! Else, we assume that it is a field from which we are copying levels
2397 else
2398 if (present(all_level_list)) then
2399 do j=1,size(all_level_list)
2400 call create_level(field, real(all_level_list(j)), idx, output_flags, &
2401 asrcname=keys(i)%cvalue, rsrclevel=real(all_level_list(j)))
2402 end do
2403 else
2404 query_field%header%field = keys(i)%cvalue ! Use same levels as source field, not level_template
2405 nullify(all_list)
2406 call storage_get_levels(query_field, all_list)
2407 if (associated(all_list)) then
2408 do j=1,size(all_list)
2409 call create_level(field, real(all_list(j)), idx, output_flags, &
2410 asrcname=keys(i)%cvalue, rsrclevel=real(all_list(j)))
2411 end do
2412 deallocate(all_list)
2414 else
2416 ! If the field doesn't have any levels (or does not exist) then we have not
2417 ! really filled all levels at this point.
2418 filled_all_lev = .false.
2419 end if
2420 end if
2422 end if
2423 end if
2425 !
2426 ! Handle individually specified levels
2427 !
2428 else
2430 read(keys(i)%ckey,*) rlevel
2432 ! See if we are filling this level with a constant
2433 call get_constant_fill_lev(keys(i)%cvalue, rfillconst, istatus)
2434 if (istatus == 0) then
2435 call create_level(field, rlevel, idx, output_flags, rfillconst=rfillconst)
2437 ! Otherwise, we are filling from another level
2438 else
2439 call get_fill_src_level(keys(i)%cvalue, asrcname, isrclevel)
2440 rsrclevel = real(isrclevel)
2441 call create_level(field, rlevel, idx, output_flags, &
2442 asrcname=asrcname, rsrclevel=rsrclevel)
2444 end if
2445 end if
2446 end do
2448 if (associated(keys)) deallocate(keys)
2450 end subroutine fill_field
2453 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2454 ! Name: create_level
2455 !
2456 ! Purpose:
2457 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2458 subroutine create_level(field_template, rlevel, idx, output_flags, &
2459 rfillconst, asrcname, rsrclevel)
2461 use storage_module
2462 use interp_option_module
2464 implicit none
2466 ! Arguments
2467 type (fg_input), intent(inout) :: field_template
2468 real, intent(in) :: rlevel
2469 integer, intent(in) :: idx
2470 character (len=128), dimension(:), intent(inout) :: output_flags
2471 real, intent(in), optional :: rfillconst, rsrclevel
2472 character (len=128), intent(in), optional :: asrcname
2474 ! Local variables
2475 integer :: i, j, istatus
2476 integer :: sm1,em1,sm2,em2
2477 type (fg_input) :: query_field
2479 !
2480 ! Check to make sure optional arguments are sane
2481 !
2482 if (present(rfillconst) .and. (present(asrcname) .or. present(rsrclevel))) then
2483 call mprintf(.true.,ERROR,'A call to create_level() cannot be given specifications '// &
2484 'for both a constant fill value and a source level.')
2486 else if ((present(asrcname) .and. .not. present(rsrclevel)) .or. &
2487 (.not. present(asrcname) .and. present(rsrclevel))) then
2488 call mprintf(.true.,ERROR,'Neither or both of optional arguments asrcname and '// &
2489 'rsrclevel must be specified to subroutine create_level().')
2491 else if (.not. present(rfillconst) .and. &
2492 .not. present(asrcname) .and. &
2493 .not. present(rsrclevel)) then
2494 call mprintf(.true.,ERROR,'A call to create_level() must be given either a specification '// &
2495 'for a constant fill value or a source level.')
2496 end if
2498 query_field%header%time_dependent = .true.
2499 query_field%header%field = field_template%header%field
2500 query_field%header%vertical_level = rlevel
2501 nullify(query_field%r_arr)
2502 nullify(query_field%valid_mask)
2503 nullify(query_field%modified_mask)
2505 call storage_query_field(query_field, istatus)
2506 if (istatus == 0) then
2507 call mprintf(.true.,INFORM,'%s at level %f already exists; leaving it alone.', &
2508 s1=field_template%header%field, f1=rlevel)
2509 return
2510 end if
2512 sm1 = field_template%header%dim1(1)
2513 em1 = field_template%header%dim1(2)
2514 sm2 = field_template%header%dim2(1)
2515 em2 = field_template%header%dim2(2)
2517 !
2518 ! Handle constant fill value case
2519 !
2520 if (present(rfillconst)) then
2522 field_template%header%vertical_level = rlevel
2523 allocate(field_template%r_arr(sm1:em1,sm2:em2))
2524 allocate(field_template%valid_mask)
2525 allocate(field_template%modified_mask)
2526 call bitarray_create(field_template%valid_mask, em1-sm1+1, em2-sm2+1)
2527 call bitarray_create(field_template%modified_mask, em1-sm1+1, em2-sm2+1)
2529 field_template%r_arr = rfillconst
2531 do j=sm2,em2
2532 do i=sm1,em1
2533 call bitarray_set(field_template%valid_mask, i-sm1+1, j-sm2+1)
2534 end do
2535 end do
2537 call storage_put_field(field_template)
2539 if (output_this_field(idx) .and. flag_in_output(idx) /= ' ') then
2540 output_flags(idx) = flag_in_output(idx)
2541 end if
2543 !
2544 ! Handle source field and source level case
2545 !
2546 else if (present(asrcname) .and. present(rsrclevel)) then
2548 query_field%header%field = ' '
2549 query_field%header%field = asrcname
2550 query_field%header%vertical_level = rsrclevel
2552 ! Check to see whether the requested source field exists at the requested level
2553 call storage_query_field(query_field, istatus)
2555 if (istatus == 0) then
2557 ! Read in requested field at requested level
2558 call storage_get_field(query_field, istatus)
2559 if ((query_field%header%dim1(1) /= field_template%header%dim1(1)) .or. &
2560 (query_field%header%dim1(2) /= field_template%header%dim1(2)) .or. &
2561 (query_field%header%dim2(1) /= field_template%header%dim2(1)) .or. &
2562 (query_field%header%dim2(2) /= field_template%header%dim2(2))) then
2563 call mprintf(.true.,ERROR,'Dimensions for %s do not match those of %s. This is '// &
2564 'probably because the staggerings of the fields do not match.', &
2565 s1=query_field%header%field, s2=field_template%header%field)
2566 end if
2568 field_template%header%vertical_level = rlevel
2569 allocate(field_template%r_arr(sm1:em1,sm2:em2))
2570 allocate(field_template%valid_mask)
2571 allocate(field_template%modified_mask)
2572 call bitarray_create(field_template%valid_mask, em1-sm1+1, em2-sm2+1)
2573 call bitarray_create(field_template%modified_mask, em1-sm1+1, em2-sm2+1)
2575 field_template%r_arr = query_field%r_arr
2577 ! We should retain information about which points in the field are valid
2578 do j=sm2,em2
2579 do i=sm1,em1
2580 if (bitarray_test(query_field%valid_mask, i-sm1+1, j-sm2+1)) then
2581 call bitarray_set(field_template%valid_mask, i-sm1+1, j-sm2+1)
2582 end if
2583 end do
2584 end do
2586 call storage_put_field(field_template)
2588 if (output_this_field(idx) .and. flag_in_output(idx) /= ' ') then
2589 output_flags(idx) = flag_in_output(idx)
2590 end if
2592 else
2593 call mprintf(.true.,INFORM,'Couldn''t find %s at level %f to fill level %f of %s.', &
2594 s1=asrcname,f1=rsrclevel,f2=rlevel,s2=field_template%header%field)
2595 end if
2597 end if
2599 end subroutine create_level
2602 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2603 ! Name: accum_continuous
2604 !
2605 ! Purpose: Sum up all of the source data points whose nearest neighbor in the
2606 ! model grid is the specified model grid point.
2607 !
2608 ! NOTE: When processing the source tile, those source points that are
2609 ! closest to a different model grid point will be added to the totals for
2610 ! such grid points; thus, an entire source tile will be processed at a time.
2611 ! This routine really processes for all model grid points that are
2612 ! within a source tile, and not just for a single grid point.
2613 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2614 subroutine accum_continuous(src_array, &
2615 src_min_x, src_max_x, src_min_y, src_max_y, src_min_z, src_max_z, bdr_width, &
2616 dst_array, n, &
2617 start_i, end_i, start_j, end_j, start_k, end_k, &
2618 istagger, &
2619 new_pts, msgval, maskval, mask_relational, mask_array, sr_x, sr_y)
2621 use bitarray_module
2622 use misc_definitions_module
2624 implicit none
2626 ! Arguments
2627 integer, intent(in) :: start_i, end_i, start_j, end_j, start_k, end_k, istagger, &
2628 src_min_x, src_max_x, src_min_y, src_max_y, src_min_z, src_max_z, bdr_width
2629 real, intent(in) :: maskval, msgval
2630 real, dimension(src_min_x:src_max_x, src_min_y:src_max_y, src_min_z:src_max_z), intent(in) :: src_array
2631 real, dimension(start_i:end_i, start_j:end_j, start_k:end_k), intent(inout) :: dst_array, n
2632 real, dimension(src_min_x:src_max_x, src_min_y:src_max_y), intent(in), optional :: mask_array
2633 integer, intent(in), optional :: sr_x, sr_y
2634 type (bitarray), intent(inout) :: new_pts
2635 character(len=1), intent(in), optional :: mask_relational
2637 ! Local variables
2638 integer :: i, j
2639 integer, pointer, dimension(:,:,:) :: where_maps_to
2640 real :: rsr_x, rsr_y
2642 rsr_x = 1.0
2643 rsr_y = 1.0
2644 if (present(sr_x)) rsr_x = real(sr_x)
2645 if (present(sr_y)) rsr_y = real(sr_y)
2647 allocate(where_maps_to(src_min_x:src_max_x,src_min_y:src_max_y,2))
2648 do i=src_min_x,src_max_x
2649 do j=src_min_y,src_max_y
2650 where_maps_to(i,j,1) = NOT_PROCESSED
2651 end do
2652 end do
2654 call process_continuous_block(src_array, where_maps_to, &
2655 src_min_x, src_min_y, src_min_z, src_max_x, src_max_y, src_max_z, &
2656 src_min_x+bdr_width, src_min_y, src_min_z, &
2657 src_max_x-bdr_width, src_max_y, src_max_z, &
2658 dst_array, n, start_i, end_i, start_j, end_j, start_k, end_k, &
2659 istagger, &
2660 new_pts, rsr_x, rsr_y, msgval, maskval, mask_relational, mask_array)
2662 deallocate(where_maps_to)
2664 end subroutine accum_continuous
2667 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2668 ! Name: process_continuous_block
2669 !
2670 ! Purpose: To recursively process a subarray of continuous data, adding the
2671 ! points in a block to the sum for their nearest grid point. The nearest
2672 ! neighbor may be estimated in some cases; for example, if the four corners
2673 ! of a subarray all have the same nearest grid point, all elements in the
2674 ! subarray are added to that grid point.
2675 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2676 recursive subroutine process_continuous_block(tile_array, where_maps_to, &
2677 src_min_x, src_min_y, src_min_z, src_max_x, src_max_y, src_max_z, &
2678 min_i, min_j, min_k, max_i, max_j, max_k, &
2679 dst_array, n, &
2680 start_x, end_x, start_y, end_y, start_z, end_z, &
2681 istagger, &
2682 new_pts, sr_x, sr_y, msgval, maskval, mask_relational, mask_array)
2684 use bitarray_module
2685 use llxy_module
2686 use misc_definitions_module
2688 implicit none
2690 ! Arguments
2691 integer, intent(in) :: min_i, min_j, min_k, max_i, max_j, max_k, &
2692 src_min_x, src_min_y, src_min_z, src_max_x, src_max_y, src_max_z, &
2693 start_x, end_x, start_y, end_y, start_z, end_z, istagger
2694 integer, dimension(src_min_x:src_max_x,src_min_y:src_max_y,2), intent(inout) :: where_maps_to
2695 real, intent(in) :: sr_x, sr_y, maskval, msgval
2696 real, dimension(src_min_x:src_max_x,src_min_y:src_max_y,src_min_z:src_max_z), intent(in) :: tile_array
2697 real, dimension(src_min_x:src_max_x,src_min_y:src_max_y), intent(in), optional :: mask_array
2698 real, dimension(start_x:end_x,start_y:end_y,start_z:end_z), intent(inout) :: dst_array, n
2699 type (bitarray), intent(inout) :: new_pts
2700 character(len=1), intent(in), optional :: mask_relational
2702 ! Local variables
2703 integer :: orig_selected_domain, x_dest, y_dest, i, j, k, center_i, center_j
2704 real :: lat_corner, lon_corner, rx, ry
2706 ! Compute the model grid point that the corners of the rectangle to be
2707 ! processed map to
2708 ! Lower-left corner
2709 if (where_maps_to(min_i,min_j,1) == NOT_PROCESSED) then
2710 orig_selected_domain = iget_selected_domain()
2711 call select_domain(SOURCE_PROJ)
2712 call xytoll(real(min_i), real(min_j), lat_corner, lon_corner, istagger)
2713 call select_domain(1)
2714 call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
2715 rx = (rx - 1.0)*sr_x + 1.0
2716 ry = (ry - 1.0)*sr_y + 1.0
2717 call select_domain(orig_selected_domain)
2718 if (real(start_x) <= rx .and. rx <= real(end_x) .and. &
2719 real(start_y) <= ry .and. ry <= real(end_y)) then
2720 where_maps_to(min_i,min_j,1) = nint(rx)
2721 where_maps_to(min_i,min_j,2) = nint(ry)
2722 else
2723 where_maps_to(min_i,min_j,1) = OUTSIDE_DOMAIN
2724 end if
2725 end if
2727 ! Upper-left corner
2728 if (where_maps_to(min_i,max_j,1) == NOT_PROCESSED) then
2729 orig_selected_domain = iget_selected_domain()
2730 call select_domain(SOURCE_PROJ)
2731 call xytoll(real(min_i), real(max_j), lat_corner, lon_corner, istagger)
2732 call select_domain(1)
2733 call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
2734 rx = (rx - 1.0)*sr_x + 1.0
2735 ry = (ry - 1.0)*sr_y + 1.0
2736 call select_domain(orig_selected_domain)
2737 if (real(start_x) <= rx .and. rx <= real(end_x) .and. &
2738 real(start_y) <= ry .and. ry <= real(end_y)) then
2739 where_maps_to(min_i,max_j,1) = nint(rx)
2740 where_maps_to(min_i,max_j,2) = nint(ry)
2741 else
2742 where_maps_to(min_i,max_j,1) = OUTSIDE_DOMAIN
2743 end if
2744 end if
2746 ! Upper-right corner
2747 if (where_maps_to(max_i,max_j,1) == NOT_PROCESSED) then
2748 orig_selected_domain = iget_selected_domain()
2749 call select_domain(SOURCE_PROJ)
2750 call xytoll(real(max_i), real(max_j), lat_corner, lon_corner, istagger)
2751 call select_domain(1)
2752 call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
2753 rx = (rx - 1.0)*sr_x + 1.0
2754 ry = (ry - 1.0)*sr_y + 1.0
2755 call select_domain(orig_selected_domain)
2756 if (real(start_x) <= rx .and. rx <= real(end_x) .and. &
2757 real(start_y) <= ry .and. ry <= real(end_y)) then
2758 where_maps_to(max_i,max_j,1) = nint(rx)
2759 where_maps_to(max_i,max_j,2) = nint(ry)
2760 else
2761 where_maps_to(max_i,max_j,1) = OUTSIDE_DOMAIN
2762 end if
2763 end if
2765 ! Lower-right corner
2766 if (where_maps_to(max_i,min_j,1) == NOT_PROCESSED) then
2767 orig_selected_domain = iget_selected_domain()
2768 call select_domain(SOURCE_PROJ)
2769 call xytoll(real(max_i), real(min_j), lat_corner, lon_corner, istagger)
2770 call select_domain(1)
2771 call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
2772 rx = (rx - 1.0)*sr_x + 1.0
2773 ry = (ry - 1.0)*sr_y + 1.0
2774 call select_domain(orig_selected_domain)
2775 if (real(start_x) <= rx .and. rx <= real(end_x) .and. &
2776 real(start_y) <= ry .and. ry <= real(end_y)) then
2777 where_maps_to(max_i,min_j,1) = nint(rx)
2778 where_maps_to(max_i,min_j,2) = nint(ry)
2779 else
2780 where_maps_to(max_i,min_j,1) = OUTSIDE_DOMAIN
2781 end if
2782 end if
2784 ! If all four corners map to same model grid point, accumulate the
2785 ! entire rectangle
2786 if (where_maps_to(min_i,min_j,1) == where_maps_to(min_i,max_j,1) .and. &
2787 where_maps_to(min_i,min_j,1) == where_maps_to(max_i,max_j,1) .and. &
2788 where_maps_to(min_i,min_j,1) == where_maps_to(max_i,min_j,1) .and. &
2789 where_maps_to(min_i,min_j,2) == where_maps_to(min_i,max_j,2) .and. &
2790 where_maps_to(min_i,min_j,2) == where_maps_to(max_i,max_j,2) .and. &
2791 where_maps_to(min_i,min_j,2) == where_maps_to(max_i,min_j,2) .and. &
2792 where_maps_to(min_i,min_j,1) /= OUTSIDE_DOMAIN) then
2793 x_dest = where_maps_to(min_i,min_j,1)
2794 y_dest = where_maps_to(min_i,min_j,2)
2796 ! If this grid point was already given a value from higher-priority source data,
2797 ! there is nothing to do.
2798 ! if (.not. bitarray_test(processed_pts, x_dest-start_x+1, y_dest-start_y+1)) then
2800 ! If this grid point has never been given a value by this level of source data,
2801 ! initialize the point
2802 if (.not. bitarray_test(new_pts, x_dest-start_x+1, y_dest-start_y+1)) then
2803 do k=min_k,max_k
2804 dst_array(x_dest,y_dest,k) = 0.
2805 end do
2806 end if
2808 ! Sum all the points whose nearest neighbor is this grid point
2809 if (present(mask_array) .and. present(mask_relational)) then
2810 do i=min_i,max_i
2811 do j=min_j,max_j
2812 do k=min_k,max_k
2813 ! Ignore masked/missing values in the source data
2814 if (tile_array(i,j,k) /= msgval) then
2815 if (mask_relational == ' ' .and. mask_array(i,j) /= maskval) then
2816 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2817 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2818 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2819 else if (mask_relational == '<' .and. mask_array(i,j) >= maskval) then
2820 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2821 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2822 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2823 else if (mask_relational == '>' .and. mask_array(i,j) <= maskval) then
2824 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2825 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2826 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2827 end if
2828 end if
2829 end do
2830 end do
2831 end do
2832 else if (present(mask_array)) then
2833 do i=min_i,max_i
2834 do j=min_j,max_j
2835 do k=min_k,max_k
2836 ! Ignore masked/missing values in the source data
2837 if ((tile_array(i,j,k) /= msgval) .and. &
2838 (mask_array(i,j) /= maskval)) then
2839 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2840 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2841 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2842 end if
2843 end do
2844 end do
2845 end do
2846 else
2847 do i=min_i,max_i
2848 do j=min_j,max_j
2849 do k=min_k,max_k
2850 ! Ignore masked/missing values in the source data
2851 if ((tile_array(i,j,k) /= msgval)) then
2852 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2853 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2854 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2855 end if
2856 end do
2857 end do
2858 end do
2859 end if
2861 ! end if
2863 ! Rectangle is a square of four points, and we can simply deal with each of the points
2864 else if (((max_i - min_i + 1) <= 2) .and. ((max_j - min_j + 1) <= 2)) then
2865 do i=min_i,max_i
2866 do j=min_j,max_j
2867 x_dest = where_maps_to(i,j,1)
2868 y_dest = where_maps_to(i,j,2)
2870 if (x_dest /= OUTSIDE_DOMAIN) then
2872 ! if (.not. bitarray_test(processed_pts, x_dest-start_x+1, y_dest-start_y+1)) then
2873 if (.not. bitarray_test(new_pts, x_dest-start_x+1, y_dest-start_y+1)) then
2874 do k=min_k,max_k
2875 dst_array(x_dest,y_dest,k) = 0.
2876 end do
2877 end if
2879 if (present(mask_array) .and. present(mask_relational)) then
2880 do k=min_k,max_k
2881 ! Ignore masked/missing values
2882 if (tile_array(i,j,k) /= msgval) then
2883 if (mask_relational == ' ' .and. mask_array(i,j) /= maskval) then
2884 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2885 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2886 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2887 else if (mask_relational == '<' .and. mask_array(i,j) >= maskval) then
2888 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2889 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2890 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2891 else if (mask_relational == '>' .and. mask_array(i,j) <= maskval) then
2892 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2893 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2894 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2895 end if
2896 end if
2897 end do
2898 else if (present(mask_array)) then
2899 do k=min_k,max_k
2900 ! Ignore masked/missing values
2901 if ((tile_array(i,j,k) /= msgval) .and. &
2902 (mask_array(i,j) /= maskval)) then
2903 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2904 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2905 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2906 end if
2907 end do
2908 else
2909 do k=min_k,max_k
2910 ! Ignore masked/missing values
2911 if ((tile_array(i,j,k) /= msgval)) then
2912 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2913 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2914 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2915 end if
2916 end do
2917 end if
2918 ! end if
2920 end if
2921 end do
2922 end do
2924 ! Not all corners map to the same grid point, and the rectangle contains more than
2925 ! four points
2926 else
2927 center_i = (max_i + min_i)/2
2928 center_j = (max_j + min_j)/2
2930 ! Recursively process lower-left rectangle
2931 call process_continuous_block(tile_array, where_maps_to, &
2932 src_min_x, src_min_y, src_min_z, &
2933 src_max_x, src_max_y, src_max_z, &
2934 min_i, min_j, min_k, &
2935 center_i, center_j, max_k, &
2936 dst_array, n, &
2937 start_x, end_x, start_y, end_y, start_z, end_z, &
2938 istagger, &
2939 new_pts, sr_x, sr_y, msgval, maskval, mask_relational, mask_array)
2941 if (center_i < max_i) then
2942 ! Recursively process lower-right rectangle
2943 call process_continuous_block(tile_array, where_maps_to, &
2944 src_min_x, src_min_y, src_min_z, &
2945 src_max_x, src_max_y, src_max_z, &
2946 center_i+1, min_j, min_k, max_i, &
2947 center_j, max_k, &
2948 dst_array, n, &
2949 start_x, end_x, start_y, &
2950 end_y, start_z, end_z, &
2951 istagger, &
2952 new_pts, sr_x, sr_y, msgval, maskval, mask_relational, mask_array)
2953 end if
2955 if (center_j < max_j) then
2956 ! Recursively process upper-left rectangle
2957 call process_continuous_block(tile_array, where_maps_to, &
2958 src_min_x, src_min_y, src_min_z, &
2959 src_max_x, src_max_y, src_max_z, &
2960 min_i, center_j+1, min_k, center_i, &
2961 max_j, max_k, &
2962 dst_array, n, &
2963 start_x, end_x, start_y, &
2964 end_y, start_z, end_z, &
2965 istagger, &
2966 new_pts, sr_x, sr_y, msgval, maskval, mask_relational, mask_array)
2967 end if
2969 if (center_i < max_i .and. center_j < max_j) then
2970 ! Recursively process upper-right rectangle
2971 call process_continuous_block(tile_array, where_maps_to, &
2972 src_min_x, src_min_y, src_min_z, &
2973 src_max_x, src_max_y, src_max_z, &
2974 center_i+1, center_j+1, min_k, max_i, &
2975 max_j, max_k, &
2976 dst_array, n, &
2977 start_x, end_x, start_y, &
2978 end_y, start_z, end_z, &
2979 istagger, &
2980 new_pts, sr_x, sr_y, msgval, maskval, mask_relational, mask_array)
2981 end if
2982 end if
2984 end subroutine process_continuous_block
2986 end module process_domain_module