diff --git a/src/biogeophys/CanopyTemperatureMod.F90 b/src/biogeophys/CanopyTemperatureMod.F90 index b1c3c6cc67..fbe76b561e 100644 --- a/src/biogeophys/CanopyTemperatureMod.F90 +++ b/src/biogeophys/CanopyTemperatureMod.F90 @@ -47,6 +47,9 @@ module CanopyTemperatureMod !------------------------------------------------------------------------------ subroutine CanopyTemperature(bounds, & num_nolakec, filter_nolakec, num_nolakep, filter_nolakep, & +!KO + num_urbanc, filter_urbanc, & +!KO clm_fates, & atm2lnd_inst, canopystate_inst, soilstate_inst, frictionvel_inst, & waterstatebulk_inst, waterdiagnosticbulk_inst, waterfluxbulk_inst, & @@ -90,6 +93,10 @@ subroutine CanopyTemperature(bounds, & integer , intent(in) :: filter_nolakec(:) ! column filter for non-lake points integer , intent(in) :: num_nolakep ! number of column non-lake points in patch filter integer , intent(in) :: filter_nolakep(:) ! patch filter for non-lake points +!KO + integer , intent(in) :: num_urbanc ! number of urban columns in clump + integer , intent(in) :: filter_urbanc(:) ! urban column filter +!KO type(hlm_fates_interface_type), intent(inout) :: clm_fates type(atm2lnd_type) , intent(in) :: atm2lnd_inst type(canopystate_type) , intent(inout) :: canopystate_inst @@ -107,6 +114,9 @@ subroutine CanopyTemperature(bounds, & integer :: j ! soil/snow level index integer :: fp ! lake filter patch index integer :: fc ! lake filter column index +!KO + integer :: nlev ! greater of nlevgrnd and nlevurb +!KO real(r8) :: qred ! soil surface relative humidity real(r8) :: avmuir ! ir inverse optical depth per unit leaf area real(r8) :: eg ! water vapor pressure at temperature T [pa] @@ -218,10 +228,14 @@ subroutine CanopyTemperature(bounds, & do j = -nlevsno+1, nlevgrnd do fc = 1,num_nolakec c = filter_nolakec(fc) - if ((col%itype(c) == icol_sunwall .or. col%itype(c) == icol_shadewall & - .or. col%itype(c) == icol_roof) .and. j > nlevurb) then - tssbef(c,j) = spval - else +!KO if ((col%itype(c) == icol_sunwall .or. col%itype(c) == icol_shadewall & +!KO .or. col%itype(c) == icol_roof) .and. j > nlevurb) then +!KO tssbef(c,j) = spval +!KO else +!KO + if (col%itype(c) /= icol_sunwall .and. col%itype(c) /= icol_shadewall & + .and. col%itype(c) /= icol_roof) then +!KO tssbef(c,j) = t_soisno(c,j) end if ! record t_h2osfc prior to updating @@ -229,6 +243,25 @@ subroutine CanopyTemperature(bounds, & end do end do +!KO + nlev = max0(nlevgrnd,nlevurb) + do j = -nlevsno+1, nlev + do fc = 1,num_urbanc + c = filter_urbanc(fc) + if (col%itype(c) == icol_sunwall .or. col%itype(c) == icol_shadewall & + .or. col%itype(c) == icol_roof) then + if (j > nlevurb) then + tssbef(c,j) = spval + else + tssbef(c,j) = t_soisno(c,j) + end if + end if + ! record t_h2osfc prior to updating + t_h2osfc_bef(c) = t_h2osfc(c) + end do + end do +!KO + ! calculate moisture stress/resistance for soil evaporation call calc_soilevap_resis(bounds, num_nolakec, filter_nolakec, & soilstate_inst, waterstatebulk_inst, waterdiagnosticbulk_inst, temperature_inst) diff --git a/src/biogeophys/HydrologyDrainageMod.F90 b/src/biogeophys/HydrologyDrainageMod.F90 index c9f2c9eba1..f4cc9175d7 100644 --- a/src/biogeophys/HydrologyDrainageMod.F90 +++ b/src/biogeophys/HydrologyDrainageMod.F90 @@ -151,14 +151,30 @@ subroutine HydrologyDrainage(bounds, & do j = 1, nlevgrnd do fc = 1, num_nolakec c = filter_nolakec(fc) - if ((ctype(c) == icol_sunwall .or. ctype(c) == icol_shadewall & - .or. ctype(c) == icol_roof) .and. j > nlevurb) then - else +!KO if ((ctype(c) == icol_sunwall .or. ctype(c) == icol_shadewall & +!KO .or. ctype(c) == icol_roof) .and. j > nlevurb) then +!KO + if (ctype(c) /= icol_sunwall .and. ctype(c) /= icol_shadewall & + .and. ctype(c) /= icol_roof) then +!KO +!KO else h2osoi_vol(c,j) = h2osoi_liq(c,j)/(dz(c,j)*denh2o) + h2osoi_ice(c,j)/(dz(c,j)*denice) end if end do end do +!KO + do j = 1, nlevurb + do fc = 1, num_urbanc + c = filter_urbanc(fc) + if (col%itype(c) == icol_sunwall .or. col%itype(c) == icol_shadewall & + .or. col%itype(c) == icol_roof) then + h2osoi_vol(c,j) = h2osoi_liq(c,j)/(dz(c,j)*denh2o) + h2osoi_ice(c,j)/(dz(c,j)*denice) + end if + end do + end do +!KO + call ComputeWaterMassNonLake(bounds, num_nolakec, filter_nolakec, & waterstatebulk_inst, waterdiagnosticbulk_inst, & subtract_dynbal_baselines = .false., & diff --git a/src/biogeophys/HydrologyNoDrainageMod.F90 b/src/biogeophys/HydrologyNoDrainageMod.F90 index 1823c826c0..a942c0eff6 100644 --- a/src/biogeophys/HydrologyNoDrainageMod.F90 +++ b/src/biogeophys/HydrologyNoDrainageMod.F90 @@ -514,13 +514,29 @@ subroutine HydrologyNoDrainage(bounds, & do j = 1, nlevgrnd do fc = 1, num_nolakec c = filter_nolakec(fc) - if ((ctype(c) == icol_sunwall .or. ctype(c) == icol_shadewall & - .or. ctype(c) == icol_roof) .and. j > nlevurb) then - else +!KO if ((ctype(c) == icol_sunwall .or. ctype(c) == icol_shadewall & +!KO .or. ctype(c) == icol_roof) .and. j > nlevurb) then +!KO + if (ctype(c) /= icol_sunwall .and. ctype(c) /= icol_shadewall & + .and. ctype(c) /= icol_roof) then +!KO +!KO else + h2osoi_vol(c,j) = h2osoi_liq(c,j)/(dz(c,j)*denh2o) + h2osoi_ice(c,j)/(dz(c,j)*denice) + end if + end do + end do + +!KO + do j = 1, nlevurb + do fc = 1, num_urbanc + c = filter_urbanc(fc) + if (col%itype(c) == icol_sunwall .or. col%itype(c) == icol_shadewall & + .or. col%itype(c) == icol_roof) then h2osoi_vol(c,j) = h2osoi_liq(c,j)/(dz(c,j)*denh2o) + h2osoi_ice(c,j)/(dz(c,j)*denice) end if end do end do +!KO ! if (use_cn) then ! Update soilpsi. diff --git a/src/biogeophys/SoilWaterMovementMod.F90 b/src/biogeophys/SoilWaterMovementMod.F90 index 05be0685f4..ffa37d55fa 100644 --- a/src/biogeophys/SoilWaterMovementMod.F90 +++ b/src/biogeophys/SoilWaterMovementMod.F90 @@ -2188,43 +2188,59 @@ subroutine TridiagonalCol (ci, lbj, ubj, jtop, a, b, c, r, u) bet = b(jtop) do j = lbj, ubj - if ((col%itype(ci) == icol_sunwall .or. col%itype(ci) == icol_shadewall & - .or. col%itype(ci) == icol_roof) .and. j <= nlevurb) then - if (j >= jtop) then - if (j == jtop) then - u(j) = r(j) / bet - else - gam(j) = c(j-1) / bet - bet = b(j) - a(j) * gam(j) - u(j) = (r(j) - a(j)*u(j-1)) / bet - end if - end if - else if (col%itype(ci) /= icol_sunwall .and. col%itype(ci) /= icol_shadewall & - .and. col%itype(ci) /= icol_roof) then - if (j >= jtop) then - if (j == jtop) then - u(j) = r(j) / bet - else - gam(j) = c(j-1) / bet - bet = b(j) - a(j) * gam(j) - u(j) = (r(j) - a(j)*u(j-1)) / bet - end if +!KO if ((col%itype(ci) == icol_sunwall .or. col%itype(ci) == icol_shadewall & +!KO .or. col%itype(ci) == icol_roof) .and. j <= nlevurb) then +!KO if (j >= jtop) then +!KO if (j == jtop) then +!KO u(j) = r(j) / bet +!KO else +!KO gam(j) = c(j-1) / bet +!KO bet = b(j) - a(j) * gam(j) +!KO u(j) = (r(j) - a(j)*u(j-1)) / bet +!KO end if +!KO end if +!KO else if (col%itype(ci) /= icol_sunwall .and. col%itype(ci) /= icol_shadewall & +!KO .and. col%itype(ci) /= icol_roof) then +!KO if (j >= jtop) then +!KO if (j == jtop) then +!KO u(j) = r(j) / bet +!KO else +!KO gam(j) = c(j-1) / bet +!KO bet = b(j) - a(j) * gam(j) +!KO u(j) = (r(j) - a(j)*u(j-1)) / bet +!KO end if +!KO end if +!KO end if +!KO + if (j >= jtop) then + if (j == jtop) then + u(j) = r(j) / bet + else + gam(j) = c(j-1) / bet + bet = b(j) - a(j) * gam(j) + u(j) = (r(j) - a(j)*u(j-1)) / bet end if end if +!KO end do do j = ubj-1,lbj,-1 - if ((col%itype(ci) == icol_sunwall .or. col%itype(ci) == icol_shadewall & - .or. col%itype(ci) == icol_roof) .and. j <= nlevurb-1) then - if (j >= jtop) then - u(j) = u(j) - gam(j+1) * u(j+1) - end if - else if (col%itype(ci) /= icol_sunwall .and. col%itype(ci) /= icol_shadewall & - .and. col%itype(ci) /= icol_roof) then - if (j >= jtop) then - u(j) = u(j) - gam(j+1) * u(j+1) - end if +!KO if ((col%itype(ci) == icol_sunwall .or. col%itype(ci) == icol_shadewall & +!KO .or. col%itype(ci) == icol_roof) .and. j <= nlevurb-1) then +!KO if (j >= jtop) then +!KO u(j) = u(j) - gam(j+1) * u(j+1) +!KO end if +!KO else if (col%itype(ci) /= icol_sunwall .and. col%itype(ci) /= icol_shadewall & +!KO .and. col%itype(ci) /= icol_roof) then +!KO if (j >= jtop) then +!KO u(j) = u(j) - gam(j+1) * u(j+1) +!KO end if +!KO end if +!KO + if (j >= jtop) then + u(j) = u(j) - gam(j+1) * u(j+1) end if +!KO end do end subroutine TridiagonalCol diff --git a/src/biogeophys/TotalWaterAndHeatMod.F90 b/src/biogeophys/TotalWaterAndHeatMod.F90 index be06681fb0..566de39e1c 100644 --- a/src/biogeophys/TotalWaterAndHeatMod.F90 +++ b/src/biogeophys/TotalWaterAndHeatMod.F90 @@ -348,6 +348,9 @@ subroutine AccumulateSoilLiqIceMassNonLake(bounds, num_c, filter_c, & ! !LOCAL VARIABLES: integer :: c, j, fc ! indices logical :: has_h2o ! whether this point potentially has water to add +!KO + integer :: nlev ! greater of nlevgrnd and nlevurb +!KO character(len=*), parameter :: subname = 'AccumulateSoilLiqIceMassNonLake' !----------------------------------------------------------------------- @@ -360,7 +363,13 @@ subroutine AccumulateSoilLiqIceMassNonLake(bounds, num_c, filter_c, & h2osoi_liq => waterstate_inst%h2osoi_liq_col & ! Input: [real(r8) (:,:) ] liquid water (kg/m2) ) - do j = 1, nlevgrnd +!KO + nlev = max0(nlevgrnd,nlevurb) +!KO +!KO do j = 1, nlevgrnd +!KO + do j = 1, nlev +!KO do fc = 1, num_c c = filter_c(fc) if (col%itype(c) == icol_sunwall .or. col%itype(c) == icol_shadewall) then @@ -372,7 +381,12 @@ subroutine AccumulateSoilLiqIceMassNonLake(bounds, num_c, filter_c, & has_h2o = .false. end if else - has_h2o = .true. +!KO has_h2o = .true. +!KO + if (j <= nlevgrnd) then + has_h2o = .true. + end if +!KO end if if (has_h2o) then @@ -719,6 +733,9 @@ subroutine AccumulateSoilHeatNonLake(bounds, num_c, filter_c, & ! !LOCAL VARIABLES: integer :: fc integer :: l, c, j +!KO + integer :: nlev ! greater of nlevgrnd and nlevurb +!KO logical :: has_h2o ! whether this point potentially has water to add real(r8) :: soil_heat_liquid(bounds%begc:bounds%endc) ! sum of heat content: liquid water in soil, excluding latent heat [J/m^2] @@ -755,7 +772,11 @@ subroutine AccumulateSoilHeatNonLake(bounds, num_c, filter_c, & soil_latent_heat_liquid(c) = 0._r8 end do - do j = 1, nlevgrnd +!KO + nlev = max0(nlevgrnd,nlevurb) + do j = 1, nlev +!KO +!KO do j = 1, nlevgrnd do fc = 1, num_c c = filter_c(fc) l = col%landunit(c) @@ -777,17 +798,23 @@ subroutine AccumulateSoilHeatNonLake(bounds, num_c, filter_c, & end if else - has_h2o = .true. +!KO + if (j <= nlevgrnd) then +!KO + has_h2o = .true. - if (col%itype(c) == icol_road_imperv .and. j <= nlev_improad(l)) then - soil_heat_dry_mass(c) = soil_heat_dry_mass(c) + & - TempToHeat(temp = t_soisno(c,j), cv = (cv_improad(l,j) * dz(c,j))) - else if (lun%itype(l) /= istwet .and. lun%itype(l) /= istice_mec) then - ! Note that this also includes impervious roads below nlev_improad (where - ! we have soil) - soil_heat_dry_mass(c) = soil_heat_dry_mass(c) + & - TempToHeat(temp = t_soisno(c,j), cv = (csol(c,j)*(1-watsat(c,j))*dz(c,j))) + if (col%itype(c) == icol_road_imperv .and. j <= nlev_improad(l)) then + soil_heat_dry_mass(c) = soil_heat_dry_mass(c) + & + TempToHeat(temp = t_soisno(c,j), cv = (cv_improad(l,j) * dz(c,j))) + else if (lun%itype(l) /= istwet .and. lun%itype(l) /= istice_mec) then + ! Note that this also includes impervious roads below nlev_improad (where + ! we have soil) + soil_heat_dry_mass(c) = soil_heat_dry_mass(c) + & + TempToHeat(temp = t_soisno(c,j), cv = (csol(c,j)*(1-watsat(c,j))*dz(c,j))) + end if +!KO end if +!KO end if if (has_h2o) then diff --git a/src/biogeophys/TridiagonalMod.F90 b/src/biogeophys/TridiagonalMod.F90 index 68dbd71cce..bfb296f50d 100644 --- a/src/biogeophys/TridiagonalMod.F90 +++ b/src/biogeophys/TridiagonalMod.F90 @@ -28,8 +28,8 @@ subroutine Tridiagonal (bounds, lbj, ubj, jtop, numf, filter, a, b, c, r, u) ! !USES: use shr_kind_mod , only : r8 => shr_kind_r8 use shr_log_mod , only : errMsg => shr_log_errMsg - use clm_varpar , only : nlevurb - use column_varcon , only : icol_roof, icol_sunwall, icol_shadewall +!KO use clm_varpar , only : nlevurb +!KO use column_varcon , only : icol_roof, icol_sunwall, icol_shadewall use clm_varctl , only : iulog use decompMod , only : bounds_type use ColumnType , only : col @@ -39,8 +39,12 @@ subroutine Tridiagonal (bounds, lbj, ubj, jtop, numf, filter, a, b, c, r, u) type(bounds_type), intent(in) :: bounds integer , intent(in) :: lbj, ubj ! lbinning and ubing level indices integer , intent(in) :: jtop( bounds%begc: ) ! top level for each column [col] - integer , intent(in) :: numf ! filter dimension - integer , intent(in) :: filter(:) ! filter +!KO integer , intent(in) :: numf ! filter dimension +!KO integer , intent(in) :: filter(:) ! filter +!KO + integer , intent(in) :: numf ! filter dimension (should not include hydrologically inactive points) + integer , intent(in) :: filter(:) ! filter (should not include hydrologically inactive points) +!KO real(r8), intent(in) :: a( bounds%begc: , lbj: ) ! "a" left off diagonal of tridiagonal matrix [col, j] real(r8), intent(in) :: b( bounds%begc: , lbj: ) ! "b" diagonal column for tridiagonal matrix [col, j] real(r8), intent(in) :: c( bounds%begc: , lbj: ) ! "c" right off diagonal tridiagonal matrix [col, j] @@ -70,46 +74,62 @@ subroutine Tridiagonal (bounds, lbj, ubj, jtop, numf, filter, a, b, c, r, u) do j = lbj, ubj do fc = 1,numf ci = filter(fc) - if ((col%itype(ci) == icol_sunwall .or. col%itype(ci) == icol_shadewall & - .or. col%itype(ci) == icol_roof) .and. j <= nlevurb) then - if (j >= jtop(ci)) then - if (j == jtop(ci)) then - u(ci,j) = r(ci,j) / bet(ci) - else - gam(ci,j) = c(ci,j-1) / bet(ci) - bet(ci) = b(ci,j) - a(ci,j) * gam(ci,j) - u(ci,j) = (r(ci,j) - a(ci,j)*u(ci,j-1)) / bet(ci) - end if - end if - else if (col%itype(ci) /= icol_sunwall .and. col%itype(ci) /= icol_shadewall & - .and. col%itype(ci) /= icol_roof) then - if (j >= jtop(ci)) then - if (j == jtop(ci)) then - u(ci,j) = r(ci,j) / bet(ci) - else - gam(ci,j) = c(ci,j-1) / bet(ci) - bet(ci) = b(ci,j) - a(ci,j) * gam(ci,j) - u(ci,j) = (r(ci,j) - a(ci,j)*u(ci,j-1)) / bet(ci) - end if +!KO if ((col%itype(ci) == icol_sunwall .or. col%itype(ci) == icol_shadewall & +!KO .or. col%itype(ci) == icol_roof) .and. j <= nlevurb) then +!KO if (j >= jtop(ci)) then +!KO if (j == jtop(ci)) then +!KO u(ci,j) = r(ci,j) / bet(ci) +!KO else +!KO gam(ci,j) = c(ci,j-1) / bet(ci) +!KO bet(ci) = b(ci,j) - a(ci,j) * gam(ci,j) +!KO u(ci,j) = (r(ci,j) - a(ci,j)*u(ci,j-1)) / bet(ci) +!KO end if +!KO end if +!KO else if (col%itype(ci) /= icol_sunwall .and. col%itype(ci) /= icol_shadewall & +!KO .and. col%itype(ci) /= icol_roof) then +!KO if (j >= jtop(ci)) then +!KO if (j == jtop(ci)) then +!KO u(ci,j) = r(ci,j) / bet(ci) +!KO else +!KO gam(ci,j) = c(ci,j-1) / bet(ci) +!KO bet(ci) = b(ci,j) - a(ci,j) * gam(ci,j) +!KO u(ci,j) = (r(ci,j) - a(ci,j)*u(ci,j-1)) / bet(ci) +!KO end if +!KO end if +!KO end if +!KO + if (j >= jtop(ci)) then + if (j == jtop(ci)) then + u(ci,j) = r(ci,j) / bet(ci) + else + gam(ci,j) = c(ci,j-1) / bet(ci) + bet(ci) = b(ci,j) - a(ci,j) * gam(ci,j) + u(ci,j) = (r(ci,j) - a(ci,j)*u(ci,j-1)) / bet(ci) end if end if +!KO end do end do do j = ubj-1,lbj,-1 do fc = 1,numf ci = filter(fc) - if ((col%itype(ci) == icol_sunwall .or. col%itype(ci) == icol_shadewall & - .or. col%itype(ci) == icol_roof) .and. j <= nlevurb-1) then - if (j >= jtop(ci)) then - u(ci,j) = u(ci,j) - gam(ci,j+1) * u(ci,j+1) - end if - else if (col%itype(ci) /= icol_sunwall .and. col%itype(ci) /= icol_shadewall & - .and. col%itype(ci) /= icol_roof) then - if (j >= jtop(ci)) then - u(ci,j) = u(ci,j) - gam(ci,j+1) * u(ci,j+1) - end if +!KO if ((col%itype(ci) == icol_sunwall .or. col%itype(ci) == icol_shadewall & +!KO .or. col%itype(ci) == icol_roof) .and. j <= nlevurb-1) then +!KO if (j >= jtop(ci)) then +!KO u(ci,j) = u(ci,j) - gam(ci,j+1) * u(ci,j+1) +!KO end if +!KO else if (col%itype(ci) /= icol_sunwall .and. col%itype(ci) /= icol_shadewall & +!KO .and. col%itype(ci) /= icol_roof) then +!KO if (j >= jtop(ci)) then +!KO u(ci,j) = u(ci,j) - gam(ci,j+1) * u(ci,j+1) +!KO end if +!KO end if +!KO + if (j >= jtop(ci)) then + u(ci,j) = u(ci,j) - gam(ci,j+1) * u(ci,j+1) end if +!KO end do end do diff --git a/src/biogeophys/UrbanParamsType.F90 b/src/biogeophys/UrbanParamsType.F90 index 4b4187c3af..2e3c7c4766 100644 --- a/src/biogeophys/UrbanParamsType.F90 +++ b/src/biogeophys/UrbanParamsType.F90 @@ -122,8 +122,11 @@ subroutine Init(this, bounds) ! ! !USES: use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=) - use clm_varpar , only : nlevcan, nlevcan, numrad, nlevgrnd, nlevurb - use clm_varpar , only : nlevsoi, nlevgrnd +!KO use clm_varpar , only : nlevcan, nlevcan, numrad, nlevgrnd, nlevurb +!KO use clm_varpar , only : nlevsoi, nlevgrnd +!KO + use clm_varpar , only : numrad, nlevurb +!KO use clm_varctl , only : use_vancouver, use_mexicocity use clm_varcon , only : vkc use column_varcon , only : icol_roof, icol_sunwall, icol_shadewall @@ -195,7 +198,7 @@ subroutine Init(this, bounds) do l = bounds%begl,bounds%endl - ! "0" refers to urban wall/roof surface and "nlevsoi" refers to urban wall/roof bottom +!KO ! "0" refers to urban wall/roof surface and "nlevsoi" refers to urban wall/roof bottom if (lun%urbpoi(l)) then g = lun%gridcell(l) diff --git a/src/main/clm_driver.F90 b/src/main/clm_driver.F90 index 4222a332a8..3df74fec7f 100644 --- a/src/main/clm_driver.F90 +++ b/src/main/clm_driver.F90 @@ -540,6 +540,9 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, ro call CanopyTemperature(bounds_clump, & filter(nc)%num_nolakec, filter(nc)%nolakec, & filter(nc)%num_nolakep, filter(nc)%nolakep, & +!KO + filter(nc)%num_urbanc, filter(nc)%urbanc, & +!KO clm_fates, & atm2lnd_inst, canopystate_inst, soilstate_inst, frictionvel_inst, & water_inst%waterstatebulk_inst, water_inst%waterdiagnosticbulk_inst, &