/WPS/metgrid/src/process_domain_module.F
FORTRAN Legacy | 2839 lines | 1753 code | 456 blank | 630 comment | 397 complexity | 2d29f8bcd4bc0d9d1f57c097777cb2a6 MD5 | raw file
Possible License(s): AGPL-1.0
- 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_lake, 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, sub_x, sub_y, num_land_cat, process_bdy_width
- real :: cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
- dom_dx, dom_dy, pole_lat, pole_lon
- 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 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)
-
- n_times = idiff / interval_seconds
-
- ! 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 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, num_land_cat, is_water, is_lake, 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, sub_x, sub_y, &
- cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
- pole_lat, pole_lon, dom_dx, dom_dy, landmask, xlat, xlon, xlat_u, xlon_u, &
- xlat_v, xlon_v, corner_lats, corner_lons, title)
- 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, num_land_cat, is_water, is_lake, 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, &
- pole_lat, pole_lon, parent_grid_ratio, sub_x, sub_y, &
- corner_lats, corner_lons, output_flags, 0)
- !
- ! 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
- if (t > 0) then
- process_bdy_width = process_only_bdy
- else
- process_bdy_width = 0
- end if
-
- 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, num_land_cat, is_water, is_lake, 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, &
- pole_lat, pole_lon, parent_grid_ratio, sub_x, sub_y, &
- corner_lats, corner_lons, td_output_flags, process_bdy_width)
-
- 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, num_land_cat, is_water, is_lake, 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, sub_x, sub_y, &
- cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
- pole_lat, pole_lon, 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
- use module_debug
- 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_lake, 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, sub_x, sub_y, num_land_cat
- real, pointer, dimension(:,:) :: landmask
- real, intent(inout) :: cen_lat, moad_cen_lat, cen_lon, stand_lon, truelat1, truelat2, &
- dom_dx, dom_dy, pole_lat, pole_lon
- 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
- integer :: we_mem_subgrid_s, we_mem_subgrid_e, &
- sn_mem_subgrid_s, sn_mem_subgrid_e
- integer :: we_patch_subgrid_s, we_patch_subgrid_e, &
- sn_patch_subgrid_s, sn_patch_subgrid_e
- 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
- logical :: is_subgrid
- ! 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, num_land_cat, &
- is_water, is_lake, 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, pole_lat, pole_lon, parent_grid_ratio, &
- corner_lats, corner_lons, sub_x, sub_y)
- 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
- 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') 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
- ! 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 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, is_subgrid, &
- 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
- we_mem_subgrid_s = (we_mem_s + HALO_WIDTH*lh_mult - 1)*sub_x - HALO_WIDTH*lh_mult + 1
- we_mem_subgrid_e = (we_mem_e + (1-rh_mult) - HALO_WIDTH*rh_mult )*sub_x + HALO_WIDTH*rh_mult
- we_patch_subgrid_s = (we_patch_s - 1)*sub_x + 1
- we_patch_subgrid_e = (we_patch_e + (1-rh_mult) )*sub_x
- sn_mem_subgrid_s = (sn_mem_s + HALO_WIDTH*bh_mult - 1)*sub_y - HALO_WIDTH*bh_mult + 1
- sn_mem_subgrid_e = (sn_mem_e + (1-th_mult) - HALO_WIDTH*th_mult )*sub_y + HALO_WIDTH*th_mult
- sn_patch_subgrid_s = (sn_patch_s - 1)*sub_y + 1
- sn_patch_subgrid_e = (sn_patch_e + (1-th_mult) )*sub_y
-
- ! 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%sr_x=sub_x
- field%header%sr_y=sub_y
- field%header%is_subgrid=is_subgrid
- 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 (is_subgrid) then
- field%map%stagger = M
- field%header%dim1(1) = we_mem_subgrid_s
- field%header%dim1(2) = we_mem_subgrid_e
- field%header%dim2(1) = sn_mem_subgrid_s
- field%header%dim2(2) = sn_mem_subgrid_e
- else 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') 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
- end if
- allocate(field%valid_mask)
- if (is_subgrid) then
- allocate(field%r_arr(we_mem_subgrid_s:we_mem_subgrid_e,&
- sn_mem_subgrid_s:sn_mem_subgrid_e))
- field%r_arr(we_patch_subgrid_s:we_patch_subgrid_e,sn_patch_subgrid_s:sn_patch_subgrid_e) = &
- real_array(sp1:ep1,sp2:ep2,k)
- call exchange_halo_r(field%r_arr, &
- we_mem_subgrid_s, we_mem_subgrid_e, sn_mem_subgrid_s, sn_mem_subgrid_e, 1, 1, &
- we_patch_subgrid_s, we_patch_subgrid_e, sn_patch_subgrid_s, sn_patch_subgrid_e, 1, 1)
- call bitarray_create(field%valid_mask, &
- (we_mem_subgrid_e-we_mem_subgrid_s)+1, &
- (sn_mem_subgrid_e-sn_mem_subgrid_s)+1)
- do j=1,(sn_mem_subgrid_e-sn_mem_subgrid_s)+1
- do i=1,(we_mem_subgrid_e-we_mem_subgrid_s)+1
- call bitarray_set(field%valid_mask, i, j)
- end do
- end do
- else 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, num_land_cat, is_water, is_lake, 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, &
- pole_lat, pole_lon, parent_grid_ratio, sub_x, sub_y, &
- corner_lats, corner_lons, output_flags, process_bdy_width)
-
- 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_lake, 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, sub_x, sub_y, num_land_cat, &
- process_bdy_width
- ! 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, pole_lat, pole_lon
- 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
- integer :: nmet_flags
- integer :: num_metgrid_soil_levs
- integer, pointer, dimension(:) :: soil_levels
- 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
- character (len=128), allocatable, dimension(:) :: met_flags
- 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 dataset
- ! Note that we do not currently support regional Gaussian grids
- if ((fg_data%iproj == PROJ_LATLON .and. abs(fg_data%nx * fg_data%deltalon - 360.) < 0.0001) &
- .or. (fg_data%iproj == PROJ_GAUSS)) 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%sr_x = sub_x
- field%header%sr_y = sub_y
- field%header%is_subgrid = .false.
- 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') 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') 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, &
- we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
- halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
- field%modified_mask, process_bdy_width)
- 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, &
- we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
- halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
- field%modified_mask, process_bdy_width)
- 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
- 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, &
- we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
- halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
- field%modified_mask, process_bdy_width)
- 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') 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, &
- we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
- halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
- field%modified_mask, process_bdy_width, landmask)
-
- else if (gridtype == 'E') 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, &
- we_domain_s, we_domain_e, sn_domain_s, sn_domain_e, &
- halo_slab, 1-bdr_wdth, fg_data%nx+bdr_wdth, 1, fg_data%ny, bdr_wdth, do_gcell_interp, &
- field%modified_mask, process_bdy_width, 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') 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)
- 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') 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)
- 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, and then we check to see if any fields are not
- ! completely covered with data.
- !
- 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 V3.3'
-
- ! Initialize the output module for this domain and time
- call mprintf(.true.,LOGFILE,'Initializing output module.')
- output_date = temp_date
- if ( .not. nocolons ) then
- 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
- else
- 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
- endif
- 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)
- ! Add in a flag to tell real that we have seven new msf fields
- nmet_flags = num_entries + 1
- allocate(met_flags(nmet_flags))
- do i=1,num_entries
- met_flags(i) = output_flags(i)
- end do
- if (gridtype == 'C') then
- met_flags(num_entries+1) = 'FLAG_MF_XY'
- else
- met_flags(num_entries+1) = ' '
- end if
- ! Find out how many soil levels or layers we have; this assumes a field named ST
- field % header % field = 'ST'
- nullify(soil_levels)
- call storage_get_levels(field, soil_levels)
- if (.not. associated(soil_levels)) then
- field % header % field = 'SOILT'
- nullify(soil_levels)
- call storage_get_levels(field, soil_levels)
- end if
- if (.not. associated(soil_levels)) then
- field % header % field = 'STC_WPS'
- nullify(soil_levels)
- call storage_get_levels(field, soil_levels)
- end if
- if (associated(soil_levels)) then
- num_metgrid_soil_levs = size(soil_levels)
- deallocate(soil_levels)
- else
- num_metgrid_soil_levs = 0
- end if
-
- ! 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, num_land_cat, &
- is_water, is_lake, 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, &
- pole_lat, pole_lon, parent_grid_ratio, sub_x, sub_y, &
- corner_lats, corner_lons, num_metgrid_soil_levs, &
- met_flags, nmet_flags, process_bdy_width)
- deallocate(met_flags)
-
- 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%sr_x = 1
- mask_field%header%sr_y = 1
- mask_field%header%is_subgrid = .false.
- 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
- ! Note that we do not currently support regional Gaussian grids
- if ((fg_data%iproj == PROJ_LATLON .and. abs(fg_data%nx * fg_data%deltalon - 360.) < 0.0001) &
- .or. (fg_data%iproj == PROJ_GAUSS)) 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, &
- sd1, ed1, sd2, ed2, &
- slab, minx, maxx, miny, maxy, bdr, do_gcell_interp, &
- new_pts, process_bdy_width, 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, &
- sd1, ed1, sd2, ed2, &
- minx, maxx, miny, maxy, bdr, &
- process_bdy_width
- 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, process_width
- 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
- if (process_bdy_width == 0) then
- process_width = max(ed1-sd1+1, ed2-sd2+1)
- else
- process_width = process_bdy_width
- ! Add two extra rows/cols to accommodate staggered fields: one extra row/col for
- ! averaging to mass points in real, and one beyond that for averaging during
- ! wind rotation
- if (ifieldstagger /= M) process_width = process_width + 2
- 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, maxx, miny, maxy, 1, 1, bdr, &
- field%r_arr, data_count, &
- sm1, em1, sm2, em2, 1, 1, &
- istagger, &
- new_pts, missing_value(idx), interp_mask_val(idx), interp_mask_relational(idx), mask_field%r_arr)
- else
- call accum_continuous(slab, &
- minx, maxx, miny, maxy, 1, 1, bdr, &
- field%r_arr, data_count, &
- sm1, em1, sm2, em2, 1, 1, &
- istagger, &
- new_pts, missing_value(idx), -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 (abs(i - sd1) >= process_width .and. (abs(i - ed1) >= process_width) .and. &
- abs(j - sd2) >= process_width .and. (abs(j - ed2) >= process_width)) then
- field%r_arr(i,j) = fill_missing(idx)
- call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
- cycle
- end if
- 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)
- call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
- 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_relational=interp_mask_relational(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)
- ! Assume that if missing fill value is other than default, then user has asked
- ! to fill in any missing values, and we can consider this point to have
- ! received a valid value
- if (fill_missing(idx) /= NAN) call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
- end if
- else
- if (data_count(i,j) > 0.) then
- field%r_arr(i,j) = field%r_arr(i,j) / data_count(i,j)
- call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
- 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_relational=interp_mask_relational(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)
- ! Assume that if missing fill value is other than default, then user has asked
- ! to fill in any missing values, and we can consider this point to have
- ! received a valid value
- if (fill_missing(idx) /= NAN) call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
- 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 (abs(i - sd1) >= process_width .and. (abs(i - ed1) >= process_width) .and. &
- abs(j - sd2) >= process_width .and. (abs(j - ed2) >= process_width)) then
- field%r_arr(i,j) = fill_missing(idx)
- call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
- cycle
- end if
- 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_relational=interp_land_mask_relational(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_relational=interp_water_mask_relational(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_relational=interp_mask_relational(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
- ! No landmask for this field
- 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_relational=interp_mask_relational(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 (present(landmask)) then
- 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
- 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)
- ! Assume that if missing fill value is other than default, then user has asked
- ! to fill in any missing values, and we can consider this point to have
- ! received a valid value
- if (fill_missing(idx) /= NAN) call bitarray_set(new_pts, i-sm1+1, j-sm2+1)
- 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_relational, 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
- character(len=1), intent(in), optional :: mask_relational
- ! 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) .and. present (mask_relational)) then
- interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, 1, 1, source_missing_value, &
- interp_method_list, 1, mask_relational, mask_val, mask_field)
- else 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, maskval=mask_val, mask_array=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) .and. present(mask_relational)) then
- interp_to_latlon = interp_sequence(rx, ry, 1, slab, minx, maxx, miny, maxy, &
- 1, 1, source_missing_value, &
- interp_method_list, 1, mask_relational, mask_val, mask_field)
- else 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, maskval=mask_val, mask_array=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%sr_x = 1
- field%header%sr_y = 1
- field%header%is_subgrid = .false.
- 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
- 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, bdr_width, &
- dst_array, n, &
- start_i, end_i, start_j, end_j, start_k, end_k, &
- istagger, &
- new_pts, msgval, maskval, mask_relational, mask_array, sr_x, sr_y)
-
- 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, bdr_width
- 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(start_i:end_i, start_j:end_j, start_k:end_k), intent(inout) :: dst_array, n
- real, dimension(src_min_x:src_max_x, src_min_y:src_max_y), intent(in), optional :: mask_array
- integer, intent(in), optional :: sr_x, sr_y
- type (bitarray), intent(inout) :: new_pts
- character(len=1), intent(in), optional :: mask_relational
-
- ! Local variables
- integer :: i, j
- integer, pointer, dimension(:,:,:) :: where_maps_to
- real :: rsr_x, rsr_y
- rsr_x = 1.0
- rsr_y = 1.0
- if (present(sr_x)) rsr_x = real(sr_x)
- if (present(sr_y)) rsr_y = real(sr_y)
-
- 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+bdr_width, src_min_y, src_min_z, &
- src_max_x-bdr_width, src_max_y, src_max_z, &
- dst_array, n, start_i, end_i, start_j, end_j, start_k, end_k, &
- istagger, &
- new_pts, rsr_x, rsr_y, msgval, maskval, mask_relational, 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, sr_x, sr_y, msgval, maskval, mask_relational, 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) :: sr_x, sr_y, 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
- character(len=1), intent(in), optional :: mask_relational
-
- ! 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)
- rx = (rx - 1.0)*sr_x + 1.0
- ry = (ry - 1.0)*sr_y + 1.0
- 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)
- rx = (rx - 1.0)*sr_x + 1.0
- ry = (ry - 1.0)*sr_y + 1.0
- 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)
- rx = (rx - 1.0)*sr_x + 1.0
- ry = (ry - 1.0)*sr_y + 1.0
- 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)
- rx = (rx - 1.0)*sr_x + 1.0
- ry = (ry - 1.0)*sr_y + 1.0
- 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) .and. present(mask_relational)) 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) then
- if (mask_relational == ' ' .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)
- else if (mask_relational == '<' .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)
- else if (mask_relational == '>' .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 if
- end do
- end do
- end do
- else 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) .and. present(mask_relational)) then
- do k=min_k,max_k
- ! Ignore masked/missing values
- if (tile_array(i,j,k) /= msgval) then
- if (mask_relational == ' ' .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)
- else if (mask_relational == '<' .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)
- else if (mask_relational == '>' .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 if
- end do
- else 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, sr_x, sr_y, msgval, maskval, mask_relational, 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, sr_x, sr_y, msgval, maskval, mask_relational, 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, sr_x, sr_y, msgval, maskval, mask_relational, 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, sr_x, sr_y, msgval, maskval, mask_relational, mask_array)
- end if
- end if
-
- end subroutine process_continuous_block
- end module process_domain_module