Context Navigation

inlandsis.F @ 5185

Last change on this file since 5185 was 4013, checked in by lguez, 3 years ago
Sync latest trunk changes to Ocean_skin
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 : 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	c #sw real PorVol,rWater !
227	c #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	c #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	c #SZ real Z0Sa_N ! Regime Snow Roughness Length
246	c #SZ real Z0SaSi ! 1.IF Rgm Snow Roughness Length
247	c #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	character3 usth_param ! Switch for ut 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	c #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	C + ! for Suspension
353	C + ! Pommeroy, Gray and Landine 1993,
354	C + ! J. Hydrology, 144(8) p.169
355	data nit / 5 / ! us(is0,uth) recursivity: Nb Iterations
356	cc#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	cc#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	C +--Aeolian erosion and Blowing Snow
400	C +==================================
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
416	C + **********
417	call SISVAT_BSn
418	endif
419	else
420	call SISVAT_BSn
421	C + **********
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	C +--Computation of threshold friction velocity for snow erosion
435	C ---------------------------------------------------------------
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	C +--Snow Properties
443	C + ~~~~~~~~~~~~~~~
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	C +--Mobility Index (Guyomarc'h & Merindol 1997, Ann.Glaciol.)
470	C + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
471	SaltM1 = -0.750e-2 * G1snSV(ikl,isn)
472	. -0.500e-2 * G2snSV(ikl,isn)+ 0.500e00 !dendritic case
473	C + CAUTION: Guyomarc'h & Merindol Dendricity Sign is +
474	C + ^^^^^^^^ 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	c 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	cc#AE FacRho = 1.25 - 0.0042 * ro__SV(ikl,isn)
483	cc#AE SaltMo = 0.34 * SaltMo + 0.66 * FacRho !needed for polar snow
484	MIN_Mo = 0.
485	c SaltMo = max(SaltMo,MIN_Mo)
486	c SaltMo = SaltOK * SaltMo + (1.-SaltOK) * min(SaltMo,SaltMx)
487	c #TUNE SaltMo = SaltOK * SaltMo - (1.-SaltOK) * 0.9500
488	SaltMo = max(SaltMo,epsi-unun)
489
490	C +--Influence of Density on Threshold Shear Stress
491	C + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
492	Por_BS = 1. - 300. / ro_Ice
493	ShearS = Por_BS / (1.-Por_BS)
494	C +... SheaBS = Arg(sqrt(shear = max shear stress in snow)):
495	C + shear = 3.420d00 * exp(-(Por_BS +Por_BS)
496	C + . /(unun -Por_BS))
497	C + SheaBS : see de Montmollin (1978),
498	C + These Univ. Sci. Medic. Grenoble, Fig. 1 p. 124
499
500	C +--Snow Drift Index (Guyomarc'h & Merindol 1997, Ann.Glaciol.)
501	C + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
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	C +--Threshold Friction Velocity
509	C + ~~~~~~~~~~~~~~~~~~~~~~~~~~~
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	C + 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
523	C +... Salt_us : Extension of Guyomarc'h & Merindol 1998 with
524	C +... 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.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 .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
634	C + ------------------------------------------------------
635	C +--Buffer Layer
636	C + -----------------------------------------------------
637
638	DO ikl=1,knonv
639	c 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	C +--Snow Density
646	C + ^^^^^^^^^^^^
647	Polair = zero
648	c #NP Polair = max(zero, !
649	c #NP. sign(unun,TaPole !
650	c #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	c #NP BufPro = max( rosMin, ! Fallen Snow Density
658	c #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	c #BS density_kotlyakov = .false. !C.Amory BS 2018
663	C + ... 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	C + ... [ A compromise between
669	C + ... 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_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	C + ... [ density derived from observations of the first 50cm of
680	C + ... snow - cf. Rajashree Datta - and multiplied by 0.8 ]
681	C + ... C. Agosta, 2016-09
682	cc #SD BufPro = 149.2 + 6.84VV10SV(ikl) + 0.48Tsrfsv(ikl)
683	cc #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	C +-- S.Falling Snow Properties (computed as in SISVAT_zAg)
720	C + ^^^^^^^^^^^^^^^^^^^^^^^
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	C +... REMARQUE: le type moyen (dendritique ou non) depend
784	C + ^^^^^^^^ de la comparaison avec le diametre optique
785	C + d'une neige recente de dendricite nulle
786	C +... REMARK: the mean type (dendritic or not) depends
787	C + ^^^^^^ on the comparaison with the optical diameter
788	C + 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	C +... REMARQUE: le type moyen (dendritique ou non) depend
850	C + ^^^^^^^^ de la comparaison avec le diametre optique
851	C + d'une neige recente de dendricite nulle
852	C +... REMARK: the mean type (dendritic or not) depends
853	C + ^^^^^^ on the comparaison with the optical diameter
854	C + 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	C +--Update of Buffer Layer Content & Decision about creating a new snow layer
872	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
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.eq.1) then
897
898	if(isnoSV(1).ge.1.or.NLaysv(1).ge.1) then
899	C + **********
900	call SISVAT_zSn
901	C + **********
902	endif
903	else
904	C + **********
905	call SISVAT_zSn
906	C + **********
907	endif
908
909	C + **********
910	! #ve call SISVAT_wEq('_zSn ',0)
911	C + **********
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_Rad0exdRad) !
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_Rad0exdRad) !
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)+EmiSnoLSnMsk ! Sol+Sno Emissivity
1037	emi_SV(ikl)= EmiSol(1-LSnMsk) + EmiSnoLSnMsk
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)) ! uT
1190	. /(rhT_SV(ikl) *cp) !
1191	uqs_SV(ikl) = (HLv_sv(ikl) +HLs_sv(ikl)) ! uq
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
1240	C + **********
1241	call SISVAT_GSn
1242	C + **********
1243	endif
1244	else
1245	C + **********
1246	call SISVAT_GSn
1247	C + **********
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
1262	C +--Roughness Length for Momentum
1263	C + -----------------------------
1264
1265	! ETIENNE WARNING: changes have been made wrt original SISVAT
1266
1267	C +--Land+Sea-Ice / Ice-free Sea Mask
1268	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
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
1288	C +--Z0 Smooth Regime over Snow (Andreas 1995, CRREL Report 95-16, p. 8)
1289	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^
1290	Z0m_nu = 5.e-5 ! z0s~(10-d)*exp(-vonkar/sqrt(1.1e-03))
1291
1292	C +--Z0 Saltat.Regime over Snow (Gallee et al., 2001, BLM 99 (19) p.11)
1293	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^
1294
1295	u2star = us__SV(ikl) *us__SV(ikl)
1296	Z0mBSn = u2star *0.536e-3 - 61.8e-6
1297	Z0mBSn = max(Z0mBS0 ,Z0mBSn)
1298
1299	C +--Z0 Smooth + Saltat. Regime
1300	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^
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
1317	C +--Rough Snow Surface Roughness Length (Variable Sastrugi Height)
1318	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
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
1361	c #SZ Z0Sa_N = (us__SV(ikl) -0.2)*param ! 0.0001=TUNING
1362	c #SZ. * max(zero,sign(unun,TfSnow-eps9
1363	c #SZ. -TsisSV(ikl , isnoSV(ikl))))
1364	!!#SZ Z0SaSi = max(zero,sign(unun,Z0Sa_N ))! 1 if erosion
1365	c #SZ Z0SaSi = max(zero,sign(unun,zero -eps9 -uss_SV(ikl)))!
1366	c #SZ Z0Sa_N = max(zero, Z0Sa_N)
1367	c #SZ Z0SaSV(ikl) =
1368	c #SZ. max(Z0SaSV(ikl) ,Z0SaSV(ikl)
1369	c #SZ. + Z0SaSi(Z0Sa_N-Z0SaSV(ikl))exp(-dt__SV/43200.))
1370	c #SZ. - min(dz0_SV(ikl) , Z0SaSV(ikl))
1371
1372	c #SZ A_Fact = Z0SaSV(ikl) * 5.0/0.15 ! A=5 if h~10cm
1373	C +... CAUTION: The influence of the sastrugi direction is not yet included
1374
1375	c #SZ Z0m_Sn = Z0SaSV(ikl) !
1376	c #SZ. - Z0m_nu !
1377
1378	C +--Z0 Saltat.Regime over Snow (Shao & Lin, 1999, BLM 91 (46) p.222)
1379	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^
1380	c #ZN sqrrZ0 = usthSV(ikl)/max( us__SV(ikl),0.001)
1381	c #ZN sqrrZ0 = min( sqrrZ0 ,0.999)
1382	c #ZN Z0mBSn = 0.55 0.55 exp(-sqrrZ0 *sqrrZ0)
1383	c #ZN. us__SV(ikl) us__SV(ikl)grvinv0.5
1384
1385	C +--Z0 Smooth + Saltat. Regime (Shao & Lin, 1999, BLM 91 (46) p.222)
1386	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^
1387	c #ZN Z0enSV(ikl) = (Z0m_nu ** sqrrZ0 )
1388	c #ZN. * (Z0mBSn **(1.-sqrrZ0))
1389	c #ZN Z0enSV(ikl) = max(Z0enSV(ikl), Z0m_nu)
1390
1391
1392	C +--Z0 Smooth Regime over Snow (Andreas etAl., 2004
1393	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^ ams.confex.com/ams/pdfpapers/68601.pdf)
1394	c #ZA Z0m_nu = 0.135*akmol / max(us__SV(ikl) , epsi)
1395
1396	C +--Z0 Saltat.Regime over Snow (Andreas etAl., 2004
1397	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^ ams.confex.com/ams/pdfpapers/68601.pdf)
1398	c #ZA Z0mBSn = 0.035u2star grvinv
1399
1400	C +--Z0 Smooth + Saltat. Regime (Andreas etAl., 2004
1401	! ( used by Erosion) ams.confex.com/ams/pdfpapers/68601.pdf)
1402	! ^^^^^^^^^^^^^^^^^^^^^^^^^^
1403	c #ZA Z0enSV(ikl) = Z0m_nu
1404	c #ZA. + Z0mBSn
1405
1406	C +--Z0 Rough Regime over Snow (Andreas etAl., 2004
1407	C + (.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_FactZ0mBSn exp(-u2star*u2star)
1411	c #ZA Z0m_90 =(10.-0.025*VVs_SV(ikl)/5.)
1412	c #ZA. exp(-0.4/sqrt(.00275+.00001max(0.,VVs_SV(ikl)-5.)))
1413	c #ZA Z0m_Sn = DDs_SV(ikl)* Z0m_90 / 45.
1414	c #ZA. - DDs_SV(ikl)DDs_SV(ikl) Z0m_90 /(90.*90.)
1415
1416
1417
1418
1419	C +--Z0 (Erosion) over Snow (instantaneous)
1420	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^
1421	Z0e_SV(ikl) = Z0enSV(ikl)
1422
1423	C +--Momentum Roughness Length (Etienne: changes wrt original SISVAT)
1424	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1425	Z0mnSV(ikl) = Z0m_nu *(1-SnoMsk) ! Ice z0
1426	. + (Z0m_Sn)*SnoMsk ! Snow Sastrugi Form and Snow Erosion
1427
1428
1429	C +--GIS Roughness Length
1430	C + ^^^^^^^^^^^^^^^^^^^^^
1431	c #GL Z0mnSV(ikl) =
1432	c #GL. (1-LSmask(ikl)) * Z0mnSV(ikl)
1433	c #GL. + LSmask(ikl) * max(Z0mnSV(ikl),max(Z0_GIM,
1434	c #GL. Z0_GIM+
1435	c #GL. (0.0032-Z0_GIM)*(ro__SV(ikl,isnoSV(ikl))-600.) !
1436	c #GL. /(920.00 -600.))) !
1437
1438	C +--Mom. Roughness Length, Instantaneous
1439	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1440	Z0m_SV(ikl) = Z0mnSV(ikl) ! Z0mnSV instant.
1441
1442
1443	C +--Roughness Length for Scalars
1444	C + ----------------------------
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+rstar1alors+rstar2alors*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
1491	c #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.

Download in other formats:

Original Format