Context Navigation

inlandsis.f90 @ 5159

Last change on this file since 5159 was 5117, checked in by abarral, 5 months ago
rename modules properly lmdz_* move some unused files to obsolete/ (lint) uppercase fortran keywords
File size: 59.5 KB

Rev	Line
[5105]	1	SUBROUTINE INLANDSIS(SnoMod,BloMod,jjtime,debut)
[3792]	2
[5105]	3	USE dimphy
[3792]	4
[5105]	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	!--------------------------------------------------------------------------+
[3792]	161
	162
	163
[5105]	164	! Global Variables
	165	! ================
[3792]	166
	167
[5105]	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
[3792]	182
[5105]	183	IMPLICIT NONE
[3792]	184
[5117]	185	LOGICAL :: SnoMod
	186	LOGICAL :: BloMod
	187	LOGICAL :: debut
[5116]	188	INTEGER :: jjtime
[3792]	189
	190
[5105]	191	! Internal Variables
	192	! ==================
[3792]	193
[5105]	194	! Non Local
	195	! ---------
[3792]	196
[5116]	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
[3792]	207
	208
	209
[5105]	210	! Local
	211	! -----
[3792]	212
[5116]	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 !
[5105]	226	! #sw real PorVol,rWater !
	227	! #sw real rusNEW,rdzNEW,etaNEW !
[5116]	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
[5105]	237	! #ZN real sqrrZ0 ! ut/u
[5116]	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
[5105]	245	! #SZ real Z0Sa_N ! Regime Snow Roughness Length
	246	! #SZ real Z0SaSi ! 1.IF Rgm Snow Roughness Length
	247	! #GL real Z0_GIM ! Mimimum GIMEX Roughness Length
[5116]	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
[5105]	285	! for Aeolian erosion and blowing snow
[5116]	286	INTEGER :: nit ,iit
	287	REAL :: Fac ! Correc. factor for drift ratio
	288	REAL :: dusuth,signus
	289	REAL :: sss__F,sss__N
	290	REAL :: sss__K,sss__G
	291	REAL :: us_127,us_227,us_327,us_427,us_527
	292	REAL :: VVa_OK, usuth0
	293	REAL :: ssstar
	294	REAL :: SblPom
	295	REAL :: rCd10n ! Square root of drag coefficient
	296	REAL :: DendOK ! Dendricity Switch
	297	REAL :: SaltOK ! Saltation Switch
	298	REAL :: MeltOK ! Saltation Switch (Melting Snow)
	299	REAL :: SnowOK ! Pack Top Switch
	300	REAL :: SaltM1,SaltM2,SaltMo,SaltMx ! Saltation Parameters
	301	REAL :: ShearX, ShearS ! Arg. Max Shear Stress
	302	REAL :: Por_BS ! Snow Porosity
	303	REAL :: Salt_us ! New thresh.friction velocity u*t
	304	REAL :: Fac_Mo,ArguSi,FacRho ! Numerical factors for u*t
	305	REAL :: SaltSI(klonv,0:nsno) ! Snow Drift Index !
	306	REAL :: MIN_Mo ! Minimum Mobility Fresh Fallen *
	307	CHARACTER(LEN=3) :: qsalt_param ! Switch for saltation flux param.
	308	CHARACTER(LEN=3) :: usth_param ! Switch for u*t param
[3792]	309
	310
[5105]	311	! Internal DATA
	312	! =============
[3900]	313
[5105]	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
[3792]	333
[3900]	334
[5105]	335	! DATA Emissivities ! Pielke, 1984, pp. 383,409
[3792]	336
[5105]	337	data Z0mBS0 / 0.5e-6/ ! MINimum Snow Roughness Length
[5113]	338	! for Momentum if Blowing Snow
	339	! Gallee et al. 2001 BLM 99 (19)
[5105]	340	data Z0m_S0/ 0.00005/ ! MINimum Snow Roughness Length
[5113]	341	! MegaDunes included
[5105]	342	data Z0m_S1/ 0.030 / ! MAXimum Snow Roughness Length
[5113]	343	! (Sastrugis)
[5105]	344	! #GL data Z0_GIM/ 0.0013/ ! Ice Min Z0 = 0.0013 m (Broeke)
[5113]	345	! Old Ice Z0 = 0.0500 m (Bruce)
	346	! 0.0500 m (Smeets)
	347	! 0.1200 m (Broeke)
[5105]	348	data Z0_ICE/ 0.0010/ ! Sea-Ice Z0 = 0.0010 m (Andreas)
[5113]	349	! (Ice Station Weddel -- ISW)
[5105]	350	! for aerolian erosion
	351	data SblPom/ 1.27/ ! Lower Boundary Height Parameter
	352	! + ! for Suspension
	353	! + ! Pommeroy, Gray and Landine 1993,
	354	! + ! J. Hydrology, 144(8) p.169
	355	data nit / 5 / ! us(is0,uth) recursivity: Nb Iterations
	356	!c#AE data qsalt_param/"bin"/ ! saltation part. conc. from Bintanja 2001 (p
	357	data qsalt_param/"pom"/ ! saltation part. conc. from Pomeroy and Gray
	358	!c#AE data usth_param/"lis"/ ! u*t from Liston et al. 2007
	359	data usth_param/"gal"/ ! u*t from Gallee et al. 2001
	360	data SaltMx/-5.83e-2/
[3792]	361
[5105]	362	vk2 = vonKrm * vonKrm ! Square of Von Karman Constant
[3792]	363
	364
[5105]	365	! BEGIN.main.
	366	! ===========================
[3792]	367
	368
	369
	370
[5105]	371	! "Soil" Humidity of Water Bodies
	372	! ===============================
[3792]	373
[5105]	374	DO ikl=1,knonv
[3792]	375
[5105]	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
[3792]	383
	384
[5105]	385	! Vertical Discretization Factor
	386	! ==============================
[3792]	387
[5105]	388	LSdzsv(ikl) = ist__s & ! Soil
	389	+ OcndSV * ist__w ! Water Body
	390	END DO
[3792]	391
	392
	393
	394
	395
[5105]	396	IF (SnoMod) THEN
[3792]	397
	398
[5105]	399	! +--Aeolian erosion and Blowing Snow
	400	! +==================================
[3792]	401
[3900]	402
	403
[5105]	404	DO ikl=1,knonv
	405	usthSV(ikl) = 1.0e+2
	406	END DO
[3900]	407
	408
[5105]	409	IF (BloMod) THEN
[3900]	410
[5117]	411	IF (klonv==1) THEN
	412	IF(isnoSV(1)>=2 .AND. &
	413	TsisSV(1,max(1,isnoSV(1)))<273. .AND. &
	414	ro__SV(1,max(1,isnoSV(1)))<500. .AND. &
[5116]	415	eta_SV(1,max(1,isnoSV(1)))<epsi) THEN
[5105]	416	! + **********
	417	CALL SISVAT_BSn
	418	endif
	419	else
	420	CALL SISVAT_BSn
	421	! + **********
	422	endif
[3900]	423
	424
	425
	426
	427
[3792]	428
	429
[5105]	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
[3900]	433
[5105]	434	! +--Computation of threshold friction velocity for snow erosion
	435	! ---------------------------------------------------------------
[3900]	436
[5105]	437	rCd10n = 1. / 26.5 ! Vt / u*t = 26.5
[5113]	438	! Budd et al. 1965, Antarct. Res. Series Fig.13
	439	! ratio developped during assumed neutral conditions
[3900]	440
	441
[5105]	442	! +--Snow Properties
	443	! + ~~~~~~~~~~~~~~~
[3900]	444
[5105]	445	DO ikl = 1,knonv
[3900]	446
[5105]	447	isn = isnoSV(ikl)
[3900]	448
	449
	450
[5105]	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
[3900]	460
[5105]	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)
[3900]	465
[5105]	466	SaltOK = min(unun, SaltOK + MeltOK) * SnowOK
[3900]	467
[5105]	468
	469	! +--Mobility Index (Guyomarc'h & Merindol 1997, Ann.Glaciol.)
	470	! + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	471	SaltM1 = -0.750e-2 * G1snSV(ikl,isn) &
	472	-0.500e-2 * G2snSV(ikl,isn)+ 0.500e00 !dendritic case
	473	! + CAUTION: Guyomarc'h & Merindol Dendricity Sign is +
	474	! + ^^^^^^^^ MAR Dendricity Sign is -
	475	SaltM2 = -0.833d-2 * G1snSV(ikl,isn) &
	476	-0.583d-2 * G2snSV(ikl,isn)+ 0.833d00 !non-dendritic case
	477
	478	! SaltMo = (DendOK * SaltM1 + (1.-DendOK) * SaltM2 )
	479	SaltMo = 0.625 !SaltMo pour d=s=0.5
	480
	481	!weighting SaltMo with surface snow density (Vionnet et al. 2012)
	482	!c#AE FacRho = 1.25 - 0.0042 * ro__SV(ikl,isn)
	483	!c#AE SaltMo = 0.34 * SaltMo + 0.66 * FacRho !needed for polar snow
	484	MIN_Mo = 0.
	485	! SaltMo = max(SaltMo,MIN_Mo)
	486	! SaltMo = SaltOK * SaltMo + (1.-SaltOK) * min(SaltMo,SaltMx)
	487	! #TUNE SaltMo = SaltOK * SaltMo - (1.-SaltOK) * 0.9500
	488	SaltMo = max(SaltMo,epsi-unun)
	489
	490	! +--Influence of Density on Threshold Shear Stress
	491	! + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	492	Por_BS = 1. - 300. / ro_Ice
	493	ShearS = Por_BS / (1.-Por_BS)
	494	! +... SheaBS = Arg(sqrt(shear = max shear stress in snow)):
	495	! + shear = 3.420d00 * exp(-(Por_BS +Por_BS)
	496	! + . /(unun -Por_BS))
	497	! + SheaBS : see de Montmollin (1978),
	498	! + These Univ. Sci. Medic. Grenoble, Fig. 1 p. 124
	499
	500	! +--Snow Drift Index (Guyomarc'h & Merindol 1997, Ann.Glaciol.)
	501	! + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	502	ArguSi = -0.085 *us__SV(ikl)/rCd10n
	503	!V=u*/sqrt(CD) eqs 2 to 4 Gallee et al. 2001
	504
	505	SaltSI(ikl,isn) = -2.868 * exp(ArguSi) + 1 + SaltMo
	506
	507
	508	! +--Threshold Friction Velocity
	509	! + ~~~~~~~~~~~~~~~~~~~~~~~~~~~
[5116]	510	IF(ro__SV(ikl,isn)>300.) THEN
[5105]	511	Por_BS = 1.000 - ro__SV(ikl,isn) /ro_Ice
	512	else
	513	Por_BS = 1.000 - 300. /ro_Ice
	514	endif
	515
	516	ShearX = Por_BS/max(epsi,1.-Por_BS)
	517	Fac_Mo = exp(-ShearX+ShearS)
	518	! + Gallee et al., 2001 eq 5, p5
	519
[5117]	520	IF (usth_param == "gal") THEN
[5105]	521	Salt_us = (log(2.868) - log(1 + SaltMo)) * rCd10n/0.085
	522	Salt_us = Salt_us * Fac_Mo
	523	! +... Salt_us : Extension of Guyomarc'h & Merindol 1998 with
	524	! +... de Montmollin (1978). Gallee et al. 2001
	525	endif
	526
[5117]	527	IF (usth_param == "lis") then !Liston et al. 2007
[5116]	528	IF(ro__SV(ikl,isn)>300.) THEN
[5105]	529	Salt_us = 0.005exp(0.013ro__SV(ikl,isn))
	530	else
	531	Salt_us = 0.01exp(0.003ro__SV(ikl,isn))
[3900]	532	endif
[5105]	533	endif
[3900]	534
[5105]	535	SnowOK = 1 -min(1,iabs(isn-isnoSV(ikl))) !Switch new vs old snow
[3900]	536
[5105]	537	usthSV(ikl) = SnowOK * (Salt_us) &
	538	+ (1.-SnowOK)* usthSV(ikl)
[3900]	539
[5105]	540	END DO
[3900]	541
[5105]	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
[5113]	556	! ______________ _____
	557	! Newton-Raphson (! Iteration, BEGIN)
	558	! ~~~~~~~~~~~~~~ ~~~~~
[5105]	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
[5113]	582	! #AE sss__F = vonkar /rCDmSV(ikl)
[5105]	583	! ENDDO
	584
[5113]	585	! ______________ ___
	586	! Newton-Raphson (! Iteration, END )
	587	! ~~~~~~~~~~~~~~ ~~~
[5105]	588
	589	! us_127 = exp( SblPom *log(us__SV(ikl)))
	590	! us_227 = us_127 * us__SV(ikl)
	591
[5113]	592	! Momentum Turbulent Scale u*: 0-Limit in case of no Blow. Snow
	593	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[5105]	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
[5117]	609	IF (qsalt_param == "pom") THEN
[5105]	610	qSalSV(ikl) = (us__SV(ikl)2 - usthSV(ikl)2) *signus &
	611	/ (hSalSV(ikl) * gravit * us__SV(ikl) * 3.25)
	612	endif
	613
[5117]	614	IF (qsalt_param == "bin") THEN
[5105]	615	qSalSV(ikl) = (us__SV(ikl) * us__SV(ikl) &
	616	-usthSV(ikl) * usthSV(ikl))*signus &
	617	* 0.535 / (hSalSV(ikl) * gravit)
	618	endif
	619
	620	qSalSV(ikl) = qSalSV(ikl)/rht_SV(ikl) ! conversion kg/m3 to kg/kg
	621
	622	ssstar = rCDmSV(ikl) * (qsnoSV(ikl) - qSalSV(ikl)) &
	623	* r_Turb !Bintanja 2000, BLM
	624	!r_Turb compensates for an overestim. of the blown snow part. fall velocity
	625
	626	uss_SV(ikl) = min(zero , us__SV(ikl) *ssstar)
	627	uss_SV(ikl) = max(-0.0001 , uss_SV(ikl))
	628
	629
	630
	631
	632	ENDIF ! BloMod
	633
	634	! + ------------------------------------------------------
	635	! +--Buffer Layer
	636	! + -----------------------------------------------------
	637
	638	DO ikl=1,knonv
	639	! BufsSV(ikl) [mm w.e.] i.e, i.e., [kg/m2]
	640	d_Bufs = max(dsn_SV(ikl) *dt__SV,0.) !
	641	dsn_SV(ikl) = 0. !
	642	Bufs_N = BufsSV(ikl) +d_Bufs !
	643
	644
	645	! +--Snow Density
	646	! + ^^^^^^^^^^^^
	647	Polair = zero
	648	! #NP Polair = max(zero, !
	649	! #NP. sign(unun,TaPole !
	650	! #NP. -TaT_SV(ikl))) !
	651	Polair = max(zero, & !
	652	sign(unun,TaPole & !
	653	-TaT_SV(ikl))) !
	654	Buf_ro = max( rosMin, & ! Fallen Snow Density
	655	roSn_1+roSn_2* (TaT_SV(ikl)-TfSnow) & ! [kg/m3]
	656	+roSn_3*sqrt( VV__SV(ikl))) ! Pahaut (CEN), Etienne: use wind speed at first model level instead of 10m wind
	657	! #NP BufPro = max( rosMin, ! Fallen Snow Density
	658	! #NP. 104. *sqrt( max( VV10SV(ikl)-6.0,0.0))) ! Kotlyakov (1961)
	659
	660	! C.Agosta option for snow density, same as for BS i.e.
	661	! is_ok_density_kotlyakov=.FALSE.
	662	! #BS density_kotlyakov = .FALSE. !C.Amory BS 2018
	663	! + ... Fallen Snow Density, Adapted for Antarctica
[5117]	664	IF (is_ok_density_kotlyakov) THEN
[5105]	665	tt_tmp = TaT_SV(ikl)-TfSnow
[5113]	666	!vv_tmp = VV10SV(ikl)
[5105]	667	vv_tmp=VV__SV(ikl) ! Etienne: use wind speed at first model level instead of 10m wind
	668	! + ... [ A compromise between
	669	! + ... Kotlyakov (1961) and Lenaerts (2012, JGR, Part1) ]
[5117]	670	IF (tt_tmp>=-10) THEN
[5105]	671	BufPro = max( rosMin, &
	672	104. *sqrt( max( vv_tmp-6.0,0.0))) ! Kotlyakov (1961)
[3792]	673	else
[5105]	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))
[3792]	677	endif
[5105]	678	else
	679	! + ... [ density derived from observations of the first 50cm of
	680	! + ... snow - cf. Rajashree Datta - and multiplied by 0.8 ]
	681	! + ... C. Agosta, 2016-09
	682	!c #SD BufPro = 149.2 + 6.84VV10SV(ikl) + 0.48Tsrfsv(ikl)
	683	!c #SD BufPro = 125 + 14VV10SV(ikl) + 0.6Tsrfsv(ikl) !MAJ CK and CAm
	684	! BufPro = 200 + 21 * VV10SV(ikl)!CK 29/07/19
	685	BufPro = 200 + 21 * VV__SV(ikl)!Etienne: use wind speed at first model level instead of 10m wind
	686	endif
[3792]	687
[5105]	688	Bros_N = (1. - Polair) * Buf_ro & ! Temperate Snow
	689	+ Polair * BufPro ! Polar Snow
[3792]	690
[5105]	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
[3900]	693
	694
[5105]	695	! Density of deposited blown snow
	696	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[3900]	697
[5117]	698	IF (BloMod) THEN
[5105]	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
[3900]	709
[3792]	710
[5105]	711	! Time averaged Density of deposited blown Snow
	712	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[3900]	713
[5105]	714	BrosSV(ikl) =(Bros_N * d_Bufs & !
	715	+BrosSV(ikl)* BufsSV(ikl)) & !
	716	/ max(epsi,Bufs_N) !
[3792]	717
	718
[5105]	719	! +-- S.Falling Snow Properties (computed as in SISVAT_zAg)
	720	! + ^^^^^^^^^^^^^^^^^^^^^^^
	721	Buf_G1 = max(-G1_dSV, & ! Temperate Snow
	722	min(Dendr1*VV__SV(ikl)-Dendr2, & ! Dendricity
	723	Dendr3 )) !
	724	Buf_G2 = min( Spher4, & ! Temperate Snow
	725	max(Spher1*VV__SV(ikl)+Spher2, & ! Sphericity
	726	Spher3 )) !
	727	! EV: now control buf_sph_pol and bug_siz_pol in physiq.def
	728	Buf_G1 = (1. - Polair) * Buf_G1 & ! Temperate Snow
	729	+ Polair * buf_sph_pol ! Polar Snow
	730	Buf_G2 = (1. - Polair) * Buf_G2 & ! Temperate Snow
	731	+ Polair * buf_siz_pol ! Polar Snow
	732	G1 = Buf_G1 ! NO Blown Snow
	733	G2 = Buf_G2 ! NO Blown Snow
[3792]	734
	735
	736
[5105]	737	IF (BloMod) THEN
[3900]	738
[5105]	739	! S.1. Meme Type de Neige / same Grain Type
	740	! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
[3792]	741
[5105]	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
[3792]	748
[5105]	749	! S.2. Types differents / differents Types
	750	! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	751	typ__1 = max(zero,sign(unun,epsi-Buf_G1)) ! =1.=> Dendritic
	752	zroNEW = typ__1 *(1.0-dsnbSV(ikl)) & ! fract.Dendr.Lay.
	753	+ (1.-typ__1) * dsnbSV(ikl) !
	754	G1_NEW = typ__1 *Buf_G1 & ! G1 of Dendr.Lay.
	755	+ (1.-typ__1) *G1_dSV !
	756	G2_NEW = typ__1 *Buf_G2 & ! G2 of Dendr.Lay.
	757	+ (1.-typ__1) *ADSdSV !
	758	zroOLD = (1.-typ__1) *(1.0-dsnbSV(ikl)) & ! fract.Spher.Lay.
	759	+ typ__1 * dsnbSV(ikl) !
	760	G1_OLD = (1.-typ__1) *Buf_G1 & ! G1 of Spher.Lay.
	761	+ typ__1 *G1_dSV !
	762	G2_OLD = (1.-typ__1) *Buf_G2 & ! G2 of Spher.Lay.
	763	+ typ__1 *ADSdSV !
	764	SizNEW = -G1_NEW *DDcdSV/G1_dSV & ! Size Dendr.Lay.
	765	+(1.+G1_NEW /G1_dSV) & !
	766	(G2_NEW DScdSV/G1_dSV & !
	767	+(1.-G2_NEW /G1_dSV)*DFcdSV) !
	768	SphNEW = G2_NEW /G1_dSV ! Spher.Dendr.Lay.
	769	SizOLD = G2_OLD ! Size Spher.Lay.
	770	SphOLD = G1_OLD /G1_dSV ! Spher.Spher.Lay.
	771	Siz_av = (zroNEWSizNEW+zroOLDSizOLD) ! Averaged Size
	772	Sph_av = min( zroNEWSphNEW+zroOLDSphOLD & !
	773	, unun) ! Averaged Sphericity
	774	Den_av = min((Siz_av -( Sph_av *DScdSV & !
	775	+(1.-Sph_av)*DFcdSV)) & !
	776	/ (DDcdSV -( Sph_av *DScdSV & !
	777	+(1.-Sph_av)*DFcdSV)) & !
	778	, unun) !
	779	DendOK = max(zero, & !
	780	sign(unun, Sph_av *DScdSV & ! Small Grains
	781	+(1.-Sph_av)*DFcdSV & ! Faceted Grains
	782	- Siz_av )) !
	783	! +... REMARQUE: le type moyen (dendritique ou non) depend
	784	! + ^^^^^^^^ de la comparaison avec le diametre optique
	785	! + d'une neige recente de dendricite nulle
	786	! +... REMARK: the mean type (dendritic or not) depends
	787	! + ^^^^^^ on the comparaison with the optical diameter
	788	! + of a recent snow having zero dendricity
[3792]	789
[5105]	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
	800
	801
[5105]	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)
[3792]	811
[5105]	812	! S.2. Types differents / differents Types
	813	! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
[3792]	814
[5105]	815	typ__1 = max(zero,sign(unun,epsi-Buf_G1)) ! =1.=> Dendritic
	816	zroNEW =( typ__1 *d_Bufs & ! fract.Dendr.Lay.
	817	+ (1.-typ__1) *BufsSV(ikl)) & !
	818	/max(epsi,Bufs_N) !
	819	G1_NEW = typ__1 *Buf_G1 & ! G1 of Dendr.Lay.
	820	+ (1.-typ__1) *BG1sSV(ikl) !
	821	G2_NEW = typ__1 *Buf_G2 & ! G2 of Dendr.Lay.
	822	+ (1.-typ__1) *BG2sSV(ikl) !
	823	zroOLD =((1.-typ__1) *d_Bufs & ! fract.Spher.Lay.
	824	+ typ__1 *BufsSV(ikl)) & !
	825	/max(epsi,Bufs_N) !
	826	G1_OLD = (1.-typ__1) *Buf_G1 & ! G1 of Spher.Lay.
	827	+ typ__1 *BG1sSV(ikl) !
	828	G2_OLD = (1.-typ__1) *Buf_G2 & ! G2 of Spher.Lay.
	829	+ typ__1 *BG2sSV(ikl) !
	830	SizNEW = -G1_NEW *DDcdSV/G1_dSV & ! Size Dendr.Lay.
	831	+(1.+G1_NEW /G1_dSV) & !
	832	(G2_NEW DScdSV/G1_dSV & !
	833	+(1.-G2_NEW /G1_dSV)*DFcdSV) !
	834	SphNEW = G2_NEW /G1_dSV ! Spher.Dendr.Lay.
	835	SizOLD = G2_OLD ! Size Spher.Lay.
	836	SphOLD = G1_OLD /G1_dSV ! Spher.Spher.Lay.
	837	Siz_av = ( zroNEW SizNEW+zroOLDSizOLD) ! Averaged Size
	838	Sph_av = min( zroNEW SphNEW+zroOLDSphOLD & !
	839	, unun ) ! Averaged Sphericity
	840	Den_av = min((Siz_av - ( Sph_av *DScdSV & !
	841	+(1.-Sph_av)*DFcdSV)) & !
	842	/ (DDcdSV - ( Sph_av *DScdSV & !
	843	+(1.-Sph_av)*DFcdSV)) & !
	844	, unun )!
	845	DendOK = max(zero, & !
	846	sign(unun, Sph_av *DScdSV & ! Small Grains
	847	+(1.-Sph_av)*DFcdSV & ! Faceted Grains
	848	- Siz_av )) !
	849	! +... REMARQUE: le type moyen (dendritique ou non) depend
	850	! + ^^^^^^^^ de la comparaison avec le diametre optique
	851	! + d'une neige recente de dendricite nulle
	852	! +... REMARK: the mean type (dendritic or not) depends
	853	! + ^^^^^^ on the comparaison with the optical diameter
	854	! + of a recent snow having zero dendricity
[3792]	855
[5105]	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
[3792]	864
[5105]	865	BG1sSV(ikl) = G1 & !
	866	* Bufs_N/max(epsi,Bufs_N) !
	867	BG2sSV(ikl) = G2 & !
	868	* Bufs_N/max(epsi,Bufs_N) !
[3792]	869
	870
[5105]	871	! +--Update of Buffer Layer Content & Decision about creating a new snow layer
	872	! + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	873	BufsSV(ikl) = Bufs_N ! [mm w.e.]
	874	NLaysv(ikl) = min(unun, & !
	875	max(zero, & ! Allows to create
	876	sign(unun,BufsSV(ikl) & ! a new snow Layer
	877	-SMndSV )) & ! if Buffer > SMndSV
	878	max(zero, & ! Except if Erosion
	879	sign(unun,0.50 & ! dominates
	880	-dsnbSV(ikl))) & !
	881	+max(zero, & ! Allows to create
	882	sign(unun,BufsSV(ikl) & ! a new snow Layer
	883	-SMndSV3.00))) ! is Buffer > SMndSV3
	884	Bdzssv(ikl) = 1.e-3BufsSV(ikl)ro_Wat & ! [mm w.e.] -> [m w.e.]
	885	/max(epsi,BrosSV(ikl))!& [m w.e.] -> [m]
[3792]	886
	887
[5105]	888	END DO
[3792]	889
	890
	891
[5105]	892	! Snow Pack Discretization(option XF in MAR)
	893	! ==========================================
[3792]	894
	895
[5117]	896	IF (discret_xf.AND.klonv==1) THEN
	897	IF(isnoSV(1)>=1.OR.NLaysv(1)>=1) THEN
[5105]	898	! + **********
	899	CALL SISVAT_zSn
	900	! + **********
	901	endif
	902	else
	903	! + **********
	904	CALL SISVAT_zSn
	905	! + **********
[5117]	906	ENDIF
[3792]	907
[5105]	908	! + **********
	909	! #ve CALL SISVAT_wEq('_zSn ',0)
	910	! + **********
[3792]	911
[5105]	912	! Add a new Snow Layer
	913	! ====================
[3792]	914
[5105]	915	DO ikl=1,knonv
[3792]	916
[5105]	917	isnoSV(ikl) = isnoSV(ikl) +NLaysv(ikl)
	918	isn = isnoSV(ikl)
	919	dzsnSV(ikl,isn) = dzsnSV(ikl,isn) * (1-NLaysv(ikl)) &
	920	+ Bdzssv(ikl) * NLaysv(ikl)
	921	TsisSV(ikl,isn) = TsisSV(ikl,isn) * (1-NLaysv(ikl)) &
	922	+ min(TaT_SV(ikl),Tf_Sno) *NLaysv(ikl)
	923	ro__SV(ikl,isn) = ro__SV(ikl,isn) * (1-NLaysv(ikl)) &
	924	+ Brossv(ikl) * NLaysv(ikl)
	925	eta_SV(ikl,isn) = eta_SV(ikl,isn) * (1-NLaysv(ikl)) ! + 0.
	926	agsnSV(ikl,isn) = agsnSV(ikl,isn) * (1-NLaysv(ikl)) ! + 0.
	927	G1snSV(ikl,isn) = G1snSV(ikl,isn) * (1-NLaysv(ikl)) &
	928	+ BG1ssv(ikl) * NLaysv(ikl)
	929	G2snSV(ikl,isn) = G2snSV(ikl,isn) * (1-NLaysv(ikl)) &
	930	+ BG2ssv(ikl) * NLaysv(ikl)
	931	istoSV(ikl,isn) = istoSV(ikl,isn) * (1-NLaysv(ikl)) &
	932	+ max(zer0,sign(un_1,TaT_SV(ikl) &
	933	-Tf_Sno-eps_21)) * istdSV(2) &
	934	* NLaysv(ikl)
	935	BufsSV(ikl) = BufsSV(ikl) * (1-NLaysv(ikl))
	936	NLaysv(ikl) = 0
[3792]	937
	938
[5105]	939	END DO
[3792]	940
	941
[5105]	942	! Snow Pack Thickness
	943	! -------------------
	944
	945	DO ikl=1,knonv
	946	z_snsv(ikl) = 0.0
	947	END DO
	948	DO isn=1,nsno
	949	DO ikl=1,knonv
	950	z_snsv(ikl) = z_snsv(ikl) + dzsnSV(ikl,isn)
	951	zzsnsv(ikl,isn) = z_snsv(ikl)
[3792]	952	END DO
[5105]	953	END DO
[3792]	954
	955
[5099]	956
[5105]	957	END IF ! SnoMod
[3792]	958
	959
	960
[5105]	961	! Soil Albedo: Soil Humidity Correction
	962	! ==========================================
[3792]	963
[5105]	964	! REFERENCE: McCumber and Pielke (1981), Pielke (1984)
	965	! ^^^^^^^^^
	966	DO ikl=1,knonv
	967	albssv(ikl) = &
	968	alb0SV(ikl) *(1.0-min(half,eta_SV( ikl,0) &
	969	/etadSV(isotSV(ikl))))
	970	! REMARK: Albedo of Water Surfaces (isotSV=0):
	971	! ^^^^^^ alb0SV := 2 X effective value, while
	972	! eta_SV := etadSV
	973	END DO
[3792]	974
	975
[5105]	976	! Snow Pack Optical Properties
	977	! ============================
[3792]	978
[5105]	979	IF (SnoMod) THEN
[3792]	980
[5105]	981	! ******
	982	CALL SnOptP(jjtime)
	983	! ******
[3792]	984
[5105]	985	ELSE
	986	DO ikl=1,knonv
	987	sEX_sv(ikl,1) = 1.0
	988	sEX_sv(ikl,0) = 0.0
	989	albisv(ikl) = albssv(ikl)
	990	END DO
	991	END IF
[3792]	992
	993
	994
[5105]	995	! Soil optical properties
	996	! =============================
	997	!Etienne: as in inlandis we do not CALL vgopt, we need to define
	998	!the albedo alb_SV and to calculate the
	999	!absorbed Solar Radiation by Surfac (Normaliz)[-] SoSosv
[3792]	1000
	1001
[5105]	1002	DO ikl=1,klonv
[3792]	1003
[5105]	1004	e_pRad = 2.5 * coszSV(ikl) ! exponential argument,
[5113]	1005	! V/nIR radiation partitioning,
	1006	! DR97, 2, eqn (2.53) & (2.54)
[5105]	1007	e1pRad = 1.-exp(-e_pRad) ! exponential, V/nIR Rad. Part.
	1008	exdRad= 1.
[3792]	1009
[5105]	1010	! Visible Part of the Solar Radiation Spectrum (V, 0.4--0.7mi.m)
	1011	A_Rad0 = 0.25 + 0.697 * e1pRad ! Absorbed Vis. Radiation
	1012	absg_V = (1.-albisv(ikl))(A_Rad0exdRad) !
[3792]	1013
[5105]	1014	! Near-IR Part of the Solar Radiation Spectrum (nIR, 0.7--2.8mi.m)
[3792]	1015
[5105]	1016	A_Rad0 = 0.80 + 0.185 * e1pRad ! Absorbed nIR. Radiation
	1017	absgnI = (1.-albisv(ikl))(A_Rad0exdRad) !
[3792]	1018
[5105]	1019	SoSosv(ikl) = (absg_V+absgnI)*0.5d0
[3792]	1020
[5105]	1021	alb_SV(ikl) = albisv(ikl)
[3792]	1022
[5105]	1023	END DO
[3792]	1024
[5105]	1025	! **********
	1026	! #ve CALL SISVAT_wEq('SnOptP',0)
	1027	! **********
[3792]	1028
	1029
[5105]	1030	! Surface Emissivity (Etienne: simplified calculation for landice)
	1031	! =============================================================
[3792]	1032
[5105]	1033	DO ikl=1,knonv
	1034	LSnMsk = min( 1,isnoSV(ikl))
	1035	Eso_sv(ikl)= EmiSol(1-LSnMsk)+EmiSnoLSnMsk ! Sol+Sno Emissivity
	1036	emi_SV(ikl)= EmiSol(1-LSnMsk) + EmiSnoLSnMsk
	1037	END DO
[3792]	1038
	1039
	1040
	1041
[5105]	1042	! Upward IR (INPUT, from previous time step)
	1043	! ===================================================================
[3792]	1044
[5105]	1045	DO ikl=1,knonv
	1046	! #e1 Enrsvd(ikl) = - IRs_SV(ikl)
	1047	IRupsv(ikl) = IRs_SV(ikl)
	1048	END DO
[3792]	1049
[5105]	1050
	1051	! Turbulence
	1052	! ==========
	1053
	1054	! Latent Heat of Vaporization/Sublimation
	1055	! ---------------------------------------
	1056
	1057	DO ikl=1,knonv
	1058	SnoWat = min(isnoSV(ikl),0)
	1059	Lx_H2O(ikl) = &
	1060	(1.-SnoWat) * LhvH2O &
	1061	+ SnoWat (LhsH2O (1.-eta_SV(ikl,isnoSV(ikl))) &
	1062	+LhvH2O * eta_SV(ikl,isnoSV(ikl)) )
	1063	END DO
	1064
	1065
	1066
	1067
	1068	! Aerodynamic Resistance (calculated from drags given by LMDZ)
	1069	! Commented because already calculated in surf_inlandsis_mod
	1070	! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	1071	! DO ikl=1,knonv
	1072	! ram_sv(ikl) = 1./(cdM_SV(ikl)*max(VV__SV(ikl),eps6))
	1073	! rah_sv(ikl) = 1./(cdH_SV(ikl)*max(VV__SV(ikl),eps6))
	1074	! END DO
	1075
	1076
	1077
	1078	! Soil Energy Balance
	1079	! =====================
	1080
	1081
[5117]	1082	IF (iflag_temp_inlandsis == 0) THEN
[5105]	1083	CALL SISVAT_TSo
	1084
	1085	else
	1086	DO ikl=1,knonv
	1087	Tsf_SV(ikl)=Tsrfsv(ikl)
	1088	END DO
	1089
	1090	CALL SISVAT_TS2
	1091
	1092	end if
	1093
	1094
	1095	! **********
	1096	! #ve CALL SISVAT_wEq('_TSo ',0)
	1097	! **********
	1098
	1099
	1100
	1101	! Soil Water Potential
	1102	! ------------------------
	1103
	1104	DO isl=-nsol,0
	1105	DO ikl=1,knonv
	1106	ist = isotSV(ikl) ! Soil Type
	1107	psi_sv(ikl,isl) = psidSV(ist) & ! DR97, Eqn.(3.34)
	1108	*(etadSV(ist) /max(eps6,eta_SV(ikl,isl))) & !
	1109	**bCHdSV(ist) !
	1110
	1111
	1112	! Soil Hydraulic Conductivity
	1113	! ---------------------------
	1114
	1115	Khydsv(ikl,isl) = s2__SV(ist) & ! DR97, Eqn.(3.35)
	1116	(eta_SV(ikl,isl)(2.bCHdSV(ist)+3.)) !
	1117	END DO
	1118	END DO
	1119
	1120
	1121	! Melting / Refreezing in the Snow Pack
	1122	! =====================================
	1123
	1124	IF (SnoMod) THEN
	1125
	1126	! **********
	1127	CALL SISVAT_qSn
	1128	! **********
	1129
	1130	! **********
	1131	! #ve CALL SISVAT_wEq('_qSn ',0)
	1132	! **********
	1133
	1134
	1135	! Snow Pack Thickness
	1136	! -------------------
	1137
	1138	DO ikl=1,knonv
	1139	z_snsv(ikl) = 0.0
[3792]	1140	END DO
[5105]	1141	DO isn=1,nsno
	1142	DO ikl=1,knonv
	1143	z_snsv(ikl) = z_snsv(ikl) + dzsnSV(ikl,isn)
	1144	zzsnsv(ikl,isn) = z_snsv(ikl)
	1145	END DO
	1146	END DO
[3792]	1147
	1148
[5105]	1149	! Energy in Excess is added to the first Soil Layer
	1150	! -------------------------------------------------
	1151	DO ikl=1,knonv
	1152	z_snsv(ikl) = max(zer0, &
	1153	sign(un_1,eps6-z_snsv(ikl)))
	1154	TsisSV(ikl,0) = TsisSV(ikl,0) + EExcsv(ikl) &
	1155	/(rocsSV(isotSV(ikl)) &
	1156	+rcwdSV*eta_SV(ikl,0))
	1157	EExcsv(ikl) = 0.
	1158	END DO
[3792]	1159
	1160
[5105]	1161	END IF
[3792]	1162
	1163
[5105]	1164	! Soil Water Balance
	1165	! =====================
[3792]	1166
[5105]	1167	! **********
	1168	CALL SISVAT_qSo
	1169	! #m0. (Wats_0,Wats_1,Wats_d)
	1170	! **********
[3792]	1171
	1172
[5105]	1173	! Surface Fluxes
	1174	! =====================
[3792]	1175
[5105]	1176	DO ikl=1,knonv
	1177	IRdwsv(ikl)=IRd_SV(ikl)*Eso_sv(ikl) ! Downward IR
	1178	! IRdwsv(ikl)=tau_sv(ikl) IRd_SV(ikl)Eso_sv(ikl) ! Downward IR
	1179	! . +(1.0-tau_sv(ikl))IRd_SV(ikl)Evg_sv(ikl) ! ! Etienne, remove vegetation component
	1180	IRupsv(ikl) = IRupsv(ikl) ! Upward IR
	1181	IRu_SV(ikl) = -IRupsv(ikl) & ! Upward IR
	1182	+IRd_SV(ikl) & ! (effective)
	1183	-IRdwsv(ikl) ! (positive)
[3792]	1184
[5105]	1185	TBr_sv(ikl) =sqrt(sqrt(IRu_SV(ikl)/StefBo)) ! Brightness
[5113]	1186	! Temperature
[5105]	1187	uts_SV(ikl) = (HSv_sv(ikl) +HSs_sv(ikl)) & ! uT
	1188	/(rhT_SV(ikl) *cp) !
	1189	uqs_SV(ikl) = (HLv_sv(ikl) +HLs_sv(ikl)) & ! uq
	1190	/(rhT_SV(ikl) *LhvH2O) !
	1191	LMO_SV(ikl) = TaT_SV(ikl)(us__SV(ikl)*3) &
	1192	/gravit/uts_SV(ikl)/vonKrm ! MO length
[3792]	1193
	1194
[5105]	1195	! Surface Temperature
	1196	! ^^^^^^^^^^^^^^^^^^^^
[3792]	1197
[5105]	1198	IF (iflag_tsurf_inlandsis == 0) THEN
[3792]	1199
[5105]	1200	Tsrfsv(ikl) =TsisSV(ikl,isnoSV(ikl))
[3792]	1201
[5105]	1202	ELSE IF (iflag_tsurf_inlandsis > 0) THEN
	1203	! Etienne: extrapolation from the two uppermost levels:
[3792]	1204
[5117]	1205	IF (isnoSV(ikl) >=2) THEN
[5105]	1206	zm1=-dzsnSV(ikl,isnoSV(ikl))/2.
	1207	zm2=-(dzsnSV(ikl,isnoSV(ikl)) + dzsnSV(ikl,isnoSV(ikl)-1)/2.)
[5117]	1208	ELSE IF (isnoSV(ikl) == 1) THEN
[5105]	1209	zm1=-dzsnSV(ikl,isnoSV(ikl))/2.
	1210	zm2=-(dzsnSV(ikl,isnoSV(ikl))+dz_dSV(0)/2.)
	1211	else
	1212	zm1=-dz_dSV(0)/2.
	1213	zm2=-(dz_dSV(0)+dz_dSV(-1)/2.)
[3792]	1214
[5105]	1215	end if
[3792]	1216
[5105]	1217	coefslope=(TsisSV(ikl,isnoSV(ikl))-TsisSV(ikl,isnoSV(ikl)-1)) &
	1218	/(zm1-zm2)
	1219	Tsrfsv(ikl)=TsisSV(ikl,isnoSV(ikl))+coefslope*(0. - zm1)
[3792]	1220
	1221
[5105]	1222	ELSE !(default)
[3900]	1223
[5105]	1224	Tsrfsv(ikl) =TsisSV(ikl,isnoSV(ikl))
[3900]	1225
[5105]	1226	END IF
[3792]	1227
	1228
[5105]	1229	END DO
[3792]	1230
[5105]	1231	! Snow Pack Properties (sphericity, dendricity, size)
	1232	! ===================================================
[3792]	1233
[5105]	1234	IF (SnoMod) THEN
[3792]	1235
[5117]	1236	IF (discret_xf .AND. klonv==1) THEN
[5116]	1237	IF(isnoSV(1)>=1) THEN
[5105]	1238	! + **********
	1239	CALL SISVAT_GSn
	1240	! + **********
[5117]	1241	ENDIF
[5105]	1242	else
	1243	! + **********
	1244	CALL SISVAT_GSn
	1245	! + **********
[5117]	1246	ENDIF
[3792]	1247
	1248
[5105]	1249	END IF
[3900]	1250
	1251
[5105]	1252	! Roughness Length for next time step
	1253	!====================================
[3900]	1254
[5105]	1255	! Note that in INLANDSIS, we treat only ice covered surfaces so calculation
	1256	! of z0 is much simpler (no subgrid fraction of ocean or land)
	1257	! old calculations are commented below
[3792]	1258
	1259
[5105]	1260	! +--Roughness Length for Momentum
	1261	! + -----------------------------
	1262
	1263	! ETIENNE WARNING: changes have been made wrt original SISVAT
	1264
	1265	! +--Land+Sea-Ice / Ice-free Sea Mask
	1266	! + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	1267	DO ikl=1,knonv
	1268	IcIndx(ikl) = 0
	1269	ENDDO
	1270	DO isn=1,nsno
	1271	DO ikl=1,knonv
	1272
	1273	IcIndx(ikl) = max(IcIndx(ikl), &
	1274	isn*max(0, &
	1275	sign(1, &
	1276	int(ro__SV(ikl,isn)-900.))))
	1277	ENDDO
	1278	ENDDO
	1279
	1280	DO ikl=1,knonv
	1281	LISmsk = 1. ! in inlandsis, land only
	1282	IceMsk = max(0,sign(1 ,IcIndx(ikl)-1) )
	1283	SnoMsk = max(min(isnoSV(ikl)-iiceSV(ikl),1),0)
	1284
	1285
	1286	! +--Z0 Smooth Regime over Snow (Andreas 1995, CRREL Report 95-16, p. 8)
	1287	! + ^^^^^^^^^^^^^^^^^^^^^^^^^^
	1288	Z0m_nu = 5.e-5 ! z0s~(10-d)*exp(-vonkar/sqrt(1.1e-03))
	1289
	1290	! +--Z0 Saltat.Regime over Snow (Gallee et al., 2001, BLM 99 (19) p.11)
	1291	! + ^^^^^^^^^^^^^^^^^^^^^^^^^^
	1292
	1293	u2star = us__SV(ikl) *us__SV(ikl)
	1294	Z0mBSn = u2star *0.536e-3 - 61.8e-6
	1295	Z0mBSn = max(Z0mBS0 ,Z0mBSn)
	1296
	1297	! +--Z0 Smooth + Saltat. Regime
	1298	! + ^^^^^^^^^^^^^^^^^^^^^^^^^^
	1299	Z0enSV(ikl) = Z0m_nu &
	1300	+ Z0mBSn
	1301
	1302
	1303	! Calculation of snow roughness length
	1304	!=====================================
	1305	IF (iflag_z0m_snow == 0) THEN
	1306
	1307	Z0m_Sn=prescribed_z0m_snow
	1308
	1309	ELSE IF (iflag_z0m_snow == 1) THEN
	1310
	1311	Z0m_Sn=Z0enSV(ikl)
	1312
	1313	ELSE IF (iflag_z0m_snow == 2) THEN
	1314
	1315	! +--Rough Snow Surface Roughness Length (Variable Sastrugi Height)
	1316	! + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	1317	A_Fact = 1.0000 ! Andreas et al., 2004, p.4
[5113]	1318	! ams.confex.com/ams/pdfpapers/68601.pdf
[5105]	1319
	1320	! Parameterization of z0 dependance on Temperature (C. Amory, 2017)
	1321	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1322	! Z0=f(T) deduced from observations, Adelie Land, dec2012-dec2013
	1323
	1324
	1325	coefa = 0.1658 !0.1862 !Ant
	1326	coefb = -50.3869 !-55.7718 !Ant
	1327	ta1 = 253.15 !255. Ant
	1328	ta2 = 273.15
	1329	ta3 = 273.15+3
	1330	z01 = exp(coefa*ta1 + coefb) !~0.2 ! ~0.25 mm
	1331	z02 = exp(coefa*ta2 + coefb) !~6 !~7 mm
	1332	z03 = z01
	1333	coefc = log(z03/z02)/(ta3-ta2)
	1334	coefd = log(z03)-coefc*ta3
	1335
[5117]	1336	IF (TaT_SV(ikl) < ta1) THEN
[5105]	1337	Z0_obs = z01
[5117]	1338	ELSE IF (TaT_SV(ikl)>=ta1 .AND. TaT_SV(ikl)<ta2) THEN
[5105]	1339	Z0_obs = exp(coefa*TaT_SV(ikl) + coefb)
[5117]	1340	ELSE IF (TaT_SV(ikl)>=ta2 .AND. TaT_SV(ikl)<ta3) THEN
[5113]	1341	! if st > 0, melting induce smooth surface
[5105]	1342	Z0_obs = exp(coefc*TaT_SV(ikl) + coefd)
[3900]	1343	else
[5105]	1344	Z0_obs = z03
[3900]	1345	endif
[3792]	1346
[5105]	1347	Z0m_Sn=Z0_obs
[3792]	1348
	1349
[5105]	1350	ELSE
[3792]	1351
[5105]	1352	Z0m_Sn=0.500e-3 ! default=0.500e-3m (tuning of MAR)
[3792]	1353
[5105]	1354	ENDIF
[3792]	1355
	1356
[3900]	1357
[5105]	1358	! param = Z0_obs/1. ! param(s) \| 1.(m/s)=TUNING
	1359	! #SZ Z0Sa_N = (us__SV(ikl) -0.2)*param ! 0.0001=TUNING
	1360	! #SZ. * max(zero,sign(unun,TfSnow-eps9
	1361	! #SZ. -TsisSV(ikl , isnoSV(ikl))))
	1362	!!#SZ Z0SaSi = max(zero,sign(unun,Z0Sa_N ))! 1 if erosion
	1363	! #SZ Z0SaSi = max(zero,sign(unun,zero -eps9 -uss_SV(ikl)))!
	1364	! #SZ Z0Sa_N = max(zero, Z0Sa_N)
	1365	! #SZ Z0SaSV(ikl) =
	1366	! #SZ. max(Z0SaSV(ikl) ,Z0SaSV(ikl)
	1367	! #SZ. + Z0SaSi(Z0Sa_N-Z0SaSV(ikl))exp(-dt__SV/43200.))
	1368	! #SZ. - min(dz0_SV(ikl) , Z0SaSV(ikl))
[3900]	1369
[5105]	1370	! #SZ A_Fact = Z0SaSV(ikl) * 5.0/0.15 ! A=5 if h~10cm
	1371	! +... CAUTION: The influence of the sastrugi direction is not yet included
[3792]	1372
[5105]	1373	! #SZ Z0m_Sn = Z0SaSV(ikl) !
	1374	! #SZ. - Z0m_nu !
[3792]	1375
[5105]	1376	! +--Z0 Saltat.Regime over Snow (Shao & Lin, 1999, BLM 91 (46) p.222)
	1377	! + ^^^^^^^^^^^^^^^^^^^^^^^^^^
	1378	! #ZN sqrrZ0 = usthSV(ikl)/max( us__SV(ikl),0.001)
	1379	! #ZN sqrrZ0 = min( sqrrZ0 ,0.999)
	1380	! #ZN Z0mBSn = 0.55 0.55 exp(-sqrrZ0 *sqrrZ0)
	1381	! #ZN. us__SV(ikl) us__SV(ikl)grvinv0.5
[3900]	1382
[5105]	1383	! +--Z0 Smooth + Saltat. Regime (Shao & Lin, 1999, BLM 91 (46) p.222)
	1384	! + ^^^^^^^^^^^^^^^^^^^^^^^^^^
	1385	! #ZN Z0enSV(ikl) = (Z0m_nu ** sqrrZ0 )
	1386	! #ZN. * (Z0mBSn **(1.-sqrrZ0))
	1387	! #ZN Z0enSV(ikl) = max(Z0enSV(ikl), Z0m_nu)
[3900]	1388
	1389
[5105]	1390	! +--Z0 Smooth Regime over Snow (Andreas etAl., 2004
	1391	! + ^^^^^^^^^^^^^^^^^^^^^^^^^^ ams.confex.com/ams/pdfpapers/68601.pdf)
	1392	! #ZA Z0m_nu = 0.135*akmol / max(us__SV(ikl) , epsi)
[3900]	1393
[5105]	1394	! +--Z0 Saltat.Regime over Snow (Andreas etAl., 2004
	1395	! + ^^^^^^^^^^^^^^^^^^^^^^^^^^ ams.confex.com/ams/pdfpapers/68601.pdf)
	1396	! #ZA Z0mBSn = 0.035u2star grvinv
[3900]	1397
[5105]	1398	! +--Z0 Smooth + Saltat. Regime (Andreas etAl., 2004
	1399	! ( used by Erosion) ams.confex.com/ams/pdfpapers/68601.pdf)
	1400	! ^^^^^^^^^^^^^^^^^^^^^^^^^^
	1401	! #ZA Z0enSV(ikl) = Z0m_nu
	1402	! #ZA. + Z0mBSn
[3900]	1403
[5105]	1404	! +--Z0 Rough Regime over Snow (Andreas etAl., 2004
	1405	! + (.NOT. used by Erosion) ams.confex.com/ams/pdfpapers/68601.pdf)
	1406	! ^^^^^^^^^^^^^^^^^^^^^^^^^^
	1407	!!#ZA u2star = (us__SV(ikl) -0.1800) / 0.1
	1408	!!#ZA Z0m_Sn =A_FactZ0mBSn exp(-u2star*u2star)
	1409	! #ZA Z0m_90 =(10.-0.025*VVs_SV(ikl)/5.)
	1410	! #ZA. exp(-0.4/sqrt(.00275+.00001max(0.,VVs_SV(ikl)-5.)))
	1411	! #ZA Z0m_Sn = DDs_SV(ikl)* Z0m_90 / 45.
	1412	! #ZA. - DDs_SV(ikl)DDs_SV(ikl) Z0m_90 /(90.*90.)
[3900]	1413
	1414
	1415
	1416
[5105]	1417	! +--Z0 (Erosion) over Snow (instantaneous)
	1418	! + ^^^^^^^^^^^^^^^^^^^^^^^^^^
	1419	Z0e_SV(ikl) = Z0enSV(ikl)
[3900]	1420
[5105]	1421	! +--Momentum Roughness Length (Etienne: changes wrt original SISVAT)
	1422	! + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	1423	Z0mnSV(ikl) = Z0m_nu *(1-SnoMsk) & ! Ice z0
	1424	+ (Z0m_Sn)*SnoMsk ! Snow Sastrugi Form and Snow Erosion
[3792]	1425
[3900]	1426
[5105]	1427	! +--GIS Roughness Length
	1428	! + ^^^^^^^^^^^^^^^^^^^^^
	1429	! #GL Z0mnSV(ikl) =
	1430	! #GL. (1-LSmask(ikl)) * Z0mnSV(ikl)
	1431	! #GL. + LSmask(ikl) * max(Z0mnSV(ikl),max(Z0_GIM,
	1432	! #GL. Z0_GIM+
	1433	! #GL. (0.0032-Z0_GIM)*(ro__SV(ikl,isnoSV(ikl))-600.) !
	1434	! #GL. /(920.00 -600.))) !
[3900]	1435
[5105]	1436	! +--Mom. Roughness Length, Instantaneous
	1437	! + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	1438	Z0m_SV(ikl) = Z0mnSV(ikl) ! Z0mnSV instant.
[3900]	1439
	1440
[5105]	1441	! +--Roughness Length for Scalars
	1442	! + ----------------------------
[3900]	1443
[5105]	1444	Z0hnSV(ikl) = Z0mnSV(ikl)/ 7.4
[3900]	1445
[5105]	1446	IF (is_ok_z0h_rn) THEN
[3900]	1447
[5105]	1448	rstar = Z0mnSV(ikl) * us__SV(ikl) / akmol
	1449	rstar = max(epsi,min(rstar,R_1000))
	1450	alors = log(rstar)
	1451	rstar0 = 1.250e0 * max(zero,sign(unun,0.135e0 - rstar)) &
	1452	+(1. - max(zero,sign(unun,0.135e0 - rstar))) &
	1453	(0.149e0 max(zero,sign(unun,2.500e0 - rstar)) &
	1454	+ 0.317e0 &
	1455	*(1. - max(zero,sign(unun,2.500e0 - rstar))))
	1456	rstar1 = 0. * max(zero,sign(unun,0.135e0 - rstar)) &
	1457	+(1. - max(zero,sign(unun,0.135e0 - rstar))) &
	1458	(-0.55e0 max(zero,sign(unun,2.500e0 - rstar)) &
	1459	- 0.565 &
	1460	*(1. - max(zero,sign(unun,2.500e0 - rstar))))
	1461	rstar2 = 0. * max(zero,sign(unun,0.135e0 - rstar)) &
	1462	+(1. - max(zero,sign(unun,0.135e0 - rstar))) &
	1463	(0. max(zero,sign(unun,2.500e0 - rstar)) &
	1464	- 0.183 &
	1465	*(unun - max(zero,sign(unun,2.500e0 - rstar))))
[3900]	1466
	1467
	1468
[5105]	1469	!XF #RN (is_ok_z0h_rn) does not work well over bare ice
	1470	!XF MAR is then too warm and not enough melt
[3792]	1471
[5116]	1472	IF(ro__SV(ikl,isnoSV(ikl))>50 &
[5117]	1473	.AND.ro__SV(ikl,isnoSV(ikl))<roSdSV)THEN
[5105]	1474	Z0hnSV(ikl) = max(zero &
	1475	, sign(unun,zzsnsv(ikl,isnoSV(ikl))-epsi)) &
	1476	* exp(rstar0+rstar1alors+rstar2alors*alors) &
	1477	* 0.001e0 + Z0hnSV(ikl) * ( 1. - max(zero &
	1478	, sign(unun,zzsnsv(ikl,isnoSV(ikl))-epsi)))
[3900]	1479
[5105]	1480	endif
[3900]	1481
	1482
[5105]	1483	ENDIF
[3792]	1484
[5105]	1485	Z0h_SV(ikl) = Z0hnSV(ikl)
[3792]	1486
	1487
[5105]	1488	! #MT Z0m_SV(ikl) = max(2.0e-6 ,Z0m_SV(ikl)) ! Min Z0_m (Garrat Scheme)
	1489	! Z0m_SV(ikl) = min(Z0m_SV(ikl),za__SV(ikl)*0.3333)
[3792]	1490
	1491
[5105]	1492	END DO
[3792]	1493
	1494
	1495
[5116]	1496	END SUBROUTINE inlandsis
[3792]	1497
	1498
	1499
	1500
	1501
[5105]	1502
	1503
	1504
	1505
	1506
	1507

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

Download in other formats:

Original Format