module process_domain_module
contains
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Name: process_domain
!
! Purpose:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine process_domain(n, extra_row, extra_col)
use date_pack
use gridinfo_module
use interp_option_module
use misc_definitions_module
use module_debug
use storage_module
implicit none
! Arguments
integer, intent(in) :: n
logical, intent(in) :: extra_row, extra_col
! Local variables
integer :: i, t, dyn_opt, &
we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
idiff, n_times, &
west_east_dim, south_north_dim, bottom_top_dim, map_proj, &
is_water, is_ice, is_urban, i_soilwater, &
grid_id, parent_id, i_parent_start, j_parent_start, &
i_parent_end, j_parent_end, parent_grid_ratio
real :: cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
dom_dx, dom_dy
real, dimension(16) :: corner_lats, corner_lons
real, pointer, dimension(:,:) :: landmask
real, pointer, dimension(:,:) :: xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v
logical, allocatable, dimension(:) :: got_this_field, got_const_field
character (len=19) :: valid_date, temp_date
character (len=128) :: title, mminlu
character (len=128), allocatable, dimension(:) :: output_flags, td_output_flags
! Compute number of times that we will process
call mprintf(.true.,LOGFILE,'Into subroutine')
call geth_idts(end_date(n), start_date(n), idiff)
call mprintf((idiff < 0),ERROR,'Ending date is earlier than starting date in namelist for domain %i.', i1=n)
! only do the initial time if initializing a nest
!mp if (n .eq. 1) then
n_times = idiff / interval_seconds
!mp else
!mp n_times = 0
!mp endif
! Check that the interval evenly divides the range of times to process
call mprintf((mod(idiff, interval_seconds) /= 0),WARN, &
'In namelist, interval_seconds does not evenly divide '// &
'(end_date - start_date) for domain %i. Only %i time periods '// &
'will be processed.', i1=n, i2=n_times)
! Initialize the storage module
call mprintf(.true.,LOGFILE,'Initializing storage module')
call storage_init()
!
! Do time-independent processing
!
call mprintf(.true.,LOGFILE,'To get_static_fields')
call get_static_fields(n, dyn_opt, west_east_dim, south_north_dim, bottom_top_dim, map_proj, &
we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
we_patch_s, we_patch_e, &
we_patch_stag_s, we_patch_stag_e, &
sn_patch_s, sn_patch_e, &
sn_patch_stag_s, sn_patch_stag_e, &
we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
mminlu, is_water, is_ice, is_urban, i_soilwater, &
grid_id, parent_id, &
i_parent_start, j_parent_start, i_parent_end, j_parent_end, &
parent_grid_ratio, cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
dom_dx, dom_dy, landmask, xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v, corner_lats, &
corner_lons, title)
call mprintf(.true.,LOGFILE,'Past get_static_fields')
allocate(output_flags(num_entries))
allocate(got_const_field(num_entries))
do i=1,num_entries
output_flags(i) = ' '
got_const_field(i) = .false.
end do
! This call is to process the constant met fields (SST or SEAICE, for example)
! That we process constant fields is indicated by the first argument
call process_single_met_time(.true., temp_date, n, extra_row, extra_col, xlat, xlon, &
xlat_u, xlon_u, xlat_v, xlon_v, landmask, &
title, dyn_opt, &
west_east_dim, south_north_dim, &
we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
got_const_field, &
map_proj, mminlu, is_water, is_ice, is_urban, i_soilwater, &
grid_id, parent_id, i_parent_start, &
j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, &
truelat2, parent_grid_ratio, corner_lats, corner_lons, output_flags)
!
! Begin time-dependent processing
!
allocate(td_output_flags(num_entries))
allocate(got_this_field (num_entries))
! Loop over all times to be processed for this domain
do t=0,n_times
call geth_newdate(valid_date, trim(start_date(n)), t*interval_seconds)
temp_date = ' '
if (mod(interval_seconds,3600) == 0) then
write(temp_date,'(a13)') valid_date(1:10)//'_'//valid_date(12:13)
else if (mod(interval_seconds,60) == 0) then
write(temp_date,'(a16)') valid_date(1:10)//'_'//valid_date(12:16)
else
write(temp_date,'(a19)') valid_date(1:10)//'_'//valid_date(12:19)
end if
call mprintf(.true.,STDOUT, ' Processing %s', s1=trim(temp_date))
call mprintf(.true.,LOGFILE, 'Preparing to process output time %s', s1=temp_date)
do i=1,num_entries
td_output_flags(i) = output_flags(i)
got_this_field(i) = got_const_field(i)
end do
call process_single_met_time(.false., temp_date, n, extra_row, extra_col, xlat, xlon, &
xlat_u, xlon_u, xlat_v, xlon_v, landmask, &
title, dyn_opt, &
west_east_dim, south_north_dim, &
we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
got_this_field, &
map_proj, mminlu, is_water, is_ice, is_urban, i_soilwater, &
grid_id, parent_id, i_parent_start, &
j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, &
truelat2, parent_grid_ratio, corner_lats, corner_lons, td_output_flags)
end do ! Loop over n_times
deallocate(td_output_flags)
deallocate(got_this_field)
deallocate(output_flags)
deallocate(got_const_field)
call storage_delete_all()
end subroutine process_domain
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Name: get_static_fields
!
! Purpose:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine get_static_fields(n, dyn_opt, west_east_dim, south_north_dim, bottom_top_dim, &
map_proj, &
we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
mminlu, is_water, is_ice, is_urban, i_soilwater, grid_id, parent_id, &
i_parent_start, j_parent_start, i_parent_end, j_parent_end, &
parent_grid_ratio, cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
dom_dx, dom_dy, landmask, xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v, corner_lats, &
corner_lons, title)
use gridinfo_module
use input_module
use llxy_module
use parallel_module
use storage_module
implicit none
! Arguments
integer, intent(in) :: n
integer, intent(inout) :: dyn_opt, west_east_dim, south_north_dim, bottom_top_dim, &
map_proj, &
we_dom_s, we_dom_e, sn_dom_s, sn_dom_e, &
we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
is_water, is_ice, is_urban, i_soilwater, grid_id, parent_id, &
i_parent_start, j_parent_start, i_parent_end, j_parent_end, &
parent_grid_ratio
real, pointer, dimension(:,:) :: landmask
real, intent(inout) :: cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
dom_dx, dom_dy
real, pointer, dimension(:,:) :: xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v
real, dimension(16), intent(out) :: corner_lats, corner_lons
character (len=128), intent(inout) :: title, mminlu
! Local variables
integer :: istatus, i, j, k, sp1, ep1, sp2, ep2, sp3, ep3, &
lh_mult, rh_mult, bh_mult, th_mult
real, pointer, dimension(:,:,:) :: real_array
character (len=3) :: memorder
character (len=128) :: grid_type, datestr, cname, stagger, cunits, cdesc
character (len=128), dimension(3) :: dimnames
type (fg_input) :: field
! Initialize the input module to read static input data for this domain
call mprintf(.true.,LOGFILE,'Opening static input file.')
call input_init(n, istatus)
call mprintf((istatus /= 0),ERROR, 'input_init(): Error opening input for domain %i.', i1=n)
! Read global attributes from the static data input file
call mprintf(.true.,LOGFILE,'Reading static global attributes.')
call read_global_attrs(title, datestr, grid_type, dyn_opt, west_east_dim, &
south_north_dim, bottom_top_dim, &
we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
map_proj, mminlu, is_water, is_ice, is_urban, i_soilwater, &
grid_id, parent_id, i_parent_start, &
j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, &
truelat2, parent_grid_ratio, corner_lats, corner_lons)
we_dom_s = 1
sn_dom_s = 1
if (grid_type(1:1) == 'C') then
we_dom_e = west_east_dim - 1
sn_dom_e = south_north_dim - 1
else if (grid_type(1:1) == 'E') then
we_dom_e = west_east_dim
sn_dom_e = south_north_dim
else if (grid_type(1:1) == 'B' .or. grid_type(1:1) == 'A') then
we_dom_e = west_east_dim
sn_dom_e = south_north_dim
end if
! Given the full dimensions of this domain, find out the range of indices
! that will be worked on by this processor. This information is given by
! my_minx, my_miny, my_maxx, my_maxy
call parallel_get_tile_dims(west_east_dim, south_north_dim)
! Must figure out patch dimensions from info in parallel module
if (nprocs > 1 .and. .not. do_tiled_input) then
we_patch_s = my_minx
we_patch_stag_s = my_minx
we_patch_e = my_maxx - 1
sn_patch_s = my_miny
sn_patch_stag_s = my_miny
sn_patch_e = my_maxy - 1
if (gridtype == 'C') then
if (my_x /= nproc_x - 1) then
we_patch_e = we_patch_e + 1
we_patch_stag_e = we_patch_e
else
we_patch_stag_e = we_patch_e + 1
end if
if (my_y /= nproc_y - 1) then
sn_patch_e = sn_patch_e + 1
sn_patch_stag_e = sn_patch_e
else
sn_patch_stag_e = sn_patch_e + 1
end if
else if (gridtype == 'E' .or. gridtype == 'B' .or. gridtype == 'A') then
we_patch_e = we_patch_e + 1
sn_patch_e = sn_patch_e + 1
we_patch_stag_e = we_patch_e
sn_patch_stag_e = sn_patch_e
end if
end if
call mprintf(.true.,LOGFILE,'my_minx is %i and my_maxx is %i.', i1=my_minx, i2=my_maxx)
call mprintf(.true.,LOGFILE,'my_miny is %i and my_maxy is %i.', i1=my_miny, i2=my_maxy)
call mprintf(.true.,LOGFILE,'we_patch_e is %i and we_patch_stag_e is %i.', i1=we_patch_e, i2=we_patch_stag_e)
! Compute multipliers for halo width; these must be 0/1
if (my_x /= 0) then
lh_mult = 1
else
lh_mult = 0
end if
if (my_x /= (nproc_x-1)) then
rh_mult = 1
else
rh_mult = 0
end if
if (my_y /= 0) then
bh_mult = 1
else
bh_mult = 0
end if
if (my_y /= (nproc_y-1)) then
th_mult = 1
else
th_mult = 0
end if
we_mem_s = we_patch_s - HALO_WIDTH*lh_mult
we_mem_e = we_patch_e + HALO_WIDTH*rh_mult
sn_mem_s = sn_patch_s - HALO_WIDTH*bh_mult
sn_mem_e = sn_patch_e + HALO_WIDTH*th_mult
we_mem_stag_s = we_patch_stag_s - HALO_WIDTH*lh_mult
we_mem_stag_e = we_patch_stag_e + HALO_WIDTH*rh_mult
sn_mem_stag_s = sn_patch_stag_s - HALO_WIDTH*bh_mult
sn_mem_stag_e = sn_patch_stag_e + HALO_WIDTH*th_mult
! Initialize a proj_info type for the destination grid projection. This will
! primarily be used for rotating Earth-relative winds to grid-relative winds
call mprintf(.true.,INFORM,' - Call set_domain_projection.')
call set_domain_projection(map_proj, stand_lon, truelat1, truelat2, &
dom_dx, dom_dy, dom_dx, dom_dy, west_east_dim, &
south_north_dim, real(west_east_dim)/2., &
real(south_north_dim)/2.,cen_lat, cen_lon, &
cen_lat, cen_lon)
! Read static fields using the input module; we know that there are no more
! fields to be read when read_next_field() returns a non-zero status.
istatus = 0
do while (istatus == 0)
call read_next_field(sp1, ep1, sp2, ep2, sp3, ep3, cname, cunits, cdesc, &
memorder, stagger, dimnames, real_array, istatus)
if (istatus == 0) then
call mprintf(.true.,LOGFILE, 'Read in static field %s.',s1=cname)
! We will also keep copies in core of the lat/lon arrays, for use in
! interpolation of the met fields to the model grid.
! For now, we assume that the lat/lon arrays will have known field names
if (index(cname, 'XLAT_M') /= 0 .and. &
len_trim(cname) == len_trim('XLAT_M')) then
allocate(xlat(we_mem_s:we_mem_e,sn_mem_s:sn_mem_e))
xlat(we_patch_s:we_patch_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
call exchange_halo_r(xlat, &
we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, 1, 1, &
we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, 1, 1)
else if (index(cname, 'XLONG_M') /= 0 .and. &
len_trim(cname) == len_trim('XLONG_M')) then
allocate(xlon(we_mem_s:we_mem_e,sn_mem_s:sn_mem_e))
xlon(we_patch_s:we_patch_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
call exchange_halo_r(xlon, &
we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, 1, 1, &
we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, 1, 1)
else if (index(cname, 'XLAT_U') /= 0 .and. &
len_trim(cname) == len_trim('XLAT_U')) then
allocate(xlat_u(we_mem_stag_s:we_mem_stag_e,sn_mem_s:sn_mem_e))
xlat_u(we_patch_stag_s:we_patch_stag_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
call exchange_halo_r(xlat_u, &
we_mem_stag_s, we_mem_stag_e, sn_mem_s, sn_mem_e, 1, 1, &
we_patch_stag_s, we_patch_stag_e, sn_patch_s, sn_patch_e, 1, 1)
else if (index(cname, 'XLONG_U') /= 0 .and. &
len_trim(cname) == len_trim('XLONG_U')) then
allocate(xlon_u(we_mem_stag_s:we_mem_stag_e,sn_mem_s:sn_mem_e))
xlon_u(we_patch_stag_s:we_patch_stag_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
call exchange_halo_r(xlon_u, &
we_mem_stag_s, we_mem_stag_e, sn_mem_s, sn_mem_e, 1, 1, &
we_patch_stag_s, we_patch_stag_e, sn_patch_s, sn_patch_e, 1, 1)
else if (index(cname, 'XLAT_V') /= 0 .and. &
len_trim(cname) == len_trim('XLAT_V')) then
allocate(xlat_v(we_mem_s:we_mem_e,sn_mem_stag_s:sn_mem_stag_e))
xlat_v(we_patch_s:we_patch_e,sn_patch_stag_s:sn_patch_stag_e) = real_array(:,:,1)
call exchange_halo_r(xlat_v, &
we_mem_s, we_mem_e, sn_mem_stag_s, sn_mem_stag_e, 1, 1, &
we_patch_s, we_patch_e, sn_patch_stag_s, sn_patch_stag_e, 1, 1)
else if (index(cname, 'XLONG_V') /= 0 .and. &
len_trim(cname) == len_trim('XLONG_V')) then
allocate(xlon_v(we_mem_s:we_mem_e,sn_mem_stag_s:sn_mem_stag_e))
xlon_v(we_patch_s:we_patch_e,sn_patch_stag_s:sn_patch_stag_e) = real_array(:,:,1)
call exchange_halo_r(xlon_v, &
we_mem_s, we_mem_e, sn_mem_stag_s, sn_mem_stag_e, 1, 1, &
we_patch_s, we_patch_e, sn_patch_stag_s, sn_patch_stag_e, 1, 1)
else if (index(cname, 'LANDMASK') /= 0 .and. &
len_trim(cname) == len_trim('LANDMASK')) then
allocate(landmask(we_mem_s:we_mem_e,sn_mem_s:sn_mem_e))
landmask(we_patch_s:we_patch_e,sn_patch_s:sn_patch_e) = real_array(:,:,1)
call exchange_halo_r(landmask, &
we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, 1, 1, &
we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, 1, 1)
end if
! Having read in a field, we write each level individually to the
! storage module; levels will be reassembled later on when they
! are written.
do k=sp3,ep3
field%header%version = 1
field%header%date = start_date(n)
field%header%time_dependent = .false.
field%header%mask_field = .false.
field%header%forecast_hour = 0.0
field%header%fg_source = 'geogrid_model'
field%header%field = cname
field%header%units = cunits
field%header%description = cdesc
field%header%vertical_coord = dimnames(3)
field%header%vertical_level = k
field%header%array_order = memorder
field%header%is_wind_grid_rel = .true.
field%header%array_has_missing_values = .false.
if (gridtype == 'C') then
if (trim(stagger) == 'M') then
field%map%stagger = M
field%header%dim1(1) = we_mem_s
field%header%dim1(2) = we_mem_e
field%header%dim2(1) = sn_mem_s
field%header%dim2(2) = sn_mem_e
else if (trim(stagger) == 'U') then
field%map%stagger = U
field%header%dim1(1) = we_mem_stag_s
field%header%dim1(2) = we_mem_stag_e
field%header%dim2(1) = sn_mem_s
field%header%dim2(2) = sn_mem_e
else if (trim(stagger) == 'V') then
field%map%stagger = V
field%header%dim1(1) = we_mem_s
field%header%dim1(2) = we_mem_e
field%header%dim2(1) = sn_mem_stag_s
field%header%dim2(2) = sn_mem_stag_e
end if
else if (gridtype == 'E' .or. gridtype == 'B') then
if (trim(stagger) == 'M') then
field%map%stagger = HH
else if (trim(stagger) == 'V') then
field%map%stagger = VV
end if
field%header%dim1(1) = we_mem_s
field%header%dim1(2) = we_mem_e
field%header%dim2(1) = sn_mem_s
field%header%dim2(2) = sn_mem_e
else if (gridtype == 'A') then
if (trim(stagger) == 'M') then
field%map%stagger = M
else if (trim(stagger) == 'V') then
print*, 'shouldnt have a V stagger for A grid'
field%map%stagger = M
end if
field%header%dim1(1) = we_mem_s
field%header%dim1(2) = we_mem_e
field%header%dim2(1) = sn_mem_s
field%header%dim2(2) = sn_mem_e
end if
allocate(field%valid_mask)
if (field%map%stagger == M .or. &
field%map%stagger == HH .or. &
field%map%stagger == VV) then
allocate(field%r_arr(we_mem_s:we_mem_e,&
sn_mem_s:sn_mem_e))
field%r_arr(we_patch_s:we_patch_e,sn_patch_s:sn_patch_e) = real_array(sp1:ep1,sp2:ep2,k)
call exchange_halo_r(field%r_arr, &
we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, 1, 1, &
we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, 1, 1)
call bitarray_create(field%valid_mask, &
(we_mem_e-we_mem_s)+1, &
(sn_mem_e-sn_mem_s)+1)
do j=1,(sn_mem_e-sn_mem_s)+1
do i=1,(we_mem_e-we_mem_s)+1
call bitarray_set(field%valid_mask, i, j)
end do
end do
else if (field%map%stagger == U) then
allocate(field%r_arr(we_mem_stag_s:we_mem_stag_e,&
sn_mem_s:sn_mem_e))
field%r_arr(we_patch_stag_s:we_patch_stag_e,sn_patch_s:sn_patch_e) = real_array(sp1:ep1,sp2:ep2,k)
call exchange_halo_r(field%r_arr, &
we_mem_stag_s, we_mem_stag_e, sn_mem_s, sn_mem_e, 1, 1, &
we_patch_stag_s, we_patch_stag_e, sn_patch_s, sn_patch_e, 1, 1)
call bitarray_create(field%valid_mask, &
(we_mem_stag_e-we_mem_stag_s)+1, &
(sn_mem_e-sn_mem_s)+1)
do j=1,(sn_mem_e-sn_mem_s)+1
do i=1,(we_mem_stag_e-we_mem_stag_s)+1
call bitarray_set(field%valid_mask, i, j)
end do
end do
else if (field%map%stagger == V) then
allocate(field%r_arr(we_mem_s:we_mem_e,&
sn_mem_stag_s:sn_mem_stag_e))
field%r_arr(we_patch_s:we_patch_e,sn_patch_stag_s:sn_patch_stag_e) = real_array(sp1:ep1,sp2:ep2,k)
call exchange_halo_r(field%r_arr, &
we_mem_s, we_mem_e, sn_mem_stag_s, sn_mem_stag_e, 1, 1, &
we_patch_s, we_patch_e, sn_patch_stag_s, sn_patch_stag_e, 1, 1)
call bitarray_create(field%valid_mask, &
(we_mem_e-we_mem_s)+1, &
(sn_mem_stag_e-sn_mem_stag_s)+1)
do j=1,(sn_mem_stag_e-sn_mem_stag_s)+1
do i=1,(we_mem_e-we_mem_s)+1
call bitarray_set(field%valid_mask, i, j)
end do
end do
end if
nullify(field%modified_mask)
call storage_put_field(field)
end do
end if
end do
! Done reading all static fields for this domain
call input_close()
end subroutine get_static_fields
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Name: process_single_met_time
!
! Purpose:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine process_single_met_time(do_const_processing, &
temp_date, n, extra_row, extra_col, xlat, xlon, &
xlat_u, xlon_u, xlat_v, xlon_v, landmask, &
title, dyn_opt, &
west_east_dim, south_north_dim, &
we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
got_this_field, &
map_proj, mminlu, is_water, is_ice, is_urban, i_soilwater, &
grid_id, parent_id, i_parent_start, &
j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, &
truelat2, parent_grid_ratio, corner_lats, corner_lons, output_flags)
use bitarray_module
use gridinfo_module
use interp_module
use interp_option_module
use llxy_module
use misc_definitions_module
use module_debug
use output_module
use parallel_module
use read_met_module
use rotate_winds_module
use storage_module
implicit none
! Arguments
logical, intent(in) :: do_const_processing
integer, intent(in) :: n, dyn_opt, west_east_dim, south_north_dim, map_proj, &
we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
we_mem_s, we_mem_e, we_mem_stag_s, we_mem_stag_e, &
sn_mem_s, sn_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
is_water, is_ice, is_urban, i_soilwater, &
grid_id, parent_id, i_parent_start, j_parent_start, &
i_parent_end, j_parent_end, parent_grid_ratio
! BUG: Should we be passing these around as pointers, or just declare them as arrays?
real, pointer, dimension(:,:) :: landmask
real, intent(in) :: dom_dx, dom_dy, cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2
real, dimension(16), intent(in) :: corner_lats, corner_lons
real, pointer, dimension(:,:) :: xlat, xlon, xlat_u, xlon_u, xlat_v, xlon_v
logical, intent(in) :: extra_row, extra_col
logical, dimension(:), intent(inout) :: got_this_field
character (len=19), intent(in) :: temp_date
character (len=128), intent(in) :: mminlu
character (len=128), dimension(:), intent(inout) :: output_flags
! BUG: Move this constant to misc_definitions_module?
integer, parameter :: BDR_WIDTH = 3
! Local variables
integer :: istatus, iqstatus, fg_idx, idx, idxt, i, j, bottom_top_dim, &
sm1, em1, sm2, em2, sm3, em3, &
sp1, ep1, sp2, ep2, sp3, ep3, &
sd1, ed1, sd2, ed2, sd3, ed3, &
u_idx, bdr_wdth
real :: rx, ry
real :: threshold
logical :: do_gcell_interp
integer, pointer, dimension(:) :: u_levels, v_levels
real, pointer, dimension(:,:) :: halo_slab
real, pointer, dimension(:,:,:) :: real_array
character (len=19) :: output_date
character (len=128) :: cname, title
character (len=MAX_FILENAME_LEN) :: input_name
type (fg_input) :: field, u_field, v_field
type (met_data) :: fg_data
! For this time, we need to process all first-guess filename roots. When we
! hit a root containing a '*', we assume we have hit the end of the list
fg_idx = 1
if (do_const_processing) then
input_name = constants_name(fg_idx)
else
input_name = fg_name(fg_idx)
end if
do while (input_name /= '*')
call mprintf(.true.,STDOUT, ' %s', s1=input_name)
call mprintf(.true.,LOGFILE, 'Getting input fields from %s', s1=input_name)
! Do a first pass through this fg source to get all mask fields used
! during interpolation
call get_interp_masks(trim(input_name), do_const_processing, temp_date)
istatus = 0
! Initialize the module for reading in the met fields
call read_met_init(trim(input_name), do_const_processing, temp_date, istatus)
if (istatus == 0) then
! Process all fields and levels from the current file; read_next_met_field()
! will return a non-zero status when there are no more fields to be read.
do while (istatus == 0)
call read_next_met_field(fg_data, istatus)
if (istatus == 0) then
! Find index into fieldname, interp_method, masked, and fill_missing
! of the current field
idxt = num_entries + 1
do idx=1,num_entries
if ((index(fieldname(idx), trim(fg_data%field)) /= 0) .and. &
(len_trim(fieldname(idx)) == len_trim(fg_data%field))) then
got_this_field(idx) = .true.
if (index(input_name,trim(from_input(idx))) /= 0 .or. &
(from_input(idx) == '*' .and. idxt == num_entries + 1)) then
idxt = idx
end if
end if
end do
idx = idxt
if (idx > num_entries) idx = num_entries ! The last entry is a default
! Do we need to rename this field?
if (output_name(idx) /= ' ') then
fg_data%field = output_name(idx)(1:9)
idxt = num_entries + 1
do idx=1,num_entries
if ((index(fieldname(idx), trim(fg_data%field)) /= 0) .and. &
(len_trim(fieldname(idx)) == len_trim(fg_data%field))) then
got_this_field(idx) = .true.
if (index(input_name,trim(from_input(idx))) /= 0 .or. &
(from_input(idx) == '*' .and. idxt == num_entries + 1)) then
idxt = idx
end if
end if
end do
idx = idxt
if (idx > num_entries) idx = num_entries ! The last entry is a default
end if
! Do a simple check to see whether this is a global lat/lon dataset
if ( (fg_data%iproj == PROJ_LATLON .or. fg_data%iproj == PROJ_GAUSS) .and. &
abs(fg_data%nx * fg_data%deltalon - 360.) < 0.0001) then
bdr_wdth = BDR_WIDTH
allocate(halo_slab(1-BDR_WIDTH:fg_data%nx+BDR_WIDTH,1:fg_data%ny))
halo_slab(1:fg_data%nx, 1:fg_data%ny) = &
fg_data%slab(1:fg_data%nx, 1:fg_data%ny)
halo_slab(1-BDR_WIDTH:0, 1:fg_data%ny) = &
fg_data%slab(fg_data%nx-BDR_WIDTH+1:fg_data%nx, 1:fg_data%ny)
halo_slab(fg_data%nx+1:fg_data%nx+BDR_WIDTH, 1:fg_data%ny) = &
fg_data%slab(1:BDR_WIDTH, 1:fg_data%ny)
deallocate(fg_data%slab)
else
bdr_wdth = 0
halo_slab => fg_data%slab
nullify(fg_data%slab)
end if
call mprintf(.true.,LOGFILE,'Processing %s at level %f.',s1=fg_data%field,f1=fg_data%xlvl)
call push_source_projection(fg_data%iproj, fg_data%xlonc, fg_data%truelat1, &
fg_data%truelat2, fg_data%dx, fg_data%dy, fg_data%deltalat, &
fg_data%deltalon, fg_data%starti, fg_data%startj, &
fg_data%startlat, fg_data%startlon, earth_radius=fg_data%earth_radius*1000.)
! Initialize fg_input structure to store the field
field%header%version = 1
field%header%date = fg_data%hdate//' '
if (do_const_processing) then
field%header%time_dependent = .false.
else
field%header%time_dependent = .true.
end if
field%header%forecast_hour = fg_data%xfcst
field%header%fg_source = 'FG'
field%header%field = ' '
field%header%field(1:9) = fg_data%field
field%header%units = ' '
field%header%units(1:25) = fg_data%units
field%header%description = ' '
field%header%description(1:46) = fg_data%desc
call get_z_dim_name(fg_data%field,field%header%vertical_coord)
field%header%vertical_level = nint(fg_data%xlvl)
field%header%array_order = 'XY '
field%header%is_wind_grid_rel = fg_data%is_wind_grid_rel
field%header%array_has_missing_values = .false.
nullify(field%r_arr)
nullify(field%valid_mask)
nullify(field%modified_mask)
if (output_this_field(idx) .and. flag_in_output(idx) /= ' ') then
output_flags(idx) = flag_in_output(idx)
end if
! If we should not output this field, just list it as a mask field
if (output_this_field(idx)) then
field%header%mask_field = .false.
else
field%header%mask_field = .true.
end if
!
! Before actually doing any interpolation to the model grid, we must check
! whether we will be using the average_gcell interpolator that averages all
! source points in each model grid cell
!
do_gcell_interp = .false.
if (index(interp_method(idx),'average_gcell') /= 0) then
call get_gcell_threshold(interp_method(idx), threshold, istatus)
if (istatus == 0) then
if (fg_data%dx == 0. .and. fg_data%dy == 0. .and. &
fg_data%deltalat /= 0. .and. fg_data%deltalon /= 0.) then
fg_data%dx = abs(fg_data%deltalon)
fg_data%dy = abs(fg_data%deltalat)
else
! BUG: Need to more correctly handle dx/dy in meters.
fg_data%dx = fg_data%dx / 111000. ! Convert meters to approximate degrees
fg_data%dy = fg_data%dy / 111000.
end if
if (gridtype == 'C' .or. gridtype == 'A') then
if (threshold*max(fg_data%dx,fg_data%dy)*111. <= max(dom_dx,dom_dy)/1000.) &
do_gcell_interp = .true.
else if (gridtype == 'E' .or. gridtype == 'B') then
if (threshold*max(fg_data%dx,fg_data%dy) <= max(dom_dx,dom_dy)) &
do_gcell_interp = .true.
end if
end if
end if
! Interpolate to U staggering
if (output_stagger(idx) == U) then
call storage_query_field(field, iqstatus)
if (iqstatus == 0) then
call storage_get_field(field, iqstatus)
call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// &
' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl))
if (associated(field%modified_mask)) then
call bitarray_destroy(field%modified_mask)
nullify(field%modified_mask)
end if
else
allocate(field%valid_mask)
call bitarray_create(field%valid_mask, we_mem_stag_e-we_mem_stag_s+1, sn_mem_e-sn_mem_s+1)
end if
! Save a copy of the fg_input structure for the U field so that we can find it later
if (is_u_field(idx)) call dup(field, u_field)
allocate(field%modified_mask)
call bitarray_create(field%modified_mask, we_mem_stag_e-we_mem_stag_s+1, sn_mem_e-sn_mem_s+1)
if (do_const_processing .or. field%header%time_dependent) then
call interp_met_field(input_name, fg_data%field, U, M, &
field, xlat_u, xlon_u, we_mem_stag_s, we_mem_stag_e, sn_mem_s, sn_mem_e, &
halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
field%modified_mask)
else
call mprintf(.true.,INFORM,' - already processed this field from constant file.')
end if
! Interpolate to V staggering
else if (output_stagger(idx) == V) then
call storage_query_field(field, iqstatus)
if (iqstatus == 0) then
call storage_get_field(field, iqstatus)
call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// &
' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl))
if (associated(field%modified_mask)) then
call bitarray_destroy(field%modified_mask)
nullify(field%modified_mask)
end if
else
allocate(field%valid_mask)
call bitarray_create(field%valid_mask, we_mem_e-we_mem_s+1, sn_mem_stag_e-sn_mem_stag_s+1)
end if
! Save a copy of the fg_input structure for the V field so that we can find it later
if (is_v_field(idx)) call dup(field, v_field)
allocate(field%modified_mask)
call bitarray_create(field%modified_mask, we_mem_e-we_mem_s+1, sn_mem_stag_e-sn_mem_stag_s+1)
if (do_const_processing .or. field%header%time_dependent) then
call interp_met_field(input_name, fg_data%field, V, M, &
field, xlat_v, xlon_v, we_mem_s, we_mem_e, sn_mem_stag_s, sn_mem_stag_e, &
halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
field%modified_mask)
else
call mprintf(.true.,INFORM,' - already processed this field from constant file.')
end if
! Interpolate to VV staggering
else if (output_stagger(idx) == VV) then
call storage_query_field(field, iqstatus)
if (iqstatus == 0) then
call storage_get_field(field, iqstatus)
call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// &
' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl))
if (associated(field%modified_mask)) then
call bitarray_destroy(field%modified_mask)
nullify(field%modified_mask)
end if
else
allocate(field%valid_mask)
call bitarray_create(field%valid_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1)
end if
! Save a copy of the fg_input structure for the U field so that we can find it later
if (is_u_field(idx)) call dup(field, u_field)
! Save a copy of the fg_input structure for the V field so that we can find it later
if (is_v_field(idx)) call dup(field, v_field)
allocate(field%modified_mask)
call bitarray_create(field%modified_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1)
if (do_const_processing .or. field%header%time_dependent) then
!here
call interp_met_field(input_name, fg_data%field, VV, M, &
field, xlat_v, xlon_v, we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
field%modified_mask,landmask)
else
call mprintf(.true.,INFORM,' - already processed this field from constant file.')
end if
! All other fields interpolated to M staggering for C grid, H staggering for E grid
else
call storage_query_field(field, iqstatus)
if (iqstatus == 0) then
call storage_get_field(field, iqstatus)
call mprintf((iqstatus /= 0),ERROR,'Queried field %s at level %i and found it,'// &
' but could not get data.',s1=fg_data%field,i1=nint(fg_data%xlvl))
if (associated(field%modified_mask)) then
call bitarray_destroy(field%modified_mask)
nullify(field%modified_mask)
end if
else
allocate(field%valid_mask)
call bitarray_create(field%valid_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1)
end if
allocate(field%modified_mask)
call bitarray_create(field%modified_mask, we_mem_e-we_mem_s+1, sn_mem_e-sn_mem_s+1)
if (do_const_processing .or. field%header%time_dependent) then
if (gridtype == 'C' .or. gridtype == 'A') then
call interp_met_field(input_name, fg_data%field, M, M, &
field, xlat, xlon, we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
field%modified_mask, landmask)
else if (gridtype == 'E' .or. gridtype == 'B') then
call interp_met_field(input_name, fg_data%field, HH, M, &
field, xlat, xlon, we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
field%modified_mask, landmask)
end if
else
call mprintf(.true.,INFORM,' - already processed this field from constant file.')
end if
end if
call bitarray_merge(field%valid_mask, field%modified_mask)
deallocate(halo_slab)
! Store the interpolated field
call storage_put_field(field)
call pop_source_projection()
end if
end do
call read_met_close()
call push_source_projection(fg_data%iproj, fg_data%xlonc, fg_data%truelat1, &
fg_data%truelat2, fg_data%dx, fg_data%dy, fg_data%deltalat, &
fg_data%deltalon, fg_data%starti, fg_data%startj, &
fg_data%startlat, fg_data%startlon, earth_radius=fg_data%earth_radius*1000.)
!
! If necessary, rotate winds to earth-relative for this fg source
!
call storage_get_levels(u_field, u_levels)
call storage_get_levels(v_field, v_levels)
if (associated(u_levels) .and. associated(v_levels)) then
u_idx = 1
do u_idx = 1, size(u_levels)
u_field%header%vertical_level = u_levels(u_idx)
call storage_get_field(u_field, istatus)
v_field%header%vertical_level = v_levels(u_idx)
call storage_get_field(v_field, istatus)
if (associated(u_field%modified_mask) .and. &
associated(v_field%modified_mask)) then
if (u_field%header%is_wind_grid_rel) then
if (gridtype == 'C') then
call map_to_met(u_field%r_arr, u_field%modified_mask, &
v_field%r_arr, v_field%modified_mask, &
we_mem_stag_s, sn_mem_s, &
we_mem_stag_e, sn_mem_e, &
we_mem_s, sn_mem_stag_s, &
we_mem_e, sn_mem_stag_e, &
xlon_u, xlon_v, xlat_u, xlat_v)
else if (gridtype == 'E' .or. gridtype == 'B') then
call map_to_met_nmm(u_field%r_arr, u_field%modified_mask, &
v_field%r_arr, v_field%modified_mask, &
we_mem_s, sn_mem_s, &
we_mem_e, sn_mem_e, &
xlat_v, xlon_v)
else if (gridtype == 'A') then
call map_to_met_nmm(u_field%r_arr, u_field%modified_mask, &
v_field%r_arr, v_field%modified_mask, &
we_mem_s, sn_mem_s, &
we_mem_e, sn_mem_e, &
xlat, xlon)
end if
end if
call bitarray_destroy(u_field%modified_mask)
call bitarray_destroy(v_field%modified_mask)
nullify(u_field%modified_mask)
nullify(v_field%modified_mask)
call storage_put_field(u_field)
call storage_put_field(v_field)
end if
end do
deallocate(u_levels)
deallocate(v_levels)
end if
call pop_source_projection()
else
if (do_const_processing) then
call mprintf(.true.,WARN,'Couldn''t open file %s for input.',s1=input_name)
else
call mprintf(.true.,WARN,'Couldn''t open file %s for input.',s1=trim(input_name)//':'//trim(temp_date))
end if
end if
fg_idx = fg_idx + 1
if (do_const_processing) then
input_name = constants_name(fg_idx)
else
input_name = fg_name(fg_idx)
end if
end do ! while (input_name /= '*')
!
! Rotate winds from earth-relative to grid-relative
!
call storage_get_levels(u_field, u_levels)
call storage_get_levels(v_field, v_levels)
if (associated(u_levels) .and. associated(v_levels)) then
u_idx = 1
do u_idx = 1, size(u_levels)
u_field%header%vertical_level = u_levels(u_idx)
call storage_get_field(u_field, istatus)
v_field%header%vertical_level = v_levels(u_idx)
call storage_get_field(v_field, istatus)
if (gridtype == 'C') then
call met_to_map(u_field%r_arr, u_field%valid_mask, &
v_field%r_arr, v_field%valid_mask, &
we_mem_stag_s, sn_mem_s, &
we_mem_stag_e, sn_mem_e, &
we_mem_s, sn_mem_stag_s, &
we_mem_e, sn_mem_stag_e, &
xlon_u, xlon_v, xlat_u, xlat_v)
else if (gridtype == 'E' .or. gridtype == 'B') then
call met_to_map_nmm(u_field%r_arr, u_field%valid_mask, &
v_field%r_arr, v_field%valid_mask, &
we_mem_s, sn_mem_s, &
we_mem_e, sn_mem_e, &
xlat_v, xlon_v)
else if (gridtype == 'A') then
call met_to_map_nmm(u_field%r_arr, u_field%valid_mask, &
v_field%r_arr, v_field%valid_mask, &
we_mem_s, sn_mem_s, &
we_mem_e, sn_mem_e, &
xlat, xlon)
end if
end do
deallocate(u_levels)
deallocate(v_levels)
end if
if (do_const_processing) return
!
! Now that we have all degribbed fields, we build a 3-d pressure field, and fill in any
! missing levels in the other 3-d fields
!
call mprintf(.true.,LOGFILE,'Filling missing levels.')
call fill_missing_levels(output_flags)
call mprintf(.true.,LOGFILE,'Creating derived fields.')
call create_derived_fields(gridtype, fg_data%hdate, fg_data%xfcst, &
we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
we_mem_stag_s, we_mem_stag_e, sn_mem_stag_s, sn_mem_stag_e, &
got_this_field, output_flags)
!
! Check that every mandatory field was found in input data
!
do i=1,num_entries
if (is_mandatory(i) .and. .not. got_this_field(i)) then
call mprintf(.true.,ERROR,'The mandatory field %s was not found in any input data.',s1=fieldname(i))
end if
end do
!
! Before we begin to write fields, if debug_level is set high enough, we
! write a table of which fields are available at which levels to the
! metgrid.log file
!
! call storage_print_fields()
! call find_missing_values()
!
! All of the processing is now done for this time period for this domain;
! now we simply output every field from the storage module.
!
title = 'OUTPUT FROM METGRID'
! Initialize the output module for this domain and time
call mprintf(.true.,LOGFILE,'Initializing output module.')
output_date = temp_date
if (len_trim(temp_date) == 13) then
output_date(14:19) = ':00:00'
else if (len_trim(temp_date) == 16) then
output_date(17:19) = ':00'
end if
call output_init(n, title, output_date, gridtype, dyn_opt, &
corner_lats, corner_lons, &
we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
we_patch_s, we_patch_e, sn_patch_s, sn_patch_e, &
we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
extra_col, extra_row)
call get_bottom_top_dim(bottom_top_dim)
! First write out global attributes
call mprintf(.true.,LOGFILE,'Writing global attributes to output.')
call write_global_attrs(title, output_date, gridtype, dyn_opt, west_east_dim, &
south_north_dim, bottom_top_dim, &
we_patch_s, we_patch_e, we_patch_stag_s, we_patch_stag_e, &
sn_patch_s, sn_patch_e, sn_patch_stag_s, sn_patch_stag_e, &
map_proj, mminlu, is_water, is_ice, is_urban, i_soilwater, &
grid_id, parent_id, i_parent_start, &
j_parent_start, i_parent_end, j_parent_end, dom_dx, dom_dy, &
cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, &
truelat2, parent_grid_ratio, corner_lats, corner_lons, output_flags, num_entries)
call reset_next_field()
istatus = 0
! Now loop over all output fields, writing each to the output module
do while (istatus == 0)
call get_next_output_field(cname, real_array, &
sm1, em1, sm2, em2, sm3, em3, istatus)
if (istatus == 0) then
call mprintf(.true.,LOGFILE,'Writing field %s to output.',s1=cname)
call write_field(sm1, em1, sm2, em2, sm3, em3, &
cname, output_date, real_array)
deallocate(real_array)
end if
end do
call mprintf(.true.,LOGFILE,'Closing output file.')
call output_close()
! Free up memory used by met fields for this valid time
call storage_delete_all_td()
end subroutine process_single_met_time
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Name: get_interp_masks
!
! Purpose:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine get_interp_masks(fg_prefix, is_constants, fg_date)
use interp_option_module
use read_met_module
use storage_module
implicit none
! Arguments
logical, intent(in) :: is_constants
character (len=*), intent(in) :: fg_prefix, fg_date
! BUG: Move this constant to misc_definitions_module?
integer, parameter :: BDR_WIDTH = 3
! Local variables
integer :: i, istatus, idx, idxt
type (fg_input) :: mask_field
type (met_data) :: fg_data
istatus = 0
call read_met_init(fg_prefix, is_constants, fg_date, istatus)
do while (istatus == 0)
call read_next_met_field(fg_data, istatus)
if (istatus == 0) then
! Find out which METGRID.TBL entry goes with this field
idxt = num_entries + 1
do idx=1,num_entries
if ((index(fieldname(idx), trim(fg_data%field)) /= 0) .and. &
(len_trim(fieldname(idx)) == len_trim(fg_data%field))) then
if (index(fg_prefix,trim(from_input(idx))) /= 0 .or. &
(from_input(idx) == '*' .and. idxt == num_entries + 1)) then
idxt = idx
end if
end if
end do
idx = idxt
if (idx > num_entries) idx = num_entries ! The last entry is a default
! Do we need to rename this field?
if (output_name(idx) /= ' ') then
fg_data%field = output_name(idx)(1:9)
idxt = num_entries + 1
do idx=1,num_entries
if ((index(fieldname(idx), trim(fg_data%field)) /= 0) .and. &
(len_trim(fieldname(idx)) == len_trim(fg_data%field))) then
if (index(fg_prefix,trim(from_input(idx))) /= 0 .or. &
(from_input(idx) == '*' .and. idxt == num_entries + 1)) then
idxt = idx
end if
end if
end do
idx = idxt
if (idx > num_entries) idx = num_entries ! The last entry is a default
end if
do i=1,num_entries
if (interp_mask(i) /= ' ' .and. (trim(interp_mask(i)) == trim(fg_data%field))) then
mask_field%header%version = 1
mask_field%header%date = ' '
mask_field%header%date = fg_date
if (is_constants) then
mask_field%header%time_dependent = .false.
else
mask_field%header%time_dependent = .true.
end if
mask_field%header%mask_field = .true.
mask_field%header%forecast_hour = 0.
mask_field%header%fg_source = 'degribbed met data'
mask_field%header%field = trim(fg_data%field)//'.mask'
mask_field%header%units = '-'
mask_field%header%description = '-'
mask_field%header%vertical_coord = 'none'
mask_field%header%vertical_level = 1
mask_field%header%array_order = 'XY'
mask_field%header%dim1(1) = 1
mask_field%header%dim1(2) = fg_data%nx
mask_field%header%dim2(1) = 1
mask_field%header%dim2(2) = fg_data%ny
mask_field%header%is_wind_grid_rel = .true.
mask_field%header%array_has_missing_values = .false.
mask_field%map%stagger = M
! Do a simple check to see whether this is a global lat/lon dataset
if ( (fg_data%iproj == PROJ_LATLON .or. fg_data%iproj == PROJ_GAUSS) .and. &
abs(fg_data%nx * fg_data%deltalon - 360.) < 0.0001) then
allocate(mask_field%r_arr(1-BDR_WIDTH:fg_data%nx+BDR_WIDTH,1:fg_data%ny))
mask_field%r_arr(1:fg_data%nx, 1:fg_data%ny) = &
fg_data%slab(1:fg_data%nx, 1:fg_data%ny)
mask_field%r_arr(1-BDR_WIDTH:0, 1:fg_data%ny) = &
fg_data%slab(fg_data%nx-BDR_WIDTH+1:fg_data%nx, 1:fg_data%ny)
mask_field%r_arr(fg_data%nx+1:fg_data%nx+BDR_WIDTH, 1:fg_data%ny) = &
fg_data%slab(1:BDR_WIDTH, 1:fg_data%ny)
else
allocate(mask_field%r_arr(1:fg_data%nx,1:fg_data%ny))
mask_field%r_arr = fg_data%slab
end if
nullify(mask_field%valid_mask)
nullify(mask_field%modified_mask)
call storage_put_field(mask_field)
exit
end if
end do
if (associated(fg_data%slab)) deallocate(fg_data%slab)
end if
end do
call read_met_close()
end subroutine get_interp_masks
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Name: interp_met_field
!
! Purpose:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine interp_met_field(input_name, short_fieldnm, ifieldstagger, istagger, &
field, xlat, xlon, sm1, em1, sm2, em2, &
slab, minx, maxx, miny, maxy, bdr, do_gcell_interp, &
new_pts, landmask)
use bitarray_module
use interp_module
use interp_option_module
use llxy_module
use misc_definitions_module
use storage_module
implicit none
! Arguments
integer, intent(in) :: ifieldstagger, istagger, &
sm1, em1, sm2, em2, &
minx, maxx, miny, maxy, bdr
real, dimension(minx:maxx,miny:maxy), intent(in) :: slab
real, dimension(sm1:em1,sm2:em2), intent(in) :: xlat, xlon
real, dimension(sm1:em1,sm2:em2), intent(in), optional :: landmask
logical, intent(in) :: do_gcell_interp
character (len=9), intent(in) :: short_fieldnm
character (len=MAX_FILENAME_LEN), intent(in) :: input_name
type (fg_input), intent(inout) :: field
type (bitarray), intent(inout) :: new_pts
! Local variables
integer :: i, j, idx, idxt, orig_selected_proj, interp_mask_status, &
interp_land_mask_status, interp_water_mask_status
integer, pointer, dimension(:) :: interp_array
real :: rx, ry, temp
real, pointer, dimension(:,:) :: data_count
type (fg_input) :: mask_field, mask_water_field, mask_land_field
! Find index into fieldname, interp_method, masked, and fill_missing
! of the current field
idxt = num_entries + 1
do idx=1,num_entries
if ((index(fieldname(idx), trim(short_fieldnm)) /= 0) .and. &
(len_trim(fieldname(idx)) == len_trim(short_fieldnm))) then
if (index(input_name,trim(from_input(idx))) /= 0 .or. &
(from_input(idx) == '*' .and. idxt == num_entries + 1)) then
idxt = idx
end if
end if
end do
idx = idxt
if (idx > num_entries) then
call mprintf(.true.,WARN,'Entry in METGRID.TBL not found for field %s. '// &
'Default options will be used for this field!', s1=short_fieldnm)
idx = num_entries ! The last entry is a default
end if
field%header%dim1(1) = sm1
field%header%dim1(2) = em1
field%header%dim2(1) = sm2
field%header%dim2(2) = em2
field%map%stagger = ifieldstagger
if (.not. associated(field%r_arr)) then
allocate(field%r_arr(sm1:em1,sm2:em2))
end if
interp_mask_status = 1
interp_land_mask_status = 1
interp_water_mask_status = 1
if (interp_mask(idx) /= ' ') then
mask_field%header%version = 1
mask_field%header%forecast_hour = 0.
mask_field%header%field = trim(interp_mask(idx))//'.mask'
mask_field%header%vertical_coord = 'none'
mask_field%header%vertical_level = 1
call storage_get_field(mask_field, interp_mask_status)
end if
if (interp_land_mask(idx) /= ' ') then
mask_land_field%header%version = 1
mask_land_field%header%forecast_hour = 0.
mask_land_field%header%field = trim(interp_land_mask(idx))//'.mask'
mask_land_field%header%vertical_coord = 'none'
mask_land_field%header%vertical_level = 1
call storage_get_field(mask_land_field, interp_land_mask_status)
end if
if (interp_water_mask(idx) /= ' ') then
mask_water_field%header%version = 1
mask_water_field%header%forecast_hour = 0.
mask_water_field%header%field = trim(interp_water_mask(idx))//'.mask'
mask_water_field%header%vertical_coord = 'none'
mask_water_field%header%vertical_level = 1
call storage_get_field(mask_water_field, interp_water_mask_status)
end if
interp_array => interp_array_from_string(interp_method(idx))
!
! Interpolate using average_gcell interpolation method
!
if (do_gcell_interp) then
allocate(data_count(sm1:em1,sm2:em2))
data_count = 0.
if (interp_mask_status == 0) then
call accum_continuous(slab, &
minx+bdr, maxx-bdr, miny, maxy, 1, 1, &
field%r_arr, data_count, &
sm1, em1, sm2, em2, 1, 1, &
istagger, &
new_pts, NAN, interp_mask_val(idx), mask_field%r_arr)
else
call accum_continuous(slab, &
minx+bdr, maxx-bdr, miny, maxy, 1, 1, &
field%r_arr, data_count, &
sm1, em1, sm2, em2, 1, 1, &
istagger, &
new_pts, NAN, -1.) ! The -1 is the maskval, but since we
! we do not give an optional mask, no
! no need to worry about -1s in data
end if
orig_selected_proj = iget_selected_domain()
call select_domain(SOURCE_PROJ)
do j=sm2,em2
do i=sm1,em1
if (present(landmask)) then
if (landmask(i,j) /= masked(idx)) then
if (data_count(i,j) > 0.) then
field%r_arr(i,j) = field%r_arr(i,j) / data_count(i,j)
else
if (interp_mask_status == 0) then
temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
minx, maxx, miny, maxy, bdr, missing_value(idx), &
mask_val=interp_mask_val(idx), mask_field=mask_field%r_arr)
else
temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
minx, maxx, miny, maxy, bdr, missing_value(idx))
end if
if (temp /= missing_value(idx)) then
field%r_arr(i,j) = temp
call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
end if
end if
else
field%r_arr(i,j) = fill_missing(idx)
call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
end if
if (.not. bitarray_test(new_pts, i-sm1+1, j-sm2+1) .and. &
.not. bitarray_test(field%valid_mask, i-sm1+1, j-sm2+1)) then
field%r_arr(i,j) = fill_missing(idx)
end if
else
if (data_count(i,j) > 0.) then
field%r_arr(i,j) = field%r_arr(i,j) / data_count(i,j)
else
if (interp_mask_status == 0) then
temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
minx, maxx, miny, maxy, bdr, missing_value(idx), &
mask_val=interp_mask_val(idx), mask_field=mask_field%r_arr)
else
temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
minx, maxx, miny, maxy, bdr, missing_value(idx))
end if
if (temp /= missing_value(idx)) then
field%r_arr(i,j) = temp
call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
end if
end if
if (.not. bitarray_test(new_pts, i-sm1+1, j-sm2+1) .and. &
.not. bitarray_test(field%valid_mask, i-sm1+1, j-sm2+1)) then
field%r_arr(i,j) = fill_missing(idx)
end if
end if
end do
end do
call select_domain(orig_selected_proj)
deallocate(data_count)
!
! No average_gcell interpolation method
!
else
orig_selected_proj = iget_selected_domain()
call select_domain(SOURCE_PROJ)
do j=sm2,em2
do i=sm1,em1
if (present(landmask)) then
if (masked(idx) == MASKED_BOTH) then
if (landmask(i,j) == 0) then ! WATER POINT
if (interp_land_mask_status == 0) then
temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
minx, maxx, miny, maxy, bdr, missing_value(idx), &
mask_val=interp_land_mask_val(idx), mask_field=mask_land_field%r_arr)
else
temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
minx, maxx, miny, maxy, bdr, missing_value(idx))
end if
else if (landmask(i,j) == 1) then ! LAND POINT
if (interp_water_mask_status == 0) then
temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
minx, maxx, miny, maxy, bdr, missing_value(idx), &
mask_val=interp_water_mask_val(idx), &
mask_field=mask_water_field%r_arr)
else
temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
minx, maxx, miny, maxy, bdr, missing_value(idx))
end if
end if
else if (landmask(i,j) /= masked(idx)) then
if (interp_mask_status == 0) then
temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
minx, maxx, miny, maxy, bdr, missing_value(idx), &
mask_val=interp_mask_val(idx), mask_field=mask_field%r_arr)
else
temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
minx, maxx, miny, maxy, bdr, missing_value(idx))
end if
else
temp = missing_value(idx)
end if
else
if (interp_mask_status == 0) then
temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
minx, maxx, miny, maxy, bdr, missing_value(idx), &
mask_val=interp_mask_val(idx), mask_field=mask_field%r_arr)
else
temp = interp_to_latlon(xlat(i,j), xlon(i,j), istagger, interp_array, slab, &
minx, maxx, miny, maxy, bdr, missing_value(idx))
end if
end if
if (temp /= missing_value(idx)) then
field%r_arr(i,j) = temp
call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
else if (landmask(i,j) == masked(idx)) then
field%r_arr(i,j) = fill_missing(idx)
call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
end if
if (.not. bitarray_test(new_pts, i-sm1+1, j-sm2+1) .and. &
.not. bitarray_test(field%valid_mask, i-sm1+1, j-sm2+1)) then
field%r_arr(i,j) = fill_missing(idx)
end if
end do
end do
call select_domain(orig_selected_proj)
end if
deallocate(interp_array)
end subroutine interp_met_field
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Name: interp_to_latlon
!
! Purpose:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function interp_to_latlon(rlat, rlon, istagger, interp_method_list, slab, &
minx, maxx, miny, maxy, bdr, source_missing_value, &
mask_field, mask_val)
use interp_module
use llxy_module
implicit none
! Arguments
integer, intent(in) :: minx, maxx, miny, maxy, bdr, istagger
integer, dimension(:), intent(in) :: interp_method_list
real, intent(in) :: rlat, rlon, source_missing_value
real, dimension(minx:maxx,miny:maxy), intent(in) :: slab
real, intent(in), optional :: mask_val
real, dimension(minx:maxx,miny:maxy), intent(in), optional :: mask_field
! Return value
real :: interp_to_latlon
! Local variables
real :: rx, ry
interp_to_latlon = source_missing_value
call lltoxy(rlat, rlon, rx, ry, istagger)
if (rx >= minx+bdr-0.5 .and. rx <= maxx-bdr+0.5) then
if (present(mask_field) .and. present(mask_val)) then
interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, 1, 1, source_missing_value, &
interp_method_list, 1, mask_val, mask_field)
else
interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, 1, 1, source_missing_value, &
interp_method_list, 1)
end if
else
interp_to_latlon = source_missing_value
end if
if (interp_to_latlon == source_missing_value) then
! Try a lon in the range 0. to 360.; all lons in the xlon
! array should be in the range -180. to 180.
if (rlon < 0.) then
call lltoxy(rlat, rlon+360., rx, ry, istagger)
if (rx >= minx+bdr-0.5 .and. rx <= maxx-bdr+0.5) then
if (present(mask_field) .and. present(mask_val)) then
interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, &
1, 1, source_missing_value, &
interp_method_list, 1, mask_val, mask_field)
else
interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, &
1, 1, source_missing_value, &
interp_method_list, 1)
end if
else
interp_to_latlon = source_missing_value
end if
end if
end if
return
end function interp_to_latlon
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Name: get_bottom_top_dim
!
! Purpose:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine get_bottom_top_dim(bottom_top_dim)
use interp_option_module
use list_module
use storage_module
implicit none
! Arguments
integer, intent(out) :: bottom_top_dim
! Local variables
integer :: i, j
integer, pointer, dimension(:) :: field_levels
character (len=32) :: z_dim
type (fg_input), pointer, dimension(:) :: headers
type (list) :: temp_levels
! Initialize a list to store levels that are found for 3-d fields
call list_init(temp_levels)
! Get a list of all time-dependent fields (given by their headers) from
! the storage module
call storage_get_td_headers(headers)
!
! Given headers of all fields, we first build a list of all possible levels
! for 3-d met fields (excluding sea-level, though).
!
do i=1,size(headers)
call get_z_dim_name(headers(i)%header%field, z_dim)
! We only want to consider 3-d met fields
if (z_dim(1:18) == 'num_metgrid_levels') then
! Find out what levels the current field has
call storage_get_levels(headers(i), field_levels)
do j=1,size(field_levels)
! If this level has not yet been encountered, add it to our list
if (.not. list_search(temp_levels, ikey=field_levels(j), ivalue=field_levels(j))) then
if (field_levels(j) /= 201300) then
call list_insert(temp_levels, ikey=field_levels(j), ivalue=field_levels(j))
end if
end if
end do
deallocate(field_levels)
end if
end do
bottom_top_dim = list_length(temp_levels)
call list_destroy(temp_levels)
deallocate(headers)
end subroutine get_bottom_top_dim
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Name: fill_missing_levels
!
! Purpose:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine fill_missing_levels(output_flags)
use interp_option_module
use list_module
use module_debug
use module_mergesort
use storage_module
implicit none
! Arguments
character (len=128), dimension(:), intent(inout) :: output_flags
! Local variables
integer :: i, ii, j, ix, jx, k, lower, upper, temp, istatus
integer, pointer, dimension(:) :: union_levels, field_levels
real, pointer, dimension(:) :: r_union_levels
character (len=128) :: clevel
type (fg_input) :: lower_field, upper_field, new_field, search_field
type (fg_input), pointer, dimension(:) :: headers, all_headers
type (list) :: temp_levels
type (list_item), pointer, dimension(:) :: keys
! Initialize a list to store levels that are found for 3-d fields
call list_init(temp_levels)
! Get a list of all fields (given by their headers) from the storage module
call storage_get_td_headers(headers)
call storage_get_all_headers(all_headers)
!
! Given headers of all fields, we first build a list of all possible levels
! for 3-d met fields (excluding sea-level, though).
!
do i=1,size(headers)
! Find out what levels the current field has
call storage_get_levels(headers(i), field_levels)
do j=1,size(field_levels)
! If this level has not yet been encountered, add it to our list
if (.not. list_search(temp_levels, ikey=field_levels(j), ivalue=field_levels(j))) then
if (field_levels(j) /= 201300) then
call list_insert(temp_levels, ikey=field_levels(j), ivalue=field_levels(j))
end if
end if
end do
deallocate(field_levels)
end do
if (list_length(temp_levels) > 0) then
!
! With all possible levels stored in a list, get an array of levels, sorted
! in decreasing order
!
i = 0
allocate(union_levels(list_length(temp_levels)))
do while (list_length(temp_levels) > 0)
i = i + 1
call list_get_first_item(temp_levels, ikey=union_levels(i), ivalue=temp)
end do
call mergesort(union_levels, 1, size(union_levels))
allocate(r_union_levels(size(union_levels)))
do i=1,size(union_levels)
r_union_levels(i) = real(union_levels(i))
end do
!
! With a sorted, complete list of levels, we need
! to go back and fill in missing levels for each 3-d field
!
do i=1,size(headers)
!
! Find entry in METGRID.TBL for this field, if one exists; if it does, then the
! entry may tell us how to get values for the current field at the missing level
!
do ii=1,num_entries
if (fieldname(ii) == headers(i)%header%field) exit
end do
if (ii <= num_entries) then
call dup(headers(i),new_field)
nullify(new_field%valid_mask)
nullify(new_field%modified_mask)
call fill_field(new_field, ii, output_flags, r_union_levels)
end if
end do
deallocate(union_levels)
deallocate(r_union_levels)
deallocate(headers)
call storage_get_td_headers(headers)
!
! Now we may need to vertically interpolate to missing values in 3-d fields
!
do i=1,size(headers)
call storage_get_levels(headers(i), field_levels)
! If this isn't a 3-d array, nothing to do
if (size(field_levels) > 1) then
do k=1,size(field_levels)
call dup(headers(i),search_field)
search_field%header%vertical_level = field_levels(k)
call storage_get_field(search_field,istatus)
if (istatus == 0) then
JLOOP: do jx=search_field%header%dim2(1),search_field%header%dim2(2)
ILOOP: do ix=search_field%header%dim1(1),search_field%header%dim1(2)
if (.not. bitarray_test(search_field%valid_mask, &
ix-search_field%header%dim1(1)+1, &
jx-search_field%header%dim2(1)+1)) then
call dup(search_field, lower_field)
do lower=k-1,1,-1
lower_field%header%vertical_level = field_levels(lower)
call storage_get_field(lower_field,istatus)
if (bitarray_test(lower_field%valid_mask, &
ix-search_field%header%dim1(1)+1, &
jx-search_field%header%dim2(1)+1)) &
exit
end do
call dup(search_field, upper_field)
do upper=k+1,size(field_levels)
upper_field%header%vertical_level = field_levels(upper)
call storage_get_field(upper_field,istatus)
if (bitarray_test(upper_field%valid_mask, &
ix-search_field%header%dim1(1)+1, &
jx-search_field%header%dim2(1)+1)) &
exit
end do
if (upper <= size(field_levels) .and. lower >= 1) then
search_field%r_arr(ix,jx) = real(abs(field_levels(upper)-field_levels(k))) &
/ real(abs(field_levels(upper)-field_levels(lower))) &
* lower_field%r_arr(ix,jx) &
+ real(abs(field_levels(k)-field_levels(lower))) &
/ real(abs(field_levels(upper)-field_levels(lower))) &
* upper_field%r_arr(ix,jx)
call bitarray_set(search_field%valid_mask, &
ix-search_field%header%dim1(1)+1, &
jx-search_field%header%dim2(1)+1)
end if
end if
end do ILOOP
end do JLOOP
else
call mprintf(.true.,ERROR, &
'This is bad, could not get %s at level %i.', &
s1=trim(search_field%header%field), i1=field_levels(k))
end if
end do
end if
deallocate(field_levels)
end do
end if
call list_destroy(temp_levels)
deallocate(all_headers)
deallocate(headers)
end subroutine fill_missing_levels
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Name: create_derived_fields
!
! Purpose:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine create_derived_fields(arg_gridtype, hdate, xfcst, &
we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
we_mem_stag_s, we_mem_stag_e, sn_mem_stag_s, sn_mem_stag_e, &
created_this_field, output_flags)
use interp_option_module
use list_module
use module_mergesort
use storage_module
implicit none
! Arguments
integer, intent(in) :: we_mem_s, we_mem_e, sn_mem_s, sn_mem_e, &
we_mem_stag_s, we_mem_stag_e, sn_mem_stag_s, sn_mem_stag_e
real, intent(in) :: xfcst
logical, dimension(:), intent(inout) :: created_this_field
character (len=1), intent(in) :: arg_gridtype
character (len=24), intent(in) :: hdate
character (len=128), dimension(:), intent(inout) :: output_flags
! Local variables
integer :: idx, i, j, istatus
type (fg_input) :: field
! Initialize fg_input structure to store the field
field%header%version = 1
field%header%date = hdate//' '
field%header%time_dependent = .true.
field%header%mask_field = .false.
field%header%forecast_hour = xfcst
field%header%fg_source = 'Derived from FG'
field%header%field = ' '
field%header%units = ' '
field%header%description = ' '
field%header%vertical_level = 0
field%header%array_order = 'XY '
field%header%is_wind_grid_rel = .true.
field%header%array_has_missing_values = .false.
nullify(field%r_arr)
nullify(field%valid_mask)
nullify(field%modified_mask)
!
! Check each entry in METGRID.TBL to see whether it is a derive field
!
do idx=1,num_entries
if (is_derived_field(idx)) then
created_this_field(idx) = .true.
call mprintf(.true.,INFORM,'Going to create the field %s',s1=fieldname(idx))
! Intialize more fields in storage structure
field%header%field = fieldname(idx)
call get_z_dim_name(fieldname(idx),field%header%vertical_coord)
field%map%stagger = output_stagger(idx)
if (arg_gridtype == 'E') then
field%header%dim1(1) = we_mem_s
field%header%dim1(2) = we_mem_e
field%header%dim2(1) = sn_mem_s
field%header%dim2(2) = sn_mem_e
elseif (arg_gridtype == 'B' .or. arg_gridtype == 'A') then
field%header%dim1(1) = we_mem_s
field%header%dim1(2) = we_mem_e
field%header%dim2(1) = sn_mem_s
field%header%dim2(2) = sn_mem_e
else if (arg_gridtype == 'C') then
if (output_stagger(idx) == M) then
field%header%dim1(1) = we_mem_s
field%header%dim1(2) = we_mem_e
field%header%dim2(1) = sn_mem_s
field%header%dim2(2) = sn_mem_e
else if (output_stagger(idx) == U) then
field%header%dim1(1) = we_mem_stag_s
field%header%dim1(2) = we_mem_stag_e
field%header%dim2(1) = sn_mem_s
field%header%dim2(2) = sn_mem_e
else if (output_stagger(idx) == V) then
field%header%dim1(1) = we_mem_s
field%header%dim1(2) = we_mem_e
field%header%dim2(1) = sn_mem_stag_s
field%header%dim2(2) = sn_mem_stag_e
end if
end if
call fill_field(field, idx, output_flags)
end if
end do
end subroutine create_derived_fields
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Name: fill_field
!
! Purpose:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine fill_field(field, idx, output_flags, all_level_list)
use interp_option_module
use list_module
use module_mergesort
use storage_module
implicit none
! Arguments
integer, intent(in) :: idx
type (fg_input), intent(inout) :: field
character (len=128), dimension(:), intent(inout) :: output_flags
real, dimension(:), intent(in), optional :: all_level_list
! Local variables
integer :: i, j, istatus, isrclevel
integer, pointer, dimension(:) :: all_list
real :: rfillconst, rlevel, rsrclevel
type (fg_input) :: query_field
type (list_item), pointer, dimension(:) :: keys
character (len=128) :: asrcname
logical :: filled_all_lev
filled_all_lev = .false.
!
! Get a list of all levels to be filled for this field
!
keys => list_get_keys(fill_lev_list(idx))
do i=1,list_length(fill_lev_list(idx))
!
! First handle a specification for levels "all"
!
if (trim(keys(i)%ckey) == 'all') then
! We only want to fill all levels if we haven't already filled "all" of them
if (.not. filled_all_lev) then
filled_all_lev = .true.
query_field%header%time_dependent = .true.
query_field%header%field = ' '
nullify(query_field%r_arr)
nullify(query_field%valid_mask)
nullify(query_field%modified_mask)
! See if we are filling this level with a constant
call get_constant_fill_lev(keys(i)%cvalue, rfillconst, istatus)
if (istatus == 0) then
if (present(all_level_list)) then
do j=1,size(all_level_list)
call create_level(field, real(all_level_list(j)), idx, output_flags, rfillconst=rfillconst)
end do
else
query_field%header%field = level_template(idx)
nullify(all_list)
call storage_get_levels(query_field, all_list)
if (associated(all_list)) then
do j=1,size(all_list)
call create_level(field, real(all_list(j)), idx, output_flags, rfillconst=rfillconst)
end do
deallocate(all_list)
end if
end if
! Else see if we are filling this level with a constant equal
! to the value of the level
else if (trim(keys(i)%cvalue) == 'vertical_index') then
if (present(all_level_list)) then
do j=1,size(all_level_list)
call create_level(field, real(all_level_list(j)), idx, output_flags, &
rfillconst=real(all_level_list(j)))
end do
else
query_field%header%field = level_template(idx)
nullify(all_list)
call storage_get_levels(query_field, all_list)
if (associated(all_list)) then
do j=1,size(all_list)
call create_level(field, real(all_list(j)), idx, output_flags, rfillconst=real(all_list(j)))
end do
deallocate(all_list)
end if
end if
! Else, we assume that it is a field from which we are copying levels
else
if (present(all_level_list)) then
do j=1,size(all_level_list)
call create_level(field, real(all_level_list(j)), idx, output_flags, &
asrcname=keys(i)%cvalue, rsrclevel=real(all_level_list(j)))
end do
else
query_field%header%field = keys(i)%cvalue ! Use same levels as source field, not level_template
nullify(all_list)
call storage_get_levels(query_field, all_list)
if (associated(all_list)) then
do j=1,size(all_list)
call create_level(field, real(all_list(j)), idx, output_flags, &
asrcname=keys(i)%cvalue, rsrclevel=real(all_list(j)))
end do
deallocate(all_list)
else
! If the field doesn't have any levels (or does not exist) then we have not
! really filled all levels at this point.
filled_all_lev = .false.
end if
end if
end if
end if
!
! Handle individually specified levels
!
else
read(keys(i)%ckey,*) rlevel
! See if we are filling this level with a constant
call get_constant_fill_lev(keys(i)%cvalue, rfillconst, istatus)
if (istatus == 0) then
call create_level(field, rlevel, idx, output_flags, rfillconst=rfillconst)
! Otherwise, we are filling from another level
else
call get_fill_src_level(keys(i)%cvalue, asrcname, isrclevel)
rsrclevel = real(isrclevel)
call create_level(field, rlevel, idx, output_flags, &
asrcname=asrcname, rsrclevel=rsrclevel)
end if
end if
end do
if (associated(keys)) deallocate(keys)
end subroutine fill_field
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Name: create_level
!
! Purpose:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine create_level(field_template, rlevel, idx, output_flags, &
rfillconst, asrcname, rsrclevel)
use storage_module
use interp_option_module
implicit none
! Arguments
type (fg_input), intent(inout) :: field_template
real, intent(in) :: rlevel
integer, intent(in) :: idx
character (len=128), dimension(:), intent(inout) :: output_flags
real, intent(in), optional :: rfillconst, rsrclevel
character (len=128), intent(in), optional :: asrcname
! Local variables
integer :: i, j, istatus
integer :: sm1,em1,sm2,em2
type (fg_input) :: query_field
!
! Check to make sure optional arguments are sane
!
if (present(rfillconst) .and. (present(asrcname) .or. present(rsrclevel))) then
call mprintf(.true.,ERROR,'A call to create_level() cannot be given specifications '// &
'for both a constant fill value and a source level.')
else if ((present(asrcname) .and. .not. present(rsrclevel)) .or. &
(.not. present(asrcname) .and. present(rsrclevel))) then
call mprintf(.true.,ERROR,'Neither or both of optional arguments asrcname and '// &
'rsrclevel must be specified to subroutine create_level().')
else if (.not. present(rfillconst) .and. &
.not. present(asrcname) .and. &
.not. present(rsrclevel)) then
call mprintf(.true.,ERROR,'A call to create_level() must be given either a specification '// &
'for a constant fill value or a source level.')
end if
query_field%header%time_dependent = .true.
query_field%header%field = field_template%header%field
query_field%header%vertical_level = rlevel
nullify(query_field%r_arr)
nullify(query_field%valid_mask)
nullify(query_field%modified_mask)
call storage_query_field(query_field, istatus)
if (istatus == 0) then
call mprintf(.true.,INFORM,'%s at level %f already exists; leaving it alone.', &
s1=field_template%header%field, f1=rlevel)
return
end if
sm1 = field_template%header%dim1(1)
em1 = field_template%header%dim1(2)
sm2 = field_template%header%dim2(1)
em2 = field_template%header%dim2(2)
!
! Handle constant fill value case
!
if (present(rfillconst)) then
field_template%header%vertical_level = rlevel
allocate(field_template%r_arr(sm1:em1,sm2:em2))
allocate(field_template%valid_mask)
allocate(field_template%modified_mask)
call bitarray_create(field_template%valid_mask, em1-sm1+1, em2-sm2+1)
call bitarray_create(field_template%modified_mask, em1-sm1+1, em2-sm2+1)
field_template%r_arr = rfillconst
do j=sm2,em2
do i=sm1,em1
call bitarray_set(field_template%valid_mask, i-sm1+1, j-sm2+1)
end do
end do
call storage_put_field(field_template)
if (output_this_field(idx) .and. flag_in_output(idx) /= ' ') then
output_flags(idx) = flag_in_output(idx)
end if
!
! Handle source field and source level case
!
else if (present(asrcname) .and. present(rsrclevel)) then
query_field%header%field = ' '
query_field%header%field = asrcname
query_field%header%vertical_level = rsrclevel
! Check to see whether the requested source field exists at the requested level
call storage_query_field(query_field, istatus)
if (istatus == 0) then
! Read in requested field at requested level
call storage_get_field(query_field, istatus)
if ((query_field%header%dim1(1) /= field_template%header%dim1(1)) .or. &
(query_field%header%dim1(2) /= field_template%header%dim1(2)) .or. &
(query_field%header%dim2(1) /= field_template%header%dim2(1)) .or. &
(query_field%header%dim2(2) /= field_template%header%dim2(2))) then
call mprintf(.true.,ERROR,'Dimensions for %s do not match those of %s. This is '// &
'probably because the staggerings of the fields do not match.', &
s1=query_field%header%field, s2=field_template%header%field)
end if
field_template%header%vertical_level = rlevel
allocate(field_template%r_arr(sm1:em1,sm2:em2))
allocate(field_template%valid_mask)
allocate(field_template%modified_mask)
call bitarray_create(field_template%valid_mask, em1-sm1+1, em2-sm2+1)
call bitarray_create(field_template%modified_mask, em1-sm1+1, em2-sm2+1)
field_template%r_arr = query_field%r_arr
! We should retain information about which points in the field are valid
do j=sm2,em2
do i=sm1,em1
if (bitarray_test(query_field%valid_mask, i-sm1+1, j-sm2+1)) then
call bitarray_set(field_template%valid_mask, i-sm1+1, j-sm2+1)
end if
end do
end do
call storage_put_field(field_template)
if (output_this_field(idx) .and. flag_in_output(idx) /= ' ') then
output_flags(idx) = flag_in_output(idx)
end if
else
call mprintf(.true.,INFORM,'Couldn''t find %s at level %f to fill level %f of %s.', &
s1=asrcname,f1=rsrclevel,f2=rlevel,s2=field_template%header%field)
end if
end if
end subroutine create_level
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Name: accum_continuous
!
! Purpose: Sum up all of the source data points whose nearest neighbor in the
! model grid is the specified model grid point.
!
! NOTE: When processing the source tile, those source points that are
! closest to a different model grid point will be added to the totals for
! such grid points; thus, an entire source tile will be processed at a time.
! This routine really processes for all model grid points that are
! within a source tile, and not just for a single grid point.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine accum_continuous(src_array, &
src_min_x, src_max_x, src_min_y, src_max_y, src_min_z, src_max_z, &
dst_array, n, &
start_i, end_i, start_j, end_j, start_k, end_k, &
istagger, &
new_pts, msgval, maskval, mask_array)
use bitarray_module
use misc_definitions_module
implicit none
! Arguments
integer, intent(in) :: start_i, end_i, start_j, end_j, start_k, end_k, istagger, &
src_min_x, src_max_x, src_min_y, src_max_y, src_min_z, src_max_z
real, intent(in) :: maskval, msgval
real, dimension(src_min_x:src_max_x, src_min_y:src_max_y, src_min_z:src_max_z), intent(in) :: src_array
real, dimension(src_min_x:src_max_x, src_min_y:src_max_y), intent(in), optional :: mask_array
real, dimension(start_i:end_i, start_j:end_j, start_k:end_k), intent(inout) :: dst_array, n
type (bitarray), intent(inout) :: new_pts
! Local variables
integer :: i, j
integer, pointer, dimension(:,:,:) :: where_maps_to
allocate(where_maps_to(src_min_x:src_max_x,src_min_y:src_max_y,2))
do i=src_min_x,src_max_x
do j=src_min_y,src_max_y
where_maps_to(i,j,1) = NOT_PROCESSED
end do
end do
call process_continuous_block(src_array, where_maps_to, &
src_min_x, src_min_y, src_min_z, src_max_x, src_max_y, src_max_z, &
src_min_x, src_min_y, src_min_z, &
src_max_x, src_max_y, src_max_z, &
dst_array, n, start_i, end_i, start_j, end_j, start_k, end_k, &
istagger, &
new_pts, msgval, maskval, mask_array)
deallocate(where_maps_to)
end subroutine accum_continuous
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Name: process_continuous_block
!
! Purpose: To recursively process a subarray of continuous data, adding the
! points in a block to the sum for their nearest grid point. The nearest
! neighbor may be estimated in some cases; for example, if the four corners
! of a subarray all have the same nearest grid point, all elements in the
! subarray are added to that grid point.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
recursive subroutine process_continuous_block(tile_array, where_maps_to, &
src_min_x, src_min_y, src_min_z, src_max_x, src_max_y, src_max_z, &
min_i, min_j, min_k, max_i, max_j, max_k, &
dst_array, n, &
start_x, end_x, start_y, end_y, start_z, end_z, &
istagger, &
new_pts, msgval, maskval, mask_array)
use bitarray_module
use llxy_module
use misc_definitions_module
implicit none
! Arguments
integer, intent(in) :: min_i, min_j, min_k, max_i, max_j, max_k, &
src_min_x, src_min_y, src_min_z, src_max_x, src_max_y, src_max_z, &
start_x, end_x, start_y, end_y, start_z, end_z, istagger
integer, dimension(src_min_x:src_max_x,src_min_y:src_max_y,2), intent(inout) :: where_maps_to
real, intent(in) :: maskval, msgval
real, dimension(src_min_x:src_max_x,src_min_y:src_max_y,src_min_z:src_max_z), intent(in) :: tile_array
real, dimension(src_min_x:src_max_x,src_min_y:src_max_y), intent(in), optional :: mask_array
real, dimension(start_x:end_x,start_y:end_y,start_z:end_z), intent(inout) :: dst_array, n
type (bitarray), intent(inout) :: new_pts
! Local variables
integer :: orig_selected_domain, x_dest, y_dest, i, j, k, center_i, center_j
real :: lat_corner, lon_corner, rx, ry
! Compute the model grid point that the corners of the rectangle to be
! processed map to
! Lower-left corner
if (where_maps_to(min_i,min_j,1) == NOT_PROCESSED) then
orig_selected_domain = iget_selected_domain()
call select_domain(SOURCE_PROJ)
call xytoll(real(min_i), real(min_j), lat_corner, lon_corner, istagger)
call select_domain(1)
call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
call select_domain(orig_selected_domain)
if (real(start_x) <= rx .and. rx <= real(end_x) .and. &
real(start_y) <= ry .and. ry <= real(end_y)) then
where_maps_to(min_i,min_j,1) = nint(rx)
where_maps_to(min_i,min_j,2) = nint(ry)
else
where_maps_to(min_i,min_j,1) = OUTSIDE_DOMAIN
end if
end if
! Upper-left corner
if (where_maps_to(min_i,max_j,1) == NOT_PROCESSED) then
orig_selected_domain = iget_selected_domain()
call select_domain(SOURCE_PROJ)
call xytoll(real(min_i), real(max_j), lat_corner, lon_corner, istagger)
call select_domain(1)
call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
call select_domain(orig_selected_domain)
if (real(start_x) <= rx .and. rx <= real(end_x) .and. &
real(start_y) <= ry .and. ry <= real(end_y)) then
where_maps_to(min_i,max_j,1) = nint(rx)
where_maps_to(min_i,max_j,2) = nint(ry)
else
where_maps_to(min_i,max_j,1) = OUTSIDE_DOMAIN
end if
end if
! Upper-right corner
if (where_maps_to(max_i,max_j,1) == NOT_PROCESSED) then
orig_selected_domain = iget_selected_domain()
call select_domain(SOURCE_PROJ)
call xytoll(real(max_i), real(max_j), lat_corner, lon_corner, istagger)
call select_domain(1)
call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
call select_domain(orig_selected_domain)
if (real(start_x) <= rx .and. rx <= real(end_x) .and. &
real(start_y) <= ry .and. ry <= real(end_y)) then
where_maps_to(max_i,max_j,1) = nint(rx)
where_maps_to(max_i,max_j,2) = nint(ry)
else
where_maps_to(max_i,max_j,1) = OUTSIDE_DOMAIN
end if
end if
! Lower-right corner
if (where_maps_to(max_i,min_j,1) == NOT_PROCESSED) then
orig_selected_domain = iget_selected_domain()
call select_domain(SOURCE_PROJ)
call xytoll(real(max_i), real(min_j), lat_corner, lon_corner, istagger)
call select_domain(1)
call lltoxy(lat_corner, lon_corner, rx, ry, istagger)
call select_domain(orig_selected_domain)
if (real(start_x) <= rx .and. rx <= real(end_x) .and. &
real(start_y) <= ry .and. ry <= real(end_y)) then
where_maps_to(max_i,min_j,1) = nint(rx)
where_maps_to(max_i,min_j,2) = nint(ry)
else
where_maps_to(max_i,min_j,1) = OUTSIDE_DOMAIN
end if
end if
! If all four corners map to same model grid point, accumulate the
! entire rectangle
if (where_maps_to(min_i,min_j,1) == where_maps_to(min_i,max_j,1) .and. &
where_maps_to(min_i,min_j,1) == where_maps_to(max_i,max_j,1) .and. &
where_maps_to(min_i,min_j,1) == where_maps_to(max_i,min_j,1) .and. &
where_maps_to(min_i,min_j,2) == where_maps_to(min_i,max_j,2) .and. &
where_maps_to(min_i,min_j,2) == where_maps_to(max_i,max_j,2) .and. &
where_maps_to(min_i,min_j,2) == where_maps_to(max_i,min_j,2) .and. &
where_maps_to(min_i,min_j,1) /= OUTSIDE_DOMAIN) then
x_dest = where_maps_to(min_i,min_j,1)
y_dest = where_maps_to(min_i,min_j,2)
! If this grid point was already given a value from higher-priority source data,
! there is nothing to do.
! if (.not. bitarray_test(processed_pts, x_dest-start_x+1, y_dest-start_y+1)) then
! If this grid point has never been given a value by this level of source data,
! initialize the point
if (.not. bitarray_test(new_pts, x_dest-start_x+1, y_dest-start_y+1)) then
do k=min_k,max_k
dst_array(x_dest,y_dest,k) = 0.
end do
end if
! Sum all the points whose nearest neighbor is this grid point
if (present(mask_array)) then
do i=min_i,max_i
do j=min_j,max_j
do k=min_k,max_k
! Ignore masked/missing values in the source data
if ((tile_array(i,j,k) /= msgval) .and. &
(mask_array(i,j) /= maskval)) then
dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
end if
end do
end do
end do
else
do i=min_i,max_i
do j=min_j,max_j
do k=min_k,max_k
! Ignore masked/missing values in the source data
if ((tile_array(i,j,k) /= msgval)) then
dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
end if
end do
end do
end do
end if
! end if
! Rectangle is a square of four points, and we can simply deal with each of the points
else if (((max_i - min_i + 1) <= 2) .and. ((max_j - min_j + 1) <= 2)) then
do i=min_i,max_i
do j=min_j,max_j
x_dest = where_maps_to(i,j,1)
y_dest = where_maps_to(i,j,2)
if (x_dest /= OUTSIDE_DOMAIN) then
! if (.not. bitarray_test(processed_pts, x_dest-start_x+1, y_dest-start_y+1)) then
if (.not. bitarray_test(new_pts, x_dest-start_x+1, y_dest-start_y+1)) then
do k=min_k,max_k
dst_array(x_dest,y_dest,k) = 0.
end do
end if
if (present(mask_array)) then
do k=min_k,max_k
! Ignore masked/missing values
if ((tile_array(i,j,k) /= msgval) .and. &
(mask_array(i,j) /= maskval)) then
dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
end if
end do
else
do k=min_k,max_k
! Ignore masked/missing values
if ((tile_array(i,j,k) /= msgval)) then
dst_array(x_dest,y_dest,k) = dst_array(x_dest,y_dest,k) + tile_array(i,j,k)
n(x_dest,y_dest,k) = n(x_dest,y_dest,k) + 1.0
call bitarray_set(new_pts, x_dest-start_x+1, y_dest-start_y+1)
end if
end do
end if
! end if
end if
end do
end do
! Not all corners map to the same grid point, and the rectangle contains more than
! four points
else
center_i = (max_i + min_i)/2
center_j = (max_j + min_j)/2
! Recursively process lower-left rectangle
call process_continuous_block(tile_array, where_maps_to, &
src_min_x, src_min_y, src_min_z, &
src_max_x, src_max_y, src_max_z, &
min_i, min_j, min_k, &
center_i, center_j, max_k, &
dst_array, n, &
start_x, end_x, start_y, end_y, start_z, end_z, &
istagger, &
new_pts, msgval, maskval, mask_array)
if (center_i < max_i) then
! Recursively process lower-right rectangle
call process_continuous_block(tile_array, where_maps_to, &
src_min_x, src_min_y, src_min_z, &
src_max_x, src_max_y, src_max_z, &
center_i+1, min_j, min_k, max_i, &
center_j, max_k, &
dst_array, n, &
start_x, end_x, start_y, &
end_y, start_z, end_z, &
istagger, &
new_pts, msgval, maskval, mask_array)
end if
if (center_j < max_j) then
! Recursively process upper-left rectangle
call process_continuous_block(tile_array, where_maps_to, &
src_min_x, src_min_y, src_min_z, &
src_max_x, src_max_y, src_max_z, &
min_i, center_j+1, min_k, center_i, &
max_j, max_k, &
dst_array, n, &
start_x, end_x, start_y, &
end_y, start_z, end_z, &
istagger, &
new_pts, msgval, maskval, mask_array)
end if
if (center_i < max_i .and. center_j < max_j) then
! Recursively process upper-right rectangle
call process_continuous_block(tile_array, where_maps_to, &
src_min_x, src_min_y, src_min_z, &
src_max_x, src_max_y, src_max_z, &
center_i+1, center_j+1, min_k, max_i, &
max_j, max_k, &
dst_array, n, &
start_x, end_x, start_y, &
end_y, start_z, end_z, &
istagger, &
new_pts, msgval, maskval, mask_array)
end if
end if
end subroutine process_continuous_block
end module process_domain_module