! #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