diff --git a/src/drivers/mct/cime_config/namelist_definition_drv.xml b/src/drivers/mct/cime_config/namelist_definition_drv.xml index 82a7d634e78..cd9aff7d39c 100644 --- a/src/drivers/mct/cime_config/namelist_definition_drv.xml +++ b/src/drivers/mct/cime_config/namelist_definition_drv.xml @@ -572,6 +572,39 @@ + + char + control + seq_infodata_inparm + on,off,on_if_glc_coupled_fluxes + + Whether to renormalize the surface mass balance (smb) sent from lnd to glc so that the + global integral on the glc grid agrees with the global integral on the lnd grid. + + Unlike most fluxes, smb is remapped with bilinear rather than conservative mapping weights, + so this option is needed for conservation. However, conservation is not required in many + cases, since we often run glc as a diagnostic (one-way-coupled) component. + + Allowable values are: + 'on': always do this renormalization + 'off': never do this renormalization (see WARNING below) + 'on_if_glc_coupled_fluxes': Determine at runtime whether to do this renormalization. + Does the renormalization if we're running a two-way-coupled glc that sends fluxes + to other components (which is the case where we need conservation). + Does NOT do the renormalization if we're running a one-way-coupled glc, or if + we're running a glc-only compset (T compsets). + (In these cases, conservation is not important.) + + Only used if running with a prognostic GLC component. + + WARNING: Setting this to 'off' will break conservation when running with an + evolving, two-way-coupled glc. + + + on_if_glc_coupled_fluxes + + + real control diff --git a/src/drivers/mct/main/CMakeLists.txt b/src/drivers/mct/main/CMakeLists.txt index 537d5834b91..53cabe62ba2 100644 --- a/src/drivers/mct/main/CMakeLists.txt +++ b/src/drivers/mct/main/CMakeLists.txt @@ -1,13 +1,9 @@ list(APPEND drv_sources component_type_mod.F90 map_glc2lnd_mod.F90 - map_lnd2glc_mod.F90 map_lnd2rof_irrig_mod.F90 seq_map_mod.F90 seq_map_type_mod.F90 - vertical_gradient_calculator_base.F90 - vertical_gradient_calculator_2nd_order.F90 - vertical_gradient_calculator_factory.F90 ) sourcelist_to_parent(drv_sources) \ No newline at end of file diff --git a/src/drivers/mct/main/map_lnd2glc_mod.F90 b/src/drivers/mct/main/map_lnd2glc_mod.F90 index da2ea515880..94c2bb899a7 100644 --- a/src/drivers/mct/main/map_lnd2glc_mod.F90 +++ b/src/drivers/mct/main/map_lnd2glc_mod.F90 @@ -8,19 +8,19 @@ module map_lnd2glc_mod ! elevation class) onto the GLC grid ! ! For high-level design, see: - ! https://docs.google.com/document/d/1sjsaiPYsPJ9A7dVGJIHGg4rVIY2qF5aRXbNzSXVAafU/edit?usp=sharing + ! https://docs.google.com/document/d/1H_SuK6SfCv1x6dK91q80dFInPbLYcOkUj_iAa6WRnqQ/edit #include "shr_assert.h" - use seq_comm_mct, only : logunit + use seq_comm_mct, only: CPLID, GLCID, logunit use shr_kind_mod, only : r8 => shr_kind_r8 + use shr_kind_mod, only : cxx => SHR_KIND_CXX use glc_elevclass_mod, only : glc_get_num_elevation_classes, glc_get_elevation_class, & - glc_elevclass_as_string, GLC_ELEVCLASS_ERR_NONE, GLC_ELEVCLASS_ERR_TOO_LOW, & + glc_elevclass_as_string, glc_all_elevclass_strings, GLC_ELEVCLASS_STRLEN, & + GLC_ELEVCLASS_ERR_NONE, GLC_ELEVCLASS_ERR_TOO_LOW, & GLC_ELEVCLASS_ERR_TOO_HIGH, glc_errcode_to_string use mct_mod use seq_map_type_mod, only : seq_map use seq_map_mod, only : seq_map_map - use vertical_gradient_calculator_base, only : vertical_gradient_calculator_base_type - use shr_log_mod, only : errMsg => shr_log_errMsg use shr_sys_mod, only : shr_sys_abort implicit none @@ -39,12 +39,12 @@ module map_lnd2glc_mod private :: get_glc_elevation_classes ! get the elevation class of each point on the glc grid private :: map_bare_land ! remap the field of interest for the bare land "elevation class" - private :: map_one_elevation_class ! remap the field of interest for one ice elevation class + private :: map_ice_covered ! remap the field of interest for all elevation classes (excluding bare land) contains !----------------------------------------------------------------------- - subroutine map_lnd2glc(l2x_l, landfrac_l, g2x_g, fieldname, gradient_calculator, & + subroutine map_lnd2glc(l2x_l, landfrac_l, g2x_g, fieldname, & mapper, l2x_g) ! ! !DESCRIPTION: @@ -83,7 +83,6 @@ subroutine map_lnd2glc(l2x_l, landfrac_l, g2x_g, fieldname, gradient_calculator, type(mct_aVect) , intent(in) :: landfrac_l ! lfrac field on the land grid type(mct_aVect) , intent(in) :: g2x_g ! glc -> cpl fields on the glc grid character(len=*) , intent(in) :: fieldname ! name of the field to map - class(vertical_gradient_calculator_base_type), intent(inout) :: gradient_calculator type(seq_map) , intent(inout) :: mapper type(mct_aVect) , intent(inout) :: l2x_g ! lnd -> cpl fields on the glc grid ! @@ -92,15 +91,16 @@ subroutine map_lnd2glc(l2x_l, landfrac_l, g2x_g, fieldname, gradient_calculator, ! fieldname with trailing blanks removed character(len=:), allocatable :: fieldname_trimmed - ! index for looping over elevation classes - integer :: elevclass - ! number of points on the GLC grid integer :: lsize_g - ! data for one elevation class on the GLC grid + ! data for bare land on the GLC grid + ! needs to be a pointer to satisfy the MCT interface + real(r8), pointer :: data_g_bareland(:) + + ! data for ice-covered regions on the GLC grid ! needs to be a pointer to satisfy the MCT interface - real(r8), pointer :: data_g_oneEC(:) + real(r8), pointer :: data_g_ice_covered(:) ! final data on the GLC grid ! needs to be a pointer to satisfy the MCT interface @@ -126,7 +126,9 @@ subroutine map_lnd2glc(l2x_l, landfrac_l, g2x_g, fieldname, gradient_calculator, ! ------------------------------------------------------------------------ lsize_g = mct_aVect_lsize(l2x_g) - allocate(data_g_oneEC(lsize_g)) + + allocate(data_g_ice_covered(lsize_g)) + allocate(data_g_bareland(lsize_g)) allocate(data_g(lsize_g)) fieldname_trimmed = trim(fieldname) @@ -148,10 +150,10 @@ subroutine map_lnd2glc(l2x_l, landfrac_l, g2x_g, fieldname, gradient_calculator, call get_glc_elevation_classes(glc_ice_covered, glc_topo, glc_elevclass) ! ------------------------------------------------------------------------ - ! Map elevation class 0 (bare land) + ! Map elevation class 0 (bare land) and ice elevation classes ! ------------------------------------------------------------------------ - call map_bare_land(l2x_l, landfrac_l, fieldname_trimmed, mapper, data_g_oneEC) + call map_bare_land(l2x_l, landfrac_l, fieldname_trimmed, mapper, data_g_bareland) ! Start by setting the output data equal to the bare land value everywhere; this will ! later get overwritten in places where we have ice @@ -161,21 +163,15 @@ subroutine map_lnd2glc(l2x_l, landfrac_l, g2x_g, fieldname, gradient_calculator, ! current glint implementation, which sets acab and artm to 0 over ocean (although ! notes that this could lead to a loss of conservation). Figure out how to handle ! this case. - data_g(:) = data_g_oneEC(:) + data_g(:) = data_g_bareland(:) - ! ------------------------------------------------------------------------ - ! Map ice elevation classes - ! ------------------------------------------------------------------------ + ! Map the SMB to ice-covered cells + call map_ice_covered(l2x_l, landfrac_l, fieldname_trimmed, & + glc_topo, mapper, data_g_ice_covered) - call gradient_calculator%calc_gradients() - do elevclass = 1, glc_get_num_elevation_classes() - call map_one_elevation_class(l2x_l, landfrac_l, fieldname_trimmed, elevclass, & - gradient_calculator, glc_topo, mapper, data_g_oneEC) - - where (glc_elevclass == elevclass) - data_g = data_g_oneEC - end where - end do + where (glc_elevclass /= 0) + data_g = data_g_ice_covered + end where ! ------------------------------------------------------------------------ ! Set field in output attribute vector @@ -187,7 +183,8 @@ subroutine map_lnd2glc(l2x_l, landfrac_l, g2x_g, fieldname, gradient_calculator, ! Clean up ! ------------------------------------------------------------------------ - deallocate(data_g_oneEC) + deallocate(data_g_ice_covered) + deallocate(data_g_bareland) deallocate(data_g) deallocate(glc_ice_covered) deallocate(glc_topo) @@ -305,147 +302,180 @@ subroutine map_bare_land(l2x_l, landfrac_l, fieldname, mapper, data_g_bare_land) end subroutine map_bare_land - !----------------------------------------------------------------------- - subroutine map_one_elevation_class(l2x_l, landfrac_l, fieldname, elevclass, & - gradient_calculator, topo_g, mapper, data_g_thisEC) + subroutine map_ice_covered(l2x_l, landfrac_l, fieldname, & + topo_g, mapper, data_g_ice_covered) + ! ! !DESCRIPTION: - ! Remaps the field of interest for a single ice elevation class. + ! Remaps the field of interest from the land grid (in multiple elevation classes) + ! to the glc grid ! - ! Puts the output in data_g_thisEC, which should already be allocated to have size + ! Puts the output in data_g_ice_covered, which should already be allocated to have size ! equal to the number of GLC points that this processor is responsible for. - ! - ! To do this remapping, we remap the field adjusted by the vertical gradient. That is, - ! rather than mapping data_l itself, we instead remap: - ! - ! data_l + (vertical_gradient_l) * (topo_g - topo_l) - ! - ! (where _l denotes quantities on the land grid, _g on the glc grid) - ! - ! However, in order to do the remapping with existing routines, we do this by - ! performing two separate remappings: - ! - ! (1) Remap (data_l - vertical_gradient_l * topo_l); put result in partial_remap_g - ! (note: in variables in the code, the parenthesized term is called partial_remap, - ! either on the land or glc grid) - ! - ! (2) Remap vertical_gradient_l; put result in vertical_gradient_g - ! - ! Then data_g = partial_remap_g + topo_g * vertical_gradient_g - ! + ! ! !USES: ! ! !ARGUMENTS: type(mct_aVect) , intent(in) :: l2x_l ! lnd -> cpl fields on the land grid type(mct_aVect) , intent(in) :: landfrac_l ! lfrac field on the land grid character(len=*) , intent(in) :: fieldname ! name of the field to map (should have NO trailing blanks) - integer , intent(in) :: elevclass ! elevation class index to map - class(vertical_gradient_calculator_base_type), intent(in) :: gradient_calculator real(r8) , intent(in) :: topo_g(:) ! topographic height for each point on the glc grid type(seq_map) , intent(inout) :: mapper - real(r8) , intent(out) :: data_g_thisEC(:) - ! + real(r8) , intent(out) :: data_g_ice_covered(:) ! field remapped to glc grid + ! !LOCAL VARIABLES: - ! Fields contained in the temporary attribute vectors: - character(len=*), parameter :: partial_remap_tag = 'partial_remap' - character(len=*), parameter :: vertical_gradient_tag = 'vertical_gradient' - character(len=*), parameter :: attr_tags = & - partial_remap_tag // ':' // vertical_gradient_tag - - ! Base name for the topo fields in l2x_l. Actual fields will have an elevation class suffix. - character(len=*), parameter :: toponame = 'Sl_topo' + character(len=*), parameter :: toponame = 'Sl_topo' ! base name for topo fields in l2x_l; + ! actual names will have elevation class suffix + character(len=GLC_ELEVCLASS_STRLEN), allocatable :: all_elevclass_strings(:) character(len=:), allocatable :: elevclass_as_string character(len=:), allocatable :: fieldname_ec character(len=:), allocatable :: toponame_ec - integer :: lsize_l ! number of points for attribute vectors on the land grid - integer :: lsize_g ! number of points for attribute vectors on the glc grid - type(mct_aVect) :: l2x_l_temp + character(len=:), allocatable :: fieldnamelist + character(len=:), allocatable :: toponamelist + character(len=:), allocatable :: totalfieldlist + + integer :: nEC ! number of elevation classes + integer :: lsize_g ! number of cells on glc grid + integer :: n, ec + integer :: strlen + + real(r8) :: elev_l, elev_u ! lower and upper elevations in interpolation range + real(r8) :: d_elev ! elev_u - elev_l + type(mct_aVect) :: l2x_g_temp ! temporary attribute vector holding the remapped fields for this elevation class - ! Note that arrays passed to MCT routines need to be pointers - ! Temporary fields on the land grid: - real(r8), pointer :: data_l(:) - real(r8), pointer :: topo_l(:) - real(r8), pointer :: vertical_gradient_l(:) - real(r8), pointer :: partial_remap_l(:) - ! Temporary fields on the glc grid: - real(r8), pointer :: vertical_gradient_g(:) - real(r8), pointer :: partial_remap_g(:) - - character(len=*), parameter :: subname = 'map_one_elevation_class' + real(r8), pointer :: tmp_field_g(:) ! must be a pointer to satisfy the MCT interface + real, pointer :: data_g_EC(:,:) ! remapped field in each glc cell, in each EC + real, pointer :: topo_g_EC(:,:) ! remapped topo in each glc cell, in each EC + + ! 1 is probably enough, but use 10 to be safe, in case the length of the delimiter + ! changes + integer, parameter :: extra_len_for_list_merge = 10 + + character(len=*), parameter :: subname = 'map_ice_covered' !----------------------------------------------------------------------- - lsize_g = size(data_g_thisEC) + lsize_g = size(data_g_ice_covered) + nEC = glc_get_num_elevation_classes() SHR_ASSERT_FL((size(topo_g) == lsize_g), __FILE__, __LINE__) - + ! ------------------------------------------------------------------------ - ! Create temporary attribute vectors + ! Create temporary vectors ! ------------------------------------------------------------------------ - - lsize_l = mct_aVect_lsize(l2x_l) - call mct_aVect_init(l2x_l_temp, rList = attr_tags, lsize = lsize_l) - call mct_aVect_init(l2x_g_temp, rList = attr_tags, lsize = lsize_g) - + + allocate(tmp_field_g(lsize_g)) + allocate(data_g_EC (lsize_g,nEC)) + allocate(topo_g_EC (lsize_g,nEC)) + ! ------------------------------------------------------------------------ - ! Create fields to remap on the source (land) grid + ! Make a string that concatenates all EC levels of field, as well as the topo + ! The resulting list will look something like this: + ! 'Flgl_qice01:Flgl_qice02: ... :Flgl_qice10:Sl_topo01:Sl_topo02: ... :Sltopo10' ! ------------------------------------------------------------------------ - allocate(data_l(lsize_l)) - allocate(topo_l(lsize_l)) - allocate(vertical_gradient_l(lsize_l)) - allocate(partial_remap_l(lsize_l)) - - elevclass_as_string = glc_elevclass_as_string(elevclass) - fieldname_ec = fieldname // elevclass_as_string - toponame_ec = toponame // elevclass_as_string - - call mct_aVect_exportRattr(l2x_l, fieldname_ec, data_l) - call mct_aVect_exportRattr(l2x_l, toponame_ec, topo_l) - call gradient_calculator%get_gradients_one_class(elevclass, vertical_gradient_l) - partial_remap_l = data_l - (vertical_gradient_l * topo_l) - - call mct_aVect_importRattr(l2x_l_temp, partial_remap_tag, partial_remap_l) - call mct_aVect_importRattr(l2x_l_temp, vertical_gradient_tag, vertical_gradient_l) - + allocate(all_elevclass_strings(1:glc_get_num_elevation_classes())) + all_elevclass_strings = glc_all_elevclass_strings(include_zero = .false.) + fieldnamelist = shr_string_listFromSuffixes( & + suffixes = all_elevclass_strings, & + strBase = fieldname) + toponamelist = shr_string_listFromSuffixes( & + suffixes = all_elevclass_strings, & + strBase = toponame) + strlen = len_trim(fieldnamelist) + len_trim(toponamelist) + extra_len_for_list_merge + allocate(character(len=strlen) :: totalfieldlist) + call shr_string_listMerge(fieldnamelist, toponamelist, totalfieldlist ) + ! ------------------------------------------------------------------------ - ! Remap to destination (glc) grid + ! Make a temporary attribute vector. + ! For each grid cell on the land grid, this attribute vector contains the field and + ! topo values for all ECs. + ! ------------------------------------------------------------------------ + call mct_aVect_init(l2x_g_temp, rList = totalfieldlist, lsize = lsize_g) + + ! ------------------------------------------------------------------------ + ! Remap all these fields from the land (source) grid to the glc (destination) grid. ! ------------------------------------------------------------------------ - call seq_map_map(mapper = mapper, av_s = l2x_l_temp, av_d = l2x_g_temp, & - norm = .true., & - avwts_s = landfrac_l, & - avwtsfld_s = 'lfrac') - + call seq_map_map(mapper = mapper, & + av_s = l2x_l, & + av_d = l2x_g_temp, & + fldlist = totalfieldlist, & + norm = .true., & + avwts_s = landfrac_l, & + avwtsfld_s = 'lfrac') + ! ------------------------------------------------------------------------ - ! Compute final field on the destination (glc) grid + ! Export all elevation classes out of attribute vector and into local 2D arrays (xy,z) ! ------------------------------------------------------------------------ - allocate(partial_remap_g(lsize_g)) - allocate(vertical_gradient_g(lsize_g)) + do ec = 1, nEC + elevclass_as_string = glc_elevclass_as_string(ec) + fieldname_ec = fieldname // elevclass_as_string + toponame_ec = toponame // elevclass_as_string + call mct_aVect_exportRattr(l2x_g_temp, fieldname_ec, tmp_field_g) + data_g_EC(:,ec) = tmp_field_g + call mct_aVect_exportRattr(l2x_g_temp, toponame_ec, tmp_field_g) + topo_g_EC(:,ec) = tmp_field_g + enddo + + ! ------------------------------------------------------------------------ + ! Perform vertical interpolation of data onto ice sheet topography + ! ------------------------------------------------------------------------ + + data_g_ice_covered(:) = 0._r8 + + do n = 1, lsize_g - call mct_aVect_exportRattr(l2x_g_temp, partial_remap_tag, partial_remap_g) - call mct_aVect_exportRattr(l2x_g_temp, vertical_gradient_tag, vertical_gradient_g) + ! For each ice sheet point, find bounding EC values... + if (topo_g(n) < topo_g_EC(n,1)) then + ! lower than lowest mean EC elevation value + data_g_ice_covered(n) = data_g_EC(n,1) - data_g_thisEC = partial_remap_g + (topo_g * vertical_gradient_g) + else if (topo_g(n) >= topo_g_EC(n,nEC)) then + ! higher than highest mean EC elevation value + data_g_ice_covered(n) = data_g_EC(n,nEC) + else + ! do linear interpolation of data in the vertical + do ec = 2, nEC + if (topo_g(n) < topo_g_EC(n, ec)) then + elev_l = topo_g_EC(n, ec-1) + elev_u = topo_g_EC(n, ec) + d_elev = elev_u - elev_l + if (d_elev <= 0) then + ! This shouldn't happen, but handle it in case it does. In this case, + ! let's arbitrarily use the mean of the two elevation classes, rather + ! than the weighted mean. + write(logunit,*) subname//' WARNING: topo diff between elevation classes <= 0' + write(logunit,*) 'n, ec, elev_l, elev_u = ', n, ec, elev_l, elev_u + write(logunit,*) 'Simply using mean of the two elevation classes,' + write(logunit,*) 'rather than the weighted mean.' + data_g_ice_covered(n) = data_g_EC(n,ec-1) * 0.5_r8 & + + data_g_EC(n,ec) * 0.5_r8 + else + data_g_ice_covered(n) = data_g_EC(n,ec-1) * (elev_u - topo_g(n)) / d_elev & + + data_g_EC(n,ec) * (topo_g(n) - elev_l) / d_elev + end if + + exit + end if + end do + end if ! topo_g(n) + end do ! lsize_g + ! ------------------------------------------------------------------------ ! Clean up ! ------------------------------------------------------------------------ + deallocate(tmp_field_g) + deallocate(data_g_EC) + deallocate(topo_g_EC) + call mct_aVect_clean(l2x_g_temp) - call mct_aVect_clean(l2x_l_temp) - deallocate(data_l) - deallocate(topo_l) - deallocate(vertical_gradient_l) - deallocate(partial_remap_l) - deallocate(vertical_gradient_g) - deallocate(partial_remap_g) - - end subroutine map_one_elevation_class - - + + end subroutine map_ice_covered end module map_lnd2glc_mod diff --git a/src/drivers/mct/main/prep_glc_mod.F90 b/src/drivers/mct/main/prep_glc_mod.F90 index 52eeb02bfa1..e14f78a24a0 100644 --- a/src/drivers/mct/main/prep_glc_mod.F90 +++ b/src/drivers/mct/main/prep_glc_mod.F90 @@ -1,5 +1,6 @@ module prep_glc_mod +#include "shr_assert.h" use shr_kind_mod , only: r8 => SHR_KIND_R8 use shr_kind_mod , only: cl => SHR_KIND_CL use shr_sys_mod , only: shr_sys_abort, shr_sys_flush @@ -14,7 +15,10 @@ module prep_glc_mod use mct_mod use perf_mod use component_type_mod, only: component_get_x2c_cx, component_get_c2x_cx + use component_type_mod, only: component_get_dom_cx use component_type_mod, only: glc, lnd + use glc_elevclass_mod, only : glc_get_num_elevation_classes, glc_elevclass_as_string + use glc_elevclass_mod, only : glc_all_elevclass_strings, GLC_ELEVCLASS_STRLEN implicit none save @@ -44,8 +48,12 @@ module prep_glc_mod ! Private interfaces !-------------------------------------------------------------------------- + private :: prep_glc_do_renormalize_smb + private :: prep_glc_set_g2x_lx_fields private :: prep_glc_merge - private :: prep_glc_map_one_field_lnd2glc + private :: prep_glc_map_one_state_field_lnd2glc + private :: prep_glc_map_qice_conservative_lnd2glc + private :: prep_glc_renormalize_smb !-------------------------------------------------------------------------- ! Private data @@ -54,6 +62,7 @@ module prep_glc_mod ! mappers type(seq_map), pointer :: mapper_Sl2g type(seq_map), pointer :: mapper_Fl2g + type(seq_map), pointer :: mapper_Fg2l ! attribute vectors type(mct_aVect), pointer :: l2x_gx(:) ! Lnd export, glc grid, cpl pes - allocated in driver @@ -64,6 +73,24 @@ module prep_glc_mod ! other module variables integer :: mpicom_CPLID ! MPI cpl communicator + + ! Whether to renormalize the SMB for conservation. + ! Should be set to true for 2-way coupled runs with evolving ice sheets. + ! Does not need to be true for 1-way coupling. + logical :: smb_renormalize + + ! Name of flux field giving surface mass balance + character(len=*), parameter :: qice_fieldname = 'Flgl_qice' + + ! Names of some other fields + character(len=*), parameter :: Sg_frac_field = 'Sg_ice_covered' + character(len=*), parameter :: Sg_topo_field = 'Sg_topo' + character(len=*), parameter :: Sg_icemask_field = 'Sg_icemask' + + ! Fields needed in the g2x_lx attribute vector used as part of mapping qice from lnd to glc + character(len=:), allocatable :: g2x_lx_fields + + !================================================================================================ contains @@ -104,6 +131,9 @@ subroutine prep_glc_init(infodata, lnd_c2_glc) allocate(mapper_Sl2g) allocate(mapper_Fl2g) + allocate(mapper_Fg2l) + + smb_renormalize = prep_glc_do_renormalize_smb(infodata) if (glc_present .and. lnd_c2_glc) then @@ -128,6 +158,7 @@ subroutine prep_glc_init(infodata, lnd_c2_glc) l2gacc_lx_cnt = 0 if (lnd_c2_glc) then + samegrid_lg = .true. if (trim(lnd_gnam) /= trim(glc_gnam)) samegrid_lg = .false. @@ -138,6 +169,7 @@ subroutine prep_glc_init(infodata, lnd_c2_glc) call seq_map_init_rcfile(mapper_Sl2g, lnd(1), glc(1), & 'seq_maps.rc', 'lnd2glc_smapname:', 'lnd2glc_smaptype:', samegrid_lg, & 'mapper_Sl2g initialization', esmf_map_flag) + if (iamroot_CPLID) then write(logunit,*) ' ' write(logunit,F00) 'Initializing mapper_Fl2g' @@ -145,6 +177,19 @@ subroutine prep_glc_init(infodata, lnd_c2_glc) call seq_map_init_rcfile(mapper_Fl2g, lnd(1), glc(1), & 'seq_maps.rc', 'lnd2glc_fmapname:', 'lnd2glc_fmaptype:', samegrid_lg, & 'mapper_Fl2g initialization', esmf_map_flag) + + ! We need to initialize our own Fg2l mapper because in some cases (particularly + ! TG compsets - dlnd forcing CISM) the system doesn't otherwise create a Fg2l + ! mapper. + if (iamroot_CPLID) then + write(logunit,*) ' ' + write(logunit,F00) 'Initializing mapper_Fg2l' + end if + call seq_map_init_rcfile(mapper_Fg2l, glc(1), lnd(1), & + 'seq_maps.rc', 'glc2lnd_fmapname:', 'glc2lnd_fmaptype:', samegrid_lg, & + 'mapper_Fg2l initialization', esmf_map_flag) + + call prep_glc_set_g2x_lx_fields() end if call shr_sys_flush(logunit) @@ -152,6 +197,96 @@ subroutine prep_glc_init(infodata, lnd_c2_glc) end subroutine prep_glc_init + !================================================================================================ + + function prep_glc_do_renormalize_smb(infodata) result(do_renormalize_smb) + ! Returns a logical saying whether we should do the smb renormalization + logical :: do_renormalize_smb ! function return value + ! + ! Arguments + type (seq_infodata_type) , intent(in) :: infodata + + ! Local variables + character(len=cl) :: glc_renormalize_smb ! namelist option saying whether to do smb renormalization + logical :: glc_coupled_fluxes ! does glc send fluxes to other components? + logical :: lnd_prognostic ! is lnd a prognostic component? + + character(len=*), parameter :: subname = '(prep_glc_do_renormalize_smb)' + !--------------------------------------------------------------- + + call seq_infodata_getdata(infodata, & + glc_renormalize_smb = glc_renormalize_smb, & + glc_coupled_fluxes = glc_coupled_fluxes, & + lnd_prognostic = lnd_prognostic) + + select case (glc_renormalize_smb) + case ('on') + do_renormalize_smb = .true. + case ('off') + do_renormalize_smb = .false. + case ('on_if_glc_coupled_fluxes') + if (.not. lnd_prognostic) then + ! Do not renormalize if we're running glc with dlnd (T compsets): In this case + ! there is no feedback from glc to lnd, and conservation is not important + do_renormalize_smb = .false. + else if (.not. glc_coupled_fluxes) then + ! Do not renormalize if glc isn't sending fluxes to other components: In this + ! case conservation is not important + do_renormalize_smb = .false. + else + ! lnd_prognostic is true and glc_coupled_fluxes is true + do_renormalize_smb = .true. + end if + case default + write(logunit,*) subname,' ERROR: unknown value for glc_renormalize_smb: ', & + trim(glc_renormalize_smb) + call shr_sys_abort(subname//' ERROR: unknown value for glc_renormalize_smb') + end select + end function prep_glc_do_renormalize_smb + + !================================================================================================ + + subroutine prep_glc_set_g2x_lx_fields() + + !--------------------------------------------------------------- + ! Description + ! Sets the module-level g2x_lx_fields variable. + ! + ! This gives the fields needed in the g2x_lx attribute vector used as part of mapping + ! qice from lnd to glc. + ! + ! Local Variables + character(len=GLC_ELEVCLASS_STRLEN), allocatable :: all_elevclass_strings(:) + character(len=:), allocatable :: frac_fields + character(len=:), allocatable :: topo_fields + integer :: strlen + + ! 1 is probably enough, but use 10 to be safe, in case the length of the delimiter + ! changes + integer, parameter :: extra_len_for_list_merge = 10 + + character(len=*), parameter :: subname = '(prep_glc_set_g2x_lx_fields)' + !--------------------------------------------------------------- + + allocate(all_elevclass_strings(0:glc_get_num_elevation_classes())) + all_elevclass_strings = glc_all_elevclass_strings(include_zero = .true.) + frac_fields = shr_string_listFromSuffixes( & + suffixes = all_elevclass_strings, & + strBase = Sg_frac_field) + ! Sg_topo is not actually needed on the land grid in + ! prep_glc_map_qice_conservative_lnd2glc, but it is required by the current interface + ! for map_glc2lnd_ec. + topo_fields = shr_string_listFromSuffixes( & + suffixes = all_elevclass_strings, & + strBase = Sg_topo_field) + + strlen = len_trim(frac_fields) + len_trim(topo_fields) + extra_len_for_list_merge + allocate(character(len=strlen) :: g2x_lx_fields) + call shr_string_listMerge(frac_fields, topo_fields, g2x_lx_fields) + + end subroutine prep_glc_set_g2x_lx_fields + + !================================================================================================ subroutine prep_glc_accum(timer) @@ -318,21 +453,35 @@ subroutine prep_glc_merge( l2x_g, fractions_g, x2g_g ) index_lfrac = mct_aVect_indexRA(fractions_g,"lfrac") do i = 1, num_flux_fields + call seq_flds_getField(field, i, seq_flds_x2g_fluxes) index_l2x = mct_aVect_indexRA(l2x_g, trim(field)) index_x2g = mct_aVect_indexRA(x2g_g, trim(field)) - if (first_time) then - mrgstr(mrgstr_index) = subname//'x2g%'//trim(field)//' =' // & - ' = lfrac*l2x%'//trim(field) - end if + if (trim(field) == qice_fieldname) then - do n = 1, lsize - lfrac = fractions_g%rAttr(index_lfrac,n) - x2g_g%rAttr(index_x2g,n) = l2x_g%rAttr(index_l2x,n) * lfrac - end do + if (first_time) then + mrgstr(mrgstr_index) = subname//'x2g%'//trim(field)//' =' // & + ' = l2x%'//trim(field) + end if + + ! treat qice as if it were a state variable, with a simple copy. + do n = 1, lsize + x2g_g%rAttr(index_x2g,n) = l2x_g%rAttr(index_l2x,n) + end do + + else + write(logunit,*) subname,' ERROR: Flux fields other than ', & + qice_fieldname, ' currently are not handled in lnd2glc remapping.' + write(logunit,*) '(Attempt to handle flux field <', trim(field), '>.)' + write(logunit,*) 'Substantial thought is needed to determine how to remap other fluxes' + write(logunit,*) 'in a smooth, conservative manner.' + call shr_sys_abort(subname//& + ' ERROR: Flux fields other than qice currently are not handled in lnd2glc remapping.') + endif ! qice_fieldname mrgstr_index = mrgstr_index + 1 + end do if (first_time) then @@ -383,35 +532,55 @@ subroutine prep_glc_calc_l2x_gx(fractions_lx, timer) do field_num = 1, num_flux_fields call seq_flds_getField(fieldname, field_num, seq_flds_x2g_fluxes) - call prep_glc_map_one_field_lnd2glc(egi=egi, eli=eli, & - fieldname = fieldname, & - fractions_lx = fractions_lx(efi), & - mapper = mapper_Fl2g) + + if (trim(fieldname) == qice_fieldname) then + + ! Use a bilinear (Sl2g) mapper, as for states. + ! The Fg2l mapper is needed to map some glc fields to the land grid + ! for purposes of conservation. + call prep_glc_map_qice_conservative_lnd2glc(egi=egi, eli=eli, & + fractions_lx = fractions_lx(efi), & + mapper_Sl2g = mapper_Sl2g, & + mapper_Fg2l = mapper_Fg2l) + + else + write(logunit,*) subname,' ERROR: Flux fields other than ', & + qice_fieldname, ' currently are not handled in lnd2glc remapping.' + write(logunit,*) '(Attempt to handle flux field <', trim(fieldname), '>.)' + write(logunit,*) 'Substantial thought is needed to determine how to remap other fluxes' + write(logunit,*) 'in a smooth, conservative manner.' + call shr_sys_abort(subname//& + ' ERROR: Flux fields other than qice currently are not handled in lnd2glc remapping.') + endif ! qice_fieldname + end do do field_num = 1, num_state_fields call seq_flds_getField(fieldname, field_num, seq_flds_x2g_states) - call prep_glc_map_one_field_lnd2glc(egi=egi, eli=eli, & + call prep_glc_map_one_state_field_lnd2glc(egi=egi, eli=eli, & fieldname = fieldname, & fractions_lx = fractions_lx(efi), & mapper = mapper_Sl2g) end do - enddo + + enddo ! egi + call t_drvstopf (trim(timer)) + end subroutine prep_glc_calc_l2x_gx !================================================================================================ - subroutine prep_glc_map_one_field_lnd2glc(egi, eli, fieldname, fractions_lx, mapper) + subroutine prep_glc_map_one_state_field_lnd2glc(egi, eli, fieldname, fractions_lx, mapper) ! Maps a single field from the land grid to the glc grid. ! - ! Note that we remap each field separately because each field needs its own - ! vertical gradient calculator. + ! This mapping is not conservative, so should only be used for state fields. + ! + ! NOTE(wjs, 2017-05-10) We used to map each field separately because each field needed + ! its own vertical gradient calculator. Now that we don't need vertical gradient + ! calculators, we may be able to change this to map multiple fields at once, at least + ! for part of map_lnd2glc. - use vertical_gradient_calculator_2nd_order, only : vertical_gradient_calculator_2nd_order_type - use vertical_gradient_calculator_factory - use glc_elevclass_mod, only : glc_get_num_elevation_classes, & - glc_get_elevclass_bounds, glc_all_elevclass_strings use map_lnd2glc_mod, only : map_lnd2glc ! Arguments @@ -422,27 +591,19 @@ subroutine prep_glc_map_one_field_lnd2glc(egi, eli, fieldname, fractions_lx, map type(seq_map), intent(inout) :: mapper ! ! Local Variables - type(mct_avect), pointer :: g2x_gx - type(vertical_gradient_calculator_2nd_order_type) :: gradient_calculator + type(mct_avect), pointer :: g2x_gx ! glc export, glc grid, cpl pes - allocated in driver !--------------------------------------------------------------- g2x_gx => component_get_c2x_cx(glc(egi)) - gradient_calculator = create_vertical_gradient_calculator_2nd_order( & - attr_vect = l2gacc_lx(eli), & - fieldname = fieldname, & - toponame = 'Sl_topo', & - elevclass_names = glc_all_elevclass_strings(), & - elevclass_bounds = glc_get_elevclass_bounds()) call map_lnd2glc(l2x_l = l2gacc_lx(eli), & landfrac_l = fractions_lx, & g2x_g = g2x_gx, & fieldname = fieldname, & - gradient_calculator = gradient_calculator, & mapper = mapper, & l2x_g = l2x_gx(eli)) - end subroutine prep_glc_map_one_field_lnd2glc + end subroutine prep_glc_map_one_state_field_lnd2glc !================================================================================================ @@ -471,6 +632,451 @@ end subroutine prep_glc_zero_fields !================================================================================================ + subroutine prep_glc_map_qice_conservative_lnd2glc(egi, eli, fractions_lx, & + mapper_Sl2g, mapper_Fg2l) + + ! Maps the surface mass balance field (qice) from the land grid to the glc grid. + ! + ! Use a smooth, non-conservative (bilinear) mapping, followed by a correction for + ! conservation. + ! + ! For high-level design, see: + ! https://docs.google.com/document/d/1H_SuK6SfCv1x6dK91q80dFInPbLYcOkUj_iAa6WRnqQ/edit + + use map_lnd2glc_mod, only : map_lnd2glc + + ! Arguments + integer, intent(in) :: egi ! glc instance index + integer, intent(in) :: eli ! lnd instance index + type(mct_aVect) , intent(in) :: fractions_lx ! fractions on the land grid, for this frac instance + type(seq_map), intent(inout) :: mapper_Sl2g ! state mapper from land to glc grid; non-conservative + type(seq_map), intent(inout) :: mapper_Fg2l ! flux mapper from glc to land grid; conservative + ! + ! Local Variables + type(mct_aVect), pointer :: g2x_gx ! glc export, glc grid + + logical :: iamroot + + !Note: The sums in this subroutine use the coupler areas aream_l and aream_g. + ! The coupler areas can differ from the native areas area_l and area_g. + ! (For CISM with a polar stereographic projection, area_g can differ from aream_g + ! by up to ~10%.) + ! If so, then the calls to subroutine mct_avect_vecmult in component_mod.F90 + ! (just before and after the call to comp_run) should adjust the SMB fluxes + ! such that in each grid cell, the native value of area*flux is equal to the + ! coupler value of aream*flux. This assumes that the SMB field is contained in + ! seq_fields l2x_fluxes and seq_fields_x2g_fluxes. + + real(r8), dimension(:), allocatable :: aream_g ! cell areas on glc grid, for mapping + real(r8), dimension(:), allocatable :: area_g ! cell areas on glc grid, according to glc model + + type(mct_ggrid), pointer :: dom_g ! glc grid info + + integer :: lsize_g ! number of points on glc grid + + integer :: n + integer :: km, ka + + real(r8), pointer :: qice_g(:) ! qice data on glc grid + + !--------------------------------------------------------------- + + call seq_comm_getdata(CPLID, iamroot=iamroot) + + if (iamroot) then + write(logunit,*) ' ' + write(logunit,*) 'In prep_glc_map_qice_conservative_lnd2glc' + write(logunit,*) 'smb_renormalize = ', smb_renormalize + endif + + ! Get attribute vector needed for mapping and conservation + g2x_gx => component_get_c2x_cx(glc(egi)) + + ! get grid size + lsize_g = mct_aVect_lsize(l2x_gx(eli)) + + ! allocate and fill area arrays on the glc grid + ! (Note that we get domain information from instance 1, following what's done in + ! other parts of the coupler.) + dom_g => component_get_dom_cx(glc(1)) + + allocate(aream_g(lsize_g)) + km = mct_aVect_indexRa(dom_g%data, "aream" ) + aream_g(:) = dom_g%data%rAttr(km,:) + + allocate(area_g(lsize_g)) + ka = mct_aVect_indexRa(dom_g%data, "area" ) + area_g(:) = dom_g%data%rAttr(ka,:) + + ! Map the SMB from the land grid to the glc grid, using a non-conservative state mapper. + call map_lnd2glc(l2x_l = l2gacc_lx(eli), & + landfrac_l = fractions_lx, & + g2x_g = g2x_gx, & + fieldname = qice_fieldname, & + mapper = mapper_Sl2g, & + l2x_g = l2x_gx(eli)) + + ! Export the remapped SMB to a local array + allocate(qice_g(lsize_g)) + call mct_aVect_exportRattr(l2x_gx(eli), trim(qice_fieldname), qice_g) + + ! Make a preemptive adjustment to qice_g to account for area differences between CISM and the coupler. + ! In component_mod.F90, there is a call to mct_avect_vecmult, which multiplies the fluxes + ! by aream_g/area_g for conservation purposes. Where CISM areas are larger (area_g > aream_g), + ! the fluxes are reduced, and where CISM areas are smaller, the fluxes are increased. + ! As a result, an SMB of 1 m/yr in CLM would be converted to an SMB ranging from + ! ~0.9 to 1.05 m/yr in CISM (with smaller values where CISM areas are larger, and larger + ! values where CISM areas are smaller). + ! Here, to keep CISM values close to the CLM values in the corresponding locations, + ! we anticipate the later correction and multiply qice_g by area_g/aream_g. + ! Then the later call to mct_avect_vecmult will bring qice back to the original values + ! obtained from bilinear remapping. + ! If Flgl_qice were changed to a state (and not included in seq_flds_x2g_fluxes), + ! then we could skip this adjustment. + ! + ! Note that we are free to do this or any other adjustments we want to qice at this + ! point in the remapping, because the conservation correction will ensure that we + ! still conserve globally despite these adjustments (and smb_renormalize = .false. + ! should only be used in cases where conservation doesn't matter anyway). + + do n = 1, lsize_g + if (aream_g(n) > 0.0_r8) then + qice_g(n) = qice_g(n) * area_g(n)/aream_g(n) + else + qice_g(n) = 0.0_r8 + endif + enddo + + if (smb_renormalize) then + call prep_glc_renormalize_smb( & + eli = eli, & + fractions_lx = fractions_lx, & + g2x_gx = g2x_gx, & + mapper_Fg2l = mapper_Fg2l, & + aream_g = aream_g, & + qice_g = qice_g) + end if + + ! Put the adjusted SMB back into l2x_gx. + ! + ! If we are doing renormalization, then this is the renormalized SMB. Whether or not + ! we are doing renormalization, this captures the preemptive adjustment to qice_g to + ! account for area differences between CISM and the coupler. + call mct_aVect_importRattr(l2x_gx(eli), qice_fieldname, qice_g) + + ! clean up + + deallocate(aream_g) + deallocate(area_g) + deallocate(qice_g) + + end subroutine prep_glc_map_qice_conservative_lnd2glc + + !================================================================================================ + + subroutine prep_glc_renormalize_smb(eli, fractions_lx, g2x_gx, mapper_Fg2l, aream_g, qice_g) + + ! Renormalizes surface mass balance (smb, here named qice_g) so that the global + ! integral on the glc grid is equal to the global integral on the land grid. + ! + ! This is required for conservation - although conservation is only necessary if we + ! are running with a fully-interactive, two-way-coupled glc. + ! + ! For high-level design, see: + ! https://docs.google.com/document/d/1H_SuK6SfCv1x6dK91q80dFInPbLYcOkUj_iAa6WRnqQ/edit + + use map_glc2lnd_mod, only : map_glc2lnd_ec + + ! Arguments + integer , intent(in) :: eli ! lnd instance index + type(mct_aVect) , intent(in) :: fractions_lx ! fractions on the land grid, for this frac instance + type(mct_aVect) , intent(in) :: g2x_gx ! glc export, glc grid + type(seq_map) , intent(inout) :: mapper_Fg2l ! flux mapper from glc to land grid; conservative + real(r8) , intent(in) :: aream_g(:) ! cell areas on glc grid, for mapping + real(r8) , intent(inout) :: qice_g(:) ! qice data on glc grid + + ! + ! Local Variables + integer :: mpicom + logical :: iamroot + + type(mct_ggrid), pointer :: dom_l ! land grid info + + integer :: lsize_l ! number of points on land grid + integer :: lsize_g ! number of points on glc grid + + real(r8), dimension(:), allocatable :: aream_l ! cell areas on land grid, for mapping + + real(r8), pointer :: qice_l(:,:) ! SMB (Flgl_qice) on land grid + real(r8), pointer :: frac_l(:,:) ! EC fractions (Sg_ice_covered) on land grid + real(r8), pointer :: tmp_field_l(:) ! temporary field on land grid + + ! The following need to be pointers to satisfy the MCT interface + ! Note: Sg_icemask defines where the ice sheet model can receive a nonzero SMB from the land model. + real(r8), pointer :: Sg_icemask_g(:) ! icemask on glc grid + real(r8), pointer :: Sg_icemask_l(:) ! icemask on land grid + real(r8), pointer :: lfrac(:) ! land fraction on land grid + + type(mct_aVect) :: g2x_lx ! glc export, lnd grid (not a pointer: created locally) + type(mct_avect) :: Sg_icemask_l_av ! temporary attribute vector holding Sg_icemask on the land grid + + integer :: nEC ! number of elevation classes + integer :: n + integer :: ec + integer :: km + + ! various strings for building field names + character(len=:), allocatable :: elevclass_as_string + character(len=:), allocatable :: qice_field + character(len=:), allocatable :: frac_field + + ! local and global sums of accumulation and ablation; used to compute renormalization factors + + real(r8) :: local_accum_on_land_grid + real(r8) :: global_accum_on_land_grid + real(r8) :: local_accum_on_glc_grid + real(r8) :: global_accum_on_glc_grid + + real(r8) :: local_ablat_on_land_grid + real(r8) :: global_ablat_on_land_grid + real(r8) :: local_ablat_on_glc_grid + real(r8) :: global_ablat_on_glc_grid + + ! renormalization factors (should be close to 1, e.g. in range 0.95 to 1.05) + real(r8) :: accum_renorm_factor ! ratio between global accumulation on the two grids + real(r8) :: ablat_renorm_factor ! ratio between global ablation on the two grids + + real(r8) :: effective_area ! grid cell area multiplied by min(lfrac,Sg_icemask_l). + ! This is the area that can contribute SMB to the ice sheet model. + + + !--------------------------------------------------------------- + + lsize_g = size(qice_g) + SHR_ASSERT_FL((size(aream_g) == lsize_g), __FILE__, __LINE__) + + call seq_comm_setptrs(CPLID, mpicom=mpicom) + call seq_comm_getdata(CPLID, iamroot=iamroot) + lsize_l = mct_aVect_lsize(l2gacc_lx(eli)) + + ! allocate and fill area arrays on the land grid + ! (Note that we get domain information from instance 1, following what's done in + ! other parts of the coupler.) + dom_l => component_get_dom_cx(lnd(1)) + + allocate(aream_l(lsize_l)) + km = mct_aVect_indexRa(dom_l%data, "aream" ) + aream_l(:) = dom_l%data%rAttr(km,:) + + ! Export land fractions from fractions_lx to a local array + allocate(lfrac(lsize_l)) + call mct_aVect_exportRattr(fractions_lx, "lfrac", lfrac) + + ! Map Sg_icemask from the glc grid to the land grid. + ! This may not be necessary, if Sg_icemask_l has already been mapped from Sg_icemask_g. + ! It is done here for two reasons: + ! (1) The mapping will *not* have been done if we are running with dlnd (e.g., a TG case). + ! (2) Because of coupler lags, the current Sg_icemask_l might not be up to date with + ! Sg_icemask_g. This probably isn't a problem in practice, but doing the mapping + ! here ensures the mask is up to date. + ! + ! This mapping uses the same options as the standard glc -> lnd mapping done in + ! prep_lnd_calc_g2x_lx. If that mapping ever changed (e.g., changing norm to + ! .false.), then we should change this mapping, too. + ! + ! BUG(wjs, 2017-05-11, #1516) I think we actually want norm = .false. here, but this + ! requires some more thought + call mct_aVect_init(Sg_icemask_l_av, rList = Sg_icemask_field, lsize = lsize_l) + call seq_map_map(mapper = mapper_Fg2l, & + av_s = g2x_gx, & + av_d = Sg_icemask_l_av, & + fldlist = Sg_icemask_field, & + norm = .true.) + + ! Export Sg_icemask_l from the temporary attribute vector to a local array + allocate(Sg_icemask_l(lsize_l)) + call mct_aVect_exportRattr(Sg_icemask_l_av, Sg_icemask_field, Sg_icemask_l) + + ! Clean the temporary attribute vector + call mct_aVect_clean(Sg_icemask_l_av) + + ! Map Sg_ice_covered from the glc grid to the land grid. + ! This gives the fields Sg_ice_covered00, Sg_ice_covered01, etc. on the land grid. + ! These fields are needed to integrate the total SMB on the land grid, for conservation purposes. + ! As above, the mapping may not be necessary, because Sg_ice_covered might already have been mapped. + ! However, the mapping will not have been done in a TG case with dlnd, and it might not + ! be up to date because of coupler lags (though the latter probably isn't a problem + ! in practice). + ! + ! Note that, for a case with full two-way coupling, we will only conserve if the + ! actual land cover used over the course of the year matches these currently-remapped + ! values. This should generally be the case with the current coupling setup. + ! + ! One could argue that it would be safer (for conservation purposes) if LND sent its + ! grid cell average SMB values, or if it sent its own notion of the area in each + ! elevation class for the purpose of creating grid cell average SMB values here. But + ! these options cause problems if we're not doing full two-way coupling (e.g., in a TG + ! case with dlnd, or in the common case where GLC is a diagnostic component that + ! doesn't cause updates in the glacier areas in LND). In these cases without full + ! two-way coupling, if we use the LND's notion of the area in each elevation class, + ! then the conservation corrections would end up correcting for discrepancies in + ! elevation class areas between LND and GLC, rather than just correcting for + ! discrepancies arising from the remapping of SMB. (And before you get worried: It + ! doesn't matter that we are not conserving in these cases without full two-way + ! coupling, because GLC isn't connected with the rest of the system in terms of energy + ! and mass in these cases. So in these cases, it's okay that the LND integral computed + ! here differs from the integral that LND itself would compute.) + + ! Create an attribute vector g2x_lx to hold the mapped fields + call mct_aVect_init(g2x_lx, rList=g2x_lx_fields, lsize=lsize_l) + + ! Map Sg_ice_covered and Sg_topo from glc to land + call map_glc2lnd_ec( & + g2x_g = g2x_gx, & + frac_field = Sg_frac_field, & + topo_field = Sg_topo_field, & + icemask_field = Sg_icemask_field, & + extra_fields = ' ', & ! no extra fields + mapper = mapper_Fg2l, & + g2x_l = g2x_lx) + + ! Export qice and Sg_ice_covered in each elevation class to local arrays. + ! Note: qice comes from l2gacc_lx; frac comes from g2x_lx. + + nEC = glc_get_num_elevation_classes() + + allocate(qice_l(lsize_l,0:nEC)) + allocate(frac_l(lsize_l,0:nEC)) + allocate(tmp_field_l(lsize_l)) + + do ec = 0, nEC + elevclass_as_string = glc_elevclass_as_string(ec) + + frac_field = Sg_frac_field // elevclass_as_string ! Sg_ice_covered01, etc. + call mct_aVect_exportRattr(g2x_lx, trim(frac_field), tmp_field_l) + frac_l(:,ec) = tmp_field_l(:) + + qice_field = qice_fieldname // elevclass_as_string ! Flgl_qice01, etc. + call mct_aVect_exportRattr(l2gacc_lx(eli), trim(qice_field), tmp_field_l) + qice_l(:,ec) = tmp_field_l(:) + + enddo + + ! clean the temporary attribute vector g2x_lx + call mct_aVect_clean(g2x_lx) + + ! Sum qice over local land grid cells + + ! initialize qice sum + local_accum_on_land_grid = 0.0_r8 + local_ablat_on_land_grid = 0.0_r8 + + do n = 1, lsize_l + + effective_area = min(lfrac(n),Sg_icemask_l(n)) * aream_l(n) + + do ec = 0, nEC + + if (qice_l(n,ec) >= 0.0_r8) then + local_accum_on_land_grid = local_accum_on_land_grid & + + effective_area * frac_l(n,ec) * qice_l(n,ec) + else + local_ablat_on_land_grid = local_ablat_on_land_grid & + + effective_area * frac_l(n,ec) * qice_l(n,ec) + endif + + enddo ! ec + + enddo ! n + + call shr_mpi_sum(local_accum_on_land_grid, & + global_accum_on_land_grid, & + mpicom, 'accum_l') + + call shr_mpi_sum(local_ablat_on_land_grid, & + global_ablat_on_land_grid, & + mpicom, 'ablat_l') + + call shr_mpi_bcast(global_accum_on_land_grid, mpicom) + call shr_mpi_bcast(global_ablat_on_land_grid, mpicom) + + ! Sum qice_g over local glc grid cells. + ! Note: This sum uses the coupler areas (aream_g), which differ from the native CISM areas. + ! But since the original qice_g (from bilinear remapping) has been multiplied by + ! area_g/aream_g above, this calculation is equivalent to multiplying the original qice_g + ! by the native CISM areas (area_g). + ! If Flgl_qice were changed to a state (and not included in seq_flds_x2g_fluxes), + ! then it would be appropriate to use the native CISM areas in this sum. + + ! Export Sg_icemask from g2x_gx to a local array + allocate(Sg_icemask_g(lsize_g)) + call mct_aVect_exportRattr(g2x_gx, Sg_icemask_field, Sg_icemask_g) + + local_accum_on_glc_grid = 0.0_r8 + local_ablat_on_glc_grid = 0.0_r8 + + do n = 1, lsize_g + + if (qice_g(n) >= 0.0_r8) then + local_accum_on_glc_grid = local_accum_on_glc_grid & + + Sg_icemask_g(n) * aream_g(n) * qice_g(n) + else + local_ablat_on_glc_grid = local_ablat_on_glc_grid & + + Sg_icemask_g(n) * aream_g(n) * qice_g(n) + endif + + enddo ! n + + call shr_mpi_sum(local_accum_on_glc_grid, & + global_accum_on_glc_grid, & + mpicom, 'accum_g') + + call shr_mpi_sum(local_ablat_on_glc_grid, & + global_ablat_on_glc_grid, & + mpicom, 'ablat_g') + + call shr_mpi_bcast(global_accum_on_glc_grid, mpicom) + call shr_mpi_bcast(global_ablat_on_glc_grid, mpicom) + + ! Renormalize + + if (global_accum_on_glc_grid > 0.0_r8) then + accum_renorm_factor = global_accum_on_land_grid / global_accum_on_glc_grid + else + accum_renorm_factor = 0.0_r8 + endif + + if (global_ablat_on_glc_grid < 0.0_r8) then ! negative by definition + ablat_renorm_factor = global_ablat_on_land_grid / global_ablat_on_glc_grid + else + ablat_renorm_factor = 0.0_r8 + endif + + if (iamroot) then + write(logunit,*) 'accum_renorm_factor = ', accum_renorm_factor + write(logunit,*) 'ablat_renorm_factor = ', ablat_renorm_factor + endif + + do n = 1, lsize_g + if (qice_g(n) >= 0.0_r8) then + qice_g(n) = qice_g(n) * accum_renorm_factor + else + qice_g(n) = qice_g(n) * ablat_renorm_factor + endif + enddo + + deallocate(aream_l) + deallocate(lfrac) + deallocate(Sg_icemask_l) + deallocate(Sg_icemask_g) + deallocate(tmp_field_l) + deallocate(qice_l) + deallocate(frac_l) + + end subroutine prep_glc_renormalize_smb + + !================================================================================================ + function prep_glc_get_l2x_gx() type(mct_aVect), pointer :: prep_glc_get_l2x_gx(:) prep_glc_get_l2x_gx => l2x_gx(:) diff --git a/src/drivers/mct/main/prep_lnd_mod.F90 b/src/drivers/mct/main/prep_lnd_mod.F90 index 204b3c17f77..74282df1356 100644 --- a/src/drivers/mct/main/prep_lnd_mod.F90 +++ b/src/drivers/mct/main/prep_lnd_mod.F90 @@ -451,6 +451,13 @@ subroutine prep_lnd_calc_g2x_lx(timer) ! These are mapped using a simple area-conservative remapping. (Note that we use ! the flux mapper even though these contain states, because we need these icemask ! fields to be mapped conservatively.) + ! + ! Note that this mapping is redone for Sg_icemask in prep_glc_mod: + ! prep_glc_map_qice_conservative_lnd2glc. If we ever change this mapping (e.g., + ! changing norm to .false.), then we should change the mapping there, too. + ! + ! BUG(wjs, 2017-05-11, #1516) I think we actually want norm = .false. here, but + ! this requires some more thought call seq_map_map(mapper_Fg2l, g2x_gx, g2x_lx(egi), & fldlist = glc2lnd_non_ec_fields, norm=.true.) diff --git a/src/drivers/mct/main/vertical_gradient_calculator_2nd_order.F90 b/src/drivers/mct/main/vertical_gradient_calculator_2nd_order.F90 deleted file mode 100644 index e53a7894d1d..00000000000 --- a/src/drivers/mct/main/vertical_gradient_calculator_2nd_order.F90 +++ /dev/null @@ -1,410 +0,0 @@ -module vertical_gradient_calculator_2nd_order - - !--------------------------------------------------------------------- - ! - ! Purpose: - ! - ! This module defines a subclass of vertical_gradient_calculator_base_type for - ! computing vertical gradients using a second-order centered difference. - ! - ! If the topo values are nearly equal across the gradient (i.e., denominator is near 0), - ! returns a gradient of 0. - ! - ! If there is only one elevation class, returns a gradient of 0. - -#include "shr_assert.h" - use seq_comm_mct, only : logunit - use vertical_gradient_calculator_base, only : vertical_gradient_calculator_base_type - use shr_kind_mod, only : r8 => shr_kind_r8 - use shr_log_mod, only : errMsg => shr_log_errMsg - use shr_sys_mod, only : shr_sys_abort - use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=) - implicit none - private - - public :: vertical_gradient_calculator_2nd_order_type - - type, extends(vertical_gradient_calculator_base_type) :: & - vertical_gradient_calculator_2nd_order_type - private - - integer :: nelev ! number of elevation classes - integer :: num_points - real(r8), allocatable :: field(:,:) ! field(i,j) is point i, elevation class j - real(r8), allocatable :: topo(:,:) ! topo(i,j) is point i, elevation class j - - ! Bounds of each elevation class. This array has one more element than the number of - ! elevation classes, since it contains lower and upper bounds for each elevation - ! class. The indices go 0:nelev. These bounds - ! are guaranteed to be monotonically increasing. - real(r8), allocatable :: elevclass_bounds(:) - - ! precomputed vertical gradients; vertical_gradient(i,j) is point i, elevation class - ! j - real(r8), allocatable :: vertical_gradient(:,:) - - logical :: calculated ! whether gradients have been calculated yet - - contains - procedure :: calc_gradients - procedure :: get_gradients_one_class - procedure :: get_gradients_one_point - - procedure, private :: check_topo ! check topographic heights - procedure, private :: limit_gradient - - end type vertical_gradient_calculator_2nd_order_type - - interface vertical_gradient_calculator_2nd_order_type - module procedure constructor - end interface vertical_gradient_calculator_2nd_order_type - -contains - - !----------------------------------------------------------------------- - function constructor(field, topo, elevclass_bounds) result(this) - ! - ! !DESCRIPTION: - ! Creates a vertical_gradient_calculator_2nd_order_type object. - ! - ! Pre-condition: elevclass_bounds must be monotonically increasing. - ! - ! Pre-condition: Topographic heights must all lie inside the bounds of their - ! respective elevation class (given by elevclass_bounds), with the possible exception - ! of the lowest elevation class (topographic heights can lie below the arbitrary lower - ! bound of the elevation class) and the highest elevation class (topographic heights - ! can lie above the arbitrary upper bound of the elevation class). (This pre-condition - ! is mainly important for the sake of calculating the limiter.) - ! - ! TODO(wjs, 2016-04-21) Currently the topographic heights pre-condition is not - ! checked: see below. - ! - ! !USES: - ! - ! !ARGUMENTS: - type(vertical_gradient_calculator_2nd_order_type) :: this ! function result - real(r8), intent(in) :: field(:,:) ! field(i,j) is point i, elevation class j - real(r8), intent(in) :: topo(:,:) ! topo(i,j) is point i, elevation class j - - ! bounds of each elevation class; this array should have one more element than the - ! number of elevation classes, since it contains lower and upper bounds for each - ! elevation class - real(r8), intent(in) :: elevclass_bounds(0:) - ! - ! !LOCAL VARIABLES: - - character(len=*), parameter :: subname = 'constructor' - !----------------------------------------------------------------------- - - this%calculated = .false. - - this%num_points = size(field, 1) - this%nelev = size(field, 2) - SHR_ASSERT_ALL_FL((ubound(topo) == (/this%num_points, this%nelev/)), __FILE__, __LINE__) - SHR_ASSERT_ALL_FL((ubound(elevclass_bounds) == (/this%nelev/)), __FILE__, __LINE__) - - allocate(this%elevclass_bounds(0:this%nelev)) - this%elevclass_bounds(:) = elevclass_bounds(:) - - ! (In principle, we could also handle monotonically decreasing elevclass_bounds, but - ! that would require generalizing some code, such as in check_topo.) - call this%check_elevclass_bounds_monotonic_increasing(this%elevclass_bounds) - - allocate(this%field(this%num_points, this%nelev)) - this%field(:,:) = field(:,:) - allocate(this%topo(this%num_points, this%nelev)) - this%topo(:,:) = topo(:,:) - - ! TODO(wjs, 2016-04-21) Uncomment this call to check_topo. It's important for - ! topographic heights to be within bounds in order for the limiter to be applied - ! correctly. However, this currently isn't the case for some of the old TG forcing - ! data. At a glance, it looks like the problems are just outside of Greenland, so this - ! should be okay. When we have new TG forcing data, we should try uncommenting this - ! call to check_topo. - ! - ! Alternatively, we could change check_topo to set a flag for each point saying - ! whether topo values are bad for that point. Then, when computing gradients, set - ! them to 0 for all points with bad topo values. (We did something similar for the - ! now-deleted vertical_gradient_calculator_continuous.) However, longer-term, an - ! abort may be more appropriate rather than silently setting gradients to 0. - - ! call this%check_topo() - - allocate(this%vertical_gradient(this%num_points, this%nelev)) - this%vertical_gradient(:,:) = nan - - end function constructor - - !----------------------------------------------------------------------- - subroutine calc_gradients(this) - ! - ! !DESCRIPTION: - ! Calculates all vertical gradients - ! - ! !USES: - ! - ! !ARGUMENTS: - class(vertical_gradient_calculator_2nd_order_type), intent(inout) :: this - ! - ! !LOCAL VARIABLES: - ! Tolerance for considering two topo values to be nearly equal - real(r8), parameter :: topo_equality_tolerance = 1.e-13_r8 - - integer :: i - integer :: elevation_class - integer :: ec_low ! elevation class index to use as the lower bound of the gradient - integer :: ec_high ! elevation class index to use as the upper bound of the gradient - logical :: two_sided ! true if we're estimating the gradient with a two-sided difference - - character(len=*), parameter :: subname = 'calc_gradients' - !----------------------------------------------------------------------- - - if (this%calculated) then - ! nothing to do - return - end if - - if (this%nelev == 1) then - this%vertical_gradient(:,:) = 0._r8 - two_sided = .false. - - else - - do elevation_class = 1, this%nelev - if (elevation_class == 1) then - ec_low = elevation_class - ec_high = elevation_class + 1 - two_sided = .false. - else if (elevation_class == this%nelev) then - ec_low = elevation_class - 1 - ec_high = elevation_class - two_sided = .false. - else - ec_low = elevation_class - 1 - ec_high = elevation_class + 1 - two_sided = .true. - end if - - do i = 1, this%num_points - if (abs(this%topo(i, ec_high) - this%topo(i, ec_low)) < topo_equality_tolerance) then - this%vertical_gradient(i, elevation_class) = 0._r8 - else - this%vertical_gradient(i, elevation_class) = & - (this%field(i, ec_high) - this%field(i, ec_low)) / & - (this%topo (i, ec_high) - this%topo (i, ec_low)) - end if - end do - - if (two_sided) then - call this%limit_gradient(elevation_class, ec_low, ec_high, & - this%vertical_gradient(:,elevation_class)) - end if - end do - - end if - - this%calculated = .true. - - end subroutine calc_gradients - - !----------------------------------------------------------------------- - subroutine get_gradients_one_class(this, elevation_class, gradients) - ! - ! !DESCRIPTION: - ! Returns the vertical gradient for all points, at a given elevation class. - ! - ! this%calc_gradients should already have been called - ! - ! !USES: - ! - ! !ARGUMENTS: - class(vertical_gradient_calculator_2nd_order_type), intent(in) :: this - integer, intent(in) :: elevation_class - - ! gradients should already be allocated to the appropriate size - real(r8), intent(out) :: gradients(:) - ! - ! !LOCAL VARIABLES: - - character(len=*), parameter :: subname = 'get_gradients_one_class' - !----------------------------------------------------------------------- - - ! Assert pre-conditions - - SHR_ASSERT_FL(this%calculated, __FILE__, __LINE__) - SHR_ASSERT_FL((size(gradients) == this%num_points), __FILE__, __LINE__) - - if (elevation_class < 1 .or. elevation_class > this%nelev) then - write(logunit,*) subname, ': ERROR: elevation class out of bounds: ', & - elevation_class, this%nelev - call shr_sys_abort(subname//': ERROR: elevation class out of bounds') - end if - - gradients(:) = this%vertical_gradient(:, elevation_class) - - end subroutine get_gradients_one_class - - !----------------------------------------------------------------------- - subroutine get_gradients_one_point(this, point, gradients) - ! - ! !DESCRIPTION: - ! Returns the vertical gradient for all elevation classes, for one point - ! - ! this%calc_gradients should already have been called - ! - ! !USES: - ! - ! !ARGUMENTS: - class(vertical_gradient_calculator_2nd_order_type), intent(in) :: this - integer, intent(in) :: point - - ! gradients should already be allocated to the appropriate size - real(r8), intent(out) :: gradients(:) - ! - ! !LOCAL VARIABLES: - - character(len=*), parameter :: subname = 'get_gradients_one_point' - !----------------------------------------------------------------------- - - SHR_ASSERT_FL(this%calculated, __FILE__, __LINE__) - SHR_ASSERT_FL(point <= this%num_points, __FILE__, __LINE__) - SHR_ASSERT_FL((size(gradients) == this%nelev), __FILE__, __LINE__) - - gradients(:) = this%vertical_gradient(point, :) - - end subroutine get_gradients_one_point - - !----------------------------------------------------------------------- - subroutine check_topo(this) - ! - ! !DESCRIPTION: - ! Check topographic heights; abort if there is a problem - ! - ! Topographic heights in the attribute vector must all lie inside the bounds of their - ! respective elevation class (given by elevclass_bounds), with the possible exception - ! of the lowest elevation class (topographic heights can lie below the arbitrary lower - ! bound of the elevation class) and the highest elevation class (topographic heights - ! can lie above the arbitrary upper bound of the elevation class) - ! - ! !USES: - ! - ! !ARGUMENTS: - class(vertical_gradient_calculator_2nd_order_type), intent(in) :: this - ! - ! !LOCAL VARIABLES: - integer :: elevclass - integer :: i - - ! Absolute tolerance for error checks. This is chosen so that it allows for - ! double-precision roundoff-level errors on values of order 10,000. - real(r8), parameter :: tol = 1.e-10_r8 - - character(len=*), parameter :: subname = 'check_topo' - !----------------------------------------------------------------------- - - do elevclass = 1, this%nelev - if (elevclass > 1) then - do i = 1, this%num_points - if (this%topo(i,elevclass) - this%elevclass_bounds(elevclass-1) < -tol) then - write(logunit,*) subname, ': ERROR: topo lower than lower bound of elevation class:' - write(logunit,*) 'i, elevclass, topo, lower_bound = ', & - i, elevclass, this%topo(i,elevclass), this%elevclass_bounds(elevclass-1) - call shr_sys_abort(subname//': ERROR: topo lower than lower bound of elevation class') - end if - end do - end if - - if (elevclass < this%nelev) then - do i = 1, this%num_points - if (this%topo(i,elevclass) - this%elevclass_bounds(elevclass) > tol) then - write(logunit,*) subname, ': ERROR: topo higher than upper bound of elevation class:' - write(logunit,*) 'i, elevclass, topo, upper_bound = ', & - i, elevclass, this%topo(i,elevclass), this%elevclass_bounds(elevclass) - call shr_sys_abort(subname//': ERROR: topo higher than upper bound of elevation class') - end if - end do - end if - end do - - end subroutine check_topo - - !----------------------------------------------------------------------- - subroutine limit_gradient(this, k, ec_low, ec_high, vertical_gradient) - ! - ! !DESCRIPTION: - ! Limit the gradient: Ensure that the interface values lie inside the range defined - ! by the max and min of the mean values in this class and its 2 adjacent neighbors. - ! - ! Should only be called for two-sided differences (ec_low < k < ec_high) - ! - ! !ARGUMENTS: - class(vertical_gradient_calculator_2nd_order_type), intent(in) :: this - integer , intent(in) :: k ! elevation class index - integer , intent(in) :: ec_low ! elevation class index used as the lower bound of the gradient - integer , intent(in) :: ec_high ! elevation class index used as the upper bound of the gradient - real(r8), intent(inout) :: vertical_gradient(:) - ! - ! !LOCAL VARIABLES: - integer :: i - real(r8) :: deviation_high - real(r8) :: deviation_low - real(r8) :: deviation_max - real(r8) :: deviation_min - real(r8) :: diff_max - real(r8) :: diff_min - real(r8) :: factor1 - real(r8) :: factor2 - real(r8) :: limiting_factor - - character(len=*), parameter :: subname = 'limit_gradient' - !----------------------------------------------------------------------- - - ! Basic idea: In 1D with a linear reconstruction, the extreme values of the data will - ! lie at the interfaces between adjacent elevation classes. The interface values - ! should not lie outside the range defined by the max and min of the mean values in - ! this class and its 2 adjacent neighbors. - - ! This code only works correctly if we're doing a two-sided difference (otherwise, - ! one of diff_min or diff_max will be 0, leading to 0 gradient - when in fact we - ! don't want to do any limiting for a one-sided difference). - SHR_ASSERT(ec_low < k, subname//': Only works for two-sided difference: must have ec_low < k') - SHR_ASSERT(ec_high > k, subname//': Only works for two-sided difference: must have ec_high > k') - - do i = 1, this%num_points - ! First compute the max and min values of the deviation of the data from its mean - ! value. With a linear gradient, the max differences must lie at the adjacent - ! interfaces. - deviation_high = vertical_gradient(i) * (this%elevclass_bounds(k) - this%topo(i,k)) - deviation_low = vertical_gradient(i) * (this%elevclass_bounds(k-1) - this%topo(i,k)) - deviation_max = max(deviation_high, deviation_low) - deviation_min = min(deviation_high, deviation_low) - - ! Now compute the max and min of the data in the cell and its nearest neighbors. - ! (Actually, the difference between this max/min value and the mean value in the - ! current class.) - diff_max = max(this%field(i,ec_high), this%field(i,k), this%field(i,ec_low)) - this%field(i,k) - diff_min = min(this%field(i,ec_high), this%field(i,k), this%field(i,ec_low)) - this%field(i,k) - - ! Now limit the gradient using the information computed above. - - if (abs(deviation_min) > 0._r8) then - factor1 = max(0._r8, diff_min/deviation_min) - else - factor1 = 1._r8 - endif - - if (abs(deviation_max) > 0._r8) then - factor2 = max(0._r8, diff_max/deviation_max) - else - factor2 = 1._r8 - endif - - ! limiting factor will lie between 0 and 1 - limiting_factor = min(1._r8, factor1, factor2) - vertical_gradient(i) = vertical_gradient(i) * limiting_factor - end do - - end subroutine limit_gradient - -end module vertical_gradient_calculator_2nd_order - diff --git a/src/drivers/mct/main/vertical_gradient_calculator_base.F90 b/src/drivers/mct/main/vertical_gradient_calculator_base.F90 deleted file mode 100644 index 27fb7a7aa41..00000000000 --- a/src/drivers/mct/main/vertical_gradient_calculator_base.F90 +++ /dev/null @@ -1,100 +0,0 @@ -module vertical_gradient_calculator_base - - !--------------------------------------------------------------------- - ! - ! Purpose: - ! - ! This module defines an abstract base class for computing the vertical gradient of a - ! field. - ! - ! Usage: - ! - ! - First call calc_gradients - ! - ! - Then can query the computed vertical gradients using the other methods - - use seq_comm_mct, only : logunit - use shr_kind_mod, only : r8 => shr_kind_r8 - use shr_sys_mod , only : shr_sys_abort - - implicit none - private - - public :: vertical_gradient_calculator_base_type - - type, abstract :: vertical_gradient_calculator_base_type - contains - ! Calculate the vertical gradients for all points and all elevation classes - procedure(calc_gradients_interface), deferred :: calc_gradients - - ! Get the vertical gradients for all points for a single elevation class - procedure(get_gradients_one_class_interface), deferred :: get_gradients_one_class - - ! Get the vertical gradients for all elevation classes for a single point - procedure(get_gradients_one_point_interface), deferred :: get_gradients_one_point - - ! These routines are utility methods for derived classes; they should not be called - ! by clients of this class. - procedure, nopass :: check_elevclass_bounds_monotonic_increasing - - - end type vertical_gradient_calculator_base_type - - abstract interface - subroutine calc_gradients_interface(this) - import :: vertical_gradient_calculator_base_type - class(vertical_gradient_calculator_base_type), intent(inout) :: this - end subroutine calc_gradients_interface - - subroutine get_gradients_one_class_interface(this, elevation_class, gradients) - import :: vertical_gradient_calculator_base_type - import :: r8 - class(vertical_gradient_calculator_base_type), intent(in) :: this - integer, intent(in) :: elevation_class - - ! vertical_gradient should already be allocated to the appropriate size - real(r8), intent(out) :: gradients(:) - end subroutine get_gradients_one_class_interface - - subroutine get_gradients_one_point_interface(this, point, gradients) - import :: vertical_gradient_calculator_base_type - import :: r8 - class(vertical_gradient_calculator_base_type), intent(in) :: this - integer, intent(in) :: point - - ! vertical_gradient should already be allocated to the appropriate size - real(r8), intent(out) :: gradients(:) - end subroutine get_gradients_one_point_interface - end interface - -contains - - !----------------------------------------------------------------------- - subroutine check_elevclass_bounds_monotonic_increasing(elevclass_bounds) - ! - ! !DESCRIPTION: - ! Ensure that elevclass_bounds are monotonically increasing; abort if there is a - ! problem - ! - ! !ARGUMENTS: - real(r8), intent(in) :: elevclass_bounds(:) - ! - ! !LOCAL VARIABLES: - integer :: i - - character(len=*), parameter :: subname = 'check_elevclass_bounds' - !----------------------------------------------------------------------- - - do i = 2, size(elevclass_bounds) - if (elevclass_bounds(i-1) >= elevclass_bounds(i)) then - write(logunit,*) subname, ': ERROR: elevclass_bounds must be monotonically increasing' - write(logunit,*) 'elevclass_bounds = ', elevclass_bounds - call shr_sys_abort(subname//': ERROR: elevclass_bounds must be monotonically increasing') - end if - end do - - end subroutine check_elevclass_bounds_monotonic_increasing - - - -end module vertical_gradient_calculator_base diff --git a/src/drivers/mct/main/vertical_gradient_calculator_factory.F90 b/src/drivers/mct/main/vertical_gradient_calculator_factory.F90 deleted file mode 100644 index c6501a87e53..00000000000 --- a/src/drivers/mct/main/vertical_gradient_calculator_factory.F90 +++ /dev/null @@ -1,128 +0,0 @@ -module vertical_gradient_calculator_factory - !--------------------------------------------------------------------- - ! - ! Purpose: - ! - ! This module creates vertical gradient objects - -#include "shr_assert.h" - use shr_kind_mod, only : r8 => shr_kind_r8 - use shr_log_mod, only : errMsg => shr_log_errMsg - use vertical_gradient_calculator_2nd_order, only : vertical_gradient_calculator_2nd_order_type - use mct_mod - - implicit none - private - - public :: create_vertical_gradient_calculator_2nd_order - - ! The following routines are public just to support unit testing, and shouldn't be - ! called from production code - public :: extract_data_from_attr_vect - -contains - - !----------------------------------------------------------------------- - function create_vertical_gradient_calculator_2nd_order( & - attr_vect, fieldname, toponame, elevclass_names, elevclass_bounds) & - result(calculator) - ! - ! !DESCRIPTION: - ! Creates and returns a vertical_gradient_calculator_2nd_order_type object. - ! - ! The attribute vector is assumed to have fields named fieldname // - ! elevclass_names(1), toponame // elevclass_names(1), etc. - ! - ! !USES: - ! - ! !ARGUMENTS: - type(vertical_gradient_calculator_2nd_order_type) :: calculator ! function result - type(mct_aVect) , intent(in) :: attr_vect ! attribute vector in which we can find the data - character(len=*) , intent(in) :: fieldname ! base name of the field of interest - character(len=*) , intent(in) :: toponame ! base name of the topographic field - character(len=*) , intent(in) :: elevclass_names(:) ! strings corresponding to each elevation class - ! bounds of each elevation class; this array should have one more element than the - ! number of elevation classes, since it contains lower and upper bounds for each - ! elevation class - real(r8) , intent(in) :: elevclass_bounds(0:) - ! - ! !LOCAL VARIABLES: - integer :: nelev - real(r8), allocatable :: field(:,:) - real(r8), allocatable :: topo(:,:) - - character(len=*), parameter :: subname = 'create_vertical_gradient_calculator_2nd_order' - !----------------------------------------------------------------------- - - nelev = size(elevclass_names) - SHR_ASSERT_ALL_FL((ubound(elevclass_bounds) == (/nelev/)), __FILE__, __LINE__) - - call extract_data_from_attr_vect(attr_vect, fieldname, toponame, elevclass_names, & - field, topo) - - calculator = vertical_gradient_calculator_2nd_order_type( & - field = field, topo = topo, elevclass_bounds = elevclass_bounds) - - end function create_vertical_gradient_calculator_2nd_order - - !----------------------------------------------------------------------- - subroutine extract_data_from_attr_vect(attr_vect, fieldname, toponame, elevclass_names, & - field_extracted, topo_extracted) - ! - ! !DESCRIPTION: - ! Extract topo and data from attribute vector. - ! - ! Allocates and sets topo_extracted and data_extracted - ! - ! The attribute vector is assumed to have fields named fieldname // - ! elevclass_names(1), toponame // elevclass_names(1), etc. - ! - ! !USES: - ! - ! !ARGUMENTS: - type(mct_aVect) , intent(in) :: attr_vect ! attribute vector in which we can find the data - character(len=*) , intent(in) :: fieldname ! base name of the field of interest - character(len=*) , intent(in) :: toponame ! base name of the topographic field - character(len=*) , intent(in) :: elevclass_names(:) ! strings corresponding to each elevation class - - ! field_extracted(i,j) is point i, elevation class j; same for topo_extracted - ! these are both allocated here - real(r8), intent(out), allocatable :: field_extracted(:,:) - real(r8), intent(out), allocatable :: topo_extracted(:,:) - ! - ! !LOCAL VARIABLES: - integer :: npts - integer :: nelev - integer :: ec - character(len=:), allocatable :: fieldname_ec - character(len=:), allocatable :: toponame_ec - - ! The following temporary array is needed because mct wants pointers - real(r8), pointer :: temp(:) - - character(len=*), parameter :: subname = 'extract_data_from_attr_vect' - !----------------------------------------------------------------------- - - nelev = size(elevclass_names) - npts = mct_aVect_lsize(attr_vect) - - allocate(field_extracted(npts, nelev)) - allocate(topo_extracted(npts, nelev)) - allocate(temp(npts)) - - do ec = 1, nelev - fieldname_ec = trim(fieldname) // trim(elevclass_names(ec)) - call mct_aVect_exportRattr(attr_vect, fieldname_ec, temp) - field_extracted(:,ec) = temp(:) - - toponame_ec = trim(toponame) // trim(elevclass_names(ec)) - call mct_aVect_exportRattr(attr_vect, toponame_ec, temp) - topo_extracted(:,ec) = temp(:) - end do - - deallocate(temp) - - end subroutine extract_data_from_attr_vect - - -end module vertical_gradient_calculator_factory diff --git a/src/drivers/mct/shr/glc_elevclass_mod.F90 b/src/drivers/mct/shr/glc_elevclass_mod.F90 index 8ae9ac652a5..7ab9f5fdc8f 100644 --- a/src/drivers/mct/shr/glc_elevclass_mod.F90 +++ b/src/drivers/mct/shr/glc_elevclass_mod.F90 @@ -46,6 +46,9 @@ module glc_elevclass_mod integer, parameter, public :: GLC_ELEVCLASS_ERR_TOO_LOW = 2 ! err_code indicating topo below lowest elevation class integer, parameter, public :: GLC_ELEVCLASS_ERR_TOO_HIGH = 3 ! err_code indicating topo above highest elevation class + ! String length for glc elevation classes represented as strings + integer, parameter, public :: GLC_ELEVCLASS_STRLEN = 2 + !-------------------------------------------------------------------------- ! Private data !-------------------------------------------------------------------------- @@ -324,38 +327,60 @@ function glc_elevclass_as_string(elevation_class) result(ec_string) ! !USES: ! ! !ARGUMENTS: - character(len=2) :: ec_string ! function result + character(len=GLC_ELEVCLASS_STRLEN) :: ec_string ! function result integer, intent(in) :: elevation_class ! ! !LOCAL VARIABLES: + character(len=16) :: format_string character(len=*), parameter :: subname = 'glc_elevclass_as_string' !----------------------------------------------------------------------- - write(ec_string,'(i2.2)') elevation_class + ! e.g., for GLC_ELEVCLASS_STRLEN = 2, format_string will be '(i2.2)' + write(format_string,'(a,i0,a,i0,a)') '(i', GLC_ELEVCLASS_STRLEN, '.', GLC_ELEVCLASS_STRLEN, ')' + + write(ec_string,trim(format_string)) elevation_class end function glc_elevclass_as_string !----------------------------------------------------------------------- - function glc_all_elevclass_strings() result(ec_strings) + function glc_all_elevclass_strings(include_zero) result(ec_strings) ! ! !DESCRIPTION: ! Returns an array of strings corresponding to all elevation classes from 1 to glc_nec ! + ! If include_zero is present and true, then includes elevation class 0 - so goes from + ! 0 to glc_nec + ! ! These strings can be used as suffixes for fields in MCT attribute vectors. ! ! !USES: ! ! !ARGUMENTS: - character(len=2), allocatable :: ec_strings(:) ! function result + character(len=GLC_ELEVCLASS_STRLEN), allocatable :: ec_strings(:) ! function result + logical, intent(in), optional :: include_zero ! if present and true, include elevation class 0 (default is false) ! ! !LOCAL VARIABLES: + logical :: l_include_zero ! local version of optional include_zero argument + integer :: lower_bound integer :: i character(len=*), parameter :: subname = 'glc_all_elevclass_strings' !----------------------------------------------------------------------- - allocate(ec_strings(1:glc_nec)) - do i = 1, glc_nec + if (present(include_zero)) then + l_include_zero = include_zero + else + l_include_zero = .false. + end if + + if (l_include_zero) then + lower_bound = 0 + else + lower_bound = 1 + end if + + allocate(ec_strings(lower_bound:glc_nec)) + do i = lower_bound, glc_nec ec_strings(i) = glc_elevclass_as_string(i) end do diff --git a/src/drivers/mct/shr/seq_infodata_mod.F90 b/src/drivers/mct/shr/seq_infodata_mod.F90 index d62014fa87f..01c6b67689f 100644 --- a/src/drivers/mct/shr/seq_infodata_mod.F90 +++ b/src/drivers/mct/shr/seq_infodata_mod.F90 @@ -124,6 +124,7 @@ MODULE seq_infodata_mod logical :: flux_albav ! T => no diurnal cycle in ocn albedos logical :: flux_diurnal ! T => diurnal cycle in atm/ocn fluxes real(SHR_KIND_R8) :: gust_fac ! wind gustiness factor + character(SHR_KIND_CL) :: glc_renormalize_smb ! Whether to renormalize smb sent from lnd -> glc real(SHR_KIND_R8) :: wall_time_limit ! force stop time limit (hours) character(SHR_KIND_CS) :: force_stop_at ! when to force a stop (month, day, etc) character(SHR_KIND_CL) :: atm_gnam ! atm grid @@ -205,6 +206,7 @@ MODULE seq_infodata_mod logical :: glcocn_present ! does glc have ocean runoff on logical :: glcice_present ! does glc have iceberg coupling on logical :: glc_prognostic ! does component model need input data from driver + logical :: glc_coupled_fluxes ! does glc send fluxes to other components (only relevant if glc_present is .true.) logical :: wav_present ! does component model exist logical :: wav_prognostic ! does component model need input data from driver logical :: esp_present ! does component model exist @@ -355,6 +357,7 @@ SUBROUTINE seq_infodata_Init( infodata, nmlfile, ID, pioid) logical :: flux_albav ! T => no diurnal cycle in ocn albedos logical :: flux_diurnal ! T => diurnal cycle in atm/ocn fluxes real(SHR_KIND_R8) :: gust_fac ! wind gustiness factor + character(SHR_KIND_CL) :: glc_renormalize_smb ! Whether to renormalize smb sent from lnd -> glc real(SHR_KIND_R8) :: wall_time_limit ! force stop time limit (hours) character(SHR_KIND_CS) :: force_stop_at ! when to force a stop (month, day, etc) character(SHR_KIND_CL) :: atm_gnam ! atm grid @@ -424,7 +427,7 @@ SUBROUTINE seq_infodata_Init( infodata, nmlfile, ID, pioid) orb_iyear, orb_obliq, orb_eccen, orb_mvelp, & wv_sat_scheme, wv_sat_transition_start, & wv_sat_use_tables, wv_sat_table_spacing, & - tfreeze_option, & + tfreeze_option, glc_renormalize_smb, & ice_gnam, rof_gnam, glc_gnam, wav_gnam, & atm_gnam, lnd_gnam, ocn_gnam, cpl_decomp, & shr_map_dopole, vect_map, aoflux_grid, do_histinit, & @@ -495,6 +498,7 @@ SUBROUTINE seq_infodata_Init( infodata, nmlfile, ID, pioid) flux_albav = .false. flux_diurnal = .false. gust_fac = huge(1.0_SHR_KIND_R8) + glc_renormalize_smb = 'on_if_glc_coupled_fluxes' wall_time_limit = -1.0 force_stop_at = 'month' atm_gnam = 'undefined' @@ -601,6 +605,7 @@ SUBROUTINE seq_infodata_Init( infodata, nmlfile, ID, pioid) infodata%flux_albav = flux_albav infodata%flux_diurnal = flux_diurnal infodata%gust_fac = gust_fac + infodata%glc_renormalize_smb = glc_renormalize_smb infodata%wall_time_limit = wall_time_limit infodata%force_stop_at = force_stop_at infodata%atm_gnam = atm_gnam @@ -680,6 +685,11 @@ SUBROUTINE seq_infodata_Init( infodata, nmlfile, ID, pioid) infodata%ocnrof_prognostic = .false. infodata%ice_prognostic = .false. infodata%glc_prognostic = .false. + ! It's safest to assume glc_coupled_fluxes = .true. if it's not set elsewhere, + ! because this is needed for conservation in some cases. Note that it is ignored + ! if glc_present is .false., so it's okay to just start out assuming it's .true. + ! in all cases. + infodata%glc_coupled_fluxes = .true. infodata%wav_prognostic = .false. infodata%iceberg_prognostic = .false. infodata%esp_prognostic = .false. @@ -907,6 +917,7 @@ SUBROUTINE seq_infodata_GetData_explicit( infodata, cime_model, case_name, case_ atm_present, atm_prognostic, lnd_present, lnd_prognostic, rof_prognostic, & rof_present, ocn_present, ocn_prognostic, ocnrof_prognostic, & ice_present, ice_prognostic, glc_present, glc_prognostic, & + glc_coupled_fluxes, & flood_present, wav_present, wav_prognostic, rofice_present, & glclnd_present, glcocn_present, glcice_present, iceberg_prognostic,& esp_present, esp_prognostic, & @@ -927,7 +938,7 @@ SUBROUTINE seq_infodata_GetData_explicit( infodata, cime_model, case_name, case_ cpl_cdf64, orb_iyear, orb_iyear_align, orb_mode, orb_mvelp, & orb_eccen, orb_obliqr, orb_lambm0, orb_mvelpp, wv_sat_scheme, & wv_sat_transition_start, wv_sat_use_tables, wv_sat_table_spacing, & - tfreeze_option, & + tfreeze_option, glc_renormalize_smb, & glc_phase, rof_phase, atm_phase, lnd_phase, ocn_phase, ice_phase, & wav_phase, esp_phase, wav_nx, wav_ny, atm_nx, atm_ny, & lnd_nx, lnd_ny, rof_nx, rof_ny, ice_nx, ice_ny, ocn_nx, ocn_ny, & @@ -986,6 +997,7 @@ SUBROUTINE seq_infodata_GetData_explicit( infodata, cime_model, case_name, case_ logical, optional, intent(OUT) :: flux_albav ! T => no diurnal cycle in ocn albedos logical, optional, intent(OUT) :: flux_diurnal ! T => diurnal cycle in atm/ocn flux real(SHR_KIND_R8), optional, intent(OUT) :: gust_fac ! wind gustiness factor + character(len=*), optional, intent(OUT) :: glc_renormalize_smb ! Whether to renormalize smb sent from lnd -> glc real(SHR_KIND_R8), optional, intent(OUT) :: wall_time_limit ! force stop wall time (hours) character(len=*), optional, intent(OUT) :: force_stop_at ! force stop at next (month, day, etc) character(len=*), optional, intent(OUT) :: atm_gnam ! atm grid @@ -1064,6 +1076,7 @@ SUBROUTINE seq_infodata_GetData_explicit( infodata, cime_model, case_name, case_ logical, optional, intent(OUT) :: glcocn_present logical, optional, intent(OUT) :: glcice_present logical, optional, intent(OUT) :: glc_prognostic + logical, optional, intent(OUT) :: glc_coupled_fluxes logical, optional, intent(OUT) :: wav_present logical, optional, intent(OUT) :: wav_prognostic logical, optional, intent(OUT) :: esp_present @@ -1156,6 +1169,7 @@ SUBROUTINE seq_infodata_GetData_explicit( infodata, cime_model, case_name, case_ if ( present(flux_albav) ) flux_albav = infodata%flux_albav if ( present(flux_diurnal) ) flux_diurnal = infodata%flux_diurnal if ( present(gust_fac) ) gust_fac = infodata%gust_fac + if ( present(glc_renormalize_smb)) glc_renormalize_smb = infodata%glc_renormalize_smb if ( present(wall_time_limit)) wall_time_limit= infodata%wall_time_limit if ( present(force_stop_at) ) force_stop_at = infodata%force_stop_at if ( present(atm_gnam) ) atm_gnam = infodata%atm_gnam @@ -1234,6 +1248,7 @@ SUBROUTINE seq_infodata_GetData_explicit( infodata, cime_model, case_name, case_ if ( present(glcocn_present) ) glcocn_present = infodata%glcocn_present if ( present(glcice_present) ) glcice_present = infodata%glcice_present if ( present(glc_prognostic) ) glc_prognostic = infodata%glc_prognostic + if ( present(glc_coupled_fluxes)) glc_coupled_fluxes = infodata%glc_coupled_fluxes if ( present(wav_present) ) wav_present = infodata%wav_present if ( present(wav_prognostic) ) wav_prognostic = infodata%wav_prognostic if ( present(esp_present) ) esp_present = infodata%esp_present @@ -1467,6 +1482,7 @@ SUBROUTINE seq_infodata_PutData_explicit( infodata, cime_model, case_name, case_ atm_present, atm_prognostic, lnd_present, lnd_prognostic, rof_prognostic, & rof_present, ocn_present, ocn_prognostic, ocnrof_prognostic, & ice_present, ice_prognostic, glc_present, glc_prognostic, & + glc_coupled_fluxes, & flood_present, wav_present, wav_prognostic, rofice_present, & glclnd_present, glcocn_present, glcice_present, iceberg_prognostic,& esp_present, esp_prognostic, & @@ -1487,7 +1503,7 @@ SUBROUTINE seq_infodata_PutData_explicit( infodata, cime_model, case_name, case_ cpl_cdf64, orb_iyear, orb_iyear_align, orb_mode, orb_mvelp, & orb_eccen, orb_obliqr, orb_lambm0, orb_mvelpp, wv_sat_scheme, & wv_sat_transition_start, wv_sat_use_tables, wv_sat_table_spacing, & - tfreeze_option, & + tfreeze_option, glc_renormalize_smb, & glc_phase, rof_phase, atm_phase, lnd_phase, ocn_phase, ice_phase, & wav_phase, esp_phase, wav_nx, wav_ny, atm_nx, atm_ny, & lnd_nx, lnd_ny, rof_nx, rof_ny, ice_nx, ice_ny, ocn_nx, ocn_ny, & @@ -1546,6 +1562,7 @@ SUBROUTINE seq_infodata_PutData_explicit( infodata, cime_model, case_name, case_ logical, optional, intent(IN) :: flux_albav ! T => no diurnal cycle in ocn albedos logical, optional, intent(IN) :: flux_diurnal ! T => diurnal cycle in atm/ocn flux real(SHR_KIND_R8), optional, intent(IN) :: gust_fac ! wind gustiness factor + character(len=*), optional, intent(IN) :: glc_renormalize_smb ! Whether to renormalize smb sent from lnd -> glc real(SHR_KIND_R8), optional, intent(IN) :: wall_time_limit ! force stop wall time (hours) character(len=*), optional, intent(IN) :: force_stop_at ! force a stop at next (month, day, etc) character(len=*), optional, intent(IN) :: atm_gnam ! atm grid @@ -1624,6 +1641,7 @@ SUBROUTINE seq_infodata_PutData_explicit( infodata, cime_model, case_name, case_ logical, optional, intent(IN) :: glcocn_present logical, optional, intent(IN) :: glcice_present logical, optional, intent(IN) :: glc_prognostic + logical, optional, intent(IN) :: glc_coupled_fluxes logical, optional, intent(IN) :: wav_present logical, optional, intent(IN) :: wav_prognostic logical, optional, intent(IN) :: esp_present @@ -1715,6 +1733,7 @@ SUBROUTINE seq_infodata_PutData_explicit( infodata, cime_model, case_name, case_ if ( present(flux_albav) ) infodata%flux_albav = flux_albav if ( present(flux_diurnal) ) infodata%flux_diurnal = flux_diurnal if ( present(gust_fac) ) infodata%gust_fac = gust_fac + if ( present(glc_renormalize_smb)) infodata%glc_renormalize_smb = glc_renormalize_smb if ( present(wall_time_limit)) infodata%wall_time_limit= wall_time_limit if ( present(force_stop_at) ) infodata%force_stop_at = force_stop_at if ( present(atm_gnam) ) infodata%atm_gnam = atm_gnam @@ -1793,6 +1812,7 @@ SUBROUTINE seq_infodata_PutData_explicit( infodata, cime_model, case_name, case_ if ( present(glcocn_present) ) infodata%glcocn_present = glcocn_present if ( present(glcice_present) ) infodata%glcice_present = glcice_present if ( present(glc_prognostic) ) infodata%glc_prognostic = glc_prognostic + if ( present(glc_coupled_fluxes)) infodata%glc_coupled_fluxes = glc_coupled_fluxes if ( present(wav_present) ) infodata%wav_present = wav_present if ( present(wav_prognostic) ) infodata%wav_prognostic = wav_prognostic if ( present(esp_present) ) infodata%esp_present = esp_present @@ -2133,6 +2153,7 @@ subroutine seq_infodata_bcast(infodata,mpicom) call shr_mpi_bcast(infodata%flux_albav, mpicom) call shr_mpi_bcast(infodata%flux_diurnal, mpicom) call shr_mpi_bcast(infodata%gust_fac, mpicom) + call shr_mpi_bcast(infodata%glc_renormalize_smb, mpicom) call shr_mpi_bcast(infodata%wall_time_limit, mpicom) call shr_mpi_bcast(infodata%force_stop_at, mpicom) call shr_mpi_bcast(infodata%atm_gnam, mpicom) @@ -2211,6 +2232,7 @@ subroutine seq_infodata_bcast(infodata,mpicom) call shr_mpi_bcast(infodata%glcocn_present, mpicom) call shr_mpi_bcast(infodata%glcice_present, mpicom) call shr_mpi_bcast(infodata%glc_prognostic, mpicom) + call shr_mpi_bcast(infodata%glc_coupled_fluxes, mpicom) call shr_mpi_bcast(infodata%wav_present, mpicom) call shr_mpi_bcast(infodata%wav_prognostic, mpicom) call shr_mpi_bcast(infodata%esp_present, mpicom) @@ -2498,6 +2520,7 @@ subroutine seq_infodata_Exchange(infodata,ID,type) call shr_mpi_bcast(infodata%glcocn_present, mpicom, pebcast=cmppe) call shr_mpi_bcast(infodata%glcice_present, mpicom, pebcast=cmppe) call shr_mpi_bcast(infodata%glc_prognostic, mpicom, pebcast=cmppe) + call shr_mpi_bcast(infodata%glc_coupled_fluxes, mpicom, pebcast=cmppe) call shr_mpi_bcast(infodata%glc_nx, mpicom, pebcast=cmppe) call shr_mpi_bcast(infodata%glc_ny, mpicom, pebcast=cmppe) ! dead_comps is true if it's ever set to true @@ -2543,6 +2566,7 @@ subroutine seq_infodata_Exchange(infodata,ID,type) call shr_mpi_bcast(infodata%glcocn_present, mpicom, pebcast=cplpe) call shr_mpi_bcast(infodata%glcice_present, mpicom, pebcast=cplpe) call shr_mpi_bcast(infodata%glc_prognostic, mpicom, pebcast=cplpe) + call shr_mpi_bcast(infodata%glc_coupled_fluxes, mpicom, pebcast=cplpe) call shr_mpi_bcast(infodata%wav_present, mpicom, pebcast=cplpe) call shr_mpi_bcast(infodata%wav_prognostic, mpicom, pebcast=cplpe) call shr_mpi_bcast(infodata%esp_present, mpicom, pebcast=cplpe) @@ -2797,6 +2821,7 @@ SUBROUTINE seq_infodata_print( infodata ) write(logunit,F0L) subname,'flux_albav = ', infodata%flux_albav write(logunit,F0L) subname,'flux_diurnal = ', infodata%flux_diurnal write(logunit,F0R) subname,'gust_fac = ', infodata%gust_fac + write(logunit,F0A) subname,'glc_renormalize_smb = ', trim(infodata%glc_renormalize_smb) write(logunit,F0R) subname,'wall_time_limit = ', infodata%wall_time_limit write(logunit,F0A) subname,'force_stop_at = ', trim(infodata%force_stop_at) write(logunit,F0A) subname,'atm_gridname = ', trim(infodata%atm_gnam) @@ -2879,6 +2904,7 @@ SUBROUTINE seq_infodata_print( infodata ) write(logunit,F0L) subname,'glcocn_present = ', infodata%glcocn_present write(logunit,F0L) subname,'glcice_present = ', infodata%glcice_present write(logunit,F0L) subname,'glc_prognostic = ', infodata%glc_prognostic + write(logunit,F0L) subname,'glc_coupled_fluxes = ', infodata%glc_coupled_fluxes write(logunit,F0L) subname,'wav_present = ', infodata%wav_present write(logunit,F0L) subname,'wav_prognostic = ', infodata%wav_prognostic write(logunit,F0L) subname,'esp_present = ', infodata%esp_present diff --git a/src/drivers/mct/unit_test/CMakeLists.txt b/src/drivers/mct/unit_test/CMakeLists.txt index a747f9e23a2..83010758c29 100644 --- a/src/drivers/mct/unit_test/CMakeLists.txt +++ b/src/drivers/mct/unit_test/CMakeLists.txt @@ -60,7 +60,5 @@ list(APPEND DRV_UNIT_TEST_LIBS ${NETCDF_LIBRARIES}) add_subdirectory(avect_wrapper_test) add_subdirectory(seq_map_test) add_subdirectory(glc_elevclass_test) -add_subdirectory(vertical_gradient_calculator_test) -add_subdirectory(map_lnd2glc_test) add_subdirectory(map_glc2lnd_test) add_subdirectory(map_lnd2rof_irrig_test) diff --git a/src/drivers/mct/unit_test/glc_elevclass_test/test_glc_elevclass.pf b/src/drivers/mct/unit_test/glc_elevclass_test/test_glc_elevclass.pf index 8c33b80599c..f0f8ebed22a 100644 --- a/src/drivers/mct/unit_test/glc_elevclass_test/test_glc_elevclass.pf +++ b/src/drivers/mct/unit_test/glc_elevclass_test/test_glc_elevclass.pf @@ -223,7 +223,7 @@ contains @Test subroutine test_glc_elevclass_as_string_0(this) class(TestGLCElevclass), intent(inout) :: this - character(len=16) :: str + character(len=GLC_ELEVCLASS_STRLEN) :: str str = glc_elevclass_as_string(0) @assertEqual('00', trim(str)) @@ -232,7 +232,7 @@ contains @Test subroutine test_glc_elevclass_as_string_1digit(this) class(TestGLCElevclass), intent(inout) :: this - character(len=16) :: str + character(len=GLC_ELEVCLASS_STRLEN) :: str str = glc_elevclass_as_string(2) @assertEqual('02', trim(str)) @@ -241,7 +241,7 @@ contains @Test subroutine test_glc_elevclass_as_string_2digits(this) class(TestGLCElevclass), intent(inout) :: this - character(len=16) :: str + character(len=GLC_ELEVCLASS_STRLEN) :: str str = glc_elevclass_as_string(12) @assertEqual('12', trim(str)) @@ -254,15 +254,33 @@ contains @Test subroutine test_glc_all_elevclass_strings(this) class(TestGLCElevclass), intent(inout) :: this - character(len=16) :: elevclass_strings(3) + character(len=GLC_ELEVCLASS_STRLEN), allocatable :: elevclass_strings(:) call glc_elevclass_init(3) elevclass_strings = glc_all_elevclass_strings() + @assertEqual(3, size(elevclass_strings)) ! There doesn't seem to be an assertEqual method for an array of strings @assertEqual('01', elevclass_strings(1)) @assertEqual('02', elevclass_strings(2)) @assertEqual('03', elevclass_strings(3)) end subroutine test_glc_all_elevclass_strings + + @Test + subroutine test_glc_all_elevclass_strings_include_zero(this) + class(TestGLCElevclass), intent(inout) :: this + character(len=GLC_ELEVCLASS_STRLEN), allocatable :: elevclass_strings(:) + + call glc_elevclass_init(3) + elevclass_strings = glc_all_elevclass_strings(include_zero=.true.) + + @assertEqual(4, size(elevclass_strings)) + ! There doesn't seem to be an assertEqual method for an array of strings + @assertEqual('00', elevclass_strings(1)) + @assertEqual('01', elevclass_strings(2)) + @assertEqual('02', elevclass_strings(3)) + @assertEqual('03', elevclass_strings(4)) + end subroutine test_glc_all_elevclass_strings_include_zero + end module test_glc_elevclass diff --git a/src/drivers/mct/unit_test/map_lnd2glc_test/CMakeLists.txt b/src/drivers/mct/unit_test/map_lnd2glc_test/CMakeLists.txt deleted file mode 100644 index 459660cbb47..00000000000 --- a/src/drivers/mct/unit_test/map_lnd2glc_test/CMakeLists.txt +++ /dev/null @@ -1,8 +0,0 @@ -set (pfunit_sources - test_map_lnd2glc.pf - ) - -create_pFUnit_test(map_lnd2glc map_lnd2glc_exe - "${pfunit_sources}" "") - -target_link_libraries(map_lnd2glc_exe ${DRV_UNIT_TEST_LIBS}) diff --git a/src/drivers/mct/unit_test/map_lnd2glc_test/test_map_lnd2glc.pf b/src/drivers/mct/unit_test/map_lnd2glc_test/test_map_lnd2glc.pf deleted file mode 100644 index f677b6b2ec6..00000000000 --- a/src/drivers/mct/unit_test/map_lnd2glc_test/test_map_lnd2glc.pf +++ /dev/null @@ -1,773 +0,0 @@ -module test_map_lnd2glc - - ! Tests of map_lnd2glc_mod - - use pfunit_mod - use map_lnd2glc_mod - use glc_elevclass_mod, only : glc_elevclass_init, glc_elevclass_clean - use mct_mod, only : mct_aVect, mct_aVect_clean - use seq_map_type_mod, only : seq_map - use mct_wrapper_mod, only : mct_init, mct_clean - use avect_wrapper_mod - use simple_map_mod - use create_mapper_mod - use vertical_gradient_calculator_base, only : vertical_gradient_calculator_base_type - use vertical_gradient_calculator_specified, only : & - vertical_gradient_calculator_specified_type, vgc_specified_ec_times_ptSquared - use shr_kind_mod, only : r8 => shr_kind_r8 - - implicit none - - real(r8), parameter :: tol = 1.e-11_r8 - - integer, parameter :: n_elev_classes = 3 - - ! Assume 3 elevation classes, with boundaries of: - ! (1) 0 - 100 m - ! (2) 100 - 200 m - ! (3) 200 - 300 m - real(r8), parameter :: elev_class_boundaries(0:n_elev_classes) = & - [0._r8, 100._r8, 200._r8, 300._r8] - - - ! This type holds the data in a single land grid cell - type :: lnd_data_type - ! Index 0 is bare land - real(r8) :: topo(0:n_elev_classes) - real(r8) :: data(0:n_elev_classes) - end type lnd_data_type - - @TestCase - type, extends(TestCase) :: TestMapLnd2glc - type(seq_map) :: mapper - type(mct_aVect) :: data_l ! data on the LND (source) grid - type(mct_aVect) :: landfrac_l ! landfrac on the LND (source) grid - type(mct_aVect) :: data_g ! data on the GLC (destination) grid - type(mct_aVect) :: g2x ! data sent from glc -> cpl - contains - procedure :: setUp - procedure :: tearDown - procedure :: setup_inputs - procedure :: run_map_lnd2glc - procedure :: verify_data_g - end type TestMapLnd2glc - -contains - - ! ======================================================================== - ! Utility routines - ! ======================================================================== - - subroutine setUp(this) - class(TestMapLnd2glc), intent(inout) :: this - - call mct_init() - - end subroutine setUp - - subroutine tearDown(this) - class(TestMapLnd2glc), intent(inout) :: this - - call clean_mapper(this%mapper) - call mct_aVect_clean(this%data_l) - call mct_aVect_clean(this%landfrac_l) - call mct_aVect_clean(this%data_g) - call mct_aVect_clean(this%g2x) - call glc_elevclass_clean() - call mct_clean() - end subroutine tearDown - - subroutine setup_inputs(this, frac_glc, topo_glc, lnd_data, my_map, landfrac) - ! This utility function sets up inputs that are needed for the map_lnd2glc call - - class(TestMapLnd2glc), intent(inout) :: this - real(r8), intent(in) :: frac_glc(:) ! ice fraction in each glc cell - real(r8), intent(in) :: topo_glc(:) ! ice topographic height in each glc cell - type(lnd_data_type), intent(in) :: lnd_data(:) ! land data in each grid cell - type(simple_map_type), intent(in) :: my_map ! mapping information from land to glc - - ! If provided, this gives the landfrac of each land cell. If absent, landfrac is - ! assumed to be 1 for every land cell - real(r8), intent(in), optional :: landfrac(:) - - integer :: npts_glc - integer :: npts_lnd - integer :: lnd_index - real(r8), allocatable :: l_landfrac(:) ! local version of landfrac - - ! ------------------------------------------------------------------------ - ! Do some initial error-checking to make sure this routine is being called properly - ! ------------------------------------------------------------------------ - npts_glc = size(frac_glc) - @assertEqual(npts_glc, size(topo_glc)) - - npts_lnd = size(lnd_data) - if (present(landfrac)) then - @assertEqual(npts_lnd, size(landfrac)) - end if - - ! ------------------------------------------------------------------------ - ! Set local version of optional arguments - ! ------------------------------------------------------------------------ - - allocate(l_landfrac(npts_lnd)) - if (present(landfrac)) then - l_landfrac = landfrac - else - l_landfrac(:) = 1._r8 - end if - - ! ------------------------------------------------------------------------ - ! Setup - ! ------------------------------------------------------------------------ - - call glc_elevclass_init(n_elev_classes, elev_class_boundaries) - - ! The following assumes that n_elev_classes is 3: - call create_aVect_with_data_rows_are_fields(this%data_l, & - attr_tags = ['data00 ', 'data01 ', 'data02 ', 'data03 ', & - 'Sl_topo00', 'Sl_topo01', 'Sl_topo02', 'Sl_topo03'], & - data = reshape([(lnd_data(lnd_index)%data, lnd_data(lnd_index)%topo, & - lnd_index = 1, npts_lnd)], & - [8,npts_lnd])) - - call create_aVect_with_data_rows_are_fields(this%landfrac_l, & - attr_tags = ['lfrac'], & - data = reshape(l_landfrac, [1, npts_lnd])) - - call create_aVect_with_data_rows_are_points(this%g2x, & - attr_tags = ['Sg_ice_covered', 'Sg_topo '], & - data = reshape([frac_glc, topo_glc], [npts_glc, 2])) - - call create_aVect_without_data(this%data_g, attr_tags = ['data'], lsize = npts_glc) - - call create_mapper(this%mapper, my_map) - - end subroutine setup_inputs - - subroutine run_map_lnd2glc(this, gradient_calculator, fieldname) - ! This utility function wraps the call to the map_lnd2glc routine - - class(TestMapLnd2glc), intent(inout) :: this - class(vertical_gradient_calculator_base_type), intent(inout) :: gradient_calculator - - ! Name of field to map. If not provided, uses 'data'. (This argument is available to - ! test particular cases, such as having trailing blanks in the fieldname; for most - ! tests, it can be omitted.) - character(len=*), intent(in), optional :: fieldname - - character(len=:), allocatable :: l_fieldname ! local version of fieldname - - l_fieldname = 'data' - if (present(fieldname)) then - l_fieldname = fieldname - end if - - call map_lnd2glc(this%data_l, this%landfrac_l, this%g2x, l_fieldname, gradient_calculator, & - this%mapper, this%data_g) - end subroutine run_map_lnd2glc - - subroutine verify_data_g(this, expected_data_glc, message) - ! Verify that the remapped data (in this%data_g) matches expected_data_glc - - class(TestMapLnd2glc), intent(in) :: this - real(r8), intent(in) :: expected_data_glc(:) ! expected outputs in each glc cell - character(len=*), intent(in) :: message ! message printed if assertion fails - - real(r8), allocatable :: actual_data_glc(:) - - actual_data_glc = aVect_exportRattr(this%data_g, 'data') - @assertEqual(expected_data_glc, actual_data_glc, message=message, tolerance=tol) - - end subroutine verify_data_g - - subroutine create_data_for_EC2_gradient2(my_map, gradient_calculator, lnd_data, & - frac_glc, topo_glc, expected_data_glc) - ! Helper routine to set up all of the data needed to run a test with a single source - ! grid cell, a single destination point in elevation class 2, and a gradient of 2*EC. - ! - ! This can be used for one-offs off of this setup that is not too simple, but not too - ! complex. - - type(simple_map_type), intent(out) :: my_map - type(vertical_gradient_calculator_specified_type), intent(out) :: gradient_calculator - type(lnd_data_type), intent(out) :: lnd_data - real(r8), intent(out) :: frac_glc(1) - real(r8), intent(out) :: topo_glc(1) - real(r8), intent(out) :: expected_data_glc(1) - - my_map = create_simple_map_with_one_source(ndest = 1) - - gradient_calculator = & - vgc_specified_ec_times_ptSquared(num_points = 1, & - nelev = n_elev_classes, gradient = 2._r8) - - ! data in elev class: 0 1 2 3 - lnd_data%topo(:) = [25._r8, 50._r8, 150._r8, 250._r8] - lnd_data%data(:) = [10._r8, 11._r8, 12._r8, 13._r8] - - frac_glc(1) = 1._r8 - topo_glc(1) = 125._r8 - - ! For expected data, note that we multiply the topo difference by 4: - ! gradient * EC = 2 * 2 = 4 - expected_data_glc(1) = lnd_data%data(2) + 4._r8*(125._r8 - lnd_data%topo(2)) - - end subroutine create_data_for_EC2_gradient2 - - ! ======================================================================== - ! Actual tests - ! ======================================================================== - - @Test - subroutine test_mapLnd2glc_with_EC0_gradient0(this) - ! Do a test of the map_lnd2glc routine with only an elevation class 0 destination - ! point, and a gradient of 0. This tests the very basic operation of map_lnd2glc. - class(TestMapLnd2glc), intent(inout) :: this - - type(lnd_data_type) :: lnd_data - - integer, parameter :: npts_glc = 1 - real(r8), parameter :: frac_glc(npts_glc) = [0._r8] - real(r8), parameter :: topo_glc(npts_glc) = [75._r8] - real(r8) :: expected_data_glc(npts_glc) - - type(simple_map_type) :: my_map - type(vertical_gradient_calculator_specified_type) :: gradient_calculator - - ! ------------------------------------------------------------------------ - ! Setup - ! ------------------------------------------------------------------------ - - my_map = create_simple_map_with_one_source(ndest = npts_glc) - - gradient_calculator = & - vgc_specified_ec_times_ptSquared(num_points = 1, & - nelev = n_elev_classes, gradient = 0._r8) - - ! data in elev class: 0 1 2 3 - lnd_data%topo(:) = [25._r8, 50._r8, 150._r8, 250._r8] - lnd_data%data(:) = [10._r8, 11._r8, 12._r8, 13._r8] - - call this%setup_inputs(frac_glc, topo_glc, [lnd_data], my_map) - - ! ------------------------------------------------------------------------ - ! Exercise - ! ------------------------------------------------------------------------ - - call this%run_map_lnd2glc(gradient_calculator) - - ! ------------------------------------------------------------------------ - ! Verify - ! ------------------------------------------------------------------------ - - expected_data_glc(1) = lnd_data%data(0) - - call this%verify_data_g(expected_data_glc, message = "test_mapLnd2glc_with_EC0_gradient0") - - end subroutine test_mapLnd2glc_with_EC0_gradient0 - - @Test - subroutine test_mapLnd2glc_with_EC2_gradient0(this) - ! Do a test of the map_lnd2glc routine with only an elevation class 2 destination - ! point, and a gradient of 0. - class(TestMapLnd2glc), intent(inout) :: this - - type(lnd_data_type) :: lnd_data - - integer, parameter :: npts_glc = 1 - real(r8), parameter :: frac_glc(npts_glc) = [1._r8] - real(r8), parameter :: topo_glc(npts_glc) = [125._r8] - real(r8) :: expected_data_glc(npts_glc) - - type(simple_map_type) :: my_map - type(vertical_gradient_calculator_specified_type) :: gradient_calculator - - ! ------------------------------------------------------------------------ - ! Setup - ! ------------------------------------------------------------------------ - - my_map = create_simple_map_with_one_source(ndest = npts_glc) - - gradient_calculator = & - vgc_specified_ec_times_ptSquared(num_points = 1, & - nelev = n_elev_classes, gradient = 0._r8) - - ! data in elev class: 0 1 2 3 - lnd_data%topo(:) = [25._r8, 50._r8, 150._r8, 250._r8] - lnd_data%data(:) = [10._r8, 11._r8, 12._r8, 13._r8] - - call this%setup_inputs(frac_glc, topo_glc, [lnd_data], my_map) - - ! ------------------------------------------------------------------------ - ! Exercise - ! ------------------------------------------------------------------------ - - call this%run_map_lnd2glc(gradient_calculator) - - ! ------------------------------------------------------------------------ - ! Verify - ! ------------------------------------------------------------------------ - - expected_data_glc(1) = lnd_data%data(2) - - call this%verify_data_g(expected_data_glc, message = "test_mapLnd2glc_with_EC2_gradient0") - - end subroutine test_mapLnd2glc_with_EC2_gradient0 - - @Test - subroutine test_mapLnd2glc_with_multipleDest_gradient0(this) - ! Do a test of the map_lnd2glc routine with multiple destination points, and a - ! gradient of 0 - i.e., not trying to correct for the vertical gradient. This tests to - ! make sure that each destination (GLC) grid cell gets data from the appropriate - ! source (LND) elevation class. - class(TestMapLnd2glc), intent(inout) :: this - - type(lnd_data_type) :: lnd_data - - ! On the glc grid, there are 4 grid cells: one bare land and one in each elevation class - integer, parameter :: npts_glc = 4 - real(r8), parameter :: frac_glc(npts_glc) = [1._r8, 1._r8, 1._r8, 0._r8] - real(r8), parameter :: topo_glc(npts_glc) = [225._r8, 125._r8, 25._r8, 75._r8] - real(r8) :: expected_data_glc(npts_glc) - - type(simple_map_type) :: my_map - type(vertical_gradient_calculator_specified_type) :: gradient_calculator - - ! ------------------------------------------------------------------------ - ! Setup - ! ------------------------------------------------------------------------ - - my_map = create_simple_map_with_one_source(ndest = npts_glc) - - gradient_calculator = & - vgc_specified_ec_times_ptSquared(num_points = 1, & - nelev = n_elev_classes, gradient = 0._r8) - - ! data in elev class: 0 1 2 3 - lnd_data%topo(:) = [25._r8, 50._r8, 150._r8, 250._r8] - lnd_data%data(:) = [10._r8, 11._r8, 12._r8, 13._r8] - - call this%setup_inputs(frac_glc, topo_glc, [lnd_data], my_map) - - ! ------------------------------------------------------------------------ - ! Exercise - ! ------------------------------------------------------------------------ - - call this%run_map_lnd2glc(gradient_calculator) - - ! ------------------------------------------------------------------------ - ! Verify - ! ------------------------------------------------------------------------ - - expected_data_glc(:) = [lnd_data%data(3), lnd_data%data(2), lnd_data%data(1), lnd_data%data(0)] - - call this%verify_data_g(expected_data_glc, message = "test_mapLnd2glc_with_multipleDest_gradient0") - - end subroutine test_mapLnd2glc_with_multipleDest_gradient0 - - @Test - subroutine test_mapLnd2glc_with_EC2_gradient2(this) - ! Do a test of the map_lnd2glc routine with only an elevation class 2 destination - ! point, and a gradient of 2*EC - class(TestMapLnd2glc), intent(inout) :: this - - type(lnd_data_type) :: lnd_data - - integer, parameter :: npts_glc = 1 - real(r8) :: frac_glc(npts_glc) - real(r8) :: topo_glc(npts_glc) - real(r8) :: expected_data_glc(npts_glc) - - type(simple_map_type) :: my_map - type(vertical_gradient_calculator_specified_type) :: gradient_calculator - - ! ------------------------------------------------------------------------ - ! Setup - ! ------------------------------------------------------------------------ - - call create_data_for_EC2_gradient2(my_map, gradient_calculator, lnd_data, & - frac_glc, topo_glc, expected_data_glc) - call this%setup_inputs(frac_glc, topo_glc, [lnd_data], my_map) - - ! ------------------------------------------------------------------------ - ! Exercise - ! ------------------------------------------------------------------------ - - call this%run_map_lnd2glc(gradient_calculator) - - ! ------------------------------------------------------------------------ - ! Verify - ! ------------------------------------------------------------------------ - - call this%verify_data_g(expected_data_glc, message = "test_mapLnd2glc_with_EC2_gradient2") - - end subroutine test_mapLnd2glc_with_EC2_gradient2 - - @Test - subroutine test_mapLnd2glc_with_trailing_blanks_in_fieldname(this) - ! Do a test with trailing blanks in the field name - class(TestMapLnd2glc), intent(inout) :: this - - type(lnd_data_type) :: lnd_data - - integer, parameter :: npts_glc = 1 - real(r8) :: frac_glc(npts_glc) - real(r8) :: topo_glc(npts_glc) - real(r8) :: expected_data_glc(npts_glc) - - type(simple_map_type) :: my_map - type(vertical_gradient_calculator_specified_type) :: gradient_calculator - - ! ------------------------------------------------------------------------ - ! Setup - ! ------------------------------------------------------------------------ - - call create_data_for_EC2_gradient2(my_map, gradient_calculator, lnd_data, & - frac_glc, topo_glc, expected_data_glc) - call this%setup_inputs(frac_glc, topo_glc, [lnd_data], my_map) - - ! ------------------------------------------------------------------------ - ! Exercise - ! ------------------------------------------------------------------------ - - call this%run_map_lnd2glc(gradient_calculator, fieldname = 'data ') - - ! ------------------------------------------------------------------------ - ! Verify - ! ------------------------------------------------------------------------ - - call this%verify_data_g(expected_data_glc, message = "test_mapLnd2glc_with_trailing_blanks_in_fieldname") - - end subroutine test_mapLnd2glc_with_trailing_blanks_in_fieldname - - @Test - subroutine test_mapLnd2glc_with_multipleSource_multipleDest(this) - ! Do a test of the map_lnd2glc routine with multiple source points and multiple - ! destination points. - ! - ! In particular, we have 2 source points and 3 dest points: d1 entirely in s1, d2 - ! entirely in s2, d3 part in s1, part in s2 - class(TestMapLnd2glc), intent(inout) :: this - - type(lnd_data_type) :: lnd_data(2) - - integer, parameter :: npts_glc = 3 - ! For simplicity, all glc points are in elevation class 2: - real(r8), parameter :: frac_glc(npts_glc) = [1._r8, 1._r8, 1._r8] - real(r8), parameter :: topo_glc(npts_glc) = [125._r8, 125._r8, 125._r8] - real(r8) :: expected_data_glc(npts_glc) - - type(simple_map_type) :: my_map - type(vertical_gradient_calculator_specified_type) :: gradient_calculator - - - ! ------------------------------------------------------------------------ - ! Setup - ! ------------------------------------------------------------------------ - - my_map = simple_map_type( & - source_indices = [1, 2, 1, 2], & - dest_indices = [1, 2, 3, 3], & - overlap_weights = [1._r8, 1._r8, 0.4_r8, 0.6_r8]) - - gradient_calculator = & - vgc_specified_ec_times_ptSquared(num_points = 2, & - nelev = n_elev_classes, gradient = 2._r8) - - ! data in elev class: 0 1 2 3 - lnd_data(1)%topo(:) = [25._r8, 50._r8, 150._r8, 250._r8] - lnd_data(1)%data(:) = [10._r8, 11._r8, 12._r8, 13._r8] - lnd_data(2)%topo(:) = [25._r8, 10._r8, 110._r8, 210._r8] - lnd_data(2)%data(:) = [14._r8, 15._r8, 16._r8, 17._r8] - - call this%setup_inputs(frac_glc, topo_glc, lnd_data, my_map) - - ! ------------------------------------------------------------------------ - ! Exercise - ! ------------------------------------------------------------------------ - - call this%run_map_lnd2glc(gradient_calculator) - - ! ------------------------------------------------------------------------ - ! Verify - ! ------------------------------------------------------------------------ - - ! gradient * EC * gridcell^2 = 2 * 2 * 1 = 4 - expected_data_glc(1) = lnd_data(1)%data(2) + 4._r8*(topo_glc(1) - lnd_data(1)%topo(2)) - ! gradient * EC * gridcell^2 = 2 * 2 * 4 = 16 - expected_data_glc(2) = lnd_data(2)%data(2) + 16._r8*(topo_glc(2) - lnd_data(2)%topo(2)) - - ! To determine the expected answer in the glc grid cell with two land source cells, - ! we do a straightforward application of equation 2.3 in the design document: - ! https://docs.google.com/a/ucar.edu/document/d/1sjsaiPYsPJ9A7dVGJIHGg4rVIY2qF5aRXbNzSXVAafU/edit# - ! - i.e., remap: data_l + (vertical_gradient_l) * (topo_g - topo_l) - expected_data_glc(3) = & - 0.4_r8 * (lnd_data(1)%data(2) + 4._r8*(topo_glc(3) - lnd_data(1)%topo(2))) + & - 0.6_r8 * (lnd_data(2)%data(2) + 16._r8*(topo_glc(3) - lnd_data(2)%topo(2))) - - call this%verify_data_g(expected_data_glc, message = "test_mapLnd2glc_with_multipleSource_multipleDest") - - end subroutine test_mapLnd2glc_with_multipleSource_multipleDest - - @Test - subroutine test_mapLnd2glc_with_landfrac_EC0_gradient0(this) - ! Do a test of the map_lnd2glc routine with landfrac < 1; make sure normalization - ! happens properly. This test uses elevation class 0, to cover the call to - ! seq_map_map for the bare land elevation class. In order to test the landfrac - ! normalization sufficiently, this uses two source points. - class(TestMapLnd2glc), intent(inout) :: this - - type(lnd_data_type) :: lnd_data(2) - real(r8), parameter :: lnd_overlaps(2) = [0.3_r8, 0.7_r8] - real(r8), parameter :: landfracs(2) = [0.4_r8, 0.9_r8] - - integer, parameter :: npts_glc = 1 - real(r8), parameter :: frac_glc(npts_glc) = [0._r8] - real(r8), parameter :: topo_glc(npts_glc) = [75._r8] - real(r8) :: expected_data_glc(npts_glc) - - type(simple_map_type) :: my_map - type(vertical_gradient_calculator_specified_type) :: gradient_calculator - - ! ------------------------------------------------------------------------ - ! Setup - ! ------------------------------------------------------------------------ - - my_map = simple_map_type( & - source_indices = [1 , 2], & - dest_indices = [1 , 1], & - overlap_weights = lnd_overlaps) - - gradient_calculator = & - vgc_specified_ec_times_ptSquared(num_points = 2, & - nelev = n_elev_classes, gradient = 0._r8) - - ! data in elev class: 0 1 2 3 - lnd_data(1)%topo(:) = [25._r8, 50._r8, 150._r8, 250._r8] - lnd_data(1)%data(:) = [10._r8, 11._r8, 12._r8, 13._r8] - lnd_data(2)%topo(:) = [25._r8, 10._r8, 110._r8, 210._r8] - lnd_data(2)%data(:) = [14._r8, 15._r8, 16._r8, 17._r8] - - call this%setup_inputs(frac_glc, topo_glc, lnd_data, my_map, & - landfrac = landfracs) - - ! ------------------------------------------------------------------------ - ! Exercise - ! ------------------------------------------------------------------------ - - call this%run_map_lnd2glc(gradient_calculator) - - ! ------------------------------------------------------------------------ - ! Verify - ! ------------------------------------------------------------------------ - - expected_data_glc(1) = (lnd_overlaps(1) * landfracs(1) * lnd_data(1)%data(0) + & - lnd_overlaps(2) * landfracs(2) * lnd_data(2)%data(0)) / & - (lnd_overlaps(1) * landfracs(1) + lnd_overlaps(2) * landfracs(2)) - - call this%verify_data_g(expected_data_glc, message = "test_mapLnd2glc_with_landfrac_EC0_gradient0") - - end subroutine test_mapLnd2glc_with_landfrac_EC0_gradient0 - - @Test - subroutine test_mapLnd2glc_with_landfrac_EC2_gradient0(this) - ! Do a test of the map_lnd2glc routine with landfrac < 1; make sure normalization - ! happens properly. This test uses elevation class 2, to cover the call to - ! seq_map_map for the ice sheet elevation classes. In order to test the landfrac - ! normalization sufficiently, this uses two source points. - class(TestMapLnd2glc), intent(inout) :: this - - type(lnd_data_type) :: lnd_data(2) - real(r8), parameter :: lnd_overlaps(2) = [0.3_r8, 0.7_r8] - real(r8), parameter :: landfracs(2) = [0.4_r8, 0.9_r8] - - integer, parameter :: npts_glc = 1 - real(r8), parameter :: frac_glc(npts_glc) = [1._r8] - real(r8), parameter :: topo_glc(npts_glc) = [125._r8] - real(r8) :: expected_data_glc(npts_glc) - - type(simple_map_type) :: my_map - type(vertical_gradient_calculator_specified_type) :: gradient_calculator - - ! ------------------------------------------------------------------------ - ! Setup - ! ------------------------------------------------------------------------ - - my_map = simple_map_type( & - source_indices = [1 , 2], & - dest_indices = [1 , 1], & - overlap_weights = lnd_overlaps) - - gradient_calculator = & - vgc_specified_ec_times_ptSquared(num_points = 2, & - nelev = n_elev_classes, gradient = 0._r8) - - ! data in elev class: 0 1 2 3 - lnd_data(1)%topo(:) = [25._r8, 50._r8, 150._r8, 250._r8] - lnd_data(1)%data(:) = [10._r8, 11._r8, 12._r8, 13._r8] - lnd_data(2)%topo(:) = [25._r8, 10._r8, 110._r8, 210._r8] - lnd_data(2)%data(:) = [14._r8, 15._r8, 16._r8, 17._r8] - - call this%setup_inputs(frac_glc, topo_glc, lnd_data, my_map, & - landfrac = landfracs) - - ! ------------------------------------------------------------------------ - ! Exercise - ! ------------------------------------------------------------------------ - - call this%run_map_lnd2glc(gradient_calculator) - - ! ------------------------------------------------------------------------ - ! Verify - ! ------------------------------------------------------------------------ - - expected_data_glc(1) = (lnd_overlaps(1) * landfracs(1) * lnd_data(1)%data(2) + & - lnd_overlaps(2) * landfracs(2) * lnd_data(2)%data(2)) / & - (lnd_overlaps(1) * landfracs(1) + lnd_overlaps(2) * landfracs(2)) - - call this%verify_data_g(expected_data_glc, message = "test_mapLnd2glc_with_landfrac_EC2_gradient0") - - end subroutine test_mapLnd2glc_with_landfrac_EC2_gradient0 - - @Test - subroutine test_mapLnd2glc_conservation(this) - ! This is a more complex test with multiple lnd points and multiple glc points, where - ! the different glc points are in different elevation classes. In addition to the - ! standard tests, this also demonstrates that the remapping is conservative. For this - ! to be true, though, we need to be a little more careful about how we set the - ! topographic heights of the land cells. - class(TestMapLnd2glc), intent(inout) :: this - - type(lnd_data_type) :: lnd_data(2) - - integer, parameter :: npts_glc = 4 - ! The different GLC points are in EC: 0 1 2 2 - real(r8), parameter :: frac_glc(npts_glc) = [0._r8, 1._r8, 1._r8, 1._r8] - real(r8), parameter :: topo_glc(npts_glc) = [5._r8, 25._r8, 125._r8, 190._r8] - ! Area of each glc point (needed for the conservation check): - real(r8), parameter :: area_glc(npts_glc) = [11._r8, 13._r8, 15._r8, 17._r8] - real(r8) :: expected_data_glc(npts_glc) - - type(simple_map_type) :: my_map - type(vertical_gradient_calculator_specified_type) :: gradient_calculator - - ! These parameters give the fraction of each glc cell that is in lnd1 and the - ! fraction in lnd2 (note that this sums to 1 for each glc point) - real(r8), parameter :: overlaps_with_lnd1(npts_glc) = [0.9_r8, 0.7_r8, 0.3_r8, 0.1_r8] - real(r8), parameter :: overlaps_with_lnd2(npts_glc) = [0.1_r8, 0.3_r8, 0.7_r8, 0.9_r8] - - real(r8) :: areas_lnd1(npts_glc) ! area of each glc cell in lnd1 - real(r8) :: areas_lnd2(npts_glc) ! area of each glc cell in lnd2 - real(r8) :: ec_areas_lnd1(0:n_elev_classes) ! area of each elevation class in lnd1 - real(r8) :: ec_areas_lnd2(0:n_elev_classes) ! area of each elevation class in lnd2 - real(r8) :: topo_ec2_lnd1, topo_ec2_lnd2 - real(r8) :: lnd_sum - real(r8) :: glc_sum - - ! ------------------------------------------------------------------------ - ! Setup - ! ------------------------------------------------------------------------ - - my_map = simple_map_type( & - source_indices = [1, 1, 1, 1, 2, 2, 2, 2], & - dest_indices = [1, 2, 3, 4, 1, 2, 3, 4], & - overlap_weights = [overlaps_with_lnd1, overlaps_with_lnd2]) - - gradient_calculator = & - vgc_specified_ec_times_ptSquared(num_points = 2, & - nelev = n_elev_classes, gradient = 2._r8) - - areas_lnd1(:) = overlaps_with_lnd1(:) * area_glc(:) - areas_lnd2(:) = overlaps_with_lnd2(:) * area_glc(:) - - ! Determine topographic height of elevation class 2 in each land grid cell. This - ! takes some work, because we need to take the weighted average of two CISM cells for - ! each of the two land cells. - ! - ! NOTE(wjs, 2015-05-17) In principle, we could replace this manual determination of - ! topographic heights with a call to map_glc2lnd_ec. But then it would make sense to - ! move this test routine somewhere else, to indicate that it's testing a combination - ! of routines. This would take more work than I'm up for right now. But eventually, - ! we may want to do that, along with introducing some additional conservation checks. - ! Then these conservation checks would become more integration-ish tests than unit - ! tests. - topo_ec2_lnd1 = (areas_lnd1(3) * topo_glc(3) + areas_lnd1(4) * topo_glc(4)) / & - (areas_lnd1(3) + areas_lnd1(4)) - topo_ec2_lnd2 = (areas_lnd2(3) * topo_glc(3) + areas_lnd2(4) * topo_glc(4)) / & - (areas_lnd2(3) + areas_lnd2(4)) - - ! Note that the topographic height in elevation class 1 comes from the single CISM - ! cell in elevation class 1. The topographic heights of elevation classes 0 and 3 are - ! arbitrary. (The topographic height of elevation class 0 doesn't matter for - ! conservation, and there is no CISM cell in elevation class 3.) - - ! data in elev class: 0 1 2 3 - lnd_data(1)%topo(:) = [25._r8, 25._r8, topo_ec2_lnd1, 250._r8] - lnd_data(1)%data(:) = [10._r8, 11._r8, 12._r8, 13._r8] - lnd_data(2)%topo(:) = [25._r8, 25._r8, topo_ec2_lnd2, 210._r8] - lnd_data(2)%data(:) = [14._r8, 15._r8, 16._r8, 17._r8] - - call this%setup_inputs(frac_glc, topo_glc, lnd_data, my_map) - - ! ------------------------------------------------------------------------ - ! Exercise - ! ------------------------------------------------------------------------ - - call this%run_map_lnd2glc(gradient_calculator) - - ! ------------------------------------------------------------------------ - ! Verify: Ensure that the actual glc values are the same as expected - ! ------------------------------------------------------------------------ - - ! glc 1 is in EC 0 - expected_data_glc(1) = overlaps_with_lnd1(1) * lnd_data(1)%data(0) + & - overlaps_with_lnd2(1) * lnd_data(2)%data(0) - ! glc 2 is in EC 1. The vertical gradient isn't important here, since the topographic - ! height of lnd EC 1 matches the topographic height of the glc cell. - expected_data_glc(2) = overlaps_with_lnd1(2) * lnd_data(1)%data(1) + & - overlaps_with_lnd2(2) * lnd_data(2)%data(1) - - ! glc 3 & 4 are in EC 2. Here we need to account for the vertical gradient. We do a - ! straightforward application of equation 2.3 in the design document: - ! https://docs.google.com/a/ucar.edu/document/d/1sjsaiPYsPJ9A7dVGJIHGg4rVIY2qF5aRXbNzSXVAafU/edit# - ! - i.e., remap: data_l + (vertical_gradient_l) * (topo_g - topo_l) - expected_data_glc(3) = & - overlaps_with_lnd1(3) * (lnd_data(1)%data(2) + 4._r8*(topo_glc(3) - lnd_data(1)%topo(2))) + & - overlaps_with_lnd2(3) * (lnd_data(2)%data(2) + 16._r8*(topo_glc(3) - lnd_data(2)%topo(2))) - - expected_data_glc(4) = & - overlaps_with_lnd1(4) * (lnd_data(1)%data(2) + 4._r8*(topo_glc(4) - lnd_data(1)%topo(2))) + & - overlaps_with_lnd2(4) * (lnd_data(2)%data(2) + 16._r8*(topo_glc(4) - lnd_data(2)%topo(2))) - - call this%verify_data_g(expected_data_glc, message = "test_mapLnd2glc_conservation") - - ! ------------------------------------------------------------------------ - ! Verify: Ensure conservation - ! ------------------------------------------------------------------------ - - ! Determine area of each elevation class (column) on the land grid - ec_areas_lnd1(0) = areas_lnd1(1) - ec_areas_lnd1(1) = areas_lnd1(2) - ec_areas_lnd1(2) = areas_lnd1(3) + areas_lnd1(4) - ec_areas_lnd1(3) = 0._r8 - - ec_areas_lnd2(0) = areas_lnd2(1) - ec_areas_lnd2(1) = areas_lnd2(2) - ec_areas_lnd2(2) = areas_lnd2(3) + areas_lnd2(4) - ec_areas_lnd2(3) = 0._r8 - - ! Determine weighted sum of field on the land grid - lnd_sum = sum(ec_areas_lnd1(:) * lnd_data(1)%data(:) + & - ec_areas_lnd2(:) * lnd_data(2)%data(:)) - - ! Determine weighted sum of EXPECTED field on the glc grid (we have already shown that - ! the actual field on the glc grid is the same as expected) - glc_sum = sum(area_glc(:) * expected_data_glc(:)) - - ! Show these are the same - @assertEqual(glc_sum, lnd_sum, message='Conservation', tolerance=tol) - - end subroutine test_mapLnd2glc_conservation - -end module test_map_lnd2glc diff --git a/src/drivers/mct/unit_test/stubs/CMakeLists.txt b/src/drivers/mct/unit_test/stubs/CMakeLists.txt index a3097917ed9..572c3166057 100644 --- a/src/drivers/mct/unit_test/stubs/CMakeLists.txt +++ b/src/drivers/mct/unit_test/stubs/CMakeLists.txt @@ -1,6 +1,5 @@ list(APPEND drv_sources seq_timemgr_mod.F90 - vertical_gradient_calculator_constant.F90 ) sourcelist_to_parent(drv_sources) \ No newline at end of file diff --git a/src/drivers/mct/unit_test/stubs/vertical_gradient_calculator_constant.F90 b/src/drivers/mct/unit_test/stubs/vertical_gradient_calculator_constant.F90 deleted file mode 100644 index fb4520f4491..00000000000 --- a/src/drivers/mct/unit_test/stubs/vertical_gradient_calculator_constant.F90 +++ /dev/null @@ -1,233 +0,0 @@ -module vertical_gradient_calculator_specified - - !--------------------------------------------------------------------- - ! - ! Purpose: - ! - ! This module defines a subclass of vertical_gradient_calculator_base_type that is - ! useful for unit testing. It returns a specified vertical gradient. - ! - ! This module also provides convenience functions for creating a - ! vertical_gradient_calculator_specified_type object with various functional forms. - - ! It computes the gradient as a constant times the elevation - ! class index times (the grid cell index squared). - -#include "shr_assert.h" - use vertical_gradient_calculator_base, only : vertical_gradient_calculator_base_type - use shr_kind_mod , only : r8 => shr_kind_r8 - use shr_sys_mod, only : shr_sys_abort - use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=) - - implicit none - private - - public :: vertical_gradient_calculator_specified_type - - type, extends(vertical_gradient_calculator_base_type) :: & - vertical_gradient_calculator_specified_type - private - integer :: num_points - integer :: nelev - real(r8), allocatable :: vertical_gradient(:,:) ! [point, elev classs] - logical :: calculated - contains - procedure :: calc_gradients - procedure :: get_gradients_one_class - procedure :: get_gradients_one_point - end type vertical_gradient_calculator_specified_type - - interface vertical_gradient_calculator_specified_type - module procedure constructor - end interface vertical_gradient_calculator_specified_type - - ! Creates a calculator where the gradient in ec i, pt j is gradient * i * j^2 - public :: vgc_specified_ec_times_ptSquared - - ! Creates a calculator where the gradient is constant for each point, set as the mean - ! slope from the lowest to highest elev class - public :: vgc_specified_mean_slope -contains - - !----------------------------------------------------------------------- - function vgc_specified_ec_times_ptSquared(num_points, nelev, gradient) & - result(calculator) - ! - ! !DESCRIPTION: - ! Creates a calculator where the gradient in ec i, pt j is gradient * i * j^2 - ! - ! num_points gives the number of points for which a gradient is needed (e.g., if - ! computing the vertical gradient on the land domain, then num_points is the number - ! of land points). - ! - ! !USES: - ! - ! !ARGUMENTS: - type(vertical_gradient_calculator_specified_type) :: calculator ! function result - integer, intent(in) :: num_points - integer, intent(in) :: nelev - real(r8), intent(in) :: gradient - ! - ! !LOCAL VARIABLES: - real(r8), allocatable :: gradients(:,:) - integer :: pt, ec - - character(len=*), parameter :: subname = 'vgc_specified_ec_times_ptSquared' - !----------------------------------------------------------------------- - - allocate(gradients(num_points, nelev)) - - do ec = 1, nelev - do pt = 1, num_points - gradients(pt, ec) = gradient * ec * pt**2 - end do - end do - - calculator = vertical_gradient_calculator_specified_type(gradients) - - end function vgc_specified_ec_times_ptSquared - - !----------------------------------------------------------------------- - function vgc_specified_mean_slope(data, topo) result(calculator) - ! - ! !DESCRIPTION: - ! Creates a calculator where the gradient is constant for all elevation classes - - ! though can differ for each point. Specifically, it is set to the mean slope from - ! the lowest to highest elev class - ! - ! !USES: - ! - ! !ARGUMENTS: - type(vertical_gradient_calculator_specified_type) :: calculator ! function result - real(r8), intent(in) :: data(:,:) ! [pt, ec] - real(r8), intent(in) :: topo(:,:) ! [pt, ec] - ! - ! !LOCAL VARIABLES: - integer :: num_points - integer :: nelev - real(r8), allocatable :: gradients(:,:) - integer pt - - character(len=*), parameter :: subname = 'vgc_specified_mean_slope' - !----------------------------------------------------------------------- - - num_points = size(data,1) - nelev = size(data,2) - SHR_ASSERT_ALL((ubound(topo) == (/num_points, nelev/)), 'bad size for topo') - - allocate(gradients(num_points, nelev)) - - do pt = 1, num_points - gradients(pt, :) = (data(pt,nelev) - data(pt,1)) / & - (topo(pt,nelev) - topo(pt,1)) - end do - - calculator = vertical_gradient_calculator_specified_type(gradients) - - end function vgc_specified_mean_slope - - !----------------------------------------------------------------------- - function constructor(gradients) result(this) - ! - ! !DESCRIPTION: - ! Create a new vertical_gradient_calculator_specified_type object. - ! - ! !USES: - ! - ! !ARGUMENTS: - type(vertical_gradient_calculator_specified_type) :: this ! function result - real(r8), intent(in) :: gradients(:,:) ! [pt, ec] - ! - ! !LOCAL VARIABLES: - - character(len=*), parameter :: subname = 'constructor' - !----------------------------------------------------------------------- - - this%calculated = .false. - this%num_points = size(gradients, 1) - this%nelev = size(gradients, 2) - - allocate(this%vertical_gradient(this%num_points, this%nelev)) - this%vertical_gradient(:,:) = gradients(:,:) - - end function constructor - - !----------------------------------------------------------------------- - subroutine calc_gradients(this) - ! - ! !DESCRIPTION: - ! Calculate the vertical gradients - ! - ! !USES: - ! - ! !ARGUMENTS: - class(vertical_gradient_calculator_specified_type), intent(inout) :: this - ! - ! !LOCAL VARIABLES: - - character(len=*), parameter :: subname = 'calc_gradients' - !----------------------------------------------------------------------- - - SHR_ASSERT(.not. this%calculated, 'gradients already calculated') - - ! Nothing to do in this stub - - this%calculated = .true. - - end subroutine calc_gradients - - - !----------------------------------------------------------------------- - subroutine get_gradients_one_class(this, elevation_class, gradients) - ! - ! !DESCRIPTION: - ! Return the vertical gradients for one elevation class, for all points - ! - ! !USES: - ! - ! !ARGUMENTS: - class(vertical_gradient_calculator_specified_type), intent(in) :: this - integer, intent(in) :: elevation_class - - ! gradients should already be allocated to the appropriate size - real(r8), intent(out) :: gradients(:) - ! - ! !LOCAL VARIABLES: - character(len=*), parameter :: subname = 'get_gradients_one_class' - !----------------------------------------------------------------------- - - SHR_ASSERT(this%calculated, 'gradients not yet calculated') - SHR_ASSERT(elevation_class <= this%nelev, subname//': elevation class exceeds bounds') - SHR_ASSERT((size(gradients) == this%num_points), subname//': wrong size for vertical gradient') - - gradients(:) = this%vertical_gradient(:, elevation_class) - end subroutine get_gradients_one_class - - !----------------------------------------------------------------------- - subroutine get_gradients_one_point(this, point, gradients) - ! - ! !DESCRIPTION: - ! Return the vertical gradient for all elevation classes, for one point - ! - ! !USES: - ! - ! !ARGUMENTS: - class(vertical_gradient_calculator_specified_type), intent(in) :: this - integer, intent(in) :: point - - ! gradients should already be allocated to the appropriate size - real(r8), intent(out) :: gradients(:) - ! - ! !LOCAL VARIABLES: - - character(len=*), parameter :: subname = 'get_gradients_one_class' - !----------------------------------------------------------------------- - - SHR_ASSERT(this%calculated, 'gradients not yet calculated') - SHR_ASSERT(point <= this%num_points, subname//': elevation class exceeds bounds') - SHR_ASSERT((size(gradients) == this%nelev), subname//': wrong size for vertical gradient') - - gradients(:) = this%vertical_gradient(point, :) - end subroutine get_gradients_one_point - -end module vertical_gradient_calculator_specified diff --git a/src/drivers/mct/unit_test/vertical_gradient_calculator_test/CMakeLists.txt b/src/drivers/mct/unit_test/vertical_gradient_calculator_test/CMakeLists.txt deleted file mode 100644 index ca7810516d5..00000000000 --- a/src/drivers/mct/unit_test/vertical_gradient_calculator_test/CMakeLists.txt +++ /dev/null @@ -1,9 +0,0 @@ -set (pfunit_sources - test_vertical_gradient_calculator_2nd_order.pf - test_vertical_gradient_calculator_factory.pf - ) - -create_pFUnit_test(vertical_gradient_calculator vertical_gradient_calculator_exe - "${pfunit_sources}" "") - -target_link_libraries(vertical_gradient_calculator_exe ${DRV_UNIT_TEST_LIBS}) diff --git a/src/drivers/mct/unit_test/vertical_gradient_calculator_test/README b/src/drivers/mct/unit_test/vertical_gradient_calculator_test/README deleted file mode 100644 index 7fa7d251bec..00000000000 --- a/src/drivers/mct/unit_test/vertical_gradient_calculator_test/README +++ /dev/null @@ -1,10 +0,0 @@ -The script plot_gradient and the example input file gradient_example.txt are not -directly related to the unit tests here. Rather, this script can be used to plot -the output from a gradient calculator, in a sort of "functional unit testing" -sense. - -If you look back at the history of this directory, you'll find some "functional -unit tests" that resulted in output files that could be plotted with this -script. However, these have been deleted because they were testing a vertical -gradient calculator implementation that no longer exists. - diff --git a/src/drivers/mct/unit_test/vertical_gradient_calculator_test/gradient_example.txt b/src/drivers/mct/unit_test/vertical_gradient_calculator_test/gradient_example.txt deleted file mode 100644 index 70d70dc9d0c..00000000000 --- a/src/drivers/mct/unit_test/vertical_gradient_calculator_test/gradient_example.txt +++ /dev/null @@ -1,5 +0,0 @@ -3 -0 10 20 30 -5 15 25 --3 7 15 -2 1 3 diff --git a/src/drivers/mct/unit_test/vertical_gradient_calculator_test/plot_gradient b/src/drivers/mct/unit_test/vertical_gradient_calculator_test/plot_gradient deleted file mode 100755 index bb549ecd6ae..00000000000 --- a/src/drivers/mct/unit_test/vertical_gradient_calculator_test/plot_gradient +++ /dev/null @@ -1,180 +0,0 @@ -#!/usr/bin/env python - -from __future__ import print_function - -import sys -import traceback -import os.path - -if sys.hexversion < 0x02070000: - print(70 * "*") - print("ERROR: {0} requires python >= 2.7.x. ".format(sys.argv[0])) - print("It appears that you are running python {0}".format( - ".".join(str(x) for x in sys.version_info[0:3]))) - print(70 * "*") - sys.exit(1) - -import argparse -#These are not generally available, avoid pylint error when not found -#pylint: disable=import-error -import matplotlib.pyplot as plt -import matplotlib.pylab as pylab - -class GradientInfo: - - def __init__(self, nelev, elevclass_bounds, topo, field, gradient): - """Create a GradientInfo object - - nelev: int - elevclass_bounds: tuple of (nelev+1) floats - topo: tuple of (nelev) floats - field: tuple of (nelev) floats - gradient: tuple of (nelev) floats - """ - - self.nelev = nelev - - if (len(elevclass_bounds) != nelev+1): - raise ValueError('elevclass_bounds should be of size nelev+1') - self.elevclass_bounds = elevclass_bounds - - if (len(topo) != nelev): - raise ValueError('topo should be of size nelev') - self.topo = topo - - if (len(field) != nelev): - raise ValueError('topo should be of size nelev') - self.field = field - - if (len(gradient) != nelev): - raise ValueError('gradients should be of size nelev') - self.gradient = gradient - - @classmethod - def from_file(cls, filename): - """Create a GradientInfo object by reading a file - - File should be formatted as: - nelev (int) - elevclass_bounds (list of floats; length nelev+1) - topo (list of floats; length nelev) - field (list of floats; length nelev) - gradient (list of floats; length nelev) - - For example: - 3 - 0. 10. 20. 30. - 5. 15. 25. - -3. 7. 15. - 2. 1. 3. - """ - - with open(filename) as f: - nelev = int(f.readline()) - elevclass_bounds = [float(x) for x in f.readline().split()] - topo = [float(x) for x in f.readline().split()] - field = [float(x) for x in f.readline().split()] - gradient = [float(x) for x in f.readline().split()] - - return cls(nelev, elevclass_bounds, topo, field, gradient) - - def draw_figure(self, output_filename): - """Draw a figure of this gradient info, and save it to - output_filename""" - - field_min = min(self.field) - field_max = max(self.field) - - plt.plot(self.topo, self.field, 'ro') - - # Limit upper bound of top elevation class - upper_bound = min(self.elevclass_bounds[self.nelev], - self.topo[self.nelev-1] + - (self.topo[self.nelev-1] - self.elevclass_bounds[self.nelev-1])) - - for ec in range(self.nelev): - if (ec < self.nelev - 1): - my_upper_bound = self.elevclass_bounds[ec+1] - else: - my_upper_bound = upper_bound - (xs, ys) = gradient_line(self.topo[ec], self.field[ec], self.gradient[ec], - self.elevclass_bounds[ec], my_upper_bound) - plt.plot(xs, ys, 'b') - - # limit x axes - plt.xlim([self.elevclass_bounds[0], upper_bound]) - - # set y axes ourselves, rather than letting them be dynamic, for easier - # comparison between figures - y_range = field_max - field_min - y_max = field_max + 0.2 * y_range - y_min = field_min - 0.2 * y_range - plt.ylim([y_min, y_max]) - - # plot elevation class bounds - vertical lines - # (don't draw upper bound of last EC) - for ec_bound in self.elevclass_bounds[:len(self.elevclass_bounds)-1]: - plt.plot([ec_bound, ec_bound], [y_min, y_max], 'k') - - pylab.savefig(output_filename) - plt.close() - - -def commandline_options(): - """Process command-line arguments""" - - parser = argparse.ArgumentParser( - description = 'Creates plots of gradients from one or more input files', - epilog = """Each file is expected to be formatted as follows: - nelev (int) - elevclass_bounds (list of floats; length nelev+1) - topo (list of floats; length nelev) - field (list of floats; length nelev) - gradient (list of floats; length nelev) - - For example: - 3 - 0. 10. 20. 30. - 5. 15. 25. - -3. 7. 15. - 2. 1. 3.""" - ) - - parser.add_argument('files', nargs='+', - help='names of file(s) containing gradients to plot') - - parser.add_argument('--backtrace', action='store_true', - help='show exception backtraces as extra debugging output') - - loptions = parser.parse_args() - return loptions - -def gradient_line(x, y, slope, x_lb, x_ub): - """Returns two tuples (x1, x2), (y1, y2) giving the end points of a line - that: - - - Has center (x, y) - - Has slope 'slope' - - Has x coordinates going from x_lb to x_ub - """ - - y_lb = y + (x_lb - x)*slope - y_ub = y + (x_ub - x)*slope - return ((x_lb, x_ub), (y_lb, y_ub)) - -def main(loptions): - for input_filename in loptions.files: - file_base = os.path.splitext(input_filename)[0] - gradient_info = GradientInfo.from_file(input_filename) - gradient_info.draw_figure(file_base + '.pdf') - -if __name__ == "__main__": - options = commandline_options() - try: - status = main(options) - sys.exit(status) - except Exception as error: - print(str(error)) - if options.backtrace: - traceback.print_exc() - sys.exit(1) diff --git a/src/drivers/mct/unit_test/vertical_gradient_calculator_test/test_vertical_gradient_calculator_2nd_order.pf b/src/drivers/mct/unit_test/vertical_gradient_calculator_test/test_vertical_gradient_calculator_2nd_order.pf deleted file mode 100644 index e34cbbe8bf3..00000000000 --- a/src/drivers/mct/unit_test/vertical_gradient_calculator_test/test_vertical_gradient_calculator_2nd_order.pf +++ /dev/null @@ -1,514 +0,0 @@ -module test_vertical_gradient_calculator_2nd_order - - ! Tests of vertical_gradient_calculator_2nd_order - - use pfunit_mod - use vertical_gradient_calculator_base - use vertical_gradient_calculator_2nd_order - use shr_kind_mod, only : r8 => shr_kind_r8 - implicit none - - real(r8), parameter :: tol = 1.e-13_r8 - - @TestCase - type, extends(TestCase) :: TestVertGradCalc2ndOrder - contains - procedure :: setUp - procedure :: tearDown - procedure :: create_calculator - procedure :: calculateAndVerifyGradient_1point_ECmid - end type TestVertGradCalc2ndOrder - -contains - - subroutine setUp(this) - class(TestVertGradCalc2ndOrder), intent(inout) :: this - - end subroutine setUp - - subroutine tearDown(this) - class(TestVertGradCalc2ndOrder), intent(inout) :: this - - end subroutine tearDown - - function create_calculator(this, topo, data, elevclass_bounds) & - result(calculator) - type(vertical_gradient_calculator_2nd_order_type) :: calculator - class(TestVertGradCalc2ndOrder), intent(inout) :: this - real(r8), intent(in) :: topo(:,:) ! topo(i,j) is point i, elevation class j - real(r8), intent(in) :: data(:,:) ! data(i,j) is point i, elevation class j - - ! bounds of each elevation class; this array should have one more element than the - ! number of elevation classes, since it contains lower and upper bounds for each - ! elevation class - real(r8), intent(in) :: elevclass_bounds(:) - - integer :: n_elev_classes - - n_elev_classes = size(data,2) - @assertEqual(size(data), size(topo)) - @assertEqual(n_elev_classes + 1, size(elevclass_bounds)) - - calculator = vertical_gradient_calculator_2nd_order_type( & - field = data, & - topo = topo, & - elevclass_bounds = elevclass_bounds) - call calculator%calc_gradients() - - end function create_calculator - - subroutine calculateAndVerifyGradient_1point_ECmid(this, & - elevclass_bounds, topo, data, expected_vertical_gradient, & - msg) - ! Parameterized test: Setup a vertical gradient calculator for a single point with 3 - ! ECs, calculate the vertical gradient for the middle EC, and verify that the - ! vertical gradient matches the expected vertical gradient - class(TestVertGradCalc2ndOrder), intent(inout) :: this - real(r8), intent(in) :: elevclass_bounds(:) ! elevation class bounds (should be size 4) - real(r8), intent(in) :: topo(:) ! topographic height for each EC (should be size 3) - real(r8), intent(in) :: data(:) ! data for each EC (should be size 3) - real(r8), intent(in) :: expected_vertical_gradient - character(len=*), intent(in) :: msg ! message to print if test fails - - type(vertical_gradient_calculator_2nd_order_type) :: calculator - real(r8) :: vertical_gradient(1) - - ! Check arguments - @assertEqual(4, size(elevclass_bounds)) - @assertEqual(3, size(topo)) - @assertEqual(3, size(data)) - - ! Setup - calculator = this%create_calculator( & - topo = reshape(topo, [1, 3]), & - data = reshape(data, [1, 3]), & - elevclass_bounds = elevclass_bounds) - - ! Exercise - call calculator%get_gradients_one_class(2, vertical_gradient) - - ! Verify - @assertEqual(expected_vertical_gradient, vertical_gradient(1), tolerance=tol, message = msg) - end subroutine calculateAndVerifyGradient_1point_ECmid - - @Test - subroutine ECmid(this) - ! Test with an elevation class in the middle of the range (standard case, not an edge - ! case). This uses a single grid cell. - class(TestVertGradCalc2ndOrder), intent(inout) :: this - real(r8), parameter :: elevclass_bounds(4) = [0._r8, 100._r8, 200._r8, 300._r8] - real(r8), parameter :: topo(3) = [50._r8, 125._r8, 275._r8] - real(r8), parameter :: data(3) = [11._r8, 12._r8, 13._r8] - real(r8) :: expected_vertical_gradient - - expected_vertical_gradient = (data(3) - data(1)) / (topo(3) - topo(1)) - call this%calculateAndVerifyGradient_1point_ECmid( & - elevclass_bounds = elevclass_bounds, & - topo = topo, & - data = data, & - expected_vertical_gradient = expected_vertical_gradient, & - msg = 'ECmid') - end subroutine ECmid - - @Test - subroutine ECmid_almostLimitedPositiveLB(this) - ! Make sure that a positive gradient that should *almost* (but not quite) be limited - ! by the limiter (due to the lower bound) isn't limited. - class(TestVertGradCalc2ndOrder), intent(inout) :: this - real(r8), parameter :: elevclass_bounds(4) = [0._r8, 100._r8, 200._r8, 300._r8] - real(r8), parameter :: topo(3) = [50._r8, 125._r8, 275._r8] - real(r8), parameter :: data(3) = [11._r8, 12._r8, 19.9999_r8] - real(r8) :: expected_vertical_gradient - - expected_vertical_gradient = (data(3) - data(1)) / (topo(3) - topo(1)) - call this%calculateAndVerifyGradient_1point_ECmid( & - elevclass_bounds = elevclass_bounds, & - topo = topo, & - data = data, & - expected_vertical_gradient = expected_vertical_gradient, & - msg = 'ECmid_almostLimitedPositiveLB') - end subroutine ECmid_almostLimitedPositiveLB - - @Test - subroutine ECmid_almostLimitedPositiveUB(this) - ! Make sure that a positive gradient that should *almost* (but not quite) be limited - ! by the limiter (due to the upper bound) isn't limited. - class(TestVertGradCalc2ndOrder), intent(inout) :: this - real(r8), parameter :: elevclass_bounds(4) = [0._r8, 100._r8, 200._r8, 300._r8] - real(r8), parameter :: topo(3) = [50._r8, 125._r8, 275._r8] - real(r8), parameter :: data(3) = [10.0001_r8, 12._r8, 13._r8] - real(r8) :: expected_vertical_gradient - - expected_vertical_gradient = (data(3) - data(1)) / (topo(3) - topo(1)) - call this%calculateAndVerifyGradient_1point_ECmid( & - elevclass_bounds = elevclass_bounds, & - topo = topo, & - data = data, & - expected_vertical_gradient = expected_vertical_gradient, & - msg = 'ECmid_almostLimitedPositiveUB') - end subroutine ECmid_almostLimitedPositiveUB - - @Test - subroutine ECmid_almostLimitedNegativeLB(this) - ! Make sure that a negative gradient that should *almost* (but not quite) be limited - ! by the limiter (due to the lower bound) isn't limited. - class(TestVertGradCalc2ndOrder), intent(inout) :: this - real(r8), parameter :: elevclass_bounds(4) = [0._r8, 100._r8, 200._r8, 300._r8] - real(r8), parameter :: topo(3) = [50._r8, 125._r8, 275._r8] - real(r8), parameter :: data(3) = [13._r8, 12._r8, 4.0001_r8] - real(r8) :: expected_vertical_gradient - - expected_vertical_gradient = (data(3) - data(1)) / (topo(3) - topo(1)) - call this%calculateAndVerifyGradient_1point_ECmid( & - elevclass_bounds = elevclass_bounds, & - topo = topo, & - data = data, & - expected_vertical_gradient = expected_vertical_gradient, & - msg = 'ECmid_almostLimitedNegativeLB') - end subroutine ECmid_almostLimitedNegativeLB - - @Test - subroutine ECmid_almostLimitedNegativeUB(this) - ! Make sure that a negative gradient that should *almost* (but not quite) be limited - ! by the limiter (due to the upper bound) isn't limited. - class(TestVertGradCalc2ndOrder), intent(inout) :: this - real(r8), parameter :: elevclass_bounds(4) = [0._r8, 100._r8, 200._r8, 300._r8] - real(r8), parameter :: topo(3) = [50._r8, 125._r8, 275._r8] - real(r8), parameter :: data(3) = [13.9999_r8, 12._r8, 11._r8] - real(r8) :: expected_vertical_gradient - - expected_vertical_gradient = (data(3) - data(1)) / (topo(3) - topo(1)) - call this%calculateAndVerifyGradient_1point_ECmid( & - elevclass_bounds = elevclass_bounds, & - topo = topo, & - data = data, & - expected_vertical_gradient = expected_vertical_gradient, & - msg = 'ECmid_almostLimitedNegativeUB') - end subroutine ECmid_almostLimitedNegativeUB - - @Test - subroutine ECbottom(this) - ! Test with an elevation class at the bottom edge. This uses a single grid cell. - class(TestVertGradCalc2ndOrder), intent(inout) :: this - type(vertical_gradient_calculator_2nd_order_type) :: calculator - real(r8), parameter :: elevclass_bounds(4) = [0._r8, 100._r8, 200._r8, 300._r8] - real(r8), parameter :: topo(1,3) = reshape([40._r8, 125._r8, 275._r8], [1,3]) - real(r8), parameter :: data(1,3) = reshape([11._r8, 12._r8, 13._r8], [1,3]) - real(r8) :: vertical_gradient(1) - real(r8) :: expected_vertical_gradient(1) - - calculator = this%create_calculator(topo=topo, data=data, & - elevclass_bounds=elevclass_bounds) - - call calculator%get_gradients_one_class(1, vertical_gradient) - - expected_vertical_gradient(1) = (data(1,2) - data(1,1)) / (topo(1,2) - topo(1,1)) - @assertEqual(expected_vertical_gradient, vertical_gradient, tolerance=tol) - - end subroutine ECbottom - - @Test - subroutine ECtop(this) - ! Test with an elevation class at the top edge. This uses a single grid cell. - class(TestVertGradCalc2ndOrder), intent(inout) :: this - type(vertical_gradient_calculator_2nd_order_type) :: calculator - real(r8), parameter :: elevclass_bounds(4) = [0._r8, 100._r8, 200._r8, 300._r8] - real(r8), parameter :: topo(1,3) = reshape([50._r8, 125._r8, 275._r8], [1,3]) - real(r8), parameter :: data(1,3) = reshape([11._r8, 12._r8, 13._r8], [1,3]) - real(r8) :: vertical_gradient(1) - real(r8) :: expected_vertical_gradient(1) - - calculator = this%create_calculator(topo=topo, data=data, & - elevclass_bounds=elevclass_bounds) - - call calculator%get_gradients_one_class(3, vertical_gradient) - - expected_vertical_gradient(1) = (data(1,3) - data(1,2)) / (topo(1,3) - topo(1,2)) - @assertEqual(expected_vertical_gradient, vertical_gradient, tolerance=tol) - - end subroutine ECtop - - @Test - subroutine OneEC(this) - ! Test with a single elevation class. This uses a single grid cell. - class(TestVertGradCalc2ndOrder), intent(inout) :: this - type(vertical_gradient_calculator_2nd_order_type) :: calculator - real(r8), parameter :: elevclass_bounds(2) = [0._r8, 100._r8] - real(r8), parameter :: topo(1,1) = reshape([50._r8], [1,1]) - real(r8), parameter :: data(1,1) = reshape([11._r8], [1,1]) - real(r8) :: vertical_gradient(1) - real(r8) :: expected_vertical_gradient(1) - - calculator = this%create_calculator(topo=topo, data=data, & - elevclass_bounds=elevclass_bounds) - - call calculator%get_gradients_one_class(1, vertical_gradient) - - expected_vertical_gradient(1) = 0._r8 - @assertEqual(expected_vertical_gradient, vertical_gradient, tolerance=tol) - - end subroutine OneEC - - @Test - subroutine toposEqual(this) - ! Test with topo values equal - make sure this edge case is handled correctly. - class(TestVertGradCalc2ndOrder), intent(inout) :: this - type(vertical_gradient_calculator_2nd_order_type) :: calculator - real(r8), parameter :: elevclass_bounds(3) = [0._r8, 100._r8, 200._r8] - real(r8), parameter :: topo(1,2) = reshape([100._r8, 100._r8], [1,2]) - real(r8), parameter :: data(1,2) = reshape([11._r8, 12._r8], [1,2]) - real(r8) :: vertical_gradient(1) - real(r8) :: expected_vertical_gradient(1) - - calculator = this%create_calculator(topo=topo, data=data, & - elevclass_bounds=elevclass_bounds) - - call calculator%get_gradients_one_class(2, vertical_gradient) - - expected_vertical_gradient(1) = 0._r8 - @assertEqual(expected_vertical_gradient, vertical_gradient, tolerance=tol) - - end subroutine toposEqual - - ! ------------------------------------------------------------------------ - ! Tests that trigger the limiter - ! ------------------------------------------------------------------------ - - @Test - subroutine ECmid_limitedLocalMaximum(this) - ! If values go low, high, low, then gradient should be 0 - class(TestVertGradCalc2ndOrder), intent(inout) :: this - real(r8), parameter :: elevclass_bounds(4) = [0._r8, 100._r8, 200._r8, 300._r8] - real(r8), parameter :: topo(3) = [50._r8, 125._r8, 275._r8] - real(r8), parameter :: data(3) = [11._r8, 12._r8, 10._r8] - real(r8), parameter :: expected_vertical_gradient = 0._r8 - - call this%calculateAndVerifyGradient_1point_ECmid( & - elevclass_bounds = elevclass_bounds, & - topo = topo, & - data = data, & - expected_vertical_gradient = expected_vertical_gradient, & - msg = 'ECmid_limitedLocalMaximum') - end subroutine ECmid_limitedLocalMaximum - - @Test - subroutine ECmid_limitedLocalMinimum(this) - ! If values go high, low, high, then gradient should be 0 - class(TestVertGradCalc2ndOrder), intent(inout) :: this - real(r8), parameter :: elevclass_bounds(4) = [0._r8, 100._r8, 200._r8, 300._r8] - real(r8), parameter :: topo(3) = [50._r8, 125._r8, 275._r8] - real(r8), parameter :: data(3) = [13._r8, 12._r8, 14._r8] - real(r8), parameter :: expected_vertical_gradient = 0._r8 - - call this%calculateAndVerifyGradient_1point_ECmid( & - elevclass_bounds = elevclass_bounds, & - topo = topo, & - data = data, & - expected_vertical_gradient = expected_vertical_gradient, & - msg = 'ECmid_limitedLocalMinimum') - end subroutine ECmid_limitedLocalMinimum - - @Test - subroutine ECmid_limitedPositiveLB(this) - ! Make sure that a positive gradient that should be limited by the limiter (due to the - ! lower bound) is in fact limited. - class(TestVertGradCalc2ndOrder), intent(inout) :: this - real(r8), parameter :: elevclass_bounds(4) = [0._r8, 100._r8, 200._r8, 300._r8] - real(r8), parameter :: topo(3) = [50._r8, 125._r8, 275._r8] - real(r8), parameter :: data(3) = [11._r8, 12._r8, 21._r8] - real(r8) :: expected_vertical_gradient - - expected_vertical_gradient = 1._r8/25._r8 - call this%calculateAndVerifyGradient_1point_ECmid( & - elevclass_bounds = elevclass_bounds, & - topo = topo, & - data = data, & - expected_vertical_gradient = expected_vertical_gradient, & - msg = 'ECmid_limitedPositiveLB') - end subroutine ECmid_limitedPositiveLB - - @Test - subroutine ECmid_limitedPositiveUB(this) - ! Make sure that a positive gradient that should be limited by the limiter (due to the - ! upper bound) is in fact limited. - class(TestVertGradCalc2ndOrder), intent(inout) :: this - real(r8), parameter :: elevclass_bounds(4) = [0._r8, 100._r8, 200._r8, 300._r8] - real(r8), parameter :: topo(3) = [50._r8, 125._r8, 275._r8] - real(r8), parameter :: data(3) = [9._r8, 12._r8, 13._r8] - real(r8) :: expected_vertical_gradient - - expected_vertical_gradient = 1._r8/75._r8 - call this%calculateAndVerifyGradient_1point_ECmid( & - elevclass_bounds = elevclass_bounds, & - topo = topo, & - data = data, & - expected_vertical_gradient = expected_vertical_gradient, & - msg = 'ECmid_limitedPositiveUB') - end subroutine ECmid_limitedPositiveUB - - @Test - subroutine ECmid_limitedNegativeLB(this) - ! Make sure that a negative gradient that should be limited by the limiter (due to the - ! lower bound) is in fact limited. - class(TestVertGradCalc2ndOrder), intent(inout) :: this - real(r8), parameter :: elevclass_bounds(4) = [0._r8, 100._r8, 200._r8, 300._r8] - real(r8), parameter :: topo(3) = [50._r8, 125._r8, 275._r8] - real(r8), parameter :: data(3) = [13._r8, 12._r8, 3._r8] - real(r8) :: expected_vertical_gradient - - expected_vertical_gradient = -1._r8/25._r8 - call this%calculateAndVerifyGradient_1point_ECmid( & - elevclass_bounds = elevclass_bounds, & - topo = topo, & - data = data, & - expected_vertical_gradient = expected_vertical_gradient, & - msg = 'ECmid_limitedNegativeLB') - end subroutine ECmid_limitedNegativeLB - - @Test - subroutine ECmid_limitedNegativeUB(this) - ! Make sure that a negative gradient that should be limited by the limiter (due to the - ! upper bound) is in fact limited. - class(TestVertGradCalc2ndOrder), intent(inout) :: this - real(r8), parameter :: elevclass_bounds(4) = [0._r8, 100._r8, 200._r8, 300._r8] - real(r8), parameter :: topo(3) = [50._r8, 125._r8, 275._r8] - real(r8), parameter :: data(3) = [15._r8, 12._r8, 11._r8] - real(r8) :: expected_vertical_gradient - - expected_vertical_gradient = -1._r8/75._r8 - call this%calculateAndVerifyGradient_1point_ECmid( & - elevclass_bounds = elevclass_bounds, & - topo = topo, & - data = data, & - expected_vertical_gradient = expected_vertical_gradient, & - msg = 'ECmid_limitedNegativeUB') - end subroutine ECmid_limitedNegativeUB - - ! ------------------------------------------------------------------------ - ! Test that demonstrates that we can still have non-monotonic behavior - ! - ! Unlike most tests, this test isn't necessarily something we want - it is just a - ! demonstration of current behavior. So this test can be removed if this behavior - ! changes. - ! ------------------------------------------------------------------------ - - @Test - subroutine evenWithLimiter_canStillBeNonMonotonic(this) - ! This test demonstrates that, even though the incoming values are monotonic, the - ! interpolated values are not. - ! - ! Unlike most tests, this test isn't necessarily something we want - it is just a - ! demonstration of current behavior. So this test can be removed if this behavior - ! changes. - class(TestVertGradCalc2ndOrder), intent(inout) :: this - type(vertical_gradient_calculator_2nd_order_type) :: calculator - real(r8), parameter :: elevclass_bounds(5) = [0._r8, 100._r8, 200._r8, 300._r8, 400._r8] - real(r8), parameter :: topo(1,4) = reshape([50._r8, 125._r8, 275._r8, 350._r8], [1,4]) - real(r8), parameter :: data(1,4) = reshape([9._r8, 12._r8, 13._r8 , 14._r8], [1,4]) - real(r8) :: vertical_gradient_ec2(1) - real(r8) :: vertical_gradient_ec3(1) - real(r8) :: value_200m_ec2 - real(r8) :: value_200m_ec3 - real(r8) :: value_199m - real(r8) :: value_201m - - calculator = this%create_calculator(topo=topo, data=data, & - elevclass_bounds=elevclass_bounds) - - call calculator%get_gradients_one_class(2, vertical_gradient_ec2) - call calculator%get_gradients_one_class(3, vertical_gradient_ec3) - - ! Show non-monotonicity in two ways: - - ! (1) value at 200m in EC2 > value at 200m in EC3 - value_200m_ec2 = data(1,2) + vertical_gradient_ec2(1) * (200._r8 - topo(1,2)) - value_200m_ec3 = data(1,3) + vertical_gradient_ec3(1) * (200._r8 - topo(1,3)) - @assertEqual(13._r8, value_200m_ec2, tolerance=tol) - ! In the following, use 12.9 rather than 13 to show that value_200m_ec3 is even less - ! than 12.9 (i.e., it's not just a roundoff problem) - @assertGreaterThan(12.9_r8, value_200m_ec3) - - ! (2) value at 199m (in EC2) > value at 201m (in EC3) - value_199m = data(1,2) + vertical_gradient_ec2(1) * (199._r8 - topo(1,2)) - value_201m = data(1,3) + vertical_gradient_ec3(1) * (201._r8 - topo(1,3)) - @assertGreaterThan(value_199m, value_201m) - - end subroutine evenWithLimiter_canStillBeNonMonotonic - - ! ------------------------------------------------------------------------ - ! Tests with multiple points - ! ------------------------------------------------------------------------ - - @Test - subroutine multiplePoints(this) - ! Test with multiple grid cells. One has topo values equal, two are normal cases. - class(TestVertGradCalc2ndOrder), intent(inout) :: this - type(vertical_gradient_calculator_2nd_order_type) :: calculator - - integer, parameter :: npts = 3 - integer, parameter :: nelev = 2 - real(r8), parameter :: elevclass_bounds(3) = [0._r8, 100._r8, 200._r8] - ! In the following, each line is one elevation class (with all points for that - ! elevation class) - real(r8), parameter :: topo(npts,nelev) = reshape( & - [50._r8, 100._r8, 99._r8, & - 125._r8, 100._r8, 101._r8], & - [npts,nelev]) - real(r8), parameter :: data(npts,nelev) = reshape( & - [11._r8, 100._r8, 1000._r8, & - 12._r8, 200._r8, 2000._r8], & - [npts,nelev]) - - real(r8) :: vertical_gradient(npts) - real(r8) :: expected_vertical_gradient(npts) - - calculator = this%create_calculator(topo=topo, data=data, & - elevclass_bounds=elevclass_bounds) - - call calculator%get_gradients_one_class(2, vertical_gradient) - - expected_vertical_gradient(1) = (data(1,2) - data(1,1)) / (topo(1,2) - topo(1,1)) - expected_vertical_gradient(2) = 0._r8 - expected_vertical_gradient(3) = (data(3,2) - data(3,1)) / (topo(3,2) - topo(3,1)) - @assertEqual(expected_vertical_gradient, vertical_gradient, tolerance=tol) - - end subroutine multiplePoints - - @Test - subroutine multiplePoints_someLimited(this) - ! Test with multiple grid cells, some (but not all) of which trigger the limiter. - class(TestVertGradCalc2ndOrder), intent(inout) :: this - type(vertical_gradient_calculator_2nd_order_type) :: calculator - - integer, parameter :: npts = 3 - integer, parameter :: nelev = 3 - real(r8), parameter :: elevclass_bounds(4) = [0._r8, 100._r8, 200._r8, 300._r8] - ! In the following, each line is one elevation class (with all points for that - ! elevation class) - real(r8), parameter :: topo(npts,nelev) = reshape( & - [50._r8, 50._r8, 50._r8, & - 125._r8, 125._r8, 125._r8, & - 275._r8, 275._r8, 275._r8], & - [npts,nelev]) - ! points are: limited by lower bound, non-limited, limited by upper bound - real(r8), parameter :: data(npts,nelev) = reshape( & - [11._r8, 11._r8, 9._r8, & - 12._r8, 12._r8, 12._r8, & - 21._r8, 13._r8, 13._r8], & - [npts,nelev]) - - real(r8) :: vertical_gradient(npts) - real(r8) :: expected_vertical_gradient(npts) - - calculator = this%create_calculator(topo=topo, data=data, & - elevclass_bounds=elevclass_bounds) - - call calculator%get_gradients_one_class(2, vertical_gradient) - - expected_vertical_gradient(1) = 1._r8/25._r8 - expected_vertical_gradient(2) = 2._r8/225._r8 - expected_vertical_gradient(3) = 1._r8/75._r8 - @assertEqual(expected_vertical_gradient, vertical_gradient, tolerance=tol) - end subroutine multiplePoints_someLimited - -end module test_vertical_gradient_calculator_2nd_order diff --git a/src/drivers/mct/unit_test/vertical_gradient_calculator_test/test_vertical_gradient_calculator_factory.pf b/src/drivers/mct/unit_test/vertical_gradient_calculator_test/test_vertical_gradient_calculator_factory.pf deleted file mode 100644 index d72dd917cea..00000000000 --- a/src/drivers/mct/unit_test/vertical_gradient_calculator_test/test_vertical_gradient_calculator_factory.pf +++ /dev/null @@ -1,135 +0,0 @@ -module test_vertical_gradient_calculator_factory - - ! Tests of vertical_gradient_calculator_factory - - use pfunit_mod - use vertical_gradient_calculator_factory - use shr_kind_mod , only : r8 => shr_kind_r8 - use mct_mod, only : mct_aVect, mct_aVect_clean - use mct_wrapper_mod, only : mct_init, mct_clean - use avect_wrapper_mod - - implicit none - - @TestCase - type, extends(TestCase) :: TestVertGradCalcFactory - type(mct_aVect) :: av - contains - procedure :: setUp - procedure :: tearDown - end type TestVertGradCalcFactory - - real(r8), parameter :: tol = 1.e-13_r8 - -contains - - subroutine setUp(this) - class(TestVertGradCalcFactory), intent(inout) :: this - - call mct_init() - end subroutine setUp - - subroutine tearDown(this) - class(TestVertGradCalcFactory), intent(inout) :: this - - call mct_aVect_clean(this%av) - call mct_clean() - end subroutine tearDown - - function two_digit_string(val) - ! Converts val to a two-digit string - character(len=2) :: two_digit_string - integer, intent(in) :: val - - write(two_digit_string, '(i2.2)') val - end function two_digit_string - - function elevclass_names(n_elev_classes) - ! Returns array of elevation class names - integer, intent(in) :: n_elev_classes - character(len=16) :: elevclass_names(n_elev_classes) - - integer :: i - - do i = 1, n_elev_classes - elevclass_names(i) = two_digit_string(i) - end do - end function elevclass_names - - subroutine create_av(topo, data, toponame, dataname, av) - ! Creates the attribute vector 'av' - real(r8), intent(in) :: topo(:,:) ! topo(i,j) is point i, elevation class j - real(r8), intent(in) :: data(:,:) ! data(i,j) is point i, elevation class j - character(len=*), intent(in) :: toponame - character(len=*), intent(in) :: dataname - type(mct_aVect), intent(out) :: av - - integer :: npts - integer :: n_elev_classes - integer :: elevclass - character(len=64), allocatable :: attr_tags(:) - - npts = size(topo, 1) - n_elev_classes = size(topo, 2) - - @assertEqual(ubound(data), [npts, n_elev_classes]) - - allocate(attr_tags(2*n_elev_classes)) - do elevclass = 1, n_elev_classes - attr_tags(elevclass) = dataname // two_digit_string(elevclass) - end do - do elevclass = 1, n_elev_classes - attr_tags(n_elev_classes + elevclass) = toponame // two_digit_string(elevclass) - end do - - call create_aVect_with_data_rows_are_points(av, & - attr_tags = attr_tags, & - data = reshape([data, topo], [npts, n_elev_classes * 2])) - - end subroutine create_av - - @Test - subroutine test_create_av(this) - ! Tests the create_av helper routine - class(TestVertGradCalcFactory), intent(inout) :: this - ! 3 points, 2 elevation classes - real(r8), parameter :: topo(3,2) = reshape( & - [1._r8, 2._r8, 3._r8, & - 4._r8, 5._r8, 6._r8], & - [3, 2]) - real(r8), parameter :: data(3,2) = reshape( & - [11._r8, 12._r8, 13._r8, & - 14._r8, 15._r8, 16._r8], & - [3, 2]) - - call create_av(topo, data, 'topo', 'data', this%av) - - @assertEqual([4._r8, 5._r8, 6._r8], aVect_exportRattr(this%av, 'topo' // two_digit_string(2))) - - @assertEqual([14._r8, 15._r8, 16._r8], aVect_exportRattr(this%av, 'data' // two_digit_string(2))) - - end subroutine test_create_av - - @Test - subroutine test_extract_data(this) - class(TestVertGradCalcFactory), intent(inout) :: this - integer, parameter :: npts = 2 - integer, parameter :: nelev = 3 - real(r8), parameter :: topo(npts,nelev) = & - reshape([1._r8, 2._r8, 3._r8, 4._r8, 5._r8, 6._r8], [npts, nelev]) - real(r8), parameter :: data(npts,nelev) = & - reshape([11._r8, 12._r8, 13._r8, 14._r8, 15._r8, 16._r8], [npts, nelev]) - real(r8), allocatable :: topo_extracted(:,:) - real(r8), allocatable :: data_extracted(:,:) - - call create_av(topo, data, 'topo', 'data', this%av) - - call extract_data_from_attr_vect(this%av, 'data', 'topo', elevclass_names(nelev), & - data_extracted, topo_extracted) - - @assertEqual(data, data_extracted) - @assertEqual(topo, topo_extracted) - end subroutine test_extract_data - -end module test_vertical_gradient_calculator_factory - diff --git a/src/share/test/unit/shr_string_test/test_shr_string.pf b/src/share/test/unit/shr_string_test/test_shr_string.pf index 3ac1fe603a6..82afb99aa75 100644 --- a/src/share/test/unit/shr_string_test/test_shr_string.pf +++ b/src/share/test/unit/shr_string_test/test_shr_string.pf @@ -118,6 +118,41 @@ contains @assertEqual('first:second:third:fourth', actual) end subroutine test_shr_string_listDiff_elementNotInList1 + ! ------------------------------------------------------------------------ + ! Tests of shr_string_listFromSuffixes + ! ------------------------------------------------------------------------ + + @Test + subroutine test_shr_string_listFromSuffixes_with_1() + ! 1 suffix -> list of length 1 + character(len=list_len) :: actual + + actual = shr_string_listFromSuffixes(suffixes = ['_s1'], strBase = 'foo') + @assertEqual('foo_s1', actual) + end subroutine test_shr_string_listFromSuffixes_with_1 + + @Test + subroutine test_shr_string_listFromSuffixes_with_3() + ! 3 suffixes -> list of length 3 + character(len=list_len) :: actual + + actual = shr_string_listFromSuffixes(suffixes = ['_s1', '_s2', '_s3'], strBase = 'foo') + @assertEqual('foo_s1:foo_s2:foo_s3', actual) + end subroutine test_shr_string_listFromSuffixes_with_3 + + ! ------------------------------------------------------------------------ + ! Tests of shr_string_listCreateField + ! ------------------------------------------------------------------------ + + @Test + subroutine test_shr_string_listCreateField_basic() + character(len=list_len) :: actual, expected + + actual = shr_string_listCreateField(numFields = 5, strBase = 'LAI') + expected = 'LAI_1:LAI_2:LAI_3:LAI_4:LAI_5' + @assertEqual(expected, actual) + end subroutine test_shr_string_listCreateField_basic + ! ------------------------------------------------------------------------ ! Tests of shr_string_listAddSuffix ! ------------------------------------------------------------------------ diff --git a/src/share/util/shr_string_mod.F90 b/src/share/util/shr_string_mod.F90 index 307f4390e33..03b5bc82906 100644 --- a/src/share/util/shr_string_mod.F90 +++ b/src/share/util/shr_string_mod.F90 @@ -61,6 +61,8 @@ module shr_string_mod public :: shr_string_listPrepend ! prepend list in front of another public :: shr_string_listSetDel ! Set field delimiter in lists public :: shr_string_listGetDel ! Get field delimiter in lists + public :: shr_string_listFromSuffixes! return colon delimited field list + ! given array of suffixes and a base string public :: shr_string_listCreateField ! return colon delimited field list ! given number of fields N and a base string public :: shr_string_listAddSuffix ! add a suffix to every field in a field list @@ -1729,6 +1731,43 @@ subroutine shr_string_listGetDel(del) end subroutine shr_string_listGetDel +!=============================================================================== +! +! shr_string_listFromSuffixes +! +! Returns a string of colon delimited fields given an array of suffixes and a base string +! +! given suffixes = ['_s1', '_s2', '_s3'] and strBase = 'foo', returns: +! 'foo_s1:foo_s2:foo_s3' +! +!=============================================================================== +function shr_string_listFromSuffixes( suffixes, strBase ) result ( retString ) + + character(len=*), intent(in) :: suffixes(:) + character(len=*), intent(in) :: strBase + character(len=:), allocatable :: retString + + integer :: nfields + integer :: i + integer(SHR_KIND_IN) :: t01 = 0 ! timer + + character(len=*), parameter :: subName = "(shr_string_listFromSuffixes) " + +!------------------------------------------------------------------------------- + + if ( debug > 1 .and. t01 < 1 ) call shr_timer_get( t01,subName ) + if ( debug > 1 ) call shr_timer_start( t01 ) + + nfields = size(suffixes) + retString = trim(strBase) // suffixes(1) + do i = 2, nfields + retString = trim(retString) // ':' // trim(strBase) // suffixes(i) + end do + + if ( debug > 1 ) call shr_timer_stop ( t01 ) + +end function shr_string_listFromSuffixes + !=============================================================================== ! ! shr_string_listCreateField