From bd4dfbfe31fc69682c66f6f982e8429369d72f8f Mon Sep 17 00:00:00 2001
From: Philippe WAUTELET <philippe.wautelet@aero.obs-mip.fr>
Date: Wed, 17 Nov 2021 12:56:43 +0100
Subject: [PATCH] Philippe 17/11/2021: ZSOLVER: port of last development to
5.5.0 version
---
src/ZSOLVER/advection_metsv.f90 | 245 +++++------
src/ZSOLVER/advection_uvw_cen.f90 | 51 +--
src/ZSOLVER/gdiv.f90 | 1 -
src/ZSOLVER/pressure.f90 | 689 ------------------------------
src/ZSOLVER/viscosity.f90 | 83 ++--
5 files changed, 202 insertions(+), 867 deletions(-)
delete mode 100644 src/ZSOLVER/pressure.f90
diff --git a/src/ZSOLVER/advection_metsv.f90 b/src/ZSOLVER/advection_metsv.f90
index f148ad521..481586dd4 100644
--- a/src/ZSOLVER/advection_metsv.f90
+++ b/src/ZSOLVER/advection_metsv.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1994-2019 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1994-2021 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.
@@ -8,7 +8,7 @@
! ###########################
!
INTERFACE
- SUBROUTINE ADVECTION_METSV (TPFILE, OCLOSE_OUT,HUVW_ADV_SCHEME, &
+ SUBROUTINE ADVECTION_METSV (TPFILE, HUVW_ADV_SCHEME, &
HMET_ADV_SCHEME,HSV_ADV_SCHEME, HCLOUD, KSPLIT, &
OSPLIT_CFL, PSPLIT_CFL, OCFL_WRIT, &
HLBCX, HLBCY, KRR, KSV, TPDTCUR, PTSTEP, &
@@ -17,11 +17,9 @@ INTERFACE
PRTHS, PRRS, PRTKES, PRSVS, &
PRTHS_CLD, PRRS_CLD, PRSVS_CLD, PRTKES_ADV )
!
-USE MODD_IO, ONLY: TFILEDATA
+USE MODD_IO, ONLY: TFILEDATA
USE MODD_TYPE_DATE, ONLY: DATE_TIME
!
-LOGICAL, INTENT(IN) :: OCLOSE_OUT ! switch for syncronous
- ! file opening
TYPE(TFILEDATA), INTENT(IN) :: TPFILE ! Output file
CHARACTER(LEN=6), INTENT(IN) :: HMET_ADV_SCHEME, & ! Control of the
HSV_ADV_SCHEME, & ! scheme applied
@@ -65,7 +63,7 @@ END INTERFACE
!
END MODULE MODI_ADVECTION_METSV
! ##########################################################################
- SUBROUTINE ADVECTION_METSV (TPFILE, OCLOSE_OUT,HUVW_ADV_SCHEME, &
+ SUBROUTINE ADVECTION_METSV (TPFILE, HUVW_ADV_SCHEME, &
HMET_ADV_SCHEME,HSV_ADV_SCHEME, HCLOUD, KSPLIT, &
OSPLIT_CFL, PSPLIT_CFL, OCFL_WRIT, &
HLBCX, HLBCY, KRR, KSV, TPDTCUR, PTSTEP, &
@@ -134,10 +132,14 @@ END MODULE MODI_ADVECTION_METSV
!! 04/2016 (C.Lac) : correction of negativity for KHKO
!! 10/2016 (C.Lac) Correction on the flag for Strang splitting
!! to insure reproducibility between START and RESTA
-!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
-!! 07/2017 (V. Vionnet) : add advection of 2D variables at
-!! the surface for the blowing snow scheme
+! V. Vionnet 07/2017: add advection of 2D variables at the surface for the blowing snow scheme
+! P. Wautelet 05/2016-04/2018: new data structures and calls for I/O
+! P. Wautelet 02/2020: use the new data structures and subroutines for budgets
+! B. Vie 03/2020: LIMA negativity checks after turbulence, advection and microphysics budgets
! P. Wautelet 20/05/2019: add name argument to ADDnFIELD_ll + new ADD4DFIELD_ll subroutine
+! P. Wautelet 11/06/2020: bugfix: correct PRSVS array indices
+! P. Wautelet + Benoît Vié 06/2020: improve removal of negative scalar variables + adapt the corresponding budgets
+! P. Wautelet 30/06/2020: move removal of negative scalar variables to Sources_neg_correct
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -148,55 +150,60 @@ USE MODD_MPIF
use modd_precision, only: MNHREAL_MPI
USE MODD_VAR_ll, ONLY: NMNH_COMM_WORLD
#endif
-USE MODD_BUDGET
-USE MODD_CST
-USE MODD_CTURB, ONLY: XTKEMIN
+USE MODD_BLOWSNOW
+USE MODD_BLOWSNOW_n
+use modd_budget, only: lbudget_th, lbudget_tke, lbudget_rv, lbudget_rc, &
+ lbudget_rr, lbudget_ri, lbudget_rs, lbudget_rg, lbudget_rh, lbudget_sv, &
+ NBUDGET_TH, NBUDGET_TKE, NBUDGET_RV, NBUDGET_RC, &
+ NBUDGET_RR, NBUDGET_RI, NBUDGET_RS, NBUDGET_RG, NBUDGET_RH, NBUDGET_SV1, &
+ tbudgets
USE MODD_CONF, ONLY: LNEUTRAL,NHALO,L1D, L2D
+USE MODD_CST
+USE MODD_CTURB, ONLY: XTKEMIN
+use modd_field, only: tfielddata, TYPEREAL
+USE MODD_IBM_PARAM_n, ONLY: LIBM,XIBM_LS,XIBM_EPSI
USE MODD_IO, ONLY: TFILEDATA
USE MODD_LUNIT_n, ONLY: TLUOUT
+USE MODD_NSV
+USE MODD_PARAMETERS
+USE MODD_PARAM_LIMA
USE MODD_PARAM_n
+USE MODD_REF_n, ONLY: XRHODJ,XRHODREF
USE MODD_TYPE_DATE, ONLY: DATE_TIME
-USE MODD_BLOWSNOW
-USE MODD_BLOWSNOW_n
-USE MODD_PARAMETERS
!
use mode_argslist_ll, only: ADD3DFIELD_ll, ADD4DFIELD_ll, CLEANLIST_ll, LIST_ll
+use mode_budget, only: Budget_store_init, Budget_store_end
+#ifdef MNH_OPENACC
+USE MODE_DEVICE
+USE MODE_MNH_ZWORK, ONLY: MNH_ALLOCATE_ZT3D, MNH_ALLOCATE_ZT4D, MNH_CHECK_IN_ZT3D, MNH_CHECK_OUT_ZT3D, &
+ MNH_GET_ZT3D, MNH_REL_ZT3D, MNH_REL_ZT4D, ZT3D
+#endif
use mode_exchange_ll, only: UPDATE_HALO_ll
-USE MODE_FIELD, ONLY: TFIELDDATA, TYPEREAL
USE MODE_IO_FIELD_WRITE, only: IO_Field_write
use mode_mppdb
USE MODE_MSG
+use mode_sources_neg_correct, only: Sources_neg_correct
+#ifndef MNH_OPENACC
use mode_sum_ll, only: MAX_ll
+#endif
use mode_tools_ll, only: GET_INDICE_ll, lnorth_ll, lsouth_ll, least_ll, lwest_ll
!
USE MODI_ADV_BOUNDARIES
-USE MODI_BUDGET
+#ifdef MNH_BITREP
+USE MODI_BITREP
+#endif
USE MODI_CONTRAV
USE MODI_GET_HALO
USE MODI_PPM_RHODJ
USE MODI_PPM_MET
USE MODI_PPM_SCALAR
!
-#ifdef MNH_OPENACC
-USE MODE_DEVICE
-USE MODE_MNH_ZWORK, ONLY : ZT3D, MNH_GET_ZT3D , MNH_REL_ZT3D
-USE MODE_MNH_ZWORK, ONLY : MNH_ALLOCATE_ZT3D , MNH_REL_ZT3D, MNH_ALLOCATE_ZT3DP , MNH_ALLOCATE_ZT2D, &
- MNH_ALLOCATE_ZT4D , MNH_REL_ZT4D, &
- MNH_CHECK_IN_ZT3D,MNH_CHECK_OUT_ZT3D
-USE MODI_GET_HALO
-#endif
-#ifdef MNH_BITREP
-USE MODI_BITREP
-#endif
-!
!-------------------------------------------------------------------------------
!
IMPLICIT NONE
!
!* 0.1 Declarations of dummy arguments :
!
-LOGICAL, INTENT(IN) :: OCLOSE_OUT ! switch for synchronous
- ! file opening
TYPE(TFILEDATA), INTENT(IN) :: TPFILE ! Output file
CHARACTER(LEN=6), INTENT(IN) :: HMET_ADV_SCHEME, & ! Control of the
HSV_ADV_SCHEME, & ! scheme applied
@@ -250,7 +257,9 @@ REAL, DIMENSION(:,:,:),pointer , contiguous :: ZCFLW
! ! CFL numbers on each direction
REAL, DIMENSION(:,:,:),pointer , contiguous :: ZCFL
! ! CFL number
+#ifdef MNH_OPENACC
INTEGER :: IZRUCPPM,IZRVCPPM,IZRWCPPM,IZCFLU,IZCFLV,IZCFLW,IZCFL
+#endif
!
REAL :: ZCFLU_MAX, ZCFLV_MAX, ZCFLW_MAX, ZCFL_MAX ! maximum CFL numbers
!
@@ -260,30 +269,40 @@ REAL, DIMENSION(:,:,:),pointer , contiguous :: ZRTHS_OTHER
REAL, DIMENSION(:,:,:),pointer , contiguous :: ZRTKES_OTHER
REAL, DIMENSION(:,:,:),pointer , contiguous :: ZRTHS_PPM
REAL, DIMENSION(:,:,:),pointer , contiguous :: ZRTKES_PPM
+#ifdef MNH_OPENACC
INTEGER :: IZTH,IZTKE,IZRTHS_OTHER,IZRTKES_OTHER,IZRTHS_PPM,IZRTKES_PPM
+#endif
REAL, DIMENSION(:,:,:,:),pointer , contiguous :: ZR
REAL, DIMENSION(:,:,:,:),pointer , contiguous :: ZSV
REAL, DIMENSION(:,:,:,:),pointer , contiguous :: ZSNWC
REAL, DIMENSION(:,:,:,:),pointer , contiguous :: ZSNWC_INIT
REAL, DIMENSION(:,:,:,:),pointer , contiguous :: ZRSNWCS
+#ifdef MNH_OPENACC
INTEGER :: IZR,IZSV,IZSNWC,IZSNWC_INIT,IZRSNWCS
+#endif
! Guess at the sub time step
REAL, DIMENSION(:,:,:,:),pointer , contiguous :: ZRRS_OTHER
REAL, DIMENSION(:,:,:,:),pointer , contiguous :: ZRSVS_OTHER
REAL, DIMENSION(:,:,:,:),pointer , contiguous :: ZRSNWCS_OTHER
+#ifdef MNH_OPENACC
INTEGER :: IZRRS_OTHER,IZRSVS_OTHER,IZRSNWCS_OTHER
+#endif
! Tendencies since the beginning of the time step
REAL, DIMENSION(:,:,:,:),pointer , contiguous :: ZRRS_PPM
REAL, DIMENSION(:,:,:,:),pointer , contiguous :: ZRSVS_PPM
REAL, DIMENSION(:,:,:,:),pointer , contiguous :: ZRSNWCS_PPM
+#ifdef MNH_OPENACC
INTEGER :: IZRRS_PPM,IZRSVS_PPM,IZRSNWCS_PPM
+#endif
! Guess at the end of the sub time step
REAL, DIMENSION(:,:,:),pointer , contiguous :: ZRHOX1,ZRHOX2
REAL, DIMENSION(:,:,:),pointer , contiguous :: ZRHOY1,ZRHOY2
REAL, DIMENSION(:,:,:),pointer , contiguous :: ZRHOZ1,ZRHOZ2
REAL, DIMENSION(:,:,:),pointer , contiguous :: ZT,ZEXN,ZLV,ZLS,ZCPH
+#ifdef MNH_OPENACC
INTEGER :: IZRHOX1,IZRHOX2,IZRHOY1,IZRHOY2,IZRHOZ1,IZRHOZ2 &
,IZT,IZEXN,IZLV,IZLS,IZCPH
+#endif
! Temporary advected rhodj for PPM routines
!
@@ -301,6 +320,7 @@ INTEGER :: ILUOUT ! logical unit
INTEGER :: ISPLIT_PPM ! temporal time splitting
INTEGER :: IIB, IIE, IJB, IJE,IKB,IKE
#ifdef MNH_OPENACC
+CHARACTER(LEN=3) :: YNUM
INTEGER :: IZ1, IZ2
#endif
TYPE(TFIELDDATA) :: TZFIELD
@@ -427,16 +447,30 @@ IZCPH = MNH_ALLOCATE_ZT3D( ZCPH , JIU,JJU,JKU )
!
!* 0. INITIALIZATION
! --------------
-!
+
+GTKE=(SIZE(PTKET)/=0)
+
+if ( lbudget_th ) call Budget_store_init( tbudgets(NBUDGET_TH ), 'ADV', prths (:, :, :) )
+if ( lbudget_tke ) call Budget_store_init( tbudgets(NBUDGET_TKE), 'ADV', prtkes(:, :, :) )
+if ( lbudget_rv ) call Budget_store_init( tbudgets(NBUDGET_RV ), 'ADV', prrs (:, :, :, 1) )
+if ( lbudget_rc ) call Budget_store_init( tbudgets(NBUDGET_RC ), 'ADV', prrs (:, :, :, 2) )
+if ( lbudget_rr ) call Budget_store_init( tbudgets(NBUDGET_RR ), 'ADV', prrs (:, :, :, 3) )
+if ( lbudget_ri ) call Budget_store_init( tbudgets(NBUDGET_RI ), 'ADV', prrs (:, :, :, 4) )
+if ( lbudget_rs ) call Budget_store_init( tbudgets(NBUDGET_RS ), 'ADV', prrs (:, :, :, 5) )
+if ( lbudget_rg ) call Budget_store_init( tbudgets(NBUDGET_RG ), 'ADV', prrs (:, :, :, 6) )
+if ( lbudget_rh ) call Budget_store_init( tbudgets(NBUDGET_RH ), 'ADV', prrs (:, :, :, 7) )
+if ( lbudget_sv) then
+ do jsv = 1, ksv
+ call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + jsv ), 'ADV', prsvs(:, :, :, jsv) )
+ end do
+end if
+
ILUOUT = TLUOUT%NLU
!
CALL GET_INDICE_ll(IIB,IJB,IIE,IJE)
IKB=1+JPVEXT
IKE=SIZE(PSVT,3) - JPVEXT
-!
-GTKE=(SIZE(PTKET)/=0)
-!
#ifdef MNH_OPENACC
CALL INIT_ON_HOST_AND_DEVICE(PRTKES_ADV,PVALUE=-1e99,HNAME='ADVECTION_METSV::PRTKES_ADV')
CALL INIT_ON_HOST_AND_DEVICE(ZRUCPPM,PVALUE=-1e90,HNAME='ADVECTION_METSV::ZRUCPPM')
@@ -516,6 +550,26 @@ IF (.NOT. L1D) THEN
ZCFLU(IIB:IIE,IJB:IJE,:) = ABS(ZRUCPPM(IIB:IIE,IJB:IJE,:) * PTSTEP)
ZCFLV(IIB:IIE,IJB:IJE,:) = ABS(ZRVCPPM(IIB:IIE,IJB:IJE,:) * PTSTEP)
ZCFLW(IIB:IIE,IJB:IJE,:) = ABS(ZRWCPPM(IIB:IIE,IJB:IJE,:) * PTSTEP)
+ IF (LIBM) THEN
+#ifndef MNH_BITREP
+ ZCFLU(IIB:IIE,IJB:IJE,:) = ZCFLU(IIB:IIE,IJB:IJE,:)*(1.-exp(-(XIBM_LS(IIB:IIE,IJB:IJE,:,2)/&
+ (XRHODJ(IIB:IIE,IJB:IJE,:)/XRHODREF(IIB:IIE,IJB:IJE,:))**(1./3.))**2.))
+ ZCFLV(IIB:IIE,IJB:IJE,:) = ZCFLV(IIB:IIE,IJB:IJE,:)*(1.-exp(-(XIBM_LS(IIB:IIE,IJB:IJE,:,3)/&
+ (XRHODJ(IIB:IIE,IJB:IJE,:)/XRHODREF(IIB:IIE,IJB:IJE,:))**(1./3.))**2.))
+ ZCFLW(IIB:IIE,IJB:IJE,:) = ZCFLW(IIB:IIE,IJB:IJE,:)*(1.-exp(-(XIBM_LS(IIB:IIE,IJB:IJE,:,4)/&
+ (XRHODJ(IIB:IIE,IJB:IJE,:)/XRHODREF(IIB:IIE,IJB:IJE,:))**(1./3.))**2.))
+#else
+ ZCFLU(IIB:IIE,IJB:IJE,:) = ZCFLU(IIB:IIE,IJB:IJE,:)*(1.-Br_exp(-Br_pow(XIBM_LS(IIB:IIE,IJB:IJE,:,2)/&
+ Br_pow(XRHODJ(IIB:IIE,IJB:IJE,:)/XRHODREF(IIB:IIE,IJB:IJE,:),1./3.),2.)))
+ ZCFLV(IIB:IIE,IJB:IJE,:) = ZCFLV(IIB:IIE,IJB:IJE,:)*(1.-Br_exp(-Br_pow(XIBM_LS(IIB:IIE,IJB:IJE,:,3)/&
+ Br_pow(XRHODJ(IIB:IIE,IJB:IJE,:)/XRHODREF(IIB:IIE,IJB:IJE,:),1./3.),2.)))
+ ZCFLW(IIB:IIE,IJB:IJE,:) = ZCFLW(IIB:IIE,IJB:IJE,:)*(1.-Br_exp(-Br_pow(XIBM_LS(IIB:IIE,IJB:IJE,:,4)/&
+ Br_pow(XRHODJ(IIB:IIE,IJB:IJE,:)/XRHODREF(IIB:IIE,IJB:IJE,:),1./3.),2.)))
+#endif
+ WHERE (XIBM_LS(IIB:IIE,IJB:IJE,:,2).GT.(-XIBM_EPSI)) ZCFLU(IIB:IIE,IJB:IJE,:)=0.
+ WHERE (XIBM_LS(IIB:IIE,IJB:IJE,:,3).GT.(-XIBM_EPSI)) ZCFLV(IIB:IIE,IJB:IJE,:)=0.
+ WHERE (XIBM_LS(IIB:IIE,IJB:IJE,:,4).GT.(-XIBM_EPSI)) ZCFLW(IIB:IIE,IJB:IJE,:)=0.
+ ENDIF
!$acc end kernels
#ifndef MNH_BITREP
IF (.NOT. L2D) THEN
@@ -562,7 +616,7 @@ END IF
!
!* prints in the file the 3D Courant numbers (one should flag this)
!
-IF (OCLOSE_OUT .AND. OCFL_WRIT .AND. (.NOT. L1D)) THEN
+IF ( tpfile%lopened .AND. OCFL_WRIT .AND. (.NOT. L1D) ) THEN
!$acc update host(ZCFLU,ZCFLV,ZCFLW,ZCFL)
TZFIELD%CMNHNAME = 'CFLU'
TZFIELD%CSTDNAME = ''
@@ -695,7 +749,7 @@ END IF
ZTSTEP_PPM = PTSTEP / REAL(KSPLIT)
!
!
-!* 2.4 normalized contravariant components for splitted PPM time-step
+!* 2.4 normalized contravariant components for split PPM time-step
!
!$acc kernels
ZRUCPPM(:,:,:) = ZRUCPPM(:,:,:)*ZTSTEP_PPM
@@ -773,14 +827,19 @@ NULLIFY(TZFIELDS0_ll)
CALL GET_HALO_D(ZRTHS_OTHER, HNAME='ADVECTION_METSV::ZRTHS_OTHER')
IF (GTKE) CALL GET_HALO_D(ZRTKES_OTHER,HNAME='ADVECTION_METSV::ZRTKES_OTHER')
DO JR=1,KRR
- CALL GET_HALO_D(ZRRS_OTHER(:,:,:,JR),HNAME='ADVECTION_METSV::ZRRS_OTHER')
+ WRITE(YNUM, '( I3.3 )' ) JR
+ CALL GET_HALO_D(ZRRS_OTHER(:,:,:,JR),HNAME='ADVECTION_METSV::ZRRS_OTHER:'//YNUM)
END DO
DO JSV = 1, KSV
- CALL GET_HALO_D(ZRSVS_OTHER(:,:,:,JSV),HNAME='ADVECTION_METSV::ZRSVS_OTHER')
- END DO
- DO JSV = 1,NBLOWSNOW_2D
- CALL GET_HALO_D(ZRSNWCS_OTHER(:,:,:,JSV),HNAME='ADVECTION_METSV::ZRSNWCS_OTHER')
+ WRITE(YNUM, '( I3.3 )' ) JSV
+ CALL GET_HALO_D(ZRSVS_OTHER(:,:,:,JSV),HNAME='ADVECTION_METSV::ZRSVS_OTHER:'//YNUM)
END DO
+ IF ( LBLOWSNOW ) THEN
+ DO JSV = 1,NBLOWSNOW_2D
+ WRITE(YNUM, '( I3.3 )' ) JSV
+ CALL GET_HALO_D(ZRSNWCS_OTHER(:,:,:,JSV),HNAME='ADVECTION_METSV::ZRSNWCS_OTHER:'//YNUM)
+ END DO
+ END IF
!PW: TODO: update only what is needed...
! acc update device(ZRTHS_OTHER,ZRTKES_OTHER,ZRRS_OTHER,ZRSVS_OTHER)
#endif
@@ -886,9 +945,6 @@ DO JSPL=1,KSPLIT
!$acc end kernels
END IF
! acc end kernels
-
-!PW: bug PGI 18.10: not necessary for PRRS,PRSVS but error with decriptor not present
-!$acc update self(PRRS,PRSVS)
!
! Top and bottom Boundaries and LBC for the guesses
!
@@ -958,9 +1014,11 @@ DO JSPL=1,KSPLIT
DO JSV = 1, KSV
CALL GET_HALO_D(ZSV(:,:,:,JSV),HNAME='ADVECTION_METSV::ZSV')
END DO
- DO JSV = 1,NBLOWSNOW_2D
- CALL GET_HALO_D(ZSNWC(:,:,:,JSV),HNAME='ADVECTION_METSV::ZSNWC')
- END DO
+ IF ( LBLOWSNOW ) THEN
+ DO JSV = 1,NBLOWSNOW_2D
+ CALL GET_HALO_D(ZSNWC(:,:,:,JSV),HNAME='ADVECTION_METSV::ZSNWC')
+ END DO
+ END IF
#endif
END IF
@@ -1005,84 +1063,27 @@ END IF
!* 5. BUDGETS
! -------
!
-IF (LBUDGET_TH) THEN
-!$acc update self(PRTHS)
- CALL BUDGET (PRTHS,4,'ADV_BU_RTH')
-END IF
-IF (LBUDGET_TKE) THEN
-!$acc update self(PRTKES)
- CALL BUDGET (PRTKES,5,'ADV_BU_RTKE')
-END IF
-IF (KRR>=1.AND.LBUDGET_RV) THEN
-!$acc update self(PRRS(:,:,:,1))
- CALL BUDGET (PRRS(:,:,:,1),6,'ADV_BU_RRV')
-END IF
-IF (KRR>=2.AND.LBUDGET_RC) THEN
-!$acc update self(PRRS(:,:,:,2))
- CALL BUDGET (PRRS(:,:,:,2),7,'ADV_BU_RRC')
-END IF
-IF (KRR>=3.AND.LBUDGET_RR) THEN
-!$acc update self(PRRS(:,:,:,3))
- CALL BUDGET (PRRS(:,:,:,3),8,'ADV_BU_RRR')
-END IF
-IF (KRR>=4.AND.LBUDGET_RI) THEN
-!$acc update self(PRRS(:,:,:,4))
- CALL BUDGET (PRRS(:,:,:,4),9,'ADV_BU_RRI')
-END IF
-IF (KRR>=5.AND.LBUDGET_RS) THEN
-!$acc update self(PRRS(:,:,:,5))
- CALL BUDGET (PRRS(:,:,:,5),10,'ADV_BU_RRS')
-END IF
-IF (KRR>=6.AND.LBUDGET_RG) THEN
-!$acc update self(PRRS(:,:,:,6))
- CALL BUDGET (PRRS(:,:,:,6),11,'ADV_BU_RRG')
-END IF
-IF (KRR>=7.AND.LBUDGET_RH) THEN
-!$acc update self(PRRS(:,:,:,7))
- CALL BUDGET (PRRS(:,:,:,7),12,'ADV_BU_RRH')
-END IF
-IF (LBUDGET_SV) THEN
-!$acc update self(PRSVS)
- DO JSV=1,KSV
- CALL BUDGET (PRSVS(:,:,:,JSV),JSV+12,'ADV_BU_RSV')
- END DO
-END IF
-!
-IF ((HCLOUD == 'KHKO') .OR. (HCLOUD == 'C2R2')) THEN
-#ifdef MNH_OPENACC
- call Print_msg( NVERB_ERROR, 'GEN', 'ADVECTION_METSV', 'OpenACC: HCLOUD=''KHKO'' OR ''C2R2'' not yet implemented' )
-#endif
- ZEXN(:,:,:)= (PPABST(:,:,:)/XP00)**(XRD/XCPD)
- ZT(:,:,:)= PTHT(:,:,:)*ZEXN(:,:,:)
- ZLV(:,:,:)=XLVTT +(XCPV-XCL) *(ZT(:,:,:)-XTT)
- ZLS(:,:,:)=XLSTT +(XCPV-XCI) *(ZT(:,:,:)-XTT)
- ZCPH(:,:,:)=XCPD +XCPV*PRT(:,:,:,1)
-! CALL GET_HALO(PRRS(:,:,:,2))
-! CALL GET_HALO(PRSVS(:,:,:,2))
-! CALL GET_HALO(PRSVS(:,:,:,3))
- WHERE (PRRS(:,:,:,2) < 0. .OR. PRSVS(:,:,:,2) < 0.)
- PRSVS(:,:,:,1) = 0.0
- END WHERE
- DO JSV = 2, 3
- WHERE (PRRS(:,:,:,JSV) < 0. .OR. PRSVS(:,:,:,JSV) < 0.)
- PRRS(:,:,:,1) = PRRS(:,:,:,1) + PRRS(:,:,:,JSV)
- PRTHS(:,:,:) = PRTHS(:,:,:) - PRRS(:,:,:,JSV) * ZLV(:,:,:) / &
- ZCPH(:,:,:) / ZEXN(:,:,:)
- PRRS(:,:,:,JSV) = 0.0
- PRSVS(:,:,:,JSV) = 0.0
- END WHERE
- END DO
-!
- IF (LBUDGET_TH) CALL BUDGET (PRTHS(:,:,:) , 4,'NEADV_BU_RTH')
- IF (LBUDGET_RV) CALL BUDGET (PRRS(:,:,:,1), 6,'NEADV_BU_RRV')
- IF (LBUDGET_RC) CALL BUDGET (PRRS(:,:,:,2), 7,'NEADV_BU_RRC')
+if ( lbudget_th ) call Budget_store_end( tbudgets(NBUDGET_TH ), 'ADV', prths (:, :, :) )
+if ( lbudget_tke ) call Budget_store_end( tbudgets(NBUDGET_TKE), 'ADV', prtkes(:, :, :) )
+if ( lbudget_rv ) call Budget_store_end( tbudgets(NBUDGET_RV ), 'ADV', prrs (:, :, :, 1) )
+if ( lbudget_rc ) call Budget_store_end( tbudgets(NBUDGET_RC ), 'ADV', prrs (:, :, :, 2) )
+if ( lbudget_rr ) call Budget_store_end( tbudgets(NBUDGET_RR ), 'ADV', prrs (:, :, :, 3) )
+if ( lbudget_ri ) call Budget_store_end( tbudgets(NBUDGET_RI ), 'ADV', prrs (:, :, :, 4) )
+if ( lbudget_rs ) call Budget_store_end( tbudgets(NBUDGET_RS ), 'ADV', prrs (:, :, :, 5) )
+if ( lbudget_rg ) call Budget_store_end( tbudgets(NBUDGET_RG ), 'ADV', prrs (:, :, :, 6) )
+if ( lbudget_rh ) call Budget_store_end( tbudgets(NBUDGET_RH ), 'ADV', prrs (:, :, :, 7) )
+if ( lbudget_sv) then
+ do jsv = 1, ksv
+ call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + jsv ), 'ADV', prsvs(:, :, :, jsv) )
+ end do
+end if
+
+! Remove non-physical negative values (unnecessary in a perfect world) + corresponding budgets
+call Sources_neg_correct( hcloud, 'NEADV', krr, ptstep, ppabst, ptht, prt, prths, prrs, prsvs )
-END IF
-!
#ifdef MNH_OPENACC
CALL MNH_REL_ZT3D(IZ1, IZ2)
#endif
-!-------------------------------------------------------------------------------
IF (MPPDB_INITIALIZED) THEN
!Check all INOUT arrays
diff --git a/src/ZSOLVER/advection_uvw_cen.f90 b/src/ZSOLVER/advection_uvw_cen.f90
index 534b5ac01..fd980ac36 100644
--- a/src/ZSOLVER/advection_uvw_cen.f90
+++ b/src/ZSOLVER/advection_uvw_cen.f90
@@ -89,38 +89,35 @@ END MODULE MODI_ADVECTION_UVW_CEN
!! Original 01/2013 (from ADVECTION routine)
!! Modif
!! J.Escobar 21/03/2013: for HALOK comment all NHALO=1 test
-! P. Wautelet 20/05/2019: add name argument to ADDnFIELD_ll + new ADD4DFIELD_ll subroutine
-!
+! P. Wautelet 20/05/2019: add name argument to ADDnFIELD_ll + new ADD4DFIELD_ll subroutine
+! P. Wautelet 02/2020: use the new data structures and subroutines for budgets
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
! ------------
!
-USE MODE_ll
-use mode_mppdb
USE MODD_ARGSLIST_ll, ONLY : LIST_ll, HALO2LIST_ll
USE MODD_CONF
-USE MODD_PARAMETERS
+use modd_budget, only: lbudget_u, lbudget_v, lbudget_w, NBUDGET_U, NBUDGET_V, NBUDGET_W, tbudgets
USE MODD_GRID_n
-!
-#ifndef MNH_OPENACC
-USE MODI_SHUMAN
-#else
-USE MODI_SHUMAN_DEVICE
-#endif
-USE MODI_CONTRAV
-USE MODI_ADVECUVW_2ND
-USE MODI_ADVECUVW_4TH
-!
-USE MODD_BUDGET
-USE MODI_BUDGET
-!
+
+use mode_budget, only: Budget_store_init, Budget_store_end
+USE MODE_ll
#ifdef MNH_OPENACC
USE MODE_DEVICE
USE MODE_MNH_ZWORK, ONLY : ZT3D, MNH_GET_ZT3D , MNH_REL_ZT3D
use mode_msg
#endif
-!
+use mode_mppdb
+
+USE MODI_ADVECUVW_2ND
+USE MODI_ADVECUVW_4TH
+USE MODI_CONTRAV
+#ifndef MNH_OPENACC
+USE MODI_SHUMAN
+#else
+USE MODI_SHUMAN_DEVICE
+#endif
!-------------------------------------------------------------------------------
!
IMPLICIT NONE
@@ -212,6 +209,10 @@ IF (MPPDB_INITIALIZED) THEN
CALL MPPDB_CHECK(PRWS,"ADVECTION_UVW_CEN beg:PRWS")
END IF
+if ( lbudget_u ) call Budget_store_init( tbudgets(NBUDGET_U), 'ADV', prus(:, :, :) )
+if ( lbudget_v ) call Budget_store_init( tbudgets(NBUDGET_V), 'ADV', prvs(:, :, :) )
+if ( lbudget_w ) call Budget_store_init( tbudgets(NBUDGET_W), 'ADV', prws(:, :, :) )
+
allocate( zus ( size( put, 1 ), size( put, 2 ), size( put, 3 ) ) )
allocate( zvs ( size( put, 1 ), size( put, 2 ), size( put, 3 ) ) )
allocate( zws ( size( put, 1 ), size( put, 2 ), size( put, 3 ) ) )
@@ -250,7 +251,7 @@ CALL INIT_ON_HOST_AND_DEVICE(ZMZM_RHODJ,-3e97,'ADVECTION_UVW_CEN::ZMZM_RHODJ')
!
CALL MNH_GET_ZT3D(IZ1, IZ2)
#endif
-!
+
#ifndef MNH_OPENACC
ZMXM_RHODJ = MXM(PRHODJ)
ZMYM_RHODJ = MYM(PRHODJ)
@@ -345,11 +346,11 @@ PDUM(:,:,:) = ZUS(:,:,:) - PUM(:,:,:)
PDVM(:,:,:) = ZVS(:,:,:) - PVM(:,:,:)
PDWM(:,:,:) = ZWS(:,:,:) - PWM(:,:,:)
!$acc end kernels
-!
-IF (LBUDGET_U) CALL BUDGET (PRUS,1,'ADV_BU_RU')
-IF (LBUDGET_V) CALL BUDGET (PRVS,2,'ADV_BU_RV')
-IF (LBUDGET_W) CALL BUDGET (PRWS,3,'ADV_BU_RW')
-!
+
+if ( lbudget_u ) call Budget_store_end( tbudgets(NBUDGET_U), 'ADV', prus(:, :, :) )
+if ( lbudget_v ) call Budget_store_end( tbudgets(NBUDGET_V), 'ADV', prvs(:, :, :) )
+if ( lbudget_w ) call Budget_store_end( tbudgets(NBUDGET_W), 'ADV', prws(:, :, :) )
+
#ifdef MNH_OPENACC
CALL MNH_REL_ZT3D(IZ1, IZ2)
#endif
diff --git a/src/ZSOLVER/gdiv.f90 b/src/ZSOLVER/gdiv.f90
index c78da2ed6..b24671946 100644
--- a/src/ZSOLVER/gdiv.f90
+++ b/src/ZSOLVER/gdiv.f90
@@ -151,7 +151,6 @@ INTEGER :: IKE ! indice K for the last inner mass point along z
INTEGER :: JI,JJ,JK ! loop indexes
!
INTEGER :: IIU,IJU,IKU
-#ifdef MNH_OPENACC
!
#ifdef MNH_OPENACC
REAL, DIMENSION(:,:,:) , POINTER , CONTIGUOUS :: ZTMP1,ZTMP2
diff --git a/src/ZSOLVER/pressure.f90 b/src/ZSOLVER/pressure.f90
deleted file mode 100644
index 7c11f6593..000000000
--- a/src/ZSOLVER/pressure.f90
+++ /dev/null
@@ -1,689 +0,0 @@
-!MNH_LIC Copyright 1994-2020 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.
-!-----------------------------------------------------------------
-!###################
-MODULE MODI_PRESSURE
-!###################
-!
-INTERFACE
-!
- SUBROUTINE PRESSURE( &
- HLBCX,HLBCY,HPRESOPT,KITR,OITRADJ,KTCOUNT,PRELAX,KMI, &
- PRHODJ,PDXX,PDYY,PDZZ,PDZX,PDZY,PDXHATM,PDYHATM,PRHOM, &
- PAF,PBF,PCF,PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,PPABSM, &
- KRR,KRRL,KRRI,PDRYMASST,PREFMASS,PMASS_O_PHI0, &
- PTHT,PRT,PRHODREF,PTHVREF,PRVREF,PEXNREF,PLINMASS, &
- PRUS,PRVS,PRWS,PPABST,PRESIDUAL)
-!
-IMPLICIT NONE
-!
-CHARACTER (LEN=*), DIMENSION(:), INTENT(IN) :: HLBCX ! x-direction LBC type
-CHARACTER (LEN=*), DIMENSION(:), INTENT(IN) :: HLBCY ! y-direction LBC type
-!
-CHARACTER (LEN=5), INTENT(IN) :: HPRESOPT ! choice of the pressure solver
-!
-INTEGER, INTENT(INOUT) :: KITR ! number of iterations for the
- ! pressure solver
-LOGICAL, INTENT(IN) :: OITRADJ ! switch to adjust or not KITR
-INTEGER, INTENT(IN) :: KTCOUNT ! counter value of the
- ! model temporal loop
-INTEGER, INTENT(IN) :: KMI ! Model index
-REAL, INTENT(IN) :: PRELAX ! relaxation coefficient for
- ! the Richardson's method
-!
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! density of reference state
- ! * J
-!
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX,PDYY,PDZZ
-REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PDZX,PDZY ! metric coefficients
-!
-REAL, INTENT(IN) :: PDXHATM ! mean grid increment in the x
- ! direction
-REAL, INTENT(IN) :: PDYHATM ! mean grid increment in the y
- ! direction
-!
-REAL, DIMENSION (:), INTENT(IN) :: PRHOM ! mean of XRHODJ on the plane x y
- ! localized at a mass level
-!
-REAL, DIMENSION(:), INTENT(IN) :: PAF,PCF ! vectors giving the nonvanishing
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PBF ! elements of the tri-diag.
- ! matrix in the pressure eq.
-!
- ! arrays of sin or cos values
- ! for the FFT :
-REAL, DIMENSION(:), INTENT(IN) :: PTRIGSX ! - along x
-REAL, DIMENSION(:), INTENT(IN) :: PTRIGSY ! - along y
-!
- ! decomposition in prime
- ! numbers for the FFT:
-INTEGER, DIMENSION(19), INTENT(IN) :: KIFAXX ! - along x
-INTEGER, DIMENSION(19), INTENT(IN) :: KIFAXY ! - along y
-!
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PPABSM ! pressure (t-dt)
-!
-INTEGER, INTENT(IN) :: KRR ! Total number of water var.
-INTEGER, INTENT(IN) :: KRRL ! Number of liquid water var.
-INTEGER, INTENT(IN) :: KRRI ! Number of ice water var.
-!
-REAL, INTENT(IN) :: PDRYMASST ! Mass of dry air and of
-REAL, INTENT(IN) :: PREFMASS ! the ref. atmosphere
-REAL, INTENT(IN) :: PMASS_O_PHI0 ! Mass / Phi0
-!
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHT ! Temperature and water
-REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRT ! variables at time t
-!
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODREF ! dry Density
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHVREF ! Virtual Temperature
- ! of the reference state
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PRVREF ! mixing ratio of the
- ! reference state
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PEXNREF ! Exner function
- ! of the reference state
-REAL, INTENT(IN) :: PLINMASS ! lineic mass through
- ! open boundaries
-!
-REAL, INTENT(INOUT) :: PRUS(:,:,:) ! source term along x
-REAL, INTENT(INOUT) :: PRVS(:,:,:) ! source term along y
-REAL, INTENT(INOUT) :: PRWS(:,:,:) ! source term along z
-!
-REAL, INTENT(INOUT) :: PPABST(:,:,:) ! pressure(t)
-!JUAN
-REAL, OPTIONAL :: PRESIDUAL
-!JUAN
-!
-END SUBROUTINE PRESSURE
-!
-END INTERFACE
-!
-END MODULE MODI_PRESSURE
-! ######################################################################
- SUBROUTINE PRESSURE( &
- HLBCX,HLBCY,HPRESOPT,KITR,OITRADJ,KTCOUNT,PRELAX,KMI, &
- PRHODJ,PDXX,PDYY,PDZZ,PDZX,PDZY,PDXHATM,PDYHATM,PRHOM, &
- PAF,PBF,PCF,PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,PPABSM, &
- KRR,KRRL,KRRI,PDRYMASST,PREFMASS,PMASS_O_PHI0, &
- PTHT,PRT,PRHODREF,PTHVREF,PRVREF,PEXNREF,PLINMASS, &
- PRUS,PRVS,PRWS,PPABST,PRESIDUAL)
-! ######################################################################
-!
-!!**** *PRESSURE * - solve the pressure equation and add the pressure term
-!! to the sources
-!!
-!! PURPOSE
-!! -------
-! The purpose of this routine is to solve the pressure equation:
-! with either the conjugate gradient method or the Richardson's method.
-! The pressure gradient is added to the sources in order
-! to nullify the divergence of the momentum* Thetavref*(1+Rvref)
-! at the time t+dt.
-!
-!!** METHOD
-!! ------
-!! The divergence of the sources ( RHS of the pressure equation ) is
-!! computed. The pressure equation is then solved by either CG method,
-!! either Richardson's method, or an exact method. Finally, the pressure
-!! gradient is added to the sources RUS, RVS, RWS.
-!! Finally, the absolute pressure is diagnozed from the total mass
-!! included in the simulation domain.
-!!
-!! EXTERNAL
-!! --------
-!! Subroutine MASS_LEAK : assures global non-divergence condition in the
-!! case of open boundaries
-!! Subroutine FLAT_INV : solve the pressure equation for the case
-!! without orography
-!! Subroutine RICHARDSON: solve the pressure equation with the
-!! Richardson's method
-!! Subroutine CONJGRAD : solve the pressure equation with the Conjugate
-!! Gradient algorithm
-!! Function GX_M_U : compute the gradient along x
-!! Function GY_M_V : compute the gradient along y
-!! Function GZ_M_W : compute the gradient along z
-!! Subroutine GDIV : compute J times the divergence of 1/J times a vector
-!! Function MXM: compute an average in the x direction for a variable
-!! at a mass localization
-!! Function MYM: compute an average in the y direction for a variable
-!! at a mass localization
-!! Function MZM: compute an average in the z direction for a variable
-!! at a mass localization
-!! Subroutine P_ABS : compute the constant for PABS and therefore, the
-!! absolute pressure function
-!!
-!! IMPLICIT ARGUMENTS
-!! ------------------
-!! Module MODD_CONF: model configuration
-!! LFLAT: logical switch for zero orography
-!! L2D : logical switch for two-dimensional configuration
-!! LCARTESIAN : logical switch for cartesian geometry
-!! Module MODD_PARAMETERS: declaration of parameter variables
-!! JPHEXT, JPVEXT: define the number of marginal points out of the
-!! physical domain along horizontal and vertical directions respectively
-!! Module MODD_CST: physical constants
-!! XCPD
-!!
-!! REFERENCE
-!! ---------
-!! Book2 of documentation (subroutine PRESSURE) + Book1 ( )
-!!
-!! AUTHOR
-!! ------
-!! P. Hereil and J. Stein * Meteo France *
-!!
-!! MODIFICATIONS
-!! -------------
-!! Original 05/07/94
-!! Modification 03/01/95 (Lafore) To add the absolute pressure diagnosis
-!! Modification 31/01/95 (Stein) Copy of the pressure function in the
-!! 2D case in the two outermost planes
-!! Modification 16/02/95 (Mallet) Add the call to MASS_LEAK
-!! Modification 16/03/95 (Stein) change the argument list of the
-!! gradient and remove R from the historical var.
-!! Modification 30/06/95 (Stein) Add a test not to compute the absolute
-!! pressure in the Boussinesq case
-!! 16/10/95 (J. Stein) change the budget calls
-!! 29/01/96 (J. Stein) call iterative resolution for
-!! non-cartessian geometry
-!! 19/12/96 (J.-P. Pinty) update the budget calls
-!! 14/01/97 (Stein,Lafore) New anelastic equations
-!! 17/12/97 ( Stein )include the case of non-vanishing
-!! orography at the lbc
-!! 26/03/98 (Stein,Jabouille) fix the value of the corner point
-!! 15/06/98 (D.Lugato, R.Guivarch) Parallelisation
-!! 25/08/99 (J.-P. Pinty) add CRESI option to CPRESOPT
-!! 06/11/02 (V. Masson) update the budget calls
-!! 24/08/2005 (J. escobar) BUG : remove IIE+1, IJE+1 out of bound
-!! references in parallel run
-!! 08/2010 (V.Masson, C.Lac) Add UPDATE_HALO
-!! 11/2010 (V.Masson, C.Lac) PPABST, must not be cyclic => add temp array
-!! to save it before UPDATE_HALO
-!! 06/2011 (J.escobar ) Bypass Bug with ifort11/12 on HLBCX,HLBCY
-!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
-! P. Wautelet 20/05/2019: add name argument to ADDnFIELD_ll + new ADD4DFIELD_ll subroutine
-!
-!-------------------------------------------------------------------------------
-!
-!* 0. DECLARATIONS
-! ------------
-!
-USE MODD_PARAMETERS
-USE MODD_BUDGET
-USE MODD_CONF
-USE MODD_CST
-USE MODD_LUNIT_n, ONLY: TLUOUT
-USE MODI_MASS_LEAK
-USE MODI_GDIV
-USE MODI_FLAT_INV
-USE MODI_RICHARDSON
-USE MODI_CONJGRAD
-USE MODI_CONRESOL
-USE MODI_GRADIENT_M
-USE MODI_SHUMAN
-USE MODI_P_ABS
-USE MODI_BUDGET
-!
-USE MODD_ARGSLIST_ll, ONLY : LIST_ll
-USE MODE_ll
-!
-IMPLICIT NONE
-!
-!* 0.1 declarations of arguments
-!
- CHARACTER (LEN=*), DIMENSION(:), INTENT(IN) :: HLBCX ! x-direction LBC type
- CHARACTER (LEN=*), DIMENSION(:), INTENT(IN) :: HLBCY ! y-direction LBC type
-!
-CHARACTER (LEN=5), INTENT(IN) :: HPRESOPT ! choice of the pressure solver
-!
-INTEGER, INTENT(INOUT) :: KITR ! number of iterations for the
- ! pressure solver
-LOGICAL, INTENT(IN) :: OITRADJ ! switch to adjust or not KITR
-INTEGER, INTENT(IN) :: KTCOUNT ! counter value of the
- ! model temporal loop
-INTEGER, INTENT(IN) :: KMI ! Model index
-REAL, INTENT(IN) :: PRELAX ! relaxation coefficient for
- ! the Richardson's method
-!
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! density of reference state
- ! * J
-!
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PDXX,PDYY,PDZZ
-REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PDZX,PDZY ! metric coefficients
-!
-REAL, INTENT(IN) :: PDXHATM ! mean grid increment in the x
- ! direction
-REAL, INTENT(IN) :: PDYHATM ! mean grid increment in the y
- ! direction
-!
-REAL, DIMENSION (:), INTENT(IN) :: PRHOM ! mean of XRHODJ on the plane x y
- ! localized at a mass level
-!
-REAL, DIMENSION(:), INTENT(IN) :: PAF,PCF ! vectors giving the nonvanishing
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PBF ! elements of the tri-diag.
- ! matrix in the pressure eq.
-!
- ! arrays of sin or cos values
- ! for the FFT :
-REAL, DIMENSION(:), INTENT(IN) :: PTRIGSX ! - along x
-REAL, DIMENSION(:), INTENT(IN) :: PTRIGSY ! - along y
-!
- ! decomposition in prime
- ! numbers for the FFT:
-INTEGER, DIMENSION(19), INTENT(IN) :: KIFAXX ! - along x
-INTEGER, DIMENSION(19), INTENT(IN) :: KIFAXY ! - along y
-!
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PPABSM ! pressure (t-dt)
-!
-INTEGER, INTENT(IN) :: KRR ! Total number of water var.
-INTEGER, INTENT(IN) :: KRRL ! Number of liquid water var.
-INTEGER, INTENT(IN) :: KRRI ! Number of ice water var.
-!
-REAL, INTENT(IN) :: PDRYMASST ! Mass of dry air and of
-REAL, INTENT(IN) :: PREFMASS ! the ref. atmosphere
-REAL, INTENT(IN) :: PMASS_O_PHI0 ! Mass / Phi0
-!
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHT ! Temperature and water
-REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRT ! variables at time t
-!
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODREF ! dry Density
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHVREF ! Virtual Temperature
- ! of the reference state
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PRVREF ! mixing ratio of the
- ! reference state
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PEXNREF ! Exner function
- ! of the reference state
-REAL, INTENT(IN) :: PLINMASS ! lineic mass through
- ! open boundaries
-!
-REAL, INTENT(INOUT) :: PRUS(:,:,:) ! source term along x
-REAL, INTENT(INOUT) :: PRVS(:,:,:) ! source term along y
-REAL, INTENT(INOUT) :: PRWS(:,:,:) ! source term along z
-!
-REAL, INTENT(INOUT) :: PPABST(:,:,:) ! pressure(t)
-!JUAN
-REAL, OPTIONAL :: PRESIDUAL
-!JUAN
-!
-!
-!* 0.2 declarations of local variables
-!
-! Metric coefficients:
-!
-REAL, DIMENSION(SIZE(PPABSM,1),SIZE(PPABSM,2),SIZE(PPABSM,3)) :: ZDV_SOURCE
-! ! divergence of the sources
-!
-INTEGER :: IIB ! indice I for the first inner mass point along x
-INTEGER :: IIE ! indice I for the last inner mass point along x
-INTEGER :: IJB ! indice J for the first inner mass point along y
-INTEGER :: IJE ! indice J for the last inner mass point along y
-INTEGER :: IKB ! indice K for the first inner mass point along z
-INTEGER :: IKE ! indice K for the last inner mass point along z
-INTEGER :: ILUOUT ! Logical unit of output listing
-!
-REAL, DIMENSION(SIZE(PPABSM,1),SIZE(PPABSM,2),SIZE(PPABSM,3)) :: ZTHETAV, &
- ! virtual potential temperature
- ZPHIT
- ! MAE + DUR => Exner function perturbation
- ! LHE => Exner function perturbation * CPD * THVREF
-!
-REAL :: ZRV_OV_RD ! XRV / XRD
-REAL :: ZMAXVAL, ZMAXRES ! for print
-INTEGER, DIMENSION(3) :: IMAXLOC ! purpose
-INTEGER :: JWATER ! loop index on water species
-INTEGER :: IIU,IJU,IKU ! array sizes in I,J,K
-INTEGER :: JK ! loop index on the vertical levels
-INTEGER :: JI,JJ
-!
-REAL, DIMENSION(SIZE(PDXX,1),SIZE(PDXX,3)) :: ZPABS_S ! local pressure on southern side
-REAL, DIMENSION(SIZE(PDXX,1),SIZE(PDXX,3)) :: ZPABS_N ! local pressure on northern side
-REAL, DIMENSION(SIZE(PDYY,2),SIZE(PDXX,3)) :: ZPABS_E ! local pressure on eastern side
-REAL, DIMENSION(SIZE(PDYY,2),SIZE(PDXX,3)) :: ZPABS_W ! local pressure on western side
-INTEGER :: IINFO_ll
-TYPE(LIST_ll), POINTER :: TZFIELDS_ll ! list of fields to exchange
-TYPE(LIST_ll), POINTER :: TZFIELDS_2_ll ! list of fields to exchange
-!
-!
-!------------------------------------------------------------------------------
-!-------------------------------------------------------------------------------
-NULLIFY(TZFIELDS_ll)
-NULLIFY(TZFIELDS_2_ll)
-!
-!* 1. PRELIMINARIES
-! -------------
-!
-ILUOUT = TLUOUT%NLU
-!
-CALL GET_PHYSICAL_ll(IIB,IJB,IIE,IJE)
-CALL GET_DIM_EXT_ll('B',IIU,IJU)
-!
-IKB= 1+JPVEXT
-IKU= SIZE(PPABSM,3)
-IKE= IKU - JPVEXT
-!
-ZPABS_S(:,:) = 0.
-ZPABS_N(:,:) = 0.
-ZPABS_E(:,:) = 0.
-ZPABS_W(:,:) = 0.
-!
-!
-!-------------------------------------------------------------------------------
-!
-!* 3. COMPUTE THE LINEIC MASS
-! -----------------------
-!
-IF ( ANY(HLBCX(:)=='OPEN') .OR. ANY(HLBCY(:)=='OPEN') ) THEN
- CALL MASS_LEAK(PDXX,PDYY,HLBCX,HLBCY,PLINMASS,PRHODJ,PRUS,PRVS)
-END IF
-!
-!-------------------------------------------------------------------------------
-!
-!* 4. COMPUTE THE FORCING TERM FOR THE PRESSURE EQUATION
-! --------------------------------------------------
-!
-!
-CALL ADD3DFIELD_ll( TZFIELDS_ll, PRUS, 'PRESSURE::PRUS' )
-CALL ADD3DFIELD_ll( TZFIELDS_ll, PRVS, 'PRESSURE::PRVS' )
-CALL ADD3DFIELD_ll( TZFIELDS_ll, PRWS, 'PRESSURE::PRWS' )
-CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll)
-CALL CLEANLIST_ll(TZFIELDS_ll)
-!
-CALL GDIV(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRUS,PRVS,PRWS,ZDV_SOURCE)
-!
-! The non-homogenous Neuman problem is transformed in an homogenous Neuman
-! problem in the non-periodic cases
-IF (HLBCX(1) /= 'CYCL') THEN
- IF (LWEST_ll()) ZDV_SOURCE(IIB-1,:,:) = 0.
- IF (LEAST_ll()) ZDV_SOURCE(IIE+1,:,:) = 0.
-ENDIF
-!
-IF (.NOT. L2D .AND. HLBCY(1) /= 'CYCL') THEN
- IF (LSOUTH_ll()) ZDV_SOURCE(:,IJB-1,:) = 0.
- IF (LNORTH_ll()) ZDV_SOURCE(:,IJE+1,:) = 0.
-ENDIF
-!
-!-------------------------------------------------------------------------------
-!
-!* 5. SOLVE THE PRESSURE EQUATION
-! ---------------------------
-!
-!
-!* 5.1 Compute the virtual theta and the pressure perturbation
-! -------------------------------------------------------
-!
-IF(CEQNSYS=='MAE' .OR. CEQNSYS=='DUR') THEN
- IF(KRR > 0) THEN
- !
- ! compute the ratio : 1 + total water mass / dry air mass
- ZRV_OV_RD = XRV / XRD
- ZTHETAV(:,:,:) = 1. + PRT(:,:,:,1)
- DO JWATER = 2 , 1+KRRL+KRRI
- ZTHETAV(:,:,:) = ZTHETAV(:,:,:) + PRT(:,:,:,JWATER)
- END DO
- ! compute the virtual potential temperature when water is present in any
- ! form
- ZTHETAV(:,:,:) = PTHT(:,:,:) * (1. + PRT(:,:,:,1) * ZRV_OV_RD) / ZTHETAV(:,:,:)
- ELSE
- ! compute the virtual potential temperature when water is absent
- ZTHETAV(:,:,:) = PTHT(:,:,:)
- END IF
- !
- ZPHIT(:,:,:)=(PPABSM(:,:,:)/XP00)**(XRD/XCPD)-PEXNREF(:,:,:)
- !
-ELSEIF(CEQNSYS=='LHE') THEN
- ZPHIT(:,:,:)= ((PPABSM(:,:,:)/XP00)**(XRD/XCPD)-PEXNREF(:,:,:)) &
- * XCPD * PTHVREF(:,:,:)
- !
-END IF
-!
-IF(CEQNSYS=='LHE'.AND. LFLAT .AND. LCARTESIAN) THEN
- ! flat cartesian LHE case -> exact solution
- !
- CALL FLAT_INV(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF, &
- PTRIGSX,PTRIGSY,KIFAXX,KIFAXY,ZDV_SOURCE,ZPHIT)
-ELSE
- SELECT CASE(HPRESOPT)
- CASE('RICHA') ! Richardson's method
-!
- CALL RICHARDSON(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,ZTHETAV, &
- PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF,PTRIGSX,PTRIGSY, &
- KIFAXX,KIFAXY,KITR,KTCOUNT,PRELAX,ZDV_SOURCE,ZPHIT)
-!
- CASE('CGRAD') ! Conjugate Gradient method
- CALL CONJGRAD(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,ZTHETAV, &
- PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF,PTRIGSX,PTRIGSY, &
- KIFAXX,KIFAXY,KITR,ZDV_SOURCE,ZPHIT)
-!
- CASE('CRESI') ! Conjugate Residual method
- CALL CONRESOL(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRHODJ,ZTHETAV, &
- PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF,PTRIGSX,PTRIGSY, &
- KIFAXX,KIFAXY,KITR,ZDV_SOURCE,ZPHIT)
- END SELECT
-END IF
-!
-!-------------------------------------------------------------------------------
-!
-!* 6. ADD THE PRESSURE GRADIENT TO THE SOURCES
-! ----------------------------------------
-!
-IF ( HLBCX(1) /= 'CYCL' ) THEN
- IF(LWEST_ll()) ZPHIT(IIB-1,:,IKB-1) = ZPHIT(IIB,:,IKB-1)
- IF(LEAST_ll()) ZPHIT(IIE+1,:,IKB-1) = ZPHIT(IIE,:,IKB-1)
-ENDIF
-IF ( HLBCY(1) /= 'CYCL' ) THEN
- IF (LSOUTH_ll()) ZPHIT(:,IJB-1,IKB-1) = ZPHIT(:,IJB,IKB-1)
- IF (LNORTH_ll()) ZPHIT(:,IJE+1,IKB-1) = ZPHIT(:,IJE,IKB-1)
-ENDIF
-!
-IF ( L2D ) THEN
- IF (LSOUTH_ll()) ZPHIT(:,IJB-1,:) = ZPHIT(:,IJB,:)
- IF (LNORTH_ll()) ZPHIT(:,IJE+1,:) = ZPHIT(:,IJB,:)
-END IF
-!
-ZDV_SOURCE = GX_M_U(1,IKU,1,ZPHIT,PDXX,PDZZ,PDZX)
-!
-IF ( HLBCX(1) /= 'CYCL' ) THEN
- IF(LWEST_ll()) THEN
-!!!!!!!!!!!!!!!! FUJI compiler directive !!!!!!!!!!
-!!!!!!!!!!!!!!!! FUJI compiler directive !!!!!!!!!!
- DO JK=2,IKU-1
- ZDV_SOURCE(2,:,JK)= &
- (ZPHIT(2,:,JK) - ZPHIT(1,:,JK) - 0.5 * ( &
- PDZX(2,:,JK) * (ZPHIT(2,:,JK)-ZPHIT(2,:,JK-1)) / PDZZ(2,:,JK) &
- +PDZX(2,:,JK+1) * (ZPHIT(2,:,JK+1)-ZPHIT(2,:,JK)) / PDZZ(2,:,JK+1) &
- ) &
- ) / PDXX(2,:,JK)
- END DO
- ENDIF
- !
- IF(LEAST_ll()) THEN
- DO JK=2,IKU-1
- ZDV_SOURCE(IIU,:,JK)= &
- (ZPHIT(IIU,:,JK) - ZPHIT(IIU-1,:,JK) - 0.5 * ( &
- PDZX(IIU,:,JK) * (ZPHIT(IIU-1,:,JK)-ZPHIT(IIU-1,:,JK-1)) &
- / PDZZ(IIU-1,:,JK) &
- +PDZX(IIU,:,JK+1) * (ZPHIT(IIU-1,:,JK+1)-ZPHIT(IIU-1,:,JK)) &
- / PDZZ(IIU-1,:,JK+1) &
- ) &
- ) / PDXX(IIU,:,JK)
- END DO
- END IF
-END IF
-!
-IF(CEQNSYS=='MAE' .OR. CEQNSYS=='DUR') THEN
- PRUS = PRUS - MXM(PRHODJ * XCPD * ZTHETAV) * ZDV_SOURCE
- PRWS = PRWS - MZM(PRHODJ * XCPD * ZTHETAV) * GZ_M_W(1,IKU,1,ZPHIT,PDZZ)
-ELSEIF(CEQNSYS=='LHE') THEN
- PRUS = PRUS - MXM(PRHODJ) * ZDV_SOURCE
- PRWS = PRWS - MZM(PRHODJ) * GZ_M_W(1,IKU,1,ZPHIT,PDZZ)
-END IF
-!
-IF(.NOT. L2D) THEN
-!
- ZDV_SOURCE = GY_M_V(1,IKU,1,ZPHIT,PDYY,PDZZ,PDZY)
-!
- IF ( HLBCY(1) /= 'CYCL' ) THEN
- IF (LSOUTH_ll()) THEN
-!!!!!!!!!!!!!!!! FUJI compiler directive !!!!!!!!!!
-!!!!!!!!!!!!!!!! FUJI compiler directive !!!!!!!!!!
- DO JK=2,IKU-1
- ZDV_SOURCE(:,2,JK)= &
- (ZPHIT(:,2,JK) - ZPHIT(:,1,JK) - 0.5 * ( &
- PDZY(:,2,JK) * (ZPHIT(:,2,JK)-ZPHIT(:,2,JK-1)) / PDZZ(:,2,JK) &
- +PDZY(:,2,JK+1) * (ZPHIT(:,2,JK+1)-ZPHIT(:,2,JK)) / PDZZ(:,2,JK+1) &
- ) &
- ) / PDYY(:,2,JK)
- END DO
- END IF
- !
- IF (LNORTH_ll()) THEN
- DO JK=2,IKU-1
- ZDV_SOURCE(:,IJU,JK)= &
- (ZPHIT(:,IJU,JK) - ZPHIT(:,IJU-1,JK) - 0.5 * ( &
- PDZY(:,IJU,JK) * (ZPHIT(:,IJU-1,JK)-ZPHIT(:,IJU-1,JK-1)) &
- / PDZZ(:,IJU-1,JK) &
- +PDZY(:,IJU,JK+1) * (ZPHIT(:,IJU-1,JK+1)-ZPHIT(:,IJU-1,JK)) &
- / PDZZ(:,IJU-1,JK+1) &
- ) &
- ) / PDYY(:,IJU,JK)
- END DO
- END IF
- END IF
-!
- IF(CEQNSYS=='MAE' .OR. CEQNSYS=='DUR') THEN
- PRVS = PRVS - MYM(PRHODJ * XCPD * ZTHETAV) * ZDV_SOURCE
- ELSEIF(CEQNSYS=='LHE') THEN
- PRVS = PRVS - MYM(PRHODJ) * ZDV_SOURCE
- END IF
-END IF
-!
-!! same boundary conditions as in gdiv ... !! (provisory coding)
-!! (necessary when NVERB=1)
-DO JJ = IJB,IJE ! copy the horizontal components under
- DO JI = IIB,IIE
- PRUS(JI,JJ,IKB-1)=PRUS(JI,JJ,IKB)
- PRUS(JI,JJ,IKE+1)=PRUS(JI,JJ,IKE)
- END DO
-END DO
-!
-DO JJ = IJB,IJE
- DO JI = IIB,IIE ! the ground and above the top
- PRVS(JI,JJ,IKB-1)=PRVS(JI,JJ,IKB)
- PRVS(JI,JJ,IKE+1)=PRVS(JI,JJ,IKE)
- END DO
-END DO
-!!
-!
-! compute the residual divergence
-CALL GDIV(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRUS,PRVS,PRWS,ZDV_SOURCE)
-!
-IF ( CEQNSYS=='DUR' ) THEN
- IF ( SIZE(PRVREF,1) == 0 ) THEN
- ZDV_SOURCE=ZDV_SOURCE/PRHODJ/XTH00*PRHODREF*PTHVREF
- ELSE
- ZDV_SOURCE=ZDV_SOURCE/PRHODJ/XTH00*PRHODREF*PTHVREF*(1.+PRVREF)
- END IF
-ELSEIF( CEQNSYS=='MAE' .OR. CEQNSYS=='LHE' ) THEN
- ZDV_SOURCE=ZDV_SOURCE/PRHODJ*PRHODREF
-END IF
-!
-ZMAXVAL=MAX_ll(ABS(ZDV_SOURCE),IINFO_ll)
-IF (PRESENT(PRESIDUAL)) PRESIDUAL = ZMAXVAL
-IMAXLOC=MAXLOC( ABS(ZDV_SOURCE(IIB:IIE,IJB:IJE,IKB:IKE))) !provisory coding one one processor only
-!
-WRITE(ILUOUT,*) 'residual divergence / 2 DT', ZMAXVAL, &
- ' located at ', IMAXLOC
-! number of iterations adjusted
-IF (LFLAT .AND. LCARTESIAN) THEN
- ZMAXRES = 1.E-12
-ELSE
- ZMAXRES = 1.E-9
-END IF
-!
-IF (OITRADJ) THEN
- IF (ZMAXVAL>10.*ZMAXRES) THEN
- KITR=KITR+2
- WRITE(ILUOUT,*) 'NITR adjusted to ', KITR
- ELSE IF (ZMAXVAL<ZMAXRES) THEN
- KITR=MAX(KITR-1,1)
- WRITE(ILUOUT,*) 'NITR adjusted to ', KITR
- ENDIF
-ENDIF
-!
-!* 7. STORAGE OF THE FIELDS IN BUDGET ARRAYS
-! --------------------------------------
-!
-IF (LBUDGET_U) CALL BUDGET (PRUS,1,'PRES_BU_RU')
-IF (LBUDGET_V) CALL BUDGET (PRVS,2,'PRES_BU_RV')
-IF (LBUDGET_W) CALL BUDGET (PRWS,3,'PRES_BU_RW')
-!
-!-------------------------------------------------------------------------------
-!
-!* 8. ABSOLUTE PRESSURE COMPUTATION
-! -----------------------------
-!
-!IF ( ABS(PRHODREF(IIB,IJB,IKB)-PRHODREF(IIB,IJB,IKE)) > 1.E-16 &
-IF ( ABS(PRHODREF(IIB,IJB,IKB)-PRHODREF(IIB,IJB,IKE)) > 1.E-12 &
- .AND. KTCOUNT >0 ) THEN
- CALL P_ABS ( KRR, KRRL, KRRI, PDRYMASST, PREFMASS, PMASS_O_PHI0, &
- PTHT, PRT, PRHODJ, PRHODREF, ZTHETAV, PTHVREF, &
- PRVREF, PEXNREF, ZPHIT )
-!
- IF(CEQNSYS=='MAE' .OR. CEQNSYS=='DUR') THEN
- PPABST(:,:,:)=XP00*(ZPHIT+PEXNREF)**(XCPD/XRD)
- ELSEIF(CEQNSYS=='LHE') THEN
- PPABST(:,:,:)=XP00*(ZPHIT/(XCPD*PTHVREF)+PEXNREF)**(XCPD/XRD)
- ENDIF
-!
- IF( HLBCX(1) == 'CYCL' ) THEN
- IF (LWEST_ll()) THEN
- ZPABS_W(:,:)= PPABST(IIB,:,:)
- END IF
-!
- IF (LEAST_ll()) THEN
- ZPABS_E(:,:)= PPABST(IIE+1,:,:)
- END IF
-!
- END IF
-!
- IF( HLBCY(1) == 'CYCL' ) THEN
- IF (LSOUTH_ll()) THEN
- ZPABS_S(:,:)= PPABST(:,IJB,:)
- END IF
-!
- IF (LNORTH_ll()) THEN
- ZPABS_N(:,:)= PPABST(:,IJE+1,:)
- END IF
-!
- END IF
-!
- CALL ADD3DFIELD_ll( TZFIELDS_2_ll, PPABST, 'PRESSURE::PPABST' )
- CALL UPDATE_HALO_ll(TZFIELDS_2_ll,IINFO_ll)
- CALL CLEANLIST_ll(TZFIELDS_2_ll)
-!
- IF( HLBCX(1) == 'CYCL' ) THEN
- IF (LWEST_ll()) THEN
- PPABST(IIB,:,:) = ZPABS_W(:,:)
- END IF
-!
- IF (LEAST_ll()) THEN
- PPABST(IIE+1,:,:) = ZPABS_E(:,:)
- END IF
-!
- END IF
-!
- IF( HLBCY(1) == 'CYCL' ) THEN
- IF (LSOUTH_ll()) THEN
- PPABST(:,IJB,:) = ZPABS_S(:,:)
- END IF
-!
- IF (LNORTH_ll()) THEN
- PPABST(:,IJE+1,:) = ZPABS_N(:,:)
- END IF
-!
- END IF
-!
-END IF
-!
-!-------------------------------------------------------------------------------
-!
-END SUBROUTINE PRESSURE
diff --git a/src/ZSOLVER/viscosity.f90 b/src/ZSOLVER/viscosity.f90
index a450ca9d7..e202cf20e 100644
--- a/src/ZSOLVER/viscosity.f90
+++ b/src/ZSOLVER/viscosity.f90
@@ -1,4 +1,4 @@
-!MNH_LIC Copyright 1994-2020 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC Copyright 1994-2021 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.
@@ -91,21 +91,29 @@ SUBROUTINE VISCOSITY(HLBCX, HLBCY, KRR, KSV, PNU, PPRANDTL, &
!! 01/18 (C.Lac) Add budgets
! P. Wautelet 20/05/2019: add name argument to ADDnFIELD_ll + new ADD4DFIELD_ll subroutine
! P. Wautelet 08/11/2019: corrected wrong budget name VISC_BU_RU -> VISC_BU_RTH
+! P. Wautelet 02/2020: use the new data structures and subroutines for budgets
+! T. Nagel 02/2021: add adhesion condition in case of an IBM-obstacle at the domain top boundary
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
! ------------
!
- USE MODI_LAP_M
- USE MODI_SHUMAN
- USE MODD_PARAMETERS
+ USE MODD_ARGSLIST_ll, ONLY: LIST_ll
+ use modd_budget, only: lbudget_u, lbudget_v, lbudget_w, lbudget_th, lbudget_rv, lbudget_rc, &
+ lbudget_rr, lbudget_ri, lbudget_rs, lbudget_rg, lbudget_rh, lbudget_sv, &
+ NBUDGET_U, NBUDGET_V, NBUDGET_W, NBUDGET_TH, NBUDGET_RV, NBUDGET_RC, &
+ NBUDGET_RR, NBUDGET_RI, NBUDGET_RS, NBUDGET_RG, NBUDGET_RH, NBUDGET_SV1, &
+ tbudgets
USE MODD_CONF
- USE MODD_VISCOSITY
USE MODD_DRAG_n
- USE MODD_BUDGET
+ USE MODD_PARAMETERS
+ USE MODD_VISCOSITY
+
+ use mode_budget, only: Budget_store_init, Budget_store_end
USE MODE_ll
- USE MODD_ARGSLIST_ll, ONLY : LIST_ll
- USE MODI_BUDGET
+
+ USE MODI_SHUMAN
+ USE MODI_LAP_M
!
!-------------------------------------------------------------------------------
!
@@ -180,6 +188,23 @@ IIU=SIZE(PWT,1)
IJU=SIZE(PWT,2)
IKU=SIZE(PWT,3)
+if ( lbudget_u .and. ovisc_uvw ) call Budget_store_init( tbudgets(NBUDGET_U ), 'VISC', prus (:, :, :) )
+if ( lbudget_v .and. ovisc_uvw ) call Budget_store_init( tbudgets(NBUDGET_V ), 'VISC', prvs (:, :, :) )
+if ( lbudget_w .and. ovisc_uvw ) call Budget_store_init( tbudgets(NBUDGET_W ), 'VISC', prws (:, :, :) )
+if ( lbudget_th .and. ovisc_th ) call Budget_store_init( tbudgets(NBUDGET_TH), 'VISC', prths(:, :, :) )
+if ( lbudget_rv .and. ovisc_r ) call Budget_store_init( tbudgets(NBUDGET_RV), 'VISC', prrs (:, :, :, 1) )
+if ( lbudget_rc .and. ovisc_r ) call Budget_store_init( tbudgets(NBUDGET_RC), 'VISC', prrs (:, :, :, 2) )
+if ( lbudget_rr .and. ovisc_r ) call Budget_store_init( tbudgets(NBUDGET_RR), 'VISC', prrs (:, :, :, 3) )
+if ( lbudget_ri .and. ovisc_r ) call Budget_store_init( tbudgets(NBUDGET_RI), 'VISC', prrs (:, :, :, 4) )
+if ( lbudget_rs .and. ovisc_r ) call Budget_store_init( tbudgets(NBUDGET_RS), 'VISC', prrs (:, :, :, 5) )
+if ( lbudget_rg .and. ovisc_r ) call Budget_store_init( tbudgets(NBUDGET_RG), 'VISC', prrs (:, :, :, 6) )
+if ( lbudget_rh .and. ovisc_r ) call Budget_store_init( tbudgets(NBUDGET_RH), 'VISC', prrs (:, :, :, 7) )
+if ( lbudget_sv .and. ovisc_sv ) then
+ do ik = 1, ksv
+ call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + ik), 'VISC', prsvs(:, :, :, ik) )
+ end do
+end if
+
!* 1. Viscous forcing for potential temperature
! -----------------------------------------
!
@@ -193,8 +218,6 @@ IF (OVISC_TH) THEN
!
END IF
!
-IF (LBUDGET_TH) CALL BUDGET (PRTHS,4,'VISC_BU_RTH')
-!
!-------------------------------------------------------------------------------
!
!* 2. Viscous forcing for moisture
@@ -211,14 +234,6 @@ IF (OVISC_R .AND. (SIZE(PRT,1) > 0)) THEN
!
END IF
!
-IF (LBUDGET_RV) CALL BUDGET (PRRS(:,:,:,1),6,'VISC_BU_RRV')
-IF (LBUDGET_RC) CALL BUDGET (PRRS(:,:,:,2),7,'VISC_BU_RRC')
-IF (LBUDGET_RR) CALL BUDGET (PRRS(:,:,:,3),8,'VISC_BU_RRR')
-IF (LBUDGET_RI) CALL BUDGET (PRRS(:,:,:,4),9,'VISC_BU_RRI')
-IF (LBUDGET_RS) CALL BUDGET (PRRS(:,:,:,5),10,'VISC_BU_RRS')
-IF (LBUDGET_RG) CALL BUDGET (PRRS(:,:,:,6),11,'VISC_BU_RRG')
-IF (LBUDGET_RH) CALL BUDGET (PRRS(:,:,:,7),12,'VISC_BU_RRH')
-!
!-------------------------------------------------------------------------------
!
!* 3. Viscous forcing for passive scalars
@@ -234,13 +249,6 @@ IF (OVISC_SV .AND. (SIZE(PSVT,1) > 0)) THEN
!
END IF
!
-IF (LBUDGET_SV) THEN
- DO IK = 1, KSV
- CALL BUDGET (PRSVS(:,:,:,IK), 12+IK, 'VISC_BU_RSV')
- END DO
-END IF
-!
-
!-------------------------------------------------------------------------------
!
!* 4. Viscous forcing for momentum
@@ -254,6 +262,8 @@ IF (OVISC_UVW) THEN
ZY1 = MXF(PUT)
IF (ODRAG) THEN
ZY1(:,:,1) = PDRAG * ZY1(:,:,2)
+!!Add adhesion condition in case of an IBM-obstacle at the domain top boundary
+! ZY1(:,:,IKU) = PDRAG * ZY1(:,:,IKE)
ENDIF
!
!
@@ -333,9 +343,22 @@ ENDIF
ENDIF
ENDIF
END IF
-!
-IF (LBUDGET_U) CALL BUDGET (PRUS,1,'VISC_BU_RU')
-IF (LBUDGET_V) CALL BUDGET (PRVS,2,'VISC_BU_RV')
-IF (LBUDGET_W) CALL BUDGET (PRWS,2,'VISC_BU_RW')
-!
+
+if ( lbudget_u .and. ovisc_uvw ) call Budget_store_end( tbudgets(NBUDGET_U ), 'VISC', prus (:, :, :) )
+if ( lbudget_v .and. ovisc_uvw ) call Budget_store_end( tbudgets(NBUDGET_V ), 'VISC', prvs (:, :, :) )
+if ( lbudget_w .and. ovisc_uvw ) call Budget_store_end( tbudgets(NBUDGET_W ), 'VISC', prws (:, :, :) )
+if ( lbudget_th .and. ovisc_th ) call Budget_store_end( tbudgets(NBUDGET_TH), 'VISC', prths(:, :, :) )
+if ( lbudget_rv .and. ovisc_r ) call Budget_store_end( tbudgets(NBUDGET_RV), 'VISC', prrs (:, :, :, 1) )
+if ( lbudget_rc .and. ovisc_r ) call Budget_store_end( tbudgets(NBUDGET_RC), 'VISC', prrs (:, :, :, 2) )
+if ( lbudget_rr .and. ovisc_r ) call Budget_store_end( tbudgets(NBUDGET_RR), 'VISC', prrs (:, :, :, 3) )
+if ( lbudget_ri .and. ovisc_r ) call Budget_store_end( tbudgets(NBUDGET_RI), 'VISC', prrs (:, :, :, 4) )
+if ( lbudget_rs .and. ovisc_r ) call Budget_store_end( tbudgets(NBUDGET_RS), 'VISC', prrs (:, :, :, 5) )
+if ( lbudget_rg .and. ovisc_r ) call Budget_store_end( tbudgets(NBUDGET_RG), 'VISC', prrs (:, :, :, 6) )
+if ( lbudget_rh .and. ovisc_r ) call Budget_store_end( tbudgets(NBUDGET_RH), 'VISC', prrs (:, :, :, 7) )
+if ( lbudget_sv .and. ovisc_sv ) then
+ do ik = 1, ksv
+ call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + ik), 'VISC', prsvs(:, :, :, ik) )
+ end do
+end if
+
END SUBROUTINE VISCOSITY
--
GitLab