source: LMDZ6/branches/Amaury_dev/libf/phylmd/inlandsis/sisvat_gsn.f90 @ 5136

Last change on this file since 5136 was 5117, checked in by abarral, 4 months ago

rename modules properly lmdz_*
move some unused files to obsolete/
(lint) uppercase fortran keywords

File size: 42.0 KB
Line 
1
2SUBROUTINE SISVAT_GSn
3
4  ! +------------------------------------------------------------------------+
5  ! | MAR          SISVAT_GSn                                20-09-2003  MAR |
6  ! |   SubRoutine SISVAT_GSn simulates SNOW Metamorphism                    |
7  ! +------------------------------------------------------------------------+
8  ! |                                                                        |
9  ! |   PARAMETERS:  knonv: Total Number of columns =                        |
10  ! |   ^^^^^^^^^^        = Total Number of continental     grid boxes       |
11  ! |                     X       Number of Mosaic Cell per grid box         |
12  ! |                                                                        |
13  ! |   INPUT /  isnoSV   = total Nb of Ice/Snow Layers                      |
14  ! |   OUTPUT:  iiceSV   = total Nb of Ice      Layers                      |
15  ! |   ^^^^^^   istoSV   = 0,...,5 :   Snow     History (see istdSV data)   |
16  ! |                                                                        |
17  ! |   INPUT:   TsisSV   : Soil/Ice Temperatures (layers -nsol,-nsol+1,..,0)|
18  ! |   ^^^^^             & Snow     Temperatures (layers  1,2,...,nsno) [K] |
19  ! |            ro__SV   : Soil/Snow Volumic Mass                   [kg/m3] |
20  ! |            eta_SV   : Soil/Snow Water   Content                [m3/m3] |
21  ! |            slopSV   : Surface Slope                                [-] |
22  ! |            dzsnSV   : Snow Layer        Thickness                  [m] |
23  ! |            dt__SV2   : Time  Step                                   [s] |
24  ! |                                                                        |
25  ! |   INPUT /  G1snSV   : Dendricity (<0) or Sphericity (>0) of Snow Layer |
26  ! |   OUTPUT:  G2snSV   : Sphericity (>0) or Size            of Snow Layer |
27  ! |   ^^^^^^                                                               |
28  ! |                                                                        |
29  ! |   Formalisme adopte pour la Representation des Grains:                 |
30  ! |   Formalism         for the Representation of  Grains:                 |
31  ! |   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                 |
32  ! |                                                                        |
33  ! |             1       - -1                 Neige Fraiche                 |
34  ! |            / \      |                    -------------                 |
35  ! |           /   \     |  Dendricite        decrite  par Dendricite       |
36  ! |          /     \    |  Dendricity                  et Sphericite       |
37  ! |         /       \   |                                                  |
38  ! |        2---------3  -  0                 described by Dendricity       |
39  ! |                                                   and Sphericity       |
40  ! |        |---------|                                                     |
41  ! |        0         1                                                     |
42  ! |        Sphericite                                                      |
43  ! |        Sphericity                                                      |
44  ! |                                                                        |
45  ! |        4---------5  -                                                  |
46  ! |        |         |  |                                                  |
47  ! |        |         |  |  Diametre (1/10eme de mm) (ou Taille)            |
48  ! |        |         |  |  Diameter (1/10th  of mm) (or Size  )            |
49  ! |        |         |  |                                                  |
50  ! |        |         |  |                    Neige non dendritique         |
51  ! |        6---------7  -                    ---------------------         |
52  ! |                                          decrite  par Sphericite       |
53  ! |                                                    et     Taille       |
54  ! |                                          described by Sphericity       |
55  ! |                                                   and       Size       |
56  ! |                                                                        |
57  ! |   Les Variables du Modele:                                             |
58  ! |   Model         Variables:                                             |
59  ! |   ^^^^^^^^^^^^^^^^^^^^^^^^                                             |
60  ! |     Cas Dendritique               Cas non Dendritique                  |
61  ! |                                                                        |
62  ! |     G1snSV        : Dendricite    G1snSV        : Sphericite           |
63  ! |     G2snSV        : Sphericite    G2snSV        : Taille (1/10e mm)    |
64  ! |                                                   Size                 |
65  ! |                                                                        |
66  ! |   Cas Dendritique/ Dendritic Case                                      |
67  ! |   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                                      |
68  ! |   Dendricite(Dendricity) G1snSV                                        |
69  ! |            varie     de -G1_dSV (-99 par defaut / etoile)          a 0 |
70  ! |            division par -G1_dSV pour obtenir des valeurs entre 1  et 0 |
71  ! |            varies  from -G1_dSV (default -99    / fresh snow)     to 0 |
72  ! |            division  by -G1_dSV to obtain values       between 1 and 0 |
73  ! |                                                                        |
74  ! |   Sphericite(Sphericity) G2snSV                                        |
75  ! |            varie     de  0         (cas completement anguleux)         |
76  ! |                       a  G1_dSV (99 par defaut, cas spherique)         |
77  ! |            division par  G1_dSV pour obtenir des valeurs entre 0  et 1 |
78  ! |            varies  from  0      (full faceted)               to G1_dSV |
79  ! |                                                                        |
80  ! |   Cas non Dendritique / non Dendritic Case                             |
81  ! |   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                             |
82  ! |   Sphericite(Sphericity) G1snSV                                        |
83  ! |            varie     de  0         (cas completement anguleux)         |
84  ! |                       a  G1_dSV (99 par defaut, cas spherique)         |
85  ! |            division par  G1_dSV pour obtenir des valeurs entre 0  et 1 |
86  ! |            varies  from  0      (full faceted)               to G1_dSV |
87  ! |                                                                        |
88  ! |   Taille    (Size)       G2snSV                                        |
89  ! |            superieure a  ADSdSV (.4 mm) et ne fait que croitre         |
90  ! |            greater than  ADSdSV (.4 mm) always increases               |
91  ! |                                                                        |
92  ! |   Exemples: Points caracteristiques des Figures ci-dessus              |
93  ! |   ^^^^^^^^^                                                            |
94  ! |                                                                        |
95  ! |               G1snSV    G2snSV     dendricite  sphericite  taille      |
96  ! |                                    dendricity  sphericity  size        |
97  ! |   ------------------------------------------------------------------   |
98  ! |                                                            [1/10 mm]   |
99  ! |     1        -G1_dSV    sph3SN          1           0.5                |
100  ! |     2           0         0             0           0                  |
101  ! |     3           0       G1_dSV          0           1                  |
102  ! |     4           0       ADSdSV                      0       4.         |
103  ! |     5         G1_dSV    ADSdSV-vsphe1               1       3.         |
104  ! |     6           0         --                        0       --         |
105  ! |     7         G1_dSV      --                        1       --         |
106  ! |                                                                        |
107  ! |   par defaut: G1_dSV=99.                                               |
108  ! |                         sph3SN=50.                                     |
109  ! |                         ADSdSV= 4.                                     |
110  ! |                                vsphe1=1.                               |
111  ! |                                                                        |
112  ! |   Methode:                                                             |
113  ! |   ^^^^^^^^                                                             |
114  ! |   1. Evolution Types de Grains selon Lois de Brun et al. (1992):       |
115  ! |      Grain metamorphism according to         Brun et al. (1992):       |
116  ! |      Plusieurs Cas sont a distiguer  /  the different Cases are:       |
117  ! |       1.1 Metamorphose Neige humide  /  wet Snow                       |
118  ! |       1.2 Metamorphose Neige seche   /  dry Snow                       |
119  ! |         1.2.1 Gradient faible        /  low      Temperature Gradient  |
120  ! |         1.2.2 Gradient moyen         /  moderate Temperature Gradient  |
121  ! |         1.2.3 Gradient fort          /  high     Temperature Gradient  |
122  ! |      Dans chaque Cas on separe Neige Dendritique et non Dendritique    |
123  ! |                           le Passage Dendritique -> non Dendritique    |
124  ! |                           se fait lorsque  G1snSV devient > 0          |
125  ! |      the Case of Dentritic or non Dendritic Snow is treated separately |
126  ! |      the Limit   Dentritic -> non Dendritic is reached when G1snSV > 0 |
127  ! |                                                                        |
128  ! |   2. Tassement: Loi de Viscosite adaptee selon le Type de Grains       |
129  ! |      Snow Settling:    Viscosity depends on the   Grain Type           |
130  ! |                                                                        |
131  ! |   3. Update Variables historiques (cas non dendritique seulement)      |
132  ! |      nhSNow defaut                                                     |
133  ! |                0    Cas normal                                         |
134  ! |      istdSV(1) 1    Grains anguleux / faceted cristal                  |
135  ! |      istdSV(2) 2    Grains ayant ete en presence d eau liquide         |
136  ! |                     mais n'ayant pas eu de caractere anguleux    /     |
137  ! |                     liquid water and no faceted cristals before        |
138  ! |      istdSV(3) 3    Grains ayant ete en presence d eau liquide         |
139  ! |                     ayant eu auparavant un caractere anguleux    /     |
140  ! |                     liquid water and    faceted cristals before        |
141  ! |                                                                        |
142  ! |   REFER. : Brun et al.      1989, J. Glaciol 35 pp. 333--342           |
143  ! |   ^^^^^^^^ Brun et al.      1992, J. Glaciol 38 pp.  13-- 22           |
144  ! |            (CROCUS Model, adapted to MAR at CEN by H.Gallee)           |
145  ! |                                                                        |
146  ! |   REFER. : Marbouty, D.     1980, J. Glaciol 26 pp. xxx--xxx           |
147  ! |   ^^^^^^^^ (CROCUS Model, adapted to MAR at CEN by H.Gallee)           |
148  ! |            (for angular shapes)                                        |
149  ! |                                                                        |
150  ! |   Preprocessing  Option: SISVAT IO (not always a standard preprocess.) |
151  ! |   ^^^^^^^^^^^^^^^^^^^^^  ^^^^^^^^^                                     |
152  ! |   FILE                 |      CONTENT                                  |
153  ! |   ~~~~~~~~~~~~~~~~~~~~~+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
154  ! | # SISVAT_GSn.vp        | #vp: OUTPUT/Verification: Snow   Properties   |
155  ! |                        |      unit 47, SubRoutines SISVAT_zSn, _GSn    |
156  ! | # stdout               | #wp: OUTPUT/Verification: Snow   Properties   |
157  ! |                        |      unit  6, SubRoutine  SISVAT_GSn          |
158  ! |                                                                        |
159  ! +------------------------------------------------------------------------+
160
161
162
163
164  ! +--Global Variables
165  ! +  ================
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
179  ! +--INPUT/OUTPUT
180  ! +  ------------
181
182
183  ! +--OUTPUT
184  ! +  ------
185
186  INTEGER :: dt__SV2
187
188
189  ! +--Local  Variables
190  ! +  ================
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
287  ! +--Snow Properties: IO
288  ! +  ~~~~~~~~~~~~~~~~~~~
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
295  ! +--DATA
296  ! +  ====
297
298  data       vector/.TRUE./               ! Vectorization  Switch
299  data       vdent1/ 0.5e8/               ! Wet Snow Metamorphism
300  !XF                      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
345  ! +--DATA (Coefficient Fonction fort Gradient Marbouty)
346  ! +  --------------------------------------------------
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
380  ! +-- 1. Metamorphoses dans les Strates
381  ! +      Metamorphism
382  ! +      ==============================
383
384  dt__SV2= dt__SV
385  frac_j = dt__SV2 / 86400.                        ! Time Step [Day]
386
387  zn4_SV = 0
388
389
390  ! +-- 1.1 Initialisation: teneur en eau liquide et gradient de temperature
391  ! +   ------------------  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
404  !!$OMP PARALLEL DO default(firstprivate)
405  !!$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))
414  ! +...    Factor 1.d-2 for Conversion K/m --> K/cm
415
416
417  ! +-- 1.2 Metamorphose humide
418  ! +       Wet Snow Metamorphism
419  ! +       ---------------------
420
421      Wet_OK = max(zero,sign(unun,eta_SV(ikl,isn)-epsi))
422
423
424  ! +--     Vitesse de diminution de la dendricite
425  ! +       Rate of the dendricity decrease
426  ! +       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
427      sWater=1.d-1*ro__SV(ikl,isn)*eta_SV(ikl,isn) &
428            /max(epsi,ro_dry(ikl,isn))
429  ! +...    sWater:Water Content [%]
430  ! +              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
435  ! +-- 1.2.1 Cas dendritique/dendritic Case
436  ! +         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
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
455  ! +-- 1.2.2 Cas non dendritique non completement spherique
456  ! +         Evolution de la Sphericite seulement.
457  ! +         Non dendritic and not completely spheric Case
458  ! +         Evolution of    Sphericity only (not size)
459  ! +         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
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
469  ! +-- 1.2.3 Cas non dendritique et spherique / non dendritic and spheric
470  ! +         Evolution de la Taille seulement / Evolution of Size only
471  ! +         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
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
479  ! +-- 1.3 Metamorposes seches / Dry Metamorphism
480  ! +       --------------------------------------
481
482
483  ! +-- 1.3.1 Calcul Metamorphose faible/low Gradient (0.00-0.05 deg/cm)
484  ! +         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
485      OKlowT=max(zero, & !
486            sign(unun, vgrat1 & !
487            -dTsndz         ))     !
488
489      facVap=exp(vvap1/TsisSV(ikl,isn))
490
491  ! +-- 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
511  ! +-- 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
532  ! +-- 1.3.2 Calcul Metamorphose Gradient Moyen/Moderate (0.05-0.15)
533  ! +         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
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
543  ! +-- 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
563  ! +-- 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
569  ! +-- 1.3.3 Calcul Metamorphose fort / high Gradient
570  ! +   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
571      facVap=vdent1*exp(vvap1/TsisSV(ikl,isn)) &
572            *   (1.e2 *dTsndz)**vvap2
573
574  ! +-- 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
584  ! +                                                   ! 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
594  ! +-- 1.3.3.2 Cas non dendritique non completement anguleux.
595  ! +           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
605  ! +-- 1.3.3.3 Cas non dendritique et anguleux
606  ! +           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
625  ! +-- Influence de la Temperature /Temperature Influence
626  ! +   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
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
636  ! +-- Influence de la Masse Volumique /Density Influence
637  ! +   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
638            * ro1_OK &
639            *(   ro2_OK*(1. - (ro_dry(ikl,isn)-vrang2) & !
640            /(vrang1-vrang2)) & !
641            +1.-ro2_OK                                ) & !
642
643  ! +-- Influence du Gradient de Temperature /Temperature Gradient Influence
644  ! +   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
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
662  ! +--New Properties
663  ! +  --------------
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
680  !XF
681  IF(G1snSV(ikl,isn)<0.1) &
682        G2_hds = G2snSV(ikl,isn) + 1.d1 *AngSno*vfi     *frac_j
683  !XF
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
698  ! +--Snow Properties: IO Set Up
699  ! +  ~~~~~~~~~~~~~~~~~~~~~~~~~~
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
738  ! +--Snow Properties: IO
739  ! +  ~~~~~~~~~~~~~~~~~~~
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
762  !!$OMP END PARALLEL DO
763
764  ! +-- 2. Mise a Jour Variables Historiques (Cas non dendritique)
765  ! +      Update of the historical Variables
766  ! +      =======================================================
767
768  IF (vector)                                                   THEN
769  !XF
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
801  ! +-- 2. Mise a Jour Variables Historiques (Cas non dendritique)
802  ! +      Update of the historical Variables
803  ! +      =======================================================
804
805    DO ikl=1,knonv
806    DO isn=iiceSV(ikl),isnoSV(ikl)
807      IF  (G1snSV(ikl,isn)>=0.)                               THEN
808        IF(G1snSV(ikl,isn)<vsphe4.AND.istoSV(ikl,isn)==0)  THEN
809               istoSV(ikl,isn)=istdSV(1)
810        ELSEIF(G1_dSV-G1snSV(ikl,isn)         <vsphe4.AND. &
811              etaSno(ikl,isn)/dzsnSV(ikl,isn)>vtelv1)       THEN
812          IF  (istoSV(ikl,isn)==0) &
813                istoSV(ikl,isn)=   istdSV(2)
814          IF  (istoSV(ikl,isn)==istdSV(1)) &
815                istoSV(ikl,isn)=   istdSV(3)
816        ELSEIF(TsisSV(ikl,isn)<TfSnow)                       THEN
817          IF  (istoSV(ikl,isn)==istdSV(2)) &
818                istoSV(ikl,isn)=   istdSV(4)
819          IF  (istoSV(ikl,isn)==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
828  ! +-- 3. Tassement mecanique /mechanical Settlement
829  ! +      ==========================================
830
831    DO ikl=1,knonv
832      SnMass(ikl) = 0.
833    END DO
834  !XF
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
852  ! +-- Changement de Viscosite si Teneur en Eau liquide
853  ! +   Change of the Viscosity if liquid Water Content
854  ! +   ------------------------------------------------
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
880  ! +-- Calcul nouvelle Epaisseur / new Thickness
881  ! +   -----------------------------------------
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
898  ! +...    ro_dry: Dry Density (g/cm3)
899  ! +
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
925END SUBROUTINE sisvat_gsn
Note: See TracBrowser for help on using the repository browser.