!
! $Header$
!
MODULE ocean_cpl_mod
!
! This module is used both for the sub-surface ocean and sea-ice for the case of a 
! coupled model configuration, ocean=couple. 
!

  USE dimphy,           ONLY : klon
  USE cpl_mod
  USE calcul_fluxs_mod, ONLY : calcul_fluxs
  USE climb_wind_mod,   ONLY : calcul_wind_flux

  IMPLICIT NONE
  PRIVATE

  PUBLIC :: ocean_cpl_init, ocean_cpl_get_vars, ocean_cpl_noice, ocean_cpl_ice

  REAL, ALLOCATABLE, DIMENSION(:), SAVE       :: tmp_flux_o
  !$OMP THREADPRIVATE(tmp_flux_o)
  REAL, ALLOCATABLE, DIMENSION(:), SAVE       :: tmp_flux_g
  !$OMP THREADPRIVATE(tmp_flux_g)

!****************************************************************************************
!
CONTAINS
!
!****************************************************************************************
!
  SUBROUTINE ocean_cpl_init(dtime, rlon, rlat)
!
! Allocate fields for this module and initailize the module mod_cpl
!
! Input arguments
!*************************************************************************************
    REAL, INTENT(IN)                  :: dtime
    REAL, DIMENSION(klon), INTENT(IN) :: rlon, rlat

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


    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

! Initialize module cpl_init
    CALL cpl_init(dtime, rlon, rlat)
    
  END SUBROUTINE ocean_cpl_init
!
!****************************************************************************************
!
  SUBROUTINE ocean_cpl_noice( &
       sollw, albedo, &
       windsp, &
       fder_old, &
       itime, dtime, knon, knindex, &
       swdown, &
       p1lay, tq_cdrag, precip_rain, precip_snow, temp_air, spechum, &
       petAcoef, peqAcoef, petBcoef, peqBcoef, &
       ps, u1_lay, v1_lay, pctsrf_in, &
       radsol, snow, qsurf, agesno, &
       evap, fluxsens, fluxlat, &
       tsurf_new, dflux_s, dflux_l, pctsrf_oce)
!
! This subroutine treats the "open ocean", all grid points that are not entierly covered
! by ice. The subroutine first receives fields from coupler, then some calculations at 
! surface is done and finally it sends some fields to the coupler.
!
    INCLUDE "indicesol.h"
    INCLUDE "YOMCST.h"
!    
! Input arguments  
!****************************************************************************************
    INTEGER, INTENT(IN)                      :: itime, knon
    INTEGER, DIMENSION(klon), INTENT(IN)     :: knindex
    REAL, INTENT(IN)                         :: dtime
    REAL, DIMENSION(klon), INTENT(IN)        :: sollw
    REAL, DIMENSION(klon), INTENT(IN)        :: albedo
    REAL, DIMENSION(klon), INTENT(IN)        :: windsp
    REAL, DIMENSION(klon), INTENT(IN)        :: fder_old
    REAL, DIMENSION(klon), INTENT(IN)        :: swdown    
    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
    REAL, DIMENSION(klon,nbsrf), INTENT(IN)  :: pctsrf_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)       :: tsurf_new
    REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l      
    REAL, DIMENSION(klon), INTENT(OUT)       :: pctsrf_oce

! Local variables
!****************************************************************************************
    INTEGER               :: i
    INTEGER, DIMENSION(1) :: iloc
    REAL, DIMENSION(klon) :: cal, beta, dif_grnd
    REAL, DIMENSION(klon) :: zx_sl
    REAL, DIMENSION(klon) :: fder_new
    REAL, DIMENSION(klon) :: tsurf_cpl
    REAL, DIMENSION(klon) :: taux, tauy
    LOGICAL               :: check=.FALSE.

! End definitions
!****************************************************************************************

    IF (check) WRITE(*,*)' Entering ocean_cpl_noice'

!****************************************************************************************
! Receive sea-surface temperature(tsurf_cpl) and new fraction of ocean surface(pctsrf_oce) 
! from coupler
!
!****************************************************************************************
    CALL cpl_receive_ocean_fields(itime, dtime, knon, knindex, pctsrf_in, &
         tsurf_cpl, pctsrf_oce)

!****************************************************************************************
! Calculate fluxes at surface
!
!****************************************************************************************
    cal = 0.
    beta = 1.
    dif_grnd = 0.
    agesno(:) = 0.
    
    CALL calcul_fluxs(knon, is_oce, dtime, &
         tsurf_cpl, 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)
    
    ! Calcultate the flux of u and v at surface
    CALL calcul_wind_flux(knon, dtime, taux, tauy)
    

!****************************************************************************************
! Calculate fder : flux derivative (sensible and latente)
!
!****************************************************************************************
    fder_new(:) = fder_old(:) + dflux_s(:) + dflux_l(:)
    
    iloc = MAXLOC(fder_new(1:klon))
    IF (check .AND. fder_new(iloc(1))> 0.) THEN
       WRITE(*,*)'**** Debug fder****'
       WRITE(*,*)'max fder(',iloc(1),') = ',fder_new(iloc(1))
       WRITE(*,*)'fder_old, dflux_s, dflux_l',fder_old(iloc(1)), &
            dflux_s(iloc(1)), dflux_l(iloc(1))
    ENDIF

!****************************************************************************************
! Flux ocean-atmosphere useful for "slab" ocean but here calculated only for printing
! usage later in physiq  
! 
!****************************************************************************************
    tmp_flux_o(:) = 0.0
    DO i=1, knon
       zx_sl(i) = RLVTT
       IF (tsurf_new(i) .LT. RTT) zx_sl(i) = RLSTT
       !IM     flux_o(i) = fluxsens(i)-evap(i)*zx_sl(i)
       !       flux_o(i) = fluxsens(i) + fluxlat(i)
       IF (pctsrf_oce(knindex(i)) .GT. epsfra) THEN
          tmp_flux_o(knindex(i)) = fluxsens(i) + fluxlat(i)
       ENDIF
    ENDDO


!****************************************************************************************
! Send and cumulate fields to the coupler
!
!****************************************************************************************

    CALL cpl_send_ocean_fields(itime, knon, knindex, &
         swdown, sollw, fluxlat, fluxsens, &
         precip_rain, precip_snow, evap, tsurf_new, fder_new, albedo, taux, tauy, windsp)
    

  END SUBROUTINE ocean_cpl_noice
!
!****************************************************************************************
!
  SUBROUTINE ocean_cpl_ice( &
       rlon, rlat, sollw, albedo, &
       fder_old, &
       itime, dtime, knon, knindex, &
       lafin, &
       swdown, &
       p1lay, tq_cdrag, precip_rain, precip_snow, temp_air, spechum, &
       petAcoef, peqAcoef, petBcoef, peqBcoef, &
       ps, u1_lay, v1_lay, pctsrf_in, &
       radsol, snow, qsurf, &
       alblw, evap, fluxsens, fluxlat, &
       tsurf_new, alb_new, dflux_s, dflux_l, pctsrf_sic)
!
! This subroutine treats the ocean where there is ice. The subroutine first receives 
! fields from coupler, then some calculations at surface is done and finally sends 
! some fields to the coupler.
!    
    INCLUDE "indicesol.h"
    INCLUDE "YOMCST.h"

! Input arguments
!****************************************************************************************
    INTEGER, INTENT(IN)                      :: itime, knon
    INTEGER, DIMENSION(klon), INTENT(IN)     :: knindex
    LOGICAL, INTENT(IN)                      :: lafin
    REAL, INTENT(IN)                         :: dtime
    REAL, DIMENSION(klon), INTENT(IN)        :: rlon, rlat
    REAL, DIMENSION(klon), INTENT(IN)        :: sollw
    REAL, DIMENSION(klon), INTENT(IN)        :: albedo
    REAL, DIMENSION(klon), INTENT(IN)        :: fder_old
    REAL, DIMENSION(klon), INTENT(IN)        :: swdown
    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
    REAL, DIMENSION(klon,nbsrf), INTENT(IN)  :: pctsrf_in

! In/output arguments
!****************************************************************************************
    REAL, DIMENSION(klon), INTENT(INOUT)     :: radsol
    REAL, DIMENSION(klon), INTENT(INOUT)     :: snow

! Output arguments
!****************************************************************************************
    REAL, DIMENSION(klon), INTENT(OUT)       :: qsurf
    REAL, DIMENSION(klon), INTENT(OUT)       :: alblw
    REAL, DIMENSION(klon), INTENT(OUT)       :: evap, fluxsens, fluxlat
    REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new, alb_new
    REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l      
    REAL, DIMENSION(klon), INTENT(OUT)       :: pctsrf_sic

! Local variables
!****************************************************************************************
    INTEGER                 :: i
    INTEGER, DIMENSION(1)   :: iloc
    LOGICAL                 :: check=.FALSE.
    REAL, PARAMETER         :: t_grnd=271.35
    REAL, DIMENSION(klon)   :: cal, beta, dif_grnd
    REAL, DIMENSION(klon)   :: tsurf_cpl, fder_new
    REAL, DIMENSION(klon)   :: taux, tauy

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

!****************************************************************************************
! Receive ocean temperature(tsurf_cpl), albedo(alb_new) and new fraction of 
! seaice(pctsrf_sic) from coupler
!
!****************************************************************************************

    CALL cpl_receive_seaice_fields(knon, knindex, &
         tsurf_cpl, alb_new, pctsrf_sic)
    
!****************************************************************************************
! Calculate fluxes at surface
!
!****************************************************************************************
    cal = 0.
    dif_grnd = 0.
    beta = 1.0
    

    CALL calcul_fluxs(knon, is_sic, dtime, &
         tsurf_cpl, 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)

    ! Calcultate the flux of u and v at surface
    CALL calcul_wind_flux(knon, dtime, taux, tauy)
    
!****************************************************************************************
! Flux ocean-atmosphere useful for "slab" ocean but here calculated only for printing
! usage later in physiq  
! 
!  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
!
!****************************************************************************************
    tmp_flux_g(:) = 0.0
    DO i = 1, knon
       IF (cal(i) .GT. 1.0e-15 .AND. pctsrf_sic(knindex(i)) .GT. epsfra) &
            tmp_flux_g(knindex(i)) = (tsurf_new(i)-t_grnd) * &
            dif_grnd(i) * RCPD/cal(i) 
    ENDDO
    
!****************************************************************************************
! Calculate fder : flux derivative (sensible and latente)
!
!****************************************************************************************
    fder_new(:) = fder_old(:) + dflux_s(:) + dflux_l(:)
    
    iloc = MAXLOC(fder_new(1:klon))
    IF (check .AND. fder_new(iloc(1))> 0.) THEN
       WRITE(*,*)'**** Debug fder ****'
       WRITE(*,*)'max fder(',iloc(1),') = ',fder_new(iloc(1))
       WRITE(*,*)'fder_old, dflux_s, dflux_l',fder_old(iloc(1)), &
            dflux_s(iloc(1)), dflux_l(iloc(1))
    ENDIF

!****************************************************************************************
! Send and cumulate fields to the coupler
!
!****************************************************************************************

    CALL cpl_send_seaice_fields(itime, dtime, knon, knindex, &
       pctsrf_in, lafin, rlon, rlat, &
       swdown, sollw, fluxlat, fluxsens, &
       precip_rain, precip_snow, evap, tsurf_new, fder_new, albedo, taux, tauy)
 

    alblw(1:knon) = alb_new(1:knon)    

  END SUBROUTINE ocean_cpl_ice
!  
!****************************************************************************************
!
  SUBROUTINE ocean_cpl_get_vars(flux_o, flux_g)

! This subroutine returns variables private in this module to an external 
! routine (physiq).

    REAL, DIMENSION(klon), INTENT(OUT) :: flux_o
    REAL, DIMENSION(klon), INTENT(OUT) :: flux_g

! Set the output variables
    flux_o(:) = tmp_flux_o(:)
    flux_g(:) = tmp_flux_g(:)

  END SUBROUTINE ocean_cpl_get_vars
!
!****************************************************************************************
!
END MODULE ocean_cpl_mod
