
subroutine SISVAT_GSn

  ! +------------------------------------------------------------------------+
  ! | MAR          SISVAT_GSn                                20-09-2003  MAR |
  ! |   SubRoutine SISVAT_GSn simulates SNOW Metamorphism                    |
  ! +------------------------------------------------------------------------+
  ! |                                                                        |
  ! |   PARAMETERS:  knonv: Total Number of columns =                        |
  ! |   ^^^^^^^^^^        = Total Number of continental     grid boxes       |
  ! |                     X       Number of Mosaic Cell per grid box         |
  ! |                                                                        |
  ! |   INPUT /  isnoSV   = total Nb of Ice/Snow Layers                      |
  ! |   OUTPUT:  iiceSV   = total Nb of Ice      Layers                      |
  ! |   ^^^^^^   istoSV   = 0,...,5 :   Snow     History (see istdSV data)   |
  ! |                                                                        |
  ! |   INPUT:   TsisSV   : Soil/Ice Temperatures (layers -nsol,-nsol+1,..,0)|
  ! |   ^^^^^             & Snow     Temperatures (layers  1,2,...,nsno) [K] |
  ! |            ro__SV   : Soil/Snow Volumic Mass                   [kg/m3] |
  ! |            eta_SV   : Soil/Snow Water   Content                [m3/m3] |
  ! |            slopSV   : Surface Slope                                [-] |
  ! |            dzsnSV   : Snow Layer        Thickness                  [m] |
  ! |            dt__SV2   : Time  Step                                   [s] |
  ! |                                                                        |
  ! |   INPUT /  G1snSV   : Dendricity (<0) or Sphericity (>0) of Snow Layer |
  ! |   OUTPUT:  G2snSV   : Sphericity (>0) or Size            of Snow Layer |
  ! |   ^^^^^^                                                               |
  ! |                                                                        |
  ! |   Formalisme adopte pour la Representation des Grains:                 |
  ! |   Formalism         for the Representation of  Grains:                 |
  ! |   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                 |
  ! |                                                                        |
  ! |             1       - -1                 Neige Fraiche                 |
  ! |            / \      |                    -------------                 |
  ! |           /   \     |  Dendricite        decrite  par Dendricite       |
  ! |          /     \    |  Dendricity                  et Sphericite       |
  ! |         /       \   |                                                  |
  ! |        2---------3  -  0                 described by Dendricity       |
  ! |                                                   and Sphericity       |
  ! |        |---------|                                                     |
  ! |        0         1                                                     |
  ! |        Sphericite                                                      |
  ! |        Sphericity                                                      |
  ! |                                                                        |
  ! |        4---------5  -                                                  |
  ! |        |         |  |                                                  |
  ! |        |         |  |  Diametre (1/10eme de mm) (ou Taille)            |
  ! |        |         |  |  Diameter (1/10th  of mm) (or Size  )            |
  ! |        |         |  |                                                  |
  ! |        |         |  |                    Neige non dendritique         |
  ! |        6---------7  -                    ---------------------         |
  ! |                                          decrite  par Sphericite       |
  ! |                                                    et     Taille       |
  ! |                                          described by Sphericity       |
  ! |                                                   and       Size       |
  ! |                                                                        |
  ! |   Les Variables du Modele:                                             |
  ! |   Model         Variables:                                             |
  ! |   ^^^^^^^^^^^^^^^^^^^^^^^^                                             |
  ! |     Cas Dendritique               Cas non Dendritique                  |
  ! |                                                                        |
  ! |     G1snSV        : Dendricite    G1snSV        : Sphericite           |
  ! |     G2snSV        : Sphericite    G2snSV        : Taille (1/10e mm)    |
  ! |                                                   Size                 |
  ! |                                                                        |
  ! |   Cas Dendritique/ Dendritic Case                                      |
  ! |   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                                      |
  ! |   Dendricite(Dendricity) G1snSV                                        |
  ! |            varie     de -G1_dSV (-99 par defaut / etoile)          a 0 |
  ! |            division par -G1_dSV pour obtenir des valeurs entre 1  et 0 |
  ! |            varies  from -G1_dSV (default -99    / fresh snow)     to 0 |
  ! |            division  by -G1_dSV to obtain values       between 1 and 0 |
  ! |                                                                        |
  ! |   Sphericite(Sphericity) G2snSV                                        |
  ! |            varie     de  0         (cas completement anguleux)         |
  ! |                       a  G1_dSV (99 par defaut, cas spherique)         |
  ! |            division par  G1_dSV pour obtenir des valeurs entre 0  et 1 |
  ! |            varies  from  0      (full faceted)               to G1_dSV |
  ! |                                                                        |
  ! |   Cas non Dendritique / non Dendritic Case                             |
  ! |   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^                             |
  ! |   Sphericite(Sphericity) G1snSV                                        |
  ! |            varie     de  0         (cas completement anguleux)         |
  ! |                       a  G1_dSV (99 par defaut, cas spherique)         |
  ! |            division par  G1_dSV pour obtenir des valeurs entre 0  et 1 |
  ! |            varies  from  0      (full faceted)               to G1_dSV |
  ! |                                                                        |
  ! |   Taille    (Size)       G2snSV                                        |
  ! |            superieure a  ADSdSV (.4 mm) et ne fait que croitre         |
  ! |            greater than  ADSdSV (.4 mm) always increases               |
  ! |                                                                        |
  ! |   Exemples: Points caracteristiques des Figures ci-dessus              |
  ! |   ^^^^^^^^^                                                            |
  ! |                                                                        |
  ! |               G1snSV    G2snSV     dendricite  sphericite  taille      |
  ! |                                    dendricity  sphericity  size        |
  ! |   ------------------------------------------------------------------   |
  ! |                                                            [1/10 mm]   |
  ! |     1        -G1_dSV    sph3SN          1           0.5                |
  ! |     2           0         0             0           0                  |
  ! |     3           0       G1_dSV          0           1                  |
  ! |     4           0       ADSdSV                      0       4.         |
  ! |     5         G1_dSV    ADSdSV-vsphe1               1       3.         |
  ! |     6           0         --                        0       --         |
  ! |     7         G1_dSV      --                        1       --         |
  ! |                                                                        |
  ! |   par defaut: G1_dSV=99.                                               |
  ! |                         sph3SN=50.                                     |
  ! |                         ADSdSV= 4.                                     |
  ! |                                vsphe1=1.                               |
  ! |                                                                        |
  ! |   Methode:                                                             |
  ! |   ^^^^^^^^                                                             |
  ! |   1. Evolution Types de Grains selon Lois de Brun et al. (1992):       |
  ! |      Grain metamorphism according to         Brun et al. (1992):       |
  ! |      Plusieurs Cas sont a distiguer  /  the different Cases are:       |
  ! |       1.1 Metamorphose Neige humide  /  wet Snow                       |
  ! |       1.2 Metamorphose Neige seche   /  dry Snow                       |
  ! |         1.2.1 Gradient faible        /  low      Temperature Gradient  |
  ! |         1.2.2 Gradient moyen         /  moderate Temperature Gradient  |
  ! |         1.2.3 Gradient fort          /  high     Temperature Gradient  |
  ! |      Dans chaque Cas on separe Neige Dendritique et non Dendritique    |
  ! |                           le Passage Dendritique -> non Dendritique    |
  ! |                           se fait lorsque  G1snSV devient > 0          |
  ! |      the Case of Dentritic or non Dendritic Snow is treated separately |
  ! |      the Limit   Dentritic -> non Dendritic is reached when G1snSV > 0 |
  ! |                                                                        |
  ! |   2. Tassement: Loi de Viscosite adaptee selon le Type de Grains       |
  ! |      Snow Settling:    Viscosity depends on the   Grain Type           |
  ! |                                                                        |
  ! |   3. Update Variables historiques (cas non dendritique seulement)      |
  ! |      nhSNow defaut                                                     |
  ! |                0    Cas normal                                         |
  ! |      istdSV(1) 1    Grains anguleux / faceted cristal                  |
  ! |      istdSV(2) 2    Grains ayant ete en presence d eau liquide         |
  ! |                     mais n'ayant pas eu de caractere anguleux    /     |
  ! |                     liquid water and no faceted cristals before        |
  ! |      istdSV(3) 3    Grains ayant ete en presence d eau liquide         |
  ! |                     ayant eu auparavant un caractere anguleux    /     |
  ! |                     liquid water and    faceted cristals before        |
  ! |                                                                        |
  ! |   REFER. : Brun et al.      1989, J. Glaciol 35 pp. 333--342           |
  ! |   ^^^^^^^^ Brun et al.      1992, J. Glaciol 38 pp.  13-- 22           |
  ! |            (CROCUS Model, adapted to MAR at CEN by H.Gallee)           |
  ! |                                                                        |
  ! |   REFER. : Marbouty, D.     1980, J. Glaciol 26 pp. xxx--xxx           |
  ! |   ^^^^^^^^ (CROCUS Model, adapted to MAR at CEN by H.Gallee)           |
  ! |            (for angular shapes)                                        |
  ! |                                                                        |
  ! |   Preprocessing  Option: SISVAT IO (not always a standard preprocess.) |
  ! |   ^^^^^^^^^^^^^^^^^^^^^  ^^^^^^^^^                                     |
  ! |   FILE                 |      CONTENT                                  |
  ! |   ~~~~~~~~~~~~~~~~~~~~~+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
  ! | # SISVAT_GSn.vp        | #vp: OUTPUT/Verification: Snow   Properties   |
  ! |                        |      unit 47, SubRoutines SISVAT_zSn, _GSn    |
  ! | # stdout               | #wp: OUTPUT/Verification: Snow   Properties   |
  ! |                        |      unit  6, SubRoutine  SISVAT_GSn          |
  ! |                                                                        |
  ! +------------------------------------------------------------------------+




  ! +--Global Variables
  ! +  ================

  use VARphy
  use VAR_SV
  use VARdSV
  use VAR0SV
  use VARxSV
  use VARtSV


  IMPLICIT NONE



  ! +--INPUT/OUTPUT
  ! +  ------------


  ! +--OUTPUT
  ! +  ------

  integer :: dt__SV2


  ! +--Local  Variables
  ! +  ================

  logical :: vector                        !
  integer :: ikl                           !
  integer :: isn   ,isnp                   !
  integer :: istoOK                        !
  real :: G1_bak,G2_bak                 ! Old Values of G1, G2
  real :: ro_dry(knonv,      nsno)      ! Dry Density            [g/cm3]
  real :: etaSno(knonv,      nsno)      ! Liquid Water Content   [g/cm2]
  real :: SnMass(knonv)                 ! Snow   Mass            [kg/m2]
  real :: dTsndz                        ! Temperature Gradient
  real :: sWater                        !        Water Content       [%]
  real :: exp1Wa                        !
  real :: dDENDR                        ! Dendricity Increment
  real :: DENDRn                        ! Normalized Dendricity
  real :: SPHERn                        ! Normalized Sphericity
  real :: Wet_OK                        ! Wet Metamorphism Switch
  real :: OK__DE                        !
  real :: OK__wd                        ! New G*, from wet Dendritic
  real :: G1__wd                        ! New G1, from wet Dendritic
  real :: G2__wd                        ! New G2, from wet Dendritic
  real :: OKlowT                        !
  real :: facVap                        !
  real :: OK_ldd                        !
  real :: G1_ldd                        !
  real :: G2_ldd                        !
  real :: DiamGx                        !
  real :: DiamOK                        !
  real :: No_Big                        !
  real :: dSPHER                        !
  real :: SPHER0                        !
  real :: SPHbig                        !
  real :: G1_lds                        !
  real :: OK_mdT                        !
  real :: OKmidT                        !
  real :: OKhigT                        !
  real :: OK_mdd                        !
  real :: G1_mdd                        !
  real :: G2_mdd                        !
  real :: G1_mds                        !
  real :: OK_hdd                        !
  real :: G1_hdd                        !
  real :: G2_hdd                        !
  real :: OK_hds                        !
  real :: G1_hds                        !
  real :: T1__OK,T2__OK                 !
  real :: T3_xOK,T3__OK,T3_nOK          !
  real :: ro1_OK,ro2_OK                 !
  real :: dT1_OK,dT2_OK,dT3xOK,dT3_OK   !
  real :: dT4xOK,dT4_OK,dT4nOK,AngSno   !
  real :: G2_hds,SphrOK,HISupd          !
  real :: H1a_OK,H1b_OK,H1__OK          !
  real :: H23aOK,H23bOK,H23_OK          !
  real :: H2__OK,H3__OK                 !
  real :: H45_OK,H4__OK,H5__OK          !
  real :: ViscSn,OK_Liq,OK_Ang,OKxLiq   !
  real :: dSnMas,dzsnew,rosnew,rosmax,smb_old,smb_new
  real :: zn_old,zn_new

  real :: epsi5                         ! Alpha ev67 single precision	
  real :: vdiam1                        ! Small Grains Min.Diam.[.0001m]
  real :: vdiam2                        ! Spher.Variat.Max Diam.    [mm]
  real :: vdiam3                        ! Min.Diam.|Limit Spher.    [mm]
  real :: vdiam4                        ! Min.Diam.|Viscosity Change
  real :: vsphe1                        ! Max Sphericity
  real :: vsphe2                        ! Low T Metamorphism  Coeff.
  real :: vsphe3                        ! Max.Sphericity (history=1)
  real :: vsphe4                        ! Min.Sphericity=>history=1
  real :: vtang1,vtang2,vtang3,vtang4   ! Temperature Contribution
  real :: vtang5,vtang6,vtang7,vtang8   !
  real :: vtang9,vtanga,vtangb,vtangc   !
  real :: vrang1,vrang2                 ! Density     Contribution
  real :: vgang1,vgang2,vgang3,vgang4   ! Grad(T)     Contribution
  real :: vgang5,vgang6,vgang7,vgang8   !
  real :: vgang9,vganga,vgangb,vgangc   !
  real :: vgran6                        ! Max.Sphericity for Settling
  real :: vtelv1                        ! Threshold | history = 2, 3
  real :: vvap1                        ! Vapor Pressure Coefficient
  real :: vvap2                        ! Vapor Pressure Exponent
  real :: vgrat1                        ! Boundary weak/mid   grad(T)
  real :: vgrat2                        ! Boundary mid/strong grad(T)
  real :: vfi                        ! PHI,         strong grad(T)
  real :: vvisc1,vvisc2,vvisc3,vvisc4   ! Viscosity Coefficients
  real :: vvisc5,vvisc6,vvisc7          ! id., wet Snow
  real :: rovisc                        ! Wet Snow Density  Influence
  real :: vdz3                        ! Maximum Layer Densification
  real :: OK__ws                        ! New G2
  real :: G1__ws                        ! New G1, from wet Spheric
  real :: G2__ws                        ! New G2, from wet Spheric
  real :: husi_0,husi_1,husi_2,husi_3   ! Constants for New G2
  real :: vtail1,vtail2                 ! Constants for New G2
  real :: frac_j                        ! Time Step            [Day]

  real :: vdent1                       ! Wet Snow Metamorphism
  integer :: nvdent1                       ! (Coefficients for
  integer :: nvdent2                       !           Dendricity)

  ! +--Snow Properties: IO
  ! +  ~~~~~~~~~~~~~~~~~~~
  ! #vp real      G_curr(18),Gcases(18)
  ! #vp common   /GSnLOC/    Gcases
  ! #wp real                 D__MAX
  ! #wp common   /GSnMAX/    D__MAX


  ! +--DATA
  ! +  ====

  data       vector/.true./               ! Vectorization  Switch
  data       vdent1/ 0.5e8/               ! Wet Snow Metamorphism
  !XF                      tuned for Greenland (2.e8=old value)
  data      nvdent1/ 3   /                ! (Coefficients for
  data      nvdent2/16   /                !           Dendricity)

  data       husi_0 /20.      /           !   10  * 2
  data       husi_1 / 0.23873 /           ! (3/4) /pi
  data       husi_2 / 4.18880 /           ! (4/3) *pi
  data       husi_3 / 0.33333 /           !  1/3
  data       vtail1 / 1.28e-08/           !  Wet Metamorphism
  data       vtail2 / 4.22e-10/           ! (NON Dendritic / Spheric)

  data       epsi5  / 1.0e-5  /           !

  data       vdiam1 / 4.0     /           ! Small Grains Min.Diameter

  data       vdiam2 / 0.5     /           ! Spher.Variat.Max Diam.[mm]
  data       vdiam3 / 3.0     /           ! Min.Diam.|Limit Spher.[mm]
  data       vdiam4 / 2.0     /           ! Min.Diam.|Viscosity Change

  data       vsphe1 / 1.0     /           ! Max Sphericity
  data       vsphe2 / 1.0e9   /           ! Low T Metamorphism  Coeff.
  data       vsphe3 / 0.5     /           ! Max.Sphericity (history=1)
  data       vsphe4 / 0.1     /           ! Min.Sphericity=>history=1

  data       vgran6 / 51.     /           ! Max.Sphericity for Settling
  data       vtelv1 / 5.e-1   /           ! Threshold | history = 2, 3

  data        vvap1 /-6.e3    /           ! Vapor Pressure Coefficient
  data        vvap2 / 0.4     /           ! Vapor Pressure Exponent

  data       vgrat1 /0.05     /           ! Boundary weak/mid   grad(T)
  data       vgrat2 /0.15     /           ! Boundary mid/strong grad(T)
  data          vfi /0.09     /           ! PHI,         strong grad(T)

  data       vvisc1 / 0.70    /           ! Viscosity Coefficients
  data       vvisc2 / 1.11e5  /           !
  data       vvisc3 /23.00    /           !
  data       vvisc4 / 0.10    /           !
  data       vvisc5 / 1.00    /           ! id., wet Snow
  data       vvisc6 / 2.00    /           !
  data       vvisc7 /10.00    /           !
  data       rovisc / 0.25    /           ! Wet Snow Density  Influence
  data         vdz3 / 0.30    /           ! Maximum Layer Densification


  ! +--DATA (Coefficient Fonction fort Gradient Marbouty)
  ! +  --------------------------------------------------

  data       vtang1 /40.0/                ! Temperature Contribution
  data       vtang2 / 6.0/                !
  data       vtang3 /22.0/                !
  data       vtang4 / 0.7/                !
  data       vtang5 / 0.3/                !
  data       vtang6 / 6.0/                !
  data       vtang7 / 1.0/                !
  data       vtang8 / 0.8/                !
  data       vtang9 /16.0/                !
  data       vtanga / 0.2/                !
  data       vtangb / 0.2/                !
  data       vtangc /18.0/                !

  data       vrang1 / 0.40/               ! Density     Contribution
  data       vrang2 / 0.15/               !

  data       vgang1 / 0.70/               ! Grad(T)     Contribution
  data       vgang2 / 0.25/               !
  data       vgang3 / 0.40/               !
  data       vgang4 / 0.50/               !
  data       vgang5 / 0.10/               !
  data       vgang6 / 0.15/               !
  data       vgang7 / 0.10/               !
  data       vgang8 / 0.55/               !
  data       vgang9 / 0.65/               !
  data       vganga / 0.20/               !
  data       vgangb / 0.85/               !
  data       vgangc / 0.15/               !

  ! #wp data       D__MAX / 4.00/               !


  ! +-- 1. Metamorphoses dans les Strates
  ! +      Metamorphism
  ! +      ==============================

  dt__SV2= dt__SV
  frac_j = dt__SV2 / 86400.                        ! Time Step [Day]

  zn4_SV = 0


  ! +-- 1.1 Initialisation: teneur en eau liquide et gradient de temperature
  ! +   ------------------  liquid water content and temperature gradient

   DO ikl=1,knonv
    DO   isn=1,isnoSV(ikl)

      ro_dry(ikl,isn) = 1.e-3 *ro__SV(ikl,isn) & ! Dry Density
            *(1.    -eta_SV(ikl,isn))   !         [g/cm3]
      etaSno(ikl,isn) = 1.e-1 *dzsnSV(ikl,isn) & ! Liquid Water
            *        ro__SV(ikl,isn) & ! Content [g/cm2]
            *        eta_SV(ikl,isn)    !
    END DO
  END DO

  !!$OMP PARALLEL DO default(firstprivate)
  !!$OMP.shared (/xSISVAT_I/,/xSISVAT_R/,/SoR0SV/,/SoI0SV/,/Sn_dSV/)
   DO ikl=1,knonv
    DO   isn=1,isnoSV(ikl)
      isnp   = min(isn+1,isnoSV(ikl))

      dTsndz = abs( (TsisSV(ikl,isnp)-TsisSV(ikl,isn-1)) *2.e-2 &
            /max(((dzsnSV(ikl,isnp)+dzsnSV(ikl,isn)  ) &
            *(           isnp -           isn) &
            +(dzsnSV(ikl,isn )+dzsnSV(ikl,isn-1))),epsi))
  ! +...    Factor 1.d-2 for Conversion K/m --> K/cm


  ! +-- 1.2 Metamorphose humide
  ! +       Wet Snow Metamorphism
  ! +       ---------------------

      Wet_OK = max(zero,sign(unun,eta_SV(ikl,isn)-epsi))


  ! +--     Vitesse de diminution de la dendricite
  ! +       Rate of the dendricity decrease
  ! +       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      sWater=1.d-1*ro__SV(ikl,isn)*eta_SV(ikl,isn) &
            /max(epsi,ro_dry(ikl,isn))
  ! +...    sWater:Water Content [%]
  ! +              1.d-1= 1.d2(1->%) * 1.d-3(ro__SV*eta_SV:kg/m3->g/cm3)

      exp1Wa=   sWater**nvdent1
      dDENDR=max(exp1Wa/nvdent2,vdent1*exp(vvap1/TfSnow))

  ! +-- 1.2.1 Cas dendritique/dendritic Case
  ! +         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      OK__wd=max(zero, & !
            sign(unun,-G1snSV(ikl,isn) & !
            -epsi           ))     !

      DENDRn=-G1snSV(ikl,isn)/G1_dSV  ! Normalized Dendricity (+)
      SPHERn= G2snSV(ikl,isn)/G1_dSV  ! Normalized Sphericity
      DENDRn= DENDRn -dDENDR *frac_j  ! New        Dendricity (+)
      SPHERn= SPHERn +dDENDR *frac_j  ! New        Sphericity

      OK__DE=max(zero, & ! IF 1.,
            sign(unun, DENDRn & ! NO change
            -epsi           ))     ! Dendr. -> Spheric

      G1__wd=OK__DE *    (      -DENDRn*G1_dSV) & ! Dendritic
            +(1.-OK__DE)* min(G1_dSV,SPHERn*G1_dSV)   ! Dendr. -> Spheric
      G2__wd=OK__DE * min(G1_dSV,SPHERn*G1_dSV) & ! Spheric
            +(1.-OK__DE)*(ADSdSV-min(SPHERn,vsphe1))  ! Spher. -> Size

  ! +-- 1.2.2 Cas non dendritique non completement spherique
  ! +         Evolution de la Sphericite seulement.
  ! +         Non dendritic and not completely spheric Case
  ! +         Evolution of    Sphericity only (not size)
  ! +         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      OK__ws=max(zero, & !
            sign(unun, G1_dSV & !
            -epsi5 & !
            -G1snSV(ikl,isn)))     !

      SPHERn= G1snSV(ikl,isn)/G1_dSV
      SPHERn= SPHERn +dDENDR *frac_j
      G1__ws=         min(G1_dSV,SPHERn*G1_dSV)

  ! +-- 1.2.3 Cas non dendritique et spherique / non dendritic and spheric
  ! +         Evolution de la Taille seulement / Evolution of Size only
  ! +         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      G2__ws =  husi_0 &
            *( husi_1 &
            *(husi_2 *( G2snSV(ikl,isn)/husi_0)**3 &
            +(vtail1 +vtail2 *exp1Wa    )*dt__SV2)) &
            ** husi_3


  ! +-- 1.3 Metamorposes seches / Dry Metamorphism
  ! +       --------------------------------------


  ! +-- 1.3.1 Calcul Metamorphose faible/low Gradient (0.00-0.05 deg/cm)
  ! +         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      OKlowT=max(zero, & !
            sign(unun, vgrat1 & !
            -dTsndz         ))     !

      facVap=exp(vvap1/TsisSV(ikl,isn))

  ! +-- 1.3.1.1 Cas dendritique / dendritic Case

      OK_ldd=max(zero, & !
            sign(unun,-G1snSV(ikl,isn) & !
            -epsi           ))     !

      DENDRn=-G1snSV(ikl,isn)     /G1_dSV
      SPHERn= G2snSV(ikl,isn)     /G1_dSV
      DENDRn= DENDRn-vdent1*facVap*frac_j
      SPHERn= SPHERn+vsphe2*facVap*frac_j

      OK__DE=max(zero, & ! IF 1.,
            sign(unun, DENDRn & ! NO change
            -epsi           ))     ! Dendr. -> Spheric

      G1_ldd= OK__DE *    (      -DENDRn*G1_dSV) & ! Dendritic
            +(1.-OK__DE)* min(G1_dSV,SPHERn*G1_dSV)  ! Dendr. -> Spheric
      G2_ldd= OK__DE * min(G1_dSV,SPHERn*G1_dSV) & ! Spheric
            +(1.-OK__DE)*(ADSdSV-min(SPHERn,vsphe1)) ! Spher. -> Size

  ! +-- 1.3.1.2 Cas non dendritique / non dendritic Case

      SPHERn=G1snSV(ikl,isn)/G1_dSV
      DiamGx=G2snSV(ikl,isn)*0.1

      istoOK=min( abs(istoSV(ikl,isn)- &
            istdSV(1)      ),1)         ! zero if istoSV = 1
      DiamOK=max(zero,  sign(unun,vdiam2-DiamGx))
      No_Big=    istoOK+DiamOK
      No_Big=min(No_Big,unun)

      dSPHER=           vsphe2*facVap*frac_j      !
      SPHER0=    SPHERn+dSPHER                    ! small grains
      SPHbig=    SPHERn+dSPHER & ! big   grains
            *exp(min(zero,vdiam3-G2snSV(ikl,isn)))  ! (history = 2 or 3)
      SPHbig=       min(vsphe3,SPHbig)            ! limited sphericity
      SPHERn= No_Big *  SPHER0 &
            + (1.-No_Big)*  SPHbig

      G1_lds=       min(G1_dSV,SPHERn*G1_dSV)

  ! +-- 1.3.2 Calcul Metamorphose Gradient Moyen/Moderate (0.05-0.15)
  ! +         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      OK_mdT=max(zero, & !
            sign(unun, vgrat2 & !
            -dTsndz))              !
      OKmidT=               OK_mdT  *(1.-OKlowT)  !
      OKhigT=          (1. -OK_mdT) *(1.-OKlowT)  !

      facVap=vdent1*exp(vvap1/TsisSV(ikl,isn)) &
            *   (1.e2 *dTsndz)**vvap2

  ! +-- 1.3.2.1 cas dendritique / dendritic case.

      OK_mdd=max(zero, & !
            sign(unun,-G1snSV(ikl,isn) & !
            -epsi           ))     !

      DENDRn=-G1snSV(ikl,isn)/G1_dSV
      SPHERn= G2snSV(ikl,isn)/G1_dSV
      DENDRn= DENDRn - facVap*frac_j
      SPHERn= SPHERn - facVap*frac_j

      OK__DE=max(zero, & ! IF 1.,
            sign(unun, DENDRn & ! NO change
            -epsi           ))     ! Dendr. -> Spheric

      G1_mdd= OK__DE *    (      -DENDRn*G1_dSV) & ! Dendritic
            +(1.-OK__DE)* max(zero  ,SPHERn*G1_dSV)  ! Dendr. -> Spheric
      G2_mdd= OK__DE * max(zero  ,SPHERn*G1_dSV) & ! Spheric
            +(1.-OK__DE)*(ADSdSV-max(SPHERn,zero  )) ! Spher. -> Size

  ! +-- 1.3.2.2 Cas non dendritique / non dendritic Case

      SPHERn=G1snSV(ikl,isn)/G1_dSV
      SPHERn=         SPHERn-facVap*frac_j
      G1_mds=max(zero,SPHERn*G1_dSV)

  ! +-- 1.3.3 Calcul Metamorphose fort / high Gradient
  ! +   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      facVap=vdent1*exp(vvap1/TsisSV(ikl,isn)) &
            *   (1.e2 *dTsndz)**vvap2

  ! +-- 1.3.3.1 Cas dendritique / dendritic Case

      OK_hdd=max(zero, & !
            sign(unun,-G1snSV(ikl,isn) & !
            -epsi           ))     !

      DENDRn=-G1snSV(ikl,isn)/G1_dSV              !
      SPHERn= G2snSV(ikl,isn)/G1_dSV              !
      DENDRn= DENDRn - facVap*frac_j              !
      SPHERn= SPHERn - facVap*frac_j              ! Non dendritic
  ! +                                                   ! and angular
      OK__DE=max(zero, & ! IF 1.,
            sign(unun, DENDRn & ! NO change
            -epsi  ))              ! Dendr. -> Spheric

      G1_hdd= OK__DE *    (      -DENDRn*G1_dSV) & ! Dendritic
            +(1.-OK__DE)* max(zero  ,SPHERn*G1_dSV)  ! Dendr. -> Spheric
      G2_hdd= OK__DE * max(zero  ,SPHERn*G1_dSV) & ! Spheric
            +(1.-OK__DE)*(ADSdSV-max(SPHERn,zero  )) ! Spher. -> Size

  ! +-- 1.3.3.2 Cas non dendritique non completement anguleux.
  ! +           non dendritic and spericity gt. 0

      OK_hds=max(zero, & !
            sign(unun, G1snSV(ikl,isn) & !
            -epsi           ))     !

      SPHERn= G1snSV(ikl,isn)/G1_dSV
      SPHERn= SPHERn - facVap*frac_j
      G1_hds= max(zero,SPHERn*G1_dSV)

  ! +-- 1.3.3.3 Cas non dendritique et anguleux
  ! +           dendritic and spericity = 0.

      T1__OK = max(zero,sign(unun,TsisSV(ikl,isn)-TfSnow+vtang1))
      T2__OK = max(zero,sign(unun,TsisSV(ikl,isn)-TfSnow+vtang2))
      T3_xOK = max(zero,sign(unun,TsisSV(ikl,isn)-TfSnow+vtang3))
      T3__OK =                    T3_xOK  * (1. - T2__OK)
      T3_nOK =              (1. - T3_xOK) * (1. - T2__OK)
      ro1_OK = max(zero,sign(unun,vrang1-ro_dry(ikl,isn)))
      ro2_OK = max(zero,sign(unun,ro_dry(ikl,isn)-vrang2))
      dT1_OK = max(zero,sign(unun,vgang1-dTsndz         ))
      dT2_OK = max(zero,sign(unun,vgang2-dTsndz         ))
      dT3xOK = max(zero,sign(unun,vgang3-dTsndz         ))
      dT3_OK =                    dT3xOK  * (1. - dT2_OK)
      dT4xOK = max(zero,sign(unun,vgang4-dTsndz         ))
      dT4_OK =                    dT4xOK  * (1. - dT3_OK) &
            * (1. - dT2_OK)
      dT4nOK =              (1. - dT4xOK) * (1. - dT3_OK) &
            * (1. - dT2_OK)

  ! +-- Influence de la Temperature /Temperature Influence
  ! +   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      AngSno = &
            T1__OK & ! 11
            *(T2__OK*(vtang4+vtang5*(TfSnow       -TsisSV(ikl,isn)) & ! 12
            /vtang6) & !
            +T3__OK*(vtang7-vtang8*(TfSnow-vtang2-TsisSV(ikl,isn)) & ! 13
            /vtang9) & !
            +T3_nOK*(vtanga-vtangb*(TfSnow-vtang3-TsisSV(ikl,isn)) & ! 14
            /vtangc)) & !

  ! +-- Influence de la Masse Volumique /Density Influence
  ! +   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            * ro1_OK &
            *(   ro2_OK*(1. - (ro_dry(ikl,isn)-vrang2) & !
            /(vrang1-vrang2)) & !
            +1.-ro2_OK                                ) & !

  ! +-- Influence du Gradient de Temperature /Temperature Gradient Influence
  ! +   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            *(   dT1_OK*(dT2_OK*vgang5*(dTsndz-vgang6) & ! 15
            /(vgang2-vgang6) & !
            +dT3_OK*vgang7 & ! 16
            +dT4_OK*vgang9 & ! 17
            +dT4nOK*vgangb                ) & ! 18
            +1.-dT1_OK                                ) & !
            + ro1_OK &
            *    dT1_OK*(dT3_OK*vgang8*(dTsndz-vgang2) &
            /(vgang3-vgang2) &
            +dT4_OK*vganga*(dTsndz-vgang3) &
            /(vgang4-vgang3) &
            +dT4nOK*vgangc*(dTsndz-vgang4) &
            /(vgang1-vgang4))

      G2_hds = G2snSV(ikl,isn) + 1.d2 *AngSno*vfi     *frac_j


  ! +--New Properties
  ! +  --------------

      G1_bak          = G1snSV(ikl,isn)
      G2_bak          = G2snSV(ikl,isn)

      G1snSV(ikl,isn) = Wet_OK * (    OK__wd             *G1__wd & !  1
            +(1.-OK__wd)*    OK__ws *G1__ws & !  2
            +(1.-OK__wd)*(1.-OK__ws)*G1_bak) & !  3
            +(1. - Wet_OK) & !
            *(    OKlowT  *(    OK_ldd             *G1_ldd & !  4
            +(1.-OK_ldd)            *G1_lds) & !  5
            + OKmidT  *(    OK_mdd             *G1_mdd & !  6
            +(1.-OK_mdd)            *G1_mds) & !  7
            + OKhigT  *(    OK_hdd             *G1_hdd & !  8
            +(1.-OK_hdd)*    OK_hds *G1_hds & !  9
            +(1.-OK_hdd)*(1.-OK_hds)*G1_bak))  ! 10

  !XF
  if(G1snSV(ikl,isn)<0.1) &
        G2_hds = G2snSV(ikl,isn) + 1.d1 *AngSno*vfi     *frac_j
  !XF


      G2snSV(ikl,isn) = Wet_OK * (    OK__wd             *G2__wd & !  1
            +(1.-OK__wd)*    OK__ws *G2_bak & !  2
            +(1.-OK__wd)*(1.-OK__ws)*G2__ws) & !  3
            +(1. - Wet_OK) & !
            *(    OKlowT  *(    OK_ldd             *G2_ldd & !  4
            +(1.-OK_ldd)            *G2_bak) & !  5
            + OKmidT  *(    OK_mdd             *G2_mdd & !  6
            +(1.-OK_mdd)            *G2_bak) & !  7
            + OKhigT  *(    OK_hdd             *G2_hdd & !  8
            +(1.-OK_hdd)*    OK_hds *G2_bak & !  9
            +(1.-OK_hdd)*(1.-OK_hds)*G2_hds))  ! 10

  ! +--Snow Properties: IO Set Up
  ! +  ~~~~~~~~~~~~~~~~~~~~~~~~~~
  ! #vp     G_curr( 1) =     Wet_OK             *    OK__wd
  ! #vp     G_curr( 2) =     Wet_OK             *(1.-OK__wd)*    OK__ws
  ! #vp     G_curr( 3) =     Wet_OK             *(1.-OK__wd)*(1.-OK__ws)
  ! #vp     G_curr( 4) = (1.-Wet_OK)*    OKlowT *    OK_ldd
  ! #vp     G_curr( 5) = (1.-Wet_OK)*    OKlowT *(1.-OK_ldd)
  ! #vp     G_curr( 6) = (1.-Wet_OK)*    OKmidT *    OK_mdd
  ! #vp     G_curr( 7) = (1.-Wet_OK)*    OKmidT *(1.-OK_mdd)
  ! #vp     G_curr( 8) = (1.-Wet_OK)*    OKhigT *    OK_hdd
  ! #vp     G_curr( 9) = (1.-Wet_OK)*    OKhigT *(1.-OK_hdd)*    OK_hds
  ! #vp     G_curr(10) = (1.-Wet_OK)*    OKhigT *(1.-OK_hdd)*(1.-OK_hds)
  ! #vp     G_curr(11) =     T1__OK                         * G_curr(10)
  ! #vp     G_curr(12) =     T2__OK                         * G_curr(10)
  ! #vp     G_curr(13) =     T3__OK                         * G_curr(10)
  ! #vp     G_curr(14) =     T3_nOK                         * G_curr(10)
  ! #vp     G_curr(15) =     ro1_OK*     dT1_OK *    dT2_OK * G_curr(10)
  ! #vp     G_curr(16) =     ro1_OK*     dT1_OK *    dT3_OK * G_curr(10)
  ! #vp     G_curr(17) =     ro1_OK*     dT1_OK *    dT4_OK * G_curr(10)
  ! #vp     G_curr(18) =     ro1_OK*     dT1_OK *    dT4nOK * G_curr(10)

  ! #vp     Gcases( 1) = max(Gcases( 1),G_curr( 1))
  ! #vp     Gcases( 2) = max(Gcases( 2),G_curr( 2))
  ! #vp     Gcases( 3) = max(Gcases( 3),G_curr( 3))
  ! #vp     Gcases( 4) = max(Gcases( 4),G_curr( 4))
  ! #vp     Gcases( 5) = max(Gcases( 5),G_curr( 5))
  ! #vp     Gcases( 6) = max(Gcases( 6),G_curr( 6))
  ! #vp     Gcases( 7) = max(Gcases( 7),G_curr( 7))
  ! #vp     Gcases( 8) = max(Gcases( 8),G_curr( 8))
  ! #vp     Gcases( 9) = max(Gcases( 9),G_curr( 9))
  ! #vp     Gcases(10) = max(Gcases(10),G_curr(10))
  ! #vp     Gcases(11) = max(Gcases(11),G_curr(11))
  ! #vp     Gcases(12) = max(Gcases(12),G_curr(12))
  ! #vp     Gcases(13) = max(Gcases(13),G_curr(13))
  ! #vp     Gcases(14) = max(Gcases(14),G_curr(14))
  ! #vp     Gcases(15) = max(Gcases(15),G_curr(15))
  ! #vp     Gcases(16) = max(Gcases(16),G_curr(16))
  ! #vp     Gcases(17) = max(Gcases(17),G_curr(17))
  ! #vp     Gcases(18) = max(Gcases(18),G_curr(18))

  ! +--Snow Properties: IO
  ! +  ~~~~~~~~~~~~~~~~~~~
  ! #vp     IF          (isn    .le.     isnoSV(ikl))
  ! #vp.    write(47,471)isn            ,isnoSV(ikl)                    ,
  ! #vp.                 TsisSV(ikl,isn),ro__SV(ikl,isn),eta_SV(ikl,isn),
  ! #vp.                 G1_bak         ,G2_bak         ,istoSV(ikl,isn),
  ! #vp.                 dTsndz,
  ! #vp.                (       k ,k=1,18),
  ! #vp.                (G_curr(k),k=1,18),
  ! #vp.                (Gcases(k),k=1,18),
  ! #vp.                 Wet_OK,OK__wd,G1__wd,G2__wd,
  ! #vp.                     1.-OK__wd,OK__ws,G1__ws,1.-OK__ws,G2__ws,
  ! #vp.              1.-Wet_OK,OKlowT,OK_ldd,G1_ldd,          G2_ldd,
  ! #vp.                            1.-OK_ldd,G1_lds,
  ! #vp.                        OKmidT,OK_mdd,G1_mdd,          G1_mdd,
  ! #vp.                            1.-OK_mdd,G1_mds,
  ! #vp.                        OKhigT,OK_hdd,G1_hdd,          G2_hdd,
  ! #vp.                            1.-OK_hdd,OK_hds,          G1_hds,
  ! #vp.                                             1.-OK_hds,G2_hds,
  ! #vp.                 G1snSV(ikl,isn),
  ! #vp.                 G2snSV(ikl,isn)

    END DO
  END DO
  !!$OMP END PARALLEL DO

  ! +-- 2. Mise a Jour Variables Historiques (Cas non dendritique)
  ! +      Update of the historical Variables
  ! +      =======================================================

  IF (vector)                                                   THEN
  !XF
     DO  ikl=1,knonv
      DO isn=1,isnoSV(ikl)
      SphrOK = max(zero,sign(unun,       G1snSV(ikl,isn)))
      H1a_OK = max(zero,sign(unun,vsphe4-G1snSV(ikl,isn)))
      H1b_OK =     1   - min(1   ,       istoSV(ikl,isn))
      H1__OK =                    H1a_OK*H1b_OK
      H23aOK = max(zero,sign(unun,vsphe4-G1_dSV &
            +G1snSV(ikl,isn)))
      H23bOK = max(zero,sign(unun,etaSno(ikl,isn) &
            /max(epsi,dzsnSV(ikl,isn)) &
            -vtelv1         ))
      H23_OK =                    H23aOK*H23bOK
      H2__OK =     1   - min(1   ,       istoSV(ikl,isn))
      H3__OK =     1   - min(1   ,   abs(istoSV(ikl,isn)-istdSV(1)))
      H45_OK = max(zero,sign(unun,TfSnow-TsisSV(ikl,isn)+epsi))
      H4__OK =     1   - min(1   ,   abs(istoSV(ikl,isn)-istdSV(2)))
      H5__OK =     1   - min(1   ,   abs(istoSV(ikl,isn)-istdSV(3)))

      HISupd          = &
            SphrOK*(H1__OK                             *istdSV(1) &
            +(1.-H1__OK)*    H23_OK         *(H2__OK*istdSV(2) &
            +H3__OK*istdSV(3)) &
            +(1.-H1__OK)*(1.-H23_OK) *H45_OK*(H4__OK*istdSV(4) &
            +H5__OK*istdSV(5)))
      istoSV(ikl,isn) =   HISupd  + &
            (1.-min(unun,HISupd))               *istoSV(ikl,isn)
    END DO
    END DO
  ELSE


  ! +-- 2. Mise a Jour Variables Historiques (Cas non dendritique)
  ! +      Update of the historical Variables
  ! +      =======================================================

    DO ikl=1,knonv
    DO isn=iiceSV(ikl),isnoSV(ikl)
      IF  (G1snSV(ikl,isn).ge.0.)                               THEN
        IF(G1snSV(ikl,isn).lt.vsphe4.and.istoSV(ikl,isn).eq.0)  THEN
               istoSV(ikl,isn)=istdSV(1)
        ELSEIF(G1_dSV-G1snSV(ikl,isn)         .lt.vsphe4.and. &
              etaSno(ikl,isn)/dzsnSV(ikl,isn).gt.vtelv1)       THEN
          IF  (istoSV(ikl,isn).eq.0) &
                istoSV(ikl,isn)=   istdSV(2)
          IF  (istoSV(ikl,isn).eq.istdSV(1)) &
                istoSV(ikl,isn)=   istdSV(3)
        ELSEIF(TsisSV(ikl,isn).lt.TfSnow)                       THEN
          IF  (istoSV(ikl,isn).eq.istdSV(2)) &
                istoSV(ikl,isn)=   istdSV(4)
          IF  (istoSV(ikl,isn).eq.istdSV(3)) &
                istoSV(ikl,isn)=   istdSV(5)
        END IF
      END IF
    END DO
    END DO
  END IF


  ! +-- 3. Tassement mecanique /mechanical Settlement
  ! +      ==========================================

    DO ikl=1,knonv
      SnMass(ikl) = 0.
    END DO
  !XF
  DO ikl=1,knonv

    smb_old = 0.
     zn_old = 0
    DO isn  = 1, isnoSV(ikl)
    smb_old = smb_old + dzsnSV(ikl,isn) *ro__SV(ikl,isn)
     zn_old = zn_old  + dzsnSV(ikl,isn)
    ENDDO

    DO   isn=isnoSV(ikl),1,-1
      dSnMas     = 100.*dzsnSV(ikl,isn)*ro_dry(ikl,isn)
      SnMass(ikl)=      SnMass(ikl)+0.5*dSnMas
      ViscSn     =      vvisc1         *vvisc2 &
            *exp(vvisc3           *ro_dry(ikl,isn) &
            +vvisc4*abs(TfSnow-TsisSV(ikl,isn))) &
            *ro_dry(ikl,isn)/rovisc

  ! +-- Changement de Viscosite si Teneur en Eau liquide
  ! +   Change of the Viscosity if liquid Water Content
  ! +   ------------------------------------------------

      OK_Liq     =    max(zero,sign(unun,etaSno(ikl,isn)-epsi))
      OK_Ang     =    max(zero,sign(unun,vgran6-G1snSV(ikl,isn))) &
            *(1-min(1   , abs(istoSV(ikl,isn)-istdSV(1))))
  ! #wp     IF (G1snSV(ikl,isn).gt.0..AND.G1snSV(ikl,isn).lt.vsphe4
  ! #wp.                             .AND.istoSV(ikl,isn).eq.     0)
  ! #wp.    THEN
  ! #wp       write(6,*) ikl,isn,' G1,G2,hist,OK_Ang  ',
  ! #wp.          G1snSV(ikl,isn), G2snSV(ikl,isn),istoSV(ikl,isn),OK_Ang
  ! #wp       stop "Grains anguleux mal d?finis"
  ! #wp     END IF
      OKxLiq     =    max(zero,sign(unun,vtelv1-etaSno(ikl,isn) &
            /max(epsi,dzsnSV(ikl,isn)))) &
            *    max(0   ,sign(1   ,istoSV(ikl,isn) &
            -istdSV(1)      ))
      ViscSn     = &
            ViscSn*(    OK_Liq/(vvisc5+vvisc6*etaSno(ikl,isn) &
            /max(epsi,dzsnSV(ikl,isn))) &
            +(1.-OK_Liq)                               ) &
            *(    OK_Ang*exp(min(ADSdSV,G2snSV(ikl,isn)-vdiam4)) &
            +(1.-OK_Ang)                                       ) &
            *(    OKxLiq*        vvisc7 &
            +(1.-OKxLiq)              )


  ! +-- Calcul nouvelle Epaisseur / new Thickness
  ! +   -----------------------------------------

      dzsnew         = &
            dzsnSV(ikl,isn) &
            *max(vdz3, &
            (unun-dt__SV2*max(SnMass(ikl)*cos(slopSV(ikl)),unun) &
            /max(ViscSn                      ,epsi)))
      rosnew         = ro__SV(ikl,isn) *dzsnSV(ikl,isn) &
            /max(1e-10,dzsnew)
      rosmax         = 1.   /( (1. -eta_SV(ikl,isn)) /ro_Ice &
            +      eta_SV(ikl,isn)  /ro_Wat)
      rosnew         = min(rosnew ,rosmax)
      dzsnew         = dzsnSV(ikl,isn) *ro__SV(ikl,isn) &
            /max(1e-10,rosnew)
      ro__SV(ikl,isn)= rosnew
      dzsnSV(ikl,isn)= dzsnew
      ro_dry(ikl,isn)= ro__SV(ikl,isn)*(1.-eta_SV(ikl,isn))*1.e-3
  ! +...    ro_dry: Dry Density (g/cm3)
  ! +
      SnMass(ikl)    = SnMass(ikl)+dSnMas*0.5
    END DO

    smb_new = 0.
    DO isn  = 1, isnoSV(ikl)
    smb_new = smb_new + dzsnSV(ikl,isn) *ro__SV(ikl,isn)
    ENDDO

    isn=1
    if (dzsnSV(ikl,isn)>0.and.ro__SV(ikl,isn)>0) then
    dzsnSV(ikl,isn) = dzsnSV(ikl,isn) +0.9999*(smb_old-smb_new) &
          / ro__SV(ikl,isn)
    endif

     zn_new = 0
    DO isn  = 1, isnoSV(ikl)
     zn_new = zn_new  + dzsnSV(ikl,isn)
    ENDDO
    zn4_SV(ikl) = zn4_SV(ikl) + (zn_new - zn_old)

  END DO



  return
end subroutine sisvat_gsn
