!
! $Header$
!
MODULE ocean_slab_mod
!
! This module is used for both surface ocean and sea-ice when using the slab ocean,
! "ocean=slab".
!
  USE surface_data,     ONLY : tau_gl, calice, calsno
  USE fonte_neige_mod,  ONLY : fonte_neige
  USE calcul_fluxs_mod, ONLY : calcul_fluxs
  USE dimphy
  
  IMPLICIT NONE

  INTEGER, PRIVATE, SAVE                           :: lmt_pas, julien, idayvrai
  !$OMP THREADPRIVATE(lmt_pas,julien,idayvrai)
  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE   :: tmp_seaice
  !$OMP THREADPRIVATE(tmp_seaice)
  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE   :: tmp_tslab_loc
  !$OMP THREADPRIVATE(tmp_tslab_loc)
  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE   :: slab_bils
  !$OMP THREADPRIVATE(slab_bils)
  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE , SAVE  :: lmt_bils
  !$OMP THREADPRIVATE(lmt_bils)
  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: tmp_pctsrf_slab
  !$OMP THREADPRIVATE(tmp_pctsrf_slab)
  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE   :: tmp_tslab
  !$OMP THREADPRIVATE(tmp_tslab)
  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE   :: tmp_radsol
  !$OMP THREADPRIVATE(tmp_radsol)
  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE   :: tmp_flux_o, tmp_flux_g
  !$OMP THREADPRIVATE(tmp_flux_o,tmp_flux_g)
  LOGICAL, PRIVATE, SAVE                           :: check = .FALSE.
  !$OMP THREADPRIVATE(check)

CONTAINS
!
!****************************************************************************************
!
  SUBROUTINE ocean_slab_init(dtime, tslab_rst, seaice_rst, pctsrf_rst)

    INCLUDE "indicesol.h"
    INCLUDE "iniprint.h"

! Input variables
!****************************************************************************************
    REAL, INTENT(IN)                         :: dtime
! Variables read from restart file
    REAL, DIMENSION(klon), INTENT(IN)        :: tslab_rst          
    REAL, DIMENSION(klon), INTENT(IN)        :: seaice_rst
    REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf_rst


! Local variables
!****************************************************************************************
    INTEGER                :: error
    CHARACTER (len = 80)   :: abort_message
    CHARACTER (len = 20)   :: modname = 'ocean_slab_intit'


    WRITE(lunout,*) '************************'
    WRITE(lunout,*) 'SLAB OCEAN est actif, prenez precautions !'
    WRITE(lunout,*) '************************'    

! Allocate variables initialize from restart fields
    ALLOCATE(tmp_tslab(klon), stat = error)
    IF (error /= 0) THEN
       abort_message='Pb allocation tmp_tslab'
       CALL abort_gcm(modname,abort_message,1)
    ENDIF
    tmp_tslab(:) = tslab_rst(:)

    ALLOCATE(tmp_tslab_loc(klon), stat = error)
    IF (error /= 0) THEN
       abort_message='Pb allocation tmp_tslab_loc'
       CALL abort_gcm(modname,abort_message,1)
    ENDIF
    tmp_tslab_loc(:) = tslab_rst(:)

    ALLOCATE(tmp_seaice(klon), stat = error)
    IF (error /= 0) THEN
       abort_message='Pb allocation tmp_seaice'
       CALL abort_gcm(modname,abort_message,1)
    ENDIF
    tmp_seaice(:) = seaice_rst(:)

    ALLOCATE(tmp_pctsrf_slab(klon,nbsrf), stat = error)
    IF (error /= 0) THEN
       abort_message='Pb allocation tmp_pctsrf_slab'
       CALL abort_gcm(modname,abort_message,1)
    ENDIF
    tmp_pctsrf_slab(:,:) = pctsrf_rst(:,:)
    
! Allocate some other variables internal in module mod_oceanslab
    ALLOCATE(tmp_radsol(klon), stat = error)
    IF (error /= 0) THEN
       abort_message='Pb allocation tmp_radsol'
       CALL abort_gcm(modname,abort_message,1)
    ENDIF

    ALLOCATE(tmp_flux_o(klon), stat = error)
    IF (error /= 0) THEN
       abort_message='Pb allocation tmp_flux_o'
       CALL abort_gcm(modname,abort_message,1)
    ENDIF
    
    ALLOCATE(tmp_flux_g(klon), stat = error)
    IF (error /= 0) THEN
       abort_message='Pb allocation tmp_flux_g'
       CALL abort_gcm(modname,abort_message,1)
    ENDIF

! a mettre un slab_bils aussi en force !!!
    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(lmt_bils(klon), stat = error)
    IF (error /= 0) THEN
       abort_message='Pb allocation lmt_bils'
       CALL abort_gcm(modname,abort_message,1)
    ENDIF


! pour une lecture une fois par jour    
    lmt_pas = NINT(86400./dtime * 1.0) 

  END SUBROUTINE ocean_slab_init
!
!****************************************************************************************
!
  SUBROUTINE ocean_slab_noice( & 
       dtime, knon, knindex, &
       p1lay, tq_cdrag, precip_rain, precip_snow, temp_air, spechum, &
       petAcoef, peqAcoef, petBcoef, peqBcoef, &
       ps, u1_lay, v1_lay, &
       radsol, snow, agesno, &
       qsurf, evap, fluxsens, fluxlat, &
       tsurf_new, dflux_s, dflux_l, pctsrf_oce)

    INCLUDE "indicesol.h"
    INCLUDE "iniprint.h"

! Input arguments
!****************************************************************************************
    INTEGER, INTENT(IN)                  :: knon
    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
    REAL, INTENT(IN)                     :: dtime
    REAL, DIMENSION(klon), INTENT(IN)    :: p1lay
    REAL, DIMENSION(klon), INTENT(IN)    :: tq_cdrag
    REAL, DIMENSION(klon), INTENT(IN)    :: precip_rain, precip_snow
    REAL, DIMENSION(klon), INTENT(IN)    :: temp_air, spechum
    REAL, DIMENSION(klon), INTENT(IN)    :: petAcoef, peqAcoef
    REAL, DIMENSION(klon), INTENT(IN)    :: petBcoef, peqBcoef
    REAL, DIMENSION(klon), INTENT(IN)    :: ps
    REAL, DIMENSION(klon), INTENT(IN)    :: u1_lay, v1_lay

! In/Output arguments
!****************************************************************************************
    REAL, DIMENSION(klon), INTENT(INOUT) :: radsol
    REAL, DIMENSION(klon), INTENT(INOUT) :: snow
    REAL, DIMENSION(klon), INTENT(INOUT) :: agesno
    
! Output arguments
!****************************************************************************************
    REAL, DIMENSION(klon), INTENT(OUT)   :: qsurf
    REAL, DIMENSION(klon), INTENT(OUT)   :: evap, fluxsens, fluxlat
    REAL, DIMENSION(klon), INTENT(OUT)   :: tsurf_new
    REAL, DIMENSION(klon), INTENT(OUT)   :: dflux_s, dflux_l      
    REAL, DIMENSION(klon), INTENT(OUT)   :: pctsrf_oce

! Local variables
!****************************************************************************************
    INTEGER                :: i
    REAL, DIMENSION(klon)  :: cal, beta, dif_grnd
    REAL, DIMENSION(klon)  :: tsurf_temp

!****************************************************************************************
    IF (check) WRITE(*,*)' Entering ocean_slab_noice'    

    tsurf_new(1:knon) = tmp_tslab(knindex(1:knon))
    pctsrf_oce(:)   = tmp_pctsrf_slab(:,is_oce)
    
    tsurf_temp(:) = tsurf_new(:)
    cal = 0.
    beta = 1.
    dif_grnd = 0.
    agesno(:) = 0.
    
    CALL calcul_fluxs(knon, is_oce, dtime, &
         tsurf_temp, p1lay, cal, beta, tq_cdrag, ps, &
         precip_rain, precip_snow, snow, qsurf,  &
         radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
         petAcoef, peqAcoef, petBcoef, peqBcoef, &
         tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
    
    tmp_flux_o(:) = 0.0
    tmp_radsol(:) = 0.0

    DO i=1, knon
       tmp_radsol(knindex(i))=radsol(i)
       
       IF (pctsrf_oce(knindex(i)) .GT. epsfra) &
            tmp_flux_o(knindex(i)) = fluxsens(i) + fluxlat(i)
    ENDDO
    
  END SUBROUTINE ocean_slab_noice
!
!****************************************************************************************
!
  SUBROUTINE ocean_slab_ice(   &
       itime, dtime, jour, knon, knindex, &
       debut, &
       tsurf, p1lay, tq_cdrag, precip_rain, precip_snow, temp_air, spechum, &
       petAcoef, peqAcoef, petBcoef, peqBcoef, &
       ps, u1_lay, v1_lay, &
       radsol, snow, qsurf, qsol, agesno, tsoil, &
       alb1_new, alb2_new, evap, fluxsens, fluxlat, &
       tsurf_new, dflux_s, dflux_l, pctsrf_sic)

    INCLUDE "indicesol.h"
    INCLUDE "dimsoil.h"
    INCLUDE "YOMCST.h"
    INCLUDE "iniprint.h"
    INCLUDE "clesphys.h"

! Input arguments  
!****************************************************************************************
    INTEGER, INTENT(IN)                  :: itime, jour, knon
    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
    REAL, INTENT(IN)                     :: dtime
    REAL, DIMENSION(klon), INTENT(IN)    :: tsurf
    REAL, DIMENSION(klon), INTENT(IN)    :: p1lay
    REAL, DIMENSION(klon), INTENT(IN)    :: tq_cdrag
    REAL, DIMENSION(klon), INTENT(IN)    :: precip_rain, precip_snow
    REAL, DIMENSION(klon), INTENT(IN)    :: temp_air, spechum
    REAL, DIMENSION(klon), INTENT(IN)    :: petAcoef, peqAcoef
    REAL, DIMENSION(klon), INTENT(IN)    :: petBcoef, peqBcoef
    REAL, DIMENSION(klon), INTENT(IN)    :: ps
    REAL, DIMENSION(klon), INTENT(IN)    :: u1_lay, v1_lay
    LOGICAL, INTENT(IN)                  :: debut

!In/Output arguments
!****************************************************************************************
    REAL, DIMENSION(klon), INTENT(INOUT)          :: radsol
    REAL, DIMENSION(klon), INTENT(INOUT)          :: snow, qsol
    REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
    REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil

! Output arguments
!****************************************************************************************
    REAL, DIMENSION(klon), INTENT(OUT)            :: qsurf
    REAL, DIMENSION(klon), INTENT(OUT)            :: alb1_new  ! new albedo in visible SW interval
    REAL, DIMENSION(klon), INTENT(OUT)            :: alb2_new  ! new albedo in near IR interval
    REAL, DIMENSION(klon), INTENT(OUT)            :: evap, fluxsens, fluxlat
    REAL, DIMENSION(klon), INTENT(OUT)            :: tsurf_new
    REAL, DIMENSION(klon), INTENT(OUT)            :: dflux_s, dflux_l      
    REAL, DIMENSION(klon), INTENT(OUT)            :: pctsrf_sic

! Local variables
!****************************************************************************************
    INTEGER                              :: i
    REAL, DIMENSION(klon)                :: cal, beta, dif_grnd, capsol
    REAL, DIMENSION(klon)                :: alb_neig, tsurf_temp
    REAL, DIMENSION(klon)                :: soilcap, soilflux
    REAL, DIMENSION(klon)                :: zfra
    REAL, PARAMETER                      :: t_grnd=271.35
    REAL                                 :: amn, amx 
    REAL, DIMENSION(klon)                :: tslab
    REAL, DIMENSION(klon)                :: seaice ! glace de mer (kg/m2)
    REAL, DIMENSION(klon,nbsrf)          :: pctsrf_new

!****************************************************************************************

    IF (check) WRITE(*,*)'Entering surface_seaice, knon=',knon 

! Initialization of output variables
    alb1_new(:) = 0.0

!****************************************************************************************
! 
!
!****************************************************************************************
    IF ( ok_slab_sicOBS) THEN   
       ! glace de mer observee, lecture conditions limites
       CALL interfoce_lim(itime, dtime, jour, & 
            knon, knindex, &
            debut, &
            tsurf_new, pctsrf_new)

       tmp_pctsrf_slab(:,:) = pctsrf_new(:,:)
       WRITE(lunout,*) 'jour lecture pctsrf_new sic =',jour

    ELSE
       pctsrf_new=tmp_pctsrf_slab
    ENDIF

    DO i = 1, knon
       tsurf_new(i) = tsurf(i)
       IF (pctsrf_new(knindex(i),is_sic) < EPSFRA) THEN
          snow(i) = 0.0
          tsurf_new(i) = RTT - 1.8
          IF (soil_model) tsoil(i,:) = RTT -1.8
       ENDIF
    ENDDO
    
    CALL calbeta(dtime, is_sic, knon, snow, qsol, beta, capsol, dif_grnd)
    
    IF (soil_model) THEN 
       CALL soil(dtime, is_sic, knon, snow, tsurf_new, tsoil, soilcap, soilflux)
       cal(1:knon) = RCPD / soilcap(1:knon)
       radsol(1:knon) = radsol(1:knon)  + soilflux(1:knon)
    ELSE 
       dif_grnd = 1.0 / tau_gl
       cal = RCPD * calice
       WHERE (snow > 0.0) cal = RCPD * calsno 
    ENDIF
    tsurf_temp = tsurf_new
    beta = 1.0

!****************************************************************************************
!
!
!****************************************************************************************
    CALL calcul_fluxs(knon, is_sic, dtime, &
         tsurf_temp, p1lay, cal, beta, tq_cdrag, ps, &
         precip_rain, precip_snow, snow, qsurf,  &
         radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
         petAcoef, peqAcoef, petBcoef, peqBcoef, &
         tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
    
    CALL fonte_neige( knon, is_sic, knindex, dtime, &
         tsurf_temp, precip_rain, precip_snow, &
         snow, qsol, tsurf_new, evap)

!****************************************************************************************
!     calcul albedo
!
!****************************************************************************************
    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)))
    alb1_new(1 : knon) = alb_neig(1 : knon) *zfra(1:knon) + & 
         0.6 * (1.0-zfra(1:knon))
    
    alb2_new(:) = alb1_new(:)

!
!IM: flux entre l'ocean et la glace de mer pour le "slab" ocean
    tmp_flux_g(:) = 0.0
    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 ((cal(i).GT.1.0e-15) .AND. (pctsrf_new(knindex(i),is_sic) .GT. epsfra)) THEN
          tmp_flux_g(knindex(i))=(tsurf_new(i)-t_grnd) &
               * dif_grnd(i) *RCPD/cal(i)
       ENDIF
!
!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
    
!
! 2eme appel a interfoce pour le cumul et le passage des flux a l'ocean
!

    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
       
       WRITE(lunout,*) ' debut avant interfoce_slab min max tmp_tslab',amn,amx
    ENDIF !(check) THEN

    tslab = tmp_tslab   

    CALL interfoce_slab(klon, debut, itime, dtime, jour, &
         tmp_radsol, tmp_flux_o, tmp_flux_g, tmp_pctsrf_slab, &
         tslab, seaice, pctsrf_new)
    
    tmp_pctsrf_slab=pctsrf_new

    DO i=1, knon
       tmp_tslab(knindex(i))=tslab(knindex(i))
    ENDDO
     
        
!****************************************************************************************
! Return the fraction of sea-ice
! NB! jg : Peut-etre un probleme, faut-il prend aussi tmp_pctsrf_slab(:,is_oce)???
!****************************************************************************************
    pctsrf_sic(:) =  tmp_pctsrf_slab(:,is_sic)


  END SUBROUTINE ocean_slab_ice
! 
!****************************************************************************************
!
  SUBROUTINE interfoce_slab(klon, debut, itap, dtime, ijour, &
       radsol, fluxo, fluxg, pctsrf, &
       tslab, seaice, pctsrf_slab)
!
! Cette routine calcule la temperature d'un slab ocean, la glace de mer 
! et les pourcentages de la maille couverte par l'ocean libre et/ou 
! la glace de mer pour un "slab" ocean de 50m
!
! Conception: Laurent Li
! Re-ecriture + adaptation LMDZ4: I. Musat
!
! input:
!   klon         nombre total de points de grille
!   debut        logical: 1er appel a la physique
!   itap         numero du pas de temps
!   dtime        pas de temps de la physique (en s)
!   ijour        jour dans l'annee en cours
!   radsol       rayonnement net au sol (LW + SW)
!   fluxo        flux turbulent (sensible + latent) sur les mailles oceaniques 
!   fluxg        flux de conduction entre la surface de la glace de mer et l'ocean
!   pctsrf       tableau des pourcentages de surface de chaque maille
! output: 
!   tslab        temperature de l'ocean libre
!   seaice       glace de mer (kg/m2)
!   pctsrf_slab  "pourcentages" (valeurs entre 0. et 1.) surfaces issus du slab

    INCLUDE "indicesol.h"
    INCLUDE "YOMCST.h"
    INCLUDE "iniprint.h"
    INCLUDE "clesphys.h"

! Input arguments
!****************************************************************************************
    INTEGER, INTENT(IN)                       :: klon
    LOGICAL, INTENT(IN)                       :: debut    ! not used
    INTEGER, INTENT(IN)                       :: itap
    REAL, INTENT(IN)                          :: dtime       ! not used
    INTEGER, INTENT(IN)                       :: ijour
    REAL, DIMENSION(klon), INTENT(IN)         :: radsol
    REAL, DIMENSION(klon), INTENT(IN)         :: fluxo
    REAL, DIMENSION(klon), INTENT(IN)         :: fluxg
    REAL, DIMENSION(klon, nbsrf), INTENT(IN)  :: pctsrf

! Output arguments
!****************************************************************************************
    REAL, DIMENSION(klon), INTENT(OUT)        :: tslab
    REAL, DIMENSION(klon), INTENT(OUT)        :: seaice ! glace de mer (kg/m2)
    REAL, DIMENSION(klon, nbsrf), INTENT(OUT) :: pctsrf_slab

! Local variables
!****************************************************************************************
    REAL                    :: amn, amx
    REAL, PARAMETER         :: unjour=86400.
    REAL, PARAMETER         :: cyang=50.0 * 4.228e+06 ! capacite calorifique volumetrique de l'eau J/(m2 K)
    REAL, PARAMETER         :: cbing=0.334e+05        ! J/kg
    REAL, DIMENSION(klon)   :: siceh !hauteur de la glace de mer (m)
    INTEGER                 :: i
    REAL                    :: zz, za, zb
!
!****************************************************************************************
!
    julien = MOD(ijour,360)

    IF (check ) THEN
       amn=MIN(tmp_tslab_loc(1),1000.)
       amx=MAX(tmp_tslab_loc(1),-1000.)
       DO i=2, klon
          amn=MIN(tmp_tslab_loc(i),amn)
          amx=MAX(tmp_tslab_loc(i),amx)
       ENDDO

       WRITE(lunout,*) ' debut min max tslab',amn,amx 
       WRITE(lunout,*) ' itap,lmt_pas unjour',itap,lmt_pas,unjour
    ENDIF

    pctsrf_slab(1:klon,1:nbsrf) = pctsrf(1:klon,1:nbsrf)
!
! lecture du bilan au sol lmt_bils issu d'une simulation forcee en debut de journee
! 
    IF (MOD(itap,lmt_pas) .EQ. 1) THEN 
       ! 1er pas de temps de la journee
       idayvrai = ijour
       CALL condsurf(julien,idayvrai, lmt_bils)
    ENDIF

    DO i = 1, klon
       IF((pctsrf_slab(i,is_oce).GT.epsfra).OR. &
            (pctsrf_slab(i,is_sic).GT.epsfra)) THEN
!
! fabriquer de la glace si congelation atteinte:
!          
          IF (tmp_tslab_loc(i).LT.(RTT-1.8)) THEN
             zz =  (RTT-1.8)-tmp_tslab_loc(i)
             tmp_seaice(i) = tmp_seaice(i) + cyang/cbing * zz
             seaice(i) = tmp_seaice(i)
             tmp_tslab_loc(i) = RTT-1.8
          ENDIF
!
! faire fondre de la glace si temperature est superieure a 0:
!
          IF ((tmp_tslab_loc(i).GT.RTT) .AND. (tmp_seaice(i).GT.0.0)) THEN
             zz = cyang/cbing * (tmp_tslab_loc(i)-RTT)
             zz = MIN(zz,tmp_seaice(i))
             tmp_seaice(i) = tmp_seaice(i) - zz
             seaice(i) = tmp_seaice(i)
             tmp_tslab_loc(i) = tmp_tslab_loc(i) - zz*cbing/cyang
          ENDIF
!
! limiter la glace de mer a 10 metres (10000 kg/m2)
!
          IF(tmp_seaice(i).GT.45.) THEN
             tmp_seaice(i) = MIN(tmp_seaice(i),10000.0)
          ELSE
             tmp_seaice(i) = 0. 
          ENDIF
          seaice(i) = tmp_seaice(i)
          siceh(i)=tmp_seaice(i)/1000. !en metres
!
! determiner la nature du sol (glace de mer ou ocean libre):
!
! on fait dependre la fraction de seaice "pctsrf(i,is_sic)" 
! de l'epaisseur de seaice :
! pctsrf(i,is_sic)=1. si l'epaisseur de la glace de mer est >= a 20cm
! 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) WRITE(lunout,*) 'cas ok_slab_sicOBS TRUE : passe sur la modif.'
          ENDIF !(.NOT.ok_slab_sicOBS) then
       ENDIF !pctsrf
    ENDDO
!
! Calculer le bilan du flux de chaleur au sol :
!
    DO i = 1, klon
       za = radsol(i) + fluxo(i)
       zb = fluxg(i)
       IF((pctsrf_slab(i,is_oce).GT.epsfra).OR. &
            (pctsrf_slab(i,is_sic).GT.epsfra)) THEN
          slab_bils(i)=slab_bils(i)+(za*pctsrf_slab(i,is_oce) &
               +zb*pctsrf_slab(i,is_sic))/ FLOAT(lmt_pas)
       ENDIF
    ENDDO !klon
!
! calcul tslab 
!
    IF (MOD(itap,lmt_pas).EQ.0) THEN !fin de journee
!
! calcul tslab 
       amn=MIN(tmp_tslab_loc(1),1000.)
       amx=MAX(tmp_tslab_loc(1),-1000.)
       DO i = 1, klon
          IF ((pctsrf_slab(i,is_oce).GT.epsfra).OR. &
               (pctsrf_slab(i,is_sic).GT.epsfra)) THEN
             tmp_tslab_loc(i) = tmp_tslab_loc(i) + &
                  (slab_bils(i)-lmt_bils(i)) &
                  /cyang*unjour

! on remet l'accumulation a 0
             slab_bils(i) = 0.
          ENDIF !pctsrf
!
          IF (check) THEN
             IF(i.EQ.1) THEN  
                WRITE(lunout,*) 'i tmp_tslab_loc MOD(itap,lmt_pas).EQ.0 sicINTER',i,itap, &
                     tmp_tslab_loc(i), tmp_seaice(i)
             ENDIF
             
             amn=MIN(tmp_tslab_loc(i),amn)
             amx=MAX(tmp_tslab_loc(i),amx)
          ENDIF
       ENDDO !klon
    ENDIF !(MOD(itap,lmt_pas).EQ.0) THEN

    IF ( check ) THEN
       WRITE(lunout,*) 'fin min max tslab',amn,amx 
    ENDIF

    tslab  = tmp_tslab_loc
    seaice = tmp_seaice

  END SUBROUTINE interfoce_slab
!
!****************************************************************************************  
!
  SUBROUTINE ocean_slab_final(tslab_rst, seaice_rst)

! This subroutine will send to phyredem the variables concerning the slab 
! ocean that should be written to restart file.

!****************************************************************************************

    REAL, DIMENSION(klon), INTENT(OUT) :: tslab_rst
    REAL, DIMENSION(klon), INTENT(OUT) :: seaice_rst

!****************************************************************************************
! Set the output variables
    tslab_rst(:)  = tmp_tslab(:)
!    tslab_rst(:)  = tmp_tslab_loc(:)
    seaice_rst(:) = tmp_seaice(:)

! Deallocation of all variables in module
    DEALLOCATE(tmp_tslab, tmp_tslab_loc, tmp_pctsrf_slab)
    DEALLOCATE(tmp_seaice, tmp_radsol, tmp_flux_o, tmp_flux_g)
    DEALLOCATE(slab_bils, lmt_bils)

  END SUBROUTINE ocean_slab_final
!
!****************************************************************************************
!
  SUBROUTINE ocean_slab_get_vars(tslab_loc, seaice_loc, flux_o_loc, flux_g_loc)
! "Get some variables from module ocean_slab_mod"
! This subroutine prints variables to a external routine

    REAL, DIMENSION(klon), INTENT(OUT) :: tslab_loc
    REAL, DIMENSION(klon), INTENT(OUT) :: seaice_loc
    REAL, DIMENSION(klon), INTENT(OUT) :: flux_o_loc
    REAL, DIMENSION(klon), INTENT(OUT) :: flux_g_loc

! Set the output variables
    tslab_loc(:)  = tmp_tslab(:)
    seaice_loc(:) = tmp_seaice(:)
    flux_o_loc(:) = tmp_flux_o(:)
    flux_g_loc(:) = tmp_flux_g(:)

  END SUBROUTINE ocean_slab_get_vars
!
!****************************************************************************************
!
END MODULE ocean_slab_mod