From b3aac7198fce8a7edb4ce8bb4ede7ee3bd6868db Mon Sep 17 00:00:00 2001 From: Ignacio Aguirre Belmar Date: Fri, 4 Jul 2025 08:45:25 -0600 Subject: [PATCH 01/15] implement mHeight Diff and wind speed Diff updated in the aStability call --- build/source/engine/vegNrgFlux.f90 | 23 ++++++++++++++++------- 1 file changed, 16 insertions(+), 7 deletions(-) diff --git a/build/source/engine/vegNrgFlux.f90 b/build/source/engine/vegNrgFlux.f90 index aa36df5ad..1c8847334 100644 --- a/build/source/engine/vegNrgFlux.f90 +++ b/build/source/engine/vegNrgFlux.f90 @@ -1509,6 +1509,8 @@ subroutine aeroResist(& real(rkind) :: singleLeafConductance ! leaf boundary layer conductance (m s-1) real(rkind) :: canopyLeafConductance ! leaf boundary layer conductance -- scaled up to the canopy (m s-1) real(rkind) :: leaf2CanopyScaleFactor ! factor to scale from the leaf to the canopy [m s-(1/2)] + real(rkind) :: mHeightDiff ! difference between measurement height and reference height (m) + real(rkind) :: windspdDiff ! wind speed at the top of the canopy under neutral conditions (m s-1) ! ----------------------------------------------------------------------------------------------------------------------------------------- ! initialize error control err=0; message='aeroResist/' @@ -1566,6 +1568,9 @@ subroutine aeroResist(& ! check measurement height if (mHeight < zeroPlaneDisplacement+z0Canopy) then; err=20; message=trim(message)//'measurement height is below the displacement height'; return; end if + ! Above the aStability call + referenceHeight = z0Canopy+zeroPlaneDisplacement + windspdCanopyRef = windspd/log((mHeight - snowDepth - zeroPlaneDisplacement)/z0Canopy) ! ----------------------------------------------------------------------------------------------------------------------------------------- ! ----------------------------------------------------------------------------------------------------------------------------------------- ! * compute resistance for the case where the canopy is exposed @@ -1574,7 +1579,7 @@ subroutine aeroResist(& ! input ixStability, & ! input: choice of stability function ! input: forcing data, diagnostic and state variables - mHeight, & ! input: measurement height (m) + mHeight, & ! input: measurement height (m) airTemp, & ! input: air temperature above the canopy (K) canairTemp, & ! input: temperature of the canopy air space (K) windspd, & ! input: wind speed above the canopy (m s-1) @@ -1591,28 +1596,32 @@ subroutine aeroResist(& err, cmessage ) ! output: error control if (err/=0) then; message=trim(message)//trim(cmessage); return; end if + + mHeightDiff = mHeight - zeroPlaneDisplacement + windspdDiff = windspd - windspdCanopyRef + + ! compute turbulent exchange coefficient (-) - canopyExNeut = (vkc**2_i4b) / ( log((mHeight - zeroPlaneDisplacement)/z0Canopy))**2_i4b ! coefficient under conditions of neutral stability + canopyExNeut = (vkc**2_i4b) / ( log((mHeightDiff - zeroPlaneDisplacement)/z0Canopy))**2_i4b ! coefficient under conditions of neutral stability sfc2AtmExchangeCoeff_canopy = canopyExNeut*canopyStabilityCorrection ! after stability corrections ! compute the friction velocity (m s-1) - frictionVelocity = windspd * sqrt(sfc2AtmExchangeCoeff_canopy) + frictionVelocity = windspdDiff * sqrt(sfc2AtmExchangeCoeff_canopy) ! compute the above-canopy resistance (s m-1) - canopyResistance = 1._rkind/(sfc2AtmExchangeCoeff_canopy*windspd) + canopyResistance = 1._rkind/(sfc2AtmExchangeCoeff_canopy*windspdDiff) if (canopyResistance < 0._rkind) then; err=20; message=trim(message)//'canopy resistance < 0'; return; end if ! compute windspeed at the top of the canopy above snow depth (m s-1) ! NOTE: stability corrections cancel out - windConvFactor_fv = log((heightCanopyTopAboveSnow - zeroPlaneDisplacement)/z0Canopy) / log((mHeight - snowDepth - zeroPlaneDisplacement)/z0Canopy) - windspdCanopyTop = windspd*windConvFactor_fv + windConvFactor_fv = log((heightCanopyTopAboveSnow - zeroPlaneDisplacement)/z0Canopy) / log((mHeightDiff - snowDepth - zeroPlaneDisplacement)/z0Canopy) + windspdCanopyTop = windspdDiff*windConvFactor_fv ! compute the windspeed reduction ! 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) - referenceHeight = z0Canopy+zeroPlaneDisplacement windConvFactor = exp(-windReductionFactor*(1._rkind - (referenceHeight/heightCanopyTopAboveSnow))) windspdRefHeight = windspdCanopyTop*windConvFactor if(heightCanopyTopAboveSnow < referenceHeight)then; err=20; message=trim(message)//'canopy top height above snow < reference height'; return; end if From d9bb8c09bf6d90b4752e7fd259bca668b8765221 Mon Sep 17 00:00:00 2001 From: Ignacio Aguirre Belmar Date: Fri, 4 Jul 2025 08:52:02 -0600 Subject: [PATCH 02/15] adding windspdCanopyRef a new variable real kind --- build/source/engine/vegNrgFlux.f90 | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/build/source/engine/vegNrgFlux.f90 b/build/source/engine/vegNrgFlux.f90 index 1c8847334..f0c388916 100644 --- a/build/source/engine/vegNrgFlux.f90 +++ b/build/source/engine/vegNrgFlux.f90 @@ -1509,8 +1509,9 @@ subroutine aeroResist(& real(rkind) :: singleLeafConductance ! leaf boundary layer conductance (m s-1) real(rkind) :: canopyLeafConductance ! leaf boundary layer conductance -- scaled up to the canopy (m s-1) real(rkind) :: leaf2CanopyScaleFactor ! factor to scale from the leaf to the canopy [m s-(1/2)] - real(rkind) :: mHeightDiff ! difference between measurement height and reference height (m) - real(rkind) :: windspdDiff ! wind speed at the top of the canopy under neutral conditions (m s-1) + real(rkind) :: mHeightDiff ! difference between measurement height and reference height (m) Fixes aStability call + real(rkind) :: windspdDiff ! wind speed at the top of the canopy under neutral conditions (m s-1) Fixes aStability call + real(rkind) :: windspdCanopyRef ! wind speed at the reference height (m s-1) Fixes aStability call ! ----------------------------------------------------------------------------------------------------------------------------------------- ! initialize error control err=0; message='aeroResist/' @@ -1570,7 +1571,7 @@ subroutine aeroResist(& ! Above the aStability call referenceHeight = z0Canopy+zeroPlaneDisplacement - windspdCanopyRef = windspd/log((mHeight - snowDepth - zeroPlaneDisplacement)/z0Canopy) + windspdCanopyRef = windspd/log((mHeight - snowDepth - zeroPlaneDisplacement)/z0Canopy) ! This is also a new variable ! ----------------------------------------------------------------------------------------------------------------------------------------- ! ----------------------------------------------------------------------------------------------------------------------------------------- ! * compute resistance for the case where the canopy is exposed From 09624fd0c6f54397b9e48dddf78752409c25ef79 Mon Sep 17 00:00:00 2001 From: Ignacio Aguirre Belmar Date: Fri, 4 Jul 2025 09:07:27 -0600 Subject: [PATCH 03/15] testing only mHeight change --- build/source/engine/vegNrgFlux.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/build/source/engine/vegNrgFlux.f90 b/build/source/engine/vegNrgFlux.f90 index f0c388916..c991b2498 100644 --- a/build/source/engine/vegNrgFlux.f90 +++ b/build/source/engine/vegNrgFlux.f90 @@ -1571,7 +1571,7 @@ subroutine aeroResist(& ! Above the aStability call referenceHeight = z0Canopy+zeroPlaneDisplacement - windspdCanopyRef = windspd/log((mHeight - snowDepth - zeroPlaneDisplacement)/z0Canopy) ! This is also a new variable + !windspdCanopyRef = windspd/log((mHeight - snowDepth - zeroPlaneDisplacement)/z0Canopy) ! This is also a new variable ! ----------------------------------------------------------------------------------------------------------------------------------------- ! ----------------------------------------------------------------------------------------------------------------------------------------- ! * compute resistance for the case where the canopy is exposed @@ -1599,7 +1599,7 @@ subroutine aeroResist(& mHeightDiff = mHeight - zeroPlaneDisplacement - windspdDiff = windspd - windspdCanopyRef + windspdDiff = windspd ! - windspdCanopyRef ! compute turbulent exchange coefficient (-) From f8341279461bb99fa73174a602a99a3773e257c1 Mon Sep 17 00:00:00 2001 From: Ignacio Aguirre Belmar Date: Fri, 4 Jul 2025 10:32:48 -0600 Subject: [PATCH 04/15] Changes following Martyn instructions: update the aStability call input function --- build/source/engine/vegNrgFlux.f90 | 23 ++++++++--------------- 1 file changed, 8 insertions(+), 15 deletions(-) diff --git a/build/source/engine/vegNrgFlux.f90 b/build/source/engine/vegNrgFlux.f90 index c991b2498..8fe4c1641 100644 --- a/build/source/engine/vegNrgFlux.f90 +++ b/build/source/engine/vegNrgFlux.f90 @@ -1509,8 +1509,6 @@ subroutine aeroResist(& real(rkind) :: singleLeafConductance ! leaf boundary layer conductance (m s-1) real(rkind) :: canopyLeafConductance ! leaf boundary layer conductance -- scaled up to the canopy (m s-1) real(rkind) :: leaf2CanopyScaleFactor ! factor to scale from the leaf to the canopy [m s-(1/2)] - real(rkind) :: mHeightDiff ! difference between measurement height and reference height (m) Fixes aStability call - real(rkind) :: windspdDiff ! wind speed at the top of the canopy under neutral conditions (m s-1) Fixes aStability call real(rkind) :: windspdCanopyRef ! wind speed at the reference height (m s-1) Fixes aStability call ! ----------------------------------------------------------------------------------------------------------------------------------------- ! initialize error control @@ -1571,7 +1569,7 @@ subroutine aeroResist(& ! Above the aStability call referenceHeight = z0Canopy+zeroPlaneDisplacement - !windspdCanopyRef = windspd/log((mHeight - snowDepth - zeroPlaneDisplacement)/z0Canopy) ! This is also a new variable + windspdCanopyRef = windspd/log((mHeight - snowDepth - zeroPlaneDisplacement)/z0Canopy) ! This is also a new variable ! ----------------------------------------------------------------------------------------------------------------------------------------- ! ----------------------------------------------------------------------------------------------------------------------------------------- ! * compute resistance for the case where the canopy is exposed @@ -1580,10 +1578,10 @@ subroutine aeroResist(& ! input ixStability, & ! input: choice of stability function ! input: forcing data, diagnostic and state variables - mHeight, & ! input: measurement height (m) + mHeight - zeroPlaneDisplacement, & ! input: measurement height difference (m) airTemp, & ! input: air temperature above the canopy (K) canairTemp, & ! input: temperature of the canopy air space (K) - windspd, & ! input: wind speed above the canopy (m s-1) + windspd - windspdCanopyRef, & ! input: wind speed difference (m s-1) ! input: stability parameters critRichNumber, & ! input: critical value for the bulk Richardson number where turbulence ceases (-) Louis79_bparam, & ! input: parameter in Louis (1979) stability function @@ -1597,26 +1595,21 @@ subroutine aeroResist(& err, cmessage ) ! output: error control if (err/=0) then; message=trim(message)//trim(cmessage); return; end if - - mHeightDiff = mHeight - zeroPlaneDisplacement - windspdDiff = windspd ! - windspdCanopyRef - - ! compute turbulent exchange coefficient (-) - canopyExNeut = (vkc**2_i4b) / ( log((mHeightDiff - zeroPlaneDisplacement)/z0Canopy))**2_i4b ! coefficient under conditions of neutral stability + canopyExNeut = (vkc**2_i4b) / ( log((mHeight - zeroPlaneDisplacement)/z0Canopy))**2_i4b ! coefficient under conditions of neutral stability sfc2AtmExchangeCoeff_canopy = canopyExNeut*canopyStabilityCorrection ! after stability corrections ! compute the friction velocity (m s-1) - frictionVelocity = windspdDiff * sqrt(sfc2AtmExchangeCoeff_canopy) + frictionVelocity = windspd * sqrt(sfc2AtmExchangeCoeff_canopy) ! compute the above-canopy resistance (s m-1) - canopyResistance = 1._rkind/(sfc2AtmExchangeCoeff_canopy*windspdDiff) + canopyResistance = 1._rkind/(sfc2AtmExchangeCoeff_canopy*windspd) if (canopyResistance < 0._rkind) then; err=20; message=trim(message)//'canopy resistance < 0'; return; end if ! compute windspeed at the top of the canopy above snow depth (m s-1) ! NOTE: stability corrections cancel out - windConvFactor_fv = log((heightCanopyTopAboveSnow - zeroPlaneDisplacement)/z0Canopy) / log((mHeightDiff - snowDepth - zeroPlaneDisplacement)/z0Canopy) - windspdCanopyTop = windspdDiff*windConvFactor_fv + windConvFactor_fv = log((heightCanopyTopAboveSnow - zeroPlaneDisplacement)/z0Canopy) / log((mHeight - snowDepth - zeroPlaneDisplacement)/z0Canopy) + 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). From 3304b4fb32856e2058c52dd55f07f197282ef8f8 Mon Sep 17 00:00:00 2001 From: Ignacio Aguirre Belmar Date: Fri, 4 Jul 2025 10:39:57 -0600 Subject: [PATCH 05/15] repeat reference Height in the code --- build/source/engine/vegNrgFlux.f90 | 1 + 1 file changed, 1 insertion(+) diff --git a/build/source/engine/vegNrgFlux.f90 b/build/source/engine/vegNrgFlux.f90 index 8fe4c1641..e2b40161a 100644 --- a/build/source/engine/vegNrgFlux.f90 +++ b/build/source/engine/vegNrgFlux.f90 @@ -1616,6 +1616,7 @@ subroutine aeroResist(& windReductionFactor = windReductionParam * exposedVAI**twoThirds * (heightCanopyTopAboveSnow - heightCanopyBottomAboveSnow)**oneThird / leafDimension**oneThird ! compute windspeed at the height z0Canopy+zeroPlaneDisplacement (m s-1) + referenceHeight = z0Canopy+zeroPlaneDisplacement windConvFactor = exp(-windReductionFactor*(1._rkind - (referenceHeight/heightCanopyTopAboveSnow))) windspdRefHeight = windspdCanopyTop*windConvFactor if(heightCanopyTopAboveSnow < referenceHeight)then; err=20; message=trim(message)//'canopy top height above snow < reference height'; return; end if From 9be4fd43daa4b59446186e635405958a469794ea Mon Sep 17 00:00:00 2001 From: Ignacio Aguirre Belmar Date: Fri, 4 Jul 2025 10:49:20 -0600 Subject: [PATCH 06/15] printing changes --- build/source/engine/vegNrgFlux.f90 | 1 + 1 file changed, 1 insertion(+) diff --git a/build/source/engine/vegNrgFlux.f90 b/build/source/engine/vegNrgFlux.f90 index e2b40161a..aa41e3d13 100644 --- a/build/source/engine/vegNrgFlux.f90 +++ b/build/source/engine/vegNrgFlux.f90 @@ -1570,6 +1570,7 @@ subroutine aeroResist(& ! Above the aStability call referenceHeight = z0Canopy+zeroPlaneDisplacement windspdCanopyRef = windspd/log((mHeight - snowDepth - zeroPlaneDisplacement)/z0Canopy) ! This is also a new variable + print*, mHeight - zeroPlaneDisplacement, windspdCanopyRef, windspd - windspdCanopyRef ! ----------------------------------------------------------------------------------------------------------------------------------------- ! ----------------------------------------------------------------------------------------------------------------------------------------- ! * compute resistance for the case where the canopy is exposed From c6ae59846f71c15f5c47eef0bd960bbe4ea2c8ef Mon Sep 17 00:00:00 2001 From: Ignacio Aguirre Belmar Date: Fri, 4 Jul 2025 10:52:44 -0600 Subject: [PATCH 07/15] also print wind speed --- build/source/engine/vegNrgFlux.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/build/source/engine/vegNrgFlux.f90 b/build/source/engine/vegNrgFlux.f90 index aa41e3d13..bfa610011 100644 --- a/build/source/engine/vegNrgFlux.f90 +++ b/build/source/engine/vegNrgFlux.f90 @@ -1570,7 +1570,7 @@ subroutine aeroResist(& ! Above the aStability call referenceHeight = z0Canopy+zeroPlaneDisplacement windspdCanopyRef = windspd/log((mHeight - snowDepth - zeroPlaneDisplacement)/z0Canopy) ! This is also a new variable - print*, mHeight - zeroPlaneDisplacement, windspdCanopyRef, windspd - windspdCanopyRef + print*, mHeight - zeroPlaneDisplacement, windspdCanopyRef, windspd - windspdCanopyRef, windspd ! ----------------------------------------------------------------------------------------------------------------------------------------- ! ----------------------------------------------------------------------------------------------------------------------------------------- ! * compute resistance for the case where the canopy is exposed From caf2e86578646e53217b638ff5770c332cb2248f Mon Sep 17 00:00:00 2001 From: Ignacio Aguirre Belmar Date: Fri, 4 Jul 2025 13:15:29 -0600 Subject: [PATCH 08/15] new variables heighDiff and windDiff, plus print of RiBulkCanopy --- build/source/engine/vegNrgFlux.f90 | 22 +++++++++++++++------- 1 file changed, 15 insertions(+), 7 deletions(-) diff --git a/build/source/engine/vegNrgFlux.f90 b/build/source/engine/vegNrgFlux.f90 index bfa610011..5532c831d 100644 --- a/build/source/engine/vegNrgFlux.f90 +++ b/build/source/engine/vegNrgFlux.f90 @@ -1510,6 +1510,9 @@ subroutine aeroResist(& real(rkind) :: canopyLeafConductance ! leaf boundary layer conductance -- scaled up to the canopy (m s-1) real(rkind) :: leaf2CanopyScaleFactor ! factor to scale from the leaf to the canopy [m s-(1/2)] real(rkind) :: windspdCanopyRef ! wind speed at the reference height (m s-1) Fixes aStability call + real(rkind) :: heightDiff ! difference between measurement height and zero plane displacement (m) + real(rkind) :: windDiff ! difference between wind speed at measurement height and + ! ----------------------------------------------------------------------------------------------------------------------------------------- ! initialize error control err=0; message='aeroResist/' @@ -1575,14 +1578,16 @@ subroutine aeroResist(& ! ----------------------------------------------------------------------------------------------------------------------------------------- ! * compute resistance for the case where the canopy is exposed ! compute the stability correction for resistance from canopy air space to air above the canopy (-) + heightDiff = mHeight - zeroPlaneDisplacement + windDiff = windspd - windspdCanopyRef call aStability(& ! input ixStability, & ! input: choice of stability function ! input: forcing data, diagnostic and state variables - mHeight - zeroPlaneDisplacement, & ! input: measurement height difference (m) + heightDiff, & ! input: measurement height difference (m) airTemp, & ! input: air temperature above the canopy (K) canairTemp, & ! input: temperature of the canopy air space (K) - windspd - windspdCanopyRef, & ! input: wind speed difference (m s-1) + windDiff, & ! input: wind speed difference (m s-1) ! input: stability parameters critRichNumber, & ! input: critical value for the bulk Richardson number where turbulence ceases (-) Louis79_bparam, & ! input: parameter in Louis (1979) stability function @@ -1597,20 +1602,20 @@ subroutine aeroResist(& if (err/=0) then; message=trim(message)//trim(cmessage); return; end if ! compute turbulent exchange coefficient (-) - canopyExNeut = (vkc**2_i4b) / ( log((mHeight - zeroPlaneDisplacement)/z0Canopy))**2_i4b ! coefficient under conditions of neutral stability + canopyExNeut = (vkc**2_i4b) / ( log((heightDiff - zeroPlaneDisplacement)/z0Canopy))**2_i4b ! coefficient under conditions of neutral stability sfc2AtmExchangeCoeff_canopy = canopyExNeut*canopyStabilityCorrection ! after stability corrections ! compute the friction velocity (m s-1) - frictionVelocity = windspd * sqrt(sfc2AtmExchangeCoeff_canopy) + frictionVelocity = windDiff * sqrt(sfc2AtmExchangeCoeff_canopy) ! compute the above-canopy resistance (s m-1) - canopyResistance = 1._rkind/(sfc2AtmExchangeCoeff_canopy*windspd) + canopyResistance = 1._rkind/(sfc2AtmExchangeCoeff_canopy*windDiff) if (canopyResistance < 0._rkind) then; err=20; message=trim(message)//'canopy resistance < 0'; return; end if ! compute windspeed at the top of the canopy above snow depth (m s-1) ! NOTE: stability corrections cancel out - windConvFactor_fv = log((heightCanopyTopAboveSnow - zeroPlaneDisplacement)/z0Canopy) / log((mHeight - snowDepth - zeroPlaneDisplacement)/z0Canopy) - windspdCanopyTop = windspd*windConvFactor_fv + windConvFactor_fv = log((heightCanopyTopAboveSnow - zeroPlaneDisplacement)/z0Canopy) / log((heightDiff - snowDepth - zeroPlaneDisplacement)/z0Canopy) + windspdCanopyTop = windDiff*windConvFactor_fv ! compute the windspeed reduction ! Refs: Norman et al. (Ag. Forest Met., 1995) -- citing Goudriaan (1977 manuscript "crop micrometeorology: a simulation study", Wageningen). @@ -1638,6 +1643,9 @@ subroutine aeroResist(& ! Note: max used to avoid dividing by zero eddyDiffusCanopyTop = max(vkc*FrictionVelocity*(heightCanopyTopAboveSnow - zeroPlaneDisplacement), mpe) + + print*, windspd, windspd, RiBulkCanopy, canopyStabilityCorrection, + ! compute the resistance between the surface and canopy air UNDER NEUTRAL CONDITIONS (s m-1) ! case 1: assume exponential profile extends from the snow depth plus surface roughness length to the displacement height plus vegetation roughness if (ixWindProfile==exponential) then From a223818130ad562f2a993ddf3c77d27e45c04ffa Mon Sep 17 00:00:00 2001 From: Ignacio Aguirre Belmar Date: Fri, 4 Jul 2025 13:18:52 -0600 Subject: [PATCH 09/15] typo --- build/source/engine/vegNrgFlux.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/build/source/engine/vegNrgFlux.f90 b/build/source/engine/vegNrgFlux.f90 index 5532c831d..a88685dd6 100644 --- a/build/source/engine/vegNrgFlux.f90 +++ b/build/source/engine/vegNrgFlux.f90 @@ -1644,7 +1644,7 @@ subroutine aeroResist(& eddyDiffusCanopyTop = max(vkc*FrictionVelocity*(heightCanopyTopAboveSnow - zeroPlaneDisplacement), mpe) - print*, windspd, windspd, RiBulkCanopy, canopyStabilityCorrection, + print*, windspd, windspd, RiBulkCanopy, canopyStabilityCorrection ! compute the resistance between the surface and canopy air UNDER NEUTRAL CONDITIONS (s m-1) ! case 1: assume exponential profile extends from the snow depth plus surface roughness length to the displacement height plus vegetation roughness From 002ba547d1ae0f5df556458e2616eaf64bc2dad6 Mon Sep 17 00:00:00 2001 From: Ignacio Aguirre Belmar Date: Fri, 4 Jul 2025 14:04:49 -0600 Subject: [PATCH 10/15] more printing --- build/source/engine/vegNrgFlux.f90 | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/build/source/engine/vegNrgFlux.f90 b/build/source/engine/vegNrgFlux.f90 index a88685dd6..f363d21a2 100644 --- a/build/source/engine/vegNrgFlux.f90 +++ b/build/source/engine/vegNrgFlux.f90 @@ -1574,6 +1574,8 @@ subroutine aeroResist(& referenceHeight = z0Canopy+zeroPlaneDisplacement windspdCanopyRef = windspd/log((mHeight - snowDepth - zeroPlaneDisplacement)/z0Canopy) ! This is also a new variable print*, mHeight - zeroPlaneDisplacement, windspdCanopyRef, windspd - windspdCanopyRef, windspd + print* , 'referenceHeight = ', referenceHeight, 'z0Canopy = ', z0Canopy, 'zeroPlaneDisplacement = ', zeroPlaneDisplacement + print*, 'mHeight = ', mHeight, 'heightCanopyTopAboveSnow = ', heightCanopyTopAboveSnow, 'heightCanopyBottomAboveSnow = ', heightCanopyBottomAboveSnow ! ----------------------------------------------------------------------------------------------------------------------------------------- ! ----------------------------------------------------------------------------------------------------------------------------------------- ! * compute resistance for the case where the canopy is exposed @@ -1646,6 +1648,13 @@ subroutine aeroResist(& print*, windspd, windspd, RiBulkCanopy, canopyStabilityCorrection + print*, 'windspd =', windspd, ', RiBulkCanopy =', RiBulkCanopy, ', canopyStabilityCorrection =', canopyStabilityCorrection + print*, 'dCanopyStabilityCorrection_dRich =', dCanopyStabilityCorrection_dRich + print*, 'dCanopyStabilityCorrection_dAirTemp =', dCanopyStabilityCorrection_dAirTemp + print*, 'dCanopyStabilityCorrection_dCasTemp =', dCanopyStabilityCorrection_dCasTemp + print*, 'canopyExNeut =', canopyExNeut, ', sfc2AtmExchangeCoeff_canopy =', sfc2AtmExchangeCoeff_canopy + + ! compute the resistance between the surface and canopy air UNDER NEUTRAL CONDITIONS (s m-1) ! case 1: assume exponential profile extends from the snow depth plus surface roughness length to the displacement height plus vegetation roughness if (ixWindProfile==exponential) then From cf592dc2cf57b1006306cee839b0c905c4bb28f2 Mon Sep 17 00:00:00 2001 From: Ignacio Aguirre Belmar Date: Sun, 13 Jul 2025 23:02:18 -0600 Subject: [PATCH 11/15] testing v8 inside the function --- build/source/engine/vegNrgFlux.f90 | 36 +++++++++++++++++------------- 1 file changed, 20 insertions(+), 16 deletions(-) diff --git a/build/source/engine/vegNrgFlux.f90 b/build/source/engine/vegNrgFlux.f90 index f363d21a2..77d2fe4ee 100644 --- a/build/source/engine/vegNrgFlux.f90 +++ b/build/source/engine/vegNrgFlux.f90 @@ -1569,27 +1569,21 @@ subroutine aeroResist(& ! check measurement height if (mHeight < zeroPlaneDisplacement+z0Canopy) then; err=20; message=trim(message)//'measurement height is below the displacement height'; return; end if - - ! Above the aStability call - referenceHeight = z0Canopy+zeroPlaneDisplacement - windspdCanopyRef = windspd/log((mHeight - snowDepth - zeroPlaneDisplacement)/z0Canopy) ! This is also a new variable - print*, mHeight - zeroPlaneDisplacement, windspdCanopyRef, windspd - windspdCanopyRef, windspd - print* , 'referenceHeight = ', referenceHeight, 'z0Canopy = ', z0Canopy, 'zeroPlaneDisplacement = ', zeroPlaneDisplacement - print*, 'mHeight = ', mHeight, 'heightCanopyTopAboveSnow = ', heightCanopyTopAboveSnow, 'heightCanopyBottomAboveSnow = ', heightCanopyBottomAboveSnow + + ! ----------------------------------------------------------------------------------------------------------------------------------------- ! ----------------------------------------------------------------------------------------------------------------------------------------- ! * compute resistance for the case where the canopy is exposed ! compute the stability correction for resistance from canopy air space to air above the canopy (-) - heightDiff = mHeight - zeroPlaneDisplacement - windDiff = windspd - windspdCanopyRef + call aStability(& ! input ixStability, & ! input: choice of stability function ! input: forcing data, diagnostic and state variables - heightDiff, & ! input: measurement height difference (m) + mHeight, & ! input: measurement height difference (m) airTemp, & ! input: air temperature above the canopy (K) canairTemp, & ! input: temperature of the canopy air space (K) - windDiff, & ! input: wind speed difference (m s-1) + windspd, & ! input: wind speed difference (m s-1) ! input: stability parameters critRichNumber, & ! input: critical value for the bulk Richardson number where turbulence ceases (-) Louis79_bparam, & ! input: parameter in Louis (1979) stability function @@ -1603,6 +1597,13 @@ subroutine aeroResist(& err, cmessage ) ! output: error control if (err/=0) then; message=trim(message)//trim(cmessage); return; end if + ! Inside the aStability call + referenceHeight = z0Canopy+zeroPlaneDisplacement + windspdCanopyRef = windspd/log((mHeight - snowDepth - zeroPlaneDisplacement)/z0Canopy) ! This is also a new variable + + heightDiff = mHeight - zeroPlaneDisplacement + windDiff = windspd - windspdCanopyRef + ! compute turbulent exchange coefficient (-) canopyExNeut = (vkc**2_i4b) / ( log((heightDiff - zeroPlaneDisplacement)/z0Canopy))**2_i4b ! coefficient under conditions of neutral stability sfc2AtmExchangeCoeff_canopy = canopyExNeut*canopyStabilityCorrection ! after stability corrections @@ -1646,13 +1647,16 @@ subroutine aeroResist(& eddyDiffusCanopyTop = max(vkc*FrictionVelocity*(heightCanopyTopAboveSnow - zeroPlaneDisplacement), mpe) - print*, windspd, windspd, RiBulkCanopy, canopyStabilityCorrection + !print*, windspd, windspd, RiBulkCanopy, canopyStabilityCorrection print*, 'windspd =', windspd, ', RiBulkCanopy =', RiBulkCanopy, ', canopyStabilityCorrection =', canopyStabilityCorrection - print*, 'dCanopyStabilityCorrection_dRich =', dCanopyStabilityCorrection_dRich - print*, 'dCanopyStabilityCorrection_dAirTemp =', dCanopyStabilityCorrection_dAirTemp - print*, 'dCanopyStabilityCorrection_dCasTemp =', dCanopyStabilityCorrection_dCasTemp - print*, 'canopyExNeut =', canopyExNeut, ', sfc2AtmExchangeCoeff_canopy =', sfc2AtmExchangeCoeff_canopy + !print*, 'dCanopyStabilityCorrection_dRich =', dCanopyStabilityCorrection_dRich + !print*, 'dCanopyStabilityCorrection_dAirTemp =', dCanopyStabilityCorrection_dAirTemp + !print*, 'dCanopyStabilityCorrection_dCasTemp =', dCanopyStabilityCorrection_dCasTemp + !print*, 'canopyExNeut =', canopyExNeut, ', sfc2AtmExchangeCoeff_canopy =', sfc2AtmExchangeCoeff_canopy + !print*, mHeight - zeroPlaneDisplacement, windspdCanopyRef, windspd - windspdCanopyRef, windspd + !print* , 'referenceHeight = ', referenceHeight, 'z0Canopy = ', z0Canopy, 'zeroPlaneDisplacement = ', zeroPlaneDisplacement + !print*, 'mHeight = ', mHeight, 'heightCanopyTopAboveSnow = ', heightCanopyTopAboveSnow, 'heightCanopyBottomAboveSnow = ', heightCanopyBottomAboveSnow ! compute the resistance between the surface and canopy air UNDER NEUTRAL CONDITIONS (s m-1) From 0b978a873cbca286c4f02fe1beb3df938c45bff4 Mon Sep 17 00:00:00 2001 From: Ignacio Aguirre Belmar Date: Sun, 13 Jul 2025 23:41:21 -0600 Subject: [PATCH 12/15] add more prints --- build/source/engine/vegNrgFlux.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/build/source/engine/vegNrgFlux.f90 b/build/source/engine/vegNrgFlux.f90 index 77d2fe4ee..19cd77e3a 100644 --- a/build/source/engine/vegNrgFlux.f90 +++ b/build/source/engine/vegNrgFlux.f90 @@ -1653,10 +1653,10 @@ subroutine aeroResist(& !print*, 'dCanopyStabilityCorrection_dRich =', dCanopyStabilityCorrection_dRich !print*, 'dCanopyStabilityCorrection_dAirTemp =', dCanopyStabilityCorrection_dAirTemp !print*, 'dCanopyStabilityCorrection_dCasTemp =', dCanopyStabilityCorrection_dCasTemp - !print*, 'canopyExNeut =', canopyExNeut, ', sfc2AtmExchangeCoeff_canopy =', sfc2AtmExchangeCoeff_canopy + print*, 'canopyExNeut =', canopyExNeut, ', sfc2AtmExchangeCoeff_canopy =', sfc2AtmExchangeCoeff_canopy !print*, mHeight - zeroPlaneDisplacement, windspdCanopyRef, windspd - windspdCanopyRef, windspd - !print* , 'referenceHeight = ', referenceHeight, 'z0Canopy = ', z0Canopy, 'zeroPlaneDisplacement = ', zeroPlaneDisplacement - !print*, 'mHeight = ', mHeight, 'heightCanopyTopAboveSnow = ', heightCanopyTopAboveSnow, 'heightCanopyBottomAboveSnow = ', heightCanopyBottomAboveSnow + print* , 'referenceHeight = ', referenceHeight, 'z0Canopy = ', z0Canopy, 'zeroPlaneDisplacement = ', zeroPlaneDisplacement + print*, 'mHeight = ', mHeight, 'heightCanopyTopAboveSnow = ', heightCanopyTopAboveSnow, 'heightCanopyBottomAboveSnow = ', heightCanopyBottomAboveSnow ! compute the resistance between the surface and canopy air UNDER NEUTRAL CONDITIONS (s m-1) From ea571bf6aa8a79c839f0a57ad9b9123329757a9b Mon Sep 17 00:00:00 2001 From: Ignacio Aguirre Belmar Date: Mon, 14 Jul 2025 00:09:06 -0600 Subject: [PATCH 13/15] stability test chaging the canopyExNeut --- build/source/engine/vegNrgFlux.f90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/build/source/engine/vegNrgFlux.f90 b/build/source/engine/vegNrgFlux.f90 index 19cd77e3a..85336c6df 100644 --- a/build/source/engine/vegNrgFlux.f90 +++ b/build/source/engine/vegNrgFlux.f90 @@ -1605,7 +1605,9 @@ subroutine aeroResist(& windDiff = windspd - windspdCanopyRef ! compute turbulent exchange coefficient (-) - canopyExNeut = (vkc**2_i4b) / ( log((heightDiff - zeroPlaneDisplacement)/z0Canopy))**2_i4b ! coefficient under conditions of neutral stability + !canopyExNeut = (vkc**2_i4b) / ( log((heightDiff - zeroPlaneDisplacement)/z0Canopy))**2_i4b ! coefficient under conditions of neutral stability ! as should be + canopyExNeut = (vkc**2_i4b) / ( log((mHeight - zeroPlaneDisplacement)/z0Canopy))**2_i4b ! coefficient under conditions of neutral stability + sfc2AtmExchangeCoeff_canopy = canopyExNeut*canopyStabilityCorrection ! after stability corrections ! compute the friction velocity (m s-1) From 6ffe2228e5efa5ef83402094b492323b43f7ab97 Mon Sep 17 00:00:00 2001 From: Ignacio Aguirre Belmar Date: Mon, 14 Jul 2025 00:14:48 -0600 Subject: [PATCH 14/15] print friction velocity --- build/source/engine/vegNrgFlux.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/build/source/engine/vegNrgFlux.f90 b/build/source/engine/vegNrgFlux.f90 index 85336c6df..cb53d04d6 100644 --- a/build/source/engine/vegNrgFlux.f90 +++ b/build/source/engine/vegNrgFlux.f90 @@ -1651,7 +1651,7 @@ subroutine aeroResist(& !print*, windspd, windspd, RiBulkCanopy, canopyStabilityCorrection - print*, 'windspd =', windspd, ', RiBulkCanopy =', RiBulkCanopy, ', canopyStabilityCorrection =', canopyStabilityCorrection + print*, 'windspd =', windspd, ', RiBulkCanopy =', RiBulkCanopy, ', canopyStabilityCorrection =', canopyStabilityCorrection, 'FrictionVelocity =', frictionVelocity !print*, 'dCanopyStabilityCorrection_dRich =', dCanopyStabilityCorrection_dRich !print*, 'dCanopyStabilityCorrection_dAirTemp =', dCanopyStabilityCorrection_dAirTemp !print*, 'dCanopyStabilityCorrection_dCasTemp =', dCanopyStabilityCorrection_dCasTemp From fb6878ccbb002d8bce67ff78ae34ba05638c376a Mon Sep 17 00:00:00 2001 From: Ignacio Aguirre Belmar Date: Mon, 14 Jul 2025 15:26:04 -0600 Subject: [PATCH 15/15] remove prints --- build/source/engine/vegNrgFlux.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/build/source/engine/vegNrgFlux.f90 b/build/source/engine/vegNrgFlux.f90 index cb53d04d6..b2a9e5198 100644 --- a/build/source/engine/vegNrgFlux.f90 +++ b/build/source/engine/vegNrgFlux.f90 @@ -1655,10 +1655,10 @@ subroutine aeroResist(& !print*, 'dCanopyStabilityCorrection_dRich =', dCanopyStabilityCorrection_dRich !print*, 'dCanopyStabilityCorrection_dAirTemp =', dCanopyStabilityCorrection_dAirTemp !print*, 'dCanopyStabilityCorrection_dCasTemp =', dCanopyStabilityCorrection_dCasTemp - print*, 'canopyExNeut =', canopyExNeut, ', sfc2AtmExchangeCoeff_canopy =', sfc2AtmExchangeCoeff_canopy + !print*, 'canopyExNeut =', canopyExNeut, ', sfc2AtmExchangeCoeff_canopy =', sfc2AtmExchangeCoeff_canopy !print*, mHeight - zeroPlaneDisplacement, windspdCanopyRef, windspd - windspdCanopyRef, windspd - print* , 'referenceHeight = ', referenceHeight, 'z0Canopy = ', z0Canopy, 'zeroPlaneDisplacement = ', zeroPlaneDisplacement - print*, 'mHeight = ', mHeight, 'heightCanopyTopAboveSnow = ', heightCanopyTopAboveSnow, 'heightCanopyBottomAboveSnow = ', heightCanopyBottomAboveSnow + !print* , 'referenceHeight = ', referenceHeight, 'z0Canopy = ', z0Canopy, 'zeroPlaneDisplacement = ', zeroPlaneDisplacement + !print*, 'mHeight = ', mHeight, 'heightCanopyTopAboveSnow = ', heightCanopyTopAboveSnow, 'heightCanopyBottomAboveSnow = ', heightCanopyBottomAboveSnow ! compute the resistance between the surface and canopy air UNDER NEUTRAL CONDITIONS (s m-1)