source: LMDZ6/trunk/libf/phylmd/inlandsis/inlandsis.F @ 4744

Last change on this file since 4744 was 3900, checked in by evignon, 3 years ago

Commission de la nouvelle interface LMDZ-SISVAT
la prochaine commission consistera a supprimer l'ancien repertoire sisvat
et a faire un peu de nettoyage.

File size: 64.0 KB
Line 
1      subroutine INLANDSIS(SnoMod,BloMod,jjtime,debut)
2
3      USE dimphy
4
5!--------------------------------------------------------------------------+
6!     INLANDSIS module                                                     |
7!     Simplified SISVAT module, containing ice and snow processes for      |
8!     ice-covered surfaces                                                 |
9!     version MARv3, november 2020                                         |
10!     SubRoutine INLANDSIS contains the fortran 77 code of the             |
11!                Soil/Ice Snow Vegetation Atmosphere Transfer Scheme       |
12!                                                                          |
13!--------------------------------------------------------------------------+
14!     PARAMETERS:  klonv: Total Number of columns =                        |
15!     ^^^^^^^^^^        = Total Number of continental     grid boxes       |
16!                       X       Number of Mosaic Cell per grid box         |
17!                                                                          |
18!     INPUT:   daHost   : Date Host Model                                  |
19!     ^^^^^                                                                |
20!                                                                          |
21!     INPUT:   LSmask   : 1:          Land       MASK                      |
22!     ^^^^^               0:          Sea        MASK                      |
23!              isotSV   = 0,...,12:   Soil       Type                      |
24!                         0:          Water,          Liquid (Sea, Lake)   |
25!                        12:          Water, Solid           (Ice)         |
26!                                                                          |
27!     INPUT:   coszSV   : Cosine of the Sun Zenithal Distance          [-] |
28!     ^^^^^    sol_SV   : Surface Downward  Solar      Radiation    [W/m2] |
29!              IRd_SV   : Surface Downward  Longwave   Radiation    [W/m2] |
30!              drr_SV   : Rain  Intensity                        [kg/m2/s] |
31!              dsn_SV   : Snow  Intensity                      [mm w.e./s] |
32!              dsnbSV   : Snow  Intensity,  Drift Fraction             [-] |
33!              dbs_SV   : Drift Amount                           [mm w.e.] |
34!              za__SV   : Surface Boundary Layer (SBL) Height          [m] |
35!              VV__SV   :(SBL Top)   Wind Velocity                   [m/s] |
36!              TaT_SV   : SBL Top    Temperature                       [K] |
37!              rhT_SV   : SBL Top    Air  Density                  [kg/m3] |
38!              QaT_SV   : SBL Top    Specific  Humidity            [kg/kg] |
39!              qsnoSV   : SBL Mean   Snow      Content             [kg/kg] |
40!              alb0SV   : Soil Basic Albedo                            [-] |
41!              slopSV   : Surface    Slope                             [-] |
42!              dt__SV   : Time  Step                                   [s] |
43!                                                                          |
44!     INPUT /  isnoSV   = total Nb of Ice/Snow Layers                      |
45!     OUTPUT:  ispiSV   = 0,...,nsno: Uppermost Superimposed Ice Layer     |
46!     ^^^^^^   iiceSV   = total Nb of Ice      Layers                      |
47!              istoSV   = 0,...,5 :   Snow     History (see istdSV data)   |
48!                                                                          |
49!     INPUT /  alb_SV   : Surface        Albedo                        [-] |
50!     OUTPUT:  emi_SV   : Surface        Emissivity                    [-] |
51!     ^^^^^^   IRs_SV   : Soil           IR Flux  (negative)        [W/m2] |
52!              LMO_SV   : Monin-Obukhov               Scale            [m] |
53!              us__SV   : Friction          Velocity                 [m/s] |
54!              uts_SV   : Temperature       Turbulent Scale          [m/s] |
55!              uqs_SV   : Specific Humidity Velocity                 [m/s] |
56!              uss_SV   : Blowing Snow      Turbulent Scale          [m/s] |
57!              usthSV   : Blowing Snow      Erosion   Threshold      [m/s] |
58!              Z0m_SV   : Momentum     Roughness Length                [m] |
59!              Z0mmSV   : Momentum     Roughness Length (time mean)    [m] |
60!              Z0mnSV   : Momentum     Roughness Length (instantaneous)[m] |
61!              Z0SaSV   : Sastrugi     Roughness Length                [m] |
62!              Z0e_SV   : Erosion Snow Roughness Length                [m] |
63!              Z0emSV   : Erosion Snow Roughness Length (time mean)    [m] |
64!              Z0enSV   : Erosion Snow Roughness Length (instantaneous)[m] |
65!              Z0roSV   : Subgrid Topo Roughness Length                [m] |
66!              Z0h_SV   : Heat         Roughness Length                [m] |
67!              TsisSV   : Soil/Ice Temperatures (layers -nsol,-nsol+1,..,0)|
68!                       & Snow     Temperatures (layers  1,2,...,nsno) [K] |
69!              ro__SV   : Soil/Snow Volumic Mass                   [kg/m3] |
70!              eta_SV   : Soil/Snow Water   Content                [m3/m3] |
71!              G1snSV   : snow dendricity/sphericity                       |
72!              G2snSV   : snow sphericity/grain size                       |
73!              dzsnSV   : Snow Layer        Thickness                  [m] |
74!              agsnSV   : Snow       Age                             [day] |
75!              BufsSV   : Snow Buffer Layer              [kg/m2] .OR. [mm] |
76!              BrosSV   : Snow Buffer Layer Density      [kg/m3]           |
77!              BG1sSV   : Snow Buffer Layer Dendricity / Sphericity    [-] |
78!              BG2sSV   : Snow Buffer Layer Sphericity / Size [-] [0.1 mm] |
79!              rusnSV   : Surficial   Water              [kg/m2] .OR. [mm] |
80!                                                                          |
81!     OUTPUT:  no__SV   : OUTPUT file Unit Number                      [-] |
82!     ^^^^^^   i___SV   : OUTPUT point   i Coordinate                  [-] |
83!              j___SV   : OUTPUT point   j Coordinate                  [-] |
84!              n___SV   : OUTPUT point   n Coordinate                  [-] |
85!              lwriSV   : OUTPUT point vec Index                       [-] |
86!                                                                          |
87!     OUTPUT:  IRu_SV   : Upward     IR Flux (+, upw., effective)      [K] |
88!     ^^^^^^   hSalSV   : Saltating Layer Height                       [m] |
89!              qSalSV   : Saltating Snow  Concentration            [kg/kg] |
90!              RnofSV   : RunOFF Intensity                       [kg/m2/s] |
91!                                                                          |
92!     Internal Variables:                                                  |
93!     ^^^^^^^^^^^^^^^^^^                                                   |
94!              NLaysv   = New            Snow Layer Switch             [-] |
95!              albisv   : Snow/Ice/Water/Soil Integrated Albedo        [-] |
96!              SoSosv   : Absorbed Solar Radiation by Surfac.(Normaliz)[-] |
97!              TBr_sv   : Brightness Temperature                       [K] |
98!              IRupsv   : Upward     IR Flux (-, upw.)              [W/m2] |
99!              ram_sv   : Aerodynamic Resistance for Momentum        [s/m] |
100!              rah_sv   : Aerodynamic Resistance for Heat            [s/m] |
101!              Evp_sv   : Evaporation                              [kg/m2] |
102!              EvT_sv   : Evapotranspiration                       [kg/m2] |
103!              HSs_sv   : Surface    Sensible Heat Flux + => absorb.[W/m2] |
104!              HLs_sv   : Surface    Latent   Heat Flux + => absorb.[W/m2] |
105!              Lx_H2O   : Latent Heat of Vaporization/Sublimation   [J/kg] |
106!              Tsrfsv   : Surface    Temperature                       [K] |
107!              sEX_sv   : Verticaly Integrated Extinction Coefficient  [-] |
108!              LSdzsv   : Vertical   Discretization Factor             [-] |
109!                       =    1. Soil                                       |
110!                       = 1000. Ocean                                      |
111!              z_snsv   : Snow Pack  Thickness                         [m] |
112!              zzsnsv   : Snow Pack  Thickness                         [m] |
113!              albssv   : Soil       Albedo                            [-] |
114!              Eso_sv   : Soil+Snow       Emissivity                   [-] |
115!              Khydsv   : Soil   Hydraulic    Conductivity           [m/s] |
116!                                                                          |
117!              ETSo_0   : Snow/Soil Energy Power, before Forcing    [W/m2] |
118!              ETSo_1   : Snow/Soil Energy Power, after  Forcing    [W/m2] |
119!              ETSo_d   : Snow/Soil Energy Power         Forcing    [W/m2] |
120!              EqSn_0   : Snow      Energy, before Phase Change     [J/m2] |
121!              EqSn_1   : Snow      Energy, after  Phase Change     [J/m2] |
122!              EqSn_d   : Snow      Energy,       net    Forcing    [J/m2] |
123!              Enrsvd   : SVAT      Energy Power         Forcing    [W/m2] |
124!              Enrbal   : SVAT      Energy Balance                  [W/m2] |
125!              Wats_0   : Soil Water,  before Forcing                 [mm] |
126!              Wats_1   : Soil Water,  after  Forcing                 [mm] |
127!              Wats_d   : Soil Water          Forcing                 [mm] |
128!              SIWm_0   : Snow        initial Mass               [mm w.e.] |
129!              SIWm_1   : Snow        final   Mass               [mm w.e.] |
130!              SIWa_i   : Snow Atmos. initial Forcing            [mm w.e.] |
131!              SIWa_f   : Snow Atmos. final   Forcing(noConsumed)[mm w.e.] |
132!              SIWe_i   : SnowErosion initial Forcing            [mm w.e.] |
133!              SIWe_f   : SnowErosion final   Forcing(noConsumed)[mm w.e.] |
134!              SIsubl   : Snow sublimed/deposed  Mass            [mm w.e.] |
135!              SImelt   : Snow Melted            Mass            [mm w.e.] |
136!              SIrnof   : Surficial Water + Run OFF Change       [mm w.e.] |
137!              SIvAcr   : Sea-Ice    vertical Acretion           [mm w.e.] |
138!              Watsvd   : SVAT Water          Forcing                 [mm] |
139!              Watbal   : SVAT Water  Balance                       [W/m2] |
140!                                                                          |
141!              vk2      : Square of Von Karman Constant                [-] |
142!              sqrCm0   : Factor of   Neutral Drag Coeffic.Momentum  [s/m] |
143!              sqrCh0   : Factor of   Neutral Drag Coeffic.Heat      [s/m] |
144!              EmiSol   : Soil          Emissivity                     [-] |
145!              EmiSno   : Snow          Emissivity                     [-] |
146!              EmiWat   : Water         Emissivity                     [-] |
147!              Z0mLnd   :          Land Roughness Length               [m] |
148!              sqrrZ0   : u*t/u*                                           |
149!              f_eff    : Marticorena & B. 1995 JGR (20)                   |
150!              A_Fact   : Fundamental * Roughness                          |
151!              Z0mBSn   :         BSnow Roughness Length               [m] |
152!              Z0mBS0   : Mimimum BSnow Roughness Length (blown* )     [m] |
153!              Z0m_Sn   :          Snow Roughness Length (surface)     [m] |
154!              Z0m_S0   : Mimimum  Snow Roughness Length               [m] |
155!              Z0m_S1   : Maximum  Snow Roughness Length               [m] |
156!              Z0_GIM   : Minimum GIMEX Roughness Length               [m] |
157!              Z0_ICE   : Sea Ice ISW   Roughness Length               [m] |
158!                                                                          |
159!                                                                          |
160!--------------------------------------------------------------------------+
161     
162
163
164! Global Variables
165! ================
166
167
168      USE VARphy
169      USE VAR_SV
170      USE VARdSV
171      USE VAR0SV
172      USE VARxSV
173      USE VARySV
174      USE VARtSV
175      USE surface_data, ONLY: is_ok_z0h_rn,
176     .                        is_ok_density_kotlyakov,
177     .                        prescribed_z0m_snow,
178     .                        iflag_z0m_snow,
179     .                        iflag_tsurf_inlandsis,
180     .                        iflag_temp_inlandsis,
181     .                        discret_xf, buf_sph_pol,buf_siz_pol     
182
183      IMPLICIT NONE
184
185      logical   SnoMod
186      logical   BloMod
187      logical   debut
188      integer   jjtime
189
190
191! Internal Variables
192! ==================
193
194! Non Local
195! ---------
196
197      real      TBr_sv(klonv)                 ! Brightness Temperature
198      real      IRdwsv(klonv)                 ! DOWNward   IR Flux
199      real      IRupsv(klonv)                 ! UPward     IR Flux
200      real      d_Bufs,Bufs_N                 ! Buffer Snow Layer Increment
201      real      Buf_ro,Bros_N                 ! Buffer Snow Layer Density
202      real      BufPro                        ! Buffer Snow Layer Density
203      real      Buf_G1,BG1__N                 ! Buffer Snow Layer Dendr/Sphe[-]
204      real      Buf_G2,BG2__N                 ! Buffer Snow Layer Spher/Size[-]
205      real      Bdzssv(klonv)                 ! Buffer Snow Layer Thickness
206      real      z_snsv(klonv)                 ! Snow-Ice, current Thickness
207
208
209
210! Local
211! -----
212
213      integer   iwr
214      integer   ikl   ,isn   ,isl   ,ist      !
215      integer   ist__s,ist__w                 ! Soil/Water Body Identifier
216      integer   growth                        ! Seasonal               Mask
217      integer   LISmsk                        ! Land+Ice / Open    Sea Mask
218      integer   LSnMsk                        ! Snow-Ice / No Snow-Ice Mask
219      integer   IceMsk,IcIndx(klonv)          !      Ice / No      Ice Mask
220      integer   SnoMsk                        ! Snow     / No Snow     Mask
221      real      roSMin,roSMax,roSn_1,roSn_2,roSn_3   ! Fallen Snow Density (PAHAUT)
222      real      Dendr1,Dendr2,Dendr3          ! Fallen Snow Dendric.(GIRAUD)
223      real      Spher1,Spher2,Spher3,Spher4   ! Fallen Snow Spheric.(GIRAUD)
224      real      Polair                        ! Polar  Snow Switch
225      real      PorSno,Salt_f,PorRef   !
226c #sw real      PorVol,rWater                 !
227c #sw real      rusNEW,rdzNEW,etaNEW          !
228      real      ro_new                        !
229      real      TaPole                        ! Maximum     Polar Temperature
230      real      T__Min                        ! Minimum realistic Temperature
231      real      EmiSol                        ! Emissivity of       Soil
232      real      EmiSno                        ! Emissivity of            Snow
233      real      EmiWat                        ! Emissivity of a Water Area
234      real      vk2                           ! Square of Von Karman Constant
235      real      u2star                        !(u*)**2
236      real      Z0mLnd                        !          Land Roughness Length
237c #ZN real      sqrrZ0                        ! u*t/u*
238      real      f_eff                         ! Marticorena & B. 1995 JGR (20)
239      real      A_Fact                        ! Fundamental * Roughness
240      real      Z0m_nu                        ! Smooth R Snow Roughness Length
241      real      Z0mBSn                        !         BSnow Roughness Length
242      real      Z0mBS0                        ! Mimimum BSnow Roughness Length
243      real      Z0m_S0                        ! Mimimum  Snow Roughness Length
244      real      Z0m_S1                        ! Maximum  Snow Roughness Length
245c #SZ real      Z0Sa_N                        ! Regime   Snow Roughness Length
246c #SZ real      Z0SaSi                        ! 1.IF Rgm Snow Roughness Length
247c #GL real      Z0_GIM                        ! Mimimum GIMEX Roughness Length
248      real      Z0_ICE                        ! Ice ISW   Roughness Length
249      real      Z0m_Sn,Z0m_90                 ! Snow  Surface Roughness Length
250      real      SnoWat                        ! Snow Layer    Switch
251      real      rstar,alors                   !
252      real      rstar0,rstar1,rstar2          !
253      real      SameOK                        ! 1. => Same Type of Grains
254      real      G1same                        ! Averaged G1,  same Grains
255      real      G2same                        ! Averaged G2,  same Grains
256      real      typ__1                        ! 1. => Lay1 Type: Dendritic
257      real      zroNEW                        ! dz X ro, if fresh Snow
258      real      G1_NEW                        ! G1,      if fresh Snow
259      real      G2_NEW                        ! G2,      if fresh Snow
260      real      zroOLD                        ! dz X ro, if old   Snow
261      real      G1_OLD                        ! G1,      if old   Snow
262      real      G2_OLD                        ! G2,      if old   Snow
263      real      SizNEW                        ! Size,    if fresh Snow
264      real      SphNEW                        ! Spheric.,if fresh Snow
265      real      SizOLD                        ! Size,    if old   Snow
266      real      SphOLD                        ! Spheric.,if old   Snow
267      real      Siz_av                        ! Averaged    Grain Size
268      real      Sph_av                        ! Averaged    Grain Spher.
269      real      Den_av                        ! Averaged    Grain Dendr.
270      real      G1diff                        ! Averaged G1, diff. Grains
271      real      G2diff                        ! Averaged G2, diff. Grains
272      real      G1                            ! Averaged G1
273      real      G2                            ! Averaged G2
274      real      param                         ! Polynomial   fit z0=f(T)
275      real      Z0_obs                        ! Fit Z0_obs=f(T) (m)
276      real      tamin                         ! min T of linear fit (K)
277      real      tamax                         ! max T of linear fit (K)
278      real      coefa,coefb,coefc,coefd       ! Coefs for z0=f(T)
279      real      ta1,ta2,ta3                   ! Air temperature thresholds
280      real      z01,z02,z03                   ! z0 thresholds
281      real      tt_c,vv_c                     ! Critical param.
282      real      tt_tmp,vv_tmp,vv_virt         ! Temporary variables
283      real      e_prad,e1pRad,A_Rad0,absg_V,absgnI,exdRad ! variables for SoSosv calculations
284      real      zm1, zm2, coefslope                    ! variables for surface temperature extrapolation
285! for Aeolian erosion and blowing snow
286      integer   nit   ,iit
287      real      Fac                           ! Correc. factor for drift ratio
288      real      dusuth,signus
289      real      sss__F,sss__N
290      real      sss__K,sss__G
291      real      us_127,us_227,us_327,us_427,us_527
292      real      VVa_OK, usuth0
293      real      ssstar
294      real      SblPom
295      real      rCd10n                        ! Square root of drag coefficient
296      real      DendOK                        ! Dendricity Switch
297      real      SaltOK                        ! Saltation  Switch
298      real      MeltOK                        ! Saltation  Switch (Melting Snow)
299      real      SnowOK                        ! Pack Top   Switch
300      real      SaltM1,SaltM2,SaltMo,SaltMx   ! Saltation  Parameters
301      real      ShearX, ShearS                ! Arg. Max Shear Stress
302      real      Por_BS                        ! Snow Porosity
303      real      Salt_us                       ! New thresh.friction velocity u*t
304      real      Fac_Mo,ArguSi,FacRho          ! Numerical factors for u*t
305      real      SaltSI(klonv,0:nsno)          ! Snow Drift Index              !
306      real      MIN_Mo                        ! Minimum Mobility Fresh Fallen *
307      character*3    qsalt_param              ! Switch for saltation flux param.
308      character*3    usth_param               ! Switch for u*t param
309
310
311! Internal DATA
312! =============
313
314      data      T__Min / 200.00/              ! Minimum realistic Temperature
315      data      TaPole / 268.15/              ! Maximum Polar Temperature (value from C. Agosta)
316      data      roSMin / 300.  /              ! Minimum Snow  Density
317      data      roSMax / 400.  /              ! Max Fresh Snow Density
318      data      tt_c   / -2.0  /              ! Critical Temp. (degC)
319      data      vv_c   / 14.3  /              ! Critical Wind speed (m/s)
320      data      roSn_1 / 109.  /              ! Fall.Sno.Density, Indep. Param.
321      data      roSn_2 /   6.  /              ! Fall.Sno.Density, Temper.Param.
322      data      roSn_3 /  26.  /              ! Fall.Sno.Density, Wind   Param.
323      data      Dendr1 /  17.12/              ! Fall.Sno.Dendric.,Wind 1/Param.
324      data      Dendr2 / 128.  /              ! Fall.Sno.Dendric.,Wind 2/Param.
325      data      Dendr3 / -20.  /              ! Fall.Sno.Dendric.,Indep. Param.
326      data      Spher1 /   7.87/              ! Fall.Sno.Spheric.,Wind 1/Param.
327      data      Spher2 /  38.  /              ! Fall.Sno.Spheric.,Wind 2/Param.
328      data      Spher3 /  50.  /              ! Fall.Sno.Spheric.,Wind 3/Param.
329      data      Spher4 /  90.  /              ! Fall.Sno.Spheric.,Indep. Param.
330      data      EmiSol /   0.99999999/        ! 0.94Emissivity of Soil
331      data      EmiWat /   0.99999999/        ! Emissivity of a Water Area
332      data      EmiSno /   0.99999999/        ! Emissivity of Snow
333
334     
335!     DATA      Emissivities                  ! Pielke, 1984, pp. 383,409
336
337      data      Z0mBS0 /   0.5e-6/            ! MINimum Snow Roughness Length
338                                              ! for Momentum if Blowing Snow
339                                              ! Gallee et al. 2001 BLM 99 (19)
340      data      Z0m_S0/    0.00005/           ! MINimum Snow Roughness Length
341                                              ! MegaDunes    included
342      data      Z0m_S1/    0.030  /           ! MAXimum Snow Roughness Length
343                                              !        (Sastrugis)
344c #GL data      Z0_GIM/    0.0013/            ! Ice Min Z0 = 0.0013 m (Broeke)
345!                                             ! Old Ice Z0 = 0.0500 m (Bruce)
346!                                             !              0.0500 m (Smeets)
347!                                             !              0.1200 m (Broeke)
348      data      Z0_ICE/    0.0010/            ! Sea-Ice Z0 = 0.0010 m (Andreas)
349!                                             !    (Ice Station Weddel -- ISW)
350! for aerolian erosion
351      data      SblPom/ 1.27/   ! Lower Boundary Height Parameter
352C +                             ! for Suspension
353C +                             ! Pommeroy, Gray and Landine 1993,
354C +                             ! J. Hydrology, 144(8) p.169
355      data      nit   / 5   /   ! us(is0,uth) recursivity: Nb Iterations
356cc#AE data      qsalt_param/"bin"/ ! saltation part. conc. from Bintanja 2001 (p
357      data      qsalt_param/"pom"/ ! saltation part. conc. from Pomeroy and Gray
358cc#AE data      usth_param/"lis"/  ! u*t from Liston et al. 2007
359      data      usth_param/"gal"/  ! u*t from Gallee et al. 2001
360      data      SaltMx/-5.83e-2/
361
362      vk2    =  vonKrm  *  vonKrm             ! Square of Von Karman Constant
363
364
365! BEGIN.main.
366! ===========================
367
368
369
370
371! "Soil" Humidity of Water Bodies
372! ===============================
373
374      DO ikl=1,knonv
375
376          ist    =      isotSV(ikl)                       ! Soil Type
377          ist__s =  min(ist, 1)                           ! 1 => Soil
378          ist__w =  1 - ist__s                            ! 1 => Water Body
379        DO isl=-nsol,0
380          eta_SV(ikl,isl) = eta_SV(ikl,isl) * ist__s      ! Soil
381     .                    + etadSV(ist)     * ist__w      ! Water Body
382        END DO
383
384
385! Vertical Discretization Factor
386! ==============================
387
388          LSdzsv(ikl)     =                   ist__s      ! Soil
389     .                    + OcndSV          * ist__w      ! Water Body
390      END DO
391
392
393
394
395
396      IF (SnoMod)                            THEN
397
398 
399C +--Aeolian erosion and Blowing Snow
400C +==================================
401
402
403
404        DO ikl=1,knonv
405            usthSV(ikl) =                     1.0e+2
406        END DO
407
408
409        IF (BloMod) THEN
410 
411        if (klonv.eq.1) then
412          if(isnoSV(1).ge.2                   .and.
413     .         TsisSV(1,max(1,isnoSV(1)))<273.  .and.
414     .         ro__SV(1,max(1,isnoSV(1)))<500.  .and.
415     .         eta_SV(1,max(1,isnoSV(1)))<epsi) then
416C +                       **********
417                     call SISVAT_BSn
418          endif
419         else
420                     call SISVAT_BSn
421C +                       **********
422        endif
423
424
425
426
427
428
429
430! Calculate threshold erosion velocity for next time step
431! Unlike in sisvat, computation is of threshold velocity made here (instead of sisvaesbl)
432! since we do not use sisvatesbl for the coupling with LMDZ
433
434C +--Computation of threshold friction velocity for snow erosion
435C ---------------------------------------------------------------
436
437        rCd10n =  1. / 26.5 ! Vt / u*t = 26.5
438                     ! Budd et al. 1965, Antarct. Res. Series Fig.13
439                     ! ratio developped during assumed neutral conditions
440 
441
442C +--Snow Properties
443C +  ~~~~~~~~~~~~~~~
444
445        DO ikl = 1,knonv
446
447          isn      =  isnoSV(ikl)
448
449
450 
451          DendOK   =  max(zero,sign(unun,epsi-G1snSV(ikl,isn)  ))  !
452          SaltOK   =  min(1   , max(istdSV(2)-istoSV(ikl,isn),0))  !
453          MeltOK   =     (unun                                     !
454     .             -max(zero,sign(unun,TfSnow-epsi                 !
455     .             -TsisSV(ikl,isn)  )))                           ! Melting Snow
456     .             *  min(unun,DendOK                              !
457     .                  +(1.-DendOK)                               !
458     .                      *sign(unun,     G2snSV(ikl,isn)-1.0))  ! 1.0 for 1mm
459          SnowOK   =  min(1   , max(isnoSV(ikl)      +1 -isn ,0))  ! Snow Switch
460 
461          G1snSV(ikl,isn) =      SnowOK *    G1snSV(ikl,isn)
462     .                  + (1.- SnowOK)*min(G1snSV(ikl,isn),G1_dSV)
463          G2snSV(ikl,isn) =      SnowOK *    G2snSV(ikl,isn)
464     .                  + (1.- SnowOK)*min(G2snSV(ikl,isn),G1_dSV)
465 
466          SaltOK   =  min(unun, SaltOK + MeltOK) * SnowOK
467 
468 
469C +--Mobility Index (Guyomarc'h & Merindol 1997, Ann.Glaciol.)
470C +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
471          SaltM1   = -0.750e-2 * G1snSV(ikl,isn)
472     .             -0.500e-2 * G2snSV(ikl,isn)+ 0.500e00 !dendritic case
473C +     CAUTION:  Guyomarc'h & Merindol Dendricity Sign is +
474C +     ^^^^^^^^                    MAR Dendricity Sign is -
475          SaltM2   = -0.833d-2 * G1snSV(ikl,isn)
476     .             -0.583d-2 * G2snSV(ikl,isn)+ 0.833d00 !non-dendritic case
477 
478c       SaltMo   = (DendOK   * SaltM1 + (1.-DendOK) *     SaltM2       )
479          SaltMo   = 0.625 !SaltMo pour d=s=0.5
480 
481!weighting SaltMo with surface snow density (Vionnet et al. 2012)
482cc#AE   FacRho   = 1.25 - 0.0042 * ro__SV(ikl,isn)
483cc#AE   SaltMo   = 0.34 * SaltMo + 0.66 * FacRho !needed for polar snow
484          MIN_Mo   =  0.
485c       SaltMo   =  max(SaltMo,MIN_Mo)
486c       SaltMo   =  SaltOK   * SaltMo + (1.-SaltOK) * min(SaltMo,SaltMx)
487c #TUNE SaltMo   =  SaltOK   * SaltMo - (1.-SaltOK) *     0.9500
488          SaltMo   =  max(SaltMo,epsi-unun)
489 
490C +--Influence of Density on Threshold Shear Stress
491C +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
492          Por_BS =  1. - 300. / ro_Ice
493          ShearS = Por_BS / (1.-Por_BS)
494C +...         SheaBS =  Arg(sqrt(shear = max shear stress in snow)):
495C +            shear  =  3.420d00 * exp(-(Por_BS      +Por_BS)
496C +  .                                  /(unun        -Por_BS))
497C +            SheaBS :  see de Montmollin         (1978),
498C +                      These Univ. Sci. Medic. Grenoble, Fig. 1 p. 124
499 
500C +--Snow Drift Index (Guyomarc'h & Merindol 1997, Ann.Glaciol.)
501C +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
502          ArguSi      =     -0.085 *us__SV(ikl)/rCd10n
503!V=u*/sqrt(CD) eqs 2 to 4 Gallee et al. 2001
504 
505          SaltSI(ikl,isn) = -2.868 * exp(ArguSi) + 1 + SaltMo
506 
507
508C +--Threshold Friction Velocity
509C +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~
510          if(ro__SV(ikl,isn)>300.) then
511             Por_BS      =  1.000       - ro__SV(ikl,isn)     /ro_Ice
512          else
513             Por_BS      =  1.000  - 300. /ro_Ice
514          endif
515 
516          ShearX =  Por_BS/max(epsi,1.-Por_BS)
517          Fac_Mo = exp(-ShearX+ShearS)
518C +     Gallee et al., 2001    eq 5, p5
519 
520          if (usth_param .eq. "gal") then
521            Salt_us   =   (log(2.868) - log(1 + SaltMo)) * rCd10n/0.085
522            Salt_us   = Salt_us * Fac_Mo
523C +...  Salt_us   :  Extension of  Guyomarc'h & Merindol 1998 with
524C +...              de Montmollin (1978). Gallee et al. 2001
525          endif
526 
527          if (usth_param .eq. "lis") then !Liston et al. 2007
528            if(ro__SV(ikl,isn)>300.) then
529              Salt_us   = 0.005*exp(0.013*ro__SV(ikl,isn))
530            else
531              Salt_us   = 0.01*exp(0.003*ro__SV(ikl,isn))
532            endif
533          endif
534 
535          SnowOK   =  1 -min(1,iabs(isn-isnoSV(ikl))) !Switch new vs old snow
536 
537          usthSV(ikl) =     SnowOK *   (Salt_us)
538     .                + (1.-SnowOK)*    usthSV(ikl)
539 
540        END DO
541
542
543 
544!  Feeback between blowing snow turbulent Scale  u* (commented here
545!  since ustar is an input variable (not in/out) of inlandsis)
546!  -----------------------------------------------------------------
547
548
549!           VVa_OK      =  max(0.000001,       VVaSBL(ikl))
550!           sss__N      =  vonkar      *       VVa_OK
551!           sss__F      = (sqrCm0(ikl) - psim_z + psim_0)
552!           usuth0      =  sss__N /sss__F                ! u* if NO Blow. Snow
553 
554!           sss__G      =  0.27417     * gravit
555 
556! !  ______________               _____
557! !  Newton-Raphson (! Iteration, BEGIN)
558! !  ~~~~~~~~~~~~~~               ~~~~~
559!           DO iit=1,nit
560!           sss__K      =  gravit      * r_Turb * A_Turb *za__SV(ikl)
561!      .                                     *rCDmSV(ikl)*rCDmSV(ikl)
562!      .                           /(1.+0.608*QaT_SV(ikl)-qsnoSV(ikl))
563!           us_127      =  exp(    SblPom *log(us__SV(ikl)))
564!           us_227      =  us_127         *    us__SV(ikl)
565!           us_327      =  us_227         *    us__SV(ikl)
566!           us_427      =  us_327         *    us__SV(ikl)
567!           us_527      =  us_427         *    us__SV(ikl)
568 
569!           us__SV(ikl) =  us__SV(ikl)
570!      .    - (  us_527     *sss__F     /sss__N
571!      .      -  us_427
572!      .      -  us_227     *qsnoSV(ikl)*sss__K
573!      .      + (us__SV(ikl)*us__SV(ikl)-usthSV(ikl)*usthSV(ikl))/sss__G)
574!      .     /(  us_427*5.27*sss__F     /sss__N
575!      .      -  us_327*4.27
576!      .      -  us_127*2.27*qsnoSV(ikl)*sss__K
577!      .      +  us__SV(ikl)*2.0                                 /sss__G)
578 
579!           us__SV(ikl)= min(us__SV(ikl),usuth0)
580!           us__SV(ikl)= max(us__SV(ikl),epsi  )
581!           rCDmSV(ikl)=     us__SV(ikl)/VVa_OK
582! ! #AE     sss__F     =     vonkar     /rCDmSV(ikl)
583!           ENDDO
584 
585! !  ______________               ___
586! !  Newton-Raphson (! Iteration, END  )
587! !  ~~~~~~~~~~~~~~               ~~~
588 
589!           us_127      =  exp(    SblPom *log(us__SV(ikl)))
590!           us_227      =  us_127         *    us__SV(ikl)
591 
592! !  Momentum            Turbulent Scale  u*: 0-Limit in case of no Blow. Snow
593! !  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
594!           dusuth      =  us__SV(ikl) - usthSV(ikl)       ! u* - uth*
595!           signus      =  max(sign(unun,dusuth),zero)     ! 1 <=> u* - uth* > 0
596!           us__SV(ikl) =                                  !
597!      .                   us__SV(ikl)  *signus  +         ! u* (_BS)
598!      .                   usuth0                          ! u* (nBS)
599!      .                            *(1.-signus)           !       
600
601
602
603
604!  Blowing Snow        Turbulent Scale ss*
605!  ---------------------------------------
606 
607        hSalSV(ikl) = 8.436e-2  * us__SV(ikl)**SblPom
608 
609        if (qsalt_param .eq. "pom") then
610          qSalSV(ikl) = (us__SV(ikl)**2 - usthSV(ikl)**2) *signus
611     .               / (hSalSV(ikl) * gravit * us__SV(ikl) * 3.25)
612        endif
613 
614        if (qsalt_param .eq. "bin") then
615          qSalSV(ikl) = (us__SV(ikl) * us__SV(ikl)
616     .                -usthSV(ikl) * usthSV(ikl))*signus
617     .                * 0.535 / (hSalSV(ikl) * gravit)
618        endif
619 
620        qSalSV(ikl) = qSalSV(ikl)/rht_SV(ikl) ! conversion kg/m3 to kg/kg
621 
622        ssstar      = rCDmSV(ikl) * (qsnoSV(ikl) - qSalSV(ikl))
623     .              * r_Turb !Bintanja 2000, BLM
624!r_Turb compensates for an overestim. of the blown snow part. fall velocity
625 
626        uss_SV(ikl) = min(zero    , us__SV(ikl) *ssstar)
627        uss_SV(ikl) = max(-0.0001 , uss_SV(ikl))   
628
629
630
631
632        ENDIF   ! BloMod
633 
634C + ------------------------------------------------------
635C +--Buffer Layer
636C +  -----------------------------------------------------
637 
638          DO ikl=1,knonv
639c  BufsSV(ikl) [mm w.e.] i.e, i.e., [kg/m2]
640            d_Bufs      =  max(dsn_SV(ikl) *dt__SV,0.)  !
641            dsn_SV(ikl) =      0.                       !
642            Bufs_N      =      BufsSV(ikl) +d_Bufs      !
643 
644 
645C +--Snow Density
646C +  ^^^^^^^^^^^^
647            Polair      =      zero
648c #NP       Polair      =  max(zero,                    !
649c #NP.                         sign(unun,TaPole         !
650c #NP.                                  -TaT_SV(ikl)))  !
651            Polair      =  max(zero,                    !
652     .                         sign(unun,TaPole         !
653     .                                  -TaT_SV(ikl)))  !
654            Buf_ro      =  max( rosMin,                 ! Fallen Snow Density
655     .      roSn_1+roSn_2*     (TaT_SV(ikl)-TfSnow)     ! [kg/m3]
656     .            +roSn_3*sqrt( VV__SV(ikl)))           ! Pahaut    (CEN), Etienne: use wind speed at first model level instead of 10m wind
657c #NP       BufPro      =  max( rosMin,                 ! Fallen Snow Density
658c #NP.         104. *sqrt( max( VV10SV(ikl)-6.0,0.0)))  ! Kotlyakov (1961)
659 
660!          C.Agosta option for snow density, same as for BS i.e.
661!          is_ok_density_kotlyakov=.false.
662c #BS       density_kotlyakov = .false.  !C.Amory BS 2018
663C + ...     Fallen Snow Density, Adapted for Antarctica
664            if (is_ok_density_kotlyakov) then
665                tt_tmp = TaT_SV(ikl)-TfSnow
666                !vv_tmp = VV10SV(ikl)
667                vv_tmp=VV__SV(ikl) ! Etienne: use wind speed at first model level instead of 10m wind
668C + ...         [ A compromise between
669C + ...           Kotlyakov (1961) and Lenaerts (2012, JGR, Part1) ]
670                if (tt_tmp.ge.-10) then
671                  BufPro   =  max( rosMin,
672     .            104. *sqrt( max( vv_tmp-6.0,0.0))) ! Kotlyakov (1961)
673                else
674                  vv_virt = (tt_c*vv_tmp+vv_c*(tt_tmp+10))
675     .                     /(tt_c+tt_tmp+10)
676                  BufPro  = 104. *sqrt( max( vv_virt-6.0,0.0))
677                endif
678            else
679C + ...         [ density derived from observations of the first 50cm of
680C + ...           snow - cf. Rajashree Datta - and multiplied by 0.8 ]
681C + ...           C. Agosta, 2016-09
682cc #SD           BufPro = 149.2 + 6.84*VV10SV(ikl) + 0.48*Tsrfsv(ikl)
683cc #SD           BufPro = 125 + 14*VV10SV(ikl) + 0.6*Tsrfsv(ikl) !MAJ CK and CAm
684!                BufPro = 200 + 21 * VV10SV(ikl)!CK 29/07/19
685                 BufPro = 200 + 21 * VV__SV(ikl)!Etienne: use wind speed at first model level instead of 10m wind
686            endif
687 
688            Bros_N      = (1. - Polair) *   Buf_ro      ! Temperate Snow
689     .                        + Polair  *   BufPro      ! Polar     Snow
690 
691            Bros_N = max( 20.,max(rosMin,  Bros_N))
692            Bros_N = min(400.,min(rosMax-1,Bros_N)) ! for dz_min in SISVAT_zSn
693 
694 
695!    Density of deposited blown snow
696!    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
697 
698         if (BloMod) then
699         Bros_N      = frsno
700         ro_new      = ro__SV(ikl,max(1,isnoSV(ikl)))
701         ro_new      = max(Bros_N,min(roBdSV,ro_new))
702         Fac         = 1-((ro__SV(ikl,max(1,isnoSV(ikl)))
703     .               -roBdSV)/(500.-roBdSV))
704         Fac         = max(0.,min(1.,Fac))
705         dsnbSV(ikl) = Fac*dsnbSV(ikl)
706         Bros_N      = Bros_N     * (1.0-dsnbSV(ikl))
707     .               + ro_new     *      dsnbSV(ikl)
708         endif
709
710 
711!    Time averaged Density of deposited blown Snow
712!    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
713           
714            BrosSV(ikl) =(Bros_N     *      d_Bufs      !
715     .                   +BrosSV(ikl)*      BufsSV(ikl))!
716     .                   /         max(epsi,Bufs_N)     !
717 
718 
719C +-- S.Falling Snow Properties (computed as in SISVAT_zAg)
720C +     ^^^^^^^^^^^^^^^^^^^^^^^
721            Buf_G1      =  max(-G1_dSV,                 ! Temperate Snow
722     .               min(Dendr1*VV__SV(ikl)-Dendr2,     !     Dendricity
723     .                   Dendr3                   ))    !
724            Buf_G2      =  min( Spher4,                 ! Temperate Snow
725     .               max(Spher1*VV__SV(ikl)+Spher2,     !     Sphericity
726     .                   Spher3                   ))    !
727! EV: now control buf_sph_pol and bug_siz_pol in physiq.def
728            Buf_G1      = (1. - Polair) *   Buf_G1      ! Temperate Snow
729     .                        + Polair  *   buf_sph_pol ! Polar Snow
730            Buf_G2      = (1. - Polair) *   Buf_G2      ! Temperate Snow
731     .                        + Polair  *   buf_siz_pol ! Polar Snow
732                G1      =                   Buf_G1      ! NO  Blown Snow
733                G2      =                   Buf_G2      ! NO  Blown Snow
734
735
736
737            IF (BloMod) THEN
738
739!     S.1. Meme  Type  de Neige  / same Grain Type
740!          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
741
742           SameOK  =  max(zero,
743     .         sign(unun,    Buf_G1             *G1_dSV
744     .                            - eps_21                    ))
745           G1same  = ((1.0-dsnbSV(ikl))*Buf_G1+dsnbSV(ikl) *G1_dSV)
746           G2same  = ((1.0-dsnbSV(ikl))*Buf_G2+dsnbSV(ikl) *ADSdSV)
747!           Blowing Snow Properties:                         G1_dSV, ADSdSV
748 
749!     S.2. Types differents / differents Types
750!          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
751           typ__1  =  max(zero,sign(unun,epsi-Buf_G1))   ! =1.=> Dendritic
752           zroNEW  =     typ__1  *(1.0-dsnbSV(ikl))      ! fract.Dendr.Lay.
753     .            + (1.-typ__1) *     dsnbSV(ikl)       !
754           G1_NEW  =     typ__1  *Buf_G1                 ! G1 of Dendr.Lay.
755     .            + (1.-typ__1) *G1_dSV                 !
756           G2_NEW  =     typ__1  *Buf_G2                 ! G2 of Dendr.Lay.
757     .            + (1.-typ__1) *ADSdSV                 !
758           zroOLD  = (1.-typ__1) *(1.0-dsnbSV(ikl))      ! fract.Spher.Lay.
759     .            +     typ__1  *     dsnbSV(ikl)       !
760           G1_OLD  = (1.-typ__1) *Buf_G1                 ! G1 of Spher.Lay.
761     .            +     typ__1  *G1_dSV                 !
762           G2_OLD  = (1.-typ__1) *Buf_G2                 ! G2 of Spher.Lay.
763     .            +     typ__1  *ADSdSV                 !
764           SizNEW  =    -G1_NEW  *DDcdSV/G1_dSV          ! Size  Dendr.Lay.
765     .            +(1.+G1_NEW         /G1_dSV)          !
766     .                  *(G2_NEW  *DScdSV/G1_dSV        !
767     .            +(1.-G2_NEW         /G1_dSV)*DFcdSV)  !
768           SphNEW  =     G2_NEW         /G1_dSV          ! Spher.Dendr.Lay.
769           SizOLD  =     G2_OLD                          ! Size  Spher.Lay.
770           SphOLD  =     G1_OLD         /G1_dSV          ! Spher.Spher.Lay.
771           Siz_av  =     (zroNEW*SizNEW+zroOLD*SizOLD)   ! Averaged Size
772           Sph_av  = min( zroNEW*SphNEW+zroOLD*SphOLD    !
773     .                 ,  unun)                         ! Averaged Sphericity
774           Den_av  = min((Siz_av -(    Sph_av *DScdSV    !
775     .            +(1.-Sph_av)*DFcdSV))                 !
776     .            / (DDcdSV -(    Sph_av *DScdSV        !
777     .            +(1.-Sph_av)*DFcdSV))                 !
778     .                   ,  unun)                       !
779           DendOK  = max(zero,                           !
780     .                    sign(unun,     Sph_av *DScdSV   ! Small   Grains
781     .                              +(1.-Sph_av)*DFcdSV   ! Faceted Grains
782     .                              -    Siz_av        )) !
783C +...      REMARQUE: le  type moyen (dendritique ou non) depend
784C +         ^^^^^^^^  de la  comparaison avec le diametre optique
785C +                   d'une neige recente de   dendricite nulle
786C +...      REMARK:   the mean type  (dendritic   or not) depends
787C +         ^^^^^^    on the comparaison with the optical diameter
788C +                   of a recent snow    having zero dendricity
789 
790           G1diff  =(   -DendOK *Den_av
791     .            +(1.-DendOK)*Sph_av) *G1_dSV
792           G2diff  =     DendOK *Sph_av  *G1_dSV
793     .            +(1.-DendOK)*Siz_av
794           G1      =     SameOK *G1same
795     .            +(1.-SameOK)*G1diff
796           G2      =     SameOK *G2same
797     .            +(1.-SameOK)*G2diff
798           ENDIF
799
800
801 
802!     S.1. Meme  Type  de Neige  / same Grain Type
803!          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
804            SameOK  =  max(zero,
805     .                     sign(unun,    Buf_G1 *BG1sSV(ikl)
806     .                                 - eps_21                    ))
807            G1same  = (d_Bufs*Buf_G1+BufsSV(ikl)*BG1sSV(ikl))
808     .                     /max(epsi,Bufs_N)
809            G2same  = (d_Bufs*Buf_G2+BufsSV(ikl)*BG2sSV(ikl))
810     .                     /max(epsi,Bufs_N)
811 
812!     S.2. Types differents / differents Types
813!          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
814
815            typ__1  =  max(zero,sign(unun,epsi-Buf_G1))   ! =1.=> Dendritic
816            zroNEW  =(    typ__1  *d_Bufs                 ! fract.Dendr.Lay.
817     .              + (1.-typ__1) *BufsSV(ikl))           !
818     .                   /max(epsi,Bufs_N)                !
819            G1_NEW  =     typ__1  *Buf_G1                 ! G1 of Dendr.Lay.
820     .              + (1.-typ__1) *BG1sSV(ikl)            !
821            G2_NEW  =     typ__1  *Buf_G2                 ! G2 of Dendr.Lay.
822     .              + (1.-typ__1) *BG2sSV(ikl)            !
823            zroOLD  =((1.-typ__1) *d_Bufs                 ! fract.Spher.Lay.
824     .              +     typ__1  *BufsSV(ikl))           !
825     .                   /max(epsi,Bufs_N)                !
826            G1_OLD  = (1.-typ__1) *Buf_G1                 ! G1 of Spher.Lay.
827     .              +     typ__1  *BG1sSV(ikl)            !
828            G2_OLD  = (1.-typ__1) *Buf_G2                 ! G2 of Spher.Lay.
829     .              +     typ__1  *BG2sSV(ikl)            !
830            SizNEW  =    -G1_NEW  *DDcdSV/G1_dSV          ! Size  Dendr.Lay.
831     .               +(1.+G1_NEW         /G1_dSV)         !
832     .                  *(G2_NEW  *DScdSV/G1_dSV          !
833     .               +(1.-G2_NEW         /G1_dSV)*DFcdSV) !
834            SphNEW  =     G2_NEW         /G1_dSV          ! Spher.Dendr.Lay.
835            SizOLD  =     G2_OLD                          ! Size  Spher.Lay.
836            SphOLD  =     G1_OLD         /G1_dSV          ! Spher.Spher.Lay.
837            Siz_av  =   ( zroNEW  *SizNEW+zroOLD*SizOLD)  ! Averaged Size
838            Sph_av = min( zroNEW  *SphNEW+zroOLD*SphOLD   !
839     .                  ,   unun                       )  ! Averaged Sphericity
840            Den_av = min((Siz_av  - (    Sph_av *DScdSV   !
841     .                              +(1.-Sph_av)*DFcdSV)) !
842     .                 / (DDcdSV  - (    Sph_av *DScdSV   !
843     .                              +(1.-Sph_av)*DFcdSV)) !
844     .                  ,   unun                         )!
845            DendOK  = max(zero,                           !
846     .                    sign(unun,     Sph_av *DScdSV   ! Small   Grains
847     .                              +(1.-Sph_av)*DFcdSV   ! Faceted Grains
848     .                              -    Siz_av        )) !
849C +...      REMARQUE: le  type moyen (dendritique ou non) depend
850C +         ^^^^^^^^  de la  comparaison avec le diametre optique
851C +                   d'une neige recente de   dendricite nulle
852C +...      REMARK:   the mean type  (dendritic   or not) depends
853C +         ^^^^^^    on the comparaison with the optical diameter
854C +                   of a recent snow    having zero dendricity
855 
856            G1diff  =(   -DendOK *Den_av
857     .               +(1.-DendOK)*Sph_av) *G1_dSV
858            G2diff  =     DendOK *Sph_av  *G1_dSV
859     .               +(1.-DendOK)*Siz_av
860            G1      =     SameOK *G1same
861     .               +(1.-SameOK)*G1diff
862            G2      =     SameOK *G2same
863     .               +(1.-SameOK)*G2diff
864 
865            BG1sSV(ikl) =                       G1      !
866     .                  *       Bufs_N/max(epsi,Bufs_N) !
867            BG2sSV(ikl) =                       G2      !
868     .                  *       Bufs_N/max(epsi,Bufs_N) !
869 
870
871C +--Update of Buffer Layer Content & Decision about creating a new snow layer
872C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
873            BufsSV(ikl) =       Bufs_N                  !     [mm w.e.]
874            NLaysv(ikl) = min(unun,                     !
875     .                    max(zero,                     ! Allows to create
876     .                        sign(unun,BufsSV(ikl)     ! a new snow Layer
877     .                                 -SMndSV     ))   ! if Buffer > SMndSV
878     .                   *max(zero,                     ! Except if * Erosion
879     .                        sign(unun,0.50            ! dominates
880     .                                 -dsnbSV(ikl)))   !
881     .                   +max(zero,                     ! Allows to create
882     .                        sign(unun,BufsSV(ikl)     ! a new snow Layer
883     .                                 -SMndSV*3.00)))  ! is Buffer > SMndSV*3
884            Bdzssv(ikl) = 1.e-3*BufsSV(ikl)*ro_Wat      ! [mm w.e.] -> [m w.e.]
885     .                            /max(epsi,BrosSV(ikl))!& [m w.e.] -> [m]
886 
887 
888          END DO
889 
890
891
892! Snow Pack Discretization(option XF in MAR)
893! ==========================================
894
895         
896      if (discret_xf.AND.klonv.eq.1) then
897
898       if(isnoSV(1).ge.1.or.NLaysv(1).ge.1) then
899C +          **********
900         call SISVAT_zSn
901C +          **********
902       endif
903      else
904C +          **********
905        call SISVAT_zSn
906C +          **********
907      endif
908 
909C +          **********
910! #ve   call SISVAT_wEq('_zSn  ',0)
911C +          **********
912
913! Add a new Snow Layer
914! ====================
915
916          DO ikl=1,knonv
917
918            isnoSV(ikl)     = isnoSV(ikl)         +NLaysv(ikl)
919            isn             = isnoSV(ikl)
920            dzsnSV(ikl,isn) = dzsnSV(ikl,isn) * (1-NLaysv(ikl))
921     .                      + Bdzssv(ikl)     *    NLaysv(ikl)
922            TsisSV(ikl,isn) = TsisSV(ikl,isn) * (1-NLaysv(ikl))
923     .                  + min(TaT_SV(ikl),Tf_Sno) *NLaysv(ikl)
924            ro__SV(ikl,isn) = ro__SV(ikl,isn) * (1-NLaysv(ikl))
925     .                      + Brossv(ikl)     *    NLaysv(ikl)
926            eta_SV(ikl,isn) = eta_SV(ikl,isn) * (1-NLaysv(ikl))   ! + 0.
927            agsnSV(ikl,isn) = agsnSV(ikl,isn) * (1-NLaysv(ikl))   ! + 0.
928            G1snSV(ikl,isn) = G1snSV(ikl,isn) * (1-NLaysv(ikl))
929     .                      + BG1ssv(ikl)     *    NLaysv(ikl)
930            G2snSV(ikl,isn) = G2snSV(ikl,isn) * (1-NLaysv(ikl))
931     .                      + BG2ssv(ikl)     *    NLaysv(ikl)
932            istoSV(ikl,isn) = istoSV(ikl,isn) * (1-NLaysv(ikl))
933     .   + max(zer0,sign(un_1,TaT_SV(ikl)
934     .                       -Tf_Sno-eps_21)) *    istdSV(2)
935     .                                        *    NLaysv(ikl)
936            BufsSV(ikl)     = BufsSV(ikl)     * (1-NLaysv(ikl))
937            NLaysv(ikl)     = 0
938
939
940          END DO
941
942
943! Snow Pack Thickness
944! -------------------
945
946        DO ikl=1,knonv
947            z_snsv(ikl)     = 0.0
948        END DO
949        DO   isn=1,nsno
950          DO ikl=1,knonv
951            z_snsv(ikl)     = z_snsv(ikl) + dzsnSV(ikl,isn)
952            zzsnsv(ikl,isn) = z_snsv(ikl)
953          END DO
954        END DO
955
956
957
958      END IF  ! SnoMod
959
960
961
962! Soil Albedo: Soil Humidity Correction
963! ==========================================
964
965!         REFERENCE: McCumber and Pielke (1981), Pielke (1984)
966!         ^^^^^^^^^
967          DO ikl=1,knonv
968            albssv(ikl) =
969     .      alb0SV(ikl) *(1.0-min(half,eta_SV(       ikl,0)
970     .                                /etadSV(isotSV(ikl))))
971!         REMARK:    Albedo of Water Surfaces (isotSV=0):
972!         ^^^^^^     alb0SV := 2  X  effective value, while
973!                    eta_SV :=          etadSV
974          END DO
975
976
977! Snow Pack Optical Properties
978! ============================
979
980      IF (SnoMod)                                                 THEN
981
982!            ******
983        call SnOptP(jjtime)
984!            ******
985
986      ELSE
987        DO ikl=1,knonv
988          sEX_sv(ikl,1) = 1.0
989          sEX_sv(ikl,0) = 0.0
990          albisv(ikl)   = albssv(ikl)
991        END DO
992      END IF
993
994
995
996! Soil optical properties
997! =============================
998!Etienne: as in inlandis we do not call vgopt, we need to define
999!the albedo alb_SV and to calculate the
1000!absorbed Solar Radiation by Surfac (Normaliz)[-] SoSosv
1001
1002 
1003      DO ikl=1,klonv
1004
1005            e_pRad = 2.5   *  coszSV(ikl)       ! exponential argument,
1006                                                ! V/nIR radiation partitioning,
1007                                                ! DR97, 2, eqn (2.53) & (2.54)
1008            e1pRad = 1.-exp(-e_pRad)            ! exponential, V/nIR Rad. Part.
1009            exdRad= 1.
1010 
1011! Visible Part of the Solar Radiation Spectrum (V,   0.4--0.7mi.m)
1012            A_Rad0 =      0.25 + 0.697 * e1pRad ! Absorbed    Vis. Radiation
1013            absg_V = (1.-albisv(ikl))*(A_Rad0*exdRad)  !
1014 
1015! Near-IR Part of the Solar Radiation Spectrum (nIR, 0.7--2.8mi.m)
1016 
1017            A_Rad0 =      0.80 + 0.185 * e1pRad ! Absorbed    nIR. Radiation
1018            absgnI = (1.-albisv(ikl))*(A_Rad0*exdRad)  !
1019
1020            SoSosv(ikl) = (absg_V+absgnI)*0.5d0
1021
1022            alb_SV(ikl) = albisv(ikl)
1023
1024      END DO
1025 
1026!            **********
1027! #ve   call SISVAT_wEq('SnOptP',0)
1028!            **********
1029
1030
1031! Surface Emissivity (Etienne: simplified calculation for landice)
1032! =============================================================
1033!
1034       DO ikl=1,knonv
1035            LSnMsk     =     min( 1,isnoSV(ikl))
1036            Eso_sv(ikl)=  EmiSol*(1-LSnMsk)+EmiSno*LSnMsk  ! Sol+Sno Emissivity
1037            emi_SV(ikl)= EmiSol*(1-LSnMsk) + EmiSno*LSnMsk
1038        END DO
1039
1040
1041
1042
1043!  Upward IR (INPUT, from previous time step)
1044! ===================================================================
1045
1046        DO ikl=1,knonv
1047! #e1     Enrsvd(ikl) =    - IRs_SV(ikl)
1048           IRupsv(ikl) =      IRs_SV(ikl)
1049        END DO
1050
1051
1052! Turbulence
1053! ==========
1054
1055! Latent Heat of Vaporization/Sublimation
1056! ---------------------------------------
1057
1058        DO ikl=1,knonv
1059          SnoWat      =                     min(isnoSV(ikl),0)
1060          Lx_H2O(ikl) =
1061     .    (1.-SnoWat) * LhvH2O
1062     .  +     SnoWat  *(LhsH2O * (1.-eta_SV(ikl,isnoSV(ikl)))
1063     .                 +LhvH2O *     eta_SV(ikl,isnoSV(ikl)) )
1064        END DO
1065
1066
1067
1068
1069! Aerodynamic Resistance (calculated from drags given by LMDZ)
1070! Commented because already calculated in surf_inlandsis_mod
1071! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1072!       DO ikl=1,knonv
1073!          ram_sv(ikl) = 1./(cdM_SV(ikl)*max(VV__SV(ikl),eps6))
1074!          rah_sv(ikl) = 1./(cdH_SV(ikl)*max(VV__SV(ikl),eps6))
1075!        END DO
1076
1077
1078
1079! Soil   Energy Balance
1080! =====================
1081
1082
1083      if (iflag_temp_inlandsis .eq. 0) then
1084
1085       call SISVAT_TSo
1086
1087      else
1088        DO ikl=1,knonv
1089        Tsf_SV(ikl)=Tsrfsv(ikl)
1090        END DO
1091
1092       call SISVAT_TS2
1093
1094      end if
1095
1096
1097!            **********
1098! #ve   call SISVAT_wEq('_TSo  ',0)
1099!            **********
1100
1101
1102
1103! Soil Water     Potential
1104! ------------------------
1105
1106      DO   isl=-nsol,0
1107        DO ikl=1,knonv
1108          ist             =     isotSV(ikl)        ! Soil Type
1109          psi_sv(ikl,isl) =     psidSV(ist)        ! DR97, Eqn.(3.34)
1110     .  *(etadSV(ist) /max(eps6,eta_SV(ikl,isl)))  !
1111     .  **bCHdSV(ist)                              !
1112
1113
1114! Soil Hydraulic Conductivity
1115! ---------------------------
1116
1117          Khydsv(ikl,isl) =    s2__SV(ist)         ! DR97, Eqn.(3.35)
1118     .  *(eta_SV(ikl,isl)**(2.*bCHdSV(ist)+3.))    ! 
1119        END DO
1120      END DO
1121
1122
1123! Melting / Refreezing in the Snow Pack
1124! =====================================
1125
1126      IF (SnoMod)                                                 THEN
1127
1128!            **********
1129        call SISVAT_qSn
1130!            **********
1131
1132!            **********
1133! #ve   call SISVAT_wEq('_qSn  ',0)
1134!            **********
1135
1136
1137! Snow Pack Thickness
1138! -------------------
1139
1140          DO ikl=1,knonv
1141            z_snsv(ikl)     = 0.0
1142          END DO
1143        DO   isn=1,nsno
1144          DO ikl=1,knonv
1145            z_snsv(ikl)     = z_snsv(ikl) + dzsnSV(ikl,isn)
1146            zzsnsv(ikl,isn) = z_snsv(ikl)
1147          END DO
1148        END DO
1149
1150
1151! Energy in Excess is added to the first Soil Layer
1152! -------------------------------------------------
1153        DO ikl=1,knonv
1154            z_snsv(ikl)   = max(zer0,
1155     .                          sign(un_1,eps6-z_snsv(ikl)))
1156            TsisSV(ikl,0) = TsisSV(ikl,0)    + EExcsv(ikl)
1157     .                                       /(rocsSV(isotSV(ikl))
1158     .                                        +rcwdSV*eta_SV(ikl,0))
1159            EExcsv(ikl)   = 0.
1160        END DO
1161
1162
1163      END IF
1164
1165
1166! Soil   Water  Balance
1167! =====================
1168
1169!            **********
1170        call SISVAT_qSo
1171! #m0.                 (Wats_0,Wats_1,Wats_d)
1172!            **********
1173
1174
1175! Surface Fluxes
1176! =====================
1177
1178        DO ikl=1,knonv
1179         IRdwsv(ikl)=IRd_SV(ikl)*Eso_sv(ikl) ! Downward IR
1180!         IRdwsv(ikl)=tau_sv(ikl) *IRd_SV(ikl)*Eso_sv(ikl) ! Downward IR
1181!    .          +(1.0-tau_sv(ikl))*IRd_SV(ikl)*Evg_sv(ikl) !  ! Etienne, remove vegetation component
1182          IRupsv(ikl) =      IRupsv(ikl)                   ! Upward   IR
1183          IRu_SV(ikl) =     -IRupsv(ikl)                   ! Upward   IR
1184     .                      +IRd_SV(ikl)                   ! (effective)
1185     .                      -IRdwsv(ikl)                   ! (positive)
1186
1187          TBr_sv(ikl) =sqrt(sqrt(IRu_SV(ikl)/StefBo))      ! Brightness
1188!                                                          ! Temperature
1189          uts_SV(ikl) =     (HSv_sv(ikl) +HSs_sv(ikl))     ! u*T*
1190     .                     /(rhT_SV(ikl) *cp)          !
1191          uqs_SV(ikl) =     (HLv_sv(ikl) +HLs_sv(ikl))     ! u*q*
1192     .                     /(rhT_SV(ikl) *LhvH2O)          !
1193          LMO_SV(ikl) = TaT_SV(ikl)*(us__SV(ikl)**3)
1194     .                     /gravit/uts_SV(ikl)/vonKrm      ! MO length
1195     
1196
1197! Surface Temperature
1198! ^^^^^^^^^^^^^^^^^^^^
1199
1200          IF (iflag_tsurf_inlandsis .EQ. 0) THEN   
1201
1202            Tsrfsv(ikl) =TsisSV(ikl,isnoSV(ikl))
1203
1204          ELSE IF (iflag_tsurf_inlandsis .GT. 0) THEN
1205! Etienne: extrapolation from the two uppermost levels:
1206
1207         if (isnoSV(ikl) >=2) then
1208           zm1=-dzsnSV(ikl,isnoSV(ikl))/2.
1209           zm2=-(dzsnSV(ikl,isnoSV(ikl)) + dzsnSV(ikl,isnoSV(ikl)-1)/2.)
1210         else if (isnoSV(ikl) .EQ. 1) then
1211           zm1=-dzsnSV(ikl,isnoSV(ikl))/2.
1212           zm2=-(dzsnSV(ikl,isnoSV(ikl))+dz_dSV(0)/2.)
1213         else
1214           zm1=-dz_dSV(0)/2.
1215           zm2=-(dz_dSV(0)+dz_dSV(-1)/2.)
1216
1217         end if
1218
1219           coefslope=(TsisSV(ikl,isnoSV(ikl))-TsisSV(ikl,isnoSV(ikl)-1))
1220     .               /(zm1-zm2)
1221           Tsrfsv(ikl)=TsisSV(ikl,isnoSV(ikl))+coefslope*(0. - zm1)
1222
1223
1224         ELSE !(default)
1225
1226           Tsrfsv(ikl) =TsisSV(ikl,isnoSV(ikl))
1227
1228         END IF
1229
1230
1231         END DO
1232
1233! Snow Pack Properties (sphericity, dendricity, size)
1234! ===================================================
1235
1236      IF (SnoMod)                                                 THEN
1237
1238      if (discret_xf .AND. klonv.eq.1) then
1239      if(isnoSV(1).ge.1) then
1240C +          **********
1241      call SISVAT_GSn
1242C +          **********
1243      endif
1244      else
1245C +          **********
1246        call SISVAT_GSn
1247C +          **********
1248      endif
1249
1250
1251      END IF
1252
1253
1254! Roughness Length for next time step
1255!====================================
1256
1257! Note that in INLANDSIS, we treat only ice covered surfaces so calculation
1258! of z0 is much simpler (no subgrid fraction of ocean or land)
1259! old calculations are commented below
1260
1261
1262C +--Roughness Length for Momentum
1263C +  -----------------------------
1264
1265! ETIENNE WARNING: changes have been made wrt original SISVAT
1266 
1267C +--Land+Sea-Ice / Ice-free Sea Mask
1268C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1269        DO ikl=1,knonv
1270          IcIndx(ikl) = 0
1271        ENDDO
1272        DO isn=1,nsno
1273        DO ikl=1,knonv
1274
1275          IcIndx(ikl) = max(IcIndx(ikl),
1276     .                  isn*max(0,
1277     .                  sign(1,
1278     .                  int(ro__SV(ikl,isn)-900.))))
1279        ENDDO
1280        ENDDO
1281 
1282        DO ikl=1,knonv
1283          LISmsk    =     1. ! in inlandsis, land only
1284          IceMsk    =     max(0,sign(1   ,IcIndx(ikl)-1)  )
1285          SnoMsk    = max(min(isnoSV(ikl)-iiceSV(ikl),1),0)
1286
1287
1288C +--Z0 Smooth Regime over Snow (Andreas 1995, CRREL Report 95-16, p. 8)
1289C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^
1290          Z0m_nu =       5.e-5 ! z0s~(10-d)*exp(-vonkar/sqrt(1.1e-03))
1291
1292C +--Z0 Saltat.Regime over Snow (Gallee  et al., 2001, BLM 99 (19) p.11)
1293C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^
1294
1295          u2star =       us__SV(ikl) *us__SV(ikl)
1296          Z0mBSn =       u2star      *0.536e-3   -  61.8e-6
1297          Z0mBSn =   max(Z0mBS0      ,Z0mBSn)
1298
1299C +--Z0 Smooth + Saltat. Regime
1300C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^
1301          Z0enSV(ikl) =  Z0m_nu
1302     .                +  Z0mBSn
1303
1304       
1305! Calculation of snow roughness length
1306!=====================================
1307          IF (iflag_z0m_snow .EQ. 0) THEN
1308
1309          Z0m_Sn=prescribed_z0m_snow
1310
1311          ELSE IF (iflag_z0m_snow .EQ. 1) THEN
1312
1313          Z0m_Sn=Z0enSV(ikl)
1314
1315          ELSE IF (iflag_z0m_snow .EQ. 2) THEN                             
1316
1317C +--Rough   Snow Surface Roughness Length (Variable Sastrugi Height)
1318C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1319          A_Fact      =  1.0000        ! Andreas et al., 2004, p.4
1320                                       ! ams.confex.com/ams/pdfpapers/68601.pdf
1321 
1322! Parameterization of z0 dependance on Temperature (C. Amory, 2017)
1323! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1324! Z0=f(T) deduced from observations, Adelie Land, dec2012-dec2013
1325
1326         
1327          coefa = 0.1658 !0.1862 !Ant
1328          coefb = -50.3869 !-55.7718 !Ant
1329          ta1 = 253.15 !255. Ant
1330          ta2 = 273.15
1331          ta3 = 273.15+3
1332          z01 = exp(coefa*ta1 + coefb) !~0.2 ! ~0.25 mm
1333          z02 = exp(coefa*ta2 + coefb) !~6  !~7 mm
1334          z03 = z01
1335          coefc = log(z03/z02)/(ta3-ta2)
1336          coefd = log(z03)-coefc*ta3
1337
1338          if (TaT_SV(ikl) .lt. ta1) then
1339            Z0_obs = z01
1340          else if (TaT_SV(ikl).ge.ta1 .and. TaT_SV(ikl).lt.ta2) then
1341            Z0_obs = exp(coefa*TaT_SV(ikl) + coefb)
1342          else if (TaT_SV(ikl).ge.ta2 .and. TaT_SV(ikl).lt.ta3) then
1343            ! if st > 0, melting induce smooth surface
1344            Z0_obs = exp(coefc*TaT_SV(ikl) + coefd)
1345          else
1346            Z0_obs = z03
1347          endif
1348 
1349          Z0m_Sn=Z0_obs
1350
1351
1352          ELSE
1353
1354          Z0m_Sn=0.500e-3  ! default=0.500e-3m (tuning of MAR)
1355
1356          ENDIF
1357 
1358
1359
1360!          param = Z0_obs/1. ! param(s) | 1.(m/s)=TUNING
1361c #SZ     Z0Sa_N =                   (us__SV(ikl) -0.2)*param   ! 0.0001=TUNING
1362c #SZ.           * max(zero,sign(unun,TfSnow-eps9
1363c #SZ.                               -TsisSV(ikl , isnoSV(ikl))))
1364!!#SZ     Z0SaSi = max(zero,sign(unun,Z0Sa_N                  ))! 1 if erosion
1365c #SZ     Z0SaSi = max(zero,sign(unun,zero  -eps9 -uss_SV(ikl)))!
1366c #SZ     Z0Sa_N = max(zero,          Z0Sa_N)
1367c #SZ     Z0SaSV(ikl) =
1368c #SZ.             max(Z0SaSV(ikl)   ,Z0SaSV(ikl)
1369c #SZ.               + Z0SaSi*(Z0Sa_N-Z0SaSV(ikl))*exp(-dt__SV/43200.))
1370c #SZ.               -            min(dz0_SV(ikl) ,     Z0SaSV(ikl))
1371 
1372c #SZ     A_Fact      =               Z0SaSV(ikl) *  5.0/0.15   ! A=5 if h~10cm
1373C +...    CAUTION: The influence of the sastrugi direction is not yet included
1374 
1375c #SZ     Z0m_Sn =                    Z0SaSV(ikl)               !
1376c #SZ.                              - Z0m_nu                    !
1377 
1378C +--Z0 Saltat.Regime over Snow (Shao & Lin, 1999, BLM 91 (46)  p.222)
1379C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^
1380c #ZN     sqrrZ0 =       usthSV(ikl)/max( us__SV(ikl),0.001)
1381c #ZN     sqrrZ0 =                   min( sqrrZ0     ,0.999)
1382c #ZN     Z0mBSn =       0.55 *0.55 *exp(-sqrrZ0     *sqrrZ0)
1383c #ZN.                  *us__SV(ikl)*     us__SV(ikl)*grvinv*0.5
1384 
1385C +--Z0 Smooth + Saltat. Regime (Shao & Lin, 1999, BLM 91 (46)  p.222)
1386C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^
1387c #ZN     Z0enSV(ikl) = (Z0m_nu     **    sqrrZ0 )
1388c #ZN.                * (Z0mBSn     **(1.-sqrrZ0))
1389c #ZN     Z0enSV(ikl) =  max(Z0enSV(ikl), Z0m_nu)
1390 
1391
1392C +--Z0 Smooth Regime over Snow (Andreas etAl., 2004
1393C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^  ams.confex.com/ams/pdfpapers/68601.pdf)
1394c #ZA     Z0m_nu = 0.135*akmol  / max(us__SV(ikl) , epsi)
1395 
1396C +--Z0 Saltat.Regime over Snow (Andreas etAl., 2004
1397C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^  ams.confex.com/ams/pdfpapers/68601.pdf)
1398c #ZA     Z0mBSn = 0.035*u2star      *grvinv
1399 
1400C +--Z0 Smooth + Saltat. Regime (Andreas etAl., 2004
1401!    (      used by Erosion)     ams.confex.com/ams/pdfpapers/68601.pdf)
1402!    ^^^^^^^^^^^^^^^^^^^^^^^^^^
1403c #ZA     Z0enSV(ikl) =  Z0m_nu
1404c #ZA.                +  Z0mBSn
1405 
1406C +--Z0 Rough  Regime over Snow (Andreas etAl., 2004
1407C +  (.NOT. used by Erosion)     ams.confex.com/ams/pdfpapers/68601.pdf)
1408!    ^^^^^^^^^^^^^^^^^^^^^^^^^^
1409!!#ZA     u2star =      (us__SV(ikl) -0.1800)     / 0.1
1410!!#ZA     Z0m_Sn =A_Fact*Z0mBSn *exp(-u2star*u2star)
1411c #ZA     Z0m_90 =(10.-0.025*VVs_SV(ikl)/5.)
1412c #ZA.            *exp(-0.4/sqrt(.00275+.00001*max(0.,VVs_SV(ikl)-5.)))
1413c #ZA     Z0m_Sn =           DDs_SV(ikl)* Z0m_90 / 45.
1414c #ZA.         - DDs_SV(ikl)*DDs_SV(ikl)* Z0m_90 /(90.*90.)
1415
1416
1417
1418
1419C +--Z0  (Erosion)    over Snow (instantaneous)
1420C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^
1421          Z0e_SV(ikl) =  Z0enSV(ikl)
1422 
1423C +--Momentum  Roughness Length (Etienne: changes wrt original SISVAT)
1424C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                             
1425          Z0mnSV(ikl) =  Z0m_nu *(1-SnoMsk)                     ! Ice z0
1426     .                + (Z0m_Sn)*SnoMsk                         ! Snow Sastrugi Form and Snow Erosion
1427 
1428
1429C +--GIS  Roughness Length
1430C +  ^^^^^^^^^^^^^^^^^^^^^
1431c #GL     Z0mnSV(ikl) =
1432c #GL.      (1-LSmask(ikl)) *     Z0mnSV(ikl)
1433c #GL.    +    LSmask(ikl)  * max(Z0mnSV(ikl),max(Z0_GIM,
1434c #GL.                                            Z0_GIM+
1435c #GL.      (0.0032-Z0_GIM)*(ro__SV(ikl,isnoSV(ikl))-600.)   !
1436c #GL.                     /(920.00                 -600.))) !
1437 
1438C +--Mom. Roughness Length, Instantaneous
1439C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1440          Z0m_SV(ikl) =  Z0mnSV(ikl)                         ! Z0mnSV  instant.
1441 
1442 
1443C +--Roughness Length for Scalars
1444C +  ----------------------------
1445 
1446          Z0hnSV(ikl) =     Z0mnSV(ikl)/  7.4
1447 
1448          IF (is_ok_z0h_rn) THEN
1449
1450          rstar       =     Z0mnSV(ikl) * us__SV(ikl) / akmol
1451          rstar       = max(epsi,min(rstar,R_1000))
1452          alors       =          log(rstar)
1453          rstar0      = 1.250e0 * max(zero,sign(unun,0.135e0 - rstar))
1454     .                +(1.      - max(zero,sign(unun,0.135e0 - rstar)))
1455     .                *(0.149e0 * max(zero,sign(unun,2.500e0 - rstar))
1456     .                + 0.317e0
1457     .                *(1.      - max(zero,sign(unun,2.500e0 - rstar))))
1458          rstar1      = 0.      * max(zero,sign(unun,0.135e0 - rstar))
1459     .                +(1.      - max(zero,sign(unun,0.135e0 - rstar)))
1460     .                *(-0.55e0 * max(zero,sign(unun,2.500e0 - rstar))
1461     .                - 0.565
1462     .                *(1.      - max(zero,sign(unun,2.500e0 - rstar))))
1463          rstar2      = 0.      * max(zero,sign(unun,0.135e0 - rstar))
1464     .                +(1.      - max(zero,sign(unun,0.135e0 - rstar)))
1465     .                *(0.      * max(zero,sign(unun,2.500e0 - rstar))
1466     .                - 0.183
1467     .                *(unun    - max(zero,sign(unun,2.500e0 - rstar))))
1468 
1469         
1470
1471!XF    #RN (is_ok_z0h_rn) does not work well over bare ice
1472!XF    MAR is then too warm and not enough melt
1473 
1474         if(ro__SV(ikl,isnoSV(ikl))>50
1475     .  .and.ro__SV(ikl,isnoSV(ikl))<roSdSV)then
1476 
1477             Z0hnSV(ikl) = max(zero
1478     .                , sign(unun,zzsnsv(ikl,isnoSV(ikl))-epsi))
1479     .                * exp(rstar0+rstar1*alors+rstar2*alors*alors)
1480     .                * 0.001e0 + Z0hnSV(ikl) * ( 1. - max(zero
1481     .                , sign(unun,zzsnsv(ikl,isnoSV(ikl))-epsi)))
1482 
1483          endif
1484
1485
1486          ENDIF
1487 
1488          Z0h_SV(ikl) =     Z0hnSV(ikl)
1489 
1490
1491c #MT     Z0m_SV(ikl) = max(2.0e-6     ,Z0m_SV(ikl)) ! Min Z0_m (Garrat Scheme)
1492!          Z0m_SV(ikl) = min(Z0m_SV(ikl),za__SV(ikl)*0.3333)
1493     
1494
1495       END DO
1496 
1497
1498       return
1499       end
1500
1501
1502
1503
1504
1505
1506
1507
1508
1509
1510 
Note: See TracBrowser for help on using the repository browser.