Index: LMDZ6/branches/Ocean_skin/libf/phylmd/surf_landice_mod.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmd/surf_landice_mod.F90	(revision 3798)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/surf_landice_mod.F90	(revision 4013)
@@ -19,5 +19,5 @@
        tsoil, z0m, z0h, SFRWL, alb_dir, alb_dif, evap, fluxsens, fluxlat, &
        tsurf_new, dflux_s, dflux_l, &
-       slope, cloudf, &
+       alt, slope, cloudf, &
        snowhgt, qsnow, to_ice, sissnow, &
        alb3, runoff, &
@@ -25,6 +25,7 @@
 
     USE dimphy
-    USE surface_data,     ONLY : type_ocean, calice, calsno, landice_opt, n_dtis
-    USE fonte_neige_mod,  ONLY : fonte_neige, run_off_lic
+    USE geometry_mod,     ONLY : longitude,latitude
+    USE surface_data,     ONLY : type_ocean, calice, calsno, landice_opt, iflag_albcalc
+    USE fonte_neige_mod,  ONLY : fonte_neige,run_off_lic,fqcalving_global,ffonte_global,fqfonte_global,runofflic_global
     USE cpl_mod,          ONLY : cpl_send_landice_fields
     USE calcul_fluxs_mod
@@ -33,7 +34,4 @@
     USE ioipsl_getin_p_mod, ONLY : getin_p
 
-#ifdef CPP_SISVAT
-    USE surf_sisvat_mod,  ONLY : surf_sisvat
-#endif
 
 #ifdef CPP_INLANDSIS
@@ -75,4 +73,5 @@
     REAL, DIMENSION(klon), INTENT(IN)             :: albedo  !mean albedo
     REAL, DIMENSION(klon), INTENT(IN)             :: pphi1    
+    REAL, DIMENSION(klon), INTENT(IN)             :: alt   !mean altitude of the grid box  
     REAL, DIMENSION(klon), INTENT(IN)             :: slope   !mean slope in grid box  
     REAL, DIMENSION(klon), INTENT(IN)             :: cloudf  !total cloud fraction
@@ -115,9 +114,10 @@
     REAL, DIMENSION(klon)    :: u0, v0, u1_lay, v1_lay, ustar
     INTEGER                  :: i,j,nt
-
+    REAL, DIMENSION(klon)    :: fqfonte,ffonte
     REAL, DIMENSION(klon)    :: emis_new                  !Emissivity
     REAL, DIMENSION(klon)    :: swdown,lwdown
-    REAL, DIMENSION(klon)    :: precip_snow_adv, snow_adv !Snow Drift precip./advection
-    REAL, DIMENSION(klon)    :: zsl_height, wind_velo      !surface layer height, wind spd
+    REAL, DIMENSION(klon)    :: precip_snow_adv, snow_adv !Snow Drift precip./advection (not used in inlandsis)
+    REAL, DIMENSION(klon)    :: erod                      !erosion of surface snow (flux, kg/m2/s like evap)
+    REAL, DIMENSION(klon)    :: zsl_height, wind_velo     !surface layer height, wind spd
     REAL, DIMENSION(klon)    :: dens_air,  snow_cont_air  !air density; snow content air
     REAL, DIMENSION(klon)    :: alb_soil                  !albedo of underlying ice
@@ -132,8 +132,6 @@
 
 
-!albedo SB >>>
-    real,dimension(klon) :: alb1,alb2
-!albedo SB <<<
-
+    REAL,DIMENSION(klon) :: alb1,alb2
+    REAL, DIMENSION (klon,6) :: alb6
 ! End definition
 !****************************************************************************************
@@ -179,82 +177,13 @@
 !****************************************************************************************
 !  landice_opt = 0 : soil_model, calcul_flux, fonte_neige, ...  
-!  landice_opt = 1  : prepare and call SISVAT snow model
-!  landice_opt = 2  : prepare and call INLANDSIS snow model
+!  landice_opt = 1  : prepare and call INterace Lmdz SISvat (INLANDSIS)
 !****************************************************************************************
 
 
     IF (landice_opt .EQ. 1) THEN
-    
-!****************************************************************************************
-! CALL to SISVAT interface
-!****************************************************************************************
-
-#ifdef CPP_SISVAT
-       ! Prepare for calling SISVAT
-       
-       ! Calculate incoming flux for SW and LW interval: swdown, lwdown
-       swdown(:)        = 0.0
-       lwdown(:)        = 0.0
-       DO i = 1, knon
-          swdown(i)        = swnet(i)/(1-albedo(i))
-          lwdown(i)        = lwdownm(i)
-       END DO
-       
-       ! Set constants and compute some input for SISVAT
-       snow_adv(:)      = 0.                          ! no snow blown in for now
-       snow_cont_air(:) = 0.       
-       alb_soil(:)      = albedo(:)
-       pref             = 100000.                     ! = 1000 hPa
-       DO i = 1, knon
-          wind_velo(i)     = u1(i)**2 + v1(i)**2
-          wind_velo(i)     = wind_velo(i)**0.5
-          pexner(i)        = (p1lay(i)/pref)**(RD/RCPD)
-          dens_air(i)      = p1lay(i)/RD/temp_air(i)  ! dry air density
-          zsl_height(i)     = pphi1(i)/RG              
-       END DO
-
-
-       ! config: compute everything with SV but temperatures afterwards with soil/calculfluxs
-       DO i = 1, knon
-          tsoil0(i,:)=tsoil(i,:)
-       END DO
-           ! Martin
-           PRINT*, 'on appelle surf_sisvat'
-           ! Martin
-       CALL surf_sisvat(knon, rlon, rlat, knindex, itime, dtime, debut, lafin, &
-            rmu0, swdown, lwdown, pexner, ps, p1lay, &
-            precip_rain, precip_snow, precip_snow_adv, snow_adv, &
-            zsl_height, wind_velo, temp_air, dens_air, spechum, tsurf, &
-            rugoro, snow_cont_air, alb_soil, slope, cloudf, &
-            radsol, qsol, tsoil0, snow, snowhgt, qsnow, to_ice,sissnow, agesno, &
-            AcoefH, AcoefQ, BcoefH, BcoefQ, cdragh, &
-            run_off_lic, evap, fluxsens, fluxlat, dflux_s, dflux_l, &        
-            tsurf_new, alb1, alb2, alb3, &
-            emis_new, z0m, qsurf)
-       z0h(1:knon)=z0m(1:knon) ! en attendant mieux
-       
-       ! Suppose zero surface speed
-       u0(:)            = 0.0
-       v0(:)            = 0.0
-       ! The calculation of heat/water fluxes, otherwise done by "CALL calcul_fluxs" is 
-       ! integrated in SISVAT, using the same method. It can be found in "sisvat.f", in the
-       ! subroutine "SISVAT_TS2". 
-       ! u0, v0=0., dif_grnd=0. and beta=1 are assumed there! 
-       
-       CALL calcul_flux_wind(knon, dtime, &
-            u0, v0, u1, v1, gustiness, cdragm, &
-            AcoefU, AcoefV, BcoefU, BcoefV, &
-            p1lay, temp_air, &
-            flux_u1, flux_v1)
-#else
-       abort_message='Pb de coherence: landice_opt = 1 mais CPP_SISVAT = .false.'
-       CALL abort_physic(modname,abort_message,1)
-#endif
-
-!****************************************************************************************
+
+!****************************************************************************************    
 ! CALL to INLANDSIS interface
 !****************************************************************************************
-
-    ELSE IF (landice_opt .EQ. 2) THEN
 #ifdef CPP_INLANDSIS
 
@@ -278,6 +207,5 @@
        swdown(:)        = 0.0
        lwdown(:)        = 0.0
-       snow_adv(:)      = 0.  ! no snow blown in for now
-       snow_cont_air(:) = 0.       
+       snow_cont_air(:) = 0.  ! the snow content in air is not a prognostic variable of the model      
        alb_soil(:)      = 0.4 ! before albedo(:) but here it is the ice albedo that we have to set
        ustar(:)         = 0.
@@ -296,33 +224,41 @@
        
 
-       ! Subtimestepping
-       
-        dtis=dtime/n_dtis
-
-        DO nt=1,n_dtis
-
-          IF (lafin .and. nt.eq.n_dtis) THEN
+
+        dtis=dtime
+
+          IF (lafin) THEN
             lafin_is=.true.
           END IF
 
-        !PRINT*,'RENTRE DANS INLANDSIS','itime',itime,'dtime',dtime,'dtis',dtis
-        CALL surf_inlandsis(knon, rlon, rlat, knindex, itime, dtis, debut_is, lafin_is, &
-            rmu0, swdown, lwdown, albedo, pexner, ps, p1lay, &
-            precip_rain, precip_snow, precip_snow_adv, snow_adv, &
-            zsl_height, wind_velo, ustar, temp_air, dens_air, spechum, tsurf, &
-            rugoro, snow_cont_air, alb_soil, slope, cloudf, &
-            radsol, qsol, tsoil0, snow, zfra, snowhgt, qsnow, to_ice,sissnow, agesno, &
+          CALL surf_inlandsis(knon, rlon, rlat, knindex, itime, dtis, debut_is, lafin_is,&
+            rmu0, swdown, lwdown, albedo, pexner, ps, p1lay, precip_rain, precip_snow,   &
+            zsl_height, wind_velo, ustar, temp_air, dens_air, spechum, tsurf,&
+            rugoro, snow_cont_air, alb_soil, alt, slope, cloudf, &
+            radsol, qsol, tsoil0, snow, zfra, snowhgt, qsnow, to_ice, sissnow,agesno,   &
             AcoefH, AcoefQ, BcoefH, BcoefQ, cdragm, cdragh, &
-            run_off_lic, evap, fluxsens, fluxlat, dflux_s, dflux_l, &        
-            tsurf_new, alb1, alb2, alb3, &
-            emis_new, z0m, z0h, qsurf)      
-
-         debut_is=.false.
-
-        END DO
+            run_off_lic, fqfonte, ffonte, evap, erod, fluxsens, fluxlat,dflux_s, dflux_l, &
+            tsurf_new, alb1, alb2, alb3, alb6, &
+            emis_new, z0m, z0h, qsurf)
+
+          debut_is=.false.
+
+
+        ! Treatment of snow melting and calving
+
+        ! for consistency with standard LMDZ, add calving to run_off_lic
+        run_off_lic(:)=run_off_lic(:) + to_ice(:)
+
+        DO i = 1, knon
+           ffonte_global(knindex(i),is_lic)    = ffonte(i)
+           fqfonte_global(knindex(i),is_lic)   = fqfonte(i)! net melting= melting - refreezing
+           fqcalving_global(knindex(i),is_lic) = to_ice(i) ! flux
+           runofflic_global(knindex(i)) = run_off_lic(i)
+        ENDDO
+        ! Here, we assume that the calving term is equal to the to_ice term
+        ! (no ice accumulation)
 
 
 #else
-       abort_message='Pb de coherence: landice_opt = 2 mais CPP_INLANDSIS = .false.'
+       abort_message='Pb de coherence: landice_opt = 1  mais CPP_INLANDSIS = .false.'
        CALL abort_physic(modname,abort_message,1)
 #endif
@@ -343,5 +279,6 @@
     ! use soil model and recalculate properly cal
     IF (soil_model) THEN 
-       CALL soil(dtime, is_lic, knon, snow, tsurf, tsoil, soilcap, soilflux)
+       CALL soil(dtime, is_lic, knon, snow, tsurf, qsol, &
+        & longitude(knindex(1:knon)), latitude(knindex(1:knon)), tsoil, soilcap, soilflux)
        cal(1:knon) = RCPD / soilcap(1:knon)
        radsol(1:knon)  = radsol(1:knon) + soilflux(1:knon)
@@ -420,7 +357,4 @@
 
 
-
-    
-
     END IF ! landice_opt
 
@@ -428,5 +362,5 @@
 !****************************************************************************************
 ! Send run-off on land-ice to coupler if coupled ocean.
-! run_off_lic has been calculated in fonte_neige or surf_sisvat
+! run_off_lic has been calculated in fonte_neige or surf_inlandsis
 !
 !****************************************************************************************
@@ -476,10 +410,20 @@
        alb_dir(1:knon,5)=alb2(1:knon)
        alb_dir(1:knon,6)=alb2(1:knon)
+
+       IF ((landice_opt .EQ. 1) .AND. (iflag_albcalc .EQ. 2)) THEN
+       alb_dir(1:knon,1)=alb6(1:knon,1)
+       alb_dir(1:knon,2)=alb6(1:knon,2)
+       alb_dir(1:knon,3)=alb6(1:knon,3)
+       alb_dir(1:knon,4)=alb6(1:knon,4)
+       alb_dir(1:knon,5)=alb6(1:knon,5)
+       alb_dir(1:knon,6)=alb6(1:knon,6)
+       ENDIF
+
      end select
 alb_dif=alb_dir
 !albedo SB <<<
 
-
-
+  
+  
 
   END SUBROUTINE surf_landice
