diff --git a/src/MNH/budget.f90 b/src/MNH/budget.f90
index 5d3d18d1188165d942c51c1238a41e2caf903581..95f3a48d9cdfd6480e4768bfe217f7bf474250ba 100644
--- a/src/MNH/budget.f90
+++ b/src/MNH/budget.f90
@@ -5,13 +5,14 @@
!-----------------------------------------------------------------
! Modifications
! P. Wautelet 28/01/2020: new subroutines: Budget_store_init, Budget_store_end and Budget_source_id_find in new module mode_budget
+! P. Wautelet 17/08/2020: treat LES budgets correctly
!-----------------------------------------------------------------
!#################
module mode_budget
!#################
-use modd_budget, only: cbutype, nbutime, tbudgetdata
+use modd_budget, only: cbutype, lbu_enable, nbutime, tbudgetdata
use modi_cart_compress, only: Cart_compress
use modi_mask_compress, only: Mask_compress
@@ -30,6 +31,8 @@ public :: Budget_store_add
contains
subroutine Budget_store_init( tpbudget, hsource, pvars )
+ use modd_les, only: lles_call
+
type(tbudgetdata), intent(inout) :: tpbudget ! Budget datastructure
character(len=*), intent(in) :: hsource ! Name of the source term
real, dimension(:,:,:), intent(in) :: pvars ! Current value to be stored
@@ -38,6 +41,18 @@ subroutine Budget_store_init( tpbudget, hsource, pvars )
call Print_msg( NVERB_DEBUG, 'BUD', 'Budget_store_init', trim( tpbudget%cname )//':'//trim( hsource ) )
+ if ( lles_call ) then
+ if ( allocated( tpbudget%xtmplesstore ) ) then
+ call Print_msg( NVERB_ERROR, 'BUD', 'Budget_store_init', 'xtmplesstore already allocated' )
+ else
+ allocate( tpbudget%xtmplesstore( Size( pvars, 1 ), Size( pvars, 2 ), Size ( pvars, 3 ) ) )
+ end if
+ tpbudget%xtmplesstore(:, :, :) = pvars(:, :, :)
+ end if
+
+ ! Nothing else to do if budgets are not enabled
+ if ( .not. lbu_enable ) return
+
call Budget_source_id_find( tpbudget, hsource, iid )
if ( tpbudget%ntmpstoresource /= 0 ) then
@@ -76,15 +91,38 @@ end subroutine Budget_store_init
subroutine Budget_store_end( tpbudget, hsource, pvars )
+ use modd_les, only: lles_call
+
+ use modi_les_budget, only: Les_budget
+
type(tbudgetdata), intent(inout) :: tpbudget ! Budget datastructure
character(len=*), intent(in) :: hsource ! Name of the source term
real, dimension(:,:,:), intent(in) :: pvars ! Current value to be stored
integer :: iid ! Reference number of the current source term
integer :: igroup ! Number of the group where to store the source term
+ real, dimension(:,:,:), allocatable :: zvars_add
call Print_msg( NVERB_DEBUG, 'BUD', 'Budget_store_end', trim( tpbudget%cname )//':'//trim( hsource ) )
+ if ( lles_call ) then
+ if ( allocated( tpbudget%xtmplesstore ) ) then
+ ! Do the call to Les_budget with oadd=.true.
+ ! This is necessary when the call to Budget_store_init was done with pvars not strictly
+ ! equal to the source term
+ Allocate( zvars_add( Size( pvars, 1 ), Size( pvars, 2 ), Size ( pvars, 3 ) ) )
+ zvars_add(:, :, :) = pvars(:, :, :) - tpbudget%xtmplesstore(:, :, :)
+ call Les_budget( zvars_add, tpbudget%nid, hsource, oadd = .true. )
+ Deallocate( zvars_add )
+ Deallocate( tpbudget%xtmplesstore )
+ else
+ call Les_budget( pvars, tpbudget%nid, hsource, oadd = .false. )
+ end if
+ end if
+
+ ! Nothing to do if budgets are not enabled
+ if ( .not. lbu_enable ) return
+
call Budget_source_id_find( tpbudget, hsource, iid )
if ( tpbudget%tsources(iid)%lenabled ) then
@@ -152,6 +190,10 @@ end subroutine Budget_store_end
subroutine Budget_store_add( tpbudget, hsource, pvars )
+ use modd_les, only: lles_call
+
+ use modi_les_budget, only: Les_budget
+
type(tbudgetdata), intent(inout) :: tpbudget ! Budget datastructure
character(len=*), intent(in) :: hsource ! Name of the source term
real, dimension(:,:,:), intent(in) :: pvars ! Current value to be stored
@@ -161,6 +203,11 @@ subroutine Budget_store_add( tpbudget, hsource, pvars )
call Print_msg( NVERB_DEBUG, 'BUD', 'Budget_store_add', trim( tpbudget%cname )//':'//trim( hsource ) )
+ if ( lles_call ) call Les_budget( pvars, tpbudget%nid, hsource, oadd = .true. )
+
+ ! Nothing to do if budgets are not enabled
+ if ( .not. lbu_enable ) return
+
call Budget_source_id_find( tpbudget, hsource, iid )
if ( tpbudget%tsources(iid)%lenabled ) then
diff --git a/src/MNH/les_budget.f90 b/src/MNH/les_budget.f90
index da9b09c2022185ae5f4e27524d17ab32eb9f08ff..a5114416fe6f7bb49bf713e308d40a3fc894540b 100644
--- a/src/MNH/les_budget.f90
+++ b/src/MNH/les_budget.f90
@@ -9,12 +9,12 @@
!
INTERFACE
!
- SUBROUTINE LES_BUDGET(PVARS,KBUDN,HBUVAR)
+ SUBROUTINE LES_BUDGET(PVARS,KBUDN,HBUVAR,oadd)
REAL, DIMENSION(:,:,:), INTENT(IN) :: PVARS ! Source
-INTEGER , INTENT(IN) :: KBUDN ! variable number
-CHARACTER (LEN=*) , INTENT(IN) :: HBUVAR ! Identifier of the Budget of the
- ! variable that is considered
+INTEGER, INTENT(IN) :: KBUDN ! variable number
+CHARACTER(LEN=*), INTENT(IN) :: HBUVAR ! Identifier of the budget of the variable that is considered
+logical, intent(in) :: oadd ! Flag to determine if source is to be added (true) or if is absolute (false)
END SUBROUTINE LES_BUDGET
@@ -23,7 +23,7 @@ END INTERFACE
END MODULE MODI_LES_BUDGET
!
! ####################################
- SUBROUTINE LES_BUDGET(PVARS,KBUDN,HBUVAR)
+ SUBROUTINE LES_BUDGET(PVARS,KBUDN,HBUVAR,oadd)
! ####################################
!
!!**** *LES_BUDGET* - stores
@@ -48,6 +48,7 @@ END MODULE MODI_LES_BUDGET
!! Original September 19, 2002
!! 25/11/2016 Q.Rodier correction bug variance u'^2 v'^2 w'^2
! P. Wautelet 20/05/2019: add name argument to ADDnFIELD_ll + new ADD4DFIELD_ll subroutine
+! P. Wautelet 17/08/2020: treat LES budgets correctly
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -73,10 +74,10 @@ IMPLICIT NONE
! ------------------------------
!
REAL, DIMENSION(:,:,:), INTENT(IN) :: PVARS ! Source
-INTEGER , INTENT(IN) :: KBUDN ! variable number
-CHARACTER (LEN=*) , INTENT(IN) :: HBUVAR ! Identifier of the Budget of the
- ! variable that is considered
-
+INTEGER, INTENT(IN) :: KBUDN ! variable number
+CHARACTER(LEN=*), INTENT(IN) :: HBUVAR ! Identifier of the budget of the variable that is considered
+logical, intent(in) :: oadd ! Flag to determine if source is to be added (true) or if is absolute (false)
+!
!* 0.2 declaration of local variables
! ------------------------------
!
@@ -124,7 +125,11 @@ SELECT CASE (KBUDN)
!* u
!
CASE( NBUDGET_U )
- CALL LES_BUDGET_ANOMALY(PVARS,'X',ZANOM)
+ if ( oadd ) then
+ CALL LES_BUDGET_ANOMALY(XCURRENT_RUS + PVARS,'X',ZANOM)
+ else
+ CALL LES_BUDGET_ANOMALY(PVARS,'X',ZANOM)
+ end if
!
!* action in KE budget
ZWORK_LES = 0.5*( ZANOM ** 2 - XU_ANOM ** 2 ) / XCURRENT_TSTEP
@@ -132,13 +137,21 @@ SELECT CASE (KBUDN)
X_LES_BU_RES_KE(:,ILES_BU) = X_LES_BU_RES_KE(:,ILES_BU) + ZLES_PROF(:)
!
!* update fields
- XCURRENT_RUS = PVARS
+ if ( oadd ) then
+ XCURRENT_RUS = XCURRENT_RUS + PVARS
+ else
+ XCURRENT_RUS = PVARS
+ end if
XU_ANOM = ZANOM
!
!* v
!
CASE( NBUDGET_V )
- CALL LES_BUDGET_ANOMALY(PVARS,'Y',ZANOM)
+ if ( oadd ) then
+ CALL LES_BUDGET_ANOMALY(XCURRENT_RVS + PVARS,'Y',ZANOM)
+ else
+ CALL LES_BUDGET_ANOMALY(PVARS,'Y',ZANOM)
+ end if
!
!* action in KE budget
ZWORK_LES = 0.5*( ZANOM ** 2 - XV_ANOM ** 2 ) / XCURRENT_TSTEP
@@ -146,13 +159,21 @@ SELECT CASE (KBUDN)
X_LES_BU_RES_KE(:,ILES_BU) = X_LES_BU_RES_KE(:,ILES_BU) + ZLES_PROF(:)
!
!* update fields
- XCURRENT_RVS = PVARS
+ if ( oadd ) then
+ XCURRENT_RVS = XCURRENT_RVS + PVARS
+ else
+ XCURRENT_RVS = PVARS
+ end if
XV_ANOM = ZANOM
!
!* w
!
CASE( NBUDGET_W )
- CALL LES_BUDGET_ANOMALY(PVARS,'Z',ZANOM)
+ if ( oadd ) then
+ CALL LES_BUDGET_ANOMALY(XCURRENT_RWS + PVARS,'Z',ZANOM)
+ else
+ CALL LES_BUDGET_ANOMALY(PVARS,'Z',ZANOM)
+ end if
!
!* action in KE budget
ZWORK_LES = 0.5*( ZANOM ** 2 - XW_ANOM ** 2 ) / XCURRENT_TSTEP
@@ -179,13 +200,21 @@ SELECT CASE (KBUDN)
END DO
!
!* update fields
- XCURRENT_RWS = PVARS
+ if ( oadd ) then
+ XCURRENT_RWS = XCURRENT_RWS + PVARS
+ else
+ XCURRENT_RWS = PVARS
+ end if
XW_ANOM = ZANOM
!
!* Th
!
CASE( NBUDGET_TH )
- XCURRENT_RTHLS = XCURRENT_RTHLS + PVARS - XCURRENT_RTHS
+ if ( oadd ) then
+ XCURRENT_RTHLS = XCURRENT_RTHLS + PVARS
+ else
+ XCURRENT_RTHLS = XCURRENT_RTHLS + PVARS - XCURRENT_RTHS
+ end if
CALL LES_BUDGET_ANOMALY(XCURRENT_RTHLS,'-',ZANOM)
!
!* action in WTHL budget
@@ -207,15 +236,24 @@ SELECT CASE (KBUDN)
END IF
!
!* update fields
- XCURRENT_RTHS = PVARS
+ if ( oadd ) then
+ XCURRENT_RTHS = XCURRENT_RTHS + PVARS
+ else
+ XCURRENT_RTHS = PVARS
+ end if
XTHL_ANOM = ZANOM
!
!* Tke
!
CASE( NBUDGET_TKE )
ALLOCATE(ZTEND(IIU,IJU,IKU))
- ZTEND(:,:,:) = (PVARS(:,:,:)-XCURRENT_RTKES(:,:,:)) / XCURRENT_RHODJ
- XCURRENT_RTKES = PVARS
+ if ( oadd ) then
+ ZTEND(:,:,:) = PVARS(:,:,:) / XCURRENT_RHODJ
+ XCURRENT_RTKES = XCURRENT_RTKES+PVARS
+ else
+ ZTEND(:,:,:) = (PVARS(:,:,:)-XCURRENT_RTKES(:,:,:)) / XCURRENT_RHODJ
+ XCURRENT_RTKES = PVARS
+ end if
CALL LES_VER_INT( ZTEND, ZWORK_LES )
DEALLOCATE(ZTEND)
CALL LES_MEAN_ll( ZWORK_LES, LLES_CURRENT_CART_MASK, ZLES_PROF)
@@ -225,7 +263,11 @@ SELECT CASE (KBUDN)
!
CASE( NBUDGET_RV, NBUDGET_RR, NBUDGET_RI, NBUDGET_RS, NBUDGET_RG, NBUDGET_RH )
!* transformation into conservative variables: RT
- XCURRENT_RRTS = XCURRENT_RRTS + PVARS(:,:,:) - XCURRENT_RRS(:,:,:,KBUDN-(NBUDGET_RV-1))
+ if ( oadd ) then
+ XCURRENT_RRTS = XCURRENT_RRTS + PVARS(:,:,:)
+ else
+ XCURRENT_RRTS = XCURRENT_RRTS + PVARS(:,:,:) - XCURRENT_RRS(:,:,:,KBUDN-(NBUDGET_RV-1))
+ end if
CALL LES_BUDGET_ANOMALY(XCURRENT_RRTS,'-',ZANOM)
!
!* action in WRT budget
@@ -245,16 +287,29 @@ SELECT CASE (KBUDN)
X_LES_BU_RES_THLRT(:,ILES_BU) = X_LES_BU_RES_THLRT(:,ILES_BU) + ZLES_PROF(:)
!
!* update fields
- XCURRENT_RRS(:,:,:,KBUDN-(NBUDGET_RV-1)) = PVARS
+ if ( oadd ) then
+ XCURRENT_RRS(:,:,:,KBUDN-(NBUDGET_RV-1)) = XCURRENT_RRS(:,:,:,KBUDN-(NBUDGET_RV-1)) + PVARS
+ else
+ XCURRENT_RRS(:,:,:,KBUDN-(NBUDGET_RV-1)) = PVARS
+ end if
XRT_ANOM = ZANOM
!
!* Rc
!
CASE( NBUDGET_RC )
!* transformation into conservative variables: theta_l; RT
- XCURRENT_RRTS = XCURRENT_RRTS + PVARS(:,:,:) - XCURRENT_RRS(:,:,:,KBUDN-(NBUDGET_RV-1))
- XCURRENT_RTHLS = XCURRENT_RTHLS - XCURRENT_L_O_EXN_CP &
- * (PVARS(:,:,:) - XCURRENT_RRS(:,:,:,KBUDN-(NBUDGET_RV-1)))
+ if ( oadd ) then
+ XCURRENT_RRTS = XCURRENT_RRTS + PVARS(:,:,:)
+ else
+ XCURRENT_RRTS = XCURRENT_RRTS + PVARS(:,:,:) - XCURRENT_RRS(:,:,:,KBUDN-(NBUDGET_RV-1))
+ end if
+ if ( oadd ) then
+ XCURRENT_RTHLS = XCURRENT_RTHLS - XCURRENT_L_O_EXN_CP &
+ * (PVARS(:,:,:) )
+ else
+ XCURRENT_RTHLS = XCURRENT_RTHLS - XCURRENT_L_O_EXN_CP &
+ * (PVARS(:,:,:) - XCURRENT_RRS(:,:,:,KBUDN-(NBUDGET_RV-1)))
+ end if
!* anomaly of THL
ALLOCATE(ZTHL_ANOM(IIU,IJU,NLES_K))
@@ -292,7 +347,11 @@ SELECT CASE (KBUDN)
!
!
!* update fields
- XCURRENT_RRS(:,:,:,KBUDN-(NBUDGET_RV-1)) = PVARS
+ if ( oadd ) then
+ XCURRENT_RRS(:,:,:,KBUDN-(NBUDGET_RV-1)) = XCURRENT_RRS(:,:,:,KBUDN-(NBUDGET_RV-1)) + PVARS
+ else
+ XCURRENT_RRS(:,:,:,KBUDN-(NBUDGET_RV-1)) = PVARS
+ end if
XRT_ANOM = ZANOM
XTHL_ANOM = ZTHL_ANOM
DEALLOCATE(ZTHL_ANOM)
@@ -300,7 +359,11 @@ SELECT CASE (KBUDN)
!* SV
!
CASE( NBUDGET_SV1: )
- CALL LES_BUDGET_ANOMALY(PVARS,'-',ZANOM)
+ if ( oadd ) then
+ CALL LES_BUDGET_ANOMALY(XCURRENT_RSVS(:,:,:,KBUDN-(NBUDGET_SV1-1)) + PVARS,'-',ZANOM)
+ else
+ CALL LES_BUDGET_ANOMALY(PVARS,'-',ZANOM)
+ end if
!
!* action in WSV budget
ZWORK_LES = ( ZANOM * XW_ANOM - XSV_ANOM(:,:,:,KBUDN-(NBUDGET_SV1-1)) * XW_ANOM ) / &
@@ -315,7 +378,12 @@ SELECT CASE (KBUDN)
X_LES_BU_RES_SV2(:,ILES_BU,KBUDN-(NBUDGET_SV1-1)) = X_LES_BU_RES_SV2(:,ILES_BU,KBUDN-(NBUDGET_SV1-1)) + ZLES_PROF(:)
!
!* update fields
- XCURRENT_RSVS(:,:,:,KBUDN-(NBUDGET_SV1-1)) = PVARS
+ if ( oadd ) then
+ XCURRENT_RSVS(:,:,:,KBUDN-(NBUDGET_SV1-1)) = XCURRENT_RSVS(:,:,:,KBUDN-(NBUDGET_SV1-1)) + PVARS
+ else
+ XCURRENT_RSVS(:,:,:,KBUDN-(NBUDGET_SV1-1)) = PVARS
+ end if
+ XV_ANOM = ZANOM
XSV_ANOM(:,:,:,KBUDN-(NBUDGET_SV1-1)) = ZANOM
END SELECT
@@ -339,37 +407,37 @@ CHARACTER (LEN=*), INTENT(IN) :: HBU ! Identifier of the Budget of the
! variable that is considered
INTEGER, INTENT(OUT) :: KLES_BU ! LES budget identifier
!
-IF (HBU(1:3)=='ADV') THEN
+IF (HBU=='ADV') THEN
KLES_BU = NLES_TOTADV
-ELSE IF (HBU(1:3)=='REL') THEN
+ELSE IF (HBU=='REL') THEN
KLES_BU = NLES_RELA
-ELSE IF (HBU(1:5)=='VTURB') THEN
+ELSE IF (HBU=='VTURB') THEN
KLES_BU = NLES_VTURB
-ELSE IF (HBU(1:5)=='HTURB') THEN
+ELSE IF (HBU=='HTURB') THEN
KLES_BU = NLES_HTURB
-ELSE IF (HBU(1:4)=='GRAV') THEN
+ELSE IF (HBU=='GRAV') THEN
KLES_BU = NLES_GRAV
-ELSE IF (HBU(1:4)=='PRES') THEN
+ELSE IF (HBU=='PRES') THEN
KLES_BU = NLES_PRES
-ELSE IF (HBU(1:4)=='PREF') THEN
+ELSE IF (HBU=='PREF') THEN
KLES_BU = NLES_PREF
-ELSE IF (HBU(1:4)=='CURV') THEN
+ELSE IF (HBU=='CURV') THEN
KLES_BU = NLES_CURV
-ELSE IF (HBU(1:3)=='COR') THEN
+ELSE IF (HBU=='COR') THEN
KLES_BU = NLES_COR
-ELSE IF (HBU(1:2)=='DP') THEN
+ELSE IF (HBU=='DP') THEN
KLES_BU = NLES_DP
-ELSE IF (HBU(1:2)=='TP') THEN
+ELSE IF (HBU=='TP') THEN
KLES_BU = NLES_TP
-ELSE IF (HBU(1:2)=='TR') THEN
+ELSE IF (HBU=='TR') THEN
KLES_BU = NLES_TR
-ELSE IF (HBU(1:4)=='DISS') THEN
+ELSE IF (HBU=='DISS') THEN
KLES_BU = NLES_DISS
-ELSE IF (HBU(1:3)=='DIF') THEN
+ELSE IF (HBU=='DIF') THEN
KLES_BU = NLES_DIFF
-ELSE IF (HBU(1:3)=='RAD') THEN
+ELSE IF (HBU=='RAD') THEN
KLES_BU = NLES_RAD
-ELSE IF (HBU(1:4)=='NEST') THEN
+ELSE IF (HBU=='NEST') THEN
KLES_BU = NLES_NEST
ELSE
KLES_BU = NLES_MISC
diff --git a/src/MNH/modd_budget.f90 b/src/MNH/modd_budget.f90
index 486a6af78f3fc2b3ec7d6d6d7ac184961163bc45..abb8a950c679c498375b45b7352f82b177e15c2e 100644
--- a/src/MNH/modd_budget.f90
+++ b/src/MNH/modd_budget.f90
@@ -54,6 +54,7 @@
! P. Wautelet 17/04/2020: set default values for budgets switch values
! P. Wautelet 23/04/2020: add nid in tbudgetdata datatype
! P. Wautelet 30/06/2020: add NNETURSV, NNEADVSV and NNECONSV variables
+! P. Wautelet 17/08/2020: add xtmplesstore in tbudgetdata datatype
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -95,6 +96,8 @@ type tbudgetdata
logical :: lenabled = .false. ! True if corresponding budget flag is set to true
real, dimension(:,:,:), allocatable :: xtmpstore ! Array to store temporary data
! (to allow to store the difference between 2 places)
+ real, dimension(:,:,:), allocatable :: xtmplesstore ! Array to store temporary data for LES budgets
+ ! (to allow to store the difference between 2 places)
type(tbusourcedata), dimension(:), allocatable :: tsources ! Full list of source terms (used or not)
type(tbugroupdata), dimension(:), allocatable :: tgroups ! Full list of groups of source terms (to be written)
type(tburhodata), pointer :: trhodj => null() ! Budget array for rhodj