!
! $Header$
!
MODULE ocean_forced_mod
!
! This module is used for both the sub-surfaces ocean and sea-ice for the case of a 
! forced ocean,  "ocean=force".
!
  USE surface_data,     ONLY : calice, calsno, tau_gl
  USE fonte_neige_mod,  ONLY : fonte_neige
  USE calcul_fluxs_mod, ONLY : calcul_fluxs
  USE dimphy

  IMPLICIT NONE

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

CONTAINS
!
!****************************************************************************************
! 
  SUBROUTINE ocean_forced_init
! Allocate fields needed for this module
!    
    INTEGER              :: error
    CHARACTER (len = 80) :: abort_message
    CHARACTER (len = 20) :: modname = 'ocean_forced_init'
!****************************************************************************************

    ALLOCATE(tmp_flux_o(1: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(1:klon), stat = error)
    IF (error /= 0) THEN
       abort_message='Pb allocation tmp_flux_g'
       CALL abort_gcm(modname,abort_message,1)
    ENDIF    


  END SUBROUTINE ocean_forced_init
!
!****************************************************************************************
!
  SUBROUTINE ocean_forced_noice(itime, dtime, jour, knon, knindex, &
       debut, &
       p1lay, tq_cdrag, precip_rain, precip_snow, &
       temp_air, spechum, &
       petAcoef, peqAcoef, petBcoef, peqBcoef, &
       ps, u1_lay, v1_lay, &
       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 routine reads data from climatologie file and does some calculations at the 
! surface. 
!
    INCLUDE "indicesol.inc"
    INCLUDE "YOMCST.inc"

! Input arguments
!****************************************************************************************
    INTEGER, INTENT(IN)                      :: itime, jour, knon
    INTEGER, DIMENSION(klon), INTENT(IN)     :: knindex
    LOGICAL, INTENT(IN)                      :: debut
    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 !? put to 0 in ocean
  
! 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)       :: alb_neig, tsurf_lim, zx_sl
    LOGICAL                     :: check=.FALSE.
    REAL, DIMENSION(klon,nbsrf) :: pctsrf_lim

!****************************************************************************************
! Start calculation
!****************************************************************************************
    IF (check) WRITE(*,*)' Entering ocean_forced_noice'

!****************************************************************************************
! 1)    
! Read from climatologie file SST and fraction of sub-surfaces
!
!****************************************************************************************
! Get from file tsurf_lim and pctsrf_lim
    CALL interfoce_lim(itime, dtime, jour, & 
         knon, knindex, &
         debut, &
         tsurf_lim, pctsrf_lim)
    
!****************************************************************************************
! 2)
! Flux calculation
!
!****************************************************************************************
! Set some variables for calcul_fluxs
    cal = 0.
    beta = 1.
    dif_grnd = 0.
    alb_neig(:) = 0.
    agesno(:) = 0.
  
! Calcul de tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l and qsurf
    CALL calcul_fluxs(knon, is_oce, dtime, &
         tsurf_lim, 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)

!****************************************************************************************
! 3)
! Calculate tmp_flux_o
!
!****************************************************************************************
!IM: flux ocean-atmosphere utile pour le "slab" ocean
! The flux are written to hist file
    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_lim(knindex(i),is_oce) .GT. epsfra) THEN
          tmp_flux_o(knindex(i)) = fluxsens(i) + fluxlat(i)
       ENDIF
    ENDDO

!****************************************************************************************
! 4)
! Return the new values for the ocean fraction
!
!****************************************************************************************

    pctsrf_oce(:) = pctsrf_lim(:,is_oce)

  END SUBROUTINE ocean_forced_noice
!
!****************************************************************************************
!
  SUBROUTINE ocean_forced_ice(itime, dtime, jour, knon, knindex, &
       debut, &
       tsurf_in, 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, alblw, evap, fluxsens, fluxlat, &
       tsurf_new, alb_new, dflux_s, dflux_l, pctsrf_sic)
!
! This subroutine treats the ocean where there is ice. 
! The routine reads data from climatologie file and does flux calculations at the 
! surface.
!    
    INCLUDE "indicesol.inc"
    INCLUDE "dimsoil.h"
    INCLUDE "YOMCST.inc"
    INCLUDE "clesphys.inc"

! Input arguments
!****************************************************************************************
    INTEGER, INTENT(IN)                  :: itime, jour, knon
    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
    LOGICAL, INTENT(IN)                  :: debut
    REAL, INTENT(IN)                     :: dtime
    REAL, DIMENSION(klon), INTENT(IN)    :: tsurf_in
    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, 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)            :: 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
!****************************************************************************************
    LOGICAL                     :: check=.FALSE.
    INTEGER                     :: i
    REAL                        :: zfra
    REAL, PARAMETER             :: t_grnd=271.35
    REAL, DIMENSION(klon)       :: cal, beta, dif_grnd, capsol
    REAL, DIMENSION(klon)       :: alb_neig, tsurf_tmp
    REAL, DIMENSION(klon)       :: soilcap, soilflux
    REAL, DIMENSION(klon,nbsrf) :: pctsrf_lim

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

!****************************************************************************************
! 1)    
! Read from climatologie file SST and fraction of sub-surfaces
!
!****************************************************************************************
    CALL interfoce_lim(itime, dtime, jour, & 
         knon, knindex, &
         debut, &
         tsurf_tmp, pctsrf_lim)
  
    DO i = 1, knon
       tsurf_tmp(i) = tsurf_in(i)
       IF (pctsrf_lim(knindex(i),is_sic) < EPSFRA) THEN
          snow(i) = 0.0
          tsurf_tmp(i) = RTT - 1.8
          IF (soil_model) tsoil(i,:) = RTT -1.8
       ENDIF
    ENDDO

!****************************************************************************************
! 2) 
! Flux calculation : tsurf_new, evap, fluxlat, fluxsens, 
!                    dflux_s, dflux_l and qsurf
!****************************************************************************************

! calculate the parameters cal, beta, capsol and dif_grnd

    CALL calbeta(dtime, is_sic, knon, snow, qsol, beta, capsol, dif_grnd)
    
    IF (soil_model) THEN 
! update tsoil and calculate soilcap and soilflux
       CALL soil(dtime, is_sic, knon,snow, tsurf_tmp, tsoil,soilcap, soilflux)
       cal(1:knon) = RCPD / soilcap(1:knon)
       radsol(1:knon) = radsol(1:knon)  + soilflux(1:knon)
       dif_grnd = 1.0 / tau_gl
    ELSE 
       dif_grnd = 1.0 / tau_gl
       cal = RCPD * calice
       WHERE (snow > 0.0) cal = RCPD * calsno 
    ENDIF

    beta = 1.0
    CALL calcul_fluxs(knon, is_sic, dtime, &
         tsurf_tmp, 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)

!****************************************************************************************
! 3)
! Calculations due to snow and runoff
!
!****************************************************************************************
    CALL fonte_neige( knon, is_sic, knindex, dtime, &
         tsurf_tmp, precip_rain, precip_snow, &
         snow, qsol, tsurf_new, evap)
    
! Calculation of albedo at snow (alb_neig) and update the age of snow (agesno)
! 
    CALL albsno(klon, knon, dtime, agesno(:), alb_neig(:), precip_snow(:))  

    WHERE (snow(1:knon) .LT. 0.0001) agesno(1:knon) = 0.

    alb_new(:) = 0.0
    DO i=1, knon
       zfra = MAX(0.0,MIN(1.0,snow(i)/(snow(i)+10.0)))
       alb_new(i) = alb_neig(i) * zfra +  0.6 * (1.0-zfra)
    ENDDO
!!      alb_new(1 : knon) = 0.6

!****************************************************************************************
! 4)
! Calculate tmp_flux_g
!
!****************************************************************************************
!
    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_lim(knindex(i),is_sic) .GT. epsfra) &
            tmp_flux_g(knindex(i)) = (tsurf_new(i)-t_grnd) * dif_grnd(i) *RCPD/cal(i)
       
    ENDDO

!!$    z0_new = 0.002
!!$    z0_new = SQRT(z0_new**2+rugoro**2)
    alblw(1:knon) = alb_new(1:knon)

!****************************************************************************************
! 5)
! Return the new values for the seaice fraction
!
!****************************************************************************************

    pctsrf_sic(:) = pctsrf_lim(:,is_sic)

  END SUBROUTINE ocean_forced_ice
!
!****************************************************************************************
!
  SUBROUTINE ocean_forced_get_vars(flux_o, flux_g)
! Get some variables from module oceanforced.
! This subroutine returns variables to a external routine

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

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

  END SUBROUTINE ocean_forced_get_vars
!
!****************************************************************************************
!
END MODULE ocean_forced_mod