coupling_seafluxn.F90

Go to the documentation of this file.

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
!     ###############################################################################
SUBROUTINE coupling_seaflux_n (CHS, DTS, DGS, O, OR, G, S, DST, SLT, &
                               HPROGRAM, HCOUPLING, PTIMEC, PTSTEP, KYEAR, KMONTH, KDAY, PTIME, &
                               KI, KSV, KSW, PTSUN, PZENITH, PZENITH2, PAZIM, PZREF, PUREF,     &
                               PU, PV, PQA, PTA, PRHOA, PSV, PCO2, HSV, PRAIN, PSNOW, PLW,      &
                               PDIR_SW, PSCA_SW, PSW_BANDS, PPS, PPA, PSFTQ, PSFTH, PSFTS,      &
                               PSFCO2, PSFU, PSFV, PTRAD, PDIR_ALB, PSCA_ALB, PEMIS, PTSURF,    &
                               PZ0, PZ0H, PQSURF, PPEW_A_COEF, PPEW_B_COEF, PPET_A_COEF,        &
                               PPEQ_A_COEF, PPET_B_COEF, PPEQ_B_COEF,  HTEST                 )  
!     ###############################################################################
!
!!****  *COUPLING_SEAFLUX_n * - Driver of the WATER_FLUX scheme for sea   
!!
!!    PURPOSE
!!    -------
!
!!**  METHOD
!!    ------
!!
!!    REFERENCE
!!    ---------
!!      
!!
!!    AUTHOR
!!    ------
!!     V. Masson 
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    01/2004
!!      Modified    01/2006 : sea flux parameterization.
!!      Modified    09/2006 : P. Tulet Introduce Sea salt aerosol Emission/Deposition
!!      Modified    03/2009 : B. Decharme SST could change during a run => ALB and EMIS 
!!      Modified    05/2009 : V. Masson : implicitation of momentum fluxes
!!      Modified    09/2009 : B. Decharme Radiative properties at time t+1
!!      Modified    01/2010 : B. Decharme Add XTTS
!!      Modified    09/2012 : B. Decharme New wind implicitation
!!      Modified    10/2012 : P. Le Moigne CMO1D update
!!      Modified    04/2013 : P. Le Moigne Wind implicitation and SST update displaced
!!      Modified    04/2013 : B. Decharme new coupling variables
!!      Modified    01/2014 : S. Senesi : handle sea-ice cover, sea-ice model interface, 
!!                               and apply to Gelato
!!      Modified    01/2014 : S. Belamari Remove MODE_THERMOS and XLVTT
!!      Modified    05/2014 : S. Belamari New ECUME : Include salinity & atm. pressure impact 
!!      Modified    01/2015 : R. Séférian interactive ocaen surface albedo
!!                                       
!!---------------------------------------------------------------------
!
!
USE modd_ch_seaflux_n, ONLY : ch_seaflux_t
USE modd_data_seaflux_n, ONLY : data_seaflux_t
USE modd_surfex_n, ONLY : seaflux_diag_t
USE modd_ocean_n, ONLY : ocean_t
USE modd_ocean_rel_n, ONLY : ocean_rel_t
USE modd_sfx_grid_n, ONLY : grid_t
USE modd_seaflux_n, ONLY : seaflux_t
!
USE modd_dst_n, ONLY : dst_t
USE modd_slt_n, ONLY : slt_t
!
USE modd_reprod_oper, ONLY : cimplicit_wind
!
USE modd_csts,       ONLY : xrd, xcpd, xp00, xtt, xtts, xttsi, xday
USE modd_surf_par,   ONLY : xundef
USE modd_sfx_oasis,  ONLY : lcpl_sea, lcpl_seaice
USE modd_water_par,  ONLY : xemiswat, xemiswatice
!
USE modd_water_par, ONLY : xalbseaice
!
!
USE modi_water_flux
USE modi_mr98
USE modi_ecume_seaflux
USE modi_coare30_seaflux
USE modi_add_forecast_to_date_surf
USE modi_mod1d_n
USE modi_diag_inline_seaflux_n
USE modi_ch_aer_dep
USE modi_ch_dep_water
USE modi_dslt_dep
USE modi_sst_update
USE modi_interpol_sst_mth
USE modi_update_rad_sea
!
USE mode_dslt_surf
USE modd_dst_surf
USE modd_slt_surf
! 
USE modd_ocean_grid,   ONLY : nockmin
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
USE modi_abor1_sfx
!
USE modi_coupling_iceflux_n
USE modi_seaice_gelato1d_n
!
USE modi_coupling_slt_n
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
!
TYPE(ch_seaflux_t), INTENT(INOUT) :: CHS
TYPE(data_seaflux_t), INTENT(INOUT) :: DTS
TYPE(seaflux_diag_t), INTENT(INOUT) :: DGS
TYPE(ocean_t), INTENT(INOUT) :: O
TYPE(ocean_rel_t), INTENT(INOUT) :: OR
TYPE(grid_t), INTENT(INOUT) :: G
TYPE(seaflux_t), INTENT(INOUT) :: S 
TYPE(dst_t), INTENT(INOUT) :: DST
TYPE(slt_t), INTENT(INOUT) :: SLT
!
 CHARACTER(LEN=6),    INTENT(IN)  :: HPROGRAM  ! program calling surf. schemes
 CHARACTER(LEN=1),    INTENT(IN)  :: HCOUPLING ! type of coupling
                                      ! 'E' : explicit
                                      ! 'I' : implicit
REAL,                INTENT(IN)  :: PTIMEC    ! current duration since start of the run (s)
INTEGER,             INTENT(IN)  :: KYEAR     ! current year (UTC)
INTEGER,             INTENT(IN)  :: KMONTH    ! current month (UTC)
INTEGER,             INTENT(IN)  :: KDAY      ! current day (UTC)
REAL,                INTENT(IN)  :: PTIME     ! current time since midnight (UTC, s)
INTEGER,             INTENT(IN)  :: KI        ! number of points
INTEGER,             INTENT(IN)  :: KSV       ! number of scalars
INTEGER,             INTENT(IN)  :: KSW       ! number of short-wave spectral bands
REAL, DIMENSION(KI), INTENT(IN)  :: PTSUN     ! solar time                    (s from midnight)
REAL,                INTENT(IN)  :: PTSTEP    ! atmospheric time-step                 (s)
REAL, DIMENSION(KI), INTENT(IN)  :: PZREF     ! height of T,q forcing                 (m)
REAL, DIMENSION(KI), INTENT(IN)  :: PUREF     ! height of wind forcing                (m)
!
REAL, DIMENSION(KI), INTENT(IN)  :: PTA       ! air temperature forcing               (K)
REAL, DIMENSION(KI), INTENT(IN)  :: PQA       ! air humidity forcing                  (kg/m3)
REAL, DIMENSION(KI), INTENT(IN)  :: PRHOA     ! air density                           (kg/m3)
REAL, DIMENSION(KI,KSV),INTENT(IN) :: PSV     ! scalar variables
!                                             ! chemistry:   first char. in HSV: '#'  (molecule/m3)
!                                             !
 CHARACTER(LEN=6), DIMENSION(KSV),INTENT(IN):: HSV  ! name of all scalar variables
REAL, DIMENSION(KI), INTENT(IN)  :: PU        ! zonal wind                            (m/s)
REAL, DIMENSION(KI), INTENT(IN)  :: PV        ! meridian wind                         (m/s)
REAL, DIMENSION(KI,KSW),INTENT(IN) :: PDIR_SW ! direct  solar radiation (on horizontal surf.)
!                                             !                                       (W/m2)
REAL, DIMENSION(KI,KSW),INTENT(IN) :: PSCA_SW ! diffuse solar radiation (on horizontal surf.)
!                                             !                                       (W/m2)
REAL, DIMENSION(KSW),INTENT(IN)  :: PSW_BANDS ! mean wavelength of each shortwave band (m)
REAL, DIMENSION(KI), INTENT(IN)  :: PZENITH   ! zenithal angle at t  (radian from the vertical)
REAL, DIMENSION(KI), INTENT(IN)  :: PZENITH2  ! zenithal angle at t+1(radian from the vertical)
REAL, DIMENSION(KI), INTENT(IN)  :: PAZIM     ! azimuthal angle      (radian from North, clockwise)
REAL, DIMENSION(KI), INTENT(IN)  :: PLW       ! longwave radiation (on horizontal surf.)
!                                             !                                       (W/m2)
REAL, DIMENSION(KI), INTENT(IN)  :: PPS       ! pressure at atmospheric model surface (Pa)
REAL, DIMENSION(KI), INTENT(IN)  :: PPA       ! pressure at forcing level             (Pa)
REAL, DIMENSION(KI), INTENT(IN)  :: PCO2      ! CO2 concentration in the air          (kg/m3)
REAL, DIMENSION(KI), INTENT(IN)  :: PSNOW     ! snow precipitation                    (kg/m2/s)
REAL, DIMENSION(KI), INTENT(IN)  :: PRAIN     ! liquid precipitation                  (kg/m2/s)
!
REAL, DIMENSION(KI), INTENT(OUT) :: PSFTH     ! flux of heat                          (W/m2)
REAL, DIMENSION(KI), INTENT(OUT) :: PSFTQ     ! flux of water vapor                   (kg/m2/s)
REAL, DIMENSION(KI), INTENT(OUT) :: PSFU      ! zonal momentum flux                   (Pa)
REAL, DIMENSION(KI), INTENT(OUT) :: PSFV      ! meridian momentum flux                (Pa)
REAL, DIMENSION(KI), INTENT(OUT) :: PSFCO2    ! flux of CO2                           (m/s*kg_CO2/kg_air)
REAL, DIMENSION(KI,KSV),INTENT(OUT):: PSFTS   ! flux of scalar var.                   (kg/m2/s)
!
REAL, DIMENSION(KI), INTENT(OUT) :: PTRAD     ! radiative temperature                 (K)
REAL, DIMENSION(KI,KSW),INTENT(OUT):: PDIR_ALB! direct albedo for each spectral band  (-)
REAL, DIMENSION(KI,KSW),INTENT(OUT):: PSCA_ALB! diffuse albedo for each spectral band (-)
REAL, DIMENSION(KI), INTENT(OUT) :: PEMIS     ! emissivity                            (-)
!
REAL, DIMENSION(KI), INTENT(OUT) :: PTSURF    ! surface effective temperature         (K)
REAL, DIMENSION(KI), INTENT(OUT) :: PZ0       ! roughness length for momentum         (m)
REAL, DIMENSION(KI), INTENT(OUT) :: PZ0H      ! roughness length for heat             (m)
REAL, DIMENSION(KI), INTENT(OUT) :: PQSURF    ! specific humidity at surface          (kg/kg)
!
REAL, DIMENSION(KI), INTENT(IN) :: PPEW_A_COEF! implicit coefficients   (m2s/kg)
REAL, DIMENSION(KI), INTENT(IN) :: PPEW_B_COEF! needed if HCOUPLING='I' (m/s)
REAL, DIMENSION(KI), INTENT(IN) :: PPET_A_COEF
REAL, DIMENSION(KI), INTENT(IN) :: PPEQ_A_COEF
REAL, DIMENSION(KI), INTENT(IN) :: PPET_B_COEF
REAL, DIMENSION(KI), INTENT(IN) :: PPEQ_B_COEF
 CHARACTER(LEN=2),    INTENT(IN) :: HTEST ! must be equal to 'OK'
!
!*      0.2    declarations of local variables
!     
REAL, DIMENSION(KI,KSW) :: ZDIR_ALB   ! Direct albedo at time t
REAL, DIMENSION(KI,KSW) :: ZSCA_ALB   ! Diffuse albedo at time t
!
REAL, DIMENSION(KI) :: ZEXNA      ! Exner function at forcing level
REAL, DIMENSION(KI) :: ZEXNS      ! Exner function at surface level
REAL, DIMENSION(KI) :: ZU         ! zonal wind
REAL, DIMENSION(KI) :: ZV         ! meridian wind
REAL, DIMENSION(KI) :: ZWIND      ! Wind
REAL, DIMENSION(KI) :: ZCD        ! Drag coefficient on open sea
REAL, DIMENSION(KI) :: ZCD_ICE    ! "     "          on seaice
REAL, DIMENSION(KI) :: ZCDN       ! Neutral Drag coefficient on open sea
REAL, DIMENSION(KI) :: ZCDN_ICE   ! "     "          on seaice
REAL, DIMENSION(KI) :: ZCH        ! Heat transfer coefficient on open sea
REAL, DIMENSION(KI) :: ZCH_ICE    ! "     "          on seaice
REAL, DIMENSION(KI) :: ZCE        ! Vaporization heat transfer coefficient on open sea
REAL, DIMENSION(KI) :: ZCE_ICE    ! "     "          on seaice
REAL, DIMENSION(KI) :: ZRI        ! Richardson number on open sea
REAL, DIMENSION(KI) :: ZRI_ICE    ! "     "          on seaice
REAL, DIMENSION(KI) :: ZRESA_SEA  ! aerodynamical resistance on open sea
REAL, DIMENSION(KI) :: ZRESA_SEA_ICE  ! "     "          on seaice
REAL, DIMENSION(KI) :: ZUSTAR     ! friction velocity (m/s) on open sea
REAL, DIMENSION(KI) :: ZUSTAR_ICE ! "     "          on seaice
REAL, DIMENSION(KI) :: ZZ0        ! roughness length over open sea
REAL, DIMENSION(KI) :: ZZ0_ICE    ! roughness length over seaice
REAL, DIMENSION(KI) :: ZZ0H       ! heat roughness length over open sea
REAL, DIMENSION(KI) :: ZZ0H_ICE   ! heat roughness length over seaice
REAL, DIMENSION(KI) :: ZZ0W       ! Work array for Z0 and Z0H computation
REAL, DIMENSION(KI) :: ZQSAT      ! humidity at saturation on open sea
REAL, DIMENSION(KI) :: ZQSAT_ICE  ! "     "          on seaice
!
REAL, DIMENSION(KI) :: ZSFTH      ! Heat flux for open sea (and for sea-ice points if merged)
REAL, DIMENSION(KI) :: ZSFTQ      ! Water vapor flux on open sea (and for sea-ice points if merged)
REAL, DIMENSION(KI) :: ZSFU       ! zonal      momentum flux on open sea (and for sea-ice points if merged)(Pa)
REAL, DIMENSION(KI) :: ZSFV       ! meridional momentum flux on open sea (and for sea-ice points if merged)(Pa)
!
REAL, DIMENSION(KI) :: ZSFTH_ICE  ! Heat flux on sea ice 
REAL, DIMENSION(KI) :: ZSFTQ_ICE  ! Sea-ice sublimation flux
REAL, DIMENSION(KI) :: ZSFU_ICE   ! zonal      momentum flux on seaice (Pa)
REAL, DIMENSION(KI) :: ZSFV_ICE   ! meridional momentum flux on seaice (Pa)

REAL, DIMENSION(KI) :: ZHU        ! Near surface relative humidity
REAL, DIMENSION(KI) :: ZQA        ! specific humidity (kg/kg)
REAL, DIMENSION(KI) :: ZEMIS      ! Emissivity at time t
REAL, DIMENSION(KI) :: ZTRAD      ! Radiative temperature at time t
!
REAL, DIMENSION(KI) :: ZSST       ! XSST corrected for anomalously low values (which actually are sea-ice temp)
REAL, DIMENSION(KI) :: ZMASK      ! A mask for diagnosing where seaice exists (or, for coupling_iceflux, may appear)
!
REAL                             :: ZCONVERTFACM0_SLT, ZCONVERTFACM0_DST
REAL                             :: ZCONVERTFACM3_SLT, ZCONVERTFACM3_DST
REAL                             :: ZCONVERTFACM6_SLT, ZCONVERTFACM6_DST
!
INTEGER                          :: ISIZE_WATER  ! number of points with some sea water 
INTEGER                          :: ISIZE_ICE    ! number of points with some sea ice
!
INTEGER                          :: ISWB       ! number of shortwave spectral bands
INTEGER                          :: JSWB       ! loop counter on shortwave spectral bands
INTEGER                          :: ISLT       ! number of sea salt variable
!
INTEGER :: IBEG, IEND
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!-------------------------------------------------------------------------------------
! Preliminaries:
!-------------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('COUPLING_SEAFLUX_N',0,zhook_handle)
IF (htest/='OK') THEN
  CALL abor1_sfx('COUPLING_SEAFLUXN: FATAL ERROR DURING ARGUMENT TRANSFER')        
END IF
!-------------------------------------------------------------------------------------
!
zexna(:) = xundef
zexns(:) = xundef
zu(:) = xundef
zv(:) = xundef
zwind(:) = xundef
zsftq(:) = xundef
zsfth(:) = xundef
zcd(:) = xundef    
zcdn(:) = xundef
zch(:) = xundef
zce(:) = xundef
zri(:) = xundef
zhu(:) = xundef
zresa_sea(:) = xundef
zustar(:) = xundef
zz0(:) = xundef
zz0h(:) = xundef
zqsat(:) = xundef
!
zsftq_ice(:) = xundef
zsfth_ice(:) = xundef
zcd_ice(:) = xundef    
zcdn_ice(:) = xundef
zch_ice(:) = xundef
zce_ice(:) = xundef
zri_ice(:) = xundef
zresa_sea_ice= xundef
zustar_ice(:) = xundef
zz0_ice(:) = xundef
zz0h_ice(:) = xundef
zqsat_ice(:) = xundef
!
zemis(:) = xundef
ztrad(:) = xundef
zdir_alb(:,:) = xundef
zsca_alb(:,:) = xundef
!
!-------------------------------------------------------------------------------------
!
zexns(:)     = (pps(:)/xp00)**(xrd/xcpd)
zexna(:)     = (ppa(:)/xp00)**(xrd/xcpd)
!
IF(lcpl_sea)THEN 
  !Sea currents are taken into account
  zu(:)=pu(:)-s%XUMER(:)
  zv(:)=pv(:)-s%XVMER(:)
ELSE
  zu(:)=pu(:)
  zv(:)=pv(:)        
ENDIF
!
zwind(:) = sqrt(zu(:)**2+zv(:)**2)
!
psfts(:,:) = 0.
!
zhu = 1.
!
zqa(:) = pqa(:) / prhoa(:)
!
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Time evolution
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
s%TTIME%TIME = s%TTIME%TIME + ptstep
 CALL add_forecast_to_date_surf(s%TTIME%TDATE%YEAR,s%TTIME%TDATE%MONTH,s%TTIME%TDATE%DAY,s%TTIME%TIME)
!
!--------------------------------------------------------------------------------------
! Fluxes over water according to Charnock formulae
!--------------------------------------------------------------------------------------
!
IF (s%LHANDLE_SIC) THEN 
  ! Flux for sea are computed everywhere
  isize_water = SIZE(zmask)
  ! Ensure freezing SST values where XSST actually has very low (sea-ice) values (old habits)
  zsst(:)=max(s%XSST(:), xttsi) 
  ! Flux over sea-ice will not be computed by next calls, but by coupling_iceflux. Hence :
  isize_ice   = 0
  ! Flux over sea-ice will be computed by coupling_iceflux anywhere sea-ice could form in one 
  ! time-step (incl. under forcing). ZMASK value is set to 1. on these points
  zmask(:)=0.
  WHERE ( s%XSIC(:) > 0. ) zmask(:)=1. 
  ! To be large, assume that seaice may form where SST is < 10C
  WHERE ( s%XSST(:) - xtts <= 10. ) zmask(:)=1.
  IF (s%LINTERPOL_SIC) WHERE (s%XFSIC(:) > 0. ) zmask(:)=1. 
  IF (s%LINTERPOL_SIT) WHERE (s%XFSIT(:) > 0. ) zmask(:)=1.
ELSE
  zsst(:) = s%XSST(:)
  zmask(:) = s%XSST(:) - xtts
  isize_water = count(zmask(:)>=0.)
  isize_ice   = SIZE(s%XSST) - isize_water
ENDIF
!
SELECT CASE (s%CSEA_FLUX)
CASE ('DIRECT')
 CALL water_flux(s%XZ0, pta, zexna, prhoa, zsst, zexns, zqa, &
                prain, psnow, xtts, zwind, pzref, puref,    &
                pps, s%LHANDLE_SIC, zqsat,  zsfth, zsftq,   &
                zustar, zcd, zcdn, zch, zri, zresa_sea, zz0h )  
CASE ('ITERAT')
 CALL mr98     (s%XZ0, pta, zexna, prhoa, s%XSST, zexns, zqa,  &
               xtts, zwind, pzref, puref, pps, zqsat,         &
               zsfth, zsftq, zustar,                          &
               zcd, zcdn, zch, zri, zresa_sea, zz0h           )  

CASE ('ECUME ','ECUME6')
 CALL ecume_seaflux(s, zmask, isize_water, isize_ice,       &
              pta, zexna ,prhoa, zsst, zexns, zqa,         &
              prain, psnow, zwind, pzref, puref, pps, ppa, &
              zqsat, zsfth, zsftq, zustar,                 &
              zcd, zcdn, zch, zce, zri, zresa_sea, zz0h    )
CASE ('COARE3')
 CALL coare30_seaflux(s, zmask, isize_water, isize_ice,        &
              pta, zexna ,prhoa, zsst, zexns, zqa, prain,     &
              psnow, zwind, pzref, puref, pps, zqsat,         &
              zsfth, zsftq, zustar,                           &
              zcd, zcdn, zch, zce, zri, zresa_sea, zz0h       )  
END SELECT
!
!-------------------------------------------------------------------------------------
!radiative properties at time t
!-------------------------------------------------------------------------------------
!
iswb = SIZE(psw_bands)
!
DO jswb=1,iswb
zdir_alb(:,jswb) = s%XDIR_ALB(:)
zsca_alb(:,jswb) = s%XSCA_ALB(:)
END DO
!
IF (s%LHANDLE_SIC) THEN 
zemis(:) =   (1 - s%XSIC(:)) * xemiswat    + s%XSIC(:) * xemiswatice
ztrad(:) = (((1 - s%XSIC(:)) * xemiswat    * s%XSST (:)**4 + &
             s%XSIC(:)  * xemiswatice * s%XTICE(:)**4)/ zemis(:)) ** 0.25
ELSE
ztrad(:) = s%XSST (:)
zemis(:) = s%XEMIS(:)
END IF
!
!-------------------------------------------------------------------------------------
!Specific fields for seaice model (when using earth system model or embedded  
!seaice scheme)
!-------------------------------------------------------------------------------------
!
IF(lcpl_seaice.OR.s%LHANDLE_SIC)THEN
 CALL coupling_iceflux_n(ki, pta, zexna, prhoa, s%XTICE, zexns, &
             zqa, prain, psnow, zwind, pzref, puref,     &
             pps, s%XSST, xtts, zsfth_ice, zsftq_ice,    &  
             s%LHANDLE_SIC, zmask, zqsat_ice, zz0_ice,   &
             zustar_ice, zcd_ice, zcdn_ice, zch_ice,   &
             zri_ice, zresa_sea_ice, zz0h_ice          )
ENDIF
!
IF (s%LHANDLE_SIC) CALL complement_each_other_flux
!
!-------------------------------------------------------------------------------------
! Momentum fluxes over sea or se-ice
!-------------------------------------------------------------------------------------
!
 CALL sea_momentum_fluxes(zcd, zsfu, zsfv)
!
! Momentum fluxes over sea-ice if embedded seaice scheme is used
!
IF (s%LHANDLE_SIC) CALL sea_momentum_fluxes(zcd_ice, zsfu_ice, zsfv_ice)
!
! CO2 flux
!
psfco2(:) = 0.0
!
!IF(LCPL_SEA.AND.CSEACO2=='NONE')THEN
!  PSFCO2(:) = XSEACO2(:)
!ELSEIF(CSEACO2=='CST ')THEN
!  PSFCO2 = E * deltapCO2 
!  According to Wanninkhof (medium hypothesis) : 
!  E = 1.13.10^-3 * WIND^2 CO2mol.m-2.yr-1.uatm-1
!    = 1.13.10^-3 * WIND^2 * Mco2.10^-3 * (1/365*24*3600)
!  deltapCO2 = -8.7 uatm (Table 1 half hypothesis)
psfco2(:) = - zwind(:)**2 * 1.13e-3 * 8.7 * 44.e-3 / ( 365*24*3600 )
!ENDIF
!
!-------------------------------------------------------------------------------------
! Scalar fluxes:
!-------------------------------------------------------------------------------------
!
IF (chs%SVS%NBEQ>0) THEN
  !
  IF (chs%CCH_DRY_DEP == "WES89") THEN
    !
    ibeg = chs%SVS%NSV_CHSBEG
    iend = chs%SVS%NSV_CHSEND
    !
    CALL ch_dep_water  (zresa_sea, zustar, pta, ztrad,psv(:,ibeg:iend), &
                        chs%SVS%CSV(ibeg:iend), chs%XDEP(:,1:chs%SVS%NBEQ) )  
    !
    psfts(:,ibeg:iend) = - psv(:,ibeg:iend) * chs%XDEP(:,1:chs%SVS%NBEQ)  
    !
    IF (chs%SVS%NAEREQ > 0 ) THEN
      !
      ibeg = chs%SVS%NSV_AERBEG
      iend = chs%SVS%NSV_AEREND
      !
      CALL ch_aer_dep(psv(:,ibeg:iend),psfts(:,ibeg:iend),zustar,zresa_sea,pta,prhoa)   
      !  
    END IF
    !
  ELSE
    !
    ibeg = chs%SVS%NSV_AERBEG
    iend = chs%SVS%NSV_AEREND
    !
    psfts(:,ibeg:iend) =0.
    IF (iend.GT.ibeg)     psfts(:,ibeg:iend) =0.
    !
  ENDIF
  !
ENDIF
!
IF (chs%SVS%NDSTEQ>0) THEN
  !
  ibeg = chs%SVS%NSV_DSTBEG
  iend = chs%SVS%NSV_DSTEND
  !
  CALL dslt_dep(psv(:,ibeg:iend), psfts(:,ibeg:iend), zustar, zresa_sea, pta, &
                prhoa, dst%XEMISSIG_DST, dst%XEMISRADIUS_DST, jpmode_dst,     &
                xdensity_dst, xmolarweight_dst, zconvertfacm0_dst, zconvertfacm6_dst, &
                zconvertfacm3_dst, lvarsig_dst, lrgfix_dst, cvermod  )  
  !
  CALL massflux2momentflux(         &
                psfts(:,ibeg:iend), & !I/O ![kg/m2/sec] In: flux of only mass, out: flux of moments
                prhoa,                          & !I [kg/m3] air density
                dst%XEMISRADIUS_DST,                &!I [um] emitted radius for the modes (max 3)
                dst%XEMISSIG_DST,                   &!I [-] emitted sigma for the different modes (max 3)
                ndstmde,                        &
                zconvertfacm0_dst,              &
                zconvertfacm6_dst,              &
                zconvertfacm3_dst,              &
                lvarsig_dst, lrgfix_dst         )  
  !
ENDIF

!
IF (chs%SVS%NSLTEQ>0) THEN
  !
  ibeg = chs%SVS%NSV_SLTBEG
  iend = chs%SVS%NSV_SLTEND
  !
  islt = iend - ibeg + 1
  !
  CALL coupling_slt_n(slt, &
                      SIZE(zustar,1),           & !I [nbr] number of sea point
                      islt,                     & !I [nbr] number of sea salt variables
                      zwind,                    & !I [m/s] wind velocity
                      psfts(:,ibeg:iend) )   
  !
  CALL dslt_dep(psv(:,ibeg:iend), psfts(:,ibeg:iend), zustar, zresa_sea, pta, &
                prhoa, slt%XEMISSIG_SLT, slt%XEMISRADIUS_SLT, jpmode_slt,     &
                xdensity_slt, xmolarweight_slt, zconvertfacm0_slt, zconvertfacm6_slt, &
                zconvertfacm3_slt, lvarsig_slt, lrgfix_slt, cvermod  )  
  !
  CALL massflux2momentflux(         &
                psfts(:,ibeg:iend), & !I/O ![kg/m2/sec] In: flux of only mass, out: flux of moments
                prhoa,                          & !I [kg/m3] air density
                slt%XEMISRADIUS_SLT,            &!I [um] emitted radius for the modes (max 3)
                slt%XEMISSIG_SLT,               &!I [-] emitted sigma for the different modes (max 3)
                nsltmde,                        &
                zconvertfacm0_slt,              &
                zconvertfacm6_slt,              &
                zconvertfacm3_slt,              &
                lvarsig_slt, lrgfix_slt         ) 
  !
ENDIF
!
!-------------------------------------------------------------------------------
! Inline diagnostics at time t for SST and TRAD
!-------------------------------------------------------------------------------
!
 CALL diag_inline_seaflux_n(dgs%O, dgs%D, dgs%DC, dgs%DI, dgs%DIC, dgs%DMI, &
                           s, ptstep, pta, zqa, ppa, pps, prhoa, pu,       &
                           pv, pzref, puref, zcd, zcdn, zch, zce, zri, zhu,&
                           zz0h, zqsat, zsfth, zsftq, zsfu, zsfv,          &
                           pdir_sw, psca_sw, plw, zdir_alb, zsca_alb,      &
                           zemis, ztrad, prain, psnow,                     &
                           zcd_ice, zcdn_ice, zch_ice, zce_ice, zri_ice,   &
                           zz0_ice, zz0h_ice, zqsat_ice, zsfth_ice,        &
                           zsftq_ice, zsfu_ice, zsfv_ice)
!
!-------------------------------------------------------------------------------
! A kind of "average_flux"
!-------------------------------------------------------------------------------
!
IF (s%LHANDLE_SIC) THEN
  psfth(:) = zsfth(:) * ( 1 - s%XSIC (:)) + zsfth_ice(:) * s%XSIC(:)
  psftq(:) = zsftq(:) * ( 1 - s%XSIC (:)) + zsftq_ice(:) * s%XSIC(:)
  psfu(:) = zsfu(:) * ( 1 - s%XSIC (:)) +  zsfu_ice(:) * s%XSIC(:)
  psfv(:) = zsfv(:) * ( 1 - s%XSIC (:)) +  zsfv_ice(:) * s%XSIC(:)
ELSE
  psfth(:) = zsfth(:) 
  psftq(:) = zsftq(:) 
  psfu(:) = zsfu(:) 
  psfv(:) = zsfv(:) 
ENDIF
!
!-------------------------------------------------------------------------------
! IMPOSED SSS OR INTERPOLATED SSS AT TIME t+1
!-------------------------------------------------------------------------------
!
! Daily update Sea surface salinity from monthly data
!
IF (s%LINTERPOL_SSS .AND. mod(s%TTIME%TIME,xday) == 0.) THEN
  CALL interpol_sst_mth(s,'S')
  IF (any(s%XSSS(:)<0.0)) THEN
    CALL abor1_sfx('COUPLING_SEAFLUX_N: XSSS should be >=0') 
  ENDIF                      
ENDIF
!
!-------------------------------------------------------------------------------
! SEA-ICE coupling at time t+1
!-------------------------------------------------------------------------------
!
IF (s%LHANDLE_SIC) THEN
  !
  IF (s%LINTERPOL_SIC) THEN
    IF ((mod(s%TTIME%TIME,xday) == 0.) .OR. (ptimec <= ptstep )) THEN
      ! Daily update Sea Ice Cover constraint from monthly data
      CALL interpol_sst_mth(s,'C')
      IF (any(s%XFSIC(:)>1.0).OR.any(s%XFSIC(:)<0.0)) THEN
        CALL abor1_sfx('COUPLING_SEAFLUX_N: FSIC should be >=0 and <=1') 
      ENDIF                    
    ENDIF
  ENDIF
  !
  IF (s%LINTERPOL_SIT) THEN
    IF ((mod(s%TTIME%TIME,xday) == 0.) .OR. (ptimec <= ptstep )) THEN
      ! Daily update Sea Ice Thickness constraint from monthly data
      CALL interpol_sst_mth(s,'H')
      IF (any(s%XFSIT(:)<0.0)) THEN
        CALL abor1_sfx('COUPLING_SEAFLUX_N: XFSIT should be >=0') 
      ENDIF  
    ENDIF
  ENDIF
  !
  IF (s%CSEAICE_SCHEME=='GELATO') THEN
    CALL seaice_gelato1d_n(s, hprogram,ptimec, ptstep)
  ENDIF
  ! Update of cell-averaged albedo, emissivity and radiative 
  ! temperature is done later
ENDIF
!
!-------------------------------------------------------------------------------
! OCEANIC COUPLING, IMPOSED SST OR INTERPOLATED SST AT TIME t+1
!-------------------------------------------------------------------------------
!
IF (o%LMERCATOR) THEN
  !
  ! Update SST reference profile for relaxation purpose
  IF (dts%LSST_DATA) THEN
    CALL sst_update(dts, s, or%XSEAT_REL(:,nockmin+1))
    !
    ! Convert to degree C for ocean model
    or%XSEAT_REL(:,nockmin+1) = or%XSEAT_REL(:,nockmin+1) - xtt
  ENDIF
  !
  CALL mod1d_n(dgs%GO, o, or, g%XLAT, s, &
               hprogram,ptime,zemis(:),zdir_alb(:,1:ksw),zsca_alb(:,1:ksw),&
               plw(:),psca_sw(:,1:ksw),pdir_sw(:,1:ksw),psfth(:),          &
               psftq(:),psfu(:),psfv(:),prain(:))
  !
ELSEIF(dts%LSST_DATA) THEN 
  !
  ! Imposed SST 
  !
  CALL sst_update(dts, s, s%XSST)
  !
ELSEIF (s%LINTERPOL_SST.AND.mod(s%TTIME%TIME,xday) == 0.) THEN
   !
   ! Imposed monthly SST 
   !
   CALL interpol_sst_mth(s,'T')
   !
ENDIF
!
!-------------------------------------------------------------------------------
!Physical properties see by the atmosphere in order to close the energy budget 
!between surfex and the atmosphere. All variables should be at t+1 but very 
!difficult to do. Maybe it will be done later. However, Ts is at time t+1
!-------------------------------------------------------------------------------
!
IF (s%LHANDLE_SIC) THEN
  IF (s%CSEAICE_SCHEME/='GELATO') THEN
     s%XTICE    = s%XSST
     s%XSIC     = s%XFSIC
     s%XICE_ALB = xalbseaice           
  ENDIF         
  ptsurf(:) = s%XSST(:) * ( 1 - s%XSIC (:)) +   s%XTICE(:) * s%XSIC(:)
  pqsurf(:) = zqsat(:) * ( 1 - s%XSIC (:)) + zqsat_ice(:) * s%XSIC(:)
  zz0w(:) = ( 1 - s%XSIC(:) ) * 1.0/(log(puref(:)/zz0(:))    **2)  +  &
                     s%XSIC(:)   * 1.0/(log(puref(:)/zz0_ice(:))**2)  
  pz0(:) = puref(:) * exp( - sqrt( 1./  zz0w(:) ))
  zz0w(:) = ( 1 - s%XSIC(:) ) * 1.0/(log(pzref(:)/zz0h(:))    **2)  +  &
                     s%XSIC(:)   * 1.0/(log(pzref(:)/zz0h_ice(:))**2)  
  pz0h(:) = pzref(:) * exp( - sqrt( 1./  zz0w(:) ))
ELSE
  ptsurf(:) = s%XSST(:) 
  pqsurf(:) = zqsat(:) 
  pz0(:) = s%XZ0 (:) 
  pz0h(:) = zz0h(:) 
ENDIF
!
!-------------------------------------------------------------------------------
!Radiative properties at time t+1 (see by the atmosphere) in order to close
!the energy budget between surfex and the atmosphere
!-------------------------------------------------------------------------------
!
 CALL update_rad_sea(s,pzenith2,xtts,pdir_alb,psca_alb,pemis,ptrad,pu,pv)
!
!=======================================================================================
!
IF (lhook) CALL dr_hook('COUPLING_SEAFLUX_N',1,zhook_handle)
!
!=======================================================================================
!
CONTAINS
!
SUBROUTINE sea_momentum_fluxes(PCD, PSFU, PSFV)
!
IMPLICIT NONE
!
REAL, DIMENSION(KI), INTENT(IN)  :: PCD       ! Drag coefficient (on open sea or seaice)
REAL, DIMENSION(KI), INTENT(OUT) :: PSFU      ! zonal momentum flux                   (Pa)
REAL, DIMENSION(KI), INTENT(OUT) :: PSFV      ! meridian momentum flux                (Pa)
!
REAL, DIMENSION(KI) :: ZUSTAR2    ! square of friction velocity (m2/s2)
REAL, DIMENSION(KI) :: ZWORK      ! Work array
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
IF (lhook) CALL dr_hook('COUPLING_SEAFLUX_N: SEA_MOMENTUM_FLUXES',0,zhook_handle)
!
zwork(:) = xundef
zustar2(:) = xundef
!
IF(cimplicit_wind=='OLD')THEN
! old implicitation (m2/s2)
  zustar2(:) = (pcd(:)*zwind(:)*ppew_b_coef(:)) /            &
              (1.0-prhoa(:)*pcd(:)*zwind(:)*ppew_a_coef(:))
ELSE
! new implicitation (m2/s2)
  zustar2(:) = (pcd(:)*zwind(:)*(2.*ppew_b_coef(:)-zwind(:))) /&
              (1.0-2.0*prhoa(:)*pcd(:)*zwind(:)*ppew_a_coef(:))
!                   
  zwork(:)  = prhoa(:)*ppew_a_coef(:)*zustar2(:) + ppew_b_coef(:)
  zwork(:) = max(zwork(:),0.)
!
  WHERE(ppew_a_coef(:)/= 0.)
        zustar2(:) = max( ( zwork(:) - ppew_b_coef(:) ) / (prhoa(:)*ppew_a_coef(:)), 0.)
  ENDWHERE
!              
ENDIF
!
psfu = 0.
psfv = 0.
WHERE (zwind(:)>0.)
  psfu(:) = - prhoa(:) * zustar2(:) * zu(:) / zwind(:)
  psfv(:) = - prhoa(:) * zustar2(:) * zv(:) / zwind(:)
END WHERE
!
IF (lhook) CALL dr_hook('COUPLING_SEAFLUX_N: SEA_MOMENTUM_FLUXES',1,zhook_handle)
!
END SUBROUTINE sea_momentum_fluxes
!
!=======================================================================================
!
SUBROUTINE complement_each_other_flux
!
! Provide dummy fluxes on places with no open-sea or no sea-ice 
! Allows a smooth computing of CLS parameters in all cases while avoiding 
! having to pack arrays (in routines PARAM_CLS and CLS_TQ)
!
IMPLICIT NONE
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
IF (lhook) CALL dr_hook('COUPLING_SEAFLUX_N: COMPLEMENT_EACH_OTHER_FLUX',0,zhook_handle)
!
  WHERE (s%XSIC(:) == 1.)
     zsfth=zsfth_ice 
     zsftq=zsftq_ice 
     zsfu=zsfu_ice
     zsfv=zsfv_ice
     zqsat=zqsat_ice
     zcd=zcd_ice
     zcdn=zcdn_ice
     zch=zch_ice
     zce=zce_ice
     zri=zri_ice 
     zz0h=zz0h_ice
  END WHERE
  WHERE (s%XSIC(:) == 0.)
     zsfth_ice=zsfth 
     zsftq_ice=zsftq 
     zsfu_ice=zsfu
     zsfv_ice=zsfv
     zqsat_ice=zqsat
     zcd_ice=zcd
     zcdn_ice=zcdn
     zch_ice=zch
     zce_ice=zce
     zri_ice=zri 
     zz0h_ice=zz0h
  END WHERE
!
IF (lhook) CALL dr_hook('COUPLING_SEAFLUX_N: COMPLEMENT_EACH_OTHER_FLUX',1,zhook_handle)
!
END SUBROUTINE complement_each_other_flux
!
!=======================================================================================
!
END SUBROUTINE coupling_seaflux_n