Index: LMDZ6/trunk/libf/phylmd/inlandsis/sisvat_sno_albedo.F
===================================================================
--- LMDZ6/trunk/libf/phylmd/inlandsis/sisvat_sno_albedo.F	(revision 3980)
+++ LMDZ6/trunk/libf/phylmd/inlandsis/sisvat_sno_albedo.F	(revision 3980)
@@ -0,0 +1,765 @@
+      subroutine SnOptP(jjtime)
+ 
+C +------------------------------------------------------------------------+
+C | MAR/SISVAT   SnOptP                                    12-08-2019  MAR |
+C |   SubRoutine SnOptP computes the Snow Pack optical Properties          |
+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 |   ^^^^^    ispiSV   = 0,...,nsno: Uppermost Superimposed Ice Layer     |
+C |                                                                        |
+C |                                                                        |
+C |   INPUT:   G1snSV   : Dendricity (<0) or Sphericity (>0) of Snow Layer |
+C |   ^^^^^    G2snSV   : Sphericity (>0) or Size            of Snow Layer |
+C |            agsnSV   : Snow       Age                             [day] |
+C |            ro__SV   : Snow/Soil  Volumic Mass                  [kg/m3] |
+C |            eta_SV   : Water      Content                       [m3/m3] |
+C |            rusnSV   : Surficial  Water   Thickness   [kg/m2] .OR. [mm] |
+C |            SWS_SV   : Surficial  Water   Status                        |
+C |            dzsnSV   : Snow       Layer   Thickness                 [m] |
+C |                                                                        |
+C |            albssv   : Soil       Albedo                            [-] |
+C |            zzsnsv   : Snow       Pack    Thickness                 [m] |
+C |                                                                        |
+C |   OUTPUT:  albisv   : Snow/Ice/Water/Soil Integrated Albedo        [-] |
+C |   ^^^^^^   sEX_sv   : Verticaly Integrated Extinction Coefficient      |
+C |            DOPsnSV   : Snow Optical diameter [m]                       |
+C |                                                                        |
+C |   Internal Variables:                                                  |
+C |   ^^^^^^^^^^^^^^^^^^                                                   |
+C |            SnOpSV   : Snow Grain optical Size                      [m] |
+C |            EX1_sv   : Integrated Snow Extinction (0.3--0.8micr.m)      |
+C |            EX2_sv   : Integrated Snow Extinction (0.8--1.5micr.m)      |
+C |            EX3_sv   : Integrated Snow Extinction (1.5--2.8micr.m)      |
+C |                                                                        |
+C |   METHODE:    Calcul de la taille optique des grains ? partir de       |
+C |   ^^^^^^^    -leur type decrit par les deux variables descriptives     |
+C |                    continues sur la plage -99/+99 passees en appel.    |
+C |              -la taille optique (1/10mm) des etoiles,                  |
+C |                                          des grains fins et            |
+C |                                          des jeunes faces planes       |
+C |                                                                        |
+C |   METHOD:     Computation of the optical diameter of the grains        |
+C |   ^^^^^^      described with the CROCUS formalism G1snSV / G2snSV      |
+C |                                                                        |
+C |   REFERENCE: Brun et al.      1989, J. Glaciol 35 pp. 333--342         |
+C |   ^^^^^^^^^  Brun et al.      1992, J. Glaciol 38 pp.  13-- 22         |
+C |              Eric Martin Sept.1996                                     |
+C |                                                                        |
+C |                                                                        |
+C +------------------------------------------------------------------------+
+ 
+ 
+ 
+ 
+C +--Global Variables
+C +  ================
+ 
+
+      use VARphy
+      use VAR_SV
+      use VARdSV
+      use VARxSV
+      use VARySV
+      use VARtSV
+      USE surface_data, only: iflag_albcalc,correc_alb
+
+      IMPLICIT NONE
+
+ 
+C + -- INPUT
+      integer jjtime
+
+C +--Internal Variables
+C +  ==================
+ 
+      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
+c #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
+ 
+C +--Local   DATA
+C +  ============
+ 
+C +--For the computation of the solar irradiance extinction in snow
+C +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      data      sbeta1/0.0192/,sbeta2/0.4000/,sbeta3/0.1098/
+      data      sbeta4/1.0000/
+      data      sbeta5/2.00e1/
+ 
+C +--Snow Age Maximum (Taiga, e.g. Col de Porte)
+C +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      data      AgeMax  /60.0/
+C +...          AgeMax:  Snow Age Maximum                              [day]
+ 
+      data      AlbMin  /0.94/
+C +...          AlbMin:  Albedo   Minimum / visible (0.3--0.8 micrometers)
+ 
+      data      HSnoSV  /0.01/
+C +...          HSnoSV:  Snow     Thickness through witch
+C +                      Albedo is interpolated to Ice  Albedo
+      data      HIceSV  /0.10/
+C +...          HIceSV:  Snow/Ice Thickness through witch
+C +                      Albedo is interpolated to Soil Albedo
+      data      doptmx  /2.3e-3/
+C +...          doptmx:  Maximum  optical Diameter    (pi * R**2)        [m]
+C +
+      data      czeMAX  /0.173648178/  ! 80.deg (Segal et al., 1991 JAS)
+
+      data      bsegal  /4.00       /  !
+      data      albmax  /0.99       /  ! Albedo max
+
+C +-- 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/
+ 
+ 
+C +--Snow Grain optical Size
+C +  =======================
+ 
+        DO ikl=1,knonv
+         DO   isn=1,max(1,isnoSV(ikl))
+ 
+          G2snSV(ikl,isn) =  max(epsi,G2snSV(ikl,isn))
+C +...    Avoid non physical Values
+ 
+          SignG1          = sign(unun,G1snSV(ikl,isn))
+          Sph_OK          =  max(zero,SignG1)
+ 
+          SnOpSV(ikl,isn) =   1.e-4 *
+C +...    SI:           (from 1/10 mm to m)
+ 
+ 
+C +--Contribution of Non Dendritic Snow
+C +  ----------------------------------
+ 
+     .    (    Sph_OK *(      G2snSV(ikl,isn)*G1snSV(ikl,isn)/G1_dSV
+     .              +max(demi*G2snSV(ikl,isn),DFcdSV)
+     .                 *(unun-G1snSV(ikl,isn)                /G1_dSV))
+ 
+ 
+C +--Contribution of     Dendritic Snow
+C +  ----------------------------------
+ 
+     .    +(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))
+
+C + --For outputs (Etienne)
+C + ------------------------
+          DOPsnSV(ikl,isn)=SnOpSV(ikl,isn)
+        END DO
+      END DO
+ 
+
+
+ 
+C +--Snow/Ice Albedo
+C +  ===============
+
+ 
+ 
+C +--Uppermost effective Snow Layer
+C +  ------------------------------
+ 
+        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))
+ 
+cCA    +--Correction of snow albedo for Antarctica/Greenland
+cCA       --------------------------------------------------
+
+          
+          albCor = correc_alb 
+c #GL     albCor = 1.01
+c #AC    albCor = 1.01
+
+
+          IF (iflag_albcalc .GE. 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))
+c #cz     dalbeS =   ((bsegal+unun)/(unun+2.0*bsegal*csegal)
+c #cz.                -       unun                          )*0.32
+c #cz.              /  bsegal
+c #cz     dalbeS = max(dalbeS,zero)
+c #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
+c #GL     albSn1 =  albSn1*max(0.9,(1.-1.5*eta_SV(ikl,isn)))
+c #GL     albSn2 =  albSn2*max(0.9,(1.-1.5*eta_SV(ikl,isn)))
+c #GL     albSn3 =  albSn3*max(0.9,(1.-1.5*eta_SV(ikl,isn)))
+
+          ENDIF
+
+ 
+          albSno =  So1dSV*albSn1
+     .           +  So2dSV*albSn2
+     .           +  So3dSV*albSn3
+ 
+cXF
+          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)
+
+ 
+c +           ro < roSdSV |          min al > aI3dSV
+c +  roSdSV < ro < rocdSV | aI2dSV < min al < aI3dSV (fct of density)
+ 
+ 
+C +--Snow/Ice Pack Thickness
+C +  -----------------------
+ 
+          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)
+ 
+C +       The value of SnownH is set to 1 in case of ice lenses above
+C +       1m of dry snow (ro<600kg/m3) for using CROCUS albedo
+ 
+c         ro_ave =  0.
+c         dz_ave =  0.
+c         SnowOK =  1.
+c      do isn    =  isnoSV(ikl),1,-1
+c         ro_ave =  ro_ave + ro__SV(ikl,isn) * dzsnSV(ikl,isn) * SnowOK
+c         dz_ave =  dz_ave +                   dzsnSV(ikl,isn) * SnowOK
+c         SnowOK =  max(zero,sign(unun,1.-dz_ave))
+c      enddo
+ 
+c         ro_ave =  ro_ave / max(dz_ave,epsi)
+c         SnowOK =  max(zero,sign(unun,600.-ro_ave))
+c         SnownH =  SnowOK + SnownH * (1. - SnowOK)
+ 
+C +--Integrated Snow/Ice Albedo: Case of Water on Bare Ice
+C +  -----------------------------------------------------
+ 
+          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)                        !
+c         albWIc = max(aI1dSV,min(aI2dSV,albWIc+slopSV(ikl)*
+c    .             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))
+ 
+c +  rocdSV < ro < ro_ice | aI2dSV< al <aI3dSV (fct of density)
+c +           ro > ro_ice | aI1dSV< al <aI2dSV (fct of superficial water content
+ 
+ 
+C +--Integrated Snow/Ice      Albedo
+C +  -------------------------------
+ 
+          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
+
+cc #AG     agesno =      min(agsnSV(ikl,isn)          ,AgeMax)
+cc #AG     a_SII1      =     a_SII1      -0.175*agesno/AgeMax
+C +...                                   Impurities: Col de Porte Parameter.
+ 
+
+ 
+C +--Elsewhere    Integrated Snow/Ice Albedo
+C +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+c #cp     ELSE
+            albSII =     So1dSV*a_SII1
+     .                 + So2dSV*a_SII2
+     .                 + So3dSV*a_SII3
+c #cp     END IF
+ 
+ 
+C +--Integrated Snow/Ice/Soil Albedo
+C +  -------------------------------
+ 
+            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
+
+ 
+C +--Integrated Snow/Ice/Soil Albedo: Clouds Correction! Greuell & all., 1994
+C +  --------------------------------------------------! 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))
+ 
+C +--Integrated Snow/Ice/Soil Albedo: Minimum snow albedo = aI1dSV
+C +  -------------------------------------------------------------
+ 
+            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)
+
+
+C +--Integrated Snow/Ice/Soil Albedo: Maximum albedo = 95% 
+C +  -----------------------------------------------------
+ 
+            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 
+	
+c #AO      if (LSmask(ikl) .eq. 0 .and. coupling_ao .eq. .true.) then 
+c #AO       if (AOmask(ikl) .eq. 0) then  
+c #AO       albisv(ikl) =  (1.-AOmask(ikl))* albAOsisv(ikl)
+c #AO.                    +(AOmask(ikl)*albisv(ikl))
+c #AO       alb1sv(ikl) =   (1.-AOmask(ikl))* albAOsisv(ikl)
+c #AO.                    +(AOmask(ikl)*alb1sv(ikl))
+c #AO       alb2sv(ikl) =   (1.-AOmask(ikl))* albAOsisv(ikl)
+c #AO.                    +(AOmask(ikl)*alb2sv(ikl))
+c #AO       alb3sv(ikl) =   (1.-AOmask(ikl))* albAOsisv(ikl)
+c #AO.                    +(AOmask(ikl)*alb3sv(ikl))
+c #AO       endif
+c #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
+ 
+ 
+C +--Extinction Coefficient: Exponential Factor
+C +  ==========================================
+ 
+        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
+ 
+cXF
+ 
+        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
+ 
+ 
+C +--Integrated Extinction of Solar Irradiance (Normalized Value)
+C +  ============================================================
+ 
+          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
+ 
+ 
+      return
+
+
+      end
+
+
+!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+      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).GE.lambdamin).AND.
+     .       (wavelengths(i).LT.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
+  
Index: LMDZ6/trunk/libf/phylmd/inlandsis/sisvat_sop.F
===================================================================
--- LMDZ6/trunk/libf/phylmd/inlandsis/sisvat_sop.F	(revision 3978)
+++ 	(revision )
@@ -1,765 +1,0 @@
-      subroutine SnOptP(jjtime)
- 
-C +------------------------------------------------------------------------+
-C | MAR/SISVAT   SnOptP                                    12-08-2019  MAR |
-C |   SubRoutine SnOptP computes the Snow Pack optical Properties          |
-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 |   ^^^^^    ispiSV   = 0,...,nsno: Uppermost Superimposed Ice Layer     |
-C |                                                                        |
-C |                                                                        |
-C |   INPUT:   G1snSV   : Dendricity (<0) or Sphericity (>0) of Snow Layer |
-C |   ^^^^^    G2snSV   : Sphericity (>0) or Size            of Snow Layer |
-C |            agsnSV   : Snow       Age                             [day] |
-C |            ro__SV   : Snow/Soil  Volumic Mass                  [kg/m3] |
-C |            eta_SV   : Water      Content                       [m3/m3] |
-C |            rusnSV   : Surficial  Water   Thickness   [kg/m2] .OR. [mm] |
-C |            SWS_SV   : Surficial  Water   Status                        |
-C |            dzsnSV   : Snow       Layer   Thickness                 [m] |
-C |                                                                        |
-C |            albssv   : Soil       Albedo                            [-] |
-C |            zzsnsv   : Snow       Pack    Thickness                 [m] |
-C |                                                                        |
-C |   OUTPUT:  albisv   : Snow/Ice/Water/Soil Integrated Albedo        [-] |
-C |   ^^^^^^   sEX_sv   : Verticaly Integrated Extinction Coefficient      |
-C |            DOPsnSV   : Snow Optical diameter [m]                       |
-C |                                                                        |
-C |   Internal Variables:                                                  |
-C |   ^^^^^^^^^^^^^^^^^^                                                   |
-C |            SnOpSV   : Snow Grain optical Size                      [m] |
-C |            EX1_sv   : Integrated Snow Extinction (0.3--0.8micr.m)      |
-C |            EX2_sv   : Integrated Snow Extinction (0.8--1.5micr.m)      |
-C |            EX3_sv   : Integrated Snow Extinction (1.5--2.8micr.m)      |
-C |                                                                        |
-C |   METHODE:    Calcul de la taille optique des grains ? partir de       |
-C |   ^^^^^^^    -leur type decrit par les deux variables descriptives     |
-C |                    continues sur la plage -99/+99 passees en appel.    |
-C |              -la taille optique (1/10mm) des etoiles,                  |
-C |                                          des grains fins et            |
-C |                                          des jeunes faces planes       |
-C |                                                                        |
-C |   METHOD:     Computation of the optical diameter of the grains        |
-C |   ^^^^^^      described with the CROCUS formalism G1snSV / G2snSV      |
-C |                                                                        |
-C |   REFERENCE: Brun et al.      1989, J. Glaciol 35 pp. 333--342         |
-C |   ^^^^^^^^^  Brun et al.      1992, J. Glaciol 38 pp.  13-- 22         |
-C |              Eric Martin Sept.1996                                     |
-C |                                                                        |
-C |                                                                        |
-C +------------------------------------------------------------------------+
- 
- 
- 
- 
-C +--Global Variables
-C +  ================
- 
-
-      use VARphy
-      use VAR_SV
-      use VARdSV
-      use VARxSV
-      use VARySV
-      use VARtSV
-      USE surface_data, only: iflag_albcalc,correc_alb
-
-      IMPLICIT NONE
-
- 
-C + -- INPUT
-      integer jjtime
-
-C +--Internal Variables
-C +  ==================
- 
-      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
-c #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
- 
-C +--Local   DATA
-C +  ============
- 
-C +--For the computation of the solar irradiance extinction in snow
-C +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-      data      sbeta1/0.0192/,sbeta2/0.4000/,sbeta3/0.1098/
-      data      sbeta4/1.0000/
-      data      sbeta5/2.00e1/
- 
-C +--Snow Age Maximum (Taiga, e.g. Col de Porte)
-C +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-      data      AgeMax  /60.0/
-C +...          AgeMax:  Snow Age Maximum                              [day]
- 
-      data      AlbMin  /0.94/
-C +...          AlbMin:  Albedo   Minimum / visible (0.3--0.8 micrometers)
- 
-      data      HSnoSV  /0.01/
-C +...          HSnoSV:  Snow     Thickness through witch
-C +                      Albedo is interpolated to Ice  Albedo
-      data      HIceSV  /0.10/
-C +...          HIceSV:  Snow/Ice Thickness through witch
-C +                      Albedo is interpolated to Soil Albedo
-      data      doptmx  /2.3e-3/
-C +...          doptmx:  Maximum  optical Diameter    (pi * R**2)        [m]
-C +
-      data      czeMAX  /0.173648178/  ! 80.deg (Segal et al., 1991 JAS)
-
-      data      bsegal  /4.00       /  !
-      data      albmax  /0.99       /  ! Albedo max
-
-C +-- 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/
- 
- 
-C +--Snow Grain optical Size
-C +  =======================
- 
-        DO ikl=1,knonv
-         DO   isn=1,max(1,isnoSV(ikl))
- 
-          G2snSV(ikl,isn) =  max(epsi,G2snSV(ikl,isn))
-C +...    Avoid non physical Values
- 
-          SignG1          = sign(unun,G1snSV(ikl,isn))
-          Sph_OK          =  max(zero,SignG1)
- 
-          SnOpSV(ikl,isn) =   1.e-4 *
-C +...    SI:           (from 1/10 mm to m)
- 
- 
-C +--Contribution of Non Dendritic Snow
-C +  ----------------------------------
- 
-     .    (    Sph_OK *(      G2snSV(ikl,isn)*G1snSV(ikl,isn)/G1_dSV
-     .              +max(demi*G2snSV(ikl,isn),DFcdSV)
-     .                 *(unun-G1snSV(ikl,isn)                /G1_dSV))
- 
- 
-C +--Contribution of     Dendritic Snow
-C +  ----------------------------------
- 
-     .    +(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))
-
-C + --For outputs (Etienne)
-C + ------------------------
-          DOPsnSV(ikl,isn)=SnOpSV(ikl,isn)
-        END DO
-      END DO
- 
-
-
- 
-C +--Snow/Ice Albedo
-C +  ===============
-
- 
- 
-C +--Uppermost effective Snow Layer
-C +  ------------------------------
- 
-        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))
- 
-cCA    +--Correction of snow albedo for Antarctica/Greenland
-cCA       --------------------------------------------------
-
-          
-          albCor = correc_alb 
-c #GL     albCor = 1.01
-c #AC    albCor = 1.01
-
-
-          IF (iflag_albcalc .GE. 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))
-c #cz     dalbeS =   ((bsegal+unun)/(unun+2.0*bsegal*csegal)
-c #cz.                -       unun                          )*0.32
-c #cz.              /  bsegal
-c #cz     dalbeS = max(dalbeS,zero)
-c #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
-c #GL     albSn1 =  albSn1*max(0.9,(1.-1.5*eta_SV(ikl,isn)))
-c #GL     albSn2 =  albSn2*max(0.9,(1.-1.5*eta_SV(ikl,isn)))
-c #GL     albSn3 =  albSn3*max(0.9,(1.-1.5*eta_SV(ikl,isn)))
-
-          ENDIF
-
- 
-          albSno =  So1dSV*albSn1
-     .           +  So2dSV*albSn2
-     .           +  So3dSV*albSn3
- 
-cXF
-          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)
-
- 
-c +           ro < roSdSV |          min al > aI3dSV
-c +  roSdSV < ro < rocdSV | aI2dSV < min al < aI3dSV (fct of density)
- 
- 
-C +--Snow/Ice Pack Thickness
-C +  -----------------------
- 
-          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)
- 
-C +       The value of SnownH is set to 1 in case of ice lenses above
-C +       1m of dry snow (ro<600kg/m3) for using CROCUS albedo
- 
-c         ro_ave =  0.
-c         dz_ave =  0.
-c         SnowOK =  1.
-c      do isn    =  isnoSV(ikl),1,-1
-c         ro_ave =  ro_ave + ro__SV(ikl,isn) * dzsnSV(ikl,isn) * SnowOK
-c         dz_ave =  dz_ave +                   dzsnSV(ikl,isn) * SnowOK
-c         SnowOK =  max(zero,sign(unun,1.-dz_ave))
-c      enddo
- 
-c         ro_ave =  ro_ave / max(dz_ave,epsi)
-c         SnowOK =  max(zero,sign(unun,600.-ro_ave))
-c         SnownH =  SnowOK + SnownH * (1. - SnowOK)
- 
-C +--Integrated Snow/Ice Albedo: Case of Water on Bare Ice
-C +  -----------------------------------------------------
- 
-          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)                        !
-c         albWIc = max(aI1dSV,min(aI2dSV,albWIc+slopSV(ikl)*
-c    .             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))
- 
-c +  rocdSV < ro < ro_ice | aI2dSV< al <aI3dSV (fct of density)
-c +           ro > ro_ice | aI1dSV< al <aI2dSV (fct of superficial water content
- 
- 
-C +--Integrated Snow/Ice      Albedo
-C +  -------------------------------
- 
-          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
-
-cc #AG     agesno =      min(agsnSV(ikl,isn)          ,AgeMax)
-cc #AG     a_SII1      =     a_SII1      -0.175*agesno/AgeMax
-C +...                                   Impurities: Col de Porte Parameter.
- 
-
- 
-C +--Elsewhere    Integrated Snow/Ice Albedo
-C +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-c #cp     ELSE
-            albSII =     So1dSV*a_SII1
-     .                 + So2dSV*a_SII2
-     .                 + So3dSV*a_SII3
-c #cp     END IF
- 
- 
-C +--Integrated Snow/Ice/Soil Albedo
-C +  -------------------------------
- 
-            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
-
- 
-C +--Integrated Snow/Ice/Soil Albedo: Clouds Correction! Greuell & all., 1994
-C +  --------------------------------------------------! 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))
- 
-C +--Integrated Snow/Ice/Soil Albedo: Minimum snow albedo = aI1dSV
-C +  -------------------------------------------------------------
- 
-            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)
-
-
-C +--Integrated Snow/Ice/Soil Albedo: Maximum albedo = 95% 
-C +  -----------------------------------------------------
- 
-            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 
-	
-c #AO      if (LSmask(ikl) .eq. 0 .and. coupling_ao .eq. .true.) then 
-c #AO       if (AOmask(ikl) .eq. 0) then  
-c #AO       albisv(ikl) =  (1.-AOmask(ikl))* albAOsisv(ikl)
-c #AO.                    +(AOmask(ikl)*albisv(ikl))
-c #AO       alb1sv(ikl) =   (1.-AOmask(ikl))* albAOsisv(ikl)
-c #AO.                    +(AOmask(ikl)*alb1sv(ikl))
-c #AO       alb2sv(ikl) =   (1.-AOmask(ikl))* albAOsisv(ikl)
-c #AO.                    +(AOmask(ikl)*alb2sv(ikl))
-c #AO       alb3sv(ikl) =   (1.-AOmask(ikl))* albAOsisv(ikl)
-c #AO.                    +(AOmask(ikl)*alb3sv(ikl))
-c #AO       endif
-c #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
- 
- 
-C +--Extinction Coefficient: Exponential Factor
-C +  ==========================================
- 
-        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
- 
-cXF
- 
-        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
- 
- 
-C +--Integrated Extinction of Solar Irradiance (Normalized Value)
-C +  ============================================================
- 
-          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
- 
- 
-      return
-
-
-      end
-
-
-!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-      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).GE.lambdamin).AND.
-     .       (wavelengths(i).LT.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
-