Index: LMDZ6/branches/Ocean_skin/libf/phylmd/surf_landice_mod.F90
===================================================================
--- LMDZ6/branches/Ocean_skin/libf/phylmd/surf_landice_mod.F90	(revision 3428)
+++ LMDZ6/branches/Ocean_skin/libf/phylmd/surf_landice_mod.F90	(revision 3798)
@@ -25,5 +25,5 @@
 
     USE dimphy
-    USE surface_data,     ONLY : type_ocean, calice, calsno, ok_snow
+    USE surface_data,     ONLY : type_ocean, calice, calsno, landice_opt, n_dtis
     USE fonte_neige_mod,  ONLY : fonte_neige, run_off_lic
     USE cpl_mod,          ONLY : cpl_send_landice_fields
@@ -36,4 +36,9 @@
     USE surf_sisvat_mod,  ONLY : surf_sisvat
 #endif
+
+#ifdef CPP_INLANDSIS
+    USE surf_inlandsis_mod,  ONLY : surf_inlandsis
+#endif
+
     USE indice_sol_mod
 
@@ -86,6 +91,6 @@
 !    REAL, DIMENSION(klon), INTENT(OUT)            :: alb1  ! new albedo in visible SW interval
 !    REAL, DIMENSION(klon), INTENT(OUT)            :: alb2  ! new albedo in near IR interval
-    REAL, DIMENSION(6), INTENT(IN)              ::SFRWL
-    REAL, DIMENSION(klon,nsw), INTENT(OUT)        ::alb_dir,alb_dif
+    REAL, DIMENSION(6), INTENT(IN)                :: SFRWL
+    REAL, DIMENSION(klon,nsw), INTENT(OUT)        :: alb_dir,alb_dif
 !albedo SB <<<
     REAL, DIMENSION(klon), INTENT(OUT)            :: evap, fluxsens, fluxlat
@@ -108,16 +113,18 @@
     REAL, DIMENSION(klon)    :: zfra, alb_neig
     REAL, DIMENSION(klon)    :: radsol
-    REAL, DIMENSION(klon)    :: u0, v0, u1_lay, v1_lay
-    INTEGER                  :: i,j
+    REAL, DIMENSION(klon)    :: u0, v0, u1_lay, v1_lay, ustar
+    INTEGER                  :: i,j,nt
 
     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)    :: bl_height, wind_velo      !height boundary layer, wind spd
+    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
     REAL, DIMENSION(klon)    :: pexner                    !Exner potential
     REAL                     :: pref
-    REAL, DIMENSION(klon,nsoilmx) :: tsoil0 !modfi
+    REAL, DIMENSION(klon,nsoilmx) :: tsoil0               !modif
+    REAL                          :: dtis                ! subtimestep
+    LOGICAL                       :: debut_is, lafin_is  ! debut and lafin for inlandsis
 
     CHARACTER (len = 20)                      :: modname = 'surf_landice'
@@ -139,7 +146,8 @@
   LOGICAL, SAVE :: firstcall = .TRUE.
   !$OMP THREADPRIVATE(firstcall)
-!FC
-
-
+
+
+!FC firtscall initializations
+!******************************************************************************************
   IF (firstcall) THEN
   alb_vis_sno_lic=0.77
@@ -149,7 +157,11 @@
   CALL getin_p('alb_nir_sno_lic',alb_nir_sno_lic)
            PRINT*, 'alb_nir_sno_lic',alb_nir_sno_lic
+  
+!  z0m=1.e-3
+!  z0h = z0m
   firstcall=.false.
   ENDIF
-!
+!******************************************************************************************
+
 ! Initialize output variables
     alb3(:) = 999999.
@@ -166,9 +178,16 @@
 
 !****************************************************************************************
-!   ok_snow = TRUE  : prepare and call SISVAT snow model
-!   ok_snow = FALSE : soil_model, calcul_flux, fonte_neige, ... 
-!
-!****************************************************************************************
-    IF (ok_snow) THEN
+!  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
+!****************************************************************************************
+
+
+    IF (landice_opt .EQ. 1) THEN
+    
+!****************************************************************************************
+! CALL to SISVAT interface
+!****************************************************************************************
+
 #ifdef CPP_SISVAT
        ! Prepare for calling SISVAT
@@ -192,11 +211,8 @@
           pexner(i)        = (p1lay(i)/pref)**(RD/RCPD)
           dens_air(i)      = p1lay(i)/RD/temp_air(i)  ! dry air density
-          bl_height(i)     = pphi1(i)/RG              
+          zsl_height(i)     = pphi1(i)/RG              
        END DO
 
-!****************************************************************************************
-! CALL to SISVAT interface
-!
-!****************************************************************************************
+
        ! config: compute everything with SV but temperatures afterwards with soil/calculfluxs
        DO i = 1, knon
@@ -209,5 +225,5 @@
             rmu0, swdown, lwdown, pexner, ps, p1lay, &
             precip_rain, precip_snow, precip_snow_adv, snow_adv, &
-            bl_height, wind_velo, temp_air, dens_air, spechum, tsurf, &
+            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, &
@@ -232,8 +248,87 @@
             flux_u1, flux_v1)
 #else
-       abort_message='Pb de coherence: ok_snow = .true. mais CPP_SISVAT = .false.'
+       abort_message='Pb de coherence: landice_opt = 1 mais CPP_SISVAT = .false.'
        CALL abort_physic(modname,abort_message,1)
 #endif
-    ELSE ! ok_snow=FALSE
+
+!****************************************************************************************
+! CALL to INLANDSIS interface
+!****************************************************************************************
+
+    ELSE IF (landice_opt .EQ. 2) THEN
+#ifdef CPP_INLANDSIS
+
+        debut_is=debut
+        lafin_is=.false.
+        ! Suppose zero surface speed
+        u0(:)            = 0.0
+        v0(:)            = 0.0
+
+
+        CALL calcul_flux_wind(knon, dtime, &
+         u0, v0, u1, v1, gustiness, cdragm, &
+         AcoefU, AcoefV, BcoefU, BcoefV, &
+         p1lay, temp_air, &
+         flux_u1, flux_v1)
+
+       
+       ! Set constants and compute some input for SISVAT
+       ! = 1000 hPa
+       ! and calculate incoming flux for SW and LW interval: swdown, lwdown
+       swdown(:)        = 0.0
+       lwdown(:)        = 0.0
+       snow_adv(:)      = 0.  ! no snow blown in for now
+       snow_cont_air(:) = 0.       
+       alb_soil(:)      = 0.4 ! before albedo(:) but here it is the ice albedo that we have to set
+       ustar(:)         = 0.
+       pref             = 100000.       
+       DO i = 1, knon
+          swdown(i)        = swnet(i)/(1-albedo(i))
+          lwdown(i)        = lwdownm(i)
+          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      
+          tsoil0(i,:)      = tsoil(i,:)  
+          ustar(i)= (cdragm(i)*(wind_velo(i)**2))**0.5   
+       END DO
+       
+
+       ! Subtimestepping
+       
+        dtis=dtime/n_dtis
+
+        DO nt=1,n_dtis
+
+          IF (lafin .and. nt.eq.n_dtis) 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, &
+            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
+
+
+#else
+       abort_message='Pb de coherence: landice_opt = 2 mais CPP_INLANDSIS = .false.'
+       CALL abort_physic(modname,abort_message,1)
+#endif
+
+
+
+    ELSE 
 
 !****************************************************************************************
@@ -241,4 +336,10 @@
 ! 
 !****************************************************************************************
+
+    ! EV: use calbeta
+    CALL calbeta(dtime, is_lic, knon, snow, qsol, beta, cal, dif_grnd)
+
+
+    ! use soil model and recalculate properly cal
     IF (soil_model) THEN 
        CALL soil(dtime, is_lic, knon, snow, tsurf, tsoil, soilcap, soilflux)
@@ -255,6 +356,6 @@
 !
 !****************************************************************************************
-    beta(:) = 1.0
-    dif_grnd(:) = 0.0
+!    beta(:) = 1.0
+!    dif_grnd(:) = 0.0
 
 ! Suppose zero surface speed
@@ -281,5 +382,5 @@
 !    
 !****************************************************************************************
-    CALL fonte_neige( knon, is_lic, knindex, dtime, &
+    CALL fonte_neige(knon, is_lic, knindex, dtime, &
          tsurf, precip_rain, precip_snow, &
          snow, qsol, tsurf_new, evap)
@@ -291,4 +392,5 @@
 !****************************************************************************************
     CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))  
+
     WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.
     zfra(1:knon) = MAX(0.0,MIN(1.0,snow(1:knon)/(snow(1:knon)+10.0)))
@@ -303,5 +405,5 @@
 !IM: KstaTER0.77 & LMD_ARMIP6    
 
-! Attantion: alb1 and alb2 are the same!
+! Attantion: alb1 and alb2 are not the same!
     alb1(1:knon)  = alb_vis_sno_lic
     alb2(1:knon)  = alb_nir_sno_lic
@@ -316,5 +418,10 @@
     z0m = SQRT(z0m**2+rugoro**2)
 
-    END IF ! ok_snow
+
+
+
+    
+
+    END IF ! landice_opt
 
 
@@ -333,12 +440,13 @@
   
 !****************************************************************************************
-       snow_o=0.
-       zfra_o = 0.
-       DO j = 1, knon
-           i = knindex(j)
-           snow_o(i) = snow(j) 
-           zfra_o(i) = zfra(j)
-       ENDDO
-
+! Etienne: comment these lines because of duplication just below
+!       snow_o=0.
+!       zfra_o = 0.
+!       DO j = 1, knon
+!           i = knindex(j)
+!           snow_o(i) = snow(j) 
+!           zfra_o(i) = zfra(j)
+!       ENDDO
+!
 !****************************************************************************************
        snow_o=0.
