From 9aae7ac91368801b3ff37df07c82dcf5850aaf5c Mon Sep 17 00:00:00 2001
From: Gaelle DELAUTIER <gaelle.delautier@meteo.fr>
Date: Mon, 12 Dec 2016 15:26:08 +0100
Subject: [PATCH] Gaelle 12/12/16 : bug CPL_WAVE

---
 src/ARCH_SRC/CPL_WAVE/MNH/ground_paramn.f90   |   2 +-
 .../CPL_WAVE/MNH/mnh_oasis_define.F90         |   2 +-
 .../CPL_WAVE/SURFEX/diag_seaflux_initn.F90    |   2 +-
 src/ARCH_SRC/CPL_WAVE/SURFEX/read_lcover.F90  |   2 +-
 .../CPL_WAVE/SURFEX/read_seafluxn.F90         |   4 +-
 src/ARCH_SRC/CPL_WAVE/SURFEX/write_lcover.F90 |   2 +-
 .../CPL_WAVE/SURFEX/writesurf_seafluxn.F90    |   2 +-
 src/MNH/compute_entr_detr.f90                 | 215 ++++++++----------
 8 files changed, 109 insertions(+), 122 deletions(-)

diff --git a/src/ARCH_SRC/CPL_WAVE/MNH/ground_paramn.f90 b/src/ARCH_SRC/CPL_WAVE/MNH/ground_paramn.f90
index a5b570de7..a4ec50b76 100644
--- a/src/ARCH_SRC/CPL_WAVE/MNH/ground_paramn.f90
+++ b/src/ARCH_SRC/CPL_WAVE/MNH/ground_paramn.f90
@@ -119,8 +119,8 @@ USE MODD_IO_SURF_MNH, ONLY : NLUOUT
 USE MODI_GET_HALO
 USE MODI_MNH_OASIS_RECV
 USE MODI_MNH_OASIS_SEND
-USE MODD_DYN, ONLY : XSEGLEN
 USE MODD_SFX_OASIS, ONLY : LOASIS
+USE MODD_DYN, ONLY : XSEGLEN
 #endif
 !
 USE MODD_CST,        ONLY : XP00, XCPD, XRD, XRV,XRHOLW, XDAY, XPI, XLVTT, XMD, XAVOGADRO
diff --git a/src/ARCH_SRC/CPL_WAVE/MNH/mnh_oasis_define.F90 b/src/ARCH_SRC/CPL_WAVE/MNH/mnh_oasis_define.F90
index 06afba192..7fca9653c 100755
--- a/src/ARCH_SRC/CPL_WAVE/MNH/mnh_oasis_define.F90
+++ b/src/ARCH_SRC/CPL_WAVE/MNH/mnh_oasis_define.F90
@@ -64,7 +64,7 @@ USE MOD_OASIS
 USE MODD_MNH_SURFEX_n
 #endif
 !
-USE MODD_CONF, ONLY : NVERB
+USE MODD_CONF, ONLY : NVERB,NHALO
 USE MODD_PARAMETERS, ONLY : JPHEXT
 !
 USE MODE_ll
diff --git a/src/ARCH_SRC/CPL_WAVE/SURFEX/diag_seaflux_initn.F90 b/src/ARCH_SRC/CPL_WAVE/SURFEX/diag_seaflux_initn.F90
index f822e8e2a..7d9588ed4 100644
--- a/src/ARCH_SRC/CPL_WAVE/SURFEX/diag_seaflux_initn.F90
+++ b/src/ARCH_SRC/CPL_WAVE/SURFEX/diag_seaflux_initn.F90
@@ -93,7 +93,7 @@ INTEGER, INTENT(IN) :: KSW   ! number of SW spectral bands
 !
 INTEGER           :: IVERSION
 INTEGER           :: IRESP          ! IRESP  : return-code if a problem appears
- CHARACTER(LEN=12) :: YREC           ! Name of the article to be read
+ CHARACTER(LEN=LEN_HREC) :: YREC           ! Name of the article to be read
 !
 REAL(KIND=JPRB)   :: ZHOOK_HANDLE
 !
diff --git a/src/ARCH_SRC/CPL_WAVE/SURFEX/read_lcover.F90 b/src/ARCH_SRC/CPL_WAVE/SURFEX/read_lcover.F90
index 6855dbefa..2a3c869bd 100644
--- a/src/ARCH_SRC/CPL_WAVE/SURFEX/read_lcover.F90
+++ b/src/ARCH_SRC/CPL_WAVE/SURFEX/read_lcover.F90
@@ -73,7 +73,7 @@ LOGICAL, DIMENSION(JPCOVER)    :: OCOVER   ! list of covers
 !              -------------------------------
 !
 INTEGER           :: IRESP          ! Error code after redding
- CHARACTER(LEN=12) :: YRECFM         ! Name of the article to be read
+ CHARACTER(LEN=LEN_HREC) :: YRECFM         ! Name of the article to be read
 INTEGER           :: IVERSION       ! version of surfex file being read
 LOGICAL, DIMENSION(:), ALLOCATABLE :: GCOVER ! cover list in the file
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
diff --git a/src/ARCH_SRC/CPL_WAVE/SURFEX/read_seafluxn.F90 b/src/ARCH_SRC/CPL_WAVE/SURFEX/read_seafluxn.F90
index c65ac933c..d77db880c 100644
--- a/src/ARCH_SRC/CPL_WAVE/SURFEX/read_seafluxn.F90
+++ b/src/ARCH_SRC/CPL_WAVE/SURFEX/read_seafluxn.F90
@@ -87,7 +87,7 @@ INTEGER           :: ILU          ! 1D physical dimension
 !
 INTEGER           :: IRESP          ! Error code after redding
 !
- CHARACTER(LEN=12) :: YRECFM         ! Name of the article to be read
+ CHARACTER(LEN=LEN_HREC) :: YRECFM         ! Name of the article to be read
 !
 INTEGER           :: IVERSION       ! surface version
 !
@@ -256,7 +256,7 @@ SUBROUTINE CHECK_SEA(HFIELD,PFIELD)
 !
 IMPLICIT NONE
 !
- CHARACTER(LEN=12),  INTENT(IN) :: HFIELD
+ CHARACTER(LEN=LEN_HREC),  INTENT(IN) :: HFIELD
 REAL, DIMENSION(:), INTENT(IN) :: PFIELD
 !
 REAL            :: ZMAX,ZMIN
diff --git a/src/ARCH_SRC/CPL_WAVE/SURFEX/write_lcover.F90 b/src/ARCH_SRC/CPL_WAVE/SURFEX/write_lcover.F90
index f4a6ec576..33061c51b 100644
--- a/src/ARCH_SRC/CPL_WAVE/SURFEX/write_lcover.F90
+++ b/src/ARCH_SRC/CPL_WAVE/SURFEX/write_lcover.F90
@@ -73,7 +73,7 @@ LOGICAL, DIMENSION(JPCOVER)    :: OCOVER   ! list of covers
 TYPE(DIAG_SURF_ATM_t), INTENT(INOUT) :: DGU
 TYPE(SURF_ATM_t), INTENT(INOUT) :: U
 INTEGER           :: IRESP          ! Error code after reading
-CHARACTER(LEN=12) :: YRECFM         ! Name of the article to be read
+CHARACTER(LEN=LEN_HREC) :: YRECFM         ! Name of the article to be read
 CHARACTER(LEN=100):: YCOMMENT       ! Comment string
 LOGICAL, DIMENSION(JPCOVER)    :: GCOVER   ! tmp list of covers
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
diff --git a/src/ARCH_SRC/CPL_WAVE/SURFEX/writesurf_seafluxn.F90 b/src/ARCH_SRC/CPL_WAVE/SURFEX/writesurf_seafluxn.F90
index 35f7469b9..eb75fefa4 100644
--- a/src/ARCH_SRC/CPL_WAVE/SURFEX/writesurf_seafluxn.F90
+++ b/src/ARCH_SRC/CPL_WAVE/SURFEX/writesurf_seafluxn.F90
@@ -90,7 +90,7 @@ INTEGER           :: JMTH, INMTH
  CHARACTER(LEN=2 ) :: YMTH
 !
 INTEGER           :: IRESP          ! IRESP  : return-code if a problem appears
- CHARACTER(LEN=12) :: YRECFM         ! Name of the article to be read
+ CHARACTER(LEN=LEN_HREC) :: YRECFM         ! Name of the article to be read
  CHARACTER(LEN=100):: YCOMMENT       ! Comment string
 REAL(KIND=JPRB) :: ZHOOK_HANDLE
 !
diff --git a/src/MNH/compute_entr_detr.f90 b/src/MNH/compute_entr_detr.f90
index 413f40c44..b7346fad9 100644
--- a/src/MNH/compute_entr_detr.f90
+++ b/src/MNH/compute_entr_detr.f90
@@ -234,33 +234,35 @@ INTEGER :: JI,JLOOP
   ZMIXRT(:)=0.1
 
 !                1.4 Estimation of PPART_DRY
-  WHERE(OTEST .AND. OTESTLCL)
+  DO JLOOP=1,SIZE(OTEST)
+    IF(OTEST(JLOOP) .AND. OTESTLCL(JLOOP)) THEN
       !No dry part when condensation level is reached
-      PPART_DRY(:)=0.
-      ZDZ_STOP(:)=0.
-      ZPRE(:)=PPRE_MINUS_HALF(:)
-  ELSEWHERE (OTEST .AND. .NOT. OTESTLCL)
+      PPART_DRY(JLOOP)=0.
+      ZDZ_STOP(JLOOP)=0.
+      ZPRE(JLOOP)=PPRE_MINUS_HALF(JLOOP)
+    ELSE IF (OTEST(JLOOP) .AND. .NOT. OTESTLCL(JLOOP)) THEN
       !Temperature at flux level KK
-      ZT(:)=PTH_UP(:)*(PPRE_MINUS_HALF(:)/XP00) ** (XRD/XCPD)
+      ZT(JLOOP)=PTH_UP(JLOOP)*(PPRE_MINUS_HALF(JLOOP)/XP00) ** (XRD/XCPD)
       !Saturating vapor pressure at flux level KK
-      ZFOESW(:) = MIN(EXP( XALPW - XBETAW/ZT(:) - XGAMW*LOG(ZT(:))  ), 0.99*PPRE_MINUS_HALF(:))
-      ZFOESI(:) = MIN(EXP( XALPI - XBETAI/ZT(:) - XGAMI*LOG(ZT(:))  ), 0.99*PPRE_MINUS_HALF(:))
+      ZFOESW(JLOOP) = MIN(EXP( XALPW - XBETAW/ZT(JLOOP) - XGAMW*LOG(ZT(JLOOP))  ), 0.99*PPRE_MINUS_HALF(JLOOP))
+      ZFOESI(JLOOP) = MIN(EXP( XALPI - XBETAI/ZT(JLOOP) - XGAMI*LOG(ZT(JLOOP))  ), 0.99*PPRE_MINUS_HALF(JLOOP))
       !Computation of d.Rsat / dP (partial derivations with respect to P and T
       !and use of T=Theta*(P/P0)**(R/Cp) to transform dT into dP with theta_up
       !constant at the vertical)
-      ZDRSATODP(:)=(XBETAW/ZT(:)-XGAMW)*(1-ZFRAC_ICE(:))+(XBETAI/ZT(:)-XGAMI)*ZFRAC_ICE(:)
-      ZDRSATODP(:)=((XRD/XCPD)*ZDRSATODP(:)-1.)*PRSAT_UP(:)/ &
-                  &(PPRE_MINUS_HALF(:)-(ZFOESW(:)*(1-ZFRAC_ICE(:)) + ZFOESI(:)*ZFRAC_ICE(:)))
+      ZDRSATODP(JLOOP)=(XBETAW/ZT(JLOOP)-XGAMW)*(1-ZFRAC_ICE(JLOOP))+(XBETAI/ZT(JLOOP)-XGAMI)*ZFRAC_ICE(JLOOP)
+      ZDRSATODP(JLOOP)=((XRD/XCPD)*ZDRSATODP(JLOOP)-1.)*PRSAT_UP(JLOOP)/ &
+                  &(PPRE_MINUS_HALF(JLOOP)-(ZFOESW(JLOOP)*(1-ZFRAC_ICE(JLOOP)) + ZFOESI(JLOOP)*ZFRAC_ICE(JLOOP)))
       !Use of d.Rsat / dP and pressure at flux level KK to find pressure (ZPRE)
       !where Rsat is equal to PRT_UP
-      ZPRE(:)=PPRE_MINUS_HALF(:)+(PRT_UP(:)-PRSAT_UP(:))/ZDRSATODP(:)
+      ZPRE(JLOOP)=PPRE_MINUS_HALF(JLOOP)+(PRT_UP(JLOOP)-PRSAT_UP(JLOOP))/ZDRSATODP(JLOOP)
       !Fraction of dry part (computed with pressure and used with heights, no
       !impact found when using log function here and for pressure on flux levels
       !computation)
-      PPART_DRY(:)=MAX(0., MIN(1., (PPRE_MINUS_HALF(:)-ZPRE(:))/(PPRE_MINUS_HALF(:)-PPRE_PLUS_HALF(:))))
+      PPART_DRY(JLOOP)=MAX(0., MIN(1., (PPRE_MINUS_HALF(JLOOP)-ZPRE(JLOOP))/(PPRE_MINUS_HALF(JLOOP)-PPRE_PLUS_HALF(JLOOP))))
       !Height above flux level KK of the cloudy part
-      ZDZ_STOP(:) = (PZZ(:,KK+KKL)-PZZ(:,KK))*PPART_DRY(:)
-  ENDWHERE
+      ZDZ_STOP(JLOOP) = (PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK))*PPART_DRY(JLOOP)
+    END IF
+  END DO
 
 !               1.5 Gradient and flux values of thetav
   IF(KK/=KKB)THEN
@@ -332,29 +334,31 @@ ENDDO
   ZTHV_UP_F2(:) = ZTHMIX(:)*(1.+ZRVORD*ZRVMIX(:))/(1.+PRT_UP(:))
 
   ! Integral buoyancy for cloudy part
-  WHERE(OTEST .AND. PPART_DRY(:)<1.)
-    !Gradient of Theta V updraft over the cloudy part, assuming that thetaV updraft don't change
-    !between flux level KK and bottom of cloudy part
-    ZCOTHVU(:)=(ZTHV_UP_F2(:)-PTHV_UP(:))/((PZZ(:,KK+KKL)-PZZ(:,KK))*(1-PPART_DRY(:)))
-
-    !Computation in two parts to use change of gradient of theta v of environment
-    !Between bottom of cloudy part (if under mass level) and mass level KK
-    ZDZ(:)=MAX(0., 0.5*(PZZ(:,KK+KKL)-PZZ(:,KK))-ZDZ_STOP(:))
-    PBUO_INTEG_CLD(:) = ZG_O_THVREF(:)*ZDZ(:)*&
-              (0.5*( ZCOTHVU(:) - ZCOEFF_MINUS_HALF(:))*ZDZ(:) &
-                - (PTHVM(:,KK)-ZDZ(:)*ZCOEFF_MINUS_HALF(:)) + PTHV_UP(:) )
-
-    !Between max(mass level, bottom of cloudy part) and flux level KK+KKL
-    ZDZ(:)=(PZZ(:,KK+KKL)-PZZ(:,KK))-MAX(ZDZ_STOP(:),0.5*(PZZ(:,KK+KKL)-PZZ(:,KK)))
-    PBUO_INTEG_CLD(:) = PBUO_INTEG_CLD(:)+ZG_O_THVREF(:)*ZDZ(:)*&
-                        (0.5*( ZCOTHVU(:) - ZCOEFF_PLUS_HALF(:))*ZDZ(:)&
-                - (PTHVM(:,KK)+(0.5*((PZZ(:,KK+KKL)-PZZ(:,KK)))-ZDZ(:))*ZCOEFF_PLUS_HALF(:)) +&
-                PTHV_UP(:) )
-
-  ELSEWHERE
-    !No cloudy part
-    PBUO_INTEG_CLD(:)=0.
-  ENDWHERE
+  DO JLOOP=1,SIZE(OTEST)
+    IF(OTEST(JLOOP) .AND. PPART_DRY(JLOOP)<1.) THEN
+      !Gradient of Theta V updraft over the cloudy part, assuming that thetaV updraft don't change
+      !between flux level KK and bottom of cloudy part
+      ZCOTHVU(JLOOP)=(ZTHV_UP_F2(JLOOP)-PTHV_UP(JLOOP))/((PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK))*(1-PPART_DRY(JLOOP)))
+
+      !Computation in two parts to use change of gradient of theta v of environment
+      !Between bottom of cloudy part (if under mass level) and mass level KK
+      ZDZ(JLOOP)=MAX(0., 0.5*(PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK))-ZDZ_STOP(JLOOP))
+      PBUO_INTEG_CLD(JLOOP) = ZG_O_THVREF(JLOOP)*ZDZ(JLOOP)*&
+              (0.5*( ZCOTHVU(JLOOP) - ZCOEFF_MINUS_HALF(JLOOP))*ZDZ(JLOOP) &
+                - (PTHVM(JLOOP,KK)-ZDZ(JLOOP)*ZCOEFF_MINUS_HALF(JLOOP)) + PTHV_UP(JLOOP) )
+
+      !Between max(mass level, bottom of cloudy part) and flux level KK+KKL
+      ZDZ(JLOOP)=(PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK))-MAX(ZDZ_STOP(JLOOP),0.5*(PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK)))
+      PBUO_INTEG_CLD(JLOOP) = PBUO_INTEG_CLD(JLOOP)+ZG_O_THVREF(JLOOP)*ZDZ(JLOOP)*&
+                        (0.5*( ZCOTHVU(JLOOP) - ZCOEFF_PLUS_HALF(JLOOP))*ZDZ(JLOOP)&
+                - (PTHVM(JLOOP,KK)+(0.5*((PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK)))-ZDZ(JLOOP))*ZCOEFF_PLUS_HALF(JLOOP)) +&
+                PTHV_UP(JLOOP) )
+
+    ELSE
+      !No cloudy part
+      PBUO_INTEG_CLD(JLOOP)=0.
+    END IF
+  END DO
 
 !               3.2 Critical mixed fraction for KK+KKL flux level (ZKIC_F2) and
 !                   for bottom of cloudy part (ZKIC), then a mean for the cloudy part
@@ -372,26 +376,27 @@ ENDDO
 
   !   JI computed to avoid KKL(KK-KKL) being < KKL*KKB
   JI=KKL*MAX(KKL*(KK-KKL),KKL*KKB)
-
-  WHERE(OTEST .AND. PPART_DRY(:)>0.5)
-    ZDZ(:)=ZDZ_STOP(:)-0.5*(PZZ(:,KK+KKL)-PZZ(:,KK))
-    ZTHV(:)= PTHVM(:,KK)+ZCOEFF_PLUS_HALF(:)*ZDZ(:)
-    ZMIXTHL(:) = ZKIC_INIT * &
-                 (PTHLM(:,KK)+ZDZ(:)*(PTHLM(:,KK+KKL)-PTHLM(:,KK))/PDZZ(:,KK+KKL)) + &
-                 (1. - ZKIC_INIT)*PTHL_UP(:)
-    ZMIXRT(:)  = ZKIC_INIT * &
-                 (PRTM(:,KK)+ZDZ(:)*(PRTM(:,KK+KKL)-PRTM(:,KK))/PDZZ(:,KK+KKL)) +   &
-                 (1. - ZKIC_INIT)*PRT_UP(:)
-  ELSEWHERE(OTEST)
-    ZDZ(:)=0.5*(PZZ(:,KK+KKL)-PZZ(:,KK))-ZDZ_STOP(:)
-    ZTHV(:)= PTHVM(:,KK)-ZCOEFF_MINUS_HALF(:)*ZDZ(:)
-    ZMIXTHL(:) = ZKIC_INIT * &
-                 (PTHLM(:,KK)-ZDZ(:)*(PTHLM(:,KK)-PTHLM(:,JI))/PDZZ(:,KK)) + &
-                 (1. - ZKIC_INIT)*PTHL_UP(:)
-    ZMIXRT(:)  = ZKIC_INIT * &
-                 (PRTM(:,KK)-ZDZ(:)*(PRTM(:,KK)-PRTM(:,JI))/PDZZ(:,KK)) + &
-                 (1. - ZKIC_INIT)*PRT_UP(:)
-  ENDWHERE
+  DO JLOOP=1,SIZE(OTEST)
+    IF(OTEST(JLOOP) .AND. PPART_DRY(JLOOP)>0.5) THEN
+      ZDZ(JLOOP)=ZDZ_STOP(JLOOP)-0.5*(PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK))
+      ZTHV(JLOOP)= PTHVM(JLOOP,KK)+ZCOEFF_PLUS_HALF(JLOOP)*ZDZ(JLOOP)
+      ZMIXTHL(JLOOP) = ZKIC_INIT * &
+                 (PTHLM(JLOOP,KK)+ZDZ(JLOOP)*(PTHLM(JLOOP,KK+KKL)-PTHLM(JLOOP,KK))/PDZZ(JLOOP,KK+KKL)) + &
+                 (1. - ZKIC_INIT)*PTHL_UP(JLOOP)
+      ZMIXRT(JLOOP)  = ZKIC_INIT * &
+                 (PRTM(JLOOP,KK)+ZDZ(JLOOP)*(PRTM(JLOOP,KK+KKL)-PRTM(JLOOP,KK))/PDZZ(JLOOP,KK+KKL)) +   &
+                 (1. - ZKIC_INIT)*PRT_UP(JLOOP)
+     ELSEIF(OTEST(JLOOP)) THEN
+      ZDZ(JLOOP)=0.5*(PZZ(JLOOP,KK+KKL)-PZZ(JLOOP,KK))-ZDZ_STOP(JLOOP)
+      ZTHV(JLOOP)= PTHVM(JLOOP,KK)-ZCOEFF_MINUS_HALF(JLOOP)*ZDZ(JLOOP)
+      ZMIXTHL(JLOOP) = ZKIC_INIT * &
+                 (PTHLM(JLOOP,KK)-ZDZ(JLOOP)*(PTHLM(JLOOP,KK)-PTHLM(JLOOP,JI))/PDZZ(JLOOP,KK)) + &
+                 (1. - ZKIC_INIT)*PTHL_UP(JLOOP)
+      ZMIXRT(JLOOP)  = ZKIC_INIT * &
+                 (PRTM(JLOOP,KK)-ZDZ(JLOOP)*(PRTM(JLOOP,KK)-PRTM(JLOOP,JI))/PDZZ(JLOOP,KK)) + &
+                 (1. - ZKIC_INIT)*PRT_UP(JLOOP)
+    ENDIF
+  ENDDO
   CALL TH_R_FROM_THL_RT_1D(HFRAC_ICE,ZFRAC_ICE,&
                ZPRE,ZMIXTHL,ZMIXRT,&
                ZTHMIX,ZRVMIX,PRC_MIX,PRI_MIX,&
@@ -408,50 +413,28 @@ ENDDO
   ZTHVMIX_F2(:) = ZTHMIX(:)*(1.+ZRVORD*ZRVMIX(:))/(1.+ZMIXRT(:))
 
   !Computation of mean ZKIC over the cloudy part
-  WHERE (OTEST)
+  DO JLOOP=1,SIZE(OTEST) 
+    IF (OTEST(JLOOP)) THEN
     ! Compute ZKIC at the bottom of cloudy part
     ! Thetav_up at bottom is equal to Thetav_up at flux level KK
-    WHERE (ABS(PTHV_UP(:)-ZTHVMIX(:))<1.E-10)
-      ZKIC(:)=1.
-    ELSEWHERE
-      ZKIC(:) = MAX(0.,PTHV_UP(:)-ZTHV(:))*ZKIC_INIT /  &  
-                   (PTHV_UP(:)-ZTHVMIX(:))
-    ENDWHERE
-    ! Compute ZKIC_F2 at flux level KK+KKL
-    WHERE (ABS(ZTHV_UP_F2(:)-ZTHVMIX_F2(:))<1.E-10)
-      ZKIC_F2(:)=1.
-    ELSEWHERE
-      ZKIC_F2(:) = MAX(0.,ZTHV_UP_F2(:)-ZTHV_PLUS_HALF(:))*ZKIC_INIT /  &  
-                   (ZTHV_UP_F2(:)-ZTHVMIX_F2(:))
-    ENDWHERE
-    !Mean ZKIC over the cloudy part
-    ZKIC(:)=MAX(MIN(0.5*(ZKIC(:)+ZKIC_F2(:)),1.),0.)
-  ENDWHERE
-
-
-
-!====
-  !Computation of mean ZKIC over the cloudy part
-!  WHERE (OTEST .AND. ABS(PTHV_UP(:)-ZTHVMIX(:))<1.E-10)
-!    ! Compute ZKIC at the bottom of cloudy part
-!    ! Thetav_up at bottom is equal to Thetav_up at flux level KK
-!      ZKIC(:)=1.
-!    ELSEWHERE (OTEST .AND. ABS(PTHV_UP(:)-ZTHVMIX(:))>=1.E-10)
-!      ZKIC(:) = MAX(0.,PTHV_UP(:)-ZTHV(:))*ZKIC_INIT /  &  
-!                   (PTHV_UP(:)-ZTHVMIX(:))
-!  ENDWHERE
-    ! Compute ZKIC_F2 at flux level KK+KKL
-!  WHERE (OTEST .AND. ABS(ZTHV_UP_F2(:)-ZTHVMIX_F2(:))<1.E-10)
-!      ZKIC_F2(:)=1.
-!  ELSEWHERE (OTEST .AND. ABS(ZTHV_UP_F2(:)-ZTHVMIX_F2(:))>=1.E-10)
-!      ZKIC_F2(:) = MAX(0.,ZTHV_UP_F2(:)-ZTHV_PLUS_HALF(:))*ZKIC_INIT /  &  
-!                   (ZTHV_UP_F2(:)-ZTHVMIX_F2(:))
-!  ENDWHERE
-!  WHERE (OTEST)
-!    !Mean ZKIC over the cloudy part
-!    ZKIC(:)=MAX(MIN(0.5*(ZKIC(:)+ZKIC_F2(:)),1.),0.)
-!  ENDWHERE
-!======
+      IF (ABS(PTHV_UP(JLOOP)-ZTHVMIX(JLOOP))<1.E-10) THEN
+        ZKIC(JLOOP)=1.
+      ELSE
+        ZKIC(JLOOP) = MAX(0.,PTHV_UP(JLOOP)-ZTHV(JLOOP))*ZKIC_INIT /  &  
+                   (PTHV_UP(JLOOP)-ZTHVMIX(JLOOP))
+      END IF
+      ! Compute ZKIC_F2 at flux level KK+KKL
+      IF (ABS(ZTHV_UP_F2(JLOOP)-ZTHVMIX_F2(JLOOP))<1.E-10) THEN
+        ZKIC_F2(JLOOP)=1.
+      ELSE
+        ZKIC_F2(JLOOP) = MAX(0.,ZTHV_UP_F2(JLOOP)-ZTHV_PLUS_HALF(JLOOP))*ZKIC_INIT /  &  
+                   (ZTHV_UP_F2(JLOOP)-ZTHVMIX_F2(JLOOP))
+      END IF
+      !Mean ZKIC over the cloudy part
+      ZKIC(JLOOP)=MAX(MIN(0.5*(ZKIC(JLOOP)+ZKIC_F2(JLOOP)),1.),0.)
+    END IF
+  END DO
+
 
 !               3.3 Integration of PDF
 !                   According to Kain and Fritsch (1990), we replace delta Mt
@@ -460,10 +443,12 @@ ENDDO
 
   !Constant PDF
   !For this PDF, eq. (5) is delta Me=0.5*delta Mt
-  WHERE(OTEST)
-    ZEPSI(:) = ZKIC(:)**2. !integration multiplied by 2
-    ZDELTA(:) = (1.-ZKIC(:))**2. !idem
-  ENDWHERE
+  DO JLOOP=1,SIZE(OTEST)
+    IF(OTEST(JLOOP)) THEN
+      ZEPSI(JLOOP) = ZKIC(JLOOP)**2. !integration multiplied by 2
+      ZDELTA(JLOOP) = (1.-ZKIC(JLOOP))**2. !idem
+    ENDIF
+  ENDDO
 
   !Triangular PDF
   !Calculus must be verified before activating this part, but in this state,
@@ -481,15 +466,17 @@ ENDDO
   !ENDWHERE
 
 !               3.4 Computation of PENTR and PDETR
-  WHERE (OTEST)
-    ZEPSI_CLOUD=MIN(ZDELTA,ZEPSI)
-    PENTR_CLD(:) = (1.-PPART_DRY(:))*ZCOEFFMF_CLOUD*PRHODREF(:)*ZEPSI_CLOUD(:)
-    PDETR_CLD(:) = (1.-PPART_DRY(:))*ZCOEFFMF_CLOUD*PRHODREF(:)*ZDELTA(:)
-    PENTR(:) = PENTR(:)+PENTR_CLD(:)
-    PDETR(:) = PDETR(:)+PDETR_CLD(:)
-  ELSEWHERE
-    PENTR_CLD(:) = 0.
-    PDETR_CLD(:) = 0.
-  ENDWHERE
+  DO JLOOP=1,SIZE(OTEST)
+    IF(OTEST(JLOOP)) THEN
+      ZEPSI_CLOUD=MIN(ZDELTA,ZEPSI)
+      PENTR_CLD(JLOOP) = (1.-PPART_DRY(JLOOP))*ZCOEFFMF_CLOUD*PRHODREF(JLOOP)*ZEPSI_CLOUD(JLOOP)
+      PDETR_CLD(JLOOP) = (1.-PPART_DRY(JLOOP))*ZCOEFFMF_CLOUD*PRHODREF(JLOOP)*ZDELTA(JLOOP)
+      PENTR(JLOOP) = PENTR(JLOOP)+PENTR_CLD(JLOOP)
+      PDETR(JLOOP) = PDETR(JLOOP)+PDETR_CLD(JLOOP)
+    ELSE
+      PENTR_CLD(JLOOP) = 0.
+      PDETR_CLD(JLOOP) = 0.
+    ENDIF
+  ENDDO 
 
 END SUBROUTINE COMPUTE_ENTR_DETR  
-- 
GitLab