ecearth3/sources/tm5mp/proj/cbm4/emission_dust.F90

!
#define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr
#define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if
#define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if
#define IF_NOTOK_MDF(action) if (status/=0) then; TRACEBACK; action; call MDF_CLose(dust_FileID,status); status=1; return; end if
!
#include "tm5.inc"
!
!-----------------------------------------------------------------------------
!                    TM5                                                     !
!-----------------------------------------------------------------------------
!BOP
!
! !MODULE:      EMISSION_DUST
!
! !DESCRIPTION: 
!\\
! 
! AEROCOM emissions
! -----------------
!
! Example of file content:
!
!    netcdf dust200001 {
!    dimensions:
!            longitude = 360 ;
!            latitude = 180 ;
!            time = 31 ;
!
!    variables:
!            float longitude(longitude) ;
!                    longitude:TITLE = "longitude" ;
!                    longitude:UNITS = "degrees_east" ;
!            float latitude(latitude) ;
!                    latitude:TITLE = "latitude" ;
!                    latitude:UNITS = "degrees_north" ;
!            float time(time) ;
!                    time:TITLE = "day in Jan" ;
!            float gridbox_area(latitude, longitude) ;
!                    gridbox_area:TITLE = "total area in gridbox" ;
!                    gridbox_area:UNITS = "m2/gridbox" ;
!            float total_flux(time, latitude, longitude) ;
!                    total_flux:TITLE = "total daily_dust_flux" ;
!                    total_flux:UNITS = "kg/gridbox" ;
!            float mode2_radius(time, latitude, longitude) ;
!                    mode2_radius:TITLE = "number mode-radius for log-normal distr. (std.dev=1.59) for accumu mode" ;
!                    mode2_radius:UNITS = "um" ;
!            float mode3_radius(time, latitude, longitude) ;
!                    mode3_radius:TITLE = "number mode-radius for log-normal distr. (std.dev=2.00) for coarse mode" ;
!                    mode3_radius:UNITS = "um" ;
!            float mode2_number(time, latitude, longitude) ;
!                    mode2_number:TITLE = "daily dust particles in accumu mode (.1-1um)" ;
!                    mode2_number:UNITS = "number/gridbox" ;
!            float mode3_number(time, latitude, longitude) ;
!                    mode3_number:TITLE = "daily dust particles in coarse mode (1-6um)" ;
!                    mode3_number:UNITS = "number/gridbox" ;
!            float mode2_flux(time, latitude, longitude) ;
!                    mode2_flux:TITLE = "daily dust flux in accumu mode (0.1-1um sizes)" ;
!                    mode2_flux:UNITS = "kg/gridbox" ;
!            float mode3_flux(time, latitude, longitude) ;
!                    mode3_flux:TITLE = "daily dust flux in coarse mode (1-6um sizes)" ;
!                    mode3_flux:UNITS = "kg/gridbox" ;
!
!    // global attributes:
!                    :filename = "ginoux2000_dust.nc" ;
!                    :title = "daily dust fluxes for Jan 2000" ;
!                    :history = "created by S.Kinne Nov/2003" ;
!                    :institution = "MPI-Meteorology, Hamburg" ;
!    }
!
!
! Dust emissions Tegen/Vignati
! ----------------------------
!
! From the 'readme.txt' that accompanies the files:
!
!  "Dust emission fields are prepared using an application of the 
!  Ina Tegen model (Tegen et al. JGR 107, D21, 2002), 
!  extended by B. Heinhold (JGR, 112, 2007)
!  and adapted by E. Vignati using the ECMWF fields as input, 
!  so they are coherent with the TM5 input.
!
!  The work of Tegen et al. (2002) is heavily based on
!  Marticorena and Bergametti (JGR, 1995) - MB95
!
!  The fields are prepared without the soil moisture because 
!  it was not trivial to use the proper fields available in the 
!  TM5 ECMWF fields.
!
!  Please contact Elisabetta Vignati (elisabetta.vignati@jrc.it) 
!  for a proper ackowlegment to use in case the fields are used 
!  for publication."
!
! Example of file content:
!
!  netcdf dust200201 {
!  dimensions:
!          lon = 360 ;
!          lat = 180 ;
!          time = 31 ;
!          nmonth = 1 ;
!
!  variables:
!          float mode2_mass(time, lat, lon) ;
!          float mode2_number(time, lat, lon) ;
!          float mode3_mass(time, lat, lon) ;
!          float mode3_number(time, lat, lon) ;
!          float lon(lon) ;
!                  lon:units = "[degrees from -180]" ;
!          float lat(lat) ;
!                  lat:units = "[degrees from -90]" ;
!          float gridbox_area(lat, lon) ;
!                  gridbox_area:units = "[square m]" ;
!          float days(time) ;
!                  days:units = "[day/month]" ;
!
!  // global attributes:
!                  :filename = "E:\\global_emissions\\dust_online\\dust200201.nc" ;
!                  :lat = 180 ;
!                  :lon = 360 ;
!                  :nmonth = 1 ;
!                  :days = 31 ;
!                  :message = "dust for each mode mass in kg/gridbox and number in number/gridbox daily fields" ;
!  }
!
!
!
! Online dust emissions based on Tegen/Vignati/Strunk
! ---------------------------------------------------
!
! Please read the section above for background information about the underlying 
! approach. An improved and modified online implementation has been accomplished 
! from which. It can be activated by setting 
!
!    input.emis.dust : ONLINE
!
! in the rc-file. An additional netcdf file is needed for some input parameters. 
! The path to which needs to be defined in the key
!
!    input.emis.dust.dir : /ms_perm/TM/TM5/emissions/other/Dust_online/onlinedust.nc
! 
! For every time step there will be particles emitted, scaled to monthly 
! amounts (both mass and numbers) in order to keep compliance with assumptions
! about the aerosol emissions in sedimentation.F90. 
!
!-----------------------------------------------------------------------
!\\
! !INTERFACE: 
!
MODULE EMISSION_DUST
  !
  ! !USES:
  !
  USE dims,          ONLY : okdebug, nlon360, nlat180
  USE GO,            ONLY : gol, goErr, goPr 
  USE global_types,  ONLY : d3_data, emis_data
  USE tm5_distgrid,  ONLY : dgrid, get_distgrid, scatter, gather
  USE partools,      ONLY : isRoot
  USE mo_aero_m7,    ONLY : ddust
  USE emission_data, ONLY : emis_input_dir_dust, emis_input_dust
  USE meteodata,     ONLY : tv_dat, sd_dat, lsmask_dat
  USE meteodata,     ONLY : u10m_dat, v10m_dat, nveg
  USE meteodata,     ONLY : spm_dat, t2m_dat

  IMPLICIT NONE
  PRIVATE
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  PUBLIC   :: calc_emission_dust    ! online computation
  PUBLIC   :: emission_dust_declare !
  PUBLIC   :: emission_dust_init
  PUBLIC   :: emission_dust_done
  PUBLIC   :: read_emission_dust
  !
  ! !PRIVATE DATA MEMBERS:
  !
  CHARACTER(len=*), PARAMETER  :: mname = 'emission_dust'
  INTEGER                      :: dust_FileID                   ! file id for input parameters

  ! parameters for online emission input file ("onlinedust.nc")
  ! fields on 1x1 deg grid
  INTEGER, PARAMETER :: nsoilph =  5, &
                        nfpar   = 12, &
                        nz0     = 13  ! number of {soilph, par, z0} fields
                                      ! entry nz0 indicates the annual mean.
  REAL, DIMENSION(:,:),   ALLOCATABLE :: soil_type, pot_source, cult
  REAL, DIMENSION(:,:,:), ALLOCATABLE :: z0, fpar, soilph

  ! local arrays during runtime
  REAL, DIMENSION(:),     ALLOCATABLE ::  Uth
  REAL, DIMENSION(:,:),   ALLOCATABLE ::  srel, srelV, su_srelV  

  REAL, DIMENSION(:,:),   ALLOCATABLE :: snowcover, desert, umin2, alpha, c_eff, area, lai_eff
  REAL, DIMENSION(:,:),   ALLOCATABLE :: fnum_ai, flux_ai, fnum_ci, flux_ci
  !REAL, DIMENSION(:),     ALLOCATABLE :: fluxtyp, dpk
  REAL, DIMENSION(:),     ALLOCATABLE :: fluxtyp
  !REAL, DIMENSION(:),     ALLOCATABLE :: dbmin, dbmax, fluxtot, fdust
  REAL, DIMENSION(:),     ALLOCATABLE :: fluxtot, fdust 
  !
  ! !REVISION HISTORY: 
  !          2005 - Elisabetta Vignati - changed for coupling with M7
  !    1 Oct 2010 - Achim Strunk       - added Tegen-Vignati option
  !      Nov 2011 - Achim Strunk       - included online dust emissions
  !   26 Jun 2012 - Ph. Le Sager       - adapted for lon-lat MPI domain decomposition
  !    1 Jul 2013 - Ph. Le Sager       - added init routine
  !    April 2015 - T. van Noije       - corrections in online dust emissions
  !
  ! !REMARKS:
  !
  !EOP
  !------------------------------------------------------------------------

CONTAINS

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_DUST_INIT
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_DUST_INIT( status )
    !
    ! !USES:
    !
    USE dims,      ONLY : iglbsfc
    USE meteo,     ONLY : Set
    !
    ! !OUTPUT PARAMETERS:
    !
    INTEGER, INTENT(out) ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Jul 2013 - Ph Le Sager - v0
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC
    integer :: iveg
    
    CALL set( lsmask_dat(iglbsfc), status, used=.TRUE. )
    CALL set(     spm_dat(iglbsfc), status, used=.TRUE. )
    CALL set(    t2m_dat(iglbsfc), status, used=.TRUE. )
    CALL set(     sd_dat(iglbsfc), status, used=.TRUE. )
    CALL set(   u10m_dat(iglbsfc), status, used=.TRUE. )
    CALL set(   v10m_dat(iglbsfc), status, used=.TRUE. )
    DO iveg = 1, nveg
       CALL set( tv_dat(iglbsfc,iveg), status, used=.TRUE. )
    END DO
    
  END SUBROUTINE EMISSION_DUST_INIT
  !EOC


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    CALC_EMISSION_DUST
  !
  ! !DESCRIPTION: online dust computation. See module !DESCRIPTION above.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE CALC_EMISSION_DUST( status )
    !
    ! !USES:
    !
    USE dims,          ONLY : newday, idate, iglbsfc, sec_month, im, jm, lm, nregions
    USE toolbox,       ONLY : coarsen_emission
    USE partools,      ONLY : localComm, MPI_INFO_NULL
    USE chem_param,    ONLY : mode_aci, mode_coi
    USE emission_data, ONLY : emis_mass, emis_number, emis_temp, msg_emis
    USE emission_data, ONLY : vd_class_name_len, emission_vdist_by_sector
    USE binas,         ONLY : rgas, xmair
    USE binas,         ONLY : grav, vkarman, pi
    USE MDF,           ONLY : MDF_Open,           MDF_Close,       MDF_Inq_VarID
    USE MDf,           ONLY : MDF_Var_Par_Access, MDF_INDEPENDENT, MDF_COLLECTIVE
    USE MDF,           ONLY : MDF_Get_Var,        MDF_READ,        MDF_NETCDF4
    USE mo_aero_m7,    ONLY : iacci, icoai, sigma
    !
    ! !OUTPUT PARAMETERS:
    !
    INTEGER, INTENT(out)        ::  status
    !
    ! !REVISION HISTORY: 
    !      Nov 2011 - Achim Strunk - v0
    !    5 Jul 2012 - Ph. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:
    !   (1) this routine is basically the "declare" routine for the ONLINE case.
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    CHARACTER(len=*), PARAMETER ::  rname = mname//'/calc_emission_dust'

    ! --- local ----------------------------------------

    REAL, DIMENSION(:,:), ALLOCATABLE    :: glb_arr
    TYPE(d3_data),   DIMENSION(nregions) :: emis3d
    TYPE(emis_data), DIMENSION(nregions) :: emis_glb
    CHARACTER(len=vd_class_name_len)     :: splittype
    INTEGER                              :: imr, jmr, lmr, imode
    !INTEGER, PARAMETER                   :: amonth = 2


    ! parameters for online dust calculations
    INTEGER, PARAMETER              :: ntraced=8                     ! number of coarse-grained bins 
                                                                     ! in the original emission model
    INTEGER, PARAMETER              :: nbin=24                       ! number of discretization points per bin
    INTEGER, PARAMETER              :: nclass=ntraced*nbin           ! total number of discretization points
    INTEGER, PARAMETER              :: nats=12                       ! number of soil types
    INTEGER, PARAMETER              :: nmode=4                       ! number of particle size distributions in soils,
                                                                     ! which distinguishes between clay, silt, 
                                                                     ! medium/fine sand, and coarse sand
    INTEGER, PARAMETER              :: nspe=nmode*3+2                ! for explanation, see below

    ! Constants used in the parameterization of the efficient friction velocity ratio,
    ! see Eqs. (17-20) in MB95: 
    REAL, PARAMETER                 :: aeff=0.35                     
    REAL, PARAMETER                 :: xeff=10.                     
    !
    !>>> TvN
    ! u1fac is a tuning parameter necessary to obtain 
    ! a reasonable global annual emission amount.
    ! u1fac < 1 is used to reduce the threshold friction velocity.
    ! In ECHAM-HAM simulations at T63 values of
    ! 0.86 and 0.56 were used by Cheng et al. (ACP, 2008).
    ! The lower value was introduced to increase emissions
    ! when surface roughness lengths were increased 
    ! from a constant value of 0.001 cm to values
    ! based on satellite measurements from Prigent et al. (JGR, 2005). 
    ! It is unclear where the value 0.66 specified below
    ! is based on.
    ! In ECHAM-HAM2 (Zhang et al., ACP, 2012)
    ! the satellite based surface roughness values
    ! were abandoned again.
    !REAL, PARAMETER                 :: u1fac=0.66                    ! scaling factor for threshold friction velocity
    ! Testing different values:
    !REAL, PARAMETER                 :: u1fac=0.56 
    REAL, PARAMETER                 :: u1fac=0.6
    !<<< TvN
    
    REAL, PARAMETER                 :: cd=1.2507E-06                 ! flux dimensioning parameter [g s^2/cm^4]
    
    !<<< TvN                                                                 ! (=roa/(grav*1.e2))
    ! ustar_min is not used:
    !REAL, PARAMETER                 :: ustar_min=5.                  ! min. fricton velocity (cm/s)
    ! minimum surface roughness length z0 (cm)
    ! The minimum value in the data set 
    ! from Prigent et al. is 1e-3 cm.
    ! but that seems very low.
    ! For instance, the minimum value in the 
    ! measurements used in the regression
    ! in that study is 2.3e-3 cm.
    ! Also, at very low z0, volume scattering
    ! of the microwave radiation will take place 
    ! that can significantly decrease the radar 
    ! backscatter coefficient (p. 8).
    ! Furthermore, using 1e-3 cm leads to 
    ! an overestimation of AOD (compared to MODIS)
    ! in the areas concerned,
    ! in particular around the dust hot spots
    ! of the Sahara (using current u1fac value).
    ! For these reasons the minimum value
    ! has been increased. 
    !REAL, PARAMETER                 :: z0_min=1.e-3
    !REAL, PARAMETER                 :: z0_min=5.e-3
    REAL, PARAMETER                 :: z0_min=1.e-2
    !REAL, PARAMETER                 :: z0_min=2.e-2
    !<<< TvN

    REAL, PARAMETER                 :: lai_lim=0.25
    REAL, PARAMETER                 :: lai_lim2=0.5

    ! d_thrsld [cm^2.5] = 0.006/(ddust * grav*1.e2) with ddust = 2.65 g/cm^3,
    ! see Eq. (4) in MB95:
    REAL, PARAMETER                 :: d_thrsld=2.31e-6           ! threshold value
    !>>> TvN
    ! There are eight coarse-grained size bins,
    ! of which only the first four are used here.
    ! According to Tegen et al., Heinold et al.,
    ! the radius boundaries of the first seven bins are 
    ! at 0.1, 0.3, 0.9, 2.6, 8.0, 24, 72, and 220 um.
    ! However, these number don't seem to be exact.
    ! Since there is a constant ratio between the right
    ! and low boundaries, it seems this ratio is 3.0.
    ! Indeed, in Laurent et al. (JGR, 2010), 
    ! 2.6 is corrected to 2.7, which would be consistent
    ! with 8.0/3.0 = 2.67.
    ! This would imply that the radius boundaries are at
    ! 0.0987654 = 72./(3.^6), 0.296296, 0.889, 2.67, 8.0, 24, 72, 216, 
    ! and 648 um.

    ! Next, each bin is discretized with 24 size points,
    ! where d(n+1) = d(n) * exp(Dstep).
    ! Thus, Dstep = ln(3.)/24 = 0.04577551202.
    ! Dmin is the diameter of the first size point,
    ! given by 2* 72./(3.^6)) * exp(0.5*Dstep) = 0.20210403762 um.
    ! Similarly, the last size point is at a diameter
    ! 2* 648. * exp(-0.5*Dstep) = 1266.67434757 um.
    ! 
    ! With the original bin settings,
    ! the number of size points is 191 not 192 (=8*24).
    !
    !REAL, PARAMETER                 :: Dmin=0.00002                  ! minimum partic. diameter (cm)
    !REAL, PARAMETER                 :: Dmax=0.130                    ! maximum partic. diameter (cm)
    !REAL, PARAMETER                 :: Dstep=0.0460517018598807      ! diameter increment
    REAL, PARAMETER                 :: Dmin=2.0210403762e-5          ! diameter (cm) at first discretization point
    REAL, PARAMETER                 :: Dmax=0.126667434757           ! diameter (cm) at last discretization point
    REAL, PARAMETER                 :: Dstep=0.04577551202           ! diameter increment in log-space
    !<<< TvN

    ! Constants in the parameterization of the Reynolds number,
    ! see Eq. (5) in MB95:
    REAL, PARAMETER                 :: a_rnolds=1331.647             ! Reynolds constant
    REAL, PARAMETER                 :: b_rnolds=0.38194              ! Reynolds constant
    REAL, PARAMETER                 :: x_rnolds=1.561228             ! Reynolds constant
    !
    ! Air density has been made variable,
    ! to account for orographic effects.
    ! Previously, a global value for the 
    ! threshold friction velocity Uth was calculated.
    ! To keep its unit the same,
    ! roa is kept as a reference value,
    ! but its exact value is not important anymore.
    REAL, PARAMETER                 :: roa=0.001227                  ! reference air density (g/cm^3)
    REAL                            :: rho_air                       ! variable air density (g/cm^3)
    REAL, PARAMETER                 :: airfac=1./rgas*xmair*1.e-6    ! factor for rho_air
    REAL                            :: airdens_ratio, airdens_ratio2
    !<<< TvN
    REAL, PARAMETER                 :: umin=13.75                    ! minimum threshold friction velocity (cm/s)
    REAL, PARAMETER                 :: ZZ=1000.                      ! wind measurement height (cm)

    ! parameters for the grouping in 2 modes
    ! The code follows the ECHAM-HAM implementation
    ! of Stier et al. (JGR, 2005),
    ! where the emission distribution is
    ! fitted onto three log-normal modes
    ! corresponding to the accumulation, coarse and super-coarse mode.
    ! (see presentation E. Vignati, TM meeting, 6 June 2008).
    !
    ! According to Heinold et al., 
    ! the three largest dust bins
    ! are less important for long-range transport,
    ! so particles with radius larger than 24 um
    ! can safely be neglected.
    ! However, a substantial part of the emitted mass
    ! is carried by particles with a radius larger than 10 um
    ! (see Tegen et al., Table 5).
    !
    ! The amounts of mass emitted in the accumulation and coarse modes
    ! are calculated from the masses emitted in the bin model,
    ! using two size ranges:
    ! r1 from 0.0987654 to 0.296296 um, and 
    ! r2 from 0.296296 to 8.0 um. 
    !
    ! Boundaries for Acc. mode
    INTEGER, PARAMETER              :: min_ai=1
    INTEGER, PARAMETER              :: max_ai=1
    ! Boundaries for Coa. mode
    INTEGER, PARAMETER              :: min_ci=2
    INTEGER, PARAMETER              :: max_ci=4
    !
    ! These size ranges include only part of
    ! the mass in the accumulation and coarse modes.
    ! The corresponding mass fractions are given by 
    ! mf(rmin,rmax) = 0.5*(
    !                erf(ln(rmax/mmr)/(sqrt(2)*ln(sigma)))-
    !                erf(ln(rmin/mmr)/(sqrt(2)*ln(sigma))) ),
    ! where mmr is the mass median radius.
    ! Applying this formula, 
    ! we find the following numbers:
    ! mf_acc(0,0.0987654)=0.00219913
    ! mf_acc_r1=mf_acc(0.0987654,0.296296)=0.313758
    ! mf_acc_r2=mf_acc(0.296296,8.0)=0.684043
    ! mf_acc(0.296296,inf)=0.684043
    !
    ! mf_coa(0,0.296296)=0.00519991
    ! mf_coa_r1=mf_coa(0.0987654,0.296296)=0.00518309 
    ! mf_coa_r2=mf_coa(0.296296,8.0)=0.980634
    ! mf_coa(8.0,inf)=0.0141665
    !
    REAL, PARAMETER :: mf_acc_r1 = 0.313758
    REAL, PARAMETER :: mf_acc_r2 = 0.684043
    REAL, PARAMETER :: mf_coa_r1 = 0.00518309
    REAL, PARAMETER :: mf_coa_r2 = 0.980634
    !
    ! Most importantly, r1 contains only about 31.4% 
    ! of the mass in the accumulation mode!
    ! This implies that we cannot just put the emissions
    ! from r1 to the accumulation mode, 
    ! and those from r2 to the coarse mode!
    !
    ! Instead, the modal emissions are determined
    ! by the following system of linear equations:
    ! mf_acc_r1 * flux_ai + mf_coa_r1 * flux_ci = flux_r1
    ! mf_acc_r2 * flux_ai + mf_coa_r2 * flux_ci = flux_r2,
    ! which relates the mass emitted in the ranges r1 and r2
    ! to the mass emitted in the accumulation and coarse modes.
    ! The solution is expressed using 
    ! the following parameters:
    !
    REAL, PARAMETER :: ratio_coa = mf_coa_r1/mf_coa_r2
    REAL, PARAMETER :: ratio_acc = mf_acc_r2/mf_acc_r1
    REAL, PARAMETER :: denom_acc_inv = 1./(mf_acc_r1-ratio_coa*mf_acc_r2)
    REAL, PARAMETER :: denom_coa_inv = 1./(mf_coa_r2-ratio_acc*mf_coa_r1)
    REAL, PARAMETER :: mf_acc_r12_inv = 1./(mf_acc_r1+mf_acc_r2)
    REAL, PARAMETER :: mf_coa_r12_inv = 1./(mf_coa_r1+mf_coa_r2) 
    ! 
    ! Source mass median radius (cm)
    ! Stier et al. (2005) uses very similar numbers 
    ! for mass median radii, 
    ! but uses 0.37 um for the accumulation mode.
    ! Thus, it seems these numbers are not mass mean,
    ! but mass median radii.
    !
    ! The super-coarse mode has 
    ! a mass median radius of 15.0 and sigma=2.0,
    ! but is not included.
    !
    ! The AeroCom recommendation of Dentener et al. (ACP, 2006)
    ! is to use a number median radius
    ! of 0.65 um for the coarse mode,
    ! which corresponds to mass median radius of 2.75 um
    ! (the conversion factor is exp(3.0*ln(sigma)^2),
    ! see Zender, Particle Size Distributions:
    ! Theory and Application to Aerosols, Clouds, and Soils, 2002).
    ! 
    !REAL, PARAMETER                 :: mmr_ai=0.35E-4
    REAL, PARAMETER                 :: mmr_ai=0.37E-4
    REAL, PARAMETER                 :: mmr_ci=1.75E-4
    !<<< TvN

    !----------------------------------------------------------------
    ! SOIL CARACTERISTICS:
    ! ZOBLER texture classes 
    !----------------------------------------------------------------

    INTEGER :: jp

    ! solspe includes for each soil type (first dimension)
    ! the mass median diameter (cm) and standard deviation (see Table 1, MB95)
    ! and the relative contribution (Table 2, MP95) for the four size populations.
    ! The two additional entries describe the saltation efficiency alpha (cm^-1),
    ! and the residual moisture, which is currently not used.
    ! Efficiencies are calculated as averages over the four populations
    ! (as in Eq. (8) in Marticorena et al. (JGR, 1997),
    ! where 1e-7, 1e-6 and 1e-5 cm^-1 is used for coarse sand, 
    ! medium/fine sand and silt, respectively,
    ! and 1e-6 for clay for soils with clay fractions below 45%
    ! and 1e-7 for clay for soils with clay fractions above 45%.
    ! (Tegen et al.).
    DOUBLE PRECISION, DIMENSION(nats,nspe) :: solspe

    !--     soil type 1 : Coarse
    DATA (solspe(1,jp),jp=1,nspe)/  &
         0.0707, 2.,  0.43 ,      &
         0.0158, 2.,  0.4 ,       &
         0.0015, 2.,  0.17 ,      &
         0.0002 ,2.,  0. ,        &
         2.1E-06,   0.2/
    !--     soil type 2 : Medium
    DATA (solspe(2,jp),jp=1,nspe)/  &
         0.0707, 2.,  0. ,            &
         0.0158, 2.,  0.37 ,          &
         0.0015, 2.,  0.33 ,          &
         0.0002, 2.,  0.3 ,           &
         4.0e-6,    0.25/
    !--     soil type 3 : Fine
    DATA (solspe(3,jp),jp=1,nspe)/  &
         0.0707, 2.,  0. ,            &
         0.0158, 2.,  0. ,            &
         0.0015, 2.,  0.33 ,          &
         0.0002, 2.,  0.67 ,          &
         !>>> TvN
         ! 33% x 1e-5 + 67% x 1e-7 = 3.367e-6 cm^-1
         !1.E-07,   0.5/
         3.4e-6,   0.5/
         !<<< TvN
    !--     soil type 4 : Coarse Medium
    DATA (solspe(4,jp),jp=1,nspe)/  &
         0.0707, 2.,  0.1 ,           &
         0.0158, 2.,  0.5 ,           &
         0.0015, 2.,  0.2 ,           &
         0.0002, 2.,  0.2 ,           &
         2.7E-06,   0.23/
    !--     soil type 5 : Coarse Fine
    DATA (solspe(5,jp),jp=1,nspe)/  &
         0.0707, 2.,  0. ,            &
         0.0158, 2.,  0.5 ,           &
         0.0015, 2.,  0.12 ,          &
         0.0002, 2.,  0.38 ,          &
         !>>> TvN
         ! 50% x 1e-6 + 12% x 1e-5 + 38% x 1e-6 = 2.08e-6 cm^-1
         !2.8E-06,   0.25/
         2.1e-6,   0.25/
         !<<< TvN
    !--     soil type 6 : Medium Fine
    DATA (solspe(6,jp),jp=1,nspe)/  &
         0.0707, 2.,  0.   ,          &
         0.0158, 2.,  0.27 ,          &
         0.0015, 2.,  0.25 ,          &
         0.0002, 2.,  0.48 ,          &
         !>>> TvN
         ! 27% x 1e-6 + 25% x 1e-5 + 48% x 1e-7 = 2.818e-6 cm^-1
         !1e-07,   0.36/
         2.8e-6,   0.36/
         !<<< TvN
    !--     soil type 7 : Coarse, Medium, Fine
    DATA (solspe(7,jp),jp=1,nspe)/  &
         0.0707, 2.,  0.23 ,          &
         0.0158, 2.,  0.23 ,          &
         0.0015, 2.,  0.19 ,          &
         0.0002, 2.,  0.35 ,          &
         2.5E-06,  0.25/
    !--     soil type 8 : Organic
    DATA (solspe(8,jp),jp=1,nspe)/  &
         0.0707, 2.,  0.25 ,          &
         0.0158, 2.,  0.25 ,          &
         0.0015, 2.,  0.25 ,          &
         0.0002, 2.,  0.25 ,          &
         0.,   0.5/
    !--     soil type 9 : Ice
    DATA (solspe(9,jp),jp=1,nspe)/  &
         0.0707,  2.,  0.25 ,         &
         0.0158,  2.,  0.25 ,         &
         0.0015,  2.,  0.25 ,         &
         0.0002,  2.,  0.25 ,         &
         0.,       0.5/
    !--     soil type 10 : Potential Lakes (additional)
    !       GENERAL CASE
    DATA (solspe(10,jp),jp=1,nspe)/  &
         0.0707,  2.,  0. ,            &
         0.0158,  2.,  0. ,            &
         0.0015,  2.,  1. ,            &
         0.0002,  2.,  0. ,            &
         1.E-05,  0.25/
    !--     soil type 11 : Potential Lakes (clay)
    !       GENERAL CASE
    DATA (solspe(11,jp),jp=1,nspe)/  &
         0.0707,  2.,  0. ,            &
         0.0158,  2.,  0. ,            &
         0.0015,  2.,  0. ,            &
         0.0002,  2.,  1. ,            &
         1.E-05,  0.25/
    !--     soil type 12 : Potential Lakes Australia
    DATA (solspe(12,jp),jp=1,nspe)/  &
         0.0707,  2.,  0. ,            &
         0.0158,  2.,  0. ,            &
         0.0027,  2.,  1. ,            &
         0.0002,  2.,  0. ,            &
         1.E-05,  0.25/


    INTEGER, PARAMETER :: add_field = 0

    REAL    :: veget, lai_max, lai_avg, lai_cur, z0s, dpd, flux_diam, cultfac1, dlast
    REAL    :: aaa, bb, ccc, ff, feff, dbstart, dp, uthp, wind10m, ustar
    REAL    :: xk, ddd, ee, stotal, stotalV, fdp1, fdp2
    REAL    :: su, suV, su_loc, su_locV, xl, xm, xn, xnv
    REAL    :: flux_r1, flux_r2
    INTEGER :: istat, nd, nsi, nm, np, ns, region, var_id, sds_id
    INTEGER :: i, j, i_s1, i_s11, i_s111, idust, lai_flag, month, kk, iveg
    INTEGER :: kkk, kfirst, kkmin, nn
    INTEGER :: i01, j01, i02, j02
    INTEGER :: i1, j1, i2, j2, access_mode

    ! for newsrun
    CHARACTER(len=200)      :: dust_filename, var_name
    REAL, DIMENSION(:),     ALLOCATABLE :: su_class, su_classv, utest
    REAL, DIMENSION(:,:,:), ALLOCATABLE :: global_3d
    REAL, DIMENSION(:,:),   ALLOCATABLE :: global_2d

    ! saving the status of being called
    LOGICAL, SAVE :: initial = .TRUE.


    ! ------------ start ------------------------------

    ! only valid for online method
    IF( TRIM( emis_input_dust ) /= 'ONLINE' ) RETURN

    ! MPI tile bounds of 1x1
    CALL get_distgrid( dgrid(iglbsfc), I_STRT=i01, I_STOP=i02, J_STRT=j01, J_STOP=j02 )
    
    
    ! =========================== INIT
    IF( initial ) THEN 

       ! global fields, which have to be available throughout the whole
       ! online emission procedure
       ALLOCATE(      uth(       nclass ) )
       ALLOCATE(    srel ( nats, nclass ) )
       ALLOCATE(    srelV( nats, nclass ) )
       ALLOCATE( su_srelV( nats, nclass ) )

       ! gridded 1x1 fields from input file
       ALLOCATE( soil_type ( i01:i02, j01:j02          ) )
       ALLOCATE( pot_source( i01:i02, j01:j02          ) )
       ALLOCATE( cult      ( i01:i02, j01:j02          ) )
       ALLOCATE( area      ( i01:i02, j01:j02          ) )
       ALLOCATE( z0        ( i01:i02, j01:j02, nz0     ) ) 
       ALLOCATE( fpar      ( i01:i02, j01:j02, nfpar   ) )
       ALLOCATE( soilph    ( i01:i02, j01:j02, nsoilph ) )

       ! additional 1x1 fields 
       ALLOCATE( snowcover( i01:i02, j01:j02 ) )
       ALLOCATE( desert   ( i01:i02, j01:j02 ) )
       ALLOCATE( umin2    ( i01:i02, j01:j02 ) )
       ALLOCATE( alpha    ( i01:i02, j01:j02 ) )
       ALLOCATE( c_eff    ( i01:i02, j01:j02 ) )
       ALLOCATE( lai_eff  ( i01:i02, j01:j02 ) )


       ! only needed within "initial"
       ALLOCATE( su_class ( nclass ) )
       ALLOCATE( su_classv( nclass ) )
       ALLOCATE( utest    ( nats   ) )


       ! ---------------------------------------------
       ! read input file 
       ! ---------------------------------------------
       dust_filename = TRIM(emis_input_dir_dust)//'/onlinedust_4.nc'
    
       WRITE(gol,'("Opening dust emission input file for ONLINE option: ",a)') dust_filename; CALL goPr
       
       IF (isRoot) THEN
          ! open file as source for information on 1x1 grid
          CALL MDF_Open( TRIM(dust_fileName), MDF_NETCDF4, MDF_READ, dust_FileID, status )
          IF_NOTOK_RETURN(status=1)
       ELSE
          ALLOCATE( global_3d(1,1,1) )
       ENDIF
       
       ! --- surface roughness (3d)
       IF (isRoot) THEN
          ALLOCATE( global_3d( nlon360, nlat180, nz0) )
          CALL MDF_Inq_VarID ( dust_FileID,   'z0',  var_id, status ) ; IF_NOTOK_MDF()
          CALL MDF_Get_Var   ( dust_FileID, var_id, global_3d, status ) ; IF_NOTOK_MDF()
       ENDIF
       
       CALL SCATTER ( dgrid(iglbsfc), z0, global_3d, 0, status)
       IF_NOTOK_RETURN(status=1)
       if (isRoot) DEALLOCATE( global_3d )

       ! --- soilph (3d)
       IF (isRoot) THEN
          ALLOCATE( global_3d( nlon360, nlat180, nsoilph) )
          CALL MDF_Inq_VarID ( dust_FileID, 'soilph', var_id, status ) ; IF_NOTOK_MDF()
          CALL MDF_Get_Var   ( dust_FileID, var_id,  global_3d, status ) ; IF_NOTOK_MDF()
       ENDIF
       
       CALL SCATTER ( dgrid(iglbsfc), soilph, global_3d, 0, status)
       IF_NOTOK_RETURN(status=1)
       if (isRoot) DEALLOCATE( global_3d )
       
       ! --- fpar (3d)
       IF (isRoot) THEN
          ALLOCATE( global_3d( nlon360, nlat180, nfpar) )
          CALL MDF_Inq_VarID ( dust_FileID, 'fpar',  var_id, status ) ; IF_NOTOK_MDF()
          CALL MDF_Get_Var   ( dust_FileID, var_id, global_3d, status ) ; IF_NOTOK_MDF()
       ENDIF
       
       CALL SCATTER ( dgrid(iglbsfc), fpar, global_3d, 0, status)
       IF_NOTOK_RETURN(status=1)
       DEALLOCATE( global_3d )  ! done for all
       
       ! -- other fields are 2D
       IF (isRoot) THEN
          ALLOCATE( global_2d( nlon360, nlat180) )
       ELSE
          ALLOCATE( global_2d(1,1) )
       ENDIF
       
       ! soiltype (2d)
       IF (isRoot) THEN
          CALL MDF_Inq_VarID ( dust_FileID, 'soiltype', var_id,  status ) ; IF_NOTOK_MDF()
          CALL MDF_Get_Var   ( dust_FileID,  var_id,   global_2d, status ) ; IF_NOTOK_MDF()
       ENDIF
       CALL SCATTER ( dgrid(iglbsfc), soil_type, global_2d, 0, status)
       IF_NOTOK_RETURN(status=1)

       ! potsrc
       IF (isRoot) THEN
          CALL MDF_Inq_VarID ( dust_FileID, 'potsrc', var_id, status ) ; IF_NOTOK_MDF()
          CALL MDF_Get_Var   ( dust_FileID, var_id, global_2d, status ) ; IF_NOTOK_MDF()
       ENDIF
       CALL SCATTER ( dgrid(iglbsfc), pot_source, global_2d, 0, status)
       IF_NOTOK_RETURN(status=1)

       ! cult
       IF (isRoot) THEN
          CALL MDF_Inq_VarID ( dust_FileID, 'cult', var_id, status ) ; IF_NOTOK_MDF()
          CALL MDF_Get_Var   ( dust_FileID, var_id, global_2d, status ) ; IF_NOTOK_MDF()
       ENDIF
       CALL SCATTER ( dgrid(iglbsfc), cult, global_2d, 0, status)
       IF_NOTOK_RETURN(status=1)

       ! grid area
       IF (isRoot) THEN
          CALL MDF_Inq_VarID ( dust_FileID, 'area', var_id, status ) ; IF_NOTOK_MDF()
          CALL MDF_Get_Var   ( dust_FileID, var_id, global_2d, status ) ; IF_NOTOK_MDF()
       ENDIF
       CALL SCATTER ( dgrid(iglbsfc), area, global_2d, 0, status)
       IF_NOTOK_RETURN(status=1)
       
       DEALLOCATE( global_2d ) ! done for all

       IF (isRoot) THEN
          CALL MDF_Close( dust_FileID, status )
          IF_NOTOK_RETURN(status=1)
       endif
          
       WRITE(gol,'("Closing dust emission input file")'); CALL goPr 

       !---------------------------------------------------------------------------------------
       !        initializations
       !---------------------------------------------------------------------------------------
       uth      = 0. 
       srel     = 0.          ! fraction of the grid area correspondent to each soil population
       srelV    = 0.          ! fraction of volume
       su_srelV = 0.
       utest    = 0.

       !---------------------------------------------------------------------------------------
       !       Uth calculation 
       !       Threshold friction velocity dependent on the particle diameter
       !       following Eqs. (3-5) in MB95.
       !---------------------------------------------------------------------------------------
       nn = 0
       dp = Dmin
       DO WHILE( dp <= Dmax + 1E-05 )
          nn  = nn + 1
          BB  = a_rnolds * (dp ** x_rnolds) + b_rnolds
          XK  = SQRT(ddust * grav*100. * dp / roa)        ! grav should be in cm s-2
          CCC = SQRT(1. + d_thrsld /(dp ** 2.5))
          IF( BB < 10. ) THEN
             DDD = SQRT(1.928 * (BB ** 0.092) - 1.)
             Uth(nn) = 0.129 * XK * CCC / DDD
          ELSE
             EE = -0.0617 * (BB - 10.) 
             FF = 1. -0.0858 * EXP(EE)
             Uth(nn) = 0.12 * XK * CCC * FF       
          END IF
          dp = dp * EXP(Dstep)  
       END DO


       !---------------------------------------------------------------------------------------
       !       surface calculation - calculation of the soil size distribution
       !       Through all soil particle diameter the calculation of the relative contribution
       !       in surface and volume of the soil population independently of the grid 
       !---------------------------------------------------------------------------------------
       DO ns = 1, nats ! soil types

          Stotal    = 0.
          StotalV   = 0.
          su_class  = 0. 
          su_classV = 0. 

          kk = 0
          dp = Dmin
          DO WHILE( dp <= Dmax + 1.0E-5 )        
             kk  = kk + 1
             su  = 0.
             suV = 0.
             DO nm = 1, Nmode            ! particle size populations in soils
                nd  = ((nm - 1) *3 ) + 1 ! index to mass median diameter
                nsi = nd + 1             ! index to standard deviation
                np  = nd + 2             ! index to relative contribution
                !
                !   based on soil type and contribution of population of the soil type the soil size
                !   distribution population is calculated
                !

                !>>> TvN
                ! Bug in the original code: nd should be np
                ! Since solspe(ns,nd) is never zero
                ! and the final result is proportional to solspe(ns,np),
                ! the bug has no impact on the results. 
                !IF (solspe(ns,nd).EQ.0.) THEN            
                IF (solspe(ns,np).EQ.0.) THEN
                !<<< TvN
                   su_loc = 0.
                   su_locV=0.
                ELSE
                   xk      = solspe(ns,np)/(SQRT(2.* pi)*LOG(solspe(ns,nsi)))
                   xl      = ( (LOG(dp) - LOG( solspe(ns,nd ) ))**2 ) / &
                        (2.*(LOG( solspe(ns,nsi) ))**2 )
                   xm      =  xk * EXP(-xl)         ! value of the lognormal mass size distribution
                                                    ! dM/dln(dp) in Eq. (29) in MB95
                                                    ! (Aerosol Sci. Technol., 1994)
                   xn      =  ddust*(2./3.)*(dp/2.) ! surface
                                                    ! cf. the denominator in Eq. (30) in MB95
                                                    ! The factor 2 difference is irrelevant,
                                                    ! since only relative contributions are used.
                   xnV     =  1. !volume
                   su_loc  = (xm*Dstep/xn)          ! Eq. (30) in MB95
                   su_locV = (xm*Dstep/xnV)        
                END IF !
                su  = su  + su_loc
                suV = suV + su_locV
             END DO !Nmode

             su_class(kk)   = su
             su_classV(kk)  = suV
             Stotal         = Stotal + su
             StotalV        = StotalV + suV
             dp             = dp * EXP(Dstep)
          END DO !dp

          DO nn = 1,Nclass
             IF (Stotal.EQ.0.)THEN
                srel (ns,nn) = 0.
                srelV(ns,nn) = 0.
             ELSE
                srel    (ns,nn) = su_class(nn)/Stotal
                srelV   (ns,nn) = su_classV(nn)/StotalV
                utest   (ns   ) = utest(ns)+srelV(ns,nn)
                su_srelV(ns,nn) = utest(ns) 
             END IF
          END DO !j=1,nclass
       END DO !ns (soil type)


       DEALLOCATE( su_class, su_classV, utest )

       ! reset initial 
       initial = .FALSE.

    END IF ! =========================== INIT

    ! only once per day 
    IF( newday ) THEN
       
       ! calculation of snow cover from snow dept
       ! Tegen et al. fraction (0-1)
       snowcover = sd_dat(iglbsfc)%data(:,:,1) / 0.015
       WHERE( snowcover > 1. ) snowcover = 1.

       !---------------------------------------------------------------------------------------
       !       Prepare the flux calculation
       !---------------------------------------------------------------------------------------
       !
       !       Calculations done on monthly fields

       ! default: no dust source due to 
       !          - vegetation
       !          - not a desert pixel or 
       !          - no pure land grid cell
       lai_eff = 0. 

       ! per grid box
       DO j=j01, j02
          DO i=i01, i02

             !---------------------------------------------------------------------------------------
             !       Selection of potential dust sources areas
             !---------------------------------------------------------------------------------------
             !      Preferential Sources = Potential lakes
 
             !>>> TvN
             ! If monthly surface roughness is not available
             ! use the annual mean value, if available.
             ! Since the annual mean is calculated
             ! based on all available months,
             ! it has a much better spatial coverage 
             ! than the individual months.
             IF( z0(i,j,idate(2)) <= 0. .AND. z0(i,j,nz0) > 0. ) THEN 
               z0(i,j,idate(2)) = z0(i,j,nz0) 
             ENDIF
             !<<< TvN
        
             IF( pot_source(i,j) > 0.5 ) THEN 
                ! if the potential lake area is > 50%, it is a pot. lake grid
                soil_type(i,j) = 10.                 
                !>>> TvN
                ! Use minimum value for roughness length.
                ! Since there are only few potential source areas
                ! where the annual mean is not available,
                ! this will only have a limited impact.
                !IF( z0(i,j,idate(2)) <= 0. ) z0(i,j,idate(2)) = 0.001 !! if z0 is not valid or missing (cm), PhD thesis Marticorena p.85
                IF( z0(i,j,idate(2)) <= 0. ) z0(i,j,idate(2)) = z0_min
                !<<< TvN
             END IF

             !---------------------------------------------------------------------------------------
             !       Calculation of the ratio: horizontal/vertical flux (alpha)
             !---------------------------------------------------------------------------------------
             !---------------------------------------------------------------------------------------
             !       Test on the vegetation type
             !---------------------------------------------------------------------------------------
             !  When cult=0, the cultivation field info is not used. Otherwise: cult(i,j)=3
!!$             cult(i,j)   = 0.

             desert(i,j) = soilph(i,j,3) + soilph(i,j,4)
             veget=0.
             DO iveg=1,nveg
                veget = veget + tv_dat(iglbsfc,iveg)%data(i,j,1)
             END DO

             ! default: no dust emissions
             idust = 0 
             ! dust emissions only when 
             ! 1) there is only land (almost)
             ! 2) 'desert' is positive or vegetation active
             IF( lsmask_dat(iglbsfc)%data(i,j,1) >= 99. .AND. (desert(i,j) > 0. .OR. veget > TINY(veget)) ) &
                  idust = 1

             ! here is dust uptake possible
             IF( idust == 1 ) THEN

                !---------------------------------------------------------------------------------------
                !--  Calculate effective surface for fpar < lai_lim (as proxy for
                !--  veg. cover), shrubby vegetation is determined by max
                !--  annual fpar, grassy by monthly fpar (Tegen et al.2002)
                !---------------------------------------------------------------------------------------

                ! so we start with no vegetation --> full area available
                lai_eff(i,j) = 1. 


                !--    get max/mean fpar of the full year --> needed for shrub land
                lai_max = MAXVAL(fpar(i,j,1:12))
                lai_avg =    SUM(fpar(i,j,1:12)) / 12. 
                lai_cur = fpar(i,j,idate(2))


                ! ---------------------------------------------
                ! 3 classes: grass, shrub, mixed{grass,shrub}
                ! ---------------------------------------------

                ! first: grass dominated (tv(2) and tv(7))
                !        current fpar determines available area
                IF( (tv_dat(iglbsfc,2)%data(i,j,1) + tv_dat(iglbsfc,7)%data(i,j,1)) > 50 ) THEN 

                   lai_eff  (i,j) = 1. - lai_cur / lai_lim

                   ! second: shrub dominated (tv(16) and tv(17))
                   !         if max(fpar) > 0.25 --> no dust 
                   !         else max(fpar) determines area
                ELSEIF( (tv_dat(iglbsfc,16)%data(i,j,1) + tv_dat(iglbsfc,17)%data(i,j,1)) > 50 ) THEN 

                   ! lai_eff is zero for lai_max > lai_min and 
                   ! [0,1] for lai_max < lai_lim
                   lai_eff  (i,j) = 1. - lai_max / lai_lim

                   ! third: mixtures of grass and shrub land
                   !        if mean(fpar) > 0.5 --> shrub dominated --> use max(fpar) for scaling
                   !        else grass dominated --> use current(fpar) for scaling
                ELSE

                   IF( lai_avg > lai_lim2 ) THEN 
                      lai_eff  (i,j) = 1. - lai_max / lai_lim
                   ELSE
                      lai_eff  (i,j) = 1. - lai_cur / lai_lim
                   END IF

                END IF

                ! limit to valid range [0,1]
                lai_eff(i,j) = MAX( 0., MIN( 1., lai_eff(i,j) ) )


!!$                ............... !!!!hier ist das AND falsch!!!! ..................
!!$                DO month = 1, 12
!!$                   IF(     ( tv_dat(i,j,16) > 50 ) .OR. ( tv_dat(i,j,17) > 50 ) .AND. ( lai_flag == 0 ) ) then
!!$                      lai_eff(i,j,month) = 1. - fpar(i,j,month) / lai_lim  !lai_lim=0.25
!!$                   ELSEIF( ( tv_dat(i,j, 2) > 50 ) .OR. ( tv_dat(i,j, 7) > 50 ) .OR.  ( desert(i,j) > 0.) ) then 
!!$                      lai_eff(i,j,month) = 1. - fpar(i,j,month) / lai_lim  !lai_lim=0.25
!!$                   ELSE
!!$                      lai_eff(i,j,month) = 1. - lai_max / lai_lim       !shrubs=1
!!$                   END IF
!!$                   ! original formulation
!!$                   !                 lai_eff(j,i,1,month)=1.-(lai(j,i,1,month)     &
!!$                   !                 *(1.-shrub(int(sp(j,i,1,5))))                   &
!!$                   !                 +lai_max*shrub(int(sp(j,i,1,5)))           &
!!$                   !                        )*1./lai_lim
!!$                   IF( lai_eff(i,j,month) <= 0 ) lai_eff(i,j,month) = 0
!!$                   IF( lai_eff(i,j,month) >  1 ) lai_eff(i,j,month) = 1
!!$                END DO  !month

             END IF    ! if idust=1

             !             print *, 'lai_eff=1 everywhere'
             !---------------------------------------------------------------------------------------
             !     Lowering the threshold friction velocity depending on the presence of cultivations
             !---------------------------------------------------------------------------------------
             !       Factors according to dsf increase seen in data **
             !---------------------------------------------------------------------------------------
             umin2(i,j) = umin
             ! 
             !---------------------------------------------------------------------------------------
             IF( cult(i,j) <= 0.5 .AND. cult(i,j) > 0.08 ) THEN
                IF( desert(i,j) > 0. .OR. tv_dat(iglbsfc,16)%data(i,j,1) > 50 .OR. tv_dat(iglbsfc,17)%data(i,j,1) > 50 ) & 
                     umin2(i,j) = umin * 0.93
                ! 
                !---------------------------------------------------------------------------------------
                IF( tv_dat(iglbsfc,2)%data(i,j,1) > 50 .OR. tv_dat(iglbsfc,7)%data(i,j,1) > 50 ) & 
                     umin2(i,j) = umin * 0.99
             END IF !cult=2

             !  
             !---------------------------------------------------------------------------------------
             IF( cult(i,j) > 0.5 ) THEN
                IF( ( desert(i,j) > 0 ) .OR. ( tv_dat(iglbsfc,16)%data(i,j,1) > 50 ) .OR. ( tv_dat(iglbsfc,17)%data(i,j,1) > 50 ) ) &
                     umin2(i,j) = umin * 0.73                         
             END IF !cult=1

             !---------------------------------------------------------------------------------------
             !       Daily z0 and efficient fraction feff
             !---------------------------------------------------------------------------------------

             i_s1 = INT( soil_type(i,j) )             ! soil type index for the calcl. of horiz. dust flux
             IF( i_s1 == 0 ) i_s1 = 9                 ! set it the same as ice if the soil type is not defined

             ! Roughness length [cm] of the surface without obstacles, i.e. of the smooth surface:
             Z0S  = 0.001 !! en cm, these Marticorena p.85    ! optimum value for the calculation of energy loss

              
             ! Soil-type dependent saltation efficiency,
             ! i.e. the ratio between vertical and horizontal fluxes,
             ! (see  Eq. (42) in MB95; Eq. (3) in Heinold et al.):
             alpha(i,j) = solspe(i_s1,nmode*3+1)

             ! for now moist is not included but when it is done then:

             !---------------------------------------------------------------------------------------
             !       Calculation of the threshold soil moisture (w')  [Fecan, F. et al., 1999] 
             !---------------------------------------------------------------------------------------
             !          when moist is included   !!!!!!!!!!!!!!!!!!
             !          w_str(j,i,1) = 0.0014*(solspe(i_s1,nmode*3)*100)**2 + 0.17*(solspe(i_s1,nmode*3)*100)
             !          W0   = 0.99           ! used by Bernd solspe(i_s1,nmode*3+2)
             feff = 0.
             !          * partition of energy between the surface and the elements of rugosity *
             !               these pp 111-112

             IF( z0(i,j,idate(2) ) <= 0. ) THEN   ! if there are no info on z0 and no potential sources
                z0(i,j,idate(2)) = 1.             ! then z0 is set to 1 and no dust can be produced
                feff = 0.
             ELSE
                !>>> TvN
                ! Use minimum value for roughness length.
                z0(i,j,idate(2)) = max(z0_min,z0(i,j,idate(2)) )
                !<<< TvN
                ! Eq. (20) in MB95:
                AAA = LOG( z0(i,j,idate(2)) / Z0S )
                BB  = LOG( aeff * (xeff / Z0S)**0.8)
                CCC = 1. - AAA/BB
                !          * partition between Z01 and Z02 * which are z0 of larger stone which cannot be mobilized
                FF = 1.    ! we do not separate roughness length between soil which
                           ! gives dust and solid material which is not mobilised
                ! total efficient friction velocity ratio:
                feff = FF * CCC
                ! restrict to [0,1]
                feff = MIN( 1., feff )
                feff = MAX( 0., feff )
             END IF

             c_eff(i,j) = feff  ! scaling parameter for the threshold friction velocity

             ! due to energy loss
             !---------------------------------------------------------------------------------------
          END DO   ! j
       END DO     ! i
       !---------------------------------------------------------------------------------------
       !      End of daily base calculations

    END IF ! newday 


    ! additional fields
    ALLOCATE( fluxtyp  (nclass  ) )
    !ALLOCATE( dpk      (nclass  ) )
    !ALLOCATE( dbmin    (ntraced ) )
    !ALLOCATE( dbmax    (ntraced ) )
    ALLOCATE( fluxtot  (ntraced ) )
    ALLOCATE( fdust    (ntraced ) )

    ! 1x1 flux mass and numbers
    ALLOCATE( fnum_ai  ( i01:i02,j01:j02 ) )
    ALLOCATE( flux_ai  ( i01:i02,j01:j02 ) )
    ALLOCATE( fnum_ci  ( i01:i02,j01:j02 ) )
    ALLOCATE( flux_ci  ( i01:i02,j01:j02 ) )

    ! reset flux masses 
    flux_ai = 0. 
    flux_ci = 0. 

!!$! start omp sharing work
!!$!$OMP PARALLEL default( private ) &
!!$!$OMP          shared ( c_eff, z0, umin2, uth, srel, srelV, su_srelV, alpha, snowcover,   & 
!!$!$OMP                   lai_eff, flux_ai, flux_ci, fnum_ai, fnum_ci, idate, area,         &
!!$!$OMP                   soil_type, lsmask_dat, sec_month, sigma, u10m_dat, v10m_dat,      &
!!$!$OMP                   spm_dat, t2m_dat )
!!$!$OMP DO 
    DO j = j01, j02
       DO i = i01, i02

          !-- initialisation of the fields
          !   size: ntraced
          fluxtot = 0. 
          fdust   = 0.

          !----- --------------------------------------------------------------------------
          !     Calculation of dust emission flux
          !     dependent on the 3 hourly wind fields
          !----------------------------------------------------------------------
          IF( c_eff(i,j) > 0. ) THEN 

             ! Calculation of ustar

             ! AS: initialise ustar (for those cases where if statement(s) are not fulfilled)
             ustar = 0. 

             IF( lsmask_dat(iglbsfc)%data(i,j,1) > 0. ) THEN 
                wind10m = SQRT(u10m_dat(iglbsfc)%data(i,j,1)**2 + &
                     v10m_dat(iglbsfc)%data(i,j,1)**2) * 100. ! cm/s
                ustar = (vKarman * wind10m) / ( alog( ZZ / z0(i,j,idate(2)) ) ) ! cm/s
             ENDIF

             IF( Ustar > 0 .AND. (Ustar > umin2(i,j) / c_eff(i,j)) ) THEN

                !>>> TvN 
                rho_air = spm_dat(iglbsfc)%data(i,j,1)/t2m_dat(iglbsfc)%data(i,j,1)*airfac ! g/cm3
                airdens_ratio  = rho_air/roa
                airdens_ratio2 = sqrt(roa/rho_air)
                !<<< TvN

                !-- initialisation of the fields
                !   size: ntraced
                !dbmin   = 0. 
                !dbmax   = 0. 
                !   size: nclass
                fluxtyp = 0.


                ! soil type index for the calcl. of horiz. dust flux
                i_s1 = INT( soil_type(i,j) )                
                ! set it the same as ice
                IF( i_s1 == 0 ) i_s1 = 9                    
                ! to separate from now on between saltation and mobilisation
                i_s11 = i_s1                                  
                ! to separate between mobilisation and saltation and dust particles
                IF( i_s1 == 10 .OR. i_s1 == 12 ) i_s11 = 11 

                kk = 0
                dp = Dmin
                DO WHILE( dp <= Dmax+1E-5)
                   kk    = kk+1
                   uthp  = uth(kk) * umin2(i,j) / umin * u1fac !reduce saltation threshold for cultivated soils
                   !>>> TvN
                   ! Include correction factor for variable air density
                   uthp = uthp * airdens_ratio2
                   !<<< TvN

                   ! See Eq. (28) in MB95; Eq. (6) in Tegen et al.; Eq. (2) in Heinold et al.
                   ! Note that (1+R)^2 * (1-R) = (1+R) * (1-R^2)
                   fdp1 = (1.-(Uthp/(c_eff(i,j) * Ustar)))          ! component of the horiz. flux
                   fdp2 = (1.+(Uthp/(c_eff(i,j) * Ustar)))**2.      !    

                   IF( fdp1 > 0 .AND. fdp2 > 0) THEN

                      ! vertical flux dust weighted by the surface area relative to each soil type
                      flux_diam = srel(i_s1,kk) * fdp1 * fdp2 * cd * Ustar**3 * alpha(i,j)
                      !>>> TvN
                      ! Include correction factor for variable air density
                      flux_diam = flux_diam * airdens_ratio
                      !<<< TvN

                      !----------------------------------------------------------------------
                      !   all particles even the small ones can be mobilised by saltation
                      !----------------------------------------------------------------------
                      dbstart = dmin

                      IF( dbstart >= dp ) THEN 
                         fluxtyp(kk) = fluxtyp(kk) + flux_diam
                      ELSE
                         !----------------------------------------------------------------------
                         !  loop over dislocated dust particle sizes
                         !----------------------------------------------------------------------
                         dpd    = dmin
                         kkk    = 0
                         kfirst = 0
                         DO WHILE( dpd <= dp+1e-5 )
                            kkk = kkk + 1
                            IF( dpd >= dbstart ) THEN                   ! the particles produced by saltation are put
                               IF( kfirst == 0 ) kkmin = kkk               ! in finer bins
                               kfirst = 1
                               !----------------------------------------------------------------------
                               !  scaling with relative contribution of dust size  fraction
                               !  we take into account the volume contribution of the particle types:
                               !  all the particles from soil type 10 are put into the 11 soil type when
                               !  we are in the production region
                               !----------------------------------------------------------------------
                               IF( kk > kkmin ) THEN
                                  ! remember: i_s11 puts the mobilised
                                  fluxtyp(kkk) = fluxtyp(kkk) + flux_diam * srelV(i_s11,kkk) / &
                                       (su_srelV(i_s11,kk) - su_srelV(i_s11,kkmin) )
                                  ! particles in smaller bins
                               END IF !kk.gt.kmin
                            END IF !dpd.gt.dbstart
                            dpd = dpd * EXP(dstep)
                         END DO !dpd
                         !----------------------------------------------------------------------
                         !  end of saltation loop
                         !----------------------------------------------------------------------
                      END IF !dbstart.lt.dp

                   END IF !fdp1

                   dp = dp * EXP(Dstep)

                END DO !dp   
                !----------------------------------------------------------------------
                !  assign fluxes to bins: flux is in g cm-2 s-1 for each bin
                !  192 sub-bins are put into 8 bins
                !----------------------------------------------------------------------
                dp    = dmin   
                dlast = dmin
                nn    = 1
                kk    = 0
                DO WHILE( dp <= dmax+1e-5 )  
                   kk = kk+1
                   ! add to total
                   IF( nn <= ntraced ) fluxtot(nn) = fluxtot(nn) + fluxtyp(kk) 

                   IF( MOD(kk,nbin) == 0 ) THEN
                      !dbmax(nn) = dp * 10000. * 0.5  !radius in um
                      !dbmin(nn) = dlast * 10000. * 0.5
                      !dpk(nn)   = SQRT( dbmax(nn) * dbmin(nn) )
                      nn        = nn+1
                      dlast     = dp
                   END IF

                   dp = dp * EXP(Dstep)

                END DO !dp      

             END IF   !ustar
          END IF   !c_eff 

          ! Masking the area covered by snow, vegetation and [...?...]
          cultfac1 = 1.

          DO nn = 1, ntraced
             !          fluxtot: g/cm2/sec 
             !    MASK: Effective area determined by cultfac1/snow
             fdust(nn) = fluxtot(nn) * cultfac1 * (1.-snowcover(i,j))

             !    MASK: Effective area determined by fpar:

             fdust(nn) = fdust(nn) * lai_eff(i,j) ! turn off vegetation limitation here!
             ! TvN: an alternative approach based on surface roughness
             ! is applied by Laurent et al. (JGR, 2006).


             !    MASK: Soil moisture threshold, using w0
             !          when moisture is included   !!!!!!!!!!!!!!!!!!
             !             IF(qrsur(i,j).GE.w0) THEN
             !               fdust(i,j,nn)=0.
             !             END IF
          END DO

          ! ------------------------------------------------------------------------------
          ! Grouping into 2 modes: 1sec accumulation
          !
          !>>> TvN
          !   Accumulation
          flux_r1 = 0.
          DO nn = min_ai, max_ai
           !flux_ai(i,j) = flux_ai(i,j) + fdust(nn)
           flux_r1 = flux_r1 + fdust(nn)
          END DO

          !   Coarse
          flux_r2 = 0.
          DO nn = min_ci, max_ci
             !flux_ci(i,j) = flux_ci(i,j) + fdust(nn)
             flux_r2 = flux_r2 + fdust(nn)
          END DO

          ! The solution of the system of linear equations
          ! (see comments above).
          ! For special conditions, 
          ! the solution can give a negative mass flux 
          ! in either the accumulation or coarse mode.
          ! In those case, all mass is put into
          ! the other mode.
          flux_ai(i,j) = flux_r1 - ratio_coa * flux_r2
          flux_ci(i,j) = flux_r2 - ratio_acc * flux_r1
          IF (flux_ai(i,j) .gt. 0. .AND. flux_ci(i,j) .gt. 0.) THEN
            flux_ai(i,j) = flux_ai(i,j) * denom_acc_inv
            flux_ci(i,j) = flux_ci(i,j) * denom_coa_inv
          ELSEIF (flux_ai(i,j) .lt. 0.) THEN
            flux_ai(i,j) = 0.
            flux_ci(i,j) = (flux_r1 + flux_r2) * mf_coa_r12_inv
          ELSEIF (flux_ci(i,j) .lt. 0.) THEN
            flux_ai(i,j) = (flux_r1 + flux_r2) * mf_acc_r12_inv
            flux_ci(i,j) = 0.
          ENDIF
          !<<< TvN

          ! now scale the emissions
          ! convert from g/cm2/s to g/grid_cell/month (area is in m2)
          flux_ai(i,j) = flux_ai(i,j) * sec_month * 1.E4 * area(i,j)
          flux_ci(i,j) = flux_ci(i,j) * sec_month * 1.E4 * area(i,j)

          !-------------------------------------------------------------------------------
          !  Calculating number flux (#/grid_cell/month)
          !
          !   Accumulation
          fnum_ai(i,j) = flux_ai(i,j) * 3. / (4.*pi*ddust*mmr_ai**3) * EXP(4.5*LOG(sigma(iacci))**2)
          !   Coarse
          fnum_ci(i,j) = flux_ci(i,j) * 3. / (4.*pi*ddust*mmr_ci**3) * EXP(4.5*LOG(sigma(icoai))**2)


          ! scale the flux from g to kg
          flux_ai(i,j) = flux_ai(i,j) * 1.E-03
          flux_ci(i,j) = flux_ci(i,j) * 1.E-03

       END DO ! j
    END DO ! i

!!$!$OMP END DO 
!!$!$OMP END PARALLEL
!!$! end omp sharing work 

    ! free memory
    !DEALLOCATE( fluxtyp, dpk, dbmin, dbmax, fluxtot, fdust )
    DEALLOCATE( fluxtyp, fluxtot, fdust )


    ! ------------------------------
    ! add sources to emission arrays
    ! ------------------------------

    ! vertical splitting is that class
    splittype = 'nearsurface'

    ! work arrays for gather-coarsen-scatter (skip if region #1 is at 1x1,
    ! assuming that no zoom region then)
    IF ( (iglbsfc /= 1) .OR. okdebug) THEN
       IF (isRoot) THEN
          ALLOCATE(glb_arr(nlon360,nlat180))
          DO region = 1, nregions
             ALLOCATE(emis_glb(region)%surf(im(region), jm(region)))
          END DO
       ELSE
          ALLOCATE(glb_arr(1,1))
          DO region = 1, nregions
             ALLOCATE(emis_glb(region)%surf(1,1))
          END DO
       END IF
    END IF

    ! work array for vertical distribution
    DO region = 1, nregions
       CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )          
       lmr = lm(region)
       ALLOCATE( emis3d(region)%d3(i1:i2, j1:j2, lm(region)) )
    END DO

    ! ------------------------------
    ! accumulation mode
    
    ! number
    CALL fill_target_array( fnum_ai, 'number acc', 'dust_number' )  ! fill emis_temp(region)
    IF_NOTOK_RETURN(status=1)

    ! vertically distribute according to sector
    DO region = 1, nregions
       emis3d(region)%d3 = 0.0
       CALL emission_vdist_by_sector( splittype, 'DUST', region, emis_temp(region), emis3d(region), status )
       emis_number(region,mode_aci)%d4(:,:,:,1)   = emis3d(region)%d3
    ENDDO

    ! mass
    CALL fill_target_array( flux_ai, 'mass acc', 'dust_mass' )  ! fill emis_temp(region)
    IF_NOTOK_RETURN(status=1)

    ! vertically distribute according to sector
    DO region = 1, nregions
       emis3d(region)%d3 = 0.0
       CALL emission_vdist_by_sector( splittype, 'DUST', region, emis_temp(region), emis3d(region), status )
       emis_mass(region,mode_aci)%d4(:,:,:,1)   = emis3d(region)%d3
    ENDDO

    ! ------------------------------
    ! coarse mode

    ! number
    CALL fill_target_array( fnum_ci, 'number coa', 'dust_number' )  ! fill emis_temp(region)
    IF_NOTOK_RETURN(status=1)

    ! vertically distribute according to sector
    DO region = 1, nregions
       emis3d(region)%d3 = 0.0
       CALL emission_vdist_by_sector( splittype, 'DUST', region, emis_temp(region), emis3d(region), status )
       emis_number(region,mode_coi)%d4(:,:,:,1)   = emis3d(region)%d3
    ENDDO

    ! mass
    CALL fill_target_array( flux_ci, 'mass coa', 'dust_mass' )  ! fill emis_temp(region)
    IF_NOTOK_RETURN(status=1)

    ! vertically distribute according to sector
    DO region = 1, nregions
       emis3d(region)%d3 = 0.0
       CALL emission_vdist_by_sector( splittype, 'DUST', region, emis_temp(region), emis3d(region), status )
       emis_mass(region,mode_coi)%d4(:,:,:,1)   = emis3d(region)%d3
    ENDDO


    ! free memory
    DEALLOCATE( fnum_ai, flux_ai, fnum_ci, flux_ci )
    DO region = 1, nregions
       IF (ASSOCIATED(emis_glb(region)%surf)) DEALLOCATE(emis_glb(region)%surf)
       DEALLOCATE(emis3d(region)%d3)
    END DO
    IF (ALLOCATED(glb_arr)) DEALLOCATE(glb_arr)

  CONTAINS

    !--------------------------------------------------------------------------
    !                    TM5                                                  !
    !--------------------------------------------------------------------------
    !BOP
    !
    ! !IROUTINE:    FILL_TARGET_ARRAY
    !
    ! !DESCRIPTION: Convenient internal program to fill EMIS_TEMP - same as routine in emission_ss.F90    
    !\\
    !\\
    ! !INTERFACE:
    !
    SUBROUTINE fill_target_array( arr_in, str1, str2 )
      !
      ! !INPUT PARAMETERS:
      !
      REAL,             INTENT(in) :: arr_in(i01:,j01:)
      CHARACTER(len=*), INTENT(in) :: str1, str2
      !
      ! !REVISION HISTORY: 
      !   25 Jun 2012 - P. Le Sager - v0
      !
      !EOP
      !------------------------------------------------------------------------
      !BOC

      ! Test for optimization: if region #1 is at 1x1, assume no zoom region
      ! and skip call to coarsen and mpi comm 
      IF (iglbsfc == 1) THEN

         emis_temp(1)%surf = arr_in

         IF (okdebug) THEN
            CALL gather(dgrid(iglbsfc), arr_in, glb_arr, 0, status)
            IF_NOTOK_RETURN(status=1)
!FIXME - really needed? too much messaging
!             IF (isRoot) THEN
!                CALL msg_emis( amonth, ' ', TRIM(str1), 'DU', 1., SUM(glb_arr) )
!             END IF
         END IF

      ELSE

         DO region = 1, nregions
            emis_temp(region)%surf = 0.0
         END DO

         ! gather on 1x1, coarsen to other grid on root, then scatter
         !-----------------------------------------------------------
         ! FIXME-PERF : Need a coarsen that works independtly on each proc, regardless of zooming... :(    
         CALL gather(dgrid(iglbsfc), arr_in, glb_arr, 0, status)
         IF_NOTOK_RETURN(status=1)

         IF (isRoot) THEN

!FIXME - really needed? too much messaging
!            CALL msg_emis( amonth, ' ', TRIM(str1), 'DU', 1., SUM(glb_arr) )
            CALL coarsen_emission(TRIM(str2), nlon360, nlat180, glb_arr, emis_glb, add_field, status)
            IF_NOTOK_RETURN(status=1)

         END IF

         DO region = 1, nregions
            CALL scatter(dgrid(region), emis_temp(region)%surf, emis_glb(region)%surf, 0, status)
            IF_NOTOK_RETURN(status=1)
         ENDDO

      ENDIF

      status=0

    END SUBROUTINE FILL_TARGET_ARRAY
    !EOC

  END SUBROUTINE CALC_EMISSION_DUST
  !EOC



  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_DUST_DECLARE
  !
  ! !DESCRIPTION: Open the input file(s), if appropriate. 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_DUST_DECLARE( status )
    !
    ! !USES:
    !
    USE io_hdf,        ONLY : DFACC_READ
    USE dims,          ONLY : idate
    !
    ! !OUTPUT PARAMETERS:
    !
    INTEGER, INTENT(out)        ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - added TEGEN-VIGNATI and ONLINE options
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    CHARACTER(len=*), PARAMETER :: rname = mname//'/emission_dust_declare'
    CHARACTER(len=256)          :: dust_filename
    INTEGER                     :: stat, sfstart, sfend

    ! --- begin -----------------------------------------

    
    SELECT CASE( TRIM( emis_input_dust ) )
    CASE( "AEROCOM" )

       WRITE (dust_filename,'(a,"/dust",i4.4,i2.2,".nc")') TRIM(emis_input_dir_dust), 2000, idate(2)
       WRITE (gol,'("WARNING - using dust emissions for 2000 ...")'); CALL goPr
       WRITE (gol,'("Filename for dust fluxes : ",a)') TRIM(dust_filename); CALL goPr

    CASE( "TEGEN-VIGNATI" ) 

       ! files 'dust200201.nc' etc.
       WRITE (dust_filename,'(a,"/dust",i4.4,i2.2,".nc")') TRIM(emis_input_dir_dust), idate(1), idate(2)

    CASE( "ONLINE" ) 

       ! stuff is done above (in *calc)
       status=0; RETURN

    CASE default 

       WRITE (gol,'("no valid dust emission method provided.")') ; CALL goErr
       TRACEBACK
       status=1; RETURN

    END SELECT

    IF (isRoot) THEN
       ! closing eventual previous dust file; ignore failure
       status = sfend( dust_FileID )

       dust_FileID = sfstart(TRIM(dust_fileName),DFACC_READ)
       IF ( dust_FileID == -1 ) THEN
          WRITE (gol,'("opening dust file : ",a)') TRIM(dust_fileName); CALL goErr
          WRITE (gol,'("in ",a)') rname; CALL goErr; status=1; RETURN
       ENDIF
    ENDIF

    
    ! ok
    status = 0

  END SUBROUTINE emission_dust_declare
  !EOC


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_DUST_DONE
  !
  ! !DESCRIPTION: Close open file(s). 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_DUST_DONE
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - v0
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    INTEGER    :: sfend

    DEALLOCATE(      uth )
    DEALLOCATE(    srel  )
    DEALLOCATE(    srelV )
    DEALLOCATE( su_srelV )

    ! gridded 1x1 fields from input file
    DEALLOCATE( soil_type  )
    DEALLOCATE( pot_source )
    DEALLOCATE( cult       )
    DEALLOCATE( area       )
    DEALLOCATE( z0         )
    DEALLOCATE( fpar       )
    DEALLOCATE( soilph     )

    ! additional 1x1 fields
    DEALLOCATE( snowcover )
    DEALLOCATE( desert    )
    DEALLOCATE( umin2     )
    DEALLOCATE( alpha     )
    DEALLOCATE( c_eff     )
    DEALLOCATE( lai_eff   )

    IF( isRoot .AND. ( TRIM( emis_input_dust ) /= 'ONLINE' )) THEN
       WRITE(gol,*) 'Closing dust emission file', isRoot; CALL goPr
       WRITE(gol,*) 'Closing dust emission file:', sfend(dust_FileID); CALL goPr
    END IF

  END SUBROUTINE EMISSION_DUST_DONE
  !EOC


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    READ_EMISSION_DUST
  !
  ! !DESCRIPTION: Opens, reads and evaluates input files (per month). 
  !               Provides emissions on 2d/3d-arrays which are then given 
  !               to emis_number and emis_mass, which are used in 
  !               sedimentation. (no *_apply!)
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE READ_EMISSION_DUST( status )
    !
    ! !USES:
    !
    USE dims,          ONLY : newday, nlon360, nlat180, idate, mlen, im, jm, lm, itau, okdebug, nregions
    USE dims,          ONLY : iglbsfc
    USE chem_param,    ONLY : sigma_lognormal, dust_density, nmodes, mode_aci, mode_coi
    USE chem_param,    ONLY : iaci_n, iduaci, icoi_n, iducoi
    USE toolbox,       ONLY : coarsen_emission
    USE Binas,         ONLY : pi
    USE datetime,      ONLY : tau2date
    USE emission_data, ONLY : emis_mass, emis_number, emis_temp
    USE emission_data, ONLY : vd_class_name_len, emission_vdist_by_sector
    !
    ! !OUTPUT PARAMETERS:
    !
    INTEGER, INTENT(out)    ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - added TEGEN-VIGNATI option
    !   26 Jun 2012 - P. Le Sager  - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:
    ! (1) now read on root, need to switch (FIXME) to MDF for parallel I/O
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    CHARACTER(len=*), PARAMETER :: rname = mname//'/read_emission_dust'

    REAL(kind=4), DIMENSION(:,:), ALLOCATABLE        :: mode_number, mode_radius, mode_mass
    INTEGER                                          :: daynumber, var_id, region, istat, sds_id
    INTEGER                                          :: sfn2index,sfselect,sfrdata,sfendacc
    REAL                                             :: ln2s
    CHARACTER(len=20)                                :: var_name
    CHARACTER(len=1 )                                :: modes
    INTEGER, PARAMETER                               :: add_field = 0
    !INTEGER, PARAMETER                               :: amonth=2
    INTEGER                                          :: imr, jmr, lmr, imode
    INTEGER, DIMENSION(6)                            :: idater

    INTEGER                              :: i, j, i1, i2, j1, j2, i01, i02, j01, j02
    TYPE(emis_data), DIMENSION(nregions) :: emis_glb
    TYPE(d3_data)                        :: emis3d
    CHARACTER(len=vd_class_name_len)     :: splittype


    ! --- begin -----------------------------------------

    ! vertical splitting is that class
    splittype = 'nearsurface'

    ! don't do anything for online method
    IF( TRIM(emis_input_dust) == 'ONLINE' ) RETURN

    IF(newday) THEN   

       !===================
       ! Work arrays
       !===================
!       CALL GET_DISTGRID( dgrid(iglbsfc), I_STRT=i01, I_STOP=i02, J_STRT=j01, J_STOP=j02 )
!
!       allocate(number(i01:i02,j01:j02))
!       allocate(mass  (i01:i02,j01:j02))

       IF (isRoot) THEN

          ALLOCATE(mode_number(nlon360,nlat180))
          ALLOCATE(mode_mass  (nlon360,nlat180))

          IF( TRIM(emis_input_dust) /= "TEGEN-VIGNATI" ) &
               ALLOCATE(mode_radius(nlon360,nlat180))

          DO region = 1, nregions
             ALLOCATE(emis_glb(region)%surf(im(region), jm(region)))
          END DO

       ELSE
             
          DO region = 1, nregions
             ALLOCATE(emis_glb(region)%surf(1,1))
          END DO
          
       END IF


       CALL tau2date(itau-3*3600, idater)
       idater(4) = 21

       daynumber = idate(3)

       !===================
       ! Read accumulation mode
       !===================
       imode = 2
       WRITE(modes,'(i1)') imode

       IF (isRoot) THEN

          IF( TRIM(emis_input_dust) == "TEGEN-VIGNATI" ) THEN 
             var_name = 'mode'//modes//'_mass'
             var_id   = sfn2index(dust_FileID,TRIM(var_name))
             sds_id   = sfselect (dust_FileID,var_id)
             istat    = sfrdata(sds_id,(/0,0,daynumber-1/), (/1,1,1/), (/nlon360,nlat180,1/),mode_mass)
          ELSE 
             var_name = 'mode'//modes//'_radius'
             var_id   = sfn2index(dust_FileID,TRIM(var_name))
             sds_id   = sfselect (dust_FileID,var_id)
             istat    = sfrdata(sds_id,(/0,0,daynumber-1/), (/1,1,1/), (/nlon360,nlat180,1/),mode_radius)
          END IF

          status=istat
          IF_NOTOK_RETURN(status=1)

          istat    = sfendacc(sds_id)
          var_name = 'mode'//modes//'_number'
          var_id   = sfn2index(dust_FileID,TRIM(var_name))
          sds_id   = sfselect (dust_FileID,var_id)
          istat    = sfrdata(sds_id,(/0,0,daynumber-1/), (/1,1,1/), (/nlon360,nlat180,1/),mode_number)

          status=istat
          IF_NOTOK_RETURN(status=1)

          istat    = sfendacc(sds_id)

          IF( TRIM(emis_input_dust) == "TEGEN-VIGNATI" ) THEN 
             mode_number = mode_number * mlen(idate(2))   ! #/gridbox/day ---> #/gridbox/month
             mode_mass   = mode_mass   * mlen(idate(2))   ! #/gridbox/day ---> #/gridbox/month
          ELSE 
             ln2s = (alog(sigma_lognormal(mode_aci)))**2
             ! transform to #/gridbox  kg/gridbox and shift...
             mode_number = CSHIFT(mode_number,nlon360/2,1)   ! shift emissions starting at 0 degrees to -180
             mode_radius = CSHIFT(mode_radius,nlon360/2,1) 
             mode_number = mode_number*mlen(idate(2))   ! #/gridbox/day ---> #/gridbox/month
             ! 1e18: um3 --> m3; exp(9/2...) is for Volume(r_geom^3) to Volume(r_mass^3)
             mode_mass=pi*4./3.*mode_radius**3.*mode_number*EXP(9./2.*ln2s) /1e18*dust_density   ! kg/gridbox/month
          END IF

       END IF
       
       ! Coarsen, scatter, vertical distribution
       ! ----------------------       
       IF (isRoot) THEN
          CALL coarsen_emission('dust_number ', nlon360, nlat180, REAL(mode_number) , emis_glb, add_field, status)
          IF_NOTOK_RETURN(status=1)
       END IF

       DO region = 1, nregions
          CALL scatter(dgrid(region), emis_temp(region)%surf, emis_glb(region)%surf, 0, status)
          IF_NOTOK_RETURN(status=1)
          CALL get_distgrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
          ALLOCATE( emis3d%d3(i1:i2,j1:j2,lm(region)) ) ; emis3d%d3 = 0.0
          CALL emission_vdist_by_sector( splittype, 'DUST', region, emis_temp(region), emis3d, status )
          emis_number(region,mode_aci)%d4(:,:,:,1)   = emis3d%d3
          DEALLOCATE(emis3d%d3)
       ENDDO

       ! Coarsen, scatter, vertical distribution
       ! ----------------------       
       IF (isRoot) THEN
          CALL coarsen_emission('dust_mass ', nlon360, nlat180, REAL(mode_mass) , emis_glb, add_field, status)
          IF_NOTOK_RETURN(status=1)
       END IF

       DO region = 1, nregions
          CALL scatter(dgrid(region), emis_temp(region)%surf, emis_glb(region)%surf, 0, status)
          IF_NOTOK_RETURN(status=1)

          CALL get_distgrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
          ALLOCATE( emis3d%d3(i1:i2,j1:j2,lm(region)) ) ; emis3d%d3 = 0.0
          CALL emission_vdist_by_sector( splittype, 'DUST', region, emis_temp(region), emis3d, status )
          emis_mass(region,mode_aci)%d4(:,:,:,1)   = emis3d%d3
          DEALLOCATE(emis3d%d3)
       ENDDO


       !===================
       ! Read coarse mode
       !===================
       imode = 3
       WRITE(modes,'(i1)') imode

       IF (isRoot) THEN 

          IF( TRIM(emis_input_dust) == "TEGEN-VIGNATI" ) THEN 
             var_name = 'mode'//modes//'_mass'
             var_id   = sfn2index(dust_FileID,TRIM(var_name))
             sds_id   = sfselect (dust_FileID,var_id)
             istat    = sfrdata(sds_id,(/0,0,daynumber-1/), (/1,1,1/), (/nlon360,nlat180,1/),mode_mass)
          ELSE 
             var_name = 'mode'//modes//'_radius'
             var_id   = sfn2index(dust_FileID,TRIM(var_name))
             sds_id   = sfselect (dust_FileID,var_id)
             istat    = sfrdata(sds_id,(/0,0,daynumber-1/), (/1,1,1/), (/nlon360,nlat180,1/),mode_radius)
          END IF

          status=istat
          IF_NOTOK_RETURN(status=1)

          istat    = sfendacc(sds_id)
          var_name = 'mode'//modes//'_number'
          var_id   = sfn2index(dust_FileID,TRIM(var_name))
          sds_id   = sfselect (dust_FileID,var_id)
          istat    = sfrdata(sds_id,(/0,0,daynumber-1/), (/1,1,1/), (/nlon360,nlat180,1/),mode_number)

          status=istat
          IF_NOTOK_RETURN(status=1)

          istat    = sfendacc(sds_id)

          IF( TRIM(emis_input_dust) == "TEGEN-VIGNATI" ) THEN 
             mode_number = mode_number*mlen(idate(2))   ! #/gridbox/day ---> #/gridbox/month
             mode_mass = mode_mass*mlen(idate(2))    ! #/gridbox/day ---> #/gridbox/month
          ELSE 
             ln2s = (alog(sigma_lognormal(mode_coi)))**2
             ! transform to #/gridbox  kg/gridbox and shift...
             mode_number = CSHIFT(mode_number,nlon360/2,1)   ! shift emissions starting at 0 degrees to -180
             mode_radius = CSHIFT(mode_radius,nlon360/2,1) 
             mode_number = mode_number*mlen(idate(2))   ! #/gridbox/day ---> #/gridbox/month
             ! 1e18: um3 --> m3; exp(9/2...) is for Volume(r_geom^3) to Volume(r_mass^3)
             mode_mass=pi*4./3.*mode_radius**3.*mode_number*EXP(9./2.*ln2s) /1e18*dust_density   ! kg/gridbox/month
          END IF

       END IF

       ! Coarsen, scatter, vertical distribution
       ! ----------------------       
       IF (isRoot) THEN
          CALL coarsen_emission('dust_number ', nlon360, nlat180, REAL(mode_number) , emis_glb, add_field, status)
          IF_NOTOK_RETURN(status=1)
       END IF

       DO region = 1, nregions
          CALL scatter(dgrid(region), emis_temp(region)%surf, emis_glb(region)%surf, 0, status)
          IF_NOTOK_RETURN(status=1)
          CALL get_distgrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
          ALLOCATE( emis3d%d3(i1:i2,j1:j2,lm(region)) ) ; emis3d%d3 = 0.0
          CALL emission_vdist_by_sector( splittype, 'DUST', region, emis_temp(region), emis3d, status )
          emis_number(region,mode_coi)%d4(:,:,:,1)   = emis3d%d3
          DEALLOCATE(emis3d%d3)
       ENDDO

       ! Coarsen, scatter, vertical distribution
       ! ----------------------       
       IF (isRoot) THEN
          CALL coarsen_emission('dust_mass ', nlon360, nlat180, REAL(mode_mass) , emis_glb, add_field, status)
          IF_NOTOK_RETURN(status=1)
       END IF

       DO region = 1, nregions
          CALL scatter(dgrid(region), emis_temp(region)%surf, emis_glb(region)%surf, 0, status)
          IF_NOTOK_RETURN(status=1)

          CALL get_distgrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
          ALLOCATE( emis3d%d3(i1:i2,j1:j2,lm(region)) ) ; emis3d%d3 = 0.0
          CALL emission_vdist_by_sector( splittype, 'DUST', region, emis_temp(region), emis3d, status )
          emis_mass(region,mode_coi)%d4(:,:,:,1)   = emis3d%d3
          DEALLOCATE(emis3d%d3)
       ENDDO

       
       !===================
       ! Done
       !===================

       DO region = 1, nregions
          DEALLOCATE(emis_glb(region)%surf)
       END DO

       IF (isRoot) THEN
          DEALLOCATE(mode_number, mode_mass)

          IF( TRIM(emis_input_dust) /= "TEGEN-VIGNATI" ) &
               DEALLOCATE(mode_radius)
       END IF
       

    ENDIF ! newday  

  END SUBROUTINE READ_EMISSION_DUST
  !EOC
  
END MODULE EMISSION_DUST