Context Navigation

sisvat_zsn.F @ 3991

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

Line
1
2
3	subroutine SISVAT_zSn
4
5	C +------------------------------------------------------------------------+
6	C \| MAR SISVAT_zSn 12-07-2019 MAR \|
7	C \| SubRoutine SISVAT_zSn manages the Snow Pack vertical Discretization \|
8	C \| \|
9	C +------------------------------------------------------------------------+
10	C \| \|
11	C \| PARAMETERS: knonv: Total Number of columns = \|
12	C \| ^^^^^^^^^^ = Total Number of continental grid boxes \|
13	C \| X Number of Mosaic Cell per grid box \|
14	C \| \|
15	C \| INPUT / NLaysv = New Snow Layer Switch \|
16	C \| OUTPUT: isnoSV = total Nb of Ice/Snow Layers \|
17	C \| ^^^^^^ ispiSV = 0,...,nsno: Uppermost Superimposed Ice Layer \|
18	C \| iiceSV = total Nb of Ice Layers \|
19	C \| istoSV = 0,...,5 : Snow History (see istdSV data) \|
20	C \| \|
21	C \| INPUT / TsisSV : Soil/Ice Temperatures (layers -nsol,-nsol+1,..,0)\|
22	C \| OUTPUT: & Snow Temperatures (layers 1,2,...,nsno) [K] \|
23	C \| ^^^^^^ ro__SV : Soil/Snow Volumic Mass [kg/m3] \|
24	C \| eta_SV : Soil/Snow Water Content [m3/m3] \|
25	C \| dzsnSV : Snow Layer Thickness [m] \|
26	C \| G1snSV : Dendricity (<0) or Sphericity (>0) of Snow Layer \|
27	C \| G2snSV : Sphericity (>0) or Size of Snow Layer \|
28	C \| agsnSV : Snow Age [day] \|
29	C \| \|
30	C \| METHOD: 1) Agregate the thinest Snow Layer \|
31	C \| ^^^^^^ if a new Snow Layer has been precipitated (NLaysv = 1) \|
32	C \| 2) Divide a too thick Snow Layer except \|
33	C \| if the maximum Number of Layer is reached \|
34	C \| in this case forces NLay_s = 1 \|
35	C \| 3) Agregate the thinest Snow Layer \|
36	C \| in order to divide a too thick Snow Layer \|
37	C \| at next Time Step when NLay_s = 1 \|
38	C \| \|
39	C +------------------------------------------------------------------------+
40
41
42
43
44	C +--Global Variables
45	C + ================
46
47
48	use VARphy
49	use VAR_SV
50	use VARdSV
51	use VAR0SV
52	use VARxSV
53	use VARySV
54
55	IMPLICIT NONE
56
57
58	C +--Internal Variables
59	C + ==================
60
61	integer ikl ,isn ,i !
62
63
64
65
66
67
68	integer NLay_s(knonv) ! Split Snow Layer Switch
69	integer isagr1(knonv) ! 1st Layer History
70	integer isagr2(knonv) ! 2nd Layer History
71	integer LstLay ! 0 ====> isnoSV = 1
72	integer isno_n ! Snow Normal.Profile
73	integer iice_n ! Ice Normal.Profile
74	integer iiceOK ! Ice Switch
75	integer icemix ! 0 ====> Agregated Snow+Ice=Snow
76	C + ! 1 Ice
77	integer isn1 (knonv) ! 1st layer to stagger
78	real staggr ! stagger Switch
79
80	real WEagre(knonv) ! Snow Water Equivalent Thickness
81	real dzthin(knonv) ! Thickness of the thinest layer
82	real OKthin ! Swich ON a new thinest layer
83	real dz_dif ! difference from ideal discret.
84	real thickL ! Thick Layer Indicator
85	real OK_ICE ! Swich ON uppermost Ice Layer
86
87	real Agrege(knonv) ! 1. when Agregation constrained
88	real dzepsi ! Min Single Snw Layer Thickness
89	real dzxmin ! Min Acceptable Layer Thickness
90	real dz_min ! Min Layer Thickness
91	real dz_max ! Max Layer Thickness
92	real dzagr1(knonv) ! 1st Layer Thickness
93	real dzagr2(knonv) ! 2nd Layer Thickness
94	real T_agr1(knonv) ! 1st Layer Temperature
95	real T_agr2(knonv) ! 2nd Layer Temperature
96	real roagr1(knonv) ! 1st Layer Density
97	real roagr2(knonv) ! 2nd Layer Density
98	real etagr1(knonv) ! 1st Layer Water Content
99	real etagr2(knonv) ! 2nd Layer Water Content
100	real G1agr1(knonv) ! 1st Layer Dendricity/Spher.
101	real G1agr2(knonv) ! 2nd Layer Dendricity/Spher.
102	real G2agr1(knonv) ! 1st Layer Sphericity/Size
103	real G2agr2(knonv) ! 2nd Layer Sphericity/Size
104	real agagr1(knonv) ! 1st Layer Age
105	real agagr2(knonv) ! 2nd Layer Age
106
107
108	C +--DATA
109	C + ====
110
111	data icemix / 0 / ! 0 ====> Agregated Snow+Ice=Snow
112	data dzepsi / 0.0020/ ! Min single Layer Thickness
113	data dzxmin / 0.0025/ ! Min accept.Layer Thickness
114	c #EU data dz_min / 0.0050/ ! Min Local Layer Thickness < SMn
115	data dz_min / 0.0040/ ! Min Local Layer Thickness < SMn
116	data dz_max / 0.0300/ ! Min Gener. Layer Thickness
117	C + CAUTION: dz_max > dz_min*2 is required ! Otherwise re-agregation is
118	C + ! activated after splitting
119
120
121
122
123
124	C +--Constrains Agregation of too thin Layers
125	C + =================================================
126
127	C +--Search the thinest non-zero Layer
128	C + ----------------------------------
129
130	DO ikl=1,knonv
131	if(isnoSV(ikl)<=2) dz_min=max(0.0050,dz_min)
132
133	dzepsi=0.0015
134	if(ro__SV(ikl,isnoSV(ikl))>920) dzepsi=0.0020
135
136	dzthin(ikl) = 0. ! Arbitrary unrealistic
137	END DO ! Layer Thickness
138	cXF
139	DO ikl=1,knonv
140	DO isn=1,isnoSV(ikl)-3 ! no agregation of 3 first snowlayers
141	! XF 04/07/2019
142
143	isno_n = isnoSV(ikl)-isn+1 ! Snow Normal.Profile
144	iice_n = iiceSV(ikl)-isn ! Ice Normal.Profile
145	iiceOK = min(1,max(0,iice_n +1)) ! Ice Switch
146	! #vz dz_ref(isn) = !
147	! #vz. dz_min ((1-iiceOK)isno_n*isno_n ! Theoretical Profile
148	! #vz. + iiceOK * 2**iice_n) !
149	! #vz. /max(1,isnoSV(ikl)) !
150	dz_dif = max(zero, ! Actual Profile
151	. dz_min !
152	. ((1-iiceOK)isno_n*isno_n ! Theoretical Profile
153	. + iiceOK 2. *iice_n) !
154	. - dzsnSV(ikl, isn) ) ! Actual Profile
155	! #vz dzwdif(isn) = dz_dif !
156	OKthin = max(zero, !
157	. sign(unun, !
158	. dz_dif-dzthin(ikl))) ! 1.=> New thinest Lay.
159	. * max(0, ! 1 => .le. isnoSV
160	. min(1, ! 1 => isn is in the
161	. isnoSV(ikl)-isn +1 )) ! Snow Pack
162	. * min(unun, !
163	!
164	! 1st additional Condition to accept OKthin
165	. max(zero, ! combination
166	. sign(unun,G1snSV(ikl, isn ) ! G1 with same
167	. *G1snSV(ikl,max(1,isn-1)))) ! sign => OK
168	!
169	! 2nd additional Condition to accept OKthin
170	. + max(zero, ! G1>0
171	. sign(unun,G1snSV(ikl, isn ))) ! =>OK
172	!
173	! 3rd additional Condition to accept OKthin
174	. + max(zero, ! dz too small
175	. sign(unun,dzxmin ! =>OK
176	. -dzsnSV(ikl, isn ))))!
177
178	i_thin(ikl) = (1. - OKthin) * i_thin(ikl) ! Update thinest Lay.
179	. + OKthin * isn ! Index
180	dzthin(ikl) = (1. - OKthin) * dzthin(ikl) !
181	. + OKthin * dz_dif !
182	END DO
183	END DO
184
185
186
187	DO ikl=1,knonv
188	DO isn=1,isnoSV(ikl)
189	OKthin = max(zero, !
190	. sign(unun, !
191	. dz_min !
192	. -dzsnSV(ikl,isn))) !
193	. * max(zero, ! ON if dz > 0
194	. sign(unun, !
195	. dzsnSV(ikl,isn)-epsi)) !
196	. *min(1,max(0, ! Multiple Snow Lay.
197	. min (1, ! Switch = 1
198	. isnoSV(ikl) ! if isno > iice + 1
199	. -iiceSV(ikl)-1)) !
200	C + !
201	. +int(max(zero, !
202	. sign(unun, !
203	. dzepsi ! Minimum accepted for
204	. -dzsnSV(ikl,isn)))) ! 1 Snow Layer over Ice
205	. *int(max(zero, ! ON if dz > 0
206	. sign(unun, !
207	. dzsnSV(ikl,isn)-epsi)))!
208	. *(1 -min (abs(isnoSV(ikl) ! Switch = 1
209	. -iiceSV(ikl)-1),1)) ! if isno = iice + 1
210	C + !
211	. +max(0, ! Ice
212	. min (1, ! Switch
213	. iiceSV(ikl)+1-isn))) !
214	. *min(unun, !
215	. max(zero, ! combination
216	. sign(unun,G1snSV(ikl, isn ) ! G1>0 + G1<0
217	. *G1snSV(ikl,max(1,isn-1)))) ! NO
218	. + max(zero, !
219	. sign(unun,G1snSV(ikl, isn ))) !
220	. + max(zero, !
221	. sign(unun,dzxmin !
222	. -dzsnSV(ikl, isn ))))!
223	i_thin(ikl) = (1. - OKthin) * i_thin(ikl) ! Update thinest Lay.
224	. + OKthin * isn ! Index
225	END DO
226	END DO
227
228
229
230
231	C + ***************
232	call SISVAT_zCr
233	C + ***************
234
235
236	C +--Assign the 2 Layers to agregate
237	C + -------------------------------
238
239	DO ikl=1,knonv
240	isn = i_thin(ikl)
241	if(LIndsv(ikl)>0) isn=min(nsno-1,isn) ! cXF
242	isagr1(ikl) = istoSV(ikl,isn)
243	isagr2(ikl) = istoSV(ikl,isn+LIndsv(ikl))
244	dzagr1(ikl) = dzsnSV(ikl,isn)
245	dzagr2(ikl) = dzsnSV(ikl,isn+LIndsv(ikl))
246	T_agr1(ikl) = TsisSV(ikl,isn)
247	T_agr2(ikl) = TsisSV(ikl,isn+LIndsv(ikl))
248	roagr1(ikl) = ro__SV(ikl,isn)
249	roagr2(ikl) = ro__SV(ikl,isn+LIndsv(ikl))
250	etagr1(ikl) = eta_SV(ikl,isn)
251	etagr2(ikl) = eta_SV(ikl,isn+LIndsv(ikl))
252	G1agr1(ikl) = G1snSV(ikl,isn)
253	G1agr2(ikl) = G1snSV(ikl,isn+LIndsv(ikl))
254	G2agr1(ikl) = G2snSV(ikl,isn)
255	G2agr2(ikl) = G2snSV(ikl,isn+LIndsv(ikl))
256	agagr1(ikl) = agsnSV(ikl,isn)
257	agagr2(ikl) = agsnSV(ikl,isn+LIndsv(ikl))
258	LstLay = min(1,max( 0,isnoSV(ikl) -1)) ! 0 if single Layer
259	isnoSV(ikl) = isnoSV(ikl) ! decrement isnoSV
260	. -(1-LstLay)* max(zero, ! if downmost Layer
261	. sign(unun,eps_21 ! < 1.e-21 m
262	. -dzsnSV(ikl,1))) !
263	isnoSV(ikl) = max( 0, isnoSV(ikl) ) !
264	Agrege(ikl) = max(zero, !
265	. sign(unun,dz_min ! No Agregation
266	. -dzagr1(ikl) )) ! if too thick Layer
267	. *LstLay ! if a single Layer
268	. * min( max(0 ,isnoSV(ikl)+1 ! if Agregation
269	. -i_thin(ikl) ! with a Layer
270	. -LIndsv(ikl) ),1) ! above the Pack
271
272	WEagre(ikl) = 0.
273	END DO
274
275
276	DO ikl=1,knonv
277	DO isn=1,isnoSV(ikl)
278	WEagre(ikl) = WEagre(ikl) + ro__SV(ikl,isn)*dzsnSV(ikl,isn)
279	. *min(1,max(0,i_thin(ikl)+1-isn))
280	ENDDO
281	ENDDO
282
283
284	C +--Agregates
285	C + ---------
286
287	C + ***************
288	call SISVAT_zAg
289	. (isagr1,isagr2,WEagre
290	. ,dzagr1,dzagr2,T_agr1,T_agr2
291	. ,roagr1,roagr2,etagr1,etagr2
292	. ,G1agr1,G1agr2,G2agr1,G2agr2
293	. ,agagr1,agagr2,Agrege
294	. )
295	C + ***************
296
297
298	C +--Rearranges the Layers
299	C + ---------------------
300
301	C +--New (agregated) Snow layer
302	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^
303	DO ikl=1,knonv
304	isn = i_thin(ikl)
305	isn = min(isn,isn+LIndsv(ikl))
306	isnoSV(ikl) = max(0.,isnoSV(ikl) -Agrege(ikl))
307	iiceSV(ikl) = iiceSV(ikl)
308	. -max(0,sign(1,iiceSV(ikl) -isn +icemix))
309	. *Agrege(ikl)
310	. *max(0,sign(1,iiceSV(ikl) -1 ))
311	istoSV(ikl,isn) = (1.-Agrege(ikl))*istoSV(ikl,isn)
312	. + Agrege(ikl) *isagr1(ikl)
313	dzsnSV(ikl,isn) = (1.-Agrege(ikl))*dzsnSV(ikl,isn)
314	. + Agrege(ikl) *dzagr1(ikl)
315	TsisSV(ikl,isn) = (1.-Agrege(ikl))*TsisSV(ikl,isn)
316	. + Agrege(ikl) *T_agr1(ikl)
317	ro__SV(ikl,isn) = (1.-Agrege(ikl))*ro__SV(ikl,isn)
318	. + Agrege(ikl) *roagr1(ikl)
319	eta_SV(ikl,isn) = (1.-Agrege(ikl))*eta_SV(ikl,isn)
320	. + Agrege(ikl) *etagr1(ikl)
321	G1snSV(ikl,isn) = (1.-Agrege(ikl))*G1snSV(ikl,isn)
322	. + Agrege(ikl) *G1agr1(ikl)
323	G2snSV(ikl,isn) = (1.-Agrege(ikl))*G2snSV(ikl,isn)
324	. + Agrege(ikl) *G2agr1(ikl)
325	agsnSV(ikl,isn) = (1.-Agrege(ikl))*agsnSV(ikl,isn)
326	. + Agrege(ikl) *agagr1(ikl)
327	END DO
328
329	C +--Above
330	C + ^^^^^
331	DO ikl=1,knonv
332	isn1(ikl)=max(i_thin(ikl),i_thin(ikl)+LIndsv(ikl))
333	END DO
334	DO i= 1,nsno-1
335	DO ikl=1,knonv
336	staggr = min(1,max(0,i +1 -isn1(ikl) ))
337	istoSV(ikl,i) = (1.-staggr )*istoSV(ikl,i )
338	. + staggr((1.-Agrege(ikl))istoSV(ikl,i )
339	. + Agrege(ikl) *istoSV(ikl,i+1))
340	dzsnSV(ikl,i) = (1.-staggr )*dzsnSV(ikl,i )
341	. + staggr((1.-Agrege(ikl))dzsnSV(ikl,i )
342	. + Agrege(ikl) *dzsnSV(ikl,i+1))
343	TsisSV(ikl,i) = (1.-staggr )*TsisSV(ikl,i )
344	. + staggr((1.-Agrege(ikl))TsisSV(ikl,i )
345	. + Agrege(ikl) *TsisSV(ikl,i+1))
346	ro__SV(ikl,i) = (1.-staggr )*ro__SV(ikl,i )
347	. + staggr((1.-Agrege(ikl))ro__SV(ikl,i )
348	. + Agrege(ikl) *ro__SV(ikl,i+1))
349	eta_SV(ikl,i) = (1.-staggr )*eta_SV(ikl,i )
350	. + staggr((1.-Agrege(ikl))eta_SV(ikl,i )
351	. + Agrege(ikl) *eta_SV(ikl,i+1))
352	G1snSV(ikl,i) = (1.-staggr )*G1snSV(ikl,i )
353	. + staggr((1.-Agrege(ikl))G1snSV(ikl,i )
354	. + Agrege(ikl) *G1snSV(ikl,i+1))
355	G2snSV(ikl,i) = (1.-staggr )*G2snSV(ikl,i )
356	. + staggr((1.-Agrege(ikl))G2snSV(ikl,i )
357	. + Agrege(ikl) *G2snSV(ikl,i+1))
358	agsnSV(ikl,i) = (1.-staggr )*agsnSV(ikl,i )
359	. + staggr((1.-Agrege(ikl))agsnSV(ikl,i )
360	. + Agrege(ikl) *agsnSV(ikl,i+1))
361	END DO
362	END DO
363
364	DO ikl=1,knonv
365	isn = min(isnoSV(ikl) +1,nsno)
366	istoSV(ikl,isn) = (1.-Agrege(ikl))*istoSV(ikl,isn)
367	dzsnSV(ikl,isn) = (1.-Agrege(ikl))*dzsnSV(ikl,isn)
368	TsisSV(ikl,isn) = (1.-Agrege(ikl))*TsisSV(ikl,isn)
369	ro__SV(ikl,isn) = (1.-Agrege(ikl))*ro__SV(ikl,isn)
370	eta_SV(ikl,isn) = (1.-Agrege(ikl))*eta_SV(ikl,isn)
371	G1snSV(ikl,isn) = (1.-Agrege(ikl))*G1snSV(ikl,isn)
372	G2snSV(ikl,isn) = (1.-Agrege(ikl))*G2snSV(ikl,isn)
373	agsnSV(ikl,isn) = (1.-Agrege(ikl))*agsnSV(ikl,isn)
374	END DO
375
376
377
378
379	C +--Constrains Splitting of too thick Layers
380	C + =================================================
381
382
383	C +--Search the thickest non-zero Layer
384	C + ----------------------------------
385
386	DO ikl=1,knonv
387	dzthin(ikl) = 0. ! Arbitrary unrealistic
388	END DO
389	DO ikl=1,knonv
390	DO isn=1,isnoSV(ikl)
391	isno_n = isnoSV(ikl)-isn+1 ! Snow Normal.Profile
392	iice_n = iiceSV(ikl)-isn ! Ice Normal.Profile
393	iiceOK = min(1,max(0,iice_n +1)) ! Ice Switch
394	dz_dif =( dzsnSV(ikl,isn) ! Actual Profile
395	. - dz_max ((1-iiceOK)isno_n*isno_n ! Theoretical Profile
396	. + iiceOK 2. *iice_n) ) !
397	. /max(dzsnSV(ikl,isn),epsi) !
398	OKthin = max(zero, !
399	. sign(unun, !
400	. dz_dif-dzthin(ikl))) ! 1.=>New thickest Lay.
401	. * max(0, ! 1 =>.le. isnoSV
402	. min(1, !
403	. isnoSV(ikl)-isn +1 )) !
404	i_thin(ikl) = (1. - OKthin) * i_thin(ikl) ! Update thickest Lay.
405	. + OKthin * isn ! Index
406	dzthin(ikl) = (1. - OKthin) * dzthin(ikl) !
407	. + OKthin * dz_dif !
408	END DO
409
410	isn=max(1,isnoSV(ikl)-3)
411	if(dzsnSV(ikl,isn)>0.30) then ! surface layer > 30cm
412	i_thin(ikl) = isn ! XF 04/07/2019
413	dzthin(ikl) = dzsnSV(ikl,isn)
414	endif
415
416	isn=max(1,isnoSV(ikl)-2)
417	if(dzsnSV(ikl,isn)>0.10) then ! surface layer > 10cm
418	i_thin(ikl) = isn ! XF 04/07/2019
419	dzthin(ikl) = dzsnSV(ikl,isn)
420	endif
421
422	isn=max(1,isnoSV(ikl)-1)
423	if(dzsnSV(ikl,isn)>0.05) then ! surface layer > 5cm
424	i_thin(ikl) = isn ! XF 04/07/2019
425	dzthin(ikl) = dzsnSV(ikl,isn)
426	endif
427
428	isn=max(1,isnoSV(ikl))
429	if(dzsnSV(ikl,isn)>0.02) then ! surface layer > 2cm
430	i_thin(ikl) = isn ! XF 04/07/2019
431	dzthin(ikl) = dzsnSV(ikl,isn)
432	endif
433
434	END DO
435
436	DO ikl=1,knonv
437	ThickL = max(zero, ! 1. => a too thick
438	. sign(unun,dzthin(ikl) ! Layer exists
439	. -epsi )) !
440	. * max(0,1-max(0 , isnoSV(ikl) ! No spliting allowed
441	. -nsno+1 )) ! if isno > nsno - 1
442	Agrege(ikl) = ThickL ! 1. => effective split
443	. * max(0,1-max(0 , NLaysv(ikl) !
444	. +isnoSV(ikl) !
445	. -nsno+1 )) !
446	NLay_s(ikl) = ThickL ! Agregation
447	. * max(0,1-max(0 , NLaysv(ikl) ! to allow Splitting
448	. +isnoSV(ikl) ! at next Time Step
449	. -nsno )) !
450	. -Agrege(ikl) !
451	NLay_s(ikl) = max(0 , NLay_s(ikl)) ! Agregation effective
452	END DO
453
454
455	C +--Rearranges the Layers
456	C + ---------------------
457
458	DO isn=nsno,2,-1
459	DO ikl=1,knonv
460	IF (Agrege(ikl).gt.0..AND.i_thin(ikl).lt.isnoSV(ikl)) THEN
461	staggr = min(1,max(0,isn-i_thin(ikl) -1))
462	. * min(1,max(0, isnoSV(ikl)-isn+2))
463	istoSV(ikl,isn) = staggr * istoSV(ikl ,isn-1)
464	. + (1. - staggr) * istoSV(ikl ,isn )
465	dzsnSV(ikl,isn) = staggr * dzsnSV(ikl ,isn-1)
466	. + (1. - staggr) * dzsnSV(ikl ,isn )
467	TsisSV(ikl,isn) = staggr * TsisSV(ikl ,isn-1)
468	. + (1. - staggr) * TsisSV(ikl ,isn )
469	ro__SV(ikl,isn) = staggr * ro__SV(ikl ,isn-1)
470	. + (1. - staggr) * ro__SV(ikl ,isn )
471	eta_SV(ikl,isn) = staggr * eta_SV(ikl ,isn-1)
472	. + (1. - staggr) * eta_SV(ikl ,isn )
473	G1snSV(ikl,isn) = staggr * G1snSV(ikl ,isn-1)
474	. + (1. - staggr) * G1snSV(ikl ,isn )
475	G2snSV(ikl,isn) = staggr * G2snSV(ikl ,isn-1)
476	. + (1. - staggr) * G2snSV(ikl ,isn )
477	agsnSV(ikl,isn) = staggr * agsnSV(ikl ,isn-1)
478	. + (1. - staggr) * agsnSV(ikl ,isn )
479	END IF
480	END DO
481	END DO
482
483	DO ikl=1,knonv
484	isn = i_thin(ikl)
485	dzsnSV(ikl,isn) = 0.5Agrege(ikl) dzsnSV(ikl,isn)
486	. + (1.-Agrege(ikl))*dzsnSV(ikl,isn)
487
488	isn = min(i_thin(ikl) +1,nsno)
489	istoSV(ikl,isn) = Agrege(ikl) *istoSV(ikl,isn-1)
490	. + (1.-Agrege(ikl))*istoSV(ikl,isn)
491	dzsnSV(ikl,isn) = Agrege(ikl) *dzsnSV(ikl,isn-1)
492	. + (1.-Agrege(ikl))*dzsnSV(ikl,isn)
493	TsisSV(ikl,isn) = Agrege(ikl) *TsisSV(ikl,isn-1)
494	. + (1.-Agrege(ikl))*TsisSV(ikl,isn)
495	ro__SV(ikl,isn) = Agrege(ikl) *ro__SV(ikl,isn-1)
496	. + (1.-Agrege(ikl))*ro__SV(ikl,isn)
497	eta_SV(ikl,isn) = Agrege(ikl) *eta_SV(ikl,isn-1)
498	. + (1.-Agrege(ikl))*eta_SV(ikl,isn)
499	G1snSV(ikl,isn) = Agrege(ikl) *G1snSV(ikl,isn-1)
500	. + (1.-Agrege(ikl))*G1snSV(ikl,isn)
501	G2snSV(ikl,isn) = Agrege(ikl) *G2snSV(ikl,isn-1)
502	. + (1.-Agrege(ikl))*G2snSV(ikl,isn)
503	agsnSV(ikl,isn) = Agrege(ikl) *agsnSV(ikl,isn-1)
504	. + (1.-Agrege(ikl))*agsnSV(ikl,isn)
505	isnoSV(ikl) = min(Agrege(ikl) +isnoSV(ikl),real(nsno))
506	iiceSV(ikl) = iiceSV(ikl)
507	. + Agrege(ikl) *max(0,sign(1,iiceSV(ikl)
508	. -isn +icemix))
509	. *max(0,sign(1,iiceSV(ikl)
510	. -1 ))
511	END DO
512
513
514	C +--Constrains Agregation in case of too much Layers
515	C + =================================================
516
517	C +--Search the thinest non-zero Layer
518	C + -----------------------------------
519
520
521
522	DO ikl=1,knonv
523	dzthin(ikl) = 0. ! Arbitrary unrealistic
524	END DO ! Layer Thickness
525	DO ikl=1,knonv
526	DO isn=1,isnoSV(ikl)-3 ! no agregation of 3 first snowlayers
527	! XF 04/07/2019
528
529	isno_n = isnoSV(ikl)-isn+1 ! Snow Normal.Profile
530	iice_n = iiceSV(ikl)-isn ! Ice Normal.Profile
531	iiceOK = min(1,max(0,iice_n +1)) ! Ice Switch
532	! #vz dz_ref(isn) = !
533	! #vz. dz_min ((1-iiceOK)isno_n*isno_n ! Theoretical Profile
534	! #vz. + iiceOK * 2**iice_n) !
535	! #vz. /max(1,isnoSV(ikl)) !
536	dz_dif = dz_min ! Actual Profile
537	. - dzsnSV(ikl ,isn) !
538	. /max(epsi,((1-iiceOK)isno_nisno_n ! Theoretical Profile
539	. + iiceOK 2. *iice_n)) !
540	! #vz dzwdif(isn) = dz_dif !
541	OKthin = max(zero, !
542	. sign(unun, !
543	. dz_dif - dzthin(ikl)))! 1.=> New thinest Lay.
544	. * max(0, ! 1 => .le. isnoSV
545	. min(1, !
546	. isnoSV(ikl)-isn +1 )) !
547	i_thin(ikl) = (1. - OKthin) * i_thin(ikl) ! Update thinest Lay.
548	. + OKthin * isn ! Index
549	dzthin(ikl) = (1. - OKthin) * dzthin(ikl) !
550	. + OKthin * dz_dif !
551
552
553	END DO
554	END DO
555
556
557
558
559
560	C +--Index of the contiguous Layer to agregate
561	C + -----------------------------------------
562
563	C + ***************
564	call SISVAT_zCr
565	C + ***************
566
567
568	C +--Assign the 2 Layers to agregate
569	C + -------------------------------
570
571	DO ikl=1,knonv
572	isn = i_thin(ikl)
573	if(LIndsv(ikl)>0) isn=min(isn, nsno-1) !cXF
574	isagr1(ikl) = istoSV(ikl,isn)
575	isagr2(ikl) = istoSV(ikl,isn+LIndsv(ikl))
576	dzagr1(ikl) = dzsnSV(ikl,isn)
577	dzagr2(ikl) = dzsnSV(ikl,isn+LIndsv(ikl))
578	T_agr1(ikl) = TsisSV(ikl,isn)
579	T_agr2(ikl) = TsisSV(ikl,isn+LIndsv(ikl))
580	roagr1(ikl) = ro__SV(ikl,isn)
581	roagr2(ikl) = ro__SV(ikl,isn+LIndsv(ikl))
582	etagr1(ikl) = eta_SV(ikl,isn)
583	etagr2(ikl) = eta_SV(ikl,isn+LIndsv(ikl))
584	G1agr1(ikl) = G1snSV(ikl,isn)
585	G1agr2(ikl) = G1snSV(ikl,isn+LIndsv(ikl))
586	G2agr1(ikl) = G2snSV(ikl,isn)
587	G2agr2(ikl) = G2snSV(ikl,isn+LIndsv(ikl))
588	agagr1(ikl) = agsnSV(ikl,isn)
589	agagr2(ikl) = agsnSV(ikl,isn+LIndsv(ikl))
590	LstLay = min(1,max( 0, isnoSV(ikl)-1 ))
591	Agrege(ikl) = min(1,
592	. max(0,
593	. NLaysv(ikl) +isnoSV(ikl)-nsno
594	. +NLay_s(ikl) )
595	. *LstLay )
596
597	C + minimum uppermost layer thickness to guarantee a correct reproduction of the snow
598	C + atmosphere coupling
599	if(dzsnSV(ikl,max(1,isnoSV(ikl)-0))>0.02 .or. ! surface layers> 2-5-10
600	. dzsnSV(ikl,max(1,isnoSV(ikl)-1))>0.05 .or. ! XF 04/07/2019
601	. dzsnSV(ikl,max(1,isnoSV(ikl)-2))>0.10 .or.
602	. dzsnSV(ikl,max(1,isnoSV(ikl)-3))>0.30 )then
603	Agrege(ikl) = min(1,
604	. max(0,
605	. NLaysv(ikl) +isnoSV(ikl)+1-nsno ! nsno-1 layers ma
606	. +NLay_s(ikl) )
607	. *LstLay )
608	endif
609
610	isnoSV(ikl) = isnoSV(ikl)
611	. -(1-LstLay)*max(zero,
612	. sign(unun, eps_21
613	. -dzsnSV(ikl,1) ))
614	isnoSV(ikl) =max( 0, isnoSV(ikl) )
615
616	WEagre(ikl) = 0.
617	END DO
618
619	DO isn=1,nsno
620	DO ikl=1,knonv
621	WEagre(ikl) = WEagre(ikl) + ro__SV(ikl,isn)*dzsnSV(ikl,isn)
622	. *min(1,max(0,i_thin(ikl)+1-isn))
623	ENDDO
624	ENDDO
625
626	C +--Agregates
627	C + ---------
628
629	C + ***************
630	call SISVAT_zAg
631	. (isagr1,isagr2,WEagre
632	. ,dzagr1,dzagr2,T_agr1,T_agr2
633	. ,roagr1,roagr2,etagr1,etagr2
634	. ,G1agr1,G1agr2,G2agr1,G2agr2
635	. ,agagr1,agagr2,Agrege
636	. )
637	C + ***************
638
639
640	C +--Rearranges the Layers
641	C + ---------------------
642
643	C +--New (agregated) Snow layer
644	C + ^^^^^^^^^^^^^^^^^^^^^^^^^^
645	DO ikl=1,knonv
646	isn = i_thin(ikl)
647	isn = min(isn,isn+LIndsv(ikl))
648	isnoSV(ikl) = max(0.,isnoSV(ikl) -Agrege(ikl))
649	iiceSV(ikl) = iiceSV(ikl)
650	. -max(0,sign(1,iiceSV(ikl) -isn +icemix))
651	. *Agrege(ikl)
652	. *max(0,sign(1,iiceSV(ikl) -1 ))
653	istoSV(ikl,isn) = (1.-Agrege(ikl))*istoSV(ikl,isn)
654	. + Agrege(ikl) *isagr1(ikl)
655	dzsnSV(ikl,isn) = (1.-Agrege(ikl))*dzsnSV(ikl,isn)
656	. + Agrege(ikl) *dzagr1(ikl)
657	TsisSV(ikl,isn) = (1.-Agrege(ikl))*TsisSV(ikl,isn)
658	. + Agrege(ikl) *T_agr1(ikl)
659	ro__SV(ikl,isn) = (1.-Agrege(ikl))*ro__SV(ikl,isn)
660	. + Agrege(ikl) *roagr1(ikl)
661	eta_SV(ikl,isn) = (1.-Agrege(ikl))*eta_SV(ikl,isn)
662	. + Agrege(ikl) *etagr1(ikl)
663	G1snSV(ikl,isn) = (1.-Agrege(ikl))*G1snSV(ikl,isn)
664	. + Agrege(ikl) *G1agr1(ikl)
665	G2snSV(ikl,isn) = (1.-Agrege(ikl))*G2snSV(ikl,isn)
666	. + Agrege(ikl) *G2agr1(ikl)
667	agsnSV(ikl,isn) = (1.-Agrege(ikl))*agsnSV(ikl,isn)
668	. + Agrege(ikl) *agagr1(ikl)
669	END DO
670
671	C +--Above
672	C + ^^^^^
673	DO ikl=1,knonv
674	isn1(ikl)=max(i_thin(ikl),i_thin(ikl)+LIndsv(ikl))
675	END DO
676	DO i= 1,nsno-1
677	DO ikl=1,knonv
678	staggr = min(1,max(0,i +1 -isn1(ikl) ))
679	istoSV(ikl,i) = (1.-staggr )*istoSV(ikl,i )
680	. + staggr((1.-Agrege(ikl))istoSV(ikl,i )
681	. + Agrege(ikl) *istoSV(ikl,i+1))
682	dzsnSV(ikl,i) = (1.-staggr )*dzsnSV(ikl,i )
683	. + staggr((1.-Agrege(ikl))dzsnSV(ikl,i )
684	. + Agrege(ikl) *dzsnSV(ikl,i+1))
685	TsisSV(ikl,i) = (1.-staggr )*TsisSV(ikl,i )
686	. + staggr((1.-Agrege(ikl))TsisSV(ikl,i )
687	. + Agrege(ikl) *TsisSV(ikl,i+1))
688	ro__SV(ikl,i) = (1.-staggr )*ro__SV(ikl,i )
689	. + staggr((1.-Agrege(ikl))ro__SV(ikl,i )
690	. + Agrege(ikl) *ro__SV(ikl,i+1))
691	eta_SV(ikl,i) = (1.-staggr )*eta_SV(ikl,i )
692	. + staggr((1.-Agrege(ikl))eta_SV(ikl,i )
693	. + Agrege(ikl) *eta_SV(ikl,i+1))
694	G1snSV(ikl,i) = (1.-staggr )*G1snSV(ikl,i )
695	. + staggr((1.-Agrege(ikl))G1snSV(ikl,i )
696	. + Agrege(ikl) *G1snSV(ikl,i+1))
697	G2snSV(ikl,i) = (1.-staggr )*G2snSV(ikl,i )
698	. + staggr((1.-Agrege(ikl))G2snSV(ikl,i )
699	. + Agrege(ikl) *G2snSV(ikl,i+1))
700	agsnSV(ikl,i) = (1.-staggr )*agsnSV(ikl,i )
701	. + staggr((1.-Agrege(ikl))agsnSV(ikl,i )
702	. + Agrege(ikl) *agsnSV(ikl,i+1))
703	END DO
704	END DO
705
706	DO ikl=1,knonv
707	isn = min(isnoSV(ikl) +1,nsno)
708	istoSV(ikl,isn) = (1.-Agrege(ikl))*istoSV(ikl,isn)
709	dzsnSV(ikl,isn) = (1.-Agrege(ikl))*dzsnSV(ikl,isn)
710	TsisSV(ikl,isn) = (1.-Agrege(ikl))*TsisSV(ikl,isn)
711	ro__SV(ikl,isn) = (1.-Agrege(ikl))*ro__SV(ikl,isn)
712	eta_SV(ikl,isn) = (1.-Agrege(ikl))*eta_SV(ikl,isn)
713	G1snSV(ikl,isn) = (1.-Agrege(ikl))*G1snSV(ikl,isn)
714	G2snSV(ikl,isn) = (1.-Agrege(ikl))*G2snSV(ikl,isn)
715	agsnSV(ikl,isn) = (1.-Agrege(ikl))*agsnSV(ikl,isn)
716	END DO
717
718
719	return
720	end

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

Download in other formats:

Original Format