#define TRACEBACK write (gol,'("in ",a," (",a,i6,")")') 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 ! #include "tm5.inc" #include "output.inc" ! !------------------------------------------------------------------------------ ! TM5 ! !------------------------------------------------------------------------------ !BOP ! ! !MODULE: SEDIMENTATION ! ! !DESCRIPTION: This module calculates sedimentation of aerosol particles. ! Also the deposition at the surface is calculated here since ! it uses similar routines. ! ! It is assumed that the distribution is described by nmodes log-normal ! distributions as described in chem_param.F90. ! ! Each mode has a number and mass and a sigma_lognormal. Number and mass are ! related and the mean aerosol radius can thus be calculated for each mode. ! ! mass=pi*4./3.*radius**3.*number*exp(9./2.*ln2s) /1e18*density, with: ! ! density = density of aerosol type ! ln2s = (alog(sigma_lognormal(mode)))**2 ! mass = mode mass ! number = mode number ! ! mode1 : accumulation ! mode2 : coarse ! mode3 : super coarse (ss) coarse ! ! For each mode a separate fall velocity is calculated according to the mass ! and number mean radii. Water take-up by seasalt particles is taken into ! account. This changes the density, radius, and sigma of the distribution. ! ! Also included is the deposition calculation. based on a lookup table ! calculated for a reference aerosol density (e.g. 1800 kg/m3) and a number of ! radii. This deposition curve is convoluted with the number/volume ! distribution of the aerosols. ! ! Again, for SS the water takeup is accounted for, and the effects on density, ! sigma and radius are calculated. The density has effect on the impaction ! term is the depotion calculation. This can be modeled by a shift in the ! radius. Thus the radii of the lookup table are adapted for density ! differences when impaction becomes important. ! ! routine Sedimentation_Init : initialisation of sedimentation routine ! routine Sedimentation_Done : free memory ! routine sedimentation_calcv : calculate velocities ! routine deposition_calcv : calculate velocities ! routine Sedimentation_Apply : do the accutual sedimentation !\\ !\\ ! !INTERFACE: ! MODULE SEDIMENTATION ! ! !USES: ! use GO, only: gol, goPr, goErr, goBug use GO, only: GO_Timer_Def, GO_Timer_End, GO_Timer_Start use Binas, only: grav, xmair, Avog use dims, only: nregions use tm5_distgrid, only: dgrid, Get_DistGrid, gather use partools, only: isRoot use global_types, only: d3_data, emis_data, aerosol_emis_data use chem_param, only: ntracet, imsa, inh4, ino3_a #ifdef with_m7 use chem_param, only: inus_n, iso4nus, isoanus, iais_n, iso4ais, ibcais, ipomais, isoaais use chem_param, only: iacs_n, iso4acs, ibcacs, ipomacs, issacs, iduacs, isoaacs use chem_param, only: icos_n, iso4cos, ibccos, ipomcos, isscos, iducos, isoacos use chem_param, only: iaii_n, ibcaii, ipomaii, isoaaii, iaci_n, iduaci, icoi_n, iducoi #else use chem_param, only: iso4 #endif use chem_param, only: nmod #ifdef with_budgets use budget_global, only: nbudg, nbud_vg, nzon_vg #endif IMPLICIT NONE PRIVATE ! ! !PUBLIC MEMBER FUNCTIONS: ! public Sedimentation_Init, Sedimentation_Done public Sedimentation_Apply public calculate_rh public aerosol_radius ! ! !PUBLIC DATA MEMBERS: ! real, dimension(nregions), public :: sum_drydep, sum_sedimentation ! budgets for tracer #1 (significant on root only) type(d3_data), dimension(nregions), public :: rh, rha ! rh (0-1), cloudfree vs. allsky #ifdef with_m7 integer, parameter, public :: nsed=33 ! number of tracers for which sedimentation is applied integer, parameter, dimension(nsed), public :: ised = (/ inh4, ino3_a, imsa, & inus_n, iso4nus, isoanus, iais_n, iso4ais, ibcais, ipomais, isoaais, & iacs_n, iso4acs, ibcacs, ipomacs, issacs, iduacs, isoaacs, & icos_n, iso4cos, ibccos, ipomcos, isscos, iducos, isoacos, & iaii_n, ibcaii, ipomaii, isoaaii, iaci_n, iduaci, icoi_n, & iducoi /) #else ! sedimentation of few gas-phase species integer, parameter, public :: nsed= 4 ! number of tracers for which sedimentation is applied integer, parameter, dimension(nsed), public :: ised = (/ iso4, inh4, ino3_a, imsa /) #endif ! TB This was made public, to allow getting indices in AerChemMIP output integer, dimension(ntracet),public :: sindex ! pointer to match ntracet array with nsed ! ! !PRIVATE DATA MEMBERS: ! ! Mean values at the surface, per mode (i,j,mode) type(aerosol_emis_data), dimension(nregions) :: vn_deposition_mean ! number mean deposition velocity (m/s) type(aerosol_emis_data), dimension(nregions) :: vn_sedimentation_mean ! number mean sedim. velocity at surface (m/s) type(aerosol_emis_data), dimension(nregions) :: vm_deposition_mean ! mass mean deposition velocity (m/s) type(aerosol_emis_data), dimension(nregions) :: vm_sedimentation_mean ! mass mean sedim. velocity at surface (m/s) type(d3_data), dimension(nregions,nmod) :: vn_sedimentation ! number sedimentation velocities (Pa/s) type(emis_data), dimension(nregions,nmod) :: vn_deposition ! number deposition velocity (Pa/s) type(d3_data), dimension(nregions,nmod) :: vm_sedimentation ! mass sedimentation velocities (Pa/s) type(emis_data), dimension(nregions,nmod) :: vm_deposition ! mass deposition velocity (Pa/s) integer, dimension(nregions) :: n_deposition_mean ! counter integer, dimension(nregions) :: n_sedimentation_mean ! counter logical, parameter :: vd_in_deposition = .false. ! is sedimentation treated in deposition routine? integer, parameter :: ndp = 3 ! limit of the velocity to ndp times the layer thickness ! iteration will account for a fall length through multiple layers real, parameter :: m_to_pa = 7.24e16*grav*xmair*1e3/Avog !factor from m/s --> Pa/s ! ! !PUBLIC TYPES: ! #ifdef with_budgets ! budget terms on grid basis: type, public :: buddep_data real, dimension(:,:,:,:), pointer :: sed ! im,jm,nbud_vg,nsed end type buddep_data type(buddep_data),dimension(nregions),target,public :: buddep_m7_dat #endif integer :: itim_appl real :: sigma, density character(len=*), parameter :: mname = 'sedimentation' ! ! !REVISION HISTORY: ! Feb 2004 by MK - ! 2 Apr 2010 - P. Le Sager - Optimize arrays de/allocation if m7 is not used. ! 19 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! ! !REMARKS: ! (1) Only relative humidity is used when not using m7. ! (2) sedimentation is defined at the bottom of the gridbox (see init) ! (3) Note the switch from emiss_data to aerosol_emiss_data type for the ! *_mean variables (simplify and speedup mpi comm in sedimentation_done) ! ! !TODO: ! !EOP !------------------------------------------------------------------------------ CONTAINS !------------------------------------------------------------------------------ ! TM5 ! !------------------------------------------------------------------------------ !BOP ! ! !IROUTINE: SEDIMENTATION_INIT ! ! !DESCRIPTION: allocate memory, map ised !\\ !\\ ! !INTERFACE: ! SUBROUTINE SEDIMENTATION_INIT ( status ) ! ! !USES: ! use dims, only : lm, iglbsfc use chem_param, only : ntracet use meteodata, only : Set, spm_dat, temper_dat, humid_dat, cc_dat, phlb_dat use meteodata, only : lsp_dat, cp_dat ! ! !OUTPUT PARAMETERS: ! integer,intent(out) :: status ! ! !REVISION HISTORY: ! 2 Apr 2010 - P. Le Sager - Added test on "with_m7" ! 20 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Sedimentation_Init' integer :: region, lmr, i integer :: imode, n, i1, i2, j1, j2 ! --- begin --------------------------------- ! meteo to be used do region = 1, nregions call Set( temper_dat(region), status, used=.true. ) call Set( phlb_dat(region), status, used=.true. ) call Set( humid_dat(region), status, used=.true. ) call Set( cc_dat(region), status, used=.true. ) call Set( lsp_dat(region), status, used=.true. ) call Set( cp_dat(region), status, used=.true. ) end do #ifdef with_m7 ! Calculate mapping of ntracer to sedimentation arrays: sindex(:) = 0 do i=1,nsed do n=1,ntracet if(ised(i) == n) sindex(n) = i enddo enddo #endif ! Allocate do region = 1, nregions call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) lmr=lm(region) #ifdef with_m7 #ifdef with_budgets sum_sedimentation(region) = 0.0 sum_drydep(region) = 0.0 allocate( buddep_m7_dat(region)%sed(i1:i2, j1:j2, nbud_vg, nsed) ) buddep_m7_dat(region)%sed = 0.0 #endif allocate(vn_deposition_mean (region)%surf(i1:i2, j1:j2, nmod)) ! defined at bottom of surface layer allocate(vn_sedimentation_mean(region)%surf(i1:i2, j1:j2, nmod)) ! defined at bottom of surface layer allocate(vm_sedimentation_mean(region)%surf(i1:i2, j1:j2, nmod)) ! defined at bottom of surface layer allocate(vm_deposition_mean (region)%surf(i1:i2, j1:j2, nmod)) ! defined at bottom of surface layer vn_deposition_mean (region)%surf = 0.0 vn_sedimentation_mean(region)%surf = 0.0 vm_sedimentation_mean(region)%surf = 0.0 vm_deposition_mean (region)%surf = 0.0 do imode = 1, nmod allocate(vn_sedimentation(region,imode)%d3 (i1:i2, j1:j2, lmr)) ! defined at bottom of layer allocate(vn_deposition (region,imode)%surf(i1:i2, j1:j2 )) ! defined at bottom of surface layer allocate(vm_sedimentation(region,imode)%d3 (i1:i2, j1:j2, lmr)) ! defined at bottom of layer allocate(vm_deposition (region,imode)%surf(i1:i2, j1:j2 )) ! defined at bottom of surface layer vn_sedimentation(region,imode)%d3 = 0.0 vn_deposition (region,imode)%surf = 0.0 vm_sedimentation(region,imode)%d3 = 0.0 vm_deposition (region,imode)%surf = 0.0 enddo #endif allocate( rh(region)%d3(i1:i2, j1:j2, lmr)) allocate(rha(region)%d3(i1:i2, j1:j2, lmr)) n_deposition_mean(region) = 0 n_sedimentation_mean(region) = 0 enddo call GO_Timer_Def( itim_appl, 'sedimentation appl', status ) IF_NOTOK_RETURN(status=1) status = 0 END SUBROUTINE SEDIMENTATION_INIT !EOC !------------------------------------------------------------------------------ ! TM5 ! !------------------------------------------------------------------------------ !BOP ! ! !IROUTINE: SEDIMENTATION_DONE ! ! !DESCRIPTION: gather/aggregate/write budgets, and free memory !\\ !\\ ! !INTERFACE: ! SUBROUTINE SEDIMENTATION_DONE( status ) ! ! !USES: ! use dims, only : im, jm #ifdef with_budgets use budget_global, only : budget_file_global, nbud_vg, budg_dat, nbudg, NHAB #ifdef with_hdf4 use file_hdf, only : THdfFile, TSds use file_hdf, only : Init, Done, WriteAttribute, WriteData, SetDim #endif use Dims, only : region_name use partools, only : par_reduce_element #endif ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Sedimentation_Done' Integer :: region, i1, i2, j1, j2 #ifdef with_m7 #ifdef with_budgets #ifdef with_hdf4 type(THdfFile) :: io type(TSds) :: sds #endif integer :: j, i, n, nzone, nzone_v real, dimension(nbudg, nbud_vg, nsed) :: budsedg real, dimension(nbudg, nbud_vg, nsed) :: budsedg_all ! for buggy MPI real, dimension(:,:), allocatable :: budget2d ! to gather budget terms on PEs real, dimension(:,:,:), allocatable :: budget3d real, dimension(:,:,:,:),allocatable :: budget4d #endif #endif integer :: imode ! --- begin ------------------------------ #ifdef with_m7 #ifdef with_budgets #ifdef with_hdf4 ! open file if( isRoot ) then call Init(io, budget_file_global, 'write', status) IF_NOTOK_RETURN(status=1) call WriteAttribute(io,'sum_drydep_m7',sum_drydep,status) IF_NOTOK_RETURN(status=1) call WriteAttribute(io,'sum_sedimentation_m7',sum_sedimentation,status) IF_NOTOK_RETURN(status=1) call WriteAttribute(io,'nsed',nsed,status) IF_NOTOK_RETURN(status=1) call WriteAttribute(io,'ised',ised,status) IF_NOTOK_RETURN(status=1) endif #endif budsedg(:,:,:) = 0.0 REG: do region=1,nregions call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) !-- bin aggregated sed budget do n=1,nsed do j=j1,j2 do i=i1,i2 nzone = budg_dat(region)%nzong(i,j) budsedg(nzone,:,n) = budsedg(nzone,:,n) + buddep_m7_dat(region)%sed(i,j,:,n) end do end do end do !-- write Non-Horizontally-Aggregated-Budgets if (NHAB) then ! gather sed budget and write to file if (isRoot )then allocate(budget4d(im(region), jm(region), nbud_vg, nsed)) else allocate(budget4d(1,1,1,1)) end if call GATHER( dgrid(region), buddep_m7_dat(region)%sed, budget4d, 0, status) IF_NOTOK_RETURN(status=1) #ifdef with_hdf4 if(isRoot )then call Init(Sds,io, 'buddep_dat_sed_'//region_name(region),(/im(region),jm(region),nbud_vg,nsed/), 'real(4)', status) IF_NOTOK_RETURN(status=1) call WriteAttribute(Sds,'ised',ised,status) IF_NOTOK_RETURN(status=1) call SetDim(Sds, 0, 'im_'//region_name(region),'longitude', (/(j, j=1,im(region))/), status) call SetDim(Sds, 1, 'jm_'//region_name(region),'latitude', (/(j, j=1,jm(region))/), status) call SetDim(Sds, 2, 'nbud_vg','vertical layer', (/(j, j=1,nbud_vg)/), status) call SetDim(Sds, 3, 'nsed','sedimentation', (/(j, j=1,nsed)/), status) call WriteData(Sds, budget4d ,status) IF_NOTOK_RETURN(status=1) call Done(Sds, status) IF_NOTOK_RETURN(status=1) endif #endif deallocate(budget4d) ! gather deposition velocities and write to file if (isRoot ) then allocate(budget3d(im(region), jm(region), nmod)) else allocate(budget3d(1,1,1)) end if call GATHER( dgrid(region), vn_deposition_mean(region)%surf, budget3d, 0, status) IF_NOTOK_RETURN(status=1) if(isRoot ) then if (n_deposition_mean(region)==0)then budget3d = 0. else budget3d = budget3d / n_deposition_mean(region) end if #ifdef with_hdf4 call Init(Sds,io, 'vd_n_'//region_name(region),(/im(region),jm(region),nmod/), 'real(4)', status) IF_NOTOK_RETURN(status=1) call SetDim(Sds, 0, 'im_'//region_name(region),'longitude', (/(j, j=1,im(region))/), status) IF_NOTOK_RETURN(status=1) call SetDim(Sds, 1, 'jm_'//region_name(region),'latitude', (/(j, j=1,jm(region))/), status) IF_NOTOK_RETURN(status=1) call SetDim(Sds, 2, 'nmodes','deposition modes', (/(j, j=1,nmod)/), status) IF_NOTOK_RETURN(status=1) call WriteData(Sds, budget3d ,status) IF_NOTOK_RETURN(status=1) call Done(Sds, status) IF_NOTOK_RETURN(status=1) #endif endif ! gather and write deposition velocities (mass tracers) call GATHER( dgrid(region), vm_deposition_mean(region)%surf, budget3d, 0, status) IF_NOTOK_RETURN(status=1) if(isRoot ) then if (n_deposition_mean(region)==0)then budget3d = 0. else budget3d = budget3d / n_deposition_mean(region) end if #ifdef with_hdf4 call Init(Sds,io, 'vd_m_'//region_name(region),(/im(region),jm(region),nmod/), 'real(4)', status) IF_NOTOK_RETURN(status=1) call SetDim(Sds, 0, 'im_'//region_name(region),'longitude', (/(j, j=1,im(region))/), status) IF_NOTOK_RETURN(status=1) call SetDim(Sds, 1, 'jm_'//region_name(region),'latitude', (/(j, j=1,jm(region))/), status) IF_NOTOK_RETURN(status=1) call SetDim(Sds, 2, 'nmodes','deposition modes', (/(j, j=1,nmod)/), status) IF_NOTOK_RETURN(status=1) call WriteData(Sds, budget3d ,status) IF_NOTOK_RETURN(status=1) call Done(Sds, status) IF_NOTOK_RETURN(status=1) #endif endif ! gather and write sedimentation velocities call GATHER( dgrid(region), vn_sedimentation_mean(region)%surf, budget3d, 0, status) IF_NOTOK_RETURN(status=1) if(isRoot ) then if (n_sedimentation_mean(region)==0)then budget3d = 0. else budget3d = budget3d / n_sedimentation_mean(region) end if #ifdef with_hdf4 call Init(Sds,io, 'vs_n_'//region_name(region),(/im(region),jm(region),nmod/), 'real(4)', status) IF_NOTOK_RETURN(status=1) call SetDim(Sds, 0, 'im_'//region_name(region),'longitude', (/(j, j=1,im(region))/), status) IF_NOTOK_RETURN(status=1) call SetDim(Sds, 1, 'jm_'//region_name(region),'latitude', (/(j, j=1,jm(region))/), status) IF_NOTOK_RETURN(status=1) call SetDim(Sds, 2, 'nmodes','deposition modes', (/(j, j=1,nmod)/), status) IF_NOTOK_RETURN(status=1) call WriteData(Sds, budget3d ,status) IF_NOTOK_RETURN(status=1) call Done(Sds, status) IF_NOTOK_RETURN(status=1) #endif endif ! gather and write sedimentation velocities (mass tracers) call GATHER( dgrid(region), vm_sedimentation_mean(region)%surf, budget3d, 0, status) IF_NOTOK_RETURN(status=1) if(isRoot) then if (n_sedimentation_mean(region)==0)then budget3d = 0. else budget3d = budget3d / n_sedimentation_mean(region) end if #ifdef with_hdf4 call Init(Sds,io, 'vs_m_'//region_name(region),(/im(region),jm(region),nmod/), 'real(4)', status) IF_NOTOK_RETURN(status=1) call SetDim(Sds, 0, 'im_'//region_name(region),'longitude', (/(j, j=1,im(region))/), status) IF_NOTOK_RETURN(status=1) call SetDim(Sds, 1, 'jm_'//region_name(region),'latitude', (/(j, j=1,jm(region))/), status) IF_NOTOK_RETURN(status=1) call SetDim(Sds, 2, 'nmodes','deposition modes', (/(j, j=1,nmod)/), status) IF_NOTOK_RETURN(status=1) call WriteData(Sds, budget3d ,status) IF_NOTOK_RETURN(status=1) call Done(Sds, status) IF_NOTOK_RETURN(status=1) #endif endif deallocate(budget3d) endif ! NHAB enddo REG ! regions !------- Gather and write aggregated budget ! Sum up contribution from each proc into root array call PAR_REDUCE_ELEMENT( budsedg, 'SUM', budsedg_all, status ) IF_NOTOK_RETURN(status=1) #ifdef with_hdf4 if(isRoot )then call Init(Sds, io, 'budsed_m7', (/nbudg,nbud_vg,nsed/), 'real(8)', status) IF_NOTOK_RETURN(status=1) call SetDim(Sds, 0, 'nbudg', 'horizontal region', (/(j, j=1,nbudg)/), status) call SetDim(Sds, 1, 'nbud_vg', 'vertical layer', (/(j, j=1,nbud_vg)/), status) call SetDim(Sds, 2, 'nsed', 'sedimentation', (/(j, j=1,nsed)/), status) IF_NOTOK_RETURN(status=1) call WriteData(Sds, budsedg_all, status) IF_NOTOK_RETURN(status=1) call Done(Sds, status) IF_NOTOK_RETURN(status=1) call Done(io, status) IF_NOTOK_RETURN(status=1) endif #endif #endif /* BUDGET */ #endif /* M7 */ !-------- Done do region = 1, nregions #ifdef with_m7 #ifdef with_budgets deallocate( buddep_m7_dat(region)%sed) #endif deallocate( vn_deposition_mean (region)%surf ) deallocate( vn_sedimentation_mean(region)%surf ) deallocate( vm_deposition_mean (region)%surf ) deallocate( vm_sedimentation_mean(region)%surf ) do imode = 1, nmod deallocate(vn_sedimentation(region,imode)%d3 ) deallocate(vn_deposition (region,imode)%surf) deallocate(vm_sedimentation(region,imode)%d3 ) deallocate(vm_deposition (region,imode)%surf) enddo #endif deallocate(rh(region)%d3) deallocate(rha(region)%d3) enddo status = 0 END SUBROUTINE SEDIMENTATION_DONE !EOC !------------------------------------------------------------------------------ ! TM5 ! !------------------------------------------------------------------------------ !BOP ! ! !IROUTINE: DEPOSITION_CALCV ! ! !DESCRIPTION: calculate velocities !\\ !\\ ! !INTERFACE: ! SUBROUTINE DEPOSITION_CALCV( region, status ) ! ! !USES: ! use global_data, only : mass_dat, region_dat use meteodata, only : spm_dat, temper_dat, cc_dat use chem_param, only : sigma_lognormal, density_ref use chem_param, only : mode_start, mode_end, names use dims, only : at, bt, nsrce, tref, im, jm, lm, okdebug use chem_param, only : nrdep, lur use dry_deposition, only : vd=>vda use binas, only : pi use partools, only : par_reduce #ifdef with_m7 use m7_data, only : dens_mode, rw_mode ! density (kg/m3) r (m) #endif ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 2 Apr 2010 - P. Le Sager - ! 21 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! ! !REMARKS: ! (1) Called only if m7 used. ! !EOP !------------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/DEPOSITION_CALCV' real, dimension(:,:,:), pointer :: p ! surface pressure (Pa) real, dimension(:,:,:), pointer :: t ! temperature (K) real, dimension(:,:,:), pointer :: cc ! cloud_cover (0-1) real, dimension(:,:,:,:), pointer :: rm ! tracer array (kg) real :: pb ! pressure at bottom of box (Pa) real :: dt ! general timestep (s) real :: dp ! pressure difference layer integer :: i,j,l, mode, n, proc, itn, itracer real :: temp ! temperature real :: to_pascal real, dimension(nmod) :: lns real :: radius real, dimension(nmod,2) :: vd_mean_all, vd_max_all real, dimension(nmod,2) :: vd_mean, vd_max integer, dimension(nmod,2) :: nvd, nvd_all real, dimension(nmod) :: r_mean_all, r_max_all real, dimension(nmod) :: r_mean, r_max integer, dimension(nmod) :: nr, nr_all real, dimension(nrdep) :: d_aer ! diameter vd loopup table (um) real, dimension(nrdep) :: nnumb,nvolume ! number and volume distribution real, dimension(nrdep) :: vdi ! for the integration real, dimension(nrdep) :: vdi_def ! for the integration real, dimension(nrdep) :: lure ! effective loo real, dimension(nrdep) :: ddpi ! integration bin-sizes real, dimension(nrdep+1) :: dlogdp,ddp ! integration edges real, dimension(nrdep) :: logdp,logde ! log(diameter) real :: dpg, ntot, vt real :: lnsigma, rho_aer, sigma integer :: ir, ir1, nshift, nrd, nglob integer :: i1, i2, j1, j2 !________________________start_________________________________ nullify(rm) nullify(p) nullify(t) nullify(cc) rm => mass_dat (region)%rm p => spm_dat (region)%data t => temper_dat(region)%data cc => cc_dat (region)%data dt = float(nsrce)/(2*tref(region)) ! timestep do mode =1,nmod lns(mode) = log(sigma_lognormal(mode)) enddo ! calculate the binsizes (um) around the radii of the pre-calculated vd's d_aer = 2.0*lur ! diameter (um) logdp = log10(d_aer) ! log(diameter) n_deposition_mean(region) = n_deposition_mean(region) + 1 call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) if (okdebug) then if (isRoot) then write(gol,*) '________________________________________________________________________________' ; call goPr write(gol,*) 'Statistics for the deposition routine:' ; call goPr write(gol,*) '________________________________________________________________________________' ; call goPr write(gol,*) 'level mode r_mean r_max(um) vd_n(cm/s) vd_max(cm/s) vd_m(cm/s) vd_mmax(cm/s) ' ; call goPr write(gol,*) '________________________________________________________________________________' ; call goPr end if vd_mean = 0.0 r_mean = 0.0 vd_max = 0.0 r_max = 0.0 nvd = 0 nr = 0 endif !$OMP PARALLEL & !$OMP default (none) & !$OMP shared (region, t, at, bt, p, dt, vd, & !$OMP rm, rh, logdp, & !$OMP vn_deposition_mean, vn_deposition, vm_deposition_mean, vm_deposition, & !$OMP lur, sigma_lognormal, mode_start, mode_end) & !$OMP shared (okdebug, vd_mean, vd_max, nvd, r_max, r_mean, nr) & #ifdef with_m7 !$OMP shared (rw_mode, dens_mode) & #endif !$OMP private (i, j, l, temp, pb, dp, to_pascal, vdi_def, & !$OMP itn, radius, vdi, sigma, lnsigma, density, & !$OMP lure, dlogdp, logde, nshift, nrd, & !$OMP ddp, ddpi, d_aer, dpg, ntot, nnumb, nvolume, vt) l = 1 do j=j1,j2 do i=i1,i2 temp = t(i,j,1) ! at surface to temp box pb = at(l) + bt(l)*p(i,j,1) ! pressure at bottom of box (Pa) dp = (at(l)-at(l+1)) + p(i,j,1)*(bt(l)-bt(l+1)) ! layer thickness to_pascal = m_to_pa*dt*pb/temp ! convert from m/s ---> Pa/timestep do n=1,nrdep vdi_def(n) = vd(region,n)%surf(i,j) enddo #ifdef with_m7 M7MODES: do mode = 1, nmod itn = mode_start(mode) ! position of number tracer ! compute radius: radius = rw_mode(region,mode)%d3(i,j,1) ! initial deposition velocities for increasing radia: vdi = vdi_def sigma = sigma_lognormal(mode) lnsigma = log(sigma) density = dens_mode(region,mode)%d3(i,j,1) if(okdebug) then if(radius > tiny(radius)) then r_mean(mode) = r_mean(mode) + radius nr(mode) = nr(mode) + 1 r_max(mode) = max(r_max(mode), radius) endif endif RADENS: if (radius > 1e-11 .and. density > 1e-2) then ! account for density different than density_ref of the look-up table (lur --> vdi): lure(:) = lur(:) logde(:) = logdp(:) do ir = 2, nrdep if(vdi(ir) > vdi(ir-1)) exit ! for bigger r's : impaction dominates (density effects) if ( ir == nrdep ) exit ! trap upper boundary enddo do ir1 = ir, nrdep lure(ir1) = lur(ir1)*sqrt(density_ref/density) logde(ir1) = log10(2*lure(ir1)) enddo ! compress look-up table such that radii are increasing monotonic: nshift = 0 ir1 = ir do if ( logde(ir1) > logde(ir-1) ) exit nshift = nshift + 1 ir = ir -1 if(ir == 1) exit enddo nrd = nrdep - nshift if (nshift > 0) then do ir1 = ir, nrd logde(ir1) = logde(ir1+nshift) lure(ir1) = lure(ir1+nshift) vdi(ir1) = vdi(ir1+nshift) enddo endif ! do the integration of the shifted lookup table: dlogdp(1) = -3.0 do ir=2,nrd dlogdp(ir) = 0.5*(logde(ir-1)+logde(ir)) ! take middle of the log scale enddo dlogdp(nrd+1) = 3.0 ! 1000 um ddp(1:nrd+1) = 10**dlogdp(1:nrd+1) ddpi(1:nrd) = ddp(2:nrd+1)-ddp(1:nrd) ! integration intervals (um) d_aer(1:nrd) = 2.0*lure(1:nrd) ! perform convolution with log-normal distribution: dpg = 2*radius*1e6 ! diameter in um ntot = rm(i,j,1,itn) ! calculate the distribution (number and mass) over the deposition bins: if(ntot > 1.0 .and. radius > tiny(radius) ) then ! you need some aerosol! do n=1,nrd nnumb(n) = ntot/(sqrt(2.*pi)*d_aer(n)*lnsigma)*exp(-(log(d_aer(n))-log(dpg))**2/(2*lnsigma**2)) nvolume(n) = nnumb(n)*(pi/6.0)*d_aer(n)**3 enddo vt = sum(nnumb(1:nrd)*ddpi(1:nrd)*vdi(1:nrd))/sum(nnumb(1:nrd)*ddpi(1:nrd)) else vt = 0.0 endif vn_deposition_mean(region)%surf(i,j,mode) = vn_deposition_mean(region)%surf(i,j, mode) + vt vn_deposition(region,mode)%surf(i,j) = min(to_pascal*vt,ndp*dp) ! in Pa/timestep downwards if(okdebug) then if ( vt > tiny(vt) ) then vd_mean(mode,1) = vd_mean(mode,1) + vt vd_max(mode,1) = max(vd_max(mode,1) , vt) nvd(mode,1) = nvd(mode,1) + 1 endif endif ! for mass: if(ntot > 1.0 .and. radius > tiny(radius) ) then ! you need some aerosol! vt = sum(nvolume(1:nrd)*ddpi(1:nrd)*vdi(1:nrd))/sum(nvolume(1:nrd)*ddpi(1:nrd)) else vt = 0.0 endif vm_deposition_mean(region)%surf(i,j, mode) = vm_deposition_mean(region)%surf(i,j, mode) + vt vm_deposition(region,mode)%surf(i,j) = min(to_pascal*vt,ndp*dp) ! in Pa/timestep downwards if(okdebug) then if ( vt > tiny(vt) ) then vd_mean(mode,2) = vd_mean(mode,2) + vt vd_max(mode,2) = max(vd_max(mode,2) , vt) nvd(mode,2) = nvd(mode,2) + 1 endif endif ! else vm_deposition(region,mode)%surf(i,j) = 0.0 vn_deposition(region,mode)%surf(i,j) = 0.0 endif RADENS end do M7MODES #endif /* M7 */ enddo !i enddo !j !$OMP END PARALLEL if(okdebug) then do mode = 1, nmod call Par_reduce( r_mean(mode) , 'sum', r_mean_all(mode) , status) IF_NOTOK_RETURN(status=1) call Par_reduce( nr(mode) , 'sum', nr_all(mode) , status) IF_NOTOK_RETURN(status=1) call Par_reduce( r_max(mode) , 'max', r_max_all(mode) , status) IF_NOTOK_RETURN(status=1) call Par_reduce( vd_mean(mode,1), 'sum', vd_mean_all(mode,1), status) IF_NOTOK_RETURN(status=1) call Par_reduce( nvd(mode,1) , 'sum', nvd_all(mode,1) , status) IF_NOTOK_RETURN(status=1) call Par_reduce( vd_max(mode,1) , 'max', vd_max_all(mode,1) , status) IF_NOTOK_RETURN(status=1) call Par_reduce( vd_mean(mode,2), 'sum', vd_mean_all(mode,2), status) IF_NOTOK_RETURN(status=1) call Par_reduce( nvd(mode,2) , 'sum', nvd_all(mode,2) , status) IF_NOTOK_RETURN(status=1) call Par_reduce( vd_max(mode,2) , 'max', vd_max_all(mode,2) , status) IF_NOTOK_RETURN(status=1) if (isRoot) then if(nr_all(mode) > 0) then r_mean_all(mode) = r_mean_all(mode)*1e6/nr_all(mode) else r_mean_all(mode) = 0. end if if(nvd_all(mode,1) > 0) then vd_mean_all(mode,1) = 1e2*vd_mean_all(mode,1)/nvd_all(mode,1) else vd_mean_all(mode,1) = 0. end if if(nvd_all(mode,2) > 0) then vd_mean_all(mode,2) = 1e2*vd_mean_all(mode,2)/nvd_all(mode,2) else vd_mean_all(mode,2) = 0. end if print '(i6,i4,2f10.4,2f12.6,2f12.6)', 0, mode, r_mean_all(mode), 1e6*r_max_all(mode), & vd_mean_all(mode,1), 1e2*vd_max_all(mode,1), vd_mean_all(mode,2), 1e2*vd_max_all(mode,2) end if enddo write(gol,*) '________________________________________________________________________________' ; call goPr endif status=0 END SUBROUTINE DEPOSITION_CALCV !EOC !------------------------------------------------------------------------------ ! TM5 ! !------------------------------------------------------------------------------ !BOP ! ! !IROUTINE: SEDIMENTATION_CALCV ! ! !DESCRIPTION: calculate velocities !\\ !\\ ! !INTERFACE: ! SUBROUTINE SEDIMENTATION_CALCV(region, status) ! ! !USES: ! use global_data, only : mass_dat, region_dat use meteodata, only : spm_dat, temper_dat, cc_dat use chem_param, only : sigma_lognormal, names, mode_start, mode_end use dims, only : at, bt, nsrce, tref, im, jm, lm, okdebug #ifdef with_m7 use m7_data, only : rw_mode, dens_mode #endif use partools, only : par_reduce ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 2 Apr 2010 - P. Le Sager - ! 21 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! ! !REMARKS: ! (1) Called only if m7 used. ! !EOP !------------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Sedimentation_calcv' real, dimension(:,:,:), pointer :: p ! surface pressure (Pa) real, dimension(:,:,:), pointer :: t ! temperature (K) real, dimension(:,:,:), pointer :: cc ! cloud cover (0-1) real, dimension(:,:,:), pointer :: q ! humidity (kg/kg) real, dimension(:,:,:,:), pointer :: rm ! tracer array (kg) real :: pb ! pressure at bottom of box (Pa) real :: dt ! general timestep (s) real :: dp ! pressure difference layer integer :: i,j,l, mode, proc, itn, n, nglob real :: temp ! temperature real :: to_pascal real :: slinnfac real :: zvis, rho_air, radius, Dpg, vt real, dimension(nmod,2) :: vd_mean_all, vd_max_all real, dimension(nmod,2) :: vd_mean, vd_max integer, dimension(nmod,2) :: nvd, nvd_all real, dimension(nmod) :: r_mean_all, r_max_all real, dimension(nmod) :: r_mean, r_max integer, dimension(nmod) :: nr, nr_all real :: lnsigma, sigma integer :: i1, j1, i2, j2 !________________________start_________________________________ nullify(rm) nullify(p) nullify(cc) nullify(t) p => spm_dat (region)%data t => temper_dat(region)%data cc => cc_dat (region)%data rm => mass_dat (region)%rm dt = float(nsrce)/(2*tref(region)) ! timestep call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) ! transfer sedimentation vel. from m/s to Pa/s (note define this as positive = falling) ! dp = -rho*g*dz ! rho is calculated from temperature and pressure.... ! 7.24e16*p/T (#/cm3) (pV=nRT--> nR/V = p/T ) ! #/cm3 ---> kg/m3 f = xmair*1e3/Avog ! Thus: rho (kg/m3) = 7.24e16*p/T * xmair*1e3/Avog ! and (dPa) = -rho*g*(vdep)*dt = -7.24e16*grav*xmair*1e3/Avog * (p/T) *dt n_sedimentation_mean(region) = n_sedimentation_mean(region) + 1 if(okdebug.and.isRoot) then write(gol,*) '________________________________________________________________________________' ; call goPr write(gol,*) 'Statistics for the sedimentation routine:' ; call goPr write(gol,*) '________________________________________________________________________________' ; call goPr write(gol,*) 'level mode r_mean r_max(um) vd_n(cm/s) vd_max(cm/s) vd_m(cm/s) vd_mmax(cm/s) ' ; call goPr write(gol,*) '________________________________________________________________________________' ; call goPr endif #ifdef with_m7 !$OMP PARALLEL & !$OMP default (none) & !$OMP shared (region, t, at, bt, p, dt, & !$OMP rm, lm, & !$OMP vn_sedimentation_mean, vn_sedimentation, vm_sedimentation_mean, vm_sedimentation, & !$OMP sigma_lognormal, mode_start, mode_end) & !$OMP shared (okdebug, vd_mean, vd_max, nvd, r_max, r_mean, nr) & !$OMP shared (rw_mode, dens_mode) & !$OMP private (i, j, l, temp, pb, dp, zvis, rho_air, & !$OMP to_pascal, mode, itn, radius, sigma, lnsigma, & !$OMP density, slinnfac, dpg, vt) LEVS: do l=1, lm(region) if(okdebug) then vd_mean = 0.0 r_mean = 0.0 vd_max = 0.0 r_max = 0.0 nvd = 0 nr = 0 endif do j=j1, j2 do i=i1, i2 if(l == 1) then temp = t(i,j,1) ! at surface to temp box else temp = 0.5*(t(i,j,l)+t(i,j,l-1)) ! interpolate to bottom of box endif pb = at(l) + bt(l)*p(i,j,1) ! pressure at bottom of box (Pa) dp = (at(l)-at(l+1)) + p(i,j,1)*(bt(l)-bt(l+1)) ! layer thickness zvis = 0.0000014963*temp*sqrt(temp)/(temp+120.) ! viscocity (kg/ms) rho_air = 7.24e16*pb/temp * xmair*1e3/Avog ! kg/m3 to_pascal = m_to_pa*dt*pb/temp ! convert from m/s ---> Pa/timestep M7MODES: do mode = 1, nmod radius = rw_mode(region,mode)%d3(i,j,l) sigma = sigma_lognormal(mode) lnsigma = log(sigma) density = dens_mode(region,mode)%d3(i,j,l) slinnfac = exp(2.0*lnsigma*lnsigma) if(okdebug) then if(radius > tiny(radius)) then r_mean(mode) = r_mean(mode) + radius nr(mode) = nr(mode) + 1 r_max(mode) = max(r_max(mode), radius) endif endif ! for number: take first mode (Seinfeld, 1986, pg 286) dpg = radius*2.0*exp(lnsigma*lnsigma) ! diameter in m vt = fall1(pb,Dpg,zvis,temp,rho_air,density) if(okdebug) then if ( vt > tiny(vt) ) then vd_mean(mode,1) = vd_mean(mode,1) + vt vd_max(mode,1) = max(vd_max(mode,1) , vt) nvd(mode,1) = nvd(mode,1) + 1 endif endif vn_sedimentation(region,mode)%d3(i,j,l) = min(to_pascal*vt,ndp*dp) ! in Pa/timestep downwards if(l == 1) then vn_sedimentation_mean(region)%surf(i,j,mode) = & vn_sedimentation_mean(region)%surf(i,j,mode) + vt endif ! for mass: the third moment dpg = radius*2.0*exp(3*lnsigma*lnsigma) ! diameter in m vt = fall1(pb,Dpg,zvis,temp,rho_air,density) if(okdebug) then if ( vt > tiny(vt) ) then vd_mean(mode,2) = vd_mean(mode,2) + vt*slinnfac vd_max(mode,2) = max(vd_max(mode,2) , vt*slinnfac) nvd(mode,2) = nvd(mode,2) + 1 endif endif vm_sedimentation(region,mode)%d3(i,j,l) = min(to_pascal*vt*slinnfac,ndp*dp) ! multiply with slinnfac if(l == 1) then vm_sedimentation_mean(region)%surf(i,j,mode) = & vm_sedimentation_mean(region)%surf(i,j,mode) + vt*slinnfac endif enddo M7MODES enddo ! i enddo ! j if(okdebug) then do mode = 1, nmod call Par_reduce( r_mean(mode) , 'sum', r_mean_all(mode) , status) IF_NOTOK_RETURN(status=1) call Par_reduce( nr(mode) , 'sum', nr_all(mode) , status) IF_NOTOK_RETURN(status=1) call Par_reduce( r_max(mode) , 'max', r_max_all(mode) , status) IF_NOTOK_RETURN(status=1) call Par_reduce( vd_mean(mode,1), 'sum', vd_mean_all(mode,1), status) IF_NOTOK_RETURN(status=1) call Par_reduce( nvd(mode,1) , 'sum', nvd_all(mode,1) , status) IF_NOTOK_RETURN(status=1) call Par_reduce( vd_max(mode,1) , 'max', vd_max_all(mode,1) , status) IF_NOTOK_RETURN(status=1) call Par_reduce( vd_mean(mode,2), 'sum', vd_mean_all(mode,2), status) IF_NOTOK_RETURN(status=1) call Par_reduce( nvd(mode,2) , 'sum', nvd_all(mode,2) , status) IF_NOTOK_RETURN(status=1) call Par_reduce( vd_max(mode,2) , 'max', vd_max_all(mode,2) , status) IF_NOTOK_RETURN(status=1) if (isRoot) then if(nr_all(mode) > 0) then r_mean_all(mode) = r_mean_all(mode)*1e6/nr_all(mode) else r_mean_all(mode) = 0. end if if(nvd_all(mode,1) > 0) then vd_mean_all(mode,1) = 1e2*vd_mean_all(mode,1)/nvd_all(mode,1) else vd_mean_all(mode,1) = 0. end if if(nvd_all(mode,2) > 0) then vd_mean_all(mode,2) = 1e2*vd_mean_all(mode,2)/nvd_all(mode,2) else vd_mean_all(mode,2) = 0. end if print '(i6,i4,2f10.4,2f12.6,2f12.6)', l, mode, r_mean_all(mode), 1e6*r_max_all(mode), & vd_mean_all(mode,1), 1e2*vd_max_all(mode,1), vd_mean_all(mode,2), 1e2*vd_max_all(mode,2) end if enddo write(gol,*) '________________________________________________________________________________' ; call goPr endif enddo LEVS !$OMP END PARALLEL #endif /* M7 */ status = 0 END SUBROUTINE SEDIMENTATION_CALCV !EOC !------------------------------------------------------------------------------ ! TM5 ! !------------------------------------------------------------------------------ !BOP ! ! !IROUTINE: SEDIMENTATION_APPLY ! ! !DESCRIPTION: Sedimentation calculation based on pre-calculated ! v_sedimentation. ! Based on dynamw routine, including the slopes. !\\ !\\ ! !INTERFACE: ! SUBROUTINE SEDIMENTATION_APPLY( region, status ) ! ! !USES: ! use global_data, only : mass_dat, region_dat use meteodata , only : spm_dat use chem_param, only : ra use dims, only : lm, at, bt, limits, sec_month, nsrce, tref, okdebug use chem_param, only : mode_nm, mode_tracers, names use partools, only : par_reduce, par_reduce_element #ifdef with_m7 use chem_param, only : mode_tracers_sed, mode_nm_sed use emission_data, only : bb_lm #ifndef without_emission use emission_data, only : emis_mass, emis_number #endif #endif #ifdef with_budgets use emission_data, only : budemi_dat, sum_emission use ebischeme, only : buddep_dat #endif ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 19 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! ! !REMARKS: ! (1) effective only if m7 used... ! !EOP !------------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Sedimentation_Apply' real,dimension(:,:,:,:),pointer :: rm,rxm,rym,rzm real,dimension(:,:,:),pointer :: vsold real,dimension(:,: ),pointer :: vdold real,dimension(:,:,:),pointer :: p ! surface pressure real,dimension(:,:,:),allocatable :: f, pf, fx, fy, vs real,dimension(:,:,:),allocatable :: emit real,dimension(:,:),allocatable :: vd integer :: i, j, l, n, lmr, mode real, parameter :: one = 1.0 real :: gamma, gam, l_gam real :: dp, dtime, month2dt real :: rmold, rmnew, rmplus, fdep, fsed integer,parameter :: iter_limit=200 integer :: n_advz, iter integer :: nzone_v integer :: ls, le, ii, inmode real :: l_sum(3), g_sum(3) integer :: i1, j1, i2, j2, iflag, jflag, lflag !________________________start_________________________________ ! Big if-condition over all routine ! start timing: call GO_Timer_Start( itim_appl, status ) IF_NOTOK_RETURN(status=1) #ifdef with_m7 ! Bounds call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) lmr = lm(region) allocate( f ( i1:i2, j1:j2, lmr )) ! flux of tracer (kg) allocate(pf ( i1:i2, j1:j2, lmr )) ! rxz moment flux allocate(fx ( i1:i2, j1:j2, lmr )) ! rmx flux allocate(fy ( i1:i2, j1:j2, lmr )) ! rmy flux allocate(vs ( i1:i2, j1:j2, lmr )) allocate(vd ( i1:i2, j1:j2 )) allocate(emit( i1:i2, j1:j2, bb_lm)) call sedimentation_calcv(region, status) ! calculate sedimentation.... IF_NOTOK_RETURN(status=1) call deposition_calcv(region, status) ! calculate deposition..... IF_NOTOK_RETURN(status=1) nullify(rm) nullify(rxm) nullify(rym) nullify(rzm) nullify(p) p => spm_dat(region)%data rm => mass_dat(region)%rm rxm => mass_dat(region)%rxm rym => mass_dat(region)%rym rzm => mass_dat(region)%rzm if(okdebug) then write(gol,*) ' call sedimentation'; call goPr end if ! If requested, limit vertical slopes such that no non-negative tracer masses should occur ls = 1 ; le = lmr if (limits) then do n = 1, ubound (rzm, 4) do l = ls, le do j = j1, j2 do i = i1, i2 rzm(i,j,l,n) = max(min(rzm(i,j,l,n),rm(i,j,l,n)), -rm(i,j,l,n)) end do end do end do end do endif ! determine timestep dtime=float(nsrce)/(2*tref(region)) ! conversion to emission per timestep month2dt=dtime/sec_month ! ================= ! Loop over tracers ! ================= do mode =1,nmod !do inmode=0,mode_nm(mode) do inmode=0,mode_nm_sed(mode) n = mode_tracers_sed(inmode,mode) !n = mode_tracers(inmode,mode) !------------- reset n_advz=1 nullify(vsold) nullify(vdold) if (inmode == 0) then ! number or mass tracer vsold => vn_sedimentation(region,mode)%d3 vdold => vn_deposition(region,mode)%surf else vsold => vm_sedimentation(region,mode)%d3 vdold => vm_deposition(region,mode)%surf endif !------------- find # iter needed for no CFL violation cfl: do vs(:,j1:j2,:) = vsold(:,j1:j2,:)/float(n_advz) vd(:,j1:j2) = vdold(:,j1:j2) /float(n_advz) advz: do iter=1,n_advz ! reset gamma l_gam = 0. do l= 1, lmr do j= j1, j2 do i= i1, i2 ! vs (Pa) always downwards ! calculate the flux at the bottom of box: ! note that we define this as negative (vs>0, f<0) ! pressure difference in layer dp (Pa) dp = (at(l)-at(l+1)) + p(i,j,1)*(bt(l)-bt(l+1)) ! positive.. if(l /= 1) then gamma=-vs(i,j,l)/dp !always negative (according to negative cm dynamw) else gamma=-(vs(i,j,l)+vd(i,j))/dp !always negative (according to negative cm dynamw) endif ! replace 'one' with 0.99999 (NB/PLS) if (abs(gamma)>=0.9999999) then l_gam=max(l_gam,abs(gamma)) ! exit advz ! commented for consistency with TM5 end if enddo enddo enddo enddo advz call Par_REDUCE( l_gam, 'MAX', gam, status, all=.true.) ! if CFL violation increase iteration counter and re-start cfl loop ! but check if the number of iterations becomes too large if(abs(gam)>=0.9999999) then ! PLS: replace "one" with 0.9999999 to be consistent with above !Commented since "exit advz" is commented, and hence not meaningfull ! if(okdebug)then ! ! iflag=min(i,i2) ! jflag=min(j,j2) ! lflag=min(l,lmr) ! ! ! i-j-l information correct if from the processor where violation happened ! print *,' --- CFL violation: gamma=',gam,vs(iflag,jflag,lflag),' in (region,i,j,l,n)= ',region,iflag,jflag,lflag,n ! print*,' #iterations:', n_advz ! endif n_advz = n_advz + 1 if ( n_advz > iter_limit ) then status=1 IF_NOTOK_RETURN(status=1) end if cycle cfl else ! situation OK, no CFL: use this n_advz exit cfl endif enddo cfl !------------- Apply number of iterations to calculate new tracer distributions #ifdef with_budgets l_sum = 0.0 #endif #ifndef without_emission if (n==ino3_a.or.n==inh4.or.n==imsa) then emit(:,j1:j2,:) = 0.0 else if(inmode == 0) then emit(:,j1:j2,:) = 0.0 do ii=1,mode_nm(mode) ! add up all number emissions in the mode... emit(:,j1:j2,:) = emit(:,j1:j2,:) + emis_number(region,mode)%d4(:,j1:j2,:,ii)*month2dt/float(n_advz) enddo else ! this is a 'mass' emission with index nmode emit(:,j1:j2,:) = emis_mass(region,mode)%d4(:,j1:j2,:,inmode)*month2dt/float(n_advz) endif #else emit = 0.0 #endif ! do the iteration do iter=1,n_advz ! compute f, pf, fx and fy do l= 1, lmr do j= j1, j2 do i= i1, i2 ! vs (Pa) always downwards ! calculate the flux at the bottom of box: ! note that we define this as negative (vs>0, f<0) ! pressure difference in layer dp (Pa) dp = (at(l)-at(l+1)) + p(i,j,1)*(bt(l)-bt(l+1)) ! positive.. if(l /= 1) then gamma=-vs(i,j,l)/dp !always negative (according to negative cm dynamw) else gamma=-(vs(i,j,l)+vd(i,j))/dp !always negative (according to negative cm dynamw) endif f(i,j,l) = gamma*(rm(i,j,l,n)-(one+gamma)*rzm(i,j,l,n)) !kg (neg.) pf(i,j,l) = -vs(i,j,l)*(gamma*gamma*rzm(i,j,l,n)-3.*f(i,j,l)) !neg. fx(i,j,l) = gamma*rxm(i,j,l,n) !kg (neg.) fy(i,j,l) = gamma*rym(i,j,l,n) !kg (neg.) ! - Cannot happen unless there is a bug (pls) !if (abs(gamma)>=one) then ! status=1 ! IF_NOTOK_RETURN(status=1) !end if enddo enddo enddo ! calculate new tracer mass, and tracer mass slopes ! update rm, rzm, rxm and rym in interior layers of the column l = lmr do j = j1, j2 do i = i1, i2 dp = (at(l)-at(l+1)) + p(i,j,1)*(bt(l)-bt(l+1)) rm (i,j,l,n)=rm(i,j,l,n) + f(i,j,l) rzm(i,j,l,n)=rzm(i,j,l,n)+ & ( pf(i,j,l) & -(-vs(i,j,l)) *rzm(i,j,l,n) & +3.*( ( -vs(i,j,l)) *rm (i,j,l,n) & -( f(i,j,l))*dp ))/dp if(limits) then rzm(i,j,l,n) = max(min(rzm(i,j,l,n), rm(i,j,l,n)), -rm(i,j,l,n)) endif rxm(i,j,l,n)=rxm(i,j,l,n)+(fx(i,j,l)) rym(i,j,l,n)=rym(i,j,l,n)+(fy(i,j,l)) #ifdef with_budgets nzone_v=nzon_vg(l) buddep_m7_dat(region)%sed(i,j,nzone_v,sindex(n)) = & buddep_m7_dat(region)%sed(i,j,nzone_v,sindex(n)) - f(i,j,l)*1e3/ra(n) ! in mole ! Downward flux is negative. And it is defined at the bottom of the box. A positive flux ! is positive and would be an income for the respective grid cell. A minus sign takes care ! We want to define the sedimentation as a cost. #endif enddo enddo do l = lmr-1, 2, -1 do j = j1, j2 do i = i1, i2 dp = (at(l)-at(l+1)) + p(i,j,1)*(bt(l)-bt(l+1)) rmold = rm(i,j,l,n) if (l .le. bb_lm) then rmplus = rmold + f(i,j,l)-f(i,j,l+1) + emit(i,j,l) else rmplus = rmold + f(i,j,l)-f(i,j,l+1) endif rm(i,j,l,n) = rmplus rzm(i,j,l,n)=rzm(i,j,l,n)+ & ( pf(i,j,l)-pf(i,j,l+1) & -(vs(i,j,l+1)-vs(i,j,l)) *rzm(i,j,l,n) & +3.*( (-vs(i,j,l+1)-vs(i,j,l)) *rm (i,j,l,n) & -(f(i,j,l+1)+ f(i,j,l)) *dp ))/dp if(limits) then rzm(i,j,l,n) = max(min(rzm(i,j,l,n), rm(i,j,l,n)), -rm(i,j,l,n)) endif rxm(i,j,l,n)=rxm(i,j,l,n)+(fx(i,j,l)-fx(i,j,l+1)) rym(i,j,l,n)=rym(i,j,l,n)+(fy(i,j,l)-fy(i,j,l+1)) #ifdef with_budgets nzone_v=nzon_vg(l) buddep_m7_dat(region)%sed(i,j,nzone_v,sindex(n)) = & buddep_m7_dat(region)%sed(i,j,nzone_v,sindex(n)) - (f(i,j,l)-f(i,j,l+1))*1e3/ra(n) ! in mole ! The sedimentation is calculated as 'income' again. With the minus sign, those are costs. #ifndef without_emission if ( l <= bb_lm ) then budemi_dat(region)%emi(i,j,nzone_v,n) = & budemi_dat(region)%emi(i,j,nzone_v,n) + emit(i,j,l)*1e3/ra(n) ! in mole end if #endif #endif enddo enddo enddo l = 1 do j = j1, j2 do i = i1, i2 dp = (at(l)-at(l+1)) + p(i,j,1)*(bt(l)-bt(l+1)) if(vd_in_deposition) then write(gol,*)' apply_sedimentation: vs in deposition disabled'; call goBug status=1 IF_NOTOK_RETURN(status=1) end if !if(vd_in_deposition) then ! rm (i,j,l,n)=rm(i,j,l,n) - f(i,j,l+1) ! rzm(i,j,l,n)=rzm(i,j,l,n)+ & ! ( -pf(i,j,l+1) & ! -(vs(i,j,l+1)) *rzm(i,j,l,n) & ! +3.*( (-vs(i,j,l+1)) *rm (i,j,l,n) & ! -(f(i,j,l+1))*dp ))/dp ! if(limits) then ! rzm(i,j,l,n) = max(min(rzm(i,j,l,n), rm(i,j,l,n)), -rm(i,j,l,n)) ! endif ! rxm(i,j,l,n)=rxm(i,j,l,n)+(-fx(i,j,l+1)) ! rym(i,j,l,n)=rym(i,j,l,n)+(-fy(i,j,l+1)) !else ! rm (i,j,l,n)=rm(i,j,l,n) + f(i,j,l)-f(i,j,l+1) ! rzm(i,j,l,n)=rzm(i,j,l,n)+ & ! ( pf(i,j,l)-pf(i,j,l+1) & ! -(vs(i,j,l+1)-vs(i,j,l)) *rzm(i,j,l,n) & ! +3.*( (-vs(i,j,l+1)-vs(i,j,l)) *rm (i,j,l,n) & ! -(f(i,j,l+1)+ f(i,j,l))*dp ))/dp ! if(limits) then ! rzm(i,j,l,n) = max(min(rzm(i,j,l,n), rm(i,j,l,n)), -rm(i,j,l,n)) ! endif ! apply to rxm, rym the fluxes: rxm(i,j,l,n)=rxm(i,j,l,n)+(fx(i,j,l)-fx(i,j,l+1)) rym(i,j,l,n)=rym(i,j,l,n)+(fy(i,j,l)-fy(i,j,l+1)) ! for rm: apply emissions, sedimentation flux from above ! deposition and sedimentation at surface: Backward Eularian rmold = rm(i,j,l,n) rmplus = (rmold - f(i,j,l+1) + emit(i,j,l)) ! note f is negative ! rmnew = rmplus/(1. + (vd(i,j) + vs(i,j,1))/dp) rm(i,j,l,n) = rmnew if(rmold > 1e-8) rzm(i,j,l,n) = rzm(i,j,l,n)*rmnew/rmold ! budget: #ifdef with_budgets if((vs(i,j,1) + vd(i,j)) > 1e-12) then fsed = -(rmplus-rmnew)*vs(i,j,1)/(vs(i,j,1) + vd(i,j)) ! is negative fdep = -(rmplus-rmnew)*vd(i,j) /(vs(i,j,1) + vd(i,j)) else fsed = 0.0 fdep = 0.0 endif if(n == 1) then l_sum(1) = l_sum(1) + f(i,j,l+1) - fsed ! goes into sum_sedimentation l_sum(2) = l_sum(2) - fdep ! goes into sum_drydep #ifndef without_emission l_sum(3) = l_sum(3) + emit(i,j,1) ! goes into sum_emission #endif endif nzone_v=nzon_vg(l) buddep_m7_dat(region)%sed(i,j,nzone_v,sindex(n)) = & buddep_m7_dat(region)%sed(i,j,nzone_v,sindex(n)) - (fsed - f(i,j,l+1) )*1e3/ra(n) ! in mole ! The comment line says that fsed is negative. As we regard sedimentation as a cost, we define ! a positive sedimentation as a loss of material. A negative fsed is loss of material, so we ! need to have -fsed. The (upward) flux on the top of the grid cell is a cost as well. buddep_dat(region)%dry(i,j,n) = & buddep_dat(region)%dry(i,j,n) - fdep*1e3/ra(n) ! in mole #ifndef without_emission budemi_dat(region)%emi(i,j,1,n) = & budemi_dat(region)%emi(i,j,1,n) + emit(i,j,1)*1e3/ra(n) ! in mole #endif #endif /* BUDGETS */ enddo !i enddo !j enddo ! iter ! #ifdef with_budgets call PAR_REDUCE_ELEMENT( l_sum, 'SUM', g_sum, status ) if(isRoot ) then sum_sedimentation(region) = sum_sedimentation(region) + g_sum(1) sum_drydep (region) = sum_drydep (region) + g_sum(2) #ifndef without_emission sum_emission (region) = sum_emission (region) + g_sum(3) #endif end if #endif enddo ! loop over tracers in mode enddo ! loop over modes deallocate( f) deallocate(pf) deallocate(fx) deallocate(fy) deallocate(vs) deallocate(vd) deallocate(emit) nullify(rm) nullify(rxm) nullify(rym) nullify(rzm) nullify(vsold) nullify(vdold) #endif /* M7 */ ! end timing: call GO_Timer_End( itim_appl, status ) IF_NOTOK_RETURN(status=1) status = 0 END SUBROUTINE SEDIMENTATION_APPLY !EOC #ifdef with_m7 !------------------------------------------------------------------------------ ! TM5 ! !------------------------------------------------------------------------------ !BOP ! ! !IROUTINE: FALL1 ! ! !DESCRIPTION: function to calculate the fall velocity from the particle ! diameter, as a function of density, temperature and pressure. !\\ !\\ ! !INTERFACE: ! REAL FUNCTION FALL1( press, zmd, zvis, t, zdensair, densaer_p) result(vt) ! ! !INPUT PARAMETERS: ! real, intent(in) :: press ! pressure in Pa (or bar?) real, intent(in) :: zmd ! median radius of aerosol (m) real, intent(in) :: zvis ! viscosity (kg/(sm)) real, intent(in) :: t ! temperature (K) real, intent(in) :: zdensair ! density air (kg/m3) real, intent(in) :: densaer_p ! density aerosol particles (kg/m3) ! ! !REVISION HISTORY: ! 2 Apr 2010 - P. Le Sager - ! ! !REMARKS: ! Called in Sedimentation_Apply, only if m7 used. ! !EOP !------------------------------------------------------------------------------ !BOC real :: zlair ! ----------- start if (zmd > tiny(zmd)) then vt=2./9.*(densaer_p-zdensair) * grav/zvis*(zmd/2.)**2. zlair=0.066*(1.01325e5/press)*t/293.15*1e-6 !--- With cunnigham slip- flow correction (S&P, Equation 8.34): vt = vt * (1.+ 1.257*zlair/zmd*2. + 0.4*zlair/zmd*2. * exp(-1.1/(zlair/zmd*2.)) ) else vt = 0.0 endif END FUNCTION FALL1 !EOC #endif /* M7 */ !------------------------------------------------------------------------------ ! TM5 ! !------------------------------------------------------------------------------ !BOP ! ! !IROUTINE: AEROSOL_RADIUS ! ! !DESCRIPTION: function that calculates the median aerosol radius (m), ! given the total mass and number of a log-normal distribution. !\\ !\\ ! !INTERFACE: ! REAL FUNCTION AEROSOL_RADIUS(mtot, ntot, sigma, rho_aer) result(radius) ! ! !USES: ! use Binas, only: Pi ! ! !INPUT PARAMETERS: ! real, intent(in) :: mtot ! total mass (kg) real, intent(in) :: ntot ! total number (#) real, intent(in) :: sigma ! the sigma of the log-normal aerosol size distribution real, intent(in) :: rho_aer ! the density of the aerosol (kg/m3) ! ! !REVISION HISTORY: ! 2 Apr 2010 - P. Le Sager - ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------------ !BOC real :: lns if(mtot > tiny(mtot) .and. ntot > tiny(ntot)) then lns = alog(sigma) radius = (mtot/(rho_aer*4.0*pi/3.0 *ntot *exp((9./2.)*lns**2)))**(1./3.) else radius = 0.0 endif END FUNCTION AEROSOL_RADIUS !EOC !------------------------------------------------------------------------------ ! TM5 ! !------------------------------------------------------------------------------ !BOP ! ! !IROUTINE: CALCULATE_RH ! ! !DESCRIPTION: calculate the rh, with 100% rh assumption in cloudy part. ! In the cloud free part, the rh thus is smaller, and the water ! uptake by aerosols will be smaller. !\\ !\\ ! !INTERFACE: ! SUBROUTINE CALCULATE_RH ! ! !USES: ! use chem_param, only : xmh2o use meteodata, only : phlb_dat use MeteoData, only : temper_dat, humid_dat, cc_dat use dims, only : nregions, im, jm ,lm ! ! !REVISION HISTORY: ! 19 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! !EOP !------------------------------------------------------------------------------ !BOC real, dimension(:,:,:), pointer :: phlb ! pressure at border (Pa) real, dimension(:,:,:), pointer :: t ! temperature (K) real, dimension(:,:,:), pointer :: q ! humidity (kg h2o / kg air) real, dimension(:,:,:), pointer :: cc ! cloud cover (0-1) real :: tr, wv, airn, h2on, rrh, s, ccs integer :: region, i,j,l, i1,i2,j1,j2 do region = 1, nregions call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) nullify(phlb) nullify(t) nullify(q) nullify(cc) phlb => phlb_dat(region)%data t => temper_dat(region)%data q => humid_dat(region)%data cc => cc_dat(region)%data do l=1, lm(region) do j=j1, j2 do i=i1, i2 tr = 1. - 373.15/t(i,j,l) wv=exp((((-.1299*tr-.6445)*tr-1.976)*tr+13.3185)*tr) airn = 7.24e16*0.5*(phlb(i,j,l) + phlb(i,j,l+1))/t(i,j,l) ! somethings seem redundent here! h2on = q(i,j,l)*airn*xmair/xmh2o ! leave it for readability! rrh = h2on*t(i,j,l)/(1013.25*wv*7.24e16) s = 0.01*max(0.0, min(rrh, 99.9 ) ) ! 0-0.999 scale! rha(region)%d3(i,j,l) = s ! scale relative humidity to cloudfree part ! assuming 100% rh in the cloudy part, but never smaller than 0.75! if (s > 0.90) then ccs = cc(i,j,l) if((1. - ccs) > tiny(ccs)) s = max(0.75, (s-ccs)/(1. - ccs)) endif rh(region)%d3(i,j,l) = s enddo enddo enddo nullify(phlb) nullify(t) nullify(q) nullify(cc) enddo END SUBROUTINE CALCULATE_RH !EOC END MODULE SEDIMENTATION