From bd1dada2f81b2a8cd9777170e8ff6053ad39b51f Mon Sep 17 00:00:00 2001
From: Philippe WAUTELET <philippe.wautelet@aero.obs-mip.fr>
Date: Fri, 19 Feb 2021 11:56:02 +0100
Subject: [PATCH] Philippe 19/02/2021: budgets: TR source term correctly
 computed

---
 src/MNH/tke_eps_sources.f90 | 13 +++++++++++--
 1 file changed, 11 insertions(+), 2 deletions(-)

diff --git a/src/MNH/tke_eps_sources.f90 b/src/MNH/tke_eps_sources.f90
index f37dff102..4efe246be 100644
--- a/src/MNH/tke_eps_sources.f90
+++ b/src/MNH/tke_eps_sources.f90
@@ -273,7 +273,6 @@ IKE=KKU-JPVEXT_TURB*KKL
 ZKEFF(:,:,:) = PLM(:,:,:) * SQRT(PTKEM(:,:,:)) 
 
 if (lbudget_th)  call Budget_store_init( tbudgets(NBUDGET_TH),  'DISSH', prthls(:, :, :) )
-if (lbudget_tke) call Budget_store_init( tbudgets(NBUDGET_TKE), 'TR',    prtkes(:, :, :) )
 
 !----------------------------------------------------------------------------
 !
@@ -303,7 +302,7 @@ PDP(:,:,IKB) = PDP(:,:,IKB) * (1. + PDZZ(:,:,IKB+KKL)/PDZZ(:,:,IKB))
 ! Compute the source terms for TKE: ( ADVECtion + NUMerical DIFFusion + ..)
 ! + (Dynamical Production) + (Thermal Production) - (dissipation) 
 ZFLX(:,:,:) = XCED * SQRT(PTKEM(:,:,:)) / PLEPS(:,:,:)
-ZSOURCE(:,:,:) = PRTKES(:,:,:) / PRHODJ(:,:,:)  +  PRTKESM(:,:,:) / PRHODJ(:,:,:) &
+ZSOURCE(:,:,:) = ( PRTKES(:,:,:) +  PRTKESM(:,:,:) ) / PRHODJ(:,:,:) &
    - PTKEM(:,:,:) / PTSTEP &
    + PDP(:,:,:) + PTP(:,:,:) + PTR(:,:,:) - PEXPL * ZFLX(:,:,:) * PTKEM(:,:,:)
 !
@@ -379,6 +378,16 @@ end if
 !
 !*       2.5  computes the final RTKE and stores the whole turbulent transport
 !              with the removal of the advection part 
+
+if (lbudget_tke) then
+  !Store the previous source terms in prtkes before initializing the next one
+  PRTKES(:,:,:) = PRTKES(:,:,:) + PRHODJ(:,:,:) *                                                           &
+                  ( PDP(:,:,:) + PTP(:,:,:)                                                                 &
+                    - XCED * SQRT(PTKEM(:,:,:)) / PLEPS(:,:,:) * ( PEXPL*PTKEM(:,:,:) + PIMPL*ZRES(:,:,:) ) )
+
+  call Budget_store_init( tbudgets(NBUDGET_TKE), 'TR', prtkes(:, :, :) )
+end if
+
 PRTKES(:,:,:) = ZRES(:,:,:) * PRHODJ(:,:,:) / PTSTEP -  PRTKESM(:,:,:)
 !
 ! stores the whole turbulent transport
-- 
GitLab