Merge branch 'hillslope_hydrology' into hillslope_hydrology-ssr

src/biogeophys/HillslopeHydrologyMod.F90

@@ -27,8 +27,6 @@ module HillslopeHydrologyMod
   public SetHillslopeSoilThickness
   public HillslopeSoilThicknessProfile
   public HillslopeSetLowlandUplandPfts
-  public HillslopeDominantPftIndex
-  public HillslopeTwoLargestPftIndices
   public HillslopeDominantPft
   public HillslopeDominantLowlandPft
   public HillslopePftFromFile
@@ -184,6 +182,7 @@ subroutine InitHillslope(bounds,fsurdat)
     use fileutils       , only : getfil
     use clm_varcon      , only : spval, ispval, grlnd 
     use landunit_varcon , only : istsoil
+    use subgridWeightsMod , only : compute_higher_order_weights
     use ncdio_pio
 
     !
@@ -577,9 +576,6 @@ subroutine InitHillslope(bounds,fsurdat)
     ! to ensure patch exists in every gridcell
     if (pft_distribution_method == pft_from_file) then
        call HillslopePftFromFile(bounds,col_pftndx)
-    else if (pft_distribution_method == pft_uniform_dominant_pft) then
-       ! Specify single dominant pft per gridcell
-       call HillslopeDominantPft(bounds)
     else if (pft_distribution_method == pft_lowland_dominant_pft) then
        ! Specify different pfts for uplands / lowlands
        call HillslopeDominantLowlandPft(bounds)
@@ -595,6 +591,9 @@ subroutine InitHillslope(bounds,fsurdat)
        deallocate(col_pftndx)
     endif
 
+    ! Update higher order weights and check that weights sum to 1
+    call compute_higher_order_weights(bounds)
+
     call ncd_pio_closefile(ncid)
 
   end subroutine InitHillslope
@@ -865,76 +864,6 @@ subroutine HillslopeSetLowlandUplandPfts(bounds,lowland_ivt,upland_ivt)
 
   end subroutine HillslopeSetLowlandUplandPfts
 
-  !------------------------------------------------------------------------
-  subroutine HillslopeDominantPftIndex(weights,lbound,pindex)
-    !
-    ! !DESCRIPTION:
-    ! Locate each gridcell's most dominant pft on the input dataset.
-    ! This is different than using n_dom_pts = 1, because it is meant
-    ! to be applied only to gridcells having active hillslope hydrology.
-    ! This routine is called from surfrd_hillslope.
-
-    !
-    ! !USES
-    !
-    ! !ARGUMENTS:
-    real(r8), intent(in)  :: weights(:)        ! array of patch weights
-    integer,  intent(in)  :: lbound            ! lower bound of weights
-    integer,  intent(out) :: pindex            ! index of largest weight
-    !
-    ! !LOCAL VARIABLES:
-    integer :: pdom(1)
-
-    !------------------------------------------------------------------------
-
-    pdom = maxloc(weights)
-    pindex = pdom(1) + (lbound - 1)
-
-  end subroutine HillslopeDominantPftIndex
-
-  !------------------------------------------------------------------------
-  subroutine HillslopeTwoLargestPftIndices(weights,lbound,pindex1,pindex2)
-    !
-    ! !DESCRIPTION:
-    ! Locate each gridcell's two most dominant patches on the input dataset.
-    ! This is different than using n_dom_pts = 2, because it is meant
-    ! to be applied only to gridcells having active hillslope hydrology.
-    ! This routine is called from surfrd_hillslope.
-
-    !
-    ! !USES
-    !
-    ! !ARGUMENTS:
-    real(r8), intent(in)  :: weights(:)        ! array of patch weights
-    integer,  intent(in)  :: lbound            ! lower bound of weights
-    integer,  intent(out) :: pindex1           ! index of largest weight
-    integer,  intent(out) :: pindex2           ! index of next largest weight
-    !
-    ! !LOCAL VARIABLES:
-    integer :: pdom(1),psubdom(1)
-    real(r8),allocatable :: mask(:)
-
-    !------------------------------------------------------------------------
-
-    pdom = maxloc(weights)
-    ! create mask to exclude pdom
-    allocate(mask(size(weights)))
-    mask(:) = 1.
-    mask(pdom(1)) = 0.
-    ! check that more than one nonzero patch weight exists,
-    ! if not return identical indices
-    if (sum(mask*weights) > 0) then
-       psubdom = maxloc(mask*weights)
-    else
-       psubdom = pdom
-    endif
-    deallocate(mask)
-
-    pindex1 = pdom(1) + (lbound - 1)
-    pindex2 = psubdom(1) + (lbound - 1)
-
-  end subroutine HillslopeTwoLargestPftIndices
-
   !------------------------------------------------------------------------
   subroutine HillslopeDominantPft(bounds)
     !
@@ -952,6 +881,7 @@ subroutine HillslopeDominantPft(bounds)
     use decompMod       , only : get_clump_bounds, get_proc_clumps
     use clm_varcon      , only : ispval
     use PatchType       , only : patch
+    use array_utils     , only : find_k_max_indices
     !
     ! !ARGUMENTS:
     type(bounds_type), intent(in) :: bounds
@@ -965,7 +895,8 @@ subroutine HillslopeDominantPft(bounds)
 
     do c = bounds%begc,bounds%endc
        if (col%is_hillslope_column(c) .and. (col%npatches(c) > 1)) then
-          pdom = maxloc(patch%wtcol(col%patchi(c):col%patchf(c)))
+
+          call find_k_max_indices(patch%wtcol(col%patchi(c):col%patchf(c)),1,1,pdom)
           pdom = pdom + (col%patchi(c) - 1)
 
           sum_wtcol = sum(patch%wtcol(col%patchi(c):col%patchf(c)))
@@ -1003,32 +934,34 @@ subroutine HillslopeDominantLowlandPft(bounds)
     use clm_varcon      , only : ispval
     use PatchType       , only : patch
     use pftconMod       , only : pftcon, ndllf_evr_tmp_tree, nc3_nonarctic_grass, nc4_grass
+    use array_utils     , only : find_k_max_indices
     !
     ! !ARGUMENTS:
     type(bounds_type), intent(in) :: bounds
     !
     ! !LOCAL VARIABLES:
-    integer :: nc,p,pu,pl,l,c    ! indices
-    integer :: pdom(1),psubdom(1)
+    integer :: p,c    ! indices
     integer :: plow, phigh
+    integer :: max_index(1)
+    integer, allocatable  :: max_indices(:)    ! largest weight pft indices
     real(r8) :: sum_wtcol, sum_wtlun, sum_wtgrc
-    real(r8),allocatable :: mask(:)
 
     !------------------------------------------------------------------------
 
+    allocate(max_indices(2))
     do c = bounds%begc,bounds%endc
        if (col%is_hillslope_column(c)) then
-          pdom = maxloc(patch%wtcol(col%patchi(c):col%patchf(c)))
-          ! create mask to exclude pdom
-          allocate(mask(col%npatches(c)))
-          mask(:) = 1.
-          mask(pdom(1)) = 0.
 
-          pdom = pdom + (col%patchi(c) - 1)
-
-          psubdom = maxloc(mask*patch%wtcol(col%patchi(c):col%patchf(c)))
-          psubdom = psubdom + (col%patchi(c) - 1)
-          deallocate(mask)
+          ! if only one pft exists, find dominant pft index and set 2nd index to the same value
+
+          if (size(patch%wtcol(col%patchi(c):col%patchf(c))) == 1) then
+             call find_k_max_indices(patch%wtcol(col%patchi(c):col%patchf(c)),1,1,max_index)
+             max_indices(1) = max_index(1) + (col%patchi(c) - 1)
+             max_indices(2) = max_indices(1)
+          else
+             call find_k_max_indices(patch%wtcol(col%patchi(c):col%patchf(c)),1,2,max_indices)
+             max_indices = max_indices + (col%patchi(c) - 1)
+          endif
 
           sum_wtcol = sum(patch%wtcol(col%patchi(c):col%patchf(c)))
           sum_wtlun = sum(patch%wtlunit(col%patchi(c):col%patchf(c)))
@@ -1037,34 +970,36 @@ subroutine HillslopeDominantLowlandPft(bounds)
           patch%wtcol(col%patchi(c):col%patchf(c)) = 0._r8
           patch%wtlunit(col%patchi(c):col%patchf(c)) = 0._r8
           patch%wtgcell(col%patchi(c):col%patchf(c)) = 0._r8
-
-          ! assumes trees are 1-8, shrubs 9-11, and grasses 12-14
-          ! and puts the lowest ivt on the lowland column
-          if ((patch%itype(pdom(1)) > patch%itype(psubdom(1)) &
-               .and. patch%itype(psubdom(1)) > 0) .or. &
-               (patch%itype(pdom(1)) == 0)) then
-             plow = psubdom(1)
-             phigh = pdom(1)
+
+          ! Put the highest stature vegetation on the lowland column
+          ! non-tree and tree       ; place tree on lowland
+          ! grass and shrub         ; place shrub on lowland
+          ! bare soil and vegetation; place vegetation on lowland
+          if ((.not. pftcon%is_tree(patch%itype(max_indices(1))) .and. pftcon%is_tree(patch%itype(max_indices(2)))) &
+               .or. (pftcon%is_grass(patch%itype(max_indices(1))) .and. pftcon%is_shrub(patch%itype(max_indices(2)))) &
+               .or. (patch%itype(max_indices(1)) == 0)) then
+             plow = max_indices(2)
+             phigh = max_indices(1)
           else
-             plow = pdom(1)
-             phigh = psubdom(1)
+             plow = max_indices(1)
+             phigh = max_indices(2)
           endif
 
           ! Special cases (subjective)
 
           ! if NET/BDT assign BDT to lowland
-          if ((patch%itype(pdom(1)) == ndllf_evr_tmp_tree) .and. pftcon%is_tree(patch%itype(psubdom(1)))) then
-             plow = psubdom(1)        
-             phigh = pdom(1)
+          if ((patch%itype(max_indices(1)) == ndllf_evr_tmp_tree) .and. pftcon%is_tree(patch%itype(max_indices(2)))) then
+             plow = max_indices(2)        
+             phigh = max_indices(1)
           endif
           ! if C3/C4 assign C4 to lowland
-          if ((patch%itype(pdom(1)) == nc4_grass) .and. (patch%itype(psubdom(1)) == nc3_nonarctic_grass)) then
-             plow = pdom(1)        
-             phigh = psubdom(1)
+          if ((patch%itype(max_indices(1)) == nc4_grass) .and. (patch%itype(max_indices(2)) == nc3_nonarctic_grass)) then
+             plow = max_indices(1)        
+             phigh = max_indices(2)
           endif
-          if ((patch%itype(pdom(1)) == nc3_nonarctic_grass) .and. (patch%itype(psubdom(1)) == nc4_grass)) then
-             plow = psubdom(1)        
-             phigh = pdom(1)
+          if ((patch%itype(max_indices(1)) == nc3_nonarctic_grass) .and. (patch%itype(max_indices(2)) == nc4_grass)) then
+             plow = max_indices(2)        
+             phigh = max_indices(1)
           endif
 
           if(col%cold(c) == ispval) then
@@ -1080,7 +1015,8 @@ subroutine HillslopeDominantLowlandPft(bounds)
           endif
        endif
     enddo    ! end loop c
-
+    deallocate(max_indices)
+
   end subroutine HillslopeDominantLowlandPft
 
   !------------------------------------------------------------------------
@@ -1095,6 +1031,7 @@ subroutine HillslopePftFromFile(bounds,col_pftndx)
     ! !USES
     use ColumnType      , only : col                
     use PatchType       , only : patch
+    use clm_varpar      , only : natpft_lb
 
     !
     ! !ARGUMENTS:
@@ -1115,6 +1052,8 @@ subroutine HillslopePftFromFile(bounds,col_pftndx)
           npatches_per_column = 0
           do p = col%patchi(c), col%patchf(c)
              patch%itype(p) = col_pftndx(c)
+             ! update mxy as is done in initSubgridMod.add_patch
+             patch%mxy(p) = patch%itype(p) + (1 - natpft_lb)
              npatches_per_column = npatches_per_column + 1
           enddo
           if (check_npatches) then

src/cpl/utils/lnd_import_export_utils.F90

@@ -163,7 +163,7 @@ subroutine check_for_nans(array, fname, begg, direction)
        write(iulog,*) 'Which are NaNs = ', isnan(array)
        do i = 1, size(array)
           if (isnan(array(i))) then
-             write(iulog,*) "NaN found in field ", trim(fname), ' at gridcell index ',begg+i-1,grc%londeg(begg+i-1),grc%latdeg(begg+i-1)
+             write(iulog,*) "NaN found in field ", trim(fname), ' at gridcell index/lon/lat: ',begg+i-1,grc%londeg(begg+i-1),grc%latdeg(begg+i-1)
           end if
        end do
        call shr_sys_abort(' ERROR: One or more of the CTSM cap '//direction//' fields are NaN ' )

src/main/atm2lndMod.F90

@@ -809,8 +809,7 @@ subroutine downscale_hillslope_solar(bounds, atm2lnd_inst, surfalb_inst)
             g = col%gridcell(c)
             do n = 1,numrad
                ! absorbed energy is solar flux x area landunit (sum_wtgcell)
-               !if(sum_solar(g,n) > 0._r8) then
-               if(sum_solar(g,n) > 0._r8.and.illum_frac(g) > illumination_threshold) then
+               if(sum_solar(g,n) > 0._r8 .and. illum_frac(g) > illumination_threshold) then
                   norm(n) = sum_wtgcell(g)*forc_solad_grc(g,n)/sum_solar(g,n)
                   forc_solad_col(c,n)  = forc_solad_col(c,n)*norm(n)
                else

src/main/initVerticalMod.F90

@@ -182,22 +182,6 @@ subroutine initVertical(bounds, glc_behavior, thick_wall, thick_roof)
     integer               :: begc, endc
     integer               :: begl, endl
     integer               :: jmin_bedrock
-    ! Possible values for levgrnd_class. The important thing is that, for a given column,
-    ! layers that are fundamentally different (e.g., soil vs bedrock) have different
-    ! values. This information is used in the vertical interpolation in init_interp.
-    !
-    ! IMPORTANT: These values should not be changed lightly. e.g., try to avoid changing
-    ! the values assigned to LEVGRND_CLASS_STANDARD, LEVGRND_CLASS_DEEP_BEDROCK, etc.  The
-    ! problem with changing these is that init_interp expects that layers with a value of
-    ! (e.g.) 1 on the source file correspond to layers with a value of 1 on the
-    ! destination file. So if you change the values of these constants, you either need to
-    ! adequately inform users of this change, or build in some translation mechanism in
-    ! init_interp (such as via adding more metadata to the restart file on the meaning of
-    ! these different values).
-    !
-    ! The distinction between "shallow" and "deep" bedrock is not made explicitly
-    ! elsewhere. But, since these classes have somewhat different behavior, they are
-    ! distinguished explicitly here.
     character(len=*), parameter :: subname = 'initVertical'
     !------------------------------------------------------------------------
 

src/main/surfrdMod.F90

@@ -899,18 +899,19 @@ subroutine surfrd_hillslope(begg, endg, ncid, ns)
     use ncdio_pio,  only : ncd_inqdid, ncd_inqdlen
     use pftconMod , only : noveg
     use HillslopeHydrologyMod, only : pft_distribution_method, pft_from_file, pft_uniform_dominant_pft, pft_lowland_dominant_pft, pft_lowland_upland
-    use HillslopeHydrologyMod, only : HillslopeDominantPftIndex,HillslopeTwoLargestPftIndices
-
+    use array_utils, only: find_k_max_indices
+    
     !
     ! !ARGUMENTS:
     integer, intent(in) :: begg, endg
     type(file_desc_t),intent(inout) :: ncid   ! netcdf id
     integer          ,intent(in)    :: ns     ! domain size
     !
     ! !LOCAL VARIABLES:
-    integer  :: g, nh, m, n                    ! index
+    integer  :: g, nh, m, n                       ! index
     integer  :: dimid,varid                    ! netCDF id's
     integer  :: ier                            ! error status	
+    integer, allocatable  :: max_indices(:)    ! largest weight pft indices
     logical  :: readvar                        ! is variable on dataset
     integer,pointer :: arrayl(:)               ! local array (needed because ncd_io expects a pointer)
     character(len=32) :: subname = 'surfrd_hillslope'  ! subroutine name
@@ -962,36 +963,50 @@ subroutine surfrd_hillslope(begg, endg, ncid, ns)
     ! pft_uniform_dominant_pft uses the patch with the
     ! largest weight for all hillslope columns in the gridcell
     if (pft_distribution_method == pft_uniform_dominant_pft) then
+       allocate(max_indices(1))
        do g = begg, endg
           ! If hillslopes will be used in a gridcell, modify wt_nat_patch,
           ! otherwise use original patch distribution
           if(ncolumns_hillslope(g) > 0) then
-             call HillslopeDominantPftIndex(wt_nat_patch(g,:),natpft_lb,m)
+
+             call find_k_max_indices(wt_nat_patch(g,:),natpft_lb,1,max_indices)
              wt_nat_patch(g,:) = 0._r8
-             wt_nat_patch(g,m) = 100._r8
+             wt_nat_patch(g,max_indices(1)) = 100._r8
+
           endif
        enddo
+       deallocate(max_indices)
     endif
 
     ! pft_lowland_dominant_pft uses the two patches with the
     ! largest weights for the hillslope columns in the gridcell
     if (pft_distribution_method == pft_lowland_dominant_pft) then
+       allocate(max_indices(2))
        do g = begg, endg
           ! If hillslopes will be used in a gridcell, modify wt_nat_patch, otherwise use original patch distribution
           if(ncolumns_hillslope(g) > 0) then
-             call HillslopeTwoLargestPftIndices(wt_nat_patch(g,:),natpft_lb,m,n)
+
              ! Preserve the relative weights of the largest and
              ! next largest weights using arbitrarily chosen values
-             ! (i.e. m should be larger than n)  This will minimize
-             ! memory usage while still allowing HillslopeDominantLowlandPft
-             ! to pick out the two largest patch types.
-             if(m /= n) then
+             ! (i.e. 1 should be larger than 2) that sum to 100.
+             ! This will minimize memory usage while still allowing
+             ! HillslopeDominantLowlandPft to pick out the two largest patch types.
+
+             call find_k_max_indices(wt_nat_patch(g,:),natpft_lb,2,max_indices)
+             ! check that 2nd index weight is non-zero
+             if (wt_nat_patch(g,max_indices(2)) > 0._r8) then
+                wt_nat_patch(g,:) = 0._r8
+                wt_nat_patch(g,max_indices(1)) = 75._r8
+                wt_nat_patch(g,max_indices(2)) = 25._r8
+             else
+                ! if only one pft exists, set its weight to 100 per cent
                 wt_nat_patch(g,:) = 0._r8
-                wt_nat_patch(g,m) = 75._r8
-                wt_nat_patch(g,n) = 25._r8
+                wt_nat_patch(g,max_indices(1)) = 100._r8
              endif
+
           endif
        enddo
+       deallocate(max_indices)
     endif
 
   end subroutine surfrd_hillslope