diff --git a/docs/TODO b/docs/TODO
index 2904850445616917db387ca0428dc3299ea2ef05..aed2597c58fc0f29a8d55cb0455ee5bcf6f8fb06 100644
--- a/docs/TODO
+++ b/docs/TODO
@@ -64,3 +64,7 @@ SPP
Gradients/shuman:
- essayer de mettre des abort dans les routines arome (shuman doit suffire)
+
+Nettoyage apl_arome non fait (pb a la compilation) ==> 4 arguments dans aro_turb_mnh supprimés (non utilisés)
+
+turb.F90 : il reste un CALL à SOURCES_NEG_CORRECT à ajouter. Besoin de récupérer CCLOUD dans apl_arome : comment ?
diff --git a/src/arome/aux/ibm_mixinglength.f90 b/src/arome/aux/ibm_mixinglength.f90
new file mode 100644
index 0000000000000000000000000000000000000000..766627888e3d49a472bd2eaa7a0b49bb58667abe
--- /dev/null
+++ b/src/arome/aux/ibm_mixinglength.f90
@@ -0,0 +1,35 @@
+!MNH_LIC Copyright 2019-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+!-----------------------------------------------------------------
+!
+! ############################
+MODULE MODI_IBM_MIXINGLENGTH
+ ! ############################
+ !
+ INTERFACE
+ !
+ SUBROUTINE IBM_MIXINGLENGTH(PLM,PLEPS,PMU,PHI,PTKE)
+ !
+ REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PLM
+ REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PLEPS
+ REAL, DIMENSION(:,:,:), INTENT(OUT) :: PMU
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PHI
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PTKE
+ !
+ END SUBROUTINE IBM_MIXINGLENGTH
+ !
+ END INTERFACE
+END MODULE MODI_IBM_MIXINGLENGTH
+ !
+ SUBROUTINE IBM_MIXINGLENGTH(PLM,PLEPS,PMU,PHI,PTKE)
+ !
+ REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PLM
+ REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PLEPS
+ REAL, DIMENSION(:,:,:), INTENT(OUT) :: PMU
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PHI
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PTKE
+ !
+ END SUBROUTINE IBM_MIXINGLENGTH
+ !
diff --git a/src/arome/aux/modd_ibm_paramn.f90 b/src/arome/aux/modd_ibm_paramn.f90
new file mode 100644
index 0000000000000000000000000000000000000000..716a93e6d6da60bf1ca8f19e2ecdffa9faf2fc32
--- /dev/null
+++ b/src/arome/aux/modd_ibm_paramn.f90
@@ -0,0 +1,15 @@
+!MNH_LIC Copyright 2019-2021 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+!-----------------------------------------------------------------
+!
+! #######################
+MODULE MODD_IBM_PARAM_n
+ ! #######################
+ IMPLICIT NONE
+ LOGICAL :: LIBM,LIBM_TROUBLE
+ REAL, DIMENSION(:,:,:,:) , POINTER :: XIBM_LS=>NULL() ! LSF for MNH
+ REAL, DIMENSION(:,:,:) , POINTER :: XIBM_XMUT=>NULL()
+END MODULE MODD_IBM_PARAM_n
+!
diff --git a/src/arome/aux/mode_argslist_ll.F90 b/src/arome/aux/mode_argslist_ll.F90
index a01b2324990bdc736fc1884ed585fbf074677ecd..f80bc15fb3ae6e2d9ebdb0c7d3c6796127648e05 100644
--- a/src/arome/aux/mode_argslist_ll.F90
+++ b/src/arome/aux/mode_argslist_ll.F90
@@ -8,6 +8,17 @@ IMPLICIT NONE
TYPE(LIST_ll), POINTER :: TPLIST ! List of fields
CALL ABORT
END SUBROUTINE CLEANLIST_ll
+!
+ SUBROUTINE ADD2DFIELD_ll(TPLIST, PFIELD, HNAME)
+IMPLICIT NONE
+
+ TYPE(LIST_ll), POINTER :: TPLIST ! list of fields
+ REAL, DIMENSION(:,:), TARGET :: PFIELD ! field to be added to the list
+ ! of fields
+ character(len=*), intent(in) :: HNAME ! Name of the field to be added
+ !
+ CALL ABORT
+END SUBROUTINE ADD2DFIELD_ll
!
SUBROUTINE ADD3DFIELD_ll(TPLIST, PFIELD, HNAME)
IMPLICIT NONE
diff --git a/src/arome/ext/aro_turb_mnh.F90 b/src/arome/ext/aro_turb_mnh.F90
index d9609464f953c8421f9a75c5ac136711125cbf2a..db3aea28115b43d27036f277ee2e28e69faf106c 100644
--- a/src/arome/ext/aro_turb_mnh.F90
+++ b/src/arome/ext/aro_turb_mnh.F90
@@ -415,15 +415,14 @@ DO JRR=1, NBUDGET_RI
YLBUDGET(JRR)%YDMDDH=>YDMDDH
ENDDO
-CL=HINST_SFU
CALL TURB (KLEV+2,1,KKL,IMI, KRR, KRRL, KRRI, HLBCX, HLBCY, ISPLIT,IMI, &
& OTURB_FLX,OTURB_DIAG,OSUBG_COND,ORMC01, &
- & HTURBDIM,HTURBLEN,'NONE','NONE', CL, &
- & HMF_UPDRAFT,ZIMPL, &
+ & HTURBDIM,HTURBLEN,'NONE','NONE', &
+ & ZIMPL, &
& 2*PTSTEP,ZTFILE, &
& ZDXX,ZDYY,ZDZZ,ZDZX,ZDZY,ZZZ, &
& ZDIRCOSXW,ZDIRCOSYW,ZDIRCOSZW,ZCOSSLOPE,ZSINSLOPE, &
- & PRHODJ,PTHVREF,PRHODREF, &
+ & PRHODJ,PTHVREF, &
& PSFTH,PSFRV,PSFSV,PSFU,PSFV, &
& PPABSM,PUM,PVM,PWM,PTKEM,ZSVM,PSRCM, &
& PLENGTHM,PLENGTHH,MFMOIST, &
@@ -431,11 +430,10 @@ CALL TURB (KLEV+2,1,KKL,IMI, KRR, KRRL, KRRI, HLBCX, HLBCY, ISPLIT,IMI, &
& ZCEI,ZCEI_MIN,ZCEI_MAX,ZCOEF_AMPL_SAT, &
& PTHM,ZRM, &
& PRUS,PRVS,PRWS,PRTHS,ZRRS,ZRSVS,PRTKES_OUT, &
- & ZHGRAD,PSIGS, &
+ & PSIGS, &
& PDRUS_TURB,PDRVS_TURB, &
& PDRTHLS_TURB,PDRRTS_TURB,ZDRSVS_TURB, &
& PFLXZTHVMF,ZWTH,ZWRC,ZWSV,PDP,PTP,PTPMF,PTDIFF,PTDISS,&
- & YDDDH,YDLDDH,YDMDDH, &
& YLBUDGET, KBUDGETS=SIZE(YLBUDGET), &
& PEDR=PEDR)
!
diff --git a/src/arome/turb/modi_update_lm.F90 b/src/arome/turb/modi_update_lm.F90
deleted file mode 100644
index d2a31786e9a15acb0f239ed06035cbba5a160ac8..0000000000000000000000000000000000000000
--- a/src/arome/turb/modi_update_lm.F90
+++ /dev/null
@@ -1,18 +0,0 @@
-! ######spl
- MODULE MODI_UPDATE_LM
-! ###################
-INTERFACE
-!
-SUBROUTINE UPDATE_LM(HLBCX,HLBCY,PLM,PLEPS)
-!
-CHARACTER(LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX ! X boundary type
-CHARACTER(LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY ! Y boundary type
-!
-REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PLM ! mixing length
-REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PLEPS ! dissipative length
-!
-END SUBROUTINE UPDATE_LM
-!
-END INTERFACE
-!
-END MODULE MODI_UPDATE_LM
diff --git a/src/common/turb/mode_rotate_wind.F90 b/src/common/turb/mode_rotate_wind.F90
new file mode 100644
index 0000000000000000000000000000000000000000..e117b95f81770e26771aee5999a8ef9aedaef688
--- /dev/null
+++ b/src/common/turb/mode_rotate_wind.F90
@@ -0,0 +1,204 @@
+!MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+! #######################
+ MODULE MODE_ROTATE_WIND
+! #######################
+IMPLICIT NONE
+CONTAINS
+! ###########################################################
+ SUBROUTINE ROTATE_WIND(PU,PV,PW, &
+ PDIRCOSXW, PDIRCOSYW, PDIRCOSZW, &
+ PCOSSLOPE,PSINSLOPE, &
+ PDXX,PDYY,PDZZ, &
+ PUSLOPE,PVSLOPE )
+! ###########################################################
+!
+!
+!!**** *ROTATE_WIND* - computes the wind components along the maximum slope
+!! direction and its normal direction in the first mass level.
+!!
+!! PURPOSE
+!! -------
+!!****
+! The purpose of this routine is to compute the wind component parallel
+! to the orography at the first mass level. The exact location where these
+! components are computed is the point of intersection between the normal
+! to the orography and the first mass-level hyper-plane at PDZZ(:,:,IKB)/2
+!
+!!** METHOD
+!! ------
+!! The values of the 3 cartesian components of the wind are determined
+!! by a bilinear interpolation between the 4 nearest points in the first
+!! mass-level hyper-plane. These points are found according to the signs of
+!! the slopes' sinus and cosinus. For each direction of interpolation, the
+!! two different localizations (mass or flux grids) are used to avoid
+!! lateral boundary problems.
+!! Then, the rotation is performed for the wind components. The rotation
+!! angle is the angle between the x axe and the maximum slope direction
+!! defined by the slope vector (dZs/dx , dZs/dy).
+!! Finally, the horizontal components are set at the marginal points
+!! according to cyclic boundary conditions because this is the only case
+!! where these points can be considered.
+!!
+!! EXTERNAL
+!! --------
+!! NONE
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!!
+!! MODD_CONF : L2D switch for 2D model version
+!!
+!!
+!! REFERENCE
+!! ---------
+!! Book 1 of documentation (Chapter: Turbulence)
+!!
+!! AUTHOR
+!! ------
+!! Joel Stein * Meteo-France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original 14/11/95
+!! Modifications: 15/05/96, (N. wood)
+!! take into account no slip conditions
+!! at the surface
+!! 14/02/01 (V. Masson)
+!! Slip condition at the surface restored
+!!
+!! --------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! ------------
+USE MODD_PARAMETERS
+!
+IMPLICIT NONE
+!
+!
+!* 0.1 declarations of arguments
+!
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PU,PV,PW ! cartesian components
+ ! of the wind
+REAL, DIMENSION(:,:), INTENT(IN) :: PDIRCOSXW, PDIRCOSYW, PDIRCOSZW
+! Director Cosinus along x, y and z directions at surface w-point
+REAL, DIMENSION(:,:), INTENT(IN) :: PCOSSLOPE ! cosinus of the angle
+ ! between i and the slope vector
+REAL, DIMENSION(:,:), INTENT(IN) :: PSINSLOPE ! sinus of the angle
+ ! between i and the slope vector
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX, PDYY, PDZZ
+ ! Metric coefficients
+REAL, DIMENSION(:,:), INTENT(OUT) :: PUSLOPE ! wind component along
+ ! the maximum slope direction
+REAL, DIMENSION(:,:), INTENT(OUT) :: PVSLOPE ! wind component along
+ ! the direction normal to the maximum slope one
+!
+!-------------------------------------------------------------------------------
+!
+! 0.2 declaration of local variables
+!
+INTEGER, DIMENSION(SIZE(PDIRCOSXW,1),SIZE(PDIRCOSXW,2)) :: ILOC,JLOC
+ ! shift index to find the 4 nearest points in x and y directions
+REAL, DIMENSION(SIZE(PDIRCOSXW,1),SIZE(PDIRCOSXW,2)) :: ZCOEFF,ZCOEFM, &
+ ! interpolation weigths for flux and mass locations
+ ZUINT,ZVINT,ZWINT, &
+ ! intermediate values of the cartesian components after x interp.
+ ZUFIN,ZVFIN,ZWFIN, &
+ ! final values of the cartesian components after the 2 interp.
+ ZWGROUND
+ ! vertical velocity at the surface
+INTEGER :: IIB,IIE,IJB,IJE,IKB
+ ! index values for the Beginning or the End of the physical
+ ! domain in x,y and z directions
+INTEGER :: IIU,IJU
+ ! arrays' sizes for i and j indices
+INTEGER :: JI,JJ
+!
+!----------------------------------------------------------------------------
+!
+!* 1. PRELIMINARIES
+! -------------
+!
+PUSLOPE=0.
+PVSLOPE=0.
+!
+IIB = 2
+IJB = 2
+IIU = SIZE(PU,1)
+IJU = SIZE(PU,2)
+IIE = IIU - 1
+IJE = IJU - 1
+IKB = 1+JPVEXT
+!
+ZWGROUND(:,:) = PW(:,:,IKB)
+!
+!* 2. INTERPOLATE THE CARTESIAN COMPONENTS
+! ------------------------------------
+!
+ILOC(:,:)=NINT(SIGN(1.,-PCOSSLOPE(:,:)))
+JLOC(:,:)=NINT(SIGN(1.,-PSINSLOPE(:,:)))
+!
+! interpolation in x direction
+!
+DO JJ = 1,IJU
+ DO JI = IIB,IIE
+ ZCOEFF(JI,JJ) = &
+ (0.5*PDXX(JI,JJ,IKB) + 0.5*PDZZ(JI,JJ,IKB)*PDIRCOSXW(JI,JJ) ) &
+ * 2. / (PDXX(JI,JJ,IKB)+PDXX(JI+1,JJ,IKB))
+ ZUINT(JI,JJ) = ZCOEFF(JI,JJ) * PU(JI+1,JJ,IKB) + &
+ (1.-ZCOEFF(JI,JJ)) * PU(JI,JJ,IKB)
+ !
+ ZCOEFM(JI,JJ) = 1. - 0.5 * PDZZ(JI,JJ,IKB) * ABS(PDIRCOSXW(JI,JJ)) &
+ / PDXX(JI+(ILOC(JI,JJ)+1)/2,JJ,IKB)
+ ZVINT(JI,JJ) = ZCOEFM(JI,JJ) * PV(JI,JJ,IKB) + &
+ (1.-ZCOEFM(JI,JJ)) * PV(JI+ILOC(JI,JJ),JJ,IKB)
+ !
+ ZWINT(JI,JJ) = ZCOEFM(JI,JJ) * (PW(JI,JJ,IKB+1)+ZWGROUND(JI,JJ)) * 0.5 &
+ + (1.-ZCOEFM(JI,JJ)) &
+ *(PW(JI+ILOC(JI,JJ),JJ,IKB+1)+ZWGROUND(JI+ILOC(JI,JJ),JJ)) * 0.5
+ END DO
+END DO
+!
+! interpolation in y direction
+!
+DO JJ = IJB,IJE
+ DO JI = IIB,IIE
+ ZCOEFF(JI,JJ) = &
+ (0.5*PDYY(JI,JJ,IKB) + 0.5*PDZZ(JI,JJ,IKB)*PDIRCOSYW(JI,JJ) ) &
+ * 2. / (PDYY(JI,JJ,IKB)+PDYY(JI+1,JJ,IKB))
+ ZVFIN(JI,JJ) = ZCOEFF(JI,JJ) * ZVINT(JI,JJ+1) + &
+ (1.-ZCOEFF(JI,JJ)) * ZVINT(JI,JJ)
+ !
+ ZCOEFM(JI,JJ) = 1. - 0.5 * PDZZ(JI,JJ,IKB) * ABS(PDIRCOSYW(JI,JJ)) &
+ / PDYY(JI,JJ+(JLOC(JI,JJ)+1)/2,IKB)
+ ZUFIN(JI,JJ) = ZCOEFM(JI,JJ) * ZUINT(JI,JJ) + &
+ (1.-ZCOEFM(JI,JJ)) * ZUINT(JI,JJ+JLOC(JI,JJ))
+ ZWFIN(JI,JJ) = ZCOEFM(JI,JJ) * ZWINT(JI,JJ) + &
+ (1.-ZCOEFM(JI,JJ)) * ZWINT(JI,JJ+JLOC(JI,JJ))
+ END DO
+END DO
+!
+!* 3. ROTATE THE WIND
+! ---------------
+!
+!
+DO JJ = IJB,IJE
+ DO JI = IIB,IIE
+ PUSLOPE(JI,JJ) = PCOSSLOPE(JI,JJ) * PDIRCOSZW(JI,JJ) * ZUFIN(JI,JJ) + &
+ PSINSLOPE(JI,JJ) * PDIRCOSZW(JI,JJ) * ZVFIN(JI,JJ) + &
+ SQRT(1.-PDIRCOSZW(JI,JJ)**2) * ZWFIN(JI,JJ)
+ !
+ PVSLOPE(JI,JJ) =-PSINSLOPE(JI,JJ) * ZUFIN(JI,JJ) + &
+ PCOSSLOPE(JI,JJ) * ZVFIN(JI,JJ)
+ !
+ END DO
+END DO
+!
+!
+!
+!----------------------------------------------------------------------------
+!
+END SUBROUTINE ROTATE_WIND
+END MODULE MODE_ROTATE_WIND
diff --git a/src/common/turb/mode_update_lm.F90 b/src/common/turb/mode_update_lm.F90
new file mode 100644
index 0000000000000000000000000000000000000000..d5ab737b0b7d15fd55afbec372f9990a33654f2a
--- /dev/null
+++ b/src/common/turb/mode_update_lm.F90
@@ -0,0 +1,119 @@
+!MNH_LIC Copyright 2006-2019 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+!-----------------------------------------------------------------
+ MODULE MODE_UPDATE_LM
+IMPLICIT NONE
+CONTAINS
+SUBROUTINE UPDATE_LM(HLBCX,HLBCY,PLM,PLEPS)
+! #################################################################
+!
+!!**** *UPDATE_LM* - routine to set external points for mixing length
+!!
+!! PURPOSE
+!! -------
+!
+!!** METHOD
+!! ------
+!!
+!! EXTERNAL
+!! --------
+!!
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!!
+!! REFERENCE
+!! ---------
+!! Book2 of documentation (routine UPDATE_LM)
+!!
+!! AUTHOR
+!! ------
+!! V. Masson * Meteo France *
+!!
+!! MODIFICATIONS
+!! -------------
+!! Original april 2006
+!! V.Masson : Exchange of East and North sides
+!! J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1
+! P. Wautelet 20/05/2019: add name argument to ADDnFIELD_ll + new ADD4DFIELD_ll subroutine
+!-------------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+!
+USE MODD_CONF
+USE MODD_PARAMETERS
+!
+USE MODE_ll
+USE MODD_ARGSLIST_ll, ONLY : LIST_ll
+!
+IMPLICIT NONE
+!
+!
+!* 0.1 declarations of arguments
+!
+CHARACTER(LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX ! X boundary type
+CHARACTER(LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY ! Y boundary type
+!
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PLM ! mixing length
+REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PLEPS ! dissipative length
+!
+!* 0.2 declarations of local variables
+!
+INTEGER :: IIB ! First physical index in x direction
+INTEGER :: IJB ! First physical index in y direction
+INTEGER :: IIE ! last physical index in x direction
+INTEGER :: IJE ! last physical index in y direction
+INTEGER :: JI ! loop index
+!
+TYPE(LIST_ll), POINTER :: TZLM_ll ! list of fields to exchange
+INTEGER :: IINFO_ll ! return code of parallel routine
+!
+!-------------------------------------------------------------------------------
+!
+!* 1. COMPUTE DIMENSIONS OF ARRAYS :
+! ----------------------------
+CALL GET_INDICE_ll (IIB,IJB,IIE,IJE)
+NULLIFY(TZLM_ll)
+!
+!-------------------------------------------------------------------------------
+!
+!* 2. UPDATE HALOs :
+! -------------
+!
+!
+!!$IF(NHALO == 1) THEN
+ CALL ADD3DFIELD_ll( TZLM_ll, PLM, 'UPDATE_LM::PLM' )
+ CALL ADD3DFIELD_ll( TZLM_ll, PLEPS, 'UPDATE_LM::PLEPS' )
+ CALL UPDATE_HALO_ll(TZLM_ll,IINFO_ll)
+ CALL CLEANLIST_ll(TZLM_ll)
+!!$END IF
+!
+!-------------------------------------------------------------------------------
+!
+!* 3. UPDATE EXTERNAL POINTS OF GLOBAL DOMAIN:
+! ---------------------------------------
+!
+IF ( HLBCX(1) /= "CYCL" .AND. LWEST_ll()) THEN
+ PLM (IIB-1,:,:) = PLM (IIB,:,:)
+ PLEPS(IIB-1,:,:) = PLEPS(IIB,:,:)
+END IF
+IF ( HLBCX(2) /= "CYCL" .AND. LEAST_ll()) THEN
+ PLM (IIE+1,:,:) = PLM (IIE,:,:)
+ PLEPS(IIE+1,:,:) = PLEPS(IIE,:,:)
+END IF
+IF ( HLBCY(1) /= "CYCL" .AND. LSOUTH_ll()) THEN
+ DO JI=1,SIZE(PLM,1)
+ PLM (JI,IJB-1,:) = PLM (JI,IJB,:)
+ PLEPS(JI,IJB-1,:) = PLEPS(JI,IJB,:)
+ END DO
+END IF
+IF ( HLBCY(2) /= "CYCL" .AND. LNORTH_ll()) THEN
+ DO JI=1,SIZE(PLM,1)
+ PLM (JI,IJE+1,:) = PLM (JI,IJE,:)
+ PLEPS(JI,IJE+1,:) = PLEPS(JI,IJE,:)
+ END DO
+END IF
+!-----------------------------------------------------------------------------
+END SUBROUTINE UPDATE_LM
+END MODULE MODE_UPDATE_LM
diff --git a/src/common/turb/modi_turb.F90 b/src/common/turb/modi_turb.F90
index 032c257502511d07f5ae84ce9133d9b93d9c9239..ae2a67c527b0e7ff218c6dae44ad35ce58e05323 100644
--- a/src/common/turb/modi_turb.F90
+++ b/src/common/turb/modi_turb.F90
@@ -7,11 +7,11 @@ INTERFACE
SUBROUTINE TURB(KKA, KKU, KKL, KMI,KRR,KRRL,KRRI,HLBCX,HLBCY, &
& KSPLIT,KMODEL_CL, &
& OTURB_FLX,OTURB_DIAG,OSUBG_COND,ORMC01, &
- & HTURBDIM,HTURBLEN,HTOM,HTURBLEN_CL,HINST_SFU, &
- & HMF_UPDRAFT,PIMPL,PTSTEP,TPFILE, &
+ & HTURBDIM,HTURBLEN,HTOM,HTURBLEN_CL, &
+ & PIMPL,PTSTEP,TPFILE, &
& PDXX,PDYY,PDZZ,PDZX,PDZY,PZZ, &
& PDIRCOSXW,PDIRCOSYW,PDIRCOSZW,PCOSSLOPE,PSINSLOPE, &
- & PRHODJ,PTHVREF,PRHODREF, &
+ & PRHODJ,PTHVREF, &
& PSFTH,PSFRV,PSFSV,PSFU,PSFV, &
& PPABST,PUT,PVT,PWT,PTKET,PSVT,PSRCT, &
& PLENGTHM,PLENGTHH,MFMOIST, &
@@ -19,17 +19,13 @@ INTERFACE
& PCEI,PCEI_MIN,PCEI_MAX,PCOEF_AMPL_SAT, &
& PTHLT,PRT, &
& PRUS,PRVS,PRWS,PRTHLS,PRRS,PRSVS,PRTKES, &
- & PHGRAD,PSIGS, &
+ & PSIGS, &
& PDRUS_TURB,PDRVS_TURB, &
& PDRTHLS_TURB,PDRRTS_TURB,PDRSVS_TURB, &
& PFLXZTHVMF,PWTH,PWRC,PWSV,PDP,PTP,PTPMF,PTDIFF,PTDISS,&
- & YDDDH,YDLDDH,YDMDDH, &
& TBUDGETS, KBUDGETS, &
& PTR,PDISS,PEDR,PLEM )
!
-USE DDH_MIX, ONLY : TYP_DDH
-USE YOMLDDH, ONLY : TLDDH
-USE YOMMDDH, ONLY : TMDDH
USE MODD_BUDGET, ONLY : TBUDGETDATA
USE MODD_IO, ONLY : TFILEDATA
!
@@ -41,7 +37,6 @@ INTEGER, INTENT(IN) :: KRR ! number of moist var.
INTEGER, INTENT(IN) :: KRRL ! number of liquid water var.
INTEGER, INTENT(IN) :: KRRI ! number of ice water var.
CHARACTER(LEN=*),DIMENSION(:),INTENT(IN):: HLBCX, HLBCY ! X- and Y-direc LBC
-CHARACTER(LEN=4),INTENT(IN) :: HMF_UPDRAFT ! Type of mass flux
INTEGER, INTENT(IN) :: KSPLIT ! number of time-splitting
INTEGER, INTENT(IN) :: KMODEL_CL ! model number for cloud mixing length
LOGICAL, INTENT(IN) :: OTURB_FLX ! switch to write the
@@ -56,9 +51,8 @@ CHARACTER(LEN=4) , INTENT(IN) :: HTURBDIM ! dimensionality of the
CHARACTER(LEN=4) , INTENT(IN) :: HTURBLEN ! kind of mixing length
CHARACTER(LEN=4) , INTENT(IN) :: HTOM ! kind of Third Order Moment
CHARACTER(LEN=4) , INTENT(IN) :: HTURBLEN_CL ! kind of cloud mixing length
-CHARACTER(LEN=1) , INTENT(IN) :: HINST_SFU ! temporal location of the
- ! surface friction flux
REAL, INTENT(IN) :: PIMPL ! degree of implicitness
+!CHARACTER (LEN=4), INTENT(IN) :: HCLOUD ! Kind of microphysical scheme
REAL, INTENT(IN) :: PTSTEP ! Timestep
TYPE(TFILEDATA), INTENT(IN) :: TPFILE ! Output file for MesoNH
!
@@ -77,8 +71,6 @@ REAL, DIMENSION(:,:,:), INTENT(IN) :: MFMOIST ! Moist mass flux DUal sch
REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHVREF ! Virtual Potential
! Temperature of the reference state
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODREF ! dry density of the
- ! reference state
!
REAL, DIMENSION(:,:), INTENT(IN) :: PSFTH,PSFRV, &
! normal surface fluxes of theta and Rv
@@ -119,7 +111,6 @@ REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRRS
REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRSVS
! Sigma_s at time t+1 : square root of the variance of the deviation to the
! saturation
-REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PHGRAD
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PSIGS
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PDRUS_TURB ! evolution of rhoJ*U by turbulence only
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PDRVS_TURB ! evolution of rhoJ*V by turbulence only
@@ -141,10 +132,6 @@ REAL, DIMENSION(:,:,:), INTENT(OUT) :: PTDISS ! Dissipation TKE term
REAL, DIMENSION(:,:,:), INTENT(IN) :: PLENGTHM
REAL, DIMENSION(:,:,:), INTENT(IN) :: PLENGTHH
!
-TYPE(TYP_DDH), INTENT(INOUT) :: YDDDH
-TYPE(TLDDH), INTENT(IN) :: YDLDDH
-TYPE(TMDDH), INTENT(IN) :: YDMDDH
-!
TYPE(TBUDGETDATA), DIMENSION(KBUDGETS), INTENT(INOUT) :: TBUDGETS
INTEGER, INTENT(IN) :: KBUDGETS
!
diff --git a/src/common/turb/turb.F90 b/src/common/turb/turb.F90
index 2e97d63c673ebe1703d8079c6907c4b0903eb740..f86eaf87ae70fc648f756f90c641b8f8e717c400 100644
--- a/src/common/turb/turb.F90
+++ b/src/common/turb/turb.F90
@@ -6,11 +6,10 @@
SUBROUTINE TURB(KKA,KKU,KKL,KMI,KRR,KRRL,KRRI,HLBCX,HLBCY, &
& KSPLIT,KMODEL_CL, &
& OTURB_FLX,OTURB_DIAG,OSUBG_COND,ORMC01, &
- & HTURBDIM,HTURBLEN,HTOM,HTURBLEN_CL,HINST_SFU, &
- & HMF_UPDRAFT,PIMPL, &
- & PTSTEP,TPFILE, PDXX,PDYY,PDZZ,PDZX,PDZY,PZZ, &
+ & HTURBDIM,HTURBLEN,HTOM,HTURBLEN_CL,PIMPL, &
+ & PTSTEP,TPFILE,PDXX,PDYY,PDZZ,PDZX,PDZY,PZZ, &
& PDIRCOSXW,PDIRCOSYW,PDIRCOSZW,PCOSSLOPE,PSINSLOPE, &
- & PRHODJ,PTHVREF,PRHODREF, &
+ & PRHODJ,PTHVREF, &
& PSFTH,PSFRV,PSFSV,PSFU,PSFV, &
& PPABST,PUT,PVT,PWT,PTKET,PSVT,PSRCT, &
& PLENGTHM,PLENGTHH,MFMOIST, &
@@ -18,11 +17,10 @@
& PCEI,PCEI_MIN,PCEI_MAX,PCOEF_AMPL_SAT, &
& PTHLT,PRT, &
& PRUS,PRVS,PRWS,PRTHLS,PRRS,PRSVS,PRTKES, &
- & PHGRAD, PSIGS, &
+ & PSIGS, &
& PDRUS_TURB,PDRVS_TURB, &
& PDRTHLS_TURB,PDRRTS_TURB,PDRSVS_TURB, &
& PFLXZTHVMF,PWTH,PWRC,PWSV,PDP,PTP,PTPMF,PTDIFF,PTDISS,&
- & YDDDH,YDLDDH,YDMDDH, &
& TBUDGETS, KBUDGETS, &
& PTR,PDISS,PEDR,PLEM )
! #################################################################
@@ -211,11 +209,24 @@
!! 06/2011 (J.escobar ) Bypass Bug with ifort11/12 on HLBCX,HLBC
!! 2012-02 Y. Seity, add possibility to run with reversed
!! vertical levels
+!! 10/2012 (J. Colin) Correct bug in DearDoff for dry simulations
+!! 10/2012 J.Escobar Bypass PGI bug , redefine some allocatable array inplace of automatic
!! 2014-11 Y. Seity, add output terms for TKE DDHs budgets
!! July 2015 (Wim de Rooy) modifications to run with RACMO
!! turbulence (LHARAT=TRUE)
-!! --------------------------------------------------------------------------
-!
+!! 04/2016 (C.Lac) correction of negativity for KHKO
+! P. Wautelet 05/2016-04/2018: new data structures and calls for I/O
+! Q. Rodier 01/2018: introduction of RM17
+! P. Wautelet 20/05/2019: add name argument to ADDnFIELD_ll + new ADD4DFIELD_ll subroutine
+! P. Wautelet 02/2020: use the new data structures and subroutines for budgets
+! B. Vie 03/2020: LIMA negativity checks after turbulence, advection and microphysics budgets
+! P. Wautelet 11/06/2020: bugfix: correct PRSVS array indices
+! P. Wautelet + Benoit Vié 06/2020: improve removal of negative scalar variables + adapt the corresponding budgets
+! P. Wautelet 30/06/2020: move removal of negative scalar variables to Sources_neg_correct
+! R. Honnert/V. Masson 02/2021: new mixing length in the grey zone
+! J.L. Redelsperger 03/2021: add Ocean LES case
+! --------------------------------------------------------------------------
+!
!* 0. DECLARATIONS
! ------------
!
@@ -240,20 +251,18 @@ USE MODD_NSV
!
USE MODE_BL89, ONLY: BL89
USE MODE_TURB_VER, ONLY : TURB_VER
-!!MODIF AROME
-!USE MODI_ROTATE_WIND
+USE MODE_ROTATE_WIND, ONLY: ROTATE_WIND
USE MODE_TURB_HOR_SPLT, ONLY: TURB_HOR_SPLT
USE MODE_TKE_EPS_SOURCES, ONLY: TKE_EPS_SOURCES
USE MODI_SHUMAN, ONLY : MZF, MXF, MYF
USE MODI_GRADIENT_M
USE MODI_GRADIENT_U
USE MODI_GRADIENT_V
-USE MODI_BUDGET_DDH
USE MODI_LES_MEAN_SUBGRID
USE MODE_RMC01, ONLY: RMC01
USE MODI_GRADIENT_W
USE MODE_TM06, ONLY: TM06
-USE MODI_UPDATE_LM
+USE MODE_UPDATE_LM, ONLY: UPDATE_LM
!
USE MODE_BUDGET, ONLY: BUDGET_STORE_INIT, BUDGET_STORE_END
USE MODE_IO_FIELD_WRITE, ONLY: IO_FIELD_WRITE
@@ -263,9 +272,8 @@ USE MODE_SOURCES_NEG_CORRECT, ONLY: SOURCES_NEG_CORRECT
USE MODE_EMOIST, ONLY: EMOIST
USE MODE_ETHETA, ONLY: ETHETA
!
-USE DDH_MIX, ONLY : TYP_DDH
-USE YOMLDDH, ONLY : TLDDH
-USE YOMMDDH, ONLY : TMDDH
+USE MODD_IBM_PARAM_n, ONLY: LIBM, XIBM_LS, XIBM_XMUT
+USE MODI_IBM_MIXINGLENGTH
!
IMPLICIT NONE
!
@@ -296,12 +304,9 @@ CHARACTER(LEN=4), INTENT(IN) :: HTURBDIM ! dimensionality of the
CHARACTER(LEN=4), INTENT(IN) :: HTURBLEN ! kind of mixing length
CHARACTER(LEN=4), INTENT(IN) :: HTOM ! kind of Third Order Moment
CHARACTER(LEN=4), INTENT(IN) :: HTURBLEN_CL ! kind of cloud mixing length
-CHARACTER(LEN=1), INTENT(IN) :: HINST_SFU ! temporal location of the
- ! surface friction flux
REAL, INTENT(IN) :: PIMPL ! degree of implicitness
+!CHARACTER (LEN=4), INTENT(IN) :: HCLOUD ! Kind of microphysical scheme
REAL, INTENT(IN) :: PTSTEP ! timestep
-!
-CHARACTER(LEN=4), INTENT(IN) :: HMF_UPDRAFT ! Type of Mass Flux Scheme
TYPE(TFILEDATA), INTENT(IN) :: TPFILE ! Output file
!
REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX,PDYY,PDZZ,PDZX,PDZY
@@ -319,8 +324,6 @@ REAL, DIMENSION(:,:,:), INTENT(IN) :: MFMOIST ! moist mass flux dual schem
REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHVREF ! Virtual Potential
! Temperature of the reference state
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODREF ! dry density of the
- ! reference state
!
REAL, DIMENSION(:,:), INTENT(IN) :: PSFTH,PSFRV, &
! normal surface fluxes of theta and Rv
@@ -362,7 +365,6 @@ REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRRS
REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRSVS
! Sigma_s at time t+1 : square root of the variance of the deviation to the
! saturation
-REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PHGRAD
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PSIGS
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PDRUS_TURB ! evolution of rhoJ*U by turbulence only
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PDRVS_TURB ! evolution of rhoJ*V by turbulence only
@@ -383,10 +385,6 @@ REAL, DIMENSION(:,:,:), INTENT(OUT) :: PDP ! Dynamic TKE production
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PTDIFF ! Diffusion TKE term
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PTDISS ! Dissipation TKE term
!
-TYPE(TYP_DDH), INTENT(INOUT) :: YDDDH
-TYPE(TLDDH), INTENT(IN) :: YDLDDH
-TYPE(TMDDH), INTENT(IN) :: YDMDDH
-!
TYPE(TBUDGETDATA), DIMENSION(KBUDGETS), INTENT(INOUT) :: TBUDGETS
INTEGER, INTENT(IN) :: KBUDGETS
!
@@ -736,9 +734,15 @@ IF (HTURBLEN=='ADAP') ZLEPS = MIN(ZLEPS,ZLMW*XCADAP)
! ----------------------------------------------------------
!
IF (HTURBDIM=="3DIM") THEN
-!****FOR AROME****
-! CALL UPDATE_LM(HLBCX,HLBCY,ZLM,ZLEPS)
+ CALL UPDATE_LM(HLBCX,HLBCY,ZLM,ZLEPS)
END IF
+!
+!* 3.9 Mixing length correction if immersed walls
+! ------------------------------------------
+!
+IF (LIBM) THEN
+ CALL IBM_MIXINGLENGTH(PLEM,ZLEPS,XIBM_XMUT,XIBM_LS(:,:,:,1),PTKET)
+ENDIF
!----------------------------------------------------------------------------
!
!* 4. GO INTO THE AXES FOLLOWING THE SURFACE
@@ -747,65 +751,30 @@ END IF
!
!* 4.1 rotate the wind at time t
!
-IF ( HINST_SFU == 'T' ) THEN
-!
-!
- IF (CPROGRAM=='AROME ') THEN
- ZUSLOPE=PUT(:,:,KKA)
- ZVSLOPE=PVT(:,:,KKA)
- ELSE
-! CALL ROTATE_WIND(PUT,PVT,PWT, &
-! PDIRCOSXW, PDIRCOSYW, PDIRCOSZW, &
-! PCOSSLOPE,PSINSLOPE, &
-! PDXX,PDYY,PDZZ, &
-! ZUSLOPE,ZVSLOPE )
-!
-! CALL UPDATE_ROTATE_WIND(ZUSLOPE,ZVSLOPE)
- END IF
-!
-!
-!* 4.2 compute the proportionality coefficient between wind and stress
-!
- ZCDUEFF(:,:) =-SQRT ( (PSFU(:,:)**2 + PSFV(:,:)**2) / &
- (1.E-60 + ZUSLOPE(:,:)**2 + ZVSLOPE(:,:)**2 ) &
- )
-!
-!* 4.3 rotate the wind at time t-delta t
!
- IF (CPROGRAM/='AROME ') THEN
-! CALL ROTATE_WIND(PUT,PVT,PWT, &
-! PDIRCOSXW, PDIRCOSYW, PDIRCOSZW, &
-! PCOSSLOPE,PSINSLOPE, &
-! PDXX,PDYY,PDZZ, &
-! ZUSLOPE,ZVSLOPE )
!
-! CALL UPDATE_ROTATE_WIND(ZUSLOPE,ZVSLOPE)
- END IF
+IF (CPROGRAM/='AROME ') THEN
+ CALL ROTATE_WIND(PUT,PVT,PWT, &
+ PDIRCOSXW, PDIRCOSYW, PDIRCOSZW, &
+ PCOSSLOPE,PSINSLOPE, &
+ PDXX,PDYY,PDZZ, &
+ ZUSLOPE,ZVSLOPE )
!
+ CALL UPDATE_ROTATE_WIND(ZUSLOPE,ZVSLOPE)
ELSE
+ ZUSLOPE=PUT(:,:,KKA)
+ ZVSLOPE=PVT(:,:,KKA)
+END IF
!
-!* 4.4 rotate the wind at time t-delta t
-!
- IF (CPROGRAM=='AROME ') THEN
- ZUSLOPE=PUT(:,:,KKA)
- ZVSLOPE=PVT(:,:,KKA)
- ELSE
-!
-! CALL ROTATE_WIND(PUT,PVT,PWT, &
-! PDIRCOSXW, PDIRCOSYW, PDIRCOSZW, &
-! PCOSSLOPE,PSINSLOPE, &
-! PDXX,PDYY,PDZZ, &
-! ZUSLOPE,ZVSLOPE )
-!
-! CALL UPDATE_ROTATE_WIND(ZUSLOPE,ZVSLOPE)
- END IF
!
-!* 4.5 compute the proportionality coefficient between wind and stress
+!* 4.2 compute the proportionality coefficient between wind and stress
!
- ZCDUEFF(:,:) =-SQRT ( (PSFU(:,:)**2 + PSFV(:,:)**2) / &
- (1.E-60 + ZUSLOPE(:,:)**2 + ZVSLOPE(:,:)**2 ) &
- )
-END IF
+ZCDUEFF(:,:) =-SQRT ( (PSFU(:,:)**2 + PSFV(:,:)**2) / &
+#ifdef REPRO48
+ (1.E-60 + ZUSLOPE(:,:)**2 + ZVSLOPE(:,:)**2 ) )
+#else
+ (XMNH_TINY + ZUSLOPE(:,:)**2 + ZVSLOPE(:,:)**2 ) )
+#endif
!
!* 4.6 compute the surface tangential fluxes
!
@@ -948,8 +917,10 @@ IF( LBUDGET_SV ) THEN
END DO
END IF
!
-#ifdef REPRO48 !Les budgets des termes horizontaux de la turb sont présents dans AROME
-#else ! alors que ces termes ne sont pas calculés
+!Les budgets des termes horizontaux de la turb sont présents dans AROME
+! alors que ces termes ne sont pas calculés
+#ifdef REPRO48
+#else
IF( HTURBDIM == '3DIM' ) THEN
#endif
IF( LBUDGET_U ) CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_U ), 'HTURB', PRUS (:, :, :) )
@@ -985,7 +956,8 @@ IF( HTURBDIM == '3DIM' ) THEN
CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_SV1 - 1 + JSV), 'HTURB', PRSVS(:, :, :, JSV) )
END DO
END IF
-#ifdef REPRO48 !à supprimer une fois le précédent ifdef REPRO48 validé
+!à supprimer une fois le précédent ifdef REPRO48 validé
+#ifdef REPRO48
#else
CALL TURB_HOR_SPLT(KSPLIT, KRR, KRRL, KRRI, PTSTEP, &
HLBCX,HLBCY,OTURB_FLX,OSUBG_COND, &
@@ -1170,6 +1142,8 @@ IF ( KRRL >= 1 ) THEN
END IF
END IF
+! Remove non-physical negative values (unnecessary in a perfect world) + corresponding budgets
+!CALL SOURCES_NEG_CORRECT(HCLOUD, 'NETUR',KRR,PTSTEP,PPABST,PTHLT,PRT,PRTHLS,PRRS,PRSVS)
!----------------------------------------------------------------------------
!
!* 9. LES averaged surface fluxes
@@ -1205,8 +1179,8 @@ IF (LLES_CALL) THEN
!
IF (HTURBDIM=="1DIM") THEN
CALL LES_MEAN_SUBGRID(2./3.*PTKET,X_LES_SUBGRID_U2)
- CALL LES_MEAN_SUBGRID(2./3.*PTKET,X_LES_SUBGRID_V2)
- CALL LES_MEAN_SUBGRID(2./3.*PTKET,X_LES_SUBGRID_W2)
+ X_LES_SUBGRID_V2 = X_LES_SUBGRID_U2
+ X_LES_SUBGRID_W2 = X_LES_SUBGRID_U2
CALL LES_MEAN_SUBGRID(2./3.*PTKET*MZF(GZ_M_W(KKA,KKU,KKL,PTHLT,PDZZ),&
KKA, KKU, KKL),X_LES_RES_ddz_Thl_SBG_W2)
IF (KRR>=1) &
@@ -1285,9 +1259,7 @@ IF (LHOOK) CALL DR_HOOK('TURB:UPDATE_ROTATE_WIND',0,ZHOOK_HANDLE)
!
NULLIFY(TZFIELDS_ll)
!
-IIU=SIZE(PUSLOPE,1)
-IJU=SIZE(PUSLOPE,2)
-CALL GET_INDICE_ll (IIB,IJB,IIE,IJE,IIU,IJU)
+CALL GET_INDICE_ll (IIB,IJB,IIE,IJE)
!
! 2 Update halo if necessary
!
diff --git a/src/mesonh/turb/turb.f90 b/src/mesonh/turb/turb.f90
index 8bd2693f0812ac0a89bb1a92ef6fb60872b8f518..3c705ba5b3288c0794eb50353e564a38267396ad 100644
--- a/src/mesonh/turb/turb.f90
+++ b/src/mesonh/turb/turb.f90
@@ -209,6 +209,9 @@
!! vertical levels
!! 10/2012 (J. Colin) Correct bug in DearDoff for dry simulations
!! 10/2012 J.Escobar Bypass PGI bug , redefine some allocatable array inplace of automatic
+!! 2014-11 Y. Seity, add output terms for TKE DDHs budgets
+!! July 2015 (Wim de Rooy) modifications to run with RACMO
+!! turbulence (LHARAT=TRUE)
!! 04/2016 (C.Lac) correction of negativity for KHKO
! P. Wautelet 05/2016-04/2018: new data structures and calls for I/O
! Q. Rodier 01/2018: introduction of RM17