source: LMDZ6/branches/blowing_snow/libf/phylmd/inlandsis/sisvat_zsn.F @ 5327

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