Index: LMDZ5/trunk/libf/phylmd/ocean_slab_mod.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/ocean_slab_mod.F90	(revision 1907)
+++ LMDZ5/trunk/libf/phylmd/ocean_slab_mod.F90	(revision 2057)
@@ -5,7 +5,26 @@
 ! "ocean=slab".
 !
+
+  USE dimphy
+  USE indice_sol_mod
+
   IMPLICIT NONE
   PRIVATE
-  PUBLIC :: ocean_slab_frac, ocean_slab_noice
+  PUBLIC :: ocean_slab_init, ocean_slab_frac, ocean_slab_noice!, ocean_slab_ice
+
+  INTEGER, PRIVATE, SAVE                           :: cpl_pas
+  !$OMP THREADPRIVATE(cpl_pas)
+  REAL, PRIVATE, SAVE                              :: cyang
+  !$OMP THREADPRIVATE(cyang)
+  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE   :: slabh
+  !$OMP THREADPRIVATE(slabh)
+  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: tslab
+  !$OMP THREADPRIVATE(tslab)
+  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE  :: pctsrf
+  !$OMP THREADPRIVATE(pctsrf)
+  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: slab_bils
+  !$OMP THREADPRIVATE(slab_bils)
+  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: bils_cum
+  !$OMP THREADPRIVATE(bils_cum)
 
 CONTAINS
@@ -13,10 +32,74 @@
 !****************************************************************************************
 !
-  SUBROUTINE ocean_slab_frac(itime, dtime, jour, pctsrf, is_modified)
-
-    USE dimphy
+  SUBROUTINE ocean_slab_init(dtime, pctsrf_rst)
+  !, seaice_rst etc
+
+    use IOIPSL
+
+    INCLUDE "iniprint.h"
+    ! For ok_xxx vars (Ekman...)
+    INCLUDE "clesphys.h"
+
+    ! Input variables
+!****************************************************************************************
+    REAL, INTENT(IN)                         :: dtime
+! Variables read from restart file
+    REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf_rst
+
+! Local variables
+!****************************************************************************************
+    INTEGER                :: error
+    CHARACTER (len = 80)   :: abort_message
+    CHARACTER (len = 20)   :: modname = 'ocean_slab_intit'
+
+!****************************************************************************************
+! Allocate surface fraction read from restart file
+!****************************************************************************************
+    ALLOCATE(pctsrf(klon,nbsrf), stat = error)
+    IF (error /= 0) THEN
+       abort_message='Pb allocation tmp_pctsrf_slab'
+       CALL abort_gcm(modname,abort_message,1)
+    ENDIF
+    pctsrf(:,:) = pctsrf_rst(:,:)
+
+!****************************************************************************************
+! Allocate local variables
+!****************************************************************************************
+    ALLOCATE(slab_bils(klon), stat = error)
+    IF (error /= 0) THEN
+       abort_message='Pb allocation slab_bils'
+       CALL abort_gcm(modname,abort_message,1)
+    ENDIF
+    slab_bils(:) = 0.0   
+    ALLOCATE(bils_cum(klon), stat = error)
+    IF (error /= 0) THEN
+       abort_message='Pb allocation slab_bils_cum'
+       CALL abort_gcm(modname,abort_message,1)
+    ENDIF
+    bils_cum(:) = 0.0   
+
+! Layer thickness
+    ALLOCATE(slabh(nslay), stat = error)
+    IF (error /= 0) THEN
+       abort_message='Pb allocation slabh'
+       CALL abort_gcm(modname,abort_message,1)
+    ENDIF
+    slabh(1)=50.
+! cyang = 1/heat capacity of top layer (rho.c.H)
+    cyang=1/(slabh(1)*4.228e+06)
+
+! cpl_pas periode de couplage avec slab (update tslab, pctsrf)
+! pour un calcul à chaque pas de temps, cpl_pas=1
+    cpl_pas = NINT(86400./dtime * 1.0) ! une fois par jour
+    CALL getin('cpl_pas',cpl_pas)
+    print *,'cpl_pas',cpl_pas
+ END SUBROUTINE ocean_slab_init
+!
+!****************************************************************************************
+!
+  SUBROUTINE ocean_slab_frac(itime, dtime, jour, pctsrf_chg, is_modified)
+
     USE limit_read_mod
     USE surface_data
-    USE indice_sol_mod
 
 !    INCLUDE "clesphys.h"
@@ -27,20 +110,15 @@
     INTEGER, INTENT(IN)                        :: jour    ! jour a lire dans l'annee
     REAL   , INTENT(IN)                        :: dtime   ! pas de temps de la physique (en s)
-    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: pctsrf  ! sub-surface fraction
+    REAL, DIMENSION(klon,nbsrf), INTENT(INOUT) :: pctsrf_chg  ! sub-surface fraction
     LOGICAL, INTENT(OUT)                       :: is_modified ! true if pctsrf is modified at this time step
 
 ! Local variables
 !****************************************************************************************
-    CHARACTER (len = 80)   :: abort_message
-    CHARACTER (len = 20)   :: modname = 'ocean_slab_frac'
-
-
-    IF (version_ocean == 'sicOBS') THEN   
-       CALL limit_read_frac(itime, dtime, jour, pctsrf, is_modified)
+
+    IF (version_ocean == 'sicOBS'.OR. version_ocean == 'sicNO') THEN   
+       CALL limit_read_frac(itime, dtime, jour, pctsrf_chg, is_modified)
     ELSE
-       abort_message='Ocean slab model without forced sea-ice fractions has to be rewritten!!!'
-       CALL abort_gcm(modname,abort_message,1)
-! Here should sea-ice/ocean fraction either be calculated or returned if saved as a module varaiable 
-! (in the case the new fractions are calculated in ocean_slab_ice or ocean_slab_noice subroutines).  
+       pctsrf_chg(:,:)=pctsrf(:,:)
+       is_modified=.TRUE.
     END IF
 
@@ -57,9 +135,8 @@
        radsol, snow, agesno, &
        qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
-       tsurf_new, dflux_s, dflux_l, lmt_bils)
+       tsurf_new, dflux_s, dflux_l, qflux)
     
-    USE dimphy
     USE calcul_fluxs_mod
-    USE indice_sol_mod
+    USE surface_data
 
     INCLUDE "iniprint.h"
@@ -95,15 +172,13 @@
     REAL, DIMENSION(klon), INTENT(OUT)   :: tsurf_new
     REAL, DIMENSION(klon), INTENT(OUT)   :: dflux_s, dflux_l      
-    REAL, DIMENSION(klon), INTENT(OUT)   :: lmt_bils
+    REAL, DIMENSION(klon), INTENT(OUT)   :: qflux
 
 ! Local variables
 !****************************************************************************************
-    INTEGER               :: i
+    INTEGER               :: i,ki
     REAL, DIMENSION(klon) :: cal, beta, dif_grnd
-    REAL, DIMENSION(klon) :: lmt_bils_oce, lmt_foce, diff_sst
+    REAL, DIMENSION(klon) :: diff_sst, lmt_bils
     REAL, DIMENSION(klon) :: u0, v0
     REAL, DIMENSION(klon) :: u1_lay, v1_lay
-    REAL                  :: calc_bils_oce, deltat
-    REAL, PARAMETER       :: cyang=50.0 * 4.228e+06 ! capacite calorifique volumetrique de l'eau J/(m2 K)
 
 !****************************************************************************************
@@ -136,28 +211,124 @@
          flux_u1, flux_v1)  
 
-!****************************************************************************************
-! 2) Get global variables lmt_bils and lmt_foce from file limit_slab.nc
-!
-!****************************************************************************************
-    CALL limit_slab(itime, dtime, jour, lmt_bils, lmt_foce, diff_sst)  ! global pour un processus
-
-    lmt_bils_oce(:) = 0.
-    WHERE (lmt_foce > 0.) 
-       lmt_bils_oce = lmt_bils / lmt_foce ! global 
-    END WHERE
+! Accumulate total fluxes locally
+    slab_bils(:)=0.
+    DO i=1,knon
+        ki=knindex(i)
+        slab_bils(ki)=(fluxlat(i)+fluxsens(i)+radsol(i))*pctsrf(ki,is_oce)/(1.-zmasq(ki))
+        bils_cum(ki)=bils_cum(ki)+slab_bils(ki)
+! Also taux, tauy, saved vars...
+    END DO
+
+!****************************************************************************************
+! 2) Get global variables lmt_bils and diff_sst from file limit_slab.nc
+!
+!****************************************************************************************
+    lmt_bils(:)=0.
+    CALL limit_slab(itime, dtime, jour, lmt_bils, diff_sst)  ! global pour un processus
+    ! lmt_bils and diff_sst saved by limit_slab
+    qflux(:)=lmt_bils(:)+diff_sst(:)/cyang/86400.
+    ! qflux = total QFlux correction (in W/m2)
 
 !****************************************************************************************
 ! 3) Recalculate new temperature
 !
-!****************************************************************************************
-    DO i = 1, knon
-       calc_bils_oce = radsol(i) + fluxsens(i) + fluxlat(i)
-       deltat        = (calc_bils_oce - lmt_bils_oce(knindex(i)))*dtime/cyang +diff_sst(knindex(i))
-       tsurf_new(i)  = tsurf_in(i) + deltat
-    END DO
-
+!***********************************************o*****************************************
+    tsurf_new=tsurf_in
+    IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab & fraction
+        ! Compute transport
+        ! Add read QFlux and SST tendency
+        tslab(:,1)=tslab(:,1)+qflux(:)*cyang*dtime*cpl_pas 
+        ! Add cumulated surface fluxes
+        tslab(:,1)=tslab(:,1)+bils_cum(:)*cyang*dtime
+        ! Update surface temperature
+        SELECT CASE(version_ocean)
+        CASE('sicNO')
+            DO i=1,knon
+                ki=knindex(i)
+                tsurf_new(i)=tslab(ki,1)
+            END DO
+        CASE('sicOBS') ! check for sea ice or tsurf below freezing
+            DO i=1,knon
+                ki=knindex(i)
+                IF ((tslab(ki,1).LT.t_freeze).OR.(pctsrf(ki,is_sic).GT.epsfra)) THEN
+                    tsurf_new(i)=t_freeze
+                    tslab(ki,1)=t_freeze
+                ELSE
+                    tsurf_new(i)=tslab(ki,1)
+                END IF
+            END DO
+        CASE('sicINT')
+            DO i=1,knon
+                ki=knindex(i)
+                IF (pctsrf(ki,is_sic).LT.epsfra) THEN ! Free of ice
+                    IF (tslab(ki,1).GT.t_freeze) THEN 
+                        tsurf_new(i)=tslab(ki,1)
+                    ELSE
+                        tsurf_new(i)=t_freeze
+                        ! Call new ice routine
+                        tslab(ki,1)=t_freeze
+                    END IF
+                ELSE ! ice present, tslab update completed in slab_ice
+                    tsurf_new(i)=t_freeze
+                END IF !ice free
+            END DO
+        END SELECT
+        bils_cum(:)=0.0! clear cumulated fluxes
+    END IF ! coupling time
   END SUBROUTINE ocean_slab_noice
 !
 !****************************************************************************************
 !
+!  SUBROUTINE ocean_slab_ice(   &
+!       itime, dtime, jour, knon, knindex, &
+!       tsurf_in, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
+!       AcoefH, AcoefQ, BcoefH, BcoefQ, &
+!       AcoefU, AcoefV, BcoefU, BcoefV, &
+!       ps, u1, v1, &
+!       radsol, snow, qsurf, qsol, agesno, tsoil, &
+!       alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
+!       tsurf_new, dflux_s, dflux_l)
+!
+!****************************************************************************************
+! 1) Flux calculation
+!****************************************************************************************
+! set beta, cal etc. depends snow / ice surf ?
+! calcul_fluxs (sens, lat etc)
+! calcul_flux_wind
+
+!****************************************************************************************
+! 2) Update surface
+!****************************************************************************************
+! neige, fonte
+! flux glace-ocean
+! update temperature
+! neige precip, evap
+! Melt snow & ice from above
+! New albedo
+
+!****************************************************************************************
+! 3) Recalculate new ocean temperature
+!    Melt / freeze from below
+!***********************************************o*****************************************
+
+
+!  END SUBROUTINE ocean_slab_ice
+!
+!****************************************************************************************
+!
+  SUBROUTINE ocean_slab_final
+  !, seaice_rst etc
+
+    ! For ok_xxx vars (Ekman...)
+    INCLUDE "clesphys.h"
+
+!****************************************************************************************
+! Deallocate module variables
+!
+!****************************************************************************************
+    IF (ALLOCATED(pctsrf)) DEALLOCATE(pctsrf)
+    IF (ALLOCATED(tslab)) DEALLOCATE(tslab)
+
+  END SUBROUTINE ocean_slab_final
+
 END MODULE ocean_slab_mod
