Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

sisvat_qsn.f90 @ 5452

Last change on this file since 5452 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

Line
1
2
3	subroutine SISVAT_qSn &
4	( &
5	! #e1. EqSn_0,EqSn_1,EqSn_d
6	! #m1. ,SIsubl,SImelt,SIrnof
7	)
8
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	! +------------------------------------------------------------------------+
50
51
52
53
54	! +--Global Variables
55	! + ================
56
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
64
65
66	IMPLICIT NONE
67
68
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
83
84
85	! +--Internal Variables
86	! + ==================
87
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
127
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