 MODULE sflxturb
!!======================================================================
!!                   ***  MODULE  sflxturb  ***
!! Surface fluxes : estimate of turbulent (latent and sensible) atmospherical fluxes
!!
!!      The upward sensible heat flux depends on (1) the wind speed,
!!      (2) the difference between the sea temperature and the air
!!      temperature and (3) the heat capacity rho_a c_p per unit volume
!!      of the air.
!!      Q_s = rho_a c_p c_H u (T_s - T_a), c_h being a dimensionless
!!      coefficient, sometimes called the Stanton number.
!!
!!      The upward latent heat flux depends on (1) the wind speed,
!!      (2) the difference between the specific humidities (mass of
!!      water vapor per unit mass of air) of the sea and the air
!!      and (3) the air density and the latent heat of vaporization
!!      of water.
!!      Q_l = rho_a Lnu c_E u (q_s - q_a), c_E being a dimensionless
!!      coefficient, sometimes called the Dalton number.
!!
!!======================================================================





























































!! * Modules used
   USE parameters,ONLY : ljmin,ljmax,ljminp1,imin,imax,valmanq,rsh
   USE cometdyn,  ONLY : cds
   USE cometdyn,  ONLY : l_sflx_rrain,rhoair
   USE cometemp,  ONLY : l_rapflux

   IMPLICIT NONE

   PRIVATE

!! * Accessibility
   PUBLIC sflx_coef_init        ! subroutine called by main.F90, Agrif_User.F90
   PUBLIC sflx_turb             ! subroutine called by sflxflux.F90
   PUBLIC cd                    ! subroutine called by sflxcalc.F90 or sflxturb.F90

!! * Interface
   INTERFACE sflx_turb

     MODULE PROCEDURE sflx_turb_default

   END INTERFACE

!! * Shared module variables
   REAL(KIND=rsh), PUBLIC, PARAMETER   :: lnu=2.5e+6           ! latent heat of vaporization

!! * Private variables
   REAL(KIND=rsh), PARAMETER           :: em_w=0.985_rsh,   &  ! water emissivity
                                          suswe=0.98_rsh,   &  ! salty / unsalty sea water effect (-2% for salty water)
                                          kelvin=273.15_rsh,&  ! shift to use kelvin unit
                                          cp=1004.0_rsh,    &  ! specific heat of air at constant pressure (J/kg/K)
                                          zt=2.0_rsh           ! air temperature measurement height  [m]
! ATTENTION si temperature de l air a 10 m et non 2m idem pour humidite relative

   INTEGER                             ::   i,j                !  Loop indices

 CONTAINS

!!======================================================================

  SUBROUTINE sflx_turb_large(ztemp,zsat,zwindmod,zpa,zrh,zrad,zirlos,zrain, &
                             zsenslos,zlatlos,zevap,ztosx,ztosy)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE sflx_turb_large  ***
!&E
!&E ** Purpose : latent and sensible heat fluxes estimates
!&E
!&E ** Description :
!&E
!&E ** Called by : sflxsurf.F90
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables :
!&E
!&E ** Reference : Large and Yeager (2004)
!&E                METEOCEAN report (2009)
!&E
!&E ** History :
!&E             issued from OPA 9.0
!&E       !  2004-08  (U. Schweckendiek)  Original code
!&E       !  2005-04  (L. Brodeau, A.M. Treguier) severals modifications:
!&E                   - new bulk routine for efficiency
!&E                   - WINDS ARE NOW ASSUMED TO BE AT T POINTS in input files
!&E       !  2006-02  (S. Masson, G. Madec)  IOM read + NEMO "style"
!&E       !  2006-12  (L. Brodeau) handle both 2m and 10m input fields
!&E       !  2009-06  (B. Levier, V. Garnier) MARS implementation
!&E
!&E Rq vl : Symphonie utilise 6 iterations et remonte la temp potentielle de l eau en surface
!&E---------------------------------------------------------------------
!! * Modules used
   USE comvars2d,  ONLY : jb,jh,ig,id,dt,nbouc,iscreenlog
   USE cometdyn,   ONLY : l_sflx_rath2m
   USE toolmpi,  ONLY : findmax_mpi_3,MASTER

!! * Arguments
   REAL(KIND=rsh), INTENT( IN )      :: ztemp,zsat,zwindmod,zpa,zrh, &
                                        zrad,zirlos,zrain
   REAL(KIND=rsh), INTENT( OUT )     :: zsenslos,zlatlos,zevap
   REAL(KIND=rsh), INTENT( INOUT )   :: ztosx,ztosy

!! * Local declarations
!   REAL(KIND=rsh)            ::   &
!         zzu , zzv              ! scalar
   REAL(KIND=rsh)                                   ::   &
         zsst,            &  ! SST                                    [K]
         zat,             &  ! air temperature                        [K]
         zssu,            &  ! zonal sea surface current              [m/s]
         zssv,            &  ! meridional sea surface current         [m/s]
         dUnormt,         &  ! relative wind (norm) on T points       [m/s]
!         ztaux, ztauy,    &  ! wind stress computed at T-point
         sh_w,            &  ! specific humidity at zSST              [kg/kg]
         sh_a,            &  ! air spec. hum.                         [kg/kg]
         cheat,           &  ! coefficient for sensible heat flux
         cevap,           &  ! coefficient for latent heat flux
         cdrag,           &  ! coefficient for wind stress
         zat_zu,          &  ! air temperature at wind speed height   [K]
         sh_a_zu             ! air spec. hum.  at wind speed height   [kg/kg]
   REAL(KIND=rsh), PARAMETER ::   &
         rhoair  = 1.22_rsh,  &  ! air density                        [kg/m3]
         cp    = 1000.5_rsh,  &  ! specific heat of air               [J/kg/K]
         zu    = 10.0_rsh        ! height for dU                      [m]
   INTEGER :: ii,jj

!!----------------------------------------------------------------------
!! * Executable part

! --------------
! Initialization
! --------------

! Temperatures in Kelvin
! ----------------------

! sea temperature [Kelvin]
! Caution set to kelvin over land, not 0 Kelvin (otherwise floating point exception in bulk computation)
   zsst = MAX( (ztemp + kelvin) , 0.5_rsh )

! air potential temperature [Kelvin]
! the air temperature increases by 0.98 degrees Celsius per 100 meters
   zat = zsat + kelvin + 0.0098_rsh*zt


! Sea surface current to estimate scalar wind ( = | U10m - SSU | )
! It is important to take into account the sea surface courant
! ----------------------------------------------------------------

! Interpolation of surface current at T-point, zssu and zssv :
!   DO j = 2, ljmaxm1
!      DO i = 1, imax
!         zssu(i,j) = ( gu(i-1,j  ) + gu(i,j) ) * 0.5_rsh
!         zssv(i,j) = ( gv(i  ,j-1) + gv(i,j) ) * 0.5_rsh
!      END DO
!   END DO

   zssu = 0.0_rsh
   zssv = 0.0_rsh



! Module of relative wind
! -----------------------

! Wind components are provided on T-points within the netcdf input file.
! Thus, the wind module is computded at T-points taking into account the sea

!      zzu = MIN(ABS(u10 - zssu),50.0_rsh)
!      zzv = MIN(ABS(v10 - zssv),50.0_rsh)
!      dUnormt = SQRT( zzu*zzu + zzv*zzv )
       dUnormt = zwindmod

! specific humidities
! -------------------
   sh_w = suswe*640380.0_rsh*exp(-5107.4_rsh/zsst)/rhoair
   sh_a = zrh*640380.0_rsh*exp(-5107.4_rsh/zat)/rhoair


! ---------------------------------------------------------------------
! CORE iterative algo for computation of cdrag, cheat, cevap at T-point
! ---------------------------------------------------------------------

   IF ( l_sflx_rath2m ) THEN   ! air temp. and spec. hum. are given at different height (2m) than wind (10m)

     IF ( nbouc == 1 )   THEN
       IF (MASTER) THEN
       WRITE(iscreenlog,*) 
       WRITE(iscreenlog,*) 'Calling coef_large_2z for bulk transfert coefficients'
       WRITE(iscreenlog,*) 
       ENDIF
     ENDIF
     CALL coef_large_2z(zt,zu,zsst, zat, sh_w, sh_a, &
         &              dUnormt, cdrag, cheat, cevap, zat_zu, sh_a_zu)

   ELSE                ! air temp. and spec. hum. are given at same height than wind (10m)

     IF( nbouc == 1 )    THEN
       IF (MASTER) THEN
       WRITE(iscreenlog,*) 
       WRITE(iscreenlog,*) 'Calling coef_large_1z for bulk transfert coefficients'
       ENDIF
     ENDIF
     CALL coef_large_1z(zu, zsst, zat, sh_w, sh_a,  &
         &              dUnormt, cdrag, cheat, cevap )
   ENDIF


! -------------------------
! Momentum over ocean
! -------------------------

   ztosx=rhoair*cdrag*ztosx
   ztosy=rhoair*cdrag*ztosy

! dUnormt is the modified module of the wind
! from which all the exchange coefficients have been estimated.
! ztosx = ||V||*Vx ; ztosy = ||V||*Vy, so dUnormt/zwindmod removes
! the module of the wind in the formula,
! the second factor dUnormt/zwindmod corrects the component of the wind

   IF (zwindmod > 0.0_rsh) THEN
     ztosx=ztosx*(dUnormt/zwindmod)**2
     ztosy=ztosy*(dUnormt/zwindmod)**2
   ENDIF




! ---------------------------
! Turbulent fluxes over ocean
! ---------------------------

   IF( l_sflx_rath2m ) THEN

! Values of temp. and hum. adjusted to 10m must be used instead of 2m values
      zsenslos = -rhoair*cp*cheat*(zsst - zat_zu)*dUnormt
      zevap    = rhoair*cevap*(sh_w - sh_a_zu)*dUnormt

   ELSE

      zsenslos = -rhoair*cp*cheat*(zsst - zat)*dUnormt
      zevap    = rhoair*cevap*(sh_w - sh_a)*dUnormt

   END IF

   zlatlos =  -lnu * zevap

  END SUBROUTINE sflx_turb_large

!!======================================================================

  SUBROUTINE coef_large_1z(zu, sst, T_a, q_w, q_a,   &
   &                     dU, cdrag, cheat, cevap   )
!!----------------------------------------------------------------------
!!                      ***  ROUTINE  turb_core  ***
!!
!! ** Purpose :   Computes turbulent transfert coefficients of surface
!!                fluxes according to Large & Yeager (2004)
!!
!! ** Method  :   I N E R T I A L   D I S S I P A T I O N   M E T H O D
!!      Momentum, Latent and sensible heat exchange coefficients
!!      Caution: this procedure should only be used in cases when air
!!      temperature (T_air), air specific humidity (q_air) and wind (dU)
!!      are provided at the same height 'zzu' (10 m)!
!!
!! References :
!!   Large & Yeager, 2004 : Diurnal to decadal global forcing
!!                          for ocean and sea-ice models: the data sets
!!                          and flux climatologies. CGD Division of the National Center
!!                          for Atmospheric Research, NCAR Technical Note: NCAR/TN-460+STR.
!!
!! History :
!!        !  XX-XX  (???     )  Original code
!!   9.0  !  05-08  (L. Brodeau) Rewriting and optimization
!!----------------------------------------------------------------------
!! * Modules used
   USE parameters, ONLY : grav

!! * Arguments
   REAL(KIND=rsh), INTENT(in) :: zu           ! altitude of wind measurement       [m]
   REAL(KIND=rsh), INTENT(in)                                   ::  &
         sst,      &       ! sea surface temperature         [Kelvin]
         T_a,      &       ! potential air temperature       [Kelvin]
         q_w,      &       ! sea surface specific humidity   [kg/kg]
         q_a,      &       ! specific air humidity           [kg/kg]
         dU                ! wind module |U(zu)-U(0)|        [m/s]
   REAL(KIND=rsh), intent(out)                                    :: &
         cdrag,    &       ! transfert coefficient for momentum       (tau)
         cheat,    &       ! transfert coefficient for temperature (Q_sens)
         cevap             ! transfert coefficient for evaporation  (Q_lat)

!! * Local declarations
   REAL(KIND=rsh)                                    ::   &
         dU10,           & ! relative wind dU                     [m/s]
         dT,             & ! air/sea temperature differeence      [K]
         dq,             & ! air/sea humidity difference          [kg/kg]
         cdrag_n10,      & ! 10m neutral drag coefficient
         cevap_n10,      & ! 10m neutral latent coefficient
         cheat_n10,      & ! 10m neutral sensible coefficient
         sqrt_cdrag_n10, & ! root square of cdrag_n10
         sqrt_cdrag,     & ! root square of cdrag
         T_vpot,         & ! virtual potential temperature        [K]
         T_star,         & ! turbulent scale of tem. fluct.
         q_star,         & ! turbulent humidity of temp. fluct.
         U_star,         & ! turb. scale of velocity fluct.
         L,              & ! Monin-Obukov length                  [m]
         zeta,           & ! stability parameter at height zu
         U_n10,          & ! neutral wind velocity at 10m         [m]
         xlogt, xct, zpsi_h, zpsi_m
   REAL(KIND=rsh), PARAMETER        :: milli = 0.001_rsh
!!
   INTEGER :: j_itt
   INTEGER, PARAMETER :: nb_itt = 3
   INTEGER                                    ::   &
         stab              ! 1st guess stability test integer

   REAL(KIND=rsh), PARAMETER ::                        &
         karman  = 0.4_rsh          ! von Karman s constant
 
!! * Start
!! Air/sea differences
   dU10 = max(0.5_rsh, dU) ! we don t want to fall under 0.5 m/s
   dT = T_a - sst          ! assuming that T_a is allready the potential temp. at zzu
   dq = q_a - q_w
!!
!! Virtual potential temperature
   T_vpot = T_a*(1.0_rsh + 0.608_rsh*q_a)
!!
!! Neutral Drag Coefficient
   stab    = 0.5_rsh + sign(0.5_rsh,dT)    ! stable : stab = 1 ; unstable : stab = 0
   cdrag_n10  = milli * ( 2.7_rsh/dU10 + 0.142_rsh + dU10/13.09_rsh )    !   L & Y eq. (6a)
   sqrt_cdrag_n10 = sqrt(cdrag_n10)
   cevap_n10  = milli * ( 34.6_rsh * sqrt_cdrag_n10 )                     !   L & Y eq. (6b)
   cheat_n10  = milli*sqrt_cdrag_n10*(18.0_rsh*stab + 32.7_rsh*(1-stab))  !   L & Y eq. (6c), (6d)
!!
!! Initializing transfert coefficients with their first guess neutral equivalents :
   cdrag = cdrag_n10 ;  cevap = cevap_n10 ;  cheat = cheat_n10 ;  sqrt_cdrag = sqrt(cdrag)

!! * Now starting iteration loop
   DO j_itt=1, nb_itt
!! Turbulent scales :
      U_star  = sqrt_cdrag*dU10                !   L & Y eq. (7a)
      T_star  = cheat/sqrt_cdrag*dT            !   L & Y eq. (7b)
      q_star  = cevap/sqrt_cdrag*dq            !   L & Y eq. (7c)

!! Estimate the Monin-Obukov length      ! L & Y eq. (8a)
      L  = (U_star**2)/( karman*grav*(T_star/T_vpot + q_star/(q_a + 1.0_rsh/0.608_rsh)) )

!! Stability parameters :
      zeta = zu/L ;  zeta   = sign( min(abs(zeta),10.0_rsh), zeta )
      zpsi_h = psi_h(zeta)
      zpsi_m = psi_m(zeta)

!! Shifting the wind speed to 10m and neutral stability :
      U_n10 = dU10*1.0_rsh/(1.0_rsh + sqrt_cdrag_n10/karman*(log(zu/10.0_rsh) - zpsi_m)) !  L & Y eq. (9a)

!! Updating the neutral 10m transfer coefficients :
      cdrag_n10  = milli * (2.7_rsh/U_n10 + 0.142_rsh + U_n10/13.09_rsh)     !  L & Y eq. (6a)
      sqrt_cdrag_n10 = sqrt(cdrag_n10)
      cevap_n10  = milli * (34.6_rsh * sqrt_cdrag_n10)                       !  L & Y eq. (6b)
      stab    = 0.5_rsh + sign(0.5_rsh,zeta)
      cheat_n10  = milli*sqrt_cdrag_n10*(18.0_rsh*stab + 32.7_rsh*(1-stab))  !  L & Y eq. (6c), (6d)

!! Shifting the neutral  10m transfer coefficients to ( zu , zeta ) :
!!
      xct = 1.0_rsh + sqrt_cdrag_n10/karman*(log(zu/10) - zpsi_m)
      cdrag  = cdrag_n10/(xct*xct) ;  sqrt_cdrag = sqrt(cdrag)     !  L & Y eq. (10a)
!!
      xlogt = log(zu/10.0_rsh) - zpsi_h
!!
      xct = 1.0_rsh + cheat_n10*xlogt/karman/sqrt_cdrag_n10
      cheat  = cheat_n10*sqrt_cdrag/sqrt_cdrag_n10/xct             !  L & Y eq. (10b)
!!
      xct = 1.0_rsh + cevap_n10*xlogt/karman/sqrt_cdrag_n10
      cevap  = cevap_n10*sqrt_cdrag/sqrt_cdrag_n10/xct             !  L & Y eq. (10c)
!!
   END DO
!!
  END SUBROUTINE coef_large_1z

!!======================================================================

  SUBROUTINE coef_large_2z(zt, zu, sst, T_zt, q_w, q_zt, dU, cdrag, cheat, cevap, T_zu, q_zu)
!!----------------------------------------------------------------------
!!                      ***  ROUTINE  turb_core  ***
!!
!! ** Purpose :   Computes turbulent transfert coefficients of surface
!!                fluxes according to Large & Yeager (2004).
!!
!! ** Method  :   I N E R T I A L   D I S S I P A T I O N   M E T H O D
!!      Momentum, Latent and sensible heat exchange coefficients
!!      Caution: this procedure should only be used in cases when air
!!      temperature (T_air) and air specific humidity (q_air) are at 2m
!!      whereas wind (dU) is at 10m.
!!
!! References :
!!      Large & Yeager, 2004 : ???
!! History :
!!   9.0  !  06-12  (L. Brodeau) Original code for 2Z
!!----------------------------------------------------------------------
!! * Modules used
   USE parameters, ONLY : grav

!! * Arguments
   REAL(KIND=rsh), INTENT(in)   :: &
         zt,      &     ! height for T_zt and q_zt                   [m]
         zu             ! height for dU                              [m]
   REAL(KIND=rsh), INTENT(in)                                   ::  &
         sst,     &     ! sea surface temperature              [Kelvin]
         T_zt,    &     ! air potential temperature            [Kelvin]
         q_w,     &     ! sea surface specific humidity         [kg/kg]
         q_zt,    &     ! specific air humidity                 [kg/kg]
         dU             ! relative wind module |U(zu)-U(0)|       [m/s]
   REAL(KIND=rsh), INTENT(out)                                    ::  &
         cdrag,   &     ! transfer coefficient for momentum         (tau)
         cheat,   &     ! transfer coefficient for sensible heat (Q_sens)
         cevap,   &     ! transfert coefficient for evaporation   (Q_lat)
         T_zu,    &     ! air temp. shifted at zu                     [K]
         q_zu           ! spec. hum.  shifted at zu               [kg/kg]

!! * Local declarations
   REAL(KIND=rsh), PARAMETER        :: milli = 0.001_rsh, &
                                       karman  = 0.4_rsh      ! von Karman s constant
   REAL(KIND=rsh)                                   ::  &
         dU10,           & ! relative wind dU                  [m/s]
         dT,             & ! air/sea temperature differeence   [K]
         dq,             & ! air/sea humidity difference       [kg/kg]
         cdrag_n10,      & ! 10m neutral drag coefficient
         cevap_n10,      & ! 10m neutral latent coefficient
         cheat_n10,      & ! 10m neutral sensible coefficient
         sqrt_cdrag_n10, & ! root square of cdrag_n10
         sqrt_cdrag,     & ! root square of cdrag
         T_vpot_u,       & ! virtual potential temperature        [K]
         T_star,         & ! turbulent scale of tem. fluct.
         q_star,         & ! turbulent humidity of temp. fluct.
         U_star,         & ! turb. scale of velocity fluct.
         L,              & ! Monin-Obukov length                  [m]
         zeta_u,         & ! stability parameter at height zu
         zeta_t,         & ! stability parameter at height zt
         U_n10,          & ! neutral wind velocity at 10m        [m]
         xlogt, xct, zpsi_hu, zpsi_ht, zpsi_m

   INTEGER :: j_itt
   INTEGER, PARAMETER :: nb_itt = 3   ! number of itterations
   INTEGER                                   :: &
           &     stab                ! 1st stability test integer

!!  * Start

!! Initial air/sea differences
   dU10 = max(0.5_rsh, dU)      !  we don t want to fall under 0.5 m/s
   dT = T_zt - sst 
   dq = q_zt - q_w

!! Neutral Drag Coefficient :
   stab = 0.5_rsh + sign(0.5_rsh,dT)                 ! stab = 1  if dT > 0  -> STABLE
   cdrag_n10  = milli*( 2.7_rsh/dU10 + 0.142_rsh + dU10/13.09_rsh )          !   L & Y eq. (6a)
   sqrt_cdrag_n10 = sqrt(cdrag_n10)
   cevap_n10  = milli*( 34.6_rsh * sqrt_cdrag_n10 )                          !   L & Y eq. (6b)
   cheat_n10  = milli*sqrt_cdrag_n10*(18.0_rsh*stab + 32.7_rsh*(1 - stab))   !   L & Y eq. (6c, 6d)

!! Initializing transf. coeff. with their first guess neutral equivalents :
   cdrag = cdrag_n10 ;  cevap = cevap_n10 ;  cheat = cheat_n10 ;  sqrt_cdrag = sqrt(cdrag)

!! Initializing z_u values with z_t values :
   T_zu = T_zt ;  q_zu = q_zt

!!  * Now starting iteration loop
   DO j_itt=1, nb_itt
      dT = T_zu - sst ;  dq = q_zu - q_w ! Updating air/sea differences
      T_vpot_u = T_zu*(1.0_rsh + 0.608_rsh*q_zu)    ! Updating virtual potential temperature at zu
      U_star = sqrt_cdrag*dU10                ! Updating turbulent scales :   (L & Y eq. (7a))
      T_star  = cheat/sqrt_cdrag*dT           !   L & Y eq. (7b)
      q_star  = cevap/sqrt_cdrag*dq           !   L & Y eq. (7c)
!!

      IF (T_star + 0.608_rsh*T_zu*q_star /= 0.0_rsh) THEN
!! Estimate the Monin-Obukov length at height zu    ! L & Y eq. (8a)
!L  = (U_star**2)/( karman*grav*(T_star/T_vpot_u + q_star/(q_zu + 1.0_rsh/0.608_rsh)) )
!FAUXL = (U_star**2) / (karman*grav/T_vpot_u*(T_star*(1.0_rsh+0.608_rsh*q_zu) + 0.608_rsh*T_zu*q_star))
        L = (U_star**2) / (karman*grav/T_vpot_u*(T_star + 0.608_rsh*T_zu*q_star))
!! Stability parameters :
        zeta_u  = zu/L  ;  zeta_u = sign( min(abs(zeta_u),10.0_rsh), zeta_u )
        zeta_t  = zt/L  ;  zeta_t = sign( min(abs(zeta_t),10.0_rsh), zeta_t )
      ELSE
! not really clean but... it happens very rarely and by luck not at the last iteration ?
        zeta_u  = 0.0_rsh
        zeta_t  = 0.0_rsh
      END IF
      zpsi_hu = psi_h(zeta_u)     !   L & Y eq. (8c, 8e)
      zpsi_ht = psi_h(zeta_t)     !   L & Y eq. (8c, 8e)
      zpsi_m  = psi_m(zeta_u)     !   L & Y eq. (8c, 8d)
!!
!! Shifting the wind speed to 10m and neutral stability : (L & Y eq.(9a))
      U_n10 = dU10/(1.0_rsh + sqrt_cdrag_n10/karman*(log(zu/10.0_rsh) - zpsi_m))
!!
!! Shifting temperature and humidity at zu :          (L & Y eq. (9b-9c))
!! Main difference with coef_large_1z !!!
      T_zu = T_zt - T_star/karman*(log(zt/zu) + zpsi_hu - zpsi_ht)
      q_zu = q_zt - q_star/karman*(log(zt/zu) + zpsi_hu - zpsi_ht)
!!
!! q_zu cannot have a negative value : forcing 0
      stab = 0.5_rsh + sign(0.5_rsh,q_zu) ;  q_zu = stab*q_zu
!!
!! Updating the neutral 10m transfer coefficients :
      cdrag_n10  = milli * (2.7_rsh/U_n10 + 0.142_rsh + U_n10/13.09_rsh)    ! L & Y eq. (6a)
      sqrt_cdrag_n10 = sqrt(cdrag_n10)
      cevap_n10  = milli * (34.6_rsh * sqrt_cdrag_n10)                 ! L & Y eq. (6b)
      stab    = 0.5_rsh + sign(0.5_rsh,zeta_u)
      cheat_n10  = milli*sqrt_cdrag_n10*(18.0_rsh*stab + 32.7_rsh*(1-stab)) ! L & Y eq. (6c-6d)
!!
!!
!! Shifting the neutral 10m transfer coefficients to (zu,zeta_u) :
      xct = 1.0_rsh + sqrt_cdrag_n10/karman*(log(zu/10.0_rsh) - zpsi_m)
      cdrag = cdrag_n10/(xct*xct) ; sqrt_cdrag = sqrt(cdrag)     !  L & Y eq. (10a)
!!
      xlogt = log(zu/10.0_rsh) - zpsi_hu
!!
      xct = 1.0_rsh + cheat_n10*xlogt/karman/sqrt_cdrag_n10
      cheat  = cheat_n10*sqrt_cdrag/sqrt_cdrag_n10/xct           !  L & Y eq. (10b)
!!
      xct = 1.0_rsh + cevap_n10*xlogt/karman/sqrt_cdrag_n10
      cevap  = cevap_n10*sqrt_cdrag/sqrt_cdrag_n10/xct           !  L & Y eq. (10c)
!!
!!
   END DO
  END SUBROUTINE coef_large_2z

!!======================================================================

  FUNCTION psi_m(zta)
!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE psi_m  ***
!&E
!&E ** Purpose : evaluate stability function for wind speed and scalars
!&E              for the momentum variables
!&E
!&E ** Description : matches Kansas and free convection forms with weighting
!&E             convective form follows Fairall et al (1996) with profile constants
!&E             from Grachev et al (2000) BLM
!&E             stable form from Beljaars and Holtslag (1991)
!&E
!&E ** Called by : sflx_turb_large
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables :
!&E
!&E ** Reference : Psis, L & Y eq. (8c), (8d), (8e)
!&E
!&E ** History :
!&E          issued from OPA 9.0
!&E       !  2009-06  (V. Garnier) MARS implementation
!&E
!&E---------------------------------------------------------------------
!! * Modules used

   REAL(KIND=rsh),                                   INTENT(in) :: zta

   REAL(KIND=rsh), PARAMETER :: pi = 3.14159_rsh
   REAL(KIND=rsh)                                               :: psi_m
   REAL(KIND=rsh)                                               :: X2, X, stabit
   X2 = sqrt(abs(1.0_rsh - 16.0_rsh*zta))  ;  X2 = max(X2 , 1.0_rsh) ;  X  = sqrt(X2)
   stabit    = 0.5_rsh + sign(0.5_rsh,zta)
   psi_m = -5.0_rsh*zta*stabit  &                                                  ! Stable
        & + (1 - stabit)*(2.0_rsh*log((1.0_rsh + X)/2.0_rsh) + log((1.0_rsh + X2)/2.0_rsh) - 2.0_rsh*atan(X) + pi/2.0_rsh)  ! Unstable
  END FUNCTION psi_m

!!======================================================================

  FUNCTION psi_h(zta)
!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE psi_h  ***
!&E
!&E ** Purpose : evaluate stability function for wind speed and scalars
!&E              for the heat fluxes variables
!&E
!&E ** Description : matches Kansas and free convection forms with weighting
!&E             convective form follows Fairall et al (1996) with profile constants
!&E             from Grachev et al (2000) BLM
!&E             stable form from Beljaars and Holtslag (1991)
!&E
!&E ** Called by : sflx_turb_large
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables :
!&E
!&E ** Reference : Psis, L & Y eq. (8c), (8d), (8e)
!&E
!&E ** History :
!&E          issued from OPA 9.0
!&E       !  2009-06  (V. Garnier) MARS implementation
!&E
!&E---------------------------------------------------------------------
!! * Modules used

   REAL(KIND=rsh),                                   INTENT(in) :: zta

   REAL(KIND=rsh)                                               :: psi_h
   REAL(KIND=rsh)                                               :: X2, X, stabit
   X2 = sqrt(abs(1.0_rsh - 16.0_rsh*zta))  ;  X2 = max(X2 , 1.0_rsh) ;  X  = sqrt(X2)
   stabit    = 0.5_rsh + sign(0.5_rsh,zta)
   psi_h = -5.0_rsh*zta*stabit  &                                       ! Stable
        & + (1 - stabit)*(2.0_rsh*log( (1.0_rsh + X2)/2.0_rsh ))        ! Unstable
    END FUNCTION psi_h


!!======================================================================

!!======================================================================

  SUBROUTINE sflx_turb_ayina(ztemp,zsat,zwindmod,zpa,zrh,zrad,zirlos,zrain, &
                             zsenslos,zlatlos,zevap,ztosx,ztosy)
                              

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE sflx_turb_ayina  ***
!&E
!&E ** Purpose : Computation of the heat fluxes and wind stress
!&E              coefficients at the ocean surface
!&E
!&E ** Description : The Wind stress and the flux of heat at
!&E              the ocean surfaces are derived from
!&E              semi-empirical ( or bulk ) formulae which
!&E              relate the flux to the properties of the surface
!&E              and of the lower atmosphere.
!&E
!&E              zpa      ! atmospherical pressure             [Pa]
!&E              zsat     ! air temperature                    [Celsius]
!&E              ztemp    ! sea surface temperature            [Celsius]
!&E              zwindmod ! module of wind at 10 m             [m/s]
!&E              zrh      ! relative humidity of air           [kg/kg]
!&E              suswe    ! relative humidity of the ocean     [kg/kg]
!&E
!&E ** Called by : sflxsurf.F90
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables :
!&E
!&E ** Reference : Ayina et Bentamy (2007), Ayina et al (2008)
!&E                METEOCEAN report (2009)
!&E
!&E ** History :
!&E       !  2003-10  (L.-H. Ayina)  Original code
!&E       !  2009-06  (V. Garnier) MARS implementation
!&E
!&E----------------------------------------------------------------------
   USE comvars2d,  ONLY : jb,jh,ig,id,nbouc
!   USE cometdyn,   ONLY : tosx, tosy, u10, v10

!! * Arguments
   REAL(KIND=rsh), INTENT( IN )    :: ztemp,zsat,zwindmod,zpa,zrh, &
                                      zrad,zirlos,zrain
   REAL(KIND=rsh), INTENT( OUT )   :: zsenslos,zlatlos,zevap
   REAL(KIND=rsh), INTENT( INOUT ) :: ztosx,ztosy

!! * Local declarations


   REAL(KIND=rsh), PARAMETER ::          &     ! wind scaling
                        cv0=0.0_rsh, cv2=2.0_rsh, cv3=3.0_rsh,       &
                        cv4=4.0_rsh, cv5=5.0_rsh, cv6=6.0_rsh,       &
                        cv7=7.0_rsh, cv8=8.0_rsh, cv9=9.0_rsh,       &
                        cv10=10.0_rsh, cv11=11.0_rsh, cv13=13.0_rsh, &
                        cv14=14.0_rsh, cv18=18.0_rsh, cv19=19.0_rsh, &
                        cv24=24.0_rsh, cv25=25.0_rsh, cv52=52.0_rsh

   REAL(KIND=rsh), PARAMETER ::    &
                        zal=0.032_rsh,     &
                        zgr= 10.0_rsh,     &
                        zga= 287.04_rsh,   &  !   []
                        zc1= 0.7859_rsh,   &
                        zc2= 0.03477_rsh,  &
                        zc3= 0.00412_rsh,  &
                        zc4= 0.622_rsh,    &
                        zc5= 1.0_rsh-zc4,  &
                        zc6= 0.61_rsh,     &
                        zc8= 2.2e-3,       &
                        zzv= 10.0_rsh,     &  ! def  [m]
                        zca= 0.4_rsh,      &
                        zm1= 13.2_rsh,     &
                        zm2= 0.0_rsh,      &
                        zm3= 1.e-06,       &
                        zco= 1.e-04_rsh,   &
                        cqs= 1.32e-3,      &
                        cql= 1.15e-3,      &
                        zcodu= 1.e-3,      &
                        zcodt= 1.e-3,      &
                        zcodq = 1.e-3,     &
                        zhlw= 2.5e+06,     &
                        zhse= 1004.0_rsh,  &
                        zzero= 0.0_rsh,    &
                        krenf= 1.2_rsh,    &
                        zdept= 10.0_rsh,   &
                        zvita= 1.0_rsh,    &
                        zvito= 0.5_rsh,    &
                        zzone= 1.0_rsh,    &
                        coest= -0.005_rsh, &
                        coein= -0.005_rsh, &
                        a1= -0.146785_rsh, &
                        b1= -0.292400_rsh, &
                        c1= -2.207748_rsh, &
                        d1=  1.612292_rsh, &
                        csss=0.527_rsh

   INTEGER zitime

   REAL(KIND=rsh)    ::   &
                        zu_star, zz_utas, zz_usta,               &
                        zu_plog, zu_rapu, zu_dist,zt_sta1,       &
                        zt_sta2,zt_dist , zq_sta1,zq_sta2,       &
                        zq_dist,zu_refer, zq_evapo,              &
                        zq_hlate,zt_hsens, zu_zotes,zlmob,       &
                        zt_artz, zt_usta, zt_plog,               &
                        zu_prlog,zq_artz,zq_usta,                &
                        zq_plog, zkav,zdvp,zkov,zkot,            &
                        zkaq,dift,difq,zwind,                    &
                        zvis,zcin,zdaq,zkoq,zdiq,zdit,           &
                        nkov,nkot,nkoq,ztax,ztay,zten,           &
                        ztqa,zesa,zaqa,zrha,ztva,ztpa,ztka,zvpa, &

                        ztqo,zeso,zaqo,ztvo,ztpo,ztko,zvpo,ztmo, &
                        zcte,zkma,zcoq,zkmi,                     &
                        zcqe,zkat,zcor,zobmi,zobms,              &

                        st_zpsiv,st_zpsit,st_zpsiq,              &
                        in_zpsiv,in_ta,in_zpsit,st_zobmi,        &
                        in_zpsiq,in_zobmi,in_zx,in_l1,in_l2,     &
                        zlma,zpsiv,zpsit,zpsiq

!!----------------------------------------------------------------------
!! * Executable part


!------------------------------------------------------------------!
!  computation of wind direction                                   !
!------------------------------------------------------------------!

         zwind  = MAX(zwindmod,0.5_rsh)          ! Warning : a minimum is required [m/s]
         zvis   = zm1 + zm2*zsat                 ! [degrees Celsius] []
         zcin   = zvis*zm3                       ! kinematic viscosity []

!------------------------------------------------------------------!
!   Atmospherical parameters                                       !
!------------------------------------------------------------------!

         ztqa   = (zc1+zc2*zsat)/(1.0_rsh+zc3*zsat)
         ztqa   = 10.0_rsh**ztqa                     !  saturation vapor pressure [mbar]
         ztqa   = 100.0_rsh*ztqa                     !  [mbr] -----> [Pa]
         ztqa   = zrh*ztqa                           !  relative humidity
         zesa   = zc4*ztqa/(zpa-zc5*ztqa)            !  specific humidity

         ztka   = zsat + kelvin           ! Air temperature [Kelvin]
         ztpa   = ztka + zt*0.0098_rsh    ! Potential temperature
! the air temperature increases by 0.98 degrees Celsius per 100 meters in adiabatic conditions

         ztva   =  ztka*(1.0_rsh+zc6*zesa)    ! Virtual Temperature

         zvpa   =  ztva                  ! Virtual Potential temperature
! normalement on doit ajouter un alpha
! qui est egale zt*0.0098. ATTENTION approximation
         zaqa   =  ztqa/(zga*ztka)       ! Absolue Humidity    []
         zrha   =  zpa/(zga*ztka)        ! air humidity density -> CLIO  []

!------------------------------------------------------------------!
!   Ocean parameters                                               !
!------------------------------------------------------------------!

         ztqo = (zc1+zc2*ztemp)/(1.0_rsh+zc3*ztemp)
         ztqo = 10.0_rsh**ztqo                  ! saturation vapor pressure [mbar]
         ztqo = 100.0_rsh*ztqo                  ! [mbr] -----> [Pa]
         ztqo = suswe*ztqo                      ! salinity effect
         zeso = zc4*ztqo/(zpa-zc5*ztqo)         ! specific humidity -> CLIO

         ztko = ztemp + kelvin                  ! Sea surface temperature [Kelvin]
         ztpo = ztko                            ! Potential temperature  ATTENTION approximation
         ztvo = ztko*(1.0_rsh+zc6*zeso)         ! Virtual Temperature
         zvpo = ztvo                            ! Virtual Potential temperature ATTENTION approximation
         zaqo = ztqo/(zga*ztko)                 ! Absolue Humidity -> CLIO ??? []

!-------------------------------------------------------------------!
!      Computation of stability criterios                           !
!                               : terme de gradient adiabatique     !
!-------------------------------------------------------------------!

! difference of humidity [] [kg/kg](au lieu
         zdiq   =   zesa - zeso          ! de prendre la dif d humidite absolue
! je vais prendre celle de l humidite
! specifique)
         zdaq   =   zeso - zesa

         zdit   =   ztpo - ztpa          ! difference of potential temperatures  [K]

         zdvp   =   zvpo - zvpa          ! difference of virtual potential temperatures [K]

         ztmo   =   ztka*(1.0_rsh + zc8*ztka*zesa)

!------------------------------------------------------------------!
!   Initialization parameters with Smith et al . 1998 (neutral)    !
!------------------------------------------------------------------!

         zkot = 1.3_rsh    ! Drag coefficient []
         zkov = 1.3_rsh

         IF ((zwind.ge.cv0).AND.(zwind.le.cv52)) THEN
           zkot = a1*EXP(b1*(zwind+c1))+(d1/zwind)+1.0_rsh    ! Drag coefficient (smith)
         ENDIF

         IF ((zwind.ge.cv0).AND.(zwind.le.cv2)) THEN
           zkov = 1.0_rsh-0.01_rsh*zwind
         ELSE IF ((zwind.ge.cv3).AND.(zwind.le.cv4)) THEN
           zkov =  0.91_rsh+0.02_rsh*zwind
         ELSE IF ((zwind.ge.cv5).AND.(zwind.le.cv6)) THEN
           zkov =  0.83_rsh+0.04_rsh*zwind
         ELSE IF ((zwind.ge.cv7).AND.(zwind.le.cv8)) THEN
           zkov =  0.71_rsh+0.06_rsh*zwind
         ELSE IF ((zwind.ge.cv9).AND.(zwind.le.cv10)) THEN
           zkov =  0.70_rsh+0.06_rsh*zwind
         ELSE IF ((zwind.ge.cv11).AND.(zwind.le.cv13)) THEN
           zkov =  0.80_rsh+0.05_rsh*zwind
         ELSE IF ((zwind.ge.cv14).AND.(zwind.le.cv18)) THEN
           zkov =  0.81_rsh+0.05_rsh*zwind
         ELSE IF ((zwind.ge.cv19).AND.(zwind.le.cv24)) THEN
           zkov =  0.80_rsh+0.05_rsh*zwind
         ELSE IF ((zwind.ge.cv25).AND.(zwind.le.cv52)) THEN
           zkov =  1.04_rsh+0.04_rsh*zwind
         END IF

         zkov  = zkov*1.e-3
         zkot  = zkot*1.e-3
         zkoq  = zkot
         zkav  = zkov
         zkat  = zkot
         zkaq  = zkot
         nkov  = zkov
         nkot  = zkot
         nkoq  = zkoq
         zkma  = 2.9e-3
         zkmi  = 1.1e-3


!------------------------------------------------------------------!
!   Obukhov  parameters                                            !
!------------------------------------------------------------------!

         IF ((ABS(zdiq).ge.zzero).AND.(ABS(zdit).ge.zzero)) THEN

           zitime = 0
           zu_dist = 1.0_rsh
           zt_dist = 1.0_rsh
           zq_dist = 1.0_rsh

!---------------------------------------------------------------------!
! Controle de la coherence entre les deux termes definisants les 3 flux!
!---------------------------------------------------------------------!
           DO WHILE ((zu_dist > zcodu) .OR. (zt_dist > zcodt) .OR. (zq_dist > zcodq))

             zitime = zitime + 1

             zcor    =   zwind*zwind*ztva      ! U*U*Tv
             zcor    =   zcor*(zkov**1.5_rsh)  ! U*U*TvCd***1.5
             zcor    =   zca*zgr*zzv*zkot/zcor ! K*Z*g*Ct/U*U*Tv*Cd***1.5

             zcte    =   zdit*zcor

             zcoq    =   zc6*ztka*zkoq         !
             zcoq    =   zcoq/(zkot)
             zcqe    =   zcoq*zdaq*zcor        !

             zobms   =   coest*(zcte + zcqe)    !
             zobmi   =   coein*(zcte + zcqe)    ! (Z/L) parametre de stabilite
             zlmob   =   10.0_rsh/coein*(zcte + zcqe)

!------------------------------------------------------------------!
!   Stable conditions                                              !
!------------------------------------------------------------------!

             st_zobmi  = MAX(zzero,zobms)
             st_zpsiv  =  -7.0_rsh * st_zobmi  ! stability function for wind
             st_zpsit  =  st_zpsiv             ! stability function for hsensible
             st_zpsiq  =  st_zpsiv             ! stability function for hlatent

!------------------------------------------------------------------!
!   Instable conditions                                            !
!------------------------------------------------------------------!

             in_zobmi  = MIN(zzero,zobmi)
             in_zx     =  (1.0_rsh - 16.0_rsh * in_zobmi)**0.25_rsh
             in_l1     =  2.0_rsh * LOG((1.0_rsh+in_zx)/2.0_rsh)
             in_l2     =  LOG((1.0_rsh+ in_zx*in_zx)/2.0_rsh)
             in_ta     =  2.0_rsh * ATAN(in_zx)
             in_zpsiv  =  in_l1 + in_l2 - in_ta + 3.1416_rsh/2.0_rsh
             in_zpsit  =  2.0_rsh*in_l2
             in_zpsiq  =  in_zpsit

!------------------------------------------------------------------!
!   Computation of stability function                              !
!------------------------------------------------------------------!

             zlma     =   MAX(zzero,SIGN(zzone,zdvp)) ! zlma = 0. instable
!      = 1. stable
             zpsiv    =   zlma*in_zpsiv + (1.0_rsh - zlma)*st_zpsiv
             zpsit    =   zlma*in_zpsit + (1.0_rsh - zlma)*st_zpsit
             zpsiq    =   zpsit

!------------------------------------------------------------------!
!  calcul de u_star    : correction sur le cd et le ustatr         !
!------------------------------------------------------------------!

             zu_star  =   (zkov**0.5_rsh)*zwind         ! old u*
             zz_utas  =   zal*zu_star*zu_star/zgr       ! z* by charnock
             zz_usta  =   zz_utas +  zcin/(zgr*zu_star) ! z*+effet de faible vent

             zu_plog  =   (LOG(zzv/zz_usta)-zpsiv)/zca  ! profil log de u
             zu_rapu  =   zwind/zu_star                 ! rapport des u
             zu_star  =   zu_star*zu_rapu/zu_plog       ! new u* corrige
             zkov  =   (zu_star/zwind)**2               ! new coef de transfere
             zu_dist  =   zu_plog - zu_rapu             ! diff of news profils
             zu_dist  =   ABS(zu_dist)

             IF (zkov > zkma)   zkov = nkov + zkmi
             IF (zkov < zkmi)   zkov = nkov*krenf

!-------------------------------------------------------------------!
!  calcul de  Q_star : correction sur la temperature et le ct       !
!-------------------------------------------------------------------!

             IF (zitime == 1) then
               zt_artz   =   zca*zca/(zkot*LOG(10.0_rsh/zz_usta))
               zt_usta   =   10.0_rsh/EXP(zt_artz)
               zt_plog   =   LOG(zt/zt_usta)-zpsit
             ENDIF

             zt_sta1   =   1.0_rsh/(LOG(zzv/zt_usta)-zpsit)
             zt_sta2   =   zkot*zwind/(zca*zu_star)
             zkot   =   zkot*zt_sta1/zt_sta2
             zt_dist   =   zdit*(zt_sta1 - zt_sta2)
             zt_dist   =   ABS(zt_dist)

             IF (zkot > zkma) zkot = nkot + zkmi
             IF (zkot < zkmi) zkot = nkot*krenf

!-------------------------------------------------------------------!
!  calcul de q_star     : correction sur lhumidité  et le c        !
!-------------------------------------------------------------------!

             IF (zitime == 1) THEN
               zq_artz   =   zca*zca/(zkoq*LOG(10.0_rsh/zz_usta))
               zq_usta   =   10.0_rsh/EXP(zq_artz)
               zq_plog   =   LOG(zt/zq_usta)-zpsiq
             ENDIF

             zq_sta1   =   1.0_rsh/(LOG(zzv/zq_usta)-zpsiq)
             zq_sta2   =   zkoq*zwind/(zca*zu_star)
             zkoq   =   zkoq*zq_sta1/zq_sta2
             zq_dist   =   zdiq*(zq_sta1 - zq_sta2)
             zq_dist   =   ABS(zq_dist)

             IF (zkoq > zkma)   zkoq = nkov + zkmi
             IF (zkoq < zkmi)   zkoq = nkoq*krenf

!-------------------------------------------------------------------!
!  cas non convergent                                               !
!-------------------------------------------------------------------!

             IF (zitime == 10)  then
               zt_dist = 0.00001_rsh    ! convergence criterium is satisfied
               zq_dist = 0.00001_rsh    ! convergence criterium is satisfied
               zu_dist = 0.00001_rsh    ! convergence criterium is satisfied
               zkot    = nkot           ! keep the values of the last iteration
               zkoq    = nkoq
               zkov    = nkov
               dift    = zdit           ! difference of temperature
               difq    = zdaq
             ENDIF

           END DO

!---------------------------------------------------------------------!
!       calcul de Long Obukhov & Monin et de U_ref                    !
!---------------------------------------------------------------------!

           zu_prlog = (LOG(zzv/zz_usta) - zpsiv)/zca
           zu_refer = zu_star*zu_prlog
           zu_refer = zwind

         ELSE                                          ! cas neutre et stable

           zkoq   =   nkoq
           zkot   =   nkot
           zkov   =   nkov
           zu_refer  =  zwind
           zu_star   =  (zkov**0.5_rsh)*zwind
           zz_utas   =  zal*zu_star*zu_star/zgr          ! z* by charnock
           zz_usta   =  zz_utas + zcin/(zgr*zu_star)     ! z*+effet de faible vent
           zt_sta2   =  zdit*zkot*zwind/(zca*zu_star)
           zq_sta2   =  zdaq*zkoq*zwind/(zca*zu_star)
           zu_plog  =   LOG(zzv/zz_usta)/zca
           zkov  =   (1.0_rsh/zu_plog)**2
           zkoq   =   nkoq
           zkot   =   nkot

         ENDIF

!---------------------------------------------------------------------!
! Calcul des trois flux : tension de vent,                            !
!                        flux de chaleur latente et sensible          !
!---------------------------------------------------------------------!

         zq_evapo  = zrha*zkoq*zu_refer*(zeso-zesa)        ! evaporation m/s
         zq_hlate  = - zhlw*zq_evapo                       ! latent heat flux

         zt_hsens  = - zrha*zhse*zkot*zu_refer*(ztko-ztka) ! sensible heat flux

         zu_zotes  = zrha*zkov                        ! wind
         zu_zotes  = zu_zotes*zu_refer
         zu_zotes  = zu_zotes*zu_refer                ! wind stress

! Wind stress
         ztosx=zrha*zkov*ztosx
         ztosy=zrha*zkov*ztosy

! zu_refer is the modified module of the wind
! from which all the exchange coefficients have been estimated.
! ztosx = ||V||*Vx ; ztosy = ||V||*Vy, so zu_refer/zwindmod removes
! the module of the wind in the formula,
! the second factor zu_refer/zwindmod corrects the component of the wind

         IF (zwindmod > 0.0_rsh) THEN
           ztosx=ztosx*(zu_refer/zwindmod)**2
           ztosy=ztosy*(zu_refer/zwindmod)**2
         ENDIF




!---------------------------------------------------------------------!
! Assign out variables                                                !
!---------------------------------------------------------------------!

         zsenslos = zt_hsens
         zlatlos  = zq_hlate
         zevap    = zq_evapo


  END SUBROUTINE sflx_turb_ayina

!!======================================================================


!!======================================================================

  SUBROUTINE sflx_turb_fairall(ztemp,zsat,zwindmod,zpa,zrh,zrad,zirlos,zrain, &
                               zsenslos,zlatlos,zevap,ztosx,ztosy)

!&E--------------------------------------------------------------------
!&E                 ***  ROUTINE sflx_turb_fairall  ***
!&E
!&E ** Purpose : Computation of the heat fluxes and wind stress
!&E              coefficients at the ocean surface
!&E
!&E ** Description :
!&E        COARE Bulk Flux Algorithm version 3.0a. See Fairall et al.,2003:
!&E        J.Climate,16,571-591
!&E
!&E        Written 05/09/2003 by Frank Bradley frank.bradley@csiro.au
!&E
!&E        Some details of changes from versions 2.5 and 2.6.
!&E        For full description see Cor30_readme.txt
!&E
!&E        Wave options for Zo - models by Taylor and Yelland
!&E        (JPO 2001,31,572) or Oost et al.(BLM 2002,103,409)
!&E        Wind speed dependant Charnock parameter above 10 m/s
!&E        Grachev and Fairall first guess z/L routine (JAM 36,406-414,1997)
!&E        to cut iterations from 20 to 3
!&E        rr-zoq,charnock relations as described in COARE 3.0 J. Climate paper
!&E        (see above) Use separate phiu,phit functions with Grachev et al.
!&E        (2000) constants Stable functions according to Beljaars and Holtslag
!&E        (JAM 30, 327-341,1991) Return Webb Wbar because Webb flux correction already
!&E        included in latent heat flux Set up to read COARE test data set, test3_0.txt
!&E        with 4-digit year tab delimited see Bradley, Moncrieff and Weller,
!&E        Proceedings of the COARE Flux Group
!&E        workshop, May 14-16 1997,c NCAR/MMM, PO Box 3000, Boulder, CO 80307-3000, USA
!&E
!&E        THIS version "f" uses Chris Fairall s seasnake for Ts nominal depth 0.05m
!&E        Another version "h" uses Mike Gregg s ASP Ts data at 6m depth, prepared
!&E        by Hemantha Wijesekera Outputs tst3_0af.out or tst3_0ah.out written from
!&E        subroutine bulk_flux
!&E        Wave option outputs tst30afo.out or tst30aft.out (Jwave = 1 or 2)
!&E
!&E ** Input data - introduced in calling program "fluxes" and passed to subroutines in COMMON
!&E        xtime (COARE convention yyyymnddhhmmss.ss GMT)
!&E        u (wind speed relative to the sea surface) m/s  i.e. add CURRENT VECTOR if available
!&E        ts=temp (sea temperature)  deg. C
!&E        t=sat (air temperature) deg. C
!&E        q (specific humidity) g/kg or RH (as decimal) - code detects which then works internally in kg/kg
!&E        rs (shortwave radiation) W/m2
!&E        rl (downwelling longwave) W/m2
!&E        zrain (average rainrate in timestep) mm/hour
!&E        xlat (latitude) degrees  [latitude north +ve, latitude south -ve]
!&E        xlon (longitude) degrees [longitude east +ve, longitude west -ve]
!&E        hwt alternative sea temperature (6m depth)
!&E        p (pressure) mb; use 1008mb if data unavailable
!&E        zi (boundary-layer depth) m; use 600m if data unavailable
!&E        zu (wind measurement height) m
!&E        zt (T measurement height) m
!&E        zq (q measurement height) m
!   !&E        zus (wind standard height) m - usually 10m
!   !&E        zts (T and q standard height) m - usually 2m
!   !&E        ts_depth (depth of ts measurement) - positive m - for warm layer correction
!&E        jcool (=1 for cool skin calculation;  =0 if SST measured by IR radiometer)
!&E        jwarm (=1 for warm layer calculation; =0 if SST measured by IR radiometer)
!&E        jwave (wave state options 0=Charnock,1=Oost et al,2=Taylor and Yelland)
!&E
!&E ** Calculated inputs:
!&E        qs (humidity at ocean interface)
!&E        tsw (ts with warm layer added)
!&E
!&E ** Outputs available - write statements in subroutine bulk_flux
!&E        hf W/m**2     (Sensible heat flux) - turbulent part only
!&E        ef W/m**2     (Latent heat flux) - turbulent part only
!&E        rf W/m**2     (Heat flux due to rainfall)
!&E        tau N/m**2   (wind stress) - turbulent part only
!&E        taur - momentum flux due to rain N/m**2
!&E        usr m/s     (M-O velocity scaling parameter u* = friction velocity)
!&E        tsr C       (M-O temperature scaling parameter t*)
!&E        qsr kg/kg   (M-O humidity scaling parameter q*)
!&E        Cdn - neutral drag coefficient
!&E        Chn - neutral transfer coefficient for heat
!&E        Cen - neutral transfer coefficient for moisture
!&E        RR - Roughness Reynolds number
!&E        RT - Roughness Reynolds number for temperature
!&E        RQ - Roughness Reynolds number for moisture
!&E        zL - stability parameter height/L where L is the Obukhov length
!&E        zo - velocity roughness length
!&E        zot - temperature roughness length
!&E        zoq - humidity roughness length
!&E        dt_wrm - total warm layer temperature difference C
!&E        tk_pwp - thickness of warm layer
!&E        dsea - (=dt_wrm*ts_depth/tk_pwp) warming above sea temperature measurement
!&E        dter - cool skin temperature difference C
!&E        tkt - cool skin thickness mm (printed out *1000)
!&E        sst - skin temperature C (sst = ts - dter + dsea)
!&E        Wg - gustiness factor m/s
!&E        Wbar - Webb mean vertical velocity m/s
!&E        ef_webb - latent flux due to the Webb correction
!&E        u_zs - wind velocity at standard height
!&E        t_zs - air temperature ditto
!&E        q_zs,qa_zs,e_zs,rh_zs - humidity ditto in a variety of units
!&E
!&E ** Called by : sflxsurf.F90
!&E
!&E ** External calls :
!&E
!&E ** Used ij-arrays :
!&E
!&E ** Modified variables : zsenslos, zlatlos, zevap, zchs, zces
!&E
!&E ** Reference : COARE Bulk Flux Algorithm version 3.0a. See Fairall et al.,2003. J.Climate,16,571-591
!&E                METEOCEAN report (2009)
!&E
!&E ** History :
!&E       !  2003/09  (Frank Bradley) Original code
!&E       !  2008-08  (L.-H. Ayina)
!&E       !  2009-06  (V. Garnier) MARS implementation
!&E
!&E----------------------------------------------------------------------
!! * Modules used
   USE parameters, ONLY : grav
   USE comvars2d, ONLY  : jb,jh,ig,id,nbouc,t,iwarnlog
   USE comfluxdo, ONLY  : rain

!! * Arguments
   REAL(KIND=rsh), INTENT( IN )      :: ztemp,zsat,zwindmod,zpa,zrh,zrain, &
                         zirlos, &   ! downward IR flux [W/m**2]
                         zrad        ! downward solar flux [W/m**2]
   REAL(KIND=rsh), INTENT( OUT )     :: zsenslos,zlatlos,zevap
   REAL(KIND=rsh), INTENT( INOUT )   :: ztosx,ztosy

!! * Local declarations
!   REAL(KIND=rsh), DIMENSION(imin:imax,ljmin:ljmax) ::   &
!         xlat



   INTEGER, PARAMETER :: jwave=0,       &   ! implement wave dependent roughness model
! Charnock formulation
                         jcool=0            ! implement cool calculation skin (0=no, 1=yes)
! SST measured by IR radiometer


!! * Fairall parameters
!! --------------------

   REAL(KIND=rsh), PARAMETER ::        &
                        zu=10.0_rsh,          &    ! wind measurement height             [m]
                        zq=2.0_rsh,           &    ! humidity (q) measurement height     [m]
!                        ts_depth=1.0_rsh,     &
                        Cpa=1004.67_rsh,      &    ! specific heat of air at constant pressure [J/kg/K]
                        zs=10.0_rsh,          &    ! standard height                     [m]
                        p=1008.0_rsh,         &    ! pressure if measurement not available [mbar]
                        zi=600.0_rsh,         &    ! boundary layer depth [m]
                        Rgas=287.1_rsh,       &
                        cpw=4000.0_rsh,       &
                        rhow=1022.0_rsh,      &
                        von=0.4_rsh,          &
                        Beta=1.2_rsh,         &
                        be=0.026_rsh,         &
                        twopi=3.141593_rsh*2.0_rsh, &    ! 2*pi
                        visw=1.e-6,           &
                        tcw=0.6_rsh

   REAL(KIND=rsh) :: expo,expa,zprs,hwave,twave,ztsw,dter,zsst,ea,es,qa, &
                     Cpv, rhoa,Al,ztka,ztko,qs,                          &
                     bigc,lnut,wetc,visa,cwave,                          &
                     lwave,dqer,zo,Wg,tkt,Du,Dt,Dq,u10,                  &
                     zo10,Cd10,Ch10,Ct10,zot10,Ct,Cd,CC,Ribcu,           &
                     Ribu,zetu,L1,psu10,pst10,psq10

   REAL(KIND=rsh) :: usr,tsr,qsr,charn,rr,zoq,zot, zL,L,psu,pst,         &
                     Bf,Rnl,hsb,hlb,qout,dels,qcool,alq, xlamx,          &
                     dwat,dtmp,dqs_dt,alfac,rf,taur,Wbar,ef_webb,        &
                     CDS,CHS,CES,Cdn,Chn,Cen,                            &
                     u_zs,t_zs,q_zs,qa_zs,e_zs,esat_zs,rh_zs,            &
                     zwind,L10,tau,ef,hf

   INTEGER :: nits,iter                              


!!----------------------------------------------------------------------
!! * Executable part


! --------------
! Initialization
! --------------


! Relative wind velocity [m/s] (module)
! prevoir u_wind-u_ocean selon la direction du vent
! -------------------------------------

         zwind = zwindmod


! default wave parameters
! -----------------------

!  ztso=hwt !un-comment to use 6m depth sea temperature; also change ts_depth above
!  default values for wave height and period for equilibrium sea
!  if known, add to data and to read statement above,then comment next two lines

         hwave=0.018_rsh*zwind*zwind
         hwave=hwave*(1.0_rsh+.015_rsh*zwind)
         twave=0.729_rsh*zwind


! Sea surface temperature [Celsius]
! ---------------------------------
 
         ztsw=ztemp
         dter=0.3_rsh*jcool
         zsst=ztemp-dter

! jump over warm layer calculation
! convert humidities and calculate fluxes in ASL
         IF (jwave == 0) ztsw=ztemp


! air and sea temperature [K]
!----------------------------

         ztka   = zsat + kelvin                 ! air temperature [K]
         ztko   = ztemp + kelvin                ! sea temperature [K]


! convert humidities and calculate fluxes in ASL
! ----------------------------------------------

! Atmospherical pressure [Pa ---> mbar]
         zprs = zpa/100.0_rsh

! The same formulation must be used for the air and the ocean
         CALL sh_fairall(zsat,zrh,zprs,ea,qa)      ! Teten s returns sat. air qa in kg/kg
         CALL sh_fairall(ztemp,suswe,zprs,es,qs)   ! Teten s returns sat. sea qs in kg/kg


!         CALL gravity(xlat(i,j),grav)                     ! gravity


! ------------------
! Fairall parameters
! ------------------


! Constants and coefficients (Stull 1988 p640).
!----------------------------------------------

         Cpv=Cpa*(1.0_rsh+0.00084_rsh*qa)                        ! qa (kg/kg)
         rhoa=zprs*100.0_rsh/(Rgas*ztka*(1.0_rsh+0.61_rsh*qa))   ! qa (kg/kg)  and zprs (mb)  (bug in core3.0)
         Al=2.1e-5*(ztemp+3.2_rsh)**0.79_rsh


! Factors :Beta=1.2            !Given as 1.25 in Fairall et al.(1996)
!          twopi=3.14159*2.
!
! Additional constants needed for cool skin
! visw=1.e-6                   !m2/s kinematic viscosity water
! tcw=0.6                      !W/m/K   Thermal conductivity water

         bigc=16.0_rsh*grav*cpw*(rhow*visw)**3
         bigc=bigc/(tcw*tcw*rhoa*rhoa)
         lnut=(2.501_rsh-0.00237_rsh*ztsw)*1e+6         ! latent heat of vaporization [J/kg]
! at ts (3C warming=0.3%)
         wetc=0.622_rsh*lnut*qs/(Rgas*(ztsw+kelvin)**2) ! Clausius-Clapeyron

         visa= 4.84e-9*zsat*zsat*zsat
         visa= 8.301e-6*zsat*zsat-visa
         visa=1.326e-5*(1.0_rsh+6.542e-3*zsat+visa)



! Wave parameters of Oost et al (2002)
! ------------------------------------

         cwave=grav*twave/twopi
         lwave=cwave*twave



! -----------
! First guess
! -----------

         dter=0.3_rsh*jcool              ! cool skin Dt
         dqer=wetc*dter                  ! cool skin Dq
         zo=0.0001_rsh
         Wg=0.5_rsh                      ! Gustiness factor initial guess
         tkt= 0.001_rsh*jcool            ! Cool skin thickness first guess


! Air-sea differences - includes warm layer in Dt and Dq
! ------------------------------------------------------

         Du=(zwind**2+Wg**2)**0.5_rsh    ! include gustiness in wind spd. difference
         Dt=ztsw-zsat-0.0098_rsh*zt      ! potential temperature difference.
         Dq=qs-qa


! neutral coefficients
! --------------------

         u10=Du*LOG(10/zo)/LOG(zu/zo)
         usr=0.035_rsh*u10
         zo10=0.011_rsh*usr*usr/grav+0.11_rsh*visa/usr
         Cd10=(von/LOG(10/zo10))**2
         Ch10=0.00115_rsh
         Ct10=Ch10/sqrt(Cd10)
         zot10=10.0_rsh/EXP(von/Ct10)
         Cd=(von/LOG(zu/zo10))**2


! Grachev and Fairall (JAM, 1997)
! -------------------------------

         Ct=von/LOG(zt/zot10)                    ! Temperature transfer coefficient
         CC=von*Ct/Cd                            ! z/L vs Rib linear coefficient
         Ribcu=-zu/(zi*0.004_rsh*Beta**3)        ! Saturation or plateau Rib
         Ribu=-grav*zu*((Dt-dter)+0.61_rsh*ztka*Dq)/(ztka*Du**2)
         IF (Ribu < 0.0_rsh) THEN
           zetu=CC*Ribu/(1.0_rsh+Ribu/Ribcu)               ! Unstable G and F
         ELSE
           zetu=CC*Ribu*(1.0_rsh+27.0_rsh/9.0_rsh*Ribu/CC) ! Stable
         END IF
         IF (zetu == 0.0_rsh) THEN
           zetu = -0.1_rsh ! cval WARNING AIE WHICH VALUE IS CORRECT ?

           WRITE(iwarnlog,*) "FAIRALL at t: ", t,"zetu == 0.0_rsh and set to -0.1_rsh"

         END IF
         L10=zu/zetu                          ! MO length
         IF (zetu > 50.0_rsh) THEN
           nits=1
         ELSE
           nits=3                       ! number of iterations
         END IF


! First guess M-O stability dependent scaling params.(u*,t*,q*) to estimate zo and z/L
! ------------------------------------------------------------------------------------

         psu10=psiu(zu/L10)
         pst10=psit(zt/L10)
         psq10=psit(zq/L10)

         usr= Du*von/(LOG(zu/zo10)-psu10)
         tsr=-(Dt-dter)*von/(LOG(zt/zot10)-pst10)
         qsr=-(Dq-dqer)*von/(LOG(zq/zot10)-psq10)


! Charnock for high wind speeds Chris  data
! -----------------------------------------

         charn=0.011_rsh     !then modify Charnock for high wind speeds Chris  data
         IF (Du > 10.0_rsh) charn=0.011_rsh+(0.018_rsh-0.011_rsh)*(Du-10.0_rsh)/(18.0_rsh-10.0_rsh)
         IF (Du > 18.0_rsh) charn=0.018_rsh



!----------------------------------------------------------------------------
! **** Iterate across u*(t*,q*),zo(zot,zoq) and z/L including cool skin ****
!----------------------------------------------------------------------------

         DO iter=1,nits


! zo estimate
! -----------

           IF (jwave == 0) THEN
             zo=charn*usr*usr/grav + 0.11_rsh*visa/usr          ! after Smith 1988
           ELSE IF (jwave == 1) THEN
             zo=(50.0_rsh/twopi)*lwave*(usr/cwave)**4.5_rsh+0.11_rsh*visa/usr ! Oost et al.
           ELSE IF (jwave == 2) THEN
             zo=1200.0_rsh*hwave*(hwave/lwave)**4.5_rsh+0.11_rsh*visa/usr     ! Taylor and Yelland
           END IF
           rr=zo*usr/visa


! zoq and zot fitted to results from several ETL cruises
! ------------------------------------------------------

           zoq=min(1.15e-4,5.5e-5/rr**0.6_rsh)
           zot=zoq

!
! -------------------------------------------------------------------------
           zL=von*grav*zu*(tsr*(1.+0.61_rsh*qa)+0.61_rsh*ztka*qsr)           &
              /(ztka*usr*usr*(1.+0.61_rsh*qa))
           L=zu/zL
           psu=psiu(zu/L)
           pst=psit(zt/L)
           dqer=wetc*dter*jcool
           usr=Du*von/(LOG(zu/zo)-psiu(zu/L))
           tsr=-(Dt-dter)*von/(LOG(zt/zot)-psit(zt/L))
           qsr=-(Dq-dqer)*von/(LOG(zq/zoq)-psit(zq/L))
           Bf=-grav/ztka*usr*(tsr+0.61_rsh*ztka*qsr)

           IF (Bf > 0.0_rsh) THEN
             Wg=Beta*(Bf*zi)**0.333_rsh
           ELSE
             Wg=0.2_rsh
           END IF

           Du=sqrt(zwind**2+Wg**2)                     ! include gustiness in wind spd.


! Cool skin
! ---------

           IF (jcool == 1) THEN

! background sensible and latent heat
             hsb=-rhoa*cpa*usr*tsr
             hlb=-rhoa*lnut*usr*qsr

! Recompute net longwave; cool skin=-2W/m2
!             Rnl= 0.97_rsh*(5.67e-8*(ztsw-dter+kelvin)**4-zrlo)
!             qout=Rnl+hsb+hlb

! thermal + sensible + latent heat fluxes
             qout=zirlos+hsb+hlb         ! ATTENTION signe et unite de irlos

! solar flux estimate
!             dels=(0.065_rsh+11.0_rsh*tkt-6.6e-5/tkt*(1.0_rsh-EXP(-tkt/8.0e-4))) !Eq.16 Ohlmann
!             dels=zrso*0.945_rsh*dels  ! the ocean sees only 94.5% of the incident solar  (penetrative flux)
!             qcool=qout-dels

! ATTENTION calcul de qout et qcool n est pas clair - signe des flux ?

! cooling flux (the penetrative solar flux -which warms the sea surface- is removed)
! apres etude de sensibilite, utiliser zrad * exp(-coext*dz) (on veut utiliser le flux penetrant)
             qcool=qout-zrad            ! ATTENTION signe et unite
             alq=Al*qcool+be*hlb*cpw/lnut                         ! Eq. 7 Buoy flux water

             IF (alq > 0.0_rsh) THEN                              ! originally (qcool.gt.0)
               xlamx=6.0_rsh/(1.0_rsh+(bigc*alq/usr**4)**0.75_rsh)**0.333_rsh     ! Eq 13 Saunders coeff.
               tkt=xlamx*visw/(sqrt(rhoa/rhow)*usr)               ! Eq.11 Sublayer thickness
             ELSE
               xlamx=6.0_rsh                                      ! prevent excessive warm skins
               tkt=MIN(0.01_rsh,xlamx*visw/(sqrt(rhoa/rhow)*usr)) ! Limit tkt
             END IF

             dter=qcool*tkt/tcw                                   ! Eq.12 Cool skin
             dqer=wetc*dter

           END IF    ! ends skin cool

         END DO                                             ! end iterations

!---------------------------
! **** End of iteration ****
!---------------------------


!--------------------------------------------
! compute surface fluxes and other parameters
!--------------------------------------------

         zsst=ztsw-dter*jcool          ! final skin temperature at this timestep
! warming at the sea surface relative to the heat lost (net flux)
! between the mid ocean surface layer and the sea surface
! dT/dz = Q

! relations 1.3.2 in METEOCEAN 2009
         tau=rhoa*usr*usr*zwind/Du     ! stress N/m2
         hf=-cpa*rhoa*usr*tsr          ! sensible W/m2
         ef=-lnut*rhoa*usr*qsr         ! latent W/m2


!----------------------------------
! compute heat flux due to rainfall
!----------------------------------

         dwat=2.11e-5*((ztka)/kelvin)**1.94_rsh                  ! water vapour diffusivity
         dtmp=(1.0_rsh+3.309e-3*zsat-1.44e-6*zsat*zsat)*0.02411_rsh
         dtmp=dtmp/(rhoa*cpa)                                    ! heat diffusivity
         dqs_dt=qa*lnut/(Rgas*(ztka)**2)                         ! Clausius-Clapeyron
         alfac= 1.0_rsh/(1.0_rsh+(wetc*lnut*dwat)/(cpa*dtmp))    ! wet bulb factor
!rf= zrain*alfac*cpw*((zsst-zsat)+(qs-qa-dqer)*lnut/cpa)/3600.0_rsh    ! ATTENTION signe et unite (/3600 car passe en mm/h dans sflxread originalement)
         rf= zrain*alfac*cpw*((zsst-zsat)+(qs-qa-dqer)*lnut/cpa)               ! ATTENTION signe et unite (stupide, on reste en mm/s dans sflxread)
         rf = 0.0_rsh                                            ! the rain is not taken into account at the moment


!--------------------------------------
! compute momentum flux due to rainfall
!--------------------------------------

!taur=0.85_rsh*zrain/3600.0_rsh*zwind                    ! ATTENTION signe et unite de la pluie
         taur=0.85_rsh*zrain*zwind                    ! ATTENTION signe et unite de la pluie (on reste en mm/s dans sflxread)
         taur=0.0_rsh                                            ! the rain is not taken into account at the moment


!-------------------------------------------------------------------------------------
! Webb correction to latent heat flux already in ef via zoq/rr function so return Wbar
!-------------------------------------------------------------------------------------

!      Wbar=-1.61_rsh*usr*qsr/(1.0_rsh+1.61_rsh*qa)-usr*tsr/(ztka)
         expo=-1.61_rsh*usr*qsr/(1.0_rsh+1.61_rsh*qa)
         expa=usr*tsr
         expa=expa/ztka
         Wbar= expo-expa

! latent heat flux due to the Webb correction
         ef_webb=lnut*rhoa*Wbar*qa
         ef_webb=0.0_rsh


!---------------------------------------------------------------------------
! compute neutral transfer coefficients and met variables at standard height
!---------------------------------------------------------------------------

         Cdn=(von/LOG(zs/zo))**2
         Chn=von*von/(LOG(zs/zo)*LOG(zs/zot))
         Cen=von*von/(LOG(zs/zo)*LOG(zs/zoq))


!------------------------------
! compute transfer coefficients
!------------------------------

         CDS=(usr/Du)**2
         IF (zsat/=zsst) THEN
           CHS=usr*tsr/(Du*(zsat-zsst+0.0098_rsh*zt))
         ELSE
           CHS=Chn
         END IF
         IF (qa/=qs) THEN
           CES=usr*qsr/(Du*(qa-qs+dqer))
         ELSE
           CES=Cen
         END IF



!---------------------------------------------------
! adjust meteorological variables to standard height
! (aboard, measurements are done at different atitudes)
!---------------------------------------------------




!--------------
! output fluxes
!--------------

! Wind stress (N/m2)

         ztosx=rhoa*CDS*ztosx
         ztosy=rhoa*CDS*ztosy

! Du is the modified module of the wind
! from which all the exchange coefficients have been estimated.
! ztosx = ||V||*Vx ; ztosy = ||V||*Vy, so Du/zwindmod removes
! the module of the wind in the formula,
! the second factor Du/zwindmod corrects the component of the wind

         IF (zwindmod > 0.0_rsh) THEN
           ztosx=ztosx*(Du/zwindmod)**2
           ztosy=ztosy*(Du/zwindmod)**2
         ENDIF



! Turbulent fluxes

         zsenslos = -(hf+rf)
         zlatlos = -(ef+ef_webb)
         zevap = -zlatlos/lnut



  END SUBROUTINE sflx_turb_fairall


!!======================================================================

  SUBROUTINE sh_fairall(ssti,rhim,pris,pvis,shim)

!&E---------------------------------------------------------------------
!&E                       ***  ROUTINE sh_fairall ***
!&E
!&E       ** Purpose : calculate specific humidity and saturated vapour pressure
!&E
!&E       ** Description : - determined as function of relative humidity
!&E                        - alternative using dew point temperature instead
!&E                          of relative humidity is possible but not included
!&E
!&E       ** Called by :
!&E
!&E       ** External calls :
!&E
!&E       ** Reference : - Teten s formula for saturation vo Buck (1981)
!&E                                                   JAM 20, 1527-1532
!&E       ** History :
!&E             !
!&E              press = 1013.25 ! standard atmospheric pressure [mb]
!&E---------------------------------------------------------------------
   REAL(KIND=rsh), INTENT( IN )  :: ssti,rhim,pris
   REAL(KIND=rsh), INTENT( OUT ) :: shim
   REAL(KIND=rsh), INTENT( OUT)  :: pvis
   REAL(KIND=rsh), PARAMETER     ::         &
                                    cons1=1.0007_rsh,    &
                                    cons2=3.46e-6,       &
                                    cons3=6.112_rsh,     &
                                    cons4=17.502_rsh,    &
                                    cons5=240.97_rsh,    &
                                    cons6=0.622_rsh,     &
                                    cons7=1.0_rsh-cons6

   REAL(KIND=rsh)                :: expo,expo1,expa

   expa = cons4*ssti/(cons5+ssti)
   expo = exp(expa)
   expo1 = (cons1+cons2*pris)*cons3*expo
!   IF (rhim < 2.0_rsh) THEN
     pvis = expo1*rhim
     shim = (cons6*pvis)/(pris-cons7*pvis)
!   ELSE
!     shim = rhim/1000.0_rsh                 !g/kg to kg/kg
!   END IF

  END SUBROUTINE sh_fairall


!!======================================================================

  SUBROUTINE humidity(T,P,esat)

!&E------------------------------------------------------------------
!&E      Tetens  formula for saturation vp Buck(1981) JAM 20, 1527-1532
!&E------------------------------------------------------------------
   REAL(KIND=rsh), INTENT( IN )  ::  T,P
   REAL(KIND=rsh), INTENT( OUT ) ::  esat

   esat = (1.0007_rsh+3.46e-6*P)*6.1121_rsh*EXP(17.502_rsh*T/(240.97_rsh+T)) ! mb

  END SUBROUTINE humidity


!!======================================================================

  FUNCTION psiu(zL)
!&E---------------------------------------------------------------------
!&Epsiu and psit evaluate stability function for wind speed and scalars
!&E  matching Kansas and free convection forms with weighting f
!&E     convective form follows Fairall et al (1996) with profile constants
!&E       from Grachev et al (2000) BLM
!&E       stable form from Beljaars and Holtslag (1991)
!&E--------------------------------------------------------------------
   REAL(KIND=rsh), INTENT( IN ) :: zL
   REAL(KIND=rsh)               :: x,y,psik,psic,f,psiu,c

   IF (zL <= 0.0_rsh) THEN
     x=(1.0_rsh-15.0_rsh*zL)**0.25_rsh              ! Kansas unstable
     psik=2.0_rsh*LOG((1.0_rsh+x)/2.0_rsh)+LOG((1.0_rsh+x*x)/2.0_rsh)  &
          - 2.0_rsh*ATAN(x)+2.0_rsh*ATAN(1.0_rsh)
     y=(1.0_rsh-10.15_rsh*zL)**0.3333_rsh                    ! Convective
     psic=1.5_rsh*LOG((1.0_rsh+y+y*y)/3.0_rsh)                       &
            - SQRT(3.0_rsh)*ATAN((1.0_rsh+2.0_rsh*y)/SQRT(3.0_rsh))  &
            + 4.0_rsh*ATAN(1.0_rsh)/SQRT(3.0_rsh)
     f=zL*zL/(1.+zL*zL)
     psiu=(1.0_rsh-f)*psik+f*psic
   ELSE
     c=min(50.0_rsh,0.35_rsh*zL)                       ! Stable
!cval       psiu=-((1.0_rsh+1.0_rsh*zL)**1.+.6667*(zL-14.28)/EXP(c)+8.525)
     psiu=-((1.0_rsh+1.0_rsh*zL)+0.6667_rsh*(zL-14.28_rsh)/EXP(c)+8.525_rsh)
   END IF

  END FUNCTION psiu


!!======================================================================

  FUNCTION psit(zL)
!&E---------------------------------------------------------------------
!&E       psiu and psit evaluate stability function for wind speed and scalars
!&E       matching Kansas and free convection forms with weighting f
!&E       convective form follows Fairall et al (1996) with profile constants
!&E       from Grachev et al (2000) BLM
!&E       stable form from Beljaars and Holtslag (1991)
!&E--------------------------------------------------------------------

   REAL(KIND=rsh), INTENT( IN )  :: zL
   REAL(KIND=rsh)                :: x,y,psik,psic,f,psit,c

   IF (zL <= 0.0_rsh) THEN
     x=(1.0_rsh-15.0_rsh*zL)**0.5_rsh                        ! Kansas unstable
     psik=2.0_rsh*LOG((1.0_rsh+x)/2.0_rsh)
     y=(1.0_rsh-34.15_rsh*zL)**0.3333_rsh                    ! Convective
     psic=1.5_rsh*LOG((1.0_rsh+y+y*y)/3.0_rsh)                       &
            - SQRT(3.0_rsh)*ATAN((1.0_rsh+2.0_rsh*y)/SQRT(3.0_rsh))  &
            + 4.0_rsh*ATAN(1.0_rsh)/SQRT(3.0_rsh)
     f=zL*zL/(1.0_rsh+zL*zL)
     psit=(1.0_rsh-f)*psik+f*psic
   ELSE
     c=min(50.0_rsh,0.35_rsh*zL)                        !Stable
     psit=-((1.0_rsh+2.0_rsh*zL/3.0_rsh)**1.5_rsh+0.6667_rsh*(zL-14.28_rsh)/EXP(c)+8.525_rsh)
   END IF

  END FUNCTION psit

!!======================================================================

  SUBROUTINE gravity(zzlat,zgra)
!&E--------------------------------------------------------------------------
!&E      calculates g as a function of latitude using the 1980 IUGG formula
!&E
!&E      Bulletin Geodesique, Vol 62, No 3, 1988 (Geodesist s Handbook)
!&E      p 356, 1980 Gravity Formula (IUGG, H. Moritz)
!&E      units are in m/sec^2 and have a relative precision of 1 part
!&E      in 10^10 (0.1 microGal)
!&E      code by M. Zumberge.
!&E
!&E      check values are:
!&E
!&E      g = 9.780326772 at latitude  0.0
!&E      g = 9.806199203 at latitude 45.0
!&E      g = 9.832186368 at latitude 90.0
!&E--------------------------------------------------------------------------
   USE parameters, ONLY : rlg
   REAL(KIND=rsh), INTENT( IN )  :: zzlat
   REAL(KIND=rsh), INTENT( OUT ) :: zgra
   REAL(KIND=rlg), PARAMETER     :: ppi=3.14159265358979_rlg
   REAL(KIND=rlg)                :: gamma, c1, c2, c3, c4, phi

   phi = REAL(zzlat,rlg) * ppi/180.0_rlg
   gamma = 9.7803267715_rlg
   c1 = 0.0052790414_rlg
   c2 = 0.0000232718_rlg
   c3 = 0.0000001262_rlg
   c4 = 0.0000000007_rlg
   zgra = REAL ( gamma * (1.0_rlg    &
          + c1 * ((sin(phi))**2)     &
          + c2 * ((sin(phi))**4)     &
          + c3 * ((sin(phi))**6)     &
          + c4 * ((sin(phi))**8)), rsh)

  END SUBROUTINE gravity


!!======================================================================

  SUBROUTINE sflx_turb_luyten(ztemp,zsat,zwindmod,zpa,zrh,zrad,zirlos,zrain, &
                              zsenslos,zlatlos,zevap,ztosx,ztosy)
 
!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE sflx_turb_luyten  ***
!&E
!&E ** Purpose : latent and sensible heat fluxes estimate
!&E
!&E ** Description :
!&E
!&E      zwindmod      ! module of wind velocity at 10 m (m/s)
!&E      zrh           ! relative humidity       (0 < rh  < 1)
!&E      atmp          ! atmospherical pressure in mb (but read in Pascal)
!&E      at            ! air temperature (Celsius)
!&E      st            ! sea surface temperature (Celsius)
!&E      roair         ! air density (kg/m3)
!&E      cp            ! specific heat of air at constant pressure (1004 J/kg/K)
!&E      cevap         ! latent heat flux coefficient
!&E      cheat         ! sensible heat flux coefficient
!&E      e             ! vapour pressure (in mb)
!&E                      e = rh*10.0**(9.4_rsh - 2353.0_rsh/temp
!&E      sh            ! specific humidity
!&E                      sh=0.622_rsh * e / (atmp-0.378_rsh*e)
!&E
!&E
!&E ** Called by : sflxsurf.F90
!&E
!&E ** External calls :
!&E
!&E ** Reference : Luyten & De Mulder (1992)
!&E                METEOCEAN report (2009)
!&E
!&E ** History :
!&E       !           (P. Lazure)  Original code
!&E       !  2006-05  (V. Garnier) 0 < rh, fcc < 1
!&E       !  2009-05  (V. Garnier) correction according to METEOCEAN report
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE sflxatm,   ONLY   : drag, exch, imeteo_exchtype, imeteo_dragtype

!! * Declaration subroutine

!! * Arguments
   REAL(KIND=rsh), INTENT( IN )  :: ztemp,zsat,zwindmod,zpa,zrh, &
                                    zrad,zirlos,zrain
   REAL(KIND=rsh), INTENT( OUT ) :: zsenslos,zlatlos,zevap
   REAL(KIND=rsh), INTENT( INOUT ) :: ztosx,ztosy

!! * Local declarations
   REAL(KIND=rsh), PARAMETER   :: roair=1.25
   REAL(KIND=rsh)              :: st,at,e_a,e_w,sh_a,sh_w,cevap,var,tdif,u,cheat,atmp
   REAL(KIND=rsh)              :: zdsenslos,zdlatlos

!!----------------------------------------------------------------------
!! * Executable part


! sea temperature (degree C)
         st=ztemp
! air temperature (degree C)
         at=zsat
! difference of potential temperature at the air-sea interface
! the air temperature increases by 0.98 degrees Celsius per 100 meters in adiabatic conditions
         tdif=at + 0.0098_rsh*zt -st
         IF(tdif> 5.0_rsh) tdif=5.0_rsh
         IF(tdif<-5.0_rsh) tdif=-5.0_rsh
! module of wind velocity
         u=zwindmod
! atmospherical pressure (Pascal to mb)
         atmp = zpa/100.0_rsh

! water
! latent heat flux coefficient
         cevap=ce(u,tdif,exch,imeteo_exchtype)
! vapour pressure (in mb)
         e_w=suswe*10**(9.4_rsh - 2353.0_rsh/(st+kelvin))
! specific humidity
         sh_w=0.622_rsh * e_w / (atmp-0.378_rsh*e_w)

! air
! vapour pressure (in mb)
         var=9.4_rsh-2353.0_rsh/(at+kelvin)
         e_a=zrh * 10**var
! specific humidity
         sh_a=0.622_rsh * e_a / (atmp-0.378_rsh*e_a)

! latent flux
         zlatlos=-roair*lnu*cevap*(sh_w-sh_a)*u

! evaporation
         IF (l_sflx_rrain) zevap=-zlatlos / lnu

! sensible heat flux coefficient
         cheat=cevap
! sensible flux
         zsenslos=cheat*rhoair*cp*u*tdif

! wind stress
         ztosx=rhoair*cd(zwindmod,tdif,drag,imeteo_dragtype,imeteo_exchtype)*ztosx
         ztosy=rhoair*cd(zwindmod,tdif,drag,imeteo_dragtype,imeteo_exchtype)*ztosy



! sea surface flux relaxation
         IF (l_rapflux) THEN
           zdsenslos=-cheat*rhoair*cp*u
           PRINT*,'sflxturb : dlatlos not ready yet :'
           PRINT*,'the simulation is stopped'
           STOP
!           zdlatlos=-roair*lnu*cevap*u*dq/dt
         END IF


  END SUBROUTINE sflx_turb_luyten

!!======================================================================

  SUBROUTINE sflx_turb_default(ztemp,zsat,zwindmod,zpa,zrh,zrad,zirlos,zrain, &
                               zsenslos,zlatlos,zevap,ztosx,ztosy)
 
!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE sflx_turb_default ***
!&E
!&E ** Purpose : latent and sensible heat fluxes estimate
!&E
!&E ** Description :
!&E
!&E      zmodwind      ! module of wind velocity at 10 m (m/s)
!&E      zrh           ! relative humidity       (0 < rh  < 1)
!&E      atmp          ! atmospherical pressure in mb (but read in Pascal)
!&E      at            ! air temperature (Celsius)
!&E      st            ! sea surface temperature (Celsius)
!&E      roair         ! air density (kg/m3)
!&E      cp            ! specific heat of air at constant pressure (1004 J/kg/K)
!&E      cevap         ! latent heat flux coefficient
!&E                    ! cevap=0.0012_rsh*(1.0_rsh+1/u)
!&E      cheat         ! sensible heat flux coefficient
!&E                    ! 1.45e-3
!&E      e             ! vapour pressure (in mb)
!&E                      e = exp(2.302*(7.5*temp/(237.3+temp)+0.7858))
!&E      sh            ! specific humidity
!&E                      sh=0.622_rsh * e / (atmp-0.378_rsh*e)
!&E      lnt           ! latent heat of vaporization depends on the sea surface temperature
!&E
!&E ** Called by : sflxsurf.F90
!&E
!&E ** External calls :
!&E
!&E ** Reference : Clark et al (1995) latent flux + private modification for roair and lnu
!&E                Elliot et Clark (1990)  sensible flux + private modification for roair
!&E                METEOCEAN report (2009)
!&E
!&E ** History :
!&E       !           (P. Lazure)  Original code
!&E       !  2006-05  (V. Garnier) 0 < rh, fcc < 1
!&E       !  2009-05  (V. Garnier) correction according to METEOCEAN report
!&E
!&E---------------------------------------------------------------------
!! * Modules used
   USE sflxatm,   ONLY   : drag, imeteo_exchtype, imeteo_dragtype

!! * Declaration subroutine

!! * Arguments
   REAL(KIND=rsh), INTENT( IN ) :: ztemp,zsat,zwindmod,zpa,zrh, &
                                   zrad,zirlos,zrain
   REAL(KIND=rsh), INTENT( OUT ) :: zsenslos,zlatlos,zevap
   REAL(KIND=rsh), INTENT( INOUT ) :: ztosx,ztosy

!! * Local declarations
   REAL(KIND=rsh), PARAMETER   :: cheat = 1.45e-3
   REAL(KIND=rsh)              :: e_a,e_w,sh_a,sh_w,lnt,st,at,cevap,var,tdif,u,roair,atmp
   REAL(KIND=rsh)              :: zdsenslos,zdlatlos

!!----------------------------------------------------------------------
!! * Executable part

! sea temperature (degree C)
         st=ztemp
! air temperature (degree C)
         at=zsat
! square of wind velocity
         u=MAX(zwindmod,1.0_rsh)
! atmospherical pressure (Pascal to mb)
         atmp = zpa/100.0_rsh

! water
! air density is function of the air temperature
         roair=1.293_rsh*kelvin/(at+kelvin)
! latent heat flux coefficient
         cevap=0.0012_rsh*(1.0_rsh+1.0_rsh/u)
! vapour pressure (in mb)
         e_w=suswe*EXP(2.302_rsh*(7.5_rsh*st/(237.3_rsh+st)+0.7858_rsh))
! specific humidity
         sh_w=0.622_rsh*e_w/(atmp-0.378_rsh*e_w)

! air
! vapour pressure (in mb)
         e_a=zrh*EXP(2.302_rsh*(7.5_rsh*at/(237.3_rsh+at)+0.7858_rsh))
! specific humidity
         sh_a=0.622_rsh*e_a/(atmp-0.378_rsh*e_a)
! latent heat of vaporization (J/kg) from Gill 1982 (t in degrees Celsius)
         lnt=(2500.8_rsh-2.3_rsh*st)*1000.0_rsh

! latent flux
         zlatlos=-lnt*roair*cevap*(sh_w-sh_a)*u

! evaporation
         IF (l_sflx_rrain) zevap=-zlatlos / lnt

! sensible flux
         zsenslos=-cheat*rhoair*cp*u*(st-at-0.0098_rsh*zt)
! difference of potential temperatures
! the air temperature increases by 0.98 degrees Celsius per 100 meters in adiabatic conditions

! wind stress
         ztosx=rhoair*cd(zwindmod,tdif,drag,imeteo_dragtype,imeteo_exchtype)*ztosx
         ztosy=rhoair*cd(zwindmod,tdif,drag,imeteo_dragtype,imeteo_exchtype)*ztosy



! sea surface flux relaxation
         IF (l_rapflux) THEN
           zdsenslos=-cheat*rhoair*cp*u
           zdlatlos=-lnt*roair*cevap*sh_w
         END IF

  END SUBROUTINE sflx_turb_default

!!======================================================================
    
  FUNCTION cd(wind,tdif,drag,imeteo_dragtype,imeteo_exchtype)
    
!§E---------------------------------------------------------------------
!§E                 ***  FUNCTION cd ***
!§E
!§E ** Purpose :  calculate surface drag coefficient
!§E
!§E ** Description :
!§E     - options selected by imeteo_dragtype (wind dependence) and imeteo_exchtype
!§E               (dependence on air-sea temperature difference)
!§E       imeteo_dragtype = 0 -> constant value (if imeteo_exchtype=0)
!§E             = 1 -> large and pond (1981)
!§E             = 2 -> smith and banke (1975)
!§E             = 3 -> geernaert et al. (1986)
!§E             = 4 -> charnock s relation
!§E       imeteo_exchtype = 0 -> no dependence on tdif
!§E             = 1 -> as function of tdif
!§E
!§E     - uses bilinear interpolation if imeteo_exchtype = 1
!§E             (dependence on air-sea temperature difference and wind)
!§E
!§E ** Called by : sflx_surf
!§E
!§E ** External calls : charno
!§E
!§E ** Reference :
!§E             - section iii-1.6.4 of the user documentation
!§E             - geernaert g.l., 1990. bulk parameterizations
!§E               for the wind stress and heat fluxes.
!§E               in: g.l. geernaert and w.j. plant (editors),
!§E               surface waves and fluxes, vol. 1 - current theory.
!§E               kluwer academic publishers, dordrecht, pp. 91-172.
!§E
!§E ** History :
!§E       !         (P. Luyten)  Original code
!§E       ! 04-12   (G. Andre)   Introduction dans MARS
!§E       ! 06-04   (V. Garnier) Fortran 90
!§E
!§E---------------------------------------------------------------------
 
!! * Declaration function
   REAL(kind=rsh)             :: cd

!! * Arguments
   INTEGER,       INTENT(in)  :: imeteo_exchtype, &   ! switch to select dependence
! on tdif
                                 imeteo_dragtype      ! type of drag coefficient
! formulation
   REAL(kind=rsh),INTENT(in)  :: tdif,  &   ! air-sea temperature difference
! (used only if imeteo_exchtype=1)
                                 wind       ! wind speed at 10 m height [m/s]
   REAL(kind=rsh),INTENT(in),DIMENSION(13,11) :: &
                                 drag       ! array used for bilinear evaluation
! of drag coefficient

!! * Local declarations
   INTEGER                    :: irow, jcol
   REAL(kind=rsh) :: f11, f12, f21, f22, z1, z2

!!----------------------------------------------------------------------
!! * Executable part
 
!! 1. neutral value
!! ----------------
   IF(imeteo_exchtype==0) THEN
     IF(imeteo_dragtype==0) THEN
           cd = cds
       ELSEIF(imeteo_dragtype==1) THEN
         IF(wind .LT. 11.0_rsh) THEN
               cd = 0.0012_rsh
         ELSE
               cd = 0.001_rsh*(0.49_rsh+0.065_rsh*wind)
         ENDIF
       ELSEIF(imeteo_dragtype==2) THEN
            cd = 0.001_rsh*(0.63_rsh+0.066_rsh*wind)
       ELSEIF(imeteo_dragtype==3) THEN
            cd = 0.001_rsh*(0.43_rsh+0.097_rsh*wind)
       ELSEIF(imeteo_dragtype==4) THEN
            cd = charno(wind)
       ENDIF

!! 2. with dependence on tdif
!! --------------------------
 
   ELSEIF(imeteo_exchtype==1) THEN
     irow = MIN(12.0_rsh,0.4_rsh*wind+1.0_rsh)
     jcol = MIN(10.0_rsh,tdif+6.0_rsh)
     jcol = max0(1,jcol)
     f11 = drag(irow,jcol)
     f12 = drag(irow,jcol+1)
     f21 = drag(irow+1,jcol)
     f22 = drag(irow+1,jcol+1)
     z1 = 0.4_rsh*wind-irow+1.0_rsh
     z2 = tdif-jcol+6.0_rsh
     cd = (1.0_rsh-z1)*(1.0_rsh-z2)*f11+z2*(1.0_rsh-z1)*f12+z1*(1.0_rsh-z2)*f21+z1*z2*f22
   ENDIF
!
  END FUNCTION cd

!!======================================================================

  FUNCTION charno(wind)

!§E---------------------------------------------------------------------
!§E                 ***  FUNCTION charno ***
!§E
!§E ** Purpose : calculate the surface drag coefficient using
!§E             charnock s relation
!§E
!§E ** Description : uses a newton-raphson iteration procedure
!§E
!§E ** Called by : cd
!§E
!§E ** External calls :
!§E
!§E
!§E ** Reference : Charnock h., 1955. wind stress on a water surface.
!§E                quarterly journal royal meteorological society,
!§E                81, 639-640.
!§E
!§E ** History :
!§E       !          (P. Luyten)  Original code
!§E       ! 2004-12  (G. Andre)   Introduction dans MARS
!§E       ! 2006-04  (V. Garnier) Fortran 90
!§E
!§E---------------------------------------------------------------------
   USE parameters, ONLY : grav

!! * Declaration function
   REAL(kind=rsh)             :: charno
 
!! * Arguments
   REAL(kind=rsh), INTENT(in) :: wind      ! wind speed at 10m height [m/s]

!! * Local declarations
   REAL(kind=rsh), PARAMETER  :: ckar=0.4_rsh, alchar=0.014_rsh, eps=1.0e-06
   REAL(kind=rsh), PARAMETER  :: ci = 10.0_rsh*grav/alchar
   REAL(kind=rsh)             :: c, df, f, x, w10

!!----------------------------------------------------------------------
!! * Executable part

   c=LOG(ci)
   df=1.0_rsh
   f=1.0_rsh
   x = 0.001_rsh
   w10 = MAX(0.1_rsh,wind)

   DO WHILE (ABS(f/df).GT.eps)
      f = LOG(x)+ckar/SQRT(x)+2.0_rsh*LOG(w10)-c
      df = 1.0_rsh/x-0.5_rsh*ckar*x**(-1.5_rsh)
      x = ABS(x-f/df)

      IF (x==0.0_rsh) THEN
        WRITE(*,*) 'WARNING : x=0 for the drag coefficient of Charnock formulation' 
        df=1.0_rsh
        f=1.0_rsh
        x = 0.0011_rsh
        DO WHILE (ABS(f/df).GT.eps)
          f = LOG(x)+ckar/SQRT(x)+2.0_rsh*LOG(w10)-c
          df = 1.0_rsh/x-0.5_rsh*ckar*x**(-1.5_rsh)
          x = ABS(x-f/df)
          IF (x==0.0_rsh) THEN
            WRITE(*,*) 'WARNING : x=0 AGAIN ! for the drag coefficient of Charnock formulation' 
            WRITE(*,*) 'YOU SHOULD HAVE A LOOK'
            WRITE(*,*) 'charnock is set to 0.0016'
            EXIT
          END IF
        END DO
      END IF
      IF (x==0.0_rsh) THEN
        x=0.0016_rsh
        EXIT
      END IF
   END DO

   charno = x

  END FUNCTION charno

!!======================================================================

  FUNCTION ce(wind,tdif,exch,imeteo_exchtype)

!§E---------------------------------------------------------------------
!§E                 ***  FUNCTION ce ***
!§E
!§E ** Purpose : calculate thermal exchange coefficient
!§E
!§E ** Description :
!§E  - options selected by imeteo_exchtype
!§E       imeteo_exchtype = 0 -> no dependence on tdif
!§E             = 1 -> as function of tdif
!§E  - uses bilinear interpolation if imeteo_exchtype = 1
!§E
!§E ** Called by : sflx_surf
!§E
!§E ** External calls :
!§E
!§E ** Reference : - section iii-1.6.4 of the user documentation
!§E                - geernaert g.l., 1990. bulk parameterizations
!§E                  for the wind stress and heat fluxes.
!§E                  in: g.l. geernaert and w.j. plant (editors),
!§E                  surface waves and fluxes, vol. 1 - current theory.
!§E                  kluwer academic publishers, dordrecht, pp. 91-172.
!§E
!§E ** History :
!§E       !         (P. Luyten)  Original code
!§E       ! 04-12   (G. Andre)   Introduction dans MARS
!§E       ! 06-04   (V. Garnier) Fortran 90
!§E
!§E---------------------------------------------------------------------

!! * Declaration function
   REAL(kind=rsh)             :: ce

!! * Arguments
   INTEGER,       INTENT(in)  :: imeteo_exchtype      ! switch to select dependence
! on tdif
   REAL(kind=rsh),INTENT(in)  :: tdif,  &   ! air-sea temperature difference
! (used only if imeteo_exchtype=1)
                                 wind       ! wind speed at 10 m height [m/s]
   REAL(kind=rsh),INTENT(in),DIMENSION(13,11) :: &
                                 exch       ! array used for bilinear evaluation
! of thermal exchange coefficient

!! * Local declarations
   INTEGER                    :: irow, jcol
   REAL(kind=rsh)             :: f11, f12, f21, f22, z1, z2

!!----------------------------------------------------------------------
!! * Executable part

!! 1. neutral value
!! ----------------
      
   IF(imeteo_exchtype==0) THEN
      ce = 0.00113_rsh

!! 2. with dependence on tdif
!! --------------------------

   ELSEIF(imeteo_exchtype==1) THEN
     irow = MIN(12.0_rsh,0.4_rsh*wind+1.0_rsh)
     jcol = MIN(10.0_rsh,tdif+6.0_rsh)
     jcol = max0(1,jcol)
     f11 = exch(irow,jcol)
     f12 = exch(irow,jcol+1)
     f21 = exch(irow+1,jcol)
     f22 = exch(irow+1,jcol+1)
     z1 = 0.4_rsh*wind-irow+1
     z2 = tdif-jcol+6.0_rsh
     ce = (1.0_rsh-z1)*(1.0_rsh-z2)*f11+z2*(1.0_rsh-z1)*f12+z1*(1.0_rsh-z2)*f21+z1*z2*f22
   ENDIF

  END FUNCTION ce


  SUBROUTINE sflx_coef_init(drag,exch,imeteo_dragtype)

!§E---------------------------------------------------------------------
!§E                 ***  ROUTINE fluxco ***
!§E
!§E ** Purpose : calculate surface drag and exchange coefficient at fixed
!§E                  intervals of wind speed and air-sea temperature difference
!§E
!§E ** Description : - the routine solves the equations described in
!§E                    section iii-1.6.4 of the user documentation,
!§E                    by a newton-raphson iteration scheme for selected
!§E                    values of wind speed and air-sea temperzature difference
!§E                  - results are stored into the arrays drag and exch which
!§E                    will be used for the evaluation of cd and ce by
!§E                    bilinear interpolation
!§E
!§E ** Called by : flx_init
!§E
!§E ** External calls : cd, ce, error, humid
!§E
!§E ** Reference : -luyten p.j. and de mulder t., 1992. a module representing
!§E                 surface fluxes of momentum and heat. technical report
!§E                 no. 9 mast-0050-c (mumm), 30 pp.
!§E                -geernaert g.l., 1990. bulk parameterizations for the wind
!§E                 stress and heat fluxes. in: g.l. geernaert and w.j. plant
!§E                 (editors), surface waves and fluxes, vol. 1-current theory.
!§E                  kluwer academic publishers, dordrecht, pp. 91-172.
!§E
!§E ** History :
!§E       !         (P. Luyten)  Original code
!§E       ! 04-12   (G. Andre)   Introduction dans MARS
!§E       ! 06-04   (V. Garnier) Fortran 90
!§E
!§E ** Local variables :
!§E a          matrix of linear system (jacobian matrix in
!§E            newton algorithm)
!§E b          r.h.s. of linear system (newton algorithm)
!§E dx         solution of linear system (newton algorithm)
!§E x          values of cd and ce at each iteration
!§E alpha      constant in monin-obukhov function phi_m
!§E            (unstable strat.)
!§E avap       vapour pressure at air temperature           [mb]
!§E beta       constant in monin-obukhov function phi_m
!§E            (stable strat.)
!§E cdn        neutral drag coefficient
!§E cen        neutral exchange coefficient
!§E dtemp      interval of air-sea temperature difference
!§E            used in interpolation
!§E dwind      interval of wind speed used in interpolation
!§E eps        precision of the algorithm
!§E phih       monin-obukhov similarity function phi_h
!§E phim       monin-obukhov similarity function phi_m
!§E psih       vertically integrated function phi_h
!§E psim       vertically integrated function phi_m
!§E qa         specific humudity of air
!§E qs         specific humidity at sea surface temperature
!§E rhref      relative humudity (reference value)
!§E svap       vapour pressure at sea surface temperature [mbar]
!§E tair       air temperature                           [deg c]
!§E tdif       air/sea temperature difference            [deg c]
!§E tkelv      air temperature                           [deg k]
!§E tsea       sea surface temperature                   [deg c]
!§E tvir       virtual temperature                       [deg k]
!§E wind       wind speed at 10 m height                   [m/s]
!§E zeta       standard height [10m] normalised by monin-obukhov
!§E            length                                        [m]
!§E---------------------------------------------------------------------

!! * Modules used
   USE parameters, ONLY : pi
   IMPLICIT NONE

!   REAL(kind=rsh) cd, ce

!! * Arguments
   INTEGER,       INTENT(in)   :: imeteo_dragtype     ! type of drag coefficient
! formulation
   REAL(kind=rsh),INTENT(out),DIMENSION(13,11) :: &
                                  exch,  &  ! array which contains values of
! thermal exchange coefficient
! at specified intervals of wind
! speed and air/sea temperature
! difference
                                  drag      ! array which contains values of
! drag coefficient at specified
! intervals of wind speed and
! air/sea temperature difference

!! * Local declarations
   INTEGER                        :: imeteo_exchtype
   REAL(kind=rsh), PARAMETER      :: alpha = 16.0_rsh, beta = 5.0_rsh, eps = 1.0e-05, &
                                     g=9.81_rsh,ckar=0.4_rsh
   REAL(kind=rsh), PARAMETER      :: dwind=2.5_rsh,dtemp=1.0_rsh,tempref=12.0_rsh
   REAL(kind=rsh), PARAMETER      :: tsea=MIN(25.0_rsh,tempref)
!  REAL(kind=rsh), PARAMETER      :: tsea=max(5.0_rsh,tempref)
   REAL(kind=rsh), DIMENSION(2,2) :: a
   REAL(kind=rsh), DIMENSION(2)   :: b, dx, x
   REAL(kind=rsh)                 :: avap, cdn, cen, det, fac,                &
                                     f1, f2, f3, phih, phim, psih, psihx,     &
                                     psihy, psim, psimx, psimy, qa, qs, rhref,&
                                     svap, tair, tdif, tkelv, tv,             &
                                     tvir, wind, zeta, zetax, zetay
!   REAL(kind=rsh),EXTERNAL        :: ce,cd

!!----------------------------------------------------------------------
!! * Executable part

!! 1. fix parameters
!! -----------------

   rhref = 0.8_rsh
   exch(:,:)=0.0_rsh
   drag(:,:)=0.0_rsh
   imeteo_exchtype=0

!! 2. calculate coefficients at specified intervals
!! ------------------------------------------------

   DO i=2,13
     wind = (i-1)*dwind
     tdif = 0.0_rsh
     cdn = cd(wind,tdif,drag,imeteo_dragtype,imeteo_exchtype)
     cen = ce(wind,tdif,exch,imeteo_exchtype)
     f1 = SQRT(cdn)/ckar
     f2 = cen/(ckar*SQRT(cdn))
     f3 = cen/(ckar*ckar)
     DO j=1,11
       x(1) = cdn
       x(2) = cen
       tdif = j-6.0_rsh
       IF( .NOT. (tdif.GT.0.0_rsh .AND. i==2) ) THEN
         tair = tsea+tdif
         rhref = 0.8_rsh
         CALL humid(tair,rhref,avap,qa)
         rhref = 1.0_rsh
         CALL humid(tsea,rhref,svap,qs)

!! 2.1 start iteration
!! -------------------

         tkelv = tair+273.15_rsh
         tvir = tkelv*(1.0_rsh+0.61_rsh*qa)
         fac = 10.0_rsh*g*ckar/(tvir*wind*wind)
         dx(:)=1.0_rsh
         DO WHILE( ABS(dx(1))>eps .OR. ABS(dx(2))>eps ) 
           IF(x(1).LE.0.0_rsh .OR. x(2).LE.0.0_rsh) STOP 'error dans fluxco 6'

!! 2.2 monin-obukhov functions
!! ---------------------------

           tv = tdif+0.61_rsh*tkelv*(qa-qs)
           zeta = fac*x(1)**(-1.5_rsh)*x(2)*tv
           zetax = -1.5_rsh*zeta/x(1)
           zetay = zeta/x(2)
           IF(zeta.LT.0.0_rsh) THEN
             phim = (1.0_rsh-alpha*zeta)**(-0.25_rsh)
             phih = phim*phim
             tv = 1.0_rsh/phim
             psim = 2.0_rsh*LOG(0.5_rsh+0.5_rsh*tv)+LOG(0.5_rsh+0.5_rsh*tv*tv)    &
                   -2.0_rsh*ATAN(tv)+0.5_rsh*pi
             psih = 2.0_rsh*LOG(0.5_rsh+0.5_rsh/phih)
           ELSE
             phim = 1.0_rsh+beta*zeta
             phih = phim
             psim = 1.0_rsh-phim
             psih = 1.0_rsh-phih
           ENDIF
           psimx = (1.0_rsh-phim)*zetax/zeta
           psimy = (1.0_rsh-phim)*zetay/zeta
           psihx = (1.0_rsh-phih)*zetax/zeta
           psihy = (1.0_rsh-phih)*zetay/zeta

!! 2.3 evaluate matrices of linear system
!! --------------------------------------

           tv = cdn**(-0.5_rsh)+(LOG(cdn/x(1))-psim)/ckar
           b(1) = x(1)*tv*tv-1.0_rsh
           b(2) = x(2)*(1.0_rsh-f1*psim-f2*psih+f3*psim*psih)-cen
           a(1,1) = tv*tv-2.0_rsh*tv*(x(1)*psimx+1.0)/ckar
           a(1,2) = -2.0_rsh*x(1)*tv*psimy/ckar
           a(2,1) = -x(2)*(f1*psimx+f2*psihx-f3*(psimx*psih+psim*psihx))
           a(2,2) = 1.0_rsh-f1*psim-f2*psih+f3*psim*psih                      &
                    -x(2)*(f1*psimy+f2*psihy-f3*(psimy*psih+psim*psihy))

!! 2.4 determine new values
!! ------------------------

           det = a(1,1)*a(2,2)-a(1,2)*a(2,1)
           dx(1) = (b(1)*a(2,2)-b(2)*a(1,2))/det
           dx(2) = (a(1,1)*b(2)-a(2,1)*b(1))/det
           x(1) = x(1)-dx(1)
           x(2) = x(2)-dx(2)
         END DO
         drag(i,j) = x(1)
         exch(i,j) = x(2)
       END IF
     END DO
   END DO

!! 3. low or zero wind speed
!! -------------------------

   drag(2,7:11) = 0.5_rsh*drag(3,7:11)
   exch(2,7:11) = 0.5_rsh*exch(3,7:11)
   drag(1,:) = 0.0_rsh
   exch(1,:) = 0.0_rsh

  END SUBROUTINE sflx_coef_init

!!======================================================================

  SUBROUTINE humid(t,rh,pvap,q)
!§E---------------------------------------------------------------------
!§E                 ***  ROUTINE humid  ***
!§E
!§E ** Purpose : calculate specific humidity and saturated vapour pressure
!§E
!§E ** Description : - determined as function of relative humidity
!§E                  - alternative using dew point temperature instead
!§E                    of relative humidity is possible but not included
!§E
!§E ** Called by : fluxco, sflx_surf
!§E
!§E ** External calls :
!§E
!§E ** Reference : - gill a.e., 1982. atmosphere-ocean dynamics.
!§E                  international geophysics series, vol.30.
!§E                  academic press, orlando, 662 pp.
!§E
!§E ** History :
!§E       !         (P. Luyten)  Original code
!§E       ! 04-12   (G. Andre)   Introduction dans MARS
!§E       ! 06-04   (V. Garnier) Fortran 90
!§E
!§E---------------------------------------------------------------------

!! * Modules used
   IMPLICIT NONE

!! * Arguments
   REAL(kind=rsh), INTENT(in)    :: rh,   & ! relative humidity
                                    t       ! temperature               [deg c]
   REAL(kind=rsh), INTENT(out)   :: pvap, & ! saturated vapour pressure [mbar]
                                    q       !specific humidity

!! * Local declarations
   REAL(kind=rsh)             :: ex
   REAL(kind=rsh), PARAMETER  :: pair = 1013.25_rsh ! standard atmospheric pressure [mb]
!!----------------------------------------------------------------------
!! * Executable part

   ex = (0.7859_rsh+0.03477_rsh*t)/(1.0_rsh+0.00412_rsh*t)
   pvap = rh*10.0_rsh**ex
   q = 0.622_rsh*pvap/(pair-0.378_rsh*pvap)

  END SUBROUTINE humid

!!======================================================================

  SUBROUTINE sflx_turb_empty(ztemp,zsat,zwindmod,zpa,zrh,zrad,zirlos,zrain, &
                                zsenslos,zlatlos,zevap,ztosx,ztosy)

!&E---------------------------------------------------------------------
!&E                 ***  ROUTINE sflx_turb_empty  ***
!&E
!&E ** Purpose : warning
!&E
!&E---------------------------------------------------------------------
!! * Modules used

!! * Declaration function

!! * Arguments
   REAL(KIND=rsh), INTENT( IN )    :: ztemp,zsat,zwindmod,zpa,zrh, &
                                      zrad,zirlos,zrain
   REAL(KIND=rsh), INTENT( OUT )   :: zsenslos,zlatlos,zevap
   REAL(KIND=rsh), INTENT( INOUT ) :: ztosx,ztosy

!! * Local declarations
   INTEGER                       :: IERR_MPI

!!----------------------------------------------------------------------
!! * Executable part

   zsenslos=0.0_rsh
   zlatlos=0.0_rsh
   zevap=0.0_rsh
   ztosx=ztosx
   ztosy=ztosy

   WRITE(*,*) 'You do not read the turbulent atmospherical fluxes (latent ans sensible)'
   WRITE(*,*) 'and you have not chosen any formulation either'
   WRITE(*,*) 'The simulation is stopped'
   CALL MPI_FINALIZE(IERR_MPI)
   STOP

  END SUBROUTINE sflx_turb_empty

!!======================================================================

 END MODULE sflxturb
