diff --git a/build/source/engine/computFlux.f90 b/build/source/engine/computFlux.f90
index 75b8fc486fe20ebab28f9ba7e95b4f8892bb62c5..1c7e5cdf0a84ba67139ce638ef0aec3338e4cdb6 100755
--- a/build/source/engine/computFlux.f90
+++ b/build/source/engine/computFlux.f90
@@ -521,6 +521,7 @@ contains
endif ! if computing energy fluxes throughout the snow+soil domain
+! print*, "After ssdNRGFlux call ", flux_data%var(iLookFLUX%iLayerLiqFluxSoil)%dat
! *****
! * CALCULATE THE LIQUID FLUX THROUGH VEGETATION...
@@ -675,6 +676,9 @@ contains
err,cmessage) ! intent(out): error control
if(err/=0)then; message=trim(message)//trim(cmessage); return; endif
+! print*, "After soil Liq call ", flux_data%var(iLookFLUX%iLayerLiqFluxSoil)%dat
+
+
! calculate net liquid water fluxes for each soil layer (s-1)
do iLayer=1,nSoil
mLayerLiqFluxSoil(iLayer) = -(iLayerLiqFluxSoil(iLayer) - iLayerLiqFluxSoil(iLayer-1))/mLayerDepth(iLayer+nSnow)
diff --git a/build/source/engine/coupled_em.f90 b/build/source/engine/coupled_em.f90
index 4fecab68ae997ebb59baba35986baae0583b7219..de184bfdf69ef4679b9053bcf18412e36fdeae80 100755
--- a/build/source/engine/coupled_em.f90
+++ b/build/source/engine/coupled_em.f90
@@ -99,7 +99,7 @@ contains
! ************************************************************************************************
! public subroutine coupled_em: run the coupled energy-mass model for one timestep
! ************************************************************************************************
- subroutine coupled_em(&
+subroutine coupled_em(&
! model control
indxHRU, & ! intent(in): hruId
dt_init, & ! intent(inout): used to initialize the size of the sub-step
@@ -763,8 +763,6 @@ contains
! check for all errors (error recovery within opSplittin)
if(err/=0)then; err=20; message=trim(message)//trim(cmessage); return; end if
- !print*, 'completed step'
- !print*, 'PAUSE: '; read(*,*)
! process the flag for too much melt
diff --git a/build/source/engine/opSplittin.f90 b/build/source/engine/opSplittin.f90
index d0612d2f277d97ac14b42fc63e0ef7c64056a6d2..38a175d7e84c7a86f83314763bbecefe6a4e0902 100755
--- a/build/source/engine/opSplittin.f90
+++ b/build/source/engine/opSplittin.f90
@@ -798,9 +798,9 @@ contains
if(err>0) return
endif ! (check for errors)
- !print*, trim(message)//'after varSubstep: scalarSnowDrainage = ', flux_data%var(iLookFLUX%scalarSnowDrainage)%dat
- !print*, trim(message)//'after varSubstep: iLayerLiqFluxSnow = ', flux_data%var(iLookFLUX%iLayerLiqFluxSnow)%dat
- !print*, trim(message)//'after varSubstep: iLayerLiqFluxSoil = ', flux_data%var(iLookFLUX%iLayerLiqFluxSoil)%dat
+ ! print*, trim(message)//'after varSubstep: scalarSnowDrainage = ', flux_data%var(iLookFLUX%scalarSnowDrainage)%dat
+ ! print*, trim(message)//'after varSubstep: iLayerLiqFluxSnow = ', flux_data%var(iLookFLUX%iLayerLiqFluxSnow)%dat
+ ! print*, trim(message)//'after varSubstep: iLayerLiqFluxSoil = ', flux_data%var(iLookFLUX%iLayerLiqFluxSoil)%dat
! check
!if(ixSolution==scalar)then
diff --git a/build/source/engine/soilLiqFlx.f90 b/build/source/engine/soilLiqFlx.f90
index 52eb06ce616794ac22edbfebc34d97c1d01e5745..e5c57fe148c4f5f38c32c4acdaa789f9ebaaebed 100755
--- a/build/source/engine/soilLiqFlx.f90
+++ b/build/source/engine/soilLiqFlx.f90
@@ -529,7 +529,7 @@ contains
! output: error control
err,cmessage) ! intent(out): error control
if(err/=0)then; message=trim(message)//trim(cmessage); return; end if
-
+ ! print*, "scalarGroundEvaporation =", scalarGroundEvaporation
! include base soil evaporation as the upper boundary flux
iLayerLiqFluxSoil(0) = scalarGroundEvaporation/iden_water + scalarSurfaceInfiltration
@@ -542,7 +542,7 @@ contains
end select
end if
- !write(*,'(a,1x,10(f30.15))') 'scalarRainPlusMelt, scalarSurfaceInfiltration = ', scalarRainPlusMelt, scalarSurfaceInfiltration
+ ! write(*,'(a,1x,10(f30.15))') 'scalarRainPlusMelt, scalarSurfaceInfiltration = ', scalarRainPlusMelt, scalarSurfaceInfiltration
end do ! (looping through different flux calculations -- one or multiple calls depending if desire for numerical or analytical derivatives)
@@ -550,7 +550,7 @@ contains
if(deriv_desired .and. ixDerivMethod==numerical)then
dq_dHydStateBelow(0) = (scalarFlux_dStateBelow - scalarFlux)/dx ! change in surface flux w.r.t. change in the soil moisture in the top soil layer (m s-1)
end if
- !print*, 'scalarSurfaceInfiltration, iLayerLiqFluxSoil(0) = ', scalarSurfaceInfiltration, iLayerLiqFluxSoil(0)
+! print*, 'scalarSurfaceInfiltration, iLayerLiqFluxSoil(0) = ', scalarSurfaceInfiltration, iLayerLiqFluxSoil(0)
!print*, '(ixDerivMethod==numerical), dq_dHydStateBelow(0) = ', (ixDerivMethod==numerical), dq_dHydStateBelow(0)
!pause
diff --git a/build/source/engine/varSubstep.f90 b/build/source/engine/varSubstep.f90
index f9d91585a69e30c2e5c626818b8ef980c92dae10..0d9c0c5cf2d221d0f4a054df1fe3af295143665f 100755
--- a/build/source/engine/varSubstep.f90
+++ b/build/source/engine/varSubstep.f90
@@ -815,7 +815,7 @@ contains
!write(*,'(a,1x,f20.10)') 'dt = ', dt
!write(*,'(a,1x,f20.10)') 'soilBalance0 = ', soilBalance0
!write(*,'(a,1x,f20.10)') 'soilBalance1 = ', soilBalance1
- !write(*,'(a,1x,f20.10)') 'vertFlux = ', vertFlux
+ ! write(*,'(a,1x,f20.10)') 'vertFlux = ', vertFlux
!write(*,'(a,1x,f20.10)') 'tranSink = ', tranSink
!write(*,'(a,1x,f20.10)') 'baseSink = ', baseSink
!write(*,'(a,1x,f20.10)') 'compSink = ', compSink
diff --git a/build/source/engine/vegNrgFlux.f90 b/build/source/engine/vegNrgFlux.f90
index 47bfba9a9ba9d8b7930d32e66dcbbdfc6d57d151..2e80729623d8fd510594274617ef2a047dbabf67 100755
--- a/build/source/engine/vegNrgFlux.f90
+++ b/build/source/engine/vegNrgFlux.f90
@@ -727,6 +727,7 @@ contains
! compute the total vegetation area index (leaf plus stem)
VAI = scalarLAI + scalarSAI ! vegetation area index
+
exposedVAI = scalarExposedLAI + scalarExposedSAI ! exposed vegetation area index
! compute emissivity of the canopy (-)
@@ -795,12 +796,6 @@ contains
dCanopyWetFraction_dT = 0._dp ! derivative in wetted fraction w.r.t. canopy temperature (K-1)
end if
!write(*,'(a,1x,L1,1x,f25.15,1x))') 'computeVegFlux, scalarCanopyWetFraction = ', computeVegFlux, scalarCanopyWetFraction
- !print*, 'dCanopyWetFraction_dWat = ', dCanopyWetFraction_dWat
- !print*, 'dCanopyWetFraction_dT = ', dCanopyWetFraction_dT
- !print*, 'canopyLiqTrial = ', canopyLiqTrial
- !print*, 'canopyIceTrial = ', canopyIceTrial
- !print*, 'scalarCanopyLiqMax = ', scalarCanopyLiqMax
- !print*, 'scalarCanopyIceMax = ', scalarCanopyIceMax
! *******************************************************************************************************************************************************************
! *******************************************************************************************************************************************************************
@@ -1153,6 +1148,7 @@ contains
! assign scalar resistances for un-perturbed cases
else
+
trialLeafResistance = scalarLeafResistance
trialGroundResistance = scalarGroundResistance
trialCanopyResistance = scalarCanopyResistance
@@ -1278,8 +1274,8 @@ contains
!write(*,'(a,f25.15)') 'scalarLatHeatCanopyEvap = ', scalarLatHeatCanopyEvap
!write(*,'(a,f25.15)') 'scalarLatHeatCanopyTrans = ', scalarLatHeatCanopyTrans
- !print*, 'scalarSenHeatGround = ', scalarSenHeatGround
- !print*, 'scalarLatHeatGround = ', scalarLatHeatGround
+ ! print*, 'scalarSenHeatGround = ', scalarSenHeatGround
+ ! print*, 'scalarLatHeatGround = ', scalarLatHeatGround
!notUsed_scalarCanopyStabilityCorrection ! stability correction for the canopy (-)
!notUsed_scalarGroundStabilityCorrection ! stability correction for the ground surface (-)
@@ -1427,9 +1423,9 @@ contains
scalarGroundEvaporation = scalarLatHeatGround/LH_vap
scalarSnowSublimation = 0._dp ! no sublimation from snow if no snow layers have formed
end if
- !print*, 'scalarSnowSublimation, scalarLatHeatGround = ', scalarSnowSublimation, scalarLatHeatGround
-
- !print*, 'canopyWetFraction, scalarCanopyEvaporation = ', canopyWetFraction, scalarCanopyEvaporation
+ ! print*, 'nsnow = ', nsnow
+ ! print*, 'scalarSnowSublimation, scalarLatHeatGround = ', scalarSnowSublimation, scalarLatHeatGround
+ ! print*, 'canopyWetFraction, scalarCanopyEvaporation = ', canopyWetFraction, scalarCanopyEvaporation
! *******************************************************************************************************************************************************************
! *******************************************************************************************************************************************************************
@@ -2221,7 +2217,7 @@ contains
windspdCanopyTop = windspd*windConvFactor_fv
! compute the windspeed reduction
- ! Refs: Norman et al. (Ag. Forest Met., 1995) -- citing Goudriaan (1977 manuscript "crop micrometeorology: a simulation study", Wageningen).
+ ! Refs: Norman et al. (Ag. Forest Met., 1995) -- citing Goudriaan (1977 manuscript "crop micrometeorology: a simulation study", Wageningen).
windReductionFactor = windReductionParam * exposedVAI**twoThirds * (heightCanopyTopAboveSnow - heightCanopyBottomAboveSnow)**oneThird / leafDimension**oneThird
! compute windspeed at the height z0Canopy+zeroPlaneDisplacement (m s-1)
@@ -2246,10 +2242,16 @@ contains
eddyDiffusCanopyTop = max(vkc*FrictionVelocity*(heightCanopyTopAboveSnow - zeroPlaneDisplacement), mpe)
! compute the resistance between the surface and canopy air UNDER NEUTRAL CONDITIONS (s m-1)
-
+ ! print*, ""
+ ! print*, "-windReductionFactor = ", -windReductionFactor
+ ! print*, " z0Ground = ", z0Ground
+ ! print*, " heightCanopyTopAboveSnow = ", heightCanopyTopAboveSnow
+ ! print*, " zeroPlaneDisplacement = ", zeroPlaneDisplacement
+ ! print*, " heightCanopyTopAboveSnow = ", heightCanopyTopAboveSnow
+ ! print*, " eddyDiffusCanopyTop = ", eddyDiffusCanopyTop
+ ! print*, ""
! case 1: assume exponential profile extends from the snow depth plus surface roughness length to the displacement height plus vegetation roughness
if(ixWindProfile==exponential .or. heightCanopyBottomAboveSnow<z0Ground+xTolerance)then
-
! compute the neutral ground resistance
tmp1 = exp(-windReductionFactor* z0Ground/heightCanopyTopAboveSnow)
tmp2 = exp(-windReductionFactor*(z0Canopy+zeroPlaneDisplacement)/heightCanopyTopAboveSnow)
@@ -2257,8 +2259,7 @@ contains
! case 2: logarithmic profile from snow depth plus roughness height to bottom of the canopy
! NOTE: heightCanopyBottomAboveSnow>z0Ground+xTolerance
- else
-
+ else
! compute the neutral ground resistance
! (first, component between heightCanopyBottomAboveSnow and z0Canopy+zeroPlaneDisplacement)
tmp1 = exp(-windReductionFactor* heightCanopyBottomAboveSnow/heightCanopyTopAboveSnow)
@@ -2266,7 +2267,7 @@ contains
groundResistanceNeutral = ( heightCanopyTopAboveSnow*exp(windReductionFactor) / (windReductionFactor*eddyDiffusCanopyTop) ) * (tmp1 - tmp2)
! (add log-below-canopy component)
groundResistanceNeutral = groundResistanceNeutral + (1._dp/(max(0.1_dp,windspdCanopyBottom)*vkc**2._dp))*(log(heightCanopyBottomAboveSnow/z0Ground))**2._dp
-
+
endif ! switch between exponential profile and log-below-canopy
! compute the stability correction for resistance from the ground to the canopy air space (-)
@@ -2797,7 +2798,6 @@ contains
end if
groundConductanceLH = 1._dp/(groundResistance + soilResistance) ! NOTE: soilResistance accounts for fractional snow, and =0 when snow cover is 100%
totalConductanceLH = evapConductance + transConductance + groundConductanceLH + canopyConductance
-
! check sensible heat conductance
if(totalConductanceSH < -tinyVal .or. groundConductanceSH < -tinyVal .or. canopyConductance < -tinyVal)then
message=trim(message)//'negative conductance for sensible heat'
@@ -2899,7 +2899,19 @@ contains
latHeatGround = -latHeatSubVapGround*latentHeatConstant*groundConductanceLH*(satVP_GroundTemp*soilRelHumidity - VPair) ! (positive downwards)
senHeatTotal = senHeatGround
end if
- !write(*,'(a,10(f25.15,1x))') 'latHeatGround = ', latHeatGround
+! print*, "-volHeatCapacitAir = ", volHeatCapacityAir
+! print*, "groundConductanceSH = ", groundConductanceSH
+! print*, "groundTemp = ", groundTemp
+! print*, "canairTemp = ", canairTemp
+! print*, "latHeatSubVapGround", -latHeatSubVapGround
+! print*, "latentHeatConstant", latentHeatConstant
+! print*, "groundConductanceLH", groundConductanceLH
+! print*, "satVP_GroundTemp", satVP_GroundTemp
+! print*, "soilRelHumidity", soilRelHumidity
+! print*, "VP_CanopyAir", VP_CanopyAir
+
+
+! write(*,'(a,10(f25.15,1x))') 'latHeatGround = ', latHeatGround
! compute latent heat flux from the canopy air space to the atmosphere
! NOTE: VP_CanopyAir is a diagnostic variable
diff --git a/build/source/noah-mp/module_sf_noahmplsm.F b/build/source/noah-mp/module_sf_noahmplsm.F
index 8cf25a894d2ac85d4ee1ed1e6d9d3661684f4718..2566d35b215ed6db86e1c1dd78472db03c7ee06a 100755
--- a/build/source/noah-mp/module_sf_noahmplsm.F
+++ b/build/source/noah-mp/module_sf_noahmplsm.F
@@ -482,23 +482,23 @@ END MODULE NOAHMP_RAD_PARAMETERS
! ==================================================================================================
MODULE NOAHMP_ROUTINES
- use nrtype
- USE NOAHMP_GLOBALS
- IMPLICIT NONE
-
- public :: noahmp_options
- public :: REDPRM
-
- public :: PHENOLOGY
- public :: RADIATION
- private :: ALBEDO
- private :: SNOW_AGE
- private :: SNOWALB_BATS
- private :: SNOWALB_CLASS
- private :: GROUNDALB
- public :: TWOSTREAM
- public :: STOMATA
- public :: CANRES
+ use nrtype
+ USE NOAHMP_GLOBALS
+ IMPLICIT NONE
+
+ public :: noahmp_options
+ public :: REDPRM
+
+ public :: PHENOLOGY
+ public :: RADIATION
+ private :: ALBEDO
+ private :: SNOW_AGE
+ private :: SNOWALB_BATS
+ private :: SNOWALB_CLASS
+ private :: GROUNDALB
+ public :: TWOSTREAM
+ public :: STOMATA
+ public :: CANRES
contains
!
@@ -546,7 +546,6 @@ contains
REAL(rkind) :: WT1,WT2 !interpolation weights
REAL(rkind) :: T !current month (1.00, ..., 12.00)
! --------------------------------------------------------------------------------------------------
-
! Interpolate monthly SAI and LAI to daily values
IF ( DVEG == 1 .or. DVEG == 3 .or. DVEG == 4 ) THEN
@@ -569,46 +568,77 @@ contains
LAI = WT1*LAIM(VEGTYP,IT1) + WT2*LAIM(VEGTYP,IT2)
SAI = WT1*SAIM(VEGTYP,IT1) + WT2*SAIM(VEGTYP,IT2)
ENDIF
-
- ! Realism check: no leaves without stems
- IF (SAI == 0.0) LAI = 0.0
+ IF (SAI < 0.05) SAI = 0.0 ! MB: SAI CHECK
+ IF (LAI < 0.05 .OR. SAI == 0.0) LAI = 0.0 ! MB: LAI CHECK
- ! Realism check: warn about no stems for vegetated land classes
- IF ( (SAI == 0.0) .and. ( VEGTYP /= ISWATER ) .and. ( VEGTYP /= ISBARREN ) .and. ( VEGTYP /= ISSNOW ) .and. ( VEGTYP /= ISURBAN) ) THEN
- write(*,'(A,I3,A)') ' WARNING: module_sf_noahmplsm/PHENOLOGY: Stem Area Index (SAI) = 0.0 may be unrealistic for vegetation type ',VEGTYP,'. Continuing.'
- ENDIF
-
- ! Realism check: no vegetation should exist on certain land classes
IF ( ( VEGTYP == ISWATER ) .OR. ( VEGTYP == ISBARREN ) .OR. ( VEGTYP == ISSNOW ) .or. ( VEGTYP == ISURBAN) ) THEN
LAI = 0.
SAI = 0.
ENDIF
+ ! Realism check: no leaves without stems
+ ! IF (SAI == 0.0) LAI = 0.0
+
+ ! ! Realism check: warn about no stems for vegetated land classes
+ ! IF ( (SAI == 0.0) .and. ( VEGTYP /= ISWATER ) .and. ( VEGTYP /= ISBARREN ) .and. ( VEGTYP /= ISSNOW ) .and. ( VEGTYP /= ISURBAN) ) THEN
+ ! write(*,'(A,I3,A)') ' WARNING: module_sf_noahmplsm/PHENOLOGY: Stem Area Index (SAI) = 0.0 may be unrealistic for vegetation type ',VEGTYP,'. Continuing.'
+ ! ENDIF
+
+ ! ! Realism check: no vegetation should exist on certain land classes
+ ! IF ( ( VEGTYP == ISWATER ) .OR. ( VEGTYP == ISBARREN ) .OR. ( VEGTYP == ISSNOW ) .or. ( VEGTYP == ISURBAN) ) THEN
+ ! LAI = 0.
+ ! SAI = 0.
+ ! ENDIF
!buried by snow
+ DB = MIN( MAX(SNOWH - HVB(VEGTYP),0.), HVT(VEGTYP)-HVB(VEGTYP) )
+ FB = DB / MAX(1.E-06,HVT(VEGTYP)-HVB(VEGTYP))
+ !print*, 'HVB(VEGTYP), HVT(VEGTYP), DB, FB = ', HVB(VEGTYP), HVT(VEGTYP), DB, FB
+
+ IF(HVT(VEGTYP)> 0. .AND. HVT(VEGTYP) <= 0.5) THEN
+ SNOWHC = HVT(VEGTYP)*EXP(-SNOWH/0.1)
+ IF (SNOWH>SNOWHC) THEN
+ FB = 1.
+ ELSE
+ FB = SNOWH/SNOWHC
+ ENDIF
+ ENDIF
- DB = MIN( MAX(SNOWH - HVB(VEGTYP),0.), HVT(VEGTYP)-HVB(VEGTYP) )
- FB = DB / MAX(1.E-06,HVT(VEGTYP)-HVB(VEGTYP))
- !print*, 'HVB(VEGTYP), HVT(VEGTYP), DB, FB = ', HVB(VEGTYP), HVT(VEGTYP), DB, FB
-
- IF(HVT(VEGTYP)> 0. .AND. HVT(VEGTYP) <= 0.5) THEN
- SNOWHC = HVT(VEGTYP)*EXP(-SNOWH/0.1)
- IF (SNOWH>SNOWHC) THEN
- FB = 1.
- ELSE
- FB = SNOWH/SNOWHC
- ENDIF
- ENDIF
+ ELAI = LAI*(1.-FB)
+ ESAI = SAI*(1.-FB)
+ IF (ESAI < 0.05) ESAI = 0.0 ! MB: ESAI CHECK
+ IF (ELAI < 0.05 .OR. ESAI == 0.0) ELAI = 0.0 ! MB: LAI CHECK
+
+ IF (TV .GT. TMIN(VEGTYP)) THEN
+ IGS = 1.
+ ELSE
+ IGS = 0.
+ ENDIF
- ELAI = LAI*(1.-FB)
- ESAI = SAI*(1.-FB)
+ HTOP = HVT(VEGTYP)
- IF (TV .GT. TMIN(VEGTYP)) THEN
- IGS = 1.
- ELSE
- IGS = 0.
- ENDIF
+ ! DB = MIN( MAX(SNOWH - HVB(VEGTYP),0.), HVT(VEGTYP)-HVB(VEGTYP) )
+ ! FB = DB / MAX(1.E-06,HVT(VEGTYP)-HVB(VEGTYP))
+ ! !print*, 'HVB(VEGTYP), HVT(VEGTYP), DB, FB = ', HVB(VEGTYP), HVT(VEGTYP), DB, FB
+
+ ! IF(HVT(VEGTYP)> 0. .AND. HVT(VEGTYP) <= 0.5) THEN
+ ! SNOWHC = HVT(VEGTYP)*EXP(-SNOWH/0.1)
+ ! IF (SNOWH>SNOWHC) THEN
+ ! FB = 1.
+ ! ELSE
+ ! FB = SNOWH/SNOWHC
+ ! ENDIF
+ ! ENDIF
+
+ ! ELAI = LAI*(1.-FB)
+ ! ESAI = SAI*(1.-FB)
+
+ ! IF (TV .GT. TMIN(VEGTYP)) THEN
+ ! IGS = 1.
+ ! ELSE
+ ! IGS = 0.
+ ! ENDIF
- HTOP = HVT(VEGTYP)
+ ! HTOP = HVT(VEGTYP)
END SUBROUTINE PHENOLOGY
! ==================================================================================================