source: trunk/LMDZ.MARS/libf/phymars/dustdevil.F90 @ 3582

MODULE dustdevil_mod

IMPLICIT NONE

INTEGER,SAVE :: dust_devil_scheme ! dust devil scheme number
! 0 (default): old Renno et al. 1998 (JAS 55, 3244-3252) scheme

!$OMP THREADPRIVATE(dust_devil_scheme)

CONTAINS

  SUBROUTINE dustdevil(ngrid,nlay,nq, pplev,pu,pv,pt, ptsurf,pq2, &
                       pdqdev,pdqs_dev)

    use tracer_mod, only: alpha_devil
    use surfdat_h, only: z0_default
    USE comcstfi_h, ONLY: g, cpp, r, rcp
    USE mod_phys_lmdz_para, ONLY: is_master, bcast
    USE ioipsl_getin_p_mod, ONLY : getin_p
    IMPLICIT NONE

!=======================================================================
!
!  Parameterization of dust devil activities
!  to estimate dust lifting
!  The dust devil activity is estimated based on 
!  Renno et al. 1998 (JAS 55, 3244-3252)  
!
!  It is proportional to (1-b)*Fs
!
!  With b= [ps**(rcp+1) - ptop**(rcp+1)] / [(ps-ptop)*(rcp+1)* ps**rcp]
!  with ptop pressure of the top of the boundary layer
!       rcp = R/cp
!
!  and Fs the surface sensible heat flux = Cd*|U|*(T(1) -Tsurf)
!       
!=======================================================================

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

    INTEGER,INTENT(IN) :: ngrid,nlay
    REAL,INTENT(IN) :: pplev(ngrid,nlay+1)
    REAL,INTENT(IN) :: pt(ngrid,nlay)
    REAL,INTENT(IN) :: pu(ngrid,nlay)
    REAL,INTENT(IN) :: pv(ngrid,nlay)
    REAL,INTENT(IN) :: pq2(ngrid,nlay+1)
    REAL,INTENT(IN) :: ptsurf(ngrid)

!    Traceurs :
    INTEGER,INTENT(IN) :: nq 
    real,intent(out) :: pdqdev(ngrid,nlay,nq) 
    real,intent(out) :: pdqs_dev(ngrid,nq) 
      
!   local:
!   ------

    INTEGER :: ig,l,iq
    real,save :: Cd
!$OMP THREADPRIVATE(Cd)
    real :: z1

    LOGICAL,SAVE :: firstcall=.true.
!$OMP THREADPRIVATE(firstcall)

    character(len=20),parameter :: rname="dustdevil"

    REAL :: devila(ngrid)
    integer :: ltop(ngrid)
    real :: b,rho,Fs,wind

    REAL,SAVE :: q2top=.5 ! value of q2 at the top of PBL
    REAL,SAVE :: seuil=.3 ! value of minimum dust devil activity
                          ! for dust lifting      
!$OMP THREADPRIVATE(q2top,seuil)

!-----------------------------------------------------------------------
!    initialisation
!    --------------

    ! AS: OK firstcall absolute

    IF (firstcall) THEN

      ! get scheme number (default=0, old scheme)
      dust_devil_scheme=0
      call getin_p("dust_devil_scheme",dust_devil_scheme)

      ! sanity check on available dust_devil_scheme values
      if (dust_devil_scheme.ne.0) then
        write(*,*) trim(rname)//" wrong value for dust_devil_scheme:",&
                   dust_devil_scheme
        call abort_physic(rname,"bad dust_devil_scheme value",1)
      endif

      if(is_master) then

        write(*,*) 'In dustdevil :'
        write(*,*) '    q2top= ',q2top,'     seuil= ', seuil 

! A rough estimation of the horizontal drag coefficient Cd
! (the scale heigh is taken to be 13 km since we are typically
! dealing with daytime temperature around 250K.
! 
        z1= -0.5*13.e3*log(pplev(1,2)/pplev(1,1))
        Cd = (0.4/log(z1/z0_default))**2

!       Temporaire
!       open(77,file='devil')
      
      endif !is_master

      ! share the info with all cores
      CALL bcast(z1)
      CALL bcast(Cd)

      firstcall=.false.

    ENDIF ! firstcall 


!-----------------------------------------------------------------------
! Initialisation of outputs
    do iq=1,nq
      do l=1,nlay
        do ig=1,ngrid
          pdqdev(ig,l,iq)= 0
        end do
      end do
    end do
    pdqs_dev(1:ngrid,1:nq)=0

    if (dust_devil_scheme==0) then
!-----------------------------------------------------------------------
!      Determining the top of the boundary layer
!      ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      do ig=1,ngrid
        do  l=2,nlay-1
          if (pq2(ig,l).lt.q2top)then
            ltop(ig)=l
            goto 99
          end if
        enddo
 99     continue

!        ***************************************
!c        if (ptsurf(ig).gt.255)then
!         write(*,*) 'tsurf, ztop (km): ', ig,
!     &   ptsurf(ig), -12.*log(pplev(ig,ltop(ig))/pplev(ig,1))
!         if ((ptsurf(ig).gt.50.).and.(
!     &      (-12.*log(pplev(ig,ltop(ig))/pplev(ig,1))).gt.60.))then
!            do l=1,nlay
!             write(*,*) l,
!     &       -12.*log(pplev(ig,l)/pplev(ig,1)),pq2(ig,l)
!            end do
!         end if
!c        end if
!        ***************************************
     
      enddo

!        ***************************************
!        do ig=100,148
!           write(*,*)'time,tsurf,ztop', localtime(ig),ptsurf(ig),
!    &      -12.*log(pplev(ig,ltop(ig))/pplev(ig,1))
!        end do
!        ***************************************


!     Calculation : dust devil intensity
!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      do ig=1,ngrid

! --------------------------------------------------
! 1) Version 1 : sensible heat flux using actual wind :
!      Wind magnitude:
!      wind = sqrt(pu(ig,1)**2+pv(ig,1)**2)
!
! --------------------------------------------------
! 2) Version 2 : sensible heat flux using  wind = 15 m/s
        wind = 15.
! ----------------------------------------------------
!      Density :
        rho=pplev(ig,1)/(R*pt(ig,1))

!        Sensible heat flux (W.m-2) (>0 if up)
        Fs= rho*cpp*Cd * wind * (ptsurf(ig) -pt(ig,1))
        b= (pplev(ig,1)**(rcp+1) - pplev(ig,ltop(ig))**(rcp+1)) / &
           ( (pplev(ig,1)-pplev(ig,ltop(ig)))*(rcp+1)*pplev(ig,1)**rcp)

!     b_diag(ig) = b     ! TEMPORAIRE (pour diagnostique)

!   Energy flux available to drive dust devil (W.m-2) : (1-b)*Fs
!   Dust devil activity : 
        devila(ig)= max( 0. , (1-b)*Fs - seuil ) 
      enddo
!   
!     lifted dust (kg m-2 s-1)  (<0 when lifting)
!     ~~~~~~~~~~  
      do iq=1,nq
        do ig=1,ngrid
          pdqs_dev(ig,iq)= - alpha_devil(iq) * devila(ig)
        end do
      end do 

  !   Injection of dust in the atmosphere (up to the top of pbl)
  !   ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      do iq=1,nq
        do ig=1,ngrid
          if (devila(ig).ne.0.) then
            do l=1,ltop(ig)
               pdqdev(ig,l,iq)=-pdqs_dev(ig,iq)*g/ &
                               (pplev(ig,1)-pplev(ig,ltop(ig)))
            end do
          end if
        end do
      end do
    
    endif ! of if (dust_devil_scheme==0)
     
  END SUBROUTINE dustdevil

END MODULE dustdevil_mod