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)>0..AND.i_thin(ikl)<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.5*Agrege(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_n*isno_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 |
---|