From fa26d08d8154124ac51514b11bc671eb312e20f8 Mon Sep 17 00:00:00 2001 From: Jonathan Beezley Date: Sat, 20 Nov 2010 23:24:16 -0700 Subject: [PATCH] cleaning subgrid code in wps and fixing a bug when sr_?=1 --- WPS/geogrid/src/output_module.F | 103 ++++++++++++------------ WPS/geogrid/src/process_tile_module.F | 80 +++++++++--------- WPS/geogrid/src/source_data_module.F | 94 +++++++++++---------- WPS/metgrid/src/datatype_module.F | 1 + WPS/metgrid/src/input_module.F | 37 ++++----- WPS/metgrid/src/process_domain_module.F | 41 +++++----- WPS/metgrid/src/storage_module.F | 36 ++++----- 7 files changed, 195 insertions(+), 197 deletions(-) diff --git a/WPS/geogrid/src/output_module.F b/WPS/geogrid/src/output_module.F index 105d0877..70b3d072 100644 --- a/WPS/geogrid/src/output_module.F +++ b/WPS/geogrid/src/output_module.F @@ -28,7 +28,7 @@ module output_module integer :: ndims, istagger integer, dimension(MAX_DIMENSIONS) :: dom_start, mem_start, patch_start integer, dimension(MAX_DIMENSIONS) :: dom_end, mem_end, patch_end - integer :: sr_x, sr_y + logical :: is_subgrid real, pointer, dimension(:,:,:) :: rdata_arr character (len=128), dimension(MAX_DIMENSIONS) :: dimnames @@ -422,6 +422,7 @@ subroutine init_output_fields(nest_num, grid_type, & use storage_module #endif use parallel_module + use gridinfo_module, only : subgrid_ratio_x, subgrid_ratio_y implicit none @@ -445,7 +446,14 @@ subroutine init_output_fields(nest_num, grid_type, & character (len=128) :: memorder, units, description character (len=128), dimension(3) :: dimnames integer :: sr_x, sr_y, istatus - + logical :: subgrid_var + + ! + ! Get local nest subgrid refinement ratio + ! + sr_x=subgrid_ratio_x(nest_num) + sr_y=subgrid_ratio_y(nest_num) + ! ! First find out how many fields there are ! @@ -457,7 +465,7 @@ subroutine init_output_fields(nest_num, grid_type, & call get_next_output_fieldname(nest_num, fieldname, ndims, & min_category, max_category, & istagger, memorder, dimnames, & - units, description, sr_x, sr_y, ifieldstatus) + units, description, subgrid_var, ifieldstatus) if (ifieldstatus == 0) nfields = nfields + 1 end do @@ -473,10 +481,6 @@ subroutine init_output_fields(nest_num, grid_type, & NUM_FIELDS = nfields+NUM_AUTOMATIC_FIELDS allocate(fields(NUM_FIELDS)) - ! Automatic fields will always be on the non-refined grid - sr_x=1 - sr_y=1 - ! ! There are some fields that will always be computed ! Initialize those fields first, followed by all user-specified fields @@ -635,8 +639,7 @@ subroutine init_output_fields(nest_num, grid_type, & fields(i)%dimnames(3) = ' ' fields(i)%mem_order = 'XY' fields(i)%stagger = 'M' - fields(i)%sr_x = 1 - fields(i)%sr_y = 1 + fields(i)%is_subgrid = .false. if (grid_type == 'C') then fields(i)%istagger = M else if (grid_type == 'E') then @@ -774,29 +777,23 @@ subroutine init_output_fields(nest_num, grid_type, & ! Perform adjustments for subgrid lat/lon fields ! do i=23,24 - - call get_subgrid_dim_default(nest_num,fields(i)%dimnames, & - fields(i)%sr_x,fields(i)%sr_y,istatus) - sr_x=fields(i)%sr_x - sr_y=fields(i)%sr_y - - if(sr_x.gt.1.and.sr_y.gt.1)then - fields(i)%mem_start(1) = (fields(i)%mem_start(1) - 1)*sr_x + 1 - fields(i)%patch_start(1) = (fields(i)%patch_start(1) - 1)*sr_x + 1 - fields(i)%dom_start(1) = (fields(i)%dom_start(1) - 1)*sr_x + 1 - - fields(i)%mem_end(1) = (fields(i)%mem_end(1) + 1)*sr_x - fields(i)%patch_end(1) = (fields(i)%patch_end(1) + 1)*sr_x - fields(i)%dom_end(1) = (fields(i)%dom_end(1) + 1)*sr_x - - fields(i)%mem_start(2) = (fields(i)%mem_start(2) - 1)*sr_y + 1 - fields(i)%patch_start(2) = (fields(i)%patch_start(2) - 1)*sr_y + 1 - fields(i)%dom_start(2) = (fields(i)%dom_start(2) - 1)*sr_y + 1 - - fields(i)%mem_end(2) = (fields(i)%mem_end(2) + 1)*sr_y - fields(i)%patch_end(2) = (fields(i)%patch_end(2) + 1)*sr_y - fields(i)%dom_end(2) = (fields(i)%dom_end(2) + 1)*sr_y - end if + fields(i)%is_subgrid=.true. + call get_subgrid_dim_name(fields(i)%dimnames(1:2)) + fields(i)%mem_start(1) = (fields(i)%mem_start(1) - 1)*sr_x + 1 + fields(i)%patch_start(1) = (fields(i)%patch_start(1) - 1)*sr_x + 1 + fields(i)%dom_start(1) = (fields(i)%dom_start(1) - 1)*sr_x + 1 + + fields(i)%mem_end(1) = (fields(i)%mem_end(1) + 1)*sr_x + fields(i)%patch_end(1) = (fields(i)%patch_end(1) + 1)*sr_x + fields(i)%dom_end(1) = (fields(i)%dom_end(1) + 1)*sr_x + + fields(i)%mem_start(2) = (fields(i)%mem_start(2) - 1)*sr_y + 1 + fields(i)%patch_start(2) = (fields(i)%patch_start(2) - 1)*sr_y + 1 + fields(i)%dom_start(2) = (fields(i)%dom_start(2) - 1)*sr_y + 1 + + fields(i)%mem_end(2) = (fields(i)%mem_end(2) + 1)*sr_y + fields(i)%patch_end(2) = (fields(i)%patch_end(2) + 1)*sr_y + fields(i)%dom_end(2) = (fields(i)%dom_end(2) + 1)*sr_y enddo else if (grid_type == 'E') then @@ -829,7 +826,7 @@ subroutine init_output_fields(nest_num, grid_type, & call get_next_output_fieldname(nest_num, fieldname, ndims, & min_category, max_category, & istagger, memorder, dimnames, & - units, description, sr_x, sr_y, ifieldstatus) + units, description, subgrid_var, ifieldstatus) if (ifieldstatus == 0) then !{ @@ -874,22 +871,21 @@ subroutine init_output_fields(nest_num, grid_type, & fields(nfields)%patch_end(2) = end_patch_2 fields(nfields)%patch_end(3) = max_category - fields(nfields)%sr_x=sr_x - fields(nfields)%sr_y=sr_y + fields(nfields)%is_subgrid=subgrid_var - if (extra_col .and. (istagger == U .or. sr_x > 1)) then !{ + if (extra_col .and. (istagger == U .or. subgrid_var)) then !{ fields(nfields)%dom_end(1) = fields(nfields)%dom_end(1) + 1 fields(nfields)%mem_end(1) = fields(nfields)%mem_end(1) + 1 fields(nfields)%patch_end(1) = fields(nfields)%patch_end(1) + 1 - else if ((istagger == U .or. sr_x > 1) .and. my_proc_id == IO_NODE .and. .not. do_tiled_output) then + else if ((istagger == U .or. subgrid_var) .and. my_proc_id == IO_NODE .and. .not. do_tiled_output) then fields(nfields)%dom_end(1)=fields(nfields)%dom_end(1) + 1 end if !} - if (extra_row .and. (istagger == V .or. sr_y > 1)) then !{ + if (extra_row .and. (istagger == V .or. subgrid_var)) then !{ fields(nfields)%dom_end(2) = fields(nfields)%dom_end(2) + 1 fields(nfields)%mem_end(2) = fields(nfields)%mem_end(2) + 1 fields(nfields)%patch_end(2) = fields(nfields)%patch_end(2) + 1 - else if ((istagger == V .or. sr_y > 1) .and. my_proc_id == IO_NODE .and. .not. do_tiled_output) then + else if ((istagger == V .or. subgrid_var) .and. my_proc_id == IO_NODE .and. .not. do_tiled_output) then fields(nfields)%dom_end(2)=fields(nfields)%dom_end(2) + 1 end if !} @@ -905,7 +901,7 @@ subroutine init_output_fields(nest_num, grid_type, & thalo_width = 0 #endif - if (sr_x > 1 .and. sr_y > 1) then + if (subgrid_var) then fields(nfields)%mem_start(1) = (fields(nfields)%mem_start(1) + lhalo_width - 1)*sr_x + 1 - lhalo_width fields(nfields)%patch_start(1) = (fields(nfields)%patch_start(1) - 1)*sr_x + 1 fields(nfields)%dom_start(1) = (fields(nfields)%dom_start(1) - 1)*sr_x + 1 @@ -962,6 +958,7 @@ subroutine write_field(start_mem_i, end_mem_i, & integer, dimension(3) :: sd, ed, sp, ep, sm, em real, pointer, dimension(:,:,:) :: real_dom_array logical :: allocated_real_locally + integer :: is_subgrid allocated_real_locally = .false. @@ -1014,7 +1011,13 @@ subroutine write_field(start_mem_i, end_mem_i, & do i=1,NUM_FIELDS if ( (index(cname, trim(fields(i)%fieldname) ) /= 0) .and. & (len_trim(cname) == len_trim(fields(i)%fieldname)) ) then - + + if (fields(i)%is_subgrid) then + is_subgrid=1 + else + is_subgrid=0 + end if + ! Here, the output array has dimensions of the full grid if it was gathered together ! from all processors if (my_proc_id == IO_NODE .and. nprocs > 1 .and. .not. do_tiled_output) then @@ -1072,10 +1075,8 @@ subroutine write_field(start_mem_i, end_mem_i, & trim(fields(i)%fieldname), trim(fields(i)%descr), istatus) call ext_int_put_var_ti_char(handle, 'stagger', & trim(fields(i)%fieldname), trim(fields(i)%stagger), istatus) - call ext_int_put_var_ti_integer(handle,'sr_x', & - trim(fields(i)%fieldname),(/fields(i)%sr_x/),1, istatus) - call ext_int_put_var_ti_integer(handle,'sr_y', & - trim(fields(i)%fieldname),(/fields(i)%sr_y/),1, istatus) + call ext_int_put_var_ti_integer(handle,'subgrid', & + trim(fields(i)%fieldname),(/is_subgrid/),1, istatus) end if #endif #ifdef IO_NETCDF @@ -1086,10 +1087,8 @@ subroutine write_field(start_mem_i, end_mem_i, & trim(fields(i)%fieldname), trim(fields(i)%descr), istatus) call ext_ncd_put_var_ti_char(handle, 'stagger', & trim(fields(i)%fieldname), trim(fields(i)%stagger), istatus) - call ext_ncd_put_var_ti_integer(handle,'sr_x', & - trim(fields(i)%fieldname),(/fields(i)%sr_x/),1, istatus) - call ext_ncd_put_var_ti_integer(handle,'sr_y', & - trim(fields(i)%fieldname),(/fields(i)%sr_y/),1, istatus) + call ext_ncd_put_var_ti_integer(handle,'subgrid', & + trim(fields(i)%fieldname),(/is_subgrid/),1, istatus) end if #endif #ifdef IO_GRIB1 @@ -1100,10 +1099,8 @@ subroutine write_field(start_mem_i, end_mem_i, & trim(fields(i)%fieldname), trim(fields(i)%descr), istatus) call ext_gr1_put_var_ti_char(handle, 'stagger', & trim(fields(i)%fieldname), trim(fields(i)%stagger), istatus) - call ext_gr1_put_var_ti_integer(handle,'sr_x', & - trim(fields(i)%fieldname),(/fields(i)%sr_x/),1, istatus) - call ext_gr1_put_var_ti_integer(handle,'sr_y', & - trim(fields(i)%fieldname),(/fields(i)%sr_y/),1, istatus) + call ext_gr1_put_var_ti_integer(handle,'subgrid', & + trim(fields(i)%fieldname),(/is_subgrid/),1, istatus) end if #endif end if diff --git a/WPS/geogrid/src/process_tile_module.F b/WPS/geogrid/src/process_tile_module.F index 34e2293d..6ac0c034 100644 --- a/WPS/geogrid/src/process_tile_module.F +++ b/WPS/geogrid/src/process_tile_module.F @@ -32,6 +32,7 @@ subroutine process_tile(which_domain, grid_type, dynopt, & use output_module use smooth_module use source_data_module + use gridinfo_module, only : subgrid_ratio_x, subgrid_ratio_y implicit none @@ -74,8 +75,12 @@ subroutine process_tile(which_domain, grid_type, dynopt, & type (bitarray) :: processed_domain type (hashtable) :: processed_fieldnames character (len=128), dimension(2) :: dimnames - integer :: sub_x, sub_y, idom - + integer :: sr_x, sr_y + logical :: subgrid_var + + sr_x=subgrid_ratio_x(which_domain) + sr_y=subgrid_ratio_y(which_domain) + datestr = '0000-00-00_00:00:00' field_count = 0 mass_flag=1.0 @@ -139,28 +144,18 @@ subroutine process_tile(which_domain, grid_type, dynopt, & end_dom_stag_j = dummy_end_dom_j end if - call get_subgrid_dim_default(which_domain,dimnames,sub_x,sub_y,istatus) - if(sub_x.gt.1)then - start_mem_sub_i = (start_mem_i - 1 ) * sub_x + 1 + start_mem_sub_i = (start_mem_i - 1 ) * sr_x + 1 if(.not.extra_col)then - end_mem_sub_i = (end_mem_i) * sub_x + end_mem_sub_i = (end_mem_i) * sr_x else - end_mem_sub_i = (end_mem_i + 1) * sub_x + end_mem_sub_i = (end_mem_i + 1) * sr_x end if - else - start_mem_sub_i=start_mem_i - end_mem_sub_i=end_mem_i - end if - if(sub_y.gt.1)then - start_mem_sub_j = (start_mem_j -1 ) * sub_y + 1 + + start_mem_sub_j = (start_mem_j -1 ) * sr_y + 1 if(.not.extra_row)then - end_mem_sub_j = (end_mem_j) * sub_y + end_mem_sub_j = (end_mem_j) * sr_y else - end_mem_sub_j = (end_mem_j + 1) * sub_y - end if - else - start_mem_sub_j=start_mem_j - end_mem_sub_j=end_mem_j + end_mem_sub_j = (end_mem_j + 1) * sr_y end if ! Allocate arrays to hold all lat/lon fields; these will persist for the duration of @@ -201,7 +196,7 @@ subroutine process_tile(which_domain, grid_type, dynopt, & call get_lat_lon_fields(clat_array, clon_array, start_mem_i, & start_mem_j, end_mem_i, end_mem_j, M, comp_ll=.true.) call get_lat_lon_fields(xlat_array_subgrid, xlon_array_subgrid, & - start_mem_sub_i, start_mem_sub_j, end_mem_sub_i, end_mem_sub_j, M, sub_x=sub_x, sub_y=sub_y) + start_mem_sub_i, start_mem_sub_j, end_mem_sub_i, end_mem_sub_j, M, sub_x=sr_x, sub_y=sr_y) call get_map_factor(xlat_array_subgrid, xlon_array_subgrid, mapfac_array_x_subgrid, & mapfac_array_y_subgrid, start_mem_sub_i, start_mem_sub_j, end_mem_sub_i, end_mem_sub_j) @@ -712,10 +707,11 @@ subroutine process_tile(which_domain, grid_type, dynopt, & call get_output_stagger(fieldname, istagger, istatus) dimnames(:) = 'null' - call get_subgrid_dim_name(which_domain, fieldname, dimnames, & - sub_x, sub_y, istatus) + + subgrid_var=is_subgrid_var(fieldname) + if( subgrid_var ) call get_subgrid_dim_name(dimnames) - if (istagger == M .and. (sub_x .le. 1) .and. (sub_y .le. 1)) then + if (istagger == M .and. .not. subgrid_var ) then sm1 = start_mem_i em1 = end_mem_i sm2 = start_mem_j @@ -724,7 +720,7 @@ subroutine process_tile(which_domain, grid_type, dynopt, & xlon_ptr => xlon_array mapfac_ptr_x => mapfac_array_m_x mapfac_ptr_y => mapfac_array_m_y - else if ( (sub_x .gt. 1) .or. (sub_y .gt. 1) ) then + else if ( subgrid_var ) then sm1 = start_mem_sub_i em1 = end_mem_sub_i sm2 = start_mem_sub_j @@ -787,14 +783,14 @@ subroutine process_tile(which_domain, grid_type, dynopt, & call mprintf(.true.,LOGFILE,'Processing field %i of %i (%s)', & i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=fieldname) - if ((sub_x > 1) .or. (sub_y > 1)) then + if ( subgrid_var ) then call calc_field(fieldname, field, xlat_ptr, xlon_ptr, istagger, & sm1, em1, sm2, em2, min_level, max_level, & - processed_domain, 1, sr_x=sub_x, sr_y=sub_y) + processed_domain, 1, sr_x=sr_x, sr_y=sr_y) else call calc_field(fieldname, field, xlat_ptr, xlon_ptr, istagger, & sm1, em1, sm2, em2, min_level, max_level, & - processed_domain, 1, landmask=landmask, sr_x=sub_x, sr_y=sub_y) + processed_domain, 1, landmask=landmask) end if ! If user wants to halt when a missing value is found in output field, check now @@ -891,10 +887,10 @@ subroutine process_tile(which_domain, grid_type, dynopt, & if (grid_type == 'C') then call calc_dfdx(field, slp_field, sm1, sm2, min_level, em1, em2, max_level, & - which_domain, mapfac_ptr_x, sub_x) + which_domain, mapfac_ptr_x, subgrid_var) else if (grid_type == 'E') then call calc_dfdx(field, slp_field, sm1, sm2, min_level, em1, em2, max_level, & - which_domain, sr_x=sub_x) + which_domain, subgrid_var=subgrid_var) end if call write_field(sm1, em1, sm2, em2, & min_level, max_level, trim(gradname), datestr, real_array=slp_field) @@ -910,10 +906,10 @@ subroutine process_tile(which_domain, grid_type, dynopt, & if (grid_type == 'C') then call calc_dfdy(field, slp_field, sm1, sm2, min_level, em1, em2, max_level, & - which_domain, mapfac_ptr_y, sub_y) + which_domain, mapfac_ptr_y, subgrid_var) else if (grid_type == 'E') then call calc_dfdy(field, slp_field, sm1, sm2, min_level, em1, em2, max_level, & - which_domain, sr_y=sub_y) + which_domain, subgrid_var=subgrid_var) end if call write_field(sm1, em1, sm2, em2, & min_level, max_level, trim(gradname), datestr, real_array=slp_field) @@ -940,14 +936,14 @@ subroutine process_tile(which_domain, grid_type, dynopt, & i1=field_count,i2=NUM_FIELDS-NUM_AUTOMATIC_FIELDS,s1=domname) end if - if ((sub_x > 1) .or. (sub_y > 1)) then + if ( subgrid_var ) then call calc_field(fieldname, field, xlat_ptr, xlon_ptr, istagger, & sm1, em1, sm2, em2, min_category, max_category, & - processed_domain, 1, sr_x=sub_x, sr_y=sub_y) + processed_domain, 1, sr_x=sr_x, sr_y=sr_y) else call calc_field(fieldname, field, xlat_ptr, xlon_ptr, istagger, & sm1, em1, sm2, em2, min_category, max_category, & - processed_domain, 1, landmask=landmask, sr_x=sub_x, sr_y=sub_y) + processed_domain, 1, landmask=landmask) end if ! Find fractions @@ -1944,7 +1940,7 @@ end subroutine process_neighbor ! the result in dst_array. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine calc_dfdy(src_arr, dst_arr, start_mem_i, start_mem_j, start_mem_k, & - end_mem_i, end_mem_j, end_mem_k, idom, mapfac, sr_y) + end_mem_i, end_mem_j, end_mem_k, idom, mapfac, subgrid_var) ! Modules use gridinfo_module @@ -1957,7 +1953,7 @@ subroutine calc_dfdy(src_arr, dst_arr, start_mem_i, start_mem_j, start_mem_k, & real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j,start_mem_k:end_mem_k), intent(out) :: dst_arr integer,intent(in) :: idom real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(in), optional :: mapfac - integer, intent(in), optional :: sr_y + logical, intent(in), optional :: subgrid_var ! Local variables integer :: i, j, k @@ -1965,7 +1961,9 @@ subroutine calc_dfdy(src_arr, dst_arr, start_mem_i, start_mem_j, start_mem_k, & ! calculate the refinement ratio from the top-level domain l_sr_y=1. - if (present(sr_y)) l_sr_y=sr_y + if (present(subgrid_var)) then + if (subgrid_var) l_sr_y=subgrid_ratio_y(idom) + end if l_sr_y=l_sr_y*moad_grid_ratio(idom) if (present(mapfac)) then @@ -2012,7 +2010,7 @@ end subroutine calc_dfdy ! the result in dst_array. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine calc_dfdx(src_arr, dst_arr, start_mem_i, start_mem_j, & - start_mem_k, end_mem_i, end_mem_j, end_mem_k, idom, mapfac, sr_x) + start_mem_k, end_mem_i, end_mem_j, end_mem_k, idom, mapfac, subgrid_var) ! Modules use gridinfo_module @@ -2025,7 +2023,7 @@ subroutine calc_dfdx(src_arr, dst_arr, start_mem_i, start_mem_j, & real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j, start_mem_k:end_mem_k), intent(out) :: dst_arr integer, intent(in) :: idom real, dimension(start_mem_i:end_mem_i, start_mem_j:end_mem_j), intent(in), optional :: mapfac - integer, intent(in), optional :: sr_x + logical, intent(in), optional :: subgrid_var ! Local variables integer :: i, j, k @@ -2033,7 +2031,9 @@ subroutine calc_dfdx(src_arr, dst_arr, start_mem_i, start_mem_j, & ! calculate the refinement ratio from the top-level domain l_sr_x=1. - if (present(sr_x)) l_sr_x=sr_x + if (present(subgrid_var)) then + if (subgrid_var) l_sr_x=subgrid_ratio_x(idom) + end if l_sr_x=l_sr_x*moad_grid_ratio(idom) if (present(mapfac)) then diff --git a/WPS/geogrid/src/source_data_module.F b/WPS/geogrid/src/source_data_module.F index eeea0f9f..6041e145 100644 --- a/WPS/geogrid/src/source_data_module.F +++ b/WPS/geogrid/src/source_data_module.F @@ -1137,7 +1137,7 @@ end subroutine get_next_fieldname recursive subroutine get_next_output_fieldname(nest_num, field_name, ndims, & min_cat, max_cat, & istagger, memorder, dimnames, units, & - description, sr_x, sr_y, istatus) + description, subgrid_var, istatus) use gridinfo_module @@ -1148,7 +1148,7 @@ recursive subroutine get_next_output_fieldname(nest_num, field_name, ndims, & ! Arguments integer, intent(in) :: nest_num integer, intent(out) :: istatus, ndims, istagger, min_cat, max_cat - integer, intent(out) :: sr_x, sr_y + logical, intent(out) :: subgrid_var character (len=128), intent(out) :: memorder, field_name, units, description character (len=128), dimension(3), intent(out) :: dimnames @@ -1175,7 +1175,7 @@ recursive subroutine get_next_output_fieldname(nest_num, field_name, ndims, & call get_next_output_fieldname(nest_num, field_name, ndims, & min_cat, max_cat, istagger, & memorder, dimnames, units, description, & - sr_x, sr_y, istatus) + subgrid_var, istatus) return else ndims = 2 @@ -1218,7 +1218,10 @@ recursive subroutine get_next_output_fieldname(nest_num, field_name, ndims, & memorder = ' ' dimnames(3) = ' ' end if - call get_subgrid_dim_name(nest_num, field_name, dimnames(1:2), sr_x, sr_y, istatus) + subgrid_var=is_subgrid_var(field_name) + if ( subgrid_var ) then + call get_subgrid_dim_name(dimnames(1:2)) + end if call get_source_units(field_name, 1, units, istatus) if (istatus /= 0) units = '-' call get_source_descr(field_name, 1, description, istatus) @@ -1231,7 +1234,7 @@ recursive subroutine get_next_output_fieldname(nest_num, field_name, ndims, & call get_next_output_fieldname(nest_num, field_name, ndims, & min_cat, max_cat, istagger, & memorder, dimnames, units, description, & - sr_x, sr_y, istatus) + subgrid_var, istatus) return end if @@ -1248,7 +1251,7 @@ recursive subroutine get_next_output_fieldname(nest_num, field_name, ndims, & call get_next_output_fieldname(nest_num, field_name, ndims, & min_cat, max_cat, istagger, & memorder, dimnames, units, description, & - sr_x, sr_y, istatus) + subgrid_var, istatus) return ! Return the fractional field @@ -1293,7 +1296,12 @@ recursive subroutine get_next_output_fieldname(nest_num, field_name, ndims, & memorder = ' ' dimnames(3) = ' ' end if - call get_subgrid_dim_name(nest_num, field_name, dimnames(1:2), sr_x, sr_y, istatus) + + subgrid_var=is_subgrid_var(field_name) + if ( subgrid_var ) then + call get_subgrid_dim_name(dimnames(1:2)) + end if + call get_source_units(field_name, 1, units, istatus) if (istatus /= 0) units = '-' call get_source_descr(field_name, 1, description, istatus) @@ -1309,7 +1317,7 @@ recursive subroutine get_next_output_fieldname(nest_num, field_name, ndims, & call get_next_output_fieldname(nest_num, field_name, ndims, & min_cat, max_cat, istagger, & memorder, dimnames, units, description, & - sr_x, sr_y, istatus) + subgrid_var, istatus) return end if @@ -1347,7 +1355,11 @@ recursive subroutine get_next_output_fieldname(nest_num, field_name, ndims, & dimnames(3) = ' ' memorder = 'XY ' - call get_subgrid_dim_name(nest_num, field_name, dimnames(1:2), sr_x, sr_y, istatus) + subgrid_var=is_subgrid_var(field_name) + if ( subgrid_var ) then + call get_subgrid_dim_name(dimnames(1:2)) + end if + field_name = domcat_name units = 'category' description = 'Dominant category' @@ -1365,7 +1377,7 @@ recursive subroutine get_next_output_fieldname(nest_num, field_name, ndims, & call get_next_output_fieldname(nest_num, field_name, ndims, & min_cat, max_cat, istagger, & memorder, dimnames, units, description, & - sr_x, sr_y, istatus) + subgrid_var, istatus) return end if else @@ -1426,7 +1438,11 @@ recursive subroutine get_next_output_fieldname(nest_num, field_name, ndims, & end if units = '-' - call get_subgrid_dim_name(nest_num, field_name, dimnames(1:2), sr_x, sr_y, istatus) + subgrid_var=is_subgrid_var(field_name) + if ( subgrid_var ) then + call get_subgrid_dim_name(dimnames(1:2)) + end if + field_name = dfdx_name description = 'df/dx' if (output_field_state == RETURN_DFDX) then @@ -1443,7 +1459,7 @@ recursive subroutine get_next_output_fieldname(nest_num, field_name, ndims, & call get_next_output_fieldname(nest_num, field_name, ndims, & min_cat, max_cat, istagger, & memorder, dimnames, units, description, & - sr_x, sr_y, istatus) + subgrid_var, istatus) return end if else @@ -1503,7 +1519,11 @@ recursive subroutine get_next_output_fieldname(nest_num, field_name, ndims, & dimnames(3) = ' ' end if - call get_subgrid_dim_name(nest_num, field_name, dimnames(1:2), sr_x, sr_y, istatus) + subgrid_var=is_subgrid_var(field_name) + if ( subgrid_var ) then + call get_subgrid_dim_name(dimnames(1:2)) + end if + field_name = dfdy_name units = '-' description = 'df/dy' @@ -1513,7 +1533,7 @@ recursive subroutine get_next_output_fieldname(nest_num, field_name, ndims, & call get_next_output_fieldname(nest_num, field_name, ndims, & min_cat, max_cat, istagger, & memorder, dimnames, units, description, & - sr_x, sr_y, istatus) + subgrid_var, istatus) return end if else @@ -2319,61 +2339,47 @@ end subroutine get_output_stagger ! ! Pupose: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine get_subgrid_dim_name(nest_num, field_name, dimnames, & - sub_x, sub_y, istatus) + logical function is_subgrid_var(field_name) use gridinfo_module implicit none - integer, intent(in) :: nest_num - integer, intent(out) :: sub_x, sub_y, istatus character(len=128), intent(in) :: field_name - character(len=128), dimension(2), intent(inout) :: dimnames integer :: idx, nlen - sub_x = 1 - sub_y = 1 + is_subgrid_var=.false. - istatus = 0 if (((index(trim(field_name),'FXLAT') /= 0) .and. & (len_trim(field_name) == len_trim('FXLAT'))) .or. & ((index(trim(field_name),'FXLONG') /= 0) .and. & (len_trim(field_name) == len_trim('FXLONG'))))then - call get_subgrid_dim_default(nest_num,dimnames,sub_x,sub_y,istatus) - endif + is_subgrid_var=.true. + else - do idx=1,num_entries - if ((index(source_fieldname(idx),trim(field_name)) /= 0) .and. & - (len_trim(source_fieldname(idx)) == len_trim(field_name))) then - if (is_subgrid(idx)) then - istatus = 0 - if (is_output_stagger(idx)) then - call mprintf(.true.,ERROR,'Cannot use subgrids on variables with staggered grids') - end if - call get_subgrid_dim_default(nest_num,dimnames,sub_x,sub_y,istatus) + do idx=1,num_entries + if ((index(source_fieldname(idx),trim(field_name)) /= 0) .and. & + (len_trim(source_fieldname(idx)) == len_trim(field_name))) then + if (is_subgrid(idx)) is_subgrid_var=.true. + goto 110 end if - end if - end do + end do - end subroutine get_subgrid_dim_name + end if + 110 continue - subroutine get_subgrid_dim_default(nest_num,dimnames,sub_x,sub_y,istatus) - use gridinfo_module + end function is_subgrid_var + + subroutine get_subgrid_dim_name(dimnames) implicit none - integer,intent(in)::nest_num - integer,intent(out)::sub_x,sub_y,istatus character(len=128),dimension(2),intent(out)::dimnames dimnames(1)='west_east_subgrid' dimnames(2)='south_north_subgrid' - sub_x=subgrid_ratio_x(nest_num) - sub_y=subgrid_ratio_y(nest_num) - istatus=0 - end subroutine get_subgrid_dim_default + end subroutine get_subgrid_dim_name !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! diff --git a/WPS/metgrid/src/datatype_module.F b/WPS/metgrid/src/datatype_module.F index 066ba9b0..ed7f3eaf 100644 --- a/WPS/metgrid/src/datatype_module.F +++ b/WPS/metgrid/src/datatype_module.F @@ -37,6 +37,7 @@ module datatype_module logical :: is_wind_grid_rel logical :: array_has_missing_values + logical :: is_subgrid integer :: sr_x, sr_y end type header_info diff --git a/WPS/metgrid/src/input_module.F b/WPS/metgrid/src/input_module.F index b4ee3e24..6e13f31e 100644 --- a/WPS/metgrid/src/input_module.F +++ b/WPS/metgrid/src/input_module.F @@ -114,18 +114,18 @@ subroutine read_next_field(start_patch_i, end_patch_i, & start_patch_j, end_patch_j, & start_patch_k, end_patch_k, & cname, cunits, cdesc, memorder, stagger, & - dimnames, sr_x, sr_y, real_array, istatus) + dimnames, is_subgrid, real_array, istatus) implicit none ! Arguments integer, intent(out) :: start_patch_i, end_patch_i, & start_patch_j, end_patch_j, & - start_patch_k, end_patch_k, & - sr_x, sr_y + start_patch_k, end_patch_k real, pointer, dimension(:,:,:) :: real_array character (len=*), intent(out) :: cname, memorder, stagger, cunits, cdesc character (len=128), dimension(3), intent(inout) :: dimnames + logical, intent(out) :: is_subgrid integer, intent(inout) :: istatus #include "wrf_io_flags.h" @@ -138,6 +138,7 @@ subroutine read_next_field(start_patch_i, end_patch_i, & real, pointer, dimension(:,:,:) :: real_domain character (len=20) :: datestr type (q_data) :: qd + integer :: sr_x, sr_y if (my_proc_id == IO_NODE .or. do_tiled_input) then @@ -175,32 +176,33 @@ subroutine read_next_field(start_patch_i, end_patch_i, & #ifdef IO_BINARY if (io_form_input == BINARY) then call ext_int_get_var_info(handle, cname, ndim, memorder, stagger, domain_start, domain_end, wrftype, istatus) - call ext_int_get_var_ti_integer(handle, 'sr_x', & + call ext_int_get_var_ti_integer(handle, 'subgrid', & trim(cname), temp(1), 1, temp(3), istatus) - call ext_int_get_var_ti_integer(handle, 'sr_y', & - trim(cname), temp(2), 1, temp(3), istatus) end if #endif #ifdef IO_NETCDF if (io_form_input == NETCDF) then call ext_ncd_get_var_info(handle, cname, ndim, memorder, stagger, domain_start, domain_end, wrftype, istatus) - call ext_ncd_get_var_ti_integer(handle, 'sr_x', & + call ext_ncd_get_var_ti_integer(handle, 'subgrid', & trim(cname), temp(1), 1, temp(3), istatus) - call ext_ncd_get_var_ti_integer(handle, 'sr_y', & - trim(cname), temp(2), 1, temp(3), istatus) end if #endif #ifdef IO_GRIB1 if (io_form_input == GRIB1) then call ext_gr1_get_var_info(handle, cname, ndim, memorder, stagger, domain_start, domain_end, wrftype, istatus) - call ext_gr1_get_var_ti_integer(handle, 'sr_x', & + call ext_gr1_get_var_ti_integer(handle, 'subgrid', & trim(cname), temp(1), 1, temp(3), istatus) - call ext_gr1_get_var_ti_integer(handle, 'sr_y', & - trim(cname), temp(2), 1, temp(3), istatus) end if #endif - sr_x = temp(1) - sr_y = temp(2) + + call ext_get_dom_ti_integer_scalar('sr_x', sr_x) + call ext_get_dom_ti_integer_scalar('sr_y', sr_y) + + if (temp(1) .eq. 1) then + is_subgrid=.true. + else + is_subgrid=.false. + end if call mprintf((istatus /= 0),ERROR,'In read_next_field(), problems with ext_pkg_get_var_info()') @@ -287,8 +289,7 @@ subroutine read_next_field(start_patch_i, end_patch_i, & call parallel_bcast_char(dimnames(3), 128) call parallel_bcast_int(domain_start(3)) call parallel_bcast_int(domain_end(3)) - call parallel_bcast_int(sr_x) - call parallel_bcast_int(sr_y) + call parallel_bcast_logical(is_subgrid) sp1 = my_minx ep1 = my_maxx - 1 @@ -298,10 +299,10 @@ subroutine read_next_field(start_patch_i, end_patch_i, & ep3 = domain_end(3) if (internal_gridtype == 'C') then - if (my_x /= nproc_x - 1 .or. stagger == 'U' .or. sr_x > 1) then + if (my_x /= nproc_x - 1 .or. stagger == 'U' .or. is_subgrid) then ep1 = ep1 + 1 end if - if (my_y /= nproc_y - 1 .or. stagger == 'V' .or. sr_y > 1) then + if (my_y /= nproc_y - 1 .or. stagger == 'V' .or. is_subgrid) then ep2 = ep2 + 1 end if else if (internal_gridtype == 'E') then diff --git a/WPS/metgrid/src/process_domain_module.F b/WPS/metgrid/src/process_domain_module.F index 79c1461e..8462e6cb 100644 --- a/WPS/metgrid/src/process_domain_module.F +++ b/WPS/metgrid/src/process_domain_module.F @@ -217,7 +217,7 @@ subroutine get_static_fields(n, dyn_opt, west_east_dim, south_north_dim, bottom_ ! Local variables integer :: istatus, i, j, k, sp1, ep1, sp2, ep2, sp3, ep3, & - lh_mult, rh_mult, bh_mult, th_mult, subx, suby + 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, & @@ -227,6 +227,7 @@ subroutine get_static_fields(n, dyn_opt, west_east_dim, south_north_dim, bottom_ 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.') @@ -338,7 +339,7 @@ subroutine get_static_fields(n, dyn_opt, west_east_dim, south_north_dim, bottom_ istatus = 0 do while (istatus == 0) call read_next_field(sp1, ep1, sp2, ep2, sp3, ep3, cname, cunits, cdesc, & - memorder, stagger, dimnames, subx, suby, & + memorder, stagger, dimnames, is_subgrid, & real_array, istatus) if (istatus == 0) then @@ -405,25 +406,22 @@ subroutine get_static_fields(n, dyn_opt, west_east_dim, south_north_dim, bottom_ 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 + 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=subx - field%header%sr_y=suby + 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. @@ -439,7 +437,7 @@ subroutine get_static_fields(n, dyn_opt, west_east_dim, south_north_dim, bottom_ 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 + if (is_subgrid) then field%map%stagger = M field%header%dim1(1) = we_mem_subgrid_s field%header%dim1(2) = we_mem_subgrid_e @@ -478,7 +476,7 @@ subroutine get_static_fields(n, dyn_opt, west_east_dim, south_north_dim, bottom_ allocate(field%valid_mask) - if (subx > 1 .or. suby > 1) then + 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) = & @@ -767,8 +765,9 @@ subroutine process_single_met_time(do_const_processing, & 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%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. @@ -1318,6 +1317,7 @@ subroutine get_interp_masks(fg_prefix, is_constants, fg_date) 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 @@ -2044,6 +2044,7 @@ subroutine create_derived_fields(arg_gridtype, hdate, xfcst, & 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. diff --git a/WPS/metgrid/src/storage_module.F b/WPS/metgrid/src/storage_module.F index a1159e78..7ca03fba 100644 --- a/WPS/metgrid/src/storage_module.F +++ b/WPS/metgrid/src/storage_module.F @@ -339,7 +339,7 @@ end subroutine storage_query_field subroutine get_next_output_fieldname(nest_num, field_name, ndims, & min_level, max_level, & istagger, mem_order, dim_names, units, description, & - sr_x, sr_y, & + subgrid_var, & istatus) implicit none @@ -347,7 +347,7 @@ subroutine get_next_output_fieldname(nest_num, field_name, ndims, & ! Arguments integer, intent(in) :: nest_num integer, intent(out) :: ndims, min_level, max_level, istagger, istatus - integer, intent(out) :: sr_x, sr_y + logical, intent(out) :: subgrid_var character (len=128), intent(out) :: field_name, mem_order, units, description character (len=128), dimension(3), intent(out) :: dim_names @@ -412,9 +412,9 @@ subroutine get_next_output_fieldname(nest_num, field_name, ndims, & dim_names(2) = 'j-dimension' end if units = get_units(next_output_field%fg_data) - description = get_description(next_output_field%fg_data) - call get_subgrid_dim_name(nest_num, field_name, dim_names(1:2), & - sr_x, sr_y, istatus) + description = get_description(next_output_field%fg_data) + subgrid_var=is_subgrid_var(field_name) + if (subgrid_var) call get_subgrid_dim_name(dim_names(1:2)) next_output_field => next_output_field%next end if @@ -427,36 +427,28 @@ end subroutine get_next_output_fieldname ! ! Purpose: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine get_subgrid_dim_name(nest_num, field_name, dimnames, & - sub_x, sub_y, istatus) + logical function is_subgrid_var(field_name) use gridinfo_module implicit none ! Arguments - integer, intent(in) :: nest_num - integer, intent(out) :: sub_x, sub_y, istatus character(len=128), intent(in) :: field_name - character(len=128), dimension(2), intent(inout) :: dimnames - ! Local variables - integer :: idx, nlen + is_subgrid_var=next_output_field%fg_data%header%is_subgrid - sub_x = next_output_field%fg_data%header%sr_x - sub_y = next_output_field%fg_data%header%sr_y + end function is_subgrid_var - if (sub_x > 1) then - dimnames(1) = trim(dimnames(1))//"_subgrid" - end if - if (sub_y > 1) then - dimnames(2) = trim(dimnames(2))//"_subgrid" - end if + subroutine get_subgrid_dim_name(dimnames) + implicit none - istatus = 0 + character(len=128),dimension(2),intent(out)::dimnames - end subroutine get_subgrid_dim_name + dimnames(1)='west_east_subgrid' + dimnames(2)='south_north_subgrid' + end subroutine get_subgrid_dim_name !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Name: get_next_output_field