Index: LMDZ6/trunk/libf/phylmd/inlandsis/sisvat_sop.F
===================================================================
--- LMDZ6/trunk/libf/phylmd/inlandsis/sisvat_sop.F	(revision 3792)
+++ LMDZ6/trunk/libf/phylmd/inlandsis/sisvat_sop.F	(revision 3900)
@@ -1,4 +1,2 @@
- 
- 
       subroutine SnOptP(jjtime)
  
@@ -30,4 +28,5 @@
 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:                                                  |
@@ -68,5 +67,5 @@
       use VARySV
       use VARtSV
-      USE surface_data, only: iflag_albzenith
+      USE surface_data, only: iflag_albcalc,correc_alb
 
       IMPLICIT NONE
@@ -89,9 +88,9 @@
 c #AG real      agesno
  
-      integer   isn   ,ikl   ,isn1  
+      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
+      real      bsegal,czemax,csegal,csza
       real      RoFrez,DiffRo,SignRo,SnowOK,OpSqrt
       real      albSn1,albIc1,a_SnI1,a_SII1
@@ -103,7 +102,7 @@
       real      albedo_old,albCor
       real      ro_ave,dz_ave,minalb
-      real      thetazs,thetazs0,aap,bbp,ccp,alb0      ! parameters for zenoth angle effect at Dome C
-
-
+      real      l1min,l1max,l2min,l2max,l3min,l3max
+      real      l6min(6), l6max(6), albSn6(6), a_SII6(6)
+      real      lmintmp,lmaxtmp,albtmp
  
 C +--Local   DATA
@@ -138,8 +137,14 @@
       data      albmax  /0.99       /  ! Albedo max
 
-C +-- For the computation of solar zenoth angle effect from Dome C data
-C +  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-      data      aap/0.00016/,bbp/-0.017/,ccp/1.2/
-
+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/
  
  
@@ -177,7 +182,13 @@
      .                                       *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
@@ -191,5 +202,5 @@
         DO ikl=1,knonv
  
-           isn   =  max(iun,isnoSV(ikl))
+          isn   =  max(iun,isnoSV(ikl))
  
           SignRo = sign(unun, rocdSV - ro__SV(ikl,isn))
@@ -200,8 +211,51 @@
 cCA    +--Correction of snow albedo for Antarctica/Greenland
 cCA       --------------------------------------------------
-          albCor = 1.
+
+          
+          albCor = correc_alb 
 c #GL     albCor = 1.01
-c #AC     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)
@@ -219,8 +273,14 @@
           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
@@ -237,4 +297,6 @@
           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
@@ -304,47 +366,15 @@
           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.
  
- 
-!    Zenith Angle Correction (Segal et al.,         1991, JAS 48, p.1025)
-!    ----------------------- (Wiscombe & Warren, dec1980, JAS   , p.2723)
-!                            (Warren,               1982,  RG   , p.  81)
-!                            (Vignon, PhD Thesis, pp 150, conditions at Dome C)
-!                            -------------------------------------------------
- 
-         
-          dalbed = 0.0
-
-          if (iflag_albzenith .GE. 0) then
-          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))
-
-
-! Vignon PhD thesis, Dome C conditions
-
-            if (iflag_albzenith .EQ. 1) then
-            thetazs0=-bbp/(2.0*aap)
-            alb0=(bbp**2)/4./aap-(bbp**2)/2./aap+ccp
-            thetazs=max(acos(coszSV(ikl))/pi*180.0,thetazs0) ! in degrees
-
-
-            dalbeW = aap*(thetazs**2)+bbp*thetazs+ccp-alb0       
-            dalbed =     dalbeW      *       min(1,isnoSV(ikl))
-            end if
-
-
-            end if
+
  
 C +--Elsewhere    Integrated Snow/Ice Albedo
@@ -368,8 +398,13 @@
             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
@@ -382,4 +417,6 @@
             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))
@@ -397,7 +434,11 @@
             alb3sv(ikl) = alb3sv(ikl) - 1.0/3.0             ! 33 %
      .                  * (albedo_old-albisv(ikl)) / So3dSV
- 
- 
-C +--Integrated Snow/Ice/Soil Albedo: Maximum albedo = 95%
+            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 +  -----------------------------------------------------
  
@@ -410,4 +451,8 @@
             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)
 
 
@@ -436,5 +481,8 @@
             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
  
@@ -519,3 +567,199 @@
  
       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
+  
