module coare_cp_mod
  IMPLICIT NONE
  PRIVATE
  public psit_30, psiuo, coare_cp

contains

  real function psit_30(zet)
    IMPLICIT NONE
    real, intent(in) :: zet

    REAL :: x, psik, psic, f, c

    IF(zet<0) THEN
      x=(1.-(15.*zet))**.5
      psik=2.*log((1.+x)/2)
      x=(1.-(34.15*zet))**.3333
      psic=1.5*log((1.+x+x*x)/3.)-sqrt(3.)*atan((1.+2.*x)/sqrt(3.))+4.*atan(1.)/sqrt(3.)
      f=zet*zet/(1.+zet*zet)
      psit_30=(1.-f)*psik+f*psic
    else
      c=min(50.,.35*zet)
      psit_30=-((1.+2./3.*zet)**1.5+.6667*(zet-14.28)/exp(c)+8.525)
    endif

  END FUNCTION psit_30

  real function psiuo(zet)
    IMPLICIT NONE
    real, intent(in) :: zet

    REAL :: x, psik, psic, f, c

    if (zet<0) THEN
       x=(1.-15.*zet)**.25
       psik=2.*log((1.+x)/2.)+log((1.+x*x)/2.)-2.*atan(x)+2.*atan(1.)
       x=(1.-10.15*zet)**.3333
       psic=1.5*log((1.+x+x*x)/3.)-sqrt(3.)*atan((1.+2.*x)/sqrt(3.))+4.*atan(1.)/sqrt(3.)
       f=zet*zet/(1+zet*zet)
       psiuo=(1-f)*psik+f*psic
    else
     c=min(50.,.35*zet)
     psiuo=-((1.+1.0*zet)**1.0+.667*(zet-14.28)/exp(c)+8.525)
  endif

END FUNCTION psiuo


SUBROUTINE coare_cp(du,dt,dq,&
     t,q,&
     zu,zt,zq,&
     p,zi,&
     nits,&
     coeffs,rugosm,rugosh)
  !     zo,tau,hsb,hlb,&
  !     var)
  !version with shortened iteration  modified Rt and Rq
  !uses wave information wave period in s and wave ht in m
  !no wave, standard coare 2.6 charnock:  jwave=0
  !Oost et al.  zo=50/2/pi L (u*/c)**4.5 if jwave=1
  !taylor and yelland  zo=1200 h*(L/h)**4.5 jwave=2

  USE MODD_CSTS,       ONLY: XKARMAN, XG, XSTEFAN, XRD, XRV, XPI, &
                            XLVTT, XCL, XCPD, XCPV, XRHOLW, XTT, &
                            XP00

  IMPLICIT NONE

  real, intent(in) :: du,dt,dq,t,q
  real, intent(in) :: zu,zt,zq,p,zi
  integer, intent(in) :: nits
  ! real, dimension (nits), intent(out) :: zo,tau,hsb,hlb
  ! real, dimension(nits,3), intent(out) :: var
  real, dimension(3), intent(out) :: coeffs
  real, intent(out) :: rugosm
  real, intent(out) :: rugosh

  real, parameter ::  beta=1.2, von=.4, fdg = 1. ,&
       tdk = 273.16, pi = 3.141593, grav = 9.82,&
       rgas = 287.1

  integer, dimension(3) :: shape_input

  real bf, cc, rhoa, visa,&
       u10, ut, uts, ut0, ug,&
       cd10, ch10, ct, ct10,&
       charn, ribu,&
       l, l10, &
       ribcu,rr,&
       zet,zetu, zom10, zoh10, zot, zoq, zot10,&
       cd, ch, le, cpa, cpv,&
       usr,qsr,tsr,&
       zom, t_c

!, ZTWAVE, ZHWAVE, ZCWAVE, ZLWAVE




  real old_usr, old_tsr, old_qsr,tmp

  real, external :: grv
  integer i,j,k
!---------------- Rajout pour prendre en compte différent Z0 --------------------------------!
!  INTEGER :: NGRVWAVES        ! Pour le choix du z0
!  NGRVWAVES = 2
!---------------------------------------------------------------------------------------------
!-------------------- Attention Modif réalisée pas SURRR -------------------------------------
!---------------------------------------------------------------------------------------------

  Ribcu=-zu/zi/.004/Beta**3
  cpa=1004.67

  t_c = t - 273.15


  Le=(2.501-.00237*(t_c-dt))*1e6


  visa=1.326e-5*(1+6.542e-3*(t_c)+8.301e-6*t_c**2-4.84e-9*t_c**3)


  cpv=cpa*(1+0.84*Q)
  rhoa=P/(Rgas*t*(1+0.61*Q))






  ug=.2

  ut=sqrt(du*du+ug*ug)
  ut=MAX(ut , 0.1 * MIN(10.,zu) )

  u10=ut*log(10/1e-4)/log(zu/1e-4)
  usr=.035*u10                              ! turbulent friction velocity (m/s), including gustiness
  zom10=0.011*usr*usr/grav+0.11*visa/usr     ! roughness length for u (smith 88)
  Cd10=(von/log(10/zom10))**2
!  zoh10=0.40*visa/usr+1.4e-5
!  Ch10=((von**2)/((log(10/zom10))*(log(10/zoh10)))) !ammener à devenir ca
  Ch10=0.00115
  Ct10=Ch10/sqrt(Cd10)
  zot10=10/exp(von/Ct10)                    ! roughness length for t
  Cd=(von/log(zu/zom10))**2
  Ct=von/log(zt/zot10)
  CC=von*Ct/Cd

  ut0 = ut


  ut = ut0
  Ribu=grav*zu/t*(dt+.61*t*dq)/ut**2

  if (Ribu < 0) THEN
     zetu=CC*Ribu/(1+Ribu/Ribcu)
  else
     zetu=CC*Ribu*(1+27/9*Ribu/CC)
  endif

  L10=zu/zetu

  ! if (zetu .GT. 50) THEN
  !    nits=1
  ! endif

  usr=ut*von/(log(zu/zom10)-psiuo(zetu))
  tsr=dt*von*fdg/(log(zt/zot10)-psit_30(zt/L10))
  qsr=dq*von*fdg/(log(zq/zot10)-psit_30(zq/L10))
  !  tkt=.001

  ! charnock constant - lin par morceau - constant
  if ( ut <= 10. ) THEN
     charn=0.011
  else
     if (ut > 18) THEN
        charn=0.018
     else
        charn=0.011+(ut-10)/(18-10)*(0.018-0.011)
     endif
  endif

!  ZHWAVE = 0.018*ut*ut*(1.+0.015*ut)
!  ZTWAVE = 0.729*ut
!  ZCWAVE = XG*ZTWAVE/(2.*pi)
!  ZLWAVE = ZTWAVE*ZCWAVE


  !***************  bulk loop ************
  do i=1, nits

     zet=von*grav*zu/t*(tsr*(1+0.61*Q)+.61*t*qsr)/(usr*usr)/(1+0.61*Q)

!     IF (NGRVWAVES==0) THEN
!       zom = charn*usr*usr/XG + 0.11*visa/usr !Smith 1988
!     ELSE IF (NGRVWAVES==1) THEN
!       zom = (50./(2.*pi))*ZLWAVE*(usr/ZCWAVE)**4.5 &
!               + 0.11*visa/usr                       !Oost et al. 2002
!     ELSE IF (NGRVWAVES==2) THEN
!       zom = 1200.*ZHWAVE*(ZHWAVE/ZLWAVE)**4.5 &
!               + 0.11*visa/usr                       !Taulor and Yelland 2001
!     ENDIF

     zom=charn*usr*usr/grav+0.11*visa/usr

     rr=zom*usr/visa
     L=zu/zet
     zoq=min(1.15e-4,5.5e-5/rr**.6)  ! a modifier
     zot=zoq                         ! a modifier
!     zot=0.40*visa/usr+1.4e-5

     old_usr = usr
     old_tsr = tsr
     old_qsr = qsr

     usr=ut*von/(log(zu/zom)-psiuo(zu/L))
     tsr=dt*von*fdg/(log(zt/zot)-psit_30(zt/L))
     qsr=dq*von*fdg/(log(zq/zoq)-psit_30(zq/L))

     Bf=-grav/t*usr*(tsr+.61*t*qsr)

     if (Bf > 0) THEN
        ug=Beta*(Bf*zi)**.333
     else
        ug=.2
     endif

     ut=sqrt(du*du+ug*ug)
     ut=MAX(ut , 0.1 * MIN(10.,zu) )


  enddo !bulk iter loop

  ! coeffs(m1,m2,m3,m4,m5,1)=rhoa*usr*usr*du/ut !stress
  ! coeffs(m1,m2,m3,m4,m5,2)=rhoa*cpa*usr*tsr
  ! coeffs(m1,m2,m3,m4,m5,3)=rhoa*Le*usr*qsr

  tmp = (von/(log(zu/zom)-psiuo(zu/L)) ) * ut / du

  coeffs = [tmp**2,&
       von*fdg/(log(zt/zot)-psit_30(zt/L)) * tmp,&
       von*fdg/(log(zq/zoq)-psit_30(zq/L)) * tmp]

  rugosm = zom
  rugosh = zot

END SUBROUTINE  coare_cp

end module coare_cp_mod