!
MODULE ocean_slab_mod
!
! This module is used for both surface ocean and sea-ice when using the slab ocean,
! "ocean=slab".
!

  USE dimphy
  USE indice_sol_mod

  IMPLICIT NONE
  PRIVATE
  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
!
!****************************************************************************************
!
  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

!    INCLUDE "clesphys.h"

! Arguments
!****************************************************************************************
    INTEGER, INTENT(IN)                        :: itime   ! numero du pas de temps courant
    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_chg  ! sub-surface fraction
    LOGICAL, INTENT(OUT)                       :: is_modified ! true if pctsrf is modified at this time step

! Local variables
!****************************************************************************************

    IF (version_ocean == 'sicOBS'.OR. version_ocean == 'sicNO') THEN   
       CALL limit_read_frac(itime, dtime, jour, pctsrf_chg, is_modified)
    ELSE
       pctsrf_chg(:,:)=pctsrf(:,:)
       is_modified=.TRUE.
    END IF

  END SUBROUTINE ocean_slab_frac
!
!****************************************************************************************
!
  SUBROUTINE ocean_slab_noice( & 
       itime, dtime, jour, knon, knindex, &
       p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
       AcoefH, AcoefQ, BcoefH, BcoefQ, &
       AcoefU, AcoefV, BcoefU, BcoefV, &
       ps, u1, v1, tsurf_in, &
       radsol, snow, agesno, &
       qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
       tsurf_new, dflux_s, dflux_l, qflux)
    
    USE calcul_fluxs_mod
    USE surface_data

    INCLUDE "iniprint.h"

! Input arguments
!****************************************************************************************
    INTEGER, INTENT(IN)                  :: itime
    INTEGER, INTENT(IN)                  :: jour
    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)    :: cdragh, cdragm
    REAL, DIMENSION(klon), INTENT(IN)    :: precip_rain, precip_snow
    REAL, DIMENSION(klon), INTENT(IN)    :: temp_air, spechum
    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefH, AcoefQ, BcoefH, BcoefQ
    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefU, AcoefV, BcoefU, BcoefV
    REAL, DIMENSION(klon), INTENT(IN)    :: ps
    REAL, DIMENSION(klon), INTENT(IN)    :: u1, v1
    REAL, DIMENSION(klon), INTENT(IN)    :: tsurf_in

! 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)   :: flux_u1, flux_v1
    REAL, DIMENSION(klon), INTENT(OUT)   :: tsurf_new
    REAL, DIMENSION(klon), INTENT(OUT)   :: dflux_s, dflux_l      
    REAL, DIMENSION(klon), INTENT(OUT)   :: qflux

! Local variables
!****************************************************************************************
    INTEGER               :: i,ki
    REAL, DIMENSION(klon) :: cal, beta, dif_grnd
    REAL, DIMENSION(klon) :: diff_sst, lmt_bils
    REAL, DIMENSION(klon) :: u0, v0
    REAL, DIMENSION(klon) :: u1_lay, v1_lay

!****************************************************************************************
! 1) Flux calculation
!
!****************************************************************************************
    cal(:)      = 0.
    beta(:)     = 1.
    dif_grnd(:) = 0.
    agesno(:)   = 0.
    
! Suppose zero surface speed
    u0(:)=0.0
    v0(:)=0.0
    u1_lay(:) = u1(:) - u0(:)
    v1_lay(:) = v1(:) - v0(:)

    CALL calcul_fluxs(knon, is_oce, dtime, &
         tsurf_in, p1lay, cal, beta, cdragh, ps, &
         precip_rain, precip_snow, snow, qsurf,  &
         radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, &
         AcoefH, AcoefQ, BcoefH, BcoefQ, &
         tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)

! - Flux calculation at first modele level for U and V
    CALL calcul_flux_wind(knon, dtime, &
         u0, v0, u1, v1, cdragm, &
         AcoefU, AcoefV, BcoefU, BcoefV, &
         p1lay, temp_air, &
         flux_u1, flux_v1)  

! 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
!
!***********************************************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