Philippe 17/09/2021: OpenACC: improve Source_neg_correct + bit-reproducibility

src/MNH/sources_neg_correct.f90

@@ -32,9 +32,16 @@ use modd_param_lima, only: lcold_lima => lcold, lrain_lima => lrain, lspro_lima
                            xctmin_lima => xctmin, xrtmin_lima => xrtmin
 
 use mode_budget,         only: Budget_store_init, Budget_store_end
+#ifdef MNH_OPENACC
+use mode_mnh_zwork,  only: Mnh_allocate_zt3d , Mnh_rel_zt3d
+#endif
 use mode_mppdb
 use mode_msg
 
+#ifdef MNH_BITREP
+use modi_bitrep
+#endif
+
 implicit none
 
 character(len=*),            intent(in)           :: hcloud   ! Kind of cloud parameterization
@@ -49,15 +56,21 @@ real, dimension(:, :, :, :), intent(inout)        :: prrs     ! Source terms
 real, dimension(:, :, :, :), intent(inout)        :: prsvs    ! Source terms
 real, dimension(:, :, :),    intent(in), optional :: prhodj   ! Dry density * jacobian
 
+integer :: jiu, jju, jku
 integer :: ji, jj, jk
 integer :: jr
 integer :: jrmax
 integer :: jsv
 integer :: isv_lima_end
-real, dimension(:, :, :), allocatable :: zt, zexn, zlv, zls, zcph, zcor
+real, dimension(:, :, :), pointer, contiguous :: zt, zexn, zlv, zls, zcph, zcor
+#ifdef MNH_OPENACC
+integer :: izt, izexn, izlv, izls, izcph, izcor
+#endif
 
 if ( krr == 0 ) return
 
+zcor => null()
+
 !$acc data present( ppabst, ptht, prt, prths, prrs, prsvs, prhodj )
 
 if ( mppdb_initialized ) then
@@ -141,18 +154,48 @@ else !NECON + NEGA
   end if
 end if
 
-allocate( zt  ( Size( prths, 1 ), Size( prths, 2 ), Size( prths, 3 ) ) )
-allocate( zexn( Size( prths, 1 ), Size( prths, 2 ), Size( prths, 3 ) ) )
-allocate( zlv ( Size( prths, 1 ), Size( prths, 2 ), Size( prths, 3 ) ) )
-allocate( zcph( Size( prths, 1 ), Size( prths, 2 ), Size( prths, 3 ) ) )
+jiu = Size(prths, 1 )
+jju = Size(prths, 2 )
+jku = Size(prths, 3 )
+
+#ifndef MNH_OPENACC
+allocate( zt  ( jiu, jju, jku ) )
+allocate( zexn( jiu, jju, jku ) )
+allocate( zlv ( jiu, jju, jku ) )
+allocate( zcph( jiu, jju, jku ) )
 if ( hcloud == 'ICE3' .or. hcloud == 'ICE4' .or. hcloud == 'LIMA' ) then
-  allocate( zls( Size( prths, 1 ), Size( prths, 2 ), Size( prths, 3 ) ) )
+  allocate( zls( jiu, jju, jku ) )
+  if ( krr > 3 ) then
+    allocate( zcor( jiu, jju, jku ) )
+  end if
 end if
+if ( .not. Associated( zcor ) ) Allocate( zcor(0, 0, 0) )
+#else
+izt   = Mnh_allocate_zt3d( zt,   jiu, jju, jku )
+izexn = Mnh_allocate_zt3d( zexn, jiu, jju, jku )
+izlv  = Mnh_allocate_zt3d( zlv,  jiu, jju, jku )
+izcph = Mnh_allocate_zt3d( zcph, jiu, jju, jku )
+if ( hcloud == 'ICE3' .or. hcloud == 'ICE4' .or. hcloud == 'LIMA' ) then
+  izls = Mnh_allocate_zt3d( zls, jiu, jju, jku )
+  if ( krr > 3 ) then
+    izcor = Mnh_allocate_zt3d( zcor, jiu, jju, jku )
+  else
+    izcor = Mnh_allocate_zt3d( zcor, 0, 0, 0 )
+  end if
+else
+  izls  = Mnh_allocate_zt3d( zls,  0, 0, 0 )
+  izcor = Mnh_allocate_zt3d( zcor, 0, 0, 0 )
+end if
+#endif
 
-!$acc data create( zt, zexn, zlv, zcph, zls )
+!$acc data create( zt, zexn, zlv, zcph, zls, zcor )
 
 !$acc kernels
+#ifndef MNH_BITREP
 zexn(:, :, :) = ( ppabst(:, :, :) / xp00 ) ** (xrd / xcpd )
+#else
+zexn(:, :, :) = Br_pow( ppabst(:, :, :) / xp00,  xrd / xcpd )
+#endif
 zt  (:, :, :) = ptht(:, :, :) * zexn(:, :, :)
 zlv (:, :, :) = xlvtt + ( xcpv - xcl ) * ( zt(:, :, :) - xtt )
 !$acc end kernels
@@ -165,7 +208,9 @@ end if
 zcph(:, :, :) = xcpd + xcpv * prt(:, :, :, 1)
 !$acc end kernels
 
+#ifndef MNH_OPENACC
 deallocate( zt )
+#endif
 
 CLOUD: select case ( hcloud )
   case ( 'KESS' )
@@ -240,8 +285,7 @@ CLOUD: select case ( hcloud )
 !$acc end kernels
 !   ice
     if ( krr > 3 ) then
-      allocate( zcor( Size( prths, 1 ), Size( prths, 2 ), Size( prths, 3 ) ) )
-!$acc kernels create( zcor )
+!$acc kernels
       where ( prrs(:, :, :, 1) < 0. .and. prrs(:, :, :, 4) > 0. )
         zcor(:, :, :) = Min( -prrs(:, :, :, 1), prrs(:, :, :, 4) )
         prrs(:, :, :, 1) = prrs(:, :, :, 1) + zcor(:, :, :)
@@ -250,7 +294,6 @@ CLOUD: select case ( hcloud )
         prrs(:, :, :, 4) = prrs(:, :, :, 4) - zcor(:, :, :)
       end where
 !$acc end kernels
-      deallocate( zcor )
     end if
 !
 !
@@ -338,8 +381,7 @@ CLOUD: select case ( hcloud )
         end do
       end if
       if(krr > 3) then
-        allocate( zcor( Size( prths, 1 ), Size( prths, 2 ), Size( prths, 3 ) ) )
-!$acc kernels create( zcor )
+!$acc kernels
         where ( prrs(:, :, :, 1) < 0. .and. prrs(:, :, :, 4) > 0. )
           zcor(:, :, :) = Min( -prrs(:, :, :, 1), prrs(:, :, :, 4) )
           prrs(:, :, :, 1) = prrs(:, :, :, 1) + zcor(:, :, :)
@@ -348,7 +390,6 @@ CLOUD: select case ( hcloud )
           prrs(:, :, :, 4) = prrs(:, :, :, 4) - zcor(:, :, :)
         end where
 !$acc end kernels
-        deallocate( zcor )
       end if
     end if
 
@@ -360,6 +401,10 @@ end select CLOUD
 
 !$acc end data
 
+#ifdef MNH_OPENACC
+call Mnh_rel_zt3d ( izt, izexn, izlv, izcph, izls, izcor )
+#endif
+
 if ( hbudname /= 'NECON' .and. hbudname /= 'NEGA' ) then
   if ( hcloud == 'KESS' .or. hcloud == 'ICE3' .or. hcloud == 'ICE4' .or. &
        hcloud == 'KHKO' .or. hcloud == 'C2R2' .or. hcloud == 'LIMA' ) then