From d1ce800868c4c01027e892acf9435fff0155026e Mon Sep 17 00:00:00 2001
From: Quentin Rodier <quentin.rodier@meteo.fr>
Date: Tue, 29 Mar 2022 18:09:27 +0200
Subject: [PATCH] Quentin 29/03/2022: Expand Array mode_tke_eps_sources.F90
---
src/common/turb/mode_tke_eps_sources.F90 | 153 ++++++++++++++++++-----
1 file changed, 123 insertions(+), 30 deletions(-)
diff --git a/src/common/turb/mode_tke_eps_sources.F90 b/src/common/turb/mode_tke_eps_sources.F90
index 0478fea7f..3d9991f8d 100644
--- a/src/common/turb/mode_tke_eps_sources.F90
+++ b/src/common/turb/mode_tke_eps_sources.F90
@@ -213,7 +213,10 @@ REAL, DIMENSION(D%NIT,D%NJT,D%NKT) :: &
ZFLX, & ! horizontal or vertical flux of the treated variable
ZSOURCE, & ! source of evolution for the treated variable
ZKEFF, & ! effectif diffusion coeff = LT * SQRT( TKE )
- ZTR ! Transport term
+ ZTR, & ! Transport term
+ ZMWORK1,ZMWORK2,& ! working var. for MZM/MZF operators (array syntax)
+ ZDWORK1,ZDWORK2 ! working var. for DZM/DZF operators (array syntax)
+
LOGICAL,DIMENSION(D%NIT,D%NJT,D%NKT) :: GTKENEG
! 3D mask .T. if TKE < CSTURB%XTKEMIN
INTEGER :: IKB ! Index values for the Beginning and End
@@ -223,6 +226,7 @@ TYPE(LIST_ll), POINTER :: TZFIELDDISS_ll ! list of fields to exchange
INTEGER :: IINFO_ll ! return code of parallel routine
TYPE(TFIELDDATA) :: TZFIELD
REAL(KIND=JPRB) :: ZHOOK_HANDLE
+INTEGER :: JI,JJ,JK
!
!----------------------------------------------------------------------------
NULLIFY(TZFIELDDISS_ll)
@@ -235,7 +239,13 @@ IF (LHOOK) CALL DR_HOOK('TKE_EPS_SOURCES',0,ZHOOK_HANDLE)
IKB=D%NKB
!
! compute the effective diffusion coefficient at the mass point
-ZKEFF(:,:,:) = PLM(:,:,:) * SQRT(PTKEM(:,:,:))
+DO JK=1,D%NKT
+ DO JJ=1,D%NJT
+ DO JI=1,D%NIT
+ZKEFF(JI,JJ,JK) = PLM(JI,JJ,JK) * SQRT(PTKEM(JI,JJ,JK))
+ ENDDO
+ ENDDO
+ENDDO
!
!----------------------------------------------------------------------------
!
@@ -258,14 +268,26 @@ END IF
!
! extrapolate the dynamic production with a 1/Z law from its value at the
! W(IKB+1) value stored in PDP(IKB) to the mass localization tke(IKB)
-PDP(:,:,IKB) = PDP(:,:,IKB) * (1. + PDZZ(:,:,IKB+D%NKL)/PDZZ(:,:,IKB))
+!
+DO JJ=1,D%NJT
+ DO JI=1,D%NIT
+PDP(JI,JJ,IKB) = PDP(JI,JJ,IKB) * (1. + PDZZ(JI,JJ,IKB+D%NKL)/PDZZ(JI,JJ,IKB))
+ ENDDO
+ENDDO
!
! Compute the source terms for TKE: ( ADVECtion + NUMerical DIFFusion + ..)
! + (Dynamical Production) + (Thermal Production) - (dissipation)
-ZFLX(:,:,:) = CSTURB%XCED * SQRT(PTKEM(:,:,:)) / PLEPS(:,:,:)
-ZSOURCE(:,:,:) = ( PRTKES(:,:,:) + PRTKEMS(:,:,:) ) / PRHODJ(:,:,:) &
- - PTKEM(:,:,:) / PTSTEP &
- + PDP(:,:,:) + PTP(:,:,:) + ZTR(:,:,:) - PEXPL * ZFLX(:,:,:) * PTKEM(:,:,:)
+!
+ZMWORK1 = MZM(ZKEFF, D%NKA, D%NKU, D%NKL)
+ZMWORK2 = MZM(PRHODJ, D%NKA, D%NKU, D%NKL)
+!
+DO JK=1,D%NKT
+ DO JJ=1,D%NJT
+ DO JI=1,D%NIT
+ZFLX(JI,JJ,JK) = CSTURB%XCED * SQRT(PTKEM(JI,JJ,JK)) / PLEPS(JI,JJ,JK)
+ZSOURCE(JI,JJ,JK) = ( PRTKES(JI,JJ,JK) + PRTKEMS(JI,JJ,JK) ) / PRHODJ(JI,JJ,JK) &
+ - PTKEM(JI,JJ,JK) / PTSTEP &
+ + PDP(JI,JJ,JK) + PTP(JI,JJ,JK) + ZTR(JI,JJ,JK) - PEXPL * ZFLX(JI,JJ,JK) * PTKEM(JI,JJ,JK)
!
!* 2.2 implicit vertical TKE transport
!
@@ -273,8 +295,10 @@ ZSOURCE(:,:,:) = ( PRTKES(:,:,:) + PRTKEMS(:,:,:) ) / PRHODJ(:,:,:) &
! Compute the vector giving the elements just under the diagonal for the
! matrix inverted in TRIDIAG
!
-ZA(:,:,:) = - PTSTEP * CSTURB%XCET * &
- MZM(ZKEFF, D%NKA, D%NKU, D%NKL) * MZM(PRHODJ, D%NKA, D%NKU, D%NKL) / PDZZ**2
+ZA(JI,JJ,JK) = - PTSTEP * CSTURB%XCET * ZMWORK1(JI,JJ,JK) * ZMWORK2(JI,JJ,JK) / PDZZ(JI,JJ,JK)**2
+ ENDDO
+ ENDDO
+ENDDO
!
! Compute TKE at time t+deltat: ( stored in ZRES )
!
@@ -285,8 +309,15 @@ CALL GET_HALO(ZRES)
!* diagnose the dissipation
!
IF (LDIAG_IN_RUN) THEN
- XCURRENT_TKE_DISS = ZFLX(:,:,:) * PTKEM(:,:,:) &
- *(PEXPL*PTKEM(:,:,:) + PIMPL*ZRES(:,:,:))
+ DO JK=1,D%NKT
+ DO JJ=1,D%NJT
+ DO JI=1,D%NIT
+ XCURRENT_TKE_DISS(JI,JJ,JK) = ZFLX(JI,JJ,JK) * PTKEM(JI,JJ,JK) &
+ *(PEXPL*PTKEM(JI,JJ,JK) + PIMPL*ZRES(JI,JJ,JK))
+ ENDDO
+ ENDDO
+ENDDO
+!
CALL ADD3DFIELD_ll( TZFIELDDISS_ll, XCURRENT_TKE_DISS, 'TKE_EPS_SOURCES::XCURRENT_TKE_DISS' )
CALL UPDATE_HALO_ll(TZFIELDDISS_ll,IINFO_ll)
CALL CLEANLIST_ll(TZFIELDDISS_ll)
@@ -295,27 +326,54 @@ ENDIF
! TKE must be greater than its minimum value
! CL : Now done at the end of the time step in ADVECTION_METSV for MesoNH
IF(HPROGRAM/='MESONH') THEN
- GTKENEG = ZRES <= CSTURB%XTKEMIN
- WHERE ( GTKENEG )
- ZRES = CSTURB%XTKEMIN
- END WHERE
+ GTKENEG = ZRES <= CSTURB%XTKEMIN
+ DO JK=1,D%NKT
+ DO JJ=1,D%NJT
+ DO JI=1,D%NIT
+ IF ( GTKENEG(JI,JJ,JK) ) THEN
+ ZRES(JI,JJ,JK) = CSTURB%XTKEMIN
+ ENDIF
+ ENDDO
+ ENDDO
+ENDDO
END IF
-PTDISS(:,:,:) = - ZFLX(:,:,:)*(PEXPL*PTKEM(:,:,:) + PIMPL*ZRES(:,:,:))
+!
+DO JK=1,D%NKT
+ DO JJ=1,D%NJT
+ DO JI=1,D%NIT
+PTDISS(JI,JJ,JK) = - ZFLX(JI,JJ,JK)*(PEXPL*PTKEM(JI,JJ,JK) + PIMPL*ZRES(JI,JJ,JK))
+ ENDDO
+ ENDDO
+ENDDO
!
IF ( OLES_CALL .OR. &
(OTURB_DIAG .AND. TPFILE%LOPENED) ) THEN
!
! Compute the cartesian vertical flux of TKE in ZFLX
!
- ZFLX(:,:,:) = - CSTURB%XCET * MZM(ZKEFF, D%NKA, D%NKU, D%NKL) * &
- DZM(PIMPL * ZRES + PEXPL * PTKEM, D%NKA, D%NKU, D%NKL) / PDZZ
+ ZMWORK1 = MZM(ZKEFF, D%NKA, D%NKU, D%NKL)
+ ZDWORK1 = DZM(PIMPL * ZRES + PEXPL * PTKEM, D%NKA, D%NKU, D%NKL)
+ DO JK=1,D%NKT
+ DO JJ=1,D%NJT
+ DO JI=1,D%NIT
+ ZFLX(JI,JJ,JK) = - CSTURB%XCET * ZMWORK1(JI,JJ,JK) * ZDWORK1(JI,JJ,JK) / PDZZ(JI,JJ,JK)
+ ENDDO
+ ENDDO
+ENDDO
!
ZFLX(:,:,IKB) = 0.
ZFLX(:,:,D%NKA) = 0.
!
! Compute the whole turbulent TRansport of TKE:
!
- ZTR(:,:,:)= ZTR - DZF(MZM(PRHODJ, D%NKA, D%NKU, D%NKL) * ZFLX / PDZZ, D%NKA, D%NKU, D%NKL) /PRHODJ
+ ZDWORK1 = DZF(MZM(PRHODJ, D%NKA, D%NKU, D%NKL) * ZFLX / PDZZ, D%NKA, D%NKU, D%NKL)
+ DO JK=1,D%NKT
+ DO JJ=1,D%NJT
+ DO JI=1,D%NIT
+ ZTR(JI,JJ,JK)= ZTR(JI,JJ,JK) - ZDWORK1(JI,JJ,JK) /PRHODJ(JI,JJ,JK)
+ ENDDO
+ ENDDO
+ENDDO
!
! Storage in the LES configuration
!
@@ -343,16 +401,28 @@ END IF
!Store the previous source terms in prtkes before initializing the next one
!Should be in IF LBUDGET_TKE only. Was removed out for a correct comput. of PTDIFF in case of LBUDGET_TKE=F in AROME
-PRTKES(:,:,:) = PRTKES(:,:,:) + PRHODJ(:,:,:) * &
- ( PDP(:,:,:) + PTP(:,:,:) &
- - CSTURB%XCED * SQRT(PTKEM(:,:,:)) / PLEPS(:,:,:) * ( PEXPL*PTKEM(:,:,:) + PIMPL*ZRES(:,:,:) ) )
+DO JK=1,D%NKT
+ DO JJ=1,D%NJT
+ DO JI=1,D%NIT
+PRTKES(JI,JJ,JK) = PRTKES(JI,JJ,JK) + PRHODJ(JI,JJ,JK) * &
+ ( PDP(JI,JJ,JK) + PTP(JI,JJ,JK) &
+ - CSTURB%XCED * SQRT(PTKEM(JI,JJ,JK)) / PLEPS(JI,JJ,JK) * ( PEXPL*PTKEM(JI,JJ,JK) + PIMPL*ZRES(JI,JJ,JK) ) )
+!
+PTDIFF(JI,JJ,JK) = ZRES(JI,JJ,JK) / PTSTEP - PRTKES(JI,JJ,JK)/PRHODJ(JI,JJ,JK) &
+ & - PDP(JI,JJ,JK)- PTP(JI,JJ,JK) - PTDISS(JI,JJ,JK)
+ ENDDO
+ ENDDO
+ENDDO
!
-PTDIFF(:,:,:) = ZRES(:,:,:) / PTSTEP - PRTKES(:,:,:)/PRHODJ(:,:,:) &
- & - PDP(:,:,:)- PTP(:,:,:) - PTDISS(:,:,:)
-
IF (BUCONF%LBUDGET_TKE) CALL BUDGET_STORE_INIT( TBUDGETS(NBUDGET_TKE), 'TR', PRTKES(:, :, :) )
!
-PRTKES(:,:,:) = ZRES(:,:,:) * PRHODJ(:,:,:) / PTSTEP - PRTKEMS(:,:,:)
+DO JK=1,D%NKT
+ DO JJ=1,D%NJT
+ DO JI=1,D%NIT
+PRTKES(JI,JJ,JK) = ZRES(JI,JJ,JK) * PRHODJ(JI,JJ,JK) / PTSTEP - PRTKEMS(JI,JJ,JK)
+ ENDDO
+ ENDDO
+ENDDO
!
! stores the whole turbulent transport
!
@@ -363,16 +433,39 @@ IF (BUCONF%LBUDGET_TKE) CALL BUDGET_STORE_END( TBUDGETS(NBUDGET_TKE), 'TR', PRTK
!* 3. COMPUTE THE DISSIPATIVE HEATING
! -------------------------------
!
-PRTHLS(:,:,:) = PRTHLS(:,:,:) + CSTURB%XCED * SQRT(PTKEM(:,:,:)) / PLEPS(:,:,:) * &
- (PEXPL*PTKEM(:,:,:) + PIMPL*ZRES(:,:,:)) * PRHODJ(:,:,:) * PCOEF_DISS(:,:,:)
+DO JK=1,D%NKT
+ DO JJ=1,D%NJT
+ DO JI=1,D%NIT
+PRTHLS(JI,JJ,JK) = PRTHLS(JI,JJ,JK) + CSTURB%XCED * SQRT(PTKEM(JI,JJ,JK)) / PLEPS(JI,JJ,JK) * &
+ (PEXPL*PTKEM(JI,JJ,JK) + PIMPL*ZRES(JI,JJ,JK)) * PRHODJ(JI,JJ,JK) * PCOEF_DISS(JI,JJ,JK)
+ ENDDO
+ ENDDO
+ENDDO
!----------------------------------------------------------------------------
!
!* 4. STORES SOME DIAGNOSTICS
! -----------------------
!
-IF(PRESENT(PDISS)) PDISS(:,:,:) = -CSTURB%XCED * (PTKEM(:,:,:)**1.5) / PLEPS(:,:,:)
IF(PRESENT(PTR)) PTR=ZTR
-IF(PRESENT(PEDR)) PEDR = CSTURB%XCED * (PTKEM(:,:,:)**1.5) / PLEPS(:,:,:)
+IF(PRESENT(PDISS)) THEN
+ DO JK=1,D%NKT
+ DO JJ=1,D%NJT
+ DO JI=1,D%NIT
+ PDISS(JI,JJ,JK) = -CSTURB%XCED * (PTKEM(JI,JJ,JK)**1.5) / PLEPS(JI,JJ,JK)
+ ENDDO
+ ENDDO
+ENDDO
+END IF
+!
+IF(PRESENT(PEDR)) THEN
+ DO JK=1,D%NKT
+ DO JJ=1,D%NJT
+ DO JI=1,D%NIT
+ PEDR(JI,JJ,JK) = CSTURB%XCED * (PTKEM(JI,JJ,JK)**1.5) / PLEPS(JI,JJ,JK)
+ ENDDO
+ ENDDO
+ENDDO
+END IF
!
IF ( OTURB_DIAG .AND. TPFILE%LOPENED ) THEN
!
--
GitLab