MODULE dust_windstress_lift_mod

IMPLICIT NONE

INTEGER,SAVE :: windstress_lift_scheme ! wind stress lifting scheme number
! 0 (default): old scheme

!$OMP THREADPRIVATE(windstress_lift_scheme)

CONTAINS

  SUBROUTINE dust_windstress_lift(ngrid,nlay,nq,rho, &
                                pcdh_true,pcdh,co2ice, &
                                dqslift)

#ifndef MESOSCALE
  use tracer_mod, only: alpha_lift, radius
#else
  use tracer_mod, only: alpha_lift, radius, &
                      igcm_dust_mass, igcm_dust_number, &
                      ref_r0,r3n_q
#endif
  use comcstfi_h, only: g
  USE ioipsl_getin_p_mod, ONLY : getin_p
  IMPLICIT NONE

!=======================================================================
!
!  Dust lifting by surface winds
!  Computing flux to the middle of the first layer
!  (Called by vdifc)
!
!=======================================================================

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

  integer,intent(in) :: ngrid, nlay, nq  
  real,intent(in) ::  rho(ngrid)  ! density (kg.m-3) at surface
  real,intent(in) ::  pcdh_true(ngrid) ! Cd 
  real,intent(in) ::  pcdh(ngrid) ! Cd * |V| 
  real,intent(in) ::  co2ice(ngrid) ! surface CO2 ice (kg/m2)

  real,intent(out) ::  dqslift(ngrid,nq) !surface dust flux to mid-layer (<0 when lifing)

!   local:
!   ------
  INTEGER :: ig,iq
  REAL :: fhoriz(ngrid)  ! Horizontal dust flux 
  REAL :: ust,us
  REAL,SAVE :: stress_seuil=0.0225 ! stress lifting threshold (N.m2)
!$OMP THREADPRIVATE(stress_seuil)

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

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

#ifdef MESOSCALE
!!!! AS: In the mesoscale model we'd like to easily set
!!!! AS: ... stress for lifting
!!!! AS: you have to compile with -DMESOSCALE to do so
  REAL :: alpha
  REAL :: r0_lift
  INTEGER :: ierr
  REAL :: ulim

    OPEN(99,file='stress.def',status='old',form='formatted',iostat=ierr)
    !!! no file => default values
    IF(ierr.EQ.0) THEN
      READ(99,*) ulim !ulim = sqrt(stress_seuil/rho) avec rho = 0.02. 
                      !prendre ulim = 1.061 m/s pour avoir stress_seuil = 0.0225 
      READ(99,*) alpha
      stress_seuil = 0.02 * ulim * ulim
      write(*,*) 'USER-DEFINED threshold: ', stress_seuil, alpha
      CLOSE(99)
      alpha_lift(igcm_dust_mass) = alpha
      r0_lift = radius(igcm_dust_mass) / ref_r0
      alpha_lift(igcm_dust_number)=r3n_q* &
                              alpha_lift(igcm_dust_mass)/r0_lift**3
      write(*,*) 'set dust number: ', alpha_lift(igcm_dust_number)
    ENDIF
#endif

    if (firstcall) then
      ! get wind stress lifting scheme number (default 0)
      windstress_lift_scheme=0 ! default
      call getin_p("windstress_lift_scheme",windstress_lift_scheme)
      
      ! sanity check on available windstress_lift_scheme values
      if (windstress_lift_scheme.ne.0) then
        write(*,*) trim(rname)//" wrong value for windstress_lift_scheme:",&
                   windstress_lift_scheme
        call abort_physic(rname,"bad windstress_lift_scheme value",1)
      endif
      
      firstcall=.false.
    endif ! of if (firstcall)


    if (windstress_lift_scheme==0) then

!   ---------------------------------
!   1. Compute horizontal flux: fhoriz
!   ---------------------------------

      do ig=1,ngrid
        fhoriz(ig) = 0.      ! initialisation

!         Selection of points where surface dust is available
!         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!         if (latid(ig).ge.80.) goto 99  ! N permanent  polar caps
!         if (latid(ig).le.-80.) goto 99 ! S polar deposits
!         if  ((longd(ig).ge.-141. .and. longd(ig).le.-127.)
!    &   .and.(latid(ig).ge.12.   .and. latid(ig).le.23.))goto 99 ! olympus
!         if  ((longd(ig).ge.-125. .and. longd(ig).le.-118.)
!    &   .and.(latid(ig).ge.-12.   .and. latid(ig).le.-6.))goto 99 ! Arsia
!         if  ((longd(ig).ge.-116. .and. longd(ig).le.-109.)
!    &   .and.(latid(ig).ge.-5.   .and. latid(ig).le. 5.))goto 99 ! pavonis
!         if  ((longd(ig).ge.-109. .and. longd(ig).le.-100.)
!    &   .and.(latid(ig).ge. 7.   .and. latid(ig).le. 16.))goto 99 ! ascraeus
!         if  ((longd(ig).ge.  61. .and. longd(ig).le.  63.)
!    &   .and.(latid(ig).ge. 63. .and. latid(ig).le. 64.))goto 99 !weird point
        if (co2ice(ig).gt.0.) goto 99


!       Is the wind strong enough ?
!       ~~~~~~~~~~~~~~~~~~~~~~~~~~~
        ust = sqrt(stress_seuil/rho(ig))
        us = pcdh(ig) /  sqrt(pcdh_true(ig)) ! ustar=cd*v /sqrt(cd)
        if (us.gt.ust) then
!         If lifting ? 
!         Compute flux according to Marticorena (in fact white (1979))

          fhoriz(ig) = 2.61*(rho(ig)/g) * (us -ust) * (us + ust)**2
         endif

 99   continue
      enddo ! of do ig=1,ngrid

!     -------------------------------------
!     2. Compute vertical flux and diffusion
!     -------------------------------------
 
      do iq=1,nq
        do ig=1,ngrid
          dqslift(ig,iq)= -alpha_lift(iq)* fhoriz(ig)

!   the vertical flux replaces the turbulent diffusion term which is set to zero
!            zb(ig,1) = 0.
!c           If surface deposition by turbulence diffusion (impaction...)
!c           if(fhoriz(ig).ne.0) then
!c           zb(ig,1) = zcdh(ig)*zb0(ig,1)
!c           AMount of Surface deposition !
!c           pdqs_dif(ig,iq)=pdqs_dif(ig,iq) +
!c    &      zb(ig,1)*zq(ig,1,iq)/ptimestep
!c          write(*,*) 'zb(1)  = ' ,  zb(ig,1),zcdh(ig),zb0(ig,1)
!c

        enddo
      enddo

    endif ! of if (windstress_lift_scheme==0) 
    
  END SUBROUTINE dust_windstress_lift

END MODULE dust_windstress_lift_mod