From 3857e9ff4edbf06d109b4ccfb57a96e8c3d2f39e Mon Sep 17 00:00:00 2001
From: Gaelle DELAUTIER <gaelle.delautier@meteo.fr>
Date: Tue, 19 Jun 2018 14:44:07 +0200
Subject: [PATCH] Gaellle 19/6/2018 : mv *.f90 *.F90
---
src/SURFEX/blowsnw_dep.f90 | 66 --
src/SURFEX/blowsnw_diffk.f90 | 69 --
src/SURFEX/blowsnw_init_names.f90 | 97 --
src/SURFEX/blowsnw_velgrav1d.f90 | 286 -----
src/SURFEX/canopy_blowsnw_subl.f90 | 276 -----
src/SURFEX/canopy_evol_blowsnw.f90 | 274 -----
src/SURFEX/gamma.f90 | 224 ----
src/SURFEX/gamma_inc_low.f90 | 130 ---
src/SURFEX/modd_blowsnw_sedim_lkt1d.f90 | 54 -
src/SURFEX/modd_blowsnw_surf.f90 | 77 --
src/SURFEX/modd_blowsnwn.f90 | 51 -
src/SURFEX/mode_blowsnw_mblutl.f90 | 667 ------------
src/SURFEX/mode_blowsnw_sedim_lkt1d.f90 | 1297 -----------------------
src/SURFEX/mode_blowsnw_surf.f90 | 118 ---
src/SURFEX/modn_blowsnw.f90 | 42 -
src/SURFEX/snowpack_evol.f90 | 744 -------------
16 files changed, 4472 deletions(-)
delete mode 100644 src/SURFEX/blowsnw_dep.f90
delete mode 100644 src/SURFEX/blowsnw_diffk.f90
delete mode 100644 src/SURFEX/blowsnw_init_names.f90
delete mode 100644 src/SURFEX/blowsnw_velgrav1d.f90
delete mode 100644 src/SURFEX/canopy_blowsnw_subl.f90
delete mode 100644 src/SURFEX/canopy_evol_blowsnw.f90
delete mode 100644 src/SURFEX/gamma.f90
delete mode 100644 src/SURFEX/gamma_inc_low.f90
delete mode 100644 src/SURFEX/modd_blowsnw_sedim_lkt1d.f90
delete mode 100644 src/SURFEX/modd_blowsnw_surf.f90
delete mode 100644 src/SURFEX/modd_blowsnwn.f90
delete mode 100644 src/SURFEX/mode_blowsnw_mblutl.f90
delete mode 100644 src/SURFEX/mode_blowsnw_sedim_lkt1d.f90
delete mode 100644 src/SURFEX/mode_blowsnw_surf.f90
delete mode 100644 src/SURFEX/modn_blowsnw.f90
delete mode 100644 src/SURFEX/snowpack_evol.f90
diff --git a/src/SURFEX/blowsnw_dep.f90 b/src/SURFEX/blowsnw_dep.f90
deleted file mode 100644
index 0d82d32e8..000000000
--- a/src/SURFEX/blowsnw_dep.f90
+++ /dev/null
@@ -1,66 +0,0 @@
-
-!###########################################################
- SUBROUTINE BLOWSNW_DEP (K2D_SNWBEG, K2D_SNWEND,PTA,PPA, &
- PRHODREF,PSVT,PSFSNW)
-!###########################################################
- !-------------------------------------------------------------------------------
- !
- !* 0. DECLARATIONS
- ! ------------
- !
- USE MODE_BLOWSNW_SURF
- USE MODI_BLOWSNW_VELGRAV1D
- USE MODD_CSTS, ONLY : XG,XRHOLI, XPI
- USE MODD_BLOWSNW_SURF
-
- !
- IMPLICIT NONE
- !
- !* 0.1 Declarations of dummy arguments :
- !
- INTEGER, INTENT(IN) :: K2D_SNWBEG, K2D_SNWEND ! index of first and last blowing snow 2D variable sent to MNH
- REAL, DIMENSION(:), INTENT(IN) :: PTA ! air temperature
- REAL, DIMENSION(:), INTENT(IN) :: PPA ! air pressure
- REAL, DIMENSION(:), INTENT(IN) :: PRHODREF ! air density
- REAL, DIMENSION(:,:), INTENT(INOUT) :: PSVT ! blowing snow concentration
- ! in surface units (__/m3)
- REAL, DIMENSION(:,:), INTENT(INOUT) :: PSFSNW ! blowing snow concentration fluxes
-
-
-
-
- !
- !
- !* 0.2 Declarations of local variables :
- !
- REAL, DIMENSION(SIZE(PSVT,1),1, SIZE(PSVT,2)) :: ZSVT
- REAL, DIMENSION(SIZE(PSVT,1),1) :: ZBET, ZRG,ZTA,ZPA
- REAL, DIMENSION(SIZE(PSVT,1),1,SIZE(PSVT,2)) :: ZVGK ! Terminal fallspeed (m/s)
-
- INTEGER :: JN
-
-
-! Save scalars in local array
-DO JN=1,SIZE(PSVT,2)
- ZSVT(:,1,JN)=MAX(PSVT(:,JN),1E-60)
-END DO
-ZTA(:,1)=PTA(:)
-ZPA(:,1)=PPA(:)
-
-! Get gamma parameter distribution : scale factor and mean radius at the first
-! level in the atmosphere
-CALL SNOWMOMENT2SIZE(ZSVT, PBETA1D=ZBET, PRG1D=ZRG )
-! Compute mass-weighted terminal fall speed based on particles distribution and
-! atmospheric conditions.
-CALL BLOWSNW_VELGRAV1D(ZBET, ZRG, ZTA, PRHODREF,ZPA, ZVGK)
-
-! Compute sedimentation fluxes of blowing snow variables
-
-PSVT(:,K2D_SNWBEG) = PSVT(:,1) * ZVGK(:,1,1)
-PSVT(:,K2D_SNWBEG+1) = PSVT(:,2) * ZVGK(:,1,2)
-
-
-PSFSNW(:,1:2) = PSVT(:,K2D_SNWBEG:(K2D_SNWBEG+1))
-
-
-END SUBROUTINE BLOWSNW_DEP
diff --git a/src/SURFEX/blowsnw_diffk.f90 b/src/SURFEX/blowsnw_diffk.f90
deleted file mode 100644
index 725d3271f..000000000
--- a/src/SURFEX/blowsnw_diffk.f90
+++ /dev/null
@@ -1,69 +0,0 @@
-!! #######################################
-SUBROUTINE BLOWSNW_DIFFK(PK, PTKE, PVGK,KSNW,PKSNW)
-
-!! Compute particle eddy diffusivity for number and momentum
-!! ZZETA is the ratio between eddy diffusivity for momentum
-!! znd eddy diffusivity for number and mass.
-!! The value 0.25 is based on a sensititivy analysis using
-!! vertical profile of blowing snow fluxes and concentration
-!! measured using SPC and mechanical snow trap at Col du Lac
-!! Blanc experimental site. More in details in the PhD thesis
-!! of V. Vionnet : Etudes du transport de la neige par le vent en
-!! conditions alpines : observations et modelisation a l'aide d'un
-!! modele couple atmosphere/manteau neigeux
-
-USE MODD_BLOWSNW_SURF
-!
-USE YOMHOOK ,ONLY : LHOOK, DR_HOOK
-USE PARKIND1 ,ONLY : JPRB
-
-IMPLICIT NONE
- !
- !* 0.1 Declarations of dummy arguments :
- !
-REAL, DIMENSION(:,:), INTENT(IN) :: PK ! flow eddy diffusivity
-REAL, DIMENSION(:,:), INTENT(IN) :: PTKE ! turbulent kineti energy
-REAL, DIMENSION(:,:,:), INTENT(IN) :: PVGK ! particle fall speed
-INTEGER, INTENT(IN) :: KSNW ! number of blowing snow variables
-!
-REAL, DIMENSION(:,:,:), INTENT(OUT) :: PKSNW ! particle eddy diffusivity
- ! including inertial effects
- !
- !
- !* 0.2 Declarations of local variables :
- !
-REAL :: ZBETAC
-REAL :: ZZETA
-INTEGER :: JSV ! Loop counter
-LOGICAL :: ONEW_K
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!
- !
- !
- !* 1.0 Initialize constants
- !
-IF (LHOOK) CALL DR_HOOK('BLOWSNW_DIFFK',0,ZHOOK_HANDLE)
-
-
-ZBETAC = 1.
-ONEW_K=.FALSE.
-
- !
- !
- !* 2.0 Compute eddy diffusivity for number and mass
- !
-IF(ONEW_K) THEN
- DO JSV =1,KSNW
- PKSNW(:,:,JSV) = PK(:,:)*(1.+(ZBETAC*PVGK(:,:,JSV))**2./(1./3.*PTKE(:,:)))**(-1.)
- ENDDO
-ELSE
-!
- ZZETA = 0.25
- DO JSV =1,KSNW
- PKSNW(:,:,JSV)=PK(:,:)/XRSNOW_SBL
- ENDDO
-ENDIF
-
-IF (LHOOK) CALL DR_HOOK('BLOWSNW_DIFFK',1,ZHOOK_HANDLE)
-
-END SUBROUTINE BLOWSNW_DIFFK
diff --git a/src/SURFEX/blowsnw_init_names.f90 b/src/SURFEX/blowsnw_init_names.f90
deleted file mode 100644
index 4d6d8746e..000000000
--- a/src/SURFEX/blowsnw_init_names.f90
+++ /dev/null
@@ -1,97 +0,0 @@
-!! ###########################################
-SUBROUTINE BLOWSNW_INIT_NAMES (KLUOUT,HSV,KSNWEQ,KSV_SNWBEG,KSV_SNWEND, &
- K2DSNWEQ,K2D_SNWBEG,K2D_SNWEND )
-!! ###########################################
-!!
-!!*** *SNW_INIT_NAMES*
-!!
-!! PURPOSE
-!! -------
-!! Read and filter all chemical species into the CSV array
-!! initialize NSV_SNWBEG and NSV_SNWEND index for the begin and the ending chemical index
-!!
-!!
-!! REFERENCE
-!! ---------
-!! Modified dst_init_names (february 2005)
-!!
-!! AUTHOR
-!! ------
-!! Vincent VIONNET <vincent.vionnet@meteo.fr>
-!!
-!! MODIFICATIONS
-!! -------------
-!!
-!! EXTERNAL
-!! --------
-!!
-!! IMPLICIT ARGUMENTS
-!! ------------------
-!-------------------------------------------------------------------------------
-!
-!* 0. DECLARATIONS
-! ------------
-USE YOMHOOK ,ONLY : LHOOK, DR_HOOK
-USE PARKIND1 ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!* 0.1 declarations of arguments
-!
-
-INTEGER, INTENT(IN) :: KLUOUT ! output listing channel
-CHARACTER(LEN=*), DIMENSION(:), INTENT(IN) :: HSV ! name of chemical species
- ! with character # for chemistry
-INTEGER, INTENT(OUT) :: KSNWEQ ! number of blowing snow related variables
-INTEGER, INTENT(OUT) :: KSV_SNWBEG ! first blowing snow related scalar
-INTEGER, INTENT(OUT) :: KSV_SNWEND ! last blowing snow related scalar
-INTEGER, INTENT(OUT) :: K2DSNWEQ ! number of 2D blowing snow related variables
-INTEGER, INTENT(OUT) :: K2D_SNWBEG ! first 2D blowing snow related scalar
-INTEGER, INTENT(OUT) :: K2D_SNWEND ! last 2D blowing snow related scalar
-!
-!* 0.2 declarations of local variables
-INTEGER :: JSV !! loop on scalar variables
-CHARACTER(LEN=4) :: YRC1
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!
-!-------------------------------------------------------------------------------
-
-!Initialize output variables
-IF (LHOOK) CALL DR_HOOK('BLOWSNW_INIT_NAMES',0,ZHOOK_HANDLE)
-
-KSNWEQ = 0
-KSV_SNWBEG=0
-KSV_SNWEND=0
-
-K2DSNWEQ = 0
-K2D_SNWBEG=0
-K2D_SNWEND=0
-
-DO JSV=1, SIZE(HSV)
- YRC1= HSV(JSV)(1:4)
- IF (YRC1 == 'SNWM') THEN
- KSNWEQ = KSNWEQ + 1
- IF (KSNWEQ == 1) THEN
- KSV_SNWBEG=JSV
- ENDIF !Check on first time
- ELSE IF (YRC1 == 'SNWC') THEN
- K2DSNWEQ = K2DSNWEQ + 1
- IF (K2DSNWEQ == 1) THEN
- K2D_SNWBEG=JSV
- ENDIF !Check on first time
- ELSE !Not snow variables
- !DO NOTHING
- ENDIF
-ENDDO
-!
-! Set the output list of scalar to the input list of scalars
-
-! Get the index of the last blowing snow relevant tracer
-KSV_SNWEND = KSV_SNWBEG + KSNWEQ -1
-! Get the index of the last 2D blowing snow relevant tracer
-K2D_SNWEND = K2D_SNWBEG + K2DSNWEQ -1
-
-!
-IF (LHOOK) CALL DR_HOOK('BLOWSNW_INIT_NAMES',1,ZHOOK_HANDLE)
-!
-END SUBROUTINE BLOWSNW_INIT_NAMES
diff --git a/src/SURFEX/blowsnw_velgrav1d.f90 b/src/SURFEX/blowsnw_velgrav1d.f90
deleted file mode 100644
index d5bc2a19d..000000000
--- a/src/SURFEX/blowsnw_velgrav1d.f90
+++ /dev/null
@@ -1,286 +0,0 @@
-!SFX_LIC Copyright 1994-2018 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 BLOWSNW_VELGRAV1D(PBETA, PRG, PTA, PRHODREF,PPABST,PVGK)
-
-!! -------
-!!
-!! Compute number- and mass-averaged settling velocity for
-!! blowing snow particles based on several methods :
-!! - Mitchell (1996) : numerical integration assuming spherical
-!! particles (expensive)
-!! - Carrier (1953) and Dover (1993) : numerical integration (expensive)
-!! - look-up table based on Carrier (1953) depending on mean radius and
-!! pressure
-!! - None : assume no settling velocity
-!!
-!! REFERENCE
-!! ---------
-!!
-!! Mitchell (1996) : Use of mass- and area-dimensionla power laws for
-!! determining precipitation particle terminal velocities, JAS,
-!! 53(12),1710-1723
-!! Carrier, C. : On Slow Viscous Flow, Tech. rep., Office of Naval Research, Contract Nonr-653(00), Brown
-!! University, Providence, RI, 1953.
-!! Dover, S. : Numerical Modelling of Blowing Snow, Ph.D. thesis, University of Leeds, U.K., 1993.
-!!
-!! AUTHOR
-!! ------
-!! V. Vionnet (CNRM/GMME/MOSAYC)
-!!
-!! MODIFICATIONS
-!! -------------
-!! Philippe Wautelet 28/05/2018: corrected truncated integer division (1*10**(-6) -> 1E-6)
-!!
-!! NB : this routine is similar to the routine implemented in Meso-NH (blowsnow_velgrav.f90)
-!
-!!
-
-USE MODD_BLOWSNW_SURF
-USE MODD_CSTS
-
-USE MODI_GAMMA_INC_LOW
-
-USE MODE_BLOWSNW_SEDIM_LKT1D
-
-USE YOMHOOK ,ONLY : LHOOK, DR_HOOK
-USE PARKIND1 ,ONLY : JPRB
-
-IMPLICIT NONE
- !
- !* 0.1 Declarations of dummy arguments :
- !
- REAL, DIMENSION(:,:), INTENT(IN) :: PBETA ! distribution scale parameter [m]
- REAL, DIMENSION(:,:), INTENT(IN) :: PRG ! mean radius [m]
- REAL, DIMENSION(:,:), INTENT(IN) :: PTA ! air temperature
- REAL, DIMENSION(:,:), INTENT(IN) :: PPABST ! air pressure
- REAL, DIMENSION(:), INTENT(IN) :: PRHODREF ! air density
- REAL, DIMENSION(:,:,:), INTENT(OUT) :: PVGK ! terminal fallspeed [m/s]
- !
- !
- !* 0.2 Declarations of local variables :
- !
-
- REAL, DIMENSION(SIZE(PTA,1),SIZE(PTA,2)) :: ZMU !air kinematic viscosity [m2/s]
- REAL, DIMENSION(SIZE(PTA,1),SIZE(PTA,2)) :: ZR1 ! first limit radius in integration
- REAL, DIMENSION(SIZE(PTA,1),SIZE(PTA,2)) :: ZR2 ! second limit radius in integration [m]
- REAL, DIMENSION(SIZE(PTA,1),SIZE(PTA,2)) :: ZAM1 ! Parameter avr in Mitchell
- REAL, DIMENSION(SIZE(PTA,1),SIZE(PTA,2)) :: ZAM2 !
- REAL, DIMENSION(SIZE(PTA,1),SIZE(PTA,2)) :: ZAM3 !
-
-REAL, DIMENSION(SIZE(PTA,1),SIZE(PTA,2)) :: ZAA
-REAL, DIMENSION(SIZE(PTA,1)) :: ZBB
-INTEGER, DIMENSION(SIZE(PTA,1),SIZE(PTA,2)) :: NMAX
-INTEGER :: ZNUM_EXP,ZMAS_EXP
-REAL :: ZGAM,ZVEL_CARRIER,ZR
-REAL :: ZW_M3,ZW_M0
-REAL :: ZSUM_VEL_M3,ZSUM_M3,ZSUM_VEL_M0,ZSUM_M0
-REAL :: ZDELTAR
-REAL :: ZGAMB,ZGAM_BM3, ZGAM_BM3B
- INTEGER JI,JJ,II,JK,ILAYER ! Loop counter
-
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-
- !
- !* 0.3 Initialization of variables :
- !
-IF (LHOOK) CALL DR_HOOK('BLOWSNW_VELGRAV1D',0,ZHOOK_HANDLE)
-
-ZDELTAR = 1e-6
-ZGAM = GAMMA(XEMIALPHA_SNW)
-ILAYER=SIZE(PTA,2)
-
-! Sutherland's equation for kinematic viscosity
-DO JK=1,ILAYER
- ZMU(:,JK)=1.8325d-5*416.16/(PTA(:,JK)+120)*(PTA(:,JK)/296.16)*SQRT(PTA(:,JK)/296.16)/PRHODREF(:)
-END DO
-
- !
- !* 1 Compute number- and mass-averaged settling velocity :
- !
-
-IF(CSNOW_SEDIM=='TABC') THEN
-!
-! Sedimentation of snow particles is computed according to Carrier's drag coofficient.
-! To reduce computational time; look-up tables are used. They depend on the
-! average radius and the pressure (interpolation)
-!
- CALL BLOWSNW_SEDIM_LKT1D(PRG,PPABST,PVGK)
-
-ELSE IF(CSNOW_SEDIM=='MITC') THEN ! Sedimentation following Mitchell (1996)
-
- ZGAMB = GAMMA(XEMIALPHA_SNW+3)
- ZGAM_BM3 = GAMMA(3*XBM3-1+XEMIALPHA_SNW)
- ZGAM_BM3B = GAMMA(3*XBM3+2+XEMIALPHA_SNW)
-
- ! Compute limit radius for integration of Mitchell's formulation
- ZR1(:,:)=RLIM(ZMU,PRHODREF,XBESTL_1)
- ZR2(:,:)=RLIM(ZMU,PRHODREF,XBESTL_2)
- ! Compute parameter avr for integration of Mitchell's formulation
- ZAM1(:,:)=AVR(XAM1,XBM1,PRHODREF,ZMU)
- ZAM2(:,:)=AVR(XAM2,XBM2,PRHODREF,ZMU)
- ZAM3(:,:)=AVR(XAM3,XBM3,PRHODREF,ZMU)
-
-DO JI=1,SIZE(PTA,1)
- DO JK=1,ILAYER
-!Number weighted terminal velocity
- PVGK(JI,JK,1)=(PBETA(JI,JK)**(3*XBM1-1)*ZAM1(JI,JK)* &
- GAMMA_INC_LOW(3*XBM1-1+XEMIALPHA_SNW,ZR1(JI,JK)/PBETA(JI,JK)) + &
- PBETA(JI,JK)**(3*XBM2-1)*ZAM2(JI,JK)* &
- (GAMMA_INC_LOW(3*XBM2-1+XEMIALPHA_SNW,ZR2(JI,JK)/PBETA(JI,JK))- &
- GAMMA_INC_LOW(3*XBM2-1+XEMIALPHA_SNW,ZR1(JI,JK)/PBETA(JI,JK)))+ &
- PBETA(JI,JK)**(3*XBM3-1)*ZAM3(JI,JK)* &
- (ZGAM_BM3-GAMMA_INC_LOW(3*XBM3-1+XEMIALPHA_SNW,ZR2(JI,JK)/PBETA(JI,JK))))/ZGAM
-!Mass weighted terminal velocity
- PVGK(JI,JK,2)=(PBETA(JI,JK)**(3*XBM1-1)*ZAM1(JI,JK)* &
- GAMMA_INC_LOW(3*XBM1+2+XEMIALPHA_SNW,ZR1(JI,JK)/PBETA(JI,JK)) + &
- PBETA(JI,JK)**(3*XBM2-1)*ZAM2(JI,JK)* &
- (GAMMA_INC_LOW(3*XBM2+2+XEMIALPHA_SNW,ZR2(JI,JK)/PBETA(JI,JK))- &
- GAMMA_INC_LOW(3*XBM2+2+XEMIALPHA_SNW,ZR1(JI,JK)/PBETA(JI,JK)))+ &
- PBETA(JI,JK)**(3*XBM3-1)*ZAM3(JI,JK)* &
- (ZGAM_BM3B-GAMMA_INC_LOW(3*XBM3+2+XEMIALPHA_SNW,ZR2(JI,JK)/PBETA(JI,JK))))/ZGAMB
- ENDDO
-ENDDO
-
-
-ELSE IF(CSNOW_SEDIM=='CARR') THEN
-! Settling velocity is computed according to Carrier's drag coofficient.
-! This method is used in other blowing snow model such as PIEKTUK or SNOWSTORM
-! We perfom a numerical integration since no analytical solution exists.
-
-ZAA(:,:) = 6.203*ZMU(:,:)/2
-ZBB(:) = 5.516*XRHOLI/(4*PRHODREF(:))*XG
-NMAX(:,:)=GET_INDEX(PBETA,ZDELTAR)
-
-! Exponent used to weight the number-averaged falling speed
-ZNUM_EXP=0.
-! Exponent used to weight the mass-averaged falling speed
-ZMAS_EXP=3.
-
-
-DO JI=1,SIZE(PTA,1)
- DO JK=1,ILAYER
- ZSUM_M3=0.
- ZSUM_M0=0.
- ZSUM_VEL_M0=0.
- ZSUM_VEL_M3=0.
- DO II=1,NMAX(JI,JK)
- ZR = 1E-6+(II-0.5)*ZDELTAR
- ZVEL_CARRIER = - ZAA(JI,JK)/ZR+((ZAA(JI,JK)/ZR)**2.+ZBB(JI)*ZR)**0.5
- ZW_M0=ZR**(XEMIALPHA_SNW-1)*exp(-ZR/PBETA(JI,JK))/(PBETA(JI,JK)**XEMIALPHA_SNW*ZGAM)
-
- ZW_M3=ZR**ZMAS_EXP*ZW_M0
- ZW_M0=ZR**ZNUM_EXP*ZW_M0
-
- ZSUM_M3 = ZSUM_M3+ZW_M3*ZDELTAR
- ZSUM_M0 = ZSUM_M0+ZW_M0*ZDELTAR
-
-
- ZSUM_VEL_M0 = ZSUM_VEL_M0+ ZW_M0*ZVEL_CARRIER*ZDELTAR
- ZSUM_VEL_M3 = ZSUM_VEL_M3+ ZW_M3*ZVEL_CARRIER*ZDELTAR
- ENDDO
- PVGK(JI,JK,1) = ZSUM_VEL_M0/ZSUM_M0
- PVGK(JI,JK,2) = ZSUM_VEL_M3/ZSUM_M3
- ENDDO
-ENDDO
-
-ELSE
-
-! Settling velocity is set equal to 0.
-PVGK(:,:,1) = 0.
-PVGK(:,:,2) = 0.
-
-
-END IF
-
-IF (LHOOK) CALL DR_HOOK('BLOWSNW_VELGRAV1D',1,ZHOOK_HANDLE)
-
-CONTAINS
-
-
-FUNCTION RLIM(PMU,PRHODREF,PBEST_LIM) RESULT(PRLIM)
-!
-!! PURPOSE
-!! -------
-! Calculate the radius of a sperical particle for a given Best Number
-!
-!
-USE MODD_CSTS, ONLY : XRHOLI,XG
-!
-!
-!* 0.1 declarations of arguments
-!
-REAL, DIMENSION(:), INTENT(IN) :: PRHODREF ! (kg/m3)
-REAL, DIMENSION(:,:), INTENT(IN) :: PMU ! (m2/s)
-REAL, INTENT(IN) :: PBEST_LIM! (-)
-
-!
-REAL, DIMENSION(SIZE(PMU,1),SIZE(PMU,2)) :: PRLIM ! (m)
-INTEGER JK
-!
-DO JK=1,SIZE(PMU,2)
- PRLIM(:,JK)=(3./32.*PRHODREF(:)/(XRHOLI*XG)*PMU(:,JK)**2.*PBEST_LIM)**0.333333333
-END DO
-
-END FUNCTION RLIM
-
-FUNCTION AVR(PARE,PBRE,PRHODEF,PMU) RESULT(PAVR)
-!
-!! PURPOSE
-!! -------
-! Calculate the parameter av_r in Michell formulation
-!
-!
-USE MODD_CSTS, ONLY : XRHOLI,XG
-
-!
-!
-!* 0.1 declarations of arguments
-!
-REAL, INTENT(IN) :: PARE ! (-)
-REAL, INTENT(IN) :: PBRE ! (-)
-REAL, DIMENSION(:), INTENT(IN) :: PRHODEF ! (kg/m3)
-REAL, DIMENSION(:,:), INTENT(IN) :: PMU ! (m2/s)
-
-!
-REAL, DIMENSION(SIZE(PMU,1),SIZE(PMU,2)) :: PAVR ! (-)
-!
-INTEGER JK
-
-DO JK=1,SIZE(PMU,2)
- PAVR(:,JK)=2.**(3.*PBRE-1.)*PARE*PMU(:,JK)**(1.-2.*PBRE)*(4./3.*XRHOLI/PRHODEF(:)*XG)**PBRE
-END DO
-
-END FUNCTION AVR
-
-FUNCTION GET_INDEX(PBETA,PDELTAR) RESULT(KMAX)
-!
-!! PURPOSE
-!! -------
-! Calculate the upper index in numerical integration of Carrier's formulation
-! Index equals to 5* mean radius
-!
-!
-USE MODD_BLOWSNW_SURF
-
-!
-IMPLICIT NONE
-!
-!* 0.1 declarations of arguments
-!
-REAL, INTENT(IN) :: PDELTAR ! (-)
-REAL, DIMENSION(:,:), INTENT(IN) :: PBETA ! (kg/m3)
-
-!
-INTEGER, DIMENSION(SIZE(PBETA,1),SIZE(PBETA,2)) :: KMAX ! (-)
-!
-
-KMAX(:,:)=int(PBETA(:,:)*XEMIALPHA_SNW*5/PDELTAR)
-
-
-END FUNCTION GET_INDEX
-
-END SUBROUTINE BLOWSNW_VELGRAV1D
diff --git a/src/SURFEX/canopy_blowsnw_subl.f90 b/src/SURFEX/canopy_blowsnw_subl.f90
deleted file mode 100644
index e6e2ad5eb..000000000
--- a/src/SURFEX/canopy_blowsnw_subl.f90
+++ /dev/null
@@ -1,276 +0,0 @@
-! #########################################
- SUBROUTINE CANOPY_BLOWSNW_SUBL(KI,KLVL,KSNW,PTSTEP,PBLOWSNW,PT,PQ,PRHOA,PP, &
- PBLOWSNWS,PDZ,PSFTQ,PSFTH,PSNW_SUBL )
-
-USE MODD_CSTS
-USE MODD_BLOWSNW_SURF
-!
-USE MODE_BLOWSNW_SURF
-!
-USE MODI_BLOWSNW_VELGRAV1D
-
-IMPLICIT NONE
-!
-!* 0.1 Declarations of arguments
-! -------------------------
-!
-INTEGER, INTENT(IN) :: KI ! number of horizontal points
-INTEGER, INTENT(IN) :: KLVL ! number of levels in canopy
-INTEGER, INTENT(IN) :: KSNW ! number of snow variables in canopy
-REAL, INTENT(IN) :: PTSTEP ! time-step
-REAL, DIMENSION(KI,KLVL), INTENT(IN) :: PT ! Temperature at canopy levels (K)
-REAL, DIMENSION(KI,KLVL),INTENT(IN) :: PP ! air presseure at canopy levels (Pa)
-REAL, DIMENSION(KI,KLVL), INTENT(IN) :: PQ ! humidity at canopy levels (kg/m3)
-REAL, DIMENSION(KI,KLVL,KSNW), INTENT(IN) :: PBLOWSNW ! blowing snow variables at canopy levels (1: #_s/kg, 2: kg_s/kg)
-REAL, DIMENSION(KI,KLVL), INTENT(IN) :: PDZ ! depth of canopy levels (m)
-
-REAL, DIMENSION(KI), INTENT(IN) :: PRHOA ! air density (kg/m3)
-
-REAL, DIMENSION(KI), INTENT(INOUT) :: PSFTH ! flux of heat (W/m2)
-REAL, DIMENSION(KI), INTENT(INOUT) :: PSFTQ ! flux of water vapor (kg/m2/s)
-
-REAL, DIMENSION(KI,KLVL,KSNW), INTENT(INOUT) :: PBLOWSNWS ! snow variables tendency due to sublimation at Canopy levels
-! (1: #_s/kg/s, 2: kg_s/kg/s)
-REAL, DIMENSION(KI,KSNW+1), INTENT(INOUT) :: PSNW_SUBL ! Diagnostic !Sublimation rate
- ! 1: Instantaneous number 2: Instantaneous mass 3: Accumulated Mass
-!
-!* 0.2 Declarations of local variables
-! -------------------------------
-!
-INTEGER JK,JI,JSV
-
-REAL, DIMENSION(KI,KLVL) :: ZBET ! Scale parameter of the gamma distribution (m)
-REAL, DIMENSION(KI,KLVL,KSNW) :: ZVGK ! Terminal fallspeed (m/s)
-REAL, DIMENSION(KI,KLVL) :: ZRG ! Mean radius (m)
-REAL, DIMENSION(KI,KLVL) :: ZEXN ! Exner function at full levels
-REAL, DIMENSION(KI,KLVL) :: ZW ! working array
-REAL, DIMENSION(KI) :: ZT ! average temprature of Canopy
-REAL, DIMENSION(KI,KLVL) :: ZLSFACT! L_s/(C_ph)
-
-REAL, DIMENSION(KI,KLVL,KSNW) :: ZSNW ! work variable : pot. temp at futur instant
-! ! (or past at the end of the routine)
-REAL, DIMENSION(KI,KLVL-1) :: ZMU ! air kinematic viscosity
-REAL, DIMENSION(KI,KLVL-1) :: ZNU ! Nusselt number (-)
-REAL, DIMENSION(KI,KLVL-1) :: ZKA ! air thermal conductivity
-REAL, DIMENSION(KI,KLVL-1) :: ZDV ! diffusivity of water vapor in the air.
-REAL, DIMENSION(KI,KLVL-1) :: ZAI ! denominator in Thorpe and Masson (66) equation
-REAL, DIMENSION(KI,KLVL-1) :: ZUSI ! undersaturation over ice
-REAL, DIMENSION(KI,KLVL-1) :: ZSUBL_NUM,ZSUBL_MASS ! sublimation rate (number and mass)
-REAL, DIMENSION(KI) :: ZCOL_SUBL ! integrated sublimation rate (kg/m2/s)
-!
-LOGICAL :: LSUBL_ALPINE3D
-REAL :: ZR_EFF ! Effective radius for the computation of sublimation
-REAL, DIMENSION(KI,KLVL) :: ZAA ! Coeff used in the computation of terminal fallspeed (m/s)
-REAL, DIMENSION(KI,KLVL) :: ZBB ! Coeff used in the computation of terminal fallspeed (m/s)
-!
-!-------------------------------------------------------------------------------
-!
-!
-!* 1. initializations
-!
-!-------------------------------------------------------------------------------
-DO JK=1,(KLVL-1)
-! Sutherland's equation for kinematic viscosity
- ZMU(:,JK)=1.8325d-5*416.16/(PT(:,JK)+120)*(PT(:,JK)/296.16)*SQRT(PT(:,JK)/296.16)/PRHOA(:)
-! Thermal conductivity of the air
- ZKA(:,JK) = 2.38E-2 + 0.0071E-2 * ( PT(:,JK) - XTT ) ! k_a
-! Diffusivity of water vapor in the air.
- ZDV(:,JK) = 0.211E-4 * (PT(:,JK)/XTT)**1.94 * (XP00/PP(:,JK)) ! D_v
-END DO
-
-!* Compute the denominator in the Thorpe and Masson (66) equation
-!
-ZAI(:,:) = EXP( XALPI - XBETAI/PT(:,1:(KLVL-1)) - XGAMI*ALOG(PT(:,1:(KLVL-1)))) ! es_i
-!
-DO JK=1,(KLVL-1)
- ! Undersaturation over ice
- ZUSI(:,JK) = PQ(:,JK)/PRHOA(:)*( PP(:,JK)-ZAI(:,JK) ) / ( (XMV/XMD) * ZAI(:,JK) ) - 1.0
- ! denominator in the Thorpe and Masson (66) equation
- ZAI(:,JK) = ( XLSTT + (XCPV-XCI)*(PT(:,JK)-XTT) ) / (ZKA(:,JK)*PT(:,JK)) &
- *( ( XLSTT + (XCPV-XCI)*(PT(:,JK)-XTT) ) / (XRV*PT(:,JK)) - 1.) &
- + (XRV*PT(:,JK)) / (ZDV(:,JK)*ZAI(:,JK))
-END DO
-
-ZUSI(:,:) =MIN(ZUSI(:,:),0.)
-
-!
-DO JK=1,(KLVL)
- ZW(:,JK) = XCPD+XCPV*PQ(:,JK)/PRHOA(:)
-END DO
-!
-ZLSFACT(:,:) = (XLSTT+(XCPV-XCI)*(PT(:,:)-XTT))/ZW(:,:) ! L_s/(*C_ph)
-!
-ZEXN = (PP/XP00)**(XRD/XCPD)
-!
-ZSUBL_NUM = 0.
-ZSUBL_MASS = 0.
-ZCOL_SUBL=0.
-!
-LSUBL_ALPINE3D = .FALSE. ! Compute sublimation using the method of reprsentative
- ! radius implemented in Alpine 3D (Groot and al, 2011)
-!
-
-IF(LSUBL_ALPINE3D) THEN
-!
-!-------------------------------------------------------------------------------
-!
-!
-!* 2. Compute settling velocity based on the effective radius
-! ---------------
-!
-!-------------------------------------------------------------------------------
- ZR_EFF = 73.5e-6 ! Effective radius computed following the Swiss
- ! method. This effective radius give the same total
- ! sublimation for a equal concentration an ensemble of
- ! gamma distributed particles with rm = 35e-6 m and
- ! alpha=3
- ZRG(:,:) = ZR_EFF
- ZBET(:,:) = ZR_EFF/XEMIALPHA_SNW
-!
-!
-!* compute gravitational velocities
-!
- CALL BLOWSNW_VELGRAV1D(ZBET, ZRG, PT, PRHOA,PP, ZVGK)
-!
-!-------------------------------------------------------------------------------
-!
-!
-!* 3. Compute sublimation rate based on formulation of Thorpe and Masson (66)
-! Note that no integration is computed on the particle spectra
-! ---------------
-!
-!-------------------------------------------------------------------------------
-!
-! Nusselt Number using effective radius
-ZNU(:,:) = NUSSELT(ZRG(:,1:(KLVL-1)),ZMU(:,1:(KLVL-1)),ZVGK(:,1:(KLVL-1),2))
-
-ZSUBL_MASS(:,:) = 3*PBLOWSNW(:,1:(KLVL-1),2)*ZNU(:,:)*ZUSI(:,:) &
- /(2*ZAI(:,:)*XRHOLI*ZR_EFF**2)
-
-
-ELSE
-
-
-!-------------------------------------------------------------------------------
-!
-!
-!* 2. Compute settling velocity based on the size distribution of the previous time
-! step: used as ventilation velocity
-! ---------------
-!
-!-------------------------------------------------------------------------------
-
-!
-! Convert canopy variables in _/m3
-!
-DO JK=1,KLVL
- DO JSV=1,KSNW
- ZSNW(:,JK,JSV)=PBLOWSNW(:,JK,JSV)*PRHOA(:)
- ENDDO
-ENDDO
-!
-!* compute BETA, RG and moments
-!
-CALL SNOWMOMENT2SIZE(ZSNW, PBETA1D=ZBET, PRG1D=ZRG )
-!
-!* compute gravitational velocities
-!
-CALL BLOWSNW_VELGRAV1D(ZBET, ZRG, PT, PRHOA,PP, ZVGK)
-!
-!-------------------------------------------------------------------------------
-!
-!
-!* 3. Compute sublimation rate based on formulation of Thorpe and Masson (66)
-! ---------------
-!
-!-------------------------------------------------------------------------------
-!
-! Nusselt Number using mean radius of particle size distribution
-ZNU(:,:) = NUSSELT(ZRG(:,1:(KLVL-1)),ZMU(:,1:(KLVL-1)),ZVGK(:,1:(KLVL-1),2))
-
-
-! mass averaged sublimation rate follows Dery and Yau (1999) and avoids
-! numerical integration over the particle spectrum
-ZSUBL_MASS(:,:) = PBLOWSNW(:,1:(KLVL-1),2)*ZNU(:,:)*ZUSI(:,:)/ &
- (ZAI(:,:)*2*XRHOLI*ZRG(:,1:(KLVL-1))**2)
-!
-!
-ENDIF
-! Restriction of ZSUM_SUBL to insure coherence between vapor and blowing snow
-! mixing ratio
-!
-DO JI=1,KI
- DO JK=1,(KLVL-1)
- ZSUBL_MASS(JI,JK) = MIN( PQ(JI,JK)/(PRHOA(JI)*PTSTEP),ZSUBL_MASS(JI,JK))* &
- (0.5+SIGN(0.5,ZSUBL_MASS(JI,JK)))- &
- MIN(PBLOWSNW(JI,JK,2)/PTSTEP,ABS(ZSUBL_MASS(JI,JK)))* &
- (0.5-SIGN(0.5,ZSUBL_MASS(JI,JK)))
- ZSUBL_MASS(JI,JK)=MIN(0.,ZSUBL_MASS(JI,JK)) ! Sink of snow
-!
-! number-averaged sublimation rate
-! Change in concentration rate Sn = Sb*N/qb (Dery and Yau,2000)
-!
- IF(ZSUBL_MASS(JI,JK)<0.) THEN
- ZSUBL_NUM(JI,JK) = ZSUBL_MASS(JI,JK)*PBLOWSNW(JI,JK,1)/PBLOWSNW(JI,JK,2)
- END IF
- ZCOL_SUBL(JI) = ZCOL_SUBL(JI)+ZSUBL_MASS(JI,JK)*PRHOA(JI)*PDZ(JI,JK)
-!* Compute modification of heat and vapor fluxes due to
-! sublimation of blowing snow particles in Canopy
-! These fluxes are then sent to MNH
- PSFTQ(JI) = PSFTQ(JI) - ZSUBL_MASS(JI,JK)*PRHOA(JI)*PDZ(JI,JK)
- PSFTH(JI) = PSFTH(JI) + ZSUBL_MASS(JI,JK)*PRHOA(JI)*PDZ(JI,JK)* &
- (XLSTT+(XCPV-XCI)*(PT(JI,JK)-XTT))/ZEXN(JI,JK)
- END DO
-! ZT(JI) = SUM(PT(JI,1:(KLVL-1)))/(KLVL-1)
-END DO
-!
-! Store number and mass sublimation rate
-PBLOWSNWS(:,1:(KLVL-1),1) = ZSUBL_NUM(:,:)
-PBLOWSNWS(:,1:(KLVL-1),2) = ZSUBL_MASS(:,:)
-
-! Store integrated sublimation rate in mmSWE/day
-PSNW_SUBL(:,2) = ZCOL_SUBL(:)*3600*24/XRHOLW*1000
-! Store accumulated sublimation rate in mmSWE
-PSNW_SUBL(:,3) = PSNW_SUBL(:,3)+ZCOL_SUBL(:)*PTSTEP
-
-CONTAINS
-
-FUNCTION NUSSELT(PR,PMU,PVEL_VENT) RESULT(PNU)
-!
-!! PURPOSE
-!! -------
-! Calculate the Nusselt number for a given particle radius
-! Formulation based on Lee (1975)
-!
-!
-IMPLICIT NONE
-!
-!* 0.1 declarations of arguments
-!
-REAL, DIMENSION(:,:), INTENT(IN) :: PR ! (m)
-REAL, DIMENSION(:,:), INTENT(IN) :: PMU ! (m2/s)
-REAL, DIMENSION(:,:), INTENT(IN) :: PVEL_VENT ! (m/s)
-!
-REAL, DIMENSION(SIZE(PMU,1),SIZE(PMU,2)) :: PNU ! (m/s)
-!
-!
-!* 0.2 declaration of local variables
-!
-REAL, DIMENSION(SIZE(PMU,1),SIZE(PMU,2)) :: ZRE
-!
-!
-!* 1 Calculate Reynolds number
-!
-ZRE(:,:) = 2*PR(:,:)*PVEL_VENT(:,:)/PMU(:,:)
-!
-!* 2 Calculate Nusselt number
-!
-WHERE(ZRE(:,:)<10)
- PNU(:,:) = 1.79+0.606*ZRE(:,:)**0.5
-ELSEWHERE
- PNU(:,:) = 1.88+0.580*ZRE(:,:)**0.5
-END WHERE
-
-END FUNCTION NUSSELT
-
-END SUBROUTINE CANOPY_BLOWSNW_SUBL
diff --git a/src/SURFEX/canopy_evol_blowsnw.f90 b/src/SURFEX/canopy_evol_blowsnw.f90
deleted file mode 100644
index 088ca2aff..000000000
--- a/src/SURFEX/canopy_evol_blowsnw.f90
+++ /dev/null
@@ -1,274 +0,0 @@
-! #########################################
- SUBROUTINE CANOPY_EVOL_BLOWSNW(KI,KLVL,KSNW,PTSTEP,PSNWA,PK, &
- PSFLUX_SNW,PZ,PDZ,PDZF,PSNW,PRHOA,PSNWDEP, &
- PTHT,PPABST,PSNW_SUBL,PTKE)
-! #########################################
-!
-!!**** *CANOPY_EVOL_BLOWSNW* - evolution of blowing snow variables in canopy
-!!
-!!
-!! PURPOSE
-!! -------
-!!
-!!** METHOD
-!! ------
-!!
-!! EXTERNAL
-!! --------
-!!
-!!
-!! IMPLICIT ARGUMENTS
-!! ------------------
-!!
-!! REFERENCE
-!! ---------
-!!
-!!
-!! AUTHOR
-!! ------
-!! V. Vionnet *Meteo France*
-!!
-!! MODIFICATIONS
-!! -------------
-!! Original 04/2012
-!-------------------------------------------------------------------------------
-!
-!* 0. DECLARATIONS
-! ------------
-!
-USE MODD_BLOWSNW_n
-USE MODD_BLOWSNW_SURF
-
-USE MODE_BLOWSNW_SURF
-
-USE MODI_TRIDIAG_GROUND
-USE MODI_BLOWSNW_VELGRAV1D
-USE MODI_BLOWSNW_DIFFK
-
-!
-IMPLICIT NONE
-!
-!* 0.1 Declarations of arguments
-! -------------------------
-!
-INTEGER, INTENT(IN) :: KI ! number of horizontal points
-INTEGER, INTENT(IN) :: KLVL ! number of levels in canopy
-INTEGER, INTENT(IN) :: KSNW ! number of drifting snow variables
-REAL, INTENT(IN) :: PTSTEP ! time-step (s)
-REAL, DIMENSION(KI,KSNW), INTENT(IN) :: PSNWA ! Blowing snow variables at forcing levels (__ /kg)
-REAL, DIMENSION(KI,KLVL), INTENT(IN) :: PK ! mixing exchange coefficient (m2/s)
-REAL, DIMENSION(KI,KSNW), INTENT(INOUT) :: PSFLUX_SNW ! surface flux w'snw' (mkg/kg/s)
-REAL, DIMENSION(KI,KLVL), INTENT(IN) :: PZ ! heights of canopy levels (m)
-REAL, DIMENSION(KI,KLVL), INTENT(IN) :: PDZ ! deltaZ between canopy half levels,
-! ! located at full levels (m)
-REAL, DIMENSION(KI,KLVL), INTENT(IN) :: PDZF ! deltaZ between canopy (full) levels,
-! ! located at half levels (m)
-REAL, DIMENSION(KI,KLVL,KSNW), INTENT(INOUT) :: PSNW ! drifting snow at canopy levels (kg/kg)
-REAL, DIMENSION(KI,KLVL,KSNW), INTENT(IN) :: PSNW_SUBL ! drifting snow sublimation rate at canopy levels (kg/kg/s)
- ! zero when sublimation is not activated
-REAL, DIMENSION(KI), INTENT(IN) :: PRHOA ! air density (kg/m3)
-REAL, DIMENSION(KI,KLVL), INTENT(IN) :: PTHT ! air temperature (K)
-REAL, DIMENSION(KI,KLVL), INTENT(IN) :: PPABST ! air pressure (Pa)
-REAL, DIMENSION(KI,KLVL), INTENT(IN) :: PTKE ! air pressure (Pa)
-REAL, DIMENSION(KI,KSNW), INTENT(OUT) :: PSNWDEP ! sedimentation flux at the bottom of Canopy
-
-!
-!
-!* 0.2 Declarations of local variables
-! -------------------------------
-!
-INTEGER :: JLAYER ! loop counter on layers
-INTEGER :: JSV ! loop counter on blowing snow variables
-INTEGER :: JI ! loop counter on grid points
-!
-REAL, DIMENSION(KI) :: ZZREF ! height of forcing level
-REAL, DIMENSION(KI,KSNW) :: ZSNWC ! Blowing Snow content in canopy layers (__ /m3)
-
-REAL, DIMENSION(KI,KLVL) :: ZBET ! Scale parameter of the gamma distribution (m)
-REAL, DIMENSION(KI,KLVL,KSNW) :: ZVGK ! Terminal fallspeed (m/s)
-REAL, DIMENSION(KI,KLVL,KSNW) :: ZKSNW ! particle eddy diffusivity
-REAL, DIMENSION(KI,KLVL) :: ZRG ! Mean radius (m)
-REAL, DIMENSION(KI,KLVL-1,KSNW) :: ZA,ZB,ZC,ZRHS ! Term of the tridiagonal matrix
-REAL, DIMENSION(KI,KLVL,KSNW) :: ZSNW ! work variable : pot. temp at futur instant
-! ! (or past at the end of the routine)
-LOGICAL, DIMENSION(KI) :: GEMIS ! Logical=TRUE if snow is emitted in the atmosphere for the
- ! first time
-INTEGER :: IPRINT
-!
-!-------------------------------------------------------------------------------
-!
-!
-!* 1. initializations
-! ---------------
-!
-IPRINT = 1
-!
-GEMIS(:)= .FALSE.
-!
-!-------------------------------------------------------------------------------
-!
-!
-!* 2. Compute settling velocity based on the size distribution of the previous time
-! step
-! ---------------
-!
-!-------------------------------------------------------------------------------
-!
-! Convert canopy variables in _/m3
-!
-DO JLAYER=1,KLVL
- DO JSV=1,KSNW
- ZSNW(:,JLAYER,JSV)=PSNW(:,JLAYER,JSV)*PRHOA(:)
- ENDDO
-ENDDO
-
-!
-!* compute BETA, RG and moments
-!
-CALL SNOWMOMENT2SIZE(ZSNW, PBETA1D=ZBET, PRG1D=ZRG )
-!
-!* Initialize profile of radius if no snow is present in the atmosphere
-! We use Pomeroy's theoretical profile.
-!
-
-DO JI=1,KI
- IF(PSFLUX_SNW(JI,2)>0 .AND. ALL(PSNW(JI,:,2)==0.)) THEN
- GEMIS(JI) = .TRUE.
- DO JLAYER=1,KLVL
- ZRG(JI,JLAYER)=XEMIRADIUS_SNW*(PZ(JI,JLAYER)/0.05)**(-0.258)
- ZBET(JI,JLAYER) = ZRG(JI,JLAYER)/XEMIALPHA_SNW
- END DO
- END IF
-END DO
-!* compute gravitational velocities
-!
-!
-CALL BLOWSNW_VELGRAV1D(ZBET, ZRG, PTHT, PRHOA,PPABST, ZVGK)
-
-!
-! Mettre en place cas spécial avec flux positif et pas de neige dans l'atmo
-! initialisation du profil de rayon et de vitesse de chute suivant la méthode
-! de Déry
-!
-!* compute particle eddy diffusivity
-!
-!
-CALL BLOWSNW_DIFFK(PK, PTKE, ZVGK,KSNW, ZKSNW)
-!
-!
-!-------------------------------------------------------------------------------
-!
-!
-!* 3. Initialize coefficients of the tridiagonal matrix
-! ---------------
-!
-!-------------------------------------------------------------------------------
-!
-! Lower boundary: turbulent flux PSFLUX_SNW at the surface
-!
-DO JSV = 1,KSNW
- ZA(:,1,JSV) = 0.
- ZB(:,1,JSV) = 1+(ZKSNW(:,2,JSV)/PDZF(:,2)+ZVGK(:,1,JSV))*PTSTEP/PDZ(:,1)
- ZC(:,1,JSV) = -(ZKSNW(:,2,JSV)/PDZF(:,2)+ZVGK(:,2,JSV))*PTSTEP/PDZ(:,1)
- ZRHS(:,1,JSV)= PSNW(:,1,JSV)+PSFLUX_SNW(:,JSV)*PTSTEP/PDZ(:,1)+ &
- PSNW_SUBL(:,1,JSV)*PTSTEP
-END DO
-!
-! Upper boundary: imposed concentration at level KLVL
-!
-DO JSV = 1,KSNW
- ZC(:,KLVL-1,JSV) = 0.
- ZRHS(:,KLVL-1,JSV)= PSNW(:,KLVL-1,JSV)+PSNW_SUBL(:,KLVL-1,JSV)*PTSTEP+ &
- PSNW(:,KLVL,JSV)*PTSTEP/ &
- PDZ(:,KLVL-1)*(ZVGK(:,KLVL,JSV)+ZKSNW(:,KLVL,JSV)/PDZF(:,KLVL))
-END DO
-!
-! Values at other levels
-!
-DO JSV = 1,KSNW
- DO JLAYER=2,(KLVL-1)
- ZA(:,JLAYER,JSV)= -(ZKSNW(:,JLAYER,JSV)/PDZF(:,JLAYER))*PTSTEP/PDZ(:,JLAYER)
- ZB(:,JLAYER,JSV)=1+PTSTEP/PDZ(:,JLAYER)*(ZKSNW(:,JLAYER+1,JSV)/PDZF(:,JLAYER+1)+ &
- ZKSNW(:,JLAYER,JSV)/PDZF(:,JLAYER)+ZVGK(:,JLAYER,JSV))
- IF(JLAYER<(KLVL-1)) THEN
- ZC(:,JLAYER,JSV) = -(ZKSNW(:,JLAYER+1,JSV)/PDZF(:,JLAYER+1)+ZVGK(:,JLAYER+1,JSV)) &
- *PTSTEP/PDZ(:,JLAYER)
- ZRHS(:,JLAYER,JSV) = PSNW(:,JLAYER,JSV)+PSNW_SUBL(:,JLAYER,JSV)*PTSTEP
- END IF
- END DO
-END DO
-!
-!-------------------------------------------------------------------------------
-!
-!
-!* 4. Solve tridiagonal system
-! ---------------
-!
-!-------------------------------------------------------------------------------
-!
-DO JSV = 1,KSNW
-
- CALL TRIDIAG_GROUND(ZA(:,:,JSV),ZB(:,:,JSV),ZC(:,:,JSV),ZRHS(:,:,JSV), &
- ZSNW(:,1:KLVL-1,JSV))
-END DO
-!
-!
-!-------------------------------------------------------------------------------
-!
-!* 5. updated fluxes
-! ----------------------
-!
-! Compute net flux sent to Meso-NH: use top of canopy instead of surface
-!
-DO JI=1,KI
- DO JSV=1,KSNW
- IF(GEMIS(JI)) THEN
- PSFLUX_SNW(JI,JSV) = MAX(-ZKSNW (JI,KLVL,JSV) * (PSNW(JI,KLVL,JSV)-ZSNW(JI,KLVL-1,JSV))/PDZF(JI,KLVL)* &
- (1-ZVGK(JI,KLVL,JSV)*PTSTEP/(2*PZ(JI,KLVL))),0.)
- ELSE
- PSFLUX_SNW(JI,JSV) = -ZKSNW (JI,KLVL,JSV) * (PSNW(JI,KLVL,JSV)-ZSNW(JI,KLVL-1,JSV))/PDZF(JI,KLVL)- &
- ZVGK(JI,KLVL,JSV)*PSNW(JI,KLVL,JSV)
- ENDIF
- PSFLUX_SNW(JI,JSV) = PSFLUX_SNW(JI,JSV)*PRHOA(JI) ! Convert flux in kg_{snow}/m2/s
- END DO
-END DO
-!
-! Update sedimentation flux send to Crocus at the next time step
-!
-
-DO JSV=1,KSNW
- PSNWDEP(:,JSV) = ZVGK(:,1,JSV)*ZSNW(:,1,JSV)*PRHOA(:)
-END DO
-!
-! Store sedimentation flux at the bottom of Canopy as a diagnostic
-!XSNSED(:,1:2) = -PSNWDEP(:,1:2) ! Instantaneous fluxes (number and mass)
-!XSNSED(:,3) = XSNSED(:,3)+XSNSED(:,2)*PTSTEP ! Accumulated flux (mass)
-
-
-!
-!-------------------------------------------------------------------------------
-!
-!* 6. New value of blowing snow variables (at full levels)
-! ----------------------------------
-!
-PSNW(:,1:(KLVL-1),:) = ZSNW(:,1:(KLVL-1),:)
-!
-!-------------------------------------------------------------------------------
-!
-!* 7. Update total Canopy variables
-! ----------------------------------
-!
-!ZSNWC =0.
-! Compute number and mass content of Canopy layer.
-!DO JI=1,KI
-! ZZREF(JI) = SUM(PDZ(JI,1:(KLVL-1)))
-! DO JSV=1,KSNW
-! DO JLAYER=1,(KLVL-1)
-! ZSNWC(JI,JSV) = ZSNWC(JI,JSV)+PSNW(JI,JLAYER,JSV)*PDZ(JI,JLAYER)*PRHOA(JI)
-! END DO
-! PSNOW_CANO(JI,JSV)=ZSNWC(JI,JSV)/ZZREF(JI)
-! END DO
-!END DO
-!
-
-END SUBROUTINE CANOPY_EVOL_BLOWSNW
diff --git a/src/SURFEX/gamma.f90 b/src/SURFEX/gamma.f90
deleted file mode 100644
index cd7e04fbd..000000000
--- a/src/SURFEX/gamma.f90
+++ /dev/null
@@ -1,224 +0,0 @@
-!MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
-!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
-!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
-!MNH_LIC for details. version 1.
-!########################
- MODULE MODI_GAMMA
-!########################
-!
-INTERFACE GAMMA
-!
-FUNCTION GAMMA_X0D(PX) RESULT(PGAMMA)
-REAL, INTENT(IN) :: PX
-REAL :: PGAMMA
-END FUNCTION GAMMA_X0D
-!
-FUNCTION GAMMA_X1D(PX) RESULT(PGAMMA)
-REAL, DIMENSION(:), INTENT(IN) :: PX
-REAL, DIMENSION(SIZE(PX)) :: PGAMMA
-END FUNCTION GAMMA_X1D
-!
-END INTERFACE
-END MODULE MODI_GAMMA
-!
-!--------------------------------------------------------------------------
-!
-!
-!* 1. FUNCTION GAMMA FOR SCALAR VARIABLE
-!
-!
-! ######################################
- FUNCTION GAMMA_X0D(PX) RESULT(PGAMMA)
-! ######################################
-!
-!
-!!**** *GAMMA * - Gamma function
-!!
-!!
-!! PURPOSE
-!! -------
-! The purpose of this function is to compute the Generalized gamma
-! function of its argument.
-!
-!
-!!** METHOD
-!! ------
-!!
-!! EXTERNAL
-!! --------
-!! NONE
-!!
-!! IMPLICIT ARGUMENTS
-!! ------------------
-!! None
-!!
-!! REFERENCE
-!! ---------
-!! Press, Teukolsky, Vetterling and Flannery: Numerical Recipes, 206-207
-!!
-!! AUTHOR
-!! ------
-!! Jean-Pierre Pinty *LA/OMP*
-!!
-!! MODIFICATIONS
-!! -------------
-!! Original 7/11/95
-!! C. Barthe 9/11/09 add a function for 1D arguments
-!
-!* 0. DECLARATIONS
-! ------------
-!
-IMPLICIT NONE
-!
-!* 0.1 declarations of arguments and result
-!
-REAL, INTENT(IN) :: PX
-REAL :: PGAMMA
-!
-!* 0.2 declarations of local variables
-!
-INTEGER :: JJ ! Loop index
-REAL :: ZSER,ZSTP,ZTMP,ZX,ZY,ZCOEF(6)
-REAL :: ZPI
-!
-!-------------------------------------------------------------------------------
-!
-!* 1. SOME CONSTANTS
-! --------------
-!
-ZCOEF(1) = 76.18009172947146
-ZCOEF(2) =-86.50532032941677
-ZCOEF(3) = 24.01409824083091
-ZCOEF(4) = -1.231739572450155
-ZCOEF(5) = 0.1208650973866179E-2
-ZCOEF(6) = -0.5395239384953E-5
-ZSTP = 2.5066282746310005
-!
-ZPI = 3.141592654
-!
-!-------------------------------------------------------------------------------
-!
-!* 2. COMPUTE GAMMA
-! -------------
-!
-IF (PX .LT. 0.) THEN
- ZX = 1. - PX
-ELSE
- ZX = PX
-END IF
-ZY = ZX
-ZTMP = ZX + 5.5
-ZTMP = (ZX + 0.5) * ALOG(ZTMP) - ZTMP
-ZSER = 1.000000000190015
-!
-DO JJ = 1, 6
- ZY = ZY + 1.0
- ZSER = ZSER + ZCOEF(JJ) / ZY
-END DO
-!
-IF (PX .LT. 0.) THEN
- PGAMMA = ZPI / SIN(ZPI*PX) / EXP(ZTMP + ALOG(ZSTP*ZSER/ZX))
-ELSE
- PGAMMA = EXP(ZTMP + ALOG(ZSTP*ZSER/ZX))
-END IF
-RETURN
-!
-END FUNCTION GAMMA_X0D
-!
-!-------------------------------------------------------------------------------
-!
-!
-!* 1. FUNCTION GAMMA FOR 1D ARRAY
-!
-!
-! ######################################
- FUNCTION GAMMA_X1D(PX) RESULT(PGAMMA)
-! ######################################
-!
-!
-!!**** *GAMMA * - Gamma function
-!!
-!!
-!! PURPOSE
-!! -------
-! The purpose of this function is to compute the Generalized gamma
-! function of its argument.
-!
-!
-!!** METHOD
-!! ------
-!!
-!! EXTERNAL
-!! --------
-!! NONE
-!!
-!! IMPLICIT ARGUMENTS
-!! ------------------
-!! None
-!!
-!! REFERENCE
-!! ---------
-!! Press, Teukolsky, Vetterling and Flannery: Numerical Recipes, 206-207
-!!
-!! AUTHOR
-!! ------
-!! Jean-Pierre Pinty *LA/OMP*
-!!
-!! MODIFICATIONS
-!! -------------
-!! Original 7/11/95
-!!
-!-------------------------------------------------------------------------------
-!
-!* 0. DECLARATIONS
-! ------------
-!
-IMPLICIT NONE
-!
-!* 0.1 declarations of arguments and result
-!
-REAL, DIMENSION(:), INTENT(IN) :: PX
-REAL, DIMENSION(SIZE(PX)) :: PGAMMA
-!
-!* 0.2 declarations of local variables
-!
-INTEGER :: JJ ! Loop index
-REAL, DIMENSION(SIZE(PX)) :: ZSER,ZSTP,ZTMP,ZX,ZY
-REAL :: ZCOEF(6)
-REAL :: ZPI
-!
-!-------------------------------------------------------------------------------
-!
-!* 1. SOME CONSTANTS
-! --------------
-!
-ZCOEF(1) = 76.18009172947146
-ZCOEF(2) =-86.50532032941677
-ZCOEF(3) = 24.01409824083091
-ZCOEF(4) = -1.231739572450155
-ZCOEF(5) = 0.1208650973866179E-2
-ZCOEF(6) = -0.5395239384953E-5
-ZSTP = 2.5066282746310005
-!
-ZPI = 3.141592654
-ZX(:) = PX(:)
-WHERE ( PX(:)<0.0 )
- ZX(:) = 1.- PX(:)
-END WHERE
-ZY(:) = ZX(:)
-ZTMP(:) = ZX(:) + 5.5
-ZTMP(:) = (ZX(:) + 0.5)*ALOG(ZTMP(:)) - ZTMP(:)
-ZSER(:) = 1.000000000190015
-!
-DO JJ = 1 , 6
- ZY(:) = ZY(:) + 1.0
- ZSER(:) = ZSER(:) + ZCOEF(JJ)/ZY(:)
-END DO
-!
-PGAMMA(:) = EXP( ZTMP(:) + ALOG( ZSTP*ZSER(:)/ZX(:) ) )
-WHERE ( PX(:)<0.0 )
- PGAMMA(:) = ZPI/SIN(ZPI*PX(:))/PGAMMA(:)
-END WHERE
-RETURN
-!
-END FUNCTION GAMMA_X1D
diff --git a/src/SURFEX/gamma_inc_low.f90 b/src/SURFEX/gamma_inc_low.f90
deleted file mode 100644
index 1329e2ce7..000000000
--- a/src/SURFEX/gamma_inc_low.f90
+++ /dev/null
@@ -1,130 +0,0 @@
-!####################
-MODULE MODI_GAMMA_INC_LOW
-!####################
-!
-INTERFACE
-!
-FUNCTION GAMMA_INC_LOW(PA,PX) RESULT(PGAMMA_INC_LOW)
-REAL, INTENT(IN) :: PA
-REAL, INTENT(IN) :: PX
-REAL :: PGAMMA_INC_LOW
-END FUNCTION GAMMA_INC_LOW
-!
-END INTERFACE
-!
-END MODULE MODI_GAMMA_INC_LOW
-! #############################################
- FUNCTION GAMMA_INC_LOW(PA,PX) RESULT(PGAMMA_INC_LOW)
-! #############################################
-!
-!
-!!**** *GAMMA_INC_LOW * - Generalized gamma function
-!!
-!!
-!! PURPOSE
-!! -------
-! The purpose of this function is to compute the generalized
-!! lower incomplete Gamma function of its argument.
-!!
-!! /X
-!! |
-!! GAMMA_INC_LOW(A,X)= ---- | Z**(A-1) EXP(-Z) dZ with A >0
-!! |
-!! /0
-!!
-!!** METHOD
-!! ------
-!!
-!! EXTERNAL
-!! --------
-!! NONE
-!!
-!! IMPLICIT ARGUMENTS
-!! ------------------
-!! MODULE MODI_GAMMA : computation of the Gamma function
-!!
-!! REFERENCE
-!! ---------
-!! U. Blahak : Efficient approximation of the incomplete gamma function for
-!! use in cloud model applications, GMD, 2010
-!!
-!!
-!! AUTHOR
-!! ------
-!! V. Vionnet (CNRM/GMME)
-!!
-!! MODIFICATIONS
-!! -------------
-!! Original 20/09/10
-!
-!* 0. DECLARATIONS
-! ------------
-!
-USE MODI_GAMMA
-!
-IMPLICIT NONE
-!
-!* 0.1 declarations of arguments and result
-!
-REAL, INTENT(IN) :: PA
-REAL, INTENT(IN) :: PX
-REAL :: PGAMMA_INC_LOW
-!
-!* 0.2 declarations of local variables
-!
-REAL :: ZP(6), ZQ(4), ZR(4), ZS(5)
-REAL :: ZC(4)
-REAL :: ZWORK
-!
-!* 0.3 initializations of local variables
-!
-ZP(1) = 9.4368392235E-3
-ZP(2) = -1.0782666481E-4
-ZP(3) = -5.8969657295E-6
-ZP(4) = 2.8939523781E-7
-ZP(5) = 1.0043326298E-1
-ZP(6) = 5.5637848465E-1
-
-ZQ(1) = 1.1464706419E-1
-ZQ(2) = 2.6963429121
-ZQ(3) = -2.9647038257
-ZQ(4) = 2.1080724954
-
-ZR(1) = 0.0
-ZR(2) = 1.1428716184
-ZR(3) = -6.6981186438E-3
-ZR(4) = 1.0480765092E-4
-
-ZS(1) = 1.0356711153
-ZS(2) = 2.3423452308
-ZS(3) = -3.6174503174E-1
-ZS(4) = -3.1376557650
-ZS(5) = 2.9092306039
-!
-!* 1 Compute coefficients
-!
-IF( (PX.LT.0.0).OR.(PA.LE.0.0) ) THEN
- PRINT *,' BAD ARGUMENTS IN GAMMA_INC_LOW'
-!callabortstop
-CALL ABORT
- STOP
-END IF
-!
-!
-ZC(1) = 1.+ZP(1)*PA+ZP(2)*PA**2+ZP(3)*PA**3+ZP(4)*PA**4+ZP(5)*(EXP(-ZP(6)*PA)-1)
-!
-ZC(2) = ZQ(1) + ZQ(2)/PA + ZQ(3)/PA**2 + ZQ(4)/PA**3
-!
-ZC(3) = ZR(1)+ZR(2)*PA+ZR(3)*PA**2+ZR(4)*PA**3
-!
-ZC(4) = ZS(1) + ZS(2)/PA + ZS(3)/PA**2 + ZS(4)/PA**3 + ZS(5)/PA**4
-!
-!* 2 Compute final results
-!
-ZWORK = 0.5+0.5*TANH(ZC(2)*(PX-ZC(3)))
-
-PGAMMA_INC_LOW = EXP(-PX)* PX**PA * (1./PA +ZC(1)*PX/(PA*(PA+1.))+(ZC(1)*PX)**2/(PA*(PA+1.)*(PA+2.))) &
- * (1.-ZWORK) + GAMMA(PA)*ZWORK*(1.-ZC(4)**(-PX))
-RETURN
-!
-END FUNCTION GAMMA_INC_LOW
diff --git a/src/SURFEX/modd_blowsnw_sedim_lkt1d.f90 b/src/SURFEX/modd_blowsnw_sedim_lkt1d.f90
deleted file mode 100644
index e2a909cbd..000000000
--- a/src/SURFEX/modd_blowsnw_sedim_lkt1d.f90
+++ /dev/null
@@ -1,54 +0,0 @@
-
-!! ########################
- MODULE MODD_BLOWSNW_SEDIM_LKT1D
-!! ########################
-!!
-!! PURPOSE
-!! -------
-!!
-!! Purpose: Contains look up tables for settling velocity of drifitng snow particles
-!! The parameters to be looked up are:
-!! 1) Number-averaged settling velocity
-!! 2) Mass-averaged settling velocity
-!!
-!! All values are pre-calculated using matlab.
-!!
-!!
-!! METHOD
-!! ------
-!!
-!!
-!! REFERENCE
-!! ---------
-!! Based on MODD_DUST_OPT_LKT (Pierre Tulet)
-!!
-!!
-!! AUTHOR
-!! ------
-!! Vincent VIONNET (CNRM)
-!!
-!!
-!! MODIFICATIONS
-!! -------------
-!!
-!!--------------------------------------------------------------------
-!! DECLARATIONS
-!! ------------
-
- IMPLICIT NONE
- PUBLIC
-
- INTEGER, PARAMETER :: NMAX_RADIUS_LKT1D=196 !Max number of radii in look up tables ()
- INTEGER, PARAMETER :: NMAX_PRESSURE_LKT1D=4 !Max number of pressure in lkt
-
- !Declaration of the look up tables
- REAL, DIMENSION(NMAX_RADIUS_LKT1D,NMAX_PRESSURE_LKT1D) :: XNUMB_SPEED_LKT1D
- REAL, DIMENSION(NMAX_RADIUS_LKT1D,NMAX_PRESSURE_LKT1D) :: XMASS_SPEED_LKT1D
-
- !Declaration of the max and min values taken into account in the tables
- REAL, PARAMETER :: XRADIUS_LKT1D_MIN = 5 ![um] smallest number median radius taken into account
- REAL, PARAMETER :: XRADIUS_LKT1D_MAX = 200 ![um] largest number median radius taken into account
- REAL, PARAMETER :: XPRESSURE_LKT1D_MIN = 45000 ![Pa] smallest pressure coefficient taken into account
- REAL, PARAMETER :: XPRESSURE_LKT1D_MAX = 105000 ![Pa] largest pressure coefficient taken into account
-
-END MODULE MODD_BLOWSNW_SEDIM_LKT1D
diff --git a/src/SURFEX/modd_blowsnw_surf.f90 b/src/SURFEX/modd_blowsnw_surf.f90
deleted file mode 100644
index 343848755..000000000
--- a/src/SURFEX/modd_blowsnw_surf.f90
+++ /dev/null
@@ -1,77 +0,0 @@
-!-----------------------------------------------------------------
- MODULE MODD_BLOWSNW_SURF
-!
-! Contains parameters for blowing snow simulation
-!
- IMPLICIT NONE
-!
-! Snow density : cf Liston et Sturm (1998)
- REAL, PARAMETER :: XRHO_DEP=280.0
-!
-! Deposited snow grains are round and small : s=1,r=0.4mm
- REAL,PARAMETER :: XSDEP_GRA1=99.0
- REAL,PARAMETER :: XSDEP_GRA2=0.0004
-!
-! Minimal mean radius (um)
- REAL,PARAMETER :: XINIRADIUS_SNW = 5.e-6
-! Minimum allowed number concentration (#/m3)
- REAL,PARAMETER :: XN0MIN_SNW = 1
-
-! Parameters used in KC02 parameterization for settling velocity
- REAL,PARAMETER :: XC0 = 0.29
- REAL,PARAMETER :: XDELTA0 = 9.06
-
-! Parameters used in Mitchell (96) parameterization for settling velocity
- REAL,PARAMETER :: XAM1 = 0.04394
- REAL,PARAMETER :: XAM2 = 0.06049
- REAL,PARAMETER :: XAM3 = 0.2072
- REAL,PARAMETER :: XBM1 = 0.970
- REAL,PARAMETER :: XBM2 = 0.831
- REAL,PARAMETER :: XBM3 = 0.638
- REAL,PARAMETER :: XBESTL_1 = 10.0
- REAL,PARAMETER :: XBESTL_2 = 585.
-
-! Grain size distribution is a two parameter gamma distribution
- REAL :: XEMIRADIUS_SNW
- REAL :: XEMIALPHA_SNW
-!
- LOGICAL :: LSNOW_SALT ! Flag to fix snow concentration in the saltation layer
- ! Snow concentration in the saltation layer (kg_{snow}/m3_{air}) if LSNOW_SALT = T
- REAL :: XCONC_SALT
- !
- LOGICAL :: LSNOW_PRECIP ! Flag to impose uniform and constant precipitation rate
- ! Fixed snow precipitation rate SWE (kg/ (m2 s)) if LSNOW_PRECIP = T
- REAL :: XSNOW_PRECIP
- LOGICAL :: LBLOWSNW_ADV ! Flag to account for advection effects on
- ! total mass and number in Canopy variables and
- ! to compute divergence of
- ! saltation flux
- LOGICAL :: LBLOWSNW_CANOSUBL ! Flag to activate sublimation of Canopy variables
-
- LOGICAL :: LBLOWSNW_CANODIAG ! Flag to get additional diagnostic at Canopy levels : mean radius
- ! and number and mass variables in _/m3
-
- LOGICAL :: LSNOW_WIND ! Flag to activate effects of snow particle on wind profile
- ! Increase in surface rougnhess due to saltation following Dover (1993) z0_s = z0_ns*(u*/uth*)²
- REAL :: XSNOW_ROUGHNESS ! = 0 not activated; =1 activated
- ! Buoyancy effect in presence of suspended particles of blowing snow.
- REAL :: XSNOW_BUOYANCY ! = 0 not activated; =1 activated
- CHARACTER(LEN=4) :: CSNOW_SALT ! Paramaterization to compute particle flux in saltation
-! 'POME': Pomeroy and Gray, 1990
-! 'SORE': Sorensen (1991) : used at SLF before Doorshot's model
-! 'MANN': Concentration in the saltation layer is computed according to Mann
-! parameterization for particle number density just above the ground (10 cm)
- CHARACTER(LEN=4) :: CSNOW_SEDIM ! Paramaterization to compute settling velocity
-! 'CARR': follow Carrier's drag coefficient
-! 'TABC': based on tabuleted values of Carrier's formulation for alpha=3
-! 'MITC': based on Mitchell's formulation for settling velocity of ice spheres
-! 'NONE': sedimentation is desactivated (for test only!)
-!
- REAL :: XRSNOW_SBL ! Ratio between diffusion coefficient for scalar
- ! variables and blowing snow variables
- ! RSNOW_SBL = KSCA/KSNOW = 4.
- ! See Vionnet (PhD, 2012, In French) and Vionnet et al (TC, 2014)
- ! for a complete dicsussion
-!
- END MODULE MODD_BLOWSNW_SURF
-
diff --git a/src/SURFEX/modd_blowsnwn.f90 b/src/SURFEX/modd_blowsnwn.f90
deleted file mode 100644
index 346e56dae..000000000
--- a/src/SURFEX/modd_blowsnwn.f90
+++ /dev/null
@@ -1,51 +0,0 @@
-MODULE MODD_BLOWSNW_n
-
-!Purpose:
-!Declare variables and constants necessary to do the blowing snow calculations
-!Here are only the variables which depend on the grid!
-
-!Author: Vincent Vionnet
-! based on modd_dstn.F90
-!
-USE YOMHOOK ,ONLY : LHOOK, DR_HOOK
-USE PARKIND1 ,ONLY : JPRB
-!
-IMPLICIT NONE
-
-TYPE BLOWSNW_t
-
- REAL, DIMENSION(:,:),POINTER :: XSFSNW ! Blowing snow fluxes
- ! 1: Number deposition flux 2: Mass deposition flux 3: Streamwise saltation flux
- REAL, DIMENSION(:,:),POINTER :: XSNW_SUBL !Sublimation rate
- ! 1: Instantaneous number 2: Instantaneous mass 3: Accumulated Mass
- REAL, DIMENSION(:,:),POINTER :: XSNW_FTURB ! Snow surface turbulent flux
- ! 1: Instantaneous number 2: Instantaneous mass 3: Accumulated Mass
- REAL, DIMENSION(:,:),POINTER :: XSNW_FSED ! Snow surface sedimentation flux
- ! 1: Instantaneous number 2: Instantaneous mass 3: Accumulated Mass
- REAL, DIMENSION(:,:),POINTER :: XSNW_FNET ! Total surface net flux
- ! (salt+ susp)
- ! 1: Instantaneous number 2: Instantaneous mass 3: Accumulated Mass
- REAL, DIMENSION(:,:),POINTER :: XSNW_FSALT ! Saltation deposition/erosion flux
- ! 1: streamwise flux 2: Instantaneous 3: Accumulated
- REAL, DIMENSION(:,:),POINTER :: XSNW_CANO_RGA ! !Blowing snow radius at canopy levels (m)
- REAL, DIMENSION(:,:,:),POINTER :: XSNW_CANO_VAR ! !Blowing snow variables at canopy levels (_/m3)
-
-
-END TYPE BLOWSNW_t
-
-CONTAINS
-
-SUBROUTINE BLOWSNW_INIT(YBLOWSNW)
-TYPE(BLOWSNW_t), INTENT(INOUT) :: YBLOWSNW
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-IF (LHOOK) CALL DR_HOOK("MODD_BLOWSNW_N:BLOWSNW_INIT",0,ZHOOK_HANDLE)
- NULLIFY(YBLOWSNW%XSFSNW)
- NULLIFY(YBLOWSNW%XSNW_SUBL)
- NULLIFY(YBLOWSNW%XSNW_FTURB)
- NULLIFY(YBLOWSNW%XSNW_FSED)
- NULLIFY(YBLOWSNW%XSNW_FNET)
- NULLIFY(YBLOWSNW%XSNW_FSALT)
-IF (LHOOK) CALL DR_HOOK("MODD_BLOWSNW_N:BLOWSNW_INIT",1,ZHOOK_HANDLE)
-END SUBROUTINE BLOWSNW_INIT
-
-END MODULE MODD_BLOWSNW_n
diff --git a/src/SURFEX/mode_blowsnw_mblutl.f90 b/src/SURFEX/mode_blowsnw_mblutl.f90
deleted file mode 100644
index f730e1bc8..000000000
--- a/src/SURFEX/mode_blowsnw_mblutl.f90
+++ /dev/null
@@ -1,667 +0,0 @@
-!-----------------------------------------------------------------
-! ##################
- MODULE MODE_BLOWSNW_MBLUTL
-! ##################
-!
-!!**** *MODE_BLOWSNW_MBLUTL * - contains subroutines to determine:
-! - threshold velocity for snow erosion as a function of snow
-! pack properties,
-! - blown snow particles concentration in the saltation layer,
-! - surface turbulent flux of blown snow particles
-!!
-!! PURPOSE
-!! -------
-!
-!!** METHOD
-!! ------
-!!
-!! EXTERNAL
-!! --------
-!!
-!! IMPLICIT ARGUMENTS
-!! ------------------
-!!
-!! REFERENCE
-!! ---------
-!!
-!! AUTHOR
-!! ------
-!! V. Vionnet * Meteo France *
-!!
-!! MODIFICATIONS
-!! -------------
-!! Original 21/03/08
-!!
-!-----------------------------------------------------------------------------
-!
-!* 0. DECLARATIONS
-!
-!
-!-------------------------------------------------------------------------------
-CONTAINS
-
-
- FUNCTION SNOW_MBL_INDEX(PSNOWGRAN1, PSNOWGRAN2) &
-RESULT(PMBLINDEX)
-!
-!! PURPOSE
-!! -------
-! Calculate the mobility index as a function of snow grain type
-! Based on G&M,1998 : PROTEON, Centre d'Etude de la Neige
-!
-USE MODD_SNOW_METAMO, ONLY : XGRAN,XEPSI
-USE YOMHOOK ,ONLY : LHOOK, DR_HOOK
-USE PARKIND1 ,ONLY : JPRB
-!
-!
-!* 0.1 Declarations of arguments
-!
-REAL, INTENT(IN) :: PSNOWGRAN1, PSNOWGRAN2
-!
-REAL :: PMBLINDEX
-!
-!* 0.2 declarations of local variables
-!
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!
-!
-! 0.3 Initialization
-!
-IF (LHOOK) CALL DR_HOOK('MODE_BLOWSNW_MBLUTL:SNOW_MBL_INDEX',0,ZHOOK_HANDLE)
-!
-PMBLINDEX=0.
-!
-! 1. Computation
-!
-IF(PSNOWGRAN1<XEPSI)THEN
- PMBLINDEX=-0.75*PSNOWGRAN1/XGRAN-0.5*PSNOWGRAN2/ &
- XGRAN+0.5
-ELSE IF(PSNOWGRAN1>=0 .AND. PSNOWGRAN2>0) THEN
- PMBLINDEX=-0.583*PSNOWGRAN2*1000-0.833*PSNOWGRAN1/ &
- XGRAN+0.833
-ELSE
- PMBLINDEX=0.
-ENDIF
-
-PMBLINDEX=MAX(MIN(PMBLINDEX,1.),-0.9999)
-
-IF (LHOOK) CALL DR_HOOK('MODE_BLOWSNW_MBLUTL:SNOW_MBL_INDEX',1,ZHOOK_HANDLE)
-
-END FUNCTION SNOW_MBL_INDEX
-
-
- SUBROUTINE WIND_RFR_THR(PWIND_RFR_THR,PSNOWGRAN1, &
- PSNOWGRAN2)
-!
-!! PURPOSE
-!! -------
-!
-! Computation of threshold wind speed at level of reference (5
-! meter from G&M)
-!
-USE YOMHOOK ,ONLY : LHOOK, DR_HOOK
-USE PARKIND1 ,ONLY : JPRB
- IMPLICIT NONE
-!
-! 0.1 declarations of arguments
-!
-!
- REAL, INTENT(IN) :: PSNOWGRAN1, PSNOWGRAN2
-
- REAL, INTENT(INOUT) ::PWIND_RFR_THR
-!
-! 0.2 declaration of local variables
-!
- REAL :: ZMBLINDEX
- REAL(KIND=JPRB) :: ZHOOK_HANDLE
-
-!
-! 1. mobility index
-!
-IF (LHOOK) CALL DR_HOOK('MODE_BLOWSNW_MBLUTL:WIND_RFR_THR',0,ZHOOK_HANDLE)
-!
- ZMBLINDEX= SNOW_MBL_INDEX(PSNOWGRAN1,PSNOWGRAN2)
-!
-! 2. threshold 5-m wind speed
-!
- PWIND_RFR_THR=(log(2.868)-log(1+ZMBLINDEX))/0.085
-!
-IF (LHOOK) CALL DR_HOOK('MODE_BLOWSNW_MBLUTL:WIND_RFR_THR',1,ZHOOK_HANDLE)
-!
- END SUBROUTINE WIND_RFR_THR
-
- SUBROUTINE WIND_FRC_THR(PWIND_RFR_THR,PZ0EFF, &
- PWIND_FRC_THR)
-!
-!! PURPOSE
-!! -------
-!
-! Compute threshold friction velocity from the threshold
-! velocity at a reference height of 5 m. We assume a log profile
-! in the SBL with z0 the same as used to compute CDSNOW
-!
-!
- USE MODD_CSTS,ONLY : XKARMAN
-!
-USE YOMHOOK ,ONLY : LHOOK, DR_HOOK
-USE PARKIND1 ,ONLY : JPRB
-!
- IMPLICIT NONE
-!
-! 0.1 declarations of arguments
-!
- REAL, INTENT(IN) :: PZ0EFF
- REAL, INTENT(IN) :: PWIND_RFR_THR
-!
- REAL, INTENT(OUT) :: PWIND_FRC_THR
-!
-! 0.2 declaration of local variables
- REAL :: ZCD
- REAL(KIND=JPRB) :: ZHOOK_HANDLE
-
-! !
-! 0.3 Initialization
-!
-IF (LHOOK) CALL DR_HOOK('MODE_BLOWSNW_MBLUTL:WIND_FRC_THR',0,ZHOOK_HANDLE)
-!
- ZCD = (XKARMAN/LOG(5/PZ0EFF))**2
-
-! 1. Threshold friction velocity
-!
- PWIND_FRC_THR=SQRT(ZCD)*PWIND_RFR_THR
-!
-IF (LHOOK) CALL DR_HOOK('MODE_BLOWSNW_MBLUTL:WIND_FRC_THR',1,ZHOOK_HANDLE)
-!
- END SUBROUTINE WIND_FRC_THR
-
- SUBROUTINE CONC_SALT(PCONC_SALT,PWIND_FRC_SALT,PWIND_FRC_THR, &
- PRHOA,PQSALT,PSNOWLIQ,PSNOWHIST)
-!
-!! PURPOSE
-!! -------
-!
-! Saltation layer parameterization : computes
-! - snow particle concentration (kg/m3)
-! - streamwise saltation flux (kg/m/s)
-!
-!
-USE MODD_CSTS, ONLY : XG,XRHOLI
-USE MODD_BLOWSNW_SURF
-USE YOMHOOK ,ONLY : LHOOK, DR_HOOK
-USE PARKIND1 ,ONLY : JPRB
-
- IMPLICIT NONE
-!
-! 0.1 declarations of arguments
-!
- REAL, INTENT(IN) :: PWIND_FRC_SALT,PRHOA, &
- PWIND_FRC_THR
-
- REAL, INTENT(IN) :: PSNOWLIQ,PSNOWHIST
- REAL, INTENT(INOUT) :: PCONC_SALT
- REAL, INTENT(OUT) :: PQSALT
-!
-! 0.2 declaration of local variables
-!
- REAL :: ZSALT_FLUX
- REAL :: ZHSALT
- REAL :: ZUPART
- REAL(KIND=JPRB) :: ZHOOK_HANDLE
-
-!
-! 0.3 Initialization
-!
-IF (LHOOK) CALL DR_HOOK('MODE_BLOWSNW_MBLUTL:CONC_SALT',0,ZHOOK_HANDLE)
-!
-! 1. Streamiwise saltation flux : horizontal saltation, transport rate under
-! equilibrium conditions
-
-!
- IF(PWIND_FRC_SALT>=PWIND_FRC_THR) THEN ! Transport in saltation
- IF(CSNOW_SALT=='POME' .OR. CSNOW_SALT=='TPOM') THEN
-!
-! Formulation of Pomeroy and Gray, 1990
-!
- ZSALT_FLUX=0.68*PRHOA*PWIND_FRC_THR* &
- (PWIND_FRC_SALT**2-PWIND_FRC_THR**2)/ &
- (PWIND_FRC_SALT*XG)
-
- ELSE IF(CSNOW_SALT=='SORE'.OR. CSNOW_SALT=='TSOR') THEN
-!
-! Formulation of Sorensen (2004) (initially developped for sand)
-! with coefficients adapted for snow based on the wind-tunnel
-! measurement of Nishimura and Hunt (2000). More details in the PhD
-! thesis of Vionnet (2012)
-!
- ZSALT_FLUX=PRHOA/XG*PWIND_FRC_SALT**3.* &
- (1.-PWIND_FRC_THR**2./PWIND_FRC_SALT**2.)* &
- (2.6+2.*PWIND_FRC_THR/PWIND_FRC_SALT+2.5* &
- PWIND_FRC_THR**2./PWIND_FRC_SALT**2.)
- ENDIF
- ELSE ! No transport
- ZSALT_FLUX=0.
- ENDIF
-!
-! Saltation does not occur when :
-! - surface layer is wet
-! - historical variable is higher than 1 : crust of non-transportable snow
- IF(PSNOWLIQ>0 .OR. PSNOWHIST>0) THEN
- ZSALT_FLUX=0
- ENDIF
-!
-! Store streamise saltation flux for future advection
- PQSALT=ZSALT_FLUX
-!
-! 2. Saltation height (Pomeroy and Male, 1992)
-!
- ZHSALT = 0.08436*PWIND_FRC_SALT**1.27
-!
-! 3. Horizontalparticle velocity Up=2,8.U*t
-! according to Pomeroy and Gray (1900)
-!
- ZUPART = 2.8*PWIND_FRC_THR
-!
-! 4. Saltating particle concentration
-!
- IF(ZSALT_FLUX>0.) THEN
- PCONC_SALT=ZSALT_FLUX/(ZUPART*ZHSALT)
- ELSE
- PCONC_SALT=0.
- ENDIF
-
- IF(CSNOW_SALT=='TPOM' .OR. CSNOW_SALT=='TSOR') THEN
-! Compute concentration at the top of the saltation layer assuming
-! an exponential flux profile in the saltation layer (cf Doorschot
-! et al., 2004 and Nishimura and Hunt, 2000).
- IF(ZSALT_FLUX>0.) THEN
- PCONC_SALT = ZSALT_FLUX*0.45*XG/PWIND_FRC_SALT**2 &
- *EXP(-0.45*ZHSALT*XG/PWIND_FRC_SALT**2) &
- /(ZUPART)
- ENDIF
- ENDIF
-
- IF(LSNOW_SALT) THEN
-! Imposed concentration in the saltation layer (only for sensitivity
-! studies)
- PCONC_SALT= XCONC_SALT
- PQSALT = XCONC_SALT*ZUPART*ZHSALT
- END IF
-!
-IF (LHOOK) CALL DR_HOOK('MODE_BLOWSNW_MBLUTL:CONC_SALT',1,ZHOOK_HANDLE)
-!
- END SUBROUTINE CONC_SALT
-
-
- SUBROUTINE SNOW_FLUX(PCONC_SALT,PVMOD,PCONC_SURF, &
- PCDSNOW,PSNOW_FLUX)
-! Compute snow turbulent vertical flux as in Gallee et al (2001)
-! Flux=u*q*=Cd.V.(Q(Surf)-Q(Salt)) (kg/(m2.s))
-! Cd : drag coefficient of momentum
-! V wind speed in the lowest level of the atmo. model (m/s)
-! Q(Surf) blown snow concentration in the lowest level of the atmo. model (kg/m3)
-! Q(Salt) blown snow concentration in the saltation layer (kg/m3)
-! Only positive flux : from the saltation layer to the model lowest level. Deposition is computed later
-!
-! NOTE: the lowest level of the atmo. model can be Canopy or directly Meso-NH if Canopy is
-! not activated
-!
-USE MODD_CSTS, ONLY : XG,XRHOLI, XPI
-USE MODD_BLOWSNW_SURF
-USE YOMHOOK ,ONLY : LHOOK, DR_HOOK
-USE PARKIND1 ,ONLY : JPRB
-
- IMPLICIT NONE
-!
-! 0.1 declarations of arguments
-!
- REAL, INTENT(IN) :: PCONC_SALT,PCDSNOW,PVMOD
- REAL, DIMENSION(:),INTENT(IN) :: PCONC_SURF
-!
- REAL, DIMENSION(:),INTENT(OUT) :: PSNOW_FLUX
-!
-!
-! 0.2 declaration of local variables
-!
- REAL :: ZNUM_SALT ! Number particle concentration in the saltation
-! layer (#.m^{-3})
- REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!
-! 1. Initialisation following the gamma distribution
-!
-IF (LHOOK) CALL DR_HOOK('MODE_BLOWSNW_MBLUTL:SNOW_FLUX',0,ZHOOK_HANDLE)
-!
- ZNUM_SALT = 3*PCONC_SALT/(4*XPI*XRHOLI*(XEMIALPHA_SNW+2.)*(XEMIALPHA_SNW+1.)* &
- XEMIALPHA_SNW*(XEMIRADIUS_SNW/XEMIALPHA_SNW)**3)
-!
-! 2. Turbulent fluxes
-!
-! Number flux
-!
- PSNOW_FLUX(1)=MAX(0.,-PCDSNOW*PVMOD*(PCONC_SURF(1)- &
- ZNUM_SALT))
-! Mass flux
- PSNOW_FLUX(2)=MAX(0.,-PCDSNOW*PVMOD*(PCONC_SURF(2)- &
- PCONC_SALT))
-!
-IF (LHOOK) CALL DR_HOOK('MODE_BLOWSNW_MBLUTL:SNOW_FLUX',1,ZHOOK_HANDLE)
-!
-END SUBROUTINE SNOW_FLUX
-
-
- SUBROUTINE SURFACE_RI_PART(PCONC_SALT, PCONC_SURF,PRHOA, &
- PZREF, PUREF, PDIRCOSZW, PVMOD, PRI_PART )
-! ######################################################################
-!
-!!**** *SURFACE_RI*
-!!
-!! PURPOSE
-!! -------
-!
-! Computes the component of richardson number near the ground associated
-! with particle-induced buoyancy.
-!
-!
-!!** METHOD
-!! ------
-!
-! Based on surface_ri.f90
-!
-!
-!! EXTERNAL
-!! --------
-!!
-!!
-!! IMPLICIT ARGUMENTS
-!! ------------------
-!!
-!! MODD_CST
-!! MODD_GROUND_PAR
-!!
-!!
-!! REFERENCE
-!! ---------
-!!
-!!
-!! AUTHOR
-!! ------
-!!
-!! V. Vionnet * Meteo-France *
-!!
-!! MODIFICATIONS
-!! -------------
-!! Original 04/10/10
-!-------------------------------------------------------------------------------
-!
-!* 0. DECLARATIONS
-! ------------
-!
-USE MODD_CSTS,ONLY : XRV, XRD, XG
-USE MODI_WIND_THRESHOLD
-USE YOMHOOK ,ONLY : LHOOK, DR_HOOK
-USE PARKIND1 ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!* 0.1 declarations of arguments
-!
-!
-REAL, INTENT(IN) :: PCONC_SALT ! Concentration of blowing snow particle
-! in the saltation layer (kg/m3)
-REAL, INTENT(IN) :: PCONC_SURF ! Concentration of blowing snow particle
-! at the 1st atmospheric level (kg/m3)
-REAL, INTENT(IN) :: PVMOD ! module of the horizontal wind
-REAL, INTENT(IN) :: PRHOA ! air density (kg/m3)
-!
-REAL, INTENT(IN) :: PZREF ! reference height of the first
- ! atmospheric level
-REAL, INTENT(IN) :: PUREF ! reference height of the wind
-! ! NOTE this is different from ZZREF
-! ! ONLY in stand-alone/forced mode,
-! ! NOT when coupled to a model (MesoNH)
-REAL, INTENT(IN) :: PDIRCOSZW! Cosine of the angle between
-! ! the normal to the surface and
-! ! the vertical
-!
-REAL, INTENT(OUT) :: PRI_PART ! "Particle" Richardson number
-!
-!* 0.2 declarations of local variables
-!
-REAL :: ZVMOD
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!-------------------------------------------------------------------------------
-!
-! 1. Richardson number
-! -----------------
-!
-IF (LHOOK) CALL DR_HOOK('MODE_BLOWSNW_MBLUTL:SURFACE_RI_PART',0,ZHOOK_HANDLE)
-!
-ZVMOD = PVMOD
-!
-!
- ! Richardson's number
-PRI_PART = XG * PDIRCOSZW/PRHOA * PUREF * PUREF &
- * (PCONC_SALT-PCONC_SURF) &
- / (ZVMOD*ZVMOD) /PZREF
-!-------------------------------------------------------------------------------
-!
-IF (LHOOK) CALL DR_HOOK('MODE_BLOWSNW_MBLUTL:SURFACE_RI_PART',1,ZHOOK_HANDLE)
-!
-END SUBROUTINE SURFACE_RI_PART
-
-! #################################################################
- SUBROUTINE SURFACE_CD_PART(PRI, PZREF, PUREF, PZ0EFF, PZ0H, &
- PCD, PCDN)
-! #################################################################
-!
-!!**** *SURFACE_CD_PART*
-!!
-!! PURPOSE
-!! -------
-!
-! Computes the drag coefficients for momentum near the ground
-! including particle-induced buoyancy.
-!
-!!** METHOD
-!! ------
-!
-! 1 and 2 : computation of relative humidity near the ground
-!
-! 3 : richardson number including particle-induced buoyancy
-!
-! 4 : the aerodynamical resistance for heat transfers is deduced
-!
-! 5 : the drag coefficient for momentum ZCD is computed
-! including particle-induced buoyancy
-!
-!! REFERENCE
-!! ---------
-!!
-!!
-!! AUTHOR
-!! ------
-!!
-!! V. Vionnet * Meteo-France *
-!! Based on surface_cd.f90 by V. Masson
-!!
-!! MODIFICATIONS
-!! -------------
-!! Original 04/10/10
-
-!-------------------------------------------------------------------------------
-!
-!* 0. DECLARATIONS
-! ------------
-!
-USE MODD_CSTS,ONLY : XKARMAN
-!
-USE MODE_THERMOS
-!
-USE YOMHOOK ,ONLY : LHOOK, DR_HOOK
-USE PARKIND1 ,ONLY : JPRB
-!
-IMPLICIT NONE
-!
-!* 0.1 declarations of arguments
-!
-!
-REAL, INTENT(IN) :: PRI ! Richardson number
-REAL, INTENT(IN) :: PZREF ! reference height of the first
- ! atmospheric level
-REAL, INTENT(IN) :: PUREF ! reference height of the wind
-! ! NOTE this is different from ZZREF
-! ! ONLY in stand-alone/forced mode,
-! ! NOT when coupled to a model (MesoNH)
-REAL, INTENT(IN) :: PZ0EFF ! roughness length for momentum
- ! with subgrid-scale orography
-REAL, INTENT(IN) :: PZ0H ! roughness length for heat
-!
-REAL, INTENT(OUT) :: PCD ! drag coefficient for momentum
-REAL, INTENT(OUT) :: PCDN ! neutral drag coefficient for momentum
-!
-!* 0.2 declarations of local variables
-!
-!
-REAL :: ZZ0EFF, ZZ0H, ZMU, &
- ZCMSTAR, ZPM, ZCM, ZFM
-INTEGER :: JJ
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-
-!-------------------------------------------------------------------------------
-!
-!* 1. Drag coefficient for momentum transfers
-! ---------------------------------------
-!
-IF (LHOOK) CALL DR_HOOK('MODE_BLOWSNW_MBLUTL:SURFACE_CD_PART',0,ZHOOK_HANDLE)
-!
- ZZ0EFF = MIN(PZ0EFF,PUREF*0.5)
- ZZ0H = MIN(ZZ0EFF,PZ0H)
-!
- ZMU = LOG( MIN(ZZ0EFF/ZZ0H,200.) )
-!
- PCDN = (XKARMAN/LOG(PUREF/ZZ0EFF))**2
-
- ZCMSTAR = CMSTAR(ZMU)
- ZPM = PM(ZMU)
-!
- ZCM = 10.*ZCMSTAR*PCDN*( PUREF/ZZ0EFF )**ZPM
-!
- IF ( PRI > 0.0 ) THEN
- ZFM = 1. + 10.*PRI / SQRT( 1.+5.*PRI )
- ZFM = 1. / ZFM
- ELSE
- ZFM = 1. - 10.*PRI / ( 1.+ZCM*SQRT(-PRI) )
- ENDIF
-!
- PCD = PCDN*ZFM
-!
-!
-!-------------------------------------------------------------------------------
-!
-!
-!
-IF (LHOOK) CALL DR_HOOK('MODE_BLOWSNW_MBLUTL:SURFACE_CD_PART',1,ZHOOK_HANDLE)
-
-!
-CONTAINS
-!
-! functions used in the calculation
-! the terms Cm
-!
- FUNCTION CMSTAR(X)
- REAL, INTENT(IN) :: X
- REAL :: CMSTAR
- !
- CMSTAR = 6.8741 + 2.6933*X - 0.3601*X*X + 0.0154*X*X*X
- !
- END FUNCTION CMSTAR
- !
- !
- FUNCTION PM(X)
- REAL, INTENT(IN) :: X
- REAL :: PM
- !
- PM = 0.5233 - 0.0815*X + 0.0135*X*X - 0.0010*X*X*X
- !
- END FUNCTION PM
-
-!
-!-------------------------------------------------------------------------------
-
-
-!
-END SUBROUTINE SURFACE_CD_PART
-
-SUBROUTINE SOLVE_CD(PUSTAR,PVMOD,PRHOA,PCONC_SURF,PZREF,PUREF,PDIRCOSZW, &
- PWIND_FRC_THR,PRI,PZ0EFF,HSNOWRES,PZ0H,PRES,PRI_TOT, &
- PCDSNOW,PSNOWLIQ,PSNOWHIST)
-
-!
-! Compute the value of the function used in the Newton's iterative algorithm.
-!
-
-
-USE MODD_BLOWSNW_SURF
-USE MODD_SNOW_PAR, ONLY : X_RI_MAX
-USE YOMHOOK ,ONLY : LHOOK, DR_HOOK
-USE PARKIND1 ,ONLY : JPRB
-
-IMPLICIT NONE
-!
-!* 0.1 declarations of arguments
-!
-REAL, INTENT(IN) :: PUSTAR
-REAL, INTENT(IN) :: PSNOWLIQ,PSNOWHIST
-REAL, INTENT(IN) :: PVMOD,PRHOA,PCONC_SURF,PZREF,PUREF
-REAL, INTENT(IN) :: PDIRCOSZW,PWIND_FRC_THR,PRI,PZ0EFF,PZ0H
-REAL, INTENT(OUT) :: PCDSNOW, PRI_TOT,PRES
-CHARACTER(LEN=*), INTENT(IN) :: HSNOWRES
-!
-!* 0.2 declarations of local variables
-!
-REAL :: ZCONC_SALT,ZRI_PART,ZCD, ZCDN, ZZ0_TEMP,ZQSALT
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-
-PRI_TOT = PRI
-ZZ0_TEMP = PZ0EFF
-
-IF (LHOOK) CALL DR_HOOK('MODE_BLOWSNW_MBLUTL:SOLVE_CD',0,ZHOOK_HANDLE)
-
-IF(XSNOW_BUOYANCY==1) THEN
-!
-! Blown snow particles concentration (kg/m3) in the saltation
-! layer
-!
-CALL CONC_SALT(ZCONC_SALT,PUSTAR,PWIND_FRC_THR,PRHOA,ZQSALT ,&
- PSNOWLIQ,PSNOWHIST)
-
-CALL SURFACE_RI_PART(ZCONC_SALT,PCONC_SURF, PRHOA, &
- PZREF, PUREF, PDIRCOSZW, PVMOD, ZRI_PART )
-
-!
-PRI_TOT = PRI + ZRI_PART
-!
-! Simple adaptation of method by Martin and Lejeune (1998)
-! to apply a lower limit to turbulent transfer coef
-! by defining a maximum Richarson number for stable
-! conditions:
-!
-IF(HSNOWRES=='RIL') PRI_TOT = MIN(X_RI_MAX, PRI_TOT)
-END IF
-
-IF(XSNOW_ROUGHNESS==1 .AND. PUSTAR > PWIND_FRC_THR ) THEN
- ! Increase in surface rougnhess due to saltation following
- ! Dover (1993) z0_s = z0_ns*(u*/uth*)�
- ! Limit increase to value observed over snow surface
- ZZ0_TEMP = MIN(PZ0EFF * (PUSTAR/PWIND_FRC_THR)**2,0.01)
-END IF
-
-CALL SURFACE_CD_PART(PRI_TOT, PZREF, PUREF, ZZ0_TEMP, PZ0H, PCDSNOW, ZCDN)
-
-PRES = PUSTAR - SQRT(PCDSNOW)* PVMOD
-
-IF (LHOOK) CALL DR_HOOK('MODE_BLOWSNW_MBLUTL:SOLVE_CD',1,ZHOOK_HANDLE)
-
-
-END SUBROUTINE SOLVE_CD
-
-END MODULE MODE_BLOWSNW_MBLUTL
diff --git a/src/SURFEX/mode_blowsnw_sedim_lkt1d.f90 b/src/SURFEX/mode_blowsnw_sedim_lkt1d.f90
deleted file mode 100644
index 095aa45ce..000000000
--- a/src/SURFEX/mode_blowsnw_sedim_lkt1d.f90
+++ /dev/null
@@ -1,1297 +0,0 @@
-
-!-----------------------------------------------------------------
-!--------------- special set of characters for RCS information
-!-----------------------------------------------------------------
-! $Source: /sxturb1/data1/mesonh/CODEMNH/mode_dustopt.f90,v $ $Revision: 1.1.2.1 $
-! masdev4_7 BUG1 2007/06/29 12:06:27
-!-----------------------------------------------------------------
-! ###################
-MODULE MODE_BLOWSNW_SEDIM_LKT1D
-! ###################
-!
-!!
-!! PURPOSE
-!! -------
-!!
-!! Sedimentation of snow particles is computed according to Carrier's drag coofficient.
-!! To reduce computational time; look-up tables are used. They depend on the
-!! average radius and the pressure (interpolation)
-!! AUTHOR
-!! ------
-!! Vincent Vionnet (CNRM/GMEI)
-!!
-!! MODIFICATIONS
-!! -------------
-!!
- IMPLICIT NONE
- PUBLIC
-
-CONTAINS
-
-
-! Sedimentation of snow particles is computed according to Carrier's drag coofficient.
-! To reduce computational time; look-up tables are used. They depend on the
-! average radius and the pressure (interpolation)
-
-SUBROUTINE BLOWSNW_SEDIM_LKT1D(PRG, & ! Mean radius (m)
- PABST,& ! Pressure
- PVGK) ! Number and mass average settling velocity
-
-USE MODD_BLOWSNW_SEDIM_LKT1D
-
-USE YOMHOOK ,ONLY : LHOOK, DR_HOOK
-USE PARKIND1 ,ONLY : JPRB
-!
-!* 0.1 declarations of arguments
-!
-REAL, DIMENSION(:,:), INTENT(IN) :: PRG
-REAL, DIMENSION(:,:), INTENT(IN) :: PABST
-REAL, DIMENSION(:,:,:), INTENT(OUT) :: PVGK
-!
-!* 0.2 declaration of local variables
-!
-REAL, DIMENSION(SIZE(PRG,1),SIZE(PRG,2)) :: ZRG, ZPRESSURE
-REAL :: ZIND,ZFACT_RG,ZFACT_PRESSURE
-REAL :: ZV1_MASS,ZV1_NUMB,ZV2_MASS,ZV2_NUMB
-INTEGER :: RG_IDX,PR_IDX
-INTEGER :: JI,JJ,JK,II !Loop counter
-REAL :: M2UM ! Convert m to um
-
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-
-IF (LHOOK) CALL DR_HOOK('BLOWSNW_SEDIM_LKT1D',0,ZHOOK_HANDLE)
-
-M2UM = 1000000
-
-
-!Remove unphysical values for rg and pressure
-ZRG(:,:) = MIN( MAX(XRADIUS_LKT1D_MIN,PRG(:,:)*M2UM), XRADIUS_LKT1D_MAX)
-ZPRESSURE(:,:) = MIN( MAX(XPRESSURE_LKT1D_MIN,PABST(:,:)), XPRESSURE_LKT1D_MAX-1.)
-ZFACT_RG = DBLE(NMAX_RADIUS_LKT1D-1)/(XRADIUS_LKT1D_MAX-XRADIUS_LKT1D_MIN)
-ZFACT_PRESSURE = DBLE(NMAX_PRESSURE_LKT1D-1)/(XPRESSURE_LKT1D_MAX-XPRESSURE_LKT1D_MIN)
-
-DO JI=1,SIZE(PRG,1)
- DO JK=1,SIZE(PRG,2)
- !Get the correct indexes for the look up tables
- RG_IDX = INT((ZRG(JI,JK)-XRADIUS_LKT1D_MIN)*ZFACT_RG + 1.5)
- ZIND = (ZPRESSURE(JI,JK)-XPRESSURE_LKT1D_MIN)*ZFACT_PRESSURE + 1.
- PR_IDX = INT(ZIND)
- ZV1_NUMB = XNUMB_SPEED_LKT1D(RG_IDX,PR_IDX)
- ZV2_NUMB = XNUMB_SPEED_LKT1D(RG_IDX,PR_IDX+1)
- ZV1_MASS = XMASS_SPEED_LKT1D(RG_IDX,PR_IDX)
- ZV2_MASS = XMASS_SPEED_LKT1D(RG_IDX,PR_IDX+1)
- ! Linear interpolation to get correct values
- PVGK(JI,JK,1) = ZV1_NUMB + (ZIND - PR_IDX) * (ZV2_NUMB-ZV1_NUMB)
- PVGK(JI,JK,2) = ZV1_MASS + (ZIND - PR_IDX) * (ZV2_MASS-ZV1_MASS)
-! PVGK(JI,JK,3) = PVGK(JI,JK,2)
- ENDDO
-ENDDO
-
-IF (LHOOK) CALL DR_HOOK('BLOWSNW_SEDIM_LKT1D',1,ZHOOK_HANDLE)
-
-END SUBROUTINE BLOWSNW_SEDIM_LKT1D
-
-SUBROUTINE BLOWSNW_SEDIM_LKT1D_SET
-
-USE MODD_BLOWSNW_SEDIM_LKT1D
-
-USE YOMHOOK ,ONLY : LHOOK, DR_HOOK
-USE PARKIND1 ,ONLY : JPRB
-!
-!* 0.2 declaration of local variables
-!
-LOGICAL :: OPIEKTUK
-
-REAL(KIND=JPRB) :: ZHOOK_HANDLE
-
-IF (LHOOK) CALL DR_HOOK('BLOWSNW_SEDIM_LKT1D_SET',0,ZHOOK_HANDLE)
-
-OPIEKTUK=.FALSE.
-
-IF(OPIEKTUK) THEN
-
-XNUMB_SPEED_LKT1D(1,1:4) = (/0.015056,0.013491,0.012659,0.012142/)
-XMASS_SPEED_LKT1D(1,1:4) = (/0.026814,0.025878,0.025363,0.025028/)
-
-XNUMB_SPEED_LKT1D(2,1:4) = (/0.018595,0.017278,0.016572,0.016129/)
-XMASS_SPEED_LKT1D(2,1:4) = (/0.036750,0.035901,0.035405,0.035061/)
-
-XNUMB_SPEED_LKT1D(3,1:4) = (/0.023150,0.021999,0.021372,0.020972/)
-XMASS_SPEED_LKT1D(3,1:4) = (/0.048615,0.047758,0.047209,0.046796/)
-
-XNUMB_SPEED_LKT1D(4,1:4) = (/0.028622,0.027576,0.026992,0.026605/)
-XMASS_SPEED_LKT1D(4,1:4) = (/0.062287,0.061322,0.060637,0.060081/)
-
-XNUMB_SPEED_LKT1D(5,1:4) = (/0.034946,0.033956,0.033378,0.032979/)
-XMASS_SPEED_LKT1D(5,1:4) = (/0.077657,0.076464,0.075545,0.074760/)
-
-XNUMB_SPEED_LKT1D(6,1:4) = (/0.042073,0.041092,0.040487,0.040046/)
-XMASS_SPEED_LKT1D(6,1:4) = (/0.094610,0.093053,0.091783,0.090669/)
-
-XNUMB_SPEED_LKT1D(7,1:4) = (/0.049963,0.048943,0.048273,0.047759/)
-XMASS_SPEED_LKT1D(7,1:4) = (/0.113031,0.110950,0.109200,0.107647/)
-
-XNUMB_SPEED_LKT1D(8,1:4) = (/0.058579,0.057469,0.056693,0.056071/)
-XMASS_SPEED_LKT1D(8,1:4) = (/0.132797,0.130017,0.127645,0.125534/)
-
-XNUMB_SPEED_LKT1D(9,1:4) = (/0.067886,0.066631,0.065704,0.064936/)
-XMASS_SPEED_LKT1D(9,1:4) = (/0.153787,0.150115,0.146970,0.144180/)
-
-XNUMB_SPEED_LKT1D(10,1:4) = (/0.077848,0.076388,0.075262,0.074306/)
-XMASS_SPEED_LKT1D(10,1:4) = (/0.175877,0.171111,0.167037,0.163443/)
-
-XNUMB_SPEED_LKT1D(11,1:4) = (/0.088431,0.086701,0.085324,0.084137/)
-XMASS_SPEED_LKT1D(11,1:4) = (/0.198946,0.192874,0.187714,0.183195/)
-
-XNUMB_SPEED_LKT1D(12,1:4) = (/0.099598,0.097531,0.095847,0.094384/)
-XMASS_SPEED_LKT1D(12,1:4) = (/0.222876,0.215284,0.208884,0.203321/)
-
-XNUMB_SPEED_LKT1D(13,1:4) = (/0.111314,0.108836,0.106789,0.105003/)
-XMASS_SPEED_LKT1D(13,1:4) = (/0.247555,0.238228,0.230436,0.223719/)
-
-XNUMB_SPEED_LKT1D(14,1:4) = (/0.123545,0.120580,0.118111,0.115954/)
-XMASS_SPEED_LKT1D(14,1:4) = (/0.272876,0.261603,0.252275,0.244297/)
-
-XNUMB_SPEED_LKT1D(15,1:4) = (/0.136256,0.132725,0.129774,0.127197/)
-XMASS_SPEED_LKT1D(15,1:4) = (/0.298740,0.285316,0.274314,0.264979/)
-
-XNUMB_SPEED_LKT1D(16,1:4) = (/0.149412,0.145235,0.141740,0.138697/)
-XMASS_SPEED_LKT1D(16,1:4) = (/0.325053,0.309281,0.296478,0.285697/)
-
-XNUMB_SPEED_LKT1D(17,1:4) = (/0.162981,0.158074,0.153975,0.150418/)
-XMASS_SPEED_LKT1D(17,1:4) = (/0.351731,0.333423,0.318701,0.306394/)
-
-XNUMB_SPEED_LKT1D(18,1:4) = (/0.176929,0.171209,0.166446,0.162329/)
-XMASS_SPEED_LKT1D(18,1:4) = (/0.378696,0.357675,0.340926,0.327022/)
-
-XNUMB_SPEED_LKT1D(19,1:4) = (/0.191225,0.184610,0.179123,0.174402/)
-XMASS_SPEED_LKT1D(19,1:4) = (/0.405877,0.381980,0.363105,0.347539/)
-
-XNUMB_SPEED_LKT1D(20,1:4) = (/0.205840,0.198246,0.191977,0.186607/)
-XMASS_SPEED_LKT1D(20,1:4) = (/0.433212,0.406284,0.385195,0.367912/)
-
-XNUMB_SPEED_LKT1D(21,1:4) = (/0.220743,0.212088,0.204981,0.198922/)
-XMASS_SPEED_LKT1D(21,1:4) = (/0.460642,0.430545,0.407161,0.388113/)
-
-XNUMB_SPEED_LKT1D(22,1:4) = (/0.235908,0.226112,0.218111,0.211323/)
-XMASS_SPEED_LKT1D(22,1:4) = (/0.488118,0.454722,0.428975,0.408121/)
-
-XNUMB_SPEED_LKT1D(23,1:4) = (/0.251307,0.240291,0.231345,0.223789/)
-XMASS_SPEED_LKT1D(23,1:4) = (/0.515594,0.478783,0.450612,0.427917/)
-
-XNUMB_SPEED_LKT1D(24,1:4) = (/0.266917,0.254604,0.244661,0.236303/)
-XMASS_SPEED_LKT1D(24,1:4) = (/0.543030,0.502700,0.472052,0.447487/)
-
-XNUMB_SPEED_LKT1D(25,1:4) = (/0.282714,0.269029,0.258041,0.248846/)
-XMASS_SPEED_LKT1D(25,1:4) = (/0.570392,0.526448,0.493278,0.466821/)
-
-XNUMB_SPEED_LKT1D(26,1:4) = (/0.298675,0.283547,0.271468,0.261404/)
-XMASS_SPEED_LKT1D(26,1:4) = (/0.597648,0.550009,0.514278,0.485909/)
-
-XNUMB_SPEED_LKT1D(27,1:4) = (/0.314780,0.298139,0.284925,0.273963/)
-XMASS_SPEED_LKT1D(27,1:4) = (/0.624772,0.573364,0.535041,0.504748/)
-
-XNUMB_SPEED_LKT1D(28,1:4) = (/0.331009,0.312788,0.298398,0.286511/)
-XMASS_SPEED_LKT1D(28,1:4) = (/0.651741,0.596502,0.555561,0.523332/)
-
-XNUMB_SPEED_LKT1D(29,1:4) = (/0.347343,0.327479,0.311875,0.299037/)
-XMASS_SPEED_LKT1D(29,1:4) = (/0.678535,0.619410,0.575831,0.541661/)
-
-XNUMB_SPEED_LKT1D(30,1:4) = (/0.363766,0.342198,0.325343,0.311530/)
-XMASS_SPEED_LKT1D(30,1:4) = (/0.705138,0.642081,0.595848,0.559734/)
-
-XNUMB_SPEED_LKT1D(31,1:4) = (/0.380262,0.356932,0.338792,0.323983/)
-XMASS_SPEED_LKT1D(31,1:4) = (/0.731535,0.664507,0.615609,0.577551/)
-
-XNUMB_SPEED_LKT1D(32,1:4) = (/0.396817,0.371670,0.352213,0.336387/)
-XMASS_SPEED_LKT1D(32,1:4) = (/0.757714,0.686683,0.635115,0.595114/)
-
-XNUMB_SPEED_LKT1D(33,1:4) = (/0.413415,0.386400,0.365596,0.348735/)
-XMASS_SPEED_LKT1D(33,1:4) = (/0.783665,0.708606,0.654364,0.612425/)
-
-XNUMB_SPEED_LKT1D(34,1:4) = (/0.430045,0.401112,0.378935,0.361022/)
-XMASS_SPEED_LKT1D(34,1:4) = (/0.809380,0.730275,0.673358,0.629488/)
-
-XNUMB_SPEED_LKT1D(35,1:4) = (/0.446695,0.415799,0.392222,0.373242/)
-XMASS_SPEED_LKT1D(35,1:4) = (/0.834854,0.751687,0.692100,0.646306/)
-
-XNUMB_SPEED_LKT1D(36,1:4) = (/0.463354,0.430450,0.405452,0.385391/)
-XMASS_SPEED_LKT1D(36,1:4) = (/0.860079,0.772843,0.710590,0.662883/)
-
-XNUMB_SPEED_LKT1D(37,1:4) = (/0.480012,0.445061,0.418618,0.397465/)
-XMASS_SPEED_LKT1D(37,1:4) = (/0.885054,0.793743,0.728833,0.679223/)
-
-XNUMB_SPEED_LKT1D(38,1:4) = (/0.496660,0.459623,0.431717,0.409459/)
-XMASS_SPEED_LKT1D(38,1:4) = (/0.909775,0.814389,0.746832,0.695331/)
-
-XNUMB_SPEED_LKT1D(39,1:4) = (/0.513288,0.474131,0.444744,0.421372/)
-XMASS_SPEED_LKT1D(39,1:4) = (/0.934241,0.834783,0.764591,0.711212/)
-
-XNUMB_SPEED_LKT1D(40,1:4) = (/0.529891,0.488579,0.457695,0.433201/)
-XMASS_SPEED_LKT1D(40,1:4) = (/0.958451,0.854928,0.782114,0.726871/)
-
-XNUMB_SPEED_LKT1D(41,1:4) = (/0.546459,0.502964,0.470567,0.444943/)
-XMASS_SPEED_LKT1D(41,1:4) = (/0.982404,0.874825,0.799404,0.742312/)
-
-XNUMB_SPEED_LKT1D(42,1:4) = (/0.562987,0.517279,0.483358,0.456598/)
-XMASS_SPEED_LKT1D(42,1:4) = (/1.006102,0.894479,0.816467,0.757540/)
-
-XNUMB_SPEED_LKT1D(43,1:4) = (/0.579469,0.531523,0.496065,0.468163/)
-XMASS_SPEED_LKT1D(43,1:4) = (/1.029545,0.913893,0.833306,0.772561/)
-
-XNUMB_SPEED_LKT1D(44,1:4) = (/0.595899,0.545691,0.508686,0.479638/)
-XMASS_SPEED_LKT1D(44,1:4) = (/1.052735,0.933071,0.849927,0.787379/)
-
-XNUMB_SPEED_LKT1D(45,1:4) = (/0.612272,0.559780,0.521219,0.491021/)
-XMASS_SPEED_LKT1D(45,1:4) = (/1.075673,0.952016,0.866335,0.802000/)
-
-XNUMB_SPEED_LKT1D(46,1:4) = (/0.628584,0.573789,0.533664,0.502312/)
-XMASS_SPEED_LKT1D(46,1:4) = (/1.098364,0.970733,0.882533,0.816428/)
-
-XNUMB_SPEED_LKT1D(47,1:4) = (/0.644831,0.587714,0.546018,0.513511/)
-XMASS_SPEED_LKT1D(47,1:4) = (/1.120808,0.989225,0.898526,0.830668/)
-
-XNUMB_SPEED_LKT1D(48,1:4) = (/0.661009,0.601554,0.558281,0.524618/)
-XMASS_SPEED_LKT1D(48,1:4) = (/1.143009,1.007497,0.914320,0.844724/)
-
-XNUMB_SPEED_LKT1D(49,1:4) = (/0.677115,0.615307,0.570453,0.535632/)
-XMASS_SPEED_LKT1D(49,1:4) = (/1.164969,1.025553,0.929918,0.858602/)
-
-XNUMB_SPEED_LKT1D(50,1:4) = (/0.693145,0.628971,0.582533,0.546554/)
-XMASS_SPEED_LKT1D(50,1:4) = (/1.186693,1.043397,0.945325,0.872306/)
-
-XNUMB_SPEED_LKT1D(51,1:4) = (/0.709098,0.642547,0.594520,0.557383/)
-XMASS_SPEED_LKT1D(51,1:4) = (/1.208184,1.061034,0.960546,0.885840/)
-
-XNUMB_SPEED_LKT1D(52,1:4) = (/0.724970,0.656031,0.606415,0.568121/)
-XMASS_SPEED_LKT1D(52,1:4) = (/1.229445,1.078468,0.975584,0.899209/)
-
-XNUMB_SPEED_LKT1D(53,1:4) = (/0.740760,0.669425,0.618217,0.578768/)
-XMASS_SPEED_LKT1D(53,1:4) = (/1.250479,1.095703,0.990445,0.912417/)
-
-XNUMB_SPEED_LKT1D(54,1:4) = (/0.756466,0.682727,0.629927,0.589324/)
-XMASS_SPEED_LKT1D(54,1:4) = (/1.271291,1.112743,1.005133,0.925468/)
-
-XNUMB_SPEED_LKT1D(55,1:4) = (/0.772086,0.695937,0.641545,0.599790/)
-XMASS_SPEED_LKT1D(55,1:4) = (/1.291884,1.129592,1.019651,0.938366/)
-
-XNUMB_SPEED_LKT1D(56,1:4) = (/0.787619,0.709055,0.653071,0.610166/)
-XMASS_SPEED_LKT1D(56,1:4) = (/1.312262,1.146255,1.034004,0.951114/)
-
-XNUMB_SPEED_LKT1D(57,1:4) = (/0.803064,0.722080,0.664505,0.620455/)
-XMASS_SPEED_LKT1D(57,1:4) = (/1.332428,1.162735,1.048195,0.963717/)
-
-XNUMB_SPEED_LKT1D(58,1:4) = (/0.818419,0.735013,0.675849,0.630655/)
-XMASS_SPEED_LKT1D(58,1:4) = (/1.352387,1.179036,1.062228,0.976179/)
-
-XNUMB_SPEED_LKT1D(59,1:4) = (/0.833684,0.747853,0.687103,0.640769/)
-XMASS_SPEED_LKT1D(59,1:4) = (/1.372142,1.195163,1.076108,0.988502/)
-
-XNUMB_SPEED_LKT1D(60,1:4) = (/0.848858,0.760601,0.698268,0.650797/)
-XMASS_SPEED_LKT1D(60,1:4) = (/1.391697,1.211118,1.089838,1.000691/)
-
-XNUMB_SPEED_LKT1D(61,1:4) = (/0.863940,0.773257,0.709343,0.660740/)
-XMASS_SPEED_LKT1D(61,1:4) = (/1.411055,1.226906,1.103420,1.012748/)
-
-XNUMB_SPEED_LKT1D(62,1:4) = (/0.878931,0.785821,0.720331,0.670599/)
-XMASS_SPEED_LKT1D(62,1:4) = (/1.430221,1.242531,1.116860,1.024677/)
-
-XNUMB_SPEED_LKT1D(63,1:4) = (/0.893829,0.798295,0.731231,0.680376/)
-XMASS_SPEED_LKT1D(63,1:4) = (/1.449198,1.257995,1.130160,1.036481/)
-
-XNUMB_SPEED_LKT1D(64,1:4) = (/0.908635,0.810678,0.742046,0.690070/)
-XMASS_SPEED_LKT1D(64,1:4) = (/1.467989,1.273303,1.143323,1.048163/)
-
-XNUMB_SPEED_LKT1D(65,1:4) = (/0.923348,0.822970,0.752774,0.699684/)
-XMASS_SPEED_LKT1D(65,1:4) = (/1.486598,1.288458,1.156352,1.059726/)
-
-XNUMB_SPEED_LKT1D(66,1:4) = (/0.937969,0.835174,0.763418,0.709218/)
-XMASS_SPEED_LKT1D(66,1:4) = (/1.505028,1.303463,1.169251,1.071173/)
-
-XNUMB_SPEED_LKT1D(67,1:4) = (/0.952498,0.847288,0.773979,0.718673/)
-XMASS_SPEED_LKT1D(67,1:4) = (/1.523284,1.318321,1.182023,1.082507/)
-
-XNUMB_SPEED_LKT1D(68,1:4) = (/0.966934,0.859315,0.784456,0.728051/)
-XMASS_SPEED_LKT1D(68,1:4) = (/1.541367,1.333036,1.194671,1.093730/)
-
-XNUMB_SPEED_LKT1D(69,1:4) = (/0.981278,0.871254,0.794852,0.737352/)
-XMASS_SPEED_LKT1D(69,1:4) = (/1.559283,1.347610,1.207196,1.104844/)
-
-XNUMB_SPEED_LKT1D(70,1:4) = (/0.995531,0.883107,0.805167,0.746578/)
-XMASS_SPEED_LKT1D(70,1:4) = (/1.577033,1.362047,1.219603,1.115853/)
-
-XNUMB_SPEED_LKT1D(71,1:4) = (/1.009692,0.894874,0.815403,0.755729/)
-XMASS_SPEED_LKT1D(71,1:4) = (/1.594621,1.376349,1.231894,1.126759/)
-
-XNUMB_SPEED_LKT1D(72,1:4) = (/1.023762,0.906556,0.825559,0.764807/)
-XMASS_SPEED_LKT1D(72,1:4) = (/1.612050,1.390520,1.244071,1.137565/)
-
-XNUMB_SPEED_LKT1D(73,1:4) = (/1.037742,0.918153,0.835638,0.773813/)
-XMASS_SPEED_LKT1D(73,1:4) = (/1.629324,1.404562,1.256136,1.148271/)
-
-XNUMB_SPEED_LKT1D(74,1:4) = (/1.051632,0.929668,0.845640,0.782747/)
-XMASS_SPEED_LKT1D(74,1:4) = (/1.646445,1.418477,1.268093,1.158882/)
-
-XNUMB_SPEED_LKT1D(75,1:4) = (/1.065432,0.941100,0.855567,0.791611/)
-XMASS_SPEED_LKT1D(75,1:4) = (/1.663415,1.432269,1.279944,1.169398/)
-
-XNUMB_SPEED_LKT1D(76,1:4) = (/1.079144,0.952451,0.865418,0.800406/)
-XMASS_SPEED_LKT1D(76,1:4) = (/1.680239,1.445941,1.291690,1.179822/)
-
-XNUMB_SPEED_LKT1D(77,1:4) = (/1.092767,0.963721,0.875196,0.809133/)
-XMASS_SPEED_LKT1D(77,1:4) = (/1.696919,1.459493,1.303335,1.190156/)
-
-XNUMB_SPEED_LKT1D(78,1:4) = (/1.106302,0.974911,0.884901,0.817793/)
-XMASS_SPEED_LKT1D(78,1:4) = (/1.713457,1.472930,1.314880,1.200402/)
-
-XNUMB_SPEED_LKT1D(79,1:4) = (/1.119751,0.986023,0.894534,0.826387/)
-XMASS_SPEED_LKT1D(79,1:4) = (/1.729856,1.486253,1.326328,1.210562/)
-
-XNUMB_SPEED_LKT1D(80,1:4) = (/1.133114,0.997057,0.904096,0.834916/)
-XMASS_SPEED_LKT1D(80,1:4) = (/1.746119,1.499464,1.337680,1.220637/)
-
-XNUMB_SPEED_LKT1D(81,1:4) = (/1.146391,1.008013,0.913588,0.843380/)
-XMASS_SPEED_LKT1D(81,1:4) = (/1.762249,1.512567,1.348938,1.230630/)
-
-XNUMB_SPEED_LKT1D(82,1:4) = (/1.159583,1.018894,0.923011,0.851782/)
-XMASS_SPEED_LKT1D(82,1:4) = (/1.778248,1.525562,1.360105,1.240542/)
-
-XNUMB_SPEED_LKT1D(83,1:4) = (/1.172691,1.029699,0.932367,0.860121/)
-XMASS_SPEED_LKT1D(83,1:4) = (/1.794117,1.538452,1.371182,1.250374/)
-
-XNUMB_SPEED_LKT1D(84,1:4) = (/1.185716,1.040431,0.941655,0.868399/)
-XMASS_SPEED_LKT1D(84,1:4) = (/1.809861,1.551240,1.382171,1.260130/)
-
-XNUMB_SPEED_LKT1D(85,1:4) = (/1.198659,1.051089,0.950877,0.876616/)
-XMASS_SPEED_LKT1D(85,1:4) = (/1.825481,1.563927,1.393074,1.269809/)
-
-XNUMB_SPEED_LKT1D(86,1:4) = (/1.211519,1.061674,0.960034,0.884774/)
-XMASS_SPEED_LKT1D(86,1:4) = (/1.840979,1.576515,1.403892,1.279414/)
-
-XNUMB_SPEED_LKT1D(87,1:4) = (/1.224299,1.072188,0.969127,0.892874/)
-XMASS_SPEED_LKT1D(87,1:4) = (/1.856357,1.589006,1.414628,1.288946/)
-
-XNUMB_SPEED_LKT1D(88,1:4) = (/1.236999,1.082631,0.978157,0.900916/)
-XMASS_SPEED_LKT1D(88,1:4) = (/1.871619,1.601402,1.425282,1.298406/)
-
-XNUMB_SPEED_LKT1D(89,1:4) = (/1.249619,1.093005,0.987124,0.908901/)
-XMASS_SPEED_LKT1D(89,1:4) = (/1.886765,1.613704,1.435857,1.307797/)
-
-XNUMB_SPEED_LKT1D(90,1:4) = (/1.262160,1.103309,0.996029,0.916830/)
-XMASS_SPEED_LKT1D(90,1:4) = (/1.901798,1.625916,1.446354,1.317118/)
-
-XNUMB_SPEED_LKT1D(91,1:4) = (/1.274624,1.113546,1.004874,0.924705/)
-XMASS_SPEED_LKT1D(91,1:4) = (/1.916720,1.638037,1.456774,1.326373/)
-
-XNUMB_SPEED_LKT1D(92,1:4) = (/1.287011,1.123716,1.013659,0.932525/)
-XMASS_SPEED_LKT1D(92,1:4) = (/1.931533,1.650070,1.467120,1.335561/)
-
-XNUMB_SPEED_LKT1D(93,1:4) = (/1.299322,1.133819,1.022385,0.940291/)
-XMASS_SPEED_LKT1D(93,1:4) = (/1.946239,1.662017,1.477391,1.344684/)
-
-XNUMB_SPEED_LKT1D(94,1:4) = (/1.311558,1.143857,1.031053,0.948005/)
-XMASS_SPEED_LKT1D(94,1:4) = (/1.960839,1.673878,1.487590,1.353744/)
-
-XNUMB_SPEED_LKT1D(95,1:4) = (/1.323719,1.153830,1.039663,0.955667/)
-XMASS_SPEED_LKT1D(95,1:4) = (/1.975336,1.685656,1.497718,1.362741/)
-
-XNUMB_SPEED_LKT1D(96,1:4) = (/1.335806,1.163739,1.048217,0.963278/)
-XMASS_SPEED_LKT1D(96,1:4) = (/1.989732,1.697353,1.507777,1.371677/)
-
-XNUMB_SPEED_LKT1D(97,1:4) = (/1.347820,1.173586,1.056715,0.970839/)
-XMASS_SPEED_LKT1D(97,1:4) = (/2.004027,1.708969,1.517766,1.380553/)
-
-XNUMB_SPEED_LKT1D(98,1:4) = (/1.359762,1.183370,1.065158,0.978350/)
-XMASS_SPEED_LKT1D(98,1:4) = (/2.018225,1.720505,1.527689,1.389370/)
-
-XNUMB_SPEED_LKT1D(99,1:4) = (/1.371633,1.193093,1.073547,0.985812/)
-XMASS_SPEED_LKT1D(99,1:4) = (/2.032326,1.731964,1.537546,1.398129/)
-
-XNUMB_SPEED_LKT1D(100,1:4) = (/1.383434,1.202756,1.081883,0.993226/)
-XMASS_SPEED_LKT1D(100,1:4) = (/2.046332,1.743347,1.547338,1.406831/)
-
-XNUMB_SPEED_LKT1D(101,1:4) = (/1.395164,1.212358,1.090166,1.000593/)
-XMASS_SPEED_LKT1D(101,1:4) = (/2.060244,1.754655,1.557066,1.415476/)
-
-XNUMB_SPEED_LKT1D(102,1:4) = (/1.406826,1.221902,1.098397,1.007913/)
-XMASS_SPEED_LKT1D(102,1:4) = (/2.074065,1.765889,1.566731,1.424067/)
-
-XNUMB_SPEED_LKT1D(103,1:4) = (/1.418419,1.231388,1.106576,1.015187/)
-XMASS_SPEED_LKT1D(103,1:4) = (/2.087796,1.777051,1.576336,1.432604/)
-
-XNUMB_SPEED_LKT1D(104,1:4) = (/1.429945,1.240816,1.114705,1.022415/)
-XMASS_SPEED_LKT1D(104,1:4) = (/2.101438,1.788142,1.585879,1.441088/)
-
-XNUMB_SPEED_LKT1D(105,1:4) = (/1.441405,1.250187,1.122785,1.029599/)
-XMASS_SPEED_LKT1D(105,1:4) = (/2.114993,1.799162,1.595363,1.449519/)
-
-XNUMB_SPEED_LKT1D(106,1:4) = (/1.452798,1.259502,1.130815,1.036739/)
-XMASS_SPEED_LKT1D(106,1:4) = (/2.128462,1.810114,1.604789,1.457899/)
-
-XNUMB_SPEED_LKT1D(107,1:4) = (/1.464126,1.268762,1.138796,1.043835/)
-XMASS_SPEED_LKT1D(107,1:4) = (/2.141846,1.820998,1.614157,1.466229/)
-
-XNUMB_SPEED_LKT1D(108,1:4) = (/1.475390,1.277968,1.146730,1.050888/)
-XMASS_SPEED_LKT1D(108,1:4) = (/2.155147,1.831815,1.623469,1.474509/)
-
-XNUMB_SPEED_LKT1D(109,1:4) = (/1.486590,1.287119,1.154617,1.057899/)
-XMASS_SPEED_LKT1D(109,1:4) = (/2.168367,1.842567,1.632725,1.482740/)
-
-XNUMB_SPEED_LKT1D(110,1:4) = (/1.497727,1.296218,1.162457,1.064869/)
-XMASS_SPEED_LKT1D(110,1:4) = (/2.181506,1.853254,1.641926,1.490923/)
-
-XNUMB_SPEED_LKT1D(111,1:4) = (/1.508802,1.305264,1.170252,1.071797/)
-XMASS_SPEED_LKT1D(111,1:4) = (/2.194566,1.863878,1.651074,1.499059/)
-
-XNUMB_SPEED_LKT1D(112,1:4) = (/1.519815,1.314258,1.178001,1.078686/)
-XMASS_SPEED_LKT1D(112,1:4) = (/2.207547,1.874440,1.660168,1.507149/)
-
-XNUMB_SPEED_LKT1D(113,1:4) = (/1.530768,1.323200,1.185705,1.085534/)
-XMASS_SPEED_LKT1D(113,1:4) = (/2.220452,1.884940,1.669211,1.515192/)
-
-XNUMB_SPEED_LKT1D(114,1:4) = (/1.541660,1.332093,1.193365,1.092342/)
-XMASS_SPEED_LKT1D(114,1:4) = (/2.233281,1.895380,1.678202,1.523191/)
-
-XNUMB_SPEED_LKT1D(115,1:4) = (/1.552493,1.340935,1.200982,1.099112/)
-XMASS_SPEED_LKT1D(115,1:4) = (/2.246036,1.905760,1.687143,1.531145/)
-
-XNUMB_SPEED_LKT1D(116,1:4) = (/1.563267,1.349728,1.208556,1.105844/)
-XMASS_SPEED_LKT1D(116,1:4) = (/2.258718,1.916081,1.696034,1.539056/)
-
-XNUMB_SPEED_LKT1D(117,1:4) = (/1.573983,1.358473,1.216088,1.112538/)
-XMASS_SPEED_LKT1D(117,1:4) = (/2.271327,1.926345,1.704876,1.546924/)
-
-XNUMB_SPEED_LKT1D(118,1:4) = (/1.584642,1.367169,1.223578,1.119195/)
-XMASS_SPEED_LKT1D(118,1:4) = (/2.283865,1.936552,1.713669,1.554749/)
-
-XNUMB_SPEED_LKT1D(119,1:4) = (/1.595243,1.375818,1.231026,1.125815/)
-XMASS_SPEED_LKT1D(119,1:4) = (/2.296333,1.946702,1.722416,1.562533/)
-
-XNUMB_SPEED_LKT1D(120,1:4) = (/1.605789,1.384420,1.238434,1.132398/)
-XMASS_SPEED_LKT1D(120,1:4) = (/2.308732,1.956798,1.731115,1.570275/)
-
-XNUMB_SPEED_LKT1D(121,1:4) = (/1.616279,1.392976,1.245802,1.138946/)
-XMASS_SPEED_LKT1D(121,1:4) = (/2.321063,1.966839,1.739769,1.577977/)
-
-XNUMB_SPEED_LKT1D(122,1:4) = (/1.626714,1.401487,1.253130,1.145459/)
-XMASS_SPEED_LKT1D(122,1:4) = (/2.333327,1.976827,1.748377,1.585640/)
-
-XNUMB_SPEED_LKT1D(123,1:4) = (/1.637095,1.409952,1.260419,1.151937/)
-XMASS_SPEED_LKT1D(123,1:4) = (/2.345524,1.986762,1.756940,1.593263/)
-
-XNUMB_SPEED_LKT1D(124,1:4) = (/1.647423,1.418372,1.267670,1.158380/)
-XMASS_SPEED_LKT1D(124,1:4) = (/2.357657,1.996645,1.765460,1.600847/)
-
-XNUMB_SPEED_LKT1D(125,1:4) = (/1.657697,1.426749,1.274882,1.164790/)
-XMASS_SPEED_LKT1D(125,1:4) = (/2.369725,2.006477,1.773936,1.608394/)
-
-XNUMB_SPEED_LKT1D(126,1:4) = (/1.667919,1.435082,1.282057,1.171166/)
-XMASS_SPEED_LKT1D(126,1:4) = (/2.381730,2.016258,1.782369,1.615902/)
-
-XNUMB_SPEED_LKT1D(127,1:4) = (/1.678089,1.443372,1.289194,1.177509/)
-XMASS_SPEED_LKT1D(127,1:4) = (/2.393673,2.025990,1.790760,1.623374/)
-
-XNUMB_SPEED_LKT1D(128,1:4) = (/1.688208,1.451619,1.296295,1.183819/)
-XMASS_SPEED_LKT1D(128,1:4) = (/2.405554,2.035672,1.799109,1.630809/)
-
-XNUMB_SPEED_LKT1D(129,1:4) = (/1.698276,1.459825,1.303360,1.190098/)
-XMASS_SPEED_LKT1D(129,1:4) = (/2.417374,2.045306,1.807417,1.638208/)
-
-XNUMB_SPEED_LKT1D(130,1:4) = (/1.708294,1.467989,1.310389,1.196344/)
-XMASS_SPEED_LKT1D(130,1:4) = (/2.429135,2.054892,1.815685,1.645572/)
-
-XNUMB_SPEED_LKT1D(131,1:4) = (/1.718262,1.476113,1.317382,1.202560/)
-XMASS_SPEED_LKT1D(131,1:4) = (/2.440836,2.064431,1.823912,1.652900/)
-
-XNUMB_SPEED_LKT1D(132,1:4) = (/1.728182,1.484195,1.324341,1.208744/)
-XMASS_SPEED_LKT1D(132,1:4) = (/2.452479,2.073924,1.832101,1.660194/)
-
-XNUMB_SPEED_LKT1D(133,1:4) = (/1.738053,1.492238,1.331265,1.214897/)
-XMASS_SPEED_LKT1D(133,1:4) = (/2.464065,2.083371,1.840251,1.667454/)
-
-XNUMB_SPEED_LKT1D(134,1:4) = (/1.747876,1.500242,1.338155,1.221021/)
-XMASS_SPEED_LKT1D(134,1:4) = (/2.475594,2.092772,1.848362,1.674680/)
-
-XNUMB_SPEED_LKT1D(135,1:4) = (/1.757651,1.508206,1.345011,1.227114/)
-XMASS_SPEED_LKT1D(135,1:4) = (/2.487067,2.102129,1.856435,1.681873/)
-
-XNUMB_SPEED_LKT1D(136,1:4) = (/1.767380,1.516132,1.351834,1.233179/)
-XMASS_SPEED_LKT1D(136,1:4) = (/2.498485,2.111442,1.864471,1.689033/)
-
-XNUMB_SPEED_LKT1D(137,1:4) = (/1.777063,1.524020,1.358625,1.239214/)
-XMASS_SPEED_LKT1D(137,1:4) = (/2.509848,2.120712,1.872471,1.696161/)
-
-XNUMB_SPEED_LKT1D(138,1:4) = (/1.786699,1.531870,1.365383,1.245220/)
-XMASS_SPEED_LKT1D(138,1:4) = (/2.521158,2.129938,1.880434,1.703257/)
-
-XNUMB_SPEED_LKT1D(139,1:4) = (/1.796290,1.539683,1.372109,1.251198/)
-XMASS_SPEED_LKT1D(139,1:4) = (/2.532415,2.139123,1.888361,1.710321/)
-
-XNUMB_SPEED_LKT1D(140,1:4) = (/1.805837,1.547459,1.378803,1.257148/)
-XMASS_SPEED_LKT1D(140,1:4) = (/2.543619,2.148265,1.896253,1.717354/)
-
-XNUMB_SPEED_LKT1D(141,1:4) = (/1.815339,1.555198,1.385466,1.263070/)
-XMASS_SPEED_LKT1D(141,1:4) = (/2.554771,2.157367,1.904109,1.724357/)
-
-XNUMB_SPEED_LKT1D(142,1:4) = (/1.824797,1.562902,1.392098,1.268965/)
-XMASS_SPEED_LKT1D(142,1:4) = (/2.565873,2.166428,1.911932,1.731329/)
-
-XNUMB_SPEED_LKT1D(143,1:4) = (/1.834212,1.570570,1.398700,1.274833/)
-XMASS_SPEED_LKT1D(143,1:4) = (/2.576924,2.175448,1.919720,1.738272/)
-
-XNUMB_SPEED_LKT1D(144,1:4) = (/1.843584,1.578203,1.405271,1.280674/)
-XMASS_SPEED_LKT1D(144,1:4) = (/2.587925,2.184429,1.927474,1.745184/)
-
-XNUMB_SPEED_LKT1D(145,1:4) = (/1.852913,1.585801,1.411812,1.286489/)
-XMASS_SPEED_LKT1D(145,1:4) = (/2.598877,2.193370,1.935196,1.752068/)
-
-XNUMB_SPEED_LKT1D(146,1:4) = (/1.862200,1.593365,1.418324,1.292277/)
-XMASS_SPEED_LKT1D(146,1:4) = (/2.609780,2.202273,1.942884,1.758922/)
-
-XNUMB_SPEED_LKT1D(147,1:4) = (/1.871446,1.600895,1.424807,1.298040/)
-XMASS_SPEED_LKT1D(147,1:4) = (/2.620636,2.211138,1.950540,1.765749/)
-
-XNUMB_SPEED_LKT1D(148,1:4) = (/1.880651,1.608391,1.431261,1.303778/)
-XMASS_SPEED_LKT1D(148,1:4) = (/2.631444,2.219964,1.958164,1.772547/)
-
-XNUMB_SPEED_LKT1D(149,1:4) = (/1.889815,1.615855,1.437687,1.309490/)
-XMASS_SPEED_LKT1D(149,1:4) = (/2.642205,2.228754,1.965756,1.779317/)
-
-XNUMB_SPEED_LKT1D(150,1:4) = (/1.898938,1.623285,1.444084,1.315177/)
-XMASS_SPEED_LKT1D(150,1:4) = (/2.652920,2.237507,1.973317,1.786059/)
-
-XNUMB_SPEED_LKT1D(151,1:4) = (/1.908022,1.630682,1.450454,1.320840/)
-XMASS_SPEED_LKT1D(151,1:4) = (/2.663589,2.246223,1.980847,1.792775/)
-
-XNUMB_SPEED_LKT1D(152,1:4) = (/1.917067,1.638048,1.456795,1.326478/)
-XMASS_SPEED_LKT1D(152,1:4) = (/2.674213,2.254903,1.988347,1.799464/)
-
-XNUMB_SPEED_LKT1D(153,1:4) = (/1.926072,1.645382,1.463110,1.332092/)
-XMASS_SPEED_LKT1D(153,1:4) = (/2.684792,2.263547,1.995816,1.806126/)
-
-XNUMB_SPEED_LKT1D(154,1:4) = (/1.935039,1.652684,1.469398,1.337683/)
-XMASS_SPEED_LKT1D(154,1:4) = (/2.695328,2.272156,2.003255,1.812761/)
-
-XNUMB_SPEED_LKT1D(155,1:4) = (/1.943967,1.659955,1.475659,1.343249/)
-XMASS_SPEED_LKT1D(155,1:4) = (/2.705819,2.280731,2.010665,1.819371/)
-
-XNUMB_SPEED_LKT1D(156,1:4) = (/1.952858,1.667196,1.481894,1.348793/)
-XMASS_SPEED_LKT1D(156,1:4) = (/2.716268,2.289271,2.018046,1.825955/)
-
-XNUMB_SPEED_LKT1D(157,1:4) = (/1.961711,1.674406,1.488102,1.354314/)
-XMASS_SPEED_LKT1D(157,1:4) = (/2.726673,2.297777,2.025398,1.832514/)
-
-XNUMB_SPEED_LKT1D(158,1:4) = (/1.970527,1.681586,1.494285,1.359812/)
-XMASS_SPEED_LKT1D(158,1:4) = (/2.737037,2.306249,2.032721,1.839048/)
-
-XNUMB_SPEED_LKT1D(159,1:4) = (/1.979307,1.688736,1.500443,1.365287/)
-XMASS_SPEED_LKT1D(159,1:4) = (/2.747359,2.314689,2.040016,1.845557/)
-
-XNUMB_SPEED_LKT1D(160,1:4) = (/1.988050,1.695857,1.506575,1.370740/)
-XMASS_SPEED_LKT1D(160,1:4) = (/2.757640,2.323095,2.047283,1.852041/)
-
-XNUMB_SPEED_LKT1D(161,1:4) = (/1.996758,1.702948,1.512682,1.376171/)
-XMASS_SPEED_LKT1D(161,1:4) = (/2.767879,2.331469,2.054523,1.858501/)
-
-XNUMB_SPEED_LKT1D(162,1:4) = (/2.005429,1.710011,1.518765,1.381580/)
-XMASS_SPEED_LKT1D(162,1:4) = (/2.778079,2.339810,2.061735,1.864937/)
-
-XNUMB_SPEED_LKT1D(163,1:4) = (/2.014066,1.717045,1.524823,1.386968/)
-XMASS_SPEED_LKT1D(163,1:4) = (/2.788238,2.348120,2.068920,1.871349/)
-
-XNUMB_SPEED_LKT1D(164,1:4) = (/2.022667,1.724051,1.530857,1.392335/)
-XMASS_SPEED_LKT1D(164,1:4) = (/2.798359,2.356398,2.076078,1.877737/)
-
-XNUMB_SPEED_LKT1D(165,1:4) = (/2.031235,1.731029,1.536867,1.397680/)
-XMASS_SPEED_LKT1D(165,1:4) = (/2.808440,2.364645,2.083210,1.884102/)
-
-XNUMB_SPEED_LKT1D(166,1:4) = (/2.039768,1.737979,1.542854,1.403004/)
-XMASS_SPEED_LKT1D(166,1:4) = (/2.818482,2.372861,2.090316,1.890445/)
-
-XNUMB_SPEED_LKT1D(167,1:4) = (/2.048267,1.744901,1.548817,1.408308/)
-XMASS_SPEED_LKT1D(167,1:4) = (/2.828486,2.381046,2.097396,1.896764/)
-
-XNUMB_SPEED_LKT1D(168,1:4) = (/2.056733,1.751797,1.554757,1.413592/)
-XMASS_SPEED_LKT1D(168,1:4) = (/2.838452,2.389202,2.104450,1.903061/)
-
-XNUMB_SPEED_LKT1D(169,1:4) = (/2.065165,1.758666,1.560673,1.418855/)
-XMASS_SPEED_LKT1D(169,1:4) = (/2.848381,2.397327,2.111478,1.909335/)
-
-XNUMB_SPEED_LKT1D(170,1:4) = (/2.073565,1.765508,1.566568,1.424098/)
-XMASS_SPEED_LKT1D(170,1:4) = (/2.858272,2.405422,2.118481,1.915587/)
-
-XNUMB_SPEED_LKT1D(171,1:4) = (/2.081932,1.772324,1.572439,1.429321/)
-XMASS_SPEED_LKT1D(171,1:4) = (/2.868127,2.413489,2.125460,1.921817/)
-
-XNUMB_SPEED_LKT1D(172,1:4) = (/2.090267,1.779113,1.578289,1.434525/)
-XMASS_SPEED_LKT1D(172,1:4) = (/2.877945,2.421526,2.132413,1.928025/)
-
-XNUMB_SPEED_LKT1D(173,1:4) = (/2.098570,1.785877,1.584116,1.439709/)
-XMASS_SPEED_LKT1D(173,1:4) = (/2.887727,2.429534,2.139343,1.934212/)
-
-XNUMB_SPEED_LKT1D(174,1:4) = (/2.106841,1.792616,1.589922,1.444874/)
-XMASS_SPEED_LKT1D(174,1:4) = (/2.897473,2.437514,2.146247,1.940378/)
-
-XNUMB_SPEED_LKT1D(175,1:4) = (/2.115081,1.799329,1.595706,1.450021/)
-XMASS_SPEED_LKT1D(175,1:4) = (/2.907184,2.445466,2.153128,1.946522/)
-
-XNUMB_SPEED_LKT1D(176,1:4) = (/2.123290,1.806017,1.601469,1.455148/)
-XMASS_SPEED_LKT1D(176,1:4) = (/2.916860,2.453389,2.159985,1.952645/)
-
-XNUMB_SPEED_LKT1D(177,1:4) = (/2.131469,1.812680,1.607211,1.460256/)
-XMASS_SPEED_LKT1D(177,1:4) = (/2.926502,2.461285,2.166819,1.958747/)
-
-XNUMB_SPEED_LKT1D(178,1:4) = (/2.139617,1.819319,1.612931,1.465347/)
-XMASS_SPEED_LKT1D(178,1:4) = (/2.936108,2.469153,2.173629,1.964829/)
-
-XNUMB_SPEED_LKT1D(179,1:4) = (/2.147735,1.825934,1.618631,1.470418/)
-XMASS_SPEED_LKT1D(179,1:4) = (/2.945681,2.476994,2.180415,1.970891/)
-
-XNUMB_SPEED_LKT1D(180,1:4) = (/2.155823,1.832524,1.624310,1.475472/)
-XMASS_SPEED_LKT1D(180,1:4) = (/2.955219,2.484808,2.187179,1.976932/)
-
-XNUMB_SPEED_LKT1D(181,1:4) = (/2.163881,1.839091,1.629969,1.480508/)
-XMASS_SPEED_LKT1D(181,1:4) = (/2.964725,2.492595,2.193920,1.982952/)
-
-XNUMB_SPEED_LKT1D(182,1:4) = (/2.171911,1.845633,1.635608,1.485526/)
-XMASS_SPEED_LKT1D(182,1:4) = (/2.974196,2.500355,2.200638,1.988953/)
-
-XNUMB_SPEED_LKT1D(183,1:4) = (/2.179911,1.852153,1.641227,1.490527/)
-XMASS_SPEED_LKT1D(183,1:4) = (/2.983635,2.508089,2.207334,1.994935/)
-
-XNUMB_SPEED_LKT1D(184,1:4) = (/2.187882,1.858649,1.646826,1.495510/)
-XMASS_SPEED_LKT1D(184,1:4) = (/2.993041,2.515797,2.214008,2.000896/)
-
-XNUMB_SPEED_LKT1D(185,1:4) = (/2.195825,1.865123,1.652406,1.500476/)
-XMASS_SPEED_LKT1D(185,1:4) = (/3.002415,2.523478,2.220659,2.006838/)
-
-XNUMB_SPEED_LKT1D(186,1:4) = (/2.203740,1.871573,1.657966,1.505425/)
-XMASS_SPEED_LKT1D(186,1:4) = (/3.011756,2.531134,2.227288,2.012761/)
-
-XNUMB_SPEED_LKT1D(187,1:4) = (/2.211626,1.878001,1.663507,1.510356/)
-XMASS_SPEED_LKT1D(187,1:4) = (/3.021065,2.538765,2.233896,2.018664/)
-
-XNUMB_SPEED_LKT1D(188,1:4) = (/2.219485,1.884407,1.669029,1.515272/)
-XMASS_SPEED_LKT1D(188,1:4) = (/3.030343,2.546369,2.240482,2.024548/)
-
-XNUMB_SPEED_LKT1D(189,1:4) = (/2.227317,1.890790,1.674532,1.520170/)
-XMASS_SPEED_LKT1D(189,1:4) = (/3.039589,2.553949,2.247046,2.030413/)
-
-XNUMB_SPEED_LKT1D(190,1:4) = (/2.235121,1.897152,1.680017,1.525052/)
-XMASS_SPEED_LKT1D(190,1:4) = (/3.048803,2.561504,2.253589,2.036260/)
-
-XNUMB_SPEED_LKT1D(191,1:4) = (/2.242899,1.903492,1.685483,1.529918/)
-XMASS_SPEED_LKT1D(191,1:4) = (/3.057987,2.569033,2.260111,2.042088/)
-
-XNUMB_SPEED_LKT1D(192,1:4) = (/2.250649,1.909810,1.690930,1.534767/)
-XMASS_SPEED_LKT1D(192,1:4) = (/3.067139,2.576538,2.266612,2.047897/)
-
-XNUMB_SPEED_LKT1D(193,1:4) = (/2.258373,1.916107,1.696360,1.539601/)
-XMASS_SPEED_LKT1D(193,1:4) = (/3.076261,2.584018,2.273091,2.053687/)
-
-XNUMB_SPEED_LKT1D(194,1:4) = (/2.266071,1.922383,1.701771,1.544419/)
-XMASS_SPEED_LKT1D(194,1:4) = (/3.085353,2.591474,2.279550,2.059459/)
-
-XNUMB_SPEED_LKT1D(195,1:4) = (/2.273743,1.928638,1.707164,1.549221/)
-XMASS_SPEED_LKT1D(195,1:4) = (/3.094414,2.598906,2.285989,2.065213/)
-
-XNUMB_SPEED_LKT1D(196,1:4) = (/2.281388,1.934872,1.712540,1.554007/)
-XMASS_SPEED_LKT1D(196,1:4) = (/3.103445,2.606313,2.292406,2.070949/)
-
-ELSE
-
-XNUMB_SPEED_LKT1D(1,1:4) = (/0.012923,0.010176,0.008721,0.007819/)
-XMASS_SPEED_LKT1D(1,1:4) = (/0.018036,0.016773,0.016098,0.015674/)
-
-XNUMB_SPEED_LKT1D(2,1:4) = (/0.013461,0.011087,0.009829,0.009049/)
-XMASS_SPEED_LKT1D(2,1:4) = (/0.023476,0.022401,0.021816,0.021441/)
-
-XNUMB_SPEED_LKT1D(3,1:4) = (/0.014584,0.012500,0.011393,0.010706/)
-XMASS_SPEED_LKT1D(3,1:4) = (/0.030187,0.029222,0.028679,0.028317/)
-
-XNUMB_SPEED_LKT1D(4,1:4) = (/0.016219,0.014360,0.013371,0.012754/)
-XMASS_SPEED_LKT1D(4,1:4) = (/0.038072,0.037155,0.036609,0.036224/)
-
-XNUMB_SPEED_LKT1D(5,1:4) = (/0.018310,0.016629,0.015730,0.015165/)
-XMASS_SPEED_LKT1D(5,1:4) = (/0.047064,0.046132,0.045537,0.045090/)
-
-XNUMB_SPEED_LKT1D(6,1:4) = (/0.020817,0.019276,0.018445,0.017918/)
-XMASS_SPEED_LKT1D(6,1:4) = (/0.057103,0.056092,0.055395,0.054841/)
-
-XNUMB_SPEED_LKT1D(7,1:4) = (/0.023711,0.022279,0.021496,0.020993/)
-XMASS_SPEED_LKT1D(7,1:4) = (/0.068133,0.066973,0.066115,0.065403/)
-
-XNUMB_SPEED_LKT1D(8,1:4) = (/0.026970,0.025618,0.024866,0.024373/)
-XMASS_SPEED_LKT1D(8,1:4) = (/0.080101,0.078711,0.077628,0.076702/)
-
-XNUMB_SPEED_LKT1D(9,1:4) = (/0.030573,0.029277,0.028539,0.028043/)
-XMASS_SPEED_LKT1D(9,1:4) = (/0.092951,0.091243,0.089862,0.088661/)
-
-XNUMB_SPEED_LKT1D(10,1:4) = (/0.034506,0.033242,0.032501,0.031988/)
-XMASS_SPEED_LKT1D(10,1:4) = (/0.106627,0.104505,0.102750,0.101209/)
-
-XNUMB_SPEED_LKT1D(11,1:4) = (/0.038754,0.037498,0.036737,0.036194/)
-XMASS_SPEED_LKT1D(11,1:4) = (/0.121072,0.118434,0.116222,0.114275/)
-
-XNUMB_SPEED_LKT1D(12,1:4) = (/0.043304,0.042033,0.041234,0.040646/)
-XMASS_SPEED_LKT1D(12,1:4) = (/0.136230,0.132967,0.130214,0.127791/)
-
-XNUMB_SPEED_LKT1D(13,1:4) = (/0.048144,0.046834,0.045980,0.045332/)
-XMASS_SPEED_LKT1D(13,1:4) = (/0.152045,0.148043,0.144661,0.141694/)
-
-XNUMB_SPEED_LKT1D(14,1:4) = (/0.053263,0.051889,0.050960,0.050238/)
-XMASS_SPEED_LKT1D(14,1:4) = (/0.168461,0.163602,0.159505,0.155926/)
-
-XNUMB_SPEED_LKT1D(15,1:4) = (/0.058650,0.057188,0.056164,0.055350/)
-XMASS_SPEED_LKT1D(15,1:4) = (/0.185425,0.179591,0.174690,0.170431/)
-
-XNUMB_SPEED_LKT1D(16,1:4) = (/0.064293,0.062716,0.061579,0.060657/)
-XMASS_SPEED_LKT1D(16,1:4) = (/0.202884,0.195954,0.190165,0.185160/)
-
-XNUMB_SPEED_LKT1D(17,1:4) = (/0.070183,0.068465,0.067192,0.066145/)
-XMASS_SPEED_LKT1D(17,1:4) = (/0.220788,0.212643,0.205881,0.200068/)
-
-XNUMB_SPEED_LKT1D(18,1:4) = (/0.076309,0.074421,0.072992,0.071802/)
-XMASS_SPEED_LKT1D(18,1:4) = (/0.239089,0.229611,0.221794,0.215114/)
-
-XNUMB_SPEED_LKT1D(19,1:4) = (/0.082660,0.080575,0.078967,0.077617/)
-XMASS_SPEED_LKT1D(19,1:4) = (/0.257743,0.246815,0.237865,0.230261/)
-
-XNUMB_SPEED_LKT1D(20,1:4) = (/0.089227,0.086915,0.085107,0.083579/)
-XMASS_SPEED_LKT1D(20,1:4) = (/0.276705,0.264215,0.254057,0.245475/)
-
-XNUMB_SPEED_LKT1D(21,1:4) = (/0.095999,0.093431,0.091400,0.089677/)
-XMASS_SPEED_LKT1D(21,1:4) = (/0.295936,0.281775,0.270337,0.260729/)
-
-XNUMB_SPEED_LKT1D(22,1:4) = (/0.102968,0.100113,0.097837,0.095901/)
-XMASS_SPEED_LKT1D(22,1:4) = (/0.315398,0.299461,0.286676,0.275995/)
-
-XNUMB_SPEED_LKT1D(23,1:4) = (/0.110123,0.106952,0.104407,0.102240/)
-XMASS_SPEED_LKT1D(23,1:4) = (/0.335056,0.317241,0.303047,0.291251/)
-
-XNUMB_SPEED_LKT1D(24,1:4) = (/0.117455,0.113937,0.111101,0.108685/)
-XMASS_SPEED_LKT1D(24,1:4) = (/0.354877,0.335090,0.319428,0.306477/)
-
-XNUMB_SPEED_LKT1D(25,1:4) = (/0.124956,0.121059,0.117911,0.115228/)
-XMASS_SPEED_LKT1D(25,1:4) = (/0.374832,0.352981,0.335796,0.321654/)
-
-XNUMB_SPEED_LKT1D(26,1:4) = (/0.132616,0.128310,0.124826,0.121860/)
-XMASS_SPEED_LKT1D(26,1:4) = (/0.394892,0.370892,0.352133,0.336769/)
-
-XNUMB_SPEED_LKT1D(27,1:4) = (/0.140427,0.135680,0.131839,0.128573/)
-XMASS_SPEED_LKT1D(27,1:4) = (/0.415033,0.388802,0.368423,0.351807/)
-
-XNUMB_SPEED_LKT1D(28,1:4) = (/0.148381,0.143163,0.138942,0.135360/)
-XMASS_SPEED_LKT1D(28,1:4) = (/0.435230,0.406693,0.384651,0.366757/)
-
-XNUMB_SPEED_LKT1D(29,1:4) = (/0.156469,0.150749,0.146127,0.142212/)
-XMASS_SPEED_LKT1D(29,1:4) = (/0.455462,0.424548,0.400805,0.381609/)
-
-XNUMB_SPEED_LKT1D(30,1:4) = (/0.164684,0.158432,0.153387,0.149124/)
-XMASS_SPEED_LKT1D(30,1:4) = (/0.475710,0.442354,0.416873,0.396356/)
-
-XNUMB_SPEED_LKT1D(31,1:4) = (/0.173019,0.166204,0.160716,0.156089/)
-XMASS_SPEED_LKT1D(31,1:4) = (/0.495955,0.460096,0.432847,0.410989/)
-
-XNUMB_SPEED_LKT1D(32,1:4) = (/0.181465,0.174058,0.168107,0.163101/)
-XMASS_SPEED_LKT1D(32,1:4) = (/0.516182,0.477765,0.448718,0.425503/)
-
-XNUMB_SPEED_LKT1D(33,1:4) = (/0.190017,0.181988,0.175554,0.170155/)
-XMASS_SPEED_LKT1D(33,1:4) = (/0.536376,0.495348,0.464478,0.439893/)
-
-XNUMB_SPEED_LKT1D(34,1:4) = (/0.198667,0.189987,0.183050,0.177245/)
-XMASS_SPEED_LKT1D(34,1:4) = (/0.556524,0.512839,0.480122,0.454155/)
-
-XNUMB_SPEED_LKT1D(35,1:4) = (/0.207409,0.198050,0.190592,0.184366/)
-XMASS_SPEED_LKT1D(35,1:4) = (/0.576614,0.530228,0.495645,0.468287/)
-
-XNUMB_SPEED_LKT1D(36,1:4) = (/0.216236,0.206171,0.198173,0.191513/)
-XMASS_SPEED_LKT1D(36,1:4) = (/0.596634,0.547509,0.511043,0.482285/)
-
-XNUMB_SPEED_LKT1D(37,1:4) = (/0.225143,0.214345,0.205789,0.198683/)
-XMASS_SPEED_LKT1D(37,1:4) = (/0.616576,0.564676,0.526312,0.496148/)
-
-XNUMB_SPEED_LKT1D(38,1:4) = (/0.234125,0.222565,0.213435,0.205871/)
-XMASS_SPEED_LKT1D(38,1:4) = (/0.636432,0.581725,0.541450,0.509875/)
-
-XNUMB_SPEED_LKT1D(39,1:4) = (/0.243175,0.230829,0.221107,0.213073/)
-XMASS_SPEED_LKT1D(39,1:4) = (/0.656192,0.598651,0.556454,0.523465/)
-
-XNUMB_SPEED_LKT1D(40,1:4) = (/0.252288,0.239130,0.228801,0.220286/)
-XMASS_SPEED_LKT1D(40,1:4) = (/0.675851,0.615450,0.571323,0.536917/)
-
-XNUMB_SPEED_LKT1D(41,1:4) = (/0.261459,0.247465,0.236513,0.227507/)
-XMASS_SPEED_LKT1D(41,1:4) = (/0.695402,0.632120,0.586055,0.550232/)
-
-XNUMB_SPEED_LKT1D(42,1:4) = (/0.270684,0.255830,0.244240,0.234732/)
-XMASS_SPEED_LKT1D(42,1:4) = (/0.714841,0.648658,0.600651,0.563410/)
-
-XNUMB_SPEED_LKT1D(43,1:4) = (/0.279958,0.264220,0.251977,0.241958/)
-XMASS_SPEED_LKT1D(43,1:4) = (/0.734162,0.665062,0.615109,0.576451/)
-
-XNUMB_SPEED_LKT1D(44,1:4) = (/0.289277,0.272632,0.259723,0.249184/)
-XMASS_SPEED_LKT1D(44,1:4) = (/0.753363,0.681331,0.629429,0.589356/)
-
-XNUMB_SPEED_LKT1D(45,1:4) = (/0.298636,0.281063,0.267475,0.256405/)
-XMASS_SPEED_LKT1D(45,1:4) = (/0.772438,0.697464,0.643612,0.602127/)
-
-XNUMB_SPEED_LKT1D(46,1:4) = (/0.308031,0.289509,0.275228,0.263621/)
-XMASS_SPEED_LKT1D(46,1:4) = (/0.791386,0.713460,0.657659,0.614765/)
-
-XNUMB_SPEED_LKT1D(47,1:4) = (/0.317459,0.297967,0.282982,0.270828/)
-XMASS_SPEED_LKT1D(47,1:4) = (/0.810204,0.729318,0.671569,0.627270/)
-
-XNUMB_SPEED_LKT1D(48,1:4) = (/0.326916,0.306435,0.290733,0.278026/)
-XMASS_SPEED_LKT1D(48,1:4) = (/0.828890,0.745039,0.685344,0.639645/)
-
-XNUMB_SPEED_LKT1D(49,1:4) = (/0.336399,0.314908,0.298479,0.285211/)
-XMASS_SPEED_LKT1D(49,1:4) = (/0.847442,0.760622,0.698986,0.651891/)
-
-XNUMB_SPEED_LKT1D(50,1:4) = (/0.345905,0.323386,0.306218,0.292382/)
-XMASS_SPEED_LKT1D(50,1:4) = (/0.865859,0.776069,0.712494,0.664010/)
-
-XNUMB_SPEED_LKT1D(51,1:4) = (/0.355429,0.331865,0.313949,0.299539/)
-XMASS_SPEED_LKT1D(51,1:4) = (/0.884139,0.791379,0.725871,0.676004/)
-
-XNUMB_SPEED_LKT1D(52,1:4) = (/0.364970,0.340344,0.321668,0.306678/)
-XMASS_SPEED_LKT1D(52,1:4) = (/0.902284,0.806553,0.739118,0.687873/)
-
-XNUMB_SPEED_LKT1D(53,1:4) = (/0.374525,0.348819,0.329376,0.313799/)
-XMASS_SPEED_LKT1D(53,1:4) = (/0.920290,0.821592,0.752236,0.699621/)
-
-XNUMB_SPEED_LKT1D(54,1:4) = (/0.384090,0.357289,0.337069,0.320900/)
-XMASS_SPEED_LKT1D(54,1:4) = (/0.938160,0.836498,0.765227,0.711250/)
-
-XNUMB_SPEED_LKT1D(55,1:4) = (/0.393664,0.365752,0.344747,0.327981/)
-XMASS_SPEED_LKT1D(55,1:4) = (/0.955892,0.851270,0.778093,0.722760/)
-
-XNUMB_SPEED_LKT1D(56,1:4) = (/0.403243,0.374207,0.352408,0.335041/)
-XMASS_SPEED_LKT1D(56,1:4) = (/0.973486,0.865911,0.790835,0.734154/)
-
-XNUMB_SPEED_LKT1D(57,1:4) = (/0.412827,0.382651,0.360051,0.342077/)
-XMASS_SPEED_LKT1D(57,1:4) = (/0.990944,0.880422,0.803455,0.745434/)
-
-XNUMB_SPEED_LKT1D(58,1:4) = (/0.422411,0.391083,0.367675,0.349091/)
-XMASS_SPEED_LKT1D(58,1:4) = (/1.008265,0.894804,0.815955,0.756602/)
-
-XNUMB_SPEED_LKT1D(59,1:4) = (/0.431996,0.399502,0.375279,0.356079/)
-XMASS_SPEED_LKT1D(59,1:4) = (/1.025450,0.909058,0.828337,0.767660/)
-
-XNUMB_SPEED_LKT1D(60,1:4) = (/0.441578,0.407905,0.382861,0.363043/)
-XMASS_SPEED_LKT1D(60,1:4) = (/1.042500,0.923186,0.840602,0.778610/)
-
-XNUMB_SPEED_LKT1D(61,1:4) = (/0.451155,0.416293,0.390420,0.369981/)
-XMASS_SPEED_LKT1D(61,1:4) = (/1.059416,0.937190,0.852753,0.789453/)
-
-XNUMB_SPEED_LKT1D(62,1:4) = (/0.460727,0.424663,0.397957,0.376892/)
-XMASS_SPEED_LKT1D(62,1:4) = (/1.076198,0.951071,0.864790,0.800192/)
-
-XNUMB_SPEED_LKT1D(63,1:4) = (/0.470291,0.433014,0.405469,0.383776/)
-XMASS_SPEED_LKT1D(63,1:4) = (/1.092848,0.964831,0.876717,0.810828/)
-
-XNUMB_SPEED_LKT1D(64,1:4) = (/0.479845,0.441346,0.412957,0.390633/)
-XMASS_SPEED_LKT1D(64,1:4) = (/1.109366,0.978471,0.888534,0.821364/)
-
-XNUMB_SPEED_LKT1D(65,1:4) = (/0.489389,0.449658,0.420420,0.397462/)
-XMASS_SPEED_LKT1D(65,1:4) = (/1.125755,0.991992,0.900243,0.831801/)
-
-XNUMB_SPEED_LKT1D(66,1:4) = (/0.498921,0.457948,0.427856,0.404263/)
-XMASS_SPEED_LKT1D(66,1:4) = (/1.142015,1.005398,0.911847,0.842140/)
-
-XNUMB_SPEED_LKT1D(67,1:4) = (/0.508440,0.466215,0.435266,0.411034/)
-XMASS_SPEED_LKT1D(67,1:4) = (/1.158147,1.018689,0.923347,0.852385/)
-
-XNUMB_SPEED_LKT1D(68,1:4) = (/0.517944,0.474460,0.442648,0.417777/)
-XMASS_SPEED_LKT1D(68,1:4) = (/1.174153,1.031866,0.934744,0.862536/)
-
-XNUMB_SPEED_LKT1D(69,1:4) = (/0.527432,0.482681,0.450004,0.424490/)
-XMASS_SPEED_LKT1D(69,1:4) = (/1.190034,1.044933,0.946041,0.872595/)
-
-XNUMB_SPEED_LKT1D(70,1:4) = (/0.536903,0.490877,0.457331,0.431174/)
-XMASS_SPEED_LKT1D(70,1:4) = (/1.205791,1.057889,0.957240,0.882565/)
-
-XNUMB_SPEED_LKT1D(71,1:4) = (/0.546356,0.499047,0.464630,0.437828/)
-XMASS_SPEED_LKT1D(71,1:4) = (/1.221426,1.070738,0.968341,0.892445/)
-
-XNUMB_SPEED_LKT1D(72,1:4) = (/0.555790,0.507193,0.471900,0.444452/)
-XMASS_SPEED_LKT1D(72,1:4) = (/1.236941,1.083480,0.979347,0.902240/)
-
-XNUMB_SPEED_LKT1D(73,1:4) = (/0.565205,0.515311,0.479141,0.451046/)
-XMASS_SPEED_LKT1D(73,1:4) = (/1.252336,1.096117,0.990259,0.911949/)
-
-XNUMB_SPEED_LKT1D(74,1:4) = (/0.574599,0.523403,0.486353,0.457610/)
-XMASS_SPEED_LKT1D(74,1:4) = (/1.267613,1.108651,1.001079,0.921574/)
-
-XNUMB_SPEED_LKT1D(75,1:4) = (/0.583971,0.531468,0.493535,0.464143/)
-XMASS_SPEED_LKT1D(75,1:4) = (/1.282774,1.121083,1.011808,0.931117/)
-
-XNUMB_SPEED_LKT1D(76,1:4) = (/0.593321,0.539505,0.500688,0.470646/)
-XMASS_SPEED_LKT1D(76,1:4) = (/1.297820,1.133415,1.022448,0.940580/)
-
-XNUMB_SPEED_LKT1D(77,1:4) = (/0.602648,0.547514,0.507811,0.477119/)
-XMASS_SPEED_LKT1D(77,1:4) = (/1.312753,1.145648,1.033001,0.949964/)
-
-XNUMB_SPEED_LKT1D(78,1:4) = (/0.611952,0.555495,0.514903,0.483560/)
-XMASS_SPEED_LKT1D(78,1:4) = (/1.327573,1.157785,1.043468,0.959271/)
-
-XNUMB_SPEED_LKT1D(79,1:4) = (/0.621231,0.563447,0.521966,0.489972/)
-XMASS_SPEED_LKT1D(79,1:4) = (/1.342283,1.169826,1.053851,0.968501/)
-
-XNUMB_SPEED_LKT1D(80,1:4) = (/0.630486,0.571370,0.528998,0.496353/)
-XMASS_SPEED_LKT1D(80,1:4) = (/1.356884,1.181774,1.064150,0.977656/)
-
-XNUMB_SPEED_LKT1D(81,1:4) = (/0.639715,0.579264,0.536000,0.502703/)
-XMASS_SPEED_LKT1D(81,1:4) = (/1.371378,1.193628,1.074368,0.986738/)
-
-XNUMB_SPEED_LKT1D(82,1:4) = (/0.648918,0.587128,0.542971,0.509023/)
-XMASS_SPEED_LKT1D(82,1:4) = (/1.385765,1.205392,1.084505,0.995748/)
-
-XNUMB_SPEED_LKT1D(83,1:4) = (/0.658095,0.594963,0.549912,0.515313/)
-XMASS_SPEED_LKT1D(83,1:4) = (/1.400047,1.217066,1.094564,1.004687/)
-
-XNUMB_SPEED_LKT1D(84,1:4) = (/0.667245,0.602768,0.556823,0.521572/)
-XMASS_SPEED_LKT1D(84,1:4) = (/1.414226,1.228653,1.104546,1.013556/)
-
-XNUMB_SPEED_LKT1D(85,1:4) = (/0.676368,0.610543,0.563703,0.527802/)
-XMASS_SPEED_LKT1D(85,1:4) = (/1.428304,1.240152,1.114451,1.022357/)
-
-XNUMB_SPEED_LKT1D(86,1:4) = (/0.685464,0.618287,0.570552,0.534001/)
-XMASS_SPEED_LKT1D(86,1:4) = (/1.442280,1.251566,1.124281,1.031092/)
-
-XNUMB_SPEED_LKT1D(87,1:4) = (/0.694531,0.626002,0.577371,0.540170/)
-XMASS_SPEED_LKT1D(87,1:4) = (/1.456158,1.262896,1.134038,1.039760/)
-
-XNUMB_SPEED_LKT1D(88,1:4) = (/0.703570,0.633686,0.584159,0.546309/)
-XMASS_SPEED_LKT1D(88,1:4) = (/1.469938,1.274143,1.143723,1.048363/)
-
-XNUMB_SPEED_LKT1D(89,1:4) = (/0.712581,0.641340,0.590918,0.552418/)
-XMASS_SPEED_LKT1D(89,1:4) = (/1.483622,1.285309,1.153336,1.056903/)
-
-XNUMB_SPEED_LKT1D(90,1:4) = (/0.721563,0.648964,0.597645,0.558498/)
-XMASS_SPEED_LKT1D(90,1:4) = (/1.497210,1.296395,1.162879,1.065381/)
-
-XNUMB_SPEED_LKT1D(91,1:4) = (/0.730516,0.656557,0.604343,0.564548/)
-XMASS_SPEED_LKT1D(91,1:4) = (/1.510705,1.307402,1.172354,1.073797/)
-
-XNUMB_SPEED_LKT1D(92,1:4) = (/0.739440,0.664119,0.611010,0.570569/)
-XMASS_SPEED_LKT1D(92,1:4) = (/1.524108,1.318331,1.181761,1.082152/)
-
-XNUMB_SPEED_LKT1D(93,1:4) = (/0.748334,0.671651,0.617648,0.576560/)
-XMASS_SPEED_LKT1D(93,1:4) = (/1.537419,1.329184,1.191101,1.090449/)
-
-XNUMB_SPEED_LKT1D(94,1:4) = (/0.757199,0.679153,0.624255,0.582523/)
-XMASS_SPEED_LKT1D(94,1:4) = (/1.550640,1.339962,1.200377,1.098687/)
-
-XNUMB_SPEED_LKT1D(95,1:4) = (/0.766034,0.686624,0.630833,0.588456/)
-XMASS_SPEED_LKT1D(95,1:4) = (/1.563773,1.350665,1.209587,1.106867/)
-
-XNUMB_SPEED_LKT1D(96,1:4) = (/0.774839,0.694064,0.637380,0.594361/)
-XMASS_SPEED_LKT1D(96,1:4) = (/1.576819,1.361296,1.218735,1.114992/)
-
-XNUMB_SPEED_LKT1D(97,1:4) = (/0.783614,0.701474,0.643898,0.600238/)
-XMASS_SPEED_LKT1D(97,1:4) = (/1.589779,1.371855,1.227820,1.123060/)
-
-XNUMB_SPEED_LKT1D(98,1:4) = (/0.792358,0.708854,0.650387,0.606085/)
-XMASS_SPEED_LKT1D(98,1:4) = (/1.602653,1.382343,1.236844,1.131075/)
-
-XNUMB_SPEED_LKT1D(99,1:4) = (/0.801073,0.716204,0.656846,0.611905/)
-XMASS_SPEED_LKT1D(99,1:4) = (/1.615444,1.392761,1.245807,1.139036/)
-
-XNUMB_SPEED_LKT1D(100,1:4) = (/0.809757,0.723523,0.663276,0.617697/)
-XMASS_SPEED_LKT1D(100,1:4) = (/1.628153,1.403111,1.254711,1.146944/)
-
-XNUMB_SPEED_LKT1D(101,1:4) = (/0.818411,0.730812,0.669677,0.623461/)
-XMASS_SPEED_LKT1D(101,1:4) = (/1.640780,1.413394,1.263557,1.154800/)
-
-XNUMB_SPEED_LKT1D(102,1:4) = (/0.827035,0.738071,0.676049,0.629197/)
-XMASS_SPEED_LKT1D(102,1:4) = (/1.653326,1.423610,1.272346,1.162605/)
-
-XNUMB_SPEED_LKT1D(103,1:4) = (/0.835628,0.745300,0.682393,0.634906/)
-XMASS_SPEED_LKT1D(103,1:4) = (/1.665794,1.433760,1.281077,1.170360/)
-
-XNUMB_SPEED_LKT1D(104,1:4) = (/0.844190,0.752499,0.688707,0.640588/)
-XMASS_SPEED_LKT1D(104,1:4) = (/1.678183,1.443846,1.289753,1.178066/)
-
-XNUMB_SPEED_LKT1D(105,1:4) = (/0.852722,0.759668,0.694994,0.646243/)
-XMASS_SPEED_LKT1D(105,1:4) = (/1.690496,1.453868,1.298375,1.185723/)
-
-XNUMB_SPEED_LKT1D(106,1:4) = (/0.861224,0.766808,0.701252,0.651870/)
-XMASS_SPEED_LKT1D(106,1:4) = (/1.702733,1.463828,1.306942,1.193332/)
-
-XNUMB_SPEED_LKT1D(107,1:4) = (/0.869695,0.773918,0.707481,0.657472/)
-XMASS_SPEED_LKT1D(107,1:4) = (/1.714894,1.473726,1.315456,1.200894/)
-
-XNUMB_SPEED_LKT1D(108,1:4) = (/0.878135,0.780999,0.713683,0.663046/)
-XMASS_SPEED_LKT1D(108,1:4) = (/1.726982,1.483564,1.323918,1.208409/)
-
-XNUMB_SPEED_LKT1D(109,1:4) = (/0.886545,0.788050,0.719858,0.668595/)
-XMASS_SPEED_LKT1D(109,1:4) = (/1.738997,1.493341,1.332329,1.215879/)
-
-XNUMB_SPEED_LKT1D(110,1:4) = (/0.894925,0.795072,0.726004,0.674117/)
-XMASS_SPEED_LKT1D(110,1:4) = (/1.750941,1.503059,1.340688,1.223305/)
-
-XNUMB_SPEED_LKT1D(111,1:4) = (/0.903274,0.802066,0.732124,0.679614/)
-XMASS_SPEED_LKT1D(111,1:4) = (/1.762813,1.512720,1.348998,1.230685/)
-
-XNUMB_SPEED_LKT1D(112,1:4) = (/0.911593,0.809030,0.738216,0.685085/)
-XMASS_SPEED_LKT1D(112,1:4) = (/1.774616,1.522323,1.357258,1.238023/)
-
-XNUMB_SPEED_LKT1D(113,1:4) = (/0.919882,0.815966,0.744281,0.690531/)
-XMASS_SPEED_LKT1D(113,1:4) = (/1.786349,1.531869,1.365470,1.245317/)
-
-XNUMB_SPEED_LKT1D(114,1:4) = (/0.928140,0.822873,0.750319,0.695951/)
-XMASS_SPEED_LKT1D(114,1:4) = (/1.798015,1.541360,1.373634,1.252569/)
-
-XNUMB_SPEED_LKT1D(115,1:4) = (/0.936369,0.829751,0.756331,0.701346/)
-XMASS_SPEED_LKT1D(115,1:4) = (/1.809613,1.550796,1.381751,1.259780/)
-
-XNUMB_SPEED_LKT1D(116,1:4) = (/0.944567,0.836601,0.762316,0.706717/)
-XMASS_SPEED_LKT1D(116,1:4) = (/1.821145,1.560178,1.389821,1.266949/)
-
-XNUMB_SPEED_LKT1D(117,1:4) = (/0.952735,0.843423,0.768275,0.712063/)
-XMASS_SPEED_LKT1D(117,1:4) = (/1.832612,1.569506,1.397846,1.274078/)
-
-XNUMB_SPEED_LKT1D(118,1:4) = (/0.960874,0.850218,0.774208,0.717385/)
-XMASS_SPEED_LKT1D(118,1:4) = (/1.844014,1.578782,1.405825,1.281167/)
-
-XNUMB_SPEED_LKT1D(119,1:4) = (/0.968982,0.856984,0.780115,0.722682/)
-XMASS_SPEED_LKT1D(119,1:4) = (/1.855353,1.588006,1.413760,1.288217/)
-
-XNUMB_SPEED_LKT1D(120,1:4) = (/0.977061,0.863722,0.785996,0.727955/)
-XMASS_SPEED_LKT1D(120,1:4) = (/1.866629,1.597178,1.421652,1.295229/)
-
-XNUMB_SPEED_LKT1D(121,1:4) = (/0.985111,0.870434,0.791852,0.733205/)
-XMASS_SPEED_LKT1D(121,1:4) = (/1.877843,1.606301,1.429500,1.302202/)
-
-XNUMB_SPEED_LKT1D(122,1:4) = (/0.993131,0.877117,0.797682,0.738431/)
-XMASS_SPEED_LKT1D(122,1:4) = (/1.888996,1.615373,1.437305,1.309137/)
-
-XNUMB_SPEED_LKT1D(123,1:4) = (/1.001121,0.883774,0.803487,0.743633/)
-XMASS_SPEED_LKT1D(123,1:4) = (/1.900088,1.624397,1.445069,1.316036/)
-
-XNUMB_SPEED_LKT1D(124,1:4) = (/1.009083,0.890404,0.809268,0.748812/)
-XMASS_SPEED_LKT1D(124,1:4) = (/1.911121,1.633372,1.452791,1.322898/)
-
-XNUMB_SPEED_LKT1D(125,1:4) = (/1.017015,0.897007,0.815023,0.753969/)
-XMASS_SPEED_LKT1D(125,1:4) = (/1.922095,1.642299,1.460472,1.329723/)
-
-XNUMB_SPEED_LKT1D(126,1:4) = (/1.024918,0.903583,0.820754,0.759102/)
-XMASS_SPEED_LKT1D(126,1:4) = (/1.933011,1.651179,1.468113,1.336514/)
-
-XNUMB_SPEED_LKT1D(127,1:4) = (/1.032792,0.910133,0.826461,0.764213/)
-XMASS_SPEED_LKT1D(127,1:4) = (/1.943870,1.660012,1.475715,1.343269/)
-
-XNUMB_SPEED_LKT1D(128,1:4) = (/1.040638,0.916656,0.832143,0.769302/)
-XMASS_SPEED_LKT1D(128,1:4) = (/1.954672,1.668800,1.483277,1.349990/)
-
-XNUMB_SPEED_LKT1D(129,1:4) = (/1.048455,0.923153,0.837802,0.774368/)
-XMASS_SPEED_LKT1D(129,1:4) = (/1.965418,1.677543,1.490800,1.356677/)
-
-XNUMB_SPEED_LKT1D(130,1:4) = (/1.056243,0.929625,0.843437,0.779412/)
-XMASS_SPEED_LKT1D(130,1:4) = (/1.976110,1.686240,1.498285,1.363330/)
-
-XNUMB_SPEED_LKT1D(131,1:4) = (/1.064003,0.936070,0.849048,0.784434/)
-XMASS_SPEED_LKT1D(131,1:4) = (/1.986746,1.694894,1.505733,1.369950/)
-
-XNUMB_SPEED_LKT1D(132,1:4) = (/1.071735,0.942490,0.854635,0.789435/)
-XMASS_SPEED_LKT1D(132,1:4) = (/1.997329,1.703504,1.513143,1.376537/)
-
-XNUMB_SPEED_LKT1D(133,1:4) = (/1.079438,0.948885,0.860199,0.794414/)
-XMASS_SPEED_LKT1D(133,1:4) = (/2.007859,1.712071,1.520517,1.383092/)
-
-XNUMB_SPEED_LKT1D(134,1:4) = (/1.087114,0.955254,0.865741,0.799372/)
-XMASS_SPEED_LKT1D(134,1:4) = (/2.018336,1.720596,1.527855,1.389615/)
-
-XNUMB_SPEED_LKT1D(135,1:4) = (/1.094762,0.961598,0.871259,0.804309/)
-XMASS_SPEED_LKT1D(135,1:4) = (/2.028761,1.729079,1.535157,1.396106/)
-
-XNUMB_SPEED_LKT1D(136,1:4) = (/1.102382,0.967917,0.876755,0.809225/)
-XMASS_SPEED_LKT1D(136,1:4) = (/2.039136,1.737520,1.542423,1.402567/)
-
-XNUMB_SPEED_LKT1D(137,1:4) = (/1.109975,0.974212,0.882228,0.814120/)
-XMASS_SPEED_LKT1D(137,1:4) = (/2.049459,1.745921,1.549655,1.408996/)
-
-XNUMB_SPEED_LKT1D(138,1:4) = (/1.117540,0.980482,0.887679,0.818994/)
-XMASS_SPEED_LKT1D(138,1:4) = (/2.059733,1.754281,1.556852,1.415396/)
-
-XNUMB_SPEED_LKT1D(139,1:4) = (/1.125079,0.986727,0.893107,0.823848/)
-XMASS_SPEED_LKT1D(139,1:4) = (/2.069957,1.762601,1.564016,1.421766/)
-
-XNUMB_SPEED_LKT1D(140,1:4) = (/1.132590,0.992948,0.898514,0.828682/)
-XMASS_SPEED_LKT1D(140,1:4) = (/2.080133,1.770882,1.571146,1.428106/)
-
-XNUMB_SPEED_LKT1D(141,1:4) = (/1.140074,0.999145,0.903899,0.833496/)
-XMASS_SPEED_LKT1D(141,1:4) = (/2.090260,1.779125,1.578243,1.434417/)
-
-XNUMB_SPEED_LKT1D(142,1:4) = (/1.147531,1.005319,0.909262,0.838290/)
-XMASS_SPEED_LKT1D(142,1:4) = (/2.100340,1.787329,1.585307,1.440699/)
-
-XNUMB_SPEED_LKT1D(143,1:4) = (/1.154962,1.011468,0.914604,0.843064/)
-XMASS_SPEED_LKT1D(143,1:4) = (/2.110373,1.795495,1.592339,1.446953/)
-
-XNUMB_SPEED_LKT1D(144,1:4) = (/1.162366,1.017594,0.919924,0.847819/)
-XMASS_SPEED_LKT1D(144,1:4) = (/2.120359,1.803623,1.599339,1.453179/)
-
-XNUMB_SPEED_LKT1D(145,1:4) = (/1.169744,1.023697,0.925223,0.852555/)
-XMASS_SPEED_LKT1D(145,1:4) = (/2.130300,1.811715,1.606307,1.459377/)
-
-XNUMB_SPEED_LKT1D(146,1:4) = (/1.177096,1.029777,0.930502,0.857271/)
-XMASS_SPEED_LKT1D(146,1:4) = (/2.140195,1.819770,1.613245,1.465548/)
-
-XNUMB_SPEED_LKT1D(147,1:4) = (/1.184422,1.035833,0.935759,0.861968/)
-XMASS_SPEED_LKT1D(147,1:4) = (/2.150045,1.827789,1.620152,1.471691/)
-
-XNUMB_SPEED_LKT1D(148,1:4) = (/1.191722,1.041867,0.940996,0.866647/)
-XMASS_SPEED_LKT1D(148,1:4) = (/2.159851,1.835772,1.627028,1.477808/)
-
-XNUMB_SPEED_LKT1D(149,1:4) = (/1.198996,1.047877,0.946213,0.871306/)
-XMASS_SPEED_LKT1D(149,1:4) = (/2.169612,1.843720,1.633874,1.483899/)
-
-XNUMB_SPEED_LKT1D(150,1:4) = (/1.206245,1.053866,0.951409,0.875947/)
-XMASS_SPEED_LKT1D(150,1:4) = (/2.179331,1.851633,1.640691,1.489964/)
-
-XNUMB_SPEED_LKT1D(151,1:4) = (/1.213468,1.059832,0.956585,0.880570/)
-XMASS_SPEED_LKT1D(151,1:4) = (/2.189007,1.859512,1.647479,1.496002/)
-
-XNUMB_SPEED_LKT1D(152,1:4) = (/1.220666,1.065775,0.961741,0.885175/)
-XMASS_SPEED_LKT1D(152,1:4) = (/2.198640,1.867357,1.654237,1.502016/)
-
-XNUMB_SPEED_LKT1D(153,1:4) = (/1.227839,1.071697,0.966877,0.889761/)
-XMASS_SPEED_LKT1D(153,1:4) = (/2.208232,1.875168,1.660967,1.508004/)
-
-XNUMB_SPEED_LKT1D(154,1:4) = (/1.234987,1.077597,0.971994,0.894330/)
-XMASS_SPEED_LKT1D(154,1:4) = (/2.217782,1.882946,1.667668,1.513967/)
-
-XNUMB_SPEED_LKT1D(155,1:4) = (/1.242110,1.083475,0.977091,0.898881/)
-XMASS_SPEED_LKT1D(155,1:4) = (/2.227291,1.890690,1.674342,1.519906/)
-
-XNUMB_SPEED_LKT1D(156,1:4) = (/1.249209,1.089331,0.982169,0.903414/)
-XMASS_SPEED_LKT1D(156,1:4) = (/2.236759,1.898403,1.680988,1.525820/)
-
-XNUMB_SPEED_LKT1D(157,1:4) = (/1.256283,1.095166,0.987228,0.907930/)
-XMASS_SPEED_LKT1D(157,1:4) = (/2.246188,1.906083,1.687607,1.531710/)
-
-XNUMB_SPEED_LKT1D(158,1:4) = (/1.263333,1.100980,0.992267,0.912428/)
-XMASS_SPEED_LKT1D(158,1:4) = (/2.255577,1.913731,1.694198,1.537577/)
-
-XNUMB_SPEED_LKT1D(159,1:4) = (/1.270359,1.106773,0.997288,0.916909/)
-XMASS_SPEED_LKT1D(159,1:4) = (/2.264927,1.921348,1.700763,1.543420/)
-
-XNUMB_SPEED_LKT1D(160,1:4) = (/1.277360,1.112544,1.002290,0.921374/)
-XMASS_SPEED_LKT1D(160,1:4) = (/2.274238,1.928934,1.707301,1.549240/)
-
-XNUMB_SPEED_LKT1D(161,1:4) = (/1.284338,1.118295,1.007274,0.925821/)
-XMASS_SPEED_LKT1D(161,1:4) = (/2.283510,1.936489,1.713814,1.555036/)
-
-XNUMB_SPEED_LKT1D(162,1:4) = (/1.291292,1.124026,1.012239,0.930252/)
-XMASS_SPEED_LKT1D(162,1:4) = (/2.292745,1.944014,1.720300,1.560810/)
-
-XNUMB_SPEED_LKT1D(163,1:4) = (/1.298222,1.129736,1.017186,0.934666/)
-XMASS_SPEED_LKT1D(163,1:4) = (/2.301942,1.951509,1.726761,1.566562/)
-
-XNUMB_SPEED_LKT1D(164,1:4) = (/1.305129,1.135425,1.022115,0.939064/)
-XMASS_SPEED_LKT1D(164,1:4) = (/2.311102,1.958974,1.733197,1.572291/)
-
-XNUMB_SPEED_LKT1D(165,1:4) = (/1.312013,1.141094,1.027026,0.943445/)
-XMASS_SPEED_LKT1D(165,1:4) = (/2.320226,1.966409,1.739607,1.577999/)
-
-XNUMB_SPEED_LKT1D(166,1:4) = (/1.318874,1.146744,1.031919,0.947810/)
-XMASS_SPEED_LKT1D(166,1:4) = (/2.329313,1.973815,1.745993,1.583684/)
-
-XNUMB_SPEED_LKT1D(167,1:4) = (/1.325711,1.152373,1.036794,0.952159/)
-XMASS_SPEED_LKT1D(167,1:4) = (/2.338364,1.981193,1.752354,1.589348/)
-
-XNUMB_SPEED_LKT1D(168,1:4) = (/1.332526,1.157983,1.041652,0.956493/)
-XMASS_SPEED_LKT1D(168,1:4) = (/2.347379,1.988542,1.758692,1.594991/)
-
-XNUMB_SPEED_LKT1D(169,1:4) = (/1.339318,1.163573,1.046492,0.960810/)
-XMASS_SPEED_LKT1D(169,1:4) = (/2.356360,1.995863,1.765005,1.600613/)
-
-XNUMB_SPEED_LKT1D(170,1:4) = (/1.346088,1.169143,1.051315,0.965112/)
-XMASS_SPEED_LKT1D(170,1:4) = (/2.365305,2.003155,1.771294,1.606214/)
-
-XNUMB_SPEED_LKT1D(171,1:4) = (/1.352835,1.174695,1.056121,0.969398/)
-XMASS_SPEED_LKT1D(171,1:4) = (/2.374216,2.010421,1.777560,1.611794/)
-
-XNUMB_SPEED_LKT1D(172,1:4) = (/1.359560,1.180227,1.060910,0.973669/)
-XMASS_SPEED_LKT1D(172,1:4) = (/2.383093,2.017658,1.783803,1.617354/)
-
-XNUMB_SPEED_LKT1D(173,1:4) = (/1.366263,1.185740,1.065682,0.977925/)
-XMASS_SPEED_LKT1D(173,1:4) = (/2.391936,2.024869,1.790023,1.622893/)
-
-XNUMB_SPEED_LKT1D(174,1:4) = (/1.372943,1.191234,1.070437,0.982166/)
-XMASS_SPEED_LKT1D(174,1:4) = (/2.400745,2.032053,1.796220,1.628413/)
-
-XNUMB_SPEED_LKT1D(175,1:4) = (/1.379602,1.196710,1.075176,0.986391/)
-XMASS_SPEED_LKT1D(175,1:4) = (/2.409522,2.039211,1.802394,1.633912/)
-
-XNUMB_SPEED_LKT1D(176,1:4) = (/1.386240,1.202167,1.079898,0.990602/)
-XMASS_SPEED_LKT1D(176,1:4) = (/2.418265,2.046342,1.808546,1.639393/)
-
-XNUMB_SPEED_LKT1D(177,1:4) = (/1.392856,1.207605,1.084604,0.994798/)
-XMASS_SPEED_LKT1D(177,1:4) = (/2.426976,2.053447,1.814676,1.644853/)
-
-XNUMB_SPEED_LKT1D(178,1:4) = (/1.399450,1.213025,1.089294,0.998979/)
-XMASS_SPEED_LKT1D(178,1:4) = (/2.435655,2.060527,1.820784,1.650295/)
-
-XNUMB_SPEED_LKT1D(179,1:4) = (/1.406023,1.218427,1.093968,1.003146/)
-XMASS_SPEED_LKT1D(179,1:4) = (/2.444302,2.067581,1.826871,1.655717/)
-
-XNUMB_SPEED_LKT1D(180,1:4) = (/1.412575,1.223811,1.098625,1.007298/)
-XMASS_SPEED_LKT1D(180,1:4) = (/2.452918,2.074610,1.832936,1.661121/)
-
-XNUMB_SPEED_LKT1D(181,1:4) = (/1.419107,1.229177,1.103267,1.011437/)
-XMASS_SPEED_LKT1D(181,1:4) = (/2.461502,2.081613,1.838980,1.666506/)
-
-XNUMB_SPEED_LKT1D(182,1:4) = (/1.425617,1.234525,1.107893,1.015560/)
-XMASS_SPEED_LKT1D(182,1:4) = (/2.470056,2.088592,1.845003,1.671872/)
-
-XNUMB_SPEED_LKT1D(183,1:4) = (/1.432106,1.239856,1.112504,1.019670/)
-XMASS_SPEED_LKT1D(183,1:4) = (/2.478579,2.095547,1.851005,1.677220/)
-
-XNUMB_SPEED_LKT1D(184,1:4) = (/1.438575,1.245168,1.117099,1.023766/)
-XMASS_SPEED_LKT1D(184,1:4) = (/2.487071,2.102477,1.856986,1.682550/)
-
-XNUMB_SPEED_LKT1D(185,1:4) = (/1.445024,1.250464,1.121678,1.027848/)
-XMASS_SPEED_LKT1D(185,1:4) = (/2.495533,2.109383,1.862948,1.687863/)
-
-XNUMB_SPEED_LKT1D(186,1:4) = (/1.451452,1.255742,1.126243,1.031917/)
-XMASS_SPEED_LKT1D(186,1:4) = (/2.503966,2.116266,1.868888,1.693157/)
-
-XNUMB_SPEED_LKT1D(187,1:4) = (/1.457860,1.261003,1.130792,1.035971/)
-XMASS_SPEED_LKT1D(187,1:4) = (/2.512369,2.123125,1.874809,1.698433/)
-
-XNUMB_SPEED_LKT1D(188,1:4) = (/1.464248,1.266247,1.135327,1.040013/)
-XMASS_SPEED_LKT1D(188,1:4) = (/2.520743,2.129960,1.880710,1.703693/)
-
-XNUMB_SPEED_LKT1D(189,1:4) = (/1.470616,1.271474,1.139846,1.044041/)
-XMASS_SPEED_LKT1D(189,1:4) = (/2.529088,2.136772,1.886591,1.708934/)
-
-XNUMB_SPEED_LKT1D(190,1:4) = (/1.476965,1.276684,1.144351,1.048055/)
-XMASS_SPEED_LKT1D(190,1:4) = (/2.537404,2.143562,1.892453,1.714159/)
-
-XNUMB_SPEED_LKT1D(191,1:4) = (/1.483293,1.281878,1.148841,1.052056/)
-XMASS_SPEED_LKT1D(191,1:4) = (/2.545692,2.150328,1.898295,1.719367/)
-
-XNUMB_SPEED_LKT1D(192,1:4) = (/1.489603,1.287055,1.153316,1.056045/)
-XMASS_SPEED_LKT1D(192,1:4) = (/2.553951,2.157072,1.904119,1.724558/)
-
-XNUMB_SPEED_LKT1D(193,1:4) = (/1.495892,1.292215,1.157777,1.060020/)
-XMASS_SPEED_LKT1D(193,1:4) = (/2.562183,2.163794,1.909923,1.729732/)
-
-XNUMB_SPEED_LKT1D(194,1:4) = (/1.502163,1.297360,1.162224,1.063982/)
-XMASS_SPEED_LKT1D(194,1:4) = (/2.570387,2.170493,1.915708,1.734889/)
-
-XNUMB_SPEED_LKT1D(195,1:4) = (/1.508414,1.302488,1.166656,1.067932/)
-XMASS_SPEED_LKT1D(195,1:4) = (/2.578563,2.177171,1.921475,1.740031/)
-
-XNUMB_SPEED_LKT1D(196,1:4) = (/1.514647,1.307599,1.171075,1.071869/)
-XMASS_SPEED_LKT1D(196,1:4) = (/2.586713,2.183826,1.927223,1.745156/)
-
-ENDIF
-
-IF (LHOOK) CALL DR_HOOK('BLOWSNW_SEDIM_LKT1D_SET',0,ZHOOK_HANDLE)
-
-END SUBROUTINE BLOWSNW_SEDIM_LKT1D_SET
-
-END MODULE MODE_BLOWSNW_SEDIM_LKT1D
diff --git a/src/SURFEX/mode_blowsnw_surf.f90 b/src/SURFEX/mode_blowsnw_surf.f90
deleted file mode 100644
index 996b5b97c..000000000
--- a/src/SURFEX/mode_blowsnw_surf.f90
+++ /dev/null
@@ -1,118 +0,0 @@
-!! ########################
-MODULE MODE_BLOWSNW_SURF
-!! ########################
-!!
-
- USE MODD_BLOWSNW_SURF
- USE MODD_CSTS
- USE MODD_BLOWSNW_n
-
-
- IMPLICIT NONE
- PUBLIC
-
-CONTAINS
-
-SUBROUTINE SNOWMOMENT2SIZE( &
- PSVT & !I [XX/m3] input scalar variables (moment of distribution)
- , PBETA1D & !O [m] scale factor of blowing snow particles distribution (gamma distribution)
- , PRG1D & !O [m] number mean radius of blowing snow particles distribution
- )
- !! ############################################################
- !!
- !!
- !! PURPOSE
- !! -------
- !! Translate the two moments M0 and M3 into
- !! Values which can be understood more easily (R, beta)
- !! At this point, M3 is in kg_{snow}/m3, M0 in #/m3
- !!
- !!
- !! REFERENCE
- !! ---------
- !! Based on mode_dst_surf.f90 (Tulet et al)
- !!
- !! AUTHOR
- !! ------
- !! Vincent Vionnet
- !!
- !! MODIFICATIONS
- !! -------------
- !!
- !! EXTERNAL
- !! --------
- !! None
- !!
- USE YOMHOOK ,ONLY : LHOOK, DR_HOOK
- USE PARKIND1 ,ONLY : JPRB
-
- IMPLICIT NONE
- !!
- !-------------------------------------------------------------------------------
- !
- !* 0. DECLARATIONS
- ! ------------
- !
- !* 0.1 declarations of arguments
- !
- !INPUT
- REAL, DIMENSION(:,:,:), INTENT(IN) :: PSVT !I [ __ /m3] moments in surface units
-
- !OUTPUT
- REAL, DIMENSION(:,:), OPTIONAL, INTENT(OUT) :: PBETA1D !O [m] scale factor deviation
- REAL, DIMENSION(:,:), OPTIONAL, INTENT(OUT) :: PRG1D !O [m] number median diameter
- !
- !
- !* 0.2 declarations local variables
- !
- REAL,DIMENSION(:,:,:), ALLOCATABLE :: ZSV ! [snow particles moment concentration]
- REAL,DIMENSION(:,:), ALLOCATABLE :: ZBETA ! [-] standard deviation
- REAL,DIMENSION(:,:), ALLOCATABLE :: ZRG ! [um] number median diameter
- INTEGER :: JK !Loop index
- REAL, DIMENSION(2) :: ZPMIN
- REAL :: ZRGMIN
- REAL(KIND=JPRB) :: ZHOOK_HANDLE
-
-IF (LHOOK) CALL DR_HOOK('MODE_BLOWSNW_SURF:SNOWMOMENT2SIZE',0,ZHOOK_HANDLE)
-
- ALLOCATE (ZBETA(SIZE(PSVT,1),SIZE(PSVT,2)))
- ALLOCATE (ZRG(SIZE(PSVT,1),SIZE(PSVT,2)))
- ALLOCATE (ZSV(SIZE(PSVT,1), SIZE(PSVT,2),SIZE(PSVT,3)))
-
- !Save the moments in a local array
- ZSV(:,:,:) = MAX(PSVT(:,:,:), 1E-80)
-
-!Get minimum values possible
- ZPMIN(1) = XN0MIN_SNW
- ZRGMIN = XINIRADIUS_SNW
- ZPMIN(2) = 4*XPI*XRHOLI/3*(ZRGMIN/XEMIALPHA_SNW)**3.*(XEMIALPHA_SNW+2)*(XEMIALPHA_SNW+1)*XEMIALPHA_SNW*XN0MIN_SNW
- ZSV(:,:,1)=MAX(PSVT(:,:,1),ZPMIN(1))
- ZSV(:,:,2)=MAX(PSVT(:,:,2),ZPMIN(2))
-
- DO JK=1,SIZE(PSVT,2)
- ZBETA(:,JK)=((3*ZSV(:,JK,2))/(4*XPI*XRHOLI*(XEMIALPHA_SNW+2)*(XEMIALPHA_SNW+1)*XEMIALPHA_SNW*ZSV(:,JK,1)))**(1./3.)
-
- ZRG(:,JK) = ZBETA(:,JK)*XEMIALPHA_SNW
-
- END DO
-
- !Give the beta-values to the passed array
- IF(PRESENT(PBETA1D))THEN
- PBETA1D(:,:) = ZBETA(:,:)
- ENDIF
-
- !Get the number median radius
- IF(PRESENT(PRG1D))THEN
- PRG1D(:,:)= ZRG(:,:)
- ENDIF
-
- DEALLOCATE(ZSV)
- DEALLOCATE(ZRG)
- DEALLOCATE(ZBETA)
-
-IF (LHOOK) CALL DR_HOOK('MODE_BLOWSNW_SURF:SNOWMOMENT2SIZE',1,ZHOOK_HANDLE)
-
-
-END SUBROUTINE SNOWMOMENT2SIZE
-
-END MODULE MODE_BLOWSNW_SURF
diff --git a/src/SURFEX/modn_blowsnw.f90 b/src/SURFEX/modn_blowsnw.f90
deleted file mode 100644
index d03644939..000000000
--- a/src/SURFEX/modn_blowsnw.f90
+++ /dev/null
@@ -1,42 +0,0 @@
-!!
-!! #####################
- MODULE MODN_BLOWSNW
-!! #####################
-!!
-!!*** *MODN_SNW*
-!!
-!! PURPOSE
-!! -------
-! Namelist for snow drift scheme
-!!
-!!** AUTHOR
-!! ------
-!! V. Vionnet *CNRM*
-!
-!! MODIFICATIONS
-!! -------------
-!! Original 03/2010
-!!
-!! IMPLICIT ARGUMENTS
-!! ------------------
-USE MODD_BLOWSNW_SURF, ONLY : LSNOW_SALT, XCONC_SALT, XSNOW_PRECIP, LSNOW_PRECIP, &
- LSNOW_WIND, XSNOW_ROUGHNESS, XSNOW_BUOYANCY, &
- XEMIRADIUS_SNW,XEMIALPHA_SNW,CSNOW_SALT, &
- CSNOW_SEDIM,LBLOWSNW_CANOSUBL,LBLOWSNW_CANODIAG, &
- XRSNOW_SBL,LBLOWSNW_ADV
-
-!!
-!-----------------------------------------------------------------------------
-!
-!* 0. DECLARATIONS
-! -----------------
-IMPLICIT NONE
-SAVE
-NAMELIST /NAM_SURF_BLOWSNW/ &
- LSNOW_SALT, XCONC_SALT, XSNOW_PRECIP, LSNOW_PRECIP, LSNOW_WIND, & !Parameterization type
- XSNOW_ROUGHNESS, XSNOW_BUOYANCY,XEMIRADIUS_SNW,XEMIALPHA_SNW , &
- CSNOW_SALT, CSNOW_SEDIM,LBLOWSNW_CANOSUBL,LBLOWSNW_CANODIAG, &
- XRSNOW_SBL,LBLOWSNW_ADV
-
-!
-END MODULE MODN_BLOWSNW
diff --git a/src/SURFEX/snowpack_evol.f90 b/src/SURFEX/snowpack_evol.f90
deleted file mode 100644
index 370879352..000000000
--- a/src/SURFEX/snowpack_evol.f90
+++ /dev/null
@@ -1,744 +0,0 @@
-SUBROUTINE SNOWPACK_EVOL(HSNOWRES,PBLOWSNW_FLUX,PSNOWHEAT,PSNOWSWE,PSNOWRHO, &
- PSNOWGRAN1,PSNOWGRAN2,PSNOWHIST,PSNOWAGE, &
- PTSTEP,PRHOA,PTA,PBLOWSNW_CONC, &
- PVMOD,PQA,PPS,PUREF,PEXNS,PDIRCOSZW,PZREF, &
- PZ0EFF,PZ0H,PBLOWSNW_DEP,PTG)
-! ######################################################################
-! ##########################################################################
-!
-!!
-!! PURPOSE
-!! -------
-!! The purpose of this routine is to simulate the mass exchange between the
-!! snowpack and the atmosphere when coupled to Meso-NH. A net mass flux flux is computed.
-!! - If erosion occurs, it is simumated by this routine.
-!! - If deposition occurs, the deposited mass flux is computed, sent to Crocus
-!! and is added to the snowfall
-!!
-!!
-!!** METHOD
-!! ------
-!!
-!! EXTERNAL
-!! --------
-!!
-!! IMPLICIT ARGUMENTS
-!! ------------------
-!!
-!! REFERENCE
-!! ---------
-!!
-!! Vionnet et al, The Cryosphere, 2014
-!!
-!! AUTHOR
-!! ------
-!! V. Vionnet * CNRM/GAME*
-!!
-!! MODIFICATIONS
-!! -------------
-!! Original 10/10/11
-!
-!* 0. DECLARATIONS
-
-USE MODD_CSTS, ONLY : XLMTT, XTT, XCI,XRHOLW
-USE MODD_SNOW_PAR, ONLY : X_RI_MAX
-USE MODD_BLOWSNW_SURF
-USE MODD_BLOWSNW_n
-USE MODD_SNOW_METAMO
-
-USE MODE_SNOW3L
-USE MODE_THERMOS
-
-USE MODI_SURFACE_RI
-USE MODI_SURFACE_CD
-
-USE YOMHOOK ,ONLY : LHOOK, DR_HOOK
-USE PARKIND1 ,ONLY : JPRB
-
-IMPLICIT NONE
-! ===============================================================
-!
-! 0.1 declarations of arguments
-!
-! ===============================================================
-!
- CHARACTER(LEN=*), INTENT(IN) :: HSNOWRES
-! HSNOWRES = ISBA-SNOW3L turbulant exchange option
-! 'DEF' = Default: Louis (ISBA: Noilhan and Mahfouf 1996)
-! 'RIL' = Limit Richarson number under very stable
-! conditions
-
- REAL, DIMENSION(:), INTENT(IN) :: PTA,PRHOA,PVMOD,PZ0EFF,PQA,PEXNS, &
- PUREF,PZREF,PDIRCOSZW,PPS,PZ0H,PTG
-! PTA = atmospheric temperature at level za (K)
-! PRHOA = air density (kg/m3)
-! PVMOD = modulus of the wind parallel to the orography (m/s)
-! PZ0EFF = roughness length for momentum
-! PQA = atmospheric specific humidity
-! at level za
-! PEXNS = Exner function
-! PUREF = reference height of the wind
-! PZREF = reference height of the first
-! PDIRCOSZW = Cosinus of the angle between the
-! normal to the surface and the vertical
-! atmospheric level
-! PPS = surface pressure
-! PZ0H = grid box average roughness length for heat
-! PTG = Surface soil temperature (effective
-! temperature of the layer lying below snowpack)
-
- REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWHEAT,PSNOWSWE,PSNOWRHO
-! PSNOWHEAT = Snow layer(s) heat content (J/m3)
-! PSNOWRHO = Snow layer(s) averaged density (kg/m3)
-! PSNOWSWE = Snow layer(s) liquid Water Equivalent (SWE:kg m-2)
-
- REAL, DIMENSION(:,:), INTENT(INOUT) :: PSNOWGRAN1,PSNOWGRAN2, &
- PSNOWHIST,PSNOWAGE
-! PSNOWGRAN1 = Snow layer(s) grain parameter 1
-! PSNOWGRAN2 = Snow layer(s) grain parameter 2
-! PSNOWHIST = Snow layer(s) grain historical parameter
-! PSNOWAGE = Snow layer(s) age
-
- REAL, DIMENSION(:,:) , INTENT(IN) :: PBLOWSNW_CONC
-! PBLOWSNW_CONC = blowing snow particle concentration
-! (1: number; 2: mass)
-
- REAL, DIMENSION(:,:) , INTENT(INOUT) :: PBLOWSNW_FLUX
-! PBLOWSNW_FLUX = Blowing snow particles flux
-! IN : sedimentation flux from previous time step (1,2)
-! contribution of saltation transport to snowpack mass balance (3)
-! 1: number (#/m2/s); 2: mass (kg/m2/s), 3: (kg/m2/s)
-! OUT: emitted turbulent flux + streamwise saltation flux
-! 1: number (#/m2/s); 2: mass (kg/m2/s), 3: kg/m/s
-!
- REAL, DIMENSION(:) , INTENT(OUT) :: PBLOWSNW_DEP
-! PBLOWSNW_DEP = deposion flux of blowing snow particles (sent to Crocus)
-!
- REAL, INTENT(IN) :: PTSTEP
-! PTSTEP = time step of the integration
-
-! ===============================================================
-!
-! 0.2 declarations of local variables
-!
-! ===============================================================
- INTEGER JWRK, JJ ! Loop counters
- INTEGER INLVLS ! maximum number of snow layers
- REAL SWE_MAX ! Maximum SWE that can be removed from a snow layer
- REAL SWE_DEP ! SWE which is deposited
- REAL ZMOB ! Mobility index of the eroded layer
-
-
- REAL, DIMENSION(SIZE(PRHOA)) :: T_EROSION ! Erosion time step
-
- INTEGER, DIMENSION(SIZE(PRHOA)) :: INLVLS_USE ! decrit nbre de couches effectives initiales
- INTEGER, DIMENSION(SIZE(PRHOA)) :: INLVLS_USE_TMP ! decrit nbre de couches effectives en cours d'�volution
-
- LOGICAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: LAYER_USE !
- LOGICAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: LAYER_USE_TMP !
-
- REAL , DIMENSION(SIZE(PRHOA),2) :: ZTURBFLUX ! blowing snow turbulent flux
- REAL , DIMENSION(SIZE(PRHOA),2) :: ZSEDFLUX ! blowing snow sedimentation flux
- REAL , DIMENSION(SIZE(PRHOA)) :: ZSALT_CONTR ! contribution of saltation transport to snowpack mass balance
- REAL , DIMENSION(SIZE(PRHOA),2) :: ZTOTFLUX_TURB ! total turbulent blowing snow flux
- REAL , DIMENSION(SIZE(PRHOA)) :: ZTOTFLUX_SALT ! total saltation mass flux
- REAL , DIMENSION(SIZE(PRHOA)) :: ZSALTFLUX ! streamwise saltation mass flux
- REAL , DIMENSION(SIZE(PRHOA)) :: ZEMIMOB ! emitted mobility index
- REAL , DIMENSION(SIZE(PRHOA)) :: ZDEPMOB ! deposited mobility index
- REAL , DIMENSION(SIZE(PRHOA)) :: ZUSTAR ! friction velocity
- REAL , DIMENSION(SIZE(PRHOA)) :: ZRISNOW ! Richardson number over snow
- REAL , DIMENSION(SIZE(PRHOA)) :: ZCDSNOW ! Drag coefficient for momentum over snow
-
-
- REAL, DIMENSION(SIZE(PRHOA)) :: ZEXNA, ZQSAT, ZCDN,ZTEMP_RI
-
- ! Temporary snow variables
- REAL , DIMENSION(SIZE(PRHOA,1)) :: ZSNOW ! total snow depth
-
- REAL,DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZSNOWGRAN1,ZSNOWGRAN2, &
- ZSNOWHIST,ZSNOWAGE
-
- REAL, DIMENSION(SIZE(PSNOWRHO,1),SIZE(PSNOWRHO,2)) :: ZSNOWHEAT,ZSNOWDZ, &
- ZSNOWSWE,ZSNOWRHO,ZSNOWTEMP,ZSNOWLIQ, &
- ZSCAP
-
- REAL(KIND=JPRB) :: ZHOOK_HANDLE
-! ===============================================================
-!
-! 1. Initialization of local variables
-!
-! ===============================================================
-!
-IF (LHOOK) CALL DR_HOOK('SNOWPACK_EVOL',0,ZHOOK_HANDLE)
-
-INLVLS = SIZE(PSNOWSWE(:,:),2)
-
-
-
-ZSEDFLUX(:,:)= -PBLOWSNW_FLUX(:,1:2)
-ZSALT_CONTR(:) = -PBLOWSNW_FLUX(:,3)
-
-ZTOTFLUX_TURB(:,:)=0.
-
-ZTOTFLUX_SALT(:)=0.
-ZSALTFLUX(:)=0.
-
-ZEMIMOB(:) =0.
-ZDEPMOB(:) =0.
-
-! Initialize snowpack characteristics
-ZSNOWGRAN1(:,:)=PSNOWGRAN1(:,:)
-ZSNOWGRAN2(:,:)=PSNOWGRAN2(:,:)
-ZSNOWHIST(:,:) =PSNOWHIST(:,:)
-ZSNOWAGE(:,:) =PSNOWAGE(:,:)
-ZSNOWHEAT(:,:) =PSNOWHEAT(:,:)
-ZSNOWRHO(:,:) =PSNOWRHO(:,:)
-ZSNOWSWE(:,:) =PSNOWSWE(:,:)
-
-
-! Initially : no deposition of blowing snow particles
-PBLOWSNW_DEP(:) = 0.
-
-
-WHERE (PSNOWSWE(:,:)>0.)
- LAYER_USE(:,:)=.TRUE.
- ZSNOWDZ(:,:) = PSNOWSWE(:,:)/PSNOWRHO(:,:)
-ELSEWHERE
- LAYER_USE(:,:)=.FALSE.
- ZSNOWDZ(:,:) =0.
-ENDWHERE
-
-ZSNOWTEMP(:,:)= 0.
-
-! Initialize temperature where snow is present
-WHERE (LAYER_USE(:,:))
-ZSCAP(:,:) = SNOW3LSCAP(PSNOWRHO)
-!
-ZSNOWTEMP(:,:) = XTT + ( ((PSNOWHEAT(:,:)/ZSNOWDZ(:,:)) &
- + XLMTT*PSNOWRHO(:,:))/ZSCAP(:,:) )
-!
-ZSNOWLIQ(:,:) = MAX(0.0,ZSNOWTEMP(:,:)-XTT)*ZSCAP(:,:)* &
- ZSNOWDZ(:,:)/(XLMTT*XRHOLW)
-!
-ZSNOWTEMP(:,:) = MIN(XTT,ZSNOWTEMP(:,:))
-END WHERE
-
-! Temperature used to compute Richardson number depending if the soil is
-! covered by snow or not.
-ZTEMP_RI(:) = 0.
-WHERE(ZSNOWTEMP(:,1)>0.)
- ZTEMP_RI(:) = ZSNOWTEMP(:,1)
-ELSEWHERE
- ZTEMP_RI(:) = PTG(:)
-END WHERE
-
-
-ZQSAT(:) = QSATI(ZTEMP_RI(:),PPS(:))
-ZEXNA(:) = PEXNS(:)
-!
-! Compute Richardson number
-!
-CALL SURFACE_RI(ZTEMP_RI(:), ZQSAT, PEXNS, ZEXNA, PTA, PQA, &
- PZREF, PUREF, PDIRCOSZW, PVMOD, ZRISNOW )
-!
-! Simple adaptation of method by Martin and Lejeune (1998)
-! to apply a lower limit to turbulent transfer coef
-! by defining a maximum Richarson number for stable
-! conditions:
-!
-IF(HSNOWRES=='RIL') ZRISNOW(:) = MIN(X_RI_MAX, ZRISNOW(:))
-!
-! Compute drag coefficient for momentum
-!
-CALL SURFACE_CD(ZRISNOW, PZREF, PUREF, PZ0EFF, PZ0H, ZCDSNOW, ZCDN)
-!
-! Compute friction velocity
-!
-ZUSTAR(:) = SQRT(ZCDSNOW(:)*PVMOD(:)**2)
-!
-!
-! Initialize snow layering
-!
-
-
-INLVLS_USE(:)=0
-DO JJ=1,INLVLS
- WHERE(LAYER_USE(:,JJ)) INLVLS_USE(:)=JJ
-ENDDO
-INLVLS_USE_TMP(:)=INLVLS_USE(:)
-LAYER_USE_TMP(:,:)=LAYER_USE(:,:)
-
-! ===============================================================
-!
-! 2. Compute net flux of blowing snow particles
-! and update of snowpack accordingly
-!
-! ===============================================================
-! Loop on snow covered points
- DO JWRK=1,SIZE(PSNOWSWE,1)
- T_EROSION(JWRK)=PTSTEP
- DO JJ=1,INLVLS_USE(JWRK)
- ZSNOW(JWRK) = SUM(ZSNOWDZ(JWRK,1:INLVLS_USE_TMP(JWRK)))
-
-! Turbulent flux
- IF(LSNOW_WIND) THEN
- ! Buoyancy effects in presence of suspended particles are taken into account.
- ! Iterations are necessary to reach a balance between u* and
- ! concentration in the saltation layer.
- CALL SNOWFLUX_GET_TURB(ZSNOWGRAN1(JWRK,1),ZSNOWGRAN2(JWRK,1), &
- ZSNOWLIQ(JWRK,1),ZSNOWHIST(JWRK,1),PRHOA(JWRK),ZTURBFLUX(JWRK,:), &
- PBLOWSNW_CONC(JWRK,:),PVMOD(JWRK),ZCDSNOW(JWRK),PZ0EFF(JWRK),ZMOB, &
- ZRISNOW(JWRK),PZREF(JWRK),PZ0H(JWRK),PUREF(JWRK), HSNOWRES, &
- PDIRCOSZW(JWRK),ZSALTFLUX(JWRK))
- ELSE
- ! Buoyancy effects in presence of suspended particles are not
- ! taken into account. Only thermal effect are included. Snow turbulent flux is directly computed.
- CALL SNOWFLUX_GET(ZSNOWGRAN1(JWRK,1),ZSNOWGRAN2(JWRK,1), &
- ZSNOWLIQ(JWRK,1),ZSNOWHIST(JWRK,1),ZUSTAR(JWRK),PRHOA(JWRK),ZTURBFLUX(JWRK,:), &
- PBLOWSNW_CONC(JWRK,:),PVMOD(JWRK),ZCDSNOW(JWRK),PZ0EFF(JWRK),ZMOB,ZSALTFLUX(JWRK))
- END IF
- ! Negative net flux : snow deposition
- IF(ZTURBFLUX(JWRK,2)+ZSEDFLUX(JWRK,2)+ZSALT_CONTR(JWRK)<=0.) THEN
- SWE_DEP=-(ZTURBFLUX(JWRK,2)+ZSEDFLUX(JWRK,2)+ZSALT_CONTR(JWRK))*T_EROSION(JWRK)
- PBLOWSNW_DEP(JWRK) = SWE_DEP / PTSTEP
- ZTOTFLUX_TURB(JWRK,:)=ZTOTFLUX_TURB(JWRK,:)+T_EROSION(JWRK)/PTSTEP*ZTURBFLUX(JWRK,:)
- ZTOTFLUX_SALT(JWRK)=ZTOTFLUX_SALT(JWRK)+T_EROSION(JWRK)/PTSTEP*ZSALTFLUX(JWRK)
- ZEMIMOB(JWRK) = ZEMIMOB(JWRK)+T_EROSION(JWRK)/PTSTEP*ZMOB
- T_EROSION(JWRK)=0.
- EXIT
-
- ! Positive net flux : snow erosion
- ELSE IF(ZTURBFLUX(JWRK,2)+ZSEDFLUX(JWRK,2)+ZSALT_CONTR(JWRK)>0.) THEN
- ! SWE potentially eroded during the T_EROSION
- SWE_MAX=(ZTURBFLUX(JWRK,2)+ZSEDFLUX(JWRK,2)+ZSALT_CONTR(JWRK))*T_EROSION(JWRK)
-
- ! Only layer JJ is partially eroded
- IF(SWE_MAX<=ZSNOWSWE(JWRK,1)) THEN ! Only layer JJ is partially eroded
- CALL SNOWEROSION(SWE_MAX,ZSNOWHEAT(JWRK,:),ZSNOWDZ(JWRK,:),ZSNOWSWE(JWRK,:), &
- ZSNOWRHO(JWRK,:),ZSNOWGRAN1(JWRK,:),ZSNOWGRAN2(JWRK,:) &
- ,ZSNOWHIST(JWRK,:),ZSNOWAGE(JWRK,:),ZSNOWLIQ(JWRK,:), &
- INLVLS_USE_TMP(JWRK),LAYER_USE_TMP(JWRK,:),PVMOD(JWRK),PTSTEP)
- ZTOTFLUX_TURB(JWRK,:)=ZTOTFLUX_TURB(JWRK,:)+T_EROSION(JWRK)/PTSTEP*ZTURBFLUX(JWRK,:)
- ZTOTFLUX_SALT(JWRK)=ZTOTFLUX_SALT(JWRK)+T_EROSION(JWRK)/PTSTEP*ZSALTFLUX(JWRK)
- ZEMIMOB(JWRK) = ZEMIMOB(JWRK)+T_EROSION(JWRK)/PTSTEP*ZMOB
- T_EROSION(JWRK)=0.
- EXIT
-
- ! Layer JJ is totally eroded and removed from the snowpack, next layer is then checked
- ! to see if blowing snow occurs
- ELSE
- T_EROSION(JWRK)=T_EROSION(JWRK)*(1-ZSNOWSWE(JWRK,1)/SWE_MAX)
- ZTOTFLUX_TURB(JWRK,:)=ZTOTFLUX_TURB(JWRK,:)+(1-T_EROSION(JWRK)/PTSTEP)*ZTURBFLUX(JWRK,:)
- ZTOTFLUX_SALT(JWRK)=ZTOTFLUX_SALT(JWRK)+(1-T_EROSION(JWRK)/PTSTEP)*ZSALTFLUX(JWRK)
- CALL SNOWEROSION(ZSNOWSWE(JWRK,1),ZSNOWHEAT(JWRK,:),ZSNOWDZ(JWRK,:), &
- ZSNOWSWE(JWRK,:), ZSNOWRHO(JWRK,:),ZSNOWGRAN1(JWRK,:), &
- ZSNOWGRAN2(JWRK,:),ZSNOWHIST(JWRK,:),ZSNOWAGE(JWRK,:),ZSNOWLIQ(JWRK,:), &
- INLVLS_USE_TMP(JWRK),LAYER_USE_TMP(JWRK,:),PVMOD(JWRK),PTSTEP)
- ZEMIMOB(JWRK) = ZEMIMOB(JWRK)+(1-T_EROSION(JWRK)/PTSTEP)*ZMOB
- END IF
- END IF
- END DO
-
-
- IF(T_EROSION(JWRK)>0.AND.(ZSEDFLUX(JWRK,2)+ZSALT_CONTR(JWRK))<-XUEPSI) THEN
- ! Take into account two cases :
- ! - Snow pack has been totally eroded during less
- ! than one time step and potential deposition is
- ! computed for the remaining part of the time step
- ! - Deposition of blown snow particles on bare ground
- !write(*,*) 'Deposition en excedent'
-
- SWE_DEP=-(ZSEDFLUX(JWRK,2)+ZSALT_CONTR(JWRK))*T_EROSION(JWRK)
- PBLOWSNW_DEP(JWRK) = SWE_DEP / PTSTEP
- write(*,*)'pb T_EROSION',T_EROSION(JWRK),JWRK
- write(*,*)'pb INLVLS_OLD',SWE_DEP
- write(*,*) ZTURBFLUX(JWRK,2)
-! STOP 'deposition en excedent'
- ENDIF
-
- ! Compute final snow depth that may have changed because of erosion.
- IF(INLVLS_USE_TMP(JWRK)>0) THEN
- ZSNOW(JWRK) = SUM(ZSNOWDZ(JWRK,1:INLVLS_USE_TMP(JWRK)))
- ELSE
- ZSNOW(JWRK) = 0.
- ENDIF
-
- END DO
-
-
-! ===============================================================
-!
-! 3. Final computation
-!
-! ===============================================================
-
-! Update final snowpack
- PSNOWGRAN1(:,:)=ZSNOWGRAN1(:,:)
- PSNOWGRAN2(:,:)=ZSNOWGRAN2(:,:)
- PSNOWHIST (:,:)=ZSNOWHIST (:,:)
- PSNOWAGE (:,:)=ZSNOWAGE (:,:)
- PSNOWHEAT (:,:)=ZSNOWHEAT (:,:)
- PSNOWRHO (:,:)=ZSNOWRHO (:,:)
- PSNOWSWE (:,:)=ZSNOWDZ(:,:)*PSNOWRHO(:,:)
-
-! Snow fluxes sent to Canopy or directly to Meso-NH if Canopy is activated
-
- PBLOWSNW_FLUX(:,1)=ZTOTFLUX_TURB(:,1) ! Save total number turbulent flux
- PBLOWSNW_FLUX(:,2)=ZTOTFLUX_TURB(:,2) ! Save total mass turbulent flux
- PBLOWSNW_FLUX(:,3)=ZTOTFLUX_SALT(:) ! Save streamwise saltation flux
-
-IF (LHOOK) CALL DR_HOOK('SNOWPACK_EVOL',1,ZHOOK_HANDLE)
-
-CONTAINS
-
- SUBROUTINE SNOWFLUX_GET(PSNOWGRAN1,PSNOWGRAN2, &
- PSNOWLIQ,PSNOWHIST,PUSTAR,PRHOA,PBLOWSNW_FLUX, &
- PBLOWSNW_CONC,PVMOD,PCDSNOW,PZ0EFF,PMOB,PQSALT)
-!
-!! PURPOSE
-!! -------
-! Calculates turbulent blowing snow flux and streamwise saltation flux as a function
-! of wind speed and snowpack properties
-!
-! Returns :
-! - emitted turbulent flux
-! - mobility index of emitted snow
-! - streamwise saltation flux
-!
-!
- USE MODE_BLOWSNW_MBLUTL
-
- IMPLICIT NONE
-!
-! 0.1 declarations of arguments
-!
-!
- REAL, INTENT(IN) :: PSNOWGRAN1,PSNOWGRAN2,PSNOWLIQ,PSNOWHIST
-!
- REAL, INTENT(IN) :: PUSTAR,PRHOA
- REAL, INTENT(IN) :: PVMOD,PCDSNOW,PZ0EFF
- REAL, DIMENSION(:) , INTENT(IN) :: PBLOWSNW_CONC
-
- REAL, INTENT(OUT) :: PMOB
- REAL, DIMENSION(:) , INTENT(OUT) :: PBLOWSNW_FLUX
- REAL, INTENT(OUT) :: PQSALT ! Streamwise saltation mass flux (kg/m/s)
-!
-! 0.2 declaration of local variables
-!
- REAL :: ZMBLINDEX
- REAL :: ZWIND_RFR_THR,ZWIND_FRC,ZWIND_FRC_THR, &
- ZWIND_FRC_SALT
-
- REAL :: ZCONC_SALT
-
- REAL(KIND=JPRB) :: ZHOOK_HANDLE
-!
-! 0.3 Initialization
-!
-IF (LHOOK) CALL DR_HOOK('SNOWPACK_EVOL:SNOWFLUX_GET',0,ZHOOK_HANDLE)
-
-!
- ZWIND_FRC = PUSTAR
- PMOB = SNOW_MBL_INDEX(PSNOWGRAN1, PSNOWGRAN2)
-!
-! 1. Threshold Wind at 5m
-!
- CALL WIND_RFR_THR(ZWIND_RFR_THR,PSNOWGRAN1,PSNOWGRAN2)
-!
-! 2. Threshold Friction velocity
-!
- CALL WIND_FRC_THR(ZWIND_RFR_THR,PZ0EFF,ZWIND_FRC_THR)
-!
-! 3. Blown snow particles concentration (kg/m3) in the saltation layer
-
-!
- CALL CONC_SALT(ZCONC_SALT,ZWIND_FRC,ZWIND_FRC_THR,PRHOA,PQSALT, &
- PSNOWLIQ,PSNOWHIST)
-!
-! 4. Turbulent flux of blowing snow particles
-!
- CALL SNOW_FLUX(ZCONC_SALT,PVMOD,PBLOWSNW_CONC,PCDSNOW,PBLOWSNW_FLUX)
-!
-! Snowdrfit does not occur when:
-! - surface layer is wet
-! - historical variable is higher than 1 : crust of non-transportable snow
- IF(PSNOWLIQ>0 .OR. PSNOWHIST>0) THEN
- PBLOWSNW_FLUX(:)=0
- ENDIF
-!
-
-IF (LHOOK) CALL DR_HOOK('SNOWPACK_EVOL:SNOWFLUX_GET',1,ZHOOK_HANDLE)
-
- END SUBROUTINE SNOWFLUX_GET
-
-
- SUBROUTINE SNOWFLUX_GET_TURB(PSNOWGRAN1,PSNOWGRAN2, &
- PSNOWLIQ,PSNOWHIST,PRHOA,PBLOWSNW_FLUX, &
- PBLOWSNW_CONC,PVMOD,PCDSNOW,PZ0EFF,PMOB,PRI,PZREF,PZ0H, &
- PUREF, HSNOWRES, PDIRCOSZW,PQSALT )
-
- USE MODE_BLOWSNW_MBLUTL
- USE MODD_BLOWSNW_SURF
-
-IMPLICIT NONE
-!
-!
-!! PURPOSE
-!! -------
-! Routine that compute emitted turbulent flux accounting for influence of
-! blown snow particles on the wind profile.
-! Return :
-! - emitted turbulent flux
-! - mobility index of emitted snow
-! - streamwise saltation flux
-! - Richardson number including particle-induced buoyancy
-! - modified drag coefficient for momentum
-!
-! Routine to be tested!!
-!
-! 0.1 declarations of arguments
-!
-!
- REAL, INTENT(IN) :: PSNOWGRAN1,PSNOWGRAN2,PSNOWLIQ,PSNOWHIST
-!
- REAL, INTENT(IN) :: PRHOA
- REAL, INTENT(IN) :: PVMOD,PZ0EFF
- REAL, INTENT(IN) :: PZREF, PZ0H, PUREF, PDIRCOSZW
- REAL, DIMENSION(:) , INTENT(IN) :: PBLOWSNW_CONC
-
- CHARACTER(LEN=*), INTENT(IN) :: HSNOWRES
-!
- REAL, INTENT(INOUT) :: PCDSNOW, PRI
-!
- REAL, INTENT(OUT) :: PMOB
- REAL, DIMENSION(:) , INTENT(OUT) :: PBLOWSNW_FLUX
- REAL, INTENT(OUT) :: PQSALT ! Streamwise saltation mass flux (kg/m/s)
-!
-! 0.2 declaration of local variables
-!
- REAL :: ZMBLINDEX
- REAL :: ZWIND_RFR_THR,ZWIND_FRC,ZWIND_FRC_THR
- REAL :: ZCD_TEMP, ZCDN
- REAL :: ZCONC_SALT
- REAL :: ZRI_TOT ! Total richardson number (thermal + particle
-! buoyancy included)
- REAL :: ZRI_PART ! Particle richardson number ( particle
-! buoyancy included)
- REAL :: ZB,ZDELTA,ZY,ZY_DRIV,ZY_TEMP
- INTEGER :: II,JJ
- INTEGER :: NITER
-
- REAL(KIND=JPRB) :: ZHOOK_HANDLE
-
-
-!
-! 0.3 Initialization
-!
-IF (LHOOK) CALL DR_HOOK('SNOWPACK_EVOL:SNOWFLUX_GET_TURB',0,ZHOOK_HANDLE)
-!
- ZWIND_FRC = SQRT(PCDSNOW*PVMOD**2)
- ZCD_TEMP = PCDSNOW ! Drag coefficient without particle-induced buoyancy
-! and increase in aerodynamic rougnhess due to saltation
- ZRI_TOT = PRI ! Richardson number without particle-induced buoyancy
-! Parameters for Newton's method
- ZDELTA = 0.01
- JJ=0
- ZB = ZWIND_FRC-999
- NITER = 5 ! Number of iterations to compute u* and Conc_salt
- !
- ! Compute mobility index based on surface layer properties
- PMOB = SNOW_MBL_INDEX(PSNOWGRAN1, PSNOWGRAN2);
-!
-! 1. Threshold Wind at 5m
-!
- CALL WIND_RFR_THR(ZWIND_RFR_THR,PSNOWGRAN1,PSNOWGRAN2)
-!
-! 2. Threshold Friction velocity
-!
- CALL WIND_FRC_THR(ZWIND_RFR_THR,PZ0EFF,ZWIND_FRC_THR)
-!
-! 3. Newton method to find u*, Ri and CD
-!
-DO WHILE(ABS(ZWIND_FRC-ZB)>0.001)
- JJ=JJ+1
- ZB=ZWIND_FRC
- IF(JJ>NITER) EXIT
- CALL SOLVE_CD(ZWIND_FRC+ZDELTA,PVMOD,PRHOA,PBLOWSNW_CONC(2),PZREF,PUREF,PDIRCOSZW, &
- ZWIND_FRC_THR,PRI,PZ0EFF,HSNOWRES,PZ0H,ZY_TEMP,ZRI_TOT,ZCD_TEMP, &
- PSNOWLIQ,PSNOWHIST)
-
- CALL SOLVE_CD(ZWIND_FRC,PVMOD,PRHOA,PBLOWSNW_CONC(2),PZREF,PUREF,PDIRCOSZW, &
- ZWIND_FRC_THR,PRI,PZ0EFF,HSNOWRES,PZ0H,ZY,ZRI_TOT,ZCD_TEMP, &
- PSNOWLIQ,PSNOWHIST)
-
- ZY_DRIV = (ZY_TEMP-ZY)/ZDELTA
- ZWIND_FRC = ZB- ZY/ZY_DRIV
-END DO
-!
-! 4. Blown snow particles concentration in saltation (kg/m3)
-!
- CALL CONC_SALT(ZCONC_SALT,ZWIND_FRC,ZWIND_FRC_THR,PRHOA,PQSALT, &
- PSNOWLIQ,PSNOWHIST)
-!
-! 5. Turbulent flux of blown snow particles
-!
- CALL SNOW_FLUX(ZCONC_SALT,PVMOD,PBLOWSNW_CONC,ZCD_TEMP,PBLOWSNW_FLUX)
-!
-! Snowdrfit does not occur when :
-! - surface layer is wet
-! - historical variable is higher than 1 : crust of non-transportable snow
-!
- IF(PSNOWLIQ>0 .OR. PSNOWHIST>0) THEN
- PBLOWSNW_FLUX(:)=0
- ENDIF
-
-! Update Richardson number and drag coefficient for momentum
- PRI = ZRI_TOT
- PCDSNOW = ZCD_TEMP
-
- IF (LHOOK) CALL DR_HOOK('SNOWPACK_EVOL:SNOWFLUX_GET_TURB',1,ZHOOK_HANDLE)
-
- END SUBROUTINE SNOWFLUX_GET_TURB
-
-
- SUBROUTINE SNOWEROSION(PSWE_REM,PSNOWHEAT,PSNOWDZ,PSNOWSWE, &
- PSNOWRHO,PSNOWGRAN1,PSNOWGRAN2,PSNOWHIST,&
- PSNOWAGE,PSNOWLIQ,INLVLS_USE,LAYER_USE, &
- PVMOD, PTSTEP)
-!
-USE MODE_SNOW3L
-USE MODD_SNOW_PAR
-USE MODD_CSTS, ONLY : XLMTT, XTT, XCI
-USE MODE_BLOWSNW_MBLUTL
-!
-IMPLICIT NONE
-!
-!! PURPOSE
-!! -------
-!
-! Erosion of the snowpack surface layer :
-! reduces heat content, SWE and thickness
-!
-! ASSUMPTION: the snowpack surface layer is never totally eroded.
-! Total erosion is insured earlier.
-!
-! 0.1 declarations of arguments
-!
-!
- REAL, INTENT(IN) :: PSWE_REM ! Snow mass to be removed
-! from the snow top layer (never exceed top layer thickness)
- REAL, INTENT(IN) :: PVMOD, PTSTEP
-!
-! Top layers properties to be updated
-!
- REAL, DIMENSION(:), INTENT(INOUT) :: PSNOWHEAT,PSNOWDZ,PSNOWSWE,PSNOWRHO, &
- PSNOWGRAN1,PSNOWGRAN2,PSNOWHIST,PSNOWAGE,PSNOWLIQ
-
- LOGICAL, DIMENSION(:), INTENT(INOUT) :: LAYER_USE
- INTEGER,INTENT(INOUT) :: INLVLS_USE
-
-! 0.2 declaration of local variables
-!
- REAL :: ZSNOW_REM
- REAL ::ZSNOWHMASS,ZSNOWTEMP
- REAL :: ZRMOB,ZRDRIFT,ZPROFEQ,ZRT,ZSNOWRHO2
- REAL :: ZDRO, ZDGR1,ZDGR2
- REAL, DIMENSION(SIZE(PSNOWRHO,1)) :: ZSCAP
- INTEGER :: JLAYER,KNLVLS
- LOGICAL :: OEVOL_GRAIN
- REAL(KIND=JPRB) :: ZHOOK_HANDLE
-
-!
-! 0.3 initialization
-!
-IF (LHOOK) CALL DR_HOOK('SNOWPACK_EVOL:SNOWEROSION',0,ZHOOK_HANDLE)
-!
-ZSNOW_REM=0.0
-ZPROFEQ=0.0
-KNLVLS = SIZE(PSNOWSWE,1)
-!
-! Option to simulated evolution of snow properties in areas where
-! erosion occurs
-!
-OEVOL_GRAIN = .FALSE.
-
-
-! 1. Remove snow from surface layer
-
- IF(PSWE_REM==PSNOWSWE(1)) THEN ! Surface layer is totally removed.
- DO JLAYER=1,INLVLS_USE-1
- PSNOWSWE(JLAYER)=PSNOWSWE(JLAYER+1)
- PSNOWHEAT(JLAYER)=PSNOWHEAT(JLAYER+1)
- PSNOWDZ(JLAYER)=PSNOWDZ(JLAYER+1)
- PSNOWRHO(JLAYER)=PSNOWRHO(JLAYER+1)
- PSNOWLIQ(JLAYER)=PSNOWLIQ(JLAYER+1)
- PSNOWGRAN1(JLAYER)=PSNOWGRAN1(JLAYER+1)
- PSNOWGRAN2(JLAYER)=PSNOWGRAN2(JLAYER+1)
- PSNOWHIST(JLAYER)=PSNOWHIST(JLAYER+1)
- PSNOWAGE(JLAYER)=PSNOWAGE(JLAYER+1)
- ENDDO
- PSNOWSWE(INLVLS_USE)=0.0
- PSNOWRHO(INLVLS_USE)=999.
- PSNOWDZ(INLVLS_USE)=0.
- PSNOWGRAN1(INLVLS_USE)=0.
- PSNOWGRAN2(INLVLS_USE)=0.
- PSNOWHIST(INLVLS_USE)=0.
- PSNOWAGE(INLVLS_USE)=0.
- PSNOWHEAT(INLVLS_USE)=0.
- PSNOWLIQ(INLVLS_USE)=0.
- LAYER_USE(INLVLS_USE)=.FALSE.
- INLVLS_USE=INLVLS_USE-1
- ELSE
- ZSCAP = SNOW3LSCAP(PSNOWRHO)
- ZSNOWTEMP = XTT + (PSNOWHEAT(1) + &
- XLMTT*PSNOWRHO(1)*PSNOWDZ(1))/ &
- (ZSCAP(1)*MAX(XSNOWDMIN/KNLVLS,PSNOWDZ(1)))
- ZSNOWTEMP = MIN(XTT, ZSNOWTEMP)
- ZSNOWHMASS = PSWE_REM*(XCI*(ZSNOWTEMP-XTT)-XLMTT)
-!
-! 2.Reduce total pack depth
- ZSNOW_REM=PSWE_REM/PSNOWRHO(1)
- PSNOWDZ(1)=PSNOWDZ(1)-ZSNOW_REM
-! 3.Reduce heat content
- PSNOWHEAT(1)=PSNOWHEAT(1)-ZSNOWHMASS
- PSNOWSWE(1)=PSNOWDZ(1)*PSNOWRHO(1)
-! 4.Evolution of the remaining layer to simulate subgrid density increase
-! and grain shape and size evolution due snowdrift.
- IF(OEVOL_GRAIN) THEN
- DO JLAYER=1,INLVLS_USE
- ZSNOWRHO2 = PSNOWRHO(JLAYER)
- ZRMOB = SNOW_MBL_INDEX(PSNOWGRAN1(JLAYER),PSNOWGRAN2(JLAYER))
- IF(PSNOWHIST(JLAYER).ge.2) ZRMOB = MIN(ZRMOB,-0.0583)
- ZRDRIFT = ZRMOB-(2.868*EXP(-0.085*PVMOD)-1.)
- ZPROFEQ = ZPROFEQ + 0.5* PSNOWDZ(JLAYER) * (3.25-ZRDRIFT)
- ZRT=MAX(0.,ZRDRIFT*EXP(-ZPROFEQ/0.1))
-
- IF(PSNOWRHO(JLAYER).lt.350.) THEN
- ZDRO=(350.-PSNOWRHO(JLAYER))*PTSTEP/12./3600.*ZRT
- PSNOWRHO(JLAYER)=MIN(350.,PSNOWRHO(JLAYER)+ZDRO)
- PSNOWDZ(JLAYER)=PSNOWDZ(JLAYER)*ZSNOWRHO2/ &
- PSNOWRHO(JLAYER)
- END IF
-! Dendritic snow
- IF (PSNOWGRAN1(JLAYER).lt.0.) THEN
- ZDGR1=-PSNOWGRAN1(JLAYER)*PTSTEP/12./3600.*ZRT*0.5
- ZDGR1=MIN(ZDGR1,-0.99*PSNOWGRAN1(JLAYER))
- ZDGR2=ZRT*PTSTEP/12./3600.*(99.-PSNOWGRAN2(JLAYER))
- PSNOWGRAN1(JLAYER)=PSNOWGRAN1(JLAYER)+ZDGR1
- PSNOWGRAN2(JLAYER)=MIN(99.,PSNOWGRAN2(JLAYER)+zdgr2)
- ELSE ! Non dendritic snow
- ZDGR1=PTSTEP/12./3600.*ZRT*(99.-PSNOWGRAN1(JLAYER))
- ZDGR2=PTSTEP/12./3600.*ZRT*5./10000.
- PSNOWGRAN1(JLAYER)=MIN(99.,PSNOWGRAN1(JLAYER)+ZDGR1)
- PSNOWGRAN2(JLAYER)=MAX(0.0003,PSNOWGRAN2(JLAYER)-ZDGR2)
- END IF
- END DO
- END IF
-
- END IF
-
- IF (LHOOK) CALL DR_HOOK('SNOWPACK_EVOL:SNOWEROSION',1,ZHOOK_HANDLE)
-
- END SUBROUTINE SNOWEROSION
-
-END SUBROUTINE SNOWPACK_EVOL
--
GitLab