From 9097519f0ef1af5a2cf5ee5e0622d74517ddf858 Mon Sep 17 00:00:00 2001
From: Quentin Rodier <quentin.rodier@meteo.fr>
Date: Tue, 22 Feb 2022 17:11:05 +0100
Subject: [PATCH] Quentin 22/02/2022: Merge from main commit 1230d9

---
 src/mesonh/turb/mf_turb.f90                   |   2 +-
 src/mesonh/turb/mf_turb_greyzone.f90          |   2 +-
 src/mesonh/turb/mode_tridiag_massflux.f90     | 280 ++++++++++++++++++
 src/mesonh/turb/mode_turb_ver_dyn_flux.f90    |  14 +-
 src/mesonh/turb/mode_turb_ver_sv_flux.f90     |   3 +-
 src/mesonh/turb/mode_turb_ver_thermo_corr.f90 |   8 -
 src/mesonh/turb/mode_turb_ver_thermo_flux.f90 |   6 +-
 src/mesonh/turb/turb.f90                      |  20 +-
 8 files changed, 306 insertions(+), 29 deletions(-)
 create mode 100644 src/mesonh/turb/mode_tridiag_massflux.f90

diff --git a/src/mesonh/turb/mf_turb.f90 b/src/mesonh/turb/mf_turb.f90
index 2a96b713a..3ac54d317 100644
--- a/src/mesonh/turb/mf_turb.f90
+++ b/src/mesonh/turb/mf_turb.f90
@@ -137,7 +137,7 @@ END MODULE MODI_MF_TURB
 USE MODD_PARAM_MFSHALL_n
 !
 USE MODI_SHUMAN_MF
-USE MODI_TRIDIAG_MASSFLUX
+USE MODE_TRIDIAG_MASSFLUX, ONLY: TRIDIAG_MASSFLUX
 !
 IMPLICIT NONE
 !
diff --git a/src/mesonh/turb/mf_turb_greyzone.f90 b/src/mesonh/turb/mf_turb_greyzone.f90
index ab28b6c61..be0394221 100644
--- a/src/mesonh/turb/mf_turb_greyzone.f90
+++ b/src/mesonh/turb/mf_turb_greyzone.f90
@@ -138,7 +138,7 @@ END MODULE MODI_MF_TURB_GREYZONE
 USE MODD_PARAM_MFSHALL_n
 !
 USE MODI_SHUMAN_MF
-USE MODI_TRIDIAG_MASSFLUX
+USE MODE_TRIDIAG_MASSFLUX, ONLY: TRIDIAG_MASSFLUX
 !
 IMPLICIT NONE
 !
diff --git a/src/mesonh/turb/mode_tridiag_massflux.f90 b/src/mesonh/turb/mode_tridiag_massflux.f90
new file mode 100644
index 000000000..3e909c2c5
--- /dev/null
+++ b/src/mesonh/turb/mode_tridiag_massflux.f90
@@ -0,0 +1,280 @@
+!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_TRIDIAG_MASSFLUX
+IMPLICIT NONE
+CONTAINS
+SUBROUTINE TRIDIAG_MASSFLUX(KKA,KKB,KKE,KKU,KKL,PVARM,PF,PDFDT,PTSTEP,PIMPL,  &
+                                 PDZZ,PRHODJ,PVARP             )
+
+       USE PARKIND1, ONLY : JPRB
+       USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+!      #################################################
+!
+!
+!!****   *TRIDIAG_MASSFLUX* - routine to solve a time implicit scheme
+!!
+!!
+!!     PURPOSE
+!!     -------
+!        The purpose of this routine is to give a field PVARP at t+1, by 
+!      solving an implicit TRIDIAGonal system obtained by the 
+!      discretization of the vertical turbulent diffusion. It should be noted 
+!      that the degree of implicitness can be varied (PIMPL parameter) and that
+!      the function of F(T) must have been linearized.
+!      PVARP is localized at a mass point.
+!
+!!**   METHOD
+!!     ------
+!!
+!!        [T(+) - T(-)]/2Dt = -d{ F + dF/dT *impl*[T(+) + T(-)] }/dz
+!!
+!!     It is discretized as follows:
+!!
+!!    PRHODJ(k)*PVARP(k)/PTSTEP
+!!              = 
+!!    PRHODJ(k)*PVARM(k)/PTSTEP 
+!!  - (PRHODJ(k+1)+PRHODJ(k)  )/2. * PF(k+1)/PDZZ(k+1)
+!!  + (PRHODJ(k)  +PRHODJ(k-1))/2. * PF(k)  /PDZZ(k)
+!!  + (PRHODJ(k+1)+PRHODJ(k)  )/2. * 0.5*PIMPL* PDFDT(k+1) * PVARM(k+1)/PDZZ(k+1)
+!!  - (PRHODJ(k+1)+PRHODJ(k)  )/2. * 0.5*PIMPL* PDFDT(k+1) * PVARP(k+1)/PDZZ(k+1)
+!!  + (PRHODJ(k+1)+PRHODJ(k)  )/2. * 0.5*PIMPL* PDFDT(k+1) * PVARM(k)  /PDZZ(k+1)
+!!  - (PRHODJ(k+1)+PRHODJ(k)  )/2. * 0.5*PIMPL* PDFDT(k+1) * PVARP(k)  /PDZZ(k+1)
+!!  - (PRHODJ(k)  +PRHODJ(k-1))/2. * 0.5*PIMPL* PDFDT(k)   * PVARM(k)  /PDZZ(k)
+!!  + (PRHODJ(k)  +PRHODJ(k-1))/2. * 0.5*PIMPL* PDFDT(k)   * PVARP(k)  /PDZZ(k)
+!!  - (PRHODJ(k)  +PRHODJ(k-1))/2. * 0.5*PIMPL* PDFDT(k)   * PVARM(k-1)/PDZZ(k)
+!!  + (PRHODJ(k)  +PRHODJ(k-1))/2. * 0.5*PIMPL* PDFDT(k)   * PVARP(k-1)/PDZZ(k)
+!!
+!!
+!!    The system to solve is:
+!!
+!!      A*PVARP(k-1) + B*PVARP(k) + C*PVARP(k+1) = Y(k)
+!!
+!!
+!!    The RHS of the linear system in PVARP writes:
+!!
+!! y(k)    = PRHODJ(k)*PVARM(k)/PTSTEP
+!!  - (PRHODJ(k+1)+PRHODJ(k)  )/2. * PF(k+1)/PDZZ(k+1)
+!!  + (PRHODJ(k)  +PRHODJ(k-1))/2. * PF(k)  /PDZZ(k)
+!!  + (PRHODJ(k+1)+PRHODJ(k)  )/2. * 0.5*PIMPL* PDFDT(k+1) * PVARM(k+1)/PDZZ(k+1)
+!!  + (PRHODJ(k+1)+PRHODJ(k)  )/2. * 0.5*PIMPL* PDFDT(k+1) * PVARM(k)  /PDZZ(k+1)
+!!  - (PRHODJ(k)  +PRHODJ(k-1))/2. * 0.5*PIMPL* PDFDT(k)   * PVARM(k)  /PDZZ(k)
+!!  - (PRHODJ(k)  +PRHODJ(k-1))/2. * 0.5*PIMPL* PDFDT(k)   * PVARM(k-1)/PDZZ(k)
+!!
+!!                      
+!!        Then, the classical TRIDIAGonal algorithm is used to invert the 
+!!     implicit operator. Its matrix is given by:
+!!
+!!     ( b(KKB)   c(KKB)      0        0        0         0        0        0  )
+!!     ( a(KKB+1) b(KKB+1) c(KKB+1)    0  ...    0        0        0        0  ) 
+!!     (   0      a(KKB+2) b(KKB+2) c(KKB+2).    0        0        0        0  ) 
+!!      .......................................................................
+!!     (   0   ...   0     a(k)     b(k)     c(k)         0   ...  0        0  ) 
+!!      .......................................................................
+!!     (   0         0        0        0        0 ...a(KKE-1) b(KKE-1) c(KKE-1))
+!!     (   0         0        0        0        0 ...     0   a(KKE)   b(KKE)  )
+!!
+!!     KKB and KKE represent the first and the last inner mass levels of the
+!!     model. The coefficients are:
+!!         
+!! a(k) = - (PRHODJ(k)  +PRHODJ(k-1))/2. * 0.5*PIMPL* PDFDT(k)  /PDZZ(k)
+!! b(k) =    PRHODJ(k) / PTSTEP
+!!        + (PRHODJ(k+1)+PRHODJ(k)  )/2. * 0.5*PIMPL* PDFDT(k+1)/PDZZ(k+1)
+!!        - (PRHODJ(k)  +PRHODJ(k-1))/2. * 0.5*PIMPL* PDFDT(k)  /PDZZ(k)
+!! c(k) = + (PRHODJ(k+1)+PRHODJ(k)  )/2. * 0.5*PIMPL* PDFDT(k+1)/PDZZ(k+1)
+!!
+!!          for all k /= KKB or KKE
+!!
+!!
+!! b(KKB) =  PRHODJ(KKB) / PTSTEP
+!!          +(PRHODJ(KKB+1)+PRHODJ(KKB))/2.*0.5*PIMPL*PDFDT(KKB+1)/PDZZ(KKB+1)
+!! c(KKB) = +(PRHODJ(KKB+1)+PRHODJ(KKB))/2.*0.5*PIMPL*PDFDT(KKB+1)/PDZZ(KKB+1)
+!!
+!! b(KKE) =  PRHODJ(KKE) / PTSTEP
+!!          -(PRHODJ(KKE)+PRHODJ(KKE-1))/2.*0.5*PIMPL*PDFDT(KKE)/PDZZ(KKE)
+!! a(KKE) = -(PRHODJ(KKE)+PRHODJ(KKE-1))/2.*0.5*PIMPL*PDFDT(KKE)/PDZZ(KKE)
+!!
+!!
+!!     EXTERNAL
+!!     --------
+!!
+!!       NONE
+!!
+!!     IMPLICIT ARGUMENTS
+!!     ------------------
+!!
+!!     REFERENCE
+!!     ---------
+!!       Press et al: Numerical recipes (1986) Cambridge Univ. Press
+!!
+!!     AUTHOR
+!!     ------
+!!       V. Masson and S. Malardel         * Meteo-France *   
+!! 
+!!     MODIFICATIONS
+!!     -------------
+!!       Original        07/2006
+!!       V.Masson : Optimization
+!!       S. Riette Jan 2012: support for both order of vertical levels
+!! ---------------------------------------------------------------------
+!
+!*       0. DECLARATIONS
+!
+USE MODD_PARAMETERS, ONLY: JPVEXT
+USE MODI_SHUMAN_MF
+!
+IMPLICIT NONE
+!
+!
+!*       0.1 declarations of arguments
+!
+INTEGER,                INTENT(IN)   :: KKA          ! near ground array index
+INTEGER,                INTENT(IN)   :: KKB          ! near ground physical index
+INTEGER,                INTENT(IN)   :: KKE          ! uppest atmosphere physical index
+INTEGER,                INTENT(IN)   :: KKU          ! uppest atmosphere array index
+INTEGER,                INTENT(IN)   :: KKL          ! +1 if grid goes from ground to atmosphere top, -1 otherwise
+REAL, DIMENSION(:,:), INTENT(IN) :: PVARM   ! variable at t-1      at mass point
+REAL, DIMENSION(:,:), INTENT(IN) :: PF      ! flux in dT/dt=-dF/dz at flux point
+REAL, DIMENSION(:,:), INTENT(IN) :: PDFDT   ! dF/dT                at flux point
+REAL,                   INTENT(IN) :: PTSTEP  ! Double time step
+REAL,                   INTENT(IN) :: PIMPL   ! implicit weight
+REAL, DIMENSION(:,:), INTENT(IN) :: PDZZ    ! Dz                   at flux point
+REAL, DIMENSION(:,:), INTENT(IN) :: PRHODJ  ! (dry rho)*J          at mass point
+!
+REAL, DIMENSION(:,:), INTENT(OUT):: PVARP   ! variable at t+1      at mass point
+!
+!
+!*       0.2 declarations of local variables
+!
+REAL, DIMENSION(SIZE(PVARM,1),SIZE(PVARM,2))  :: ZRHODJ_DFDT_O_DZ
+REAL, DIMENSION(SIZE(PVARM,1),SIZE(PVARM,2))  :: ZMZM_RHODJ
+REAL, DIMENSION(SIZE(PVARM,1),SIZE(PVARM,2))  :: ZA, ZB, ZC
+REAL, DIMENSION(SIZE(PVARM,1),SIZE(PVARM,2))  :: ZY ,ZGAM 
+                                         ! RHS of the equation, 3D work array
+REAL, DIMENSION(SIZE(PVARM,1))                :: ZBET
+                                         ! 2D work array
+INTEGER                              :: JK            ! loop counter
+!
+! ---------------------------------------------------------------------------
+!                                              
+!*      1.  Preliminaries
+!           -------------
+!
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('TRIDIAG_MASSFLUX',0,ZHOOK_HANDLE)
+ZMZM_RHODJ  = MZM_MF(KKA, KKU, KKL, PRHODJ)
+ZRHODJ_DFDT_O_DZ = ZMZM_RHODJ*PDFDT/PDZZ
+!
+ZA=0.
+ZB=0.
+ZC=0.
+ZY=0.
+!
+!
+!*      2.  COMPUTE THE RIGHT HAND SIDE
+!           ---------------------------
+!
+ZY(:,KKB) = PRHODJ(:,KKB)*PVARM(:,KKB)/PTSTEP             &
+    - ZMZM_RHODJ(:,KKB+KKL) * PF(:,KKB+KKL)/PDZZ(:,KKB+KKL)     &
+    + ZMZM_RHODJ(:,KKB  ) * PF(:,KKB  )/PDZZ(:,KKB  )     &
+    + ZRHODJ_DFDT_O_DZ(:,KKB+KKL) * 0.5*PIMPL * PVARM(:,KKB+KKL)    &
+    + ZRHODJ_DFDT_O_DZ(:,KKB+KKL) * 0.5*PIMPL * PVARM(:,KKB  )
+!
+DO JK=2+JPVEXT,SIZE(ZY,2)-JPVEXT-1
+  ZY(:,JK) = PRHODJ(:,JK)*PVARM(:,JK)/PTSTEP          &
+    - ZMZM_RHODJ(:,JK+KKL) * PF(:,JK+KKL)/PDZZ(:,JK+KKL)    &
+    + ZMZM_RHODJ(:,JK  ) * PF(:,JK  )/PDZZ(:,JK  )    &
+    + ZRHODJ_DFDT_O_DZ(:,JK+KKL) * 0.5*PIMPL * PVARM(:,JK+KKL)  &
+    + ZRHODJ_DFDT_O_DZ(:,JK+KKL) * 0.5*PIMPL * PVARM(:,JK  )  &
+    - ZRHODJ_DFDT_O_DZ(:,JK  ) * 0.5*PIMPL * PVARM(:,JK  )  &
+    - ZRHODJ_DFDT_O_DZ(:,JK  ) * 0.5*PIMPL * PVARM(:,JK-KKL)
+END DO
+! 
+IF (JPVEXT==0) THEN
+  ZY(:,KKE) = PRHODJ(:,KKE)*PVARM(:,KKE)/PTSTEP 
+ELSE
+  ZY(:,KKE) = PRHODJ(:,KKE)*PVARM(:,KKE)/PTSTEP &
+   - ZMZM_RHODJ(:,KKE+KKL) * PF(:,KKE+KKL)/PDZZ(:,KKE+KKL) &
+   + ZMZM_RHODJ(:,KKE  ) * PF(:,KKE  )/PDZZ(:,KKE  ) &
+   - ZRHODJ_DFDT_O_DZ(:,KKE ) * 0.5*PIMPL * PVARM(:,KKE  ) &
+   - ZRHODJ_DFDT_O_DZ(:,KKE ) * 0.5*PIMPL * PVARM(:,KKE-KKL)
+ENDIF
+!
+!
+!*       3.  INVERSION OF THE TRIDIAGONAL SYSTEM
+!            -----------------------------------
+!
+IF ( PIMPL > 1.E-10 ) THEN
+!
+!*       3.1 arrays A, B, C
+!            --------------
+!
+  ZB(:,KKB) =   PRHODJ(:,KKB)/PTSTEP                   &
+                + ZRHODJ_DFDT_O_DZ(:,KKB+KKL) * 0.5*PIMPL
+  ZC(:,KKB) =   ZRHODJ_DFDT_O_DZ(:,KKB+KKL) * 0.5*PIMPL
+
+  DO JK=2+JPVEXT,SIZE(ZY,2)-JPVEXT-1
+    ZA(:,JK) = - ZRHODJ_DFDT_O_DZ(:,JK  ) * 0.5*PIMPL
+    ZB(:,JK) =   PRHODJ(:,JK)/PTSTEP                   &
+                 + ZRHODJ_DFDT_O_DZ(:,JK+KKL) * 0.5*PIMPL &
+                 - ZRHODJ_DFDT_O_DZ(:,JK  ) * 0.5*PIMPL
+    ZC(:,JK) =   ZRHODJ_DFDT_O_DZ(:,JK+KKL) * 0.5*PIMPL
+  END DO
+
+  ZA(:,KKE) = - ZRHODJ_DFDT_O_DZ(:,KKE  ) * 0.5*PIMPL
+  ZB(:,KKE) =   PRHODJ(:,KKE)/PTSTEP                   &
+                - ZRHODJ_DFDT_O_DZ(:,KKE  ) * 0.5*PIMPL
+!
+!*       3.2 going up
+!            --------
+!
+  ZBET(:) = ZB(:,KKB)  ! bet = b(KKB)
+  PVARP(:,KKB) = ZY(:,KKB) / ZBET(:)
+
+  !
+  DO JK = KKB+KKL,KKE-KKL,KKL
+    ZGAM(:,JK) = ZC(:,JK-KKL) / ZBET(:)
+                                                    ! gam(k) = c(k-1) / bet
+    ZBET(:)    = ZB(:,JK) - ZA(:,JK) * ZGAM(:,JK)
+                                                    ! bet = b(k) - a(k)* gam(k)  
+    PVARP(:,JK)= ( ZY(:,JK) - ZA(:,JK) * PVARP(:,JK-KKL) ) / ZBET(:)
+                                        ! res(k) = (y(k) -a(k)*res(k-1))/ bet 
+  END DO 
+  ! special treatment for the last level
+  ZGAM(:,KKE) = ZC(:,KKE-KKL) / ZBET(:)
+                                                    ! gam(k) = c(k-1) / bet
+  ZBET(:)     = ZB(:,KKE) - ZA(:,KKE) * ZGAM(:,KKE)
+                                                    ! bet = b(k) - a(k)* gam(k)  
+  PVARP(:,KKE)= ( ZY(:,KKE) - ZA(:,KKE) * PVARP(:,KKE-KKL) ) / ZBET(:)
+                                       ! res(k) = (y(k) -a(k)*res(k-1))/ bet 
+!
+!*       3.3 going down
+!            ----------
+!
+  DO JK = KKE-KKL,KKB,-KKL
+    PVARP(:,JK) = PVARP(:,JK) - ZGAM(:,JK+KKL) * PVARP(:,JK+KKL)
+  END DO
+!
+!
+ELSE
+  !!! EXPLICIT FORMULATION
+  !
+  DO JK=1+JPVEXT,SIZE(PVARP,2)-JPVEXT
+    PVARP(:,JK) = ZY(:,JK) * PTSTEP / PRHODJ(:,JK)
+  ENDDO
+  !
+END IF 
+!
+!
+!*       4.  FILL THE UPPER AND LOWER EXTERNAL VALUES
+!            ----------------------------------------
+!
+PVARP(:,KKA)=PVARP(:,KKB)
+PVARP(:,KKU)=PVARP(:,KKE)
+!
+!-------------------------------------------------------------------------------
+!
+IF (LHOOK) CALL DR_HOOK('TRIDIAG_MASSFLUX',1,ZHOOK_HANDLE)
+END SUBROUTINE TRIDIAG_MASSFLUX
+END MODULE MODE_TRIDIAG_MASSFLUX
diff --git a/src/mesonh/turb/mode_turb_ver_dyn_flux.f90 b/src/mesonh/turb/mode_turb_ver_dyn_flux.f90
index 6caeab9b0..4e01503a6 100644
--- a/src/mesonh/turb/mode_turb_ver_dyn_flux.f90
+++ b/src/mesonh/turb/mode_turb_ver_dyn_flux.f90
@@ -115,9 +115,6 @@ SUBROUTINE TURB_VER_DYN_FLUX(KKA,KKU,KKL,                     &
 !!      DXF,DYF,DZF,DZM
 !!                             :  Shuman functions (difference operators)     
 !!                               
-!!      SUBROUTINE TRIDIAG     : to compute the split implicit evolution
-!!                               of a variable located at a mass point
-!!
 !!      SUBROUTINE TRIDIAG_WIND: to compute the split implicit evolution
 !!                               of a variable located at a wind point
 !!
@@ -230,8 +227,7 @@ USE MODI_GRADIENT_M
 USE MODI_SECOND_MNH
 USE MODI_SHUMAN , ONLY: MZM, MZF, MXM, MXF, MYM, MYF,&
                       & DZM, DXF, DXM, DYF, DYM
-USE MODI_TRIDIAG 
-USE MODI_TRIDIAG_WIND 
+USE MODE_TRIDIAG_WIND, ONLY: TRIDIAG_WIND
 USE MODI_LES_MEAN_SUBGRID
 !
 USE MODE_IO_FIELD_WRITE, only: IO_Field_write
@@ -331,16 +327,16 @@ INTEGER             :: IIU,IJU      ! size of array in x,y,z directions
 !
 REAL :: ZTIME1, ZTIME2, ZCMFS
 TYPE(TFIELDDATA) :: TZFIELD
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
 !----------------------------------------------------------------------------
 !
 !*       1.   PRELIMINARIES
 !             -------------
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('TURB_VER_DYN_FLUX',0,ZHOOK_HANDLE)
+!
 ZA=XUNDEF
 PDP=XUNDEF
 !
-IF (LHOOK) CALL DR_HOOK('TURB_VER_DYN_FLUX',0,ZHOOK_HANDLE)
-!
 IIU=SIZE(PUM,1)
 IJU=SIZE(PUM,2)
 CALL GET_INDICE_ll (IIB,IJB,IIE,IJE,IIU,IJU)
@@ -351,8 +347,6 @@ IKE=KKU-JPVEXT_TURB*KKL
 IKT=SIZE(PUM,3)          
 IKTB=1+JPVEXT_TURB              
 IKTE=IKT-JPVEXT_TURB
-
-
 !
 ZSOURCE = 0.
 ZFLXZ   = 0.
diff --git a/src/mesonh/turb/mode_turb_ver_sv_flux.f90 b/src/mesonh/turb/mode_turb_ver_sv_flux.f90
index e9ae7372c..624f9b12b 100644
--- a/src/mesonh/turb/mode_turb_ver_sv_flux.f90
+++ b/src/mesonh/turb/mode_turb_ver_sv_flux.f90
@@ -226,8 +226,7 @@ USE MODI_GRADIENT_V
 USE MODI_GRADIENT_W
 USE MODI_GRADIENT_M
 USE MODI_SHUMAN , ONLY : DZM, MZM, MZF
-USE MODI_TRIDIAG 
-USE MODI_TRIDIAG_WIND 
+USE MODE_TRIDIAG, ONLY: TRIDIAG
 USE MODI_EMOIST
 USE MODI_ETHETA
 USE MODI_LES_MEAN_SUBGRID
diff --git a/src/mesonh/turb/mode_turb_ver_thermo_corr.f90 b/src/mesonh/turb/mode_turb_ver_thermo_corr.f90
index 7913b4e6a..db08f0ce6 100644
--- a/src/mesonh/turb/mode_turb_ver_thermo_corr.f90
+++ b/src/mesonh/turb/mode_turb_ver_thermo_corr.f90
@@ -117,12 +117,6 @@ SUBROUTINE TURB_VER_THERMO_CORR(KKA,KKU,KKL,KRR,KRRL,KRRI,    &
 !!      DXF,DYF,DZF,DZM
 !!                             :  Shuman functions (difference operators)     
 !!                               
-!!      SUBROUTINE TRIDIAG     : to compute the split implicit evolution
-!!                               of a variable located at a mass point
-!!
-!!      SUBROUTINE TRIDIAG_WIND: to compute the split implicit evolution
-!!                               of a variable located at a wind point
-!!
 !!      FUNCTIONs ETHETA and EMOIST  :  
 !!            allows to compute:
 !!            - the coefficients for the turbulent correlation between
@@ -223,9 +217,7 @@ USE MODI_GRADIENT_V
 USE MODI_GRADIENT_W
 USE MODI_GRADIENT_M
 USE MODI_SHUMAN , ONLY : DZM, MZM, MZF
-USE MODI_TRIDIAG 
 USE MODI_LES_MEAN_SUBGRID
-USE MODI_TRIDIAG_THERMO
 !
 USE MODE_IO_FIELD_WRITE, only: IO_Field_write
 USE MODE_PRANDTL
diff --git a/src/mesonh/turb/mode_turb_ver_thermo_flux.f90 b/src/mesonh/turb/mode_turb_ver_thermo_flux.f90
index 78a188588..6d0f1a009 100644
--- a/src/mesonh/turb/mode_turb_ver_thermo_flux.f90
+++ b/src/mesonh/turb/mode_turb_ver_thermo_flux.f90
@@ -250,10 +250,10 @@ USE MODI_GRADIENT_V
 USE MODI_GRADIENT_W
 USE MODI_GRADIENT_M
 USE MODI_SHUMAN , ONLY : DZF, DZM, MZF, MZM, MYF, MXF
-USE MODI_TRIDIAG 
+USE MODE_TRIDIAG, ONLY: TRIDIAG
 USE MODI_LES_MEAN_SUBGRID
-USE MODI_TRIDIAG_THERMO
-USE MODI_TM06_H
+USE MODE_TRIDIAG_THERMO, ONLY: TRIDIAG_THERMO
+USE MODE_TM06_H, ONLY: TM06_H
 !
 USE MODE_IO_FIELD_WRITE, ONLY: IO_FIELD_WRITE
 USE MODE_PRANDTL
diff --git a/src/mesonh/turb/turb.f90 b/src/mesonh/turb/turb.f90
index c255d2dd5..f1b3f9e11 100644
--- a/src/mesonh/turb/turb.f90
+++ b/src/mesonh/turb/turb.f90
@@ -360,6 +360,9 @@ END MODULE MODI_TURB
 !*      0. DECLARATIONS
 !          ------------
 !
+USE PARKIND1, ONLY : JPRB
+USE YOMHOOK , ONLY : LHOOK, DR_HOOK
+!
 use modd_budget,      only: lbudget_u,  lbudget_v,  lbudget_w,  lbudget_th, lbudget_rv, lbudget_rc,  &
                             lbudget_rr, lbudget_ri, lbudget_rs, lbudget_rg, lbudget_rh, lbudget_sv,  &
                             NBUDGET_U,  NBUDGET_V,  NBUDGET_W,  NBUDGET_TH, NBUDGET_RV, NBUDGET_RC,  &
@@ -369,7 +372,7 @@ USE MODD_CONF
 USE MODD_CST
 USE MODD_CTURB
 USE MODD_DYN_n, ONLY : LOCEAN
-use modd_field,          only: tfielddata, TYPEREAL
+USE MODD_FIELD, ONLY: TFIELDDATA,TYPEREAL
 USE MODD_IO, ONLY: TFILEDATA
 USE MODD_LES
 USE MODD_NSV
@@ -390,7 +393,7 @@ USE MODI_GRADIENT_M
 USE MODI_LES_MEAN_SUBGRID
 USE MODI_RMC01
 USE MODI_GRADIENT_W
-USE MODI_TM06
+USE MODE_TM06, ONLY: TM06
 USE MODI_UPDATE_LM
 USE MODI_GET_HALO
 !
@@ -1341,15 +1344,19 @@ REAL, DIMENSION(:,:), INTENT(INOUT) :: PUSLOPE,PVSLOPE
 !
 !*       0.2   Declarations of local variables :
 !
-INTEGER             :: IIB,IIE,IJB,IJE ! index values for the physical subdomain
+INTEGER             :: IIB,IIE,IJB,IJE,IIU,IJU ! index values for the physical subdomain
 TYPE(LIST_ll), POINTER :: TZFIELDS_ll  ! list of fields to exchange
 INTEGER                :: IINFO_ll     ! return code of parallel routine
+REAL(KIND=JPRB) :: ZHOOK_HANDLE
+IF (LHOOK) CALL DR_HOOK('TURB:UPDATE_ROTATE_WIND',0,ZHOOK_HANDLE)
 !
 !*        1  PROLOGUE
 !
 NULLIFY(TZFIELDS_ll)
 !
-CALL GET_INDICE_ll (IIB,IJB,IIE,IJE)
+IIU=SIZE(PUSLOPE,1)
+IJU=SIZE(PUSLOPE,2)
+CALL GET_INDICE_ll (IIB,IJB,IIE,IJE,IIU,IJU)
 !
 !         2 Update halo if necessary
 !
@@ -1379,6 +1386,8 @@ IF(  HLBCY(2) /= "CYCL" .AND. LNORTH_ll()) THEN
   PVSLOPE(:,IJE+1)=PVSLOPE(:,IJE)
 END IF
 !
+IF (LHOOK) CALL DR_HOOK('TURB:UPDATE_ROTATE_WIND',1,ZHOOK_HANDLE)
+!
 END SUBROUTINE UPDATE_ROTATE_WIND
 !
 !     ########################################################################
@@ -1421,6 +1430,8 @@ REAL, DIMENSION(SIZE(PEXN,1),SIZE(PEXN,2),SIZE(PEXN,3)) :: ZDRVSATDT
 !
 !-------------------------------------------------------------------------------
 !
+  REAL(KIND=JPRB) :: ZHOOK_HANDLE
+  IF (LHOOK) CALL DR_HOOK('TURB:COMPUTE_FUNCTION_THERMO',0,ZHOOK_HANDLE)
   ZEPS = XMV / XMD
 !
 !*       1.1 Lv/Cph at  t
@@ -1462,6 +1473,7 @@ REAL, DIMENSION(SIZE(PEXN,1),SIZE(PEXN,2),SIZE(PEXN,3)) :: ZDRVSATDT
 !
   PLOCPEXN(:,:,:) = PLOCPEXN(:,:,:) / PEXN(:,:,:)
 !
+IF (LHOOK) CALL DR_HOOK('TURB:COMPUTE_FUNCTION_THERMO',1,ZHOOK_HANDLE)
 END SUBROUTINE COMPUTE_FUNCTION_THERMO
 !
 !     ####################
-- 
GitLab