diff --git a/src/PHYEX/turb/mode_tke_eps_sources.f90 b/src/PHYEX/turb/mode_tke_eps_sources.f90
index f7cb314bc30fa0fa2ccbf4b15d6b9b8b34d8b3ce..e9fc927b002ae295ba978c402c26f7cb7e3c7455 100644
--- a/src/PHYEX/turb/mode_tke_eps_sources.f90
+++ b/src/PHYEX/turb/mode_tke_eps_sources.f90
@@ -238,6 +238,7 @@ IKT=D%NKT
IKA=D%NKA
IKL=D%NKL
!
+!$acc kernels
! compute the effective diffusion coefficient at the mass point
!$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
ZKEFF(:,:) = PLM(:,:) * SQRT(PTKEM(:,:))
@@ -278,9 +279,6 @@ END IF
! Compute the source terms for TKE: ( ADVECtion + NUMerical DIFFusion + ..)
! + (Dynamical Production) + (Thermal Production) - (dissipation)
!
-CALL MZM_PHY(D,ZKEFF,ZMWORK1)
-CALL MZM_PHY(D,PRHODJ,ZMWORK2)
-!
!$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
ZFLX(:,:) = CSTURB%XCED * SQRT(PTKEM(:,:)) / PLEPS(:,:)
ZSOURCE(:,:) = ( PRTKES(:,:) + PRTKEMS(:,:) ) &
@@ -300,18 +298,18 @@ IF (OOCEAN) THEN
ZSOURCE(:,IKE)=ZSOURCE(:,IKE)-1.E2*((PSFUM(:)**2 + PSFVM(:)**2)**1.5) /PDZZ(:,IKE)
!$mnh_end_expand_array(JIJ=IIJB:IIJE)
END IF
+!$acc end kernels
! Compute the vector giving the elements just under the diagonal for the
! matrix inverted in TRIDIAG
!
-!$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
-ZA(:,:) = - PTSTEP * CSTURB%XCET * ZMWORK1(:,:) &
- * ZMWORK2(:,:) / PDZZ(:,:)**2
-!$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
+!
+ZA(:,:) = - PTSTEP * CSTURB%XCET * MZM(ZKEFF) * MZM(PRHODJ) / PDZZ(:,:)**2
!
! Compute TKE at time t+deltat: ( stored in ZRES )
!
CALL TRIDIAG_TKE(D,PTKEM,ZA,PTSTEP,PEXPL,TURBN%XIMPL,PRHODJ,ZSOURCE,PTSTEP*ZFLX,ZRES)
CALL GET_HALO_PHY(D,ZRES)
+!$acc update device(ZRES)
!
!* diagnose the dissipation
!
@@ -337,41 +335,29 @@ IF(HPROGRAM/='MESONH') THEN
!$mnh_end_expand_where(JIJ=IIJB:IIJE,JK=1:IKT)
END IF
!
+!$acc kernels
!$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
PTDISS(:,:) = - ZFLX(:,:)*(PEXPL*PTKEM(:,:) &
+ TURBN%XIMPL*ZRES(:,:))
!$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
+!$acc end kernels
!
IF ( TLES%LLES_CALL .OR. &
(TURBN%LTURB_DIAG .AND. TPFILE%LOPENED) ) THEN
!
! Compute the cartesian vertical flux of TKE in ZFLX
!
- CALL MZM_PHY(D,ZKEFF,ZMWORK1)
- !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
- ZDWORK1(:,:) = TURBN%XIMPL * ZRES(:,:) + PEXPL * PTKEM(:,:)
- !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
- CALL DZM_PHY(D,ZDWORK1,ZDWORK2)
- !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
- ZFLX(:,:) = - CSTURB%XCET * ZMWORK1(:,:) &
- * ZDWORK2(:,:) / PDZZ(:,:)
- !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
+ ZFLX(:,:) = - CSTURB%XCET * MZM(ZKEFF(:,:)) * &
+ DZM(TURBN%XIMPL * ZRES(:,:) + PEXPL * PTKEM(:,:) ) / PDZZ(:,:)
!
+!$acc kernels
ZFLX(:,IKB) = 0.
ZFLX(:,IKA) = 0.
+!$acc end kernels
!
! Compute the whole turbulent TRansport of TKE:
!
- CALL MZM_PHY(D,PRHODJ,ZMWORK1)
- !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
- ZMWORK2(:,:) = ZMWORK1(:,:) * ZFLX(:,:) &
- / PDZZ(:,:)
- !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
- CALL DZF_PHY(D,ZMWORK2,ZDWORK1)
- !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
- ZTR(:,:)= ZTR(:,:) - ZDWORK1(:,:) &
- /PRHODJ(:,:)
- !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
+ PTR(:,:)= PTR(:,:) - DZF( MZM(PRHODJ(:,:)) * ZFLX(:,:) / PDZZ(:,:) ) /PRHODJ(:,:)
!
! Storage in the LES configuration
!
@@ -387,22 +373,28 @@ END IF
!
IF (BUCONF%LBUDGET_TKE) THEN
! Dynamical production
+ !$acc kernels
!$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
ZMWORK1(:,:) = PDP(:,:) * PRHODJ(:,:)
!$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
+ !$acc end kernels
CALL BUDGET_STORE_ADD_PHY(D, TBUDGETS(NBUDGET_TKE), 'DP', ZMWORK1)
!
! Thermal production
+ !$acc kernels
!$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
ZMWORK1(:,:) = PTP(:,:) * PRHODJ(:,:)
!$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
+ !$acc end kernels
CALL BUDGET_STORE_ADD_PHY(D, TBUDGETS(NBUDGET_TKE), 'TP', ZMWORK1)
!
! Dissipation
+ !$acc kernels
!$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
ZMWORK1(:,:) = -CSTURB%XCED * SQRT(PTKEM(:,:))/PLEPS(:,:) * &
(PEXPL*PTKEM(:,:) + TURBN%XIMPL*ZRES(:,:))*PRHODJ(:,:)
!$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
+ !$acc end kernels
CALL BUDGET_STORE_ADD_PHY(D, TBUDGETS(NBUDGET_TKE), 'DISS',ZMWORK1)
END IF
!
@@ -410,6 +402,7 @@ END IF
! with the removal of the advection part for MesoNH
!Should be in IF LBUDGET_TKE only. Was removed out for a correct comput. of PTDIFF in case of LBUDGET_TKE=F in AROME
+!$acc kernels present_cr(ZRES)
!$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
#ifdef REPRO48
IF (BUCONF%LBUDGET_TKE) THEN
@@ -429,13 +422,16 @@ PTDIFF(:,:) = ZRES(:,:) / PTSTEP - PRTKES(:,:)&
/PRHODJ(:,:) &
& - PDP(:,:)- PTP(:,:) - PTDISS(:,:)
!$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
+!$acc end kernels
!
IF (BUCONF%LBUDGET_TKE) CALL BUDGET_STORE_INIT_PHY(D, TBUDGETS(NBUDGET_TKE), 'TR', PRTKES)
!
+!$acc kernels
!$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
PRTKES(:,:) = ZRES(:,:) * PRHODJ(:,:) / PTSTEP &
- PRTKEMS(:,:)
!$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
+!$acc end kernels
!
! stores the whole turbulent transport
!
@@ -446,31 +442,44 @@ IF (BUCONF%LBUDGET_TKE) CALL BUDGET_STORE_END_PHY(D, TBUDGETS(NBUDGET_TKE), 'TR'
!* 3. COMPUTE THE DISSIPATIVE HEATING
! -------------------------------
!
+!$acc kernels
!$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
PRTHLS(:,:) = PRTHLS(:,:) + &
CSTURB%XCED * SQRT(PTKEM(:,:)) / PLEPS(:,:) * &
(PEXPL*PTKEM(:,:) + TURBN%XIMPL*ZRES(:,:)) &
* PRHODJ(:,:) * PCOEF_DISS(:,:)
!$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
+!$acc end kernels
!----------------------------------------------------------------------------
!
!* 4. STORES SOME DIAGNOSTICS
! -----------------------
!
-IF(PRESENT(PTR)) PTR=ZTR
+IF(PRESENT(PTR)) THEN
+!$acc kernels
+ !$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
+ PTR(:,:)=ZTR(:,:)
+ !$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
+!$acc end kernels
+END IF
IF(PRESENT(PDISS)) THEN
+!$acc kernels
!$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
PDISS(:,:) = -CSTURB%XCED * (PTKEM(:,:)**1.5) / PLEPS(:,:)
!$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
+!$acc end kernels
END IF
!
IF(PRESENT(PEDR)) THEN
+!$acc kernels
!$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
PEDR(:,:) = CSTURB%XCED * (PTKEM(:,:)**1.5) / PLEPS(:,:)
!$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
+!$acc end kernels
END IF
!
IF ( TURBN%LTURB_DIAG .AND. TPFILE%LOPENED ) THEN
+!$acc update self(PDP,PTP,PTR,PDISS)
!
! stores the dynamic production
!
diff --git a/src/PHYEX/turb/mode_turb_ver.f90 b/src/PHYEX/turb/mode_turb_ver.f90
index 0bd0644af7dbb3c72d7830a30f6ca7543aafd1f1..feb115323b9230a73704e75bfd770851c261d3f2 100644
--- a/src/PHYEX/turb/mode_turb_ver.f90
+++ b/src/PHYEX/turb/mode_turb_ver.f90
@@ -203,6 +203,7 @@ SUBROUTINE TURB_VER(D,CST,CSTURB,TURBN,TLES,KRR,KRRL,KRRI,KGRADIENTS,&
!! 10/2012 (J.Escobar) Bypass PGI bug , redefine some allocatable array inplace of automatic
!! 08/2014 (J.Escobar) Bypass PGI memory leak bug , replace IF statement with IF THEN ENDIF
!! Modifications: July, 2015 (Wim de Rooy) switch for HARATU (Racmo turbulence scheme)
+!! 04/2016 (M.Moge) Use openACC directives to port the TURB part of Meso-NH on GPU
!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
!! JL Redelsperger 03/2021 : add Ocean LES case
!!--------------------------------------------------------------------------
@@ -422,6 +423,7 @@ CALL PRANDTL(D,CST,CSTURB,KRR,KSV,KRRI,TURBN%LTURB_FLX, &
!
! Buoyancy coefficient
!
+!$acc kernels
IF (OOCEAN) THEN
!$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
ZBETA(:,:) = CST%XG*CST%XALPHAOC
@@ -437,16 +439,20 @@ END IF
!$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
ZSQRT_TKE(:,:) = SQRT(PTKEM(:,:))
!$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
+!$acc end kernels
!
! gradients of mean quantities at previous time-step
!
CALL GZ_M_W_PHY(D,PTHLM,PDZZ,ZDTH_DZ)
+!$acc kernels
ZDR_DZ(:,:) = 0.
+!$acc end kernels
IF (KRR>0) CALL GZ_M_W_PHY(D,PRM(:,:,1),PDZZ,ZDR_DZ)
!
!
! Denominator factor in 3rd order terms
!
+!$acc kernels
IF (.NOT. TURBN%LHARAT) THEN
!$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
ZD(:,:) = (1.+ZREDTH1(:,:)+ZREDR1(:,:)) * &
@@ -457,6 +463,7 @@ ELSE
ZD(:,:) = 1.
!$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
ENDIF
+!$acc end kernels
!
! Phi3 and Psi3 Prandtl numbers
!
@@ -495,12 +502,13 @@ END IF
!* 4. TURBULENT CORRELATIONS : <w Rc>, <THl THl>, <THl Rnp>, <Rnp Rnp>
! ----------------------------------------------------------------
!
-
+!$acc kernels
IF (TURBN%LHARAT) THEN
ZLM(:,:)=PLENGTHH(:,:)
ELSE
ZLM(:,:)=PLM(:,:)
ENDIF
+!$acc end kernels
!
CALL TURB_VER_THERMO_FLUX(D,CST,CSTURB,TURBN,TLES, &
KRR,KRRL,KRRI,KSV,KGRADIENTS, &
@@ -626,6 +634,7 @@ IF ( TURBN%LTURB_FLX .AND. TPFILE%LOPENED .AND. .NOT. TURBN%LHARAT) THEN
NTYPE = TYPEREAL, &
NDIMS = 3, &
LTIMEDEP = .TRUE. )
+ !$acc update self(ZPHI3)
CALL IO_FIELD_WRITE_PHY(D,TPFILE,TZFIELD,ZPHI3)
!
! stores the Turbulent Schmidt number
@@ -641,6 +650,7 @@ IF ( TURBN%LTURB_FLX .AND. TPFILE%LOPENED .AND. .NOT. TURBN%LHARAT) THEN
NTYPE = TYPEREAL, &
NDIMS = 3, &
LTIMEDEP = .TRUE. )
+ !$acc update self(ZPSI3)
CALL IO_FIELD_WRITE_PHY(D,TPFILE,TZFIELD,ZPSI3)
!
!
@@ -655,6 +665,7 @@ IF ( TURBN%LTURB_FLX .AND. TPFILE%LOPENED .AND. .NOT. TURBN%LHARAT) THEN
NTYPE = TYPEREAL, &
NDIMS = 3, &
LTIMEDEP = .TRUE. )
+ !$acc update self(ZPSI_SV)
DO JSV=1,KSV
WRITE(TZFIELD%CMNHNAME, '("PSI_SV_",I3.3)') JSV
TZFIELD%CLONGNAME = TRIM(TZFIELD%CMNHNAME)
diff --git a/src/Rules.LXnvhpc2202.mk b/src/Rules.LXnvhpc2202.mk
index 496b4e005bc8a7b8003670f58f116843ff1668b0..861f8953acf3bfa664d0b9171c7eec662558d0a5 100644
--- a/src/Rules.LXnvhpc2202.mk
+++ b/src/Rules.LXnvhpc2202.mk
@@ -334,6 +334,12 @@ $(OBJS_I4) : OPT = $(OPT_BASE_I4)
endif
SPLL = spll_new
+
+PHYEX_OPTDEFAULT = --addMPPDB_CHECKS --addStack MESONH --stopScopes toto
+
PHYEX_LIST = $(notdir $(shell find PHYEX/micro PHYEX/turb -follow -type f -name "*.f*" | sed -e 's/\(.*\)\(\.\).*/\1.D/g' ))
-$(PHYEX_LIST) : PYFT = pyft_tool.py --addMPPDB_CHECKS --shumanFUNCtoCALL --addStack MESONH --stopScopes toto
+$(PHYEX_LIST) : PYFT = pyft_tool.py $(PHYEX_OPTDEFAULT)
+
+PHYEX_SHUMAN = mode_tke_eps_sources.D
+$(PHYEX_SHUMAN) : PYFT = pyft_tool.py --shumanFUNCtoCALL $(PHYEX_OPTDEFAULT)