source: LMDZ6/trunk/libf/phylmd/inlandsis/sisvat_zsn.f90 @ 5452

Last change on this file since 5452 was 5246, checked in by abarral, 2 months 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
3subroutine 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.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  ! +--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_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  ! +--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
741end subroutine sisvat_zsn
Note: See TracBrowser for help on using the repository browser.