Noah-MP LSM M-O stability paramaters

[apattant]

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

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?

 

[Ming Chen]

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.

 

[Ming Chen]

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