Context Navigation

inlandsis.f90 @ 5151

Last change on this file since 5151 was 5117, checked in by abarral, 4 months ago
rename modules properly lmdz_* move some unused files to obsolete/ (lint) uppercase fortran keywords
File size: 59.5 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 : ut/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 !
226	! #sw real PorVol,rWater !
227	! #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
237	! #ZN real sqrrZ0 ! ut/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
245	! #SZ real Z0Sa_N ! Regime Snow Roughness Length
246	! #SZ real Z0SaSi ! 1.IF Rgm Snow Roughness Length
247	! #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(LEN=3) :: qsalt_param ! Switch for saltation flux param.
308	CHARACTER(LEN=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)
344	! #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
352	! + ! for Suspension
353	! + ! Pommeroy, Gray and Landine 1993,
354	! + ! J. Hydrology, 144(8) p.169
355	data nit / 5 / ! us(is0,uth) recursivity: Nb Iterations
356	!c#AE data qsalt_param/"bin"/ ! saltation part. conc. from Bintanja 2001 (p
357	data qsalt_param/"pom"/ ! saltation part. conc. from Pomeroy and Gray
358	!c#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
399	! +--Aeolian erosion and Blowing Snow
400	! +==================================
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==1) THEN
412	IF(isnoSV(1)>=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
416	! + **********
417	CALL SISVAT_BSn
418	endif
419	else
420	CALL SISVAT_BSn
421	! + **********
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
434	! +--Computation of threshold friction velocity for snow erosion
435	! ---------------------------------------------------------------
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
442	! +--Snow Properties
443	! + ~~~~~~~~~~~~~~~
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
469	! +--Mobility Index (Guyomarc'h & Merindol 1997, Ann.Glaciol.)
470	! + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
471	SaltM1 = -0.750e-2 * G1snSV(ikl,isn) &
472	-0.500e-2 * G2snSV(ikl,isn)+ 0.500e00 !dendritic case
473	! + CAUTION: Guyomarc'h & Merindol Dendricity Sign is +
474	! + ^^^^^^^^ 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
478	! 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)
482	!c#AE FacRho = 1.25 - 0.0042 * ro__SV(ikl,isn)
483	!c#AE SaltMo = 0.34 * SaltMo + 0.66 * FacRho !needed for polar snow
484	MIN_Mo = 0.
485	! SaltMo = max(SaltMo,MIN_Mo)
486	! SaltMo = SaltOK * SaltMo + (1.-SaltOK) * min(SaltMo,SaltMx)
487	! #TUNE SaltMo = SaltOK * SaltMo - (1.-SaltOK) * 0.9500
488	SaltMo = max(SaltMo,epsi-unun)
489
490	! +--Influence of Density on Threshold Shear Stress
491	! + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
492	Por_BS = 1. - 300. / ro_Ice
493	ShearS = Por_BS / (1.-Por_BS)
494	! +... SheaBS = Arg(sqrt(shear = max shear stress in snow)):
495	! + shear = 3.420d00 * exp(-(Por_BS +Por_BS)
496	! + . /(unun -Por_BS))
497	! + SheaBS : see de Montmollin (1978),
498	! + These Univ. Sci. Medic. Grenoble, Fig. 1 p. 124
499
500	! +--Snow Drift Index (Guyomarc'h & Merindol 1997, Ann.Glaciol.)
501	! + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
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
508	! +--Threshold Friction Velocity
509	! + ~~~~~~~~~~~~~~~~~~~~~~~~~~~
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)
518	! + Gallee et al., 2001 eq 5, p5
519
520	IF (usth_param == "gal") THEN
521	Salt_us = (log(2.868) - log(1 + SaltMo)) * rCd10n/0.085
522	Salt_us = Salt_us * Fac_Mo
523	! +... Salt_us : Extension of Guyomarc'h & Merindol 1998 with
524	! +... de Montmollin (1978). Gallee et al. 2001
525	endif
526
527	IF (usth_param == "lis") then !Liston et al. 2007
528	IF(ro__SV(ikl,isn)>300.) THEN
529	Salt_us = 0.005exp(0.013ro__SV(ikl,isn))
530	else
531	Salt_us = 0.01exp(0.003ro__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_4275.27sss__F /sss__N
575	! . - us_327*4.27
576	! . - us_1272.27qsnoSV(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 == "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 == "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
634	! + ------------------------------------------------------
635	! +--Buffer Layer
636	! + -----------------------------------------------------
637
638	DO ikl=1,knonv
639	! 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
645	! +--Snow Density
646	! + ^^^^^^^^^^^^
647	Polair = zero
648	! #NP Polair = max(zero, !
649	! #NP. sign(unun,TaPole !
650	! #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
657	! #NP BufPro = max( rosMin, ! Fallen Snow Density
658	! #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.
662	! #BS density_kotlyakov = .FALSE. !C.Amory BS 2018
663	! + ... 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
668	! + ... [ A compromise between
669	! + ... Kotlyakov (1961) and Lenaerts (2012, JGR, Part1) ]
670	IF (tt_tmp>=-10) THEN
671	BufPro = max( rosMin, &
672	104. *sqrt( max( vv_tmp-6.0,0.0))) ! Kotlyakov (1961)
673	else
674	vv_virt = (tt_cvv_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
679	! + ... [ density derived from observations of the first 50cm of
680	! + ... snow - cf. Rajashree Datta - and multiplied by 0.8 ]
681	! + ... C. Agosta, 2016-09
682	!c #SD BufPro = 149.2 + 6.84VV10SV(ikl) + 0.48Tsrfsv(ikl)
683	!c #SD BufPro = 125 + 14VV10SV(ikl) + 0.6Tsrfsv(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
719	! +-- S.Falling Snow Properties (computed as in SISVAT_zAg)
720	! + ^^^^^^^^^^^^^^^^^^^^^^^
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 = (zroNEWSizNEW+zroOLDSizOLD) ! Averaged Size
772	Sph_av = min( zroNEWSphNEW+zroOLDSphOLD & !
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 )) !
783	! +... REMARQUE: le type moyen (dendritique ou non) depend
784	! + ^^^^^^^^ de la comparaison avec le diametre optique
785	! + d'une neige recente de dendricite nulle
786	! +... REMARK: the mean type (dendritic or not) depends
787	! + ^^^^^^ on the comparaison with the optical diameter
788	! + 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_BufsBuf_G1+BufsSV(ikl)BG1sSV(ikl)) &
808	/max(epsi,Bufs_N)
809	G2same = (d_BufsBuf_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+zroOLDSizOLD) ! Averaged Size
838	Sph_av = min( zroNEW SphNEW+zroOLDSphOLD & !
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 )) !
849	! +... REMARQUE: le type moyen (dendritique ou non) depend
850	! + ^^^^^^^^ de la comparaison avec le diametre optique
851	! + d'une neige recente de dendricite nulle
852	! +... REMARK: the mean type (dendritic or not) depends
853	! + ^^^^^^ on the comparaison with the optical diameter
854	! + 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
871	! +--Update of Buffer Layer Content & Decision about creating a new snow layer
872	! + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
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	-SMndSV3.00))) ! is Buffer > SMndSV3
884	Bdzssv(ikl) = 1.e-3BufsSV(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==1) THEN
897	IF(isnoSV(1)>=1.OR.NLaysv(1)>=1) THEN
898	! + **********
899	CALL SISVAT_zSn
900	! + **********
901	endif
902	else
903	! + **********
904	CALL SISVAT_zSn
905	! + **********
906	ENDIF
907
908	! + **********
909	! #ve CALL SISVAT_wEq('_zSn ',0)
910	! + **********
911
912	! Add a new Snow Layer
913	! ====================
914
915	DO ikl=1,knonv
916
917	isnoSV(ikl) = isnoSV(ikl) +NLaysv(ikl)
918	isn = isnoSV(ikl)
919	dzsnSV(ikl,isn) = dzsnSV(ikl,isn) * (1-NLaysv(ikl)) &
920	+ Bdzssv(ikl) * NLaysv(ikl)
921	TsisSV(ikl,isn) = TsisSV(ikl,isn) * (1-NLaysv(ikl)) &
922	+ min(TaT_SV(ikl),Tf_Sno) *NLaysv(ikl)
923	ro__SV(ikl,isn) = ro__SV(ikl,isn) * (1-NLaysv(ikl)) &
924	+ Brossv(ikl) * NLaysv(ikl)
925	eta_SV(ikl,isn) = eta_SV(ikl,isn) * (1-NLaysv(ikl)) ! + 0.
926	agsnSV(ikl,isn) = agsnSV(ikl,isn) * (1-NLaysv(ikl)) ! + 0.
927	G1snSV(ikl,isn) = G1snSV(ikl,isn) * (1-NLaysv(ikl)) &
928	+ BG1ssv(ikl) * NLaysv(ikl)
929	G2snSV(ikl,isn) = G2snSV(ikl,isn) * (1-NLaysv(ikl)) &
930	+ BG2ssv(ikl) * NLaysv(ikl)
931	istoSV(ikl,isn) = istoSV(ikl,isn) * (1-NLaysv(ikl)) &
932	+ max(zer0,sign(un_1,TaT_SV(ikl) &
933	-Tf_Sno-eps_21)) * istdSV(2) &
934	* NLaysv(ikl)
935	BufsSV(ikl) = BufsSV(ikl) * (1-NLaysv(ikl))
936	NLaysv(ikl) = 0
937
938
939	END DO
940
941
942	! Snow Pack Thickness
943	! -------------------
944
945	DO ikl=1,knonv
946	z_snsv(ikl) = 0.0
947	END DO
948	DO isn=1,nsno
949	DO ikl=1,knonv
950	z_snsv(ikl) = z_snsv(ikl) + dzsnSV(ikl,isn)
951	zzsnsv(ikl,isn) = z_snsv(ikl)
952	END DO
953	END DO
954
955
956
957	END IF ! SnoMod
958
959
960
961	! Soil Albedo: Soil Humidity Correction
962	! ==========================================
963
964	! REFERENCE: McCumber and Pielke (1981), Pielke (1984)
965	! ^^^^^^^^^
966	DO ikl=1,knonv
967	albssv(ikl) = &
968	alb0SV(ikl) *(1.0-min(half,eta_SV( ikl,0) &
969	/etadSV(isotSV(ikl))))
970	! REMARK: Albedo of Water Surfaces (isotSV=0):
971	! ^^^^^^ alb0SV := 2 X effective value, while
972	! eta_SV := etadSV
973	END DO
974
975
976	! Snow Pack Optical Properties
977	! ============================
978
979	IF (SnoMod) THEN
980
981	! ******
982	CALL SnOptP(jjtime)
983	! ******
984
985	ELSE
986	DO ikl=1,knonv
987	sEX_sv(ikl,1) = 1.0
988	sEX_sv(ikl,0) = 0.0
989	albisv(ikl) = albssv(ikl)
990	END DO
991	END IF
992
993
994
995	! Soil optical properties
996	! =============================
997	!Etienne: as in inlandis we do not CALL vgopt, we need to define
998	!the albedo alb_SV and to calculate the
999	!absorbed Solar Radiation by Surfac (Normaliz)[-] SoSosv
1000
1001
1002	DO ikl=1,klonv
1003
1004	e_pRad = 2.5 * coszSV(ikl) ! exponential argument,
1005	! V/nIR radiation partitioning,
1006	! DR97, 2, eqn (2.53) & (2.54)
1007	e1pRad = 1.-exp(-e_pRad) ! exponential, V/nIR Rad. Part.
1008	exdRad= 1.
1009
1010	! Visible Part of the Solar Radiation Spectrum (V, 0.4--0.7mi.m)
1011	A_Rad0 = 0.25 + 0.697 * e1pRad ! Absorbed Vis. Radiation
1012	absg_V = (1.-albisv(ikl))(A_Rad0exdRad) !
1013
1014	! Near-IR Part of the Solar Radiation Spectrum (nIR, 0.7--2.8mi.m)
1015
1016	A_Rad0 = 0.80 + 0.185 * e1pRad ! Absorbed nIR. Radiation
1017	absgnI = (1.-albisv(ikl))(A_Rad0exdRad) !
1018
1019	SoSosv(ikl) = (absg_V+absgnI)*0.5d0
1020
1021	alb_SV(ikl) = albisv(ikl)
1022
1023	END DO
1024
1025	! **********
1026	! #ve CALL SISVAT_wEq('SnOptP',0)
1027	! **********
1028
1029
1030	! Surface Emissivity (Etienne: simplified calculation for landice)
1031	! =============================================================
1032
1033	DO ikl=1,knonv
1034	LSnMsk = min( 1,isnoSV(ikl))
1035	Eso_sv(ikl)= EmiSol(1-LSnMsk)+EmiSnoLSnMsk ! Sol+Sno Emissivity
1036	emi_SV(ikl)= EmiSol(1-LSnMsk) + EmiSnoLSnMsk
1037	END DO
1038
1039
1040
1041
1042	! Upward IR (INPUT, from previous time step)
1043	! ===================================================================
1044
1045	DO ikl=1,knonv
1046	! #e1 Enrsvd(ikl) = - IRs_SV(ikl)
1047	IRupsv(ikl) = IRs_SV(ikl)
1048	END DO
1049
1050
1051	! Turbulence
1052	! ==========
1053
1054	! Latent Heat of Vaporization/Sublimation
1055	! ---------------------------------------
1056
1057	DO ikl=1,knonv
1058	SnoWat = min(isnoSV(ikl),0)
1059	Lx_H2O(ikl) = &
1060	(1.-SnoWat) * LhvH2O &
1061	+ SnoWat (LhsH2O (1.-eta_SV(ikl,isnoSV(ikl))) &
1062	+LhvH2O * eta_SV(ikl,isnoSV(ikl)) )
1063	END DO
1064
1065
1066
1067
1068	! Aerodynamic Resistance (calculated from drags given by LMDZ)
1069	! Commented because already calculated in surf_inlandsis_mod
1070	! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1071	! DO ikl=1,knonv
1072	! ram_sv(ikl) = 1./(cdM_SV(ikl)*max(VV__SV(ikl),eps6))
1073	! rah_sv(ikl) = 1./(cdH_SV(ikl)*max(VV__SV(ikl),eps6))
1074	! END DO
1075
1076
1077
1078	! Soil Energy Balance
1079	! =====================
1080
1081
1082	IF (iflag_temp_inlandsis == 0) THEN
1083	CALL SISVAT_TSo
1084
1085	else
1086	DO ikl=1,knonv
1087	Tsf_SV(ikl)=Tsrfsv(ikl)
1088	END DO
1089
1090	CALL SISVAT_TS2
1091
1092	end if
1093
1094
1095	! **********
1096	! #ve CALL SISVAT_wEq('_TSo ',0)
1097	! **********
1098
1099
1100
1101	! Soil Water Potential
1102	! ------------------------
1103
1104	DO isl=-nsol,0
1105	DO ikl=1,knonv
1106	ist = isotSV(ikl) ! Soil Type
1107	psi_sv(ikl,isl) = psidSV(ist) & ! DR97, Eqn.(3.34)
1108	*(etadSV(ist) /max(eps6,eta_SV(ikl,isl))) & !
1109	**bCHdSV(ist) !
1110
1111
1112	! Soil Hydraulic Conductivity
1113	! ---------------------------
1114
1115	Khydsv(ikl,isl) = s2__SV(ist) & ! DR97, Eqn.(3.35)
1116	(eta_SV(ikl,isl)(2.bCHdSV(ist)+3.)) !
1117	END DO
1118	END DO
1119
1120
1121	! Melting / Refreezing in the Snow Pack
1122	! =====================================
1123
1124	IF (SnoMod) THEN
1125
1126	! **********
1127	CALL SISVAT_qSn
1128	! **********
1129
1130	! **********
1131	! #ve CALL SISVAT_wEq('_qSn ',0)
1132	! **********
1133
1134
1135	! Snow Pack Thickness
1136	! -------------------
1137
1138	DO ikl=1,knonv
1139	z_snsv(ikl) = 0.0
1140	END DO
1141	DO isn=1,nsno
1142	DO ikl=1,knonv
1143	z_snsv(ikl) = z_snsv(ikl) + dzsnSV(ikl,isn)
1144	zzsnsv(ikl,isn) = z_snsv(ikl)
1145	END DO
1146	END DO
1147
1148
1149	! Energy in Excess is added to the first Soil Layer
1150	! -------------------------------------------------
1151	DO ikl=1,knonv
1152	z_snsv(ikl) = max(zer0, &
1153	sign(un_1,eps6-z_snsv(ikl)))
1154	TsisSV(ikl,0) = TsisSV(ikl,0) + EExcsv(ikl) &
1155	/(rocsSV(isotSV(ikl)) &
1156	+rcwdSV*eta_SV(ikl,0))
1157	EExcsv(ikl) = 0.
1158	END DO
1159
1160
1161	END IF
1162
1163
1164	! Soil Water Balance
1165	! =====================
1166
1167	! **********
1168	CALL SISVAT_qSo
1169	! #m0. (Wats_0,Wats_1,Wats_d)
1170	! **********
1171
1172
1173	! Surface Fluxes
1174	! =====================
1175
1176	DO ikl=1,knonv
1177	IRdwsv(ikl)=IRd_SV(ikl)*Eso_sv(ikl) ! Downward IR
1178	! IRdwsv(ikl)=tau_sv(ikl) IRd_SV(ikl)Eso_sv(ikl) ! Downward IR
1179	! . +(1.0-tau_sv(ikl))IRd_SV(ikl)Evg_sv(ikl) ! ! Etienne, remove vegetation component
1180	IRupsv(ikl) = IRupsv(ikl) ! Upward IR
1181	IRu_SV(ikl) = -IRupsv(ikl) & ! Upward IR
1182	+IRd_SV(ikl) & ! (effective)
1183	-IRdwsv(ikl) ! (positive)
1184
1185	TBr_sv(ikl) =sqrt(sqrt(IRu_SV(ikl)/StefBo)) ! Brightness
1186	! Temperature
1187	uts_SV(ikl) = (HSv_sv(ikl) +HSs_sv(ikl)) & ! uT
1188	/(rhT_SV(ikl) *cp) !
1189	uqs_SV(ikl) = (HLv_sv(ikl) +HLs_sv(ikl)) & ! uq
1190	/(rhT_SV(ikl) *LhvH2O) !
1191	LMO_SV(ikl) = TaT_SV(ikl)(us__SV(ikl)*3) &
1192	/gravit/uts_SV(ikl)/vonKrm ! MO length
1193
1194
1195	! Surface Temperature
1196	! ^^^^^^^^^^^^^^^^^^^^
1197
1198	IF (iflag_tsurf_inlandsis == 0) THEN
1199
1200	Tsrfsv(ikl) =TsisSV(ikl,isnoSV(ikl))
1201
1202	ELSE IF (iflag_tsurf_inlandsis > 0) THEN
1203	! Etienne: extrapolation from the two uppermost levels:
1204
1205	IF (isnoSV(ikl) >=2) THEN
1206	zm1=-dzsnSV(ikl,isnoSV(ikl))/2.
1207	zm2=-(dzsnSV(ikl,isnoSV(ikl)) + dzsnSV(ikl,isnoSV(ikl)-1)/2.)
1208	ELSE IF (isnoSV(ikl) == 1) THEN
1209	zm1=-dzsnSV(ikl,isnoSV(ikl))/2.
1210	zm2=-(dzsnSV(ikl,isnoSV(ikl))+dz_dSV(0)/2.)
1211	else
1212	zm1=-dz_dSV(0)/2.
1213	zm2=-(dz_dSV(0)+dz_dSV(-1)/2.)
1214
1215	end if
1216
1217	coefslope=(TsisSV(ikl,isnoSV(ikl))-TsisSV(ikl,isnoSV(ikl)-1)) &
1218	/(zm1-zm2)
1219	Tsrfsv(ikl)=TsisSV(ikl,isnoSV(ikl))+coefslope*(0. - zm1)
1220
1221
1222	ELSE !(default)
1223
1224	Tsrfsv(ikl) =TsisSV(ikl,isnoSV(ikl))
1225
1226	END IF
1227
1228
1229	END DO
1230
1231	! Snow Pack Properties (sphericity, dendricity, size)
1232	! ===================================================
1233
1234	IF (SnoMod) THEN
1235
1236	IF (discret_xf .AND. klonv==1) THEN
1237	IF(isnoSV(1)>=1) THEN
1238	! + **********
1239	CALL SISVAT_GSn
1240	! + **********
1241	ENDIF
1242	else
1243	! + **********
1244	CALL SISVAT_GSn
1245	! + **********
1246	ENDIF
1247
1248
1249	END IF
1250
1251
1252	! Roughness Length for next time step
1253	!====================================
1254
1255	! Note that in INLANDSIS, we treat only ice covered surfaces so calculation
1256	! of z0 is much simpler (no subgrid fraction of ocean or land)
1257	! old calculations are commented below
1258
1259
1260	! +--Roughness Length for Momentum
1261	! + -----------------------------
1262
1263	! ETIENNE WARNING: changes have been made wrt original SISVAT
1264
1265	! +--Land+Sea-Ice / Ice-free Sea Mask
1266	! + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1267	DO ikl=1,knonv
1268	IcIndx(ikl) = 0
1269	ENDDO
1270	DO isn=1,nsno
1271	DO ikl=1,knonv
1272
1273	IcIndx(ikl) = max(IcIndx(ikl), &
1274	isn*max(0, &
1275	sign(1, &
1276	int(ro__SV(ikl,isn)-900.))))
1277	ENDDO
1278	ENDDO
1279
1280	DO ikl=1,knonv
1281	LISmsk = 1. ! in inlandsis, land only
1282	IceMsk = max(0,sign(1 ,IcIndx(ikl)-1) )
1283	SnoMsk = max(min(isnoSV(ikl)-iiceSV(ikl),1),0)
1284
1285
1286	! +--Z0 Smooth Regime over Snow (Andreas 1995, CRREL Report 95-16, p. 8)
1287	! + ^^^^^^^^^^^^^^^^^^^^^^^^^^
1288	Z0m_nu = 5.e-5 ! z0s~(10-d)*exp(-vonkar/sqrt(1.1e-03))
1289
1290	! +--Z0 Saltat.Regime over Snow (Gallee et al., 2001, BLM 99 (19) p.11)
1291	! + ^^^^^^^^^^^^^^^^^^^^^^^^^^
1292
1293	u2star = us__SV(ikl) *us__SV(ikl)
1294	Z0mBSn = u2star *0.536e-3 - 61.8e-6
1295	Z0mBSn = max(Z0mBS0 ,Z0mBSn)
1296
1297	! +--Z0 Smooth + Saltat. Regime
1298	! + ^^^^^^^^^^^^^^^^^^^^^^^^^^
1299	Z0enSV(ikl) = Z0m_nu &
1300	+ Z0mBSn
1301
1302
1303	! Calculation of snow roughness length
1304	!=====================================
1305	IF (iflag_z0m_snow == 0) THEN
1306
1307	Z0m_Sn=prescribed_z0m_snow
1308
1309	ELSE IF (iflag_z0m_snow == 1) THEN
1310
1311	Z0m_Sn=Z0enSV(ikl)
1312
1313	ELSE IF (iflag_z0m_snow == 2) THEN
1314
1315	! +--Rough Snow Surface Roughness Length (Variable Sastrugi Height)
1316	! + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1317	A_Fact = 1.0000 ! Andreas et al., 2004, p.4
1318	! ams.confex.com/ams/pdfpapers/68601.pdf
1319
1320	! Parameterization of z0 dependance on Temperature (C. Amory, 2017)
1321	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1322	! Z0=f(T) deduced from observations, Adelie Land, dec2012-dec2013
1323
1324
1325	coefa = 0.1658 !0.1862 !Ant
1326	coefb = -50.3869 !-55.7718 !Ant
1327	ta1 = 253.15 !255. Ant
1328	ta2 = 273.15
1329	ta3 = 273.15+3
1330	z01 = exp(coefa*ta1 + coefb) !~0.2 ! ~0.25 mm
1331	z02 = exp(coefa*ta2 + coefb) !~6 !~7 mm
1332	z03 = z01
1333	coefc = log(z03/z02)/(ta3-ta2)
1334	coefd = log(z03)-coefc*ta3
1335
1336	IF (TaT_SV(ikl) < ta1) THEN
1337	Z0_obs = z01
1338	ELSE IF (TaT_SV(ikl)>=ta1 .AND. TaT_SV(ikl)<ta2) THEN
1339	Z0_obs = exp(coefa*TaT_SV(ikl) + coefb)
1340	ELSE IF (TaT_SV(ikl)>=ta2 .AND. TaT_SV(ikl)<ta3) THEN
1341	! if st > 0, melting induce smooth surface
1342	Z0_obs = exp(coefc*TaT_SV(ikl) + coefd)
1343	else
1344	Z0_obs = z03
1345	endif
1346
1347	Z0m_Sn=Z0_obs
1348
1349
1350	ELSE
1351
1352	Z0m_Sn=0.500e-3 ! default=0.500e-3m (tuning of MAR)
1353
1354	ENDIF
1355
1356
1357
1358	! param = Z0_obs/1. ! param(s) \| 1.(m/s)=TUNING
1359	! #SZ Z0Sa_N = (us__SV(ikl) -0.2)*param ! 0.0001=TUNING
1360	! #SZ. * max(zero,sign(unun,TfSnow-eps9
1361	! #SZ. -TsisSV(ikl , isnoSV(ikl))))
1362	!!#SZ Z0SaSi = max(zero,sign(unun,Z0Sa_N ))! 1 if erosion
1363	! #SZ Z0SaSi = max(zero,sign(unun,zero -eps9 -uss_SV(ikl)))!
1364	! #SZ Z0Sa_N = max(zero, Z0Sa_N)
1365	! #SZ Z0SaSV(ikl) =
1366	! #SZ. max(Z0SaSV(ikl) ,Z0SaSV(ikl)
1367	! #SZ. + Z0SaSi(Z0Sa_N-Z0SaSV(ikl))exp(-dt__SV/43200.))
1368	! #SZ. - min(dz0_SV(ikl) , Z0SaSV(ikl))
1369
1370	! #SZ A_Fact = Z0SaSV(ikl) * 5.0/0.15 ! A=5 if h~10cm
1371	! +... CAUTION: The influence of the sastrugi direction is not yet included
1372
1373	! #SZ Z0m_Sn = Z0SaSV(ikl) !
1374	! #SZ. - Z0m_nu !
1375
1376	! +--Z0 Saltat.Regime over Snow (Shao & Lin, 1999, BLM 91 (46) p.222)
1377	! + ^^^^^^^^^^^^^^^^^^^^^^^^^^
1378	! #ZN sqrrZ0 = usthSV(ikl)/max( us__SV(ikl),0.001)
1379	! #ZN sqrrZ0 = min( sqrrZ0 ,0.999)
1380	! #ZN Z0mBSn = 0.55 0.55 exp(-sqrrZ0 *sqrrZ0)
1381	! #ZN. us__SV(ikl) us__SV(ikl)grvinv0.5
1382
1383	! +--Z0 Smooth + Saltat. Regime (Shao & Lin, 1999, BLM 91 (46) p.222)
1384	! + ^^^^^^^^^^^^^^^^^^^^^^^^^^
1385	! #ZN Z0enSV(ikl) = (Z0m_nu ** sqrrZ0 )
1386	! #ZN. * (Z0mBSn **(1.-sqrrZ0))
1387	! #ZN Z0enSV(ikl) = max(Z0enSV(ikl), Z0m_nu)
1388
1389
1390	! +--Z0 Smooth Regime over Snow (Andreas etAl., 2004
1391	! + ^^^^^^^^^^^^^^^^^^^^^^^^^^ ams.confex.com/ams/pdfpapers/68601.pdf)
1392	! #ZA Z0m_nu = 0.135*akmol / max(us__SV(ikl) , epsi)
1393
1394	! +--Z0 Saltat.Regime over Snow (Andreas etAl., 2004
1395	! + ^^^^^^^^^^^^^^^^^^^^^^^^^^ ams.confex.com/ams/pdfpapers/68601.pdf)
1396	! #ZA Z0mBSn = 0.035u2star grvinv
1397
1398	! +--Z0 Smooth + Saltat. Regime (Andreas etAl., 2004
1399	! ( used by Erosion) ams.confex.com/ams/pdfpapers/68601.pdf)
1400	! ^^^^^^^^^^^^^^^^^^^^^^^^^^
1401	! #ZA Z0enSV(ikl) = Z0m_nu
1402	! #ZA. + Z0mBSn
1403
1404	! +--Z0 Rough Regime over Snow (Andreas etAl., 2004
1405	! + (.NOT. used by Erosion) ams.confex.com/ams/pdfpapers/68601.pdf)
1406	! ^^^^^^^^^^^^^^^^^^^^^^^^^^
1407	!!#ZA u2star = (us__SV(ikl) -0.1800) / 0.1
1408	!!#ZA Z0m_Sn =A_FactZ0mBSn exp(-u2star*u2star)
1409	! #ZA Z0m_90 =(10.-0.025*VVs_SV(ikl)/5.)
1410	! #ZA. exp(-0.4/sqrt(.00275+.00001max(0.,VVs_SV(ikl)-5.)))
1411	! #ZA Z0m_Sn = DDs_SV(ikl)* Z0m_90 / 45.
1412	! #ZA. - DDs_SV(ikl)DDs_SV(ikl) Z0m_90 /(90.*90.)
1413
1414
1415
1416
1417	! +--Z0 (Erosion) over Snow (instantaneous)
1418	! + ^^^^^^^^^^^^^^^^^^^^^^^^^^
1419	Z0e_SV(ikl) = Z0enSV(ikl)
1420
1421	! +--Momentum Roughness Length (Etienne: changes wrt original SISVAT)
1422	! + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1423	Z0mnSV(ikl) = Z0m_nu *(1-SnoMsk) & ! Ice z0
1424	+ (Z0m_Sn)*SnoMsk ! Snow Sastrugi Form and Snow Erosion
1425
1426
1427	! +--GIS Roughness Length
1428	! + ^^^^^^^^^^^^^^^^^^^^^
1429	! #GL Z0mnSV(ikl) =
1430	! #GL. (1-LSmask(ikl)) * Z0mnSV(ikl)
1431	! #GL. + LSmask(ikl) * max(Z0mnSV(ikl),max(Z0_GIM,
1432	! #GL. Z0_GIM+
1433	! #GL. (0.0032-Z0_GIM)*(ro__SV(ikl,isnoSV(ikl))-600.) !
1434	! #GL. /(920.00 -600.))) !
1435
1436	! +--Mom. Roughness Length, Instantaneous
1437	! + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1438	Z0m_SV(ikl) = Z0mnSV(ikl) ! Z0mnSV instant.
1439
1440
1441	! +--Roughness Length for Scalars
1442	! + ----------------------------
1443
1444	Z0hnSV(ikl) = Z0mnSV(ikl)/ 7.4
1445
1446	IF (is_ok_z0h_rn) THEN
1447
1448	rstar = Z0mnSV(ikl) * us__SV(ikl) / akmol
1449	rstar = max(epsi,min(rstar,R_1000))
1450	alors = log(rstar)
1451	rstar0 = 1.250e0 * max(zero,sign(unun,0.135e0 - rstar)) &
1452	+(1. - max(zero,sign(unun,0.135e0 - rstar))) &
1453	(0.149e0 max(zero,sign(unun,2.500e0 - rstar)) &
1454	+ 0.317e0 &
1455	*(1. - max(zero,sign(unun,2.500e0 - rstar))))
1456	rstar1 = 0. * max(zero,sign(unun,0.135e0 - rstar)) &
1457	+(1. - max(zero,sign(unun,0.135e0 - rstar))) &
1458	(-0.55e0 max(zero,sign(unun,2.500e0 - rstar)) &
1459	- 0.565 &
1460	*(1. - max(zero,sign(unun,2.500e0 - rstar))))
1461	rstar2 = 0. * max(zero,sign(unun,0.135e0 - rstar)) &
1462	+(1. - max(zero,sign(unun,0.135e0 - rstar))) &
1463	(0. max(zero,sign(unun,2.500e0 - rstar)) &
1464	- 0.183 &
1465	*(unun - max(zero,sign(unun,2.500e0 - rstar))))
1466
1467
1468
1469	!XF #RN (is_ok_z0h_rn) does not work well over bare ice
1470	!XF MAR is then too warm and not enough melt
1471
1472	IF(ro__SV(ikl,isnoSV(ikl))>50 &
1473	.AND.ro__SV(ikl,isnoSV(ikl))<roSdSV)THEN
1474	Z0hnSV(ikl) = max(zero &
1475	, sign(unun,zzsnsv(ikl,isnoSV(ikl))-epsi)) &
1476	* exp(rstar0+rstar1alors+rstar2alors*alors) &
1477	* 0.001e0 + Z0hnSV(ikl) * ( 1. - max(zero &
1478	, sign(unun,zzsnsv(ikl,isnoSV(ikl))-epsi)))
1479
1480	endif
1481
1482
1483	ENDIF
1484
1485	Z0h_SV(ikl) = Z0hnSV(ikl)
1486
1487
1488	! #MT Z0m_SV(ikl) = max(2.0e-6 ,Z0m_SV(ikl)) ! Min Z0_m (Garrat Scheme)
1489	! Z0m_SV(ikl) = min(Z0m_SV(ikl),za__SV(ikl)*0.3333)
1490
1491
1492	END DO
1493
1494
1495
1496	END SUBROUTINE inlandsis
1497
1498
1499
1500
1501
1502
1503
1504
1505
1506
1507

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format