!-------------------------------------------------------------------------------------- ! Module for mineral dust fluxes calculations with various options ! The saltation/sandblasting scheme: ! ifluxv=0: [Marticorena and Bergametti, 1995] ! ifluxv=1: [Alfaro and Gomes, 2001] ! ifluxv=2: [Alfaro and Gomes, 2001] + optimization of [Menut et al., 2005] ! ifluxv=3: [Kok, 2014] ! ! uth: the smooth ustar threshold [Shao and Lu, 2000] or [Iversen and White 1982] ! uth2: the smooth ustar threshold, with fine resolution ! emufac: factor to convert g/cm2/s into molec/cm2/s ! emiprof: normalized function to convert modes to bins ! ! Dimensions for the wind: all parameters have to be in [cm] ! w10m(nz,nm,ks) : in [cm/s], already converted in met2chim ! wsta(nz,nm,ks) : in [cm/s], already converted in met2chim ! utmin: fit of (Gillette, 1998, JGR) ! ! locU_ms, locW_ms: are loc U and W in m/s ! ! Weibull distribution: ! lambda=locU_ms/gamma(1+1/kref) ! Gamma function estimated with the Stirling formula ! The distribution may have a variable shape kref, depending on the sub-grid scale topography ! ikref=1: Constant value k=4 ! =2: Grini et al 2004 k=0.94 sqrt(U) ! =3: sub-grid scale orography test ! =4: max(4.,U**0.8) ! ! coheng: Binding energies of the 3 aerosol populations (g.cm2/s2) ! Initial values proposed by [Alfaro et al. 1998] are e_ = (/ 3.61, 3.52, 3.46 /) ! They are often corrected using a correction factor 'div' as e=e_/div ! For div increase, cohesion energies decrease, then emissions increase. !-------------------------------------------------------------------------------------- module worker_dust use chimere_common use gridvars, only : xlong,xlati, nzonal, nmerid use interp_tools implicit none ! Width of the Weibull distribution integer, parameter :: ikref=2 ! [Alfaro and Gomes, 2001] parameters real(kind=iprec), parameter :: div = 3d0 real(kind=iprec), dimension(nmode), parameter :: coheng = (/3.61d0/div, 3.52d0/div, 3.46d0/div/) ! Parameters for the soil size distribution integer, parameter :: nsoildmax=20000 ! number of bins for the soil size distribution real(kind=iprec), parameter :: soil_dmin = 1.d-4 ! soil size distribution min (cm) real(kind=iprec), parameter :: soil_dmax = 2.d-1 ! soil size distribution max (cm) real(kind=iprec), parameter :: z1000cm = 1.d+3 ! 10m height in [cm] real(kind=iprec), parameter :: rop = 2.65d0 ! particle density (g cm-3) i.e quartz real(kind=iprec), parameter :: snow_th=1.d-3 ! Snow cover if more that [snow_th] m real(kind=iprec), parameter :: minsurf_th=0.d0 ! No flux if more than (minsurf_th x 100) % of water real(kind=iprec), parameter :: maxemisdust=5.d-8 ! max dust flux in g/cm2/s (unit used in scheme) real(kind=iprec), dimension(nmineral), parameter :: mineral_dp = (/0.0690, 0.021, 0.0125, 0.0002, 0.052/) real(kind=iprec), dimension(nmineral), parameter :: mineral_sig = (/1.6, 1.8, 1.6, 2., 1.5/) real(kind=iprec), allocatable, dimension(:) :: size1, size2, uth, uth2 !------------ contains !------------------------------------------------------------------------------------ subroutine inidust_worker implicit none ! print settings if (rank.eq.0.and.chimverb.ge.2) then if(ifluxv.eq.0) print*,'- Init dust emissions [Marticorena Bergametti 1995]' if(ifluxv.eq.1) print*,'- Init dust emissions [Alfaro Gomes 2001]' if(ifluxv.eq.2) print*,'- Init dust emissions [Alfaro Gomes 2001] optimized' if(ifluxv.eq.3) print*,'- Init dust emissions [Kok et al 2014]' endif ! u*t and relative soil surface distributions call uthsoil ! Pre-calc constants before adjustement for the modified [Alfaro and Gomes, 2001] parameterization if(ifluxv.eq.2) call prepag01modif end subroutine inidust_worker !====================================================================================== ! Calculation of dust fluxes for one hour noted 'ks' ! Result is in emisd(ks) in molecules/cm2/s ! maxvoro: maximum of topography variance on a whole domain. Here fixed. ! aodthreshold: the MMD size (in um) for the calculation of the fine/coarse AOD !------------------------------------------------------------- subroutine diagdust_hour(ks) implicit none integer,intent(in) :: ks integer :: sp, sm, si, no, isoil, nz, nm, kfin, ilupc, kwb, k, nb, nd integer :: im1, im2, icurdate, itestwater, iaodt ! Parameters for the u* threshold real(kind=iprec), parameter :: umin = 21.d0 ! Min wind speed for the 'Owen' effect (cm.s-1) ! Weibull distribution integer, parameter :: nwb=15 ! Number of values for the distribution real(kind=iprec), parameter :: dwindmax=15d0 ! Maximum of wind speed (m/s) for use of Weibull ! Parameters for the partition function between clay and silt real(kind=iprec), parameter :: ca=20.d0 real(kind=iprec), parameter :: cb=2.d0 real(kind=iprec), parameter :: cc=0.8d0 ! Parameter to refine the AOD calculation real(kind=iprec), parameter :: aodthreshold=20. ! constants for the sub-grid variance of orography (m) real(kind=iprec), parameter :: kmin=0.8 real(kind=iprec), parameter :: kmax=1.2 real(kind=iprec), parameter :: ka=20.d0 real(kind=iprec), parameter :: kb=10.d0 real(kind=iprec), parameter :: kc=1.d0 real(kind=iprec) :: oroslope, zratio, w1, w2, kreft, uth_eff, factflux real(kind=iprec) :: locemidhr(nbinshr) real(kind=iprec), dimension(ms) :: locemisd real(kind=iprec) :: emicumul, emisolub real(kind=iprec) :: aodlr(nzonal,nmerid), aodhr(nzonal,nmerid) real(kind=iprec) :: conctau(ntau), ratioaod real(kind=iprec) :: locU_ms, locW_ms, modu10cm, modu10m, ft1, ft2, ft3, locsoim, locrelh real(kind=iprec) :: cdnms, soilpcent, cerod, merod, lupcent real(kind=iprec) :: z0dyn, pdfu, auxreal, weilambda real(kind=iprec) :: probu(nwb), wind10ms(nwb), wind10cms(nwb) real(kind=iprec) :: factrs, factr, facts, r_ratio, r_alpha, r_tau real(kind=iprec) :: num, den, dvd_dlnD, clay real(kind=iprec) :: ustarns, utmin, kref real(kind=iprec) :: flux1,flux2,flux3 real(kind=iprec), allocatable,dimension(:,:,:) :: fluxegrid, emisdtot real(kind=iprec), allocatable,dimension(:) :: fsilt, fclay, cumul real(kind=iprec), allocatable,dimension(:,:,:) :: emisd_clay, emisd_silt !--------------------------------------------------------------- if (.not. allocated(fluxegrid))allocate(fluxegrid(nmode,nzonal,nmerid)) if (.not. allocated(emisdtot))allocate(emisdtot(ms,nzonal,nmerid)) if(iusedust.ge.2)then if (.not. allocated(fsilt))allocate(fsilt(ms)) if (.not. allocated(fclay))allocate(fclay(ms)) if (.not. allocated(cumul))allocate(cumul(ms)) if (.not. allocated(emisd_clay))allocate(emisd_clay(ms,nzonal,nmerid)) if (.not. allocated(emisd_silt))allocate(emisd_silt(ms,nzonal,nmerid)) endif ! Searching for the time in the month for the MODIS erodibility call reldat(idatestart,ihourrun,icurdate) call hourinmonth(icurdate,w1,w2,im1,im2) ! Saltation z0 and soil efficiency call calc_z0s_effsoil emisd(:,:,:,:,ks)=dzero fluxegrid = dzero emisdtot= dzero do nm=1,nmerid do nz=1,nzonal call check_water(nz,nm,itestwater) if(itestwater.eq.1)cycle locU_ms=abs(w10m(nz,nm,ks)/100.) locW_ms=abs(wsta(nz,nm,ks)/100.) modu10m=sqrt(locU_ms**2d0+(1.2d0*locW_ms)**2d0) locsoim=soim(nz,nm,ks) locrelh=sreh(nz,nm,ks) if(topc(nz,nm,ks).gt.dzero.or.locrelh.gt.0.99999d+00)locrelh=0.99999d+0 merod=w1*dataerod(nz,nm,im1)+w2*dataerod(nz,nm,im2) ! Constant value if(ikref.eq.1)kref=4. ! Grini et al 2004 (relation to use with U10m in m/s) if(ikref.eq.2)kref=0.94 * sqrt(modu10m) ! sub-grid scale orography test if(ikref.eq.3)then zratio=varsso(nz,nm)/1000. oroslope=1.-(kc)/(1.+ka*exp(-kb*zratio)) kreft=kmin+(kmax-kmin)*oroslope kref=kreft * (0.94 * sqrt(modu10m)) endif ! Weibull wind speed distribution auxreal=dun+dun/kref weilambda = modu10m/exp(auxreal*log(auxreal)-auxreal & & - 0.5d0*log(auxreal/(2.d0*pi))+dun/(12.d0*auxreal) & & - dun/(360.d0*(auxreal**3.d0))+dun/(1260.d0*(auxreal**5.e0))) do kwb=1,nwb if (nwb.gt.1) then wind10ms(kwb) = kwb*2d0*modu10m/nwb pdfu = (kref/weilambda)*(wind10ms(kwb)/weilambda)**(kref-dun) & * exp(-(wind10ms(kwb)/weilambda)**kref) probu(kwb) = pdfu*2.d0*modu10m/nwb if(probu(kwb).le.1.d-2) cycle else wind10ms(kwb) = modu10m probu(kwb) = dun endif ! Limit to avoid unrealistic high values of the wind speed if(wind10ms(kwb).ge.dwindmax)wind10ms(kwb)=dwindmax enddo do kwb=1,nwb wind10cms(kwb)=wind10ms(kwb)*100d0 ! Loop over the soil and surface characteristics ! The index 4 of table 'mineral_frac' represents the clay's % of the soil in the current cell flux1=dzero flux2=dzero flux3=dzero do isoil=1,nsoils do ilupc=1,nlupcell ceff= feff(nz,nm,ilupc,isoil) z0dyn= surf_z01(nz,nm,ilupc,isoil) call calc_dust_erod(nz,nm,ilupc,isoil,merod,cerod) soilpcent= SOILrelsurf(nz,nm,isoil) lupcent= ratiol(nz,nm,ilupc,isoil) clay=mineral_frac(isoil,4) ustarsalt=dzero if(ceff.le.0.) cycle if(z0dyn.eq.0) cycle if(cerod.eq.0.) cycle if(soilpcent.eq.0) cycle if(lupcent.eq.0) cycle ! in dust schemes, U10m is now in cm/s cdnms= vkarm/(log(z1000cm/z0dyn)) ustarns=cdnms * wind10cms(kwb) utmin=umin/(cdnms*ceff) ustarsalt=ustarns if(wind10cms(kwb).ge.utmin)ustarsalt=ustarns+(0.3*(wind10cms(kwb)/m2cm-utmin/m2cm)**2.) if(ustarsalt.lt.umin/ceff)cycle call copyncl(kfin,isoil) if(ifecan.eq.1)call fecan1999(kfin,locsoim,clay) ft1=dzero ft2=dzero ft3=dzero if(ifluxv == 0) call mb1995(kfin,ft1,ft2,ft3,uth_eff,isoil) if(ifluxv == 1) call ag2001(kfin,ft1,ft2,ft3,uth_eff,isoil) if(ifluxv == 2) call ag2001(kfin,ft1,ft2,ft3,uth_eff,isoil) if(ifluxv == 3) call kok2014(nz,nm,kfin,ks,clay,ft1,isoil) factflux=soilpcent*cerod*lupcent*(1.-locrelh) flux1=flux1+ft1*factflux flux2=flux2+ft2*factflux flux3=flux3+ft3*factflux enddo ! end of loop over nlupcell enddo ! end of loop over nsoils fluxegrid(1,nz,nm) = fluxegrid(1,nz,nm)+flux1*probu(kwb) fluxegrid(2,nz,nm) = fluxegrid(2,nz,nm)+flux2*probu(kwb) fluxegrid(3,nz,nm) = fluxegrid(3,nz,nm)+flux3*probu(kwb) enddo ! of loop over wind speed loop enddo ! of nz enddo ! of nm 110 format(10(1x,f12.6)) !---------------------------------------------------------- ! Attenuation of mineral dust fluxes due to rain and/or snow ! In case of snow at the surface, no dust flux ! In case of rain during the last hours, attenuated dust flux estimation do nm=1,nmerid do nz=1,nzonal factr=dun facts=dun factrs=dun if(snowh_ok.eq.1.and.snowh(nz,nm,ks).gt.snow_th)facts=dzero r_tau=12.d0 r_ratio=-2.d0*pi*lastrain(nz,nm)/r_tau factr=dun-exp(r_ratio) if(facts.lt.dun.or.factr.lt.dun)factrs=min(facts,factr) fluxegrid(:,nz,nm)=fluxegrid(:,nz,nm)*factrs enddo enddo !---------------------------------------------------------- ! Limitation of the mineral dust flux to avoid unrealistic high values ! The upper limit is defined in chimere_params as: maxemisdust (g/m2/s) ! Default value: 5.10^-4 g/m2/s (i.e 5.10^-8 g/cm2/s) ! do nm=1,nmerid ! do nz=1,nzonal ! do k=1,nmode ! fluxegrid(k,nz,nm)=min(fluxegrid(k,nz,nm),maxemisdust) ! enddo ! enddo ! enddo !---------------------------------------------------------- ! Redistribute flux in three (nmode) modes into (ms) model bins. ! For MB1995 and AG2001, fluxes are emitted in three modes ! For Kok2014, fluxes are in one single mode (put in the first one of fluxegrid) ! Check AOD to unbias the 'bins' effect ! itvhr is a diameter in meter, need to be converted as radius in um. ! mmd_m: the mean mass median diameter (m) aodhr=0. aodlr=0. ! looking for aodthreshold in the bintau distribution ! aodthreshold diameter in um ! bintau: radius in um do no=1,ntau-1 if(bintau(no).le.aodthreshold/2..and.bintau(no+1).ge.aodthreshold/2.)iaodt=no enddo do nz=1,nzonal do nm=1,nmerid locemidhr=0. if(ifluxv.eq.0.or.ifluxv.eq.1.or.ifluxv.eq.2)then do k=1,nmode do nb=1,nbinshr locemidhr(nb)=locemidhr(nb)+fluxegrid(k,nz,nm)*emufac*emiprofhr(idmineral(1,k),nb) enddo enddo endif if(ifluxv.eq.3)then do nb=1,nbinshr locemidhr(nb)=fluxegrid(1,nz,nm)*emufac*emikokprof(nb) enddo endif call itpconcAPP(locemidhr,nbinshr,conctau,ntau,ntau,itvlog_hr,ilut_hrntau) do no=1,iaodt aodhr(nz,nm)=aodhr(nz,nm)+conctau(no)*lutau(no) enddo call cuminbins(mmdhr,locemidhr,nbinshr,ds,locemisd,ms) call itpconcAPP(locemisd,ms,conctau,ntau,ntau,itvlog_ms,ilut_msntau) do no=1,iaodt aodlr(nz,nm)=aodlr(nz,nm)+conctau(no)*lutau(no) enddo ratioaod=1. ! if(aodlr(nz,nm).ne.0.)ratioaod=aodhr(nz,nm)/aodlr(nz,nm) do nb=1,ms if(mmd_m(nb)*1.d6.le.aodthreshold)then emisdtot(nb,nz,nm)=locemisd(nb)*ratioaod else emisdtot(nb,nz,nm)=locemisd(nb) endif enddo enddo enddo !----------------------------------------------------------------------------------- ! If only one dust species is calculated if(iusedust.eq.1)then do nd=1,nemisd do nz=1,nzonal do nm=1,nmerid do nb=1,ms emisd(nd,nb,nz,nm,ks)=emisdtot(nb,nz,nm) enddo enddo enddo enddo endif !----------------------------------------------------------------------------------- ! If calculation of the mineral dust mineralogy. ! ! The vertical flux is modified to account for the clay and silt fraction. ! (hypothesis: the total mass of Fv is only composed of clay and silt) if(iusedust.eq.2.or.iusedust.eq.3)then emisd_clay= dzero emisd_silt= dzero do nb=1,ms fclay(nb)=1.d0-cc/(1.d0+ca*exp(-cb*mmd_m(nb)*1.d+6)) fsilt(nb)=1.d0-fclay(nb) enddo do nz=1,nzonal do nm=1,nmerid ! Fraction of silt and clay, assuming that the whole vertical flux is composed only of these two do nb=1,ms emisd_clay(nb,nz,nm)=emisdtot(nb,nz,nm)*fclay(nb) emisd_silt(nb,nz,nm)=emisdtot(nb,nz,nm)*fsilt(nb) enddo ! Estimation of the chemical elements fraction and calculation of soluble and insoluble parts ! clay_ratio: [0:1] ! silt_ratio: [0:1] ! mineral_ratio: Ratio of chemical element in mineral ! iusedust=2: Case of emission of chemical elements if(iusedust.eq.2)then do nd=1,nemisd if(idustother(nd).eq.0)then sm=idustchem(nd) si=idustsolub(nd) do sp=1,nmineralspec do nb=1,ms emicumul=emisd_clay(nb,nz,nm)*clay_ratio(nz,nm,sp)*mineral_ratio(sp,sm)/100.d0 & & +emisd_silt(nb,nz,nm)*silt_ratio(nz,nm,sp)*mineral_ratio(sp,sm)/100.d0 if(si.eq.1)then emisolub=mineral_solubility(sp,sm)/100.d0 else emisolub=(1d0-mineral_solubility(sp,sm)/100.d0) endif emisd(nd,nb,nz,nm,ks)=emisd(nd,nb,nz,nm,ks)+emicumul*emisolub enddo enddo endif enddo endif ! iusedust=3: Case of emission of mineral species if(iusedust.eq.3)then do nd=1,nemisd if(idustother(nd).eq.0)then sp=idustspec(nd) do nb=1,ms emisd(nd,nb,nz,nm,ks)=emisd_clay(nb,nz,nm)*clay_ratio(nz,nm,sp) & & +emisd_silt(nb,nz,nm)*silt_ratio(nz,nm,sp) enddo endif enddo endif ! The OTher species DuOT is always created (if iusedust>2) but is not in the mineral tables cumul=0d0 do nd=1,nemisd if(idustother(nd).eq.0)then do nb=1,ms cumul(nb)=cumul(nb)+emisd(nd,nb,nz,nm,ks) enddo endif enddo do nd=1,nemisd if(idustother(nd).eq.1)then do nb=1,ms emisd(nd,nb,nz,nm,ks)=emisdtot(nb,nz,nm)-cumul(nb) enddo endif enddo enddo enddo endif !----------------------------------------------- if (allocated(fluxegrid)) deallocate(fluxegrid) if (allocated(emisdtot)) deallocate(emisdtot) if(iusedust.ge.2)then if (allocated(fsilt)) deallocate(fsilt) if (allocated(fclay)) deallocate(fclay) if (allocated(cumul)) deallocate(cumul) if (allocated(emisd_clay)) deallocate(emisd_clay) if (allocated(emisd_silt)) deallocate(emisd_silt) endif end subroutine diagdust_hour !///////////////////////////////////////////////////////////////////////////////////////////// ! Initialization routines !================================================================================ !subroutine read_dust_inputdata !end subroutine read_dust_inputdata !================================================================================ ! uthsoil: ! Calculations of the smooth threshold friction velocity ! This routine should be surface air density dependent. ! Since the coefficients were fitted with a constant ROa value, ! this model version used fixed roa. ! icuth=0 Iversen and White 1982 ! icuth=1 Shao and Lu 2000 ! ! All units must be [g] and [cm] ! ! Relative soil surface distributions ! 3 modes ('nmodes') described by three log-normal distributions ! for each type of soil ('nsoils') ! for each distrib: ! - mineral_dp: diameter (cm) ! - mineral_sig: sigma of distrib ! - mineral_frac: relative fraction of this mode [0:1] ! ! After [Shao and Lu, 200, JGR]: ! sig_shao=rop/roa: the particle to air density ratio ! with rop = 2.65d0 (g cm-3) ! roa = 0.001227d0 (g cm-3) the air density (may be replaced by airm later!) ! gam_shao: 3.e-4 kg.s-2 then 3.e-1 g.s-2 ! an_shao=0.0123 [ad.] !-------------------------------------------------------------------------------- subroutine uthsoil implicit none integer :: ns, i, j real(kind=iprec), parameter :: adust = 1331.647d0 real(kind=iprec), parameter :: bdust = 0.38194d0 real(kind=iprec), parameter :: xdust = 1.561228d0 real(kind=iprec), parameter :: ddust=6d-3/(rop*gravitcm) real(kind=iprec), parameter :: soil_eps=1.d-5 real(kind=iprec), parameter :: dstep=2.d-3 ! test au lieu de 5d-3 real(kind=iprec) :: an_shao, gam_shao, sig_shao, x1, x2 real(kind=iprec) :: bb, cc, xk, dd, ee, ff, xl, xm, xn real(kind=iprec) :: stotale, su, su_loc, dp2 real(kind=iprec), allocatable, dimension(:) :: subsoildist real(kind=iprec), allocatable, dimension(:) :: newstep !-------------------------------------------------------- ! Looking for the number of soil classes dp2=soil_dmin do i=1,nsoildmax dp2=dp2*exp(dstep) if (dp2.ge.(soil_dmax+soil_eps).or.ncl.eq.nsoildmax) then ncl=i exit endif enddo ! Allocate table knowing the number of soil classes if(.not. allocated(srel)) allocate(srel(nsoils,ncl)) if(.not. allocated(srel2)) allocate(srel2(nsoils,ncl)) if(.not. allocated(size1)) allocate(size1(ncl)) if(.not. allocated(size2)) allocate(size2(ncl)) if(.not. allocated(subsoildist)) allocate(subsoildist(ncl)) if(.not. allocated(newstep)) allocate(newstep(ncl)) if(.not. allocated(uth)) allocate(uth(ncl)) if(.not. allocated(uth2)) allocate(uth2(ncl)) ! Calculate the soil size distribution dp2=soil_dmin size1(1) = soil_dmin newstep(1) = dstep do i=2,ncl dp2=dp2*exp(dstep) size1(i) = dp2 newstep(i) = size1(i)-size1(i-1) enddo !-------------------------------------------------------- ! icuth=0: Iversen and White 1982 if (icuth.eq.0) then do i=1,ncl bb=adust*(size1(i)**xdust)+bdust cc=sqrt(dun+ddust*(size1(i)**(-2.5d0))) xk=sqrt(abs(rop*gravitcm*size1(i)/roa)) if (bb < 10d0) then dd=sqrt(1.928d0*(bb**0.092d0)-dun) uth(i) = 0.129d0*xk*cc/dd else ee=-0.0617d0*(bb-10d0) ff=dun-0.0858d0*exp(ee) uth(i) = 0.12d0*xk*cc*ff endif enddo endif !-------------------------------------------------------- ! icuth=1: Shao and Lu 2000 if(icuth.eq.1) then an_shao=0.0123d0 gam_shao=0.3d0 do i=1,ncl sig_shao=rop/roa x1=sig_shao*gravitcm*size1(i) x2=gam_shao/(roa*size1(i)) uth(i) = sqrt(an_shao*(x1+x2)) enddo endif !-------------------------------------------------------- ! Relative soil surface distributions ! Calculation could be on 'nmineral' but is limited to the first three populations: ! coarse sand, fine medium sand, silt do ns=1,nsoils stotale=dzero subsoildist=dzero do i=1,ncl su=dzero do j=1,3 if (mineral_dp(j).eq.0.) then su_loc=0. else xk=mineral_frac(ns,j)/(sqrt(2.d0*pi)*log(mineral_sig(j))) xl = ((log(size1(i))-log(mineral_dp(j)))**2d0)/(2.d0*(log(mineral_sig(j)))**2d0) xm=xk*exp(-xl) xn=rop*(2.d0/3.d0)*(size1(i)/2.d0) su_loc=xm*newstep(i)/xn endif su=su+su_loc enddo subsoildist(i) = su stotale=stotale+su enddo do i=1,ncl if (subsoildist(i) > 0. .and. stotale > 0.) then srel(ns,i) = subsoildist(i)/stotale else srel(ns,i) = dzero endif enddo enddo end subroutine uthsoil !========================================================================= subroutine prepag01modif ! pre-calculation of constants before adjustement ! only needed for the modified [Alfaro and Gomes, 2001] parameterization ! Original values: ! e1=0.376 kg.m-2.s-2 ! e2=0.366 kg.m-2.s-2 ! e3=0.346 kg.m-2.s-2 ! Kinetic energy values are defined in chimere_params and expressed in g.cm-2.s-2 ! Optimization following [Menut et al., JAOTech, 2005] implicit none real(kind=iprec) :: et real(kind=iprec), parameter :: p3t=1.d-3 !---------------------- ! check binding energy values if (coheng(3) >= coheng(2) .or. coheng(2) >= coheng(1) .and. rank == 0) then print*,'Incorrect energy values : ',coheng(3),'<',coheng(2),'<',coheng(1) call exit1('') endif var3a=((3*coheng(3))/(100d0*rop*pi))**(1.d0/3.d0) et=coheng(3)+sqrt(((coheng(2)-coheng(3))*(coheng(1)-coheng(3)))/p3t) var3b=((3*et)/(100d0*rop*pi))**(1.d0/3.d0) end subroutine prepag01modif !========================================================================= subroutine copyncl(kfin,isoil) ! Estimation of the reduced soil size distribution ! Calculation depends on the u*th value. u*th is different for each time-step, each cell and each ! soil type. ! if the scheme is NOT AG2001 optmized, we just copy the soil size distribution ! if the scheme IS AG2001, we reduce the soil size distribution to speed up calculation implicit none integer :: i, ideb, ifin, k, kfin, isoil real(kind=iprec) :: dsmin, dsmax, srelmean integer, parameter :: istep=5 srel2(isoil,:)=0. if(ifluxv.ne.2)then do i=1,ncl size2(i) = size1(i) uth2(i) = uth(i) srel2(isoil,i) = srel(isoil,i) enddo kfin=ncl endif if(ifluxv.eq.2) then dsmin=var3a*(ustarsalt**(-2.d0/3.d0)) dsmax=var3b*(ustarsalt**(-2.d0/3.d0)) ! check consistency of values if(dsmin.lt.size1(1))dsmin=size1(1) if(dsmax.gt.size1(ncl))dsmax=size1(ncl) ! Looking for the position of dsmin and dsmax in the [1:ncl] soil distribution ideb=0 ifin=0 call dichotomy(size1,ncl,1,ncl,dsmin,ideb) call dichotomy(size1,ncl,1,ncl,dsmax,ifin) kfin=0 do i=ideb,ifin kfin=kfin+1 size2(kfin) = size1(i) uth2(kfin) = uth(i) srel2(isoil,kfin) = srel(isoil,i) enddo do i=ifin,ncl-istep,istep kfin=kfin+1 size2(kfin) = size1(i) uth2(kfin) = uth(i) srelmean=0d0 do k=1,istep srelmean=srelmean+srel(isoil,i+k-1) enddo if(kfin.gt.ncl) exit srel2(isoil,kfin) = srelmean enddo endif end subroutine copyncl !========================================================================= subroutine dichotomy(siz,ncl,i1,i2,size2find,iout) ! siz: array to scan the dimension 'nsoildmax' and reduced size of interest [i1:i2] ! size2find: value to find in 'siz' ! iout: index in the table 'siz' corresponding to location of 'size2find' implicit none integer, intent(in) :: ncl real(kind=iprec), intent(in) :: siz(ncl),size2find integer, intent(out) :: iout integer :: i1,i2,ismin,ismax,k2,ihalf,idiff iout=0 ismin=i1 ismax=i2 ihalf=int((ismax+ismin)/2.) do k2=1,1000000 if(size2find.ge.siz(ihalf))then ismin=ihalf else ismax=ihalf endif ihalf=int((ismax+ismin)/2.) idiff=ismax-ismin if(idiff.le.1)then iout=ismin exit endif enddo end subroutine dichotomy !========================================================================== subroutine fecan1999(kfin,locsoim,clay) ! Change the ustar threshold value accounting for soil moisture. ! w (noted gsoim) and w' (noted soimt) have to be in % [0:100] ! soim in [m3/m3], but gsoim must be in [%] ! clay_pc: The clay content in % (already the case with clay) ! clay_th: a threshold: Attenuation is not applied for cells with less than 10% of clay implicit none integer :: i,kfin real(kind=iprec) :: gsoim, soimt, fact_fecan, locsoim, clay, clay_pc real(kind=iprec), parameter :: clay_th=10d0 gsoim=100.d0*locsoim soimt=dzero fact_fecan=dzero clay_pc=clay if (clay_pc.gt.clay_th)then soimt = 0.0014d0 * (clay_pc)**2d0 + 0.17d0 * (clay_pc) if(gsoim.gt.soimt)then fact_fecan=sqrt(dun+1.21d0*((gsoim-soimt)**0.68d0)) endif endif if(fact_fecan.gt.dun)then do i=1,kfin uth2(i) = uth2(i)*fact_fecan enddo endif end subroutine fecan1999 !========================================================================== subroutine kok2014(nz,nm,kfin,ks,clay,ft1,isoil) ! [Kok et al., 2014] mineral dust emission scheme ! Values and constants in [kg], [m] and [s] ! ustarsalt must be converted here in [m/s] ! srel is the relative fraction of soil [0:1] ! At the end, the flux must be converted from [kg/m2] to [g/cm2] ! airmdust: conversion of airm to have kg/m3 ! ! Compared to [Kok 2014 ACP] ! fbare: fraction of the surface that consists of bare soil ! cerod is used here (the erodibility fraction) ! fclay: fraction of clay in the soil. We use here clay(%)/100 ! ! Original parameters: ! kok_ustarST0=0.16d0 ! m/s ! kok_Ce=2d0 ! [ad] ! kok_Cd0=4.4d-5 ! [ad] ! kok_alpha=2.7d0 ! [ad] ! airmdust0=1.225d0 ! kg/m3 implicit none integer :: i integer,intent(in) :: ks, nz, nm, kfin, isoil real(kind=iprec) :: airmdust, kok_Cd, ustarsalt_ms, ft1 real(kind=iprec) :: kokflux, fact1, fact2, fact3, ustarT_ms, clay, ustarST_ms real(kind=iprec), parameter :: kok_ustarST0=0.16d0 ! m/s real(kind=iprec), parameter :: kok_Ce=1.5d0 ! [ad] real(kind=iprec), parameter :: kok_Cd0=1.0d-4 ! [ad] real(kind=iprec), parameter :: kok_alpha=2.7d0 ! [ad] original 2.7 real(kind=iprec), parameter :: airmdust0=1.225d0 ! kg/m3 !------------------------------------------------------------------- ft1=dzero airmdust=airm(nz,nm,1,ks)/avogadro*29d0*1.d3 ustarsalt_ms=ustarsalt/m2cm do i=1,kfin kokflux=dzero ustarT_ms=uth2(i)/ceff/m2cm ustarST_ms=ustarT_ms*sqrt(airmdust/airmdust0) if(srel2(isoil,i).gt.0..and.ustarsalt_ms.gt.ustarT_ms)then kok_Cd=kok_Cd0*exp(-kok_Ce*(ustarST_ms-kok_ustarST0)/kok_ustarST0) fact1=airmdust*(ustarsalt_ms**2d0-ustarT_ms**2d0)/ustarST_ms fact2=ustarsalt_ms/ustarST_ms fact3=kok_alpha*(ustarST_ms-kok_ustarST0)/kok_ustarST0 kokflux=kok_Cd * srel2(isoil,i) * fact1*(fact2**fact3) * clay/100. ft1=ft1+kokflux endif enddo ! from [kg/m2] to [g/cm2] -> *1e3 / 1e4 => /10. ft1=ft1/10d0 end subroutine kok2014 !========================================================================== subroutine mb1995(kfin,ft1,ft2,ft3,uth_eff,isoil) ! [Marticorena and Bergametti, 1995] parameters ! ! cad: percentage of flux in each of the three modes ! uth and uth2: the smooth threshold friction velocity ! ustar: friction velocity under neutral conditions ! ceff: drag efficiency implicit none integer :: i integer,intent(in) :: kfin, isoil real(kind=iprec), dimension(nmode), parameter :: cad = (/0.2, 0.6, 0.2/) real(kind=iprec) :: mbflux, fdp1, fdp2, uth_eff, ft1, ft2, ft3 do i=1,kfin uth_eff=uth2(i)/ceff if(srel2(isoil,i).gt.0..and.ustarsalt.gt.uth_eff)then fdp1=dun-(uth_eff/ustarsalt) fdp2=(dun+(uth_eff/ustarsalt))**2d0 mbflux=cd*fdp1*fdp2*(ustarsalt**3d0)*alpha_dust(isoil)*srel2(isoil,i) ft1=ft1+cad(1)*mbflux ft2=ft2+cad(2)*mbflux ft3=ft3+cad(3)*mbflux endif enddo end subroutine mb1995 !============================================================================== subroutine ag2001(kfin,ft1,ft2,ft3,uth_eff,isoil) ! ! [Alfaro and Gomes 2001] parameters ! ! diam_aer must be in meters (and is converted in cm in the equation) ! coheng: cohesion energy between sand grains ! beta_dust: tunable parameter found in [AG01] as beta=16300 cm/s2 implicit none integer :: i real(kind=iprec), parameter :: beta_dust=16300.d0 integer,intent(in) :: kfin, isoil real(kind=iprec) :: t1, t2, t3, p1, p2, p3, dfec1, dfec2, dfec3, uth_eff real(kind=iprec) :: hflux, dec, fdp1, fdp2, ad1, ad2, ad3, ec, ft1, ft2, ft3 do i=1,kfin uth_eff=uth2(i)/ceff ad1=dzero ad2=dzero ad3=dzero if(srel2(isoil,i).gt.0..and.ustarsalt.gt.uth_eff)then fdp1=dun-(uth_eff/ustarsalt) fdp2=(dun+(uth_eff/ustarsalt))**2d0 hflux=cd*fdp1*fdp2*(ustarsalt**3.d0)*srel2(isoil,i) dec=hflux*beta_dust ec = (pi/3.d0)*1d+2*rop*(size2(i)**3.d0) *(ustarsalt**2.d0) dfec1 = (ec-coheng(1)) dfec2 = (ec-coheng(2)) dfec3 = (ec-coheng(3)) t1=dzero t2=dzero t3=dzero if(ec >= coheng(1))t1=1d0 if(ec >= coheng(2))t2=1d0 if(ec >= coheng(3))t3=1d0 if(dfec3.ne.dzero)then p1=t1*dfec1/dfec3 p2=t2*(dun-p1)*dfec2/dfec3 else p1=dzero p2=dzero endif p3=t3*(dun-p1-p2) ad1 = (pi/6.d0)*dec*rop*p1*((diam_aer(idmineral(1,1))*1d+2)**3.d0)/coheng(1) ad2 = (pi/6.d0)*dec*rop*p2*((diam_aer(idmineral(1,2))*1d+2)**3.d0)/coheng(2) ad3 = (pi/6.d0)*dec*rop*p3*((diam_aer(idmineral(1,3))*1d+2)**3.d0)/coheng(3) endif ft1=ft1+ad1 ft2=ft2+ad2 ft3=ft3+ad3 enddo end subroutine ag2001 !------------------------------------------------------------------------ subroutine calc_z0s_effsoil ! DRAG PARTITION CALCULATION, FACTOR FEFF ranging [0:1] ! Parametrization after [Marticorena and Bergametti, JGR, 1995] ! ! z0s must be calculated only for coarse sand, fine-medium sand and silt ! because clay is too fine and does not represent erodible material. ! ! xeff is a constant (mean distance between non erodible elements), ! equal to 10cm (MB95) but may be 12250 cm (MacKinnon, 2004) ! aeff is a constant: (0.35 in MB95 p.16,420) but 0.7 in other studies ! mineral_dp(): in [cm] ! surf_z0s(): in [cm] implicit none integer :: nz,nm,i,j real(kind=iprec), parameter :: aeff=0.7d0 real(kind=iprec), parameter :: xeff=12.25d0 real(kind=iprec) :: z0s_num, z0s_den, diam_d real(kind=iprec) :: surf_z0s(nsoils) do i=1,nsoils diam_d=dzero do j=1,nmineral if(mineral_frac(i,j).gt.dzero .and. mineral_dp(j).gt.diam_d .and. j.le.3) diam_d=mineral_dp(j) enddo surf_z0s(i) = diam_d/30d0 enddo do nm=1,nmerid do nz=1,nzonal do i=1,nlupcell do j=1,nsoils if (surf_z0s(j).eq.0. .or. surf_z01(nz,nm,i,j).eq.0.) then feff(nz,nm,i,j) = dzero else z0s_num=log(surf_z01(nz,nm,i,j)/surf_z0s(j)) z0s_den=log(aeff*(xeff/surf_z0s(j))**0.8d0) feff(nz,nm,i,j) = dun-z0s_num/z0s_den if (feff(nz,nm,i,j) < dzero) feff(nz,nm,i,j) = dzero if (feff(nz,nm,i,j) > dun) feff(nz,nm,i,j) = dun endif enddo enddo enddo enddo end subroutine calc_z0s_effsoil !------------------------------------------------------------------------ subroutine calc_dust_erod(nz,nm,i,j,merod,cerod) !----------------------- ! Value of the erodibility % per cell, depending on the vegetation ! Landuse type from z0 file : ! USGS landuse: "erodible areas" are limited to Cropland -> lsurf=2 to 7 ! Shrubs -> lsurf=8 and 9 ! Barren -> lsurf=19 ! In all cases, the water mask is applied. ! merod: MODIS erodibility value, not dependent on soil and landuse. ! Output: cerod, erodibility factor between 0 and 1 !----------------------- implicit none integer :: nz, nm, i, j, lsurf real(kind=iprec) :: cerod, merod, nogreen nogreen=(1d0-greenfrac(nz,nm,imonth(ihourrun))) cerod=0d0 lsurf=surf_code(nz,nm,i,j) if (lsurf.ge.2.and.lsurf.le.9)cerod=nogreen if (lsurf.eq.19)cerod=nogreen end subroutine calc_dust_erod !------------------------------------------------------------------------ subroutine check_water(nz,nm,itestwater) ! Check if there is water in the cell. ! To avoid unrealistic fluxes along the coast ! (high mean wind speed over low percentage of erobible surface) ! There is no flux ig the % of water surface is higher than minsurf_th ! ! The indices of water in 'landuse' and 'soil' data are here hardcoded for the USGS database integer :: isoil, ilupc, nz, nm, itestwater, lsurf itestwater=0 ! First check in the soil relative surfaces if(SOILrelsurf(nz,nm,isoilwater).gt.minsurf_th)itestwater=1 if(landusef(nz,nm,iswater).gt.minsurf_th)itestwater=1 !if(landusef(nz,nm,islake).gt.minsurf_th)itestwater=1 ! Second check in the landuse database ! surf_code(nz,nm,i,j): code of the landuse for a specific soil in a cell ! ratiol: % of the landuse for a specific soil in a cell if(itestwater.ne.1)then do isoil=1,nsoils do ilupc=1,nlupcell lsurf=surf_code(nz,nm,ilupc,isoil) if(lsurf.eq.iswater.or.lsurf.ge.24)then if(ratiol(nz,nm,ilupc,isoil).gt.minsurf_th)then itestwater=1 goto 100 endif endif enddo enddo endif 100 continue end subroutine check_water !------------------------------------------------------------------------------------ end module worker_dust