#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: CHEMISTRY ! ! !DESCRIPTION: perform CB05 chemistry simulation in TM5 + M7 !\\ !\\ ! !INTERFACE: ! MODULE CHEMISTRY ! ! !USES: ! use GO, only : gol, goErr, goPr, goBug use tm5_distgrid, only : dgrid, Get_DistGrid, gather use dims, only : nregions #ifdef with_budgets use budget_global, only : nbudg, nbud_vg use chem_param, only : ntrace, ntracet #endif IMPLICIT NONE PRIVATE ! ! !PUBLIC MEMBER FUNCTIONS: ! public :: Chemistry_Init, Chemistry_Done public :: Chemie ! ! !PRIVATE DATA MEMBERS: ! #ifdef with_budgets real, dimension(:,:,:), allocatable :: budchem real, Dimension(:,:,:), Allocatable :: budeqsam real :: eminox #ifdef with_m7 real, dimension(:,:,:), allocatable :: budm7proc real, dimension(:,:,:), allocatable :: budm7 real, dimension(:,:,:),allocatable,public::d_nucle, growth1_2, coag_sink_nuc, cond1_soa, cond1_su, m_nucle_su real, dimension(:,:,:),allocatable,public::cond2_soa, cond3_soa, cond4_soa, cond5_soa, m_nucle_soa #endif #endif character(len=*), parameter :: mname = 'chemistry' ! ! minimal concentration required to apply scaling for slopes ! (to avoid floating point overflows): real, parameter :: eps = 1.0e-20 ! ! !REVISION HISTORY: ! 2005 - Elisabetta Vignati - changed for coupling with M7 ! 22 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ CONTAINS !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: CHEMISTRY_INIT ! ! !DESCRIPTION: allocate and initialize budget arrays. Init M7 if needed. !\\ !\\ ! !INTERFACE: ! SUBROUTINE CHEMISTRY_INIT( status ) ! ! !USES: ! use Dims , only : nregions, im, jm, lm #ifdef with_budgets use Chem_param, only : ntrace, ntracet use ebischeme , only : sum_deposition, sum_wetchem, sum_chemistry use ebischeme , only : buddep_dat, budrrg, budrjg, budrwg,diag_prod,AC_diag_prod,AC_O3_lp,AC_loss,nprod,nprod_AC,nprod_AC_o3,n_loss !PLS use chem_param, only : ch4oh, so4pg, so4pa, o3p, o3l use chem_param, only : o3p, o3l #ifdef with_m7 use m7_data, only : nm7procs #endif #endif #ifdef with_m7 use mo_aero_m7, only : m7_initialize use m7_data, only : init_m7_data #endif use meteoData , only : Set use MeteoData , only : phlb_dat, m_dat, temper_dat, humid_dat, gph_dat use MeteoData , only : lwc_dat, iwc_dat, cc_dat ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 22 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! !EOP !------------------------------------------------------------------------ !BOC integer :: region, i1, i2, j1, j2 character(len=*), parameter :: rname = mname//'/Chemistry_Init' ! --- begin -------------------------------- ! select meteo to be used: do region = 1, nregions call Set( phlb_dat(region), status, used=.true. ) call Set( m_dat(region), status, used=.true. ) call Set( temper_dat(region), status, used=.true. ) call Set( humid_dat(region), status, used=.true. ) call Set( gph_dat(region), status, used=.true. ) call Set( lwc_dat(region), status, used=.true. ) call Set( iwc_dat(region), status, used=.true. ) call Set( cc_dat(region), status, used=.true. ) enddo #ifdef with_budgets do region = 1, nregions sum_chemistry (region) = 0.0 sum_deposition(region) = 0.0 sum_wetchem (region) = 0.0 call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) allocate(buddep_dat(region)%dry(i1:i2, j1:j2, ntrace)) buddep_dat(region)%dry = 0.0 !nprod=4 #ifdef with_m7 allocate(diag_prod (region)%prod(i1:i2, j1:j2, lm(region), nprod)) allocate(AC_diag_prod(region)%prod(i1:i2, j1:j2, lm(region), nprod_AC)) allocate(AC_O3_lp (region)%prod(i1:i2, j1:j2, lm(region), nprod_AC_o3)) allocate(AC_loss (region)%prod(i1:i2, j1:j2, lm(region), n_loss)) diag_prod(region)%prod = 0.0 AC_diag_prod(region)%prod = 0.0 AC_O3_lp(region)%prod = 0.0 AC_loss(region)%prod = 0.0 allocate(d_nucle(i1:i2,j1:j2,lm(region))) allocate(m_nucle_su(i1:i2,j1:j2,lm(region))) allocate(m_nucle_soa(i1:i2,j1:j2,lm(region))) allocate(growth1_2(i1:i2,j1:j2,lm(region))) allocate(coag_sink_nuc(i1:i2,j1:j2,lm(region))) allocate(cond1_soa(i1:i2,j1:j2,lm(region))) allocate(cond2_soa(i1:i2,j1:j2,lm(region))) allocate(cond3_soa(i1:i2,j1:j2,lm(region))) allocate(cond4_soa(i1:i2,j1:j2,lm(region))) allocate(cond5_soa(i1:i2,j1:j2,lm(region))) allocate(cond1_su(i1:i2,j1:j2,lm(region))) d_nucle=0.0 m_nucle_soa=0.0 m_nucle_su=0.0 growth1_2=0.0 coag_sink_nuc=0.0 cond1_soa=0.0 cond2_soa=0.0 cond3_soa=0.0 cond4_soa=0.0 cond5_soa=0.0 cond1_su=0.0 #endif allocate(o3p(region)%d3(i1:i2, j1:j2, lm(region))); o3p(region)%d3 = 0.0 allocate(o3l(region)%d3(i1:i2, j1:j2, lm(region))); o3l(region)%d3 = 0.0 ! >>> ! TvN: still in sources_sinks.F90 ! allocate(ch4oh(region)%d3(i1:i2, j1:j2, lm(region))); ch4oh(region)%d3 = 0.0 ! allocate(so4pg(region)%d3(i1:i2, j1:j2, lm(region))); so4pg(region)%d3 = 0.0 ! allocate(so4pa(region)%d3(i1:i2, j1:j2, lm(region))); so4pa(region)%d3 = 0.0 ! <<< enddo budrrg = 0.0 budrjg = 0.0 budrwg = 0.0 Allocate(budchem(nbudg,nbud_vg,ntrace)) budchem = 0.0 ! Two acid-base reactions Allocate(budeqsam(nbudg,nbud_vg,2)) budeqsam = 0.0 #ifdef with_m7 ! M7 Bugets Allocate(budm7proc(nbudg,nbud_vg,nm7procs)) Allocate(budm7(nbudg,nbud_vg,ntracet)) budm7proc = 0.0 budm7 = 0.0 #endif ! global lightning NOx eminox = 0. #endif /* BUDGETS */ #ifdef with_m7 CALL M7_INITIALIZE CALL INIT_M7_DATA( status ) ! declare M7 output data structure IF_NOTOK_RETURN(status=1) #endif ! ok status = 0 END SUBROUTINE CHEMISTRY_INIT !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: CHEMISTRY_DONE ! ! !DESCRIPTION: gather and write budgets !\\ !\\ ! !INTERFACE: ! SUBROUTINE CHEMISTRY_DONE( status ) ! ! !USES: ! use dims, only : nregions, im, jm, lm #ifdef with_budgets use dims, only : region_name #ifdef with_hdf4 use file_hdf, only : Init, Done, WriteAttribute, WriteData, SetDim, THdfFile, TSds #endif use budget_global, only : budget_file_global, budg_dat, NHAB #ifndef without_boundary use boundary, only : budstratg #endif #ifndef without_emission use emission, only : budemig_all #endif use partools, only : isRoot, par_reduce_element, par_reduce use ebischeme, only : sum_deposition, sum_wetchem, sum_chemistry use ebischeme, only : buddep_dat, budrrg, budrjg, budrwg, budmarkg, diag_prod, AC_diag_prod, AC_O3_lp, AC_loss !PLS use chem_param, only : ch4oh, so4pg, so4pa, o3p, o3l, marknam, nmark use chem_param, only : o3p, o3l, marknam, nmark use chem_param, only : marknam, ntrace use chem_param, only : inox use reaction_data, only : nreac, ratnam, nreacw, rwnam use photolysis_data, only : nj, jnam #ifdef with_cariolle use chem_param, only : ncar, carnam use chemistry_cariolle, only : budcarg #endif #ifdef with_m7 use m7_data, only : nm7procs #endif #endif /* BUDGETS */ #ifdef with_m7 use m7_data, only : free_m7_data #endif ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 22 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/Chemistry_Done' integer :: region #ifdef with_budgets #ifdef with_hdf4 type(THdfFile) :: io type(TSds) :: sds #endif real, dimension(:,:,:), allocatable :: collectb real,dimension(nbudg,nbud_vg,ntrace) :: buddry real,dimension(nbudg,nbud_vg,ntrace) :: budchem_result integer :: i, j, l, nzone, n, nsend, i1, i2, j1, j2 ! for buggy MPI (see comment in budget_global.F90 for details) real, dimension(nregions) :: sum_chemistry_all, sum_wetchem_all, sum_deposition_all real, dimension(nbudg,nbud_vg,nreac) :: budrrg_all real, dimension(nbudg,nbud_vg,nj) :: budrjg_all real, dimension(nbudg,nbud_vg,nreacw) :: budrwg_all real, dimension(nbudg,nbud_vg,nmark) :: budmarkg_all real, dimension(nbudg,nbud_vg,ntrace) :: buddry_all real, dimension(nbudg,nbud_vg,ntrace) :: budchem_all real, dimension(nbudg,nbud_vg,2) :: budeqsam_all real :: eminox_all #ifndef without_boundary real, dimension(nbudg,nbud_vg,ntracet) :: budstratg_all #endif #ifdef with_cariolle real, dimension(nbudg,nbud_vg,ncar) :: budcarg_all #endif #ifdef with_m7 real, dimension(nbudg,nbud_vg,ntracet) :: budm7_all real, dimension(nbudg,nbud_vg,nm7procs) :: budm7proc_all #endif #ifdef with_m7 ! Do not change this order. And if you do, look through the M7 code to pprocess (parameter process) where it is written with a fixed paramter. You can, however, change this order if it is made dynamic. Character(len=8), Dimension(nm7procs), Parameter :: m7procnames = (/'NucN ','NucSU ','Coag11N ','Coag12N ','Coag12SU',& 'Coag13N ','Coag13SU','Coag14N ','Coag14SU','Coag15N ','Coag15SU','Coag16N ','Coag16SU','Coag17N ','Coag17SU',& 'Coag22N ','Coag23N ','Coag23SU','Coag23BC','Coag23OC','Coag24N ','Coag24SU','Coag24BC','Coag24OC','Coag25N ',& 'Coag25BC','Coag25OC','Coag26N ','Coag26SU','Coag26BC','Coag26OC','Coag27N ','Coag27SU','Coag27BC','Coag27OC',& 'Coag33N ','Coag35N ','Coag35BC','Coag35OC','Coag55N ','Cond1SU ','Cond2SU ','Cond3SU ','Cond4SU ','Cond5SU ',& 'Cond6SU ','Cond7SU ','Age5N ','Age5BC ','Age5OC ','Age6N ','Age6DU ','Age7N ','Age7DU ','Grow1N ',& 'Grow1SU ','Grow2N ','Grow2SU ','Grow2BC ','Grow2OC ','Grow3N ','Grow3SU ','Grow3BC ','Grow3OC ','Grow3SS ',& 'Grow3DU ','Coa12SOA','Coa13SOA','Coa14SOA','Coa15SOA','Coa16SOA','Coa17SOA','Coa23SOA','Coa24SOA','Coa25SOA',& 'Coa26SOA','Coa27SOA','Coa35SOA','Cond1SOA','Cond2SOA','Cond3SOA','Cond4SOA','Cond5SOA','Age5SOA ','Grow1SOA',& 'Grow2SOA','Grow3SOA','NucSOA '/) ! Array for left hand sides, static attribute for the M7 process budget data set. Character(len=8), Dimension(nm7procs), Parameter :: m7proclhs = (/'void ','H2SO4 ','NUS_N ','NUS_N ','SO4NUS ',& 'NUS_N ','SO4NUS ','NUS_N ','SO4NUS ','NUS_N ','SO4NUS ','NUS_N ','SO4NUS ','NUS_N ','SO4NUS ',& 'AIS_N ','AIS_N ','SO4AIS ','BCAIS ','POMAIS ','AIS_N ','SO4AIS ','BCAIS ','POMAIS ','AII_N ',& 'BCAII ','POMAII ','AIS_N ','SO4AIS ','BCAIS ','POMAIS ','AIS_N ','SO4AIS ','BCAIS ','POMAIS ',& 'ACS_N ','AII_N ','BCAII ','POMAII ','AII_N ','H2SO4 ','H2SO4 ','H2SO4 ','H2SO4 ','H2SO4 ',& 'H2SO4 ','H2SO4 ','AII_N ','BCAII ','POMAII ','ACI_N ','DUACI ','COI_N ','DUCOI ','NUS_N ',& 'SO4NUS ','AIS_N ','SO4AIS ','BCAIS ','POMAIS ','ACS_N ','SO4ACS ','BCACS ','POMACS ','SSACS ',& 'DUACS ','SOANUS ','SOANUS ','SOANUS ','SOANUS ','SOANUS ','SOANUS ','SOAAIS ','SOAAIS ','SOAAII ',& 'SOAAIS ','SOAAIS ','SOAAII ','EL+SVOC ','EL+SVOC ','EL+SVOC ','EL+SVOC ','EL+SVOC ','SOAAII ','SOANUS ',& 'SOAAIS ','SOAACS ','ELVOC '/) ! Array for right hand sides, static attribute for the M7 process budget data set. Character(len=8), Dimension(nm7procs), Parameter :: m7procrhs = (/'NUS_N ','SO4NUS ','void ','void ','SO4AIS ',& 'void ','SO4ACS ','void ','SO4COS ','void ','SO4AIS ','void ','SO4ACS ','void ','SO4COS ',& 'void ','void ','SO4ACS ','BCACS ','POMACS ','void ','SO4COS ','BCCOS ','POMCOS ','void ',& 'BCAIS ','POMAIS ','void ','SO4ACS ','BCACS ','POMACS ','void ','SO4COS ','BCCOS ','POMCOS ',& 'void ','void ','BCACS ','POMACS ','void ','SO4NUS ','SO4AIS ','SO4ACS ','SO4COS ','SO4AIS ',& 'SO4ACS ','SO4COS ','AIS_N ','BCAIS ','POMAIS ','ACS_N ','DUACS ','COS_N ','DUCOS ','AIS_N ',& 'SO4AIS ','ACS_N ','SO4ACS ','BCACS ','POMACS ','COS_N ','SO4COS ','BCCOS ','POMCOS ','SSCOS ',& 'DUCOS ','SOAAIS ','SOAACS ','SOACOS ','SOAAIS ','SOAACS ','SOACOS ','SOAACS ','SOAACS ','SOAAIS ',& 'SOAACS ','SOACOS ','SOAACS ','SOANUS ','SOAAIS ','SOAACS ','SOACOS ','SOAAII ','SOAAIS ','SOAAIS ',& 'SOAACS ','SOACOS ','SOANUS '/) #endif #endif /* BUDGETS */ ! --- begin -------------------------------- #ifdef with_budgets eminox_all=0.0 ! Sum up contribution from each proc into root array call PAR_REDUCE_ELEMENT(budrrg, 'SUM', budrrg_all, status ) do region = 1, nregions CALL PAR_REDUCE(sum_chemistry(region), 'SUM', sum_chemistry_all(region), status) IF_NOTOK_RETURN(status=1) CALL PAR_REDUCE(sum_deposition(region), 'SUM', sum_deposition_all(region), status) IF_NOTOK_RETURN(status=1) CALL PAR_REDUCE(sum_wetchem(region), 'SUM', sum_wetchem_all(region), status) IF_NOTOK_RETURN(status=1) end do ! Sum up contribution from each proc into root array for "total LiNOx" call PAR_REDUCE( eminox, 'SUM', eminox_all, status ) eminox_all = eminox_all*1.e-9 ! kg -> Tg #ifdef with_hdf4 if ( isRoot ) then call Init(io, budget_file_global, 'write', status) IF_NOTOK_RETURN(status=1) call WriteAttribute(io, 'sum_chemistry' , sum_chemistry_all, status) IF_NOTOK_RETURN(status=1) call WriteAttribute(io, 'sum_wetchem' , sum_wetchem_all, status) IF_NOTOK_RETURN(status=1) call WriteAttribute(io, 'sum_deposition' , sum_deposition_all, status) IF_NOTOK_RETURN(status=1) call WriteAttribute(io, 'total_lightning_nox_emissions_Tg_N' , eminox_all, status) IF_NOTOK_RETURN(status=1) call Init(Sds, io, 'budrr', (/nbudg,nbud_vg,nreac/), 'real(8)', status) IF_NOTOK_RETURN(status=1) call WriteAttribute(Sds, 'ratnam', ratnam, 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, 'nreac', 'reaction number', (/(j, j=1,nreac)/), status) IF_NOTOK_RETURN(status=1) call WriteData(Sds, BUDRRG_all, status) IF_NOTOK_RETURN(status=1) call Done(Sds, status) IF_NOTOK_RETURN(status=1) end if #endif ! Sum up contribution from each proc into root array call PAR_REDUCE_ELEMENT(budrjg, 'SUM', budrjg_all, status ) #ifdef with_hdf4 if ( isRoot ) then call Init(Sds,io, 'budrj',(/nbudg,nbud_vg,nj/), 'real(8)', status) IF_NOTOK_RETURN(status=1) call WriteAttribute(Sds, 'jnam', jnam, 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, 'nj','reaction number', (/(j, j=1,nj)/), status) IF_NOTOK_RETURN(status=1) call WriteData(Sds, BUDRJG_all, status) IF_NOTOK_RETURN(status=1) call Done(Sds, status) IF_NOTOK_RETURN(status=1) end if #endif ! Lets us write down the WET CHEMISTRY call PAR_REDUCE_ELEMENT(budrwg, 'SUM', budrwg_all, status ) #ifdef with_hdf4 if ( isRoot ) then call Init(Sds,io, 'budrw',(/nbudg,nbud_vg,nreacw/), 'real(8)', status) IF_NOTOK_RETURN(status=1) call WriteAttribute(Sds, 'rwnam', rwnam, 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, 'nreacw','wet reaction number', (/(j, j=1,nreacw)/), status) IF_NOTOK_RETURN(status=1) call WriteData(Sds, BUDRWG_all, status) IF_NOTOK_RETURN(status=1) call Done(Sds, status) IF_NOTOK_RETURN(status=1) end if #endif ! Sum up contribution from each proc into root array call PAR_REDUCE_ELEMENT(budmarkg, 'SUM', budmarkg_all, status ) #ifdef with_hdf4 if ( isRoot ) then call Init(Sds,io, 'budmark',(/nbudg,nbud_vg,nmark/), 'real(8)', status) IF_NOTOK_RETURN(status=1) call WriteAttribute(Sds, 'marknam', marknam, 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, 'nmark','marked tracer number', (/(j, j=1,nmark)/), status) IF_NOTOK_RETURN(status=1) call WriteData(Sds, BUDMARKG_all, status) IF_NOTOK_RETURN(status=1) call Done(Sds, status) IF_NOTOK_RETURN(status=1) end if #endif ! DRYDEP budget buddry = 0.0 do region=1,nregions ! aggregates horizontally budget call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) do n=1,ntrace do j=j1,j2 do i=i1,i2 nzone = budg_dat(region)%nzong(i,j) buddry(nzone,1,n) = buddry(nzone,1,n) + buddep_dat(region)%dry(i,j,n) end do end do end do call PAR_REDUCE_ELEMENT(buddry, 'SUM', buddry_all, status ) ! contribution from all procs ! gather and write Non-Horizontally-Aggregated-Budgets if (NHAB) then if ( isRoot ) then allocate( collectb(im(region),jm(region),ntrace )) else allocate( collectb(1,1,1)) end if call GATHER( dgrid(region), buddep_dat(region)%dry, collectb, 0, status) #ifdef with_hdf4 if(isRoot) then call Init(Sds,io, 'buddep_dat_dry_'//region_name(region),(/im(region),jm(region),ntrace/), 'real(4)', status) 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, 'ntrace','tracer index', (/(j, j=1,ntrace)/), status) IF_NOTOK_RETURN(status=1) call WriteData(Sds, collectb, status) IF_NOTOK_RETURN(status=1) call Done(Sds, status) IF_NOTOK_RETURN(status=1) endif #endif deallocate(collectb) endif ! NHAB enddo ! regions !-- write more aggregated budgets call PAR_REDUCE_ELEMENT(budchem, 'SUM', budchem_all, status ) call PAR_REDUCE_ELEMENT(budeqsam, 'SUM', budeqsam_all, status ) #ifndef without_boundary call PAR_REDUCE_ELEMENT(budstratg,'SUM', budstratg_all, status ) #endif #ifdef with_m7 call PAR_REDUCE_ELEMENT(budm7, 'SUM', budm7_all, status ) call PAR_REDUCE_ELEMENT(budm7proc,'SUM', budm7proc_all, status ) #endif if(isRoot) then #ifdef with_hdf4 call Init(Sds,io, 'buddry',(/nbudg,nbud_vg,ntrace/), '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, 'ntrace','tracer index', (/(j, j=1,ntrace)/), status) IF_NOTOK_RETURN(status=1) call WriteData(Sds,buddry_all,status) IF_NOTOK_RETURN(status=1) call Done(Sds, status) IF_NOTOK_RETURN(status=1) Call Init(Sds,io,"budchem",(/nbudg,nbud_vg,ntrace/),'real(8)',status) call SetDim(Sds, 0, 'nbudg','horizontal budget', (/(j, j=1,nbudg)/), status) call SetDim(Sds, 1, 'nbud_vg','vertical budget', (/(j, j=1,nbud_vg)/), status) call SetDim(Sds, 2, 'ntrace','tracer number', (/(j, j=1,ntrace)/), status) IF_NOTOK_RETURN(status=1) call WriteAttribute(Sds, 'comment', 'EBI budget stuffed together. Hopefully it does not explode', status) #endif ! correct chemistry budget budchem_result = budchem_all+buddry_all #ifndef without_emission budchem_result(:,:,inox) = budchem_result(:,:,inox) - budemig_all(:,:,inox) #endif #ifdef with_hdf4 call WriteData(Sds,budchem_result,status) IF_NOTOK_RETURN(status=1) call Done(Sds, status) IF_NOTOK_RETURN(status=1) Call Init(Sds,io,"budeqsam",(/nbudg,nbud_vg,2/),'real(8)',status) call SetDim(Sds, 0, 'nbudg','horizontal budget', (/(j, j=1,nbudg)/), status) call SetDim(Sds, 1, 'nbud_vg','vertical budget', (/(j, j=1,nbud_vg)/), status) call SetDim(Sds, 2, 'AB_reac','Acid-base reaction number', (/(j, j=1,2)/), status) IF_NOTOK_RETURN(status=1) call WriteAttribute(Sds, 'comment', 'EQSAM evil acid base reaction. Positive forms aerosol material.', status) call WriteData(Sds,budeqsam_all,status) IF_NOTOK_RETURN(status=1) call Done(Sds, status) IF_NOTOK_RETURN(status=1) #ifndef without_boundary Call Init(Sds,io,"budstrat",(/nbudg,nbud_vg,ntracet/),'real(8)',status) call SetDim(Sds, 0, 'nbudg','horizontal budget', (/(j, j=1,nbudg)/), status) call SetDim(Sds, 1, 'nbud_vg','vertical budget', (/(j, j=1,nbud_vg)/), status) call SetDim(Sds, 2, 'ntracet','tracer number', (/(j, j=1,ntracet)/), status) IF_NOTOK_RETURN(status=1) call WriteAttribute(Sds, 'comment', 'Stratosphere (from boundary.F90)', status) call WriteData(Sds,budstratg_all,status) IF_NOTOK_RETURN(status=1) call Done(Sds, status) IF_NOTOK_RETURN(status=1) #endif #endif endif #ifdef with_cariolle call PAR_REDUCE_ELEMENT(budcarg, 'SUM', budcarg_all, status ) #ifdef with_hdf4 if ( isRoot ) then call Init(Sds,io, 'budcar',(/nbudg,nbud_vg,ncar/), 'real(8)', status) IF_NOTOK_RETURN(status=1) call WriteAttribute(Sds, 'carnam', carnam, 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, 'ncar','Cariolle tracer number', (/(j, j=1,ncar)/), status) IF_NOTOK_RETURN(status=1) call WriteData(Sds,budcarg_all,status) IF_NOTOK_RETURN(status=1) call Done(Sds, status) IF_NOTOK_RETURN(status=1) end if #endif #endif #ifdef with_m7 #ifdef with_hdf4 if ( isRoot ) then ! M7 process budget Call Init(Sds,io,"budm7proc",(/nbudg,nbud_vg,nm7procs/),'real(8)',status) call SetDim(Sds, 0, 'nbudg','horizontal budget', (/(j, j=1,nbudg)/), status) call SetDim(Sds, 1, 'nbud_vg','vertical budget', (/(j, j=1,nbud_vg)/), status) call SetDim(Sds, 2, 'nm7proc','M7 process number', (/(j, j=1,nm7procs)/), status) IF_NOTOK_RETURN(status=1) call WriteAttribute(Sds, 'procnames',m7procnames, status) call WriteAttribute(Sds, 'lefthandsides',m7proclhs, status) call WriteAttribute(Sds, 'righthandsides',m7procrhs, status) call WriteData(Sds,budm7proc_all,status) IF_NOTOK_RETURN(status=1) call Done(Sds, status) IF_NOTOK_RETURN(status=1) ! M7 budget Call Init(Sds,io,"budm7",(/nbudg,nbud_vg,ntracet/),'real(8)',status) call SetDim(Sds, 0, 'nbudg','horizontal budget', (/(j, j=1,nbudg)/), status) call SetDim(Sds, 1, 'nbud_vg','vertical budget', (/(j, j=1,nbud_vg)/), status) call SetDim(Sds, 2, 'ntracet','tracer number', (/(j, j=1,ntracet)/), status) IF_NOTOK_RETURN(status=1) call WriteAttribute(Sds, 'comment', 'M7 budget stuffed together. Hopefully the M7 tracers have now correct budgets', status) call WriteData(Sds,budm7_all,status) IF_NOTOK_RETURN(status=1) call Done(Sds, status) IF_NOTOK_RETURN(status=1) endif if ( isRoot ) then call Done(io, status) IF_NOTOK_RETURN(status=1) endif #endif #endif do region = 1, nregions #ifdef with_m7 deallocate(diag_prod (region)%prod) deallocate(AC_diag_prod(region)%prod) deallocate(AC_O3_lp (region)%prod) deallocate(AC_loss (region)%prod) #endif deallocate(buddep_dat(region)%dry) deallocate(o3p(region)%d3) deallocate(o3l(region)%d3) enddo deallocate(budchem) deallocate(budeqsam) #ifdef with_m7 deallocate(budm7proc) deallocate(budm7) deallocate(d_nucle) deallocate(m_nucle_su) deallocate(m_nucle_soa) deallocate(growth1_2) deallocate(coag_sink_nuc) deallocate(cond1_soa) deallocate(cond2_soa) deallocate(cond3_soa) deallocate(cond4_soa) deallocate(cond5_soa) deallocate(cond1_su) #endif #endif /* BUDGETS */ #ifdef with_m7 call free_m7_data #endif ! ok status = 0 END SUBROUTINE CHEMISTRY_DONE !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: CHEMIE ! ! !DESCRIPTION: performs chemistry transformations of the tracers !\\ !\\ ! !INTERFACE: ! SUBROUTINE CHEMIE( region, tr, status ) ! ! !USES: ! use GO, only : TDate use binas, only : Avog, pi, xmair, Rgas_air use dims, only : itaur, nchem, tref, sec_month, revert, okdebug use dims, only : dx, xref, dy, yref, im, jm, lm use dims, only : nregions, ndyn, ndyn_max use dims, only : xbeg, ybeg, adv_scheme, mlen use chem_param, only : xmbc, xmpom, xmnacl, xmdust !EV use chem_param, only : maxtrace, n_extra, iso4, ino3_a, inh4, ih2opart use chem_param, only : irinc, ino, xmn, inox, ieno, iairm, iairn, ntrace use chem_param, only : ntracet, fscale, imsa, iacid, iair, ih2o, io2, ih2on use chem_param, only : iairm,i_temp, i_pres, inh3, ihno3, xmh2o, names, iterp, xmelvoc, xmsvoc #ifdef with_m7 use chem_param, only : iso4nus, iso4ais, iso4acs, iso4cos, ibcais, ibcacs, ibccos, ibcaii, isoanus, ipomais, ipomacs use chem_param, only : ipomcos, ipomaii, issacs, isscos, iduacs, iducos, iduaci,iducoi use chem_param, only : isoaais, isoaacs, isoacos, isoaaii use chem_param, only : inus_n, xmnumb, iais_n, iacs_n, icos_n, iaii_n, iaci_n, icoi_n use chem_param, only : nh4no3_factor, nh4no3_density, msa_density use mo_aero_m7, only : cmr2ram, ram2cmr use binas, only : rol use chem_param, only : ielvoc, isvoc #endif use chem_param, only : idz use meteodata, only : m_dat, gph_dat use reaction_data, only : nreac, aerdens use chem_rates, only : calrates use eqsam_param, only : eqsam_v03d use toolbox, only : dumpfield use datetime, only : tau2date, get_day use global_data, only : region_dat, mass_dat, chem_dat #ifndef without_photolysis #ifdef with_optics use photolysis, only : aerosol_info_m7 #else use photolysis, only : aerosol_info #endif use photolysis, only : daysim, get_sza, photo_flux use photolysis_data #endif use ebischeme, only : ebi #ifdef with_budgets #ifdef with_m7 use ebischeme, only :diag_prod,nprod,AC_diag_prod,iprod_soa3dmon,iprod_soa2dhour #endif #endif #ifndef without_emission use dims, only : at, bt use emission_nox, only : eminox_lightning, nox_emis_3d use emission_nox, only : nox_emis_3d_bb_app, emission_nox_bb_daily_cycle #endif #ifdef with_budgets use emission_data, only : budemi_dat, sum_emission use budget_global, only : budg_dat, nzon_vg #endif #ifndef without_sedimentation use sedimentation, only : rh #endif #ifdef with_cariolle use chem_param, only : io3 use chemistry_cariolle, only : l_cariolle, with_trop_force #endif ! use tracer_data, only : tracer_print ! for debugging #ifdef with_m7 use mo_aero, only : nsoa !RM use mo_aero_m7, only : naermod, nmod, nsol use m7_data, only : rw_mode, rwd_mode, dens_mode, h2o_mode !use emission_data, only : oc2pom !factor for conversion of OC mass to POM #ifdef with_budgets Use m7_data, only : nm7procs #endif !#ifdef with_optics ! Use Optics, only : nwl, aod_count, AOD, calculate_aop ! Use chem_param, only : xmso4, xmno3 !#endif external :: m7 #endif /* M7 */ ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region type(TDate), intent(in) :: tr(2) ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 22 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/chemie' ! --- local ------------------------------ real,dimension(:,:,:),pointer :: m real,dimension(:,:,:,:),pointer :: rm #ifdef slopes real,dimension(:,:,:,:),pointer :: rxm, rym, rzm #ifdef secmom real,dimension(:,:,:,:),pointer :: rxxm, rxym, rxzm, ryym, ryzm, rzzm #endif #endif real,dimension(:,:,:,:),pointer :: rmc real,dimension(:,:,:),pointer :: r_cloud ! store per level actinic fluxes, photolysis rates, and reaction rates #ifndef without_photolysis real,dimension(:,:,:),allocatable :: fact #endif real,dimension(:,:,:),allocatable :: rj, rr ! new photolysis array real,dimension(:,:,:,:),allocatable :: rj_new ! store concentrations for chemistry .... real,dimension(:,:,:),allocatable :: y, y0, yhelp ! store information that is transferred ! between chemistry modules real,dimension(:,:,:),allocatable :: ye real,dimension(:,:), allocatable :: th, tha, sza, mu, reff real :: starttime,deltat ! eqsam: real,dimension(:,:), allocatable :: yi, yo ! sad: real,dimension(:,:,:), allocatable :: sad_aer_out, sad_cld_out, sad_ice_out #ifdef with_m7 ! m7: real, dimension(:), pointer :: dxyp real, dimension(:,:), allocatable :: presm7, rhm7, tempm7, dryairm7 ! meteo parameters real, dimension(:,:), allocatable :: h2so4m7 ! h2so4 gas phase real, dimension(:,:,:), allocatable :: aemassm7 ! aerosol mass real, dimension(:,:,:), allocatable :: aenumm7 ! aerosol number real, dimension(:,:), allocatable :: elvocm7, svocm7 ! ELVOC & SVOC gas phase !RM !#ifdef with_optics ! Integer :: iband ! real,dimension(:,:,:),allocatable :: Ext ! extinction coefficient [1/m] ! real,dimension(:,:,:),allocatable :: a ! single-scattering albedo [unitless] ! real,dimension(:,:,:),allocatable :: g ! asymmetry parameter [unitless] !#endif real, dimension(:,:,:), allocatable :: aedensm7 ! aerosol density real, dimension(:,:,:), allocatable :: aewatm7 ! aerosol water content real, dimension(:,:,:), allocatable :: aeradm7 ! aerosol radius in equilibrium with water real, dimension(:,:,:), allocatable :: aeradrm7 ! aerosol dry radius ! ! To convert units from M7 to TM5 units real, parameter :: convnumb = 1.e3/xmnumb*Avog ! This value is needed directly real, parameter :: convsu = 1. ! Nothing happens for sulfate. real, parameter :: convbc = 1e-12/xmbc*Avog ! BC !>>> TvN ! I have the impression that OC in M7 actually refers to POM. ! doc seems to be the density of POM not OC, and is equal to the density of BC (dbc). ! Compare with GLOMAP, where the densities of POM and BC are the same (Mann et al., GMD, 2010). ! This would imply that also the mass of POM and not OC should be used in M7, ! and that "OC" in the paper by Vignati et al. (JGR, 2004) is actually referring to POM, ! consistent with the terminology used by Stier et al. (ACP, 2005). !real, parameter :: convoc = 1e-12/xmpom*Avog*oc2pom ! OC (M7) POM (TM5) real, parameter :: convoc = 1e-12/xmpom*Avog ! POM (in both M7 and TM5) !<<< TvN !RM real, parameter :: convelvoc = 1e-12/xmelvoc*Avog ! ELVOC real, parameter :: convsvoc = 1e-12/xmsvoc*Avog ! S/LVOC real, parameter :: convss = 1e-12/xmnacl*Avog ! SS real, parameter :: convdu = 1e-12/xmdust*Avog ! DU #ifdef with_budgets real,dimension(:,:,:), allocatable :: processm7 ! M7 processes in M7 units. ! Array that convertts M7 processes from M7 to TM5 units. real, dimension(nm7procs), parameter :: processm7totm5 = (/ & ! 1 2 3 4 5 6 7 8 9 10 convnumb,convsu ,convnumb,convnumb,convsu ,convnumb,convsu ,convnumb,convsu,convnumb,& ! 0? convsu ,convnumb,convsu ,convnumb,convsu ,convnumb,convnumb,convsu ,convbc,convoc ,& ! 1? convnumb,convsu ,convbc ,convoc ,convnumb,convbc ,convoc ,convnumb,convsu,convbc ,& ! 2? convoc ,convnumb,convsu ,convbc ,convoc ,convnumb,convnumb,convbc ,convoc,convnumb,& ! 3? convsu ,convsu ,convsu ,convsu ,convsu ,convsu ,convsu ,convnumb,convbc,convoc ,& ! 4? convnumb,convdu ,convnumb,convdu ,convnumb,convsu ,convnumb,convsu ,convbc,convoc ,& ! 5? convnumb,convsu ,convbc ,convoc ,convss ,convdu ,convoc ,convoc ,convoc,convoc ,& ! 6? convoc ,convoc ,convoc ,convoc ,convoc ,convoc ,convoc ,convoc ,convoc,convoc ,& ! 7? convoc ,convoc ,convoc ,convoc ,convoc ,convoc ,convoc ,convoc/) ! 8? #endif real :: dryvol_m7, vol_nh4no3, vol_msa, vol_h2o real :: number_conc integer :: imode #endif /* M7 */ real :: dtime, gmt, dxx, dyy, lat, lon real :: px, skf, x, month2dt integer :: day, level, nzone, nzone_v integer,dimension(6) :: idater, idate_temp integer :: i, j, l, n integer :: ispec integer :: imax, nca, nco, iloop integer :: eqsam_option real :: dry_aerosol_mass, dry_aerosol_volume real :: water_aerosol_mass, water_aerosol_volume, ccs, s integer :: i1, i2, j1, j2, lmr, is3, ie3, klm ! --- begin --------------------------------- if ( okdebug ) write (*,*) 'start of chemistry' m => m_dat (region)%data rmc => chem_dat(region)%rm rm => mass_dat(region)%rm #ifdef slopes rxm => mass_dat(region)%rxm rym => mass_dat(region)%rym rzm => mass_dat(region)%rzm #ifdef secmom rxxm => mass_dat(region)%rxxm rxym => mass_dat(region)%rxym rxzm => mass_dat(region)%rxzm ryym => mass_dat(region)%ryym ryzm => mass_dat(region)%ryzm rzzm => mass_dat(region)%rzzm #endif #endif #ifdef with_m7 dxyp => region_dat(region)%dxyp #endif r_cloud => phot_dat(region)%cloud_reff ! ------- Time stuff call tau2date(itaur(region), idater) ! current time dtime = nchem/(2*tref(region)) ! time step month2dt = dtime/sec_month ! conversion to emission per timestep call tau2date(itaur(region)+revert*nint(dtime*0.5),idate_temp) ! get date @ halfway interval. gmt = idate_temp(4) + idate_temp(5)/60.0 + idate_temp(6)/3600.0 ! GMT in hours day = get_day(idate_temp(2),idate_temp(3),mlen) if ( okdebug ) then write(gol,*)'chemistry: region ',region, ' at date: ',idater,' with time step ',dtime call goPr write(gol,*)'chemistry: day, gmt, date@halfway interval ', day, gmt, idate_temp call goPr end if ! ------- Resolution and bounds dxx = dx/xref(region) ! gridsize at current region dyy = dy/yref(region) ! gridsize at current region call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) lmr = lm(region) is3=(i1-1)*3+1 ! to cover three times oversampled ie3=i2*3 ! allocate(tha (is3:ie3,j1:j2)) allocate(tha (im(region)*3,j1:j2)) allocate(sza (i1:i2, j1:j2 )) allocate(mu (i1:i2, j1:j2 )) allocate(reff (i1:i2, j1:j2 )) allocate(fact (i1:i2, j1:j2, nbands_trop )) allocate(rj (i1:i2, j1:j2, nj )) allocate(rr (i1:i2, j1:j2, nreac )) allocate(rj_new(i1:i2, j1:j2, lmr, nj )) allocate(y (i1:i2, j1:j2, maxtrace )) allocate(y0 (i1:i2, j1:j2, maxtrace )) allocate(yhelp (i1:i2, j1:j2, ntrace )) allocate(ye (i1:i2, j1:j2, n_extra )) ! Eqsam allocate(sad_cld_out(i1:i2, j1:j2, lmr )) ; sad_cld_out = 0. allocate(sad_ice_out(i1:i2, j1:j2, lmr )) ; sad_ice_out = 0. allocate(sad_aer_out(i1:i2, j1:j2, lmr )) ; sad_aer_out = 0. imax = (i2-i1+1)*(j2-j1+1) !nb of grid box per layer (was im(region)*jm(region)) nca = 11 nco = 37 ! was 36, JadB added one for the water associated with just the nitrate. allocate(yi(imax,nca)) ; yi = 0.0 allocate(yo(imax,nco)) ; yo = 0.0 #ifdef with_m7 allocate( presm7 (imax, 1 )) allocate( rhm7 (imax, 1 )) allocate( tempm7 (imax, 1 )) allocate( dryairm7(imax, 1 )) allocate( h2so4m7 (imax, 1 )) allocate( aemassm7(imax, 1, naermod)) allocate( aenumm7 (imax, 1, nmod )) allocate( aedensm7(imax, 1, nmod )) allocate( aewatm7 (imax, 1, nmod )) allocate( aeradm7 (imax, 1, nmod )) allocate( aeradrm7(imax, 1, nsol )) allocate( elvocm7 (imax, 1 )) allocate( svocm7 (imax, 1 )) presm7 =0. rhm7 =0. tempm7 =0. dryairm7=0. h2so4m7 =0. aemassm7=0. aenumm7 =0. aedensm7=0. aewatm7 =0. aeradm7 =0. aeradrm7=0. elvocm7 =0. svocm7 =0. #ifdef with_budgets allocate(processm7(imax,1,nm7procs)) #endif #endif tha = 0.0 sza = 0.0 mu = 0.0 fact = 0.0 rj = 0.0 rr = 0.0 rj_new = 0.0 y = 0.0 y0 = 0.0 ye = 0.0 ! alternative calculation of zenith angle with higher sampling (time and space). ! yields factor of two lower angles in the tropics! ! Difference is resolution dependent. The coarser the larger the differences starttime = real(idate_temp(4)) + real(idate_temp(5))/60.0 + real(idate_temp(6))/3600. deltat = nchem/3600.0 call daysim(region, day, is3, i1,i2,j1, j2, tha) call get_sza(region, i1, j1, i2, j2, starttime, deltat, tha, sza) ! ascribe mu do j=j1,j2 do i=i1,i2 mu(i,j)=max(1.e-5,COS(sza(i,j)*pi/180.)) enddo enddo #ifdef with_optics call aerosol_info_m7(region, status) IF_NOTOK_RETURN(status=1) #else call aerosol_info(region) #endif ! calculate actinic flux outside level loop for online fluxes call photo_flux(region, i1, j1, sza, rj_new) ! scale NOx BMB if needed #ifndef without_emission call emission_nox_bb_daily_cycle(status) IF_NOTOK_RETURN(status=1) #endif ! compute global LiNox budget #ifdef with_budgets #ifndef without_emission #ifndef without_convection eminox=eminox+sum(eminox_lightning(region)%d3)*month2dt ! kg [N] month-1 * dt /(sec_month) = kg [N] #endif #endif #endif LEVS: do level = 1, lmr rj(:,:,:) = rj_new (:,:,level,:) reff(:,:) = r_cloud(:,:,level) #ifdef with_budgets nzone_v=nzon_vg(level) #endif ! get/calculate all the fields needed in chemistry call GET_EXTRA( region, level, ye, i1, j1) ! calculate the rinc: relative radius growth of aerosols due to water ! Eqsam wants to calculate aerosol radii ! irinc is used by chem_rates in calchet2 ! calchet2 is used by no one ! So, this part is not so important ! ----------------------------------- ! To get this code back, take a stone-age copy of chemistry.F90 ! But, to prevent even more errors. ! ----------------------------------- ye(:,:,irinc) = 1.0 call CALRATES(region, level, i1, j1, rj, rr, reff, ye, dtime, & sad_cld_out(:,:,level), sad_ice_out(:,:,level), sad_aer_out(:,:,level), status) ! also returns rh in ye IF_ERROR_RETURN(status=1) ! CALCULATE NOX EMISSIONS ! --------------------------- do j=j1,j2 do i=i1,i2 ! nox emis in this cell; init to zero: x = 0.0 ! kg(N)/month #ifndef without_emission ! add total 3d emissions if any if(associated(nox_emis_3d(region)%d3)) then x = x + nox_emis_3d(region)%d3(i,j,level) end if if(associated(nox_emis_3d_bb_app(region)%d3)) then x = x + nox_emis_3d_bb_app(region)%d3(i,j,level) end if #ifndef without_convection ! add lightning emissions: x = x + eminox_lightning(region)%d3(i,j,level) #endif #endif /* EMISS */ ! emission budget (since nox emissions are added in chemistry) #ifdef with_budgets #ifndef without_emission budemi_dat(region)%emi(i,j,nzone_v,inox) = & budemi_dat(region)%emi(i,j,nzone_v,inox) + x/xmn*month2dt*1e3 ! g if (inox == 1) sum_emission(region) = sum_emission(region) + x*month2dt ! kg N #endif #endif ! put emissions in the units #/s/cm3 used in chemistry ye(i,j,ieno) = x/ye(i,j,iairm)*ye(i,j,iairn)* & (xmair/xmn)/sec_month !kg N/month ----> #(NO)/cm3/s end do end do ! INITIAL CONCENTRATIONS ! --------------------------- do n=1,ntracet do j=j1,j2 do i=i1,i2 y(i,j,n) = rm(i,j,level,n)/ye(i,j,iairm)*ye(i,j,iairn)* fscale(n) !kg ----> #/cm3 end do end do end do do n=ntracet+1,ntrace do j=j1,j2 do i=i1,i2 ! remember nontransported tracers are in rmc y(i,j,n) =rmc(i,j,level,n)/ye(i,j,iairm)*ye(i,j,iairn)* fscale(n) !kg ----> #/cm3 end do end do end do ! FILL EXTRA SPECIES : (from ntrace+1 to maxtrace) ! --------------------------- do j=j1,j2 do i=i1,i2 ! these are in #/cm3 and used in budget....should be in y.... y(i,j,iair) = ye(i,j,iairn) y(i,j,ih2o) = ye(i,j,ih2on) ! introduced to calculate correct budget for JO2 y(i,j,io2) = 1. #ifdef with_m7 !RM y(i,j,ielvoc)=0. y(i,j,isvoc)=0. #endif end do end do ! BACKUP INIT CONCENTRATIONS ! --------------------------- y0 = y do klm=1,ntrace yhelp(:,:,klm) = y(:,:,klm)*ye(:,:,iairm)*1000./xmair/y(:,:,iair) end do ! EBI ! --------------------------- CALL EBI( region, level, i1, j1, rj, rr, y, ye, status) IF_ERROR_RETURN(status=1) ! EBI scheme budgets #ifdef with_budgets do j=j1,j2 do i=i1,i2 nzone=budg_dat(region)%nzong(i,j) ! global budget zone budchem(nzone,nzone_v,1:ntrace) = budchem(nzone,nzone_v,1:ntrace) + & y(i,j,1:ntrace)*ye(i,j,iairm)*1000./xmair/y(i,j,iair)-yhelp(i,j,:) !production data !d_gas_prod_so4=1 !d_liq_prod_so4=2 !d_prod_elvoc=3 !d_prod_svoc=4 ! cannot be done here, or would have also the existing so4-> ebischeme !diag_prod(region)%prod(i,j,1)=y(i,j,iso4) !diag_prod(region)%prod(i,j,1)=y(i,j,ielvoc) !=y(i,j,n)/y(i,j,iair)*ye(i,j,iairm)/fscale(n) ! y in #cm-3 ! y -> #cm-3s-1 : /y(air)-> #/#(air) : *y(airm)/xmair*xmelvoc -> kg(elvoc) #ifdef with_m7 diag_prod(region)%prod(i,j,level,3)=diag_prod(region)%prod(i,j,level,3)+y(i,j,ielvoc)/y(i,j,iair)*ye(i,j,iairm)/xmair*xmelvoc diag_prod(region)%prod(i,j,level,4)=diag_prod(region)%prod(i,j,level,4)+y(i,j,isvoc)/y(i,j,iair)*ye(i,j,iairm)/xmair*xmsvoc !AERCHEMMIP AC_diag_prod(region)%prod(i,j,level,iprod_soa3dmon)=AC_diag_prod(region)%prod(i,j,level,iprod_soa3dmon)+y(i,j,ielvoc)/y(i,j,iair)*ye(i,j,iairm)/xmair*xmelvoc+y(i,j,isvoc)/y(i,j,iair)*ye(i,j,iairm)/xmair*xmsvoc !instantaneous 2d hour/6hour, make it hourly mean AC_diag_prod(region)%prod(i,j,level,iprod_soa2dhour)=AC_diag_prod(region)%prod(i,j,level,iprod_soa2dhour)+y(i,j,ielvoc)/y(i,j,iair)*ye(i,j,iairm)/xmair*xmelvoc+y(i,j,isvoc)/y(i,j,iair)*ye(i,j,iairm)/xmair*xmsvoc #endif end do end do #endif !--- Begin of M7: ------------------------------------------------------ #ifdef with_m7 ! M7 iloop = 0 do j=j1,j2 do i=i1,i2 #ifdef with_budgets ! For the M7 budgets, we first subtract the old concentrations. nzone=budg_dat(region)%nzong(i,j) ! global budget zone budm7(nzone,nzone_v,:)=budm7(nzone,nzone_v,:)-y(i,j,1:ntracet)*ye(i,j,iairm)*1000./xmair/y(i,j,iair) #endif iloop = iloop + 1 ! meteo parameters presm7 (iloop,1) = ye(i,j,i_pres) ! Pa rhm7 (iloop,1) = rh(region)%d3(i,j,level) ! fraction, in cloud free part! tempm7 (iloop,1) = ye(i,j,i_temp) ! K dryairm7(iloop,1) = presm7(iloop,1)/(Rgas_air*tempm7(iloop,1)) ! density of dry air ! H2SO4 gas concentrations molecule/cm3 h2so4m7(iloop,1) = y(i,j,iso4)/convsu ! molecule/cm3 ----> molecules/cm3 if (nsoa==2) then !RM ELVOC+SVOC elvocm7(iloop,1) = y(i,j,ielvoc)/convelvoc svocm7(iloop,1) = y(i,j,isvoc)/convsvoc end if ! sulphate mass in the 4 soluble modes molecule/cm3 aemassm7(iloop,1,1) = y(i,j,iso4nus)/convsu aemassm7(iloop,1,2) = y(i,j,iso4ais)/convsu aemassm7(iloop,1,3) = y(i,j,iso4acs)/convsu aemassm7(iloop,1,4) = y(i,j,iso4cos)/convsu ! molecule/cm3 ----> molecules/cm3 ! BC mass in the 3 soluble modes and 1 insoluble aemassm7(iloop,1,5) = y(i,j,ibcais)/convbc aemassm7(iloop,1,6) = y(i,j,ibcacs)/convbc aemassm7(iloop,1,7) = y(i,j,ibccos)/convbc aemassm7(iloop,1,8) = y(i,j,ibcaii)/convbc ! molecule/cm3 -----> ug/m3 ! was: POM (tm5) OC (m7) mass in the 3 soluble modes and 1 insoluble ! POM mass in the 3 soluble modes and 1 insoluble aemassm7(iloop,1,9) = y(i,j,ipomais)/convoc aemassm7(iloop,1,10) = y(i,j,ipomacs)/convoc aemassm7(iloop,1,11) = y(i,j,ipomcos)/convoc aemassm7(iloop,1,12) = y(i,j,ipomaii)/convoc ! molecule/cm3 -----> ug/m3 ! SOA mass in the 3 soluble modes and 1 insoluble aemassm7(iloop,1,19) = y(i,j,isoanus)/convoc aemassm7(iloop,1,20) = y(i,j,isoaais)/convoc aemassm7(iloop,1,21) = y(i,j,isoaacs)/convoc aemassm7(iloop,1,22) = y(i,j,isoacos)/convoc aemassm7(iloop,1,23) = y(i,j,isoaaii)/convoc ! molecule/cm3 -----> ug/m3 ! SS mass in the 2 soluble modes aemassm7(iloop,1,13) = y(i,j,issacs)/convss aemassm7(iloop,1,14) = y(i,j,isscos)/convss ! molecule/cm3 -----> ug/m3 ! DUST mass in the 2 soluble modes and 2 insoluble aemassm7(iloop,1,15) = y(i,j,iduacs)/convdu aemassm7(iloop,1,16) = y(i,j,iducos)/convdu aemassm7(iloop,1,17) = y(i,j,iduaci)/convdu aemassm7(iloop,1,18) = y(i,j,iducoi)/convdu ! molecule/cm3 -----> ug/m3 ! number concentrations of the 7 modes aenumm7(iloop,1,1) = y(i,j,inus_n)/convnumb aenumm7(iloop,1,2) = y(i,j,iais_n)/convnumb aenumm7(iloop,1,3) = y(i,j,iacs_n)/convnumb aenumm7(iloop,1,4) = y(i,j,icos_n)/convnumb aenumm7(iloop,1,5) = y(i,j,iaii_n)/convnumb aenumm7(iloop,1,6) = y(i,j,iaci_n)/convnumb aenumm7(iloop,1,7) = y(i,j,icoi_n)/convnumb ! ! molecule/cm3 -----> particle/cm3: ! 1 # air !rm in total aerosol number (#aer) ----> #aer * ------ * ------ is y ! kg(air) cm-3 !#aer # air 1 !---- * ------- ----> 1 kg(air) = 10^3 g air = 10^3 * ----- * Avog !cm-3 kg(air) xmair ! end do end do ! The M7 processes will be calculated when budgets are used. The extra ! argument will appear when the compiler flag with_budgets is set. #ifdef with_budgets call M7(iloop, imax, 1, & ! TM5 indices presm7, rhm7, tempm7, & ! " thermodynamics h2so4m7, elvocm7, svocm7, aemassm7, aenumm7, & ! M7 tracers aedensm7, aewatm7, aeradm7, aeradrm7, dtime, & ! M7 aerosol propertis processm7) ! The extra argument, in which the rates of each process is stored in M7 units. #else call M7(iloop, imax, 1, & ! TM5 indices presm7, rhm7, tempm7, & ! " thermodynamics h2so4m7, elvocm7, svocm7, aemassm7, aenumm7, & ! M7 tracers aedensm7, aewatm7, aeradm7, aeradrm7, dtime) ! M7 aerosol propertis #endif iloop = 0 do j=j1,j2 do i=i1,i2 iloop = iloop + 1 ! H2SO4 gas concentrations in molecule/cm3 y(i,j,iso4) = h2so4m7(iloop,1)*convsu ! molecule/cm3 ----> molecule/cm3 if (nsoa==2)then !RM ELVOC+SVOC y(i,j,ielvoc) = elvocm7(iloop,1)*convelvoc y(i,j,isvoc) = svocm7(iloop,1)*convsvoc end if ! sulphate mass in the 4 soluble modes y(i,j,iso4nus) = aemassm7(iloop,1,1)*convsu y(i,j,iso4ais) = aemassm7(iloop,1,2)*convsu y(i,j,iso4acs) = aemassm7(iloop,1,3)*convsu y(i,j,iso4cos) = aemassm7(iloop,1,4)*convsu ! molecule/cm3 ----> molecule/cm3 ! BC mass in the 3 soluble modes and 1 insoluble y(i,j,ibcais) = aemassm7(iloop,1,5)*convbc y(i,j,ibcacs) = aemassm7(iloop,1,6)*convbc y(i,j,ibccos) = aemassm7(iloop,1,7)*convbc y(i,j,ibcaii) = aemassm7(iloop,1,8)*convbc ! ug/m3 ----> molecule/cm3 ! was: POM (tm5) OC (m7) mass in the 3 soluble modes and 1 insoluble ! POM mass in the 3 soluble modes and 1 insoluble y(i,j,ipomais) = aemassm7(iloop,1,9)*convoc y(i,j,ipomacs) = aemassm7(iloop,1,10)*convoc y(i,j,ipomcos) = aemassm7(iloop,1,11)*convoc y(i,j,ipomaii) = aemassm7(iloop,1,12)*convoc ! ug/m3 ----> molecule/cm3 ! POM mass in the 3 soluble modes and 1 insoluble y(i,j,isoanus) = aemassm7(iloop,1,19)*convoc y(i,j,isoaais) = aemassm7(iloop,1,20)*convoc y(i,j,isoaacs) = aemassm7(iloop,1,21)*convoc y(i,j,isoacos) = aemassm7(iloop,1,22)*convoc y(i,j,isoaaii) = aemassm7(iloop,1,23)*convoc ! ug/m3 ----> molecule/cm3 ! SS mass in the 2 soluble modes y(i,j,issacs) = aemassm7(iloop,1,13)*convss y(i,j,isscos) = aemassm7(iloop,1,14)*convss ! ug/m3 ----> molecule/cm3 ! DUST mass in the 2 soluble modes and 2 insoluble y(i,j,iduacs) = aemassm7(iloop,1,15)*convdu y(i,j,iducos) = aemassm7(iloop,1,16)*convdu y(i,j,iduaci) = aemassm7(iloop,1,17)*convdu y(i,j,iducoi) = aemassm7(iloop,1,18)*convdu ! ug/m3 ----> molecule/cm3 ! number concentrations of the 7 modes y(i,j,inus_n) = aenumm7(iloop,1,1)*convnumb y(i,j,iais_n) = aenumm7(iloop,1,2)*convnumb y(i,j,iacs_n) = aenumm7(iloop,1,3)*convnumb y(i,j,icos_n) = aenumm7(iloop,1,4)*convnumb y(i,j,iaii_n) = aenumm7(iloop,1,5)*convnumb y(i,j,iaci_n) = aenumm7(iloop,1,6)*convnumb y(i,j,icoi_n) = aenumm7(iloop,1,7)*convnumb !number of particles/cm3 ----> molecule/cm3 (see comment above) #ifdef with_budgets ! For the M7 budgets, we now add the new concentrations. nzone=budg_dat(region)%nzong(i,j) ! global budget zone budm7(nzone,nzone_v,:)=budm7(nzone,nzone_v,:)+y(i,j,1:ntracet)*ye(i,j,iairm)*1000./xmair/y(i,j,iair) budm7proc(nzone,nzone_v,:)=budm7proc(nzone,nzone_v,:)+processm7(iloop,1,:)*processm7totm5*ye(i,j,iairm)*& 1000./xmair/y(i,j,iair) ! processm7totm5 is a 88-element array with unit conversion factors. !nucleation in processm7 is saved as #/cm3 within one timestep d_nucle(i,j,level)= d_nucle(i,j,level)+processm7(iloop,1,1) m_nucle_su(i,j,level)= m_nucle_su(i,j,level)+processm7(iloop,1,2) m_nucle_soa(i,j,level)= m_nucle_soa(i,j,level)+processm7(iloop,1,88) growth1_2(i,j,level)= growth1_2(i,j,level)+processm7(iloop,1,55) coag_sink_nuc(i,j,level)=coag_sink_nuc(i,j,level)+(processm7(iloop,1,3)+processm7(iloop,1,4)+processm7(iloop,1,6)& +processm7(iloop,1,8)+processm7(iloop,1,10)+processm7(iloop,1,12)& +processm7(iloop,1,14)) cond1_soa(i,j,level)= cond1_soa(i,j,level)+processm7(iloop,1,79) cond2_soa(i,j,level)= cond2_soa(i,j,level)+processm7(iloop,1,80) cond3_soa(i,j,level)= cond3_soa(i,j,level)+processm7(iloop,1,81) cond4_soa(i,j,level)= cond4_soa(i,j,level)+processm7(iloop,1,82) cond5_soa(i,j,level)= cond5_soa(i,j,level)+processm7(iloop,1,83) cond1_su(i,j,level)= cond1_su(i,j,level)+processm7(iloop,1,41) !1000./xmair/y(i,j,iair) #endif ! output m7 ! kgH2O/m3 ---> kg water per gridbox do imode=1,nsol h2o_mode(region,imode)%d3(i,j,level) = aewatm7(iloop,1,imode) *dxyp(j)*ye(i,j,idz) rwd_mode (region,imode)%d3(i,j,level) = aeradrm7(iloop,1,imode)*1e-2 ! from cm to meter end Do do imode=1,nmod dens_mode(region,imode)%d3(i,j,level) = aedensm7(iloop,1,imode)*1e3 ! from g/cm^3 --> kg/m^3 rw_mode (region,imode)%d3(i,j,level) = aeradm7(iloop,1,imode)*1e-2 ! from cm to meter end do end do end do #endif /* M7 */ ! ====================== EQSAM ! Using HNO3 / NO3_A, NH4 / NH2 and SO4. yi = 0.0 yo = 0.0 iloop = 0 do j=j1,j2 do i=i1,i2 iloop = iloop + 1 yi(iloop,1)=ye(i,j,i_temp) ! K #ifndef without_sedimentation #ifdef with_m7 s = rh(region)%d3(i,j,level) ! in cloud free part! #else s = 0.0 #endif #else s = 0.0 #endif yi(iloop,2)=max(0.001,min(0.95,s)) yi(iloop,11)=ye(i,j,i_pres)/100. ! Pa -> hPa ! yi(iloop,6)= c(nv2,ina)/Avog*1.e12 yi(iloop,6)=0. #ifdef with_m7 yi(iloop,4)=(y(i,j,iso4)+y(i,j,iso4nus)+y(i,j,iso4ais)+y(i,j,iso4acs)+y(i,j,iso4cos))/Avog*1.e12 ! molec/cm3 -> umol/m3 #else yi(iloop,4)=(y(i,j,iso4))/Avog*1.e12 ! molec/cm3 -> umol/m3 #endif yi(iloop,3)=(y(i,j,inh3)+y(i,j,inh4))/Avog*1.e12 yi(iloop,5)=(y(i,j,ihno3)+y(i,j,ino3_a))/Avog*1.e12 !yi(iloop,7)=(c(nv2,ihcl)+c(nv2,icl))/Avog*1.e12 yi(iloop,7)=0. !yi(iloop,8)=c(nv2,ik)/Avog*1.e12 yi(iloop,8)=0. !yi(iloop,9)=c(nv2,ica)/Avog*1.e12 yi(iloop,9)=0. !yi(iloop,10)=c(nv2,img)/Avog*1.e12 yi(iloop,10)=0. end do end do eqsam_option = 1 call EQSAM_V03D( yi, yo, nca, nco, eqsam_option, iloop, imax, 0, (/0,0,0,0,0,0/)) iloop = 0 do j=j1,j2 do i=i1,i2 iloop = iloop + 1 #ifdef with_budgets nzone=budg_dat(region)%nzong(i,j) ! global budget zone budeqsam(nzone,nzone_v,1) = budeqsam(nzone,nzone_v,1) + & (yo(iloop,19)*Avog*1.e-12-y(i,j,inh4)) *ye(i,j,iairm)*1000./xmair/y(i,j,iair) ! Ammonium budeqsam(nzone,nzone_v,2) = budeqsam(nzone,nzone_v,2) + & (yo(iloop,20)*Avog*1.e-12-y(i,j,ino3_a))*ye(i,j,iairm)*1000./xmair/y(i,j,iair) ! Nitrate #endif y(i,j,ihno3) = yo(iloop, 9)*Avog*1.e-12 ! umol/m3 -> molec/cm3 y(i,j,inh3) = yo(iloop,10)*Avog*1.e-12 ! umol/m3 -> molec/cm3 y(i,j,inh4) = yo(iloop,19)*Avog*1.e-12 ! umol/m3 -> molec/cm3 y(i,j,ino3_a) = yo(iloop,20)*Avog*1.e-12 ! umol/m3 -> molec/cm3 ! Sulfuric acid is unchanged: ! y(i,j,iso4) = yo(iloop,21)*Avog*1.e-12 ! The NO3_A should also get some water. M7 takes the water for the ! anoins: SO4-- and Cl- (Sea Salt). Cations are neglected. ! Just the NO3- water is used apart from the M7 water. ! Ignore the direct eqsam water ! ug/m3 -> molec/cm3 ! y(i,j,ih2opart) = yo(iloop,12)*Avog*1.e-12/xmh2o #ifdef with_m7 ! The water associated with nitrate aerosol is put into ! the soluble accumulation mode. h2o_mode(region,3)%d3(i,j,level) = h2o_mode(region,3)%d3(i,j,level) + yo(iloop,37)*dxyp(j)*ye(i,j,idz)*1.e-9 ! from ug/m^3 (according to eqsam), then take the M7 conversion (from kg/m^3) times 1.e-9. ! This EQSAM water is added to the optics water at the call of optics. #endif end do end do ! END EQSAM ! Check for negative concentrations. do n=1,ntrace do j=j1,j2 do i=i1,i2 if ( y(i,j,n) < 0. ) then !print *,' chemistry: negatives appeared in chemie' !write(6,*)i,j,level,idater !write(6,907)(y0(i,j,ispec),ispec=1,ntrace) !write(6,907)(y(i,j,ispec),ispec=1,ntrace) if ( n > ntracet ) then do ispec=ntracet+1,ntrace y(i,j,ispec)=max(0.,y(i,j,ispec)) end do else write (gol,'(a,a8,a,4i4,e12.4)') & 'chemistry: negatives for species:',names(n), ' at ', i,j,level,region, & y(i,j,n)/y(i,j,iair); call goPr y(i,j,n)=max(0.,y(i,j,n)) end if end if end do !i end do !j end do !n ! Write the tracers back to the mass_dat (rm) do n=1,ntrace #ifdef with_cariolle if ( (n==io3) .and. ( (level.ge.l_cariolle(region)) .or. with_trop_force ) ) then ! O3 is changed by Cariolle chemistry else #endif do j=j1,j2 do i=i1,i2 if ( n <= ntracet ) then rm(i,j,level,n)=y(i,j,n)/y(i,j,iair)*ye(i,j,iairm)/fscale(n) else rmc(i,j,level,n)=y(i,j,n)/y(i,j,iair)*ye(i,j,iairm)/fscale(n) end if #ifdef slopes if ( n <= ntracet ) then skf=0. if ( y(i,j,n) > eps .and. y0(i,j,n) > eps ) skf=y(i,j,n)/y0(i,j,n) rxm(i,j,level,n)=rxm(i,j,level,n)*skf rym(i,j,level,n)=rym(i,j,level,n)*skf rzm(i,j,level,n)=rzm(i,j,level,n)*skf #ifdef secmom rxxm(i,j,level,n)=rxxm(i,j,level,n)*skf rxym(i,j,level,n)=rxym(i,j,level,n)*skf rxzm(i,j,level,n)=rxzm(i,j,level,n)*skf ryym(i,j,level,n)=ryym(i,j,level,n)*skf ryzm(i,j,level,n)=ryzm(i,j,level,n)*skf rzzm(i,j,level,n)=rzzm(i,j,level,n)*skf #endif endif #endif end do end do !j,i #ifdef with_cariolle endif #endif end do !n ! For the optics calculations, ! it is assumed that nitrate aerosol is formed ! by condensation onto existing particles ! in the soluble accumulation mode (Adams et al., JGR, 2001). ! Therefore, the presence of nitrate does not increase ! the number of particles, nor does it have an impact ! on the processes described by M7. ! Above, the water associated with nitrate aerosol ! is added to the total aerosol water (h2o_mode) in the accumulation mode. ! However, this only affects the calculation of the refractive index, ! which also accounts for the nitrate mass itself. ! In order to properly account for the additional ammonium nitrate ! and the associated water in the optics routine, ! the corresponding dry and wet particle radii have to be increased as well; ! the aerosol density (dens_mode) can be left unchanged, ! because it is not used in the optics routine. ! We use the same refractive index for ammonium nitrate as for sulfate. ! This is a reasonable approximation, since, ! according to Tang (JGR, 1997), for sulfate and nitrate aerosols ! the effect of chemical composition on light scattering ! is outweighted by the size effect of the particles. ! Lowenthal et al. (Atmos. Environ., 2000) ! give a refractive index of 1.44 for sulfuric acid ! and 1.55 for ammonium nitrate. ! If desired, the refractive index of ammonium nitrate ! could be estimated at different wavelengths ! by scaling the OPAC data for sulfate, ! based on these values. ! In M7 and in the optics calculations, ! sulfate aerosols are assumed to be pure sulfuric acid, ! and we do not try to account for the presence of ammonium (bi)sulfate. ! The reason is that the presence of ammonium (bi)sulfate ! has competing effects on scattering, ! which cannot be easily accounted for in a model based on M7. ! Sulfuric acid particles are extremely hygroscopic ! and will draw significant water mass into the aerosol phase ! at any relative humidity (RH). ! If these particles are partially or completely neutralized ! by drawing ammonia from the gas phase, ! there will be an increase in particle mass due to the added ammonium ! but a decrease in particle hygroscopicity at low to moderate RH, ! and thus a decrease in particulate water mass. ! On the other hand, ammonium bisulfate and ammonium sulfate both ! have a higher refractive index than sulfuric acid (see e.g. Boucher and Anderson, 1995). ! The net result of these competing factors is that ! one mole of sulfuric acid scatters about 25% more sunlight ! than one mole of ammonium bisulfate at 80% relative humidity. ! Ammonium sulfate is intermediate between the two. ! For more information see Boucher and Anderson (1995) and Boucher et al. (JGR, 1998). ! An excellent discussion of the optical properties of ammonium (bi)sulfate ! versus sulfuric acid is also given by Adams et al. (JGR, 2001). ! ! In remote regions, especially outside the tropics, ! the contribution of methane sulfonate (MSA-) aerosol ! to the total particulate sulfure burden ! cannot be neglected. ! Like nitrate, MSA- aerosols are mainly formed by ! condensation onto exisiting particles, ! in the soluble accumulation and coarse modes ! (Saltzmann et al., JGR, 1983; Pzenny et al., J. Atmos. Chem., 1992). ! For simplicity, we assume in the model ! that MSA- is formed only in the accumulation mode. ! Moreover, it seems reasonable to use the refractive index of ! sulfuric acid (H2SO4) for MSA, as is done by Kinne et al. (2003). ! The similarity between the refractive index for MSA and H2SO4 ! is also discussed by Myhre et al. (Applied Optics, 2004). ! #ifdef with_m7 do j=j1,j2 do i=i1,i2 ! prevent division for zero for extremely low aerosol concentrations ! calculate corresponding number concentration in #/cm3 ! and test if this is larger than 1.e-10 ! the same criterion is applied in m7_equil number_conc = rm(i,j,level,iacs_n) / ( 1.e6 * dxyp(j) * ye(i,j,idz) ) if ( number_conc > 1.e-10 ) then dryvol_m7=((rwd_mode(region,3)%d3(i,j,level)*cmr2ram(3))**3.0)*pi/0.75 ! dry particle volume [m3] ! in the accumulation mode from M7 vol_nh4no3=rm(i,j,level,ino3_a)*nh4no3_factor/ & (rm(i,j,level,iacs_n)*nh4no3_density) ! ammonium-nitrate aerosol volume per particle [m3] vol_msa=rm(i,j,level,imsa)/(rm(i,j,level,iacs_n)*msa_density) ! MSA- aerosol volume per particle [m3] vol_h2o=h2o_mode(region,3)%d3(i,j,level)/(rm(i,j,level,iacs_n)*rol) ! total aerosol water per particle [m3] rw_mode(region,3)%d3(i,j,level) = ((dryvol_m7+vol_nh4no3+vol_msa+vol_h2o)*0.75/pi)**(1./3.)*ram2cmr(3) rwd_mode(region,3)%d3(i,j,level) = ((dryvol_m7+vol_nh4no3+vol_msa)*0.75/pi)**(1./3.)*ram2cmr(3) else ! don't change rwd_mode and set rw_mode equal to rwd_mode rw_mode(region,3)%d3(i,j,level) = rwd_mode(region,3)%d3(i,j,level) endif enddo enddo #endif end do LEVS ! loop over vertical levels... ! SAD phot_dat(region)%sad_cld=sad_cld_out phot_dat(region)%sad_ice=sad_ice_out phot_dat(region)%sad_aer=sad_aer_out ! perform averaging... phot_dat(region)%nsad_av = phot_dat(region)%nsad_av + float(ndyn)/float(ndyn_max) phot_dat(region)%sad_cld_av = phot_dat(region)%sad_cld_av + (float(ndyn)/float(ndyn_max))*phot_dat(region)%sad_cld phot_dat(region)%sad_ice_av = phot_dat(region)%sad_ice_av + (float(ndyn)/float(ndyn_max))*phot_dat(region)%sad_ice phot_dat(region)%sad_aer_av = phot_dat(region)%sad_aer_av + (float(ndyn)/float(ndyn_max))*phot_dat(region)%sad_aer ! cleanup deallocate(tha) deallocate(mu) deallocate(sza) deallocate(reff) deallocate(fact) deallocate(rj) deallocate(rj_new) deallocate(rr) deallocate(y) deallocate(y0) deallocate(ye, yhelp) ! eqsam deallocate(yo) deallocate(yi) ! sad deallocate(sad_cld_out,sad_ice_out,sad_aer_out) #ifdef with_m7 deallocate(presm7) deallocate(rhm7) deallocate(tempm7) deallocate(dryairm7) deallocate(h2so4m7) deallocate(aemassm7) deallocate(aenumm7) deallocate(aedensm7) deallocate(aewatm7) deallocate(aeradm7) deallocate(aeradrm7) deallocate(elvocm7) deallocate(svocm7 ) #ifdef with_budgets deallocate(processm7) #endif nullify(dxyp) #endif nullify(m) nullify(rm) nullify(rmc) #ifdef slopes nullify(rxm) nullify(rym) nullify(rzm) #ifdef secmom nullify(rxxm) nullify(rxym) nullify(rxzm) nullify(ryym) nullify(ryzm) nullify(rzzm) #endif #endif nullify(r_cloud) if ( okdebug ) write (*,*) 'end of chemistry' ! ok status = 0 END SUBROUTINE CHEMIE !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: GET_EXTRA ! ! !DESCRIPTION: ! calculate / get 2D fields needed in the chemistry routines... ! some fields (like pH) are calculated later on... !\\ !\\ ! !INTERFACE: ! SUBROUTINE GET_EXTRA( region, level, ye, is, js) ! ! !USES: ! use Binas, only : Avog, grav use dims, only : im, jm, lm use Dims, only : CheckShape use global_data, only : region_dat, mass_dat use meteodata, only : phlb_dat, m_dat, temper_dat, humid_dat, gph_dat, cc_dat, iwc_dat, lwc_dat use chem_param, only : n_extra, i_pres, i_temp, iairn, ih2on use chem_param, only : xmair, xmh2o, iairm, ilwc, iiwc, icc, iph, idz,iciwc,iclwc ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region ! region # integer, intent(in) :: level, is, js ! level and start indices ! ! !OUTPUT PARAMETERS: ! real, intent(out) :: ye(is:,js:,:) ! ! !REVISION HISTORY: ! 22 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC real,dimension(:,:,:),pointer :: phlb,m,gph,t,q,lwc,iwc,cc integer :: i,j,l real :: px,x1,x2,xice,xliq,aird integer :: imr, jmr, lmr, i1, i2, j1, j2, status ! --- begin -------------------------------- call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) ! check arguments ... call CheckShape( (/i2-i1+1, j2-j1+1, n_extra/), shape(ye ), status ) lmr = lm(region) phlb => phlb_dat(region)%data m => m_dat(region)%data t => temper_dat(region)%data q => humid_dat(region)%data gph => gph_dat(region)%data ! (:,:,1:lm(region)+1) lwc => lwc_dat(region)%data iwc => iwc_dat(region)%data cc => cc_dat(region)%data do j=j1,j2 do i=i1,i2 px = 0.5*(phlb(i,j,level)+phlb(i,j,level+1)) ye(i,j,i_pres) = px ye(i,j,i_temp) = t(i,j,level) ye(i,j,iairn ) = 7.24e16*px/t(i,j,level) ye(i,j,ih2on ) = q(i,j,level)*ye(i,j,iairn)*xmair/xmh2o ye(i,j,iairm ) = m(i,j,level) ye(i,j,iclwc ) = lwc(i,j,level) ye(i,j,iciwc ) = iwc(i,j,level) ! x1: kg water to m3 water (m3 water/ kg air) x1=lwc(i,j,level)*1.e-3 aird = ye(i,j,iairn) x2=xmair*1e3*aird/Avog ! kg/m3 (air) xliq = x1/x2 ! dimensionless number (m^3/m^3) ! avoid negatives and artificial values(1e-10 is about 0.0001 g/m3) if ( xliq < 1e-10 ) xliq=0. ye(i,j,ilwc) = xliq x1=iwc(i,j,level)*1.e-3 ! kg ice to m3 ice xice=x1/x2 ! avoid negatives and artificial values) if ( xice < 1e-10 ) xice=0. ye(i,j,iiwc) = xice ye(i,j,icc) = cc(i,j,level) ye(i,j,iph) = 0.0 ! only set in wetS when cloudy!!! ye(i,j,idz) = gph(i,j,level+1)-gph(i,j,level) !dz end do end do nullify(phlb) nullify(m) nullify(gph) nullify(t) nullify(q) nullify(lwc) nullify(iwc) nullify(cc) END SUBROUTINE GET_EXTRA !EOC END MODULE CHEMISTRY