From dc9ef33218cc2507520e1ec64f93f500fd5e708d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9bastien=20Riette?= <sebastien.riette@meteo.fr>
Date: Thu, 14 Apr 2022 18:03:04 +0200
Subject: [PATCH] S. Riette 14/4/2022 bugfixes around mnh_expand

Execution time is now equivalent with or without array-syntax / do loops conversion
---
 src/common/turb/mode_prandtl.F90              |  8 ++++----
 src/common/turb/mode_tke_eps_sources.F90      |  2 +-
 src/common/turb/mode_tridiag_thermo.F90       | 16 +++++++++++++---
 src/common/turb/mode_tridiag_tke.F90          |  4 ++--
 src/common/turb/mode_tridiag_wind.F90         |  4 ++--
 src/common/turb/mode_turb_ver.F90             |  6 +++---
 src/common/turb/mode_turb_ver_sv_flux.F90     |  4 ++++
 src/common/turb/mode_turb_ver_thermo_corr.F90 |  2 +-
 src/common/turb/mode_turb_ver_thermo_flux.F90 | 16 +++++++++++-----
 9 files changed, 41 insertions(+), 21 deletions(-)

diff --git a/src/common/turb/mode_prandtl.F90 b/src/common/turb/mode_prandtl.F90
index 9d27bd513..2ba414d71 100644
--- a/src/common/turb/mode_prandtl.F90
+++ b/src/common/turb/mode_prandtl.F90
@@ -488,9 +488,9 @@ DO JSV=1,KSV
 !
     !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)    
     IF (KRR /= 0) THEN
-      ZWORK1 = ZW1(:,:,:)*PETHETA(:,:,:)
+      ZWORK1(:,:,:) = ZW1(:,:,:)*PETHETA(:,:,:)
     ELSE
-      ZWORK1 = ZW1(:,:,:)
+      ZWORK1(:,:,:) = ZW1(:,:,:)
     END IF
     PRED2THS3(:,:,:,JSV) = PREDTH1(:,:,:) * PREDS1(:,:,:,JSV)   +        &
                        ZWORK1(:,:,:) * ZWORK2(:,:,:)
@@ -519,9 +519,9 @@ DO JSV=1,KSV
 
     !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)    
     IF (KRR /= 0) THEN
-      ZWORK1 = ZW1(:,:,:)*PETHETA(:,:,:)
+      ZWORK1(:,:,:) = ZW1(:,:,:)*PETHETA(:,:,:)
     ELSE
-      ZWORK1 = ZW1(:,:,:)
+      ZWORK1(:,:,:) = ZW1(:,:,:)
     END IF
     PRED2THS3(:,:,:,JSV) = PREDTH1(:,:,:) * PREDS1(:,:,:,JSV)   +        &
                        ZWORK1(:,:,:)*ZWORK2(:,:,:)
diff --git a/src/common/turb/mode_tke_eps_sources.F90 b/src/common/turb/mode_tke_eps_sources.F90
index b399fb1f5..526b14d70 100644
--- a/src/common/turb/mode_tke_eps_sources.F90
+++ b/src/common/turb/mode_tke_eps_sources.F90
@@ -311,8 +311,8 @@ ENDIF
 ! TKE must be greater than its minimum value
 ! CL : Now done at the end of the time step in ADVECTION_METSV for MesoNH
 IF(HPROGRAM/='MESONH') THEN
- GTKENEG =  ZRES <= CSTURB%XTKEMIN
  !$mnh_expand_where(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
+ GTKENEG(:,:,:) =  ZRES(:,:,:) <= CSTURB%XTKEMIN
  WHERE ( GTKENEG(:,:,:) ) 
    ZRES(:,:,:) = CSTURB%XTKEMIN
  END WHERE
diff --git a/src/common/turb/mode_tridiag_thermo.F90 b/src/common/turb/mode_tridiag_thermo.F90
index a16d374ec..df45e1944 100644
--- a/src/common/turb/mode_tridiag_thermo.F90
+++ b/src/common/turb/mode_tridiag_thermo.F90
@@ -223,15 +223,19 @@ IF ( PIMPL > 1.E-10 ) THEN
   ZB(:,:,IKB) =   PRHODJ(:,:,IKB)/PTSTEP                   &
                 - ZRHODJ_DFDDTDZ_O_DZ2(:,:,IKB+D%NKL) * PIMPL
   ZC(:,:,IKB) =   ZRHODJ_DFDDTDZ_O_DZ2(:,:,IKB+D%NKL) * PIMPL
+  !$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)
     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
+    !$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)
   ZA(:,:,IKE) =   ZRHODJ_DFDDTDZ_O_DZ2(:,:,IKE  ) * PIMPL
   ZB(:,:,IKE) =   PRHODJ(:,:,IKE)/PTSTEP                   &
                 - ZRHODJ_DFDDTDZ_O_DZ2(:,:,IKE  ) * PIMPL
@@ -241,39 +245,45 @@ IF ( PIMPL > 1.E-10 ) THEN
 !
   ZBET(:,:) = ZB(:,:,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
+    !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
     ZGAM(:,:,JK) = ZC(:,:,JK-D%NKL) / 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-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 
   ! special treatment for the last level
+  !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
   ZGAM(:,:,IKE) = ZC(:,:,IKE-D%NKL) / ZBET(:,:) 
                                                     ! gam(k) = c(k-1) / bet
   ZBET(:,:)     = ZB(:,:,IKE) - ZA(:,:,IKE) * ZGAM(:,:,IKE)
                                                     ! bet = b(k) - a(k)* gam(k)  
   PVARP(:,:,IKE)= ( ZY(:,:,IKE) - ZA(:,:,IKE) * PVARP(:,:,IKE-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)
 !
 !*       3.3 going down
 !            ----------
 !
   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
 !
-  !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
 ELSE
 ! 
-  !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
   DO JK=IKTB,IKTE
+    !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
     PVARP(:,:,JK) = ZY(:,:,JK) * PTSTEP / PRHODJ(:,:,JK)
+    !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
   END DO
-  !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
 !
 END IF 
 !
diff --git a/src/common/turb/mode_tridiag_tke.F90 b/src/common/turb/mode_tridiag_tke.F90
index 81364c2f1..e69a64f73 100644
--- a/src/common/turb/mode_tridiag_tke.F90
+++ b/src/common/turb/mode_tridiag_tke.F90
@@ -232,11 +232,11 @@ IF ( PIMPL > 1.E-10 ) THEN
 !
 ELSE
 ! 
-  !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
   DO JK=IKTB,IKTE
+    !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
     PVARP(:,:,JK) = ZY(:,:,JK)
+    !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
   END DO
-  !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
 !
 END IF 
 !
diff --git a/src/common/turb/mode_tridiag_wind.F90 b/src/common/turb/mode_tridiag_wind.F90
index ca73dd9e1..27faf80c3 100644
--- a/src/common/turb/mode_tridiag_wind.F90
+++ b/src/common/turb/mode_tridiag_wind.F90
@@ -235,11 +235,11 @@ IF ( PIMPL > 1.E-10 ) THEN
 !
 ELSE
 ! 
-  !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
   DO JK=IKTB,IKTE
+    !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
     PVARP(:,:,JK) = ZY(:,:,JK)
+    !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
   END DO
-  !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
 !
 END IF 
 !
diff --git a/src/common/turb/mode_turb_ver.F90 b/src/common/turb/mode_turb_ver.F90
index 4cbdcbecb..5b8945071 100644
--- a/src/common/turb/mode_turb_ver.F90
+++ b/src/common/turb/mode_turb_ver.F90
@@ -436,7 +436,7 @@ END IF
 ! Square root of Tke
 !
 !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-ZSQRT_TKE = SQRT(PTKEM)
+ZSQRT_TKE(:,:,:) = SQRT(PTKEM(:,:,:))
 !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
 !
 ! gradients of mean quantities at previous time-step
@@ -448,13 +448,13 @@ IF (KRR>0) ZDR_DZ  = GZ_M_W(D%NKA, D%NKU, D%NKL,PRM(:,:,:,1),PDZZ)
 !
 ! Denominator factor in 3rd order terms
 !
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
 IF (.NOT. OHARAT) THEN
+  !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
   ZD(:,:,:) = (1.+ZREDTH1(:,:,:)+ZREDR1(:,:,:)) * (1.+0.5*(ZREDTH1(:,:,:)+ZREDR1(:,:,:)))
+  !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
 ELSE
   ZD(:,:,:) = 1.
 ENDIF
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
 !
 ! Phi3 and Psi3 Prandtl numbers
 !
diff --git a/src/common/turb/mode_turb_ver_sv_flux.F90 b/src/common/turb/mode_turb_ver_sv_flux.F90
index e4e1b20f9..e960864b4 100644
--- a/src/common/turb/mode_turb_ver_sv_flux.F90
+++ b/src/common/turb/mode_turb_ver_sv_flux.F90
@@ -418,9 +418,13 @@ DO JSV=1,KSV
     !
     !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
     ZFLXZ(:,:,D%NKA) = ZFLXZ(:,:,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)
       PWSV(:,:,JK,JSV)=0.5*(ZFLXZ(:,:,JK)+ZFLXZ(:,:,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)
     PWSV(:,:,IKB,JSV)=0.5*(ZFLXZ(:,:,IKB)+ZFLXZ(:,:,IKB+D%NKL))
     PWSV(:,:,IKE,JSV)=PWSV(:,:,IKE-D%NKL,JSV)
     !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
diff --git a/src/common/turb/mode_turb_ver_thermo_corr.F90 b/src/common/turb/mode_turb_ver_thermo_corr.F90
index 7e2badb3c..529bae238 100644
--- a/src/common/turb/mode_turb_ver_thermo_corr.F90
+++ b/src/common/turb/mode_turb_ver_thermo_corr.F90
@@ -472,7 +472,7 @@ ENDIF
       !
       !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
       ZF(:,:,:) = ZF(:,:,:) + ZWORK1(:,:,:) * PFR2(:,:,:)
-      ZDFDDTDZ = ZDFDDTDZ + ZWORK2(:,:,:) * PFR2(:,:,:)
+      ZDFDDTDZ(:,:,:) = ZDFDDTDZ(:,:,:) + ZWORK2(:,:,:) * PFR2(:,:,:)
       !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
     END IF
     !
diff --git a/src/common/turb/mode_turb_ver_thermo_flux.F90 b/src/common/turb/mode_turb_ver_thermo_flux.F90
index 6bf75778d..29d2228ec 100644
--- a/src/common/turb/mode_turb_ver_thermo_flux.F90
+++ b/src/common/turb/mode_turb_ver_thermo_flux.F90
@@ -708,11 +708,13 @@ IF (OOCEAN) THEN
   !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
 END IF
 !
-!$mnh_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)
   PWTH(:,:,JK)=0.5*(ZFLXZ(:,:,JK)+ZFLXZ(:,:,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)
 PWTH(:,:,IKB)=0.5*(ZFLXZ(:,:,IKB)+ZFLXZ(:,:,IKB+D%NKL)) 
 !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)    
 !
@@ -772,7 +774,7 @@ END IF
 !
 ZWORK1 = MZM(PETHETA, D%NKA, D%NKU, D%NKL)
 !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
-PWTHV = ZWORK1(:,:,:) * ZFLXZ
+PWTHV(:,:,:) = ZWORK1(:,:,:) * ZFLXZ(:,:,:)
 !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT,JK=1:D%NKT)
 !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
 PWTHV(:,:,IKB) = PETHETA(:,:,IKB) * ZFLXZ(:,:,IKB)
@@ -780,9 +782,9 @@ PWTHV(:,:,IKB) = PETHETA(:,:,IKB) * ZFLXZ(:,:,IKB)
 !
 IF (OOCEAN) THEN
   ! temperature contribution to Buy flux
-!$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT) 
+  !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT) 
   PWTHV(:,:,IKE) = PETHETA(:,:,IKE) * ZFLXZ(:,:,IKE)
-!$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+  !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
 END IF
 !*       2.3  Partial vertical divergence of the < Rc w > flux
 ! Correction for qc and qi negative in AROME 
@@ -1044,10 +1046,14 @@ IF (KRR /= 0) THEN
   !
   !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
   ZFLXZ(:,:,D%NKA) = ZFLXZ(:,:,IKB) 
+  !$mnh_end_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
   !
   DO JK=IKTB+1,IKTE-1
-   PWRC(:,:,JK)=0.5*(ZFLXZ(:,:,JK)+ZFLXZ(:,:,JK+D%NKL))
+    !$mnh_expand_array(JI=1:D%NIT,JJ=1:D%NJT)
+    PWRC(:,:,JK)=0.5*(ZFLXZ(:,:,JK)+ZFLXZ(:,:,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)
   PWRC(:,:,IKB)=0.5*(ZFLXZ(:,:,IKB)+ZFLXZ(:,:,IKB+D%NKL))
   PWRC(:,:,D%NKA)=0.5*(ZFLXZ(:,:,D%NKA)+ZFLXZ(:,:,D%NKA+D%NKL))
   PWRC(:,:,IKE)=PWRC(:,:,IKE-D%NKL)
-- 
GitLab