Context Navigation

sisvat_qsn.F @ 4490

Last change on this file since 4490 was 3900, checked in by evignon, 4 years ago
Commission de la nouvelle interface LMDZ-SISVAT la prochaine commission consistera a supprimer l'ancien repertoire sisvat et a faire un peu de nettoyage.
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	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	C +--Internal Variables
86	C + ==================
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	C +--OUTPUT of SISVAT Trace Statistics (see assignation in PHY_SISVAT)
129	C + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
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	C +--Energy and Mass Budget
138	C + ~~~~~~~~~~~~~~~~~~~~~~
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	C +--Slush Diagnostic: IO
149	C + ~~~~~~~~~~~~~~~~~~~~
150	! #vu logical su_opn ! IO Switch
151	! #vu common/SI_qSn_L/su_opn !
152
153
154	C +--DATA
155	C + ====
156
157	data dzepsi/0.0001/ ! Minim. Snow Layer Thickness (!)
158	c #?? data dz_Min/0.005/ ! Minim. Snow Layer Thickness
159	c ... 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	C +--Energy Budget (IN)
165	C + ==================
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	C +--Water Budget (IN)
180	C + ==================
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	C +--Snow Melt Budget
194	C + ================
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	C +--Initialization
203	C + ==============
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	C +--Melting/Freezing Energy
217	C + =======================
218
219	C +...REMARK: Snow liquid Water Temperature assumed = TfSnow
220	C + ^^^^^^
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	C +--Surficial Water Status
231	C + ----------------------
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	C +--Energy, store Previous Content
243	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
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	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	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	C +--Snow Water Content
323	C + ==================
324
325	C +--Percolation Velocity
326	C + ^^^^^^^^^^^^^^^^^^^^
327	c #PW SGDiam = 1.6d-4
328	c #PW. + 1.1d-13 (ro__SV(ikl,isn)ro__SV(ikl,isn)
329	c #PW. ro__SV(ikl,isn)ro__SV(ikl,isn))
330
331	C +--Pore Volume [-]
332	C + ^^^^^^^^^^^^^^^^^
333	rosDry =(1. - eta_SV(ikl,isn))* ro__SV(ikl,isn) !
334	PorVol = 1. - rosDry / ro_Ice !
335	PorVol = max(PorVol , zero ) !
336
337	C +--Water Retention
338	C + ^^^^^^^^^^^^^^^^
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	C +--Pore Hole Close OFF
349	C + ^^^^^^^^^^^^^^^^^^^
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	cXF
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	c if (isn>1.and.isn<nsno .and.
381	c . ro__SV(ikl,isn-1)>900 .and.
382	c . ro__SV(ikl,isn) >roCdSV .and.
383	c . ro__SV(ikl,isn) <900 .and.
384	c . TsisSV(ikl,isn) >273.14 .and.
385	c . TsisSV(ikl,isn+1)<273.15 .and.
386	c . drr_SV(ikl) >0) then
387	c TsisSV(ikl,isn)=273.14
388	c PClose = 1
389	c endif
390
391	cXF
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	C +--Remove Zero-Thickness Layers
404	C + ============================
405
406	1000 CONTINUE
407	isnitr = 0
408	DO ikl=1,knonv
409	isnUpD = 0
410	isinew = 0
411	cXF
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	C +--New upper Limit of the non erodible Snow (istoSV .GT. 1)
444	C + ========================================
445
446	DO ikl=1,knonv
447	nh = 0
448	cXF
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	cXF
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	C +--Energy Budget (OUT)
467	C + ===================
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	C +--"Negative Heat" from supercooled rain
483	C + ------------------------------------
484
485	DO ikl=1,knonv
486	EExcsv(ikl) = EExcsv(ikl) + EExdum(ikl)
487
488
489	C +--Surficial Water Run OFF
490	C + -----------------------
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	c #EU rusnew = 0.
498	c #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	C +--Percolation down the Continental Ice Pack
511	C + -----------------------------------------
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	cXF 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	C +--Slush Formation (Activated. CAUTION: ADD RunOff Possibility before Activation)
542	C + --------------- ^^^^^^^ ^^^
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	C +--Available Additional Pore Volume [-]
551	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
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	C +--Impact of the Sublimation/Deposition on the Surface Mass Balance
577	C + ================================================================
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	C +--Additional Energy
589	C + -----------------
590
591	c #VH AdEnrg = dzVap0 * ro__SV(ikl,isnoSV(ikl)) ! Water Vapor
592	c #VH. C__Wat (TsisSV(ikl,isnoSV(ikl)) -TfSnow) ! Sensible Heat
593
594	c #aH B_Enrg =(Cn_dSV *(TsisSV(ikl,isn) -TfSnow )
595	c #aH. -Lf_H2O *(1. -eta_SV(ikl,isn)))
596	c #aH. /(1. + dzVap0 /max(epsi,dzsnSV(ikl,isn)))
597	c #aH eta_SV(ikl,isn) =
598	c #aH. max(zero,unun +(B_Enrg
599	c #aH. -(TsisSV(ikl,isn) -TfSnow)*Cn_dSV)
600	c #aH. /Lf_H2O )
601	c #aH TsisSV(ikl,isn) = ( B_Enrg
602	c #aH. +(1. -eta_SV(ikl,isn))
603	c #aH. *Lf_H2O )
604	c #aH. / Cn_dSV
605	c #aH. + TfSnow
606
607	! #e1 STOP "PLEASE add Energy (#aH) from deposition/sublimation"
608
609
610	C +--Update of the upper Snow layer Thickness
611	C + ----------------------------------------
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	C +--Energy Budget (OUT)
624	C + ===================
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	C +--Snow/I Budget
640	C + -------------
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	C +--Anticipated Disappearance of a rapidly Melting too thin Last Snow Layer
651	C + =======================================================================
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	C +--Water Production
667	C + ^^^^^^^^^^^^^^^^^
668	drr_SV(ikl) = drr_SV(ikl)
669	. + ro__SV(ikl,1) * z_Melt /dt__SV
670	END DO
671
672
673	C +--Update Nb of Layers
674	C + ===================
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	C +--Water Budget (OUT)
698	C + ===================
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

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

Download in other formats:

Original Format