MODULE clmain_ideal_mod
      
      IMPLICIT NONE
      
      CONTAINS
c
c
      SUBROUTINE clmain_ideal(dtime,itap,
     .                  t,u,v,
     .                  rmu0, 
     .                  ts,
     .                  ftsoil,
     .                  paprs,pplay,ppk,radsol,albe,
     .                  solsw, sollw, sollwdown, fder,
     .                  rlon, rlat, cufi, cvfi, 
     .                  debut, lafin, 
     .                  d_t,d_u,d_v,d_ts,
     .                  flux_t,flux_u,flux_v,cdragh,cdragm,
     .                  dflux_t,
     .                  zcoefh,zu1,zv1) 

c---------------------------------------------------------------
c VENUS
c Routine for a very simple idealized Planetary Boundary layer scheme:
c  - Rayleigh friction in the lowest atmospheric layer, tau=3Ed=2.6e5s
c  - Kedd=0.15 m^2/s

c S Lebonnois, 10/11/08
c---------------------------------------------------------------

      use dimphy, only: klon, klev
      use soil_mod, only: nsoilmx
      use clmain_mod, only: clvent

      IMPLICIT none
c======================================================================
c
      REAL,INTENT(IN) :: dtime ! physics time step (s)
      integer,intent(in) :: itap ! physics time step counter
      REAL,INTENT(IN) :: t(klon,klev) ! atmospheric temperature (K)
      REAL,INTENT(IN) :: u(klon,klev) ! zonal wind (m/s)
      REAL,INTENT(IN) :: v(klon,klev) ! meridional wind (m/s)
      REAL,INTENT(IN) :: paprs(klon,klev+1) ! pressure at layer boundaries (Pa)
      REAL,INTENT(IN) :: pplay(klon,klev) ! pressure at mid-layer (Pa)
      REAL,INTENT(IN) :: radsol(klon) ! Net radiative flux (positive downwards) in W/m2
! ADAPTATION GCM FOR CP(T)
      real,intent(in) :: ppk(klon,klev)
      real,intent(in) :: rlon(klon) ! longitudes (deg)
      real,intent(in) :: rlat(klon) ! latitudes (deg)
      real,intent(in) :: cufi(klon) ! mesh resolution (m)
      real,intent(in) :: cvfi(klon) ! mesh resolution (m)
      REAL,INTENT(OUT) :: d_t(klon, klev) ! temperature increment (K)
      REAL,INTENT(OUT) :: d_u(klon, klev) ! zonal wind increment (m/s)
      REAL,INTENT(OUT) :: d_v(klon, klev) ! meridional wind increment (m/s)
      REAL,INTENT(OUT) :: flux_t(klon,klev) ! latent heat flux (CpT) J/m**2/s (W/m**2)
                                            ! (positive when downwards)
      REAL,INTENT(OUT) :: dflux_t(klon) ! derivative of sensible heat flux

      REAL,INTENT(OUT) :: flux_u(klon,klev) ! zonal wind stress (kg m/s)/(m**2 s) or Pa
      REAL,INTENT(OUT) :: flux_v(klon,klev) ! meridional wind stress (kg m/s)/(m**2 s) or Pa
      REAL,INTENT(OUT) :: cdragh(klon)
      REAL,INTENt(OUT) :: cdragm(klon)
      real,intent(in) :: rmu0(klon)         ! cosine of solar zenithal angle
      LOGICAL,INTENT(IN) :: debut ! .true. if first call to physics
      LOGICAL,INTENT(IN) :: lafin ! .true. if last call to physics
c
      REAL,INTENT(IN) :: ts(klon) ! surface temperature (K)
      REAL,INTENT(OUT) :: d_ts(klon) ! surface temperature increment (K)
      REAL,INTENT(INOUT) :: albe(klon) ! albedo of the surface
C
      REAL,INTENT(IN) :: fder(klon)
      REAL,INTENT(IN) :: sollw(klon), solsw(klon), sollwdown(klon)
cAA
      REAL,INTENT(OUT) :: zcoefh(klon,klev)
      REAL,INTENT(OUT) :: zu1(klon) ! zonal wind in 1st layer (m/s)
      REAL,INTENT(OUT) :: zv1(klon) ! meridional wind in 1st layer (m/s)
cAA
c$$$ PB ajout pour soil
      REAL,INTENT(INOUT) :: ftsoil(klon,nsoilmx) ! subsurface temperatures (K)

      REAL ytsoil(klon,nsoilmx)
c======================================================================
      REAL yts(klon)
      REAL yalb(klon)
      REAL yu1(klon), yv1(klon)
      real ysollw(klon), ysolsw(klon), ysollwdown(klon)
      real yfder(klon), ytaux(klon), ytauy(klon)
      REAL yrads(klon)
C
      REAL y_d_ts(klon)
      REAL y_d_t(klon, klev)
      REAL y_d_u(klon, klev), y_d_v(klon, klev)
      REAL y_flux_t(klon,klev)
      REAL y_flux_u(klon,klev), y_flux_v(klon,klev)
      REAL y_dflux_t(klon)
      REAL ycoefh(klon,klev), ycoefm(klon,klev)
      REAL yu(klon,klev), yv(klon,klev)
      REAL yt(klon,klev)
      REAL ypaprs(klon,klev+1), ypplay(klon,klev), ydelp(klon,klev)
c
      REAL ycoefm0(klon,klev), ycoefh0(klon,klev)

      real yzlay(klon,klev),yzlev(klon,klev+1)
      real yteta(klon,klev)
      real ykmm(klon,klev+1),ykmn(klon,klev+1)
      real ykmq(klon,klev+1)
      real yustar(klon),y_cd_m(klon),y_cd_h(klon)
c
      REAL u1lay(klon), v1lay(klon)
      REAL delp(klon,klev)
      INTEGER i, k
      INTEGER ni(klon), knon, j
      
c======================================================================
      REAL zx_alf1, zx_alf2 ! ambient values used for extrapolation
c======================================================================
c
      REAL zt, zdelta, zcor
C
      real taurelax

c=========================================================
c DEBUT
c=========================================================
          
      DO k = 1, klev   ! thickness of atmospheric layers
      DO i = 1, klon
         delp(i,k) = paprs(i,k)-paprs(i,k+1)
      ENDDO
      ENDDO
      DO i = 1, klon  ! wind in the first layer
ccc         zx_alf1 = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
         zx_alf1 = 1.0
         zx_alf2 = 1.0 - zx_alf1
         u1lay(i) = u(i,1)*zx_alf1 + u(i,2)*zx_alf2
         v1lay(i) = v(i,1)*zx_alf1 + v(i,2)*zx_alf2
      ENDDO
c
c initialisation:
c
      DO i = 1, klon
         cdragh(i) = 0.0
         cdragm(i) = 0.0
         dflux_t(i) = 0.0
         zu1(i) = 0.0
         zv1(i) = 0.0
      ENDDO
      yts = 0.0
      yalb = 0.0
      yfder = 0.0
      ytaux = 0.0
      ytauy = 0.0
      ysolsw = 0.0
      ysollw = 0.0
      ysollwdown = 0.0
      yu1 = 0.0
      yv1 = 0.0
      yrads = 0.0
      ypaprs = 0.0
      ypplay = 0.0
      ydelp = 0.0
      yu = 0.0
      yv = 0.0
      yt = 0.0
      y_flux_u = 0.0
      y_flux_v = 0.0
      y_d_ts = 0.0
      y_d_t = 0.0
      y_d_u = 0.0 
      y_d_v = 0.0
      y_flux_t = 0.0
C$$ PB
      y_dflux_t = 0.0
      ytsoil = 999999.
      DO i = 1, klon
         d_ts(i) = 0.0
      ENDDO
      flux_t = 0.
      flux_u = 0.
      flux_v = 0.
      DO k = 1, klev
      DO i = 1, klon
         d_t(i,k) = 0.0
         d_u(i,k) = 0.0
         d_v(i,k) = 0.0
         zcoefh(i,k) = 0.0
      ENDDO
      ENDDO
c
c identify indexes:
      DO j = 1, klon
         ni(j) = j
      ENDDO
      knon = klon

      DO j = 1, knon
      i = ni(j)
        yts(j) = ts(i)
        yalb(j) = albe(i)
        yfder(j) = fder(i)
        ytaux(j) = flux_u(i,1)
        ytauy(j) = flux_v(i,1)
        ysolsw(j) = solsw(i)
        ysollw(j) = sollw(i)
        ysollwdown(j) = sollwdown(i)
        yu1(j) = u1lay(i)
        yv1(j) = v1lay(i)
        yrads(j) =  ysolsw(j)+ ysollw(j)
        ypaprs(j,klev+1) = paprs(i,klev+1)
      END DO
C
c$$$ PB ajour pour soil
      DO k = 1, nsoilmx
        DO j = 1, knon
          i = ni(j)
          ytsoil(j,k) = ftsoil(i,k)
        END DO  
      END DO 
      DO k = 1, klev
      DO j = 1, knon
      i = ni(j)
        ypaprs(j,k) = paprs(i,k)
        ypplay(j,k) = pplay(i,k)
        ydelp(j,k) = delp(i,k)
        yu(j,k) = u(i,k)
        yv(j,k) = v(i,k)
        yt(j,k) = t(i,k)
      ENDDO
      ENDDO
c
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c RAYLEIGH FRICTION (implicit scheme) in the first layer
c Ref: PhD of C. Lee Oxford 2006
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      taurelax = 2.6e5
      yu1 = yu1 / (1+dtime/taurelax)
      yv1 = yv1 / (1+dtime/taurelax)
      yu(:,1) = yu(:,1) / (1+dtime/taurelax)
      yv(:,1) = yv(:,1) / (1+dtime/taurelax)

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c Coefficient for vertical diffusion
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      ycoefm = 0.15

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c compute diffusion for winds "u" and "v"
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      CALL clvent(knon,dtime,yu1,yv1,ycoefm,yt,yu,ypaprs,ypplay,ydelp,
     s            y_d_u,y_flux_u)
      CALL clvent(knon,dtime,yu1,yv1,ycoefm,yt,yv,ypaprs,ypplay,ydelp,
     s            y_d_v,y_flux_v)

c for the coupling
      ytaux = y_flux_u(:,1)
      ytauy = y_flux_v(:,1)

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c no diffusion for "q" and "h"
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      ycoefh = 0.

c=========================
c END: compute tendencies
c=========================

      DO j = 1, knon
         i = ni(j)
         d_ts(i) = y_d_ts(j)
         albe(i) = yalb(j)
         cdragh(i) = cdragh(i) + ycoefh(j,1)
         cdragm(i) = cdragm(i) + ycoefm(j,1)
         dflux_t(i) = dflux_t(i) + y_dflux_t(j)
         zu1(i) = zu1(i) + yu1(j)
         zv1(i) = zv1(i) + yv1(j)
      END DO

c$$$ PB ajout pour soil
      DO k = 1, nsoilmx
        DO j = 1, knon
         i = ni(j)
         ftsoil(i, k) = ytsoil(j,k)
        ENDDO
      END DO
      
      DO k = 1, klev
        DO j = 1, knon
         i = ni(j)
         flux_t(i,k) = y_flux_t(j,k)
         flux_u(i,k) = y_flux_u(j,k)
         flux_v(i,k) = y_flux_v(j,k)
         d_t(i,k) = d_t(i,k) + y_d_t(j,k)
         d_u(i,k) = d_u(i,k) + y_d_u(j,k)
         d_v(i,k) = d_v(i,k) + y_d_v(j,k)
         zcoefh(i,k) = zcoefh(i,k) + ycoefh(j,k)
        ENDDO
      ENDDO

c --------------------
c TEST!!!!! PAS DE MELANGE PAR TURBULENCE !!!
c       d_u = 0. 
c       d_v = 0.
c       flux_u = 0.
c       flux_v = 0.
c --------------------

c     print*,"y_d_t apres clqh=",y_d_t(klon/2,:)

      END SUBROUTINE clmain_ideal

      END MODULE clmain_ideal_mod