diff --git a/src/ZSOLVER/SURCOUCHE/mode_exchange2_ll.f90 b/src/ZSOLVER/SURCOUCHE/mode_exchange2_ll.f90 index 8d2affca44ef628c6003ff27fbe700901b5bba7e..d4942e8ce3a03d5e30b30bbf6a4488eaeb794cf9 100644 --- a/src/ZSOLVER/SURCOUCHE/mode_exchange2_ll.f90 +++ b/src/ZSOLVER/SURCOUCHE/mode_exchange2_ll.f90 @@ -1,4 +1,4 @@ -!MNH_LIC Copyright 1998-2019 CNRS, Meteo-France and Universite Paul Sabatier +!MNH_LIC Copyright 1998-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. @@ -134,8 +134,7 @@ implicit none ! !* 1.1 Allocate the current HALO2_ll ! - ALLOCATE(TZHALO2LIST%HALO2) - + ALLOCATE(TZHALO2LIST%HALO2) ALLOCATE(TZHALO2LIST%HALO2%WEST(KDIMY, KDIMZ)) ALLOCATE(TZHALO2LIST%HALO2%EAST(KDIMY, KDIMZ)) ALLOCATE(TZHALO2LIST%HALO2%SOUTH(KDIMX, KDIMZ)) @@ -494,8 +493,6 @@ implicit none USE MODD_MPIF ! IMPLICIT NONE - -! INCLUDE 'mpif.h' ! !* 0.1 declarations of arguments ! @@ -728,8 +725,6 @@ implicit none ! IMPLICIT NONE ! -! INCLUDE 'mpif.h' -! !* 0.1 declarations of arguments ! TYPE(CRSPD_ll), POINTER :: TPCRSPDSEND, TPCRSPDRECV diff --git a/src/ZSOLVER/advection_uvw.f90 b/src/ZSOLVER/advection_uvw.f90 index 5c92dba01d1c9910880d4c3c759bf75d12817e6c..e6f09679fd4ce76d6fd181abd87dc3df1bfff146 100644 --- a/src/ZSOLVER/advection_uvw.f90 +++ b/src/ZSOLVER/advection_uvw.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. @@ -92,34 +92,35 @@ END MODULE MODI_ADVECTION_UVW !! J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1 !! C.LAC 10/2016 : Add OSPLIT_WENO ! 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 MODD_ARGSLIST_ll, ONLY : LIST_ll, HALO2LIST_ll -USE MODD_PARAMETERS, ONLY : JPVEXT +use modd_budget, only: lbudget_u, lbudget_v, lbudget_w, NBUDGET_U, NBUDGET_V, NBUDGET_W, tbudgets USE MODD_CONF, ONLY : NHALO -USE MODD_BUDGET -! +USE MODD_PARAMETERS, ONLY : JPVEXT + +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_ALLOCATE_ZT3D, MNH_GET_ZT3D , MNH_REL_ZT3D, MNH_GET_ZT4D , MNH_REL_ZT4D +#endif +use mode_mppdb + +USE MODI_ADV_BOUNDARIES +USE MODI_ADVECUVW_RK +USE MODI_CONTRAV #ifndef MNH_OPENACC USE MODI_SHUMAN #else +USE MODI_GET_HALO, ONLY: GET_HALO_D USE MODI_SHUMAN_DEVICE #endif -USE MODI_CONTRAV -USE MODI_ADVECUVW_RK -USE MODI_ADV_BOUNDARIES -USE MODI_BUDGET -USE MODI_GET_HALO ! -#ifdef MNH_OPENACC -USE MODE_DEVICE -USE MODE_MNH_ZWORK, ONLY : ZT3D, MNH_GET_ZT3D , MNH_REL_ZT3D, MNH_GET_ZT4D , MNH_REL_ZT4D, & - MNH_ALLOCATE_ZT3D -#endif -use mode_mppdb ! !------------------------------------------------------------------------------- ! @@ -158,7 +159,9 @@ INTEGER :: IKE ! indice K End in z direction REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZRUT REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZRVT REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZRWT +#ifdef MNH_OPENACC INTEGER :: IZRUT,IZRVT,IZRWT +#endif ! cartesian ! components of ! momentum @@ -166,18 +169,24 @@ INTEGER :: IZRUT,IZRVT,IZRWT REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZRUCT REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZRVCT REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZRWCT +#ifdef MNH_OPENACC INTEGER :: IZRUCT,IZRVCT,IZRWCT +#endif ! contravariant ! components ! of momentum ! REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZU, ZV, ZW +#ifdef MNH_OPENACC INTEGER :: IZU, IZV, IZW +#endif ! Guesses at the end of the sub time step REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZRUS_OTHER REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZRVS_OTHER REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZRWS_OTHER +#ifdef MNH_OPENACC INTEGER :: IZRUS_OTHER,IZRVS_OTHER,IZRWS_OTHER +#endif ! Contribution of the RK time step REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZRUS_ADV REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZRVS_ADV @@ -185,8 +194,10 @@ REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZRWS_ADV REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZMXM_RHODJ REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZMYM_RHODJ REAL, DIMENSION(:,:,:), POINTER, CONTIGUOUS :: ZMZM_RHODJ +#ifdef MNH_OPENACC INTEGER :: IZRUS_ADV,IZRVS_ADV,IZRWS_ADV INTEGER :: IZMXM_RHODJ,IZMYM_RHODJ,IZMZM_RHODJ +#endif ! ! Momentum tendencies due to advection INTEGER :: ISPLIT ! Number of splitting loops @@ -200,7 +211,7 @@ TYPE(LIST_ll), POINTER :: TZFIELDS0_ll ! list of fields to exchange ! #ifdef MNH_OPENACC INTEGER :: ISPL, IZUT, IZVT, IZWT, IZ1, IZ2 -INTEGER :: IZRUSB, IZRUSE, IZRVSB, IZRVSE, IZRWSB, IZRWSE +INTEGER :: IZRUSB, IZRUSE, IZRVSB, IZRVSE, IZRWSB, IZRWSE, IIBMS, IIBME #endif ! INTEGER :: IIU,IJU,IKU @@ -332,6 +343,7 @@ CALL MNH_GET_ZT3D(IZUT, IZVT, IZWT, IZ1, IZ2) CALL MNH_GET_ZT4D(ISPL, IZRUSB, IZRUSE) CALL MNH_GET_ZT4D(ISPL, IZRVSB, IZRVSE) CALL MNH_GET_ZT4D(ISPL, IZRWSB, IZRWSE) +CALL MNH_GET_ZT4D(3, IIBMS, IIBME) #endif ! IKE = SIZE(PWT,3) - JPVEXT @@ -345,7 +357,11 @@ CALL MXM_DEVICE(PRHODJ,ZMXM_RHODJ) CALL MYM_DEVICE(PRHODJ,ZMYM_RHODJ) CALL MZM_DEVICE(PRHODJ,ZMZM_RHODJ) #endif -! + +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(:, :, :) ) + !------------------------------------------------------------------------------- ! !* 1. COMPUTES THE CONTRAVARIANT COMPONENTS @@ -489,7 +505,7 @@ DO JSPL=1,ISPLIT ) #else ,ZT3D(:,:,:,IZUT), ZT3D(:,:,:,IZVT), ZT3D(:,:,:,IZWT), & - ZT3D(:,:,:,IZRUSB:IZRUSE), ZT3D(:,:,:,IZRVSB:IZRVSE), ZT3D(:,:,:,IZRWSB:IZRWSE) ) + ZT3D(:,:,:,IZRUSB:IZRUSE), ZT3D(:,:,:,IZRVSB:IZRVSE), ZT3D(:,:,:,IZRWSB:IZRWSE), ZT3D(:,:,:,IIBMS:IIBME) ) #endif ! ! Tendencies on wind @@ -539,21 +555,14 @@ END DO !* 4. BUDGETS ! ------- ! -IF (LBUDGET_U) THEN -!$acc update self(PRUS) - CALL BUDGET (PRUS,1,'ADV_BU_RU') -END IF -IF (LBUDGET_V) THEN -!$acc update self(PRVS) - CALL BUDGET (PRVS,2,'ADV_BU_RV') -END IF -IF (LBUDGET_W) THEN -!$acc update self(PRWS) - CALL BUDGET (PRWS,3,'ADV_BU_RW') -END IF +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_ZT4D(3, IIBMS ) CALL MNH_REL_ZT4D(ISPL, IZRWSB) CALL MNH_REL_ZT4D(ISPL, IZRVSB) CALL MNH_REL_ZT4D(ISPL, IZRUSB) diff --git a/src/ZSOLVER/advecuvw_rk.f90 b/src/ZSOLVER/advecuvw_rk.f90 index 6bdfae2fe4342040acd6d061ac9fa988b5f72be6..1352e75c9ea4ac02ae9c7f82c1979c9382fd3845 100644 --- a/src/ZSOLVER/advecuvw_rk.f90 +++ b/src/ZSOLVER/advecuvw_rk.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. @@ -20,7 +20,7 @@ INTERFACE #ifndef MNH_OPENACC ) #else - ,ZUT, ZVT, ZWT, ZRUS, ZRVS, ZRWS) + , ZUT, ZVT, ZWT, ZRUS, ZRVS, ZRWS, ZIBM ) #endif ! CHARACTER(LEN=6), INTENT(IN) :: HUVW_ADV_SCHEME! to the selected @@ -49,8 +49,8 @@ REAL, DIMENSION(:,:,:), INTENT(IN) :: PRUS_OTHER , PRVS_OTHER, PRWS_OTHER ! ! tendencies from other processes #ifdef MNH_OPENACC ! Work arrays -REAL, DIMENSION(:,:,:) :: ZUT, ZVT, ZWT -REAL, DIMENSION(:,:,:,:) :: ZRUS,ZRVS,ZRWS +REAL, DIMENSION(:,:,:) :: ZUT, ZVT, ZWT +REAL, DIMENSION(:,:,:,:) :: ZRUS, ZRVS, ZRWS, ZIBM #endif ! ! @@ -72,7 +72,7 @@ END MODULE MODI_ADVECUVW_RK #ifndef MNH_OPENACC ) #else - ,ZUT, ZVT, ZWT, ZRUS, ZRVS, ZRWS) + , ZUT, ZVT, ZWT, ZRUS, ZRVS, ZRWS, ZIBM ) #endif ! ########################################################################## ! @@ -121,6 +121,7 @@ END MODULE MODI_ADVECUVW_RK !! C.Lac 10/16 : Correction on RK loop ! P. Wautelet 10/04/2019: replace ABORT and STOP calls by Print_msg ! P. Wautelet 20/05/2019: add name argument to ADDnFIELD_ll + new ADD4DFIELD_ll subroutine +! T. Nagel,F.Auguste 06/2021: add IBM !! !------------------------------------------------------------------------------- ! @@ -129,16 +130,21 @@ END MODULE MODI_ADVECUVW_RK ! USE MODD_ARGSLIST_ll, ONLY: LIST_ll, HALO2LIST_ll USE MODD_CONF, ONLY: NHALO +USE MODD_IBM_PARAM_n, ONLY: LIBM, CIBM_ADV, XIBM_LS, XIBM_EPSI USE MODD_PARAMETERS, ONLY: JPVEXT +USE MODD_SUB_MODEL_n, ONLY: XT_IBM_FORC ! USE MODE_ll USE MODE_MPPDB use mode_msg ! USE MODI_ADV_BOUNDARIES +USE MODI_ADVECUVW_2ND USE MODI_ADVECUVW_4TH USE MODI_ADVECUVW_WENO_K USE MODI_GET_HALO +USE MODI_IBM_FORCING_ADV +USE MODI_SECOND_MNH USE MODI_SHUMAN ! ! @@ -178,8 +184,8 @@ REAL, DIMENSION(:,:,:), INTENT(OUT) :: PRUS_ADV , PRVS_ADV, PRWS_ADV REAL, DIMENSION(:,:,:), INTENT(IN) :: PRUS_OTHER , PRVS_OTHER, PRWS_OTHER ! ! tendencies from other processes #ifdef MNH_OPENACC -REAL, DIMENSION(:,:,:) :: ZUT, ZVT, ZWT -REAL, DIMENSION(:,:,:,:) :: ZRUS,ZRVS,ZRWS +REAL, DIMENSION(:,:,:) :: ZUT, ZVT, ZWT +REAL, DIMENSION(:,:,:,:) :: ZRUS, ZRVS, ZRWS, ZIBM #endif ! ! @@ -195,7 +201,7 @@ INTEGER :: IKE ! indice K End in z direction REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZUT, ZVT, ZWT ! Intermediate Guesses inside the RK loop ! -REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZRUS,ZRVS,ZRWS +REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZRUS,ZRVS,ZRWS,ZIBM #endif ! Momentum tendencies due to advection REAL, DIMENSION(:,:), ALLOCATABLE :: ZBUT ! Butcher array coefficients @@ -224,10 +230,12 @@ TYPE(LIST_ll), POINTER :: TZFIELDS_ll ! list of fields to exchange TYPE(LIST_ll), POINTER :: TZFIELDS0_ll ! list of fields to exchange TYPE(LIST_ll), POINTER :: TZFIELDS4_ll ! list of fields to exchange ! +REAL :: ZTIME1,ZTIME2 !------------------------------------------------------------------------------- -!$acc data present( PU, PV, PW, PUT, PVT, PWT, PMXM_RHODJ, PMYM_RHODJ, PMZM_RHODJ, & -!$acc & PRUCT, PRVCT, PRWCT, PRUS_ADV, PRVS_ADV, PRWS_ADV, & -!$acc & PRUS_OTHER, PRVS_OTHER, PRWS_OTHER, ZUT, ZVT, ZWT, ZRUS, ZRVS, ZRWS ) +!$acc data present( PU, PV, PW, PUT, PVT, PWT, PMXM_RHODJ, PMYM_RHODJ, PMZM_RHODJ, & +!$acc & PRUCT, PRVCT, PRWCT, PRUS_ADV, PRVS_ADV, PRWS_ADV, & +!$acc & PRUS_OTHER, PRVS_OTHER, PRWS_OTHER, ZUT, ZVT, ZWT, & +!$acc & ZRUS, ZRVS, ZRWS, ZIBM ) IF (MPPDB_INITIALIZED) THEN !Check all IN arrays @@ -247,9 +255,13 @@ IF (MPPDB_INITIALIZED) THEN CALL MPPDB_CHECK(PRVS_OTHER,"ADVECUVW_RK beg:PRVS_OTHER") CALL MPPDB_CHECK(PRWS_OTHER,"ADVECUVW_RK beg:PRWS_OTHER") END IF + +#ifdef MNH_OPENACC +if ( LIBM ) call Print_msg( NVERB_FATAL, 'GEN', 'ADVECUVW_RK', 'OpenACC: LIBM=T not yet implemented' ) +#endif ! -!* 0. INITIALIZATION -! -------------- +!* 0. INITIALIZATION +! --------------------- ! #ifndef MNH_OPENACC allocate(ZUT(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3))) @@ -391,8 +403,25 @@ ALLOCATE(ZRUS(SIZE(PUT,1),SIZE(PUT,2),SIZE(PWT,3),ISPL)) ALLOCATE(ZRVS(SIZE(PUT,1),SIZE(PUT,2),SIZE(PWT,3),ISPL)) ALLOCATE(ZRWS(SIZE(PUT,1),SIZE(PUT,2),SIZE(PWT,3),ISPL)) #endif + +IF ( LIBM ) THEN +#ifndef MNH_OPENACC + ALLOCATE( ZIBM(SIZE(PUT,1), SIZE(PUT,2), SIZE(PWT,3), 3) ) +#endif +!$acc kernels + ZIBM(:,:,:,:) = 1. +!$acc end kernels +END IF ! !$acc kernels present(PRUS_ADV,PRVS_ADV,PRWS_ADV) present(ZUT,ZVT,ZWT) present(PU,PV,PW) +IF (LIBM .AND. CIBM_ADV=='FREEZE') THEN + + WHERE (XIBM_LS(:,:,:,2).GT.-XIBM_EPSI) ZIBM(:,:,:,1) = 0. + WHERE (XIBM_LS(:,:,:,3).GT.-XIBM_EPSI) ZIBM(:,:,:,2) = 0. + WHERE (XIBM_LS(:,:,:,4).GT.-XIBM_EPSI) ZIBM(:,:,:,3) = 0. + +ENDIF +! PRUS_ADV = 0. PRVS_ADV = 0. PRWS_ADV = 0. @@ -400,7 +429,7 @@ PRWS_ADV = 0. !------------------------------------------------------------------------------- ! !* 2. Wind guess before RK loop -! ------------------------- +! -------------------------------- ! ZUT = PU ZVT = PV @@ -430,64 +459,68 @@ ZRVS(:, :, :, : ) = 0. ZRWS(:, :, :, : ) = 0. !$acc end kernels -!Necessary to work around a PGI bug (19.10) -!because following update ZRUS(:,:,:,JS) are done on the WHOLE array -!$acc update self(ZRUS,ZRVS,ZRWS) - !------------------------------------------------------------------------------- ! !* 3. BEGINNING of Runge-Kutta loop -! ----------------------------- +! ------------------------------------ ! - DO JS = 1, ISPL +RKLOOP: DO JS = 1, ISPL ! #ifndef MNH_OPENACC - CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZUT, PUT, 'U' ) - CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZVT, PVT, 'V' ) - CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZWT, PWT, 'W' ) + CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZUT, PUT, 'U' ) + CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZVT, PVT, 'V' ) + CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZWT, PWT, 'W' ) #else - CALL ADV_BOUNDARIES_DEVICE (HLBCX, HLBCY, ZUT, PUT, 'U' ) - CALL ADV_BOUNDARIES_DEVICE (HLBCX, HLBCY, ZVT, PVT, 'V' ) - CALL ADV_BOUNDARIES_DEVICE (HLBCX, HLBCY, ZWT, PWT, 'W' ) + CALL ADV_BOUNDARIES_DEVICE (HLBCX, HLBCY, ZUT, PUT, 'U' ) + CALL ADV_BOUNDARIES_DEVICE (HLBCX, HLBCY, ZVT, PVT, 'V' ) + CALL ADV_BOUNDARIES_DEVICE (HLBCX, HLBCY, ZWT, PWT, 'W' ) #endif - ! + ! #ifndef MNH_OPENACC - CALL UPDATE_HALO_ll(TZFIELDMT_ll,IINFO_ll) - CALL UPDATE_HALO2_ll(TZFIELDMT_ll, TZHALO2MT_ll, IINFO_ll) + CALL UPDATE_HALO_ll(TZFIELDMT_ll,IINFO_ll) + CALL UPDATE_HALO2_ll(TZFIELDMT_ll, TZHALO2MT_ll, IINFO_ll) #else ! acc update self(ZUT,ZVT,ZWT) !!$ CALL GET_HALO_D(ZUT,HNAME='ZUT') !!$ CALL GET_HALO_D(ZVT,HNAME='ZVT') !!$ CALL GET_HALO_D(ZWT,HNAME='ZWT') !!$ CALL UPDATE_HALO2_ll(TZFIELDMT_ll, TZHALO2MT_ll, IINFO_ll) -IF (GFIRST_CALL_ADVECUVW_RK) THEN - GFIRST_CALL_ADVECUVW_RK = .FALSE. - NULLIFY(TZHALO2_UT,TZHALO2_VT,TZHALO2_WT) - CALL INIT_HALO2_ll(TZHALO2_UT,1,IIU,IJU,IKU) - CALL INIT_HALO2_ll(TZHALO2_VT,1,IIU,IJU,IKU) - CALL INIT_HALO2_ll(TZHALO2_WT,1,IIU,IJU,IKU) -END IF - CALL GET_HALO2_DF(ZUT,TZHALO2_UT,HNAME='ZUT') - CALL GET_HALO2_DF(ZVT,TZHALO2_VT,HNAME='ZVT') - CALL GET_HALO2_DF(ZWT,TZHALO2_WT,HNAME='ZWT') -! acc update device(ZUT,ZVT,ZWT) -#endif - + IF (GFIRST_CALL_ADVECUVW_RK) THEN + GFIRST_CALL_ADVECUVW_RK = .FALSE. + NULLIFY(TZHALO2_UT,TZHALO2_VT,TZHALO2_WT) + CALL INIT_HALO2_ll(TZHALO2_UT,1,IIU,IJU,IKU) + CALL INIT_HALO2_ll(TZHALO2_VT,1,IIU,IJU,IKU) + CALL INIT_HALO2_ll(TZHALO2_WT,1,IIU,IJU,IKU) + END IF - + CALL GET_HALO2_DF(ZUT,TZHALO2_UT,HNAME='ZUT') + CALL GET_HALO2_DF(ZVT,TZHALO2_VT,HNAME='ZVT') + CALL GET_HALO2_DF(ZWT,TZHALO2_WT,HNAME='ZWT') +! acc update device(ZUT,ZVT,ZWT) +#endif ! !* 4. Advection with WENO -! ------------------- +! -------------------------- ! -!!$TZHALO2_UT => TZHALO2MT_ll ! 1rst add3dfield in model_n -!!$TZHALO2_VT => TZHALO2MT_ll%NEXT ! 2nd add3dfield in model_n -!!$TZHALO2_WT => TZHALO2MT_ll%NEXT%NEXT ! 3rst add3dfield in model_n - +!$acc kernels + IF (LIBM .AND. CIBM_ADV=='LOWORD') THEN + ZIBM(:,:,:,1)=ZRUS(:,:,:,JS) + ZIBM(:,:,:,2)=ZRVS(:,:,:,JS) + ZIBM(:,:,:,3)=ZRWS(:,:,:,JS) + ENDIF +!$acc end kernels +! +#ifndef MNH_OPENACC +TZHALO2_UT => TZHALO2MT_ll ! 1rst add3dfield in model_n +TZHALO2_VT => TZHALO2MT_ll%NEXT ! 2nd add3dfield in model_n +TZHALO2_WT => TZHALO2MT_ll%NEXT%NEXT ! 3rst add3dfield in model_n +#endif + IF (HUVW_ADV_SCHEME=='WENO_K') THEN CALL ADVECUVW_WENO_K (HLBCX, HLBCY, KWENO_ORDER, ZUT, ZVT, ZWT, & PRUCT, PRVCT, PRWCT, & ZRUS(:,:,:,JS), ZRVS(:,:,:,JS), ZRWS(:,:,:,JS), & - TZHALO2_UT,TZHALO2_VT,TZHALO2_WT & + TZHALO2_UT, TZHALO2_VT, TZHALO2_WT & #ifndef MNH_OPENACC ) #else @@ -500,30 +533,68 @@ END IF TZHALO2MT_ll ) ENDIF ! - + IF (LIBM .AND. CIBM_ADV=='LOWORD') THEN + IF (HUVW_ADV_SCHEME=='WENO_K') THEN + CALL ADVECUVW_WENO_K (HLBCX, HLBCY, 3, ZUT, ZVT, ZWT, & + PRUCT, PRVCT, PRWCT, & + ZIBM(:,:,:,1), ZIBM(:,:,:,2), ZIBM(:,:,:,3) ,& + TZHALO2_UT, TZHALO2_VT, TZHALO2_WT & +#ifndef MNH_OPENACC + ) +#else + , ZT3D(:,:,:,IZMEAN), ZT3D(:,:,:,IZWORK) ) +#endif + ENDIF + IF (HUVW_ADV_SCHEME=='CEN4TH') THEN + CALL ADVECUVW_2ND (ZUT, ZVT, ZWT, PRUCT, PRVCT, PRWCT, & + ZIBM(:,:,:,1), ZIBM(:,:,:,2), ZIBM(:,:,:,3)) + ENDIF + WHERE(XIBM_LS(:,:,:,2).GT.-XIBM_EPSI) ZRUS(:,:,:,JS)=ZIBM(:,:,:,1) + WHERE(XIBM_LS(:,:,:,3).GT.-XIBM_EPSI) ZRVS(:,:,:,JS)=ZIBM(:,:,:,2) + WHERE(XIBM_LS(:,:,:,4).GT.-XIBM_EPSI) ZRWS(:,:,:,JS)=ZIBM(:,:,:,3) + ZIBM(:,:,:,:)=1. + ENDIF ! - write ( ynum, '( I3 )' ) js + write ( ynum, '( I3 )' ) js #ifndef MNH_OPENACC - NULLIFY(TZFIELDS4_ll) - CALL ADD3DFIELD_ll( TZFIELDS4_ll, ZRUS(:,:,:,JS), 'ADVECUVW_RK::ZRUS(:,:,:,'//trim( adjustl( ynum ) )//')' ) - CALL ADD3DFIELD_ll( TZFIELDS4_ll, ZRVS(:,:,:,JS), 'ADVECUVW_RK::ZRVS(:,:,:,'//trim( adjustl( ynum ) )//')' ) - CALL ADD3DFIELD_ll( TZFIELDS4_ll, ZRWS(:,:,:,JS), 'ADVECUVW_RK::ZRWS(:,:,:,'//trim( adjustl( ynum ) )//')' ) - CALL UPDATE_HALO_ll(TZFIELDS4_ll,IINFO_ll) - CALL CLEANLIST_ll(TZFIELDS4_ll) + NULLIFY(TZFIELDS4_ll) + CALL ADD3DFIELD_ll( TZFIELDS4_ll, ZRUS(:,:,:,JS), 'ADVECUVW_RK::ZRUS(:,:,:,'//trim( adjustl( ynum ) )//')' ) + CALL ADD3DFIELD_ll( TZFIELDS4_ll, ZRVS(:,:,:,JS), 'ADVECUVW_RK::ZRVS(:,:,:,'//trim( adjustl( ynum ) )//')' ) + CALL ADD3DFIELD_ll( TZFIELDS4_ll, ZRWS(:,:,:,JS), 'ADVECUVW_RK::ZRWS(:,:,:,'//trim( adjustl( ynum ) )//')' ) + CALL UPDATE_HALO_ll(TZFIELDS4_ll,IINFO_ll) + CALL CLEANLIST_ll(TZFIELDS4_ll) #else - CALL GET_HALO_D(ZRUS(:,:,:,JS),HNAME='ADVECUVW_RK::ZRUS(:,:,:,'//trim( adjustl( ynum ) )//')' ) - CALL GET_HALO_D(ZRVS(:,:,:,JS),HNAME='ADVECUVW_RK::ZRVS(:,:,:,'//trim( adjustl( ynum ) )//')' ) - CALL GET_HALO_D(ZRWS(:,:,:,JS),HNAME='ADVECUVW_RK::ZRWS(:,:,:,'//trim( adjustl( ynum ) )//')' ) + CALL GET_HALO_D(ZRUS(:,:,:,JS),HNAME='ADVECUVW_RK::ZRUS(:,:,:,'//trim( adjustl( ynum ) )//')' ) + CALL GET_HALO_D(ZRVS(:,:,:,JS),HNAME='ADVECUVW_RK::ZRVS(:,:,:,'//trim( adjustl( ynum ) )//')' ) + CALL GET_HALO_D(ZRWS(:,:,:,JS),HNAME='ADVECUVW_RK::ZRWS(:,:,:,'//trim( adjustl( ynum ) )//')' ) ! acc update device(ZRUS(:,:,:,JS),ZRVS(:,:,:,JS),ZRWS(:,:,:,JS)) -#endif +#endif ! + IF (LIBM .AND. CIBM_ADV=='FREEZE') THEN + WHERE(XIBM_LS(:,:,:,2).GT.-XIBM_EPSI) ZRUS(:,:,:,JS)=ZUT(:,:,:)*PMXM_RHODJ(:,:,:)/PTSTEP + WHERE(XIBM_LS(:,:,:,3).GT.-XIBM_EPSI) ZRVS(:,:,:,JS)=ZVT(:,:,:)*PMYM_RHODJ(:,:,:)/PTSTEP + WHERE(XIBM_LS(:,:,:,4).GT.-XIBM_EPSI) ZRWS(:,:,:,JS)=ZWT(:,:,:)*PMZM_RHODJ(:,:,:)/PTSTEP + ENDIF + + IF (LIBM .AND. CIBM_ADV=='FORCIN') THEN + CALL SECOND_MNH(ZTIME1) + CALL IBM_FORCING_ADV(ZRUS(:,:,:,JS),ZRVS(:,:,:,JS),ZRWS(:,:,:,JS)) + CALL SECOND_MNH(ZTIME2) + XT_IBM_FORC = XT_IBM_FORC + ZTIME2 - ZTIME1 + ENDIF ! ! Guesses at the end of the RK loop ! -!$acc kernels present(PRUS_ADV,PRVS_ADV,PRWS_ADV,ZBUTS) present(ZRUS,ZRVS,ZRWS) - PRUS_ADV(:,:,:) = PRUS_ADV(:,:,:) + ZBUTS(JS) * ZRUS(:,:,:,JS) - PRVS_ADV(:,:,:) = PRVS_ADV(:,:,:) + ZBUTS(JS) * ZRVS(:,:,:,JS) - PRWS_ADV(:,:,:) = PRWS_ADV(:,:,:) + ZBUTS(JS) * ZRWS(:,:,:,JS) +!$acc kernels present(PRUS_ADV,PRVS_ADV,PRWS_ADV,ZBUTS) present(ZRUS,ZRVS,ZRWS,ZIBM) + IF ( .NOT. LIBM ) THEN + PRUS_ADV(:,:,:) = PRUS_ADV(:,:,:) + ZBUTS(JS) * ZRUS(:,:,:,JS) + PRVS_ADV(:,:,:) = PRVS_ADV(:,:,:) + ZBUTS(JS) * ZRVS(:,:,:,JS) + PRWS_ADV(:,:,:) = PRWS_ADV(:,:,:) + ZBUTS(JS) * ZRWS(:,:,:,JS) + ELSE + PRUS_ADV(:,:,:) = PRUS_ADV(:,:,:) + ZBUTS(JI) * ZRUS(:,:,:,JI) * ZIBM(:,:,:,1) + PRVS_ADV(:,:,:) = PRVS_ADV(:,:,:) + ZBUTS(JI) * ZRVS(:,:,:,JI) * ZIBM(:,:,:,2) + PRWS_ADV(:,:,:) = PRWS_ADV(:,:,:) + ZBUTS(JI) * ZRWS(:,:,:,JI) * ZIBM(:,:,:,3) + END IF !$acc end kernels ! IF ( JS < ISPL ) THEN @@ -531,20 +602,29 @@ END IF !$acc & present(ZRUS,ZRVS,ZRWS) present(PRUS_OTHER,PRVS_OTHER,PRWS_OTHER) & !$acc & present(PMXM_RHODJ,PMYM_RHODJ,PMZM_RHODJ) ! - ZUT = PU - ZVT = PV - ZWT = PW + ZUT(:,:,:) = PU(:,:,:) + ZVT(:,:,:) = PV(:,:,:) + ZWT(:,:,:) = PW(:,:,:) ! DO JI = 1, JS ! ! Intermediate guesses inside the RK loop ! - ZUT(:,:,:) = ZUT(:,:,:) + ZBUT(JS,JI) * PTSTEP * & - ( ZRUS(:,:,:,JI) + PRUS_OTHER(:,:,:) ) / PMXM_RHODJ - ZVT(:,:,:) = ZVT(:,:,:) + ZBUT(JS,JI) * PTSTEP * & - ( ZRVS(:,:,:,JI) + PRVS_OTHER(:,:,:) ) / PMYM_RHODJ - ZWT(:,:,:) = ZWT(:,:,:) + ZBUT(JS,JI) * PTSTEP * & - ( ZRWS(:,:,:,JI) + PRWS_OTHER(:,:,:) ) / PMZM_RHODJ + IF ( .NOT. LIBM ) THEN + ZUT(:,:,:) = ZUT(:,:,:) + ZBUT(JS,JI) * PTSTEP * & + ( ZRUS(:,:,:,JI) + PRUS_OTHER(:,:,:) ) / PMXM_RHODJ(:,:,:) + ZVT(:,:,:) = ZVT(:,:,:) + ZBUT(JS,JI) * PTSTEP * & + ( ZRVS(:,:,:,JI) + PRVS_OTHER(:,:,:) ) / PMYM_RHODJ(:,:,:) + ZWT(:,:,:) = ZWT(:,:,:) + ZBUT(JS,JI) * PTSTEP * & + ( ZRWS(:,:,:,JI) + PRWS_OTHER(:,:,:) ) / PMZM_RHODJ(:,:,:) + ELSE + ZUT(:,:,:) = ZUT(:,:,:) + ZBUT(JS,JI) * PTSTEP * & + ( ZRUS(:,:,:,JI) + PRUS_OTHER(:,:,:) ) / PMXM_RHODJ(:,:,:) * ZIBM(:,:,:,1) + ZVT(:,:,:) = ZVT(:,:,:) + ZBUT(JS,JI) * PTSTEP * & + ( ZRVS(:,:,:,JI) + PRVS_OTHER(:,:,:) ) / PMYM_RHODJ(:,:,:) * ZIBM(:,:,:,2) + ZWT(:,:,:) = ZWT(:,:,:) + ZBUT(JS,JI) * PTSTEP * & + ( ZRWS(:,:,:,JI) + PRWS_OTHER(:,:,:) ) / PMZM_RHODJ(:,:,:) * ZIBM(:,:,:,3) + END IF ! END DO !$acc end kernels @@ -552,7 +632,7 @@ END IF END IF ! ! End of the RK loop - END DO + END DO RKLOOP ! ! #ifdef MNH_OPENACC diff --git a/src/ZSOLVER/advecuvw_weno_k.f90 b/src/ZSOLVER/advecuvw_weno_k.f90 index f41009ae4ef63792834dff1b2f025e0808f1a2e2..fff7b03cddebe3b8ce237ee2a3c8430e20625559 100644 --- a/src/ZSOLVER/advecuvw_weno_k.f90 +++ b/src/ZSOLVER/advecuvw_weno_k.f90 @@ -73,30 +73,27 @@ END MODULE MODI_ADVECUVW_WENO_K !* 0. DECLARATIONS ! ------------ ! -USE MODE_ll -! -USE MODD_PARAMETERS -USE MODD_CONF USE MODD_ARGSLIST_ll, ONLY : HALO2LIST_ll -! -#ifndef MNH_OPENACC -USE MODI_SHUMAN -#else -USE MODI_SHUMAN_DEVICE +USE MODD_CONF, ONLY : L2D, NHALO +USE MODD_PARAMETERS + +USE MODE_ll +#ifdef MNH_OPENACC +USE MODE_DEVICE +USE MODE_MNH_ZWORK, ONLY : ZT3D, MNH_GET_ZT3D , MNH_REL_ZT3D #endif +USE MODE_MPPDB + USE MODI_ADVEC_WENO_K_1_AUX USE MODI_ADVEC_WENO_K_2_AUX USE MODI_ADVEC_WENO_K_3_AUX -! -USE MODD_CONF, ONLY : NHALO -USE MODE_MPPDB USE MODI_GET_HALO -! -#ifdef MNH_OPENACC -USE MODE_DEVICE -USE MODE_MNH_ZWORK, ONLY : ZT3D, MNH_GET_ZT3D , MNH_REL_ZT3D +#ifndef MNH_OPENACC +USE MODI_SHUMAN +#else +USE MODI_SHUMAN_DEVICE #endif -! + IMPLICIT NONE ! !* 0.1 Declarations of dummy arguments : @@ -210,7 +207,7 @@ CASE(1) ! WENO 1 CALL UP_UX_DEVICE(PUT,ZWORK,ZMEAN) CALL DXM_DEVICE(ZMEAN,ZWORK) !$acc kernels - PRUS = PRUS - ZWORK + PRUS(:,:,:) = PRUS(:,:,:) - ZWORK(:,:,:) !$acc end kernels ! !PRUS = PRUS - DYF(UP_MY(PUT,MXM(PRVCT))) @@ -218,7 +215,7 @@ CASE(1) ! WENO 1 CALL UP_MY_DEVICE(PUT,ZWORK,ZMEAN) CALL DYF_DEVICE(ZMEAN,ZWORK) !$acc kernels - PRUS = PRUS - ZWORK + PRUS(:,:,:) = PRUS(:,:,:) - ZWORK(:,:,:) !$acc end kernels ! !PRUS = PRUS - DZF(1,IKU,1,UP_MZ(PUT,MXM(PRWCT))) @@ -226,7 +223,7 @@ CASE(1) ! WENO 1 CALL UP_MZ_DEVICE(PUT,ZWORK,ZMEAN) CALL DZF_DEVICE(1,IKU,1,ZMEAN,ZWORK) !$acc kernels - PRUS = PRUS - ZWORK + PRUS(:,:,:) = PRUS(:,:,:) - ZWORK(:,:,:) !$acc end kernels ! ! V component @@ -236,7 +233,7 @@ CASE(1) ! WENO 1 CALL UP_MX_DEVICE(PVT,ZWORK,ZMEAN) CALL DXF_DEVICE(ZMEAN,ZWORK) !$acc kernels - PRVS = PRVS - ZWORK + PRVS(:,:,:) = PRVS(:,:,:) - ZWORK(:,:,:) !$acc end kernels ! !PRVS = PRVS - DYM(UP_VY(PVT,MYF(PRVCT))) @@ -244,7 +241,7 @@ CASE(1) ! WENO 1 CALL UP_VY_DEVICE(PVT,ZWORK,ZMEAN) CALL DYM_DEVICE(ZMEAN,ZWORK) !$acc kernels - PRVS = PRVS - ZWORK + PRVS(:,:,:) = PRVS(:,:,:) - ZWORK(:,:,:) !$acc end kernels ! !PRVS = PRVS - DZF(1,IKU,1,UP_MZ(PVT,MYM(PRWCT))) @@ -252,7 +249,7 @@ CASE(1) ! WENO 1 CALL UP_MZ_DEVICE(PVT,ZWORK,ZMEAN) CALL DZF_DEVICE(1,IKU,1,ZMEAN,ZWORK) !$acc kernels - PRVS = PRVS - ZWORK + PRVS(:,:,:) = PRVS(:,:,:) - ZWORK(:,:,:) !$acc end kernels ! ! W component @@ -262,7 +259,7 @@ CASE(1) ! WENO 1 CALL UP_MX_DEVICE(PWT,ZWORK,ZMEAN) CALL DXF_DEVICE(ZMEAN,ZWORK) !$acc kernels - PRWS = PRWS - ZWORK + PRWS(:,:,:) = PRWS(:,:,:) - ZWORK(:,:,:) !$acc end kernels ! !PRWS = PRWS - DYF(UP_MY(PWT,MZM(1,IKU,1,PRVCT))) @@ -270,7 +267,7 @@ CASE(1) ! WENO 1 CALL UP_MY_DEVICE(PWT,ZWORK,ZMEAN) CALL DYF_DEVICE(ZMEAN,ZWORK) !$acc kernels - PRWS = PRWS - ZWORK + PRWS(:,:,:) = PRWS(:,:,:) - ZWORK(:,:,:) !$acc end kernels ! !PRWS = PRWS - DZM(1,IKU,1,UP_WZ(PWT,MZF(1,IKU,1,PRWCT))) @@ -278,7 +275,7 @@ CASE(1) ! WENO 1 CALL UP_WZ_DEVICE(PWT,ZWORK,ZMEAN) CALL DZM_DEVICE(1,IKU,1,ZMEAN,ZWORK) !$acc kernels - PRWS = PRWS - ZWORK + PRWS(:,:,:) = PRWS(:,:,:) - ZWORK(:,:,:) !$acc end kernels #endif ! @@ -340,7 +337,7 @@ CASE(3) ! WENO 3 ZT3D(:,:,:,IZOMP1), ZT3D(:,:,:,IZOMP2), ZT3D(:,:,:,IZOMN1), ZT3D(:,:,:,IZOMN2) ) CALL DXM_DEVICE(ZMEAN,ZWORK) !$acc kernels - PRUS = PRUS - ZWORK + PRUS(:,:,:) = PRUS(:,:,:) - ZWORK(:,:,:) !$acc end kernels ! IF (.NOT.L2D) THEN @@ -351,7 +348,7 @@ CASE(3) ! WENO 3 ZT3D(:,:,:,IZOMP1), ZT3D(:,:,:,IZOMP2), ZT3D(:,:,:,IZOMN1), ZT3D(:,:,:,IZOMN2) ) CALL DYF_DEVICE(ZMEAN,ZWORK) !$acc kernels - PRUS = PRUS - ZWORK + PRUS(:,:,:) = PRUS(:,:,:) - ZWORK(:,:,:) !$acc end kernels END IF ! @@ -363,7 +360,7 @@ CASE(3) ! WENO 3 ZT3D(:,:,:,IZOMP1), ZT3D(:,:,:,IZOMP2), ZT3D(:,:,:,IZOMN1), ZT3D(:,:,:,IZOMN2) ) CALL DZF_DEVICE(1,IKU,1,ZMEAN,ZWORK) !$acc kernels - PRUS = PRUS - ZWORK + PRUS(:,:,:) = PRUS(:,:,:) - ZWORK(:,:,:) !$acc end kernels ! ! V component @@ -376,7 +373,7 @@ CASE(3) ! WENO 3 ZT3D(:,:,:,IZOMP1), ZT3D(:,:,:,IZOMP2), ZT3D(:,:,:,IZOMN1), ZT3D(:,:,:,IZOMN2) ) CALL DXF_DEVICE(ZMEAN,ZWORK) !$acc kernels - PRVS = PRVS - ZWORK + PRVS(:,:,:) = PRVS(:,:,:) - ZWORK(:,:,:) !$acc end kernels ! CALL MYF_DEVICE(PRVCT,ZWORK) @@ -386,7 +383,7 @@ CASE(3) ! WENO 3 ZT3D(:,:,:,IZOMP1), ZT3D(:,:,:,IZOMP2), ZT3D(:,:,:,IZOMN1), ZT3D(:,:,:,IZOMN2) ) CALL DYM_DEVICE(ZMEAN,ZWORK) !$acc kernels - PRVS = PRVS - ZWORK + PRVS(:,:,:) = PRVS(:,:,:) - ZWORK(:,:,:) !$acc end kernels ! ! PRVS = PRVS - DZF(1,IKU,1,WENO_K_2_MZ(PVT, MYM(PRWCT))) @@ -397,7 +394,7 @@ CASE(3) ! WENO 3 ZT3D(:,:,:,IZOMP1), ZT3D(:,:,:,IZOMP2), ZT3D(:,:,:,IZOMN1), ZT3D(:,:,:,IZOMN2) ) CALL DZF_DEVICE(1,IKU,1,ZMEAN,ZWORK) !$acc kernels - PRVS = PRVS - ZWORK + PRVS(:,:,:) = PRVS(:,:,:) - ZWORK(:,:,:) !$acc end kernels END IF ! @@ -411,7 +408,7 @@ CASE(3) ! WENO 3 ZT3D(:,:,:,IZOMP1), ZT3D(:,:,:,IZOMP2), ZT3D(:,:,:,IZOMN1), ZT3D(:,:,:,IZOMN2) ) CALL DXF_DEVICE(ZMEAN,ZWORK) !$acc kernels - PRWS = PRWS - ZWORK + PRWS(:,:,:) = PRWS(:,:,:) - ZWORK(:,:,:) !$acc end kernels ! IF (.NOT.L2D) THEN @@ -423,7 +420,7 @@ CASE(3) ! WENO 3 ZT3D(:,:,:,IZOMP1), ZT3D(:,:,:,IZOMP2), ZT3D(:,:,:,IZOMN1), ZT3D(:,:,:,IZOMN2) ) CALL DYF_DEVICE(ZMEAN,ZWORK) !$acc kernels - PRWS = PRWS - ZWORK + PRWS(:,:,:) = PRWS(:,:,:) - ZWORK(:,:,:) !$acc end kernels END IF ! @@ -435,7 +432,7 @@ CASE(3) ! WENO 3 ZT3D(:,:,:,IZOMP1), ZT3D(:,:,:,IZOMP2), ZT3D(:,:,:,IZOMN1), ZT3D(:,:,:,IZOMN2) ) CALL DZM_DEVICE(1,IKU,1,ZMEAN,ZWORK) !$acc kernels - PRWS = PRWS - ZWORK + PRWS(:,:,:) = PRWS(:,:,:) - ZWORK(:,:,:) !$acc end kernels ! CALL MNH_REL_ZT3D(IZFPOS1,IZFPOS2,IZFNEG1,IZFNEG2,IZBPOS1,IZBPOS2,IZBNEG1,IZBNEG2,IZOMP1,IZOMP2,IZOMN1,IZOMN2) @@ -517,9 +514,9 @@ CASE(5) ! WENO 5 CALL GET_HALO_D(ZMEAN)! Update HALO CALL DXM_DEVICE(ZMEAN,ZWORK) !$acc kernels - PRUS = PRUS - ZWORK + PRUS(:,:,:) = PRUS(:,:,:) - ZWORK(:,:,:) !$acc end kernels -! +! IF (.NOT.L2D) THEN CALL MXM_DEVICE(PRVCT,ZWORK) CALL ADVEC_WENO_K_3_MY(HLBCY, PUT, ZWORK, ZMEAN, & @@ -532,7 +529,7 @@ CASE(5) ! WENO 5 CALL GET_HALO_D(ZMEAN)! Update HALO CALL DYF_DEVICE(ZMEAN,ZWORK) !$acc kernels - PRUS = PRUS - ZWORK + PRUS(:,:,:) = PRUS(:,:,:) - ZWORK(:,:,:) !$acc end kernels END IF ! @@ -547,7 +544,7 @@ CASE(5) ! WENO 5 CALL GET_HALO_D(ZMEAN)! Update HALO - maybe not necessary (T.Lunet) CALL DZF_DEVICE(1,IKU,1,ZMEAN,ZWORK) !$acc kernels - PRUS = PRUS - ZWORK + PRUS(:,:,:) = PRUS(:,:,:) - ZWORK(:,:,:) !$acc end kernels ! ! V component @@ -564,7 +561,7 @@ CASE(5) ! WENO 5 CALL GET_HALO_D(ZMEAN)! Update HALO CALL DXF_DEVICE(ZMEAN,ZWORK) !$acc kernels - PRVS = PRVS - ZWORK + PRVS(:,:,:) = PRVS(:,:,:) - ZWORK(:,:,:) !$acc end kernels ! CALL MYF_DEVICE(PRVCT,ZWORK) @@ -578,7 +575,7 @@ CASE(5) ! WENO 5 CALL GET_HALO_D(ZMEAN)! Update HALO CALL DYM_DEVICE(ZMEAN,ZWORK) !$acc kernels - PRVS = PRVS - ZWORK + PRVS(:,:,:) = PRVS(:,:,:) - ZWORK(:,:,:) !$acc end kernels ! CALL MYM_DEVICE(PRWCT,ZWORK) @@ -592,7 +589,7 @@ CASE(5) ! WENO 5 CALL GET_HALO_D(ZMEAN)! Update HALO - maybe not necessary (T.Lunet) CALL DZF_DEVICE(1,IKU,1,ZMEAN,ZWORK) !$acc kernels - PRVS = PRVS - ZWORK + PRVS(:,:,:) = PRVS(:,:,:) - ZWORK(:,:,:) !$acc end kernels END IF ! @@ -609,7 +606,7 @@ CASE(5) ! WENO 5 CALL GET_HALO_D(ZMEAN)! Update HALO CALL DXF_DEVICE(ZMEAN,ZWORK) !$acc kernels - PRWS = PRWS - ZWORK + PRWS(:,:,:) = PRWS(:,:,:) - ZWORK(:,:,:) !$acc end kernels ! IF (.NOT.L2D) THEN @@ -624,7 +621,7 @@ CASE(5) ! WENO 5 CALL GET_HALO_D(ZMEAN)! Update HALO CALL DYF_DEVICE(ZMEAN,ZWORK) !$acc kernels - PRWS = PRWS - ZWORK + PRWS(:,:,:) = PRWS(:,:,:) - ZWORK(:,:,:) !$acc end kernels END IF ! @@ -639,7 +636,7 @@ CASE(5) ! WENO 5 CALL GET_HALO_D(ZMEAN)! Update HALO - maybe not necessary (T.Lunet) CALL DZM_DEVICE(1,IKU,1,ZMEAN,ZWORK) !$acc kernels - PRWS = PRWS - ZWORK + PRWS(:,:,:) = PRWS(:,:,:) - ZWORK(:,:,:) !$acc end kernels ! CALL MNH_REL_ZT3D(IZFPOS1,IZFPOS2,IZFPOS3,IZFNEG1,IZFNEG2,IZFNEG3,IZBPOS1, & diff --git a/src/ZSOLVER/anel_balancen.f90 b/src/ZSOLVER/anel_balancen.f90 index c1659356c3616b260ddda416eb328e07ddb3faca..27c6a691d385fb74d3c90db4da5031e2201fda85 100644 --- a/src/ZSOLVER/anel_balancen.f90 +++ b/src/ZSOLVER/anel_balancen.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. @@ -169,8 +169,6 @@ INTEGER :: IRRI ! Number of solid water variables REAL :: ZDRYMASST ! Mass of dry air Md REAL :: ZREFMASS ! Total mass of the ref. atmosphere REAL :: ZMASS_O_PHI0 ! Mass / Phi0 -LOGICAL :: GCLOSE_OUT ! switch for the LFI writing -CHARACTER (LEN= 28) :: YFMFILE ! virtual FM file INTEGER :: IMI ! model index !JUAN INTEGER :: IIU_B,IJU_B,IKU @@ -287,8 +285,6 @@ ZMASS_O_PHI0 = 1. ! | which is here not needed | IRR = 0 ! | | IRRL = 0 ! | | IRRI = 0 ! ============================================== -GCLOSE_OUT=.FALSE. -YFMFILE='UNUSED' ! IMI = GET_CURRENT_MODEL_INDEX() CALL PRESSUREZ(CLBCX,CLBCY,CPRESOPT,NITR,LITRADJ,ITCOUNT,XRELAX,IMI, & diff --git a/src/ZSOLVER/contrav.f90 b/src/ZSOLVER/contrav.f90 index 97df82096f4d57928bc287dce6eaf863ea163fa5..e14bb869b8d22be39a1db8f232069003c9a16489 100644 --- a/src/ZSOLVER/contrav.f90 +++ b/src/ZSOLVER/contrav.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. @@ -691,7 +691,7 @@ ELSE ! ------------------------------------ ! ! Z1(:,:,:) = 0. -! Z2(:,:,:) = 0. +! Z2(:,:,:) = 0. ! IF (KADV_ORDER == 2 ) THEN #ifdef MNH_OPENACC diff --git a/src/ZSOLVER/gdiv.f90 b/src/ZSOLVER/gdiv.f90 index 5bfd218bf519a7bf8bbbcae521aeffcf7e33ff0e..55c94571a3feaaa9c307900c45f3625e5b1c5456 100644 --- a/src/ZSOLVER/gdiv.f90 +++ b/src/ZSOLVER/gdiv.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. @@ -150,9 +150,9 @@ INTEGER :: IKE ! indice K for the last inner mass point along z ! INTEGER :: JI,JJ,JK ! loop indexes ! -#ifdef MNH_OPENACC INTEGER :: IIU,IJU,IKU ! +#ifdef MNH_OPENACC REAL, DIMENSION(:,:,:) , POINTER , CONTIGUOUS :: ZTMP1,ZTMP2 INTEGER :: IZTMP1,IZTMP2 #endif diff --git a/src/ZSOLVER/get_halo.f90 b/src/ZSOLVER/get_halo.f90 index 7dec4456cba761f8600a66e3bdf78a86b3c860de..e9950df6a4c015088369b1a146405741d140dc2a 100644 --- a/src/ZSOLVER/get_halo.f90 +++ b/src/ZSOLVER/get_halo.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. @@ -173,7 +173,7 @@ USE MODD_ARGSLIST_ll, ONLY : LIST_ll, HALO2LIST_ll ! IMPLICIT NONE ! -REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PSRC ! variable at t +REAL, DIMENSION(:,:,:), INTENT(IN) :: PSRC ! variable at t TYPE(HALO2LIST_ll), POINTER :: TP_PSRC_HALO2_ll ! halo2 for SRC character(len=*), optional, intent(in) :: HNAME ! Name of the field to be added ! @@ -215,7 +215,7 @@ USE MODD_ARGSLIST_ll, ONLY : LIST_ll ! IMPLICIT NONE ! -REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PSRC ! variable at t +REAL, DIMENSION(:,:,:), INTENT(IN) :: PSRC ! variable at t CHARACTER(len=4), OPTIONAL :: HDIR ! to send only halo on X or Y direction character(len=*), optional, intent(in) :: HNAME ! Name of the field to be added ! diff --git a/src/ZSOLVER/ini_modeln.f90 b/src/ZSOLVER/ini_modeln.f90 index 7fa50a868f56b4f03f92f506a153227281417929..4b63fd6b58eb24b958392ab3310afa91f59c7973 100644 --- a/src/ZSOLVER/ini_modeln.f90 +++ b/src/ZSOLVER/ini_modeln.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. @@ -287,6 +287,12 @@ END MODULE MODI_INI_MODEL_n ! P. Wautelet 10/04/2019: replace ABORT and STOP calls by Print_msg ! P. Wautelet 19/04/2019: removed unused dummy arguments and variables ! P. Wautelet 07/06/2019: allocate lookup tables for optical properties only when needed +! P. Wautelet 13/09/2019: budget: simplify and modernize date/time management +! C. Lac 11/2019: correction in the drag formula and application to building in addition to tree +! S. Riette 04/2020: XHL* fields +! F. Auguste 02/2021: add IBM +! T.Nigel 02/2021: add turbulence recycling +! J.L.Redelsperger 06/2011: OCEAN case !--------------------------------------------------------------------------------- ! !* 0. DECLARATIONS @@ -326,7 +332,7 @@ USE MODD_DEF_EDDYUV_FLUX_n ! FOR UV USE MODD_DIAG_FLAG, only: LCHEMDIAG, CSPEC_BU_DIAG USE MODD_DIM_n USE MODD_DRAG_n -USE MODD_DRAGTREE +USE MODD_DRAGTREE_n USE MODD_DUST use MODD_DUST_OPT_LKT, only: NMAX_RADIUS_LKT_DUST=>NMAX_RADIUS_LKT, NMAX_SIGMA_LKT_DUST=>NMAX_SIGMA_LKT, & NMAX_WVL_SW_DUST=>NMAX_WVL_SW, & @@ -337,6 +343,7 @@ USE MODD_DYN_n USE MODD_DYNZD USE MODD_DYNZD_n USE MODD_ELEC_n, only: XCION_POS_FW, XCION_NEG_FW +USE MODD_EOL_MAIN USE MODD_FIELD_n #ifdef MNH_FOREFIRE USE MODD_FOREFIRE @@ -347,10 +354,12 @@ USE MODD_FRC_n USE MODD_GET_n USE MODD_GRID_n USE MODD_GRID, only: XLONORI,XLATORI +USE MODD_IBM_PARAM_n, only: LIBM, XIBM_IEPS, XIBM_LS, XIBM_XMUT USE MODD_IO, only: CIO_DIR, TFILEDATA, TFILE_DUMMY USE MODD_IO_SURF_MNH, only: IO_SURF_MNH_MODEL USE MODD_LATZ_EDFLX USE MODD_LBC_n, only: CLBCX, CLBCY +use modd_les USE MODD_LSFIELD_n USE MODD_LUNIT_n USE MODD_MEAN_FIELD @@ -361,6 +370,7 @@ USE MODD_NESTING, only: CDAD_NAME, NDAD, NDT_2_WAY, NDTRATIO, NDXRATIO USE MODD_NSV USE MODD_NSV USE MODD_NUDGING_n, only: LNUDGING +USE MODD_OCEANH USE MODD_OUT_n USE MODD_PARAMETERS USE MODD_PARAM_KAFR_n @@ -372,6 +382,7 @@ USE MODD_PASPOL_n USE MODD_PAST_FIELD_n use modd_precision, only: LFIINT USE MODD_RADIATIONS_n +USE MODD_RECYCL_PARAM_n USE MODD_REF USE MODD_REF_n USE MODD_RELFRC_n @@ -390,6 +401,7 @@ USE MODD_TURB_n USE MODD_VAR_ll, only: IP USE MODE_GATHER_ll +use mode_ini_budget, only: Budget_preallocate, Ini_budget USE MODE_INI_ONE_WAY_n USE MODE_IO USE MODE_IO_FIELD_READ, only: IO_Field_read @@ -418,12 +430,13 @@ USE MODI_INI_AEROSET6 USE MODI_INI_AIRCRAFT_BALLOON USE MODI_INI_AIRCRAFT_BALLOON USE MODI_INI_BIKHARDT_n -USE MODI_INI_BUDGET USE MODI_INI_CPL USE MODI_INI_DEEP_CONVECTION USE MODI_INI_DRAG USE MODI_INI_DYNAMICS USE MODI_INI_ELEC_n +USE MODI_INI_EOL_ADNR +USE MODI_INI_EOL_ALM USE MODI_INI_LES_N USE MODI_INI_LG USE MODI_INI_LW_SETUP @@ -457,6 +470,11 @@ USE MODI_SUNPOS_n USE MODI_SURF_SOLAR_GEOM USE MODI_UPDATE_METRICS USE MODI_UPDATE_NSV +#ifdef MNH_ECRAD +#if ( VER_ECRAD == 140 ) +USE YOERDI , ONLY :RCCO2 +#endif +#endif ! IMPLICIT NONE ! @@ -484,6 +502,7 @@ LOGICAL :: GINIDCONV ! logical switch for the deep convection ! initialization LOGICAL :: GINIRAD ! logical switch for the radiation ! initialization +logical :: gles ! Logical to determine if LES diagnostics are enabled ! ! TYPE(LIST_ll), POINTER :: TZINITHALO2D_ll ! pointer for the list of 2D fields @@ -509,11 +528,12 @@ REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZDIR_ALB ! direct albedo REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZSCA_ALB ! diffuse albedo REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZEMIS ! emissivity REAL, DIMENSION(:,:), ALLOCATABLE :: ZTSRAD ! surface temperature +REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZIBM_LS ! LevelSet IBM ! ! INTEGER, DIMENSION(:,:),ALLOCATABLE :: IINDEX ! indices of non-zero terms INTEGER, DIMENSION(:),ALLOCATABLE :: IIND -INTEGER :: JM +INTEGER :: JM, JT ! !------------------------------------------ ! Dummy pointers needed to correct an ifort Bug @@ -551,6 +571,12 @@ IF (CCLOUD == 'LIMA') THEN LHORELAX_SVC1R3=LHORELAX_SVLIMA END IF ! +! UPDATE CONSTANTS FOR OCEAN MODEL +IF (LOCEAN) THEN + XP00=XP00OCEAN + XTH00=XTH00OCEAN +END IF +! ! NULLIFY(TZINITHALO2D_ll) NULLIFY(TZINITHALO3D_ll) @@ -691,6 +717,41 @@ CALL UPDATE_NSV(KMI) ! !* 3. ALLOCATE MEMORY ! ----------------- +! * Module RECYCL +! +IF (LRECYCL) THEN +! + NR_COUNT = 0 +! + ALLOCATE(XUMEANW(IJU,IKU,INT(XNUMBELT))) ; XUMEANW = 0.0 + ALLOCATE(XVMEANW(IJU,IKU,INT(XNUMBELT))) ; XVMEANW = 0.0 + ALLOCATE(XWMEANW(IJU,IKU,INT(XNUMBELT))) ; XWMEANW = 0.0 + ALLOCATE(XUMEANN(IIU,IKU,INT(XNUMBELT))) ; XUMEANN = 0.0 + ALLOCATE(XVMEANN(IIU,IKU,INT(XNUMBELT))) ; XVMEANN = 0.0 + ALLOCATE(XWMEANN(IIU,IKU,INT(XNUMBELT))) ; XWMEANN = 0.0 + ALLOCATE(XUMEANE(IJU,IKU,INT(XNUMBELT))) ; XUMEANE = 0.0 + ALLOCATE(XVMEANE(IJU,IKU,INT(XNUMBELT))) ; XVMEANE = 0.0 + ALLOCATE(XWMEANE(IJU,IKU,INT(XNUMBELT))) ; XWMEANE = 0.0 + ALLOCATE(XUMEANS(IIU,IKU,INT(XNUMBELT))) ; XUMEANS = 0.0 + ALLOCATE(XVMEANS(IIU,IKU,INT(XNUMBELT))) ; XVMEANS = 0.0 + ALLOCATE(XWMEANS(IIU,IKU,INT(XNUMBELT))) ; XWMEANS = 0.0 + ALLOCATE(XTBV(IIU,IJU,IKU)) ; XTBV = 0.0 +ELSE + ALLOCATE(XUMEANW(0,0,0)) + ALLOCATE(XVMEANW(0,0,0)) + ALLOCATE(XWMEANW(0,0,0)) + ALLOCATE(XUMEANN(0,0,0)) + ALLOCATE(XVMEANN(0,0,0)) + ALLOCATE(XWMEANN(0,0,0)) + ALLOCATE(XUMEANE(0,0,0)) + ALLOCATE(XVMEANE(0,0,0)) + ALLOCATE(XWMEANE(0,0,0)) + ALLOCATE(XUMEANS(0,0,0)) + ALLOCATE(XVMEANS(0,0,0)) + ALLOCATE(XWMEANS(0,0,0)) + ALLOCATE(XTBV (0,0,0)) +END IF +! ! !* 3.1 Module MODD_FIELD_n ! @@ -703,6 +764,7 @@ IF (LMEAN_FIELD) THEN ALLOCATE(XWM_MEAN(IIU,IJU,IKU)) ; XWM_MEAN = 0.0 ALLOCATE(XTHM_MEAN(IIU,IJU,IKU)) ; XTHM_MEAN = 0.0 ALLOCATE(XTEMPM_MEAN(IIU,IJU,IKU)) ; XTEMPM_MEAN = 0.0 + ALLOCATE(XSVT_MEAN(IIU,IJU,IKU)) ; XSVT_MEAN = 0.0 IF (CTURB/='NONE') THEN ALLOCATE(XTKEM_MEAN(IIU,IJU,IKU)) XTKEM_MEAN = 0.0 @@ -714,6 +776,7 @@ IF (LMEAN_FIELD) THEN ALLOCATE(XU2_MEAN(IIU,IJU,IKU)) ; XU2_MEAN = 0.0 ALLOCATE(XV2_MEAN(IIU,IJU,IKU)) ; XV2_MEAN = 0.0 ALLOCATE(XW2_MEAN(IIU,IJU,IKU)) ; XW2_MEAN = 0.0 + ALLOCATE(XUW_MEAN(IIU,IJU,IKU)) ; XUW_MEAN = 0.0 ALLOCATE(XTH2_MEAN(IIU,IJU,IKU)) ; XTH2_MEAN = 0.0 ALLOCATE(XTEMP2_MEAN(IIU,IJU,IKU)) ; XTEMP2_MEAN = 0.0 ALLOCATE(XPABS2_MEAN(IIU,IJU,IKU)) ; XPABS2_MEAN = 0.0 @@ -736,12 +799,14 @@ ELSE ALLOCATE(XWM_MEAN(0,0,0)) ALLOCATE(XTHM_MEAN(0,0,0)) ALLOCATE(XTEMPM_MEAN(0,0,0)) + ALLOCATE(XSVT_MEAN(0,0,0)) ALLOCATE(XTKEM_MEAN(0,0,0)) ALLOCATE(XPABSM_MEAN(0,0,0)) ! ALLOCATE(XU2_MEAN(0,0,0)) ALLOCATE(XV2_MEAN(0,0,0)) ALLOCATE(XW2_MEAN(0,0,0)) + ALLOCATE(XUW_MEAN(0,0,0)) ALLOCATE(XTH2_MEAN(0,0,0)) ALLOCATE(XTEMP2_MEAN(0,0,0)) ALLOCATE(XPABS2_MEAN(0,0,0)) @@ -793,6 +858,43 @@ ALLOCATE(XRWS_PRES(IIU,IJU,IKU)); XRWS_PRES = 0.0 ALLOCATE(XRTHS(IIU,IJU,IKU)) ; XRTHS = 0.0 !$acc enter data copyin(XRTHS) ALLOCATE(XRTHS_CLD(IIU,IJU,IKU)); XRTHS_CLD = 0.0 + +IF ( LIBM ) THEN + ALLOCATE(ZIBM_LS(IIU,IJU,IKU)) ; ZIBM_LS = 0.0 + ALLOCATE(XIBM_XMUT(IIU,IJU,IKU)); XIBM_XMUT = 0.0 +ELSE + ALLOCATE(ZIBM_LS (0,0,0)) + ALLOCATE(XIBM_XMUT(0,0,0)) +END IF + +IF ( LRECYCL ) THEN + ALLOCATE(XFLUCTUNW(IJU,IKU)) ; XFLUCTUNW = 0.0 + ALLOCATE(XFLUCTVNN(IIU,IKU)) ; XFLUCTVNN = 0.0 + ALLOCATE(XFLUCTUTN(IIU,IKU)) ; XFLUCTUTN = 0.0 + ALLOCATE(XFLUCTVTW(IJU,IKU)) ; XFLUCTVTW = 0.0 + ALLOCATE(XFLUCTUNE(IJU,IKU)) ; XFLUCTUNE = 0.0 + ALLOCATE(XFLUCTVNS(IIU,IKU)) ; XFLUCTVNS = 0.0 + ALLOCATE(XFLUCTUTS(IIU,IKU)) ; XFLUCTUTS = 0.0 + ALLOCATE(XFLUCTVTE(IJU,IKU)) ; XFLUCTVTE = 0.0 + ALLOCATE(XFLUCTWTW(IJU,IKU)) ; XFLUCTWTW = 0.0 + ALLOCATE(XFLUCTWTN(IIU,IKU)) ; XFLUCTWTN = 0.0 + ALLOCATE(XFLUCTWTE(IJU,IKU)) ; XFLUCTWTE = 0.0 + ALLOCATE(XFLUCTWTS(IIU,IKU)) ; XFLUCTWTS = 0.0 +ELSE + ALLOCATE(XFLUCTUNW(0,0)) + ALLOCATE(XFLUCTVNN(0,0)) + ALLOCATE(XFLUCTUTN(0,0)) + ALLOCATE(XFLUCTVTW(0,0)) + ALLOCATE(XFLUCTUNE(0,0)) + ALLOCATE(XFLUCTVNS(0,0)) + ALLOCATE(XFLUCTUTS(0,0)) + ALLOCATE(XFLUCTVTE(0,0)) + ALLOCATE(XFLUCTWTW(0,0)) + ALLOCATE(XFLUCTWTN(0,0)) + ALLOCATE(XFLUCTWTE(0,0)) + ALLOCATE(XFLUCTWTS(0,0)) +END IF +! IF (CTURB /= 'NONE') THEN ALLOCATE(XTKET(IIU,IJU,IKU)) ALLOCATE(XRTKES(IIU,IJU,IKU)) @@ -841,10 +943,25 @@ ELSE ALLOCATE(XSRCT(0,0,0)) ALLOCATE(XSIGS(0,0,0)) END IF +IF (CCLOUD == 'ICE3'.OR.CCLOUD == 'ICE4') THEN + ALLOCATE(XHLC_HRC(IIU,IJU,IKU)) + ALLOCATE(XHLC_HCF(IIU,IJU,IKU)) + ALLOCATE(XHLI_HRI(IIU,IJU,IKU)) + ALLOCATE(XHLI_HCF(IIU,IJU,IKU)) + XHLC_HRC(:,:,:)=0. + XHLC_HCF(:,:,:)=0. + XHLI_HRI(:,:,:)=0. + XHLI_HCF(:,:,:)=0. +ELSE + ALLOCATE(XHLC_HRC(0,0,0)) + ALLOCATE(XHLC_HCF(0,0,0)) + ALLOCATE(XHLI_HRI(0,0,0)) + ALLOCATE(XHLI_HCF(0,0,0)) +END IF ! IF (NRR>1) THEN - ALLOCATE(XCLDFR(IIU,IJU,IKU)) - ALLOCATE(XRAINFR(IIU,IJU,IKU)) + ALLOCATE(XCLDFR(IIU,IJU,IKU)); XCLDFR (:, :, :) = 0. + ALLOCATE(XRAINFR(IIU,IJU,IKU)); XRAINFR(:, :, :) = 0. ELSE ALLOCATE(XCLDFR(0,0,0)) ALLOCATE(XRAINFR(0,0,0)) @@ -906,11 +1023,19 @@ ALLOCATE(XDZZ(IIU,IJU,IKU)) ! !* 3.3 Modules MODD_REF and MODD_REF_n ! -IF (KMI == 1) THEN +! Different reference states for Ocean and Atmosphere models +! For the moment, same reference states for O and A +!IF ((KMI == 1).OR.LCOUPLES) THEN +IF (KMI==1) THEN ALLOCATE(XRHODREFZ(IKU),XTHVREFZ(IKU)) +ELSE IF (LCOUPLES) THEN +! in coupled O-A case, need different variables for ocean + ALLOCATE(XRHODREFZO(IKU),XTHVREFZO(IKU)) ELSE !Do not allocate XRHODREFZ and XTHVREFZ because they are the same on all grids (not 'n' variables) END IF +! +ALLOCATE(XPHIT(IIU,IJU,IKU)) ALLOCATE(XRHODREF(IIU,IJU,IKU)) ALLOCATE(XTHVREF(IIU,IJU,IKU)) ALLOCATE(XEXNREF(IIU,IJU,IKU)) @@ -1339,15 +1464,31 @@ END IF ! Initialization of SW bands NSWB_OLD = 6 ! Number of bands in ECMWF original scheme (from Fouquart et Bonnel (1980)) ! then modified through INI_RADIATIONS_ECMWF but remains equal to 6 practically + +#ifdef MNH_ECRAD +#if ( VER_ECRAD == 140 ) +NLWB_OLD = 16 ! For XEMIS initialization (should be spectral in the future) +#endif +#endif + +NLWB_MNH = 16 ! For XEMIS initialization (should be spectral in the future) + IF (CRAD == 'ECRA') THEN NSWB_MNH = 14 +#ifdef MNH_ECRAD +#if ( VER_ECRAD == 140 ) + NLWB_MNH = 16 +#endif +#endif ELSE NSWB_MNH = NSWB_OLD +#ifdef MNH_ECRAD +#if ( VER_ECRAD == 140 ) + NLWB_MNH = NLWB_OLD +#endif +#endif END IF -NLWB_MNH = 16 ! For XEMIS initialization (should be spectral in the future) - - ALLOCATE(XSW_BANDS (NSWB_MNH)) ALLOCATE(XLW_BANDS (NLWB_MNH)) ALLOCATE(XZENITH (IIU,IJU)) @@ -1498,7 +1639,17 @@ ALLOCATE(XRI_MF(IIU,IJU,IKU)) ; XRI_MF=0.0 ! ALLOCATE(ZJ(IIU,IJU,IKU)) ! -!* 3.10 Forcing variables (Module MODD_FRC) +!* 3.10 Forcing variables (Module MODD_FRC and MODD_FRCn) +! +IF ( LFORCING ) THEN + ALLOCATE(XWTFRC(IIU,IJU,IKU)) ; XWTFRC = XUNDEF + ALLOCATE(XUFRC_PAST(IIU,IJU,IKU)) ; XUFRC_PAST = XUNDEF + ALLOCATE(XVFRC_PAST(IIU,IJU,IKU)) ; XVFRC_PAST = XUNDEF +ELSE + ALLOCATE(XWTFRC(0,0,0)) + ALLOCATE(XUFRC_PAST(0,0,0)) + ALLOCATE(XVFRC_PAST(0,0,0)) +END IF ! IF (KMI == 1) THEN IF ( LFORCING ) THEN @@ -1530,15 +1681,6 @@ IF (KMI == 1) THEN ALLOCATE(XTENDUFRC(0,0)) ALLOCATE(XTENDVFRC(0,0)) END IF - IF ( LFORCING ) THEN - ALLOCATE(XWTFRC(IIU,IJU,IKU)) - ALLOCATE(XUFRC_PAST(IIU,IJU,IKU)) ; XUFRC_PAST = XUNDEF - ALLOCATE(XVFRC_PAST(IIU,IJU,IKU)) ; XVFRC_PAST = XUNDEF - ELSE - ALLOCATE(XWTFRC(0,0,0)) - ALLOCATE(XUFRC_PAST(0,0,0)) - ALLOCATE(XVFRC_PAST(0,0,0)) - END IF ELSE !Do not allocate because they are the same on all grids (not 'n' variables) END IF @@ -1667,12 +1809,21 @@ ENDIF !* 4. INITIALIZE BUDGET VARIABLES ! --------------------------- ! +gles = lles_mean .or. lles_resolved .or. lles_subgrid .or. lles_updraft & + .or. lles_downdraft .or. lles_spectra +!Called if budgets are enabled via NAM_BUDGET +!or if LES budgets are enabled via NAM_LES (condition on kmi==1 to call it max once) +if ( ( cbutype /= "NONE" .and. nbumod == kmi ) .or. ( gles .and. kmi == 1 ) ) THEN + call Budget_preallocate() +end if + IF ( CBUTYPE /= "NONE" .AND. NBUMOD == KMI ) THEN - CALL INI_BUDGET(ILUOUT,XTSTEP,NSV,NRR, & + CALL Ini_budget(ILUOUT,XTSTEP,NSV,NRR, & LNUMDIFU,LNUMDIFTH,LNUMDIFSV, & LHORELAX_UVWTH,LHORELAX_RV, LHORELAX_RC,LHORELAX_RR, & LHORELAX_RI,LHORELAX_RS,LHORELAX_RG, LHORELAX_RH,LHORELAX_TKE, & - LHORELAX_SV,LVE_RELAX,LCHTRANS,LNUDGING,LDRAGTREE,LDEPOTREE, & + LHORELAX_SV, LVE_RELAX, LVE_RELAX_GRD, & + LCHTRANS,LNUDGING,LDRAGTREE,LDEPOTREE,LMAIN_EOL, & CRAD,CDCONV,CSCONV,CTURB,CTURBDIM,CCLOUD ) END IF ! @@ -1710,13 +1861,6 @@ IF (KMI == 1) THEN IF ( NDAD(KMI) == 7) CDAD_NAME(KMI) = CEXP//'.7.'//CSEG IF ( NDAD(KMI) == 8) CDAD_NAME(KMI) = CEXP//'.8.'//CSEG END IF -ELSE - ALLOCATE(XUM(0,0,0)) - ALLOCATE(XVM(0,0,0)) - ALLOCATE(XWM(0,0,0)) - ALLOCATE(XDUM(0,0,0)) - ALLOCATE(XDVM(0,0,0)) - ALLOCATE(XDWM(0,0,0)) END IF ! !------------------------------------------------------------------------------- @@ -1756,8 +1900,8 @@ NDT_2_WAY(KMI)=4 ! IF ( LUSECHEM .OR. LCHEMDIAG ) THEN IF ((KMI==1).AND.(CPROGRAM == "MESONH".OR.CPROGRAM == "DIAG ")) & - CALL CH_INIT_JVALUES(TDTCUR%TDATE%DAY, TDTCUR%TDATE%MONTH, & - TDTCUR%TDATE%YEAR, ILUOUT, XCH_TUV_DOBNEW) + CALL CH_INIT_JVALUES(TDTCUR%nday, TDTCUR%nmonth, & + TDTCUR%nyear, ILUOUT, XCH_TUV_DOBNEW) ! IF (LORILAM) THEN CALL CH_AER_MOD_INIT @@ -1774,11 +1918,12 @@ IF (CCLOUD=='LIMA') CALL INIT_AEROSOL_PROPERTIES ! -------------------------------- ! CALL MPPDB_CHECK3D(XUT,"INI_MODEL_N-before read_field::XUT",PRECISION) -CALL READ_FIELD(TPINIFILE,IIU,IJU,IKU, & +CALL READ_FIELD(KMI,TPINIFILE,IIU,IJU,IKU, & CGETTKET,CGETRVT,CGETRCT,CGETRRT,CGETRIT,CGETCIT,CGETZWS, & CGETRST,CGETRGT,CGETRHT,CGETSVT,CGETSRCT,CGETSIGS,CGETCLDFR, & - CGETBL_DEPTH,CGETSBL_DEPTH,CGETPHC,CGETPHR,CUVW_ADV_SCHEME, & - CTEMP_SCHEME,NSIZELBX_ll,NSIZELBXU_ll,NSIZELBY_ll,NSIZELBYV_ll,& + CGETBL_DEPTH,CGETSBL_DEPTH,CGETPHC,CGETPHR, & + CUVW_ADV_SCHEME, CTEMP_SCHEME, & + NSIZELBX_ll, NSIZELBXU_ll, NSIZELBY_ll, NSIZELBYV_ll, & NSIZELBXTKE_ll,NSIZELBYTKE_ll, & NSIZELBXR_ll,NSIZELBYR_ll,NSIZELBXSV_ll,NSIZELBYSV_ll, & XUM,XVM,XWM,XDUM,XDVM,XDWM, & @@ -1797,7 +1942,10 @@ CALL READ_FIELD(TPINIFILE,IIU,IJU,IKU, & NADVFRC,TDTADVFRC,XDTHFRC,XDRVFRC, & NRELFRC,TDTRELFRC,XTHREL,XRVREL, & XVTH_FLUX_M,XWTH_FLUX_M,XVU_FLUX_M, & - XRUS_PRES,XRVS_PRES,XRWS_PRES,XRTHS_CLD,XRRS_CLD,XRSVS_CLD ) + XRUS_PRES,XRVS_PRES,XRWS_PRES,XRTHS_CLD,XRRS_CLD,XRSVS_CLD, & + ZIBM_LS,XIBM_XMUT,XUMEANW,XVMEANW,XWMEANW,XUMEANN,XVMEANN, & + XWMEANN,XUMEANE,XVMEANE,XWMEANE,XUMEANS,XVMEANS,XWMEANS ) + ! !------------------------------------------------------------------------------- ! @@ -1949,7 +2097,6 @@ IF ((KMI==1).AND.(.NOT. LSTEADYLS)) THEN ENDDO ! END IF -ALLOCATE(XLSZWSM(IIU,IJU)) ; XLSZWSM = -1. ! IF ( KMI > 1) THEN ! Use dummy pointers to correct an ifort BUG @@ -2093,6 +2240,13 @@ CALL INI_DYNAMICS(XLON,XLAT,XRHODJ,XTHVREF,XMAP,XZZ,XDXHAT,XDYHAT, & IF (LDRAG) THEN CALL INI_DRAG(LMOUNT,XZS,XHSTART,NSTART,XDRAG) ENDIF +!* 16.2 Initialize the LevelSet function +! ------------- +IF (LIBM) THEN + ALLOCATE(XIBM_LS(IIU,IJU,IKU,4)) ; XIBM_LS = -XIBM_IEPS + XIBM_LS(:,:,:,1)=ZIBM_LS(:,:,:) + DEALLOCATE(ZIBM_LS) +ENDIF !------------------------------------------------------------------------------- ! !* 17. SURFACE FIELDS @@ -2167,10 +2321,14 @@ ALLOCATE(ZSCA_ALB(IIU,IJU,NSWB_MNH)) ALLOCATE(ZEMIS (IIU,IJU,NLWB_MNH)) ALLOCATE(ZTSRAD (IIU,IJU)) ! -IF ((TPINIFILE%NMNHVERSION(1)==4 .AND. TPINIFILE%NMNHVERSION(2)>=6) .OR. TPINIFILE%NMNHVERSION(1)>4) THEN - CALL IO_Field_read(TPINIFILE,'SURF',CSURF) +IF (LCOUPLES.AND.(KMI>1))THEN + CSURF ="NONE" ELSE - CSURF = "EXTE" + IF ((TPINIFILE%NMNHVERSION(1)==4 .AND. TPINIFILE%NMNHVERSION(2)>=6) .OR. TPINIFILE%NMNHVERSION(1)>4) THEN + CALL IO_Field_read(TPINIFILE,'SURF',CSURF) + ELSE + CSURF = "EXTE" + END IF END IF ! ! @@ -2217,7 +2375,6 @@ ELSE XTSRAD = XTT XSEA = 1. END IF - END IF IF (CSURF=='EXTE' .AND. (CPROGRAM=='SPAWN ')) THEN ! ouverture du fichier PGD @@ -2279,14 +2436,14 @@ IF (CRAD == 'ECMW') THEN ZBARE(:,:) = 0. END IF ! - IF ( CAOP=='EXPL' .AND. LDUST ) THEN + IF ( CAOP=='EXPL' .AND. LDUST .AND. KMI==1) THEN ALLOCATE( XEXT_COEFF_WVL_LKT_DUST( NMAX_RADIUS_LKT_DUST, NMAX_SIGMA_LKT_DUST, NMAX_WVL_SW_DUST ) ) ALLOCATE( XEXT_COEFF_550_LKT_DUST( NMAX_RADIUS_LKT_DUST, NMAX_SIGMA_LKT_DUST ) ) ALLOCATE( XPIZA_LKT_DUST ( NMAX_RADIUS_LKT_DUST, NMAX_SIGMA_LKT_DUST, NMAX_WVL_SW_DUST ) ) ALLOCATE( XCGA_LKT_DUST ( NMAX_RADIUS_LKT_DUST, NMAX_SIGMA_LKT_DUST, NMAX_WVL_SW_DUST ) ) END IF ! - IF ( CAOP=='EXPL' .AND. LSALT ) THEN + IF ( CAOP=='EXPL' .AND. LSALT .AND. KMI==1) THEN ALLOCATE( XEXT_COEFF_WVL_LKT_SALT( NMAX_RADIUS_LKT_SALT, NMAX_SIGMA_LKT_SALT, NMAX_WVL_SW_SALT ) ) ALLOCATE( XEXT_COEFF_550_LKT_SALT( NMAX_RADIUS_LKT_SALT, NMAX_SIGMA_LKT_SALT ) ) ALLOCATE( XPIZA_LKT_SALT ( NMAX_RADIUS_LKT_SALT, NMAX_SIGMA_LKT_SALT, NMAX_WVL_SW_SALT ) ) @@ -2424,18 +2581,18 @@ CALL INI_AIRCRAFT_BALLOON(TPINIFILE,XTSTEP, TDTSEG, XSEGLEN, NRR, NSV, & !* 24. STATION initializations ! ----------------------- ! -CALL INI_SURFSTATION_n(XTSTEP, TDTSEG, XSEGLEN, NRR, NSV, & - CTURB=="TKEL" , & - XLATORI, XLONORI ) +CALL INI_SURFSTATION_n(XTSTEP, XSEGLEN, NRR, NSV, & + CTURB=="TKEL" , KMI, & + XLATORI, XLONORI ) ! !------------------------------------------------------------------------------- ! !* 25. PROFILER initializations ! ------------------------ ! -CALL INI_POSPROFILER_n(XTSTEP, TDTSEG, XSEGLEN, NRR, NSV, & - CTURB=="TKEL", & - XLATORI, XLONORI ) +CALL INI_POSPROFILER_n(XTSTEP, XSEGLEN, NRR, NSV, & + CTURB=="TKEL", & + XLATORI, XLONORI ) ! !------------------------------------------------------------------------------- ! @@ -2473,7 +2630,7 @@ ENDIF !------------------------------------ IF ( LFOREFIRE ) THEN CALL INIT_FOREFIRE_n(KMI, ILUOUT, IP & - , TDTCUR%TDATE%YEAR, TDTCUR%TDATE%MONTH, TDTCUR%TDATE%DAY, TDTCUR%TIME, XTSTEP) + , TDTCUR%nyear, TDTCUR%nmonth, TDTCUR%nday, TDTCUR%xtime, XTSTEP) END IF #endif @@ -2521,6 +2678,38 @@ IF (LCHEMDIAG) THEN ELSE ALLOCATE(XTCHEM(0)) END IF - +!------------------------------------------------------------------------------- +! +!* 32. Wind turbine +! +IF (LMAIN_EOL .AND. KMI == NMODEL_EOL) THEN + ALLOCATE(XFX_RG(IIU,IJU,IKU)) + ALLOCATE(XFY_RG(IIU,IJU,IKU)) + ALLOCATE(XFZ_RG(IIU,IJU,IKU)) + ALLOCATE(XFX_SMR_RG(IIU,IJU,IKU)) + ALLOCATE(XFY_SMR_RG(IIU,IJU,IKU)) + ALLOCATE(XFZ_SMR_RG(IIU,IJU,IKU)) + SELECT CASE(CMETH_EOL) + CASE('ADNR') + CALL INI_EOL_ADNR + CASE('ALM') + CALL INI_EOL_ALM(XDXX,XDYY) + END SELECT +END IF +! +!* 33. Auto-coupling Atmos-Ocean LES NH +! +IF (LCOUPLES) THEN + ALLOCATE(XSSUFL_C(IIU,IJU,1)); XSSUFL_C=0.0 + ALLOCATE(XSSVFL_C(IIU,IJU,1)); XSSVFL_C=0.0 + ALLOCATE(XSSTFL_C(IIU,IJU,1)); XSSTFL_C=0.0 + ALLOCATE(XSSRFL_C(IIU,IJU,1)); XSSRFL_C=0. +ELSE + ALLOCATE(XSSUFL_C(0,0,0)) + ALLOCATE(XSSVFL_C(0,0,0)) + ALLOCATE(XSSTFL_C(0,0,0)) + ALLOCATE(XSSRFL_C(0,0,0)) +END IF +! END SUBROUTINE INI_MODEL_n diff --git a/src/ZSOLVER/ini_spectren.f90 b/src/ZSOLVER/ini_spectren.f90 index 4783443207ecf52aa183b97da09c2f4014f94b58..2385b3f6d1a573768bb008d1d811a600de531292 100644 --- a/src/ZSOLVER/ini_spectren.f90 +++ b/src/ZSOLVER/ini_spectren.f90 @@ -58,7 +58,7 @@ USE MODD_CTURB USE MODD_CURVCOR_n USE MODD_DEEP_CONVECTION_n USE MODD_DIM_n -USE MODD_DRAGTREE +USE MODD_DRAGTREE_n USE MODD_DUST USE MODD_DYN USE MODD_DYN_n diff --git a/src/ZSOLVER/modd_dynn.f90 b/src/ZSOLVER/modd_dynn.f90 index 445105291501c6f96285c0eb1654124d801036a5..4ac45311f08f8fe70a71f607168f27529ee1f0ee 100644 --- a/src/ZSOLVER/modd_dynn.f90 +++ b/src/ZSOLVER/modd_dynn.f90 @@ -1,6 +1,6 @@ -!MNH_LIC Copyright 1994-2018 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 version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt !MNH_LIC for details. version 1. !----------------------------------------------------------------- ! ################# @@ -43,6 +43,7 @@ !! Modification 01/2016 (JP Pinty) Add LIMA !! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O !! Modification 07/2017 (V. Vionnet) Add blowing snow variable +!! Modification 03/2021 (JL Redelsperger) Add logical LOCEAN !------------------------------------------------------------------------------- ! !* 0. DECLARATIONS @@ -55,6 +56,7 @@ TYPE DYN_t ! INTEGER :: NSTOP ! Number of time step REAL :: XTSTEP ! Time step + LOGICAL :: LOCEAN ! !++++++++++++++++++++++++++++++++++ !PART USED BY THE PRESSURE SOLVER @@ -202,6 +204,7 @@ CHARACTER(LEN=5), POINTER :: CPRESOPT=>NULL() INTEGER, POINTER :: NITR=>NULL() LOGICAL, POINTER :: LITRADJ=>NULL() LOGICAL, POINTER :: LRES=>NULL() +LOGICAL, POINTER :: LOCEAN=>NULL() REAL, POINTER :: XRES=>NULL() REAL, POINTER :: XRELAX=>NULL() REAL, POINTER :: XDXHATM=>NULL() @@ -331,6 +334,7 @@ CPRESOPT=>DYN_MODEL(KTO)%CPRESOPT NITR=>DYN_MODEL(KTO)%NITR LITRADJ=>DYN_MODEL(KTO)%LITRADJ LRES=>DYN_MODEL(KTO)%LRES +LOCEAN=>DYN_MODEL(KTO)%LOCEAN XRES=>DYN_MODEL(KTO)%XRES XRELAX=>DYN_MODEL(KTO)%XRELAX XDXHATM=>DYN_MODEL(KTO)%XDXHATM diff --git a/src/ZSOLVER/modeln.f90 b/src/ZSOLVER/modeln.f90 index a65e4e58abd7369dc942b9d7956f844df1cfa6dc..f06c77e4c4af9ef274c078717b1bf9e8e7d3803c 100644 --- a/src/ZSOLVER/modeln.f90 +++ b/src/ZSOLVER/modeln.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. @@ -168,7 +168,7 @@ END MODULE MODI_MODEL_n !! July 29,1996 (Lafore) nesting introduction !! Aug. 1,1996 (Lafore) synchronization between models !! Sept. 4,1996 (Lafore) modification of call to routine SET_COUPLING -!! now splitted in 2 routines +!! now split in 2 routines !! (UVW_LS_COUPLING and SCALAR_LS_COUPLING) !! Sept 5,1996 (V.Masson) print of loop index for debugging !! purposes @@ -265,7 +265,14 @@ END MODULE MODI_MODEL_n ! J. Escobar 09/07/2019: norme Doctor -> Rename Module Type variable TZ -> T ! J. Escobar 09/07/2019: for bug in management of XLSZWSM variable, add/use specific 2D TLSFIELD2D_ll pointer ! P. Wautelet 23/07/2019: OpenACC: move data creations from resolved_cloud to modeln and optimize updates +! P. Wautelet 13/09/2019: budget: simplify and modernize date/time management ! J. Escobar 27/09/2019: add missing report timing of RESOLVED_ELEC +! P. Wautelet 02-03/2020: use the new data structures and subroutines for budgets +! P. Wautelet 12/10/2020: Write_les_n: remove HLES_AVG dummy argument and group all 4 calls +! F. Auguste 01/02/2021: add IBM +! T. Nagel 01/02/2021: add turbulence recycling +! P. Wautelet 19/02/2021: add NEGA2 term for SV budgets +! J.L. Redelsperger 03/2021: add Call NHOA_COUPLN (coupling O & A LES version) !!------------------------------------------------------------------------------- ! !* 0. DECLARATIONS @@ -274,34 +281,41 @@ END MODULE MODI_MODEL_n USE MODD_2D_FRC USE MODD_ADV_n USE MODD_AIRCRAFT_BALLOON +USE MODD_ARGSLIST_ll, ONLY : LIST_ll USE MODD_BAKOUT USE MODD_BIKHARDT_n -USE MODD_BLANK -USE MODD_BUDGET +USE MODD_BLANK_n +USE MODD_BLOWSNOW +USE MODD_BLOWSNOW_n +use modd_budget, only: cbutype, lbu_ru, lbu_rv, lbu_rw, lbudget_u, lbudget_v, lbudget_w, lbudget_sv, lbu_enable, & + NBUDGET_U, NBUDGET_V, NBUDGET_W, NBUDGET_SV1, nbumod, nbutime, & + tbudgets, tburhodj, & + xtime_bu, xtime_bu_process USE MODD_CH_AERO_n, ONLY: XSOLORG, XMI USE MODD_CH_MNHC_n, ONLY: LUSECHEM,LCH_CONV_LINOX,LUSECHAQ,LUSECHIC, & LCH_INIT_FIELD USE MODD_CLOUD_MF_n -USE MODD_VISCOSITY -USE MODD_DRAG_n USE MODD_CLOUDPAR_n USE MODD_CONF USE MODD_CONF_n USE MODD_CURVCOR_n USE MODD_DEEP_CONVECTION_n USE MODD_DIM_n +USE MODD_DRAG_n USE MODD_DUST, ONLY: LDUST USE MODD_DYN USE MODD_DYN_n USE MODD_DYNZD USE MODD_DYNZD_n USE MODD_ELEC_DESCR +USE MODD_EOL_MAIN USE MODD_FIELD_n USE MODD_FRC USE MODD_FRC_n USE MODD_GET_n USE MODD_GRID, ONLY: XLONORI,XLATORI USE MODD_GRID_n +USE MODD_IBM_PARAM_n, ONLY: CIBM_ADV, LIBM, LIBM_TROUBLE, XIBM_LS USE MODD_ICE_C1R3_DESCR, ONLY: XRTMIN_C1R3=>XRTMIN USE MODD_IO, ONLY: LIO_NO_WRITE, TFILEDATA, TFILE_SURFEX, TFILE_DUMMY USE MODD_LBC_n @@ -309,7 +323,7 @@ USE MODD_LES USE MODD_LES_BUDGET USE MODD_LIMA_PRECIP_SCAVENGING_n USE MODD_LSFIELD_n -USE MODD_LUNIT, ONLY: TLUOUT0,TOUTDATAFILE +USE MODD_LUNIT, ONLY: TOUTDATAFILE USE MODD_LUNIT_n, ONLY: TDIAFILE,TINIFILE,TINIFILEPGD,TLUOUT USE MODD_MEAN_FIELD USE MODD_MEAN_FIELD_n @@ -327,8 +341,6 @@ USE MODD_PARAM_LIMA, ONLY: MSEDC => LSEDC, MWARM => LWARM, MRAIN => LRAIN, & MACTIT => LACTIT, LSCAV, LCOLD, & MSEDI => LSEDI, MHHONI => LHHONI, LHAIL, & XRTMIN_LIMA=>XRTMIN, MACTTKE=>LACTTKE -USE MODD_BLOWSNOW_n -USE MODD_BLOWSNOW USE MODD_PARAM_MFSHALL_n USE MODD_PARAM_n USE MODD_PAST_FIELD_n @@ -337,6 +349,8 @@ use modd_precision, only: MNHTIME USE MODD_PROFILER_n USE MODD_RADIATIONS_n, ONLY: XTSRAD,XSCAFLASWD,XDIRFLASWD,XDIRSRFSWD, XAER, XDTHRAD USE MODD_RAIN_ICE_DESCR, ONLY: XRTMIN +USE MODD_RECYCL_PARAM_n, ONLY: LRECYCL +USE MODD_REF, ONLY: LCOUPLES USE MODD_REF_n USE MODD_SALT, ONLY: LSALT USE MODD_SERIES, ONLY: LSERIES @@ -348,7 +362,9 @@ USE MODD_TIME_n USE MODD_TIMEZ USE MODD_TURB_CLOUD, ONLY: NMODEL_CLOUD,CTURBLEN_CLOUD,XCEI USE MODD_TURB_n +USE MODD_VISCOSITY ! +use mode_budget, only: Budget_store_init, Budget_store_end USE MODE_DATETIME USE MODE_ELEC_ll USE MODE_GRIDCART @@ -357,11 +373,19 @@ USE MODE_IO_FIELD_WRITE, only: IO_Field_user_write, IO_Fieldlist_write, IO_Hea USE MODE_IO_FILE, only: IO_File_close, IO_File_open USE MODE_IO_MANAGE_STRUCT, only: IO_File_add2list USE MODE_ll +#ifdef MNH_IOLFI +use mode_menu_diachro, only: MENU_DIACHRO +#endif USE MODE_MNH_TIMING USE MODE_MODELN_HANDLER USE MODE_MPPDB +USE MODE_MSG USE MODE_ONE_WAY_n +use mode_write_les_n, only: Write_les_n +use mode_write_lfifmn_fordiachro_n, only: WRITE_LFIFMN_FORDIACHRO_n +USE MODE_WRITE_PROFILER_n, ONLY: WRITE_PROFILER_n ! +USE MODI_ADDFLUCTUATIONS USE MODI_ADVECTION_METSV USE MODI_ADVECTION_UVW USE MODI_ADVECTION_UVW_CEN @@ -384,18 +408,20 @@ USE MODI_FORC_SQUALL_LINE USE MODI_FORC_WIND USE MODI_GET_HALO USE MODI_GRAVITY_IMPL +USE MODI_IBM_INIT +USE MODI_IBM_FORCING +USE MODI_IBM_FORCING_TR +USE MODI_IBM_FORCING_ADV USE MODI_INI_DIAG_IN_RUN USE MODI_INI_LG USE MODI_INI_MEAN_FIELD USE MODI_INITIAL_GUESS USE MODI_LES_INI_TIMESTEP_n USE MODI_LES_N -USE MODI_VISCOSITY USE MODI_LIMA_PRECIP_SCAVENGING USE MODI_LS_COUPLING USE MODI_MASK_COMPRESS USE MODI_MEAN_FIELD -USE MODI_MENU_DIACHRO USE MODI_MNHGET_SURF_PARAM_n USE MODI_MNHWRITE_ZS_DUMMY_n USE MODI_NUDGING @@ -404,6 +430,7 @@ USE MODI_PHYS_PARAM_n USE MODI_PRESSUREZ USE MODI_PROFILER_n USE MODI_RAD_BOUND +USE MODI_RECYCLING USE MODI_RELAX2FW_ION USE MODI_RELAXATION USE MODI_REL_FORCING_n @@ -418,13 +445,11 @@ USE MODI_STATION_n USE MODI_TURB_CLOUD_INDEX USE MODI_TWO_WAY USE MODI_UPDATE_NSV +USE MODI_VISCOSITY USE MODI_WRITE_AIRCRAFT_BALLOON USE MODI_WRITE_DESFM_n USE MODI_WRITE_DIAG_SURF_ATM_N -USE MODI_WRITE_LES_n USE MODI_WRITE_LFIFM_n -USE MODI_WRITE_LFIFMN_FORDIACHRO_n -USE MODI_WRITE_PROFILER_n USE MODI_WRITE_SERIES_n USE MODI_WRITE_STATION_n USE MODI_WRITE_SURF_ATM_N @@ -612,10 +637,6 @@ IF (KTCOUNT == 1) THEN ! ALLOCATE(XWT_ACT_NUC(SIZE(XWT,1),SIZE(XWT,2),SIZE(XWT,3))) ALLOCATE(GMASKkids(SIZE(XWT,1),SIZE(XWT,2))) -! -! initialization of the FM file backup/output number - IBAK=0 - IOUT=0 ! IF ( .NOT. LIO_NO_WRITE ) THEN CALL IO_File_open(TDIAFILE) @@ -719,12 +740,15 @@ IF (KTCOUNT == 1) THEN XT_TURB = 0.0_MNHTIME XT_MAFL = 0.0_MNHTIME XT_DRAG = 0.0_MNHTIME + XT_EOL = 0.0_MNHTIME XT_TRACER = 0.0_MNHTIME XT_SHADOWS = 0.0_MNHTIME XT_ELEC = 0.0_MNHTIME XT_CHEM = 0.0_MNHTIME XT_2WAY = 0.0_MNHTIME ! + XT_IBM_FORC = 0.0_MNHTIME + ! END IF ! !* 1.7 Allocation of arrays for observation diagnostics @@ -744,11 +768,14 @@ CALL SECOND_MNH2(ZTIME1) ! ISYNCHRO = MODULO (KTCOUNT, NDTRATIO(IMI) ) ! test of synchronisation ! - - -IF (IMI/=1 .AND. NDAD(IMI)/=IMI .AND. (ISYNCHRO==1 .OR. NDTRATIO(IMI) == 1) ) THEN -! - ! Use dummy pointers to correct an ifort BUG +! +IF (LCOUPLES.AND.LOCEAN) THEN + CALL NHOA_COUPL_n(NDAD(IMI),XTSTEP,IMI,KTCOUNT,IKU) +END IF +! No Gridnest in coupled OA LES for now +IF (.NOT. LCOUPLES .AND. IMI/=1 .AND. NDAD(IMI)/=IMI .AND. (ISYNCHRO==1 .OR. NDTRATIO(IMI) == 1) ) THEN +! +! Use dummy pointers to correct an ifort BUG DPTR_XBMX1=>XBMX1 DPTR_XBMX2=>XBMX2 DPTR_XBMX3=>XBMX3 @@ -875,7 +902,28 @@ IF (IMI/=1 .AND. NDAD(IMI)/=IMI .AND. (ISYNCHRO==1 .OR. NDTRATIO(IMI) == 1) ) TH END IF ! CALL SECOND_MNH2(ZTIME2) -XT_1WAY = XT_1WAY + ZTIME2 - ZTIME1 +XT_1WAY = XT_1WAY + ZTIME2 - ZTIME1 +! +!* 2.1 RECYCLING TURBULENCE +! ---- +IF (CTURB /= 'NONE' .AND. LRECYCL) THEN + CALL RECYCLING(XFLUCTUNW,XFLUCTVNN,XFLUCTUTN,XFLUCTVTW,XFLUCTWTW,XFLUCTWTN, & + XFLUCTUNE,XFLUCTVNS,XFLUCTUTS,XFLUCTVTE,XFLUCTWTE,XFLUCTWTS, & + KTCOUNT) +ENDIF +! +!* 2.2 IBM +! ---- +! +IF (LIBM .AND. KTCOUNT==1) THEN + ! + IF (.NOT.LCARTESIAN) THEN + CALL PRINT_MSG(NVERB_FATAL, 'GEN', 'MODELN', 'IBM can only be used in combination with cartesian coordinates') + ENDIF + ! + CALL IBM_INIT(XIBM_LS) + ! +ENDIF ! !------------------------------------------------------------------------------- ! @@ -895,7 +943,7 @@ CALL BOUNDARIES ( & XLBYUM,XLBYVM,XLBYWM,XLBYTHM,XLBYTKEM,XLBYRM,XLBYSVM, & XLBXUS,XLBXVS,XLBXWS,XLBXTHS,XLBXTKES,XLBXRS,XLBXSVS, & XLBYUS,XLBYVS,XLBYWS,XLBYTHS,XLBYTKES,XLBYRS,XLBYSVS, & - XRHODJ, & + XRHODJ,XRHODREF, & XUT, XVT, XWT, XTHT, XTKET, XRT, XSVT, XSRCT ) END IF ! @@ -913,19 +961,18 @@ IF (CSURF=='EXTE') CALL GOTO_SURFEX(IMI) ! ZTIME1 = ZTIME2 ! -IF (IBAK < NBAK_NUMB ) THEN - IF (KTCOUNT == TBACKUPN(IBAK+1)%NSTEP) THEN - IBAK=IBAK+1 - GCLOSE_OUT=.TRUE. +IF ( nfile_backup_current < NBAK_NUMB ) THEN + IF ( KTCOUNT == TBACKUPN(nfile_backup_current + 1)%NSTEP ) THEN + nfile_backup_current = nfile_backup_current + 1 ! - TZBAKFILE => TBACKUPN(IBAK)%TFILE + TZBAKFILE => TBACKUPN(nfile_backup_current)%TFILE IVERB = TZBAKFILE%NLFIVERB ! CALL IO_File_open(TZBAKFILE) ! CALL WRITE_DESFM_n(IMI,TZBAKFILE) - CALL IO_Header_write(TBACKUPN(IBAK)%TFILE) - CALL WRITE_LFIFM_n(TBACKUPN(IBAK)%TFILE,TBACKUPN(IBAK)%TFILE%TDADFILE%CNAME) + CALL IO_Header_write( TBACKUPN(nfile_backup_current)%TFILE ) + CALL WRITE_LFIFM_n( TBACKUPN(nfile_backup_current)%TFILE, TBACKUPN(nfile_backup_current)%TFILE%TDADFILE%CNAME ) TOUTDATAFILE => TZBAKFILE CALL MNHWRITE_ZS_DUMMY_n(TZBAKFILE) IF (CSURF=='EXTE') THEN @@ -958,17 +1005,17 @@ ELSE TZBAKFILE => TFILE_DUMMY END IF ! -IF (IOUT < NOUT_NUMB ) THEN - IF (KTCOUNT == TOUTPUTN(IOUT+1)%NSTEP) THEN - IOUT=IOUT+1 +IF ( nfile_output_current < NOUT_NUMB ) THEN + IF ( KTCOUNT == TOUTPUTN(nfile_output_current + 1)%NSTEP ) THEN + nfile_output_current = nfile_output_current + 1 ! - TZOUTFILE => TOUTPUTN(IOUT)%TFILE + TZOUTFILE => TOUTPUTN(nfile_output_current)%TFILE ! CALL IO_File_open(TZOUTFILE) ! CALL IO_Header_write(TZOUTFILE) - CALL IO_Fieldlist_write(TOUTPUTN(IOUT)) - CALL IO_Field_user_write(TOUTPUTN(IOUT)) + CALL IO_Fieldlist_write( TOUTPUTN(nfile_output_current) ) + CALL IO_Field_user_write( TOUTPUTN(nfile_output_current) ) ! CALL IO_File_close(TZOUTFILE) ! @@ -981,6 +1028,42 @@ XT_STORE = XT_STORE + ZTIME2 - ZTIME1 ! !------------------------------------------------------------------------------- ! +!* 4.BIS IBM and Fluctuations application +! ----------------------------- +! +!* 4.B1 Add fluctuations at the domain boundaries +! +IF (LRECYCL) THEN + CALL ADDFLUCTUATIONS ( & + CLBCX,CLBCY, & + XUT, XVT, XWT, XTHT, XTKET, XRT, XSVT, XSRCT, & + XFLUCTUTN,XFLUCTVTW,XFLUCTUTS,XFLUCTVTE, & + XFLUCTWTW,XFLUCTWTN,XFLUCTWTS,XFLUCTWTE ) +ENDIF +! +!* 4.B2 Immersed boundaries +! +IF (LIBM) THEN + ! + ZTIME1=ZTIME2 + ! + IF (.NOT.LCARTESIAN) THEN + CALL PRINT_MSG(NVERB_FATAL, 'GEN', 'MODELN', 'IBM can only be used in combination with cartesian coordinates') + ENDIF + ! + CALL IBM_FORCING(XUT,XVT,XWT,XTHT,XRT,XSVT,XTKET) + ! + IF (LIBM_TROUBLE) THEN + CALL IBM_FORCING_TR(XUT,XVT,XWT,XTHT,XRT,XSVT,XTKET) + ENDIF + ! + CALL SECOND_MNH2(ZTIME2) + ! + XT_IBM_FORC = XT_IBM_FORC + ZTIME2 - ZTIME1 + ! +ENDIF +!------------------------------------------------------------------------------- +! !* 5. INITIALIZATION OF THE BUDGET VARIABLES ! -------------------------------------- ! @@ -991,28 +1074,34 @@ ELSE END IF ! IF (NBUMOD==IMI .AND. CBUTYPE=='MASK' ) THEN - CALL SET_MASK - IF (LBU_RU) XBURHODJU(:,NBUTIME,:) = XBURHODJU(:,NBUTIME,:) & - + MASK_COMPRESS(MXM(XRHODJ)) - IF (LBU_RV) XBURHODJV(:,NBUTIME,:) = XBURHODJV(:,NBUTIME,:) & - + MASK_COMPRESS(MYM(XRHODJ)) - IF (LBU_RW) XBURHODJW(:,NBUTIME,:) = XBURHODJW(:,NBUTIME,:) & - + MASK_COMPRESS(MZM(XRHODJ)) - IF (ALLOCATED(XBURHODJ)) & - XBURHODJ (:,NBUTIME,:) = XBURHODJ (:,NBUTIME,:) & - + MASK_COMPRESS(XRHODJ) + CALL SET_MASK() + if ( lbu_ru ) then + tbudgets(NBUDGET_U)%trhodj%xdata(:, nbutime, :) = tbudgets(NBUDGET_U)%trhodj%xdata(:, nbutime, :) & + + Mask_compress( Mxm( xrhodj(:, :, :) ) ) + end if + if ( lbu_rv ) then + tbudgets(NBUDGET_V)%trhodj%xdata(:, nbutime, :) = tbudgets(NBUDGET_V)%trhodj%xdata(:, nbutime, :) & + + Mask_compress( Mym( xrhodj(:, :, :) ) ) + end if + if ( lbu_rw ) then + tbudgets(NBUDGET_W)%trhodj%xdata(:, nbutime, :) = tbudgets(NBUDGET_W)%trhodj%xdata(:, nbutime, :) & + + Mask_compress( Mzm( xrhodj(:, :, :) ) ) + end if + if ( associated( tburhodj ) ) tburhodj%xdata(:, nbutime, :) = tburhodj%xdata(:, nbutime, :) + Mask_compress( xrhodj(:, :, :) ) END IF ! IF (NBUMOD==IMI .AND. CBUTYPE=='CART' ) THEN - IF (LBU_RU) XBURHODJU(:,:,:) = XBURHODJU(:,:,:) & - + CART_COMPRESS(MXM(XRHODJ)) - IF (LBU_RV) XBURHODJV(:,:,:) = XBURHODJV(:,:,:) & - + CART_COMPRESS(MYM(XRHODJ)) - IF (LBU_RW) XBURHODJW(:,:,:) = XBURHODJW(:,:,:) & - + CART_COMPRESS(MZM(XRHODJ)) - IF (ALLOCATED(XBURHODJ)) & - XBURHODJ (:,:,:) = XBURHODJ (:,:,:) & - + CART_COMPRESS(XRHODJ) + if ( lbu_ru ) then + tbudgets(NBUDGET_U)%trhodj%xdata(:, :, :) = tbudgets(NBUDGET_U)%trhodj%xdata(:, :, :) + Cart_compress( Mxm( xrhodj(:, :, :) ) ) + end if + if ( lbu_rv ) then + tbudgets(NBUDGET_V)%trhodj%xdata(:, :, :) = tbudgets(NBUDGET_V)%trhodj%xdata(:, :, :) + Cart_compress( Mym( xrhodj(:, :, :) ) ) + end if + if ( lbu_rw ) then + tbudgets(NBUDGET_W)%trhodj%xdata(:, :, :) = tbudgets(NBUDGET_W)%trhodj%xdata(:, :, :) & + + Cart_compress( Mzm( xrhodj(:, :, :) ) ) + end if + if ( associated( tburhodj ) ) tburhodj%xdata(:, :, :) = tburhodj%xdata(:, :, :) + Cart_compress( xrhodj(:, :, :) ) END IF ! CALL BUDGET_FLAGS(LUSERV, LUSERC, LUSERR, & @@ -1116,7 +1205,7 @@ END IF ! IF ( LFORCING ) THEN CALL FORCING(XTSTEP,LUSERV,XRHODJ,XCORIOZ,XZHAT,XZZ,TDTCUR,& - XUFRC_PAST, XVFRC_PAST, & + XUFRC_PAST, XVFRC_PAST,XWTFRC, & XUT,XVT,XWT,XTHT,XTKET,XRT,XSVT, & XRUS,XRVS,XRWS,XRTHS,XRTKES,XRRS,XRSVS,IMI,ZJ) END IF @@ -1212,7 +1301,13 @@ IF ( LNUMDIFU .OR. LNUMDIFTH .OR. LNUMDIFSV ) THEN LZDIFFU,LNUMDIFU, LNUMDIFTH, LNUMDIFSV, & THALO2T_ll, TLSHALO2_ll,XZDIFFU_HALO2 ) END IF -! + +if ( lbudget_sv ) then + do jsv = 1, nsv + call Budget_store_init( tbudgets(jsv + NBUDGET_SV1 - 1), 'NEGA2', xrsvs(:, :, :, jsv) ) + end do +end if + DO JSV = NSV_CHEMBEG,NSV_CHEMEND XRSVS(:,:,:,JSV) = MAX(XRSVS(:,:,:,JSV),0.) END DO @@ -1258,6 +1353,12 @@ IF (CELEC .NE. 'NONE') THEN XRSVS(:,:,:,NSV_ELECBEG) = MAX(XRSVS(:,:,:,NSV_ELECBEG),0.) XRSVS(:,:,:,NSV_ELECEND) = MAX(XRSVS(:,:,:,NSV_ELECEND),0.) END IF + +if ( lbudget_sv ) then + do jsv = 1, nsv + call Budget_store_end( tbudgets(jsv + NBUDGET_SV1 - 1), 'NEGA2', xrsvs(:, :, :, jsv) ) + end do +end if ! CALL SECOND_MNH2(ZTIME2) ! @@ -1319,10 +1420,10 @@ XT_RELAX = XT_RELAX + ZTIME2 - ZTIME1 & ! ZTIME1 = ZTIME2 ! -CALL PHYS_PARAM_n(KTCOUNT,TZBAKFILE, GCLOSE_OUT, & - XT_RAD,XT_SHADOWS,XT_DCONV,XT_GROUND,XT_MAFL, & - XT_DRAG,XT_TURB,XT_TRACER, & - ZTIME,ZWETDEPAER,GMASKkids,GCLOUD_ONLY) +CALL PHYS_PARAM_n( KTCOUNT, TZBAKFILE, & + XT_RAD, XT_SHADOWS, XT_DCONV, XT_GROUND, & + XT_MAFL, XT_DRAG, XT_EOL, XT_TURB, XT_TRACER, & + ZTIME, ZWETDEPAER, GMASKkids, GCLOUD_ONLY ) ! IF (CDCONV/='NONE') THEN XPACCONV = XPACCONV + XPRCONV * XTSTEP @@ -1332,8 +1433,8 @@ IF (CDCONV/='NONE') THEN END IF END IF ! -IF (IBAK>0 .AND. IBAK <= NBAK_NUMB ) THEN - IF (KTCOUNT == TBACKUPN(IBAK)%NSTEP) THEN +IF ( nfile_backup_current > 0 .AND. nfile_backup_current <= NBAK_NUMB ) THEN + IF ( KTCOUNT == TBACKUPN(nfile_backup_current)%NSTEP ) THEN IF (CSURF=='EXTE') THEN CALL GOTO_SURFEX(IMI) CALL DIAG_SURF_ATM_n(YSURF_CUR,'MESONH') @@ -1464,9 +1565,9 @@ XT_COUPL = XT_COUPL + ZTIME2 - ZTIME1 !* 8 Bis . Blowing snow scheme ! --------- ! -IF(LBLOWSNOW) THEN - CALL BLOWSNOW(CLBCX,CLBCY,XTSTEP,NRR,XPABST,XTHT,XRT,XZZ,XRHODREF, & - XRHODJ,XEXNREF,XRRS,XRTHS,XSVT,XRSVS,XSNWSUBL3D ) +IF ( LBLOWSNOW ) THEN + CALL BLOWSNOW( XTSTEP, NRR, XPABST, XTHT, XRT, XZZ, XRHODREF, & + XRHODJ, XEXNREF, XRRS, XRTHS, XSVT, XRSVS, XSNWSUBL3D ) ENDIF ! !----------------------------------------------------------------------- @@ -1516,7 +1617,7 @@ XTIME_LES_BU_PROCESS = 0. !$acc data copyin (XTKET, XRSVS_CLD) & !$acc & copy (XRTKES, XRSVS) & !$acc & copyout(XRTKEMS) -CALL ADVECTION_METSV ( TZBAKFILE, GCLOSE_OUT,CUVW_ADV_SCHEME, & +CALL ADVECTION_METSV ( TZBAKFILE, CUVW_ADV_SCHEME, & CMET_ADV_SCHEME, CSV_ADV_SCHEME, CCLOUD, NSPLIT, & LSPLIT_CFL, XSPLIT_CFL, LCFL_WRIT, & CLBCX, CLBCY, NRR, NSV, TDTCUR, XTSTEP, & @@ -1559,6 +1660,18 @@ CALL SECOND_MNH2(ZTIME2) ! XT_GRAV = XT_GRAV + ZTIME2 - ZTIME1 - XTIME_LES_BU_PROCESS - XTIME_BU_PROCESS ! +IF ( LIBM .AND. CIBM_ADV=='FORCIN' ) THEN + ! + ZTIME1=ZTIME2 + ! + CALL IBM_FORCING_ADV (XRUS,XRVS,XRWS) + ! + CALL SECOND_MNH2(ZTIME2) + ! + XT_IBM_FORC = XT_IBM_FORC + ZTIME2 - ZTIME1 + ! +ENDIF +! ZTIME1 = ZTIME2 XTIME_BU_PROCESS = 0. XTIME_LES_BU_PROCESS = 0. @@ -1614,10 +1727,10 @@ XT_ADVUVW = XT_ADVUVW + ZTIME2 - ZTIME1 - XTIME_LES_BU_PROCESS - XTIME_BU_PROCES !------------------------------------------------------------------------------- ! IF (NMODEL_CLOUD==IMI .AND. CTURBLEN_CLOUD/='NONE') THEN - CALL TURB_CLOUD_INDEX(XTSTEP,TZBAKFILE, & - LTURB_DIAG,GCLOSE_OUT,NRRI, & - XRRS,XRT,XRHODJ,XDXX,XDYY,XDZZ,XDZX,XDZY, & - XCEI ) + CALL TURB_CLOUD_INDEX( XTSTEP, TZBAKFILE, & + LTURB_DIAG, NRRI, & + XRRS, XRT, XRHODJ, XDXX, XDYY, XDZZ, XDZX, XDZY, & + XCEI ) END IF ! !------------------------------------------------------------------------------- @@ -1630,12 +1743,19 @@ ZTIME1 = ZTIME2 ZRUS=XRUS ZRVS=XRVS ZRWS=XRWS -! - CALL RAD_BOUND (CLBCX,CLBCY,CTURB,XCARPKMAX, & + +if ( .not. l1d ) then + if ( lbudget_u ) call Budget_store_init( tbudgets(NBUDGET_U), 'PRES', xrus(:, :, :) ) + if ( lbudget_v ) call Budget_store_init( tbudgets(NBUDGET_V), 'PRES', xrvs(:, :, :) ) + if ( lbudget_w ) call Budget_store_init( tbudgets(NBUDGET_W), 'PRES', xrws(:, :, :) ) +end if + +CALL RAD_BOUND (CLBCX,CLBCY,CTURB,XCARPKMAX, & XTSTEP, & XDXHAT, XDYHAT, XZHAT, & XUT, XVT, & XLBXUM, XLBYVM, XLBXUS, XLBYVS, & + XFLUCTUNW,XFLUCTVNN,XFLUCTUNE,XFLUCTVNS, & XCPHASE, XCPHASE_PBL, XRHODJ, & XTKET,XRUS, XRVS, XRWS ) ZRUS=XRUS-ZRUS @@ -1801,7 +1921,7 @@ IF (CCLOUD /= 'NONE' .AND. CELEC == 'NONE') THEN CALL RESOLVED_CLOUD ( CCLOUD, CACTCCN, CSCONV, CMF_CLOUD, NRR, NSPLITR, & NSPLITG, IMI, KTCOUNT, & CLBCX,CLBCY,TZBAKFILE, CRAD, CTURBDIM, & - GCLOSE_OUT, LSUBG_COND,LSIGMAS,CSUBG_AUCV,XTSTEP, & + LSUBG_COND,LSIGMAS,CSUBG_AUCV,XTSTEP, & XZZ, XRHODJ, XRHODREF, XEXNREF, & ZPABST, XTHT,XRT,XSIGS,VSIGQSAT,XMFCONV,XTHM,XRCM, & XPABSM, XWT_ACT_NUC,XDTHRAD, XRTHS, XRRS, & @@ -1809,13 +1929,14 @@ IF (CCLOUD /= 'NONE' .AND. CELEC == 'NONE') THEN XSRCT, XCLDFR,XCIT, & LSEDIC,KACTIT, KSEDC, KSEDI, KRAIN, KWARM, KHHONI, & LCONVHG, XCF_MF,XRC_MF, XRI_MF, & -! XINPRC,ZINPRC3D,XINPRR, XINPRR3D, XEVAP3D, & +! XINPRC,ZINPRC3D,XINPRR, XINPRR3D, XEVAP3D, & ! XINPRS,ZINPRS3D, XINPRG,ZINPRG3D, XINPRH,ZINPRH3D, & - XINPRC,XINPRR, XINPRR3D, XEVAP3D, & - XINPRS, XINPRG,XINPRH, & ! XSOLORG, XMI,ZSPEEDC, ZSPEEDR, ZSPEEDS, ZSPEEDG, ZSPEEDH, & + XINPRC,XINPRR, XINPRR3D, XEVAP3D, & + XINPRS, XINPRG, XINPRH, & XSOLORG, XMI, & XINDEP, XSUPSAT, XNACT, XNPRO,XSSPRO, XRAINFR, & + XHLC_HRC, XHLC_HCF, XHLI_HRI, XHLI_HCF, & ZSEA, ZTOWN ) !$acc end data DEALLOCATE(ZTOWN) @@ -1823,7 +1944,7 @@ IF (CCLOUD /= 'NONE' .AND. CELEC == 'NONE') THEN CALL RESOLVED_CLOUD ( CCLOUD, CACTCCN, CSCONV, CMF_CLOUD, NRR, NSPLITR, & NSPLITG, IMI, KTCOUNT, & CLBCX,CLBCY,TZBAKFILE, CRAD, CTURBDIM, & - GCLOSE_OUT, LSUBG_COND,LSIGMAS,CSUBG_AUCV, & + LSUBG_COND,LSIGMAS,CSUBG_AUCV, & XTSTEP,XZZ, XRHODJ, XRHODREF, XEXNREF, & ZPABST, XTHT,XRT,XSIGS,VSIGQSAT,XMFCONV,XTHM,XRCM, & XPABSM, XWT_ACT_NUC,XDTHRAD, XRTHS, XRRS, & @@ -1831,13 +1952,14 @@ IF (CCLOUD /= 'NONE' .AND. CELEC == 'NONE') THEN XSRCT, XCLDFR,XCIT, & LSEDIC,KACTIT, KSEDC, KSEDI, KRAIN, KWARM, KHHONI, & LCONVHG, XCF_MF,XRC_MF, XRI_MF, & -! XINPRC,ZINPRC3D,XINPRR, XINPRR3D, XEVAP3D, & +! XINPRC,ZINPRC3D,XINPRR, XINPRR3D, XEVAP3D, & ! XINPRS,ZINPRS3D, XINPRG,ZINPRG3D, XINPRH,ZINPRH3D, & - XINPRC,XINPRR, XINPRR3D, XEVAP3D, & - XINPRS, XINPRG,XINPRH, & ! XSOLORG, XMI,ZSPEEDC, ZSPEEDR, ZSPEEDS, ZSPEEDG, ZSPEEDH, & + XINPRC,XINPRR, XINPRR3D, XEVAP3D, & + XINPRS, XINPRG, XINPRH, & XSOLORG, XMI, & - XINDEP, XSUPSAT, XNACT, XNPRO,XSSPRO, XRAINFR ) + XINDEP, XSUPSAT, XNACT, XNPRO,XSSPRO, XRAINFR, & + XHLC_HRC, XHLC_HCF, XHLI_HRI, XHLI_HCF ) END IF !$acc end data @@ -1968,7 +2090,7 @@ XT_SPECTRA = XT_SPECTRA + ZTIME2 - ZTIME1 + XTIME_LES_BU + XTIME_LES ! -------------------- ! IF (LMEAN_FIELD) THEN - CALL MEAN_FIELD(XUT, XVT, XWT, XTHT, XTKET, XPABST) + CALL MEAN_FIELD(XUT, XVT, XWT, XTHT, XTKET, XPABST, XSVT(:,:,:,1)) END IF ! !------------------------------------------------------------------------------- @@ -2026,7 +2148,6 @@ ZTIME1 = ZTIME2 ! IF (LFLYER) & CALL AIRCRAFT_BALLOON(XTSTEP, & - TDTEXP, TDTMOD, TDTSEG, TDTCUR, & XXHAT, XYHAT, XZZ, XMAP, XLONORI, XLATORI, & XUT, XVT, XWT, XPABST, XTHT, XRT, XSVT, XTKET, XTSRAD, & XRHODREF,XCIT,PSEA=ZSEA(:,:)) @@ -2039,7 +2160,6 @@ IF (LFLYER) & ! IF (LSTATION) & CALL STATION_n(XTSTEP, & - TDTEXP, TDTMOD, TDTSEG, TDTCUR, & XXHAT, XYHAT, XZZ, & XUT, XVT, XWT, XTHT, XRT, XSVT, XTKET, XTSRAD, XPABST ) ! @@ -2050,10 +2170,9 @@ IF (LSTATION) & ! IF (LPROFILER) & CALL PROFILER_n(XTSTEP, & - TDTEXP, TDTMOD, TDTSEG, TDTCUR, & XXHAT, XYHAT, XZZ,XRHODREF, & XUT, XVT, XWT, XTHT, XRT, XSVT, XTKET, XTSRAD, XPABST, & - XAER, XCLDFR, XCIT) + XAER, XCLDFR, XCIT,PSEA=ZSEA(:,:)) ! ! CALL SECOND_MNH2(ZTIME2) @@ -2077,7 +2196,7 @@ ZTIME1 = ZTIME2 ! IF ( .NOT. LIO_NO_WRITE ) THEN IF (NBUMOD==IMI .AND. CBUTYPE/='NONE') THEN - CALL ENDSTEP_BUDGET(TDIAFILE,KTCOUNT,TDTCUR,TDTMOD,XTSTEP,NSV) + CALL ENDSTEP_BUDGET(TDIAFILE,KTCOUNT,TDTCUR,XTSTEP,NSV) END IF END IF ! @@ -2090,8 +2209,7 @@ XT_STEP_BUD = XT_STEP_BUD + ZTIME2 - ZTIME1 + XTIME_BU !* 26. FM FILE CLOSURE ! --------------- ! -IF (GCLOSE_OUT) THEN - GCLOSE_OUT=.FALSE. +IF ( tzbakfile%lopened ) THEN CALL IO_File_close(TZBAKFILE) END IF ! @@ -2100,7 +2218,7 @@ END IF !* 27. CURRENT TIME REFRESH ! -------------------- ! -TDTCUR%TIME=TDTCUR%TIME + XTSTEP +TDTCUR%xtime=TDTCUR%xtime + XTSTEP CALL DATETIME_CORRECTDATE(TDTCUR) ! !------------------------------------------------------------------------------- @@ -2123,11 +2241,10 @@ IF (OEXIT) THEN CALL WRITE_AIRCRAFT_BALLOON(TDIAFILE) CALL WRITE_STATION_n(TDIAFILE) CALL WRITE_PROFILER_n(TDIAFILE) - CALL WRITE_LES_n(TDIAFILE,' ') - CALL WRITE_LES_n(TDIAFILE,'A') - CALL WRITE_LES_n(TDIAFILE,'E') - CALL WRITE_LES_n(TDIAFILE,'H') + call Write_les_n( tdiafile ) +#ifdef MNH_IOLFI CALL MENU_DIACHRO(TDIAFILE,'END') +#endif CALL IO_File_close(TDIAFILE) END IF ! @@ -2163,6 +2280,7 @@ IF (OEXIT) THEN CALL TIME_STAT_ll(XT_ADVUVW,ZTOT, ' ADVECTION UVW','=') CALL TIME_STAT_ll(XT_GRAV,ZTOT, ' GRAVITY','=') CALL TIME_STAT_ll(XT_FORCING,ZTOT, ' FORCING','=') + CALL TIME_STAT_ll(XT_IBM_FORC,ZTOT, ' IBM','=') CALL TIME_STAT_ll(XT_NUDGING,ZTOT, ' NUDGING','=') CALL TIME_STAT_ll(XT_SOURCES,ZTOT, ' DYN_SOURCES','=') CALL TIME_STAT_ll(XT_DIFF,ZTOT, ' NUM_DIFF','=') @@ -2178,6 +2296,7 @@ IF (OEXIT) THEN CALL TIME_STAT_ll(XT_TURB,ZTOT, ' TURB = '//CTURB ,'-') CALL TIME_STAT_ll(XT_MAFL,ZTOT, ' MAFL = '//CSCONV,'-') CALL TIME_STAT_ll(XT_CHEM,ZTOT, ' CHIMIE' ,'-') + CALL TIME_STAT_ll(XT_EOL,ZTOT, ' WIND TURBINE' ,'-') CALL TIMING_LEGEND() CALL TIME_STAT_ll(XT_COUPL,ZTOT, ' SET_COUPLING','=') CALL TIME_STAT_ll(XT_RAD_BOUND,ZTOT, ' RAD_BOUND','=') @@ -2202,12 +2321,13 @@ IF (OEXIT) THEN CALL TIME_STAT_ll(XT_STEP_BUD,ZTOT, ' BUDGETS','=') CALL TIME_STAT_ll(XT_SPECTRA,ZTOT, ' LES','=') CALL TIME_STAT_ll(XT_STEP_MISC,ZTOT, ' MISCELLANEOUS','=') + IF (LIBM) CALL TIME_STAT_ll(XT_IBM_FORC,ZTOT,' IBM FORCING','=') ! ! sum of call subroutine ! ZALL = XT_1WAY + XT_BOUND + XT_STORE + XT_GUESS + XT_2WAY + & XT_ADV + XT_FORCING + XT_NUDGING + XT_SOURCES + XT_DIFF + & - XT_ADVUVW + XT_GRAV + & + XT_ADVUVW + XT_GRAV + XT_IBM_FORC + & XT_RELAX+ XT_PARAM + XT_COUPL + XT_RAD_BOUND+XT_PRESS + & XT_CLOUD+ XT_ELEC + XT_HALO + XT_SPECTRA + XT_STEP_SWA + & XT_STEP_MISC+ XT_STEP_BUD @@ -2246,10 +2366,6 @@ IF (OEXIT) THEN ! CALL TIMING_SEPARATOR('=') ! - ! - ! - CALL IO_File_close(TLUOUT) - IF (IMI==NMODEL) CALL IO_File_close(TLUOUT0) END IF !$acc end data diff --git a/src/ZSOLVER/p_abs.f90 b/src/ZSOLVER/p_abs.f90 index adfa745cb1665ad351177b31e3b00abef6d583b3..cd97868883559d76bd3bd53060b0b899123fb271 100644 --- a/src/ZSOLVER/p_abs.f90 +++ b/src/ZSOLVER/p_abs.f90 @@ -1,13 +1,8 @@ -!MNH_LIC Copyright 1994-2014 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 version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt !MNH_LIC for details. version 1. !----------------------------------------------------------------- -!--------------- special set of characters for RCS information -!----------------------------------------------------------------- -! $Source$ $Revision$ -! MASDEV4_7 solver 2006/05/18 13:07:25 -!----------------------------------------------------------------- ! ################# MODULE MODI_P_ABS ! ################# @@ -16,7 +11,7 @@ INTERFACE ! SUBROUTINE P_ABS (KRR, KRRL, KRRI, PDRYMASST, PREFMASS, PMASS_O_PHI0, & PTHT, PRT, PRHODJ, PRHODREF, PTHETAV, PTHVREF, & - PRVREF, PEXNREF, PPHIT ) + PRVREF, PEXNREF, PPHIT, PPHI0) ! IMPLICIT NONE ! @@ -44,7 +39,7 @@ REAL, DIMENSION(:,:,:), INTENT(IN) :: PEXNREF! Exner function of the ! REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PPHIT ! Perturbation of ! either the Exner function Pi or Pi * Cpd * THvref -! +REAL, INTENT(INOUT) :: PPHI0 ! Phi0 at time t ! ! END SUBROUTINE P_ABS ! @@ -54,7 +49,7 @@ END MODULE MODI_P_ABS ! ####################################################################### SUBROUTINE P_ABS (KRR, KRRL, KRRI, PDRYMASST, PREFMASS, PMASS_O_PHI0, & PTHT, PRT, PRHODJ, PRHODREF, PTHETAV, PTHVREF, & - PRVREF, PEXNREF, PPHIT ) + PRVREF, PEXNREF, PPHIT, PPHI0 ) ! ####################################################################### ! !!**** *P_ABS * - routine to compute the absolute Exner pressure deviation PHI @@ -108,6 +103,8 @@ END MODULE MODI_P_ABS !! from Durran (1989), MAE and DUR respectively !! 15/06/98 (D.Lugato, R.Guivarch) Parallelisation !! J. Colin 07/13 Add LBOUSS +!! J.L Redelsperger 03/2021 Change of one step to pressure computation +!! in order to perform Ocean runs (equivalent to LHE shallow convection) !------------------------------------------------------------------------------- ! !* 0. DECLARATIONS @@ -115,13 +112,13 @@ END MODULE MODI_P_ABS ! USE MODD_CST USE MODD_CONF +USE MODD_DYN_n, ONLY: LOCEAN +USE MODD_IBM_PARAM_n, ONLY: XIBM_LS, LIBM, XIBM_EPSI USE MODD_PARAMETERS -USE MODD_REF, ONLY : LBOUSS +USE MODD_REF, ONLY: LBOUSS ! USE MODE_ll -!JUAN USE MODE_REPRO_SUM -!JUAN ! #ifdef MNH_BITREP USE MODI_BITREP @@ -157,6 +154,7 @@ REAL, DIMENSION(:,:,:), INTENT(IN) :: PRVREF ! vapor mixing ratio REAL, DIMENSION(:,:,:), INTENT(IN) :: PEXNREF! Exner function of the ! reference state ! +REAL, INTENT(INOUT) :: PPHI0 ! PHI0 at t REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PPHIT ! Perturbation of ! either the Exner function Pi or Pi * Cpd * THvref ! @@ -188,7 +186,9 @@ REAL :: ZPI0 ! constant to retrieve the absolute Exner pressure INTEGER :: JWATER ! loop index on the different types of water REAL, DIMENSION(:,:,:) , POINTER , CONTIGUOUS & :: ZRTOT, ZRHOREF, ZWORK +#ifdef MNH_OPENACC INTEGER :: IZRTOT, IZRHOREF, IZWORK +#endif REAL :: ZPHI0 ! INTEGER :: IINFO_ll @@ -276,7 +276,7 @@ IF ( CEQNSYS=='DUR' .OR. CEQNSYS=='MAE' ) THEN DO JJ = IJB,IJE DO JI = IIB,IIE ZMASSGUESS_2D(JI,JJ) = ZMASSGUESS_2D(JI,JJ) + & -#ifndef MNH_OPENACC +#ifndef MNH_BITREP (PEXNREF(JI,JJ,JK)+PPHIT(JI,JJ,JK))**ZCVD_O_RD & #else BR_POW((PEXNREF(JI,JJ,JK)+PPHIT(JI,JJ,JK)),ZCVD_O_RD) & @@ -401,8 +401,13 @@ ELSEIF( CEQNSYS == 'LHE' ) THEN ZRTOT(:,:,:) = 0. END IF ! + IF (LIBM) THEN + WHERE (XIBM_LS(:,:,:,1).GT.-XIBM_EPSI) + ZWORK(:,:,:) = PTHVREF(:,:,:) + ENDWHERE + ENDIF ! - ! compute the absolute pressure function + ! compute the absolute pressure function (LHE equation system case) ! ! ! @@ -424,8 +429,16 @@ ELSEIF( CEQNSYS == 'LHE' ) THEN ZMASSGUESS = SUM_DD_R2_ll(ZMASSGUESS_2D) ZWATERMASST = SUM_DD_R2_ll(ZWATERMASST_2D) ! - ZPHI0 = (PDRYMASST + ZWATERMASST - 2. * PREFMASS + ZMASSGUESS ) / PMASS_O_PHI0 - PPHIT(:,:,:) = PPHIT(:,:,:) + ZPHI0 + ! case shallow bouss : to get the real pressure fluctuation + ! Eq 2.40 p15 : constant not resolved in poisson equation + IF (.NOT. LOCEAN) THEN + PPHI0 = (PDRYMASST + ZWATERMASST - 2. * PREFMASS + ZMASSGUESS ) / PMASS_O_PHI0 + ELSE + ! PPHI0 = 0. => to be possibly modified for ocean LES case + PPHI0=0. + END IF + ! following computation moved in PRESSURE routine (Eq 2.40 bis p15: Phi_total) + ! PPHIT(:,:,:) = PPHIT(:,:,:) + ZPHI0 ! END IF ! diff --git a/src/ZSOLVER/ppm.f90 b/src/ZSOLVER/ppm.f90 index 6c80487d465aecf2a6ddbc78f07daae49da0017f..a713d19be5dcf58bfefef06a0123b576e5f8b9a1 100644 --- a/src/ZSOLVER/ppm.f90 +++ b/src/ZSOLVER/ppm.f90 @@ -2314,47 +2314,6 @@ END WHERE ! ZDQ = ZQR - ZQL #else -!!! -!!! initialize final parabolae parameters -!!! -!!ZQL = ZQL0 -!!ZQR = ZQR0 -!!ZQ6 = ZQ60 -!!! -!!! eliminate over and undershoots and create qL and qR as in Lin96 -!!! -!!!PW: BUG: done like that because if using PGI (tested up to 16.10) -!!! will cause crashes at run (address not mapped) -!!! and problems at compilation -!!$WHERE ( ZDMQ == 0.0 ) -!!$ ZQL = PSRC -!!$ ZQR = PSRC -!!$ ZQ6 = 0.0 -!!$END WHERE -!!#ifndef MNH_BITREP -!!$WHERE ( ( ZDMQ /= 0.0 ) .AND. ( ZQ60*ZDQ < -ZDQ**2 ) ) -!!#else -!!$WHERE ( ( ZDMQ /= 0.0 ) .AND. ( ZQ60*ZDQ < -BR_P2(ZDQ) ) ) -!!#endif -!!$WHERE ( ( ZDMQ /= 0.0 ) .AND. ( ZQ60*ZDQ < -BR_P2(ZDQ) ) ) -!!$ ZQ6 = 3.0*(ZQL0 - PSRC) -!!$ ZQR = ZQL0 - ZQ6 -!!$ ZQL = ZQL0 -!!$END WHERE -!!#ifndef MNH_BITREP -!!$WHERE ( ( ZDMQ /= 0.0 ) .AND. ( ZQ60*ZDQ > ZDQ**2 ) ) -!!#else -!!$WHERE ( ( ZDMQ /= 0.0 ) .AND. ( ZQ60*ZDQ > BR_P2(ZDQ) ) ) -!!#endif -!!$WHERE ( ( ZDMQ /= 0.0 ) .AND. ( ZQ60*ZDQ > BR_P2(ZDQ) ) ) -!!$ ZQ6 = 3.0*(ZQR0 - PSRC) -!!$ ZQL = ZQR0 - ZQ6 -!!$ ZQR = ZQR0 -!!$END WHERE -!!! -!!! recalculate coefficients of the parabolae -!!! -!!ZDQ = ZQR - ZQL !$acc loop independent collapse(3) DO K=1,IKU DO J=1,IJU @@ -2735,6 +2694,7 @@ REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR INTEGER:: IIB,IJB ! Begining useful area in x,y,z directions INTEGER:: IIE,IJE ! End useful area in x,y,z directions INTEGER :: IJS,IJN +INTEGER :: IIU, IJU, IKU ! LOGICAL :: GWEST, GEAST @@ -2755,11 +2715,11 @@ REAL, DIMENSION(:,:,:) :: ZPHAT REAL, DIMENSION(:,:,:) :: ZRHO_MXM, ZCR_MXM , ZCR_DXF INTEGER :: I,J,K ! +LOGICAL, SAVE :: GFIRST_CALL_PPM_S0_X = .TRUE. REAL, DIMENSION(:,:) :: ZPSRC_HALO2_WEST #endif TYPE(HALO2LIST_ll), SAVE , POINTER :: TZ_PSRC_HALO2_ll ! halo2 for PSRC -LOGICAL, SAVE :: GFIRST_CALL_PPM_S0_X = .TRUE. ! REAL , POINTER , CONTIGUOUS , DIMENSION(:,:) :: ZWEST @@ -2788,21 +2748,28 @@ IJN=IJE ! GWEST = LWEST_ll() GEAST = LEAST_ll() +! +IIU = SIZE( PSRC, 1 ) +IJU = SIZE( PSRC, 2 ) +IKU = SIZE( PSRC, 3 ) #endif ! !BEG JUAN PPM_LL ! !* initialise & update halo & halo2 for PSRC ! -!!$CALL GET_HALO2(PSRC, TZ_PSRC_HALO2_ll, HNAME='PSRC') +#ifndef MNH_OPENACC +CALL GET_HALO2(PSRC, TZ_PSRC_HALO2_ll, HNAME='PSRC') !!$ZPSRC_HALO2_WEST(:,:) = TZ_PSRC_HALO2_ll%HALO2%WEST(:,:) !!$!$acc update device (ZPSRC_HALO2_WEST) +#else IF (GFIRST_CALL_PPM_S0_X) THEN GFIRST_CALL_PPM_S0_X = .FALSE. NULLIFY(TZ_PSRC_HALO2_ll) CALL INIT_HALO2_ll(TZ_PSRC_HALO2_ll,1,IIU,IJU,IKU) END IF CALL GET_HALO2_DF(PSRC, TZ_PSRC_HALO2_ll, HNAME='PSRC') +#endif ZWEST => TZ_PSRC_HALO2_ll%HALO2%WEST !$acc kernels ZPSRC_HALO2_WEST(:,:) = ZWEST(:,:) @@ -3182,6 +3149,7 @@ INTEGER:: IIB,IJB ! Begining useful area in x,y,z directions INTEGER:: IIE,IJE ! End useful area in x,y,z directions INTEGER :: IJS,IJN INTEGER :: IIW,IIA +INTEGER :: IIU, IJU, IKU ! LOGICAL :: GNORTH, GSOUTH ! @@ -3204,7 +3172,6 @@ REAL, DIMENSION(:,:,:) :: ZFPOS, ZFNEG REAL, DIMENSION(:,:,:) :: ZPHAT ! TYPE(HALO2LIST_ll), SAVE ,POINTER :: TZ_PSRC_HALO2_ll ! halo2 for PSRC -LOGICAL, SAVE :: GFIRST_CALL_PPM_S0_Y = .TRUE. TYPE(HALO2LIST_ll), POINTER :: TZ_PHAT_HALO2_ll ! halo2 for ZPHAT ! @@ -3212,6 +3179,7 @@ REAL, DIMENSION(:,:,:) :: ZRHO_MYM , ZCR_MYM , ZCR_DYF ! INTEGER :: I,J,K ! +LOGICAL, SAVE :: GFIRST_CALL_PPM_S0_Y = .TRUE. REAL, DIMENSION(:,:) :: ZPSRC_HALO2_SOUTH #endif ! @@ -3242,6 +3210,10 @@ IIA=IIE ! GNORTH = LNORTH_ll() GSOUTH = LSOUTH_ll() +! +IIU = SIZE( PSRC, 1 ) +IJU = SIZE( PSRC, 2 ) +IKU = SIZE( PSRC, 3 ) #endif ! !------------------------------------------------------------------------------- @@ -3256,15 +3228,18 @@ IF ( L2D ) THEN ! RETURN ELSE !not L2D ! -!!$CALL GET_HALO2(PSRC, TZ_PSRC_HALO2_ll, HNAME='PSRC') +#ifndef MNH_OPENACC +CALL GET_HALO2(PSRC, TZ_PSRC_HALO2_ll, HNAME='PSRC') !!$ZPSRC_HALO2_SOUTH(:,:) = TZ_PSRC_HALO2_ll%HALO2%SOUTH(:,:) !!$!$acc update device (ZPSRC_HALO2_SOUTH) +#else IF (GFIRST_CALL_PPM_S0_Y) THEN GFIRST_CALL_PPM_S0_Y = .FALSE. NULLIFY(TZ_PSRC_HALO2_ll) CALL INIT_HALO2_ll(TZ_PSRC_HALO2_ll,1,IIU,IJU,IKU) END IF CALL GET_HALO2_DF(PSRC, TZ_PSRC_HALO2_ll, HNAME='PSRC') +#endif ZSOUTH => TZ_PSRC_HALO2_ll%HALO2%SOUTH(:,:) ! ! Initialize with relalistic value all work array diff --git a/src/ZSOLVER/pressurez.f90 b/src/ZSOLVER/pressurez.f90 index b7221fa4acb1ab350c0658a531032e553d8f22b5..ea71d1adf9919117e5263b2192bdb44f4b109592 100644 --- a/src/ZSOLVER/pressurez.f90 +++ b/src/ZSOLVER/pressurez.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. @@ -231,16 +231,21 @@ END MODULE MODI_PRESSUREZ !! Philippe Wautelet: 22/01/2019: use standard FLUSH statement instead of non standard intrinsics ! P. Wautelet 10/04/2019: replace ABORT and STOP calls by Print_msg ! P. Wautelet 20/05/2019: add name argument to ADDnFIELD_ll + new ADD4DFIELD_ll subroutine +! P. Wautelet 28/01/2020: use the new data structures and subroutines for budgets for U +!! JL Redelsperger 03/2021 : Shallow convection case added in LHE case: +!! working for both atmosphere and ocean cases !------------------------------------------------------------------------------- ! !* 0. DECLARATIONS ! ------------ ! USE MODD_ARGSLIST_ll, ONLY: LIST_ll -USE MODD_BUDGET +use modd_budget, only: lbudget_u, lbudget_v, lbudget_w, NBUDGET_U, NBUDGET_V, NBUDGET_W, tbudgets USE MODD_CST USE MODD_CONF -USE MODD_DYN_n, ONLY: LRES, XRES +USE MODD_DYN_n, ONLY: LRES, XRES,LOCEAN +USE MODD_FIELD_n, ONLY: XPHIT +USE MODD_IBM_PARAM_n, ONLY: LIBM, NIBM_ITR, XIBM_EPSI, XIBM_LS, XIBM_SU USE MODD_LUNIT_n, ONLY: TLUOUT USE MODD_MPIF USE MODD_PARAMETERS @@ -248,35 +253,39 @@ use modd_precision, only: MNHREAL_MPI USE MODD_REF, ONLY: LBOUSS USE MODD_VAR_ll, ONLY: NMNH_COMM_WORLD , NPROC ! +use mode_budget, only: Budget_store_end USE MODE_ll +#ifdef MNH_OPENACC +USE MODE_MNH_ZWORK, ONLY: MNH_ALLOCATE_ZT3D, MNH_REL_ZT3D +#endif USE MODE_MPPDB USE MODE_MSG USE MODE_SUM2_ll, ONLY: GMAXLOC_ll ! -USE MODI_BUDGET +#ifdef MNH_BITREP +USE MODI_BITREP +#endif USE MODI_CONJGRAD -USE MODI_ZCONJGRAD USE MODI_CONRESOL USE MODI_CONRESOLZ USE MODI_FLAT_INV -USE MODI_ZSOLVER USE MODI_FLAT_INVZ USE MODI_GDIV +#ifdef MNH_OPENACC +USE MODI_GET_HALO +#endif USE MODI_GRADIENT_M +USE MODI_IBM_BALANCE USE MODI_MASS_LEAK USE MODI_P_ABS USE MODI_RICHARDSON USE MODI_SHUMAN #ifdef MNH_OPENACC USE MODI_SHUMAN_DEVICE -USE MODI_GET_HALO -#endif -! -#ifdef MNH_BITREP -USE MODI_BITREP -USE MODE_MNH_ZWORK , ONLY : MNH_ALLOCATE_ZT3D , MNH_REL_ZT3D #endif -! +USE MODI_ZCONJGRAD +USE MODI_ZSOLVER + IMPLICIT NONE ! !* 0.1 declarations of arguments @@ -368,7 +377,9 @@ REAL, OPTIONAL :: PRESIDUAL ! REAL, DIMENSION(:,:,:) , POINTER , CONTIGUOUS :: ZDV_SOURCE ! ! divergence of the sources +#ifdef MNH_OPENACC INTEGER :: IZDV_SOURCE +#endif ! INTEGER :: IIB ! indice I for the first inner mass point along x INTEGER :: IIE ! indice I for the last inner mass point along x @@ -384,8 +395,11 @@ REAL, DIMENSION(:,:,:) , POINTER , CONTIGUOUS :: ZTHETAV, & ZPHIT ! MAE + DUR => Exner function perturbation ! LHE => Exner function perturbation * CPD * THVREF -INTEGER :: IZTHETAV,IZPHIT +#ifdef MNH_OPENACC +INTEGER :: IZTHETAV, IZPHIT +#endif ! +REAL :: ZPHI0 REAL :: ZRV_OV_RD ! XRV / XRD REAL :: ZMAXVAL, ZMAXRES, ZMAX,ZMAX_ll ! for print INTEGER, DIMENSION(3) :: IMAXLOC ! purpose @@ -404,10 +418,12 @@ TYPE(LIST_ll), POINTER :: TZFIELDS_ll, TZFIELDS2_ll ! list of fields to exchang INTEGER :: IIB_I,IIE_I,IJB_I,IJE_I INTEGER :: IIMAX_ll,IJMAX_ll ! +#ifdef MNH_OPENACC REAL, DIMENSION(:,:,:) , POINTER , CONTIGUOUS :: ZPRHODJ,ZMXM_PRHODJ,ZMZM_PRHODJ,ZGZ_M_W -INTEGER :: IZPRHODJ,IZMXM_PRHODJ,IZMZM_PRHODJ,IZGZ_M_W +INTEGER :: IZPRHODJ, IZMXM_PRHODJ, IZMZM_PRHODJ, IZGZ_M_W REAL, DIMENSION(:,:,:) , POINTER , CONTIGUOUS :: ZMYM_PRHODJ INTEGER :: IZMYM_PRHODJ +#endif ! ! LOGICAL :: GWEST,GEAST,GSOUTH,GNORTH @@ -425,7 +441,7 @@ ILUOUT = TLUOUT%NLU ! CALL GET_GLOBALDIMS_ll (IIMAX_ll,IJMAX_ll) IF ( ( MIN(IIMAX_ll,IJMAX_ll) < NPROC ) .AND. ( HPRESOPT /= 'ZRESI' ) .AND. ( HPRESOPT /= 'ZSOLV' )) THEN - WRITE(UNIT=ILUOUT,FMT=*) 'ERROR IN PRESSUREZ:: YOU WANT TO USE TO MANY PROCESSOR WITHOUT CPRESOPT="ZRESI/ZSOLV" ' + WRITE(UNIT=ILUOUT,FMT=*) 'ERROR IN PRESSUREZ:: YOU WANT TO USE TO MANY PROCESSES WITHOUT CPRESOPT="ZRESI/ZSOLV" ' WRITE(UNIT=ILUOUT,FMT=*) 'MIN(IIMAX_ll,IJMAX_ll)=',MIN(IIMAX_ll,IJMAX_ll),' < NPROC =', NPROC WRITE(UNIT=ILUOUT,FMT=*) 'YOU HAVE TO SET CPRESOPT="ZRESI => JOB ABORTED ' CALL PRINT_MSG(NVERB_FATAL,'GEN','PRESSUREZ','') @@ -446,6 +462,11 @@ GNORTH2D = ( L2D .AND. LNORTH_ll() ) ! GPRVREF0 = ( SIZE(PRVREF,1) == 0 ) ! +#ifndef MNH_OPENACC +ALLOCATE( ZDV_SOURCE(IIU, IJU, IKU) ) +ALLOCATE( ZTHETAV (IIU, IJU, IKU) ) +ALLOCATE( ZPHIT (IIU, IJU, IKU) ) +#else IZDV_SOURCE = MNH_ALLOCATE_ZT3D(ZDV_SOURCE ,IIU,IJU,IKU ) IZTHETAV = MNH_ALLOCATE_ZT3D(ZTHETAV ,IIU,IJU,IKU ) IZPHIT = MNH_ALLOCATE_ZT3D(ZPHIT ,IIU,IJU,IKU ) @@ -455,6 +476,7 @@ IZMXM_PRHODJ = MNH_ALLOCATE_ZT3D( ZMXM_PRHODJ,IIU,IJU,IKU ) IZMZM_PRHODJ = MNH_ALLOCATE_ZT3D( ZMZM_PRHODJ,IIU,IJU,IKU ) IZGZ_M_W = MNH_ALLOCATE_ZT3D( ZGZ_M_W,IIU,IJU,IKU ) IZMYM_PRHODJ = MNH_ALLOCATE_ZT3D( ZMYM_PRHODJ,IIU,IJU,IKU ) +#endif ! IF (GFIRST_CALL_PRESSUREZ) THEN GFIRST_CALL_PRESSUREZ = .FALSE. @@ -466,8 +488,12 @@ ZPABS_N(:,:) = 0. ZPABS_E(:,:) = 0. ZPABS_W(:,:) = 0. !$acc end kernels -! -! + +! Done in model_n before call to Rad_bound +! if ( lbudget_u ) call Budget_store_init( tbudgets(NBUDGET_U), 'PRES', prus(:, :, :) ) +! if ( lbudget_v ) call Budget_store_init( tbudgets(NBUDGET_V), 'PRES', prvs(:, :, :) ) +! if ( lbudget_w ) call Budget_store_init( tbudgets(NBUDGET_W), 'PRES', prws(:, :, :) ) + !------------------------------------------------------------------------------- ! !* 3. COMPUTE THE LINEIC MASS @@ -482,6 +508,11 @@ END IF !* 4. COMPUTE THE FORCING TERM FOR THE PRESSURE EQUATION ! -------------------------------------------------- ! +IF (LIBM) THEN + WHERE(XIBM_LS(:,:,:,2).GT.-XIBM_EPSI) PRUS(:,:,:) = 0. + WHERE(XIBM_LS(:,:,:,3).GT.-XIBM_EPSI) PRVS(:,:,:) = 0. + WHERE(XIBM_LS(:,:,:,4).GT.-XIBM_EPSI) PRWS(:,:,:) = 0. +ENDIF ! CALL MPPDB_CHECK3D(PRUS,"pressurez 4-before update_halo_ll::PRUS",PRECISION) CALL MPPDB_CHECK3D(PRVS,"pressurez 4-before update_halo_ll::PRVS",PRECISION) @@ -504,6 +535,10 @@ CALL MPPDB_CHECK3D(PRWS,"pressurez 4-after update_halo_ll::PRWS",PRECISION) ! CALL GDIV(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRUS,PRVS,PRWS,ZDV_SOURCE) ! +IF (LIBM) THEN + CALL IBM_BALANCE(XIBM_LS,XIBM_SU,PRUS,PRVS,PRWS,ZDV_SOURCE) +ENDIF +! ! The non-homogenous Neuman problem is transformed in an homogenous Neuman ! problem in the non-periodic cases IF ( GWEST ) THEN @@ -527,6 +562,23 @@ IF ( GNORTH ) THEN !$acc end kernels END IF !$acc wait + +IF (LIBM) THEN + ! + IF (HLBCX(1) == 'CYCL') THEN + IF (LWEST_ll()) ZDV_SOURCE(IIB-1,:,:) = ZDV_SOURCE(IIB-1,:,:)*XIBM_SU(IIB,:,:,1) + IF (LEAST_ll()) ZDV_SOURCE(IIE+1,:,:) = ZDV_SOURCE(IIE+1,:,:)*XIBM_SU(IIE,:,:,1) + ENDIF + ! + IF (HLBCY(1) == 'CYCL') THEN + IF (LSOUTH_ll()) ZDV_SOURCE(:,IJB-1,:) = ZDV_SOURCE(:,IJB-1,:)*XIBM_SU(:,IJB,:,1) + IF (LNORTH_ll()) ZDV_SOURCE(:,IJE+1,:) = ZDV_SOURCE(:,IJE+1,:)*XIBM_SU(:,IJE,:,1) + ENDIF + ! + ZDV_SOURCE(:,:,IKB-1) = ZDV_SOURCE(:,:,IKB-1)*XIBM_SU(:,:,IKB,1) + ZDV_SOURCE(:,:,IKE+1) = ZDV_SOURCE(:,:,IKE+1)*XIBM_SU(:,:,IKE,1) + ! +ENDIF ! !------------------------------------------------------------------------------- ! @@ -557,23 +609,43 @@ IF(CEQNSYS=='MAE' .OR. CEQNSYS=='DUR') THEN END IF !$acc end kernels ! + IF (LIBM) THEN + WHERE (XIBM_LS(:,:,:,1).GT.-XIBM_EPSI) + ZTHETAV(:,:,:) = PTHVREF(:,:,:) + ENDWHERE + ENDIF + ! #ifndef MNH_OPENACC +#ifndef MNH_BITREP ZPHIT(:,:,:)=(PPABST(:,:,:)/XP00)**(XRD/XCPD)-PEXNREF(:,:,:) +#else + ZPHIT(:,:,:)=BR_POW(PPABST(:,:,:)/XP00,XRD/XCPD)-PEXNREF(:,:,:) +#endif #else !$acc kernels DO CONCURRENT ( JI=1:IIU,JJ=1:IJU,JK=1:IKU ) - ZPHIT(JI,JJ,JK)=BR_POW((PPABST(JI,JJ,JK)/XP00),(XRD/XCPD))-PEXNREF(JI,JJ,JK) +#ifndef MNH_BITREP + ZPHIT(JI,JJ,JK)=(PPABST(JI,JJ,JK)/XP00)**(XRD/XCPD)-PEXNREF(JI,JJ,JK) +#else + ZPHIT(JI,JJ,JK)=BR_POW((PPABST(JI,JJ,JK)/XP00),(XRD/XCPD))-PEXNREF(JI,JJ,JK) +#endif END DO !$acc end kernels #endif ! ELSEIF(CEQNSYS=='LHE') THEN - ZPHIT(:,:,:)= ((PPABST(:,:,:)/XP00)**(XRD/XCPD)-PEXNREF(:,:,:)) & - * XCPD * PTHVREF(:,:,:) - ! + IF ( .NOT. LOCEAN) THEN + ZPHIT(:,:,:)= ((PPABST(:,:,:)/XP00)**(XRD/XCPD)-PEXNREF(:,:,:)) & + * XCPD * PTHVREF(:,:,:) + ELSE + ! Field at T- DT for LHE anelastic approx + ! not used in ocean case (flat LHE) + ZPHIT(:,:,:)=0. + END IF +! END IF ! -IF(CEQNSYS=='LHE'.AND. LFLAT .AND. LCARTESIAN) THEN +IF(CEQNSYS=='LHE'.AND. LFLAT .AND. LCARTESIAN .AND. .NOT. LIBM) THEN ! flat cartesian LHE case -> exact solution IF ( HPRESOPT /= "ZRESI" ) THEN CALL FLAT_INV(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOT,PAF,PBF,PCF, & @@ -660,6 +732,32 @@ IF ( GNORTH2D ) THEN END IF !$acc wait ! +IF (LIBM) THEN + ! + IF ( HLBCX(1) == 'CYCL' ) THEN + IF (LWEST_ll()) THEN + ZPHIT(IIB-1,:,:) = ZPHIT(IIB,:,:)*(1.-XIBM_SU(IIB,:,:,1))+XIBM_SU(IIB,:,:,1)*ZPHIT(IIB-1,:,:) + ENDIF + IF (LEAST_ll()) THEN + ZPHIT(IIE+1,:,:) = ZPHIT(IIE,:,:)*(1.-XIBM_SU(IIE,:,:,1))+XIBM_SU(IIE,:,:,1)*ZPHIT(IIE+1,:,:) + ENDIF + ENDIF + ! + IF ( HLBCY(1) == 'CYCL' ) THEN + IF (LSOUTH_ll()) THEN + ZPHIT(:,IJB-1,:) = ZPHIT(:,IJB,:)*(1.-XIBM_SU(:,IJB,:,1))+XIBM_SU(:,IJB,:,1)*ZPHIT(:,IJB-1,:) + ENDIF + IF (LNORTH_ll()) THEN + ZPHIT(:,IJE+1,:) = ZPHIT(:,IJE,:)*(1.-XIBM_SU(:,IJE,:,1))+XIBM_SU(:,IJE,:,1)*ZPHIT(:,IJE+1,:) + ENDIF + ENDIF + ! + !-------------Bottom Boundary conditions + ZPHIT(:,:,IKB-1) = ZPHIT(:,:,IKB-1)*XIBM_SU(:,:,IKB,1)+(1.-XIBM_SU(:,:,IKB,1))*ZPHIT(:,:,IKB) + ZPHIT(:,:,IKE+1) = ZPHIT(:,:,IKE+1)*XIBM_SU(:,:,IKE,1)+(1.-XIBM_SU(:,:,IKE,1))*ZPHIT(:,:,IKE) + ! +ENDIF +! #ifndef MNH_OPENACC ZDV_SOURCE = GX_M_U(1,IKU,1,ZPHIT,PDXX,PDZZ,PDZX) #else @@ -699,8 +797,10 @@ END IF ! CALL MPPDB_CHECK3DM("before MXM PRESSUREZ :PRU/V/WS",PRECISION,PRUS,PRVS,PRWS) 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) +#ifndef MNH_OPENACC + PRUS = PRUS - MXM(PRHODJ * XCPD * ZTHETAV) * ZDV_SOURCE + PRWS = PRWS - MZM(PRHODJ * XCPD * ZTHETAV) * GZ_M_W(1,IKU,1,ZPHIT,PDZZ) +#else !$acc kernels ZPRHODJ = PRHODJ * XCPD * ZTHETAV !$acc end kernels @@ -711,6 +811,7 @@ IF(CEQNSYS=='MAE' .OR. CEQNSYS=='DUR') THEN PRUS = PRUS - ZMXM_PRHODJ * ZDV_SOURCE PRWS = PRWS - ZMZM_PRHODJ * ZGZ_M_W !$acc end kernels +#endif ELSEIF(CEQNSYS=='LHE') THEN PRUS = PRUS - MXM(PRHODJ) * ZDV_SOURCE PRWS = PRWS - MZM(PRHODJ) * GZ_M_W(1,IKU,1,ZPHIT,PDZZ) @@ -757,11 +858,14 @@ IF(.NOT. L2D) THEN ! CALL MPPDB_CHECK3DM("before MYM PRESSUREZ :PRU/V/WS",PRECISION,PRUS,PRVS,PRWS) IF(CEQNSYS=='MAE' .OR. CEQNSYS=='DUR') THEN -!!$ PRVS = PRVS - MYM(PRHODJ * XCPD * ZTHETAV) * ZDV_SOURCE +#ifndef MNH_OPENACC + PRVS = PRVS - MYM(PRHODJ * XCPD * ZTHETAV) * ZDV_SOURCE +#else CALL MYM_DEVICE(ZPRHODJ,ZMYM_PRHODJ) !$acc kernels PRVS = PRVS - ZMYM_PRHODJ * ZDV_SOURCE !$acc end kernels +#endif ELSEIF(CEQNSYS=='LHE') THEN PRVS = PRVS - MYM(PRHODJ) * ZDV_SOURCE END IF @@ -793,12 +897,16 @@ CALL GET_HALO_D(PRWS,HNAME='PRESSUREZ::PRWS' ) ! compute the residual divergence CALL GDIV(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,PRUS,PRVS,PRWS,ZDV_SOURCE) ! +IF (LIBM) THEN + ZDV_SOURCE(:,:,:)=ZDV_SOURCE(:,:,:)*XIBM_SU(:,:,:,2) +ENDIF +! IF ( CEQNSYS=='DUR' ) THEN - IF ( GPRVREF0 ) THEN + IF ( GPRVREF0 ) THEN !$acc kernels ZDV_SOURCE=ZDV_SOURCE/PRHODJ/XTH00*PRHODREF*PTHVREF !$acc end kernels - ELSE + ELSE !$acc kernels ZDV_SOURCE=ZDV_SOURCE/PRHODJ/XTH00*PRHODREF*PTHVREF*(1.+PRVREF) !$acc end kernels @@ -838,12 +946,16 @@ IF (OITRADJ) THEN ENDIF ENDIF ! +IF (LIBM) THEN + KITR=MIN(NIBM_ITR,KITR) +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') +if ( lbudget_u ) call Budget_store_end( tbudgets(NBUDGET_U), 'PRES', prus(:, :, :) ) +if ( lbudget_v ) call Budget_store_end( tbudgets(NBUDGET_V), 'PRES', prvs(:, :, :) ) +if ( lbudget_w ) call Budget_store_end( tbudgets(NBUDGET_W), 'PRES', prws(:, :, :) ) ! !------------------------------------------------------------------------------- ! @@ -859,18 +971,29 @@ IF ((ZMAX_ll > 1.E-12) .AND. KTCOUNT >0 ) THEN !IF ( KTCOUNT >0 .AND. .NOT.LBOUSS ) THEN CALL P_ABS ( KRR, KRRL, KRRI, PDRYMASST, PREFMASS, PMASS_O_PHI0, & PTHT, PRT, PRHODJ, PRHODREF, ZTHETAV, PTHVREF, & - PRVREF, PEXNREF, ZPHIT ) + PRVREF, PEXNREF, ZPHIT, ZPHI0 ) ! IF(CEQNSYS=='MAE' .OR. CEQNSYS=='DUR') THEN -#ifndef MNH_OPENACC + !$acc kernels +#ifndef MNH_BITREP PPABST(:,:,:)=XP00*(ZPHIT+PEXNREF)**(XCPD/XRD) #else - !$acc kernels PPABST(:,:,:)=XP00*BR_POW((ZPHIT+PEXNREF),(XCPD/XRD)) - !$acc end kernels #endif + !$acc end kernels ELSEIF(CEQNSYS=='LHE') THEN - PPABST(:,:,:)=XP00*(ZPHIT/(XCPD*PTHVREF)+PEXNREF)**(XCPD/XRD) + IF (.NOT. LOCEAN) THEN + ! Deep atmosphere case : computing of PI fluctuation ; ZPHI0 (computed in P_ABS routine) is added + XPHIT(:,:,:) = (ZPHIT+ZPHI0)/(XCPD*PTHVREF) + ! Absolute Pressure + PPABST(:,:,:)=XP00*(XPHIT(:,:,:)+PEXNREF)**(XCPD/XRD) + ! Computing press fluctuation + XPHIT(:,:,:) = PPABST(:,:,:) - XP00*PEXNREF**(XCPD/XRD) + ELSE +! Shallow atmosphere ou ocean + XPHIT(:,:,:) = (ZPHIT+ZPHI0)*PRHODREF + PPABST(:,:,:)=XPHIT(:,:,:) + XP00*PEXNREF**(XCPD/XRD) + END IF ENDIF ! IF( GWEST ) THEN diff --git a/src/ZSOLVER/qlap.f90 b/src/ZSOLVER/qlap.f90 index 5f76e50d9b5f77db9717497307716ce2d7e2fef0..050ad22829f4b3840b6b2298cab4807f1d4d3976 100644 --- a/src/ZSOLVER/qlap.f90 +++ b/src/ZSOLVER/qlap.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. @@ -137,6 +137,7 @@ END MODULE MODI_QLAP !! 06/12 V.Masson : update_halo due to CONTRAV changes !! J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1 ! P. Wautelet 20/05/2019: add name argument to ADDnFIELD_ll + new ADD4DFIELD_ll subroutine +! F. Auguste 02/21: add IBM !------------------------------------------------------------------------------- ! !* 0. DECLARATIONS @@ -153,6 +154,9 @@ USE MODI_SHUMAN ! USE MODE_MPPDB ! +USE MODD_IBM_PARAM_n, ONLY: XIBM_LS, LIBM, XIBM_SU +USE MODI_IBM_BALANCE +! IMPLICIT NONE ! !* 0.1 declarations of arguments @@ -187,7 +191,7 @@ INTEGER :: IIU,IJU,IKU ! I,J,K array sizes INTEGER :: JK,JJ,JI ! vertical loop index TYPE(LIST_ll), POINTER :: TZFIELDS_ll ! list of fields to exchange INTEGER :: IINFO_ll -INTEGER :: IIB,IIE,IJB,IJE +INTEGER :: IIB,IIE,IJB,IJE,IKB,IKE !------------------------------------------------------------------------------- ! ! @@ -197,6 +201,8 @@ INTEGER :: IIB,IIE,IJB,IJE CALL GET_DIM_EXT_ll('B',IIU,IJU) CALL GET_INDICE_ll (IIB,IJB,IIE,IJE) IKU=SIZE(PY,3) +IKE = IKU - JPVEXT +IKB = 1 + JPVEXT ! ZU = GX_M_U(1,IKU,1,PY,PDXX,PDZZ,PDZX) CALL MPPDB_CHECK3D(ZU,'QLAP::ZU',PRECISION) @@ -284,7 +290,26 @@ CALL ADD3DFIELD_ll( TZFIELDS_ll, ZW, 'QLAP::ZW' ) CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll) CALL CLEANLIST_ll(TZFIELDS_ll) ! -CALL GDIV(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,ZU,ZV,ZW,PQLAP) +CALL GDIV(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,ZU,ZV,ZW,PQLAP) +! +IF (LIBM) THEN + ! + CALL IBM_BALANCE(XIBM_LS,XIBM_SU,ZU,ZV,ZW,PQLAP) + ! + PQLAP(:,:,IKB-1) = PQLAP(:,:,IKB-1)*XIBM_SU(:,:,IKB,1) + PQLAP(:,:,IKE+1) = PQLAP(:,:,IKE+1)*XIBM_SU(:,:,IKE,1) + ! + IF ( HLBCX(1) /= 'CYCL' ) THEN + IF(LWEST_ll()) PQLAP(IIB-1,:,:) = PQLAP(IIB-1,:,:)*XIBM_SU(IIB,:,:,1) + IF(LEAST_ll()) PQLAP(IIE+1,:,:) = PQLAP(IIE+1,:,:)*XIBM_SU(IIE,:,:,1) + ENDIF + ! + IF ( HLBCY(1) /= 'CYCL' ) THEN + IF (LSOUTH_ll()) PQLAP(:,IJB-1,:) = PQLAP(:,IJB-1,:)*XIBM_SU(:,IJB,:,1) + IF (LNORTH_ll()) PQLAP(:,IJE+1,:) = PQLAP(:,IJE+1,:)*XIBM_SU(:,IJE,:,1) + ENDIF + ! +ENDIF ! !------------------------------------------------------------------------------- ! @@ -361,6 +386,7 @@ END FUNCTION QLAP !! 06/12 V.Masson : update_halo due to CONTRAV changes !! J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1 ! P. Wautelet 20/05/2019: add name argument to ADDnFIELD_ll + new ADD4DFIELD_ll subroutine +! F. Auguste 02/21: add IBM !------------------------------------------------------------------------------- ! !* 0. DECLARATIONS @@ -381,6 +407,8 @@ USE MODE_MPPDB USE MODE_MNH_ZWORK , ONLY : MNH_ALLOCATE_ZT3D , MNH_REL_ZT3D ! USE MODI_GET_HALO +USE MODD_IBM_PARAM_n, ONLY: XIBM_LS, LIBM, XIBM_SU +USE MODI_IBM_BALANCE ! IMPLICIT NONE ! @@ -418,7 +446,7 @@ INTEGER :: IIU,IJU,IKU ! I,J,K array sizes INTEGER :: JK,JJ,JI ! vertical loop index TYPE(LIST_ll), POINTER :: TZFIELDS_ll ! list of fields to exchange INTEGER :: IINFO_ll -INTEGER :: IIB,IIE,IJB,IJE +INTEGER :: IIB,IIE,IJB,IJE,IKB,IKE ! REAL, DIMENSION(:,:,:), pointer , contiguous :: ZMXM,ZMYM,ZMZM,ZRHODJ,ZGZMW INTEGER :: IZMXM,IZMYM,IZMZM,IZRHODJ,IZGZMW @@ -433,6 +461,8 @@ LOGICAL :: GWEST,GEAST,GSOUTH,GNORTH CALL GET_DIM_EXT_ll('B',IIU,IJU) CALL GET_INDICE_ll (IIB,IJB,IIE,IJE) IKU=SIZE(PY,3) +IKE = IKU - JPVEXT +IKB = 1 + JPVEXT ! GWEST = ( HLBCX(1) /= 'CYCL' .AND. LWEST_ll() ) GEAST = ( HLBCX(2) /= 'CYCL' .AND. LEAST_ll() ) @@ -559,7 +589,26 @@ CALL GET_HALO_D(ZU,HNAME='QLAP::ZU') CALL GET_HALO_D(ZV,HNAME='QLAP::ZV') CALL GET_HALO_D(ZW,HNAME='QLAP::ZW') ! -CALL GDIV(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,ZU,ZV,ZW,PQLAP) +CALL GDIV(HLBCX,HLBCY,PDXX,PDYY,PDZX,PDZY,PDZZ,ZU,ZV,ZW,PQLAP) +! +IF (LIBM) THEN + ! + CALL IBM_BALANCE(XIBM_LS,XIBM_SU,ZU,ZV,ZW,PQLAP) + ! + PQLAP(:,:,IKB-1) = PQLAP(:,:,IKB-1)*XIBM_SU(:,:,IKB,1) + PQLAP(:,:,IKE+1) = PQLAP(:,:,IKE+1)*XIBM_SU(:,:,IKE,1) + ! + IF ( HLBCX(1) /= 'CYCL' ) THEN + IF(LWEST_ll()) PQLAP(IIB-1,:,:) = PQLAP(IIB-1,:,:)*XIBM_SU(IIB,:,:,1) + IF(LEAST_ll()) PQLAP(IIE+1,:,:) = PQLAP(IIE+1,:,:)*XIBM_SU(IIE,:,:,1) + ENDIF + ! + IF ( HLBCY(1) /= 'CYCL' ) THEN + IF (LSOUTH_ll()) PQLAP(:,IJB-1,:) = PQLAP(:,IJB-1,:)*XIBM_SU(:,IJB,:,1) + IF (LNORTH_ll()) PQLAP(:,IJE+1,:) = PQLAP(:,IJE+1,:)*XIBM_SU(:,IJE,:,1) + ENDIF + ! +ENDIF ! CALL MNH_REL_ZT3D ( IZU,IZV,IZW, IZMXM,IZMYM,IZMZM,IZRHODJ,IZGZMW ) !------------------------------------------------------------------------------- diff --git a/src/ZSOLVER/read_exsegn.f90 b/src/ZSOLVER/read_exsegn.f90 index 85e867132c14957b8d40a881ad9fe9d337879228..53e7a74fa142fefbad6ffa478f7b10688f6b1e91 100644 --- a/src/ZSOLVER/read_exsegn.f90 +++ b/src/ZSOLVER/read_exsegn.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. @@ -293,95 +293,114 @@ END MODULE MODI_READ_EXSEG_n !! Modification 01/2019 (R. Honnert) remove SURF in CMF_UPDRAFT !! Bielli S. 02/2019 Sea salt : significant sea wave height influences salt emission; 5 salt modes ! P. Wautelet 10/04/2019: replace ABORT and STOP calls by Print_msg -!!------------------------------------------------------------------------------ +! C. Lac 11/2019: correction in the drag formula and application to building in addition to tree +! Q. Rodier 03/2020: add abort if use of any LHORELAX and cyclic conditions +! F.Auguste 02/2021: add IBM +! T.Nagel 02/2021: add turbulence recycling +! E.Jezequel 02/2021: add stations read from CSV file +! P. Wautelet 09/03/2021: simplify allocation of scalar variable names +! P. Wautelet 09/03/2021: move some chemistry initializations to ini_nsv +! P. Wautelet 10/03/2021: move scalar variable name initializations to ini_nsv +! R. Honnert 23/04/2021: add ADAP mixing length and delete HRIO and BOUT from CMF_UPDRAFT +! S. Riette 11/05/2021 HighLow cloud +!------------------------------------------------------------------------------ ! !* 0. DECLARATIONS ! ------------ -USE MODD_PARAMETERS +USE MODD_BLOWSNOW +USE MODD_BUDGET +USE MODD_CH_AEROSOL +USE MODD_CH_M9_n, ONLY : NEQ +USE MODD_CONDSAMP USE MODD_CONF -USE MODD_CONFZ USE MODD_CONF_n, ONLY: CSTORAGE_TYPE +USE MODD_CONFZ +! USE MODD_DRAG_n +USE MODD_DUST +USE MODD_DYN +USE MODD_DYN_n, ONLY : LHORELAX_SVLIMA +#ifdef MNH_FOREFIRE +USE MODD_FOREFIRE +#endif +USE MODD_GET_n +USE MODD_GR_FIELD_n USE MODD_IO, ONLY: TFILEDATA USE MODD_LUNIT_n, ONLY: TLUOUT +USE MODD_NSV,NSV_USER_n=>NSV_USER +USE MODD_PARAMETERS +USE MODD_PASPOL +USE MODD_SALT USE MODD_VAR_ll, ONLY: NPROC -! +USE MODD_VISCOSITY + +USE MODE_MSG +USE MODE_POS + +USE MODI_INI_NSV +USE MODI_TEST_NAM_VAR + +USE MODN_2D_FRC +USE MODN_ADV_n ! The final filling of these modules for the model n is USE MODN_BACKUP +USE MODN_BLANK_n +USE MODN_BLOWSNOW +USE MODN_BLOWSNOW_n USE MODN_BUDGET -USE MODN_LES +USE MODN_CH_MNHC_n +USE MODN_CH_ORILAM +USE MODN_CH_SOLVER_n +USE MODN_CONDSAMP USE MODN_CONF +USE MODN_CONF_n USE MODN_CONFZ -USE MODN_FRC +USE MODN_DRAGBLDG_n +USE MODN_DRAG_n +USE MODN_DRAGTREE_n +USE MODN_DUST USE MODN_DYN -USE MODN_NESTING -USE MODN_OUTPUT -USE MODN_CONF_n -USE MODN_LBC_n ! routine is used for each nested model. This has been done USE MODN_DYN_n ! to avoid the duplication of this routine for each model. -USE MODN_ADV_n ! The final filling of these modules for the model n is -USE MODN_PARAM_n ! realized in subroutine ini_model n -USE MODN_PARAM_RAD_n -USE MODN_PARAM_ECRAD_n -USE MODN_PARAM_KAFR_n -USE MODN_PARAM_MFSHALL_n -USE MODN_PARAM_ICE +USE MODN_ELEC +USE MODN_EOL +USE MODN_EOL_ADNR +USE MODN_EOL_ALM +#ifdef MNH_FOREFIRE +USE MODN_FOREFIRE +#endif +USE MODN_FRC +USE MODN_IBM_PARAM_n +USE MODN_LATZ_EDFLX +USE MODN_LBC_n ! routine is used for each nested model. This has been done +USE MODN_LES USE MODN_LUNIT_n +USE MODN_MEAN +USE MODN_NESTING USE MODN_NUDGING_n -USE MODN_TURB_n -USE MODN_DRAG_n -USE MODN_BLANK -USE MODN_CH_MNHC_n -USE MODN_CH_SOLVER_n -USE MODN_PARAM_C2R2, ONLY : EPARAM_CCN=>HPARAM_CCN, EINI_CCN=>HINI_CCN, & - WNUC=>XNUC, WALPHAC=>XALPHAC, NAM_PARAM_C2R2 +USE MODN_OUTPUT USE MODN_PARAM_C1R3, ONLY : NAM_PARAM_C1R3, CPRISTINE_ICE_C1R3, & CHEVRIMED_ICE_C1R3 +USE MODN_PARAM_C2R2, ONLY : EPARAM_CCN=>HPARAM_CCN, EINI_CCN=>HINI_CCN, & + WNUC=>XNUC, WALPHAC=>XALPHAC, NAM_PARAM_C2R2 +USE MODN_PARAM_ECRAD_n +USE MODN_PARAM_ICE +USE MODN_PARAM_KAFR_n USE MODN_PARAM_LIMA, ONLY : FINI_CCN=>HINI_CCN,NAM_PARAM_LIMA,NMOD_CCN,LSCAV, & CPRISTINE_ICE_LIMA, CHEVRIMED_ICE_LIMA, NMOD_IFN, & - LCOLD, LACTI, LNUCL, XALPHAC, XNUC, LMEYERS, LHAIL -USE MODN_ELEC + LCOLD, LACTI, LNUCL, XALPHAC, XNUC, LMEYERS, LHAIL,& + LPTSPLIT +USE MODN_PARAM_MFSHALL_n +USE MODN_PARAM_n ! realized in subroutine ini_model n +USE MODN_PARAM_RAD_n +USE MODN_PASPOL +USE MODN_RECYCL_PARAM_n +USE MODN_SALT USE MODN_SERIES -USE MODN_SERIES_n -USE MODN_TURB_CLOUD +USE MODN_SERIES_n +USE MODN_STATION_n USE MODN_TURB -USE MODN_MEAN -USE MODN_DRAGTREE -USE MODN_LATZ_EDFLX -! -USE MODD_NSV,NSV_USER_n=>NSV_USER -USE MODD_DYN -USE MODD_DYN_n, ONLY : LHORELAX_SVLIMA -USE MODD_GET_n -USE MODD_GR_FIELD_n -! -USE MODE_POS -USE MODE_MSG -! -USE MODI_TEST_NAM_VAR -USE MODI_INI_NSV -USE MODI_CH_INIT_SCHEME_n -USE MODN_CH_ORILAM -USE MODD_CH_AEROSOL -USE MODD_DUST -USE MODD_SALT -USE MODD_PASPOL -#ifdef MNH_FOREFIRE -USE MODD_FOREFIRE -USE MODN_FOREFIRE -#endif -USE MODD_CONDSAMP -USE MODD_BLOWSNOW -USE MODN_DUST -USE MODN_SALT -USE MODD_CH_M9_n, ONLY : NEQ -USE MODN_PASPOL -USE MODN_CONDSAMP -USE MODN_BLOWSNOW -USE MODN_BLOWSNOW_n -USE MODN_2D_FRC +USE MODN_TURB_CLOUD +USE MODN_TURB_n USE MODN_VISCOSITY -USE MODD_VISCOSITY -USE MODD_DRAG_n -! + IMPLICIT NONE ! !* 0.1 declarations of arguments @@ -444,8 +463,6 @@ INTEGER :: JS,JCI,JI,JSV ! Loop indexes LOGICAL :: GRELAX LOGICAL :: GFOUND ! Return code when searching namelist ! -INTEGER :: IMOMENTS, JMODE, IMODEIDX, JMOM, JSV_NAME, JMOD, I -! !------------------------------------------------------------------------------- ! !* 1. READ EXSEG FILE @@ -461,6 +478,8 @@ CCPLFILE(:)=" " CALL INIT_NAM_CONFN CALL INIT_NAM_DYNN CALL INIT_NAM_ADVN +CALL INIT_NAM_DRAGTREEN +CALL INIT_NAM_DRAGBLDGN CALL INIT_NAM_PARAMN CALL INIT_NAM_PARAM_RADN #ifdef MNH_ECRAD @@ -471,11 +490,15 @@ CALL INIT_NAM_PARAM_MFSHALLN CALL INIT_NAM_LBCN CALL INIT_NAM_NUDGINGN CALL INIT_NAM_TURBN +CALL INIT_NAM_BLANKN CALL INIT_NAM_DRAGN +CALL INIT_NAM_IBM_PARAMN +CALL INIT_NAM_RECYCL_PARAMN CALL INIT_NAM_CH_MNHCN CALL INIT_NAM_CH_SOLVERN CALL INIT_NAM_SERIESN CALL INIT_NAM_BLOWSNOWN +CALL INIT_NAM_STATIONn ! WRITE(UNIT=ILUOUT,FMT="(/,'READING THE EXSEG.NAM FILE')") CALL POSNAM(ILUSEG,'NAM_LUNITN',GFOUND,ILUOUT) @@ -506,14 +529,32 @@ CALL POSNAM(ILUSEG,'NAM_TURBN',GFOUND,ILUOUT) IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_TURBn) CALL POSNAM(ILUSEG,'NAM_DRAGN',GFOUND,ILUOUT) IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_DRAGn) +CALL POSNAM(ILUSEG,'NAM_IBM_PARAMN',GFOUND,ILUOUT) +IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_IBM_PARAMn) +CALL POSNAM(ILUSEG,'NAM_RECYCL_PARAMN',GFOUND,ILUOUT) +IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_RECYCL_PARAMn) CALL POSNAM(ILUSEG,'NAM_CH_MNHCN',GFOUND,ILUOUT) IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_CH_MNHCn) CALL POSNAM(ILUSEG,'NAM_CH_SOLVERN',GFOUND,ILUOUT) IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_CH_SOLVERn) CALL POSNAM(ILUSEG,'NAM_SERIESN',GFOUND,ILUOUT) IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_SERIESn) +CALL POSNAM(ILUSEG,'NAM_BLANKN',GFOUND,ILUOUT) +IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_BLANKn) CALL POSNAM(ILUSEG,'NAM_BLOWSNOWN',GFOUND,ILUOUT) IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_BLOWSNOWn) +CALL POSNAM(ILUSEG,'NAM_DRAGTREEN',GFOUND,ILUOUT) +IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_DRAGTREEn) +CALL POSNAM(ILUSEG,'NAM_DRAGBLDGN',GFOUND,ILUOUT) +IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_DRAGBLDGn) +CALL POSNAM(ILUSEG,'NAM_EOL',GFOUND,ILUOUT) +IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_EOL) +CALL POSNAM(ILUSEG,'NAM_EOL_ADNR',GFOUND,ILUOUT) +IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_EOL_ADNR) +CALL POSNAM(ILUSEG,'NAM_EOL_ALM',GFOUND,ILUOUT) +IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_EOL_ALM) +CALL POSNAM(ILUSEG,'NAM_STATIONN',GFOUND,ILUOUT) +IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_STATIONn) ! IF (KMI == 1) THEN WRITE(UNIT=ILUOUT,FMT="(' namelists common to all the models ')") @@ -584,40 +625,182 @@ IF (KMI == 1) THEN END IF CALL POSNAM(ILUSEG,'NAM_BUDGET',GFOUND,ILUOUT) IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_BUDGET) + CALL POSNAM(ILUSEG,'NAM_BU_RU',GFOUND,ILUOUT) - IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_BU_RU) + IF (GFOUND) THEN + IF ( ALLOCATED( CBULIST_RU ) ) THEN + CALL Print_msg( NVERB_WARNING, 'IO', 'READ_EXSEG_n', 'unexpected: CBULIST_RU was already allocated' ) + DEALLOCATE( CBULIST_RU ) + END IF + ALLOCATE( CHARACTER(LEN=NBULISTMAXLEN) :: CBULIST_RU(NBULISTMAXLINES) ) + CBULIST_RU(:) = '' + READ(UNIT=ILUSEG,NML=NAM_BU_RU) + ELSE + ALLOCATE( CHARACTER(LEN=NBULISTMAXLEN) :: CBULIST_RU(0) ) + END IF + CALL POSNAM(ILUSEG,'NAM_BU_RV',GFOUND,ILUOUT) - IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_BU_RV) + IF (GFOUND) THEN + IF ( ALLOCATED( CBULIST_RV ) ) THEN + CALL Print_msg( NVERB_WARNING, 'IO', 'READ_EXSEG_n', 'unexpected: CBULIST_RV was already allocated' ) + DEALLOCATE( CBULIST_RV ) + END IF + ALLOCATE( CHARACTER(LEN=NBULISTMAXLEN) :: CBULIST_RV(NBULISTMAXLINES) ) + CBULIST_RV(:) = '' + READ(UNIT=ILUSEG,NML=NAM_BU_RV) + ELSE + ALLOCATE( CHARACTER(LEN=NBULISTMAXLEN) :: CBULIST_RV(0) ) + END IF + CALL POSNAM(ILUSEG,'NAM_BU_RW',GFOUND,ILUOUT) - IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_BU_RW) + IF (GFOUND) THEN + IF ( ALLOCATED( CBULIST_RW ) ) THEN + CALL Print_msg( NVERB_WARNING, 'IO', 'READ_EXSEG_n', 'unexpected: CBULIST_RW was already allocated' ) + DEALLOCATE( CBULIST_RW ) + END IF + ALLOCATE( CHARACTER(LEN=NBULISTMAXLEN) :: CBULIST_RW(NBULISTMAXLINES) ) + CBULIST_RW(:) = '' + READ(UNIT=ILUSEG,NML=NAM_BU_RW) + ELSE + ALLOCATE( CHARACTER(LEN=NBULISTMAXLEN) :: CBULIST_RW(0) ) + END IF + CALL POSNAM(ILUSEG,'NAM_BU_RTH',GFOUND,ILUOUT) - IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_BU_RTH) + IF (GFOUND) THEN + IF ( ALLOCATED( CBULIST_RTH ) ) THEN + CALL Print_msg( NVERB_WARNING, 'IO', 'READ_EXSEG_n', 'unexpected: CBULIST_RTH was already allocated' ) + DEALLOCATE( CBULIST_RTH ) + END IF + ALLOCATE( CHARACTER(LEN=NBULISTMAXLEN) :: CBULIST_RTH(NBULISTMAXLINES) ) + CBULIST_RTH(:) = '' + READ(UNIT=ILUSEG,NML=NAM_BU_RTH) + ELSE + ALLOCATE( CHARACTER(LEN=NBULISTMAXLEN) :: CBULIST_RTH(0) ) + END IF + CALL POSNAM(ILUSEG,'NAM_BU_RTKE',GFOUND,ILUOUT) - IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_BU_RTKE) + IF (GFOUND) THEN + IF ( ALLOCATED( CBULIST_RTKE ) ) THEN + CALL Print_msg( NVERB_WARNING, 'IO', 'READ_EXSEG_n', 'unexpected: CBULIST_RTKE was already allocated' ) + DEALLOCATE( CBULIST_RTKE ) + END IF + ALLOCATE( CHARACTER(LEN=NBULISTMAXLEN) :: CBULIST_RTKE(NBULISTMAXLINES) ) + CBULIST_RTKE(:) = '' + READ(UNIT=ILUSEG,NML=NAM_BU_RTKE) + ELSE + ALLOCATE( CHARACTER(LEN=NBULISTMAXLEN) :: CBULIST_RTKE(0) ) + END IF + CALL POSNAM(ILUSEG,'NAM_BU_RRV',GFOUND,ILUOUT) - IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_BU_RRV) + IF (GFOUND) THEN + IF ( ALLOCATED( CBULIST_RRV ) ) THEN + CALL Print_msg( NVERB_WARNING, 'IO', 'READ_EXSEG_n', 'unexpected: CBULIST_RRV was already allocated' ) + DEALLOCATE( CBULIST_RRV ) + END IF + ALLOCATE( CHARACTER(LEN=NBULISTMAXLEN) :: CBULIST_RRV(NBULISTMAXLINES) ) + CBULIST_RRV(:) = '' + READ(UNIT=ILUSEG,NML=NAM_BU_RRV) + ELSE + ALLOCATE( CHARACTER(LEN=NBULISTMAXLEN) :: CBULIST_RRV(0) ) + END IF + CALL POSNAM(ILUSEG,'NAM_BU_RRC',GFOUND,ILUOUT) - IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_BU_RRC) + IF (GFOUND) THEN + IF ( ALLOCATED( CBULIST_RRC ) ) THEN + CALL Print_msg( NVERB_WARNING, 'IO', 'READ_EXSEG_n', 'unexpected: CBULIST_RRC was already allocated' ) + DEALLOCATE( CBULIST_RRC ) + END IF + ALLOCATE( CHARACTER(LEN=NBULISTMAXLEN) :: CBULIST_RRC(NBULISTMAXLINES) ) + CBULIST_RRC(:) = '' + READ(UNIT=ILUSEG,NML=NAM_BU_RRC) + ELSE + ALLOCATE( CHARACTER(LEN=NBULISTMAXLEN) :: CBULIST_RRC(0) ) + END IF + CALL POSNAM(ILUSEG,'NAM_BU_RRR',GFOUND,ILUOUT) - IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_BU_RRR) + IF (GFOUND) THEN + IF ( ALLOCATED( CBULIST_RRR ) ) THEN + CALL Print_msg( NVERB_WARNING, 'IO', 'READ_EXSEG_n', 'unexpected: CBULIST_RRR was already allocated' ) + DEALLOCATE( CBULIST_RRR ) + END IF + ALLOCATE( CHARACTER(LEN=NBULISTMAXLEN) :: CBULIST_RRR(NBULISTMAXLINES) ) + CBULIST_RRR(:) = '' + READ(UNIT=ILUSEG,NML=NAM_BU_RRR) + ELSE + ALLOCATE( CHARACTER(LEN=NBULISTMAXLEN) :: CBULIST_RRR(0) ) + END IF + CALL POSNAM(ILUSEG,'NAM_BU_RRI',GFOUND,ILUOUT) - IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_BU_RRI) + IF (GFOUND) THEN + IF ( ALLOCATED( CBULIST_RRI ) ) THEN + CALL Print_msg( NVERB_WARNING, 'IO', 'READ_EXSEG_n', 'unexpected: CBULIST_RRI was already allocated' ) + DEALLOCATE( CBULIST_RRI ) + END IF + ALLOCATE( CHARACTER(LEN=NBULISTMAXLEN) :: CBULIST_RRI(NBULISTMAXLINES) ) + CBULIST_RRI(:) = '' + READ(UNIT=ILUSEG,NML=NAM_BU_RRI) + ELSE + ALLOCATE( CHARACTER(LEN=NBULISTMAXLEN) :: CBULIST_RRI(0) ) + END IF + CALL POSNAM(ILUSEG,'NAM_BU_RRS',GFOUND,ILUOUT) - IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_BU_RRS) + IF (GFOUND) THEN + IF ( ALLOCATED( CBULIST_RRS ) ) THEN + CALL Print_msg( NVERB_WARNING, 'IO', 'READ_EXSEG_n', 'unexpected: CBULIST_RRS was already allocated' ) + DEALLOCATE( CBULIST_RRS ) + END IF + ALLOCATE( CHARACTER(LEN=NBULISTMAXLEN) :: CBULIST_RRS(NBULISTMAXLINES) ) + CBULIST_RRS(:) = '' + READ(UNIT=ILUSEG,NML=NAM_BU_RRS) + ELSE + ALLOCATE( CHARACTER(LEN=NBULISTMAXLEN) :: CBULIST_RRS(0) ) + END IF + CALL POSNAM(ILUSEG,'NAM_BU_RRG',GFOUND,ILUOUT) - IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_BU_RRG) + IF (GFOUND) THEN + IF ( ALLOCATED( CBULIST_RRG ) ) THEN + CALL Print_msg( NVERB_WARNING, 'IO', 'READ_EXSEG_n', 'unexpected: CBULIST_RRG was already allocated' ) + DEALLOCATE( CBULIST_RRG ) + END IF + ALLOCATE( CHARACTER(LEN=NBULISTMAXLEN) :: CBULIST_RRG(NBULISTMAXLINES) ) + CBULIST_RRG(:) = '' + READ(UNIT=ILUSEG,NML=NAM_BU_RRG) + ELSE + ALLOCATE( CHARACTER(LEN=NBULISTMAXLEN) :: CBULIST_RRG(0) ) + END IF + CALL POSNAM(ILUSEG,'NAM_BU_RRH',GFOUND,ILUOUT) - IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_BU_RRH) + IF (GFOUND) THEN + IF ( ALLOCATED( CBULIST_RRH ) ) THEN + CALL Print_msg( NVERB_WARNING, 'IO', 'READ_EXSEG_n', 'unexpected: CBULIST_RRH was already allocated' ) + DEALLOCATE( CBULIST_RRH ) + END IF + ALLOCATE( CHARACTER(LEN=NBULISTMAXLEN) :: CBULIST_RRH(NBULISTMAXLINES) ) + CBULIST_RRH(:) = '' + READ(UNIT=ILUSEG,NML=NAM_BU_RRH) + ELSE + ALLOCATE( CHARACTER(LEN=NBULISTMAXLEN) :: CBULIST_RRH(0) ) + END IF + CALL POSNAM(ILUSEG,'NAM_BU_RSV',GFOUND,ILUOUT) - IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_BU_RSV) + IF (GFOUND) THEN + IF ( ALLOCATED( CBULIST_RSV ) ) THEN + CALL Print_msg( NVERB_WARNING, 'IO', 'READ_EXSEG_n', 'unexpected: CBULIST_RSV was already allocated' ) + DEALLOCATE( CBULIST_RSV ) + END IF + ALLOCATE( CHARACTER(LEN=NBULISTMAXLEN) :: CBULIST_RSV(NBULISTMAXLINES) ) + CBULIST_RSV(:) = '' + READ(UNIT=ILUSEG,NML=NAM_BU_RSV) + ELSE + ALLOCATE( CHARACTER(LEN=NBULISTMAXLEN) :: CBULIST_RSV(0) ) + END IF + CALL POSNAM(ILUSEG,'NAM_LES',GFOUND,ILUOUT) IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_LES) CALL POSNAM(ILUSEG,'NAM_MEAN',GFOUND,ILUOUT) IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_MEAN) CALL POSNAM(ILUSEG,'NAM_PDF',GFOUND,ILUOUT) IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_PDF) - CALL POSNAM(ILUSEG,'NAM_BLANK',GFOUND,ILUOUT) - IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_BLANK) CALL POSNAM(ILUSEG,'NAM_FRC',GFOUND,ILUOUT) IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_FRC) CALL POSNAM(ILUSEG,'NAM_PARAM_ICE',GFOUND,ILUOUT) @@ -650,8 +833,6 @@ IF (KMI == 1) THEN #endif CALL POSNAM(ILUSEG,'NAM_CONDSAMP',GFOUND,ILUOUT) IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_CONDSAMP) - CALL POSNAM(ILUSEG,'NAM_DRAGTREE',GFOUND,ILUOUT) - IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_DRAGTREE) CALL POSNAM(ILUSEG,'NAM_2D_FRC',GFOUND,ILUOUT) IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_2D_FRC) CALL POSNAM(ILUSEG,'NAM_LATZ_EDFLX',GFOUND) @@ -704,9 +885,13 @@ CALL TEST_NAM_VAR(ILUOUT,'CLBCY(1)',CLBCY(1),'CYCL','WALL','OPEN') CALL TEST_NAM_VAR(ILUOUT,'CLBCY(2)',CLBCY(2),'CYCL','WALL','OPEN') ! CALL TEST_NAM_VAR(ILUOUT,'CTURBDIM',CTURBDIM,'1DIM','3DIM') -CALL TEST_NAM_VAR(ILUOUT,'CTURBLEN',CTURBLEN,'DELT','BL89','RM17','DEAR','BLKR') +CALL TEST_NAM_VAR(ILUOUT,'CTURBLEN',CTURBLEN,'DELT','BL89','RM17','DEAR','BLKR','ADAP') CALL TEST_NAM_VAR(ILUOUT,'CTOM',CTOM,'NONE','TM06') -CALL TEST_NAM_VAR(ILUOUT,'CSUBG_AUCV',CSUBG_AUCV,'NONE','CLFR','SIGM','PDF') +CALL TEST_NAM_VAR(ILUOUT,'CSUBG_AUCV',CSUBG_AUCV,'NONE','CLFR','SIGM','PDF','ADJU') +CALL TEST_NAM_VAR(ILUOUT,'CSUBG_AUCV_RI',CSUBG_AUCV_RI,'NONE','CLFR','ADJU') +CALL TEST_NAM_VAR(ILUOUT,'CCONDENS',CCONDENS,'CB02','GAUS') +CALL TEST_NAM_VAR(ILUOUT,'CLAMBDA3',CLAMBDA3,'CB','NONE') +CALL TEST_NAM_VAR(ILUOUT,'CSUBG_MF_PDF',CSUBG_MF_PDF,'NONE','TRIANGLE') ! CALL TEST_NAM_VAR(ILUOUT,'CCH_TDISCRETIZATION',CCH_TDISCRETIZATION, & 'SPLIT ','CENTER ','LAGGED ') @@ -725,8 +910,7 @@ CALL TEST_NAM_VAR(ILUOUT,'CTURBLEN_CLOUD',CTURBLEN_CLOUD,'NONE','DEAR','DELT','B ! ! The test on the mass flux scheme for shallow convection ! -CALL TEST_NAM_VAR(ILUOUT,'CMF_UPDRAFT',CMF_UPDRAFT,'NONE','EDKF','RHCJ',& - 'HRIO','BOUT') +CALL TEST_NAM_VAR(ILUOUT,'CMF_UPDRAFT',CMF_UPDRAFT,'NONE','EDKF','RHCJ') CALL TEST_NAM_VAR(ILUOUT,'CMF_CLOUD',CMF_CLOUD,'NONE','STAT','DIRE') ! ! The test on the CSOLVER name is made elsewhere @@ -1073,16 +1257,23 @@ SELECT CASE ( CCLOUD ) LUSERH=LHAIL END IF ! - IF (LSUBG_COND .AND. LCOLD) THEN - WRITE(UNIT=ILUOUT,FMT=9003) KMI - WRITE(UNIT=ILUOUT,FMT=*) 'YOU WANT TO USE BOTH THE SIMPLE MIXED PHASE' - WRITE(UNIT=ILUOUT,FMT=*) 'MICROPHYS. SCHEME AND THE SUBGRID COND. SCHEME.' - WRITE(UNIT=ILUOUT,FMT=*) 'THIS IS NOT YET AVAILABLE. SET LSUBG_COND ' - WRITE(UNIT=ILUOUT,FMT=*) 'TO FALSE OR CCLOUD TO "REVE", "KESS" ' - !callabortstop - CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_EXSEG_n','') - END IF -! +!!$ IF (LSUBG_COND .AND. LCOLD) THEN +!!$ WRITE(UNIT=ILUOUT,FMT=9003) KMI +!!$ WRITE(UNIT=ILUOUT,FMT=*) 'YOU WANT TO USE BOTH THE SIMPLE MIXED PHASE' +!!$ WRITE(UNIT=ILUOUT,FMT=*) 'MICROPHYS. SCHEME AND THE SUBGRID COND. SCHEME.' +!!$ WRITE(UNIT=ILUOUT,FMT=*) 'THIS IS NOT YET AVAILABLE. SET LSUBG_COND ' +!!$ WRITE(UNIT=ILUOUT,FMT=*) 'TO FALSE OR CCLOUD TO "REVE", "KESS" ' +!!$ !callabortstop +!!$ CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_EXSEG_n','') +!!$ END IF +! + IF (CCLOUD == 'LIMA' .AND. LSUBG_COND .AND. (.NOT. LPTSPLIT)) THEN + WRITE(UNIT=ILUOUT,FMT=9001) KMI + WRITE(UNIT=ILUOUT,FMT=*) 'YOU MUST USE LPTSPLIT=T with CCLOUD=LIMA' + WRITE(UNIT=ILUOUT,FMT=*) 'AND LSUBG_COND ' + CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_EXSEG_n','use LPTSPLIT=T with LIMA and LSUBG_COND=T') + END IF +! IF ( XALPHAC /= 3.0 .OR. XNUC /= 2.0) THEN WRITE(UNIT=ILUOUT,FMT=9001) KMI WRITE(UNIT=ILUOUT,FMT=*) 'IT IS ADVISED TO USE XALPHAC=3. and XNUC=2.' @@ -1133,6 +1324,8 @@ LUSETKE(KMI) = (CTURB /= 'NONE') ! !* 2.3 Chemical and NSV_* variables initializations ! +CALL UPDATE_NAM_IBM_PARAMN +CALL UPDATE_NAM_RECYCL_PARAMN CALL UPDATE_NAM_PARAMN CALL UPDATE_NAM_DYNN CALL UPDATE_NAM_CONFN @@ -1209,12 +1402,6 @@ IF ((LUSECHIC).AND.(LCH_RET_ICE)) THEN WRITE(UNIT=ILUOUT,FMT=*) 'TO FALSE IE NO CHEMICAL SPECIES IN ICE' ENDIF ! -IF (LUSECHEM) THEN - CALL CH_INIT_SCHEME_n(KMI,LUSECHAQ,LUSECHIC,LCH_PH,ILUOUT,NVERB) - IF (LORILAM) CALL CH_AER_INIT_SOA(ILUOUT, NVERB) -END IF -! - CALL UPDATE_NAM_CH_MNHCN CALL INI_NSV(KMI) ! @@ -1420,7 +1607,7 @@ ELSE END IF END IF ! -IF(CTURBLEN=='RM17') THEN +IF(CTURBLEN=='RM17' .OR. CTURBLEN=='ADAP') THEN XCEDIS=0.34 ELSE XCEDIS=0.84 @@ -1703,44 +1890,7 @@ IF (LDUST) THEN !callabortstop CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_EXSEG_n','') END IF - IF(.NOT.ALLOCATED(CDUSTNAMES)) THEN - IMOMENTS = (NSV_DSTEND - NSV_DSTBEG +1 )/NMODE_DST - ALLOCATE(CDUSTNAMES(IMOMENTS*NMODE_DST)) - !Loop on all dust modes - IF (IMOMENTS == 1) THEN - DO JMODE=1,NMODE_DST - IMODEIDX=JPDUSTORDER(JMODE) - JSV_NAME = (IMODEIDX - 1)*3 + 2 - CDUSTNAMES(JMODE) = YPDUST_INI(JSV_NAME) - END DO - ELSE - DO JMODE=1,NMODE_DST - !Find which mode we are dealing with - IMODEIDX=JPDUSTORDER(JMODE) - DO JMOM=1,IMOMENTS - !Find which number this is of the list of scalars - JSV = (JMODE-1)*IMOMENTS + JMOM - !Find what name this corresponds to, always 3 moments assumed in YPDUST_INI - JSV_NAME = (IMODEIDX - 1)*3 + JMOM - !Get the right CDUSTNAMES which should follow the list of scalars transported in XSVM/XSVT - CDUSTNAMES(JSV) = YPDUST_INI(JSV_NAME) - ENDDO ! Loop on moments - ENDDO ! Loop on dust modes - END IF - END IF - ! Initialization of deposition scheme - IF (LDEPOS_DST(KMI)) THEN - IF(.NOT.ALLOCATED(CDEDSTNAMES)) THEN - ALLOCATE(CDEDSTNAMES(NMODE_DST*2)) - DO JMODE=1,NMODE_DST - IMODEIDX=JPDUSTORDER(JMODE) - CDEDSTNAMES(JMODE) = YPDEDST_INI(IMODEIDX) - CDEDSTNAMES(NMODE_DST+JMODE) = YPDEDST_INI(NMODE_DST+IMODEIDX) - ENDDO - ENDIF - ENDIF - -END IF +END IF ! ! Sea Salt case ! @@ -1786,42 +1936,6 @@ IF (LSALT) THEN !callabortstop CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_EXSEG_n','') END IF - IF(.NOT.ALLOCATED(CSALTNAMES)) THEN - IMOMENTS = (NSV_SLTEND - NSV_SLTBEG +1 )/NMODE_SLT - ALLOCATE(CSALTNAMES(IMOMENTS*NMODE_SLT)) - !Loop on all dust modes - IF (IMOMENTS == 1) THEN - DO JMODE=1,NMODE_SLT - IMODEIDX=JPSALTORDER(JMODE) - JSV_NAME = (IMODEIDX - 1)*3 + 2 - CSALTNAMES(JMODE) = YPSALT_INI(JSV_NAME) - END DO - ELSE - DO JMODE=1,NMODE_SLT - !Find which mode we are dealing with - IMODEIDX=JPSALTORDER(JMODE) - DO JMOM=1,IMOMENTS - !Find which number this is of the list of scalars - JSV = (JMODE-1)*IMOMENTS + JMOM - !Find what name this corresponds to, always 3 moments assumed in YPSALT_INI - JSV_NAME = (IMODEIDX - 1)*3 + JMOM - !Get the right CSALTNAMES which should follow the list of scalars transported in XSVM/XSVT - CSALTNAMES(JSV) = YPSALT_INI(JSV_NAME) - ENDDO ! Loop on moments - ENDDO ! Loop on dust modes - END IF - END IF - ! Initialization of deposition scheme - IF (LDEPOS_SLT(KMI)) THEN - IF(.NOT.ALLOCATED(CDESLTNAMES)) THEN - ALLOCATE(CDESLTNAMES(NMODE_SLT*2)) - DO JMODE=1,NMODE_SLT - IMODEIDX=JPSALTORDER(JMODE) - CDESLTNAMES(JMODE) = YPDESLT_INI(IMODEIDX) - CDESLTNAMES(NMODE_SLT+JMODE) = YPDESLT_INI(NMODE_SLT+IMODEIDX) - ENDDO - ENDIF - ENDIF END IF ! ! Orilam SV case @@ -1860,13 +1974,6 @@ IF (LORILAM) THEN CGETSVT(NSV_AERDEPBEG:NSV_AERDEPEND)='INIT' END IF END IF -! Initialization of deposition scheme - IF (LDEPOS_AER(KMI)) THEN - IF(.NOT.ALLOCATED(CDEAERNAMES)) THEN - ALLOCATE(CDEAERNAMES(JPMODE*2)) - CDEAERNAMES(:) = YPDEAER_INI(:) - ENDIF - ENDIF END IF ! ! Lagrangian variables @@ -1964,13 +2071,6 @@ IF (LBLOWSNOW) THEN & "THE BLOWING SNOW VARIABLES HAVE BEEN INITIALIZED TO ZERO ")') CGETSVT(NSV_SNWBEG:NSV_SNWEND)='INIT' END IF - IF(.NOT.ALLOCATED(CSNOWNAMES)) THEN - IMOMENTS = (NSV_SNWEND - NSV_SNWBEG +1 ) - ALLOCATE(CSNOWNAMES(IMOMENTS)) - DO JMOM=1,IMOMENTS - CSNOWNAMES(JMOM) = YPSNOW_INI(JMOM) - ENDDO ! Loop on moments - END IF END IF ! ! @@ -2264,15 +2364,24 @@ IF ( LFORCING ) THEN WRITE(ILUOUT,FMT=*) 'YOU CHOSE A TEMPERATURE AND HUMIDITY RELAXATION' WRITE(ILUOUT,FMT=*) 'TOGETHER WITH TENDENCY OR GEOSTROPHIC FORCING' WRITE(ILUOUT,FMT=*) & - 'YOU MIGHT CHECK YOUR SWITCHES: LRELAX_THRV_FRC, LTEND_THRV_FRC, AND' + 'YOU MIGHT CHECK YOUR SWITCHES: LRELAX_THRV_FRC, LTEND_THRV_FRC, AND' WRITE(ILUOUT,FMT=*) 'LGEOST_TH_FRC' END IF ! - IF ( LRELAX_UV_FRC .AND. LGEOST_UV_FRC ) THEN + IF ( LRELAX_UV_FRC .AND. LRELAX_UVMEAN_FRC) THEN + WRITE(UNIT=ILUOUT,FMT=9003) KMI + WRITE(ILUOUT,FMT=*) 'YOU MUST CHOOSE BETWEEN A RELAXATION APPLIED TO' + WRITE(ILUOUT,FMT=*) 'THE 3D FULL WIND FIELD (LRELAX_UV_FRC) OR' + WRITE(ILUOUT,FMT=*) 'THE HORIZONTAL MEAN WIND (LRELAX_UVMEAN_FRC)' + !callabortstop + CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_EXSEG_n','') + END IF +! + IF ( (LRELAX_UV_FRC .OR. LRELAX_UVMEAN_FRC) .AND. LGEOST_UV_FRC ) THEN WRITE(UNIT=ILUOUT,FMT=9003) KMI WRITE(ILUOUT,FMT=*) 'YOU MUST NOT USE A WIND RELAXATION' WRITE(ILUOUT,FMT=*) 'TOGETHER WITH A GEOSTROPHIC FORCING' - WRITE(ILUOUT,FMT=*) 'CHECK SWITCHES: LRELAX_UV_FRC, LGEOST_UV_FRC' + WRITE(ILUOUT,FMT=*) 'CHECK SWITCHES: LRELAX_UV_FRC, LRELAX_UVMEAN_FRC, LGEOST_UV_FRC' !callabortstop CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_EXSEG_n','') END IF @@ -2584,7 +2693,31 @@ IF ((LHORELAX_UVWTH .OR. LHORELAX_SVPP .OR. & !callabortstop CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_EXSEG_n','') END IF -! +! +IF ((LHORELAX_UVWTH .OR. LHORELAX_SVPP .OR. & + LHORELAX_SVCS .OR. & +#ifdef MNH_FOREFIRE + LHORELAX_SVFF .OR. & +#endif + LHORELAX_SVC2R2 .OR. LHORELAX_SVC1R3 .OR. & + LHORELAX_SVLIMA .OR. & + LHORELAX_SVELEC .OR. LHORELAX_SVCHEM .OR. & + LHORELAX_SVLG .OR. ANY(LHORELAX_SV) .OR. & + LHORELAX_RV .OR. LHORELAX_RC .OR. & + LHORELAX_RR .OR. LHORELAX_RI .OR. & + LHORELAX_RG .OR. LHORELAX_RS .OR. & + LHORELAX_RH .OR. LHORELAX_TKE.OR. & + LHORELAX_SVCHIC ) & + .AND. (CLBCX(1)=='CYCL'.OR.CLBCX(2)=='CYCL' & + .OR.CLBCY(1)=='CYCL'.OR.CLBCY(2)=='CYCL')) THEN + WRITE(UNIT=ILUOUT,FMT=9003) KMI + WRITE(ILUOUT,FMT=*) 'YOU WANT TO USE THE HORIZONTAL RELAXATION ' + WRITE(ILUOUT,FMT=*) 'FOR CYCLIC CLBCX OR CLBCY VALUES' + WRITE(ILUOUT,FMT=*) 'CHANGE LHORELAX TO FALSE' + !callabortstop + CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_EXSEG_n','') +END IF +! IF (KMI==1) THEN GRELAX = .NOT.(OUSERV) .AND. LUSERV .AND. LHORELAX_RV ELSE @@ -2816,6 +2949,8 @@ END IF ! CALL UPDATE_NAM_LUNITN CALL UPDATE_NAM_CONFN +CALL UPDATE_NAM_DRAGTREEN +CALL UPDATE_NAM_DRAGBLDGN CALL UPDATE_NAM_DYNN CALL UPDATE_NAM_ADVN CALL UPDATE_NAM_PARAMN @@ -2828,10 +2963,12 @@ CALL UPDATE_NAM_PARAM_MFSHALLN CALL UPDATE_NAM_LBCN CALL UPDATE_NAM_NUDGINGN CALL UPDATE_NAM_TURBN +CALL UPDATE_NAM_BLANKN CALL UPDATE_NAM_CH_MNHCN CALL UPDATE_NAM_CH_SOLVERN CALL UPDATE_NAM_SERIESN CALL UPDATE_NAM_BLOWSNOWN +CALL UPDATE_NAM_STATIONn !------------------------------------------------------------------------------- WRITE(UNIT=ILUOUT,FMT='(/)') !------------------------------------------------------------------------------- diff --git a/src/ZSOLVER/spectre.f90 b/src/ZSOLVER/spectre.f90 deleted file mode 100644 index 545c60629410ef29c5a16bfd5289f96038c41e28..0000000000000000000000000000000000000000 --- a/src/ZSOLVER/spectre.f90 +++ /dev/null @@ -1,217 +0,0 @@ -!MNH_LIC Copyright 2011-2019 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. -!----------------------------------------------------------------- -! ######spl - PROGRAM SPECTRE -! ############ -! -!!**** -!! -!! PURPOSE -!! ------- -!! compute energy spectra from a MESONH file -!! -!! -!! -!! -!! Modifications: -!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O -!------------------------------------------------------------------------------- -! -!* 0. DECLARATIONS -! ------------ -! -! -USE MODD_CONF -USE MODD_IO, ONLY: NIO_VERB,NVERB_DEBUG,TFILEDATA -USE MODD_LUNIT -USE MODD_LUNIT_n -USE MODD_TIME_n -USE MODD_DIM_ll -USE MODD_SPECTRE -! -USE MODI_SPECTRE_MESONH -USE MODI_SPECTRE_AROME -! -USE MODE_MSG -USE MODE_POS -USE MODE_IO, only: IO_Config_set, IO_Init -USE MODE_IO_FILE, only: IO_File_close, IO_File_open -USE MODE_IO_MANAGE_STRUCT, only: IO_File_add2list, IO_Filelist_print -use mode_init_ll, only: END_PARA_ll -USE MODE_MODELN_HANDLER -!USE MODD_TYPE_DATE -USE MODI_VERSION -! -USE MODN_CONFZ -USE MODN_CONFIO, ONLY : NAM_CONFIO -! -IMPLICIT NONE -! -!* 0.1 declarations of local variables -! -CHARACTER (LEN=28), DIMENSION(1) :: YINIFILE ! names of the INPUT FM-file -CHARACTER (LEN=50) :: YOUTFILE ! names of the OUTPUT FM-file -INTEGER :: IRESP ! return code in FM routines -INTEGER :: ILUOUT0 ! Logical unit number for the output listing -INTEGER :: ILUNAM ! Logical unit numbers for the namelist file - ! and for output_listing file -LOGICAL :: GFOUND ! Return code when searching namelist -! -INTEGER :: IINFO_ll ! return code for _ll routines -! -REAL,DIMENSION(:,:,:),ALLOCATABLE:: ZWORK ! work array -REAL,DIMENSION(:,:,:),ALLOCATABLE:: ZWORKAROME ! work array -INTEGER :: NI,NJ,NK -REAL ::XDELTAX,XDELTAY -TYPE(TFILEDATA),POINTER :: TZNMLFILE => NULL() -! -NAMELIST/NAM_SPECTRE/ LSPECTRE_U,LSPECTRE_V,LSPECTRE_W,LSPECTRE_TH,LSPECTRE_RV,& - LSPECTRE_LSU,LSPECTRE_LSV,LSPECTRE_LSW,LSPECTRE_LSTH,LSPECTRE_LSRV,LSMOOTH -! -NAMELIST/NAM_SPECTRE_FILE/ YINIFILE,CTYPEFILE,YOUTFILE,LSTAT -NAMELIST/NAM_ZOOM_SPECTRE/ LZOOM,NITOT,NJTOT,NXDEB,NYDEB -NAMELIST/NAM_DOMAIN_AROME/ NI,NJ,NK,XDELTAX,XDELTAY -! -!------------------------------------------------------------------------------- -! -!* 0.0 Initializations -! --------------- -! -! -CALL GOTO_MODEL(1) -! -CALL VERSION -CPROGRAM='SPEC ' -! -CALL IO_Init() -! -! initialization -YINIFILE(:) = ' ' -CTYPEFILE = 'MESONH' -LSPECTRE_U = .FALSE. -LSPECTRE_V = .FALSE. -LSPECTRE_W = .FALSE. -LSPECTRE_TH = .FALSE. -LSPECTRE_RV = .FALSE. -LSPECTRE_LSU = .FALSE. -LSPECTRE_LSV = .FALSE. -LSPECTRE_LSW = .FALSE. -LSPECTRE_LSTH = .FALSE. -LSPECTRE_LSRV = .FALSE. -LSMOOTH = .FALSE. -LZOOM = .FALSE. -YOUTFILE = ' ' -LSTAT = .FALSE. -NI=750 -NJ=720 -NK=60 -XDELTAX=2500. -XDELTAY=2500. -! -! -!------------------------------------------------------------------------------- -! -!* 1.0 Namelist reading -! ---------------- -! -PRINT*, ' ' -PRINT*, '*********************************************************************' -PRINT*, '*********************************************************************' -PRINT*, ' ' -! -CALL IO_File_add2list(TZNMLFILE,'SPEC1.nam','NML','READ') -CALL IO_File_open(TZNMLFILE) -ILUNAM = TZNMLFILE%NLU -! -PRINT*, 'READ THE SPEC1.NAM FILE' -! -CALL POSNAM(ILUNAM,'NAM_SPECTRE',GFOUND) -IF (GFOUND) THEN - READ(UNIT=ILUNAM,NML=NAM_SPECTRE) - PRINT*, ' namelist NAM_SPECTRE read' -END IF -! -! -CALL POSNAM(ILUNAM,'NAM_SPECTRE_FILE',GFOUND) -IF (GFOUND) THEN - READ(UNIT=ILUNAM,NML=NAM_SPECTRE_FILE) - PRINT*, ' namelist NAM_SPECTRE_FILE read' -END IF -! -CALL POSNAM(ILUNAM,'NAM_ZOOM_SPECTRE',GFOUND) -IF (GFOUND) THEN - READ(UNIT=ILUNAM,NML=NAM_ZOOM_SPECTRE) - PRINT*, ' namelist NAM_ZOOM_SPECTRE read' -END IF -! -CALL POSNAM(ILUNAM,'NAM_DOMAIN_AROME',GFOUND) -IF (GFOUND) THEN - READ(UNIT=ILUNAM,NML=NAM_DOMAIN_AROME) - PRINT*, ' namelist NAM_DOMAIN_AROME read' -END IF -! -CALL POSNAM(ILUNAM,'NAM_CONFZ',GFOUND) -IF (GFOUND) THEN - READ(UNIT=ILUNAM,NML=NAM_CONFZ) - PRINT*, ' namelist NAM_CONFZ read' -END IF -! -CALL POSNAM(ILUNAM,'NAM_CONFIO',GFOUND) -IF (GFOUND) THEN - READ(UNIT=ILUNAM,NML=NAM_CONFIO) - PRINT*, ' namelist NAM_CONFIO read' -END IF -CALL IO_Config_set() -! -CALL IO_File_close(TZNMLFILE) -! -CINIFILE = YINIFILE(1) -! -!------------------------------------------------------------------------------- -! -!* 2.0 file -! ----------- -! -IF ( LEN_TRIM(CINIFILE)==0 ) THEN - CALL PRINT_MSG(NVERB_FATAL,'GEN','SPECTRE','LEN_TRIM(CINIFILE)==0') -ENDIF -! -IF ( LEN_TRIM(YOUTFILE)==0 ) THEN - WRITE(YOUTFILE,FMT='(A,A)') "spectra_",TRIM(ADJUSTL(CINIFILE)) -ENDIF -! -!------------------------------------------------------------------------------- -! -!* 3.0 Fields initialization and spectra computation -! -IF (CTYPEFILE=='MESONH') THEN - CALL SPECTRE_MESONH(YOUTFILE) - ! - CALL IO_File_close(LUNIT_MODEL(1)%TINIFILE) - IF(NIO_VERB>=NVERB_DEBUG) CALL IO_Filelist_print() - CALL IO_File_close(TLUOUT0) - CALL IO_File_close(TLUOUT) -ELSEIF (CTYPEFILE=='AROME ')THEN - CALL SPECTRE_AROME(CINIFILE,YOUTFILE,XDELTAX,XDELTAY,NI,NJ,NK) -ELSE - print*,"This type of file is not accept for SPECTRE PROGRAM" -ENDIF -! -!------------------------------------------------------------------------------- -! -!* 4. FINALIZE THE PARALLEL SESSION -! ----------------------------- -! -CALL END_PARA_ll(IINFO_ll) -! -PRINT*, ' ' -PRINT*, '****************************************************' -PRINT*, '* EXIT SPECTRE CORRECTLY *' -PRINT*, '****************************************************' -PRINT*, ' ' -!------------------------------------------------------------------------------- -END PROGRAM SPECTRE - diff --git a/src/ZSOLVER/tridiag_thermo.f90 b/src/ZSOLVER/tridiag_thermo.f90 index 7ba28585cf889c05a029487228b15e8dcdc5a908..8eea372766e0639d923a9c90e9674087b25d154a 100644 --- a/src/ZSOLVER/tridiag_thermo.f90 +++ b/src/ZSOLVER/tridiag_thermo.f90 @@ -187,9 +187,13 @@ REAL, DIMENSION(:,:,:), pointer , contiguous :: ZMZM_RHODJ REAL, DIMENSION(:,:,:), pointer , contiguous :: ZA, ZB, ZC REAL, DIMENSION(:,:,:), pointer , contiguous :: ZY ,ZGAM ! RHS of the equation, 3D work array +#ifdef MNH_OPENACC INTEGER :: IZRHODJ_DFDDTDZ_O_DZ2,IZMZM_RHODJ,IZA,IZB,IZC,IZY,IZGAM +#endif REAL, DIMENSION(:,:), pointer , contiguous :: ZBET +#ifdef MNH_OPENACC INTEGER :: IZBET +#endif ! 2D work array INTEGER :: JI,JJ,JK ! loop counter INTEGER :: IKB,IKE ! inner vertical limits diff --git a/src/ZSOLVER/turb.f90 b/src/ZSOLVER/turb.f90 index 9d5f3bfa99c5c9e6b916e628f882e8bb9f8e6d7c..174af509ceeeed91564857d103dce1a96abf37ae 100644 --- a/src/ZSOLVER/turb.f90 +++ b/src/ZSOLVER/turb.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. @@ -22,20 +22,21 @@ implicit none private -public :: turb +public :: Turb contains ! ################################################################# -SUBROUTINE TURB(KKA,KKU,KKL,KMI,KRR,KRRL,KRRI,HLBCX,HLBCY,KSPLIT,KMODEL_CL, & - OCLOSE_OUT,OTURB_FLX,OTURB_DIAG,OSUBG_COND,ORMC01, & + SUBROUTINE TURB(KKA, KKU, KKL, KMI,KRR,KRRL,KRRI,HLBCX,HLBCY, & + KSPLIT,KMODEL_CL, & + OTURB_FLX,OTURB_DIAG,OSUBG_COND,ORMC01, & HTURBDIM,HTURBLEN,HTOM,HTURBLEN_CL,HCLOUD,PIMPL, & PTSTEP,TPFILE,PDXX,PDYY,PDZZ,PDZX,PDZY,PZZ, & PDIRCOSXW,PDIRCOSYW,PDIRCOSZW,PCOSSLOPE,PSINSLOPE, & PRHODJ,PTHVREF, & PSFTH,PSFRV,PSFSV,PSFU,PSFV, & PPABST,PUT,PVT,PWT,PTKET,PSVT,PSRCT, & - PBL_DEPTH,PSBL_DEPTH, & + PBL_DEPTH, PSBL_DEPTH, & PCEI,PCEI_MIN,PCEI_MAX,PCOEF_AMPL_SAT, & PTHLT,PRT, & PRUS,PRVS,PRWS,PRTHLS,PRRS,PRSVS,PRTKES,PRTKEMS,PSIGS,& @@ -142,7 +143,6 @@ SUBROUTINE TURB(KKA,KKU,KKL,KMI,KRR,KRRL,KRRI,HLBCX,HLBCY,KSPLIT,KMODEL_CL, & !! Module MODD_BUDGET: !! NBUMOD !! CBUTYPE -!! NBUPROCCTR !! LBU_RU !! LBU_RV !! LBU_RW @@ -230,24 +230,41 @@ SUBROUTINE TURB(KKA,KKU,KKL,KMI,KRR,KRRL,KRRI,HLBCX,HLBCY,KSPLIT,KMODEL_CL, & !! 10/2012 (J. Colin) Correct bug in DearDoff for dry simulations !! 10/2012 J.Escobar Bypass PGI bug , redefine some allocatable array inplace of automatic !! 04/2016 (C.Lac) correction of negativity for KHKO -!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O -!! 01/2018 (Q.Rodier) Introduction of RM17 +! P. Wautelet 05/2016-04/2018: new data structures and calls for I/O +! Q. Rodier 01/2018: introduction of RM17 ! P. Wautelet 20/05/2019: add name argument to ADDnFIELD_ll + new ADD4DFIELD_ll subroutine ! P. Wautelet 20/06/2019: take DELT and DEAR subroutines out of the TURB one (PGI compiler bug workaround) + transform into a mode_ module -!! -------------------------------------------------------------------------- -! +! 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 11/06/2020: bugfix: correct PRSVS array indices +! P. Wautelet + Benoit 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 +! R. Honnert/V. Masson 02/2021: new mixing length in the grey zone +! J.L. Redelsperger 03/2021: add Ocean LES case +! -------------------------------------------------------------------------- +! !* 0. DECLARATIONS ! ------------ ! -USE MODD_PARAMETERS, ONLY: JPVEXT_TURB +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_CST USE MODD_CTURB -USE MODD_CONF -USE MODD_BUDGET +USE MODD_DYN_n, ONLY : LOCEAN +use modd_field, only: tfielddata, TYPEREAL +USE MODD_IBM_PARAM_n, ONLY: LIBM, XIBM_LS, XIBM_XMUT USE MODD_IO, ONLY: TFILEDATA USE MODD_LES USE MODD_NSV +USE MODD_PARAMETERS, ONLY: JPVEXT_TURB +USE MODD_PARAM_LIMA +USE MODD_TURB_n, ONLY: XCADAP ! +USE MODI_IBM_MIXINGLENGTH USE MODI_GRADIENT_M USE MODI_GRADIENT_U USE MODI_GRADIENT_V @@ -264,7 +281,6 @@ USE MODI_SHUMAN USE MODI_SHUMAN_DEVICE #endif USE MODI_GRADIENT_M -USE MODI_BUDGET USE MODI_LES_MEAN_SUBGRID USE MODI_RMC01 USE MODI_GRADIENT_W @@ -272,18 +288,17 @@ USE MODI_TM06 USE MODI_UPDATE_LM USE MODI_GET_HALO ! -USE MODE_FIELD, ONLY: TFIELDDATA, TYPEREAL +use mode_budget, only: Budget_store_init, Budget_store_end USE MODE_IO_FIELD_WRITE, only: IO_Field_write +USE MODE_MPPDB USE MODE_SBL +use mode_sources_neg_correct, only: Sources_neg_correct ! USE MODI_EMOIST USE MODI_ETHETA ! USE MODI_SECOND_MNH ! -USE MODE_MPPDB -! -!! use, intrinsic :: ISO_C_BINDING ! IMPLICIT NONE ! @@ -302,8 +317,6 @@ INTEGER, INTENT(IN) :: KRRI ! number of ice water var. CHARACTER(LEN=*),DIMENSION(:),INTENT(IN):: HLBCX, HLBCY ! X- and Y-direc LBC INTEGER, INTENT(IN) :: KSPLIT ! number of time-splitting INTEGER, INTENT(IN) :: KMODEL_CL ! model number for cloud mixing length -LOGICAL, INTENT(IN) :: OCLOSE_OUT ! switch for syncronous - ! file opening LOGICAL, INTENT(IN) :: OTURB_FLX ! switch to write the ! turbulent fluxes in the syncronous FM-file LOGICAL, INTENT(IN) :: OTURB_DIAG ! switch to write some @@ -364,6 +377,7 @@ REAL, INTENT(IN) :: PCOEF_AMPL_SAT ! saturation of the amplification coeff REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PTHLT ! conservative pot. temp. REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRT ! water var. where ! PRT(:,:,:,1) is the conservative mixing ratio +! ! sources of momentum, conservative potential temperature, Turb. Kin. Energy, ! TKE dissipation REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRUS,PRVS,PRWS,PRTHLS,PRTKES @@ -397,6 +411,7 @@ REAL, POINTER , CONTIGUOUS, DIMENSION(:,:,:) ::& ZEXN, & ! EXN at t-1 ZT, & ! T at t-1 ZLOCPEXNM, & ! Lv/Cp/EXNREF at t-1 + ZLMW, & ! Turbulent mixing length (work array) ZLEPS, & ! Dissipative length ZTRH, & ! Dynamic and Thermal Production of TKE ZATHETA,ZAMOIST, & ! coefficients for s = f (Thetal,Rnp) @@ -405,9 +420,11 @@ REAL, POINTER , CONTIGUOUS, DIMENSION(:,:,:) ::& ZMWTH,ZMWR,ZMTH2,ZMR2,ZMTHR,& ! 3rd order moments ZFWTH,ZFWR,ZFTH2,ZFR2,ZFTHR,& ! opposite of verticale derivate of 3rd order moments ZTHLM ! initial potential temp. +#ifdef MNH_OPENACC INTEGER :: IZCP,IZEXN,IZT,IZLOCPEXNM,IZLEPS,IZTRH,IZATHETA,IZAMOIST & ,IZCOEF_DISS,IZFRAC_ICE,IZMWTH,IZMWR,IZMTH2,IZMR2,IZMTHR & ,IZFWTH,IZFWR,IZFTH2,IZFR2,IZFTHR,IZTHLM +#endif ! REAL, POINTER , CONTIGUOUS, DIMENSION(:,:,:,:) :: & ZRM ! initial mixing ratio @@ -426,12 +443,16 @@ REAL, POINTER , CONTIGUOUS, DIMENSION(:,:) :: ZTAU11M,ZTAU12M, & ! ! Virtual Potential Temp. used ! in the Deardorff mixing length computation +#ifdef MNH_OPENACC INTEGER :: IZRM,IZTAU11M,IZTAU12M,IZTAU22M,IZTAU33M,IZUSLOPE,IZVSLOPE & ,IZCDUEFF,IZUSTAR,IZLMO,IZRVM,IZSFRV +#endif REAL, DIMENSION(:,:,:), POINTER , CONTIGUOUS :: & ZLVOCPEXNM,ZLSOCPEXNM, & ! Lv/Cp/EXNREF and Ls/Cp/EXNREF at t-1 ZATHETA_ICE,ZAMOIST_ICE ! coefficients for s = f (Thetal,Rnp) +#ifdef MNH_OPENACC INTEGER :: IZLVOCPEXNM,IZLSOCPEXNM,IZATHETA_ICE,IZAMOIST_ICE +#endif ! REAL :: ZEXPL ! 1-PIMPL deg of expl. REAL :: ZRVORD ! RV/RD @@ -443,13 +464,16 @@ INTEGER :: IKTB,IKTE ! start, end of k loops in physical domain INTEGER :: JRR,JK,JSV ! loop counters INTEGER :: JI,JJ ! loop counters REAL :: ZL0 ! Max. Mixing Length in Blakadar formula -REAL :: ZALPHA ! proportionnality constant between Dz/2 and -! ! BL89 mixing length near the surface +REAL :: ZALPHA ! work coefficient : + ! - proportionnality constant between Dz/2 and +! ! BL89 mixing length near the surface ! REAL :: ZTIME1, ZTIME2 REAL, DIMENSION(:,:,:), pointer , contiguous :: ZTT,ZEXNE,ZLV,ZCPH REAL, DIMENSION(:,:,:), pointer , contiguous :: ZSHEAR, ZDUDZ, ZDVDZ +#ifdef MNH_OPENACC INTEGER :: IZTT,IZEXNE,IZLV,IZCPH,IZSHEAR, IZDUDZ, IZDVDZ +#endif TYPE(TFIELDDATA) :: TZFIELD ! #ifdef MNH_OPENACC @@ -695,7 +719,7 @@ ZRVORD= XRV / XRD !$acc update device(PTHLT,PRT) !$acc kernels !Copy data into ZTHLM and ZRM only if needed -IF (HTURBLEN=='BL89' .OR. HTURBLEN=='RM17' .OR. ORMC01) THEN +IF (HTURBLEN=='BL89' .OR. HTURBLEN=='RM17' .OR. HTURBLEN == 'ADAP' .OR. ORMC01) THEN ZTHLM(:,:,:) = PTHLT(:,:,:) ZRM(:,:,:,:) = PRT(:,:,:,:) END IF @@ -725,14 +749,18 @@ END DO ! !* 2.2 Exner function at t ! +IF (LOCEAN) THEN + ZEXN(:,:,:) = 1. +ELSE !PW: "BUG" PGI : results different on CPU and GPU due to the power function !See http://www.pgroup.com/userforum/viewtopic.php?t=5364&sid=ec7b762b17fb9bb3332a47f0db57af55 !Use of own functions allows bit-reproducible results #ifndef MNH_BITREP -ZEXN(:,:,:) = (PPABST(:,:,:)/XP00) ** (XRD/XCPD) + ZEXN(:,:,:) = (PPABST(:,:,:)/XP00) ** (XRD/XCPD) #else -ZEXN(:,:,:) = BR_POW(PPABST(:,:,:)/XP00,XRD/XCPD) + ZEXN(:,:,:) = BR_POW(PPABST(:,:,:)/XP00,XRD/XCPD) #endif +END IF ! !* 2.3 dissipative heating coeff a t ! @@ -790,7 +818,7 @@ IF (KRRL >=1) THEN END IF ! ! - IF (OCLOSE_OUT .AND. OTURB_DIAG) THEN + IF ( tpfile%lopened .AND. OTURB_DIAG ) THEN !$acc update self(ZAMOIST,ZATHETA) TZFIELD%CMNHNAME = 'ATHETA' TZFIELD%CSTDNAME = '' @@ -877,20 +905,41 @@ SELECT CASE (HTURBLEN) ZSHEAR = SQRT(ZDUDZ*ZDUDZ + ZDVDZ*ZDVDZ) CALL BL89(KKA,KKU,KKL,PZZ,PDZZ,PTHVREF,ZTHLM,KRR,ZRM,PTKET,ZSHEAR,PLEM) ! -!* 3.3 Delta mixing length +!* 3.3 Grey-zone combined RM17 & Deardorff mixing lengths +! -------------------------------------------------- + + CASE ('ADAP') +#ifdef MNH_OPENACC + call Print_msg( NVERB_FATAL, 'GEN', 'TURB', 'OpenACC: HTURBLEN=ADAP not yet implemented' ) +#endif + ZDUDZ = MXF(MZF(GZ_U_UW(PUT,PDZZ))) + ZDVDZ = MYF(MZF(GZ_V_VW(PVT,PDZZ))) + ZSHEAR = SQRT(ZDUDZ*ZDUDZ + ZDVDZ*ZDVDZ) + CALL BL89(KKA,KKU,KKL,PZZ,PDZZ,PTHVREF,ZTHLM,KRR,ZRM,PTKET,ZSHEAR,PLEM) + + CALL DELT(KKA,KKU,KKL,IKB, IKE,IKTB, IKTE,ORMC01,HTURBDIM,PDXX, PDYY,PZZ,PDIRCOSZW,ZLMW,ODZ=.FALSE.) + ! The minimum mixing length is chosen between Horizontal grid mesh (not taking into account the vertical grid mesh) and RM17. + ! For large horizontal grid meshes, this is equal to RM17 + ! For LES grid meshes, this is equivalent to Deardorff : the base mixing lentgh is the horizontal grid mesh, + ! and it is limited by a stability-based length (RM17), as was done in Deardorff length (but taking into account shear as well) + ! For grid meshes in the grey zone, then this is the smaller of the two. + PLEM = MIN(PLEM,XCADAP*ZLMW) +! +!* 3.4 Delta mixing length ! ------------------- ! CASE ('DELT') - CALL DELT(KKA,KKU,KKL,IKB, IKE,IKTB, IKTE,ORMC01,HTURBDIM,PDXX, PDYY,PZZ,PDIRCOSZW,PLEM) + CALL DELT(KKA,KKU,KKL,IKB, IKE,IKTB, IKTE,ORMC01,HTURBDIM,PDXX, PDYY,PZZ,PDIRCOSZW,PLEM,ODZ=.TRUE.) ! -!* 3.4 Deardorff mixing length +!* 3.5 Deardorff mixing length ! ----------------------- ! CASE ('DEAR') CALL DEAR(KKA,KKU,KKL,KRR, KRRI, IKB, IKE,IKTB, IKTE, & ORMC01,HTURBDIM,PDXX, PDYY, PDZZ,PZZ,PDIRCOSZW,PTHLT,PTHVREF,PTKET,PSRCT,PRT,& ZLOCPEXNM,ZATHETA, ZAMOIST, PLEM) -!* 3.5 Blackadar mixing length +! +!* 3.6 Blackadar mixing length ! ----------------------- ! CASE ('BLKR') @@ -944,27 +993,37 @@ IF (ORMC01) THEN call Print_msg( NVERB_FATAL, 'GEN', 'TURB', 'OpenACC: ORMC01 not yet implemented' ) #endif #ifndef MNH_BITREP - ZUSTAR=(PSFU**2+PSFV**2)**(0.25) + ZUSTAR(:,:) = (PSFU(:,:)**2+PSFV(:,:)**2)**(0.25) #else - ZUSTAR=BR_POW(BR_P2(PSFU)+BR_P2(PSFV),0.25) + ZUSTAR(:,:) = BR_POW( BR_P2( PSFU(:,:) ) + BR_P2( PSFV(:,:) ), 0.25 ) #endif IF (KRR>0) THEN - ZLMO=LMO(ZUSTAR,ZTHLM(:,:,IKB),ZRM(:,:,IKB,1),PSFTH,PSFRV) + ZLMO(:,:) = LMO( ZUSTAR(:,:), ZTHLM(:,:,IKB), ZRM(:,:,IKB,1), PSFTH(:,:), PSFRV(:,:) ) ELSE - ZRVM=0. - ZSFRV=0. - ZLMO=LMO(ZUSTAR,ZTHLM(:,:,IKB),ZRVM,PSFTH,ZSFRV) + ZRVM (:,:) = 0. + ZSFRV(:,:) = 0. + ZLMO(:,:) = LMO( ZUSTAR(:,:), ZTHLM(:,:,IKB), ZRVM(:,:), PSFTH(:,:), ZSFRV(:,:) ) END IF - CALL RMC01(HTURBLEN,KKA,KKU,KKL,PZZ,PDXX,PDYY,PDZZ,PDIRCOSZW,PSBL_DEPTH,ZLMO,PLEM,ZLEPS) + CALL RMC01( HTURBLEN, KKA, KKU, KKL, PZZ, PDXX, PDYY, PDZZ, PDIRCOSZW, PSBL_DEPTH, ZLMO, PLEM, ZLEPS ) !$acc update device(PLEM,ZLEPS) END IF ! +!RMC01 is only applied on RM17 in ADAP +IF (HTURBLEN=='ADAP') ZLEPS = MIN(ZLEPS,ZLMW*XCADAP) +! !* 3.8 Mixing length in external points (used if HTURBDIM="3DIM") ! ---------------------------------------------------------- ! IF (HTURBDIM=="3DIM") THEN CALL UPDATE_LM(HLBCX,HLBCY,PLEM,ZLEPS) END IF +! +!* 3.9 Mixing length correction if immersed walls +! ------------------------------------------ +! +IF (LIBM) THEN + CALL IBM_MIXINGLENGTH(PLEM,ZLEPS,XIBM_XMUT,XIBM_LS(:,:,:,1),PTKET) +ENDIF !---------------------------------------------------------------------------- ! !* 4. GO INTO THE AXES FOLLOWING THE SURFACE @@ -1084,10 +1143,44 @@ ENDIF !* 5. TURBULENT SOURCES ! ----------------- ! +if ( lbudget_u ) call Budget_store_init( tbudgets(NBUDGET_U ), 'VTURB', prus (:, :, :) ) +if ( lbudget_v ) call Budget_store_init( tbudgets(NBUDGET_V ), 'VTURB', prvs (:, :, :) ) +if ( lbudget_w ) call Budget_store_init( tbudgets(NBUDGET_W ), 'VTURB', prws (:, :, :) ) + +if ( lbudget_th ) then + if ( krri >= 1 .and. krrl >= 1 ) then + call Budget_store_init( tbudgets(NBUDGET_TH), 'VTURB', prthls(:, :, :) + zlvocpexnm(:, :, :) * prrs(:, :, :, 2) & + + zlsocpexnm(:, :, :) * prrs(:, :, :, 4) ) + else if ( krrl >= 1 ) then + call Budget_store_init( tbudgets(NBUDGET_TH), 'VTURB', prthls(:, :, :) + zlocpexnm(:, :, :) * prrs(:, :, :, 2) ) + else + call Budget_store_init( tbudgets(NBUDGET_TH), 'VTURB', prthls(:, :, :) ) + end if +end if + +if ( lbudget_rv ) then + if ( krri >= 1 .and. krrl >= 1 ) then + call Budget_store_init( tbudgets(NBUDGET_RV), 'VTURB', prrs(:, :, :, 1) - prrs(:, :, :, 2) - prrs(:, :, :, 4) ) + else if ( krrl >= 1 ) then + call Budget_store_init( tbudgets(NBUDGET_RV), 'VTURB', prrs(:, :, :, 1) - prrs(:, :, :, 2) ) + else + call Budget_store_init( tbudgets(NBUDGET_RV), 'VTURB', prrs(:, :, :, 1) ) + end if +end if + +if ( lbudget_rc ) call Budget_store_init( tbudgets(NBUDGET_RC), 'VTURB', prrs (:, :, :, 2) ) +if ( lbudget_ri ) call Budget_store_init( tbudgets(NBUDGET_RI), 'VTURB', prrs (:, :, :, 4) ) + +if ( lbudget_sv ) then + do jsv = 1, nsv + call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + jsv), 'VTURB', prsvs(:, :, :, jsv) ) + end do +end if + !$acc update device(PRHODJ) !$acc update device(PRUS,PRVS,PRWS,PRSVS) CALL TURB_VER(KKA,KKU,KKL,KRR, KRRL, KRRI, & - OCLOSE_OUT,OTURB_FLX, & + OTURB_FLX, & HTURBDIM,HTOM,PIMPL,ZEXPL, & PTSTEP,TPFILE, & PDXX,PDYY,PDZZ,PDZX,PDZY,PDIRCOSZW,PZZ, & @@ -1103,60 +1196,78 @@ CALL TURB_VER(KKA,KKU,KKL,KRR, KRRL, KRRI, & PRUS,PRVS,PRWS,PRTHLS,PRRS,PRSVS, & PDYP,PTHP,PSIGS,PWTH,PWRC,PWSV ) !$acc update self(PWTH,PWRC,PWSV) + +if ( lbudget_u ) call Budget_store_end( tbudgets(NBUDGET_U), 'VTURB', prus(:, :, :) ) +if ( lbudget_v ) call Budget_store_end( tbudgets(NBUDGET_V), 'VTURB', prvs(:, :, :) ) +if ( lbudget_w ) call Budget_store_end( tbudgets(NBUDGET_W), 'VTURB', prws(:, :, :) ) + +if ( lbudget_th ) then + if ( krri >= 1 .and. krrl >= 1 ) then + call Budget_store_end( tbudgets(NBUDGET_TH), 'VTURB', prthls(:, :, :) + zlvocpexnm(:, :, :) * prrs(:, :, :, 2) & + + zlsocpexnm(:, :, :) * prrs(:, :, :, 4) ) + else if ( krrl >= 1 ) then + call Budget_store_end( tbudgets(NBUDGET_TH), 'VTURB', prthls(:, :, :) + zlocpexnm(:, :, :) * prrs(:, :, :, 2) ) + else + call Budget_store_end( tbudgets(NBUDGET_TH), 'VTURB', prthls(:, :, :) ) + end if +end if + +if ( lbudget_rv ) then + if ( krri >= 1 .and. krrl >= 1 ) then + call Budget_store_end( tbudgets(NBUDGET_RV), 'VTURB', prrs(:, :, :, 1) - prrs(:, :, :, 2) - prrs(:, :, :, 4) ) + else if ( krrl >= 1 ) then + call Budget_store_end( tbudgets(NBUDGET_RV), 'VTURB', prrs(:, :, :, 1) - prrs(:, :, :, 2) ) + else + call Budget_store_end( tbudgets(NBUDGET_RV), 'VTURB', prrs(:, :, :, 1) ) + end if +end if + +if ( lbudget_rc ) call Budget_store_end( tbudgets(NBUDGET_RC), 'VTURB', prrs(:, :, :, 2) ) +if ( lbudget_ri ) call Budget_store_end( tbudgets(NBUDGET_RI), 'VTURB', prrs(:, :, :, 4) ) + +if ( lbudget_sv ) then + do jsv = 1, nsv + call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + jsv), 'VTURB', prsvs(:, :, :, jsv) ) + end do +end if ! -#ifdef MNH_OPENACC -IF ( ( LBUDGET_TH .AND. ( ( KRRI >= 1 .AND. KRRL >= 1 ) .OR. ( KRRL >= 1 ) ) ) .OR. & - LBUDGET_RV .OR. LBUDGET_RC .OR. LBUDGET_RI ) THEN -!$acc update self(PRRS) -ENDIF -#endif -! -IF (LBUDGET_U) THEN -!$acc update self(PRUS) - CALL BUDGET (PRUS,1,'VTURB_BU_RU') -END IF -IF (LBUDGET_V) THEN -!$acc update self(PRVS) - CALL BUDGET (PRVS,2,'VTURB_BU_RV') -END IF -IF (LBUDGET_W) THEN -!$acc update self(PRWS) - CALL BUDGET (PRWS,3,'VTURB_BU_RW') -END IF -IF (LBUDGET_TH) THEN -!$acc update self(PRTHLS) - IF ( KRRI >= 1 .AND. KRRL >= 1 ) THEN -!$acc update self(ZLVOCPEXNM,ZLSOCPEXNM) - CALL BUDGET (PRTHLS+ ZLVOCPEXNM * PRRS(:,:,:,2) + ZLSOCPEXNM * PRRS(:,:,:,4),4,'VTURB_BU_RTH') - ELSE IF ( KRRL >= 1 ) THEN -!$acc update self(ZLOCPEXNM) - CALL BUDGET (PRTHLS+ ZLOCPEXNM * PRRS(:,:,:,2),4,'VTURB_BU_RTH') - ELSE - CALL BUDGET (PRTHLS,4,'VTURB_BU_RTH') - END IF -END IF -IF (LBUDGET_SV) THEN -!$acc update self(PRSVS) - DO JSV = 1,NSV - CALL BUDGET (PRSVS(:,:,:,JSV),JSV+12,'VTURB_BU_RSV') - END DO -END IF -IF (LBUDGET_RV) THEN - IF ( KRRI >= 1 .AND. KRRL >= 1) THEN - CALL BUDGET (PRRS(:,:,:,1)-PRRS(:,:,:,2)-PRRS(:,:,:,4),6,'VTURB_BU_RRV') - ELSE IF ( KRRL >= 1 ) THEN - CALL BUDGET (PRRS(:,:,:,1)-PRRS(:,:,:,2),6,'VTURB_BU_RRV') - ELSE - CALL BUDGET (PRRS(:,:,:,1),6,'VTURB_BU_RRV') - END IF -END IF -IF (LBUDGET_RC) CALL BUDGET (PRRS(:,:,:,2),7,'VTURB_BU_RRC') -IF (LBUDGET_RI) CALL BUDGET (PRRS(:,:,:,4),9,'VTURB_BU_RRI') -! -! -IF (HTURBDIM=='3DIM') THEN +if ( hturbdim == '3DIM' ) then + if ( lbudget_u ) call Budget_store_init( tbudgets(NBUDGET_U ), 'HTURB', prus (:, :, :) ) + if ( lbudget_v ) call Budget_store_init( tbudgets(NBUDGET_V ), 'HTURB', prvs (:, :, :) ) + if ( lbudget_w ) call Budget_store_init( tbudgets(NBUDGET_W ), 'HTURB', prws (:, :, :) ) + + if (lbudget_th) then + if ( krri >= 1 .and. krrl >= 1 ) then + call Budget_store_init( tbudgets(NBUDGET_TH), 'HTURB', prthls(:, :, :) + zlvocpexnm(:, :, :) * prrs(:, :, :, 2) & + + zlsocpexnm(:, :, :) * prrs(:, :, :, 4) ) + else if ( krrl >= 1 ) then + call Budget_store_init( tbudgets(NBUDGET_TH), 'HTURB', prthls(:, :, :) + zlocpexnm(:, :, :) * prrs(:, :, :, 2) ) + else + call Budget_store_init( tbudgets(NBUDGET_TH), 'HTURB', prthls(:, :, :) ) + end if + end if + + if ( lbudget_rv ) then + if ( krri >= 1 .and. krrl >= 1 ) then + call Budget_store_init( tbudgets(NBUDGET_RV), 'HTURB', prrs(:, :, :, 1) - prrs(:, :, :, 2) - prrs(:, :, :, 4) ) + else if ( krrl >= 1 ) then + call Budget_store_init( tbudgets(NBUDGET_RV), 'HTURB', prrs(:, :, :, 1) - prrs(:, :, :, 2) ) + else + call Budget_store_init( tbudgets(NBUDGET_RV), 'HTURB', prrs(:, :, :, 1) ) + end if + end if + + if ( lbudget_rc ) call Budget_store_init( tbudgets(NBUDGET_RC), 'HTURB', prrs(:, :, :, 2) ) + if ( lbudget_ri ) call Budget_store_init( tbudgets(NBUDGET_RI), 'HTURB', prrs(:, :, :, 4) ) + + if ( lbudget_sv ) then + do jsv = 1, nsv + call Budget_store_init( tbudgets(NBUDGET_SV1 - 1 + jsv), 'HTURB', prsvs(:, :, :, jsv) ) + end do + end if + CALL TURB_HOR_SPLT(KSPLIT, KRR, KRRL, KRRI, PTSTEP, & - HLBCX,HLBCY,OCLOSE_OUT,OTURB_FLX,OSUBG_COND, & + HLBCX,HLBCY,OTURB_FLX,OSUBG_COND, & TPFILE, & PDXX,PDYY,PDZZ,PDZX,PDZY,PZZ, & PDIRCOSXW,PDIRCOSYW,PDIRCOSZW, & @@ -1170,48 +1281,44 @@ IF (HTURBDIM=='3DIM') THEN PDYP,PTHP,PSIGS, & ZTRH, & PRUS,PRVS,PRWS,PRTHLS,PRRS,PRSVS ) -END IF + + if ( lbudget_u ) call Budget_store_end( tbudgets(NBUDGET_U), 'HTURB', prus(:, :, :) ) + if ( lbudget_v ) call Budget_store_end( tbudgets(NBUDGET_V), 'HTURB', prvs(:, :, :) ) + if ( lbudget_w ) call Budget_store_end( tbudgets(NBUDGET_W), 'HTURB', prws(:, :, :) ) + + if ( lbudget_th ) then + if ( krri >= 1 .and. krrl >= 1 ) then + call Budget_store_end( tbudgets(NBUDGET_TH), 'HTURB', prthls(:, :, :) + zlvocpexnm(:, :, :) * prrs(:, :, :, 2) & + + zlsocpexnm(:, :, :) * prrs(:, :, :, 4) ) + else if ( krrl >= 1 ) then + call Budget_store_end( tbudgets(NBUDGET_TH), 'HTURB', prthls(:, :, :) + zlocpexnm(:, :, :) * prrs(:, :, :, 2) ) + else + call Budget_store_end( tbudgets(NBUDGET_TH), 'HTURB', prthls(:, :, :) ) + end if + end if + + if ( lbudget_rv ) then + if ( krri >= 1 .and. krrl >= 1 ) then + call Budget_store_end( tbudgets(NBUDGET_RV), 'HTURB', prrs(:, :, :, 1) - prrs(:, :, :, 2) - prrs(:, :, :, 4) ) + else if ( krrl >= 1 ) then + call Budget_store_end( tbudgets(NBUDGET_RV), 'HTURB', prrs(:, :, :, 1) - prrs(:, :, :, 2) ) + else + call Budget_store_end( tbudgets(NBUDGET_RV), 'HTURB', prrs(:, :, :, 1) ) + end if + end if + + if ( lbudget_rc ) call Budget_store_end( tbudgets(NBUDGET_RC), 'HTURB', prrs(:, :, :, 2) ) + if ( lbudget_ri ) call Budget_store_end( tbudgets(NBUDGET_RI), 'HTURB', prrs(:, :, :, 4) ) + + if ( lbudget_sv ) then + do jsv = 1, nsv + call Budget_store_end( tbudgets(NBUDGET_SV1 - 1 + jsv), 'HTURB', prsvs(:, :, :, jsv) ) + end do + end if +end if + !$acc update self(PSIGS,PRUS,PRVS,PRWS,PRSVS) -! -! -#ifdef MNH_OPENACC -IF ( ( LBUDGET_TH .AND. ( ( KRRI >= 1 .AND. KRRL >= 1 ) .OR. ( KRRL >= 1 ) ) ) .OR. & - LBUDGET_RV .OR. LBUDGET_RC .OR. LBUDGET_RI ) THEN -!$acc update self(PRRS) -ENDIF -#endif -! -IF (LBUDGET_U) CALL BUDGET (PRUS,1,'HTURB_BU_RU') -IF (LBUDGET_V) CALL BUDGET (PRVS,2,'HTURB_BU_RV') -IF (LBUDGET_W) CALL BUDGET (PRWS,3,'HTURB_BU_RW') -IF (LBUDGET_TH) THEN -!$acc update self(PRTHLS) - IF ( KRRI >= 1 .AND. KRRL >= 1 ) THEN - CALL BUDGET (PRTHLS+ZLVOCPEXNM*PRRS(:,:,:,2)+ZLSOCPEXNM*PRRS(:,:,:,4) & - ,4,'HTURB_BU_RTH') - ELSE IF ( KRRL >= 1 ) THEN - CALL BUDGET (PRTHLS+ ZLOCPEXNM * PRRS(:,:,:,2),4,'HTURB_BU_RTH') - ELSE - CALL BUDGET (PRTHLS,4,'HTURB_BU_RTH') - END IF -END IF -IF (LBUDGET_SV) THEN - DO JSV = 1,NSV - CALL BUDGET (PRSVS(:,:,:,JSV),JSV+12,'HTURB_BU_RSV') - END DO -END IF -IF (LBUDGET_RV) THEN - IF ( KRRI >= 1 .AND. KRRL >= 1) THEN - CALL BUDGET (PRRS(:,:,:,1)-PRRS(:,:,:,2)-PRRS(:,:,:,4),6,'HTURB_BU_RRV') - ELSE IF ( KRRL >= 1 ) THEN - CALL BUDGET (PRRS(:,:,:,1)-PRRS(:,:,:,2),6,'HTURB_BU_RRV') - ELSE - CALL BUDGET (PRRS(:,:,:,1),6,'HTURB_BU_RRV') - END IF -END IF -IF (LBUDGET_RC) CALL BUDGET (PRRS(:,:,:,2),7,'HTURB_BU_RRC') -IF (LBUDGET_RI) CALL BUDGET (PRRS(:,:,:,4),9,'HTURB_BU_RRI') -! + !---------------------------------------------------------------------------- ! !* 6. EVOLUTION OF THE TKE AND ITS DISSIPATION @@ -1234,33 +1341,21 @@ CALL TKE_EPS_SOURCES(KKA,KKU,KKL,KMI,PTKET,PLEM,ZLEPS,PDYP,ZTRH, & PRHODJ,PDZZ,PDXX,PDYY,PDZX,PDZY,PZZ, & PTSTEP,PIMPL,ZEXPL, & HTURBLEN,HTURBDIM, & - TPFILE,OCLOSE_OUT,OTURB_DIAG, & + TPFILE,OTURB_DIAG, & PTHP,PRTKES,PRTKEMS,PRTHLS,ZCOEF_DISS,PTR,PDISS ) ! !$acc update self(PTR,PDISS) !$acc update self(PDYP,PTHP) ! !$acc update self(PRTKES) -IF (LBUDGET_TH) THEN -!$acc update self(PRTHLS) - IF ( KRRI >= 1 .AND. KRRL >= 1 ) THEN - CALL BUDGET (PRTHLS+ZLVOCPEXNM*PRRS(:,:,:,2)+ZLSOCPEXNM*PRRS(:,:,:,4) & - ,4,'DISSH_BU_RTH') - ELSE IF ( KRRL >= 1 ) THEN - CALL BUDGET (PRTHLS+ZLOCPEXNM* PRRS(:,:,:,2),4,'DISSH_BU_RTH') - ELSE - CALL BUDGET (PRTHLS,4,'DISSH_BU_RTH') - END IF -END IF -! + !---------------------------------------------------------------------------- ! !* 7. STORES SOME INFORMATIONS RELATED TO THE TURBULENCE SCHEME ! --------------------------------------------------------- ! !$acc update self(PLEM) -! -IF ( OTURB_DIAG .AND. OCLOSE_OUT ) THEN +IF ( OTURB_DIAG .AND. tpfile%lopened ) THEN ! ! stores the mixing length ! @@ -1340,37 +1435,12 @@ IF ( KRRL >= 1 ) THEN !$acc update self(PRT(:,:,:,1)) END IF END IF + !$acc update self(PRRS,PTHLT,PRTHLS) -! -IF ((HCLOUD == 'KHKO') .OR. (HCLOUD == 'C2R2')) THEN -#ifdef MNH_OPENACC - call Print_msg( NVERB_FATAL, 'GEN', 'TURB', 'OpenACC: HCLOUD=KHKO or C2R2 not yet implemented' ) -#endif - ZEXNE(:,:,:)= (PPABST(:,:,:)/XP00)**(XRD/XCPD) - ZTT(:,:,:)= PTHLT(:,:,:)*ZEXNE(:,:,:) - ZLV(:,:,:)=XLVTT +(XCPV-XCL) *(ZTT(:,:,:)-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) - PRTHLS(:,:,:) = PRTHLS(:,:,:) - PRRS(:,:,:,JSV) * ZLV(:,:,:) / & - ZCPH(:,:,:) / ZEXNE(:,:,:) - PRRS(:,:,:,JSV) = 0.0 - PRSVS(:,:,:,JSV) = 0.0 - END WHERE - END DO -! - IF (LBUDGET_TH) CALL BUDGET (PRTHLS(:,:,:), 4,'NETUR_BU_RTH') - IF (LBUDGET_RV) CALL BUDGET (PRRS(:,:,:,1), 6,'NETUR_BU_RRV') - IF (LBUDGET_RC) CALL BUDGET (PRRS(:,:,:,2), 7,'NETUR_BU_RRC') -END IF -! + +! Remove non-physical negative values (unnecessary in a perfect world) + corresponding budgets +call Sources_neg_correct( hcloud, 'NETUR', krr, ptstep, ppabst, pthlt, prt, prthls, prrs, prsvs ) + !---------------------------------------------------------------------------- ! !* 9. LES averaged surface fluxes @@ -1613,19 +1683,22 @@ CALL GET_INDICE_ll (IIB,IJB,IIE,IJE) ! ! 2 Update halo if necessary ! +#ifndef MNH_OPENACC !!$IF (NHALO == 1) THEN !!$!$acc update self(PUSLOPE,PVSLOPE) -!!$ CALL ADD2DFIELD_ll( TZFIELDS_ll, PUSLOPE, 'UPDATE_ROTATE_WIND::PUSLOPE' ) -!!$ CALL ADD2DFIELD_ll( TZFIELDS_ll, PVSLOPE, 'UPDATE_ROTATE_WIND::PVSLOPE' ) -!!$ CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll) -!!$ CALL CLEANLIST_ll(TZFIELDS_ll) + CALL ADD2DFIELD_ll( TZFIELDS_ll, PUSLOPE, 'UPDATE_ROTATE_WIND::PUSLOPE' ) + CALL ADD2DFIELD_ll( TZFIELDS_ll, PVSLOPE, 'UPDATE_ROTATE_WIND::PVSLOPE' ) + CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll) + CALL CLEANLIST_ll(TZFIELDS_ll) !!$!$acc update device(PUSLOPE,PVSLOPE) +!!$ENDIF +#else ! ! /!\ warning conner needed -> GET_HALO...C ! CALL GET_2D_HALO_DDC( PUSLOPE, HNAME='UPDATE_ROTATE_WIND::PUSLOPE' ) CALL GET_2D_HALO_DDC( PVSLOPE, HNAME='UPDATE_ROTATE_WIND::PVSLOPE' ) -!!$ENDIF +#endif ! ! 3 Boundary conditions for non cyclic case ! @@ -1700,7 +1773,9 @@ REAL, DIMENSION(:,:,:), INTENT(OUT) :: PAMOIST,PATHETA REAL :: ZEPS ! XMV / XMD real, dimension(:,:,:), pointer , contiguous :: zrvsat real, dimension(:,:,:), pointer , contiguous :: zdrvsatdt +#ifdef MNH_OPENACC INTEGER :: izrvsat, izdrvsatdt +#endif ! !------------------------------------------------------------------------------- @@ -1911,7 +1986,7 @@ ELSE !* 3.2 Delta mixing length ! ------------------- CASE ('DELT') - CALL DELT(KKA,KKU,KKL,IKB, IKE,IKTB, IKTE,ORMC01,HTURBDIM,PDXX, PDYY,PZZ,PDIRCOSZW,ZLM_CLOUD) + CALL DELT(KKA,KKU,KKL,IKB, IKE,IKTB, IKTE,ORMC01,HTURBDIM,PDXX, PDYY,PZZ,PDIRCOSZW,ZLM_CLOUD,ODZ=.TRUE.) ! !* 3.3 Deardorff mixing length ! ----------------------- @@ -1927,7 +2002,7 @@ ENDIF ! ----------------------------------------------- ! ! Impression before modification of the mixing length -IF ( OTURB_DIAG .AND. OCLOSE_OUT ) THEN +IF ( OTURB_DIAG .AND. tpfile%lopened ) THEN TZFIELD%CMNHNAME = 'LM_CLEAR_SKY' TZFIELD%CSTDNAME = '' TZFIELD%CLONGNAME = 'LM_CLEAR_SKY' @@ -1953,7 +2028,7 @@ WHERE (PCEI(:,:,:) == -1.) PLEM(:,:,:) = ZLM_CLOUD(:,:,:) !* 5. IMPRESSION ! ---------- ! -IF ( OTURB_DIAG .AND. OCLOSE_OUT ) THEN +IF ( OTURB_DIAG .AND. tpfile%lopened ) THEN TZFIELD%CMNHNAME = 'COEF_AMPL' TZFIELD%CSTDNAME = '' TZFIELD%CLONGNAME = 'COEF_AMPL' @@ -1988,7 +2063,7 @@ END SUBROUTINE TURB !################### -SUBROUTINE DELT(KKA,KKU,KKL,KKB, KKE,KKTB, KKTE,ORMC01,HTURBDIM,PDXX, PDYY,PZZ,PDIRCOSZW,PLM) +SUBROUTINE DELT(KKA,KKU,KKL,KKB, KKE,KKTB, KKTE,ORMC01,HTURBDIM,PDXX, PDYY,PZZ,PDIRCOSZW,PLM, ODZ) !################### !! !!**** *DELT* routine to compute mixing length for DELT case @@ -2007,7 +2082,8 @@ SUBROUTINE DELT(KKA,KKU,KKL,KKB, KKE,KKTB, KKTE,ORMC01,HTURBDIM,PDXX, PDYY,PZZ,P ! !* 0. DECLARATIONS ! ------------ -use modd_conf, only: l2d +use modd_conf, only: l2d +USE MODD_DYN_n, ONLY: LOCEAN use mode_mppdb @@ -2037,6 +2113,7 @@ REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! physical distance REAL, DIMENSION(:,:), INTENT(IN) :: PDIRCOSZW ! Director Cosinus along x, y and z directions at surface w-point REAL, DIMENSION(:,:,:), INTENT(OUT) :: PLM +LOGICAL, INTENT(IN) :: ODZ ! ! !* 0.2 Declarations of local variables @@ -2069,43 +2146,79 @@ allocate( ztmp2_device( size( pdxx, 1 ), size( pdxx, 2 ), size( pdxx, 3 ) ) ) !$acc data create( ztmp1_device, ztmp2_device ) +IF (ODZ) THEN !$acc kernels -DO JK = KKTB,KKTE ! 1D turbulence scheme - PLM(:,:,JK) = PZZ(:,:,JK+KKL) - PZZ(:,:,JK) -END DO -PLM(:,:,KKU) = PLM(:,:,KKE) -PLM(:,:,KKA) = PZZ(:,:,KKB) - PZZ(:,:,KKA) + ! Dz is take into account in the computation + DO JK = KKTB,KKTE ! 1D turbulence scheme + PLM(:,:,JK) = PZZ(:,:,JK+KKL) - PZZ(:,:,JK) + END DO + PLM(:,:,KKU) = PLM(:,:,KKE) + PLM(:,:,KKA) = PZZ(:,:,KKB) - PZZ(:,:,KKA) !$acc end kernels -IF ( HTURBDIM /= '1DIM' ) THEN ! 3D turbulence scheme - IF ( L2D ) THEN + IF ( HTURBDIM /= '1DIM' ) THEN ! 3D turbulence scheme + IF ( L2D ) THEN #ifndef MNH_OPENACC - PLM(:,:,:) = SQRT( PLM(:,:,:)*MXF(PDXX(:,:,:)) ) + PLM(:,:,:) = SQRT( PLM(:,:,:)*MXF(PDXX(:,:,:)) ) #else - CALL MXF_DEVICE( PDXX, ZTMP1_DEVICE ) + CALL MXF_DEVICE( PDXX, ZTMP1_DEVICE ) !$acc kernels - PLM(:,:,:) = SQRT( PLM(:,:,:) * ZTMP1_DEVICE ) + PLM(:,:,:) = SQRT( PLM(:,:,:) * ZTMP1_DEVICE ) !$acc end kernels #endif - ELSE + ELSE #ifndef MNH_OPENACC #ifndef MNH_BITREP - PLM(:,:,:) = (PLM(:,:,:)*MXF(PDXX(:,:,:))*MYF(PDYY(:,:,:)) ) ** (1./3.) + PLM(:,:,:) = (PLM(:,:,:)*MXF(PDXX(:,:,:))*MYF(PDYY(:,:,:)) ) ** (1./3.) +#else + PLM(:,:,:) = BR_POW( PLM(:, :, : ) * MXF( PDXX(:, :, : ) ) * MYF( PDYY(:, :, : ) ), 1. / 3. ) +#endif +#else + CALL MXF_DEVICE( PDXX, ZTMP1_DEVICE ) + CALL MYF_DEVICE( PDYY, ZTMP2_DEVICE ) +!$acc kernels +#ifndef MNH_BITREP + PLM(:,:,:) = ( PLM(:,:,:) * ZTMP1_DEVICE * ZTMP2_DEVICE ) ** (1./3.) +#else + PLM(:,:,:) = BR_POW( PLM(:,:,:) * ZTMP1_DEVICE * ZTMP2_DEVICE, 1./3. ) +#endif +!$acc end kernels +#endif + END IF + END IF +ELSE + ! Dz not taken into account in computation to assure invariability with vertical grid mesh +!$acc kernels + PLM=1.E10 +!$acc end kernels + IF ( HTURBDIM /= '1DIM' ) THEN ! 3D turbulence scheme + IF ( L2D) THEN +#ifndef MNH_OPENACC + PLM(:,:,:) = MXF(PDXX(:,:,:)) +#else + CALL MXF_DEVICE( PDXX, PLM ) +#endif + ELSE +#ifndef MNH_OPENACC +#ifndef MNH_BITREP + PLM(:,:,:) = (MXF(PDXX(:,:,:))*MYF(PDYY(:,:,:)) ) ** (1./2.) #else - PLM(:,:,:) = BR_POW( PLM(:, :, : ) * MXF( PDXX(:, :, : ) ) * MYF( PDYY(:, :, : ) ), 1. / 3. ) + PLM(:,:,:) = BR_POW(MXF(PDXX(:,:,:))*MYF(PDYY(:,:,:)), 1. / 2. ) #endif #else - CALL MXF_DEVICE( PDXX, ZTMP1_DEVICE ) - CALL MYF_DEVICE( PDYY, ZTMP2_DEVICE ) + CALL MXF_DEVICE( PDXX, ZTMP1_DEVICE ) + CALL MYF_DEVICE( PDYY, ZTMP2_DEVICE ) !$acc kernels #ifndef MNH_BITREP - PLM(:,:,:) = ( PLM(:,:,:) * ZTMP1_DEVICE * ZTMP2_DEVICE ) ** (1./3.) + PLM(:,:,:) = ( ZTMP1_DEVICE * ZTMP2_DEVICE ) ** (1./2.) #else - PLM(:,:,:) = BR_POW( PLM(:,:,:) * ZTMP1_DEVICE * ZTMP2_DEVICE, 1./3. ) + PLM(:,:,:) = BR_POW( ZTMP1_DEVICE * ZTMP2_DEVICE, 1. / 2. ) #endif !$acc end kernels #endif + END IF END IF END IF + ! ! mixing length limited by the distance normal to the surface ! (with the same factor as for BL89) @@ -2120,14 +2233,25 @@ IF (.NOT. ORMC01) THEN ! DO JJ=1,SIZE(PLM,2) DO JI=1,SIZE(PLM,1) - DO JK=KKTB,KKTE - ZD = ZALPHA * ( 0.5 * ( PZZ(JI, JJ, JK ) + PZZ(JI, JJ, JK+KKL ) ) - PZZ(JI, JJ, KKB ) ) * PDIRCOSZW(JI, JJ ) - IF ( PLM(JI,JJ,JK) > ZD ) THEN - PLM(JI,JJ,JK) = ZD - ELSE - EXIT - ENDIF - END DO + IF (LOCEAN) THEN + DO JK=KKTE,KKTB,-1 + ZD=ZALPHA*(PZZ(JI,JJ,KKTE+1)-PZZ(JI,JJ,JK)) + IF ( PLM(JI,JJ,JK)>ZD) THEN + PLM(JI,JJ,JK)=ZD + ELSE + EXIT + ENDIF + END DO + ELSE + DO JK=KKTB,KKTE + ZD = ZALPHA * ( 0.5 * ( PZZ(JI, JJ, JK) + PZZ(JI, JJ, JK+KKL) ) - PZZ(JI, JJ, KKB) ) * PDIRCOSZW(JI, JJ) + IF ( PLM(JI,JJ,JK) > ZD ) THEN + PLM(JI,JJ,JK) = ZD + ELSE + EXIT + ENDIF + END DO + ENDIF END DO END DO END IF @@ -2172,8 +2296,9 @@ SUBROUTINE DEAR(KKA,KKU,KKL,KRR, KRRI, KKB, KKE,KKTB, KKTE, & ! !* 0. DECLARATIONS ! ------------ -use modd_conf, only: l2d -use modd_cst, only: XG, XMNH_EPSILON +use modd_conf, only: l2d +use modd_cst, only: XALPHAOC, XBETAOC, XG, XMNH_EPSILON +USE MODD_DYN_n, ONLY: LOCEAN use mode_mppdb @@ -2231,7 +2356,9 @@ REAL, DIMENSION(:,:,:), POINTER , CONTIGUOUS :: & ZDTHLDZ,ZDRTDZ, &!dtheta_l/dz, drt_dz used for computing the stablity ! ! criterion ZETHETA,ZEMOIST !coef ETHETA and EMOIST +#ifdef MNH_OPENACC INTEGER :: IZWORK2D,IZDTHLDZ,IZDRTDZ,IZETHETA,IZEMOIST +#endif ! #ifdef MNH_OPENACC REAL, DIMENSION(:,:,:), POINTER , CONTIGUOUS :: ZTMP1_DEVICE,ZTMP2_DEVICE @@ -2362,7 +2489,6 @@ CALL EMOIST(KRR,KRRI,PTHLT,PRT,PLOCPEXNM,PAMOIST,PSRCT,ZEMOIST) #endif ! !$acc kernels -! For dry simulations IF (KRR>0) THEN !$acc loop independent collapse(3) private(ZVAR) DO JK = KKTB+1,KKTE-1 @@ -2372,8 +2498,12 @@ IF (KRR>0) THEN (PTHLT(JI,JJ,JK )-PTHLT(JI,JJ,JK-KKL))/PDZZ(JI,JJ,JK )) ZDRTDZ(JI,JJ,JK) = 0.5*((PRT(JI,JJ,JK+KKL,1)-PRT(JI,JJ,JK ,1))/PDZZ(JI,JJ,JK+KKL)+ & (PRT(JI,JJ,JK ,1)-PRT(JI,JJ,JK-KKL,1))/PDZZ(JI,JJ,JK )) - ZVAR=XG/PTHVREF(JI,JJ,JK)* & + IF (LOCEAN) THEN + ZVAR=XG*(XALPHAOC*ZDTHLDZ(JI,JJ,JK)-XBETAOC*ZDRTDZ(JI,JJ,JK)) + ELSE + ZVAR=XG/PTHVREF(JI,JJ,JK)* & (ZETHETA(JI,JJ,JK)*ZDTHLDZ(JI,JJ,JK)+ZEMOIST(JI,JJ,JK)*ZDRTDZ(JI,JJ,JK)) + END IF ! IF (ZVAR>0.) THEN PLM(JI,JJ,JK)=MAX(XMNH_EPSILON,MIN(PLM(JI,JJ,JK), & @@ -2382,15 +2512,19 @@ IF (KRR>0) THEN END DO END DO END DO -ELSE +ELSE! For dry atmos or unsalted ocean runs !$acc loop independent collapse(3) private(ZVAR) DO JK = KKTB+1,KKTE-1 DO JJ=1,SIZE(PLM,2) DO JI=1,SIZE(PLM,1) ZDTHLDZ(JI,JJ,JK)= 0.5*((PTHLT(JI,JJ,JK+KKL)-PTHLT(JI,JJ,JK ))/PDZZ(JI,JJ,JK+KKL)+ & (PTHLT(JI,JJ,JK )-PTHLT(JI,JJ,JK-KKL))/PDZZ(JI,JJ,JK )) - ZVAR=XG/PTHVREF(JI,JJ,JK)*ZETHETA(JI,JJ,JK)*ZDTHLDZ(JI,JJ,JK) - ! + IF (LOCEAN) THEN + ZVAR= XG*XALPHAOC*ZDTHLDZ(JI,JJ,JK) + ELSE + ZVAR= XG/PTHVREF(JI,JJ,JK)*ZETHETA(JI,JJ,JK)*ZDTHLDZ(JI,JJ,JK) + END IF +! IF (ZVAR>0.) THEN PLM(JI,JJ,JK)=MAX(XMNH_EPSILON,MIN(PLM(JI,JJ,JK), & 0.76* SQRT(PTKET(JI,JJ,JK)/ZVAR))) @@ -2408,8 +2542,12 @@ ELSE ZDRTDZ(:,:,KKB)=0 ENDIF ! -ZWORK2D(:,:)=XG/PTHVREF(:,:,KKB)* & - (ZETHETA(:,:,KKB)*ZDTHLDZ(:,:,KKB)+ZEMOIST(:,:,KKB)*ZDRTDZ(:,:,KKB)) +IF (LOCEAN) THEN + ZWORK2D(:,:)=XG*(XALPHAOC*ZDTHLDZ(:,:,KKB)-XBETAOC*ZDRTDZ(:,:,KKB)) +ELSE + ZWORK2D(:,:)=XG/PTHVREF(:,:,KKB)* & + (ZETHETA(:,:,KKB)*ZDTHLDZ(:,:,KKB)+ZEMOIST(:,:,KKB)*ZDRTDZ(:,:,KKB)) +END IF WHERE(ZWORK2D(:,:)>0.) PLM(:,:,KKB)=MAX(XMNH_EPSILON,MIN( PLM(:,:,KKB), & 0.76* SQRT(PTKET(:,:,KKB)/ZWORK2D(:,:)))) @@ -2422,15 +2560,26 @@ IF (.NOT. ORMC01) THEN ! DO JJ=1,SIZE(PLM,2) DO JI=1,SIZE(PLM,1) - DO JK=KKTB,KKTE - ZD=ZALPHA*(0.5*(PZZ(JI,JJ,JK)+PZZ(JI,JJ,JK+KKL))-PZZ(JI,JJ,KKB)) & - *PDIRCOSZW(JI,JJ) - IF ( PLM(JI,JJ,JK)>ZD) THEN - PLM(JI,JJ,JK)=ZD - ELSE - EXIT - ENDIF - END DO + IF (LOCEAN) THEN + DO JK=KKTE,KKTB,-1 + ZD=ZALPHA*(PZZ(JI,JJ,KKTE+1)-PZZ(JI,JJ,JK)) + IF ( PLM(JI,JJ,JK)>ZD) THEN + PLM(JI,JJ,JK)=ZD + ELSE + EXIT + ENDIF + END DO + ELSE + DO JK=KKTB,KKTE + ZD=ZALPHA*(0.5*(PZZ(JI,JJ,JK)+PZZ(JI,JJ,JK+KKL))-PZZ(JI,JJ,KKB)) & + *PDIRCOSZW(JI,JJ) + IF ( PLM(JI,JJ,JK)>ZD) THEN + PLM(JI,JJ,JK)=ZD + ELSE + EXIT + ENDIF + END DO + ENDIF END DO END DO END IF diff --git a/src/ZSOLVER/turb_hor_dyn_corr.f90 b/src/ZSOLVER/turb_hor_dyn_corr.f90 index b59a3dfaa26314a1a3360157b626ff755290bc70..84f5cea6c398fa05be5283cf2085aa994c821f77 100644 --- a/src/ZSOLVER/turb_hor_dyn_corr.f90 +++ b/src/ZSOLVER/turb_hor_dyn_corr.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. @@ -8,7 +8,7 @@ MODULE MODI_TURB_HOR_DYN_CORR INTERFACE ! SUBROUTINE TURB_HOR_DYN_CORR(KSPLT, PTSTEP, & - OCLOSE_OUT,OTURB_FLX,KRR, & + OTURB_FLX,KRR, & TPFILE, & PK,PINV_PDZZ, & PDXX,PDYY,PDZZ,PDZX,PDZY,PZZ, & @@ -26,8 +26,6 @@ USE MODD_IO, ONLY: TFILEDATA ! INTEGER, INTENT(IN) :: KSPLT ! split process index REAL, INTENT(IN) :: PTSTEP ! timestep -LOGICAL, INTENT(IN) :: OCLOSE_OUT ! switch for syncronous - ! file opening LOGICAL, INTENT(IN) :: OTURB_FLX ! switch to write the ! turbulent fluxes in the syncronous FM-file INTEGER, INTENT(IN) :: KRR ! number of moist var. @@ -77,7 +75,7 @@ END INTERFACE END MODULE MODI_TURB_HOR_DYN_CORR ! ################################################################ SUBROUTINE TURB_HOR_DYN_CORR(KSPLT, PTSTEP, & - OCLOSE_OUT,OTURB_FLX,KRR, & + OTURB_FLX,KRR, & TPFILE, & PK,PINV_PDZZ, & PDXX,PDYY,PDZZ,PDZX,PDZY,PZZ, & @@ -146,16 +144,19 @@ USE MODD_ARGSLIST_ll, ONLY: LIST_ll USE MODD_CST USE MODD_CONF USE MODD_CTURB +use modd_field, only: tfielddata, TYPEREAL USE MODD_IO, ONLY: TFILEDATA USE MODD_PARAMETERS USE MODD_LES USE MODD_NSV ! USE MODE_ll -USE MODE_FIELD, ONLY: TFIELDDATA, TYPEREAL USE MODE_IO_FIELD_WRITE, only: IO_Field_write use mode_mppdb ! +#ifdef MNH_OPENACC +USE MODI_GET_HALO +#endif USE MODI_GRADIENT_M USE MODI_GRADIENT_U USE MODI_GRADIENT_V @@ -177,7 +178,6 @@ USE MODI_BITREP #ifdef MNH_OPENACC USE MODE_MNH_ZWORK , ONLY : MNH_ALLOCATE_ZT3D , MNH_REL_ZT3D, MNH_ALLOCATE_ZT3DP , MNH_ALLOCATE_ZT2D #endif -USE MODI_GET_HALO ! IMPLICIT NONE ! @@ -188,8 +188,6 @@ IMPLICIT NONE ! INTEGER, INTENT(IN) :: KSPLT ! split process index REAL, INTENT(IN) :: PTSTEP ! timestep -LOGICAL, INTENT(IN) :: OCLOSE_OUT ! switch for syncronous - ! file opening LOGICAL, INTENT(IN) :: OTURB_FLX ! switch to write the ! turbulent fluxes in the syncronous FM-file INTEGER, INTENT(IN) :: KRR ! number of moist var. @@ -237,8 +235,10 @@ REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PDP,PTP ! TKE production terms REAL, DIMENSION(:,:,:), pointer , contiguous :: ZFLX,ZWORK ! work arrays, PK is the turb. mixing coef. ! REAL, DIMENSION(:,:), pointer , contiguous :: ZDIRSINZW -INTEGER :: IZFLX,IZWORK,IZDIRSINZW ! sinus of the angle between the vertical and the normal to the orography +#ifdef MNH_OPENACC +INTEGER :: IZFLX,IZWORK,IZDIRSINZW +#endif INTEGER :: IKB,IKE ! Index values for the Beginning and End ! mass points of the domain @@ -249,7 +249,9 @@ REAL, DIMENSION(:,:,:), pointer , contiguous :: GX_U_M_PUM REAL, DIMENSION(:,:,:), pointer , contiguous :: GY_V_M_PVM REAL, DIMENSION(:,:,:), pointer , contiguous :: GZ_W_M_PWM REAL, DIMENSION(:,:,:), pointer , contiguous :: GZ_W_M_ZWP +#ifdef MNH_OPENACC INTEGER :: IGX_U_M_PUM,IGY_V_M_PVM,IGZ_W_M_PWM,IGZ_W_M_ZWP +#endif REAL, DIMENSION(:,:,:), pointer , contiguous :: ZMZF_DZZ ! MZF(PDZZ) REAL, DIMENSION(:,:,:), pointer , contiguous :: ZDFDDWDZ ! formal derivative of the ! ! flux (variable: dW/dz) @@ -260,8 +262,10 @@ REAL, DIMENSION(:,:,:), pointer , contiguous :: ZDV_DZ_DZS_DY ! dv/dz*dzs/dy sur REAL, DIMENSION(:,:,:), pointer , contiguous :: ZDU_DX ! du/dx surf REAL, DIMENSION(:,:,:), pointer , contiguous :: ZDV_DY ! dv/dy surf REAL, DIMENSION(:,:,:), pointer , contiguous :: ZDW_DZ ! dw/dz surf +#ifdef MNH_OPENACC INTEGER :: IZMZF_DZZ,IZDFDDWDZ,IZWP,IZDU_DZ_DZS_DX,IZDV_DZ_DZS_DY & ,IZDU_DX,IZDV_DY,IZDW_DZ +#endif ! INTEGER :: IINFO_ll ! return code of parallel routine TYPE(LIST_ll), POINTER :: TZFIELDS_ll ! list of fields to exchange @@ -272,7 +276,9 @@ REAL :: ZTIME1, ZTIME2 REAL, DIMENSION(:,:,:), pointer , contiguous :: ZCOEFF , ZDZZ ! coefficients for the uncentred gradient ! computation near the ground +#ifdef MNH_OPENACC INTEGER :: IZCOEFF , IZDZZ +#endif ! #ifdef MNH_OPENACC REAL, DIMENSION(:,:,:), pointer , contiguous :: ZTMP1_DEVICE, ZTMP2_DEVICE, ZTMP3_DEVICE, ZTMP4_DEVICE @@ -660,12 +666,16 @@ ZFLX(:,:,IKB-1) = 2. * ZFLX(:,:,IKB-1) - ZFLX(:,:,IKB) !!! to be absolutely sure, we do a wait !$acc wait ! +#ifndef MNH_OPENACC +CALL UPDATE_HALO_ll(TZFIELDS_ll, IINFO_ll) +#else !!$!$acc update self(ZFLX) !!$CALL UPDATE_HALO_ll(TZFIELDS_ll, IINFO_ll) !!$!$acc update device(ZFLX) async(10) CALL GET_HALO_D(ZFLX) +#endif ! -IF ( OCLOSE_OUT .AND. OTURB_FLX ) THEN +IF ( tpfile%lopened .AND. OTURB_FLX ) THEN ! stores <U U> TZFIELD%CMNHNAME = 'U_VAR' TZFIELD%CSTDNAME = '' @@ -869,12 +879,16 @@ ZFLX(:,:,IKB-1) = 2. * ZFLX(:,:,IKB-1) - ZFLX(:,:,IKB) ! ! !$acc wait(3) ! !$acc wait(3) +#ifndef MNH_OPENACC +CALL UPDATE_HALO_ll(TZFIELDS_ll, IINFO_ll) +#else !!$!$acc update self(ZFLX) !!$CALL UPDATE_HALO_ll(TZFIELDS_ll, IINFO_ll) !!$!$acc update device(ZFLX) async(10) CALL GET_HALO_D(ZFLX) +#endif ! -IF ( OCLOSE_OUT .AND. OTURB_FLX ) THEN +IF ( tpfile%lopened .AND. OTURB_FLX ) THEN ! stores <V V> TZFIELD%CMNHNAME = 'V_VAR' TZFIELD%CSTDNAME = '' @@ -1079,7 +1093,7 @@ END DO ! CONCURRENT ZFLX(:,:,IKB-1) = 2. * ZFLX(:,:,IKB-1) - ZFLX(:,:,IKB) !$acc end kernels ! -IF ( OCLOSE_OUT .AND. OTURB_FLX ) THEN +IF ( tpfile%lopened .AND. OTURB_FLX ) THEN !$acc wait(3) !$acc update self(ZFLX) ! stores <W W> diff --git a/src/ZSOLVER/turb_hor_thermo_flux.f90 b/src/ZSOLVER/turb_hor_thermo_flux.f90 index 76a1d60c284b3896f10722e6bf4f6a0112c2f37c..9c02d72a17d3f2d5136e0cec6a6b92a342789d4e 100644 --- a/src/ZSOLVER/turb_hor_thermo_flux.f90 +++ b/src/ZSOLVER/turb_hor_thermo_flux.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. @@ -10,7 +10,7 @@ INTERFACE ! SUBROUTINE TURB_HOR_THERMO_FLUX(KSPLT, KRR, KRRL, KRRI, & - OCLOSE_OUT,OTURB_FLX,OSUBG_COND, & + OTURB_FLX,OSUBG_COND, & TPFILE, & PK,PINV_PDXX,PINV_PDYY,PINV_PDZZ,PMZM_PRHODJ, & PDXX,PDYY,PDZZ,PDZX,PDZY, & @@ -27,8 +27,6 @@ INTEGER, INTENT(IN) :: KSPLT ! split process index INTEGER, INTENT(IN) :: KRR ! number of moist var. INTEGER, INTENT(IN) :: KRRL ! number of liquid water var. INTEGER, INTENT(IN) :: KRRI ! number of ice water var. -LOGICAL, INTENT(IN) :: OCLOSE_OUT ! switch for syncronous - ! file opening LOGICAL, INTENT(IN) :: OTURB_FLX ! switch to write the ! turbulent fluxes in the syncronous FM-file LOGICAL, INTENT(IN) :: OSUBG_COND ! Switch for sub-grid @@ -74,7 +72,7 @@ END INTERFACE END MODULE MODI_TURB_HOR_THERMO_FLUX ! ################################################################ SUBROUTINE TURB_HOR_THERMO_FLUX(KSPLT, KRR, KRRL, KRRI, & - OCLOSE_OUT,OTURB_FLX,OSUBG_COND, & + OTURB_FLX,OSUBG_COND, & TPFILE, & PK,PINV_PDXX,PINV_PDYY,PINV_PDZZ,PMZM_PRHODJ, & PDXX,PDYY,PDZZ,PDZX,PDZY, & @@ -132,11 +130,11 @@ END MODULE MODI_TURB_HOR_THERMO_FLUX USE MODD_CST USE MODD_CONF USE MODD_CTURB +use modd_field, only: tfielddata, TYPEREAL USE MODD_IO, ONLY: TFILEDATA USE MODD_PARAMETERS USE MODD_LES ! -USE MODE_FIELD, ONLY: TFIELDDATA, TYPEREAL USE MODE_IO_FIELD_WRITE, only: IO_Field_write use mode_mppdb ! @@ -170,8 +168,6 @@ INTEGER, INTENT(IN) :: KSPLT ! split process index INTEGER, INTENT(IN) :: KRR ! number of moist var. INTEGER, INTENT(IN) :: KRRL ! number of liquid water var. INTEGER, INTENT(IN) :: KRRI ! number of ice water var. -LOGICAL, INTENT(IN) :: OCLOSE_OUT ! switch for syncronous - ! file opening LOGICAL, INTENT(IN) :: OTURB_FLX ! switch to write the ! turbulent fluxes in the syncronous FM-file LOGICAL, INTENT(IN) :: OSUBG_COND ! Switch for sub-grid @@ -221,7 +217,9 @@ INTEGER :: IKB,IKE,IKU REAL, DIMENSION(:,:,:), pointer , contiguous :: ZCOEFF ! coefficients for the uncentred gradient ! computation near the ground +#ifdef MNH_OPENACC INTEGER :: IZFLX,IZFLXC,IZCOEFF +#endif ! REAL :: ZTIME1, ZTIME2 ! @@ -592,7 +590,7 @@ END IF !!ZWORK(:,:,:) = ZFLX(:,:,:) ! ! stores the horizontal <U THl> -IF ( OCLOSE_OUT .AND. OTURB_FLX ) THEN +IF ( tpfile%lopened .AND. OTURB_FLX ) THEN TZFIELD%CMNHNAME = 'UTHL_FLX' TZFIELD%CSTDNAME = '' TZFIELD%CLONGNAME = 'UTHL_FLX' @@ -735,7 +733,7 @@ IF (KRR/=0) THEN END IF ! ! stores the horizontal <U Rnp> - IF ( OCLOSE_OUT .AND. OTURB_FLX ) THEN + IF ( tpfile%lopened .AND. OTURB_FLX ) THEN TZFIELD%CMNHNAME = 'UR_FLX' TZFIELD%CSTDNAME = '' TZFIELD%CLONGNAME = 'UR_FLX' @@ -928,7 +926,7 @@ END DO END IF ! ! stores the horizontal <U Rnp> - IF ( OCLOSE_OUT .AND. OTURB_FLX ) THEN + IF ( tpfile%lopened .AND. OTURB_FLX ) THEN !$acc update self(ZFLX) TZFIELD%CMNHNAME = 'UR_FLX' TZFIELD%CSTDNAME = '' @@ -1011,7 +1009,7 @@ END IF !! ZFLX(:,:,:)*MXM(EMOIST(KRR,KRRI,PTHLT,PEXNREF,PRT,PLOCPT,PSRCM)) !! ! !! ! stores the horizontal <U VPT> -!! IF ( OCLOSE_OUT .AND. OTURB_FLX ) THEN +!! IF ( tpfile%lopened .AND. OTURB_FLX ) THEN !! TZFIELD%CMNHNAME = 'UVPT_FLX' !! TZFIELD%CSTDNAME = '' !! TZFIELD%CLONGNAME = 'UVPT_FLX' @@ -1118,7 +1116,7 @@ END IF !!ZWORK(:,:,:) = ZFLX(:,:,:) ! ! stores the horizontal <V THl> -IF ( OCLOSE_OUT .AND. OTURB_FLX ) THEN +IF ( tpfile%lopened .AND. OTURB_FLX ) THEN TZFIELD%CMNHNAME = 'VTHL_FLX' TZFIELD%CSTDNAME = '' TZFIELD%CLONGNAME = 'VTHL_FLX' @@ -1330,7 +1328,7 @@ END IF !!ZWORK(:,:,:) = ZFLX(:,:,:) ! ! stores the horizontal <V THl> -IF ( OCLOSE_OUT .AND. OTURB_FLX ) THEN +IF ( tpfile%lopened .AND. OTURB_FLX ) THEN !$acc update self(ZFLX) TZFIELD%CMNHNAME = 'VTHL_FLX' TZFIELD%CSTDNAME = '' @@ -1474,7 +1472,7 @@ IF (KRR/=0) THEN END IF ! ! stores the horizontal <V Rnp> - IF ( OCLOSE_OUT .AND. OTURB_FLX ) THEN + IF ( tpfile%lopened .AND. OTURB_FLX ) THEN TZFIELD%CMNHNAME = 'VR_FLX' TZFIELD%CSTDNAME = '' TZFIELD%CLONGNAME = 'VR_FLX' @@ -1687,7 +1685,7 @@ IF (KRR/=0) THEN END IF ! ! stores the horizontal <V Rnp> - IF ( OCLOSE_OUT .AND. OTURB_FLX ) THEN + IF ( tpfile%lopened .AND. OTURB_FLX ) THEN !$acc update self(ZFLX) TZFIELD%CMNHNAME = 'VR_FLX' TZFIELD%CSTDNAME = '' @@ -1774,7 +1772,7 @@ END IF !! END IF !! ! !! ! stores the horizontal <V VPT> -!! IF ( OCLOSE_OUT .AND. OTURB_FLX ) THEN +!! IF ( tpfile%lopened .AND. OTURB_FLX ) THEN !! TZFIELD%CMNHNAME = 'VVPT_FLX' !! TZFIELD%CSTDNAME = '' !! TZFIELD%CLONGNAME = 'VVPT_FLX' diff --git a/src/ZSOLVER/update_lm.f90 b/src/ZSOLVER/update_lm.f90 index 7b82c64c7d766f33f34af358abd3c5c5cb3b37ae..1d5b5069faac3d60df15af68f0a3886307b61749 100644 --- a/src/ZSOLVER/update_lm.f90 +++ b/src/ZSOLVER/update_lm.f90 @@ -66,6 +66,7 @@ USE MODD_PARAMETERS ! USE MODE_ll use mode_mppdb +! USE MODI_GET_HALO ! IMPLICIT NONE @@ -111,19 +112,22 @@ CALL GET_INDICE_ll (IIB,IJB,IIE,IJE) ! ------------- ! ! +#ifndef MNH_OPENACC !!$IF(NHALO == 1) THEN !!$!$acc update self(PLM,PLEPS) -!!$ CALL ADD3DFIELD_ll( TZLM_ll, PLM, 'UPDATE_LM::PLM' ) -!!$ CALL ADD3DFIELD_ll( TZLM_ll, PLEPS, 'UPDATE_LM::PLEPS' ) -!!$ CALL UPDATE_HALO_ll(TZLM_ll,IINFO_ll) -!!$ CALL CLEANLIST_ll(TZLM_ll) + CALL ADD3DFIELD_ll( TZLM_ll, PLM, 'UPDATE_LM::PLM' ) + CALL ADD3DFIELD_ll( TZLM_ll, PLEPS, 'UPDATE_LM::PLEPS' ) + CALL UPDATE_HALO_ll(TZLM_ll,IINFO_ll) + CALL CLEANLIST_ll(TZLM_ll) !!$!$acc update device(PLM,PLEPS) +!!$END IF +#else ! ! /!\ Corner update needed ! -> GET_HALO_DDC ! CALL GET_HALO_DDC( PLM, HNAME='UPDATE_LM::PLM' ) CALL GET_HALO_DDC( PLEPS, HNAME='UPDATE_LM::PLEPS' ) -!!$END IF +#endif ! !------------------------------------------------------------------------------- ! diff --git a/src/ZSOLVER/zsolver_inv.f90 b/src/ZSOLVER/zsolver_inv.f90 index 29965951e8ee7546dad9e2e27d6f93e6a1cd3e1e..4c85d3aecad9394aaf80c7d0df8836cdd97d86cb 100644 --- a/src/ZSOLVER/zsolver_inv.f90 +++ b/src/ZSOLVER/zsolver_inv.f90 @@ -312,7 +312,11 @@ SUBROUTINE ZSOLVER_INV(HLBCX,HLBCY,PDXHATM,PDYHATM,PRHOM,PAF,PBF,PCF, & ! ! WARNING WITH GET_HALO_D not BITREPROD !!! ! +#ifndef MNH_OPENACC + CALL GET_HALO( PF_1_Y ) +#else CALL GET_HALO_DDC(PF_1_Y) +#endif ! CALL PF_1_Y_BOUND(PF_1_Y) !-------------------------------------------------------------------------------