source: CONFIG/publications/ICOLMDZORINCA_CO2_Transport_GMD_2023/INCA/src/INCA_SRC/dustsource.F90 @ 6610

!$Id: dustsource.F90 104 2008-12-23 10:28:51Z acosce $
!! =========================================================================
!! INCA - INteraction with Chemistry and Aerosols
!! 
!! Copyright Laboratoire des Sciences du Climat et de l'Environnement (LSCE)
!!           Unite mixte CEA-CNRS-UVSQ
!! 
!! Contributors to this INCA subroutine:
!! 
!! Michael Schulz, LSCE, Michael.Schulz@cea.fr
!! 
!! Anne Cozic, LSCE, anne.cozic@cea.fr
!! Yann Meurdesoif, LSCE, yann.meurdesoif@cea.fr
!!
!! Olivier Boucher, LMD
!! introduction of a Weibull subgrid scale wind distribution
!!
!! This software is a computer program whose purpose is to simulate the 
!! atmospheric gas phase and aerosol composition. The model is designed to be
!! used within a transport model or a general circulation model. This version 
!! of INCA was designed to be coupled to the LMDz GCM. LMDz-INCA accounts
!! for emissions, transport (resolved and sub-grid scale), photochemical 
!! transformations, and scavenging (dry deposition and washout) of chemical 
!! species and aerosols interactively in the GCM. Several versions of the INCA
!! model are currently used depending on the envisaged applications with the 
!! chemistry-climate model. 
!! 
!! This software is governed by the CeCILL  license under French law and
!! abiding by the rules of distribution of free software.  You can  use, 
!! modify and/ or redistribute the software under the terms of the CeCILL
!! license as circulated by CEA, CNRS and INRIA at the following URL
!! "http://www.cecill.info". 
!! 
!! As a counterpart to the access to the source code and  rights to copy,
!! modify and redistribute granted by the license, users are provided only
!! with a limited warranty  and the software's author,  the holder of the
!! economic rights,  and the successive licensors  have only  limited
!! liability. 
!! 
!! In this respect, the user's attention is drawn to the risks associated
!! with loading,  using,  modifying and/or developing or reproducing the
!! software by the user in light of its specific status of free software,
!! that may mean  that it is complicated to manipulate,  and  that  also
!! therefore means  that it is reserved for developers  and  experienced
!! professionals having in-depth computer knowledge. Users are therefore
!! encouraged to load and test the software's suitability as regards their
!! requirements in conditions enabling the security of their systems and/or 
!! data to be ensured and,  more generally, to use and operate it in the 
!! same conditions as regards security. 
!! 
!! The fact that you are presently reading this means that you have had
!! knowledge of the CeCILL license and that you accept its terms.
!! =========================================================================

#include <inca_define.h>


#ifdef AER
SUBROUTINE DUSTSOURCE(delt)
  !     -----------------------------------------------------------------------
  !
  !     Author
  !     ------
  !     Michael Schulz  
  !       / Laboratoire des Sciences du Climat et de l'Environnement / Saclay
  !       9. June 2002
  !
  !     purpose
  !     -------
  !     compute interactively emission flux of mineral dust
  !
  !     interface
  !     ---------
  !      input 
  !       delt           time step duration       [sec]
  !
  !      from aerosol modules
  !       zspeed         wind speed at 10 m                 [m/s]
  !       landmask       sea land mask                      []
  !       ztempc         surface air temperature            [°C]
  !       zprecip        precipitation amount               [kg m-2 s-1]
  !       zprecipinsoil  proxy for soil humidity            [mm]
  !       rop            particle density
  !       rhv            source strength factor             [kg s2 m-5]
  !       wth            threshold velocity                 [m/s]
  !       cly            clay content                       [%]
  !
  !       srcsigma       sigma of aerosol mode 
  !       srcsigmaln
  !
  !      variables saved
  !       entered        input fields read flag
  !
  !      local variables
  !       fdist          number/mass factor                 [#/kg]
  !       nbdays         soil drying rate                   [days/mm]
  !       mask           dust emission area mask            [] 
  !
  !      output 
  !       eflux          emission of dust per tracer        [kg m-2 s-1]
  !       
  !     method
  !     ------
  !     
  !     The dust aerosol emitted is described as a lognormal size
  !     distribution with a mass median radius of 2.5 um and a standard
  !     deviation of s = 2. Such a source size distribution was found to
  !     provide a consistent mineral aerosol transport model, which
  !     provides enough mineral dust for long range transport from the
  !     Sahara to the Atlantic and the Mediterranean. Dust plumes
  !     simulated with such a model provide aerosol optical thickness
  !     values very close to those observed with satellites [Schulz et
  !     al., 1998] [Guelle et al., 2000]. Using a fixed lognormal
  !     distribution at the source implies that number concentrations can
  !     be derived from mass emission flux and imposed source size
  !     distribution [Schulz et al., 1998].
  !
  !     While surface wind speed is produced interactively or read from
  !     ECMWF winds (see routine aerosol_meteo_calc), the source function
  !     incorporates also spatial information on three parameters a) the
  !     source strength factor (rhv); b) the threshold velocity (wth) and
  !     c) the clay fraction of the uppermost soil layer (cly). 
  !
  !         eflux = max(rhv*(landwind-wth)*landwind**2,0.)*landmask
  !
  !     This work is based on an offline dust source developped by
  !     [Claquin, 1999].  He used a global TM3 dust simulation with a
  !     dustsource in all arid regions to obtain an approximate dust
  !     field. In a second step he compared regional averages of TOMS
  !     observations of aerosol index to this simulation, to obtain a
  !     source strength factor for a given region [see also Balkanski et
  !     al. [submitted , 2002]. The threshold velocity was obtained by
  !     comparing the FAO soil database to the maps of Marticorena [1995]
  !     to derive by this procedure a soil specific threshold velocity for
  !     desert soils. Combined with the high resolution FAO dataset this
  !     results in a global map of threshold velocities for dust rise in
  !     arid and semi-arid regions. The offline source formulation was
  !     converted to an interactive dust source for this work [Schulz und
  !     Timmreck, in prep] including prognostic variables produced by the
  !     GCM, such as wind speed, precipitation and surface temperature.
  !
  !     Dust emission is thought to be inhibited, when the surface soil is
  !     wettened or frozen. Drying of desert soils can be parameterised as
  !     a function of recent accumulated precipitation, surface
  !     temperature. The clay content of the soil is a key parameter for
  !     this process, since clay retains water for a longer period. While
  !     Claquin [1999] used an offline calculation to define the periods
  !     when the soil is available for emission of dust, we prefer an
  !     interactive approach here. Soil wettening and drying rate are
  !     computed from GCM calculated preceipitation and surface
  !     temperatures. The clay content is assembled from the FAO soil
  !     data. During freezing conditions, no drying of the soil and no
  !     emissions are assumed.  The parameterisation uses the conditions
  !     given by Claquin [1999] and computes an evaporation rate acting on
  !     the accumulated precipitation amount [Schulz und Timmreck, in
  !     prep]. This might overestimate the time which is needed to dry off
  !     the top soil, since no run-off or losses to the ground water
  !     reservoirs are assumed.
  !     (parameterised M. Schulz LSCE 13.2.2001)
  !     -----------------------------------------------------------------------


  USE SPECIES_NAMES
  USE AEROSOL_METEO, only : landwind,ztempc,landmask,zprecip,zprecipinsoil
  USE AEROSOL_MOD, ONLY : srcsigma,srcsigmaln,rop,rhv,wth,cly,cimode
  USE SFLX, ONLY : eflux
  USE AEROSOL_DIAG, ONLY : DUSTE
  USE INCA_DIM
  USE MOD_INCA_PARA
  USE CHEM_CONTROLS, ONLY : lafin
  USE CONST_MOD, ONLY : pi
  USE PRINT_INCA
  USE PARAM_CHEM

  IMPLICIT NONE


  REAL, INTENT(in) ::     delt


  ! source mass median diameter [m]
  REAL ::  srcmmd_CIDUSTM = 2.5e-6

  ! temp variables
  REAL    :: fdist          ! Number/Mass factor to compute number flux
  REAL    :: nbday(PLON)
  REAL    :: clyfac(PLON)
  REAL    :: avgdryrate(PLON),drying(PLON)
  INTEGER :: mask(PLON) 
  REAL    :: efluxECP_MPI(PLON_MPI)
  REAL    :: efluxECP_MPI2(PLON_MPI)
  REAL    :: efluxECP(PLON)
  REAL,ALLOCATABLE :: wth_glo(:)
  REAL,ALLOCATABLE :: rhv_glo(:)
  REAL,ALLOCATABLE :: cly_glo(:)
  REAL,ALLOCATABLE, SAVE :: zprecipinsoil_glo(:)
!$OMP THREADPRIVATE(zprecipinsoil_glo) 
  INTEGER :: i,j
  LOGICAL, SAVE   :: entered = .FALSE.
  LOGICAL, SAVE   :: ECM_HR  = .TRUE.   ! Interpolation from 1°x1°
!$OMP THREADPRIVATE (entered, ECM_HR)

 ! Weibull parameters and variables
 !--number of bins in weibull distribution (nwb=1 reduces to having no distribution)
  INTEGER, PARAMETER :: nwb=12
  INTEGER kwb
  REAL auxreal, wind10ms, zlandwind, weilambda, pdfu, probu, pdfcum
  REAL, PARAMETER :: rpi = 4.*ATAN(1.)


  mask(:)=0


  IF( .not. LMDZ_10m_winds ) THEN
!$OMP MASTER
     CALL DUSTECMWF_XIOS(efluxECP_MPI)
!$OMP END MASTER

     CALL scatter_omp(efluxECP_MPI,efluxECP)

  ENDIF

  ! read of geographical source info only once
  IF( .NOT. entered ) THEN

     CALL read_threshold_map() 
     ALLOCATE(zprecipinsoil_glo(nbp_glo))

!$OMP MASTER
     IF (is_mpi_root) THEN
          OPEN(53,file='precipinsoil.dat',status='old',form='formatted',err=999)
          READ(53,'(G18.10)') (zprecipinsoil_glo(i),i=1,nbp_glo)
          GOTO 1000
999       zprecipinsoil_glo=0.
1000      CLOSE(53)
      ENDIF
!$OMP END MASTER

      CALL scatter(zprecipinsoil_glo,zprecipinsoil)
      DEALLOCATE(zprecipinsoil_glo)

      entered = .TRUE.

  END IF

  WHERE (abs(cly).LT.9990. .AND. landmask.LT.0.5 .AND. abs(wth).LT.9990.)
      mask = 1
  ENDWHERE

  ! precipitation is accumulated everywhere, where dust emission is in
  ! principal possible
  zprecipinsoil = (zprecipinsoil + zprecip*delt) * mask


  ! surface soil wettening computations      
  ! old function, did accumulate zprecipinsoil at the edges of the arid zones
  ! thus avoiding dust emissions there
  !      where (ztempc(:,1).gt.0. .and. mask.eq.1)
  !        clyfac=MIN(16.,cly*0.4+8.)
  !        nbday=MIN(16.,CLYFAC/CLYFAC**((ZTEMPC(:,1)-10.)/15.))
  !        zprecipinsoil=max(0.,zprecipinsoil - delt/86400./nbday)
  !      endwhere
  !
  ! new soil wettening formulation f(clay content and soil temperature)
  ! Dependent on clay content and interpolated to obtain 
  ! clyfac, which is the maximim number of dry days after a rain event and corresponds to a 
  ! maximum amount of water hold in top soil [mm] at the lowest temperatures [5°C].
  ! The avg drying rate scales the drying to keep the area an arid area with an aridity limit
  ! assumed to be 300 mm precipitation per year.
  ! More rapid drying only depends on temperature and 
  ! fits the 'dry days table' of Claquin to a function with two fitting parameters.

  avgdryrate=300./365.*delt/86400. ! [mm/timestep] average drying rate at lowest temperature        
  WHERE (mask.EQ.1)
      clyfac=MIN(16.,cly*0.4+8.)                    ! [days] number of days before dry
      !  = max amount of water hold in top soil [mm]
      drying=avgdryrate*EXP(0.03905491*EXP(0.17446*ztempc(:,1))) ! [mm]
      zprecipinsoil=MIN(MAX(0.,zprecipinsoil - drying),clyfac)   ! [mm] remaining water in top soil
  ENDWHERE

  IF (lafin) THEN
     ALLOCATE(zprecipinsoil_glo(nbp_glo))
     zprecipinsoil_glo(:) = 0. 
     CALL gather(zprecipinsoil,zprecipinsoil_glo)
!$OMP MASTER
     IF (is_mpi_root) THEN
          OPEN(54,file='reprecipinsoil.dat',status='new',form='formatted')
          WRITE(54,'(G18.10)') (zprecipinsoil_glo(i),i=1,nbp_glo)
          CLOSE(54)
      ENDIF
!$OMP END MASTER
  DEALLOCATE(zprecipinsoil_glo)
  ENDIF

  !  conversion of mass median diameter**3 to diameter of average mass**3,
  !  => exp(1.5*srcsigmaln(asmode)**2)**3 = exp(4.5*srcsigmaln(asmode)**2)
  ! Source size distribution is assumed to be constant
  fdist=1./pi*6./rop(id_CIDUSTM) &
     /srcmmd_CIDUSTM**3 *EXP(4.5*srcsigmaln(cimode)**2)



  !****************source function        
  !--------unit kg m**-2 s-1 
  ! computation of flux 
  IF( LMDZ_10m_winds ) THEN
      WHERE (ztempc(:,1).GT.0. .AND. mask.EQ.1 .AND. zprecipinsoil.LT.1e-10) 
          eflux(:,id_CIDUSTM) = MAX(rhv*(landwind-wth)*landwind**2,0.)
      ENDWHERE
  ENDIF


  IF( (.not. LMDZ_10m_winds) .and. ECM_HR ) THEN
      WHERE (ztempc(:,1).GT.0. .AND. mask.EQ.1 .AND. zprecipinsoil.LT.1e-10) 
          DUSTE = MAX(rhv*(landwind-wth)*landwind**2,0.)
          eflux(:,id_CIDUSTM) = efluxECP
      ENDWHERE

  endif


  ! and for the number flux, allowing other aerosol species sources to add
  ! to the number flux
  eflux(:,id_CIN) = eflux(:,id_CIN)+ eflux(:,id_CIDUSTM)*fdist
  
  
  RETURN
END SUBROUTINE dustsource

#endif