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

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

  
      IMPLICIT NONE
  
     
 
C +--INPUT/OUTPUT
C +  ------------
 
 
C +--OUTPUT
C +  ------
 
      integer   dt__SV2
 
 
C +--Local  Variables
C +  ================
 
      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)
 
C +--Snow Properties: IO
C +  ~~~~~~~~~~~~~~~~~~~
! #vp real      G_curr(18),Gcases(18)
! #vp common   /GSnLOC/    Gcases
! #wp real                 D__MAX
! #wp common   /GSnMAX/    D__MAX
 
 
C +--DATA
C +  ====
 
      data       vector/.TRUE./               ! Vectorization  Switch
      data       vdent1/ 0.5e8/               ! Wet Snow Metamorphism
cXF                      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
 
 
C +--DATA (Coefficient Fonction fort Gradient Marbouty)
C +  --------------------------------------------------
 
      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/               !
 
 
C +-- 1. Metamorphoses dans les Strates
C +      Metamorphism
C +      ==============================
 
      dt__SV2= dt__SV
      frac_j = dt__SV2 / 86400.                        ! Time Step [Day]
 
      zn4_SV = 0
 
 
C +-- 1.1 Initialisation: teneur en eau liquide et gradient de temperature
C +   ------------------  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
 
c!$OMP PARALLEL DO default(firstprivate)
c!$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))
C +...    Factor 1.d-2 for Conversion K/m --> K/cm
 
 
C +-- 1.2 Metamorphose humide
C +       Wet Snow Metamorphism
C +       ---------------------
 
          Wet_OK = max(zero,sign(unun,eta_SV(ikl,isn)-epsi))
 
 
C +--     Vitesse de diminution de la dendricite
C +       Rate of the dendricity decrease
C +       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          sWater=1.d-1*ro__SV(ikl,isn)*eta_SV(ikl,isn)
     .       /max(epsi,ro_dry(ikl,isn))
C +...    sWater:Water Content [%]
C +              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))
 
C +-- 1.2.1 Cas dendritique/dendritic Case
C +         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          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
 
C +-- 1.2.2 Cas non dendritique non completement spherique
C +         Evolution de la Sphericite seulement.
C +         Non dendritic and not completely spheric Case
C +         Evolution of    Sphericity only (not size)
C +         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          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)
 
C +-- 1.2.3 Cas non dendritique et spherique / non dendritic and spheric
C +         Evolution de la Taille seulement / Evolution of Size only
C +         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          G2__ws =  husi_0
     .           *( husi_1
     .            *(husi_2 *( G2snSV(ikl,isn)/husi_0)**3
     .                      +(vtail1 +vtail2 *exp1Wa    )*dt__SV2))
     .           ** husi_3
 
 
C +-- 1.3 Metamorposes seches / Dry Metamorphism
C +       --------------------------------------
 
 
C +-- 1.3.1 Calcul Metamorphose faible/low Gradient (0.00-0.05 deg/cm)
C +         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          OKlowT=max(zero,                            !
     .               sign(unun, vgrat1                !
     .                         -dTsndz         ))     !
 
          facVap=exp(vvap1/TsisSV(ikl,isn))
 
C +-- 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
 
C +-- 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)
 
C +-- 1.3.2 Calcul Metamorphose Gradient Moyen/Moderate (0.05-0.15)
C +         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          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
 
C +-- 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
 
C +-- 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)
 
C +-- 1.3.3 Calcul Metamorphose fort / high Gradient
C +   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          facVap=vdent1*exp(vvap1/TsisSV(ikl,isn))
     .                 *   (1.e2 *dTsndz)**vvap2
 
C +-- 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
C +                                                   ! 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
 
C +-- 1.3.3.2 Cas non dendritique non completement anguleux.
C +           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)
 
C +-- 1.3.3.3 Cas non dendritique et anguleux
C +           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)
 
C +-- Influence de la Temperature /Temperature Influence
C +   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          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))                                      !
 
C +-- Influence de la Masse Volumique /Density Influence
C +   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     .    * ro1_OK
     .        *(   ro2_OK*(1. - (ro_dry(ikl,isn)-vrang2)                 !
     .                                  /(vrang1-vrang2))                !
     .         +1.-ro2_OK                                )               !
 
C +-- Influence du Gradient de Temperature /Temperature Gradient Influence
C +   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     .        *(   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
 
 
C +--New Properties
C +  --------------
 
          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
 
cXF
      if(G1snSV(ikl,isn)<0.1)
     .    G2_hds = G2snSV(ikl,isn) + 1.d1 *AngSno*vfi     *frac_j
cXF
 
 
          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
 
C +--Snow Properties: IO Set Up
C +  ~~~~~~~~~~~~~~~~~~~~~~~~~~
! #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))
 
C +--Snow Properties: IO
C +  ~~~~~~~~~~~~~~~~~~~
! #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
c!$OMP END PARALLEL DO
 
C +-- 2. Mise a Jour Variables Historiques (Cas non dendritique)
C +      Update of the historical Variables
C +      =======================================================
 
      IF (vector)                                                   THEN
cXF
         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
 
 
C +-- 2. Mise a Jour Variables Historiques (Cas non dendritique)
C +      Update of the historical Variables
C +      =======================================================
 
        DO ikl=1,knonv
        DO isn=iiceSV(ikl),isnoSV(ikl)
          IF  (G1snSV(ikl,isn)>=0.)                               THEN
            IF(G1snSV(ikl,isn)<vsphe4.and.istoSV(ikl,isn)==0)  THEN
                   istoSV(ikl,isn)=istdSV(1)
            ELSEIF(G1_dSV-G1snSV(ikl,isn)         <vsphe4.and.
     .             etaSno(ikl,isn)/dzsnSV(ikl,isn)>vtelv1)       THEN
              IF  (istoSV(ikl,isn)==0)
     .             istoSV(ikl,isn)=   istdSV(2)
              IF  (istoSV(ikl,isn)==istdSV(1))
     .             istoSV(ikl,isn)=   istdSV(3)
            ELSEIF(TsisSV(ikl,isn)<TfSnow)                       THEN
              IF  (istoSV(ikl,isn)==istdSV(2))
     .             istoSV(ikl,isn)=   istdSV(4)
              IF  (istoSV(ikl,isn)==istdSV(3))
     .             istoSV(ikl,isn)=   istdSV(5)
            END IF
          END IF
        END DO
        END DO
      END IF
 
 
C +-- 3. Tassement mecanique /mechanical Settlement
C +      ==========================================
 
        DO ikl=1,knonv
          SnMass(ikl) = 0.
        END DO
cXF
      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
 
C +-- Changement de Viscosite si Teneur en Eau liquide
C +   Change of the Viscosity if liquid Water Content
C +   ------------------------------------------------
 
          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)              )
 
 
C +-- Calcul nouvelle Epaisseur / new Thickness
C +   -----------------------------------------
 
          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
C +...    ro_dry: Dry Density (g/cm3)
C +
          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
