Limiting MSTAR by MSTAR_CAP in MOM_energetic_PBL.F90

[ThomasNeumann2]

The limit for mstar seems to work only in the OM4 scheme but not in Reichl_h18. I don't know whether this is intended. I made a correction in code line 3255.

MOM_energetics_PBL.F90:

3226 elseif (CS%mstar_scheme == MStar_from_Ekman) then
3227
3228 if (CS%answer_date < 20190101) then
3229 ! The limit for the balance of rotation and stabilizing is f(L_Ekman,L_Obukhov)
3230 MStar_S = CS%MStar_coefsqrt(max(0.0,Buoyancy_Flux) / UStar**2 / &
3231 (Abs_Coriolis + 1.e-10US%T_to_s) )
3232 ! The limit for rotation (Ekman length) limited mixing
3233 MStar_N = CS%C_Ek * log( max( 1., UStar / (Abs_Coriolis + 1.e-10US%T_to_s) / BLD ) )
3234 else
3235 ! The limit for the balance of rotation and stabilizing is f(L_Ekman,L_Obukhov)
3236 MStar_S = CS%MSTAR_COEFsqrt(max(0.0, Buoyancy_Flux) / (UStar**2 * max(Abs_Coriolis, 1.e-20US%T_to_s)))
3237 ! The limit for rotation (Ekman length) limited mixing
3238 MStar_N = 0.0
3239 if (UStar > Abs_Coriolis * BLD) Mstar_N = CS%C_EK * log(UStar / (Abs_Coriolis * BLD))
3240 endif
3241
3242 ! Here 1.25 is about .5/von Karman, which gives the Obukhov limit.
3243 MStar = max(MStar_S, min(1.25, MStar_N))
3244 if (CS%MStar_Cap > 0.0) MStar = min( CS%MStar_Cap,MStar )
3245 elseif ( CS%mstar_scheme == MStar_from_RH18 ) then
3246 if (CS%answer_date < 20190101) then
3247 MStar_N = CS%RH18_MStar_cn1 * ( 1.0 - 1.0 / ( 1. + CS%RH18_MStar_cn2 * &
3248 exp( CS%RH18_mstar_CN3 * BLD * Abs_Coriolis / UStar) ) )
3249 else
3250 MSN_term = CS%RH18_MStar_cn2 * exp( CS%RH18_mstar_CN3 * BLD * Abs_Coriolis / UStar)
3251 MStar_N = (CS%RH18_MStar_cn1 * MSN_term) / ( 1. + MSN_term)
3252 endif
3253 MStar_S = CS%RH18_MStar_CS1 * ( max(0.0, Buoyancy_Flux)2 * BLD / &
3254 ( UStar5 * max(Abs_Coriolis,1.e-20US%T_to_s) ) )**CS%RH18_mstar_cs2
3255 #ifdef IOW
3256 if (CS%MStar_Cap > 0.0) MStar = min( CS%MStar_Cap,MStar_N + MStar_S)
3257 #else
3258 MStar = MStar_N + MStar_S
3259 #endif
3260 endif