From 315f0351bcfca8a56a2ec305fd062fa7ed12778c Mon Sep 17 00:00:00 2001
From: Quentin Rodier <quentin.rodier@meteo.fr>
Date: Tue, 15 Jun 2021 18:18:51 +0200
Subject: [PATCH] =?UTF-8?q?P.Marquet=2015/06/2021:=20correction=20on=20the?=
 =?UTF-8?q?=20exponent=20of=20final=20BL89=20mixing=20length=20according?=
 =?UTF-8?q?=20to=20Lemari=C3=A9=20et=20al.=202021?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 src/MNH/bl89.f90 | 11 ++++++++---
 1 file changed, 8 insertions(+), 3 deletions(-)

diff --git a/src/MNH/bl89.f90 b/src/MNH/bl89.f90
index 31db77267..665e76c1a 100644
--- a/src/MNH/bl89.f90
+++ b/src/MNH/bl89.f90
@@ -73,6 +73,7 @@ END MODULE MODI_BL89
 !!  Philippe 13/02/2018: use ifdef MNH_REAL to prevent problems with intrinsics on Blue Gene/Q
 !!                  01/2019 (Q. Rodier) support for RM17 mixing length
 !!                  03/2021 (JL Redelsperger) Ocean model case 
+!!                  06/2021 (P. Marquet) correction of exponent on final length according to Lemarié et al. 2021
 !-------------------------------------------------------------------------------
 !
 !*       0.    DECLARATIONS
@@ -138,7 +139,7 @@ INTEGER :: JRR        ! moist loop counter
 REAL    :: ZRVORD     ! Rv/Rd
 REAL    :: ZPOTE,ZLWORK1,ZLWORK2
 REAL    :: ZTEST,ZTEST0,ZTESTM ! test for vectorization
-REAL    :: Z2SQRT2
+REAL    :: Z2SQRT2,ZUSRBL89,ZBL89EXP
 !-------------------------------------------------------------------------------
 !
 Z2SQRT2=2.*SQRT(2.)
@@ -187,6 +188,10 @@ ELSE
 END IF
 !
 ZSQRT_TKE = SQRT(ZTKEM)
+!
+!ZBL89EXP is defined here because (and not in ini_cturb) because XCED is defined in read_exseg (depending on BL89/RM17)
+ZBL89EXP = LOG(16.)/(4.*LOG(XKARMAN)+LOG(XCED)-3.*LOG(XCMFS))
+ZUSRBL89 = 1./ZBL89EXP
 !-------------------------------------------------------------------------------
 !
 !*       2.    Virtual potential temperature on the model grid
@@ -354,8 +359,8 @@ DO JK=IKTB,IKTE
     ZLWORK1=MAX(ZLMDN(J1D,JK),1.E-10_MNHREAL)
     ZLWORK2=MAX(ZLWORK(J1D),1.E-10_MNHREAL)
     ZPOTE = ZLWORK1 / ZLWORK2
-    ZLWORK2=1.d0 + ZPOTE**(2./3.)
-    ZLM(J1D,JK) = Z2SQRT2*ZLWORK1/(ZLWORK2*SQRT(ZLWORK2))
+    ZLWORK2=1.d0 + ZPOTE**ZBL89EXP
+    ZLM(J1D,JK) = ZLWORK1*(2./ZLWORK2)**ZUSRBL89
   END DO 
 
 ZLM(:,JK)=MAX(ZLM(:,JK),XLINI)
-- 
GitLab