From 139eb28569b5367d00262a8d39e8a29f023987bf Mon Sep 17 00:00:00 2001
From: Quentin Rodier <quentin.rodier@meteo.fr>
Date: Wed, 27 Nov 2019 18:44:18 +0100
Subject: [PATCH] Christine 27/11/2019: correction in the drag formula and
application to building in addition to tree
---
src/MNH/default_desfmn.f90 | 10 +-
src/MNH/drag_bld.f90 | 235 ++++++++++++
src/MNH/drag_veg.f90 | 578 +++++++++++++++---------------
src/MNH/goto_model_wrapper.f90 | 5 +
src/MNH/ini_modeln.f90 | 3 +-
src/MNH/mnhget_surf_paramn.f90 | 43 ++-
src/MNH/modd_dragbldgn.f90 | 53 +++
src/MNH/modd_dragtreen.f90 | 64 ++++
src/MNH/modn_dragbldgn.f90 | 51 +++
src/MNH/modn_dragtreen.f90 | 61 ++++
src/MNH/phys_paramn.f90 | 34 +-
src/MNH/read_exsegn.f90 | 14 +-
src/SURFEX/allocate_physio.F90 | 7 +-
src/SURFEX/convert_patch_isba.F90 | 3 +-
src/SURFEX/get_surf_varn.F90 | 33 +-
src/SURFEX/get_var_townn.F90 | 42 ++-
src/SURFEX/get_vegn.F90 | 9 +-
17 files changed, 916 insertions(+), 329 deletions(-)
create mode 100644 src/MNH/drag_bld.f90
create mode 100644 src/MNH/modd_dragbldgn.f90
create mode 100644 src/MNH/modd_dragtreen.f90
create mode 100644 src/MNH/modn_dragbldgn.f90
create mode 100644 src/MNH/modn_dragtreen.f90
diff --git a/src/MNH/default_desfmn.f90 b/src/MNH/default_desfmn.f90
index 38c7bba55..934b3447a 100644
--- a/src/MNH/default_desfmn.f90
+++ b/src/MNH/default_desfmn.f90
@@ -229,6 +229,7 @@ END MODULE MODI_DEFAULT_DESFM_n
!! 01/2019 (R. Honnert) add reduction of the mass-flux surface closure with the resolution
!! Bielli S. 02/2019 Sea salt : significant sea wave height influences salt emission; 5 salt modes
!! 05/2019 F.Brient add tracer emission from the top of the boundary-layer
+!! 11/2019 C.Lac correction in the drag formula and application to building in addition to tree
!!
!-------------------------------------------------------------------------------
!
@@ -272,7 +273,8 @@ USE MODD_SALT
USE MODD_PASPOL
USE MODD_CONDSAMP
USE MODD_MEAN_FIELD
-USE MODD_DRAGTREE
+USE MODD_DRAGTREE_n
+USE MODD_DRAGBLDG_n
!
!
USE MODD_PARAM_LIMA, ONLY : LCOLD, LNUCL, LSEDI, LHHONI, LSNOW, LHAIL, LMEYERS,&
@@ -544,6 +546,12 @@ VSIGQSAT = 0.02
LDRAGTREE = .FALSE.
LDEPOTREE = .FALSE.
XVDEPOTREE = 0.02 ! 2 cm/s
+!------------------------------------------------------------------------------
+!
+!* 10c. SET DEFAULT VALUES FOR MODD_DRAGB
+! ----------------------------------
+!
+LDRAGBLDG = .FALSE.
!
!-------------------------------------------------------------------------------
!
diff --git a/src/MNH/drag_bld.f90 b/src/MNH/drag_bld.f90
new file mode 100644
index 000000000..3d3e4cd3e
--- /dev/null
+++ b/src/MNH/drag_bld.f90
@@ -0,0 +1,235 @@
+!MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+! #######################
+MODULE MODI_DRAG_BLD
+ ! #######################
+ !
+ INTERFACE
+ !
+ 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) :: PRHODJ ! dry Density * Jacobian
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! Height (z)
+ !
+ REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRUS, PRVS ! Sources of Momentum
+ REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRTKES ! Sources of Tke
+ !
+ END SUBROUTINE DRAG_BLD
+
+ END INTERFACE
+
+END MODULE MODI_DRAG_BLD
+!
+! ###################################################################
+SUBROUTINE DRAG_BLD(PTSTEP, PUT, PVT, PTKET, PRHODJ, PZZ, PRUS, PRVS, PRTKES )
+ ! ###################################################################
+ !
+ !!**** *DRAG_BLD_n * -
+ !!
+ !! PURPOSE
+ !! -------
+ !
+ ! Drag force due to buildings
+ !
+ !!** METHOD
+ !! ------
+ !!
+ !! REFERENCE
+ !! ---------
+ !!
+ !! AUTHOR
+ !! ------
+ !! R. Schoetter
+ !!
+ !! MODIFICATIONS
+ !! -------------
+ !! Original 09/2019
+ !!---------------------------------------------------------------
+ !
+ !* 0. DECLARATIONS
+ ! ------------
+ !
+ USE MODD_CONF
+ USE MODD_CST
+ USE MODD_DRAGBLDG_n
+ USE MODD_DYN
+ USE MODD_DYN_n
+ USE MODD_GROUND_PAR
+ USE MODD_PGDFIELDS
+ USE MODD_NSV
+ USE MODI_MNHGET_SURF_PARAM_n
+ USE MODI_SHUMAN
+ !
+ IMPLICIT NONE
+ !
+ !* 0.1 Declarations of dummy arguments :
+ !
+ REAL, INTENT(IN) :: PTSTEP ! Time step
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PUT, PVT ! variables
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PTKET ! at t
+ !
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PRHODJ ! dry Density * Jacobian
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PZZ ! Height (z)
+ !
+ REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRUS, PRVS ! Sources of Momentum
+ REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRTKES ! Sources of Tke
+ !
+ !* 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
+ !
+ REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: &
+ ZWORK1, ZWORK2, ZWORK3, ZUT_SCAL, ZVT_SCAL, &
+ ZUS, ZVS, ZTKES, ZTKET
+ 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
+ !
+ IIU = SIZE(PUT,1)
+ IJU = SIZE(PUT,2)
+ IKU = SIZE(PUT,3)
+ !
+ ZUS (:,:,:) = 0.0
+ 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
+ ENDWHERE
+ !
+ ! For buildings, the frontal wall area density is calculated [m^2(walls facing the wind)/m^2]
+ ! The division by PI means that cylindrical buildings are assumed (circle perimeter = PI*D, circle frontal area = D)
+ ! [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
+ !
+ DEALLOCATE(ZFRAC_TOWN_PGD)
+ DEALLOCATE(ZH_BUILD_PGD)
+ DEALLOCATE(ZWALL_O_HOR_PGD)
+ !
+ !-------------------------------------------------------------------------------
+ !
+ !
+ !* 1. COMPUTES THE TRUE VELOCITY COMPONENTS
+ ! -------------------------------------
+ !
+ ZUT_SCAL(:,:,:) = MXF(PUT(:,:,:))
+ ZVT_SCAL(:,:,:) = MYF(PVT(:,:,:))
+ ZTKET(:,:,:) = PTKET(:,:,:)
+ !-------------------------------------------------------------------------------
+ !
+ !* 1. Computations of wind tendency due to canopy drag
+ ! ------------------------------------------------
+ !
+ ! Ext = - Cdrag * u- * u- * Sv tree canopy drag
+ ! - u'w'(ground) * Sh horizontal surfaces (ground)
+ ! ------------------------------
+ !
+ DO JJ=2,(IJU-1)
+ DO JI=2,(IIU-1)
+ !
+ ! Set density and drag coefficient for buildings
+ !
+ IF (ZH_BLD(JI,JJ) /= 0) THEN
+ !
+ DO JK=2,(IKU-1)
+ !
+ IF ( (PZZ(JI,JJ,JK)-PZZ(JI,JJ,2)) .LT. ZH_BLD(JI,JJ) ) THEN
+ !
+ ! FIXME: Check literature values for the drag coefficient here
+ !
+ ZCDRAG(JI,JJ,JK) = 0.2
+ !
+ ! A uniform distribution of building heights is assumed
+ !
+ ZDENSITY(JI,JJ,JK) = ZF_BLD(JI,JJ) / ZH_BLD(JI,JJ)
+ !
+ ENDIF
+ !
+ ENDDO
+ ENDIF
+ !
+ ENDDO
+ ENDDO
+ !
+ !* 1.2 Drag force by wall surfaces
+ ! ---------------------------
+ !
+ !* drag force by vertical surfaces
+ !
+ ZUS(:,:,:) = PUT(:,:,:)/( 1.0 + MXM ( ZCDRAG(:,:,:) * ZDENSITY(:,:,:) &
+ * PTSTEP * SQRT(ZUT_SCAL(:,:,:)**2+ZVT_SCAL(:,:,:)**2) ) )
+ !
+ ZVS(:,:,:) = PVT(:,:,:)/( 1.0 + MYM ( ZCDRAG(:,:,:) * ZDENSITY(:,:,:) &
+ * 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)
+ !
+ ! Ext = + Cd * u+^3 * Sv/Vair vertical surfaces or trees
+ !
+ ! 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
+ !
+ PRTKES(:,:,:) = PRTKES(:,:,:) + (ZTKES(:,:,:)-ZTKET(:,:,:))*PRHODJ(:,:,:)/PTSTEP
+ !
+END SUBROUTINE DRAG_BLD
diff --git a/src/MNH/drag_veg.f90 b/src/MNH/drag_veg.f90
index ce157a154..ec75de6d9 100644
--- a/src/MNH/drag_veg.f90
+++ b/src/MNH/drag_veg.f90
@@ -3,303 +3,311 @@
!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!MNH_LIC for details. version 1.
! #######################
- MODULE MODI_DRAG_VEG
+ MODULE MODI_DRAG_VEG
! #######################
!
INTERFACE
-SUBROUTINE DRAG_VEG(PTSTEP,PUT,PVT,PTKET,ODEPOTREE, PVDEPOTREE, &
- HCLOUD,PPABST,PTHT,PRT,PSVT, &
- PRHODJ,PZZ,PRUS, PRVS, PRTKES, &
- PTHS,PRRS,PSVS)
-!
-REAL, INTENT(IN) :: PTSTEP ! Time step
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PUT, PVT ! variables
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PTKET ! at t
-LOGICAL, INTENT(IN) :: ODEPOTREE ! Droplet deposition on tree
-REAL, INTENT(IN) :: PVDEPOTREE! Velocity deposition on tree
-CHARACTER (LEN=4), INTENT(IN) :: HCLOUD ! Kind of microphysical scheme
-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(INOUT) :: PRUS, PRVS ! Sources of Momentum
-REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRTKES ! Sources of Tke
-REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRRS
-REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PSVS
-REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PTHS
-!
-!
-
-END SUBROUTINE DRAG_VEG
+ SUBROUTINE DRAG_VEG(PTSTEP,PUT,PVT,PTKET,ODEPOTREE, PVDEPOTREE, &
+ HCLOUD,PPABST,PTHT,PRT,PSVT, &
+ PRHODJ,PZZ,PRUS, PRVS, PRTKES, &
+ PTHS,PRRS,PSVS)
+!
+ REAL, INTENT(IN) :: PTSTEP ! Time step
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PUT, PVT ! variables
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PTKET ! at t
+ LOGICAL, INTENT(IN) :: ODEPOTREE ! Droplet deposition on tree
+ REAL, INTENT(IN) :: PVDEPOTREE! Velocity deposition on tree
+ CHARACTER (LEN=4), INTENT(IN) :: HCLOUD ! Kind of microphysical scheme
+ 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(INOUT) :: PRUS, PRVS ! Sources of Momentum
+ REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRTKES ! Sources of Tke
+ REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRRS
+ REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PSVS
+ REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PTHS
+ !
+ END SUBROUTINE DRAG_VEG
-END INTERFACE
+ END INTERFACE
END MODULE MODI_DRAG_VEG
!
! ###################################################################
SUBROUTINE DRAG_VEG(PTSTEP,PUT,PVT,PTKET,ODEPOTREE, PVDEPOTREE, &
- HCLOUD,PPABST,PTHT,PRT,PSVT, &
- PRHODJ,PZZ,PRUS, PRVS, PRTKES, &
- PTHS,PRRS,PSVS)
-! ###################################################################
-!
-!!**** *DRAG_VEG_n * -
-!!
-!! PURPOSE
-!! -------
-!
-!!** METHOD
-!! ------
-!!
-!! REFERENCE
-!! ---------
-!!
-!!
-!! AUTHOR
-!! ------
-!! P. Aumond
-!!
-!! MODIFICATIONS
-!! -------------
-!! Original 07/2009
-!! C.Lac 07/2011 : Add budgets
-!! S. Donier 06/2015 : bug surface aerosols
-!! C.Lac 07/2016 : Add droplet deposition
-!! C.Lac 10/2017 : Correction on deposition
-!!---------------------------------------------------------------
-!
-!
-!* 0. DECLARATIONS
-! ------------
-!
-USE MODD_CONF
-USE MODD_CST
-USE MODD_DYN
-USE MODD_DYN_n
-USE MODD_VEG_n
-USE MODD_BUDGET
-USE MODD_PARAM_C2R2
-USE MODD_NSV
+ HCLOUD,PPABST,PTHT,PRT,PSVT, &
+ PRHODJ,PZZ,PRUS, PRVS, PRTKES, &
+ PTHS,PRRS,PSVS)
+ ! ###################################################################
+ !
+ !!**** *DRAG_VEG_n * -
+ !!
+ !! PURPOSE
+ !! -------
+ !
+ !!** METHOD
+ !! ------
+ !!
+ !! REFERENCE
+ !! ---------
+ !!
+ !!
+ !! AUTHOR
+ !! ------
+ !! P. Aumond
+ !!
+ !! MODIFICATIONS
+ !! -------------
+ !! Original 07/2009
+ !! C.Lac 07/2011 : Add budgets
+ !! S. Donier 06/2015 : bug surface aerosols
+ !! C.Lac 07/2016 : Add droplet deposition
+ !! C.Lac 10/2017 : Correction on deposition
+ !! C.Lac 11/2019 : Correction in the drag formula and application to building in addition to tree
+ !!---------------------------------------------------------------
+ !
+ !
+ !* 0. DECLARATIONS
+ ! ------------
+ !
+ USE MODD_CONF
+ USE MODD_CST
+ USE MODD_DYN
+ USE MODD_DYN_n
+ USE MODD_VEG_n
+ USE MODD_BUDGET
+ USE MODD_PARAM_C2R2
+ USE MODD_NSV
-!
-USE MODI_SHUMAN
-USE MODD_PGDFIELDS
-USE MODD_GROUND_PAR
-USE MODI_MNHGET_SURF_PARAM_n
-USE MODI_BUDGET
+ !
+ USE MODI_SHUMAN
+ USE MODD_PGDFIELDS
+ USE MODD_GROUND_PAR
+ USE MODI_MNHGET_SURF_PARAM_n
+ USE MODI_BUDGET
-!
-IMPLICIT NONE
-!
-!* 0.1 Declarations of dummy arguments :
-!
-REAL, INTENT(IN) :: PTSTEP ! Time step
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PUT, PVT ! variables
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PTKET ! at t
-LOGICAL, INTENT(IN) :: ODEPOTREE ! Droplet deposition on tree
-REAL, INTENT(IN) :: PVDEPOTREE! Velocity deposition on tree
-CHARACTER (LEN=4), INTENT(IN) :: HCLOUD ! Kind of microphysical scheme
-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(INOUT) :: PRUS, PRVS ! Sources of Momentum
-REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRTKES ! Sources of Tke
-REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRRS
-REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PSVS
-REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PTHS
-!
-!
-!* 0.2 Declarations of local variables :
-!
-INTEGER :: IIU,IJU,IKU,IKV ! array size along the k direction
-INTEGER :: JI, JJ, JK ! loop index
-!
-!
-REAL, DIMENSION(:,:), ALLOCATABLE :: ZH_TREE_PGD ! surface cover types
-REAL, DIMENSION(:,:), ALLOCATABLE :: ZLAI_PGD ! surface cover types
-REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: &
- ZWORK1, ZWORK2, ZWORK3, ZUT, ZVT, &
- ZUS, ZVS, ZTKES, ZTKET
-REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: &
- ZCDRAG, ZDENSITY
-REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2)) :: &
- ZVH,ZLAI ! LAI, Vegetation height
-REAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2),SIZE(PZZ,3)):: ZT,ZEXN,ZLV,ZCPH
-LOGICAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) &
- :: GDEP
-REAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2),SIZE(PZZ,3)):: ZWDEPR,ZWDEPS
-
-!
-!
-IIU = SIZE(PUT,1)
-IJU = SIZE(PUT,2)
-IKU = SIZE(PUT,3)
-!
-!* 0.3 Initialisation de kelkes variables
-!
-ZVH(:,:)=0.
-ZLAI(:,:)=0.
-ZCDRAG(:,:,:)=0.
-ZDENSITY(:,:,:)=0.
-!
-ALLOCATE(ZH_TREE_PGD(IIU,IJU))
-ALLOCATE(ZLAI_PGD(IIU,IJU))
-!
-CALL MNHGET_SURF_PARAM_n(PH_TREE=ZH_TREE_PGD,PLAI_TREE=ZLAI_PGD)
-!
-ZVH(:,:)=ZH_TREE_PGD(:,:)
-ZLAI(:,:)=ZLAI_PGD(:,:)
-!
-DEALLOCATE(ZH_TREE_PGD)
-DEALLOCATE(ZLAI_PGD)
-!
-!-------------------------------------------------------------------------------
-!
-!
-!* 1. COMPUTES THE TRUE VELOCITY COMPONENTS
-! -------------------------------------
-!
-ZUT(:,:,:) = PUT(:,:,:)
-ZVT(:,:,:) = PVT(:,:,:)
-ZTKET(:,:,:) = PTKET(:,:,:)
-!-------------------------------------------------------------------------------
-!
-!* 1. Computations of wind tendency due to canopy drag
-! ------------------------------------------------
-!
-!
-!
-! Ext = - Cdrag * u- * u- * Sv tree canopy drag
-! - u'w'(ground) * Sh horizontal surfaces (ground)
-!
-!* 1.1 Drag coefficient by vegetation (Patton et al 2001)
-! ------------------------------
-!
-GDEP(:,:,:) = .FALSE.
-DO JJ=2,(IJU-1)
- DO JI=2,(IIU-1)
- IF (ZVH(JI,JJ) /= 0) THEN
- DO JK=2,(IKU-1)
- IF ((ZVH(JI,JJ)+PZZ(JI,JJ,2))<PZZ(JI,JJ,JK)) EXIT
- IF ((HCLOUD=='C2R2') .OR. (HCLOUD=='KHKO')) THEN
- IF ((PRRS(JI,JJ,JK,2) >0.) .AND. (PSVS(JI,JJ,JK,NSV_C2R2BEG+1) >0.)) &
- GDEP(JI,JJ,JK) = .TRUE.
- ELSE IF (HCLOUD /= 'NONE' .AND. HCLOUD /= 'REVE') THEN
- IF (PRRS(JI,JJ,JK,2) >0.) GDEP(JI,JJ,JK) = .TRUE.
- END IF
- ZCDRAG(JI,JJ,JK) = 0.2 !0.075
- ZDENSITY(JI,JJ,JK) = MAX((4 * (ZLAI(JI,JJ) *&
- (PZZ(JI,JJ,JK)-PZZ(JI,JJ,2)) *&
- (PZZ(JI,JJ,JK)-PZZ(JI,JJ,2)) *&
- (ZVH(JI,JJ)-(PZZ(JI,JJ,JK)-PZZ(JI,JJ,2)))/&
- ZVH(JI,JJ)**3)-&
- (0.30*((ZLAI(JI,JJ) *&
- (PZZ(JI,JJ,JK)-PZZ(JI,JJ,2)) *&
- (PZZ(JI,JJ,JK)-PZZ(JI,JJ,2)) *&
- (PZZ(JI,JJ,JK)-PZZ(JI,JJ,2)) /&
- (ZVH(JI,JJ)**3))-ZLAI(JI,JJ))))/&
- ZVH(JI,JJ), 0.)
-
-
- END DO
- END IF
- END DO
-END DO
-! To exclude the first vertical level already dealt in rain_ice or rain_c2r2_khko
-GDEP(:,:,2) = .FALSE.
-!
-!* 1.2 Drag force by wall surfaces
-! ---------------------------
-!
-!* drag force by vertical surfaces
-!
-ZUS(:,:,:)= ZUT(:,:,:)/(1 + ZCDRAG(:,:,:)* ZDENSITY(:,:,:)*PTSTEP &
- *SQRT(ZUT(:,:,:)**2+ZVT(:,:,:)**2))
-!
-ZVS(:,:,:)= ZVT(:,:,:)/(1 + ZCDRAG(:,:,:)* ZDENSITY(:,:,:)*PTSTEP &
- *SQRT(ZUT(:,:,:)**2+ZVT(:,:,:)**2))
-!
-PRUS(:,:,:)=PRUS(:,:,:)+((ZUS(:,:,:)-ZUT(:,:,:))*PRHODJ(:,:,:))/PTSTEP
-!
-PRVS(:,:,:)=PRVS(:,:,:)+((ZVS(:,:,:)-ZVT(:,:,:))*PRHODJ(:,:,:))/PTSTEP
-!
-IF (ODEPOTREE) THEN
- ZEXN(:,:,:)= (PPABST(:,:,:)/XP00)**(XRD/XCPD)
- ZT(:,:,:)= PTHT(:,:,:)*ZEXN(:,:,:)
- ZLV(:,:,:)=XLVTT +(XCPV-XCL) *(ZT(:,:,:)-XTT)
- ZCPH(:,:,:)=XCPD +XCPV*PRT(:,:,:,1)
- ZWDEPR(:,:,:)= 0.
- ZWDEPS(:,:,:)= 0.
- WHERE (GDEP)
- ZWDEPR(:,:,:)= PVDEPOTREE * PRT(:,:,:,2) * PRHODJ(:,:,:)
- END WHERE
- IF ((HCLOUD=='C2R2') .OR. (HCLOUD=='KHKO') .OR. (HCLOUD=='LIMA')) THEN
- WHERE (GDEP)
- ZWDEPS(:,:,:)= PVDEPOTREE * PSVT(:,:,:,NSV_C2R2BEG+1) * PRHODJ(:,:,:)
- END WHERE
- END IF
+ !
+ IMPLICIT NONE
+ !
+ !* 0.1 Declarations of dummy arguments :
+ !
+ REAL, INTENT(IN) :: PTSTEP ! Time step
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PUT, PVT ! variables
+ REAL, DIMENSION(:,:,:), INTENT(IN) :: PTKET ! at t
+ LOGICAL, INTENT(IN) :: ODEPOTREE ! Droplet deposition on tree
+ REAL, INTENT(IN) :: PVDEPOTREE! Velocity deposition on tree
+ CHARACTER (LEN=4), INTENT(IN) :: HCLOUD ! Kind of microphysical scheme
+ 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(INOUT) :: PRUS, PRVS ! Sources of Momentum
+ REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PRTKES ! Sources of Tke
+ REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PRRS
+ REAL, DIMENSION(:,:,:,:), INTENT(INOUT) :: PSVS
+ REAL, DIMENSION(:,:,:), INTENT(INOUT) :: PTHS
+ !
+ !* 0.2 Declarations of local variables :
+ !
+ INTEGER :: IIU,IJU,IKU,IKV ! array size along the k direction
+ INTEGER :: JI, JJ, JK ! loop index
+ !
+ REAL, DIMENSION(:,:), ALLOCATABLE :: ZH_TREE_PGD ! surface cover types
+ REAL, DIMENSION(:,:), ALLOCATABLE :: ZLAI_PGD ! surface cover types
+ REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: &
+ ZWORK1, ZWORK2, ZWORK3, ZUT_SCAL, ZVT_SCAL, &
+ ZUS, ZVS, ZTKES, ZTKET
+ REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: &
+ ZCDRAG, ZDENSITY
+ !
+ REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2)) :: &
+ ZH,ZLAI ! LAI, Vegetation height
+ REAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2),SIZE(PZZ,3)):: ZT,ZEXN,ZLV,ZCPH
+ LOGICAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) &
+ :: GDEP
+ REAL, DIMENSION(SIZE(PZZ,1),SIZE(PZZ,2),SIZE(PZZ,3)):: ZWDEPR,ZWDEPS
+ !
+ !* 0.3 Initialisation de kelkes variables
+ !
+ IIU = SIZE(PUT,1)
+ IJU = SIZE(PUT,2)
+ IKU = SIZE(PUT,3)
+ !
+ ZUS (:,:,:) = 0.0
+ ZVS (:,:,:) = 0.0
+ ZTKES (:,:,:) = 0.0
+ !
+ ZH (:,:) = 0.
+ ZLAI (:,:) = 0.
+ !
+ ZCDRAG (:,:,:) = 0.
+ ZDENSITY (:,:,:) = 0.
+ !
+ ALLOCATE(ZH_TREE_PGD(IIU,IJU))
+ ALLOCATE(ZLAI_PGD(IIU,IJU))
+ !
+ ZH_TREE_PGD (:,:) = XUNDEF
+ ZLAI_PGD (:,:) = XUNDEF
+ !
+ CALL MNHGET_SURF_PARAM_n( PH_TREE = ZH_TREE_PGD, PLAI_TREE = ZLAI_PGD )
+ !
+ ZH (:,:) = ZH_TREE_PGD(:,:)
+ ZLAI(:,:) = ZLAI_PGD(:,:)
+ !
+ WHERE ( ZH (:,:).GT.998.0 ) ZH (:,:) = 0.0
+ WHERE ( ZLAI (:,:).GT.998.0 ) ZLAI (:,:) = 0.0
+ !
+ DEALLOCATE(ZH_TREE_PGD)
+ DEALLOCATE(ZLAI_PGD)
+ !
+ !-------------------------------------------------------------------------------
+ !
+ !
+ !* 1. COMPUTES THE TRUE VELOCITY COMPONENTS
+ ! -------------------------------------
+ !
+ ZUT_SCAL(:,:,:) = MXF(PUT(:,:,:))
+ ZVT_SCAL(:,:,:) = MYF(PVT(:,:,:))
+ ZTKET(:,:,:) = PTKET(:,:,:)
+ !-------------------------------------------------------------------------------
+ !
+ !* 1. Computations of wind tendency due to canopy drag
+ ! ------------------------------------------------
+ !
+ !
+ !
+ ! Ext = - Cdrag * u- * u- * Sv tree canopy drag
+ ! - u'w'(ground) * Sh horizontal surfaces (ground)
+ !
+ !* 1.1 Drag coefficient by vegetation (Patton et al 2001)
+ ! ------------------------------
+ !
+ GDEP(:,:,:) = .FALSE.
+ !
DO JJ=2,(IJU-1)
- DO JI=2,(IIU-1)
- DO JK=2,(IKU-2)
- IF (GDEP(JI,JJ,JK)) THEN
- PRRS(JI,JJ,JK,2) = PRRS(JI,JJ,JK,2) + (ZWDEPR(JI,JJ,JK+1)-ZWDEPR(JI,JJ,JK))/ &
- (PZZ(JI,JJ,JK+1)-PZZ(JI,JJ,JK))
- IF ((HCLOUD=='C2R2') .OR. (HCLOUD=='KHKO').OR. (HCLOUD=='LIMA')) THEN
- PSVS(JI,JJ,JK,NSV_C2R2BEG+1) = PSVS(JI,JJ,JK,NSV_C2R2BEG+1) + &
- (ZWDEPS(JI,JJ,JK+1)-ZWDEPS(JI,JJ,JK))/(PZZ(JI,JJ,JK+1)-PZZ(JI,JJ,JK))
- END IF
- END IF
+ DO JI=2,(IIU-1)
+ !
+ ! Set density and drag coefficient for vegetation
+ !
+ IF (ZH(JI,JJ) /= 0) THEN
+ !
+ DO JK=2,(IKU-1)
+ !
+ IF ( (PZZ(JI,JJ,JK)-PZZ(JI,JJ,2)) .LT. ZH(JI,JJ) ) THEN
+ !
+ IF ((HCLOUD=='C2R2') .OR. (HCLOUD=='KHKO')) THEN
+ IF ((PRRS(JI,JJ,JK,2) >0.) .AND. (PSVS(JI,JJ,JK,NSV_C2R2BEG+1) >0.)) &
+ GDEP(JI,JJ,JK) = .TRUE.
+ ELSE IF (HCLOUD /= 'NONE' .AND. HCLOUD /= 'REVE') THEN
+ IF (PRRS(JI,JJ,JK,2) >0.) GDEP(JI,JJ,JK) = .TRUE.
+ ENDIF
+ !
+ ZCDRAG(JI,JJ,JK) = 0.2 !0.075
+ ZDENSITY(JI,JJ,JK) = MAX((4 * (ZLAI(JI,JJ) *&
+ (PZZ(JI,JJ,JK)-PZZ(JI,JJ,2)) *&
+ (PZZ(JI,JJ,JK)-PZZ(JI,JJ,2)) *&
+ (ZH(JI,JJ)-(PZZ(JI,JJ,JK)-PZZ(JI,JJ,2)))/&
+ ZH(JI,JJ)**3)-&
+ (0.30*((ZLAI(JI,JJ) *&
+ (PZZ(JI,JJ,JK)-PZZ(JI,JJ,2)) *&
+ (PZZ(JI,JJ,JK)-PZZ(JI,JJ,2)) *&
+ (PZZ(JI,JJ,JK)-PZZ(JI,JJ,2)) /&
+ (ZH(JI,JJ)**3))-ZLAI(JI,JJ))))/&
+ ZH(JI,JJ), 0.)
+ !
+ ENDIF
+ !
+ ENDDO
+ ENDIF
+ !
+ ENDDO
+ ENDDO
+ !
+ ! To exclude the first vertical level already dealt in rain_ice or rain_c2r2_khko
+ GDEP(:,:,2) = .FALSE.
+ !
+ !* 1.2 Drag force by wall surfaces
+ ! ---------------------------
+ !
+ !* drag force by vertical surfaces
+ !
+ ZUS(:,:,:) = PUT(:,:,:)/( 1.0 + MXM ( ZCDRAG(:,:,:) * ZDENSITY(:,:,:) &
+ * PTSTEP * SQRT(ZUT_SCAL(:,:,:)**2+ZVT_SCAL(:,:,:)**2) ) )
+ !
+ ZVS(:,:,:) = PVT(:,:,:)/( 1.0 + MYM ( ZCDRAG(:,:,:) * ZDENSITY(:,:,:) &
+ * PTSTEP * SQRT(ZUT_SCAL(:,:,:)**2+ZVT_SCAL(:,:,:)**2) ) )
+ !
+ PRUS(:,:,:) = PRUS(:,:,:) + (ZUS(:,:,:)-PUT(:,:,:)) * MXM(PRHODJ(:,:,:)) / PTSTEP
+ !
+ PRVS(:,:,:) = PRVS(:,:,:) + (ZVS(:,:,:)-PVT(:,:,:)) * MYM(PRHODJ(:,:,:)) / PTSTEP
+ !
+ IF (ODEPOTREE) THEN
+ ZEXN(:,:,:)= (PPABST(:,:,:)/XP00)**(XRD/XCPD)
+ ZT(:,:,:)= PTHT(:,:,:)*ZEXN(:,:,:)
+ ZLV(:,:,:)=XLVTT +(XCPV-XCL) *(ZT(:,:,:)-XTT)
+ ZCPH(:,:,:)=XCPD +XCPV*PRT(:,:,:,1)
+ ZWDEPR(:,:,:)= 0.
+ ZWDEPS(:,:,:)= 0.
+ WHERE (GDEP)
+ ZWDEPR(:,:,:)= PVDEPOTREE * PRT(:,:,:,2) * PRHODJ(:,:,:)
+ END WHERE
+ IF ((HCLOUD=='C2R2') .OR. (HCLOUD=='KHKO') .OR. (HCLOUD=='LIMA')) THEN
+ WHERE (GDEP)
+ ZWDEPS(:,:,:)= PVDEPOTREE * PSVT(:,:,:,NSV_C2R2BEG+1) * PRHODJ(:,:,:)
+ END WHERE
+ END IF
+ DO JJ=2,(IJU-1)
+ DO JI=2,(IIU-1)
+ DO JK=2,(IKU-2)
+ IF (GDEP(JI,JJ,JK)) THEN
+ PRRS(JI,JJ,JK,2) = PRRS(JI,JJ,JK,2) + (ZWDEPR(JI,JJ,JK+1)-ZWDEPR(JI,JJ,JK))/ &
+ (PZZ(JI,JJ,JK+1)-PZZ(JI,JJ,JK))
+ IF ((HCLOUD=='C2R2') .OR. (HCLOUD=='KHKO').OR. (HCLOUD=='LIMA')) THEN
+ PSVS(JI,JJ,JK,NSV_C2R2BEG+1) = PSVS(JI,JJ,JK,NSV_C2R2BEG+1) + &
+ (ZWDEPS(JI,JJ,JK+1)-ZWDEPS(JI,JJ,JK))/(PZZ(JI,JJ,JK+1)-PZZ(JI,JJ,JK))
+ END IF
+ END IF
+ END DO
+ END DO
END DO
- END DO
- END DO
-!
-!
-END IF
-!
-IF (LBUDGET_U) CALL BUDGET (PRUS,1,'DRAG_BU_RU')
-IF (LBUDGET_V) CALL BUDGET (PRVS,2,'DRAG_BU_RV')
-IF (LBUDGET_RC) CALL BUDGET (PRRS(:,:,:,2),7,'DEPOTR_BU_RRC')
-IF (LBUDGET_SV) CALL BUDGET (PSVS(:,:,:,NSV_C2R2BEG+1),14+(NSV_C2R2BEG-1),'DEPOTR_BU_RSV')
-!
-!
-!* 3. Computations of TKE tendency due to canopy drag
-! ------------------------------------------------
-
-!* 3.1 Creation of TKE by wake
-! -----------------------
-!
-! from Kanda and Hino (1994)
-!
-! Ext = + Cd * u+^3 * Sv/Vair vertical surfaces or trees
-! Ext = - Cd * e * u * Sv trees Destruction of TKE due to
-! small-scale motions forced by leaves from Kanda and Hino (1994)
-!
-! with Vair = Vair/Vtot * Vtot = (Vair/Vtot) * Stot * Dz
-! and Sv/Vair = (Sv/Stot) * Stot/Vair = (Sv/Stot) / (Vair/Vtot) / Dz
-!
-!ZTKES(:,:,:)= (ZTKET(:,:,:) + (ZCDRAG(:,:,:)* ZDENSITY(:,:,:) &
-! *(SQRT(ZUT(:,:,:)**2+ZVT(:,:,:)**2))**3)) /&
-! (1.+(2.*ZCDRAG(:,:,:)* ZDENSITY(:,:,:)*SQRT(ZUT(:,:,:)**2+ZVT(:,:,:)**2)))
-ZTKES(:,:,:)= (ZTKET(:,:,:) + (ZCDRAG(:,:,:)* ZDENSITY(:,:,:) &
- *(SQRT(ZUT(:,:,:)**2+ZVT(:,:,:)**2))**3))*PTSTEP /&
- (1.+PTSTEP*ZCDRAG(:,:,:)* ZDENSITY(:,:,:)*SQRT(ZUT(:,:,:)**2+ZVT(:,:,:)**2))
-!
-PRTKES(:,:,:)=PRTKES(:,:,:)+((ZTKES(:,:,:)-ZTKET(:,:,:))*PRHODJ(:,:,:)/PTSTEP)
-!
-IF (LBUDGET_TKE) CALL BUDGET (PRTKES(:,:,:),5,'DRAG_BU_RTKE')
-!
+ !
+ END IF
+ !
+ IF (LBUDGET_U) CALL BUDGET (PRUS,1,'DRAG_BU_RU')
+ IF (LBUDGET_V) CALL BUDGET (PRVS,2,'DRAG_BU_RV')
+ IF (LBUDGET_RC) CALL BUDGET (PRRS(:,:,:,2),7,'DEPOTR_BU_RRC')
+ IF (LBUDGET_SV) CALL BUDGET (PSVS(:,:,:,NSV_C2R2BEG+1),14+(NSV_C2R2BEG-1),'DEPOTR_BU_RSV')
+ !
+ !* 3. Computations of TKE tendency due to canopy drag
+ ! ------------------------------------------------
+ !* 3.1 Creation of TKE by wake
+ ! -----------------------
+ !
+ ! from Kanda and Hino (1994)
+ !
+ ! Ext = + Cd * u+^3 * Sv/Vair vertical surfaces or trees
+ ! Ext = - Cd * e * u * Sv trees Destruction of TKE due to
+ ! small-scale motions forced by leaves from Kanda and Hino (1994)
+ !
+ ! 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 ) / &
+ ( 1. + PTSTEP * ZCDRAG(:,:,:) * ZDENSITY(:,:,:) * SQRT(ZUT_SCAL(:,:,:)**2+ZVT_SCAL(:,:,:)**2))
+ !
+ PRTKES(:,:,:) = PRTKES(:,:,:) + (ZTKES(:,:,:)-ZTKET(:,:,:))*PRHODJ(:,:,:)/PTSTEP
+ !
+ IF (LBUDGET_TKE) CALL BUDGET (PRTKES(:,:,:),5,'DRAG_BU_RTKE')
+ !
END SUBROUTINE DRAG_VEG
diff --git a/src/MNH/goto_model_wrapper.f90 b/src/MNH/goto_model_wrapper.f90
index 1d3e057e3..e58e193dc 100644
--- a/src/MNH/goto_model_wrapper.f90
+++ b/src/MNH/goto_model_wrapper.f90
@@ -15,6 +15,7 @@
! 02/2018 Q.Libois ECRAD
!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
! 2017 V.Vionnet blow snow
+! 11/2019 C.Lac correction in the drag formula and application to building in addition to tree
!-----------------------------------------------------------------
MODULE MODI_GOTO_MODEL_WRAPPER
@@ -42,6 +43,8 @@ USE MODD_CONF_n
USE MODD_CURVCOR_n
!USE MODD_DEEP_CONVECTION_n
USE MODD_DIM_n
+USE MODD_DRAGTREE_n
+USE MODD_DRAGBLDG_n
USE MODD_DUMMY_GR_FIELD_n
USE MODD_DYN_n
USE MODD_DYNZD_n
@@ -151,6 +154,8 @@ CALL CONF_GOTO_MODEL(KFROM, KTO)
CALL CURVCOR_GOTO_MODEL(KFROM, KTO)
!CALL DEEP_CONVECTION_GOTO_MODEL(KFROM, KTO)
CALL DIM_GOTO_MODEL(KFROM, KTO)
+CALL DRAGTREE_GOTO_MODEL(KFROM, KTO)
+CALL DRAGBLDG_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/ini_modeln.f90 b/src/MNH/ini_modeln.f90
index 74312b192..63ef283bf 100644
--- a/src/MNH/ini_modeln.f90
+++ b/src/MNH/ini_modeln.f90
@@ -284,6 +284,7 @@ END MODULE MODI_INI_MODEL_n
!! 02/2019 C.Lac add rain fraction as an output field
!! Bielli S. 02/2019 Sea salt : significant sea wave height influences salt emission; 5 salt modes
! P. Wautelet 14/03/2019: correct ZWS when variable not present in file (set to XZWS_DEFAULT)
+! 11/2019 C.Lac correction in the drag formula and application to building in addition to tree
!---------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -317,7 +318,7 @@ USE MODD_TIME
USE MODD_TURB_CLOUD, ONLY: NMODEL_CLOUD, CTURBLEN_CLOUD,XCEI
USE MODD_NESTING
USE MODD_PASPOL
-USE MODD_DRAGTREE
+USE MODD_DRAGTREE_n
USE MODD_METRICS_n
USE MODD_DYN_n
USE MODD_DYNZD_n
diff --git a/src/MNH/mnhget_surf_paramn.f90 b/src/MNH/mnhget_surf_paramn.f90
index 238c9d4ae..dbe76b4f3 100644
--- a/src/MNH/mnhget_surf_paramn.f90
+++ b/src/MNH/mnhget_surf_paramn.f90
@@ -14,15 +14,19 @@
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 )
+ PBARE, PLAI_TREE, PH_TREE, PWALL_O_HOR, &
+ PBUILD_HEIGHT,PNATURE )
!
REAL, DIMENSION(:,:,:), INTENT(OUT), OPTIONAL :: PCOVER ! cover types
REAL, DIMENSION(:,:), INTENT(OUT), OPTIONAL :: PSEA ! sea fraction
REAL, DIMENSION(:,:), INTENT(OUT), OPTIONAL :: PTOWN ! town fraction
+REAL, DIMENSION(:,:), INTENT(OUT), OPTIONAL :: PNATURE ! nature fraction
INTEGER, INTENT(OUT), OPTIONAL :: KCOVER ! number of cover types
REAL, DIMENSION(:,:), INTENT(OUT), OPTIONAL :: PBARE ! Bare soil fraction
-REAL, DIMENSION(:,:), INTENT(OUT), OPTIONAL :: PLAI_TREE !
-REAL, DIMENSION(:,:), INTENT(OUT), OPTIONAL :: PH_TREE
+REAL, DIMENSION(:,:), INTENT(OUT), OPTIONAL :: PLAI_TREE ! Tree leaf area index [m^2(leaf)/m^2(nature)]
+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(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)
@@ -43,7 +47,8 @@ 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 )
+ PBARE, PLAI_TREE, PH_TREE, PWALL_O_HOR, &
+ PBUILD_HEIGHT,PNATURE )
! ########################################
!
!!**** *MNHGET_SURF_PARAM_n* - gets some surface fields on MESONH grid
@@ -79,6 +84,7 @@ END MODULE MODI_MNHGET_SURF_PARAM_n
!! S. Donier 06/2015 : bug surface aerosols
!! 06/2016 (G.Delautier) phasage surfex 8
!! 01/2018 (G.Delautier) SURFEX 8.1
+!! 11/2019 C.Lac correction in the drag formula and application to building in addition to tree
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -104,10 +110,13 @@ IMPLICIT NONE
REAL, DIMENSION(:,:,:), INTENT(OUT), OPTIONAL :: PCOVER ! cover types
REAL, DIMENSION(:,:), INTENT(OUT), OPTIONAL :: PSEA ! sea fraction
REAL, DIMENSION(:,:), INTENT(OUT), OPTIONAL :: PTOWN ! town fraction
+REAL, DIMENSION(:,:), INTENT(OUT), OPTIONAL :: PNATURE ! nature fraction
INTEGER, INTENT(OUT), OPTIONAL :: KCOVER ! number of cover types
REAL, DIMENSION(:,:), INTENT(OUT), OPTIONAL :: PBARE ! Bare soil fraction
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(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)
@@ -140,6 +149,8 @@ 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 :: ZBUILD_HEIGHT ! Building height [m]
REAL, DIMENSION(:), ALLOCATABLE :: ZBARE ! bare soil fraction
REAL, DIMENSION(:), ALLOCATABLE :: ZZS ! orography
REAL, DIMENSION(:), ALLOCATABLE :: ZRN ! net radiation at surface (W/m2)
@@ -173,8 +184,9 @@ IF (PRESENT(PCOVER)) THEN
DEALLOCATE(ZCOVER)
END IF
!
-IF (PRESENT(PSEA) .OR. PRESENT(PTOWN) .OR. &
- PRESENT(PBARE) .OR. PRESENT(PLAI_TREE) .OR. PRESENT(PH_TREE)) 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) ) THEN
ALLOCATE(ZSEA ( ILU ))
ALLOCATE(ZWATER ( ILU ))
ALLOCATE(ZNATURE( ILU ))
@@ -186,6 +198,9 @@ IF (PRESENT(PSEA) .OR. PRESENT(PTOWN) .OR. &
IF (PRESENT(PTOWN)) THEN
CALL REMOVE_HALO(ZTOWN,PTOWN)
END IF
+ IF (PRESENT(PNATURE)) THEN
+ CALL REMOVE_HALO(ZNATURE,PNATURE)
+ END IF
END IF
!
IF (PRESENT(PBARE)) THEN
@@ -251,6 +266,22 @@ IF (PRESENT(PH_TREE) .OR.PRESENT(PLAI_TREE)) THEN
DEALLOCATE(ZLAI)
END IF
!
+IF (PRESENT(PWALL_O_HOR) .OR. PRESENT(PBUILD_HEIGHT)) THEN
+ PBUILD_HEIGHT(:,:) = XUNDEF
+ 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%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 )
+ CALL REMOVE_HALO(ZBUILD_HEIGHT,PBUILD_HEIGHT)
+ CALL REMOVE_HALO(ZWALL_O_HOR,PWALL_O_HOR)
+ DEALLOCATE(ZBUILD_HEIGHT)
+ DEALLOCATE(ZWALL_O_HOR)
+END IF
+!
IF (ALLOCATED(ZSEA)) THEN
DEALLOCATE(ZSEA )
DEALLOCATE(ZWATER )
diff --git a/src/MNH/modd_dragbldgn.f90 b/src/MNH/modd_dragbldgn.f90
new file mode 100644
index 000000000..3bcd2e590
--- /dev/null
+++ b/src/MNH/modd_dragbldgn.f90
@@ -0,0 +1,53 @@
+!MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+!!
+!! #####################
+ MODULE MODD_DRAGBLDG_n
+!! #####################
+!!
+!!*** *MODD_DRAGBLDG*
+!!
+!! PURPOSE
+!! -------
+! Declaration to take into account building drag in Meso-NH instead of SURFEX.
+!!
+!!** AUTHOR
+!! ------
+!! R.Schoetter *CNRM*
+!
+!! MODIFICATIONS
+!! -------------
+!! Original 09/2019
+!-----------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! -----------------
+!
+USE MODD_PARAMETERS, ONLY: JPMODELMAX
+!
+IMPLICIT NONE
+!
+TYPE DRAGBLDG_t
+ !
+ LOGICAL :: LDRAGBLDG ! flag used to take into account building drag in
+ ! ! the atmospheric model instead of SURFEX.
+ !
+END TYPE DRAGBLDG_t
+!
+TYPE(DRAGBLDG_t), DIMENSION(JPMODELMAX), TARGET, SAVE :: DRAGBLDG_MODEL
+!
+LOGICAL, POINTER :: LDRAGBLDG=>NULL()
+!
+CONTAINS
+!
+SUBROUTINE DRAGBLDG_GOTO_MODEL(KFROM, KTO)
+ !
+ INTEGER, INTENT(IN) :: KFROM, KTO
+ !
+ LDRAGBLDG=>DRAGBLDG_MODEL(KTO)%LDRAGBLDG
+ !
+END SUBROUTINE DRAGBLDG_GOTO_MODEL
+!
+END MODULE MODD_DRAGBLDG_n
diff --git a/src/MNH/modd_dragtreen.f90 b/src/MNH/modd_dragtreen.f90
new file mode 100644
index 000000000..6f957b4a0
--- /dev/null
+++ b/src/MNH/modd_dragtreen.f90
@@ -0,0 +1,64 @@
+!MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+!!
+!! #####################
+ MODULE MODD_DRAGTREE_n
+!! #####################
+!!
+!!*** *MODD_DRAGTREE*
+!!
+!! PURPOSE
+!! -------
+! Declaration to take into account tree drag in Meso-NH
+! instead of SURFEX.
+!!
+!!** AUTHOR
+!! ------
+!! C.Lac *CNRM*
+!
+!! MODIFICATIONS
+!! -------------
+!! Original 30/06/11
+!! 06/16 (C.Lac) Add droplet deposition
+!! 11/2019 C.Lac correction in the drag formula and application to building in addition to tree
+
+!-----------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! -----------------
+!
+USE MODD_PARAMETERS, ONLY: JPMODELMAX
+!
+IMPLICIT NONE
+!
+TYPE DRAGTREE_t
+ !
+ LOGICAL :: LDRAGTREE ! flag used to take into account tree drag in
+ ! ! the atmospheric model instead of SURFEX.
+ LOGICAL :: LDEPOTREE ! flag for droplet deposition on trees
+ !
+ REAL :: XVDEPOTREE ! Droplet deposition velocity
+ !
+END TYPE DRAGTREE_t
+!
+TYPE(DRAGTREE_t), DIMENSION(JPMODELMAX), TARGET, SAVE :: DRAGTREE_MODEL
+!
+LOGICAL,POINTER :: LDRAGTREE=>NULL()
+LOGICAL,POINTER :: LDEPOTREE=>NULL()
+REAL ,POINTER :: XVDEPOTREE=>NULL()
+!
+CONTAINS
+!
+SUBROUTINE DRAGTREE_GOTO_MODEL(KFROM, KTO)
+ !
+ INTEGER, INTENT(IN) :: KFROM, KTO
+ !
+ LDRAGTREE=>DRAGTREE_MODEL(KTO)%LDRAGTREE
+ LDEPOTREE=>DRAGTREE_MODEL(KTO)%LDEPOTREE
+ XVDEPOTREE=>DRAGTREE_MODEL(KTO)%XVDEPOTREE
+ !
+END SUBROUTINE DRAGTREE_GOTO_MODEL
+!
+END MODULE MODD_DRAGTREE_n
diff --git a/src/MNH/modn_dragbldgn.f90 b/src/MNH/modn_dragbldgn.f90
new file mode 100644
index 000000000..d60daa31e
--- /dev/null
+++ b/src/MNH/modn_dragbldgn.f90
@@ -0,0 +1,51 @@
+!MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+!!
+!! #####################
+ 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
+!
+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
+!
+END MODULE MODN_DRAGBLDG_n
diff --git a/src/MNH/modn_dragtreen.f90 b/src/MNH/modn_dragtreen.f90
new file mode 100644
index 000000000..b0654cc6f
--- /dev/null
+++ b/src/MNH/modn_dragtreen.f90
@@ -0,0 +1,61 @@
+!MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+!!
+!! #####################
+ MODULE MODN_DRAGTREE_n
+!! #####################
+!!
+!!*** *MODN_DRAGTREE*
+!!
+!! PURPOSE
+!! -------
+! Namelist to take into account tree drag in the atmospheric model
+! instead of SURFEX.
+!!
+!!** AUTHOR
+!! ------
+!! C.Lac *CNRM*
+!
+!! MODIFICATIONS
+!! -------------
+!! Original 30/06/11
+!!
+!! 10/2016 : (C.Lac) Add droplet deposition on trees
+!! 11/2019 C.Lac correction in the drag formula and application to building in addition to tree
+!! IMPLICIT ARGUMENTS
+!! ------------------
+!
+USE MODD_DRAGTREE_n, ONLY : &
+ LDRAGTREE_n => LDRAGTREE, &
+ LDEPOTREE_n => LDEPOTREE, &
+ XVDEPOTREE_n => XVDEPOTREE
+!
+!-----------------------------------------------------------------------------
+!
+!* 0. DECLARATIONS
+! -----------------
+IMPLICIT NONE
+!
+LOGICAL, SAVE :: LDRAGTREE
+LOGICAL, SAVE :: LDEPOTREE
+REAL, SAVE :: XVDEPOTREE
+!
+NAMELIST /NAM_DRAGTREEn/LDRAGTREE,LDEPOTREE,XVDEPOTREE
+!
+CONTAINS
+!
+SUBROUTINE INIT_NAM_DRAGTREEn
+ LDRAGTREE = LDRAGTREE_n
+ LDEPOTREE = LDEPOTREE_n
+ XVDEPOTREE = XVDEPOTREE_n
+END SUBROUTINE INIT_NAM_DRAGTREEn
+!
+SUBROUTINE UPDATE_NAM_DRAGTREEn
+ LDRAGTREE_n = LDRAGTREE
+ LDEPOTREE_n = LDEPOTREE
+ XVDEPOTREE_n = XVDEPOTREE
+END SUBROUTINE UPDATE_NAM_DRAGTREEn
+!
+END MODULE MODN_DRAGTREE_n
diff --git a/src/MNH/phys_paramn.f90 b/src/MNH/phys_paramn.f90
index 6cea6a86e..9854052e4 100644
--- a/src/MNH/phys_paramn.f90
+++ b/src/MNH/phys_paramn.f90
@@ -233,6 +233,7 @@ END MODULE MODI_PHYS_PARAM_n
! P. Wautelet 28/03/2018: replace TEMPORAL_DIST by DATETIME_DISTANCE
! P. Wautelet 05/2016-04/2018: new data structures and calls for I/O
! P. Wautelet 21/11/2019: ZRG_HOUR and ZRAT_HOUR are now parameter arrays
+!! 11/2019 C.Lac correction in the drag formula and application to building in addition to tree
!!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -311,6 +312,7 @@ USE MODI_SEDIM_SALT
USE MODI_DUST_FILTER
USE MODI_SALT_FILTER
USE MODI_DRAG_VEG
+USE MODI_DRAG_BLD
USE MODD_DUST
USE MODD_SALT
USE MODD_PASPOL
@@ -321,7 +323,8 @@ USE MODE_SALT_PSD
USE MODE_AERO_PSD
USE MODE_MNH_TIMING
USE MODD_TURB_FLUX_AIRCRAFT_BALLOON, ONLY : XTHW_FLUX, XRCW_FLUX, XSVW_FLUX
-USE MODD_DRAGTREE
+USE MODD_DRAGTREE_n
+USE MODD_DRAGBLDG_n
!
USE MODD_TIME, ONLY : TDTEXP ! Ajout PP
USE MODI_AEROZON ! Ajout PP
@@ -737,6 +740,23 @@ CALL SUNPOS_n ( XZENITH, ZCOSZEN, ZSINZEN, ZAZIMSOL )
WRITE(UNIT=ILUOUT,FMT='(" RADIATIONS called for KTCOUNT=",I6, &
& "with the CLOUD_ONLY option set ",L2)') KTCOUNT,OCLOUD_ONLY
!
+ !
+ WHERE (XDIRFLASWD.LT.0.0)
+ XDIRFLASWD=0.0
+ ENDWHERE
+ !
+ WHERE (XDIRFLASWD.GT.1500.0)
+ XDIRFLASWD=1500.0
+ ENDWHERE
+ !
+ WHERE (XSCAFLASWD.LT.0.0)
+ XSCAFLASWD=0.0
+ ENDWHERE
+ !
+ WHERE (XSCAFLASWD.GT.1500.0)
+ XSCAFLASWD=1500.0
+ ENDWHERE
+ !
WHERE( XDIRFLASWD(:,:,1) + XSCAFLASWD(:,:,1) >0. )
XALBUV(:,:) = ( XDIR_ALB(:,:,1) * XDIRFLASWD(:,:,1) &
+ XSCA_ALB(:,:,1) * XSCAFLASWD(:,:,1) ) &
@@ -1234,10 +1254,11 @@ ZTIME1 = ZTIME2
XTIME_BU_PROCESS = 0.
XTIME_LES_BU_PROCESS = 0.
!
-IF (LDRAGTREE) CALL DRAG_VEG(XTSTEP,XUT,XVT,XTKET,LDEPOTREE,XVDEPOTREE, &
- CCLOUD, XPABST,XTHT,XRT,XSVT, &
- XRHODJ,XZZ,XRUS, XRVS, &
- XRTKES,XRTHS, XRRS,XRSVS)
+IF (LDRAGTREE) CALL DRAG_VEG( XTSTEP, XUT, XVT, XTKET, LDEPOTREE, XVDEPOTREE, &
+ CCLOUD, XPABST, XTHT, XRT, XSVT, XRHODJ, XZZ, &
+ XRUS, XRVS, XRTKES, XRTHS, XRRS, XRSVS )
+!
+IF (LDRAGBLDG) CALL DRAG_BLD ( XTSTEP, XUT, XVT, XTKET, XRHODJ, XZZ, XRUS, XRVS, XRTKES )
!
CALL SECOND_MNH2(ZTIME2)
!
@@ -1299,6 +1320,7 @@ IF ( CTURB == 'TKEL' ) THEN
ENDIF
ZSFCO2(IIB-1,:)=ZSFCO2(IIB,:)
END IF
+ !
IF ( CLBCX(2) /= "CYCL" .AND. LEAST_ll()) THEN
ZSFTH(IIE+1,:)=ZSFTH(IIE,:)
ZSFRV(IIE+1,:)=ZSFRV(IIE,:)
@@ -1312,6 +1334,7 @@ IF ( CTURB == 'TKEL' ) THEN
ENDIF
ZSFCO2(IIE+1,:)=ZSFCO2(IIE,:)
END IF
+ !
IF ( CLBCY(1) /= "CYCL" .AND. LSOUTH_ll()) THEN
ZSFTH(:,IJB-1)=ZSFTH(:,IJB)
ZSFRV(:,IJB-1)=ZSFRV(:,IJB)
@@ -1325,6 +1348,7 @@ IF ( CTURB == 'TKEL' ) THEN
ENDIF
ZSFCO2(:,IJB-1)=ZSFCO2(:,IJB)
END IF
+ !
IF ( CLBCY(2) /= "CYCL" .AND. LNORTH_ll()) THEN
ZSFTH(:,IJE+1)=ZSFTH(:,IJE)
ZSFRV(:,IJE+1)=ZSFRV(:,IJE)
diff --git a/src/MNH/read_exsegn.f90 b/src/MNH/read_exsegn.f90
index d09c3323b..dbff7eddb 100644
--- a/src/MNH/read_exsegn.f90
+++ b/src/MNH/read_exsegn.f90
@@ -292,6 +292,7 @@ END MODULE MODI_READ_EXSEG_n
!! Modification 01/2019 (P. Wautelet) bugs correction: incorrect writes
!! Modification 01/2019 (R. Honnert) remove SURF in CMF_UPDRAFT
!! Bielli S. 02/2019 Sea salt : significant sea wave height influences salt emission; 5 salt modes
+!! C.Lac 11/2019 correction in the drag formula and application to building in addition to tree
!!------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -343,7 +344,8 @@ USE MODN_SERIES_n
USE MODN_TURB_CLOUD
USE MODN_TURB
USE MODN_MEAN
-USE MODN_DRAGTREE
+USE MODN_DRAGTREE_n
+USE MODN_DRAGBLDG_n
USE MODN_LATZ_EDFLX
!
USE MODD_NSV,NSV_USER_n=>NSV_USER
@@ -461,6 +463,8 @@ CCPLFILE(:)=" "
CALL INIT_NAM_CONFN
CALL INIT_NAM_DYNN
CALL INIT_NAM_ADVN
+CALL INIT_NAM_DRAGTREEN
+CALL INIT_NAM_DRAGBLDGN
CALL INIT_NAM_PARAMN
CALL INIT_NAM_PARAM_RADN
#ifdef MNH_ECRAD
@@ -514,6 +518,10 @@ CALL POSNAM(ILUSEG,'NAM_SERIESN',GFOUND,ILUOUT)
IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_SERIESn)
CALL POSNAM(ILUSEG,'NAM_BLOWSNOWN',GFOUND,ILUOUT)
IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_BLOWSNOWn)
+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)
!
IF (KMI == 1) THEN
WRITE(UNIT=ILUOUT,FMT="(' namelists common to all the models ')")
@@ -650,8 +658,6 @@ IF (KMI == 1) THEN
#endif
CALL POSNAM(ILUSEG,'NAM_CONDSAMP',GFOUND,ILUOUT)
IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_CONDSAMP)
- CALL POSNAM(ILUSEG,'NAM_DRAGTREE',GFOUND,ILUOUT)
- IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_DRAGTREE)
CALL POSNAM(ILUSEG,'NAM_2D_FRC',GFOUND,ILUOUT)
IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_2D_FRC)
CALL POSNAM(ILUSEG,'NAM_LATZ_EDFLX',GFOUND)
@@ -2815,6 +2821,8 @@ END IF
!
CALL UPDATE_NAM_LUNITN
CALL UPDATE_NAM_CONFN
+CALL UPDATE_NAM_DRAGTREEN
+CALL UPDATE_NAM_DRAGBLDGN
CALL UPDATE_NAM_DYNN
CALL UPDATE_NAM_ADVN
CALL UPDATE_NAM_PARAMN
diff --git a/src/SURFEX/allocate_physio.F90 b/src/SURFEX/allocate_physio.F90
index 1028fa249..745a7d007 100644
--- a/src/SURFEX/allocate_physio.F90
+++ b/src/SURFEX/allocate_physio.F90
@@ -33,6 +33,7 @@
!! -------------
!! Original xx/xxxx
!! Modified 10/2014 P. Samuelsson MEB
+!! 11/2019 C.Lac correction in the drag formula and application to building in addition to tree
!
!
USE MODD_ISBA_OPTIONS_n, ONLY : ISBA_OPTIONS_t
@@ -99,11 +100,7 @@ ELSE
ENDIF
! - vegetation: Ags parameters ('AGS', 'LAI', 'AST', 'LST', 'NIT' options)
!
-IF (IO%CPHOTO/='NON'.OR.LTREEDRAG) THEN
- ALLOCATE(PK%XH_TREE (ISIZE ))
-ELSE
- ALLOCATE(PK%XH_TREE (0 ))
-ENDIF
+ALLOCATE(PK%XH_TREE (ISIZE ))
!
IF (IO%CPHOTO/='NON') THEN
ALLOCATE(PK%XRE25 (ISIZE ))
diff --git a/src/SURFEX/convert_patch_isba.F90 b/src/SURFEX/convert_patch_isba.F90
index cad666663..b04278c20 100644
--- a/src/SURFEX/convert_patch_isba.F90
+++ b/src/SURFEX/convert_patch_isba.F90
@@ -46,6 +46,7 @@
!! coupled to atmosphere)
!! P Samuelsson 10/2014 MEB
! P. Wautelet 15/02/2019: bugfix: allocate ZSTRESS only when its size has a meaning
+!! 11/2019 C.Lac correction in the drag formula and application to building in addition to tree
!
!----------------------------------------------------------------------------
!
@@ -215,7 +216,6 @@ IF (OFIX) THEN
PK%NR_P,IO%NPATCH,KPATCH,KDECADE=KDEC)
ENDIF
!
- IF (IO%CPHOTO/='NON'.OR.LTREEDRAG) THEN
IF (GDATA .AND. ANY(DTV%LDATA_H_TREE)) THEN
CALL AV_PGD_PARAM(DTV%XPAR_LAI, DTV%XPAR_VEG, &
PK%XH_TREE,DTV%XPAR_VEGTYPE,DTV%XPAR_H_TREE,YTREE,'ARI',PK%NR_P,IO%NPATCH,KPATCH)
@@ -223,7 +223,6 @@ IF (OFIX) THEN
CALL AV_PGD_1P(DTCO, PK%XH_TREE,PCOVER,XDATA_H_TREE(:,:),YTREE,'ARI',OCOVER,&
PK%NR_P,IO%NPATCH,KPATCH,KDECADE=KDEC)
ENDIF
- ENDIF
!
IF (IO%CPHOTO/='NON') THEN
!
diff --git a/src/SURFEX/get_surf_varn.F90 b/src/SURFEX/get_surf_varn.F90
index f732f0316..de14d16e0 100644
--- a/src/SURFEX/get_surf_varn.F90
+++ b/src/SURFEX/get_surf_varn.F90
@@ -10,7 +10,8 @@
PZ0H_WATER, PZ0H_NATURE, PZ0H_TOWN, PQS_SEA, &
PQS_WATER, PQS_NATURE, PQS_TOWN, PPSNG, PPSNV, &
PZS, PSERIES, PTWSNOW, PSSO_STDEV, PLON, PLAT, &
- PBARE, PLAI_TREE, PH_TREE )
+ PBARE, PLAI_TREE, PH_TREE, &
+ PWALL_O_HOR, PBUILD_HEIGHT )
! #######################################################################
!
!!**** *GET_SURF_VAR_n* - gets some surface fields on atmospheric grid
@@ -47,6 +48,7 @@
! S. Riette 06/2010 PSSO_STDEV and PTWSNOW added
! B. Decharme 09/2012 Argument added in GET_FLUX_n
! B. Decharme 05/2013 Argument added in GET_FLUX_n for debug in ARP/AL/AR
+!! C. Lac 11/2019 correction in the drag formula and application to building in addition to tree
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -144,8 +146,11 @@ REAL, DIMENSION(:), INTENT(OUT), OPTIONAL :: PLON ! longitude
REAL, DIMENSION(:), INTENT(OUT), OPTIONAL :: PLAT ! latitude
!
REAL, DIMENSION(:), INTENT(OUT), OPTIONAL :: PBARE ! bare soil fraction on grid mesh (-)
-REAL, DIMENSION(:), INTENT(OUT), OPTIONAL :: PLAI_TREE ! Leaf Area Index on grid mesh (-)
-REAL, DIMENSION(:), INTENT(OUT), OPTIONAL :: PH_TREE ! Height of trees on grid mesh (-)
+REAL, DIMENSION(:), INTENT(OUT), OPTIONAL :: PLAI_TREE ! Leaf Area Index on grid mesh (-)
+REAL, DIMENSION(:), INTENT(OUT), OPTIONAL :: PH_TREE ! Height of trees on grid mesh (-)
+!
+REAL, DIMENSION(:), INTENT(OUT), OPTIONAL :: PWALL_O_HOR ! Facade area density on grid mesh [m^2(fac.)/m^2(town)]
+REAL, DIMENSION(:), INTENT(OUT), OPTIONAL :: PBUILD_HEIGHT ! Building height on grid mesh [m]
!
!-------------------------------------------------------------------------------
!
@@ -425,7 +430,8 @@ ENDIF
!
!-------------------------------------------------------------------------------
!
-IF ( PRESENT(PQS_TOWN) .OR. PRESENT(PZ0_TOWN) .OR. PRESENT(PZ0H_TOWN) ) THEN
+IF ( PRESENT(PQS_TOWN) .OR. PRESENT(PZ0_TOWN) .OR. PRESENT(PZ0H_TOWN) .OR. &
+ PRESENT(PWALL_O_HOR) .OR. PRESENT (PBUILD_HEIGHT) ) THEN
!
! Get parameters over town tile
!
@@ -440,8 +446,9 @@ IF ( PRESENT(PQS_TOWN) .OR. PRESENT(PZ0_TOWN) .OR. PRESENT(PZ0H_TOWN) ) THEN
IMASK(:)=0
CALL GET_1D_MASK(KI_TOWN, KI, PTOWN, IMASK(1:KI_TOWN))
!
- CALL GET_VAR_TOWN_n(TM%TD%O, TM%TD%D, HPROGRAM, KI_TOWN, &
- ZFIELD1(1:KI_TOWN), ZFIELD2(1:KI_TOWN), ZFIELD3(1:KI_TOWN))
+ CALL GET_VAR_TOWN_n(TM%TOP, TM%TD%O, TM%TD%D, TM%NT,HPROGRAM, KI_TOWN, &
+ ZFIELD1(1:KI_TOWN), ZFIELD2(1:KI_TOWN), ZFIELD3(1:KI_TOWN),&
+ ZFIELD4(1:KI_TOWN), ZFIELD5(1:KI_TOWN) )
!
IF(PRESENT(PQS_TOWN))THEN
PQS_TOWN (:) = XUNDEF
@@ -464,6 +471,20 @@ IF ( PRESENT(PQS_TOWN) .OR. PRESENT(PZ0_TOWN) .OR. PRESENT(PZ0H_TOWN) ) THEN
END DO
ENDIF
!
+ IF(PRESENT(PWALL_O_HOR))THEN
+ PWALL_O_HOR (:) = XUNDEF
+ DO JI = 1, KI_TOWN
+ PWALL_O_HOR(IMASK(JI)) = ZFIELD4(JI)
+ END DO
+ ENDIF
+ !
+ IF(PRESENT(PBUILD_HEIGHT))THEN
+ PBUILD_HEIGHT (:) = XUNDEF
+ DO JI = 1, KI_TOWN
+ PBUILD_HEIGHT(IMASK(JI)) = ZFIELD5(JI)
+ END DO
+ ENDIF
+ !
END IF
!
!* 5. Orography
diff --git a/src/SURFEX/get_var_townn.F90 b/src/SURFEX/get_var_townn.F90
index 0f2e9c9d7..382f32216 100644
--- a/src/SURFEX/get_var_townn.F90
+++ b/src/SURFEX/get_var_townn.F90
@@ -3,7 +3,8 @@
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!SFX_LIC for details. version 1.
! #########
- SUBROUTINE GET_VAR_TOWN_n (DGO, D, HPROGRAM,KI,PQS,PZ0,PZ0H)
+ SUBROUTINE GET_VAR_TOWN_n (TOP, DGO, D, NT, HPROGRAM,KI,PQS,PZ0,PZ0H, &
+ PWALL_O_HOR,PBUILD_HEIGHT )
! ###################################################
!
!!**** *GET_VAR_TOWN_n* - routine to get variables defined only over town
@@ -33,6 +34,7 @@
!! -------------
!! Original 02/2006
! M. Jidane 08/2008 Z0 and Z0H recovery from town tiles
+!! C.Lac 11/2019 correction in the drag formula and application to building in addition to tree
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -40,6 +42,8 @@
!
!
USE MODD_DIAG_n, ONLY : DIAG_t, DIAG_OPTIONS_t
+USE MODD_TEB_n, ONLY : TEB_NP_t
+USE MODD_TEB_OPTION_n, ONLY : TEB_OPTIONS_t
!
USE MODI_GET_LUOUT
USE MODD_SURF_PAR, ONLY : XUNDEF
@@ -54,19 +58,24 @@ IMPLICIT NONE
!* 0.1 Declarations of arguments
! -------------------------
!
+TYPE(TEB_OPTIONS_t), INTENT(IN) :: TOP
TYPE(DIAG_OPTIONS_t), INTENT(IN) :: DGO
-TYPE(DIAG_t), INTENT(IN) :: D
+TYPE(DIAG_t), INTENT(IN) :: D
+TYPE(TEB_NP_t), INTENT(IN) :: NT
!
CHARACTER(LEN=6), INTENT(IN) :: HPROGRAM
INTEGER, INTENT(IN) :: KI ! Number of points
REAL, DIMENSION(KI), INTENT(OUT) :: PQS ! surface humidity
REAL, DIMENSION(KI), INTENT(OUT) :: PZ0 ! surface roughness length
REAL, DIMENSION(KI), INTENT(OUT) :: PZ0H ! surface roughness length for heat
+REAL, DIMENSION(KI), INTENT(OUT) :: PWALL_O_HOR ! Facade surface density [m^2(fac.)/m^2(town)]
+REAL, DIMENSION(KI), INTENT(OUT) :: PBUILD_HEIGHT ! Building height [m]
!
!
!* 0.2 Declarations of local variables
! -------------------------------
!
+INTEGER :: JP ! loop counter on TEB patches
INTEGER :: ILUOUT
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
@@ -76,17 +85,28 @@ IF (LHOOK) CALL DR_HOOK('GET_VAR_TOWN_N',0,ZHOOK_HANDLE)
!-------------------------------------------------------------------------------
!
IF (DGO%LSURF_VARS) THEN
- PQS = D%XQS
- ELSE
- PQS = XUNDEF
+ PQS = D%XQS
+ELSE
+ PQS = XUNDEF
ENDIF
IF (DGO%LCOEF) THEN
- PZ0 = D%XZ0
- PZ0H = D%XZ0H
- ELSE
- PZ0 = XUNDEF
- PZ0H = XUNDEF
-ENDIF
+ PZ0 = D%XZ0
+ PZ0H = D%XZ0H
+ELSE
+ PZ0 = XUNDEF
+ PZ0H = XUNDEF
+ENDIF
+!
+!* building height and external wall coverage fraction
+!
+PWALL_O_HOR = 0.
+PBUILD_HEIGHT = 0.
+
+DO JP=1,TOP%NTEB_PATCH
+ PWALL_O_HOR (:) = PWALL_O_HOR (:) + TOP%XTEB_PATCH(:,JP) * NT%AL(JP)%XWALL_O_HOR(:)
+ PBUILD_HEIGHT(:) = PBUILD_HEIGHT(:) + TOP%XTEB_PATCH(:,JP) * NT%AL(JP)%XBLD_HEIGHT(:)
+END DO
+
IF (LHOOK) CALL DR_HOOK('GET_VAR_TOWN_N',1,ZHOOK_HANDLE)
!
!==============================================================================
diff --git a/src/SURFEX/get_vegn.F90 b/src/SURFEX/get_vegn.F90
index a41238194..d1e9c583d 100644
--- a/src/SURFEX/get_vegn.F90
+++ b/src/SURFEX/get_vegn.F90
@@ -29,6 +29,7 @@
!! MODIFICATIONS
!! -------------
!! Original 07/2009
+!! 11/2019 C.Lac correction in the drag formula and application to building in addition to tree
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -100,16 +101,16 @@ IPATCH_BONE = VEGTYPE_TO_PATCH(NVT_BONE, IO%NPATCH)
IPATCH_BOND = VEGTYPE_TO_PATCH(NVT_BOND, IO%NPATCH)
-!ZWORK(:) = S%XVEGTYPE(:,NVT_TRBE) + S%XVEGTYPE(:,NVT_TRBD) + S%XVEGTYPE(:,NVT_TEBE) + &
-! S%XVEGTYPE(:,NVT_TEBD) + S%XVEGTYPE(:,NVT_TENE) + S%XVEGTYPE(:,NVT_BOBD) + &
-! S%XVEGTYPE(:,NVT_BONE) + S%XVEGTYPE(:,NVT_BOND)
-
ZH_TREE(:) = 0.
ZLAI(:) = 0.
ZWORK(:) = 0.
!
DO JP = 1,IO%NPATCH
+ PK => NP%AL(JP)
+ PEK => NPE%AL(JP)
+END DO
!
+DO JP = 1,IO%NPATCH
IF (JP==IPATCH_TRBE .OR. JP==IPATCH_TRBD .OR. JP==IPATCH_TEBE .OR. JP==IPATCH_TEBD .OR. &
JP==IPATCH_TENE .OR. JP==IPATCH_BOBD .OR. JP==IPATCH_BONE .OR. JP==IPATCH_BOND) THEN
!
--
GitLab