| blob | 69fed556fabe99bdcea9ce735678740bbbe5e349 |
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_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
37 real :: cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
38 dom_dx, dom_dy
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 character (len=19) :: valid_date, temp_date
43 character (len=128) :: title, mminlu
44 character (len=128), pointer, dimension(:) :: output_flags, td_output_flags
46 ! Compute number of times that we will process
47 call geth_idts(end_date(n), start_date(n), idiff)
48 call mprintf((idiff < 0),ERROR,'Ending date is earlier than starting date in namelist for domain %i.', i1=n)
50 n_times = idiff / interval_seconds
52 ! Check that the interval evenly divides the range of times to process
53 call mprintf((mod(idiff, interval_seconds) /= 0),WARN, &
54 'In namelist, interval_seconds does not evenly divide '// &
55 '(end_date - start_date) for domain %i. Only %i time periods '// &
56 'will be processed.', i1=n, i2=n_times)
58 ! Initialize the storage module
59 call mprintf(.true.,LOGFILE,'Initializing storage module')
60 call storage_init()
62 !
63 ! Do time-independent processing
64 !
65 call get_static_fields(n, dyn_opt, west_east_dim, south_north_dim, bottom_top_dim, map_proj, &
66 we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
67 we_patch_s, we_patch_e, &
68 we_patch_stag_s, we_patch_stag_e, &
69 sn_patch_s, sn_patch_e, &
70 sn_patch_stag_s, sn_patch_stag_e, &
71 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
72 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
73 mminlu, is_water, is_ice, is_urban, i_soilwater, &
74 grid_id, parent_id, &
75 i_parent_start, j_parent_start, i_parent_end, j_parent_end, &
76 parent_grid_ratio, cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
77 dom_dx, dom_dy, landmask, xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v, corner_lats, &
78 corner_lons, title)
80 allocate(output_flags(num_entries))
81 do i=1,num_entries
82 output_flags(i) = ' '
83 end do
85 ! This call is to process the constant met fields (SST or SEAICE, for example)
86 ! That we process constant fields is indicated by the first argument
87 call process_single_met_time(.true., temp_date, n, extra_row, extra_col, xlat, xlon, &
88 xlat_u, xlon_u, xlat_v, xlon_v, landmask, &
89 title, dyn_opt, &
90 west_east_dim, south_north_dim, &
91 we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
92 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
93 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
94 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
95 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
96 map_proj, mminlu, is_water, is_ice, is_urban, i_soilwater, &
97 grid_id, parent_id, i_parent_start, &
98 j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
99 cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, &
100 truelat2, parent_grid_ratio, corner_lats, corner_lons, output_flags)
102 !
103 ! Begin time-dependent processing
104 !
106 ! Loop over all times to be processed for this domain
107 do t=0,n_times
109 call geth_newdate(valid_date, trim(start_date(n)), t*interval_seconds)
110 temp_date = ' '
112 if (mod(interval_seconds,3600) == 0) then
113 write(temp_date,'(a13)') valid_date(1:10)//'_'//valid_date(12:13)
114 else if (mod(interval_seconds,60) == 0) then
115 write(temp_date,'(a16)') valid_date(1:10)//'_'//valid_date(12:16)
116 else
117 write(temp_date,'(a19)') valid_date(1:10)//'_'//valid_date(12:19)
118 end if
120 call mprintf(.true.,STDOUT, ' Processing %s', s1=trim(temp_date))
121 call mprintf(.true.,LOGFILE, 'Preparing to process output time %s', s1=temp_date)
123 allocate(td_output_flags(num_entries))
124 do i=1,num_entries
125 td_output_flags(i) = output_flags(i)
126 end do
128 call process_single_met_time(.false., temp_date, n, extra_row, extra_col, xlat, xlon, &
129 xlat_u, xlon_u, xlat_v, xlon_v, landmask, &
130 title, dyn_opt, &
131 west_east_dim, south_north_dim, &
132 we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
133 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
134 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
135 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
136 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
137 map_proj, mminlu, is_water, is_ice, is_urban, i_soilwater, &
138 grid_id, parent_id, i_parent_start, &
139 j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
140 cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, &
141 truelat2, parent_grid_ratio, corner_lats, corner_lons, td_output_flags)
143 deallocate(td_output_flags)
145 end do ! Loop over n_times
147 deallocate(output_flags)
149 call storage_delete_all()
151 end subroutine process_domain
154 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
155 ! Name: get_static_fields
156 !
157 ! Purpose:
158 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
159 subroutine get_static_fields(n, dyn_opt, west_east_dim, south_north_dim, bottom_top_dim, &
160 map_proj, &
161 we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
162 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
163 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
164 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
165 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
166 mminlu, is_water, is_ice, is_urban, i_soilwater, grid_id, parent_id, &
167 i_parent_start, j_parent_start, i_parent_end, j_parent_end, &
168 parent_grid_ratio, cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
169 dom_dx, dom_dy, landmask, xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v, corner_lats, &
170 corner_lons, title)
172 use gridinfo_module
173 use input_module
174 use llxy_module
175 use parallel_module
176 use storage_module
178 implicit none
180 ! Arguments
181 integer, intent(in) :: n
182 integer, intent(inout) :: dyn_opt, west_east_dim, south_north_dim, bottom_top_dim, &
183 map_proj, &
184 we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
185 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
186 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
187 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
188 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
189 is_water, is_ice, is_urban, i_soilwater, grid_id, parent_id, &
190 i_parent_start, j_parent_start, i_parent_end, j_parent_end, &
191 parent_grid_ratio
192 real, pointer, dimension(:,:) :: landmask
193 real, intent(inout) :: cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
194 dom_dx, dom_dy
195 real, pointer, dimension(:,:) :: xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v
196 real, dimension(16), intent(out) :: corner_lats, corner_lons
197 character (len=128), intent(inout) :: title, mminlu
199 ! Local variables
200 integer :: istatus, i, j, k, sp1, ep1, sp2, ep2, sp3, ep3, &
201 lh_mult, rh_mult, bh_mult, th_mult
202 real, pointer, dimension(:,:,:) :: real_array
203 character (len=3) :: memorder
204 character (len=128) :: grid_type, datestr, cname, stagger, cunits, cdesc
205 character (len=128), dimension(3) :: dimnames
206 type (fg_input) :: field
208 ! Initialize the input module to read static input data for this domain
209 call mprintf(.true.,LOGFILE,'Opening static input file.')
210 call input_init(n, istatus)
211 call mprintf((istatus /= 0),ERROR, 'input_init(): Error opening input for domain %i.', i1=n)
213 ! Read global attributes from the static data input file
214 call mprintf(.true.,LOGFILE,'Reading static global attributes.')
215 call read_global_attrs(title, datestr, grid_type, dyn_opt, west_east_dim, &
216 south_north_dim, bottom_top_dim, &
217 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
218 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
219 map_proj, mminlu, is_water, is_ice, is_urban, i_soilwater, &
220 grid_id, parent_id, i_parent_start, &
221 j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
222 cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, &
223 truelat2, parent_grid_ratio, corner_lats, corner_lons)
225 we_dom_s = 1
226 sn_dom_s = 1
227 if (grid_type(1:1) == 'C') then
228 we_dom_e = west_east_dim - 1
229 sn_dom_e = south_north_dim - 1
230 else if (grid_type(1:1) == 'E') then
231 we_dom_e = west_east_dim
232 sn_dom_e = south_north_dim
233 end if
235 ! Given the full dimensions of this domain, find out the range of indices
236 ! that will be worked on by this processor. This information is given by
237 ! my_minx, my_miny, my_maxx, my_maxy
238 call parallel_get_tile_dims(west_east_dim, south_north_dim)
240 ! Must figure out patch dimensions from info in parallel module
241 if (nprocs > 1 .and. .not. do_tiled_input) then
243 we_patch_s = my_minx
244 we_patch_stag_s = my_minx
245 we_patch_e = my_maxx - 1
246 sn_patch_s = my_miny
247 sn_patch_stag_s = my_miny
248 sn_patch_e = my_maxy - 1
250 if (gridtype == 'C') then
251 if (my_x /= nproc_x - 1) then
252 we_patch_e = we_patch_e + 1
253 we_patch_stag_e = we_patch_e
254 else
255 we_patch_stag_e = we_patch_e + 1
256 end if
257 if (my_y /= nproc_y - 1) then
258 sn_patch_e = sn_patch_e + 1
259 sn_patch_stag_e = sn_patch_e
260 else
261 sn_patch_stag_e = sn_patch_e + 1
262 end if
263 else if (gridtype == 'E') then
264 we_patch_e = we_patch_e + 1
265 sn_patch_e = sn_patch_e + 1
266 we_patch_stag_e = we_patch_e
267 sn_patch_stag_e = sn_patch_e
268 end if
270 end if
272 ! Compute multipliers for halo width; these must be 0/1
273 if (my_x /= 0) then
274 lh_mult = 1
275 else
276 lh_mult = 0
277 end if
278 if (my_x /= (nproc_x-1)) then
279 rh_mult = 1
280 else
281 rh_mult = 0
282 end if
283 if (my_y /= 0) then
284 bh_mult = 1
285 else
286 bh_mult = 0
287 end if
288 if (my_y /= (nproc_y-1)) then
289 th_mult = 1
290 else
291 th_mult = 0
292 end if
294 we_mem_s = we_patch_s - HALO_WIDTH*lh_mult
295 we_mem_e = we_patch_e + HALO_WIDTH*rh_mult
296 sn_mem_s = sn_patch_s - HALO_WIDTH*bh_mult
297 sn_mem_e = sn_patch_e + HALO_WIDTH*th_mult
298 we_mem_stag_s = we_patch_stag_s - HALO_WIDTH*lh_mult
299 we_mem_stag_e = we_patch_stag_e + HALO_WIDTH*rh_mult
300 sn_mem_stag_s = sn_patch_stag_s - HALO_WIDTH*bh_mult
301 sn_mem_stag_e = sn_patch_stag_e + HALO_WIDTH*th_mult
303 ! Initialize a proj_info type for the destination grid projection. This will
304 ! primarily be used for rotating Earth-relative winds to grid-relative winds
305 call set_domain_projection(map_proj, stand_lon, truelat1, truelat2, &
306 dom_dx, dom_dy, dom_dx, dom_dy, west_east_dim, &
307 south_north_dim, real(west_east_dim)/2., &
308 real(south_north_dim)/2.,cen_lat, cen_lon, &
309 cen_lat, cen_lon)
311 ! Read static fields using the input module; we know that there are no more
312 ! fields to be read when read_next_field() returns a non-zero status.
313 istatus = 0
314 do while (istatus == 0)
315 call read_next_field(sp1, ep1, sp2, ep2, sp3, ep3, cname, cunits, cdesc, &
316 memorder, stagger, dimnames, real_array, istatus)
317 if (istatus == 0) then
319 call mprintf(.true.,LOGFILE, 'Read in static field %s.',s1=cname)
321 ! We will also keep copies in core of the lat/lon arrays, for use in
322 ! interpolation of the met fields to the model grid.
323 ! For now, we assume that the lat/lon arrays will have known field names
324 if (index(cname, 'XLAT_M') /= 0 .and. &
325 len_trim(cname) == len_trim('XLAT_M')) then
326 allocate(xlat(we_mem_s:we_mem_e,sn_mem_s:sn_mem_e))
327 xlat(we_patch_s:we_patch_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
328 call exchange_halo_r(xlat, &
329 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, 1, 1, &
330 we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, 1, 1)
332 else if (index(cname, 'XLONG_M') /= 0 .and. &
333 len_trim(cname) == len_trim('XLONG_M')) then
334 allocate(xlon(we_mem_s:we_mem_e,sn_mem_s:sn_mem_e))
335 xlon(we_patch_s:we_patch_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
336 call exchange_halo_r(xlon, &
337 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, 1, 1, &
338 we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, 1, 1)
340 else if (index(cname, 'XLAT_U') /= 0 .and. &
341 len_trim(cname) == len_trim('XLAT_U')) then
342 allocate(xlat_u(we_mem_stag_s:we_mem_stag_e,sn_mem_s:sn_mem_e))
343 xlat_u(we_patch_stag_s:we_patch_stag_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
344 call exchange_halo_r(xlat_u, &
345 we_mem_stag_s, we_mem_stag_e, sn_mem_s, sn_mem_e, 1, 1, &
346 we_patch_stag_s, we_patch_stag_e, sn_patch_s, sn_patch_e, 1, 1)
348 else if (index(cname, 'XLONG_U') /= 0 .and. &
349 len_trim(cname) == len_trim('XLONG_U')) then
350 allocate(xlon_u(we_mem_stag_s:we_mem_stag_e,sn_mem_s:sn_mem_e))
351 xlon_u(we_patch_stag_s:we_patch_stag_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
352 call exchange_halo_r(xlon_u, &
353 we_mem_stag_s, we_mem_stag_e, sn_mem_s, sn_mem_e, 1, 1, &
354 we_patch_stag_s, we_patch_stag_e, sn_patch_s, sn_patch_e, 1, 1)
356 else if (index(cname, 'XLAT_V') /= 0 .and. &
357 len_trim(cname) == len_trim('XLAT_V')) then
358 allocate(xlat_v(we_mem_s:we_mem_e,sn_mem_stag_s:sn_mem_stag_e))
359 xlat_v(we_patch_s:we_patch_e,sn_patch_stag_s:sn_patch_stag_e) = real_array(:,:,1)
360 call exchange_halo_r(xlat_v, &
361 we_mem_s, we_mem_e, sn_mem_stag_s, sn_mem_stag_e, 1, 1, &
362 we_patch_s, we_patch_e, sn_patch_stag_s, sn_patch_stag_e, 1, 1)
364 else if (index(cname, 'XLONG_V') /= 0 .and. &
365 len_trim(cname) == len_trim('XLONG_V')) then
366 allocate(xlon_v(we_mem_s:we_mem_e,sn_mem_stag_s:sn_mem_stag_e))
367 xlon_v(we_patch_s:we_patch_e,sn_patch_stag_s:sn_patch_stag_e) = real_array(:,:,1)
368 call exchange_halo_r(xlon_v, &
369 we_mem_s, we_mem_e, sn_mem_stag_s, sn_mem_stag_e, 1, 1, &
370 we_patch_s, we_patch_e, sn_patch_stag_s, sn_patch_stag_e, 1, 1)
372 else if (index(cname, 'LANDMASK') /= 0 .and. &
373 len_trim(cname) == len_trim('LANDMASK')) then
374 allocate(landmask(we_mem_s:we_mem_e,sn_mem_s:sn_mem_e))
375 landmask(we_patch_s:we_patch_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
376 call exchange_halo_r(landmask, &
377 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, 1, 1, &
378 we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, 1, 1)
380 end if
382 ! Having read in a field, we write each level individually to the
383 ! storage module; levels will be reassembled later on when they
384 ! are written.
385 do k=sp3,ep3
386 field%header%version = 1
387 field%header%date = start_date(n)
388 field%header%time_dependent = .false.
389 field%header%mask_field = .false.
390 field%header%forecast_hour = 0.0
391 field%header%fg_source = 'geogrid_model'
392 field%header%field = cname
393 field%header%units = cunits
394 field%header%description = cdesc
395 field%header%vertical_coord = dimnames(3)
396 field%header%vertical_level = k
397 field%header%array_order = memorder
398 field%header%is_wind_grid_rel = .true.
399 field%header%array_has_missing_values = .false.
400 if (gridtype == 'C') then
401 if (trim(stagger) == 'M') then
402 field%map%stagger = M
403 field%header%dim1(1) = we_mem_s
404 field%header%dim1(2) = we_mem_e
405 field%header%dim2(1) = sn_mem_s
406 field%header%dim2(2) = sn_mem_e
407 else if (trim(stagger) == 'U') then
408 field%map%stagger = U
409 field%header%dim1(1) = we_mem_stag_s
410 field%header%dim1(2) = we_mem_stag_e
411 field%header%dim2(1) = sn_mem_s
412 field%header%dim2(2) = sn_mem_e
413 else if (trim(stagger) == 'V') then
414 field%map%stagger = V
415 field%header%dim1(1) = we_mem_s
416 field%header%dim1(2) = we_mem_e
417 field%header%dim2(1) = sn_mem_stag_s
418 field%header%dim2(2) = sn_mem_stag_e
419 end if
420 else if (gridtype == 'E') then
421 if (trim(stagger) == 'M') then
422 field%map%stagger = HH
423 else if (trim(stagger) == 'V') then
424 field%map%stagger = VV
425 end if
426 field%header%dim1(1) = we_mem_s
427 field%header%dim1(2) = we_mem_e
428 field%header%dim2(1) = sn_mem_s
429 field%header%dim2(2) = sn_mem_e
430 end if
432 allocate(field%valid_mask)
433 if (field%map%stagger == M .or. &
434 field%map%stagger == HH .or. &
435 field%map%stagger == VV) then
436 allocate(field%r_arr(we_mem_s:we_mem_e,&
437 sn_mem_s:sn_mem_e))
438 field%r_arr(we_patch_s:we_patch_e,sn_patch_s:sn_patch_e) = real_array(sp1:ep1,sp2:ep2,k)
439 call exchange_halo_r(field%r_arr, &
440 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, 1, 1, &
441 we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, 1, 1)
442 call bitarray_create(field%valid_mask, &
443 (we_mem_e-we_mem_s)+1, &
444 (sn_mem_e-sn_mem_s)+1)
445 do j=1,(sn_mem_e-sn_mem_s)+1
446 do i=1,(we_mem_e-we_mem_s)+1
447 call bitarray_set(field%valid_mask, i, j)
448 end do
449 end do
450 else if (field%map%stagger == U) then
451 allocate(field%r_arr(we_mem_stag_s:we_mem_stag_e,&
452 sn_mem_s:sn_mem_e))
453 field%r_arr(we_patch_stag_s:we_patch_stag_e,sn_patch_s:sn_patch_e) = real_array(sp1:ep1,sp2:ep2,k)
454 call exchange_halo_r(field%r_arr, &
455 we_mem_stag_s, we_mem_stag_e, sn_mem_s, sn_mem_e, 1, 1, &
456 we_patch_stag_s, we_patch_stag_e, sn_patch_s, sn_patch_e, 1, 1)
457 call bitarray_create(field%valid_mask, &
458 (we_mem_stag_e-we_mem_stag_s)+1, &
459 (sn_mem_e-sn_mem_s)+1)
460 do j=1,(sn_mem_e-sn_mem_s)+1
461 do i=1,(we_mem_stag_e-we_mem_stag_s)+1
462 call bitarray_set(field%valid_mask, i, j)
463 end do
464 end do
465 else if (field%map%stagger == V) then
466 allocate(field%r_arr(we_mem_s:we_mem_e,&
467 sn_mem_stag_s:sn_mem_stag_e))
468 field%r_arr(we_patch_s:we_patch_e,sn_patch_stag_s:sn_patch_stag_e) = real_array(sp1:ep1,sp2:ep2,k)
469 call exchange_halo_r(field%r_arr, &
470 we_mem_s, we_mem_e, sn_mem_stag_s, sn_mem_stag_e, 1, 1, &
471 we_patch_s, we_patch_e, sn_patch_stag_s, sn_patch_stag_e, 1, 1)
472 call bitarray_create(field%valid_mask, &
473 (we_mem_e-we_mem_s)+1, &
474 (sn_mem_stag_e-sn_mem_stag_s)+1)
475 do j=1,(sn_mem_stag_e-sn_mem_stag_s)+1
476 do i=1,(we_mem_e-we_mem_s)+1
477 call bitarray_set(field%valid_mask, i, j)
478 end do
479 end do
480 end if
482 nullify(field%modified_mask)
484 call storage_put_field(field)
486 end do
488 end if
489 end do
491 ! Done reading all static fields for this domain
492 call input_close()
494 end subroutine get_static_fields
497 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
498 ! Name: process_single_met_time
499 !
500 ! Purpose:
501 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
502 subroutine process_single_met_time(do_const_processing, &
503 temp_date, n, extra_row, extra_col, xlat, xlon, &
504 xlat_u, xlon_u, xlat_v, xlon_v, landmask, &
505 title, dyn_opt, &
506 west_east_dim, south_north_dim, &
507 we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
508 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
509 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
510 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
511 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
512 map_proj, mminlu, is_water, is_ice, is_urban, i_soilwater, &
513 grid_id, parent_id, i_parent_start, &
514 j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
515 cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, &
516 truelat2, parent_grid_ratio, corner_lats, corner_lons, output_flags)
518 use bitarray_module
519 use gridinfo_module
520 use interp_module
521 use interp_option_module
522 use llxy_module
523 use misc_definitions_module
524 use module_debug
525 use output_module
526 use parallel_module
527 use read_met_module
528 use rotate_winds_module
529 use storage_module
531 implicit none
533 ! Arguments
534 logical, intent(in) :: do_const_processing
535 integer, intent(in) :: n, dyn_opt, west_east_dim, south_north_dim, map_proj, &
536 we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
537 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
538 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
539 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
540 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
541 is_water, is_ice, is_urban, i_soilwater, &
542 grid_id, parent_id, i_parent_start, j_parent_start, &
543 i_parent_end, j_parent_end, parent_grid_ratio
544 ! BUG: Should we be passing these around as pointers, or just declare them as arrays?
545 real, pointer, dimension(:,:) :: landmask
546 real, intent(in) :: dom_dx, dom_dy, cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2
547 real, dimension(16), intent(in) :: corner_lats, corner_lons
548 real, pointer, dimension(:,:) :: xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v
549 logical, intent(in) :: extra_row, extra_col
550 character (len=19), intent(in) :: temp_date
551 character (len=128), intent(in) :: mminlu
552 character (len=128), pointer, dimension(:) :: output_flags
554 ! BUG: Move this constant to misc_definitions_module?
555 integer, parameter :: BDR_WIDTH = 3
557 ! Local variables
558 integer :: istatus, iqstatus, fg_idx, idx, idxt, i, j, bottom_top_dim, &
559 sm1, em1, sm2, em2, sm3, em3, &
560 sp1, ep1, sp2, ep2, sp3, ep3, &
561 sd1, ed1, sd2, ed2, sd3, ed3, met_map_proj, &
562 version, nx, ny, u_idx, bdr_wdth
563 real :: rx, ry, xfcst, xlvl, startlat, startlon, starti, startj, deltalat, deltalon, earth_radius
564 real :: threshold, met_dx, met_dy
565 real :: met_cen_lon, met_truelat1, met_truelat2
566 logical :: is_wind_grid_rel, do_gcell_interp
567 integer, pointer, dimension(:) :: u_levels, v_levels
568 real, pointer, dimension(:,:) :: slab, halo_slab
569 real, pointer, dimension(:,:,:) :: real_array
570 logical, pointer, dimension(:) :: got_this_field
571 character (len=9) :: short_fieldnm
572 character (len=19) :: output_date
573 character (len=24) :: hdate
574 character (len=25) :: units
575 character (len=32) :: map_src
576 character (len=46) :: desc
577 character (len=128) :: cname, title, input_name
578 type (fg_input) :: field, u_field, v_field
580 allocate(got_this_field(num_entries))
581 do i=1,num_entries
582 got_this_field(i) = .false.
583 end do
585 ! For this time, we need to process all first-guess filename roots. When we
586 ! hit a root containing a '*', we assume we have hit the end of the list
587 fg_idx = 1
588 if (do_const_processing) then
589 input_name = constants_name(fg_idx)
590 else
591 input_name = fg_name(fg_idx)
592 end if
593 do while (input_name /= '*')
595 call mprintf(.true.,STDOUT, ' %s', s1=input_name)
596 call mprintf(.true.,LOGFILE, 'Getting input fields from %s', s1=input_name)
598 ! Do a first pass through this fg source to get all mask fields used
599 ! during interpolation
600 call get_interp_masks(trim(input_name), do_const_processing, temp_date)
602 istatus = 0
604 ! Initialize the module for reading in the met fields
605 call read_met_init(trim(input_name), do_const_processing, temp_date, istatus)
607 if (istatus == 0) then
609 ! Process all fields and levels from the current file; read_next_met_field()
610 ! will return a non-zero status when there are no more fields to be read.
611 do while (istatus == 0)
613 call read_next_met_field(version, short_fieldnm, hdate, xfcst, xlvl, units, desc, &
614 met_map_proj, startlat, startlon, starti, startj, deltalat, &
615 deltalon, met_dx, met_dy, met_cen_lon, &
616 met_truelat1, met_truelat2, earth_radius, nx, ny, &
617 map_src, slab, is_wind_grid_rel, istatus)
619 if (istatus == 0) then
621 ! Find index into fieldname, interp_method, masked, and fill_missing
622 ! of the current field
623 idxt = num_entries + 1
624 do idx=1,num_entries
625 if ((index(fieldname(idx), trim(short_fieldnm)) /= 0) .and. &
626 (len_trim(fieldname(idx)) == len_trim(short_fieldnm))) then
628 got_this_field(idx) = .true.
630 if (index(input_name,trim(from_input(idx))) /= 0 .or. &
631 (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
632 idxt = idx
633 end if
635 end if
636 end do
637 idx = idxt
638 if (idx > num_entries) idx = num_entries ! The last entry is a default
640 ! Do we need to rename this field?
641 if (output_name(idx) /= ' ') then
642 short_fieldnm = output_name(idx)(1:9)
644 idxt = num_entries + 1
645 do idx=1,num_entries
646 if ((index(fieldname(idx), trim(short_fieldnm)) /= 0) .and. &
647 (len_trim(fieldname(idx)) == len_trim(short_fieldnm))) then
649 got_this_field(idx) = .true.
651 if (index(input_name,trim(from_input(idx))) /= 0 .or. &
652 (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
653 idxt = idx
654 end if
656 end if
657 end do
658 idx = idxt
659 if (idx > num_entries) idx = num_entries ! The last entry is a default
660 end if
662 ! Do a simple check to see whether this is a global lat/lon dataset
663 if (met_map_proj == PROJ_LATLON .and. &
664 nx * deltalon == 360.) then
665 bdr_wdth = BDR_WIDTH
666 allocate(halo_slab(1-BDR_WIDTH:nx+BDR_WIDTH,1:ny))
667 halo_slab(1:nx, 1:ny) = slab(1:nx, 1:ny)
668 halo_slab(1-BDR_WIDTH:0, 1:ny) = slab(nx-BDR_WIDTH+1:nx, 1:ny)
669 halo_slab(nx+1:nx+BDR_WIDTH, 1:ny) = slab(1:BDR_WIDTH, 1:ny)
670 deallocate(slab)
671 else
672 bdr_wdth = 0
673 halo_slab => slab
674 nullify(slab)
675 end if
677 call mprintf(.true.,LOGFILE,'Processing %s at level %f.',s1=short_fieldnm,f1=xlvl)
679 call push_source_projection(met_map_proj, met_cen_lon, met_truelat1, &
680 met_truelat2, met_dx, met_dy, deltalat, deltalon, starti, startj, &
681 startlat, startlon, earth_radius=earth_radius*1000.)
683 ! Initialize fg_input structure to store the field
684 field%header%version = 1
685 field%header%date = hdate//' '
686 if (do_const_processing) then
687 field%header%time_dependent = .false.
688 else
689 field%header%time_dependent = .true.
690 end if
691 field%header%forecast_hour = xfcst
692 field%header%fg_source = 'FG'
693 field%header%field = ' '
694 field%header%field(1:9) = short_fieldnm
695 field%header%units = ' '
696 field%header%units(1:25) = units
697 field%header%description = ' '
698 field%header%description(1:46) = desc
699 call get_z_dim_name(short_fieldnm,field%header%vertical_coord)
700 field%header%vertical_level = nint(xlvl)
701 field%header%array_order = 'XY '
702 field%header%is_wind_grid_rel = is_wind_grid_rel
703 field%header%array_has_missing_values = .false.
704 nullify(field%r_arr)
705 nullify(field%valid_mask)
706 nullify(field%modified_mask)
708 if (output_this_field(idx) .and. flag_in_output(idx) /= ' ') then
709 output_flags(idx) = flag_in_output(idx)
710 end if
712 ! If we should not output this field, just list it as a mask field
713 if (output_this_field(idx)) then
714 field%header%mask_field = .false.
715 else
716 field%header%mask_field = .true.
717 end if
719 !
720 ! Before actually doing any interpolation to the model grid, we must check
721 ! whether we will be using the average_gcell interpolator that averages all
722 ! source points in each model grid cell
723 !
724 do_gcell_interp = .false.
725 if (index(interp_method(idx),'average_gcell') /= 0) then
727 call get_gcell_threshold(interp_method(idx), threshold, istatus)
728 if (istatus == 0) then
729 if (met_dx == 0. .and. met_dy == 0. .and. &
730 deltalat /= 0. .and. deltalon /= 0.) then
731 met_dx = abs(deltalon)
732 met_dy = abs(deltalat)
733 else
734 ! BUG: Need to more correctly handle dx/dy in meters.
735 met_dx = met_dx / 111000. ! Convert meters to approximate degrees
736 met_dy = met_dy / 111000.
737 end if
738 if (gridtype == 'C') then
739 if (threshold*max(met_dx,met_dy)*111. <= max(dom_dx,dom_dy)/1000.) &
740 do_gcell_interp = .true.
741 else if (gridtype == 'E') then
742 if (threshold*max(met_dx,met_dy) <= max(dom_dx,dom_dy)) &
743 do_gcell_interp = .true.
744 end if
745 end if
746 end if
748 ! Interpolate to U staggering
749 if (output_stagger(idx) == U) then
751 call storage_query_field(field, iqstatus)
752 if (iqstatus == 0) then
753 call storage_get_field(field, iqstatus)
754 call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// &
755 ' but could not get data.',s1=short_fieldnm,i1=nint(xlvl))
756 if (associated(field%modified_mask)) then
757 call bitarray_destroy(field%modified_mask)
758 nullify(field%modified_mask)
759 end if
760 else
761 allocate(field%valid_mask)
762 call bitarray_create(field%valid_mask, we_mem_stag_e-we_mem_stag_s+1, sn_mem_e-sn_mem_s+1)
763 end if
765 ! Save a copy of the fg_input structure for the U field so that we can find it later
766 if (is_u_field(idx)) call dup(field, u_field)
768 allocate(field%modified_mask)
769 call bitarray_create(field%modified_mask, we_mem_stag_e-we_mem_stag_s+1, sn_mem_e-sn_mem_s+1)
771 if (do_const_processing .or. field%header%time_dependent) then
772 call interp_met_field(input_name, short_fieldnm, U, M, &
773 field, xlat_u, xlon_u, we_mem_stag_s, we_mem_stag_e, sn_mem_s, sn_mem_e, &
774 halo_slab, 1-bdr_wdth, nx+bdr_wdth, 1, ny, bdr_wdth, do_gcell_interp, &
775 field%modified_mask)
776 else
777 call mprintf(.true.,INFORM,' - already processed this field from constant file.')
778 end if
780 ! Interpolate to V staggering
781 else if (output_stagger(idx) == V) then
783 call storage_query_field(field, iqstatus)
784 if (iqstatus == 0) then
785 call storage_get_field(field, iqstatus)
786 call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// &
787 ' but could not get data.',s1=short_fieldnm,i1=nint(xlvl))
788 if (associated(field%modified_mask)) then
789 call bitarray_destroy(field%modified_mask)
790 nullify(field%modified_mask)
791 end if
792 else
793 allocate(field%valid_mask)
794 call bitarray_create(field%valid_mask, we_mem_e-we_mem_s+1, sn_mem_stag_e-sn_mem_stag_s+1)
795 end if
797 ! Save a copy of the fg_input structure for the V field so that we can find it later
798 if (is_v_field(idx)) call dup(field, v_field)
800 allocate(field%modified_mask)
801 call bitarray_create(field%modified_mask, we_mem_e-we_mem_s+1, sn_mem_stag_e-sn_mem_stag_s+1)
803 if (do_const_processing .or. field%header%time_dependent) then
804 call interp_met_field(input_name, short_fieldnm, V, M, &
805 field, xlat_v, xlon_v, we_mem_s, we_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
806 halo_slab, 1-bdr_wdth, nx+bdr_wdth, 1, ny, bdr_wdth, do_gcell_interp, &
807 field%modified_mask)
808 else
809 call mprintf(.true.,INFORM,' - already processed this field from constant file.')
810 end if
812 ! Interpolate to VV staggering
813 else if (output_stagger(idx) == VV) then
815 call storage_query_field(field, iqstatus)
816 if (iqstatus == 0) then
817 call storage_get_field(field, iqstatus)
818 call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// &
819 ' but could not get data.',s1=short_fieldnm,i1=nint(xlvl))
820 if (associated(field%modified_mask)) then
821 call bitarray_destroy(field%modified_mask)
822 nullify(field%modified_mask)
823 end if
824 else
825 allocate(field%valid_mask)
826 call bitarray_create(field%valid_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1)
827 end if
829 ! Save a copy of the fg_input structure for the U field so that we can find it later
830 if (is_u_field(idx)) call dup(field, u_field)
832 ! Save a copy of the fg_input structure for the V field so that we can find it later
833 if (is_v_field(idx)) call dup(field, v_field)
835 allocate(field%modified_mask)
836 call bitarray_create(field%modified_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1)
838 if (do_const_processing .or. field%header%time_dependent) then
839 call interp_met_field(input_name, short_fieldnm, VV, M, &
840 field, xlat_v, xlon_v, we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
841 halo_slab, 1-bdr_wdth, nx+bdr_wdth, 1, ny, bdr_wdth, do_gcell_interp, &
842 field%modified_mask)
843 else
844 call mprintf(.true.,INFORM,' - already processed this field from constant file.')
845 end if
847 ! All other fields interpolated to M staggering for C grid, H staggering for E grid
848 else
850 call storage_query_field(field, iqstatus)
851 if (iqstatus == 0) then
852 call storage_get_field(field, iqstatus)
853 call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// &
854 ' but could not get data.',s1=short_fieldnm,i1=nint(xlvl))
855 if (associated(field%modified_mask)) then
856 call bitarray_destroy(field%modified_mask)
857 nullify(field%modified_mask)
858 end if
859 else
860 allocate(field%valid_mask)
861 call bitarray_create(field%valid_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1)
862 end if
864 allocate(field%modified_mask)
865 call bitarray_create(field%modified_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1)
867 if (do_const_processing .or. field%header%time_dependent) then
868 if (gridtype == 'C') then
869 call interp_met_field(input_name, short_fieldnm, M, M, &
870 field, xlat, xlon, we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
871 halo_slab, 1-bdr_wdth, nx+bdr_wdth, 1, ny, bdr_wdth, do_gcell_interp, &
872 field%modified_mask, landmask)
874 else if (gridtype == 'E') then
875 call interp_met_field(input_name, short_fieldnm, HH, M, &
876 field, xlat, xlon, we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
877 halo_slab, 1-bdr_wdth, nx+bdr_wdth, 1, ny, bdr_wdth, do_gcell_interp, &
878 field%modified_mask, landmask)
879 end if
880 else
881 call mprintf(.true.,INFORM,' - already processed this field from constant file.')
882 end if
884 end if
886 call bitarray_merge(field%valid_mask, field%modified_mask)
888 deallocate(halo_slab)
890 ! Store the interpolated field
891 call storage_put_field(field)
893 call pop_source_projection()
895 end if
896 end do
898 call read_met_close()
900 call push_source_projection(met_map_proj, met_cen_lon, met_truelat1, &
901 met_truelat2, met_dx, met_dy, deltalat, deltalon, starti, startj, &
902 startlat, startlon, earth_radius=earth_radius*1000.)
904 !
905 ! If necessary, rotate winds to earth-relative for this fg source
906 !
908 call storage_get_levels(u_field, u_levels)
909 call storage_get_levels(v_field, v_levels)
911 if (associated(u_levels) .and. associated(v_levels)) then
912 u_idx = 1
913 do u_idx = 1, size(u_levels)
914 u_field%header%vertical_level = u_levels(u_idx)
915 call storage_get_field(u_field, istatus)
916 v_field%header%vertical_level = v_levels(u_idx)
917 call storage_get_field(v_field, istatus)
919 if (associated(u_field%modified_mask) .and. &
920 associated(v_field%modified_mask)) then
922 ! BUG: Need to consider that E grid doesn't have u_s and u_e, since no U staggering.
923 ! BUG: Perhaps we need entirely separate map rotation routines for E staggering, where U and V are colocated.
924 if (u_field%header%is_wind_grid_rel) then
925 if (gridtype == 'C') then
926 call map_to_met(u_field%r_arr, u_field%modified_mask, &
927 v_field%r_arr, v_field%modified_mask, &
928 we_mem_stag_s, sn_mem_s, &
929 we_mem_stag_e, sn_mem_e, &
930 we_mem_s, sn_mem_stag_s, &
931 we_mem_e, sn_mem_stag_e, &
932 xlon_u, xlon_v, xlat_u, xlat_v)
933 else if (gridtype == 'E') then
934 call map_to_met_nmm(u_field%r_arr, u_field%modified_mask, &
935 v_field%r_arr, v_field%modified_mask, &
936 we_mem_s, sn_mem_s, &
937 we_mem_e, sn_mem_e, &
938 xlat_v, xlon_v)
939 end if
940 end if
942 call bitarray_destroy(u_field%modified_mask)
943 call bitarray_destroy(v_field%modified_mask)
944 nullify(u_field%modified_mask)
945 nullify(v_field%modified_mask)
946 call storage_put_field(u_field)
947 call storage_put_field(v_field)
948 end if
950 end do
952 deallocate(u_levels)
953 deallocate(v_levels)
955 end if
957 call pop_source_projection()
959 else
960 if (do_const_processing) then
961 call mprintf(.true.,WARN,'Couldn''t open file %s for input.',s1=input_name)
962 else
963 call mprintf(.true.,WARN,'Couldn''t open file %s for input.',s1=trim(input_name)//':'//trim(temp_date))
964 end if
965 end if
967 fg_idx = fg_idx + 1
968 if (do_const_processing) then
969 input_name = constants_name(fg_idx)
970 else
971 input_name = fg_name(fg_idx)
972 end if
973 end do ! while (input_name /= '*')
975 !
976 ! Rotate winds from earth-relative to grid-relative
977 !
979 call storage_get_levels(u_field, u_levels)
980 call storage_get_levels(v_field, v_levels)
982 if (associated(u_levels) .and. associated(v_levels)) then
983 u_idx = 1
984 do u_idx = 1, size(u_levels)
985 u_field%header%vertical_level = u_levels(u_idx)
986 call storage_get_field(u_field, istatus)
987 v_field%header%vertical_level = v_levels(u_idx)
988 call storage_get_field(v_field, istatus)
990 ! BUG: Need to consider that E grid doesn't have u_s and u_e, since no U staggering.
991 ! BUG: Perhaps we need entirely separate map rotation routines for E staggering, where U and V are colocated.
992 if (gridtype == 'C') then
993 call met_to_map(u_field%r_arr, u_field%valid_mask, &
994 v_field%r_arr, v_field%valid_mask, &
995 we_mem_stag_s, sn_mem_s, &
996 we_mem_stag_e, sn_mem_e, &
997 we_mem_s, sn_mem_stag_s, &
998 we_mem_e, sn_mem_stag_e, &
999 xlon_u, xlon_v, xlat_u, xlat_v)
1000 else if (gridtype == 'E') then
1001 call met_to_map_nmm(u_field%r_arr, u_field%valid_mask, &
1002 v_field%r_arr, v_field%valid_mask, &
1003 we_mem_s, sn_mem_s, &
1004 we_mem_e, sn_mem_e, &
1005 xlat_v, xlon_v)
1006 end if
1008 end do
1010 deallocate(u_levels)
1011 deallocate(v_levels)
1013 end if
1015 if (do_const_processing) return
1017 !
1018 ! Now that we have all degribbed fields, we build a 3-d pressure field, and fill in any
1019 ! missing levels in the other 3-d fields
1020 !
1021 call mprintf(.true.,LOGFILE,'Filling missing levels.')
1022 call fill_missing_levels(output_flags)
1024 call mprintf(.true.,LOGFILE,'Creating derived fields.')
1025 call create_derived_fields(gridtype, hdate, xfcst, &
1026 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
1027 we_mem_stag_s, we_mem_stag_e, sn_mem_stag_s, sn_mem_stag_e, &
1028 got_this_field, output_flags)
1030 !
1031 ! Check that every mandatory field was found in input data
1032 !
1033 do i=1,num_entries
1034 if (is_mandatory(i) .and. .not. got_this_field(i)) then
1035 call mprintf(.true.,ERROR,'The mandatory field %s was not found in any input data.',s1=fieldname(i))
1036 end if
1037 end do
1038 deallocate(got_this_field)
1040 !
1041 ! Before we begin to write fields, if debug_level is set high enough, we
1042 ! write a table of which fields are available at which levels to the
1043 ! metgrid.log file
1044 !
1045 ! call storage_print_fields()
1047 !
1048 ! All of the processing is now done for this time period for this domain;
1049 ! now we simply output every field from the storage module.
1050 !
1052 title = 'OUTPUT FROM METGRID'
1054 ! Initialize the output module for this domain and time
1055 call mprintf(.true.,LOGFILE,'Initializing output module.')
1056 output_date = temp_date
1057 if (len_trim(temp_date) == 13) then
1058 output_date(14:19) = ':00:00'
1059 else if (len_trim(temp_date) == 16) then
1060 output_date(17:19) = ':00'
1061 end if
1062 call output_init(n, title, output_date, gridtype, dyn_opt, &
1063 corner_lats, corner_lons, &
1064 we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
1065 we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, &
1066 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
1067 extra_col, extra_row)
1069 call get_bottom_top_dim(bottom_top_dim)
1071 ! First write out global attributes
1072 call mprintf(.true.,LOGFILE,'Writing global attributes to output.')
1073 call write_global_attrs(title, output_date, gridtype, dyn_opt, west_east_dim, &
1074 south_north_dim, bottom_top_dim, &
1075 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
1076 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
1077 map_proj, mminlu, is_water, is_ice, is_urban, i_soilwater, &
1078 grid_id, parent_id, i_parent_start, &
1079 j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
1080 cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, &
1081 truelat2, parent_grid_ratio, corner_lats, corner_lons, output_flags, num_entries)
1083 call reset_next_field()
1085 istatus = 0
1087 ! Now loop over all output fields, writing each to the output module
1088 do while (istatus == 0)
1089 call get_next_output_field(cname, real_array, &
1090 sm1, em1, sm2, em2, sm3, em3, istatus)
1091 if (istatus == 0) then
1093 call mprintf(.true.,LOGFILE,'Writing field %s to output.',s1=cname)
1094 call write_field(sm1, em1, sm2, em2, sm3, em3, &
1095 cname, output_date, real_array)
1096 deallocate(real_array)
1098 end if
1100 end do
1102 call mprintf(.true.,LOGFILE,'Closing output file.')
1103 call output_close()
1105 ! Free up memory used by met fields for this valid time
1106 call storage_delete_all_td()
1108 end subroutine process_single_met_time
1111 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1112 ! Name: get_interp_masks
1113 !
1114 ! Purpose:
1115 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1116 subroutine get_interp_masks(fg_prefix, is_constants, fg_date)
1118 use interp_option_module
1119 use read_met_module
1120 use storage_module
1122 implicit none
1124 ! Arguments
1125 logical, intent(in) :: is_constants
1126 character (len=*), intent(in) :: fg_prefix, fg_date
1128 ! BUG: Move this constant to misc_definitions_module?
1129 integer, parameter :: BDR_WIDTH = 3
1131 ! Local variables
1132 integer :: i, istatus, version, nx, ny, iproj, idx, idxt
1133 real :: xfcst, xlvl, startlat, startlon, starti, startj, &
1134 deltalat, deltalon, dx, dy, xlonc, truelat1, truelat2, &
1135 earth_radius
1136 real, pointer, dimension(:,:) :: slab
1137 logical :: is_wind_grid_rel
1138 character (len=9) :: field
1139 character (len=24) :: hdate
1140 character (len=25) :: units
1141 character (len=32) :: map_source
1142 character (len=46) :: desc
1143 type (fg_input) :: mask_field
1145 istatus = 0
1147 call read_met_init(fg_prefix, is_constants, fg_date, istatus)
1149 do while (istatus == 0)
1151 call read_next_met_field(version, field, hdate, xfcst, xlvl, units, desc, &
1152 iproj, startlat, startlon, starti, startj, deltalat, &
1153 deltalon, dx, dy, xlonc, truelat1, truelat2, &
1154 earth_radius, nx, ny, map_source, &
1155 slab, is_wind_grid_rel, istatus)
1157 if (istatus == 0) then
1159 ! Find out which METGRID.TBL entry goes with this field
1160 idxt = num_entries + 1
1161 do idx=1,num_entries
1162 if ((index(fieldname(idx), trim(field)) /= 0) .and. &
1163 (len_trim(fieldname(idx)) == len_trim(field))) then
1165 if (index(fg_prefix,trim(from_input(idx))) /= 0 .or. &
1166 (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
1167 idxt = idx
1168 end if
1170 end if
1171 end do
1172 idx = idxt
1173 if (idx > num_entries) idx = num_entries ! The last entry is a default
1175 ! Do we need to rename this field?
1176 if (output_name(idx) /= ' ') then
1177 field = output_name(idx)(1:9)
1179 idxt = num_entries + 1
1180 do idx=1,num_entries
1181 if ((index(fieldname(idx), trim(field)) /= 0) .and. &
1182 (len_trim(fieldname(idx)) == len_trim(field))) then
1184 if (index(fg_prefix,trim(from_input(idx))) /= 0 .or. &
1185 (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
1186 idxt = idx
1187 end if
1189 end if
1190 end do
1191 idx = idxt
1192 if (idx > num_entries) idx = num_entries ! The last entry is a default
1193 end if
1195 do i=1,num_entries
1196 if (interp_mask(i) /= ' ' .and. (trim(interp_mask(i)) == trim(field))) then
1198 mask_field%header%version = 1
1199 mask_field%header%date = ' '
1200 mask_field%header%date = fg_date
1201 if (is_constants) then
1202 mask_field%header%time_dependent = .false.
1203 else
1204 mask_field%header%time_dependent = .true.
1205 end if
1206 mask_field%header%mask_field = .true.
1207 mask_field%header%forecast_hour = 0.
1208 mask_field%header%fg_source = 'degribbed met data'
1209 mask_field%header%field = trim(field)//'.mask'
1210 mask_field%header%units = '-'
1211 mask_field%header%description = '-'
1212 mask_field%header%vertical_coord = 'none'
1213 mask_field%header%vertical_level = 1
1214 mask_field%header%array_order = 'XY'
1215 mask_field%header%dim1(1) = 1
1216 mask_field%header%dim1(2) = nx
1217 mask_field%header%dim2(1) = 1
1218 mask_field%header%dim2(2) = ny
1219 mask_field%header%is_wind_grid_rel = .true.
1220 mask_field%header%array_has_missing_values = .false.
1221 mask_field%map%stagger = M
1223 ! Do a simple check to see whether this is a global lat/lon dataset
1224 if (iproj == PROJ_LATLON .and. &
1225 nx * deltalon == 360.) then
1226 allocate(mask_field%r_arr(1-BDR_WIDTH:nx+BDR_WIDTH,1:ny))
1227 mask_field%r_arr(1:nx, 1:ny) = slab(1:nx, 1:ny)
1228 mask_field%r_arr(1-BDR_WIDTH:0, 1:ny) = slab(nx-BDR_WIDTH+1:nx, 1:ny)
1229 mask_field%r_arr(nx+1:nx+BDR_WIDTH, 1:ny) = slab(1:BDR_WIDTH, 1:ny)
1230 else
1231 allocate(mask_field%r_arr(1:nx,1:ny))
1232 mask_field%r_arr = slab
1233 end if
1235 nullify(mask_field%valid_mask)
1236 nullify(mask_field%modified_mask)
1238 call storage_put_field(mask_field)
1240 exit
1242 end if
1243 end do
1245 if (associated(slab)) deallocate(slab)
1247 end if
1249 end do
1251 call read_met_close()
1253 end subroutine get_interp_masks
1256 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1257 ! Name: interp_met_field
1258 !
1259 ! Purpose:
1260 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1261 subroutine interp_met_field(input_name, short_fieldnm, ifieldstagger, istagger, &
1262 field, xlat, xlon, sm1, em1, sm2, em2, &
1263 slab, minx, maxx, miny, maxy, bdr, do_gcell_interp, &
1264 new_pts, landmask)
1266 use bitarray_module
1267 use interp_module
1268 use interp_option_module
1269 use llxy_module
1270 use misc_definitions_module
1271 use storage_module
1273 implicit none
1275 ! Arguments
1276 integer, intent(in) :: ifieldstagger, istagger, &
1277 sm1, em1, sm2, em2, &
1278 minx, maxx, miny, maxy, bdr
1279 real, dimension(minx:maxx,miny:maxy), intent(in) :: slab
1280 real, dimension(sm1:em1,sm2:em2), intent(in) :: xlat, xlon
1281 real, dimension(sm1:em1,sm2:em2), intent(in), optional :: landmask
1282 logical, intent(in) :: do_gcell_interp
1283 character (len=9), intent(in) :: short_fieldnm
1284 character (len=128), intent(in) :: input_name
1285 type (fg_input), intent(inout) :: field
1286 type (bitarray), intent(inout) :: new_pts
1288 ! Local variables
1289 integer :: i, j, idx, idxt, orig_selected_proj, interp_mask_status, &
1290 interp_land_mask_status, interp_water_mask_status
1291 integer, pointer, dimension(:) :: interp_array
1292 real :: rx, ry, temp
1293 real, pointer, dimension(:,:) :: data_count
1294 type (fg_input) :: mask_field, mask_water_field, mask_land_field
1296 ! Find index into fieldname, interp_method, masked, and fill_missing
1297 ! of the current field
1298 idxt = num_entries + 1
1299 do idx=1,num_entries
1300 if ((index(fieldname(idx), trim(short_fieldnm)) /= 0) .and. &
1301 (len_trim(fieldname(idx)) == len_trim(short_fieldnm))) then
1302 if (index(input_name,trim(from_input(idx))) /= 0 .or. &
1303 (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
1304 idxt = idx
1305 end if
1306 end if
1307 end do
1308 idx = idxt
1309 if (idx > num_entries) then
1310 call mprintf(.true.,WARN,'Entry in METGRID.TBL not found for field %s. '// &
1311 'Default options will be used for this field!', s1=short_fieldnm)
1312 idx = num_entries ! The last entry is a default
1313 end if
1315 field%header%dim1(1) = sm1
1316 field%header%dim1(2) = em1
1317 field%header%dim2(1) = sm2
1318 field%header%dim2(2) = em2
1319 field%map%stagger = ifieldstagger
1320 if (.not. associated(field%r_arr)) then
1321 allocate(field%r_arr(sm1:em1,sm2:em2))
1322 end if
1324 interp_mask_status = 1
1325 interp_land_mask_status = 1
1326 interp_water_mask_status = 1
1328 if (interp_mask(idx) /= ' ') then
1329 mask_field%header%version = 1
1330 mask_field%header%forecast_hour = 0.
1331 mask_field%header%field = trim(interp_mask(idx))//'.mask'
1332 mask_field%header%vertical_coord = 'none'
1333 mask_field%header%vertical_level = 1
1335 call storage_get_field(mask_field, interp_mask_status)
1337 end if
1338 if (interp_land_mask(idx) /= ' ') then
1339 mask_land_field%header%version = 1
1340 mask_land_field%header%forecast_hour = 0.
1341 mask_land_field%header%field = trim(interp_land_mask(idx))//'.mask'
1342 mask_land_field%header%vertical_coord = 'none'
1343 mask_land_field%header%vertical_level = 1
1345 call storage_get_field(mask_land_field, interp_land_mask_status)
1347 end if
1348 if (interp_water_mask(idx) /= ' ') then
1349 mask_water_field%header%version = 1
1350 mask_water_field%header%forecast_hour = 0.
1351 mask_water_field%header%field = trim(interp_water_mask(idx))//'.mask'
1352 mask_water_field%header%vertical_coord = 'none'
1353 mask_water_field%header%vertical_level = 1
1355 call storage_get_field(mask_water_field, interp_water_mask_status)
1357 end if
1359 interp_array => interp_array_from_string(interp_method(idx))
1362 !
1363 ! Interpolate using average_gcell interpolation method
1364 !
1365 if (do_gcell_interp) then
1366 allocate(data_count(sm1:em1,sm2:em2))
1367 data_count = 0.
1369 if (interp_mask_status == 0) then
1370 call accum_continuous(slab, &
1371 minx+bdr, maxx-bdr, miny, maxy, 1, 1, &
1372 field%r_arr, data_count, &
1373 sm1, em1, sm2, em2, 1, 1, &
1374 istagger, &
1375 new_pts, NAN, interp_mask_val(idx), mask_field%r_arr)
1376 else
1377 call accum_continuous(slab, &
1378 minx+bdr, maxx-bdr, miny, maxy, 1, 1, &
1379 field%r_arr, data_count, &
1380 sm1, em1, sm2, em2, 1, 1, &
1381 istagger, &
1382 new_pts, NAN, -1.) ! The -1 is the maskval, but since we
1383 ! we do not give an optional mask, no
1384 ! no need to worry about -1s in data
1385 end if
1387 orig_selected_proj = iget_selected_domain()
1388 call select_domain(SOURCE_PROJ)
1389 do j=sm2,em2
1390 do i=sm1,em1
1392 if (present(landmask)) then
1394 if (landmask(i,j) /= masked(idx)) then
1395 if (data_count(i,j) > 0.) then
1396 field%r_arr(i,j) = field%r_arr(i,j) / data_count(i,j)
1397 else
1399 if (interp_mask_status == 0) then
1400 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1401 minx, maxx, miny, maxy, bdr, missing_value(idx), &
1402 mask_val=interp_mask_val(idx), mask_field=mask_field%r_arr)
1403 else
1404 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1405 minx, maxx, miny, maxy, bdr, missing_value(idx))
1406 end if
1408 if (temp /= missing_value(idx)) then
1409 field%r_arr(i,j) = temp
1410 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1411 end if
1413 end if
1414 else
1415 field%r_arr(i,j) = fill_missing(idx)
1416 end if
1418 if (.not. bitarray_test(new_pts, i-sm1+1, j-sm2+1) .and. &
1419 .not. bitarray_test(field%valid_mask, i-sm1+1, j-sm2+1)) then
1420 field%r_arr(i,j) = fill_missing(idx)
1421 end if
1423 else
1425 if (data_count(i,j) > 0.) then
1426 field%r_arr(i,j) = field%r_arr(i,j) / data_count(i,j)
1427 else
1429 if (interp_mask_status == 0) then
1430 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1431 minx, maxx, miny, maxy, bdr, missing_value(idx), &
1432 mask_val=interp_mask_val(idx), mask_field=mask_field%r_arr)
1433 else
1434 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1435 minx, maxx, miny, maxy, bdr, missing_value(idx))
1436 end if
1438 if (temp /= missing_value(idx)) then
1439 field%r_arr(i,j) = temp
1440 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1441 end if
1443 end if
1445 if (.not. bitarray_test(new_pts, i-sm1+1, j-sm2+1) .and. &
1446 .not. bitarray_test(field%valid_mask, i-sm1+1, j-sm2+1)) then
1447 field%r_arr(i,j) = fill_missing(idx)
1448 end if
1450 end if
1452 end do
1453 end do
1454 call select_domain(orig_selected_proj)
1455 deallocate(data_count)
1457 !
1458 ! No average_gcell interpolation method
1459 !
1460 else
1462 orig_selected_proj = iget_selected_domain()
1463 call select_domain(SOURCE_PROJ)
1464 do j=sm2,em2
1465 do i=sm1,em1
1466 if (present(landmask)) then
1468 if (masked(idx) == MASKED_BOTH) then
1470 if (landmask(i,j) == 0) then ! WATER POINT
1472 if (interp_land_mask_status == 0) then
1473 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1474 minx, maxx, miny, maxy, bdr, missing_value(idx), &
1475 mask_val=interp_land_mask_val(idx), mask_field=mask_land_field%r_arr)
1476 else
1477 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1478 minx, maxx, miny, maxy, bdr, missing_value(idx))
1479 end if
1481 else if (landmask(i,j) == 1) then ! LAND POINT
1483 if (interp_water_mask_status == 0) then
1484 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1485 minx, maxx, miny, maxy, bdr, missing_value(idx), &
1486 mask_val=interp_water_mask_val(idx), &
1487 mask_field=mask_water_field%r_arr)
1488 else
1489 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1490 minx, maxx, miny, maxy, bdr, missing_value(idx))
1491 end if
1493 end if
1495 else if (landmask(i,j) /= masked(idx)) then
1497 if (interp_mask_status == 0) then
1498 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1499 minx, maxx, miny, maxy, bdr, missing_value(idx), &
1500 mask_val=interp_mask_val(idx), mask_field=mask_field%r_arr)
1501 else
1502 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1503 minx, maxx, miny, maxy, bdr, missing_value(idx))
1504 end if
1506 else
1507 temp = missing_value(idx)
1508 end if
1510 else
1512 if (interp_mask_status == 0) then
1513 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1514 minx, maxx, miny, maxy, bdr, missing_value(idx), &
1515 mask_val=interp_mask_val(idx), mask_field=mask_field%r_arr)
1516 else
1517 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1518 minx, maxx, miny, maxy, bdr, missing_value(idx))
1519 end if
1521 end if
1523 if (temp /= missing_value(idx)) then
1524 field%r_arr(i,j) = temp
1525 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1526 end if
1528 if (.not. bitarray_test(new_pts, i-sm1+1, j-sm2+1) .and. &
1529 .not. bitarray_test(field%valid_mask, i-sm1+1, j-sm2+1)) then
1530 field%r_arr(i,j) = fill_missing(idx)
1531 end if
1533 end do
1534 end do
1535 call select_domain(orig_selected_proj)
1536 end if
1538 deallocate(interp_array)
1540 end subroutine interp_met_field
1543 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1544 ! Name: interp_to_latlon
1545 !
1546 ! Purpose:
1547 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1548 function interp_to_latlon(rlat, rlon, istagger, interp_method_list, slab, &
1549 minx, maxx, miny, maxy, bdr, source_missing_value, &
1550 mask_field, mask_val)
1552 use interp_module
1553 use llxy_module
1555 implicit none
1557 ! Arguments
1558 integer, intent(in) :: minx, maxx, miny, maxy, bdr, istagger
1559 integer, dimension(:), intent(in) :: interp_method_list
1560 real, intent(in) :: rlat, rlon, source_missing_value
1561 real, dimension(minx:maxx,miny:maxy), intent(in) :: slab
1562 real, intent(in), optional :: mask_val
1563 real, dimension(minx:maxx,miny:maxy), intent(in), optional :: mask_field
1565 ! Return value
1566 real :: interp_to_latlon
1568 ! Local variables
1569 real :: rx, ry
1571 interp_to_latlon = source_missing_value
1573 call lltoxy(rlat, rlon, rx, ry, istagger)
1574 if (rx >= minx+bdr-0.5 .and. rx <= maxx-bdr+0.5) then
1575 if (present(mask_field) .and. present(mask_val)) then
1576 interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, 1, 1, source_missing_value, &
1577 interp_method_list, 1, mask_val, mask_field)
1578 else
1579 interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, 1, 1, source_missing_value, &
1580 interp_method_list, 1)
1581 end if
1582 else
1583 interp_to_latlon = source_missing_value
1584 end if
1586 if (interp_to_latlon == source_missing_value) then
1588 ! Try a lon in the range 0. to 360.; all lons in the xlon
1589 ! array should be in the range -180. to 180.
1590 if (rlon < 0.) then
1591 call lltoxy(rlat, rlon+360., rx, ry, istagger)
1592 if (rx >= minx+bdr-0.5 .and. rx <= maxx-bdr+0.5) then
1593 if (present(mask_field) .and. present(mask_val)) then
1594 interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, &
1595 1, 1, source_missing_value, &
1596 interp_method_list, 1, mask_val, mask_field)
1597 else
1598 interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, &
1599 1, 1, source_missing_value, &
1600 interp_method_list, 1)
1601 end if
1602 else
1603 interp_to_latlon = source_missing_value
1604 end if
1606 end if
1608 end if
1610 return
1612 end function interp_to_latlon
1615 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1616 ! Name: get_bottom_top_dim
1617 !
1618 ! Purpose:
1619 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1620 subroutine get_bottom_top_dim(bottom_top_dim)
1622 use interp_option_module
1623 use list_module
1624 use storage_module
1626 implicit none
1628 ! Arguments
1629 integer, intent(out) :: bottom_top_dim
1631 ! Local variables
1632 integer :: i, j
1633 integer, pointer, dimension(:) :: field_levels
1634 character (len=32) :: z_dim
1635 type (fg_input), pointer, dimension(:) :: headers
1636 type (list) :: temp_levels
1638 ! Initialize a list to store levels that are found for 3-d fields
1639 call list_init(temp_levels)
1641 ! Get a list of all time-dependent fields (given by their headers) from
1642 ! the storage module
1643 call storage_get_td_headers(headers)
1645 !
1646 ! Given headers of all fields, we first build a list of all possible levels
1647 ! for 3-d met fields (excluding sea-level, though).
1648 !
1649 do i=1,size(headers)
1650 call get_z_dim_name(headers(i)%header%field, z_dim)
1652 ! We only want to consider 3-d met fields
1653 if (z_dim(1:18) == 'num_metgrid_levels') then
1655 ! Find out what levels the current field has
1656 call storage_get_levels(headers(i), field_levels)
1657 do j=1,size(field_levels)
1659 ! If this level has not yet been encountered, add it to our list
1660 if (.not. list_search(temp_levels, ikey=field_levels(j), ivalue=field_levels(j))) then
1661 if (field_levels(j) /= 201300) then
1662 call list_insert(temp_levels, ikey=field_levels(j), ivalue=field_levels(j))
1663 end if
1664 end if
1666 end do
1668 deallocate(field_levels)
1670 end if
1672 end do
1674 bottom_top_dim = list_length(temp_levels)
1676 call list_destroy(temp_levels)
1677 deallocate(headers)
1679 end subroutine get_bottom_top_dim
1682 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1683 ! Name: fill_missing_levels
1684 !
1685 ! Purpose:
1686 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1687 subroutine fill_missing_levels(output_flags)
1689 use interp_option_module
1690 use list_module
1691 use module_debug
1692 use module_mergesort
1693 use storage_module
1695 implicit none
1697 ! Arguments
1698 character (len=128), dimension(:), intent(inout) :: output_flags
1700 ! Local variables
1701 integer :: i, ii, j, ix, jx, k, lower, upper, temp, istatus
1702 integer, pointer, dimension(:) :: union_levels, field_levels
1703 real, pointer, dimension(:) :: r_union_levels
1704 character (len=128) :: clevel
1705 type (fg_input) :: lower_field, upper_field, new_field, search_field
1706 type (fg_input), pointer, dimension(:) :: headers, all_headers
1707 type (list) :: temp_levels
1708 type (list_item), pointer, dimension(:) :: keys
1710 ! Initialize a list to store levels that are found for 3-d fields
1711 call list_init(temp_levels)
1713 ! Get a list of all fields (given by their headers) from the storage module
1714 call storage_get_td_headers(headers)
1715 call storage_get_all_headers(all_headers)
1717 !
1718 ! Given headers of all fields, we first build a list of all possible levels
1719 ! for 3-d met fields (excluding sea-level, though).
1720 !
1721 do i=1,size(headers)
1723 ! Find out what levels the current field has
1724 call storage_get_levels(headers(i), field_levels)
1725 do j=1,size(field_levels)
1727 ! If this level has not yet been encountered, add it to our list
1728 if (.not. list_search(temp_levels, ikey=field_levels(j), ivalue=field_levels(j))) then
1729 if (field_levels(j) /= 201300) then
1730 call list_insert(temp_levels, ikey=field_levels(j), ivalue=field_levels(j))
1731 end if
1732 end if
1734 end do
1736 deallocate(field_levels)
1738 end do
1740 if (list_length(temp_levels) > 0) then
1742 !
1743 ! With all possible levels stored in a list, get an array of levels, sorted
1744 ! in decreasing order
1745 !
1746 i = 0
1747 allocate(union_levels(list_length(temp_levels)))
1748 do while (list_length(temp_levels) > 0)
1749 i = i + 1
1750 call list_get_first_item(temp_levels, ikey=union_levels(i), ivalue=temp)
1751 end do
1752 call mergesort(union_levels, 1, size(union_levels))
1754 allocate(r_union_levels(size(union_levels)))
1755 do i=1,size(union_levels)
1756 r_union_levels(i) = real(union_levels(i))
1757 end do
1759 !
1760 ! With a sorted, complete list of levels, we need
1761 ! to go back and fill in missing levels for each 3-d field
1762 !
1763 do i=1,size(headers)
1765 !
1766 ! Find entry in METGRID.TBL for this field, if one exists; if it does, then the
1767 ! entry may tell us how to get values for the current field at the missing level
1768 !
1769 do ii=1,num_entries
1770 if (fieldname(ii) == headers(i)%header%field) exit
1771 end do
1772 if (ii <= num_entries) then
1773 call dup(headers(i),new_field)
1774 nullify(new_field%valid_mask)
1775 nullify(new_field%modified_mask)
1776 call fill_field(new_field, ii, output_flags, r_union_levels)
1777 end if
1779 end do
1781 deallocate(union_levels)
1782 deallocate(r_union_levels)
1783 deallocate(headers)
1785 call storage_get_td_headers(headers)
1787 !
1788 ! Now we may need to vertically interpolate to missing values in 3-d fields
1789 !
1790 do i=1,size(headers)
1792 call storage_get_levels(headers(i), field_levels)
1794 ! If this isn't a 3-d array, nothing to do
1795 if (size(field_levels) > 1) then
1797 do k=1,size(field_levels)
1798 call dup(headers(i),search_field)
1799 search_field%header%vertical_level = field_levels(k)
1800 call storage_get_field(search_field,istatus)
1801 if (istatus == 0) then
1802 JLOOP: do jx=search_field%header%dim2(1),search_field%header%dim2(2)
1803 ILOOP: do ix=search_field%header%dim1(1),search_field%header%dim1(2)
1804 if (.not. bitarray_test(search_field%valid_mask, &
1805 ix-search_field%header%dim1(1)+1, &
1806 jx-search_field%header%dim2(1)+1)) then
1807 call bitarray_set(search_field%valid_mask, &
1808 ix-search_field%header%dim1(1)+1, &
1809 jx-search_field%header%dim2(1)+1)
1811 call dup(search_field, lower_field)
1812 do lower=k-1,1,-1
1813 lower_field%header%vertical_level = field_levels(lower)
1814 call storage_get_field(lower_field,istatus)
1815 if (bitarray_test(lower_field%valid_mask, &
1816 ix-search_field%header%dim1(1)+1, &
1817 jx-search_field%header%dim2(1)+1)) &
1818 exit
1820 end do
1822 call dup(search_field, upper_field)
1823 do upper=k+1,size(field_levels)
1824 upper_field%header%vertical_level = field_levels(upper)
1825 call storage_get_field(upper_field,istatus)
1826 if (bitarray_test(upper_field%valid_mask, &
1827 ix-search_field%header%dim1(1)+1, &
1828 jx-search_field%header%dim2(1)+1)) &
1829 exit
1831 end do
1832 if (upper <= size(field_levels) .and. lower >= 1) then
1833 search_field%r_arr(ix,jx) = real(abs(field_levels(upper)-field_levels(k))) &
1834 / real(abs(field_levels(upper)-field_levels(lower))) &
1835 * lower_field%r_arr(ix,jx) &
1836 + real(abs(field_levels(k)-field_levels(lower))) &
1837 / real(abs(field_levels(upper)-field_levels(lower))) &
1838 * upper_field%r_arr(ix,jx)
1839 end if
1840 end if
1841 end do ILOOP
1842 end do JLOOP
1843 else
1844 call mprintf(.true.,ERROR, &
1845 'This is bad, could not get %s at level %i.', &
1846 s1=trim(search_field%header%field), i1=field_levels(k))
1847 end if
1848 end do
1850 end if
1852 deallocate(field_levels)
1854 end do
1856 end if
1858 call list_destroy(temp_levels)
1859 deallocate(all_headers)
1860 deallocate(headers)
1862 end subroutine fill_missing_levels
1865 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1866 ! Name: create_derived_fields
1867 !
1868 ! Purpose:
1869 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1870 subroutine create_derived_fields(arg_gridtype, hdate, xfcst, &
1871 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
1872 we_mem_stag_s, we_mem_stag_e, sn_mem_stag_s, sn_mem_stag_e, &
1873 created_this_field, output_flags)
1875 use interp_option_module
1876 use list_module
1877 use module_mergesort
1878 use storage_module
1880 implicit none
1882 ! Arguments
1883 integer, intent(in) :: we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
1884 we_mem_stag_s, we_mem_stag_e, sn_mem_stag_s, sn_mem_stag_e
1885 real, intent(in) :: xfcst
1886 logical, dimension(:), intent(inout) :: created_this_field
1887 character (len=1), intent(in) :: arg_gridtype
1888 character (len=24), intent(in) :: hdate
1889 character (len=128), dimension(:), intent(inout) :: output_flags
1891 ! Local variables
1892 integer :: idx, i, j, istatus
1893 type (fg_input) :: field
1895 ! Initialize fg_input structure to store the field
1896 field%header%version = 1
1897 field%header%date = hdate//' '
1898 field%header%time_dependent = .true.
1899 field%header%mask_field = .false.
1900 field%header%forecast_hour = xfcst
1901 field%header%fg_source = 'Derived from FG'
1902 field%header%field = ' '
1903 field%header%units = ' '
1904 field%header%description = ' '
1905 field%header%vertical_level = 0
1906 field%header%array_order = 'XY '
1907 field%header%is_wind_grid_rel = .true.
1908 field%header%array_has_missing_values = .false.
1909 nullify(field%r_arr)
1910 nullify(field%valid_mask)
1911 nullify(field%modified_mask)
1913 !
1914 ! Check each entry in METGRID.TBL to see whether it is a derive field
1915 !
1916 do idx=1,num_entries
1917 if (is_derived_field(idx)) then
1919 created_this_field(idx) = .true.
1921 call mprintf(.true.,INFORM,'Going to create the field %s',s1=fieldname(idx))
1923 ! Intialize more fields in storage structure
1924 field%header%field = fieldname(idx)
1925 call get_z_dim_name(fieldname(idx),field%header%vertical_coord)
1926 field%map%stagger = output_stagger(idx)
1927 if (arg_gridtype == 'E') then
1928 field%header%dim1(1) = we_mem_s
1929 field%header%dim1(2) = we_mem_e
1930 field%header%dim2(1) = sn_mem_s
1931 field%header%dim2(2) = sn_mem_e
1932 else if (arg_gridtype == 'C') then
1933 if (output_stagger(idx) == M) then
1934 field%header%dim1(1) = we_mem_s
1935 field%header%dim1(2) = we_mem_e
1936 field%header%dim2(1) = sn_mem_s
1937 field%header%dim2(2) = sn_mem_e
1938 else if (output_stagger(idx) == U) then
1939 field%header%dim1(1) = we_mem_stag_s
1940 field%header%dim1(2) = we_mem_stag_e
1941 field%header%dim2(1) = sn_mem_s
1942 field%header%dim2(2) = sn_mem_e
1943 else if (output_stagger(idx) == V) then
1944 field%header%dim1(1) = we_mem_s
1945 field%header%dim1(2) = we_mem_e
1946 field%header%dim2(1) = sn_mem_stag_s
1947 field%header%dim2(2) = sn_mem_stag_e
1948 end if
1949 end if
1951 call fill_field(field, idx, output_flags)
1953 end if
1954 end do
1957 end subroutine create_derived_fields
1960 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1961 ! Name: fill_field
1962 !
1963 ! Purpose:
1964 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1965 subroutine fill_field(field, idx, output_flags, all_level_list)
1967 use interp_option_module
1968 use list_module
1969 use module_mergesort
1970 use storage_module
1972 implicit none
1974 ! Arguments
1975 integer, intent(in) :: idx
1976 type (fg_input), intent(inout) :: field
1977 character (len=128), dimension(:), intent(inout) :: output_flags
1978 real, dimension(:), intent(in), optional :: all_level_list
1980 ! Local variables
1981 integer :: i, j, istatus, isrclevel
1982 integer, pointer, dimension(:) :: all_list
1983 real :: rfillconst, rlevel, rsrclevel
1984 type (fg_input) :: query_field
1985 type (list_item), pointer, dimension(:) :: keys
1986 character (len=128) :: asrcname
1987 logical :: filled_all_lev
1989 filled_all_lev = .false.
1991 !
1992 ! Get a list of all levels to be filled for this field
1993 !
1994 keys => list_get_keys(fill_lev_list(idx))
1996 do i=1,list_length(fill_lev_list(idx))
1998 !
1999 ! First handle a specification for levels "all"
2000 !
2001 if (trim(keys(i)%ckey) == 'all') then
2003 ! We only want to fill all levels if we haven't already filled "all" of them
2004 if (.not. filled_all_lev) then
2006 filled_all_lev = .true.
2008 query_field%header%time_dependent = .true.
2009 query_field%header%field = ' '
2010 nullify(query_field%r_arr)
2011 nullify(query_field%valid_mask)
2012 nullify(query_field%modified_mask)
2014 ! See if we are filling this level with a constant
2015 call get_constant_fill_lev(keys(i)%cvalue, rfillconst, istatus)
2016 if (istatus == 0) then
2017 if (present(all_level_list)) then
2018 do j=1,size(all_level_list)
2019 call create_level(field, real(all_level_list(j)), idx, output_flags, rfillconst=rfillconst)
2020 end do
2021 else
2022 query_field%header%field = level_template(idx)
2023 nullify(all_list)
2024 call storage_get_levels(query_field, all_list)
2025 if (associated(all_list)) then
2026 do j=1,size(all_list)
2027 call create_level(field, real(all_list(j)), idx, output_flags, rfillconst=rfillconst)
2028 end do
2029 deallocate(all_list)
2030 end if
2031 end if
2033 ! Else see if we are filling this level with a constant equal
2034 ! to the value of the level
2035 else if (trim(keys(i)%cvalue) == 'vertical_index') then
2036 if (present(all_level_list)) then
2037 do j=1,size(all_level_list)
2038 call create_level(field, real(all_level_list(j)), idx, output_flags, &
2039 rfillconst=real(all_level_list(j)))
2040 end do
2041 else
2042 query_field%header%field = level_template(idx)
2043 nullify(all_list)
2044 call storage_get_levels(query_field, all_list)
2045 if (associated(all_list)) then
2046 do j=1,size(all_list)
2047 call create_level(field, real(all_list(j)), idx, output_flags, rfillconst=real(all_list(j)))
2048 end do
2049 deallocate(all_list)
2050 end if
2051 end if
2053 ! Else, we assume that it is a field from which we are copying levels
2054 else
2055 if (present(all_level_list)) then
2056 do j=1,size(all_level_list)
2057 call create_level(field, real(all_level_list(j)), idx, output_flags, &
2058 asrcname=keys(i)%cvalue, rsrclevel=real(all_level_list(j)))
2059 end do
2060 else
2061 query_field%header%field = keys(i)%cvalue ! Use same levels as source field, not level_template
2062 nullify(all_list)
2063 call storage_get_levels(query_field, all_list)
2064 if (associated(all_list)) then
2065 do j=1,size(all_list)
2066 call create_level(field, real(all_list(j)), idx, output_flags, &
2067 asrcname=keys(i)%cvalue, rsrclevel=real(all_list(j)))
2068 end do
2069 deallocate(all_list)
2071 else
2073 ! If the field doesn't have any levels (or does not exist) then we have not
2074 ! really filled all levels at this point.
2075 filled_all_lev = .false.
2076 end if
2077 end if
2079 end if
2080 end if
2082 !
2083 ! Handle individually specified levels
2084 !
2085 else
2087 read(keys(i)%ckey,*) rlevel
2089 ! See if we are filling this level with a constant
2090 call get_constant_fill_lev(keys(i)%cvalue, rfillconst, istatus)
2091 if (istatus == 0) then
2092 call create_level(field, rlevel, idx, output_flags, rfillconst=rfillconst)
2094 ! Otherwise, we are filling from another level
2095 else
2096 call get_fill_src_level(keys(i)%cvalue, asrcname, isrclevel)
2097 rsrclevel = real(isrclevel)
2098 call create_level(field, rlevel, idx, output_flags, &
2099 asrcname=asrcname, rsrclevel=rsrclevel)
2101 end if
2102 end if
2103 end do
2105 if (associated(keys)) deallocate(keys)
2107 end subroutine fill_field
2110 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2111 ! Name: create_level
2112 !
2113 ! Purpose:
2114 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2115 subroutine create_level(field_template, rlevel, idx, output_flags, &
2116 rfillconst, asrcname, rsrclevel)
2118 use storage_module
2119 use interp_option_module
2121 implicit none
2123 ! Arguments
2124 type (fg_input), intent(inout) :: field_template
2125 real, intent(in) :: rlevel
2126 integer, intent(in) :: idx
2127 character (len=128), dimension(:), intent(inout) :: output_flags
2128 real, intent(in), optional :: rfillconst, rsrclevel
2129 character (len=128), intent(in), optional :: asrcname
2131 ! Local variables
2132 integer :: i, j, istatus
2133 integer :: sm1,em1,sm2,em2
2134 type (fg_input) :: query_field
2136 !
2137 ! Check to make sure optional arguments are sane
2138 !
2139 if (present(rfillconst) .and. (present(asrcname) .or. present(rsrclevel))) then
2140 call mprintf(.true.,ERROR,'A call to create_level() cannot be given specifications '// &
2141 'for both a constant fill value and a source level.')
2143 else if ((present(asrcname) .and. .not. present(rsrclevel)) .or. &
2144 (.not. present(asrcname) .and. present(rsrclevel))) then
2145 call mprintf(.true.,ERROR,'Neither or both of optional arguments asrcname and '// &
2146 'rsrclevel must be specified to subroutine create_level().')
2148 else if (.not. present(rfillconst) .and. &
2149 .not. present(asrcname) .and. &
2150 .not. present(rsrclevel)) then
2151 call mprintf(.true.,ERROR,'A call to create_level() must be given either a specification '// &
2152 'for a constant fill value or a source level.')
2153 end if
2155 query_field%header%time_dependent = .true.
2156 query_field%header%field = field_template%header%field
2157 query_field%header%vertical_level = rlevel
2158 nullify(query_field%r_arr)
2159 nullify(query_field%valid_mask)
2160 nullify(query_field%modified_mask)
2162 call storage_query_field(query_field, istatus)
2163 if (istatus == 0) then
2164 call mprintf(.true.,INFORM,'%s at level %f already exists; leaving it alone.', &
2165 s1=field_template%header%field, f1=rlevel)
2166 return
2167 end if
2169 sm1 = field_template%header%dim1(1)
2170 em1 = field_template%header%dim1(2)
2171 sm2 = field_template%header%dim2(1)
2172 em2 = field_template%header%dim2(2)
2174 !
2175 ! Handle constant fill value case
2176 !
2177 if (present(rfillconst)) then
2179 field_template%header%vertical_level = rlevel
2180 allocate(field_template%r_arr(sm1:em1,sm2:em2))
2181 allocate(field_template%valid_mask)
2182 allocate(field_template%modified_mask)
2183 call bitarray_create(field_template%valid_mask, em1-sm1+1, em2-sm2+1)
2184 call bitarray_create(field_template%modified_mask, em1-sm1+1, em2-sm2+1)
2186 field_template%r_arr = rfillconst
2188 do j=sm2,em2
2189 do i=sm1,em1
2190 call bitarray_set(field_template%valid_mask, i-sm1+1, j-sm2+1)
2191 end do
2192 end do
2194 call storage_put_field(field_template)
2196 if (output_this_field(idx) .and. flag_in_output(idx) /= ' ') then
2197 output_flags(idx) = flag_in_output(idx)
2198 end if
2200 !
2201 ! Handle source field and source level case
2202 !
2203 else if (present(asrcname) .and. present(rsrclevel)) then
2205 query_field%header%field = ' '
2206 query_field%header%field = asrcname
2207 query_field%header%vertical_level = rsrclevel
2209 ! Check to see whether the requested source field exists at the requested level
2210 call storage_query_field(query_field, istatus)
2212 if (istatus == 0) then
2214 ! Read in requested field at requested level
2215 call storage_get_field(query_field, istatus)
2216 if ((query_field%header%dim1(1) /= field_template%header%dim1(1)) .or. &
2217 (query_field%header%dim1(2) /= field_template%header%dim1(2)) .or. &
2218 (query_field%header%dim2(1) /= field_template%header%dim2(1)) .or. &
2219 (query_field%header%dim2(2) /= field_template%header%dim2(2))) then
2220 call mprintf(.true.,ERROR,'Dimensions for %s do not match those of %s. This is '// &
2221 'probably because the staggerings of the fields do not match.', &
2222 s1=query_field%header%field, s2=field_template%header%field)
2223 end if
2225 field_template%header%vertical_level = rlevel
2226 allocate(field_template%r_arr(sm1:em1,sm2:em2))
2227 allocate(field_template%valid_mask)
2228 allocate(field_template%modified_mask)
2229 call bitarray_create(field_template%valid_mask, em1-sm1+1, em2-sm2+1)
2230 call bitarray_create(field_template%modified_mask, em1-sm1+1, em2-sm2+1)
2232 field_template%r_arr = query_field%r_arr
2234 ! We should retain information about which points in the field are valid
2235 do j=sm2,em2
2236 do i=sm1,em1
2237 if (bitarray_test(query_field%valid_mask, i-sm1+1, j-sm2+1)) then
2238 call bitarray_set(field_template%valid_mask, i-sm1+1, j-sm2+1)
2239 end if
2240 end do
2241 end do
2243 call storage_put_field(field_template)
2245 if (output_this_field(idx) .and. flag_in_output(idx) /= ' ') then
2246 output_flags(idx) = flag_in_output(idx)
2247 end if
2249 else
2250 call mprintf(.true.,INFORM,'Couldn''t find %s at level %f to fill level %f of %s.', &
2251 s1=asrcname,f1=rsrclevel,f2=rlevel,s2=field_template%header%field)
2252 end if
2254 end if
2256 end subroutine create_level
2259 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2260 ! Name: accum_continuous
2261 !
2262 ! Purpose: Sum up all of the source data points whose nearest neighbor in the
2263 ! model grid is the specified model grid point.
2264 !
2265 ! NOTE: When processing the source tile, those source points that are
2266 ! closest to a different model grid point will be added to the totals for
2267 ! such grid points; thus, an entire source tile will be processed at a time.
2268 ! This routine really processes for all model grid points that are
2269 ! within a source tile, and not just for a single grid point.
2270 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2271 subroutine accum_continuous(src_array, &
2272 src_min_x, src_max_x, src_min_y, src_max_y, src_min_z, src_max_z, &
2273 dst_array, n, &
2274 start_i, end_i, start_j, end_j, start_k, end_k, &
2275 istagger, &
2276 new_pts, msgval, maskval, mask_array)
2278 use bitarray_module
2279 use misc_definitions_module
2281 implicit none
2283 ! Arguments
2284 integer, intent(in) :: start_i, end_i, start_j, end_j, start_k, end_k, istagger, &
2285 src_min_x, src_max_x, src_min_y, src_max_y, src_min_z, src_max_z
2286 real, intent(in) :: maskval, msgval
2287 real, dimension(src_min_x:src_max_x, src_min_y:src_max_y, src_min_z:src_max_z), intent(in) :: src_array
2288 real, dimension(src_min_x:src_max_x, src_min_y:src_max_y), intent(in), optional :: mask_array
2289 real, dimension(start_i:end_i, start_j:end_j, start_k:end_k), intent(inout) :: dst_array, n
2290 type (bitarray), intent(inout) :: new_pts
2292 ! Local variables
2293 integer :: i, j
2294 integer, pointer, dimension(:,:,:) :: where_maps_to
2296 allocate(where_maps_to(src_min_x:src_max_x,src_min_y:src_max_y,2))
2297 do i=src_min_x,src_max_x
2298 do j=src_min_y,src_max_y
2299 where_maps_to(i,j,1) = NOT_PROCESSED
2300 end do
2301 end do
2303 call process_continuous_block(src_array, where_maps_to, &
2304 src_min_x, src_min_y, src_min_z, src_max_x, src_max_y, src_max_z, &
2305 src_min_x, src_min_y, src_min_z, &
2306 src_max_x, src_max_y, src_max_z, &
2307 dst_array, n, start_i, end_i, start_j, end_j, start_k, end_k, &
2308 istagger, &
2309 new_pts, msgval, maskval, mask_array)
2311 deallocate(where_maps_to)
2313 end subroutine accum_continuous
2316 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2317 ! Name: process_continuous_block
2318 !
2319 ! Purpose: To recursively process a subarray of continuous data, adding the
2320 ! points in a block to the sum for their nearest grid point. The nearest
2321 ! neighbor may be estimated in some cases; for example, if the four corners
2322 ! of a subarray all have the same nearest grid point, all elements in the
2323 ! subarray are added to that grid point.
2324 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2325 recursive subroutine process_continuous_block(tile_array, where_maps_to, &
2326 src_min_x, src_min_y, src_min_z, src_max_x, src_max_y, src_max_z, &
2327 min_i, min_j, min_k, max_i, max_j, max_k, &
2328 dst_array, n, &
2329 start_x, end_x, start_y, end_y, start_z, end_z, &
2330 istagger, &
2331 new_pts, msgval, maskval, mask_array)
2333 use bitarray_module
2334 use llxy_module
2335 use misc_definitions_module
2337 implicit none
2339 ! Arguments
2340 integer, intent(in) :: min_i, min_j, min_k, max_i, max_j, max_k, &
2341 src_min_x, src_min_y, src_min_z, src_max_x, src_max_y, src_max_z, &
2342 start_x, end_x, start_y, end_y, start_z, end_z, istagger
2343 integer, dimension(src_min_x:src_max_x,src_min_y:src_max_y,2), intent(inout) :: where_maps_to
2344 real, intent(in) :: maskval, msgval
2345 real, dimension(src_min_x:src_max_x,src_min_y:src_max_y,src_min_z:src_max_z), intent(in) :: tile_array
2346 real, dimension(src_min_x:src_max_x,src_min_y:src_max_y), intent(in), optional :: mask_array
2347 real, dimension(start_x:end_x,start_y:end_y,start_z:end_z), intent(inout) :: dst_array, n
2348 type (bitarray), intent(inout) :: new_pts
2350 ! Local variables
2351 integer :: orig_selected_domain, x_dest, y_dest, i, j, k, center_i, center_j
2352 real :: lat_corner, lon_corner, rx, ry
2354 ! Compute the model grid point that the corners of the rectangle to be
2355 ! processed map to
2356 ! Lower-left corner
2357 if (where_maps_to(min_i,min_j,1) == NOT_PROCESSED) then
2358 orig_selected_domain = iget_selected_domain()
2359 call select_domain(SOURCE_PROJ)
2360 call xytoll(real(min_i), real(min_j), lat_corner, lon_corner, istagger)
2361 call select_domain(1)
2362 call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
2363 call select_domain(orig_selected_domain)
2364 if (real(start_x) <= rx .and. rx <= real(end_x) .and. &
2365 real(start_y) <= ry .and. ry <= real(end_y)) then
2366 where_maps_to(min_i,min_j,1) = nint(rx)
2367 where_maps_to(min_i,min_j,2) = nint(ry)
2368 else
2369 where_maps_to(min_i,min_j,1) = OUTSIDE_DOMAIN
2370 end if
2371 end if
2373 ! Upper-left corner
2374 if (where_maps_to(min_i,max_j,1) == NOT_PROCESSED) then
2375 orig_selected_domain = iget_selected_domain()
2376 call select_domain(SOURCE_PROJ)
2377 call xytoll(real(min_i), real(max_j), lat_corner, lon_corner, istagger)
2378 call select_domain(1)
2379 call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
2380 call select_domain(orig_selected_domain)
2381 if (real(start_x) <= rx .and. rx <= real(end_x) .and. &
2382 real(start_y) <= ry .and. ry <= real(end_y)) then
2383 where_maps_to(min_i,max_j,1) = nint(rx)
2384 where_maps_to(min_i,max_j,2) = nint(ry)
2385 else
2386 where_maps_to(min_i,max_j,1) = OUTSIDE_DOMAIN
2387 end if
2388 end if
2390 ! Upper-right corner
2391 if (where_maps_to(max_i,max_j,1) == NOT_PROCESSED) then
2392 orig_selected_domain = iget_selected_domain()
2393 call select_domain(SOURCE_PROJ)
2394 call xytoll(real(max_i), real(max_j), lat_corner, lon_corner, istagger)
2395 call select_domain(1)
2396 call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
2397 call select_domain(orig_selected_domain)
2398 if (real(start_x) <= rx .and. rx <= real(end_x) .and. &
2399 real(start_y) <= ry .and. ry <= real(end_y)) then
2400 where_maps_to(max_i,max_j,1) = nint(rx)
2401 where_maps_to(max_i,max_j,2) = nint(ry)
2402 else
2403 where_maps_to(max_i,max_j,1) = OUTSIDE_DOMAIN
2404 end if
2405 end if
2407 ! Lower-right corner
2408 if (where_maps_to(max_i,min_j,1) == NOT_PROCESSED) then
2409 orig_selected_domain = iget_selected_domain()
2410 call select_domain(SOURCE_PROJ)
2411 call xytoll(real(max_i), real(min_j), lat_corner, lon_corner, istagger)
2412 call select_domain(1)
2413 call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
2414 call select_domain(orig_selected_domain)
2415 if (real(start_x) <= rx .and. rx <= real(end_x) .and. &
2416 real(start_y) <= ry .and. ry <= real(end_y)) then
2417 where_maps_to(max_i,min_j,1) = nint(rx)
2418 where_maps_to(max_i,min_j,2) = nint(ry)
2419 else
2420 where_maps_to(max_i,min_j,1) = OUTSIDE_DOMAIN
2421 end if
2422 end if
2424 ! If all four corners map to same model grid point, accumulate the
2425 ! entire rectangle
2426 if (where_maps_to(min_i,min_j,1) == where_maps_to(min_i,max_j,1) .and. &
2427 where_maps_to(min_i,min_j,1) == where_maps_to(max_i,max_j,1) .and. &
2428 where_maps_to(min_i,min_j,1) == where_maps_to(max_i,min_j,1) .and. &
2429 where_maps_to(min_i,min_j,2) == where_maps_to(min_i,max_j,2) .and. &
2430 where_maps_to(min_i,min_j,2) == where_maps_to(max_i,max_j,2) .and. &
2431 where_maps_to(min_i,min_j,2) == where_maps_to(max_i,min_j,2) .and. &
2432 where_maps_to(min_i,min_j,1) /= OUTSIDE_DOMAIN) then
2433 x_dest = where_maps_to(min_i,min_j,1)
2434 y_dest = where_maps_to(min_i,min_j,2)
2436 ! If this grid point was already given a value from higher-priority source data,
2437 ! there is nothing to do.
2438 ! if (.not. bitarray_test(processed_pts, x_dest-start_x+1, y_dest-start_y+1)) then
2440 ! If this grid point has never been given a value by this level of source data,
2441 ! initialize the point
2442 if (.not. bitarray_test(new_pts, x_dest-start_x+1, y_dest-start_y+1)) then
2443 do k=min_k,max_k
2444 dst_array(x_dest,y_dest,k) = 0.
2445 end do
2446 end if
2448 ! Sum all the points whose nearest neighbor is this grid point
2449 if (present(mask_array)) then
2450 do i=min_i,max_i
2451 do j=min_j,max_j
2452 do k=min_k,max_k
2453 ! Ignore masked/missing values in the source data
2454 if ((tile_array(i,j,k) /= msgval) .and. &
2455 (mask_array(i,j) /= maskval)) then
2456 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2457 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2458 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2459 end if
2460 end do
2461 end do
2462 end do
2463 else
2464 do i=min_i,max_i
2465 do j=min_j,max_j
2466 do k=min_k,max_k
2467 ! Ignore masked/missing values in the source data
2468 if ((tile_array(i,j,k) /= msgval)) then
2469 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2470 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2471 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2472 end if
2473 end do
2474 end do
2475 end do
2476 end if
2478 ! end if
2480 ! Rectangle is a square of four points, and we can simply deal with each of the points
2481 else if (((max_i - min_i + 1) <= 2) .and. ((max_j - min_j + 1) <= 2)) then
2482 do i=min_i,max_i
2483 do j=min_j,max_j
2484 x_dest = where_maps_to(i,j,1)
2485 y_dest = where_maps_to(i,j,2)
2487 if (x_dest /= OUTSIDE_DOMAIN) then
2489 ! if (.not. bitarray_test(processed_pts, x_dest-start_x+1, y_dest-start_y+1)) then
2490 if (.not. bitarray_test(new_pts, x_dest-start_x+1, y_dest-start_y+1)) then
2491 do k=min_k,max_k
2492 dst_array(x_dest,y_dest,k) = 0.
2493 end do
2494 end if
2496 if (present(mask_array)) then
2497 do k=min_k,max_k
2498 ! Ignore masked/missing values
2499 if ((tile_array(i,j,k) /= msgval) .and. &
2500 (mask_array(i,j) /= maskval)) then
2501 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2502 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2503 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2504 end if
2505 end do
2506 else
2507 do k=min_k,max_k
2508 ! Ignore masked/missing values
2509 if ((tile_array(i,j,k) /= msgval)) then
2510 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2511 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2512 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2513 end if
2514 end do
2515 end if
2516 ! end if
2518 end if
2519 end do
2520 end do
2522 ! Not all corners map to the same grid point, and the rectangle contains more than
2523 ! four points
2524 else
2525 center_i = (max_i + min_i)/2
2526 center_j = (max_j + min_j)/2
2528 ! Recursively process lower-left rectangle
2529 call process_continuous_block(tile_array, where_maps_to, &
2530 src_min_x, src_min_y, src_min_z, &
2531 src_max_x, src_max_y, src_max_z, &
2532 min_i, min_j, min_k, &
2533 center_i, center_j, max_k, &
2534 dst_array, n, &
2535 start_x, end_x, start_y, end_y, start_z, end_z, &
2536 istagger, &
2537 new_pts, msgval, maskval, mask_array)
2539 if (center_i < max_i) then
2540 ! Recursively process lower-right rectangle
2541 call process_continuous_block(tile_array, where_maps_to, &
2542 src_min_x, src_min_y, src_min_z, &
2543 src_max_x, src_max_y, src_max_z, &
2544 center_i+1, min_j, min_k, max_i, &
2545 center_j, max_k, &
2546 dst_array, n, &
2547 start_x, end_x, start_y, &
2548 end_y, start_z, end_z, &
2549 istagger, &
2550 new_pts, msgval, maskval, mask_array)
2551 end if
2553 if (center_j < max_j) then
2554 ! Recursively process upper-left rectangle
2555 call process_continuous_block(tile_array, where_maps_to, &
2556 src_min_x, src_min_y, src_min_z, &
2557 src_max_x, src_max_y, src_max_z, &
2558 min_i, center_j+1, min_k, center_i, &
2559 max_j, max_k, &
2560 dst_array, n, &
2561 start_x, end_x, start_y, &
2562 end_y, start_z, end_z, &
2563 istagger, &
2564 new_pts, msgval, maskval, mask_array)
2565 end if
2567 if (center_i < max_i .and. center_j < max_j) then
2568 ! Recursively process upper-right rectangle
2569 call process_continuous_block(tile_array, where_maps_to, &
2570 src_min_x, src_min_y, src_min_z, &
2571 src_max_x, src_max_y, src_max_z, &
2572 center_i+1, center_j+1, min_k, max_i, &
2573 max_j, max_k, &
2574 dst_array, n, &
2575 start_x, end_x, start_y, &
2576 end_y, start_z, end_z, &
2577 istagger, &
2578 new_pts, msgval, maskval, mask_array)
2579 end if
2580 end if
2582 end subroutine process_continuous_block
2584 end module process_domain_module