| blob | 79c1461e5467a827f2c47fcf49b7f6094dee74bb |
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
37 real :: cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
38 dom_dx, dom_dy, pole_lat, pole_lon
39 real, dimension(16) :: corner_lats, corner_lons
40 real, pointer, dimension(:,:) :: landmask
41 real, pointer, dimension(:,:) :: xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v
42 logical, allocatable, dimension(:) :: got_this_field, got_const_field
43 character (len=19) :: valid_date, temp_date
44 character (len=128) :: title, mminlu
45 character (len=128), allocatable, dimension(:) :: output_flags, td_output_flags
47 ! Compute number of times that we will process
48 call geth_idts(end_date(n), start_date(n), idiff)
49 call mprintf((idiff < 0),ERROR,'Ending date is earlier than starting date in namelist for domain %i.', i1=n)
51 n_times = idiff / interval_seconds
53 ! Check that the interval evenly divides the range of times to process
54 call mprintf((mod(idiff, interval_seconds) /= 0),WARN, &
55 'In namelist, interval_seconds does not evenly divide '// &
56 '(end_date - start_date) for domain %i. Only %i time periods '// &
57 'will be processed.', i1=n, i2=n_times)
59 ! Initialize the storage module
60 call mprintf(.true.,LOGFILE,'Initializing storage module')
61 call storage_init()
63 !
64 ! Do time-independent processing
65 !
66 call get_static_fields(n, dyn_opt, west_east_dim, south_north_dim, bottom_top_dim, map_proj, &
67 we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
68 we_patch_s, we_patch_e, &
69 we_patch_stag_s, we_patch_stag_e, &
70 sn_patch_s, sn_patch_e, &
71 sn_patch_stag_s, sn_patch_stag_e, &
72 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
73 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
74 mminlu, num_land_cat, is_water, is_lake, is_ice, is_urban, i_soilwater, &
75 grid_id, parent_id, &
76 i_parent_start, j_parent_start, i_parent_end, j_parent_end, &
77 parent_grid_ratio, sub_x, sub_y, &
78 cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
79 pole_lat, pole_lon, dom_dx, dom_dy, landmask, xlat, xlon, xlat_u, xlon_u, &
80 xlat_v, xlon_v, corner_lats, corner_lons, title)
83 allocate(output_flags(num_entries))
84 allocate(got_const_field(num_entries))
86 do i=1,num_entries
87 output_flags(i) = ' '
88 got_const_field(i) = .false.
89 end do
91 ! This call is to process the constant met fields (SST or SEAICE, for example)
92 ! That we process constant fields is indicated by the first argument
93 call process_single_met_time(.true., temp_date, n, extra_row, extra_col, xlat, xlon, &
94 xlat_u, xlon_u, xlat_v, xlon_v, landmask, &
95 title, dyn_opt, &
96 west_east_dim, south_north_dim, &
97 we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
98 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
99 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
100 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
101 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
102 got_const_field, &
103 map_proj, mminlu, num_land_cat, is_water, is_lake, is_ice, is_urban, i_soilwater, &
104 grid_id, parent_id, i_parent_start, &
105 j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
106 cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
107 pole_lat, pole_lon, parent_grid_ratio, sub_x, sub_y, &
108 corner_lats, corner_lons, output_flags)
110 !
111 ! Begin time-dependent processing
112 !
114 allocate(td_output_flags(num_entries))
115 allocate(got_this_field (num_entries))
117 ! Loop over all times to be processed for this domain
118 do t=0,n_times
120 call geth_newdate(valid_date, trim(start_date(n)), t*interval_seconds)
121 temp_date = ' '
123 if (mod(interval_seconds,3600) == 0) then
124 write(temp_date,'(a13)') valid_date(1:10)//'_'//valid_date(12:13)
125 else if (mod(interval_seconds,60) == 0) then
126 write(temp_date,'(a16)') valid_date(1:10)//'_'//valid_date(12:16)
127 else
128 write(temp_date,'(a19)') valid_date(1:10)//'_'//valid_date(12:19)
129 end if
131 call mprintf(.true.,STDOUT, ' Processing %s', s1=trim(temp_date))
132 call mprintf(.true.,LOGFILE, 'Preparing to process output time %s', s1=temp_date)
134 do i=1,num_entries
135 td_output_flags(i) = output_flags(i)
136 got_this_field(i) = got_const_field(i)
137 end do
139 call process_single_met_time(.false., temp_date, n, extra_row, extra_col, xlat, xlon, &
140 xlat_u, xlon_u, xlat_v, xlon_v, landmask, &
141 title, dyn_opt, &
142 west_east_dim, south_north_dim, &
143 we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
144 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
145 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
146 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
147 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
148 got_this_field, &
149 map_proj, mminlu, num_land_cat, is_water, is_lake, is_ice, is_urban, i_soilwater, &
150 grid_id, parent_id, i_parent_start, &
151 j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
152 cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
153 pole_lat, pole_lon, parent_grid_ratio, sub_x, sub_y, &
154 corner_lats, corner_lons, td_output_flags)
156 end do ! Loop over n_times
159 deallocate(td_output_flags)
160 deallocate(got_this_field)
162 deallocate(output_flags)
163 deallocate(got_const_field)
165 call storage_delete_all()
167 end subroutine process_domain
170 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
171 ! Name: get_static_fields
172 !
173 ! Purpose:
174 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
175 subroutine get_static_fields(n, dyn_opt, west_east_dim, south_north_dim, bottom_top_dim, &
176 map_proj, &
177 we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
178 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
179 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
180 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
181 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
182 mminlu, num_land_cat, is_water, is_lake, is_ice, is_urban, i_soilwater, &
183 grid_id, parent_id, &
184 i_parent_start, j_parent_start, i_parent_end, j_parent_end, &
185 parent_grid_ratio, sub_x, sub_y, &
186 cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
187 pole_lat, pole_lon, dom_dx, dom_dy, landmask, xlat, xlon, xlat_u, xlon_u, &
188 xlat_v, xlon_v, corner_lats, corner_lons, title)
190 use gridinfo_module
191 use input_module
192 use llxy_module
193 use parallel_module
194 use storage_module
195 use module_debug
197 implicit none
199 ! Arguments
200 integer, intent(in) :: n
201 integer, intent(inout) :: dyn_opt, west_east_dim, south_north_dim, bottom_top_dim, &
202 map_proj, &
203 we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
204 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
205 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
206 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
207 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
208 is_water, is_lake, is_ice, is_urban, i_soilwater, grid_id, parent_id, &
209 i_parent_start, j_parent_start, i_parent_end, j_parent_end, &
210 parent_grid_ratio, sub_x, sub_y, num_land_cat
211 real, pointer, dimension(:,:) :: landmask
212 real, intent(inout) :: cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
213 dom_dx, dom_dy, pole_lat, pole_lon
214 real, pointer, dimension(:,:) :: xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v
215 real, dimension(16), intent(out) :: corner_lats, corner_lons
216 character (len=128), intent(inout) :: title, mminlu
218 ! Local variables
219 integer :: istatus, i, j, k, sp1, ep1, sp2, ep2, sp3, ep3, &
220 lh_mult, rh_mult, bh_mult, th_mult, subx, suby
221 integer :: we_mem_subgrid_s, we_mem_subgrid_e, &
222 sn_mem_subgrid_s, sn_mem_subgrid_e
223 integer :: we_patch_subgrid_s, we_patch_subgrid_e, &
224 sn_patch_subgrid_s, sn_patch_subgrid_e
225 real, pointer, dimension(:,:,:) :: real_array
226 character (len=3) :: memorder
227 character (len=128) :: grid_type, datestr, cname, stagger, cunits, cdesc
228 character (len=128), dimension(3) :: dimnames
229 type (fg_input) :: field
231 ! Initialize the input module to read static input data for this domain
232 call mprintf(.true.,LOGFILE,'Opening static input file.')
233 call input_init(n, istatus)
234 call mprintf((istatus /= 0),ERROR, 'input_init(): Error opening input for domain %i.', i1=n)
236 ! Read global attributes from the static data input file
237 call mprintf(.true.,LOGFILE,'Reading static global attributes.')
238 call read_global_attrs(title, datestr, grid_type, dyn_opt, west_east_dim, &
239 south_north_dim, bottom_top_dim, &
240 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
241 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
242 map_proj, mminlu, num_land_cat, &
243 is_water, is_lake, is_ice, is_urban, i_soilwater, &
244 grid_id, parent_id, i_parent_start, &
245 j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
246 cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, &
247 truelat2, pole_lat, pole_lon, parent_grid_ratio, &
248 corner_lats, corner_lons, sub_x, sub_y)
250 we_dom_s = 1
251 sn_dom_s = 1
252 if (grid_type(1:1) == 'C') then
253 we_dom_e = west_east_dim - 1
254 sn_dom_e = south_north_dim - 1
255 else if (grid_type(1:1) == 'E') then
256 we_dom_e = west_east_dim
257 sn_dom_e = south_north_dim
258 end if
260 ! Given the full dimensions of this domain, find out the range of indices
261 ! that will be worked on by this processor. This information is given by
262 ! my_minx, my_miny, my_maxx, my_maxy
263 call parallel_get_tile_dims(west_east_dim, south_north_dim)
265 ! Must figure out patch dimensions from info in parallel module
266 if (nprocs > 1 .and. .not. do_tiled_input) then
268 we_patch_s = my_minx
269 we_patch_stag_s = my_minx
270 we_patch_e = my_maxx - 1
271 sn_patch_s = my_miny
272 sn_patch_stag_s = my_miny
273 sn_patch_e = my_maxy - 1
275 if (gridtype == 'C') then
276 if (my_x /= nproc_x - 1) then
277 we_patch_e = we_patch_e + 1
278 we_patch_stag_e = we_patch_e
279 else
280 we_patch_stag_e = we_patch_e + 1
281 end if
282 if (my_y /= nproc_y - 1) then
283 sn_patch_e = sn_patch_e + 1
284 sn_patch_stag_e = sn_patch_e
285 else
286 sn_patch_stag_e = sn_patch_e + 1
287 end if
288 else if (gridtype == 'E') then
289 we_patch_e = we_patch_e + 1
290 sn_patch_e = sn_patch_e + 1
291 we_patch_stag_e = we_patch_e
292 sn_patch_stag_e = sn_patch_e
293 end if
295 end if
297 ! Compute multipliers for halo width; these must be 0/1
298 if (my_x /= 0) then
299 lh_mult = 1
300 else
301 lh_mult = 0
302 end if
303 if (my_x /= (nproc_x-1)) then
304 rh_mult = 1
305 else
306 rh_mult = 0
307 end if
308 if (my_y /= 0) then
309 bh_mult = 1
310 else
311 bh_mult = 0
312 end if
313 if (my_y /= (nproc_y-1)) then
314 th_mult = 1
315 else
316 th_mult = 0
317 end if
319 we_mem_s = we_patch_s - HALO_WIDTH*lh_mult
320 we_mem_e = we_patch_e + HALO_WIDTH*rh_mult
321 sn_mem_s = sn_patch_s - HALO_WIDTH*bh_mult
322 sn_mem_e = sn_patch_e + HALO_WIDTH*th_mult
323 we_mem_stag_s = we_patch_stag_s - HALO_WIDTH*lh_mult
324 we_mem_stag_e = we_patch_stag_e + HALO_WIDTH*rh_mult
325 sn_mem_stag_s = sn_patch_stag_s - HALO_WIDTH*bh_mult
326 sn_mem_stag_e = sn_patch_stag_e + HALO_WIDTH*th_mult
328 ! Initialize a proj_info type for the destination grid projection. This will
329 ! primarily be used for rotating Earth-relative winds to grid-relative winds
330 call set_domain_projection(map_proj, stand_lon, truelat1, truelat2, &
331 dom_dx, dom_dy, dom_dx, dom_dy, west_east_dim, &
332 south_north_dim, real(west_east_dim)/2., &
333 real(south_north_dim)/2.,cen_lat, cen_lon, &
334 cen_lat, cen_lon)
336 ! Read static fields using the input module; we know that there are no more
337 ! fields to be read when read_next_field() returns a non-zero status.
338 istatus = 0
339 do while (istatus == 0)
340 call read_next_field(sp1, ep1, sp2, ep2, sp3, ep3, cname, cunits, cdesc, &
341 memorder, stagger, dimnames, subx, suby, &
342 real_array, istatus)
343 if (istatus == 0) then
345 call mprintf(.true.,LOGFILE, 'Read in static field %s.',s1=cname)
347 ! We will also keep copies in core of the lat/lon arrays, for use in
348 ! interpolation of the met fields to the model grid.
349 ! For now, we assume that the lat/lon arrays will have known field names
350 if (index(cname, 'XLAT_M') /= 0 .and. &
351 len_trim(cname) == len_trim('XLAT_M')) then
352 allocate(xlat(we_mem_s:we_mem_e,sn_mem_s:sn_mem_e))
353 xlat(we_patch_s:we_patch_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
354 call exchange_halo_r(xlat, &
355 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, 1, 1, &
356 we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, 1, 1)
358 else if (index(cname, 'XLONG_M') /= 0 .and. &
359 len_trim(cname) == len_trim('XLONG_M')) then
360 allocate(xlon(we_mem_s:we_mem_e,sn_mem_s:sn_mem_e))
361 xlon(we_patch_s:we_patch_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
362 call exchange_halo_r(xlon, &
363 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, 1, 1, &
364 we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, 1, 1)
366 else if (index(cname, 'XLAT_U') /= 0 .and. &
367 len_trim(cname) == len_trim('XLAT_U')) then
368 allocate(xlat_u(we_mem_stag_s:we_mem_stag_e,sn_mem_s:sn_mem_e))
369 xlat_u(we_patch_stag_s:we_patch_stag_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
370 call exchange_halo_r(xlat_u, &
371 we_mem_stag_s, we_mem_stag_e, sn_mem_s, sn_mem_e, 1, 1, &
372 we_patch_stag_s, we_patch_stag_e, sn_patch_s, sn_patch_e, 1, 1)
374 else if (index(cname, 'XLONG_U') /= 0 .and. &
375 len_trim(cname) == len_trim('XLONG_U')) then
376 allocate(xlon_u(we_mem_stag_s:we_mem_stag_e,sn_mem_s:sn_mem_e))
377 xlon_u(we_patch_stag_s:we_patch_stag_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
378 call exchange_halo_r(xlon_u, &
379 we_mem_stag_s, we_mem_stag_e, sn_mem_s, sn_mem_e, 1, 1, &
380 we_patch_stag_s, we_patch_stag_e, sn_patch_s, sn_patch_e, 1, 1)
382 else if (index(cname, 'XLAT_V') /= 0 .and. &
383 len_trim(cname) == len_trim('XLAT_V')) then
384 allocate(xlat_v(we_mem_s:we_mem_e,sn_mem_stag_s:sn_mem_stag_e))
385 xlat_v(we_patch_s:we_patch_e,sn_patch_stag_s:sn_patch_stag_e) = real_array(:,:,1)
386 call exchange_halo_r(xlat_v, &
387 we_mem_s, we_mem_e, sn_mem_stag_s, sn_mem_stag_e, 1, 1, &
388 we_patch_s, we_patch_e, sn_patch_stag_s, sn_patch_stag_e, 1, 1)
390 else if (index(cname, 'XLONG_V') /= 0 .and. &
391 len_trim(cname) == len_trim('XLONG_V')) then
392 allocate(xlon_v(we_mem_s:we_mem_e,sn_mem_stag_s:sn_mem_stag_e))
393 xlon_v(we_patch_s:we_patch_e,sn_patch_stag_s:sn_patch_stag_e) = real_array(:,:,1)
394 call exchange_halo_r(xlon_v, &
395 we_mem_s, we_mem_e, sn_mem_stag_s, sn_mem_stag_e, 1, 1, &
396 we_patch_s, we_patch_e, sn_patch_stag_s, sn_patch_stag_e, 1, 1)
398 else if (index(cname, 'LANDMASK') /= 0 .and. &
399 len_trim(cname) == len_trim('LANDMASK')) then
400 allocate(landmask(we_mem_s:we_mem_e,sn_mem_s:sn_mem_e))
401 landmask(we_patch_s:we_patch_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
402 call exchange_halo_r(landmask, &
403 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, 1, 1, &
404 we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, 1, 1)
406 end if
408 if (subx > 1) then
409 we_mem_subgrid_s = (we_mem_s + HALO_WIDTH*lh_mult - 1)*subx - HALO_WIDTH*lh_mult + 1
410 we_mem_subgrid_e = (we_mem_e + (1-rh_mult) - HALO_WIDTH*rh_mult )*subx + HALO_WIDTH*rh_mult
411 we_patch_subgrid_s = (we_patch_s - 1)*subx + 1
412 we_patch_subgrid_e = (we_patch_e + (1-rh_mult) )*subx
413 end if
414 if (suby > 1) then
415 sn_mem_subgrid_s = (sn_mem_s + HALO_WIDTH*bh_mult - 1)*suby - HALO_WIDTH*bh_mult + 1
416 sn_mem_subgrid_e = (sn_mem_e + (1-th_mult) - HALO_WIDTH*th_mult )*suby + HALO_WIDTH*th_mult
417 sn_patch_subgrid_s = (sn_patch_s - 1)*suby + 1
418 sn_patch_subgrid_e = (sn_patch_e + (1-th_mult) )*suby
419 end if
421 ! Having read in a field, we write each level individually to the
422 ! storage module; levels will be reassembled later on when they
423 ! are written.
424 do k=sp3,ep3
425 field%header%sr_x=subx
426 field%header%sr_y=suby
427 field%header%version = 1
428 field%header%date = start_date(n)
429 field%header%time_dependent = .false.
430 field%header%mask_field = .false.
431 field%header%forecast_hour = 0.0
432 field%header%fg_source = 'geogrid_model'
433 field%header%field = cname
434 field%header%units = cunits
435 field%header%description = cdesc
436 field%header%vertical_coord = dimnames(3)
437 field%header%vertical_level = k
438 field%header%array_order = memorder
439 field%header%is_wind_grid_rel = .true.
440 field%header%array_has_missing_values = .false.
441 if (gridtype == 'C') then
442 if (subx > 1 .or. suby > 1) then
443 field%map%stagger = M
444 field%header%dim1(1) = we_mem_subgrid_s
445 field%header%dim1(2) = we_mem_subgrid_e
446 field%header%dim2(1) = sn_mem_subgrid_s
447 field%header%dim2(2) = sn_mem_subgrid_e
448 else if (trim(stagger) == 'M') then
449 field%map%stagger = M
450 field%header%dim1(1) = we_mem_s
451 field%header%dim1(2) = we_mem_e
452 field%header%dim2(1) = sn_mem_s
453 field%header%dim2(2) = sn_mem_e
454 else if (trim(stagger) == 'U') then
455 field%map%stagger = U
456 field%header%dim1(1) = we_mem_stag_s
457 field%header%dim1(2) = we_mem_stag_e
458 field%header%dim2(1) = sn_mem_s
459 field%header%dim2(2) = sn_mem_e
460 else if (trim(stagger) == 'V') then
461 field%map%stagger = V
462 field%header%dim1(1) = we_mem_s
463 field%header%dim1(2) = we_mem_e
464 field%header%dim2(1) = sn_mem_stag_s
465 field%header%dim2(2) = sn_mem_stag_e
466 end if
467 else if (gridtype == 'E') then
468 if (trim(stagger) == 'M') then
469 field%map%stagger = HH
470 else if (trim(stagger) == 'V') then
471 field%map%stagger = VV
472 end if
473 field%header%dim1(1) = we_mem_s
474 field%header%dim1(2) = we_mem_e
475 field%header%dim2(1) = sn_mem_s
476 field%header%dim2(2) = sn_mem_e
477 end if
479 allocate(field%valid_mask)
481 if (subx > 1 .or. suby > 1) then
482 allocate(field%r_arr(we_mem_subgrid_s:we_mem_subgrid_e,&
483 sn_mem_subgrid_s:sn_mem_subgrid_e))
484 field%r_arr(we_patch_subgrid_s:we_patch_subgrid_e,sn_patch_subgrid_s:sn_patch_subgrid_e) = &
485 real_array(sp1:ep1,sp2:ep2,k)
486 call exchange_halo_r(field%r_arr, &
487 we_mem_subgrid_s, we_mem_subgrid_e, sn_mem_subgrid_s, sn_mem_subgrid_e, 1, 1, &
488 we_patch_subgrid_s, we_patch_subgrid_e, sn_patch_subgrid_s, sn_patch_subgrid_e, 1, 1)
489 call bitarray_create(field%valid_mask, &
490 (we_mem_subgrid_e-we_mem_subgrid_s)+1, &
491 (sn_mem_subgrid_e-sn_mem_subgrid_s)+1)
492 do j=1,(sn_mem_subgrid_e-sn_mem_subgrid_s)+1
493 do i=1,(we_mem_subgrid_e-we_mem_subgrid_s)+1
494 call bitarray_set(field%valid_mask, i, j)
495 end do
496 end do
498 else if (field%map%stagger == M .or. &
499 field%map%stagger == HH .or. &
500 field%map%stagger == VV) then
501 allocate(field%r_arr(we_mem_s:we_mem_e,&
502 sn_mem_s:sn_mem_e))
503 field%r_arr(we_patch_s:we_patch_e,sn_patch_s:sn_patch_e) = real_array(sp1:ep1,sp2:ep2,k)
504 call exchange_halo_r(field%r_arr, &
505 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, 1, 1, &
506 we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, 1, 1)
507 call bitarray_create(field%valid_mask, &
508 (we_mem_e-we_mem_s)+1, &
509 (sn_mem_e-sn_mem_s)+1)
510 do j=1,(sn_mem_e-sn_mem_s)+1
511 do i=1,(we_mem_e-we_mem_s)+1
512 call bitarray_set(field%valid_mask, i, j)
513 end do
514 end do
515 else if (field%map%stagger == U) then
516 allocate(field%r_arr(we_mem_stag_s:we_mem_stag_e,&
517 sn_mem_s:sn_mem_e))
518 field%r_arr(we_patch_stag_s:we_patch_stag_e,sn_patch_s:sn_patch_e) = real_array(sp1:ep1,sp2:ep2,k)
519 call exchange_halo_r(field%r_arr, &
520 we_mem_stag_s, we_mem_stag_e, sn_mem_s, sn_mem_e, 1, 1, &
521 we_patch_stag_s, we_patch_stag_e, sn_patch_s, sn_patch_e, 1, 1)
522 call bitarray_create(field%valid_mask, &
523 (we_mem_stag_e-we_mem_stag_s)+1, &
524 (sn_mem_e-sn_mem_s)+1)
525 do j=1,(sn_mem_e-sn_mem_s)+1
526 do i=1,(we_mem_stag_e-we_mem_stag_s)+1
527 call bitarray_set(field%valid_mask, i, j)
528 end do
529 end do
530 else if (field%map%stagger == V) then
531 allocate(field%r_arr(we_mem_s:we_mem_e,&
532 sn_mem_stag_s:sn_mem_stag_e))
533 field%r_arr(we_patch_s:we_patch_e,sn_patch_stag_s:sn_patch_stag_e) = real_array(sp1:ep1,sp2:ep2,k)
534 call exchange_halo_r(field%r_arr, &
535 we_mem_s, we_mem_e, sn_mem_stag_s, sn_mem_stag_e, 1, 1, &
536 we_patch_s, we_patch_e, sn_patch_stag_s, sn_patch_stag_e, 1, 1)
537 call bitarray_create(field%valid_mask, &
538 (we_mem_e-we_mem_s)+1, &
539 (sn_mem_stag_e-sn_mem_stag_s)+1)
540 do j=1,(sn_mem_stag_e-sn_mem_stag_s)+1
541 do i=1,(we_mem_e-we_mem_s)+1
542 call bitarray_set(field%valid_mask, i, j)
543 end do
544 end do
545 end if
547 nullify(field%modified_mask)
549 call storage_put_field(field)
551 end do
553 end if
554 end do
556 ! Done reading all static fields for this domain
557 call input_close()
559 end subroutine get_static_fields
562 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
563 ! Name: process_single_met_time
564 !
565 ! Purpose:
566 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
567 subroutine process_single_met_time(do_const_processing, &
568 temp_date, n, extra_row, extra_col, xlat, xlon, &
569 xlat_u, xlon_u, xlat_v, xlon_v, landmask, &
570 title, dyn_opt, &
571 west_east_dim, south_north_dim, &
572 we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
573 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
574 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
575 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
576 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
577 got_this_field, &
578 map_proj, mminlu, num_land_cat, is_water, is_lake, is_ice, is_urban, i_soilwater, &
579 grid_id, parent_id, i_parent_start, &
580 j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
581 cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
582 pole_lat, pole_lon, parent_grid_ratio, sub_x, sub_y, &
583 corner_lats, corner_lons, output_flags)
585 use bitarray_module
586 use gridinfo_module
587 use interp_module
588 use interp_option_module
589 use llxy_module
590 use misc_definitions_module
591 use module_debug
592 use output_module
593 use parallel_module
594 use read_met_module
595 use rotate_winds_module
596 use storage_module
598 implicit none
600 ! Arguments
601 logical, intent(in) :: do_const_processing
602 integer, intent(in) :: n, dyn_opt, west_east_dim, south_north_dim, map_proj, &
603 we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
604 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
605 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
606 we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
607 sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
608 is_water, is_lake, is_ice, is_urban, i_soilwater, &
609 grid_id, parent_id, i_parent_start, j_parent_start, &
610 i_parent_end, j_parent_end, parent_grid_ratio, sub_x, sub_y, num_land_cat
611 ! BUG: Should we be passing these around as pointers, or just declare them as arrays?
612 real, pointer, dimension(:,:) :: landmask
613 real, intent(in) :: dom_dx, dom_dy, cen_lat, moad_cen_lat, cen_lon, stand_lon, &
614 truelat1, truelat2, pole_lat, pole_lon
615 real, dimension(16), intent(in) :: corner_lats, corner_lons
616 real, pointer, dimension(:,:) :: xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v
617 logical, intent(in) :: extra_row, extra_col
618 logical, dimension(:), intent(inout) :: got_this_field
619 character (len=19), intent(in) :: temp_date
620 character (len=128), intent(in) :: mminlu
621 character (len=128), dimension(:), intent(inout) :: output_flags
623 ! BUG: Move this constant to misc_definitions_module?
624 integer, parameter :: BDR_WIDTH = 3
626 ! Local variables
627 integer :: istatus, iqstatus, fg_idx, idx, idxt, i, j, bottom_top_dim, &
628 sm1, em1, sm2, em2, sm3, em3, &
629 sp1, ep1, sp2, ep2, sp3, ep3, &
630 sd1, ed1, sd2, ed2, sd3, ed3, &
631 u_idx, bdr_wdth
632 integer :: nmet_flags
633 integer :: num_metgrid_soil_levs
634 integer, pointer, dimension(:) :: soil_levels
635 real :: rx, ry
636 real :: threshold
637 logical :: do_gcell_interp
638 integer, pointer, dimension(:) :: u_levels, v_levels
639 real, pointer, dimension(:,:) :: halo_slab
640 real, pointer, dimension(:,:,:) :: real_array
641 character (len=19) :: output_date
642 character (len=128) :: cname, title
643 character (len=MAX_FILENAME_LEN) :: input_name
644 character (len=128), allocatable, dimension(:) :: met_flags
645 type (fg_input) :: field, u_field, v_field
646 type (met_data) :: fg_data
649 ! For this time, we need to process all first-guess filename roots. When we
650 ! hit a root containing a '*', we assume we have hit the end of the list
651 fg_idx = 1
652 if (do_const_processing) then
653 input_name = constants_name(fg_idx)
654 else
655 input_name = fg_name(fg_idx)
656 end if
657 do while (input_name /= '*')
659 call mprintf(.true.,STDOUT, ' %s', s1=input_name)
660 call mprintf(.true.,LOGFILE, 'Getting input fields from %s', s1=input_name)
662 ! Do a first pass through this fg source to get all mask fields used
663 ! during interpolation
664 call get_interp_masks(trim(input_name), do_const_processing, temp_date)
666 istatus = 0
668 ! Initialize the module for reading in the met fields
669 call read_met_init(trim(input_name), do_const_processing, temp_date, istatus)
671 if (istatus == 0) then
673 ! Process all fields and levels from the current file; read_next_met_field()
674 ! will return a non-zero status when there are no more fields to be read.
675 do while (istatus == 0)
677 call read_next_met_field(fg_data, istatus)
679 if (istatus == 0) then
681 ! Find index into fieldname, interp_method, masked, and fill_missing
682 ! of the current field
683 idxt = num_entries + 1
684 do idx=1,num_entries
685 if ((index(fieldname(idx), trim(fg_data%field)) /= 0) .and. &
686 (len_trim(fieldname(idx)) == len_trim(fg_data%field))) then
688 got_this_field(idx) = .true.
690 if (index(input_name,trim(from_input(idx))) /= 0 .or. &
691 (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
692 idxt = idx
693 end if
695 end if
696 end do
697 idx = idxt
698 if (idx > num_entries) idx = num_entries ! The last entry is a default
700 ! Do we need to rename this field?
701 if (output_name(idx) /= ' ') then
702 fg_data%field = output_name(idx)(1:9)
704 idxt = num_entries + 1
705 do idx=1,num_entries
706 if ((index(fieldname(idx), trim(fg_data%field)) /= 0) .and. &
707 (len_trim(fieldname(idx)) == len_trim(fg_data%field))) then
709 got_this_field(idx) = .true.
711 if (index(input_name,trim(from_input(idx))) /= 0 .or. &
712 (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
713 idxt = idx
714 end if
716 end if
717 end do
718 idx = idxt
719 if (idx > num_entries) idx = num_entries ! The last entry is a default
720 end if
722 ! Do a simple check to see whether this is a global dataset
723 ! Note that we do not currently support regional Gaussian grids
724 if ((fg_data%iproj == PROJ_LATLON .and. abs(fg_data%nx * fg_data%deltalon - 360.) < 0.0001) &
725 .or. (fg_data%iproj == PROJ_GAUSS)) then
726 bdr_wdth = BDR_WIDTH
727 allocate(halo_slab(1-BDR_WIDTH:fg_data%nx+BDR_WIDTH,1:fg_data%ny))
729 halo_slab(1:fg_data%nx, 1:fg_data%ny) = &
730 fg_data%slab(1:fg_data%nx, 1:fg_data%ny)
732 halo_slab(1-BDR_WIDTH:0, 1:fg_data%ny) = &
733 fg_data%slab(fg_data%nx-BDR_WIDTH+1:fg_data%nx, 1:fg_data%ny)
735 halo_slab(fg_data%nx+1:fg_data%nx+BDR_WIDTH, 1:fg_data%ny) = &
736 fg_data%slab(1:BDR_WIDTH, 1:fg_data%ny)
738 deallocate(fg_data%slab)
739 else
740 bdr_wdth = 0
741 halo_slab => fg_data%slab
742 nullify(fg_data%slab)
743 end if
745 call mprintf(.true.,LOGFILE,'Processing %s at level %f.',s1=fg_data%field,f1=fg_data%xlvl)
747 call push_source_projection(fg_data%iproj, fg_data%xlonc, fg_data%truelat1, &
748 fg_data%truelat2, fg_data%dx, fg_data%dy, fg_data%deltalat, &
749 fg_data%deltalon, fg_data%starti, fg_data%startj, &
750 fg_data%startlat, fg_data%startlon, earth_radius=fg_data%earth_radius*1000.)
752 ! Initialize fg_input structure to store the field
753 field%header%version = 1
754 field%header%date = fg_data%hdate//' '
755 if (do_const_processing) then
756 field%header%time_dependent = .false.
757 else
758 field%header%time_dependent = .true.
759 end if
760 field%header%forecast_hour = fg_data%xfcst
761 field%header%fg_source = 'FG'
762 field%header%field = ' '
763 field%header%field(1:9) = fg_data%field
764 field%header%units = ' '
765 field%header%units(1:25) = fg_data%units
766 field%header%description = ' '
767 field%header%description(1:46) = fg_data%desc
768 call get_z_dim_name(fg_data%field,field%header%vertical_coord)
769 field%header%vertical_level = nint(fg_data%xlvl)
770 field%header%sr_x = 1
771 field%header%sr_y = 1
772 field%header%array_order = 'XY '
773 field%header%is_wind_grid_rel = fg_data%is_wind_grid_rel
774 field%header%array_has_missing_values = .false.
775 nullify(field%r_arr)
776 nullify(field%valid_mask)
777 nullify(field%modified_mask)
779 if (output_this_field(idx) .and. flag_in_output(idx) /= ' ') then
780 output_flags(idx) = flag_in_output(idx)
781 end if
783 ! If we should not output this field, just list it as a mask field
784 if (output_this_field(idx)) then
785 field%header%mask_field = .false.
786 else
787 field%header%mask_field = .true.
788 end if
790 !
791 ! Before actually doing any interpolation to the model grid, we must check
792 ! whether we will be using the average_gcell interpolator that averages all
793 ! source points in each model grid cell
794 !
795 do_gcell_interp = .false.
796 if (index(interp_method(idx),'average_gcell') /= 0) then
798 call get_gcell_threshold(interp_method(idx), threshold, istatus)
799 if (istatus == 0) then
800 if (fg_data%dx == 0. .and. fg_data%dy == 0. .and. &
801 fg_data%deltalat /= 0. .and. fg_data%deltalon /= 0.) then
802 fg_data%dx = abs(fg_data%deltalon)
803 fg_data%dy = abs(fg_data%deltalat)
804 else
805 ! BUG: Need to more correctly handle dx/dy in meters.
806 fg_data%dx = fg_data%dx / 111000. ! Convert meters to approximate degrees
807 fg_data%dy = fg_data%dy / 111000.
808 end if
809 if (gridtype == 'C') then
810 if (threshold*max(fg_data%dx,fg_data%dy)*111. <= max(dom_dx,dom_dy)/1000.) &
811 do_gcell_interp = .true.
812 else if (gridtype == 'E') then
813 if (threshold*max(fg_data%dx,fg_data%dy) <= max(dom_dx,dom_dy)) &
814 do_gcell_interp = .true.
815 end if
816 end if
817 end if
819 ! Interpolate to U staggering
820 if (output_stagger(idx) == U) then
822 call storage_query_field(field, iqstatus)
823 if (iqstatus == 0) then
824 call storage_get_field(field, iqstatus)
825 call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// &
826 ' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl))
827 if (associated(field%modified_mask)) then
828 call bitarray_destroy(field%modified_mask)
829 nullify(field%modified_mask)
830 end if
831 else
832 allocate(field%valid_mask)
833 call bitarray_create(field%valid_mask, we_mem_stag_e-we_mem_stag_s+1, sn_mem_e-sn_mem_s+1)
834 end if
836 ! Save a copy of the fg_input structure for the U field so that we can find it later
837 if (is_u_field(idx)) call dup(field, u_field)
839 allocate(field%modified_mask)
840 call bitarray_create(field%modified_mask, we_mem_stag_e-we_mem_stag_s+1, sn_mem_e-sn_mem_s+1)
842 if (do_const_processing .or. field%header%time_dependent) then
843 call interp_met_field(input_name, fg_data%field, U, M, &
844 field, xlat_u, xlon_u, we_mem_stag_s, we_mem_stag_e, sn_mem_s, sn_mem_e, &
845 halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
846 field%modified_mask)
847 else
848 call mprintf(.true.,INFORM,' - already processed this field from constant file.')
849 end if
851 ! Interpolate to V staggering
852 else if (output_stagger(idx) == V) then
854 call storage_query_field(field, iqstatus)
855 if (iqstatus == 0) then
856 call storage_get_field(field, iqstatus)
857 call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// &
858 ' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl))
859 if (associated(field%modified_mask)) then
860 call bitarray_destroy(field%modified_mask)
861 nullify(field%modified_mask)
862 end if
863 else
864 allocate(field%valid_mask)
865 call bitarray_create(field%valid_mask, we_mem_e-we_mem_s+1, sn_mem_stag_e-sn_mem_stag_s+1)
866 end if
868 ! Save a copy of the fg_input structure for the V field so that we can find it later
869 if (is_v_field(idx)) call dup(field, v_field)
871 allocate(field%modified_mask)
872 call bitarray_create(field%modified_mask, we_mem_e-we_mem_s+1, sn_mem_stag_e-sn_mem_stag_s+1)
874 if (do_const_processing .or. field%header%time_dependent) then
875 call interp_met_field(input_name, fg_data%field, V, M, &
876 field, xlat_v, xlon_v, we_mem_s, we_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
877 halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
878 field%modified_mask)
879 else
880 call mprintf(.true.,INFORM,' - already processed this field from constant file.')
881 end if
883 ! Interpolate to VV staggering
884 else if (output_stagger(idx) == VV) then
886 call storage_query_field(field, iqstatus)
887 if (iqstatus == 0) then
888 call storage_get_field(field, iqstatus)
889 call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// &
890 ' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl))
891 if (associated(field%modified_mask)) then
892 call bitarray_destroy(field%modified_mask)
893 nullify(field%modified_mask)
894 end if
895 else
896 allocate(field%valid_mask)
897 call bitarray_create(field%valid_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1)
898 end if
900 ! Save a copy of the fg_input structure for the U field so that we can find it later
901 if (is_u_field(idx)) call dup(field, u_field)
903 ! Save a copy of the fg_input structure for the V field so that we can find it later
904 if (is_v_field(idx)) call dup(field, v_field)
906 allocate(field%modified_mask)
907 call bitarray_create(field%modified_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1)
909 if (do_const_processing .or. field%header%time_dependent) then
910 call interp_met_field(input_name, fg_data%field, VV, M, &
911 field, xlat_v, xlon_v, we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
912 halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
913 field%modified_mask)
914 else
915 call mprintf(.true.,INFORM,' - already processed this field from constant file.')
916 end if
918 ! All other fields interpolated to M staggering for C grid, H staggering for E grid
919 else
921 call storage_query_field(field, iqstatus)
922 if (iqstatus == 0) then
923 call storage_get_field(field, iqstatus)
924 call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// &
925 ' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl))
926 if (associated(field%modified_mask)) then
927 call bitarray_destroy(field%modified_mask)
928 nullify(field%modified_mask)
929 end if
930 else
931 allocate(field%valid_mask)
932 call bitarray_create(field%valid_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1)
933 end if
935 allocate(field%modified_mask)
936 call bitarray_create(field%modified_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1)
938 if (do_const_processing .or. field%header%time_dependent) then
939 if (gridtype == 'C') then
940 call interp_met_field(input_name, fg_data%field, M, M, &
941 field, xlat, xlon, we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
942 halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
943 field%modified_mask, landmask)
945 else if (gridtype == 'E') then
946 call interp_met_field(input_name, fg_data%field, HH, M, &
947 field, xlat, xlon, we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
948 halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
949 field%modified_mask, landmask)
950 end if
951 else
952 call mprintf(.true.,INFORM,' - already processed this field from constant file.')
953 end if
955 end if
957 call bitarray_merge(field%valid_mask, field%modified_mask)
959 deallocate(halo_slab)
961 ! Store the interpolated field
962 call storage_put_field(field)
964 call pop_source_projection()
966 end if
967 end do
969 call read_met_close()
971 call push_source_projection(fg_data%iproj, fg_data%xlonc, fg_data%truelat1, &
972 fg_data%truelat2, fg_data%dx, fg_data%dy, fg_data%deltalat, &
973 fg_data%deltalon, fg_data%starti, fg_data%startj, &
974 fg_data%startlat, fg_data%startlon, earth_radius=fg_data%earth_radius*1000.)
976 !
977 ! If necessary, rotate winds to earth-relative for this fg source
978 !
980 call storage_get_levels(u_field, u_levels)
981 call storage_get_levels(v_field, v_levels)
983 if (associated(u_levels) .and. associated(v_levels)) then
984 u_idx = 1
985 do u_idx = 1, size(u_levels)
986 u_field%header%vertical_level = u_levels(u_idx)
987 call storage_get_field(u_field, istatus)
988 v_field%header%vertical_level = v_levels(u_idx)
989 call storage_get_field(v_field, istatus)
991 if (associated(u_field%modified_mask) .and. &
992 associated(v_field%modified_mask)) then
994 if (u_field%header%is_wind_grid_rel) then
995 if (gridtype == 'C') then
996 call map_to_met(u_field%r_arr, u_field%modified_mask, &
997 v_field%r_arr, v_field%modified_mask, &
998 we_mem_stag_s, sn_mem_s, &
999 we_mem_stag_e, sn_mem_e, &
1000 we_mem_s, sn_mem_stag_s, &
1001 we_mem_e, sn_mem_stag_e, &
1002 xlon_u, xlon_v, xlat_u, xlat_v)
1003 else if (gridtype == 'E') then
1004 call map_to_met_nmm(u_field%r_arr, u_field%modified_mask, &
1005 v_field%r_arr, v_field%modified_mask, &
1006 we_mem_s, sn_mem_s, &
1007 we_mem_e, sn_mem_e, &
1008 xlat_v, xlon_v)
1009 end if
1010 end if
1012 call bitarray_destroy(u_field%modified_mask)
1013 call bitarray_destroy(v_field%modified_mask)
1014 nullify(u_field%modified_mask)
1015 nullify(v_field%modified_mask)
1016 call storage_put_field(u_field)
1017 call storage_put_field(v_field)
1018 end if
1020 end do
1022 deallocate(u_levels)
1023 deallocate(v_levels)
1025 end if
1027 call pop_source_projection()
1029 else
1030 if (do_const_processing) then
1031 call mprintf(.true.,WARN,'Couldn''t open file %s for input.',s1=input_name)
1032 else
1033 call mprintf(.true.,WARN,'Couldn''t open file %s for input.',s1=trim(input_name)//':'//trim(temp_date))
1034 end if
1035 end if
1037 fg_idx = fg_idx + 1
1038 if (do_const_processing) then
1039 input_name = constants_name(fg_idx)
1040 else
1041 input_name = fg_name(fg_idx)
1042 end if
1043 end do ! while (input_name /= '*')
1045 !
1046 ! Rotate winds from earth-relative to grid-relative
1047 !
1049 call storage_get_levels(u_field, u_levels)
1050 call storage_get_levels(v_field, v_levels)
1052 if (associated(u_levels) .and. associated(v_levels)) then
1053 u_idx = 1
1054 do u_idx = 1, size(u_levels)
1055 u_field%header%vertical_level = u_levels(u_idx)
1056 call storage_get_field(u_field, istatus)
1057 v_field%header%vertical_level = v_levels(u_idx)
1058 call storage_get_field(v_field, istatus)
1060 if (gridtype == 'C') then
1061 call met_to_map(u_field%r_arr, u_field%valid_mask, &
1062 v_field%r_arr, v_field%valid_mask, &
1063 we_mem_stag_s, sn_mem_s, &
1064 we_mem_stag_e, sn_mem_e, &
1065 we_mem_s, sn_mem_stag_s, &
1066 we_mem_e, sn_mem_stag_e, &
1067 xlon_u, xlon_v, xlat_u, xlat_v)
1068 else if (gridtype == 'E') then
1069 call met_to_map_nmm(u_field%r_arr, u_field%valid_mask, &
1070 v_field%r_arr, v_field%valid_mask, &
1071 we_mem_s, sn_mem_s, &
1072 we_mem_e, sn_mem_e, &
1073 xlat_v, xlon_v)
1074 end if
1076 end do
1078 deallocate(u_levels)
1079 deallocate(v_levels)
1081 end if
1083 if (do_const_processing) return
1085 !
1086 ! Now that we have all degribbed fields, we build a 3-d pressure field, and fill in any
1087 ! missing levels in the other 3-d fields
1088 !
1089 call mprintf(.true.,LOGFILE,'Filling missing levels.')
1090 call fill_missing_levels(output_flags)
1092 call mprintf(.true.,LOGFILE,'Creating derived fields.')
1093 call create_derived_fields(gridtype, fg_data%hdate, fg_data%xfcst, &
1094 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
1095 we_mem_stag_s, we_mem_stag_e, sn_mem_stag_s, sn_mem_stag_e, &
1096 got_this_field, output_flags)
1098 !
1099 ! Check that every mandatory field was found in input data
1100 !
1101 do i=1,num_entries
1102 if (is_mandatory(i) .and. .not. got_this_field(i)) then
1103 call mprintf(.true.,ERROR,'The mandatory field %s was not found in any input data.',s1=fieldname(i))
1104 end if
1105 end do
1107 !
1108 ! Before we begin to write fields, if debug_level is set high enough, we
1109 ! write a table of which fields are available at which levels to the
1110 ! metgrid.log file, and then we check to see if any fields are not
1111 ! completely covered with data.
1112 !
1113 call storage_print_fields()
1114 call find_missing_values()
1116 !
1117 ! All of the processing is now done for this time period for this domain;
1118 ! now we simply output every field from the storage module.
1119 !
1121 title = 'OUTPUT FROM METGRID V3.2.1'
1123 ! Initialize the output module for this domain and time
1124 call mprintf(.true.,LOGFILE,'Initializing output module.')
1125 output_date = temp_date
1126 if ( .not. nocolons ) then
1127 if (len_trim(temp_date) == 13) then
1128 output_date(14:19) = ':00:00'
1129 else if (len_trim(temp_date) == 16) then
1130 output_date(17:19) = ':00'
1131 end if
1132 else
1133 if (len_trim(temp_date) == 13) then
1134 output_date(14:19) = '_00_00'
1135 else if (len_trim(temp_date) == 16) then
1136 output_date(17:19) = '_00'
1137 end if
1138 endif
1140 call output_init(n, title, output_date, gridtype, dyn_opt, &
1141 corner_lats, corner_lons, &
1142 we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
1143 we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, &
1144 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
1145 extra_col, extra_row)
1147 call get_bottom_top_dim(bottom_top_dim)
1149 ! Add in a flag to tell real that we have seven new msf fields
1150 nmet_flags = num_entries + 1
1151 allocate(met_flags(nmet_flags))
1152 do i=1,num_entries
1153 met_flags(i) = output_flags(i)
1154 end do
1155 if (gridtype == 'C') then
1156 met_flags(num_entries+1) = 'FLAG_MF_XY'
1157 else
1158 met_flags(num_entries+1) = ' '
1159 end if
1161 ! Find out how many soil levels or layers we have; this assumes a field named ST
1162 field % header % field = 'ST'
1163 nullify(soil_levels)
1164 call storage_get_levels(field, soil_levels)
1166 if (.not. associated(soil_levels)) then
1167 field % header % field = 'SOILT'
1168 nullify(soil_levels)
1169 call storage_get_levels(field, soil_levels)
1170 end if
1172 if (.not. associated(soil_levels)) then
1173 field % header % field = 'STC_WPS'
1174 nullify(soil_levels)
1175 call storage_get_levels(field, soil_levels)
1176 end if
1178 if (associated(soil_levels)) then
1179 num_metgrid_soil_levs = size(soil_levels)
1180 deallocate(soil_levels)
1181 else
1182 num_metgrid_soil_levs = 0
1183 end if
1185 ! First write out global attributes
1186 call mprintf(.true.,LOGFILE,'Writing global attributes to output.')
1187 call write_global_attrs(title, output_date, gridtype, dyn_opt, west_east_dim, &
1188 south_north_dim, bottom_top_dim, &
1189 we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
1190 sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
1191 map_proj, mminlu, num_land_cat, &
1192 is_water, is_lake, is_ice, is_urban, i_soilwater, &
1193 grid_id, parent_id, i_parent_start, &
1194 j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
1195 cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
1196 pole_lat, pole_lon, parent_grid_ratio, sub_x, sub_y, &
1197 corner_lats, corner_lons, num_metgrid_soil_levs, met_flags, nmet_flags)
1199 deallocate(met_flags)
1201 call reset_next_field()
1203 istatus = 0
1205 ! Now loop over all output fields, writing each to the output module
1206 do while (istatus == 0)
1207 call get_next_output_field(cname, real_array, &
1208 sm1, em1, sm2, em2, sm3, em3, istatus)
1209 if (istatus == 0) then
1211 call mprintf(.true.,LOGFILE,'Writing field %s to output.',s1=cname)
1212 call write_field(sm1, em1, sm2, em2, sm3, em3, &
1213 cname, output_date, real_array)
1214 deallocate(real_array)
1216 end if
1218 end do
1220 call mprintf(.true.,LOGFILE,'Closing output file.')
1221 call output_close()
1223 ! Free up memory used by met fields for this valid time
1224 call storage_delete_all_td()
1226 end subroutine process_single_met_time
1229 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1230 ! Name: get_interp_masks
1231 !
1232 ! Purpose:
1233 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1234 subroutine get_interp_masks(fg_prefix, is_constants, fg_date)
1236 use interp_option_module
1237 use read_met_module
1238 use storage_module
1240 implicit none
1242 ! Arguments
1243 logical, intent(in) :: is_constants
1244 character (len=*), intent(in) :: fg_prefix, fg_date
1246 ! BUG: Move this constant to misc_definitions_module?
1247 integer, parameter :: BDR_WIDTH = 3
1249 ! Local variables
1250 integer :: i, istatus, idx, idxt
1251 type (fg_input) :: mask_field
1252 type (met_data) :: fg_data
1254 istatus = 0
1256 call read_met_init(fg_prefix, is_constants, fg_date, istatus)
1258 do while (istatus == 0)
1260 call read_next_met_field(fg_data, istatus)
1262 if (istatus == 0) then
1264 ! Find out which METGRID.TBL entry goes with this field
1265 idxt = num_entries + 1
1266 do idx=1,num_entries
1267 if ((index(fieldname(idx), trim(fg_data%field)) /= 0) .and. &
1268 (len_trim(fieldname(idx)) == len_trim(fg_data%field))) then
1270 if (index(fg_prefix,trim(from_input(idx))) /= 0 .or. &
1271 (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
1272 idxt = idx
1273 end if
1275 end if
1276 end do
1277 idx = idxt
1278 if (idx > num_entries) idx = num_entries ! The last entry is a default
1280 ! Do we need to rename this field?
1281 if (output_name(idx) /= ' ') then
1282 fg_data%field = output_name(idx)(1:9)
1284 idxt = num_entries + 1
1285 do idx=1,num_entries
1286 if ((index(fieldname(idx), trim(fg_data%field)) /= 0) .and. &
1287 (len_trim(fieldname(idx)) == len_trim(fg_data%field))) then
1289 if (index(fg_prefix,trim(from_input(idx))) /= 0 .or. &
1290 (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
1291 idxt = idx
1292 end if
1294 end if
1295 end do
1296 idx = idxt
1297 if (idx > num_entries) idx = num_entries ! The last entry is a default
1298 end if
1300 do i=1,num_entries
1301 if (interp_mask(i) /= ' ' .and. (trim(interp_mask(i)) == trim(fg_data%field))) then
1303 mask_field%header%version = 1
1304 mask_field%header%date = ' '
1305 mask_field%header%date = fg_date
1306 if (is_constants) then
1307 mask_field%header%time_dependent = .false.
1308 else
1309 mask_field%header%time_dependent = .true.
1310 end if
1311 mask_field%header%mask_field = .true.
1312 mask_field%header%forecast_hour = 0.
1313 mask_field%header%fg_source = 'degribbed met data'
1314 mask_field%header%field = trim(fg_data%field)//'.mask'
1315 mask_field%header%units = '-'
1316 mask_field%header%description = '-'
1317 mask_field%header%vertical_coord = 'none'
1318 mask_field%header%vertical_level = 1
1319 mask_field%header%sr_x = 1
1320 mask_field%header%sr_y = 1
1321 mask_field%header%array_order = 'XY'
1322 mask_field%header%dim1(1) = 1
1323 mask_field%header%dim1(2) = fg_data%nx
1324 mask_field%header%dim2(1) = 1
1325 mask_field%header%dim2(2) = fg_data%ny
1326 mask_field%header%is_wind_grid_rel = .true.
1327 mask_field%header%array_has_missing_values = .false.
1328 mask_field%map%stagger = M
1330 ! Do a simple check to see whether this is a global lat/lon dataset
1331 ! Note that we do not currently support regional Gaussian grids
1332 if ((fg_data%iproj == PROJ_LATLON .and. abs(fg_data%nx * fg_data%deltalon - 360.) < 0.0001) &
1333 .or. (fg_data%iproj == PROJ_GAUSS)) then
1334 allocate(mask_field%r_arr(1-BDR_WIDTH:fg_data%nx+BDR_WIDTH,1:fg_data%ny))
1336 mask_field%r_arr(1:fg_data%nx, 1:fg_data%ny) = &
1337 fg_data%slab(1:fg_data%nx, 1:fg_data%ny)
1339 mask_field%r_arr(1-BDR_WIDTH:0, 1:fg_data%ny) = &
1340 fg_data%slab(fg_data%nx-BDR_WIDTH+1:fg_data%nx, 1:fg_data%ny)
1342 mask_field%r_arr(fg_data%nx+1:fg_data%nx+BDR_WIDTH, 1:fg_data%ny) = &
1343 fg_data%slab(1:BDR_WIDTH, 1:fg_data%ny)
1345 else
1346 allocate(mask_field%r_arr(1:fg_data%nx,1:fg_data%ny))
1347 mask_field%r_arr = fg_data%slab
1348 end if
1350 nullify(mask_field%valid_mask)
1351 nullify(mask_field%modified_mask)
1353 call storage_put_field(mask_field)
1355 exit
1357 end if
1358 end do
1360 if (associated(fg_data%slab)) deallocate(fg_data%slab)
1362 end if
1364 end do
1366 call read_met_close()
1368 end subroutine get_interp_masks
1371 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1372 ! Name: interp_met_field
1373 !
1374 ! Purpose:
1375 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1376 subroutine interp_met_field(input_name, short_fieldnm, ifieldstagger, istagger, &
1377 field, xlat, xlon, sm1, em1, sm2, em2, &
1378 slab, minx, maxx, miny, maxy, bdr, do_gcell_interp, &
1379 new_pts, landmask)
1381 use bitarray_module
1382 use interp_module
1383 use interp_option_module
1384 use llxy_module
1385 use misc_definitions_module
1386 use storage_module
1388 implicit none
1390 ! Arguments
1391 integer, intent(in) :: ifieldstagger, istagger, &
1392 sm1, em1, sm2, em2, &
1393 minx, maxx, miny, maxy, bdr
1394 real, dimension(minx:maxx,miny:maxy), intent(in) :: slab
1395 real, dimension(sm1:em1,sm2:em2), intent(in) :: xlat, xlon
1396 real, dimension(sm1:em1,sm2:em2), intent(in), optional :: landmask
1397 logical, intent(in) :: do_gcell_interp
1398 character (len=9), intent(in) :: short_fieldnm
1399 character (len=MAX_FILENAME_LEN), intent(in) :: input_name
1400 type (fg_input), intent(inout) :: field
1401 type (bitarray), intent(inout) :: new_pts
1403 ! Local variables
1404 integer :: i, j, idx, idxt, orig_selected_proj, interp_mask_status, &
1405 interp_land_mask_status, interp_water_mask_status
1406 integer, pointer, dimension(:) :: interp_array
1407 real :: rx, ry, temp
1408 real, pointer, dimension(:,:) :: data_count
1409 type (fg_input) :: mask_field, mask_water_field, mask_land_field
1411 ! Find index into fieldname, interp_method, masked, and fill_missing
1412 ! of the current field
1413 idxt = num_entries + 1
1414 do idx=1,num_entries
1415 if ((index(fieldname(idx), trim(short_fieldnm)) /= 0) .and. &
1416 (len_trim(fieldname(idx)) == len_trim(short_fieldnm))) then
1417 if (index(input_name,trim(from_input(idx))) /= 0 .or. &
1418 (from_input(idx) == '*' .and. idxt == num_entries + 1)) then
1419 idxt = idx
1420 end if
1421 end if
1422 end do
1423 idx = idxt
1424 if (idx > num_entries) then
1425 call mprintf(.true.,WARN,'Entry in METGRID.TBL not found for field %s. '// &
1426 'Default options will be used for this field!', s1=short_fieldnm)
1427 idx = num_entries ! The last entry is a default
1428 end if
1430 field%header%dim1(1) = sm1
1431 field%header%dim1(2) = em1
1432 field%header%dim2(1) = sm2
1433 field%header%dim2(2) = em2
1434 field%map%stagger = ifieldstagger
1435 if (.not. associated(field%r_arr)) then
1436 allocate(field%r_arr(sm1:em1,sm2:em2))
1437 end if
1439 interp_mask_status = 1
1440 interp_land_mask_status = 1
1441 interp_water_mask_status = 1
1443 if (interp_mask(idx) /= ' ') then
1444 mask_field%header%version = 1
1445 mask_field%header%forecast_hour = 0.
1446 mask_field%header%field = trim(interp_mask(idx))//'.mask'
1447 mask_field%header%vertical_coord = 'none'
1448 mask_field%header%vertical_level = 1
1450 call storage_get_field(mask_field, interp_mask_status)
1452 end if
1453 if (interp_land_mask(idx) /= ' ') then
1454 mask_land_field%header%version = 1
1455 mask_land_field%header%forecast_hour = 0.
1456 mask_land_field%header%field = trim(interp_land_mask(idx))//'.mask'
1457 mask_land_field%header%vertical_coord = 'none'
1458 mask_land_field%header%vertical_level = 1
1460 call storage_get_field(mask_land_field, interp_land_mask_status)
1462 end if
1463 if (interp_water_mask(idx) /= ' ') then
1464 mask_water_field%header%version = 1
1465 mask_water_field%header%forecast_hour = 0.
1466 mask_water_field%header%field = trim(interp_water_mask(idx))//'.mask'
1467 mask_water_field%header%vertical_coord = 'none'
1468 mask_water_field%header%vertical_level = 1
1470 call storage_get_field(mask_water_field, interp_water_mask_status)
1472 end if
1474 interp_array => interp_array_from_string(interp_method(idx))
1477 !
1478 ! Interpolate using average_gcell interpolation method
1479 !
1480 if (do_gcell_interp) then
1481 allocate(data_count(sm1:em1,sm2:em2))
1482 data_count = 0.
1484 if (interp_mask_status == 0) then
1485 call accum_continuous(slab, &
1486 minx, maxx, miny, maxy, 1, 1, bdr, &
1487 field%r_arr, data_count, &
1488 sm1, em1, sm2, em2, 1, 1, &
1489 istagger, &
1490 new_pts, missing_value(idx), interp_mask_val(idx), mask_field%r_arr)
1491 else
1492 call accum_continuous(slab, &
1493 minx, maxx, miny, maxy, 1, 1, bdr, &
1494 field%r_arr, data_count, &
1495 sm1, em1, sm2, em2, 1, 1, &
1496 istagger, &
1497 new_pts, missing_value(idx), -1.) ! The -1 is the maskval, but since we
1498 ! we do not give an optional mask, no
1499 ! no need to worry about -1s in data
1500 end if
1502 orig_selected_proj = iget_selected_domain()
1503 call select_domain(SOURCE_PROJ)
1504 do j=sm2,em2
1505 do i=sm1,em1
1507 if (present(landmask)) then
1509 if (landmask(i,j) /= masked(idx)) then
1510 if (data_count(i,j) > 0.) then
1511 field%r_arr(i,j) = field%r_arr(i,j) / data_count(i,j)
1512 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1513 else
1515 if (interp_mask_status == 0) then
1516 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1517 minx, maxx, miny, maxy, bdr, missing_value(idx), &
1518 mask_val=interp_mask_val(idx), mask_field=mask_field%r_arr)
1519 else
1520 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1521 minx, maxx, miny, maxy, bdr, missing_value(idx))
1522 end if
1524 if (temp /= missing_value(idx)) then
1525 field%r_arr(i,j) = temp
1526 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1527 end if
1529 end if
1530 else
1531 field%r_arr(i,j) = fill_missing(idx)
1532 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1533 end if
1535 if (.not. bitarray_test(new_pts, i-sm1+1, j-sm2+1) .and. &
1536 .not. bitarray_test(field%valid_mask, i-sm1+1, j-sm2+1)) then
1537 field%r_arr(i,j) = fill_missing(idx)
1539 ! Assume that if missing fill value is other than default, then user has asked
1540 ! to fill in any missing values, and we can consider this point to have
1541 ! received a valid value
1542 if (fill_missing(idx) /= NAN) call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1543 end if
1545 else
1547 if (data_count(i,j) > 0.) then
1548 field%r_arr(i,j) = field%r_arr(i,j) / data_count(i,j)
1549 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1550 else
1552 if (interp_mask_status == 0) then
1553 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1554 minx, maxx, miny, maxy, bdr, missing_value(idx), &
1555 mask_val=interp_mask_val(idx), mask_field=mask_field%r_arr)
1556 else
1557 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1558 minx, maxx, miny, maxy, bdr, missing_value(idx))
1559 end if
1561 if (temp /= missing_value(idx)) then
1562 field%r_arr(i,j) = temp
1563 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1564 end if
1566 end if
1568 if (.not. bitarray_test(new_pts, i-sm1+1, j-sm2+1) .and. &
1569 .not. bitarray_test(field%valid_mask, i-sm1+1, j-sm2+1)) then
1570 field%r_arr(i,j) = fill_missing(idx)
1572 ! Assume that if missing fill value is other than default, then user has asked
1573 ! to fill in any missing values, and we can consider this point to have
1574 ! received a valid value
1575 if (fill_missing(idx) /= NAN) call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1576 end if
1578 end if
1580 end do
1581 end do
1582 call select_domain(orig_selected_proj)
1583 deallocate(data_count)
1585 !
1586 ! No average_gcell interpolation method
1587 !
1588 else
1590 orig_selected_proj = iget_selected_domain()
1591 call select_domain(SOURCE_PROJ)
1592 do j=sm2,em2
1593 do i=sm1,em1
1594 if (present(landmask)) then
1596 if (masked(idx) == MASKED_BOTH) then
1598 if (landmask(i,j) == 0) then ! WATER POINT
1600 if (interp_land_mask_status == 0) then
1601 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1602 minx, maxx, miny, maxy, bdr, missing_value(idx), &
1603 mask_val=interp_land_mask_val(idx), mask_field=mask_land_field%r_arr)
1604 else
1605 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1606 minx, maxx, miny, maxy, bdr, missing_value(idx))
1607 end if
1609 else if (landmask(i,j) == 1) then ! LAND POINT
1611 if (interp_water_mask_status == 0) then
1612 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1613 minx, maxx, miny, maxy, bdr, missing_value(idx), &
1614 mask_val=interp_water_mask_val(idx), &
1615 mask_field=mask_water_field%r_arr)
1616 else
1617 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1618 minx, maxx, miny, maxy, bdr, missing_value(idx))
1619 end if
1621 end if
1623 else if (landmask(i,j) /= masked(idx)) then
1625 if (interp_mask_status == 0) then
1626 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1627 minx, maxx, miny, maxy, bdr, missing_value(idx), &
1628 mask_val=interp_mask_val(idx), mask_field=mask_field%r_arr)
1629 else
1630 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1631 minx, maxx, miny, maxy, bdr, missing_value(idx))
1632 end if
1634 else
1635 temp = missing_value(idx)
1636 end if
1638 ! No landmask for this field
1639 else
1641 if (interp_mask_status == 0) then
1642 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1643 minx, maxx, miny, maxy, bdr, missing_value(idx), &
1644 mask_val=interp_mask_val(idx), mask_field=mask_field%r_arr)
1645 else
1646 temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
1647 minx, maxx, miny, maxy, bdr, missing_value(idx))
1648 end if
1650 end if
1652 if (temp /= missing_value(idx)) then
1653 field%r_arr(i,j) = temp
1654 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1655 else if (present(landmask)) then
1656 if (landmask(i,j) == masked(idx)) then
1657 field%r_arr(i,j) = fill_missing(idx)
1658 call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1659 end if
1660 end if
1662 if (.not. bitarray_test(new_pts, i-sm1+1, j-sm2+1) .and. &
1663 .not. bitarray_test(field%valid_mask, i-sm1+1, j-sm2+1)) then
1664 field%r_arr(i,j) = fill_missing(idx)
1666 ! Assume that if missing fill value is other than default, then user has asked
1667 ! to fill in any missing values, and we can consider this point to have
1668 ! received a valid value
1669 if (fill_missing(idx) /= NAN) call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
1670 end if
1672 end do
1673 end do
1674 call select_domain(orig_selected_proj)
1675 end if
1677 deallocate(interp_array)
1679 end subroutine interp_met_field
1682 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1683 ! Name: interp_to_latlon
1684 !
1685 ! Purpose:
1686 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1687 function interp_to_latlon(rlat, rlon, istagger, interp_method_list, slab, &
1688 minx, maxx, miny, maxy, bdr, source_missing_value, &
1689 mask_field, mask_val)
1691 use interp_module
1692 use llxy_module
1694 implicit none
1696 ! Arguments
1697 integer, intent(in) :: minx, maxx, miny, maxy, bdr, istagger
1698 integer, dimension(:), intent(in) :: interp_method_list
1699 real, intent(in) :: rlat, rlon, source_missing_value
1700 real, dimension(minx:maxx,miny:maxy), intent(in) :: slab
1701 real, intent(in), optional :: mask_val
1702 real, dimension(minx:maxx,miny:maxy), intent(in), optional :: mask_field
1704 ! Return value
1705 real :: interp_to_latlon
1707 ! Local variables
1708 real :: rx, ry
1710 interp_to_latlon = source_missing_value
1712 call lltoxy(rlat, rlon, rx, ry, istagger)
1713 if (rx >= minx+bdr-0.5 .and. rx <= maxx-bdr+0.5) then
1714 if (present(mask_field) .and. present(mask_val)) then
1715 interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, 1, 1, source_missing_value, &
1716 interp_method_list, 1, mask_val, mask_field)
1717 else
1718 interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, 1, 1, source_missing_value, &
1719 interp_method_list, 1)
1720 end if
1721 else
1722 interp_to_latlon = source_missing_value
1723 end if
1725 if (interp_to_latlon == source_missing_value) then
1727 ! Try a lon in the range 0. to 360.; all lons in the xlon
1728 ! array should be in the range -180. to 180.
1729 if (rlon < 0.) then
1730 call lltoxy(rlat, rlon+360., rx, ry, istagger)
1731 if (rx >= minx+bdr-0.5 .and. rx <= maxx-bdr+0.5) then
1732 if (present(mask_field) .and. present(mask_val)) then
1733 interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, &
1734 1, 1, source_missing_value, &
1735 interp_method_list, 1, mask_val, mask_field)
1736 else
1737 interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, &
1738 1, 1, source_missing_value, &
1739 interp_method_list, 1)
1740 end if
1741 else
1742 interp_to_latlon = source_missing_value
1743 end if
1745 end if
1747 end if
1749 return
1751 end function interp_to_latlon
1754 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1755 ! Name: get_bottom_top_dim
1756 !
1757 ! Purpose:
1758 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1759 subroutine get_bottom_top_dim(bottom_top_dim)
1761 use interp_option_module
1762 use list_module
1763 use storage_module
1765 implicit none
1767 ! Arguments
1768 integer, intent(out) :: bottom_top_dim
1770 ! Local variables
1771 integer :: i, j
1772 integer, pointer, dimension(:) :: field_levels
1773 character (len=32) :: z_dim
1774 type (fg_input), pointer, dimension(:) :: headers
1775 type (list) :: temp_levels
1777 ! Initialize a list to store levels that are found for 3-d fields
1778 call list_init(temp_levels)
1780 ! Get a list of all time-dependent fields (given by their headers) from
1781 ! the storage module
1782 call storage_get_td_headers(headers)
1784 !
1785 ! Given headers of all fields, we first build a list of all possible levels
1786 ! for 3-d met fields (excluding sea-level, though).
1787 !
1788 do i=1,size(headers)
1789 call get_z_dim_name(headers(i)%header%field, z_dim)
1791 ! We only want to consider 3-d met fields
1792 if (z_dim(1:18) == 'num_metgrid_levels') then
1794 ! Find out what levels the current field has
1795 call storage_get_levels(headers(i), field_levels)
1796 do j=1,size(field_levels)
1798 ! If this level has not yet been encountered, add it to our list
1799 if (.not. list_search(temp_levels, ikey=field_levels(j), ivalue=field_levels(j))) then
1800 if (field_levels(j) /= 201300) then
1801 call list_insert(temp_levels, ikey=field_levels(j), ivalue=field_levels(j))
1802 end if
1803 end if
1805 end do
1807 deallocate(field_levels)
1809 end if
1811 end do
1813 bottom_top_dim = list_length(temp_levels)
1815 call list_destroy(temp_levels)
1816 deallocate(headers)
1818 end subroutine get_bottom_top_dim
1821 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1822 ! Name: fill_missing_levels
1823 !
1824 ! Purpose:
1825 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1826 subroutine fill_missing_levels(output_flags)
1828 use interp_option_module
1829 use list_module
1830 use module_debug
1831 use module_mergesort
1832 use storage_module
1834 implicit none
1836 ! Arguments
1837 character (len=128), dimension(:), intent(inout) :: output_flags
1839 ! Local variables
1840 integer :: i, ii, j, ix, jx, k, lower, upper, temp, istatus
1841 integer, pointer, dimension(:) :: union_levels, field_levels
1842 real, pointer, dimension(:) :: r_union_levels
1843 character (len=128) :: clevel
1844 type (fg_input) :: lower_field, upper_field, new_field, search_field
1845 type (fg_input), pointer, dimension(:) :: headers, all_headers
1846 type (list) :: temp_levels
1847 type (list_item), pointer, dimension(:) :: keys
1849 ! Initialize a list to store levels that are found for 3-d fields
1850 call list_init(temp_levels)
1852 ! Get a list of all fields (given by their headers) from the storage module
1853 call storage_get_td_headers(headers)
1854 call storage_get_all_headers(all_headers)
1856 !
1857 ! Given headers of all fields, we first build a list of all possible levels
1858 ! for 3-d met fields (excluding sea-level, though).
1859 !
1860 do i=1,size(headers)
1862 ! Find out what levels the current field has
1863 call storage_get_levels(headers(i), field_levels)
1864 do j=1,size(field_levels)
1866 ! If this level has not yet been encountered, add it to our list
1867 if (.not. list_search(temp_levels, ikey=field_levels(j), ivalue=field_levels(j))) then
1868 if (field_levels(j) /= 201300) then
1869 call list_insert(temp_levels, ikey=field_levels(j), ivalue=field_levels(j))
1870 end if
1871 end if
1873 end do
1875 deallocate(field_levels)
1877 end do
1879 if (list_length(temp_levels) > 0) then
1881 !
1882 ! With all possible levels stored in a list, get an array of levels, sorted
1883 ! in decreasing order
1884 !
1885 i = 0
1886 allocate(union_levels(list_length(temp_levels)))
1887 do while (list_length(temp_levels) > 0)
1888 i = i + 1
1889 call list_get_first_item(temp_levels, ikey=union_levels(i), ivalue=temp)
1890 end do
1891 call mergesort(union_levels, 1, size(union_levels))
1893 allocate(r_union_levels(size(union_levels)))
1894 do i=1,size(union_levels)
1895 r_union_levels(i) = real(union_levels(i))
1896 end do
1898 !
1899 ! With a sorted, complete list of levels, we need
1900 ! to go back and fill in missing levels for each 3-d field
1901 !
1902 do i=1,size(headers)
1904 !
1905 ! Find entry in METGRID.TBL for this field, if one exists; if it does, then the
1906 ! entry may tell us how to get values for the current field at the missing level
1907 !
1908 do ii=1,num_entries
1909 if (fieldname(ii) == headers(i)%header%field) exit
1910 end do
1911 if (ii <= num_entries) then
1912 call dup(headers(i),new_field)
1913 nullify(new_field%valid_mask)
1914 nullify(new_field%modified_mask)
1915 call fill_field(new_field, ii, output_flags, r_union_levels)
1916 end if
1918 end do
1920 deallocate(union_levels)
1921 deallocate(r_union_levels)
1922 deallocate(headers)
1924 call storage_get_td_headers(headers)
1926 !
1927 ! Now we may need to vertically interpolate to missing values in 3-d fields
1928 !
1929 do i=1,size(headers)
1931 call storage_get_levels(headers(i), field_levels)
1933 ! If this isn't a 3-d array, nothing to do
1934 if (size(field_levels) > 1) then
1936 do k=1,size(field_levels)
1937 call dup(headers(i),search_field)
1938 search_field%header%vertical_level = field_levels(k)
1939 call storage_get_field(search_field,istatus)
1940 if (istatus == 0) then
1941 JLOOP: do jx=search_field%header%dim2(1),search_field%header%dim2(2)
1942 ILOOP: do ix=search_field%header%dim1(1),search_field%header%dim1(2)
1943 if (.not. bitarray_test(search_field%valid_mask, &
1944 ix-search_field%header%dim1(1)+1, &
1945 jx-search_field%header%dim2(1)+1)) then
1947 call dup(search_field, lower_field)
1948 do lower=k-1,1,-1
1949 lower_field%header%vertical_level = field_levels(lower)
1950 call storage_get_field(lower_field,istatus)
1951 if (bitarray_test(lower_field%valid_mask, &
1952 ix-search_field%header%dim1(1)+1, &
1953 jx-search_field%header%dim2(1)+1)) &
1954 exit
1956 end do
1958 call dup(search_field, upper_field)
1959 do upper=k+1,size(field_levels)
1960 upper_field%header%vertical_level = field_levels(upper)
1961 call storage_get_field(upper_field,istatus)
1962 if (bitarray_test(upper_field%valid_mask, &
1963 ix-search_field%header%dim1(1)+1, &
1964 jx-search_field%header%dim2(1)+1)) &
1965 exit
1967 end do
1968 if (upper <= size(field_levels) .and. lower >= 1) then
1969 search_field%r_arr(ix,jx) = real(abs(field_levels(upper)-field_levels(k))) &
1970 / real(abs(field_levels(upper)-field_levels(lower))) &
1971 * lower_field%r_arr(ix,jx) &
1972 + real(abs(field_levels(k)-field_levels(lower))) &
1973 / real(abs(field_levels(upper)-field_levels(lower))) &
1974 * upper_field%r_arr(ix,jx)
1975 call bitarray_set(search_field%valid_mask, &
1976 ix-search_field%header%dim1(1)+1, &
1977 jx-search_field%header%dim2(1)+1)
1978 end if
1979 end if
1980 end do ILOOP
1981 end do JLOOP
1982 else
1983 call mprintf(.true.,ERROR, &
1984 'This is bad, could not get %s at level %i.', &
1985 s1=trim(search_field%header%field), i1=field_levels(k))
1986 end if
1987 end do
1989 end if
1991 deallocate(field_levels)
1993 end do
1995 end if
1997 call list_destroy(temp_levels)
1998 deallocate(all_headers)
1999 deallocate(headers)
2001 end subroutine fill_missing_levels
2004 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2005 ! Name: create_derived_fields
2006 !
2007 ! Purpose:
2008 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2009 subroutine create_derived_fields(arg_gridtype, hdate, xfcst, &
2010 we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
2011 we_mem_stag_s, we_mem_stag_e, sn_mem_stag_s, sn_mem_stag_e, &
2012 created_this_field, output_flags)
2014 use interp_option_module
2015 use list_module
2016 use module_mergesort
2017 use storage_module
2019 implicit none
2021 ! Arguments
2022 integer, intent(in) :: we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
2023 we_mem_stag_s, we_mem_stag_e, sn_mem_stag_s, sn_mem_stag_e
2024 real, intent(in) :: xfcst
2025 logical, dimension(:), intent(inout) :: created_this_field
2026 character (len=1), intent(in) :: arg_gridtype
2027 character (len=24), intent(in) :: hdate
2028 character (len=128), dimension(:), intent(inout) :: output_flags
2030 ! Local variables
2031 integer :: idx, i, j, istatus
2032 type (fg_input) :: field
2034 ! Initialize fg_input structure to store the field
2035 field%header%version = 1
2036 field%header%date = hdate//' '
2037 field%header%time_dependent = .true.
2038 field%header%mask_field = .false.
2039 field%header%forecast_hour = xfcst
2040 field%header%fg_source = 'Derived from FG'
2041 field%header%field = ' '
2042 field%header%units = ' '
2043 field%header%description = ' '
2044 field%header%vertical_level = 0
2045 field%header%sr_x = 1
2046 field%header%sr_y = 1
2047 field%header%array_order = 'XY '
2048 field%header%is_wind_grid_rel = .true.
2049 field%header%array_has_missing_values = .false.
2050 nullify(field%r_arr)
2051 nullify(field%valid_mask)
2052 nullify(field%modified_mask)
2054 !
2055 ! Check each entry in METGRID.TBL to see whether it is a derive field
2056 !
2057 do idx=1,num_entries
2058 if (is_derived_field(idx)) then
2060 created_this_field(idx) = .true.
2062 call mprintf(.true.,INFORM,'Going to create the field %s',s1=fieldname(idx))
2064 ! Intialize more fields in storage structure
2065 field%header%field = fieldname(idx)
2066 call get_z_dim_name(fieldname(idx),field%header%vertical_coord)
2067 field%map%stagger = output_stagger(idx)
2068 if (arg_gridtype == 'E') then
2069 field%header%dim1(1) = we_mem_s
2070 field%header%dim1(2) = we_mem_e
2071 field%header%dim2(1) = sn_mem_s
2072 field%header%dim2(2) = sn_mem_e
2073 else if (arg_gridtype == 'C') then
2074 if (output_stagger(idx) == M) then
2075 field%header%dim1(1) = we_mem_s
2076 field%header%dim1(2) = we_mem_e
2077 field%header%dim2(1) = sn_mem_s
2078 field%header%dim2(2) = sn_mem_e
2079 else if (output_stagger(idx) == U) then
2080 field%header%dim1(1) = we_mem_stag_s
2081 field%header%dim1(2) = we_mem_stag_e
2082 field%header%dim2(1) = sn_mem_s
2083 field%header%dim2(2) = sn_mem_e
2084 else if (output_stagger(idx) == V) then
2085 field%header%dim1(1) = we_mem_s
2086 field%header%dim1(2) = we_mem_e
2087 field%header%dim2(1) = sn_mem_stag_s
2088 field%header%dim2(2) = sn_mem_stag_e
2089 end if
2090 end if
2092 call fill_field(field, idx, output_flags)
2094 end if
2095 end do
2098 end subroutine create_derived_fields
2101 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2102 ! Name: fill_field
2103 !
2104 ! Purpose:
2105 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2106 subroutine fill_field(field, idx, output_flags, all_level_list)
2108 use interp_option_module
2109 use list_module
2110 use module_mergesort
2111 use storage_module
2113 implicit none
2115 ! Arguments
2116 integer, intent(in) :: idx
2117 type (fg_input), intent(inout) :: field
2118 character (len=128), dimension(:), intent(inout) :: output_flags
2119 real, dimension(:), intent(in), optional :: all_level_list
2121 ! Local variables
2122 integer :: i, j, istatus, isrclevel
2123 integer, pointer, dimension(:) :: all_list
2124 real :: rfillconst, rlevel, rsrclevel
2125 type (fg_input) :: query_field
2126 type (list_item), pointer, dimension(:) :: keys
2127 character (len=128) :: asrcname
2128 logical :: filled_all_lev
2130 filled_all_lev = .false.
2132 !
2133 ! Get a list of all levels to be filled for this field
2134 !
2135 keys => list_get_keys(fill_lev_list(idx))
2137 do i=1,list_length(fill_lev_list(idx))
2139 !
2140 ! First handle a specification for levels "all"
2141 !
2142 if (trim(keys(i)%ckey) == 'all') then
2144 ! We only want to fill all levels if we haven't already filled "all" of them
2145 if (.not. filled_all_lev) then
2147 filled_all_lev = .true.
2149 query_field%header%time_dependent = .true.
2150 query_field%header%field = ' '
2151 nullify(query_field%r_arr)
2152 nullify(query_field%valid_mask)
2153 nullify(query_field%modified_mask)
2155 ! See if we are filling this level with a constant
2156 call get_constant_fill_lev(keys(i)%cvalue, rfillconst, istatus)
2157 if (istatus == 0) then
2158 if (present(all_level_list)) then
2159 do j=1,size(all_level_list)
2160 call create_level(field, real(all_level_list(j)), idx, output_flags, rfillconst=rfillconst)
2161 end do
2162 else
2163 query_field%header%field = level_template(idx)
2164 nullify(all_list)
2165 call storage_get_levels(query_field, all_list)
2166 if (associated(all_list)) then
2167 do j=1,size(all_list)
2168 call create_level(field, real(all_list(j)), idx, output_flags, rfillconst=rfillconst)
2169 end do
2170 deallocate(all_list)
2171 end if
2172 end if
2174 ! Else see if we are filling this level with a constant equal
2175 ! to the value of the level
2176 else if (trim(keys(i)%cvalue) == 'vertical_index') then
2177 if (present(all_level_list)) then
2178 do j=1,size(all_level_list)
2179 call create_level(field, real(all_level_list(j)), idx, output_flags, &
2180 rfillconst=real(all_level_list(j)))
2181 end do
2182 else
2183 query_field%header%field = level_template(idx)
2184 nullify(all_list)
2185 call storage_get_levels(query_field, all_list)
2186 if (associated(all_list)) then
2187 do j=1,size(all_list)
2188 call create_level(field, real(all_list(j)), idx, output_flags, rfillconst=real(all_list(j)))
2189 end do
2190 deallocate(all_list)
2191 end if
2192 end if
2194 ! Else, we assume that it is a field from which we are copying levels
2195 else
2196 if (present(all_level_list)) then
2197 do j=1,size(all_level_list)
2198 call create_level(field, real(all_level_list(j)), idx, output_flags, &
2199 asrcname=keys(i)%cvalue, rsrclevel=real(all_level_list(j)))
2200 end do
2201 else
2202 query_field%header%field = keys(i)%cvalue ! Use same levels as source field, not level_template
2203 nullify(all_list)
2204 call storage_get_levels(query_field, all_list)
2205 if (associated(all_list)) then
2206 do j=1,size(all_list)
2207 call create_level(field, real(all_list(j)), idx, output_flags, &
2208 asrcname=keys(i)%cvalue, rsrclevel=real(all_list(j)))
2209 end do
2210 deallocate(all_list)
2212 else
2214 ! If the field doesn't have any levels (or does not exist) then we have not
2215 ! really filled all levels at this point.
2216 filled_all_lev = .false.
2217 end if
2218 end if
2220 end if
2221 end if
2223 !
2224 ! Handle individually specified levels
2225 !
2226 else
2228 read(keys(i)%ckey,*) rlevel
2230 ! See if we are filling this level with a constant
2231 call get_constant_fill_lev(keys(i)%cvalue, rfillconst, istatus)
2232 if (istatus == 0) then
2233 call create_level(field, rlevel, idx, output_flags, rfillconst=rfillconst)
2235 ! Otherwise, we are filling from another level
2236 else
2237 call get_fill_src_level(keys(i)%cvalue, asrcname, isrclevel)
2238 rsrclevel = real(isrclevel)
2239 call create_level(field, rlevel, idx, output_flags, &
2240 asrcname=asrcname, rsrclevel=rsrclevel)
2242 end if
2243 end if
2244 end do
2246 if (associated(keys)) deallocate(keys)
2248 end subroutine fill_field
2251 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2252 ! Name: create_level
2253 !
2254 ! Purpose:
2255 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2256 subroutine create_level(field_template, rlevel, idx, output_flags, &
2257 rfillconst, asrcname, rsrclevel)
2259 use storage_module
2260 use interp_option_module
2262 implicit none
2264 ! Arguments
2265 type (fg_input), intent(inout) :: field_template
2266 real, intent(in) :: rlevel
2267 integer, intent(in) :: idx
2268 character (len=128), dimension(:), intent(inout) :: output_flags
2269 real, intent(in), optional :: rfillconst, rsrclevel
2270 character (len=128), intent(in), optional :: asrcname
2272 ! Local variables
2273 integer :: i, j, istatus
2274 integer :: sm1,em1,sm2,em2
2275 type (fg_input) :: query_field
2277 !
2278 ! Check to make sure optional arguments are sane
2279 !
2280 if (present(rfillconst) .and. (present(asrcname) .or. present(rsrclevel))) then
2281 call mprintf(.true.,ERROR,'A call to create_level() cannot be given specifications '// &
2282 'for both a constant fill value and a source level.')
2284 else if ((present(asrcname) .and. .not. present(rsrclevel)) .or. &
2285 (.not. present(asrcname) .and. present(rsrclevel))) then
2286 call mprintf(.true.,ERROR,'Neither or both of optional arguments asrcname and '// &
2287 'rsrclevel must be specified to subroutine create_level().')
2289 else if (.not. present(rfillconst) .and. &
2290 .not. present(asrcname) .and. &
2291 .not. present(rsrclevel)) then
2292 call mprintf(.true.,ERROR,'A call to create_level() must be given either a specification '// &
2293 'for a constant fill value or a source level.')
2294 end if
2296 query_field%header%time_dependent = .true.
2297 query_field%header%field = field_template%header%field
2298 query_field%header%vertical_level = rlevel
2299 nullify(query_field%r_arr)
2300 nullify(query_field%valid_mask)
2301 nullify(query_field%modified_mask)
2303 call storage_query_field(query_field, istatus)
2304 if (istatus == 0) then
2305 call mprintf(.true.,INFORM,'%s at level %f already exists; leaving it alone.', &
2306 s1=field_template%header%field, f1=rlevel)
2307 return
2308 end if
2310 sm1 = field_template%header%dim1(1)
2311 em1 = field_template%header%dim1(2)
2312 sm2 = field_template%header%dim2(1)
2313 em2 = field_template%header%dim2(2)
2315 !
2316 ! Handle constant fill value case
2317 !
2318 if (present(rfillconst)) then
2320 field_template%header%vertical_level = rlevel
2321 allocate(field_template%r_arr(sm1:em1,sm2:em2))
2322 allocate(field_template%valid_mask)
2323 allocate(field_template%modified_mask)
2324 call bitarray_create(field_template%valid_mask, em1-sm1+1, em2-sm2+1)
2325 call bitarray_create(field_template%modified_mask, em1-sm1+1, em2-sm2+1)
2327 field_template%r_arr = rfillconst
2329 do j=sm2,em2
2330 do i=sm1,em1
2331 call bitarray_set(field_template%valid_mask, i-sm1+1, j-sm2+1)
2332 end do
2333 end do
2335 call storage_put_field(field_template)
2337 if (output_this_field(idx) .and. flag_in_output(idx) /= ' ') then
2338 output_flags(idx) = flag_in_output(idx)
2339 end if
2341 !
2342 ! Handle source field and source level case
2343 !
2344 else if (present(asrcname) .and. present(rsrclevel)) then
2346 query_field%header%field = ' '
2347 query_field%header%field = asrcname
2348 query_field%header%vertical_level = rsrclevel
2350 ! Check to see whether the requested source field exists at the requested level
2351 call storage_query_field(query_field, istatus)
2353 if (istatus == 0) then
2355 ! Read in requested field at requested level
2356 call storage_get_field(query_field, istatus)
2357 if ((query_field%header%dim1(1) /= field_template%header%dim1(1)) .or. &
2358 (query_field%header%dim1(2) /= field_template%header%dim1(2)) .or. &
2359 (query_field%header%dim2(1) /= field_template%header%dim2(1)) .or. &
2360 (query_field%header%dim2(2) /= field_template%header%dim2(2))) then
2361 call mprintf(.true.,ERROR,'Dimensions for %s do not match those of %s. This is '// &
2362 'probably because the staggerings of the fields do not match.', &
2363 s1=query_field%header%field, s2=field_template%header%field)
2364 end if
2366 field_template%header%vertical_level = rlevel
2367 allocate(field_template%r_arr(sm1:em1,sm2:em2))
2368 allocate(field_template%valid_mask)
2369 allocate(field_template%modified_mask)
2370 call bitarray_create(field_template%valid_mask, em1-sm1+1, em2-sm2+1)
2371 call bitarray_create(field_template%modified_mask, em1-sm1+1, em2-sm2+1)
2373 field_template%r_arr = query_field%r_arr
2375 ! We should retain information about which points in the field are valid
2376 do j=sm2,em2
2377 do i=sm1,em1
2378 if (bitarray_test(query_field%valid_mask, i-sm1+1, j-sm2+1)) then
2379 call bitarray_set(field_template%valid_mask, i-sm1+1, j-sm2+1)
2380 end if
2381 end do
2382 end do
2384 call storage_put_field(field_template)
2386 if (output_this_field(idx) .and. flag_in_output(idx) /= ' ') then
2387 output_flags(idx) = flag_in_output(idx)
2388 end if
2390 else
2391 call mprintf(.true.,INFORM,'Couldn''t find %s at level %f to fill level %f of %s.', &
2392 s1=asrcname,f1=rsrclevel,f2=rlevel,s2=field_template%header%field)
2393 end if
2395 end if
2397 end subroutine create_level
2400 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2401 ! Name: accum_continuous
2402 !
2403 ! Purpose: Sum up all of the source data points whose nearest neighbor in the
2404 ! model grid is the specified model grid point.
2405 !
2406 ! NOTE: When processing the source tile, those source points that are
2407 ! closest to a different model grid point will be added to the totals for
2408 ! such grid points; thus, an entire source tile will be processed at a time.
2409 ! This routine really processes for all model grid points that are
2410 ! within a source tile, and not just for a single grid point.
2411 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2412 subroutine accum_continuous(src_array, &
2413 src_min_x, src_max_x, src_min_y, src_max_y, src_min_z, src_max_z, bdr_width, &
2414 dst_array, n, &
2415 start_i, end_i, start_j, end_j, start_k, end_k, &
2416 istagger, &
2417 new_pts, msgval, maskval, mask_array, sr_x, sr_y)
2419 use bitarray_module
2420 use misc_definitions_module
2422 implicit none
2424 ! Arguments
2425 integer, intent(in) :: start_i, end_i, start_j, end_j, start_k, end_k, istagger, &
2426 src_min_x, src_max_x, src_min_y, src_max_y, src_min_z, src_max_z, bdr_width
2427 real, intent(in) :: maskval, msgval
2428 real, dimension(src_min_x:src_max_x, src_min_y:src_max_y, src_min_z:src_max_z), intent(in) :: src_array
2429 real, dimension(start_i:end_i, start_j:end_j, start_k:end_k), intent(inout) :: dst_array, n
2430 real, dimension(src_min_x:src_max_x, src_min_y:src_max_y), intent(in), optional :: mask_array
2431 integer, intent(in), optional :: sr_x, sr_y
2432 type (bitarray), intent(inout) :: new_pts
2434 ! Local variables
2435 integer :: i, j
2436 integer, pointer, dimension(:,:,:) :: where_maps_to
2437 real :: rsr_x, rsr_y
2439 rsr_x = 1.0
2440 rsr_y = 1.0
2441 if (present(sr_x)) rsr_x = real(sr_x)
2442 if (present(sr_y)) rsr_y = real(sr_y)
2444 allocate(where_maps_to(src_min_x:src_max_x,src_min_y:src_max_y,2))
2445 do i=src_min_x,src_max_x
2446 do j=src_min_y,src_max_y
2447 where_maps_to(i,j,1) = NOT_PROCESSED
2448 end do
2449 end do
2451 call process_continuous_block(src_array, where_maps_to, &
2452 src_min_x, src_min_y, src_min_z, src_max_x, src_max_y, src_max_z, &
2453 src_min_x+bdr_width, src_min_y, src_min_z, &
2454 src_max_x-bdr_width, src_max_y, src_max_z, &
2455 dst_array, n, start_i, end_i, start_j, end_j, start_k, end_k, &
2456 istagger, &
2457 new_pts, rsr_x, rsr_y, msgval, maskval, mask_array)
2459 deallocate(where_maps_to)
2461 end subroutine accum_continuous
2464 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2465 ! Name: process_continuous_block
2466 !
2467 ! Purpose: To recursively process a subarray of continuous data, adding the
2468 ! points in a block to the sum for their nearest grid point. The nearest
2469 ! neighbor may be estimated in some cases; for example, if the four corners
2470 ! of a subarray all have the same nearest grid point, all elements in the
2471 ! subarray are added to that grid point.
2472 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2473 recursive subroutine process_continuous_block(tile_array, where_maps_to, &
2474 src_min_x, src_min_y, src_min_z, src_max_x, src_max_y, src_max_z, &
2475 min_i, min_j, min_k, max_i, max_j, max_k, &
2476 dst_array, n, &
2477 start_x, end_x, start_y, end_y, start_z, end_z, &
2478 istagger, &
2479 new_pts, sr_x, sr_y, msgval, maskval, mask_array)
2481 use bitarray_module
2482 use llxy_module
2483 use misc_definitions_module
2485 implicit none
2487 ! Arguments
2488 integer, intent(in) :: min_i, min_j, min_k, max_i, max_j, max_k, &
2489 src_min_x, src_min_y, src_min_z, src_max_x, src_max_y, src_max_z, &
2490 start_x, end_x, start_y, end_y, start_z, end_z, istagger
2491 integer, dimension(src_min_x:src_max_x,src_min_y:src_max_y,2), intent(inout) :: where_maps_to
2492 real, intent(in) :: sr_x, sr_y, maskval, msgval
2493 real, dimension(src_min_x:src_max_x,src_min_y:src_max_y,src_min_z:src_max_z), intent(in) :: tile_array
2494 real, dimension(src_min_x:src_max_x,src_min_y:src_max_y), intent(in), optional :: mask_array
2495 real, dimension(start_x:end_x,start_y:end_y,start_z:end_z), intent(inout) :: dst_array, n
2496 type (bitarray), intent(inout) :: new_pts
2498 ! Local variables
2499 integer :: orig_selected_domain, x_dest, y_dest, i, j, k, center_i, center_j
2500 real :: lat_corner, lon_corner, rx, ry
2502 ! Compute the model grid point that the corners of the rectangle to be
2503 ! processed map to
2504 ! Lower-left corner
2505 if (where_maps_to(min_i,min_j,1) == NOT_PROCESSED) then
2506 orig_selected_domain = iget_selected_domain()
2507 call select_domain(SOURCE_PROJ)
2508 call xytoll(real(min_i), real(min_j), lat_corner, lon_corner, istagger)
2509 call select_domain(1)
2510 call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
2511 rx = (rx - 1.0)*sr_x + 1.0
2512 ry = (ry - 1.0)*sr_y + 1.0
2513 call select_domain(orig_selected_domain)
2514 if (real(start_x) <= rx .and. rx <= real(end_x) .and. &
2515 real(start_y) <= ry .and. ry <= real(end_y)) then
2516 where_maps_to(min_i,min_j,1) = nint(rx)
2517 where_maps_to(min_i,min_j,2) = nint(ry)
2518 else
2519 where_maps_to(min_i,min_j,1) = OUTSIDE_DOMAIN
2520 end if
2521 end if
2523 ! Upper-left corner
2524 if (where_maps_to(min_i,max_j,1) == NOT_PROCESSED) then
2525 orig_selected_domain = iget_selected_domain()
2526 call select_domain(SOURCE_PROJ)
2527 call xytoll(real(min_i), real(max_j), lat_corner, lon_corner, istagger)
2528 call select_domain(1)
2529 call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
2530 rx = (rx - 1.0)*sr_x + 1.0
2531 ry = (ry - 1.0)*sr_y + 1.0
2532 call select_domain(orig_selected_domain)
2533 if (real(start_x) <= rx .and. rx <= real(end_x) .and. &
2534 real(start_y) <= ry .and. ry <= real(end_y)) then
2535 where_maps_to(min_i,max_j,1) = nint(rx)
2536 where_maps_to(min_i,max_j,2) = nint(ry)
2537 else
2538 where_maps_to(min_i,max_j,1) = OUTSIDE_DOMAIN
2539 end if
2540 end if
2542 ! Upper-right corner
2543 if (where_maps_to(max_i,max_j,1) == NOT_PROCESSED) then
2544 orig_selected_domain = iget_selected_domain()
2545 call select_domain(SOURCE_PROJ)
2546 call xytoll(real(max_i), real(max_j), lat_corner, lon_corner, istagger)
2547 call select_domain(1)
2548 call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
2549 rx = (rx - 1.0)*sr_x + 1.0
2550 ry = (ry - 1.0)*sr_y + 1.0
2551 call select_domain(orig_selected_domain)
2552 if (real(start_x) <= rx .and. rx <= real(end_x) .and. &
2553 real(start_y) <= ry .and. ry <= real(end_y)) then
2554 where_maps_to(max_i,max_j,1) = nint(rx)
2555 where_maps_to(max_i,max_j,2) = nint(ry)
2556 else
2557 where_maps_to(max_i,max_j,1) = OUTSIDE_DOMAIN
2558 end if
2559 end if
2561 ! Lower-right corner
2562 if (where_maps_to(max_i,min_j,1) == NOT_PROCESSED) then
2563 orig_selected_domain = iget_selected_domain()
2564 call select_domain(SOURCE_PROJ)
2565 call xytoll(real(max_i), real(min_j), lat_corner, lon_corner, istagger)
2566 call select_domain(1)
2567 call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
2568 rx = (rx - 1.0)*sr_x + 1.0
2569 ry = (ry - 1.0)*sr_y + 1.0
2570 call select_domain(orig_selected_domain)
2571 if (real(start_x) <= rx .and. rx <= real(end_x) .and. &
2572 real(start_y) <= ry .and. ry <= real(end_y)) then
2573 where_maps_to(max_i,min_j,1) = nint(rx)
2574 where_maps_to(max_i,min_j,2) = nint(ry)
2575 else
2576 where_maps_to(max_i,min_j,1) = OUTSIDE_DOMAIN
2577 end if
2578 end if
2580 ! If all four corners map to same model grid point, accumulate the
2581 ! entire rectangle
2582 if (where_maps_to(min_i,min_j,1) == where_maps_to(min_i,max_j,1) .and. &
2583 where_maps_to(min_i,min_j,1) == where_maps_to(max_i,max_j,1) .and. &
2584 where_maps_to(min_i,min_j,1) == where_maps_to(max_i,min_j,1) .and. &
2585 where_maps_to(min_i,min_j,2) == where_maps_to(min_i,max_j,2) .and. &
2586 where_maps_to(min_i,min_j,2) == where_maps_to(max_i,max_j,2) .and. &
2587 where_maps_to(min_i,min_j,2) == where_maps_to(max_i,min_j,2) .and. &
2588 where_maps_to(min_i,min_j,1) /= OUTSIDE_DOMAIN) then
2589 x_dest = where_maps_to(min_i,min_j,1)
2590 y_dest = where_maps_to(min_i,min_j,2)
2592 ! If this grid point was already given a value from higher-priority source data,
2593 ! there is nothing to do.
2594 ! if (.not. bitarray_test(processed_pts, x_dest-start_x+1, y_dest-start_y+1)) then
2596 ! If this grid point has never been given a value by this level of source data,
2597 ! initialize the point
2598 if (.not. bitarray_test(new_pts, x_dest-start_x+1, y_dest-start_y+1)) then
2599 do k=min_k,max_k
2600 dst_array(x_dest,y_dest,k) = 0.
2601 end do
2602 end if
2604 ! Sum all the points whose nearest neighbor is this grid point
2605 if (present(mask_array)) then
2606 do i=min_i,max_i
2607 do j=min_j,max_j
2608 do k=min_k,max_k
2609 ! Ignore masked/missing values in the source data
2610 if ((tile_array(i,j,k) /= msgval) .and. &
2611 (mask_array(i,j) /= maskval)) then
2612 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2613 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2614 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2615 end if
2616 end do
2617 end do
2618 end do
2619 else
2620 do i=min_i,max_i
2621 do j=min_j,max_j
2622 do k=min_k,max_k
2623 ! Ignore masked/missing values in the source data
2624 if ((tile_array(i,j,k) /= msgval)) then
2625 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2626 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2627 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2628 end if
2629 end do
2630 end do
2631 end do
2632 end if
2634 ! end if
2636 ! Rectangle is a square of four points, and we can simply deal with each of the points
2637 else if (((max_i - min_i + 1) <= 2) .and. ((max_j - min_j + 1) <= 2)) then
2638 do i=min_i,max_i
2639 do j=min_j,max_j
2640 x_dest = where_maps_to(i,j,1)
2641 y_dest = where_maps_to(i,j,2)
2643 if (x_dest /= OUTSIDE_DOMAIN) then
2645 ! if (.not. bitarray_test(processed_pts, x_dest-start_x+1, y_dest-start_y+1)) then
2646 if (.not. bitarray_test(new_pts, x_dest-start_x+1, y_dest-start_y+1)) then
2647 do k=min_k,max_k
2648 dst_array(x_dest,y_dest,k) = 0.
2649 end do
2650 end if
2652 if (present(mask_array)) then
2653 do k=min_k,max_k
2654 ! Ignore masked/missing values
2655 if ((tile_array(i,j,k) /= msgval) .and. &
2656 (mask_array(i,j) /= maskval)) then
2657 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2658 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2659 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2660 end if
2661 end do
2662 else
2663 do k=min_k,max_k
2664 ! Ignore masked/missing values
2665 if ((tile_array(i,j,k) /= msgval)) then
2666 dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
2667 n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
2668 call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
2669 end if
2670 end do
2671 end if
2672 ! end if
2674 end if
2675 end do
2676 end do
2678 ! Not all corners map to the same grid point, and the rectangle contains more than
2679 ! four points
2680 else
2681 center_i = (max_i + min_i)/2
2682 center_j = (max_j + min_j)/2
2684 ! Recursively process lower-left rectangle
2685 call process_continuous_block(tile_array, where_maps_to, &
2686 src_min_x, src_min_y, src_min_z, &
2687 src_max_x, src_max_y, src_max_z, &
2688 min_i, min_j, min_k, &
2689 center_i, center_j, max_k, &
2690 dst_array, n, &
2691 start_x, end_x, start_y, end_y, start_z, end_z, &
2692 istagger, &
2693 new_pts, sr_x, sr_y, msgval, maskval, mask_array)
2695 if (center_i < max_i) then
2696 ! Recursively process lower-right rectangle
2697 call process_continuous_block(tile_array, where_maps_to, &
2698 src_min_x, src_min_y, src_min_z, &
2699 src_max_x, src_max_y, src_max_z, &
2700 center_i+1, min_j, min_k, max_i, &
2701 center_j, max_k, &
2702 dst_array, n, &
2703 start_x, end_x, start_y, &
2704 end_y, start_z, end_z, &
2705 istagger, &
2706 new_pts, sr_x, sr_y, msgval, maskval, mask_array)
2707 end if
2709 if (center_j < max_j) then
2710 ! Recursively process upper-left rectangle
2711 call process_continuous_block(tile_array, where_maps_to, &
2712 src_min_x, src_min_y, src_min_z, &
2713 src_max_x, src_max_y, src_max_z, &
2714 min_i, center_j+1, min_k, center_i, &
2715 max_j, max_k, &
2716 dst_array, n, &
2717 start_x, end_x, start_y, &
2718 end_y, start_z, end_z, &
2719 istagger, &
2720 new_pts, sr_x, sr_y, msgval, maskval, mask_array)
2721 end if
2723 if (center_i < max_i .and. center_j < max_j) then
2724 ! Recursively process upper-right rectangle
2725 call process_continuous_block(tile_array, where_maps_to, &
2726 src_min_x, src_min_y, src_min_z, &
2727 src_max_x, src_max_y, src_max_z, &
2728 center_i+1, center_j+1, min_k, max_i, &
2729 max_j, max_k, &
2730 dst_array, n, &
2731 start_x, end_x, start_y, &
2732 end_y, start_z, end_z, &
2733 istagger, &
2734 new_pts, sr_x, sr_y, msgval, maskval, mask_array)
2735 end if
2736 end if
2738 end subroutine process_continuous_block
2740 end module process_domain_module