Skip to content

[DO NOT MERGE] stability corrections #602

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
ijaguirre wants to merge 15 commits into
base: develop_sundials
Choose a base branch
from
Open

[DO NOT MERGE] stability corrections #602

Changes from all commits
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
b3aac71
implement mHeight Diff and wind speed Diff updated in the aStability …
ijaguirre Jul 4, 2025
d9bb8c0
adding windspdCanopyRef a new variable real kind
ijaguirre Jul 4, 2025
09624fd
testing only mHeight change
ijaguirre Jul 4, 2025
f834127
Changes following Martyn instructions: update the aStability call inp…
ijaguirre Jul 4, 2025
3304b4f
repeat reference Height in the code
ijaguirre Jul 4, 2025
9be4fd4
printing changes
ijaguirre Jul 4, 2025
c6ae598
also print wind speed
ijaguirre Jul 4, 2025
caf2e86
new variables heighDiff and windDiff, plus print of RiBulkCanopy
ijaguirre Jul 4, 2025
a223818
typo
ijaguirre Jul 4, 2025
002ba54
more printing
ijaguirre Jul 4, 2025
cf592dc
testing v8 inside the function
ijaguirre Jul 14, 2025
0b978a8
add more prints
ijaguirre Jul 14, 2025
ea571bf
stability test chaging the canopyExNeut
ijaguirre Jul 14, 2025
6ffe222
print friction velocity
ijaguirre Jul 14, 2025
fb6878c
remove prints
ijaguirre Jul 14, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
42 changes: 35 additions & 7 deletions build/source/engine/vegNrgFlux.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1509,6 +1509,10 @@ 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) :: 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/'
Expand Down Expand Up @@ -1565,19 +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



! -----------------------------------------------------------------------------------------------------------------------------------------
! -----------------------------------------------------------------------------------------------------------------------------------------
! * 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 (-)

call aStability(&
! input
ixStability, & ! input: choice of stability function
! input: forcing data, diagnostic and state variables
mHeight, & ! input: measurement height (m)
mHeight, & ! 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, & ! 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
Expand All @@ -1591,21 +1597,30 @@ 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 ! 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)
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).
Expand Down Expand Up @@ -1633,6 +1648,19 @@ subroutine aeroResist(&
! Note: max used to avoid dividing by zero
eddyDiffusCanopyTop = max(vkc*FrictionVelocity*(heightCanopyTopAboveSnow - zeroPlaneDisplacement), mpe)


!print*, windspd, windspd, RiBulkCanopy, 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
!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)
! 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
Expand Down