
! $Id$

      SUBROUTINE clc_core_cp(du, dt, dq, t1, q1, &
     				zu, zt, zq, &
     					P,&
     				n_it, mixte, &
    			 coeffs, rugosm, rugosh)

      IMPLICIT NONE


  !  INTEGER     :: i,j,k
  REAL, INTENT(IN) :: du,dt,dq,t1,q1
  REAL, INTENT(IN) :: zu,zt,zq, P
  INTEGER, INTENT(IN):: n_it
  LOGICAL, INTENT(IN) :: mixte
  REAL, DIMENSION(3), INTENT(OUT) :: coeffs
  real, intent(out) :: rugosm
  real, intent(out) :: rugosh


!  REAL :: du,dt,t1,dq,q1
  REAL :: dt_u, dq_u



  REAL :: cd, ch, cd_rt, cd_n10, ce_n10, ce,&
       u_N, X,&
       rho,cpa,le

  REAL :: usr,tsr,qsr
  REAL, DIMENSION(3) :: zeta                           ! dans l'ordre: it courant, it precedente, var
  REAL, DIMENSION(2,3) :: psi
  REAL, PARAMETER :: g = 9.81, k=.41, sqrt_k = .640312 !avec g= gravité, k = cste von karman, sqrt_k = racine(k)
  REAL :: logzu_10,logzt_10,logzq_10,logzt_zu,logzq_zu,&
       u,cd_n10_rt,ch_n10_rt,ch_n10,chi,z_0m,z_0h,tv
  INTEGER :: i,j,m1,m2,m3,m4,m5


  REAL :: phih,phid,phie

  
  REAL, parameter :: alpha=2.7e-3,&           !alpha = a1 ds formulation smith 
       beta = 1.42e-4,&                       !beta = a2 ds la formulation smith
       gamma = 7.64e-5,&        !Large and yaeger 2006              !gamma = a3 ds formulation smith
!       gamma = 13.09e-3,&         !Large and yaeger 2004
       q0 = 1.64474                           ! q0 dzfinis mais pas reutilisée ds la suite ??

  REAL, dimension(3) :: z
  integer, dimension(2) :: tmp_shape
  z = [zu, zt, zq]
 


  logzu_10 = LOG(zu/10.)
  logzt_10 = LOG(zt/10.)
  logzq_10 = LOG(zq/10.)
  logzt_zu = logzt_10 - logzu_10
  logzq_zu = logzq_10 - logzu_10

  cpa = 1004.67                               !!! pas reutilise ds nvelle formulation...

  tv = t1 * (1+.608*q1)                       !!!!! tv pas reutiliser ds nvelle formulation...
  rho = p / (287.1 * t1 * (1.+.61*q1))        !!!! rho pas reutiliser ds nvelle formulation....


  le = (2.501-.00237*(t1-273.15-dt))*1.e6     !!!!!! le pas reutiliser ds nvelle formulation....
  
  u = max(du,.5)                              !!!! u = vitesse du vent 1 er niveau
  u_N = u                                     !!!!!!! u_N utilisée pour le calcul des premier calcul...
  

  !!!!! initialisation des parametres et premier calcul de z_0 cd et ch


  if (mixte) then 
     z_0m = 1.e-4                             !!!!! pour l'initialisation on fixe z_0 a 10-4
     cd_n10 = k**2 / (log(10/z_0m))**2        !!!!! on calcul ensuite cd a 10 m d'après smith 
     cd_n10_rt = sqrt(cd_n10)   
  else
! Large and yaeger 2006
     cd_n10 = alpha / u_n + beta + gamma * u_n !!!!! calcul de cd pour core pur
! Large and yaeger 2004
!     cd_n10 = 1E-3 * (2.7/U_n + 0.142 + U_n/13.09) 
     cd_n10_rt = sqrt(cd_n10)                  !!!!!!
     z_0m = 10. * exp(-k/cd_n10_rt)            !!!!!!! on calcul le premier z_0
  endif


  if ( dt > 0. ) then
     ch_n10 = 0.018 * cd_n10_rt                !!!!!!! si regime stable calcul du ch
!     z_0h = 10. * exp(-k/ch_n10)
  else
     ch_n10 = 0.0327 * cd_n10_rt               !!!!!!! si regime instable calcul du ch
!     z_0h = 10. * exp(-k/ch_n10)
  endif

  ce_n10 = 3.46e-2 *cd_n10_rt                  !!!!!! calcul ensuite de ce


  dt_u = dt                                    !!! def
  dq_u = dq                                    !!! def

  cd = cd_n10                                  !!! def
  ch = ch_n10                                  !!! def
  ce = ce_n10                                  !!! def

  cd_rt = sqrt(cd)                             !!! def

  !  open(7,file="err_rel.dat",position="append")


!!!!!!!!! début des itérations....

  do i = 1,n_it

     usr = cd_rt * u 
     tsr = ch / cd_rt *  dt_u
     qsr = ce / cd_rt * dq_u


     zeta = k*g / (usr**2) * tsr / t1&
          * z

     zeta = sign( min(abs(zeta),10.0), zeta )


!!!! calcul du zeta pour connaitre le signe et donc en deduire le regime...




     IF ( zeta(1) > 0 ) THEN !!! si regime stable alors valeurs pr chi et psi

        chi = 0.018 


        psi(1,1) = -5.*zeta(1)
        psi(2,1) = psi(1,1)

     ELSE !!!!! si regime instable alors valeurs pr psi et chi !!!

        chi = 0.0327

        X = sqrt( sqrt( 1-16.*zeta(1) ))
        psi(1,1) =  3.1415695 / 2 + 2 * LOG( ( 1 + X ) / 2 ) + LOG( (1 + X**2) / 2 ) - 2*ATAN(X)
        psi(2,1) =  2.*LOG( (1 + X**2) / 2 )

     END IF

     do j =2,3
        IF ( zeta(j) > 0 ) THEN
           psi(1,j) = -5.*zeta(j)
           psi(2,j) = psi(1,1)

        ELSE
           X = sqrt( sqrt( 1-16.*zeta(j) ))
           psi(1,j) =  3.1415695 / 2 + 2 * LOG( ( 1 + X ) / 2 ) + LOG( (1 + X**2) / 2 ) - 2*ATAN(X)
           psi(2,j) =  2.*LOG( (1 + X**2) / 2 )

        END IF
     enddo


     dt_u = dt - tsr / k * &
          ( logzt_zu + psi(2,1) - psi(2,2) )

     dq_u = dq - qsr / k * &
          ( logzq_zu + psi(2,1) - psi(2,3) )



     u_N = u / ( 1 + cd_n10_rt / k * ( logzu_10 - psi(1,1)) )

     if (mixte) then
        z_0m=0.018*usr**2/9.81 + 0.11*14e-6 / (usr)
!        z_0h=(0.11*14e-6 / (usr)) + 1.4e-5 
        cd_n10 = k**2 / (log(10/z_0m))**2
        cd_n10_rt = sqrt(cd_n10)
!        ch_n10 = k**2 / ((log(10/z_0m))*(log(10/z_0h)))
!        ch_n10_rt = sqrt(ch_n10)
        ch_n10 = 1.e-3
     else
! Large and yaeger 2006
        cd_n10 = alpha / u_n + beta + gamma * u_n !!!!! calcul de cd pour core pur
! Large and yaeger 2004
!        cd_n10 = 1E-3 * (2.7/U_n + 0.142 + U_n/13.09) 
        !        cd_n10_dun = -alpha / u_n**2 + gamma
        cd_n10_rt = sqrt(cd_n10)
        z_0m = 10. * exp(-k/cd_n10_rt) !!!!! c'est sur cette ligne la qu'il y a un bug... ou sur la ligne du u_N
!        ch_n10 = chi * cd_n10_rt
!        z_0h = 10. * exp(-k/ch_n10)
!        ch_n10_rt = sqrt(ch_n10)
        ch_n10 = chi * cd_n10_rt
     endif


     ce_n10 = 3.46e-2*cd_n10_rt



     phid = 1+cd_n10_rt/k * (logzu_10 - psi(1,1))
     phih = 1+chi/k * (logzu_10 - psi(2,1))
     phie = 1+3.46e-2/k * (logzu_10 - psi(2,1))


     cd_rt = cd_n10_rt / abs( phid )
     ch = ch_n10 / ( phih * phid )
     ce = ce_n10 / ( phie * phid )

     cd=cd_rt**2

  END DO

  usr = cd_rt * u
  tsr = ch / cd_rt *  dt_u
  qsr = ce / cd_rt * dq_u


  ! coeffs(:,m) = [ &
  !      rho * usr**2,&
  !      rho * cpa * usr *tsr,&
  !      rho * le * usr *qsr]

  coeffs = [ cd_rt**2, ch , ce ]

  rugosm =  z_0m
  rugosh =  z_0h

      RETURN
      END subroutine clc_core_cp
