From 46cb33c1e2ddac0f07fba5cb92a36c7cab26c79e Mon Sep 17 00:00:00 2001
From: Quentin Rodier <quentin.rodier@meteo.fr>
Date: Tue, 7 May 2024 10:24:49 +0200
Subject: [PATCH] Quentin 07/05/2024: bugfix for use of CBL_HEIGHT_DEF=FRI:
bl_depth_diag can be called with arrays in the full vertical dimensions or
the les vertical dimensions. The 1D version was also wrong. The FRI option
was missing in the user guide (option forgotten since many years)
---
src/MNH/lesn.f90 | 2 +-
src/PHYEX/turb/mode_bl_depth_diag.f90 | 47 +++++++++++++++------------
src/PHYEX/turb/mode_sbl_depth.f90 | 4 +--
3 files changed, 30 insertions(+), 23 deletions(-)
diff --git a/src/MNH/lesn.f90 b/src/MNH/lesn.f90
index 86b86a3e6..7f1de616e 100644
--- a/src/MNH/lesn.f90
+++ b/src/MNH/lesn.f90
@@ -3277,7 +3277,7 @@ ELSE IF (CBL_HEIGHT_DEF=='FRI') THEN
+( XLES_SUBGRID_WV (:,NLES_CURRENT_TCOUNT,1) &
+XLES_RESOLVED_WV(:,NLES_CURRENT_TCOUNT,1))**2 )
ZFRIC_SURF = XLES_USTAR(NLES_CURRENT_TCOUNT)**2
- CALL BL_DEPTH_DIAG(YLDIMPHYEX,ZFRIC_SURF, XLES_ZS, &
+ CALL BL_DEPTH_DIAG(YLDIMPHYEX,SIZE(ZFRIC_LES),ZFRIC_SURF, XLES_ZS, &
ZFRIC_LES, XLES_Z, &
XFTOP_O_FSURF,XLES_BL_HEIGHT(NLES_CURRENT_TCOUNT))
END IF
diff --git a/src/PHYEX/turb/mode_bl_depth_diag.f90 b/src/PHYEX/turb/mode_bl_depth_diag.f90
index 2cae3a3fe..daafb106e 100644
--- a/src/PHYEX/turb/mode_bl_depth_diag.f90
+++ b/src/PHYEX/turb/mode_bl_depth_diag.f90
@@ -12,7 +12,7 @@ END INTERFACE
!
CONTAINS
!
-SUBROUTINE BL_DEPTH_DIAG_3D(D,PSURF,PZS,PFLUX,PZZ,PFTOP_O_FSURF,BL_DEPTH_DIAG3D)
+SUBROUTINE BL_DEPTH_DIAG_3D(D,KDIM, PSURF,PZS,PFLUX,PZZ,PFTOP_O_FSURF,BL_DEPTH_DIAG3D)
USE YOMHOOK , ONLY : LHOOK, DR_HOOK, JPHOOK
!
!
@@ -56,10 +56,11 @@ USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t
IMPLICIT NONE
!
TYPE(DIMPHYEX_t), INTENT(IN) :: D
+INTEGER, INTENT(IN) :: KDIM ! Vertical dimensions (LES grid vs full grid)
REAL, DIMENSION(D%NIJT), INTENT(IN) :: PSURF ! surface flux
REAL, DIMENSION(D%NIJT), INTENT(IN) :: PZS ! orography
-REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PFLUX ! flux
-REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PZZ ! altitude of flux points
+REAL, DIMENSION(D%NIJT,KDIM), INTENT(IN) :: PFLUX ! flux
+REAL, DIMENSION(D%NIJT,KDIM), INTENT(IN) :: PZZ ! altitude of flux points
REAL, INTENT(IN) :: PFTOP_O_FSURF! Flux at BL top / Surface flux
REAL, DIMENSION(D%NIJT), INTENT(OUT) :: BL_DEPTH_DIAG3D
!
@@ -106,7 +107,7 @@ BL_DEPTH_DIAG3D(:) = BL_DEPTH_DIAG3D(:) / (1. - PFTOP_O_FSURF)
IF (LHOOK) CALL DR_HOOK('BL_DEPTH_DIAG_3D',1,ZHOOK_HANDLE)
END SUBROUTINE BL_DEPTH_DIAG_3D
!
-SUBROUTINE BL_DEPTH_DIAG_1D(D,PSURF,PZS,PFLUX,PZZ,PFTOP_O_FSURF,BL_DEPTH_DIAG1D)
+SUBROUTINE BL_DEPTH_DIAG_1D(D,KDIM, PSURF,PZS,PFLUX,PZZ,PFTOP_O_FSURF,BL_DEPTH_DIAG1D)
USE YOMHOOK , ONLY : LHOOK, DR_HOOK, JPHOOK
!
USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t
@@ -114,31 +115,37 @@ USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t
IMPLICIT NONE
!
TYPE(DIMPHYEX_t), INTENT(IN) :: D
+INTEGER, INTENT(IN) :: KDIM ! Vertical dimensions (LES grid vs full grid)
REAL, INTENT(IN) :: PSURF ! surface flux
REAL, INTENT(IN) :: PZS ! orography
-REAL, DIMENSION(D%NKT), INTENT(IN) :: PFLUX ! flux
-REAL, DIMENSION(D%NKT), INTENT(IN) :: PZZ ! altitude of flux points
+REAL, DIMENSION(KDIM), INTENT(IN) :: PFLUX ! flux
+REAL, DIMENSION(KDIM), INTENT(IN) :: PZZ ! altitude of flux points
REAL, INTENT(IN) :: PFTOP_O_FSURF! Flux at BL top / Surface flux
REAL, INTENT(OUT) :: BL_DEPTH_DIAG1D
+INTEGER :: JK ! loop counters
+INTEGER :: IKB,IIJB,IIJE,IKL
+REAL :: ZFLX ! flux at top of BL
!
-REAL, DIMENSION(1,1) :: ZSURF
-REAL, DIMENSION(1,1) :: ZZS
-REAL, DIMENSION(1,1,D%NKT) :: ZFLUX
-REAL, DIMENSION(1,1,D%NKT) :: ZZZ
-REAL, DIMENSION(1,1) :: ZBL_DEPTH_DIAG
-!
-INTEGER :: IKT
REAL(KIND=JPHOOK) :: ZHOOK_HANDLE
IF (LHOOK) CALL DR_HOOK('BL_DEPTH_DIAG_1D',0,ZHOOK_HANDLE)
-IKT=D%NKT
-ZSURF = PSURF
-ZZS = PZS
-ZFLUX(1,1,:) = PFLUX(:)
-ZZZ (1,1,:) = PZZ (:)
!
-CALL BL_DEPTH_DIAG_3D(D,ZSURF,ZZS,ZFLUX,ZZZ,PFTOP_O_FSURF,ZBL_DEPTH_DIAG)
+IKB=D%NKTB
+IKL=D%NKL
!
-BL_DEPTH_DIAG1D = ZBL_DEPTH_DIAG(1,1)
+BL_DEPTH_DIAG1D = 0.
+IF (PSURF/=0.) THEN
+ DO JK=IKB,KDIM-1,IKL
+ IF (PZZ(JK-IKL)>PZS) THEN
+ ZFLX = PSURF * PFTOP_O_FSURF
+ IF ( (PFLUX(JK)-ZFLX)*(PFLUX(JK-IKL)-ZFLX) <= 0. ) THEN
+ BL_DEPTH_DIAG1D = (PZZ (JK-IKL) - PZS) &
+ + (PZZ (JK) - PZZ (JK-IKL)) &
+ * (ZFLX - PFLUX(JK-IKL) ) &
+ / (PFLUX(JK) - PFLUX(JK-IKL) )
+ END IF
+ END IF
+ END DO
+END IF
!
!-------------------------------------------------------------------------------
!
diff --git a/src/PHYEX/turb/mode_sbl_depth.f90 b/src/PHYEX/turb/mode_sbl_depth.f90
index 8da65ac56..948948c07 100644
--- a/src/PHYEX/turb/mode_sbl_depth.f90
+++ b/src/PHYEX/turb/mode_sbl_depth.f90
@@ -109,7 +109,7 @@ ZUSTAR2(:) = SQRT(ZWU(:)**2+ZWV(:)**2)
!$mnh_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
ZWIND(:,:)=SQRT(PFLXU(:,:)**2+PFLXV(:,:)**2)
!$mnh_end_expand_array(JIJ=IIJB:IIJE,JK=1:IKT)
-CALL BL_DEPTH_DIAG(D,ZUSTAR2,PZZ(:,IKB),ZWIND,PZZ,CSTURB%XFTOP_O_FSURF,ZSBL_DYN)
+CALL BL_DEPTH_DIAG(D,D%NKT,ZUSTAR2,PZZ(:,IKB),ZWIND,PZZ,CSTURB%XFTOP_O_FSURF,ZSBL_DYN)
!$mnh_expand_array(JIJ=IIJB:IIJE)
ZSBL_DYN(:) = CSTURB%XSBL_O_BL * ZSBL_DYN(:)
!$mnh_end_expand_array(JIJ=IIJB:IIJE)
@@ -118,7 +118,7 @@ ZSBL_DYN(:) = CSTURB%XSBL_O_BL * ZSBL_DYN(:)
!
!* BL and SBL diagnosed with buoyancy flux criteria
!
-CALL BL_DEPTH_DIAG(D,ZQ0,PZZ(:,IKB),PWTHV,PZZ,CSTURB%XFTOP_O_FSURF,ZSBL_THER)
+CALL BL_DEPTH_DIAG(D,D%NKT,ZQ0,PZZ(:,IKB),PWTHV,PZZ,CSTURB%XFTOP_O_FSURF,ZSBL_THER)
!$mnh_expand_array(JIJ=IIJB:IIJE)
ZSBL_THER(:)= CSTURB%XSBL_O_BL * ZSBL_THER(:)
!$mnh_end_expand_array(JIJ=IIJB:IIJE)
--
GitLab