diff --git a/src/MNH/default_desfmn.f90 b/src/MNH/default_desfmn.f90
index c656a19af214c8c32027fdd699cc8e3d6842b52a..b61a13e5e73f9637f38684089677fbdadf0038cc 100644
--- a/src/MNH/default_desfmn.f90
+++ b/src/MNH/default_desfmn.f90
@@ -219,6 +219,7 @@ END MODULE MODI_DEFAULT_DESFM_n
! Q. Rodier 06/2021: modify default value to LGZ=F (grey-zone corr.), LSEDI and OSEDC=T (LIMA sedimentation)
! F. Couvreux 06/2021: add LRELAX_UVMEAN_FRC
! Q. Rodier 07/2021: modify XPOND=1
+! R. Schoetter 12/2021 multi-level coupling between MesoNH and SURFEX
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -266,6 +267,7 @@ USE MODD_CONDSAMP
USE MODD_MEAN_FIELD
USE MODD_DRAGTREE_n
USE MODD_DRAGBLDG_n
+USE MODD_COUPLING_LEVELS_n
USE MODD_EOL_MAIN
USE MODD_EOL_ADNR
USE MODD_EOL_ALM
@@ -560,6 +562,19 @@ LDEPOTREE = .FALSE.
XVDEPOTREE = 0.02 ! 2 cm/s
!------------------------------------------------------------------------------
!
+!* 10b. SET DEFAULT VALUES FOR MODD_DRAGBLDG_n :
+! ----------------------------------
+!
+LDRAGBLDG = .FALSE.
+LFLUXBLDG = .FALSE.
+LDRAGURBVEG = .FALSE.
+!
+!* 10c. SET DEFAULT VALUES FOR MODD_COUPLING_LEVELS_n :
+! ----------------------------------
+!
+NLEV_COUPLE = 1
+!------------------------------------------------------------------------------
+!
!* 10c. SET DEFAULT VALUES FOR MODD_DRAGB
! ----------------------------------
!
diff --git a/src/MNH/drag_bld.f90 b/src/MNH/drag_bld.f90
index fbb25dc13bdcb4837a23d39130417db344f58a48..cb94923b9b3568736dd901f5a0c9a57b9bc4ac1e 100644
--- a/src/MNH/drag_bld.f90
+++ b/src/MNH/drag_bld.f90
@@ -9,17 +9,31 @@ MODULE MODI_DRAG_BLD
!
INTERFACE
!
- SUBROUTINE DRAG_BLD(PTSTEP, PUT, PVT, PTKET, PRHODJ, PZZ, PRUS, PRVS, PRTKES )
+ SUBROUTINE DRAG_BLD(PTSTEP, PUT, PVT, PTKET, PPABST, PTHT, &
+ PRT, PSVT, PRHODJ, PZZ, PRUS, PRVS, PRTKES, PRTHS, PRRS, &
+ PSFTH_WALL, PSFTH_ROOF, PCD_ROOF, PSFRV_WALL, PSFRV_ROOF )
!
REAL, INTENT(IN) :: PTSTEP ! Time step
REAL, DIMENSION(:,:,:), INTENT(IN) :: PUT, PVT ! variables
REAL, DIMENSION(:,:,:), INTENT(IN) :: PTKET ! at t
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PPABST ! at t
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHT ! at t
+ REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRT ! at t
+ REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PSVT ! at t
!
REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! dry Density * Jacobian
REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! Height (z)
!
+ REAL, DIMENSION(:,:), INTENT(IN) :: PSFTH_WALL ! Wall flux of theta
+ REAL, DIMENSION(:,:), INTENT(IN) :: PSFTH_ROOF ! Roof flux of theta
+ REAL, DIMENSION(:,:), INTENT(IN) :: PCD_ROOF ! Drag coefficient due to roofs
+ REAL, DIMENSION(:,:), INTENT(IN) :: PSFRV_WALL ! Wall flux of vapor
+ REAL, DIMENSION(:,:), INTENT(IN) :: PSFRV_ROOF ! Roof flux of vapor
+ !
REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRUS, PRVS ! Sources of Momentum
REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRTKES ! Sources of Tke
+ REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRRS
+ REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRTHS
!
END SUBROUTINE DRAG_BLD
@@ -28,7 +42,9 @@ MODULE MODI_DRAG_BLD
END MODULE MODI_DRAG_BLD
!
! ###################################################################
-SUBROUTINE DRAG_BLD(PTSTEP, PUT, PVT, PTKET, PRHODJ, PZZ, PRUS, PRVS, PRTKES )
+SUBROUTINE DRAG_BLD(PTSTEP, PUT, PVT, PTKET, PPABST, PTHT, PRT, &
+ PSVT, PRHODJ, PZZ, PRUS, PRVS, PRTKES, PRTHS, PRRS, &
+ PSFTH_WALL, PSFTH_ROOF, PCD_ROOF, PSFRV_WALL, PSFRV_ROOF )
! ###################################################################
!
!!**** *DRAG_BLD_n * -
@@ -52,13 +68,14 @@ SUBROUTINE DRAG_BLD(PTSTEP, PUT, PVT, PTKET, PRHODJ, PZZ, PRUS, PRVS, PRTKES )
!! -------------
!! Original 09/2019
! P. Wautelet 04/03/2021: budgets: add DRAGB source term
+ ! R. Schoetter 12/2021 multi-level coupling between MesoNH and SURFEX
!!---------------------------------------------------------------
!
!* 0. DECLARATIONS
! ------------
!
- use modd_budget, only: lbudget_u, lbudget_v, lbudget_tke, &
- NBUDGET_U, NBUDGET_V, NBUDGET_TKE, &
+ use modd_budget, only: lbudget_u, lbudget_v, lbudget_tke, lbudget_th, lbudget_rv, &
+ NBUDGET_U, NBUDGET_V, NBUDGET_TKE, NBUDGET_TH, NBUDGET_RV, &
tbudgets
USE MODD_CONF
USE MODD_CST
@@ -68,9 +85,9 @@ SUBROUTINE DRAG_BLD(PTSTEP, PUT, PVT, PTKET, PRHODJ, PZZ, PRUS, PRVS, PRTKES )
USE MODD_GROUND_PAR
USE MODD_NSV
USE MODD_PGDFIELDS
-
+
use mode_budget, only: Budget_store_init, Budget_store_end
-
+ !
USE MODI_MNHGET_SURF_PARAM_n
USE MODI_SHUMAN
!
@@ -81,18 +98,31 @@ SUBROUTINE DRAG_BLD(PTSTEP, PUT, PVT, PTKET, PRHODJ, PZZ, PRUS, PRVS, PRTKES )
REAL, INTENT(IN) :: PTSTEP ! Time step
REAL, DIMENSION(:,:,:), INTENT(IN) :: PUT, PVT ! variables
REAL, DIMENSION(:,:,:), INTENT(IN) :: PTKET ! at t
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PPABST ! at t
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHT ! at t
+ REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PRT ! at t
+ REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PSVT ! at t
!
REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! dry Density * Jacobian
REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! Height (z)
!
+ REAL, DIMENSION(:,:), INTENT(IN) :: PSFTH_WALL ! Wall flux of theta
+ REAL, DIMENSION(:,:), INTENT(IN) :: PSFTH_ROOF ! Roof flux of theta
+ REAL, DIMENSION(:,:), INTENT(IN) :: PCD_ROOF ! Drag coefficient due to roofs
+ REAL, DIMENSION(:,:), INTENT(IN) :: PSFRV_WALL ! Wall flux of vapor
+ REAL, DIMENSION(:,:), INTENT(IN) :: PSFRV_ROOF ! Roof flux of vapor
+ !
REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRUS, PRVS ! Sources of Momentum
REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRTKES ! Sources of Tke
+ REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRRS
+ REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRTHS
!
!* 0.2 Declarations of local variables :
!
INTEGER :: IIU,IJU,IKU,IKV ! array size along the k direction
INTEGER :: JI, JJ, JK ! loop index
INTEGER :: INFO_ll
+ INTEGER :: ICHECK
!
REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: &
ZWORK1, ZWORK2, ZWORK3, ZUT_SCAL, ZVT_SCAL, &
@@ -100,19 +130,37 @@ SUBROUTINE DRAG_BLD(PTSTEP, PUT, PVT, PTKET, PRHODJ, PZZ, PRUS, PRVS, PRTKES )
REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: &
ZCDRAG, ZDENSITY
!
- REAL, DIMENSION(:,:), ALLOCATABLE :: ZH_BUILD_PGD
- REAL, DIMENSION(:,:), ALLOCATABLE :: ZWALL_O_HOR_PGD
- REAL, DIMENSION(:,:), ALLOCATABLE :: ZFRAC_TOWN_PGD
- !
- REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2)) :: &
- ZH_BLD,ZF_BLD ! Building height, frontal density
- REAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2),SIZE(PZZ,3)):: ZT,ZEXN,ZLV,ZCPH
- !
- !* 0.3 Initialization
+ REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZCDRAG_BUILDG
+ REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZCDRAG_ROOF
+ REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZCDRAG_URBVEG
+ REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZDENSITY_BUILDG
+ REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZDENSITY_ROOF
+ REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZDENSITY_URBVEG
+ !
+ REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZDELTA_T_WALL
+ REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZDELTA_R_WALL
+ REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZDELTA_T_ROOF
+ REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZDELTA_R_ROOF
+ !
+ REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2)) :: ZFRAC_TOWN ! Town Fraction
+ REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2)) :: ZWALL_O_HOR ! Wall surface density
+ REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2)) :: ZH_BLD ! Building height
+ REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2)) :: ZF_BLD ! Wall frontal density
+ REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2)) :: ZLAI_URBVEG ! LAI of high urban vegetation
+ REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2)) :: ZH_URBTREE ! Height of trees in urban vegetation
+ REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2)) :: ZH_URBTRUN ! Height of trunks in urban vegetation
+ REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2)) :: ZFRAC_HVEG ! Fraction of high in urban vegetation
+ REAL, DIMENSION(SIZE(PUT,3)) :: ZLAD_CAN
+ REAL :: ZSUM_LAD_CAN, ZSUM_BLD_DENSITY, ZLEV_K0, ZLEV_K1
+ REAL :: ZSUM_SFTH_WALL, ZSUM_SFTH_ROOF, ZSUM_SFRV_WALL, ZSUM_SFRV_ROOF
+ !
+ !* 0.3 Initialisation de kelkes variables
!
if ( lbudget_u ) call Budget_store_init( tbudgets(NBUDGET_U ), 'DRAGB', prus (:, :, :) )
if ( lbudget_v ) call Budget_store_init( tbudgets(NBUDGET_V ), 'DRAGB', prvs (:, :, :) )
if ( lbudget_tke ) call Budget_store_init( tbudgets(NBUDGET_TKE), 'DRAGB', prtkes(:, :, :) )
+ if ( lbudget_th ) call Budget_store_init( tbudgets(NBUDGET_TH ), 'DRAGB', prths(:, :, :) )
+ if ( lbudget_rv ) call Budget_store_init( tbudgets(NBUDGET_RV ), 'DRAGB', prrs (:, :, :,1) )
IIU = SIZE(PUT,1)
IJU = SIZE(PUT,2)
@@ -122,37 +170,41 @@ SUBROUTINE DRAG_BLD(PTSTEP, PUT, PVT, PTKET, PRHODJ, PZZ, PRUS, PRVS, PRTKES )
ZVS (:,:,:) = 0.0
ZTKES (:,:,:) = 0.0
!
- ZH_BLD (:,:) = 0.
- ZF_BLD (:,:) = 0.
- !
- ZCDRAG (:,:,:) = 0.
- ZDENSITY (:,:,:) = 0.
- !
- ALLOCATE(ZFRAC_TOWN_PGD(IIU,IJU))
- ALLOCATE(ZH_BUILD_PGD(IIU,IJU))
- ALLOCATE(ZWALL_O_HOR_PGD(IIU,IJU))
- !
- ZFRAC_TOWN_PGD (:,:) = XUNDEF
- ZH_BUILD_PGD (:,:) = XUNDEF
- ZWALL_O_HOR_PGD (:,:) = XUNDEF
- !
- CALL MNHGET_SURF_PARAM_n( PTOWN=ZFRAC_TOWN_PGD, &
- PBUILD_HEIGHT=ZH_BUILD_PGD, &
- PWALL_O_HOR=ZWALL_O_HOR_PGD )
- !
- ! FIXME: Some values of ZFRAC_TOWN_PGD are 999. This is a bit strange since the
- ! TOWN fraction should be defined everywhere.
- ! It is set to 0.0 provisionally
- !
- WHERE(ZFRAC_TOWN_PGD(:,:).GT.1.0) ZFRAC_TOWN_PGD(:,:)=0.0
- !
- ! The values for wall density and building height are set to 0.0 where the Town fraction is 0.0
- ! For the wall density this would not be necessary since it will be multiplied by the town
- ! fraction anyway.
- !
- WHERE(ZFRAC_TOWN_PGD(:,:).EQ.0.0)
- ZWALL_O_HOR_PGD(:,:) = 0.0
- ZH_BUILD_PGD(:,:) = 0.0
+ ZFRAC_TOWN (:,:) = XUNDEF
+ ZWALL_O_HOR (:,:) = XUNDEF
+ ZH_BLD (:,:) = XUNDEF
+ ZF_BLD (:,:) = XUNDEF
+ ZLAI_URBVEG (:,:) = XUNDEF
+ ZH_URBTREE (:,:) = XUNDEF
+ ZH_URBTRUN (:,:) = XUNDEF
+ ZFRAC_HVEG (:,:) = XUNDEF
+ !
+ ZCDRAG_BUILDG (:,:,:) = 0.
+ ZCDRAG_ROOF (:,:,:) = 0.
+ ZCDRAG_URBVEG (:,:,:) = 0.
+ ZDENSITY_BUILDG (:,:,:) = 0.
+ ZDENSITY_ROOF (:,:,:) = 0.
+ ZDENSITY_URBVEG (:,:,:) = 0.
+ !
+ ZDELTA_T_WALL (:,:,:) = 0.
+ ZDELTA_R_WALL (:,:,:) = 0.
+ ZDELTA_T_ROOF (:,:,:) = 0.
+ ZDELTA_R_ROOF (:,:,:) = 0.
+ !
+ CALL MNHGET_SURF_PARAM_n( PTOWN = ZFRAC_TOWN, PBUILD_HEIGHT = ZH_BLD, &
+ PWALL_O_HOR = ZWALL_O_HOR, PLAI_HVEG = ZLAI_URBVEG, &
+ PH_URBTREE = ZH_URBTREE, PHTRUNK_HVEG = ZH_URBTRUN, &
+ PFRAC_HVEG=ZFRAC_HVEG )
+ !
+ WHERE(ZFRAC_TOWN(:,:).GT.1.0) ZFRAC_TOWN(:,:)=0.0
+ !
+ WHERE(ZFRAC_TOWN (:,:).EQ.0.0)
+ ZWALL_O_HOR (:,:) = 0.0
+ ZH_BLD (:,:) = 0.0
+ ZLAI_URBVEG (:,:) = 0.0
+ ZH_URBTREE (:,:) = 0.0
+ ZH_URBTRUN (:,:) = 0.0
+ ZFRAC_HVEG (:,:) = 0.0
ENDWHERE
!
! For buildings, the frontal wall area density is calculated [m^2(walls facing the wind)/m^2]
@@ -160,15 +212,32 @@ SUBROUTINE DRAG_BLD(PTSTEP, PUT, PVT, PTKET, PRHODJ, PZZ, PRUS, PRVS, PRTKES )
! [m^2(walls facing the wind)/m^2]=[m^2(wall)/m^2(town)*m^2(town)/m^2/PI]
! It will be assumed that the frontal wall area is equally distributed with height (all buildings same height)
!
- ZH_BLD(:,:) = ZH_BUILD_PGD(:,:)
- ZF_BLD(:,:) = ZFRAC_TOWN_PGD(:,:)*ZWALL_O_HOR_PGD(:,:)/XPI
+ ZF_BLD(:,:) = ZFRAC_TOWN (:,:) * ZWALL_O_HOR(:,:) / XPI
!
- DEALLOCATE(ZFRAC_TOWN_PGD)
- DEALLOCATE(ZH_BUILD_PGD)
- DEALLOCATE(ZWALL_O_HOR_PGD)
+ ! Set urban vegetation parameters to 0.0 in the case the drag shall not be accounted for
!
- !-------------------------------------------------------------------------------
+ IF (.NOT.LDRAGURBVEG) THEN
+ ZH_URBTREE (:,:) = 0.0
+ ZLAI_URBVEG (:,:) = 0.0
+ ZH_URBTRUN (:,:) = 0.0
+ ZFRAC_HVEG (:,:) = 0.0
+ ENDIF
!
+ ! Limit trunk height to 0.3 times the tree height and check afterwards
+ !
+ DO JJ=2,(IJU-1)
+ DO JI=2,(IIU-1)
+ !
+ ZH_URBTRUN(JI,JJ) = MIN(ZH_URBTRUN(JI,JJ),0.3*ZH_URBTREE(JI,JJ))
+ !
+ IF (ZH_URBTRUN(JI,JJ).GT.ZH_URBTREE(JI,JJ)) THEN
+ STOP ("Trunk higher than tree")
+ ENDIF
+ !
+ ENDDO
+ ENDDO
+ !
+ !-------------------------------------------------------------------------------
!
!* 1. COMPUTES THE TRUE VELOCITY COMPONENTS
! -------------------------------------
@@ -176,6 +245,7 @@ SUBROUTINE DRAG_BLD(PTSTEP, PUT, PVT, PTKET, PRHODJ, PZZ, PRUS, PRVS, PRTKES )
ZUT_SCAL(:,:,:) = MXF(PUT(:,:,:))
ZVT_SCAL(:,:,:) = MYF(PVT(:,:,:))
ZTKET(:,:,:) = PTKET(:,:,:)
+ !
!-------------------------------------------------------------------------------
!
!* 1. Computations of wind tendency due to canopy drag
@@ -189,24 +259,123 @@ SUBROUTINE DRAG_BLD(PTSTEP, PUT, PVT, PTKET, PRHODJ, PZZ, PRUS, PRVS, PRTKES )
DO JI=2,(IIU-1)
!
! Set density and drag coefficient for buildings
+ ! The drag coefficient is set to 0.4 although studies like Santiago and Martilli (2010)
+ ! provide hints that better results can be achieved with more sophisticated approaches.
!
IF (ZH_BLD(JI,JJ) /= 0) THEN
+ !
+ ZSUM_BLD_DENSITY = 0.0
+ ICHECK=0
!
DO JK=2,(IKU-1)
!
- IF ( (PZZ(JI,JJ,JK)-PZZ(JI,JJ,2)) .LT. ZH_BLD(JI,JJ) ) THEN
+ ! Height above ground of w-model levels (grid boundaries for density calculation)
+ !
+ ZLEV_K0 = PZZ(JI,JJ,JK ) - PZZ(JI,JJ,2)
+ ZLEV_K1 = PZZ(JI,JJ,JK+1) - PZZ(JI,JJ,2)
+ !
+ IF ( (ZLEV_K0.LT.ZH_BLD(JI,JJ)).AND.(ZLEV_K1.LT.ZH_BLD(JI,JJ)) ) THEN
+ !
+ ZCDRAG_BUILDG(JI,JJ,JK) = 0.4
+ !
+ ZDENSITY_BUILDG(JI,JJ,JK) = ZF_BLD(JI,JJ) / ZH_BLD(JI,JJ)
+ !
+ ELSE IF ( (ZLEV_K0.LT.ZH_BLD(JI,JJ)).AND.(ZLEV_K1.GE.ZH_BLD(JI,JJ)) ) THEN
!
- ! FIXME: Check literature values for the drag coefficient here
+ ICHECK=ICHECK+1
!
- ZCDRAG(JI,JJ,JK) = 0.2
+ ZCDRAG_BUILDG(JI,JJ,JK) = 0.4
!
- ! A uniform distribution of building heights is assumed
+ ZDENSITY_BUILDG(JI,JJ,JK) = ((ZH_BLD(JI,JJ)-ZLEV_K0)/(ZLEV_K1-ZLEV_K0)) * &
+ (ZF_BLD(JI,JJ) / ZH_BLD(JI,JJ))
!
- ZDENSITY(JI,JJ,JK) = ZF_BLD(JI,JJ) / ZH_BLD(JI,JJ)
+ ZCDRAG_ROOF (JI,JJ,JK) = PCD_ROOF(JI,JJ)
+ ZDENSITY_ROOF(JI,JJ,JK) = 1.0 / ( ZLEV_K1 - ZLEV_K0 )
!
ENDIF
!
+ ZSUM_BLD_DENSITY = ZSUM_BLD_DENSITY + ZDENSITY_BUILDG(JI,JJ,JK) * (ZLEV_K1-ZLEV_K0)
+ !
ENDDO
+ !
+ IF ( ICHECK .NE. 1 ) THEN
+ STOP ("Roof level not attributed")
+ ENDIF
+ !
+ IF ( ABS(ZSUM_BLD_DENSITY-ZF_BLD(JI,JJ)) .GT. 1.0E-6 ) THEN
+ STOP ("Wrong normalisation of frontal area density")
+ ENDIF
+ !
+ ENDIF
+ !
+ ! Set density and drag coefficient for urban vegetation
+ !
+ IF (ZH_URBTREE(JI,JJ) /= 0) THEN
+ !
+ ZLAD_CAN(:) = 0.0
+ ZSUM_LAD_CAN = 0.0
+ !
+ DO JK=2,(IKU-1)
+ !
+ ! Height above ground of model levels
+ !
+ ZLEV_K0 = PZZ(JI,JJ,JK ) - PZZ(JI,JJ,2)
+ ZLEV_K1 = PZZ(JI,JJ,JK+1) - PZZ(JI,JJ,2)
+ !
+ ! A) Model level entirely in between trunk and tree height
+ !
+ IF ( (ZLEV_K0.LT.ZH_URBTREE(JI,JJ)).AND.(ZLEV_K1.LT.ZH_URBTREE(JI,JJ)) .AND. &
+ (ZLEV_K0.GT.ZH_URBTRUN(JI,JJ))) THEN
+ !
+ ZLAD_CAN(JK) = ZLAI_URBVEG(JI,JJ) / ( ZH_URBTREE(JI,JJ)-ZH_URBTRUN(JI,JJ) )
+ !
+ !
+ ! B) Model level intersects tree height, but not trunk height
+ !
+ ELSE IF ( (ZLEV_K0.LT.ZH_URBTREE(JI,JJ)).AND.(ZLEV_K1.GT.ZH_URBTREE(JI,JJ)).AND. &
+ (ZLEV_K0.GT.ZH_URBTRUN(JI,JJ)) ) THEN
+ !
+ ZLAD_CAN(JK) = (ZH_URBTREE(JI,JJ)-ZLEV_K0)/(ZLEV_K1-ZLEV_K0) * &
+ ( ZLAI_URBVEG(JI,JJ) / (ZH_URBTREE(JI,JJ)-ZH_URBTRUN(JI,JJ) ) )
+ !
+ !
+ ! C) Model level intersects trunk height, but not tree height
+ !
+ ELSE IF ( (ZLEV_K1.LT.ZH_URBTREE(JI,JJ)).AND.(ZLEV_K0.LT.ZH_URBTRUN(JI,JJ)).AND. &
+ (ZLEV_K1.GT.ZH_URBTRUN(JI,JJ)) ) THEN
+ !
+ ZLAD_CAN(JK) = (ZLEV_K1-ZH_URBTRUN(JI,JJ))/(ZLEV_K1-ZLEV_K0) * &
+ ( ZLAI_URBVEG(JI,JJ) / (ZH_URBTREE(JI,JJ)-ZH_URBTRUN(JI,JJ) ) )
+ !
+ !
+ ! D) Model level intersects both tree and trunk height
+ !
+ ELSE IF ( (ZLEV_K0.LT.ZH_URBTREE(JI,JJ)).AND.(ZLEV_K1.GT.ZH_URBTREE(JI,JJ)).AND. &
+ (ZLEV_K0.LT.ZH_URBTRUN(JI,JJ)) ) THEN
+ !
+ ZLAD_CAN(JK) = (ZH_URBTREE(JI,JJ)-ZH_URBTRUN(JI,JJ))/(ZLEV_K1-ZLEV_K0) * &
+ ( ZLAI_URBVEG(JI,JJ) / (ZH_URBTREE(JI,JJ)-ZH_URBTRUN(JI,JJ) ) )
+ !
+ ! Remark: This equation simplifies to ZLAD_CAN(JK) = ZLAI_URBVEG(JI,JJ) /(ZLEV_K1-ZLEV_K0)
+ !
+ ENDIF
+ !
+ ZSUM_LAD_CAN = ZSUM_LAD_CAN + ZLAD_CAN(JK)*(ZLEV_K1-ZLEV_K0)
+ !
+ ! The division by PI assumes isotropic orientation of leaves.
+ ! To me it is not fully clear wether this is also considered in drag_veg.
+ !
+ ZDENSITY_URBVEG(JI,JJ,JK) = ZFRAC_TOWN(JI,JJ) * ZFRAC_HVEG(JI,JJ) * ZLAD_CAN(JK) / XPI
+ ZCDRAG_URBVEG(JI,JJ,JK) = 0.20
+ !
+ ENDDO
+ !
+ ! Check for correct normalisation of PLAD_CAN
+ !
+ IF ( ABS(ZSUM_LAD_CAN-ZLAI_URBVEG(JI,JJ)).GT.1.0E-6 ) THEN
+ STOP ("Wrong normalisation of vegetation density")
+ ENDIF
+ !
ENDIF
!
ENDDO
@@ -217,35 +386,115 @@ SUBROUTINE DRAG_BLD(PTSTEP, PUT, PVT, PTKET, PRHODJ, PZZ, PRUS, PRVS, PRTKES )
!
!* drag force by vertical surfaces
!
- ZUS(:,:,:) = PUT(:,:,:)/( 1.0 + MXM ( ZCDRAG(:,:,:) * ZDENSITY(:,:,:) &
- * PTSTEP * SQRT(ZUT_SCAL(:,:,:)**2+ZVT_SCAL(:,:,:)**2) ) )
+ ZUS(:,:,:) = PUT(:,:,:)/( 1.0 + MXM ( ( ZCDRAG_BUILDG(:,:,:) * ZDENSITY_BUILDG(:,:,:) + &
+ ZCDRAG_ROOF (:,:,:) * ZDENSITY_ROOF (:,:,:) + &
+ ZCDRAG_URBVEG(:,:,:) * ZDENSITY_URBVEG(:,:,:) ) &
+ * PTSTEP * SQRT(ZUT_SCAL(:,:,:)**2+ZVT_SCAL(:,:,:)**2) ) )
!
- ZVS(:,:,:) = PVT(:,:,:)/( 1.0 + MYM ( ZCDRAG(:,:,:) * ZDENSITY(:,:,:) &
- * PTSTEP * SQRT(ZUT_SCAL(:,:,:)**2+ZVT_SCAL(:,:,:)**2) ) )
+ ZVS(:,:,:) = PVT(:,:,:)/( 1.0 + MYM ( ( ZCDRAG_BUILDG(:,:,:) * ZDENSITY_BUILDG(:,:,:) + &
+ ZCDRAG_ROOF (:,:,:) * ZDENSITY_ROOF (:,:,:) + &
+ ZCDRAG_URBVEG(:,:,:) * ZDENSITY_URBVEG(:,:,:) ) &
+ * PTSTEP * SQRT(ZUT_SCAL(:,:,:)**2+ZVT_SCAL(:,:,:)**2) ) )
!
PRUS(:,:,:) = PRUS(:,:,:) + (ZUS(:,:,:)-PUT(:,:,:)) * MXM(PRHODJ(:,:,:)) / PTSTEP
!
PRVS(:,:,:) = PRVS(:,:,:) + (ZVS(:,:,:)-PVT(:,:,:)) * MYM(PRHODJ(:,:,:)) / PTSTEP
!
+ !
!* 3. Computations of TKE tendency due to canopy drag
! ------------------------------------------------
!* 3.1 Creation of TKE by wake
! -----------------------
!
- ! from Kanda and Hino (1994)
+ ZTKES(:,:,:) = ( ZTKET(:,:,:) + PTSTEP * ( ZCDRAG_URBVEG(:,:,:) * ZDENSITY_URBVEG(:,:,:) + &
+ ZCDRAG_ROOF (:,:,:) * ZDENSITY_ROOF (:,:,:) + &
+ ZCDRAG_BUILDG(:,:,:) * ZDENSITY_BUILDG(:,:,:) ) &
+ * (SQRT( ZUT_SCAL(:,:,:)**2 + ZVT_SCAL(:,:,:)**2 ))**3 ) / &
+ ( 1. + PTSTEP * ZCDRAG_URBVEG(:,:,:) * ZDENSITY_URBVEG(:,:,:) * &
+ SQRT(ZUT_SCAL(:,:,:)**2+ZVT_SCAL(:,:,:)**2) )
!
- ! Ext = + Cd * u+^3 * Sv/Vair vertical surfaces or trees
+ PRTKES(:,:,:) = PRTKES(:,:,:) + ( ZTKES(:,:,:) - ZTKET(:,:,:) ) * PRHODJ(:,:,:) / PTSTEP
!
! with Vair = Vair/Vtot * Vtot = (Vair/Vtot) * Stot * Dz
! and Sv/Vair = (Sv/Stot) * Stot/Vair = (Sv/Stot) / (Vair/Vtot) / Dz
!
- ZTKES(:,:,:) = ZTKET(:,:,:) + &
- PTSTEP * ZCDRAG(:,:,:) * ZDENSITY(:,:,:) * (SQRT( ZUT_SCAL(:,:,:)**2 + ZVT_SCAL(:,:,:)**2 ))**3
+ !* 4. Computations of temperature and mixing ratio tendency due to wall and roof heat fluxes
+ ! ------------------------------------------------
+ !
+ DO JJ=2,(IJU-1)
+ DO JI=2,(IIU-1)
+ !
+ IF (ZH_BLD(JI,JJ) /= 0) THEN
+ !
+ ZSUM_SFTH_WALL = 0.0
+ ZSUM_SFTH_ROOF = 0.0
+ ZSUM_SFRV_WALL = 0.0
+ ZSUM_SFRV_ROOF = 0.0
+ !
+ DO JK=2,(IKU-1)
+ !
+ ! Height above ground of model levels
+ !
+ ZLEV_K0 = PZZ(JI,JJ,JK ) - PZZ(JI,JJ,2)
+ ZLEV_K1 = PZZ(JI,JJ,JK+1) - PZZ(JI,JJ,2)
+ !
+ ! Fluxes from walls distributed equally over entire building height
+ !
+ IF ( (ZLEV_K0.LT.ZH_BLD(JI,JJ)).AND.(ZLEV_K1.LT.ZH_BLD(JI,JJ)) ) THEN
+ ZDELTA_T_WALL(JI,JJ,JK) = PSFTH_WALL(JI,JJ)*PTSTEP/ZH_BLD(JI,JJ)
+ ZDELTA_R_WALL(JI,JJ,JK) = PSFRV_WALL(JI,JJ)*PTSTEP/ZH_BLD(JI,JJ)
+ ELSE IF ( (ZLEV_K0.LT.ZH_BLD(JI,JJ)).AND.(ZLEV_K1.GE.ZH_BLD(JI,JJ)) ) THEN
+ ZDELTA_T_WALL(JI,JJ,JK) = ((ZH_BLD(JI,JJ)-ZLEV_K0)/(ZLEV_K1-ZLEV_K0)) * &
+ (PSFTH_WALL(JI,JJ)*PTSTEP/ZH_BLD(JI,JJ))
+ ZDELTA_R_WALL(JI,JJ,JK) = ((ZH_BLD(JI,JJ)-ZLEV_K0)/(ZLEV_K1-ZLEV_K0)) * &
+ (PSFRV_WALL(JI,JJ)*PTSTEP/ZH_BLD(JI,JJ))
+ ENDIF
+ !
+ ZSUM_SFTH_WALL = ZSUM_SFTH_WALL + (ZLEV_K1-ZLEV_K0) * ZDELTA_T_WALL(JI,JJ,JK) / PTSTEP
+ ZSUM_SFRV_WALL = ZSUM_SFRV_WALL + (ZLEV_K1-ZLEV_K0) * ZDELTA_R_WALL(JI,JJ,JK) / PTSTEP
+ !
+ ! Fluxes from roofs are applied at roof level
+ !
+ IF ( (ZLEV_K0.LT.ZH_BLD(JI,JJ)).AND.(ZLEV_K1.GE.ZH_BLD(JI,JJ)) ) THEN
+ ZDELTA_T_ROOF(JI,JJ,JK) = PSFTH_ROOF(JI,JJ)*PTSTEP/(ZLEV_K1-ZLEV_K0)
+ ZDELTA_R_ROOF(JI,JJ,JK) = PSFRV_ROOF(JI,JJ)*PTSTEP/(ZLEV_K1-ZLEV_K0)
+ ENDIF
+ !
+ ZSUM_SFTH_ROOF = ZSUM_SFTH_ROOF + (ZLEV_K1-ZLEV_K0) * ZDELTA_T_ROOF(JI,JJ,JK) / PTSTEP
+ ZSUM_SFRV_ROOF = ZSUM_SFRV_ROOF + (ZLEV_K1-ZLEV_K0) * ZDELTA_R_ROOF(JI,JJ,JK) / PTSTEP
+ !
+ ENDDO
+ !
+ ! Check for correct normalisation of fluxes
+ !
+ IF ( ABS(ZSUM_SFTH_WALL-PSFTH_WALL(JI,JJ)).GT.1.0E-6 ) THEN
+ STOP ("Wrong normalisation of wall heat flux")
+ ENDIF
+ !
+ IF ( ABS(ZSUM_SFRV_WALL-PSFRV_WALL(JI,JJ)).GT.1.0E-6 ) THEN
+ STOP ("Wrong normalisation of wall evaporative flux")
+ ENDIF
+ !
+ IF ( ABS(ZSUM_SFTH_ROOF-PSFTH_ROOF(JI,JJ)).GT.1.0E-6 ) THEN
+ STOP ("Wrong normalisation of roof heat flux")
+ ENDIF
+ !
+ IF ( ABS(ZSUM_SFRV_ROOF-PSFRV_ROOF(JI,JJ)).GT.1.0E-6 ) THEN
+ STOP ("Wrong normalisation of roof evaporative flux")
+ ENDIF
+ !
+ ENDIF
+ !
+ ENDDO
+ ENDDO
+ !
+ PRTHS(:,:,:) = PRTHS(:,:,:) + ( ZDELTA_T_WALL(:,:,:) + ZDELTA_T_ROOF(:,:,:) ) * PRHODJ(:,:,:) / PTSTEP
+ PRRS(:,:,:,1) = PRRS(:,:,:,1) + ( ZDELTA_R_WALL(:,:,:) + ZDELTA_R_ROOF(:,:,:) ) * PRHODJ(:,:,:) / PTSTEP
!
- PRTKES(:,:,:) = PRTKES(:,:,:) + (ZTKES(:,:,:)-ZTKET(:,:,:))*PRHODJ(:,:,:)/PTSTEP
-
if ( lbudget_u ) call Budget_store_end( tbudgets(NBUDGET_U ), 'DRAGB', prus (:, :, :) )
if ( lbudget_v ) call Budget_store_end( tbudgets(NBUDGET_V ), 'DRAGB', prvs (:, :, :) )
if ( lbudget_tke ) call Budget_store_end( tbudgets(NBUDGET_TKE), 'DRAGB', prtkes(:, :, :) )
-
+ if ( lbudget_th ) call Budget_store_end( tbudgets(NBUDGET_TH ), 'DRAGB', prths (:, :, :) )
+ if ( lbudget_rv ) call Budget_store_end( tbudgets(NBUDGET_RV ), 'DRAGB', prrs (:, :, :,1) )
+ !
END SUBROUTINE DRAG_BLD
diff --git a/src/MNH/goto_model_wrapper.f90 b/src/MNH/goto_model_wrapper.f90
index 831cb2028c3c0ea54da0090bb1329b014437bb4c..11b0d2c9d02c5d4955d277d1d1e4b7f0268cdd0e 100644
--- a/src/MNH/goto_model_wrapper.f90
+++ b/src/MNH/goto_model_wrapper.f90
@@ -18,6 +18,8 @@
! 11/2019 C.Lac correction in the drag formula and application to building in addition to tree
! F. Auguste 02/21: add IBM
! T. Nagel 02/21: add turbulence recycling
+
+ ! R. Schoetter 12/2021 multi-level coupling between MesoNH and SURFEX
!-----------------------------------------------------------------
MODULE MODI_GOTO_MODEL_WRAPPER
@@ -57,6 +59,7 @@ USE MODD_DIM_n
USE MODD_DRAG_n
USE MODD_DRAGTREE_n
USE MODD_DRAGBLDG_n
+USE MODD_COUPLING_LEVELS_n
USE MODD_DUMMY_GR_FIELD_n
USE MODD_DYN_n
USE MODD_DYNZD_n
@@ -152,6 +155,7 @@ CALL CURVCOR_GOTO_MODEL(KFROM, KTO)
CALL DIM_GOTO_MODEL(KFROM, KTO)
CALL DRAGTREE_GOTO_MODEL(KFROM, KTO)
CALL DRAGBLDG_GOTO_MODEL(KFROM, KTO)
+CALL COUPLING_MULT_GOTO_MODEL(KFROM, KTO)
CALL DUMMY_GR_FIELD_GOTO_MODEL(KFROM, KTO)
CALL DYN_GOTO_MODEL(KFROM, KTO)
CALL DYNZD_GOTO_MODEL(KFROM,KTO)
diff --git a/src/MNH/ground_paramn.f90 b/src/MNH/ground_paramn.f90
index 5ea7136cad6d2cbef0e8f39a2609c3cd0b5f5331..46aae3199a509e935245b2ac28f19a81fb777767 100644
--- a/src/MNH/ground_paramn.f90
+++ b/src/MNH/ground_paramn.f90
@@ -9,14 +9,20 @@ MODULE MODI_GROUND_PARAM_n
!
INTERFACE
!
- SUBROUTINE GROUND_PARAM_n( PSFTH, PSFRV, PSFSV, PSFCO2, PSFU, PSFV, &
- PDIR_ALB, PSCA_ALB, PEMIS, PTSRAD )
+ SUBROUTINE GROUND_PARAM_n( PSFTH, PSFTH_WALL, PSFTH_ROOF, PCD_ROOF, PSFRV, PSFRV_WALL, &
+ PSFRV_ROOF, PSFSV, PSFCO2, PSFU, PSFV, PDIR_ALB, PSCA_ALB, &
+ PEMIS, PTSRAD )
!
!* surface fluxes
! --------------
!
-REAL, DIMENSION(:,:), INTENT(OUT) :: PSFTH ! surface flux of potential temperature (Km/s)
-REAL, DIMENSION(:,:), INTENT(OUT) :: PSFRV ! surface flux of water vapor (m/s*kg/kg)
+REAL, DIMENSION(:,:), INTENT(OUT) :: PSFTH ! Total surface flux of potential temperature (Km/s)
+REAL, DIMENSION(:,:), INTENT(OUT) :: PSFTH_WALL ! Wall surface flux of potential temperature (Km/s)
+REAL, DIMENSION(:,:), INTENT(OUT) :: PSFTH_ROOF ! Roof surface flux of potential temperature (Km/s)
+REAL, DIMENSION(:,:), INTENT(OUT) :: PCD_ROOF ! Drag coefficient for roofs (-)
+REAL, DIMENSION(:,:), INTENT(OUT) :: PSFRV ! Total surface flux of water vapor (m/s*kg/kg)
+REAL, DIMENSION(:,:), INTENT(OUT) :: PSFRV_WALL ! Wall surface flux of water vapor (m/s*kg/kg)
+REAL, DIMENSION(:,:), INTENT(OUT) :: PSFRV_ROOF ! Roof surface flux of water vapor (m/s*kg/kg)
REAL, DIMENSION(:,:,:),INTENT(OUT):: PSFSV ! surface flux of scalar (m/s*kg/kg)
! flux of chemical var. (ppp.m/s)
REAL, DIMENSION(:,:), INTENT(OUT) :: PSFCO2! surface flux of CO2 (m/s*kg/kg)
@@ -38,8 +44,9 @@ END INTERFACE
END MODULE MODI_GROUND_PARAM_n
!
! ######################################################################
- SUBROUTINE GROUND_PARAM_n( PSFTH, PSFRV, PSFSV, PSFCO2, PSFU, PSFV, &
- PDIR_ALB, PSCA_ALB, PEMIS, PTSRAD )
+ SUBROUTINE GROUND_PARAM_n( PSFTH, PSFTH_WALL, PSFTH_ROOF, PCD_ROOF, PSFRV, &
+ PSFRV_WALL, PSFRV_ROOF, PSFSV, PSFCO2, PSFU, &
+ PSFV, PDIR_ALB, PSCA_ALB, PEMIS, PTSRAD )
! #######################################################################
!
!
@@ -111,6 +118,7 @@ END MODULE MODI_GROUND_PARAM_n
!! (V. Vionnet) 18/07/2017 add coupling for blowing snow module
!! (Bielli S.) 02/2019 Sea salt : significant sea wave height influences salt emission; 5 salt modes
! P. Wautelet 20/05/2019: add name argument to ADDnFIELD_ll + new ADD4DFIELD_ll subroutine
+! R. Schoetter 12/2021 multi-level coupling between MesoNH and SURFEX
! P. Wautelet 09/02/2022: bugfix: add missing XCURRENT_LEI computation
!-------------------------------------------------------------------------------
!
@@ -162,6 +170,8 @@ USE MODD_CSTS_DUST, ONLY : XMOLARWEIGHT_DUST
USE MODD_CSTS_SALT, ONLY : XMOLARWEIGHT_SALT
USE MODD_CH_FLX_n, ONLY : XCHFLX
USE MODD_DIAG_FLAG, ONLY : LCHEMDIAG
+USE MODD_DRAGBLDG_n, ONLY : LFLUXBLDG
+USE MODD_COUPLING_LEVELS_n
!
USE MODI_NORMAL_INTERPOL
USE MODI_ROTATE_WIND
@@ -195,8 +205,13 @@ IMPLICIT NONE
!* surface fluxes
! --------------
!
-REAL, DIMENSION(:,:), INTENT(OUT) :: PSFTH ! surface flux of potential temperature (Km/s)
-REAL, DIMENSION(:,:), INTENT(OUT) :: PSFRV ! surface flux of water vapor (m/s*kg/kg)
+REAL, DIMENSION(:,:), INTENT(OUT) :: PSFTH ! Total surface flux of potential temperature (Km/s)
+REAL, DIMENSION(:,:), INTENT(OUT) :: PSFTH_WALL ! Wall surface flux of potential temperature (Km/s)
+REAL, DIMENSION(:,:), INTENT(OUT) :: PSFTH_ROOF ! Roof surface flux of potential temperature (Km/s)
+REAL, DIMENSION(:,:), INTENT(OUT) :: PCD_ROOF ! Drag coefficient for roofs (-)
+REAL, DIMENSION(:,:), INTENT(OUT) :: PSFRV ! Total surface flux of water vapor (m/s*kg/kg)
+REAL, DIMENSION(:,:), INTENT(OUT) :: PSFRV_WALL ! Wall surface flux of water vapor (m/s*kg/kg)
+REAL, DIMENSION(:,:), INTENT(OUT) :: PSFRV_ROOF ! Roof surface flux of water vapor (m/s*kg/kg)
REAL, DIMENSION(:,:,:),INTENT(OUT):: PSFSV ! surface flux of scalar (m/s*kg/kg)
! flux of chemical var. (ppp.m/s)
REAL, DIMENSION(:,:), INTENT(OUT) :: PSFCO2! surface flux of CO2 (m/s*kg/kg)
@@ -227,46 +242,65 @@ REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZRV ! vapor mixing ratio
!
! suffix 'A' stands for atmospheric variable at first model level
!
-REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZZREF ! Forcing height
-REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZTA ! Temperature
-REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZRVA ! vapor mixing ratio
-REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZQA ! humidity (kg/m3)
-REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZPA ! Pressure
-REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZPS ! Pressure
-REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZEXNA ! Exner function
-REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZEXNS ! Exner function
-REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZTHA ! potential temperature
REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZRAIN ! liquid precipitation (kg/m2/s)
REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZSNOW ! solid precipitation (kg/m2/s)
REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZTSUN ! solar time (s since midnight)
-!
-REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZUA ! u component of the wind
-! ! parallel to the orography
-REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZVA ! v component of the wind
-! ! parallel to the orography
-REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZU ! zonal wind
-REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZV ! meridian wind
-REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZWIND ! wind parallel to the orography
-REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZRHOA ! air density
-REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZDIR ! wind direction (rad from N clockwise)
-REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZSFU ! zonal momentum flux
-REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZSFV ! meridian momentum flux
+REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZPS ! Surface pressure
+REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZEXNS ! Surface Exner function
REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZCO2 ! CO2 concentration (kg/kg)
-REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZALFA ! angle between the wind
-! ! and the x axis
-REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2),1):: ZU2D ! u and v component of the
-REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2),1):: ZV2D ! wind at mass point
-REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZSFTH ! Turbulent flux of heat
-REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZSFTQ ! Turbulent flux of water
-REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZSFCO2 ! Turbulent flux of CO2
-REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2),NSV):: ZSFTS! Turbulent flux of scalar
!
+! Variables for which multiple levels are sent to SURFEX and related ancilliary variables
+!
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZZREF ! Forcing height
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZTA ! Temperature
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZRVA ! vapor mixing ratio
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZQA ! humidity (kg/m3)
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZPA ! Pressure
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZEXNA ! Exner function
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZTHA ! potential temperature
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZUA ! u component of the wind parallel to the orography
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZVA ! v component of the wind parallel to the orography
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZU ! zonal wind
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZV ! meridian wind
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZWIND ! wind parallel to the orography
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZRHOA ! air density
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZTKE ! Subgrid turbulent kinetic energy
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZDIR ! wind direction (rad from N clockwise)
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZALFA ! angle between the wind and the x axis
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZU2D ! u and v component of the
+REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZV2D ! wind at mass point
+!
+! SURFEX output fluxes
+!
+REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZSFU ! zonal momentum flux
+REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZSFV ! meridian momentum flux
+REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZSFTH ! Total turbulent flux of heat
+REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZSFTH_SURF ! Surface turbulent flux of heat
+REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZSFTH_WALL ! Wall turbulent flux of heat
+REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZSFTH_ROOF ! Roof turbulent flux of heat
+REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZCD_ROOF ! Drag coefficient for roofs
+REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZSFTQ ! Total turbulent flux of water
+REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZSFTQ_SURF ! Surface turbulent flux of water
+REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZSFTQ_WALL ! Wall turbulent flux of water
+REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZSFTQ_ROOF ! Roof turbulent flux of water
+REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZSFCO2 ! Turbulent flux of CO2
+REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2),NSV):: ZSFTS ! Turbulent flux of scalar
REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2),NBLOWSNOW_2D) :: ZBLOWSNOW_2D ! 2D blowing snow variables
! after advection
! They refer to the 2D fields advected by MNH including:
! - total number concentration in Canopy
! - total mass concentration in Canopy
! - equivalent concentration in the saltation layer
+
+!
+! Anxiliary variables
+!
+REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZZREF_DIST
+REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZZREF_VERT
+REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZWEIGHT_VERT
+REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZAGLW_ILEV
+REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZAGLW_ILEVP1
+REAL, DIMENSION(SIZE(PSFTH,1),SIZE(PSFTH,2)) :: ZAGLSCAL_ILEV
!
!* Dimensions
! ----------
@@ -293,29 +327,44 @@ INTEGER :: KSV_SURF ! Number of scalar variables sent to SURFEX
!* Arrays put in 1D vectors
! ------------------------
!
-REAL, DIMENSION(:), ALLOCATABLE :: ZP_TSUN ! solar time
-REAL, DIMENSION(:), ALLOCATABLE :: ZP_ZENITH ! zenithal angle
-REAL, DIMENSION(:), ALLOCATABLE :: ZP_AZIM ! azimuthal angle
-REAL, DIMENSION(:), ALLOCATABLE :: ZP_ZREF ! forcing height
-REAL, DIMENSION(:), ALLOCATABLE :: ZP_ZS ! orography
-REAL, DIMENSION(:), ALLOCATABLE :: ZP_U ! zonal wind
-REAL, DIMENSION(:), ALLOCATABLE :: ZP_V ! meridian wind
-REAL, DIMENSION(:), ALLOCATABLE :: ZP_QA ! air humidity (kg/m3)
-REAL, DIMENSION(:), ALLOCATABLE :: ZP_TA ! air temperature
-REAL, DIMENSION(:), ALLOCATABLE :: ZP_RHOA ! air density
-REAL, DIMENSION(:,:), ALLOCATABLE :: ZP_SV ! scalar at first atmospheric level
+! Pure surface variables or variables forced at only one level
+!
REAL, DIMENSION(:), ALLOCATABLE :: ZP_CO2 ! air CO2 concentration
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZP_SV ! scalar at first atmospheric level
REAL, DIMENSION(:), ALLOCATABLE :: ZP_RAIN ! liquid precipitation
REAL, DIMENSION(:), ALLOCATABLE :: ZP_SNOW ! solid precipitation
REAL, DIMENSION(:), ALLOCATABLE :: ZP_LW ! incoming longwave
REAL, DIMENSION(:,:), ALLOCATABLE :: ZP_DIR_SW ! direct incoming shortwave
REAL, DIMENSION(:,:), ALLOCATABLE :: ZP_SCA_SW ! diffuse incoming shortwave
-REAL, DIMENSION(:), ALLOCATABLE :: ZP_PS ! surface pressure
-REAL, DIMENSION(:), ALLOCATABLE :: ZP_PA ! pressure at first atmospheric level
REAL, DIMENSION(:), ALLOCATABLE :: ZP_ZWS ! significant wave height (m)
-
-REAL, DIMENSION(:), ALLOCATABLE :: ZP_SFTQ ! water vapor flux
-REAL, DIMENSION(:), ALLOCATABLE :: ZP_SFTH ! potential temperature flux
+REAL, DIMENSION(:), ALLOCATABLE :: ZP_PS ! surface pressure
+REAL, DIMENSION(:), ALLOCATABLE :: ZP_TSUN ! solar time
+REAL, DIMENSION(:), ALLOCATABLE :: ZP_ZENITH ! zenithal angle
+REAL, DIMENSION(:), ALLOCATABLE :: ZP_AZIM ! azimuthal angle
+REAL, DIMENSION(:), ALLOCATABLE :: ZP_ZS ! orography
+!
+! Variables that are forced at multiple levels
+!
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZP_ZREF ! forcing height
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZP_U ! zonal wind
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZP_V ! meridian wind
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZP_QA ! air humidity (kg/m3)
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZP_TA ! air temperature
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZP_RHOA ! air density
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZP_PA ! pressure at first atmospheric level
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZP_TKE ! Subgrid turbulent kinetic energy
+!
+! SURFEX output variables
+!
+REAL, DIMENSION(:), ALLOCATABLE :: ZP_SFTQ ! Total water vapor flux
+REAL, DIMENSION(:), ALLOCATABLE :: ZP_SFTQ_SURF ! Surface water vapor flux
+REAL, DIMENSION(:), ALLOCATABLE :: ZP_SFTQ_WALL ! Wall water vapor flux
+REAL, DIMENSION(:), ALLOCATABLE :: ZP_SFTQ_ROOF ! Roof water vapor flux
+REAL, DIMENSION(:), ALLOCATABLE :: ZP_SFTH ! Total potential temperature flux
+REAL, DIMENSION(:), ALLOCATABLE :: ZP_SFTH_SURF ! Surface potential temperature flux
+REAL, DIMENSION(:), ALLOCATABLE :: ZP_SFTH_WALL ! Wall potential temperature flux
+REAL, DIMENSION(:), ALLOCATABLE :: ZP_SFTH_ROOF ! Roof potential temperature flux
+REAL, DIMENSION(:), ALLOCATABLE :: ZP_CD_ROOF ! Drag coefficient for roofs
REAL, DIMENSION(:,:), ALLOCATABLE :: ZP_SFTS ! scalar flux
REAL, DIMENSION(:), ALLOCATABLE :: ZP_SFCO2 ! CO2 flux
REAL, DIMENSION(:), ALLOCATABLE :: ZP_SFU ! zonal momentum flux
@@ -324,12 +373,11 @@ REAL, DIMENSION(:), ALLOCATABLE :: ZP_TSRAD ! radiative surface temperature
REAL, DIMENSION(:,:), ALLOCATABLE :: ZP_DIR_ALB ! direct albedo
REAL, DIMENSION(:,:), ALLOCATABLE :: ZP_SCA_ALB ! diffuse albedo
REAL, DIMENSION(:), ALLOCATABLE :: ZP_EMIS ! emissivity
-
REAL, DIMENSION(:), ALLOCATABLE :: ZP_TSURF
REAL, DIMENSION(:), ALLOCATABLE :: ZP_Z0
REAL, DIMENSION(:), ALLOCATABLE :: ZP_Z0H
REAL, DIMENSION(:), ALLOCATABLE :: ZP_QSURF
-
+!
REAL, DIMENSION(:), ALLOCATABLE :: ZP_PEW_A_COEF ! coefficients for
REAL, DIMENSION(:), ALLOCATABLE :: ZP_PEW_B_COEF ! implicit coupling
REAL, DIMENSION(:), ALLOCATABLE :: ZP_PET_A_COEF
@@ -359,6 +407,13 @@ CHARACTER(LEN=6), DIMENSION(:), ALLOCATABLE :: YSV_SURF ! name of the scalar var
REAL :: ZTIMEC
INTEGER :: ILUOUT ! logical unit
!
+! New variables for coupling at several levels
+!
+REAL :: ZAGLW_JK
+REAL :: ZAGLW_JKP1
+REAL :: ZAGLSCAL_JK
+INTEGER :: ICOUNT, ILEV
+!
!-------------------------------------------------------------------------------
!
!
@@ -369,8 +424,14 @@ IKE=IKU-JPVEXT
!
CALL GET_INDICE_ll (IIB,IJB,IIE,IJE)
!
-PSFTH = XUNDEF
-PSFRV = XUNDEF
+PSFTH = XUNDEF
+PSFTH_WALL = XUNDEF
+PSFTH_ROOF = XUNDEF
+PCD_ROOF = XUNDEF
+PSFRV = XUNDEF
+PSFRV_WALL = XUNDEF
+PSFRV_ROOF = XUNDEF
+!
PSFSV = XUNDEF
PSFCO2 = XUNDEF
PSFU = XUNDEF
@@ -380,6 +441,26 @@ PSCA_ALB = XUNDEF
PEMIS = XUNDEF
PTSRAD = XUNDEF
!
+! Allocation of the local variables
+!
+ALLOCATE(ZZREF(SIZE(PSFTH,1),SIZE(PSFTH,2),NLEV_COUPLE))
+ALLOCATE(ZTA(SIZE(PSFTH,1),SIZE(PSFTH,2),NLEV_COUPLE))
+ALLOCATE(ZRVA(SIZE(PSFTH,1),SIZE(PSFTH,2),NLEV_COUPLE))
+ALLOCATE(ZQA(SIZE(PSFTH,1),SIZE(PSFTH,2),NLEV_COUPLE))
+ALLOCATE(ZPA(SIZE(PSFTH,1),SIZE(PSFTH,2),NLEV_COUPLE))
+ALLOCATE(ZEXNA(SIZE(PSFTH,1),SIZE(PSFTH,2),NLEV_COUPLE))
+ALLOCATE(ZTHA(SIZE(PSFTH,1),SIZE(PSFTH,2),NLEV_COUPLE))
+ALLOCATE(ZUA(SIZE(PSFTH,1),SIZE(PSFTH,2),NLEV_COUPLE))
+ALLOCATE(ZVA(SIZE(PSFTH,1),SIZE(PSFTH,2),NLEV_COUPLE))
+ALLOCATE(ZU(SIZE(PSFTH,1),SIZE(PSFTH,2),NLEV_COUPLE))
+ALLOCATE(ZV(SIZE(PSFTH,1),SIZE(PSFTH,2),NLEV_COUPLE))
+ALLOCATE(ZWIND(SIZE(PSFTH,1),SIZE(PSFTH,2),NLEV_COUPLE))
+ALLOCATE(ZRHOA(SIZE(PSFTH,1),SIZE(PSFTH,2),NLEV_COUPLE))
+ALLOCATE(ZTKE(SIZE(PSFTH,1),SIZE(PSFTH,2),NLEV_COUPLE))
+ALLOCATE(ZDIR(SIZE(PSFTH,1),SIZE(PSFTH,2),NLEV_COUPLE))
+ALLOCATE(ZALFA(SIZE(PSFTH,1),SIZE(PSFTH,2),NLEV_COUPLE))
+ALLOCATE(ZU2D(SIZE(PSFTH,1),SIZE(PSFTH,2),NLEV_COUPLE))
+ALLOCATE(ZV2D(SIZE(PSFTH,1),SIZE(PSFTH,2),NLEV_COUPLE))
!
!-------------------------------------------------------------------------------
!
@@ -401,52 +482,78 @@ END IF
! 1.2 Horizontal wind direction (rad from N clockwise)
! -------------------------
!
-ZU2D(:,:,:)=MXF(XUT(:,:,IKB:IKB))
-ZV2D(:,:,:)=MYF(XVT(:,:,IKB:IKB))
+ZU2D(:,:,:)=MXF(XUT(:,:,IKB:(IKB+NLEV_COUPLE-1)))
+ZV2D(:,:,:)=MYF(XVT(:,:,IKB:(IKB+NLEV_COUPLE-1)))
!
!* angle between Y axis and wind (rad., clockwise)
!
ZALFA = 0.
-WHERE(ZU2D(:,:,1)/=0. .OR. ZV2D(:,:,1)/=0.)
- ZALFA(:,:)=ATAN2(ZU2D(:,:,1),ZV2D(:,:,1))
-END WHERE
-WHERE(ZALFA(:,:)<0.) ZALFA(:,:) = ZALFA(:,:) + 2. * XPI
!
-!* angle between North and wind (rad., clockwise)
-!
-IF (.NOT. LCARTESIAN) THEN
- ZDIR = ( (XRPK*(XLON(:,:)-XLON0)) - XBETA ) * XPI/180. + ZALFA
+DO ILEV=1,NLEV_COUPLE
+ !
+ WHERE(ZU2D(:,:,ILEV)/=0. .OR. ZV2D(:,:,ILEV)/=0.)
+ ZALFA(:,:,ILEV)=ATAN2(ZU2D(:,:,ILEV),ZV2D(:,:,ILEV))
+ END WHERE
+ !
+ WHERE(ZALFA(:,:,ILEV)<0.) ZALFA(:,:,ILEV) = ZALFA(:,:,ILEV) + 2. * XPI
+ !
+ !* angle between North and wind (rad., clockwise)
+ !
+ IF (.NOT. LCARTESIAN) THEN
+ ZDIR(:,:,ILEV) = ( (XRPK*(XLON(:,:)-XLON0)) - XBETA ) * XPI/180. + ZALFA(:,:,ILEV)
+ ELSE
+ ZDIR(:,:,ILEV) = - XBETA * XPI/180. + ZALFA(:,:,ILEV)
+ ENDIF
+ !
+ ! 1.3 Rotate the wind
+ ! Only for the first forcing level, used for friction force direction.
+ ! ---------------
+ !
+ IF (ILEV.EQ.1) THEN
+ !
+ CALL ROTATE_WIND(XUT,XVT,XWT, &
+ XDIRCOSXW, XDIRCOSYW, XDIRCOSZW, &
+ XCOSSLOPE,XSINSLOPE, &
+ XDXX,XDYY,XDZZ, &
+ ZUA(:,:,ILEV),ZVA(:,:,ILEV) )
+ !
+ ELSE
+ !
+ ZUA(:,:,ILEV) = XUT(:,:,IKB+ILEV-1)
+ ZVA(:,:,ILEV) = XVT(:,:,IKB+ILEV-1)
+ !
+ ENDIF
+ !
+ ! 1.4 zonal and meridian components of the wind parallel to the slope
+ ! ---------------------------------------------------------------
+ !
+ ZWIND(:,:,ILEV) = SQRT( ZUA(:,:,ILEV)**2 + ZVA(:,:,ILEV)**2 )
+ !
+ ZU(:,:,ILEV) = ZWIND(:,:,ILEV) * SIN(ZDIR(:,:,ILEV))
+ ZV(:,:,ILEV) = ZWIND(:,:,ILEV) * COS(ZDIR(:,:,ILEV))
+ !
+ENDDO
+ !
+ ! 1.5 Horizontal interpolation of the thermodynamic fields
+ ! -------------------------------------------------
+ !
+ ! This horizontal interpolation is only made if the forcing is located at the first level
+ !
+IF (NLEV_COUPLE.EQ.1) THEN
+ !
+ CALL NORMAL_INTERPOL(XTHT,ZRV,XPABST, &
+ XDIRCOSXW, XDIRCOSYW, XDIRCOSZW, &
+ XCOSSLOPE,XSINSLOPE, &
+ XDXX,XDYY,XDZZ, &
+ ZTHA(:,:,1),ZRVA(:,:,1),ZEXNA(:,:,1) )
+ !
ELSE
- ZDIR = - XBETA * XPI/180. + ZALFA
-END IF
-!
-!
-! 1.3 Rotate the wind
-! ---------------
-!
-CALL ROTATE_WIND(XUT,XVT,XWT, &
- XDIRCOSXW, XDIRCOSYW, XDIRCOSZW, &
- XCOSSLOPE,XSINSLOPE, &
- XDXX,XDYY,XDZZ, &
- ZUA,ZVA )
-
-!
-! 1.4 zonal and meridian components of the wind parallel to the slope
-! ---------------------------------------------------------------
-!
-ZWIND(:,:) = SQRT( ZUA**2 + ZVA**2 )
-!
-ZU(:,:) = ZWIND(:,:) * SIN(ZDIR)
-ZV(:,:) = ZWIND(:,:) * COS(ZDIR)
-!
-! 1.5 Horizontal interpolation the thermodynamic fields
-! -------------------------------------------------
-!
-CALL NORMAL_INTERPOL(XTHT,ZRV,XPABST, &
- XDIRCOSXW, XDIRCOSYW, XDIRCOSZW, &
- XCOSSLOPE,XSINSLOPE, &
- XDXX,XDYY,XDZZ, &
- ZTHA,ZRVA,ZEXNA )
+ !
+ ZEXNA (:,:,1:NLEV_COUPLE) = (XPABST(:,:,IKB:(IKB+NLEV_COUPLE-1))/XP00) ** (XRD/XCPD)
+ ZTHA (:,:,1:NLEV_COUPLE) = XTHT(:,:,IKB:(IKB+NLEV_COUPLE-1))
+ ZRVA (:,:,1:NLEV_COUPLE) = ZRV (:,:,IKB:(IKB+NLEV_COUPLE-1))
+ !
+ENDIF
!
DEALLOCATE(ZRV)
!
@@ -454,8 +561,7 @@ DEALLOCATE(ZRV)
! 1.6 Pressure and Exner function
! ---------------------------
!
-!
-ZPA(:,:) = XP00 * ZEXNA(:,:) **(XCPD/XRD)
+ZPA(:,:,:) = XP00 * ZEXNA(:,:,:) ** (XCPD/XRD)
!
ZEXNS(:,:) = 0.5 * ( (XPABST(:,:,IKB-1)/XP00)**(XRD/XCPD) &
+(XPABST(:,:,IKB )/XP00)**(XRD/XCPD) &
@@ -465,23 +571,22 @@ ZPS(:,:) = XP00 * ZEXNS(:,:) **(XCPD/XRD)
! 1.7 humidity in kg/m3 from the mixing ratio
! ---------------------------------------
!
-!
-ZQA(:,:) = ZRVA(:,:) * XRHODREF(:,:,IKB)
-!
+ZQA(:,:,:) = ZRVA(:,:,:) * XRHODREF(:,:,IKB:(IKB+NLEV_COUPLE-1))
!
! 1.8 Temperature from the potential temperature
! ------------------------------------------
!
-!
-ZTA(:,:) = ZTHA(:,:) * ZEXNA(:,:)
-!
+ZTA(:,:,:) = ZTHA(:,:,:) * ZEXNA(:,:,:)
!
! 1.9 Air density
! -----------
!
-ZRHOA(:,:) = ZPA(:,:)/(XRD * ZTA(:,:) * ((1. + (XRD/XRV)*ZRVA(:,:))/ &
- (1. + ZRVA(:,:))))
+ZRHOA(:,:,:) = ZPA(:,:,:)/(XRD * ZTA(:,:,:) * &
+ ((1. + (XRD/XRV)*ZRVA(:,:,:)) / (1. + ZRVA(:,:,:))))
+!
+! Subgrid turbulent kinetic energy
!
+ZTKE(:,:,:) = XTKET(:,:,IKB:(IKB+NLEV_COUPLE-1))
!
! 1.10 Precipitations
! --------------
@@ -518,8 +623,39 @@ END IF
! 1.12 Forcing level
! -------------
!
-ZZREF(:,:) = 0.5*( XZZ(:,:,IKB+1)-XZZ(:,:,IKB) )*XDIRCOSZW(:,:)
-!
+! A smooth transition between vertical height above ground and
+! distance to the surface is implemented here.
+! We assume that for katabatic winds located in the first meters above
+! ground, the distance to the surface is the most relevant whereas
+! for most other processes it will be the vertical distance to the surface
+!
+DO ILEV=1,NLEV_COUPLE
+ !
+ ! Height above ground of w-levels
+ !
+ ZAGLW_ILEV (:,:) = XZZ(:,:,JPVEXT+ILEV ) - XZZ(:,:,1+JPVEXT)
+ ZAGLW_ILEVP1 (:,:) = XZZ(:,:,JPVEXT+ILEV+1) - XZZ(:,:,1+JPVEXT)
+ !
+ ! Height above ground of scalar variables and (u,v)
+ !
+ ZAGLSCAL_ILEV(:,:) = 0.5 * ( ZAGLW_ILEV(:,:) + ZAGLW_ILEVP1(:,:) )
+ !
+ ! Distance to the inclined surface and vertical distance
+ !
+ ZZREF_DIST(:,:) = ZAGLSCAL_ILEV(:,:) * XDIRCOSZW(:,:)
+ !
+ ZZREF_VERT(:,:) = ZAGLSCAL_ILEV(:,:)
+ !
+ ! Scaling between 5 m and 20 m height
+ !
+ ZWEIGHT_VERT(:,:) = MIN(1.0,MAX(ZZREF_VERT(:,:)-5.0,0.0)/15.0)
+ !
+ IF (MAXVAL(ZWEIGHT_VERT).GT.1.0) STOP ("Wrong weight")
+ IF (MINVAL(ZWEIGHT_VERT).LT.0.0) STOP ("Wrong weight")
+ !
+ ZZREF(:,:,ILEV) = ZWEIGHT_VERT(:,:) * ZZREF_VERT(:,:) + (1.0 - ZWEIGHT_VERT(:,:)) * ZZREF_DIST(:,:)
+ !
+ENDDO
!
! 1.13 CO2 concentration (kg/m3)
! -----------------
@@ -595,19 +731,43 @@ END IF
#endif
!
! Call to surface schemes
-!
-CALL COUPLING_SURF_ATM_n(YSURF_CUR,'MESONH', 'E',ZTIMEC, &
- XTSTEP, TDTCUR%nyear, TDTCUR%nmonth, TDTCUR%nday, TDTCUR%xtime, &
- IDIM1D,KSV_SURF,SIZE(XSW_BANDS), &
- ZP_TSUN, ZP_ZENITH,ZP_ZENITH, ZP_AZIM, &
- ZP_ZREF, ZP_ZREF, ZP_ZS, ZP_U, ZP_V, ZP_QA, ZP_TA, ZP_RHOA, ZP_SV, ZP_CO2, &
- ZP_ZIMPWET, ZP_ZIMPDRY, YSV_SURF, &
- ZP_RAIN, ZP_SNOW, ZP_LW, ZP_DIR_SW, ZP_SCA_SW, XSW_BANDS, ZP_PS, ZP_PA, &
- ZP_SFTQ, ZP_SFTH, ZP_SFTS, ZP_SFCO2, ZP_SFU, ZP_SFV, &
- ZP_TSRAD, ZP_DIR_ALB, ZP_SCA_ALB, ZP_EMIS, ZP_TSURF, ZP_Z0, ZP_Z0H, ZP_QSURF, &
- ZP_PEW_A_COEF, ZP_PEW_B_COEF, &
- ZP_PET_A_COEF, ZP_PEQ_A_COEF, ZP_PET_B_COEF, ZP_PEQ_B_COEF,ZP_ZWS, &
- 'OK' )
+!
+CALL COUPLING_SURF_ATM_MULTI_LEVEL_n(YSURF_CUR,'MESONH', 'E',ZTIMEC, XTSTEP, &
+ TDTCUR%nyear, TDTCUR%nmonth, TDTCUR%nday, TDTCUR%xtime, &
+ IDIM1D,KSV_SURF,SIZE(XSW_BANDS), NLEV_COUPLE, ZP_TSUN, ZP_ZENITH,ZP_ZENITH, &
+ ZP_AZIM, ZP_ZREF, ZP_ZREF, ZP_ZS, ZP_U, ZP_V, ZP_QA, ZP_TA, ZP_RHOA, ZP_SV, &
+ ZP_CO2, ZP_ZIMPWET, ZP_ZIMPDRY, YSV_SURF, &
+ ZP_RAIN, ZP_SNOW, ZP_LW, ZP_DIR_SW, ZP_SCA_SW, XSW_BANDS, &
+ ZP_PS, ZP_PA, ZP_TKE, ZP_SFTQ, ZP_SFTQ_SURF, ZP_SFTQ_WALL, ZP_SFTQ_ROOF, &
+ ZP_SFTH, ZP_SFTH_SURF, ZP_SFTH_WALL, ZP_SFTH_ROOF, ZP_CD_ROOF, ZP_SFTS, &
+ ZP_SFCO2, ZP_SFU, ZP_SFV, ZP_TSRAD, ZP_DIR_ALB, ZP_SCA_ALB, ZP_EMIS, ZP_TSURF, &
+ ZP_Z0, ZP_Z0H, ZP_QSURF, ZP_PEW_A_COEF, ZP_PEW_B_COEF, ZP_PET_A_COEF, &
+ ZP_PEQ_A_COEF, ZP_PET_B_COEF, ZP_PEQ_B_COEF, ZP_ZWS, 'OK' )
+
+
+
+!Jean Wurtz
+!In some cases LE and H are not realistic at the beginning
+!Hard fix, source problem may have to be corrected
+!Energy is hence not conservated at the beginning...
+
+ZP_SFTH(:)=MAX(MIN(ZP_SFTH(:),1500.),-1500.)
+ZP_SFTH_WALL(:)=MAX(MIN(ZP_SFTH_WALL(:),1500.),-1500.)
+ZP_SFTH_ROOF(:)=MAX(MIN(ZP_SFTH_ROOF(:),1500.),-1500.)
+ZP_SFTH_SURF(:)=MAX(MIN(ZP_SFTH_SURF(:),1500.),-1500.)
+
+
+! Enthalpie de vaporisation : 2.26476E6
+! on divise 1500 par cette enthalpie
+
+ZP_SFTQ(:)=MAX(MIN(ZP_SFTQ(:),0.000662322),-0.000662322)
+ZP_SFTQ_WALL(:)=MAX(MIN(ZP_SFTQ_WALL(:),0.000662322),-0.000662322)
+ZP_SFTQ_ROOF(:)=MAX(MIN(ZP_SFTQ_ROOF(:),0.000662322),-0.000662322)
+ZP_SFTQ_SURF(:)=MAX(MIN(ZP_SFTQ_SURF(:),0.000662322),-0.000662322)
+
+
+
+
!
#ifdef CPLOASIS
IF (LOASIS) THEN
@@ -625,9 +785,9 @@ END IF
!
IF (CPROGRAM=='DIAG ' .OR. LDIAG_IN_RUN) THEN
CALL DIAG_SURF_ATM_n(YSURF_CUR,'MESONH')
- CALL MNHGET_SURF_PARAM_n( PRN = ZP_RN, PH = ZP_H, PLE = ZP_LE, PLEI = ZP_LEI, &
- PGFLUX = ZP_GFLUX, PT2M = ZP_T2M, PQ2M = ZP_Q2M, PHU2M = ZP_HU2M, &
- PZON10M = ZP_ZON10M, PMER10M = ZP_MER10M )
+ CALL MNHGET_SURF_PARAM_n(PRN=ZP_RN,PH=ZP_H,PLE=ZP_LE,PLEI=ZP_LEI, &
+ PGFLUX=ZP_GFLUX,PT2M=ZP_T2M,PQ2M=ZP_Q2M,PHU2M=ZP_HU2M, &
+ PZON10M=ZP_ZON10M,PMER10M=ZP_MER10M )
END IF
!
! Transform 1D output fields into 2D:
@@ -637,7 +797,7 @@ CALL UNSHAPE_SURF(IDIM1,IDIM2)
!------------------------!
! COUPLING WITH FOREFIRE !
!------------------------!
-
+
IF ( LFOREFIRE ) THEN
CALL FOREFIRE_DUMP_FIELDS_n(XUT, XVT, XWT, XSVT&
, XTHT, XRT(:,:,:,1), XPABST, XTKET&
@@ -661,24 +821,59 @@ FF_TIME = FF_TIME + XTSTEP
!
! Friction of components along slope axes (U: largest local slope axis, V: zero slope axis)
!
-!
PSFU(:,:) = 0.
PSFV(:,:) = 0.
!
-WHERE (ZSFU(:,:)/=XUNDEF .AND. ZWIND(:,:)>0.)
- PSFU(:,:) = - SQRT(ZSFU**2+ZSFV**2) * ZUA(:,:) / ZWIND(:,:) / XRHODREF(:,:,IKB)
- PSFV(:,:) = - SQRT(ZSFU**2+ZSFV**2) * ZVA(:,:) / ZWIND(:,:) / XRHODREF(:,:,IKB)
+WHERE (ZSFU(:,:)/=XUNDEF .AND. ZWIND(:,:,1)>0.)
+ PSFU(:,:) = - SQRT(ZSFU**2+ZSFV**2) * ZUA(:,:,1) / ZWIND(:,:,1) / XRHODREF(:,:,IKB)
+ PSFV(:,:) = - SQRT(ZSFU**2+ZSFV**2) * ZVA(:,:,1) / ZWIND(:,:,1) / XRHODREF(:,:,IKB)
END WHERE
!
-!* conversion from H (W/m2) to w'Theta'
-!
-PSFTH(:,:) = ZSFTH(:,:) / XCPD / XRHODREF(:,:,IKB)
-!
-!
-!* conversion from water flux (kg/m2/s) to w'rv'
-!
-PSFRV(:,:) = ZSFTQ(:,:) / XRHODREF(:,:,IKB)
-!
+PCD_ROOF(:,:) = ZCD_ROOF(:,:)
+!
+! Unit conversions:
+!
+!* H: (W/m2) to w'Theta'
+!
+!* Water flux: (kg/m2/s) to w'rv'
+!
+IF (LFLUXBLDG) THEN
+ !
+ ! Robert: Here the wall and roof fluxes are substracted from the surface fluxes
+ ! since they will be applied in drag_bld.F90
+ !
+ PSFTH(:,:) = ( ZSFTH(:,:) - ZSFTH_WALL(:,:) - ZSFTH_ROOF(:,:) ) / XCPD / XRHODREF(:,:,IKB)
+ PSFRV(:,:) = ( ZSFTQ(:,:) - ZSFTQ_WALL(:,:) - ZSFTQ_ROOF(:,:) ) / XRHODREF(:,:,IKB)
+ !
+ ! Wall and roof fluxes are written on separate variables
+ !
+ PSFTH_WALL(:,:) = ZSFTH_WALL(:,:) / XCPD / XRHODREF(:,:,IKB)
+ PSFTH_ROOF(:,:) = ZSFTH_ROOF(:,:) / XCPD / XRHODREF(:,:,IKB)
+ !
+ PSFRV_WALL(:,:) = ZSFTQ_WALL(:,:) / XRHODREF(:,:,IKB)
+ PSFRV_ROOF(:,:) = ZSFTQ_ROOF(:,:) / XRHODREF(:,:,IKB)
+ !
+ ! Test conservation of fluxes
+ !
+ IF (MAXVAL(ABS(ZSFTH(:,:)/XCPD/XRHODREF(:,:,IKB) - PSFTH(:,:) - PSFTH_WALL(:,:)&
+ - PSFTH_ROOF(:,:))).GT.1.0E-6) STOP ("Wrong H flux partition")
+ IF (MAXVAL(ABS(ZSFTQ(:,:)/XRHODREF(:,:,IKB) - PSFRV(:,:) - PSFRV_WALL(:,:)&
+ - PSFRV_ROOF(:,:))).GT.1.0E-6) STOP ("Wrong Q flux partition")
+ !
+ELSE
+ !
+ ! Otherwise the full surface fluxes are taken
+ !
+ PSFTH(:,:) = ZSFTH(:,:) / XCPD / XRHODREF(:,:,IKB)
+ PSFRV(:,:) = ZSFTQ(:,:) / XRHODREF(:,:,IKB)
+ !
+ PSFTH_WALL(:,:) = 0.0
+ PSFTH_ROOF(:,:) = 0.0
+ !
+ PSFRV_WALL(:,:) = 0.0
+ PSFRV_ROOF(:,:) = 0.0
+ !
+ENDIF
!
!* conversion from scalar flux (kg/m2/s) to w'rsv'
!
@@ -733,11 +928,11 @@ END IF
!
IF (LBLOWSNOW) THEN
DO JSV=NSV_SNWBEG,NSV_SNWEND
- PSFSV(:,:,JSV) = ZSFTS(:,:,JSV)/ (ZRHOA(:,:))
+ PSFSV(:,:,JSV) = ZSFTS(:,:,JSV)/ (ZRHOA(:,:,1))
END DO
!* Update tendency for blowing snow 2D fields
DO JSV=1,(NBLOWSNOW_2D)
- XRSNWCANOS(:,:,JSV) = ZBLOWSNOW_2D(:,:,JSV)*XRHODJ(:,:,IKB)/(XTSTEP*ZRHOA(:,:))
+ XRSNWCANOS(:,:,JSV) = ZBLOWSNOW_2D(:,:,JSV)*XRHODJ(:,:,IKB)/(XTSTEP*ZRHOA(:,:,1))
END DO
ELSE
@@ -775,7 +970,7 @@ IF (LDIAG_IN_RUN) THEN
ENDDO
ENDDO
END IF
-!
+ !
NULLIFY(TZFIELDSURF_ll)
CALL ADD2DFIELD_ll( TZFIELDSURF_ll,XCURRENT_RN, 'GROUND_PARAM_n::XCURRENT_RN' )
CALL ADD2DFIELD_ll( TZFIELDSURF_ll,XCURRENT_H, 'GROUND_PARAM_n::XCURRENT_H' )
@@ -794,9 +989,10 @@ IF (LDIAG_IN_RUN) THEN
CALL ADD2DFIELD_ll( TZFIELDSURF_ll,XCURRENT_SLTAOD, 'GROUND_PARAM_n::XCURRENT_SLTAOD' )
CALL ADD2DFIELD_ll( TZFIELDSURF_ll,XCURRENT_ZWS, 'GROUND_PARAM_n::XCURRENT_ZWS' )
CALL ADD2DFIELD_ll( TZFIELDSURF_ll,XCURRENT_SFCO2, 'GROUND_PARAM_n::XCURRENT_SFCO2' )
-
+ !
CALL UPDATE_HALO_ll(TZFIELDSURF_ll,IINFO_ll)
CALL CLEANLIST_ll(TZFIELDSURF_ll)
+ !
END IF
!
!==================================================================================
@@ -812,16 +1008,23 @@ INTEGER, DIMENSION(1) :: ISHAPE_1
!
ISHAPE_1 = (/KDIM1D/)
!
+! Variables that are coupled at multiple levels
+!
+ALLOCATE(ZP_ZREF (KDIM1D,NLEV_COUPLE))
+ALLOCATE(ZP_U (KDIM1D,NLEV_COUPLE))
+ALLOCATE(ZP_V (KDIM1D,NLEV_COUPLE))
+ALLOCATE(ZP_QA (KDIM1D,NLEV_COUPLE))
+ALLOCATE(ZP_TA (KDIM1D,NLEV_COUPLE))
+ALLOCATE(ZP_PA (KDIM1D,NLEV_COUPLE))
+ALLOCATE(ZP_RHOA (KDIM1D,NLEV_COUPLE))
+ALLOCATE(ZP_TKE (KDIM1D,NLEV_COUPLE))
+!
+! 2D Variables and variables that are coupled at the surface only
+!
ALLOCATE(ZP_TSUN (KDIM1D))
ALLOCATE(ZP_ZENITH (KDIM1D))
ALLOCATE(ZP_AZIM (KDIM1D))
-ALLOCATE(ZP_ZREF (KDIM1D))
ALLOCATE(ZP_ZS (KDIM1D))
-ALLOCATE(ZP_U (KDIM1D))
-ALLOCATE(ZP_V (KDIM1D))
-ALLOCATE(ZP_QA (KDIM1D))
-ALLOCATE(ZP_TA (KDIM1D))
-ALLOCATE(ZP_RHOA (KDIM1D))
ALLOCATE(ZP_SV (KDIM1D,KSV_SURF))
ALLOCATE(ZP_CO2 (KDIM1D))
ALLOCATE(ZP_RAIN (KDIM1D))
@@ -830,11 +1033,19 @@ ALLOCATE(ZP_LW (KDIM1D))
ALLOCATE(ZP_DIR_SW (KDIM1D,SIZE(XDIRSRFSWD,3)))
ALLOCATE(ZP_SCA_SW (KDIM1D,SIZE(XSCAFLASWD,3)))
ALLOCATE(ZP_PS (KDIM1D))
-ALLOCATE(ZP_PA (KDIM1D))
ALLOCATE(ZP_ZWS (KDIM1D))
-
-ALLOCATE(ZP_SFTQ (KDIM1D))
-ALLOCATE(ZP_SFTH (KDIM1D))
+!
+! 2D SURFEX output fields
+!
+ALLOCATE(ZP_SFTQ (KDIM1D))
+ALLOCATE(ZP_SFTQ_SURF (KDIM1D))
+ALLOCATE(ZP_SFTQ_WALL (KDIM1D))
+ALLOCATE(ZP_SFTQ_ROOF (KDIM1D))
+ALLOCATE(ZP_SFTH (KDIM1D))
+ALLOCATE(ZP_SFTH_SURF (KDIM1D))
+ALLOCATE(ZP_SFTH_WALL (KDIM1D))
+ALLOCATE(ZP_SFTH_ROOF (KDIM1D))
+ALLOCATE(ZP_CD_ROOF (KDIM1D))
ALLOCATE(ZP_SFU (KDIM1D))
ALLOCATE(ZP_SFV (KDIM1D))
ALLOCATE(ZP_SFTS (KDIM1D,KSV_SURF))
@@ -857,7 +1068,7 @@ ALLOCATE(ZP_Q2M (KDIM1D))
ALLOCATE(ZP_HU2M (KDIM1D))
ALLOCATE(ZP_ZON10M (KDIM1D))
ALLOCATE(ZP_MER10M (KDIM1D))
-
+!
!* explicit coupling only
ALLOCATE(ZP_PEW_A_COEF (KDIM1D))
ALLOCATE(ZP_PEW_B_COEF (KDIM1D))
@@ -865,22 +1076,30 @@ ALLOCATE(ZP_PET_A_COEF (KDIM1D))
ALLOCATE(ZP_PEQ_A_COEF (KDIM1D))
ALLOCATE(ZP_PET_B_COEF (KDIM1D))
ALLOCATE(ZP_PEQ_B_COEF (KDIM1D))
-
-ZP_TSUN(:) = RESHAPE(ZTSUN(IIB:IIE,IJB:IJE), ISHAPE_1)
-ZP_TA(:) = RESHAPE(ZTA(IIB:IIE,IJB:IJE), ISHAPE_1)
-ZP_QA(:) = RESHAPE(ZQA(IIB:IIE,IJB:IJE), ISHAPE_1)
-ZP_RHOA(:) = RESHAPE(ZRHOA(IIB:IIE,IJB:IJE), ISHAPE_1)
-ZP_U(:) = RESHAPE(ZU(IIB:IIE,IJB:IJE), ISHAPE_1)
-ZP_V(:) = RESHAPE(ZV(IIB:IIE,IJB:IJE), ISHAPE_1)
-ZP_PS(:) = RESHAPE(ZPS(IIB:IIE,IJB:IJE), ISHAPE_1)
-ZP_PA(:) = RESHAPE(ZPA(IIB:IIE,IJB:IJE), ISHAPE_1)
-ZP_ZS(:) = RESHAPE(XZS(IIB:IIE,IJB:IJE), ISHAPE_1)
-ZP_CO2(:) = RESHAPE(ZCO2(IIB:IIE,IJB:IJE), ISHAPE_1)
-ZP_SNOW(:) = RESHAPE(ZSNOW(IIB:IIE,IJB:IJE), ISHAPE_1)
-ZP_RAIN(:) = RESHAPE(ZRAIN(IIB:IIE,IJB:IJE), ISHAPE_1)
-ZP_ZREF(:) = RESHAPE(ZZREF(IIB:IIE,IJB:IJE), ISHAPE_1)
-ZP_ZWS(:) = RESHAPE(XZWS(IIB:IIE,IJB:IJE), ISHAPE_1)
-
+!
+! 2D variables or surface only
+!
+ZP_TSUN(:) = RESHAPE(ZTSUN(IIB:IIE,IJB:IJE), ISHAPE_1)
+ZP_PS(:) = RESHAPE(ZPS(IIB:IIE,IJB:IJE), ISHAPE_1)
+ZP_ZS(:) = RESHAPE(XZS(IIB:IIE,IJB:IJE), ISHAPE_1)
+ZP_CO2(:) = RESHAPE(ZCO2(IIB:IIE,IJB:IJE), ISHAPE_1)
+ZP_SNOW(:) = RESHAPE(ZSNOW(IIB:IIE,IJB:IJE), ISHAPE_1)
+ZP_RAIN(:) = RESHAPE(ZRAIN(IIB:IIE,IJB:IJE), ISHAPE_1)
+ZP_ZWS(:) = RESHAPE(XZWS(IIB:IIE,IJB:IJE), ISHAPE_1)
+!
+! Variables that are coupled on multiple levels
+!
+DO JLAYER=1,NLEV_COUPLE
+ ZP_ZREF(:,JLAYER) = RESHAPE(ZZREF(IIB:IIE,IJB:IJE,JLAYER), ISHAPE_1)
+ ZP_PA(:,JLAYER) = RESHAPE(ZPA(IIB:IIE,IJB:IJE,JLAYER), ISHAPE_1)
+ ZP_TA(:,JLAYER) = RESHAPE(ZTA(IIB:IIE,IJB:IJE,JLAYER), ISHAPE_1)
+ ZP_QA(:,JLAYER) = RESHAPE(ZQA(IIB:IIE,IJB:IJE,JLAYER), ISHAPE_1)
+ ZP_RHOA(:,JLAYER) = RESHAPE(ZRHOA(IIB:IIE,IJB:IJE,JLAYER), ISHAPE_1)
+ ZP_TKE(:,JLAYER) = RESHAPE(ZTKE(IIB:IIE,IJB:IJE,JLAYER), ISHAPE_1)
+ ZP_U(:,JLAYER) = RESHAPE(ZU(IIB:IIE,IJB:IJE,JLAYER), ISHAPE_1)
+ ZP_V(:,JLAYER) = RESHAPE(ZV(IIB:IIE,IJB:IJE,JLAYER), ISHAPE_1)
+END DO
+!
DO JLAYER=1,NSV
ZP_SV(:,JLAYER) = RESHAPE(XSVT(IIB:IIE,IJB:IJE,IKB,JLAYER), ISHAPE_1)
END DO
@@ -893,29 +1112,29 @@ END IF
!
!chemical conversion : from part/part to molec./m3
DO JLAYER=NSV_CHEMBEG,NSV_CHEMEND
- ZP_SV(:,JLAYER) = ZP_SV(:,JLAYER) * XAVOGADRO * ZP_RHOA(:) / XMD
+ ZP_SV(:,JLAYER) = ZP_SV(:,JLAYER) * XAVOGADRO * ZP_RHOA(:,1) / XMD
END DO
DO JLAYER=NSV_AERBEG,NSV_AEREND
- ZP_SV(:,JLAYER) = ZP_SV(:,JLAYER) * XAVOGADRO * ZP_RHOA(:) / XMD
+ ZP_SV(:,JLAYER) = ZP_SV(:,JLAYER) * XAVOGADRO * ZP_RHOA(:,1) / XMD
END DO
!dust conversion : from part/part to kg/m3
DO JLAYER=NSV_DSTBEG,NSV_DSTEND
- ZP_SV(:,JLAYER) = ZP_SV(:,JLAYER) * XMOLARWEIGHT_DUST* ZP_RHOA(:) / XMD
+ ZP_SV(:,JLAYER) = ZP_SV(:,JLAYER) * XMOLARWEIGHT_DUST* ZP_RHOA(:,1) / XMD
END DO
!sea salt conversion : from part/part to kg/m3
DO JLAYER=NSV_SLTBEG,NSV_SLTEND
- ZP_SV(:,JLAYER) = ZP_SV(:,JLAYER) * XMOLARWEIGHT_SALT* ZP_RHOA(:) / XMD
+ ZP_SV(:,JLAYER) = ZP_SV(:,JLAYER) * XMOLARWEIGHT_SALT* ZP_RHOA(:,1) / XMD
END DO
!
!blowing snow conversion : from kg(snow)/kg(dry air) to kg(snow)/m3
DO JLAYER=NSV_SNWBEG,NSV_SNWEND
- ZP_SV(:,JLAYER) = ZP_SV(:,JLAYER) * ZP_RHOA(:)
+ ZP_SV(:,JLAYER) = ZP_SV(:,JLAYER) * ZP_RHOA(:,1)
END DO
IF(LBLOWSNOW) THEN ! Convert 2D blowing snow fields
! from kg(snow)/kg(dry air) to kg(snow)/m3
DO JLAYER=(NSV+1),KSV_SURF
- ZP_SV(:,JLAYER) = ZP_SV(:,JLAYER) * ZP_RHOA(:)
+ ZP_SV(:,JLAYER) = ZP_SV(:,JLAYER) * ZP_RHOA(:,1)
END DO
END IF
!
@@ -945,18 +1164,35 @@ ISHAPE_2 = (/KDIM1,KDIM2/)
!
! Arguments in call to surface:
!
-ZSFTH = XUNDEF
-ZSFTQ = XUNDEF
+ZSFTH = XUNDEF
+ZSFTH_SURF = XUNDEF
+ZSFTH_WALL = XUNDEF
+ZSFTH_ROOF = XUNDEF
+ZCD_ROOF = XUNDEF
+ZSFTQ = XUNDEF
+ZSFTQ_SURF = XUNDEF
+ZSFTQ_WALL = XUNDEF
+ZSFTQ_ROOF = XUNDEF
+!
IF (NSV>0) ZSFTS = XUNDEF
ZSFCO2 = XUNDEF
ZSFU = XUNDEF
ZSFV = XUNDEF
!
-ZSFTH (IIB:IIE,IJB:IJE) = RESHAPE(ZP_SFTH(:), ISHAPE_2)
-ZSFTQ (IIB:IIE,IJB:IJE) = RESHAPE(ZP_SFTQ(:), ISHAPE_2)
+ZSFTH (IIB:IIE,IJB:IJE) = RESHAPE(ZP_SFTH(:), ISHAPE_2)
+ZSFTH_SURF (IIB:IIE,IJB:IJE) = RESHAPE(ZP_SFTH_SURF(:), ISHAPE_2)
+ZSFTH_WALL (IIB:IIE,IJB:IJE) = RESHAPE(ZP_SFTH_WALL(:), ISHAPE_2)
+ZSFTH_ROOF (IIB:IIE,IJB:IJE) = RESHAPE(ZP_SFTH_ROOF(:), ISHAPE_2)
+ZCD_ROOF (IIB:IIE,IJB:IJE) = RESHAPE(ZP_CD_ROOF(:), ISHAPE_2)
+ZSFTQ (IIB:IIE,IJB:IJE) = RESHAPE(ZP_SFTQ(:), ISHAPE_2)
+ZSFTQ_SURF (IIB:IIE,IJB:IJE) = RESHAPE(ZP_SFTQ_SURF(:), ISHAPE_2)
+ZSFTQ_WALL (IIB:IIE,IJB:IJE) = RESHAPE(ZP_SFTQ_WALL(:), ISHAPE_2)
+ZSFTQ_ROOF (IIB:IIE,IJB:IJE) = RESHAPE(ZP_SFTQ_ROOF(:), ISHAPE_2)
+!
DO JLAYER=1,SIZE(PSFSV,3)
ZSFTS (IIB:IIE,IJB:IJE,JLAYER) = RESHAPE(ZP_SFTS(:,JLAYER), ISHAPE_2)
END DO
+!
ZSFCO2 (IIB:IIE,IJB:IJE) = RESHAPE(ZP_SFCO2(:), ISHAPE_2)
ZSFU (IIB:IIE,IJB:IJE) = RESHAPE(ZP_SFU(:), ISHAPE_2)
ZSFV (IIB:IIE,IJB:IJE) = RESHAPE(ZP_SFV(:), ISHAPE_2)
@@ -999,6 +1235,7 @@ DEALLOCATE(ZP_V )
DEALLOCATE(ZP_QA )
DEALLOCATE(ZP_TA )
DEALLOCATE(ZP_RHOA )
+DEALLOCATE(ZP_TKE )
DEALLOCATE(ZP_SV )
DEALLOCATE(ZP_CO2 )
DEALLOCATE(ZP_RAIN )
@@ -1009,9 +1246,16 @@ DEALLOCATE(ZP_SCA_SW )
DEALLOCATE(ZP_PS )
DEALLOCATE(ZP_PA )
DEALLOCATE(ZP_ZWS )
-
-DEALLOCATE(ZP_SFTQ )
-DEALLOCATE(ZP_SFTH )
+!
+DEALLOCATE(ZP_SFTQ )
+DEALLOCATE(ZP_SFTQ_SURF)
+DEALLOCATE(ZP_SFTQ_WALL)
+DEALLOCATE(ZP_SFTQ_ROOF)
+DEALLOCATE(ZP_SFTH )
+DEALLOCATE(ZP_SFTH_SURF)
+DEALLOCATE(ZP_SFTH_WALL)
+DEALLOCATE(ZP_SFTH_ROOF)
+DEALLOCATE(ZP_CD_ROOF)
DEALLOCATE(ZP_SFTS )
DEALLOCATE(ZP_SFCO2 )
DEALLOCATE(ZP_SFU )
diff --git a/src/MNH/ini_budget.f90 b/src/MNH/ini_budget.f90
index 3152cb6e5e17022fc128393ad806dff7b7a08fa6..f239b3719419e9d80cad23cac2de56b2faa148ab 100644
--- a/src/MNH/ini_budget.f90
+++ b/src/MNH/ini_budget.f90
@@ -106,7 +106,7 @@ end subroutine Budget_preallocate
OHORELAX_UVWTH,OHORELAX_RV,OHORELAX_RC,OHORELAX_RR, &
OHORELAX_RI,OHORELAX_RS, OHORELAX_RG, OHORELAX_RH,OHORELAX_TKE, &
OHORELAX_SV, OVE_RELAX, ove_relax_grd, OCHTRANS, &
- ONUDGING,ODRAGTREE,ODEPOTREE, OAERO_EOL, &
+ ONUDGING,ODRAGTREE,ODEPOTREE, ODRAGBLDG, OAERO_EOL, &
HRAD,HDCONV,HSCONV,HTURB,HTURBDIM,HCLOUD )
! #################################################################
!
@@ -208,6 +208,7 @@ end subroutine Budget_preallocate
! P. Wautelet 02/03/2021: budgets: add terms for blowing snow
! P. Wautelet 04/03/2021: budgets: add terms for drag due to buildings
! P. Wautelet 17/03/2021: choose source terms for budgets with character strings instead of multiple integer variables
+ ! R. Schoetter 12/2021 multi-level coupling between MesoNH and SURFEX
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -299,6 +300,7 @@ LOGICAL, INTENT(IN) :: OCHTRANS ! switch to activate convective
LOGICAL, INTENT(IN) :: ONUDGING ! switch to activate nudging
LOGICAL, INTENT(IN) :: ODRAGTREE ! switch to activate vegetation drag
LOGICAL, INTENT(IN) :: ODEPOTREE ! switch to activate droplet deposition on tree
+LOGICAL, INTENT(IN) :: ODRAGBLDG ! switch to activate building drag
LOGICAL, INTENT(IN) :: OAERO_EOL ! switch to activate wind turbine wake
CHARACTER (LEN=*), INTENT(IN) :: HRAD ! type of the radiation scheme
CHARACTER (LEN=*), INTENT(IN) :: HDCONV ! type of the deep convection scheme
@@ -1029,6 +1031,11 @@ if ( lbu_rth ) then
tzsource%lavailable = hdconv == 'KAFR' .OR. hsconv == 'KAFR'
call Budget_source_add( tbudgets(NBUDGET_TH), tzsource )
+ tzsource%cmnhname = 'DRAGB'
+ tzsource%clongname = 'heat released by buildings'
+ tzsource%lavailable = ldragbldg
+ call Budget_source_add( tbudgets(NBUDGET_TH), tzsource )
+
tzsource%cmnhname = 'VTURB'
tzsource%clongname = 'vertical turbulent diffusion'
tzsource%lavailable = hturb == 'TKEL'
@@ -1442,6 +1449,11 @@ if ( tbudgets(NBUDGET_RV)%lenabled ) then
tzsource%lavailable = hdconv == 'KAFR' .OR. hsconv == 'KAFR'
call Budget_source_add( tbudgets(NBUDGET_RV), tzsource )
+ tzsource%cmnhname = 'DRAGB'
+ tzsource%clongname = 'vapor released by buildings'
+ tzsource%lavailable = ldragbldg
+ call Budget_source_add( tbudgets(NBUDGET_RV), tzsource )
+
tzsource%cmnhname = 'VTURB'
tzsource%clongname = 'vertical turbulent diffusion'
tzsource%lavailable = hturb == 'TKEL'
diff --git a/src/MNH/ini_mean_field.f90 b/src/MNH/ini_mean_field.f90
index 32ddfbb8847adcfe2cbe1887368abee74b2370b9..2e1da698c6c6f86742ba1877fff7db554d992fe0 100644
--- a/src/MNH/ini_mean_field.f90
+++ b/src/MNH/ini_mean_field.f90
@@ -50,6 +50,7 @@ END MODULE MODI_INI_MEAN_FIELD
!! 10/2016 (C.Lac) Add max values
!! 02/2021 (T.Nagel) add passive scalar (XSVT) and UW wind component
!! 05/2021 (PA.Joulin) add wind turbine variables
+!! 12/2021 (R. Schoetter) adds humidity and other mean diagnostics
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -79,6 +80,20 @@ IF (CTURB /= 'NONE') XTKEM_MEAN = 0.0
XPABSM_MEAN = 0.0
XSVT_MEAN = 0.0
!
+XQ_MEAN = 0.0
+XRH_W_MEAN = 0.0
+XRH_I_MEAN = 0.0
+XRH_P_MEAN = 0.0
+XRH_W_MAXCOL_MEAN = 0.0
+XRH_I_MAXCOL_MEAN = 0.0
+XRH_P_MAXCOL_MEAN = 0.0
+!
+XWIFF_MEAN = 0.0
+XWIDD_MEAN = 0.0
+!
+XWIFF_MAX = 0.0
+XWIDD_MAX = 0.0
+!
XU2_MEAN = 0.0
XV2_MEAN = 0.0
XW2_MEAN = 0.0
diff --git a/src/MNH/ini_modeln.f90 b/src/MNH/ini_modeln.f90
index d5a6ffa05b054a5ebe8d8ef848ec4cf44b3cef0e..1e03ebd940d87d95adf7c316c8fa5946bbefaa7f 100644
--- a/src/MNH/ini_modeln.f90
+++ b/src/MNH/ini_modeln.f90
@@ -293,6 +293,8 @@ END MODULE MODI_INI_MODEL_n
! F. Auguste 02/2021: add IBM
! T.Nigel 02/2021: add turbulence recycling
! J.L.Redelsperger 06/2011: OCEAN case
+ ! R. Schoetter 12/2021 multi-level coupling between MesoNH and SURFEX
+!! 12/2021 (R. Schoetter) adds humidity and other mean diagnostics
!---------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -333,6 +335,7 @@ USE MODD_DIAG_FLAG, only: LCHEMDIAG, CSPEC_BU_DIAG
USE MODD_DIM_n
USE MODD_DRAG_n
USE MODD_DRAGTREE_n
+USE MODD_DRAGBLDG_n
USE MODD_DUST
use MODD_DUST_OPT_LKT, only: NMAX_RADIUS_LKT_DUST=>NMAX_RADIUS_LKT, NMAX_SIGMA_LKT_DUST=>NMAX_SIGMA_LKT, &
NMAX_WVL_SW_DUST=>NMAX_WVL_SW, &
@@ -772,6 +775,17 @@ IF (LMEAN_FIELD) THEN
ALLOCATE(XTKEM_MEAN(0,0,0))
END IF
ALLOCATE(XPABSM_MEAN(IIU,IJU,IKU)) ; XPABSM_MEAN = 0.0
+ ALLOCATE(XQ_MEAN(IIU,IJU,IKU)) ; XQ_MEAN = 0.0
+ ALLOCATE(XRH_W_MEAN(IIU,IJU,IKU)) ; XRH_W_MEAN = 0.0
+ ALLOCATE(XRH_I_MEAN(IIU,IJU,IKU)) ; XRH_I_MEAN = 0.0
+ ALLOCATE(XRH_P_MEAN(IIU,IJU,IKU)) ; XRH_P_MEAN = 0.0
+ ALLOCATE(XRH_W_MAXCOL_MEAN(IIU,IJU)) ; XRH_W_MAXCOL_MEAN = 0.0
+ ALLOCATE(XRH_I_MAXCOL_MEAN(IIU,IJU)) ; XRH_I_MAXCOL_MEAN = 0.0
+ ALLOCATE(XRH_P_MAXCOL_MEAN(IIU,IJU)) ; XRH_P_MAXCOL_MEAN = 0.0
+ ALLOCATE(XWIFF_MEAN(IIU,IJU,IKU)) ; XWIFF_MEAN = 0.0
+ ALLOCATE(XWIDD_MEAN(IIU,IJU,IKU)) ; XWIDD_MEAN = 0.0
+ ALLOCATE(XWIFF_MAX (IIU,IJU,IKU)) ; XWIFF_MAX = 0.0
+ ALLOCATE(XWIDD_MAX (IIU,IJU,IKU)) ; XWIDD_MAX = 0.0
!
ALLOCATE(XU2_MEAN(IIU,IJU,IKU)) ; XU2_MEAN = 0.0
ALLOCATE(XV2_MEAN(IIU,IJU,IKU)) ; XV2_MEAN = 0.0
@@ -1808,7 +1822,7 @@ IF ( CBUTYPE /= "NONE" .AND. NBUMOD == KMI ) THEN
LHORELAX_UVWTH,LHORELAX_RV, LHORELAX_RC,LHORELAX_RR, &
LHORELAX_RI,LHORELAX_RS,LHORELAX_RG, LHORELAX_RH,LHORELAX_TKE, &
LHORELAX_SV, LVE_RELAX, LVE_RELAX_GRD, &
- LCHTRANS,LNUDGING,LDRAGTREE,LDEPOTREE,LMAIN_EOL, &
+ LCHTRANS,LNUDGING,LDRAGTREE,LDEPOTREE,LDRAGBLDG,LMAIN_EOL, &
CRAD,CDCONV,CSCONV,CTURB,CTURBDIM,CCLOUD )
END IF
!
diff --git a/src/MNH/mean_field.f90 b/src/MNH/mean_field.f90
index 048eb689dfa04f798b16e555bb6b220e896b3484..85c5c54989488114260799ce84ef6e27e71553cf 100644
--- a/src/MNH/mean_field.f90
+++ b/src/MNH/mean_field.f90
@@ -10,11 +10,12 @@
!
INTERFACE
- SUBROUTINE MEAN_FIELD( PUT, PVT, PWT, PTHT, PTKET, PPABST, PSVT )
+ SUBROUTINE MEAN_FIELD( PUT, PVT, PWT, PTHT, PTKET, PPABST, PRT, PSVT )
REAL, DIMENSION(:,:,:), INTENT(IN) :: PUT, PVT, PWT ! variables
REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHT, PTKET ! variables
REAL, DIMENSION(:,:,:), INTENT(IN) :: PPABST ! variables
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRT ! variables
REAL, DIMENSION(:,:,:), INTENT(IN) :: PSVT ! Passive scalar variables
END SUBROUTINE MEAN_FIELD
@@ -23,9 +24,9 @@ END INTERFACE
END MODULE MODI_MEAN_FIELD
!
-! #################################################################
- SUBROUTINE MEAN_FIELD( PUT, PVT, PWT, PTHT, PTKET, PPABST, PSVT )
-! #################################################################
+! ######################################################################
+ SUBROUTINE MEAN_FIELD( PUT, PVT, PWT, PTHT, PTKET, PPABST, PRT, PSVT )
+! ######################################################################
!
!!**** *MEAN_FIELD * -
!!
@@ -48,6 +49,7 @@ END MODULE MODI_MEAN_FIELD
!! Original 07/2009
!! (C.Lac) 09/2016 Max values
!! (PA.Joulin) 12/2020 Wind turbine variables
+!! 12/2021 (R. Schoetter) adds humidity and other mean diagnostics
!!---------------------------------------------------------------
!
!
@@ -60,6 +62,8 @@ USE MODD_PARAM_n
USE MODD_MEAN_FIELD
USE MODD_CST
USE MODD_PASPOL
+USE MODE_THERMO
+USE MODI_SHUMAN
!
USE MODD_EOL_MAIN, ONLY: LMAIN_EOL, CMETH_EOL, NMODEL_EOL
USE MODD_EOL_SHARED_IO, ONLY: XTHRUT, XTORQT, XPOWT
@@ -75,6 +79,7 @@ IMPLICIT NONE
REAL, DIMENSION(:,:,:), INTENT(IN) :: PUT, PVT, PWT ! variables
REAL, DIMENSION(:,:,:), INTENT(IN) :: PTHT, PTKET ! variables
REAL, DIMENSION(:,:,:), INTENT(IN) :: PPABST ! variables
+REAL, DIMENSION(:,:,:), INTENT(IN) :: PRT ! variables
REAL, DIMENSION(:,:,:), INTENT(IN) :: PSVT
!
@@ -83,6 +88,20 @@ REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZTEMPT
INTEGER :: IIU,IJU,IKU,IIB,IJB,IKB,IIE,IJE,IKE ! Arrays bounds
INTEGER :: JI,JJ,JK ! Loop indexes
!
+REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZRT
+REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZQSAT_W
+REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZQSAT_I
+REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZQACT
+REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZRH_W
+REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZRH_I
+REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZRH_P
+REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZFRAC
+REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZAUX_WIFF
+REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZAUX_WIDD
+REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2)) :: ZRH_W_MAXCOL
+REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2)) :: ZRH_I_MAXCOL
+REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2)) :: ZRH_P_MAXCOL
+!
INTEGER :: IMI !Current model index
!
!
@@ -102,6 +121,44 @@ IKE=IKU-JPVEXT
!* 1. MEAN
!
ZTEMPT = PTHT*(((PPABST)/XP00)**(XRD/XCPD))
+!
+!
+! Calculation of saturation specific humidity over water/ice
+!
+ZQSAT_W = QSAT (ZTEMPT, PPABST)
+ZQSAT_I = QSATI(ZTEMPT, PPABST)
+!
+! Conversion mixing ratio -> specfic humidity
+!
+ZRT(:,:,:) = -9999.0
+WHERE (PRT(:,:,:).LT.1.0E-6)
+ ZRT(:,:,:) = 1.0E-6
+ELSEWHERE
+ ZRT(:,:,:) = PRT(:,:,:)
+ENDWHERE
+!
+ZQACT(:,:,:) = 1.0 / ( 1.0 + 1.0/ZRT(:,:,:) )
+!
+! Calculation of relative humidity with respect to water/ice
+!
+ZRH_W(:,:,:) = 100.0*ZQACT(:,:,:)/ZQSAT_W(:,:,:)
+ZRH_I(:,:,:) = 100.0*ZQACT(:,:,:)/ZQSAT_I(:,:,:)
+!
+! Fractional partitioning between liquid and solid cloud water
+! as assumed in condensations
+!
+ZFRAC(:,:,:) = ( XTT - ZTEMPT(:,:,:) ) / 20.
+ZFRAC(:,:,:) = MAX( 0., MIN(1., ZFRAC(:,:,:) ) )
+!
+! Calculation of weighted average between water and ice value
+!
+ZRH_P(:,:,:) = ZFRAC(:,:,:) * ZRH_I(:,:,:) + (1.0-ZFRAC(:,:,:)) * ZRH_W(:,:,:)
+!
+! Calculation of the column maximum of relative humidity
+!
+ZRH_W_MAXCOL(:,:) = MAXVAL(ZRH_W(:,:,:),DIM=3)
+ZRH_I_MAXCOL(:,:) = MAXVAL(ZRH_I(:,:,:),DIM=3)
+ZRH_P_MAXCOL(:,:) = MAXVAL(ZRH_P(:,:,:),DIM=3)
!
XUM_MEAN = PUT + XUM_MEAN
XVM_MEAN = PVT + XVM_MEAN
@@ -111,6 +168,13 @@ IKE=IKU-JPVEXT
IF (LPASPOL) XSVT_MEAN = PSVT + XSVT_MEAN
IF (CTURB/='NONE') XTKEM_MEAN = PTKET + XTKEM_MEAN
XPABSM_MEAN = PPABST + XPABSM_MEAN
+ XQ_MEAN = XQ_MEAN + ZQACT
+ XRH_W_MEAN = XRH_W_MEAN + ZRH_W
+ XRH_I_MEAN = XRH_I_MEAN + ZRH_I
+ XRH_P_MEAN = XRH_P_MEAN + ZRH_P
+ XRH_W_MAXCOL_MEAN = XRH_W_MAXCOL_MEAN + ZRH_W_MAXCOL
+ XRH_I_MAXCOL_MEAN = XRH_I_MAXCOL_MEAN + ZRH_I_MAXCOL
+ XRH_P_MAXCOL_MEAN = XRH_P_MAXCOL_MEAN + ZRH_P_MAXCOL
!
XU2_MEAN = PUT**2 + XU2_MEAN
XV2_MEAN = PVT**2 + XV2_MEAN
@@ -142,19 +206,36 @@ IKE=IKU-JPVEXT
!
!* 2. MAX
!
+ !
+ ! Calculation of horizontal wind speed for maximum wind speed diagnostics
+ !
+ ZAUX_WIFF(:,:,:) = SQRT(MXF(PUT(:,:,:))**2 + MYF(PVT(:,:,:))**2)
+ ZAUX_WIDD(:,:,:) = 180.0 + (90.0 - 180.0*ATAN2(MYF(PVT(:,:,:)),MXF(PUT(:,:,:)))/XPI)
+ !
+ WHERE (ZAUX_WIDD(:,:,:).GT.360.0) ZAUX_WIDD(:,:,:) = ZAUX_WIDD(:,:,:) - 360.0
+ !
+ ! Get maximum diagnostics
+ !
DO JK=IKB,IKE
- DO JJ=IJB,IJE
- DO JI=IIB,IIE
- XUM_MAX(JI,JJ,JK) = MAX(XUM_MAX(JI,JJ,JK),PUT(JI,JJ,JK))
- XVM_MAX(JI,JJ,JK) = MAX(XVM_MAX(JI,JJ,JK),PVT(JI,JJ,JK))
- XWM_MAX(JI,JJ,JK) = MAX(XWM_MAX(JI,JJ,JK),PWT(JI,JJ,JK))
- XTHM_MAX(JI,JJ,JK) = MAX(XTHM_MAX(JI,JJ,JK),PTHT(JI,JJ,JK))
- XTEMPM_MAX(JI,JJ,JK) = MAX(XTEMPM_MAX(JI,JJ,JK),ZTEMPT(JI,JJ,JK))
- IF (CTURB/='NONE') XTKEM_MAX(JI,JJ,JK) = &
- MAX(XTKEM_MAX(JI,JJ,JK),PTKET(JI,JJ,JK))
- XPABSM_MAX(JI,JJ,JK) = MAX(XPABSM_MAX(JI,JJ,JK),PPABST(JI,JJ,JK))
- END DO
- END DO
- END DO
+ DO JJ=IJB,IJE
+ DO JI=IIB,IIE
+ !
+ XUM_MAX(JI,JJ,JK) = MAX(XUM_MAX(JI,JJ,JK),PUT(JI,JJ,JK))
+ XVM_MAX(JI,JJ,JK) = MAX(XVM_MAX(JI,JJ,JK),PVT(JI,JJ,JK))
+ XWM_MAX(JI,JJ,JK) = MAX(XWM_MAX(JI,JJ,JK),PWT(JI,JJ,JK))
+ XTHM_MAX(JI,JJ,JK) = MAX(XTHM_MAX(JI,JJ,JK),PTHT(JI,JJ,JK))
+ XTEMPM_MAX(JI,JJ,JK) = MAX(XTEMPM_MAX(JI,JJ,JK),ZTEMPT(JI,JJ,JK))
+ IF (CTURB/='NONE') XTKEM_MAX(JI,JJ,JK) = &
+ MAX(XTKEM_MAX(JI,JJ,JK),PTKET(JI,JJ,JK))
+ XPABSM_MAX(JI,JJ,JK) = MAX(XPABSM_MAX(JI,JJ,JK),PPABST(JI,JJ,JK))
+ !
+ IF (ZAUX_WIFF(JI,JJ,JK).GE.XWIFF_MAX(JI,JJ,JK)) THEN
+ XWIFF_MAX(JI,JJ,JK) = ZAUX_WIFF(JI,JJ,JK)
+ XWIDD_MAX(JI,JJ,JK) = ZAUX_WIDD(JI,JJ,JK)
+ ENDIF
+ !
+ ENDDO
+ ENDDO
+ ENDDO
!
END SUBROUTINE MEAN_FIELD
diff --git a/src/MNH/mnhget_surf_paramn.f90 b/src/MNH/mnhget_surf_paramn.f90
index 72f042df0482817b41343f423f4fd683ca9e09b9..31d8e8820126007aaba14a56114ed1c2807c70b8 100644
--- a/src/MNH/mnhget_surf_paramn.f90
+++ b/src/MNH/mnhget_surf_paramn.f90
@@ -7,10 +7,11 @@
MODULE MODI_MNHGET_SURF_PARAM_n
! #######################
INTERFACE
- SUBROUTINE MNHGET_SURF_PARAM_n(PCOVER,PSEA,KCOVER,PRN,PH,PLE,PLEI,PGFLUX, &
- PT2M,PQ2M,PHU2M,PZON10M,PMER10M,PZS,PTOWN,&
- PBARE, PLAI_TREE, PH_TREE, PWALL_O_HOR, &
- PBUILD_HEIGHT,PNATURE )
+ SUBROUTINE MNHGET_SURF_PARAM_n(PCOVER,PSEA,KCOVER,PRN,PH,PLE,PLEI,PGFLUX, &
+ PT2M,PQ2M,PHU2M,PZON10M,PMER10M,PZS,PTOWN, &
+ PBARE,PLAI_TREE,PH_TREE,PWALL_O_HOR, &
+ PBUILD_HEIGHT,PNATURE,PLAI_HVEG,PH_URBTREE, &
+ PHTRUNK_HVEG,PFRAC_HVEG)
!
REAL, DIMENSION(:,:,:), INTENT(OUT), OPTIONAL :: PCOVER ! cover types
REAL, DIMENSION(:,:), INTENT(OUT), OPTIONAL :: PSEA ! sea fraction
@@ -22,6 +23,10 @@ REAL, DIMENSION(:,:), INTENT(OUT), OPTIONAL :: PLAI_TREE ! Tree leaf are
REAL, DIMENSION(:,:), INTENT(OUT), OPTIONAL :: PH_TREE ! Tree height [m]
REAL, DIMENSION(:,:), INTENT(OUT), OPTIONAL :: PWALL_O_HOR ! Facade area density [m^2(fac.)/m^2(town)]
REAL, DIMENSION(:,:), INTENT(OUT), OPTIONAL :: PBUILD_HEIGHT ! Building height [m]
+REAL, DIMENSION(:,:), INTENT(OUT), OPTIONAL :: PLAI_HVEG ! LAI of urban vegetation
+REAL, DIMENSION(:,:), INTENT(OUT), OPTIONAL :: PH_URBTREE ! Height of urban vegetation
+REAL, DIMENSION(:,:), INTENT(OUT), OPTIONAL :: PHTRUNK_HVEG ! Trunk height of urban vegetation
+REAL, DIMENSION(:,:), INTENT(OUT), OPTIONAL :: PFRAC_HVEG ! Fraction of high vegetation
REAL, DIMENSION(:), INTENT(INOUT), OPTIONAL :: PRN ! Net radiation at surface (W/m2)
REAL, DIMENSION(:), INTENT(INOUT), OPTIONAL :: PH ! Sensible heat flux (W/m2)
REAL, DIMENSION(:), INTENT(INOUT), OPTIONAL :: PLE ! Total Latent heat flux (W/m2)
@@ -40,10 +45,11 @@ END INTERFACE
END MODULE MODI_MNHGET_SURF_PARAM_n
!
! ########################################
- SUBROUTINE MNHGET_SURF_PARAM_n(PCOVER,PSEA,KCOVER,PRN,PH,PLE,PLEI,PGFLUX, &
- PT2M,PQ2M,PHU2M,PZON10M,PMER10M,PZS,PTOWN,&
- PBARE, PLAI_TREE, PH_TREE, PWALL_O_HOR, &
- PBUILD_HEIGHT,PNATURE )
+ SUBROUTINE MNHGET_SURF_PARAM_n(PCOVER,PSEA,KCOVER,PRN,PH,PLE,PLEI,PGFLUX, &
+ PT2M,PQ2M,PHU2M,PZON10M,PMER10M,PZS,PTOWN, &
+ PBARE, PLAI_TREE, PH_TREE, PWALL_O_HOR, &
+ PBUILD_HEIGHT,PNATURE,PLAI_HVEG,PH_URBTREE, &
+ PHTRUNK_HVEG,PFRAC_HVEG )
! ########################################
!
!!**** *MNHGET_SURF_PARAM_n* - gets some surface fields on MESONH grid
@@ -81,6 +87,7 @@ END MODULE MODI_MNHGET_SURF_PARAM_n
!! 01/2018 (G.Delautier) SURFEX 8.1
! C. Lac 11/2019: correction in the drag formula and application to building in addition to tree
! P. Wautelet 11/03/2020: bugfix: add present checks before working on optional arrays
+ ! R. Schoetter 12/2021 multi-level coupling between MesoNH and SURFEX
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -115,6 +122,10 @@ REAL, DIMENSION(:,:), INTENT(OUT), OPTIONAL :: PLAI_TREE !
REAL, DIMENSION(:,:), INTENT(OUT), OPTIONAL :: PH_TREE !
REAL, DIMENSION(:,:), INTENT(OUT), OPTIONAL :: PWALL_O_HOR ! Facade area density [m^2(fac.)/m^2(town)]
REAL, DIMENSION(:,:), INTENT(OUT), OPTIONAL :: PBUILD_HEIGHT ! Building height [m]
+REAL, DIMENSION(:,:), INTENT(OUT), OPTIONAL :: PLAI_HVEG ! LAI of urban vegetation
+REAL, DIMENSION(:,:), INTENT(OUT), OPTIONAL :: PH_URBTREE ! Height of urban vegetation
+REAL, DIMENSION(:,:), INTENT(OUT), OPTIONAL :: PHTRUNK_HVEG ! Trunk height of urban vegetation
+REAL, DIMENSION(:,:), INTENT(OUT), OPTIONAL :: PFRAC_HVEG ! Fraction of high vegetation
REAL, DIMENSION(:), INTENT(INOUT), OPTIONAL :: PRN ! Net radiation at surface (W/m2)
REAL, DIMENSION(:), INTENT(INOUT), OPTIONAL :: PH ! Sensible heat flux (W/m2)
REAL, DIMENSION(:), INTENT(INOUT), OPTIONAL :: PLE ! Total Latent heat flux (W/m2)
@@ -141,13 +152,17 @@ INTEGER :: ILU ! total number of physical points (including halo)
INTEGER :: ILM ! total number of physical points
!
REAL, DIMENSION(:), ALLOCATABLE :: ZCOVER ! cover field
-REAL, DIMENSION(:), ALLOCATABLE :: ZSEA ! sea fraction
-REAL, DIMENSION(:), ALLOCATABLE :: ZWATER ! lake fraction
-REAL, DIMENSION(:), ALLOCATABLE :: ZNATURE! nature fraction
-REAL, DIMENSION(:), ALLOCATABLE :: ZTOWN ! town fraction
-REAL, DIMENSION(:), ALLOCATABLE :: ZVH
-REAL, DIMENSION(:), ALLOCATABLE :: ZLAI
-REAL, DIMENSION(:), ALLOCATABLE :: ZWALL_O_HOR ! Facade surface density [m^2(fac.)/m^2(town)]
+REAL, DIMENSION(:), ALLOCATABLE :: ZSEA ! sea fraction
+REAL, DIMENSION(:), ALLOCATABLE :: ZWATER ! lake fraction
+REAL, DIMENSION(:), ALLOCATABLE :: ZNATURE! nature fraction
+REAL, DIMENSION(:), ALLOCATABLE :: ZTOWN ! town fraction
+REAL, DIMENSION(:), ALLOCATABLE :: ZVH
+REAL, DIMENSION(:), ALLOCATABLE :: ZLAI
+REAL, DIMENSION(:), ALLOCATABLE :: ZWALL_O_HOR ! Facade surface density [m^2(fac.)/m^2(town)]
+REAL, DIMENSION(:), ALLOCATABLE :: ZLAI_HVEG ! LAI of urban vegetation
+REAL, DIMENSION(:), ALLOCATABLE :: ZH_URBTREE ! Height of urban vegetation
+REAL, DIMENSION(:), ALLOCATABLE :: ZHTRUNK_HVEG ! Trunk height of urban vegetation
+REAL, DIMENSION(:), ALLOCATABLE :: ZFRAC_HVEG ! Fraction of high vegetation
REAL, DIMENSION(:), ALLOCATABLE :: ZBUILD_HEIGHT ! Building height [m]
REAL, DIMENSION(:), ALLOCATABLE :: ZBARE ! bare soil fraction
REAL, DIMENSION(:), ALLOCATABLE :: ZZS ! orography
@@ -182,9 +197,12 @@ IF (PRESENT(PCOVER)) THEN
DEALLOCATE(ZCOVER)
END IF
!
-IF (PRESENT(PSEA) .OR. PRESENT(PTOWN) .OR. PRESENT(PNATURE) .OR. &
- PRESENT(PBARE) .OR. PRESENT(PLAI_TREE) .OR. PRESENT(PH_TREE) .OR. &
- PRESENT(PWALL_O_HOR) .OR. PRESENT(PBUILD_HEIGHT) ) THEN
+IF ( PRESENT(PSEA) .OR. PRESENT(PTOWN) .OR. PRESENT(PNATURE) .OR. &
+ PRESENT(PBARE) .OR. PRESENT(PLAI_TREE) .OR. PRESENT(PH_TREE) .OR. &
+ PRESENT(PWALL_O_HOR) .OR. PRESENT(PBUILD_HEIGHT) .OR. &
+ PRESENT(PLAI_HVEG) .OR. PRESENT(PH_URBTREE) .OR. &
+ PRESENT(PHTRUNK_HVEG) .OR. PRESENT(PFRAC_HVEG) ) THEN
+ !
ALLOCATE(ZSEA ( ILU ))
ALLOCATE(ZWATER ( ILU ))
ALLOCATE(ZNATURE( ILU ))
@@ -203,10 +221,10 @@ END IF
!
IF (PRESENT(PBARE)) THEN
ALLOCATE(ZBARE ( ILU ))
- CALL GET_SURF_VAR_n(YSURF_CUR%FM,YSURF_CUR%IM,YSURF_CUR%SM,YSURF_CUR%TM, &
- YSURF_CUR%GDM,YSURF_CUR%WM,YSURF_CUR%DUO,YSURF_CUR%DU,&
- YSURF_CUR%UG, YSURF_CUR%U,YSURF_CUR%USS,&
- 'MESONH', ILU, 1, PNATURE=ZNATURE, PBARE=ZBARE)
+ CALL GET_SURF_VAR_n(YSURF_CUR%FM, YSURF_CUR%IM, YSURF_CUR%SM, YSURF_CUR%TM, &
+ YSURF_CUR%GDM, YSURF_CUR%WM, YSURF_CUR%DUO, YSURF_CUR%DU, &
+ YSURF_CUR%UG, YSURF_CUR%U, YSURF_CUR%USS, 'MESONH', ILU, &
+ 1, PNATURE=ZNATURE, PBARE=ZBARE )
CALL REMOVE_HALO(ZBARE,PBARE)
DEALLOCATE(ZBARE)
END IF
@@ -254,30 +272,55 @@ IF (PRESENT(PH_TREE) .OR.PRESENT(PLAI_TREE)) THEN
PLAI_TREE(:,:) = XUNDEF
ALLOCATE(ZVH ( ILU ))
ALLOCATE(ZLAI ( ILU ))
- CALL GET_SURF_VAR_n(YSURF_CUR%FM,YSURF_CUR%IM,YSURF_CUR%SM,YSURF_CUR%TM, &
- YSURF_CUR%GDM,YSURF_CUR%WM,YSURF_CUR%DUO,YSURF_CUR%DU,&
- YSURF_CUR%UG,YSURF_CUR%U,YSURF_CUR%USS,&
- 'MESONH',ILU,1,PNATURE=ZNATURE,PLAI_TREE=ZLAI,PH_TREE=ZVH)
+ CALL GET_SURF_VAR_n(YSURF_CUR%FM,YSURF_CUR%IM,YSURF_CUR%SM,YSURF_CUR%TM, &
+ YSURF_CUR%GDM, YSURF_CUR%WM,YSURF_CUR%DUO,YSURF_CUR%DU, &
+ YSURF_CUR%UG, YSURF_CUR%U,YSURF_CUR%USS, 'MESONH',ILU,1,&
+ PNATURE=ZNATURE, PLAI_TREE=ZLAI, PH_TREE=ZVH )
IF ( PRESENT( PLAI_TREE ) ) CALL REMOVE_HALO(ZLAI,PLAI_TREE)
IF ( PRESENT( PH_TREE ) ) CALL REMOVE_HALO(ZVH,PH_TREE)
DEALLOCATE(ZVH)
DEALLOCATE(ZLAI)
END IF
!
-IF (PRESENT(PWALL_O_HOR) .OR. PRESENT(PBUILD_HEIGHT)) THEN
- IF ( PRESENT ( PBUILD_HEIGHT ) ) PBUILD_HEIGHT(:,:) = XUNDEF
- IF ( PRESENT ( PWALL_O_HOR ) ) PWALL_O_HOR(:,:) = XUNDEF
- ALLOCATE(ZBUILD_HEIGHT ( ILU ))
- ALLOCATE(ZWALL_O_HOR ( ILU ))
- CALL GET_SURF_VAR_n(YSURF_CUR%FM,YSURF_CUR%IM,YSURF_CUR%SM,YSURF_CUR%TM, &
- YSURF_CUR%GDM,YSURF_CUR%WM,YSURF_CUR%DUO,YSURF_CUR%DU,&
- YSURF_CUR%UG,YSURF_CUR%U,YSURF_CUR%USS,&
- 'MESONH',ILU,1,PTOWN=ZTOWN, &
- PWALL_O_HOR=ZWALL_O_HOR,PBUILD_HEIGHT=ZBUILD_HEIGHT )
- IF ( PRESENT ( PBUILD_HEIGHT ) ) CALL REMOVE_HALO(ZBUILD_HEIGHT,PBUILD_HEIGHT)
- IF ( PRESENT ( PWALL_O_HOR ) ) CALL REMOVE_HALO(ZWALL_O_HOR,PWALL_O_HOR)
- DEALLOCATE(ZBUILD_HEIGHT)
- DEALLOCATE(ZWALL_O_HOR)
+IF ( PRESENT(PWALL_O_HOR) .OR. PRESENT(PBUILD_HEIGHT) .OR. PRESENT(PLAI_HVEG) .OR. &
+ PRESENT(PH_URBTREE) .OR. PRESENT(PHTRUNK_HVEG) .OR. PRESENT(PFRAC_HVEG) ) THEN
+ !
+ IF (PRESENT(PBUILD_HEIGHT)) PBUILD_HEIGHT (:,:) = XUNDEF
+ IF (PRESENT(PWALL_O_HOR)) PWALL_O_HOR (:,:) = XUNDEF
+ IF (PRESENT(PLAI_HVEG)) PLAI_HVEG (:,:) = XUNDEF
+ IF (PRESENT(PH_URBTREE)) PH_URBTREE (:,:) = XUNDEF
+ IF (PRESENT(PHTRUNK_HVEG)) PHTRUNK_HVEG (:,:) = XUNDEF
+ IF (PRESENT(PFRAC_HVEG)) PFRAC_HVEG (:,:) = XUNDEF
+ !
+ ALLOCATE(ZBUILD_HEIGHT ( ILU ))
+ ALLOCATE(ZWALL_O_HOR ( ILU ))
+ ALLOCATE(ZLAI_HVEG ( ILU ))
+ ALLOCATE(ZH_URBTREE ( ILU ))
+ ALLOCATE(ZHTRUNK_HVEG ( ILU ))
+ ALLOCATE(ZFRAC_HVEG ( ILU ))
+ !
+ CALL GET_SURF_VAR_n(YSURF_CUR%FM,YSURF_CUR%IM,YSURF_CUR%SM,YSURF_CUR%TM, &
+ YSURF_CUR%GDM, YSURF_CUR%WM,YSURF_CUR%DUO,YSURF_CUR%DU,&
+ YSURF_CUR%UG, YSURF_CUR%U,YSURF_CUR%USS, 'MESONH',ILU, &
+ 1,PTOWN=ZTOWN, PWALL_O_HOR=ZWALL_O_HOR, &
+ PBUILD_HEIGHT=ZBUILD_HEIGHT, PLAI_HVEG=ZLAI_HVEG, &
+ PH_URBTREE=ZH_URBTREE, PHTRUNK_HVEG=ZHTRUNK_HVEG, &
+ PFRAC_HVEG=ZFRAC_HVEG )
+ !
+ IF (PRESENT(PBUILD_HEIGHT)) CALL REMOVE_HALO(ZBUILD_HEIGHT,PBUILD_HEIGHT)
+ IF (PRESENT(PWALL_O_HOR)) CALL REMOVE_HALO(ZWALL_O_HOR,PWALL_O_HOR)
+ IF (PRESENT(PLAI_HVEG)) CALL REMOVE_HALO(ZLAI_HVEG,PLAI_HVEG)
+ IF (PRESENT(PH_URBTREE)) CALL REMOVE_HALO(ZH_URBTREE,PH_URBTREE)
+ IF (PRESENT(PHTRUNK_HVEG)) CALL REMOVE_HALO(ZHTRUNK_HVEG,PHTRUNK_HVEG)
+ IF (PRESENT(PFRAC_HVEG)) CALL REMOVE_HALO(ZFRAC_HVEG,PFRAC_HVEG)
+ !
+ DEALLOCATE(ZBUILD_HEIGHT)
+ DEALLOCATE(ZWALL_O_HOR)
+ DEALLOCATE(ZLAI_HVEG)
+ DEALLOCATE(ZH_URBTREE)
+ DEALLOCATE(ZHTRUNK_HVEG)
+ DEALLOCATE(ZFRAC_HVEG)
+ !
END IF
!
IF (ALLOCATED(ZSEA)) THEN
diff --git a/src/MNH/modd_budget.f90 b/src/MNH/modd_budget.f90
index 7442dfd3e0184598cd4ee7687456a5e31b900e5b..2cfdc2a04093c7d0ad22de76f196cdfd93e030df 100644
--- a/src/MNH/modd_budget.f90
+++ b/src/MNH/modd_budget.f90
@@ -48,6 +48,7 @@
! P. Wautelet 03/03/2021: add tbudiachrometadata type (useful to pass more information to Write_diachro)
! P. Wautelet 17/03/2021: choose source terms for budgets with character strings instead of multiple integer variables
! P. Wautelet 30/03/2021: budgets: cartesian subdomain limits are defined in the physical domain
+ ! R. Schoetter 12/2021 multi-level coupling between MesoNH and SURFEX
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
diff --git a/src/MNH/modd_dragbldgn.f90 b/src/MNH/modd_dragbldgn.f90
index f89fc0d09d34143191714cca3d89bbf9a55ab580..b560a689d3e1cbc76912b6babbc517522183cc5c 100644
--- a/src/MNH/modd_dragbldgn.f90
+++ b/src/MNH/modd_dragbldgn.f90
@@ -21,6 +21,7 @@
!! MODIFICATIONS
!! -------------
!! Original 09/2019
+ ! R. Schoetter 12/2021 multi-level coupling between MesoNH and SURFEX
!-----------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -35,11 +36,17 @@ TYPE DRAGBLDG_t
LOGICAL :: LDRAGBLDG ! flag used to take into account building drag in
! ! the atmospheric model instead of SURFEX.
!
+ LOGICAL :: LFLUXBLDG ! Flag used to take into account heat and moisture fluxes in
+ ! ! the atmospheric model instead of SURFEX.
+ LOGICAL :: LDRAGURBVEG ! Flag used to take into account drag of urban vegetation in
+ ! ! the atmospheric model instead of SURFEX.
END TYPE DRAGBLDG_t
!
TYPE(DRAGBLDG_t), DIMENSION(JPMODELMAX), TARGET, SAVE :: DRAGBLDG_MODEL
!
LOGICAL, POINTER :: LDRAGBLDG=>NULL()
+LOGICAL, POINTER :: LFLUXBLDG=>NULL()
+LOGICAL, POINTER :: LDRAGURBVEG=>NULL()
!
CONTAINS
!
@@ -48,6 +55,8 @@ SUBROUTINE DRAGBLDG_GOTO_MODEL(KFROM, KTO)
INTEGER, INTENT(IN) :: KFROM, KTO
!
LDRAGBLDG=>DRAGBLDG_MODEL(KTO)%LDRAGBLDG
+ LFLUXBLDG=>DRAGBLDG_MODEL(KTO)%LFLUXBLDG
+ LDRAGURBVEG=>DRAGBLDG_MODEL(KTO)%LDRAGURBVEG
!
END SUBROUTINE DRAGBLDG_GOTO_MODEL
!
diff --git a/src/MNH/modd_mean_fieldn.f90 b/src/MNH/modd_mean_fieldn.f90
index 38572bc0bb2cb8f50ea1bc4cd4352100af779c84..933368deab7528d4a451e8c9e39a9f4e86989fb9 100644
--- a/src/MNH/modd_mean_fieldn.f90
+++ b/src/MNH/modd_mean_fieldn.f90
@@ -31,6 +31,7 @@
!! -------------
!! Original 01/07/11
!! 10/2016 (C.Lac) Add max values
+!! 12/2021 (R. Schoetter) adds humidity and other mean diagnostics
!!
!-------------------------------------------------------------------------------
!
@@ -46,6 +47,17 @@ TYPE MEAN_FIELD_t
REAL, DIMENSION(:,:,:), POINTER :: XTEMPM_MEAN=>NULL()
REAL, DIMENSION(:,:,:), POINTER :: XTKEM_MEAN=>NULL()
REAL, DIMENSION(:,:,:), POINTER :: XPABSM_MEAN=>NULL()
+ REAL, DIMENSION(:,:,:), POINTER :: XQ_MEAN=>NULL()
+ REAL, DIMENSION(:,:,:), POINTER :: XRH_W_MEAN=>NULL()
+ REAL, DIMENSION(:,:,:), POINTER :: XRH_I_MEAN=>NULL()
+ REAL, DIMENSION(:,:,:), POINTER :: XRH_P_MEAN=>NULL()
+ REAL, DIMENSION(:,:) , POINTER :: XRH_W_MAXCOL_MEAN=>NULL()
+ REAL, DIMENSION(:,:) , POINTER :: XRH_I_MAXCOL_MEAN=>NULL()
+ REAL, DIMENSION(:,:) , POINTER :: XRH_P_MAXCOL_MEAN=>NULL()
+ REAL, DIMENSION(:,:,:), POINTER :: XWIFF_MEAN=>NULL()
+ REAL, DIMENSION(:,:,:), POINTER :: XWIDD_MEAN=>NULL()
+ REAL, DIMENSION(:,:,:), POINTER :: XWIFF_MAX=>NULL()
+ REAL, DIMENSION(:,:,:), POINTER :: XWIDD_MAX=>NULL()
REAL, DIMENSION(:,:,:), POINTER :: XSVT_MEAN=>NULL()
REAL, DIMENSION(:,:,:), POINTER :: XU2_MEAN=>NULL(),XV2_MEAN=>NULL(),XW2_MEAN=>NULL(),XUW_MEAN=>NULL()
@@ -72,6 +84,17 @@ REAL, DIMENSION(:,:,:), POINTER :: XTHM_MEAN=>NULL()
REAL, DIMENSION(:,:,:), POINTER :: XTEMPM_MEAN=>NULL()
REAL, DIMENSION(:,:,:), POINTER :: XTKEM_MEAN=>NULL()
REAL, DIMENSION(:,:,:), POINTER :: XPABSM_MEAN=>NULL()
+REAL, DIMENSION(:,:,:), POINTER :: XQ_MEAN=>NULL()
+REAL, DIMENSION(:,:,:), POINTER :: XRH_W_MEAN=>NULL()
+REAL, DIMENSION(:,:,:), POINTER :: XRH_I_MEAN=>NULL()
+REAL, DIMENSION(:,:,:), POINTER :: XRH_P_MEAN=>NULL()
+REAL, DIMENSION(:,:) , POINTER :: XRH_W_MAXCOL_MEAN=>NULL()
+REAL, DIMENSION(:,:) , POINTER :: XRH_I_MAXCOL_MEAN=>NULL()
+REAL, DIMENSION(:,:) , POINTER :: XRH_P_MAXCOL_MEAN=>NULL()
+REAL, DIMENSION(:,:,:), POINTER :: XWIFF_MEAN=>NULL()
+REAL, DIMENSION(:,:,:), POINTER :: XWIDD_MEAN=>NULL()
+REAL, DIMENSION(:,:,:), POINTER :: XWIFF_MAX=>NULL()
+REAL, DIMENSION(:,:,:), POINTER :: XWIDD_MAX=>NULL()
REAL, DIMENSION(:,:,:), POINTER :: XSVT_MEAN=>NULL()
REAL, DIMENSION(:,:,:), POINTER :: XU2_MEAN=>NULL(),XV2_MEAN=>NULL(),XW2_MEAN=>NULL(),XUW_MEAN=>NULL()
@@ -100,6 +123,17 @@ MEAN_FIELD_MODEL(KFROM)%XTHM_MEAN=>XTHM_MEAN
MEAN_FIELD_MODEL(KFROM)%XTEMPM_MEAN=>XTEMPM_MEAN
MEAN_FIELD_MODEL(KFROM)%XTKEM_MEAN=>XTKEM_MEAN
MEAN_FIELD_MODEL(KFROM)%XPABSM_MEAN=>XPABSM_MEAN
+MEAN_FIELD_MODEL(KFROM)%XQ_MEAN=>XQ_MEAN
+MEAN_FIELD_MODEL(KFROM)%XRH_W_MEAN=>XRH_W_MEAN
+MEAN_FIELD_MODEL(KFROM)%XRH_I_MEAN=>XRH_I_MEAN
+MEAN_FIELD_MODEL(KFROM)%XRH_P_MEAN=>XRH_P_MEAN
+MEAN_FIELD_MODEL(KFROM)%XRH_W_MAXCOL_MEAN=>XRH_W_MAXCOL_MEAN
+MEAN_FIELD_MODEL(KFROM)%XRH_I_MAXCOL_MEAN=>XRH_I_MAXCOL_MEAN
+MEAN_FIELD_MODEL(KFROM)%XRH_P_MAXCOL_MEAN=>XRH_P_MAXCOL_MEAN
+MEAN_FIELD_MODEL(KFROM)%XWIFF_MEAN=>XWIFF_MEAN
+MEAN_FIELD_MODEL(KFROM)%XWIDD_MEAN=>XWIDD_MEAN
+MEAN_FIELD_MODEL(KFROM)%XWIFF_MAX=>XWIFF_MAX
+MEAN_FIELD_MODEL(KFROM)%XWIDD_MAX=>XWIDD_MAX
MEAN_FIELD_MODEL(KFROM)%XSVT_MEAN=>XSVT_MEAN
MEAN_FIELD_MODEL(KFROM)%XUM_MAX=>XUM_MAX
@@ -127,6 +161,17 @@ XTHM_MEAN=>MEAN_FIELD_MODEL(KTO)%XTHM_MEAN
XTEMPM_MEAN=>MEAN_FIELD_MODEL(KTO)%XTEMPM_MEAN
XTKEM_MEAN=>MEAN_FIELD_MODEL(KTO)%XTKEM_MEAN
XPABSM_MEAN=>MEAN_FIELD_MODEL(KTO)%XPABSM_MEAN
+XQ_MEAN=>MEAN_FIELD_MODEL(KTO)%XQ_MEAN
+XRH_W_MEAN=>MEAN_FIELD_MODEL(KTO)%XRH_W_MEAN
+XRH_I_MEAN=>MEAN_FIELD_MODEL(KTO)%XRH_I_MEAN
+XRH_P_MEAN=>MEAN_FIELD_MODEL(KTO)%XRH_P_MEAN
+XRH_W_MAXCOL_MEAN=>MEAN_FIELD_MODEL(KTO)%XRH_W_MAXCOL_MEAN
+XRH_I_MAXCOL_MEAN=>MEAN_FIELD_MODEL(KTO)%XRH_I_MAXCOL_MEAN
+XRH_P_MAXCOL_MEAN=>MEAN_FIELD_MODEL(KTO)%XRH_P_MAXCOL_MEAN
+XWIFF_MEAN=>MEAN_FIELD_MODEL(KTO)%XWIFF_MEAN
+XWIDD_MEAN=>MEAN_FIELD_MODEL(KTO)%XWIDD_MEAN
+XWIFF_MAX=>MEAN_FIELD_MODEL(KTO)%XWIFF_MAX
+XWIDD_MAX=>MEAN_FIELD_MODEL(KTO)%XWIDD_MAX
XSVT_MEAN=>MEAN_FIELD_MODEL(KTO)%XSVT_MEAN
XUM_MAX=>MEAN_FIELD_MODEL(KTO)%XUM_MAX
diff --git a/src/MNH/modn_dragbldgn.f90 b/src/MNH/modn_dragbldgn.f90
index bf2a7c7458fb91bbd436c5a17c39e9e6b63e8dcc..eae14620b1c36f6a38b7c8be27296dc87ac2e581 100644
--- a/src/MNH/modn_dragbldgn.f90
+++ b/src/MNH/modn_dragbldgn.f90
@@ -5,48 +5,57 @@
!-----------------------------------------------------------------------------
!!
!! #####################
- MODULE MODN_DRAGBLDG_n
-!! #####################
-!!
-!!*** *MODN_DRAGBLDG*
-!!
-!! PURPOSE
-!! -------
-! Namelist to take into account building drag in the atmospheric model
-! instead of SURFEX.
-!!
-!!** AUTHOR
-!! ------
-!! R.Schoetter *CNRM*
-!
-!! MODIFICATIONS
-!! -------------
-!! Original 09/2019
-!!
-!! IMPLICIT ARGUMENTS
-!! ------------------
-!
-USE MODD_DRAGBLDG_n, ONLY : &
- LDRAGBLDG_n => LDRAGBLDG
-!
-!-----------------------------------------------------------------------------
-!
-!* 0. DECLARATIONS
-! -----------------
-IMPLICIT NONE
-!
-LOGICAL, SAVE :: LDRAGBLDG
-!
-NAMELIST /NAM_DRAGBLDGn/LDRAGBLDG
-!
+MODULE MODN_DRAGBLDG_n
+ !! #####################
+ !!
+ !!*** *MODN_DRAGBLDG*
+ !!
+ !! PURPOSE
+ !! -------
+ ! Namelist to take into account building drag in the atmospheric model
+ ! instead of SURFEX.
+ !!
+ !!** AUTHOR
+ !! ------
+ !! R.Schoetter *CNRM*
+ !
+ !! MODIFICATIONS
+ !! -------------
+ !! Original 09/2019
+ ! R. Schoetter 12/2021 multi-level coupling between MesoNH and SURFEX
+ !!
+ !! IMPLICIT ARGUMENTS
+ !! ------------------
+ !
+ USE MODD_DRAGBLDG_n, ONLY : &
+ LDRAGBLDG_n => LDRAGBLDG, &
+ LFLUXBLDG_n => LFLUXBLDG, &
+ LDRAGURBVEG_n => LDRAGURBVEG
+ !
+ !-----------------------------------------------------------------------------
+ !
+ !* 0. DECLARATIONS
+ ! -----------------
+ IMPLICIT NONE
+ !
+ LOGICAL, SAVE :: LDRAGBLDG
+ LOGICAL, SAVE :: LFLUXBLDG
+ LOGICAL, SAVE :: LDRAGURBVEG
+ !
+ NAMELIST /NAM_DRAGBLDGn/LDRAGBLDG,LFLUXBLDG,LDRAGURBVEG
+ !
CONTAINS
-!
-SUBROUTINE INIT_NAM_DRAGBLDGn
- LDRAGBLDG = LDRAGBLDG_n
-END SUBROUTINE INIT_NAM_DRAGBLDGn
-!
-SUBROUTINE UPDATE_NAM_DRAGBLDGn
- LDRAGBLDG_n = LDRAGBLDG
-END SUBROUTINE UPDATE_NAM_DRAGBLDGn
-!
+ !
+ SUBROUTINE INIT_NAM_DRAGBLDGn
+ LDRAGBLDG = LDRAGBLDG_n
+ LFLUXBLDG = LFLUXBLDG_n
+ LDRAGURBVEG = LDRAGURBVEG_n
+ END SUBROUTINE INIT_NAM_DRAGBLDGn
+ !
+ SUBROUTINE UPDATE_NAM_DRAGBLDGn
+ LDRAGBLDG_n = LDRAGBLDG
+ LFLUXBLDG_n = LFLUXBLDG
+ LDRAGURBVEG_n = LDRAGURBVEG
+ END SUBROUTINE UPDATE_NAM_DRAGBLDGn
+ !
END MODULE MODN_DRAGBLDG_n
diff --git a/src/MNH/phys_paramn.f90 b/src/MNH/phys_paramn.f90
index be5d47e97df870f165baa47c3b3efb1d37f95961..ef8e3543349356e72a6b1cd9a419ceadc8e7e0cc 100644
--- a/src/MNH/phys_paramn.f90
+++ b/src/MNH/phys_paramn.f90
@@ -237,6 +237,7 @@ END MODULE MODI_PHYS_PARAM_n
! C. Lac 11/2019: correction in the drag formula and application to building in addition to tree
! F. Auguste 02/2021: add IBM
! JL Redelsperger 03/2021: add the SW flux penetration for Ocean model case
+ ! R. Schoetter 12/2021 multi-level coupling between MesoNH and SURFEX
!!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -371,6 +372,12 @@ REAL, DIMENSION(:,:), ALLOCATABLE :: ZSFRV ! surface flux of vapor
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZSFSV ! surface flux of scalars
REAL, DIMENSION(:,:), ALLOCATABLE :: ZSFCO2! surface flux of CO2
!
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZSFTH_WALL
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZSFTH_ROOF
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZCD_ROOF
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZSFRV_WALL
+REAL, DIMENSION(:,:), ALLOCATABLE :: ZSFRV_ROOF
+!
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZDIR_ALB ! direct albedo
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZSCA_ALB ! diffuse albedo
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZEMIS ! emissivity
@@ -482,6 +489,12 @@ ALLOCATE(ZSFRV (IIU,IJU))
ALLOCATE(ZSFSV (IIU,IJU,NSV))
ALLOCATE(ZSFCO2(IIU,IJU))
!
+ALLOCATE(ZSFTH_WALL (IIU,IJU))
+ALLOCATE(ZSFTH_ROOF (IIU,IJU))
+ALLOCATE(ZCD_ROOF (IIU,IJU))
+ALLOCATE(ZSFRV_WALL (IIU,IJU))
+ALLOCATE(ZSFRV_ROOF (IIU,IJU))
+!
!* if XWAY(son)=2 save surface fields before radiation or convective scheme
! calls
!
@@ -1222,8 +1235,8 @@ IF (CSURF=='EXTE') THEN
DEALLOCATE( ZSAVE_INPRC,ZSAVE_PRCONV,ZSAVE_PRSCONV)
DEALLOCATE( ZSAVE_DIRFLASWD,ZSAVE_SCAFLASWD,ZSAVE_DIRSRFSWD)
END IF
- CALL GROUND_PARAM_n(ZSFTH, ZSFRV, ZSFSV, ZSFCO2, ZSFU, ZSFV, &
- ZDIR_ALB, ZSCA_ALB, ZEMIS, ZTSRAD )
+ CALL GROUND_PARAM_n(ZSFTH, ZSFTH_WALL, ZSFTH_ROOF, ZCD_ROOF, ZSFRV, ZSFRV_WALL, ZSFRV_ROOF, &
+ ZSFSV, ZSFCO2, ZSFU, ZSFV, ZDIR_ALB, ZSCA_ALB, ZEMIS, ZTSRAD )
!
IF (LIBM) THEN
WHERE(XIBM_LS(:,:,IKB,1).GT.-XIBM_EPSI)
@@ -1264,6 +1277,11 @@ ELSE
ZSFCO2 = 0.
ZSFU = 0.
ZSFV = 0.
+ ZSFTH_WALL = 0.
+ ZSFTH_ROOF = 0.
+ ZCD_ROOF = 0.
+ ZSFRV_WALL = 0.
+ ZSFRV_ROOF = 0.
END IF
!
CALL SECOND_MNH2(ZTIME2)
@@ -1321,11 +1339,53 @@ ZTIME1 = ZTIME2
XTIME_BU_PROCESS = 0.
XTIME_LES_BU_PROCESS = 0.
!
+CALL ADD2DFIELD_ll(TZFIELDS_ll,ZSFTH_WALL, 'PHYS_PARAM_n::ZSFTH_WALL')
+CALL ADD2DFIELD_ll(TZFIELDS_ll,ZSFTH_ROOF, 'PHYS_PARAM_n::ZSFTH_ROOF')
+CALL ADD2DFIELD_ll(TZFIELDS_ll,ZCD_ROOF, 'PHYS_PARAM_n::ZCD_ROOF')
+CALL ADD2DFIELD_ll(TZFIELDS_ll,ZSFRV_WALL, 'PHYS_PARAM_n::ZSFRV_WALL')
+CALL ADD2DFIELD_ll(TZFIELDS_ll,ZSFRV_ROOF, 'PHYS_PARAM_n::ZSFRV_ROOF')
+!
+IF ( CLBCX(1) /= "CYCL" .AND. LWEST_ll()) THEN
+ ZSFTH_WALL(IIB-1,:)=ZSFTH_WALL(IIB,:)
+ ZSFTH_ROOF(IIB-1,:)=ZSFTH_ROOF(IIB,:)
+ ZCD_ROOF (IIB-1,:)=ZCD_ROOF(IIB,:)
+ ZSFRV_WALL(IIB-1,:)=ZSFRV_WALL(IIB,:)
+ ZSFRV_ROOF(IIB-1,:)=ZSFRV_ROOF(IIB,:)
+ENDIF
+!
+IF ( CLBCX(2) /= "CYCL" .AND. LEAST_ll()) THEN
+ ZSFTH_WALL(IIE+1,:)=ZSFTH_WALL(IIE,:)
+ ZSFTH_ROOF(IIE+1,:)=ZSFTH_ROOF(IIE,:)
+ ZCD_ROOF(IIE+1,:) =ZCD_ROOF(IIE,:)
+ ZSFRV_WALL(IIE+1,:)=ZSFRV_WALL(IIE,:)
+ ZSFRV_ROOF(IIE+1,:)=ZSFRV_ROOF(IIE,:)
+ENDIF
+!
+IF ( CLBCY(1) /= "CYCL" .AND. LSOUTH_ll()) THEN
+ ZSFTH_WALL(:,IJB-1)=ZSFTH_WALL(:,IJB)
+ ZSFTH_ROOF(:,IJB-1)=ZSFTH_ROOF(:,IJB)
+ ZCD_ROOF(:,IJB-1) =ZCD_ROOF(:,IJB)
+ ZSFRV_WALL(:,IJB-1)=ZSFRV_WALL(:,IJB)
+ ZSFRV_ROOF(:,IJB-1)=ZSFRV_ROOF(:,IJB)
+ENDIF
+!
+IF ( CLBCY(2) /= "CYCL" .AND. LNORTH_ll()) THEN
+ ZSFTH_WALL(:,IJE+1)=ZSFTH_WALL(:,IJE)
+ ZSFTH_ROOF(:,IJE+1)=ZSFTH_ROOF(:,IJE)
+ ZCD_ROOF(:,IJE+1)=ZCD_ROOF(:,IJE)
+ ZSFRV_WALL(:,IJE+1)=ZSFRV_WALL(:,IJE)
+ ZSFRV_ROOF(:,IJE+1)=ZSFRV_ROOF(:,IJE)
+ENDIF
+!
+!
IF (LDRAGTREE) CALL DRAG_VEG( XTSTEP, XUT, XVT, XTKET, LDEPOTREE, XVDEPOTREE, &
CCLOUD, XPABST, XTHT, XRT, XSVT, XRHODJ, XZZ, &
XRUS, XRVS, XRTKES, XRRS, XRSVS )
!
-IF (LDRAGBLDG) CALL DRAG_BLD( XTSTEP, XUT, XVT, XTKET, XRHODJ, XZZ, XRUS, XRVS, XRTKES )
+IF (LDRAGBLDG) CALL DRAG_BLD ( XTSTEP, XUT, XVT, XTKET, XPABST, XTHT, XRT, XSVT, &
+ XRHODJ, XZZ, XRUS, XRVS, XRTKES, XRTHS, XRRS, &
+ ZSFTH_WALL, ZSFTH_ROOF, ZCD_ROOF, ZSFRV_WALL, &
+ ZSFRV_ROOF )
!
CALL SECOND_MNH2(ZTIME2)
!
@@ -1596,6 +1656,11 @@ DEALLOCATE(ZSFRV )
DEALLOCATE(ZSFSV )
DEALLOCATE(ZSFCO2)
!
+DEALLOCATE(ZSFTH_WALL )
+DEALLOCATE(ZSFTH_ROOF )
+DEALLOCATE(ZCD_ROOF )
+DEALLOCATE(ZSFRV_WALL )
+DEALLOCATE(ZSFRV_ROOF )
!-------------------------------------------------------------------------------
!
END SUBROUTINE PHYS_PARAM_n
diff --git a/src/MNH/read_exsegn.f90 b/src/MNH/read_exsegn.f90
index 4d52929004a0bf4506d4a90b8ca46a503e7a5f59..772e12e9f696debadf804b2278c361847bec9c0d 100644
--- a/src/MNH/read_exsegn.f90
+++ b/src/MNH/read_exsegn.f90
@@ -21,7 +21,7 @@ INTERFACE
OCONDSAMP,OBLOWSNOW, &
KRIMX,KRIMY, KSV_USER, &
HTURB,HTOM,ORMC01,HRAD,HDCONV,HSCONV,HCLOUD,HELEC, &
- HEQNSYS,PTSTEP_ALL,HINIFILEPGD )
+ HEQNSYS,PTSTEP_ALL,HSTORAGE_TYPE,HINIFILEPGD )
!
USE MODD_IO, ONLY: TFILEDATA
!
@@ -71,6 +71,7 @@ CHARACTER (LEN=4), INTENT(IN) :: HCLOUD ! Kind of microphysical scheme
CHARACTER (LEN=4), INTENT(IN) :: HELEC ! Kind of electrical scheme
CHARACTER (LEN=*), INTENT(IN) :: HEQNSYS! type of equations' system
REAL,DIMENSION(:), INTENT(INOUT):: PTSTEP_ALL ! Time STEP of ALL models
+CHARACTER (LEN=*), INTENT(IN) :: HSTORAGE_TYPE ! type of initial file
CHARACTER (LEN=*), INTENT(IN) :: HINIFILEPGD ! name of PGD file
!
END SUBROUTINE READ_EXSEG_n
@@ -93,7 +94,7 @@ END MODULE MODI_READ_EXSEG_n
OCONDSAMP, OBLOWSNOW, &
KRIMX,KRIMY, KSV_USER, &
HTURB,HTOM,ORMC01,HRAD,HDCONV,HSCONV,HCLOUD,HELEC, &
- HEQNSYS,PTSTEP_ALL,HINIFILEPGD )
+ HEQNSYS,PTSTEP_ALL,HSTORAGE_TYPE,HINIFILEPGD )
! #########################################################################
!
!!**** *READ_EXSEG_n * - routine to read the descriptor file EXSEG
@@ -302,6 +303,7 @@ END MODULE MODI_READ_EXSEG_n
! P. Wautelet 10/03/2021: move scalar variable name initializations to ini_nsv
! R. Honnert 23/04/2021: add ADAP mixing length and delete HRIO and BOUT from CMF_UPDRAFT
! S. Riette 11/05/2021: HighLow cloud
+! R. Schoetter 12/2021 multi-level coupling between MesoNH and SURFEX
! P. Wautelet 24/06/2022: remove check on CSTORAGE_TYPE for restart of ForeFire variables
!------------------------------------------------------------------------------
!
@@ -313,6 +315,7 @@ USE MODD_CH_AEROSOL
USE MODD_CH_M9_n, ONLY : NEQ
USE MODD_CONDSAMP
USE MODD_CONF
+USE MODD_CONF_n, ONLY: CSTORAGE_TYPE
USE MODD_CONFZ
! USE MODD_DRAG_n
USE MODD_DUST
@@ -353,6 +356,7 @@ USE MODN_CONF
USE MODN_CONF_n
USE MODN_CONFZ
USE MODN_DRAGBLDG_n
+USE MODN_COUPLING_LEVELS_n
USE MODN_DRAG_n
USE MODN_DRAGTREE_n
USE MODN_DUST
@@ -452,6 +456,7 @@ CHARACTER (LEN=4), INTENT(IN) :: HCLOUD ! Kind of microphysical scheme
CHARACTER (LEN=4), INTENT(IN) :: HELEC ! Kind of electrical scheme
CHARACTER (LEN=*), INTENT(IN) :: HEQNSYS! type of equations' system
REAL,DIMENSION(:), INTENT(INOUT):: PTSTEP_ALL ! Time STEP of ALL models
+CHARACTER (LEN=*), INTENT(IN) :: HSTORAGE_TYPE ! type of initial file
CHARACTER (LEN=*), INTENT(IN) :: HINIFILEPGD ! name of PGD file
!
!* 0.2 declarations of local variables
@@ -478,6 +483,7 @@ CALL INIT_NAM_DYNN
CALL INIT_NAM_ADVN
CALL INIT_NAM_DRAGTREEN
CALL INIT_NAM_DRAGBLDGN
+CALL INIT_NAM_COUPLING_LEVELSN
CALL INIT_NAM_PARAMN
CALL INIT_NAM_PARAM_RADN
#ifdef MNH_ECRAD
@@ -545,6 +551,8 @@ CALL POSNAM(ILUSEG,'NAM_DRAGTREEN',GFOUND,ILUOUT)
IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_DRAGTREEn)
CALL POSNAM(ILUSEG,'NAM_DRAGBLDGN',GFOUND,ILUOUT)
IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_DRAGBLDGn)
+CALL POSNAM(ILUSEG,'NAM_COUPLING_LEVELSN',GFOUND,ILUOUT)
+IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_COUPLING_LEVELSn)
CALL POSNAM(ILUSEG,'NAM_EOL',GFOUND,ILUOUT)
IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_EOL)
CALL POSNAM(ILUSEG,'NAM_EOL_ADNR',GFOUND,ILUOUT)
@@ -2951,6 +2959,7 @@ CALL UPDATE_NAM_LUNITN
CALL UPDATE_NAM_CONFN
CALL UPDATE_NAM_DRAGTREEN
CALL UPDATE_NAM_DRAGBLDGN
+CALL UPDATE_NAM_COUPLING_LEVELSN
CALL UPDATE_NAM_DYNN
CALL UPDATE_NAM_ADVN
CALL UPDATE_NAM_PARAMN
diff --git a/src/MNH/seriesn.f90 b/src/MNH/seriesn.f90
index d49f4805b781b5959425f0ec47bb81e125446bf7..a38d2ccdfb7891ddd4c7ad8be0185c7c293cf18b 100644
--- a/src/MNH/seriesn.f90
+++ b/src/MNH/seriesn.f90
@@ -41,6 +41,8 @@
!! 03/2018 (P.Wautelet) replace TEMPORAL_DIST by DATETIME_DISTANCE
!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
! P. Wautelet 13/09/2019: budget: simplify and modernize date/time management
+!! 12/2021 (R. Schoetter) adds humidity and other mean diagnostics
+
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -127,8 +129,8 @@ IF (LSURF) THEN
KI=(IIE-IIB+1)*(IJE-IJB+1)
ALLOCATE(ZSERIES(KI,5))
CALL GET_SURF_VAR_n(YSURF_CUR%FM,YSURF_CUR%IM,YSURF_CUR%SM,YSURF_CUR%TM, &
- YSURF_CUR%GDM,YSURF_CUR%WM,YSURF_CUR%DUO,YSURF_CUR%DU,&
- YSURF_CUR%UG,YSURF_CUR%U, YSURF_CUR%USS,&
+ YSURF_CUR%GDM, YSURF_CUR%WM,YSURF_CUR%DUO, &
+ YSURF_CUR%DU,YSURF_CUR%UG,YSURF_CUR%U, YSURF_CUR%USS,&
'MESONH',KI,5,PSERIES=ZSERIES)
ZTS(:,:)=0.
ZTMNW(:,:)=0.
diff --git a/src/MNH/write_lfin.f90 b/src/MNH/write_lfin.f90
index e9bab27c0a30fccb9d0b8fbb61d500a8fa9e84a9..dcbe5aef54e4b2c36090ba7e2a949ce30f023fcb 100644
--- a/src/MNH/write_lfin.f90
+++ b/src/MNH/write_lfin.f90
@@ -180,6 +180,7 @@ END MODULE MODI_WRITE_LFIFM_n
! P. Wautelet 10/03/2021: use scalar variable names for dust and salt
! P. Wautelet 11/03/2021: bugfix: correct name for NSV_LIMA_IMM_NUCL
! J.L. Redelsperger 03/2021: add OCEAN and auto-coupled O-A LES cases
+!! R. Schoetter 12/2021 adds humidity and other mean diagnostics
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -268,6 +269,7 @@ USE MODI_SALTLFI_n
USE MODI_CH_AER_REALLFI_n
USE MODI_SALT_FILTER
USE MODI_DUST_FILTER
+USE MODI_SHUMAN
!
!20131128
USE MODE_MPPDB
@@ -765,6 +767,45 @@ IF (MEAN_COUNT /= 0) THEN
TZFIELD%CCOMMENT = 'X_Y_Z_vertical max wind'
CALL IO_Field_write(TPFILE,TZFIELD,XWM_MAX)
!
+ !
+ ! Calculation of mean horizontal wind speed and
+ ! wind direction based on average components
+ !
+ XWIFF_MEAN = SQRT((MXF(XUM_MEAN)/MEAN_COUNT)**2 + (MYF(XVM_MEAN)/MEAN_COUNT)**2)
+ XWIDD_MEAN = 180.0 + (90.0 - 180.0*ATAN2(MYF(XVM_MEAN)/MEAN_COUNT,MXF(XUM_MEAN)/MEAN_COUNT)/XPI)
+ !
+ WHERE (XWIDD_MEAN(:,:,:).GT.360.0)
+ XWIDD_MEAN(:,:,:) = XWIDD_MEAN(:,:,:) - 360.0
+ ENDWHERE
+ !
+ IF ((MINVAL(XWIDD_MEAN).LT.0.0).OR.(MAXVAL(XWIDD_MEAN).GT.360.0)) STOP ("Wrong wind direction")
+ !
+ TZFIELD%CMNHNAME = 'WIFFME'
+ TZFIELD%CLONGNAME = 'WIFFME'
+ TZFIELD%CUNITS = 'm s-1'
+ TZFIELD%CCOMMENT = 'X_Y_Z_horizontal mean wind speed'
+ CALL IO_Field_write(TPFILE,TZFIELD,XWIFF_MEAN)
+ !
+ TZFIELD%CMNHNAME = 'WIDDME'
+ TZFIELD%CLONGNAME = 'WIDDME'
+ TZFIELD%CUNITS = 'm s-1'
+ TZFIELD%CCOMMENT = 'X_Y_Z_horizontal mean wind direction (degrees from north)'
+ CALL IO_Field_write(TPFILE,TZFIELD,XWIDD_MEAN)
+ !
+ TZFIELD%CMNHNAME = 'WIFFMAX'
+ TZFIELD%CLONGNAME = 'WIFFMAX'
+ TZFIELD%CUNITS = 'm s-1'
+ TZFIELD%CCOMMENT = 'X_Y_Z_horizontal maximum wind speed'
+ CALL IO_Field_write(TPFILE,TZFIELD,XWIFF_MAX)
+ XWIFF_MAX(:,:,:)=-XUNDEF
+ !
+ TZFIELD%CMNHNAME = 'WIDDMAX'
+ TZFIELD%CLONGNAME = 'WIDDMAX'
+ TZFIELD%CUNITS = 'm s-1'
+ TZFIELD%CCOMMENT = 'X_Y_Z_horizontal maximum wind direction'
+ CALL IO_Field_write(TPFILE,TZFIELD,XWIDD_MAX)
+ XWIDD_MAX(:,:,:)=-XUNDEF
+!
TZFIELD%NGRID = 1
!
TZFIELD%CMNHNAME = 'CMME'
@@ -813,6 +854,55 @@ IF (MEAN_COUNT /= 0) THEN
TZFIELD%CUNITS = 'K'
TZFIELD%CCOMMENT = 'X_Y_Z_max temperature'
CALL IO_Field_write(TPFILE,TZFIELD,XTEMPM_MAX)
+!
+ TZFIELD%CMNHNAME = 'QSPECME'
+ TZFIELD%CLONGNAME = 'QSPECME'
+ TZFIELD%CUNITS = 'kg/kg'
+ TZFIELD%CCOMMENT = 'X_Y_Z_mean specific humidity'
+ ZWORK3D = XQ_MEAN/MEAN_COUNT
+ CALL IO_Field_write(TPFILE,TZFIELD,ZWORK3D)
+!
+ TZFIELD%CMNHNAME = 'RHME'
+ TZFIELD%CLONGNAME = 'RHME'
+ TZFIELD%CUNITS = '%'
+ TZFIELD%CCOMMENT = 'X_Y_Z_mean relative humidity, water'
+ ZWORK3D = XRH_W_MEAN/MEAN_COUNT
+ CALL IO_Field_write(TPFILE,TZFIELD,ZWORK3D)
+!
+ TZFIELD%CMNHNAME = 'RHME_ICE'
+ TZFIELD%CLONGNAME = 'RHME_ICE'
+ TZFIELD%CUNITS = '%'
+ TZFIELD%CCOMMENT = 'X_Y_Z_mean relative humidity, ice'
+ ZWORK3D = XRH_I_MEAN/MEAN_COUNT
+ CALL IO_Field_write(TPFILE,TZFIELD,ZWORK3D)
+!
+ TZFIELD%CMNHNAME = 'RHME_WEIG'
+ TZFIELD%CLONGNAME = 'RHME_WEIG'
+ TZFIELD%CUNITS = '%'
+ TZFIELD%CCOMMENT = 'X_Y_Z_mean relative humidity, weighted'
+ ZWORK3D = XRH_P_MEAN/MEAN_COUNT
+ CALL IO_Field_write(TPFILE,TZFIELD,ZWORK3D)
+!
+ TZFIELD%CMNHNAME = 'RHMAXME'
+ TZFIELD%CLONGNAME = 'RHMAXME'
+ TZFIELD%CUNITS = '%'
+ TZFIELD%CCOMMENT = 'X_Y_Z_mean column maximum relative humidity, water'
+ ZWORK2D = XRH_W_MAXCOL_MEAN/MEAN_COUNT
+ CALL IO_Field_write(TPFILE,TZFIELD,ZWORK3D)
+!
+ TZFIELD%CMNHNAME = 'RHMAXME_ICE'
+ TZFIELD%CLONGNAME = 'RHMAXME_ICE'
+ TZFIELD%CUNITS = '%'
+ TZFIELD%CCOMMENT = 'X_Y_Z_mean column maximum relative humidity, ice'
+ ZWORK2D = XRH_I_MAXCOL_MEAN/MEAN_COUNT
+ CALL IO_Field_write(TPFILE,TZFIELD,ZWORK3D)
+!
+ TZFIELD%CMNHNAME = 'RHMAXME_WEIG'
+ TZFIELD%CLONGNAME = 'RHMAXME_WEIG'
+ TZFIELD%CUNITS = '%'
+ TZFIELD%CCOMMENT = 'X_Y_Z_mean column maximum relative humidity, weighted'
+ ZWORK2D = XRH_P_MAXCOL_MEAN/MEAN_COUNT
+ CALL IO_Field_write(TPFILE,TZFIELD,ZWORK3D)
!
TZFIELD%CMNHNAME = 'PABSMME'
TZFIELD%CLONGNAME = 'PABSMME'