source: LMDZ6/trunk/libf/phylmd/inlandsis/sisvat_qsn.F @ 4743

Last change on this file since 4743 was 3900, checked in by evignon, 3 years ago

Commission de la nouvelle interface LMDZ-SISVAT
la prochaine commission consistera a supprimer l'ancien repertoire sisvat
et a faire un peu de nettoyage.

File size: 27.8 KB
Line 
1 
2 
3      subroutine SISVAT_qSn
4     .                     (
5! #e1.                      EqSn_0,EqSn_1,EqSn_d
6! #m1.                     ,SIsubl,SImelt,SIrnof
7     .                     )
8 
9C +------------------------------------------------------------------------+
10C | MAR          SISVAT_qSn                           Fri 29-Jul-2011  MAR |
11C |   SubRoutine SISVAT_qSn updates  the Snow Water Content                |
12C +------------------------------------------------------------------------+
13C |                                                                        |
14C |   PARAMETERS:  knonv: Total Number of columns =                        |
15C |   ^^^^^^^^^^        = Total Number of continental     grid boxes       |
16C |                     X       Number of Mosaic Cell per grid box         |
17C |                                                                        |
18C |   INPUT:   isnoSV   = total Nb of Ice/Snow Layers                      |
19C |   ^^^^^                                                                |
20C |                                                                        |
21C |   INPUT:   TaT_SV   : SBL Top    Temperature                       [K] |
22C |   ^^^^^    dt__SV   : Time Step                                    [s] |
23C |                                                                        |
24C |   INPUT /  drr_SV   : Rain Intensity                         [kg/m2/s] |
25C |   OUTPUT:  dzsnSV   : Snow Layer Thickness                         [m] |
26C |   ^^^^^^   eta_SV   : Snow Water Content                       [m3/m3] |
27C |            ro__SV   : Snow/Soil  Volumic Mass                  [kg/m3] |
28C |            TsisSV   : Soil/Ice Temperatures (layers -nsol,-nsol+1,..,0)|
29C |                     & Snow     Temperatures (layers  1,2,...,nsno) [K] |
30C |                                                                        |
31C |   OUTPUT:  SWS_SV   : Surficial Water Status                           |
32C |   ^^^^^^                                                               |
33C |            EExcsv   : Snow Energy in Excess, initial Forcing    [J/m2] |
34C |            EqSn_d   : Snow Energy in Excess, remaining          [J/m2] |
35C |            EqSn_0   : Snow Energy, before Phase Change          [J/m2] |
36C |            EqSn_1   : Snow Energy, after  Phase Change          [J/m2] |
37C |            SIsubl   : Snow sublimed/deposed Mass             [mm w.e.] |
38C |            SImelt   : Snow Melted           Mass             [mm w.e.] |
39C |            SIrnof   : Surficial Water + Run OFF Change       [mm w.e.] |
40C |                                                                        |
41C |   Internal Variables:                                                  |
42C |   ^^^^^^^^^^^^^^^^^^                                                   |
43C |                                                                        |
44C | # OPTIONS: #E0: IO for Verification: Energy       Budget               |
45C | # ^^^^^^^                                                              |
46C | #          #su: IO for Verification: Slush        Diagnostic           |
47C |                                                                        |
48C |                                                                        |
49C +------------------------------------------------------------------------+
50 
51 
52 
53 
54C +--Global Variables
55C +  ================
56
57      use VARphy
58      use VAR_SV
59      use VARdSV
60      use VAR0SV
61      use VARxSV
62      use VARySV
63      use surface_data, only: is_ok_slush,opt_runoff_ac
64
65
66      IMPLICIT NONE
67
68 
69! Energy          Budget
70! ~~~~~~~~~~~~~~~~~~~~~~
71! #e1 real     EqSn_d(knonv)                 ! Energy in Excess, initial
72! #e1 real     EqSn_0(knonv)                 ! Snow Energy, befor Phase Change
73! #vm real     EqSn01(knonv)                 ! Snow Energy, after Phase Change
74! #vm real     EqSn02(knonv)                 ! Snow Energy, after Phase Change
75                                             !              .AND. Last Melting
76! #e1 real     EqSn_1(knonv)                 ! Snow Energy, after Phase Change
77                                             !              .AND. Mass Redistr.
78! Snow/Ice (Mass) Budget
79! ~~~~~~~~~~~~~~~~~~~~~~
80! #m1 real     SIsubl(knonv)                 ! Snow Deposed Mass
81! #m1 real     SImelt(knonv)                 ! Snow Melted  Mass
82! #m1 real     SIrnof(knonv)                 ! Local Surficial Water + Run OFF
83 
84 
85C +--Internal Variables
86C +  ==================
87 
88      integer  ikl   ,isn                    !
89      integer  nh                            ! Non erodible Snow: up.lay.Index
90      integer  LayrOK                        ! 1 (0)  if In(Above) Snow Pack
91      integer  k_face                        ! 1 (0)  if Crystal(no) faceted
92      integer  LastOK                        ! 1 ==>  1! Snow Layer
93      integer  NOLayr                        ! 1     Layer  Update
94      integer  noSnow(knonv)                 ! Nb of Layers Updater
95      integer  kSlush                        ! Slush Switch
96      real     dTSnow                        ! Temperature                  [C]
97      real     EExdum(knonv)                 ! Energy in Excess when no Snow
98      real     OKmelt                        ! 1 (0)  if        (no) Melting
99      real     EnMelt                        ! Energy in excess, for Melting
100      real     SnHLat                        ! Energy consumed   in  Melting
101      real     AdEnrg,B_Enrg                 ! Additional Energy from  Vapor
102      real     dzVap0,dzVap1                 ! Vaporized Thickness          [m]
103      real     dzMelt(knonv)                 ! Melted    Thickness          [m]
104      real     rosDry                        ! Snow volumic Mass if no Water in
105      real     PorVol                        ! Pore volume
106      real     PClose                        ! Pore Hole Close OFF Switch
107      real     SGDiam                        !      Snow Grain Diameter
108      real     SGDmax                        ! Max. Snow Grain Diameter
109      real     rWater                        ! Retained Water           [kg/m2]
110      real     drrNEW                        ! New available Water      [kg/m2]
111      real     rdzNEW                        ! Snow          Mass       [kg/m2]
112      real     rdzsno                        ! Snow          Mass       [kg/m2]
113      real     EnFrez                        ! Energy Release    in  Freezing
114      real     WaFrez                        ! Water  consumed   in  Melting
115      real     RapdOK                        ! 1. ==> Snow melts rapidly
116      real     ThinOK                        ! 1. ==> Snow Layer is thin
117      real     dzepsi                        ! Minim. Snow Layer Thickness (!)
118      real     dz_Min                        ! Minim. Snow Layer Thickness
119      real     z_Melt                        ! Last (thin) Layer Melting
120      real     rusnew                        ! Surficial Water Thickness   [mm]
121      real     zWater                        ! Max Slush Water Thickness   [mm]
122      real     zSlush                        !     Slush Water Thickness   [mm]
123      real     ro_new                        ! New Snow/ice  Density    [kg/m3]
124      real     zc,zt                         ! Non erod.Snow Thickness[mm w.e.]
125      real     rusnSV0(knonv)
126      real     Tsave
127 
128C +--OUTPUT of SISVAT Trace Statistics (see assignation in PHY_SISVAT)
129C +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
130      integer             isnnew,isinew,isnUpD,isnitr
131 
132! OUTPUT in SISVAT at specified i,j,k,n (see assignation in PHY_SISVAT)
133! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
134! #wx integer             iSV_v1,jSV_v1,nSV_v1,kSV_v1,lSV_v1
135! #wx common/SISVAT_EV/   iSV_v1,jSV_v1,nSV_v1,kSV_v1,lSV_v1
136 
137C +--Energy and Mass Budget
138C +  ~~~~~~~~~~~~~~~~~~~~~~
139! #vm real     WqSn_0(knonv)                 ! Snow Water+Forcing  Initial
140! #vm real     WqSn_1(knonv)                 ! Snow Water+Forcing, Final
141! #vm logical         emopen                 ! IO   Switch
142! #vm common/Se_qSn_L/emopen                 !
143! #vm integer         no_err                 !
144! #vm common/Se_qSn_I/no_err                 !
145! #vm real     hourer,timeer                 !
146! #vm common/Se_qSn_R/timeer                 !
147 
148C +--Slush Diagnostic: IO
149C +  ~~~~~~~~~~~~~~~~~~~~
150! #vu logical         su_opn                 ! IO   Switch
151! #vu common/SI_qSn_L/su_opn                 !
152 
153 
154C +--DATA
155C +  ====
156 
157      data      dzepsi/0.0001/                ! Minim. Snow Layer Thickness (!)
158c #?? data      dz_Min/0.005/                 ! Minim. Snow Layer Thickness
159c ... Warning: Too high for Col de Porte: precludes 1st snow (layer) apparition
160      data      dz_Min/2.5e-3/                ! Minim. Snow Layer Thickness
161      data      SGDmax/0.003/                 ! Maxim. Snow Grain Diameter  [m]
162                                              ! (Rowe et al. 1995, JGR p.16268)
163 
164C +--Energy Budget (IN)
165C +  ==================
166 
167! #e1   DO ikl=1,knonv
168! #e1     EqSn_0(ikl) = 0.
169! #e1   END DO
170! #e1 DO   isn=nsno,1,-1
171! #e1   DO ikl=1,knonv
172! #e1     EqSn_0(ikl) = EqSn_0(ikl) + ro__SV(ikl,isn) *dzsnSV(ikl,isn)
173! #e1.                *(Cn_dSV      *(TsisSV(ikl,isn) -TfSnow         )
174! #e1.                 -Lf_H2O      *(1.              -eta_SV(ikl,isn)))
175! #e1   END DO
176! #e1 END DO
177 
178 
179C +--Water  Budget (IN)
180C +  ==================
181 
182! #vm   DO ikl=1,knonv
183! #vm     WqSn_0(ikl) = drr_SV(ikl) * dt__SV
184! #vm.                 +rusnSV(ikl)
185! #vm   END DO
186! #vm DO   isn=nsno,1,-1
187! #vm   DO ikl=1,knonv
188! #vm     WqSn_0(ikl) = WqSn_0(ikl) + ro__SV(ikl,isn) *dzsnSV(ikl,isn)
189! #vm   END DO
190! #vm END DO
191 
192 
193C +--Snow Melt Budget
194C +  ================
195 
196! #m1   DO ikl=1,knonv
197! #m1     SImelt(ikl) = 0.
198! #m1     SIrnof(ikl) = rusnSV(ikl) + RnofSV(ikl) * dt__SV
199! #m1   END DO
200 
201 
202C +--Initialization
203C +  ==============
204 
205      DO ikl=1,knonv
206        noSnow(ikl)   = 0                   ! Nb of Layers Updater
207        ispiSV(ikl)   = 0                   ! Pore Hole Close OFF Index
208                                            ! (assumed to be the Top of
209                                            !  the surimposed Ice Layer)
210        zn5_SV(ikl)   = 0.
211        rusnSV0(ikl)  = 0.
212
213      END DO
214 
215 
216C +--Melting/Freezing Energy
217C +  =======================
218 
219C +...REMARK: Snow liquid Water Temperature assumed = TfSnow
220C +   ^^^^^^
221        DO ikl=1,knonv
222          EExdum(ikl) = drr_SV(ikl)     * C__Wat *(TaT_SV(ikl)-TfSnow)
223     .                                  * dt__SV
224          EExcsv(ikl) = EExdum(ikl)     *    min(1,isnoSV(ikl)) ! Snow exists
225          EExdum(ikl) = EExdum(ikl)     -          EExcsv(ikl)  !
226! #e1     EqSn_d(ikl) = EExcsv(ikl)                             !
227        END DO
228 
229 
230C +--Surficial Water Status
231C +  ----------------------
232 
233        DO ikl=1,knonv
234          SWS_SV(ikl) = max(zero,sign(unun,TfSnow
235     .                                    -TsisSV(ikl,isnoSV(ikl))))
236        END DO
237 
238      DO ikl=1,knonv
239
240      DO isn=min(nsno,isnoSV(ikl)+1),1,-1
241! EV          DO isn=nsno,1,-1
242C +--Energy, store Previous Content
243C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
244          dTSnow      = TsisSV(ikl,isn) -          TfSnow
245          EExcsv(ikl) = EExcsv(ikl)
246     .                + ro__SV(ikl,isn) * Cn_dSV * dTSnow
247     .                                           * dzsnSV(ikl,isn)
248          TsisSV(ikl,isn) =                        TfSnow
249 
250C +--Water,  store Previous Content
251C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
252          drr_SV(ikl) = drr_SV(ikl)
253     .                + ro__SV(ikl,isn)          * eta_SV(ikl,isn)
254     .                                           * dzsnSV(ikl,isn)
255     .                / dt__SV
256          ro__SV(ikl,isn) =
257     .                  ro__SV(ikl,isn) *(1.     - eta_SV(ikl,isn))
258          eta_SV(ikl,isn) =  0.
259 
260 
261C +--Melting  if EExcsv > 0
262C +  ======================
263 
264          EnMelt      =    max(zero,          EExcsv(ikl) )
265 
266C +--Energy Consumption
267C +  ^^^^^^^^^^^^^^^^^^
268          SnHLat      = ro__SV(ikl,isn) * Lf_H2O
269          dzMelt(ikl) = EnMelt      / max(SnHLat,    epsi )
270          noSnow(ikl) = noSnow(ikl)
271     .      + max(zero  ,sign(unun,dzMelt(ikl)                !
272     .                            -dzsnSV(ikl ,isn)))         ! 1 if full Melt
273     .       *min(1     , max(0 ,1+isnoSV(ikl)-isn))          ! 1 in the  Pack
274          dzMelt(ikl) =
275     .              min(dzsnSV(ikl, isn),dzMelt(ikl))
276          dzsnSV(ikl,isn) =
277     .                  dzsnSV(ikl,isn) -dzMelt(ikl)
278          zn5_SV(ikl) = zn5_SV(ikl)     +dzMelt(ikl)
279          EExcsv(ikl) = EExcsv(ikl)     -dzMelt(ikl)*SnHLat
280          wem_SV(ikl) = wem_SV(ikl)     -dzMelt(ikl)*ro__SV(ikl,isn)
281 
282C +--Water  Production
283C +  ^^^^^^^^^^^^^^^^^
284          drr_SV(ikl) = drr_SV(ikl)
285     .                + ro__SV(ikl,isn) * dzMelt(ikl)/dt__SV
286! #m1     SImelt(ikl) = SImelt(ikl)
287! #m1.                + ro__SV(ikl,isn) * dzMelt(ikl)
288          OKmelt      =max(zero,sign(unun,drr_SV(ikl)-epsi))
289 
290C +--Snow History
291C +  ^^^^^^^^^^^^
292          k_face          =       min(    istoSV(ikl,isn),istdSV(1)) ! = 1  if
293     .                           *max(0,2-istoSV(ikl,isn)          ) ! faceted
294          istoSV(ikl,isn) =                                          !
295     .        (1.-OKmelt) *               istoSV(ikl,isn)            !
296     .      +     OKmelt  *((1-k_face) *  istdSV(2)                  !
297     .                     +   k_face  *  istdSV(3)      )           !
298 
299 
300C +--Freezing if EExcsv < 0
301C +  ======================
302 
303          rdzsno      =          ro__SV(ikl,isn) * dzsnSV(ikl ,isn)
304          LayrOK      = min(   1, max(0          , isnoSV(ikl)-isn+1))
305          EnFrez      = min(zero,                  EExcsv(ikl))
306          WaFrez      =   -(     EnFrez          * LayrOK / Lf_H2O)
307          drrNEW      = max(zero,drr_SV(ikl)     - WaFrez / dt__SV)
308          WaFrez      =    (     drr_SV(ikl)     - drrNEW)* dt__SV
309          drr_SV(ikl) =          drrNEW
310          EExcsv(ikl) =          EExcsv(ikl)     + WaFrez * Lf_H2O
311          EnFrez      = min(zero,EExcsv(ikl))    * LayrOK
312          rdzNEW      = WaFrez + rdzsno
313          ro__SV(ikl,isn) =      rdzNEW /max(epsi, dzsnSV(ikl,isn))
314          TsisSV(ikl,isn) =      TfSnow
315     .                + EnFrez /(Cn_dSV *max(epsi, rdzNEW)        )
316          EExcsv(ikl) =          EExcsv(ikl)     - EnFrez
317          wer_SV(ikl) = WaFrez
318     .                + wer_SV(ikl)
319 
320 
321
322C +--Snow Water Content
323C +  ==================
324 
325C +--Percolation Velocity
326C +  ^^^^^^^^^^^^^^^^^^^^
327c #PW     SGDiam    = 1.6d-4
328c #PW.              + 1.1d-13 *(ro__SV(ikl,isn)*ro__SV(ikl,isn)
329c #PW.                         *ro__SV(ikl,isn)*ro__SV(ikl,isn))
330 
331C +--Pore   Volume [-]
332C +  ^^^^^^^^^^^^^^^^^
333          rosDry      =(1.     - eta_SV(ikl,isn))* ro__SV(ikl,isn) !
334          PorVol      = 1.     - rosDry          / ro_Ice          !
335          PorVol      =      max(PorVol          , zero  )         !
336 
337C +--Water  Retention
338C +  ^^^^^^^^^^^^^^^^
339          rWater      = ws0dSV * PorVol * ro_Wat * dzsnSV(ikl,isn)
340          drrNEW      = max(zero,drr_SV(ikl)     - rWater /dt__SV)
341          rWater      =    (     drr_SV(ikl)     - drrNEW)*dt__SV
342          drr_SV(ikl)     =      drrNEW
343          rdzNEW          =      rWater
344     .                         + rosDry          * dzsnSV(ikl,isn)
345          eta_SV(ikl,isn) =      rWater / max(epsi,rdzNEW)
346          ro__SV(ikl,isn) =      rdzNEW / max(epsi,dzsnSV(ikl,isn))
347 
348C +--Pore Hole Close OFF
349C +  ^^^^^^^^^^^^^^^^^^^
350          PClose = max(zero,
351     .                 sign(unun,ro__SV(ikl,isn)
352     .                          -roCdSV         ))
353          ispiSV(ikl) =          ispiSV(ikl)      *(1.-PClose)
354     .                +      max(ispiSV(ikl),isn)    * Pclose
355          PClose = max(0   ,                         ! Water under SuPer.Ice
356     .                 min (1   ,ispiSV(ikl)         ! contributes to
357     .                          -isn            ))   ! Surficial   Water
358 
359cXF
360          if(ro__SV(ikl,isn) >= roCdSV.and.ro__SV(ikl,1)<900)
361     .    PClose = min(0.50,PClose *
362     .    (1.-(ro_ice-ro__SV(ikl,isn))/(ro_ice-roCdSV)))
363 
364          PClose = max(0.,min(1.,PClose))
365 
366          if(isn==1) then
367               PClose = 1
368           ispiSV(ikl)= max(ispiSV(ikl),1)
369          endif
370 
371          if(drr_SV(ikl)    >0  .and.TsisSV(ikl,isn)>273.14) then
372           if((ro__SV(ikl,isn)>900.and.ro__SV(ikl,isn)<920).or.
373     .         ro__SV(ikl,isn)>950) then
374             dzsnSV(ikl,isn) = dzsnSV(ikl,isn)*ro__SV(ikl,isn)/ro_ice
375             ro__SV(ikl,isn) = ro_ice
376             PClose          = 1
377           endif
378          endif
379 
380c         if (isn>1.and.isn<nsno     .and.
381c    .      ro__SV(ikl,isn-1)>900    .and.
382c    .      ro__SV(ikl,isn)  >roCdSV .and.
383c    .      ro__SV(ikl,isn)  <900    .and.
384c    .      TsisSV(ikl,isn)  >273.14 .and.
385c    .      TsisSV(ikl,isn+1)<273.15 .and.
386c    .      drr_SV(ikl)      >0)     then
387c          TsisSV(ikl,isn)=273.14
388c          PClose = 1
389c         endif
390 
391cXF
392          rusnSV(ikl) =          rusnSV(ikl)
393     .                +          drr_SV(ikl) *dt__SV * PClose
394          rusnSV0(ikl)=         rusnSV0(ikl)
395     .                +          drr_SV(ikl) *dt__SV * PClose
396          drr_SV(ikl) =          drr_SV(ikl)      *(1.-PClose)
397 
398        END DO
399
400      END DO
401 
402 
403C +--Remove Zero-Thickness Layers
404C +  ============================
405 
406 1000 CONTINUE
407           isnitr =          0
408      DO   ikl=1,knonv
409           isnUpD =          0
410           isinew =          0
411cXF
412
413
414        DO isn=1,min(nsno-1,isnoSV(ikl))
415           isnnew =(unun-max(zero  ,sign(unun,dzsnSV(ikl,isn)-dzepsi)))
416     .             *     max(0     , min(1   ,isnoSV(ikl) +1 -isn ))
417           isnUpD =      max(isnUpD,          isnnew)
418           isnitr =      max(isnitr,          isnnew)
419           isinew =      isn*isnUpD *max(0, 1-isinew)      ! LowerMost  0-Layer
420     .                                       +isinew       ! Index
421           dzsnSV(ikl,isn) =                  dzsnSV(ikl,isn+isnnew)
422           ro__SV(ikl,isn) =                  ro__SV(ikl,isn+isnnew)
423           TsisSV(ikl,isn) =                  TsisSV(ikl,isn+isnnew)
424           eta_SV(ikl,isn) =                  eta_SV(ikl,isn+isnnew)
425           G1snSV(ikl,isn) =                  G1snSV(ikl,isn+isnnew)
426           G2snSV(ikl,isn) =                  G2snSV(ikl,isn+isnnew)
427           dzsnSV(ikl,isn+isnnew) =(1-isnnew)*dzsnSV(ikl,isn+isnnew)
428           ro__SV(ikl,isn+isnnew) =(1-isnnew)*ro__SV(ikl,isn+isnnew)
429           eta_SV(ikl,isn+isnnew) =(1-isnnew)*eta_SV(ikl,isn+isnnew)
430           G1snSV(ikl,isn+isnnew) =(1-isnnew)*G1snSV(ikl,isn+isnnew)
431           G2snSV(ikl,isn+isnnew) =(1-isnnew)*G2snSV(ikl,isn+isnnew)
432
433        END DO
434           isnoSV(ikl)   =   isnoSV(ikl)-isnUpD            ! Nb of Snow   Layer
435           ispiSV(ikl)   =   ispiSV(ikl)                   ! Nb of SuperI Layer
436     .    -isnUpD *max(0,min(ispiSV(ikl)-isinew,1))        ! Update  if I=0
437 
438      END DO
439
440      IF  (isnitr.GT.0)                                       GO TO 1000
441 
442 
443C +--New upper Limit of the non erodible Snow (istoSV .GT. 1)
444C +  ========================================
445 
446      DO   ikl=1,knonv
447           nh =     0
448cXF
449        DO isn=  isnoSV(ikl),1,-1
450           nh =    nh + isn* min(istoSV(ikl,isn)-1,1)*max(0,1-nh)
451        ENDDO
452           zc =     0.
453           zt =     0.
454cXF
455        DO isn=1,isnoSV(ikl)
456           zc =    zc +          dzsnSV(ikl,isn) *ro__SV(ikl,isn)
457     .                     * max(0,min(1,nh+1-isn))
458           zt =    zt +          dzsnSV(ikl,isn) *ro__SV(ikl,isn)
459        END DO
460           zWE_SV(ikl) =                 zt
461           zWEcSV(ikl) = min(zWEcSV(ikl),zt)
462           zWEcSV(ikl) = max(zWEcSV(ikl),zc)
463      END DO
464 
465 
466C +--Energy Budget (OUT)
467C +  ===================
468 
469! #vm   DO ikl=1,knonv
470! #vm     EqSn01(ikl) =-EqSn_0(ikl)
471! #vm.                 -EExcsv(ikl)
472! #vm   END DO
473! #vm DO   isn=nsno,1,-1
474! #vm   DO ikl=1,knonv
475! #vm     EqSn01(ikl) = EqSn01(ikl) + ro__SV(ikl,isn) *dzsnSV(ikl,isn)
476! #vm.                *(Cn_dSV      *(TsisSV(ikl,isn) -TfSnow         )
477! #vm.                 -Lf_H2O      *(1.              -eta_SV(ikl,isn)))
478! #vm   END DO
479! #vm END DO
480 
481 
482C +--"Negative Heat" from supercooled rain
483C +   ------------------------------------
484 
485      DO ikl=1,knonv
486          EExcsv(ikl) = EExcsv(ikl) + EExdum(ikl)
487 
488 
489C +--Surficial Water Run OFF
490C +  -----------------------
491 
492          rusnew      = rusnSV(ikl) * SWf_SV(ikl)
493 
494          if(isnoSV(ikl)<=1 .OR. opt_runoff_ac) rusnew = 0.
495          !if(ivgtSV(ikl)>=1) rusnew = 0.
496 
497c #EU                        rusnew = 0.
498c #AC               rusnew = 0.
499
500          RnofSV(ikl) = RnofSV(ikl)
501     .                +(rusnSV(ikl) - rusnew     ) / dt__SV
502          RuofSV(ikl,1) = RuofSV(ikl,1)
503     .                +(rusnSV(ikl) - rusnew     ) / dt__SV
504          RuofSV(ikl,4) = RuofSV(ikl,4)
505     .                +(rusnSV0(ikl)             ) / dt__SV
506          rusnSV(ikl) = rusnew
507      END DO
508 
509 
510C +--Percolation down the Continental Ice Pack
511C +  -----------------------------------------
512 
513        DO ikl=1,knonv
514          drr_SV(ikl) = drr_SV(ikl) + rusnSV(ikl)
515     .                     * (1-min(1,ispiSV(ikl)))/ dt__SV
516          rusnSV(ikl) = rusnSV(ikl)
517     .                     *    min(1,ispiSV(ikl))
518        END DO
519 
520cXF removal of too thin snowlayers if TT> 275.15 + bug if TT>> 273.15
521        DO ikl=1,knonv
522         zt=0.
523         DO isn=1,isnoSV(ikl)
524          zt=zt+dzsnSV(ikl,isn)
525         ENDDO
526         
527         if(zt<0.005+(TaT_SV(ikl)-TfSnow)/1000..and.
528     .      isnoSV(ikl)             >0         .and.
529     .      TaT_SV(ikl)             >=TfSnow   .and.
530     .      istoSV(ikl,isnoSV(ikl)) >1       ) then
531          DO isn=1,isnoSV(ikl)
532           drr_SV(ikl)    = drr_SV(ikl)
533     .                    + dzsnSV(ikl,isn)*ro__SV(ikl,isn) /dt__SV
534           dzsnSV(ikl,isn)= 0.
535
536          ENDDO
537          isnoSV(ikl)     = 0
538         endif
539        ENDDO
540 
541C +--Slush Formation (Activated. CAUTION: ADD RunOff Possibility before Activation)
542C +  ---------------  ^^^^^^^  ^^^
543 
544      IF (is_ok_slush) THEN
545
546      DO  ikl=1,knonv
547       DO isn=1,isnoSV(ikl)
548          kSlush = min(1,max(0,isn+1-ispiSV(ikl)))        ! Slush Switch
549 
550C +--Available Additional Pore   Volume [-]
551C +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
552          PorVol = 1. - ro__SV(ikl,isn)                    ! [--]
553     .           *(1. - eta_SV(ikl,isn))/ ro_Ice           !
554     .           -      eta_SV(ikl,isn)                    !
555     .                 *ro__SV(ikl,isn) / ro_Wat           !
556          PorVol =  max(PorVol          , zero  )          !
557          zWater =      dzsnSV(ikl,isn) * PorVol * 1000.   ! [mm] OR [kg/m2]
558     .           * (1. -SWS_SV(ikl)                        ! 0 <=> freezing
559     .                *(1 -min(1,iabs(isn-isnoSV(ikl)))))  ! 1 <=> isn=isnoSV
560          zSlush =  min(rusnSV(ikl)     , zWater)          ! [mm] OR [kg/m2]
561          ro_new      =(dzsnSV(ikl,isn) * ro__SV(ikl,isn)  !
562     .                 +zSlush                           ) !
563     .            / max(dzsnSV(ikl,isn) , epsi           ) !
564          if(ro_new<ro_Ice+20) then ! MAX 940kg/m3         !
565           rusnSV(ikl)  = rusnSV(ikl)          - zSlush    ! [mm] OR [kg/m2]
566           RuofSV(ikl,4)= max(0.,RuofSV(ikl,4) - zSlush/dt__SV)
567           eta_SV(ikl,isn) =(ro_new - ro__SV(ikl,isn)      !
568     .                     *(1.     - eta_SV(ikl,isn)))    !
569     .                / max (ro_new , epsi            )    !
570           ro__SV(ikl,isn) =      ro_new                   !
571          endif
572        END DO
573      END DO
574      END IF
575 
576C +--Impact of the Sublimation/Deposition on the Surface Mass Balance
577C +  ================================================================
578 
579      DO ikl=1,knonv
580         isn                     = isnoSV(ikl)
581          dzVap0                  =                   dt__SV
582     .  * HLs_sv(ikl)         * min(isn             , 1   )
583     .  /(Lx_H2O(ikl)         * max(ro__SV(ikl,isn) , epsi))
584          NOLayr=min(zero,sign(unun,dzsnSV(ikl,isn) + dzVap0))
585          dzVap1=min(zero,          dzsnSV(ikl,isn) + dzVap0)
586 
587 
588C +--Additional Energy
589C +  -----------------
590 
591c #VH     AdEnrg = dzVap0 * ro__SV(ikl,isnoSV(ikl))           ! Water   Vapor
592c #VH.            *C__Wat *(TsisSV(ikl,isnoSV(ikl)) -TfSnow)  ! Sensible Heat
593 
594c #aH     B_Enrg =(Cn_dSV      *(TsisSV(ikl,isn) -TfSnow         )
595c #aH.            -Lf_H2O      *(1.              -eta_SV(ikl,isn)))
596c #aH.           /(1.          + dzVap0 /max(epsi,dzsnSV(ikl,isn)))
597c #aH     eta_SV(ikl,isn) =
598c #aH.           max(zero,unun +(B_Enrg
599c #aH.                         -(TsisSV(ikl,isn) -TfSnow)*Cn_dSV)
600c #aH.                          /Lf_H2O                          )
601c #aH     TsisSV(ikl,isn) =    ( B_Enrg
602c #aH.                         +(1.              -eta_SV(ikl,isn))
603c #aH.                          *Lf_H2O                          )
604c #aH.                         / Cn_dSV
605c #aH.                         + TfSnow
606 
607! #e1     STOP "PLEASE add Energy (#aH) from deposition/sublimation"
608 
609 
610C +--Update of the upper Snow layer Thickness
611C +  ----------------------------------------
612 
613          dzsnSV(ikl,isn) =
614     .           max(zero,  dzsnSV(ikl,isnoSV(ikl)) + dzVap0)
615          isnoSV(ikl)     = isnoSV(ikl)             + NOLayr
616          isn             = isnoSV(ikl)
617          dzsnSV(ikl,isn) = dzsnSV(ikl,isn) + dzVap1
618          wes_SV(ikl)     = ro__SV(ikl,isn) * dzVap0
619
620      END DO
621 
622 
623C +--Energy Budget (OUT)
624C +  ===================
625 
626! #vm   DO ikl=1,knonv
627! #vm     EqSn02(ikl) =-EqSn_0(ikl)
628! #vm.                 -EExcsv(ikl)
629! #vm   END DO
630! #vm DO   isn=nsno,1,-1
631! #vm   DO ikl=1,knonv
632! #vm     EqSn02(ikl) = EqSn01(ikl) + ro__SV(ikl,isn) *dzsnSV(ikl,isn)
633! #vm.                *(Cn_dSV      *(TsisSV(ikl,isn) -TfSnow         )
634! #vm.                 -Lf_H2O      *(1.              -eta_SV(ikl,isn)))
635! #vm   END DO
636! #vm END DO
637 
638 
639C +--Snow/I Budget
640C +  -------------
641 
642! #m1   DO ikl=1,knonv
643! #m1     SIsubl(ikl) = dt__SV*HLs_sv(ikl)*min(isnoSV(ikl),1)
644! #m1.                        /Lx_H2O(ikl)
645! #m1     SIrnof(ikl) = rusnSV(ikl) + RnofSV(ikl) * dt__SV
646! #m1.                - SIrnof(ikl)
647! #m1   END DO
648 
649 
650C +--Anticipated Disappearance of a rapidly Melting too thin Last Snow Layer
651C +  =======================================================================
652 
653      DO ikl=1,knonv
654        LastOK = min(1   , max(0   ,iiceSV(ikl)-isnoSV(ikl)+2)
655     .                    *min(1   ,isnoSV(ikl)-iiceSV(ikl))
656     .                    +min(1   ,isnoSV(ikl))              )
657        RapdOK = max(zero,sign(unun,dzMelt(ikl)-epsi         ))
658        ThinOK = max(zero,sign(unun,dz_Min     -dzsnSV(ikl,1)))
659        z_Melt = LastOK     *RapdOK*ThinOK
660        noSnow(ikl)   = noSnow(ikl)   + z_Melt
661        z_Melt        =                 z_Melt *dzsnSV(ikl,1)
662        dzsnSV(ikl,1) = dzsnSV(ikl,1) - z_Melt
663        EExcsv(ikl)   = EExcsv(ikl)   - z_Melt *ro__SV(ikl,1)
664     .                                *(1.     -eta_SV(ikl,1))*Lf_H2O
665 
666C +--Water  Production
667C +  ^^^^^^^^^^^^^^^^^
668        drr_SV(ikl)   = drr_SV(ikl)
669     .                + ro__SV(ikl,1) * z_Melt /dt__SV
670      END DO
671 
672 
673C +--Update Nb of Layers
674C +  ===================
675 
676      DO ikl=1,knonv
677        isnoSV(ikl)   = isnoSV(ikl)
678     .     * min(1,iabs(isnoSV(ikl)-noSnow(ikl)))
679      END DO
680 
681 
682! Energy Budget (OUT)
683! ===================
684 
685! #e1   DO ikl=1,knonv
686! #e1     EqSn_1(ikl) = 0.
687! #e1   END DO
688! #e1 DO   isn=nsno,1,-1
689! #e1   DO ikl=1,knonv
690! #e1     EqSn_1(ikl) = EqSn_1(ikl) + ro__SV(ikl,isn) *dzsnSV(ikl,isn)
691! #e1.                *(Cn_dSV      *(TsisSV(ikl,isn) -TfSnow         )
692! #e1.                 -Lf_H2O      *(1.              -eta_SV(ikl,isn)))
693! #e1   END DO
694! #e1 END DO
695 
696 
697C +--Water  Budget (OUT)
698C +  ===================
699 
700! #vm   DO ikl=1,knonv
701! #vm     WqSn_0(ikl) = WqSn_0(ikl)
702! #vm.                + HLs_sv(ikl)    * dt__SV
703! #vm.             *min(isnoSV(ikl),1) / Lx_H2O(ikl)
704! #vm     WqSn_1(ikl) = drr_SV(ikl)    * dt__SV
705! #vm.                + rusnSV(ikl)
706! #vm.                + RnofSV(ikl)    * dt__SV
707! #vm   END DO
708! #vm DO   isn=nsno,1,-1
709! #vm   DO ikl=1,knonv
710! #vm     WqSn_1(ikl) = WqSn_1(ikl)
711! #vm.                + ro__SV(ikl,isn)* dzsnSV(ikl,isn)
712! #vm   END DO
713! #vm END DO
714 
715 
716      return
717      end
Note: See TracBrowser for help on using the repository browser.