From 994c2999c0e96bf043b81430155a0086de166b03 Mon Sep 17 00:00:00 2001 From: Philippe Blain Date: Sun, 23 Apr 2023 12:03:07 -0400 Subject: [PATCH] icepack_atmo: use separate stability factors for momentum and scalars In 37f2a17 (Add support for staggered atmospheric levels (#364), 2021-05-28), we introduced the subroutine compute_stability_function, which computes the stability function during the iteration in 'icepack_amto::atmo_boundary_layer'. This subroutine also returns the unit step function 'stable' at 'hol', since it is used at the end of the computation to compute the diagnostic temperature and humidity. Note that we overwrite 'stable' in the second call to compute_stability_function, but this is OK since we call it first for momentum and then for scalars, and it is the value for scalars that we want to use later for diagnostic temperature and humdity. In a following commit we will want to compute the diagnostic _velocity_ in a different way than what is done now, and we will need the unit step function computed at the momentum level (which is currently overwritten). In preparation for that change, add a new variable 'stablem' for the unit step function at the momentum level, and rename the existing 'stable' variable to 'stables' (for scalars). Adjust the computation of 'psix2' accordingly. --- columnphysics/icepack_atmo.F90 | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/columnphysics/icepack_atmo.F90 b/columnphysics/icepack_atmo.F90 index f1a661ac..02970031 100644 --- a/columnphysics/icepack_atmo.F90 +++ b/columnphysics/icepack_atmo.F90 @@ -166,7 +166,8 @@ subroutine atmo_boundary_layer (sfctype, & cp , & ! specific heat of moist air holm , & ! H (at zlvl ) over L hols , & ! H (at zlvs ) over L (if zlvs present) - stable, & ! stability factor + stablem, & ! stability factor (momentum) + stables, & ! stability factor (scalars) cpvir , & ! defined as cp_wv/cp_air - 1. psixh ! stability function at zlvl (at zlvs if present) (heat and water) @@ -297,8 +298,8 @@ subroutine atmo_boundary_layer (sfctype, & hols = holm endif - call compute_stability_function('momentum', holm, stable, psimh) - call compute_stability_function('scalar' , hols, stable, psixh) + call compute_stability_function('momentum', holm, stablem, psimh) + call compute_stability_function('scalar' , hols, stables, psixh) ! shift all coeffs to measurement height and stability rd = rdn / (c1+rdn/vonkar*(alzm-psimh)) @@ -381,7 +382,7 @@ subroutine atmo_boundary_layer (sfctype, & else hols = hols*zTrf/zlvl endif - psix2 = -c5*hols*stable + (c1-stable)*psi_scalar_unstable(hols) + psix2 = -c5*hols*stables + (c1-stables)*psi_scalar_unstable(hols) fac = (rh/vonkar) & * (alzs + al2 - psixh + psix2) Tref = potT - delt*fac