Index: /trunk/LMDZ.GENERIC/README
===================================================================
--- /trunk/LMDZ.GENERIC/README	(revision 1481)
+++ /trunk/LMDZ.GENERIC/README	(revision 1482)
@@ -1085,2 +1085,8 @@
 == 21/09/2015 == MT
 - Cleanup of physiq.F90 and changed condense_cloud.F90 to condense_co2.F90
+
+== 14/10/2015 == MT
+- Implementation of the Spectral Albedo. Albedo(ngrid) is now Albedo(ngrid,L_NSPECTV) and has a value for each visible band.
+- Albedodat/Albedo0 have been removed. We now use Albedo_bareground.
+- CO2 Ice Albedo North/South dichotomy has been removed. CO2 ice Albedo is no longer stocked in start file. You can now
+prescribe CO2 ice albedo in callphys.def with "albedoco2ice". Its default value is 0.5.
Index: /trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90	(revision 1481)
+++ /trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90	(revision 1482)
@@ -1,7 +1,8 @@
       subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf,           &
-          albedo,emis,mu0,pplev,pplay,pt,                      & 
+          albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt,    & 
           tsurf,fract,dist_star,aerosol,muvar,                 &
           dtlw,dtsw,fluxsurf_lw,                               &
-          fluxsurf_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn,        &
+          fluxsurf_sw,fluxsurfabs_sw,fluxtop_lw,               &
+          fluxabs_sw,fluxtop_dn,                               &
           OLR_nu,OSR_nu,                                       &
           tau_col,cloudfrac,totcloudfrac,                      &
@@ -49,5 +50,5 @@
       integer,intent(in) :: nq ! number of tracers
       REAL,INTENT(IN) :: qsurf(ngrid,nq) ! tracer on surface (kg.m-2)
-      REAL,INTENT(IN) :: albedo(ngrid)   ! SW albedo
+      REAL,INTENT(IN) :: albedo(ngrid,L_NSPECTV)   ! Spectral SW albedo.MT2015.
       REAL,INTENT(IN) :: emis(ngrid)     ! LW emissivity
       real,intent(in) :: mu0(ngrid) ! cosine of sun incident angle
@@ -64,4 +65,5 @@
       REAL,INTENT(OUT) :: fluxsurf_lw(ngrid)   ! incident LW flux to surf (W/m2)
       REAL,INTENT(OUT) :: fluxsurf_sw(ngrid)   ! incident SW flux to surf (W/m2)
+      REAL,INTENT(OUT) :: fluxsurfabs_sw(ngrid)   ! absorbed SW flux by the surface (W/m2). By MT2015.
       REAL,INTENT(OUT) :: fluxtop_lw(ngrid)    ! outgoing LW flux to space (W/m2)
       REAL,INTENT(OUT) :: fluxabs_sw(ngrid)    ! SW flux absorbed by planet (W/m2)
@@ -70,4 +72,5 @@
       REAL,INTENT(OUT) :: OSR_nu(ngrid,L_NSPECTV)! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1)
       REAL,INTENT(OUT) :: tau_col(ngrid) ! diagnostic from aeropacity
+      REAL,INTENT(OUT) :: albedo_equivalent(ngrid) ! Spectrally Integrated Albedo. For Diagnostic. By MT2015
 !     for H2O cloud fraction in aeropacity
       real,intent(in) :: cloudfrac(ngrid,nlayer)
@@ -121,5 +124,6 @@
       REAL*8 fluxupv(L_NLAYRAD),fluxupi(L_NLAYRAD)
       REAL*8 fluxdnv(L_NLAYRAD),fluxdni(L_NLAYRAD)
-      REAL*8 albi,albv,acosz
+      REAL*8 albi,acosz
+      REAL*8 albv(L_NSPECTV) ! Spectral Visible Albedo.
 
       INTEGER ig,l,k,nw,iaer,irad
@@ -182,4 +186,7 @@
       real vtmp(nlayer)
       REAL*8 muvarrad(L_LEVELS)
+      
+!     included by MT for albedo calculations.      
+      REAL*8 albedo_temp(L_NSPECTV) ! For equivalent albedo calculation.
 
 !===============================================================
@@ -212,6 +219,4 @@
          call su_aer_radii(ngrid,nlayer,reffrad,nueffrad)
 	 
-	 
-
 !--------------------------------------------------
 !     set up correlated k
@@ -448,12 +453,18 @@
 
 !     Albedo and emissivity
-         albi=1-emis(ig)        ! longwave
-         albv=albedo(ig)        ! shortwave 
-
-      if(nosurf.and.(albv.gt.0.0))then
-         print*,'For open lower boundary in callcorrk must'
-         print*,'have surface albedo set to zero!'
-         call abort
-      endif
+         albi=1-emis(ig)   ! Longwave.
+         DO nw=1,L_NSPECTV ! Shortwave loop.
+            albv(nw)=albedo(ig,nw)
+	 ENDDO
+
+         if (nosurf) then
+            DO nw=1,L_NSPECTV
+               if(albv(nw).gt.0.0) then
+	          print*,'For open lower boundary in callcorrk must'
+                  print*,'have spectral surface band albedos all set to zero!'
+                  call abort
+	       endif
+	    ENDDO         
+         endif
 
       if ((ngrid.eq.1).and.(global1d)) then       ! fixed zenith angle 'szangle' in 1D simulations w/ globally-averaged sunlight
@@ -764,4 +775,17 @@
          end if
 
+
+         ! Equivalent Albedo Calculation. MT2015
+         if(fract(ig) .ge. 1.0e-4) then ! equivalent albedo makes sense only during daylight.            
+	    DO nw=1,L_NSPECTV	          
+	       albedo_temp(nw)=albedo(ig,nw)*nfluxgndv_nu(nw)
+            ENDDO
+	    albedo_temp(1:L_NSPECTV)=albedo_temp(1:L_NSPECTV)/sum(nfluxgndv_nu(1:L_NSPECTV))
+            albedo_equivalent(ig)=sum(albedo_temp(1:L_NSPECTV))
+	 else
+	    albedo_equivalent(ig)=0.0 ! Spectrally Integrated Albedo not defined for non-irradiated grid points. So we arbitrary set the equivalent albedo to 0.
+	 endif
+
+
 !-----------------------------------------------------------------------
 !     Longwave
@@ -786,4 +810,7 @@
          fluxsurf_lw(ig) = real(fluxdni(L_NLAYRAD))
          fluxsurf_sw(ig) = real(fluxdnv(L_NLAYRAD))
+         
+!        Flux absorbed by the surface. By MT2015.         
+         fluxsurfabs_sw(ig) = fluxsurf_sw(ig)*(1.-albedo_equivalent(ig))
 
          if(fluxtop_dn(ig).lt.0.0)then
Index: /trunk/LMDZ.GENERIC/libf/phystd/callkeys_mod.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/callkeys_mod.F90	(revision 1481)
+++ /trunk/LMDZ.GENERIC/libf/phystd/callkeys_mod.F90	(revision 1482)
@@ -66,4 +66,5 @@
       real pceil
       real albedosnow
+      real albedoco2ice
       real maxicethick
       real Tsaldiff
Index: /trunk/LMDZ.GENERIC/libf/phystd/condense_co2.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/condense_co2.F90	(revision 1481)
+++ /trunk/LMDZ.GENERIC/libf/phystd/condense_co2.F90	(revision 1482)
@@ -1,14 +1,15 @@
       subroutine condense_co2(ngrid,nlayer,nq,ptimestep, &
-          pcapcal,pplay,pplev,ptsrf,pt, &
-          pphi,pdt,pdu,pdv,pdtsrf,pu,pv,pq,pdq, &
-          piceco2,psolaralb,pemisurf, &
-          pdtc,pdtsrfc,pdpsrf,pduc,pdvc, &
+          pcapcal,pplay,pplev,ptsrf,pt,                  &
+          pphi,pdt,pdu,pdv,pdtsrf,pu,pv,pq,pdq,          &
+          piceco2,albedo,pemisurf,                       &
+          albedo_bareground,albedo_co2_ice_SPECTV,       &
+          pdtc,pdtsrfc,pdpsrf,pduc,pdvc,                 &
           pdqc)
 
-      use radinc_h, only : naerkind
+      use radinc_h, only : L_NSPECTV, naerkind
       use gases_h, only: gfrac, igas_co2
       use radii_mod, only : co2_reffrad
       use aerosol_mod, only : iaero_co2
-      USE surfdat_h, only: albedodat, albedice, emisice, emissiv
+      USE surfdat_h, only: emisice, emissiv
       USE comgeomfi_h, only: lati
       USE tracer_h, only: noms, rho_co2
@@ -32,5 +33,5 @@
 !     ptsrf(ngrid)          Surface temperature
 !     
-!     pdt(ngrid,nlayer)   Time derivative before condensation/sublimation of pt
+!     pdt(ngrid,nlayer)     Time derivative before condensation/sublimation of pt
 !     pdtsrf(ngrid)         Time derivative before condensation/sublimation of ptsrf
 !     pqsurf(ngrid,nq)      Sedimentation flux at the surface (kg.m-2.s-1)
@@ -44,7 +45,7 @@
 !     Both
 !     ----
-!     piceco2(ngrid)        CO2 ice at the surface (kg/m2)
-!     psolaralb(ngrid)      Albedo at the surface
-!     pemisurf(ngrid)       Emissivity of the surface
+!     piceco2(ngrid)               CO2 ice at the surface (kg/m2)
+!     albedo(ngrid,L_NSPECTV)      Spectral albedo at the surface
+!     pemisurf(ngrid)              Emissivity of the surface
 !
 !     Authors
@@ -77,6 +78,8 @@
       REAL,INTENT(IN) :: pq(ngrid,nlayer,nq)
       REAL,INTENT(IN) :: pdq(ngrid,nlayer,nq)
+      REAL,INTENT(IN) :: albedo_bareground(ngrid)
+      REAL,INTENT(IN) :: albedo_co2_ice_SPECTV(L_NSPECTV)
       REAL,INTENT(INOUT) :: piceco2(ngrid)
-      REAL,INTENT(OUT) :: psolaralb(ngrid)
+      REAL,INTENT(OUT) :: albedo(ngrid,L_NSPECTV)
       REAL,INTENT(OUT) :: pemisurf(ngrid)
       REAL,INTENT(OUT) :: pdtc(ngrid,nlayer)
@@ -90,5 +93,5 @@
 !     Local variables
 
-      INTEGER l,ig,icap,ilay,it,iq
+      INTEGER l,ig,icap,ilay,it,iq,nw
 
       REAL reffrad(ngrid,nlayer) ! radius (m) of the co2 ice particles
@@ -461,9 +464,13 @@
             piceco2(ig)=0.
          endif
-         if (piceco2(ig).gt.0) then
-            psolaralb(ig) = albedice(icap)
+         if (piceco2(ig) .gt. 1.) then  ! CO2 Albedo condition changed to 1 mm coverage. Change by MT2015.
+	    DO nw=1,L_NSPECTV
+               albedo(ig,nw) = albedo_co2_ice_SPECTV(nw)
+	    ENDDO
             emisref(ig)   = emisice(icap)
          else
-            psolaralb(ig) = albedodat(ig)
+            DO nw=1,L_NSPECTV
+               albedo(ig,nw) = albedo_bareground(ig) ! Note : If you have some water, it will be taken into account in the "hydrol" routine.
+	    ENDDO
             emisref(ig)   = emissiv
             pemisurf(ig)  = emissiv
Index: /trunk/LMDZ.GENERIC/libf/phystd/hydrol.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/hydrol.F90	(revision 1481)
+++ /trunk/LMDZ.GENERIC/libf/phystd/hydrol.F90	(revision 1482)
@@ -1,5 +1,7 @@
-subroutine hydrol(ngrid,nq,ptimestep,rnat,tsurf,          &
+subroutine hydrol(ngrid,nq,ptimestep,rnat,tsurf,  &
      qsurf,dqsurf,dqs_hyd,pcapcal,                &
-     albedo0,albedo,mu0,pdtsurf,pdtsurf_hyd,hice, &
+     albedo,albedo_bareground,                    &
+     albedo_snow_SPECTV,albedo_co2_ice_SPECTV,    &
+     mu0,pdtsurf,pdtsurf_hyd,hice,                &
      pctsrf_sic,sea_ice)
 
@@ -12,5 +14,6 @@
   USE tracer_h
   use slab_ice_h
-  use callkeys_mod, only: albedosnow,ok_slab_ocean,Tsaldiff,maxicethick,co2cond
+  use callkeys_mod, only: albedosnow,albedoco2ice,ok_slab_ocean,Tsaldiff,maxicethick,co2cond
+  use radinc_h, only : L_NSPECTV
 
   implicit none
@@ -25,5 +28,6 @@
 !     ------- 
 !     Adapted from LMDTERRE by B. Charnay (2010). Further 
-!     modifications by R. Wordsworth (2010).
+!     Modifications by R. Wordsworth (2010).
+!     Spectral albedo by M. Turbet (2015).
 !     
 !     Called by
@@ -45,8 +49,4 @@
 !     Inputs
 !     ------
-      real albedoice
-      save albedoice
-!$OMP THREADPRIVATE(albedoice)
-
       real snowlayer
       parameter (snowlayer=33.0)        ! 33 kg/m^2 of snow, equal to a layer of 3.3 cm 
@@ -72,5 +72,8 @@
       real dqsurf(ngrid,nq), pdtsurf(ngrid)
       real hice(ngrid)
-      real albedo0(ngrid), albedo(ngrid)
+      real albedo(ngrid,L_NSPECTV)
+      real albedo_bareground(ngrid)
+      real albedo_snow_SPECTV(L_NSPECTV)
+      real albedo_co2_ice_SPECTV(L_NSPECTV)
       real pctsrf_sic(ngrid), sea_ice(ngrid)
 
@@ -85,5 +88,5 @@
 !     -----
       real a,b,E
-      integer ig,iq, icap ! wld like to remove icap
+      integer ig,iq, nw
       real fsnoi, subli, fauxo
       real twater(ngrid)
@@ -130,9 +133,5 @@
 !                      iliq is used in place of igcm_h2o_vap ONLY on the 
 !                      surface.
-
 !                      Soon to be extended to the entire water cycle...
-
-!     Ice albedo = snow albedo for now
-         albedoice=albedosnow
 
 !     Total ocean surface area
@@ -188,61 +187,67 @@
 !               albedo(ig) = MAX(MIN(albedo(ig),0.60),0.04)
 !            else
-               albedo(ig) = alb_ocean ! modif Benjamin
+               do nw=1,L_NSPECTV
+                  albedo(ig,nw) = alb_ocean ! For now, alb_ocean is defined in slab_ice_h.F90. Later we could introduce spectral dependency for alb_ocean.
+	       enddo
 !            end if
 
 
-          if(ok_slab_ocean) then
-
-! Fraction neige (hauteur critique 45kg/m2~15cm)
-            zfra = MAX(0.0,MIN(1.0,zqsurf(ig,iice)/45.0))
-! Albedo glace
-            alb_ice=alb_ice_max-(alb_ice_max-alb_ice_min) &
-            *exp(-sea_ice(ig)/h_alb_ice)
-! Albedo glace+neige
-            albedo_sic= albedosnow*zfra + alb_ice*(1.0-zfra)
-
-! Albedo final
-            albedo(ig) = pctsrf_sic(ig)*albedo_sic + (1.-pctsrf_sic(ig))*alb_ocean
-!     oceanic ice height, just for diagnostics
-            hice(ig)    = MIN(10.,sea_ice(ig)/rhowater)
-          else !ok_slab_ocean
+            if(ok_slab_ocean) then
+          
+               zfra = MAX(0.0,MIN(1.0,zqsurf(ig,iice)/45.0))     ! Snow Fraction (Critical height 45kg/m2~15cm)
+               alb_ice=alb_ice_max-(alb_ice_max-alb_ice_min) &   ! Ice Albedo
+               *exp(-sea_ice(ig)/h_alb_ice)
+               ! Albedo final calculation :
+               do nw=1,L_NSPECTV
+                  albedo(ig,nw) = pctsrf_sic(ig)*                                        &
+                                 (albedo_snow_SPECTV(nw)*zfra + alb_ice*(1.0-zfra))      &
+                               + (1.-pctsrf_sic(ig))*alb_ocean
+	       enddo
+
+               ! Oceanic ice height, just for diagnostics
+               hice(ig)    = MIN(10.,sea_ice(ig)/rhowater)
+            else !ok_slab_ocean
 
 
 !     calculate oceanic ice height including the latent heat of ice formation
 !     hice is the height of oceanic ice with a maximum of maxicethick.
-            hice(ig)    = zqsurf(ig,iice)/rhowater ! update hice to include recent snowfall
-
-!            twater(ig)  = tsurf(ig) + ptimestep*zdtsurf(ig) &
-            twater(ig)  = ztsurf(ig) - hice(ig)*RLFTT*rhowater/pcapcal(ig)
-            ! this is temperature water would have if we melted entire ocean ice layer
-            hicebis(ig) = hice(ig)
-            hice(ig)    = 0.
-
-            if(twater(ig) .lt. T_h2O_ice_liq)then
-               E=min((T_h2O_ice_liq+Tsaldiff-twater(ig))*pcapcal(ig),RLFTT*rhowater*maxicethick)
-               hice(ig)        = E/(RLFTT*rhowater)
-               hice(ig)        = max(hice(ig),0.0)
-               hice(ig)        = min(hice(ig),maxicethick)
-               pdtsurf_hyd(ig) = (hice(ig) - hicebis(ig))*RLFTT*rhowater/pcapcal(ig)/ptimestep              
-               albedo(ig)      = albedoice
-
-!               if (zqsurf(ig,iice).ge.snowlayer) then
-!                  albedo(ig) = albedoice
-!               else
-!                  albedo(ig) = albedoocean &
-!                       + (albedosnow - albedoocean)*zqsurf(ig,iice)/snowlayer
-!               endif
-
-            else
-
-               pdtsurf_hyd(ig) = -hicebis(ig)*RLFTT*rhowater/pcapcal(ig)/ptimestep                
-               albedo(ig)      = alb_ocean
-
-            endif
-
-            zqsurf(ig,iliq) = zqsurf(ig,iliq)-(hice(ig)*rhowater-zqsurf(ig,iice))
-            zqsurf(ig,iice) = hice(ig)*rhowater
-
-          endif!(ok_slab_ocean)
+               hice(ig)    = zqsurf(ig,iice)/rhowater ! update hice to include recent snowfall
+
+!              twater(ig)  = tsurf(ig) + ptimestep*zdtsurf(ig) &
+               twater(ig)  = ztsurf(ig) - hice(ig)*RLFTT*rhowater/pcapcal(ig)
+               ! this is temperature water would have if we melted entire ocean ice layer
+               hicebis(ig) = hice(ig)
+               hice(ig)    = 0.
+
+               if(twater(ig) .lt. T_h2O_ice_liq)then
+                  E=min((T_h2O_ice_liq+Tsaldiff-twater(ig))*pcapcal(ig),RLFTT*rhowater*maxicethick)
+                  hice(ig)        = E/(RLFTT*rhowater)
+                  hice(ig)        = max(hice(ig),0.0)
+                  hice(ig)        = min(hice(ig),maxicethick)
+                  pdtsurf_hyd(ig) = (hice(ig) - hicebis(ig))*RLFTT*rhowater/pcapcal(ig)/ptimestep              
+                  do nw=1,L_NSPECTV
+		     albedo(ig,nw) = albedo_snow_SPECTV(nw) ! Albedo of ice has been replaced by albedo of snow here. MT2015.
+	          enddo  
+
+!                 if (zqsurf(ig,iice).ge.snowlayer) then
+!                    albedo(ig) = albedoice
+!                 else
+!                    albedo(ig) = albedoocean &
+!                    + (albedosnow - albedoocean)*zqsurf(ig,iice)/snowlayer
+!                 endif
+
+               else
+
+                  pdtsurf_hyd(ig) = -hicebis(ig)*RLFTT*rhowater/pcapcal(ig)/ptimestep
+                  DO nw=1,L_NSPECTV
+                     albedo(ig,nw) = alb_ocean
+	          ENDDO               
+
+               endif
+
+               zqsurf(ig,iliq) = zqsurf(ig,iliq)-(hice(ig)*rhowater-zqsurf(ig,iice))
+               zqsurf(ig,iice) = hice(ig)*rhowater
+
+            endif!(ok_slab_ocean)
 
 
@@ -299,10 +304,17 @@
 
 !     re-calculate continental albedo
-            albedo(ig) = albedo0(ig) ! albedo0 = base values
+            DO nw=1,L_NSPECTV
+	       albedo(ig,nw) = albedo_bareground(ig)
+	    ENDDO
             if (zqsurf(ig,iice).ge.snowlayer) then
-               albedo(ig) = albedosnow
+               DO nw=1,L_NSPECTV
+                  albedo(ig,nw) = albedo_snow_SPECTV(nw)
+	       ENDDO
             else
-               albedo(ig) = albedo0(ig) &
-                    + (albedosnow - albedo0(ig))*zqsurf(ig,iice)/snowlayer
+               DO nw=1,L_NSPECTV
+                  albedo(ig,nw) = albedo_bareground(ig)                            & 
+                               + (albedo_snow_SPECTV(nw) - albedo_bareground(ig))  &
+                                 *zqsurf(ig,iice)/snowlayer
+	       ENDDO
             endif
 
@@ -370,12 +382,13 @@
       if(co2cond)then
 
-         icap=1
          do ig=1,ngrid
-            if (qsurf(ig,igcm_co2_ice).gt.0) then
-               albedo(ig) = albedice(icap)
-            endif
-         enddo
+            if (qsurf(ig,igcm_co2_ice).gt.1.) then ! Condition changed - Need now ~1 mm CO2 ice coverage. MT2015
+	       DO nw=1,L_NSPECTV
+                  albedo(ig,nw) = albedo_co2_ice_SPECTV(nw)
+	       ENDDO
+            endif
+         enddo ! ngrid
          
-      endif
+      endif ! co2cond
 
 
Index: /trunk/LMDZ.GENERIC/libf/phystd/inifis.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/inifis.F	(revision 1481)
+++ /trunk/LMDZ.GENERIC/libf/phystd/inifis.F	(revision 1482)
@@ -619,4 +619,9 @@
          call getin_p("albedosnow",albedosnow)
          write(*,*) " albedosnow = ",albedosnow
+         
+         write(*,*) "CO2 ice albedo ?"
+         albedoco2ice=0.5       ! default value
+         call getin_p("albedoco2ice",albedoco2ice)
+         write(*,*) " albedoco2ice = ",albedoco2ice
 
          write(*,*) "Maximum ice thickness ?"
Index: /trunk/LMDZ.GENERIC/libf/phystd/phyredem.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/phyredem.F90	(revision 1481)
+++ /trunk/LMDZ.GENERIC/libf/phystd/phyredem.F90	(revision 1482)
@@ -12,5 +12,5 @@
   use comgeomfi_h, only: area
   use surfdat_h, only: zmea, zstd, zsig, zgam, zthe, &
-                       albedice, emisice, emissiv, &
+                       emisice, emissiv,             &
                        iceradius, dtemisice, phisfi
   use iostart, only : open_restartphy, close_restartphy, & 
@@ -83,6 +83,4 @@
 
 ! Optical properties of polar caps and ground emissivity
-  tab_cntrl(22) = albedice(1)  ! Albedo of northern cap ~0.5
-  tab_cntrl(23) = albedice(2)  ! Albedo of southern cap ~0.5
   tab_cntrl(24) = emisice(1)   ! Emissivity of northern cap ~0.95
   tab_cntrl(25) = emisice(2)   ! Emissivity of southern cap ~0.95
Index: /trunk/LMDZ.GENERIC/libf/phystd/physiq.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 1481)
+++ /trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 1482)
@@ -11,8 +11,8 @@
       use watercommon_h, only : RLVTT, Psat_water,epsi
       use gases_h, only: gnom, gfrac
-      use radcommon_h, only: sigma, glat, grav
+      use radcommon_h, only: sigma, glat, grav, BWNV
       use radii_mod, only: h2o_reffrad, co2_reffrad
       use aerosol_mod, only: iaero_co2, iaero_h2o
-      use surfdat_h, only: phisfi, albedodat, zmea, zstd, zsig, zgam, zthe, &
+      use surfdat_h, only: phisfi, zmea, zstd, zsig, zgam, zthe, &
                            dryness, watercaptag
       use comdiurn_h, only: coslat, sinlat, coslon, sinlon
@@ -190,16 +190,19 @@
 ! ----------------------
 
-      integer day_ini                                  ! Initial date of the run (sol since Ls=0).
-      integer icount                                   ! Counter of calls to physiq during the run.
-      real, dimension(:),allocatable,save ::  tsurf    ! Surface temperature (K).
-      real, dimension(:,:),allocatable,save ::  tsoil  ! Sub-surface temperatures (K).
-      real, dimension(:),allocatable,save :: albedo    ! Surface albedo.
-      
-!$OMP THREADPRIVATE(tsurf,tsoil,albedo)
-
-      real,dimension(:),allocatable,save :: albedo0 ! Initial surface albedo.
-      real,dimension(:),allocatable,save :: rnat    ! Defines the type of the grid (ocean,continent,...). By BC.
-      
-!$OMP THREADPRIVATE(albedo0,rnat)
+      integer day_ini                                              ! Initial date of the run (sol since Ls=0).
+      integer icount                                               ! Counter of calls to physiq during the run.
+      real, dimension(:),allocatable,save ::  tsurf                ! Surface temperature (K).
+      real, dimension(:,:),allocatable,save ::  tsoil              ! Sub-surface temperatures (K).
+      real, dimension(:,:),allocatable,save :: albedo              ! Surface Spectral albedo. By MT2015.
+      real, dimension(:),allocatable,save :: albedo_equivalent     ! Spectral Mean albedo.
+      real, dimension(:),allocatable,save :: albedo_snow_SPECTV    ! Snow Spectral albedo.
+      real, dimension(:),allocatable,save :: albedo_co2_ice_SPECTV ! CO2 ice Spectral albedo.
+      
+!$OMP THREADPRIVATE(tsurf,tsoil,albedo,albedo_equivalent,albedo_snow_SPECTV,albedo_co2_ice_SPECTV)
+
+      real,dimension(:),allocatable,save :: albedo_bareground ! Bare Ground Albedo. By MT 2015.
+      real,dimension(:),allocatable,save :: rnat              ! Defines the type of the grid (ocean,continent,...). By BC.
+      
+!$OMP THREADPRIVATE(albedo_bareground,rnat)
 
       real,dimension(:),allocatable,save :: emis        ! Thermal IR surface emissivity.
@@ -228,21 +231,22 @@
       real pw(ngrid,nlayer)               ! Vertical velocity (m/s). (NOTE : >0 WHEN DOWNWARDS !!)
 
-      integer l,ig,ierr,iq
+      integer l,ig,ierr,iq,nw
       
       ! FOR DIAGNOSTIC :
       
-      real,dimension(:),allocatable,save :: fluxsurf_lw  ! Incident Long Wave (IR) surface flux (W.m-2).
-      real,dimension(:),allocatable,save :: fluxsurf_sw  ! incident Short Wave (stellar) surface flux (W.m-2).
-      real,dimension(:),allocatable,save :: fluxtop_lw   ! Outgoing LW (IR) flux to space (W.m-2).
-      real,dimension(:),allocatable,save :: fluxabs_sw   ! Absorbed SW (stellar) flux (W.m-2).
-      real,dimension(:),allocatable,save :: fluxtop_dn   ! Incoming SW (stellar) radiation at the top of the atmosphere (W.m-2).
-      real,dimension(:),allocatable,save :: fluxdyn      ! Horizontal heat transport by dynamics (W.m-2).
-      real,dimension(:,:),allocatable,save :: OLR_nu     ! Outgoing LW radiation in each band (Normalized to the band width (W/m2/cm-1)).
-      real,dimension(:,:),allocatable,save :: OSR_nu     ! Outgoing SW radiation in each band (Normalized to the band width (W/m2/cm-1)).
-      real,dimension(:,:),allocatable,save :: zdtlw      ! LW heating tendencies (K/s).
-      real,dimension(:,:),allocatable,save :: zdtsw      ! SW heating tendencies (K/s).
-      real,dimension(:),allocatable,save :: sensibFlux   ! Turbulent flux given by the atmosphere to the surface (W.m-2).
-      
-!$OMP THREADPRIVATE(fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn,fluxdyn,OLR_nu,OSR_nu,&
+      real,dimension(:),allocatable,save :: fluxsurf_lw     ! Incident Long Wave (IR) surface flux (W.m-2).
+      real,dimension(:),allocatable,save :: fluxsurf_sw     ! Incident Short Wave (stellar) surface flux (W.m-2).
+      real,dimension(:),allocatable,save :: fluxsurfabs_sw  ! Absorbed Short Wave (stellar) flux by the surface (W.m-2).
+      real,dimension(:),allocatable,save :: fluxtop_lw      ! Outgoing LW (IR) flux to space (W.m-2).
+      real,dimension(:),allocatable,save :: fluxabs_sw      ! Absorbed SW (stellar) flux (W.m-2).
+      real,dimension(:),allocatable,save :: fluxtop_dn      ! Incoming SW (stellar) radiation at the top of the atmosphere (W.m-2).
+      real,dimension(:),allocatable,save :: fluxdyn         ! Horizontal heat transport by dynamics (W.m-2).
+      real,dimension(:,:),allocatable,save :: OLR_nu        ! Outgoing LW radiation in each band (Normalized to the band width (W/m2/cm-1)).
+      real,dimension(:,:),allocatable,save :: OSR_nu        ! Outgoing SW radiation in each band (Normalized to the band width (W/m2/cm-1)).
+      real,dimension(:,:),allocatable,save :: zdtlw         ! LW heating tendencies (K/s).
+      real,dimension(:,:),allocatable,save :: zdtsw         ! SW heating tendencies (K/s).
+      real,dimension(:),allocatable,save :: sensibFlux      ! Turbulent flux given by the atmosphere to the surface (W.m-2).
+      
+!$OMP THREADPRIVATE(fluxsurf_lw,fluxsurf_sw,fluxsurfabs_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn,fluxdyn,OLR_nu,OSR_nu,&
 	
         !$OMP zdtlw,zdtsw,sensibFlux)
@@ -376,7 +380,10 @@
       logical clearsky ! For double radiative transfer call. By BC
       
-      real fluxsurf_lw1(ngrid), fluxsurf_sw1(ngrid), fluxtop_lw1(ngrid), fluxabs_sw1(ngrid) ! For SW/LW flux diagnostics.
-      real tau_col1(ngrid)                                                                  ! For aerosol optical depth diagnostic.
-      real OLR_nu1(ngrid,L_NSPECTI), OSR_nu1(ngrid,L_NSPECTV)                               ! For Outgoing Radiation diagnostics.
+      ! For Clear Sky Case.
+      real fluxsurf_lw1(ngrid), fluxsurf_sw1(ngrid), fluxsurfabs_sw1(ngrid)  ! For SW/LW flux.
+      real fluxtop_lw1(ngrid), fluxabs_sw1(ngrid)                            ! For SW/LW flux.
+      real albedo_equivalent1(ngrid)                                         ! For Equivalent albedo calculation.
+      real tau_col1(ngrid)                                                   ! For aerosol optical depth diagnostic.
+      real OLR_nu1(ngrid,L_NSPECTI), OSR_nu1(ngrid,L_NSPECTV)                ! For Outgoing Radiation diagnostics.
       real tf, ntf
 
@@ -435,6 +442,9 @@
         ALLOCATE(tsurf(ngrid))
         ALLOCATE(tsoil(ngrid,nsoilmx))    
-        ALLOCATE(albedo(ngrid))          
-        ALLOCATE(albedo0(ngrid))          
+        ALLOCATE(albedo(ngrid,L_NSPECTV))
+ 	ALLOCATE(albedo_equivalent(ngrid))	
+ 	ALLOCATE(albedo_snow_SPECTV(L_NSPECTV))
+ 	ALLOCATE(albedo_co2_ice_SPECTV(L_NSPECTV))
+ 	ALLOCATE(albedo_bareground(ngrid))           
         ALLOCATE(rnat(ngrid))         
         ALLOCATE(emis(ngrid))   
@@ -457,4 +467,5 @@
         ALLOCATE(fluxsurf_lw(ngrid))
         ALLOCATE(fluxsurf_sw(ngrid))
+        ALLOCATE(fluxsurfabs_sw(ngrid))
         ALLOCATE(fluxtop_lw(ngrid))
         ALLOCATE(fluxabs_sw(ngrid))
@@ -521,10 +532,17 @@
          write (*,*) 'In physiq day_ini =', day_ini
 
-!        Initialize albedo and orbital calculation.
-!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         call surfini(ngrid,nq,qsurf,albedo0)
+!        Initialize albedo calculation. 
+!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+         albedo(:,:)=0.0
+ 	 albedo_bareground(:)=0.0
+ 	 albedo_snow_SPECTV(:)=0.0
+ 	 albedo_co2_ice_SPECTV(:)=0.0
+         call surfini(ngrid,nq,qsurf,albedo,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV) 
+         
+         
+!        Initialize orbital calculation. 
+!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
 
-         albedo(:)=albedo0(:)
 
          if(tlocked)then
@@ -677,8 +695,9 @@
          if (ngrid.ne.1) then ! Note : no need to create a restart file in 1d.
             call physdem0("restartfi.nc",long,lati,nsoilmx,ngrid,nlayer,nq, &
-                          ptimestep,pday+nday,time_phys,area, &
-                          albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
+                          ptimestep,pday+nday,time_phys,area,               &
+                          albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
          endif
          
+               
       endif ! end of 'firstcall'
 
@@ -848,37 +867,39 @@
                ! standard callcorrk
                clearsky=.false.
-               call callcorrk(ngrid,nlayer,pq,nq,qsurf,                          &
-                              albedo,emis,mu0,pplev,pplay,pt,                    &
-                              tsurf,fract,dist_star,aerosol,muvar,               &
-                              zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,    &
-                              fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,               &
-                              tau_col,cloudfrac,totcloudfrac,                    &
-                              clearsky,firstcall,lastcall)      
-
-                ! Option to call scheme once more for clear regions
+               call callcorrk(ngrid,nlayer,pq,nq,qsurf,                           &
+                              albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt,   &
+                              tsurf,fract,dist_star,aerosol,muvar,                &
+                              zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,                &
+                              fluxsurfabs_sw,fluxtop_lw,                          &
+                              fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,                &
+                              tau_col,cloudfrac,totcloudfrac,                     &
+                              clearsky,firstcall,lastcall)     
+
+                ! Option to call scheme once more for clear regions  
                if(CLFvarying)then
 
                   ! ---> PROBLEMS WITH ALLOCATED ARRAYS : temporary solution in callcorrk: do not deallocate if CLFvarying ...
                   clearsky=.true.
-                  call callcorrk(ngrid,nlayer,pq,nq,qsurf,                            &
-                                 albedo,emis,mu0,pplev,pplay,pt,                      &
-                                 tsurf,fract,dist_star,aerosol,muvar,                 &
-                                 zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1,fluxtop_lw1, &
-                                 fluxabs_sw1,fluxtop_dn,OLR_nu1,OSR_nu1,              &
-                                 tau_col1,cloudfrac,totcloudfrac,                     &
+                  call callcorrk(ngrid,nlayer,pq,nq,qsurf,                           &
+                                 albedo,albedo_equivalent1,emis,mu0,pplev,pplay,pt,   &
+                                 tsurf,fract,dist_star,aerosol,muvar,                &
+                                 zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1,            &
+                                 fluxsurfabs_sw1,fluxtop_lw1,                        &
+                                 fluxabs_sw1,fluxtop_dn,OLR_nu1,OSR_nu1,             &
+                                 tau_col1,cloudfrac,totcloudfrac,                    &
                                  clearsky,firstcall,lastcall)
                   clearsky = .false.  ! just in case.     
-
 
                   ! Sum the fluxes and heating rates from cloudy/clear cases
                   do ig=1,ngrid
                      tf=totcloudfrac(ig)
-                     ntf=1.-tf
-                    
-                     fluxsurf_lw(ig) = ntf*fluxsurf_lw1(ig) + tf*fluxsurf_lw(ig)
-                     fluxsurf_sw(ig) = ntf*fluxsurf_sw1(ig) + tf*fluxsurf_sw(ig)
-                     fluxtop_lw(ig)  = ntf*fluxtop_lw1(ig)  + tf*fluxtop_lw(ig)
-                     fluxabs_sw(ig)  = ntf*fluxabs_sw1(ig)  + tf*fluxabs_sw(ig)
-                     tau_col(ig)     = ntf*tau_col1(ig)     + tf*tau_col(ig)
+                     ntf=1.-tf                    
+                     fluxsurf_lw(ig)       = ntf*fluxsurf_lw1(ig)       + tf*fluxsurf_lw(ig)
+                     fluxsurf_sw(ig)       = ntf*fluxsurf_sw1(ig)       + tf*fluxsurf_sw(ig)
+                     albedo_equivalent(ig) = ntf*albedo_equivalent1(ig) + tf*albedo_equivalent(ig)
+                     fluxsurfabs_sw(ig)    = ntf*fluxsurfabs_sw1(ig)    + tf*fluxsurfabs_sw(ig)
+                     fluxtop_lw(ig)        = ntf*fluxtop_lw1(ig)        + tf*fluxtop_lw(ig)
+                     fluxabs_sw(ig)        = ntf*fluxabs_sw1(ig)        + tf*fluxabs_sw(ig)
+                     tau_col(ig)           = ntf*tau_col1(ig)           + tf*tau_col(ig)
                     
                      zdtlw(ig,1:nlayer) = ntf*zdtlw1(ig,1:nlayer) + tf*zdtlw(ig,1:nlayer)
@@ -887,5 +908,5 @@
                      OSR_nu(ig,1:L_NSPECTV) = ntf*OSR_nu1(ig,1:L_NSPECTV) + tf*OSR_nu(ig,1:L_NSPECTV)		       
                      OLR_nu(ig,1:L_NSPECTI) = ntf*OLR_nu1(ig,1:L_NSPECTI) + tf*OLR_nu(ig,1:L_NSPECTI)		       
-                 enddo
+                  enddo                               
 
                endif ! end of CLFvarying.
@@ -896,7 +917,7 @@
 
 
-               ! Radiative flux from the sky absorbed by the surface (W.m-2)
+               ! Radiative flux from the sky absorbed by the surface (W.m-2).
                GSR=0.0
-               fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)+fluxsurf_sw(1:ngrid)*(1.-albedo(1:ngrid))
+               fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)+fluxsurfabs_sw(1:ngrid)
 
                             !if(noradsurf)then ! no lower surface; SW flux just disappears
@@ -910,7 +931,8 @@
 
             elseif(newtonian)then
-
-               ! II.b Call Newtonian cooling scheme
-               ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+            
+! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+! II.b Call Newtonian cooling scheme
+! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
                call newtrelax(ngrid,nlayer,mu0,sinlat,zpopsk,pt,pplay,pplev,dtrad,firstcall)
 
@@ -928,5 +950,17 @@
                endif
                fluxsurf_sw(1:ngrid) = fluxtop_dn(1:ngrid)
-               fluxrad_sky(1:ngrid) = fluxtop_dn(1:ngrid)*(1.-albedo(1:ngrid))
+               print*,'------------WARNING---WARNING------------' ! by MT2015.
+               print*,'You are in corrk=false mode, '
+               print*,'and the surface albedo is taken equal to its value at 0.5 micrometers ...'
+               if ( (10000./BWNV(1) .le. 0.5) .or. (10000./BWNV(L_NSPECTV) .ge. 0.5) ) then
+                  print*,'0.5 microns band doesnt match your visible bands ! Abort !'
+                  call abort
+               endif
+               DO nw=1, L_NSPECTV-1
+                  if ( (10000./BWNV(nw) .ge. 0.5) .and. (10000./BWNV(nw+1) .lt. 0.5) ) then
+                     fluxsurfabs_sw(1:ngrid) = fluxtop_dn(1:ngrid)*(1.-albedo(1:ngrid,nw))
+                     fluxrad_sky(1:ngrid)    = fluxsurfabs_sw(1:ngrid)
+                  endif
+               ENDDO
                fluxtop_lw(1:ngrid)  = emis(1:ngrid)*sigma*tsurf(1:ngrid)**4
 
@@ -950,5 +984,6 @@
 	    call planetwide_sumval(cpp*massarea(:,:)*zdtsw(:,:)/totarea_planet,dEtotSW)
 	    call planetwide_sumval(cpp*massarea(:,:)*zdtlw(:,:)/totarea_planet,dEtotLW)
-	    call planetwide_sumval(fluxsurf_sw(:)*(1.-albedo(:))*area(:)/totarea_planet,dEtotsSW) !JL13 carefull, albedo can have changed since the last time we called corrk
+	    !call planetwide_sumval(fluxsurf_sw(:)*(1.-albedo_equivalent(:))*area(:)/totarea_planet,dEtotsSW) !JL13 carefull, albedo can have changed since the last time we called corrk
+            call planetwide_sumval(fluxsurfabs_sw(:)*area(:)/totarea_planet,dEtotsSW) !JL13 carefull, albedo can have changed since the last time we called corrk
 	    call planetwide_sumval((fluxsurf_lw(:)*emis(:)-zplanck(:))*area(:)/totarea_planet,dEtotsLW)
 	    dEzRadsw(:,:)=cpp*mass(:,:)*zdtsw(:,:)
@@ -1166,4 +1201,5 @@
                            pphi,pdt,pdu,pdv,zdtsurf,pu,pv,pq,pdq,        &
                            qsurf(1,igcm_co2_ice),albedo,emis,            &
+                           albedo_bareground,albedo_co2_ice_SPECTV,      &
                            zdtc,zdtsurfc,pdpsrf,zduc,zdvc,               &
                            zdqc)
@@ -1461,6 +1497,8 @@
          if(hydrology)then
             
-            call hydrol(ngrid,nq,ptimestep,rnat,tsurf,qsurf,dqsurf,dqs_hyd,             &
-                        capcal,albedo0,albedo,mu0,zdtsurf,zdtsurf_hyd,hice,pctsrf_sic,  &
+            call hydrol(ngrid,nq,ptimestep,rnat,tsurf,qsurf,dqsurf,dqs_hyd, &
+                        capcal,albedo,albedo_bareground,                    &
+	                albedo_snow_SPECTV,albedo_co2_ice_SPECTV,           &
+                        mu0,zdtsurf,zdtsurf_hyd,hice,pctsrf_sic,            &
                         sea_ice)
 
@@ -1848,5 +1886,6 @@
          call wstats(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
          call wstats(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
-         call wstats(ngrid,"ALB","Surface albedo"," ",2,albedo)
+         !call wstats(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent)
+         !call wstats(ngrid,"ALB_1st","First Band Surface albedo"," ",2,albedo(:,1))
          call wstats(ngrid,"p","Pressure","Pa",3,pplay)
          call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
@@ -1955,5 +1994,6 @@
       if(callrad.and.(.not.newtonian))then
       
-         call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo)
+         !call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent)
+         !call writediagfi(ngrid,"ALB_1st","First Band Surface albedo"," ",2,albedo(:,1))
          call writediagfi(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
          call writediagfi(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
Index: /trunk/LMDZ.GENERIC/libf/phystd/sfluxv.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/sfluxv.F	(revision 1481)
+++ /trunk/LMDZ.GENERIC/libf/phystd/sfluxv.F	(revision 1482)
@@ -24,5 +24,6 @@
 
       integer L, NG, NW, NG1,k
-      real*8  rsfv, ubar0, f0pi, btop, bsurf, taumax, eterm
+      real*8 ubar0, f0pi, btop, bsurf, taumax, eterm
+      real*8 rsfv(L_NSPECTV) ! Spectral dependency added by MT2015.
       real*8 FZEROV(L_NSPECTV)
 
@@ -83,5 +84,5 @@
 
           ETERM = MIN(TAUV(L_NLEVRAD,NW,NG),TAUMAX)
-          BSURF = RSFV*UBAR0*STEL(NW)*EXP(-ETERM/UBAR0)
+          BSURF = RSFV(NW)*UBAR0*STEL(NW)*EXP(-ETERM/UBAR0)
 
 C         WE CAN NOW SOLVE FOR THE COEFFICIENTS OF THE TWO STREAM
@@ -94,5 +95,5 @@
 
           CALL GFLUXV(DTAUV(1,NW,NG),TAUV(1,NW,NG),TAUCUMV(1,NW,NG),
-     *                WBARV(1,NW,NG),COSBV(1,NW,NG),UBAR0,F0PI,RSFV,
+     *                WBARV(1,NW,NG),COSBV(1,NW,NG),UBAR0,F0PI,RSFV(NW),    
      *                BTOP,BSURF,FMUPV,FMDV,DIFFV,FLUXUP,FLUXDN)
 
@@ -141,5 +142,5 @@
  
         ETERM = MIN(TAUV(L_NLEVRAD,NW,NG),TAUMAX)
-        BSURF = RSFV*UBAR0*STEL(NW)*EXP(-ETERM/UBAR0)
+        BSURF = RSFV(NW)*UBAR0*STEL(NW)*EXP(-ETERM/UBAR0)
 
 
@@ -152,5 +153,5 @@
 
         CALL GFLUXV(DTAUV(1,NW,NG),TAUV(1,NW,NG),TAUCUMV(1,NW,NG),
-     *              WBARV(1,NW,NG),COSBV(1,NW,NG),UBAR0,F0PI,RSFV,
+     *              WBARV(1,NW,NG),COSBV(1,NW,NG),UBAR0,F0PI,RSFV(NW),
      *              BTOP,BSURF,FMUPV,FMDV,DIFFV,FLUXUP,FLUXDN)
 
Index: /trunk/LMDZ.GENERIC/libf/phystd/surfdat_h.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/surfdat_h.F90	(revision 1481)
+++ /trunk/LMDZ.GENERIC/libf/phystd/surfdat_h.F90	(revision 1482)
@@ -4,5 +4,5 @@
        implicit none
 
-       real,allocatable,dimension(:) :: albedodat ! albedo of bare ground
+       real,allocatable,dimension(:) :: albedodat ! albedo of bare ground stocked in startfi.nc file.
 !$OMP THREADPRIVATE(albedodat)
        ! Ehouarn: moved inertiedat to comsoil.h
@@ -10,9 +10,8 @@
        real,allocatable,dimension(:) :: phisfi ! geopotential at ground level
 !$OMP THREADPRIVATE(phisfi)
-       real,dimension(2) :: albedice
        real,dimension(2) :: emisice ! ice emissivity; 1:Northern hemisphere 2:Southern hemisphere
        real emissiv
        real,dimension(2) :: iceradius, dtemisice
-!$OMP THREADPRIVATE(albedice,emisice,emissiv,iceradius,dtemisice)
+!$OMP THREADPRIVATE(emisice,emissiv,iceradius,dtemisice)
        real,allocatable,dimension(:) :: zmea,zstd,zsig,zgam,zthe
 !$OMP THREADPRIVATE(zmea,zstd,zsig,zgam,zthe)
Index: /trunk/LMDZ.GENERIC/libf/phystd/surfini.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/surfini.F	(revision 1481)
+++ /trunk/LMDZ.GENERIC/libf/phystd/surfini.F	(revision 1482)
@@ -1,62 +1,76 @@
-      SUBROUTINE surfini(ngrid,nq,qsurf,psolaralb)
+      SUBROUTINE surfini(ngrid,nq,qsurf,albedo,albedo_bareground,
+     &                   albedo_snow_SPECTV,albedo_co2_ice_SPECTV)
 
-      USE surfdat_h, only: albedodat, albedice
+      USE surfdat_h, only: albedodat
       USE tracer_h, only: igcm_co2_ice
-      use comgeomfi_h, only: lati
       use planetwide_mod, only: planetwide_maxval, planetwide_minval
+      use radinc_h, only : L_NSPECTV
+      use callkeys_mod, only : albedosnow, albedoco2ice
 
       IMPLICIT NONE
-c=======================================================================
-c
-c   creation des calottes pour l'etat initial
-c
-c=======================================================================
-c-----------------------------------------------------------------------
+      
+      
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+cccccccccccccc                                                                 cccccccccccccc
+cccccccccccccc   Spectral Albedo Initialisation - Routine modified by MT2015.  cccccccccccccc
+cccccccccccccc                                                                 cccccccccccccc 
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+
+c--------------------
 c   Declarations:
-c   -------------
-!#include "dimensions.h"
-!#include "dimphys.h"
-c
+c--------------------
+
       INTEGER,INTENT(IN) :: ngrid
       INTEGER,INTENT(IN) :: nq
-      REAL,INTENT(OUT) :: psolaralb(ngrid)
-      REAL,INTENT(IN) :: qsurf(ngrid,nq) !tracer on surface (kg/m2)
+      REAL,INTENT(OUT) :: albedo(ngrid,L_NSPECTV)
+      REAL,INTENT(OUT) :: albedo_bareground(ngrid)
+      REAL,INTENT(OUT) :: albedo_snow_SPECTV(L_NSPECTV)
+      REAL,INTENT(OUT) :: albedo_co2_ice_SPECTV(L_NSPECTV)
+      REAL,INTENT(IN) :: qsurf(ngrid,nq) ! tracer on surface (kg/m2)
 
-      INTEGER :: ig,icap
+      INTEGER :: ig,nw
       REAL :: min_albedo,max_albedo
-c
+
 c=======================================================================
 
-
-      DO ig=1,ngrid
-         psolaralb(ig)=albedodat(ig)
+      ! Step 1 : Initialisation of the Spectral Albedos.
+      DO nw=1,L_NSPECTV
+         albedo_snow_SPECTV(nw)=albedosnow
+         albedo_co2_ice_SPECTV(nw)=albedoco2ice
       ENDDO
 
-      call planetwide_minval(albedodat,min_albedo)
-      call planetwide_maxval(albedodat,max_albedo)
-      write(*,*) 'surfini: minimum corrected albedo',min_albedo
-      write(*,*) 'surfini: maximum corrected albedo',max_albedo
 
+      ! Step 2 : We get the bare ground albedo from the start files.
+      DO ig=1,ngrid
+         albedo_bareground(ig)=albedodat(ig)
+	 DO nw=1,L_NSPECTV
+	    albedo(ig,nw)=albedo_bareground(ig)
+	 ENDDO
+      ENDDO
+      call planetwide_minval(albedo_bareground,min_albedo)
+      call planetwide_maxval(albedo_bareground,max_albedo)
+      write(*,*) 'surfini: minimum bare ground albedo',min_albedo
+      write(*,*) 'surfini: maximum bare ground albedo',max_albedo
+
+
+      ! Step 3 : We modify the albedo considering some CO2 at the surface. We dont take into account water ice (this is made in hydrol after the first timestep) ...
       if (igcm_co2_ice.ne.0) then
-! Change Albedo if there is CO2 ice on the surface
-        DO ig=1,ngrid
-          IF (qsurf(ig,igcm_co2_ice) .GT. 0.) THEN
-             IF(lati(ig).LT.0.) THEN
-                icap=2 ! Southern hemisphere
-             ELSE
-                icap=1 ! Northern hemisphere
-             ENDIF
-             psolaralb(ig) = albedice(icap)
-          END IF
-        ENDDO ! of DO ig=1,ngrid     
+         DO ig=1,ngrid
+            IF (qsurf(ig,igcm_co2_ice) .GT. 1.) THEN ! This was changed by MT2015. Condition for ~1mm of CO2 ice deposit.
+	       DO nw=1,L_NSPECTV
+	          albedo(ig,nw)=albedo_co2_ice_SPECTV(nw)
+	       ENDDO
+            END IF   
+         ENDDO   
       else
-        write(*,*) "surfini: No CO2 ice tracer on surface  ..."
-        write(*,*) "         and therefore no albedo change."
-      endif
+         write(*,*) "surfini: No CO2 ice tracer on surface  ..."
+         write(*,*) "         and therefore no albedo change."
+      endif     
+      call planetwide_minval(albedo,min_albedo)
+      call planetwide_maxval(albedo,max_albedo)
+      write(*,*) 'surfini: minimum corrected initial albedo',min_albedo
+      write(*,*) 'surfini: maximum corrected initial albedo',max_albedo
 
-      call planetwide_minval(psolaralb,min_albedo)
-      call planetwide_maxval(psolaralb,max_albedo)
-      write(*,*) 'surfini: minimum corrected albedo',min_albedo
-      write(*,*) 'surfini: maximum corrected albedo',max_albedo
 
       END
Index: /trunk/LMDZ.GENERIC/libf/phystd/tabfi.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/tabfi.F	(revision 1481)
+++ /trunk/LMDZ.GENERIC/libf/phystd/tabfi.F	(revision 1482)
@@ -46,5 +46,5 @@
       use ioipsl_getincom , only: getin
 
-      use surfdat_h, only: albedice, emisice, iceradius, dtemisice,
+      use surfdat_h, only: emisice, iceradius, dtemisice,
      &                     emissiv
       use comsoil_h, only: volcapa
@@ -133,6 +133,4 @@
       emin_turb = tab_cntrl(tab0+21)
 c optical properties of polar caps and ground emissivity
-      albedice(1)= tab_cntrl(tab0+22)
-      albedice(2)= tab_cntrl(tab0+23)
       emisice(1) = tab_cntrl(tab0+24)
       emisice(2) = tab_cntrl(tab0+25)
@@ -190,6 +188,4 @@
       write(*,5) '(24)     emisice(1)',tab_cntrl(tab0+24),emisice(1)
       write(*,5) '(25)     emisice(2)',tab_cntrl(tab0+25),emisice(2)
-      write(*,5) '(22)    albedice(1)',tab_cntrl(tab0+22),albedice(1)
-      write(*,5) '(23)    albedice(2)',tab_cntrl(tab0+23),albedice(2)
       write(*,6) '(31)   iceradius(1)',tab_cntrl(tab0+31),iceradius(1)
       write(*,6) '(32)   iceradius(2)',tab_cntrl(tab0+32),iceradius(2)
@@ -225,5 +221,4 @@
       write(*,*) '(26)         emissiv : ground emissivity'
       write(*,*) '(24 et 25)   emisice : CO2 ice max emissivity '
-      write(*,*) '(22 et 23)  albedice : CO2 ice cap albedos'
       write(*,*) '(31 et 32) iceradius : mean scat radius of CO2 snow'
       write(*,*) '(33 et 34) dtemisice : time scale for snow',
@@ -315,20 +310,4 @@
           write(*,*) ' emisice(2) (new value):',emisice(2)
 
-        else if (modif(1:len_trim(modif)) .eq. 'albedice') then
-          write(*,*) 'current value albedice(1) North:',albedice(1)
-          write(*,*) 'enter new value:'
- 108      read(*,*,iostat=ierr) albedice(1)
-          if(ierr.ne.0) goto 108
-          write(*,*) 
-          write(*,*) ' albedice(1) (new value):',albedice(1)
-          write(*,*)
-
-          write(*,*) 'current value albedice(2) South:',albedice(2)
-          write(*,*) 'enter new value:'
- 109      read(*,*,iostat=ierr) albedice(2)
-          if(ierr.ne.0) goto 109
-          write(*,*) 
-          write(*,*) ' albedice(2) (new value):',albedice(2)
-
         else if (modif(1:len_trim(modif)) .eq. 'iceradius') then
           write(*,*) 'current value iceradius(1) North:',iceradius(1)
@@ -518,6 +497,4 @@
       write(*,5) '(24)     emisice(1)',tab_cntrl(tab0+24),emisice(1)
       write(*,5) '(25)     emisice(2)',tab_cntrl(tab0+25),emisice(2)
-      write(*,5) '(22)    albedice(1)',tab_cntrl(tab0+22),albedice(1)
-      write(*,5) '(23)    albedice(2)',tab_cntrl(tab0+23),albedice(2)
       write(*,6) '(31)   iceradius(1)',tab_cntrl(tab0+31),iceradius(1)
       write(*,6) '(32)   iceradius(2)',tab_cntrl(tab0+32),iceradius(2)
