From 7bf66b8988df923eae518676aec38724c67addbd Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Thu, 19 Apr 2018 17:01:26 -0600 Subject: [PATCH 01/19] Pass endrun msg to shr_sys_abort I want to write some expected-exception tests for my init_interp changes. For these to work robustly, I want to check the error message returned against the expected error message (to make sure it's not aborting for some completely different reason). In doing so, I noticed that our endrun calls don't actually pass the abort message on to shr_sys_abort. This needs to be done in order to assert against the expected message. Also: add a new optional argument to endrun, additional_msg, which is written to iulog but not passed along to shr_sys_abort. This is to support expected-exception unit testing, particularly thinking about the common practice of writing out the line number: this is too volatile to test against in unit tests. So a call could look like this: call endrun(msg='Informative message', additional_msg=errmsg(__FILE__, __LINE__)) Fixes #347 --- src/main/abortutils.F90 | 38 ++++++++++++++++++++++++++------------ 1 file changed, 26 insertions(+), 12 deletions(-) diff --git a/src/main/abortutils.F90 b/src/main/abortutils.F90 index cd91e53a7d8..9f9ce562727 100644 --- a/src/main/abortutils.F90 +++ b/src/main/abortutils.F90 @@ -20,7 +20,7 @@ module abortutils CONTAINS !----------------------------------------------------------------------- - subroutine endrun_vanilla(msg) + subroutine endrun_vanilla(msg, additional_msg) !----------------------------------------------------------------------- ! !DESCRIPTION: @@ -31,21 +31,28 @@ subroutine endrun_vanilla(msg) ! ! !ARGUMENTS: implicit none - character(len=*), intent(in), optional :: msg ! string to be printed + + ! Generally you want to at least provide msg. The main reason to separate msg from + ! additional_msg is to supported expected-exception unit testing: you can put + ! volatile stuff in additional_msg, as in: + ! call endrun(msg='Informative message', additional_msg=errmsg(__FILE__, __LINE__)) + ! and then just assert against msg. + character(len=*), intent(in), optional :: msg ! string to be passed to shr_sys_abort + character(len=*), intent(in), optional :: additional_msg ! string to be printed, but not passed to shr_sys_abort !----------------------------------------------------------------------- - if (present (msg)) then - write(iulog,*)'ENDRUN:', msg + if (present (additional_msg)) then + write(iulog,*)'ENDRUN: ', additional_msg else - write(iulog,*)'ENDRUN: called without a message string' + write(iulog,*)'ENDRUN:' end if - call shr_sys_abort() + call shr_sys_abort(msg) end subroutine endrun_vanilla !----------------------------------------------------------------------- - subroutine endrun_globalindex(decomp_index, clmlevel, msg) + subroutine endrun_globalindex(decomp_index, clmlevel, msg, additional_msg) !----------------------------------------------------------------------- ! Description: @@ -59,7 +66,14 @@ subroutine endrun_globalindex(decomp_index, clmlevel, msg) implicit none integer , intent(in) :: decomp_index character(len=*) , intent(in) :: clmlevel - character(len=*) , intent(in), optional :: msg ! string to be printed + + ! Generally you want to at least provide msg. The main reason to separate msg from + ! additional_msg is to supported expected-exception unit testing: you can put + ! volatile stuff in additional_msg, as in: + ! call endrun(msg='Informative message', additional_msg=errmsg(__FILE__, __LINE__)) + ! and then just assert against msg. + character(len=*), intent(in), optional :: msg ! string to be passed to shr_sys_abort + character(len=*), intent(in), optional :: additional_msg ! string to be printed, but not passed to shr_sys_abort ! ! Local Variables: integer :: igrc, ilun, icol @@ -68,13 +82,13 @@ subroutine endrun_globalindex(decomp_index, clmlevel, msg) write(6,*)'calling getglobalwrite with decomp_index= ',decomp_index,' and clmlevel= ',trim(clmlevel) call GetGlobalWrite(decomp_index, clmlevel) - if (present (msg)) then - write(iulog,*)'ENDRUN:', msg + if (present (additional_msg)) then + write(iulog,*)'ENDRUN: ', additional_msg else - write(iulog,*)'ENDRUN: called without a message string' + write(iulog,*)'ENDRUN:' end if - call shr_sys_abort() + call shr_sys_abort(msg) end subroutine endrun_globalindex From ef23c72c93cffb8df5d60e3000967b5c7465ef99 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Fri, 20 Apr 2018 19:44:17 -0600 Subject: [PATCH 02/19] Add a new interpolation method, set_single_match This is an alternative to set_mindist that can work when there is a single matching input point for each output point (or possibly 0 matching inputs for an inactive output point). In this method, we only look for points that share the same lat/lon as the output point. This method is intended to be used for the case where we want to copy areas from the input file to the output file (rather than maintaining areas from the surface dataset); some of the logic here is set up the way it is in order to support this copying of areas. --- src/init_interp/initInterpMindist.F90 | 257 +++++++++++-- .../initInterpMindist_test/CMakeLists.txt | 7 +- .../initInterpMindistTestUtils.pf | 136 +++++++ ..._interp_mindist.pf => test_set_mindist.pf} | 163 ++------ .../test_set_single_match.pf | 359 ++++++++++++++++++ src/unit_test_shr/CMakeLists.txt | 1 + src/unit_test_shr/unittestUtils.F90 | 33 ++ 7 files changed, 774 insertions(+), 182 deletions(-) create mode 100644 src/init_interp/test/initInterpMindist_test/initInterpMindistTestUtils.pf rename src/init_interp/test/initInterpMindist_test/{test_init_interp_mindist.pf => test_set_mindist.pf} (85%) create mode 100644 src/init_interp/test/initInterpMindist_test/test_set_single_match.pf create mode 100644 src/unit_test_shr/unittestUtils.F90 diff --git a/src/init_interp/initInterpMindist.F90 b/src/init_interp/initInterpMindist.F90 index 025ce287e86..504affe3d11 100644 --- a/src/init_interp/initInterpMindist.F90 +++ b/src/init_interp/initInterpMindist.F90 @@ -25,6 +25,7 @@ module initInterpMindist ! Public methods public :: set_mindist + public :: set_single_match ! Public types @@ -53,6 +54,8 @@ module initInterpMindist ! Private methods + private :: set_glcmec_must_be_same_type + private :: set_icemec_adjustable_type private :: do_fill_missing_with_natveg private :: is_sametype private :: is_baresoil @@ -135,14 +138,10 @@ subroutine set_mindist(begi, endi, bego, endo, activei, activeo, subgridi, subgr logical :: topoglc_present logical :: closest - ! Whether each output icemec point has adjustable column type; only valid for icemec - ! points, and only valid for subgrid name = pft or column - logical :: icemec_adjustable_type_o(bego:endo) - ! Whether two glcmec columns/patches must be the same column/patch type to be ! considered the same type. This is only valid for glcmec points, and is only valid ! for subgrid name = 'pft' or 'column'. - logical :: glcmec_must_be_same_type + logical :: glcmec_must_be_same_type_o(bego:endo) ! -------------------------------------------------------------------- if (associated(subgridi%topoglc) .and. associated(subgrido%topoglc)) then @@ -154,10 +153,11 @@ subroutine set_mindist(begi, endi, bego, endo, activei, activeo, subgridi, subgr mindist_index(bego:endo) = 0 distmin = spval - call set_icemec_adjustable_type(bego=bego, endo=endo, dimname=subgrido%name, & - glc_behavior=glc_behavior, icemec_adjustable_type_o=icemec_adjustable_type_o) + call set_glcmec_must_be_same_type(bego=bego, endo=endo, dimname=subgrido%name, & + glc_elevclasses_same = glc_elevclasses_same, glc_behavior=glc_behavior, & + glcmec_must_be_same_type_o=glcmec_must_be_same_type_o) -!$OMP PARALLEL DO PRIVATE (ni,no,n,nmin,distmin,dx,dy,dist,closest,hgtdiffmin,hgtdiff,glcmec_must_be_same_type) +!$OMP PARALLEL DO PRIVATE (ni,no,n,nmin,distmin,dx,dy,dist,closest,hgtdiffmin,hgtdiff) do no = bego,endo ! Only interpolate onto active points. Otherwise, the mere act of running @@ -167,33 +167,6 @@ subroutine set_mindist(begi, endi, bego, endo, activei, activeo, subgridi, subgr ! could potentially lead to different behavior in a transient run, if those points ! later became active; that's undesirable. if (activeo(no)) then - - ! Set glcmec_must_be_same_type. Note that this only matters for glcmec types, - ! but we set it for all types for simplicity. - if (.not. glc_elevclasses_same) then - ! If the number or bounds of the elevation classes differ between input and - ! output, then we ignore the column and patch types of glcmec points when - ! looking for the same type - instead using closest topographic height as a - ! tie-breaker between equidistant columns/patches. - glcmec_must_be_same_type = .false. - else if (icemec_adjustable_type_o(no)) then - ! If glcmec points in this output cell have adjustable type, then we ignore - ! the column and patch types of glcmec points when looking for the same - ! type: we want to find the closest glcmec point in lat-lon space without - ! regards for column/patch type, because the column/patch types may change - ! at runtime. - glcmec_must_be_same_type = .false. - else - ! Otherwise, we require the column and patch types to be the same between - ! input and output for this glcmec output point, as is the case for most - ! other landunits. This is important for a case with interpolation to give - ! bit-for-bit answers with a case without interpolation (since glcmec - ! topographic heights can change after initialization, so we can't always - ! rely on the point with closest topographic height to be the "right" one to - ! pick as the source for interpolation). - glcmec_must_be_same_type = .true. - end if - nmin = 0 distmin = spval hgtdiffmin = spval @@ -202,7 +175,8 @@ subroutine set_mindist(begi, endi, bego, endo, activei, activeo, subgridi, subgr if (is_sametype(ni = ni, no = no, & subgridi = subgridi, subgrido = subgrido, & subgrid_special_indices = subgrid_special_indices, & - glcmec_must_be_same_type = glcmec_must_be_same_type)) then + glcmec_must_be_same_type = glcmec_must_be_same_type_o(no), & + veg_patch_just_considers_ptype = .true.)) then dy = abs(subgrido%lat(no)-subgridi%lat(ni))*re dx = abs(subgrido%lon(no)-subgridi%lon(ni))*re * & 0.5_r8*(subgrido%coslat(no)+subgridi%coslat(ni)) @@ -285,6 +259,195 @@ subroutine set_mindist(begi, endi, bego, endo, activei, activeo, subgridi, subgr end subroutine set_mindist + !----------------------------------------------------------------------- + subroutine set_single_match(begi, endi, bego, endo, activeo, subgridi, subgrido, & + subgrid_special_indices, glc_behavior, glc_elevclasses_same, & + mindist_index) + ! + ! !DESCRIPTION: + ! This is an alternative to set_mindist that can work when there is a single matching + ! input point for each output point (or possibly 0 matching inputs for an inactive + ! output point). In this method, we only look for points that share the same lat/lon + ! as the output point. + ! + ! This method is intended to be used for the case where we want to copy areas from the + ! input file to the output file (rather than maintaining areas from the surface + ! dataset); some of the logic here is set up the way it is in order to support this + ! copying of areas. + ! + ! !ARGUMENTS: + integer , intent(in) :: begi, endi + integer , intent(in) :: bego, endo + logical , intent(in) :: activeo(bego:endo) + type(subgrid_type) , intent(in) :: subgridi + type(subgrid_type) , intent(in) :: subgrido + type(subgrid_special_indices_type), intent(in) :: subgrid_special_indices + type(glc_behavior_type), intent(in) :: glc_behavior + + ! True if the number and bounds of glacier elevation classes are the same in the + ! input and output; false otherwise. + ! + ! Must be true for this method! (Otherwise we may not be able to find a single unique + ! input point for each output point.) + logical , intent(in) :: glc_elevclasses_same + + integer , intent(out) :: mindist_index(bego:endo) + ! + ! !LOCAL VARIABLES: + integer :: ni, no, ni_match + logical :: found + real(r8) :: dx, dy + logical :: ni_sametype + + ! Whether two glcmec columns/patches must be the same column/patch type to be + ! considered the same type. This is only valid for glcmec points, and is only valid + ! for subgrid name = 'pft' or 'column'. + logical :: glcmec_must_be_same_type_o(bego:endo) + + ! Tolerance in lat/lon for considering a point to be at the same location + real(r8) :: same_point_tol = 1.e-14_r8 + + character(len=*), parameter :: subname = 'set_single_match' + !----------------------------------------------------------------------- + + if (.not. glc_elevclasses_same) then + write(iulog,*) subname// & + ' ERROR: For this init_interp method, the number and bounds of' + write(iulog,*) 'glacier elevation classes must be the same in input and output' + call endrun(msg=subname//' ERROR: glc_elevclasses_same must be true for this method') + end if + + call set_glcmec_must_be_same_type(bego=bego, endo=endo, dimname=subgrido%name, & + glc_elevclasses_same = glc_elevclasses_same, glc_behavior=glc_behavior, & + glcmec_must_be_same_type_o=glcmec_must_be_same_type_o) + +!$OMP PARALLEL DO PRIVATE (ni,no,ni_match,found,dx,dy,ni_sametype) + do no = bego, endo + ni_match = 0 + found = .false. + do ni = begi, endi + dx = abs(subgrido%lon(no)-subgridi%lon(ni)) + dy = abs(subgrido%lat(no)-subgridi%lat(ni)) + if (dx < same_point_tol .and. dy < same_point_tol) then + ni_sametype = is_sametype(ni = ni, no = no, & + subgridi = subgridi, subgrido = subgrido, & + subgrid_special_indices = subgrid_special_indices, & + glcmec_must_be_same_type = glcmec_must_be_same_type_o(no), & + veg_patch_just_considers_ptype = .false.) + if (ni_sametype) then + if (found) then + write(iulog,*) subname// & + ' ERROR: found multiple input points matching output point' + write(iulog,*) 'For this init_interp mode: for a given output point,' + write(iulog,*) 'we expect to find exactly one input point at the same' + write(iulog,*) 'location with the same type.' + write(iulog,*) '(This most likely indicates a problem in CTSM, such as' + write(iulog,*) 'having multiple columns with the same column type on a' + write(iulog,*) 'given gridcell.)' + call endrun(msg=subname// & + ' ERROR: found multiple input points matching output point') + else + found = .true. + ni_match = ni + end if + end if + end if + end do + + if (found) then + mindist_index(no) = ni_match + else + mindist_index(no) = 0 + + if (activeo(no)) then + write(iulog,*) subname// & + ' ERROR: cannot find any input points matching output point' + write(iulog,*) 'For this init_interp mode: for a given active output point,' + write(iulog,*) 'we expect to find exactly one input point at the same' + write(iulog,*) 'location with the same type.' + call endrun(msg=subname//' ERROR: cannot find any input points matching output point') + end if + end if + end do +!$OMP END PARALLEL DO + + end subroutine set_single_match + + !----------------------------------------------------------------------- + subroutine set_glcmec_must_be_same_type(bego, endo, dimname, & + glc_elevclasses_same, glc_behavior, & + glcmec_must_be_same_type_o) + ! + ! !DESCRIPTION: + ! Sets the glcmec_must_be_same_type_o array for each output icemec point + ! + ! This array will be set to true for icemec output columns/patches for which the + ! column/patch type must match the input column/patch type to be considered the same + ! type. This is only valid for icemec points, and is only valid for dimname = 'pft' or + ! 'column' - for others, the value is undefined. + ! + ! This assumes that bego and endo match the bounds that are used elsewhere in the + ! model - i.e., for the decomposition within init_interp to match the model's + ! decomposition. + ! + ! !ARGUMENTS: + integer , intent(in) :: bego ! beginning index for output points + integer , intent(in) :: endo ! ending index for output points + character(len=*) , intent(in) :: dimname ! 'pft', 'column', etc. + + ! True if the number and bounds of glacier elevation classes are the same in the + ! input and output; false otherwise. + logical, intent(in) :: glc_elevclasses_same + + type(glc_behavior_type), intent(in) :: glc_behavior + + logical, intent(out) :: glcmec_must_be_same_type_o( bego: ) ! see description above + + ! + ! !LOCAL VARIABLES: + integer :: no + + ! Whether each output icemec point has adjustable column type; only valid for icemec + ! points, and only valid for subgrid name = pft or column + logical :: icemec_adjustable_type_o(bego:endo) + + character(len=*), parameter :: subname = 'set_glcmec_must_be_same_type' + !----------------------------------------------------------------------- + + SHR_ASSERT_ALL((ubound(glcmec_must_be_same_type_o) == (/endo/)), errMsg(sourcefile, __LINE__)) + + if (.not. glc_elevclasses_same) then + ! If the number or bounds of the elevation classes differ between input and + ! output, then we ignore the column and patch types of glcmec points when + ! looking for the same type - instead using closest topographic height as a + ! tie-breaker between equidistant columns/patches. + glcmec_must_be_same_type_o(bego:endo) = .false. + else + call set_icemec_adjustable_type(bego=bego, endo=endo, dimname=dimname, & + glc_behavior=glc_behavior, icemec_adjustable_type_o=icemec_adjustable_type_o) + do no = bego, endo + if (icemec_adjustable_type_o(no)) then + ! If glcmec points in this output cell have adjustable type, then we ignore + ! the column and patch types of glcmec points when looking for the same + ! type: we want to find the closest glcmec point in lat-lon space without + ! regards for column/patch type, because the column/patch types may change + ! at runtime. + glcmec_must_be_same_type_o(no) = .false. + else + ! Otherwise, we require the column and patch types to be the same between + ! input and output for this glcmec output point, as is the case for most + ! other landunits. This is important for a case with interpolation to give + ! bit-for-bit answers with a case without interpolation (since glcmec + ! topographic heights can change after initialization, so we can't always + ! rely on the point with closest topographic height to be the "right" one to + ! pick as the source for interpolation). + glcmec_must_be_same_type_o(no) = .true. + end if + end do + end if + + end subroutine set_glcmec_must_be_same_type + !----------------------------------------------------------------------- subroutine set_icemec_adjustable_type(bego, endo, dimname, glc_behavior, & icemec_adjustable_type_o) @@ -380,7 +543,7 @@ end function do_fill_missing_with_natveg !======================================================================= logical function is_sametype (ni, no, subgridi, subgrido, subgrid_special_indices, & - glcmec_must_be_same_type) + glcmec_must_be_same_type, veg_patch_just_considers_ptype) ! -------------------------------------------------------------------- ! arguments @@ -394,6 +557,17 @@ logical function is_sametype (ni, no, subgridi, subgrido, subgrid_special_indice ! considered the same type. This is only valid for glcmec points, and is only valid ! for subgrid name = 'pft' or 'column'. logical, intent(in) :: glcmec_must_be_same_type + + ! For vegetated patches (natural veg or crop): if veg_patch_just_considers_ptype is + ! true, then we consider two vegetated patches to be the same type if they have the + ! same ptype, without regard to their column and landunit types (as long as both the + ! input and the output are on either the natural veg or crop landunit). This is + ! needed to handle the generic crop properly when interpolating from non-crop to + ! crop, or vice versa. + ! + ! If false, then they need to have the same column and landunit types, too (as is the + ! general case). + logical, intent(in) :: veg_patch_just_considers_ptype ! -------------------------------------------------------------------- is_sametype = .false. @@ -403,11 +577,10 @@ logical function is_sametype (ni, no, subgridi, subgrido, subgrid_special_indice subgridi%ltype(ni) == subgrid_special_indices%ilun_landice_multiple_elevation_classes .and. & subgrido%ltype(no) == subgrid_special_indices%ilun_landice_multiple_elevation_classes) then is_sametype = .true. - else if (subgrid_special_indices%is_vegetated_landunit(subgrido%ltype(no))) then - ! If the output type is natural veg or crop, then just look for the correct PFT, - ! without regard for what column or landunit it's on (as long as it's on either - ! the natural veg or crop landunit). This is needed to handle the generic crop - ! properly when interpolating from non-crop to crop, or vice versa. + else if (veg_patch_just_considers_ptype .and. & + subgrid_special_indices%is_vegetated_landunit(subgrido%ltype(no))) then + ! See comment attached to the declaration of veg_patch_just_considers_ptype for + ! rationale for this logic. ! ! TODO(wjs, 2015-09-15) If we ever allow the same PFT to appear on multiple ! columns within a given grid cell, then this logic will need to be made diff --git a/src/init_interp/test/initInterpMindist_test/CMakeLists.txt b/src/init_interp/test/initInterpMindist_test/CMakeLists.txt index 3e1bf5f9c83..03d669ef343 100644 --- a/src/init_interp/test/initInterpMindist_test/CMakeLists.txt +++ b/src/init_interp/test/initInterpMindist_test/CMakeLists.txt @@ -1,4 +1,9 @@ +set (pfunit_sources + test_set_mindist.pf + test_set_single_match.pf + initInterpMindistTestUtils.pf) + create_pFUnit_test(initInterpMindist test_initInterpMindist_exe - "test_init_interp_mindist.pf" "") + "${pfunit_sources}" "") target_link_libraries(test_initInterpMindist_exe clm csm_share) \ No newline at end of file diff --git a/src/init_interp/test/initInterpMindist_test/initInterpMindistTestUtils.pf b/src/init_interp/test/initInterpMindist_test/initInterpMindistTestUtils.pf new file mode 100644 index 00000000000..009f4438bc7 --- /dev/null +++ b/src/init_interp/test/initInterpMindist_test/initInterpMindistTestUtils.pf @@ -0,0 +1,136 @@ +module initInterpMindistTestUtils + + ! Utilities to aid the testing of initInterpMindist + + use pfunit_mod + use shr_kind_mod , only : r8 => shr_kind_r8 + use initInterpMindist, only : subgrid_type, subgrid_special_indices_type + use glcBehaviorMod, only: glc_behavior_type + use unittestSubgridMod, only : bounds + use unittestArrayMod, only: grc_array + + implicit none + private + + public :: create_subgrid_info + public :: create_glc_behavior + + type(subgrid_special_indices_type), parameter, public :: subgrid_special_indices = & + subgrid_special_indices_type( & + ipft_not_vegetated = 0, & + icol_vegetated_or_bare_soil = 10, & + ilun_vegetated_or_bare_soil = 3, & + ilun_crop = 4, & + ilun_landice_multiple_elevation_classes = 5) + + ! value we can use for a special landunit; note that this just needs to differ from + ! ilun_vegetated_or_bare_soil and from ilun_crop + integer, parameter, public :: ilun_special = 6 + +contains + + !----------------------------------------------------------------------- + function create_subgrid_info(npts, name, lat, lon, & + beg, ptype, ctype, ltype, topoglc) & + result(subgrid_info) + ! + ! !ARGUMENTS: + type(subgrid_type) :: subgrid_info ! function result + integer, intent(in) :: npts + character(len=*), intent(in) :: name + real(r8), intent(in) :: lat(:) + real(r8), intent(in) :: lon(:) + integer, intent(in), optional :: beg ! beginning index; if not provided, assumed to be 1 (should be provided for output, not needed for input) + integer, intent(in), optional :: ptype(:) + integer, intent(in), optional :: ctype(:) + integer, intent(in), optional :: ltype(:) + real(r8), intent(in), optional :: topoglc(:) + ! + ! !LOCAL VARIABLES: + integer :: l_beg ! local version of beg + integer :: l_end ! ending index + + character(len=*), parameter :: subname = 'create_subgrid_info' + !----------------------------------------------------------------------- + + if (present(beg)) then + l_beg = beg + else + l_beg = 1 + end if + l_end = l_beg + npts - 1 + + ! Check array lengths + @assertEqual(npts, size(lat)) + @assertEqual(npts, size(lon)) + if (present(ptype)) then + @assertEqual(npts, size(ptype)) + end if + if (present(ctype)) then + @assertEqual(npts, size(ctype)) + end if + if (present(ltype)) then + @assertEqual(npts, size(ltype)) + end if + if (present(topoglc)) then + @assertEqual(npts, size(topoglc)) + end if + + ! Set subgrid_info components + + subgrid_info%name = name + + allocate(subgrid_info%lat(l_beg:l_end)) + subgrid_info%lat = lat + allocate(subgrid_info%lon(l_beg:l_end)) + subgrid_info%lon = lon + allocate(subgrid_info%coslat(l_beg:l_end)) + subgrid_info%coslat = cos(subgrid_info%lat) + + if (present(ptype)) then + allocate(subgrid_info%ptype(l_beg:l_end)) + subgrid_info%ptype = ptype + end if + if (present(ctype)) then + allocate(subgrid_info%ctype(l_beg:l_end)) + subgrid_info%ctype = ctype + end if + if (present(ltype)) then + allocate(subgrid_info%ltype(l_beg:l_end)) + subgrid_info%ltype = ltype + end if + if (present(topoglc)) then + allocate(subgrid_info%topoglc(l_beg:l_end)) + subgrid_info%topoglc = topoglc + end if + + end function create_subgrid_info + + !----------------------------------------------------------------------- + function create_glc_behavior(collapse_to_atm_topo) result(glc_behavior) + ! + ! !DESCRIPTION: + ! Creates a glc_behavior instance with the given collapse_to_atm_topo behavior set + ! for all grid cells. + ! + ! Must be called *after* setting up the subgrid structure. + ! + ! Note that collapse_to_atm_topo is the only aspect of glc_behavior that is relevant + ! for the unit tests in this module. + ! + ! !ARGUMENTS: + type(glc_behavior_type) :: glc_behavior ! function result + logical, intent(in) :: collapse_to_atm_topo + ! + ! !LOCAL VARIABLES: + + character(len=*), parameter :: subname = 'create_glc_behavior' + !----------------------------------------------------------------------- + + call glc_behavior%InitSetDirectly(bounds%begg, bounds%endg, & + has_virtual_columns = grc_array(.false.), & + collapse_to_atm_topo = grc_array(collapse_to_atm_topo)) + + end function create_glc_behavior + +end module initInterpMindistTestUtils diff --git a/src/init_interp/test/initInterpMindist_test/test_init_interp_mindist.pf b/src/init_interp/test/initInterpMindist_test/test_set_mindist.pf similarity index 85% rename from src/init_interp/test/initInterpMindist_test/test_init_interp_mindist.pf rename to src/init_interp/test/initInterpMindist_test/test_set_mindist.pf index 7fcccb857fb..9b291a8994a 100644 --- a/src/init_interp/test/initInterpMindist_test/test_init_interp_mindist.pf +++ b/src/init_interp/test/initInterpMindist_test/test_set_mindist.pf @@ -1,39 +1,28 @@ -module test_init_interp_mindist +module test_set_mindist - ! Tests of initInterpMindist + ! Tests of initInterpMindist: set_mindist use pfunit_mod use initInterpMindist + use initInterpMindistTestUtils, only : create_subgrid_info, create_glc_behavior + use initInterpMindistTestUtils, only : subgrid_special_indices, ilun_special use shr_kind_mod , only : r8 => shr_kind_r8 use clm_varcon , only: spval use unittestSimpleSubgridSetupsMod use unittestSubgridMod - use unittestArrayMod, only: grc_array use glcBehaviorMod, only: glc_behavior_type implicit none @TestCase - type, extends(TestCase) :: TestInitInterpMindist + type, extends(TestCase) :: TestSetMindist contains procedure :: setUp procedure :: tearDown - end type TestInitInterpMindist + end type TestSetMindist real(r8), parameter :: tol = 1.e-13_r8 - type(subgrid_special_indices_type), parameter :: subgrid_special_indices = & - subgrid_special_indices_type( & - ipft_not_vegetated = 0, & - icol_vegetated_or_bare_soil = 10, & - ilun_vegetated_or_bare_soil = 3, & - ilun_crop = 4, & - ilun_landice_multiple_elevation_classes = 5) - - ! value we can use for a special landunit; note that this just needs to differ from - ! ilun_vegetated_or_bare_soil and from ilun_crop - integer, parameter :: ilun_special = 6 - contains ! ======================================================================== @@ -42,119 +31,15 @@ contains subroutine setUp(this) - class(TestInitInterpMindist), intent(inout) :: this + class(TestSetMindist), intent(inout) :: this end subroutine setUp subroutine tearDown(this) - class(TestInitInterpMindist), intent(inout) :: this + class(TestSetMindist), intent(inout) :: this call unittest_subgrid_teardown() end subroutine tearDown - !----------------------------------------------------------------------- - function create_subgrid_info(npts, name, lat, lon, & - beg, ptype, ctype, ltype, topoglc) & - result(subgrid_info) - ! - ! !ARGUMENTS: - type(subgrid_type) :: subgrid_info ! function result - integer, intent(in) :: npts - character(len=*), intent(in) :: name - real(r8), intent(in) :: lat(:) - real(r8), intent(in) :: lon(:) - integer, intent(in), optional :: beg ! beginning index; if not provided, assumed to be 1 (should be provided for output, not needed for input) - integer, intent(in), optional :: ptype(:) - integer, intent(in), optional :: ctype(:) - integer, intent(in), optional :: ltype(:) - real(r8), intent(in), optional :: topoglc(:) - ! - ! !LOCAL VARIABLES: - integer :: l_beg ! local version of beg - integer :: l_end ! ending index - - character(len=*), parameter :: subname = 'create_subgrid_info' - !----------------------------------------------------------------------- - - if (present(beg)) then - l_beg = beg - else - l_beg = 1 - end if - l_end = l_beg + npts - 1 - - ! Check array lengths - @assertEqual(npts, size(lat)) - @assertEqual(npts, size(lon)) - if (present(ptype)) then - @assertEqual(npts, size(ptype)) - end if - if (present(ctype)) then - @assertEqual(npts, size(ctype)) - end if - if (present(ltype)) then - @assertEqual(npts, size(ltype)) - end if - if (present(topoglc)) then - @assertEqual(npts, size(topoglc)) - end if - - ! Set subgrid_info components - - subgrid_info%name = name - - allocate(subgrid_info%lat(l_beg:l_end)) - subgrid_info%lat = lat - allocate(subgrid_info%lon(l_beg:l_end)) - subgrid_info%lon = lon - allocate(subgrid_info%coslat(l_beg:l_end)) - subgrid_info%coslat = cos(subgrid_info%lat) - - if (present(ptype)) then - allocate(subgrid_info%ptype(l_beg:l_end)) - subgrid_info%ptype = ptype - end if - if (present(ctype)) then - allocate(subgrid_info%ctype(l_beg:l_end)) - subgrid_info%ctype = ctype - end if - if (present(ltype)) then - allocate(subgrid_info%ltype(l_beg:l_end)) - subgrid_info%ltype = ltype - end if - if (present(topoglc)) then - allocate(subgrid_info%topoglc(l_beg:l_end)) - subgrid_info%topoglc = topoglc - end if - - end function create_subgrid_info - - !----------------------------------------------------------------------- - function create_glc_behavior(collapse_to_atm_topo) result(glc_behavior) - ! - ! !DESCRIPTION: - ! Creates a glc_behavior instance with the given collapse_to_atm_topo behavior set - ! for all grid cells. - ! - ! Must be called *after* setting up the subgrid structure. - ! - ! Note that collapse_to_atm_topo is the only aspect of glc_behavior that is relevant - ! for the unit tests in this module. - ! - ! !ARGUMENTS: - type(glc_behavior_type) :: glc_behavior ! function result - logical, intent(in) :: collapse_to_atm_topo - ! - ! !LOCAL VARIABLES: - - character(len=*), parameter :: subname = 'create_glc_behavior' - !----------------------------------------------------------------------- - - call glc_behavior%InitSetDirectly(bounds%begg, bounds%endg, & - has_virtual_columns = grc_array(.false.), & - collapse_to_atm_topo = grc_array(collapse_to_atm_topo)) - - end function create_glc_behavior - subroutine wrap_set_mindist(subgridi, subgrido, mindist_index, activei, activeo, & glc_behavior, glc_elevclasses_same, fill_missing_with_natveg) ! Wrap the call to set_mindist. @@ -380,7 +265,7 @@ contains ! then make sure we find the point with the same type. ! ! This tests a column-level point - class(TestInitInterpMindist), intent(inout) :: this + class(TestSetMindist), intent(inout) :: this type(subgrid_type) :: subgridi, subgrido integer, parameter :: my_ctype = 15 integer, parameter :: my_ltype = 8 @@ -424,7 +309,7 @@ contains ! topographic height. ! ! This tests both column-level and patch-level - class(TestInitInterpMindist), intent(inout) :: this + class(TestSetMindist), intent(inout) :: this integer, parameter :: my_ctype = 15 integer, parameter :: my_ptype = 25 real(r8), parameter :: my_topo = 1000._r8 @@ -457,7 +342,7 @@ contains ! the same patch type, even if it isn't the one with the closest topographic height. ! ! This tests just patch-level - class(TestInitInterpMindist), intent(inout) :: this + class(TestSetMindist), intent(inout) :: this integer, parameter :: my_ctype = 15 integer, parameter :: my_ptype = 25 real(r8), parameter :: my_topo = 1000._r8 @@ -496,7 +381,7 @@ contains ! one with the closest height ! ! This tests both column-level and patch-level - class(TestInitInterpMindist), intent(inout) :: this + class(TestSetMindist), intent(inout) :: this integer, parameter :: my_ctype = 15 integer, parameter :: my_ptype = 25 real(r8), parameter :: my_topo = 1000._r8 @@ -530,7 +415,7 @@ contains ! the one with the closest height ! ! This tests both column-level and patch-level - class(TestInitInterpMindist), intent(inout) :: this + class(TestSetMindist), intent(inout) :: this integer, parameter :: my_ctype = 15 integer, parameter :: my_ptype = 25 real(r8), parameter :: my_topo = 1000._r8 @@ -562,7 +447,7 @@ contains ! For glc_mec, if we have some points closer in topographic height, but others closer ! in x-y space, pick the closer point in x-y space - class(TestInitInterpMindist), intent(inout) :: this + class(TestSetMindist), intent(inout) :: this type(subgrid_type) :: subgridi, subgrido type(glc_behavior_type) :: glc_behavior integer, parameter :: my_ctype = 15 @@ -609,7 +494,7 @@ contains ! For interpolation from a non-crop case to a crop case, ensure that a patch-level ! variable takes its input from the correct pft. This simulates what happens to the ! generic crop type in this case. - class(TestInitInterpMindist), intent(inout) :: this + class(TestSetMindist), intent(inout) :: this type(subgrid_type) :: subgridi, subgrido integer, parameter :: my_ptype = 5 integer :: my_ctype @@ -664,7 +549,7 @@ contains ! For interpolation from a crop case to a non-crop case, ensure that a patch-level ! variable takes its input from the correct pft. This simulates what happens to the ! generic crop type in this case. - class(TestInitInterpMindist), intent(inout) :: this + class(TestSetMindist), intent(inout) :: this type(subgrid_type) :: subgridi, subgrido integer, parameter :: my_ptype = 5 integer :: input_ctype @@ -718,7 +603,7 @@ contains subroutine noncropToCrop_specificCropFromNatveg(this) ! For interpolation from a non-crop case to a crop case, ensure that crop columns ! take their info from the natural veg column in the input. - class(TestInitInterpMindist), intent(inout) :: this + class(TestSetMindist), intent(inout) :: this type(subgrid_type) :: subgridi, subgrido integer :: my_ctype real(r8), parameter :: my_lat = 31._r8 @@ -772,9 +657,9 @@ contains ! bare soil point. ! ! In terms of code coverage, this is largely redundant with - ! noncrop_to_crop_specific_crop_from_natveg. But I'm keeping it for now since it - ! tests a different use case. - class(TestInitInterpMindist), intent(inout) :: this + ! noncropToCrop_specificCropFromNatveg. But I'm keeping it for now since it tests a + ! different use case. + class(TestSetMindist), intent(inout) :: this type(subgrid_type) :: subgridi, subgrido integer, parameter :: my_ptype = 5 real(r8), parameter :: my_lat = 31._r8 @@ -828,7 +713,7 @@ contains ! This test ensures that, when finding a match for a bare soil patch, we ignore ! special landunits. This is important because special landunits also have pft type = ! noveg. - class(TestInitInterpMindist), intent(inout) :: this + class(TestSetMindist), intent(inout) :: this type(subgrid_type) :: subgridi, subgrido real(r8), parameter :: my_lat = 31._r8 real(r8), parameter :: my_lon = 41._r8 @@ -881,7 +766,7 @@ contains ! from the natural veg landunit. This is especially important to check for ! patch-level variables, for which special landunits also use the noveg (0) pft type. ! So here we just check a patch-level variable. - class(TestInitInterpMindist), intent(inout) :: this + class(TestSetMindist), intent(inout) :: this type(subgrid_type) :: subgridi, subgrido integer, parameter :: my_ptype = 5 integer, parameter :: my_ctype = 15 @@ -934,7 +819,7 @@ contains @Test subroutine gridcell_findsClosestLatlon(this) ! For gridcell-level variables, should find closest gridcell in lat-lon space - class(TestInitInterpMindist), intent(inout) :: this + class(TestSetMindist), intent(inout) :: this type(subgrid_type) :: subgridi, subgrido integer :: i integer :: mindist_index(1) @@ -961,4 +846,4 @@ contains @assertEqual(2, mindist_index(1)) end subroutine gridcell_findsClosestLatlon -end module test_init_interp_mindist +end module test_set_mindist diff --git a/src/init_interp/test/initInterpMindist_test/test_set_single_match.pf b/src/init_interp/test/initInterpMindist_test/test_set_single_match.pf new file mode 100644 index 00000000000..f951a49432d --- /dev/null +++ b/src/init_interp/test/initInterpMindist_test/test_set_single_match.pf @@ -0,0 +1,359 @@ +module test_set_single_match + + ! Tests of initInterpMindist: set_single_match + + use pfunit_mod + use initInterpMindist + use initInterpMindistTestUtils, only : create_subgrid_info, create_glc_behavior + use initInterpMindistTestUtils, only : subgrid_special_indices, ilun_special + use shr_kind_mod , only : r8 => shr_kind_r8 + use unittestSimpleSubgridSetupsMod + use unittestSubgridMod + use unittestUtils, only : endrun_msg + use glcBehaviorMod, only: glc_behavior_type + + implicit none + + @TestCase + type, extends(TestCase) :: TestSetSingleMatch + contains + procedure :: setUp + procedure :: tearDown + end type TestSetSingleMatch + + real(r8), parameter :: tol = 1.e-13_r8 + +contains + + ! ======================================================================== + ! Utility routines + ! ======================================================================== + + subroutine setUp(this) + class(TestSetSingleMatch), intent(inout) :: this + end subroutine setUp + + subroutine tearDown(this) + class(TestSetSingleMatch), intent(inout) :: this + + call unittest_subgrid_teardown() + end subroutine tearDown + + subroutine wrap_set_single_match(subgridi, subgrido, mindist_index, activeo, & + glc_behavior) + ! Wrap the call to set_single_match + ! + ! If activeo is not provided, it is assumed to be .true. for all points + ! + ! If glc_behavior is not present, it is assumed to have collapse_to_atm_topo false + ! for all grid cells. + + ! Arguments: + type(subgrid_type), intent(in) :: subgridi + type(subgrid_type), intent(in) :: subgrido + integer, intent(out) :: mindist_index(:) + logical, intent(in), optional :: activeo(:) + type(glc_behavior_type), intent(in), optional :: glc_behavior + + ! Local variables: + integer :: npts_i, npts_o + integer :: bego, endo + logical, allocatable :: l_activeo(:) + type(glc_behavior_type) :: l_glc_behavior + + !----------------------------------------------------------------------- + + npts_i = size(subgridi%lon) + npts_o = size(subgrido%lon) + bego = lbound(subgrido%lon, 1) + endo = ubound(subgrido%lon, 1) + + @assertEqual(npts_o, size(mindist_index)) + + if (present(activeo)) then + @assertEqual(npts_o, size(activeo)) + l_activeo = activeo + else + allocate(l_activeo(npts_o)) + l_activeo = .true. + end if + + if (present(glc_behavior)) then + l_glc_behavior = glc_behavior + else + l_glc_behavior = create_glc_behavior(collapse_to_atm_topo = .false.) + end if + + call set_single_match(begi = 1, endi = npts_i, bego = bego, endo = endo, & + activeo = l_activeo, subgridi = subgridi, subgrido = subgrido, & + subgrid_special_indices = subgrid_special_indices, & + glc_behavior = l_glc_behavior, & + glc_elevclasses_same = .true., & + mindist_index = mindist_index) + + end subroutine wrap_set_single_match + + ! ======================================================================== + ! Tests + ! ======================================================================== + + @Test + subroutine singleMatch_findsMatch(this) + class(TestSetSingleMatch), intent(inout) :: this + type(subgrid_type) :: subgridi, subgrido + integer, parameter :: my_ctype = 15 + integer, parameter :: my_ltype = 8 + real(r8), parameter :: my_lat = 31._r8 + real(r8), parameter :: my_lon = 41._r8 + integer :: mindist_index(1) + integer :: i + + call setup_single_veg_patch(pft_type=1) + + subgrido = create_subgrid_info( & + npts = 1, & + beg = bounds%begc, & + name = 'column', & + ctype = [my_ctype], & + ltype = [my_ltype], & + lat = [my_lat], & + lon = [my_lon]) + + ! The target point is point 3. Both before and after the target point there are + ! points with (1) same type but different location, and (2) same location but + ! different type. + subgridi = create_subgrid_info( & + npts = 5, & + name = 'column', & + ctype = [my_ctype-1, my_ctype, my_ctype, my_ctype, my_ctype+1], & + ltype = [(my_ltype, i=1,5)], & + lat = [my_lat, my_lat+1, my_lat, my_lat-1, my_lat], & + lon = [(my_lon, i=1,5)]) + + call wrap_set_single_match(subgridi, subgrido, mindist_index) + + @assertEqual(3, mindist_index(1)) + end subroutine singleMatch_findsMatch + + @Test + subroutine singleMatch_inactive_findsMatch(this) + ! Even if the output point is inactive, it still finds a match for this point. (This + ! is in contrast to set_mindist). + class(TestSetSingleMatch), intent(inout) :: this + type(subgrid_type) :: subgridi, subgrido + integer, parameter :: my_ctype = 15 + integer, parameter :: my_ltype = 8 + real(r8), parameter :: my_lat = 31._r8 + real(r8), parameter :: my_lon = 41._r8 + integer :: mindist_index(1) + integer :: i + + call setup_single_veg_patch(pft_type=1) + + subgrido = create_subgrid_info( & + npts = 1, & + beg = bounds%begc, & + name = 'column', & + ctype = [my_ctype], & + ltype = [my_ltype], & + lat = [my_lat], & + lon = [my_lon]) + + subgridi = create_subgrid_info( & + npts = 3, & + name = 'column', & + ctype = [my_ctype-1, my_ctype, my_ctype+1], & + ltype = [(my_ltype, i=1,3)], & + lat = [(my_lat, i=1,3)], & + lon = [(my_lon, i=1,3)]) + + call wrap_set_single_match(subgridi, subgrido, mindist_index, & + activeo = [.false.]) + + @assertEqual(2, mindist_index(1)) + end subroutine singleMatch_inactive_findsMatch + + @Test + subroutine noMatches_inactive_returns0(this) + ! For an inactive point: it's okay if there are no candidate source points + class(TestSetSingleMatch), intent(inout) :: this + type(subgrid_type) :: subgridi, subgrido + integer, parameter :: my_ctype = 15 + integer, parameter :: my_ltype = 8 + real(r8), parameter :: my_lat = 31._r8 + real(r8), parameter :: my_lon = 41._r8 + integer :: mindist_index(1) + + call setup_single_veg_patch(pft_type=1) + + subgrido = create_subgrid_info( & + npts = 1, & + beg = bounds%begc, & + name = 'column', & + ctype = [my_ctype], & + ltype = [my_ltype], & + lat = [my_lat], & + lon = [my_lon]) + + subgridi = create_subgrid_info( & + npts = 1, & + name = 'column', & + ctype = [my_ctype+1], & + ltype = [my_ltype], & + lat = [my_lat], & + lon = [my_lon]) + + call wrap_set_single_match(subgridi, subgrido, mindist_index, & + activeo = [.false.]) + + @assertEqual(0, mindist_index(1)) + end subroutine noMatches_inactive_returns0 + + @Test + subroutine twoMatches_aborts(this) + ! If there are two matches for the given output point, aborts + class(TestSetSingleMatch), intent(inout) :: this + type(subgrid_type) :: subgridi, subgrido + integer, parameter :: my_ctype = 15 + integer, parameter :: my_ltype = 8 + real(r8), parameter :: my_lat = 31._r8 + real(r8), parameter :: my_lon = 41._r8 + integer :: mindist_index(1) + integer :: i + character(len=:), allocatable :: expected_msg + + call setup_single_veg_patch(pft_type=1) + + subgrido = create_subgrid_info( & + npts = 1, & + beg = bounds%begc, & + name = 'column', & + ctype = [my_ctype], & + ltype = [my_ltype], & + lat = [my_lat], & + lon = [my_lon]) + + subgridi = create_subgrid_info( & + npts = 2, & + name = 'column', & + ctype = [(my_ctype, i=1,2)], & + ltype = [(my_ltype, i=1,2)], & + lat = [(my_lat, i=1,2)], & + lon = [(my_lon, i=1,2)]) + + call wrap_set_single_match(subgridi, subgrido, mindist_index) + + expected_msg = endrun_msg( & + 'set_single_match ERROR: found multiple input points matching output point') + @assertExceptionRaised(expected_msg) + end subroutine twoMatches_aborts + + @Test + subroutine noMatches_aborts(this) + ! For an active point: aborts if there are no candidate source points + ! + ! Note for the future: We could probably relax this requirement if we added code that + ! set subgrid areas to zero for any active point in the output for which there are no + ! matching input points. + class(TestSetSingleMatch), intent(inout) :: this + type(subgrid_type) :: subgridi, subgrido + integer, parameter :: my_ctype = 15 + integer, parameter :: my_ltype = 8 + real(r8), parameter :: my_lat = 31._r8 + real(r8), parameter :: my_lon = 41._r8 + integer :: mindist_index(1) + character(len=:), allocatable :: expected_msg + + call setup_single_veg_patch(pft_type=1) + + subgrido = create_subgrid_info( & + npts = 1, & + beg = bounds%begc, & + name = 'column', & + ctype = [my_ctype], & + ltype = [my_ltype], & + lat = [my_lat], & + lon = [my_lon]) + + ! One point differs in lat, one point differs in lon, one point differs in ctype + subgridi = create_subgrid_info( & + npts = 3, & + name = 'column', & + ctype = [my_ctype, my_ctype, my_ctype+1], & + ltype = [my_ltype, my_ltype, my_ltype], & + lat = [my_lat+1, my_lat , my_lat], & + lon = [my_lon , my_lon+1, my_lon]) + + call wrap_set_single_match(subgridi, subgrido, mindist_index) + + expected_msg = endrun_msg( & + 'set_single_match ERROR: cannot find any input points matching output point') + @assertExceptionRaised(expected_msg) + + end subroutine noMatches_aborts + + @Test + subroutine noncropToCrop_patchVariable_aborts(this) + ! This test provides a contrast with the analogous test of set_mindist, + ! noncropToCrop_patchVariable_usesCorrectPft. For set_single_patch, in contrast to + ! set_mindist, it is an error if we try to go from a non-crop case to a crop case. + ! This is because areas won't add to 1 properly if we try to copy a crop patch from + ! the natural veg column onto a crop patch from its own crop column. + ! + ! This test also covers similar use cases, such as going from a crop case to a + ! non-crop case. + class(TestSetSingleMatch), intent(inout) :: this + type(subgrid_type) :: subgridi, subgrido + integer, parameter :: my_ptype = 5 + integer :: my_ctype + real(r8), parameter :: my_lat = 31._r8 + real(r8), parameter :: my_lon = 41._r8 + integer :: mindist_index(1) + integer :: i + character(len=:), allocatable :: expected_msg + + associate( & + icol_natveg => subgrid_special_indices%icol_vegetated_or_bare_soil, & + ilun_natveg => subgrid_special_indices%ilun_vegetated_or_bare_soil, & + ilun_crop => subgrid_special_indices%ilun_crop & + ) + + my_ctype = icol_natveg + 1 ! arbitrary; we just want this to differ from icol_natveg + + call setup_landunit_ncols(ltype=ilun_crop, & + ctypes=[my_ctype], & + cweights=[1._r8], & + ptype=my_ptype) + + subgrido = create_subgrid_info( & + npts = 1, & + beg = bounds%begp, & + name = 'pft', & + ptype = [my_ptype], & + ctype = [my_ctype], & + ltype = [ilun_crop], & + lat = [my_lat], & + lon = [my_lon]) + + ! Input point #2 has the same ptype, but a different ctype and ltype. Other input + ! points differ in ptype. With set_mindist, we would choose #2, but with + ! set_single_match we should abort. + subgridi = create_subgrid_info( & + npts = 3, & + name = 'pft', & + ptype = [my_ptype - 1, my_ptype, my_ptype + 1], & + ctype = [icol_natveg, icol_natveg, icol_natveg], & + ltype = [ilun_natveg, ilun_natveg, ilun_natveg], & + lat = [(my_lat, i=1,3)], & + lon = [(my_lon, i=1,3)]) + + call wrap_set_single_match(subgridi, subgrido, mindist_index) + + expected_msg = endrun_msg( & + 'set_single_match ERROR: cannot find any input points matching output point') + @assertExceptionRaised(expected_msg) + + end associate + end subroutine noncropToCrop_patchVariable_aborts + +end module test_set_single_match diff --git a/src/unit_test_shr/CMakeLists.txt b/src/unit_test_shr/CMakeLists.txt index a6ba7a87889..8e3494f8c75 100644 --- a/src/unit_test_shr/CMakeLists.txt +++ b/src/unit_test_shr/CMakeLists.txt @@ -14,6 +14,7 @@ list(APPEND clm_sources unittestSimpleSubgridSetupsMod.F90 unittestFilterBuilderMod.F90 unittestGlcMec.F90 + unittestUtils.F90 ) sourcelist_to_parent(clm_sources) diff --git a/src/unit_test_shr/unittestUtils.F90 b/src/unit_test_shr/unittestUtils.F90 new file mode 100644 index 00000000000..9121a2e9eb4 --- /dev/null +++ b/src/unit_test_shr/unittestUtils.F90 @@ -0,0 +1,33 @@ +module unittestUtils + + ! Miscellaneous utilities to aid unit testing + + implicit none + private + + public :: endrun_msg ! Gives the message thrown by shr_abort_abort, given a call to endrun(msg) + +contains + + !----------------------------------------------------------------------- + function endrun_msg(msg) + ! + ! !DESCRIPTION: + ! Gives the message thrown by shr_abort_abort, given a call to endrun(msg) + ! + ! !USES: + ! + ! !ARGUMENTS: + character(len=:), allocatable :: endrun_msg ! function result + character(len=*), intent(in) :: msg + ! + ! !LOCAL VARIABLES: + + character(len=*), parameter :: subname = 'endrun_msg' + !----------------------------------------------------------------------- + + endrun_msg = 'ABORTED: '//trim(msg) + + end function endrun_msg + +end module unittestUtils From da5171d5baf6ded1170343b1f0bc8bfb2d892254 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Tue, 24 Apr 2018 09:00:04 -0600 Subject: [PATCH 03/19] Add logic controlling init_interp method Add a namelist flag controlling which init_interp method to use. This controls two aspects of the behavior: (1) which initInterpMindist subroutine to call (2) whether to interpolate area-related fields --- bld/CLMBuildNamelist.pm | 17 +++ .../namelist_defaults_clm4_5.xml | 6 + .../namelist_definition_clm4_5.xml | 8 ++ src/init_interp/initInterp.F90 | 131 ++++++++++++++++-- src/main/TopoMod.F90 | 5 +- src/main/subgridRestMod.F90 | 12 +- src/utils/restUtilMod.F90.in | 89 ++++++++---- 7 files changed, 217 insertions(+), 51 deletions(-) diff --git a/bld/CLMBuildNamelist.pm b/bld/CLMBuildNamelist.pm index 00807d049c4..568e663b6c1 100755 --- a/bld/CLMBuildNamelist.pm +++ b/bld/CLMBuildNamelist.pm @@ -1745,6 +1745,10 @@ sub process_namelist_inline_logic { ####################################################################### setup_logic_hydrology_switches($nl, $physv); + ######################################### + # namelist group: clm_initinterp_inparm # + ######################################### + setup_logic_initinterp($opts, $nl_flags, $definition, $defaults, $nl, $physv); } #------------------------------------------------------------------------------- @@ -3723,6 +3727,19 @@ sub setup_logic_lnd2atm { #------------------------------------------------------------------------------- +sub setup_logic_initinterp { + # + # Options related to init_interp + # + my ($opts, $nl_flags, $definition, $defaults, $nl, $physv) = @_; + + if ($physv->as_long() >= $physv->as_long("clm4_5")) { + add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'init_interp_method'); + } +} + +#------------------------------------------------------------------------------- + sub setup_logic_fates { # # Set some default options related to Ecosystem Demography diff --git a/bld/namelist_files/namelist_defaults_clm4_5.xml b/bld/namelist_files/namelist_defaults_clm4_5.xml index bf45c8f8fae..6a165a9f8fc 100644 --- a/bld/namelist_files/namelist_defaults_clm4_5.xml +++ b/bld/namelist_files/namelist_defaults_clm4_5.xml @@ -2413,4 +2413,10 @@ lnd/clm2/surfdata_map/surfdata_ne120np4_78pfts_CMIP6_simyr1850_c170824.nc.false. .false. + + + + +general + diff --git a/bld/namelist_files/namelist_definition_clm4_5.xml b/bld/namelist_files/namelist_definition_clm4_5.xml index 44b98a18b43..92614eced90 100644 --- a/bld/namelist_files/namelist_definition_clm4_5.xml +++ b/bld/namelist_files/namelist_definition_clm4_5.xml @@ -2416,4 +2416,12 @@ always fill with the closest natural veg patch / column, regardless of the value flag. So interpolation from non-crop to crop cases can be done without setting this flag. + +Method to use for init_interp + +Only applies when use_init_interp = .true. + + diff --git a/src/init_interp/initInterp.F90 b/src/init_interp/initInterp.F90 index d979c6529de..a5f31584bc6 100644 --- a/src/init_interp/initInterp.F90 +++ b/src/init_interp/initInterp.F90 @@ -7,7 +7,8 @@ module initInterpMod #include "shr_assert.h" use initInterpBounds, only : interp_bounds_type - use initInterpMindist, only: set_mindist, subgrid_type, subgrid_special_indices_type + use initInterpMindist, only: set_mindist, set_single_match + use initInterpMindist, only: subgrid_type, subgrid_special_indices_type use initInterp1dData, only : interp_1d_data use initInterp2dvar, only: interp_2dvar_type use initInterpMultilevelBase, only : interp_multilevel_type @@ -21,7 +22,7 @@ module initInterpMod use clm_varctl , only: iulog use abortutils , only: endrun use spmdMod , only: masterproc - use restUtilMod , only: iflag_interp, iflag_copy, iflag_skip + use restUtilMod , only: iflag_interp, iflag_copy, iflag_skip, iflag_area use glcBehaviorMod , only: glc_behavior_type use ncdio_utils , only: find_var_on_file use ncdio_pio @@ -40,6 +41,7 @@ module initInterpMod private :: check_dim_subgrid private :: check_dim_level + private :: skip_var private :: findMinDist private :: set_subgrid_info private :: interp_0d_copy @@ -52,6 +54,20 @@ module initInterpMod character(len=8) :: created_glacier_mec_landunits + ! Allowable interpolation methods + ! - interp_method_general: The general-purpose method + ! - interp_method_finidat_areas: Take areas and other related fields from the finidat + ! file, rather than maintaining them at the values from the surface dataset. This + ! only works for interpolation from one resolution to the same resolution, and with + ! basically the same configuration. This also triggers different logic for finding + ! matching points, assuring that we find exactly one match in the input for each + ! output point. + integer, parameter :: interp_method_general = 1 + integer, parameter :: interp_method_finidat_areas = 2 + + ! One of the above methods + integer :: interp_method + ! If true, fill missing types with closest natural veg column (using bare soil for ! patch-level variables) logical :: init_interp_fill_missing_with_natveg @@ -81,14 +97,16 @@ subroutine initInterp_readnl(NLFilename) ! !LOCAL VARIABLES: integer :: ierr ! error code integer :: unitn ! unit for namelist file + character(len=64) :: init_interp_method character(len=*), parameter :: subname = 'initInterp_readnl' !----------------------------------------------------------------------- namelist /clm_initinterp_inparm/ & - init_interp_fill_missing_with_natveg + init_interp_method, init_interp_fill_missing_with_natveg ! Initialize options to default values, in case they are not specified in the namelist + init_interp_method = ' ' init_interp_fill_missing_with_natveg = .false. if (masterproc) then @@ -107,6 +125,7 @@ subroutine initInterp_readnl(NLFilename) call relavu( unitn ) end if + call shr_mpi_bcast (init_interp_method, mpicom) call shr_mpi_bcast (init_interp_fill_missing_with_natveg, mpicom) if (masterproc) then @@ -116,6 +135,32 @@ subroutine initInterp_readnl(NLFilename) write(iulog,*) ' ' end if + call translate_options(init_interp_method) + call check_nl_consistency() + + contains + subroutine translate_options(init_interp_method) + ! Translate string options into integer parameters + character(len=*), intent(in) :: init_interp_method + + select case (init_interp_method) + case ('general') + interp_method = interp_method_general + case ('use_finidat_areas') + interp_method = interp_method_finidat_areas + case default + write(iulog,*) 'Unknown value for init_interp_method: ', trim(init_interp_method) + call endrun('Unknown value for init_interp_method') + end select + end subroutine translate_options + + subroutine check_nl_consistency() + if (interp_method == interp_method_finidat_areas .and. & + init_interp_fill_missing_with_natveg) then + call endrun(msg='init_interp_method = use_finidat_areas is incompatible with init_interp_fill_missing_with_natveg') + end if + end subroutine check_nl_consistency + end subroutine initInterp_readnl @@ -396,7 +441,7 @@ subroutine initInterp (filei, fileo, bounds, glc_behavior) status = pio_inq_varid (ncido, trim(varname), vardesc) status = pio_get_att(ncido, vardesc, 'interpinic_flag', iflag_interpinic) - if (iflag_interpinic == iflag_skip) then + if (skip_var(iflag_interpinic)) then if (masterproc) then write (iulog,*) 'Skipping : ', trim(varname) end if @@ -507,6 +552,14 @@ subroutine initInterp (filei, fileo, bounds, glc_behavior) write(iulog,*) 'Skipping : ',trim(varname) end if + else + + if (masterproc) then + write(iulog,*) 'Bad interpinic flag for scalar variable ', trim(varname), & + ': ', iflag_interpinic + call endrun(msg='Bad interpinic flag for scalar variable'//errMsg(sourcefile, __LINE__)) + end if + end if !--------------------------------------------------- @@ -640,6 +693,40 @@ subroutine initInterp (filei, fileo, bounds, glc_behavior) end subroutine initInterp + !----------------------------------------------------------------------- + function skip_var(iflag_interpinic) + ! + ! !DESCRIPTION: + ! Logical function saying whether we should skip a given variable given + ! iflag_interpinic and interp_method + ! + ! !ARGUMENTS: + logical :: skip_var ! function result + integer, intent(in) :: iflag_interpinic + ! + ! !LOCAL VARIABLES: + + character(len=*), parameter :: subname = 'skip_var' + !----------------------------------------------------------------------- + + if (iflag_interpinic == iflag_skip) then + skip_var = .true. + else if (iflag_interpinic == iflag_area) then + select case (interp_method) + case (interp_method_general) + skip_var = .true. + case (interp_method_finidat_areas) + skip_var = .false. + case default + call endrun(msg='Unhandled interp_method'//errMsg(sourcefile, __LINE__)) + end select + else + skip_var = .false. + end if + + end function skip_var + + !======================================================================= subroutine findMinDist( dimname, begi, endi, bego, endo, ncidi, ncido, & @@ -680,16 +767,32 @@ subroutine findMinDist( dimname, begi, endi, bego, endo, ncidi, ncido, & call set_subgrid_info(beg=bego, end=endo, dimname=dimname, use_glob=.false., & ncid=ncido, active=activeo, subgrid=subgrido) - if (masterproc) then - write(iulog,*)'calling set_mindist for ',trim(dimname) - end if - call set_mindist(begi=begi, endi=endi, bego=bego, endo=endo, & - activei=activei, activeo=activeo, subgridi=subgridi, subgrido=subgrido, & - subgrid_special_indices=subgrid_special_indices, & - glc_behavior=glc_behavior, & - glc_elevclasses_same=glc_elevclasses_same, & - fill_missing_with_natveg=init_interp_fill_missing_with_natveg, & - mindist_index=minindx) + select case (interp_method) + case (interp_method_general) + if (masterproc) then + write(iulog,*) 'calling set_mindist for ',trim(dimname) + end if + call set_mindist(begi=begi, endi=endi, bego=bego, endo=endo, & + activei=activei, activeo=activeo, subgridi=subgridi, subgrido=subgrido, & + subgrid_special_indices=subgrid_special_indices, & + glc_behavior=glc_behavior, & + glc_elevclasses_same=glc_elevclasses_same, & + fill_missing_with_natveg=init_interp_fill_missing_with_natveg, & + mindist_index=minindx) + case (interp_method_finidat_areas) + if (masterproc) then + write(iulog,*) 'calling set_single_match for ',trim(dimname) + end if + call set_single_match(begi=begi, endi=endi, bego=bego, endo=endo, & + activeo=activeo, subgridi=subgridi, subgrido=subgrido, & + subgrid_special_indices=subgrid_special_indices, & + glc_behavior=glc_behavior, & + glc_elevclasses_same=glc_elevclasses_same, & + mindist_index=minindx) + case default + write(iulog,*) 'findMinDist: unhandled interp_method: ', interp_method + call endrun('Unhandled interp_method'//errMsg(sourcefile, __LINE__)) + end select deallocate(subgridi%lat, subgridi%lon, subgridi%coslat) deallocate(subgrido%lat, subgrido%lon, subgrido%coslat) diff --git a/src/main/TopoMod.F90 b/src/main/TopoMod.F90 index c3fa7ba1e90..43b935bed00 100644 --- a/src/main/TopoMod.F90 +++ b/src/main/TopoMod.F90 @@ -171,16 +171,19 @@ subroutine Restart(this, bounds, ncid, flag) ! TODO(wjs, 2016-04-05) Rename these restart variables to get rid of 'glc' in their ! names. However, this will require some changes to init_interp, too. + ! This one is not actually an area, but has interpinic_flag='area' because we want to + ! interpolate it under the same conditions under which we interpolate areas. call restartvar(ncid=ncid, flag=flag, varname='cols1d_topoglc', xtype=ncd_double, & dim1name='column', & long_name='mean elevation on glacier elevation classes', units='m', & - interpinic_flag='skip', readvar=readvar, data=this%topo_col) + interpinic_flag='area', readvar=readvar, data=this%topo_col) if (flag /= 'read') then do p=bounds%begp,bounds%endp c = patch%column(p) rparr(p) = this%topo_col(c) enddo + ! This one has interpinic_flag = 'skip' because it isn't read back in call restartvar(ncid=ncid, flag=flag, varname='pfts1d_topoglc', xtype=ncd_double, & dim1name='pft', & long_name='mean elevation on glacier elevation classes', units='m', & diff --git a/src/main/subgridRestMod.F90 b/src/main/subgridRestMod.F90 index 757a960cfe1..62370613937 100644 --- a/src/main/subgridRestMod.F90 +++ b/src/main/subgridRestMod.F90 @@ -478,32 +478,32 @@ subroutine subgridRest_write_and_read(bounds, ncid, flag) call restartvar(ncid=ncid, flag=flag, varname='land1d_wtxy', xtype=ncd_double, & dim1name='landunit', & long_name='landunit weight relative to corresponding gridcell', & - interpinic_flag='skip', readvar=readvar, data=lun%wtgcell) + interpinic_flag='area', readvar=readvar, data=lun%wtgcell) call restartvar(ncid=ncid, flag=flag, varname='cols1d_wtxy', xtype=ncd_double, & dim1name='column', & long_name='column weight relative to corresponding gridcell', units=' ', & - interpinic_flag='skip', readvar=readvar, data=col%wtgcell) + interpinic_flag='area', readvar=readvar, data=col%wtgcell) call restartvar(ncid=ncid, flag=flag, varname='cols1d_wtlnd', xtype=ncd_double, & dim1name='column', & long_name='column weight relative to corresponding landunit', units=' ', & - interpinic_flag='skip', readvar=readvar, data=col%wtlunit) + interpinic_flag='area', readvar=readvar, data=col%wtlunit) call restartvar(ncid=ncid, flag=flag, varname='pfts1d_wtxy', xtype=ncd_double, & dim1name='pft', & long_name='pft weight relative to corresponding gridcell', units='', & - interpinic_flag='skip', readvar=readvar, data=patch%wtgcell) + interpinic_flag='area', readvar=readvar, data=patch%wtgcell) call restartvar(ncid=ncid, flag=flag, varname='pfts1d_wtlnd', xtype=ncd_double, & dim1name='pft', & long_name='pft weight relative to corresponding landunit', units='', & - interpinic_flag='skip', readvar=readvar, data=patch%wtlunit) + interpinic_flag='area', readvar=readvar, data=patch%wtlunit) call restartvar(ncid=ncid, flag=flag, varname='pfts1d_wtcol', xtype=ncd_double, & dim1name='pft', & long_name='pft weight relative to corresponding column', units='', & - interpinic_flag='skip', readvar=readvar, data=patch%wtcol) + interpinic_flag='area', readvar=readvar, data=patch%wtcol) ! Snow column variables diff --git a/src/utils/restUtilMod.F90.in b/src/utils/restUtilMod.F90.in index fa7337b5f6c..c3e8bafc9c0 100644 --- a/src/utils/restUtilMod.F90.in +++ b/src/utils/restUtilMod.F90.in @@ -11,6 +11,7 @@ module restUtilMod use clm_varctl, only: iulog, nsrest, nsrContinue, nsrBranch use clm_varcon, only: spval, ispval use decompMod, only: bounds_type + use abortutils, only: endrun use ncdio_pio use pio use ncdio_utils, only: find_var_on_file @@ -31,9 +32,19 @@ module restUtilMod module procedure restartvar_2d_{TYPE}_bounds end interface restartvar + ! iflag_interp => interpolate variable + ! iflag_copy => copy variable + ! iflag_skip => skip variable: maintain at cold start value on output file + ! iflag_area => area-related variable: skip if we're taking areas from the surface + ! dataset, interp if we're taking areas from the input finidat file. (The latter only + ! works if we're running at the same resolution with a similar configuration.) This + ! can apply to variables that are not actually areas, but for which we want this same + ! conditional behavior. integer,parameter, public :: iflag_interp = 1 integer,parameter, public :: iflag_copy = 2 integer,parameter, public :: iflag_skip = 3 + integer,parameter, public :: iflag_area = 4 + integer,parameter, public :: iflag_noswitchdim = 0 integer,parameter, public :: iflag_switchdim = 1 @@ -54,6 +65,7 @@ module restUtilMod public :: set_missing_vals_to_constant private :: missing_field_possibly_abort + private :: write_interpinic_flag character(len=*), parameter, private :: sourcefile = & __FILE__ @@ -84,7 +96,7 @@ contains character(len=*) , intent(in) :: varname ! variable name (or colon-delimited list: see above) integer , intent(in) :: xtype ! netcdf data type character(len=*) , intent(in) :: long_name ! long name for variable - character(len=*) , intent(in) :: interpinic_flag ! interpolate variable using interpinic + character(len=*) , intent(in) :: interpinic_flag ! interpolate variable using interpinic ('interp', 'copy', 'skip' or 'area': see iflag meanings defined at the top of this module) {VTYPE} , intent(inout) :: data{DIMSTR} logical , intent(out) :: readvar ! was var read? character(len=*) , intent(in), optional :: units ! long name for variable @@ -132,15 +144,8 @@ contains status = PIO_inq_varid(ncid, trim(my_varname), vardesc) varid = vardesc%varid - if (trim(interpinic_flag) == 'interp') then - status = PIO_put_att(ncid, varid, 'interpinic_flag', iflag_interp) - else if (trim(interpinic_flag) == 'copy') then - status = PIO_put_att(ncid, varid, 'interpinic_flag', iflag_copy) - else if (trim(interpinic_flag) == 'skip') then - status = PIO_put_att(ncid, varid, 'interpinic_flag', iflag_skip) - end if - status = PIO_put_att(ncid, varid, 'interpinic_flag_meanings', & - "1=nearest neighbor, 2=copy directly, 3=skip") + call write_interpinic_flag(ncid=ncid, varid=varid, varname=varname, & + interpinic_flag=interpinic_flag) ! This attribute is written in order to communicate this metadata to initInterp call ncd_putatt(ncid, varid, 'varnames_on_old_files', trim(varname)) @@ -217,7 +222,7 @@ contains character(len=*) , intent(in) :: varname ! variable name (or colon-delimited list: see above) integer , intent(in) :: xtype ! netcdf data type character(len=*) , intent(in) :: long_name ! long name for variable - character(len=*) , intent(in) :: interpinic_flag ! interpolate variable using interpinic + character(len=*) , intent(in) :: interpinic_flag ! interpolate variable using interpinic ('interp', 'copy', 'skip' or 'area': see iflag meanings defined at the top of this module) {VTYPE} , pointer :: data{DIMSTR} logical , intent(inout) :: readvar ! was var read? character(len=*) , intent(in), optional :: dim1name ! dimension name @@ -285,15 +290,8 @@ contains status = PIO_inq_varid(ncid, trim(my_varname), vardesc) varid = vardesc%varid - if (trim(interpinic_flag) == 'interp') then - status = PIO_put_att(ncid, varid, 'interpinic_flag', iflag_interp) - else if (trim(interpinic_flag) == 'copy') then - status = PIO_put_att(ncid, varid, 'interpinic_flag', iflag_copy) - else if (trim(interpinic_flag) == 'skip') then - status = PIO_put_att(ncid, varid, 'interpinic_flag', iflag_skip) - end if - status = PIO_put_att(ncid, varid, 'interpinic_flag_meanings', & - "1=nearest neighbor, 2=copy directly, 3=skip") + call write_interpinic_flag(ncid=ncid, varid=varid, varname=varname, & + interpinic_flag=interpinic_flag) ! This attribute is written in order to communicate this metadata to initInterp call ncd_putatt(ncid, varid, 'varnames_on_old_files', trim(varname)) @@ -381,7 +379,7 @@ contains character(len=*) , intent(in) :: dim2name ! dimension name logical , intent(in) :: switchdim character(len=*) , intent(in) :: long_name ! long name for variable - character(len=*) , intent(in) :: interpinic_flag ! interpolate variable using interpinic + character(len=*) , intent(in) :: interpinic_flag ! interpolate variable using interpinic ('interp', 'copy', 'skip' or 'area': see iflag meanings defined at the top of this module) {VTYPE} , pointer :: data(:,:) ! raw data logical , intent(out) :: readvar ! was var read? integer , intent(in), optional :: lowerb2 @@ -438,15 +436,8 @@ contains varid = vardesc%varid - if (trim(interpinic_flag) == 'interp') then - status = PIO_put_att(ncid, varid, 'interpinic_flag', iflag_interp) - else if (trim(interpinic_flag) == 'copy') then - status = PIO_put_att(ncid, varid, 'interpinic_flag', iflag_copy) - else if (trim(interpinic_flag) == 'skip') then - status = PIO_put_att(ncid, varid, 'interpinic_flag', iflag_skip) - end if - status = PIO_put_att(ncid, varid, 'interpinic_flag_meanings', & - "1=>nearest_neighbor 2=>copy 3=>skip") + call write_interpinic_flag(ncid=ncid, varid=varid, varname=varname, & + interpinic_flag=interpinic_flag) ! This attribute is written in order to communicate this metadata to initInterp call ncd_putatt(ncid, varid, 'varnames_on_old_files', trim(varname)) @@ -669,4 +660,42 @@ contains end if end subroutine missing_field_possibly_abort + !----------------------------------------------------------------------- + subroutine write_interpinic_flag(ncid, varid, varname, interpinic_flag) + ! + ! !DESCRIPTION: + ! Write interpinic_flag metadata for the given variable + ! + ! !ARGUMENTS: + type(file_desc_t) , intent(inout) :: ncid ! netcdf file id + integer , intent(in) :: varid ! variable id + character(len=*) , intent(in) :: varname ! variable name (just used for output on error) + character(len=*) , intent(in) :: interpinic_flag ! interpolate variable using interpinic ('interp', 'copy', 'skip' or 'area': see iflag meanings defined at the top of this module) + ! + ! !LOCAL VARIABLES: + integer :: status ! return error code + + character(len=*), parameter :: subname = 'write_interpinic_flag' + !----------------------------------------------------------------------- + + if (interpinic_flag == 'interp') then + status = PIO_put_att(ncid, varid, 'interpinic_flag', iflag_interp) + else if (interpinic_flag == 'copy') then + status = PIO_put_att(ncid, varid, 'interpinic_flag', iflag_copy) + else if (interpinic_flag == 'skip') then + status = PIO_put_att(ncid, varid, 'interpinic_flag', iflag_skip) + else if (interpinic_flag == 'area') then + status = PIO_put_att(ncid, varid, 'interpinic_flag', iflag_area) + else + write(iulog,*) 'ERROR in restartvar for ', trim(varname) + write(iulog,*) 'Unknown interpinic_flag: ', trim(interpinic_flag) + write(iulog,*) 'Allowed values are: interp, copy, skip, area' + call endrun(msg='Unknown interpinic_flag') + end if + status = PIO_put_att(ncid, varid, 'interpinic_flag_meanings', & + "1=nearest neighbor, 2=copy directly, 3=skip, 4=area") + + end subroutine write_interpinic_flag + + end module restUtilMod From d7d071003c5735185f1ee3781077a48362765976 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Tue, 24 Apr 2018 10:41:37 -0600 Subject: [PATCH 04/19] Improve error messages --- src/init_interp/initInterpMindist.F90 | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/src/init_interp/initInterpMindist.F90 b/src/init_interp/initInterpMindist.F90 index 504affe3d11..4f4dba1194f 100644 --- a/src/init_interp/initInterpMindist.F90 +++ b/src/init_interp/initInterpMindist.F90 @@ -85,6 +85,7 @@ subroutine print_point(this, index, unit) write(unit,*) 'subgrid level, index = ',& this%name, index + write(unit,*) 'lat, lon = ', this%lat(index), ', ', this%lon(index) if (associated(this%ltype)) then write(unit,*) 'ltype: ', this%ltype(index) end if @@ -338,7 +339,9 @@ subroutine set_single_match(begi, endi, bego, endo, activeo, subgridi, subgrido, if (found) then write(iulog,*) subname// & ' ERROR: found multiple input points matching output point' - write(iulog,*) 'For this init_interp mode: for a given output point,' + call subgrido%print_point(no, iulog) + write(iulog,*) ' ' + write(iulog,*) 'For this init_interp_method: for a given output point,' write(iulog,*) 'we expect to find exactly one input point at the same' write(iulog,*) 'location with the same type.' write(iulog,*) '(This most likely indicates a problem in CTSM, such as' @@ -362,9 +365,14 @@ subroutine set_single_match(begi, endi, bego, endo, activeo, subgridi, subgrido, if (activeo(no)) then write(iulog,*) subname// & ' ERROR: cannot find any input points matching output point' - write(iulog,*) 'For this init_interp mode: for a given active output point,' + call subgrido%print_point(no, iulog) + write(iulog,*) ' ' + write(iulog,*) 'For this init_interp_method: for a given active output point,' write(iulog,*) 'we expect to find exactly one input point at the same' write(iulog,*) 'location with the same type.' + write(iulog,*) 'Note that this requires the input and output grids to be identical.' + write(iulog,*) '(If you need to interpolate to a different resolution, then' + write(iulog,*) 'you need to use a different init_interp_method.)' call endrun(msg=subname//' ERROR: cannot find any input points matching output point') end if end if From 7e955cab7d3b9beb335b7c8682b947e4ef23b901 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Tue, 24 Apr 2018 11:03:30 -0600 Subject: [PATCH 05/19] Add some notes --- .../testdefs/testmods_dirs/clm/glcMEC_spunup_1way/README | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/cime_config/testdefs/testmods_dirs/clm/glcMEC_spunup_1way/README b/cime_config/testdefs/testmods_dirs/clm/glcMEC_spunup_1way/README index 062dc84800a..ed34c2a28aa 100644 --- a/cime_config/testdefs/testmods_dirs/clm/glcMEC_spunup_1way/README +++ b/cime_config/testdefs/testmods_dirs/clm/glcMEC_spunup_1way/README @@ -12,7 +12,14 @@ glacier areas are taken from the restart file rather than from CISM. This is partly so that changes in CISM don't affect the test results, but more importantly so that LII tests can pass even if CISM changes. (GLC two-way coupling was off in the generation of the -initial conditions file used here, too.) +initial conditions file used here, too.) (4-24-18: I *think* we'll be +able to turn two-way coupling on once +https://github.com/ESCOMP/ctsm/issues/340 is resolved, though I'm not +positive about that. However, if +https://github.com/ESMCI/cime/issues/2484 is done, then we'd need this +testmods to turn off comparison of cpl initial hist files if we have +two-way coupling on, since areas will differ in initialization between +the two cases.) To update the initial conditions (finidat) file for this test: From 3ac75771ac9e6f4e1d5fe3ac1b5ee7da30988578 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Tue, 24 Apr 2018 15:52:18 -0600 Subject: [PATCH 06/19] Add an LII test covering the new init_interp method It may not be totally necessary to have a system test covering this init_interp_method: between unit tests and inline checks, we have quite a bit of testing of this code. It's important to have at least one LII test to make sure we're interpolating all the fields we're supposed to interpolate, but we don't necessarily need a LII test with this mode of operation as long as we have one with the general mode of operation. However, this test is a good check that every point has a unique type (because this mode of operation will fail if that's not the case); this is something needed even for the 'general' init_interp method to work correctly, but the 'general' mode won't catch problems in this regard. One other reason why this test is useful is to cover the threading directives in the set_single_match routine, since those aren't covered by unit tests. So this test mod should be used in an LII test with threading. --- cime_config/testdefs/testlist_clm.xml | 9 ++++ .../clm/glcMEC_spunup_1way/README | 2 +- .../clm/interp_finidat_areas/README | 41 +++++++++++++++++++ .../interp_finidat_areas/include_user_mods | 1 + .../clm/interp_finidat_areas/user_nl_clm | 1 + 5 files changed, 53 insertions(+), 1 deletion(-) create mode 100644 cime_config/testdefs/testmods_dirs/clm/interp_finidat_areas/README create mode 100644 cime_config/testdefs/testmods_dirs/clm/interp_finidat_areas/include_user_mods create mode 100644 cime_config/testdefs/testmods_dirs/clm/interp_finidat_areas/user_nl_clm diff --git a/cime_config/testdefs/testlist_clm.xml b/cime_config/testdefs/testlist_clm.xml index c8cc29ab512..9f97f220acf 100644 --- a/cime_config/testdefs/testlist_clm.xml +++ b/cime_config/testdefs/testlist_clm.xml @@ -1065,6 +1065,15 @@ + + + + + + + + + diff --git a/cime_config/testdefs/testmods_dirs/clm/glcMEC_spunup_1way/README b/cime_config/testdefs/testmods_dirs/clm/glcMEC_spunup_1way/README index ed34c2a28aa..ec74e38b4ce 100644 --- a/cime_config/testdefs/testmods_dirs/clm/glcMEC_spunup_1way/README +++ b/cime_config/testdefs/testmods_dirs/clm/glcMEC_spunup_1way/README @@ -3,7 +3,7 @@ initial condition interpolation is on (use_init_interp=T) to a case with it's of they are exact. For the interpolated result to match uninterpolation, it needs to be a case that essentially needs no interpolation so it's at the same resolution as the initial condition file (finidat file). When surface datasets are changed, or the land-mask is changed, or an -imporant change is made to model physics (for example where new fields are added to the restart +important change is made to model physics (for example where new fields are added to the restart file) -- you'll need to update the initial conditions file in this test (finidat file in the user_nl_clm file). diff --git a/cime_config/testdefs/testmods_dirs/clm/interp_finidat_areas/README b/cime_config/testdefs/testmods_dirs/clm/interp_finidat_areas/README new file mode 100644 index 00000000000..8773f4ce179 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/interp_finidat_areas/README @@ -0,0 +1,41 @@ +This testmod is designed to test init_interp_method='use_finidat_areas'. + +This should be used in an LII test. So it must be used in a +configuration for which we have a compatible out-of-the-box finidat file +(so that the run with use_init_interp = .false. runs successfully). In +constrast to our standard LII test (which uses glcMEC_spunup_1way), this +one can use a standard CISM2%NOEVOLVE configuration: we do *not* need to +set GLC_TWO_WAY_COUPLING=FALSE; in fact, it's a better test if we have +GLC_TWO_WAY_COUPLING=TRUE: with this mode of operation, areas should +match between the two runs. + +It may not be totally necessary to have a system test covering this +init_interp_method: between unit tests and inline checks, we have quite +a bit of testing of this code. It's important to have at least one LII +test to make sure we're interpolating all the fields we're supposed to +interpolate, but we don't necessarily need a LII test with this mode of +operation as long as we have one with the general mode of operation. + +However, this test is a good check that every point has a unique type +(because this mode of operation will fail if that's not the case); this +is something needed even for the 'general' init_interp method to work +correctly, but the 'general' mode won't catch problems in this regard. + +One other reason why this test is useful is to cover the threading +directives in the set_single_match routine, since those aren't covered +by unit tests. So this test mod should be used in an LII test with +threading. + +To update the initial conditions (finidat) file for this test: + +(1) Run the LII test; the 'base' case should run to completion even if the +no_interp test fails. (If the 'base' case fails, you may need to retry +with init_interp_method='general'.) + +(2) Copy the finidat_interp_dest.nc file from the 'base' case to the inputdata +space. Rename this to be similar to the out-of-the-box finidat file +currently used by this test, but with a new creation date. + +(3) Update namelist defaults to point to the new finidat file. If +updating the out-of-the-box file is not desired, then you could instead +point to this new finidat file with a user_nl_clm file in this testmod. diff --git a/cime_config/testdefs/testmods_dirs/clm/interp_finidat_areas/include_user_mods b/cime_config/testdefs/testmods_dirs/clm/interp_finidat_areas/include_user_mods new file mode 100644 index 00000000000..fe0e18cf889 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/interp_finidat_areas/include_user_mods @@ -0,0 +1 @@ +../default diff --git a/cime_config/testdefs/testmods_dirs/clm/interp_finidat_areas/user_nl_clm b/cime_config/testdefs/testmods_dirs/clm/interp_finidat_areas/user_nl_clm new file mode 100644 index 00000000000..53b8e98df1b --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/interp_finidat_areas/user_nl_clm @@ -0,0 +1 @@ +init_interp_method='use_finidat_areas' From 7eaa12f43804e00ee35c7a04ee0f699c57838728 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Tue, 24 Apr 2018 19:30:58 -0600 Subject: [PATCH 07/19] Point to new initial conditions file Point to a new initial conditions file that is compatible with the current memory allocation. This is needed for the new LII test (LII_D_P360x2_Ld1.f09_g16_gl4.I1850Clm50BgcCrop.cheyenne_intel.clm-interp_finidat_areas). It was created using the finidat_interp_dest file created by that test. --- bld/namelist_files/namelist_defaults_clm4_5.xml | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/bld/namelist_files/namelist_defaults_clm4_5.xml b/bld/namelist_files/namelist_defaults_clm4_5.xml index 6a165a9f8fc..9f97863b7e6 100644 --- a/bld/namelist_files/namelist_defaults_clm4_5.xml +++ b/bld/namelist_files/namelist_defaults_clm4_5.xml @@ -515,6 +515,10 @@ attributes from the config_cache.xml file (with keys converted to upper-case). >lnd/clm2/initdata_map/clmi.I1850Clm45BgcGs.0901-01-01.0.9x1.25_gx1v6_simyr1850_c180204.nc + lnd/clm2/initdata_map/clmi.I1850Clm50BgcCrop.1366-01-01.0.9x1.25_gx1v6_simyr1850_c171213.nc +>lnd/clm2/initdata_map/clmi.I1850Clm50BgcCrop.1366-01-01.0.9x1.25_gx1v6_simyr1850_c180424.nc From d414d83ea0851025e9f2a8bf28e402a3b29e200d Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Tue, 24 Apr 2018 20:03:25 -0600 Subject: [PATCH 08/19] Add a note --- bld/CLMBuildNamelist.pm | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/bld/CLMBuildNamelist.pm b/bld/CLMBuildNamelist.pm index 568e663b6c1..4fe4d3fd744 100755 --- a/bld/CLMBuildNamelist.pm +++ b/bld/CLMBuildNamelist.pm @@ -3734,6 +3734,11 @@ sub setup_logic_initinterp { my ($opts, $nl_flags, $definition, $defaults, $nl, $physv) = @_; if ($physv->as_long() >= $physv->as_long("clm4_5")) { + # Ideally we would only let this be added to the namelist if + # use_init_interp = .true. (e.g., as is done for variables that + # can only be set if repartition_rain_snow is true in + # setup_logic_atm_forcing). However, that doesn't work with the + # way we have set up the LII test that tests this option. add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'init_interp_method'); } } From fd2735b3e5943dedceabeede899802e4b02c4536 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Tue, 24 Apr 2018 20:19:50 -0600 Subject: [PATCH 09/19] Add some documentation --- .../namelist_definition_clm4_5.xml | 27 +++++++++++++++++-- src/main/subgridWeightsMod.F90 | 8 ++++++ 2 files changed, 33 insertions(+), 2 deletions(-) diff --git a/bld/namelist_files/namelist_definition_clm4_5.xml b/bld/namelist_files/namelist_definition_clm4_5.xml index 92614eced90..38984b2374f 100644 --- a/bld/namelist_files/namelist_definition_clm4_5.xml +++ b/bld/namelist_files/namelist_definition_clm4_5.xml @@ -2419,9 +2419,32 @@ flag. So interpolation from non-crop to crop cases can be done without setting t -Method to use for init_interp +Method to use for init_interp. Only applies when use_init_interp = .true. + +'general': The general-purpose method that can be used when changing +grids, configurations, etc. This starts off with subgrid areas taken +from the surface dataset. + +'use_finidat_areas': This starts off with subgrid areas taken from the +input finidat file. This is important to achieve bit-for-bit results in +a coupled case (where areas in initialization impact initial fields sent +to the atmosphere). For this method to work, the input finidat file +needs to be at the same resolution as the current configuration. So this +is a less general form of init_interp. However, it can be used when +transitioning from a spinup run to a transient run, or in other cases +where the only difference is in internal memory allocation. In order to +catch possible problems, this uses a different algorithm for finding the +input point for each output point, which ensures that each active output +point is associated with exactly one input point with the same latitude, +longitude and type. This method requires (a) the same grid for input and +output, within roundoff; (b) any non-zero-weight point in the input must +have memory allocated for it in this grid cell in the output (this will +be satisfied if the point is non-zero-weight on the surface dataset or +if it's a point for which we allocate memory even for zero-weight +points); (c) any active point in the output (based on the surface +dataset and rules for determining active points) must have a matching +point in this grid cell in the input. -Only applies when use_init_interp = .true. diff --git a/src/main/subgridWeightsMod.F90 b/src/main/subgridWeightsMod.F90 index 6a59f376aaf..b2c882a478c 100644 --- a/src/main/subgridWeightsMod.F90 +++ b/src/main/subgridWeightsMod.F90 @@ -673,6 +673,14 @@ subroutine check_weights (bounds, active_only) deallocate(sumwtcol, sumwtlunit, sumwtgcell) if (error_found) then + write(iulog,*) ' ' + write(iulog,*) 'If you are seeing this message at the beginning of a run with' + write(iulog,*) 'use_init_interp = .true. and init_interp_method = "use_finidat_areas",' + write(iulog,*) 'and you are seeing weights less than 1, then a likely cause is:' + write(iulog,*) 'For the above-mentioned grid cell(s):' + write(iulog,*) 'The matching input grid cell had some non-zero-weight subgrid type' + write(iulog,*) 'that is not present in memory in the new run.' + write(iulog,*) ' ' call endrun(msg=errMsg(sourcefile, __LINE__)) end if From e34b7072d45aa800eb0cc1a1a839e8505428f28d Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Thu, 26 Apr 2018 10:17:54 -0600 Subject: [PATCH 10/19] Add a usermods directory that turns on the ICE_MODEL_FRACTION field Also, change ICE_MODEL_FRACTION to time-average rather than instantaneous: I think this will more often do the right thing for various history output frequencies. For example, if someone adds this to an every-10-year hist file, they'd probably want the 10-year-average ICE_MODEL_FRACTION, not the value at the end of the 10-year period. Also: The thing that originally got me worried about having it instantaneous was not being completely sure whether the value written at the end of the year would apply to the current year or the previous year. I think it would apply to the previous year (which I think is what you'd want), but I felt it's safer to just write the time-averaged value (where this year boundary problem would, at worst, lead to errors of 1/17520 for a 30-min time step, rather than leading to a complete year offset -- i.e., if we average just the year-boundary value into the wrong year). Fixes #345 --- cime_config/testdefs/testmods_dirs/clm/cmip6/README | 6 ++++++ .../testdefs/testmods_dirs/clm/cmip6/include_user_mods | 2 +- .../usermods_dirs/cmip6_evolving_icesheet/include_user_mods | 1 + .../usermods_dirs/cmip6_evolving_icesheet/user_nl_clm | 1 + src/main/glc2lndMod.F90 | 2 +- 5 files changed, 10 insertions(+), 2 deletions(-) create mode 100644 cime_config/testdefs/testmods_dirs/clm/cmip6/README create mode 100644 cime_config/usermods_dirs/cmip6_evolving_icesheet/include_user_mods create mode 100644 cime_config/usermods_dirs/cmip6_evolving_icesheet/user_nl_clm diff --git a/cime_config/testdefs/testmods_dirs/clm/cmip6/README b/cime_config/testdefs/testmods_dirs/clm/cmip6/README new file mode 100644 index 00000000000..40b90f9ce07 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/cmip6/README @@ -0,0 +1,6 @@ +Note that this testmod directory includes cmip6_evolving_icesheet rather +than simply cmip6: cmip6_evolving_icesheet includes cmip6 and adds to +it. So, rather than having tests of the cmip6 usermods and separate +tests of cmip6_evolving_icesheet, it's more efficient to just have tests +of the cmip6_evolving_icesheet usermods directory, which then also cover +the cmip6 usermods directory. diff --git a/cime_config/testdefs/testmods_dirs/clm/cmip6/include_user_mods b/cime_config/testdefs/testmods_dirs/clm/cmip6/include_user_mods index 0caefd70b41..7e8fa5ec833 100644 --- a/cime_config/testdefs/testmods_dirs/clm/cmip6/include_user_mods +++ b/cime_config/testdefs/testmods_dirs/clm/cmip6/include_user_mods @@ -1 +1 @@ -../../../../usermods_dirs/cmip6 +../../../../usermods_dirs/cmip6_evolving_icesheet diff --git a/cime_config/usermods_dirs/cmip6_evolving_icesheet/include_user_mods b/cime_config/usermods_dirs/cmip6_evolving_icesheet/include_user_mods new file mode 100644 index 00000000000..d8d7dacfae4 --- /dev/null +++ b/cime_config/usermods_dirs/cmip6_evolving_icesheet/include_user_mods @@ -0,0 +1 @@ +../cmip6 diff --git a/cime_config/usermods_dirs/cmip6_evolving_icesheet/user_nl_clm b/cime_config/usermods_dirs/cmip6_evolving_icesheet/user_nl_clm new file mode 100644 index 00000000000..193f0624e7e --- /dev/null +++ b/cime_config/usermods_dirs/cmip6_evolving_icesheet/user_nl_clm @@ -0,0 +1 @@ +hist_fincl4 += 'ICE_MODEL_FRACTION' diff --git a/src/main/glc2lndMod.F90 b/src/main/glc2lndMod.F90 index 9df55b5e70f..c7e4402ac9b 100644 --- a/src/main/glc2lndMod.F90 +++ b/src/main/glc2lndMod.F90 @@ -170,7 +170,7 @@ subroutine InitHistory(this, bounds) this%icemask_grc(begg:endg) = spval call hist_addfld1d (fname='ICE_MODEL_FRACTION', units='unitless', & - avgflag='I', long_name='Ice sheet model fractional coverage', & + avgflag='A', long_name='Ice sheet model fractional coverage', & ptr_gcell=this%icemask_grc, default='inactive') end subroutine InitHistory From 29d9b4652ee6f9d4eb072747bd67b5f61327cd1f Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Fri, 27 Apr 2018 14:49:07 -0600 Subject: [PATCH 11/19] Tweak documentation of init_interp_method --- .../namelist_definition_clm4_5.xml | 38 ++++++++++--------- 1 file changed, 20 insertions(+), 18 deletions(-) diff --git a/bld/namelist_files/namelist_definition_clm4_5.xml b/bld/namelist_files/namelist_definition_clm4_5.xml index 38984b2374f..b72089197d8 100644 --- a/bld/namelist_files/namelist_definition_clm4_5.xml +++ b/bld/namelist_files/namelist_definition_clm4_5.xml @@ -2426,24 +2426,26 @@ grids, configurations, etc. This starts off with subgrid areas taken from the surface dataset. 'use_finidat_areas': This starts off with subgrid areas taken from the -input finidat file. This is important to achieve bit-for-bit results in -a coupled case (where areas in initialization impact initial fields sent -to the atmosphere). For this method to work, the input finidat file -needs to be at the same resolution as the current configuration. So this -is a less general form of init_interp. However, it can be used when -transitioning from a spinup run to a transient run, or in other cases -where the only difference is in internal memory allocation. In order to -catch possible problems, this uses a different algorithm for finding the -input point for each output point, which ensures that each active output -point is associated with exactly one input point with the same latitude, -longitude and type. This method requires (a) the same grid for input and -output, within roundoff; (b) any non-zero-weight point in the input must -have memory allocated for it in this grid cell in the output (this will -be satisfied if the point is non-zero-weight on the surface dataset or -if it's a point for which we allocate memory even for zero-weight -points); (c) any active point in the output (based on the surface -dataset and rules for determining active points) must have a matching -point in this grid cell in the input. +input finidat file. This is needed to achieve bit-for-bit results in a +coupled case (where areas in initialization impact initial fields sent +to the atmosphere) (but using the 'general' method will typically have +only a very minor impact on results in this case). For this method to +work, the input finidat file needs to be at the same resolution as the +current configuration. So this is a less general form of +init_interp. However, it can be used when transitioning from a spinup +run to a transient run, or in other cases where the only difference is +in internal memory allocation. In order to catch possible problems, this +uses a different algorithm for finding the input point for each output +point, which ensures that each active output point is associated with +exactly one input point with the same latitude, longitude and type. This +method requires (a) the same grid for input and output, within roundoff; +(b) any non-zero-weight point in the input must have memory allocated +for it in this grid cell in the output (this will be satisfied if the +point is non-zero-weight on the surface dataset or if it's a point for +which we allocate memory even for zero-weight points); (c) any active +point in the output (based on the surface dataset and rules for +determining active points) must have a matching point in this grid cell +in the input. From a9caf9925340340ef515e779aa0be272e10d013c Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Tue, 8 May 2018 14:27:20 -0600 Subject: [PATCH 12/19] Don't allow user to set init_interp_method if use_init_interp is false Erik wants this check in place. For this error check to work, I needed to rework the LII test that tests this option: Rather than adding init_interp_method via a testmod, I needed to make a new test type that adds this option. --- bld/CLMBuildNamelist.pm | 14 ++-- cime_config/SystemTests/lii2finidatareas.py | 76 +++++++++++++++++++ cime_config/config_tests.xml | 10 +++ cime_config/testdefs/testlist_clm.xml | 4 +- .../clm/interp_finidat_areas/README | 41 ---------- .../interp_finidat_areas/include_user_mods | 1 - .../clm/interp_finidat_areas/user_nl_clm | 1 - 7 files changed, 96 insertions(+), 51 deletions(-) create mode 100644 cime_config/SystemTests/lii2finidatareas.py delete mode 100644 cime_config/testdefs/testmods_dirs/clm/interp_finidat_areas/README delete mode 100644 cime_config/testdefs/testmods_dirs/clm/interp_finidat_areas/include_user_mods delete mode 100644 cime_config/testdefs/testmods_dirs/clm/interp_finidat_areas/user_nl_clm diff --git a/bld/CLMBuildNamelist.pm b/bld/CLMBuildNamelist.pm index 4fe4d3fd744..3595ec74615 100755 --- a/bld/CLMBuildNamelist.pm +++ b/bld/CLMBuildNamelist.pm @@ -3734,12 +3734,14 @@ sub setup_logic_initinterp { my ($opts, $nl_flags, $definition, $defaults, $nl, $physv) = @_; if ($physv->as_long() >= $physv->as_long("clm4_5")) { - # Ideally we would only let this be added to the namelist if - # use_init_interp = .true. (e.g., as is done for variables that - # can only be set if repartition_rain_snow is true in - # setup_logic_atm_forcing). However, that doesn't work with the - # way we have set up the LII test that tests this option. - add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'init_interp_method'); + my $var = 'init_interp_method'; + if ( &value_is_true($nl->get_value("use_init_interp"))) { + add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, $var); + } else { + if (defined($nl->get_value($var))) { + $log->fatal_error("$var can only be set if use_init_interp is true"); + } + } } } diff --git a/cime_config/SystemTests/lii2finidatareas.py b/cime_config/SystemTests/lii2finidatareas.py new file mode 100644 index 00000000000..21ba082fffe --- /dev/null +++ b/cime_config/SystemTests/lii2finidatareas.py @@ -0,0 +1,76 @@ +""" +Implementation of the LII2FINIDATAREAS test. + +This is similar to the LII test, but tests init_interp with mode +'use_finidat_areas'. + +As with the standard LII test, this must be used in a configuration for +which we have a compatible out-of-the-box finidat file (so that the run +with use_init_interp = .false. runs successfully). In constrast to our +standard LII test (which uses glcMEC_spunup_1way), this one can use a +standard CISM2%NOEVOLVE configuration: we do *not* need to set +GLC_TWO_WAY_COUPLING=FALSE; in fact, it's a better test if we have +GLC_TWO_WAY_COUPLING=TRUE: with this mode of operation, areas should +match between the two runs. + +It may not be totally necessary to have a system test covering this +init_interp_method: between unit tests and inline checks, we have quite +a bit of testing of this code. It's important to have at least one LII +test to make sure we're interpolating all the fields we're supposed to +interpolate, but we don't necessarily need a LII test with this mode of +operation as long as we have one with the general mode of operation. + +However, this test is a good check that every point has a unique type +(because this mode of operation will fail if that's not the case); this +is something needed even for the 'general' init_interp method to work +correctly, but the 'general' mode won't catch problems in this regard. + +One other reason why this test is useful is to cover the threading +directives in the set_single_match routine, since those aren't covered +by unit tests. So this test mod should be used in a test with threading. + +To update the initial conditions (finidat) file for this test: + +(1) Run the test; the 'base' case should run to completion even if the +no_interp test fails. (If the 'base' case fails, you may need to retry +with init_interp_method='general'.) + +(2) Copy the finidat_interp_dest.nc file from the 'base' case to the inputdata +space. Rename this to be similar to the out-of-the-box finidat file +currently used by this test, but with a new creation date. + +(3) Update namelist defaults to point to the new finidat file. If +updating the out-of-the-box file is not desired, then you could instead +point to this new finidat file with a user_nl_clm file in this testmod. +""" + +from CIME.XML.standard_module_setup import * +from CIME.SystemTests.test_utils.user_nl_utils import append_to_user_nl_files + +# We can import lii directly because the SystemTests directory has been +# added to sys.path. +# +# A cleaner and more general way to accomplish this would be: In cime: +# Rather than adding the SystemTests directory to sys.path and then +# importing an individual module from there: instead, allow each +# component to have a COMPNAME_pylib directory within its cime_config, +# and then have cime add each component's cime_config directory to +# sys.path and then do something like: +# import_module("COMPNAME_pylib.SystemTests.TESTNAME"). Then, for +# example, ctsm could access its own modules via "import +# ctsm_pylib.foo", or (in this case) "from ctsm_pylib.SystemTests.lii +# import LII". +from lii import LII + +logger = logging.getLogger(__name__) + +class LII2FINIDATAREAS(LII): + + def __init__(self, case): + super(LII2FINIDATAREAS, self).__init__(case) + + def _case_one_setup(self): + super(LII2FINIDATAREAS, self)._case_one_setup() + append_to_user_nl_files(caseroot = self._get_caseroot(), + component = "clm", + contents = "init_interp_method = 'use_finidat_areas'") diff --git a/cime_config/config_tests.xml b/cime_config/config_tests.xml index c829c8421f6..33f8b4a692a 100644 --- a/cime_config/config_tests.xml +++ b/cime_config/config_tests.xml @@ -28,6 +28,16 @@ SSP smoke CLM spinup test (only valid for CLM compsets with CLM45) $STOP_N + + CLM initial condition interpolation test using finidat_areas (requires configuration with non-blank finidat) + 1 + FALSE + FALSE + never + $STOP_OPTION + $STOP_N + + CLM test: Verify that adding virtual glacier columns doesn't change answers 1 diff --git a/cime_config/testdefs/testlist_clm.xml b/cime_config/testdefs/testlist_clm.xml index b08009587af..defc564ed77 100644 --- a/cime_config/testdefs/testlist_clm.xml +++ b/cime_config/testdefs/testlist_clm.xml @@ -1065,13 +1065,13 @@ - + - + diff --git a/cime_config/testdefs/testmods_dirs/clm/interp_finidat_areas/README b/cime_config/testdefs/testmods_dirs/clm/interp_finidat_areas/README deleted file mode 100644 index 8773f4ce179..00000000000 --- a/cime_config/testdefs/testmods_dirs/clm/interp_finidat_areas/README +++ /dev/null @@ -1,41 +0,0 @@ -This testmod is designed to test init_interp_method='use_finidat_areas'. - -This should be used in an LII test. So it must be used in a -configuration for which we have a compatible out-of-the-box finidat file -(so that the run with use_init_interp = .false. runs successfully). In -constrast to our standard LII test (which uses glcMEC_spunup_1way), this -one can use a standard CISM2%NOEVOLVE configuration: we do *not* need to -set GLC_TWO_WAY_COUPLING=FALSE; in fact, it's a better test if we have -GLC_TWO_WAY_COUPLING=TRUE: with this mode of operation, areas should -match between the two runs. - -It may not be totally necessary to have a system test covering this -init_interp_method: between unit tests and inline checks, we have quite -a bit of testing of this code. It's important to have at least one LII -test to make sure we're interpolating all the fields we're supposed to -interpolate, but we don't necessarily need a LII test with this mode of -operation as long as we have one with the general mode of operation. - -However, this test is a good check that every point has a unique type -(because this mode of operation will fail if that's not the case); this -is something needed even for the 'general' init_interp method to work -correctly, but the 'general' mode won't catch problems in this regard. - -One other reason why this test is useful is to cover the threading -directives in the set_single_match routine, since those aren't covered -by unit tests. So this test mod should be used in an LII test with -threading. - -To update the initial conditions (finidat) file for this test: - -(1) Run the LII test; the 'base' case should run to completion even if the -no_interp test fails. (If the 'base' case fails, you may need to retry -with init_interp_method='general'.) - -(2) Copy the finidat_interp_dest.nc file from the 'base' case to the inputdata -space. Rename this to be similar to the out-of-the-box finidat file -currently used by this test, but with a new creation date. - -(3) Update namelist defaults to point to the new finidat file. If -updating the out-of-the-box file is not desired, then you could instead -point to this new finidat file with a user_nl_clm file in this testmod. diff --git a/cime_config/testdefs/testmods_dirs/clm/interp_finidat_areas/include_user_mods b/cime_config/testdefs/testmods_dirs/clm/interp_finidat_areas/include_user_mods deleted file mode 100644 index fe0e18cf889..00000000000 --- a/cime_config/testdefs/testmods_dirs/clm/interp_finidat_areas/include_user_mods +++ /dev/null @@ -1 +0,0 @@ -../default diff --git a/cime_config/testdefs/testmods_dirs/clm/interp_finidat_areas/user_nl_clm b/cime_config/testdefs/testmods_dirs/clm/interp_finidat_areas/user_nl_clm deleted file mode 100644 index 53b8e98df1b..00000000000 --- a/cime_config/testdefs/testmods_dirs/clm/interp_finidat_areas/user_nl_clm +++ /dev/null @@ -1 +0,0 @@ -init_interp_method='use_finidat_areas' From 0fb1de4b78f5bf3c84518147707c18107dfd36d1 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Tue, 8 May 2018 17:10:58 -0600 Subject: [PATCH 13/19] Delete the interp testmods and associated tests Now that we're doing interpolation pretty much all the time, we are covering the cases that these originally tested. Two examples are: - Crop to non-crop: covered by ERP_Ld5.f10_f10_musgs.I1850Clm50Bgc.cheyenne_gnu.clm-default1850 (uses finidat clmi.I1850Clm50BgcCrop.1366-01-01.0.9x1.25_gx1v6_simyr1850_c180424.nc) - Non-crop to crop: covered by SMS_Ld5.f10_f10_musgs.I1850Clm45BgcCrop.cheyenne_gnu.clm-crop1850 (uses finidat clmi.I1850Clm45BgcGs.0901-01-01.0.9x1.25_gx1v6_simyr1850_c180204.nc) In addition, we cover the mindist logic via unit test. What's not covered by unit test is determining if all variables are initialized reasonably - particularly important for the non-crop to crop case. But I'm not sure how we can really be confident of that even if we have a system test covering this (other than the standard way of seeing that baseline comparisons haven't changed... but the issue is: if new behavior is introduced that changes answers for crops, then we could expect baseline comparisons to change, and so could miss some issue). That is to say: I'm not sure it's worth having system tests that explicitly test these cases: it's probably enough to have the unit tests plus incidental coverage via existing system tests. Resolves #369 --- cime_config/testdefs/testlist_clm.xml | 22 ++----------------- .../testmods_dirs/clm/interp_f19_crop/README | 4 ---- .../clm/interp_f19_crop/include_user_mods | 1 - .../clm/interp_f19_crop/user_nl_clm | 1 - .../clm/interp_f19_crop/user_nl_cpl | 2 -- .../clm/interp_f19_noncrop/README | 12 ---------- .../clm/interp_f19_noncrop/include_user_mods | 1 - .../clm/interp_f19_noncrop/user_nl_clm | 2 -- .../interp_f19_noncrop1850/include_user_mods | 1 - .../clm/interp_f19_noncrop1850/user_nl_cpl | 2 -- 10 files changed, 2 insertions(+), 46 deletions(-) delete mode 100644 cime_config/testdefs/testmods_dirs/clm/interp_f19_crop/README delete mode 100644 cime_config/testdefs/testmods_dirs/clm/interp_f19_crop/include_user_mods delete mode 100644 cime_config/testdefs/testmods_dirs/clm/interp_f19_crop/user_nl_clm delete mode 100644 cime_config/testdefs/testmods_dirs/clm/interp_f19_crop/user_nl_cpl delete mode 100644 cime_config/testdefs/testmods_dirs/clm/interp_f19_noncrop/README delete mode 100644 cime_config/testdefs/testmods_dirs/clm/interp_f19_noncrop/include_user_mods delete mode 100644 cime_config/testdefs/testmods_dirs/clm/interp_f19_noncrop/user_nl_clm delete mode 100644 cime_config/testdefs/testmods_dirs/clm/interp_f19_noncrop1850/include_user_mods delete mode 100644 cime_config/testdefs/testmods_dirs/clm/interp_f19_noncrop1850/user_nl_cpl diff --git a/cime_config/testdefs/testlist_clm.xml b/cime_config/testdefs/testlist_clm.xml index defc564ed77..a6c3b3b5476 100644 --- a/cime_config/testdefs/testlist_clm.xml +++ b/cime_config/testdefs/testlist_clm.xml @@ -1115,13 +1115,13 @@ - + - + @@ -1240,15 +1240,6 @@ - - - - - - - - - @@ -1269,15 +1260,6 @@ - - - - - - - - - diff --git a/cime_config/testdefs/testmods_dirs/clm/interp_f19_crop/README b/cime_config/testdefs/testmods_dirs/clm/interp_f19_crop/README deleted file mode 100644 index 9b2a8378d7a..00000000000 --- a/cime_config/testdefs/testmods_dirs/clm/interp_f19_crop/README +++ /dev/null @@ -1,4 +0,0 @@ -This testmods directory allows testing the interpolation of initial conditions -from one resolution (or configuration) to another, using the out-of-the-box -initial conditions file for f19 with crop. - diff --git a/cime_config/testdefs/testmods_dirs/clm/interp_f19_crop/include_user_mods b/cime_config/testdefs/testmods_dirs/clm/interp_f19_crop/include_user_mods deleted file mode 100644 index fe0e18cf889..00000000000 --- a/cime_config/testdefs/testmods_dirs/clm/interp_f19_crop/include_user_mods +++ /dev/null @@ -1 +0,0 @@ -../default diff --git a/cime_config/testdefs/testmods_dirs/clm/interp_f19_crop/user_nl_clm b/cime_config/testdefs/testmods_dirs/clm/interp_f19_crop/user_nl_clm deleted file mode 100644 index 05c8a39983b..00000000000 --- a/cime_config/testdefs/testmods_dirs/clm/interp_f19_crop/user_nl_clm +++ /dev/null @@ -1 +0,0 @@ -use_init_interp = .true. diff --git a/cime_config/testdefs/testmods_dirs/clm/interp_f19_crop/user_nl_cpl b/cime_config/testdefs/testmods_dirs/clm/interp_f19_crop/user_nl_cpl deleted file mode 100644 index 542adaaf265..00000000000 --- a/cime_config/testdefs/testmods_dirs/clm/interp_f19_crop/user_nl_cpl +++ /dev/null @@ -1,2 +0,0 @@ -orb_iyear = 1850 -orb_iyear_align = 1850 diff --git a/cime_config/testdefs/testmods_dirs/clm/interp_f19_noncrop/README b/cime_config/testdefs/testmods_dirs/clm/interp_f19_noncrop/README deleted file mode 100644 index 4578d53e8c7..00000000000 --- a/cime_config/testdefs/testmods_dirs/clm/interp_f19_noncrop/README +++ /dev/null @@ -1,12 +0,0 @@ -This testmods directory allows testing the interpolation of initial conditions -from one resolution (or configuration) to another, using the out-of-the-box -initial conditions file for f19 without prognostic crop. - -Note: if we stop supporting out-of-the-box f19 initial conditions files, this -could be changed to f09 (or whatever resolution we provide initial conditions -for). - -Note: If this test fails, it may be because you need to update the -initial conditions file it uses (finidat). finidat should be updated -to match the out-of-the-box initial conditions file for f19 without -prognostic crop. diff --git a/cime_config/testdefs/testmods_dirs/clm/interp_f19_noncrop/include_user_mods b/cime_config/testdefs/testmods_dirs/clm/interp_f19_noncrop/include_user_mods deleted file mode 100644 index fe0e18cf889..00000000000 --- a/cime_config/testdefs/testmods_dirs/clm/interp_f19_noncrop/include_user_mods +++ /dev/null @@ -1 +0,0 @@ -../default diff --git a/cime_config/testdefs/testmods_dirs/clm/interp_f19_noncrop/user_nl_clm b/cime_config/testdefs/testmods_dirs/clm/interp_f19_noncrop/user_nl_clm deleted file mode 100644 index 928df99e541..00000000000 --- a/cime_config/testdefs/testmods_dirs/clm/interp_f19_noncrop/user_nl_clm +++ /dev/null @@ -1,2 +0,0 @@ -finidat = '$DIN_LOC_ROOT/lnd/clm2/initdata_map/clmi.IGM1850GSWCLM50BGCCROP.0481-01-01.1.9x2.5_gx1v6_gl5_simyr1850_c170419.nc' -use_init_interp = .true. diff --git a/cime_config/testdefs/testmods_dirs/clm/interp_f19_noncrop1850/include_user_mods b/cime_config/testdefs/testmods_dirs/clm/interp_f19_noncrop1850/include_user_mods deleted file mode 100644 index b786856cdce..00000000000 --- a/cime_config/testdefs/testmods_dirs/clm/interp_f19_noncrop1850/include_user_mods +++ /dev/null @@ -1 +0,0 @@ -../interp_f19_noncrop diff --git a/cime_config/testdefs/testmods_dirs/clm/interp_f19_noncrop1850/user_nl_cpl b/cime_config/testdefs/testmods_dirs/clm/interp_f19_noncrop1850/user_nl_cpl deleted file mode 100644 index 542adaaf265..00000000000 --- a/cime_config/testdefs/testmods_dirs/clm/interp_f19_noncrop1850/user_nl_cpl +++ /dev/null @@ -1,2 +0,0 @@ -orb_iyear = 1850 -orb_iyear_align = 1850 From 21ca821fe65e5654e940f3cee53c798a24de26b6 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Tue, 8 May 2018 19:45:28 -0600 Subject: [PATCH 14/19] Move comment to the correct finidat file --- bld/namelist_files/namelist_defaults_clm4_5.xml | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/bld/namelist_files/namelist_defaults_clm4_5.xml b/bld/namelist_files/namelist_defaults_clm4_5.xml index 7a80c369dbe..e39f041fb39 100644 --- a/bld/namelist_files/namelist_defaults_clm4_5.xml +++ b/bld/namelist_files/namelist_defaults_clm4_5.xml @@ -517,10 +517,6 @@ attributes from the config_cache.xml file (with keys converted to upper-case). >lnd/clm2/initdata_map/clmi.I1850Clm45BgcGs.0901-01-01.0.9x1.25_gx1v6_simyr1850_c180204.nc - + Date: Wed, 9 May 2018 11:59:22 -0600 Subject: [PATCH 15/19] Only read init_interp namelist if use_init_interp is true Otherwise, when use_init_interp is false, we get an error that init_interp_method has an unknown value --- src/main/controlMod.F90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/main/controlMod.F90 b/src/main/controlMod.F90 index 9cc442d026b..ed935937825 100644 --- a/src/main/controlMod.F90 +++ b/src/main/controlMod.F90 @@ -438,7 +438,9 @@ subroutine control_init( ) ! Read in other namelists for other modules ! ---------------------------------------------------------------------- - call initInterp_readnl( NLFilename ) + if (use_init_interp) then + call initInterp_readnl( NLFilename ) + end if !I call init_hydrology to set up default hydrology sub-module methods. !For future version, I suggest to put the following two calls inside their From c061ef6187c37f669cb7ebea5604db818c59ee6b Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Wed, 9 May 2018 12:29:52 -0600 Subject: [PATCH 16/19] Revert "Point to new initial conditions file" This reverts commit 21ca821fe65e5654e940f3cee53c798a24de26b6. This reverts commit 7eaa12f43804e00ee35c7a04ee0f699c57838728. There were answer changes for the new initial conditions file. One problem with the new finidat file is that it doesn't have C isotope information on it... that's a problem, but not the (only) problem with these baseline failures, I think. Another difference is that two variables have changed from double to int. I'm not sure if these completely explain the problem, but I'm going to back out the change to finidat, just using the new file for the test that needs it. --- bld/namelist_files/namelist_defaults_clm4_5.xml | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/bld/namelist_files/namelist_defaults_clm4_5.xml b/bld/namelist_files/namelist_defaults_clm4_5.xml index e39f041fb39..42f4048e9b5 100644 --- a/bld/namelist_files/namelist_defaults_clm4_5.xml +++ b/bld/namelist_files/namelist_defaults_clm4_5.xml @@ -540,15 +540,11 @@ attributes from the config_cache.xml file (with keys converted to upper-case). - lnd/clm2/initdata_map/clmi.I1850Clm50BgcCrop.1366-01-01.0.9x1.25_gx1v6_simyr1850_c180424.nc +>lnd/clm2/initdata_map/clmi.I1850Clm50BgcCrop.1366-01-01.0.9x1.25_gx1v6_simyr1850_c171213.nc From c39540615ebc9185fac9e2b5bbe6c9f2ba992b17 Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Wed, 9 May 2018 12:41:35 -0600 Subject: [PATCH 17/19] Change new LII test to use a testmods that points to compatible finidat This is needed now that we've reverted to having the out-of-the-box be incompatible with the current configuration. --- cime_config/testdefs/testlist_clm.xml | 2 +- .../clm/compatible_finidat_f09/README | 23 +++++++++++++++++++ .../compatible_finidat_f09/include_user_mods | 1 + .../clm/compatible_finidat_f09/user_nl_clm | 1 + 4 files changed, 26 insertions(+), 1 deletion(-) create mode 100644 cime_config/testdefs/testmods_dirs/clm/compatible_finidat_f09/README create mode 100644 cime_config/testdefs/testmods_dirs/clm/compatible_finidat_f09/include_user_mods create mode 100644 cime_config/testdefs/testmods_dirs/clm/compatible_finidat_f09/user_nl_clm diff --git a/cime_config/testdefs/testlist_clm.xml b/cime_config/testdefs/testlist_clm.xml index a6c3b3b5476..4630921b4cb 100644 --- a/cime_config/testdefs/testlist_clm.xml +++ b/cime_config/testdefs/testlist_clm.xml @@ -1065,7 +1065,7 @@ - + diff --git a/cime_config/testdefs/testmods_dirs/clm/compatible_finidat_f09/README b/cime_config/testdefs/testmods_dirs/clm/compatible_finidat_f09/README new file mode 100644 index 00000000000..566826ac8ea --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/compatible_finidat_f09/README @@ -0,0 +1,23 @@ +This testmod provides an finidat file compatible with the current model +configuration, at f09 resolution. + +This testmods directory is for running LII tests which compare the answers for a case where +initial condition interpolation is on (use_init_interp=T) to a case with it's off and ensures +they are exact. For the interpolated result to match uninterpolation, it needs to be a case that +essentially needs no interpolation so it's at the same resolution as the initial condition +file (finidat file). When surface datasets are changed, or the land-mask is changed, or an +important change is made to model physics (for example where new fields are added to the restart +file) -- you'll need to update the initial conditions file in this test (finidat file in +the user_nl_clm file). + +To update the initial conditions (finidat) file for this test: + +(1) Run the LII test; the 'base' case should run to completion even if the +no_interp test fails. + +(2) Copy the finidat_interp_dest.nc file from the 'base' case to the inputdata +space. Rename this to be similar to the name of the file pointed to in this +test's user_nl_clm file, but with a new creation date. + +(3) Update this the user_nl_clm file in this directory to point to the new +finidat file. diff --git a/cime_config/testdefs/testmods_dirs/clm/compatible_finidat_f09/include_user_mods b/cime_config/testdefs/testmods_dirs/clm/compatible_finidat_f09/include_user_mods new file mode 100644 index 00000000000..fe0e18cf889 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/compatible_finidat_f09/include_user_mods @@ -0,0 +1 @@ +../default diff --git a/cime_config/testdefs/testmods_dirs/clm/compatible_finidat_f09/user_nl_clm b/cime_config/testdefs/testmods_dirs/clm/compatible_finidat_f09/user_nl_clm new file mode 100644 index 00000000000..dc3dedf6c49 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/clm/compatible_finidat_f09/user_nl_clm @@ -0,0 +1 @@ +finidat = '$DIN_LOC_ROOT/lnd/clm2/initdata_map/clmi.I1850Clm50BgcCrop.1366-01-01.0.9x1.25_gx1v6_simyr1850_c180424.nc' From ba35b252915f83e8f1660a248fc5228279c7caaa Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Wed, 9 May 2018 16:03:58 -0600 Subject: [PATCH 18/19] Broadcast use_init_interp Without this change, processors other than masterproc don't have the right value of use_init_interp and so do the wrong thing with respect to calling initInterp_readnl. --- src/main/controlMod.F90 | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/main/controlMod.F90 b/src/main/controlMod.F90 index ed935937825..3cbae8ed680 100644 --- a/src/main/controlMod.F90 +++ b/src/main/controlMod.F90 @@ -15,7 +15,8 @@ module controlMod use shr_const_mod , only: SHR_CONST_CDAY use shr_log_mod , only: errMsg => shr_log_errMsg use abortutils , only: endrun - use spmdMod , only: masterproc + use spmdMod , only: masterproc, mpicom + use spmdMod , only: MPI_CHARACTER, MPI_INTEGER, MPI_LOGICAL, MPI_REAL8 use decompMod , only: clump_pproc use clm_varcon , only: h2osno_max, int_snow_max, n_melt_glcmec use clm_varpar , only: maxpatch_pft, maxpatch_glcmec, numrad, nlevsno @@ -438,6 +439,7 @@ subroutine control_init( ) ! Read in other namelists for other modules ! ---------------------------------------------------------------------- + call mpi_bcast (use_init_interp, 1, MPI_LOGICAL, 0, mpicom, ierr) if (use_init_interp) then call initInterp_readnl( NLFilename ) end if @@ -550,9 +552,6 @@ subroutine control_spmd() ! it to all processors, or collects data from ! all processors and writes it to disk. ! - ! !USES: - use spmdMod, only : mpicom, MPI_CHARACTER, MPI_INTEGER, MPI_LOGICAL, MPI_REAL8 - ! ! !ARGUMENTS: ! ! !LOCAL VARIABLES: From 191cca920c16e1c256a85fb531c9acde00e2c5bd Mon Sep 17 00:00:00 2001 From: Bill Sacks Date: Thu, 10 May 2018 13:48:48 -0600 Subject: [PATCH 19/19] Update ChangeLog --- doc/ChangeLog | 153 ++++++++++++++++++++++++++++++++++++++++++++++++++ doc/ChangeSum | 1 + 2 files changed, 154 insertions(+) diff --git a/doc/ChangeLog b/doc/ChangeLog index b7ff6410be0..c966d8142b2 100644 --- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -1,4 +1,157 @@ =============================================================== +Tag name: clm5.0.dev009 +Originator(s): sacks (Bill Sacks) +Date: Thu May 10 13:48:32 MDT 2018 +One-line Summary: New init_interp method + +Purpose of changes +------------------ + +The main change in this tag is to introduce a new method for +init_interp. This is needed in order for subgrid areas to be the same in +initialization in runs with and without init_interp. + +The choice of init_interp method is controlled by a new namelist +variable, init_interp_method, whose documentation appears below: + + Method to use for init_interp. Only applies when use_init_interp = .true. + + 'general': The general-purpose method that can be used when changing + grids, configurations, etc. This starts off with subgrid areas taken + from the surface dataset. + + 'use_finidat_areas': This starts off with subgrid areas taken from the + input finidat file. This is needed to achieve bit-for-bit results in a + coupled case (where areas in initialization impact initial fields sent + to the atmosphere) (but using the 'general' method will typically have + only a very minor impact on results in this case). For this method to + work, the input finidat file needs to be at the same resolution as the + current configuration. So this is a less general form of + init_interp. However, it can be used when transitioning from a spinup + run to a transient run, or in other cases where the only difference is + in internal memory allocation. In order to catch possible problems, this + uses a different algorithm for finding the input point for each output + point, which ensures that each active output point is associated with + exactly one input point with the same latitude, longitude and type. This + method requires (a) the same grid for input and output, within roundoff; + (b) any non-zero-weight point in the input must have memory allocated + for it in this grid cell in the output (this will be satisfied if the + point is non-zero-weight on the surface dataset or if it's a point for + which we allocate memory even for zero-weight points); (c) any active + point in the output (based on the surface dataset and rules for + determining active points) must have a matching point in this grid cell + in the input. + +This tag also has some other changes: + +- Fixes #347 - endrun message behavior: needed for some of the new unit + tests + +- Fixes #345 - Add a cmip6_evolving_icesheet usermods directory: + unrelated, but folded in here for convenience of testing) + + +Bugs fixed or introduced +------------------------ + +Issues fixed (include CTSM Issue #): +- #346: New mode of operation for init_interp: Copy subgrid areas, too +- #347: endrun message behavior +- #345: Add a cmip6_evolving_icesheet usermods directory + + +Notes of particular relevance for users +--------------------------------------- + +Caveats for users (e.g., need to interpolate initial conditions): none + +Changes to CLM's user interface (e.g., new/renamed XML or namelist variables): +- New namelist variable - init_interp_method - documented above. Default + behavior is the same as before ('general') + +Changes made to namelist defaults (e.g., changed parameter values): none + +Changes to the datasets (e.g., parameter, surface or initial files): none + +Substantial timing or memory changes: none + +Notes of particular relevance for developers: (including Code reviews and testing) +--------------------------------------------- + +Caveats for developers (e.g., code that is duplicated that requires double maintenance): none + +Changes to tests or testing: +- Added a new test type, similar to LII, that covers the new init_interp behavior +- Removed some no-longer-needed tests with 'interp' testmods + +Code reviewed by: self + +Did you follow the steps in .CLMTrunkChecklist: yes + +CLM testing: + + [PASS means all tests PASS and OK means tests PASS other than expected fails.] + + build-namelist tests: + + cheyenne - ok + + tests pass, namelists differ as expected + + unit-tests (components/clm/src): + + cheyenne - pass + + tools-tests (components/clm/test/tools): + + cheyenne - not run + + PTCLM testing (components/clm/tools/shared/PTCLM/test): + + cheyenne - not run + + regular tests (aux_clm): + + cheyenne_intel ---- pass + cheyenne_gnu ------ pass + hobart_nag -------- pass + hobart_pgi -------- pass + hobart_intel ------ pass + + Unexpected baseline failures for these two tests that are in the + expected fail list for other reasons: + FAIL ERP_Ld5.f10_f10_musgs.I2000Clm50Vic.cheyenne_gnu.clm-decStart COMPARE_base_rest + FAIL ERP_D.f10_f10_musgs.IHistClm50Bgc.cheyenne_gnu.clm-decStart COMPARE_base_rest + + Differences are just in some cism fields that I’m not concerned + about: roundoff-level differences in internal_time; big differences + in tempstag, uvel and vvel, but those fields don’t really make sense + in this configuration (and I have removed them in the latest CISM + tag). (Note that I did two runs of each of these tests from my + branch, and for both tests, the two runs had identical cism hist + files, so I don’t think this is a reproducibility problem.) + + +CLM tag used for the baseline comparisons: clm5.0.dev008 + + +Answer changes +-------------- + +Changes answers relative to baseline: NO + + +Detailed list of changes +------------------------ + +List any externals directories updated (cime, rtm, mosart, cism, fates, etc.): none + +Pull Requests that document the changes (include PR ids): +(https://github.com/ESCOMP/ctsm/pull) +- #352 + +=============================================================== +=============================================================== Tag name: clm5.0.dev008 Originator(s): erik (Erik Kluzek,UCAR/TSS,303-497-1326) Date: Fri Apr 27 13:28:41 MDT 2018 diff --git a/doc/ChangeSum b/doc/ChangeSum index b54f022d196..c6409e49128 100644 --- a/doc/ChangeSum +++ b/doc/ChangeSum @@ -1,5 +1,6 @@ Tag Who Date Summary ============================================================================================================================ + clm5.0.dev009 sacks 05/10/2018 New init_interp method clm5.0.dev008 erik 04/27/2018 With FUN subtract out soil nitrification flux of plant uptake of soil NH3 and NO3 clm5.0.dev007 erik 04/24/2018 Bring in a few answer changing things: FATES, cism updates, IC file fix, testing 1850 compset use 1850 orbit