Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

sisvat_qsn.f90 @ 5413

Last change on this file since 5413 was 5246, checked in by abarral, 2 months ago
Convert fixed-form to free-form sources .F -> .{f,F}90 (WIP: some .F remain, will be handled in subsequent commits)
File size: 26.4 KB

Rev	Line
[3792]	1
	2
[5246]	3	subroutine SISVAT_qSn &
	4	( &
	5	! #e1. EqSn_0,EqSn_1,EqSn_d
	6	! #m1. ,SIsubl,SImelt,SIrnof
	7	)
[3900]	8
[5246]	9	! +------------------------------------------------------------------------+
	10	! \| MAR SISVAT_qSn Fri 29-Jul-2011 MAR \|
	11	! \| SubRoutine SISVAT_qSn updates the Snow Water Content \|
	12	! +------------------------------------------------------------------------+
	13	! \| \|
	14	! \| PARAMETERS: knonv: Total Number of columns = \|
	15	! \| ^^^^^^^^^^ = Total Number of continental grid boxes \|
	16	! \| X Number of Mosaic Cell per grid box \|
	17	! \| \|
	18	! \| INPUT: isnoSV = total Nb of Ice/Snow Layers \|
	19	! \| ^^^^^ \|
	20	! \| \|
	21	! \| INPUT: TaT_SV : SBL Top Temperature [K] \|
	22	! \| ^^^^^ dt__SV : Time Step [s] \|
	23	! \| \|
	24	! \| INPUT / drr_SV : Rain Intensity [kg/m2/s] \|
	25	! \| OUTPUT: dzsnSV : Snow Layer Thickness [m] \|
	26	! \| ^^^^^^ eta_SV : Snow Water Content [m3/m3] \|
	27	! \| ro__SV : Snow/Soil Volumic Mass [kg/m3] \|
	28	! \| TsisSV : Soil/Ice Temperatures (layers -nsol,-nsol+1,..,0)\|
	29	! \| & Snow Temperatures (layers 1,2,...,nsno) [K] \|
	30	! \| \|
	31	! \| OUTPUT: SWS_SV : Surficial Water Status \|
	32	! \| ^^^^^^ \|
	33	! \| EExcsv : Snow Energy in Excess, initial Forcing [J/m2] \|
	34	! \| EqSn_d : Snow Energy in Excess, remaining [J/m2] \|
	35	! \| EqSn_0 : Snow Energy, before Phase Change [J/m2] \|
	36	! \| EqSn_1 : Snow Energy, after Phase Change [J/m2] \|
	37	! \| SIsubl : Snow sublimed/deposed Mass [mm w.e.] \|
	38	! \| SImelt : Snow Melted Mass [mm w.e.] \|
	39	! \| SIrnof : Surficial Water + Run OFF Change [mm w.e.] \|
	40	! \| \|
	41	! \| Internal Variables: \|
	42	! \| ^^^^^^^^^^^^^^^^^^ \|
	43	! \| \|
	44	! \| # OPTIONS: #E0: IO for Verification: Energy Budget \|
	45	! \| # ^^^^^^^ \|
	46	! \| # #su: IO for Verification: Slush Diagnostic \|
	47	! \| \|
	48	! \| \|
	49	! +------------------------------------------------------------------------+
[3792]	50
	51
[3900]	52
[3792]	53
[5246]	54	! +--Global Variables
	55	! + ================
[3792]	56
[5246]	57	use VARphy
	58	use VAR_SV
	59	use VARdSV
	60	use VAR0SV
	61	use VARxSV
	62	use VARySV
	63	use surface_data, only: is_ok_slush,opt_runoff_ac
[3792]	64
	65
[5246]	66	IMPLICIT NONE
[3792]	67
	68
[5246]	69	! Energy Budget
	70	! ~~~~~~~~~~~~~~~~~~~~~~
	71	! #e1 real EqSn_d(knonv) ! Energy in Excess, initial
	72	! #e1 real EqSn_0(knonv) ! Snow Energy, befor Phase Change
	73	! #vm real EqSn01(knonv) ! Snow Energy, after Phase Change
	74	! #vm real EqSn02(knonv) ! Snow Energy, after Phase Change
	75	! ! .AND. Last Melting
	76	! #e1 real EqSn_1(knonv) ! Snow Energy, after Phase Change
	77	! ! .AND. Mass Redistr.
	78	! Snow/Ice (Mass) Budget
	79	! ~~~~~~~~~~~~~~~~~~~~~~
	80	! #m1 real SIsubl(knonv) ! Snow Deposed Mass
	81	! #m1 real SImelt(knonv) ! Snow Melted Mass
	82	! #m1 real SIrnof(knonv) ! Local Surficial Water + Run OFF
[3900]	83
[3792]	84
[5246]	85	! +--Internal Variables
	86	! + ==================
[3900]	87
[5246]	88	integer :: ikl ,isn !
	89	integer :: nh ! Non erodible Snow: up.lay.Index
	90	integer :: LayrOK ! 1 (0) if In(Above) Snow Pack
	91	integer :: k_face ! 1 (0) if Crystal(no) faceted
	92	integer :: LastOK ! 1 ==> 1! Snow Layer
	93	integer :: NOLayr ! 1 Layer Update
	94	integer :: noSnow(knonv) ! Nb of Layers Updater
	95	integer :: kSlush ! Slush Switch
	96	real :: dTSnow ! Temperature [C]
	97	real :: EExdum(knonv) ! Energy in Excess when no Snow
	98	real :: OKmelt ! 1 (0) if (no) Melting
	99	real :: EnMelt ! Energy in excess, for Melting
	100	real :: SnHLat ! Energy consumed in Melting
	101	real :: AdEnrg,B_Enrg ! Additional Energy from Vapor
	102	real :: dzVap0,dzVap1 ! Vaporized Thickness [m]
	103	real :: dzMelt(knonv) ! Melted Thickness [m]
	104	real :: rosDry ! Snow volumic Mass if no Water in
	105	real :: PorVol ! Pore volume
	106	real :: PClose ! Pore Hole Close OFF Switch
	107	real :: SGDiam ! Snow Grain Diameter
	108	real :: SGDmax ! Max. Snow Grain Diameter
	109	real :: rWater ! Retained Water [kg/m2]
	110	real :: drrNEW ! New available Water [kg/m2]
	111	real :: rdzNEW ! Snow Mass [kg/m2]
	112	real :: rdzsno ! Snow Mass [kg/m2]
	113	real :: EnFrez ! Energy Release in Freezing
	114	real :: WaFrez ! Water consumed in Melting
	115	real :: RapdOK ! 1. ==> Snow melts rapidly
	116	real :: ThinOK ! 1. ==> Snow Layer is thin
	117	real :: dzepsi ! Minim. Snow Layer Thickness (!)
	118	real :: dz_Min ! Minim. Snow Layer Thickness
	119	real :: z_Melt ! Last (thin) Layer Melting
	120	real :: rusnew ! Surficial Water Thickness [mm]
	121	real :: zWater ! Max Slush Water Thickness [mm]
	122	real :: zSlush ! Slush Water Thickness [mm]
	123	real :: ro_new ! New Snow/ice Density [kg/m3]
	124	real :: zc,zt ! Non erod.Snow Thickness[mm w.e.]
	125	real :: rusnSV0(knonv)
	126	real :: Tsave
[3792]	127
[5246]	128	! +--OUTPUT of SISVAT Trace Statistics (see assignation in PHY_SISVAT)
	129	! + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	130	integer :: isnnew,isinew,isnUpD,isnitr
	131
	132	! OUTPUT in SISVAT at specified i,j,k,n (see assignation in PHY_SISVAT)
	133	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	134	! #wx integer iSV_v1,jSV_v1,nSV_v1,kSV_v1,lSV_v1
	135	! #wx common/SISVAT_EV/ iSV_v1,jSV_v1,nSV_v1,kSV_v1,lSV_v1
	136
	137	! +--Energy and Mass Budget
	138	! + ~~~~~~~~~~~~~~~~~~~~~~
	139	! #vm real WqSn_0(knonv) ! Snow Water+Forcing Initial
	140	! #vm real WqSn_1(knonv) ! Snow Water+Forcing, Final
	141	! #vm logical emopen ! IO Switch
	142	! #vm common/Se_qSn_L/emopen !
	143	! #vm integer no_err !
	144	! #vm common/Se_qSn_I/no_err !
	145	! #vm real hourer,timeer !
	146	! #vm common/Se_qSn_R/timeer !
	147
	148	! +--Slush Diagnostic: IO
	149	! + ~~~~~~~~~~~~~~~~~~~~
	150	! #vu logical su_opn ! IO Switch
	151	! #vu common/SI_qSn_L/su_opn !
	152
	153
	154	! +--DATA
	155	! + ====
	156
	157	data dzepsi/0.0001/ ! Minim. Snow Layer Thickness (!)
	158	! #?? data dz_Min/0.005/ ! Minim. Snow Layer Thickness
	159	! ... Warning: Too high for Col de Porte: precludes 1st snow (layer) apparition
	160	data dz_Min/2.5e-3/ ! Minim. Snow Layer Thickness
	161	data SGDmax/0.003/ ! Maxim. Snow Grain Diameter [m]
	162	! ! (Rowe et al. 1995, JGR p.16268)
	163
	164	! +--Energy Budget (IN)
	165	! + ==================
	166
	167	! #e1 DO ikl=1,knonv
	168	! #e1 EqSn_0(ikl) = 0.
	169	! #e1 END DO
	170	! #e1 DO isn=nsno,1,-1
	171	! #e1 DO ikl=1,knonv
	172	! #e1 EqSn_0(ikl) = EqSn_0(ikl) + ro__SV(ikl,isn) *dzsnSV(ikl,isn)
	173	! #e1. (Cn_dSV (TsisSV(ikl,isn) -TfSnow )
	174	! #e1. -Lf_H2O *(1. -eta_SV(ikl,isn)))
	175	! #e1 END DO
	176	! #e1 END DO
	177
	178
	179	! +--Water Budget (IN)
	180	! + ==================
	181
	182	! #vm DO ikl=1,knonv
	183	! #vm WqSn_0(ikl) = drr_SV(ikl) * dt__SV
	184	! #vm. +rusnSV(ikl)
	185	! #vm END DO
	186	! #vm DO isn=nsno,1,-1
	187	! #vm DO ikl=1,knonv
	188	! #vm WqSn_0(ikl) = WqSn_0(ikl) + ro__SV(ikl,isn) *dzsnSV(ikl,isn)
	189	! #vm END DO
	190	! #vm END DO
	191
	192
	193	! +--Snow Melt Budget
	194	! + ================
	195
	196	! #m1 DO ikl=1,knonv
	197	! #m1 SImelt(ikl) = 0.
	198	! #m1 SIrnof(ikl) = rusnSV(ikl) + RnofSV(ikl) * dt__SV
	199	! #m1 END DO
	200
	201
	202	! +--Initialization
	203	! + ==============
	204
	205	DO ikl=1,knonv
	206	noSnow(ikl) = 0 ! Nb of Layers Updater
	207	ispiSV(ikl) = 0 ! Pore Hole Close OFF Index
	208	! ! (assumed to be the Top of
	209	! ! the surimposed Ice Layer)
	210	zn5_SV(ikl) = 0.
	211	rusnSV0(ikl) = 0.
	212
	213	END DO
	214
	215
	216	! +--Melting/Freezing Energy
	217	! + =======================
	218
	219	! +...REMARK: Snow liquid Water Temperature assumed = TfSnow
	220	! + ^^^^^^
	221	DO ikl=1,knonv
	222	EExdum(ikl) = drr_SV(ikl) * C__Wat *(TaT_SV(ikl)-TfSnow) &
	223	* dt__SV
	224	EExcsv(ikl) = EExdum(ikl) * min(1,isnoSV(ikl)) ! Snow exists
	225	EExdum(ikl) = EExdum(ikl) - EExcsv(ikl) !
	226	! #e1 EqSn_d(ikl) = EExcsv(ikl) !
	227	END DO
	228
	229
	230	! +--Surficial Water Status
	231	! + ----------------------
	232
	233	DO ikl=1,knonv
	234	SWS_SV(ikl) = max(zero,sign(unun,TfSnow &
	235	-TsisSV(ikl,isnoSV(ikl))))
	236	END DO
	237
	238	DO ikl=1,knonv
	239
	240	DO isn=min(nsno,isnoSV(ikl)+1),1,-1
	241	! EV DO isn=nsno,1,-1
	242	! +--Energy, store Previous Content
	243	! + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	244	dTSnow = TsisSV(ikl,isn) - TfSnow
	245	EExcsv(ikl) = EExcsv(ikl) &
	246	+ ro__SV(ikl,isn) * Cn_dSV * dTSnow &
	247	* dzsnSV(ikl,isn)
	248	TsisSV(ikl,isn) = TfSnow
	249
	250	! +--Water, store Previous Content
	251	! + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	252	drr_SV(ikl) = drr_SV(ikl) &
	253	+ ro__SV(ikl,isn) * eta_SV(ikl,isn) &
	254	* dzsnSV(ikl,isn) &
	255	/ dt__SV
	256	ro__SV(ikl,isn) = &
	257	ro__SV(ikl,isn) *(1. - eta_SV(ikl,isn))
	258	eta_SV(ikl,isn) = 0.
	259
	260
	261	! +--Melting if EExcsv > 0
	262	! + ======================
	263
	264	EnMelt = max(zero, EExcsv(ikl) )
	265
	266	! +--Energy Consumption
	267	! + ^^^^^^^^^^^^^^^^^^
	268	SnHLat = ro__SV(ikl,isn) * Lf_H2O
	269	dzMelt(ikl) = EnMelt / max(SnHLat, epsi )
	270	noSnow(ikl) = noSnow(ikl) &
	271	+ max(zero ,sign(unun,dzMelt(ikl) & !
	272	-dzsnSV(ikl ,isn))) & ! 1 if full Melt
	273	*min(1 , max(0 ,1+isnoSV(ikl)-isn)) ! 1 in the Pack
	274	dzMelt(ikl) = &
	275	min(dzsnSV(ikl, isn),dzMelt(ikl))
	276	dzsnSV(ikl,isn) = &
	277	dzsnSV(ikl,isn) -dzMelt(ikl)
	278	zn5_SV(ikl) = zn5_SV(ikl) +dzMelt(ikl)
	279	EExcsv(ikl) = EExcsv(ikl) -dzMelt(ikl)*SnHLat
	280	wem_SV(ikl) = wem_SV(ikl) -dzMelt(ikl)*ro__SV(ikl,isn)
	281
	282	! +--Water Production
	283	! + ^^^^^^^^^^^^^^^^^
	284	drr_SV(ikl) = drr_SV(ikl) &
	285	+ ro__SV(ikl,isn) * dzMelt(ikl)/dt__SV
	286	! #m1 SImelt(ikl) = SImelt(ikl)
	287	! #m1. + ro__SV(ikl,isn) * dzMelt(ikl)
	288	OKmelt =max(zero,sign(unun,drr_SV(ikl)-epsi))
	289
	290	! +--Snow History
	291	! + ^^^^^^^^^^^^
	292	k_face = min( istoSV(ikl,isn),istdSV(1)) & ! = 1 if
	293	*max(0,2-istoSV(ikl,isn) ) ! faceted
	294	istoSV(ikl,isn) = & !
	295	(1.-OKmelt) * istoSV(ikl,isn) & !
	296	+ OKmelt ((1-k_face) istdSV(2) & !
	297	+ k_face * istdSV(3) ) !
	298
	299
	300	! +--Freezing if EExcsv < 0
	301	! + ======================
	302
	303	rdzsno = ro__SV(ikl,isn) * dzsnSV(ikl ,isn)
	304	LayrOK = min( 1, max(0 , isnoSV(ikl)-isn+1))
	305	EnFrez = min(zero, EExcsv(ikl))
	306	WaFrez = -( EnFrez * LayrOK / Lf_H2O)
	307	drrNEW = max(zero,drr_SV(ikl) - WaFrez / dt__SV)
	308	WaFrez = ( drr_SV(ikl) - drrNEW)* dt__SV
	309	drr_SV(ikl) = drrNEW
	310	EExcsv(ikl) = EExcsv(ikl) + WaFrez * Lf_H2O
	311	EnFrez = min(zero,EExcsv(ikl)) * LayrOK
	312	rdzNEW = WaFrez + rdzsno
	313	ro__SV(ikl,isn) = rdzNEW /max(epsi, dzsnSV(ikl,isn))
	314	TsisSV(ikl,isn) = TfSnow &
	315	+ EnFrez /(Cn_dSV *max(epsi, rdzNEW) )
	316	EExcsv(ikl) = EExcsv(ikl) - EnFrez
	317	wer_SV(ikl) = WaFrez &
	318	+ wer_SV(ikl)
	319
	320
	321
	322	! +--Snow Water Content
	323	! + ==================
	324
	325	! +--Percolation Velocity
	326	! + ^^^^^^^^^^^^^^^^^^^^
	327	! #PW SGDiam = 1.6d-4
	328	! #PW. + 1.1d-13 (ro__SV(ikl,isn)ro__SV(ikl,isn)
	329	! #PW. ro__SV(ikl,isn)ro__SV(ikl,isn))
	330
	331	! +--Pore Volume [-]
	332	! + ^^^^^^^^^^^^^^^^^
	333	rosDry =(1. - eta_SV(ikl,isn))* ro__SV(ikl,isn) !
	334	PorVol = 1. - rosDry / ro_Ice !
	335	PorVol = max(PorVol , zero ) !
	336
	337	! +--Water Retention
	338	! + ^^^^^^^^^^^^^^^^
	339	rWater = ws0dSV * PorVol * ro_Wat * dzsnSV(ikl,isn)
	340	drrNEW = max(zero,drr_SV(ikl) - rWater /dt__SV)
	341	rWater = ( drr_SV(ikl) - drrNEW)*dt__SV
	342	drr_SV(ikl) = drrNEW
	343	rdzNEW = rWater &
	344	+ rosDry * dzsnSV(ikl,isn)
	345	eta_SV(ikl,isn) = rWater / max(epsi,rdzNEW)
	346	ro__SV(ikl,isn) = rdzNEW / max(epsi,dzsnSV(ikl,isn))
	347
	348	! +--Pore Hole Close OFF
	349	! + ^^^^^^^^^^^^^^^^^^^
	350	PClose = max(zero, &
	351	sign(unun,ro__SV(ikl,isn) &
	352	-roCdSV ))
	353	ispiSV(ikl) = ispiSV(ikl) *(1.-PClose) &
	354	+ max(ispiSV(ikl),isn) * Pclose
	355	PClose = max(0 , & ! Water under SuPer.Ice
	356	min (1 ,ispiSV(ikl) & ! contributes to
	357	-isn )) ! Surficial Water
	358
	359	!XF
	360	if(ro__SV(ikl,isn) >= roCdSV.and.ro__SV(ikl,1)<900) &
	361	PClose = min(0.50,PClose * &
	362	(1.-(ro_ice-ro__SV(ikl,isn))/(ro_ice-roCdSV)))
	363
	364	PClose = max(0.,min(1.,PClose))
	365
	366	if(isn==1) then
	367	PClose = 1
	368	ispiSV(ikl)= max(ispiSV(ikl),1)
	369	endif
	370
	371	if(drr_SV(ikl) >0 .and.TsisSV(ikl,isn)>273.14) then
	372	if((ro__SV(ikl,isn)>900.and.ro__SV(ikl,isn)<920).or. &
	373	ro__SV(ikl,isn)>950) then
	374	dzsnSV(ikl,isn) = dzsnSV(ikl,isn)*ro__SV(ikl,isn)/ro_ice
	375	ro__SV(ikl,isn) = ro_ice
	376	PClose = 1
	377	endif
	378	endif
	379
	380	! if (isn>1.and.isn<nsno .and.
	381	! . ro__SV(ikl,isn-1)>900 .and.
	382	! . ro__SV(ikl,isn) >roCdSV .and.
	383	! . ro__SV(ikl,isn) <900 .and.
	384	! . TsisSV(ikl,isn) >273.14 .and.
	385	! . TsisSV(ikl,isn+1)<273.15 .and.
	386	! . drr_SV(ikl) >0) then
	387	! TsisSV(ikl,isn)=273.14
	388	! PClose = 1
	389	! endif
	390
	391	!XF
	392	rusnSV(ikl) = rusnSV(ikl) &
	393	+ drr_SV(ikl) dt__SV PClose
	394	rusnSV0(ikl)= rusnSV0(ikl) &
	395	+ drr_SV(ikl) dt__SV PClose
	396	drr_SV(ikl) = drr_SV(ikl) *(1.-PClose)
	397
	398	END DO
	399
	400	END DO
	401
	402
	403	! +--Remove Zero-Thickness Layers
	404	! + ============================
	405
	406	1000 CONTINUE
	407	isnitr = 0
	408	DO ikl=1,knonv
	409	isnUpD = 0
	410	isinew = 0
	411	!XF
	412
	413
	414	DO isn=1,min(nsno-1,isnoSV(ikl))
	415	isnnew =(unun-max(zero ,sign(unun,dzsnSV(ikl,isn)-dzepsi))) &
	416	* max(0 , min(1 ,isnoSV(ikl) +1 -isn ))
	417	isnUpD = max(isnUpD, isnnew)
	418	isnitr = max(isnitr, isnnew)
	419	isinew = isnisnUpD max(0, 1-isinew) & ! LowerMost 0-Layer
	420	+isinew ! Index
	421	dzsnSV(ikl,isn) = dzsnSV(ikl,isn+isnnew)
	422	ro__SV(ikl,isn) = ro__SV(ikl,isn+isnnew)
	423	TsisSV(ikl,isn) = TsisSV(ikl,isn+isnnew)
	424	eta_SV(ikl,isn) = eta_SV(ikl,isn+isnnew)
	425	G1snSV(ikl,isn) = G1snSV(ikl,isn+isnnew)
	426	G2snSV(ikl,isn) = G2snSV(ikl,isn+isnnew)
	427	dzsnSV(ikl,isn+isnnew) =(1-isnnew)*dzsnSV(ikl,isn+isnnew)
	428	ro__SV(ikl,isn+isnnew) =(1-isnnew)*ro__SV(ikl,isn+isnnew)
	429	eta_SV(ikl,isn+isnnew) =(1-isnnew)*eta_SV(ikl,isn+isnnew)
	430	G1snSV(ikl,isn+isnnew) =(1-isnnew)*G1snSV(ikl,isn+isnnew)
	431	G2snSV(ikl,isn+isnnew) =(1-isnnew)*G2snSV(ikl,isn+isnnew)
	432
	433	END DO
	434	isnoSV(ikl) = isnoSV(ikl)-isnUpD ! Nb of Snow Layer
	435	ispiSV(ikl) = ispiSV(ikl) & ! Nb of SuperI Layer
	436	-isnUpD *max(0,min(ispiSV(ikl)-isinew,1)) ! Update if I=0
	437
	438	END DO
	439
	440	IF (isnitr.GT.0) GO TO 1000
	441
	442
	443	! +--New upper Limit of the non erodible Snow (istoSV .GT. 1)
	444	! + ========================================
	445
	446	DO ikl=1,knonv
	447	nh = 0
	448	!XF
	449	DO isn= isnoSV(ikl),1,-1
	450	nh = nh + isn* min(istoSV(ikl,isn)-1,1)*max(0,1-nh)
	451	ENDDO
	452	zc = 0.
	453	zt = 0.
	454	!XF
	455	DO isn=1,isnoSV(ikl)
	456	zc = zc + dzsnSV(ikl,isn) *ro__SV(ikl,isn) &
	457	* max(0,min(1,nh+1-isn))
	458	zt = zt + dzsnSV(ikl,isn) *ro__SV(ikl,isn)
	459	END DO
	460	zWE_SV(ikl) = zt
	461	zWEcSV(ikl) = min(zWEcSV(ikl),zt)
	462	zWEcSV(ikl) = max(zWEcSV(ikl),zc)
	463	END DO
	464
	465
	466	! +--Energy Budget (OUT)
	467	! + ===================
	468
	469	! #vm DO ikl=1,knonv
	470	! #vm EqSn01(ikl) =-EqSn_0(ikl)
	471	! #vm. -EExcsv(ikl)
	472	! #vm END DO
	473	! #vm DO isn=nsno,1,-1
	474	! #vm DO ikl=1,knonv
	475	! #vm EqSn01(ikl) = EqSn01(ikl) + ro__SV(ikl,isn) *dzsnSV(ikl,isn)
	476	! #vm. (Cn_dSV (TsisSV(ikl,isn) -TfSnow )
	477	! #vm. -Lf_H2O *(1. -eta_SV(ikl,isn)))
	478	! #vm END DO
	479	! #vm END DO
	480
	481
	482	! +--"Negative Heat" from supercooled rain
	483	! + ------------------------------------
	484
	485	DO ikl=1,knonv
	486	EExcsv(ikl) = EExcsv(ikl) + EExdum(ikl)
	487
	488
	489	! +--Surficial Water Run OFF
	490	! + -----------------------
	491
	492	rusnew = rusnSV(ikl) * SWf_SV(ikl)
	493
	494	if(isnoSV(ikl)<=1 .OR. opt_runoff_ac) rusnew = 0.
	495	! !if(ivgtSV(ikl)>=1) rusnew = 0.
	496
	497	! #EU rusnew = 0.
	498	! #AC rusnew = 0.
	499
	500	RnofSV(ikl) = RnofSV(ikl) &
	501	+(rusnSV(ikl) - rusnew ) / dt__SV
	502	RuofSV(ikl,1) = RuofSV(ikl,1) &
	503	+(rusnSV(ikl) - rusnew ) / dt__SV
	504	RuofSV(ikl,4) = RuofSV(ikl,4) &
	505	+(rusnSV0(ikl) ) / dt__SV
	506	rusnSV(ikl) = rusnew
	507	END DO
	508
	509
	510	! +--Percolation down the Continental Ice Pack
	511	! + -----------------------------------------
	512
	513	DO ikl=1,knonv
	514	drr_SV(ikl) = drr_SV(ikl) + rusnSV(ikl) &
	515	* (1-min(1,ispiSV(ikl)))/ dt__SV
	516	rusnSV(ikl) = rusnSV(ikl) &
	517	* min(1,ispiSV(ikl))
	518	END DO
	519
	520	!XF removal of too thin snowlayers if TT> 275.15 + bug if TT>> 273.15
	521	DO ikl=1,knonv
	522	zt=0.
	523	DO isn=1,isnoSV(ikl)
	524	zt=zt+dzsnSV(ikl,isn)
	525	ENDDO
	526
	527	if(zt<0.005+(TaT_SV(ikl)-TfSnow)/1000..and. &
	528	isnoSV(ikl) >0 .and. &
	529	TaT_SV(ikl) >=TfSnow .and. &
	530	istoSV(ikl,isnoSV(ikl)) >1 ) then
	531	DO isn=1,isnoSV(ikl)
	532	drr_SV(ikl) = drr_SV(ikl) &
	533	+ dzsnSV(ikl,isn)*ro__SV(ikl,isn) /dt__SV
	534	dzsnSV(ikl,isn)= 0.
	535
	536	ENDDO
	537	isnoSV(ikl) = 0
	538	endif
	539	ENDDO
	540
	541	! +--Slush Formation (Activated. CAUTION: ADD RunOff Possibility before Activation)
	542	! + --------------- ^^^^^^^ ^^^
	543
	544	IF (is_ok_slush) THEN
	545
	546	DO ikl=1,knonv
	547	DO isn=1,isnoSV(ikl)
	548	kSlush = min(1,max(0,isn+1-ispiSV(ikl))) ! Slush Switch
	549
	550	! +--Available Additional Pore Volume [-]
	551	! + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
	552	PorVol = 1. - ro__SV(ikl,isn) & ! [--]
	553	*(1. - eta_SV(ikl,isn))/ ro_Ice & !
	554	- eta_SV(ikl,isn) & !
	555	*ro__SV(ikl,isn) / ro_Wat !
	556	PorVol = max(PorVol , zero ) !
	557	zWater = dzsnSV(ikl,isn) * PorVol * 1000. & ! [mm] OR [kg/m2]
	558	* (1. -SWS_SV(ikl) & ! 0 <=> freezing
	559	*(1 -min(1,iabs(isn-isnoSV(ikl))))) ! 1 <=> isn=isnoSV
	560	zSlush = min(rusnSV(ikl) , zWater) ! [mm] OR [kg/m2]
	561	ro_new =(dzsnSV(ikl,isn) * ro__SV(ikl,isn) & !
	562	+zSlush ) & !
	563	/ max(dzsnSV(ikl,isn) , epsi ) !
	564	if(ro_new<ro_Ice+20) then ! MAX 940kg/m3 !
	565	rusnSV(ikl) = rusnSV(ikl) - zSlush ! [mm] OR [kg/m2]
	566	RuofSV(ikl,4)= max(0.,RuofSV(ikl,4) - zSlush/dt__SV)
	567	eta_SV(ikl,isn) =(ro_new - ro__SV(ikl,isn) & !
	568	*(1. - eta_SV(ikl,isn))) & !
	569	/ max (ro_new , epsi ) !
	570	ro__SV(ikl,isn) = ro_new !
	571	endif
	572	END DO
	573	END DO
	574	END IF
	575
	576	! +--Impact of the Sublimation/Deposition on the Surface Mass Balance
	577	! + ================================================================
	578
	579	DO ikl=1,knonv
	580	isn = isnoSV(ikl)
	581	dzVap0 = dt__SV &
	582	* HLs_sv(ikl) * min(isn , 1 ) &
	583	/(Lx_H2O(ikl) * max(ro__SV(ikl,isn) , epsi))
	584	NOLayr=min(zero,sign(unun,dzsnSV(ikl,isn) + dzVap0))
	585	dzVap1=min(zero, dzsnSV(ikl,isn) + dzVap0)
	586
	587
	588	! +--Additional Energy
	589	! + -----------------
	590
	591	! #VH AdEnrg = dzVap0 * ro__SV(ikl,isnoSV(ikl)) ! Water Vapor
	592	! #VH. C__Wat (TsisSV(ikl,isnoSV(ikl)) -TfSnow) ! Sensible Heat
	593
	594	! #aH B_Enrg =(Cn_dSV *(TsisSV(ikl,isn) -TfSnow )
	595	! #aH. -Lf_H2O *(1. -eta_SV(ikl,isn)))
	596	! #aH. /(1. + dzVap0 /max(epsi,dzsnSV(ikl,isn)))
	597	! #aH eta_SV(ikl,isn) =
	598	! #aH. max(zero,unun +(B_Enrg
	599	! #aH. -(TsisSV(ikl,isn) -TfSnow)*Cn_dSV)
	600	! #aH. /Lf_H2O )
	601	! #aH TsisSV(ikl,isn) = ( B_Enrg
	602	! #aH. +(1. -eta_SV(ikl,isn))
	603	! #aH. *Lf_H2O )
	604	! #aH. / Cn_dSV
	605	! #aH. + TfSnow
	606
	607	! #e1 STOP "PLEASE add Energy (#aH) from deposition/sublimation"
	608
	609
	610	! +--Update of the upper Snow layer Thickness
	611	! + ----------------------------------------
	612
	613	dzsnSV(ikl,isn) = &
	614	max(zero, dzsnSV(ikl,isnoSV(ikl)) + dzVap0)
	615	isnoSV(ikl) = isnoSV(ikl) + NOLayr
	616	isn = isnoSV(ikl)
	617	dzsnSV(ikl,isn) = dzsnSV(ikl,isn) + dzVap1
	618	wes_SV(ikl) = ro__SV(ikl,isn) * dzVap0
	619
	620	END DO
	621
	622
	623	! +--Energy Budget (OUT)
	624	! + ===================
	625
	626	! #vm DO ikl=1,knonv
	627	! #vm EqSn02(ikl) =-EqSn_0(ikl)
	628	! #vm. -EExcsv(ikl)
	629	! #vm END DO
	630	! #vm DO isn=nsno,1,-1
	631	! #vm DO ikl=1,knonv
	632	! #vm EqSn02(ikl) = EqSn01(ikl) + ro__SV(ikl,isn) *dzsnSV(ikl,isn)
	633	! #vm. (Cn_dSV (TsisSV(ikl,isn) -TfSnow )
	634	! #vm. -Lf_H2O *(1. -eta_SV(ikl,isn)))
	635	! #vm END DO
	636	! #vm END DO
	637
	638
	639	! +--Snow/I Budget
	640	! + -------------
	641
	642	! #m1 DO ikl=1,knonv
	643	! #m1 SIsubl(ikl) = dt__SVHLs_sv(ikl)min(isnoSV(ikl),1)
	644	! #m1. /Lx_H2O(ikl)
	645	! #m1 SIrnof(ikl) = rusnSV(ikl) + RnofSV(ikl) * dt__SV
	646	! #m1. - SIrnof(ikl)
	647	! #m1 END DO
	648
	649
	650	! +--Anticipated Disappearance of a rapidly Melting too thin Last Snow Layer
	651	! + =======================================================================
	652
	653	DO ikl=1,knonv
	654	LastOK = min(1 , max(0 ,iiceSV(ikl)-isnoSV(ikl)+2) &
	655	*min(1 ,isnoSV(ikl)-iiceSV(ikl)) &
	656	+min(1 ,isnoSV(ikl)) )
	657	RapdOK = max(zero,sign(unun,dzMelt(ikl)-epsi ))
	658	ThinOK = max(zero,sign(unun,dz_Min -dzsnSV(ikl,1)))
	659	z_Melt = LastOK RapdOKThinOK
	660	noSnow(ikl) = noSnow(ikl) + z_Melt
	661	z_Melt = z_Melt *dzsnSV(ikl,1)
	662	dzsnSV(ikl,1) = dzsnSV(ikl,1) - z_Melt
	663	EExcsv(ikl) = EExcsv(ikl) - z_Melt *ro__SV(ikl,1) &
	664	(1. -eta_SV(ikl,1))Lf_H2O
	665
	666	! +--Water Production
	667	! + ^^^^^^^^^^^^^^^^^
	668	drr_SV(ikl) = drr_SV(ikl) &
	669	+ ro__SV(ikl,1) * z_Melt /dt__SV
	670	END DO
	671
	672
	673	! +--Update Nb of Layers
	674	! + ===================
	675
	676	DO ikl=1,knonv
	677	isnoSV(ikl) = isnoSV(ikl) &
	678	* min(1,iabs(isnoSV(ikl)-noSnow(ikl)))
	679	END DO
	680
	681
	682	! Energy Budget (OUT)
	683	! ===================
	684
	685	! #e1 DO ikl=1,knonv
	686	! #e1 EqSn_1(ikl) = 0.
	687	! #e1 END DO
	688	! #e1 DO isn=nsno,1,-1
	689	! #e1 DO ikl=1,knonv
	690	! #e1 EqSn_1(ikl) = EqSn_1(ikl) + ro__SV(ikl,isn) *dzsnSV(ikl,isn)
	691	! #e1. (Cn_dSV (TsisSV(ikl,isn) -TfSnow )
	692	! #e1. -Lf_H2O *(1. -eta_SV(ikl,isn)))
	693	! #e1 END DO
	694	! #e1 END DO
	695
	696
	697	! +--Water Budget (OUT)
	698	! + ===================
	699
	700	! #vm DO ikl=1,knonv
	701	! #vm WqSn_0(ikl) = WqSn_0(ikl)
	702	! #vm. + HLs_sv(ikl) * dt__SV
	703	! #vm. *min(isnoSV(ikl),1) / Lx_H2O(ikl)
	704	! #vm WqSn_1(ikl) = drr_SV(ikl) * dt__SV
	705	! #vm. + rusnSV(ikl)
	706	! #vm. + RnofSV(ikl) * dt__SV
	707	! #vm END DO
	708	! #vm DO isn=nsno,1,-1
	709	! #vm DO ikl=1,knonv
	710	! #vm WqSn_1(ikl) = WqSn_1(ikl)
	711	! #vm. + ro__SV(ikl,isn)* dzsnSV(ikl,isn)
	712	! #vm END DO
	713	! #vm END DO
	714
	715
	716	return
	717	end subroutine sisvat_qsn

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

Download in other formats:

Original Format