Context Navigation

inlandsis.F @ 3819

Last change on this file since 3819 was 3792, checked in by evignon, 4 years ago
Ajout de INLANDSIS, nouvelle interface entre LMDZ et la neige de SISVAT Etienne, 04/01/2021
File size: 54.6 KB

Rev	Line
[3792]	1	subroutine INLANDSIS(SnoMod,BloMod,jjtime)
	2
	3	USE dimphy
	4
	5	!--------------------------------------------------------------------------+
	6	! INLANDSIS module \|
	7	! Simplified SISVAT module, containing ice and snow processes for \|
	8	! ice-covered surfaces \|
	9	! version MARv3, november 2020 \|
	10	! SubRoutine INLANDSIS contains the fortran 77 code of the \|
	11	! Soil/Ice Snow Vegetation Atmosphere Transfer Scheme \|
	12	! \|
	13	!--------------------------------------------------------------------------+
	14	! PARAMETERS: klonv: Total Number of columns = \|
	15	! ^^^^^^^^^^ = Total Number of continental grid boxes \|
	16	! X Number of Mosaic Cell per grid box \|
	17	! \|
	18	! INPUT: daHost : Date Host Model \|
	19	! ^^^^^ \|
	20	! \|
	21	! INPUT: LSmask : 1: Land MASK \|
	22	! ^^^^^ 0: Sea MASK \|
	23	! isotSV = 0,...,12: Soil Type \|
	24	! 0: Water, Liquid (Sea, Lake) \|
	25	! 12: Water, Solid (Ice) \|
	26	! \|
	27	! INPUT: coszSV : Cosine of the Sun Zenithal Distance [-] \|
	28	! ^^^^^ sol_SV : Surface Downward Solar Radiation [W/m2] \|
	29	! IRd_SV : Surface Downward Longwave Radiation [W/m2] \|
	30	! drr_SV : Rain Intensity [kg/m2/s] \|
	31	! dsn_SV : Snow Intensity [mm w.e./s] \|
	32	! dsnbSV : Snow Intensity, Drift Fraction [-] \|
	33	! dbs_SV : Drift Amount [mm w.e.] \|
	34	! za__SV : Surface Boundary Layer (SBL) Height [m] \|
	35	! VV__SV :(SBL Top) Wind Velocity [m/s] \|
	36	! TaT_SV : SBL Top Temperature [K] \|
	37	! rhT_SV : SBL Top Air Density [kg/m3] \|
	38	! QaT_SV : SBL Top Specific Humidity [kg/kg] \|
	39	! qsnoSV : SBL Mean Snow Content [kg/kg] \|
	40	! alb0SV : Soil Basic Albedo [-] \|
	41	! slopSV : Surface Slope [-] \|
	42	! dt__SV : Time Step [s] \|
	43	! \|
	44	! INPUT / isnoSV = total Nb of Ice/Snow Layers \|
	45	! OUTPUT: ispiSV = 0,...,nsno: Uppermost Superimposed Ice Layer \|
	46	! ^^^^^^ iiceSV = total Nb of Ice Layers \|
	47	! istoSV = 0,...,5 : Snow History (see istdSV data) \|
	48	! \|
	49	! INPUT / alb_SV : Surface Albedo [-] \|
	50	! OUTPUT: emi_SV : Surface Emissivity [-] \|
	51	! ^^^^^^ IRs_SV : Soil IR Flux (negative) [W/m2] \|
	52	! LMO_SV : Monin-Obukhov Scale [m] \|
	53	! us__SV : Friction Velocity [m/s] \|
	54	! uts_SV : Temperature Turbulent Scale [m/s] \|
	55	! uqs_SV : Specific Humidity Velocity [m/s] \|
	56	! uss_SV : Blowing Snow Turbulent Scale [m/s] \|
	57	! usthSV : Blowing Snow Erosion Threshold [m/s] \|
	58	! Z0m_SV : Momentum Roughness Length [m] \|
	59	! Z0mmSV : Momentum Roughness Length (time mean) [m] \|
	60	! Z0mnSV : Momentum Roughness Length (instantaneous)[m] \|
	61	! Z0SaSV : Sastrugi Roughness Length [m] \|
	62	! Z0e_SV : Erosion Snow Roughness Length [m] \|
	63	! Z0emSV : Erosion Snow Roughness Length (time mean) [m] \|
	64	! Z0enSV : Erosion Snow Roughness Length (instantaneous)[m] \|
	65	! Z0roSV : Subgrid Topo Roughness Length [m] \|
	66	! Z0h_SV : Heat Roughness Length [m] \|
	67	! TsisSV : Soil/Ice Temperatures (layers -nsol,-nsol+1,..,0)\|
	68	! & Snow Temperatures (layers 1,2,...,nsno) [K] \|
	69	! ro__SV : Soil/Snow Volumic Mass [kg/m3] \|
	70	! eta_SV : Soil/Snow Water Content [m3/m3] \|
	71	! G1snSV : snow dendricity/sphericity \|
	72	! G2snSV : snow sphericity/grain size \|
	73	! dzsnSV : Snow Layer Thickness [m] \|
	74	! agsnSV : Snow Age [day] \|
	75	! BufsSV : Snow Buffer Layer [kg/m2] .OR. [mm] \|
	76	! BrosSV : Snow Buffer Layer Density [kg/m3] \|
	77	! BG1sSV : Snow Buffer Layer Dendricity / Sphericity [-] \|
	78	! BG2sSV : Snow Buffer Layer Sphericity / Size [-] [0.1 mm] \|
	79	! rusnSV : Surficial Water [kg/m2] .OR. [mm] \|
	80	! \|
	81	! OUTPUT: no__SV : OUTPUT file Unit Number [-] \|
	82	! ^^^^^^ i___SV : OUTPUT point i Coordinate [-] \|
	83	! j___SV : OUTPUT point j Coordinate [-] \|
	84	! n___SV : OUTPUT point n Coordinate [-] \|
	85	! lwriSV : OUTPUT point vec Index [-] \|
	86	! \|
	87	! OUTPUT: IRu_SV : Upward IR Flux (+, upw., effective) [K] \|
	88	! ^^^^^^ hSalSV : Saltating Layer Height [m] \|
	89	! qSalSV : Saltating Snow Concentration [kg/kg] \|
	90	! RnofSV : RunOFF Intensity [kg/m2/s] \|
	91	! \|
	92	! Internal Variables: \|
	93	! ^^^^^^^^^^^^^^^^^^ \|
	94	! NLaysv = New Snow Layer Switch [-] \|
	95	! albisv : Snow/Ice/Water/Soil Integrated Albedo [-] \|
	96	! SoSosv : Absorbed Solar Radiation by Surfac.(Normaliz)[-] \|
	97	! TBr_sv : Brightness Temperature [K] \|
	98	! IRupsv : Upward IR Flux (-, upw.) [W/m2] \|
	99	! ram_sv : Aerodynamic Resistance for Momentum [s/m] \|
	100	! rah_sv : Aerodynamic Resistance for Heat [s/m] \|
	101	! Evp_sv : Evaporation [kg/m2] \|
	102	! EvT_sv : Evapotranspiration [kg/m2] \|
	103	! HSs_sv : Surface Sensible Heat Flux + => absorb.[W/m2] \|
	104	! HLs_sv : Surface Latent Heat Flux + => absorb.[W/m2] \|
	105	! Lx_H2O : Latent Heat of Vaporization/Sublimation [J/kg] \|
	106	! Tsrfsv : Surface Temperature [K] \|
	107	! sEX_sv : Verticaly Integrated Extinction Coefficient [-] \|
	108	! LSdzsv : Vertical Discretization Factor [-] \|
	109	! = 1. Soil \|
	110	! = 1000. Ocean \|
	111	! z_snsv : Snow Pack Thickness [m] \|
	112	! zzsnsv : Snow Pack Thickness [m] \|
	113	! albssv : Soil Albedo [-] \|
	114	! Eso_sv : Soil+Snow Emissivity [-] \|
	115	! Khydsv : Soil Hydraulic Conductivity [m/s] \|
	116	! \|
	117	! ETSo_0 : Snow/Soil Energy Power, before Forcing [W/m2] \|
	118	! ETSo_1 : Snow/Soil Energy Power, after Forcing [W/m2] \|
	119	! ETSo_d : Snow/Soil Energy Power Forcing [W/m2] \|
	120	! EqSn_0 : Snow Energy, before Phase Change [J/m2] \|
	121	! EqSn_1 : Snow Energy, after Phase Change [J/m2] \|
	122	! EqSn_d : Snow Energy, net Forcing [J/m2] \|
	123	! Enrsvd : SVAT Energy Power Forcing [W/m2] \|
	124	! Enrbal : SVAT Energy Balance [W/m2] \|
	125	! Wats_0 : Soil Water, before Forcing [mm] \|
	126	! Wats_1 : Soil Water, after Forcing [mm] \|
	127	! Wats_d : Soil Water Forcing [mm] \|
	128	! SIWm_0 : Snow initial Mass [mm w.e.] \|
	129	! SIWm_1 : Snow final Mass [mm w.e.] \|
	130	! SIWa_i : Snow Atmos. initial Forcing [mm w.e.] \|
	131	! SIWa_f : Snow Atmos. final Forcing(noConsumed)[mm w.e.] \|
	132	! SIWe_i : SnowErosion initial Forcing [mm w.e.] \|
	133	! SIWe_f : SnowErosion final Forcing(noConsumed)[mm w.e.] \|
	134	! SIsubl : Snow sublimed/deposed Mass [mm w.e.] \|
	135	! SImelt : Snow Melted Mass [mm w.e.] \|
	136	! SIrnof : Surficial Water + Run OFF Change [mm w.e.] \|
	137	! SIvAcr : Sea-Ice vertical Acretion [mm w.e.] \|
	138	! Watsvd : SVAT Water Forcing [mm] \|
	139	! Watbal : SVAT Water Balance [W/m2] \|
	140	! \|
	141	! vk2 : Square of Von Karman Constant [-] \|
	142	! sqrCm0 : Factor of Neutral Drag Coeffic.Momentum [s/m] \|
	143	! sqrCh0 : Factor of Neutral Drag Coeffic.Heat [s/m] \|
	144	! EmiSol : Soil Emissivity [-] \|
	145	! EmiSno : Snow Emissivity [-] \|
	146	! EmiWat : Water Emissivity [-] \|
	147	! Z0mLnd : Land Roughness Length [m] \|
	148	! sqrrZ0 : ut/u \|
	149	! f_eff : Marticorena & B. 1995 JGR (20) \|
	150	! A_Fact : Fundamental * Roughness \|
	151	! Z0mBSn : BSnow Roughness Length [m] \|
	152	! Z0mBS0 : Mimimum BSnow Roughness Length (blown* ) [m] \|
	153	! Z0m_Sn : Snow Roughness Length (surface) [m] \|
	154	! Z0m_S0 : Mimimum Snow Roughness Length [m] \|
	155	! Z0m_S1 : Maximum Snow Roughness Length [m] \|
	156	! Z0_GIM : Minimum GIMEX Roughness Length [m] \|
	157	! Z0_ICE : Sea Ice ISW Roughness Length [m] \|
	158	! \|
	159	! \|
	160	!--------------------------------------------------------------------------+
	161
	162
	163
	164	! Global Variables
	165	! ================
	166
	167
	168	USE VARphy
	169	USE VAR_SV
	170	USE VARdSV
	171	USE VAR0SV
	172	USE VARxSV
	173	USE VARySV
	174	USE VARtSV
	175	USE surface_data, only: iflag_tsurf_inlandsis
	176
	177
	178	IMPLICIT NONE
	179
	180	logical SnoMod
	181	logical BloMod
	182	integer jjtime
	183
	184
	185	! Internal Variables
	186	! ==================
	187
	188	! Non Local
	189	! ---------
	190
	191	real TBr_sv(klonv) ! Brightness Temperature
	192	real IRdwsv(klonv) ! DOWNward IR Flux
	193	real IRupsv(klonv) ! UPward IR Flux
	194	real d_Bufs,Bufs_N ! Buffer Snow Layer Increment
	195	real Buf_ro,Bros_N ! Buffer Snow Layer Density
	196	real BufPro ! Buffer Snow Layer Density
	197	real Buf_G1,BG1__N ! Buffer Snow Layer Dendr/Sphe[-]
	198	real Buf_G2,BG2__N ! Buffer Snow Layer Spher/Size[-]
	199	real Bdzssv(klonv) ! Buffer Snow Layer Thickness
	200	real z_snsv(klonv) ! Snow-Ice, current Thickness
	201
	202
	203
	204	! Local
	205	! -----
	206
	207	integer iwr
	208	integer ikl ,isn ,isl ,ist !
	209	integer ist__s,ist__w ! Soil/Water Body Identifier
	210	integer growth ! Seasonal Mask
	211	integer LISmsk ! Land+Ice / Open Sea Mask
	212	integer LSnMsk ! Snow-Ice / No Snow-Ice Mask
	213	integer IceMsk,IcIndx(klonv) ! Ice / No Ice Mask
	214	integer SnoMsk ! Snow / No Snow Mask
	215
	216	real roSMin,roSMax,roSn_1,roSn_2,roSn_3 ! Fallen Snow Density (PAHAUT)
	217	real Dendr1,Dendr2,Dendr3 ! Fallen Snow Dendric.(GIRAUD)
	218	real Spher1,Spher2,Spher3,Spher4 ! Fallen Snow Spheric.(GIRAUD)
	219	real Polair ! Polar Snow Switch
	220	real PorSno,Por_BS,Salt_f,PorRef !
	221	c #sw real PorVol,rWater !
	222	c #sw real rusNEW,rdzNEW,etaNEW !
	223	real ro_new !
	224	real TaPole ! Maximum Polar Temperature
	225	real T__Min ! Minimum realistic Temperature
	226	real EmiSol ! Emissivity of Soil
	227	real EmiSno ! Emissivity of Snow
	228	real EmiWat ! Emissivity of a Water Area
	229	real vk2 ! Square of Von Karman Constant
	230	real u2star !(u)*2
	231	real Z0mLnd ! Land Roughness Length
	232	c #ZN real sqrrZ0 ! ut/u
	233	real f_eff ! Marticorena & B. 1995 JGR (20)
	234	real A_Fact ! Fundamental * Roughness
	235	real Z0m_nu ! Smooth R Snow Roughness Length
	236	real Z0mBSn ! BSnow Roughness Length
	237	real Z0mBS0 ! Mimimum BSnow Roughness Length
	238	real Z0m_S0 ! Mimimum Snow Roughness Length
	239	real Z0m_S1 ! Maximum Snow Roughness Length
	240	c #SZ real Z0Sa_N ! Regime Snow Roughness Length
	241	c #SZ real Z0SaSi ! 1.IF Rgm Snow Roughness Length
	242	c #GL real Z0_GIM ! Mimimum GIMEX Roughness Length
	243	real Z0_ICE ! Ice ISW Roughness Length
	244	real Z0m_Sn,Z0m_90 ! Snow Surface Roughness Length
	245	real SnoWat ! Snow Layer Switch
	246	c #RN real rstar,alors !
	247	c #RN real rstar0,rstar1,rstar2 !
	248	real SameOK ! 1. => Same Type of Grains
	249	real G1same ! Averaged G1, same Grains
	250	real G2same ! Averaged G2, same Grains
	251	real typ__1 ! 1. => Lay1 Type: Dendritic
	252	real zroNEW ! dz X ro, if fresh Snow
	253	real G1_NEW ! G1, if fresh Snow
	254	real G2_NEW ! G2, if fresh Snow
	255	real zroOLD ! dz X ro, if old Snow
	256	real G1_OLD ! G1, if old Snow
	257	real G2_OLD ! G2, if old Snow
	258	real SizNEW ! Size, if fresh Snow
	259	real SphNEW ! Spheric.,if fresh Snow
	260	real SizOLD ! Size, if old Snow
	261	real SphOLD ! Spheric.,if old Snow
	262	real Siz_av ! Averaged Grain Size
	263	real Sph_av ! Averaged Grain Spher.
	264	real Den_av ! Averaged Grain Dendr.
	265	real DendOK ! 1. => Average is Dendr.
	266	real G1diff ! Averaged G1, diff. Grains
	267	real G2diff ! Averaged G2, diff. Grains
	268	real G1 ! Averaged G1
	269	real G2 ! Averaged G2
	270	real param ! Polynomial fit z0=f(T)
	271	real Z0_obs ! Fit Z0_obs=f(T) (m)
	272	real tamin ! min T of linear fit (K)
	273	real tamax ! max T of linear fit (K)
	274	real coefa,coefb,coefc,coefd ! Coefs for z0=f(T)
	275	real ta1,ta2,ta3 ! Air temperature thresholds
	276	real z01,z02,z03 ! z0 thresholds
	277	real tt_c,vv_c ! Critical param.
	278	real tt_tmp,vv_tmp,vv_virt ! Temporary variables
	279	logical density_kotlyakov ! .true. if Kotlyakov 1961
	280	real e_prad,e1pRad,A_Rad0,absg_V,absgnI,exdRad ! variables for SoSosv calculations
	281	real zm1, zm2, coefslope ! variables for surface temperature extrapolation
	282
	283
	284
	285	! Internal DATA
	286	! =============
	287
	288	data T__Min / 200.00/ ! Minimum realistic Temperature
	289	data TaPole / 263.15/ ! Maximum Polar Temperature
	290	data roSMin / 30. / ! Minimum Snow Density
	291	data roSMax / 400. / ! Max Fresh Snow Density
	292	data tt_c / -2.0 / ! Critical Temp. (degC)
	293	data vv_c / 14.3 / ! Critical Wind speed (m/s)
	294	data roSn_1 / 109. / ! Fall.Sno.Density, Indep. Param.
	295	data roSn_2 / 6. / ! Fall.Sno.Density, Temper.Param.
	296	data roSn_3 / 26. / ! Fall.Sno.Density, Wind Param.
	297	data Dendr1 / 17.12/ ! Fall.Sno.Dendric.,Wind 1/Param.
	298	data Dendr2 / 128. / ! Fall.Sno.Dendric.,Wind 2/Param.
	299	data Dendr3 / -20. / ! Fall.Sno.Dendric.,Indep. Param.
	300	data Spher1 / 7.87/ ! Fall.Sno.Spheric.,Wind 1/Param.
	301	data Spher2 / 38. / ! Fall.Sno.Spheric.,Wind 2/Param.
	302	data Spher3 / 50. / ! Fall.Sno.Spheric.,Wind 3/Param.
	303	data Spher4 / 90. / ! Fall.Sno.Spheric.,Indep. Param.
	304	data EmiSol / 0.99999999/ ! 0.94Emissivity of Soil
	305	data EmiWat / 0.99999999/ ! Emissivity of a Water Area
	306	data EmiSno / 0.99999999/ ! Emissivity of Snow
	307
	308	! DATA Emissivities ! Pielke, 1984, pp. 383,409
	309
	310	data Z0mBS0 / 0.5e-6/ ! MINimum Snow Roughness Length
	311	! for Momentum if Blowing Snow
	312	! Gallee et al. 2001 BLM 99 (19)
	313	data Z0m_S0/ 0.00005/ ! MINimum Snow Roughness Length
	314	! MegaDunes included
	315	data Z0m_S1/ 0.030 / ! MAXimum Snow Roughness Length
	316	! (Sastrugis)
	317	c #GL data Z0_GIM/ 0.0013/ ! Ice Min Z0 = 0.0013 m (Broeke)
	318	! ! Old Ice Z0 = 0.0500 m (Bruce)
	319	! ! 0.0500 m (Smeets)
	320	! ! 0.1200 m (Broeke)
	321	data Z0_ICE/ 0.0010/ ! Sea-Ice Z0 = 0.0010 m (Andreas)
	322	! ! (Ice Station Weddel -- ISW)
	323	vk2 = vonKrm * vonKrm ! Square of Von Karman Constant
	324
	325
	326	! BEGIN.main.
	327	! ===========================
	328
	329
	330
	331
	332	! "Soil" Humidity of Water Bodies
	333	! ===============================
	334
	335	DO ikl=1,knonv
	336
	337	ist = isotSV(ikl) ! Soil Type
	338	ist__s = min(ist, 1) ! 1 => Soil
	339	ist__w = 1 - ist__s ! 1 => Water Body
	340	DO isl=-nsol,0
	341	eta_SV(ikl,isl) = eta_SV(ikl,isl) * ist__s ! Soil
	342	. + etadSV(ist) * ist__w ! Water Body
	343	END DO
	344
	345
	346	! Vertical Discretization Factor
	347	! ==============================
	348
	349	LSdzsv(ikl) = ist__s ! Soil
	350	. + OcndSV * ist__w ! Water Body
	351	END DO
	352
	353
	354	! Blowing Particles Threshold Friction velocity
	355	! =============================================
	356
	357	c #AE usthSV(ikl) = 1.0e+2
	358	! END DO
	359	!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
	360
	361
	362
	363
	364	! Contribution of Snow to the Surface Snow Pack
	365	! =============================================
	366
	367	IF (SnoMod) THEN
	368
	369
	370
	371	C +--Blowing Snow
	372	C + ------------
	373
	374	IF (BloMod) then
	375	if (klonv.eq.1) then
	376	if(isnoSV(1).ge.2 .and.
	377	. TsisSV(1,max(1,isnoSV(1)))<273. .and.
	378	. ro__SV(1,max(1,isnoSV(1)))<500. .and.
	379	. eta_SV(1,max(1,isnoSV(1)))<epsi) then
	380	C + **********
	381	call SISVAT_BSn
	382	endif
	383	else
	384	call SISVAT_BSn
	385	C + **********
	386	endif
	387	ENDIF
	388
	389
	390
	391	C +--Buffer Layer
	392	C + ------------
	393
	394	DO ikl=1,knonv
	395	c BufsSV(ikl) [mm w.e.] i.e, i.e., [kg/m2]
	396	d_Bufs = max(dsn_SV(ikl) *dt__SV,0.) !
	397	dsn_SV(ikl) = 0. !
	398	Bufs_N = BufsSV(ikl) +d_Bufs !
	399
	400
	401	C +--Snow Density
	402	C + ^^^^^^^^^^^^
	403	Polair = zero
	404	c #NP Polair = max(zero, !
	405	c #NP. sign(unun,TaPole !
	406	c #NP. -TaT_SV(ikl))) !
	407	Polair = max(zero, !
	408	. sign(unun,TaPole !
	409	. -TaT_SV(ikl))) !
	410	Buf_ro = max( rosMin, ! Fallen Snow Density
	411	. roSn_1+roSn_2* (TaT_SV(ikl)-TfSnow) ! [kg/m3]
	412	. +roSn_3*sqrt( VV__SV(ikl))) ! Pahaut (CEN), Etienne: use wind speed at first model level instead of 10m wind
	413	c #NP BufPro = max( rosMin, ! Fallen Snow Density
	414	c #NP. 104. *sqrt( max( VV10SV(ikl)-6.0,0.0))) ! Kotlyakov (1961)
	415
	416	density_kotlyakov = .true.
	417	c #AC density_kotlyakov = .false. !C.Agosta snow densisty as if BS is on b
	418	c #BS density_kotlyakov = .false. !C.Amory BS 2018
	419	C + ... Fallen Snow Density, Adapted for Antarctica
	420	if (density_kotlyakov) then
	421	tt_tmp = TaT_SV(ikl)-TfSnow
	422	!vv_tmp = VV10SV(ikl)
	423	vv_tmp=VV__SV(ikl) ! Etienne: use wind speed at first model level instead of 10m wind
	424	C + ... [ A compromise between
	425	C + ... Kotlyakov (1961) and Lenaerts (2012, JGR, Part1) ]
	426	if (tt_tmp.ge.-10) then
	427	BufPro = max( rosMin,
	428	. 104. *sqrt( max( vv_tmp-6.0,0.0))) ! Kotlyakov (1961)
	429	else
	430	vv_virt = (tt_cvv_tmp+vv_c(tt_tmp+10))
	431	. /(tt_c+tt_tmp+10)
	432	BufPro = 104. *sqrt( max( vv_virt-6.0,0.0))
	433	endif
	434	else
	435	C + ... [ density derived from observations of the first 50cm of
	436	C + ... snow - cf. Rajashree Datta - and multiplied by 0.8 ]
	437	C + ... C. Agosta, 2016-09
	438	cc #SD BufPro = 149.2 + 6.84VV10SV(ikl) + 0.48Tsrfsv(ikl)
	439	cc #SD BufPro = 125 + 14VV10SV(ikl) + 0.6Tsrfsv(ikl) !MAJ CK and CAm
	440	! BufPro = 200 + 21 * VV10SV(ikl)!CK 29/07/19
	441	BufPro = 200 + 21 * VV__SV(ikl)!Etienne: use wind speed at first model level instead of 10m wind
	442	endif
	443
	444	Bros_N = (1. - Polair) * Buf_ro ! Temperate Snow
	445	. + Polair * BufPro ! Polar Snow
	446
	447	Bros_N = max( 20.,max(rosMin, Bros_N))
	448	Bros_N = min(400.,min(rosMax-1,Bros_N)) ! for dz_min in SISVAT_zSn
	449
	450
	451	! Density of deposited blown snow
	452	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	453
	454	c #BS Bros_N = frsno
	455	c #BS ro_new = ro__SV(ikl,max(1,isnoSV(ikl)))
	456	c #BS ro_new = max(Bros_N,min(roBdSV,ro_new))
	457	c #BS Fac = 1-((ro__SV(ikl,max(1,isnoSV(ikl)))
	458	c #BS. -roBdSV)/(500.-roBdSV))
	459	c #BS Fac = max(0.,min(1.,Fac))
	460	c #BS dsnbSV(ikl) = Fac*dsnbSV(ikl)
	461	c #BS Bros_N = Bros_N * (1.0-dsnbSV(ikl))
	462	c #BS. + ro_new * dsnbSV(ikl)
	463
	464
	465
	466	! Time averaged Density of deposited blown Snow
	467	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	468
	469	BrosSV(ikl) =(Bros_N * d_Bufs !
	470	. +BrosSV(ikl)* BufsSV(ikl))!
	471	. / max(epsi,Bufs_N) !
	472
	473
	474	C +-- S.Falling Snow Properties (computed as in SISVAT_zAg)
	475	C + ^^^^^^^^^^^^^^^^^^^^^^^
	476	Buf_G1 = max(-G1_dSV, ! Temperate Snow
	477	. min(Dendr1*VV__SV(ikl)-Dendr2, ! Dendricity
	478	. Dendr3 )) !
	479	Buf_G2 = min( Spher4, ! Temperate Snow
	480	. max(Spher1*VV__SV(ikl)+Spher2, ! Sphericity
	481	. Spher3 )) !
	482	Buf_G1 = (1. - Polair) * Buf_G1 ! Temperate Snow
	483	. + Polair * G1_dSV ! Polar Snow
	484	Buf_G2 = (1. - Polair) * Buf_G2 ! Temperate Snow
	485	. + Polair * ADSdSV ! Polar Snow
	486	G1 = Buf_G1 ! NO Blown Snow
	487	G2 = Buf_G2 ! NO Blown Snow
	488
	489
	490	! S.1. Meme Type de Neige / same Grain Type
	491	! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	492	c #BS SameOK = max(zero,
	493	c #BS. sign(unun, Buf_G1 *G1_dSV
	494	c #BS. - eps_21 ))
	495	c #BS G1same = ((1.0-dsnbSV(ikl))Buf_G1+dsnbSV(ikl) G1_dSV)
	496	c #BS G2same = ((1.0-dsnbSV(ikl))Buf_G2+dsnbSV(ikl) ADSdSV)
	497	! Blowing Snow Properties: G1_dSV, ADSdSV
	498
	499	! S.2. Types differents / differents Types
	500	! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	501	c #BS typ__1 = max(zero,sign(unun,epsi-Buf_G1)) ! =1.=> Dendritic
	502	c #BS zroNEW = typ__1 *(1.0-dsnbSV(ikl)) ! fract.Dendr.Lay.
	503	c #BS. + (1.-typ__1) * dsnbSV(ikl) !
	504	c #BS G1_NEW = typ__1 *Buf_G1 ! G1 of Dendr.Lay.
	505	c #BS. + (1.-typ__1) *G1_dSV !
	506	c #BS G2_NEW = typ__1 *Buf_G2 ! G2 of Dendr.Lay.
	507	c #BS. + (1.-typ__1) *ADSdSV !
	508	c #BS zroOLD = (1.-typ__1) *(1.0-dsnbSV(ikl)) ! fract.Spher.Lay.
	509	c #BS. + typ__1 * dsnbSV(ikl) !
	510	c #BS G1_OLD = (1.-typ__1) *Buf_G1 ! G1 of Spher.Lay.
	511	c #BS. + typ__1 *G1_dSV !
	512	c #BS G2_OLD = (1.-typ__1) *Buf_G2 ! G2 of Spher.Lay.
	513	c #BS. + typ__1 *ADSdSV !
	514	c #BS SizNEW = -G1_NEW *DDcdSV/G1_dSV ! Size Dendr.Lay.
	515	c #BS. +(1.+G1_NEW /G1_dSV) !
	516	c #BS. (G2_NEW DScdSV/G1_dSV !
	517	c #BS. +(1.-G2_NEW /G1_dSV)*DFcdSV) !
	518	c #BS SphNEW = G2_NEW /G1_dSV ! Spher.Dendr.Lay.
	519	c #BS SizOLD = G2_OLD ! Size Spher.Lay.
	520	c #BS SphOLD = G1_OLD /G1_dSV ! Spher.Spher.Lay.
	521	c #BS Siz_av = (zroNEWSizNEW+zroOLDSizOLD) ! Averaged Size
	522	c #BS Sph_av = min( zroNEWSphNEW+zroOLDSphOLD !
	523	c #BS. , unun) ! Averaged Sphericity
	524	c #BS Den_av = min((Siz_av -( Sph_av *DScdSV !
	525	c #BS. +(1.-Sph_av)*DFcdSV)) !
	526	c #BS. / (DDcdSV -( Sph_av *DScdSV !
	527	c #BS. +(1.-Sph_av)*DFcdSV)) !
	528	c #BS. , unun) !
	529	c #BS DendOK = max(zero, !
	530	c #BS. sign(unun, Sph_av *DScdSV ! Small Grains
	531	c #BS. +(1.-Sph_av)*DFcdSV ! Faceted Grains
	532	c #BS. - Siz_av )) !
	533	C +... REMARQUE: le type moyen (dendritique ou non) depend
	534	C + ^^^^^^^^ de la comparaison avec le diametre optique
	535	C + d'une neige recente de dendricite nulle
	536	C +... REMARK: the mean type (dendritic or not) depends
	537	C + ^^^^^^ on the comparaison with the optical diameter
	538	C + of a recent snow having zero dendricity
	539
	540	c #BS G1diff =( -DendOK *Den_av
	541	c #BS. +(1.-DendOK)Sph_av) G1_dSV
	542	c #BS G2diff = DendOK Sph_av G1_dSV
	543	c #BS. +(1.-DendOK)*Siz_av
	544	c #BS G1 = SameOK *G1same
	545	c #BS. +(1.-SameOK)*G1diff
	546	c #BS G2 = SameOK *G2same
	547	c #BS. +(1.-SameOK)*G2diff
	548
	549
	550
	551	! S.1. Meme Type de Neige / same Grain Type
	552	! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	553	SameOK = max(zero,
	554	. sign(unun, Buf_G1 *BG1sSV(ikl)
	555	. - eps_21 ))
	556	G1same = (d_BufsBuf_G1+BufsSV(ikl)BG1sSV(ikl))
	557	. /max(epsi,Bufs_N)
	558	G2same = (d_BufsBuf_G2+BufsSV(ikl)BG2sSV(ikl))
	559	. /max(epsi,Bufs_N)
	560
	561	! S.2. Types differents / differents Types
	562	! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	563
	564	typ__1 = max(zero,sign(unun,epsi-Buf_G1)) ! =1.=> Dendritic
	565	zroNEW =( typ__1 *d_Bufs ! fract.Dendr.Lay.
	566	. + (1.-typ__1) *BufsSV(ikl)) !
	567	. /max(epsi,Bufs_N) !
	568	G1_NEW = typ__1 *Buf_G1 ! G1 of Dendr.Lay.
	569	. + (1.-typ__1) *BG1sSV(ikl) !
	570	G2_NEW = typ__1 *Buf_G2 ! G2 of Dendr.Lay.
	571	. + (1.-typ__1) *BG2sSV(ikl) !
	572	zroOLD =((1.-typ__1) *d_Bufs ! fract.Spher.Lay.
	573	. + typ__1 *BufsSV(ikl)) !
	574	. /max(epsi,Bufs_N) !
	575	G1_OLD = (1.-typ__1) *Buf_G1 ! G1 of Spher.Lay.
	576	. + typ__1 *BG1sSV(ikl) !
	577	G2_OLD = (1.-typ__1) *Buf_G2 ! G2 of Spher.Lay.
	578	. + typ__1 *BG2sSV(ikl) !
	579	SizNEW = -G1_NEW *DDcdSV/G1_dSV ! Size Dendr.Lay.
	580	. +(1.+G1_NEW /G1_dSV) !
	581	. (G2_NEW DScdSV/G1_dSV !
	582	. +(1.-G2_NEW /G1_dSV)*DFcdSV) !
	583	SphNEW = G2_NEW /G1_dSV ! Spher.Dendr.Lay.
	584	SizOLD = G2_OLD ! Size Spher.Lay.
	585	SphOLD = G1_OLD /G1_dSV ! Spher.Spher.Lay.
	586	Siz_av = ( zroNEW SizNEW+zroOLDSizOLD) ! Averaged Size
	587	Sph_av = min( zroNEW SphNEW+zroOLDSphOLD !
	588	. , unun ) ! Averaged Sphericity
	589	Den_av = min((Siz_av - ( Sph_av *DScdSV !
	590	. +(1.-Sph_av)*DFcdSV)) !
	591	. / (DDcdSV - ( Sph_av *DScdSV !
	592	. +(1.-Sph_av)*DFcdSV)) !
	593	. , unun )!
	594	DendOK = max(zero, !
	595	. sign(unun, Sph_av *DScdSV ! Small Grains
	596	. +(1.-Sph_av)*DFcdSV ! Faceted Grains
	597	. - Siz_av )) !
	598	C +... REMARQUE: le type moyen (dendritique ou non) depend
	599	C + ^^^^^^^^ de la comparaison avec le diametre optique
	600	C + d'une neige recente de dendricite nulle
	601	C +... REMARK: the mean type (dendritic or not) depends
	602	C + ^^^^^^ on the comparaison with the optical diameter
	603	C + of a recent snow having zero dendricity
	604
	605	G1diff =( -DendOK *Den_av
	606	. +(1.-DendOK)Sph_av) G1_dSV
	607	G2diff = DendOK Sph_av G1_dSV
	608	. +(1.-DendOK)*Siz_av
	609	G1 = SameOK *G1same
	610	. +(1.-SameOK)*G1diff
	611	G2 = SameOK *G2same
	612	. +(1.-SameOK)*G2diff
	613
	614	BG1sSV(ikl) = G1 !
	615	. * Bufs_N/max(epsi,Bufs_N) !
	616	BG2sSV(ikl) = G2 !
	617	. * Bufs_N/max(epsi,Bufs_N) !
	618
	619
	620	C +--Update of Buffer Layer Content & Decision about creating a new snow layer
	621	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	622	BufsSV(ikl) = Bufs_N ! [mm w.e.]
	623	NLaysv(ikl) = min(unun, !
	624	. max(zero, ! Allows to create
	625	. sign(unun,BufsSV(ikl) ! a new snow Layer
	626	. -SMndSV )) ! if Buffer > SMndSV
	627	. max(zero, ! Except if Erosion
	628	. sign(unun,0.50 ! dominates
	629	. -dsnbSV(ikl))) !
	630	. +max(zero, ! Allows to create
	631	. sign(unun,BufsSV(ikl) ! a new snow Layer
	632	. -SMndSV3.00))) ! is Buffer > SMndSV3
	633	Bdzssv(ikl) = 1.e-3BufsSV(ikl)ro_Wat ! [mm w.e.] -> [m w.e.]
	634	. /max(epsi,BrosSV(ikl))!& [m w.e.] -> [m]
	635
	636
	637
	638	END DO
	639
	640
	641
	642	! Snow Pack Discretization
	643	! ========================
	644
	645	! **********
	646	call SISVAT_zSn
	647	! **********
	648
	649	! **********
	650	! #ve call SISVAT_wEq('_zSn ',0)
	651	! **********
	652
	653
	654
	655	! Add a new Snow Layer
	656	! ====================
	657
	658	DO ikl=1,knonv
	659
	660	isnoSV(ikl) = isnoSV(ikl) +NLaysv(ikl)
	661	isn = isnoSV(ikl)
	662	dzsnSV(ikl,isn) = dzsnSV(ikl,isn) * (1-NLaysv(ikl))
	663	. + Bdzssv(ikl) * NLaysv(ikl)
	664	TsisSV(ikl,isn) = TsisSV(ikl,isn) * (1-NLaysv(ikl))
	665	. + min(TaT_SV(ikl),Tf_Sno) *NLaysv(ikl)
	666
	667	ro__SV(ikl,isn) = ro__SV(ikl,isn) * (1-NLaysv(ikl))
	668	. + Brossv(ikl) * NLaysv(ikl)
	669	eta_SV(ikl,isn) = eta_SV(ikl,isn) * (1-NLaysv(ikl)) ! + 0.
	670	agsnSV(ikl,isn) = agsnSV(ikl,isn) * (1-NLaysv(ikl)) ! + 0.
	671	G1snSV(ikl,isn) = G1snSV(ikl,isn) * (1-NLaysv(ikl))
	672	. + BG1ssv(ikl) * NLaysv(ikl)
	673	G2snSV(ikl,isn) = G2snSV(ikl,isn) * (1-NLaysv(ikl))
	674	. + BG2ssv(ikl) * NLaysv(ikl)
	675	istoSV(ikl,isn) = istoSV(ikl,isn) * (1-NLaysv(ikl))
	676	. + max(zer0,sign(un_1,TaT_SV(ikl)
	677	. -Tf_Sno-eps_21)) * istdSV(2)
	678	. * NLaysv(ikl)
	679	BufsSV(ikl) = BufsSV(ikl) * (1-NLaysv(ikl))
	680	NLaysv(ikl) = 0
	681
	682
	683	END DO
	684
	685
	686	! Snow Pack Thickness
	687	! -------------------
	688
	689	DO ikl=1,knonv
	690	z_snsv(ikl) = 0.0
	691	END DO
	692	DO isn=1,nsno
	693	DO ikl=1,knonv
	694	z_snsv(ikl) = z_snsv(ikl) + dzsnSV(ikl,isn)
	695	zzsnsv(ikl,isn) = z_snsv(ikl)
	696	END DO
	697	END DO
	698
	699
	700
	701	END IF
	702
	703
	704
	705	! Soil Albedo: Soil Humidity Correction
	706	! ==========================================
	707
	708	! REFERENCE: McCumber and Pielke (1981), Pielke (1984)
	709	! ^^^^^^^^^
	710	DO ikl=1,knonv
	711	albssv(ikl) =
	712	. alb0SV(ikl) *(1.0-min(half,eta_SV( ikl,0)
	713	. /etadSV(isotSV(ikl))))
	714	! REMARK: Albedo of Water Surfaces (isotSV=0):
	715	! ^^^^^^ alb0SV := 2 X effective value, while
	716	! eta_SV := etadSV
	717	END DO
	718
	719
	720	! Snow Pack Optical Properties
	721	! ============================
	722
	723	IF (SnoMod) THEN
	724
	725	! ******
	726	call SnOptP(jjtime)
	727	! ******
	728
	729	ELSE
	730	DO ikl=1,knonv
	731	sEX_sv(ikl,1) = 1.0
	732	sEX_sv(ikl,0) = 0.0
	733	albisv(ikl) = albssv(ikl)
	734	END DO
	735	END IF
	736
	737
	738
	739	! Soil optical properties
	740	! =============================
	741	!Etienne: as in inlandis we do not call vgopt, we need to define
	742	!the albedo alb_SV and to calculate the
	743	!absorbed Solar Radiation by Surfac (Normaliz)[-] SoSosv
	744
	745
	746	DO ikl=1,klonv
	747
	748	e_pRad = 2.5 * coszSV(ikl) ! exponential argument,
	749	! V/nIR radiation partitioning,
	750	! DR97, 2, eqn (2.53) & (2.54)
	751	e1pRad = 1.-exp(-e_pRad) ! exponential, V/nIR Rad. Part.
	752	exdRad= 1.
	753
	754	! Visible Part of the Solar Radiation Spectrum (V, 0.4--0.7mi.m)
	755	A_Rad0 = 0.25 + 0.697 * e1pRad ! Absorbed Vis. Radiation
	756	absg_V = (1.-albisv(ikl))(A_Rad0exdRad) !
	757
	758	! Near-IR Part of the Solar Radiation Spectrum (nIR, 0.7--2.8mi.m)
	759
	760	A_Rad0 = 0.80 + 0.185 * e1pRad ! Absorbed nIR. Radiation
	761	absgnI = (1.-albisv(ikl))(A_Rad0exdRad) !
	762
	763	SoSosv(ikl) = (absg_V+absgnI)*0.5d0
	764
	765	alb_SV(ikl) = albisv(ikl)
	766
	767	END DO
	768
	769	! **********
	770	! #ve call SISVAT_wEq('SnOptP',0)
	771	! **********
	772
	773
	774	! Surface Emissivity (Etienne: simplified calculation for landice)
	775	! =============================================================
	776	!
	777	DO ikl=1,knonv
	778	LSnMsk = min( 1,isnoSV(ikl))
	779	Eso_sv(ikl)= EmiSol(1-LSnMsk)+EmiSnoLSnMsk ! Sol+Sno Emissivity
	780	emi_SV(ikl)= EmiSol(1-LSnMsk) + EmiSnoLSnMsk
	781	END DO
	782
	783
	784
	785
	786	! Upward IR (INPUT, from previous time step)
	787	! ===================================================================
	788
	789	DO ikl=1,knonv
	790	! #e1 Enrsvd(ikl) = - IRs_SV(ikl)
	791	IRupsv(ikl) = IRs_SV(ikl)
	792	END DO
	793
	794
	795	! Turbulence
	796	! ==========
	797
	798	! Latent Heat of Vaporization/Sublimation
	799	! ---------------------------------------
	800
	801	DO ikl=1,knonv
	802	SnoWat = min(isnoSV(ikl),0)
	803	Lx_H2O(ikl) =
	804	. (1.-SnoWat) * LhvH2O
	805	. + SnoWat (LhsH2O (1.-eta_SV(ikl,isnoSV(ikl)))
	806	. +LhvH2O * eta_SV(ikl,isnoSV(ikl)) )
	807	END DO
	808
	809
	810
	811
	812	! Aerodynamic Resistance
	813	! ^^^^^^^^^^^^^^^^^^^^^^
	814
	815
	816	DO ikl=1,knonv
	817	ram_sv(ikl) = 1./(cdM_SV(ikl)*max(VV__SV(ikl),eps6))
	818	rah_sv(ikl) = 1./(cdH_SV(ikl)*max(VV__SV(ikl),eps6))
	819	END DO
	820
	821
	822
	823	! Soil Energy Balance
	824	! =====================
	825
	826
	827	if (iflag_tsurf_inlandsis .eq. 0) then
	828
	829	call SISVAT_TSo
	830
	831	else
	832
	833	call SISVAT_TS2
	834
	835	end if
	836
	837
	838	! **********
	839	! #ve call SISVAT_wEq('_TSo ',0)
	840	! **********
	841
	842
	843
	844	! Soil Water Potential
	845	! ------------------------
	846
	847	DO isl=-nsol,0
	848	DO ikl=1,knonv
	849	ist = isotSV(ikl) ! Soil Type
	850	psi_sv(ikl,isl) = psidSV(ist) ! DR97, Eqn.(3.34)
	851	. *(etadSV(ist) /max(eps6,eta_SV(ikl,isl))) !
	852	. **bCHdSV(ist) !
	853
	854
	855	! Soil Hydraulic Conductivity
	856	! ---------------------------
	857
	858	Khydsv(ikl,isl) = s2__SV(ist) ! DR97, Eqn.(3.35)
	859	. (eta_SV(ikl,isl)(2.bCHdSV(ist)+3.)) !
	860	END DO
	861	END DO
	862
	863
	864	! Melting / Refreezing in the Snow Pack
	865	! =====================================
	866
	867	IF (SnoMod) THEN
	868
	869	! **********
	870	call SISVAT_qSn
	871	! **********
	872
	873	! **********
	874	! #ve call SISVAT_wEq('_qSn ',0)
	875	! **********
	876
	877
	878	! Snow Pack Thickness
	879	! -------------------
	880
	881	DO ikl=1,knonv
	882	z_snsv(ikl) = 0.0
	883	END DO
	884	DO isn=1,nsno
	885	DO ikl=1,knonv
	886	z_snsv(ikl) = z_snsv(ikl) + dzsnSV(ikl,isn)
	887	zzsnsv(ikl,isn) = z_snsv(ikl)
	888	END DO
	889	END DO
	890
	891
	892	! Energy in Excess is added to the first Soil Layer
	893	! -------------------------------------------------
	894	DO ikl=1,knonv
	895	z_snsv(ikl) = max(zer0,
	896	. sign(un_1,eps6-z_snsv(ikl)))
	897	TsisSV(ikl,0) = TsisSV(ikl,0) + EExcsv(ikl)
	898	. /(rocsSV(isotSV(ikl))
	899	. +rcwdSV*eta_SV(ikl,0))
	900	EExcsv(ikl) = 0.
	901	END DO
	902
	903
	904	END IF
	905
	906
	907	! Soil Water Balance
	908	! =====================
	909
	910	! **********
	911	call SISVAT_qSo
	912	! #m0. (Wats_0,Wats_1,Wats_d)
	913	! **********
	914
	915
	916	! Surface Fluxes
	917	! =====================
	918
	919	DO ikl=1,knonv
	920	IRdwsv(ikl)=IRd_SV(ikl)*Eso_sv(ikl) ! Downward IR
	921	! IRdwsv(ikl)=tau_sv(ikl) IRd_SV(ikl)Eso_sv(ikl) ! Downward IR
	922	! . +(1.0-tau_sv(ikl))IRd_SV(ikl)Evg_sv(ikl) ! ! Etienne, remove vegetation component
	923	IRupsv(ikl) = IRupsv(ikl) ! Upward IR
	924	IRu_SV(ikl) = -IRupsv(ikl) ! Upward IR
	925	. +IRd_SV(ikl) ! (effective)
	926	. -IRdwsv(ikl) ! (positive)
	927
	928	TBr_sv(ikl) =sqrt(sqrt(IRu_SV(ikl)/StefBo)) ! Brightness
	929	! ! Temperature
	930	uts_SV(ikl) = (HSv_sv(ikl) +HSs_sv(ikl)) ! uT
	931	. /(rhT_SV(ikl) *cp) !
	932	uqs_SV(ikl) = (HLv_sv(ikl) +HLs_sv(ikl)) ! uq
	933	. /(rhT_SV(ikl) *LhvH2O) !
	934	LMO_SV(ikl) = TaT_SV(ikl)(us__SV(ikl)*3)
	935	. /gravit/uts_SV(ikl)/vonKrm ! MO length
	936
	937
	938	! Surface Temperature
	939	! ^^^^^^^^^^^^^^^^^^^^
	940	! Tsrfsv(ikl) =TsisSV(ikl,isnoSV(ikl))
	941
	942	! Etienne: extrapolation from the two uppermost levels:
	943
	944	if (isnoSV(ikl) >=2) then
	945	zm1=-dzsnSV(ikl,isnoSV(ikl))/2.
	946	zm2=-(dzsnSV(ikl,isnoSV(ikl)) + dzsnSV(ikl,isnoSV(ikl)-1)/2.)
	947	else if (isnoSV(ikl) .EQ. 1) then
	948	zm1=-dzsnSV(ikl,isnoSV(ikl))/2.
	949	zm2=-(dzsnSV(ikl,isnoSV(ikl))+dz_dSV(0)/2.)
	950	else
	951	zm1=-dz_dSV(0)/2.
	952	zm2=-(dz_dSV(0)+dz_dSV(-1)/2.)
	953
	954	end if
	955
	956	coefslope=(TsisSV(ikl,isnoSV(ikl))-TsisSV(ikl,isnoSV(ikl)-1))
	957	. /(zm1-zm2)
	958	Tsrfsv(ikl)=TsisSV(ikl,isnoSV(ikl))+coefslope*(0. - zm1)
	959
	960
	961	END DO
	962
	963
	964	! Snow Pack Properties (sphericity, dendricity, size)
	965	! ===================================================
	966
	967	IF (SnoMod) THEN
	968
	969	! **********
	970	call SISVAT_GSn
	971	! **********
	972
	973	! **********
	974	! #ve call SISVAT_wEq('_GSn ',0)
	975	! **********
	976
	977
	978
	979	END IF
	980
	981
	982	! Roughness Length for next time step
	983	!====================================
	984
	985	! Note that in INLANDSIS, we treat only ice covered surfaces so calculation
	986	! of z0 is much simpler (no subgrid fraction of ocean or land)
	987	! old calculations are commented below
	988
	989
	990	C +--Roughness Length for Momentum
	991	C + -----------------------------
	992
	993	C +--Land+Sea-Ice / Ice-free Sea Mask
	994	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	995	DO ikl=1,klonv
	996	IcIndx(ikl) = 0
	997	ENDDO
	998	DO isn=1,nsno
	999	DO ikl=1,klonv
	1000	IcIndx(ikl) = max(IcIndx(ikl),
	1001	. isn*max(0,
	1002	. sign(1,
	1003	. int(ro__SV(ikl,isn)-900.))))
	1004	ENDDO
	1005	ENDDO
	1006
	1007	DO ikl=1,klonv
	1008	LISmsk = 1. ! in inlandsis, land only
	1009	IceMsk = max(0,sign(1 ,IcIndx(ikl)-1) )
	1010	SnoMsk = max(min(isnoSV(ikl)-iiceSV(ikl),1),0)
	1011
	1012
	1013
	1014	Z0mLnd =max( Z0_ICE , 5.e-5 ) ! Min set := Z0 on *
	1015
	1016	C +--Z0 Smooth Regime over Snow (Andreas 1995, CRREL Report 95-16, p. 8)
	1017	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^
	1018	Z0m_nu = 5.e-5 ! z0s~(10-d)*exp(-vonkar/sqrt(1.1e-03))
	1019
	1020	C +--Z0 Saltat.Regime over Snow (Gallee et al., 2001, BLM 99 (19) p.11)
	1021	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^
	1022	u2star = us__SV(ikl) *us__SV(ikl)
	1023	Z0mBSn = u2star *0.536e-3 - 61.8e-6
	1024	Z0mBSn = max(Z0mBS0 ,Z0mBSn)
	1025
	1026	C +--Z0 Smooth + Saltat. Regime
	1027	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^
	1028	Z0enSV(ikl) = Z0m_nu
	1029	. + Z0mBSn
	1030
	1031	C +--Rough Snow Surface Roughness Length (Typical Value)
	1032	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	1033	c #tz Z0m_Sn = 0.250e-3 ! Andreas 1995, CRREL Report 95-16, fig.1&p.2
	1034	! z0r~(10-d)*exp(-vonkar/sqrt(1.5e-03))-5.e-5
	1035	Z0m_Sn = 2.000e-3 ! Calibration of MAR
	1036	c #TZ Z0m_Sn = 1.000e-3 ! Exemple Tuning in RACMO
	1037	c #TZ Z0m_Sn = 0.500e-3 ! Exemple Tuning in MAR
	1038
	1039	C +--Rough Snow Surface Roughness Length (Variable Sastrugi Height)
	1040	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	1041	A_Fact = 1.0000 ! Andreas et al., 2004, p.4
	1042	! ams.confex.com/ams/pdfpapers/68601.pdf
	1043
	1044	! Parameterization of z0 dependance on Temperature (C. Amory, 2017)
	1045	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	1046	! Z0=f(T) deduced from observations, Adelie Land, dec2012-dec2013
	1047	coefa = 0.1658 !0.1862 !Ant
	1048	coefb = -50.3869 !-55.7718 !Ant
	1049	ta1 = 253.15 !255. Ant
	1050	ta2 = 273.15
	1051	ta3 = 273.15+3
	1052	z01 = exp(coefa*ta1 + coefb) !~0.2 ! ~0.25 mm
	1053	z02 = exp(coefa*ta2 + coefb) !~6 !~7 mm
	1054	z03 = z01
	1055	coefc = log(z03/z02)/(ta3-ta2)
	1056	coefd = log(z03)-coefc*ta3
	1057	if (TaT_SV(ikl) .lt. ta1) then
	1058	Z0_obs = z01
	1059	else if (TaT_SV(ikl).ge.ta1 .and. TaT_SV(ikl).lt.ta2) then
	1060	Z0_obs = exp(coefa*TaT_SV(ikl) + coefb)
	1061	else if (TaT_SV(ikl).ge.ta2 .and. TaT_SV(ikl).lt.ta3) then
	1062	! if st > 0, melting induce smooth surface
	1063	Z0_obs = exp(coefc*TaT_SV(ikl) + coefd)
	1064	else
	1065	Z0_obs = z03
	1066	endif
	1067
	1068
	1069	! pour le moment, on choisit une valeur fixe
	1070	Z0_obs = 1.000e-3
	1071
	1072	cCA Snow roughness lenght deduced from observations
	1073	cCA (parametrization if no Blowing Snow)
	1074	cCA ----------------------------------- C. Agosta 09-2016 -----
	1075	cCA Substract Z0enSV(ikl) because re-added later in Z0mnSV(ikl)
	1076	Z0m_Sn = Z0_obs - Z0enSV(ikl)
	1077	cCA -----------------------------------------------------------
	1078
	1079	param = Z0_obs/1. ! param(s) \| 1.(m/s)=TUNING
	1080
	1081	c #SZ Z0Sa_N = (us__SV(ikl) -0.2)*param ! 0.0001=TUNING
	1082	c #SZ. * max(zero,sign(unun,TfSnow-eps9
	1083	c #SZ. -TsisSV(ikl , isnoSV(ikl))))
	1084	!!#SZ Z0SaSi = max(zero,sign(unun,Z0Sa_N ))! 1 if erosion
	1085	c #SZ Z0SaSi = max(zero,sign(unun,zero -eps9 -uss_SV(ikl)))!
	1086	c #SZ Z0Sa_N = max(zero, Z0Sa_N)
	1087	c #SZ Z0SaSV(ikl) =
	1088	c #SZ. max(Z0SaSV(ikl) ,Z0SaSV(ikl)
	1089	c #SZ. + Z0SaSi(Z0Sa_N-Z0SaSV(ikl))exp(-dt__SV/43200.))
	1090	c #SZ. - min(dz0_SV(ikl) , Z0SaSV(ikl))
	1091
	1092	c #SZ A_Fact = Z0SaSV(ikl) * 5.0/0.15 ! A=5 if h~10cm
	1093	C +... CAUTION: The influence of the sastrugi direction is not yet included
	1094
	1095	c #SZ Z0m_Sn = Z0SaSV(ikl) !
	1096	c #SZ. - Z0m_nu !
	1097
	1098	C +--Z0 Saltat.Regime over Snow (Shao & Lin, 1999, BLM 91 (46) p.222)
	1099	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^
	1100	c #ZN sqrrZ0 = usthSV(ikl)/max( us__SV(ikl),0.001)
	1101	c #ZN sqrrZ0 = min( sqrrZ0 ,0.999)
	1102	c #ZN Z0mBSn = 0.55 0.55 exp(-sqrrZ0 *sqrrZ0)
	1103	c #ZN. us__SV(ikl) us__SV(ikl)grvinv0.5
	1104
	1105	C +--Z0 Smooth + Saltat. Regime (Shao & Lin, 1999, BLM 91 (46) p.222)
	1106	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^
	1107	c #ZN Z0enSV(ikl) = (Z0m_nu ** sqrrZ0 )
	1108	c #ZN. * (Z0mBSn **(1.-sqrrZ0))
	1109	c #ZN Z0enSV(ikl) = max(Z0enSV(ikl), Z0m_nu)
	1110
	1111	C +--Z0 Smooth Regime over Snow (Andreas etAl., 2004
	1112	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^ ams.confex.com/ams/pdfpapers/68601.pdf)
	1113	c #ZA Z0m_nu = 0.135*akmol / max(us__SV(ikl) , epsi)
	1114
	1115	C +--Z0 Saltat.Regime over Snow (Andreas etAl., 2004
	1116	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^ ams.confex.com/ams/pdfpapers/68601.pdf)
	1117	c #ZA Z0mBSn = 0.035u2star grvinv
	1118
	1119	C +--Z0 Smooth + Saltat. Regime (Andreas etAl., 2004
	1120	! ( used by Erosion) ams.confex.com/ams/pdfpapers/68601.pdf)
	1121	! ^^^^^^^^^^^^^^^^^^^^^^^^^^
	1122	c #ZA Z0enSV(ikl) = Z0m_nu
	1123	c #ZA. + Z0mBSn
	1124
	1125	C +--Z0 Rough Regime over Snow (Andreas etAl., 2004
	1126	C + (.NOT. used by Erosion) ams.confex.com/ams/pdfpapers/68601.pdf)
	1127	! ^^^^^^^^^^^^^^^^^^^^^^^^^^
	1128	!!#ZA u2star = (us__SV(ikl) -0.1800) / 0.1
	1129	!!#ZA Z0m_Sn =A_FactZ0mBSn exp(-u2star*u2star)
	1130	c #ZA Z0m_90 =(10.-0.025*VVs_SV(ikl)/5.)
	1131	c #ZA. exp(-0.4/sqrt(.00275+.00001max(0.,VVs_SV(ikl)-5.)))
	1132	c #ZA Z0m_Sn = DDs_SV(ikl)* Z0m_90 / 45.
	1133	c #ZA. - DDs_SV(ikl)DDs_SV(ikl) Z0m_90 /(90.*90.)
	1134
	1135	C +--Z0 (Erosion) over Snow (instantaneous or time average)
	1136	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^
	1137	Z0e_SV(ikl) = Z0enSV(ikl)
	1138	Z0e_SV(ikl) = Z0emSV(ikl)
	1139
	1140	C +--Momentum Roughness Length
	1141	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^ ! Contribution of
	1142	Z0mnSV(ikl) = Z0mLnd ! land Form
	1143	. + (Z0m_Sn ! Sastrugi Form
	1144	. + Z0enSV(ikl)) *SnoMsk ! Snow Erosion
	1145
	1146
	1147	C +--GIS Roughness Length
	1148	C + ^^^^^^^^^^^^^^^^^^^^^
	1149	c #GL Z0mnSV(ikl) =
	1150	c #GL. (1-LSmask(ikl)) * Z0mnSV(ikl)
	1151	c #GL. + LSmask(ikl) * max(Z0mnSV(ikl),max(Z0_GIM,
	1152	c #GL. Z0_GIM+
	1153	c #GL. (0.0032-Z0_GIM)*(ro__SV(ikl,isnoSV(ikl))-600.) !
	1154	c #GL. /(920.00 -600.))) !
	1155
	1156	C +--Mom. Roughness Length, Instantaneous OR Box Moving Average in Time
	1157	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	1158	Z0m_SV(ikl) = Z0mnSV(ikl) ! Z0mnSV instant.
	1159	! Z0m_SV(ikl) = Z0mmSV(ikl) ! Z0mnSV Average
	1160
	1161	C +--Corrected Threshold Friction Velocity before Erosion ! Marticorena and
	1162	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ! Bergametti 1995
	1163	! not used anymore since Marticorena and Bergametti disabled !CK 18/07/2018
	1164	cc #BS Z0e_SV(ikl) = min(Z0m_SV(ikl),Z0e_SV(ikl)) !
	1165	cc #MB f_eff= log(0.35(0.1 /Z0e_SV(ikl))*0.8) ! JGR 100
	1166	cc #MB f_eff=1.-(log( Z0m_SV(ikl)/Z0e_SV(ikl) ))! (20) p. 16420
	1167	cc #MB. /(max( f_eff ,epsi ))! p.16426 2nd ?
	1168	cc #MB f_eff= max( f_eff ,epsi )! CONTROL
	1169	cc #MB f_eff=1.0 -(1.0 - f_eff) /5.00 ! TUNING
	1170	cc #MB f_eff= min( f_eff ,1.00 )!
	1171	cc #MB usthSV(ikl) = usthSV(ikl)/f_eff !
	1172
	1173
	1174
	1175	C +--Roughness Length for Scalars
	1176	C + ----------------------------
	1177
	1178	Z0hnSV(ikl) = Z0mnSV(ikl)/ 7.4
	1179	c #SH Z0hnSV(ikl) = Z0mnSV(ikl)/100.0
	1180	C + Z0h = Z0m /100.0 over the Sahel
	1181	C + (Taylor & Clark, QJRMS 127,p864)
	1182
	1183	c #RN rstar = Z0mnSV(ikl) * us__SV(ikl) / akmol
	1184	c #RN rstar = max(epsi,min(rstar,thous))
	1185	c #RN alors = log(rstar)
	1186	c #RN rstar0 = 1.250e0 * max(zero,sign(unun,0.135e0 - rstar))
	1187	c #RN. +(1. - max(zero,sign(unun,0.135e0 - rstar)))
	1188	c #RN. (0.149e0 max(zero,sign(unun,2.500e0 - rstar))
	1189	c #RN. + 0.317e0
	1190	c #RN. *(1. - max(zero,sign(unun,2.500e0 - rstar))))
	1191	c #RN rstar1 = 0. * max(zero,sign(unun,0.135e0 - rstar))
	1192	c #RN. +(1. - max(zero,sign(unun,0.135e0 - rstar)))
	1193	c #RN. (-0.55e0 max(zero,sign(unun,2.500e0 - rstar))
	1194	c #RN. - 0.565
	1195	c #RN. *(1. - max(zero,sign(unun,2.500e0 - rstar))))
	1196	c #RN rstar2 = 0. * max(zero,sign(unun,0.135e0 - rstar))
	1197	c #RN. +(1. - max(zero,sign(unun,0.135e0 - rstar)))
	1198	c #RN. (0. max(zero,sign(unun,2.500e0 - rstar))
	1199	c #RN. - 0.183
	1200	c #RN. *(unun - max(zero,sign(unun,2.500e0 - rstar))))
	1201
	1202	cXF #RN does not work over bare ice
	1203	cXF MAR is then too warm and not enough melt
	1204
	1205	c #RN if(ro__SV(ikl,isnoSV(ikl))>50
	1206	c #RN. .and.ro__SV(ikl,isnoSV(ikl))<roSdSV)then
	1207
	1208	c #RN Z0hnSV(ikl) = max(zero
	1209	c #RN. , sign(unun,zzsnsv(ikl,isnoSV(ikl))-epsi))
	1210	c #RN. * exp(rstar0+rstar1alors+rstar2alors*alors)
	1211	c #RN. * 0.001e0 + Z0hnSV(ikl) * ( 1. - max(zero
	1212	c #RN. , sign(unun,zzsnsv(ikl,isnoSV(ikl))-epsi)))
	1213
	1214	c #RN endif
	1215
	1216	Z0h_SV(ikl) = Z0hnSV(ikl)
	1217	! Z0h_SV(ikl) = Z0hmSV(ikl)
	1218
	1219
	1220	c #MT Z0m_SV(ikl) = max(2.0e-6 ,Z0m_SV(ikl)) ! Min Z0_m (Garrat Scheme)
	1221	! Z0m_SV(ikl) = min(Z0m_SV(ikl),za__SV(ikl)*0.3333)
	1222
	1223
	1224	END DO
	1225
	1226
	1227	return
	1228	end
	1229
	1230
	1231
	1232
	1233
	1234
	1235
	1236
	1237
	1238
	1239

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

Download in other formats:

Original Format