From e7c040a6220ced53766c51f489c8a07c5df85249 Mon Sep 17 00:00:00 2001
From: Gaelle DELAUTIER <gaelle.delautier@meteo.fr>
Date: Thu, 26 Apr 2018 11:29:10 +0200
Subject: [PATCH] V.Vionnet + Gaelle (24/04/2018) : add blowing snow scheme
---
src/MNH/advection_metsv.f90 | 97 ++++++++++++++++++-
src/MNH/boundaries.f90 | 122 +++++++++++++++++++++++-
src/MNH/default_desfmn.f90 | 22 ++++-
src/MNH/diag.f90 | 4 +-
src/MNH/endstep.f90 | 26 ++++-
src/MNH/exchange.f90 | 15 +++
src/MNH/goto_model_wrapper.f90 | 3 +
src/MNH/ground_paramn.f90 | 98 +++++++++++++++++--
src/MNH/ini_dynamics.f90 | 18 ++--
src/MNH/ini_micron.f90 | 16 ++++
src/MNH/ini_modeln.f90 | 15 ++-
src/MNH/ini_nsv.f90 | 26 ++++-
src/MNH/ini_segn.f90 | 6 +-
src/MNH/ini_spectren.f90 | 4 +-
src/MNH/init_ground_paramn.f90 | 26 ++++-
src/MNH/initial_guess.f90 | 8 ++
src/MNH/modd_diag_flag.f90 | 1 +
src/MNH/modd_dynn.f90 | 4 +
src/MNH/modd_nsv.f90 | 9 ++
src/MNH/modeln.f90 | 29 +++++-
src/MNH/modn_dynn.f90 | 6 +-
src/MNH/read_desfmn.f90 | 27 +++++-
src/MNH/read_exsegn.f90 | 59 +++++++++++-
src/MNH/read_field.f90 | 47 +++++++++
src/MNH/relaxation.f90 | 15 ++-
src/MNH/relaxdef.f90 | 11 ++-
src/MNH/turb_hor_sv_corr.f90 | 22 +++--
src/MNH/turb_hor_sv_flux.f90 | 17 +++-
src/MNH/turb_ver_sv_corr.f90 | 17 +++-
src/MNH/turb_ver_sv_flux.f90 | 14 ++-
src/MNH/update_nsv.f90 | 4 +
src/MNH/write_desfmn.f90 | 18 +++-
src/MNH/write_lfifm1_for_diag.f90 | 115 ++++++++++++++++++++++
src/MNH/write_lfin.f90 | 72 +++++++++++++-
src/SURFEX/canopy_grid_update.F90 | 59 ++++++++----
src/SURFEX/compute_isba_parameters.F90 | 36 ++++++-
src/SURFEX/coupling_isba_canopyn.F90 | 102 ++++++++++++++++++--
src/SURFEX/coupling_isba_orographyn.F90 | 6 +-
src/SURFEX/coupling_isba_svatn.F90 | 6 +-
src/SURFEX/coupling_isban.F90 | 34 ++++++-
src/SURFEX/coupling_naturen.F90 | 10 +-
src/SURFEX/coupling_surf_atmn.F90 | 5 +-
src/SURFEX/coupling_tsz0n.F90 | 6 +-
src/SURFEX/init_chemicaln.F90 | 23 ++++-
src/SURFEX/init_flaken.F90 | 3 +-
src/SURFEX/init_isban.F90 | 10 +-
src/SURFEX/init_naturen.F90 | 6 +-
src/SURFEX/init_seafluxn.F90 | 3 +-
src/SURFEX/init_surf_atmn.F90 | 2 +-
src/SURFEX/init_watfluxn.F90 | 3 +-
src/SURFEX/isba.F90 | 18 +++-
src/SURFEX/isba_meb.F90 | 15 ++-
src/SURFEX/modd_canopyn.F90 | 3 +
src/SURFEX/modd_ch_flaken.F90 | 3 +-
src/SURFEX/modd_ch_isban.F90 | 2 +
src/SURFEX/modd_ch_seafluxn.F90 | 4 +-
src/SURFEX/modd_ch_watfluxn.F90 | 2 +
src/SURFEX/modd_surfexn.F90 | 2 +
src/SURFEX/modd_svn.F90 | 13 ++-
src/SURFEX/read_all_namelists.F90 | 2 +
src/SURFEX/read_sbln.F90 | 40 +++++++-
src/SURFEX/snow3L_isba.F90 | 118 +++++++++++++++++++++--
src/SURFEX/snowcro.F90 | 54 +++++++++--
src/SURFEX/write_diag_isban.F90 | 6 +-
src/SURFEX/write_diag_naturen.F90 | 6 +-
src/SURFEX/write_diag_seb_isban.F90 | 94 +++++++++++++++++-
src/SURFEX/write_diag_surf_atmn.F90 | 2 +-
src/SURFEX/write_isban.F90 | 2 +-
src/SURFEX/writesurf_sbln.F90 | 20 +++-
69 files changed, 1566 insertions(+), 147 deletions(-)
diff --git a/src/MNH/advection_metsv.f90 b/src/MNH/advection_metsv.f90
index 80d4edff8..b125cb434 100644
--- a/src/MNH/advection_metsv.f90
+++ b/src/MNH/advection_metsv.f90
@@ -134,6 +134,8 @@ END MODULE MODI_ADVECTION_METSV
!! 10/2016 (C.Lac) Correction on the flag for Strang splitting
!! to insure reproducibility between START and RESTA
!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!! 07/2017 (V. Vionnet) : add advection of 2D variables at
+!! the surface for the blowing snow scheme
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -147,6 +149,9 @@ USE MODD_IO_ll, ONLY: TFILEDATA
USE MODD_LUNIT_n, ONLY: TLUOUT
USE MODD_PARAM_n
USE MODD_TYPE_DATE, ONLY: DATE_TIME
+USE MODD_BLOWSNOW
+USE MODD_BLOWSNOW_n
+USE MODD_PARAMETERS
!
USE MODE_FIELD, ONLY: TFIELDDATA, TYPEREAL
USE MODE_FMWRIT
@@ -233,12 +238,17 @@ REAL, DIMENSION(SIZE(PTHT,1), SIZE(PTHT,2), SIZE(PTHT,3) ) :: ZRTHS_PPM
REAL, DIMENSION(SIZE(PTKET,1),SIZE(PTKET,2),SIZE(PTKET,3)) :: ZRTKES_PPM
REAL, DIMENSION(SIZE(PRT,1), SIZE(PRT,2), SIZE(PRT,3), SIZE(PRT,4) ) :: ZR
REAL, DIMENSION(SIZE(PSVT,1),SIZE(PSVT,2),SIZE(PSVT,3),SIZE(PSVT,4)) :: ZSV
+REAL, DIMENSION(SIZE(PSVT,1),SIZE(PSVT,2),SIZE(PSVT,3), NBLOWSNOW_2D) :: ZSNWC
+REAL, DIMENSION(SIZE(PSVT,1),SIZE(PSVT,2),SIZE(PSVT,3), NBLOWSNOW_2D) :: ZSNWC_INIT
+REAL, DIMENSION(SIZE(PSVT,1),SIZE(PSVT,2),SIZE(PSVT,3), NBLOWSNOW_2D) :: ZRSNWCS
! Guess at the sub time step
REAL, DIMENSION(SIZE(PRT,1), SIZE(PRT,2), SIZE(PRT,3), SIZE(PRT,4) ) :: ZRRS_OTHER
REAL, DIMENSION(SIZE(PSVT,1),SIZE(PSVT,2),SIZE(PSVT,3),SIZE(PSVT,4)) :: ZRSVS_OTHER
+REAL, DIMENSION(SIZE(PSVT,1),SIZE(PSVT,2),SIZE(PSVT,3),NBLOWSNOW_2D) :: ZRSNWCS_OTHER
! Tendencies since the beginning of the time step
REAL, DIMENSION(SIZE(PRT,1), SIZE(PRT,2), SIZE(PRT,3), SIZE(PRT,4) ) :: ZRRS_PPM
REAL, DIMENSION(SIZE(PSVT,1),SIZE(PSVT,2),SIZE(PSVT,3),SIZE(PSVT,4)) :: ZRSVS_PPM
+REAL, DIMENSION(SIZE(PSVT,1),SIZE(PSVT,2),SIZE(PSVT,3),NBLOWSNOW_2D) :: ZRSNWCS_PPM
! Guess at the end of the sub time step
REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZRHOX1,ZRHOX2
REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZRHOY1,ZRHOY2
@@ -246,7 +256,7 @@ REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZRHOZ1,ZRHOZ2
REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)):: ZT,ZEXN,ZLV,ZLS,ZCPH
! Temporary advected rhodj for PPM routines
!
-INTEGER :: JS,JR,JSV,JSPL ! Loop index
+INTEGER :: JS,JR,JSV,JSPL, JI, JJ ! Loop index
REAL :: ZTSTEP_PPM ! Sub Time step
LOGICAL :: GTKE
!
@@ -258,7 +268,7 @@ TYPE(LIST_ll), POINTER :: TZFIELDS1_ll ! list of fields to exchange
INTEGER :: IRESP ! Return code of FM routines
INTEGER :: ILUOUT ! logical unit
INTEGER :: ISPLIT_PPM ! temporal time splitting
-INTEGER :: IIB, IIE, IJB, IJE
+INTEGER :: IIB, IIE, IJB, IJE,IKB,IKE
TYPE(TFIELDDATA) :: TZFIELD
!-------------------------------------------------------------------------------
!
@@ -268,9 +278,24 @@ TYPE(TFIELDDATA) :: TZFIELD
ILUOUT = TLUOUT%NLU
!
CALL GET_INDICE_ll(IIB,IJB,IIE,IJE)
+IKB=1+JPVEXT
+IKE=SIZE(PSVT,3) - JPVEXT
+
!
GTKE=(SIZE(PTKET)/=0)
!
+!
+IF(LBLOWSNOW) THEN ! Put 2D Canopy blowing snow variables into a 3D array for advection
+ ZSNWC_INIT = 0.
+ ZRSNWCS = 0.
+
+ DO JSV=1,(NBLOWSNOW_2D)
+ ZSNWC_INIT(:,:,IKB,JSV) = XSNWCANO(:,:,JSV)
+ ZRSNWCS(:,:,IKB,JSV) = XRSNWCANOS(:,:,JSV)
+ END DO
+ENDIF
+!
+!
!-------------------------------------------------------------------------------
!
!* 2. COMPUTES THE CONTRAVARIANT COMPONENTS (FOR PPM ONLY)
@@ -442,6 +467,11 @@ END DO
DO JSV = 1, KSV
ZRSVS_OTHER(:,:,:,JSV) = PRSVS(:,:,:,JSV) - PSVT(:,:,:,JSV) * PRHODJ / PTSTEP
END DO
+IF(LBLOWSNOW) THEN
+ DO JSV = 1, (NBLOWSNOW_2D)
+ ZRSNWCS_OTHER(:,:,:,JSV) = ZRSNWCS(:,:,:,JSV) - ZSNWC_INIT(:,:,:,JSV) * PRHODJ / PTSTEP
+ END DO
+ENDIF
!
! Top and bottom Boundaries
!
@@ -453,6 +483,11 @@ END DO
DO JSV = 1, KSV
CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZRSVS_OTHER(:,:,:,JSV))
END DO
+IF(LBLOWSNOW) THEN
+ DO JSV = 1, (NBLOWSNOW_2D)
+ CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZRSNWCS_OTHER(:,:,:,JSV))
+ END DO
+END IF
!
! Exchanges on processors
!
@@ -466,6 +501,11 @@ NULLIFY(TZFIELDS0_ll)
DO JSV=1,KSV
CALL ADD3DFIELD_ll(TZFIELDS0_ll, ZRSVS_OTHER(:,:,:,JSV))
END DO
+ IF(LBLOWSNOW) THEN
+ DO JSV = 1, (NBLOWSNOW_2D)
+ CALL ADD3DFIELD_ll(TZFIELDS0_ll, ZRSNWCS_OTHER(:,:,:,JSV))
+ END DO
+ END IF
CALL UPDATE_HALO_ll(TZFIELDS0_ll,IINFO_ll)
CALL CLEANLIST_ll(TZFIELDS0_ll)
!!$END IF
@@ -486,6 +526,13 @@ ZTH = PTHT
ZTKE = PTKET
IF (KRR /=0 ) ZR = PRT
IF (KSV /=0 ) ZSV = PSVT
+IF(LBLOWSNOW) THEN
+ DO JSV = 1, (NBLOWSNOW_2D)
+ ZSNWC(:,:,:,JSV) = ZRSNWCS(:,:,:,JSV)* PTSTEP/ PRHODJ
+ CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZSNWC(:,:,:,JSV))
+ END DO
+ ZSNWC_INIT=ZSNWC
+ENDIF
!
IF (GTKE) PRTKES_ADV(:,:,:) = 0.
!
@@ -540,6 +587,28 @@ DO JSPL=1,KSPLIT
DO JSV = 1, KSV
CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZSV(:,:,:,JSV), PSVT(:,:,:,JSV))
END DO
+
+ IF(LBLOWSNOW) THEN ! Advection of Canopy mass at the 1st atmospheric level
+ ZRSNWCS_PPM(:,:,:,:) = 0.
+ !
+
+ CALL PPM_SCALAR (HLBCX,HLBCY, NBLOWSNOW_2D, TPDTCUR, ZRUCPPM, ZRVCPPM, ZRWCPPM,PTSTEP, &
+ ZTSTEP_PPM, PRHODJ, ZRHOX1, ZRHOX2, ZRHOY1, ZRHOY2, ZRHOZ1, ZRHOZ2, &
+ ZSNWC, ZRSNWCS_PPM, HSV_ADV_SCHEME)
+
+
+! Tendencies of PPM
+ ZRSNWCS(:,:,:,:) = ZRSNWCS(:,:,:,:) + ZRSNWCS_PPM (:,:,:,:) / KSPLIT
+! Guesses of the field inside the time splitting loop
+ DO JSV = 1, ( NBLOWSNOW_2D)
+ ZSNWC(:,:,:,JSV) = ZSNWC(:,:,:,JSV) + ZRSNWCS_PPM(:,:,:,JSV)*ZTSTEP_PPM/ PRHODJ(:,:,:)
+ END DO
+
+! Top and bottom Boundaries and LBC for the guesses
+ DO JSV = 1, (NBLOWSNOW_2D)
+ CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZSNWC(:,:,:,JSV), ZSNWC_INIT(:,:,:,JSV))
+ END DO
+ END IF
!
! Exchanges fields between processors
!
@@ -553,6 +622,11 @@ DO JSPL=1,KSPLIT
DO JSV=1,KSV
CALL ADD3DFIELD_ll(TZFIELDS1_ll, ZSV(:,:,:,JSV))
END DO
+ IF(LBLOWSNOW) THEN
+ DO JSV=1,(NBLOWSNOW_2D)
+ CALL ADD3DFIELD_ll(TZFIELDS1_ll, ZSNWC(:,:,:,JSV))
+ END DO
+ END IF
CALL UPDATE_HALO_ll(TZFIELDS1_ll,IINFO_ll)
CALL CLEANLIST_ll(TZFIELDS1_ll)
!!$ END IF
@@ -572,6 +646,25 @@ IF (GTKE) THEN
PRTKES(:,:,:) = MAX (PRTKES(:,:,:) , XTKEMIN * PRHODJ(:,:,:) / PTSTEP )
END IF
!
+!
+!-------------------------------------------------------------------------------
+! Update tendency for cano variables : from 3D to 2D
+!
+IF(LBLOWSNOW) THEN
+
+ DO JSV=1,(NBLOWSNOW_2D)
+ DO JI=1,SIZE(PSVT,1)
+ DO JJ=1,SIZE(PSVT,2)
+ XRSNWCANOS(JI,JJ,JSV) = SUM(ZRSNWCS(JI,JJ,IKB:IKE,JSV))
+ END DO
+ END DO
+ END DO
+IF(LWEST_ll()) XRSNWCANOS(IIB,:,:) = ZRSNWCS(IIB,:,IKB,:)
+IF(LEAST_ll()) XRSNWCANOS(IIE,:,:) = ZRSNWCS(IIE,:,IKB,:)
+IF(LSOUTH_ll()) XRSNWCANOS(:,IJB,:) = ZRSNWCS(:,IJB,IKB,:)
+IF(LNORTH_ll()) XRSNWCANOS(:,IJE,:) = ZRSNWCS(:,IJE,IKB,:)
+
+END IF
!-------------------------------------------------------------------------------
!
!* 5. BUDGETS
diff --git a/src/MNH/boundaries.f90 b/src/MNH/boundaries.f90
index 46fa4bc36..52a91fb29 100644
--- a/src/MNH/boundaries.f90
+++ b/src/MNH/boundaries.f90
@@ -172,6 +172,7 @@ END MODULE MODI_BOUNDARIES
!! J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1
!! Redelsperger & Pianezze : 08/2015 : add XPOND coefficient
!! Modification 01/2016 (JP Pinty) Add LIMA that is LBC for CCN and IFN
+!! Modification 18/07/17 (Vionnet) Add blowing snow variables
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -187,7 +188,9 @@ USE MODD_SALT, ONLY : LSALT
USE MODD_PASPOL, ONLY : LPASPOL
USE MODD_CONDSAMP, ONLY : LCONDSAMP
USE MODD_ELEC_DESCR
-USE MODD_ELEC_n
+USE MODD_ELEC_n
+USE MODD_BLOWSNOW, ONLY : LBLOWSNOW,NBLOWSNOW_2D
+USE MODD_BLOWSNOW_n
USE MODD_REF_n
USE MODD_PARAM_n, ONLY : CELEC,CCLOUD
USE MODD_LBC_n, ONLY : XPOND
@@ -266,7 +269,7 @@ REAL :: ZPOND ! Coeff PONDERATION LS
INTEGER :: ILBX,ILBY ! size of LB fields' arrays
LOGICAL, SAVE, DIMENSION(:), ALLOCATABLE :: GCHBOUNDARY, GAERBOUNDARY,&
GDSTBOUNDARY, GSLTBOUNDARY, GPPBOUNDARY, &
- GCSBOUNDARY, GICBOUNDARY, GLIMABOUNDARY
+ GCSBOUNDARY, GICBOUNDARY, GLIMABOUNDARY,GSNWBOUNDARY
LOGICAL, SAVE :: GFIRSTCALL1 = .TRUE.
LOGICAL, SAVE :: GFIRSTCALL2 = .TRUE.
LOGICAL, SAVE :: GFIRSTCALL3 = .TRUE.
@@ -473,6 +476,7 @@ SELECT CASE ( HLBCX(1) )
IF(SIZE(PRT) /= 0) PRT (IIB-JEXT,:,:,:) = PRT (IIB-1+JEXT,:,:,:)
IF(SIZE(PSVT) /= 0) PSVT(IIB-JEXT,:,:,:) = PSVT (IIB-1+JEXT,:,:,:)
IF(SIZE(PSRCT) /= 0) PSRCT (IIB-JEXT,:,:) = PSRCT (IIB-1+JEXT,:,:)
+ IF(LBLOWSNOW) XSNWCANO(IIB-JEXT,:,:) = XSNWCANO(IIB-1+JEXT,:,:)
!
END DO
!
@@ -548,6 +552,30 @@ SELECT CASE ( HLBCX(1) )
END IF
!
END DO
+ !
+ IF(LBLOWSNOW) THEN
+ DO JSV=1 ,NBLOWSNOW_2D
+ WHERE ( PUT(IIB,:,IKB) <= 0. ) ! OUTFLOW condition
+ XSNWCANO(IIB-1,:,JSV) = MAX(0.,2.*XSNWCANO(IIB,:,JSV) - &
+ XSNWCANO(IIB+1,:,JSV))
+ ELSEWHERE ! INFLOW condition
+ XSNWCANO(IIB-1,:,JSV) = 0. ! Assume no snow enter throug
+ ! boundaries
+ END WHERE
+ END DO
+ DO JSV=NSV_SNWBEG ,NSV_SNWEND
+ IF(SIZE(PUT) /= 0) THEN
+ WHERE ( PUT(IIB,:,:) <= 0. ) ! OUTFLOW condition
+ PSVT(IIB-1,:,:,JSV) = MAX(0.,2.*PSVT(IIB,:,:,JSV) - &
+ PSVT(IIB+1,:,:,JSV))
+ ELSEWHERE ! INFLOW condition
+ PSVT(IIB-1,:,:,JSV) = 0. ! Assume no snow enter throug
+ ! boundaries
+ END WHERE
+ END IF
+ !
+ END DO
+ ENDIF
!
!
END SELECT
@@ -576,6 +604,7 @@ SELECT CASE ( HLBCX(2) )
IF(SIZE(PRT) /= 0) PRT (IIE+JEXT,:,:,:) = PRT (IIE+1-JEXT,:,:,:)
IF(SIZE(PSVT) /= 0) PSVT(IIE+JEXT,:,:,:) = PSVT (IIE+1-JEXT,:,:,:)
IF(SIZE(PSRCT) /= 0) PSRCT (IIE+JEXT,:,:)= PSRCT (IIE+1-JEXT,:,:)
+ IF(LBLOWSNOW) XSNWCANO(IIE+JEXT,:,:) = XSNWCANO(IIE+1-JEXT,:,:)
!
END DO
!
@@ -653,6 +682,29 @@ SELECT CASE ( HLBCX(2) )
!
END DO
!
+ IF(LBLOWSNOW) THEN
+ DO JSV=1 ,3
+ WHERE ( PUT(IIE+1,:,IKB) >= 0. ) ! OUTFLOW condition
+ XSNWCANO(IIE+1,:,JSV) = MAX(0.,2.*XSNWCANO(IIE,:,JSV) - &
+ XSNWCANO(IIE-1,:,JSV))
+ ELSEWHERE ! INFLOW condition
+ XSNWCANO(IIE+1,:,JSV) = 0. ! Assume no snow enter throug
+ ! boundaries
+ END WHERE
+ END DO
+ DO JSV=NSV_SNWBEG ,NSV_SNWEND
+ IF(SIZE(PUT) /= 0) THEN
+ WHERE ( PUT(IIE+1,:,:) >= 0. ) ! OUTFLOW condition
+ PSVT(IIE+1,:,:,JSV) = MAX(0.,2.*PSVT(IIE,:,:,JSV) - &
+ PSVT(IIE-1,:,:,JSV))
+ ELSEWHERE ! INFLOW condition
+ PSVT(IIE+1,:,:,JSV) = 0. ! Assume no snow enter throug
+ ! boundaries
+ END WHERE
+ END IF
+ !
+ END DO
+ END IF
!
END SELECT
!
@@ -679,6 +731,7 @@ SELECT CASE ( HLBCY(1) )
IF(SIZE(PRT) /= 0) PRT (:,IJB-JEXT,:,:) = PRT (:,IJB-1+JEXT,:,:)
IF(SIZE(PSVT) /= 0) PSVT (:,IJB-JEXT,:,:)= PSVT (:,IJB-1+JEXT,:,:)
IF(SIZE(PSRCT) /= 0) PSRCT(:,IJB-JEXT,:) = PSRCT(:,IJB-1+JEXT,:)
+ IF(LBLOWSNOW) XSNWCANO(:,IJB-JEXT,:) = XSNWCANO(:,IJB-1+JEXT,:)
!
END DO
!
@@ -752,6 +805,30 @@ SELECT CASE ( HLBCY(1) )
END IF
!
END DO
+!
+ IF(LBLOWSNOW) THEN
+ DO JSV=1 ,3
+ WHERE ( PVT(:,IJB,IKB) <= 0. ) ! OUTFLOW condition
+ XSNWCANO(:,IJB-1,JSV) = MAX(0.,2.*XSNWCANO(:,IJB,JSV) - &
+ XSNWCANO(:,IJB+1,JSV))
+ ELSEWHERE ! INFLOW condition
+ XSNWCANO(:,IJB-1,JSV) = 0. ! Assume no snow enter throug
+ ! boundaries
+ END WHERE
+ END DO
+ DO JSV=NSV_SNWBEG ,NSV_SNWEND
+ IF(SIZE(PVT) /= 0) THEN
+ WHERE ( PVT(:,IJB,:) <= 0. ) ! OUTFLOW condition
+ PSVT(:,IJB-1,:,JSV) = MAX(0.,2.*PSVT(:,IJB,:,JSV) - &
+ PSVT(:,IJB+1,:,JSV))
+ ELSEWHERE ! INFLOW condition
+ PSVT(:,IJB-1,:,JSV) = 0. ! Assume no snow enter throug
+ ! boundaries
+ END WHERE
+ END IF
+ !
+ END DO
+ END IF
!
!
END SELECT
@@ -780,6 +857,7 @@ SELECT CASE ( HLBCY(2) )
IF(SIZE(PRT) /= 0) PRT (:,IJE+JEXT,:,:) = PRT (:,IJE+1-JEXT,:,:)
IF(SIZE(PSVT) /= 0) PSVT (:,IJE+JEXT,:,:)= PSVT (:,IJE+1-JEXT,:,:)
IF(SIZE(PSRCT) /= 0) PSRCT(:,IJE+JEXT,:) = PSRCT(:,IJE+1-JEXT,:)
+ IF(LBLOWSNOW) XSNWCANO(:,IJE+JEXT,:) = XSNWCANO(:,IJE+1-JEXT,:)
!
END DO
!
@@ -858,6 +936,30 @@ SELECT CASE ( HLBCY(2) )
!
END DO
!
+ IF(LBLOWSNOW) THEN
+ DO JSV=1 ,3
+ WHERE ( PVT(:,IJE+1,IKB) >= 0. ) ! OUTFLOW condition
+ XSNWCANO(:,IJE+1,JSV) = MAX(0.,2.*XSNWCANO(:,IJE,JSV) - &
+ XSNWCANO(:,IJE-1,JSV))
+ ELSEWHERE ! INFLOW condition
+ XSNWCANO(:,IJE+1,JSV) = 0. ! Assume no snow enter throug
+ ! boundaries
+ END WHERE
+ END DO
+ DO JSV=NSV_SNWBEG ,NSV_SNWEND
+ !
+ IF(SIZE(PVT) /= 0) THEN
+ WHERE ( PVT(:,IJE+1,:) >= 0. ) ! OUTFLOW condition
+ PSVT(:,IJE+1,:,JSV) = MAX(0.,2.*PSVT(:,IJE,:,JSV) - &
+ PSVT(:,IJE-1,:,JSV))
+ ELSEWHERE ! INFLOW condition
+ PSVT(:,IJE+1,:,JSV) = 0. ! Assume no snow enter throug
+ ! boundaries
+ END WHERE
+ END IF
+ !
+ END DO
+ ENDIF
!
END SELECT
END IF
@@ -1052,6 +1154,22 @@ IF ( LCONDSAMP .AND. IMI == 1) THEN
ENDIF
ENDDO
ENDIF
+
+IF (LBLOWSNOW .AND. IMI == 1) THEN
+ IF (GFIRSTCALL3) THEN
+ ALLOCATE(GSNWBOUNDARY(NSV_SNW))
+ GFIRSTCALL3 = .FALSE.
+ DO JSV=NSV_SNWBEG,NSV_SNWEND
+ GCHTMP = .FALSE.
+ IF (LWEST_ll().AND.HLBCX(1)=='OPEN') GCHTMP = GCHTMP .OR. ALL(PLBXSVM(1,:,:,JSV)==0)
+ IF (LEAST_ll().AND.HLBCX(2)=='OPEN') GCHTMP = GCHTMP .OR. ALL(PLBXSVM(ILBX,:,:,JSV)==0)
+ IF (LSOUTH_ll().AND.HLBCY(1)=='OPEN') GCHTMP = GCHTMP .OR. ALL(PLBYSVM(:,1,:,JSV)==0)
+ IF (LNORTH_ll().AND.HLBCY(2)=='OPEN') GCHTMP = GCHTMP .OR. ALL(PLBYSVM(:,ILBY,:,JSV)==0)
+ GSNWBOUNDARY(JSV-NSV_SNWBEG+1) = GCHTMP
+ ENDDO
+ ENDIF
+ENDIF
+
#ifdef MNH_FOREFIRE
!ForeFire
IF ( LFOREFIRE .AND. IMI == 1) THEN
diff --git a/src/MNH/default_desfmn.f90 b/src/MNH/default_desfmn.f90
index bb139c807..6d4539257 100644
--- a/src/MNH/default_desfmn.f90
+++ b/src/MNH/default_desfmn.f90
@@ -1,7 +1,8 @@
-!MNH_LIC Copyright 1994-2018 CNRS, Meteo-France and Universite Paul Sabatier
+!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.
+! _Source: /home/cvsroot/MNH-VX-Y-Z/src/MNH/default_desfmn.f90,v _
!-----------------------------------------------------------------
! ###########################
MODULE MODI_DEFAULT_DESFM_n
@@ -218,9 +219,11 @@ END MODULE MODI_DEFAULT_DESFM_n
!! 10/2016 (C.Lac) Add droplet deposition
!! 10/2016 (R.Honnert and S.Riette) : Improvement of EDKF and adaptation to the grey zone
!! 10/2016 (F Brosse) add prod/loss terms computation for chemistry
+!! 07/2017 (V. Masson) adds time step for output files writing.
!! 09/2017 Q.Rodier add LTEND_UV_FRC
!! 02/2018 Q.Libois ECRAD
! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!! 07/2017 (V. Vionnet) add blowing snow variables
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -284,6 +287,8 @@ USE MODD_PARAM_LIMA, ONLY : LCOLD, LNUCL, LSEDI, LHHONI, LSNOW, LHAIL, LMEYERS,&
!
USE MODD_LATZ_EDFLX
USE MODD_2D_FRC
+USE MODD_BLOWSNOW
+USE MODD_BLOWSNOW_n
#ifdef MNH_FOREFIRE
USE MODD_FOREFIRE
#endif
@@ -425,6 +430,8 @@ LHORELAX_SVLIMA = .FALSE.
#ifdef MNH_FOREFIRE
LHORELAX_SVFF = .FALSE.
#endif
+LHORELAX_SVSNW = .FALSE.
+!
!
!-------------------------------------------------------------------------------
!
@@ -1379,4 +1386,17 @@ ENDIF
#endif
!-------------------------------------------------------------------------------
!
+!* 28. SET DEFAULT VALUES FOR MODD_BLOWSNOW AND MODD_BLOWSNOW_n
+! ----------------------------------------
+!
+IF (KMI == 1) THEN
+ LBLOWSNOW = .FALSE.
+ XALPHA_SNOW = 3.
+ XRSNOW = 4.
+ CSNOWSEDIM = 'TABC'
+END IF
+LSNOWSUBL = .FALSE.
+!
+!
+!
END SUBROUTINE DEFAULT_DESFM_n
diff --git a/src/MNH/diag.f90 b/src/MNH/diag.f90
index 632b7c6ea..ec3f8e4bf 100644
--- a/src/MNH/diag.f90
+++ b/src/MNH/diag.f90
@@ -86,6 +86,7 @@
!! 03/2018 (P.Wautelet) replace SUBTRACT_TO_DATE and ADD_FORECAST_TO_DATE
!! by DATETIME_CORRECTDATE
!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!! V.Vionnet 07/2017 add LWIND_CONTRAV
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -208,7 +209,7 @@ NAMELIST/NAM_DIAG/ CISO, LVAR_RS, LVAR_LS, &
LVAR_MRW, LVAR_MRSV, LVAR_FRC, &
LTPZH, LMOIST_V, LMOIST_E,LMOIST_ES, LCOREF, &
LVORT, LDIV, LMEAN_POVO, XMEAN_POVO, &
- LGEO, LAGEO, LWIND_ZM, LMSLP, LTHW, &
+ LGEO, LAGEO, LWIND_ZM,LWIND_CONTRAV, LMSLP, LTHW, &
LCLD_COV, LVAR_PR, LTOTAL_PR, LMEAN_PR, XMEAN_PR, &
NCAPE, LBV_FR, LRADAR, CBLTOP, LTRAJ, &
LDIAG,XDIAG,LCHEMDIAG,LCHAQDIAG,XCHEMLAT,XCHEMLON,&
@@ -275,6 +276,7 @@ XMEAN_POVO(2)=50000
LGEO=.FALSE.
LAGEO=.FALSE.
LWIND_ZM=.FALSE.
+LWIND_CONTRAV=.FALSE.
LMSLP=.FALSE.
LTHW=.FALSE.
LCLD_COV=.FALSE.
diff --git a/src/MNH/endstep.f90 b/src/MNH/endstep.f90
index ce319106c..c5086331e 100644
--- a/src/MNH/endstep.f90
+++ b/src/MNH/endstep.f90
@@ -190,7 +190,7 @@ END MODULE MODI_ENDSTEP
!! 10/2009 (C.Lac) Correction on FIT temporal scheme for variables
!! advected with PPM
!! 04/2013 (C.Lac) FIT for all the variables
-!! 04/2014 (C.Lac) Check on the positivity of XSVT
+!! 04/2014 (C.Lac) Check on the positivity of PSVT
!! J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1
!!
!------------------------------------------------------------------------------
@@ -205,12 +205,14 @@ USE MODD_GRID_n
USE MODD_BUDGET
USE MODD_NSV, ONLY : XSVMIN, NSV_CHEMBEG, NSV_CHEMEND, &
NSV_AERBEG, NSV_AEREND,&
- NSV_DSTBEG, NSV_DSTEND
-
+ NSV_DSTBEG, NSV_DSTEND,&
+ NSV_SNWBEG, NSV_SNWEND
USE MODD_CH_AEROSOL, ONLY : LORILAM
USE MODD_DUST, ONLY : LDUST
USE MODD_PARAM_C2R2, ONLY : LACTIT
USE MODD_LBC_n, ONLY : CLBCX, CLBCY
+USE MODD_BLOWSNOW
+USE MODD_BLOWSNOW_n
USE MODI_BUDGET
USE MODI_SHUMAN
!
@@ -316,6 +318,24 @@ IF (SIZE(PTKET,1) /= 0) PTKET(:,:,:)=PTKES(:,:,:)
!
PSVT(:,:,:,1:KSV)=PSVS(:,:,:,1:KSV)
!
+IF(LBLOWSNOW) THEN
+ DO JSV=1,(NBLOWSNOW_2D)
+ XSNWCANO(:,:,JSV) = XRSNWCANOS(:,:,JSV)
+ END DO
+!* MINIMUM VALUE FOR BLOWING SNOW
+!
+ WHERE(XSNWCANO(:,:,:)<1.E-20)
+ XSNWCANO(:,:,:)=0.
+ END WHERE
+
+ IF (SIZE(PSVT,4) > 1) THEN
+ WHERE(PSVT(:,:,:,NSV_SNWBEG:NSV_SNWEND)<1.E-20)
+ PSVT(:,:,:,NSV_SNWBEG:NSV_SNWEND)=0.
+ END WHERE
+ END IF
+!
+END IF
+!
IF (LWEST_ll( ) .AND. CLBCX(1)=='OPEN') THEN
DO JSV=1,KSV
PSVT(IIB,:,:,JSV)=MAX(PSVT(IIB,:,:,JSV),XSVMIN(JSV))
diff --git a/src/MNH/exchange.f90 b/src/MNH/exchange.f90
index 7d161d6aa..fbf47a060 100644
--- a/src/MNH/exchange.f90
+++ b/src/MNH/exchange.f90
@@ -3,6 +3,12 @@
!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!MNH_LIC for details. version 1.
!-----------------------------------------------------------------
+!--------------- special set of characters for RCS information
+!-----------------------------------------------------------------
+! $Source: /home/cvsroot/MNH-VX-Y-Z/src/MNH/exchange.f90,v $ $Revision: 1.2.2.2.2.2.16.1.2.5.2.1 $ $Date: 2015/12/01 15:26:23 $
+!-----------------------------------------------------------------
+!-----------------------------------------------------------------
+!-----------------------------------------------------------------
! ####################
MODULE MODI_EXCHANGE
! ####################
@@ -106,6 +112,8 @@ USE MODI_SUM_ll
USE MODI_BUDGET
USE MODD_CH_MNHC_n, ONLY : LUSECHEM, LUSECHAQ, LUSECHIC
USE MODD_CH_AEROSOL, ONLY : LORILAM, NM6_AER
+USE MODD_BLOWSNOW
+USE MODD_BLOWSNOW_n
!
IMPLICIT NONE
!
@@ -258,6 +266,13 @@ DO JSV=1,KSV
PRSVS(:,:,:,JSV) = PRSVS(:,:,:,JSV)*PTSTEP/PRHODJ
END DO
!
+IF(LBLOWSNOW) THEN
+ DO JSV=1,(NBLOWSNOW_2D)
+ XRSNWCANOS(:,:,JSV) = XRSNWCANOS(:,:,JSV)*PTSTEP/PRHODJ(:,:,1)
+ END DO
+END IF
+!
+!
!------------------------------------------------------------------------------
!
!* 2 UPDATE THE FIRST LAYER OF THE HALO
diff --git a/src/MNH/goto_model_wrapper.f90 b/src/MNH/goto_model_wrapper.f90
index 0a355f505..ffa3d19bf 100644
--- a/src/MNH/goto_model_wrapper.f90
+++ b/src/MNH/goto_model_wrapper.f90
@@ -14,6 +14,7 @@
!! 07/2017 (M.Leriche) Add DIAG chimical surface fluxes
! 02/2018 Q.Libois ECRAD
!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+! 2017 V.Vionnet blow snow
!-----------------------------------------------------------------
MODULE MODI_GOTO_MODEL_WRAPPER
@@ -87,6 +88,7 @@ USE MODD_SERIES_n
USE MODD_STATION_n
!USE MODD_TIME_n
USE MODD_TURB_n
+USE MODD_BLOWSNOW_n
!
USE MODD_SUB_CH_FIELD_VALUE_n
USE MODD_SUB_CH_MONITOR_n
@@ -221,6 +223,7 @@ CALL ADVFRC_GOTO_MODEL(KFROM, KTO)
CALL RELFRC_GOTO_MODEL(KFROM, KTO)
CALL CH_PRODLOSSTOT_GOTO_MODEL(KFROM,KTO)
CALL CH_BUDGET_GOTO_MODEL(KFROM,KTO)
+CALL BLOWSNOW_GOTO_MODEL(KFROM, KTO)
!
IF (.NOT.GNOFIELDLIST) CALL FIELDLIST_GOTO_MODEL(KFROM, KTO)
!
diff --git a/src/MNH/ground_paramn.f90 b/src/MNH/ground_paramn.f90
index 52bd7b52f..d6c91b20f 100644
--- a/src/MNH/ground_paramn.f90
+++ b/src/MNH/ground_paramn.f90
@@ -106,6 +106,8 @@ END MODULE MODI_GROUND_PARAM_n
!! 01/2018 (G.Delautier) SURFEX 8.1
!! 02/2018 Q.Libois ECRAD
!! (P.Wautelet) 28/03/2018 replace TEMPORAL_DIST by DATETIME_DISTANCE
+
+!! (V. Vionnet) 18/07/2017 add coupling for blowing snow module
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -120,7 +122,7 @@ USE MODD_METRICS_n, ONLY : XDXX, XDYY, XDZZ
USE MODD_DIM_n, ONLY : NKMAX
USE MODD_GRID_n, ONLY : XLON, XZZ, XDIRCOSXW, XDIRCOSYW, XDIRCOSZW, &
XCOSSLOPE, XSINSLOPE, XZS
-USE MODD_REF_n, ONLY : XRHODREF
+USE MODD_REF_n, ONLY : XRHODREF,XRHODJ
USE MODD_CONF_n, ONLY : NRR
USE MODD_PARAM_n, ONLY : CDCONV,CCLOUD, CRAD
USE MODD_PRECIP_n, ONLY : XINPRC, XINPRR, XINPRS, XINPRG, XINPRH
@@ -137,6 +139,8 @@ USE MODD_PARAM_C2R2, ONLY : LSEDC
USE MODD_DIAG_IN_RUN
USE MODD_DUST, ONLY : LDUST
USE MODD_SALT, ONLY : LSALT
+USE MODD_BLOWSNOW
+USE MODD_BLOWSNOW_n
USE MODD_CH_AEROSOL, ONLY : LORILAM
USE MODD_CSTS_DUST, ONLY : XMOLARWEIGHT_DUST
USE MODD_CSTS_SALT, ONLY : XMOLARWEIGHT_SALT
@@ -188,7 +192,7 @@ REAL, DIMENSION(:,:), INTENT(OUT) :: PSFV ! momentum in x and y directions
!
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PDIR_ALB ! direct albedo for each spectral band (-)
REAL, DIMENSION(:,:,:), INTENT(OUT) :: PSCA_ALB ! diffuse albedo for each spectral band (-)
-REAL, DIMENSION(:,:,:), INTENT(OUT) :: PEMIS ! surface emissivity (-)
+REAL, DIMENSION(:,:,:), INTENT(OUT) :: PEMIS ! surface emissivity (-)
REAL, DIMENSION(:,:), INTENT(OUT) :: PTSRAD ! surface radiative temperature (K)
!
!
@@ -240,7 +244,13 @@ 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
-
+!
+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
!
!* Dimensions
! ----------
@@ -262,6 +272,8 @@ INTEGER :: IDIM2 ! Y physical dimension
INTEGER :: IDIM1D! total physical dimension
INTEGER :: IKRAD
!
+INTEGER :: KSV_SURF ! Number of scalar variables sent to SURFEX
+!
!* Arrays put in 1D vectors
! ------------------------
!
@@ -319,6 +331,10 @@ REAL, DIMENSION(:), ALLOCATABLE :: ZP_MER10M ! meridian Wind at 10 meters
TYPE(LIST_ll), POINTER :: TZFIELDSURF_ll ! list of fields to exchange
INTEGER :: IINFO_ll ! return code of parallel routine
!
+!
+CHARACTER(LEN=6), DIMENSION(:), ALLOCATABLE :: YSV_SURF ! name of the scalar variables
+ ! sent to SURFEX
+!
REAL :: ZTIMEC
!
!-------------------------------------------------------------------------------
@@ -487,6 +503,36 @@ ZZREF(:,:) = 0.5*( XZZ(:,:,IKB+1)-XZZ(:,:,IKB) )*XDIRCOSZW(:,:)
!
ZCO2(:,:) = XCCO2 * XRHODREF(:,:,IKB)
!
+!
+!
+! 1.14 Blowing snow scheme (optional)
+! -----------------
+!
+ZBLOWSNOW_2D=0.
+
+IF(LBLOWSNOW) THEN
+ KSV_SURF = NSV+NBLOWSNOW_2D ! When blowing snow scheme is used
+ ! NBLOWSN0W_2D variables are sent to SURFEX through ZP_SV.
+ ! 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
+ ! Initialize array of scalar to be sent to SURFEX including 2D blowing snow fields
+ ALLOCATE(YSV_SURF(KSV_SURF))
+ YSV_SURF(1:NSV) = CSV(:)
+ YSV_SURF(NSV+1:KSV_SURF) = YPBLOWSNOW_2D(:)
+
+
+ DO JSV=1,NBLOWSNOW_2D
+ ZBLOWSNOW_2D(:,:,JSV) = XRSNWCANOS(:,:,JSV)*XTSTEP/XRHODJ(:,:,IKB)
+ END DO
+
+ELSE
+ KSV_SURF = NSV
+ ALLOCATE(YSV_SURF(KSV_SURF))
+ YSV_SURF(:) = CSV(:)
+ENDIF
+!
!-------------------------------------------------------------------------------
!
!* 2. Call to surface monitor with 2D variables
@@ -512,9 +558,9 @@ CALL DATETIME_DISTANCE(TDTSEG,TDTCUR,ZTIMEC)
!
CALL COUPLING_SURF_ATM_n(YSURF_CUR,'MESONH', 'E',ZTIMEC, &
XTSTEP, TDTCUR%TDATE%YEAR, TDTCUR%TDATE%MONTH, TDTCUR%TDATE%DAY, TDTCUR%TIME, &
- IDIM1D,NSV,SIZE(XSW_BANDS), &
+ 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, CSV,&
+ ZP_ZREF, ZP_ZREF, ZP_ZS, ZP_U, ZP_V, ZP_QA, ZP_TA, ZP_RHOA, ZP_SV, ZP_CO2, 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, &
@@ -629,6 +675,21 @@ ELSE
PSFSV(:,:,NSV_AERBEG:NSV_AEREND) = 0.
END IF
!
+!* conversion from blowing snow flux (kg/m2/s) to [kg(snow)/kg(dry air).m.s-1]
+!
+IF (LBLOWSNOW) THEN
+ DO JSV=NSV_SNWBEG,NSV_SNWEND
+ PSFSV(:,:,JSV) = ZSFTS(:,:,JSV)/ (ZRHOA(:,:))
+ END DO
+ !* Update tendency for blowing snow 2D fields
+ DO JSV=1,(NBLOWSNOW_2D)
+ XRSNWCANOS(:,:,JSV) = ZBLOWSNOW_2D(:,:,JSV)*XRHODJ(:,:,IKB)/(XTSTEP*ZRHOA(:,:))
+ END DO
+
+ELSE
+ PSFSV(:,:,NSV_SNWBEG:NSV_SNWEND) = 0.
+END IF
+!
!* conversion from CO2 flux (kg/m2/s) to w'CO2'
!
PSFCO2(:,:) = ZSFCO2(:,:) / XRHODREF(:,:,IKB)
@@ -699,7 +760,7 @@ ALLOCATE(ZP_V (KDIM1D))
ALLOCATE(ZP_QA (KDIM1D))
ALLOCATE(ZP_TA (KDIM1D))
ALLOCATE(ZP_RHOA (KDIM1D))
-ALLOCATE(ZP_SV (KDIM1D,NSV))
+ALLOCATE(ZP_SV (KDIM1D,KSV_SURF))
ALLOCATE(ZP_CO2 (KDIM1D))
ALLOCATE(ZP_RAIN (KDIM1D))
ALLOCATE(ZP_SNOW (KDIM1D))
@@ -713,7 +774,7 @@ ALLOCATE(ZP_SFTQ (KDIM1D))
ALLOCATE(ZP_SFTH (KDIM1D))
ALLOCATE(ZP_SFU (KDIM1D))
ALLOCATE(ZP_SFV (KDIM1D))
-ALLOCATE(ZP_SFTS (KDIM1D,NSV))
+ALLOCATE(ZP_SFTS (KDIM1D,KSV_SURF))
ALLOCATE(ZP_SFCO2 (KDIM1D))
ALLOCATE(ZP_TSRAD (KDIM1D))
ALLOCATE(ZP_DIR_ALB (KDIM1D,SIZE(PDIR_ALB,3)))
@@ -759,6 +820,12 @@ DO JLAYER=1,NSV
ZP_SV(:,JLAYER) = RESHAPE(XSVT(IIB:IIE,IJB:IJE,IKB,JLAYER), ISHAPE_1)
END DO
!
+IF(LBLOWSNOW) THEN
+ DO JLAYER=1,NBLOWSNOW_2D
+ ZP_SV(:,NSV+JLAYER) = RESHAPE(ZBLOWSNOW_2D(IIB:IIE,IJB:IJE,JLAYER), ISHAPE_1)
+ END DO
+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
@@ -775,6 +842,18 @@ DO JLAYER=NSV_SLTBEG,NSV_SLTEND
ZP_SV(:,JLAYER) = ZP_SV(:,JLAYER) * XMOLARWEIGHT_SALT* ZP_RHOA(:) / 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(:)
+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(:)
+ END DO
+END IF
+!
ZP_ZENITH(:) = RESHAPE(XZENITH(IIB:IIE,IJB:IJE), ISHAPE_1)
ZP_AZIM (:) = RESHAPE(XAZIM (IIB:IIE,IJB:IJE), ISHAPE_1)
ZP_LW(:) = RESHAPE(XFLALWD(IIB:IIE,IJB:IJE), ISHAPE_1)
@@ -820,6 +899,11 @@ DO JLAYER=1,SIZE(PEMIS,3)
PEMIS (IIB:IIE,IJB:IJE,JLAYER) = RESHAPE(ZP_EMIS(:), ISHAPE_2)
END DO
PTSRAD (IIB:IIE,IJB:IJE) = RESHAPE(ZP_TSRAD(:), ISHAPE_2)
+IF(LBLOWSNOW) THEN
+ DO JLAYER=1,NBLOWSNOW_2D
+ ZBLOWSNOW_2D(IIB:IIE,IJB:IJE,JLAYER) = RESHAPE(ZP_SFTS(:,NSV+JLAYER), ISHAPE_2)
+ END DO
+END IF
!
IF (LDIAG_IN_RUN) THEN
XCURRENT_RN (IIB:IIE,IJB:IJE) = RESHAPE(ZP_RN(:), ISHAPE_2)
diff --git a/src/MNH/ini_dynamics.f90 b/src/MNH/ini_dynamics.f90
index 6dba7d3d7..e4d00f5bb 100644
--- a/src/MNH/ini_dynamics.f90
+++ b/src/MNH/ini_dynamics.f90
@@ -14,7 +14,7 @@ SUBROUTINE INI_DYNAMICS(PLON,PLAT,PRHODJ,PTHVREF,PMAP,PZZ, &
OHORELAX_RH,OHORELAX_TKE,OHORELAX_SV, &
OHORELAX_SVC2R2,OHORELAX_SVC1R3,OHORELAX_SVELEC,OHORELAX_SVLG,&
OHORELAX_SVCHEM,OHORELAX_SVAER,OHORELAX_SVDST,OHORELAX_SVSLT, &
- OHORELAX_SVPP,OHORELAX_SVCS, OHORELAX_SVCHIC, &
+ OHORELAX_SVPP,OHORELAX_SVCS, OHORELAX_SVCHIC,OHORELAX_SVSNW, &
#ifdef MNH_FOREFIRE
OHORELAX_SVFF, &
#endif
@@ -88,6 +88,8 @@ LOGICAL, INTENT(IN):: OHORELAX_SVSLT ! switch for the
! horizontal relaxation for slt variables
LOGICAL, INTENT(IN):: OHORELAX_SVPP ! switch for the
! horizontal relaxation for passive pollutants
+LOGICAL, INTENT(IN):: OHORELAX_SVSNW ! switch for the
+ ! horizontal relaxation for blowing snow variables
#ifdef MNH_FOREFIRE
LOGICAL, INTENT(IN):: OHORELAX_SVFF ! switch for the
! horizontal relaxation for ForeFire variables
@@ -183,7 +185,7 @@ SUBROUTINE INI_DYNAMICS(PLON,PLAT,PRHODJ,PTHVREF,PMAP,PZZ, &
OHORELAX_RH,OHORELAX_TKE,OHORELAX_SV, &
OHORELAX_SVC2R2,OHORELAX_SVC1R3,OHORELAX_SVELEC,OHORELAX_SVLG,&
OHORELAX_SVCHEM,OHORELAX_SVAER,OHORELAX_SVDST,OHORELAX_SVSLT, &
- OHORELAX_SVPP,OHORELAX_SVCS, OHORELAX_SVCHIC, &
+ OHORELAX_SVPP,OHORELAX_SVCS, OHORELAX_SVCHIC,OHORELAX_SVSNW, &
#ifdef MNH_FOREFIRE
OHORELAX_SVFF, &
#endif
@@ -277,6 +279,7 @@ SUBROUTINE INI_DYNAMICS(PLON,PLAT,PRHODJ,PTHVREF,PMAP,PZZ, &
!! Modification 20/05/06 Remove KEPS
!! Modification 07/2013 (Bosseur & Filippi) Adds Forefire
!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!! Vionnet 07/2017 : blow snow
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -355,6 +358,8 @@ LOGICAL, INTENT(IN):: OHORELAX_SVSLT ! switch for the
! horizontal relaxation for slt variables
LOGICAL, INTENT(IN):: OHORELAX_SVPP ! switch for the
! horizontal relaxation for passive pollutants
+LOGICAL, INTENT(IN):: OHORELAX_SVSNW ! switch for the
+ ! horizontal relaxation for blowing snow variables
#ifdef MNH_FOREFIRE
LOGICAL, INTENT(IN):: OHORELAX_SVFF ! switch for the
! horizontal relaxation for ForeFire variables
@@ -447,9 +452,9 @@ INTEGER :: IIU,IJU ! Upper bounds in x,y dir
LOGICAL :: GHORELAX
LOGICAL, DIMENSION(7) :: GHORELAXR ! local array of logical
#ifdef MNH_FOREFIRE
-LOGICAL, DIMENSION(12):: GHORELAXSV! local array of logical
+LOGICAL, DIMENSION(13):: GHORELAXSV! local array of logical
#else
-LOGICAL, DIMENSION(11):: GHORELAXSV! local array of logical
+LOGICAL, DIMENSION(12):: GHORELAXSV! local array of logical
#endif
!
!-------------------------------------------------------------------------------
@@ -519,8 +524,9 @@ GHORELAXSV(8) = OHORELAX_SVSLT
GHORELAXSV(9) = OHORELAX_SVPP
GHORELAXSV(10)= OHORELAX_SVCS
GHORELAXSV(11) = OHORELAX_SVCHIC
+GHORELAXSV(12) = OHORELAX_SVSNW
#ifdef MNH_FOREFIRE
-GHORELAXSV(12) = OHORELAX_SVFF
+GHORELAXSV(13) = OHORELAX_SVFF
#endif
!
GHORELAX=ANY(GHORELAXR) .OR. ANY(GHORELAXSV) .OR. ANY(OHORELAX_SV) &
@@ -532,7 +538,7 @@ IF (GHORELAX .OR. OVE_RELAX.OR.OVE_RELAX_GRD) THEN
OHORELAX_RH,OHORELAX_TKE,OHORELAX_SV, &
OHORELAX_SVC2R2,OHORELAX_SVC1R3,OHORELAX_SVELEC,OHORELAX_SVLG, &
OHORELAX_SVCHEM, OHORELAX_SVAER, OHORELAX_SVDST, OHORELAX_SVSLT, &
- OHORELAX_SVPP, OHORELAX_SVCS, OHORELAX_SVCHIC, &
+ OHORELAX_SVPP, OHORELAX_SVCS, OHORELAX_SVCHIC,OHORELAX_SVSNW, &
PALKTOP,PALKGRD, PALZBOT,PALZBAS, &
PZZ, PZHAT, PTSTEP, &
PRIMKMAX,KRIMX,KRIMY, &
diff --git a/src/MNH/ini_micron.f90 b/src/MNH/ini_micron.f90
index ff9761e5f..67bce21b6 100644
--- a/src/MNH/ini_micron.f90
+++ b/src/MNH/ini_micron.f90
@@ -77,6 +77,8 @@ USE MODD_CLOUDPAR_n, ONLY : NSPLITR, NSPLITG
USE MODD_PARAM_n, ONLY : CELEC
USE MODD_PARAM_ICE, ONLY : LSEDIC, LDEPOSC
USE MODD_PARAM_C2R2, ONLY : LSEDC, LACTIT, LDEPOC
+USE MODD_BLOWSNOW
+USE MODD_BLOWSNOW_n
!
USE MODI_READ_PRECIP_FIELD
USE MODI_INI_CLOUD
@@ -90,6 +92,7 @@ USE MODI_SET_CONC_ICE_C1R3
!
USE MODE_ll
USE MODE_MODELN_HANDLER
+USE MODE_BLOWSNOW_SEDIM_LKT
!
USE MODD_NSV, ONLY : NSV,NSV_CHEM,NSV_C2R2BEG,NSV_C2R2END, &
NSV_C1R3BEG,NSV_C1R3END, &
@@ -211,6 +214,19 @@ ELSE
ALLOCATE(XACPRH(0,0))
END IF
!
+IF(LBLOWSNOW) THEN
+ ALLOCATE(XSNWSUBL3D(IIU,IJU,IKU))
+ XSNWSUBL3D(:,:,:) = 0.0
+ IF(CSNOWSEDIM=='TABC') THEN
+!Read in look up tables of snow particles properties
+!No arguments, all look up tables are defined in module
+!mode_snowdrift_sedim_lkt
+ CALL BLOWSNOW_SEDIM_LKT_SET
+ END IF
+ELSE
+ ALLOCATE(XSNWSUBL3D(0,0,0))
+END IF
+!
IF(SIZE(XINPRR) == 0) RETURN
!
!* 2b. ALLOCATION for Radiative cooling
diff --git a/src/MNH/ini_modeln.f90 b/src/MNH/ini_modeln.f90
index a8dd7420b..b2c466e82 100644
--- a/src/MNH/ini_modeln.f90
+++ b/src/MNH/ini_modeln.f90
@@ -275,6 +275,7 @@ END MODULE MODI_INI_MODEL_n
!! 09/2017 Q.Rodier add LTEND_UV_FRC
!! 02/2018 Q.Libois ECRAD
!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!! V. Vionnet : 18/07/2017 : add blowing snow scheme
!---------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -353,6 +354,8 @@ USE MODD_TURB_n
USE MODD_CTURB
USE MODD_LBC_n
USE MODD_PASPOL_n
+USE MODD_BLOWSNOW
+USE MODD_BLOWSNOW_n
!
!
USE MODI_INI_BUDGET
@@ -801,6 +804,16 @@ IF (LPASPOL) THEN
XATC = 0.
END IF
!
+IF(LBLOWSNOW) THEN
+ ALLOCATE(XSNWCANO(IIU,IJU,NBLOWSNOW_2D))
+ ALLOCATE(XRSNWCANOS(IIU,IJU,NBLOWSNOW_2D))
+ XSNWCANO(:,:,:) = 0.0
+ XRSNWCANOS(:,:,:) = 0.0
+ELSE
+ ALLOCATE(XSNWCANO(0,0,0))
+ ALLOCATE(XRSNWCANOS(0,0,0))
+END IF
+!
!* 3.2 Module MODD_GRID_n and MODD_METRICS_n
!
IF (LCARTESIAN) THEN
@@ -1923,7 +1936,7 @@ CALL INI_DYNAMICS(XLON,XLAT,XRHODJ,XTHVREF,XMAP,XZZ,XDXHAT,XDYHAT, &
LHORELAX_RH,LHORELAX_TKE,LHORELAX_SV, &
LHORELAX_SVC2R2,LHORELAX_SVC1R3,LHORELAX_SVELEC,LHORELAX_SVLG, &
LHORELAX_SVCHEM,LHORELAX_SVAER,LHORELAX_SVDST,LHORELAX_SVSLT, &
- LHORELAX_SVPP,LHORELAX_SVCS,LHORELAX_SVCHIC, &
+ LHORELAX_SVPP,LHORELAX_SVCS,LHORELAX_SVCHIC,LHORELAX_SVSNW, &
#ifdef MNH_FOREFIRE
LHORELAX_SVFF, &
#endif
diff --git a/src/MNH/ini_nsv.f90 b/src/MNH/ini_nsv.f90
index 00407e247..5299ffb8a 100644
--- a/src/MNH/ini_nsv.f90
+++ b/src/MNH/ini_nsv.f90
@@ -63,6 +63,7 @@ END MODULE MODI_INI_NSV
!! the 4th C2R2 scalar variable
!! J.escobar 04/08/2015 suit Pb with writ_lfin JSA increment , modif in ini_nsv to have good order initialization
!! Modification 01/2016 (JP Pinty) Add LIMA and LUSECHEM condition
+!! Modification 07/2017 (V. Vionnet) Add blowing snow condition
!!
!-------------------------------------------------------------------------------
!
@@ -79,7 +80,8 @@ USE MODD_DYN_n, ONLY : LHORELAX_SV,LHORELAX_SVC2R2,LHORELAX_SVC1R3, &
LHORELAX_SVLIMA, &
LHORELAX_SVELEC,LHORELAX_SVCHEM,LHORELAX_SVLG, &
LHORELAX_SVDST,LHORELAX_SVAER, LHORELAX_SVSLT, &
- LHORELAX_SVPP,LHORELAX_SVCS, LHORELAX_SVCHIC
+ LHORELAX_SVPP,LHORELAX_SVCS, LHORELAX_SVCHIC, &
+ LHORELAX_SVSNW
#ifdef MNH_FOREFIRE
USE MODD_DYN_n, ONLY : LHORELAX_SVFF
USE MODD_FOREFIRE
@@ -89,6 +91,7 @@ USE MODD_LG
USE MODD_DUST
USE MODD_SALT
USE MODD_PASPOL
+USE MODD_BLOWSNOW
USE MODD_CONDSAMP
USE MODD_CH_AEROSOL
USE MODD_PREP_REAL, ONLY: XT_LS
@@ -482,6 +485,21 @@ ELSE
! in order to create a null section
END IF
!
+! scalar variables used in blowing snow model
+!
+IF (LBLOWSNOW) THEN
+ NSV_SNW_A(KMI) = NBLOWSNOW3D
+ NSV_SNWBEG_A(KMI)= ISV+1
+ NSV_SNWEND_A(KMI)= ISV+NSV_SNW_A(KMI)
+ ISV = NSV_SNWEND_A(KMI)
+ELSE
+ NSV_SNW_A(KMI) = 0
+! force First index to be superior to last index
+! in order to create a null section
+ NSV_SNWBEG_A(KMI)= 1
+ NSV_SNWEND_A(KMI)= 0
+END IF
+!
! scalar variables used as LiNOX passive tracer
!
! In case without chemistry
@@ -552,6 +570,9 @@ LHORELAX_SV(NSV_FFBEG_A(KMI):NSV_FFEND_A(KMI))=LHORELAX_SVFF
! Conditional sampling
IF (LCONDSAMP) &
LHORELAX_SV(NSV_CSBEG_A(KMI):NSV_CSEND_A(KMI))=LHORELAX_SVCS
+! Blowing snow case
+IF (LBLOWSNOW) &
+LHORELAX_SV(NSV_SNWBEG_A(KMI):NSV_SNWEND_A(KMI))=LHORELAX_SVSNW
! Update NSV* variables for model KMI
CALL UPDATE_NSV(KMI)
!
@@ -591,7 +612,8 @@ IF (LPASPOL) XSVMIN(NSV_PPBEG_A(KMI):NSV_PPEND_A(KMI))=0.
#ifdef MNH_FOREFIRE
IF (LFOREFIRE) XSVMIN(NSV_FFBEG_A(KMI):NSV_FFEND_A(KMI))=0.
#endif
-IF (LCONDSAMP) XSVMIN(NSV_CSBEG_A(KMI):NSV_CSEND_A(KMI))=0.
+IF (LCONDSAMP) XSVMIN(NSV_CSBEG_A(KMI):NSV_CSEND_A(KMI))=0.
+IF (LBLOWSNOW) XSVMIN(NSV_SNWBEG_A(KMI):NSV_SNWEND_A(KMI))=XMNH_TINY
!
! NAME OF THE SCALAR VARIABLES IN THE DIFFERENT SV GROUPS
!
diff --git a/src/MNH/ini_segn.f90 b/src/MNH/ini_segn.f90
index 4e86a1c01..48e2a16bd 100644
--- a/src/MNH/ini_segn.f90
+++ b/src/MNH/ini_segn.f90
@@ -162,6 +162,7 @@ END MODULE MODI_INI_SEG_n
!! 01/2015 add GLNOX_EXPLICIT (C. Barthe)
!! 04/2016 add ABORT if CINIFILEPGD is not specified (G.Delautier)
!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!! 07/2017 add GBLOWSNOW (V. Vionnet)
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -232,6 +233,7 @@ LOGICAL :: GPASPOL
LOGICAL :: GFOREFIRE
#endif
LOGICAL :: GCONDSAMP
+LOGICAL :: GBLOWSNOW
LOGICAL :: GCHTRANS
LOGICAL :: GLNOX_EXPLICIT ! flag for LNOx
! These variables
@@ -358,7 +360,7 @@ CALL READ_DESFM_n(KMI,TPINIFILE,YCONF,GFLAT,GUSERV,GUSERC, &
GFOREFIRE, &
#endif
GLNOX_EXPLICIT, &
- GCONDSAMP, IRIMX,IRIMY,ISV, &
+ GCONDSAMP,GBLOWSNOW, IRIMX,IRIMY,ISV, &
YTURB,YTOM,GRMC01,YRAD,YDCONV,YSCONV,YCLOUD,YELEC,YEQNSYS )
!
!-------------------------------------------------------------------------------
@@ -448,7 +450,7 @@ CALL READ_EXSEG_n(KMI,TZFILE_DES,YCONF,GFLAT,GUSERV,GUSERC, &
GFOREFIRE, &
#endif
GLNOX_EXPLICIT, &
- GCONDSAMP, IRIMX,IRIMY,ISV, &
+ GCONDSAMP,GBLOWSNOW, IRIMX,IRIMY,ISV, &
YTURB,YTOM,GRMC01,YRAD,YDCONV,YSCONV,YCLOUD,YELEC,YEQNSYS, &
PTSTEP_ALL,CSTORAGE_TYPE,CINIFILEPGD_n )
!
diff --git a/src/MNH/ini_spectren.f90 b/src/MNH/ini_spectren.f90
index 31363e2e7..d2564f410 100644
--- a/src/MNH/ini_spectren.f90
+++ b/src/MNH/ini_spectren.f90
@@ -90,6 +90,8 @@ USE MODD_PARAM_n
USE MODD_PARAM_RAD_n, ONLY: CLW, CAER
USE MODD_PASPOL
USE MODD_PASPOL_n
+USE MODD_BLOWSNOW
+USE MODD_BLOWSNOW_n
USE MODD_PAST_FIELD_n
USE MODD_RADIATIONS_n
USE MODD_REF
@@ -923,7 +925,7 @@ CALL INI_DYNAMICS(XLON,XLAT,XRHODJ,XTHVREF,XMAP,XZZ,XDXHAT,XDYHAT, &
LHORELAX_RH,LHORELAX_TKE,LHORELAX_SV, &
LHORELAX_SVC2R2,LHORELAX_SVC1R3,LHORELAX_SVELEC,LHORELAX_SVLG, &
LHORELAX_SVCHEM,LHORELAX_SVAER,LHORELAX_SVDST,LHORELAX_SVSLT, &
- LHORELAX_SVPP,LHORELAX_SVCS,LHORELAX_SVCHIC, &
+ LHORELAX_SVPP,LHORELAX_SVCS,LHORELAX_SVCHIC,LHORELAX_SVSNW, &
#ifdef MNH_FOREFIRE
LHORELAX_SVFF, &
#endif
diff --git a/src/MNH/init_ground_paramn.f90 b/src/MNH/init_ground_paramn.f90
index d705c8ef6..95a5bcee7 100644
--- a/src/MNH/init_ground_paramn.f90
+++ b/src/MNH/init_ground_paramn.f90
@@ -67,6 +67,7 @@ END MODULE MODI_INIT_GROUND_PARAM_n
!! 01/2018 (G.Delautier) SURFEX 8.1
!! 02/2018 Q.Libois ECRAD
!! 03/2018 (P.Wautelet) replace ADD_FORECAST_TO_DATE_SURF by DATETIME_CORRECTDATE
+!! 2017 V.Vionnet Blow snow
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -84,6 +85,7 @@ USE MODD_NSV
USE MODD_DUST, ONLY : CDUSTNAMES
USE MODD_SALT, ONLY : CSALTNAMES
USE MODD_CH_AEROSOL, ONLY : CAERONAMES
+USE MODD_BLOWSNOW
USE MODD_TYPE_DATE, ONLY : DATE_TIME
!
USE MODD_TYPE_DATE_SURF, ONLY : DATE_SURF=>DATE
@@ -144,6 +146,8 @@ INTEGER :: IID,IRESP
TYPE (DATE_TIME), POINTER :: TZTCUR=>NULL()
TYPE (DATE_TIME) :: TZDATE
!
+CHARACTER(LEN=6), DIMENSION(:), ALLOCATABLE :: YSV_SURF ! name of the scalar variables
+ ! sent to SURFEX
!-------------------------------------------------------------------------------
!
!
@@ -194,10 +198,29 @@ TDATE_END%YEAR = TZDATE%TDATE%YEAR
TDATE_END%MONTH = TZDATE%TDATE%MONTH
TDATE_END%DAY = TZDATE%TDATE%DAY
!
+DO JLAYER=NSV_SNWBEG,NSV_SNWEND
+ HSV(JLAYER) = TRIM(CSNOWNAMES(JLAYER-NSV_SNWBEG+1))
+END DO
+!
ISV = SIZE(HSV)
+IF(LBLOWSNOW) THEN
+ ISV = ISV+NBLOWSNOW_2D ! When blowing snow scheme is used
+ ! NBLOWSN0W_2D variables are sent to SURFEX through ZP_SV.
+ ! 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
+
+ ALLOCATE(YSV_SURF(ISV))
+ YSV_SURF(1:KSV) = HSV(:)
+ YSV_SURF(NSV+1:ISV) = YPBLOWSNOW_2D(:)
+ELSE
+ ALLOCATE(YSV_SURF(ISV))
+ YSV_SURF(:) = HSV(:)
+ENDIF
CALL INIT_SURF_ATM_n(YSURF_CUR,'MESONH',HINIT,.FALSE., &
ILU,ISV,SIZE(PSW_BANDS), &
- HSV,ZCO2,ZRHODREF, &
+ YSV_SURF,ZCO2,ZRHODREF, &
ZZENITH,ZAZIM,PSW_BANDS,ZDIR_ALB,ZSCA_ALB, &
ZEMIS,ZTSRAD,ZTSURF, &
TZTCUR%TDATE%YEAR, TZTCUR%TDATE%MONTH, &
@@ -225,6 +248,7 @@ DEALLOCATE(ZDIR_ALB)
DEALLOCATE(ZSCA_ALB)
DEALLOCATE(ZEMIS )
DEALLOCATE(ZTSRAD )
+DEALLOCATE(YSV_SURF )
!
!-------------------------------------------------------------------------------
!
diff --git a/src/MNH/initial_guess.f90 b/src/MNH/initial_guess.f90
index 32a32f530..b999796d6 100644
--- a/src/MNH/initial_guess.f90
+++ b/src/MNH/initial_guess.f90
@@ -153,6 +153,8 @@ END MODULE MODI_INITIAL_GUESS
USE MODD_CONF
USE MODD_GRID_n
USE MODD_BUDGET
+USE MODD_BLOWSNOW
+USE MODD_BLOWSNOW_n
!
USE MODI_SHUMAN
USE MODI_BUDGET
@@ -224,6 +226,12 @@ DO JSV=1,KSV
PRSVS(:,:,:,JSV) = PSVT(:,:,:,JSV) * ZINVTSTEP * PRHODJ(:,:,:)
END DO
!
+IF(LBLOWSNOW) THEN
+ DO JSV=1,(NBLOWSNOW_2D)
+ XRSNWCANOS(:,:,JSV) = XSNWCANO(:,:,JSV) * ZINVTSTEP * PRHODJ(:,:,1)
+ END DO
+END IF
+!
IF (LBU_ENABLE) THEN
IF (LBU_BEG) THEN
NBUPROCCTR(:)=1
diff --git a/src/MNH/modd_diag_flag.f90 b/src/MNH/modd_diag_flag.f90
index f5ac8fbcf..4947fdfa9 100644
--- a/src/MNH/modd_diag_flag.f90
+++ b/src/MNH/modd_diag_flag.f90
@@ -85,6 +85,7 @@ REAL, DIMENSION(2) :: XMEAN_POVO
LOGICAL :: LGEO ! Geostrophic wind (m/s)
LOGICAL :: LAGEO ! Ageostrophic wind (m/s)
LOGICAL :: LWIND_ZM ! Zonal and meridien components of wind (m/s)
+LOGICAL :: LWIND_CONTRAV ! Contravariant component of wind (m/s)
LOGICAL :: LMSLP ! Mean Sea Level Pressure (hPa)
LOGICAL :: LTHW
LOGICAL :: LCLD_COV
diff --git a/src/MNH/modd_dynn.f90 b/src/MNH/modd_dynn.f90
index 821b2129a..c3b64335c 100644
--- a/src/MNH/modd_dynn.f90
+++ b/src/MNH/modd_dynn.f90
@@ -42,6 +42,7 @@
!! Modifications 07/10 (M.Leriche) Add relaxation for ice phase chemical
!! Modification 01/2016 (JP Pinty) Add LIMA
!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!! Modification 07/2017 (V. Vionnet) Add blowing snow variable
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -135,6 +136,7 @@ TYPE DYN_t
LOGICAL :: LHORELAX_SVFF
#endif
LOGICAL :: LHORELAX_SVCS
+ LOGICAL :: LHORELAX_SVSNW
LOGICAL, DIMENSION(:),POINTER :: LHORELAX_SV =>NULL()
!
REAL :: XRIMKMAX ! Max. value of the horiz. relaxation coeff.
@@ -224,6 +226,7 @@ LOGICAL, POINTER :: LHORELAX_SVPP=>NULL()
LOGICAL, POINTER :: LHORELAX_SVFF=>NULL()
#endif
LOGICAL, POINTER :: LHORELAX_SVCS=>NULL()
+LOGICAL, POINTER :: LHORELAX_SVSNW=>NULL()
LOGICAL, DIMENSION(:), POINTER :: LHORELAX_SV=>NULL()
REAL, POINTER :: XRIMKMAX=>NULL()
INTEGER, POINTER :: NRIMX=>NULL(),NRIMY=>NULL()
@@ -331,6 +334,7 @@ LHORELAX_SVPP=>DYN_MODEL(KTO)%LHORELAX_SVPP
LHORELAX_SVFF=>DYN_MODEL(KTO)%LHORELAX_SVFF
#endif
LHORELAX_SVCS=>DYN_MODEL(KTO)%LHORELAX_SVCS
+LHORELAX_SVSNW=>DYN_MODEL(KTO)%LHORELAX_SVSNW
LHORELAX_SV=>DYN_MODEL(KTO)%LHORELAX_SV
XRIMKMAX=>DYN_MODEL(KTO)%XRIMKMAX
!NRIMX=>DYN_MODEL(KTO)%NRIMX !Done in FIELDLIST_GOTO_MODEL
diff --git a/src/MNH/modd_nsv.f90 b/src/MNH/modd_nsv.f90
index 893242925..2417f728a 100644
--- a/src/MNH/modd_nsv.f90
+++ b/src/MNH/modd_nsv.f90
@@ -26,6 +26,7 @@
!! C.Lac 07/11 add conditional sampling
!! Pialat/Tulet 15/02/12 add ForeFire
!! Modification 01/2016 (JP Pinty) Add LIMA
+!! V. Vionnet 07/17 add blowing snow
!!
!-------------------------------------------------------------------------------
!
@@ -135,6 +136,10 @@ INTEGER,DIMENSION(JPMODELMAX)::NSV_FFBEG_A = 0 ! with indices in the range :
INTEGER,DIMENSION(JPMODELMAX)::NSV_FFEND_A = 0 ! NSV_FFBEG_A...NSV_FFEND_A
#endif
!
+INTEGER,DIMENSION(JPMODELMAX)::NSV_SNW_A = 0 ! number of blowing snow scalar
+INTEGER,DIMENSION(JPMODELMAX)::NSV_SNWBEG_A = 0 ! with indices in the range :
+INTEGER,DIMENSION(JPMODELMAX)::NSV_SNWEND_A = 0 ! NSV_SNWBEG_A...NSV_SNWEND_A
+!
!###############################################################################
!
! variables updated for the current model
@@ -231,5 +236,9 @@ INTEGER :: NSV_FF = 0 ! number of ForeFire scalar variables
INTEGER :: NSV_FFBEG = 0 ! with indices in the range :
INTEGER :: NSV_FFEND = 0 ! NSV_FFBEG...NSV_FFEND
#endif
+!
+INTEGER :: NSV_SNW = 0 ! number of blowing snow scalar variables
+INTEGER :: NSV_SNWBEG = 0 ! with indices in the range :
+INTEGER :: NSV_SNWEND = 0 ! NSV_SNWBEG...NSV_SNWEND
END MODULE MODD_NSV
diff --git a/src/MNH/modeln.f90 b/src/MNH/modeln.f90
index b844e57a8..cb56e15d0 100644
--- a/src/MNH/modeln.f90
+++ b/src/MNH/modeln.f90
@@ -248,6 +248,7 @@ END MODULE MODI_MODEL_n
!! 01/2018 (G.Delautier) SURFEX 8.1
!! 03/2018 (P.Wautelet) replace ADD_FORECAST_TO_DATE by DATETIME_CORRECTDATE
!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!! 07/2017 (V. Vionnet) : Add blowing snow scheme
!!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -308,6 +309,8 @@ USE MODD_PARAM_LIMA, ONLY: MSEDC => LSEDC, MWARM => LWARM, MRAIN => LRAIN, L
MACTIT => LACTIT, LSCAV, NMOD_CCN, LCOLD, &
MSEDI => LSEDI, MHHONI => LHHONI, NMOD_IFN, LHAIL, &
XRTMIN_LIMA=>XRTMIN
+USE MODD_BLOWSNOW_n
+USE MODD_BLOWSNOW
USE MODD_PARAM_MFSHALL_n
USE MODD_PARAM_n
USE MODD_PAST_FIELD_n
@@ -404,6 +407,8 @@ USE MODI_WRITE_LFIFMN_FORDIACHRO_n
USE MODI_WRITE_PROFILER_n
USE MODI_WRITE_SERIES_n
USE MODI_WRITE_STATION_n
+USE MODI_BLOWSNOW
+
USE MODI_WRITE_SURF_ATM_N
!
IMPLICIT NONE
@@ -1247,7 +1252,9 @@ END DO
DO JSV = NSV_AERDEPBEG,NSV_AERDEPEND
XRSVS(:,:,:,JSV) = MAX(XRSVS(:,:,:,JSV),0.)
END DO
-
+DO JSV = NSV_SNWBEG,NSV_SNWEND
+ XRSVS(:,:,:,JSV) = MAX(XRSVS(:,:,:,JSV),0.)
+END DO
IF (CELEC .NE. 'NONE') THEN
XRSVS(:,:,:,NSV_ELECBEG) = MAX(XRSVS(:,:,:,NSV_ELECBEG),0.)
XRSVS(:,:,:,NSV_ELECEND) = MAX(XRSVS(:,:,:,NSV_ELECEND),0.)
@@ -1278,7 +1285,7 @@ IF(LVE_RELAX .OR. LVE_RELAX_GRD .OR. LHORELAX_UVWTH .OR. LHORELAX_RV .OR.&
LHORELAX_SVELEC,LHORELAX_SVLG, &
LHORELAX_SVCHEM,LHORELAX_SVCHIC,LHORELAX_SVAER, &
LHORELAX_SVDST,LHORELAX_SVSLT,LHORELAX_SVPP, &
- LHORELAX_SVCS, &
+ LHORELAX_SVCS,LHORELAX_SVSNW, &
#ifdef MNH_FOREFIRE
LHORELAX_SVFF, &
#endif
@@ -1437,6 +1444,12 @@ IF (.NOT. LSTEADYLS) THEN
XLBXSVS(:,:,:,JSV)=MAX(XLBXSVS(:,:,:,JSV),0.)
XLBYSVS(:,:,:,JSV)=MAX(XLBYSVS(:,:,:,JSV),0.)
ENDDO
+ !
+ DO JSV=NSV_SNWBEG,NSV_SNWEND
+ XLBXSVS(:,:,:,JSV)=MAX(XLBXSVS(:,:,:,JSV),0.)
+ XLBYSVS(:,:,:,JSV)=MAX(XLBYSVS(:,:,:,JSV),0.)
+ ENDDO
+ !
END IF
END IF
END IF
@@ -1448,6 +1461,18 @@ XT_COUPL = XT_COUPL + ZTIME2 - ZTIME1
!-------------------------------------------------------------------------------
!
!
+!
+!* 8 Bis . Blowing snow scheme
+! ---------
+!
+IF(LBLOWSNOW) THEN
+ CALL BLOWSNOW(CLBCX,CLBCY,XTSTEP,NRR,XPABST,XTHT,XRT,XZZ,XRHODREF, &
+ XRHODJ,XEXNREF,XRRS,XRTHS,XSVT,XRSVS,XSNWSUBL3D )
+ENDIF
+!
+!
+!-------------------------------------------------------------------------------
+!
!* 9. ADVECTION
! ---------
!
diff --git a/src/MNH/modn_dynn.f90 b/src/MNH/modn_dynn.f90
index f43e058b1..0f6585f99 100644
--- a/src/MNH/modn_dynn.f90
+++ b/src/MNH/modn_dynn.f90
@@ -92,6 +92,7 @@ USE MODD_DYN_n, ONLY : &
LHORELAX_SVDST_n => LHORELAX_SVDST, &
LHORELAX_SVSLT_n => LHORELAX_SVSLT, &
LHORELAX_SVAER_n => LHORELAX_SVAER, &
+ LHORELAX_SVSNW_n => LHORELAX_SVSNW, &
LHORELAX_SV_n => LHORELAX_SV, &
LVE_RELAX_n => LVE_RELAX, &
LVE_RELAX_GRD_n => LVE_RELAX_GRD, &
@@ -134,6 +135,7 @@ LOGICAL, SAVE :: LHORELAX_SVCS
LOGICAL, SAVE :: LHORELAX_SVDST
LOGICAL, SAVE :: LHORELAX_SVSLT
LOGICAL, SAVE :: LHORELAX_SVAER
+LOGICAL, SAVE :: LHORELAX_SVSNW
LOGICAL, SAVE, DIMENSION(JPSVMAX) :: LHORELAX_SV
LOGICAL,SAVE :: LVE_RELAX
LOGICAL,SAVE :: LVE_RELAX_GRD
@@ -149,7 +151,7 @@ NAMELIST/NAM_DYNn/XTSTEP,CPRESOPT,NITR,LITRADJ,LRES,XRES,XRELAX,LHORELAX_UVWTH,
LHORELAX_RS,LHORELAX_RG,LHORELAX_RH,LHORELAX_TKE,&
LHORELAX_SVC2R2,LHORELAX_SVC1R3, LHORELAX_SVELEC,&
LHORELAX_SVCHEM, LHORELAX_SVCHIC, LHORELAX_SVLG, LHORELAX_SVDST, &
- LHORELAX_SVSLT, LHORELAX_SVAER, LHORELAX_SVPP, &
+ LHORELAX_SVSLT, LHORELAX_SVAER, LHORELAX_SVPP,LHORELAX_SVSNW, &
LHORELAX_SVCS, LHORELAX_SV,LVE_RELAX,LVE_RELAX_GRD,&
#ifdef MNH_FOREFIRE
LHORELAX_SVFF, &
@@ -190,6 +192,7 @@ SUBROUTINE INIT_NAM_DYNn
LHORELAX_SVDST = LHORELAX_SVDST_n
LHORELAX_SVSLT = LHORELAX_SVSLT_n
LHORELAX_SVAER = LHORELAX_SVAER_n
+ LHORELAX_SVSNW = LHORELAX_SVSNW_n
LHORELAX_SV = LHORELAX_SV_n
LVE_RELAX = LVE_RELAX_n
LVE_RELAX_GRD = LVE_RELAX_GRD_n
@@ -232,6 +235,7 @@ SUBROUTINE UPDATE_NAM_DYNn
LHORELAX_SVDST_n = LHORELAX_SVDST
LHORELAX_SVSLT_n = LHORELAX_SVSLT
LHORELAX_SVAER_n = LHORELAX_SVAER
+ LHORELAX_SVSNW_n = LHORELAX_SVSNW
LHORELAX_SV_n = LHORELAX_SV
LVE_RELAX_n = LVE_RELAX
LVE_RELAX_GRD_n = LVE_RELAX_GRD
diff --git a/src/MNH/read_desfmn.f90 b/src/MNH/read_desfmn.f90
index 2f9fea8e4..98abbf158 100644
--- a/src/MNH/read_desfmn.f90
+++ b/src/MNH/read_desfmn.f90
@@ -18,7 +18,7 @@ INTERFACE
OFOREFIRE, &
#endif
OLNOX_EXPLICIT, &
- OCONDSAMP, &
+ OCONDSAMP,OBLOWSNOW, &
KRIMX,KRIMY,KSV_USER, &
HTURB,HTOM,ORMC01,HRAD,HDCONV,HSCONV,HCLOUD,HELEC,HEQNSYS )
!
@@ -51,6 +51,7 @@ LOGICAL, INTENT(OUT) :: OFOREFIRE! ForeFire flag
#endif
LOGICAL, INTENT(OUT) :: OLNOX_EXPLICIT ! explicit LNOx flag
LOGICAL, INTENT(OUT) :: OCONDSAMP! Conditional sampling flag
+LOGICAL, INTENT(OUT) :: OBLOWSNOW ! Blowing snow flag
LOGICAL, INTENT(OUT) :: OORILAM ! Orilam flag
LOGICAL, INTENT(OUT) :: OCHTRANS ! Deep convection on scalar
LOGICAL,DIMENSION(JPMODELMAX),INTENT(OUT) :: ODEPOS_DST ! Dust Wet Deposition flag
@@ -85,7 +86,7 @@ END MODULE MODI_READ_DESFM_n
OFOREFIRE, &
#endif
OLNOX_EXPLICIT, &
- OCONDSAMP, &
+ OCONDSAMP,OBLOWSNOW, &
KRIMX,KRIMY,KSV_USER, &
HTURB,HTOM,ORMC01,HRAD,HDCONV,HSCONV,HCLOUD,HELEC,HEQNSYS )
! #########################################################################
@@ -190,6 +191,7 @@ END MODULE MODI_READ_DESFM_n
!! Modification 11/2016 (Ph. Wautelet) Allocate/initialise some output/backup structures
!! Modification 02/2018 (Q.Libois) ECRAD
!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!! Modification 07/2017 (V. Vionnet) Add blowing snow scheme
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -254,6 +256,8 @@ USE MODN_FOREFIRE
USE MODN_CONDSAMP
USE MODN_LATZ_EDFLX
USE MODN_2D_FRC
+USE MODN_BLOWSNOW_n
+USE MODN_BLOWSNOW
!
USE MODN_PARAM_LIMA
!
@@ -305,6 +309,7 @@ LOGICAL, INTENT(OUT) :: OFOREFIRE ! ForeFire flag
#endif
LOGICAL, INTENT(OUT) :: OLNOX_EXPLICIT ! explicit LNOx flag
LOGICAL, INTENT(OUT) :: OCONDSAMP! Conditional sampling flag
+LOGICAL, INTENT(OUT) :: OBLOWSNOW! Blowing snow flag
LOGICAL, INTENT(OUT) :: ODUST ! Dust flag
LOGICAL, INTENT(OUT) :: OORILAM ! Dust flag
LOGICAL, INTENT(OUT) :: OCHTRANS ! Deep convection on scalar
@@ -428,6 +433,13 @@ IF (GFOUND) THEN
READ(UNIT=ILUDES,NML=NAM_SERIESn)
CALL UPDATE_NAM_SERIESn
END IF
+
+CALL POSNAM(ILUDES,'NAM_BLOWSNOWn',GFOUND,ILUOUT)
+CALL INIT_NAM_BLOWSNOWn
+IF (GFOUND) THEN
+ READ(UNIT=ILUDES,NML=NAM_BLOWSNOWn)
+ CALL UPDATE_NAM_BLOWSNOWn
+END IF
!
!
IF (KMI == 1) THEN
@@ -554,6 +566,8 @@ IF (KMI == 1) THEN
#endif
CALL POSNAM(ILUDES,'NAM_CONDSAMP',GFOUND,ILUOUT)
IF (GFOUND) READ(UNIT=ILUDES,NML=NAM_CONDSAMP)
+ CALL POSNAM(ILUDES,'NAM_BLOWSNOW',GFOUND,ILUOUT)
+ IF (GFOUND) READ(UNIT=ILUDES,NML=NAM_BLOWSNOW)
CALL POSNAM(ILUDES,'NAM_2D_FRC',GFOUND,ILUOUT)
IF (GFOUND) READ(UNIT=ILUDES,NML=NAM_2D_FRC)
LTEMPDEPOS_DST(:) = LDEPOS_DST(:)
@@ -596,6 +610,9 @@ OFOREFIRE = LFOREFIRE
#endif
OLNOX_EXPLICIT = LLNOX_EXPLICIT
OCONDSAMP= LCONDSAMP
+OBLOWSNOW= LBLOWSNOW
+! Initially atmosphere free of blowing snow particles
+IF(KMI>1) OBLOWSNOW=.FALSE.
KRIMX = NRIMX
KRIMY = NRIMY
KSV_USER = NSV_USER
@@ -661,6 +678,9 @@ IF (NVERB >= 10) THEN
!
WRITE(UNIT=ILUOUT,FMT="('************ CHEMICAL SOLVER ******************')")
WRITE(UNIT=ILUOUT,NML=NAM_CH_SOLVERn)
+!
+ WRITE(UNIT=ILUOUT,FMT="('********** BLOWSNOWn *******************')")
+ WRITE(UNIT=ILUOUT,NML=NAM_BLOWSNOWn)
!
IF (KMI==1) THEN
WRITE(UNIT=ILUOUT,FMT="(/,'PART OF INITIAL FILE COMMON TO ALL THE MODELS')")
@@ -748,6 +768,9 @@ IF (NVERB >= 10) THEN
#endif
WRITE(UNIT=ILUOUT,FMT="('************ CONDITIONAL SAMPLING *************')")
WRITE(UNIT=ILUOUT,NML=NAM_CONDSAMP)
+ !
+ WRITE(UNIT=ILUOUT,FMT="('********** BLOWING SNOW SCHEME******************')")
+ WRITE(UNIT=ILUOUT,NML=NAM_BLOWSNOW)
!
IF( CCLOUD == 'C2R2' ) THEN
WRITE(UNIT=ILUOUT,FMT="('************ C2R2 SCHEME **********************')")
diff --git a/src/MNH/read_exsegn.f90 b/src/MNH/read_exsegn.f90
index 8f19edb80..09005b315 100644
--- a/src/MNH/read_exsegn.f90
+++ b/src/MNH/read_exsegn.f90
@@ -18,7 +18,7 @@ INTERFACE
OFOREFIRE, &
#endif
OLNOX_EXPLICIT, &
- OCONDSAMP, &
+ OCONDSAMP,OBLOWSNOW, &
KRIMX,KRIMY, KSV_USER, &
HTURB,HTOM,ORMC01,HRAD,HDCONV,HSCONV,HCLOUD,HELEC, &
HEQNSYS,PTSTEP_ALL,HSTORAGE_TYPE,HINIFILEPGD )
@@ -52,6 +52,7 @@ LOGICAL, INTENT(IN) :: OFOREFIRE ! ForeFire FLAG in FMFILE
#endif
LOGICAL, INTENT(IN) :: OLNOX_EXPLICIT ! explicit LNOx FLAG in FMFILE
LOGICAL, INTENT(IN) :: OCONDSAMP ! Conditional sampling FLAG in FMFILE
+LOGICAL, INTENT(IN) :: OBLOWSNOW ! Blowing snow FLAG in FMFILE
LOGICAL, INTENT(IN) :: OCHTRANS ! LCHTRANS FLAG in FMFILE
LOGICAL, INTENT(IN) :: OLG ! lagrangian FLAG in FMFILE
@@ -90,7 +91,7 @@ END MODULE MODI_READ_EXSEG_n
OFOREFIRE, &
#endif
OLNOX_EXPLICIT, &
- OCONDSAMP, &
+ OCONDSAMP, OBLOWSNOW, &
KRIMX,KRIMY, KSV_USER, &
HTURB,HTOM,ORMC01,HRAD,HDCONV,HSCONV,HCLOUD,HELEC, &
HEQNSYS,PTSTEP_ALL,HSTORAGE_TYPE,HINIFILEPGD )
@@ -286,6 +287,7 @@ END MODULE MODI_READ_EXSEG_n
!! aerosol and no cloud scheme defined
!! Q.Libois 02/2018 ECRAD
!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!! Modification 07/2017 (V. Vionnet) add blowing snow scheme
!!------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -362,11 +364,14 @@ USE MODD_FOREFIRE
USE MODN_FOREFIRE
#endif
USE MODD_CONDSAMP
+USE MODD_BLOWSNOW
USE MODN_DUST
USE MODN_SALT
USE MODD_CH_M9_n, ONLY : NEQ
USE MODN_PASPOL
USE MODN_CONDSAMP
+USE MODN_BLOWSNOW
+USE MODN_BLOWSNOW_n
USE MODN_2D_FRC
!
IMPLICIT NONE
@@ -403,6 +408,7 @@ LOGICAL, INTENT(IN) :: OFOREFIRE ! ForeFire FLAG in FMFILE
LOGICAL, INTENT(IN) :: OLNOX_EXPLICIT ! explicit LNOx FLAG in FMFILE
LOGICAL, INTENT(IN) :: OCONDSAMP ! Conditional sampling FLAG in FMFILE
LOGICAL, INTENT(IN) :: OCHTRANS ! LCHTRANS FLAG in FMFILE
+LOGICAL, INTENT(IN) :: OBLOWSNOW ! Blowing snow FLAG in FMFILE
LOGICAL, INTENT(IN) :: OLG ! lagrangian FLAG in FMFILE
INTEGER, INTENT(IN) :: KRIMX, KRIMY ! number of points for the
@@ -460,6 +466,7 @@ CALL INIT_NAM_TURBN
CALL INIT_NAM_CH_MNHCN
CALL INIT_NAM_CH_SOLVERN
CALL INIT_NAM_SERIESN
+CALL INIT_NAM_BLOWSNOWN
!
WRITE(UNIT=ILUOUT,FMT="(/,'READING THE EXSEG.NAM FILE')")
CALL POSNAM(ILUSEG,'NAM_LUNITN',GFOUND,ILUOUT)
@@ -494,6 +501,8 @@ CALL POSNAM(ILUSEG,'NAM_CH_SOLVERN',GFOUND,ILUOUT)
IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_CH_SOLVERn)
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)
!
IF (KMI == 1) THEN
WRITE(UNIT=ILUOUT,FMT="(' namelists common to all the models ')")
@@ -632,6 +641,8 @@ IF (KMI == 1) THEN
IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_2D_FRC)
CALL POSNAM(ILUSEG,'NAM_LATZ_EDFLX',GFOUND)
IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_LATZ_EDFLX)
+ CALL POSNAM(ILUSEG,'NAM_BLOWSNOW',GFOUND,ILUOUT)
+ IF (GFOUND) READ(UNIT=ILUSEG,NML=NAM_BLOWSNOW)
END IF
!
!-------------------------------------------------------------------------------
@@ -717,6 +728,20 @@ IF( CCLOUD == 'LIMA' ) THEN
CALL TEST_NAM_VAR(ILUOUT,'CHEVRIMED_ICE_LIMA',CHEVRIMED_ICE_LIMA, &
'GRAU','HAIL')
END IF
+IF(LBLOWSNOW) THEN
+ CALL TEST_NAM_VAR(ILUOUT,'CSNOWSEDIM',CSNOWSEDIM,'NONE','MITC','CARR','TABC')
+ IF (XALPHA_SNOW .NE. 3 .AND. CSNOWSEDIM=='TABC') THEN
+ WRITE(ILUOUT,*) '*****************************************'
+ WRITE(ILUOUT,*) '* XALPHA_SNW must be set to 3 when '
+ WRITE(ILUOUT,*) '* CSNOWSEDIM = TABC '
+ WRITE(ILUOUT,*) '* Update the look-up table in BLOWSNOW_SEDIM_LKT1D '
+ WRITE(ILUOUT,*) '* to use TABC with a different value of XEMIALPHA_SNW'
+ WRITE(ILUOUT,*) '*****************************************'
+ !callabortstop
+ CALL PRINT_MSG(NVERB_FATAL,'GEN','READ_EXSEG_n','')
+ ENDIF
+END IF
+!
!-------------------------------------------------------------------------------!
!* 2. FIRST INITIALIZATIONS
! ---------------------
@@ -1902,6 +1927,28 @@ IF (LCONDSAMP) THEN
END IF
END IF
!
+! Blowing snow scheme
+!
+IF (LBLOWSNOW) THEN
+ IF (OBLOWSNOW) THEN
+ CGETSVT(NSV_SNWBEG:NSV_SNWEND)='READ'
+ ELSE
+ WRITE(UNIT=ILUOUT,FMT=9001) KMI
+ WRITE(UNIT=ILUOUT,FMT='("THERE IS NO SCALAR VARIABLES FOR BLOWING SNOW &
+ &SCHEME IN INITIAL FMFILE",/,&
+ & "THE BLOWING SNOW VARIABLES HAVE BEEN INITIALIZED TO ZERO ")')
+ CGETSVT(NSV_SNWBEG:NSV_SNWEND)='INIT'
+ END IF
+ IF(.NOT.ALLOCATED(CSNOWNAMES)) THEN
+ IMOMENTS = (NSV_SNWEND - NSV_SNWBEG +1 )
+ ALLOCATE(CSNOWNAMES(IMOMENTS))
+ DO JMOM=1,IMOMENTS
+ CSNOWNAMES(JMOM) = YPSNOW_INI(JMOM)
+ ENDDO ! Loop on moments
+ END IF
+END IF
+!
+!
!
!* 3.5 Check coherence between the radiation control parameters
!
@@ -2381,6 +2428,13 @@ IF (.NOT. LCONDSAMP .AND. LHORELAX_SVCS) THEN
WRITE(ILUOUT,FMT=*) 'THEREFORE LHORELAX_SVCS=FALSE'
END IF
+IF (.NOT. LBLOWSNOW .AND. LHORELAX_SVSNW) THEN
+ LHORELAX_SVSNW=.FALSE.
+ WRITE(UNIT=ILUOUT,FMT=9002) KMI
+ WRITE(ILUOUT,FMT=*) 'YOU WANT TO RELAX BLOWING SNOW FIELD BUT IT DOES NOT EXIST.'
+ WRITE(ILUOUT,FMT=*) 'THEREFORE LHORELAX_SVSNW=FALSE'
+END IF
+
IF (ANY(LHORELAX_SV(NSV+1:))) THEN
LHORELAX_SV(NSV+1:)=.FALSE.
WRITE(UNIT=ILUOUT,FMT=9002) KMI
@@ -2752,6 +2806,7 @@ CALL UPDATE_NAM_TURBN
CALL UPDATE_NAM_CH_MNHCN
CALL UPDATE_NAM_CH_SOLVERN
CALL UPDATE_NAM_SERIESN
+CALL UPDATE_NAM_BLOWSNOWN
!-------------------------------------------------------------------------------
WRITE(UNIT=ILUOUT,FMT='(/)')
!-------------------------------------------------------------------------------
diff --git a/src/MNH/read_field.f90 b/src/MNH/read_field.f90
index b8e0e4f69..12e8b2194 100644
--- a/src/MNH/read_field.f90
+++ b/src/MNH/read_field.f90
@@ -238,6 +238,7 @@ END MODULE MODI_READ_FIELD
!! C.Lac 10/16 CEN4TH with RKC4 + Correction on RK loop
!! 09/2017 Q.Rodier add LTEND_UV_FRC
!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!! V. Vionnet 07/17 add blowing snow scheme
!!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -255,6 +256,9 @@ USE MODD_ELEC_DESCR, ONLY: CELECNAMES
#ifdef MNH_FOREFIRE
USE MODD_FOREFIRE
#endif
+USE MODD_BLOWSNOW
+USE MODD_BLOWSNOW_n
+
USE MODD_ICE_C1R3_DESCR, ONLY: C1R3NAMES
USE MODD_IO_ll, ONLY: TFILEDATA
USE MODD_LATZ_EDFLX
@@ -1070,6 +1074,49 @@ IF (NSV_LNOXEND>=NSV_LNOXBEG) THEN
END SELECT
END DO
END IF
+!
+IF (NSV_SNWEND>=NSV_SNWBEG) THEN
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CUNITS = 'kg/kg'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%NGRID = 1
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 3
+ TZFIELD%LTIMEDEP = .TRUE.
+ DO JSV = NSV_SNWBEG,NSV_SNWEND
+ SELECT CASE(HGETSVT(JSV))
+ CASE ('READ')
+ TZFIELD%CMNHNAME = TRIM(CSNOWNAMES(JSV-NSV_SNWBEG+1))//'T'
+ TZFIELD%CLONGNAME = TRIM(TZFIELD%CMNHNAME)
+ WRITE(TZFIELD%CCOMMENT,'(A6,A3,I3.3)') 'X_Y_Z_','SVT',JSV
+ CALL IO_READ_FIELD(TPINIFILE,TZFIELD,PSVT(:,:,:,JSV))
+ CASE ('INIT')
+ PSVT(:,:,:,JSV) = 0.
+ END SELECT
+ END DO
+END IF
+IF (NSV_SNW>=1) THEN
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CUNITS = 'kg/kg'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%NGRID = 1
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 2
+ TZFIELD%LTIMEDEP = .TRUE.
+ DO JSV = 1,NSV_SNW
+ SELECT CASE(HGETSVT(JSV))
+ CASE ('READ')
+ WRITE(TZFIELD%CMNHNAME,'(A10,I3.3)')'SNOWCANO_M',JSV
+ TZFIELD%CLONGNAME = TRIM(TZFIELD%CMNHNAME)
+ WRITE(TZFIELD%CCOMMENT,'(A6,A8,I3.3)') 'X_Y_Z_','SNOWCANO',JSV
+ CALL IO_READ_FIELD(TPINIFILE,TZFIELD,XSNWCANO(:,:,JSV))
+ CASE ('INIT')
+ XSNWCANO(:,:,JSV) = 0.
+ END SELECT
+ END DO
+
+END IF
+
!
IF (CCONF == 'RESTA') THEN
CALL IO_READ_FIELD(TPINIFILE,'US_PRES',PRUS_PRES)
diff --git a/src/MNH/relaxation.f90 b/src/MNH/relaxation.f90
index 423b77e89..13e62a4ab 100644
--- a/src/MNH/relaxation.f90
+++ b/src/MNH/relaxation.f90
@@ -16,7 +16,7 @@ INTERFACE
OHORELAX_SVELEC,OHORELAX_SVLG, &
OHORELAX_SVCHEM,OHORELAX_SVCHIC, OHORELAX_SVAER, &
OHORELAX_SVDST, OHORELAX_SVSLT, OHORELAX_SVPP, &
- OHORELAX_SVCS, &
+ OHORELAX_SVCS,OHORELAX_SVSNW, &
#ifdef MNH_FOREFIRE
OHORELAX_SVFF, &
#endif
@@ -83,6 +83,8 @@ LOGICAL, INTENT(IN):: OHORELAX_SVFF ! switch for the
#endif
LOGICAL, INTENT(IN):: OHORELAX_SVCS ! switch for the
! horizontal relaxation for conditional sampling
+LOGICAL, INTENT(IN):: OHORELAX_SVSNW ! switch for the
+ ! horizontal relaxation for blowing snow variables
INTEGER, INTENT(IN) :: KTCOUNT! Temporal loop counter
REAL, INTENT(IN) :: PTSTEP ! Time step
INTEGER, INTENT(IN) :: KRR ! Number of moist variables
@@ -156,7 +158,7 @@ END MODULE MODI_RELAXATION
OHORELAX_SVELEC,OHORELAX_SVLG, &
OHORELAX_SVCHEM,OHORELAX_SVCHIC, OHORELAX_SVAER, &
OHORELAX_SVDST, OHORELAX_SVSLT, OHORELAX_SVPP, &
- OHORELAX_SVCS, &
+ OHORELAX_SVCS,OHORELAX_SVSNW, &
#ifdef MNH_FOREFIRE
OHORELAX_SVFF, &
#endif
@@ -329,6 +331,8 @@ LOGICAL, INTENT(IN):: OHORELAX_SVFF ! switch for the
#endif
LOGICAL, INTENT(IN):: OHORELAX_SVCS ! switch for the
! horizontal relaxation for conditional sampling
+LOGICAL, INTENT(IN):: OHORELAX_SVSNW ! switch for the
+ ! horizontal relaxation for blowing snow variables
INTEGER, INTENT(IN) :: KTCOUNT! Temporal loop counter
REAL, INTENT(IN) :: PTSTEP ! Time step
INTEGER, INTENT(IN) :: KRR ! Number of moist variables
@@ -419,9 +423,9 @@ LOGICAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: GMASK3D_RELAX ! 3D
! mask for hor. relax.
LOGICAL, DIMENSION(7) :: GHORELAXR ! local array of logical
#ifdef MNH_FOREFIRE
-LOGICAL, DIMENSION(12) :: GHORELAXSV! local array of logical
+LOGICAL, DIMENSION(13) :: GHORELAXSV! local array of logical
#else
-LOGICAL, DIMENSION(11) :: GHORELAXSV! local array of logical
+LOGICAL, DIMENSION(12) :: GHORELAXSV! local array of logical
#endif
!
!-------------------------------------------------------------------------------
@@ -459,8 +463,9 @@ GHORELAXSV(8) = OHORELAX_SVSLT
GHORELAXSV(9) = OHORELAX_SVPP
GHORELAXSV(10) = OHORELAX_SVCS
GHORELAXSV(11) = OHORELAX_SVCHIC
+GHORELAXSV(12) = OHORELAX_SVSNW
#ifdef MNH_FOREFIRE
-GHORELAXSV(12) = OHORELAX_SVFF
+GHORELAXSV(13) = OHORELAX_SVFF
#endif
!-------------------------------------------------------------------------------
!
diff --git a/src/MNH/relaxdef.f90 b/src/MNH/relaxdef.f90
index c87a5fe6c..d4579df48 100644
--- a/src/MNH/relaxdef.f90
+++ b/src/MNH/relaxdef.f90
@@ -14,7 +14,7 @@ INTERFACE
OHORELAX_RH,OHORELAX_TKE,OHORELAX_SV, &
OHORELAX_SVC2R2,OHORELAX_SVC1R3,OHORELAX_SVELEC,OHORELAX_SVLG, &
OHORELAX_SVCHEM, OHORELAX_SVAER, OHORELAX_SVDST, OHORELAX_SVSLT, &
- OHORELAX_SVPP,OHORELAX_SVCS, OHORELAX_SVCHIC, &
+ OHORELAX_SVPP,OHORELAX_SVCS, OHORELAX_SVCHIC,OHORELAX_SVSNW, &
PALKTOP,PALKGRD, PALZBOT,PALZBAS, &
PZZ, PZHAT, PTSTEP, &
PRIMKMAX,KRIMX, KRIMY, &
@@ -67,6 +67,8 @@ LOGICAL, INTENT(IN):: OHORELAX_SVPP ! switch for the
! horizontal relaxation for passive pollutants
LOGICAL, INTENT(IN):: OHORELAX_SVCS ! switch for the
! horizontal relaxation for conditional sampling
+LOGICAL, INTENT(IN):: OHORELAX_SVSNW ! switch for the
+ ! horizontal relaxation for blowing snow variables
REAL, INTENT(IN) :: PALKTOP ! Damping coef. at the
! top of the abs. layer
REAL, INTENT(IN) :: PALKGRD ! Damping coef. at the
@@ -117,7 +119,7 @@ END MODULE MODI_RELAXDEF
OHORELAX_RH,OHORELAX_TKE,OHORELAX_SV, &
OHORELAX_SVC2R2,OHORELAX_SVC1R3,OHORELAX_SVELEC,OHORELAX_SVLG, &
OHORELAX_SVCHEM, OHORELAX_SVAER, OHORELAX_SVDST, OHORELAX_SVSLT, &
- OHORELAX_SVPP,OHORELAX_SVCS, OHORELAX_SVCHIC, &
+ OHORELAX_SVPP,OHORELAX_SVCS, OHORELAX_SVCHIC,OHORELAX_SVSNW, &
PALKTOP,PALKGRD, PALZBOT,PALZBAS, &
PZZ, PZHAT, PTSTEP, &
PRIMKMAX,KRIMX, KRIMY, &
@@ -290,6 +292,8 @@ LOGICAL, INTENT(IN):: OHORELAX_SVPP ! switch for the
! horizontal relaxation for passive pollutants
LOGICAL, INTENT(IN):: OHORELAX_SVCS ! switch for the
! horizontal relaxation for conditional sampling
+LOGICAL, INTENT(IN):: OHORELAX_SVSNW ! switch for the
+ ! horizontal relaxation for blowing snow variables
REAL, INTENT(IN) :: PALKTOP ! Damping coef. at the
! top of the abs. layer
REAL, INTENT(IN) :: PALKGRD ! Damping coef. at the
@@ -380,7 +384,7 @@ REAL :: ZYWPOS ! lateral boundary for "w" locations
REAL :: ZPOS ! Normalized distance
! to the inner boundary
LOGICAL, DIMENSION(7) :: GHORELAXR ! local array of logical
-LOGICAL, DIMENSION(11) :: GHORELAXSV ! local array of logical
+LOGICAL, DIMENSION(12) :: GHORELAXSV ! local array of logical
INTEGER :: IINFO_ll ! return code of parallel routine
INTEGER :: IIMAX_ll,IJMAX_ll ! Number of points of
! Global physical domain
@@ -493,6 +497,7 @@ GHORELAXSV(8) = OHORELAX_SVSLT
GHORELAXSV(9) = OHORELAX_SVPP
GHORELAXSV(10)= OHORELAX_SVCS
GHORELAXSV(11) = OHORELAX_SVCHIC
+GHORELAXSV(12) = OHORELAX_SVSNW
!
IF ( ANY(GHORELAXR) .OR. ANY(GHORELAXSV) .OR. ANY(OHORELAX_SV) &
.OR. OHORELAX_UVWTH .OR. OHORELAX_TKE ) THEN
diff --git a/src/MNH/turb_hor_sv_corr.f90 b/src/MNH/turb_hor_sv_corr.f90
index 715f6df35..77489e8c6 100644
--- a/src/MNH/turb_hor_sv_corr.f90
+++ b/src/MNH/turb_hor_sv_corr.f90
@@ -95,6 +95,7 @@ USE MODD_CTURB
USE MODD_PARAMETERS
USE MODD_NSV, ONLY : NSV,NSV_LGBEG,NSV_LGEND
USE MODD_LES
+USE MODD_BLOWSNOW
!
USE MODI_GRADIENT_M
USE MODI_GRADIENT_U
@@ -148,11 +149,20 @@ REAL :: ZTIME1, ZTIME2
REAL :: ZCSVD = 1.2 ! constant for scalar variance dissipation
REAL :: ZCTSVD = 2.4 ! constant for temperature - scalar covariance dissipation
REAL :: ZCQSVD = 2.4 ! constant for humidity - scalar covariance dissipation
+!
+REAL :: ZCSV !constant for the scalar flux
! ---------------------------------------------------------------------------
!
IKU=SIZE(PTKEM,3)
CALL SECOND_MNH(ZTIME1)
!
+IF(LBLOWSNOW) THEN
+! See Vionnet (PhD, 2012) for a complete discussion around the value of the Schmidt number for blowing snow variables
+ ZCSV= XCHF/XRSNOW
+ELSE
+ ZCSV= XCHF
+ENDIF
+!
DO JSV=1,NSV
!
IF (LNOMIXLG .AND. JSV >= NSV_LGBEG .AND. JSV<= NSV_LGEND) CYCLE
@@ -161,11 +171,11 @@ DO JSV=1,NSV
!
IF (LLES_CALL) THEN
IF (.NOT. L2D) THEN
- ZFLX(:,:,:) = XCHF / ZCSVD * PLM(:,:,:) * PLEPS(:,:,:) * &
+ ZFLX(:,:,:) = ZCSV / ZCSVD * PLM(:,:,:) * PLEPS(:,:,:) * &
( GX_M_M(1,IKU,1,PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX)**2 &
+ GY_M_M(1,IKU,1,PSVM(:,:,:,JSV),PDYY,PDZZ,PDZY)**2 )
ELSE
- ZFLX(:,:,:) = XCHF / ZCSVD * PLM(:,:,:) * PLEPS(:,:,:) * &
+ ZFLX(:,:,:) = ZCSV / ZCSVD * PLM(:,:,:) * PLEPS(:,:,:) * &
GX_M_M(1,IKU,1,PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX)**2
END IF
CALL LES_MEAN_SUBGRID( -2.*ZCSVD*SQRT(PTKEM)*ZFLX/PLEPS, &
@@ -181,11 +191,11 @@ DO JSV=1,NSV
ZFLX(:,:,:)= PLM(:,:,:) * PLEPS(:,:,:) &
* ( GX_M_M(1,IKU,1,PTHLM,PDXX,PDZZ,PDZX) * GX_M_M(1,IKU,1,PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX) &
+ GY_M_M(1,IKU,1,PTHLM,PDYY,PDZZ,PDZY) * GY_M_M(1,IKU,1,PSVM(:,:,:,JSV),PDYY,PDZZ,PDZY) &
- ) * (XCSHF+XCHF) / (2.*ZCTSVD)
+ ) * (XCSHF+ZCSV) / (2.*ZCTSVD)
ELSE
ZFLX(:,:,:)= PLM(:,:,:) * PLEPS(:,:,:) &
* GX_M_M(1,IKU,1,PTHLM,PDXX,PDZZ,PDZX) * GX_M_M(1,IKU,1,PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX) &
- * (XCSHF+XCHF) / (2.*ZCTSVD)
+ * (XCSHF+ZCSV) / (2.*ZCTSVD)
END IF
CALL LES_MEAN_SUBGRID( ZA*ZFLX, X_LES_SUBGRID_SvThv(:,:,:,JSV) , .TRUE.)
CALL LES_MEAN_SUBGRID( -XG/PTHVREF/3.*ZA*ZFLX, X_LES_SUBGRID_SvPz(:,:,:,JSV), .TRUE. )
@@ -196,11 +206,11 @@ DO JSV=1,NSV
ZFLX(:,:,:)= PLM(:,:,:) * PLEPS(:,:,:) &
* ( GX_M_M(1,IKU,1,PRM(:,:,:,1),PDXX,PDZZ,PDZX) * GX_M_M(1,IKU,1,PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX) &
+ GY_M_M(1,IKU,1,PRM(:,:,:,1),PDYY,PDZZ,PDZY) * GY_M_M(1,IKU,1,PSVM(:,:,:,JSV),PDYY,PDZZ,PDZY) &
- ) * (XCHF+XCHF) / (2.*ZCQSVD)
+ ) * (XCHF+ZCSV) / (2.*ZCQSVD)
ELSE
ZFLX(:,:,:)= PLM(:,:,:) * PLEPS(:,:,:) &
* GX_M_M(1,IKU,1,PRM(:,:,:,1),PDXX,PDZZ,PDZX) * GX_M_M(1,IKU,1,PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX) &
- * (XCHF+XCHF) / (2.*ZCQSVD)
+ * (XCHF+ZCSV) / (2.*ZCQSVD)
END IF
CALL LES_MEAN_SUBGRID( ZA*ZFLX, X_LES_SUBGRID_SvThv(:,:,:,JSV) , .TRUE.)
CALL LES_MEAN_SUBGRID( -XG/PTHVREF/3.*ZA*ZFLX, X_LES_SUBGRID_SvPz(:,:,:,JSV), .TRUE. )
diff --git a/src/MNH/turb_hor_sv_flux.f90 b/src/MNH/turb_hor_sv_flux.f90
index 55f2b3947..acb3da871 100644
--- a/src/MNH/turb_hor_sv_flux.f90
+++ b/src/MNH/turb_hor_sv_flux.f90
@@ -119,6 +119,7 @@ USE MODD_IO_ll, ONLY: TFILEDATA
USE MODD_PARAMETERS
USE MODD_NSV, ONLY : NSV_LGBEG,NSV_LGEND
USE MODD_LES
+USE MODD_BLOWSNOW
!
USE MODE_FIELD, ONLY: TFIELDDATA,TYPEREAL
USE MODE_FMWRIT
@@ -177,6 +178,8 @@ REAL, DIMENSION(SIZE(PSVM,1),SIZE(PSVM,2),SIZE(PSVM,3)) &
! work arrays
REAL, DIMENSION(SIZE(PSVM,1),SIZE(PSVM,2),1) :: ZWORK2D
!
+REAL :: ZCSV !constant for the scalar flux
+
INTEGER :: IKB,IKE
! Index values for the Beginning and End
! mass points of the domain
@@ -200,6 +203,12 @@ IKU = SIZE(PSVM,3)
!
ISV = SIZE(PSVM,4)
!
+IF(LBLOWSNOW) THEN
+! See Vionnet (PhD, 2012) for a complete discussion around the value of the Schmidt number for blowing snow variables
+ ZCSV= XCHF/XRSNOW
+ELSE
+ ZCSV= XCHF
+ENDIF
!
! compute the coefficients for the uncentred gradient computation near the
! ground
@@ -223,12 +232,12 @@ DO JSV=1,ISV
! ----------
!
! Computes the flux in the X direction
- ZFLXX(:,:,:) = -XCHF * MXM(PK) * GX_M_U(1,IKU,1,PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX)
+ ZFLXX(:,:,:) = -ZCSV * MXM(PK) * GX_M_U(1,IKU,1,PSVM(:,:,:,JSV),PDXX,PDZZ,PDZX)
ZFLXX(:,:,IKE+1) = ZFLXX(:,:,IKE)
!
! Compute the flux at the first inner U-point with an uncentred vertical
! gradient
- ZFLXX(:,:,IKB:IKB) = -XCHF * MXM( PK(:,:,IKB:IKB) ) * &
+ ZFLXX(:,:,IKB:IKB) = -ZCSV * MXM( PK(:,:,IKB:IKB) ) * &
( DXM(PSVM(:,:,IKB:IKB,JSV)) * PINV_PDXX(:,:,IKB:IKB) &
-MXM ( ZCOEFF(:,:,IKB+2:IKB+2)*PSVM(:,:,IKB+2:IKB+2,JSV) &
+ZCOEFF(:,:,IKB+1:IKB+1)*PSVM(:,:,IKB+1:IKB+1,JSV) &
@@ -273,13 +282,13 @@ DO JSV=1,ISV
IF (.NOT. L2D) THEN
!
! Computes the flux in the Y direction
- ZFLXY(:,:,:)=-XCHF * MYM(PK) * GY_M_V(1,IKU,1,PSVM(:,:,:,JSV),PDYY,PDZZ,PDZY)
+ ZFLXY(:,:,:)=-ZCSV * MYM(PK) * GY_M_V(1,IKU,1,PSVM(:,:,:,JSV),PDYY,PDZZ,PDZY)
ZFLXY(:,:,IKE+1) = ZFLXY(:,:,IKE)
!
! Compute the flux at the first inner V-point with an uncentred vertical
! gradient
!
- ZFLXY(:,:,IKB:IKB) = -XCHF * MYM( PK(:,:,IKB:IKB) ) * &
+ ZFLXY(:,:,IKB:IKB) = -ZCSV * MYM( PK(:,:,IKB:IKB) ) * &
( DYM(PSVM(:,:,IKB:IKB,JSV)) * PINV_PDYY(:,:,IKB:IKB) &
-MYM ( ZCOEFF(:,:,IKB+2:IKB+2)*PSVM(:,:,IKB+2:IKB+2,JSV) &
+ZCOEFF(:,:,IKB+1:IKB+1)*PSVM(:,:,IKB+1:IKB+1,JSV) &
diff --git a/src/MNH/turb_ver_sv_corr.f90 b/src/MNH/turb_ver_sv_corr.f90
index 75f632213..280f94392 100644
--- a/src/MNH/turb_ver_sv_corr.f90
+++ b/src/MNH/turb_ver_sv_corr.f90
@@ -108,6 +108,7 @@ USE MODD_PARAMETERS
USE MODD_LES
USE MODD_CONF
USE MODD_NSV, ONLY : NSV,NSV_LGBEG,NSV_LGEND
+USE MODD_BLOWSNOW
!
!
USE MODI_GRADIENT_U
@@ -160,6 +161,9 @@ REAL, DIMENSION(:,:,:,:), INTENT(IN) :: PPSI_SV ! Inv.Turb.Sch.for scalars
!
REAL, DIMENSION(SIZE(PSVM,1),SIZE(PSVM,2),SIZE(PSVM,3)) :: &
ZA, ZFLXZ
+!
+REAL :: ZCSV !constant for the scalar flux
+!
INTEGER :: JSV ! loop counters
!
REAL :: ZTIME1, ZTIME2
@@ -171,6 +175,13 @@ REAL :: ZCQSVD = 2.4 ! constant for humidity - scalar covariance dissipation
!
CALL SECOND_MNH(ZTIME1)
!
+IF(LBLOWSNOW) THEN
+! See Vionnet (PhD, 2012) for a complete discussion around the value of the Schmidt number for blowing snow variables
+ ZCSV= XCHF/XRSNOW
+ELSE
+ ZCSV= XCHF
+ENDIF
+!
DO JSV=1,NSV
!
IF (LNOMIXLG .AND. JSV >= NSV_LGBEG .AND. JSV<= NSV_LGEND) CYCLE
@@ -180,7 +191,7 @@ DO JSV=1,NSV
IF (LLES_CALL) THEN
! approximation: diagnosed explicitely (without implicit term)
ZFLXZ(:,:,:) = PPSI_SV(:,:,:,JSV)*GZ_M_W(KKA,KKU,KKL,PSVM(:,:,:,JSV),PDZZ)**2
- ZFLXZ(:,:,:) = XCHF / ZCSVD * PLM * PLEPS * MZF(KKA,KKU,KKL,ZFLXZ(:,:,:) )
+ ZFLXZ(:,:,:) = ZCSV / ZCSVD * PLM * PLEPS * MZF(KKA,KKU,KKL,ZFLXZ(:,:,:) )
CALL LES_MEAN_SUBGRID( -2.*ZCSVD*SQRT(PTKEM)*ZFLXZ/PLEPS, X_LES_SUBGRID_DISS_Sv2(:,:,:,JSV) )
CALL LES_MEAN_SUBGRID( MZF(KKA,KKU,KKL,PWM)*ZFLXZ, X_LES_RES_W_SBG_Sv2(:,:,:,JSV) )
END IF
@@ -190,7 +201,7 @@ DO JSV=1,NSV
IF (LLES_CALL) THEN
! approximation: diagnosed explicitely (without implicit term)
ZA(:,:,:) = ETHETA(KRR,KRRI,PTHLM,PRM,PLOCPEXNM,PATHETA,PSRCM)
- ZFLXZ(:,:,:)= ( XCSHF * PPHI3 + XCHF * PPSI_SV(:,:,:,JSV) ) &
+ ZFLXZ(:,:,:)= ( XCSHF * PPHI3 + ZCSV * PPSI_SV(:,:,:,JSV) ) &
* GZ_M_W(KKA,KKU,KKL,PTHLM,PDZZ) &
* GZ_M_W(KKA,KKU,KKL,PSVM(:,:,:,JSV),PDZZ)
ZFLXZ(:,:,:)= PLM * PLEPS / (2.*ZCTSVD) * MZF(KKA,KKU,KKL,ZFLXZ)
@@ -199,7 +210,7 @@ DO JSV=1,NSV
!
IF (KRR>=1) THEN
ZA(:,:,:) = EMOIST(KRR,KRRI,PTHLM,PRM,PLOCPEXNM,PAMOIST,PSRCM)
- ZFLXZ(:,:,:)= ( XCHF * PPSI3 + XCHF * PPSI_SV(:,:,:,JSV) ) &
+ ZFLXZ(:,:,:)= ( XCHF * PPSI3 + ZCSV * PPSI_SV(:,:,:,JSV) ) &
* GZ_M_W(KKA,KKU,KKL,PRM(:,:,:,1),PDZZ) &
* GZ_M_W(KKA,KKU,KKL,PSVM(:,:,:,JSV),PDZZ)
ZFLXZ(:,:,:)= PLM * PLEPS / (2.*ZCQSVD) * MZF(KKA,KKU,KKL,ZFLXZ)
diff --git a/src/MNH/turb_ver_sv_flux.f90 b/src/MNH/turb_ver_sv_flux.f90
index 21c4bc13f..26b102f73 100644
--- a/src/MNH/turb_ver_sv_flux.f90
+++ b/src/MNH/turb_ver_sv_flux.f90
@@ -277,7 +277,7 @@ USE MODD_PARAMETERS
USE MODD_LES
USE MODD_CONF
USE MODD_NSV, ONLY : XSVMIN,NSV_LGBEG,NSV_LGEND
-!
+USE MODD_BLOWSNOW
USE MODE_FIELD, ONLY: TFIELDDATA,TYPEREAL
USE MODE_FMWRIT
!
@@ -361,6 +361,8 @@ INTEGER :: ISV ! number of scalar var.
REAL :: ZTIME1, ZTIME2
REAL :: ZCSVP = 4.0 ! constant for scalar flux presso-correlation (RS81)
+REAL :: ZCSV !constant for the scalar flux
+!
TYPE(TFIELDDATA) :: TZFIELD
!----------------------------------------------------------------------------
!
@@ -377,6 +379,12 @@ ISV=SIZE(PSVM,4)
!
ZKEFF(:,:,:) = MZM(KKA,KKU,KKL, PLM(:,:,:) * SQRT(PTKEM(:,:,:)) )
!
+IF(LBLOWSNOW) THEN
+! See Vionnet (PhD, 2012) for a complete discussion around the value of the Schmidt number for blowing snow variables
+ ZCSV= XCHF/XRSNOW
+ELSE
+ ZCSV= XCHF
+ENDIF
!----------------------------------------------------------------------------
!
!* 8. SOURCES OF PASSIVE SCALAR VARIABLES
@@ -387,7 +395,7 @@ DO JSV=1,ISV
IF (LNOMIXLG .AND. JSV >= NSV_LGBEG .AND. JSV<= NSV_LGEND) CYCLE
!
! Preparation of the arguments for TRIDIAG
- ZA(:,:,:) = -PTSTEP*XCHF*PPSI_SV(:,:,:,JSV) * &
+ ZA(:,:,:) = -PTSTEP*ZCSV*PPSI_SV(:,:,:,JSV) * &
ZKEFF * MZM(KKA,KKU,KKL,PRHODJ) / &
PDZZ**2
ZSOURCE(:,:,:) = 0.
@@ -423,7 +431,7 @@ DO JSV=1,ISV
IF ( (OTURB_FLX .AND. OCLOSE_OUT) .OR. LLES_CALL ) THEN
! Diagnostic of the cartesian vertical flux
!
- ZFLXZ(:,:,:) = -XCHF * PPSI_SV(:,:,:,JSV) * MZM(KKA,KKU,KKL,PLM*SQRT(PTKEM)) / PDZZ * &
+ ZFLXZ(:,:,:) = -ZCSV * PPSI_SV(:,:,:,JSV) * MZM(KKA,KKU,KKL,PLM*SQRT(PTKEM)) / PDZZ * &
DZM(KKA,KKU,KKL, PIMPL*ZRES(:,:,:) + PEXPL*PSVM(:,:,:,JSV) )
! surface flux
!* in 3DIM case, a part of the flux goes vertically, and another goes horizontally
diff --git a/src/MNH/update_nsv.f90 b/src/MNH/update_nsv.f90
index e5555891f..0677a414a 100644
--- a/src/MNH/update_nsv.f90
+++ b/src/MNH/update_nsv.f90
@@ -23,6 +23,7 @@ END MODULE MODI_UPDATE_NSV
!! the NSV_* variables.
!! Modify (Escobar ) 2/2014 : add Forefire var
!! Modify (Vie) 2016 : add LIMA
+!! V. Vionnet 7/2017 : add blowing snow var
USE MODD_CONF, ONLY : NVERB
USE MODD_NSV
IMPLICIT NONE
@@ -110,6 +111,9 @@ NSV_FFEND = NSV_FFEND_A(KMI)
NSV_CS = NSV_CS_A(KMI)
NSV_CSBEG = NSV_CSBEG_A(KMI)
NSV_CSEND = NSV_CSEND_A(KMI)
+NSV_SNW = NSV_SNW_A(KMI)
+NSV_SNWBEG = NSV_SNWBEG_A(KMI)
+NSV_SNWEND = NSV_SNWEND_A(KMI)
!
END SUBROUTINE UPDATE_NSV
diff --git a/src/MNH/write_desfmn.f90 b/src/MNH/write_desfmn.f90
index 0711827fb..223153f59 100644
--- a/src/MNH/write_desfmn.f90
+++ b/src/MNH/write_desfmn.f90
@@ -143,6 +143,7 @@ END MODULE MODI_WRITE_DESFM_n
!! Modification 01/2016 (JP Pinty) Add LIMA
!! 02/2018 Q.Libois ECRAD
!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!! Modification V. Vionnet 07/2017 add blowing snow variables
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
@@ -198,6 +199,8 @@ USE MODN_LATZ_EDFLX
USE MODN_FOREFIRE
USE MODD_FOREFIRE_n, ONLY : FFCOUPLING
#endif
+USE MODN_BLOWSNOW_n
+USE MODN_BLOWSNOW
!
IMPLICIT NONE
!
@@ -220,7 +223,7 @@ LOGICAL :: GHORELAX_UVWTH, &
GHORELAX_SVFF, &
#endif
GHORELAX_SVCHEM, GHORELAX_SVC1R3, &
- GHORELAX_SVELEC, GHORELAX_SVLIMA
+ GHORELAX_SVELEC, GHORELAX_SVLIMA,GHORELAX_SVSNW
LOGICAL :: GHORELAX_SVDST, GHORELAX_SVSLT, GHORELAX_SVAER
LOGICAL, DIMENSION(JPSVMAX) :: GHORELAX_SV
!
@@ -271,6 +274,7 @@ IF (CPROGRAM/='MESONH') THEN ! impose default value for next simulation
#endif
GHORELAX_SVCS = LHORELAX_SVCS
GHORELAX_SVAER = LHORELAX_SVAER
+ GHORELAX_SVSNW = LHORELAX_SVSNW
!
LHORELAX_UVWTH = .FALSE.
LHORELAX_RV = .FALSE.
@@ -296,6 +300,7 @@ IF (CPROGRAM/='MESONH') THEN ! impose default value for next simulation
LHORELAX_SVDST= .FALSE.
LHORELAX_SVSLT= .FALSE.
LHORELAX_SVAER= .FALSE.
+ LHORELAX_SVSNW= .FALSE.
ELSE !return to namelist meaning of LHORELAX_SV
GHORELAX_SV(:) = LHORELAX_SV(:)
LHORELAX_SV(NSV_USER+1:)=.FALSE.
@@ -357,6 +362,10 @@ IF(LUSECHEM .OR. LCH_CONV_LINOX .OR. LCH_CONV_SCAV) &
CALL INIT_NAM_CH_SOLVERn
IF(LUSECHEM) WRITE(UNIT=ILUSEG,NML=NAM_CH_SOLVERn)
!
+CALL INIT_NAM_BLOWSNOWn
+IF(LBLOWSNOW) WRITE(UNIT=ILUSEG,NML=NAM_BLOWSNOWn)
+IF(LBLOWSNOW) WRITE(UNIT=ILUSEG,NML=NAM_BLOWSNOW)
+!
IF(LDUST) WRITE(UNIT=ILUSEG,NML=NAM_DUST)
IF(LSALT) WRITE(UNIT=ILUSEG,NML=NAM_SALT)
IF(LPASPOL) WRITE(UNIT=ILUSEG,NML=NAM_PASPOL)
@@ -472,6 +481,9 @@ IF (NVERB >= 5) THEN
WRITE(UNIT=ILUOUT,FMT="('************ TEMPORAL SERIESn ******************')")
WRITE(UNIT=ILUOUT,NML=NAM_SERIESn)
!
+ WRITE(UNIT=ILUOUT,FMT="('********** BLOWING SNOW SCHEME ****************')")
+ WRITE(UNIT=ILUOUT,NML=NAM_BLOWSNOWn)
+!
IF (KMI==1) THEN
WRITE(UNIT=ILUOUT,FMT="(/,'PART OF SEGMENT FILE COMMON TO ALL THE MODELS')")
WRITE(UNIT=ILUOUT,FMT="( '---------------------------------------------')")
@@ -561,6 +573,9 @@ IF (NVERB >= 5) THEN
!
WRITE(UNIT=ILUOUT,FMT="('********** SALT SCHEME ************************')")
WRITE(UNIT=ILUOUT,NML=NAM_SALT)
+!
+ WRITE(UNIT=ILUOUT,FMT="('********** BLOWING SNOW SCHEME ****************')")
+ WRITE(UNIT=ILUOUT,NML=NAM_BLOWSNOW)
!
WRITE(UNIT=ILUOUT,FMT="('************ ORILAM SCHEME ********************')")
WRITE(UNIT=ILUOUT,NML=NAM_CH_ORILAM)
@@ -624,6 +639,7 @@ IF (CPROGRAM /='MESONH') THEN !return to previous LHORELAX_
#endif
LHORELAX_SVCS = GHORELAX_SVCS
LHORELAX_SVAER = GHORELAX_SVAER
+ LHORELAX_SVSNW = GHORELAX_SVSNW
ELSE
LHORELAX_SV(:) = GHORELAX_SV(:)
ENDIF
diff --git a/src/MNH/write_lfifm1_for_diag.f90 b/src/MNH/write_lfifm1_for_diag.f90
index 824d2834d..58f88436d 100644
--- a/src/MNH/write_lfifm1_for_diag.f90
+++ b/src/MNH/write_lfifm1_for_diag.f90
@@ -206,11 +206,13 @@ USE MODI_RADAR_SIMULATOR
USE MODD_DUST
USE MODD_CSTS_DUST
USE MODD_SALT
+USE MODD_BLOWSNOW
USE MODD_CH_AEROSOL
USE MODD_CH_AERO_n
USE MODD_CH_MNHC_n
USE MODE_DUST_PSD
USE MODE_SALT_PSD
+USE MODE_BLOWSNOW_PSD
USE MODE_AERO_PSD
USE MODI_GRADIENT_M
USE MODI_GRADIENT_W
@@ -224,6 +226,7 @@ USE MODI_UV_TO_ZONAL_AND_MERID
USE MODI_CALCSOUND
USE MODI_FREE_ATM_PROFILE
USE MODI_GPS_ZENITH
+USE MODI_CONTRAV
!
USE MODE_FM, ONLY: IO_FILE_CLOSE_ll,IO_FILE_OPEN_ll
USE MODE_GRIDPROJ
@@ -293,6 +296,8 @@ REAL,DIMENSION(SIZE(XSVT,1),SIZE(XSVT,2),SIZE(XSVT,3),NMODE_DST*2):: ZSDSTDEP
REAL,DIMENSION(SIZE(XSVT,1),SIZE(XSVT,2),SIZE(XSVT,3),NMODE_SLT*2):: ZSSLTDEP
REAL,DIMENSION(:,:,:,:), ALLOCATABLE :: ZSIG_DST, ZRG_DST, ZN0_DST
REAL,DIMENSION(:,:,:,:), ALLOCATABLE :: ZSIG_SLT, ZRG_SLT, ZN0_SLT
+REAL,DIMENSION(:,:,:), ALLOCATABLE :: ZBET_SNW, ZRG_SNW
+REAL,DIMENSION(:,:,:,:), ALLOCATABLE :: ZMA_SNW
REAL,DIMENSION(:,:,:), ALLOCATABLE :: ZRHOT, ZTMP ! work array
!
! GBOTUP = True does clustering from bottom up to top, False top down to surface
@@ -1293,8 +1298,96 @@ IF (LCONDSAMP) THEN
CALL IO_WRITE_FIELD(TPFILE,TZFIELD,XSVT(:,:,:,JSV))
END DO
END IF
+!
+! Blowing snow variables
+!
+IF(LBLOWSNOW) THEN
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%NGRID = 1
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 3
+ TZFIELD%LTIMEDEP = .TRUE.
+ !
+ TZFIELD%CUNITS = 'KG/M3/S'
+ TZFIELD%CMNHNAME = 'SNWSUBL3D'
+ TZFIELD%CLONGNAME = TRIM(TZFIELD%CMNHNAME)
+ TZFIELD%CCOMMENT = 'X_Y_INstantaneous 3D Drifting snow sublimation flux (KG/M3/S)'
+ CALL IO_WRITE_FIELD(TPFILE,TZFIELD,XSNWSUBL3D(:,:,:))
+!
+ ZWORK21(:,:) = 0.
+ DO JK = IKB,IKE
+ ZWORK21(:,:) = ZWORK21(:,:)+XSNWSUBL3D(:,:,JK) * &
+ (XZZ(:,:,JK+1)-XZZ(:,:,JK))/XRHOLW*3600*24
+ END DO
+ ZWORK21(:,:) = ZWORK21(:,:)*1000. ! vapor water in mm unit
+ TZFIELD%CUNITS = 'mmSWE/day'
+ TZFIELD%CMNHNAME = 'COL_SNWSUBL'
+ TZFIELD%CLONGNAME = TRIM(TZFIELD%CMNHNAME)
+ TZFIELD%NDIMS = 2
+ TZFIELD%CCOMMENT = 'X_Y_Column Sublimation Rate (mmSWE/day)'
+ CALL IO_WRITE_FIELD(TPFILE,TZFIELD,ZWORK21(:,:))
+!
+ IF(.NOT.ALLOCATED(ZBET_SNW)) &
+ ALLOCATE(ZBET_SNW(SIZE(XSVT,1), SIZE(XSVT,2), SIZE(XSVT,3)))
+ IF(.NOT.ALLOCATED(ZRG_SNW)) &
+ ALLOCATE(ZRG_SNW(SIZE(XSVT,1), SIZE(XSVT,2), SIZE(XSVT,3)))
+ IF(.NOT.ALLOCATED(ZMA_SNW)) &
+ ALLOCATE(ZMA_SNW(SIZE(XSVT,1), SIZE(XSVT,2), SIZE(XSVT,3),NBLOWSNOW3D))
+!
+ CALL PPP2SNOW(XSVT(:,:,:,NSV_SNWBEG:NSV_SNWEND),XRHODREF,&
+ PBET3D=ZBET_SNW, PRG3D=ZRG_SNW, PM3D=ZMA_SNW)
+!
+
+ TZFIELD%CUNITS = 'm'
+ TZFIELD%CMNHNAME = 'SNWRGA'
+ TZFIELD%CLONGNAME = TRIM(TZFIELD%CMNHNAME)
+ TZFIELD%CCOMMENT = 'RG (mean) SNOW (m)'
+ TZFIELD%NDIMS = 3
+ CALL IO_WRITE_FIELD(TPFILE,TZFIELD,ZRG_SNW(:,:,:))
+
+
+ TZFIELD%CUNITS = 'm'
+ TZFIELD%CMNHNAME = 'SNWBETA'
+ TZFIELD%CLONGNAME = TRIM(TZFIELD%CMNHNAME)
+ TZFIELD%CCOMMENT = 'BETA SNOW (m)'
+ TZFIELD%NDIMS = 3
+ CALL IO_WRITE_FIELD(TPFILE,TZFIELD,ZBET_SNW(:,:,:))
+
+ TZFIELD%CUNITS = '#/m3'
+ TZFIELD%CMNHNAME = 'SNWNOA'
+ TZFIELD%CLONGNAME = TRIM(TZFIELD%CMNHNAME)
+ TZFIELD%CCOMMENT = 'NUM CONC SNOW (#/m3)'
+ TZFIELD%NDIMS = 3
+ CALL IO_WRITE_FIELD(TPFILE,TZFIELD,ZMA_SNW(:,:,:,1))
+
+ TZFIELD%CUNITS = 'kg/m3'
+ TZFIELD%CMNHNAME = 'SNWMASS'
+ TZFIELD%CLONGNAME = TRIM(TZFIELD%CMNHNAME)
+ TZFIELD%CCOMMENT = 'MASS CONC SNOW (kg/m3)'
+ TZFIELD%NDIMS = 3
+ CALL IO_WRITE_FIELD(TPFILE,TZFIELD,ZMA_SNW(:,:,:,2))
+
+ ZWORK21(:,:) = 0.
+ DO JK = IKB,IKE
+ ZWORK21(:,:) = ZWORK21(:,:)+ZMA_SNW(:,:,JK,2) * &
+ (XZZ(:,:,JK+1)-XZZ(:,:,JK))/XRHOLW
+ END DO
+ ZWORK21(:,:) = ZWORK21(:,:)*1000. ! vapor water in mm unit
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%NGRID = 1
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 2
+ TZFIELD%LTIMEDEP = .TRUE.
+ TZFIELD%CUNITS = 'MM SWE'
+ TZFIELD%CMNHNAME = 'THDS'
+ TZFIELD%CLONGNAME = TRIM(TZFIELD%CMNHNAME)
+ TZFIELD%CCOMMENT = 'X_Y_THickness of Drifting Snow (MM SWE)'
+ CALL IO_WRITE_FIELD(TPFILE,TZFIELD,ZWORK21(:,:))
+END IF
! Lagrangian variables
IF (LTRAJ) THEN
TZFIELD%CSTDNAME = ''
@@ -2826,6 +2919,28 @@ END IF
!
!-------------------------------------------------------------------------------
!
+!* Contravariant wind field
+!
+!
+IF(LWIND_CONTRAV) THEN!$
+ CALL CONTRAV ((/"TEST","TEST"/),(/"TEST","TEST"/),XUT,XVT,XWT,XDXX,XDYY,XDZZ,XDZX,XDZY, &
+ ZWORK31,ZWORK32,ZWORK33,2)
+
+ TZFIELD%CMNHNAME = 'WNORM'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'WNORM'
+ TZFIELD%CUNITS = 'm/s'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'X_Y_Z_W surface normal wind (m/s)'
+ TZFIELD%NGRID = 2
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 3
+ TZFIELD%LTIMEDEP = .TRUE.
+ CALL IO_WRITE_FIELD(TPFILE,TZFIELD,ZWORK33)
+
+END IF
+!-------------------------------------------------------------------------------
+!
!* Mean Sea Level Pressure in hPa
!
IF (LMSLP) THEN
diff --git a/src/MNH/write_lfin.f90 b/src/MNH/write_lfin.f90
index 8cc4fabc5..3918143c9 100644
--- a/src/MNH/write_lfin.f90
+++ b/src/MNH/write_lfin.f90
@@ -167,6 +167,7 @@ END MODULE MODI_WRITE_LFIFM_n
!! JP Chaboureau 27/11/2017 add wind tendency forcing
!! 02/2018 Q.Libois move Diagnostic related to the radiations in radiations.f90
!! Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
+!! V. Vionnet 07/2017, add blowing snow variables
!!
!-------------------------------------------------------------------------------
!
@@ -224,6 +225,8 @@ USE MODD_FOREFIRE
#endif
USE MODD_CONDSAMP
USE MODD_CH_AEROSOL
+USE MODD_BLOWSNOW
+USE MODD_BLOWSNOW_n
USE MODD_PAST_FIELD_n
USE MODD_ADV_n, ONLY: CUVW_ADV_SCHEME,XRTKEMS,CTEMP_SCHEME
USE MODD_ELEC_FLASH
@@ -948,6 +951,37 @@ IF (NSV >=1) THEN
END DO
END IF
#endif
+! Blowing snow variables
+ IF (LBLOWSNOW) THEN
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CUNITS = 'kg kg-1'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%NGRID = 1
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 3
+ TZFIELD%LTIMEDEP = .TRUE.
+ DO JSV = NSV_SNWBEG,NSV_SNWEND
+ TZFIELD%CMNHNAME=TRIM(CSNOWNAMES(JSV-NSV_SNWBEG+1))//'T'
+ TZFIELD%CLONGNAME = TRIM(TZFIELD%CMNHNAME)
+ WRITE(TZFIELD%CCOMMENT,'(A6,A3,I3.3,A8)')'X_Y_Z_','SVT',JSV,' (KG/KG)'
+ CALL IO_WRITE_FIELD(TPFILE,TZFIELD,XSVT(:,:,:,JSV))
+ JSA=JSA+1
+ END DO
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CUNITS = 'kg kg-1'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%NGRID = 1
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 2
+ TZFIELD%LTIMEDEP = .TRUE.
+ DO JSV = 1,(NSV_SNW)
+ WRITE(TZFIELD%CMNHNAME,'(A3,I3.3)')'SNOWCANO_M',JSV
+ TZFIELD%CLONGNAME = TRIM(TZFIELD%CMNHNAME)
+ WRITE(TZFIELD%CCOMMENT,'(A6,A8,I3.3,A8)')'X_Y_Z_','SNOWCANO',JSV,' (KG/KG)'
+ CALL IO_WRITE_FIELD(TPFILE,TZFIELD,XSNWCANO(:,:,JSV))
+ JSA=JSA+1
+ END DO
+ ENDIF
! Conditional sampling variables
IF (LCONDSAMP) THEN
TZFIELD%CSTDNAME = ''
@@ -1728,7 +1762,43 @@ IF (CPROGRAM /= 'IDEAL') THEN
END IF
!
END IF
-!
+
+ IF(LBLOWSNOW) THEN
+ IF (ASSOCIATED(XSNWSUBL3D)) THEN
+ IF (SIZE(XSNWSUBL3D) /= 0 ) THEN
+
+ TZFIELD%CMNHNAME = 'SNWSUBL3D'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'SNWSUBL3D'
+ TZFIELD%CUNITS = 'KG/M3/S'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'X_Y_INstantaneous 3D Drifting snow sublimation flux (KG/M3/S)'
+ TZFIELD%NGRID = 1
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 3
+ TZFIELD%LTIMEDEP = .TRUE.
+ CALL IO_WRITE_FIELD(TPFILE,TZFIELD,XSNWSUBL3D(:,:,:))
+ ZWORK2D(:,:) = 0.
+ DO JK = IKB,IKE
+ ZWORK2D(:,:) = ZWORK2D(:,:)+XSNWSUBL3D(:,:,JK) * &
+ (XZZ(:,:,JK+1)-XZZ(:,:,JK))/XRHOLW*3600*24
+ END DO
+ ZWORK2D(:,:) = ZWORK2D(:,:)*1000. ! vapor water in mm unit
+
+ TZFIELD%CMNHNAME = 'COL_SNWSUBL'
+ TZFIELD%CSTDNAME = ''
+ TZFIELD%CLONGNAME = 'COL_SNWSUBL'
+ TZFIELD%CUNITS = 'mmSWE/day'
+ TZFIELD%CDIR = 'XY'
+ TZFIELD%CCOMMENT = 'X_Y_Column Sublimation Rate (mmSWE/day)'
+ TZFIELD%NGRID = 1
+ TZFIELD%NTYPE = TYPEREAL
+ TZFIELD%NDIMS = 2
+ TZFIELD%LTIMEDEP = .TRUE.
+ CALL IO_WRITE_FIELD(TPFILE,TZFIELD,ZWORK2D(:,:))
+ END IF
+ END IF
+ ENDIF
!
!* 1.11 Forcing variables
!
diff --git a/src/SURFEX/canopy_grid_update.F90 b/src/SURFEX/canopy_grid_update.F90
index 0ad827233..33f5127cb 100644
--- a/src/SURFEX/canopy_grid_update.F90
+++ b/src/SURFEX/canopy_grid_update.F90
@@ -3,7 +3,7 @@
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!SFX_LIC for details. version 1.
! #########################################
- SUBROUTINE CANOPY_GRID_UPDATE(KI,PH,PZFORC,SB)
+ SUBROUTINE CANOPY_GRID_UPDATE(KI,PH,PZFORC,SB,OLOG_GRID)
! #########################################
!
!!**** *CANOPY_GRID_UPDATE* - set the upper levels at and just below forcing level
@@ -58,6 +58,8 @@ REAL, DIMENSION(KI), INTENT(IN) :: PZFORC ! height of wind forcing
!
TYPE(CANOPY_t), INTENT(INOUT) :: SB
!
+LOGICAL, OPTIONAL, INTENT(IN) :: OLOG_GRID ! true if logarithmic grid is used (blowing snow scheme)
+!
!* 0.2 Declarations of local variables
! -------------------------------
!
@@ -70,6 +72,7 @@ INTEGER :: JI ! loop counter on points
REAL :: ZZTOP ! altitude of top of the grid of the initial level
! ! just below forcing height
REAL :: ZDZ ! difference of height between new levels
+LOGICAL :: GLOG_GRID
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!-------------------------------------------------------------------------------
@@ -79,6 +82,12 @@ IF (LHOOK) CALL DR_HOOK('CANOPY_GRID_UPDATE',0,ZHOOK_HANDLE)
IF(ALL(SB%XZ(:,SB%NLVL)==PZFORC(:)) .AND. LHOOK) CALL DR_HOOK('CANOPY_GRID_UPDATE',1,ZHOOK_HANDLE)
IF(ALL(SB%XZ(:,SB%NLVL)==PZFORC(:))) RETURN
!
+IF (PRESENT(OLOG_GRID)) THEN
+ GLOG_GRID = OLOG_GRID
+ELSE
+ GLOG_GRID = .FALSE.
+END IF
+!
!-------------------------------------------------------------------------------
!
!* 1. set upper level to forcing height
@@ -91,24 +100,38 @@ SB%XZ(:,SB%NLVL) = PZFORC(:)
!
! determination of levels below forcing height, low enough
!
-ILEVEL=0
-DO JI=1,KI
- DO JLAYER=1,SB%NLVL-1
- IF( PZFORC(JI) > SB%XZF(JI,JLAYER+1) + 0.25 * SB%XDZ(JI,JLAYER) .AND. &
- SB%XZ(JI,JLAYER) < PZFORC(JI) ) ILEVEL(JI,JLAYER) = JLAYER
- ENDDO
- ! determination of latest level from the ones selected before
- IL(JI)=MAXVAL(ILEVEL(JI,1:SB%NLVL-1))
- !
- ICOUNT = SB%NLVL-IL(JI)-1
- !
- !* determination grid top of this level
- ZZTOP = SB%XZF(JI,IL(JI)+1) ! ZZTOP=0 for IL=0
- ZDZ = 2. * ( SB%XZ(JI,SB%NLVL)-ZZTOP ) / ( 2*ICOUNT+1 )
- DO JLAYER=1,ICOUNT
- SB%XZ(JI,JLAYER+IL(JI)) = ZZTOP + (JLAYER-0.5) * ZDZ
+IF(GLOG_GRID) THEN
+ ! Logarithmic grid when blowing snow scheme is used
+ ! Lower height for flux is taken at 30 cm above the snow surface
+ SB%XZF(:,2) = 0.3
+ SB%XZF(:,SB%NLVL) = 0.95*PZFORC(:)
+ DO JI=1,KI
+ ZDZ = (LOG(SB%XZF(JI,SB%NLVL))-LOG(SB%XZF(JI,2)))/(SB%NLVL-2)
+ DO JLAYER=2,SB%NLVL
+ SB%XZF(JI,JLAYER) = EXP((JLAYER-2)*ZDZ+LOG(SB%XZF(JI,2)))
+ SB%XZ(JI,JLAYER-1) = (SB%XZF(JI,JLAYER)+SB%XZF(JI,JLAYER-1))/2
+ END DO
END DO
-END DO
+ELSE
+ ILEVEL=0
+ DO JI=1,KI
+ DO JLAYER=1,SB%NLVL-1
+ IF( PZFORC(JI) > SB%XZF(JI,JLAYER+1) + 0.25 * SB%XDZ(JI,JLAYER) .AND. &
+ SB%XZ(JI,JLAYER) < PZFORC(JI) ) ILEVEL(JI,JLAYER) = JLAYER
+ ENDDO
+ ! determination of latest level from the ones selected before
+ IL(JI)=MAXVAL(ILEVEL(JI,1:SB%NLVL-1))
+ !
+ ICOUNT = SB%NLVL-IL(JI)-1
+ !
+ !* determination grid top of this level
+ ZZTOP = SB%XZF(JI,IL(JI)+1) ! ZZTOP=0 for IL=0
+ ZDZ = 2. * ( SB%XZ(JI,SB%NLVL)-ZZTOP ) / ( 2*ICOUNT+1 )
+ DO JLAYER=1,ICOUNT
+ SB%XZ(JI,JLAYER+IL(JI)) = ZZTOP + (JLAYER-0.5) * ZDZ
+ END DO
+ END DO
+ENDIF
!
!* 3. New grid characteristics
! ------------------------
diff --git a/src/SURFEX/compute_isba_parameters.F90 b/src/SURFEX/compute_isba_parameters.F90
index d17f22f1a..16e2a1949 100644
--- a/src/SURFEX/compute_isba_parameters.F90
+++ b/src/SURFEX/compute_isba_parameters.F90
@@ -6,7 +6,7 @@
SUBROUTINE COMPUTE_ISBA_PARAMETERS (DTCO, OREAD_BUDGETC, UG, U, &
IO, DTI, SB, S, IG, K, NK, NIG, NP, NPE, &
NAG, NISS, ISS, NCHI, CHI, ID, GB, NGB, &
- NDST, SLT, SV, HPROGRAM,HINIT,OLAND_USE, &
+ NDST, SLT, BLOWSNW, SV, HPROGRAM,HINIT,OLAND_USE, &
KI,KSV,KSW,HSV,PCO2,PRHOA, &
PZENITH,PSW_BANDS,PDIR_ALB,PSCA_ALB, &
PEMIS,PTSRAD,PTSURF, HTEST )
@@ -83,6 +83,7 @@ USE MODD_SURF_ATM_n, ONLY : SURF_ATM_t
USE MODD_DST_n, ONLY : DST_NP_t, DST_t
USE MODD_SLT_n, ONLY : SLT_t
USE MODD_SV_n, ONLY : SV_t
+USE MODD_BLOWSNW_n, ONLY : BLOWSNW_t
!
USE MODD_SFX_OASIS, ONLY : LCPL_LAND, LCPL_FLOOD, LCPL_GW, LCPL_CALVING
!
@@ -107,6 +108,7 @@ USE MODD_TOPD_PAR, ONLY : NUNIT
USE MODD_TOPODYN, ONLY : NNCAT, NMESHT
!
USE MODE_RANDOM
+USE MODE_BLOWSNW_SEDIM_LKT1D
!
USE MODI_GET_1D_MASK
USE MODI_GET_Z0REL
@@ -184,6 +186,7 @@ TYPE(GR_BIOG_NP_t), INTENT(INOUT) :: NGB
TYPE(DST_NP_t), INTENT(INOUT) :: NDST
TYPE(SLT_t), INTENT(INOUT) :: SLT
TYPE(SV_t), INTENT(INOUT) :: SV
+TYPE(BLOWSNW_t), INTENT(INOUT) :: BLOWSNW
!
CHARACTER(LEN=6), INTENT(IN) :: HPROGRAM ! program calling surf. schemes
CHARACTER(LEN=3), INTENT(IN) :: HINIT ! choice of fields to initialize
@@ -488,7 +491,8 @@ ALLOCATE(ISS%XZ0REL(KI))
! contains explicitely modules from ISBAn. It should be cleaned in a future
! version.
CALL INIT_CHEMICAL_n(ILUOUT, KSV, HSV, CHI%SVI, CHI%CCH_NAMES, CHI%CAER_NAMES, &
- HDSTNAMES=CHI%CDSTNAMES, HSLTNAMES=CHI%CSLTNAMES )
+ HDSTNAMES=CHI%CDSTNAMES, HSLTNAMES=CHI%CSLTNAMES, &
+ HSNWNAMES=CHI%CSNWNAMES )
!
IF (KSV /= 0) THEN
!
@@ -524,6 +528,32 @@ IF (KSV /= 0) THEN
CALL INIT_SLT(SLT, HPROGRAM)
END IF
!
+ IF (CHI%SVI%NSNWEQ >=1) THEN
+ ALLOCATE (BLOWSNW%XSNW_FSED(KI,CHI%SVI%NSNWEQ+1)) !Output array
+ ALLOCATE (BLOWSNW%XSNW_FTURB(KI,CHI%SVI%NSNWEQ+1)) !Output array
+ ALLOCATE (BLOWSNW%XSNW_FNET(KI,CHI%SVI%NSNWEQ+1)) !Output array
+ ALLOCATE (BLOWSNW%XSNW_FSALT(KI,CHI%SVI%NSNWEQ+1)) !Output array
+ ALLOCATE (BLOWSNW%XSFSNW(KI,CHI%SVI%NSNWEQ+1)) !Output array
+ ALLOCATE (BLOWSNW%XSNW_SUBL(KI,CHI%SVI%NSNWEQ+1)) !Output array
+ BLOWSNW%XSNW_FSED (:,:) = 0.
+ BLOWSNW%XSNW_FTURB(:,:) = 0.
+ BLOWSNW%XSNW_FNET (:,:) = 0.
+ BLOWSNW%XSNW_FSALT(:,:) = 0.
+ BLOWSNW%XSNW_SUBL (:,:) = 0.
+ BLOWSNW%XSFSNW (:,:) = 0.
+ !Read in look up tables of snow particles properties
+ !No arguments, all look up tables are defined in module
+ !mode_snowdrift_sedim_lkt
+ CALL BLOWSNW_SEDIM_LKT1D_SET
+ ELSE
+ ALLOCATE(BLOWSNW%XSNW_FSED(0,0))
+ ALLOCATE(BLOWSNW%XSNW_FTURB(0,0))
+ ALLOCATE(BLOWSNW%XSNW_FSALT(0,0))
+ ALLOCATE(BLOWSNW%XSNW_FNET(0,0))
+ ALLOCATE(BLOWSNW%XSNW_SUBL(0,0))
+ ALLOCATE(BLOWSNW%XSFSNW(0,0))
+ END IF
+
ENDIF
!
!-------------------------------------------------------------------------------
@@ -963,7 +993,7 @@ END IF
!* 12. Canopy air fields:
! -----------------
!
- CALL READ_SBL_n(DTCO, U, SB, IO%LCANOPY, HPROGRAM, "NATURE")
+ CALL READ_SBL_n(DTCO, U, SB, IO%LCANOPY,HPROGRAM, "NATURE", SV=CHI%SVI,BLOWSNW=BLOWSNW)
!
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
diff --git a/src/SURFEX/coupling_isba_canopyn.F90 b/src/SURFEX/coupling_isba_canopyn.F90
index a483689c8..737cee787 100644
--- a/src/SURFEX/coupling_isba_canopyn.F90
+++ b/src/SURFEX/coupling_isba_canopyn.F90
@@ -5,7 +5,7 @@
! ###############################################################################
SUBROUTINE COUPLING_ISBA_CANOPY_n (DTCO, UG, U, USS, SB, NAG, CHI, NCHI, DTV, ID, NGB, GB, &
ISS, NISS, IG, NIG, IO, S, K, NK, NP, NPE, NDST, SLT, &
- HPROGRAM, HCOUPLING, PTSTEP, &
+ BLOWSNW, HPROGRAM, HCOUPLING, PTSTEP, &
KYEAR, KMONTH, KDAY, PTIME, KI, KSV, KSW, PTSUN, &
PZENITH, PZENITH2, PAZIM, PZREF, PUREF, PZS, PU, PV, &
PQA, PTA, PRHOA, PSV, PCO2, HSV, PRAIN, PSNOW, PLW, &
@@ -61,6 +61,7 @@ USE MODD_SURF_ATM_n, ONLY : SURF_ATM_t
USE MODD_SSO_n, ONLY : SSO_t
USE MODD_DATA_ISBA_n, ONLY : DATA_ISBA_t
USE MODD_SLT_n, ONLY : SLT_t
+USE MODD_BLOWSNW_n, ONLY : BLOWSNW_t
!
!
USE MODD_CSTS, ONLY : XCPD
@@ -68,11 +69,14 @@ USE MODD_SURF_PAR, ONLY : XUNDEF
USE MODD_CANOPY_TURB, ONLY : XALPSBL
!
USE MODE_COUPLING_CANOPY
+USE MODE_BLOWSNW_CANOPY
!
USE MODI_INIT_ISBA_SBL
!
USE MODI_CANOPY_EVOL
USE MODI_CANOPY_GRID_UPDATE
+USE MODI_BLOWSNW_DEP
+USE MODI_CANOPY_BLOWSNW
!
USE MODI_COUPLING_ISBA_n
!
@@ -113,6 +117,7 @@ TYPE(SURF_ATM_GRID_t), INTENT(INOUT) :: UG
TYPE(SURF_ATM_t), INTENT(INOUT) :: U
TYPE(SSO_t), INTENT(INOUT) :: USS
TYPE(SLT_t), INTENT(INOUT) :: SLT
+TYPE(BLOWSNW_t), INTENT(INOUT) :: BLOWSNW
!
CHARACTER(LEN=6), INTENT(IN) :: HPROGRAM ! program calling surf. schemes
CHARACTER(LEN=1), INTENT(IN) :: HCOUPLING ! type of coupling
@@ -198,10 +203,17 @@ REAL, DIMENSION(KI) :: ZUREF ! wind forcing level
REAL, DIMENSION(KI) :: ZU ! zonal wind (m/s)
REAL, DIMENSION(KI) :: ZV ! meridian wind (m/s)
REAL, DIMENSION(KI) :: ZQA ! specific humidity (kg/m3)
+REAL, DIMENSION(KI,CHI%SVI%NSNWEQ) :: ZBLOWSNWA ! lowest atmospheric blowing snow variable
+! 1: #/kg(air) 2: kg/kg(air)
REAL, DIMENSION(KI) :: ZPEQ_A_COEF ! specific humidity implicit
REAL, DIMENSION(KI) :: ZPEQ_B_COEF ! coefficients (hum. in kg/kg)
!
!
+REAL, DIMENSION(KI,KSV) :: ZSV ! scalar variables
+! ! chemistry: first char. in HSV: '#' (molecule/m3)!
+REAL, DIMENSION(KI,KSV) :: ZSFTS ! flux of scalar var. (kg/m2/s)
+!
+!
! canopy turbulence scheme
!
REAL, DIMENSION(KI) :: ZCANOPY ! height of canopy (m)
@@ -244,7 +256,9 @@ REAL, DIMENSION(KI) ::ZCANOPY_DENSITY
REAL, DIMENSION(KI) ::ZUW_GROUND
REAL, DIMENSION(KI) ::ZDUWDU_GROUND
!
-INTEGER :: JJ, JLAYER, IMASK, JP, JI
+LOGICAL :: LLOG_GRID
+!
+INTEGER :: JJ, JLAYER, IMASK, JP, JI ,JSV
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!-------------------------------------------------------------------------------------
@@ -255,6 +269,11 @@ REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
IF (LHOOK) CALL DR_HOOK('COUPLING_ISBA_CANOPY_N',0,ZHOOK_HANDLE)
!
+!* store scalar variables$
+ZSV(:,:) = PSV(:,:)
+ZSFTS(:,:) = 0.
+!
+!
IF (IO%LCANOPY) THEN
!
!* 1.1 Updates canopy vertical grid as a function of forcing height
@@ -263,7 +282,14 @@ IF (IO%LCANOPY) THEN
!* determines where is the forcing level and modifies the upper levels of the canopy grid
!
ZCANOPY = 0.
- CALL CANOPY_GRID_UPDATE(KI,ZCANOPY,PUREF,SB)
+!
+ IF(CHI%SVI%NSNWEQ>0) THEN
+ LLOG_GRID = .TRUE.
+ ELSE
+ LLOG_GRID = .FALSE.
+ ENDIF
+!
+ CALL CANOPY_GRID_UPDATE(KI,ZCANOPY,PUREF,SB,LLOG_GRID)
!
!
!
@@ -279,6 +305,17 @@ IF (IO%LCANOPY) THEN
PSCA_SW, PSW_BANDS, PRAIN, PSNOW, PZREF, PUREF, ISS%XSSO_SLOPE )
ENDIF
!
+!
+!
+! 1.2.2 Blowing snow scheme: initialize and udpate blowing snow variables
+! before calling Canopy
+!
+!
+IF(CHI%SVI%NSNWEQ>0) THEN
+ CALL INIT_BLOWSNW_SBL(KI,SB%NLVL,CHI%SVI%NSV_SNWBEG,CHI%SVI%NSV_SNWEND,CHI%SVI%NSNWEQ, &
+ CHI%SVI%N2D_SNWBEG,CHI%SVI%N2D_SNWEND,CHI%SVI%N2DSNWEQ,PTSTEP, &
+ BLOWSNW,SB%XDZ,SB%XU,PUREF,PRHOA,ZSV,ZBLOWSNWA,SB%XBLOWSNW)
+ENDIF
!* 1.3 Allocations
! -----------
!
@@ -371,7 +408,15 @@ ELSE
ZPEW_B_COEF, ZPET_A_COEF, &
ZPET_B_COEF, ZPEQ_A_COEF, &
ZPEQ_B_COEF )
-!
+!
+IF(CHI%SVI%NSNWEQ>0) THEN
+! Compute sedimentation flux if blowing snow scheme is used without Canopy
+! See example in Vionnet et al (TC, 2014)
+! Use value of M0 and M3 at 1st atmospheric level to compute sedimentation fluxes.
+ CALL BLOWSNW_DEP( CHI%SVI%N2D_SNWBEG,CHI%SVI%N2D_SNWEND,PTA,PPA,PRHOA,ZSV,BLOWSNW%XSFSNW)
+!
+END IF
+!
END IF
!
!-------------------------------------------------------------------------------------
@@ -382,9 +427,9 @@ END IF
CALL COUPLING_ISBA_n(DTCO, UG, U, USS, NAG, CHI, NCHI, DTV, ID, NGB, GB, ISS,NISS, IG, &
NIG, IO, S, K, NK, NP, NPE, NDST, SLT, HPROGRAM, GCOUPLING, &
PTSTEP, KYEAR, KMONTH, KDAY, PTIME, KI, KSV, KSW, PTSUN, PZENITH, &
- PZENITH2, ZZREF, ZUREF, PZS, ZU, ZV, ZQA, ZTA, PRHOA, PSV, PCO2, &
+ PZENITH2, ZZREF, ZUREF, PZS, ZU, ZV, ZQA, ZTA, PRHOA, ZSV, PCO2, &
HSV, PRAIN, PSNOW, PLW, PDIR_SW, PSCA_SW, PSW_BANDS, PPS, ZPA, &
- PSFTQ, PSFTH, PSFTS, PSFCO2, PSFU, PSFV, PTRAD, PDIR_ALB, &
+ PSFTQ, PSFTH, ZSFTS, PSFCO2, PSFU, PSFV, PTRAD, PDIR_ALB, &
PSCA_ALB, PEMIS, PTSURF, PZ0, PZ0H, PQSURF, ZPEW_A_COEF, &
ZPEW_B_COEF, ZPET_A_COEF, ZPEQ_A_COEF, ZPET_B_COEF, ZPEQ_B_COEF, &
'OK' )
@@ -395,7 +440,17 @@ END IF
! ------------------------
!
IF (.NOT. IO%LCANOPY .AND. LHOOK) CALL DR_HOOK('COUPLING_ISBA_CANOPY_N',1,ZHOOK_HANDLE)
-IF (.NOT. IO%LCANOPY) RETURN
+IF (.NOT. IO%LCANOPY) THEN
+ IF(CHI%SVI%NSNWEQ>0) THEN
+ ! Compute net fluxes blowing snow fluxes sent to MNH : Turb-Sed
+ DO JSV=CHI%SVI%NSV_SNWBEG,CHI%SVI%NSV_SNWEND
+ PSFTS(:,JSV) = ZSFTS(:,JSV)-BLOWSNW%XSFSNW(:,JSV)
+ END DO
+ ELSE
+ PSFTS(:,:) = ZSFTS(:,:)
+ ENDIF
+ RETURN
+ENDIF
!
!-------------------------------------------------------------------------------------
!
@@ -430,6 +485,20 @@ IF (IO%LCANOPY_DRAG) THEN
!
END IF
!
+!* 5.bis Evolution of blowing snow variables in Canopy
+! ------------------------------------
+!
+IF(CHI%SVI%NSNWEQ>0) THEN
+
+ CALL CANOPY_BLOWSNW(SB,KI,CHI%SVI%NSNWEQ,PTSTEP,BLOWSNW,SB%XZ,PRHOA,ZBLOWSNWA, &
+ ZSFTS(:,CHI%SVI%NSV_SNWBEG:CHI%SVI%NSV_SNWEND), &
+ PSFTH,PSFTQ)
+
+ ZSFLUX_T(:) = PSFTH(:) / XCPD * ZEXNA(:) / PRHOA(:)
+ ZSFLUX_Q(:) = PSFTQ(:)
+!
+ENDIF
+!
!-------------------------------------------------------------------------------------
!
!* 6. Evolution of canopy air due to these impacts
@@ -470,6 +539,25 @@ END IF
! ----------------------------------------
!
IF (ID%O%N2M>=1) CALL INIT_2M_10M(SB, ID%D, PU, PV, ZWIND, PRHOA )
+
+!-------------------------------------------------------------------------------------
+!
+!* 8. Update blowing snow variable send to MNH
+! -------------------------------------
+!
+IF(CHI%SVI%NSNWEQ>0) THEN
+ CALL UPDATE_BLOWSNW_SBL(KI,SB%NLVL,CHI%SVI%NSV_SNWBEG,CHI%SVI%NSV_SNWEND,CHI%SVI%NSNWEQ, &
+ CHI%SVI%N2D_SNWBEG,CHI%SVI%N2D_SNWEND,CHI%SVI%N2DSNWEQ, &
+ BLOWSNW,SB%XDZ,SB%XU,PRHOA,PUREF,SB%XBLOWSNW,ZSFTS)
+ENDIF
+!
+!
+!-------------------------------------------------------------------------------------
+!
+!* 9. Update scalar flux
+! -------------------------------------
+!
+PSFTS(:,:) = ZSFTS(:,:)
!
IF (LHOOK) CALL DR_HOOK('COUPLING_ISBA_CANOPY_N',1,ZHOOK_HANDLE)
!
diff --git a/src/SURFEX/coupling_isba_orographyn.F90 b/src/SURFEX/coupling_isba_orographyn.F90
index 4c17a9c4a..df3e016ac 100644
--- a/src/SURFEX/coupling_isba_orographyn.F90
+++ b/src/SURFEX/coupling_isba_orographyn.F90
@@ -5,7 +5,7 @@
! ###############################################################################
SUBROUTINE COUPLING_ISBA_OROGRAPHY_n (DTCO, UG, U, USS, SB, NAG, CHI, NCHI, DTV, ID, NGB, GB, &
ISS, NISS, IG, NIG, IO, S, K, NK, NP, NPE, NDST, SLT, &
- HPROGRAM, HCOUPLING, PTSTEP, &
+ BLOWSNW,HPROGRAM, HCOUPLING, PTSTEP, &
KYEAR, KMONTH, KDAY, PTIME, KI, KSV, KSW, PTSUN, &
PZENITH, PZENITH2, PAZIM, PZREF, PUREF, PZS, PU, PV, &
PQA, PTA, PRHOA, PSV, PCO2, HSV, PRAIN, PSNOW, PLW, &
@@ -63,6 +63,7 @@ USE MODD_SURF_ATM_n, ONLY : SURF_ATM_t
USE MODD_SSO_n, ONLY : SSO_t, SSO_NP_t
USE MODD_DATA_ISBA_n, ONLY : DATA_ISBA_t
USE MODD_SLT_n, ONLY : SLT_t
+USE MODD_BLOWSNW_n, ONLY : BLOWSNW_t
!
USE MODD_SURF_PAR,ONLY : XUNDEF
USE MODD_CSTS, ONLY : XSTEFAN, XCPD, XRD, XP00
@@ -106,6 +107,7 @@ TYPE(SURF_ATM_GRID_t), INTENT(INOUT) :: UG
TYPE(SURF_ATM_t), INTENT(INOUT) :: U
TYPE(SSO_t), INTENT(INOUT) :: USS
TYPE(SLT_t), INTENT(INOUT) :: SLT
+TYPE(BLOWSNW_t), INTENT(INOUT) :: BLOWSNW
!
CHARACTER(LEN=6), INTENT(IN) :: HPROGRAM ! program calling surf. schemes
CHARACTER(LEN=1), INTENT(IN) :: HCOUPLING ! type of coupling
@@ -322,7 +324,7 @@ ENDIF
!
CALL COUPLING_ISBA_CANOPY_n(DTCO, UG, U, USS, SB, NAG, CHI, NCHI, DTV, ID, NGB, GB, &
ISS, NISS, IG, NIG, IO, S, K, NK, NP, NPE, NDST, SLT, &
- HPROGRAM, HCOUPLING, PTSTEP, &
+ BLOWSNW,HPROGRAM, HCOUPLING, PTSTEP, &
KYEAR, KMONTH, KDAY, PTIME, KI, KSV, KSW, PTSUN, PZENITH, &
PZENITH2, PAZIM, PZREF, PUREF, PZS, PU, PV, ZQA, ZTA, &
ZRHOA, PSV, PCO2, HSV, ZRAIN, ZSNOW, ZLW, ZDIR_SW, &
diff --git a/src/SURFEX/coupling_isba_svatn.F90 b/src/SURFEX/coupling_isba_svatn.F90
index 645616d2e..a01fe18e4 100644
--- a/src/SURFEX/coupling_isba_svatn.F90
+++ b/src/SURFEX/coupling_isba_svatn.F90
@@ -3,7 +3,7 @@
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!SFX_LIC for details. version 1.
! ###############################################################################
-SUBROUTINE COUPLING_ISBA_SVAT_n (DTCO, UG, U, USS, IM, NDST, SLT, HPROGRAM, HCOUPLING, PTSTEP, &
+SUBROUTINE COUPLING_ISBA_SVAT_n (DTCO, UG, U, USS, IM, NDST, SLT, BLOWSNW, HPROGRAM, HCOUPLING, PTSTEP, &
KYEAR, KMONTH, KDAY, PTIME, KI, KSV, KSW, PTSUN, PZENITH, &
PZENITH2, PAZIM, PZREF, PUREF, PZS, PU, PV, PQA, PTA, PRHOA, &
PSV, PCO2, HSV, PRAIN, PSNOW, PLW, PDIR_SW, PSCA_SW, &
@@ -48,6 +48,7 @@ USE MODD_SURF_ATM_n, ONLY : SURF_ATM_t
USE MODD_SSO_n, ONLY : SSO_t
USE MODD_DST_n, ONLY : DST_NP_t
USE MODD_SLT_n, ONLY : SLT_t
+USE MODD_BLOWSNW_n, ONLY : BLOWSNW_t
!
USE MODD_SURF_PAR, ONLY : XUNDEF
!
@@ -67,6 +68,7 @@ TYPE(SURF_ATM_t), INTENT(INOUT) :: U
TYPE(SSO_t), INTENT(INOUT) :: USS
TYPE(DST_NP_t), INTENT(INOUT) :: NDST
TYPE(SLT_t), INTENT(INOUT) :: SLT
+TYPE(BLOWSNW_t), INTENT(INOUT) :: BLOWSNW
!
CHARACTER(LEN=6), INTENT(IN) :: HPROGRAM ! program calling surf. schemes
CHARACTER(LEN=1), INTENT(IN) :: HCOUPLING ! type of coupling
@@ -232,7 +234,7 @@ DO JT=1,IT
!
CALL COUPLING_ISBA_OROGRAPHY_n(DTCO, UG, U, USS, IM%SB, IM%NAG, IM%CHI, IM%NCHI, IM%DTV, &
IM%ID, IM%NGB, IM%GB, IM%ISS, IM%NISS, IM%G, IM%NG, IM%O, &
- IM%S, IM%K, IM%NK, IM%NP, IM%NPE, NDST, SLT, &
+ IM%S, IM%K, IM%NK, IM%NP, IM%NPE, NDST, SLT, BLOWSNW, &
HPROGRAM, HCOUPLING, ZTSTEP, KYEAR, KMONTH, KDAY, PTIME, KI, &
KSV, KSW, PTSUN, PZENITH, PZENITH2, PAZIM, PZREF, PUREF, PZS, &
PU, PV, PQA, PTA, PRHOA, PSV, PCO2, HSV, PRAIN, PSNOW, PLW, &
diff --git a/src/SURFEX/coupling_isban.F90 b/src/SURFEX/coupling_isban.F90
index f2475b18f..7f4ede4b7 100644
--- a/src/SURFEX/coupling_isban.F90
+++ b/src/SURFEX/coupling_isban.F90
@@ -306,6 +306,11 @@ REAL :: ZCONVERTFACM0_SLT, ZCONVERTFACM0_DST
REAL :: ZCONVERTFACM3_SLT, ZCONVERTFACM3_DST
REAL :: ZCONVERTFACM6_SLT, ZCONVERTFACM6_DST
!
+! for blowing snow scheme
+!
+REAL, DIMENSION(KI,CHI%SVI%N2DSNWEQ) :: ZP_BLOWSNW_FLUX ! blowing snow fluxes
+REAL, DIMENSION(KI,CHI%SVI%NSNWEQ) :: ZP_BLOWSNW_CONC ! blowing snow concentration
+!
! dimensions and loop counters
!
INTEGER :: ISWB ! number of spectral shortwave bands
@@ -860,6 +865,24 @@ IF(LNOSOF) ZP_SLOPE_COS(:) = 1.0
! (see update_frac_alb_emis_isban.f90) in order to close the energy budget
! between surfex and the atmosphere. This fact do not change the offline runs.
!
+! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+! Blowing snow scheme
+! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+IF(CHI%SVI%NSNWEQ>0.) THEN
+ ZP_BLOWSNW_CONC(:,:) = ZP_SV(:,CHI%SVI%NSV_SNWBEG:CHI%SVI%NSV_SNWEND)
+ ZP_BLOWSNW_FLUX(:,:) = ZP_SV(:,CHI%SVI%N2D_SNWBEG:CHI%SVI%N2D_SNWEND)
+! ZP_BLOWSNW_FLUX IN : fluxes sent from Canopy:
+! [1] number sedim. flux (#/m2/s)
+! [2] mass sedim flux (kg{snow}/m2/s)
+! [3] contrib. saltation (kg{snow}/m2/s)
+! OUT : fluxes towards Canopy:
+! [1] number turbulent flux (#/m2/s)
+! [2] mass turbulent flux (kg{snow}/m2/s)
+! [3] updated streamwise saltation flux (kg{snow}/m2/s)
+ELSE
+ ZP_BLOWSNW_CONC(:,:) = XUNDEF
+ ZP_BLOWSNW_FLUX(:,:) = XUNDEF
+END IF
!
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! No implicitation of Tdeep
@@ -939,7 +962,7 @@ ZIRRIG_GR(:)= 0.
ZP_ALBVIS_TSOIL, ZPALPHAN, ZZ0G_WITHOUT_SNOW, ZZ0_MEBV, ZZ0H_MEBV, ZZ0EFF_MEBV, &
ZZ0_MEBN, ZZ0H_MEBN, ZZ0EFF_MEBN, ZP_TDEEP_A, ZP_CO2, ZP_FFGNOS, ZP_FFVNOS, &
ZP_EMIS, ZP_USTAR, ZP_AC_AGG, ZP_HU_AGG, ZP_RESP_BIOMASS_INST, ZP_DEEP_FLUX, &
- ZIRRIG_GR )
+ ZIRRIG_GR, ZP_BLOWSNW_FLUX, ZP_BLOWSNW_CONC )
!
ZP_TRAD = DK%XTSRAD
DK%XLE = PEK%XLE
@@ -1027,7 +1050,14 @@ END WHERE
ZP_SFTS(:,:) = 0.
!
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-!
+! Blowing snow scheme
+! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+IF(CHI%SVI%NSNWEQ>0) THEN
+! Store emitted turbulent flux 1: number (#/m2/s); 2: mass (kg/m2/s)
+ ZP_SFTS(:,CHI%SVI%NSV_SNWBEG:CHI%SVI%NSV_SNWEND) = ZP_BLOWSNW_FLUX(:,1:CHI%SVI%NSNWEQ)
+! Store streamwise saltation flux (kg/m/s)
+ ZP_SFTS(:,CHI%SVI%N2D_SNWEND) = ZP_BLOWSNW_FLUX(:,CHI%SVI%NSNWEQ+1)
+END IF
! --------------------------------------------------------------------------------------
! Chemical dry deposition :
! --------------------------------------------------------------------------------------
diff --git a/src/SURFEX/coupling_naturen.F90 b/src/SURFEX/coupling_naturen.F90
index b811555be..a19911e45 100644
--- a/src/SURFEX/coupling_naturen.F90
+++ b/src/SURFEX/coupling_naturen.F90
@@ -3,7 +3,7 @@
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!SFX_LIC for details. version 1.
! ###############################################################################
-SUBROUTINE COUPLING_NATURE_n (DTCO, UG, U, USS, IM, DTZ, DGO, DL, DLC, NDST, SLT, &
+SUBROUTINE COUPLING_NATURE_n (DTCO, UG, U, USS, IM, DTZ, DGO, DL, DLC, NDST, SLT, BLOWSNW, &
HPROGRAM, HCOUPLING, PTIMEC, PTSTEP, KYEAR, KMONTH, KDAY, PTIME, &
KI, KSV, KSW, PTSUN, PZENITH, PZENITH2, PAZIM, PZREF, PUREF, PZS,&
PU, PV, PQA, PTA, PRHOA, PSV, PCO2, HSV, PRAIN, PSNOW, PLW, &
@@ -46,6 +46,7 @@ USE MODD_DATA_TSZ0_n, ONLY : DATA_TSZ0_t
USE MODD_DIAG_n, ONLY : DIAG_t, DIAG_OPTIONS_t
USE MODD_DST_n, ONLY : DST_NP_t
USE MODD_SLT_n, ONLY : SLT_t
+USE MODD_BLOWSNW_n, ONLY : BLOWSNW_t
!
USE MODD_CSTS, ONLY : XTT
!
@@ -73,6 +74,7 @@ TYPE(DIAG_t), INTENT(INOUT) :: DL
TYPE(DIAG_t), INTENT(INOUT) :: DLC
TYPE(DST_NP_t), INTENT(INOUT) :: NDST
TYPE(SLT_t), INTENT(INOUT) :: SLT
+TYPE(BLOWSNW_t), INTENT(INOUT) :: BLOWSNW
!
CHARACTER(LEN=6), INTENT(IN) :: HPROGRAM ! program calling surf. schemes
CHARACTER(LEN=1), INTENT(IN) :: HCOUPLING ! type of coupling
@@ -150,15 +152,15 @@ REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
IF (LHOOK) CALL DR_HOOK('COUPLING_NATURE_N',0,ZHOOK_HANDLE)
IF (U%CNATURE=='ISBA ') THEN
- CALL COUPLING_ISBA_SVAT_n(DTCO, UG, U, USS, IM, NDST, SLT, HPROGRAM, HCOUPLING, PTSTEP, &
- KYEAR, KMONTH, KDAY, PTIME, KI, KSV, KSW, PTSUN, PZENITH, &
+ CALL COUPLING_ISBA_SVAT_n(DTCO, UG, U, USS, IM, NDST, SLT, BLOWSNW, HPROGRAM, HCOUPLING, PTSTEP, &
+ KYEAR, KMONTH, KDAY, PTIME, KI, KSV, KSW, PTSUN, PZENITH,&
PZENITH2, PAZIM, PZREF, PUREF, PZS, PU, PV, PQA, PTA, PRHOA, PSV,&
PCO2, HSV, PRAIN, PSNOW, PLW, PDIR_SW, PSCA_SW, PSW_BANDS, PPS, &
PPA, PSFTQ, PSFTH, PSFTS, PSFCO2, PSFU, PSFV, PTRAD, PDIR_ALB, &
PSCA_ALB, PEMIS, PTSURF, PZ0, PZ0H, PQSURF, PPEW_A_COEF, &
PPEW_B_COEF, PPET_A_COEF, PPEQ_A_COEF, PPET_B_COEF, PPEQ_B_COEF,'OK')
ELSE IF (U%CNATURE=='TSZ0 ') THEN
- CALL COUPLING_TSZ0_n(DTCO, UG, U, USS, IM, DTZ, NDST, SLT, &
+ CALL COUPLING_TSZ0_n(DTCO, UG, U, USS, IM, DTZ, NDST, SLT, BLOWSNW, &
HPROGRAM, HCOUPLING, PTSTEP, KYEAR, KMONTH, KDAY, PTIMEC, KI, &
KSV, KSW, PTSUN, PZENITH, PZENITH2, PAZIM, PZREF, PUREF, PZS, &
PU, PV, PQA, PTA, PRHOA, PSV, PCO2, HSV, PRAIN, PSNOW, PLW, &
diff --git a/src/SURFEX/coupling_surf_atmn.F90 b/src/SURFEX/coupling_surf_atmn.F90
index 112179e74..f9fdb2f8e 100644
--- a/src/SURFEX/coupling_surf_atmn.F90
+++ b/src/SURFEX/coupling_surf_atmn.F90
@@ -601,8 +601,9 @@ ELSEIF (KTILE==2) THEN
!
ELSEIF (KTILE==3) THEN
!
- CALL COUPLING_NATURE_n(YSC%DTCO, YSC%UG, YSC%U, YSC%USS, YSC%IM, YSC%DTZ, YSC%DLO, YSC%DL, &
- YSC%DLC, YSC%NDST, YSC%SLT, HPROGRAM, HCOUPLING, PTIMEC, PTSTEP, &
+ CALL COUPLING_NATURE_n(YSC%DTCO, YSC%UG, YSC%U, YSC%USS, YSC%IM, YSC%DTZ, YSC%DLO, YSC%DL, &
+ YSC%DLC, YSC%NDST, YSC%SLT, YSC%BLOWSNW, &
+ HPROGRAM, HCOUPLING, PTIMEC, PTSTEP, &
KYEAR, KMONTH, KDAY, PTIME, YSC%U%NSIZE_NATURE, KSV, KSW, ZP_TSUN, &
ZP_ZENITH, ZP_ZENITH2, ZP_AZIM, ZP_ZREF, ZP_UREF, ZP_ZS, ZP_U, ZP_V, &
ZP_QA, ZP_TA, ZP_RHOA, ZP_SV, ZP_CO2, HSV, ZP_RAIN, ZP_SNOW, ZP_LW, &
diff --git a/src/SURFEX/coupling_tsz0n.F90 b/src/SURFEX/coupling_tsz0n.F90
index 1596f2718..712403dbf 100644
--- a/src/SURFEX/coupling_tsz0n.F90
+++ b/src/SURFEX/coupling_tsz0n.F90
@@ -3,7 +3,7 @@
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!SFX_LIC for details. version 1.
! ###############################################################################
-SUBROUTINE COUPLING_TSZ0_n (DTCO, UG, U, USS, IM, DTZ, NDST, SLT, HPROGRAM, HCOUPLING, &
+SUBROUTINE COUPLING_TSZ0_n (DTCO, UG, U, USS, IM, DTZ, NDST, SLT,BLOWSNW, HPROGRAM, HCOUPLING, &
PTSTEP, KYEAR, KMONTH, KDAY, PTIME, KI, KSV, KSW, PTSUN,&
PZENITH, PZENITH2, PAZIM, PZREF, PUREF, PZS, PU, PV, &
PQA, PTA, PRHOA, PSV, PCO2, HSV, PRAIN, PSNOW, PLW, &
@@ -50,6 +50,7 @@ USE MODD_DATA_TSZ0_n, ONLY : DATA_TSZ0_t
USE MODD_DATA_ISBA_n, ONLY : DATA_ISBA_t
USE MODD_DST_n, ONLY : DST_NP_t
USE MODD_SLT_n, ONLY : SLT_t
+USE MODD_BLOWSNW_n, ONLY : BLOWSNW_t
!
!
USE MODD_SURF_PAR, ONLY : XUNDEF
@@ -74,6 +75,7 @@ TYPE(SSO_t), INTENT(INOUT) :: USS
TYPE(DATA_TSZ0_t), INTENT(INOUT) :: DTZ
TYPE(DST_NP_t), INTENT(INOUT) :: NDST
TYPE(SLT_t), INTENT(INOUT) :: SLT
+TYPE(BLOWSNW_t), INTENT(INOUT) :: BLOWSNW
!
CHARACTER(LEN=6), INTENT(IN) :: HPROGRAM ! program calling surf. schemes
CHARACTER(LEN=1), INTENT(IN) :: HCOUPLING ! type of coupling
@@ -207,7 +209,7 @@ ENDDO
!
CALL COUPLING_ISBA_OROGRAPHY_n(DTCO, UG, U, USS, IM%SB, IM%NAG, IM%CHI, IM%NCHI, IM%DTV, IM%ID, &
IM%NGB, IM%GB, IM%ISS, IM%NISS, IM%G, IM%NG, IM%O, IM%S, IM%K, IM%NK, &
- IM%NP, IM%NPE, NDST, SLT, HPROGRAM, 'E', 0.001, KYEAR, &
+ IM%NP, IM%NPE, NDST, SLT,BLOWSNW, HPROGRAM, 'E', 0.001, KYEAR, &
KMONTH, KDAY, PTIME, KI, KSV, KSW, PTSUN, PZENITH, &
PZENITH2, PAZIM, PZREF, PUREF, PZS, PU, PV, PQA, PTA, &
PRHOA, PSV, PCO2, HSV, PRAIN, PSNOW, PLW, PDIR_SW, &
diff --git a/src/SURFEX/init_chemicaln.F90 b/src/SURFEX/init_chemicaln.F90
index 653515325..2b082d470 100644
--- a/src/SURFEX/init_chemicaln.F90
+++ b/src/SURFEX/init_chemicaln.F90
@@ -4,7 +4,7 @@
!SFX_LIC for details. version 1.
!#############################################################
SUBROUTINE INIT_CHEMICAL_n(KLUOUT, KSV, HSV, SV, HCH_NAMES, HAER_NAMES, &
- HDSTNAMES, HSLTNAMES )
+ HDSTNAMES, HSLTNAMES,HSNWNAMES )
!#############################################################
!
!!**** *INIT_CHEMICAL_n* - routine to initialize CHEMICAL SPECIES
@@ -48,6 +48,7 @@ USE MODD_SLT_SURF, ONLY: LVARSIG_SLT, NSLTMDE, NSLT_MDEBEG, LRGFIX_SLT, J
USE MODI_CH_INIT_NAMES
USE MODI_DSLT_INIT_NAMES
USE MODI_DSLT_INIT_MODES
+USE MODI_BLOWSNW_INIT_NAMES
!
USE YOMHOOK ,ONLY : LHOOK, DR_HOOK
USE PARKIND1 ,ONLY : JPRB
@@ -66,6 +67,7 @@ INTEGER, INTENT(IN) :: KSV ! number of scalars
CHARACTER(LEN=6), DIMENSION(:), POINTER, OPTIONAL :: HDSTNAMES
CHARACTER(LEN=6), DIMENSION(:), POINTER, OPTIONAL :: HSLTNAMES
+ CHARACTER(LEN=6), DIMENSION(:), POINTER, OPTIONAL :: HSNWNAMES
!
!* 0.2 Declarations of local variables
! -------------------------------
@@ -156,10 +158,29 @@ IF (KSV /= 0) THEN
ENDIF
END IF
+ CALL BLOWSNW_INIT_NAMES( &
+ KLUOUT, & !I [idx] index of writing unit
+ HSV, & !I [char] list of scalar variables
+ SV%NSNWEQ, & !O [nbr] number of blowing snow related tracers
+ SV%NSV_SNWBEG, & !O [idx] first blowing snow related scalar variable
+ SV%NSV_SNWEND, & !O [idx] last blowing snow related scalar variable
+ SV%N2DSNWEQ, & !O [nbr] number of 2D blowing snow related tracers
+ SV%N2D_SNWBEG, & !O [idx] first 2D blowing snow related scalar variable
+ SV%N2D_SNWEND & !O [idx] last 2D blowing snow related scalar variable
+ )
+
+ IF (PRESENT(HSNWNAMES)) THEN
+ IF (SV%NSNWEQ >=1) THEN
+ IF(.NOT. ASSOCIATED(HSNWNAMES)) ALLOCATE (HSNWNAMES(SV%NSNWEQ))
+ HSNWNAMES(:) = SV%CSV(SV%NSV_SNWBEG:SV%NSV_SNWEND)
+ ENDIF
+ END IF
+
ELSE
ALLOCATE(SV%CSV (0))
IF (PRESENT(HDSTNAMES)) ALLOCATE(HDSTNAMES(0))
IF (PRESENT(HSLTNAMES)) ALLOCATE(HSLTNAMES(0))
+ IF (PRESENT(HSNWNAMES)) ALLOCATE(HSNWNAMES(0))
ENDIF
!
IF (LHOOK) CALL DR_HOOK('INIT_CHEMICAL_n',1,ZHOOK_HANDLE)
diff --git a/src/SURFEX/init_flaken.F90 b/src/SURFEX/init_flaken.F90
index 57102a984..c3b1de86a 100644
--- a/src/SURFEX/init_flaken.F90
+++ b/src/SURFEX/init_flaken.F90
@@ -309,7 +309,8 @@ PTSURF(:) = FM%F%XTS(:)
!
CALL INIT_CHEMICAL_n(ILUOUT, KSV, HSV, FM%CHF%SVF, &
FM%CHF%CCH_NAMES, FM%CHF%CAER_NAMES, &
- HDSTNAMES=FM%CHF%CDSTNAMES, HSLTNAMES=FM%CHF%CSLTNAMES )
+ HDSTNAMES=FM%CHF%CDSTNAMES, HSLTNAMES=FM%CHF%CSLTNAMES , &
+ HSNWNAMES=FM%CHF%CSNWNAMES )
!
!* depositiion scheme
!
diff --git a/src/SURFEX/init_isban.F90 b/src/SURFEX/init_isban.F90
index 0cb1e95de..cac8c7a7a 100644
--- a/src/SURFEX/init_isban.F90
+++ b/src/SURFEX/init_isban.F90
@@ -1,10 +1,11 @@
-!SFX_LIC Copyright 2004-2018 CNRS, Meteo-France and Universite Paul Sabatier
+!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!SFX_LIC for details. version 1.
!#############################################################
SUBROUTINE INIT_ISBA_n (DTCO, OREAD_BUDGETC, UG, U, USS, GCP, IM, DTZ,&
- NDST, SLT, SV, HPROGRAM, HINIT, OLAND_USE, &
+ NDST, SLT, BLOWSNW, SV, HPROGRAM, HINIT, &
+ OLAND_USE, &
KI, KSV, KSW, HSV, PCO2, PRHOA, PZENITH, &
PAZIM, PSW_BANDS, PDIR_ALB, PSCA_ALB, PEMIS, &
PTSRAD, PTSURF, KYEAR, KMONTH, KDAY, PTIME, &
@@ -71,6 +72,7 @@ USE MODD_DATA_TSZ0_n, ONLY : DATA_TSZ0_t
USE MODD_DST_n, ONLY : DST_NP_t
USE MODD_SLT_n, ONLY : SLT_t
USE MODD_SV_n, ONLY : SV_t
+USE MODD_BLOWSNW_n, ONLY : BLOWSNW_t
!
USE MODD_TYPE_DATE_SURF, ONLY : DATE
!
@@ -147,6 +149,8 @@ TYPE(ISBA_MODEL_t), INTENT(INOUT) :: IM
TYPE(DATA_TSZ0_t), INTENT(INOUT) :: DTZ
TYPE(DST_NP_t), INTENT(INOUT) :: NDST
TYPE(SLT_t), INTENT(INOUT) :: SLT
+TYPE(BLOWSNW_t), INTENT(INOUT) :: BLOWSNW
+
TYPE(SV_t), INTENT(INOUT) :: SV
!
CHARACTER(LEN=6), INTENT(IN) :: HPROGRAM ! program calling surf. schemes
@@ -468,7 +472,7 @@ CALL COMPUTE_ISBA_PARAMETERS(DTCO, OREAD_BUDGETC, UG, U, &
IM%O, IM%DTV, IM%SB, IM%S, IM%G, IM%K, IM%NK, &
IM%NG, IM%NP, IM%NPE, IM%NAG, IM%NISS, IM%ISS, &
IM%NCHI, IM%CHI, IM%ID, IM%GB, IM%NGB, &
- NDST, SLT, SV, HPROGRAM,HINIT,OLAND_USE, &
+ NDST, SLT,BLOWSNW, SV,HPROGRAM,HINIT,OLAND_USE,&
KI,KSV,KSW, HSV,ZCO2,PRHOA, &
PZENITH,PSW_BANDS,PDIR_ALB,PSCA_ALB, &
PEMIS,PTSRAD,PTSURF,HTEST )
diff --git a/src/SURFEX/init_naturen.F90 b/src/SURFEX/init_naturen.F90
index 5c0b1aa95..372e883ea 100644
--- a/src/SURFEX/init_naturen.F90
+++ b/src/SURFEX/init_naturen.F90
@@ -4,7 +4,7 @@
!SFX_LIC for details. version 1.
! #############################################################
SUBROUTINE INIT_NATURE_n (DTCO, OREAD_BUDGETC, UG, U, USS, GCP, IM, &
- DTZ, DGO, DL, DLC, NDST, SLT, SV, &
+ DTZ, DGO, DL, DLC, NDST, SLT, BLOWSNW,SV, &
HPROGRAM,HINIT,OLAND_USE, KI,KSV,KSW, &
HSV,PCO2,PRHOA,PZENITH,PAZIM,PSW_BANDS, &
PDIR_ALB,PSCA_ALB,PEMIS,PTSRAD,PTSURF, &
@@ -62,6 +62,7 @@ USE MODD_DATA_TSZ0_n, ONLY : DATA_TSZ0_t
USE MODD_DST_n, ONLY : DST_NP_t
USE MODD_SLT_n, ONLY : SLT_t
USE MODD_SV_n, ONLY : SV_t
+USE MODD_BLOWSNW_n, ONLY : BLOWSNW_t
!
USE MODD_CSTS, ONLY : XTT
!
@@ -91,6 +92,7 @@ TYPE(DIAG_t), INTENT(INOUT) :: DLC
TYPE(DST_NP_t), INTENT(INOUT) :: NDST
TYPE(SLT_t), INTENT(INOUT) :: SLT
TYPE(SV_t), INTENT(INOUT) :: SV
+TYPE(BLOWSNW_t), INTENT(INOUT) :: BLOWSNW
!
CHARACTER(LEN=6), INTENT(IN) :: HPROGRAM ! program calling surf. schemes
CHARACTER(LEN=3), INTENT(IN) :: HINIT ! choice of fields to initialize
@@ -143,7 +145,7 @@ ELSE IF (U%CNATURE=='FLUX ') THEN
PTSRAD, PTSURF, 'OK' )
ELSE IF (U%CNATURE=='ISBA ' .OR. U%CNATURE=='TSZ0') THEN
CALL INIT_ISBA_n(DTCO, OREAD_BUDGETC, UG, U, USS, GCP, &
- IM, DTZ, NDST, SLT, SV, &
+ IM, DTZ, NDST, SLT, BLOWSNW, SV, &
HPROGRAM, HINIT, OLAND_USE, KI, KSV, KSW, HSV, &
PCO2, PRHOA, PZENITH, PAZIM, PSW_BANDS, &
PDIR_ALB, PSCA_ALB, PEMIS, PTSRAD, PTSURF, &
diff --git a/src/SURFEX/init_seafluxn.F90 b/src/SURFEX/init_seafluxn.F90
index 26920994c..e951f0d5c 100644
--- a/src/SURFEX/init_seafluxn.F90
+++ b/src/SURFEX/init_seafluxn.F90
@@ -405,7 +405,8 @@ ENDIF
!
CALL INIT_CHEMICAL_n(ILUOUT, KSV, HSV, SM%CHS%SVS, &
SM%CHS%CCH_NAMES, SM%CHS%CAER_NAMES, &
- HDSTNAMES=SM%CHS%CDSTNAMES, HSLTNAMES=SM%CHS%CSLTNAMES )
+ HDSTNAMES=SM%CHS%CDSTNAMES, HSLTNAMES=SM%CHS%CSLTNAMES, &
+ HSNWNAMES=SM%CHS%CSNWNAMES )
!
!* deposition scheme
!
diff --git a/src/SURFEX/init_surf_atmn.F90 b/src/SURFEX/init_surf_atmn.F90
index b27a6b3ca..809e785cd 100644
--- a/src/SURFEX/init_surf_atmn.F90
+++ b/src/SURFEX/init_surf_atmn.F90
@@ -597,7 +597,7 @@ ZFRAC_TILE(:,JTILE) = YSC%U%XNATURE(:)
IF (YSC%U%NDIM_NATURE>0) &
CALL INIT_NATURE_n(YSC%DTCO, YSC%DUO%LREAD_BUDGETC, YSC%UG, YSC%U, &
YSC%USS, YSC%GCP, YSC%IM, YSC%DTZ, YSC%DLO, YSC%DL,&
- YSC%DLC, YSC%NDST, YSC%SLT, YSC%SV, &
+ YSC%DLC, YSC%NDST, YSC%SLT,YSC%BLOWSNW, YSC%SV, &
HPROGRAM,HINIT,OLAND_USE,YSC%U%NSIZE_NATURE, &
KSV,KSW, HSV,ZP_CO2,ZP_RHOA, &
ZP_ZENITH,ZP_AZIM,PSW_BANDS,ZP_DIR_ALB,ZP_SCA_ALB, &
diff --git a/src/SURFEX/init_watfluxn.F90 b/src/SURFEX/init_watfluxn.F90
index cf54026e3..4eea7ccee 100644
--- a/src/SURFEX/init_watfluxn.F90
+++ b/src/SURFEX/init_watfluxn.F90
@@ -311,7 +311,8 @@ PTSURF(:) = WM%W%XTS(:)
!
CALL INIT_CHEMICAL_n(ILUOUT, KSV, HSV, WM%CHW%SVW, &
WM%CHW%CCH_NAMES, WM%CHW%CAER_NAMES, &
- HDSTNAMES=WM%CHW%CDSTNAMES, HSLTNAMES=WM%CHW%CSLTNAMES )
+ HDSTNAMES=WM%CHW%CDSTNAMES, HSLTNAMES=WM%CHW%CSLTNAMES, &
+ HSNWNAMES=WM%CHW%CSNWNAMES)
!
!* depositiion scheme
!
diff --git a/src/SURFEX/isba.F90 b/src/SURFEX/isba.F90
index 68d30f185..6705b544e 100644
--- a/src/SURFEX/isba.F90
+++ b/src/SURFEX/isba.F90
@@ -12,7 +12,8 @@
PPALPHAN, PZ0G_WITHOUT_SNOW, PZ0_MEBV, PZ0H_MEBV, &
PZ0EFF_MEBV, PZ0_MEBN, PZ0H_MEBN, PZ0EFF_MEBN, PTDEEP_A, &
PCSP, PFFG_NOSNOW, PFFV_NOSNOW, PEMIST, PUSTAR, PAC_AGG, &
- PHU_AGG, PRESP_BIOMASS_INST, PDEEP_FLUX, PIRRIG_GR )
+ PHU_AGG, PRESP_BIOMASS_INST, PDEEP_FLUX, PIRRIG_GR, &
+ PBLOWSNW_FLUX, PBLOWSNW_CONC )
! ##########################################################################
!
!
@@ -291,6 +292,16 @@ REAL, DIMENSION(:), INTENT(OUT) :: PDEEP_FLUX ! Heat flux at bottom of ISBA
!
REAL ,DIMENSION(:),INTENT(IN) :: PIRRIG_GR ! ground irrigation rate (kg/m2/s)
!
+!* Blowing snow variables
+! ----------------------
+!
+REAL, DIMENSION(:,:),OPTIONAL, INTENT(INOUT) :: PBLOWSNW_FLUX! Blowing snow particles flux:
+! 1: Number (#/m2/s) 2: Mass (kg/m2/s)
+! IN : contains sedimentation flux
+! OUT : contains emitted turbulent flux towards the atmosphere
+REAL, DIMENSION(:,:),OPTIONAL, INTENT(IN) :: PBLOWSNW_CONC ! Blowing snow particles concentration:
+! 1: Number (#/m3) 2: Mass (kg/m3)
+!
!
!* 0.2 declarations of local variables
!
@@ -482,7 +493,7 @@ IF(OMEB)THEN
PHU_AGG, PAC_AGG, ZDELHEATV_SFC, ZDELHEATG_SFC, ZDELHEATG, &
ZDELHEATN, ZDELHEATN_SFC, ZGSFCSNOW, PTDEEP_A, PDEEP_FLUX, &
ZRI3L, ZSNOW_THRUFAL, ZSNOW_THRUFAL_SOIL, ZEVAPCOR, ZSUBVCOR, &
- ZLITCOR, ZSNOWSFCH, ZQS3L )
+ ZLITCOR, ZSNOWSFCH, ZQS3L,PBLOWSNW_FLUX,PBLOWSNW_CONC )
ELSE
!
@@ -519,7 +530,8 @@ ELSE
PK%XDG, PK%XDZG, PPEW_A_COEF, PPEW_B_COEF, PPET_A_COEF, PPEQ_A_COEF, &
PPET_B_COEF, PPEQ_B_COEF, ZSNOW_THRUFAL_SOIL, ZGRNDFLUX, ZFLSN_COR, &
ZGSFCSNOW, ZEVAPCOR, ZLES3L, ZLEL3L, ZEVAP3L, ZSNOWSFCH, ZDELHEATN, &
- ZDELHEATN_SFC, ZRI3L, PZENITH, ZDELHEATG, ZDELHEATG_SFC, ZQS3L )
+ ZDELHEATN_SFC, ZRI3L, PZENITH, ZDELHEATG, ZDELHEATG_SFC, ZQS3L, &
+ PBLOWSNW_FLUX,PBLOWSNW_CONC )
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
!* 8.0 Plant stress, stomatal resistance and, possibly, CO2 assimilation
diff --git a/src/SURFEX/isba_meb.F90 b/src/SURFEX/isba_meb.F90
index b65e81214..78a0fb358 100644
--- a/src/SURFEX/isba_meb.F90
+++ b/src/SURFEX/isba_meb.F90
@@ -17,7 +17,7 @@
PHU_AGG, PAC_AGG, PDELHEATV_SFC, PDELHEATG_SFC, PDELHEATG,&
PDELHEATN, PDELHEATN_SFC, PRESTOREN, PTDEEP_A, PDEEP_FLUX,&
PRISNOW, PSNOW_THRUFAL, PSNOW_THRUFAL_SOIL, PEVAPCOR, &
- PSUBVCOR, PLITCOR, PSNOWSFCH, PQSNOW)
+ PSUBVCOR, PLITCOR, PSNOWSFCH, PQSNOW, PBLOWSNW_FLUX, PBLOWSNW_CONC )
! ##########################################################################
!
!
@@ -250,6 +250,17 @@ REAL, DIMENSION(:), INTENT(OUT) :: PQSNOW ! snow surface specific hum
!
REAL, DIMENSION(:,:), INTENT(OUT) :: PRESP_BIOMASS_INST ! instantaneous biomass respiration (kgCO2/kgair m/s)
!
+!
+!* Blowing snow variables
+! ----------------------
+!
+REAL, DIMENSION(:,:),OPTIONAL, INTENT(INOUT) :: PBLOWSNW_FLUX! Blowing snow particles flux:
+! 1: Number (#/m2/s) 2: Mass (kg/m2/s)
+! IN : contains sedimentation flux
+! OUT : contains emitted turbulent flux towards the atmosphere
+REAL, DIMENSION(:,:),OPTIONAL, INTENT(IN) :: PBLOWSNW_CONC ! Blowing snow particles concentration:
+! 1: Number (#/m3) 2: Mass (kg/m3)
+!
!* 0.2 declarations of local variables
!
!
@@ -922,7 +933,7 @@ ZVEGFACT(:) = ZSIGMA_F(:)*(1.0-PPALPHAN(:)*PEK%XPSN(:))
PPEQ_B_COEF, PSNOW_THRUFAL, PGRNDFLUX, PFLSN_COR, &
PRESTOREN, PEVAPCOR, DEK%XLES, DEK%XLESL, ZEVAP3L, PSNOWSFCH, &
PDELHEATN, PDELHEATN_SFC, PRISNOW, PZENITH, PDELHEATG, &
- PDELHEATG_SFC, PQSNOW )
+ PDELHEATG_SFC, PQSNOW,PBLOWSNW_FLUX,PBLOWSNW_CONC )
!
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!* 11.0 Litter layer hydrology:
diff --git a/src/SURFEX/modd_canopyn.F90 b/src/SURFEX/modd_canopyn.F90
index 0338174ec..abe67ec75 100644
--- a/src/SURFEX/modd_canopyn.F90
+++ b/src/SURFEX/modd_canopyn.F90
@@ -49,6 +49,8 @@ TYPE CANOPY_t
REAL, POINTER, DIMENSION(:,:) :: XLM ! Mixing lentgh (m)
REAL, POINTER, DIMENSION(:,:) :: XLEPS! Dissipative length (m)
REAL, POINTER, DIMENSION(:,:) :: XP ! pressure at each level in canopy (kg/m3)
+ REAL, POINTER, DIMENSION(:,:,:):: XBLOWSNW ! Blowing snow moments at each level in canopy
+ ! 1:( #/m3); 2:(kg/m3)
!
REAL, POINTER, DIMENSION(:,:) :: XDZ ! depth of each level in canopy (m)
REAL, POINTER, DIMENSION(:,:) :: XZF ! height of bottom of each level grid (m)
@@ -74,6 +76,7 @@ IF (LHOOK) CALL DR_HOOK("MODD_CANOPY_N:CANOPY_INIT",0,ZHOOK_HANDLE)
NULLIFY(SB%XDZ)
NULLIFY(SB%XZF)
NULLIFY(SB%XDZF)
+ NULLIFY(SB%XBLOWSNW)
SB%NLVL=0
IF (LHOOK) CALL DR_HOOK("MODD_CANOPY_N:CANOPY_INIT",1,ZHOOK_HANDLE)
END SUBROUTINE CANOPY_INIT
diff --git a/src/SURFEX/modd_ch_flaken.F90 b/src/SURFEX/modd_ch_flaken.F90
index da058239e..1797450b8 100644
--- a/src/SURFEX/modd_ch_flaken.F90
+++ b/src/SURFEX/modd_ch_flaken.F90
@@ -50,7 +50,7 @@ TYPE CH_FLAKE_t
CHARACTER(LEN=6), DIMENSION(:), POINTER :: CDSTNAMES
CHARACTER(LEN=6), DIMENSION(:), POINTER :: CSLTNAMES
CHARACTER(LEN=6), DIMENSION(:), POINTER :: CAER_NAMES
-
+ CHARACTER(LEN=6), DIMENSION(:), POINTER :: CSNWNAMES
!
END TYPE CH_FLAKE_t
@@ -68,6 +68,7 @@ NULLIFY(YCH_FLAKE%CCH_NAMES)
NULLIFY(YCH_FLAKE%CAER_NAMES)
NULLIFY(YCH_FLAKE%CDSTNAMES)
NULLIFY(YCH_FLAKE%CSLTNAMES)
+NULLIFY(YCH_FLAKE%CSNWNAMES)
YCH_FLAKE%CCH_DRY_DEP=' '
CALL SV_INIT(YCH_FLAKE%SVF)
IF (LHOOK) CALL DR_HOOK("MODD_CH_FLAKE_N:CH_FLAKE_INIT",1,ZHOOK_HANDLE)
diff --git a/src/SURFEX/modd_ch_isban.F90 b/src/SURFEX/modd_ch_isban.F90
index 7db03a80a..5d9eb5d3a 100644
--- a/src/SURFEX/modd_ch_isban.F90
+++ b/src/SURFEX/modd_ch_isban.F90
@@ -54,6 +54,7 @@ TYPE CH_ISBA_t
CHARACTER(LEN=6), DIMENSION(:), POINTER :: CCH_NAMES ! NAME OF CHEMICAL SPECIES
! (FOR DIAG ONLY)
CHARACTER(LEN=6), DIMENSION(:), POINTER :: CAER_NAMES ! NAME OF CHEMICAL SPECIES
+ CHARACTER(LEN=6), DIMENSION(:), POINTER :: CSNWNAMES ! NAME OF CHEMICAL SPECIES
CHARACTER(LEN=6), DIMENSION(:), POINTER :: CDSTNAMES ! NAME OF CHEMICAL SPECIES
CHARACTER(LEN=6), DIMENSION(:), POINTER :: CSLTNAMES ! NAME OF CHEMICAL SPECIES
!
@@ -78,6 +79,7 @@ NULLIFY(YCH_ISBA%CCH_NAMES)
NULLIFY(YCH_ISBA%CAER_NAMES)
NULLIFY(YCH_ISBA%CDSTNAMES)
NULLIFY(YCH_ISBA%CSLTNAMES)
+NULLIFY(YCH_ISBA%CSNWNAMES)
YCH_ISBA%CCHEM_SURF_FILE=' '
YCH_ISBA%CCH_DRY_DEP=' '
YCH_ISBA%LCH_BIO_FLUX=.FALSE.
diff --git a/src/SURFEX/modd_ch_seafluxn.F90 b/src/SURFEX/modd_ch_seafluxn.F90
index 8df308f63..177a206a4 100644
--- a/src/SURFEX/modd_ch_seafluxn.F90
+++ b/src/SURFEX/modd_ch_seafluxn.F90
@@ -48,7 +48,8 @@ TYPE CH_SEAFLUX_t
! (FOR DIAG ONLY)
CHARACTER(LEN=6), DIMENSION(:), POINTER :: CDSTNAMES
CHARACTER(LEN=6), DIMENSION(:), POINTER :: CSLTNAMES
- CHARACTER(LEN=6), DIMENSION(:), POINTER :: CAER_NAMES
+ CHARACTER(LEN=6), DIMENSION(:), POINTER :: CAER_NAMES
+ CHARACTER(LEN=6), DIMENSION(:), POINTER :: CSNWNAMES
!
END TYPE CH_SEAFLUX_t
@@ -70,6 +71,7 @@ NULLIFY(YCH_SEAFLUX%CCH_NAMES)
NULLIFY(YCH_SEAFLUX%CAER_NAMES)
NULLIFY(YCH_SEAFLUX%CDSTNAMES)
NULLIFY(YCH_SEAFLUX%CSLTNAMES)
+NULLIFY(YCH_SEAFLUX%CSNWNAMES)
YCH_SEAFLUX%CCH_DRY_DEP=' '
CALL SV_INIT(YCH_SEAFLUX%SVS)
IF (LHOOK) CALL DR_HOOK("MODD_CH_SEAFLUX_N:CH_SEAFLUX_INIT",1,ZHOOK_HANDLE)
diff --git a/src/SURFEX/modd_ch_watfluxn.F90 b/src/SURFEX/modd_ch_watfluxn.F90
index cfb2777d0..64cc89ed8 100644
--- a/src/SURFEX/modd_ch_watfluxn.F90
+++ b/src/SURFEX/modd_ch_watfluxn.F90
@@ -49,6 +49,7 @@ TYPE CH_WATFLUX_t
CHARACTER(LEN=6), DIMENSION(:), POINTER :: CDSTNAMES
CHARACTER(LEN=6), DIMENSION(:), POINTER :: CSLTNAMES
CHARACTER(LEN=6), DIMENSION(:), POINTER :: CAER_NAMES
+ CHARACTER(LEN=6), DIMENSION(:), POINTER :: CSNWNAMES
!
END TYPE CH_WATFLUX_t
!
@@ -63,6 +64,7 @@ NULLIFY(YCH_WATFLUX%CCH_NAMES)
NULLIFY(YCH_WATFLUX%CAER_NAMES)
NULLIFY(YCH_WATFLUX%CDSTNAMES)
NULLIFY(YCH_WATFLUX%CSLTNAMES)
+NULLIFY(YCH_WATFLUX%CSNWNAMES)
YCH_WATFLUX%CCH_DRY_DEP=' '
CALL SV_INIT(YCH_WATFLUX%SVW)
IF (LHOOK) CALL DR_HOOK("MODD_CH_WATFLUX_N:CH_WATFLUX_INIT",1,ZHOOK_HANDLE)
diff --git a/src/SURFEX/modd_surfexn.F90 b/src/SURFEX/modd_surfexn.F90
index 56b37e74f..e76ee83a5 100644
--- a/src/SURFEX/modd_surfexn.F90
+++ b/src/SURFEX/modd_surfexn.F90
@@ -7,6 +7,7 @@ MODULE MODD_SURFEX_n
USE MODD_AGRI_n, ONLY : AGRI_NP_t
USE MODD_BEM_OPTION_n, ONLY : BEM_OPTIONS_t
USE MODD_BLD_DESCRIPTION_n, ONLY : BLD_DESC_t
+USE MODD_BLOWSNW_n, ONLY : BLOWSNW_t
USE MODD_CH_EMIS_FIELD_n, ONLY : CH_EMIS_FIELD_t
USE MODD_CH_FLAKE_n, ONLY : CH_FLAKE_t
USE MODD_CH_ISBA_n, ONLY : CH_ISBA_t, CH_ISBA_NP_t
@@ -259,6 +260,7 @@ TYPE(CH_EMIS_SNAP_t) :: CHN
TYPE(EMIS_GR_FIELD_t) :: EGF
TYPE(DST_NP_t) :: NDST
TYPE(SLT_t) :: SLT
+TYPE(BLOWSNW_t) :: BLOWSNW
!
TYPE(FLAKE_MODEL_t) :: FM
TYPE(WATFLUX_MODEL_t) :: WM
diff --git a/src/SURFEX/modd_svn.F90 b/src/SURFEX/modd_svn.F90
index ab64b2e34..a554e0bcd 100644
--- a/src/SURFEX/modd_svn.F90
+++ b/src/SURFEX/modd_svn.F90
@@ -46,8 +46,11 @@ TYPE SV_t
INTEGER :: NSV_SLTBEG, NSV_SLTEND ! index of first and last sea salt related scalar variable
INTEGER :: NSLTEQ ! number of sea salt related species in scalar variables list
INTEGER :: NSV_AERBEG, NSV_AEREND ! index of first and last aerosol related scalar variabl
- INTEGER :: NAEREQ ! number of aerosols variables
-
+ INTEGER :: NAEREQ ! number of aerosols variables$
+ INTEGER :: NSV_SNWBEG, NSV_SNWEND ! index of first and last blowing snow related scalar variable
+ INTEGER :: NSNWEQ ! number of blowing snow related species in scalar variables list
+ INTEGER :: N2D_SNWBEG, N2D_SNWEND ! index of first and last blowing snow 2D variable sent to MNH
+ INTEGER :: N2DSNWEQ ! number of blowing snow 2D related species in scalar variables list
!
!
END TYPE SV_t
@@ -78,6 +81,12 @@ YSV%NSLTEQ=0
YSV%NSV_AERBEG=0
YSV%NSV_AEREND=0
YSV%NAEREQ=0
+YSV%NSV_SNWBEG=0
+YSV%NSV_SNWEND=0
+YSV%NSNWEQ=0
+YSV%N2D_SNWBEG=0
+YSV%N2D_SNWEND=0
+YSV%N2DSNWEQ=0
IF (LHOOK) CALL DR_HOOK("MODD_SV_N:SV_INIT",1,ZHOOK_HANDLE)
END SUBROUTINE SV_INIT
diff --git a/src/SURFEX/read_all_namelists.F90 b/src/SURFEX/read_all_namelists.F90
index e1f0a10fb..e2a4dd4fd 100644
--- a/src/SURFEX/read_all_namelists.F90
+++ b/src/SURFEX/read_all_namelists.F90
@@ -25,6 +25,7 @@ USE MODI_READ_NAMELISTS_IDEAL_n
USE MODI_READ_NAMELISTS_IDEAL
USE MODI_READ_NAMELISTS_SLT
USE MODI_READ_NAMELISTS_DST
+USE MODI_READ_NAMELISTS_BLOWSNW
USE MODI_READ_NAMELISTS_ASSIM
USE MODI_READ_NAMELISTS_TOPD
#ifdef SFX_ARO
@@ -68,6 +69,7 @@ LNAM_READ=.TRUE.
CALL READ_NAMELISTS_ISBA(HPROGRAM)
CALL READ_NAMELISTS_SLT(HPROGRAM)
CALL READ_NAMELISTS_DST(HPROGRAM)
+ CALL READ_NAMELISTS_BLOWSNW(HPROGRAM)
CALL READ_NAMELISTS_IDEAL(HPROGRAM)
CALL READ_NAMELISTS_ASSIM(HPROGRAM)
CALL READ_NAMELISTS_TOPD(HPROGRAM)
diff --git a/src/SURFEX/read_sbln.F90 b/src/SURFEX/read_sbln.F90
index 500e64ae7..07776ff4a 100644
--- a/src/SURFEX/read_sbln.F90
+++ b/src/SURFEX/read_sbln.F90
@@ -3,7 +3,7 @@
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!SFX_LIC for details. version 1.
! #########
- SUBROUTINE READ_SBL_n (DTCO, U, SB, OSBL, HPROGRAM, HSURF)
+ SUBROUTINE READ_SBL_n (DTCO, U, SB, OSBL, HPROGRAM, HSURF, SV, BLOWSNW)
! #########################################
!
!!**** *READ_SBL_n* - reads TEB fields
@@ -45,6 +45,8 @@
USE MODD_DATA_COVER_n, ONLY : DATA_COVER_t
USE MODD_SURF_ATM_n, ONLY : SURF_ATM_t
USE MODD_CANOPY_n, ONLY : CANOPY_t
+USE MODD_SV_n, ONLY : SV_t
+USE MODD_BLOWSNW_n, ONLY : BLOWSNW_t
!
USE MODD_SURF_PAR, ONLY : XUNDEF
!
@@ -64,6 +66,8 @@ IMPLICIT NONE
TYPE(DATA_COVER_t), INTENT(INOUT) :: DTCO
TYPE(SURF_ATM_t), INTENT(INOUT) :: U
TYPE(CANOPY_t), INTENT(INOUT) :: SB
+TYPE(SV_t), INTENT(IN), OPTIONAL :: SV
+TYPE(BLOWSNW_t), INTENT(INOUT),OPTIONAL :: BLOWSNW
LOGICAL, INTENT(INOUT) :: OSBL
!
CHARACTER(LEN=6), INTENT(IN) :: HPROGRAM ! calling program
@@ -79,6 +83,7 @@ LOGICAL, INTENT(INOUT) :: OSBL
CHARACTER(LEN=13) :: YFORMAT
CHARACTER(LEN=3) :: YREAD
INTEGER :: JLAYER ! loop counter on layers
+INTEGER :: JBLOWSNW ! loop counter on blowing snow variables
INTEGER :: ILU ! 1D physical dimension
INTEGER :: IRESP ! Error code after redding
INTEGER :: IVERSION, IBUGFIX ! surface version
@@ -129,6 +134,9 @@ IF (.NOT.OSBL) THEN
ALLOCATE(SB%XDZ (0,0))
ALLOCATE(SB%XZF (0,0))
ALLOCATE(SB%XDZF(0,0))
+ IF (HSURF=="NATURE") THEN
+ ALLOCATE(SB%XBLOWSNW(0,0,0))
+ ENDIF
IF (LHOOK) CALL DR_HOOK('READ_SBL_N',1,ZHOOK_HANDLE)
RETURN
ENDIF
@@ -173,6 +181,19 @@ ALLOCATE(SB%XTKE(ILU,SB%NLVL))
ALLOCATE(SB%XLMO(ILU,SB%NLVL))
ALLOCATE(SB%XP (ILU,SB%NLVL))
!
+IF (HSURF=="NATURE") THEN
+IF(SV%NSNWEQ>0) THEN
+ ALLOCATE(SB%XBLOWSNW(ILU,SB%NLVL,SV%NSNWEQ))
+ ! Allocate diagnostic
+ ALLOCATE (BLOWSNW%XSNW_CANO_RGA(ILU,SB%NLVL)) !Output array
+ ALLOCATE (BLOWSNW%XSNW_CANO_VAR(ILU,SB%NLVL,SV%NSNWEQ))
+ELSE
+ ALLOCATE(SB%XBLOWSNW (0,0,0))
+ ALLOCATE(BLOWSNW%XSNW_CANO_RGA(0,0))
+ ALLOCATE(BLOWSNW%XSNW_CANO_VAR(0,0,0))
+END IF
+ENDIF
+!
IF (IVERSION>7 .OR. IVERSION==7 .AND.IBUGFIX>=2) THEN
YRECFM='STORAGETYPE'
CALL READ_SURF(HPROGRAM,YRECFM,YREAD,IRESP)
@@ -226,6 +247,20 @@ IF(YREAD=='ALL') THEN
CALL READ_SURF(HPROGRAM,YRECFM,SB%XP(:,JLAYER),IRESP)
END DO
!
+ IF (HSURF=="NATURE") THEN
+ !* Blowing snow number and mass concentration
+ IF(SV%NSNWEQ>0) THEN
+ DO JBLOWSNW=1,(SV%NSNWEQ)
+ DO JLAYER=1,SB%NLVL
+ WRITE(YRECFM,'(A8,I1.1,A1,I2.2,A4)') 'CANSNW_M',JBLOWSNW,'L',JLAYER,' '
+ CALL READ_SURF(HPROGRAM,YRECFM,SB%XBLOWSNW(:,JLAYER,JBLOWSNW),IRESP)
+ END DO
+ END DO
+ BLOWSNW%XSNW_CANO_RGA (:,:) = 0.
+ BLOWSNW%XSNW_CANO_VAR (:,:,:) = 0.
+ ENDIF
+ ENDIF
+ !
ELSE
SB%XU (:,:) = XUNDEF
SB%XT (:,:) = XUNDEF
@@ -233,6 +268,9 @@ ELSE
SB%XTKE(:,:) = XUNDEF
SB%XLMO(:,:) = XUNDEF
SB%XP (:,:) = XUNDEF
+ IF (HSURF=="NATURE") THEN
+ SB%XBLOWSNW (:,:,:) = XUNDEF
+ ENDIF
ENDIF
!
IF (HSURF=="TOWN ") THEN
diff --git a/src/SURFEX/snow3L_isba.F90 b/src/SURFEX/snow3L_isba.F90
index adc51bec5..72d786834 100644
--- a/src/SURFEX/snow3L_isba.F90
+++ b/src/SURFEX/snow3L_isba.F90
@@ -11,7 +11,8 @@ SUBROUTINE SNOW3L_ISBA(IO, G, PK, PEK, DK, DEK, DMK, OMEB, HIMPLICIT_WIND,
PPEQ_A_COEF, PPET_B_COEF, PPEQ_B_COEF, PTHRUFAL, &
PGRNDFLUX, PFLSN_COR, PGSFCSNOW, PEVAPCOR, PLES3L, PLEL3L,&
PEVAP, PSNOWSFCH, PDELHEATN, PDELHEATN_SFC, PRI, PZENITH, &
- PDELHEATG, PDELHEATG_SFC, PQS )
+ PDELHEATG, PDELHEATG_SFC, PQS , &
+ PBLOWSNW_FLUX,PBLOWSNW_CONC)
! ######################################################################################
!
!!**** *SNOW3L_ISBA*
@@ -84,9 +85,12 @@ USE MODD_DATA_COVER_PAR, ONLY : NVT_SNOW, &
NVT_TRBD, NVT_TEBE, NVT_TENE, &
NVT_BOBD, NVT_BOND, NVT_SHRB
!
+USE MODD_BLOWSNW_SURF
+!
USE MODI_SNOW3L
USE MODI_SNOWCRO
USE MODI_SNOWCRO_DIAG
+USE MODI_SNOWPACK_EVOL
!
#ifdef SFX_OL
USE MODN_IO_OFFLINE, ONLY : XTSTEP_OUTPUT
@@ -207,6 +211,15 @@ REAL, DIMENSION(:), INTENT(OUT) :: PQS
! puis plus tard dans LALB
REAL, DIMENSION(:), INTENT(IN) :: PZENITH ! solar zenith angle
!
+REAL, DIMENSION(:,:), INTENT(INOUT) :: PBLOWSNW_FLUX
+! PBLOWSNW_FLUX = Blowing snow particles flux:
+! 1: Number (#/m2/s) 2: Mass (kg/m2/s)
+! IN : contains sedimentation flux
+! OUT : contains emitted turbulent flux towards the atmosphere
+REAL, DIMENSION(:,:), INTENT(IN) :: PBLOWSNW_CONC
+! PBLOWSNW_CONC = Blowing snow particles concentration:
+! 1: Number (#/m3) 2: Mass (kg/m3)
+!
!* 0.2 declarations of local variables
!
REAL, PARAMETER :: ZCHECK_TEMP = 50.0
@@ -216,6 +229,7 @@ INTEGER :: JWRK, JJ ! Loop control
!
INTEGER :: INLVLS ! maximum number of snow layers
INTEGER :: INLVLG ! number of ground layers
+INTEGER :: IBLOWSNW ! number of blowing snow variables
!
REAL, DIMENSION(SIZE(PTA)) :: ZRRSNOW, ZSOILCOND, ZSNOW, ZSNOWFALL, &
ZSNOWABLAT_DELTA, ZSNOWSWE_1D, ZSNOWD, &
@@ -247,6 +261,11 @@ REAL, DIMENSION(SIZE(PTA)) :: ZRRSNOW, ZSOILCOND, ZSNOW, ZSNOWFALL, &
! of total ablation during a timestep (-).
! ZWORK = local working variable (*)
! ZC2 = sub-surface heat capacity [(K m2)/J]
+REAL, DIMENSION(SIZE(PTA)) :: ZBLOWSNW_ACC,ZBLOWSNW_DEPFLUX
+! ZBLOWSNW_ACC = minimum equivalent snow depth
+! for deposition of blown snow particles
+! during the current time step (m)
+! ZBLOWSNW_DEPFLUX = deposition flux of blowing snow (kg/m2/s)
!
!* 0.3 declarations of packed variables
!
@@ -299,6 +318,8 @@ ZSOILCOND(:) = 0.0
ZRRSNOW(:) = 0.0
ZSNOWFALL(:) = 0.0
ZSNOWABLAT_DELTA(:) = 0.0
+ZBLOWSNW_ACC(:) = 0.0
+ZBLOWSNW_DEPFLUX(:) = 0.0
!
ZWGHT(:) = 0.0
ZWORK(:) = 0.0
@@ -309,6 +330,7 @@ DMK%XSNOWDZ(:,:) = 0.0
!
INLVLS = SIZE(PEK%TSNOW%WSNOW(:,:),2)
INLVLG = MIN(SIZE(PD_G(:,:),2),SIZE(PTG(:,:),2))
+IBLOWSNW = SIZE(PBLOWSNW_FLUX(:,:),2)
!
!
IF(.NOT.OMEB)THEN
@@ -359,6 +381,19 @@ IF (PEK%TSNOW%SCHEME=='3-L' .OR. IO%CISBA == 'DIF' .OR. PEK%TSNOW%SCHEME == 'CRO
ZSNOWFALL(JJ) = PSR(JJ)*PTSTEP/XRHOSMAX_ES ! maximum possible snowfall depth (m)
ENDDO
!
+! Calculate maximum deposited snow depth (m) of blown snow particles
+!
+!
+ IF(IBLOWSNW>0.)THEN
+ DO JJ=1,SIZE(PSR)
+ ZBLOWSNW_ACC(JJ)=(PBLOWSNW_FLUX(JJ,2)+PBLOWSNW_FLUX(JJ,3))*PTSTEP/XRHO_DEP
+ IF (PBLOWSNW_FLUX(JJ,2)+PBLOWSNW_FLUX(JJ,3)>0.) THEN
+ ZBLOWSNW_DEPFLUX(JJ) = (PBLOWSNW_FLUX(JJ,2)+PBLOWSNW_FLUX(JJ,3))
+ ENDIF
+ ENDDO
+ ENDIF
+!
+!
! Calculate preliminary snow depth (m)
ZSNOW(:) =0.
@@ -400,13 +435,13 @@ IF (PEK%TSNOW%SCHEME=='3-L' .OR. IO%CISBA == 'DIF' .OR. PEK%TSNOW%SCHEME == 'CRO
NMASK(:) = 0
!
DO JJ=1,SIZE(ZSNOW)
- IF (ZSNOW(JJ) >= XSNOWDMIN .OR. ZSNOWFALL(JJ) >= XSNOWDMIN) THEN
+ IF (ZSNOW(JJ) >= XSNOWDMIN .OR. ZSNOWFALL(JJ) >= XSNOWDMIN .OR. ZBLOWSNW_ACC(JJ) >= XSNOWDMIN) THEN
ISIZE_SNOW = ISIZE_SNOW + 1
NMASK(ISIZE_SNOW) = JJ
ENDIF
ENDDO
!
- IF (ISIZE_SNOW>0) CALL CALL_MODEL(ISIZE_SNOW,INLVLS,INLVLG,NMASK)
+ IF (ISIZE_SNOW>0) CALL CALL_MODEL(ISIZE_SNOW,INLVLS,INLVLG,IBLOWSNW,NMASK)
!
! ===============================================================
!
@@ -448,7 +483,7 @@ IF (PEK%TSNOW%SCHEME=='3-L' .OR. IO%CISBA == 'DIF' .OR. PEK%TSNOW%SCHEME == 'CRO
WHERE(LREMOVE_SNOW(:))
!
ZSNOWSWE_OUT(:) = 0.0
- PLES3L(:) = MIN(PLES3L(:), XLSTT*(ZSNOWSWE_1D(:)/PTSTEP + PSR(:)))
+ PLES3L(:) = MIN(PLES3L(:), XLSTT*(ZSNOWSWE_1D(:)/PTSTEP + PSR(:)+ZBLOWSNW_DEPFLUX(:)))
PLEL3L(:) = 0.0
PEVAP(:) = PLES3L(:)/PK%XLSTT(:)
PTHRUFAL(:) = MAX(0.0, ZSNOWSWE_1D(:)/PTSTEP + PSR(:) - PEVAP(:)*ZPSN(:) + ZRRSNOW(:)) ! kg m-2 s-1
@@ -584,13 +619,14 @@ IF (LHOOK) CALL DR_HOOK('SNOW3L_ISBA',1,ZHOOK_HANDLE)
CONTAINS
!
!================================================================
-SUBROUTINE CALL_MODEL(KSIZE1,KSIZE2,KSIZE3,KMASK)
+SUBROUTINE CALL_MODEL(KSIZE1,KSIZE2,KSIZE3,KSIZE4,KMASK)
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: KSIZE1
INTEGER, INTENT(IN) :: KSIZE2
INTEGER, INTENT(IN) :: KSIZE3
+INTEGER, INTENT(IN) :: KSIZE4
INTEGER, DIMENSION(:), INTENT(IN) :: KMASK
!
REAL, DIMENSION(KSIZE1,KSIZE2) :: ZP_SNOWSWE
@@ -603,6 +639,8 @@ REAL, DIMENSION(KSIZE1,KSIZE2) :: ZP_SNOWGRAN1
REAL, DIMENSION(KSIZE1,KSIZE2) :: ZP_SNOWGRAN2
REAL, DIMENSION(KSIZE1,KSIZE2) :: ZP_SNOWHIST
REAL, DIMENSION(KSIZE1,KSIZE2) :: ZP_SNOWAGE
+REAL, DIMENSION(KSIZE1,KSIZE4) :: ZP_BLOWSNW_FLUX
+REAL, DIMENSION(KSIZE1,KSIZE4-1):: ZP_BLOWSNW_CONC
REAL, DIMENSION(KSIZE1) :: ZP_SNOWALB
REAL, DIMENSION(KSIZE1) :: ZP_SWNETSNOW
REAL, DIMENSION(KSIZE1) :: ZP_SWNETSNOWS
@@ -675,6 +713,7 @@ REAL, DIMENSION(KSIZE1) :: ZP_LAT,ZP_LON
REAL, DIMENSION(KSIZE1) :: ZP_PSN_INV
REAL, DIMENSION(KSIZE1) :: ZP_PSN
REAL, DIMENSION(KSIZE1) :: ZP_PSN_GFLXCOR
+REAL, DIMENSION(KSIZE1) :: ZP_BLOWSNW_DEP
REAL, DIMENSION(KSIZE1) :: ZP_WORK
!
REAL, DIMENSION(KSIZE1,KSIZE2) :: ZP_SNOWDEND
@@ -733,6 +772,36 @@ IF (PEK%TSNOW%SCHEME=='CRO') THEN
ZP_SNOWHIST (JJ,JWRK) = PEK%TSNOW%HIST (JI,JWRK)
ENDDO
ENDDO
+
+ IF(KSIZE4>0) THEN
+ DO JWRK=1,KSIZE4
+ DO JJ=1,KSIZE1
+ JI = KMASK(JJ)
+ ZP_BLOWSNW_FLUX(JJ,JWRK) = PBLOWSNW_FLUX(JI,JWRK)
+ ENDDO
+ ENDDO
+ DO JWRK=1,KSIZE4-1
+ DO JJ=1,KSIZE1
+ JI = KMASK(JJ)
+ ZP_BLOWSNW_CONC(JJ,JWRK) = PBLOWSNW_CONC(JI,JWRK)
+ ENDDO
+ ENDDO
+ ELSE
+ DO JWRK=1,KSIZE4
+ DO JJ=1,KSIZE1
+ ZP_BLOWSNW_FLUX(JJ,JWRK) = XUNDEF
+ ENDDO
+ ENDDO
+ DO JWRK=1,KSIZE4-1
+ DO JJ=1,KSIZE1
+ ZP_BLOWSNW_CONC(JJ,JWRK) = XUNDEF
+ ENDDO
+ ENDDO
+ DO JJ=1,KSIZE1
+ ZP_BLOWSNW_DEP(JJ) = 0.
+ ENDDO
+ ENDIF
+
ELSE
DO JWRK=1,KSIZE2
DO JJ=1,KSIZE1
@@ -741,6 +810,16 @@ ELSE
ZP_SNOWHIST (JJ,JWRK) = XUNDEF
ENDDO
ENDDO
+
+ IF(KSIZE4>0) THEN
+ DO JWRK=1,KSIZE4
+ DO JJ=1,KSIZE1
+ ZP_BLOWSNW_CONC(JJ,JWRK) = XUNDEF
+ ZP_BLOWSNW_FLUX(JJ,JWRK) = XUNDEF
+ ENDDO
+ ENDDO
+ ENDIF
+
ENDIF
!
DO JWRK=1,KSIZE3
@@ -883,6 +962,24 @@ ENDIF
! Call ISBA-SNOW3L model:
!
IF (PEK%TSNOW%SCHEME=='CRO') THEN
+!
+! ------------------------
+! Blowing snow scheme (when coupled to Meso-NH)
+! Handel snow Erosion and deposition
+!
+ IF(KSIZE4>0)THEN
+ CALL SNOWPACK_EVOL(IO%CSNOWRES,ZP_BLOWSNW_FLUX,ZP_SNOWHEAT, &
+ ZP_SNOWSWE,ZP_SNOWRHO, &
+ ZP_SNOWGRAN1,ZP_SNOWGRAN2,ZP_SNOWHIST, &
+ ZP_SNOWAGE, PTSTEP,ZP_RHOA,ZP_TA, &
+ ZP_BLOWSNW_CONC,ZP_VMOD, ZP_QA,ZP_PS, &
+ ZP_UREF,ZP_EXNS,ZP_DIRCOSZW, &
+ ZP_ZREF,ZP_Z0EFF,ZP_Z0HNAT,ZP_BLOWSNW_DEP, &
+ ZP_TG(:,1))
+ ENDIF
+!
+! ------------------------
+! Main call to Crocus
CALL SNOWCRO(IO%CSNOWRES, TPTIME, IO%LGLACIER, HIMPLICIT_WIND, &
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_SNOWSWE, ZP_SNOWRHO, &
@@ -896,7 +993,7 @@ IF (PEK%TSNOW%SCHEME=='CRO') THEN
ZP_HSNOW, ZP_GFLUXSNOW, ZP_HPSNOW, ZP_LES3L, ZP_LEL3L, &
ZP_EVAP, ZP_SNDRIFT, ZP_RI,ZP_EMISNOW, ZP_CDSNOW, &
ZP_USTARSNOW, ZP_CHSNOW, ZP_SNOWHMASS, ZP_QS, ZP_VEGTYPE, &
- ZP_ZENITH, ZP_LAT, ZP_LON, IO%LSNOWDRIFT, &
+ ZP_ZENITH, ZP_LAT, ZP_LON,ZP_BLOWSNW_DEP, IO%LSNOWDRIFT, &
IO%LSNOWDRIFT_SUBLIM, IO%LSNOW_ABS_ZENITH, IO%CSNOWMETAMO, &
IO%CSNOWRAD )
!
@@ -1072,6 +1169,15 @@ IF (PEK%TSNOW%SCHEME=='CRO') THEN
ENDDO
ENDDO
+ IF(KSIZE4>0.)THEN
+ DO JWRK=1,KSIZE4
+ DO JJ=1,KSIZE1
+ JI = KMASK(JJ)
+ PBLOWSNW_FLUX(JI,JWRK) = ZP_BLOWSNW_FLUX(JJ,JWRK)
+ ENDDO
+ ENDDO
+ ENDIF
+
IF (SIZE(DMK%XSNOWDEND)>0) THEN
! This is equivalent to test the value of DGMI%LPROSNOW which does not enter in ISBA
DO JWRK = 1,KSIZE2
diff --git a/src/SURFEX/snowcro.F90 b/src/SURFEX/snowcro.F90
index 7c52c9cd7..73fec9b16 100644
--- a/src/SURFEX/snowcro.F90
+++ b/src/SURFEX/snowcro.F90
@@ -17,7 +17,7 @@
PTHRUFAL,PGRNDFLUX,PEVAPCOR,PRNSNOW,PHSNOW,PGFLUXSNOW, &
PHPSNOW,PLES3L,PLEL3L,PEVAP,PSNDRIFT,PRI, &
PEMISNOW,PCDSNOW,PUSTAR,PCHSNOW,PSNOWHMASS,PQS, &
- PPERMSNOWFRAC,PZENITH,PXLAT,PXLON, &
+ PPERMSNOWFRAC,PZENITH,PXLAT,PXLON,PBLOWSNW_DEP, &
OSNOWDRIFT,OSNOWDRIFT_SUBLIM,OSNOW_ABS_ZENITH, &
HSNOWMETAMO,HSNOWRAD)
! ##########################################################################
@@ -277,6 +277,9 @@ REAL, DIMENSION(:), INTENT(OUT) :: PRI, PQS
REAL, DIMENSION(:), INTENT(IN) :: PZENITH ! solar zenith angle
REAL, DIMENSION(:), INTENT(IN) :: PXLAT,PXLON ! LAT/LON after packing
!
+REAL, DIMENSION(:), INTENT(IN) :: PBLOWSNW_DEP ! Blowing snow deposition rate (kg/m2/s)
+! (when coupled with MNH)
+!
LOGICAL, INTENT(IN) :: OSNOWDRIFT, OSNOWDRIFT_SUBLIM ! activate snowdrift, sublimation during drift
LOGICAL, INTENT(IN) :: OSNOW_ABS_ZENITH ! activate parametrization of solar absorption for polar regions
CHARACTER(3), INTENT(IN) :: HSNOWMETAMO, HSNOWRAD
@@ -548,7 +551,7 @@ ZSNOWBIS(:) = ZSNOW(:)
PSNOWGRAN1,PSNOWGRAN2,GSNOWFALL,ZSNOWDZN, &
ZSNOWRHOF,ZSNOWDZF,ZSNOWGRAN1F,ZSNOWGRAN2F, ZSNOWHISTF, &
ZSNOWAGEF,GMODIF_MAILLAGE,INLVLS_USE,OSNOWDRIFT,PZ0EFF,PUREF,&
- HSNOWMETAMO)
+ HSNOWMETAMO,PBLOWSNW_DEP)
!
!***************************************DEBUG IN**********************************************
IF (GCRODEBUGDETAILSPRINT) THEN
@@ -3595,7 +3598,7 @@ SUBROUTINE SNOWNLFALL_UPGRID(TPTIME, OGLACIER,PTSTEP,PSR,PTA,PVMOD, &
GSNOWFALL,PSNOWDZN,PSNOWRHOF,PSNOWDZF, &
PSNOWGRAN1F,PSNOWGRAN2F,PSNOWHISTF,PSNOWAGEF, &
OMODIF_GRID,KNLVLS_USE,OSNOWDRIFT,PZ0EFF,PUREF,&
- HSNOWMETAMO)
+ HSNOWMETAMO,PBLOWSNW_DEP)
!
!! PURPOSE
!! -------
@@ -3630,6 +3633,8 @@ USE MODD_SNOW_PAR, ONLY : XRHOSMIN_ES, XSNOWDMIN, XANSMAX, XAGLAMAX, XSNOWCRITD,
XDEPTH_SURFACE, XDIFF_1, XDIFF_MAX, XSCALE_DIFF, &
XSNOWFALL_A_SN, XSNOWFALL_B_SN, XSNOWFALL_C_SN
!
+USE MODD_BLOWSNW_SURF
+!
USE MODE_SNOW3L
!
IMPLICIT NONE
@@ -3671,6 +3676,9 @@ INTEGER, DIMENSION(:), INTENT(INOUT) :: KNLVLS_USE
LOGICAL,INTENT(IN) :: OSNOWDRIFT ! if snowdrift then grain types are not modified by wind
CHARACTER(3), INTENT(IN) :: HSNOWMETAMO ! metamorphism scheme
+!
+REAL, DIMENSION(:), INTENT(IN) :: PBLOWSNW_DEP ! Blowing snow deposition flux (kg/m2/s)
+!
!* 0.2 declarations of local variables
!
!
@@ -3710,6 +3718,10 @@ REAL, PARAMETER :: PPHREF_WIND_MIN = MIN(PPHREF_WIND_RHO,PPHR
REAL, DIMENSION(SIZE(PTA)) :: ZWIND_RHO
REAL, DIMENSION(SIZE(PTA)) :: ZWIND_GRAIN
!
+REAL, DIMENSION(SIZE(PTA)) :: ZRHO_BS ! Density of deposited blowing snow
+!
+REAL :: ZMOB, ZG1_BS, ZG2_BS
+!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!* 1.0 Initialization and snowage calculation for the present date
@@ -3916,7 +3928,7 @@ END DO
!!
DO JJ = 1,SIZE(PSNOW(:))
!
- IF ( PSR(JJ)>0.0 ) THEN
+ IF ( PSR(JJ)>0.0 .OR. PBLOWSNW_DEP(JJ) > XUEPSI ) THEN
!
! newly fallen snow characteristics:
IF ( KNLVLS_USE(JJ)>0 ) THEN !Case of new snowfall on a previously snow-free surface
@@ -3941,12 +3953,21 @@ DO JJ = 1,SIZE(PSNOW(:))
ZWIND_GRAIN(JJ) = PVMOD(JJ)*LOG(PPHREF_WIND_GRAIN/ZZ0EFF)/ &
LOG(PUREF(JJ)/ZZ0EFF)
- PSNOWHMASS(JJ) = PSR(JJ) * ( XCI * ( ZSNOWTEMP(JJ)-XTT ) - XLMTT ) * PTSTEP
+ PSNOWHMASS(JJ) = (PSR(JJ)+PBLOWSNW_DEP(JJ)) * ( XCI * ( ZSNOWTEMP(JJ)-XTT ) - XLMTT ) * PTSTEP
!
PSNOWRHOF (JJ) = MAX( XRHOSMIN_ES, XSNOWFALL_A_SN + &
XSNOWFALL_B_SN * ( PTA(JJ)-XTT ) + &
XSNOWFALL_C_SN * MIN( PVMOD(JJ), SQRT(ZWIND_RHO(JJ) ) ) )
- ZSNOWFALL (JJ) = PSR(JJ) * PTSTEP / PSNOWRHOF(JJ) ! snowfall thickness (m)
+!
+! Density of deposited blown snow particles (can be adjusted depending on model application)
+! 1st option : same as fresh snow
+ ZRHO_BS(JJ) = PSNOWRHOF(JJ)
+
+! Density of accumulated snow (falling+blowing snow) : weighted average of PSNOWRHOF and ZRHO_BS
+ PSNOWRHOF(JJ) = (PSNOWRHOF(JJ)*PSR(JJ) + ZRHO_BS(JJ) * PBLOWSNW_DEP(JJ))/ &
+ (PSR(JJ)+PBLOWSNW_DEP(JJ))
+!
+ ZSNOWFALL (JJ) = (PSR(JJ)+PBLOWSNW_DEP(JJ)) * PTSTEP / PSNOWRHOF(JJ) ! snowfall thickness (m)
PSNOW (JJ) = PSNOW(JJ) + ZSNOWFALL(JJ)
PSNOWDZF (JJ) = ZSNOWFALL(JJ)
!
@@ -3960,6 +3981,25 @@ DO JJ = 1,SIZE(PSNOW(:))
PSNOWGRAN2F(JJ) = MIN( MAX( XNSPH1*ZWIND_GRAIN(JJ)+XNSPH2, XNSPH3 ), XNSPH4 )
END IF
!
+! --------------------------------------------------------------------
+! Characteristic of deposited blown snow particles (dendricity and sphericity)
+! ----------------------------------------------------------------------
+! 1st Option : same as falling snow flakes
+ ZG1_BS = PSNOWGRAN1F(JJ)
+ ZG2_BS = PSNOWGRAN2F(JJ)
+!
+! 2nd Option : fixed characteristics
+! ZG1_BS = XSDEP_GRA1
+! ZG2_BS = XSDEP_GRA2
+!
+! --------------------------------------------------------------------
+! Characteristic of accumulated snow grains (falling snow + deposited blowing
+! snow) : weighted average
+! ----------------------------------------------------------------------
+ PSNOWGRAN1F(JJ) = (PSNOWGRAN1F(JJ)*PSR(JJ) + ZG1_BS * PBLOWSNW_DEP(JJ))/ &
+ (PSR(JJ)+PBLOWSNW_DEP(JJ))
+ PSNOWGRAN2F(JJ) = (PSNOWGRAN2F(JJ)*PSR(JJ) + ZG2_BS * PBLOWSNW_DEP(JJ))/ &
+ (PSR(JJ)+PBLOWSNW_DEP(JJ))
ELSE
!
IF ( OSNOWDRIFT ) THEN
@@ -4024,7 +4064,7 @@ DO JJ=1,SIZE(PSNOW(:)) ! grid point loop
PSNOWGRAN2(JJ,1),PSNOWGRAN2F(JJ),HSNOWMETAMO )
!
IF ( ( ZDIFTYPE_SUP<XDIFF_1 .AND. PSNOWDZ(JJ,1)< ZDZOPT(JJ,1) ) .OR. &
- ( PSR(JJ)<XSNOWFALL_THRESHOLD .AND. PSNOWDZ(JJ,1)<2.*ZDZOPT(JJ,1) ) .OR. &
+ ( (PSR(JJ)+PBLOWSNW_DEP(JJ))<XSNOWFALL_THRESHOLD .AND. PSNOWDZ(JJ,1)<2.*ZDZOPT(JJ,1) ) .OR. &
PSNOWDZ(JJ,1)<XDZMIN_TOP_EXTREM ) THEN
!
! Fresh snow is similar to a shallow surface layer (< ZDZOPT)
diff --git a/src/SURFEX/write_diag_isban.F90 b/src/SURFEX/write_diag_isban.F90
index e389c8e79..12f99458a 100644
--- a/src/SURFEX/write_diag_isban.F90
+++ b/src/SURFEX/write_diag_isban.F90
@@ -3,7 +3,7 @@
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!SFX_LIC for details. version 1.
! #########
-SUBROUTINE WRITE_DIAG_ISBA_n (DTCO, DUO, U, IM, NDST, HPROGRAM, HWRITE)
+SUBROUTINE WRITE_DIAG_ISBA_n (DTCO, DUO, U, IM, NDST,BLOWSNW, HPROGRAM, HWRITE)
! ###############################################################################
!
!!**** *WRITE_DIAG_ISBA_n * - Stores ISBA diagnostics
@@ -34,6 +34,7 @@ USE MODD_DIAG_n, ONLY : DIAG_OPTIONS_t
USE MODD_SURF_ATM_n, ONLY : SURF_ATM_t
USE MODD_SURFEX_n, ONLY : ISBA_MODEL_t
USE MODD_DST_n, ONLY : DST_NP_t
+USE MODD_BLOWSNW_n, ONLY : BLOWSNW_t
!
USE MODD_SURF_PAR, ONLY : XUNDEF
!
@@ -55,6 +56,7 @@ TYPE(DIAG_OPTIONS_t), INTENT(INOUT) :: DUO
TYPE(SURF_ATM_t), INTENT(INOUT) :: U
TYPE(ISBA_MODEL_t), INTENT(INOUT) :: IM
TYPE(DST_NP_t), INTENT(INOUT) :: NDST
+TYPE(BLOWSNW_t), INTENT(INOUT) :: BLOWSNW
!
CHARACTER(LEN=6), INTENT(IN) :: HPROGRAM ! program calling surf. schemes
CHARACTER(LEN=3), INTENT(IN) :: HWRITE ! 'PGD' : only physiographic fields are written
@@ -69,7 +71,7 @@ IF (LHOOK) CALL DR_HOOK('WRITE_DIAG_ISBA_N',0,ZHOOK_HANDLE)
IF (HWRITE/='PGD') THEN
IF (IM%ID%O%XDIAG_TSTEP==XUNDEF .OR. &
ABS(NINT(IM%S%TTIME%TIME/IM%ID%O%XDIAG_TSTEP)*IM%ID%O%XDIAG_TSTEP-IM%S%TTIME%TIME)<1.E-3 ) THEN
- CALL WRITE_DIAG_SEB_ISBA_n(DTCO, DUO, U, IM%NCHI, IM%CHI, IM%ID, NDST, IM%GB, &
+ CALL WRITE_DIAG_SEB_ISBA_n(DTCO, DUO, U, IM%NCHI, IM%CHI, IM%ID, NDST, BLOWSNW, IM%GB, &
IM%O, IM%S, IM%NP, IM%NPE, HPROGRAM)
CALL WRITE_DIAG_MISC_ISBA_n(DTCO, DUO%CSELECT, DUO%LSNOWDIMNC, U, IM%ID%O%LPATCH_BUDGET, &
IM%ID%D, IM%ID%ND, IM%ID%DM, IM%ID%NDM, IM%O, IM%S, IM%K, &
diff --git a/src/SURFEX/write_diag_naturen.F90 b/src/SURFEX/write_diag_naturen.F90
index 0294d464f..85e85ff41 100644
--- a/src/SURFEX/write_diag_naturen.F90
+++ b/src/SURFEX/write_diag_naturen.F90
@@ -3,7 +3,7 @@
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!SFX_LIC for details. version 1.
! #########
-SUBROUTINE WRITE_DIAG_NATURE_n (DTCO, DUO, U, IM, NDST, HPROGRAM,HWRITE)
+SUBROUTINE WRITE_DIAG_NATURE_n (DTCO, DUO, U, IM, NDST, BLOWSNW, HPROGRAM,HWRITE)
! ###############################################################################
!
!!**** *WRITE_DIAG_NATURE_n * - Chooses the surface schemes for diagnostics over
@@ -34,6 +34,7 @@ USE MODD_DIAG_n, ONLY : DIAG_OPTIONS_t
USE MODD_SURF_ATM_n, ONLY : SURF_ATM_t
USE MODD_SURFEX_n, ONLY : ISBA_MODEL_t
USE MODD_DST_n, ONLY : DST_NP_t
+USE MODD_BLOWSNW_n, ONLY : BLOWSNW_t
!
USE MODD_SURF_PAR, ONLY : XUNDEF
!
@@ -53,6 +54,7 @@ TYPE(DIAG_OPTIONS_t), INTENT(INOUT) :: DUO
TYPE(SURF_ATM_t), INTENT(INOUT) :: U
TYPE(ISBA_MODEL_t), INTENT(INOUT) :: IM
TYPE(DST_NP_t), INTENT(INOUT) :: NDST
+TYPE(BLOWSNW_t), INTENT(INOUT) :: BLOWSNW
!
CHARACTER(LEN=6), INTENT(IN) :: HPROGRAM ! program calling surf. schemes
CHARACTER(LEN=3), INTENT(IN) :: HWRITE ! 'PGD' : only physiographic fields are written
@@ -66,7 +68,7 @@ REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
IF (LHOOK) CALL DR_HOOK('WRITE_DIAG_NATURE_N',0,ZHOOK_HANDLE)
IF (U%CNATURE=='ISBA ') THEN
- CALL WRITE_DIAG_ISBA_n(DTCO, DUO, U, IM, NDST, HPROGRAM,HWRITE)
+ CALL WRITE_DIAG_ISBA_n(DTCO, DUO, U, IM, NDST, BLOWSNW, HPROGRAM,HWRITE)
END IF
IF (LHOOK) CALL DR_HOOK('WRITE_DIAG_NATURE_N',1,ZHOOK_HANDLE)
!
diff --git a/src/SURFEX/write_diag_seb_isban.F90 b/src/SURFEX/write_diag_seb_isban.F90
index 11f7ee406..578764bfa 100644
--- a/src/SURFEX/write_diag_seb_isban.F90
+++ b/src/SURFEX/write_diag_seb_isban.F90
@@ -3,7 +3,7 @@
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!SFX_LIC for details. version 1.
! #########
- SUBROUTINE WRITE_DIAG_SEB_ISBA_n ( DTCO, DUO, U, NCHI, CHI, ID, NDST, GB, &
+ SUBROUTINE WRITE_DIAG_SEB_ISBA_n ( DTCO, DUO, U, NCHI, CHI, ID, NDST, BLOWSNW, GB, &
IO, S, NP, NPE, HPROGRAM)
! #################################
!
@@ -56,6 +56,7 @@ USE MODD_DST_n, ONLY : DST_NP_t
USE MODD_GR_BIOG_n, ONLY : GR_BIOG_t
USE MODD_ISBA_OPTIONS_n, ONLY : ISBA_OPTIONS_t
USE MODD_ISBA_n, ONLY : ISBA_NP_t, ISBA_P_t, ISBA_NPE_t, ISBA_PE_t, ISBA_S_t
+USE MODD_BLOWSNW_n, ONLY : BLOWSNW_t
!
#ifdef SFX_ARO
USE MODD_IO_SURF_ARO, ONLY : NBLOCK
@@ -72,6 +73,7 @@ USE MODD_CSTS, ONLY : XRHOLW, XTT, XLMTT
USE MODD_DST_SURF
!
USE MODD_AGRI, ONLY : LAGRIP
+USE MODD_BLOWSNW_SURF, ONLY : LBLOWSNW_CANODIAG
!
USE MODE_DIAG
!
@@ -106,6 +108,7 @@ TYPE(ISBA_OPTIONS_t), INTENT(INOUT) :: IO
TYPE(ISBA_S_t), INTENT(INOUT) :: S
TYPE(ISBA_NP_t), INTENT(INOUT) :: NP
TYPE(ISBA_NPE_t), INTENT(INOUT) :: NPE
+TYPE(BLOWSNW_t), INTENT(INOUT) :: BLOWSNW
!
CHARACTER(LEN=6), INTENT(IN) :: HPROGRAM ! program calling
!
@@ -121,7 +124,7 @@ CHARACTER(LEN=100):: YCOMMENT ! Comment string
CHARACTER(LEN=2) :: YNUM
!
LOGICAL :: GRESET
-INTEGER :: JSV, JSW, JP, ISIZE
+INTEGER :: JSV, JSW, JP, ISIZE, JLAYER
INTEGER :: ISIZE_LMEB_PATCH ! Number of patches where multi-energy balance should be applied
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
@@ -671,6 +674,93 @@ IF (CHI%SVI%NDSTEQ > 0)THEN
!
ENDIF
!
+!
+! Blowing snow variables
+!
+IF (CHI%SVI%NSNWEQ > 0)THEN
+
+ YRECFM='SNOW_SALT'
+ YCOMMENT='streamwise snow saltation flux (kg/m/s)'
+ CALL WRITE_SURF(DUO%CSELECT, &
+ HPROGRAM,YRECFM,BLOWSNW%XSNW_FSALT(:,1),IRESP,HCOMMENT=YCOMMENT)
+
+DO JSV=1,2
+
+ WRITE(YRECFM,'(A8,I1.1,A3)') 'SNW_FTUR',JSV,' '
+ YCOMMENT='Ins. surface turbulent snow flux (__ /m2/s)'
+ CALL WRITE_SURF(DUO%CSELECT, &
+ HPROGRAM,YRECFM,BLOWSNW%XSNW_FTURB(:,JSV),IRESP,HCOMMENT=YCOMMENT)
+
+ WRITE(YRECFM,'(A8,I1.1,A3)') 'SNW_FSED',JSV,' '
+ YCOMMENT='Ins. surface sedimentation snow flux (__ /m2/s)'
+ CALL WRITE_SURF(DUO%CSELECT, &
+ HPROGRAM,YRECFM,BLOWSNW%XSNW_FSED(:,JSV),IRESP,HCOMMENT=YCOMMENT)
+
+ WRITE(YRECFM,'(A8,I1.1,A3)') 'SNW_FNET',JSV,' '
+ YCOMMENT='Ins. surface net snow flux (__ /m2/s)'
+ CALL WRITE_SURF(DUO%CSELECT, &
+ HPROGRAM,YRECFM,BLOWSNW%XSNW_FNET(:,JSV),IRESP,HCOMMENT=YCOMMENT)
+
+ENDDO
+
+ YRECFM='SNW_FTUR_ACC'
+ YCOMMENT='Acc. surface turbulent snow flux (kg/m2)'
+ CALL WRITE_SURF(DUO%CSELECT, &
+ HPROGRAM,YRECFM,BLOWSNW%XSNW_FTURB(:,3),IRESP,HCOMMENT=YCOMMENT)
+
+ YRECFM='SNW_FSED_ACC'
+ YCOMMENT='Acc. surface sedimentation snow flux (kg/m2)'
+ CALL WRITE_SURF(DUO%CSELECT, &
+ HPROGRAM,YRECFM,BLOWSNW%XSNW_FSED(:,3),IRESP,HCOMMENT=YCOMMENT)
+
+ YRECFM='SNW_FNET_ACC'
+ YCOMMENT='Acc. surface net snow flux (kg/m2)'
+ CALL WRITE_SURF(DUO%CSELECT, &
+ HPROGRAM,YRECFM,BLOWSNW%XSNW_FNET(:,3),IRESP,HCOMMENT=YCOMMENT)
+
+ YRECFM='SNW_FSAL_ACC'
+ YCOMMENT='Acc. surface saltation flux (kg/m2)'
+ CALL WRITE_SURF(DUO%CSELECT, &
+ HPROGRAM,YRECFM,BLOWSNW%XSNW_FSALT(:,3),IRESP,HCOMMENT=YCOMMENT)
+
+ YRECFM='SNW_FSAL_INS'
+ YCOMMENT='Ins. surface saltatio flux (kg/m2)'
+ CALL WRITE_SURF(DUO%CSELECT, &
+ HPROGRAM,YRECFM,BLOWSNW%XSNW_FSALT(:,2),IRESP,HCOMMENT=YCOMMENT)
+
+ YRECFM='SNW_SUBL_ACC'
+ YCOMMENT='Canopy Acc. sublimation (kg/m2)'
+ CALL WRITE_SURF(DUO%CSELECT, &
+ HPROGRAM,YRECFM,BLOWSNW%XSNW_SUBL(:,3),IRESP,HCOMMENT=YCOMMENT)
+
+ YRECFM='SNW_SUBL_INS'
+ YCOMMENT='Canopy Sublimation Rate (mmSWE/day)'
+ CALL WRITE_SURF(DUO%CSELECT,&
+ HPROGRAM,YRECFM,BLOWSNW%XSNW_SUBL(:,2),IRESP,HCOMMENT=YCOMMENT)
+
+
+ IF(LBLOWSNW_CANODIAG) THEN
+ DO JLAYER=1,SIZE(BLOWSNW%XSNW_CANO_RGA,2)
+
+ WRITE(YRECFM,'(A10,I2.2)') 'CANSNW_RGL',JLAYER
+ YCOMMENT='Blowing snow radius at canopy level (m)'
+ CALL WRITE_SURF(DUO%CSELECT, &
+ HPROGRAM,YRECFM,BLOWSNW%XSNW_CANO_RGA(:,JLAYER),IRESP,HCOMMENT=YCOMMENT)
+
+ WRITE(YRECFM,'(A10,I2.2)') 'CANSNW_MAS',JLAYER
+ YCOMMENT='Blowing snow mass at canopy level (kg/m3)'
+ CALL WRITE_SURF(DUO%CSELECT, &
+ HPROGRAM,YRECFM,BLOWSNW%XSNW_CANO_VAR(:,JLAYER,2),IRESP,HCOMMENT=YCOMMENT)
+
+ WRITE(YRECFM,'(A10,I2.2)') 'CANSNW_NUM',JLAYER
+ YCOMMENT='Blowing snow number at canopy level (#/m3)'
+ CALL WRITE_SURF(DUO%CSELECT, &
+ HPROGRAM,YRECFM,BLOWSNW%XSNW_CANO_VAR(:,JLAYER,1),IRESP,HCOMMENT=YCOMMENT)
+ ENDDO
+ ENDIF
+
+ENDIF
+
!----------------------------------------------------------------------------
!
!* 5. Cumulated Energy fluxes
diff --git a/src/SURFEX/write_diag_surf_atmn.F90 b/src/SURFEX/write_diag_surf_atmn.F90
index e24dbc09e..9494df2cf 100644
--- a/src/SURFEX/write_diag_surf_atmn.F90
+++ b/src/SURFEX/write_diag_surf_atmn.F90
@@ -79,7 +79,7 @@ IF (YSC%U%NDIM_SEA >0) CALL WRITE_DIAG_SEA_n(YSC%DTCO, YSC%DUO, YSC%U, YSC%SM
IF (YSC%U%NDIM_WATER >0) CALL WRITE_DIAG_INLAND_WATER_n(YSC%DTCO, YSC%DUO, YSC%U, &
YSC%WM, YSC%FM, HPROGRAM,HWRITE)
IF (YSC%U%NDIM_NATURE >0) CALL WRITE_DIAG_NATURE_n(YSC%DTCO, YSC%DUO, YSC%U, YSC%IM, &
- YSC%NDST, HPROGRAM,HWRITE)
+ YSC%NDST,YSC%BLOWSNW, HPROGRAM,HWRITE)
IF (YSC%U%NDIM_TOWN >0) CALL WRITE_DIAG_TOWN_n(YSC%DTCO, YSC%DUO%CSELECT, YSC%U, YSC%TM, &
YSC%GDM, YSC%GRM, HPROGRAM,HWRITE)
!
diff --git a/src/SURFEX/write_isban.F90 b/src/SURFEX/write_isban.F90
index 914b3d949..0c816bedb 100644
--- a/src/SURFEX/write_isban.F90
+++ b/src/SURFEX/write_isban.F90
@@ -95,7 +95,7 @@ IF (LHOOK) CALL DR_HOOK('WRITE_ISBA_N',0,ZHOOK_HANDLE)
U%NSIZE_NATURE, HPROGRAM,OLAND_USE)
!
IF ((.NOT.LNOWRITE_CANOPY).OR.SIZE(HSELECT)>0) THEN
- CALL WRITESURF_SBL_n(HSELECT, IM%O%LCANOPY, IM%SB, HPROGRAM, HWRITE, "NATURE")
+ CALL WRITESURF_SBL_n(HSELECT, IM%O%LCANOPY, IM%SB, HPROGRAM, HWRITE, "NATURE",SV=IM%CHI%SVI)
ENDIF
!
CALL END_IO_SURF_n(HPROGRAM)
diff --git a/src/SURFEX/writesurf_sbln.F90 b/src/SURFEX/writesurf_sbln.F90
index cc44a36bd..131ef78c0 100644
--- a/src/SURFEX/writesurf_sbln.F90
+++ b/src/SURFEX/writesurf_sbln.F90
@@ -3,7 +3,7 @@
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!SFX_LIC for details. version 1.
! #########
- SUBROUTINE WRITESURF_SBL_n (HSELECT, OSBL, SB, HPROGRAM, HWRITE, HSURF)
+ SUBROUTINE WRITESURF_SBL_n (HSELECT, OSBL, SB, HPROGRAM, HWRITE, HSURF,SV)
! ####################################
!
!!**** *WRITE_FLAKE_n* - writes FLAKE fields
@@ -39,6 +39,7 @@
! ------------
!
USE MODD_CANOPY_n, ONLY : CANOPY_t
+USE MODD_SV_n, ONLY : SV_t
!
USE MODI_WRITE_SURF
USE MODI_END_IO_SURF_n
@@ -56,6 +57,7 @@ IMPLICIT NONE
LOGICAL, INTENT(IN) :: OSBL
!
TYPE(CANOPY_t), INTENT(INOUT) :: SB
+TYPE(SV_t), INTENT(IN),OPTIONAL :: SV
!
CHARACTER(LEN=6), INTENT(IN) :: HPROGRAM ! program calling
CHARACTER(LEN=3), INTENT(IN) :: HWRITE ! 'PREP' : does not write SBL XUNDEF fields
@@ -71,7 +73,7 @@ INTEGER :: IRESP ! IRESP : return-code if a problem appears
CHARACTER(LEN=13) :: YFORMAT
CHARACTER(LEN=100):: YCOMMENT ! Comment string
!
-INTEGER :: JL ! loop counter on layers
+INTEGER :: JL,JN ! loop counter on layers
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!-------------------------------------------------------------------------------
!
@@ -204,6 +206,20 @@ IF (HWRITE/='PRE') THEN
CALL WRITE_SURF(HSELECT,HPROGRAM,YRECFM,SB%XP(:,JL),IRESP,HCOMMENT=YCOMMENT)
END DO
!
+ IF (HSURF=="NATURE") THEN
+ IF(SV%NSNWEQ>0) THEN
+ DO JN=1,SV%NSNWEQ
+ !DO JN=1,2
+ DO JL=1,SB%NLVL
+ WRITE(YRECFM,'(A8,I1.1,A1,I2.2)') 'CANSNW_M',JN,'L',JL
+ YCOMMENT='Blown snow variables at canopy levels (__ /kg)'
+ CALL WRITE_SURF(HSELECT,HPROGRAM,YRECFM,SB%XBLOWSNW(:,JL,JN),IRESP,HCOMMENT=YCOMMENT)
+ END DO
+ END DO
+ ENDIF
+ ENDIF
+
+
ENDIF
!
IF (LHOOK) CALL DR_HOOK('WRITESURF_SBL_N',1,ZHOOK_HANDLE)
--
GitLab