SUBROUTINE SnOptP(jjtime)

  ! +------------------------------------------------------------------------+
  ! | MAR/SISVAT   SnOptP                                    12-08-2019  MAR |
  ! |   SubRoutine SnOptP computes the Snow Pack optical Properties          |
  ! +------------------------------------------------------------------------+
  ! |                                                                        |
  ! |   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                      |
  ! |   ^^^^^    ispiSV   = 0,...,nsno: Uppermost Superimposed Ice Layer     |
  ! |                                                                        |
  ! |                                                                        |
  ! |   INPUT:   G1snSV   : Dendricity (<0) or Sphericity (>0) of Snow Layer |
  ! |   ^^^^^    G2snSV   : Sphericity (>0) or Size            of Snow Layer |
  ! |            agsnSV   : Snow       Age                             [day] |
  ! |            ro__SV   : Snow/Soil  Volumic Mass                  [kg/m3] |
  ! |            eta_SV   : Water      Content                       [m3/m3] |
  ! |            rusnSV   : Surficial  Water   Thickness   [kg/m2] .OR. [mm] |
  ! |            SWS_SV   : Surficial  Water   Status                        |
  ! |            dzsnSV   : Snow       Layer   Thickness                 [m] |
  ! |                                                                        |
  ! |            albssv   : Soil       Albedo                            [-] |
  ! |            zzsnsv   : Snow       Pack    Thickness                 [m] |
  ! |                                                                        |
  ! |   OUTPUT:  albisv   : Snow/Ice/Water/Soil Integrated Albedo        [-] |
  ! |   ^^^^^^   sEX_sv   : Verticaly Integrated Extinction Coefficient      |
  ! |            DOPsnSV   : Snow Optical diameter [m]                       |
  ! |                                                                        |
  ! |   Internal Variables:                                                  |
  ! |   ^^^^^^^^^^^^^^^^^^                                                   |
  ! |            SnOpSV   : Snow Grain optical Size                      [m] |
  ! |            EX1_sv   : Integrated Snow Extinction (0.3--0.8micr.m)      |
  ! |            EX2_sv   : Integrated Snow Extinction (0.8--1.5micr.m)      |
  ! |            EX3_sv   : Integrated Snow Extinction (1.5--2.8micr.m)      |
  ! |                                                                        |
  ! |   METHODE:    Calcul de la taille optique des grains ? partir de       |
  ! |   ^^^^^^^    -leur type decrit par les deux variables descriptives     |
  ! |                    continues sur la plage -99/+99 passees en appel.    |
  ! |              -la taille optique (1/10mm) des etoiles,                  |
  ! |                                          des grains fins et            |
  ! |                                          des jeunes faces planes       |
  ! |                                                                        |
  ! |   METHOD:     Computation of the optical diameter of the grains        |
  ! |   ^^^^^^      described with the CROCUS formalism G1snSV / G2snSV      |
  ! |                                                                        |
  ! |   REFERENCE: Brun et al.      1989, J. Glaciol 35 pp. 333--342         |
  ! |   ^^^^^^^^^  Brun et al.      1992, J. Glaciol 38 pp.  13-- 22         |
  ! |              Eric Martin Sept.1996                                     |
  ! |                                                                        |
  ! |                                                                        |
  ! +------------------------------------------------------------------------+




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


  use VARphy
  use VAR_SV
  use VARdSV
  use VARxSV
  use VARySV
  use VARtSV
  USE surface_data, only: iflag_albcalc,correc_alb

  IMPLICIT NONE


  ! + -- INPUT
  integer :: jjtime

  ! +--Internal Variables
  ! +  ==================

  real :: coalb1(knonv)                      ! weighted Coalbedo, Vis.
  real :: coalb2(knonv)                      ! weighted Coalbedo, nIR 1
  real :: coalb3(knonv)                      ! weighted Coalbedo, nIR 2
  real :: coalbm                             ! weighted Coalbedo, mean
  real :: sExt_1(knonv)                      ! Extinction Coeff., Vis.
  real :: sExt_2(knonv)                      ! Extinction Coeff., nIR 1
  real :: sExt_3(knonv)                      ! Extinction Coeff., nIR 2
  real :: SnOpSV(knonv,      nsno)           ! Snow Grain optical Size
  ! #AG real      agesno

  integer :: isn   ,ikl   ,isn1, i
  real :: sbeta1,sbeta2,sbeta3,sbeta4,sbeta5
  real :: AgeMax,AlbMin,HSnoSV,HIceSV,doptmx,SignG1,Sph_OK
  real :: dalbed,dalbeS,dalbeW
  real :: bsegal,czemax,csegal,csza
  real :: RoFrez,DiffRo,SignRo,SnowOK,OpSqrt
  real :: albSn1,albIc1,a_SnI1,a_SII1
  real :: albSn2,albIc2,a_SnI2,a_SII2
  real :: albSn3,albIc3,a_SnI3,a_SII3
  real :: albSno,albIce,albSnI,albSII,albWIc,albmax
  real :: doptic,Snow_H,SIce_H,SnownH,SIcenH
  real :: exarg1,exarg2,exarg3,sign_0,sExt_0
  real :: albedo_old,albCor
  real :: ro_ave,dz_ave,minalb
  real :: l1min,l1max,l2min,l2max,l3min,l3max
  real :: l6min(6), l6max(6), albSn6(6), a_SII6(6)
  real :: lmintmp,lmaxtmp,albtmp

  ! +--Local   DATA
  ! +  ============

  ! +--For the computation of the solar irradiance extinction in snow
  ! +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  data      sbeta1/0.0192/,sbeta2/0.4000/,sbeta3/0.1098/
  data      sbeta4/1.0000/
  data      sbeta5/2.00e1/

  ! +--Snow Age Maximum (Taiga, e.g. Col de Porte)
  ! +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  data      AgeMax  /60.0/
  ! +...          AgeMax:  Snow Age Maximum                              [day]

  data      AlbMin  /0.94/
  ! +...          AlbMin:  Albedo   Minimum / visible (0.3--0.8 micrometers)

  data      HSnoSV  /0.01/
  ! +...          HSnoSV:  Snow     Thickness through witch
  ! +                      Albedo is interpolated to Ice  Albedo
  data      HIceSV  /0.10/
  ! +...          HIceSV:  Snow/Ice Thickness through witch
  ! +                      Albedo is interpolated to Soil Albedo
  data      doptmx  /2.3e-3/
  ! +...          doptmx:  Maximum  optical Diameter    (pi * R**2)        [m]
  ! +
  data      czeMAX  /0.173648178/  ! 80.deg (Segal et al., 1991 JAS)

  data      bsegal  /4.00       /  !
  data      albmax  /0.99       /  ! Albedo max

  ! +-- wavelength limits [m] for each broad band

  data      l1min/400.0e-9/,l1max/800.0e-9/
  data      l2min/800.0e-9/,l2max/1500.0e-9/
  data      l3min/1500.0e-9/,l3max/2800.0e-9/

  data      l6min/185.0e-9,250.0e-9,400.0e-9, &
        690.0e-9,1190.0e-9,2380.0e-9/
  data      l6max/250.0e-9,400.0e-9, &
        690.0e-9,1190.0e-9,2380.0e-9,4000.0e-9/


  ! +--Snow Grain optical Size
  ! +  =======================

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

      G2snSV(ikl,isn) =  max(epsi,G2snSV(ikl,isn))
  ! +...    Avoid non physical Values

      SignG1          = sign(unun,G1snSV(ikl,isn))
      Sph_OK          =  max(zero,SignG1)

      SnOpSV(ikl,isn) =   1.e-4 * &
  ! +...    SI:           (from 1/10 mm to m)


  ! +--Contribution of Non Dendritic Snow
  ! +  ----------------------------------

            (    Sph_OK *(      G2snSV(ikl,isn)*G1snSV(ikl,isn)/G1_dSV &
            +max(demi*G2snSV(ikl,isn),DFcdSV) &
            *(unun-G1snSV(ikl,isn)                /G1_dSV)) &


  ! +--Contribution of     Dendritic Snow
  ! +  ----------------------------------

            +(1.-Sph_OK)*(     -G1snSV(ikl,isn)*DDcdSV         /G1_dSV &
            +(unun+G1snSV(ikl,isn)                /G1_dSV) &
            *    (G2snSV(ikl,isn)*DScdSV         /G1_dSV &
            +(unun-G2snSV(ikl,isn)                /G1_dSV) &
            *DFcdSV                 )))
      SnOpSV(ikl,isn) =  max(zero,SnOpSV(ikl,isn))

  ! + --For outputs (Etienne)
  ! + ------------------------
      DOPsnSV(ikl,isn)=SnOpSV(ikl,isn)
    END DO
  END DO




  ! +--Snow/Ice Albedo
  ! +  ===============



  ! +--Uppermost effective Snow Layer
  ! +  ------------------------------

    DO ikl=1,knonv

      isn   =  max(iun,isnoSV(ikl))

      SignRo = sign(unun, rocdSV - ro__SV(ikl,isn))
      SnowOK =  max(zero,SignRo) ! Ice Density Threshold

      OpSqrt = sqrt(SnOpSV(ikl,isn))

  !CA    +--Correction of snow albedo for Antarctica/Greenland
  !CA       --------------------------------------------------


      albCor = correc_alb
  ! #GL     albCor = 1.01
  ! #AC    albCor = 1.01


      IF (iflag_albcalc >= 1) THEN  ! Albedo calculation according to Kokhanovsky and Zege 2004

      dalbed = 0.0
      doptic=SnOpSV(ikl,isn)
      csza=coszSV(ikl)

      CALL albedo_kokhanovsky(l1min,l1max,csza,doptic,albSn1)
      CALL albedo_kokhanovsky(l2min,l2max,csza,doptic,albSn2)
      CALL albedo_kokhanovsky(l3min,l3max,csza,doptic,albSn3)

      DO i=1,6
         lmintmp=l6min(i)
         lmaxtmp=l6max(i)
         CALL albedo_kokhanovsky(lmintmp,lmaxtmp,csza,doptic,albtmp)
         albSn6(i)=albtmp
      ENDDO


      ELSE ! Default calculation in SISVAT

  !    Zenith Angle Correction (Segal et al.,         1991, JAS 48, p.1025)
  !--------------------------- (Wiscombe & Warren, dec1980, JAS   , p.2723)
  !                            (Warren,               1982,  RG   , p.  81)
                         ! -------------------------------------------------

      dalbed = 0.0

      csegal = max(czemax                   ,coszSV(ikl))
  ! #cz     dalbeS =   ((bsegal+unun)/(unun+2.0*bsegal*csegal)
  ! #cz.                -       unun                          )*0.32
  ! #cz.              /  bsegal
  ! #cz     dalbeS = max(dalbeS,zero)
  ! #cz     dalbed =     dalbeS      *       min(1,isnoSV(ikl))

      dalbeW =(0.64 - csegal  )*0.0625  ! Warren 1982, RevGeo, fig.12b
                                        ! ! 0.0625 = 5% * 1/0.8,   p.81
                                        ! ! 0.64   = cos(50)
      dalbed =     dalbeW      *       min(1,isnoSV(ikl))
  !-------------------------------------------------------------------------

      albSn1 =  0.96-1.580*OpSqrt
      albSn1 =  max(albSn1,AlbMin)

      albSn1 =  max(albSn1,zero)
      albSn1 =  min(albSn1*albCor,unun)

      albSn2 =  0.95-15.40*OpSqrt
      albSn2 =  max(albSn2,zero)
      albSn2 =  min(albSn2*albCor,unun)

      doptic =  min(SnOpSV(ikl,isn),doptmx)
      albSn3 =  346.3*doptic -32.31*OpSqrt +0.88
      albSn3 =  max(albSn3,zero)
      albSn3 =  min(albSn3*albCor,unun)

      albSn6(1:3)=albSn1
      albSn6(4:6)=albSn2

      ! !snow albedo corection if wetsnow
  ! #GL     albSn1 =  albSn1*max(0.9,(1.-1.5*eta_SV(ikl,isn)))
  ! #GL     albSn2 =  albSn2*max(0.9,(1.-1.5*eta_SV(ikl,isn)))
  ! #GL     albSn3 =  albSn3*max(0.9,(1.-1.5*eta_SV(ikl,isn)))

      ENDIF


      albSno =  So1dSV*albSn1 &
            +  So2dSV*albSn2 &
            +  So3dSV*albSn3

  !XF
      minalb =  (aI2dSV + (aI3dSV -aI2dSV) &
            *  (ro__SV(ikl,isn)-ro_Ice)/(roSdSV-ro_Ice))
      minalb =  min(aI3dSV,max(aI2dSV,minalb)) ! pure/firn albedo

      SnowOK =  SnowOK*max(zero,sign(unun,     albSno-minalb))
      albSn1 =  SnowOK*albSn1+(1.0-SnowOK)*max(albSno,minalb)
      albSn2 =  SnowOK*albSn2+(1.0-SnowOK)*max(albSno,minalb)
      albSn3 =  SnowOK*albSn3+(1.0-SnowOK)*max(albSno,minalb)
      albSn6(:) =  SnowOK*albSn6(:)+(1.0-SnowOK)*max(albSno,minalb)


  ! +           ro < roSdSV |          min al > aI3dSV
  ! +  roSdSV < ro < rocdSV | aI2dSV < min al < aI3dSV (fct of density)


  ! +--Snow/Ice Pack Thickness
  ! +  -----------------------

      isn    =    max(min(isnoSV(ikl) ,ispiSV(ikl)),0)
      Snow_H =  zzsnsv(ikl,isnoSV(ikl))-zzsnsv(ikl,isn)
      SIce_H =  zzsnsv(ikl,isnoSV(ikl))
      SnownH =  Snow_H  /  HSnoSV
      SnownH =  min(unun,  SnownH)
      SIcenH =  SIce_H  / (HIceSV)
      SIcenH =  min(unun,  SIcenH)

  ! +       The value of SnownH is set to 1 in case of ice lenses above
  ! +       1m of dry snow (ro<600kg/m3) for using CROCUS albedo

      ! ro_ave =  0.
      ! dz_ave =  0.
      ! SnowOK =  1.
   ! do isn    =  isnoSV(ikl),1,-1
   !    ro_ave =  ro_ave + ro__SV(ikl,isn) * dzsnSV(ikl,isn) * SnowOK
   !    dz_ave =  dz_ave +                   dzsnSV(ikl,isn) * SnowOK
   !    SnowOK =  max(zero,sign(unun,1.-dz_ave))
   ! enddo

   !    ro_ave =  ro_ave / max(dz_ave,epsi)
   !    SnowOK =  max(zero,sign(unun,600.-ro_ave))
   !    SnownH =  SnowOK + SnownH * (1. - SnowOK)

  ! +--Integrated Snow/Ice Albedo: Case of Water on Bare Ice
  ! +  -----------------------------------------------------

      isn    =  max(min(isnoSV(ikl) ,ispiSV(ikl)),0)

      albWIc =  aI1dSV-(aI1dSV-aI2dSV) &
            *  exp(-(rusnSV(ikl) & !
            *  (1. -SWS_SV(ikl) & ! 0 <=> freezing
            *  (1  -min(1,iabs(isn-isnoSV(ikl))))) & ! 1 <=> isn=isnoSV
            /   ru_dSV)**0.50)                        !
      ! albWIc = max(aI1dSV,min(aI2dSV,albWIc+slopSV(ikl)*
  !    .             min(5.,max(1.,dx/10000.))))

      SignRo = sign(unun,ro_Ice-5.-ro__SV(ikl,isn))    ! RoSN<920kg/m3
      SnowOK =  max(zero,SignRo)

      albWIc = (1. - SnowOK) * albWIc + SnowOK &
            * (aI2dSV + (aI3dSV -aI2dSV) &
            * (ro__SV(ikl,isn)-ro_Ice)/(roSdSV-ro_Ice))

  ! +  rocdSV < ro < ro_ice | aI2dSV< al <aI3dSV (fct of density)
  ! +           ro > ro_ice | aI1dSV< al <aI2dSV (fct of superficial water content


  ! +--Integrated Snow/Ice      Albedo
  ! +  -------------------------------

      a_SII1      =     albWIc      +(albSn1-albWIc)     *SnownH
      a_SII1      = min(a_SII1       ,albSn1)

      a_SII2      =     albWIc      +(albSn2-albWIc)     *SnownH
      a_SII2      = min(a_SII2       ,albSn2)

      a_SII3      =     albWIc      +(albSn3-albWIc)     *SnownH
      a_SII3      = min(a_SII3       ,albSn3)

      DO i=1,6
      a_SII6(i)   = albWIc      +(albSn6(i)-albWIc)     *SnownH
      a_SII6(i)   = min(a_SII6(i)       ,albSn6(i))
      ENDDO

  !c #AG     agesno =      min(agsnSV(ikl,isn)          ,AgeMax)
  !c #AG     a_SII1      =     a_SII1      -0.175*agesno/AgeMax
  ! +...                                   Impurities: Col de Porte Parameter.



  ! +--Elsewhere    Integrated Snow/Ice Albedo
  ! +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  ! #cp     ELSE
        albSII =     So1dSV*a_SII1 &
              + So2dSV*a_SII2 &
              + So3dSV*a_SII3
  ! #cp     END IF


  ! +--Integrated Snow/Ice/Soil Albedo
  ! +  -------------------------------

        alb1sv(ikl) =     albssv(ikl) +(a_SII1-albssv(ikl))*SIcenH
        alb1sv(ikl) = min(alb1sv(ikl)  ,a_SII1)

        alb2sv(ikl) =     albssv(ikl) +(a_SII2-albssv(ikl))*SIcenH
        alb2sv(ikl) = min(alb2sv(ikl)  ,a_SII2)

        alb3sv(ikl) =     albssv(ikl) +(a_SII3-albssv(ikl))*SIcenH
        alb3sv(ikl) = min(alb3sv(ikl)  ,a_SII3)

        albisv(ikl) =     albssv(ikl) +(albSII-albssv(ikl))*SIcenH
        albisv(ikl) = min(albisv(ikl)  ,albSII)

        DO i=1,6
        alb6sv(ikl,i) = albssv(ikl) +(a_SII6(i)-albssv(ikl))*SIcenH
        alb6sv(ikl,i) = min(alb6sv(ikl,i)  ,a_SII6(i))
        ENDDO


  ! +--Integrated Snow/Ice/Soil Albedo: Clouds Correction! Greuell & all., 1994
  ! +  --------------------------------------------------! Glob.&t Planet.Change
                                                   ! ! (9):91-114
        alb1sv(ikl) = alb1sv(ikl) + 0.05 *(cld_SV(ikl)-0.5)*SIcenH &
              + dalbed      *    (1.-cld_SV(ikl))
        alb2sv(ikl) = alb2sv(ikl) + 0.05 *(cld_SV(ikl)-0.5)*SIcenH &
              + dalbed      *    (1.-cld_SV(ikl))
        alb3sv(ikl) = alb3sv(ikl) + 0.05 *(cld_SV(ikl)-0.5)*SIcenH &
              + dalbed      *    (1.-cld_SV(ikl))
        alb6sv(ikl,:) =alb6sv(ikl,:)+0.05 *(cld_SV(ikl)-0.5)*SIcenH &
              + dalbed      *    (1.-cld_SV(ikl))
        albisv(ikl) = albisv(ikl) + 0.05 *(cld_SV(ikl)-0.5)*SIcenH &
              + dalbed      *    (1.-cld_SV(ikl))

  ! +--Integrated Snow/Ice/Soil Albedo: Minimum snow albedo = aI1dSV
  ! +  -------------------------------------------------------------

        albedo_old  = albisv(ikl)
        albisv(ikl) = max(albisv(ikl),aI1dSV   * SIcenH &
              + albssv(ikl) *(1.0        - SIcenH))
        alb1sv(ikl) = alb1sv(ikl) - 1.0/3.0 & ! 33 %
              * (albedo_old-albisv(ikl)) / So1dSV
        alb2sv(ikl) = alb2sv(ikl) - 1.0/3.0 & ! 33 %
              * (albedo_old-albisv(ikl)) / So2dSV
        alb3sv(ikl) = alb3sv(ikl) - 1.0/3.0 & ! 33 %
              * (albedo_old-albisv(ikl)) / So3dSV
        alb6sv(ikl,1:3) = alb6sv(ikl,1:3) - 1.0/6.0 & ! 16 %
              * (albedo_old-albisv(ikl)) / (So1dSV/3)
        alb6sv(ikl,4:6) = alb6sv(ikl,4:6) - 1.0/6.0 & ! 16 %
              * (albedo_old-albisv(ikl)) / (So2dSV/3)


  ! +--Integrated Snow/Ice/Soil Albedo: Maximum albedo = 95%
  ! +  -----------------------------------------------------

        albedo_old  = albisv(ikl)
        albisv(ikl) = min(albisv(ikl),0.95)
        alb1sv(ikl) = alb1sv(ikl) - 1.0/3.0 & ! 33 %
              * (albedo_old-albisv(ikl)) / So1dSV
        alb2sv(ikl) = alb2sv(ikl) - 1.0/3.0 & ! 33 %
              * (albedo_old-albisv(ikl)) / So2dSV
        alb3sv(ikl) = alb3sv(ikl) - 1.0/3.0 & ! 33 %
              * (albedo_old-albisv(ikl)) / So3dSV
        alb6sv(ikl,1:3) = alb6sv(ikl,1:3) - 1.0/6.0 & ! 16 %
              * (albedo_old-albisv(ikl)) / (So1dSV/3)
        alb6sv(ikl,4:6) = alb6sv(ikl,4:6) - 1.0/6.0 & ! 16 %
              * (albedo_old-albisv(ikl)) / (So2dSV/3)


    !Sea Ice/snow permanent-interractive prescription from Nemo
    !AO_CK 20/02/2020

    ! ! No check if coupling update since MAR and NEMO albedo are too different
    ! and since MAR albedo is computed on properties that are not in NEMO
    ! ! prescription for each time step with NEMO values

  ! #AO      if (LSmask(ikl) .eq. 0 .and. coupling_ao .eq. .TRUE.) then
  ! #AO       if (AOmask(ikl) .eq. 0) then
  ! #AO       albisv(ikl) =  (1.-AOmask(ikl))* albAOsisv(ikl)
  ! #AO.                    +(AOmask(ikl)*albisv(ikl))
  ! #AO       alb1sv(ikl) =   (1.-AOmask(ikl))* albAOsisv(ikl)
  ! #AO.                    +(AOmask(ikl)*alb1sv(ikl))
  ! #AO       alb2sv(ikl) =   (1.-AOmask(ikl))* albAOsisv(ikl)
  ! #AO.                    +(AOmask(ikl)*alb2sv(ikl))
  ! #AO       alb3sv(ikl) =   (1.-AOmask(ikl))* albAOsisv(ikl)
  ! #AO.                    +(AOmask(ikl)*alb3sv(ikl))
  ! #AO       endif
  ! #AO      endif


        alb1sv(ikl) = min(max(zero,alb1sv(ikl)),albmax)
        alb2sv(ikl) = min(max(zero,alb2sv(ikl)),albmax)
        alb3sv(ikl) = min(max(zero,alb3sv(ikl)),albmax)

        DO i=1,6
            alb6sv(ikl,i) = min(max(zero,alb6sv(ikl,i)),albmax)
        ENDDO
    END DO


  ! +--Extinction Coefficient: Exponential Factor
  ! +  ==========================================

    DO ikl=1,knonv
      sExt_1(ikl)        = 1.
      sExt_2(ikl)        = 1.
      sExt_3(ikl)        = 1.
      sEX_sv(ikl,nsno+1) = 1.

      coalb1(ikl) = (1.          -alb1sv(ikl))*So1dSV
      coalb2(ikl) = (1.          -alb2sv(ikl))*So2dSV
      coalb3(ikl) = (1.          -alb3sv(ikl))*So3dSV
      coalbm      =  coalb1(ikl) +coalb2(ikl) +coalb3(ikl)
      coalb1(ikl) =  coalb1(ikl)              /coalbm
      coalb2(ikl) =  coalb2(ikl)              /coalbm
      coalb3(ikl) =  coalb3(ikl)              /coalbm
    END DO

  !XF

    DO   isn=nsno,1,-1
      DO ikl=1,knonv
        sEX_sv(ikl,isn) = 1.0
       ! !sEX_sv(ikl,isn) = 0.95 ! if MAR is too warm in summer
      END DO
    END DO

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

      SignRo = sign(unun, rocdSV - ro__SV(ikl,isn))
      SnowOK =  max(zero,SignRo) ! Ice Density Threshold

      RoFrez =  1.e-3      * ro__SV(ikl,isn) * (1.0-eta_SV(ikl,isn))

      OpSqrt = sqrt(max(epsi,SnOpSV(ikl,isn)))
      exarg1 =      SnowOK  *1.e2 *max(sbeta1*RoFrez/OpSqrt,sbeta2) &
            +(1.0-SnowOK)           *sbeta5
      exarg2 =      SnowOK  *1.e2 *max(sbeta3*RoFrez/OpSqrt,sbeta4) &
            +(1.0-SnowOK)           *sbeta5
      exarg3 =      SnowOK  *1.e2     *sbeta5 &
            +(1.0-SnowOK)           *sbeta5


  ! +--Integrated Extinction of Solar Irradiance (Normalized Value)
  ! +  ============================================================

      sExt_1(ikl) = sExt_1(ikl) &
            * exp(min(0.0,-exarg1 *dzsnSV(ikl,isn)))
      sign_0      =              sign(unun,eps9   -sExt_1(ikl))
      sExt_0      =               max(zero,sign_0)*sExt_1(ikl)
      sExt_1(ikl) = sExt_1(ikl)                   -sExt_0

      sExt_2(ikl) = sExt_2(ikl) &
            * exp(min(0.0,-exarg2 *dzsnSV(ikl,isn)))
      sign_0      =              sign(unun,eps9   -sExt_2(ikl))
      sExt_0      =               max(zero,sign_0)*sExt_2(ikl)
      sExt_2(ikl) = sExt_2(ikl)                   -sExt_0

      sExt_3(ikl) = sExt_3(ikl) &
            * exp(min(0.0,-exarg3 *dzsnSV(ikl,isn)))
      sign_0      =              sign(unun,eps9   -sExt_3(ikl))
      sExt_0      =               max(zero,sign_0)*sExt_3(ikl)
      sExt_3(ikl) = sExt_3(ikl)                   -sExt_0

      sEX_sv(ikl,isn) = coalb1(ikl) * sExt_1(ikl) &
            + coalb2(ikl) * sExt_2(ikl) &
            + coalb3(ikl) * sExt_3(ikl)
    END DO
  END DO

  DO   isn=0,-nsol,-1
    DO ikl=1,knonv
      sEX_sv(ikl,isn) = 0.0
    END DO
  END DO





end subroutine snoptp


!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE albedo_kokhanovsky(lambdamin,lambdamax, &
        cossza,dopt,albint)
  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  ! Authors: Hajar El Habchi El Fenniri, Etienne Vignon, Cecile Agosta
       ! Ghislain Picard
  ! Routine that calculates the  integrated snow spectral albedo between
  ! lambdamin and lambdamax following Kokhanisky and Zege 2004,
  ! Scattering optics of snow, Applied Optics, Vol 43, No7
  ! Code inspired from the snowoptics package of Ghislain Picard:
  ! https://github.com/ghislainp/snowoptics
  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


  USE VARphy

  IMPLICIT NONE

  ! Inputs
  !--------
  REAL :: lambdamin   ! minimum wavelength for integration [m]
  REAL :: lambdamax   ! maximum wavelength for integration [m]
  REAL :: cossza     ! solar zenith angle cosinus
  REAL :: dopt        ! optical diameter [m]

  !Outputs
  !-------
  REAL :: albint

  ! Local Variables
  !-----------------

  REAL :: ropt,cosalb,norm,Pas
  REAL :: SSA,alpha,gamm,R,cos30,alb30
  INTEGER :: i


  REAL :: B_amp                       ! amplification factor
  PARAMETER (B_amp=1.6)

  REAL :: g_asy                       ! asymetry factor
  PARAMETER (g_asy=0.845)

  INTEGER :: nlambda                  ! length of wavelength vector
  PARAMETER(nlambda=200)

  REAL :: lmin
  PARAMETER(lmin=185.0E-9)

  REAL :: lmax
  PARAMETER(lmax=4000.0E-9)

  REAL :: albmax
  PARAMETER(albmax=0.99)

  REAL :: wavelengths(nlambda)
  REAL :: ni(nlambda)

  DATA wavelengths / 1.85000000e-07, 2.04170854e-07, &
        2.23341709e-07, 2.42512563e-07, &
        2.61683417e-07, 2.80854271e-07, 3.00025126e-07, 3.19195980e-07, &
        3.38366834e-07, 3.57537688e-07, 3.76708543e-07, 3.95879397e-07, &
        4.15050251e-07, 4.34221106e-07, 4.53391960e-07, 4.72562814e-07, &
        4.91733668e-07, 5.10904523e-07, 5.30075377e-07, 5.49246231e-07, &
        5.68417085e-07, 5.87587940e-07, 6.06758794e-07, 6.25929648e-07, &
        6.45100503e-07, 6.64271357e-07, 6.83442211e-07, 7.02613065e-07, &
        7.21783920e-07, 7.40954774e-07, 7.60125628e-07, 7.79296482e-07, &
        7.98467337e-07, 8.17638191e-07, 8.36809045e-07, 8.55979899e-07, &
        8.75150754e-07, 8.94321608e-07, 9.13492462e-07, 9.32663317e-07, &
        9.51834171e-07, 9.71005025e-07, 9.90175879e-07, 1.00934673e-06, &
        1.02851759e-06, 1.04768844e-06, 1.06685930e-06, 1.08603015e-06, &
        1.10520101e-06, 1.12437186e-06, 1.14354271e-06, 1.16271357e-06, &
        1.18188442e-06, 1.20105528e-06, 1.22022613e-06, 1.23939698e-06, &
        1.25856784e-06, 1.27773869e-06, 1.29690955e-06, 1.31608040e-06, &
        1.33525126e-06, 1.35442211e-06, 1.37359296e-06, 1.39276382e-06, &
        1.41193467e-06, 1.43110553e-06, 1.45027638e-06, 1.46944724e-06, &
        1.48861809e-06, 1.50778894e-06, 1.52695980e-06, 1.54613065e-06, &
        1.56530151e-06, 1.58447236e-06, 1.60364322e-06, 1.62281407e-06, &
        1.64198492e-06, 1.66115578e-06, 1.68032663e-06, 1.69949749e-06, &
        1.71866834e-06, 1.73783920e-06, 1.75701005e-06, 1.77618090e-06, &
        1.79535176e-06, 1.81452261e-06, 1.83369347e-06, 1.85286432e-06, &
        1.87203518e-06, 1.89120603e-06, 1.91037688e-06, 1.92954774e-06, &
        1.94871859e-06, 1.96788945e-06, 1.98706030e-06, 2.00623116e-06, &
        2.02540201e-06, 2.04457286e-06, 2.06374372e-06, 2.08291457e-06, &
        2.10208543e-06, 2.12125628e-06, 2.14042714e-06, 2.15959799e-06, &
        2.17876884e-06, 2.19793970e-06, 2.21711055e-06, 2.23628141e-06, &
        2.25545226e-06, 2.27462312e-06, 2.29379397e-06, 2.31296482e-06, &
        2.33213568e-06, 2.35130653e-06, 2.37047739e-06, 2.38964824e-06, &
        2.40881910e-06, 2.42798995e-06, 2.44716080e-06, 2.46633166e-06, &
        2.48550251e-06, 2.50467337e-06, 2.52384422e-06, 2.54301508e-06, &
        2.56218593e-06, 2.58135678e-06, 2.60052764e-06, 2.61969849e-06, &
        2.63886935e-06, 2.65804020e-06, 2.67721106e-06, 2.69638191e-06, &
        2.71555276e-06, 2.73472362e-06, 2.75389447e-06, 2.77306533e-06, &
        2.79223618e-06, 2.81140704e-06, 2.83057789e-06, 2.84974874e-06, &
        2.86891960e-06, 2.88809045e-06, 2.90726131e-06, 2.92643216e-06, &
        2.94560302e-06, 2.96477387e-06, 2.98394472e-06, 3.00311558e-06, &
        3.02228643e-06, 3.04145729e-06, 3.06062814e-06, 3.07979899e-06, &
        3.09896985e-06, 3.11814070e-06, 3.13731156e-06, 3.15648241e-06, &
        3.17565327e-06, 3.19482412e-06, 3.21399497e-06, 3.23316583e-06, &
        3.25233668e-06, 3.27150754e-06, 3.29067839e-06, 3.30984925e-06, &
        3.32902010e-06, 3.34819095e-06, 3.36736181e-06, 3.38653266e-06, &
        3.40570352e-06, 3.42487437e-06, 3.44404523e-06, 3.46321608e-06, &
        3.48238693e-06, 3.50155779e-06, 3.52072864e-06, 3.53989950e-06, &
        3.55907035e-06, 3.57824121e-06, 3.59741206e-06, 3.61658291e-06, &
        3.63575377e-06, 3.65492462e-06, 3.67409548e-06, 3.69326633e-06, &
        3.71243719e-06, 3.73160804e-06, 3.75077889e-06, 3.76994975e-06, &
        3.78912060e-06, 3.80829146e-06, 3.82746231e-06, 3.84663317e-06, &
        3.86580402e-06, 3.88497487e-06, 3.90414573e-06, 3.92331658e-06, &
        3.94248744e-06, 3.96165829e-06, 3.98082915e-06, 4.00000000e-06/


  DATA ni /7.74508407e-10, 7.74508407e-10, &
        7.74508407e-10, 7.74508407e-10, &
        7.74508407e-10, 7.74508407e-10, 7.74508407e-10, 7.74508407e-10, &
        6.98381122e-10, 6.23170274e-10, 5.97655992e-10, 5.84106004e-10, &
        5.44327597e-10, 5.71923510e-10, 6.59723827e-10, 8.05183870e-10, &
        1.03110161e-09, 1.36680386e-09, 1.85161253e-09, 2.56487751e-09, &
        3.56462855e-09, 4.89450926e-09, 6.49252022e-09, 9.62029335e-09, &
        1.32335015e-08, 1.75502184e-08, 2.19240625e-08, 3.03304156e-08, &
        4.07715972e-08, 5.00414911e-08, 7.09722331e-08, 1.00773751e-07, &
        1.31427409e-07, 1.42289041e-07, 1.49066787e-07, 2.01558515e-07, &
        2.99106105e-07, 4.03902086e-07, 4.54292169e-07, 5.21906983e-07, &
        6.27643362e-07, 9.43955678e-07, 1.33464494e-06, 1.97278315e-06, &
        2.31801329e-06, 2.20584676e-06, 1.85568138e-06, 1.73395193e-06, &
        1.73101406e-06, 1.91333936e-06, 2.26413019e-06, 3.23959718e-06, &
        4.94316963e-06, 6.89378896e-06, 1.02237444e-05, 1.21439656e-05, &
        1.31567585e-05, 1.33448288e-05, 1.32000000e-05, 1.31608040e-05, &
        1.33048369e-05, 1.40321464e-05, 1.51526244e-05, 1.80342858e-05, &
        3.82875736e-05, 1.07325259e-04, 2.11961637e-04, 3.82008054e-04, &
        5.30897470e-04, 5.29244735e-04, 4.90876605e-04, 4.33905427e-04, &
        3.77795349e-04, 3.17633099e-04, 2.81078564e-04, 2.57579485e-04, &
        2.42203100e-04, 2.23789060e-04, 2.04306870e-04, 1.87909255e-04, &
        1.73117146e-04, 1.61533186e-04, 1.53420328e-04, 1.47578033e-04, &
        1.42334776e-04, 1.35691466e-04, 1.30495414e-04, 1.36065123e-04, &
        1.70928821e-04, 2.66389730e-04, 4.80957955e-04, 8.25041961e-04, &
        1.21654792e-03, 1.50232875e-03, 1.62316078e-03, 1.61649750e-03, &
        1.53736801e-03, 1.42343711e-03, 1.24459117e-03, 1.02388611e-03, &
        7.89112523e-04, 5.97204264e-04, 4.57152413e-04, 3.62341259e-04, &
        2.99128332e-04, 2.57035569e-04, 2.26992203e-04, 2.07110355e-04, &
        2.05835688e-04, 2.25108810e-04, 2.64262893e-04, 3.23594011e-04, &
        3.93061117e-04, 4.62789970e-04, 5.19664416e-04, 5.59739628e-04, &
        5.93476084e-04, 6.22797885e-04, 6.57484833e-04, 6.92849600e-04, &
        7.26584901e-04, 7.56604648e-04, 7.68009488e-04, 7.65579073e-04, &
        7.50526164e-04, 7.39809972e-04, 7.55622847e-04, 8.05099514e-04, &
        9.67279246e-04, 1.16281559e-03, 1.42570247e-03, 2.04986585e-03, &
        2.93971170e-03, 4.49827711e-03, 7.32537532e-03, 1.18889734e-02, &
        1.85851805e-02, 2.86242532e-02, 4.34131035e-02, 6.37828307e-02, &
        9.24145850e-02, 1.35547945e-01, 1.94143245e-01, 2.54542814e-01, &
        3.02282024e-01, 3.42214181e-01, 3.85475065e-01, 4.38000000e-01, &
        4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, &
        4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, &
        4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, &
        4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, &
        4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, &
        4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, &
        4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, &
        4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, &
        4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, &
        4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, &
        4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, &
        4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, &
        4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01/


  Pas     =(lmax-lmin)/nlambda
  ropt    = dopt/2.0
  SSA     = 3.0/(rhoIce*ropt)
  cos30   = cos(30.0/180.0*pi)


  albint=0.
  norm=0.

  DO i = 1,nlambda
      gamm = 4.0 * pi * ni(i) / wavelengths(i)
      cosalb = 2.0 / (SSA * rhoice) * B_amp * gamm
      alpha = 16. / 3 * cosalb / (1.0 - g_asy)
      R = exp(-(alpha**0.5) * 3.0 / 7.0 * (1.0 + 2.0 * cossza))
      alb30 = exp(-(alpha**0.5)* 3.0 / 7.0 * (1.0 + 20 * cos30))

      IF ((wavelengths(i)>=lambdamin).AND. &
            (wavelengths(i)<lambdamax)) THEN
      albint=albint+R*Pas  ! rectangle integration -> can be improved with trapezintegration
      norm=norm+Pas
      ENDIF

  END DO

  albint=max(0.,min(albint/max(norm,1E-10),albmax))

END SUBROUTINE albedo_kokhanovsky