Context Navigation

inlandsis.F @ 5097

Last change on this file since 5097 was 5082, checked in by abarral, 4 months ago
(lint) Fix obsolete boolean operators
File size: 63.9 KB

Rev	Line
[3900]	1	subroutine INLANDSIS(SnoMod,BloMod,jjtime,debut)
[3792]	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
[3900]	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
[3792]	182
	183	IMPLICIT NONE
	184
	185	logical SnoMod
	186	logical BloMod
[3900]	187	logical debut
[3792]	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
[3900]	225	real PorSno,Salt_f,PorRef !
[3792]	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
[3900]	251	real rstar,alors !
	252	real rstar0,rstar1,rstar2 !
[3792]	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
[3900]	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
[3792]	309
	310
	311	! Internal DATA
	312	! =============
	313
	314	data T__Min / 200.00/ ! Minimum realistic Temperature
[3900]	315	data TaPole / 268.15/ ! Maximum Polar Temperature (value from C. Agosta)
	316	data roSMin / 300. / ! Minimum Snow Density
[3792]	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
[3900]	333
[3792]	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)
[3900]	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
[3792]	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
[3900]	396	IF (SnoMod) THEN
[3792]	397
[3900]	398
	399	C +--Aeolian erosion and Blowing Snow
	400	C +==================================
[3792]	401
	402
	403
[3900]	404	DO ikl=1,knonv
	405	usthSV(ikl) = 1.0e+2
	406	END DO
[3792]	407
[3900]	408
	409	IF (BloMod) THEN
[3792]	410
[5082]	411	if (klonv==1) then
	412	if(isnoSV(1)>=2 .and.
[3900]	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
[3792]	416	C + **********
	417	call SISVAT_BSn
	418	endif
	419	else
	420	call SISVAT_BSn
	421	C + **********
[3900]	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
[3792]	440
	441
[3900]	442	C +--Snow Properties
	443	C + ~~~~~~~~~~~~~~~
[3792]	444
[3900]	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
[5082]	520	if (usth_param == "gal") then
[3900]	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
[5082]	527	if (usth_param == "lis") then !Liston et al. 2007
[3900]	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
[5082]	609	if (qsalt_param == "pom") then
[3900]	610	qSalSV(ikl) = (us__SV(ikl)2 - usthSV(ikl)2) *signus
	611	. / (hSalSV(ikl) * gravit * us__SV(ikl) * 3.25)
	612	endif
	613
[5082]	614	if (qsalt_param == "bin") then
[3900]	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 + ------------------------------------------------------
[3792]	635	C +--Buffer Layer
[3900]	636	C + -----------------------------------------------------
[3792]	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
[3900]	660	! C.Agosta option for snow density, same as for BS i.e.
	661	! is_ok_density_kotlyakov=.false.
[3792]	662	c #BS density_kotlyakov = .false. !C.Amory BS 2018
	663	C + ... Fallen Snow Density, Adapted for Antarctica
[3900]	664	if (is_ok_density_kotlyakov) then
[3792]	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) ]
[5082]	670	if (tt_tmp>=-10) then
[3792]	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
[3900]	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
[3792]	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 )) !
[3900]	727	! EV: now control buf_sph_pol and bug_siz_pol in physiq.def
[3792]	728	Buf_G1 = (1. - Polair) * Buf_G1 ! Temperate Snow
[3900]	729	. + Polair * buf_sph_pol ! Polar Snow
[3792]	730	Buf_G2 = (1. - Polair) * Buf_G2 ! Temperate Snow
[3900]	731	. + Polair * buf_siz_pol ! Polar Snow
[3792]	732	G1 = Buf_G1 ! NO Blown Snow
	733	G2 = Buf_G2 ! NO Blown Snow
	734
[3900]	735
	736
	737	IF (BloMod) THEN
	738
[3792]	739	! S.1. Meme Type de Neige / same Grain Type
	740	! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
[3900]	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)
[3792]	747	! Blowing Snow Properties: G1_dSV, ADSdSV
	748
	749	! S.2. Types differents / differents Types
	750	! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
[3900]	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 )) !
[3792]	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
[3900]	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
[3792]	799
[3900]	800
[3792]	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
[3900]	892	! Snow Pack Discretization(option XF in MAR)
	893	! ==========================================
[3792]	894
[3900]	895
[5082]	896	if (discret_xf.AND.klonv==1) then
[3900]	897
[5082]	898	if(isnoSV(1)>=1.or.NLaysv(1)>=1) then
[3900]	899	C + **********
	900	call SISVAT_zSn
	901	C + **********
	902	endif
	903	else
	904	C + **********
[3792]	905	call SISVAT_zSn
[3900]	906	C + **********
	907	endif
	908
	909	C + **********
[3792]	910	! #ve call SISVAT_wEq('_zSn ',0)
[3900]	911	C + **********
[3792]	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
[3900]	958	END IF ! SnoMod
[3792]	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
[3900]	999	!the albedo alb_SV and to calculate the
[3792]	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
[3900]	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
[3792]	1076
	1077
	1078
	1079	! Soil Energy Balance
	1080	! =====================
	1081
	1082
[5082]	1083	if (iflag_temp_inlandsis == 0) then
[3792]	1084
	1085	call SISVAT_TSo
	1086
	1087	else
[3900]	1088	DO ikl=1,knonv
	1089	Tsf_SV(ikl)=Tsrfsv(ikl)
	1090	END DO
[3792]	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
[5082]	1200	IF (iflag_tsurf_inlandsis == 0) THEN
[3900]	1201
	1202	Tsrfsv(ikl) =TsisSV(ikl,isnoSV(ikl))
	1203
[5082]	1204	ELSE IF (iflag_tsurf_inlandsis > 0) THEN
[3792]	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.)
[5082]	1210	else if (isnoSV(ikl) == 1) then
[3792]	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
[3900]	1224	ELSE !(default)
[3792]	1225
[3900]	1226	Tsrfsv(ikl) =TsisSV(ikl,isnoSV(ikl))
[3792]	1227
[3900]	1228	END IF
	1229
	1230
	1231	END DO
	1232
[3792]	1233	! Snow Pack Properties (sphericity, dendricity, size)
	1234	! ===================================================
	1235
	1236	IF (SnoMod) THEN
	1237
[5082]	1238	if (discret_xf .AND. klonv==1) then
	1239	if(isnoSV(1)>=1) then
[3900]	1240	C + **********
	1241	call SISVAT_GSn
	1242	C + **********
	1243	endif
	1244	else
	1245	C + **********
[3792]	1246	call SISVAT_GSn
[3900]	1247	C + **********
	1248	endif
[3792]	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 + -----------------------------
[3900]	1264
	1265	! ETIENNE WARNING: changes have been made wrt original SISVAT
[3792]	1266
	1267	C +--Land+Sea-Ice / Ice-free Sea Mask
	1268	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
[3900]	1269	DO ikl=1,knonv
[3792]	1270	IcIndx(ikl) = 0
	1271	ENDDO
	1272	DO isn=1,nsno
[3900]	1273	DO ikl=1,knonv
	1274
[3792]	1275	IcIndx(ikl) = max(IcIndx(ikl),
[3900]	1276	. isn*max(0,
	1277	. sign(1,
	1278	. int(ro__SV(ikl,isn)-900.))))
[3792]	1279	ENDDO
	1280	ENDDO
	1281
[3900]	1282	DO ikl=1,knonv
[3792]	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))
[3900]	1291
[3792]	1292	C +--Z0 Saltat.Regime over Snow (Gallee et al., 2001, BLM 99 (19) p.11)
	1293	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^
[3900]	1294
[3792]	1295	u2star = us__SV(ikl) *us__SV(ikl)
	1296	Z0mBSn = u2star *0.536e-3 - 61.8e-6
	1297	Z0mBSn = max(Z0mBS0 ,Z0mBSn)
[3900]	1298
[3792]	1299	C +--Z0 Smooth + Saltat. Regime
	1300	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^
	1301	Z0enSV(ikl) = Z0m_nu
	1302	. + Z0mBSn
[3900]	1303
	1304
	1305	! Calculation of snow roughness length
	1306	!=====================================
[5082]	1307	IF (iflag_z0m_snow == 0) THEN
[3900]	1308
	1309	Z0m_Sn=prescribed_z0m_snow
	1310
[5082]	1311	ELSE IF (iflag_z0m_snow == 1) THEN
[3900]	1312
	1313	Z0m_Sn=Z0enSV(ikl)
	1314
[5082]	1315	ELSE IF (iflag_z0m_snow == 2) THEN
[3900]	1316
[3792]	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
[3900]	1325
	1326
[3792]	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
[3900]	1337
[5082]	1338	if (TaT_SV(ikl) < ta1) then
[3792]	1339	Z0_obs = z01
[5082]	1340	else if (TaT_SV(ikl)>=ta1 .and. TaT_SV(ikl)<ta2) then
[3792]	1341	Z0_obs = exp(coefa*TaT_SV(ikl) + coefb)
[5082]	1342	else if (TaT_SV(ikl)>=ta2 .and. TaT_SV(ikl)<ta3) then
[3792]	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
[3900]	1349	Z0m_Sn=Z0_obs
[3792]	1350
[3900]	1351
	1352	ELSE
	1353
	1354	Z0m_Sn=0.500e-3 ! default=0.500e-3m (tuning of MAR)
	1355
	1356	ENDIF
[3792]	1357
[3900]	1358
	1359
	1360	! param = Z0_obs/1. ! param(s) \| 1.(m/s)=TUNING
[3792]	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
[3900]	1391
[3792]	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.)
[3900]	1415
	1416
	1417
	1418
	1419	C +--Z0 (Erosion) over Snow (instantaneous)
[3792]	1420	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^
	1421	Z0e_SV(ikl) = Z0enSV(ikl)
	1422
[3900]	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
[3792]	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
[3900]	1438	C +--Mom. Roughness Length, Instantaneous
	1439	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
[3792]	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
[3900]	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))))
[3792]	1468
[3900]	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
[3792]	1473
[3900]	1474	if(ro__SV(ikl,isnoSV(ikl))>50
	1475	. .and.ro__SV(ikl,isnoSV(ikl))<roSdSV)then
[3792]	1476
[3900]	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)))
[3792]	1482
[3900]	1483	endif
	1484
	1485
	1486	ENDIF
[3792]	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