Index: LMDZ4/trunk/libf/phylmd/interface_surf.F90
===================================================================
--- LMDZ4/trunk/libf/phylmd/interface_surf.F90	(revision 644)
+++ LMDZ4/trunk/libf/phylmd/interface_surf.F90	(revision 686)
@@ -39,5 +39,6 @@
 #include "YOMCST.inc"
 #include "indicesol.inc"
-
+!IM
+#include "clesphys.inc"
 
 ! run_off      ruissellement total
@@ -185,8 +186,10 @@
   real, dimension(klon), intent(IN) :: tsurf, p1lay
 !IM: "slab" ocean
+  real :: amn, amx 
   real, dimension(klon), intent(INOUT) :: tslab
   real, allocatable, dimension(:), save :: tmp_tslab
   real, dimension(klon), intent(OUT) :: flux_o, flux_g
   real, dimension(klon), intent(INOUT)        :: seaice ! glace de mer (kg/m2)
+  real, dimension(klon)                       :: siceh  ! hauteur glace de mer (m)
   REAL, DIMENSION(klon), INTENT(INOUT) :: radsol,fder
   real, dimension(klon), intent(IN) :: zmasq
@@ -255,4 +258,7 @@
   real, dimension(klon):: fder_prev
   REAL, dimension(klon) :: bidule
+!
+!IM ?? quelques variables pour netcdf
+#include "netcdf.inc"
 
   if (check) write(*,*) 'Entree ', modname
@@ -339,7 +345,4 @@
       endif
     endif
-    DO i=1, knon
-     tmp_radsol(knindex(i))=radsol(i)
-    ENDDO
     if (.not. allocated(tmp_pctsrf_slab)) then
       allocate(tmp_pctsrf_slab(klon,nbsrf), stat = error)
@@ -359,7 +362,7 @@
         call abort_gcm(modname,abort_message,1)
       endif
-    DO i=1, klon
-     tmp_seaice(i)=seaice(i)
-    ENDDO
+    IF(check) THEN
+     PRINT*,'allocation tmp_seaice nisurf itime',nisurf, itime
+    ENDIF
     endif
 !
@@ -515,7 +518,5 @@
   else if (nisurf == is_oce) then
 
-    if (check) write(*,*)'ocean, nisurf = ',nisurf 
-
-
+    if (check) write(*,*)'ocean, nisurf = ',nisurf,'knon=',knon 
 !
 ! Surface "ocean" appel a l'interface avec l'ocean
@@ -555,5 +556,7 @@
 !IM: "slab" ocean
      else if (ocean == 'slab  ') then
-      tsurf_new = tsurf
+      DO i=1, knon
+       tsurf_new(i) = tmp_tslab(knindex(i))
+      ENDDO
       pctsrf_new = tmp_pctsrf_slab
 !
@@ -599,5 +602,6 @@
      zx_sl(i) = RLVTT
      if (tsurf_new(i) .LT. RTT) zx_sl(i) = RLSTT
-     flux_o(i) = fluxsens(i)-evap(i)*zx_sl(i)
+!IM     flux_o(i) = fluxsens(i)-evap(i)*zx_sl(i)
+     flux_o(i) = fluxsens(i) + fluxlat(i)
      tmp_flux_o(knindex(i)) = flux_o(i)
      tmp_radsol(knindex(i))=radsol(i)
@@ -623,18 +627,4 @@
       & pctsrf_new)
 
-!IM: "slab" ocean 
-    else if (ocean == 'slab  ') then
-!
-       seaice=tmp_seaice
-       cumul = .true.
-       call interfoce(klon, debut, itime, dtime, jour, &
-     & tmp_radsol, tmp_flux_o, tmp_flux_g, pctsrf, &
-     & tslab, seaice, pctsrf_new)
-!
-       tmp_pctsrf_slab=pctsrf_new
-       DO i=1, knon
-        tsurf_new(i)=tslab(knindex(i))
-       ENDDO !i
-!
     endif
 
@@ -658,6 +648,5 @@
   else if (nisurf == is_sic) then
 
-    if (check) write(*,*)'sea ice, nisurf = ',nisurf 
-
+    if (check) write(*,*)'sea ice, nisurf = ',nisurf,'knon=',knon 
 !
 ! Surface "glace de mer" appel a l'interface avec l'ocean
@@ -697,5 +686,19 @@
 !IM: "slab" ocean
     else if (ocean == 'slab  ') then
+!IM ajout sicOBSERVE BEG
+   IF ( ok_slab_sicOBS) THEN 
+!                              ! lecture conditions limites
+      CALL interfoce(itime, dtime, jour, & 
+             &  klon, nisurf, knon, knindex, &
+             &  debut, &
+             &  tsurf_new, pctsrf_new)
+!
+    tmp_pctsrf_slab=pctsrf_new
+    print*,'jour lecture pctsrf_new sic =',jour
+!
+   ELSE !ok_slab_sicOBS 
      pctsrf_new=tmp_pctsrf_slab
+   ENDIF
+!IM ajout sicOBSERVE END
 !
      DO ii = 1, knon
@@ -732,5 +735,4 @@
       DO ii = 1, knon
        tsurf_new(ii) = tsurf(ii)
-!IMbad IF (pctsrf_new(ii,nisurf) < EPSFRA) then
        IF (pctsrf_new(knindex(ii),nisurf) < EPSFRA) then
           snow(ii) = 0.0
@@ -748,5 +750,5 @@
          cal(1:knon) = RCPD / soilcap(1:knon)
          radsol(1:knon) = radsol(1:knon)  + soilflux(1:knon)
-         dif_grnd = 0.
+         dif_grnd = 1.0 / tau_gl
       ELSE 
          dif_grnd = 1.0 / tau_gl
@@ -757,5 +759,5 @@
       tsurf_temp = tsurf_new
       beta = 1.0
-    ENDIF
+    ENDIF !ocean ==
 
     CALL calcul_fluxs( klon, knon, nisurf, dtime, &
@@ -766,27 +768,4 @@
          &   tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
 !
-!IM: flux entre l'ocean et la glace de mer pour le "slab" ocean
-    DO i = 1, knon
-     flux_g(i) = 0.0
-!
-!IM: faire dependre le coefficient de conduction de la glace de mer
-!    de l'epaisseur de la glace de mer, dans l'hypothese ou le coeff.
-!    actuel correspond a 3m de glace de mer, cf. L.Li
-!
-!      IF(1.EQ.0) THEN
-!       IF(siceh(i).GT.0.) THEN
-!        new_dif_grnd(i) = dif_grnd(i)*3./siceh(i)
-!       ELSE
-!        new_dif_grnd(i) = 0.
-!       ENDIF
-!      ENDIF !(1.EQ.0) THEN
-!
-     IF (cal(i).GT.1.0e-15) flux_g(i)=(tsurf_new(i)-t_grnd) &
-    &                          * dif_grnd(i) *RCPD/cal(i)
-!   &                          * new_dif_grnd(i) *RCPD/cal(i)
-      tmp_flux_g(knindex(i))=flux_g(i)
-      tmp_radsol(knindex(i))=radsol(i)
-    ENDDO
-
     IF (ocean /= 'couple') THEN
       CALL fonte_neige( klon, knon, nisurf, dtime, &
@@ -807,4 +786,33 @@
 !!      alb_new(1 : knon) = 0.6
     ENDIF
+!
+!IM: flux entre l'ocean et la glace de mer pour le "slab" ocean
+!
+    DO i = 1, knon
+!
+!IM: faire dependre le coefficient de conduction de la glace de mer
+!    de l'epaisseur de la glace de mer, dans l'hypothese ou le coeff.
+!    actuel correspond a 3m de glace de mer, cf. L.Li
+!
+       IF(1.EQ.0) THEN
+        IF(siceh(i).GT.0.) THEN
+         new_dif_grnd(i) = dif_grnd(i)*3./siceh(i)
+        ELSE
+         new_dif_grnd(i) = 0.
+        ENDIF
+       ENDIF !(1.EQ.0) THEN
+!
+     IF (cal(i).GT.1.0e-15) THEN
+       flux_g(i)=(tsurf_new(i)-t_grnd) &
+    &                          * dif_grnd(i) *RCPD/cal(i)
+!   &                          * new_dif_grnd(i) *RCPD/cal(i)
+     ENDIF
+     tmp_flux_g(knindex(i))=flux_g(i)
+!
+!IM: Attention: ne pas initialiser le tmp_radsol puisque c'est deja fait sur is_oce;
+!IM:            tmp_radsol doit etre le flux solaire qui arrive sur l'ocean 
+!IM:            et non pas celui qui arrive sur la glace de mer
+!
+    ENDDO
 
     fder_prev = fder    
@@ -841,6 +849,30 @@
       & pctsrf_new)
 
-!    else if (ocean == 'slab  ') then
-!      call interfoce(nisurf)
+!IM: "slab" ocean 
+    else if (ocean == 'slab  ') then
+!
+   IF (check) THEN
+    amn=MIN(tmp_tslab(1),1000.)
+    amx=MAX(tmp_tslab(1),-1000.)
+    DO i=2, klon
+     amn=MIN(tmp_tslab(i),amn)
+     amx=MAX(tmp_tslab(i),amx)
+    ENDDO
+!
+    PRINT*,' debut avant interfoce_slab min max tmp_tslab',amn,amx
+   ENDIF !(check) THEN
+!
+       cumul = .true.
+       tslab = tmp_tslab   
+       call interfoce(klon, debut, itime, dtime, jour, &
+     & tmp_radsol, tmp_flux_o, tmp_flux_g, tmp_pctsrf_slab, &
+     & tslab, seaice, pctsrf_new)
+!
+       tmp_seaice=seaice
+       tmp_pctsrf_slab=pctsrf_new
+       DO i=1, knon
+         tmp_tslab(knindex(i))=tslab(knindex(i))
+       ENDDO !i
+!
 
     endif
@@ -850,4 +882,5 @@
     z0_new = SQRT(z0_new**2+rugoro**2)
     alblw(1:knon) = alb_new(1:knon)
+
 
   else if (nisurf == is_lic) then
@@ -2065,5 +2098,6 @@
 ! la glace de mer pour un "slab" ocean de 50m
 !
-! I. Musat 04.02.2005
+! Conception: Laurent Li
+! Re-ecriture + adaptation LMDZ4: I. Musat
 !
 ! input:
@@ -2100,4 +2134,5 @@
 !
 ! Variables locales :
+  REAL :: amn, amx
   INTEGER, save :: lmt_pas, julien, idayvrai
   REAL, parameter :: unjour=86400.
@@ -2149,4 +2184,20 @@
 !
    ENDIF !debut
+!
+   IF (check ) THEN
+    amn=MIN(tmp_tslab(1),1000.)
+    amx=MAX(tmp_tslab(1),-1000.)
+    DO i=2, klon
+     amn=MIN(tmp_tslab(i),amn)
+     amx=MAX(tmp_tslab(i),amx)
+    ENDDO
+!
+    PRINT*,' debut min max tslab',amn,amx 
+!
+!!
+    PRINT*,' itap,lmt_pas unjour',itap,lmt_pas,unjour
+   ENDIF 
+!!
+!
    pctsrf_slab(1:klon,1:nbsrf) = pctsrf(1:klon,1:nbsrf)
 !
@@ -2198,8 +2249,13 @@
 ! et pctsrf(i,is_sic) croit lineairement avec seaice de 0. a 20cm d'epaisseur
 !
+
+     IF(.NOT.ok_slab_sicOBS) then
       pctsrf_slab(i,is_sic)=MIN(siceh(i)/0.20, &
      &                      1.-(pctsrf_slab(i,is_ter)+pctsrf_slab(i,is_lic)))
       pctsrf_slab(i,is_oce)=1.0 - &
      &      (pctsrf_slab(i,is_ter)+pctsrf_slab(i,is_lic)+pctsrf_slab(i,is_sic))
+     ELSE
+      IF (i.EQ.1) print*,'cas ok_slab_sicOBS TRUE : passe sur la modif.'
+     ENDIF !(.NOT.ok_slab_sicOBS) then
      ENDIF !pctsrf
    ENDDO
@@ -2220,4 +2276,8 @@
 !
    IF (MOD(itap,lmt_pas).EQ.0) THEN !fin de journee
+!
+! calcul tslab 
+    amn=MIN(tmp_tslab(1),1000.)
+    amx=MAX(tmp_tslab(1),-1000.)
     DO i = 1, klon
       IF ((pctsrf_slab(i,is_oce).GT.epsfra).OR. &
@@ -2229,6 +2289,20 @@
        slab_bils(i) = 0.
       ENDIF !pctsrf
+!
+      IF (check) THEN
+       IF(i.EQ.1) THEN  
+        PRINT*,'i tmp_tslab MOD(itap,lmt_pas).EQ.0 sicINTER',i,itap, &
+      & tmp_tslab(i), tmp_seaice(i)
+       ENDIF
+!
+       amn=MIN(tmp_tslab(i),amn)
+       amx=MAX(tmp_tslab(i),amx)
+      ENDIF
     ENDDO !klon
    ENDIF !(MOD(itap,lmt_pas).EQ.0) THEN
+!
+   IF ( check ) THEN
+    PRINT*,'fin min max tslab',amn,amx 
+   ENDIF 
 !
    tslab = tmp_tslab
