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 ! CWH Initialize local pointer variables nullify(landmask) nullify(xlat) nullify(xlon) nullify(xlat_u) nullify(xlon_u) nullify(xlat_v) nullify(xlon_v) ! 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, subx, suby 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 ! CWH Initialize local pointer variables nullify(real_array) ! 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, subx, suby, & 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 if (subx > 1) then we_mem_subgrid_s = (we_mem_s + HALO_WIDTH*lh_mult - 1)*subx - HALO_WIDTH*lh_mult + 1 we_mem_subgrid_e = (we_mem_e + (1-rh_mult) - HALO_WIDTH*rh_mult )*subx + HALO_WIDTH*rh_mult we_patch_subgrid_s = (we_patch_s - 1)*subx + 1 we_patch_subgrid_e = (we_patch_e + (1-rh_mult) )*subx end if if (suby > 1) then sn_mem_subgrid_s = (sn_mem_s + HALO_WIDTH*bh_mult - 1)*suby - HALO_WIDTH*bh_mult + 1 sn_mem_subgrid_e = (sn_mem_e + (1-th_mult) - HALO_WIDTH*th_mult )*suby + HALO_WIDTH*th_mult sn_patch_subgrid_s = (sn_patch_s - 1)*suby + 1 sn_patch_subgrid_e = (sn_patch_e + (1-th_mult) )*suby end if ! Having read in a field, we write each level individually to the ! storage module; levels will be reassembled later on when they ! are written. do k=sp3,ep3 field%header%sr_x=subx field%header%sr_y=suby 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 (subx > 1 .or. suby > 1) 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 (subx > 1 .or. suby > 1) 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 ! CWH Initialize local pointer variables nullify(soil_levels) nullify(u_levels) nullify(v_levels) nullify(halo_slab) nullify(real_array) ! 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 = 1 field%header%sr_y = 1 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.5.1' ! 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%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 ! CWH Initialize local pointer variables nullify(interp_array) nullify(data_count) ! 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 !CWH Initialize local pointer variables nullify(field_levels) nullify(headers) ! 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 ! CWH Initialize local pointer variables nullify(union_levels) nullify(field_levels) nullify(r_union_levels) nullify(headers) nullify(all_headers) nullify(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%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 !CWH Initialize local pointer variables nullify(all_list) nullify(keys) 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