source: LMDZ6/branches/Amaury_dev/libf/phylmd/inlandsis/sisvat_qsn.f90 @ 5158

Last change on this file since 5158 was 5117, checked in by abarral, 2 months ago

rename modules properly lmdz_*
move some unused files to obsolete/
(lint) uppercase fortran keywords
File size: 26.4 KB
Line 
1
2
3SUBROUTINE SISVAT_qSn &
4        ( &
5  ! #e1.                      EqSn_0,EqSn_1,EqSn_d
6  ! #m1.                     ,SIsubl,SImelt,SIrnof
7        )
8
9  ! +------------------------------------------------------------------------+
10  ! | MAR          SISVAT_qSn                           Fri 29-Jul-2011  MAR |
11  ! |   SubRoutine SISVAT_qSn updates  the Snow Water Content                |
12  ! +------------------------------------------------------------------------+
13  ! |                                                                        |
14  ! |   PARAMETERS:  knonv: Total Number of columns =                        |
15  ! |   ^^^^^^^^^^        = Total Number of continental     grid boxes       |
16  ! |                     X       Number of Mosaic Cell per grid box         |
17  ! |                                                                        |
18  ! |   INPUT:   isnoSV   = total Nb of Ice/Snow Layers                      |
19  ! |   ^^^^^                                                                |
20  ! |                                                                        |
21  ! |   INPUT:   TaT_SV   : SBL Top    Temperature                       [K] |
22  ! |   ^^^^^    dt__SV   : Time Step                                    [s] |
23  ! |                                                                        |
24  ! |   INPUT /  drr_SV   : Rain Intensity                         [kg/m2/s] |
25  ! |   OUTPUT:  dzsnSV   : Snow Layer Thickness                         [m] |
26  ! |   ^^^^^^   eta_SV   : Snow Water Content                       [m3/m3] |
27  ! |            ro__SV   : Snow/Soil  Volumic Mass                  [kg/m3] |
28  ! |            TsisSV   : Soil/Ice Temperatures (layers -nsol,-nsol+1,..,0)|
29  ! |                     & Snow     Temperatures (layers  1,2,...,nsno) [K] |
30  ! |                                                                        |
31  ! |   OUTPUT:  SWS_SV   : Surficial Water Status                           |
32  ! |   ^^^^^^                                                               |
33  ! |            EExcsv   : Snow Energy in Excess, initial Forcing    [J/m2] |
34  ! |            EqSn_d   : Snow Energy in Excess, remaining          [J/m2] |
35  ! |            EqSn_0   : Snow Energy, before Phase Change          [J/m2] |
36  ! |            EqSn_1   : Snow Energy, after  Phase Change          [J/m2] |
37  ! |            SIsubl   : Snow sublimed/deposed Mass             [mm w.e.] |
38  ! |            SImelt   : Snow Melted           Mass             [mm w.e.] |
39  ! |            SIrnof   : Surficial Water + Run OFF Change       [mm w.e.] |
40  ! |                                                                        |
41  ! |   Internal Variables:                                                  |
42  ! |   ^^^^^^^^^^^^^^^^^^                                                   |
43  ! |                                                                        |
44  ! | # OPTIONS: #E0: IO for Verification: Energy       Budget               |
45  ! | # ^^^^^^^                                                              |
46  ! | #          #su: IO for Verification: Slush        Diagnostic           |
47  ! |                                                                        |
48  ! |                                                                        |
49  ! +------------------------------------------------------------------------+
50
51
52
53
54  ! +--Global Variables
55  ! +  ================
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
85  ! +--Internal Variables
86  ! +  ==================
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
128  ! +--OUTPUT of SISVAT Trace Statistics (see assignation in PHY_SISVAT)
129  ! +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
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
137  ! +--Energy and Mass Budget
138  ! +  ~~~~~~~~~~~~~~~~~~~~~~
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
148  ! +--Slush Diagnostic: IO
149  ! +  ~~~~~~~~~~~~~~~~~~~~
150  ! #vu logical         su_opn                 ! IO   Switch
151  ! #vu common/SI_qSn_L/su_opn                 !
152
153
154  ! +--DATA
155  ! +  ====
156
157  data      dzepsi/0.0001/                ! Minim. Snow Layer Thickness (!)
158  ! #?? data      dz_Min/0.005/                 ! Minim. Snow Layer Thickness
159  ! ... 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
164  ! +--Energy Budget (IN)
165  ! +  ==================
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
179  ! +--Water  Budget (IN)
180  ! +  ==================
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
193  ! +--Snow Melt Budget
194  ! +  ================
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
202  ! +--Initialization
203  ! +  ==============
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
216  ! +--Melting/Freezing Energy
217  ! +  =======================
218
219  ! +...REMARK: Snow liquid Water Temperature assumed = TfSnow
220  ! +   ^^^^^^
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
230  ! +--Surficial Water Status
231  ! +  ----------------------
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
242  ! +--Energy, store Previous Content
243  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
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
250  ! +--Water,  store Previous Content
251  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
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
261  ! +--Melting  if EExcsv > 0
262  ! +  ======================
263
264      EnMelt      =    max(zero,          EExcsv(ikl) )
265
266  ! +--Energy Consumption
267  ! +  ^^^^^^^^^^^^^^^^^^
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
282  ! +--Water  Production
283  ! +  ^^^^^^^^^^^^^^^^^
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
290  ! +--Snow History
291  ! +  ^^^^^^^^^^^^
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
300  ! +--Freezing if EExcsv < 0
301  ! +  ======================
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
322  ! +--Snow Water Content
323  ! +  ==================
324
325  ! +--Percolation Velocity
326  ! +  ^^^^^^^^^^^^^^^^^^^^
327  ! #PW     SGDiam    = 1.6d-4
328  ! #PW.              + 1.1d-13 *(ro__SV(ikl,isn)*ro__SV(ikl,isn)
329  ! #PW.                         *ro__SV(ikl,isn)*ro__SV(ikl,isn))
330
331  ! +--Pore   Volume [-]
332  ! +  ^^^^^^^^^^^^^^^^^
333      rosDry      =(1.     - eta_SV(ikl,isn))* ro__SV(ikl,isn) !
334      PorVol      = 1.     - rosDry          / ro_Ice          !
335      PorVol      =      max(PorVol          , zero  )         !
336
337  ! +--Water  Retention
338  ! +  ^^^^^^^^^^^^^^^^
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
348  ! +--Pore Hole Close OFF
349  ! +  ^^^^^^^^^^^^^^^^^^^
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
359  !XF
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
380      ! if (isn>1.AND.isn<nsno     .AND.
381  !    .      ro__SV(ikl,isn-1)>900    .AND.
382  !    .      ro__SV(ikl,isn)  >roCdSV .AND.
383  !    .      ro__SV(ikl,isn)  <900    .AND.
384  !    .      TsisSV(ikl,isn)  >273.14 .AND.
385  !    .      TsisSV(ikl,isn+1)<273.15 .AND.
386  !    .      drr_SV(ikl)      >0)     THEN
387      !  TsisSV(ikl,isn)=273.14
388      !  PClose = 1
389      ! END IF
390
391  !XF
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
403  ! +--Remove Zero-Thickness Layers
404  ! +  ============================
405
406 1000   CONTINUE
407       isnitr =          0
408  DO   ikl=1,knonv
409       isnUpD =          0
410       isinew =          0
411  !XF
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>0)                                       GO TO 1000
441
442
443  ! +--New upper Limit of the non erodible Snow (istoSV .GT. 1)
444  ! +  ========================================
445
446  DO   ikl=1,knonv
447       nh =     0
448  !XF
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.
454  !XF
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
466  ! +--Energy Budget (OUT)
467  ! +  ===================
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
482  ! +--"Negative Heat" from supercooled rain
483  ! +   ------------------------------------
484
485  DO ikl=1,knonv
486      EExcsv(ikl) = EExcsv(ikl) + EExdum(ikl)
487
488
489  ! +--Surficial Water Run OFF
490  ! +  -----------------------
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
497  ! #EU                        rusnew = 0.
498  ! #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
510  ! +--Percolation down the Continental Ice Pack
511  ! +  -----------------------------------------
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
520  !XF 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
541  ! +--Slush Formation (Activated. CAUTION: ADD RunOff Possibility before Activation)
542  ! +  ---------------  ^^^^^^^  ^^^
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
550  ! +--Available Additional Pore   Volume [-]
551  ! +  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
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
576  ! +--Impact of the Sublimation/Deposition on the Surface Mass Balance
577  ! +  ================================================================
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
588  ! +--Additional Energy
589  ! +  -----------------
590
591  ! #VH     AdEnrg = dzVap0 * ro__SV(ikl,isnoSV(ikl))           ! Water   Vapor
592  ! #VH.            *C__Wat *(TsisSV(ikl,isnoSV(ikl)) -TfSnow)  ! Sensible Heat
593
594  ! #aH     B_Enrg =(Cn_dSV      *(TsisSV(ikl,isn) -TfSnow         )
595  ! #aH.            -Lf_H2O      *(1.              -eta_SV(ikl,isn)))
596  ! #aH.           /(1.          + dzVap0 /max(epsi,dzsnSV(ikl,isn)))
597  ! #aH     eta_SV(ikl,isn) =
598  ! #aH.           max(zero,unun +(B_Enrg
599  ! #aH.                         -(TsisSV(ikl,isn) -TfSnow)*Cn_dSV)
600  ! #aH.                          /Lf_H2O                          )
601  ! #aH     TsisSV(ikl,isn) =    ( B_Enrg
602  ! #aH.                         +(1.              -eta_SV(ikl,isn))
603  ! #aH.                          *Lf_H2O                          )
604  ! #aH.                         / Cn_dSV
605  ! #aH.                         + TfSnow
606
607  ! #e1     STOP "PLEASE add Energy (#aH) from deposition/sublimation"
608
609
610  ! +--Update of the upper Snow layer Thickness
611  ! +  ----------------------------------------
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
623  ! +--Energy Budget (OUT)
624  ! +  ===================
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
639  ! +--Snow/I Budget
640  ! +  -------------
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
650  ! +--Anticipated Disappearance of a rapidly Melting too thin Last Snow Layer
651  ! +  =======================================================================
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
666  ! +--Water  Production
667  ! +  ^^^^^^^^^^^^^^^^^
668    drr_SV(ikl)   = drr_SV(ikl) &
669          + ro__SV(ikl,1) * z_Melt /dt__SV
670  END DO
671
672
673  ! +--Update Nb of Layers
674  ! +  ===================
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
697  ! +--Water  Budget (OUT)
698  ! +  ===================
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
717END SUBROUTINE sisvat_qsn
Note: See TracBrowser for help on using the repository browser.