 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.
   !!
   !!======================================================================

#include "toolcpp.h"

   !! * 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
#ifdef key_sflx_turb_luyten
     MODULE PROCEDURE sflx_turb_luyten
#elif key_sflx_turb_default
     MODULE PROCEDURE sflx_turb_default
#elif key_sflx_turb_large
     MODULE PROCEDURE sflx_turb_large
#elif key_sflx_turb_ayina
     MODULE PROCEDURE sflx_turb_ayina
#elif key_sflx_turb_fairall
     MODULE PROCEDURE sflx_turb_fairall
#else
     MODULE PROCEDURE sflx_turb_empty
#endif
   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

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

#     include "sflx_turb_large.h"

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

#     include "sflx_turb_ayina.h"

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

#     include "sflx_turb_fairall.h"

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

  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

#ifdef key_windstress_module
         ztosx=rhoair*cd(zwindmod,tdif,drag,imeteo_dragtype,imeteo_exchtype)*zwindmod*zwindmod
#endif

   ! 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

#ifdef key_windstress_module
         ztosx=rhoair*cd(zwindmod,tdif,drag,imeteo_dragtype,imeteo_exchtype)*zwindmod*zwindmod
#endif

   ! 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 MPI_FINALIZE(IERR_MPI)
   STOP

  END SUBROUTINE sflx_turb_empty

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

 END MODULE sflxturb
