source: LMDZ6/trunk/libf/phylmd/inlandsis/sisvat_gsn.F @ 3792

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

Ajout de INLANDSIS, nouvelle interface entre LMDZ et la neige de SISVAT
Etienne, 04/01/2021

File size: 45.3 KB
Line 
1 
2      subroutine SISVAT_GSn
3 
4C +------------------------------------------------------------------------+
5C | MAR          SISVAT_GSn                                20-09-2003  MAR |
6C |   SubRoutine SISVAT_GSn simulates SNOW Metamorphism                    |
7C +------------------------------------------------------------------------+
8C |                                                                        |
9C |   PARAMETERS:  knonv: Total Number of columns =                        |
10C |   ^^^^^^^^^^        = Total Number of continental     grid boxes       |
11C |                     X       Number of Mosaic Cell per grid box         |
12C |                                                                        |
13C |   INPUT /  isnoSV   = total Nb of Ice/Snow Layers                      |
14C |   OUTPUT:  iiceSV   = total Nb of Ice      Layers                      |
15C |   ^^^^^^   istoSV   = 0,...,5 :   Snow     History (see istdSV data)   |
16C |                                                                        |
17C |   INPUT:   TsisSV   : Soil/Ice Temperatures (layers -nsol,-nsol+1,..,0)|
18C |   ^^^^^             & Snow     Temperatures (layers  1,2,...,nsno) [K] |
19C |            ro__SV   : Soil/Snow Volumic Mass                   [kg/m3] |
20C |            eta_SV   : Soil/Snow Water   Content                [m3/m3] |
21C |            slopSV   : Surface Slope                                [-] |
22C |            dzsnSV   : Snow Layer        Thickness                  [m] |
23C |            dt__SV2   : Time  Step                                   [s] |
24C |                                                                        |
25C |   INPUT /  G1snSV   : Dendricity (<0) or Sphericity (>0) of Snow Layer |
26C |   OUTPUT:  G2snSV   : Sphericity (>0) or Size            of Snow Layer |
27C |   ^^^^^^                                                               |
28C |                                                                        |
29C |   Formalisme adopte pour la Representation des Grains:                 |
30C |   Formalism         for the Representation of  Grains:                 |
31C |   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                 |
32C |                                                                        |
33C |             1       - -1                 Neige Fraiche                 |
34C |            / \      |                    -------------                 |
35C |           /   \     |  Dendricite        decrite  par Dendricite       |
36C |          /     \    |  Dendricity                  et Sphericite       |
37C |         /       \   |                                                  |
38C |        2---------3  -  0                 described by Dendricity       |
39C |                                                   and Sphericity       |
40C |        |---------|                                                     |
41C |        0         1                                                     |
42C |        Sphericite                                                      |
43C |        Sphericity                                                      |
44C |                                                                        |
45C |        4---------5  -                                                  |
46C |        |         |  |                                                  |
47C |        |         |  |  Diametre (1/10eme de mm) (ou Taille)            |
48C |        |         |  |  Diameter (1/10th  of mm) (or Size  )            |
49C |        |         |  |                                                  |
50C |        |         |  |                    Neige non dendritique         |
51C |        6---------7  -                    ---------------------         |
52C |                                          decrite  par Sphericite       |
53C |                                                    et     Taille       |
54C |                                          described by Sphericity       |
55C |                                                   and       Size       |
56C |                                                                        |
57C |   Les Variables du Modele:                                             |
58C |   Model         Variables:                                             |
59C |   ^^^^^^^^^^^^^^^^^^^^^^^^                                             |
60C |     Cas Dendritique               Cas non Dendritique                  |
61C |                                                                        |
62C |     G1snSV        : Dendricite    G1snSV        : Sphericite           |
63C |     G2snSV        : Sphericite    G2snSV        : Taille (1/10e mm)    |
64C |                                                   Size                 |
65C |                                                                        |
66C |   Cas Dendritique/ Dendritic Case                                      |
67C |   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                                      |
68C |   Dendricite(Dendricity) G1snSV                                        |
69C |            varie     de -G1_dSV (-99 par defaut / etoile)          a 0 |
70C |            division par -G1_dSV pour obtenir des valeurs entre 1  et 0 |
71C |            varies  from -G1_dSV (default -99    / fresh snow)     to 0 |
72C |            division  by -G1_dSV to obtain values       between 1 and 0 |
73C |                                                                        |
74C |   Sphericite(Sphericity) G2snSV                                        |
75C |            varie     de  0         (cas completement anguleux)         |
76C |                       a  G1_dSV (99 par defaut, cas spherique)         |
77C |            division par  G1_dSV pour obtenir des valeurs entre 0  et 1 |
78C |            varies  from  0      (full faceted)               to G1_dSV |
79C |                                                                        |
80C |   Cas non Dendritique / non Dendritic Case                             |
81C |   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                             |
82C |   Sphericite(Sphericity) G1snSV                                        |
83C |            varie     de  0         (cas completement anguleux)         |
84C |                       a  G1_dSV (99 par defaut, cas spherique)         |
85C |            division par  G1_dSV pour obtenir des valeurs entre 0  et 1 |
86C |            varies  from  0      (full faceted)               to G1_dSV |
87C |                                                                        |
88C |   Taille    (Size)       G2snSV                                        |
89C |            superieure a  ADSdSV (.4 mm) et ne fait que croitre         |
90C |            greater than  ADSdSV (.4 mm) always increases               |
91C |                                                                        |
92C |   Exemples: Points caracteristiques des Figures ci-dessus              |
93C |   ^^^^^^^^^                                                            |
94C |                                                                        |
95C |               G1snSV    G2snSV     dendricite  sphericite  taille      |
96C |                                    dendricity  sphericity  size        |
97C |   ------------------------------------------------------------------   |
98C |                                                            [1/10 mm]   |
99C |     1        -G1_dSV    sph3SN          1           0.5                |
100C |     2           0         0             0           0                  |
101C |     3           0       G1_dSV          0           1                  |
102C |     4           0       ADSdSV                      0       4.         |
103C |     5         G1_dSV    ADSdSV-vsphe1               1       3.         |
104C |     6           0         --                        0       --         |
105C |     7         G1_dSV      --                        1       --         |
106C |                                                                        |
107C |   par defaut: G1_dSV=99.                                               |
108C |                         sph3SN=50.                                     |
109C |                         ADSdSV= 4.                                     |
110C |                                vsphe1=1.                               |
111C |                                                                        |
112C |   Methode:                                                             |
113C |   ^^^^^^^^                                                             |
114C |   1. Evolution Types de Grains selon Lois de Brun et al. (1992):       |
115C |      Grain metamorphism according to         Brun et al. (1992):       |
116C |      Plusieurs Cas sont a distiguer  /  the different Cases are:       |
117C |       1.1 Metamorphose Neige humide  /  wet Snow                       |
118C |       1.2 Metamorphose Neige seche   /  dry Snow                       |
119C |         1.2.1 Gradient faible        /  low      Temperature Gradient  |
120C |         1.2.2 Gradient moyen         /  moderate Temperature Gradient  |
121C |         1.2.3 Gradient fort          /  high     Temperature Gradient  |
122C |      Dans chaque Cas on separe Neige Dendritique et non Dendritique    |
123C |                           le Passage Dendritique -> non Dendritique    |
124C |                           se fait lorsque  G1snSV devient > 0          |
125C |      the Case of Dentritic or non Dendritic Snow is treated separately |
126C |      the Limit   Dentritic -> non Dendritic is reached when G1snSV > 0 |
127C |                                                                        |
128C |   2. Tassement: Loi de Viscosite adaptee selon le Type de Grains       |
129C |      Snow Settling:    Viscosity depends on the   Grain Type           |
130C |                                                                        |
131C |   3. Update Variables historiques (cas non dendritique seulement)      |
132C |      nhSNow defaut                                                     |
133C |                0    Cas normal                                         |
134C |      istdSV(1) 1    Grains anguleux / faceted cristal                  |
135C |      istdSV(2) 2    Grains ayant ete en presence d eau liquide         |
136C |                     mais n'ayant pas eu de caractere anguleux    /     |
137C |                     liquid water and no faceted cristals before        |
138C |      istdSV(3) 3    Grains ayant ete en presence d eau liquide         |
139C |                     ayant eu auparavant un caractere anguleux    /     |
140C |                     liquid water and    faceted cristals before        |
141C |                                                                        |
142C |   REFER. : Brun et al.      1989, J. Glaciol 35 pp. 333--342           |
143C |   ^^^^^^^^ Brun et al.      1992, J. Glaciol 38 pp.  13-- 22           |
144C |            (CROCUS Model, adapted to MAR at CEN by H.Gallee)           |
145C |                                                                        |
146C |   REFER. : Marbouty, D.     1980, J. Glaciol 26 pp. xxx--xxx           |
147C |   ^^^^^^^^ (CROCUS Model, adapted to MAR at CEN by H.Gallee)           |
148C |            (for angular shapes)                                        |
149C |                                                                        |
150C |   Preprocessing  Option: SISVAT IO (not always a standard preprocess.) |
151C |   ^^^^^^^^^^^^^^^^^^^^^  ^^^^^^^^^                                     |
152C |   FILE                 |      CONTENT                                  |
153C |   ~~~~~~~~~~~~~~~~~~~~~+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
154C | # SISVAT_GSn.vp        | #vp: OUTPUT/Verification: Snow   Properties   |
155C |                        |      unit 47, SubRoutines SISVAT_zSn, _GSn    |
156C | # stdout               | #wp: OUTPUT/Verification: Snow   Properties   |
157C |                        |      unit  6, SubRoutine  SISVAT_GSn          |
158C |                                                                        |
159C +------------------------------------------------------------------------+
160 
161 
162 
163 
164C +--Global Variables
165C +  ================
166
167      use VARphy
168      use VAR_SV
169      use VARdSV
170      use VAR0SV
171      use VARxSV
172      use VARtSV
173
174 
175      IMPLICIT NONE
176 
177     
178 
179C +--INPUT/OUTPUT
180C +  ------------
181 
182 
183C +--OUTPUT
184C +  ------
185 
186      integer   dt__SV2
187 
188 
189C +--Local  Variables
190C +  ================
191 
192      logical   vector                        !
193      integer   ikl                           !
194      integer   isn   ,isnp                   !
195      integer   istoOK                        !
196      real      G1_bak,G2_bak                 ! Old Values of G1, G2
197      real      ro_dry(knonv,      nsno)      ! Dry Density            [g/cm3]
198      real      etaSno(knonv,      nsno)      ! Liquid Water Content   [g/cm2]
199      real      SnMass(knonv)                 ! Snow   Mass            [kg/m2]
200      real      dTsndz                        ! Temperature Gradient
201      real      sWater                        !        Water Content       [%]
202      real      exp1Wa                        !
203      real      dDENDR                        ! Dendricity Increment
204      real      DENDRn                        ! Normalized Dendricity
205      real      SPHERn                        ! Normalized Sphericity
206      real      Wet_OK                        ! Wet Metamorphism Switch
207      real      OK__DE                        !
208      real      OK__wd                        ! New G*, from wet Dendritic
209      real      G1__wd                        ! New G1, from wet Dendritic
210      real      G2__wd                        ! New G2, from wet Dendritic
211      real      OKlowT                        !
212      real      facVap                        !
213      real      OK_ldd                        !
214      real      G1_ldd                        !
215      real      G2_ldd                        !
216      real      DiamGx                        !
217      real      DiamOK                        !
218      real      No_Big                        !
219      real      dSPHER                        !
220      real      SPHER0                        !
221      real      SPHbig                        !
222      real      G1_lds                        !
223      real      OK_mdT                        !
224      real      OKmidT                        !
225      real      OKhigT                        !
226      real      OK_mdd                        !
227      real      G1_mdd                        !
228      real      G2_mdd                        !
229      real      G1_mds                        !
230      real      OK_hdd                        !
231      real      G1_hdd                        !
232      real      G2_hdd                        !
233      real      OK_hds                        !
234      real      G1_hds                        !
235      real      T1__OK,T2__OK                 !
236      real      T3_xOK,T3__OK,T3_nOK          !
237      real      ro1_OK,ro2_OK                 !
238      real      dT1_OK,dT2_OK,dT3xOK,dT3_OK   !
239      real      dT4xOK,dT4_OK,dT4nOK,AngSno   !
240      real      G2_hds,SphrOK,HISupd          !
241      real      H1a_OK,H1b_OK,H1__OK          !
242      real      H23aOK,H23bOK,H23_OK          !
243      real      H2__OK,H3__OK                 !
244      real      H45_OK,H4__OK,H5__OK          !
245      real      ViscSn,OK_Liq,OK_Ang,OKxLiq   !
246      real      dSnMas,dzsnew,rosnew,rosmax,smb_old,smb_new
247      real      zn_old,zn_new
248 
249      real      epsi5                         ! Alpha ev67 single precision     
250      real      vdiam1                        ! Small Grains Min.Diam.[.0001m]
251      real      vdiam2                        ! Spher.Variat.Max Diam.    [mm]
252      real      vdiam3                        ! Min.Diam.|Limit Spher.    [mm]
253      real      vdiam4                        ! Min.Diam.|Viscosity Change
254      real      vsphe1                        ! Max Sphericity
255      real      vsphe2                        ! Low T Metamorphism  Coeff.
256      real      vsphe3                        ! Max.Sphericity (history=1)
257      real      vsphe4                        ! Min.Sphericity=>history=1
258      real      vtang1,vtang2,vtang3,vtang4   ! Temperature Contribution
259      real      vtang5,vtang6,vtang7,vtang8   !
260      real      vtang9,vtanga,vtangb,vtangc   !
261      real      vrang1,vrang2                 ! Density     Contribution
262      real      vgang1,vgang2,vgang3,vgang4   ! Grad(T)     Contribution
263      real      vgang5,vgang6,vgang7,vgang8   !
264      real      vgang9,vganga,vgangb,vgangc   !
265      real      vgran6                        ! Max.Sphericity for Settling
266      real      vtelv1                        ! Threshold | history = 2, 3
267      real       vvap1                        ! Vapor Pressure Coefficient
268      real       vvap2                        ! Vapor Pressure Exponent
269      real      vgrat1                        ! Boundary weak/mid   grad(T)
270      real      vgrat2                        ! Boundary mid/strong grad(T)
271      real         vfi                        ! PHI,         strong grad(T)
272      real      vvisc1,vvisc2,vvisc3,vvisc4   ! Viscosity Coefficients
273      real      vvisc5,vvisc6,vvisc7          ! id., wet Snow
274      real      rovisc                        ! Wet Snow Density  Influence
275      real        vdz3                        ! Maximum Layer Densification
276      real      OK__ws                        ! New G2
277      real      G1__ws                        ! New G1, from wet Spheric
278      real      G2__ws                        ! New G2, from wet Spheric
279      real      husi_0,husi_1,husi_2,husi_3   ! Constants for New G2
280      real      vtail1,vtail2                 ! Constants for New G2
281      real      frac_j                        ! Time Step            [Day]
282 
283      real       vdent1                       ! Wet Snow Metamorphism
284      integer   nvdent1                       ! (Coefficients for
285      integer   nvdent2                       !           Dendricity)
286 
287C +--Snow Properties: IO
288C +  ~~~~~~~~~~~~~~~~~~~
289! #vp real      G_curr(18),Gcases(18)
290! #vp common   /GSnLOC/    Gcases
291! #wp real                 D__MAX
292! #wp common   /GSnMAX/    D__MAX
293 
294 
295C +--DATA
296C +  ====
297 
298      data       vector/.true./               ! Vectorization  Switch
299      data       vdent1/ 0.5e8/               ! Wet Snow Metamorphism
300cXF                      tuned for Greenland (2.e8=old value)
301      data      nvdent1/ 3   /                ! (Coefficients for
302      data      nvdent2/16   /                !           Dendricity)
303 
304      data       husi_0 /20.      /           !   10  * 2
305      data       husi_1 / 0.23873 /           ! (3/4) /pi
306      data       husi_2 / 4.18880 /           ! (4/3) *pi
307      data       husi_3 / 0.33333 /           !  1/3
308      data       vtail1 / 1.28e-08/           !  Wet Metamorphism
309      data       vtail2 / 4.22e-10/           ! (NON Dendritic / Spheric)
310 
311      data       epsi5  / 1.0e-5  /           !
312 
313      data       vdiam1 / 4.0     /           ! Small Grains Min.Diameter
314 
315      data       vdiam2 / 0.5     /           ! Spher.Variat.Max Diam.[mm]
316      data       vdiam3 / 3.0     /           ! Min.Diam.|Limit Spher.[mm]
317      data       vdiam4 / 2.0     /           ! Min.Diam.|Viscosity Change
318 
319      data       vsphe1 / 1.0     /           ! Max Sphericity
320      data       vsphe2 / 1.0e9   /           ! Low T Metamorphism  Coeff.
321      data       vsphe3 / 0.5     /           ! Max.Sphericity (history=1)
322      data       vsphe4 / 0.1     /           ! Min.Sphericity=>history=1
323 
324      data       vgran6 / 51.     /           ! Max.Sphericity for Settling
325      data       vtelv1 / 5.e-1   /           ! Threshold | history = 2, 3
326 
327      data        vvap1 /-6.e3    /           ! Vapor Pressure Coefficient
328      data        vvap2 / 0.4     /           ! Vapor Pressure Exponent
329 
330      data       vgrat1 /0.05     /           ! Boundary weak/mid   grad(T)
331      data       vgrat2 /0.15     /           ! Boundary mid/strong grad(T)
332      data          vfi /0.09     /           ! PHI,         strong grad(T)
333 
334      data       vvisc1 / 0.70    /           ! Viscosity Coefficients
335      data       vvisc2 / 1.11e5  /           !
336      data       vvisc3 /23.00    /           !
337      data       vvisc4 / 0.10    /           !
338      data       vvisc5 / 1.00    /           ! id., wet Snow
339      data       vvisc6 / 2.00    /           !
340      data       vvisc7 /10.00    /           !
341      data       rovisc / 0.25    /           ! Wet Snow Density  Influence
342      data         vdz3 / 0.30    /           ! Maximum Layer Densification
343 
344 
345C +--DATA (Coefficient Fonction fort Gradient Marbouty)
346C +  --------------------------------------------------
347 
348      data       vtang1 /40.0/                ! Temperature Contribution
349      data       vtang2 / 6.0/                !
350      data       vtang3 /22.0/                !
351      data       vtang4 / 0.7/                !
352      data       vtang5 / 0.3/                !
353      data       vtang6 / 6.0/                !
354      data       vtang7 / 1.0/                !
355      data       vtang8 / 0.8/                !
356      data       vtang9 /16.0/                !
357      data       vtanga / 0.2/                !
358      data       vtangb / 0.2/                !
359      data       vtangc /18.0/                !
360 
361      data       vrang1 / 0.40/               ! Density     Contribution
362      data       vrang2 / 0.15/               !
363 
364      data       vgang1 / 0.70/               ! Grad(T)     Contribution
365      data       vgang2 / 0.25/               !
366      data       vgang3 / 0.40/               !
367      data       vgang4 / 0.50/               !
368      data       vgang5 / 0.10/               !
369      data       vgang6 / 0.15/               !
370      data       vgang7 / 0.10/               !
371      data       vgang8 / 0.55/               !
372      data       vgang9 / 0.65/               !
373      data       vganga / 0.20/               !
374      data       vgangb / 0.85/               !
375      data       vgangc / 0.15/               !
376 
377! #wp data       D__MAX / 4.00/               !
378 
379 
380C +-- 1. Metamorphoses dans les Strates
381C +      Metamorphism
382C +      ==============================
383 
384      dt__SV2= dt__SV
385      frac_j = dt__SV2 / 86400.                        ! Time Step [Day]
386 
387      zn4_SV = 0
388 
389 
390C +-- 1.1 Initialisation: teneur en eau liquide et gradient de temperature
391C +   ------------------  liquid water content and temperature gradient
392 
393       DO ikl=1,knonv
394        DO   isn=1,isnoSV(ikl)
395 
396          ro_dry(ikl,isn) = 1.e-3 *ro__SV(ikl,isn)    ! Dry Density
397     .                    *(1.    -eta_SV(ikl,isn))   !         [g/cm3]
398          etaSno(ikl,isn) = 1.e-1 *dzsnSV(ikl,isn)    ! Liquid Water
399     .                    *        ro__SV(ikl,isn)    ! Content [g/cm2]
400     .                    *        eta_SV(ikl,isn)    !
401        END DO
402      END DO
403 
404c!$OMP PARALLEL DO default(firstprivate)
405c!$OMP.shared (/xSISVAT_I/,/xSISVAT_R/,/SoR0SV/,/SoI0SV/,/Sn_dSV/)
406       DO ikl=1,knonv
407        DO   isn=1,isnoSV(ikl)
408          isnp   = min(isn+1,isnoSV(ikl))
409 
410          dTsndz = abs( (TsisSV(ikl,isnp)-TsisSV(ikl,isn-1)) *2.e-2
411     .            /max(((dzsnSV(ikl,isnp)+dzsnSV(ikl,isn)  )
412     .                 *(           isnp -           isn)
413     .                 +(dzsnSV(ikl,isn )+dzsnSV(ikl,isn-1))),epsi))
414C +...    Factor 1.d-2 for Conversion K/m --> K/cm
415 
416 
417C +-- 1.2 Metamorphose humide
418C +       Wet Snow Metamorphism
419C +       ---------------------
420 
421          Wet_OK = max(zero,sign(unun,eta_SV(ikl,isn)-epsi))
422 
423 
424C +--     Vitesse de diminution de la dendricite
425C +       Rate of the dendricity decrease
426C +       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
427          sWater=1.d-1*ro__SV(ikl,isn)*eta_SV(ikl,isn)
428     .       /max(epsi,ro_dry(ikl,isn))
429C +...    sWater:Water Content [%]
430C +              1.d-1= 1.d2(1->%) * 1.d-3(ro__SV*eta_SV:kg/m3->g/cm3)
431 
432          exp1Wa=   sWater**nvdent1
433          dDENDR=max(exp1Wa/nvdent2,vdent1*exp(vvap1/TfSnow))
434 
435C +-- 1.2.1 Cas dendritique/dendritic Case
436C +         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
437          OK__wd=max(zero,                            !
438     .               sign(unun,-G1snSV(ikl,isn)       !
439     .                         -epsi           ))     !
440 
441          DENDRn=-G1snSV(ikl,isn)/G1_dSV  ! Normalized Dendricity (+)
442          SPHERn= G2snSV(ikl,isn)/G1_dSV  ! Normalized Sphericity
443          DENDRn= DENDRn -dDENDR *frac_j  ! New        Dendricity (+)
444          SPHERn= SPHERn +dDENDR *frac_j  ! New        Sphericity
445 
446          OK__DE=max(zero,                            ! IF 1.,
447     .               sign(unun, DENDRn                ! NO change
448     .                         -epsi           ))     ! Dendr. -> Spheric
449 
450          G1__wd=OK__DE *    (      -DENDRn*G1_dSV)   ! Dendritic
451     .      +(1.-OK__DE)* min(G1_dSV,SPHERn*G1_dSV)   ! Dendr. -> Spheric
452          G2__wd=OK__DE * min(G1_dSV,SPHERn*G1_dSV)   ! Spheric
453     .      +(1.-OK__DE)*(ADSdSV-min(SPHERn,vsphe1))  ! Spher. -> Size
454 
455C +-- 1.2.2 Cas non dendritique non completement spherique
456C +         Evolution de la Sphericite seulement.
457C +         Non dendritic and not completely spheric Case
458C +         Evolution of    Sphericity only (not size)
459C +         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
460          OK__ws=max(zero,                            !
461     .               sign(unun, G1_dSV                !
462     .                         -epsi5                 !
463     .                         -G1snSV(ikl,isn)))     !
464 
465          SPHERn= G1snSV(ikl,isn)/G1_dSV
466          SPHERn= SPHERn +dDENDR *frac_j
467          G1__ws=         min(G1_dSV,SPHERn*G1_dSV)
468 
469C +-- 1.2.3 Cas non dendritique et spherique / non dendritic and spheric
470C +         Evolution de la Taille seulement / Evolution of Size only
471C +         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
472          G2__ws =  husi_0
473     .           *( husi_1
474     .            *(husi_2 *( G2snSV(ikl,isn)/husi_0)**3
475     .                      +(vtail1 +vtail2 *exp1Wa    )*dt__SV2))
476     .           ** husi_3
477 
478 
479C +-- 1.3 Metamorposes seches / Dry Metamorphism
480C +       --------------------------------------
481 
482 
483C +-- 1.3.1 Calcul Metamorphose faible/low Gradient (0.00-0.05 deg/cm)
484C +         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
485          OKlowT=max(zero,                            !
486     .               sign(unun, vgrat1                !
487     .                         -dTsndz         ))     !
488 
489          facVap=exp(vvap1/TsisSV(ikl,isn))
490 
491C +-- 1.3.1.1 Cas dendritique / dendritic Case
492 
493          OK_ldd=max(zero,                            !
494     .               sign(unun,-G1snSV(ikl,isn)       !
495     .                         -epsi           ))     !
496 
497          DENDRn=-G1snSV(ikl,isn)     /G1_dSV
498          SPHERn= G2snSV(ikl,isn)     /G1_dSV
499          DENDRn= DENDRn-vdent1*facVap*frac_j
500          SPHERn= SPHERn+vsphe2*facVap*frac_j
501 
502          OK__DE=max(zero,                            ! IF 1.,
503     .               sign(unun, DENDRn                ! NO change
504     .                         -epsi           ))     ! Dendr. -> Spheric
505 
506          G1_ldd= OK__DE *    (      -DENDRn*G1_dSV)  ! Dendritic
507     .       +(1.-OK__DE)* min(G1_dSV,SPHERn*G1_dSV)  ! Dendr. -> Spheric
508          G2_ldd= OK__DE * min(G1_dSV,SPHERn*G1_dSV)  ! Spheric
509     .       +(1.-OK__DE)*(ADSdSV-min(SPHERn,vsphe1)) ! Spher. -> Size
510 
511C +-- 1.3.1.2 Cas non dendritique / non dendritic Case
512 
513          SPHERn=G1snSV(ikl,isn)/G1_dSV
514          DiamGx=G2snSV(ikl,isn)*0.1
515 
516          istoOK=min( abs(istoSV(ikl,isn)-
517     .                    istdSV(1)      ),1)         ! zero if istoSV = 1
518          DiamOK=max(zero,  sign(unun,vdiam2-DiamGx))
519          No_Big=    istoOK+DiamOK
520          No_Big=min(No_Big,unun)
521 
522          dSPHER=           vsphe2*facVap*frac_j      !
523          SPHER0=    SPHERn+dSPHER                    ! small grains
524          SPHbig=    SPHERn+dSPHER                    ! big   grains
525     .        *exp(min(zero,vdiam3-G2snSV(ikl,isn)))  ! (history = 2 or 3)
526          SPHbig=       min(vsphe3,SPHbig)            ! limited sphericity
527          SPHERn= No_Big *  SPHER0
528     .      + (1.-No_Big)*  SPHbig
529 
530          G1_lds=       min(G1_dSV,SPHERn*G1_dSV)
531 
532C +-- 1.3.2 Calcul Metamorphose Gradient Moyen/Moderate (0.05-0.15)
533C +         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
534          OK_mdT=max(zero,                            !
535     .               sign(unun, vgrat2                !
536     .                         -dTsndz))              !
537          OKmidT=               OK_mdT  *(1.-OKlowT)  !
538          OKhigT=          (1. -OK_mdT) *(1.-OKlowT)  !
539 
540          facVap=vdent1*exp(vvap1/TsisSV(ikl,isn))
541     .                 *   (1.e2 *dTsndz)**vvap2
542 
543C +-- 1.3.2.1 cas dendritique / dendritic case.
544 
545          OK_mdd=max(zero,                            !
546     .               sign(unun,-G1snSV(ikl,isn)       !
547     .                         -epsi           ))     !
548 
549          DENDRn=-G1snSV(ikl,isn)/G1_dSV
550          SPHERn= G2snSV(ikl,isn)/G1_dSV
551          DENDRn= DENDRn - facVap*frac_j
552          SPHERn= SPHERn - facVap*frac_j
553 
554          OK__DE=max(zero,                            ! IF 1.,
555     .               sign(unun, DENDRn                ! NO change
556     .                         -epsi           ))     ! Dendr. -> Spheric
557 
558          G1_mdd= OK__DE *    (      -DENDRn*G1_dSV)  ! Dendritic
559     .       +(1.-OK__DE)* max(zero  ,SPHERn*G1_dSV)  ! Dendr. -> Spheric
560          G2_mdd= OK__DE * max(zero  ,SPHERn*G1_dSV)  ! Spheric
561     .       +(1.-OK__DE)*(ADSdSV-max(SPHERn,zero  )) ! Spher. -> Size
562 
563C +-- 1.3.2.2 Cas non dendritique / non dendritic Case
564 
565          SPHERn=G1snSV(ikl,isn)/G1_dSV
566          SPHERn=         SPHERn-facVap*frac_j
567          G1_mds=max(zero,SPHERn*G1_dSV)
568 
569C +-- 1.3.3 Calcul Metamorphose fort / high Gradient
570C +   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
571          facVap=vdent1*exp(vvap1/TsisSV(ikl,isn))
572     .                 *   (1.e2 *dTsndz)**vvap2
573 
574C +-- 1.3.3.1 Cas dendritique / dendritic Case
575 
576          OK_hdd=max(zero,                            !
577     .               sign(unun,-G1snSV(ikl,isn)       !
578     .                         -epsi           ))     !
579 
580          DENDRn=-G1snSV(ikl,isn)/G1_dSV              !
581          SPHERn= G2snSV(ikl,isn)/G1_dSV              !
582          DENDRn= DENDRn - facVap*frac_j              !
583          SPHERn= SPHERn - facVap*frac_j              ! Non dendritic
584C +                                                   ! and angular
585          OK__DE=max(zero,                            ! IF 1.,
586     .               sign(unun, DENDRn                ! NO change
587     .                         -epsi  ))              ! Dendr. -> Spheric
588 
589          G1_hdd= OK__DE *    (      -DENDRn*G1_dSV)  ! Dendritic
590     .       +(1.-OK__DE)* max(zero  ,SPHERn*G1_dSV)  ! Dendr. -> Spheric
591          G2_hdd= OK__DE * max(zero  ,SPHERn*G1_dSV)  ! Spheric
592     .       +(1.-OK__DE)*(ADSdSV-max(SPHERn,zero  )) ! Spher. -> Size
593 
594C +-- 1.3.3.2 Cas non dendritique non completement anguleux.
595C +           non dendritic and spericity gt. 0
596 
597          OK_hds=max(zero,                            !
598     .               sign(unun, G1snSV(ikl,isn)       !
599     .                         -epsi           ))     !
600 
601          SPHERn= G1snSV(ikl,isn)/G1_dSV
602          SPHERn= SPHERn - facVap*frac_j
603          G1_hds= max(zero,SPHERn*G1_dSV)
604 
605C +-- 1.3.3.3 Cas non dendritique et anguleux
606C +           dendritic and spericity = 0.
607 
608          T1__OK = max(zero,sign(unun,TsisSV(ikl,isn)-TfSnow+vtang1))
609          T2__OK = max(zero,sign(unun,TsisSV(ikl,isn)-TfSnow+vtang2))
610          T3_xOK = max(zero,sign(unun,TsisSV(ikl,isn)-TfSnow+vtang3))
611          T3__OK =                    T3_xOK  * (1. - T2__OK)
612          T3_nOK =              (1. - T3_xOK) * (1. - T2__OK)
613          ro1_OK = max(zero,sign(unun,vrang1-ro_dry(ikl,isn)))
614          ro2_OK = max(zero,sign(unun,ro_dry(ikl,isn)-vrang2))
615          dT1_OK = max(zero,sign(unun,vgang1-dTsndz         ))
616          dT2_OK = max(zero,sign(unun,vgang2-dTsndz         ))
617          dT3xOK = max(zero,sign(unun,vgang3-dTsndz         ))
618          dT3_OK =                    dT3xOK  * (1. - dT2_OK)
619          dT4xOK = max(zero,sign(unun,vgang4-dTsndz         ))
620          dT4_OK =                    dT4xOK  * (1. - dT3_OK)
621     .                                        * (1. - dT2_OK)
622          dT4nOK =              (1. - dT4xOK) * (1. - dT3_OK)
623     .                                        * (1. - dT2_OK)
624 
625C +-- Influence de la Temperature /Temperature Influence
626C +   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
627          AngSno =
628     .      T1__OK                                                       ! 11
629     .    *(T2__OK*(vtang4+vtang5*(TfSnow       -TsisSV(ikl,isn))        ! 12
630     .                    /vtang6)                                       !
631     .     +T3__OK*(vtang7-vtang8*(TfSnow-vtang2-TsisSV(ikl,isn))        ! 13
632     .                    /vtang9)                                       !
633     .     +T3_nOK*(vtanga-vtangb*(TfSnow-vtang3-TsisSV(ikl,isn))        ! 14
634     .                    /vtangc))                                      !
635 
636C +-- Influence de la Masse Volumique /Density Influence
637C +   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
638     .    * ro1_OK
639     .        *(   ro2_OK*(1. - (ro_dry(ikl,isn)-vrang2)                 !
640     .                                  /(vrang1-vrang2))                !
641     .         +1.-ro2_OK                                )               !
642 
643C +-- Influence du Gradient de Temperature /Temperature Gradient Influence
644C +   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
645     .        *(   dT1_OK*(dT2_OK*vgang5*(dTsndz-vgang6)                 ! 15
646     .                                  /(vgang2-vgang6)                 !
647     .                    +dT3_OK*vgang7                                 ! 16
648     .                    +dT4_OK*vgang9                                 ! 17
649     .                    +dT4nOK*vgangb                )                ! 18
650     .         +1.-dT1_OK                                )               !
651     .    + ro1_OK
652     .        *    dT1_OK*(dT3_OK*vgang8*(dTsndz-vgang2)
653     .                                  /(vgang3-vgang2)
654     .                    +dT4_OK*vganga*(dTsndz-vgang3)
655     .                                  /(vgang4-vgang3)
656     .                    +dT4nOK*vgangc*(dTsndz-vgang4)
657     .                                  /(vgang1-vgang4))
658 
659          G2_hds = G2snSV(ikl,isn) + 1.d2 *AngSno*vfi     *frac_j
660 
661 
662C +--New Properties
663C +  --------------
664 
665          G1_bak          = G1snSV(ikl,isn)
666          G2_bak          = G2snSV(ikl,isn)
667 
668          G1snSV(ikl,isn) = Wet_OK * (    OK__wd             *G1__wd    !  1
669     .                               +(1.-OK__wd)*    OK__ws *G1__ws    !  2
670     .                               +(1.-OK__wd)*(1.-OK__ws)*G1_bak)   !  3
671     .               +(1. - Wet_OK)                                     !
672     .                *(    OKlowT  *(    OK_ldd             *G1_ldd    !  4
673     .                               +(1.-OK_ldd)            *G1_lds)   !  5
674     .                    + OKmidT  *(    OK_mdd             *G1_mdd    !  6
675     .                               +(1.-OK_mdd)            *G1_mds)   !  7
676     .                    + OKhigT  *(    OK_hdd             *G1_hdd    !  8
677     .                               +(1.-OK_hdd)*    OK_hds *G1_hds    !  9
678     .                               +(1.-OK_hdd)*(1.-OK_hds)*G1_bak))  ! 10
679 
680cXF
681      if(G1snSV(ikl,isn)<0.1)
682     .    G2_hds = G2snSV(ikl,isn) + 1.d1 *AngSno*vfi     *frac_j
683cXF
684 
685 
686          G2snSV(ikl,isn) = Wet_OK * (    OK__wd             *G2__wd    !  1
687     .                               +(1.-OK__wd)*    OK__ws *G2_bak    !  2
688     .                               +(1.-OK__wd)*(1.-OK__ws)*G2__ws)   !  3
689     .               +(1. - Wet_OK)                                     !
690     .                *(    OKlowT  *(    OK_ldd             *G2_ldd    !  4
691     .                               +(1.-OK_ldd)            *G2_bak)   !  5
692     .                    + OKmidT  *(    OK_mdd             *G2_mdd    !  6
693     .                               +(1.-OK_mdd)            *G2_bak)   !  7
694     .                    + OKhigT  *(    OK_hdd             *G2_hdd    !  8
695     .                               +(1.-OK_hdd)*    OK_hds *G2_bak    !  9
696     .                               +(1.-OK_hdd)*(1.-OK_hds)*G2_hds))  ! 10
697 
698C +--Snow Properties: IO Set Up
699C +  ~~~~~~~~~~~~~~~~~~~~~~~~~~
700! #vp     G_curr( 1) =     Wet_OK             *    OK__wd
701! #vp     G_curr( 2) =     Wet_OK             *(1.-OK__wd)*    OK__ws
702! #vp     G_curr( 3) =     Wet_OK             *(1.-OK__wd)*(1.-OK__ws)
703! #vp     G_curr( 4) = (1.-Wet_OK)*    OKlowT *    OK_ldd
704! #vp     G_curr( 5) = (1.-Wet_OK)*    OKlowT *(1.-OK_ldd)
705! #vp     G_curr( 6) = (1.-Wet_OK)*    OKmidT *    OK_mdd
706! #vp     G_curr( 7) = (1.-Wet_OK)*    OKmidT *(1.-OK_mdd)
707! #vp     G_curr( 8) = (1.-Wet_OK)*    OKhigT *    OK_hdd
708! #vp     G_curr( 9) = (1.-Wet_OK)*    OKhigT *(1.-OK_hdd)*    OK_hds
709! #vp     G_curr(10) = (1.-Wet_OK)*    OKhigT *(1.-OK_hdd)*(1.-OK_hds)
710! #vp     G_curr(11) =     T1__OK                         * G_curr(10)
711! #vp     G_curr(12) =     T2__OK                         * G_curr(10)
712! #vp     G_curr(13) =     T3__OK                         * G_curr(10)
713! #vp     G_curr(14) =     T3_nOK                         * G_curr(10)
714! #vp     G_curr(15) =     ro1_OK*     dT1_OK *    dT2_OK * G_curr(10)
715! #vp     G_curr(16) =     ro1_OK*     dT1_OK *    dT3_OK * G_curr(10)
716! #vp     G_curr(17) =     ro1_OK*     dT1_OK *    dT4_OK * G_curr(10)
717! #vp     G_curr(18) =     ro1_OK*     dT1_OK *    dT4nOK * G_curr(10)
718 
719! #vp     Gcases( 1) = max(Gcases( 1),G_curr( 1))
720! #vp     Gcases( 2) = max(Gcases( 2),G_curr( 2))
721! #vp     Gcases( 3) = max(Gcases( 3),G_curr( 3))
722! #vp     Gcases( 4) = max(Gcases( 4),G_curr( 4))
723! #vp     Gcases( 5) = max(Gcases( 5),G_curr( 5))
724! #vp     Gcases( 6) = max(Gcases( 6),G_curr( 6))
725! #vp     Gcases( 7) = max(Gcases( 7),G_curr( 7))
726! #vp     Gcases( 8) = max(Gcases( 8),G_curr( 8))
727! #vp     Gcases( 9) = max(Gcases( 9),G_curr( 9))
728! #vp     Gcases(10) = max(Gcases(10),G_curr(10))
729! #vp     Gcases(11) = max(Gcases(11),G_curr(11))
730! #vp     Gcases(12) = max(Gcases(12),G_curr(12))
731! #vp     Gcases(13) = max(Gcases(13),G_curr(13))
732! #vp     Gcases(14) = max(Gcases(14),G_curr(14))
733! #vp     Gcases(15) = max(Gcases(15),G_curr(15))
734! #vp     Gcases(16) = max(Gcases(16),G_curr(16))
735! #vp     Gcases(17) = max(Gcases(17),G_curr(17))
736! #vp     Gcases(18) = max(Gcases(18),G_curr(18))
737 
738C +--Snow Properties: IO
739C +  ~~~~~~~~~~~~~~~~~~~
740! #vp     IF          (isn    .le.     isnoSV(ikl))
741! #vp.    write(47,471)isn            ,isnoSV(ikl)                    ,
742! #vp.                 TsisSV(ikl,isn),ro__SV(ikl,isn),eta_SV(ikl,isn),
743! #vp.                 G1_bak         ,G2_bak         ,istoSV(ikl,isn),
744! #vp.                 dTsndz,
745! #vp.                (       k ,k=1,18),
746! #vp.                (G_curr(k),k=1,18),
747! #vp.                (Gcases(k),k=1,18),
748! #vp.                 Wet_OK,OK__wd,G1__wd,G2__wd,
749! #vp.                     1.-OK__wd,OK__ws,G1__ws,1.-OK__ws,G2__ws,
750! #vp.              1.-Wet_OK,OKlowT,OK_ldd,G1_ldd,          G2_ldd,
751! #vp.                            1.-OK_ldd,G1_lds,
752! #vp.                        OKmidT,OK_mdd,G1_mdd,          G1_mdd,
753! #vp.                            1.-OK_mdd,G1_mds,
754! #vp.                        OKhigT,OK_hdd,G1_hdd,          G2_hdd,
755! #vp.                            1.-OK_hdd,OK_hds,          G1_hds,
756! #vp.                                             1.-OK_hds,G2_hds,
757! #vp.                 G1snSV(ikl,isn),
758! #vp.                 G2snSV(ikl,isn)
759
760        END DO
761      END DO
762c!$OMP END PARALLEL DO
763 
764C +-- 2. Mise a Jour Variables Historiques (Cas non dendritique)
765C +      Update of the historical Variables
766C +      =======================================================
767 
768      IF (vector)                                                   THEN
769cXF
770         DO  ikl=1,knonv
771          DO isn=1,isnoSV(ikl)
772          SphrOK = max(zero,sign(unun,       G1snSV(ikl,isn)))
773          H1a_OK = max(zero,sign(unun,vsphe4-G1snSV(ikl,isn)))
774          H1b_OK =     1   - min(1   ,       istoSV(ikl,isn))
775          H1__OK =                    H1a_OK*H1b_OK
776          H23aOK = max(zero,sign(unun,vsphe4-G1_dSV
777     .                                      +G1snSV(ikl,isn)))
778          H23bOK = max(zero,sign(unun,etaSno(ikl,isn)
779     .                      /max(epsi,dzsnSV(ikl,isn))
780     .                                      -vtelv1         ))
781          H23_OK =                    H23aOK*H23bOK
782          H2__OK =     1   - min(1   ,       istoSV(ikl,isn))
783          H3__OK =     1   - min(1   ,   abs(istoSV(ikl,isn)-istdSV(1)))
784          H45_OK = max(zero,sign(unun,TfSnow-TsisSV(ikl,isn)+epsi))
785          H4__OK =     1   - min(1   ,   abs(istoSV(ikl,isn)-istdSV(2)))
786          H5__OK =     1   - min(1   ,   abs(istoSV(ikl,isn)-istdSV(3)))
787 
788          HISupd          =
789     .    SphrOK*(H1__OK                             *istdSV(1)
790     .       +(1.-H1__OK)*    H23_OK         *(H2__OK*istdSV(2)
791     .                                        +H3__OK*istdSV(3))
792     .       +(1.-H1__OK)*(1.-H23_OK) *H45_OK*(H4__OK*istdSV(4)
793     .                                        +H5__OK*istdSV(5)))
794          istoSV(ikl,isn) =   HISupd  +
795     .           (1.-min(unun,HISupd))               *istoSV(ikl,isn)
796        END DO
797        END DO
798      ELSE
799 
800 
801C +-- 2. Mise a Jour Variables Historiques (Cas non dendritique)
802C +      Update of the historical Variables
803C +      =======================================================
804 
805        DO ikl=1,knonv
806        DO isn=iiceSV(ikl),isnoSV(ikl)
807          IF  (G1snSV(ikl,isn).ge.0.)                               THEN
808            IF(G1snSV(ikl,isn).lt.vsphe4.and.istoSV(ikl,isn).eq.0)  THEN
809                   istoSV(ikl,isn)=istdSV(1)
810            ELSEIF(G1_dSV-G1snSV(ikl,isn)         .lt.vsphe4.and.
811     .             etaSno(ikl,isn)/dzsnSV(ikl,isn).gt.vtelv1)       THEN
812              IF  (istoSV(ikl,isn).eq.0)
813     .             istoSV(ikl,isn)=   istdSV(2)
814              IF  (istoSV(ikl,isn).eq.istdSV(1))
815     .             istoSV(ikl,isn)=   istdSV(3)
816            ELSEIF(TsisSV(ikl,isn).lt.TfSnow)                       THEN
817              IF  (istoSV(ikl,isn).eq.istdSV(2))
818     .             istoSV(ikl,isn)=   istdSV(4)
819              IF  (istoSV(ikl,isn).eq.istdSV(3))
820     .             istoSV(ikl,isn)=   istdSV(5)
821            END IF
822          END IF
823        END DO
824        END DO
825      END IF
826 
827 
828C +-- 3. Tassement mecanique /mechanical Settlement
829C +      ==========================================
830 
831        DO ikl=1,knonv
832          SnMass(ikl) = 0.
833        END DO
834cXF
835      DO ikl=1,knonv
836 
837        smb_old = 0.
838         zn_old = 0
839        DO isn  = 1, isnoSV(ikl)
840        smb_old = smb_old + dzsnSV(ikl,isn) *ro__SV(ikl,isn)
841         zn_old = zn_old  + dzsnSV(ikl,isn)
842        ENDDO
843 
844        DO   isn=isnoSV(ikl),1,-1
845          dSnMas     = 100.*dzsnSV(ikl,isn)*ro_dry(ikl,isn)
846          SnMass(ikl)=      SnMass(ikl)+0.5*dSnMas
847          ViscSn     =      vvisc1         *vvisc2
848     .               *exp(vvisc3           *ro_dry(ikl,isn)
849     .                   +vvisc4*abs(TfSnow-TsisSV(ikl,isn)))
850     .                                     *ro_dry(ikl,isn)/rovisc
851 
852C +-- Changement de Viscosite si Teneur en Eau liquide
853C +   Change of the Viscosity if liquid Water Content
854C +   ------------------------------------------------
855 
856          OK_Liq     =    max(zero,sign(unun,etaSno(ikl,isn)-epsi))
857          OK_Ang     =    max(zero,sign(unun,vgran6-G1snSV(ikl,isn)))
858     .                *(1-min(1   , abs(istoSV(ikl,isn)-istdSV(1))))
859! #wp     IF (G1snSV(ikl,isn).gt.0..AND.G1snSV(ikl,isn).lt.vsphe4
860! #wp.                             .AND.istoSV(ikl,isn).eq.     0)
861! #wp.    THEN
862! #wp       write(6,*) ikl,isn,' G1,G2,hist,OK_Ang  ',
863! #wp.          G1snSV(ikl,isn), G2snSV(ikl,isn),istoSV(ikl,isn),OK_Ang
864! #wp       stop "Grains anguleux mal d?finis"
865! #wp     END IF
866          OKxLiq     =    max(zero,sign(unun,vtelv1-etaSno(ikl,isn)
867     .                                    /max(epsi,dzsnSV(ikl,isn))))
868     .               *    max(0   ,sign(1   ,istoSV(ikl,isn)
869     .                                      -istdSV(1)      ))
870          ViscSn     =
871     .    ViscSn*(    OK_Liq/(vvisc5+vvisc6*etaSno(ikl,isn)
872     .                            /max(epsi,dzsnSV(ikl,isn)))
873     .           +(1.-OK_Liq)                               )
874     .          *(    OK_Ang*exp(min(ADSdSV,G2snSV(ikl,isn)-vdiam4))
875     .           +(1.-OK_Ang)                                       )
876     .          *(    OKxLiq*        vvisc7
877     .           +(1.-OKxLiq)              )
878 
879 
880C +-- Calcul nouvelle Epaisseur / new Thickness
881C +   -----------------------------------------
882 
883          dzsnew         =
884     .    dzsnSV(ikl,isn)
885     .     *max(vdz3,
886     .         (unun-dt__SV2*max(SnMass(ikl)*cos(slopSV(ikl)),unun)
887     .                     /max(ViscSn                      ,epsi)))
888          rosnew         = ro__SV(ikl,isn) *dzsnSV(ikl,isn)
889     .                            /max(1e-10,dzsnew)
890          rosmax         = 1.   /( (1. -eta_SV(ikl,isn)) /ro_Ice
891     .                           +      eta_SV(ikl,isn)  /ro_Wat)
892          rosnew         = min(rosnew ,rosmax)
893          dzsnew         = dzsnSV(ikl,isn) *ro__SV(ikl,isn)
894     .                            /max(1e-10,rosnew)
895          ro__SV(ikl,isn)= rosnew
896          dzsnSV(ikl,isn)= dzsnew
897          ro_dry(ikl,isn)= ro__SV(ikl,isn)*(1.-eta_SV(ikl,isn))*1.e-3
898C +...    ro_dry: Dry Density (g/cm3)
899C +
900          SnMass(ikl)    = SnMass(ikl)+dSnMas*0.5
901        END DO
902 
903        smb_new = 0.
904        DO isn  = 1, isnoSV(ikl)
905        smb_new = smb_new + dzsnSV(ikl,isn) *ro__SV(ikl,isn)
906        ENDDO
907 
908        isn=1
909        if (dzsnSV(ikl,isn)>0.and.ro__SV(ikl,isn)>0) then
910        dzsnSV(ikl,isn) = dzsnSV(ikl,isn) +0.9999*(smb_old-smb_new)
911     .                  / ro__SV(ikl,isn)
912        endif
913 
914         zn_new = 0
915        DO isn  = 1, isnoSV(ikl)
916         zn_new = zn_new  + dzsnSV(ikl,isn)
917        ENDDO
918        zn4_SV(ikl) = zn4_SV(ikl) + (zn_new - zn_old)
919 
920      END DO
921 
922
923 
924      return
925      end
Note: See TracBrowser for help on using the repository browser.