From ce7f82fb4972a2533e4f7285766b823f4b32829f Mon Sep 17 00:00:00 2001
From: Philippe WAUTELET <philippe.wautelet@aero.obs-mip.fr>
Date: Tue, 2 Jul 2019 10:06:50 +0200
Subject: [PATCH] Philippe 01/07/2019: OpenACC: optimisation of ppm_s0_x/y/z_d
 for GPU

---
 src/MNH/ppm.f90 | 440 ++++++++++++++++++++++++++++--------------------
 1 file changed, 261 insertions(+), 179 deletions(-)

diff --git a/src/MNH/ppm.f90 b/src/MNH/ppm.f90
index c4dbf881c..6323a8f9e 100644
--- a/src/MNH/ppm.f90
+++ b/src/MNH/ppm.f90
@@ -6,6 +6,7 @@
 ! Modifications:
 !  P. Wautelet 05/2016-04/2018: new data structures and calls for I/O
 !  P. Wautelet 20/06/2019: OpenACC: correct intent of some dummy variables
+!  P. Wautelet 01/07/2019: OpenACC: optimisation of ppm_s0_x/y/z_d for GPU
 !-----------------------------------------------------------------
 #ifdef _OPENACC
 !
@@ -427,6 +428,7 @@ CONTAINS
 !-------------------------------------------------------------------------------
 !
 USE MODE_ll
+use mode_msg
 #ifndef _OPENACC
 USE MODI_SHUMAN
 #else
@@ -489,6 +491,7 @@ INTEGER                          :: IJS,IJN
 !END JUAN PPM_LL
 #else
 INTEGER                          :: I,J,K 
+integer                          :: ji, jj, jk
 !
 !!$!
 !!$! terms used in parabolic interpolation, dmq, qL, qR, dq, q6
@@ -524,33 +527,40 @@ GEAST = LEAST_ll()
 #ifndef _OPENACC
 CALL GET_HALO(PSRC)
 !
-PR=PSRC
-ZQL=PSRC
-ZQR=PSRC
-ZDQ=PSRC
-ZQ6=PSRC
-ZDMQ=PSRC
-ZQL0=PSRC
-ZQR0=PSRC
-ZQ60=PSRC
-ZFPOS=PSRC
-ZFNEG=PSRC
+PR   (:,:,:) = PSRC(:,:,:)
+ZQL  (:,:,:) = PSRC(:,:,:)
+ZQR  (:,:,:) = PSRC(:,:,:)
+ZDQ  (:,:,:) = PSRC(:,:,:)
+ZQ6  (:,:,:) = PSRC(:,:,:)
+ZDMQ (:,:,:) = PSRC(:,:,:)
+ZQL0 (:,:,:) = PSRC(:,:,:)
+ZQR0 (:,:,:) = PSRC(:,:,:)
+ZQ60 (:,:,:) = PSRC(:,:,:)
+ZFPOS(:,:,:) = PSRC(:,:,:)
+ZFNEG(:,:,:) = PSRC(:,:,:)
 #else
 CALL GET_HALO_D(PSRC,HDIR="01_X")
 !
 !$acc kernels 
-PR=PSRC
-ZQL=PSRC
-ZQR=PSRC
-ZDQ=PSRC
-ZQ6=PSRC
-ZDMQ=PSRC
-ZQL0=PSRC
-ZQR0=PSRC
-ZQ60=PSRC
+!$acc loop independent collapse(3)
+  do jk = 1, iku
+    do jj = 1, iju
+      do ji = 1, iiu
+        PR   (ji, jj, jk ) = PSRC(ji, jj, jk )
+        ZQL  (ji, jj, jk ) = PSRC(ji, jj, jk )
+        ZQR  (ji, jj, jk ) = PSRC(ji, jj, jk )
+        ZDQ  (ji, jj, jk ) = PSRC(ji, jj, jk )
+        ZQ6  (ji, jj, jk ) = PSRC(ji, jj, jk )
+        ZDMQ (ji, jj, jk ) = PSRC(ji, jj, jk )
+        ZQL0 (ji, jj, jk ) = PSRC(ji, jj, jk )
+        ZQR0 (ji, jj, jk ) = PSRC(ji, jj, jk )
+        ZQ60 (ji, jj, jk ) = PSRC(ji, jj, jk )
+    end do
+  end do
+end do
+!
 ZFPOS(:,1:IJS,:)=PSRC(:,1:IJS,:)
 ZFNEG(:,1:IJS,:)=PSRC(:,1:IJS,:)
-!
 ZFPOS(:,IJN:,:)=PSRC(:,IJN:,:)
 ZFNEG(:,IJN:,:)=PSRC(:,IJN:,:)
 !$acc end kernels
@@ -565,7 +575,7 @@ SELECT CASE ( HLBCX(1) ) ! X direction LBC type: (1) for left side
 !
 CASE ('CYCL','WALL')          ! In that case one must have HLBCX(1) == HLBCX(2)
 #ifdef _OPENACC
-PRINT *,'OPENACC: ppm::PPM_01_X CYCL/WALL boundaries not yet implemented'
+  call Print_msg( NVERB_ERROR, 'GEN', 'PPM_01_X', 'OpenACC: CYCL/WALL boundaries not yet implemented' )
 #endif
 !
 ! calculate dmq
@@ -671,7 +681,7 @@ PRINT *,'OPENACC: ppm::PPM_01_X CYCL/WALL boundaries not yet implemented'
 !
 !  ZFPOS(i) = Fct[ ZQR(i-1),PCR(i),ZDQ(i-1),ZQ6(i-1) ]
 !
-   ZFPOS(IIB:IIE+1,IJS:IJN,:) = ZQR(IIB-1:IIE,IJS:IJN,:) - 0.5*PCR(IIB:IIE+1,IJS:IJN,:) * &            
+   ZFPOS(IIB:IIE+1,IJS:IJN,:) = ZQR(IIB-1:IIE,IJS:IJN,:) - 0.5*PCR(IIB:IIE+1,IJS:IJN,:) * &
         (ZDQ(IIB-1:IIE,IJS:IJN,:) - (1.0 - 2.0*PCR(IIB:IIE+1,IJS:IJN,:)/3.0)        &
         * ZQ6(IIB-1:IIE,IJS:IJN,:))
 !
@@ -694,8 +704,7 @@ PRINT *,'OPENACC: ppm::PPM_01_X CYCL/WALL boundaries not yet implemented'
    PR = DXF( PCR*MXM(PRHO)*( ZFPOS*(0.5+SIGN(0.5,PCR)) + & 
                              ZFNEG*(0.5-SIGN(0.5,PCR)) ) )
 #else
-PRINT *,'not yet implemented'
-STOP
+  call Print_msg( NVERB_ERROR, 'GEN', 'PPM_01_X', 'OpenACC: CYCL/WALL boundaries not yet implemented' )
 #endif
    CALL GET_HALO(PR)   
 !
@@ -733,11 +742,20 @@ CASE('OPEN')
 !
 !  ZDMQ(i) = Fct[ ZDMQ(i),PSRC(i),PSRC(i-1),PSRC(i+1) ]
 !
-   ZDMQ(IIB:IIE,IJS:IJN,:) = &
-        SIGN( (MIN( ABS(ZDMQ(IIB:IIE,IJS:IJN,:)),2.0*(PSRC(IIB:IIE,IJS:IJN,:) - &
-        MIN(PSRC(IIB-1:IIE-1,IJS:IJN,:),PSRC(IIB:IIE,IJS:IJN,:),PSRC(IIB+1:IIE+1,IJS:IJN,:))),    &
-        2.0*(MAX(PSRC(IIB-1:IIE-1,IJS:IJN,:),PSRC(IIB:IIE,IJS:IJN,:),PSRC(IIB+1:IIE+1,IJS:IJN,:)) - &
-        PSRC(IIB:IIE,IJS:IJN,:)) )), ZDMQ(IIB:IIE,IJS:IJN,:) )
+!$acc loop independent collapse(3)
+  do jk = 1, iku
+    do jj = ijs, ijn
+      do ji = iib, iie
+   ZDMQ(ji, jj, jk ) = SIGN( &
+            MIN( ABS(ZDMQ(ji, jj, jk )),                                                                       &
+                 2.0 * ( PSRC(ji, jj, jk )                                                                     &
+                         - MIN(PSRC(ji - 1, jj, jk ),PSRC(ji, jj, jk ),PSRC(ji + 1, jj, jk )) ),   &
+                 2.0 * (-PSRC(ji, jj, jk )                                                                     &
+                        + MAX(PSRC(ji - 1, jj, jk ),PSRC(ji, jj, jk ),PSRC(ji + 1, jj, jk )) )  ), &
+            ZDMQ(ji, jj, jk ) )
+      end do
+    end do
+  end do
 !
 !  WEST BOUND
 !
@@ -770,8 +788,14 @@ CASE('OPEN')
 !
 !  ZQL0(i) = Fct[ PSRC(i),PSRC(i-1),ZDMQ(i),ZDMQ(i-1) ]
 !
-   ZQL0(IIB:IIE+1,IJS:IJN,:) = 0.5*(PSRC(IIB:IIE+1,IJS:IJN,:) + PSRC(IIB-1:IIE,IJS:IJN,:)) - &
-        (ZDMQ(IIB:IIE+1,IJS:IJN,:) - ZDMQ(IIB-1:IIE,IJS:IJN,:))/6.0
+!$acc loop independent collapse(3)
+  do jk = 1, iku
+    do jj = ijs, ijn
+      do ji = iib, iie + 1
+        ZQL0(ji, jj, jk ) = 0.5 * ( PSRC(ji, jj, jk ) + PSRC(ji-1, jj, jk ) ) - ( ZDMQ(ji, jj, jk ) - ZDMQ(ji-1, jj, jk ) ) / 6.0
+      end do
+    end do
+  end do
 !
 #ifndef _OPENACC
    CALL  GET_HALO(ZQL0)
@@ -821,19 +845,11 @@ CASE('OPEN')
       ZQL(:,IJS:IJN,:) = PSRC(:,IJS:IJN,:)
       ZQR(:,IJS:IJN,:) = PSRC(:,IJS:IJN,:)
       ZQ6(:,IJS:IJN,:) = 0.0
-#ifndef MNH_BITREP
-   ELSEWHERE ( ZQ60(:,IJS:IJN,:)*ZDQ(:,IJS:IJN,:) < -(ZDQ(:,IJS:IJN,:))**2 )
-#else
-   ELSEWHERE ( ZQ60(:,IJS:IJN,:)*ZDQ(:,IJS:IJN,:) < -BR_P2(ZDQ(:,IJS:IJN,:)) )
-#endif
+   ELSEWHERE ( ZQ60(:,IJS:IJN,:)*ZDQ(:,IJS:IJN,:) < -ZDQ(:,IJS:IJN,:)*ZDQ(:,IJS:IJN,:) )
       ZQ6(:,IJS:IJN,:) = 3.0*(ZQL0(:,IJS:IJN,:) - PSRC(:,IJS:IJN,:))
       ZQR(:,IJS:IJN,:) = ZQL0(:,IJS:IJN,:) - ZQ6(:,IJS:IJN,:)
       ZQL(:,IJS:IJN,:) = ZQL0(:,IJS:IJN,:)
-#ifndef MNH_BITREP
-   ELSEWHERE ( ZQ60(:,IJS:IJN,:)*ZDQ(:,IJS:IJN,:) > (ZDQ(:,IJS:IJN,:))**2 )
-#else
-   ELSEWHERE ( ZQ60(:,IJS:IJN,:)*ZDQ(:,IJS:IJN,:) > BR_P2(ZDQ(:,IJS:IJN,:)) )
-#endif
+   ELSEWHERE ( ZQ60(:,IJS:IJN,:)*ZDQ(:,IJS:IJN,:) > ZDQ(:,IJS:IJN,:)*ZDQ(:,IJS:IJN,:) )
       ZQ6(:,IJS:IJN,:) = 3.0*(ZQR0(:,IJS:IJN,:) - PSRC(:,IJS:IJN,:))
       ZQL(:,IJS:IJN,:) = ZQR0(:,IJS:IJN,:) - ZQ6(:,IJS:IJN,:)
       ZQR(:,IJS:IJN,:) = ZQR0(:,IJS:IJN,:)
@@ -843,6 +859,7 @@ CASE('OPEN')
 !
    ZDQ(:,IJS:IJN,:) = ZQR(:,IJS:IJN,:) - ZQL(:,IJS:IJN,:)
 #else
+!$acc loop independent collapse(3)
 DO K=1,IKU 
  DO J = IJS,IJN
     ! acc loop vector(24)
@@ -857,7 +874,7 @@ DO K=1,IKU
 !
    ZQL(I,J,K) = ZQL0(I,J,K)
    ZQR(I,J,K) = ZQR0(I,J,K)
-   ZQ6(I,J,K) = ZQ60(I,J,K) 
+   ZQ6(I,J,K) = ZQ60(I,J,K)
 !
 ! eliminate over and undershoots and create qL and qR as in Lin96
 !
@@ -865,19 +882,11 @@ DO K=1,IKU
       ZQL(I,J,K) = PSRC(I,J,K)
       ZQR(I,J,K) = PSRC(I,J,K)
       ZQ6(I,J,K) = 0.0
-#ifndef MNH_BITREP
-   ELSEIF ( ZQ60(I,J,K)*ZDQ(I,J,K) < -(ZDQ(I,J,K))**2 ) THEN
-#else
-   ELSEIF ( ZQ60(I,J,K)*ZDQ(I,J,K) < -BR_P2(ZDQ(I,J,K)) ) THEN
-#endif
+   ELSEIF ( ZQ60(I,J,K)*ZDQ(I,J,K) < -ZDQ(I,J,K)*ZDQ(I,J,K) ) THEN
       ZQ6(I,J,K) = 3.0*(ZQL0(I,J,K) - PSRC(I,J,K))
       ZQR(I,J,K) = ZQL0(I,J,K) - ZQ6(I,J,K)
       ZQL(I,J,K) = ZQL0(I,J,K)
-#ifndef MNH_BITREP
-   ELSEIF ( ZQ60(I,J,K)*ZDQ(I,J,K) > (ZDQ(I,J,K))**2 ) THEN
-#else
-   ELSEIF ( ZQ60(I,J,K)*ZDQ(I,J,K) > BR_P2(ZDQ(I,J,K)) ) THEN
-#endif
+   ELSEIF ( ZQ60(I,J,K)*ZDQ(I,J,K) > ZDQ(I,J,K)*ZDQ(I,J,K) ) THEN
       ZQ6(I,J,K) = 3.0*(ZQR0(I,J,K) - PSRC(I,J,K))
       ZQL(I,J,K) = ZQR0(I,J,K) - ZQ6(I,J,K)
       ZQR(I,J,K) = ZQR0(I,J,K)
@@ -897,10 +906,17 @@ ENDDO ; ENDDO ; ENDDO
 !!$   ZFPOS(IIB+1:IIE+1,:,:) = ZQR(IIB:IIE,:,:) - 0.5*PCR(IIB+1:IIE+1,:,:) * &            
 !!$        (ZDQ(IIB:IIE,:,:) - (1.0 - 2.0*PCR(IIB+1:IIE+1,:,:)/3.0)          &
 !!$        * ZQ6(IIB:IIE,:,:))
-
-  ZFPOS(IIB:IIE+1,IJS:IJN,:) = ZQR(IIB-1:IIE,IJS:IJN,:) - 0.5*PCR(IIB:IIE+1,IJS:IJN,:) * &
-        (ZDQ(IIB-1:IIE,IJS:IJN,:) - (1.0 - 2.0*PCR(IIB:IIE+1,IJS:IJN,:)/3.0)        &
-        * ZQ6(IIB-1:IIE,IJS:IJN,:))
+!$acc loop independent collapse(3)
+  do jk = 1, iku
+    do jj = ijs, ijn
+      do ji = iib, iie + 1
+        ZFPOS(ji, jj, jk ) = ZQR(ji - 1, jj, jk ) - 0.5 * PCR(ji, jj, jk )                  &
+                            * ( ZDQ(ji - 1, jj, jk) - (1.0 - 2.0 * PCR(ji, jj, jk ) / 3.0 ) &
+                            * ZQ6(ji - 1, jj, jk) )
+      end do
+    end do
+  end do
+!
 !
 #ifndef _OPENACC
    CALL GET_HALO(ZFPOS)
@@ -926,8 +942,15 @@ ENDDO ; ENDDO ; ENDDO
 !!$   ZFNEG(IIB-1:IIE,:,:) = ZQL(IIB-1:IIE,:,:) - 0.5*PCR(IIB-1:IIE,:,:) *  &            
 !!$        (ZDQ(IIB-1:IIE,:,:) + (1.0 + 2.0*PCR(IIB-1:IIE,:,:)/3.0)         &
 !!$        * ZQ6(IIB-1:IIE,:,:))
-   ZFNEG(:,IJS:IJN,:) = ZQL(:,IJS:IJN,:) - 0.5*PCR(:,IJS:IJN,:) *      &            
-        ( ZDQ(:,IJS:IJN,:) + (1.0 + 2.0*PCR(:,IJS:IJN,:)/3.0) * ZQ6(:,IJS:IJN,:) )
+!$acc loop independent collapse(3)
+  do jk = 1, iku
+    do jj = ijs, ijn
+      do ji = 1, iiu
+   ZFNEG(ji, jj, jk ) = ZQL(ji, jj, jk ) - 0.5*PCR(ji, jj, jk ) *      &
+        ( ZDQ(ji, jj, jk ) + (1.0 + 2.0*PCR(ji, jj, jk )/3.0) * ZQ6(ji, jj, jk ) )
+      end do
+    end do
+  end do
 !
 #ifndef _OPENACC
    CALL GET_HALO(ZFNEG)
@@ -956,7 +979,11 @@ ENDDO ; ENDDO ; ENDDO
 !$acc end kernels
    CALL MXM_DEVICE(PRHO,ZQL)
 !$acc kernels
-   ZQR =  PCR* ZQL *( ZFPOS*(0.5+SIGN(0.5,PCR)) + ZFNEG*(0.5-SIGN(0.5,PCR)) )  
+   where ( PCR(:,:,:) > 0. )
+     ZQR(:,:,:) = PCR(:,:,:) * ZQL(:,:,:) * ZFPOS(:,:,:)
+   elsewhere
+     ZQR(:,:,:) = PCR(:,:,:) * ZQL(:,:,:) * ZFNEG(:,:,:)
+   end where
     !dxf(PR,ZQR)
 !$acc end kernels
    CALL DXF_DEVICE(ZQR,PR)
@@ -1108,6 +1135,7 @@ CONTAINS
 !-------------------------------------------------------------------------------
 !
 USE MODE_ll
+use mode_msg
 #ifndef _OPENACC
 USE MODI_SHUMAN
 #else
@@ -1184,6 +1212,7 @@ INTEGER                          :: IKB,IKE
 INTEGER                          :: IJN,IJS
 !JUAN ACC
 #endif
+integer :: ji, jj, jk
 !-------------------------------------------------------------------------------
 !
 !*       0.3.     COMPUTES THE DOMAIN DIMENSIONS
@@ -1218,15 +1247,22 @@ CALL GET_HALO_D(PSRC,HDIR="01_Y")
 !
 !
 !$acc kernels
-PR=PSRC
-ZQL=PSRC
-ZQR=PSRC
-ZDQ=PSRC
-ZQ6=PSRC
-ZDMQ=PSRC
-ZQL0=PSRC
-ZQR0=PSRC
-ZQ60=PSRC
+!$acc loop independent collapse(3)
+  do jk = 1, iku
+    do jj = 1, iju
+      do ji = 1, iiu
+        PR   (ji, jj, jk ) = PSRC(ji, jj, jk )
+        ZQL  (ji, jj, jk ) = PSRC(ji, jj, jk )
+        ZQR  (ji, jj, jk ) = PSRC(ji, jj, jk )
+        ZDQ  (ji, jj, jk ) = PSRC(ji, jj, jk )
+        ZQ6  (ji, jj, jk ) = PSRC(ji, jj, jk )
+        ZDMQ (ji, jj, jk ) = PSRC(ji, jj, jk )
+        ZQL0 (ji, jj, jk ) = PSRC(ji, jj, jk )
+        ZQR0 (ji, jj, jk ) = PSRC(ji, jj, jk )
+        ZQ60 (ji, jj, jk ) = PSRC(ji, jj, jk )
+    end do
+  end do
+end do
 #ifndef _OPENACC
 ZFPOS=PSRC
 ZFNEG=PSRC
@@ -1245,7 +1281,7 @@ SELECT CASE ( HLBCY(1) ) ! Y direction LBC type: (1) for left side
 !
 CASE ('CYCL','WALL')          ! In that case one must have HLBCY(1) == HLBCY(2)
 #ifdef _OPENACC
-PRINT *,'OPENACC: ppm::PPM_01_Y CYCL/WALL boundaries not yet implemented'
+  call Print_msg( NVERB_ERROR, 'GEN', 'PPM_01_Y', 'OpenACC: CYCL/WALL boundaries not yet implemented' )
 #endif
 !
 ! calculate dmq
@@ -1322,19 +1358,11 @@ PRINT *,'OPENACC: ppm::PPM_01_Y CYCL/WALL boundaries not yet implemented'
       ZQL(IIW:IIA,:,:) = PSRC(IIW:IIA,:,:)
       ZQR(IIW:IIA,:,:) = PSRC(IIW:IIA,:,:)
       ZQ6(IIW:IIA,:,:) = 0.0
-#ifndef MNH_BITREP
-   ELSEWHERE ( ZQ60(IIW:IIA,:,:)*ZDQ(IIW:IIA,:,:) < -(ZDQ(IIW:IIA,:,:))**2 )
-#else
-   ELSEWHERE ( ZQ60(IIW:IIA,:,:)*ZDQ(IIW:IIA,:,:) < -BR_P2(ZDQ(IIW:IIA,:,:)) )
-#endif
+   ELSEWHERE ( ZQ60(IIW:IIA,:,:)*ZDQ(IIW:IIA,:,:) < -ZDQ(IIW:IIA,:,:)*ZDQ(IIW:IIA,:,:) )
       ZQ6(IIW:IIA,:,:) = 3.0*(ZQL0(IIW:IIA,:,:) - PSRC(IIW:IIA,:,:))
       ZQR(IIW:IIA,:,:) = ZQL0(IIW:IIA,:,:) - ZQ6(IIW:IIA,:,:)
       ZQL(IIW:IIA,:,:) = ZQL0(IIW:IIA,:,:)
-#ifndef MNH_BITREP
-   ELSEWHERE ( ZQ60(IIW:IIA,:,:)*ZDQ(IIW:IIA,:,:) > (ZDQ(IIW:IIA,:,:))**2 )
-#else
-   ELSEWHERE ( ZQ60(IIW:IIA,:,:)*ZDQ(IIW:IIA,:,:) > BR_P2(ZDQ(IIW:IIA,:,:)) )
-#endif
+   ELSEWHERE ( ZQ60(IIW:IIA,:,:)*ZDQ(IIW:IIA,:,:) > ZDQ(IIW:IIA,:,:)*ZDQ(IIW:IIA,:,:) )
       ZQ6(IIW:IIA,:,:) = 3.0*(ZQR0(IIW:IIA,:,:) - PSRC(IIW:IIA,:,:))
       ZQL(IIW:IIA,:,:) = ZQR0(IIW:IIA,:,:) - ZQ6(IIW:IIA,:,:)
       ZQR(IIW:IIA,:,:) = ZQR0(IIW:IIA,:,:)
@@ -1371,8 +1399,7 @@ PRINT *,'OPENACC: ppm::PPM_01_Y CYCL/WALL boundaries not yet implemented'
    PR = DYF( PCR*MYM(PRHO)*( ZFPOS*(0.5+SIGN(0.5,PCR)) + & 
                              ZFNEG*(0.5-SIGN(0.5,PCR)) ) )
 #else
-PRINT *,'not yet implemented'
-STOP
+  call Print_msg( NVERB_ERROR, 'GEN', 'PPM_01_Y', 'OpenACC: CYCL/WALL boundaries not yet implemented' )
 #endif
    CALL GET_HALO(PR) 
 !
@@ -1404,12 +1431,19 @@ CASE('OPEN')
    ENDIF
 !
 ! monotonize the difference followinq eq. 5 in Lin94
-   ZDMQ(IIW:IIA,IJB:IJE,:) = &
-        SIGN( (MIN( ABS(ZDMQ(IIW:IIA,IJB:IJE,:)),2.0*(PSRC(IIW:IIA,IJB:IJE,:) - &
-        MIN(PSRC(IIW:IIA,IJB-1:IJE-1,:),PSRC(IIW:IIA,IJB:IJE,:),PSRC(IIW:IIA,IJB+1:IJE+1,:))),    &
-        2.0*(MAX(PSRC(IIW:IIA,IJB-1:IJE-1,:),PSRC(IIW:IIA,IJB:IJE,:),PSRC(IIW:IIA,IJB+1:IJE+1,:)) - &
-        PSRC(IIW:IIA,IJB:IJE,:)) )), ZDMQ(IIW:IIA,IJB:IJE,:) )
-!!$   ZDMQ(:,IJB-1,:) = & 
+!$acc loop independent collapse(3)
+  do jk = 1, iku
+    do jj = ijb, ije
+      do ji = iiw, iia
+        ZDMQ(ji, jj, jk ) = SIGN( &
+        MIN( ABS( ZDMQ(ji, jj, jk ) ),                                                                             &
+             2.0 * (   PSRC(ji, jj, jk ) - MIN( PSRC(ji, jj-1, jk ), PSRC(ji, jj, jk ), PSRC(ji, jj+1, jk ) ) ),   &
+             2.0 * ( - PSRC(ji, jj, jk ) + MAX( PSRC(ji, jj-1, jk ), PSRC(ji, jj, jk ), PSRC(ji, jj+1, jk ) ) ) ), &
+        ZDMQ(ji, jj, jk ) )
+    end do
+  end do
+end do
+!!$   ZDMQ(:,IJB-1,:) = &
 !!$        SIGN( (MIN( ABS(ZDMQ(:,IJB-1,:)), 2.0*(PSRC(:,IJB-1,:) - &
 !!$        MIN(PSRC(:,IJE-1,:),PSRC(:,IJB-1,:),PSRC(:,IJB,:))),   &
 !!$        2.0*(MAX(PSRC(:,IJE-1,:),PSRC(:,IJB-1,:),PSRC(:,IJB,:)) - &
@@ -1432,8 +1466,14 @@ CASE('OPEN')
 !
 ! calculate qL and qR with the modified dmq
 !
-   ZQL0(IIW:IIA,IJB:IJE+1,:) = 0.5*(PSRC(IIW:IIA,IJB:IJE+1,:) + PSRC(IIW:IIA,IJB-1:IJE,:)) - &
-        (ZDMQ(IIW:IIA,IJB:IJE+1,:) - ZDMQ(IIW:IIA,IJB-1:IJE,:))/6.0
+!$acc loop independent collapse(3)
+  do jk = 1, iku
+    do jj = ijb, ije + 1
+      do ji = iiw, iia
+        ZQL0(ji, jj, jk ) = 0.5 * ( PSRC(ji, jj, jk ) + PSRC(ji, jj-1, jk )) - ( ZDMQ(ji, jj, jk ) - ZDMQ(ji, jj-1, jk ) ) / 6.0
+      end do
+    end do
+  end do
 !
 #ifndef _OPENACC
    CALL  GET_HALO(ZQL0)
@@ -1449,7 +1489,14 @@ CALL  GET_HALO_D(ZQL0,HDIR="01_Y")
     ZQL0(IIW:IIA,IJB-1,:) = ZQL0(IIW:IIA,IJB,:)
    ENDIF
 !
-   ZQR0(IIW:IIA,IJB-1:IJE,:) = ZQL0(IIW:IIA,IJB:IJE+1,:)
+!$acc loop independent collapse(3)
+  do jk = 1, iku
+    do jj = ijb - 1, ije
+      do ji = iiw, iia
+        ZQR0(ji, jj, jk ) = ZQL0(ji, jj+1, jk )
+      end do
+    end do
+  end do
 !
 !  NORTH BOUND
 !
@@ -1475,19 +1522,11 @@ CALL  GET_HALO_D(ZQL0,HDIR="01_Y")
       ZQL(IIW:IIA,:,:) = PSRC(IIW:IIA,:,:)
       ZQR(IIW:IIA,:,:) = PSRC(IIW:IIA,:,:)
       ZQ6(IIW:IIA,:,:) = 0.0
-#ifndef MNH_BITREP
-   ELSEWHERE ( ZQ60(IIW:IIA,:,:)*ZDQ(IIW:IIA,:,:) < -(ZDQ(IIW:IIA,:,:))**2 )
-#else
-   ELSEWHERE ( ZQ60(IIW:IIA,:,:)*ZDQ(IIW:IIA,:,:) < -BR_P2(ZDQ(IIW:IIA,:,:)) )
-#endif
+   ELSEWHERE ( ZQ60(IIW:IIA,:,:)*ZDQ(IIW:IIA,:,:) < -ZDQ(IIW:IIA,:,:)*ZDQ(IIW:IIA,:,:) )
       ZQ6(IIW:IIA,:,:) = 3.0*(ZQL0(IIW:IIA,:,:) - PSRC(IIW:IIA,:,:))
       ZQR(IIW:IIA,:,:) = ZQL0(IIW:IIA,:,:) - ZQ6(IIW:IIA,:,:)
       ZQL(IIW:IIA,:,:) = ZQL0(IIW:IIA,:,:)
-#ifndef MNH_BITREP
-   ELSEWHERE ( ZQ60(IIW:IIA,:,:)*ZDQ(IIW:IIA,:,:) > (ZDQ(IIW:IIA,:,:))**2 )
-#else
-   ELSEWHERE ( ZQ60(IIW:IIA,:,:)*ZDQ(IIW:IIA,:,:) > BR_P2(ZDQ(IIW:IIA,:,:)) )
-#endif
+   ELSEWHERE ( ZQ60(IIW:IIA,:,:)*ZDQ(IIW:IIA,:,:) > ZDQ(IIW:IIA,:,:)*ZDQ(IIW:IIA,:,:) )
       ZQ6(IIW:IIA,:,:) = 3.0*(ZQR0(IIW:IIA,:,:) - PSRC(IIW:IIA,:,:))
       ZQL(IIW:IIA,:,:) = ZQR0(IIW:IIA,:,:) - ZQ6(IIW:IIA,:,:)
       ZQR(IIW:IIA,:,:) = ZQR0(IIW:IIA,:,:)
@@ -1498,12 +1537,10 @@ CALL  GET_HALO_D(ZQL0,HDIR="01_Y")
    ZDQ(IIW:IIA,:,:) = ZQR(IIW:IIA,:,:) - ZQL(IIW:IIA,:,:)
 !
 #else
-   !$acc loop gang vector
+!$acc loop independent collapse(3)
    DO K=IKB,IKE
-      !$acc loop gang vector
       DO J=IJS,IJN
-         !$acc loop gang vector
-         DO I=1,IIU 
+         DO I=1,IIU
             !
             ! determine initial coefficients of the parabolae
             !
@@ -1522,21 +1559,11 @@ CALL  GET_HALO_D(ZQL0,HDIR="01_Y")
                ZQL(I,J,K) = PSRC(I,J,K)
                ZQR(I,J,K) = PSRC(I,J,K)
                ZQ6(I,J,K) = 0.0
-            END IF
-#ifndef MNH_BITREP
-            IF ( ( ZDMQ(I,J,K) /= 0.0 ) .AND. ( ZQ60(I,J,K)*ZDQ(I,J,K) < -(ZDQ(I,J,K))**2 ) ) THEN
-#else
-            IF ( ( ZDMQ(I,J,K) /= 0.0 ) .AND. ( ZQ60(I,J,K)*ZDQ(I,J,K) < -BR_P2(ZDQ(I,J,K)) ) ) THEN
-#endif
+            ELSE IF ( ZQ60(I,J,K)*ZDQ(I,J,K) < -ZDQ(I,J,K)*ZDQ(I,J,K) ) THEN
                ZQ6(I,J,K) = 3.0*(ZQL0(I,J,K) - PSRC(I,J,K))
                ZQR(I,J,K) = ZQL0(I,J,K) - ZQ6(I,J,K)
                ZQL(I,J,K) = ZQL0(I,J,K)
-            END IF
-#ifndef MNH_BITREP
-            IF ( ( ZDMQ(I,J,K) /= 0.0 ) .AND. ( ZQ60(I,J,K)*ZDQ(I,J,K) > (ZDQ(I,J,K))**2 ) ) THEN
-#else
-            IF ( ( ZDMQ(I,J,K) /= 0.0 ) .AND. ( ZQ60(I,J,K)*ZDQ(I,J,K) > BR_P2(ZDQ(I,J,K)) ) ) THEN
-#endif
+            ELSE IF ( ZQ60(I,J,K)*ZDQ(I,J,K) > ZDQ(I,J,K)*ZDQ(I,J,K) ) THEN
                ZQ6(I,J,K) = 3.0*(ZQR0(I,J,K) - PSRC(I,J,K))
                ZQL(I,J,K) = ZQR0(I,J,K) - ZQ6(I,J,K)
                ZQR(I,J,K) = ZQR0(I,J,K)
@@ -1556,9 +1583,15 @@ CALL  GET_HALO_D(ZQL0,HDIR="01_Y")
 !!$   ZFPOS(:,IJB+1:IJE+1,:) = ZQR(:,IJB:IJE,:) - 0.5*PCR(:,IJB+1:IJE+1,:) * &            
 !!$        (ZDQ(:,IJB:IJE,:) - (1.0 - 2.0*PCR(:,IJB+1:IJE+1,:)/3.0)        &
 !!$        * ZQ6(:,IJB:IJE,:))
-   ZFPOS(IIW:IIA,IJB:IJE+1,:) = ZQR(IIW:IIA,IJB-1:IJE,:) - 0.5*PCR(IIW:IIA,IJB:IJE+1,:) * &            
-        (ZDQ(IIW:IIA,IJB-1:IJE,:) - (1.0 - 2.0*PCR(IIW:IIA,IJB:IJE+1,:)/3.0)        &
-        * ZQ6(IIW:IIA,IJB-1:IJE,:))
+!$acc loop independent collapse(3)
+  do jk = 1, iku
+    do jj = ijb, ije + 1
+      do ji = iiw, iia
+        ZFPOS(ji, jj, jk ) = ZQR(ji, jj-1, jk ) - 0.5 * PCR(ji, jj, jk )   &
+                                      * ( ZDQ(ji, jj-1, jk ) - ( 1.0 - 2.0 * PCR(ji, jj, jk ) / 3.0 ) * ZQ6(ji, jj-1, jk ) )
+      end do
+    end do
+  end do
 !
 #ifndef _OPENACC
    CALL GET_HALO(ZFPOS)
@@ -1587,15 +1620,22 @@ CALL GET_HALO_D(ZFPOS,HDIR="01_Y")
 !!$   ZFNEG(:,IJB-1:IJE,:) = ZQL(:,IJB-1:IJE,:) - 0.5*PCR(:,IJB-1:IJE,:) * &            
 !!$        ( ZDQ(:,IJB-1:IJE,:) + (1.0 + 2.0*PCR(:,IJB-1:IJE,:)/3.0) * &
 !!$        ZQ6(:,IJB-1:IJE,:) )
-   ZFNEG(IIW:IIA,:,:) = ZQL(IIW:IIA,:,:) - 0.5*PCR(IIW:IIA,:,:) *      &            
-        ( ZDQ(IIW:IIA,:,:) + (1.0 + 2.0*PCR(IIW:IIA,:,:)/3.0) * ZQ6(IIW:IIA,:,:) )
+!$acc loop independent collapse(3)
+  do jk = 1, iku
+    do jj = 1, iju
+      do ji = iiw, iia
+        ZFNEG(ji, jj, jk ) = ZQL(ji, jj, jk ) - 0.5 * PCR(ji, jj, jk )      &
+                                     * ( ZDQ(ji, jj, jk ) + ( 1.0 + 2.0 * PCR(ji, jj, jk ) / 3.0 ) * ZQ6(ji, jj, jk ) )
+      end do
+    end do
+  end do
 !
 #ifndef _OPENACC
    CALL GET_HALO(ZFNEG)
 #else
-!$acc end kernels   
+!$acc end kernels
    CALL GET_HALO_D(ZFNEG,HDIR="01_Y")
-!$acc kernels   
+!$acc kernels
 #endif
 !
 ! advection flux at open boundary when u(IJE+1) < 0
@@ -1614,19 +1654,22 @@ CALL GET_HALO_D(ZFPOS,HDIR="01_Y")
                              ZFNEG*(0.5-SIGN(0.5,PCR)) ) )
 !
 #else
-!$acc end kernels 
+!$acc end kernels
    CALL MYM_DEVICE(PRHO,ZQL)
-!$acc kernels 
-   DO K=IKB,IKE 
+!$acc kernels
+!$acc loop independent collapse(3)
+   DO K=IKB,IKE
       DO J=IJS,IJN
-         DO I=1,IIU 
-            ZQR(I,J,K) =  PCR(I,J,K)* ZQL(I,J,K) &
-                 * ( ZFPOS(I,J,K)*(0.5+SIGN(0.5,PCR(I,J,K))) &
-                 + ZFNEG(I,J,K)*(0.5-SIGN(0.5,PCR(I,J,K))) ) 
+         DO I=1,IIU
+          if ( PCR(I,J,K) > 0. ) then
+            ZQR(I,J,K) =  PCR(I,J,K)* ZQL(I,J,K) * ZFPOS(I,J,K)
+          else
+            ZQR(I,J,K) =  PCR(I,J,K)* ZQL(I,J,K) * ZFNEG(I,J,K)
+          end if
          END DO
       END DO
    END DO
-!$acc end kernels 
+!$acc end kernels
    CALL DYF_DEVICE(ZQR,PR)
 #endif
 !
@@ -1855,6 +1898,7 @@ INTEGER                          :: IIU,IJU,IKU
 INTEGER                          :: I,J,K 
 !JUAN ACC
 #endif
+integer                          :: ji, jj, jk
 !
 !-------------------------------------------------------------------------------
 !
@@ -1876,19 +1920,26 @@ IKU = SIZE(PSRC,3)
 #ifndef _OPENACC
 ZDMQ = DIF2Z(PSRC)
 #else
-!$acc kernels  
+!$acc kernels
   dif2z(ZDMQ,PSRC)
 #endif
 !
 ! monotonize the difference followinq eq. 5 in Lin94
 ! use the periodic BC here, it doesn't matter for vertical (hopefully) 
 !
-ZDMQ(:,:,IKB:IKE) = &
-     SIGN( (MIN( ABS(ZDMQ(:,:,IKB:IKE)),2.0*(PSRC(:,:,IKB:IKE) - &
-     MIN(PSRC(:,:,IKB-1:IKE-1),PSRC(:,:,IKB:IKE),PSRC(:,:,IKB+1:IKE+1))),    &
-     2.0*(MAX(PSRC(:,:,IKB-1:IKE-1),PSRC(:,:,IKB:IKE),PSRC(:,:,IKB+1:IKE+1)) - &
-     PSRC(:,:,IKB:IKE)) )), ZDMQ(:,:,IKB:IKE) )
-ZDMQ(:,:,IKB-1) = & 
+!$acc loop independent collapse(3)
+do jk = ikb, ike
+  do jj = 1, iju
+    do ji = 1, iiu
+      ZDMQ(ji, jj, jk ) = SIGN(                                                                                    &
+        MIN( ABS( ZDMQ(ji, jj, jk ) ),                                                                             &
+             2.0 * (   PSRC(ji, jj, jk ) - MIN( PSRC(ji, jj, jk-1 ), PSRC(ji, jj, jk ), PSRC(ji, jj, jk+1 ) ) ) ,  &
+             2.0 * ( - PSRC(ji, jj, jk ) + MAX( PSRC(ji, jj, jk-1 ), PSRC(ji, jj, jk ), PSRC(ji, jj, jk+1 ) ) ) ), &
+        ZDMQ(ji, jj, jk ) )
+    end do
+  end do
+end do
+ZDMQ(:,:,IKB-1) = &
      SIGN( (MIN( ABS(ZDMQ(:,:,IKB-1)), 2.0*(PSRC(:,:,IKB-1) - &
      MIN(PSRC(:,:,IKE-1),PSRC(:,:,IKB-1),PSRC(:,:,IKB))),   &
      2.0*(MAX(PSRC(:,:,IKE-1),PSRC(:,:,IKB-1),PSRC(:,:,IKB)) - &
@@ -1901,17 +1952,37 @@ ZDMQ(:,:,IKE+1) = &
 !
 ! calculate qL and qR with the modified dmq
 !
-ZQL0(:,:,IKB:IKE+1) = 0.5*(PSRC(:,:,IKB:IKE+1) + PSRC(:,:,IKB-1:IKE)) - &
-     (ZDMQ(:,:,IKB:IKE+1) - ZDMQ(:,:,IKB-1:IKE))/6.0
+!$acc loop independent collapse(3)
+do jk = ikb, ike + 1
+  do jj = 1, iju
+    do ji = 1, iiu
+      ZQL0(ji, jj, jk ) = 0.5 * ( PSRC(ji, jj, jk ) + PSRC(ji, jj, jk-1 ) ) - ( ZDMQ(ji, jj, jk ) - ZDMQ(ji, jj, jk-1 ) ) / 6.0
+    end do
+  end do
+end do
 ZQL0(:,:,IKB-1) = ZQL0(:,:,IKE)
 !
-ZQR0(:,:,IKB-1:IKE) = ZQL0(:,:,IKB:IKE+1)
+!$acc loop independent collapse(3)
+do jk = ikb - 1, ike
+  do jj = 1, iju
+    do ji = 1, iiu
+      ZQR0(ji, jj, jk ) = ZQL0(ji, jj, jk+1 )
+    end do
+  end do
+end do
 ZQR0(:,:,IKE+1) = ZQR0(:,:,IKB)
 !
 ! determine initial coefficients of the parabolae
 !
-ZDQ = ZQR0 - ZQL0
-ZQ60 = 6.0*(PSRC - 0.5*(ZQL0 + ZQR0))
+!$acc loop independent collapse(3)
+do jk = ikb - 1, ike
+  do jj = 1, iju
+    do ji = 1, iiu
+      ZDQ (ji, jj, jk ) = ZQR0(ji, jj, jk ) - ZQL0(ji, jj, jk )
+      ZQ60(ji, jj, jk ) = 6.0 * ( PSRC(ji, jj, jk ) - 0.5 * ( ZQL0(ji, jj, jk ) + ZQR0(ji, jj, jk ) ) )
+    end do
+  end do
+end do
 #ifndef _OPENACC
 !
 ! initialize final parabolae parameters
@@ -1989,12 +2060,10 @@ ZDQ = ZQR - ZQL
 !!! recalculate coefficients of the parabolae
 !!!
 !!ZDQ = ZQR - ZQL
-   !$acc loop gang vector
+!$acc loop independent collapse(3)
    DO K=1,IKU
-      !$acc loop gang vector
       DO J=1,IJU
-         !$acc loop gang vector
-         DO I=1,IIU 
+         DO I=1,IIU
             !
             ! determine initial coefficients of the parabolae
             !
@@ -2013,21 +2082,11 @@ ZDQ = ZQR - ZQL
                ZQL(I,J,K) = PSRC(I,J,K)
                ZQR(I,J,K) = PSRC(I,J,K)
                ZQ6(I,J,K) = 0.0
-            END IF
-#ifndef MNH_BITREP
-            IF ( ( ZDMQ(I,J,K) /= 0.0 ) .AND. ( ZQ60(I,J,K)*ZDQ(I,J,K) < -ZDQ(I,J,K)**2 ) ) THEN
-#else
-            IF ( ( ZDMQ(I,J,K) /= 0.0 ) .AND. ( ZQ60(I,J,K)*ZDQ(I,J,K) < -BR_P2(ZDQ(I,J,K)) ) ) THEN
-#endif
+            ELSE IF ( ZQ60(I,J,K)*ZDQ(I,J,K) < -ZDQ(I,J,K)*ZDQ(I,J,K) ) THEN
                ZQ6(I,J,K) = 3.0*(ZQL0(I,J,K) - PSRC(I,J,K))
                ZQR(I,J,K) = ZQL0(I,J,K) - ZQ6(I,J,K)
                ZQL(I,J,K) = ZQL0(I,J,K)
-            END IF
-#ifndef MNH_BITREP
-            IF ( ( ZDMQ(I,J,K) /= 0.0 ) .AND. ( ZQ60(I,J,K)*ZDQ(I,J,K) > ZDQ(I,J,K)**2 ) ) THEN
-#else
-            IF ( ( ZDMQ(I,J,K) /= 0.0 ) .AND. ( ZQ60(I,J,K)*ZDQ(I,J,K) > BR_P2(ZDQ(I,J,K)) ) ) THEN
-#endif
+            ELSE IF ( ZQ60(I,J,K)*ZDQ(I,J,K) > ZDQ(I,J,K)*ZDQ(I,J,K) ) THEN
                ZQ6(I,J,K) = 3.0*(ZQR0(I,J,K) - PSRC(I,J,K))
                ZQL(I,J,K) = ZQR0(I,J,K) - ZQ6(I,J,K)
                ZQR(I,J,K) = ZQR0(I,J,K)
@@ -2044,9 +2103,15 @@ ZDQ = ZQR - ZQL
 !
 ! and finally calculate fluxes for the advection
 !
-ZFPOS(:,:,IKB+1:IKE+1) = ZQR(:,:,IKB:IKE) - 0.5*PCR(:,:,IKB+1:IKE+1) * &            
-     (ZDQ(:,:,IKB:IKE) - (1.0 - 2.0*PCR(:,:,IKB+1:IKE+1)/3.0)        &
-     * ZQ6(:,:,IKB:IKE))
+!$acc loop independent collapse(3)
+do jk = ikb + 1, ike + 1
+  do jj = 1, iju
+    do ji = 1, iiu
+      ZFPOS(ji, jj, jk ) = ZQR(ji, jj, jk-1 ) - 0.5 * PCR(ji, jj, jk ) &
+                                   * ( ZDQ(ji, jj, jk-1 ) - ( 1.0 - 2.0 * PCR(ji, jj, jk ) / 3.0) * ZQ6(ji, jj, jk-1 ) )
+    end do
+  end do
+end do
 !
 ! advection flux at open boundary when u(IKB) > 0
 ZFPOS(:,:,IKB) = (PSRC(:,:,IKB-1) - ZQR(:,:,IKB-1))*PCR(:,:,IKB) + &
@@ -2056,9 +2121,15 @@ ZFPOS(:,:,IKB) = (PSRC(:,:,IKB-1) - ZQR(:,:,IKB-1))*PCR(:,:,IKB) + &
 ! we set it to 0
 ZFPOS(:,:,IKB-1) = 0.0
 !
-ZFNEG(:,:,IKB-1:IKE) = ZQL(:,:,IKB-1:IKE) - 0.5*PCR(:,:,IKB-1:IKE) *      &            
-     ( ZDQ(:,:,IKB-1:IKE) + (1.0 + 2.0*PCR(:,:,IKB-1:IKE)/3.0) * &
-       ZQ6(:,:,IKB-1:IKE) )
+!$acc loop independent collapse(3)
+do jk = ikb - 1, ike
+  do jj = 1, iju
+    do ji = 1, iiu
+      ZFNEG(ji, jj, jk ) = ZQL(ji, jj, jk ) - 0.5 * PCR(ji, jj, jk ) &
+                                   * ( ZDQ(ji, jj, jk ) + ( 1.0 + 2.0 * PCR(ji, jj, jk ) / 3.0) * ZQ6(ji, jj, jk ) )
+    end do
+  end do
+end do
 !
 ! advection flux at open boundary when u(IKE+1) < 0
 ZFNEG(:,:,IKE+1) = (ZQR(:,:,IKE)-PSRC(:,:,IKE+1))*PCR(:,:,IKE+1) + &
@@ -2073,7 +2144,18 @@ PR = DZF(1,IKU,1, PCR*MZM(1,IKU,1,PRHO)*( ZFPOS*(0.5+SIGN(0.5,PCR)) + &
 !$acc end kernels 
     CALL MZM_DEVICE(PRHO,ZQL)
 !$acc kernels 
-    ZQR =  PCR* ZQL*( ZFPOS*(0.5+SIGN(0.5,PCR)) + ZFNEG*(0.5-SIGN(0.5,PCR)) ) 
+!$acc loop independent collapse(3)
+do jk = 1, iku
+  do jj = 1, iju
+    do ji = 1, iiu
+      if ( PCR(ji, jj, jk ) > 0. ) then
+        ZQR(ji, jj, jk ) =  PCR(ji, jj, jk ) * ZQL(ji, jj, jk ) * ZFPOS(ji, jj, jk )
+      else
+        ZQR(ji, jj, jk ) =  PCR(ji, jj, jk ) * ZQL(ji, jj, jk ) * ZFNEG(ji, jj, jk )
+      end if
+    end do
+  end do
+end do
     !dzf(PR,ZQR)
 !$acc end kernels
     CALL DZF_DEVICE(1,1,1,ZQR,PR)
@@ -2228,6 +2310,7 @@ CONTAINS
 !-------------------------------------------------------------------------------
 !
 USE MODE_ll
+use mode_msg
 #ifndef _OPENACC
 USE MODI_SHUMAN
 #else
@@ -2352,7 +2435,7 @@ ZPHAT(IIB+1:IIE,IJS:IJN,:) = ( 7.0 * &
 SELECT CASE ( HLBCX(1) ) ! X direction LBC type: (1) for left side
 CASE ('CYCL','WALL')            ! In that case one must have HLBCX(1) == HLBCX(2)
 #ifdef _OPENACC
-PRINT *,'OPENACC: ppm::PPM_S0_X CYCL/WALL boundaries not yet implemented'
+  call Print_msg( NVERB_ERROR, 'GEN', 'PPM_S0_X', 'OpenACC: CYCL/WALL boundaries not yet implemented' )
 #endif
 !
 !!$   ZPHAT(IIB,:,:) = (7.0 * &
@@ -2410,8 +2493,7 @@ CALL GET_HALO(ZFNEG) ! JUAN
         DXF( PCR*MXM(PRHO)*( ZFPOS*(0.5+SIGN(0.5,PCR)) + & 
                              ZFNEG*(0.5-SIGN(0.5,PCR)) ) )
 #else
-PRINT *,'not yet implemented'
-STOP
+  call Print_msg( NVERB_ERROR, 'GEN', 'PPM_S0_X', 'OpenACC: CYCL/WALL boundaries not yet implemented' )
 #endif
    CALL GET_HALO(PR) ! JUAN
 !
@@ -2645,6 +2727,7 @@ CONTAINS
 !-------------------------------------------------------------------------------
 !
 USE MODE_ll
+use mode_msg
 #ifndef _OPENACC
 USE MODI_SHUMAN
 #else
@@ -2767,7 +2850,7 @@ ZPHAT(IIW:IIA,IJB+1:IJE,:) = (7.0 * &
 SELECT CASE ( HLBCY(1) ) ! Y direction LBC type: (1) for left side
 CASE ('CYCL','WALL')            ! In that case one must have HLBCY(1) == HLBCY(2)
 #ifdef _OPENACC
-PRINT *,'OPENACC: ppm::PPM_S0_Y CYCL/WALL boundaries not yet implemented'
+  call Print_msg( NVERB_ERROR, 'GEN', 'PPM_S0_Y', 'OpenACC: CYCL/WALL boundaries not yet implemented' )
 #endif
 !
 !!$   ZPHAT(:,IJB,:) = (7.0 * &
@@ -2826,8 +2909,7 @@ CALL GET_HALO(ZFNEG) ! JUAN
         DYF( PCR*MYM(PRHO)*( ZFPOS*(0.5+SIGN(0.5,PCR)) + & 
                              ZFNEG*(0.5-SIGN(0.5,PCR)) ) )
 #else
-PRINT *,'not yet implemented'
-STOP
+  call Print_msg( NVERB_ERROR, 'GEN', 'PPM_S0_Y', 'OpenACC: CYCL/WALL boundaries not yet implemented' )
 #endif
 !
 CASE ('OPEN')
@@ -3239,6 +3321,7 @@ END FUNCTION PPM_S0_Z
                           PTSTEP, PR)
 !     ########################################################################
 USE MODE_ll
+use mode_msg
 USE MODE_IO
 #ifndef _OPENACC
 USE MODI_SHUMAN
@@ -3278,8 +3361,7 @@ INTEGER :: IZPHAT,IZRUT,IZFUP,IZFCOR,IZRPOS,IZRNEG
 
 
 #ifdef _OPENACC
-PRINT *,'OPENACC: ppm::PPM_S1_X not yet implemented'
-CALL ABORT
+  call Print_msg( NVERB_ERROR, 'GEN', 'PPM_S1_X', 'OpenACC: not yet implemented' )
 #endif
 
     CALL  MNH_GET_ZT3D(IZPHAT,IZRUT,IZFUP,IZFCOR,IZRPOS,IZRNEG)
@@ -3585,6 +3667,7 @@ END FUNCTION PPM_S1_X
 !     ########################################################################
 USE MODE_ll
 USE MODE_IO
+use mode_msg
 #ifndef _OPENACC
 USE MODI_SHUMAN
 #else
@@ -3623,8 +3706,7 @@ INTEGER :: IZPHAT,IZRVT,IZFUP,IZFCOR,IZRPOS,IZRNEG
 
 
 #ifdef _OPENACC
-PRINT *,'OPENACC: ppm::PPM_S1_Y not yet implemented'
-CALL ABORT
+  call Print_msg( NVERB_ERROR, 'GEN', 'PPM_S1_Y', 'OpenACC: not yet implemented' )
 #endif
 
     CALL  MNH_GET_ZT3D(IZPHAT,IZRVT,IZFUP,IZFCOR,IZRPOS,IZRNEG)
@@ -3937,6 +4019,7 @@ END FUNCTION PPM_S1_Y
 !     ########################################################################
 USE MODE_ll
 USE MODE_IO
+use mode_msg
 #ifndef _OPENACC
 USE MODI_SHUMAN
 #else
@@ -3973,8 +4056,7 @@ INTEGER :: IZPHAT,IZRVT,IZFUP,IZFCOR,IZRPOS,IZRNEG
 
 
 #ifdef _OPENACC
-PRINT *,'OPENACC: ppm::PPM_S1_Z not yet implemented'
-CALL ABORT
+  call Print_msg( NVERB_ERROR, 'GEN', 'PPM_S1_Z', 'OpenACC: not yet implemented' )
 #endif
 
     CALL  MNH_GET_ZT3D(IZPHAT,IZRVT,IZFUP,IZFCOR,IZRPOS,IZRNEG)
-- 
GitLab