From 6f727c3cfc3bbd4303e67512775363f6b7519516 Mon Sep 17 00:00:00 2001
From: ESCOBAR Juan <escj@nuwa>
Date: Tue, 15 Mar 2022 11:08:47 +0100
Subject: [PATCH] Juan 15/01/2022:ZSOLVER/sources_neg_correct.f90, add orig for
 modif for GPU

---
 src/ZSOLVER/sources_neg_correct.f90 | 483 ++++++++++++++++++++++++++++
 1 file changed, 483 insertions(+)
 create mode 100644 src/ZSOLVER/sources_neg_correct.f90

diff --git a/src/ZSOLVER/sources_neg_correct.f90 b/src/ZSOLVER/sources_neg_correct.f90
new file mode 100644
index 000000000..5eb676fc2
--- /dev/null
+++ b/src/ZSOLVER/sources_neg_correct.f90
@@ -0,0 +1,483 @@
+!MNH_LIC Copyright 2020-2022 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+!-----------------------------------------------------------------
+! Author: P. Wautelet 25/06/2020 (deduplication of code from advection_metsv, resolved_cloud and turb)
+! Modifications:
+!  P. Wautelet 30/06/2020: remove non-local corrections in resolved_cloud for NEGA => new local corrections here
+!  J. Escobar  21/07/2020: bug <-> array of size(:,:,:,0) => return if krr=0
+!  P. Wautelet 10/02/2021: budgets: add missing sources for NSV_C2R2BEG+3 budget
+!-----------------------------------------------------------------
+module mode_sources_neg_correct
+
+implicit none
+
+private
+
+public :: Sources_neg_correct
+
+contains
+
+subroutine Sources_neg_correct( hcloud, hbudname, krr, ptstep, ppabst, ptht, prt, prths, prrs, prsvs, prhodj )
+
+use modd_budget,     only: lbudget_th, lbudget_rv, lbudget_rc, lbudget_rr, lbudget_ri, &
+                           lbudget_rs, lbudget_rg, lbudget_rh, lbudget_sv,             &
+                           NBUDGET_TH, NBUDGET_RV, NBUDGET_RC, NBUDGET_RR, NBUDGET_RI, &
+                           NBUDGET_RS, NBUDGET_RG, NBUDGET_RH, NBUDGET_SV1,            &
+                           tbudgets
+use modd_cst,        only: xci, xcl, xcpd, xcpv, xlstt, xlvtt, xp00, xrd, xtt
+use modd_nsv,        only: nsv_c2r2beg, nsv_c2r2end, nsv_lima_beg, nsv_lima_end, nsv_lima_nc, nsv_lima_nr, nsv_lima_ni
+use modd_param_lima, only: lcold_lima => lcold, lrain_lima => lrain, lspro_lima => lspro, lwarm_lima => lwarm, &
+                           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_mem_get, Mnh_mem_position_pin, Mnh_mem_release
+#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
+character(len=*),            intent(in)           :: hbudname ! Budget name
+integer,                     intent(in)           :: krr      ! Number of moist variables
+real,                        intent(in)           :: ptstep   ! Timestep
+real, dimension(:, :, :),    intent(in)           :: ppabst   ! Absolute pressure at time t
+real, dimension(:, :, :),    intent(in)           :: ptht     ! Theta at time t
+real, dimension(:, :, :, :), intent(in)           :: prt      ! Moist variables at time t
+real, dimension(:, :, :),    intent(inout)        :: prths    ! Source terms
+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(:, :, :), pointer, contiguous :: zt, zexn, zlv, zls, zcph, zcor
+
+if ( krr == 0 ) return
+
+zcor => null()
+
+!$acc data present( ppabst, ptht, prt, prths, prrs, prsvs, prhodj )
+
+if ( mppdb_initialized ) then
+  !Check all IN arrays
+  call Mppdb_check( ppabst, "Sources_neg_correct beg:ppabst")
+  call Mppdb_check( ptht,   "Sources_neg_correct beg:ptht")
+  call Mppdb_check( prt,    "Sources_neg_correct beg:prt")
+  if ( Present( prhodj ) ) call Mppdb_check( prhodj, "Sources_neg_correct beg:prhodj")
+  !Check all INOUT arrays
+  call Mppdb_check( prths, "Sources_neg_correct beg:prths")
+  call Mppdb_check( prrs,  "Sources_neg_correct beg:prrs")
+  call Mppdb_check( prsvs, "Sources_neg_correct beg:prsvs")
+end if
+
+if ( hbudname /= 'NEADV' .and. hbudname /= 'NECON' .and. hbudname /= 'NEGA' .and. hbudname /= 'NETUR' ) &
+  call Print_msg( NVERB_WARNING, 'GEN', 'Sources_neg_correct', 'budget '//hbudname//' not yet tested' )
+
+if ( hcloud == 'LIMA' ) then
+  ! The negativity correction does not apply to the SPRO (supersaturation) variable which may be naturally negative
+  if ( lspro_lima ) then
+    isv_lima_end = nsv_lima_end - 1
+  else
+    isv_lima_end = nsv_lima_end
+  end if
+end if
+
+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
+    if ( lbudget_th ) call Budget_store_init( tbudgets(NBUDGET_TH), Trim( hbudname ), prths(:, :, :) )
+    if ( lbudget_rv ) call Budget_store_init( tbudgets(NBUDGET_RV), Trim( hbudname ), prrs (:, :, :, 1) )
+    if ( lbudget_rc ) call Budget_store_init( tbudgets(NBUDGET_RC), Trim( hbudname ), prrs (:, :, :, 2) )
+    if ( lbudget_rr .and.                                                                                   &
+       (   hbudname /= 'NETUR' .or.                                                                         &
+         ( hbudname == 'NETUR' .and. ( hcloud == 'C2R2' .or. hcloud == 'KHKO' .or. hcloud == 'LIMA' ) ) ) ) &
+                    call Budget_store_init( tbudgets(NBUDGET_RR), Trim( hbudname ), prrs (:, :, :, 3) )
+        IF (lbudget_ri .and.                                                                                    &
+       (   hbudname /= 'NETUR' .or.                                                                         &
+         ( hbudname == 'NETUR' .and. ( hcloud == 'ICE3' .or. hcloud == 'ICE4' .or. hcloud == 'LIMA' ) ) ) ) &
+                    call Budget_store_init( tbudgets(NBUDGET_RI), Trim( hbudname ), prrs (:, :, :, 4) )
+    if ( lbudget_rs .and. hbudname /= 'NETUR' ) call Budget_store_init( tbudgets(NBUDGET_RS), Trim( hbudname ), prrs (:, :, :, 5) )
+    if ( lbudget_rg .and. hbudname /= 'NETUR' ) call Budget_store_init( tbudgets(NBUDGET_RG), Trim( hbudname ), prrs (:, :, :, 6) )
+    if ( lbudget_rh .and. hbudname /= 'NETUR' ) call Budget_store_init( tbudgets(NBUDGET_RH), Trim( hbudname ), prrs (:, :, :, 7) )
+  end if
+
+  if ( lbudget_sv .and. ( hcloud == 'C2R2' .or. hcloud == 'KHKO' ) ) then
+    do ji = nsv_c2r2beg, nsv_c2r2end
+      call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + ji), Trim( hbudname ), prsvs(:, :, :, ji) )
+    end do
+  end if
+  if ( lbudget_sv .and. hcloud == 'LIMA' ) then
+    do ji = nsv_lima_beg, isv_lima_end
+      call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + ji), Trim( hbudname ), prsvs(:, :, :, ji) )
+    end do
+  end if
+else !NECON + NEGA
+  if ( .not. present( prhodj ) ) &
+    call Print_msg( NVERB_FATAL, 'GEN', 'Sources_neg_correct', 'optional argument prhodj not present' )
+
+  if ( hcloud == 'KESS' .or. hcloud == 'ICE3' .or. hcloud == 'ICE4' .or. &
+       hcloud == 'KHKO' .or. hcloud == 'C2R2' .or. hcloud == 'LIMA' ) then
+    if ( lbudget_th) call Budget_store_init( tbudgets(NBUDGET_TH), Trim( hbudname ), prths(:, :, :)    * prhodj(:, :, :) )
+    if ( lbudget_rv) call Budget_store_init( tbudgets(NBUDGET_RV), Trim( hbudname ), prrs (:, :, :, 1) * prhodj(:, :, :) )
+    if ( lbudget_rc) call Budget_store_init( tbudgets(NBUDGET_RC), Trim( hbudname ), prrs (:, :, :, 2) * prhodj(:, :, :) )
+    if ( lbudget_rr) call Budget_store_init( tbudgets(NBUDGET_RR), Trim( hbudname ), prrs (:, :, :, 3) * prhodj(:, :, :) )
+    if ( lbudget_ri) call Budget_store_init( tbudgets(NBUDGET_RI), Trim( hbudname ), prrs (:, :, :, 4) * prhodj(:, :, :) )
+    if ( lbudget_rs) call Budget_store_init( tbudgets(NBUDGET_RS), Trim( hbudname ), prrs (:, :, :, 5) * prhodj(:, :, :) )
+    if ( lbudget_rg) call Budget_store_init( tbudgets(NBUDGET_RG), Trim( hbudname ), prrs (:, :, :, 6) * prhodj(:, :, :) )
+    if ( lbudget_rh) call Budget_store_init( tbudgets(NBUDGET_RH), Trim( hbudname ), prrs (:, :, :, 7) * prhodj(:, :, :) )
+  end if
+
+  if ( lbudget_sv .and. ( hcloud == 'C2R2' .or. hcloud == 'KHKO' ) ) then
+    do ji = nsv_c2r2beg, nsv_c2r2end
+      call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + ji), Trim( hbudname ), prsvs(:, :, :, ji) * prhodj(:, :, :) )
+    end do
+  end if
+  if ( lbudget_sv .and. hcloud == 'LIMA' ) then
+    do ji = nsv_lima_beg, isv_lima_end
+      call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + ji), Trim( hbudname ), prsvs(:, :, :, ji) * prhodj(:, :, :) )
+    end do
+  end if
+end if
+
+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( jiu, jju, jku ) )
+  if ( krr > 3 ) then
+    allocate( zcor( jiu, jju, jku ) )
+  end if
+else
+  allocate( zls(0, 0, 0) )
+end if
+if ( .not. Associated( zcor ) ) Allocate( zcor(0, 0, 0) )
+#else
+!Pin positions in the pools of MNH memory
+call Mnh_mem_position_pin()
+
+call Mnh_mem_get( zt,   jiu, jju, jku )
+call Mnh_mem_get( zexn, jiu, jju, jku )
+call Mnh_mem_get( zlv,  jiu, jju, jku )
+call Mnh_mem_get( zcph, jiu, jju, jku )
+if ( hcloud == 'ICE3' .or. hcloud == 'ICE4' .or. hcloud == 'LIMA' ) then
+  call Mnh_mem_get( zls, jiu, jju, jku )
+  if ( krr > 3 ) then
+    call Mnh_mem_get( zcor, jiu, jju, jku )
+  else
+    call Mnh_mem_get( zcor, 0, 0, 0 )
+  end if
+else
+  call Mnh_mem_get( zls,  0, 0, 0 )
+  call Mnh_mem_get( zcor, 0, 0, 0 )
+end if
+#endif
+
+!$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
+if ( hcloud == 'ICE3' .or. hcloud == 'ICE4' .or. hcloud == 'LIMA' ) then
+!$acc kernels
+  zls(:, :, :) = xlstt + ( xcpv - xci ) * ( zt(:, :, :) - xtt )
+!$acc end kernels
+end if
+!$acc kernels
+zcph(:, :, :) = xcpd + xcpv * prt(:, :, :, 1)
+!$acc end kernels
+
+#ifndef MNH_OPENACC
+deallocate( zt )
+#endif
+
+CLOUD: select case ( hcloud )
+  case ( 'KESS' )
+    jrmax = Size( prrs, 4 )
+    do jr = 2, jrmax
+!PW: kernels directive inside do loop on jr because compiler bug... (NVHPC 21.7)
+!$acc kernels
+      where ( prrs(:, :, :, jr) < 0. )
+        prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, jr)
+        prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, jr) * zlv(:, :, :) /  &
+           ( zcph(:, :, :) * zexn(:, :, :) )
+        prrs(:, :, :, jr) = 0.
+      end where
+!$acc end kernels
+    end do
+
+!$acc kernels
+    where ( prrs(:, :, :, 1) < 0. .and. prrs(:, :, :, 2) > 0. )
+      prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, 2)
+      prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, 2) * zlv(:, :, :) /  &
+           ( zcph(:, :, :) * zexn(:, :, :) )
+      prrs(:, :, :, 2) = 0.
+    end where
+!$acc end kernels
+
+
+  case( 'ICE3', 'ICE4' )
+    if ( hbudname == 'NETUR' ) then
+      jrmax = 4
+    else
+      jrmax = Size( prrs, 4 )
+    end if
+!$acc kernels
+    do jr = 4, jrmax
+      do concurrent( ji = 1 : jiu, jj = 1 : jju, jk = 1 : jku )
+        if ( prrs(ji, jj, jk, jr) < 0. ) then
+          prrs(ji, jj, jk, 1) = prrs(ji, jj, jk, 1) + prrs(ji, jj, jk, jr)
+          prths(ji, jj, jk)   = prths(ji, jj, jk) - prrs(ji, jj, jk, jr) * zls(ji, jj, jk) /  &
+                                ( zcph(ji, jj, jk) * zexn(ji, jj, jk) )
+          prrs(ji, jj, jk, jr) = 0.
+        end if
+      end do
+    end do
+!$acc end kernels
+!
+!   cloud
+    if ( hbudname == 'NETUR' ) then
+      jrmax = 2
+    else
+      jrmax = 3
+    end if
+!$acc kernels
+    do jr = 2, jrmax
+      do concurrent( ji = 1 : jiu, jj = 1 : jju, jk = 1 : jku )
+        if ( prrs(ji, jj, jk, jr) < 0. ) then
+          prrs(ji, jj, jk, 1) = prrs(ji, jj, jk, 1) + prrs(ji, jj, jk, jr)
+          prths(ji, jj, jk)   = prths(ji, jj, jk) - prrs(ji, jj, jk, jr) * zlv(ji, jj, jk) /  &
+                                ( zcph(ji, jj, jk) * zexn(ji, jj, jk) )
+          prrs(ji, jj, jk, jr) = 0.
+        end if
+      end do
+    end do
+!
+! if rc or ri are positive, we can correct negative rv
+!   cloud
+    do concurrent( ji = 1 : jiu, jj = 1 : jju, jk = 1 : jku )
+      if ( prrs(ji, jj, jk, 1) < 0. .and. prrs(ji, jj, jk, 2) > 0. ) then
+        prrs(ji, jj, jk, 1) = prrs(ji, jj, jk, 1) + prrs(ji, jj, jk, 2)
+        prths(ji, jj, jk)   = prths(ji, jj, jk) - prrs(ji, jj, jk, 2) * zlv(ji, jj, jk) /  &
+                              ( zcph(ji, jj, jk) * zexn(ji, jj, jk) )
+        prrs(ji, jj, jk, 2) = 0.
+      end if
+    end do
+!   ice
+    if ( krr > 3 ) then
+#ifdef MNH_COMPILER_NVHPC
+!$acc loop independent collapse(3)
+#endif
+      do concurrent( ji = 1 : jiu, jj = 1 : jju, jk = 1 : jku )
+        if ( prrs(ji, jj, jk, 1) < 0. .and. prrs(ji, jj, jk, 4) > 0. ) then
+          zcor(ji, jj, jk)    = Min( -prrs(ji, jj, jk, 1), prrs(ji, jj, jk, 4) )
+          prrs(ji, jj, jk, 1) = prrs(ji, jj, jk, 1) + zcor(ji, jj, jk)
+          prths(ji, jj, jk)   = prths(ji, jj, jk) - zcor(ji, jj, jk) * zls(ji, jj, jk) / &
+                                ( zcph(ji, jj, jk) * zexn(ji, jj, jk) )
+          prrs(ji, jj, jk, 4) = prrs(ji, jj, jk, 4) - zcor(ji, jj, jk)
+        end if
+      end do
+    end if
+!$acc end kernels
+!
+!
+  case( 'C2R2', 'KHKO' )
+!$acc kernels
+    where ( prrs(:, :, :, 2) < 0. .or. prsvs(:, :, :, nsv_c2r2beg + 1) < 0. )
+      prsvs(:, :, :, nsv_c2r2beg) = 0.
+    end where
+!$acc end kernels
+    do jsv = 2, 3
+!$acc kernels
+!PW: kernels directive inside do loop on jr because compiler bug... (NVHPC 21.7)
+      where ( prrs(:, :, :, jsv) < 0. .or. prsvs(:, :, :, nsv_c2r2beg - 1 + jsv) < 0. )
+        prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, jsv)
+        prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, jsv) * zlv(:, :, :) /  &
+                ( zcph(:, :, :) * zexn(:, :, :) )
+        prrs(:, :, :, jsv)  = 0.
+        prsvs(:, :, :, nsv_c2r2beg - 1 + jsv) = 0.
+      end where
+!$acc end kernels
+    end do
+!$acc kernels
+    where ( prrs(:, :, :, 1) < 0. .and. prrs(:, :, :, 2) > 0. )
+      prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, 2)
+      prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, 2) * zlv(:, :, :) /  &
+           ( zcph(:, :, :) * zexn(:, :, :) )
+      prrs(:, :, :, 2) = 0.
+      prsvs(:, :, :, nsv_c2r2beg + 1) = 0.
+    end where
+!$acc end kernels
+!
+!
+  case( 'LIMA' )
+!$acc kernels
+! Correction where rc<0 or Nc<0
+    if ( lwarm_lima ) then
+      where ( prrs(:, :, :, 2) < xrtmin_lima(2) / ptstep .or. prsvs(:, :, :, nsv_lima_nc) < xctmin_lima(2) / ptstep )
+        prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, 2)
+        prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, 2) * zlv(:, :, :) /  &
+                 ( zcph(:, :, :) * zexn(:, :, :) )
+        prrs(:, :, :, 2)  = 0.
+        prsvs(:, :, :, nsv_lima_nc) = 0.
+      end where
+      where ( prrs(:, :, :, 1) < 0. .and. prrs(:, :, :, 2) > 0. )
+        prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, 2)
+        prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, 2) * zlv(:, :, :) /  &
+           ( zcph(:, :, :) * zexn(:, :, :) )
+        prrs(:, :, :, 2) = 0.
+        prsvs(:, :, :, nsv_lima_nc) = 0.
+      end where
+    end if
+! Correction where rr<0 or Nr<0
+    if ( lwarm_lima .and. lrain_lima ) then
+      where ( prrs(:, :, :, 3) < xrtmin_lima(3) / ptstep .or. prsvs(:, :, :, nsv_lima_nr) < xctmin_lima(3) / ptstep )
+        prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, 3)
+        prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, 3) * zlv(:, :, :) /  &
+                 ( zcph(:, :, :) * zexn(:, :, :) )
+        prrs(:, :, :, 3)  = 0.
+        prsvs(:, :, :, nsv_lima_nr) = 0.
+      end where
+    end if
+!$acc end kernels
+! Correction where ri<0 or Ni<0
+    if ( lcold_lima ) then
+!$acc kernels
+      where ( prrs(:, :, :, 4) < xrtmin_lima(4) / ptstep .or. prsvs(:, :, :, nsv_lima_ni) < xctmin_lima(4) / ptstep )
+        prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, 4)
+        prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, 4) * zls(:, :, :) /  &
+                 ( zcph(:, :, :) * zexn(:, :, :) )
+        prrs(:, :, :, 4)  = 0.
+        prsvs(:, :, :, nsv_lima_ni) = 0.
+      end where
+!$acc end kernels
+      if ( hbudname /= 'NETUR' ) then
+        do jr = 5, Size( prrs, 4 )
+!PW: kernels directive inside do loop on jr because compiler bug... (NVHPC 21.7)
+!$acc kernels
+          where ( prrs(:, :, :, jr) < 0. )
+            prrs(:, :, :, 1) = prrs(:, :, :, 1) + prrs(:, :, :, jr)
+            prths(:, :, :) = prths(:, :, :) - prrs(:, :, :, jr) * zls(:, :, :) /  &
+                    ( zcph(:, :, :) * zexn(:, :, :) )
+            prrs(:, :, :, jr) = 0.
+          end where
+!$acc end kernels
+        end do
+      end if
+      if(krr > 3) then
+!$acc kernels
+        where ( prrs(:, :, :, 1) < 0. .and. prrs(:, :, :, 4) > 0. )
+          zcor(:, :, :) = Min( -prrs(:, :, :, 1), prrs(:, :, :, 4) )
+          prrs(:, :, :, 1) = prrs(:, :, :, 1) + zcor(:, :, :)
+          prths(:, :, :) = prths(:, :, :) - zcor(:, :, :) * zls(:, :, :) /  &
+             ( zcph(:, :, :) * zexn(:, :, :) )
+          prrs(:, :, :, 4) = prrs(:, :, :, 4) - zcor(:, :, :)
+        end where
+!$acc end kernels
+      end if
+    end if
+
+!$acc kernels
+    prsvs(:, :, :, nsv_lima_beg : isv_lima_end) = Max( 0.0, prsvs(:, :, :, nsv_lima_beg : isv_lima_end) )
+!$acc end kernels
+
+end select CLOUD
+
+!$acc end data
+
+#ifndef MNH_OPENACC
+deallocate( zexn, zlv, zcph, zls, zcor )
+#else
+!Release all memory allocated with Mnh_mem_get calls since last call to Mnh_mem_position_pin
+call Mnh_mem_release()
+#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
+    if ( lbudget_th ) call Budget_store_end( tbudgets(NBUDGET_TH), Trim( hbudname ), prths(:, :, :) )
+    if ( lbudget_rv ) call Budget_store_end( tbudgets(NBUDGET_RV), Trim( hbudname ), prrs (:, :, :, 1) )
+    if ( lbudget_rc ) call Budget_store_end( tbudgets(NBUDGET_RC), Trim( hbudname ), prrs (:, :, :, 2) )
+    if ( lbudget_rr .and.                                                                                   &
+       (   hbudname /= 'NETUR' .or.                                                                         &
+         ( hbudname == 'NETUR' .and. ( hcloud == 'C2R2' .or. hcloud == 'KHKO' .or. hcloud == 'LIMA' ) ) ) ) &
+                    call Budget_store_end( tbudgets(NBUDGET_RR), Trim( hbudname ), prrs (:, :, :, 3) )
+        IF (lbudget_ri .and.                                                                                    &
+       (   hbudname /= 'NETUR' .or.                                                                         &
+         ( hbudname == 'NETUR' .and. ( hcloud == 'ICE3' .or. hcloud == 'ICE4' .or. hcloud == 'LIMA' ) ) ) ) &
+                    call Budget_store_end( tbudgets(NBUDGET_RI), Trim( hbudname ), prrs (:, :, :, 4) )
+    if ( lbudget_rs .and. hbudname /= 'NETUR' ) call Budget_store_end( tbudgets(NBUDGET_RS), Trim( hbudname ), prrs (:, :, :, 5) )
+    if ( lbudget_rg .and. hbudname /= 'NETUR' ) call Budget_store_end( tbudgets(NBUDGET_RG), Trim( hbudname ), prrs (:, :, :, 6) )
+    if ( lbudget_rh .and. hbudname /= 'NETUR' ) call Budget_store_end( tbudgets(NBUDGET_RH), Trim( hbudname ), prrs (:, :, :, 7) )
+  end if
+
+  if ( lbudget_sv .and. ( hcloud == 'C2R2' .or. hcloud == 'KHKO' ) ) then
+    do ji = nsv_c2r2beg, nsv_c2r2end
+      call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + ji), Trim( hbudname ), prsvs(:, :, :, ji) )
+    end do
+  end if
+  if ( lbudget_sv .and. hcloud == 'LIMA' ) then
+    do ji = nsv_lima_beg, isv_lima_end
+      call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + ji), Trim( hbudname ), prsvs(:, :, :, ji) )
+    end do
+  end if
+else !NECON + NEGA
+  if ( hcloud == 'KESS' .or. hcloud == 'ICE3' .or. hcloud == 'ICE4' .or. &
+       hcloud == 'KHKO' .or. hcloud == 'C2R2' .or. hcloud == 'LIMA' ) then
+    if ( lbudget_th) call Budget_store_end( tbudgets(NBUDGET_TH), Trim( hbudname ), prths(:, :, :)    * prhodj(:, :, :) )
+    if ( lbudget_rv) call Budget_store_end( tbudgets(NBUDGET_RV), Trim( hbudname ), prrs (:, :, :, 1) * prhodj(:, :, :) )
+    if ( lbudget_rc) call Budget_store_end( tbudgets(NBUDGET_RC), Trim( hbudname ), prrs (:, :, :, 2) * prhodj(:, :, :) )
+    if ( lbudget_rr) call Budget_store_end( tbudgets(NBUDGET_RR), Trim( hbudname ), prrs (:, :, :, 3) * prhodj(:, :, :) )
+    if ( lbudget_ri) call Budget_store_end( tbudgets(NBUDGET_RI), Trim( hbudname ), prrs (:, :, :, 4) * prhodj(:, :, :) )
+    if ( lbudget_rs) call Budget_store_end( tbudgets(NBUDGET_RS), Trim( hbudname ), prrs (:, :, :, 5) * prhodj(:, :, :) )
+    if ( lbudget_rg) call Budget_store_end( tbudgets(NBUDGET_RG), Trim( hbudname ), prrs (:, :, :, 6) * prhodj(:, :, :) )
+    if ( lbudget_rh) call Budget_store_end( tbudgets(NBUDGET_RH), Trim( hbudname ), prrs (:, :, :, 7) * prhodj(:, :, :) )
+  end if
+
+  if ( lbudget_sv .and. ( hcloud == 'C2R2' .or. hcloud == 'KHKO' ) ) then
+    do ji = nsv_c2r2beg, nsv_c2r2end
+      call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + ji), Trim( hbudname ), prsvs(:, :, :, ji) * prhodj(:, :, :) )
+    end do
+  end if
+  if ( lbudget_sv .and. hcloud == 'LIMA' ) then
+    do ji = nsv_lima_beg, isv_lima_end
+      call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + ji), Trim( hbudname ), prsvs(:, :, :, ji) * prhodj(:, :, :) )
+    end do
+  end if
+end if
+
+if ( mppdb_initialized ) then
+  !Check all INOUT arrays
+  call Mppdb_check( prths, "Sources_neg_correct end:prths")
+  call Mppdb_check( prrs,  "Sources_neg_correct end:prrs")
+  call Mppdb_check( prsvs, "Sources_neg_correct end:prsvs")
+end if
+
+!$acc end data
+
+end subroutine Sources_neg_correct
+
+end module mode_sources_neg_correct
-- 
GitLab