Context Navigation

sisvat_qsn.F @ 3851

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

Line
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