Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

sisvat_zsn.f90 @ 5282

Last change on this file since 5282 was 5246, checked in by abarral, 10 days ago
Convert fixed-form to free-form sources .F -> .{f,F}90 (WIP: some .F remain, will be handled in subsequent commits)
File size: 28.9 KB

Line
1
2
3	subroutine SISVAT_zSn
4
5	! +------------------------------------------------------------------------+
6	! \| MAR SISVAT_zSn 12-07-2019 MAR \|
7	! \| SubRoutine SISVAT_zSn manages the Snow Pack vertical Discretization \|
8	! \| \|
9	! +------------------------------------------------------------------------+
10	! \| \|
11	! \| PARAMETERS: knonv: Total Number of columns = \|
12	! \| ^^^^^^^^^^ = Total Number of continental grid boxes \|
13	! \| X Number of Mosaic Cell per grid box \|
14	! \| \|
15	! \| INPUT / NLaysv = New Snow Layer Switch \|
16	! \| OUTPUT: isnoSV = total Nb of Ice/Snow Layers \|
17	! \| ^^^^^^ ispiSV = 0,...,nsno: Uppermost Superimposed Ice Layer \|
18	! \| iiceSV = total Nb of Ice Layers \|
19	! \| istoSV = 0,...,5 : Snow History (see istdSV data) \|
20	! \| \|
21	! \| INPUT / TsisSV : Soil/Ice Temperatures (layers -nsol,-nsol+1,..,0)\|
22	! \| OUTPUT: & Snow Temperatures (layers 1,2,...,nsno) [K] \|
23	! \| ^^^^^^ ro__SV : Soil/Snow Volumic Mass [kg/m3] \|
24	! \| eta_SV : Soil/Snow Water Content [m3/m3] \|
25	! \| dzsnSV : Snow Layer Thickness [m] \|
26	! \| G1snSV : Dendricity (<0) or Sphericity (>0) of Snow Layer \|
27	! \| G2snSV : Sphericity (>0) or Size of Snow Layer \|
28	! \| agsnSV : Snow Age [day] \|
29	! \| \|
30	! \| METHOD: 1) Agregate the thinest Snow Layer \|
31	! \| ^^^^^^ if a new Snow Layer has been precipitated (NLaysv = 1) \|
32	! \| 2) Divide a too thick Snow Layer except \|
33	! \| if the maximum Number of Layer is reached \|
34	! \| in this case forces NLay_s = 1 \|
35	! \| 3) Agregate the thinest Snow Layer \|
36	! \| in order to divide a too thick Snow Layer \|
37	! \| at next Time Step when NLay_s = 1 \|
38	! \| \|
39	! +------------------------------------------------------------------------+
40
41
42
43
44	! +--Global Variables
45	! + ================
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	! +--Internal Variables
60	! + ==================
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	! + ! 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	! +--DATA
110	! + ====
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	! #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	! + CAUTION: dz_max > dz_min*2 is required ! Otherwise re-agregation is
119	! + ! activated after splitting
120
121
122
123
124
125	! +--Constrains Agregation of too thin Layers
126	! + =================================================
127
128	! +--Search the thinest non-zero Layer
129	! + ----------------------------------
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	!XF
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	! + !
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	! + !
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	! + ***************
233	call SISVAT_zCr
234	! + ***************
235
236
237	! +--Assign the 2 Layers to agregate
238	! + -------------------------------
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	! +--Agregates
286	! + ---------
287
288	! + ***************
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	! + ***************
297
298
299	! +--Rearranges the Layers
300	! + ---------------------
301
302	! +--New (agregated) Snow layer
303	! + ^^^^^^^^^^^^^^^^^^^^^^^^^^
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	! +--Above
331	! + ^^^^^
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	! +--Constrains Splitting of too thick Layers
381	! + =================================================
382
383
384	! +--Search the thickest non-zero Layer
385	! + ----------------------------------
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	! +--Rearranges the Layers
457	! + ---------------------
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	! +--Constrains Agregation in case of too much Layers
516	! + =================================================
517
518	! +--Search the thinest non-zero Layer
519	! + -----------------------------------
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	! +--Index of the contiguous Layer to agregate
562	! + -----------------------------------------
563
564	! + ***************
565	call SISVAT_zCr
566	! + ***************
567
568
569	! +--Assign the 2 Layers to agregate
570	! + -------------------------------
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	! + minimum uppermost layer thickness to guarantee a correct reproduction of the snow
599	! + 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	! +--Agregates
628	! + ---------
629
630	! + ***************
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	! + ***************
639
640
641	! +--Rearranges the Layers
642	! + ---------------------
643
644	! +--New (agregated) Snow layer
645	! + ^^^^^^^^^^^^^^^^^^^^^^^^^^
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	! +--Above
673	! + ^^^^^
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	! +--Search new Ice/Snow Interface (option II in MAR)
721	! + ===============================================
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 subroutine sisvat_zsn

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

Download in other formats:

Original Format