Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

sisvat_zsn.F

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

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

Download in other formats:

Original Format