Noah-MP LSM M-O stability paramaters

This post was from a previous version of the WRF&MPAS-A Support Forum. New replies have been disabled and if you have follow up questions related to this post, then please start a new thread from the forum home page.

apattant

New member
Hello,

I am looking at the M-O functions in subroutine SFCDIF1 of Noah-MP lsm and it looks to me like the businger dyer stability functions are always negative

Code:
``````IF (MOZ .LT. 0.) THEN
TMP1 = (1. - 16.*MOZ)**0.25
TMP2 = LOG((1.+TMP1*TMP1)/2.)
TMP3 = LOG((1.+TMP1)/2.)
FMNEW = 2.*TMP3 + TMP2 - 2.*ATAN(TMP1) + 1.5707963
FHNEW = 2*TMP2

! 2-meter
TMP12 = (1. - 16.*MOZ2)**0.25
TMP22 = LOG((1.+TMP12*TMP12)/2.)
TMP32 = LOG((1.+TMP12)/2.)
FM2NEW = 2.*TMP32 + TMP22 - 2.*ATAN(TMP12) + 1.5707963
FH2NEW = 2*TMP22
ELSE
FMNEW = -5.*MOZ
FHNEW = FMNEW
FM2NEW = -5.*MOZ2
FH2NEW = FM2NEW
ENDIF

! except for first iteration, weight stability factors for previous
! iteration to help avoid flip-flops from one iteration to the next

IF (ITER == 1) THEN
FM = FMNEW
FH = FHNEW
FM2 = FM2NEW
FH2 = FH2NEW
ELSE
FM = 0.5 * (FM+FMNEW)
FH = 0.5 * (FH+FHNEW)
FM2 = 0.5 * (FM2+FM2NEW)
FH2 = 0.5 * (FH2+FH2NEW)
ENDIF

! exchange coefficients

FH = MIN(FH,0.9*TMPCH)
FM = MIN(FM,0.9*TMPCM)
FH2 = MIN(FH2,0.9*TMPCH2)
FM2 = MIN(FM2,0.9*TMPCM2)

CMFM = TMPCM-FM
CHFH = TMPCH-FH
CM2FM2 = TMPCM2-FM2
CH2FH2 = TMPCH2-FH2
IF(ABS(CMFM) <= MPE) CMFM = MPE
IF(ABS(CHFH) <= MPE) CHFH = MPE
IF(ABS(CM2FM2) <= MPE) CM2FM2 = MPE
IF(ABS(CH2FH2) <= MPE) CH2FH2 = MPE
CM  = VKC*VKC/(CMFM*CMFM)
CH  = VKC*VKC/(CMFM*CHFH)
CH2  = VKC*VKC/(CM2FM2*CH2FH2)``````

It tests for stability parameter, MOZ < 0 but if it is positive or neutral it multiplies by -5. so it is essentially forcing stability parameter to be negative and making the exchange coefficients larger but never smaller. According to the Businger-Dyer relationships it should be 4.7*MOZ if MOZ >= 0

Is this a gross error or am I missing something about this subroutine for stable/nuetral conditions?

Hi,
I agree with you that this is not correct. I have contact our experts to confirm this is indeed an issue we need to address. I will get back to you once I know for sure the opinions of our experts. Thanks for raising this problem.

Hi,
I re-examine the code based on the paper below:
Businger J. A., Wyngaard J. C., Izumi Y., and Bradley E. F. , 1971: Flux profile relationship in the atmospheric surface layer, J. Atmos. Sci., 28, 181–189.

I believe the code is correct.

IF (MOZ .LT. 0.) THEN ===> this is for unstable condition
TMP1 = (1. - 16.*MOZ)**0.25
TMP2 = LOG((1.+TMP1*TMP1)/2.)
TMP3 = LOG((1.+TMP1)/2.)
FMNEW = 2.*TMP3 + TMP2 - 2.*ATAN(TMP1) + 1.5707963
FHNEW = 2*TMP2

! 2-meter
TMP12 = (1. - 16.*MOZ2)**0.25
TMP22 = LOG((1.+TMP12*TMP12)/2.)
TMP32 = LOG((1.+TMP12)/2.)
FM2NEW = 2.*TMP32 + TMP22 - 2.*ATAN(TMP12) + 1.5707963
FH2NEW = 2*TMP22
ELSE ===> this is for stable condition
FMNEW = -5.*MOZ
FHNEW = FMNEW
FM2NEW = -5.*MOZ2
FH2NEW = FM2NEW
ENDIF