diff --git a/src/biogeophys/HillslopeHydrologyMod.F90 b/src/biogeophys/HillslopeHydrologyMod.F90 index ce015e25ae..55fa5c53bd 100644 --- a/src/biogeophys/HillslopeHydrologyMod.F90 +++ b/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 diff --git a/src/cpl/utils/lnd_import_export_utils.F90 b/src/cpl/utils/lnd_import_export_utils.F90 index 788afc578a..1b40cb0e6c 100644 --- a/src/cpl/utils/lnd_import_export_utils.F90 +++ b/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 ' ) diff --git a/src/main/atm2lndMod.F90 b/src/main/atm2lndMod.F90 index c3a4766ddb..e1904125af 100644 --- a/src/main/atm2lndMod.F90 +++ b/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 diff --git a/src/main/initVerticalMod.F90 b/src/main/initVerticalMod.F90 index b492708a14..92788a4609 100644 --- a/src/main/initVerticalMod.F90 +++ b/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' !------------------------------------------------------------------------ diff --git a/src/main/surfrdMod.F90 b/src/main/surfrdMod.F90 index 6fad760f01..3935b170d6 100644 --- a/src/main/surfrdMod.F90 +++ b/src/main/surfrdMod.F90 @@ -899,8 +899,8 @@ 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 @@ -908,9 +908,10 @@ subroutine surfrd_hillslope(begg, endg, ncid, ns) 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