From 56fa7224dfb139ac56daaf0cf264f1edac603bf7 Mon Sep 17 00:00:00 2001
From: Quentin Rodier <quentin.rodier@meteo.fr>
Date: Tue, 3 May 2022 23:01:54 +0200
Subject: [PATCH] Quentin 03/05/2022: turb GPU reduce number of caracters by
 line (max 130) after expand in physical points

---
 src/common/aux/gradient_u_phy.F90             |   3 +-
 src/common/aux/gradient_v_phy.F90             |   3 +-
 src/common/turb/mode_emoist.F90               |  10 +-
 src/common/turb/mode_etheta.F90               |   6 +-
 src/common/turb/mode_prandtl.F90              | 246 +++++++++++-------
 src/common/turb/mode_tke_eps_sources.F90      |  42 +--
 src/common/turb/mode_tridiag.F90              |   3 +-
 src/common/turb/mode_tridiag_thermo.F90       |   9 +-
 src/common/turb/mode_tridiag_tke.F90          |   3 +-
 src/common/turb/mode_tridiag_wind.F90         |   3 +-
 src/common/turb/mode_turb_ver_dyn_flux.F90    |  21 +-
 src/common/turb/mode_turb_ver_sv_flux.F90     |   9 +-
 src/common/turb/mode_turb_ver_thermo_corr.F90 |  89 ++++---
 src/common/turb/mode_turb_ver_thermo_flux.F90 |  51 ++--
 src/common/turb/turb.F90                      |  56 ++--
 15 files changed, 342 insertions(+), 212 deletions(-)

diff --git a/src/common/aux/gradient_u_phy.F90 b/src/common/aux/gradient_u_phy.F90
index 7059cbec8..5db52ebda 100644
--- a/src/common/aux/gradient_u_phy.F90
+++ b/src/common/aux/gradient_u_phy.F90
@@ -85,7 +85,8 @@ CALL DZM_PHY(D,PA,PA_WORK)
 CALL MXM_PHY(D,PDZZ,PDZZ_WORK)
 !
 !$mnh__expand_array(JI=D%NIBC:D%NIEC,JJ=D%NJBC:D%NJEC,JK=1:D%NKT)
-PGZ_U_UW(D%NIBC:D%NIEC,D%NJBC:D%NJEC,1:D%NKT)= PA_WORK(D%NIBC:D%NIEC,D%NJBC:D%NJEC,1:D%NKT) / PDZZ_WORK(D%NIBC:D%NIEC,D%NJBC:D%NJEC,1:D%NKT)
+PGZ_U_UW(D%NIBC:D%NIEC,D%NJBC:D%NJEC,1:D%NKT)= PA_WORK(D%NIBC:D%NIEC,D%NJBC:D%NJEC,1:D%NKT) &
+                                               / PDZZ_WORK(D%NIBC:D%NIEC,D%NJBC:D%NJEC,1:D%NKT)
 !$mnh_end_expand_array(JI=D%NIBC:D%NIEC,JJ=D%NJBC:D%NJEC,JK=1:D%NKT)
 !
 !----------------------------------------------------------------------------
diff --git a/src/common/aux/gradient_v_phy.F90 b/src/common/aux/gradient_v_phy.F90
index 0ed5a5044..b6cc6cabb 100644
--- a/src/common/aux/gradient_v_phy.F90
+++ b/src/common/aux/gradient_v_phy.F90
@@ -85,7 +85,8 @@ CALL DZM_PHY(D,PA,PA_WORK)
 CALL MYM_PHY(D,PDZZ,PDZZ_WORK)
 !
 !$mnh__expand_array(JI=D%NIBC:D%NIEC,JJ=D%NJBC:D%NJEC,JK=1:D%NKT)
-PGZ_V_VW(D%NIBC:D%NIEC,D%NJBC:D%NJEC,1:D%NKT)= PA_WORK(D%NIBC:D%NIEC,D%NJBC:D%NJEC,1:D%NKT) / PDZZ_WORK(D%NIBC:D%NIEC,D%NJBC:D%NJEC,1:D%NKT)
+PGZ_V_VW(D%NIBC:D%NIEC,D%NJBC:D%NJEC,1:D%NKT)= PA_WORK(D%NIBC:D%NIEC,D%NJBC:D%NJEC,1:D%NKT) &
+                                               / PDZZ_WORK(D%NIBC:D%NIEC,D%NJBC:D%NJEC,1:D%NKT)
 !$mnh_end_expand_array(JI=D%NIBC:D%NIEC,JJ=D%NJBC:D%NJEC,JK=1:D%NKT)
 !----------------------------------------------------------------------------
 !
diff --git a/src/common/turb/mode_emoist.F90 b/src/common/turb/mode_emoist.F90
index 8f56f0d01..e57be179e 100644
--- a/src/common/turb/mode_emoist.F90
+++ b/src/common/turb/mode_emoist.F90
@@ -150,7 +150,7 @@ ELSE
   !   ZC is computed from line 3 to line 5
   !   Amoist* 2 * SRC is computed at line 6
   !
-    PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT) = ZDELTA * (PTHLM(IIB:IIE,IJB:IJE,1:D%NKT) + PLOCPEXNM(IIB:IIE,IJB:IJE,1:D%NKT)*(               &
+    PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT) = ZDELTA * (PTHLM(IIB:IIE,IJB:IJE,1:D%NKT) + PLOCPEXNM(IIB:IIE,IJB:IJE,1:D%NKT)*(           &
                                                     PRM(IIB:IIE,IJB:IJE,1:D%NKT,2)+PRM(IIB:IIE,IJB:IJE,1:D%NKT,4)))&
                             / (1. + ZRW(IIB:IIE,IJB:IJE,1:D%NKT))                                &
         +( PLOCPEXNM(IIB:IIE,IJB:IJE,1:D%NKT) * ZA(IIB:IIE,IJB:IJE,1:D%NKT)                                        &
@@ -176,11 +176,11 @@ ELSE
   !   ZC is computed from line 3 to line 5
   !   Amoist* 2 * SRC is computed at line 6
   !
-    PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT) = ZDELTA * (PTHLM(IIB:IIE,IJB:IJE,1:D%NKT) + PLOCPEXNM(IIB:IIE,IJB:IJE,1:D%NKT)*PRM(IIB:IIE,IJB:IJE,1:D%NKT,2))   &
-                            / (1. + ZRW(IIB:IIE,IJB:IJE,1:D%NKT))                                &
+    PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT) = ZDELTA * (PTHLM(IIB:IIE,IJB:IJE,1:D%NKT) + PLOCPEXNM(IIB:IIE,IJB:IJE,1:D%NKT)* &
+                                       PRM(IIB:IIE,IJB:IJE,1:D%NKT,2)) / (1. + ZRW(IIB:IIE,IJB:IJE,1:D%NKT))          &
         +( PLOCPEXNM(IIB:IIE,IJB:IJE,1:D%NKT) * ZA(IIB:IIE,IJB:IJE,1:D%NKT)                                        &
-               -(1.+ZDELTA) * (PTHLM(IIB:IIE,IJB:IJE,1:D%NKT) + PLOCPEXNM(IIB:IIE,IJB:IJE,1:D%NKT)*PRM(IIB:IIE,IJB:IJE,1:D%NKT,2))   &
-                            / (1. + ZRW(IIB:IIE,IJB:IJE,1:D%NKT))                                &
+               -(1.+ZDELTA) * (PTHLM(IIB:IIE,IJB:IJE,1:D%NKT) + PLOCPEXNM(IIB:IIE,IJB:IJE,1:D%NKT)* &
+               PRM(IIB:IIE,IJB:IJE,1:D%NKT,2)) / (1. + ZRW(IIB:IIE,IJB:IJE,1:D%NKT))                                &
          ) * PAMOIST(IIB:IIE,IJB:IJE,1:D%NKT) * 2. * PSRCM(IIB:IIE,IJB:IJE,1:D%NKT)
     !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
   END IF
diff --git a/src/common/turb/mode_etheta.F90 b/src/common/turb/mode_etheta.F90
index 764b8ce2d..34c2b6349 100644
--- a/src/common/turb/mode_etheta.F90
+++ b/src/common/turb/mode_etheta.F90
@@ -174,9 +174,9 @@ ELSE
   !   - Atheta * 2. * SRC is computed at line 6 
   !
     PETHETA(IIB:IIE,IJB:IJE,1:D%NKT) = ZA(IIB:IIE,IJB:IJE,1:D%NKT)                                                 &
-        +( PLOCPEXNM(IIB:IIE,IJB:IJE,1:D%NKT) * ZA(IIB:IIE,IJB:IJE,1:D%NKT)                                        &
-               -(1.+ZDELTA) * (PTHLM(IIB:IIE,IJB:IJE,1:D%NKT) + PLOCPEXNM(IIB:IIE,IJB:IJE,1:D%NKT)*PRM(IIB:IIE,IJB:IJE,1:D%NKT,2))   &
-                            / (1. + ZRW(IIB:IIE,IJB:IJE,1:D%NKT))                                &
+        +( PLOCPEXNM(IIB:IIE,IJB:IJE,1:D%NKT) * ZA(IIB:IIE,IJB:IJE,1:D%NKT) -(1.+ZDELTA) * (PTHLM(IIB:IIE,IJB:IJE,1:D%NKT) &
+        + PLOCPEXNM(IIB:IIE,IJB:IJE,1:D%NKT)*PRM(IIB:IIE,IJB:IJE,1:D%NKT,2))   &
+         / (1. + ZRW(IIB:IIE,IJB:IJE,1:D%NKT))                                 &
          ) * PATHETA(IIB:IIE,IJB:IJE,1:D%NKT) * 2. * PSRCM(IIB:IIE,IJB:IJE,1:D%NKT)
     !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
   END IF
diff --git a/src/common/turb/mode_prandtl.F90 b/src/common/turb/mode_prandtl.F90
index 6c391103f..3ff30addf 100644
--- a/src/common/turb/mode_prandtl.F90
+++ b/src/common/turb/mode_prandtl.F90
@@ -272,7 +272,8 @@ IF (.NOT. OHARAT) THEN
 !          1.3 1D Redelsperger numbers
 !
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = CST%XG / PTHVREF(IIB:IIE,IJB:IJE,1:D%NKT) * PLM(IIB:IIE,IJB:IJE,1:D%NKT) * PLEPS(IIB:IIE,IJB:IJE,1:D%NKT) / PTKEM(IIB:IIE,IJB:IJE,1:D%NKT)
+ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = CST%XG / PTHVREF(IIB:IIE,IJB:IJE,1:D%NKT) * PLM(IIB:IIE,IJB:IJE,1:D%NKT) & 
+                                  * PLEPS(IIB:IIE,IJB:IJE,1:D%NKT) / PTKEM(IIB:IIE,IJB:IJE,1:D%NKT)
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 CALL MZM_PHY(D,ZWORK1,PBLL_O_E)
 !
@@ -280,8 +281,10 @@ CALL GZ_M_W_PHY(D,PTHLM,PDZZ,ZWORK1)
 IF (KRR /= 0) THEN                ! moist case
   CALL GZ_M_W_PHY(D,PRM(:,:,:,1),PDZZ,ZWORK2)
   !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-  PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)= CSTURB%XCTV*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT) * PETHETA(IIB:IIE,IJB:IJE,1:D%NKT) * ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)
-  PREDR1(IIB:IIE,IJB:IJE,1:D%NKT) = CSTURB%XCTV*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT) * PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT) * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
+  PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)= CSTURB%XCTV*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT) * PETHETA(IIB:IIE,IJB:IJE,1:D%NKT) &
+                                    * ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)
+  PREDR1(IIB:IIE,IJB:IJE,1:D%NKT) = CSTURB%XCTV*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT) * PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT) &
+                                    * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
   !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 ELSE                              ! dry case
   !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
@@ -708,7 +711,8 @@ IF (HTURBDIM=='3DIM') THEN
     ZW2(IIB:IIE,IJB:IJE,1:D%NKT) = 0.5* PRED2TH3(IIB:IIE,IJB:IJE,1:D%NKT)
 
     PHI3(IIB:IIE,IJB:IJE,1:D%NKT)= 1. -                                       &
-            (PRED2TH3(IIB:IIE,IJB:IJE,1:D%NKT) / PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT) + ZW2(IIB:IIE,IJB:IJE,1:D%NKT)) / ZW1(IIB:IIE,IJB:IJE,1:D%NKT)
+            (PRED2TH3(IIB:IIE,IJB:IJE,1:D%NKT) / PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT) + ZW2(IIB:IIE,IJB:IJE,1:D%NKT)) &
+            / ZW1(IIB:IIE,IJB:IJE,1:D%NKT)
   END IF
   !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
 
@@ -762,15 +766,19 @@ IJB=D%NJBC
 DO JSV=1,KSV
   !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
   PSI_SV(IIB:IIE,IJB:IJE,1:D%NKT,JSV) = ( 1.                                             &
-    - (CSTURB%XCPR3+CSTURB%XCPR5) * (PRED2THS(IIB:IIE,IJB:IJE,1:D%NKT,JSV)/PREDS1(IIB:IIE,IJB:IJE,1:D%NKT,JSV)-PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)) &
-    - (CSTURB%XCPR4+CSTURB%XCPR5) * (PRED2RS(IIB:IIE,IJB:IJE,1:D%NKT,JSV)/PREDS1(IIB:IIE,IJB:IJE,1:D%NKT,JSV)-PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)) &
-    - CSTURB%XCPR3 * PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT) * PPHI3(IIB:IIE,IJB:IJE,1:D%NKT) - CSTURB%XCPR4 * PREDR1(IIB:IIE,IJB:IJE,1:D%NKT) * PPSI3(IIB:IIE,IJB:IJE,1:D%NKT)                 &
+    - (CSTURB%XCPR3+CSTURB%XCPR5) * &
+    (PRED2THS(IIB:IIE,IJB:IJE,1:D%NKT,JSV)/PREDS1(IIB:IIE,IJB:IJE,1:D%NKT,JSV)-PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)) &
+    - (CSTURB%XCPR4+CSTURB%XCPR5) * &
+    (PRED2RS(IIB:IIE,IJB:IJE,1:D%NKT,JSV)/PREDS1(IIB:IIE,IJB:IJE,1:D%NKT,JSV)-PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)) &
+    - CSTURB%XCPR3 * &
+    PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT) * PPHI3(IIB:IIE,IJB:IJE,1:D%NKT) &
+    - CSTURB%XCPR4 * PREDR1(IIB:IIE,IJB:IJE,1:D%NKT) * PPSI3(IIB:IIE,IJB:IJE,1:D%NKT)                 &
                 ) / (  1. + CSTURB%XCPR5 * ( PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT) + PREDR1(IIB:IIE,IJB:IJE,1:D%NKT) ) )           
   
   !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
 !        control of the PSI_SV positivity
   !$mnh_expand_where(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
-  WHERE ( (PSI_SV(IIB:IIE,IJB:IJE,1:D%NKT,JSV) <=0.).AND. (PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)) <= 0. )
+  WHERE ( (PSI_SV(IIB:IIE,IJB:IJE,1:D%NKT,JSV) <=0.).AND. (PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))<=0.)
     PSI_SV(IIB:IIE,IJB:IJE,1:D%NKT,JSV)=CSTURB%XPHI_LIM
   END WHERE
   !$mnh_end_expand_where(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
@@ -818,13 +826,16 @@ IF (HTURBDIM=='3DIM') THEN
     WHERE (PPHI3(IIB:IIE,IJB:IJE,1:D%NKT)<=CSTURB%XPHI_LIM)
 #endif
       D_PHI3DTDZ_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = PPHI3(IIB:IIE,IJB:IJE,1:D%NKT)                       &
-          * (1. - PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT) * (3./2.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))          &
-               /((1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))*(1.+1./2.*(PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))))) &
-          + (1.+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))*(PRED2THR3(IIB:IIE,IJB:IJE,1:D%NKT)+PRED2TH3(IIB:IIE,IJB:IJE,1:D%NKT))                       &
+          * (1. - PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT) * (3./2.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))   &
+               /((1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)) &
+               *(1.+1./2.*(PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))))) &
+          + (1.+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))*(PRED2THR3(IIB:IIE,IJB:IJE,1:D%NKT)+PRED2TH3(IIB:IIE,IJB:IJE,1:D%NKT))         &
                / (PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)*(1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))* &
                  (1.+1./2.*(PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)))) &
-          - (1./2.*PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT) * (1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)))           &
-               / ((1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))*(1.+1./2.*(PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))))
+          - (1./2.*PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT) &
+          * (1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)))           &
+               / ((1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))&
+               *(1.+1./2.*(PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))))
     ELSEWHERE
       D_PHI3DTDZ_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = PPHI3(IIB:IIE,IJB:IJE,1:D%NKT)
     ENDWHERE
@@ -841,8 +852,10 @@ IF (HTURBDIM=='3DIM') THEN
     D_PHI3DTDZ_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = PPHI3(IIB:IIE,IJB:IJE,1:D%NKT)             &
           * (1. - PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT) * (3./2.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT))      &
                /((1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT))*(1.+1./2.*PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT))))        &
-          + PRED2TH3(IIB:IIE,IJB:IJE,1:D%NKT) / (PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)*(1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT))*(1.+1./2.*PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT))) &
-          - 1./2.*PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT) / ((1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT))*(1.+1./2.*PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)))
+          + PRED2TH3(IIB:IIE,IJB:IJE,1:D%NKT) &
+          / (PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)*(1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT))*(1.+1./2.*PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT))) &
+          - 1./2.*PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT) &
+          / ((1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT))*(1.+1./2.*PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)))
     ELSEWHERE
       D_PHI3DTDZ_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = PPHI3(IIB:IIE,IJB:IJE,1:D%NKT)
     ENDWHERE
@@ -909,13 +922,17 @@ IF (HTURBDIM=='3DIM') THEN
 #else
     WHERE (PPHI3(IIB:IIE,IJB:IJE,1:D%NKT)<=CSTURB%XPHI_LIM)
 #endif
-      D_PHI3DRDZ_O_DDRDZ(IIB:IIE,IJB:IJE,1:D%NKT) =          &
-                PPHI3(IIB:IIE,IJB:IJE,1:D%NKT) * (1.-PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)*(3./2.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)) &
-                  / ((1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))*(1.+1./2.*(PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)))))  &
-              - PREDR1(IIB:IIE,IJB:IJE,1:D%NKT) * (PRED2THR3(IIB:IIE,IJB:IJE,1:D%NKT)+PRED2TH3(IIB:IIE,IJB:IJE,1:D%NKT)) / (PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)         &
-                  * (1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))*(1.+1./2.*(PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))))    &
-              + PREDR1(IIB:IIE,IJB:IJE,1:D%NKT) * (1./2.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))                  &
-                  / ((1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))*(1.+1./2.*(PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))))
+      D_PHI3DRDZ_O_DDRDZ(IIB:IIE,IJB:IJE,1:D%NKT) = PPHI3(IIB:IIE,IJB:IJE,1:D%NKT) &
+      * (1.-PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)*(3./2.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)) &
+      / ((1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)) & 
+      *(1.+1./2.*(PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)))))  &
+      - PREDR1(IIB:IIE,IJB:IJE,1:D%NKT) &
+      * (PRED2THR3(IIB:IIE,IJB:IJE,1:D%NKT)+PRED2TH3(IIB:IIE,IJB:IJE,1:D%NKT)) / (PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)         &
+      * (1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))*&
+      (1.+1./2.*(PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))))    &
+      + PREDR1(IIB:IIE,IJB:IJE,1:D%NKT) * (1./2.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))         &
+      / ((1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))&
+      *(1.+1./2.*(PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))))
     ELSEWHERE
       D_PHI3DRDZ_O_DDRDZ(IIB:IIE,IJB:IJE,1:D%NKT) = PPHI3(IIB:IIE,IJB:IJE,1:D%NKT)
     END WHERE
@@ -980,7 +997,8 @@ IF (HTURBDIM=='3DIM') THEN
    ! by derivation of (phi3 dtdz) * dtdz according to dtdz we obtain:
    ZWORK1 = D_PHI3DTDZ_O_DDTDZ(D,CSTURB,PPHI3,PREDTH1,PREDR1,PRED2TH3,PRED2THR3,HTURBDIM,OUSERV) 
    !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-   D_PHI3DTDZ2_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = PDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) * (PPHI3(IIB:IIE,IJB:IJE,1:D%NKT) +  ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT))
+   D_PHI3DTDZ2_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = PDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) &
+   * (PPHI3(IIB:IIE,IJB:IJE,1:D%NKT) +  ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT))
    !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 ELSE
         !* 1DIM case
@@ -1032,7 +1050,8 @@ IJE=D%NJEC
 IJB=D%NJBC
 
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-M3_WTH_WTH2(IIB:IIE,IJB:IJE,1:D%NKT) = CSTURB%XCSHF*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*PETHETA(IIB:IIE,IJB:IJE,1:D%NKT)*0.5/CSTURB%XCTD        &
+M3_WTH_WTH2(IIB:IIE,IJB:IJE,1:D%NKT) = CSTURB%XCSHF*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)&
+                   * PETHETA(IIB:IIE,IJB:IJE,1:D%NKT)*0.5/CSTURB%XCTD        &
                    * (1.+0.5*PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)) / PD(IIB:IIE,IJB:IJE,1:D%NKT)
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 M3_WTH_WTH2(IIB:IIE,IJB:IJE,IKB-1)=M3_WTH_WTH2(IIB:IIE,IJB:IJE,IKB)
@@ -1063,9 +1082,11 @@ IJE=D%NJEC
 IJB=D%NJBC
 
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-D_M3_WTH_WTH2_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = (  0.5*CSTURB%XCSHF*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*PETHETA(IIB:IIE,IJB:IJE,1:D%NKT)*0.5/CSTURB%XCTD/PD(IIB:IIE,IJB:IJE,1:D%NKT) &
-                                - PM3_WTH_WTH2(IIB:IIE,IJB:IJE,1:D%NKT)/PD(IIB:IIE,IJB:IJE,1:D%NKT)*(1.5+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))  )&
-                             * PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT) * PETHETA(IIB:IIE,IJB:IJE,1:D%NKT) * CSTURB%XCTV
+D_M3_WTH_WTH2_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = &
+(0.5*CSTURB%XCSHF*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*PETHETA(IIB:IIE,IJB:IJE,1:D%NKT)*0.5/CSTURB%XCTD/PD(IIB:IIE,IJB:IJE,1:D%NKT) &
+- PM3_WTH_WTH2(IIB:IIE,IJB:IJE,1:D%NKT)/PD(IIB:IIE,IJB:IJE,1:D%NKT)&
+*(1.5+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))  )&
+* PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT) * PETHETA(IIB:IIE,IJB:IJE,1:D%NKT) * CSTURB%XCTV
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 !
 D_M3_WTH_WTH2_O_DDTDZ(IIB:IIE,IJB:IJE,IKB-1)=D_M3_WTH_WTH2_O_DDTDZ(IIB:IIE,IJB:IJE,IKB)
@@ -1096,8 +1117,9 @@ IJE=D%NJEC
 IJB=D%NJBC
 CALL MZM_PHY(D,PTKE,ZWORK1)
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-M3_WTH_W2TH(IIB:IIE,IJB:IJE,1:D%NKT) = CSTURB%XCSHF*PKEFF(IIB:IIE,IJB:IJE,1:D%NKT)*1.5/ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)              &
-  * (1. - 0.5*PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)*(1.+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))/PD(IIB:IIE,IJB:IJE,1:D%NKT) ) / (1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT))
+M3_WTH_W2TH(IIB:IIE,IJB:IJE,1:D%NKT) = CSTURB%XCSHF*PKEFF(IIB:IIE,IJB:IJE,1:D%NKT)*1.5/ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) &
+  * (1. - 0.5*PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)*(1.+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))/PD(IIB:IIE,IJB:IJE,1:D%NKT) ) &
+  / (1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT))
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 !
 M3_WTH_W2TH(IIB:IIE,IJB:IJE,IKB-1)=M3_WTH_W2TH(IIB:IIE,IJB:IJE,IKB)
@@ -1135,7 +1157,8 @@ D_M3_WTH_W2TH_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = &
  - CSTURB%XCSHF*PKEFF(IIB:IIE,IJB:IJE,1:D%NKT)*1.5/ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)/(1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT))**2 &
  * CSTURB%XCTV*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*PETHETA(IIB:IIE,IJB:IJE,1:D%NKT)  &
  * (1. - 0.5*PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)*(1.+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))/PD(IIB:IIE,IJB:IJE,1:D%NKT)* &
-   ( 1.+(1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT))*(1.5+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT))/PD(IIB:IIE,IJB:IJE,1:D%NKT)) )
+   ( 1.+(1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT))*(1.5+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT))&
+   /PD(IIB:IIE,IJB:IJE,1:D%NKT)) )
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 !
 D_M3_WTH_W2TH_O_DDTDZ(IIB:IIE,IJB:IJE,IKB-1)=D_M3_WTH_W2TH_O_DDTDZ(IIB:IIE,IJB:IJE,IKB)
@@ -1168,8 +1191,9 @@ IJB=D%NJBC
 
 CALL MZM_PHY(D,PTKE,ZWORK1)
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-M3_WTH_W2R(IIB:IIE,IJB:IJE,1:D%NKT) = - CSTURB%XCSHF*PKEFF(IIB:IIE,IJB:IJE,1:D%NKT)*0.75*CSTURB%XCTV*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT) &
-                    /ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)*PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT)*PDTDZ(IIB:IIE,IJB:IJE,1:D%NKT)/PD(IIB:IIE,IJB:IJE,1:D%NKT)
+M3_WTH_W2R(IIB:IIE,IJB:IJE,1:D%NKT) = &
+  - CSTURB%XCSHF*PKEFF(IIB:IIE,IJB:IJE,1:D%NKT)*0.75*CSTURB%XCTV*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT) &
+  /ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)*PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT)*PDTDZ(IIB:IIE,IJB:IJE,1:D%NKT)/PD(IIB:IIE,IJB:IJE,1:D%NKT)
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 !
 M3_WTH_W2R(IIB:IIE,IJB:IJE,IKB-1)=M3_WTH_W2R(IIB:IIE,IJB:IJE,IKB)
@@ -1203,9 +1227,11 @@ IJB=D%NJBC
 
 CALL MZM_PHY(D,PTKE,ZWORK1)
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-D_M3_WTH_W2R_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = - CSTURB%XCSHF*PKEFF(IIB:IIE,IJB:IJE,1:D%NKT)*0.75*CSTURB%XCTV*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT) &
+D_M3_WTH_W2R_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = &
+- CSTURB%XCSHF*PKEFF(IIB:IIE,IJB:IJE,1:D%NKT)*0.75*CSTURB%XCTV*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT) &
                                /ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)*PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT)/PD(IIB:IIE,IJB:IJE,1:D%NKT) &
-                                     * (1. -  PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)*(1.5+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))/PD(IIB:IIE,IJB:IJE,1:D%NKT))
+                                     * (1. -  PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)*(1.5+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)& 
+                                     +PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))/PD(IIB:IIE,IJB:IJE,1:D%NKT))
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 !
 D_M3_WTH_W2R_O_DDTDZ(IIB:IIE,IJB:IJE,IKB-1)=D_M3_WTH_W2R_O_DDTDZ(IIB:IIE,IJB:IJE,IKB)
@@ -1239,11 +1265,13 @@ IIB=D%NIBC
 IJE=D%NJEC
 IJB=D%NJBC
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = PBETA(IIB:IIE,IJB:IJE,1:D%NKT)*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/(PSQRT_TKE(IIB:IIE,IJB:IJE,1:D%NKT)*PTKE(IIB:IIE,IJB:IJE,1:D%NKT))
+ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = PBETA(IIB:IIE,IJB:IJE,1:D%NKT)*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT) &
+                                 /(PSQRT_TKE(IIB:IIE,IJB:IJE,1:D%NKT)*PTKE(IIB:IIE,IJB:IJE,1:D%NKT))
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 CALL MZM_PHY(D,ZWORK1,ZWORK2)
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-M3_WTH_WR2(IIB:IIE,IJB:IJE,1:D%NKT) = - CSTURB%XCSHF*PKEFF(IIB:IIE,IJB:IJE,1:D%NKT)*0.25*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*CSTURB%XCTV*PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT)**2       &
+M3_WTH_WR2(IIB:IIE,IJB:IJE,1:D%NKT) = - CSTURB%XCSHF*PKEFF(IIB:IIE,IJB:IJE,1:D%NKT)& 
+                           *0.25*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*CSTURB%XCTV*PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT)**2 &
                            *ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)/CSTURB%XCTD*PDTDZ(IIB:IIE,IJB:IJE,1:D%NKT)/PD(IIB:IIE,IJB:IJE,1:D%NKT)
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 !
@@ -1279,13 +1307,16 @@ IIB=D%NIBC
 IJE=D%NJEC
 IJB=D%NJBC
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = PBETA(IIB:IIE,IJB:IJE,1:D%NKT)*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/(PSQRT_TKE(IIB:IIE,IJB:IJE,1:D%NKT)*PTKE(IIB:IIE,IJB:IJE,1:D%NKT))
+ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = PBETA(IIB:IIE,IJB:IJE,1:D%NKT)*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)&
+                                  /(PSQRT_TKE(IIB:IIE,IJB:IJE,1:D%NKT)*PTKE(IIB:IIE,IJB:IJE,1:D%NKT))
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 CALL MZM_PHY(D,ZWORK1,ZWORK2)
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-D_M3_WTH_WR2_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = - CSTURB%XCSHF*PKEFF(IIB:IIE,IJB:IJE,1:D%NKT)*0.25*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*CSTURB%XCTV*PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT)**2 &
+D_M3_WTH_WR2_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = - CSTURB%XCSHF*PKEFF(IIB:IIE,IJB:IJE,1:D%NKT)& 
+                           *0.25*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*CSTURB%XCTV*PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT)**2 &
                            *ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)/CSTURB%XCTD/PD(IIB:IIE,IJB:IJE,1:D%NKT)     &
-                           * (1. -  PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)*(1.5+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))/PD(IIB:IIE,IJB:IJE,1:D%NKT))
+                           * (1. -  PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)* &
+                           (1.5+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))/PD(IIB:IIE,IJB:IJE,1:D%NKT))
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 !
 D_M3_WTH_WR2_O_DDTDZ(IIB:IIE,IJB:IJE,IKB-1)=D_M3_WTH_WR2_O_DDTDZ(IIB:IIE,IJB:IJE,IKB)
@@ -1319,12 +1350,14 @@ IJE=D%NJEC
 IJB=D%NJBC
 
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = PBETA(IIB:IIE,IJB:IJE,1:D%NKT)*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/(PSQRT_TKE(IIB:IIE,IJB:IJE,1:D%NKT)*PTKE(IIB:IIE,IJB:IJE,1:D%NKT))
+ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = PBETA(IIB:IIE,IJB:IJE,1:D%NKT)*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)&
+                                  /(PSQRT_TKE(IIB:IIE,IJB:IJE,1:D%NKT)*PTKE(IIB:IIE,IJB:IJE,1:D%NKT))
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 CALL MZM_PHY(D,ZWORK1,ZWORK2)
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-M3_WTH_WTHR(IIB:IIE,IJB:IJE,1:D%NKT) = CSTURB%XCSHF*PKEFF(IIB:IIE,IJB:IJE,1:D%NKT)*PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT)*ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT) &
-                         *0.5*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/CSTURB%XCTD*(1+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))/PD(IIB:IIE,IJB:IJE,1:D%NKT)
+M3_WTH_WTHR(IIB:IIE,IJB:IJE,1:D%NKT) = &
+                   CSTURB%XCSHF*PKEFF(IIB:IIE,IJB:IJE,1:D%NKT)*PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT)*ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT) &
+                   *0.5*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/CSTURB%XCTD*(1+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))/PD(IIB:IIE,IJB:IJE,1:D%NKT)
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 !
 M3_WTH_WTHR(IIB:IIE,IJB:IJE,IKB-1)=M3_WTH_WTHR(IIB:IIE,IJB:IJE,IKB)
@@ -1355,8 +1388,9 @@ IJE=D%NJEC
 IJB=D%NJBC
 
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-D_M3_WTH_WTHR_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = - PM3_WTH_WTHR(IIB:IIE,IJB:IJE,1:D%NKT) * (1.5+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))&
-                               /PD(IIB:IIE,IJB:IJE,1:D%NKT)*CSTURB%XCTV*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*PETHETA(IIB:IIE,IJB:IJE,1:D%NKT)
+D_M3_WTH_WTHR_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = &
+                - PM3_WTH_WTHR(IIB:IIE,IJB:IJE,1:D%NKT) * (1.5+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))&
+                /PD(IIB:IIE,IJB:IJE,1:D%NKT)*CSTURB%XCTV*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*PETHETA(IIB:IIE,IJB:IJE,1:D%NKT)
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 !
 D_M3_WTH_WTHR_O_DDTDZ(IIB:IIE,IJB:IJE,IKB-1)=D_M3_WTH_WTHR_O_DDTDZ(IIB:IIE,IJB:IJE,IKB)
@@ -1388,7 +1422,8 @@ IIB=D%NIBC
 IJE=D%NJEC
 IJB=D%NJBC
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = (1.-0.5*PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)*(1.+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))/PD(IIB:IIE,IJB:IJE,1:D%NKT))/(1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT))*PDTDZ(IIB:IIE,IJB:IJE,1:D%NKT)
+ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = (1.-0.5*PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)*(1.+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))& 
+                                /PD(IIB:IIE,IJB:IJE,1:D%NKT))/(1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT))*PDTDZ(IIB:IIE,IJB:IJE,1:D%NKT)
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 CALL MZF_PHY(D,ZWORK1,ZWORK2)
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
@@ -1430,13 +1465,15 @@ IF (OUSERV) THEN
 !          (1.-0.5*PREDR1*(1.+PREDR1)/PD)*(1.-(1.5+PREDTH1+PREDR1)*(1.+PREDTH1)/PD )  &
 !        / (1.+PREDTH1)**2, D%NKA, D%NKU, D%NKL)
   !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-  ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = (1.-0.5*PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)*(1.+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))/PD(IIB:IIE,IJB:IJE,1:D%NKT))*(1.-(1.5+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))*   &
-             PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)*(1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT))/PD(IIB:IIE,IJB:IJE,1:D%NKT) ) / (1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT))**2
+  ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = (1.-0.5*PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)*(1.+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))&
+             / PD(IIB:IIE,IJB:IJE,1:D%NKT))*(1.-(1.5+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))   &
+             * PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)*(1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT))/PD(IIB:IIE,IJB:IJE,1:D%NKT) ) &
+             / (1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT))**2
   !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
   CALL MZF_PHY(D,ZWORK1,ZWORK2)
   !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-  D_M3_TH2_W2TH_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = - 1.5*PLM(IIB:IIE,IJB:IJE,1:D%NKT)*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/PTKE(IIB:IIE,IJB:IJE,1:D%NKT)*CSTURB%XCTV * &
-                                 ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
+  D_M3_TH2_W2TH_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = - 1.5*PLM(IIB:IIE,IJB:IJE,1:D%NKT)*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT) &
+                                                   /PTKE(IIB:IIE,IJB:IJE,1:D%NKT)*CSTURB%XCTV * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
   !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 ELSE
   !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
@@ -1444,8 +1481,8 @@ ELSE
   !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
   CALL MZF_PHY(D,ZWORK1,ZWORK2)
   !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-  D_M3_TH2_W2TH_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = - 1.5*PLM(IIB:IIE,IJB:IJE,1:D%NKT)*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/PTKE(IIB:IIE,IJB:IJE,1:D%NKT)*CSTURB%XCTV &
-                                 * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
+  D_M3_TH2_W2TH_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = - 1.5*PLM(IIB:IIE,IJB:IJE,1:D%NKT)*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT) & 
+                                                   /PTKE(IIB:IIE,IJB:IJE,1:D%NKT)*CSTURB%XCTV * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
   !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 END IF
 !
@@ -1476,11 +1513,12 @@ IIB=D%NIBC
 IJE=D%NJEC
 IJB=D%NJBC
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = (1.+0.5*PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+1.5*PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)+0.5*PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)**2)/PD(IIB:IIE,IJB:IJE,1:D%NKT)
+ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = (1.+0.5*PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT) &
+                         +1.5*PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)+0.5*PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)**2)/PD(IIB:IIE,IJB:IJE,1:D%NKT)
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 CALL MZF_PHY(D,ZWORK1,ZWORK2)
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-M3_TH2_WTH2(IIB:IIE,IJB:IJE,1:D%NKT) = PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)*0.5/CSTURB%XCTD/PSQRT_TKE(IIB:IIE,IJB:IJE,1:D%NKT)          &
+M3_TH2_WTH2(IIB:IIE,IJB:IJE,1:D%NKT) = PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)*0.5/CSTURB%XCTD/PSQRT_TKE(IIB:IIE,IJB:IJE,1:D%NKT) &
                      * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 !
@@ -1513,13 +1551,15 @@ IIB=D%NIBC
 IJE=D%NJEC
 IJB=D%NJBC
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*PETHETA(IIB:IIE,IJB:IJE,1:D%NKT)* (0.5/PD(IIB:IIE,IJB:IJE,1:D%NKT)                                                   &
-             - (1.5+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))*(1.+0.5*PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+1.5*PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)+0.5*PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)**2)/PD(IIB:IIE,IJB:IJE,1:D%NKT)**2)
+ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*PETHETA(IIB:IIE,IJB:IJE,1:D%NKT) &
+             * (0.5/PD(IIB:IIE,IJB:IJE,1:D%NKT) - (1.5+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))& 
+             *(1.+0.5*PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+1.5*PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)& 
+             +0.5*PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)**2)/PD(IIB:IIE,IJB:IJE,1:D%NKT)**2)
  !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 CALL MZF_PHY(D,ZWORK1,ZWORK2)
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-D_M3_TH2_WTH2_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)*0.5/CSTURB%XCTD/PSQRT_TKE(IIB:IIE,IJB:IJE,1:D%NKT)*CSTURB%XCTV &
-                               * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
+D_M3_TH2_WTH2_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = PLEPS(IIB:IIE,IJB:IJE,1:D%NKT) & 
+                                 *0.5/CSTURB%XCTD/PSQRT_TKE(IIB:IIE,IJB:IJE,1:D%NKT)*CSTURB%XCTV * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 !
 D_M3_TH2_WTH2_O_DDTDZ(IIB:IIE,IJB:IJE,IKB-1)=D_M3_TH2_WTH2_O_DDTDZ(IIB:IIE,IJB:IJE,IKB)
@@ -1551,7 +1591,8 @@ IIB=D%NIBC
 IJE=D%NJEC
 IJB=D%NJBC
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT)/PD(IIB:IIE,IJB:IJE,1:D%NKT)*PDTDZ(IIB:IIE,IJB:IJE,1:D%NKT)**2
+ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT) & 
+                                 /PD(IIB:IIE,IJB:IJE,1:D%NKT)*PDTDZ(IIB:IIE,IJB:IJE,1:D%NKT)**2
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 CALL MZF_PHY(D,ZWORK1,ZWORK2)
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
@@ -1590,12 +1631,14 @@ IIB=D%NIBC
 IJE=D%NJEC
 IJB=D%NJBC
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) =  PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT)/PD(IIB:IIE,IJB:IJE,1:D%NKT)*PDTDZ(IIB:IIE,IJB:IJE,1:D%NKT)*(2.-PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)*(1.5+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))/PD(IIB:IIE,IJB:IJE,1:D%NKT))
+ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) =  PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT)& 
+ /PD(IIB:IIE,IJB:IJE,1:D%NKT)*PDTDZ(IIB:IIE,IJB:IJE,1:D%NKT)*(2.-PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)* & 
+ (1.5+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))/PD(IIB:IIE,IJB:IJE,1:D%NKT))
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 CALL MZF_PHY(D,ZWORK1,ZWORK2)
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-D_M3_TH2_W2R_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = 0.75*CSTURB%XCTV**2*PLM(IIB:IIE,IJB:IJE,1:D%NKT)*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/PTKE(IIB:IIE,IJB:IJE,1:D%NKT) &
-                              * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
+D_M3_TH2_W2R_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = 0.75*CSTURB%XCTV**2*PLM(IIB:IIE,IJB:IJE,1:D%NKT) *PLEPS(IIB:IIE,IJB:IJE,1:D%NKT) &
+                                                /PTKE(IIB:IIE,IJB:IJE,1:D%NKT) * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 !
 D_M3_TH2_W2R_O_DDTDZ(IIB:IIE,IJB:IJE,IKB-1)=D_M3_TH2_W2R_O_DDTDZ(IIB:IIE,IJB:IJE,IKB)
@@ -1626,7 +1669,8 @@ IIB=D%NIBC
 IJE=D%NJEC
 IJB=D%NJBC
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = (PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT)*PDTDZ(IIB:IIE,IJB:IJE,1:D%NKT))**2/PD(IIB:IIE,IJB:IJE,1:D%NKT)
+ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = (PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT)& 
+                                  *PDTDZ(IIB:IIE,IJB:IJE,1:D%NKT))**2/PD(IIB:IIE,IJB:IJE,1:D%NKT)
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 CALL MZF_PHY(D,ZWORK1,ZWORK2)
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
@@ -1664,12 +1708,14 @@ IIB=D%NIBC
 IJE=D%NJEC
 IJB=D%NJBC
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = (PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT))**2*PDTDZ(IIB:IIE,IJB:IJE,1:D%NKT)/PD(IIB:IIE,IJB:IJE,1:D%NKT)*(2.-PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)*(1.5+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))/PD(IIB:IIE,IJB:IJE,1:D%NKT))
+ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = (PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT))**2 & 
+*PDTDZ(IIB:IIE,IJB:IJE,1:D%NKT)/PD(IIB:IIE,IJB:IJE,1:D%NKT)*(2.-PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT) & 
+*(1.5+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))/PD(IIB:IIE,IJB:IJE,1:D%NKT))
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 CALL MZF_PHY(D,ZWORK1,ZWORK2)
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-D_M3_TH2_WR2_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = 0.25*CSTURB%XCTV**2*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/PSQRT_TKE(IIB:IIE,IJB:IJE,1:D%NKT)/CSTURB%XCTD &
-                              *  ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
+D_M3_TH2_WR2_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = 0.25*CSTURB%XCTV**2*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT) & 
+                                               / PSQRT_TKE(IIB:IIE,IJB:IJE,1:D%NKT)/CSTURB%XCTD * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 !
 D_M3_TH2_WR2_O_DDTDZ(IIB:IIE,IJB:IJE,IKB-1)=D_M3_TH2_WR2_O_DDTDZ(IIB:IIE,IJB:IJE,IKB)
@@ -1701,12 +1747,13 @@ IIB=D%NIBC
 IJE=D%NJEC
 IJB=D%NJBC
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT)*PDTDZ(IIB:IIE,IJB:IJE,1:D%NKT)*(1.+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))/PD(IIB:IIE,IJB:IJE,1:D%NKT)
+ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT) & 
+                                * PDTDZ(IIB:IIE,IJB:IJE,1:D%NKT)*(1.+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))/PD(IIB:IIE,IJB:IJE,1:D%NKT)
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 CALL MZF_PHY(D,ZWORK1,ZWORK2)
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-M3_TH2_WTHR(IIB:IIE,IJB:IJE,1:D%NKT) = - 0.5*CSTURB%XCTV*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/PSQRT_TKE(IIB:IIE,IJB:IJE,1:D%NKT)/CSTURB%XCTD &
-                       *ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
+M3_TH2_WTHR(IIB:IIE,IJB:IJE,1:D%NKT) = - 0.5*CSTURB%XCTV*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT) & 
+                                       / PSQRT_TKE(IIB:IIE,IJB:IJE,1:D%NKT)/CSTURB%XCTD * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 !
 M3_TH2_WTHR(IIB:IIE,IJB:IJE,IKB-1)=M3_TH2_WTHR(IIB:IIE,IJB:IJE,IKB)
@@ -1739,12 +1786,14 @@ IIB=D%NIBC
 IJE=D%NJEC
 IJB=D%NJBC
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT)*(1.+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))/PD(IIB:IIE,IJB:IJE,1:D%NKT) * (1. -PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)*(1.5+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))/PD(IIB:IIE,IJB:IJE,1:D%NKT))
+ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*PEMOIST(IIB:IIE,IJB:IJE,1:D%NKT)* & 
+                 (1.+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))/PD(IIB:IIE,IJB:IJE,1:D%NKT) * (1. -PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)*& 
+                 (1.5+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))/PD(IIB:IIE,IJB:IJE,1:D%NKT))
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 CALL MZF_PHY(D,ZWORK1,ZWORK2)
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-D_M3_TH2_WTHR_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = - 0.5*CSTURB%XCTV*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/PSQRT_TKE(IIB:IIE,IJB:IJE,1:D%NKT)/CSTURB%XCTD &
-                               * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT) 
+D_M3_TH2_WTHR_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = - 0.5*CSTURB%XCTV*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT) & 
+                                                / PSQRT_TKE(IIB:IIE,IJB:IJE,1:D%NKT)/CSTURB%XCTD * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT) 
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 !
 D_M3_TH2_WTHR_O_DDTDZ(IIB:IIE,IJB:IJE,IKB-1)=D_M3_TH2_WTHR_O_DDTDZ(IIB:IIE,IJB:IJE,IKB)
@@ -1774,7 +1823,8 @@ IIB=D%NIBC
 IJE=D%NJEC
 IJB=D%NJBC
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) =  (1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT))*(1.+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))/PD(IIB:IIE,IJB:IJE,1:D%NKT)
+ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) =  (1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT))* & 
+                                   (1.+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))/PD(IIB:IIE,IJB:IJE,1:D%NKT)
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 CALL MZF_PHY(D,ZWORK1,ZWORK2)
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
@@ -1811,12 +1861,14 @@ IIB=D%NIBC
 IJE=D%NJEC
 IJB=D%NJBC
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = PETHETA(IIB:IIE,IJB:IJE,1:D%NKT)*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)/PD(IIB:IIE,IJB:IJE,1:D%NKT)*(1.+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))*(1.-(1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT))*(1.5+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))/PD(IIB:IIE,IJB:IJE,1:D%NKT))
+ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = PETHETA(IIB:IIE,IJB:IJE,1:D%NKT)*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)/PD(IIB:IIE,IJB:IJE,1:D%NKT) & 
+                             *(1.+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))*(1.-(1.+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)) & 
+                             *(1.5+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))/PD(IIB:IIE,IJB:IJE,1:D%NKT))
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 CALL MZF_PHY(D,ZWORK1,ZWORK2)
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-D_M3_THR_WTHR_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = 0.5*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/PSQRT_TKE(IIB:IIE,IJB:IJE,1:D%NKT)/CSTURB%XCTD * CSTURB%XCTV &
-                               * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
+D_M3_THR_WTHR_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = 0.5*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/PSQRT_TKE(IIB:IIE,IJB:IJE,1:D%NKT) & 
+                                                / CSTURB%XCTD * CSTURB%XCTV * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 !
 D_M3_THR_WTHR_O_DDTDZ(IIB:IIE,IJB:IJE,IKB-1)=D_M3_THR_WTHR_O_DDTDZ(IIB:IIE,IJB:IJE,IKB)
@@ -1848,12 +1900,13 @@ IIB=D%NIBC
 IJE=D%NJEC
 IJB=D%NJBC
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = (1.+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*PETHETA(IIB:IIE,IJB:IJE,1:D%NKT)*PDRDZ(IIB:IIE,IJB:IJE,1:D%NKT)/PD(IIB:IIE,IJB:IJE,1:D%NKT)
+ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = (1.+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)* & 
+                                  PETHETA(IIB:IIE,IJB:IJE,1:D%NKT)*PDRDZ(IIB:IIE,IJB:IJE,1:D%NKT)/PD(IIB:IIE,IJB:IJE,1:D%NKT)
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 CALL MZF_PHY(D,ZWORK1,ZWORK2)
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-M3_THR_WTH2(IIB:IIE,IJB:IJE,1:D%NKT) = - 0.25*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/PSQRT_TKE(IIB:IIE,IJB:IJE,1:D%NKT)/CSTURB%XCTD*CSTURB%XCTV &
-                     * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
+M3_THR_WTH2(IIB:IIE,IJB:IJE,1:D%NKT) = - 0.25*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT) & 
+                                    / PSQRT_TKE(IIB:IIE,IJB:IJE,1:D%NKT)/CSTURB%XCTD*CSTURB%XCTV * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 !
 M3_THR_WTH2(IIB:IIE,IJB:IJE,IKB-1)=M3_THR_WTH2(IIB:IIE,IJB:IJE,IKB)
@@ -1886,12 +1939,15 @@ IIB=D%NIBC
 IJE=D%NJEC
 IJB=D%NJBC
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = -(1.+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))*(PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*PETHETA(IIB:IIE,IJB:IJE,1:D%NKT)/PD(IIB:IIE,IJB:IJE,1:D%NKT))**2*PDRDZ(IIB:IIE,IJB:IJE,1:D%NKT)*(1.5+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))
+ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = -(1.+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))*(PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT) & 
+                                 *PETHETA(IIB:IIE,IJB:IJE,1:D%NKT)/PD(IIB:IIE,IJB:IJE,1:D%NKT))**2&
+                                 *PDRDZ(IIB:IIE,IJB:IJE,1:D%NKT)&
+                                 *(1.5+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 CALL MZF_PHY(D,ZWORK1,ZWORK2)
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-D_M3_THR_WTH2_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = - 0.25*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/PSQRT_TKE(IIB:IIE,IJB:IJE,1:D%NKT)/CSTURB%XCTD*CSTURB%XCTV**2 &
-                               * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT) 
+D_M3_THR_WTH2_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = - 0.25*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT) &
+                                 /PSQRT_TKE(IIB:IIE,IJB:IJE,1:D%NKT)/CSTURB%XCTD*CSTURB%XCTV**2 * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT) 
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 !
 D_M3_THR_WTH2_O_DDTDZ(IIB:IIE,IJB:IJE,IKB-1)=D_M3_THR_WTH2_O_DDTDZ(IIB:IIE,IJB:IJE,IKB)
@@ -1923,13 +1979,14 @@ IIB=D%NIBC
 IJE=D%NJEC
 IJB=D%NJBC
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*PETHETA(IIB:IIE,IJB:IJE,1:D%NKT)/PD(IIB:IIE,IJB:IJE,1:D%NKT)                                              &
-       *(-(1.+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))*PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)/PD(IIB:IIE,IJB:IJE,1:D%NKT)*(1.5+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))+(1.+2.*PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)))
+ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*PETHETA(IIB:IIE,IJB:IJE,1:D%NKT)/PD(IIB:IIE,IJB:IJE,1:D%NKT)&
+       *(-(1.+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))*PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)/PD(IIB:IIE,IJB:IJE,1:D%NKT)&
+       *(1.5+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))+(1.+2.*PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)))
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 CALL MZF_PHY(D,ZWORK1,ZWORK2)
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-D_M3_THR_WTH2_O_DDRDZ(IIB:IIE,IJB:IJE,1:D%NKT) = - 0.25*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/PSQRT_TKE(IIB:IIE,IJB:IJE,1:D%NKT)/CSTURB%XCTD*CSTURB%XCTV          &
-                               * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
+D_M3_THR_WTH2_O_DDRDZ(IIB:IIE,IJB:IJE,1:D%NKT) = - 0.25*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/PSQRT_TKE(IIB:IIE,IJB:IJE,1:D%NKT)& 
+                                                / CSTURB%XCTD*CSTURB%XCTV * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 !
 D_M3_THR_WTH2_O_DDRDZ(IIB:IIE,IJB:IJE,IKB-1)=D_M3_THR_WTH2_O_DDRDZ(IIB:IIE,IJB:IJE,IKB)
@@ -1964,8 +2021,8 @@ ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = (1.+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))*PDRDZ(IIB
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 CALL MZF_PHY(D,ZWORK1,ZWORK2)
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-M3_THR_W2TH(IIB:IIE,IJB:IJE,1:D%NKT) = - 0.75*PLM(IIB:IIE,IJB:IJE,1:D%NKT)*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/PTKE(IIB:IIE,IJB:IJE,1:D%NKT) * CSTURB%XCTV      &
-                     * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
+M3_THR_W2TH(IIB:IIE,IJB:IJE,1:D%NKT) = - 0.75*PLM(IIB:IIE,IJB:IJE,1:D%NKT)*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)& 
+                                      / PTKE(IIB:IIE,IJB:IJE,1:D%NKT) * CSTURB%XCTV * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 !
 M3_THR_W2TH(IIB:IIE,IJB:IJE,IKB-1)=M3_THR_W2TH(IIB:IIE,IJB:IJE,IKB)
@@ -1999,12 +2056,14 @@ IIB=D%NIBC
 IJE=D%NJEC
 IJB=D%NJBC
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) =  -PETHETA(IIB:IIE,IJB:IJE,1:D%NKT)*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*(1.+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))*PDRDZ(IIB:IIE,IJB:IJE,1:D%NKT)*(1.5+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))/PD(IIB:IIE,IJB:IJE,1:D%NKT)**2
+ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) =  -PETHETA(IIB:IIE,IJB:IJE,1:D%NKT)*PBLL_O_E(IIB:IIE,IJB:IJE,1:D%NKT)*& 
+(1.+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))*PDRDZ(IIB:IIE,IJB:IJE,1:D%NKT)& 
+*(1.5+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))/PD(IIB:IIE,IJB:IJE,1:D%NKT)**2
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 CALL MZF_PHY(D,ZWORK1,ZWORK2)
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-D_M3_THR_W2TH_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = - 0.75*PLM(IIB:IIE,IJB:IJE,1:D%NKT)*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/PTKE(IIB:IIE,IJB:IJE,1:D%NKT) * CSTURB%XCTV**2    &
-                               * ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)
+D_M3_THR_W2TH_O_DDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = - 0.75*PLM(IIB:IIE,IJB:IJE,1:D%NKT)*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)&
+                                                / PTKE(IIB:IIE,IJB:IJE,1:D%NKT) * CSTURB%XCTV**2 * ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 !
 D_M3_THR_W2TH_O_DDTDZ(IIB:IIE,IJB:IJE,IKB-1)=D_M3_THR_W2TH_O_DDTDZ(IIB:IIE,IJB:IJE,IKB)
@@ -2035,13 +2094,14 @@ IIB=D%NIBC
 IJE=D%NJEC
 IJB=D%NJBC
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) =  -(1.+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))*PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)*(1.5+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))/PD(IIB:IIE,IJB:IJE,1:D%NKT)**2          &
+ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) =  -(1.+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))*PREDR1(IIB:IIE,IJB:IJE,1:D%NKT)&
+* (1.5+PREDTH1(IIB:IIE,IJB:IJE,1:D%NKT)+PREDR1(IIB:IIE,IJB:IJE,1:D%NKT))/PD(IIB:IIE,IJB:IJE,1:D%NKT)**2          &
         +(1.+2.*PREDR1(:,:,:))/PD(:,:,:)
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 CALL MZF_PHY(D,ZWORK1,ZWORK2)
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-D_M3_THR_W2TH_O_DDRDZ(IIB:IIE,IJB:IJE,1:D%NKT) = - 0.75*PLM(IIB:IIE,IJB:IJE,1:D%NKT)*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)/PTKE(IIB:IIE,IJB:IJE,1:D%NKT) * CSTURB%XCTV     &
-                               * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
+D_M3_THR_W2TH_O_DDRDZ(IIB:IIE,IJB:IJE,1:D%NKT) = - 0.75*PLM(IIB:IIE,IJB:IJE,1:D%NKT)*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)&
+                                                / PTKE(IIB:IIE,IJB:IJE,1:D%NKT) * CSTURB%XCTV * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 !
 D_M3_THR_W2TH_O_DDRDZ(IIB:IIE,IJB:IJE,IKB-1)=D_M3_THR_W2TH_O_DDRDZ(IIB:IIE,IJB:IJE,IKB)
diff --git a/src/common/turb/mode_tke_eps_sources.F90 b/src/common/turb/mode_tke_eps_sources.F90
index ce4bc8008..6836aa2fd 100644
--- a/src/common/turb/mode_tke_eps_sources.F90
+++ b/src/common/turb/mode_tke_eps_sources.F90
@@ -285,9 +285,10 @@ CALL MZM_PHY(D,PRHODJ,ZMWORK2)
 !
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=IKTB:IKTE)
 ZFLX(IIB:IIE,IJB:IJE,IKTB:IKTE) = CSTURB%XCED * SQRT(PTKEM(IIB:IIE,IJB:IJE,IKTB:IKTE)) / PLEPS(IIB:IIE,IJB:IJE,IKTB:IKTE)
-ZSOURCE(IIB:IIE,IJB:IJE,IKTB:IKTE) = ( PRTKES(IIB:IIE,IJB:IJE,IKTB:IKTE) +  PRTKEMS(IIB:IIE,IJB:IJE,IKTB:IKTE) ) / PRHODJ(IIB:IIE,IJB:IJE,IKTB:IKTE) &
-   - PTKEM(IIB:IIE,IJB:IJE,IKTB:IKTE) / PTSTEP &
-   + PDP(IIB:IIE,IJB:IJE,IKTB:IKTE) + PTP(IIB:IIE,IJB:IJE,IKTB:IKTE) + ZTR(IIB:IIE,IJB:IJE,IKTB:IKTE) - PEXPL * ZFLX(IIB:IIE,IJB:IJE,IKTB:IKTE) * PTKEM(IIB:IIE,IJB:IJE,IKTB:IKTE)
+ZSOURCE(IIB:IIE,IJB:IJE,IKTB:IKTE) = ( PRTKES(IIB:IIE,IJB:IJE,IKTB:IKTE) +  PRTKEMS(IIB:IIE,IJB:IJE,IKTB:IKTE) ) &
+                                     / PRHODJ(IIB:IIE,IJB:IJE,IKTB:IKTE) - PTKEM(IIB:IIE,IJB:IJE,IKTB:IKTE) / PTSTEP &
+   + PDP(IIB:IIE,IJB:IJE,IKTB:IKTE) + PTP(IIB:IIE,IJB:IJE,IKTB:IKTE) + ZTR(IIB:IIE,IJB:IJE,IKTB:IKTE) & 
+   - PEXPL * ZFLX(IIB:IIE,IJB:IJE,IKTB:IKTE) * PTKEM(IIB:IIE,IJB:IJE,IKTB:IKTE)
 !
 !*       2.2  implicit vertical TKE transport
 !
@@ -295,7 +296,8 @@ ZSOURCE(IIB:IIE,IJB:IJE,IKTB:IKTE) = ( PRTKES(IIB:IIE,IJB:IJE,IKTB:IKTE) +  PRTK
 ! Compute the vector giving the elements just under the diagonal for the 
 ! matrix inverted in TRIDIAG 
 !
-ZA(IIB:IIE,IJB:IJE,IKTB:IKTE)     = - PTSTEP * CSTURB%XCET * ZMWORK1(IIB:IIE,IJB:IJE,IKTB:IKTE) * ZMWORK2(IIB:IIE,IJB:IJE,IKTB:IKTE) / PDZZ(IIB:IIE,IJB:IJE,IKTB:IKTE)**2
+ZA(IIB:IIE,IJB:IJE,IKTB:IKTE)     = - PTSTEP * CSTURB%XCET * ZMWORK1(IIB:IIE,IJB:IJE,IKTB:IKTE) & 
+                                    * ZMWORK2(IIB:IIE,IJB:IJE,IKTB:IKTE) / PDZZ(IIB:IIE,IJB:IJE,IKTB:IKTE)**2
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=IKTB:IKTE)
 !
 ! Compute TKE at time t+deltat: ( stored in ZRES )
@@ -328,7 +330,8 @@ IF(HPROGRAM/='MESONH') THEN
 END IF
 !
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=IKTB:IKTE)
-PTDISS(IIB:IIE,IJB:IJE,IKTB:IKTE) = - ZFLX(IIB:IIE,IJB:IJE,IKTB:IKTE)*(PEXPL*PTKEM(IIB:IIE,IJB:IJE,IKTB:IKTE) + PIMPL*ZRES(IIB:IIE,IJB:IJE,IKTB:IKTE))
+PTDISS(IIB:IIE,IJB:IJE,IKTB:IKTE) = - ZFLX(IIB:IIE,IJB:IJE,IKTB:IKTE)*(PEXPL*PTKEM(IIB:IIE,IJB:IJE,IKTB:IKTE) &
+                                      + PIMPL*ZRES(IIB:IIE,IJB:IJE,IKTB:IKTE))
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=IKTB:IKTE)
 !
 IF ( OLES_CALL .OR.                         &
@@ -339,7 +342,8 @@ IF ( OLES_CALL .OR.                         &
   ZMWORK1 = MZM(ZKEFF, D%NKA, D%NKU, D%NKL)
   ZDWORK1 = DZM(PIMPL * ZRES + PEXPL * PTKEM, D%NKA, D%NKU, D%NKL)
   !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=IKTB:IKTE)
-  ZFLX(IIB:IIE,IJB:IJE,IKTB:IKTE)  = - CSTURB%XCET * ZMWORK1(IIB:IIE,IJB:IJE,IKTB:IKTE) * ZDWORK1(IIB:IIE,IJB:IJE,IKTB:IKTE) / PDZZ(IIB:IIE,IJB:IJE,IKTB:IKTE)
+  ZFLX(IIB:IIE,IJB:IJE,IKTB:IKTE)  = - CSTURB%XCET * ZMWORK1(IIB:IIE,IJB:IJE,IKTB:IKTE) &
+                                     * ZDWORK1(IIB:IIE,IJB:IJE,IKTB:IKTE) / PDZZ(IIB:IIE,IJB:IJE,IKTB:IKTE)
   !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=IKTB:IKTE)
 !
   ZFLX(:,:,IKB) = 0.
@@ -349,7 +353,8 @@ IF ( OLES_CALL .OR.                         &
 !
   ZDWORK1 = DZF(MZM(PRHODJ, D%NKA, D%NKU, D%NKL) * ZFLX / PDZZ, D%NKA, D%NKU, D%NKL)
   !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=IKTB:IKTE)
-  ZTR(IIB:IIE,IJB:IJE,IKTB:IKTE)= ZTR(IIB:IIE,IJB:IJE,IKTB:IKTE) - ZDWORK1(IIB:IIE,IJB:IJE,IKTB:IKTE) /PRHODJ(IIB:IIE,IJB:IJE,IKTB:IKTE)
+  ZTR(IIB:IIE,IJB:IJE,IKTB:IKTE)= ZTR(IIB:IIE,IJB:IJE,IKTB:IKTE) - ZDWORK1(IIB:IIE,IJB:IJE,IKTB:IKTE) &
+                                  /PRHODJ(IIB:IIE,IJB:IJE,IKTB:IKTE)
   !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=IKTB:IKTE)
 !
 ! Storage in the LES configuration
@@ -379,18 +384,21 @@ END IF
 !Store the previous source terms in prtkes before initializing the next one
 !Should be in IF LBUDGET_TKE only. Was removed out for a correct comput. of PTDIFF in case of LBUDGET_TKE=F in AROME
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=IKTB:IKTE)
-PRTKES(IIB:IIE,IJB:IJE,IKTB:IKTE) = PRTKES(IIB:IIE,IJB:IJE,IKTB:IKTE) + PRHODJ(IIB:IIE,IJB:IJE,IKTB:IKTE) *                                                           &
-                ( PDP(IIB:IIE,IJB:IJE,IKTB:IKTE) + PTP(IIB:IIE,IJB:IJE,IKTB:IKTE)                                                                 &
-                  - CSTURB%XCED * SQRT(PTKEM(IIB:IIE,IJB:IJE,IKTB:IKTE)) / PLEPS(IIB:IIE,IJB:IJE,IKTB:IKTE) * ( PEXPL*PTKEM(IIB:IIE,IJB:IJE,IKTB:IKTE) + PIMPL*ZRES(IIB:IIE,IJB:IJE,IKTB:IKTE) ) )
-!
-PTDIFF(IIB:IIE,IJB:IJE,IKTB:IKTE) =  ZRES(IIB:IIE,IJB:IJE,IKTB:IKTE) / PTSTEP - PRTKES(IIB:IIE,IJB:IJE,IKTB:IKTE)/PRHODJ(IIB:IIE,IJB:IJE,IKTB:IKTE) &
- & - PDP(IIB:IIE,IJB:IJE,IKTB:IKTE)- PTP(IIB:IIE,IJB:IJE,IKTB:IKTE) - PTDISS(IIB:IIE,IJB:IJE,IKTB:IKTE)
+PRTKES(IIB:IIE,IJB:IJE,IKTB:IKTE) = PRTKES(IIB:IIE,IJB:IJE,IKTB:IKTE) + PRHODJ(IIB:IIE,IJB:IJE,IKTB:IKTE) * &
+                ( PDP(IIB:IIE,IJB:IJE,IKTB:IKTE) + PTP(IIB:IIE,IJB:IJE,IKTB:IKTE)                           &
+                  - CSTURB%XCED * SQRT(PTKEM(IIB:IIE,IJB:IJE,IKTB:IKTE)) / PLEPS(IIB:IIE,IJB:IJE,IKTB:IKTE) &
+                  * ( PEXPL*PTKEM(IIB:IIE,IJB:IJE,IKTB:IKTE) + PIMPL*ZRES(IIB:IIE,IJB:IJE,IKTB:IKTE) ) )
+!
+PTDIFF(IIB:IIE,IJB:IJE,IKTB:IKTE) =  ZRES(IIB:IIE,IJB:IJE,IKTB:IKTE) / PTSTEP - PRTKES(IIB:IIE,IJB:IJE,IKTB:IKTE)&
+                                     /PRHODJ(IIB:IIE,IJB:IJE,IKTB:IKTE) &
+                           & - PDP(IIB:IIE,IJB:IJE,IKTB:IKTE)- PTP(IIB:IIE,IJB:IJE,IKTB:IKTE) - PTDISS(IIB:IIE,IJB:IJE,IKTB:IKTE)
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=IKTB:IKTE)
 !
 IF (BUCONF%LBUDGET_TKE) CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_TKE), 'TR', PRTKES(:, :, :) )
 !
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=IKTB:IKTE)
-PRTKES(IIB:IIE,IJB:IJE,IKTB:IKTE) = ZRES(IIB:IIE,IJB:IJE,IKTB:IKTE) * PRHODJ(IIB:IIE,IJB:IJE,IKTB:IKTE) / PTSTEP -  PRTKEMS(IIB:IIE,IJB:IJE,IKTB:IKTE)
+PRTKES(IIB:IIE,IJB:IJE,IKTB:IKTE) = ZRES(IIB:IIE,IJB:IJE,IKTB:IKTE) * PRHODJ(IIB:IIE,IJB:IJE,IKTB:IKTE) / PTSTEP &
+                                    -  PRTKEMS(IIB:IIE,IJB:IJE,IKTB:IKTE)
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=IKTB:IKTE)
 !
 ! stores the whole turbulent transport
@@ -403,8 +411,10 @@ IF (BUCONF%LBUDGET_TKE) CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_TKE), 'TR', PRTK
 !             -------------------------------
 !
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=IKTB:IKTE)
-PRTHLS(IIB:IIE,IJB:IJE,IKTB:IKTE) = PRTHLS(IIB:IIE,IJB:IJE,IKTB:IKTE) + CSTURB%XCED * SQRT(PTKEM(IIB:IIE,IJB:IJE,IKTB:IKTE)) / PLEPS(IIB:IIE,IJB:IJE,IKTB:IKTE) * &
-                (PEXPL*PTKEM(IIB:IIE,IJB:IJE,IKTB:IKTE) + PIMPL*ZRES(IIB:IIE,IJB:IJE,IKTB:IKTE)) * PRHODJ(IIB:IIE,IJB:IJE,IKTB:IKTE) * PCOEF_DISS(IIB:IIE,IJB:IJE,IKTB:IKTE)
+PRTHLS(IIB:IIE,IJB:IJE,IKTB:IKTE) = PRTHLS(IIB:IIE,IJB:IJE,IKTB:IKTE) + &
+                                    CSTURB%XCED * SQRT(PTKEM(IIB:IIE,IJB:IJE,IKTB:IKTE)) / PLEPS(IIB:IIE,IJB:IJE,IKTB:IKTE) * &
+                (PEXPL*PTKEM(IIB:IIE,IJB:IJE,IKTB:IKTE) + PIMPL*ZRES(IIB:IIE,IJB:IJE,IKTB:IKTE)) &
+                * PRHODJ(IIB:IIE,IJB:IJE,IKTB:IKTE) * PCOEF_DISS(IIB:IIE,IJB:IJE,IKTB:IKTE)
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=IKTB:IKTE)
 !----------------------------------------------------------------------------
 !
diff --git a/src/common/turb/mode_tridiag.F90 b/src/common/turb/mode_tridiag.F90
index a3fc17662..e12da996e 100644
--- a/src/common/turb/mode_tridiag.F90
+++ b/src/common/turb/mode_tridiag.F90
@@ -164,7 +164,8 @@ IJB=D%NJBC
 !
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
 ZY(IIB:IIE,IJB:IJE,IKB) = PVARM(IIB:IIE,IJB:IJE,IKB)  + PTSTEP*PSOURCE(IIB:IIE,IJB:IJE,IKB) -   &
-  PEXPL / PRHODJ(IIB:IIE,IJB:IJE,IKB) * PA(IIB:IIE,IJB:IJE,IKB+D%NKL) * (PVARM(IIB:IIE,IJB:IJE,IKB+D%NKL) - PVARM(IIB:IIE,IJB:IJE,IKB))
+  PEXPL / PRHODJ(IIB:IIE,IJB:IJE,IKB) * PA(IIB:IIE,IJB:IJE,IKB+D%NKL) * &
+  (PVARM(IIB:IIE,IJB:IJE,IKB+D%NKL) - PVARM(IIB:IIE,IJB:IJE,IKB))
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
 !
 DO JK=IKTB+1,IKTE-1
diff --git a/src/common/turb/mode_tridiag_thermo.F90 b/src/common/turb/mode_tridiag_thermo.F90
index 9309ac5dc..e0ed28e13 100644
--- a/src/common/turb/mode_tridiag_thermo.F90
+++ b/src/common/turb/mode_tridiag_thermo.F90
@@ -175,7 +175,8 @@ IJB=D%NJBC
 !
 CALL MZM_PHY(D,PRHODJ,ZMZM_RHODJ)
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-ZRHODJ_DFDDTDZ_O_DZ2(IIB:IIE,IJB:IJE,1:D%NKT) = ZMZM_RHODJ(IIB:IIE,IJB:IJE,1:D%NKT)*PDFDDTDZ(IIB:IIE,IJB:IJE,1:D%NKT)/PDZZ(IIB:IIE,IJB:IJE,1:D%NKT)**2
+ZRHODJ_DFDDTDZ_O_DZ2(IIB:IIE,IJB:IJE,1:D%NKT) = ZMZM_RHODJ(IIB:IIE,IJB:IJE,1:D%NKT)*PDFDDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) &
+                                                /PDZZ(IIB:IIE,IJB:IJE,1:D%NKT)**2
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 !
 ZA=0.
@@ -259,7 +260,8 @@ IF ( PIMPL > 1.E-10 ) THEN
                                                     ! gam(k) = c(k-1) / bet
     ZBET(IIB:IIE,IJB:IJE)    = ZB(IIB:IIE,IJB:IJE,JK) - ZA(IIB:IIE,IJB:IJE,JK) * ZGAM(IIB:IIE,IJB:IJE,JK)
                                                     ! bet = b(k) - a(k)* gam(k)  
-    PVARP(IIB:IIE,IJB:IJE,JK)= ( ZY(IIB:IIE,IJB:IJE,JK) - ZA(IIB:IIE,IJB:IJE,JK) * PVARP(IIB:IIE,IJB:IJE,JK-D%NKL) ) / ZBET(IIB:IIE,IJB:IJE)
+    PVARP(IIB:IIE,IJB:IJE,JK)= ( ZY(IIB:IIE,IJB:IJE,JK) - ZA(IIB:IIE,IJB:IJE,JK) * PVARP(IIB:IIE,IJB:IJE,JK-D%NKL) ) &
+                               / ZBET(IIB:IIE,IJB:IJE)
                                         ! res(k) = (y(k) -a(k)*res(k-1))/ bet 
     !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
   END DO 
@@ -269,7 +271,8 @@ IF ( PIMPL > 1.E-10 ) THEN
                                                     ! gam(k) = c(k-1) / bet
   ZBET(IIB:IIE,IJB:IJE)     = ZB(IIB:IIE,IJB:IJE,IKE) - ZA(IIB:IIE,IJB:IJE,IKE) * ZGAM(IIB:IIE,IJB:IJE,IKE)
                                                     ! bet = b(k) - a(k)* gam(k)  
-  PVARP(IIB:IIE,IJB:IJE,IKE)= ( ZY(IIB:IIE,IJB:IJE,IKE) - ZA(IIB:IIE,IJB:IJE,IKE) * PVARP(IIB:IIE,IJB:IJE,IKE-D%NKL) ) / ZBET(IIB:IIE,IJB:IJE)
+  PVARP(IIB:IIE,IJB:IJE,IKE)= ( ZY(IIB:IIE,IJB:IJE,IKE) - ZA(IIB:IIE,IJB:IJE,IKE) * PVARP(IIB:IIE,IJB:IJE,IKE-D%NKL) ) &
+                              / ZBET(IIB:IIE,IJB:IJE)
                                        ! res(k) = (y(k) -a(k)*res(k-1))/ bet 
   !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
 !
diff --git a/src/common/turb/mode_tridiag_tke.F90 b/src/common/turb/mode_tridiag_tke.F90
index 0832e1060..cd2fad11d 100644
--- a/src/common/turb/mode_tridiag_tke.F90
+++ b/src/common/turb/mode_tridiag_tke.F90
@@ -164,7 +164,8 @@ IJB=D%NJBC
 !
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
 ZY(IIB:IIE,IJB:IJE,IKB) = PVARM(IIB:IIE,IJB:IJE,IKB)  + PTSTEP*PSOURCE(IIB:IIE,IJB:IJE,IKB) -   &
-  PEXPL / PRHODJ(IIB:IIE,IJB:IJE,IKB) * PA(IIB:IIE,IJB:IJE,IKB+D%NKL) * (PVARM(IIB:IIE,IJB:IJE,IKB+D%NKL) - PVARM(IIB:IIE,IJB:IJE,IKB))
+  PEXPL / PRHODJ(IIB:IIE,IJB:IJE,IKB) * PA(IIB:IIE,IJB:IJE,IKB+D%NKL) *  & 
+  (PVARM(IIB:IIE,IJB:IJE,IKB+D%NKL) - PVARM(IIB:IIE,IJB:IJE,IKB))
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
 !
 DO JK=IKTB+1,IKTE-1
diff --git a/src/common/turb/mode_tridiag_wind.F90 b/src/common/turb/mode_tridiag_wind.F90
index 29c06958b..fa0816c19 100644
--- a/src/common/turb/mode_tridiag_wind.F90
+++ b/src/common/turb/mode_tridiag_wind.F90
@@ -169,7 +169,8 @@ IJB=D%NJBC
 !
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
 ZY(IIB:IIE,IJB:IJE,IKB) = PVARM(IIB:IIE,IJB:IJE,IKB)  + PTSTEP*PSOURCE(IIB:IIE,IJB:IJE,IKB) -   &
-  PEXPL / PRHODJA(IIB:IIE,IJB:IJE,IKB) * PA(IIB:IIE,IJB:IJE,IKB+D%NKL) * (PVARM(IIB:IIE,IJB:IJE,IKB+D%NKL) - PVARM(IIB:IIE,IJB:IJE,IKB))
+  PEXPL / PRHODJA(IIB:IIE,IJB:IJE,IKB) * PA(IIB:IIE,IJB:IJE,IKB+D%NKL) * &
+  (PVARM(IIB:IIE,IJB:IJE,IKB+D%NKL) - PVARM(IIB:IIE,IJB:IJE,IKB))
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
 !
 DO JK=IKTB+1,IKTE-1
diff --git a/src/common/turb/mode_turb_ver_dyn_flux.F90 b/src/common/turb/mode_turb_ver_dyn_flux.F90
index a737845d2..dfd5ecb2e 100644
--- a/src/common/turb/mode_turb_ver_dyn_flux.F90
+++ b/src/common/turb/mode_turb_ver_dyn_flux.F90
@@ -371,7 +371,8 @@ ZDIRSINZW(IIB:IIE,IJB:IJE) = SQRT(1.-PDIRCOSZW(IIB:IIE,IJB:IJE)**2)
 !
 IF (OHARAT) THEN
   !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-  ZKEFF(IIB:IIE,IJB:IJE,1:D%NKT) =  PLM(IIB:IIE,IJB:IJE,1:D%NKT) * SQRT(PTKEM(IIB:IIE,IJB:IJE,1:D%NKT)) + 50*MFMOIST(IIB:IIE,IJB:IJE,1:D%NKT)
+  ZKEFF(IIB:IIE,IJB:IJE,1:D%NKT) =  PLM(IIB:IIE,IJB:IJE,1:D%NKT) * SQRT(PTKEM(IIB:IIE,IJB:IJE,1:D%NKT)) + & 
+                                    50*MFMOIST(IIB:IIE,IJB:IJE,1:D%NKT)
   !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 ELSE
   !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
@@ -398,7 +399,8 @@ CALL MXM_PHY(D,PDZZ,ZWORK2)
 CALL MZM_PHY(D,PRHODJ,ZWORK3)
 CALL MXM_PHY(D,ZWORK3,ZWORK4)
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-ZA(IIB:IIE,IJB:IJE,1:D%NKT) = -PTSTEP * ZCMFS * ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)* ZWORK4(IIB:IIE,IJB:IJE,1:D%NKT) / ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)**2
+ZA(IIB:IIE,IJB:IJE,1:D%NKT) = -PTSTEP * ZCMFS * ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)* ZWORK4(IIB:IIE,IJB:IJE,1:D%NKT) &
+                              / ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)**2
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 !
 !
@@ -493,14 +495,16 @@ ZWORK3(IIB:IIE,IJB:IJE,1:D%NKT)=PIMPL*ZRES(IIB:IIE,IJB:IJE,1:D%NKT) + PEXPL*PUM(
 CALL DZM_PHY(D,ZWORK3,ZWORK4)
 CALL MXM_PHY(D,PDZZ,ZWORK5)
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-PRUS(IIB:IIE,IJB:IJE,1:D%NKT)=PRUS(IIB:IIE,IJB:IJE,1:D%NKT)+ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)*(ZRES(IIB:IIE,IJB:IJE,1:D%NKT)-PUM(IIB:IIE,IJB:IJE,1:D%NKT))/PTSTEP
+PRUS(IIB:IIE,IJB:IJE,1:D%NKT)= PRUS(IIB:IIE,IJB:IJE,1:D%NKT)+ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)*(ZRES(IIB:IIE,IJB:IJE,1:D%NKT) & 
+                              - PUM(IIB:IIE,IJB:IJE,1:D%NKT))/PTSTEP
 !
 !
 !*       5.2  Partial Dynamic Production
 !
 ! vertical flux of the U wind component
 !
-ZFLXZ(IIB:IIE,IJB:IJE,1:D%NKT)     = -ZCMFS * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT) * ZWORK4(IIB:IIE,IJB:IJE,1:D%NKT) / ZWORK5(IIB:IIE,IJB:IJE,1:D%NKT)
+ZFLXZ(IIB:IIE,IJB:IJE,1:D%NKT)     = -ZCMFS * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT) * ZWORK4(IIB:IIE,IJB:IJE,1:D%NKT) &
+                                   / ZWORK5(IIB:IIE,IJB:IJE,1:D%NKT)
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 !
 ! surface flux
@@ -692,7 +696,8 @@ CALL MYM_PHY(D,PDZZ,ZWORK2)
 CALL MZM_PHY(D,PRHODJ,ZWORK3)
 CALL MYM_PHY(D,ZWORK3,ZWORK4)
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-ZA(IIB:IIE,IJB:IJE,1:D%NKT) = -PTSTEP * ZCMFS * ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)* ZWORK4(IIB:IIE,IJB:IJE,1:D%NKT) / ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)**2
+ZA(IIB:IIE,IJB:IJE,1:D%NKT) = -PTSTEP * ZCMFS * ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)* ZWORK4(IIB:IIE,IJB:IJE,1:D%NKT) & 
+                              / ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)**2
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 !
 !
@@ -785,14 +790,16 @@ ZWORK3(IIB:IIE,IJB:IJE,1:D%NKT)=PIMPL*ZRES(IIB:IIE,IJB:IJE,1:D%NKT) + PEXPL*PUM(
 CALL DZM_PHY(D,ZWORK3,ZWORK4)
 CALL MYM_PHY(D,PDZZ,ZWORK5)
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-PRVS(IIB:IIE,IJB:IJE,1:D%NKT)=PRVS(IIB:IIE,IJB:IJE,1:D%NKT)+ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)*(ZRES(IIB:IIE,IJB:IJE,1:D%NKT)-PVM(IIB:IIE,IJB:IJE,1:D%NKT))/PTSTEP
+PRVS(IIB:IIE,IJB:IJE,1:D%NKT) = PRVS(IIB:IIE,IJB:IJE,1:D%NKT)+ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)*(ZRES(IIB:IIE,IJB:IJE,1:D%NKT)& 
+                               - PVM(IIB:IIE,IJB:IJE,1:D%NKT))/PTSTEP
 !
 !
 !*       6.2  Complete 1D dynamic Production
 !
 !  vertical flux of the V wind component
 !
-ZFLXZ(IIB:IIE,IJB:IJE,1:D%NKT)   = -ZCMFS * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT) * ZWORK4(IIB:IIE,IJB:IJE,1:D%NKT) / ZWORK5(IIB:IIE,IJB:IJE,1:D%NKT)
+ZFLXZ(IIB:IIE,IJB:IJE,1:D%NKT)   = -ZCMFS * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT) * ZWORK4(IIB:IIE,IJB:IJE,1:D%NKT) &
+                                   / ZWORK5(IIB:IIE,IJB:IJE,1:D%NKT)
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 !
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
diff --git a/src/common/turb/mode_turb_ver_sv_flux.F90 b/src/common/turb/mode_turb_ver_sv_flux.F90
index 73c08e7d9..95a83cca6 100644
--- a/src/common/turb/mode_turb_ver_sv_flux.F90
+++ b/src/common/turb/mode_turb_ver_sv_flux.F90
@@ -327,7 +327,8 @@ IJB=D%NJB
 !
 IF (OHARAT) THEN
   !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-  ZKEFF(IIB:IIE,IJB:IJE,IKB:IKE) =  PLM(IIB:IIE,IJB:IJE,IKB:IKE) * SQRT(PTKEM(IIB:IIE,IJB:IJE,IKB:IKE)) + 50.*MFMOIST(IIB:IIE,IJB:IJE,IKB:IKE)
+  ZKEFF(IIB:IIE,IJB:IJE,IKB:IKE) =  PLM(IIB:IIE,IJB:IJE,IKB:IKE) * SQRT(PTKEM(IIB:IIE,IJB:IJE,IKB:IKE)) & 
+                                   + 50.*MFMOIST(IIB:IIE,IJB:IJE,IKB:IKE)
   !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 ELSE
   !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
@@ -355,7 +356,8 @@ DO JSV=1,KSV
 ! Preparation of the arguments for TRIDIAG 
     IF (OHARAT) THEN
       !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)  
-      ZA(IIB:IIE,IJB:IJE,IKB:IKE) = -PTSTEP * ZKEFF(IIB:IIE,IJB:IJE,IKB:IKE) * ZWORK1(IIB:IIE,IJB:IJE,IKB:IKE) / PDZZ(IIB:IIE,IJB:IJE,IKB:IKE)**2
+      ZA(IIB:IIE,IJB:IJE,IKB:IKE) = -PTSTEP * ZKEFF(IIB:IIE,IJB:IJE,IKB:IKE) * ZWORK1(IIB:IIE,IJB:IJE,IKB:IKE) &
+                                   / PDZZ(IIB:IIE,IJB:IJE,IKB:IKE)**2
       !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
     ELSE
       !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
@@ -406,7 +408,8 @@ DO JSV=1,KSV
     CALL MZM_PHY(D,ZWORK1,ZWORK3)
     CALL DZM_PHY(D,ZWORK2,ZWORK4)
     !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-    ZFLXZ(IIB:IIE,IJB:IJE,IKB:IKE) = -ZCSV * PPSI_SV(IIB:IIE,IJB:IJE,IKB:IKE,JSV) * ZWORK3(IIB:IIE,IJB:IJE,IKB:IKE) / PDZZ(IIB:IIE,IJB:IJE,IKB:IKE) * &
+    ZFLXZ(IIB:IIE,IJB:IJE,IKB:IKE) = -ZCSV * PPSI_SV(IIB:IIE,IJB:IJE,IKB:IKE,JSV) * ZWORK3(IIB:IIE,IJB:IJE,IKB:IKE) & 
+                                    / PDZZ(IIB:IIE,IJB:IJE,IKB:IKE) * &
                   ZWORK4(IIB:IIE,IJB:IJE,IKB:IKE)
     !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
     ! surface flux
diff --git a/src/common/turb/mode_turb_ver_thermo_corr.F90 b/src/common/turb/mode_turb_ver_thermo_corr.F90
index 89f8a76a5..2a34f1e14 100644
--- a/src/common/turb/mode_turb_ver_thermo_corr.F90
+++ b/src/common/turb/mode_turb_ver_thermo_corr.F90
@@ -397,7 +397,8 @@ IF (OHARAT) THEN
   PLEPSF(IIB:IIE,IJB:IJE,D%NKT)=0.001
   !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
   !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-  ZKEFF(IIB:IIE,IJB:IJE,1:D%NKT) = PLM(IIB:IIE,IJB:IJE,1:D%NKT) * SQRT(PTKEM(IIB:IIE,IJB:IJE,1:D%NKT)) + 50*MFMOIST(IIB:IIE,IJB:IJE,1:D%NKT)
+  ZKEFF(IIB:IIE,IJB:IJE,1:D%NKT) = PLM(IIB:IIE,IJB:IJE,1:D%NKT) * SQRT(PTKEM(IIB:IIE,IJB:IJE,1:D%NKT)) & 
+                                 + 50*MFMOIST(IIB:IIE,IJB:IJE,1:D%NKT)
   !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 ELSE
   !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
@@ -448,7 +449,8 @@ END IF
     !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
     CALL MZF_PHY(D,ZWORK1,ZWORK2)
     !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-    ZF(IIB:IIE,IJB:IJE,1:D%NKT) = CSTURB%XCTV*PLM(IIB:IIE,IJB:IJE,1:D%NKT)*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)*ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
+    ZF(IIB:IIE,IJB:IJE,1:D%NKT) = CSTURB%XCTV*PLM(IIB:IIE,IJB:IJE,1:D%NKT)*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)& 
+                                 * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
     !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
   ENDIF
   ZDFDDTDZ(:,:,:) = 0.     ! this term, because of discretization, is treated separately
@@ -534,7 +536,8 @@ END IF
   CALL MZF_PHY(D,ZWORK3,ZWORK4)
   !
   !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-  ZFLXZ(IIB:IIE,IJB:IJE,1:D%NKT)   = ZF(IIB:IIE,IJB:IJE,1:D%NKT) + PIMPL * ZDFDDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) * ZWORK4(IIB:IIE,IJB:IJE,1:D%NKT)
+  ZFLXZ(IIB:IIE,IJB:IJE,1:D%NKT)   = ZF(IIB:IIE,IJB:IJE,1:D%NKT) + PIMPL * ZDFDDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) &
+                                    * ZWORK4(IIB:IIE,IJB:IJE,1:D%NKT)
   !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
   !
   ! special case near the ground ( uncentred gradient )
@@ -627,11 +630,13 @@ END IF
     !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
   ELSE
     !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-    ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = 0.5*(PPHI3(IIB:IIE,IJB:IJE,1:D%NKT)+PPSI3(IIB:IIE,IJB:IJE,1:D%NKT))*PDTH_DZ(IIB:IIE,IJB:IJE,1:D%NKT)*PDR_DZ(IIB:IIE,IJB:IJE,1:D%NKT)
+    ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = 0.5*(PPHI3(IIB:IIE,IJB:IJE,1:D%NKT)+PPSI3(IIB:IIE,IJB:IJE,1:D%NKT))& 
+                                      *PDTH_DZ(IIB:IIE,IJB:IJE,1:D%NKT)*PDR_DZ(IIB:IIE,IJB:IJE,1:D%NKT)
     !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
     CALL MZF_PHY(D,ZWORK1,ZWORK2)
     !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-    ZF(IIB:IIE,IJB:IJE,1:D%NKT) = CSTURB%XCTV*PLM(IIB:IIE,IJB:IJE,1:D%NKT)*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)*ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
+    ZF(IIB:IIE,IJB:IJE,1:D%NKT) = CSTURB%XCTV*PLM(IIB:IIE,IJB:IJE,1:D%NKT)*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)& 
+                                 * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
     !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
   ENDIF
     ZDFDDTDZ(:,:,:) = 0.     ! this term, because of discretization, is treated separately
@@ -734,7 +739,8 @@ END IF
     CALL MZF_PHY(D,ZWORK1,ZWORK7)
     CALL MZF_PHY(D,ZWORK2,ZWORK8)
     IF (OHARAT) THEN    
-      ZWORK5(IIB:IIE,IJB:IJE,1:D%NKT) = 2. *PDR_DZ(IIB:IIE,IJB:IJE,1:D%NKT)  *ZWORK3(IIB:IIE,IJB:IJE,1:D%NKT) / PDZZ(IIB:IIE,IJB:IJE,1:D%NKT)               &
+      ZWORK5(IIB:IIE,IJB:IJE,1:D%NKT) = 2. *PDR_DZ(IIB:IIE,IJB:IJE,1:D%NKT)  *ZWORK3(IIB:IIE,IJB:IJE,1:D%NKT) &
+                                       / PDZZ(IIB:IIE,IJB:IJE,1:D%NKT)               &
                + 2. *PDTH_DZ(IIB:IIE,IJB:IJE,1:D%NKT) *ZWORK4(IIB:IIE,IJB:IJE,1:D%NKT) / PDZZ(IIB:IIE,IJB:IJE,1:D%NKT)
       !
       CALL MZF_PHY(D,ZWORK5,ZWORK6)      
@@ -746,21 +752,27 @@ END IF
         + PIMPL * ZDFDDRDZ(IIB:IIE,IJB:IJE,1:D%NKT) * ZWORK8(IIB:IIE,IJB:IJE,1:D%NKT)
        !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
     ELSE
-      ZWKPHIPSI1(:,:,:) = D_PHI3DTDZ_O_DDTDZ(D,CSTURB,PPHI3,PREDTH1,PREDR1,PRED2TH3,PRED2THR3,HTURBDIM,GUSERV) ! d(phi3*dthdz)/ddthdz term
-      ZWKPHIPSI2(:,:,:) = D_PSI3DTDZ_O_DDTDZ(D,CSTURB,PPSI3,PREDR1,PREDTH1,PRED2R3,PRED2THR3,HTURBDIM,GUSERV)  ! d(psi3*dthdz)/ddthdz term
-      ZWKPHIPSI3(:,:,:) = D_PHI3DRDZ_O_DDRDZ(D,CSTURB,PPHI3,PREDTH1,PREDR1,PRED2TH3,PRED2THR3,HTURBDIM,GUSERV) ! d(phi3*drdz )/ddrdz term
-      ZWKPHIPSI4(:,:,:) = D_PSI3DRDZ_O_DDRDZ(D,CSTURB,PPSI3,PREDR1,PREDTH1,PRED2R3,PRED2THR3,HTURBDIM,GUSERV)  ! d(psi3*drdz )/ddrdz term
+      ZWKPHIPSI1(:,:,:) = D_PHI3DTDZ_O_DDTDZ(D,CSTURB,PPHI3,PREDTH1,PREDR1,PRED2TH3,PRED2THR3,HTURBDIM,GUSERV) 
+      ! d(phi3*dthdz)/ddthdz term
+      ZWKPHIPSI2(:,:,:) = D_PSI3DTDZ_O_DDTDZ(D,CSTURB,PPSI3,PREDR1,PREDTH1,PRED2R3,PRED2THR3,HTURBDIM,GUSERV) 
+      ! d(psi3*dthdz)/ddthdz term
+      ZWKPHIPSI3(:,:,:) = D_PHI3DRDZ_O_DDRDZ(D,CSTURB,PPHI3,PREDTH1,PREDR1,PRED2TH3,PRED2THR3,HTURBDIM,GUSERV)
+      ! d(phi3*drdz )/ddrdz term
+      ZWKPHIPSI4(:,:,:) = D_PSI3DRDZ_O_DDRDZ(D,CSTURB,PPSI3,PREDR1,PREDTH1,PRED2R3,PRED2THR3,HTURBDIM,GUSERV)  
+      ! d(psi3*drdz )/ddrdz term
  
       !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-      ZWORK5(IIB:IIE,IJB:IJE,1:D%NKT) = (ZWKPHIPSI1(IIB:IIE,IJB:IJE,1:D%NKT)+ZWKPHIPSI2(IIB:IIE,IJB:IJE,1:D%NKT))*PDR_DZ(IIB:IIE,IJB:IJE,1:D%NKT)*ZWORK3(IIB:IIE,IJB:IJE,1:D%NKT)/PDZZ(IIB:IIE,IJB:IJE,1:D%NKT) &
-                    + (ZWKPHIPSI3(IIB:IIE,IJB:IJE,1:D%NKT) + ZWKPHIPSI4(IIB:IIE,IJB:IJE,1:D%NKT))*PDTH_DZ(IIB:IIE,IJB:IJE,:D%NKT)*ZWORK4(IIB:IIE,IJB:IJE,1:D%NKT)/PDZZ(IIB:IIE,IJB:IJE,1:D%NKT)
+      ZWORK5(IIB:IIE,IJB:IJE,1:D%NKT) = (ZWKPHIPSI1(IIB:IIE,IJB:IJE,1:D%NKT)+ZWKPHIPSI2(IIB:IIE,IJB:IJE,1:D%NKT))& 
+      *PDR_DZ(IIB:IIE,IJB:IJE,1:D%NKT)*ZWORK3(IIB:IIE,IJB:IJE,1:D%NKT)/PDZZ(IIB:IIE,IJB:IJE,1:D%NKT) &
+                    + (ZWKPHIPSI3(IIB:IIE,IJB:IJE,1:D%NKT) + ZWKPHIPSI4(IIB:IIE,IJB:IJE,1:D%NKT)) & 
+                    *PDTH_DZ(IIB:IIE,IJB:IJE,:D%NKT)*ZWORK4(IIB:IIE,IJB:IJE,1:D%NKT)/PDZZ(IIB:IIE,IJB:IJE,1:D%NKT)
       !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
       CALL MZF_PHY(D,ZWORK5,ZWORK6)      
       
       !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-      ZFLXZ(IIB:IIE,IJB:IJE,1:D%NKT)   = ZF(IIB:IIE,IJB:IJE,1:D%NKT)                                                     &
-        + PIMPL * CSTURB%XCTV*PLM(IIB:IIE,IJB:IJE,1:D%NKT)*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)*0.5                                        &
-          * ZWORK6(IIB:IIE,IJB:IJE,1:D%NKT)                                                           &
+      ZFLXZ(IIB:IIE,IJB:IJE,1:D%NKT)   = ZF(IIB:IIE,IJB:IJE,1:D%NKT)                          &
+        + PIMPL * CSTURB%XCTV*PLM(IIB:IIE,IJB:IJE,1:D%NKT)*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)*0.5 &
+          * ZWORK6(IIB:IIE,IJB:IJE,1:D%NKT)                                                   &
         + PIMPL * ZDFDDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) * ZWORK7(IIB:IIE,IJB:IJE,1:D%NKT)         &
         + PIMPL * ZDFDDRDZ(IIB:IIE,IJB:IJE,1:D%NKT) * ZWORK8(IIB:IIE,IJB:IJE,1:D%NKT)
       !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
@@ -769,42 +781,42 @@ END IF
     ! special case near the ground ( uncentred gradient )
     IF (OHARAT) THEN
       !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
-      ZFLXZ(IIB:IIE,IJB:IJE,IKB) =                                            & 
-      (1. )   &
-      *( PEXPL *                                                  &
+      ZFLXZ(IIB:IIE,IJB:IJE,IKB) =                                                            & 
+      (1. )                                                                                   &
+      *( PEXPL *                                                                              &
          ( ZCOEFF(IIB:IIE,IJB:IJE,IKB+2*D%NKL)*PTHLM(IIB:IIE,IJB:IJE,IKB+2*D%NKL)             &
           +ZCOEFF(IIB:IIE,IJB:IJE,IKB+D%NKL  )*PTHLM(IIB:IIE,IJB:IJE,IKB+D%NKL  )             & 
-          +ZCOEFF(IIB:IIE,IJB:IJE,IKB      )*PTHLM(IIB:IIE,IJB:IJE,IKB      ))            &
+          +ZCOEFF(IIB:IIE,IJB:IJE,IKB      )*PTHLM(IIB:IIE,IJB:IJE,IKB      ))                &
         *( ZCOEFF(IIB:IIE,IJB:IJE,IKB+2*D%NKL)*PRM(IIB:IIE,IJB:IJE,IKB+2*D%NKL,1)             &
           +ZCOEFF(IIB:IIE,IJB:IJE,IKB+D%NKL  )*PRM(IIB:IIE,IJB:IJE,IKB+D%NKL,1  )             & 
-          +ZCOEFF(IIB:IIE,IJB:IJE,IKB      )*PRM(IIB:IIE,IJB:IJE,IKB  ,1    ))            &
-        +PIMPL *                                                  &
+          +ZCOEFF(IIB:IIE,IJB:IJE,IKB      )*PRM(IIB:IIE,IJB:IJE,IKB  ,1    ))                &
+        +PIMPL *                                                                              &
          ( ZCOEFF(IIB:IIE,IJB:IJE,IKB+2*D%NKL)*PTHLP(IIB:IIE,IJB:IJE,IKB+2*D%NKL)             &
           +ZCOEFF(IIB:IIE,IJB:IJE,IKB+D%NKL  )*PTHLP(IIB:IIE,IJB:IJE,IKB+D%NKL  )             &
-          +ZCOEFF(IIB:IIE,IJB:IJE,IKB      )*PTHLP(IIB:IIE,IJB:IJE,IKB      ))            &
+          +ZCOEFF(IIB:IIE,IJB:IJE,IKB      )*PTHLP(IIB:IIE,IJB:IJE,IKB      ))                &
         *( ZCOEFF(IIB:IIE,IJB:IJE,IKB+2*D%NKL)*PRP(IIB:IIE,IJB:IJE,IKB+2*D%NKL  )             &
           +ZCOEFF(IIB:IIE,IJB:IJE,IKB+D%NKL  )*PRP(IIB:IIE,IJB:IJE,IKB+D%NKL    )             & 
-          +ZCOEFF(IIB:IIE,IJB:IJE,IKB      )*PRP(IIB:IIE,IJB:IJE,IKB        ))            &
+          +ZCOEFF(IIB:IIE,IJB:IJE,IKB      )*PRP(IIB:IIE,IJB:IJE,IKB        ))                &
        )
       !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
     ELSE 
       !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
-      ZFLXZ(IIB:IIE,IJB:IJE,IKB) =                                            & 
+      ZFLXZ(IIB:IIE,IJB:IJE,IKB) =                                                            & 
       (CSTURB%XCHT1 * PPHI3(IIB:IIE,IJB:IJE,IKB+D%NKL) + CSTURB%XCHT2 * PPSI3(IIB:IIE,IJB:IJE,IKB+D%NKL))   &
-      *( PEXPL *                                                  &
+      *( PEXPL *                                                                              &
          ( ZCOEFF(IIB:IIE,IJB:IJE,IKB+2*D%NKL)*PTHLM(IIB:IIE,IJB:IJE,IKB+2*D%NKL)             &
           +ZCOEFF(IIB:IIE,IJB:IJE,IKB+D%NKL  )*PTHLM(IIB:IIE,IJB:IJE,IKB+D%NKL  )             & 
-          +ZCOEFF(IIB:IIE,IJB:IJE,IKB      )*PTHLM(IIB:IIE,IJB:IJE,IKB      ))            &
+          +ZCOEFF(IIB:IIE,IJB:IJE,IKB      )*PTHLM(IIB:IIE,IJB:IJE,IKB      ))                &
         *( ZCOEFF(IIB:IIE,IJB:IJE,IKB+2*D%NKL)*PRM(IIB:IIE,IJB:IJE,IKB+2*D%NKL,1)             &
           +ZCOEFF(IIB:IIE,IJB:IJE,IKB+D%NKL  )*PRM(IIB:IIE,IJB:IJE,IKB+D%NKL,1  )             & 
-          +ZCOEFF(IIB:IIE,IJB:IJE,IKB      )*PRM(IIB:IIE,IJB:IJE,IKB  ,1    ))            &
-        +PIMPL *                                                  &
+          +ZCOEFF(IIB:IIE,IJB:IJE,IKB      )*PRM(IIB:IIE,IJB:IJE,IKB  ,1    ))                &
+        +PIMPL *                                                                              &
          ( ZCOEFF(IIB:IIE,IJB:IJE,IKB+2*D%NKL)*PTHLP(IIB:IIE,IJB:IJE,IKB+2*D%NKL)             &
           +ZCOEFF(IIB:IIE,IJB:IJE,IKB+D%NKL  )*PTHLP(IIB:IIE,IJB:IJE,IKB+D%NKL  )             &
-          +ZCOEFF(IIB:IIE,IJB:IJE,IKB      )*PTHLP(IIB:IIE,IJB:IJE,IKB      ))            &
+          +ZCOEFF(IIB:IIE,IJB:IJE,IKB      )*PTHLP(IIB:IIE,IJB:IJE,IKB      ))                &
         *( ZCOEFF(IIB:IIE,IJB:IJE,IKB+2*D%NKL)*PRP(IIB:IIE,IJB:IJE,IKB+2*D%NKL  )             &
           +ZCOEFF(IIB:IIE,IJB:IJE,IKB+D%NKL  )*PRP(IIB:IIE,IJB:IJE,IKB+D%NKL    )             & 
-          +ZCOEFF(IIB:IIE,IJB:IJE,IKB      )*PRP(IIB:IIE,IJB:IJE,IKB        ))            &
+          +ZCOEFF(IIB:IIE,IJB:IJE,IKB      )*PRP(IIB:IIE,IJB:IJE,IKB        ))                &
        )
        !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE) 
     ENDIF
@@ -871,7 +883,8 @@ ELSE
   !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
   CALL MZF_PHY(D,ZWORK1,ZWORK2)
   !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-  ZF(IIB:IIE,IJB:IJE,1:D%NKT) = CSTURB%XCTV*PLM(IIB:IIE,IJB:IJE,1:D%NKT)*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)*ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
+  ZF(IIB:IIE,IJB:IJE,1:D%NKT) = CSTURB%XCTV*PLM(IIB:IIE,IJB:IJE,1:D%NKT)*PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)& 
+                                *ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
   !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 ENDIF
     ZDFDDRDZ(:,:,:) = 0.     ! this term, because of discretization, is treated separately
@@ -961,15 +974,16 @@ ENDIF
     CALL MZF_PHY(D,ZWORK5,ZWORK6)
     !
     !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-    ZFLXZ(IIB:IIE,IJB:IJE,1:D%NKT) = ZF(IIB:IIE,IJB:IJE,1:D%NKT)                                                             &
-          + PIMPL * PLMF(IIB:IIE,IJB:IJE,1:D%NKT) *PLEPSF(IIB:IIE,IJB:IJE,1:D%NKT)                                                   &
+    ZFLXZ(IIB:IIE,IJB:IJE,1:D%NKT) = ZF(IIB:IIE,IJB:IJE,1:D%NKT)                   &
+          + PIMPL * PLMF(IIB:IIE,IJB:IJE,1:D%NKT) *PLEPSF(IIB:IIE,IJB:IJE,1:D%NKT) &
             * ZWORK4(IIB:IIE,IJB:IJE,1:D%NKT) &
           + PIMPL * ZDFDDRDZ(IIB:IIE,IJB:IJE,1:D%NKT) * ZWORK6(IIB:IIE,IJB:IJE,1:D%NKT)
     !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
    ELSE
     ZWKPHIPSI1(:,:,:) = D_PSI3DRDZ2_O_DDRDZ(D,CSTURB,PPSI3,PREDR1,PREDTH1,PRED2R3,PRED2THR3,PDR_DZ,HTURBDIM,GUSERV)  
     !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-    ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = ZWKPHIPSI1(IIB:IIE,IJB:IJE,1:D%NKT)*ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT) / PDZZ(IIB:IIE,IJB:IJE,1:D%NKT)
+    ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) = ZWKPHIPSI1(IIB:IIE,IJB:IJE,1:D%NKT)*ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT) &
+                                      / PDZZ(IIB:IIE,IJB:IJE,1:D%NKT)
     !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
     CALL MZF_PHY(D,ZWORK1,ZWORK3)
     !
@@ -979,8 +993,8 @@ ENDIF
     CALL MZF_PHY(D,ZWORK4,ZWORK5)
     !
     !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-    ZFLXZ(IIB:IIE,IJB:IJE,1:D%NKT) = ZF(IIB:IIE,IJB:IJE,1:D%NKT)                                                             &
-          + PIMPL * CSTURB%XCTV*PLM(IIB:IIE,IJB:IJE,1:D%NKT) *PLEPS(IIB:IIE,IJB:IJE,1:D%NKT)                                                 &
+    ZFLXZ(IIB:IIE,IJB:IJE,1:D%NKT) = ZF(IIB:IIE,IJB:IJE,1:D%NKT)                             &
+          + PIMPL * CSTURB%XCTV*PLM(IIB:IIE,IJB:IJE,1:D%NKT) *PLEPS(IIB:IIE,IJB:IJE,1:D%NKT) &
             * ZWORK3(IIB:IIE,IJB:IJE,1:D%NKT) &
           + PIMPL * ZDFDDRDZ(IIB:IIE,IJB:IJE,1:D%NKT) * ZWORK5(IIB:IIE,IJB:IJE,1:D%NKT)
     !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)  
@@ -1023,7 +1037,8 @@ ENDIF
     !
     IF ( KRRL > 0 ) THEN
       !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-      PSIGS(IIB:IIE,IJB:IJE,1:D%NKT) = PSIGS(IIB:IIE,IJB:IJE,1:D%NKT) + PAMOIST(IIB:IIE,IJB:IJE,1:D%NKT) **2 * ZFLXZ(IIB:IIE,IJB:IJE,1:D%NKT)
+      PSIGS(IIB:IIE,IJB:IJE,1:D%NKT) = PSIGS(IIB:IIE,IJB:IJE,1:D%NKT) + PAMOIST(IIB:IIE,IJB:IJE,1:D%NKT) **2 &
+                                     * ZFLXZ(IIB:IIE,IJB:IJE,1:D%NKT)
       !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
     END IF
     ! stores <Rnp Rnp>    
diff --git a/src/common/turb/mode_turb_ver_thermo_flux.F90 b/src/common/turb/mode_turb_ver_thermo_flux.F90
index 591877adb..3c239da9d 100644
--- a/src/common/turb/mode_turb_ver_thermo_flux.F90
+++ b/src/common/turb/mode_turb_ver_thermo_flux.F90
@@ -476,7 +476,8 @@ GUSERV = (KRR/=0)
 IF (OHARAT) THEN
 ! OHARAT so TKE and length scales at half levels!
   !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-  ZKEFF(IIB:IIE,IJB:IJE,1:D%NKT) =  PLM(IIB:IIE,IJB:IJE,1:D%NKT) * SQRT(PTKEM(IIB:IIE,IJB:IJE,1:D%NKT)) +50.*MFMOIST(IIB:IIE,IJB:IJE,1:D%NKT)
+  ZKEFF(IIB:IIE,IJB:IJE,1:D%NKT) =  PLM(IIB:IIE,IJB:IJE,1:D%NKT) * SQRT(PTKEM(IIB:IIE,IJB:IJE,1:D%NKT)) & 
+                                   +50.*MFMOIST(IIB:IIE,IJB:IJE,1:D%NKT)
   !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 ELSE
   !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
@@ -535,7 +536,8 @@ IF (OHARAT) THEN
   !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 ELSE
   !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-  ZF(IIB:IIE,IJB:IJE,1:D%NKT) = -CSTURB%XCSHF*PPHI3(IIB:IIE,IJB:IJE,1:D%NKT)*ZKEFF(IIB:IIE,IJB:IJE,1:D%NKT)*ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)/PDZZ(IIB:IIE,IJB:IJE,1:D%NKT)
+  ZF(IIB:IIE,IJB:IJE,1:D%NKT) = -CSTURB%XCSHF*PPHI3(IIB:IIE,IJB:IJE,1:D%NKT)*ZKEFF(IIB:IIE,IJB:IJE,1:D%NKT)& 
+                                *ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)/PDZZ(IIB:IIE,IJB:IJE,1:D%NKT)
   ZDFDDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) = -CSTURB%XCSHF*ZKEFF(IIB:IIE,IJB:IJE,1:D%NKT)*ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
   !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 END IF
@@ -666,7 +668,8 @@ CALL TRIDIAG_THERMO(D,PTHLM,ZF,ZDFDDTDZ,PTSTEP,PIMPL,PDZZ,&
 ! Compute the equivalent tendency for the conservative potential temperature
 !
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
-ZRWTHL(IIB:IIE,IJB:IJE,1:D%NKT)= PRHODJ(IIB:IIE,IJB:IJE,1:D%NKT)*(PTHLP(IIB:IIE,IJB:IJE,1:D%NKT)-PTHLM(IIB:IIE,IJB:IJE,1:D%NKT))/PTSTEP
+ZRWTHL(IIB:IIE,IJB:IJE,1:D%NKT)= PRHODJ(IIB:IIE,IJB:IJE,1:D%NKT)*(PTHLP(IIB:IIE,IJB:IJE,1:D%NKT)-PTHLM(IIB:IIE,IJB:IJE,1:D%NKT))& 
+                                 /PTSTEP
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
 ! replace the flux by the Leonard terms above ZALT and ZCLD_THOLD
 IF (TURBN%LHGRAD) THEN
@@ -697,7 +700,8 @@ PRTHLS(IIB:IIE,IJB:IJE,1:D%NKT)= PRTHLS(IIB:IIE,IJB:IJE,1:D%NKT)  + ZRWTHL(IIB:I
 !  Conservative potential temperature flux : 
 !
 !
-ZFLXZ(IIB:IIE,IJB:IJE,1:D%NKT)   = ZF(IIB:IIE,IJB:IJE,1:D%NKT) + PIMPL * ZDFDDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)/ PDZZ(IIB:IIE,IJB:IJE,1:D%NKT)
+ZFLXZ(IIB:IIE,IJB:IJE,1:D%NKT)   = ZF(IIB:IIE,IJB:IJE,1:D%NKT) + PIMPL * ZDFDDTDZ(IIB:IIE,IJB:IJE,1:D%NKT) * & 
+                                   ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)/ PDZZ(IIB:IIE,IJB:IJE,1:D%NKT)
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 !
 ! replace the flux by the Leonard terms
@@ -810,16 +814,17 @@ IF(HPROGRAM/='AROME  ') THEN
    IF ( KRRI >= 1 ) THEN
      !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
      PRRS(IIB:IIE,IJB:IJE,1:D%NKT,2) = PRRS(IIB:IIE,IJB:IJE,1:D%NKT,2) -                                        &
-                     PRHODJ(IIB:IIE,IJB:IJE,1:D%NKT)*PATHETA(IIB:IIE,IJB:IJE,1:D%NKT)*2.*PSRCM(IIB:IIE,IJB:IJE,1:D%NKT)*ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT) &
-                     *(1.0-PFRAC_ICE(IIB:IIE,IJB:IJE,1:D%NKT))
+                     PRHODJ(IIB:IIE,IJB:IJE,1:D%NKT)*PATHETA(IIB:IIE,IJB:IJE,1:D%NKT)*2.*PSRCM(IIB:IIE,IJB:IJE,1:D%NKT)& 
+                     *ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT) *(1.0-PFRAC_ICE(IIB:IIE,IJB:IJE,1:D%NKT))
      PRRS(IIB:IIE,IJB:IJE,1:D%NKT,4) = PRRS(IIB:IIE,IJB:IJE,1:D%NKT,4) -                                        &
-                     PRHODJ(IIB:IIE,IJB:IJE,1:D%NKT)*PATHETA(IIB:IIE,IJB:IJE,1:D%NKT)*2.*PSRCM(IIB:IIE,IJB:IJE,1:D%NKT)* ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)&
-                     *PFRAC_ICE(IIB:IIE,IJB:IJE,1:D%NKT)
+                     PRHODJ(IIB:IIE,IJB:IJE,1:D%NKT)*PATHETA(IIB:IIE,IJB:IJE,1:D%NKT)*2.*PSRCM(IIB:IIE,IJB:IJE,1:D%NKT)&
+                     * ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)*PFRAC_ICE(IIB:IIE,IJB:IJE,1:D%NKT)
      !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
    ELSE
      !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
      PRRS(IIB:IIE,IJB:IJE,1:D%NKT,2) = PRRS(IIB:IIE,IJB:IJE,1:D%NKT,2) -                                        &
-                     PRHODJ(IIB:IIE,IJB:IJE,1:D%NKT)*PATHETA(IIB:IIE,IJB:IJE,1:D%NKT)*2.*PSRCM(IIB:IIE,IJB:IJE,1:D%NKT)*ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
+                     PRHODJ(IIB:IIE,IJB:IJE,1:D%NKT)*PATHETA(IIB:IIE,IJB:IJE,1:D%NKT)*2.*PSRCM(IIB:IIE,IJB:IJE,1:D%NKT)&
+                     *ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
      !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
    END IF
  END IF
@@ -888,7 +893,8 @@ IF (KRR /= 0) THEN
  ELSE
   ZWORK2 = D_PSI3DRDZ_O_DDRDZ(D,CSTURB,PPSI3,PREDR1,PREDTH1,PRED2R3,PRED2THR3,HTURBDIM,GUSERV)
   !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
-  ZF(IIB:IIE,IJB:IJE,1:D%NKT) = -CSTURB%XCSHF*PPSI3(IIB:IIE,IJB:IJE,1:D%NKT)*ZKEFF(IIB:IIE,IJB:IJE,1:D%NKT)*ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)/PDZZ(IIB:IIE,IJB:IJE,1:D%NKT)
+  ZF(IIB:IIE,IJB:IJE,1:D%NKT) = -CSTURB%XCSHF*PPSI3(IIB:IIE,IJB:IJE,1:D%NKT)*ZKEFF(IIB:IIE,IJB:IJE,1:D%NKT)& 
+                                *ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)/PDZZ(IIB:IIE,IJB:IJE,1:D%NKT)
   ZDFDDRDZ(IIB:IIE,IJB:IJE,1:D%NKT) = -CSTURB%XCSHF*ZKEFF(IIB:IIE,IJB:IJE,1:D%NKT)*ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)
   !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
  ENDIF
@@ -1021,7 +1027,8 @@ IF (KRR /= 0) THEN
   ! Compute the equivalent tendency for the conservative mixing ratio
   !
   !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
-  ZRWRNP(IIB:IIE,IJB:IJE,1:D%NKT) = PRHODJ(IIB:IIE,IJB:IJE,1:D%NKT)*(PRP(IIB:IIE,IJB:IJE,1:D%NKT)-PRM(IIB:IIE,IJB:IJE,1:D%NKT,1))/PTSTEP
+  ZRWRNP(IIB:IIE,IJB:IJE,1:D%NKT) = PRHODJ(IIB:IIE,IJB:IJE,1:D%NKT)*(PRP(IIB:IIE,IJB:IJE,1:D%NKT)-PRM(IIB:IIE,IJB:IJE,1:D%NKT,1))& 
+                                   /PTSTEP
   !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)    
   !
   ! replace the flux by the Leonard terms above ZALT and ZCLD_THOLD
@@ -1123,7 +1130,7 @@ IF (KRR /= 0) THEN
   !
   CALL MZM_PHY(D,PEMOIST,ZWORK1)
   !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-  PWTHV(IIB:IIE,IJB:IJE,1:D%NKT)   = PWTHV(IIB:IIE,IJB:IJE,1:D%NKT) + ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) * ZFLXZ(IIB:IIE,IJB:IJE,1:D%NKT)
+  PWTHV(IIB:IIE,IJB:IJE,1:D%NKT)=PWTHV(IIB:IIE,IJB:IJE,1:D%NKT) + ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) * ZFLXZ(IIB:IIE,IJB:IJE,1:D%NKT)
   !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
   !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
   PWTHV(IIB:IIE,IJB:IJE,IKB) = PWTHV(IIB:IIE,IJB:IJE,IKB) + PEMOIST(IIB:IIE,IJB:IJE,IKB) * ZFLXZ(IIB:IIE,IJB:IJE,IKB)
@@ -1142,16 +1149,17 @@ IF(HPROGRAM/='AROME  ') THEN
      IF ( KRRI >= 1 ) THEN
        !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
        PRRS(IIB:IIE,IJB:IJE,1:D%NKT,2) = PRRS(IIB:IIE,IJB:IJE,1:D%NKT,2) -                                        &
-                       PRHODJ(IIB:IIE,IJB:IJE,1:D%NKT)*PAMOIST(IIB:IIE,IJB:IJE,1:D%NKT)*2.*PSRCM(IIB:IIE,IJB:IJE,1:D%NKT)*ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)       &
-                       *(1.0-PFRAC_ICE(IIB:IIE,IJB:IJE,1:D%NKT))
+                       PRHODJ(IIB:IIE,IJB:IJE,1:D%NKT)*PAMOIST(IIB:IIE,IJB:IJE,1:D%NKT)*2.*PSRCM(IIB:IIE,IJB:IJE,1:D%NKT)& 
+                       *ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) *(1.0-PFRAC_ICE(IIB:IIE,IJB:IJE,1:D%NKT))
        PRRS(IIB:IIE,IJB:IJE,1:D%NKT,4) = PRRS(IIB:IIE,IJB:IJE,1:D%NKT,4) -                                        &
-                       PRHODJ(IIB:IIE,IJB:IJE,1:D%NKT)*PAMOIST(IIB:IIE,IJB:IJE,1:D%NKT)*2.*PSRCM(IIB:IIE,IJB:IJE,1:D%NKT)*ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)       &
-                       *PFRAC_ICE(IIB:IIE,IJB:IJE,1:D%NKT)
+                       PRHODJ(IIB:IIE,IJB:IJE,1:D%NKT)*PAMOIST(IIB:IIE,IJB:IJE,1:D%NKT)*2.*PSRCM(IIB:IIE,IJB:IJE,1:D%NKT)&
+                       *ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT) *PFRAC_ICE(IIB:IIE,IJB:IJE,1:D%NKT)
        !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
      ELSE
        !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
        PRRS(IIB:IIE,IJB:IJE,1:D%NKT,2) = PRRS(IIB:IIE,IJB:IJE,1:D%NKT,2) -                                        &
-                       PRHODJ(IIB:IIE,IJB:IJE,1:D%NKT)*PAMOIST(IIB:IIE,IJB:IJE,1:D%NKT)*2.*PSRCM(IIB:IIE,IJB:IJE,1:D%NKT)*ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)
+                       PRHODJ(IIB:IIE,IJB:IJE,1:D%NKT)*PAMOIST(IIB:IIE,IJB:IJE,1:D%NKT)*2.*PSRCM(IIB:IIE,IJB:IJE,1:D%NKT)&
+                       *ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)
        !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
      END IF
    END IF
@@ -1195,7 +1203,8 @@ IF ( ((OTURB_FLX .AND. TPFILE%LOPENED) .OR. OLES_CALL) .AND. (KRRL > 0) ) THEN
  CALL DZM_PHY(D,ZWORK1,ZWORK2)
  IF (OHARAT) THEN
   !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-  ZA(IIB:IIE,IJB:IJE,1:D%NKT)   = ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)/ PDZZ(IIB:IIE,IJB:IJE,1:D%NKT) * (-PLM(IIB:IIE,IJB:IJE,1:D%NKT)*PSQRT_TKE(IIB:IIE,IJB:IJE,1:D%NKT))
+  ZA(IIB:IIE,IJB:IJE,1:D%NKT)   = ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)/ PDZZ(IIB:IIE,IJB:IJE,1:D%NKT) * &
+                                 (-PLM(IIB:IIE,IJB:IJE,1:D%NKT)*PSQRT_TKE(IIB:IIE,IJB:IJE,1:D%NKT))
   !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
  ELSE
   !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
@@ -1203,7 +1212,8 @@ IF ( ((OTURB_FLX .AND. TPFILE%LOPENED) .OR. OLES_CALL) .AND. (KRRL > 0) ) THEN
   !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
   CALL MZM_PHY(D,ZWORK1,ZWORK3)
   !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-  ZA(IIB:IIE,IJB:IJE,1:D%NKT)   = ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)/ PDZZ(IIB:IIE,IJB:IJE,1:D%NKT) * (-PPHI3(IIB:IIE,IJB:IJE,1:D%NKT)*ZWORK3(IIB:IIE,IJB:IJE,1:D%NKT)) * CSTURB%XCSHF 
+  ZA(IIB:IIE,IJB:IJE,1:D%NKT)   = ZWORK2(IIB:IIE,IJB:IJE,1:D%NKT)/ PDZZ(IIB:IIE,IJB:IJE,1:D%NKT) * &
+                                  (-PPHI3(IIB:IIE,IJB:IJE,1:D%NKT)*ZWORK3(IIB:IIE,IJB:IJE,1:D%NKT)) * CSTURB%XCSHF 
   !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
  ENDIF
   !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
@@ -1218,7 +1228,8 @@ IF ( ((OTURB_FLX .AND. TPFILE%LOPENED) .OR. OLES_CALL) .AND. (KRRL > 0) ) THEN
   CALL MZM_PHY(D,ZWORK1,ZWORK3)
   CALL MZM_PHY(D,ZWORK2,ZWORK4)
   !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-  ZFLXZ(IIB:IIE,IJB:IJE,1:D%NKT) = ZWORK3(IIB:IIE,IJB:IJE,1:D%NKT)* ZFLXZ(IIB:IIE,IJB:IJE,1:D%NKT) + ZWORK4(IIB:IIE,IJB:IJE,1:D%NKT)* ZA(IIB:IIE,IJB:IJE,1:D%NKT)
+  ZFLXZ(IIB:IIE,IJB:IJE,1:D%NKT) = ZWORK3(IIB:IIE,IJB:IJE,1:D%NKT)* ZFLXZ(IIB:IIE,IJB:IJE,1:D%NKT) &
+                                 + ZWORK4(IIB:IIE,IJB:IJE,1:D%NKT)* ZA(IIB:IIE,IJB:IJE,1:D%NKT)
   !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
   !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE)
   ZFLXZ(IIB:IIE,IJB:IJE,D%NKA) = ZFLXZ(IIB:IIE,IJB:IJE,IKB)
diff --git a/src/common/turb/turb.F90 b/src/common/turb/turb.F90
index f1aaa42a8..eb1fa12b6 100644
--- a/src/common/turb/turb.F90
+++ b/src/common/turb/turb.F90
@@ -575,7 +575,8 @@ IF (KRRL >=1) THEN
 !
     !$mnh_expand_where(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
     WHERE(PRT(IIB:IIE,IJB:IJE,1:D%NKT,2)+PRT(IIB:IIE,IJB:IJE,1:D%NKT,4)>0.0)
-      ZFRAC_ICE(IIB:IIE,IJB:IJE,1:D%NKT) = PRT(IIB:IIE,IJB:IJE,1:D%NKT,4) / ( PRT(IIB:IIE,IJB:IJE,1:D%NKT,2)+PRT(IIB:IIE,IJB:IJE,1:D%NKT,4) )
+      ZFRAC_ICE(IIB:IIE,IJB:IJE,1:D%NKT) = PRT(IIB:IIE,IJB:IJE,1:D%NKT,4) / ( PRT(IIB:IIE,IJB:IJE,1:D%NKT,2) & 
+                                          +PRT(IIB:IIE,IJB:IJE,1:D%NKT,4) )
     END WHERE
     !$mnh_end_expand_where(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 !
@@ -629,12 +630,16 @@ IF ( KRRL >= 1 ) THEN
   IF ( KRRI >= 1 ) THEN
     !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
     ! Rnp at t
-    PRT(IIB:IIE,IJB:IJE,1:D%NKT,1)  = PRT(IIB:IIE,IJB:IJE,1:D%NKT,1)  + PRT(IIB:IIE,IJB:IJE,1:D%NKT,2)  + PRT(IIB:IIE,IJB:IJE,1:D%NKT,4)
-    PRRS(IIB:IIE,IJB:IJE,1:D%NKT,1) = PRRS(IIB:IIE,IJB:IJE,1:D%NKT,1) + PRRS(IIB:IIE,IJB:IJE,1:D%NKT,2) + PRRS(IIB:IIE,IJB:IJE,1:D%NKT,4)
+    PRT(IIB:IIE,IJB:IJE,1:D%NKT,1)  = PRT(IIB:IIE,IJB:IJE,1:D%NKT,1)  + PRT(IIB:IIE,IJB:IJE,1:D%NKT,2)  & 
+                                    + PRT(IIB:IIE,IJB:IJE,1:D%NKT,4)
+    PRRS(IIB:IIE,IJB:IJE,1:D%NKT,1) = PRRS(IIB:IIE,IJB:IJE,1:D%NKT,1) + PRRS(IIB:IIE,IJB:IJE,1:D%NKT,2) & 
+                                    + PRRS(IIB:IIE,IJB:IJE,1:D%NKT,4)
     ! Theta_l at t
-    PTHLT(IIB:IIE,IJB:IJE,1:D%NKT)  = PTHLT(IIB:IIE,IJB:IJE,1:D%NKT)  - ZLVOCPEXNM(IIB:IIE,IJB:IJE,1:D%NKT) * PRT(IIB:IIE,IJB:IJE,1:D%NKT,2) &
+    PTHLT(IIB:IIE,IJB:IJE,1:D%NKT)  = PTHLT(IIB:IIE,IJB:IJE,1:D%NKT)  - ZLVOCPEXNM(IIB:IIE,IJB:IJE,1:D%NKT) &
+                                    * PRT(IIB:IIE,IJB:IJE,1:D%NKT,2) &
                                   - ZLSOCPEXNM(IIB:IIE,IJB:IJE,1:D%NKT) * PRT(IIB:IIE,IJB:IJE,1:D%NKT,4)
-    PRTHLS(IIB:IIE,IJB:IJE,1:D%NKT) = PRTHLS(IIB:IIE,IJB:IJE,1:D%NKT) - ZLVOCPEXNM(IIB:IIE,IJB:IJE,1:D%NKT) * PRRS(IIB:IIE,IJB:IJE,1:D%NKT,2) &
+    PRTHLS(IIB:IIE,IJB:IJE,1:D%NKT) = PRTHLS(IIB:IIE,IJB:IJE,1:D%NKT) - ZLVOCPEXNM(IIB:IIE,IJB:IJE,1:D%NKT) &
+                                    * PRRS(IIB:IIE,IJB:IJE,1:D%NKT,2) &
                                   - ZLSOCPEXNM(IIB:IIE,IJB:IJE,1:D%NKT) * PRRS(IIB:IIE,IJB:IJE,1:D%NKT,4)
     !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
   ELSE
@@ -643,8 +648,10 @@ IF ( KRRL >= 1 ) THEN
     PRT(IIB:IIE,IJB:IJE,1:D%NKT,1)  = PRT(IIB:IIE,IJB:IJE,1:D%NKT,1)  + PRT(IIB:IIE,IJB:IJE,1:D%NKT,2) 
     PRRS(IIB:IIE,IJB:IJE,1:D%NKT,1) = PRRS(IIB:IIE,IJB:IJE,1:D%NKT,1) + PRRS(IIB:IIE,IJB:IJE,1:D%NKT,2)
     ! Theta_l at t
-    PTHLT(IIB:IIE,IJB:IJE,1:D%NKT)  = PTHLT(IIB:IIE,IJB:IJE,1:D%NKT)  - ZLOCPEXNM(IIB:IIE,IJB:IJE,1:D%NKT) * PRT(IIB:IIE,IJB:IJE,1:D%NKT,2)
-    PRTHLS(IIB:IIE,IJB:IJE,1:D%NKT) = PRTHLS(IIB:IIE,IJB:IJE,1:D%NKT) - ZLOCPEXNM(IIB:IIE,IJB:IJE,1:D%NKT) * PRRS(IIB:IIE,IJB:IJE,1:D%NKT,2)
+    PTHLT(IIB:IIE,IJB:IJE,1:D%NKT)  = PTHLT(IIB:IIE,IJB:IJE,1:D%NKT)  - ZLOCPEXNM(IIB:IIE,IJB:IJE,1:D%NKT) &
+                                    * PRT(IIB:IIE,IJB:IJE,1:D%NKT,2)
+    PRTHLS(IIB:IIE,IJB:IJE,1:D%NKT) = PRTHLS(IIB:IIE,IJB:IJE,1:D%NKT) - ZLOCPEXNM(IIB:IIE,IJB:IJE,1:D%NKT) &
+                                    * PRRS(IIB:IIE,IJB:IJE,1:D%NKT,2)
     !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
   END IF
 END IF
@@ -1084,7 +1091,8 @@ END IF
 !      cloud computation is not statistical 
 CALL MZF_PHY(D,PFLXZTHVMF,ZWORK1)
 !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-PTP(IIB:IIE,IJB:IJE,1:D%NKT) = PTP(IIB:IIE,IJB:IJE,1:D%NKT) + CST%XG / PTHVREF(IIB:IIE,IJB:IJE,1:D%NKT) * ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)
+PTP(IIB:IIE,IJB:IJE,1:D%NKT) = PTP(IIB:IIE,IJB:IJE,1:D%NKT) &
+                             + CST%XG / PTHVREF(IIB:IIE,IJB:IJE,1:D%NKT) * ZWORK1(IIB:IIE,IJB:IJE,1:D%NKT)
 !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 
 IF(PRESENT(PTPMF))  THEN
@@ -1207,20 +1215,26 @@ END IF
 IF ( KRRL >= 1 ) THEN
   IF ( KRRI >= 1 ) THEN
     !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-    PRT(IIB:IIE,IJB:IJE,1:D%NKT,1)  = PRT(IIB:IIE,IJB:IJE,1:D%NKT,1)  - PRT(IIB:IIE,IJB:IJE,1:D%NKT,2)  - PRT(IIB:IIE,IJB:IJE,1:D%NKT,4)
-    PRRS(IIB:IIE,IJB:IJE,1:D%NKT,1) = PRRS(IIB:IIE,IJB:IJE,1:D%NKT,1) - PRRS(IIB:IIE,IJB:IJE,1:D%NKT,2) - PRRS(IIB:IIE,IJB:IJE,1:D%NKT,4)
-    PTHLT(IIB:IIE,IJB:IJE,1:D%NKT)  = PTHLT(IIB:IIE,IJB:IJE,1:D%NKT)  + ZLVOCPEXNM(IIB:IIE,IJB:IJE,1:D%NKT) * PRT(IIB:IIE,IJB:IJE,1:D%NKT,2) &
-                                  + ZLSOCPEXNM(IIB:IIE,IJB:IJE,1:D%NKT) * PRT(IIB:IIE,IJB:IJE,1:D%NKT,4)
-    PRTHLS(IIB:IIE,IJB:IJE,1:D%NKT) = PRTHLS(IIB:IIE,IJB:IJE,1:D%NKT) + ZLVOCPEXNM(IIB:IIE,IJB:IJE,1:D%NKT) * PRRS(IIB:IIE,IJB:IJE,1:D%NKT,2) &
-                                  + ZLSOCPEXNM(IIB:IIE,IJB:IJE,1:D%NKT) * PRRS(IIB:IIE,IJB:IJE,1:D%NKT,4)
+    PRT(IIB:IIE,IJB:IJE,1:D%NKT,1)  = PRT(IIB:IIE,IJB:IJE,1:D%NKT,1)  - PRT(IIB:IIE,IJB:IJE,1:D%NKT,2)  &
+                                    - PRT(IIB:IIE,IJB:IJE,1:D%NKT,4)
+    PRRS(IIB:IIE,IJB:IJE,1:D%NKT,1) = PRRS(IIB:IIE,IJB:IJE,1:D%NKT,1) - PRRS(IIB:IIE,IJB:IJE,1:D%NKT,2) &
+                                    - PRRS(IIB:IIE,IJB:IJE,1:D%NKT,4)
+    PTHLT(IIB:IIE,IJB:IJE,1:D%NKT)  = PTHLT(IIB:IIE,IJB:IJE,1:D%NKT)  + ZLVOCPEXNM(IIB:IIE,IJB:IJE,1:D%NKT) &
+                                    * PRT(IIB:IIE,IJB:IJE,1:D%NKT,2) &
+                                    + ZLSOCPEXNM(IIB:IIE,IJB:IJE,1:D%NKT) * PRT(IIB:IIE,IJB:IJE,1:D%NKT,4)
+    PRTHLS(IIB:IIE,IJB:IJE,1:D%NKT) = PRTHLS(IIB:IIE,IJB:IJE,1:D%NKT) + ZLVOCPEXNM(IIB:IIE,IJB:IJE,1:D%NKT) &
+                                    * PRRS(IIB:IIE,IJB:IJE,1:D%NKT,2) &
+                                    + ZLSOCPEXNM(IIB:IIE,IJB:IJE,1:D%NKT) * PRRS(IIB:IIE,IJB:IJE,1:D%NKT,4)
     !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
 !
   ELSE
     !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
     PRT(IIB:IIE,IJB:IJE,1:D%NKT,1)  = PRT(IIB:IIE,IJB:IJE,1:D%NKT,1)  - PRT(IIB:IIE,IJB:IJE,1:D%NKT,2) 
     PRRS(IIB:IIE,IJB:IJE,1:D%NKT,1) = PRRS(IIB:IIE,IJB:IJE,1:D%NKT,1) - PRRS(IIB:IIE,IJB:IJE,1:D%NKT,2)
-    PTHLT(IIB:IIE,IJB:IJE,1:D%NKT)  = PTHLT(IIB:IIE,IJB:IJE,1:D%NKT)  + ZLOCPEXNM(IIB:IIE,IJB:IJE,1:D%NKT) * PRT(IIB:IIE,IJB:IJE,1:D%NKT,2)
-    PRTHLS(IIB:IIE,IJB:IJE,1:D%NKT) = PRTHLS(IIB:IIE,IJB:IJE,1:D%NKT) + ZLOCPEXNM(IIB:IIE,IJB:IJE,1:D%NKT) * PRRS(IIB:IIE,IJB:IJE,1:D%NKT,2)
+    PTHLT(IIB:IIE,IJB:IJE,1:D%NKT)  = PTHLT(IIB:IIE,IJB:IJE,1:D%NKT)  + ZLOCPEXNM(IIB:IIE,IJB:IJE,1:D%NKT) &
+                                    * PRT(IIB:IIE,IJB:IJE,1:D%NKT,2)
+    PRTHLS(IIB:IIE,IJB:IJE,1:D%NKT) = PRTHLS(IIB:IIE,IJB:IJE,1:D%NKT) + ZLOCPEXNM(IIB:IIE,IJB:IJE,1:D%NKT) &
+                                    * PRRS(IIB:IIE,IJB:IJE,1:D%NKT,2)
     !$mnh_end_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
   END IF
 END IF
@@ -1425,7 +1439,8 @@ INTEGER :: JI,JJ,JK
 !*       1.1 Lv/Cph at  t
 !
   !$mnh_expand_array(JI=IIB:IIE,JJ=IJB:IJE,JK=1:D%NKT)
-  PLOCPEXN(IIB:IIE,IJB:IJE,1:D%NKT) = ( PLTT + (CST%XCPV-PC) *  (PT(IIB:IIE,IJB:IJE,1:D%NKT)-CST%XTT) ) / PCP(IIB:IIE,IJB:IJE,1:D%NKT)
+  PLOCPEXN(IIB:IIE,IJB:IJE,1:D%NKT) = ( PLTT + (CST%XCPV-PC) *  (PT(IIB:IIE,IJB:IJE,1:D%NKT)-CST%XTT) ) &
+                                     / PCP(IIB:IIE,IJB:IJE,1:D%NKT)
 !
 !*      1.2 Saturation vapor pressure at t
 !
@@ -1433,7 +1448,8 @@ INTEGER :: JI,JJ,JK
 !
 !*      1.3 saturation  mixing ratio at t
 !
-  ZRVSAT(IIB:IIE,IJB:IJE,1:D%NKT) =  ZRVSAT(IIB:IIE,IJB:IJE,1:D%NKT) * ZEPS / ( PPABST(IIB:IIE,IJB:IJE,1:D%NKT) - ZRVSAT(IIB:IIE,IJB:IJE,1:D%NKT) )
+  ZRVSAT(IIB:IIE,IJB:IJE,1:D%NKT) =  ZRVSAT(IIB:IIE,IJB:IJE,1:D%NKT) &
+                                    * ZEPS / ( PPABST(IIB:IIE,IJB:IJE,1:D%NKT) - ZRVSAT(IIB:IIE,IJB:IJE,1:D%NKT) )
 !
 !*      1.4 compute the saturation mixing ratio derivative (rvs')
 !
@@ -1446,8 +1462,8 @@ INTEGER :: JI,JJ,JK
 !
 !*      1.6 compute Atheta
 !
-  PATHETA(IIB:IIE,IJB:IJE,1:D%NKT)= PAMOIST(IIB:IIE,IJB:IJE,1:D%NKT) * PEXN(IIB:IIE,IJB:IJE,1:D%NKT) *                             &
-        ( ( ZRVSAT(IIB:IIE,IJB:IJE,1:D%NKT) - PRT(IIB:IIE,IJB:IJE,1:D%NKT,1) ) * PLOCPEXN(IIB:IIE,IJB:IJE,1:D%NKT) /               &
+  PATHETA(IIB:IIE,IJB:IJE,1:D%NKT)= PAMOIST(IIB:IIE,IJB:IJE,1:D%NKT) * PEXN(IIB:IIE,IJB:IJE,1:D%NKT) *               &
+        ( ( ZRVSAT(IIB:IIE,IJB:IJE,1:D%NKT) - PRT(IIB:IIE,IJB:IJE,1:D%NKT,1) ) * PLOCPEXN(IIB:IIE,IJB:IJE,1:D%NKT) / &
           ( 1. + ZDRVSATDT(IIB:IIE,IJB:IJE,1:D%NKT) * PLOCPEXN(IIB:IIE,IJB:IJE,1:D%NKT) )        *               &
           (                                                                  &
            ZRVSAT(IIB:IIE,IJB:IJE,1:D%NKT) * (1. + ZRVSAT(IIB:IIE,IJB:IJE,1:D%NKT)/ZEPS)                         &
-- 
GitLab