From 896d62dac92842df76384277f1e6bda3f7a42c1f Mon Sep 17 00:00:00 2001
From: Quentin Rodier <quentin.rodier@meteo.fr>
Date: Mon, 4 Apr 2022 18:36:50 +0200
Subject: [PATCH] Quentin 04/04/2022: Expand Array tridiag_* except
 tridiag_massflux

---
 src/common/turb/mode_tridiag.F90        | 36 +++++++++++++----
 src/common/turb/mode_tridiag_thermo.F90 | 52 +++++++++++++++++--------
 src/common/turb/mode_tridiag_tke.F90    | 30 +++++++++++---
 src/common/turb/mode_tridiag_wind.F90   | 47 +++++++++++++++-------
 4 files changed, 123 insertions(+), 42 deletions(-)

diff --git a/src/common/turb/mode_tridiag.F90 b/src/common/turb/mode_tridiag.F90
index 960c8b241..eba6cbdb0 100644
--- a/src/common/turb/mode_tridiag.F90
+++ b/src/common/turb/mode_tridiag.F90
@@ -138,7 +138,7 @@ REAL, DIMENSION(D%NIT,D%NJT,D%NKT)  :: ZY ,ZGAM
                                          ! RHS of the equation, 3D work array
 REAL, DIMENSION(D%NIT,D%NJT)         :: ZBET
                                          ! 2D work array
-INTEGER                              :: JK            ! loop counter
+INTEGER                              :: JI,JJ,JK            ! loop counter
 INTEGER                              :: IKB,IKE       ! inner vertical limits
 INTEGER                              :: IKT           ! array size in k direction
 INTEGER                              :: IKTB,IKTE     ! start, end of k loops in physical domain 
@@ -152,26 +152,32 @@ INTEGER                              :: IKTB,IKTE     ! start, end of k loops in
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 IF (LHOOK) CALL DR_HOOK('TRIDIAG',0,ZHOOK_HANDLE)
 !
-IKT=D%NKT  
-IKTB=D%NKTB          
+IKT=D%NKT
+IKTB=D%NKTB
 IKTE=D%NKTE
 IKB=D%NKB
-IKE=D%NKE  
+IKE=D%NKE
 !
+!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
 ZY(:,:,IKB) = PVARM(:,:,IKB)  + PTSTEP*PSOURCE(:,:,IKB) -   &
   PEXPL / PRHODJ(:,:,IKB) * PA(:,:,IKB+D%NKL) * (PVARM(:,:,IKB+D%NKL) - PVARM(:,:,IKB))
+!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
 !
 DO JK=IKTB+1,IKTE-1
+  !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
   ZY(:,:,JK)= PVARM(:,:,JK)  + PTSTEP*PSOURCE(:,:,JK) -               &
       PEXPL / PRHODJ(:,:,JK) *                                           &
                              ( PVARM(:,:,JK-D%NKL)*PA(:,:,JK)                &
                               -PVARM(:,:,JK)*(PA(:,:,JK)+PA(:,:,JK+D%NKL))   &
                               +PVARM(:,:,JK+D%NKL)*PA(:,:,JK+D%NKL)              &
                              ) 
+  !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
 END DO
 ! 
+!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
 ZY(:,:,IKE)= PVARM(:,:,IKE) + PTSTEP*PSOURCE(:,:,IKE) +               &
   PEXPL / PRHODJ(:,:,IKE) * PA(:,:,IKE) * (PVARM(:,:,IKE)-PVARM(:,:,IKE-D%NKL))
+!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
 !
 !
 !*       2.  INVERSION OF THE TRIDIAGONAL SYSTEM
@@ -182,11 +188,14 @@ IF ( PIMPL > 1.E-10 ) THEN
   !
   !  going up
   !
+  !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
   ZBET(:,:) = 1. - PIMPL * PA(:,:,IKB+D%NKL) / PRHODJ(:,:,IKB)  ! bet = b(ikb)
-  PVARP(:,:,IKB) = ZY(:,:,IKB) / ZBET(:,:)                
+  PVARP(:,:,IKB) = ZY(:,:,IKB) / ZBET(:,:)
+  !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)               
   !
   DO JK = IKB+D%NKL,IKE-D%NKL,D%NKL
-      ZGAM(:,:,JK) = PIMPL * PA(:,:,JK) / PRHODJ(:,:,JK-D%NKL) / ZBET(:,:)  
+    !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+    ZGAM(:,:,JK) = PIMPL * PA(:,:,JK) / PRHODJ(:,:,JK-D%NKL) / ZBET(:,:)  
                                                     ! gam(k) = c(k-1) / bet
     ZBET(:,:)    = 1. - PIMPL * (  PA(:,:,JK) * (1. + ZGAM(:,:,JK))  &
                                  + PA(:,:,JK+D%NKL)                      &
@@ -196,7 +205,9 @@ IF ( PIMPL > 1.E-10 ) THEN
                     * PVARP(:,:,JK-D%NKL)                                 &
                    ) / ZBET(:,:)
                                         ! res(k) = (y(k) -a(k)*res(k-1))/ bet 
-  END DO 
+    !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+  END DO
+  !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
   ! special treatment for the last level
   ZGAM(:,:,IKE) = PIMPL * PA(:,:,IKE) / PRHODJ(:,:,IKE-D%NKL) / ZBET(:,:) 
                                                     ! gam(k) = c(k-1) / bet
@@ -210,13 +221,20 @@ IF ( PIMPL > 1.E-10 ) THEN
   !
   !  going down
   !
+  !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
   DO JK = IKE-D%NKL,IKB,-1*D%NKL
+  !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
     PVARP(:,:,JK) = PVARP(:,:,JK) - ZGAM(:,:,JK+D%NKL) * PVARP(:,:,JK+D%NKL) 
+  !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
   END DO
 !
 ELSE
 ! 
-  PVARP(:,:,IKTB:IKTE) = ZY(:,:,IKTB:IKTE)
+  !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+  DO JK=IKTB,IKTE
+    PVARP(:,:,JK) = ZY(:,:,JK)
+  END DO
+  !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
 !
 END IF 
 !
@@ -224,8 +242,10 @@ END IF
 !*       3.  FILL THE UPPER AND LOWER EXTERNAL VALUES
 !            ----------------------------------------
 !
+!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
 PVARP(:,:,D%NKA)=PVARP(:,:,IKB)
 PVARP(:,:,D%NKU)=PVARP(:,:,IKE)
+!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
 !
 !-------------------------------------------------------------------------------
 !
diff --git a/src/common/turb/mode_tridiag_thermo.F90 b/src/common/turb/mode_tridiag_thermo.F90
index c87a9c702..a16d374ec 100644
--- a/src/common/turb/mode_tridiag_thermo.F90
+++ b/src/common/turb/mode_tridiag_thermo.F90
@@ -149,7 +149,7 @@ REAL, DIMENSION(D%NIT,D%NJT,D%NKT)  :: ZY ,ZGAM
                                          ! RHS of the equation, 3D work array
 REAL, DIMENSION(D%NIT,D%NJT)                :: ZBET
                                          ! 2D work array
-INTEGER             :: JK            ! loop counter
+INTEGER             :: JI,JJ,JK            ! loop counter
 INTEGER             :: IKB,IKE       ! inner vertical limits
 INTEGER             :: IKT          ! array size in k direction
 INTEGER             :: IKTB,IKTE    ! start, end of k loops in physical domain 
@@ -169,7 +169,9 @@ IKE=D%NKE
 
 !
 ZMZM_RHODJ  = MZM(PRHODJ,D%NKA,D%NKU,D%NKL)
-ZRHODJ_DFDDTDZ_O_DZ2 = ZMZM_RHODJ*PDFDDTDZ/PDZZ**2
+!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
+ZRHODJ_DFDDTDZ_O_DZ2(:,:,:) = ZMZM_RHODJ(:,:,:)*PDFDDTDZ(:,:,:)/PDZZ(:,:,:)**2
+!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
 !
 ZA=0.
 ZB=0.
@@ -180,25 +182,33 @@ ZY=0.
 !*      2.  COMPUTE THE RIGHT HAND SIDE
 !           ---------------------------
 !
+!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
 ZY(:,:,IKB) = PRHODJ(:,:,IKB)*PVARM(:,:,IKB)/PTSTEP                  &
     - ZMZM_RHODJ(:,:,IKB+D%NKL) * PF(:,:,IKB+D%NKL)/PDZZ(:,:,IKB+D%NKL)    &
     + ZMZM_RHODJ(:,:,IKB  ) * PF(:,:,IKB  )/PDZZ(:,:,IKB  )          &
     + ZRHODJ_DFDDTDZ_O_DZ2(:,:,IKB+D%NKL) * PIMPL * PVARM(:,:,IKB+D%NKL) &
     - ZRHODJ_DFDDTDZ_O_DZ2(:,:,IKB+D%NKL) * PIMPL * PVARM(:,:,IKB  )
-!
-  ZY(:,:,IKTB+1:IKTE-1) = PRHODJ(:,:,IKTB+1:IKTE-1)*PVARM(:,:,IKTB+1:IKTE-1)/PTSTEP                 &
-    - ZMZM_RHODJ(:,:,IKTB+1+D%NKL:IKTE-1+D%NKL) * PF(:,:,IKTB+1+D%NKL:IKTE-1+D%NKL)/PDZZ(:,:,IKTB+1+D%NKL:IKTE-1+D%NKL)     &
-    + ZMZM_RHODJ(:,:,IKTB+1:IKTE-1  ) * PF(:,:,IKTB+1:IKTE-1  )/PDZZ(:,:,IKTB+1:IKTE-1  )           &
-    + ZRHODJ_DFDDTDZ_O_DZ2(:,:,IKTB+1+D%NKL:IKTE-1+D%NKL) * PIMPL * PVARM(:,:,IKTB+1+D%NKL:IKTE-1+D%NKL) &
-    - ZRHODJ_DFDDTDZ_O_DZ2(:,:,IKTB+1+D%NKL:IKTE-1+D%NKL) * PIMPL * PVARM(:,:,IKTB+1:IKTE-1  )   &
-    - ZRHODJ_DFDDTDZ_O_DZ2(:,:,IKTB+1:IKTE-1    ) * PIMPL * PVARM(:,:,IKTB+1:IKTE-1  )   &
-    + ZRHODJ_DFDDTDZ_O_DZ2(:,:,IKTB+1:IKTE-1    ) * PIMPL * PVARM(:,:,IKTB+1-D%NKL:IKTE-1-D%NKL)
+!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+!
+!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+DO JK=IKTB+1,IKTE-1
+  ZY(:,:,JK) = PRHODJ(:,:,JK)*PVARM(:,:,JK)/PTSTEP                 &
+    - ZMZM_RHODJ(:,:,JK+D%NKL) * PF(:,:,JK+D%NKL)/PDZZ(:,:,JK+D%NKL)     &
+    + ZMZM_RHODJ(:,:,JK  ) * PF(:,:,JK  )/PDZZ(:,:,JK  )           &
+    + ZRHODJ_DFDDTDZ_O_DZ2(:,:,JK+D%NKL) * PIMPL * PVARM(:,:,JK+D%NKL) &
+    - ZRHODJ_DFDDTDZ_O_DZ2(:,:,JK+D%NKL) * PIMPL * PVARM(:,:,JK  )   &
+    - ZRHODJ_DFDDTDZ_O_DZ2(:,:,JK    ) * PIMPL * PVARM(:,:,JK  )   &
+    + ZRHODJ_DFDDTDZ_O_DZ2(:,:,JK    ) * PIMPL * PVARM(:,:,JK-D%NKL)
+END DO
+!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
 ! 
+!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
 ZY(:,:,IKE) = PRHODJ(:,:,IKE)*PVARM(:,:,IKE)/PTSTEP               &
     - ZMZM_RHODJ(:,:,IKE+D%NKL) * PF(:,:,IKE+D%NKL)/PDZZ(:,:,IKE+D%NKL) &
     + ZMZM_RHODJ(:,:,IKE  ) * PF(:,:,IKE  )/PDZZ(:,:,IKE  )       &
     - ZRHODJ_DFDDTDZ_O_DZ2(:,:,IKE ) * PIMPL * PVARM(:,:,IKE  )   &
     + ZRHODJ_DFDDTDZ_O_DZ2(:,:,IKE ) * PIMPL * PVARM(:,:,IKE-D%NKL)
+!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
 !
 !
 !*       3.  INVERSION OF THE TRIDIAGONAL SYSTEM
@@ -209,15 +219,18 @@ IF ( PIMPL > 1.E-10 ) THEN
 !*       3.1 arrays A, B, C
 !            --------------
 !
+  !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
   ZB(:,:,IKB) =   PRHODJ(:,:,IKB)/PTSTEP                   &
                 - ZRHODJ_DFDDTDZ_O_DZ2(:,:,IKB+D%NKL) * PIMPL
   ZC(:,:,IKB) =   ZRHODJ_DFDDTDZ_O_DZ2(:,:,IKB+D%NKL) * PIMPL
 !
-  ZA(:,:,IKTB+1:IKTE-1) =   ZRHODJ_DFDDTDZ_O_DZ2(:,:,IKTB+1:IKTE-1) * PIMPL
-  ZB(:,:,IKTB+1:IKTE-1) =   PRHODJ(:,:,IKTB+1:IKTE-1)/PTSTEP                        &
-                          - ZRHODJ_DFDDTDZ_O_DZ2(:,:,IKTB+1+D%NKL:IKTE-1+D%NKL) * PIMPL &
-                          - ZRHODJ_DFDDTDZ_O_DZ2(:,:,IKTB+1:IKTE-1) * PIMPL
-  ZC(:,:,IKTB+1:IKTE-1) =   ZRHODJ_DFDDTDZ_O_DZ2(:,:,IKTB+1+D%NKL:IKTE-1+D%NKL) * PIMPL
+  DO JK=IKTB+1,IKTE-1
+    ZA(:,:,JK) =   ZRHODJ_DFDDTDZ_O_DZ2(:,:,JK) * PIMPL
+    ZB(:,:,JK) =   PRHODJ(:,:,JK)/PTSTEP                        &
+                            - ZRHODJ_DFDDTDZ_O_DZ2(:,:,JK+D%NKL) * PIMPL &
+                            - ZRHODJ_DFDDTDZ_O_DZ2(:,:,JK) * PIMPL
+    ZC(:,:,JK) =   ZRHODJ_DFDDTDZ_O_DZ2(:,:,JK+D%NKL) * PIMPL
+  END DO
 !
   ZA(:,:,IKE) =   ZRHODJ_DFDDTDZ_O_DZ2(:,:,IKE  ) * PIMPL
   ZB(:,:,IKE) =   PRHODJ(:,:,IKE)/PTSTEP                   &
@@ -253,9 +266,14 @@ IF ( PIMPL > 1.E-10 ) THEN
     PVARP(:,:,JK) = PVARP(:,:,JK) - ZGAM(:,:,JK+D%NKL) * PVARP(:,:,JK+D%NKL)
   END DO
 !
+  !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
 ELSE
 ! 
-  PVARP(:,:,IKTB:IKTE) = ZY(:,:,IKTB:IKTE) * PTSTEP / PRHODJ(:,:,IKTB:IKTE)
+  !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+  DO JK=IKTB,IKTE
+    PVARP(:,:,JK) = ZY(:,:,JK) * PTSTEP / PRHODJ(:,:,JK)
+  END DO
+  !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
 !
 END IF 
 !
@@ -263,8 +281,10 @@ END IF
 !*       4.  FILL THE UPPER AND LOWER EXTERNAL VALUES
 !            ----------------------------------------
 !
+!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
 PVARP(:,:,D%NKA)=PVARP(:,:,IKB)
 PVARP(:,:,D%NKU)=PVARP(:,:,IKE)
+!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
 !
 !-------------------------------------------------------------------------------
 !
diff --git a/src/common/turb/mode_tridiag_tke.F90 b/src/common/turb/mode_tridiag_tke.F90
index 9c9fce379..81364c2f1 100644
--- a/src/common/turb/mode_tridiag_tke.F90
+++ b/src/common/turb/mode_tridiag_tke.F90
@@ -138,8 +138,8 @@ REAL, DIMENSION(D%NIT,D%NJT,D%NKT)  :: ZY ,ZGAM
                                          ! RHS of the equation, 3D work array
 REAL, DIMENSION(D%NIT,D%NJT)        :: ZBET
                                          ! 2D work array
-INTEGER             :: JK            ! loop counter
-INTEGER             :: IKB,IKE       ! inner vertical limits
+INTEGER             :: JI,JJ,JK     ! loop counter
+INTEGER             :: IKB,IKE      ! inner vertical limits
 INTEGER             :: IKT          ! array size in k direction
 INTEGER             :: IKTB,IKTE    ! start, end of k loops in physical domain
 !
@@ -157,20 +157,26 @@ IKTE=D%NKTE
 IKB=D%NKB
 IKE=D%NKE
 !
+!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
 ZY(:,:,IKB) = PVARM(:,:,IKB)  + PTSTEP*PSOURCE(:,:,IKB) -   &
   PEXPL / PRHODJ(:,:,IKB) * PA(:,:,IKB+D%NKL) * (PVARM(:,:,IKB+D%NKL) - PVARM(:,:,IKB))
+!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
 !
 DO JK=IKTB+1,IKTE-1
+  !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
   ZY(:,:,JK)= PVARM(:,:,JK)  + PTSTEP*PSOURCE(:,:,JK) -               &
       PEXPL / PRHODJ(:,:,JK) *                                           &
                              ( PVARM(:,:,JK-D%NKL)*PA(:,:,JK)                &
                               -PVARM(:,:,JK)*(PA(:,:,JK)+PA(:,:,JK+D%NKL))   &
                               +PVARM(:,:,JK+D%NKL)*PA(:,:,JK+D%NKL)              &
                              ) 
+  !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
 END DO
 ! 
+!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
 ZY(:,:,IKE)= PVARM(:,:,IKE) + PTSTEP*PSOURCE(:,:,IKE) +               &
   PEXPL / PRHODJ(:,:,IKE) * PA(:,:,IKE) * (PVARM(:,:,IKE)-PVARM(:,:,IKE-D%NKL))
+!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
 !
 !
 !*       2.  INVERSION OF THE TRIDIAGONAL SYSTEM
@@ -181,12 +187,15 @@ IF ( PIMPL > 1.E-10 ) THEN
   !
   !  going up
   !
+  !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
   ZBET(:,:) = 1. + PIMPL * (PDIAG(:,:,IKB)-PA(:,:,IKB+D%NKL) / PRHODJ(:,:,IKB))
                                                     ! bet = b(ikb)
   PVARP(:,:,IKB) = ZY(:,:,IKB) / ZBET(:,:)                
+  !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)               
   !
   DO JK = IKB+D%NKL,IKE-D%NKL,D%NKL
-      ZGAM(:,:,JK) = PIMPL * PA(:,:,JK) / PRHODJ(:,:,JK-D%NKL) / ZBET(:,:)  
+    !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+    ZGAM(:,:,JK) = PIMPL * PA(:,:,JK) / PRHODJ(:,:,JK-D%NKL) / ZBET(:,:)  
                                                     ! gam(k) = c(k-1) / bet
     ZBET(:,:)    = 1. + PIMPL * ( PDIAG(:,:,JK) -                     &
                                  (  PA(:,:,JK) * (1. + ZGAM(:,:,JK))  &
@@ -197,7 +206,9 @@ IF ( PIMPL > 1.E-10 ) THEN
                     * PVARP(:,:,JK-D%NKL)                                 &
                    ) / ZBET(:,:)
                                         ! res(k) = (y(k) -a(k)*res(k-1))/ bet 
-  END DO 
+    !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+  END DO
+  !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
   ! special treatment for the last level
   ZGAM(:,:,IKE) = PIMPL * PA(:,:,IKE) / PRHODJ(:,:,IKE-D%NKL) / ZBET(:,:) 
                                                     ! gam(k) = c(k-1) / bet
@@ -212,13 +223,20 @@ IF ( PIMPL > 1.E-10 ) THEN
   !
   !  going down
   !
+  !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
   DO JK = IKE-D%NKL,IKB,-1*D%NKL
+  !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
     PVARP(:,:,JK) = PVARP(:,:,JK) - ZGAM(:,:,JK+D%NKL) * PVARP(:,:,JK+D%NKL) 
+  !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
   END DO
 !
 ELSE
 ! 
-  PVARP(:,:,IKTB:IKTE) = ZY(:,:,IKTB:IKTE)
+  !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+  DO JK=IKTB,IKTE
+    PVARP(:,:,JK) = ZY(:,:,JK)
+  END DO
+  !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
 !
 END IF 
 !
@@ -226,8 +244,10 @@ END IF
 !*       3.  FILL THE UPPER AND LOWER EXTERNAL VALUES
 !            ----------------------------------------
 !
+!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
 PVARP(:,:,D%NKA)=PVARP(:,:,IKB)
 PVARP(:,:,D%NKU)=PVARP(:,:,IKE)
+!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
 !
 !-------------------------------------------------------------------------------
 !
diff --git a/src/common/turb/mode_tridiag_wind.F90 b/src/common/turb/mode_tridiag_wind.F90
index 9a13f1c35..ca73dd9e1 100644
--- a/src/common/turb/mode_tridiag_wind.F90
+++ b/src/common/turb/mode_tridiag_wind.F90
@@ -143,8 +143,8 @@ REAL, DIMENSION(SIZE(PVARM,1),SIZE(PVARM,2),SIZE(PVARM,3))  :: ZY ,ZGAM
                                          ! RHS of the equation, 3D work array
 REAL, DIMENSION(SIZE(PVARM,1),SIZE(PVARM,2))                :: ZBET
                                          ! 2D work array
-INTEGER             :: JK            ! loop counter
-INTEGER             :: IKB,IKE       ! inner vertical limits
+INTEGER             :: JI,JJ,JK     ! loop counter
+INTEGER             :: IKB,IKE      ! inner vertical limits
 INTEGER             :: IKT          ! array size in k direction
 INTEGER             :: IKTB,IKTE    ! start, end of k loops in physical domain 
 !
@@ -162,19 +162,26 @@ IKTE=D%NKTE
 IKB=D%NKB
 IKE=D%NKE
 !
-! 
+!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
 ZY(:,:,IKB) = PVARM(:,:,IKB)  + PTSTEP*PSOURCE(:,:,IKB) -   &
   PEXPL / PRHODJA(:,:,IKB) * PA(:,:,IKB+D%NKL) * (PVARM(:,:,IKB+D%NKL) - PVARM(:,:,IKB))
+!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
 !
-  ZY(:,:,IKTB+1:IKTE-1)= PVARM(:,:,IKTB+1:IKTE-1)  + PTSTEP*PSOURCE(:,:,IKTB+1:IKTE-1) -               &
-      PEXPL / PRHODJA(:,:,IKTB+1:IKTE-1) *                                          &
-                             ( PVARM(:,:,IKTB+1-D%NKL:IKTE-1-D%NKL)*PA(:,:,IKTB+1:IKTE-1)                &
-                              -PVARM(:,:,IKTB+1:IKTE-1)*(PA(:,:,IKTB+1:IKTE-1)+PA(:,:,IKTB+1+D%NKL:IKTE-1+D%NKL))   &
-                              +PVARM(:,:,IKTB+1+D%NKL:IKTE-1+D%NKL)*PA(:,:,IKTB+1+D%NKL:IKTE-1+D%NKL)              &
+DO JK=IKTB+1,IKTE-1
+  !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+  ZY(:,:,JK)= PVARM(:,:,JK)  + PTSTEP*PSOURCE(:,:,JK) -               &
+      PEXPL / PRHODJA(:,:,JK) *                                          &
+                             ( PVARM(:,:,JK-D%NKL)*PA(:,:,JK)                &
+                              -PVARM(:,:,JK)*(PA(:,:,JK)+PA(:,:,JK+D%NKL))   &
+                              +PVARM(:,:,JK+D%NKL)*PA(:,:,JK+D%NKL)              &
                              ) 
+  !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+END DO
 ! 
+!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
 ZY(:,:,IKE)= PVARM(:,:,IKE) + PTSTEP*PSOURCE(:,:,IKE) +               &
   PEXPL / PRHODJA(:,:,IKE) * PA(:,:,IKE) * (PVARM(:,:,IKE)-PVARM(:,:,IKE-D%NKL))
+!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
 !
 !
 !*       2.  INVERSION OF THE TRIDIAGONAL SYSTEM
@@ -185,11 +192,14 @@ IF ( PIMPL > 1.E-10 ) THEN
   !
   !  going up
   !
+  !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
   ZBET(:,:) = 1. - PIMPL * (  PA(:,:,IKB+D%NKL) / PRHODJA(:,:,IKB) &  
                             + PCOEFS(:,:) *  PTSTEP        )   ! bet = b(ikb)
-  PVARP(:,:,IKB) = ZY(:,:,IKB) / ZBET(:,:)                
+  PVARP(:,:,IKB) = ZY(:,:,IKB) / ZBET(:,:)
+  !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)               
   !
   DO JK = IKB+D%NKL,IKE-D%NKL,D%NKL
+    !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
     ZGAM(:,:,JK) = PIMPL * PA(:,:,JK) / PRHODJA(:,:,JK-D%NKL) / ZBET(:,:)  
                                                     ! gam(k) = c(k-1) / bet
     ZBET(:,:)    = 1. - PIMPL * (  PA(:,:,JK) * (1. + ZGAM(:,:,JK))  &
@@ -197,10 +207,12 @@ IF ( PIMPL > 1.E-10 ) THEN
                                 ) / PRHODJA(:,:,JK)  
                                                     ! bet = b(k) - a(k)* gam(k)  
     PVARP(:,:,JK)= ( ZY(:,:,JK) - PIMPL * PA(:,:,JK) / PRHODJA(:,:,JK) &
-                    * PVARP(:,:,JK-D%NKL)                                  &
+                    * PVARP(:,:,JK-D%NKL)                                 &
                    ) / ZBET(:,:)
                                         ! res(k) = (y(k) -a(k)*res(k-1))/ bet 
+    !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
   END DO
+  !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
   ! special treatment for the last level
   ZGAM(:,:,IKE) = PIMPL * PA(:,:,IKE) / PRHODJA(:,:,IKE-D%NKL) / ZBET(:,:) 
                                                     ! gam(k) = c(k-1) / bet
@@ -212,15 +224,22 @@ IF ( PIMPL > 1.E-10 ) THEN
                   ) / ZBET(:,:)
                                         ! res(k) = (y(k) -a(k)*res(k-1))/ bet 
   !
-  !  going down 
+  !  going down
   !
+  !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
   DO JK = IKE-D%NKL,IKB,-1*D%NKL
+  !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
     PVARP(:,:,JK) = PVARP(:,:,JK) - ZGAM(:,:,JK+D%NKL) * PVARP(:,:,JK+D%NKL) 
+  !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
   END DO
 !
 ELSE
 ! 
-  PVARP(:,:,IKTB:IKTE) = ZY(:,:,IKTB:IKTE)
+  !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+  DO JK=IKTB,IKTE
+    PVARP(:,:,JK) = ZY(:,:,JK)
+  END DO
+  !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
 !
 END IF 
 !
@@ -228,11 +247,13 @@ END IF
 !*       3.  FILL THE UPPER AND LOWER EXTERNAL VALUES
 !            ----------------------------------------
 !
+!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
 PVARP(:,:,D%NKA)=PVARP(:,:,IKB)
 PVARP(:,:,D%NKU)=PVARP(:,:,IKE)
+!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
 !
 !-------------------------------------------------------------------------------
-! 
+!
 IF (LHOOK) CALL DR_HOOK('TRIDIAG_WIND',1,ZHOOK_HANDLE)
 END SUBROUTINE TRIDIAG_WIND
 END MODULE MODE_TRIDIAG_WIND
-- 
GitLab