MODULE vdif_mod
  IMPLICIT NONE
  SAVE
  PRIVATE
  
  REAL, PARAMETER :: karman=0.4
  REAL :: lmixmin=100., emin_turb=1e-8  

  PUBLIC :: vdif, lmixmin, emin_turb
  
CONTAINS
  
  PURE SUBROUTINE vdif_cd( ngrid,pz0,pg,pz,pu,pv,pts,ph,pcdv,pcdh)
    !=======================================================================
    !
    !   Subject: computation of the surface drag coefficient using the
    !   -------  approch developed by Loui for ECMWF.
    !
    !   Author: Frederic Hourdin  15 /10 /93
    !   -------
    !
    !   Arguments:
    !   ----------
    !
    !   inputs:
    !   ------
    !     ngrid            size of the horizontal grid
    !     pz0(ngrid)       roughness length ?
    !     pg               gravity (m s -2)
    !     pz(ngrid)        height of the first atmospheric layer
    !     pu(ngrid)        u component of the wind in that layer
    !     pv(ngrid)        v component of the wind in that layer
    !     pts(ngrid)       surfacte temperature
    !     ph(ngrid)        potential temperature T*(p/ps)^kappa
    !
    !   outputs:
    !   --------
    !     pcdv(ngrid)      Cd for the wind
    !     pcdh(ngrid)      Cd for potential temperature
    !
    !=======================================================================
    !
    !-----------------------------------------------------------------------
    !   Declarations:
    !   -------------
    
    !   Arguments:
    !   ----------
    
    INTEGER, INTENT(IN) ::  ngrid
    REAL, INTENT(IN)  ::  pg, pz0(ngrid), pz(ngrid), &
         pu(ngrid),pv(ngrid), pts(ngrid),ph(ngrid)
    REAL, INTENT(OUT) :: pcdv(ngrid),pcdh(ngrid)
    
    !   Local:
    !   ------
    
    REAL, PARAMETER :: b=5., c=5., d=5., umin2=1e-12, &
         c2b=2.*b, c3bc=3.*b*c, c3b=3.*b
    INTEGER ig
    REAL zu2,z1,zri,zcd0,zz    
    
    !-----------------------------------------------------------------------
    !   couche de surface:
    !   ------------------
    
    !     DO ig=1,ngrid
    !        zu2=pu(ig)*pu(ig)+pv(ig)*pv(ig)+umin2
    !        pcdv(ig)=pz0(ig)*(1.+sqrt(zu2))
    !        pcdh(ig)=pcdv(ig)
    !     ENDDO
    !     RETURN
    
    
!!!! WARNING, verifier la formule originale de Louis!
    DO ig=1,ngrid
       zu2=pu(ig)*pu(ig)+pv(ig)*pv(ig)+umin2
       zri=pg*pz(ig)*(ph(ig)-pts(ig))/(ph(ig)*zu2)
       z1=1.+pz(ig)/pz0(ig)
       zcd0=karman/log(z1)
       zcd0=zcd0*zcd0*sqrt(zu2)
       IF(zri.LT.0.) THEN
          z1=b*zri/(1.+c3bc*zcd0*sqrt(-z1*zri))
          pcdv(ig)=zcd0*(1.-2.*z1)
          pcdh(ig)=zcd0*(1.-3.*z1)
       ELSE
          zz=sqrt(1.+d*zri)
          pcdv(ig)=zcd0/(1.+c2b*zri/zz)
          pcdh(ig)=zcd0/(1.+c3b*zri*zz)
       ENDIF
    ENDDO
    
    !-----------------------------------------------------------------------
    
  END SUBROUTINE vdif_cd
  
  PURE SUBROUTINE vdif_k(ngrid,nlay,   &
       ptimestep,pg,pzlev,pzlay,pz0,pu,pv,ph,pkv,pkh)
    ! FIXME : pkh := pkv
    USE planet
    INTEGER, INTENT(IN) :: ngrid,nlay
    REAL, INTENT(IN) ::  ptimestep, pg, &
         pzlay(ngrid,nlay),pzlev(ngrid,nlay+1), pz0(ngrid), &
         pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay)
    REAL, INTENT(OUT) :: pkv(ngrid,nlay+1),pkh(ngrid,nlay+1)
    
    INTEGER ig,il
    REAL zdu,zdv,zri,zdvodz2,zdz,z1,lmix
    
    DO ig=1,ngrid
       pkv(ig,1)=0.
       pkh(ig,1)=0.
       pkv(ig,nlay+1)=0.
       pkh(ig,nlay+1)=0.
    ENDDO
    
    DO il=2,nlay
       DO ig=1,ngrid
          z1=pzlev(ig,il)+pz0(ig)
          lmix=karman*z1/(1.+karman*z1/lmixmin)
          !           lmix=lmixmin
          ! WARNING test lmix=lmixmin
          zdu=pu(ig,il)-pu(ig,il-1)
          zdv=pv(ig,il)-pv(ig,il-1)
          zdz=pzlay(ig,il)-pzlay(ig,il-1)
          zdvodz2=(zdu*zdu+zdv*zdv)/(zdz*zdz)
          IF(zdvodz2.LT.1.e-5) THEN
             pkv(ig,il)=lmix*sqrt(emin_turb)
          ELSE
             zri=2.*pg*(ph(ig,il)-ph(ig,il-1)) &
                  / (zdz* (ph(ig,il)+ph(ig,il-1)) *zdvodz2  )
             pkv(ig,il)= &
                  lmix*sqrt(MAX(lmix*lmix*zdvodz2*(1-zri/.4),emin_turb))
          ENDIF
          pkh(ig,il)=pkv(ig,il)
          !           IF(ig.EQ.ngrid/2+1) PRINT*,il,lmix,pkv(ig,il),
          !    s      zdu,zdv,zdz,zdvodz2,ph(ig,il)+ph(ig,il-1),
          !    s      lmix*lmix*zdvodz2*(1-zri/.4),emin_turb,zri,ph(ig,il)-ph(ig,il-1),
          !    s      ph(ig,il),ph(ig,il-1)
       ENDDO
    ENDDO
    
  END SUBROUTINE vdif_k

  
  SUBROUTINE multipl(n,x1,x2,y)
!=======================================================================
!   multiplication de deux vecteurs
!=======================================================================
    INTEGER n,i
    REAL x1(n),x2(n),y(n)

    DO i=1,n
       y(i)=x1(i)*x2(i)
    END DO
  END SUBROUTINE multipl

  SUBROUTINE vdif(ngrid,nlay,ptime, &
       ptimestep,pcapcal,pz0, &
       pplay,pplev,pzlay,pzlev, &
       pu,pv,ph,ptsrf,pemis, &
       pdufi,pdvfi,pdhfi,pfluxsrf, &
       pdudif,pdvdif,pdhdif,pdtsrf,pq2,pq2l, &
       lwrite)
    USE phys_const
    
!=======================================================================
!
!   Diffusion verticale
!   Shema implicite
!   On commence par rajouter au variables x la tendance physique
!   et on resoult en fait:
!      x(t+1) =  x(t) + dt * (dx/dt)phys(t)  +  dt * (dx/dt)difv(t+1)
!
!   arguments:
!   ----------
!
!   entree:
!   -------
!
!
!=======================================================================

!
!   arguments:
!   ----------

    INTEGER ngrid,nlay
    REAL ptime,ptimestep
    REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
    REAL pzlay(ngrid,nlay),pzlev(ngrid,nlay+1)
    REAL pu(ngrid,nlay),pv(ngrid,nlay),ph(ngrid,nlay)
    REAL ptsrf(ngrid),pemis(ngrid)
    REAL pdufi(ngrid,nlay),pdvfi(ngrid,nlay),pdhfi(ngrid,nlay)
    REAL pfluxsrf(ngrid)
    REAL pdudif(ngrid,nlay),pdvdif(ngrid,nlay),pdhdif(ngrid,nlay)
    REAL pdtsrf(ngrid),pcapcal(ngrid),pz0(ngrid)
    REAL pq2(ngrid,nlay+1),pq2l(ngrid,nlay+1)
    LOGICAL lwrite
    !
    !   local:
    !   ------
    
    INTEGER ilev,ig,ilay,nlev
    INTEGER unit,ierr,it1,it2
    INTEGER cluvdb,putdat,putvdim,setname,setvdim
    REAL z4st,zdplanck(ngrid),zu2
    REAL zkv(ngrid,nlay+1),zkh(ngrid,nlay+1)
    REAL zcdv(ngrid),zcdh(ngrid)
    REAL zu(ngrid,nlay),zv(ngrid,nlay)
    REAL zh(ngrid,nlay)
    REAL ztsrf2(ngrid)
    REAL z1(ngrid),z2(ngrid)
    REAL za(ngrid,nlay),zb(ngrid,nlay)
    REAL zb0(ngrid,nlay)
    REAL zc(ngrid,nlay),zd(ngrid,nlay)
    REAL zcst1
    
    !
    !-----------------------------------------------------------------------
    !   initialisations:
    !   ----------------
    
    nlev=nlay+1
    
    !   computation of rho*dz and dt*rho/dz=dt*rho**2 g/dp:
    !   with rho=p/RT=p/ (R Theta) (p/ps)**kappa
    !   ---------------------------------
    
    DO ilay=1,nlay
       DO ig=1,ngrid
          za(ig,ilay) = (pplev(ig,ilay)-pplev(ig,ilay+1))/g
       ENDDO
    ENDDO
    
    zcst1=4.*g*ptimestep/(r*r)
    DO ilev=2,nlev-1
       DO ig=1,ngrid
          zb0(ig,ilev)=pplev(ig,ilev) &
               *(pplev(ig,1)/pplev(ig,ilev))**rcp &
               /(ph(ig,ilev-1)+ph(ig,ilev))
          zb0(ig,ilev)=zcst1*zb0(ig,ilev)*zb0(ig,ilev) &
               / (pplay(ig,ilev-1)-pplay(ig,ilev))
       ENDDO
    ENDDO
    DO ig=1,ngrid
       zb0(ig,1)=ptimestep*pplev(ig,1)/(r*ptsrf(ig))
    ENDDO
    IF(lwrite) THEN
       ig=ngrid/2+1
       PRINT*,'Pression (mbar) ,altitude (km),u,v,theta, rho dz'
       DO ilay=1,nlay
          WRITE(*,*) .01*pplay(ig,ilay),.001*pzlay(ig,ilay), &
               pu(ig,ilay),pv(ig,ilay),ph(ig,ilay),za(ig,ilay)
       ENDDO
       PRINT*,'Pression (mbar) ,altitude (km),zb'
       DO ilev=1,nlay
          WRITE(*,*) .01*pplev(ig,ilev),.001*pzlev(ig,ilev), &
               zb0(ig,ilev)
       ENDDO
    ENDIF
    
    !-----------------------------------------------------------------------
    !   2. ajout des tendances physiques:
    !   ------------------------------
    
    DO ilev=1,nlay
       DO ig=1,ngrid
          zu(ig,ilev)=pu(ig,ilev)+pdufi(ig,ilev)*ptimestep
          zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep
          zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep
       ENDDO
    ENDDO
    
    !-----------------------------------------------------------------------
    !   3. calcul de  cd :
    !   ----------------
    !
    CALL vdif_cd( ngrid,pz0,g,pzlay,pu,pv,ptsrf,ph,zcdv,zcdh)
    CALL vdif_k(ngrid,nlay,ptimestep,g,pzlev,pzlay,pz0,pu,pv,ph,zkv,zkh)
    
    DO ig=1,ngrid
       zu2=pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1)
       zcdv(ig)=zcdv(ig)*sqrt(zu2)
       zcdh(ig)=zcdh(ig)*sqrt(zu2)
    ENDDO
    
    IF(lwrite) THEN
       PRINT*
       PRINT*,'Diagnostique diffusion verticale'
       print*,'LMIXMIN',lmixmin
       PRINT*,'coefficients Cd pour v et h'
       PRINT*,zcdv(ngrid/2+1),zcdh(ngrid/2+1)
       PRINT*,'coefficients K pour v et h'
       DO ilev=1,nlay
          PRINT*,zkv(ngrid/2+1,ilev),zkh(ngrid/2+1,ilev)
       ENDDO
    ENDIF
    
    !-----------------------------------------------------------------------
    !   integration verticale pour u:
    !   -----------------------------
    
    CALL multipl((nlay-1)*ngrid,zkv(1,2),zb0(1,2),zb(1,2))
    CALL multipl(ngrid,zcdv,zb0,zb)
    DO ig=1,ngrid
       z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
       zc(ig,nlay)=za(ig,nlay)*zu(ig,nlay)*z1(ig)
       zd(ig,nlay)=zb(ig,nlay)*z1(ig)
    ENDDO
    
    DO ilay=nlay-1,1,-1
       DO ig=1,ngrid
          z1(ig) = 1./(za(ig,ilay)+zb(ig,ilay) &
               + zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
          zc(ig,ilay) = (za(ig,ilay)*zu(ig,ilay) &
               + zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
          zd(ig,ilay) = zb(ig,ilay)*z1(ig)
       ENDDO
    ENDDO
    
    DO ig=1,ngrid
       zu(ig,1)=zc(ig,1)
    ENDDO
    DO ilay=2,nlay
       DO ig=1,ngrid
          zu(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zu(ig,ilay-1)
       ENDDO
    ENDDO
    
    !-----------------------------------------------------------------------
    !   integration verticale pour v:
    !   -----------------------------
    !
    DO ig=1,ngrid
       z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
       zc(ig,nlay)=za(ig,nlay)*zv(ig,nlay)*z1(ig)
       zd(ig,nlay)=zb(ig,nlay)*z1(ig)
    ENDDO
    
    DO ilay=nlay-1,1,-1
       DO ig=1,ngrid
          z1(ig)=1./(za(ig,ilay)+zb(ig,ilay) &
               + zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
          zc(ig,ilay)=(za(ig,ilay)*zv(ig,ilay) &
               + zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
          zd(ig,ilay)=zb(ig,ilay)*z1(ig)
       ENDDO
    ENDDO
    
    DO ig=1,ngrid
       zv(ig,1)=zc(ig,1)
    ENDDO
    DO ilay=2,nlay
       DO ig=1,ngrid
          zv(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zv(ig,ilay-1)
       ENDDO
    ENDDO
    
    !-----------------------------------------------------------------------
    !   integration verticale pour h:
    !   -----------------------------
    !
    CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2))
    CALL multipl(ngrid,zcdh,zb0,zb)
    DO ig=1,ngrid
       z1(ig)=1./(za(ig,nlay)+zb(ig,nlay))
       zc(ig,nlay)=za(ig,nlay)*zh(ig,nlay)*z1(ig)
       zd(ig,nlay)=zb(ig,nlay)*z1(ig)
    ENDDO
    
    DO ilay=nlay-1,1,-1
       DO ig=1,ngrid
          z1(ig)=1./(za(ig,ilay)+zb(ig,ilay) &
               + zb(ig,ilay+1)*(1.-zd(ig,ilay+1)))
          zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay) &
               + zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig)
          zd(ig,ilay)=zb(ig,ilay)*z1(ig)
       ENDDO
    ENDDO
    
    !-----------------------------------------------------------------------
    !   rajout eventuel de planck dans le shema implicite:
    !   --------------------------------------------------
    
    z4st=4.*5.67e-8*ptimestep
    DO ig=1,ngrid
       zdplanck(ig)=z4st*pemis(ig)*ptsrf(ig)*ptsrf(ig)*ptsrf(ig)
    ENDDO
    
    !-----------------------------------------------------------------------
    !   calcul le l'evolution de la temperature du sol:
    !   -----------------------------------------------
    
    DO ig=1,ngrid
       z1(ig) = pcapcal(ig)*ptsrf(ig)+cpp*zb(ig,1)*zc(ig,1) &
            + zdplanck(ig)*ptsrf(ig)+ pfluxsrf(ig)*ptimestep
       z2(ig) = pcapcal(ig)+cpp*zb(ig,1)*(1.-zd(ig,1))+zdplanck(ig)
       ztsrf2(ig) = z1(ig)/z2(ig)
       zh(ig,1) = zc(ig,1)+zd(ig,1)*ztsrf2(ig)
       pdtsrf(ig) = (ztsrf2(ig)-ptsrf(ig))/ptimestep
    ENDDO
    
    !-----------------------------------------------------------------------
    !   integration verticale finale:
    !   -----------------------------
    
    DO ilay=2,nlay
       DO ig=1,ngrid
          zh(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zh(ig,ilay-1)
       ENDDO
    ENDDO
    
    !-----------------------------------------------------------------------
    !   calcul final des tendances de la diffusion verticale:
    !   -----------------------------------------------------
    
    DO ilev = 1, nlay
       DO ig=1,ngrid
          pdudif(ig,ilev) = ( zu(ig,ilev) &
               - (pu(ig,ilev)+pdufi(ig,ilev)*ptimestep) )/ptimestep
          pdvdif(ig,ilev) = ( zv(ig,ilev) &
               - (pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep) )/ptimestep
          pdhdif(ig,ilev) = ( zh(ig,ilev) &
               - (ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep) )/ptimestep
       ENDDO
    ENDDO
    
    IF(lwrite) THEN
       PRINT*
       PRINT*,'Diagnostique de la diffusion verticale'
       PRINT*,'h avant et apres diffusion verticale'
       PRINT*,ptsrf(ngrid/2+1),ztsrf2(ngrid/2+1)
       DO  ilev=1,nlay
          PRINT*,ph(ngrid/2+1,ilev),zh(ngrid/2+1,ilev)
       ENDDO
    END IF

  END SUBROUTINE vdif
  
END MODULE vdif_mod