Context Navigation

sisvat_qsn.F @ 3817

Last change on this file since 3817 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: 27.8 KB

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

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

Download in other formats:

Original Format