#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 #define IF_NOTOK_MDF(action) if (status/=0) then; TRACEBACK; action; call MDF_CLose(fid,status); status=1; return; end if ! #include "tm5.inc" ! !----------------------------------------------------------------------------- ! TM5 ! !----------------------------------------------------------------------------- !BOP ! ! !MODULE: USER_OUTPUT_GENERAL ! ! !DESCRIPTION: ! ! This module provides objects and methods for the AEROCOM-2 related ! diagnostics package. It is strongly related to the setup of the ! optics module. ! A list of parameters are defined for which daily averages are output ! (so far), result is one file per day containing 50 2d-fields and ! 4 3d-fields. ! ! Life cycle(s): ! 1) output_aerocom_init: ! - if( newsrun ): - allocation and initialisation of output container ! (mixf%..., drydepos/wetdepos/emission, wdep_out) ! - initialisation of optics module ! - lookuptable ! - wavelengths ! - open output file (for today), initialise it and write attributes ! 2) output_aerocom_step: ! - gather tracer ! - add current state of model to containers ! (weighted by time since last call, --> period of validity) ! - do fluxes, surface concentrations and loads ! - do optics ! - aerocom_aopio_init: initialise and fill input parameters ! - optics::optics_calculate_aop: ! calculate optical properties (per wavelength) ! (take into account additional splitting, facilitated in ! optics::optics_calculate_aop) ! - fill containers for optical properties ! - free optics input parameters ! - fill aerocom output containers ! 3) output_aerocom_write: ! - collect_budgets: ! - use current budget containers (budemi, buddep, buddry) ! - compare with saved values and infer by that the increment ! - add to temporary output containers (wetdepos,drydepos,emission) ! - write output containers to file ! - close file ! 4) output_aerocom_done: ! - free allocated containers ! (mixf%..., drydepos/wetdepos/emission, wdep_out) ! !\\ !\\ ! !INTERFACE: ! MODULE USER_OUTPUT_GENERAL ! ! !USES: ! use go, only : gol, goErr, goPr, goLabel use dims, only : nregions use optics, only : wavelendep use meteodata, only : global_lli, levi use MDF use TM5_DISTGRID, only : dgrid, Get_DistGrid, update_halo, update_halo_iband use mo_aero_m7 , only: ncomp,nmod,naermod implicit none private ! ! !PUBLIC MEMBER FUNCTIONS: ! public :: output_general_init public :: output_general_step public :: output_general_write public :: output_general_done public :: wdep_out ! Control variables character(len=20), public :: gen_exper ! general experiment name character(len=20), public :: gen_freq ! general output frequency logical,public :: all_chemistry=.false. logical, parameter :: aai_output=.false. ! ! diagnostic CCN logical, public :: nCCNdiag = .false. ! diagnostic CCN under given supersaturations for comparison to CCNC measurements ! nCCNdiag = .false. no diagnostic CCN ! nCCNdiag = .true. with diagnostic CCN integer,public :: nSat = 0 ! Number of supersaturations for which number of CCN will be calculated real, public, allocatable :: SuperSat(:) integer,public :: n_rad = 3 ! bacchus output at n diam integer, public :: HG_rad(3)=(/50,80, 120/) integer, public :: N_diam(3)=(/50,80, 120/) ! !PRIVATE DATA MEMBERS: ! character(len=*), parameter :: mname = 'user_output_general' ! Changed naming convections to AeroCom Control 2015 experiment character(len=20), parameter :: f_dataid='general', f_model='TM5' character(len=20), parameter :: f_dimmix='global', f_dimstat='stations' ! ! !PRIVATE TYPES: ! type metafields integer :: itm5 character(len=20) :: vname character(len=64) :: lname character(len=10) :: unit character(len=10) :: positive character(len=130) :: standard_name end type metafields type metafields2 integer :: itm5 character(len=17) :: vname character(len=64) :: lname character(len=10) :: unit character(len=130) :: standard_name end type metafields2 type field3d type( metafields ) :: mf real, dimension(:,:,:), pointer :: field integer :: varid logical :: writeout end type field3d type field2d type( metafields ) :: mf real, dimension(:,:), pointer :: field integer :: varid logical :: writeout end type field2d type field1d type( metafields2 ) :: mf real, dimension(:,:), pointer :: field integer :: varid end type field1d type field0d type( metafields2 ) :: mf real, dimension(:), pointer :: field integer :: varid end type field0d type mixfile type( field3d ), dimension(:), pointer :: f3d type( field2d ), dimension(:), pointer :: f2d character(len=200) :: fname integer :: acct ! accumulation time integer :: funit ! unit number integer :: nlon ! x dimension of requested field integer :: nlat ! y dimension of requested field integer :: nlev ! z dimension of requested field real,dimension(:),pointer :: lon ! x dimension of requested field real,dimension(:),pointer :: lat ! y dimension of requested field real,dimension(:),pointer :: lev ! z dimension of requested field integer :: lonid ! x dimension id in nc integer :: latid ! y dimension id in nc integer :: levid ! z dimension id in nc integer :: timeid ! time dimension id in nc integer :: time_varid end type mixfile type(mixfile), dimension(nregions), save :: mixf integer::lon_varid integer::lat_varid integer::lev_varid ! diagnostic CCN integer, allocatable :: ind_CCN(:) integer, allocatable :: ind_HG(:) integer, allocatable :: ind_N(:) type budgetstore real, dimension(:,:,:), allocatable :: f2dslast integer :: lasttime end type budgetstore type(budgetstore), dimension(nregions), save :: drydepos, wetdepos, emission ! wavelength information type(wavelendep), dimension(:), allocatable :: wdep_out integer :: fid ! file id for IF_NOTOK_MDF macro integer :: access_mode ! netcdf-4 parallel access mode for distributed data ! generate output if true: logical :: okdebug = .true. ! reference time: integer :: time_reftime6(6) = (/2000,01,01,00,00,00/) character(len=*), parameter :: time_units = 'hours since 2000-01-01 00:00' ! sizes of arrays integer, save :: ntracer_3d, ntracer_2d integer, save :: ntracer_1dstat, ntracer_0dstat, nstations integer, save :: nextra_1dstat integer, save :: n_2d_vars=0 integer, save :: n_3d_vars=0 ! index pointers for 1d/2d/3d fields in mixf integer, save :: temp, hus, airmass, pressure integer, save :: p_gas_so4,p_liq_so4,p_elvoc,p_svoc,d_nuc,m_nuc_su,m_nuc_soa,p_svoc2d,p_elvoc2d,p_gas_so42d,p_liq_so42d integer, save :: p_elohterp,p_elo3terp,p_elohisop,p_elo3isop,p_svohterp,p_svo3terp,p_svohisop,p_svo3isop integer, save :: coag1, gr1_2, co1_soa, co1_su, co2_soa, co3_soa, co4_soa, co5_soa,gph integer, save :: so4pm1,bcpm1,orgpm1,soapm1,dupm1,sspm1,nh4pm1,no3pm1 integer, save :: so4,bc,org,du,ss,nh4,no3,n3,mmf1,mmf2,mmf3,mmf4 integer, save :: ec5503Daer, abs5503Daer, ec3503Daer, abs3503Daer integer, save :: ps , precip , sconcoa , sconcsoa , sconcbc , sconcso4 , sconcdust integer, save :: sconcss , sconcno3 , loadoa , loadsoa , loadbc , loadso4 , loaddust integer, save :: loadss , loadno3 , emioa , emibc , emiso2 , emiso4 integer, save :: emidust , emidms , emiss , dryso2 , drybc , dryso4 integer, save :: drydust , drydms , dryss , wetoa , wetsoa , wetbc , wetso2 integer, save :: emiterp , emiisop integer, save :: wetso4 , wetdust , wetdms , wetss , od550aer , od550so4, od550soa integer, save :: od550bc , od550oa , od550ss , od550dust , od550no3 , od550aerh2o integer, save :: od550lt1aer , od550lt1dust, od550lt1ss , abs550aer , ec550aer , asyaer integer, save :: ec550dryaer , abs550dryaer, asydryaer , ec550drylt1aer, abs550drylt1aer integer, save :: od440aer , od870aer , sconcmsa , dryoa , drysoa , sconcnh4 integer, save :: abs440aer , ec440dryaer , abs440dryaer integer, save :: abs870aer , ec870dryaer , abs870dryaer integer, save :: od350aer , abs350aer,ccn02 integer, save :: bso4nus, bso4ais, bso4acs, bso4cos, bbcais , bbcacs integer, save :: bbccos , bbcaii , bpomais, bpomacs, bpomcos, bpomaii integer, save :: bssacs , bsscos , bduacs , bducos , bduaci , bducoi, bno3_a integer, save :: bnus_n , bais_n , bacs_n , bcos_n , baii_n , baci_n integer, save :: bcoi_n , bh2onus, bh2oais, bh2oacs, bh2ocos integer, save :: tr2d_1, tr2d_2, tr2d_3, tr2d_4, tr2d_5, tr2d_6, tr2d_7, tr2d_8, tr2d_9, tr2d_10 integer, save :: tr2d_11, tr2d_12, tr2d_13, tr2d_14, tr2d_15, tr2d_16, tr2d_17, tr2d_18, tr2d_19 integer, save :: tr2d_20, tr2d_21, tr2d_22, tr2d_23 integer, save :: cc01, cc02, cc03, cc04, cc05, cc06, cc07 integer, save :: cc3d01, cc3d02, cc3d03, cc3d04, cc3d05, cc3d06, cc3d07 integer, save :: rw01, rw02, rw03, rw04, rw05, rw06, rw07 integer, save :: rd01, rd02, rd03, rd04, rd05, rd06, rd07 integer, save :: h2o1, h2o2, h2o3, h2o4 integer, save :: tr3d_1, tr3d_2, tr3d_3, tr3d_4, tr3d_5, tr3d_6, tr3d_7, tr3d_8, tr3d_9, tr3d_10 integer, save :: tr3d_11, tr3d_12, tr3d_13, tr3d_14, tr3d_15, tr3d_16, tr3d_17, tr3d_18, tr3d_19 integer, save :: rw3d01, rw3d02, rw3d03, rw3d04, rw3d05, rw3d06, rw3d07 integer, save :: rd3d01, rd3d02, rd3d03, rd3d04, rd3d05, rd3d06, rd3d07 integer, save :: h2o3d1, h2o3d2, h2o3d3, h2o3d4 integer, save ,dimension(2*nmod) :: radius integer, save ,dimension(ncomp) :: load integer, save ,dimension(ncomp) :: drydep integer, save ,dimension(ncomp) :: emi integer, save ,dimension(ncomp) :: wetdep integer, save ,dimension(ncomp) :: sed integer, save ,dimension(nmod) :: number integer, save ,dimension(naermod) :: masses integer, save :: ngas_output=10 integer, save, dimension(:),allocatable :: gas_output integer::itim_opt_out ! ! !REVISION HISTORY: ! 16 Nov 2010 - Achim Strunk ! 8 Dec 2016 - Tommi Bergman - Modified from AEROCOM-output ! ! !REMARKS: ! (1) compiled only if with_m7 is used. ! !EOP !------------------------------------------------------------------------ contains !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: OUTPUT_GENERAL_INIT ! ! !DESCRIPTION: Initialise various parameters, eg., ! wavelengths for optics and output containers. ! This routine should be called once per day, since fields ! are daily averages and one file per day is written. ! Additional tasks are done for newsrun==.true. . !\\ !\\ ! !INTERFACE: ! subroutine output_general_init(status, iregion) ! ! !USES: ! !use MeteoData, only : sp_dat, set use GO, only: GO_Timer_Def, GO_Timer_End, GO_Timer_Start use chem_param, only : ntracet,ntrace,names,mode_tracers,mode_start,mode_nm use chem_param, only : inus_n, iso4nus, isoanus,nmod use chem_param, only : iais_n, iso4ais, ibcais, ipomais, isoaais use chem_param, only : iacs_n, iso4acs, ibcacs, ipomacs, isoaacs, issacs, iduacs use chem_param, only : icos_n, iso4cos, ibccos, ipomcos, isoacos, isscos, iducos use chem_param, only : iaii_n, ibcaii, ipomaii, isoaaii use chem_param, only : iaci_n, icoi_n, iduaci, iducoi use chem_param, only : ino3_a, inh4, imsa use chem_param, only : Kap_su, Kap_pom, Kap_soa, Kap_bc, Kap_du, Kap_ss, Kap_na2so4, Kap_no3, Kap_msa !use chem_param, only : io3, idms, imsa, inh4, iterp, ioh, ino3, ielvoc, isvoc, iisop,iso2,iso4 use chem_param, only : io3, ih2o2, ieth, irooh, idms, imsa, imgly, inox, ico, ipar, ipan use chem_param, only : ich2o, iald2, iisop, iso4, irn222, ino3_a, ich4, ich3o2h, iole, iorgntr, inh3 use chem_param, only : ic2h6, iethoh, ic3h8, iterp, iacet, iispd, ihono, ich3o2no2, ino, iho2, ich3o2 use chem_param, only : in2o5, ihno4, ic2o3, iror, irxpar, ixo2, ixo2n, inh2, ih2opart,ic3h7o2, ihypro2 use chem_param, only : isvoc, iso2, inh4, ipb210,io3s, ich3oh, ihcooh , ioh, ino2, ino3, iaco2, inh2o2, ielvoc USE mo_aero_m7, ONLY : nmod, nsol use dims, only : nregions, newsrun, idatee, idatei ,sdate_simulation use dims, only : region_name, parent, itau use dims, only : xbeg, xend, ybeg, yend, dx, dy, xref, yref use dims, only : zbeg, zend, dz, zref use global_data, only : region_dat, outdir use datetime, only : tau2date, date2tau use budget_global, only : nbud_vg use partools, only : MPI_INFO_NULL, localComm use optics, only : Optics_Init use datetime, only : date2tau use dims, only : mlen, idatei use NetCDF, only : NF90_Def_Var,NF90_put_var ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !INPUT PARAMETERS: ! integer, intent(in), optional :: iregion ! ! !REVISION HISTORY: ! Nov 2010 - Achim Strunk ! 16 Dec 2016 - Tommi Bergman KNMI ! - Modified from AEROCOM output ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/output_general_init' ! --- local ------------------------------ integer :: imr, jmr, lmr, access_mode_sta integer :: lat_dimid,lon_dimid,lev_dimid,kappa_dimid integer :: year_varid,month_varid,kso4_varid, kss_varid,kpom_varid,kbc_varid,kdu_varid,kno3_varid,kna2so4_varid,ksoa_varid,kmsa_varid integer :: region, varid integer :: io, istat, i, j, n, sc, iSat integer :: i1, i2, j1, j2 integer, dimension(6) :: idater character(len=10) :: idates character(len=16) :: lidates character(len=3) :: cwavel real :: reftime integer :: itaucur, itauref, time_shift real :: dlat, dlon integer :: mlength integer :: itracer,index_aer,imass,iout character(len=14) :: str, str2 character(len=28) :: str_reftime ! --- begin ------------------------------- call goLabel(rname) ! access mode for distributed data (2d and 3d) #ifdef MPI #ifdef with_netcdf4_par access_mode = MDF_COLLECTIVE #else write(gol,'("General aerosol output requires netcdf4 with parallel access enabled")') ; call goErr TRACEBACK status=1; return #endif #else access_mode = MDF_INDEPENDENT #endif ! for station access_mode_sta = MDF_INDEPENDENT ! initialize (only once) if( newsrun ) then call GO_Timer_Def( itim_opt_out, 'optics-output', status ) IF_NOTOK_RETURN(status=1) ! ensure that required meteo is loaded: ! call Set( sp_dat(region), status, used=.true. ) ! set wavelength information ! wl: wavelength in microns ! split: whether to split into contributions from ! M7 constituents (incl. water) allocate( wdep_out( 4 ) ) wdep_out(1)%wl = 0.550 ; wdep_out(1)%split = .true. ; wdep_out(1)%insitu = .true. wdep_out(2)%wl = 0.440 ; wdep_out(2)%split = .false. ; wdep_out(2)%insitu = .true. wdep_out(3)%wl = 0.870 ; wdep_out(3)%split = .false. ; wdep_out(3)%insitu = .true. wdep_out(4)%wl = 0.350 ; wdep_out(4)%split = .false. ; wdep_out(4)%insitu = .false. ! get the optics code prepared call Optics_Init(size(wdep_out), wdep_out, status ) IF_NOTOK_RETURN(status=1) ! ----------------------- ! parameters needed to reference the different 1d/2d/3d-fields ! (in order to avoid errors in referencing) ! finally, this list here simply determines the order in the output files ntracer_3d = 160+nSat ! list removed, not needed anymore. ! 2d-output list needed for the moment ! optics needs to be changed first, since it relies on having optical variables ! using indices 38-> ntracer_2d = 122 ps = 1 ! 2d precip = 2 ! 2d sconcoa = 3 ! 2d sconcbc = 4 ! 2d sconcso4 = 5 ! 2d sconcdust = 6 ! 2d sconcss = 7 ! 2d sconcno3 = 8 ! 2d sconcnh4 = 9 ! 2d sconcmsa = 10 ! 2d loadoa = 11 ! 2d loadbc = 12 ! 2d loadso4 = 13 ! 2d loaddust = 14 ! 2d loadss = 15 ! 2d loadno3 = 16 ! 2d emioa = 17 ! 2d emibc = 18 ! 2d emiso2 = 19 ! 2d emiso4 = 20 ! 2d emidust = 21 ! 2d emidms = 22 ! 2d emiss = 23 ! 2d dryso2 = 24 ! 2d dryoa = 25 ! 2d drybc = 26 ! 2d dryso4 = 27 ! 2d drydust = 28 ! 2d drydms = 29 ! 2d dryss = 30 ! 2d wetoa = 31 ! 2d wetbc = 32 ! 2d wetso2 = 33 ! 2d wetso4 = 34 ! 2d wetdust = 35 ! 2d wetdms = 36 ! 2d wetss = 37 ! 2d ! --- from here onwards keep consistent with order in optics (optics.F90) ! --- begin split order od550aer = 38 ! 2d od550so4 = 39 ! 2d od550bc = 40 ! 2d od550oa = 41 ! 2d od550soa = 42 ! 2d od550ss = 43 ! 2d od550dust = 44 ! 2d od550no3 = 45 ! 2d od550aerh2o = 46 ! 2d od550lt1aer = 47 ! 2d od550lt1dust = 48 ! 2d od550lt1ss = 49 ! 2d ! --- end split order abs550aer = 50 ! 2d asyaer = 51 ! 2d ec550aer = 52 ! 2d ! --- begin in-situ data order ec550dryaer = 53 ! 2d abs550dryaer = 54 ! 2d asydryaer = 55 ! 2d ec550drylt1aer = 56 ! 2d abs550drylt1aer = 57 ! 2d ! --- end in-situ data order ! od440aer = 58 ! 2d abs440aer = 59 ! 2d ec440dryaer = 60 ! 2d abs440dryaer = 61 ! 2d ! od870aer = 62 ! 2d abs870aer = 63 ! 2d ec870dryaer = 64 ! 2d abs870dryaer = 65 ! 2d ! od350aer = 66 ! 2d abs350aer = 67 ! 2d ! tr2d_1=68 tr2d_2=69 tr2d_3=70 tr2d_4=71 tr2d_5=72 tr2d_6=73 tr2d_7=74 tr2d_8=75 tr2d_9=76 tr2d_10=77 tr2d_11=78 tr2d_12=79 tr2d_13=80 tr2d_14=81 tr2d_15=82 tr2d_16=83 tr2d_17=84 tr2d_18=85 tr2d_19=86 tr2d_20=87 tr2d_21=88 tr2d_22=89 tr2d_23=90 cc01=91 cc02=92 cc03=93 cc04=94 cc05=95 cc06=96 cc07=97 h2o1=98 h2o2=99 h2o3=100 h2o4=101 rw01=102 rw02=103 rw03=104 rw04=105 rw05=106 rw06=107 rw07=108 rd01=109 rd02=110 rd03=111 rd04=112 sconcsoa=113 loadsoa=114 drysoa=115 wetsoa=116 emiterp=117 emiisop=118 p_svoc2d=119 p_elvoc2d=120 p_liq_so42d=121 p_gas_so42d=122 end if regionloop: do region = 1, nregions ! if region given, cycle if other region! if (present(iregion)) then if(iregion /= region) cycle regionloop endif imr = global_lli(region)%nlon jmr = global_lli(region)%nlat lmr = levi%nlev if( nbud_vg /= lmr ) then write(gol,*)'output_general_init: nbud_vg /= lmr'; call goErr write(gol,*)'output_general_init: YOU MUST ADD THE PROJ/BUDGET10 TO SRC CODE !!!!!'; call goErr TRACEBACK status=1; return end if ! initialize the output: if( newsrun ) then if (all_chemistry) then ngas_output=57 allocate(gas_output(ngas_output)) gas_output=(/& io3, ih2o2, ieth, irooh, idms, imsa, imgly, inox, ico, ipar, ipan, iso2, inh4, ipb210, & ich2o, iald2, iisop, iso4, irn222, ino3_a, ich4, ich3o2h, iole, iorgntr, inh3, io3s, ich3oh, ihcooh, & ic2h6, iethoh, ic3h8, iterp, iacet, iispd, ihono, ich3o2no2, ino, iho2, ich3o2, ioh, ino2, ino3, & in2o5, ihno4, ic2o3, iror, irxpar, ixo2, ixo2n, inh2, ih2opart,ic3h7o2, ihypro2, iaco2, inh2o2, ielvoc, & isvoc/) else allocate(gas_output(ngas_output)) gas_output=(/io3,idms,iso4,iso2,iterp,ioh,ino3,ielvoc,isvoc,iisop /) end if call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) !timeidx=1 ! target structures for output allocate(mixf (region)%f3d(ntracer_3d)) allocate(mixf (region)%f2d(ntracer_2d)) ! allocate data holder for maximum amount of tracers do i=1,ntracer_3d allocate( mixf(region)%f3d(i)%field(i1:i2,j1:j2,lmr) ) end do do i=1,ntracer_2d allocate( mixf(region)%f2d(i)%field(i1:i2,j1:j2)) end do ! METEO variables call add_Variable(region,-1 , 'temp ', 'Temperature ' , 'K ',3,temp,status ) call add_Variable(region,-1 , 'hus ', 'Specific Humidity ' , '1 ',3,hus,status ) call add_Variable(region,-1 , 'airmass ', 'Air Mass ' , 'kg m-2 ',3,airmass,status ) call add_Variable(region,-1 ,'pres ', 'Pressure ' , 'Pa ',3,pressure,status ) !call add_Variable(region,-1 ,'prod_gas_so4 ', 'production of so4 ' , 'kg m-2s-1',3,p_gas_so4,status ) !call add_Variable(region,-1 ,'prod_liq_so4 ', 'production of so4 ' , 'kg m-2s-1',3,p_liq_so4,status ) !call add_Variable(region,-1 ,'prod_elvoc ', 'production of elvoc ' , 'kg m-2s-1',3,p_elvoc,status ) !call add_Variable(region,-1 ,'prod_svoc ', 'production of svoc ' , 'kg m-2s-1',3,p_svoc,status ) !call add_Variable(region,-1 ,'p_el_ohterp ', 'production of elvoc, ohterp ' , 'kg m-2s-1',3,p_elohterp,status ) !call add_Variable(region,-1 ,'p_el_o3terp ', 'production of elvoc, o3terp ' , 'kg m-2s-1',3,p_elo3terp,status ) !call add_Variable(region,-1 ,'p_el_ohisop ', 'production of elvoc, ohisop ' , 'kg m-2s-1',3,p_elohisop,status ) !call add_Variable(region,-1 ,'p_el_o3isop ', 'production of elvoc, o3isop ' , 'kg m-2s-1',3,p_elo3isop,status ) !call add_Variable(region,-1 ,'p_sv_ohterp ', 'production of svoc, ohterp ' , 'kg m-2s-1',3,p_svohterp,status ) !call add_Variable(region,-1 ,'p_sv_o3terp ', 'production of svoc, o3terp ' , 'kg m-2s-1',3,p_svo3terp,status ) !call add_Variable(region,-1 ,'p_sv_ohisop ', 'production of svoc, ohisop ' , 'kg m-2s-1',3,p_svohisop,status ) !call add_Variable(region,-1 ,'p_sv_o3isop ', 'production of svoc, o3isop ' , 'kg m-2s-1',3,p_svo3isop,status ) !call add_Variable(region,-1 ,'d_nuc ', 'production of particles ' , ' cm-3s-1',3,d_nuc,status ) !call add_Variable(region,-1 ,'m_nuc_su ', 'nucleation sulfate ' , ' molec cm-3s-1',3,m_nuc_su,status ) !call add_Variable(region,-1 ,'m_nuc_soa ', 'nucleation ELVOC ' , ' µg m-3s-1',3,m_nuc_soa,status ) !call add_Variable(region,-1 ,'gr1_2 ', 'growth of mode 1 into 2 ' , ' cm-3s-1',3,gr1_2,status ) !call add_Variable(region,-1 ,'coag1 ', 'coagulation sink of mode 1 ' , ' cm-3s-1',3,coag1,status ) !call add_Variable(region,-1 ,'co1_soa ', 'condensation of ELVOC to mode 1 ' , ' µg m-3s-1',3,co1_soa,status ) !call add_Variable(region,-1 ,'co2_soa ', 'condensation of ELVOC to mode 2 ' , ' µg m-3s-1',3,co2_soa,status ) !call add_Variable(region,-1 ,'co3_soa ', 'condensation of ELVOC to mode 3 ' , ' µg m-3s-1',3,co3_soa,status ) !call add_Variable(region,-1 ,'co4_soa ', 'condensation of ELVOC to mode 4 ' , ' µg m-3s-1',3,co4_soa,status ) !call add_Variable(region,-1 ,'co5_soa ', 'condensation of ELVOC to mode 5 ' , ' µg m-3s-1',3,co5_soa,status ) !call add_Variable(region,-1 ,'co1_su ', 'condensation of sulfate to mode 1 ' , ' molec cm-3s-1',3,co1_su,status ) call add_Variable(region,-1 ,'mf1 ', 'mass of so4 in pm1 ' , ' ug m-3',3,mmf1,status ) call add_Variable(region,-1 ,'mf2 ', 'mass of BC in pm1 ' , ' ug m-3',3,mmf2,status ) call add_Variable(region,-1 ,'mf3 ', 'mass of ORG in pm1 ' , ' ug m-3',3,mmf3,status ) call add_Variable(region,-1 ,'mf4 ', 'mass of SOA in pm1 ' , ' ug m-3',3,mmf4,status ) call add_Variable(region,-1 ,'SO4pm1 ', 'mass of so4 in pm1 ' , ' ug m-3',3,so4pm1,status ) call add_Variable(region,-1 ,'BCpm1 ', 'mass of BC in pm1 ' , ' ug m-3',3,bcpm1,status ) call add_Variable(region,-1 ,'ORGpm1 ', 'mass of ORG in pm1 ' , ' ug m-3',3,orgpm1,status ) call add_Variable(region,-1 ,'SOApm1 ', 'mass of SOA in pm1 ' , ' ug m-3',3,soapm1,status ) call add_Variable(region,-1 ,'SSpm1 ', 'mass of SS in pm1 ' , ' ug m-3',3,sspm1,status ) call add_Variable(region,-1 ,'DUpm1 ', 'mass of DU in pm1 ' , ' ug m-3',3,dupm1,status ) call add_Variable(region,-1 ,'NH4pm1 ', 'mass of NH4 in pm1 ' , ' ug m-3',3,nh4pm1,status ) call add_Variable(region,-1 ,'NO3pm1 ', 'mass of NO3 in pm1 ' , ' ug m-3',3,no3pm1,status ) call add_Variable(region,-1 ,'SO4 ', 'mass of SO4 ' , ' ug S m-3',3,so4,status ) call add_Variable(region,-1 ,'BC ', 'mass of BC ' , ' ug C m-3',3,bc,status ) call add_Variable(region,-1 ,'ORG ', 'mass of ORG ' , ' ug C m-3',3,org,status ) call add_Variable(region,-1 ,'SS ', 'mass of SS ' , ' ug m-3',3,ss,status ) call add_Variable(region,-1 ,'DU ', 'mass of DU ' , ' ug m-3',3,du,status ) call add_Variable(region,-1 ,'NH4 ', 'mass of NH4 ' , ' ug m-3',3,nh4,status ) call add_Variable(region,-1 ,'NO3 ', 'mass of NO3 ' , ' ug m-3',3,no3,status ) IF (nCCNdiag) THEN ALLOCATE(ind_CCN(nSat)) DO iSat=1,nSat !write(1111,*)SuperSat(iSat)*1000.e0 if (supersat(isat)<0.0006)then !WRITE(str,'("CCN",I0.3)') SuperSat(iSat)*10000.e0 !CCN005 (0.0005*1e4=5) WRITE(str,'("CCN005")') !SuperSat(iSat)*10000.e0 !CCN005 (0.0005*1e4=5) else if (supersat(isat)<0.0008)then WRITE(str,'("CCN0075")') !SuperSat(iSat)*100000.e0!CCN0075 (0.00075*1e5=75 else if (supersat(isat)<0.0012)then WRITE(str,'("CCN01")') !SuperSat(iSat)*1000.e0 !CCN01 (0.001*1e3=1) else if (supersat(isat)<0.0016)then WRITE(str,'("CCN015")') !SuperSat(iSat)*10000.e0!CCN015 0.0015*1e4=15) else !write(1112,*)SuperSat(iSat)*1000.e0 WRITE(str,'("CCN",I0.2)') nint(SuperSat(iSat)*1000.e0) !CCN0[23456789] & CCN10 end if !write(1112,*) str , int(SuperSat(iSat)*1000.e0) , nint(SuperSat(iSat)*1000.e0), (SuperSat(iSat)*1000.e0) WRITE(str2,'("CCN at ",F4.2,"%")') SuperSat(iSat)*100.e0 call add_Variable(region,-1 ,str, str2, ' cm-3',3,ind_CCN(iSat),status ) ENDDO allocate(ind_HG(n_rad)) allocate(ind_N(n_rad)) do isat=1,n_rad if (HG_rad(isat)>99) then WRITE(str,'("HG",I3)') HG_rad(isat) else WRITE(str,'("HG",I2)') HG_rad(isat) end if WRITE(str2,'("HG at ",I3,"nm")') HG_rad(isat) call add_variable(region,-1 ,str, str2, '1',3,ind_HG(iSat),status ) if (N_diam(isat)>99)then WRITE(str,'("N",I3)') N_diam(isat) else WRITE(str,'("N",I2)') N_diam(isat) end if !write(1000,*)str WRITE(str2,'("N at ",I3,"nm")') N_diam(isat) call add_variable(region,-1 ,str, str2, ' cm-3',3,ind_N(iSat),status ) end do ENDIF call add_variable(region,-1 ,'N3', 'N at 3nm', ' cm-3',3,n3,status ) call add_Variable(region,-1 ,'gph', 'Geopotential height', 'm',3,gph,status ) call add_Variable(region,-1 ,'CCN', 'CCN at 0.2%', '#/cm3',3,ccn02,status ) ! mass concentrations of different modes imass=0 do i=1,nmod do j=1,mode_nm(i)!ncomp !get tracer id from mode_tracers array itracer=mode_tracers(j,i) ! if there is a tracer, add variable to output array !if (itracer>0)then imass=imass+1 call add_Variable(region,itracer ,'M_'//names(itracer), 'mass concentration of '//names(itracer), 'kg m-3',3,index_aer,status ) !end if ! save output array id for masses masses(imass)=index_aer end do end do call add_Variable(region,imsa ,'M_MSA', 'mass concentration of MSA (ACS)', 'kg m-3',3,index_aer,status ) call add_Variable(region,ino3_a ,'M_NO3', 'mass concentration of NO3 (ACS)', 'kg m-3',3,index_aer,status ) call add_Variable(region,inh4 ,'M_NH4', 'mass concentration of NH4 (ACS)', 'kg m-3',3,index_aer,status ) ! Mode number concentrations do i=1,nmod !get tracer id from mode_tracers array ! 0-index of mode_tracers has the number for each mode itracer=mode_tracers(0,i) ! if there is a tracer, add variable to output array if (itracer>0)then ! call add_Variable(region,itracer ,'N_'//names(itracer)(1:3), 'number concentration of '//names(itracer), '1 m-3',3,index_aer,status ) end if ! save output array id for number number(i)=index_aer end do ! gas-phase output !!$ do i=1,ngas_output !!$ itracer=gas_output(i) !!$ call add_Variable(region,itracer ,'GAS_'//names(itracer), 'gas phase concentration of '//names(itracer), 'kg m-3',3,index_aer,status ) !!$ end do !!$ !!$ ! particle water !!$ do i=1,nsol !!$ !write(iout,'(I2.2)') i !!$ itracer=mode_start(i) !!$ call add_Variable(region,-1 , 'aerh2o3d_'//names(itracer)(1:3), 'concentration of aerosol water of mode '//names(itracer)(1:3) , 'kg m-3',3,index_aer,status ) !!$ !!$ end do ! Radius of the particles ! wet radius for all !!$ do i=1,nmod !!$ itracer=mode_start(i) !!$ call add_Variable(region,-1 ,'RWET_'//names(itracer)(1:3), 'wet radius of '//names(itracer)(1:3), 'm',3,index_aer,status ) !!$ end do !!$ ! dry radius only for soluble, as rwet=rdry for insoluble !!$ do i=1,nsol !!$ itracer=mode_start(i) !!$ !!$ call add_Variable(region,-1 ,'RDRY_'//names(itracer)(1:3), 'dry radius of '//names(itracer)(1:3), 'm',3,index_aer,status ) !!$ !!$ end do ! 2d data ! For now it is exactly the same as in AEROCOM-output, mostly due to optics depending on ! output indices now in place mixf(region)%f2d(ps )%mf = metafields( -1 , 'ps ', 'Surface Air Pressure ' , 'Pa ', '', 'surface_air_pressure' ) mixf(region)%f2d(precip )%mf = metafields( -1 , 'precip ', 'Precipitation ' , 'kg m-2 s-1', '', 'precipitation_flux' ) mixf(region)%f2d(sconcoa )%mf = metafields( -1 , 'sconcoa ', 'Surface Concentration POM ' , 'kg m-3 ', '', 'mass_concentration_of_particulate_organic_matter_dry_aerosol_in_air' ) mixf(region)%f2d(sconcbc )%mf = metafields( -1 , 'sconcbc ', 'Surface Concentration BC ' , 'kg m-3 ', '', 'mass_concentration_of_black_carbon_dry_aerosol_in_air' ) mixf(region)%f2d(sconcso4 )%mf = metafields( -1 , 'sconcso4 ', 'Surface Concentration SO4 ' , 'kg m-3 ', '', 'mass_concentration_of_sulfate_dry_aerosol_in_air' ) mixf(region)%f2d(sconcdust )%mf = metafields( -1 , 'sconcdust ', 'Surface Concentration Dust ' , 'kg m-3 ', '', 'mass_concentration_of_dust_dry_aerosol_in_air' ) mixf(region)%f2d(sconcss )%mf = metafields( -1 , 'sconcss ', 'Surface Concentration SS ' , 'kg m-3 ', '', 'mass_concentration_of_seasalt_dry_aerosol_in_air' ) mixf(region)%f2d(sconcno3 )%mf = metafields( -1 , 'sconcno3 ', 'Surface Concentration NO3 ' , 'kg m-3 ', '', 'mass_concentration_of_nitrate_dry_aerosol_in_air' ) mixf(region)%f2d(sconcnh4 )%mf = metafields( -1 , 'sconcnh4 ', 'Surface Concentration NH4 ' , 'kg m-3 ', '', 'mass_concentration_of_ammonium_dry_aerosol_in_air' ) mixf(region)%f2d(sconcmsa )%mf = metafields( -1 , 'sconcmsa ', 'Surface Concentration MSA ' , 'kg m-3 ', '', 'mass_concentration_of_methane_sulfonic_acid_in_air' ) mixf(region)%f2d(loadoa )%mf = metafields( -1 , 'loadoa ', 'Load of POM ' , 'kg m-2 ', '', 'atmosphere_mass_content_of_particulate_organic_matter_dry_aerosol' ) mixf(region)%f2d(loadbc )%mf = metafields( -1 , 'loadbc ', 'Load of BC ' , 'kg m-2 ', '', 'atmosphere_mass_content_of_black_carbon_dry_aerosol' ) mixf(region)%f2d(loadso4 )%mf = metafields( -1 , 'loadso4 ', 'Load of SO4 ' , 'kg m-2 ', '', 'atmosphere_mass_content_of_sulfate_dry_aerosol' ) mixf(region)%f2d(loaddust )%mf = metafields( -1 , 'loaddust ', 'Load of Dust ' , 'kg m-2 ', '', 'atmosphere_mass_content_of_dust_dry_aerosol' ) mixf(region)%f2d(loadss )%mf = metafields( -1 , 'loadss ', 'Load of SS ' , 'kg m-2 ', '', 'atmosphere_mass_content_of_seasalt_dry_aerosol' ) mixf(region)%f2d(loadno3 )%mf = metafields( -1 , 'loadno3 ', 'Load of NO3 ' , 'kg m-2 ', '', 'atmosphere_mass_content_of_nitrate_dry_aerosol' ) mixf(region)%f2d(emioa )%mf = metafields( -1 , 'emioa ', 'Total Emission of POM ' , 'kg m-2 s-1', 'up', 'tendency_of_atmosphere_mass_content_of_primary_particulate_organic_matter_dry_aerosol_due_to_emission' ) mixf(region)%f2d(emibc )%mf = metafields( -1 , 'emibc ', 'Total Emission of BC ' , 'kg m-2 s-1', 'up', 'tendency_of_atmosphere_mass_content_of_black_carbon_dry_aerosol_due_to_emission' ) mixf(region)%f2d(emiso2 )%mf = metafields( -1 , 'emiso2 ', 'Total Emission of SO2 ' , 'kg m-2 s-1', 'up', 'tendency_of_atmosphere_mass_content_of_sulfur_dioxide_due_to_emission' ) mixf(region)%f2d(emiso4 )%mf = metafields( -1 , 'emiso4 ', 'Total Direct Emission of SO4 ' , 'kg m-2 s-1', 'up', 'tendency_of_atmosphere_mass_content_of_sulfate_dry_aerosol_due_to_emission' ) mixf(region)%f2d(emidust )%mf = metafields( -1 , 'emidust ', 'Total Emission of Dust ' , 'kg m-2 s-1', 'up', 'tendency_of_atmosphere_mass_content_of_dust_dry_aerosol_due_to_emission' ) mixf(region)%f2d(emidms )%mf = metafields( -1 , 'emidms ', 'Total Emission of DMS ' , 'kg m-2 s-1', 'up', 'tendency_of_atmosphere_mass_content_of_dimethyl_sulfide_due_to_emission' ) mixf(region)%f2d(emiss )%mf = metafields( -1 , 'emiss ', 'Total Emission of SeaSalt ' , 'kg m-2 s-1', 'up', 'tendency_of_atmosphere_mass_content_of_seasalt_dry_aerosol_due_to_emission' ) mixf(region)%f2d(dryso2 )%mf = metafields( -1 , 'dryso2 ', 'Dry Deposition of SO2 ' , 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_sulfur_dioxide_due_to_dry_deposition' ) mixf(region)%f2d(dryoa )%mf = metafields( -1 , 'dryoa ', 'Dry Deposition of POM ' , 'kg m-2 s-1', 'down', & 'tendency_of_atmosphere_mass_content_of_primary_particulate_organic_matter_dry_aerosol_due_to_dry_deposition' ) mixf(region)%f2d(drybc )%mf = metafields( -1 , 'drybc ', 'Dry Deposition of BC ' , 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_black_carbon_dry_aerosol_due_to_dry_deposition' ) mixf(region)%f2d(dryso4 )%mf = metafields( -1 , 'dryso4 ', 'Dry Deposition of SO4 ' , 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_sulfate_dry_aerosol_due_to_dry_deposition' ) mixf(region)%f2d(drydust )%mf = metafields( -1 , 'drydust ', 'Dry Deposition of Dust ' , 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_dust_dry_aerosol_due_to_dry_deposition' ) mixf(region)%f2d(drydms )%mf = metafields( -1 , 'drydms ', 'Dry Deposition of DMS ' , 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_dimethyl_sulfide_due_to_dry_deposition' ) mixf(region)%f2d(dryss )%mf = metafields( -1 , 'dryss ', 'Dry Deposition of SeaSalt ' , 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_seasalt_dry_aerosol_due_to_dry_deposition' ) mixf(region)%f2d(wetoa )%mf = metafields( -1 , 'wetoa ', 'Wet Deposition of POM ' , 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_particulate_organic_matter_dry_aerosol_due_to_wet_deposition' ) mixf(region)%f2d(wetbc )%mf = metafields( -1 , 'wetbc ', 'Wet Deposition of BC ' , 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_black_carbon_dry_aerosol_due_to_wet_deposition' ) mixf(region)%f2d(wetso2 )%mf = metafields( -1 , 'wetso2 ', 'Wet Deposition of SO2 ' , 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_sulfur_dioxide_due_to_wet_deposition' ) mixf(region)%f2d(wetso4 )%mf = metafields( -1 , 'wetso4 ', 'Wet Deposition of SO4 ' , 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_sulfate_dry_aerosol_due_to_wet_deposition' ) mixf(region)%f2d(wetdust )%mf = metafields( -1 , 'wetdust ', 'Wet Deposition of Dust ' , 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_dust_dry_aerosol_due_to_wet_deposition' ) mixf(region)%f2d(wetdms )%mf = metafields( -1 , 'wetdms ', 'Wet Deposition of DMS ' , 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_dimethyl_sulfide_due_to_wet_deposition' ) mixf(region)%f2d(wetss )%mf = metafields( -1 , 'wetss ', 'Wet Deposition of Seasalt ' , 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_seasalt_dry_aerosol_due_to_wet_deposition' ) mixf(region)%f2d(od550aer )%mf = metafields( -1 , 'od550aer ', 'AOD@550nm ' , '1 ', '', 'atmosphere_optical_thickness_due_to_ambient_aerosol' ) mixf(region)%f2d(od550so4 )%mf = metafields( -1 , 'od550so4 ', 'SO4 AOD@550nm ' , '1 ', '', 'atmosphere_optical_thickness_due_to_sulfate_ambient_aerosol' ) mixf(region)%f2d(od550bc )%mf = metafields( -1 , 'od550bc ', 'Black Carbon AOD@550nm ' , '1 ', '', 'atmosphere_optical_thickness_due_to_black_carbon_ambient_aerosol' ) mixf(region)%f2d(od550oa )%mf = metafields( -1 , 'od550oa ', 'POM AOD@550nm ' , '1 ', '', 'atmosphere_optical_thickness_due_to_particulate_organic_matter_ambient_aerosol' ) mixf(region)%f2d(od550soa )%mf = metafields( -1 , 'od550soa ', 'SOA AOD@550nm ' , '1 ', '', 'atmosphere_optical_thickness_due_to_particulate_secondary_organic_aerosol_ambient_aerosol' ) mixf(region)%f2d(od550ss )%mf = metafields( -1 , 'od550ss ', 'SeaSalt AOD@550nm ' , '1 ', '', 'atmosphere_optical_thickness_due_to_seasalt_ambient_aerosol' ) mixf(region)%f2d(od550dust )%mf = metafields( -1 , 'od550dust ', 'Dust AOD@550nm ' , '1 ', '', 'atmosphere_optical_thickness_due_to_dust_ambient_aerosol' ) mixf(region)%f2d(od550no3 )%mf = metafields( -1 , 'od550no3 ', 'Nitrate AOD@550nm ' , '1 ', '', 'atmosphere_optical_thickness_due_to_nitrate_ambient_aerosol' ) mixf(region)%f2d(od550aerh2o )%mf = metafields( -1 , 'od550aerh2o ', 'Aerosol Water AOD@550nm ' , '1 ', '', 'atmosphere_optical_thickness_due_to_water_in_ambient_aerosol' ) mixf(region)%f2d(od550lt1aer )%mf = metafields( -1 , 'od550lt1aer ', 'Fine Mode AOD@550nm ' , '1 ', '', 'atmosphere_optical_thickness_due_to_pm1_ambient_aerosol' ) mixf(region)%f2d(od550lt1dust)%mf = metafields( -1 , 'od550lt1dust', 'Fine Mode Dust AOD@550nm ' , '1 ', '', 'atmosphere_optical_thickness_due_to_dust_pm1_ambient_aerosol' ) mixf(region)%f2d(od550lt1ss)%mf = metafields( -1 , 'od550lt1ss ', 'Fine Mode SeaSalt AOD@550nm ' , '1 ', '', 'atmosphere_optical_thickness_due_to_seasalt_pm1_ambient_aerosol' ) mixf(region)%f2d(abs550aer )%mf = metafields( -1 , 'abs550aer ', 'Absorption AOD@550nm ' , '1 ', '', 'atmosphere_absorption_optical_thickness_due_to_ambient_aerosol' ) mixf(region)%f2d(asyaer )%mf = metafields( -1 , 'asyaer ', 'Asymmetry Parameter ' , '1 ', '', 'atmosphere_asymmetry_parameter_ambient_aerosol' ) mixf(region)%f2d(ec550aer )%mf = metafields( -1 , 'ec550aer ', 'Surface Ambient Aerosol Extinction@550nm' , 'm-1 ', '', 'surface_extinction_due_to_ambient_aerosol' ) mixf(region)%f2d(ec550dryaer )%mf = metafields( -1 , 'ec550dryaer ', 'Surface Dry Aerosol Extinction@550 nm','m-1 ', '', 'surface_extinction_due_to_dry_aerosol' ) mixf(region)%f2d(abs550dryaer)%mf = metafields( -1 , 'abs550dryaer', 'Surface Dry Aerosol Absorption@550 nm','m-1 ', '', 'surface_absorption_due_to_dry_aerosol' ) mixf(region)%f2d(asydryaer )%mf = metafields( -1 , 'asydryaer ', 'Surface Dry Aersol Asymmetry Parameter' , '1 ', '', 'surface_asymmetry_parameter_dry_aerosol' ) mixf(region)%f2d(ec550drylt1aer )%mf = metafields( -1 , 'ec550drylt1aer ', 'Surface Fine Mode Dry Aerosol Extinction@550 nm' , 'm-1', '', 'surface_extinction_due_to_pm1_dry_aerosol' ) mixf(region)%f2d(abs550drylt1aer)%mf = metafields( -1 , 'abs550drylt1aer', 'Surface Fine Mode Dry Aerosol Absorption@550 nm' , 'm-1', '', 'surface_absorption_due_to_pm1_dry_aerosol' ) mixf(region)%f2d(od440aer )%mf = metafields( -1 , 'od440aer ', 'AOD@440nm ' , '1 ', '', 'atmosphere_optical_thickness_due_to_ambient_aerosol' ) mixf(region)%f2d(abs440aer )%mf = metafields( -1 , 'abs440aer ', 'Absorption AOD@440nm ' , '1 ', '', 'atmosphere_absorption_optical_thickness_due_to_ambient_aerosol' ) mixf(region)%f2d(ec440dryaer )%mf = metafields( -1 , 'ec440dryaer ', 'Surface Dry Aerosol Extinction@440nm' , 'm-1 ', '', 'surface_extinction_due_to_dry_aerosol' ) mixf(region)%f2d(abs440dryaer)%mf = metafields( -1 , 'abs440dryaer', 'Surface Dry Aerosol Absorption@440nm' , 'm-1 ', '', 'surface_absorption_due_to_dry_aerosol' ) mixf(region)%f2d(od870aer )%mf = metafields( -1 , 'od870aer ', 'AOD@870nm ' , '1 ', '', 'atmosphere_optical_thickness_due_to_ambient_aerosol' ) mixf(region)%f2d(abs870aer )%mf = metafields( -1 , 'abs870aer ', 'Absorption AOD@870nm ' , '1 ', '', 'atmosphere_absorption_optical_thickness_due_to_ambient_aerosol' ) mixf(region)%f2d(ec870dryaer )%mf = metafields( -1 , 'ec870dryaer ', 'Surface Dry Aerosol Extinction@870nm' , 'm-1 ', '', 'surface_extinction_due_to_dry_aerosol' ) mixf(region)%f2d(abs870dryaer)%mf = metafields( -1 , 'abs870dryaer', 'Surface Dry Aerosol Absorption@870nm' , 'm-1 ', '', 'surface_absorption_due_to_dry_aerosol' ) mixf(region)%f2d(od350aer )%mf = metafields( -1 , 'od350aer ', 'AOD@350nm ' , '1 ', '', 'atmosphere_optical_thickness_due_to_ambient_aerosol' ) mixf(region)%f2d(abs350aer )%mf = metafields( -1 , 'abs350aer ', 'Absorption AOD@350nm ' , '1 ', '', 'atmosphere_absorption_optical_thickness_due_to_ambient_aerosol' ) !!$ !!$ !!$ mixf(region)%f2d(52)%mf = metafields( 'od550soa ', 'AOD@550nm SOA ' , '1 ' ) !!$ mixf(region)%f2d(53)%mf = metafields( 'od550bb ', 'AOD@550nm BB ' , '1 ' ) !!$ mixf(region)%f2d(54)%mf = metafields( 'gf90aer ', 'GF @ 90 % RH ' , '1 ' ) mixf(region)%f2d( tr2d_1)%mf = metafields( iso4nus, 'SO4NUS ', 'concentration of tracer 01' , 'kg m-3', '', 'concentration_of_tracer_dry_aerosol_in_air' ) mixf(region)%f2d( tr2d_2)%mf = metafields( iso4ais, 'SO4AIS ', 'concentration of tracer 02' , 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air' ) mixf(region)%f2d( tr2d_3)%mf = metafields( iso4acs, 'SO4ACS ', 'concentration of tracer 03' , 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air' ) mixf(region)%f2d( tr2d_4)%mf = metafields( iso4cos, 'SO4COS ', 'concentration of tracer 04' , 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air' ) mixf(region)%f2d( tr2d_5)%mf = metafields( ibcais , 'BCAIS ', 'concentration of tracer 05' , 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air' ) mixf(region)%f2d( tr2d_6)%mf = metafields( ibcacs , 'BCACS ', 'concentration of tracer 06' , 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air' ) mixf(region)%f2d( tr2d_7)%mf = metafields( ibccos , 'BCCOS ', 'concentration of tracer 07' , 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air' ) mixf(region)%f2d( tr2d_8)%mf = metafields( ibcaii , 'BCAII ', 'concentration of tracer 08' , 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air' ) mixf(region)%f2d( tr2d_9)%mf = metafields( ipomais, 'POMAIS ', 'concentration of tracer 09' , 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air' ) mixf(region)%f2d(tr2d_10)%mf = metafields( ipomacs, 'POMACS ', 'concentration of tracer 10' , 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air' ) mixf(region)%f2d(tr2d_11)%mf = metafields( ipomcos, 'POMCOS ', 'concentration of tracer 11' , 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air' ) mixf(region)%f2d(tr2d_12)%mf = metafields( ipomaii, 'POMAII ', 'concentration of tracer 12' , 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air' ) mixf(region)%f2d(tr2d_13)%mf = metafields( issacs , 'SSACS ', 'concentration of tracer 13' , 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air' ) mixf(region)%f2d(tr2d_14)%mf = metafields( isscos , 'SSCOS ', 'concentration of tracer 14' , 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air' ) mixf(region)%f2d(tr2d_15)%mf = metafields( iduacs , 'DUACS ', 'concentration of tracer 15' , 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air' ) mixf(region)%f2d(tr2d_16)%mf = metafields( iducos , 'DUCOS ', 'concentration of tracer 16' , 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air' ) mixf(region)%f2d(tr2d_17)%mf = metafields( iduaci , 'DUACI ', 'concentration of tracer 17' , 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air' ) mixf(region)%f2d(tr2d_18)%mf = metafields( iducoi , 'DUCOI ', 'concentration of tracer 18' , 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air' ) mixf(region)%f2d(tr2d_19)%mf = metafields( isoanus, 'SOANUS ', 'concentration of tracer 19' , 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air' ) mixf(region)%f2d(tr2d_20)%mf = metafields( isoaais, 'SOAAIS ', 'concentration of tracer 20' , 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air' ) mixf(region)%f2d(tr2d_21)%mf = metafields( isoaacs, 'SOAACS ', 'concentration of tracer 21' , 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air' ) mixf(region)%f2d(tr2d_22)%mf = metafields( isoacos, 'SOACOS ', 'concentration of tracer 22' , 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air' ) mixf(region)%f2d(tr2d_23)%mf = metafields( isoaaii, 'SOAAII ', 'concentration of tracer 23' , 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air' ) !!$ mixf(region)%f2d(cc01)%mf = metafields( inus_n , 'conc_NUS', 'number concentration of mode 01' , 'm-3' ,'', 'number_concentration_of_ambient_aerosol_in_air' ) mixf(region)%f2d(cc02)%mf = metafields( iais_n , 'conc_AIS', 'number concentration of mode 02' , 'm-3' ,'', 'number_concentration_of_ambient_aerosol_in_air' ) mixf(region)%f2d(cc03)%mf = metafields( iacs_n , 'conc_ACS', 'number concentration of mode 03' , 'm-3' ,'', 'number_concentration_of_ambient_aerosol_in_air' ) mixf(region)%f2d(cc04)%mf = metafields( icos_n , 'conc_COS', 'number concentration of mode 04' , 'm-3' ,'', 'number_concentration_of_ambient_aerosol_in_air' ) mixf(region)%f2d(cc05)%mf = metafields( iaii_n , 'conc_AII', 'number concentration of mode 05' , 'm-3' ,'', 'number_concentration_of_ambient_aerosol_in_air' ) mixf(region)%f2d(cc06)%mf = metafields( iaci_n , 'conc_ACI', 'number concentration of mode 06' , 'm-3' ,'', 'number_concentration_of_ambient_aerosol_in_air' ) mixf(region)%f2d(cc07)%mf = metafields( icoi_n , 'conc_COI', 'number concentration of mode 07' , 'm-3' ,'', 'number_concentration_of_ambient_aerosol_in_air' ) mixf(region)%f2d(h2o1)%mf = metafields( -1 , 'aerh2o01 ', 'conc of aerosol water of mode 01' , 'kg m-3','', 'concentration_of_water_in_ambient_aerosol_in_air' ) mixf(region)%f2d(h2o2)%mf = metafields( -1 , 'aerh2o02 ', 'conc of aerosol water of mode 02' , 'kg m-3','', 'concentration_of_water_in_ambient_aerosol_in_air' ) mixf(region)%f2d(h2o3)%mf = metafields( -1 , 'aerh2o03 ', 'conc of aerosol water of mode 03' , 'kg m-3','', 'concentration_of_water_in_ambient_aerosol_in_air' ) mixf(region)%f2d(h2o4)%mf = metafields( -1 , 'aerh2o04 ', 'conc of aerosol water of mode 04' , 'kg m-3','', 'concentration_of_water_in_ambient_aerosol_in_air' ) mixf(region)%f2d(rw01)%mf = metafields( -1 , 'RWETnus ', 'rwet nus ' , 'm','', 'wet_radius' ) mixf(region)%f2d(rw02)%mf = metafields( -1 , 'RWETais ', 'rwet ais ' , 'm','', 'wet_radius' ) mixf(region)%f2d(rw03)%mf = metafields( -1 , 'RWETacs ', 'rwet acs ' , 'm','', 'wet_radius' ) mixf(region)%f2d(rw04)%mf = metafields( -1 , 'RWETcos ', 'rwet cos ' , 'm','', 'wet_radius' ) mixf(region)%f2d(rw05)%mf = metafields( -1 , 'RWETaii ', 'rwet aii ' , 'm','', 'wet_radius' ) mixf(region)%f2d(rw06)%mf = metafields( -1 , 'RWETaci ', 'rwet aci ' , 'm','', 'wet_radius' ) mixf(region)%f2d(rw07)%mf = metafields( -1 , 'RWETcoi ', 'rwet coi ' , 'm','', 'wet_radius' ) mixf(region)%f2d(rd01)%mf = metafields( -1 , 'RDRYnus ', 'rdry nus ' , 'm','', 'dry radius' ) mixf(region)%f2d(rd02)%mf = metafields( -1 , 'RDRYais ', 'rdry ais ' , 'm','', 'dry radius' ) mixf(region)%f2d(rd03)%mf = metafields( -1 , 'RDRYacs ', 'rdry acs ' , 'm','', 'dry radius' ) mixf(region)%f2d(rd04)%mf = metafields( -1 , 'RDRYcos ', 'rdry cos ' , 'm','', 'dry radius' ) mixf(region)%f2d(emiterp)%mf = metafields( -1 , 'emiterp ', 'emiterp ' , 'kgm-2s-1','', 'emission_of_terpene' ) mixf(region)%f2d(emiisop)%mf = metafields( -1 , 'emiisop ', 'emiisop ' , 'kgm-2s-1','', 'emission_of_isoprene' ) !!$ mixf(region)%f2d(rd05)%mf = metafields( -1 , 'RDRYaii ', 'rdry ' , 'm','', 'precipitation_flux' ) !!$ mixf(region)%f2d(rd06)%mf = metafields( -1 , 'RDRYaci ', 'rdry ' , 'm','', 'precipitation_flux' ) !!$ mixf(region)%f2d(rd07)%mf = metafields( -1 , 'RDRYcoi ', 'rdry ' , 'm','', 'precipitation_flux' ) !!$ !mixf(region)%f2d(31)%mf = metafields( ino3_a , 'mmrno3_a ', 'mmr of nitrate aerosol' , 'kg m-3', 'mass_fraction_of_nitrate_dry_aerosol_in_air' ) !!$ !mixf(region)%f2d(32)%mf = metafields( inh4 , 'mmrnh4 ', 'mmr of ammonium' , 'kg m-3', 'mass_fraction_of_ammonium_dry_aerosol_in_air' ) !!$ !mixf(region)%f2d(33)%mf = metafields( imsa , 'mmrmsa ', 'mmr of methane sulfonic acid' , 'kg m-3', 'mass_fraction_of_methane_sulfonic_acid_in_air' ) mixf(region)%f2d(sconcsoa )%mf = metafields( -1 , 'sconcsoa ', 'Surface Concentration SOA ' , 'kg m-3 ', '', 'mass_concentration_of_secondary_particulate_organic_matter_dry_aerosol_in_air' ) mixf(region)%f2d(loadsoa )%mf = metafields( -1 , 'loadsoa ', 'Load of SOA ' , 'kg m-2 ', '', 'atmosphere_mass_content_of_secondary_particulate_organic_matter_dry_aerosol' ) mixf(region)%f2d(drysoa )%mf = metafields( -1 , 'drysoa ', 'Dry Deposition of SOA ' , 'kg m-2 s-1', 'down', & 'tendency_of_atmosphere_mass_content_of_secondary_particulate_organic_matter_dry_aerosol_due_to_dry_deposition' ) mixf(region)%f2d(wetsoa )%mf = metafields( -1, 'wetsoa ', 'Wet Deposition of SOA ', 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_secondary_particulate_organic_matter_dry_aerosol_due_to_wet_deposition') mixf(region)%f2d(p_svoc2d )%mf = metafields( -1, 'p_svoc2D ', 'Column integral of SVOC production ', 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_secondary_particulate_organic_matter_dry_aerosol_due_to_wet_deposition') mixf(region)%f2d(p_elvoc2d )%mf = metafields( -1, 'p_elvoc2D ', 'Column integral of ELVOC production ', 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_secondary_particulate_organic_matter_dry_aerosol_due_to_wet_deposition') mixf(region)%f2d(p_gas_so42d)%mf = metafields( -1, 'prod_gas_so42D ', 'Column integral of gas SO4 production', 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_gas_phase_sulfate' ) mixf(region)%f2d(p_liq_so42d)%mf = metafields( -1, 'prod_liq_so42D ', 'Column integral of liq SO4 production', 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_liquid_phase_sulfate' ) ! set global dimensions of fields for netcdf definitions mixf (region)%nlon = imr mixf (region)%nlat = jmr mixf (region)%nlev = lmr !allocate space for lon, lat, lev information allocate(mixf (region)%lon(imr)) allocate(mixf (region)%lat(jmr)) allocate(mixf (region)%lev(lmr)) ! save the lat and lon data for use in output mixf (region)%lat=global_lli(region)%lat_deg mixf (region)%lon=global_lli(region)%lon_deg mixf (region)%lev=(/1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34/) ! intermediate structures for budgets allocate(drydepos(region)%f2dslast(i1:i2,j1:j2,8)) allocate(wetdepos(region)%f2dslast(i1:i2,j1:j2,8)) allocate(emission(region)%f2dslast(i1:i2,j1:j2,9))! pom, bc, su, so2, dms, ss, dust, terp, isop ! these here are the initial budgets (monthly): 0.0 drydepos(region)%f2dslast = 0.0 wetdepos(region)%f2dslast = 0.0 emission(region)%f2dslast = 0.0 endif ! newsrun ! reset counters and accumulators mixf (region)%acct = 0 do i = 1, ntracer_3d mixf(region)%f3d(i)%field = 0.0 end do do i = 1, ntracer_2d mixf(region)%f2d(i)%field = 0.0 end do ! ---------------- ! open NetCDF file (mixf) ! ---------------- call tau2date(itau,idater) ! time (in days since 2001-01-01 00:00) !Change to from beginning of a run (idatei) !time_reftime6=sdate_simulation call date2tau( time_reftime6, itauref ) select case (trim(gen_freq)) case ('hourly') write(idates, '(i4,3i2.2)') idater(1), idater(2), idater(3), idater(4)+1 !write(idates, '(i4,2i2.2)') idater(1), idater(2), idater(3) write(lidates, '(i4,"-",i2.2,"-",i2.2," ",i2.2,":30")') idater(1), idater(2), idater(3), idater(4) idater(5:6) = (/30,0/) case ('daily') write(idates, '(i4,2i2.2)') idater(1), idater(2), idater(3) write(lidates, '(i4,"-",i2.2,"-",i2.2," 12:00")') idater(1), idater(2), idater(3) ! set noon (recommendations) idater(4:6) = (/12,0,0/) case ('monthly') ! for monthly files, set time to middle of the month write(idates, '(i4,i2.2)') idater(1), idater(2) mlength=mlen(idater(2)) if ( mod(mlength,2) == 0 ) then idater(3:6) = (/mlength/2 + 1,0,0,0/) write(lidates, '(i4,"-",i2.2,"-",i2.2," 00:00")') idater(1), idater(2), idater(3) else idater(3:6) = (/(mlength+1)/2,12,0,0/) write(lidates, '(i4,"-",i2.2,"-",i2.2," 12:00")') idater(1), idater(2), idater(3) endif case default write (gol,'("Unknown General output frequency;")'); call goErr end select call date2tau( idater, itaucur ) !move the timestamp in middle of the average period if (trim(gen_freq)=='hourly') then ! if average period is hour move the timestamp 30 min back time_shift=30*60 else if (trim(gen_freq)=='daily') then ! if average period is hour move the timestamp 30 min back time_shift = 12*60*60 else if (trim(gen_freq)=='monthly') then ! assume 30 day months time_shift=15*60*60*24 end if reftime = (itaucur - itauref + time_shift) / 3600.00 !86400. !call date2tau( !WRITE(str_reftime,'("days since ",i4,"-",i2,"-",i2," ",i2,":",i2)') idatei(1),idatei(2),idatei(3),idatei(4),idatei(5) !WRITE(str_reftime,'("hours since ",i4,"-",i2,"-",i2," ",i2,":",i2)') idatei(1),idatei(2),idatei(3),idatei(4),idatei(5) str_reftime=time_units ! Changed file name convention to General Control 2015 mixf(region)%fname = trim(outdir)//'/'//& trim(f_dataid)//'_'//& trim(f_model) //'_'//& trim(gen_exper)//'_'//& trim(idates) //'_'//& trim(gen_freq) //'.nc' #ifdef MPI ! overwrite existing files (clobber), provide MPI stuff: call MDF_Create( mixf(region)%fname, MDF_NETCDF4, MDF_REPLACE, mixf(region)%funit, status, & mpi_comm=localComm, mpi_info=MPI_INFO_NULL ) if (status/=0) then write (gol,'("from creating NetCDF4 file for writing in parallel;")'); call goErr write (gol,'("MDF module not compiled with netcdf4_par support ?")'); call goErr TRACEBACK; status=1; return end if #else ! overwrite existing files (clobber) call MDF_Create( mixf(region)%fname, MDF_NETCDF4, MDF_REPLACE, mixf(region)%funit, status ) IF_NOTOK_RETURN(status=1) #endif if(okdebug) then write(gol,*) 'output_general_init: File ', trim(mixf(region)%fname), ' opened ' ; call goPr endif ! global attributes call MDF_Put_Att( mixf(region)%funit, MDF_GLOBAL, 'title', 'Model output for BACCHUS', status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit, MDF_GLOBAL, 'source', 'TM5-MP revision xxx in SVN repository, including SOA description', status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit, MDF_GLOBAL, 'institution', 'Royal Netherlands Meteorological Institute (KNMI), De Bilt, The Netherlands)', status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit, MDF_GLOBAL, 'contact' , 'Tommi Bergman: tommi.bergman@knmi.nl ; Twan van Noije; noije@knmi.nl', status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit, MDF_GLOBAL, 'project_id', 'BACCHUS', status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit, MDF_GLOBAL, 'conventions', 'CF-1.0 - HTAP', status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit, MDF_GLOBAL, 'date', lidates, status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit, MDF_GLOBAL, 'time_units', time_units, status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit, MDF_GLOBAL, 'references', 'http://tm5.sourceforge.net/', status ) IF_NOTOK_MDF(fid=mixf(region)%funit) ! define dimensions call MDF_Def_Dim( mixf(region)%funit, 'lon', imr, lon_dimid, status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Def_Dim( mixf(region)%funit, 'lat', jmr, lat_dimid, status ) IF_NOTOK_MDF(fid=mixf(region)%funit) !write(3000,*)lmr call MDF_Def_Dim( mixf(region)%funit, 'lev', 34, lev_dimid, status ) IF_NOTOK_MDF(fid=mixf(region)%funit) !Unlimited time causes a crash in the parallel NETCDF, when unlimited dimension is increased in the file !4.4.0 may correct this, but for now netcdf v4.4.0 on cca is not working call MDF_Def_Dim( mixf(region)%funit, 'time', 1, mixf(region)%timeid, status ) IF_NOTOK_MDF(fid=mixf(region)%funit) ! define dimension variables ! longitude call MDF_Def_Var( mixf(region)%funit, 'lon', MDF_DOUBLE, & (/lon_dimid/), mixf(region)%lonid, status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit,mixf(region)%lonid , 'units', 'degrees_east', status) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit,mixf(region)%lonid , 'axis', 'X', status) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit,mixf(region)%lonid , 'long_name', 'longitude', status) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit,mixf(region)%lonid , 'standard_name', 'longitude', status) IF_NOTOK_MDF(fid=mixf(region)%funit) ! Write out the longitudes call MDF_Put_Var( mixf(region)%funit, mixf(region)%lonid, mixf(region)%lon, status) IF_NOTOK_MDF(fid=mixf(region)%funit) ! define latitude variable call MDF_Def_Var( mixf(region)%funit, 'lat', MDF_DOUBLE, & (/lat_dimid/), mixf(region)%latid, status ) IF_NOTOK_MDF(fid=mixf(region)%funit) !write out the latitude to variable call MDF_Put_Var( mixf(region)%funit, mixf(region)%latid, mixf(region)%lat, status) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit,mixf(region)%latid , 'units', 'degrees_north', status) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit,mixf(region)%latid , 'axis', 'Y', status) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit,mixf(region)%latid , 'long_name', 'latitude', status) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit,mixf(region)%latid , 'standard_name', 'latitude', status) IF_NOTOK_MDF(fid=mixf(region)%funit) ! lev call MDF_Def_Var( mixf(region)%funit, 'lev' , MDF_DOUBLE, & (/lev_dimid/), mixf(region)%levid, status ) IF_NOTOK_MDF(fid=mixf(region)%funit) !call MDF_Put_Var( mixf(region)%funit, mixf(region)%levid, (/1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34/), status) call MDF_Put_Var( mixf(region)%funit, mixf(region)%levid, mixf(region)%lev, status) IF_NOTOK_MDF(fid=mixf(region)%funit) ! time call MDF_Def_Var( mixf(region)%funit, 'time', MDF_DOUBLE, & (/mixf(region)%timeid/), mixf(region)%time_varid, status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit,mixf(region)%time_varid , 'units', str_reftime, status) IF_NOTOK_MDF(fid=mixf(region)%funit) !call MDF_Put_Att( mixf(region)%funit,mixf(region)%timeid , 'long_name', 'time', status) !IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit,mixf(region)%time_varid , 'standard_name', 'time', status) IF_NOTOK_MDF(fid=mixf(region)%funit) ! define variables !do i = 1, ntrace!r_3d do i = 1, n_3d_vars!ntracer_3d !write(1234,*)i,trim(mixf(region)%f3d(i)%mf%vname) call MDF_Def_Var( mixf(region)%funit, trim(mixf(region)%f3d(i)%mf%vname), MDF_FLOAT, & (/mixf(region)%lonid, mixf(region)%latid, mixf(region)%levid, mixf(region)%timeid/), varid, status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Var_Par_Access( mixf(region)%funit, varid, access_mode, status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit, varid, 'long_name', trim(mixf(region)%f3d(i)%mf%lname), status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit, varid, 'standard_name', trim(mixf(region)%f3d(i)%mf%standard_name), status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit, varid, 'units', trim(mixf(region)%f3d(i)%mf%unit), status ) IF_NOTOK_MDF(fid=mixf(region)%funit) ! call MDF_Put_Att( mixf(region)%funit, varid, 'time', reftime, status ) ! IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit, varid, 'cell_methods', 'time: mean', status ) IF_NOTOK_MDF(fid=mixf(region)%funit) if( mixf(region)%f3d(i)%mf%positive /= '' ) then call MDF_Put_Att( mixf(region)%funit, varid, 'positive', trim(mixf(region)%f3d(i)%mf%positive), status ) IF_NOTOK_MDF(fid=mixf(region)%funit) end if mixf(region)%f3d(i)%varid = varid end do do i = 1, ntracer_2d ! write(1234,*)i,trim(mixf(region)%f2d(i)%mf%vname) call MDF_Def_Var( mixf(region)%funit, trim(mixf(region)%f2d(i)%mf%vname), MDF_FLOAT, & (/mixf(region)%lonid, mixf(region)%latid, mixf(region)%timeid/), varid, status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Var_Par_Access( mixf(region)%funit, varid, access_mode, status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit, varid, 'long_name', trim(mixf(region)%f2d(i)%mf%lname), status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit, varid, 'standard_name', trim(mixf(region)%f2d(i)%mf%standard_name), status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit, varid, 'units', trim(mixf(region)%f2d(i)%mf%unit), status ) IF_NOTOK_MDF(fid=mixf(region)%funit) ! call MDF_Put_Att( mixf(region)%funit, varid, 'time', reftime, status ) ! IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit, varid, 'cell_methods', 'time: mean', status ) IF_NOTOK_MDF(fid=mixf(region)%funit) if( mixf(region)%f2d(i)%mf%positive /= '' ) then call MDF_Put_Att( mixf(region)%funit, varid, 'positive', trim(mixf(region)%f2d(i)%mf%positive), status ) IF_NOTOK_MDF(fid=mixf(region)%funit) end if mixf(region)%f2d(i)%varid = varid end do call MDF_Def_Dim( mixf(region)%funit, 'kappa', 1, kappa_dimid, status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Def_Var( mixf(region)%funit, 'year', MDF_FLOAT, & (/kappa_dimid/),year_varid, status ) call MDF_Def_Var( mixf(region)%funit, 'month', MDF_FLOAT, & (/kappa_dimid/),month_varid, status ) call MDF_Def_Var( mixf(region)%funit, 'kSO4', MDF_FLOAT, & (/kappa_dimid/),kso4_varid, status ) !status = NF90_Def_Var( mixf(region)%funit, 'kSO4', MDF_FLOAT, & ! kso4_varid ) call MDF_Def_Var( mixf(region)%funit, 'kPOM', MDF_FLOAT, & (/kappa_dimid/), kpom_varid, status ) call MDF_Def_Var( mixf(region)%funit, 'kSOA', MDF_FLOAT, & (/kappa_dimid/), ksoa_varid, status ) call MDF_Def_Var( mixf(region)%funit, 'kBC', MDF_FLOAT, & (/kappa_dimid/), kbc_varid, status ) call MDF_Def_Var( mixf(region)%funit, 'kDU', MDF_FLOAT, & (/kappa_dimid/), kdu_varid, status ) call MDF_Def_Var( mixf(region)%funit, 'kSS', MDF_FLOAT, & (/kappa_dimid/), kss_varid, status ) call MDF_Def_Var( mixf(region)%funit, 'kNO3', MDF_FLOAT, & (/kappa_dimid/), kno3_varid, status ) call MDF_Def_Var( mixf(region)%funit, 'kNa2SO4', MDF_FLOAT, & (/kappa_dimid/), kna2so4_varid, status ) call MDF_Def_Var( mixf(region)%funit, 'kMSA', MDF_FLOAT, & (/kappa_dimid/), kmsa_varid, status ) call MDF_EndDef( mixf(region)%funit, status ) IF_NOTOK_MDF(fid=mixf(region)%funit) ! KAPPA values written only once ! needed for BACCHUS intercomparison... call MDF_Put_Var( mixf(region)%funit, year_varid, (/idater(1)/), status) call MDF_Put_Var( mixf(region)%funit, month_varid, (/idater(2)/), status) call MDF_Put_Var( mixf(region)%funit, kso4_varid, (/kap_su/), status) !status= NF90_Put_Var( mixf(region)%funit, kso4_varid, kap_su) !IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Var( mixf(region)%funit, kpom_varid, (/kap_pom/), status) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Var( mixf(region)%funit, kbc_varid, (/kap_bc/), status) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Var( mixf(region)%funit, ksoa_varid, (/kap_soa/), status) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Var( mixf(region)%funit, kss_varid, (/kap_ss/), status) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Var( mixf(region)%funit, kdu_varid, (/kap_du/), status) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Var( mixf(region)%funit, kno3_varid, (/kap_no3/), status) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Var( mixf(region)%funit, kna2so4_varid, (/kap_na2so4/), status) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Var( mixf(region)%funit, kmsa_varid, (/kap_msa/), status) IF_NOTOK_MDF(fid=mixf(region)%funit) end do regionloop ! nregions call goLabel() ; status = 0 end subroutine output_general_init !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: OUTPUT_GENERAL_DONE ! ! !DESCRIPTION: Free parameters. !\\ !\\ ! !INTERFACE: ! subroutine output_general_done(status, iregion) ! ! !USES: use chem_param, only : ntracet,ntrace ! ! !INPUT PARAMETERS: ! !logical, intent(in) :: stat_output ! true if stations output is requested integer, intent(in), optional :: iregion ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 29 Nov 2010 - Achim Strunk - ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/output_general_done' integer :: i, region ! --- begin ------------------------------------- call goLabel(rname) deallocate( wdep_out ) regionloop: do region = 1, nregions ! if region given, cycle if other region! if (present(iregion)) then if(iregion /= region) cycle regionloop endif do i=1,ntracer_3d deallocate( mixf(region)%f3d(i)%field ) end do do i=1,ntracer_2d deallocate( mixf(region)%f2d(i)%field ) end do deallocate( gas_output) deallocate( mixf (region)%f3d ) deallocate( mixf (region)%f2d ) deallocate( drydepos(region)%f2dslast ) deallocate( wetdepos(region)%f2dslast ) deallocate( emission(region)%f2dslast ) end do regionloop if (nCCNdiag) then deallocate(ind_CCN) ! deallocate(N_rad) ! deallocate(HG_rad) end if call goLabel() ; status = 0 end subroutine output_general_done !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: output_general_write ! ! !DESCRIPTION: This routines builds the average by dividing accumulated ! data by stack counter. The results are written to file. !\\ !\\ ! !INTERFACE: ! subroutine output_general_write(region, status) ! ! !USES: ! use chem_param, only : ntracet,ntrace use dims,only:itau use datetime, only : tau2date, date2tau ! !INPUT PARAMETERS: ! integer, intent(in) :: region !logical, intent(in) :: stat_output ! true if stations output is requested ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 29 Nov 2010 - Achim Strunk - ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/output_general_write' integer :: i, imr, jmr, lmr integer :: i1, i2, j1, j2, ilev integer :: istat ! Time definitions for General real :: reftime integer :: itauref, time_shift ! --- begin ------------------------------------- call goLabel(rname) ! grid size call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) imr = i2-i1+1 jmr = j2-j1+1 lmr = levi%nlev ! this here is already accumulated over the time period (day) call collect_budgets( region, status ) IF_NOTOK_RETURN(status=1) ! --------------------- ! divide by accumulator ! --------------------- do i = 1, n_3d_vars!ntracer_3d mixf(region)%f3d(i)%field = mixf(region)%f3d(i)%field / real( mixf(region)%acct ) end do do i = 1, ntracer_2d mixf(region)%f2d(i)%field = mixf(region)%f2d(i)%field / real( mixf(region)%acct ) end do select case (trim(gen_freq)) case ('hourly') write(gol,'("---> GENERAL diagnostics: write file for previous hour")'); call goPr case ('daily') write(gol,'("---> GENERAL diagnostics: write file for previous day")'); call goPr case ('monthly') write(gol,'("---> GENERAL diagnostics: write file for previous month")'); call goPr end select ! ------------- ! write to file ! ------------- ! Ncfile needs a timestep ! define the reference time ! time (for now in days since 2001-01-01 00:00) call date2tau( time_reftime6, itauref ) ! calculate reftime as fractional days from itauref, hence division by 86400 (24h*3600s) !call date2tau( idater, itaucur ) if (trim(gen_freq)=='hourly') then ! if average period is hour move the timestamp 30 min forward because current itau is ! at the beginning of the timestep time_shift=30.0*60.0 else if (trim(gen_freq)=='daily') then ! if average period is hour move the timestamp 30 min back time_shift = 12*60*60 else if (trim(gen_freq)=='monthly') then ! assume 30 day months time_shift=15*60*60*24 end if !reftime = (itau - itauref-time_shift) / 3600 !86400. in hours instead of days reftime = (itau - itauref+time_shift) / 3600 !86400. in hours instead of days !write (1111,*)reftime,itau,itauref,itau-itauref, itau-itauref-time_shift,time_shift reftime =reftime-0.5 !reftime = (itau - itauref) / 86400. !Add time stamp to current !Currently only allows 1 step per file, needs to be extended do i=1,ntracer_2d call MDF_Put_Var( mixf(region)%funit, mixf(region)%f2d(i)%varid, mixf(region)%f2d(i)%field(i1:i2,j1:j2), status, start=(/i1,j1,1/), count=(/imr,jmr,1/) ) IF_NOTOK_MDF(fid=mixf(region)%funit) end do do i=1,n_3d_vars!ntracer_3d !write(1111,*) mixf(region)%f3d(i)%mf%vname,imr,jmr,lmr, mixf(region)%nlev,mixf(region)%nlat,mixf(region)%nlon call MDF_Put_Var( mixf(region)%funit, mixf(region)%f3d(i)%varid, mixf(region)%f3d(i)%field(i1:i2,j1:j2,1:lmr), status, start=(/i1,j1,1,1/), count=(/imr,jmr,lmr,1/) ) IF_NOTOK_MDF(fid=mixf(region)%funit) end do call MDF_Var_Par_Access( mixf(region)%funit, mixf(region)%time_varid, MDF_INDEPENDENT, status ) IF_NOTOK_MDF(fid=mixf(region)%funit) ! Write current reftime call MDF_Put_Var( mixf(region)%funit, mixf(region)%time_varid,(/reftime/), status, start=(/1/),count=(/1/)) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Close( mixf(region)%funit, status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call goLabel() ; status = 0 end subroutine output_general_write !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: OUTPUT_GENERAL_STEP ! ! !DESCRIPTION: This is the accumulation routine. It is called following ! the user specification general.dhour in rc-file. It ! evaluates the various diagnostics and does summing. !\\ !\\ ! !INTERFACE: ! subroutine output_general_step( region, dhour, status ) ! ! !USES: ! use GO , only : TDate, NewDate, rTotal, operator(-), GO_Timer_Def, GO_Timer_End, GO_Timer_Start use Grid , only : FPressure use binas, only : rgas, rol,pi use datetime, only : tau2date use MeteoData, only : sp_dat, lsp_dat, cp_dat, m_dat use MeteoData, only : temper_dat, humid_dat, gph_dat use global_data, only : mass_dat, region_dat use tracer_data, only : chem_dat use dims, only : itaur use chemistry, only : d_nucle, m_nucle_su, m_nucle_soa, growth1_2, coag_sink_nuc, cond1_soa, cond1_su use chemistry, only : cond2_soa, cond3_soa, cond4_soa, cond5_soa use chem_param, only : iso4acs, iso4cos, iduacs, iso4ais, ibccos, ibcaii, xmair use chem_param, only : iso4nus, isscos, ino3_a, issacs, iducos, iduaci, nmod use chem_param, only : iducoi, ibcacs, ipomcos, ipomaii, ibcais, ipomacs, ipomais,isoanus use chem_param, only : isoaais, isoaacs, isoacos, isoaaii use chem_param, only : imsa, inh4 use chem_param, only : inus_n,iacs_n,icos_n,iais_n,iaii_n,iaci_n,icoi_n use chem_param, only : ntracet,ntrace,xms,xmso4,xmc,xmbc,xmpom use chem_param, only : mode_end_pom,mode_end_so4,mode_end_ss,mode_end_bc,mode_end_dust,mode_end_soa use mo_aero_m7, only : nsol,nmod,dnacl,doc,dh2so4,dbc,ddust use optics, only : optics_aop_get use m7_data, only : h2o_mode, rw_mode, rwd_mode use ebischeme, only : diag_prod use mo_aero_m7, only : cmedr2mmedr,sigma ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region integer, intent(in) :: dhour ! this is general.dhour [hours] !logical, intent(in) :: stat_output ! true if stations output is requested ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 29 Nov 2010 - Achim Strunk - Creation ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/output_general_step' ! MPI arrays to gather fields. real, dimension(:,:,:,:), pointer :: rm, rmc real, dimension(:,:,:), allocatable, target :: mg real, dimension(:), pointer :: dxyp real, dimension(:,:,:), allocatable :: pres3d integer :: i, j, k, imr, jmr, lmr, lwl, dtime, iSat integer :: i1, i2, j1, j2 integer :: ie, je ! indices for subdomain extended with halo cells integer, parameter :: l_halo=1 logical :: polar integer :: istat, imode real :: dens, load_tmp integer :: rwetcounter,h2ocounter,rdrycounter Real, Dimension(:,:,:), Allocatable :: aop_output_ext ! extinctions Real, Dimension(:,:), Allocatable :: aop_output_a ! single scattering albedo Real, Dimension(:,:), Allocatable :: aop_output_g ! assymetry factor Real, Dimension(:,:,:), Allocatable :: aop_output_add ! additional parameters real, dimension(:,:,:,:,:), allocatable :: opt_ext real, dimension(:,:,:,:), allocatable :: opt_a real, dimension(:,:,:,:), allocatable :: opt_g real, dimension(:,:,:,:,:), allocatable :: opt_add real, dimension(:,:,:), allocatable :: volume real, dimension(:,:), allocatable :: temp2d real, dimension(:,:,:,:), allocatable :: CCN real, dimension(:,:,:,:), allocatable :: dry_rad real, dimension(:,:,:,:), allocatable :: Kappa_r real, dimension(:,:,:,:), allocatable :: N_r integer :: ncontr, lvec, lct, l, indoffset, nwl,nn,irw,irwd integer :: nadd, nadd_max, indoffset_stat, indabs integer :: iadd real :: z,sizelimit,hr2,rad real :: no3fac, wght, dwght, tabs real,dimension(7) :: modfrac,lnsigma integer,dimension(6) :: idate_f type(TDate) :: t, t0 real :: time real :: perm3_to_percm3,kg_to_ug integer :: d_ind ! --- begin ------------------------------- call goLabel(rname) ! grid size call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2, hasNorthPole=polar ) imr = i2-i1+1 jmr = j2-j1+1 lmr = levi%nlev ! link the datastructures for transported and non transported variables and ! Gridbox area nullify(rmc) rmc => chem_dat(region)%rm dxyp => region_dat(region)%dxyp rm => mass_dat(region)%rm ! get helping pressure field in 3d ie=i2; je=j2 allocate( pres3d(i1:ie,j1:je,lmr) ) ! fill mid level pressure call FPressure( levi, sp_dat(region)%data(i1:ie,j1:je,1), pres3d, status ) IF_NOTOK_RETURN(status=1) allocate( CCN(nSat,i1:ie,j1:je,lmr)) allocate( dry_rad(nsol,i1:ie,j1:je,lmr)) allocate( Kappa_r(3,i1:ie,j1:je,lmr)) allocate( N_r(3,i1:ie,j1:je,lmr)) call CCN_Diag(CCN, rm(i1:ie,j1:je,1:lmr,1:ntracet), dry_rad, temper_dat(region)%data(i1:i2,j1:j2,1:lmr)) call Kappa_diag(CCN, rm(i1:ie,j1:je,1:lmr,1:ntracet), dry_rad, temper_dat(region)%data(i1:i2,j1:j2,1:lmr),N_r,Kappa_r) ! dtime is seconds per step ! dtime will be taken as weight for summing up ! (makes it fit for varying step lengths) dtime = dhour * 3600. ! sum up for later averaging mixf (region)%acct = mixf (region)%acct + dtime ! ---------------------- ! 3D meteo fields ! ---------------------- !temp=1 !hus=2 !airmass=3 !pressure=4 ! temperature mixf(region)%f3d(temp)%field = mixf(region)%f3d(temp)%field + dtime * temper_dat(region)%data(i1:i2,j1:j2,:) ! specific humidity mixf(region)%f3d(hus)%field = mixf(region)%f3d(hus)%field + dtime * humid_dat(region)%data(i1:i2,j1:j2,:) ! air mass (kg -> kg/m2) do j = j1, j2 mixf(region)%f3d(airmass)%field(:,j,:) = mixf(region)%f3d(airmass)%field(:,j,:) + & dtime * m_dat(region)%data(i1:i2,j,:) / dxyp(j) end do ! pressure (already filled above) mixf(region)%f3d(pressure)%field = mixf(region)%f3d(pressure)%field + dtime * pres3d do k=1,lmr do j=j1,j2 do i=i1,i2 dens = pres3d(i,j,k) / temper_dat(region)%data(i,j,k) * xmair * 1.E-3 / (m_dat(region)%data(i,j,k) * Rgas) !!$ mixf(region)%f3d(p_gas_so4)%field(i,j,k)= mixf(region)%f3d(p_gas_so4)%field(i,j,k)+diag_prod(region)%prod(i,j,k,1) / dxyp(j) !!$ mixf(region)%f3d(p_liq_so4)%field(i,j,k)= mixf(region)%f3d(p_liq_so4)%field(i,j,k)+diag_prod(region)%prod(i,j,k,2) / dxyp(j) mixf(region)%f2d(p_gas_so42d)%field(i,j)= mixf(region)%f2d(p_gas_so42d)%field(i,j)+diag_prod(region)%prod(i,j,k,1) / dxyp(j) mixf(region)%f2d(p_liq_so42d)%field(i,j)= mixf(region)%f2d(p_liq_so42d)%field(i,j)+diag_prod(region)%prod(i,j,k,2) / dxyp(j) !!$ !!$ mixf(region)%f3d(p_elvoc)%field(i,j,k)= mixf(region)%f3d(p_elvoc)%field(i,j,k)+diag_prod(region)%prod(i,j,k,3) / dxyp(j)! divide by 2 to calculate mean over 2 chemistry timesteps in one model timestep !!$ mixf(region)%f3d(p_svoc)%field(i,j,k)= mixf(region)%f3d(p_svoc)%field(i,j,k)+diag_prod(region)%prod(i,j,k,4) / dxyp(j) !!$ !!$ mixf(region)%f3d(p_elohterp)%field(i,j,k)= mixf(region)%f3d(p_elohterp)%field(i,j,k)+diag_prod(region)%prod(i,j,k,5) / dxyp(j) !!$ mixf(region)%f3d(p_elo3terp)%field(i,j,k)= mixf(region)%f3d(p_elo3terp)%field(i,j,k)+diag_prod(region)%prod(i,j,k,6) / dxyp(j) !!$ mixf(region)%f3d(p_elohisop)%field(i,j,k)= mixf(region)%f3d(p_elohisop)%field(i,j,k)+diag_prod(region)%prod(i,j,k,7) / dxyp(j) !!$ mixf(region)%f3d(p_elo3isop)%field(i,j,k)= mixf(region)%f3d(p_elo3isop)%field(i,j,k)+diag_prod(region)%prod(i,j,k,8) / dxyp(j) !!$ mixf(region)%f3d(p_svohterp)%field(i,j,k)= mixf(region)%f3d(p_svohterp)%field(i,j,k)+diag_prod(region)%prod(i,j,k,9) / dxyp(j) !!$ mixf(region)%f3d(p_svo3terp)%field(i,j,k)= mixf(region)%f3d(p_svo3terp)%field(i,j,k)+diag_prod(region)%prod(i,j,k,10) / dxyp(j) !!$ mixf(region)%f3d(p_svohisop)%field(i,j,k)= mixf(region)%f3d(p_svohisop)%field(i,j,k)+diag_prod(region)%prod(i,j,k,11) / dxyp(j) !!$ mixf(region)%f3d(p_svo3isop)%field(i,j,k)= mixf(region)%f3d(p_svo3isop)%field(i,j,k)+diag_prod(region)%prod(i,j,k,12) / dxyp(j) mixf(region)%f2d(p_elvoc2d)%field(i,j)= mixf(region)%f2d(p_elvoc2d)%field(i,j)+diag_prod(region)%prod(i,j,k,3) / dxyp(j) mixf(region)%f2d(p_svoc2d)%field(i,j)= mixf(region)%f2d(p_svoc2d)%field(i,j)+diag_prod(region)%prod(i,j,k,4) / dxyp(j) !!$ mixf(region)%f3d(d_nuc)%field(i,j,k)= mixf(region)%f3d(d_nuc)%field(i,j,k)+d_nucle(i,j,k)*dtime !value is saved #cm-3s-1 !!$ mixf(region)%f3d(m_nuc_su)%field(i,j,k)= mixf(region)%f3d(m_nuc_su)%field(i,j,k)+m_nucle_su(i,j,k)*dtime !value is saved #molec cm-3s-1 !!$ mixf(region)%f3d(m_nuc_soa)%field(i,j,k)= mixf(region)%f3d(m_nuc_soa)%field(i,j,k)+m_nucle_soa(i,j,k)*dtime !value is saved #µg m-3s-1 !!$ mixf(region)%f3d(gr1_2)%field(i,j,k)= mixf(region)%f3d(gr1_2)%field(i,j,k)+growth1_2(i,j,k)*dtime !value is saved #cm-3s-1 !!$ mixf(region)%f3d(coag1)%field(i,j,k)= mixf(region)%f3d(coag1)%field(i,j,k)+coag_sink_nuc(i,j,k)*dtime !value is saved #cm-3s-1 !!$ mixf(region)%f3d(co1_soa)%field(i,j,k)= mixf(region)%f3d(co1_soa)%field(i,j,k)+cond1_soa(i,j,k)*dtime !value is saved # unit? !!$ mixf(region)%f3d(co2_soa)%field(i,j,k)= mixf(region)%f3d(co2_soa)%field(i,j,k)+cond2_soa(i,j,k)*dtime !value is saved # unit? !!$ mixf(region)%f3d(co3_soa)%field(i,j,k)= mixf(region)%f3d(co3_soa)%field(i,j,k)+cond3_soa(i,j,k)*dtime !value is saved # unit? !!$ mixf(region)%f3d(co4_soa)%field(i,j,k)= mixf(region)%f3d(co4_soa)%field(i,j,k)+cond4_soa(i,j,k)*dtime !value is saved # unit? !!$ mixf(region)%f3d(co5_soa)%field(i,j,k)= mixf(region)%f3d(co5_soa)%field(i,j,k)+cond5_soa(i,j,k)*dtime !value is saved # unit? !!$ mixf(region)%f3d(co1_su)%field(i,j,k)= mixf(region)%f3d(co1_su)%field(i,j,k)+cond1_su(i,j,k)*dtime !value is saved # unit? !BACCHUS pm1 vars lnsigma=log(sigma) hr2= 0.5*sqrt(2.0) sizelimit=1.0*5.e-7! pm1(1.0 um in diam) in radius (m) Do imode = 1, nmod ! Calculate the median radius of the volume weighted distribution. rad = rw_mode(region,imode)%d3(i,j,k) * cmedr2mmedr(imode) !! if( rad == 0.0 ) then ! changed to a precision depending value if( rad < 100*TINY(1.0) ) then modfrac(imode) = 1.0 ! Should not matter, because if the radius is zero, ! there are no aerosols. But in principle, it would ! mean that all aerosols are infinitely small, so ! they would fit in any PM-class. else z = ( log(sizelimit) - log(rad) ) / lnsigma(imode) modfrac(imode) = 0.5 + 0.5 * ERF(z*hr2) end if end do perm3_to_percm3=1.0e-6 kg_to_ug=1.0e9 mixf(region)%f3d(mmf1)%field(i,j,k) = mixf(region)%f3d(mmf1)%field(i,j,k) + dtime * & modfrac(1) mixf(region)%f3d(mmf2)%field(i,j,k) = mixf(region)%f3d(mmf2)%field(i,j,k) + dtime * & modfrac(2) mixf(region)%f3d(mmf3)%field(i,j,k) = mixf(region)%f3d(mmf3)%field(i,j,k) + dtime * & modfrac(3) mixf(region)%f3d(mmf4)%field(i,j,k) = mixf(region)%f3d(mmf4)%field(i,j,k) + dtime * & modfrac(4) mixf(region)%f3d(so4pm1)%field(i,j,k) = mixf(region)%f3d(so4pm1)%field(i,j,k) + dtime * & dens * (rm(i,j,k,iso4nus)*modfrac(1)+rm(i,j,k,iso4ais)*modfrac(2)+& rm(i,j,k,iso4acs)*modfrac(3)+rm(i,j,k,iso4cos)*modfrac(4))*(xms/xmso4)*kg_to_ug mixf(region)%f3d(bcpm1)%field(i,j,k) = mixf(region)%f3d(bcpm1)%field(i,j,k) + dtime * & dens * (rm(i,j,k,ibcais)*modfrac(2)+& rm(i,j,k,ibcacs)*modfrac(3)+rm(i,j,k,ibccos)*modfrac(4)+rm(i,j,k,ibcaii)*modfrac(5))*kg_to_ug*(xmc/xmbc) mixf(region)%f3d(orgpm1)%field(i,j,k) = mixf(region)%f3d(orgpm1)%field(i,j,k) + dtime * & dens * (rm(i,j,k,isoanus)*modfrac(1)+(rm(i,j,k,isoaais)+rm(i,j,k,ipomais))*modfrac(2)+& (rm(i,j,k,isoaacs)+rm(i,j,k,ipomacs))*modfrac(3)+(rm(i,j,k,isoacos)+rm(i,j,k,ipomcos))*modfrac(4)+& (rm(i,j,k,isoaaii)+rm(i,j,k,ipomaii))*modfrac(5))*kg_to_ug*(xmc/xmpom) mixf(region)%f3d(soapm1)%field(i,j,k) = mixf(region)%f3d(soapm1)%field(i,j,k) + dtime * & dens * (rm(i,j,k,isoanus)*modfrac(1)+(rm(i,j,k,isoaais))*modfrac(2)+& (rm(i,j,k,isoaacs))*modfrac(3)+(rm(i,j,k,isoacos))*modfrac(4)+& (rm(i,j,k,isoaaii))*modfrac(5))*kg_to_ug*(xmc/xmpom) mixf(region)%f3d(sspm1)%field(i,j,k) = mixf(region)%f3d(sspm1)%field(i,j,k) + dtime * & dens * ( rm(i,j,k,issacs)*modfrac(3)+rm(i,j,k,isscos)*modfrac(4))*kg_to_ug mixf(region)%f3d(dupm1)%field(i,j,k) = mixf(region)%f3d(dupm1)%field(i,j,k) + dtime * & dens * (rm(i,j,k,iduacs)*modfrac(3)+& rm(i,j,k,iducos)*modfrac(4)+rm(i,j,k,iduaci)*modfrac(6)+rm(i,j,k,iducoi)*modfrac(7))*kg_to_ug mixf(region)%f3d(nh4pm1)%field(i,j,k) = mixf(region)%f3d(nh4pm1)%field(i,j,k) + dtime * & dens * (rm(i,j,k,inh4)*modfrac(3))*kg_to_ug mixf(region)%f3d(no3pm1)%field(i,j,k) = mixf(region)%f3d(no3pm1)%field(i,j,k) + dtime * & dens * (rm(i,j,k,ino3_a)*modfrac(3))*kg_to_ug mixf(region)%f3d(so4)%field(i,j,k) = mixf(region)%f3d(so4)%field(i,j,k) + dtime * & dens * (rm(i,j,k,iso4nus)+rm(i,j,k,iso4ais)+& rm(i,j,k,iso4acs)+rm(i,j,k,iso4cos))*(xms/xmso4)*kg_to_ug mixf(region)%f3d(bc)%field(i,j,k) = mixf(region)%f3d(bc)%field(i,j,k) + dtime * & dens * (rm(i,j,k,ibcais)+& rm(i,j,k,ibcacs)+rm(i,j,k,ibccos)+rm(i,j,k,ibcaii))*(xmc/xmbc)*kg_to_ug mixf(region)%f3d(org)%field(i,j,k) = mixf(region)%f3d(org)%field(i,j,k) + dtime * & dens * (rm(i,j,k,isoanus)+(rm(i,j,k,isoaais)+rm(i,j,k,ipomais))+& (rm(i,j,k,isoaacs)+rm(i,j,k,ipomacs))+(rm(i,j,k,isoacos)+rm(i,j,k,ipomcos))+& (rm(i,j,k,isoaaii)+rm(i,j,k,ipomaii)))*(xmc/xmpom)*kg_to_ug mixf(region)%f3d(ss)%field(i,j,k) = mixf(region)%f3d(ss)%field(i,j,k) + dtime * & dens * ( rm(i,j,k,issacs)+rm(i,j,k,isscos))*kg_to_ug mixf(region)%f3d(du)%field(i,j,k) = mixf(region)%f3d(du)%field(i,j,k) + dtime * & dens * (rm(i,j,k,iduacs)+& rm(i,j,k,iducos)+rm(i,j,k,iduaci)+rm(i,j,k,iducoi))*kg_to_ug mixf(region)%f3d(nh4)%field(i,j,k) = mixf(region)%f3d(nh4)%field(i,j,k) + dtime * & dens * (rm(i,j,k,inh4))*kg_to_ug mixf(region)%f3d(no3)%field(i,j,k) = mixf(region)%f3d(no3)%field(i,j,k) + dtime * & dens * (rm(i,j,k,ino3_a))*kg_to_ug sizelimit=3*5.e-9! pm1(1.0 um in diam) in radius (m) Do imode = 1, nmod ! Calculate the median radius of the volume weighted distribution. rad = rw_mode(region,imode)%d3(i,j,k) if( rad < 100*TINY(1.0) ) then modfrac(imode) = 1.0 ! Should not matter, because if the radius is zero, ! there are no aerosols. But in principle, it would ! mean that all aerosols are infinitely small, so ! they would fit in any PM-class. else z = ( log(sizelimit) - log(rad) ) / lnsigma(imode) ! now calculating larger than sizelimit modfrac(imode) = 1-(0.5 + 0.5 * ERF(z*hr2)) end if end do mixf(region)%f3d(n3)%field(i,j,k) = mixf(region)%f3d(n3)%field(i,j,k) + dtime * & dens *( rm(i,j,k,inus_n)*modfrac(1)+ rm(i,j,k,iais_n)*modfrac(2)+& rm(i,j,k,iacs_n)*modfrac(3)+ rm(i,j,k,icos_n)*modfrac(4)+& rm(i,j,k,iaii_n)*modfrac(5)+ rm(i,j,k,iaci_n)*modfrac(6)+& rm(i,j,k,icoi_n)*modfrac(7))*perm3_to_percm3 !!$ !50nm !!$ do d_ind=1,3 !!$ sizelimit=N_diam(d_ind)*5.e-9! pm1(1.0 um in diam) in radius (m) !!$ Do imode = 1, nmod !!$ ! Calculate the median radius of the volume weighted distribution. !!$ rad = rw_mode(region,imode)%d3(i,j,k) !!$ if( rad < 100*TINY(1.0) ) then !!$ modfrac(imode) = 1.0 ! Should not matter, because if the radius is zero, !!$ ! there are no aerosols. But in principle, it would !!$ ! mean that all aerosols are infinitely small, so !!$ ! they would fit in any PM-class. !!$ else !!$ z = ( log(sizelimit) - log(rad) ) / lnsigma(imode) !!$ !!$ ! now calculating larger than sizelimit !!$ modfrac(imode) = 1-(0.5 + 0.5 * ERF(z*hr2)) !!$ end if !!$ end do !!$ mixf(region)%f3d(ind_N(d_ind))%field(i,j,k) = mixf(region)%f3d(n3)%field(i,j,k) + dtime * & !!$ dens *( rm(i,j,k,inus_n)*modfrac(1)+ rm(i,j,k,iais_n)*modfrac(2)+& !!$ rm(i,j,k,iacs_n)*modfrac(3)+ rm(i,j,k,icos_n)*modfrac(4)+& !!$ rm(i,j,k,iaii_n)*modfrac(5)+ rm(i,j,k,iaci_n)*modfrac(6)+& !!$ rm(i,j,k,icoi_n)*modfrac(7))*perm3_to_percm3 !!$ end do ! CCN under given supersaturation IF (nCCNdiag) THEN DO iSat=1,nSat mixf(region)%f3d(ind_CCN(iSat))%field(i,j,k)= mixf(region)%f3d(ind_CCN(iSat))%field(i,j,k) + dtime * & dens * 1.e-06 * CCN(iSat,i,j,k) ! cm-3 ENDDO DO iSat=1,3 mixf(region)%f3d(ind_HG(iSat))%field(i,j,k)= mixf(region)%f3d(ind_HG(iSat))%field(i,j,k) + dtime * & Kappa_r(iSat,i,j,k) ! cm-3 mixf(region)%f3d(ind_N(iSat))%field(i,j,k)= mixf(region)%f3d(ind_N(iSat))%field(i,j,k) + dtime * & dens * 1.e-06 * N_r(iSat,i,j,k) ! cm-3 ENDDO ENDIF mixf(region)%f3d(ccn02)%field(i,j,k) = mixf(region)%f3d(ccn02)%field(i,j,k) + dtime * & dens * 1.e-06 * CCN(5,i,j,k) ! cm-3 CCN0.2 mixf(region)%f3d(gph)%field(i,j,k) = mixf(region)%f3d(gph)%field(i,j,k) + gph_dat(region)%data(i,j,k)*dtime ! rwet, rdry, h2o data holder indices go from 1-4/7 ! so use counters, could also add modeindex in mixf rwetcounter=1 rdrycounter=1 h2ocounter=1 do nn=1,n_3d_vars if (mixf(region)%f3d(nn)%mf%itm5>0)then ! transported variables in rm if (mixf(region)%f3d(nn)%mf%itm5<=ntracet) then mixf(region)%f3d(nn)%field(i,j,k) = mixf(region)%f3d(nn)%field(i,j,k) + dtime * & dens * rm(i,j,k,mixf(region)%f3d(nn)%mf%itm5) else !non transported variables in rmc mixf(region)%f3d(nn)%field(i,j,k) = mixf(region)%f3d(nn)%field(i,j,k) + dtime * & dens * rmc(i,j,k,mixf(region)%f3d(nn)%mf%itm5) end if else ! not a tracer variable ! counter to go through 7 modes when outputing rwet ! requires that rwet variables are in right order ! Wet radius of particles if ( mixf(region)%f3d(nn)%mf%vname(1:4)=='RWET')then mixf(region)%f3d(nn)%field(i,j,k) = mixf(region)%f3d(nn)%field(i,j,k) + dtime * & rw_mode(region,rwetcounter)%d3(i,j,k) rwetcounter=rwetcounter+1 end if ! aerosol water if ( mixf(region)%f3d(nn)%mf%vname(1:6)=='aerh2o')then mixf(region)%f3d(nn)%field(i,j,k) = mixf(region)%f3d(nn)%field(i,j,k) + dtime * & dens*h2o_mode(region,h2ocounter)%d3(i,j,k) h2ocounter=h2ocounter+1 end if !dry radius of soluble modes if ( mixf(region)%f3d(nn)%mf%vname(1:4)=='RDRY')then mixf(region)%f3d(nn)%field(i,j,k) = mixf(region)%f3d(nn)%field(i,j,k) + dtime * & rwd_mode(region,rdrycounter)%d3(i,j,k) rdrycounter=rdrycounter+1 end if end if end do end do end do end do ! set fields summed in chemistry to 0 after writing d_nucle = 0. m_nucle_su = 0. m_nucle_soa = 0. diag_prod(region)%prod(i1:i2,j1:j2,:,1)=0.0 diag_prod(region)%prod(i1:i2,j1:j2,:,2)=0.0 diag_prod(region)%prod(i1:i2,j1:j2,:,3)=0.0 diag_prod(region)%prod(i1:i2,j1:j2,:,4)=0.0 diag_prod(region)%prod(i1:i2,j1:j2,:,5)=0.0 diag_prod(region)%prod(i1:i2,j1:j2,:,6)=0.0 diag_prod(region)%prod(i1:i2,j1:j2,:,7)=0.0 diag_prod(region)%prod(i1:i2,j1:j2,:,8)=0.0 diag_prod(region)%prod(i1:i2,j1:j2,:,9)=0.0 diag_prod(region)%prod(i1:i2,j1:j2,:,10)=0.0 diag_prod(region)%prod(i1:i2,j1:j2,:,11)=0.0 diag_prod(region)%prod(i1:i2,j1:j2,:,12)=0.0 coag_sink_nuc = 0. growth1_2 = 0. cond1_soa = 0. cond2_soa = 0. cond3_soa = 0. cond4_soa = 0. cond5_soa = 0. cond1_su = 0. ! ---------------------- ! cycle over horizontal domain ! ---------------------- do j = j1, j2 do i = i1, i2 ! ----------------------- ! SURFACE FIELDS ! ----------------------- k = 1 ! ---------------------- ! Physical Parameters ! ---------------------- ! surface pressure [Pa] mixf(region)%f2d(ps)%field(i,j) = mixf(region)%f2d(ps)%field(i,j) + dtime * sp_dat(region)%data(i,j,k) ! precipitation ([m s-1] --> [kg m-2 s-1] with density of water `rol`) mixf(region)%f2d(precip)%field(i,j) = mixf(region)%f2d(precip)%field(i,j) + dtime * (cp_dat(region)%data(i,j,k) + lsp_dat(region)%data(i,j,k)) * rol ! density for conversion of aerosol mass densities ( ---> 1/m3 ) ! (conversion factor 1.E-03 is for g --> kg) dens = pres3d(i,j,k) / temper_dat(region)%data(i,j,k) * xmair * 1.E-3 / (m_dat(region)%data(i,j,k) * Rgas) ! ---------------------- ! Surface Concentrations in [kg m-3] ! ---------------------- !TB added nucleation mode oc ! POM: (AIS + ACS + COS + AII) mixf(region)%f2d(sconcoa)%field(i,j) = mixf(region)%f2d(sconcoa)%field(i,j) + dtime * & dens * (rm(i,j,k,iPOMais) + rm(i,j,k,iPOMacs) + rm(i,j,k,iPOMcos) + rm(i,j,k,iPOMaii)) ! SOA: (NUS + AIS + ACS + COS + AII) mixf(region)%f2d(sconcsoa)%field(i,j) = mixf(region)%f2d(sconcsoa)%field(i,j) + dtime * & dens * (rm(i,j,k,iSOAnus) + rm(i,j,k,iSOAais) + rm(i,j,k,iSOAacs) + rm(i,j,k,iSOAcos) + rm(i,j,k,iSOAaii)) ! BC: (AIS + ACS + COS + AII) mixf(region)%f2d(sconcbc)%field(i,j) = mixf(region)%f2d(sconcbc)%field(i,j) + dtime * & dens * (rm(i,j,k,iBCais) + rm(i,j,k,iBCacs) + rm(i,j,k,iBCcos) + rm(i,j,k,iBCaii)) ! SO4: (NUS + AIS + ACS + COS) mixf(region)%f2d(sconcso4)%field(i,j) = mixf(region)%f2d(sconcso4)%field(i,j) + dtime * & dens * (rm(i,j,k,iSO4nus) + rm(i,j,k,iSO4ais) + rm(i,j,k,iSO4acs) + rm(i,j,k,iSO4cos)) ! Dust: mixf(region)%f2d(sconcdust)%field(i,j) = mixf(region)%f2d(sconcdust)%field(i,j) + dtime * & dens * (rm(i,j,k,iDUacs) + rm(i,j,k,iDUcos) + rm(i,j,k,iDUaci) + rm(i,j,k,iDUcoi)) ! Seasalt: mixf(region)%f2d(sconcss)%field(i,j) = mixf(region)%f2d(sconcss)%field(i,j) + dtime * & dens * (rm(i,j,k,iSSacs) + rm(i,j,k,iSScos)) ! NO3: mixf(region)%f2d(sconcno3)%field(i,j) = mixf(region)%f2d(sconcno3)%field(i,j) + dtime * & dens * rm(i,j,k,iNO3_a) ! NH4: mixf(region)%f2d(sconcnh4)%field(i,j) = mixf(region)%f2d(sconcnh4)%field(i,j) + dtime * & dens * rm(i,j,k,iNH4) ! MSA: mixf(region)%f2d(sconcmsa)%field(i,j) = mixf(region)%f2d(sconcmsa)%field(i,j) + dtime * & dens * rm(i,j,k,iMSA) ! Mass concentrations do nn=68,90 mixf(region)%f2d(nn)%field(i,j) = mixf(region)%f2d(nn)%field(i,j) + dtime * & dens * rm(i,j,k,mixf(region)%f2d(nn)%mf%itm5) end do !number concentrations do nn=91,97 mixf(region)%f2d(nn)%field(i,j) = mixf(region)%f2d(nn)%field(i,j) + dtime * & dens * rm(i,j,k,mixf(region)%f2d(nn)%mf%itm5) end do !h2o h2ocounter=0 do nn=98,101 h2ocounter=h2ocounter+1 mixf(region)%f2d(nn)%field(i,j) = mixf(region)%f2d(nn)%field(i,j) + dtime * & dens * h2o_mode(region,h2ocounter)%d3(i,j,k)!rm(i,j,k,mixf(region)%f2d(nn)%mf%itm5) end do ! RWET for all modes (RDRY=RWET for non-soluble) ! irw=0 do nn=102,108 ! irw=irw+1 mixf(region)%f2d(nn)%field(i,j) = mixf(region)%f2d(nn)%field(i,j) +dtime *rw_mode(region,irw)%d3(i,j,k) end do ! RDRY for soluble modes irwd=0 do nn=109,112 irwd=irwd+1 mixf(region)%f2d(nn)%field(i,j) = mixf(region)%f2d(nn)%field(i,j) +dtime *rwd_mode(region,irwd)%d3(i,j,k) end do ! ---------------------- ! Load in [kg m-2] ! ---------------------- do k = 1, lmr !do ll=1,len(load) ! POM: (AIS + ACS + COS + AII) load_tmp=0.0 do nn=1,7 if (mode_end_pom(nn)>0) then load_tmp = load_tmp+ rm(i,j,k,mode_end_pom(nn)) !(rm(i,j,k,iPOMnus) + rm(i,j,k,iPOMais) + rm(i,j,k,iPOMacs) + rm(i,j,k,iPOMcos) + rm(i,j,k,iPOMaii)) end if end do mixf(region)%f2d(loadoa)%field(i,j) = mixf(region)%f2d(loadoa)%field(i,j) + load_tmp * dtime / dxyp(j) ! SOA: (NUS + AIS + ACS + COS + AII) load_tmp=0.0 do nn=1,7 if (mode_end_soa(nn)>0) then load_tmp = load_tmp+ rm(i,j,k,mode_end_soa(nn)) ENDIF end do mixf(region)%f2d(loadsoa)%field(i,j) = mixf(region)%f2d(loadsoa)%field(i,j) + load_tmp * dtime / dxyp(j) ! BC: (AIS + ACS + COS + AII) load_tmp=0.0 do nn=1,7 if (mode_end_bc(nn)>0) then load_tmp = load_tmp+ rm(i,j,k,mode_end_bc(nn)) !load_tmp = (rm(i,j,k,iBCais) + rm(i,j,k,iBCacs) + rm(i,j,k,iBCcos) + rm(i,j,k,iBCaii)) end if end do mixf(region)%f2d(loadbc)%field(i,j) = mixf(region)%f2d(loadbc)%field(i,j) + load_tmp * dtime / dxyp(j) ! SO4: (NUS + AIS + ACS + COS) load_tmp=0.0 do nn=1,7 if (mode_end_so4(nn)>0) then load_tmp = load_tmp+ rm(i,j,k,mode_end_so4(nn)) !load_tmp = (rm(i,j,k,iSO4nus) + rm(i,j,k,iSO4ais) + rm(i,j,k,iSO4acs) + rm(i,j,k,iSO4cos)) end if end do mixf(region)%f2d(loadso4)%field(i,j) = mixf(region)%f2d(loadso4)%field(i,j) + load_tmp * dtime / dxyp(j) ! Dust: (ACS + COS + ACI + COI) load_tmp=0.0 do nn=1,7 if (mode_end_dust(nn)>0) then load_tmp = load_tmp+ rm(i,j,k,mode_end_dust(nn)) !load_tmp = (rm(i,j,k,iDUacs) + rm(i,j,k,iDUcos) + rm(i,j,k,iDUaci) + rm(i,j,k,iDUcoi)) end if end do mixf(region)%f2d(loaddust)%field(i,j) = mixf(region)%f2d(loaddust)%field(i,j) + load_tmp * dtime / dxyp(j) ! Seasalt: (ACS + COS) load_tmp=0.0 do nn=1,7 if (mode_end_ss(nn)>0) then load_tmp = load_tmp+ rm(i,j,k,mode_end_ss(nn)) !load_tmp = (rm(i,j,k,iSSacs) + rm(i,j,k,iSScos)) end if end do mixf(region)%f2d(loadss)%field(i,j) = mixf(region)%f2d(loadss)%field(i,j) + load_tmp * dtime / dxyp(j) !!$ ! NO3: load_tmp = rm(i,j,k,iNO3_a) mixf(region)%f2d(loadno3)%field(i,j) = mixf(region)%f2d(loadno3)%field(i,j) + load_tmp * dtime / dxyp(j) !!$ end do ! k end do ! i end do ! j !!$ call GO_Timer_Start( itim_opt_out, status ) IF_NOTOK_RETURN(status=1) ! --------------------- ! DO OPTICS IN PARALLEL ! --------------------- ! steps needed for that: ! 1) according to the parallelisation there is different container sizes for ! the results of the optics; these have to be allocated correctly ! (aop_output_ext/a/g/add) ! 2) the optics code is called (init/calculate_aop/done), see documentation ! in the optics module ! 3) the distributed fields (levels/tracers) are reshaped to global fields ! (opt_ext/a/g/add), according to parallelisation ! 4) fields are communicated such that the result contains the right ! total extinctions, albedos, asymmetry factors etc. ! 5) post-calculations (AOD etc.) are done and results are written ! to mixf%../statf% containers as well as to the AOD dump files ! 6) done... ! ------------ REMARKS ! IMPOSSIBLE / TOO EXPENSIVE (from the GENERAL list of parameters for QUICKLOOK) ! - gf @ 90% RH ! - "AOD@550nm SOA", "AOD@550nm BB" ! --------------------------------- ! fill current M7 state into an appropriate container ! and allocate output fields (ext,a,g) ! --------------------------------- ! ----- A T T E N T I O N --------- ! in case for a 'split', we need a field big enough to contain ! various contributions ! (to be synchronously changed with optics_aop_calculate!!) ! ncontr is here number of contributors ! Total, SO4, BC, OC, SS, DU, NO3, Water, Fine, Fine Dust, Fine SS !ncontr = 10 !ncontr = 11 ncontr = 12 ! Additional fields for surface dry aerosol ! to be used for validation against in-situ data: ! extinction, absorption, asymmetry factor ! fine-mode contribution to extinction, absorption nadd = 5 lvec = imr*jmr*lmr ! allocate output fields for optics_calculate_aop Allocate( aop_output_ext(lvec, size(wdep_out), ncontr ) ) ! extinction Allocate( aop_output_a (lvec, size(wdep_out) ) ) ! single scattering albedo Allocate( aop_output_g (lvec, size(wdep_out) ) ) ! asymmetry factor Allocate( aop_output_add(lvec, size(wdep_out), nadd ) ) ! extinction/absorption due to dry aerosol at surface call optics_aop_get(lvec, region, size(wdep_out), wdep_out, ncontr, aop_output_ext, aop_output_a, aop_output_g, & status, aop_output_add ) IF_NOTOK_RETURN(status=1) ! allocate here result arrays for ext, abs, ssa, asymmetry parameter, additional parameters (PM10) ! ncontr is number of contributors for 'split' allocate( opt_ext( i1:i2, j1:j2, lmr, size(wdep_out), ncontr ) ) ; opt_ext = 0.0 allocate( opt_a ( i1:i2, j1:j2, lmr, size(wdep_out) ) ) ; opt_a = 0.0 allocate( opt_g ( i1:i2, j1:j2, lmr, size(wdep_out) ) ) ; opt_g = 0.0 allocate( opt_add( i1:i2, j1:j2, lmr, size(wdep_out), nadd ) ) ; opt_add = 0.0 ! --------------------------------- ! unpack results from calculate_aop ! --------------------------------- do lwl = 1, size(wdep_out) if( wdep_out(lwl)%split ) then ! fill the array for the extinction coefficients of contributors do lct = 1, ncontr opt_ext(:,:,:,lwl,lct) = reshape( aop_output_ext(:,lwl,lct), (/imr,jmr,lmr/) ) end do else ! fill only total extinction coefficients opt_ext(:,:,:,lwl,1) = reshape( aop_output_ext(:,lwl,1), (/imr,jmr,lmr/) ) endif opt_a(:,:,:,lwl) = reshape( aop_output_a(:,lwl),(/imr,jmr,lmr/) ) opt_g(:,:,:,lwl) = reshape( aop_output_g(:,lwl),(/imr,jmr,lmr/) ) if ( wdep_out(lwl)%insitu ) then ! additional fields for split do lct = 1, nadd opt_add(i1:i2,j1:j2,:,lwl,lct) = reshape( aop_output_add(:,lwl,lct),(/imr,jmr,lmr/) ) end do endif end do ! lwl ! free temporary arrays for results from calculate_aop deallocate( aop_output_ext ) deallocate( aop_output_a ) deallocate( aop_output_g ) deallocate( aop_output_add ) ! the following sections assume that for 550nm there is splitted information ! available and that there is *NO* split for other wavelengths! if( count( (wdep_out(:)%wl == 0.550) .and. wdep_out(:)%split ) /= 1 ) then write(gol,*) 'user_output_general: wrong setup with splitted AOD information.'; call goErr TRACEBACK status=1; return end if ! ------------------------------------- ! here additional calculations and ! file writing begin ! ------------------------------------- ! cycle over wavelengths do lwl = 1, size(wdep_out) if( wdep_out(lwl)%split ) then ! for 550nm: ! 1) AODs (35 - 44) ! fill for contributors (total, SO4, BC, POM, SS, DU, NO3, H2O, fine, fine dust, fine SS) ! 2) Absorption for 550nm (45) ! Extinction is here the sum of scattering and absorption and ! the scattering part is given by the albedo (SSA). Thus the ! remainder is absorption: Abs = Ext * (1 - SSA) indoffset = od550aer allocate(temp2d(i1:i2,j1:j2)) do lct = 1, ncontr temp2d = 0.0 do j = j1, j2 do i = i1, i2 ! multiply with height elements and sum up tabs = 0.0 do l = 1, lmr temp2d(i,j) = temp2d(i,j) + opt_ext(i,j,l,lwl,lct) * & (gph_dat(region)%data(i,j,l+1) - & gph_dat(region)%data(i,j,l )) if( lct == 1 ) then ! Absorption: do this only once for the totals tabs = tabs + opt_ext(i,j,l,lwl,lct) * (1.0 - opt_a(i,j,l,lwl)) * & (gph_dat(region)%data(i,j,l+1) - gph_dat(region)%data(i,j,l)) end if end do ! Absorption: do this only once for the totals if( lct == 1 ) then mixf(region)%f2d(abs550aer)%field(i,j) = & mixf(region)%f2d(abs550aer)%field(i,j) + tabs * dtime end if end do end do ! this here is AODs for different contributors, splitted by volumes mixf(region)%f2d(indoffset+lct-1)%field = & mixf(region)%f2d(indoffset+lct-1)%field + temp2d * dtime end do deallocate( temp2d ) ! Asymmetry Parameter, weighted by AOD and single scattering albedo at each layer allocate(temp2d(i1:i2,j1:j2)) ; temp2d = 0.0 do j = j1, j2 do i = i1, i2 wght = 0.0 do l = 1, lmr dwght = opt_ext(i,j,l,lwl,1) * (gph_dat(region)%data(i,j,l+1) - gph_dat(region)%data(i,j,l )) * opt_a(i,j,l,lwl) temp2d(i,j) = temp2d(i,j) + opt_g(i,j,l,lwl) * dwght wght = wght + dwght end do temp2d(i,j) = temp2d(i,j) / wght end do end do mixf(region)%f2d(asyaer)%field = mixf(region)%f2d(asyaer)%field + temp2d * dtime deallocate(temp2d) ! Surface Ambient Aerosol Extinction mixf(region)%f2d(ec550aer)%field = & mixf(region)%f2d(ec550aer)%field + opt_ext(:,:,1,lwl,1) * dtime else ! for 440/870/350 nm: ! 1) fill total AOD and AAOD and write to containers ! 2) dump AOD fields if ( wdep_out(lwl)%wl == 0.440 ) then indoffset = od440aer indabs = abs440aer elseif ( wdep_out(lwl)%wl == 0.550 ) then indoffset = od550aer indabs = abs550aer elseif ( wdep_out(lwl)%wl == 0.870 ) then indoffset = od870aer indabs = abs870aer elseif ( wdep_out(lwl)%wl == 0.350 ) then indoffset = od350aer indabs = abs350aer else write(gol,*) 'user_output_general: wrong wavelength given for general diagnostics' ; call goErr TRACEBACK status=1; return end if ! fill AODs allocate(temp2d(i1:i2,j1:j2)) temp2d = 0.0 do j = j1, j2 do i = i1, i2 ! multiply with height elements and sum up tabs = 0.0 do l = 1, lmr temp2d(i,j) = temp2d(i,j) + opt_ext(i,j,l,lwl,1) * & (gph_dat(region)%data(i,j,l+1) - & gph_dat(region)%data(i,j,l )) tabs = tabs + opt_ext(i,j,l,lwl,1) * (1.0 - opt_a(i,j,l,lwl)) * & (gph_dat(region)%data(i,j,l+1) - gph_dat(region)%data(i,j,l)) end do mixf(region)%f2d(indabs)%field(i,j) = & mixf(region)%f2d(indabs)%field(i,j) + tabs * dtime end do end do ! write to container mixf(region)%f2d(indoffset)%field = & mixf(region)%f2d(indoffset)%field + temp2d * dtime deallocate(temp2d) endif !!$ ! 3-D extinction and absorption (m-1) if ( aai_output ) then if ( wdep_out(lwl)%wl == 0.550 .or. wdep_out(lwl)%wl == 0.350 ) then if ( wdep_out(lwl)%wl == 0.550 ) then indoffset = ec5503Daer elseif ( wdep_out(lwl)%wl == 0.350 ) then indoffset = ec3503Daer endif mixf(region)%f3d(indoffset)%field = & mixf(region)%f3d(indoffset)%field + opt_ext(:,:,:,lwl,1) * dtime mixf(region)%f3d(indoffset+1)%field = & mixf(region)%f3d(indoffset+1)%field + opt_ext(:,:,:,lwl,1) * (1.0 - opt_a(:,:,:,lwl)) * dtime endif endif end do ! done deallocate( opt_ext, opt_a, opt_g, opt_add ) nullify(rm) nullify(dxyp) deallocate( pres3d ) deallocate( CCN ) deallocate( dry_rad) deallocate( Kappa_r) deallocate( N_r) call GO_Timer_End( itim_opt_out, status ) IF_NOTOK_RETURN(status=1) call goLabel() ; status = 0 end subroutine output_general_step !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: COLLECT_BUDGETS ! ! !DESCRIPTION: This routine does book keeping of the budget values in ! order to extract daily contributions to ! emissions / dry deposition / wet deposition. !\\ !\\ ! !INTERFACE: ! subroutine collect_budgets(region, status) ! ! !USES: ! #ifndef without_chemistry use ebischeme, only : buddry_dat => buddep_dat #endif use wet_deposition, only : buddep_dat use emission_data, only : budemi_dat use budget_global, only : nbud_vg use global_data, only : region_dat use chem_param, only : iducoi, iduacs, iducos, iduaci, isscos, issacs use chem_param, only : ibccos, ibcaii, ibcacs, ibcais, ino3_a, xmair use chem_param, only : iso4nus, iso4acs, iso4cos, iso4ais use chem_param, only : ipomcos, ipomaii, ipomacs, ipomais use chem_param, only : isoacos, isoaaii, isoaacs, isoaais, isoanus use chem_param, only : idms, iso2, ntracet, ntrace,iterp,iisop use chem_param, only : xmso2, xmso4, xmdms, xmpom, xmbc, xmdust, xmnacl,xmterp,xmisop #ifndef without_sedimentation use sedimentation, only : buddep_m7_dat use sedimentation, only : nsed, ised #endif ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 29 Nov 2010 - Achim Strunk - ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/collect_budgets' real, dimension(:,:,:,:), allocatable :: collect4d real, dimension(:,:,:), allocatable :: collect3d real, dimension(:,:), allocatable :: tmpbud2d real, dimension(:), pointer :: dxyp integer :: imr, jmr, lmr, j integer :: i1, i2, j1, j2 !>>> TvN integer :: l, n !<<< TvN ! --- begin ------------------------------- call goLabel(rname) ! grid size call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) imr=i2-i1+1 jmr=j2-j1+1 lmr = levi%nlev ! this is area element dxyp => region_dat(region)%dxyp ! -------------- ! EMISSIONS ! -------------- ! collect emission budget allocate( collect4d(i1:i2,j1:j2,nbud_vg,ntracet) ); collect4d = 0.0 ! emissions are in [mole] ! convert first to kilomole per area [k-mole/m2] #ifndef without_emission do j = j1, j2 collect4d(:,j,:,:) = budemi_dat(region)%emi(i1:i2,j,:,:) / dxyp(j) * 1.E-03 end do #endif allocate( tmpbud2d(i1:i2,j1:j2) ) ! on the way: convert from kilomole/m2 to kg/m2 via molar mass [g/mole] ! POM (AIS + ACS + COS + AII) tmpbud2d = ( sum(collect4d(:,:,:,iPOMais),3) + sum(collect4d(:,:,:,iPOMacs),3) + & sum(collect4d(:,:,:,iPOMcos),3) + sum(collect4d(:,:,:,iPOMaii),3) ) * xmpom mixf(region)%f2d(emioa)%field = tmpbud2d - emission(region)%f2dslast(:,:,1) emission(region)%f2dslast(:,:,1) = tmpbud2d ! BC (AIS + ACS + COS + AII) tmpbud2d = ( sum(collect4d(:,:,:,iBCais ),3) + sum(collect4d(:,:,:,iBCacs ),3) + & sum(collect4d(:,:,:,iBCcos ),3) + sum(collect4d(:,:,:,iBCaii ),3) ) * xmbc mixf(region)%f2d(emibc)%field = tmpbud2d - emission(region)%f2dslast(:,:,2) emission(region)%f2dslast(:,:,2) = tmpbud2d ! SO2 tmpbud2d = sum(collect4d(:,:,:,iSO2 ),3) * xmso2 mixf(region)%f2d(emiso2)%field = tmpbud2d - emission(region)%f2dslast(:,:,3) emission(region)%f2dslast(:,:,3) = tmpbud2d ! SO4 (NUS + AIS + ACS + COS) tmpbud2d = ( sum(collect4d(:,:,:,iSO4nus),3) + sum(collect4d(:,:,:,iSO4ais),3) + & sum(collect4d(:,:,:,iSO4acs),3) + sum(collect4d(:,:,:,iSO4cos),3) ) * xmso4 mixf(region)%f2d(emiso4)%field = tmpbud2d - emission(region)%f2dslast(:,:,4) emission(region)%f2dslast(:,:,4) = tmpbud2d ! Dust (ACS + COS + ACI + COI) tmpbud2d = ( sum(collect4d(:,:,:,iDUacs ),3) + sum(collect4d(:,:,:,iDUcos ),3) + & sum(collect4d(:,:,:,iDUaci ),3) + sum(collect4d(:,:,:,iDUcoi ),3) ) * xmdust mixf(region)%f2d(emidust)%field = tmpbud2d - emission(region)%f2dslast(:,:,5) emission(region)%f2dslast(:,:,5) = tmpbud2d ! DMS tmpbud2d = sum(collect4d(:,:,:,iDMS ),3) * xmdms mixf(region)%f2d(emidms)%field = tmpbud2d - emission(region)%f2dslast(:,:,6) emission(region)%f2dslast(:,:,6) = tmpbud2d ! Seasalt: (ACS + COS) tmpbud2d = ( sum(collect4d(:,:,:,iSSacs ),3) + sum(collect4d(:,:,:,iSScos ),3) ) * xmnacl mixf(region)%f2d(emiss)%field = tmpbud2d - emission(region)%f2dslast(:,:,7) emission(region)%f2dslast(:,:,7) = tmpbud2d ! terpene: tmpbud2d = sum(collect4d(:,:,:,iterp ),3) * xmterp mixf(region)%f2d(emiterp)%field = tmpbud2d - emission(region)%f2dslast(:,:,8) emission(region)%f2dslast(:,:,8) = tmpbud2d ! isoprene: tmpbud2d = sum(collect4d(:,:,:,iisop ),3) * xmisop mixf(region)%f2d(emiisop)%field = tmpbud2d - emission(region)%f2dslast(:,:,9) emission(region)%f2dslast(:,:,9) = tmpbud2d deallocate( tmpbud2d ) deallocate( collect4d ) ! -------------- ! DRY DEPOSITION ! -------------- allocate( collect3d(i1:i2,j1:j2,ntrace) ); collect3d = 0.0 ! deposition is in [mole] ! convert first to kilomole per area [k-mole/m2] do j = j1, j2 #ifndef without_chemistry collect3d(:,j,:) = buddry_dat(region)%dry(i1:i2,j,:) / dxyp(j) * 1.E-3 #endif ! Add sedimentation at the surface to dry depostion ! Sedimentation budgets have to be summed in the vertical ! to get the surface contribution. #ifndef without_sedimentation do l = 1, nbud_vg do n = 1, nsed collect3d(:,j,ised(n)) = collect3d(:,j,ised(n)) + & buddep_m7_dat(region)%sed(i1:i2,j,l,n) / dxyp(j) * 1.E-3 end do end do #endif end do allocate( tmpbud2d(i1:i2,j1:j2) ) ! on the way: convert from kilomole/m2 to kg/m2 via molar mass [g/mole] ! SO2 tmpbud2d = collect3d(:,:,iSO2) * xmso2 mixf(region)%f2d(dryso2)%field = tmpbud2d - drydepos(region)%f2dslast(:,:,1) drydepos(region)%f2dslast(:,:,1) = tmpbud2d ! POM (AIS + ACS + COS + AII) tmpbud2d = ( collect3d(:,:,iPOMais ) + collect3d(:,:,iPOMacs ) + & collect3d(:,:,iPOMcos ) + collect3d(:,:,iPOMaii )) * xmpom mixf(region)%f2d(dryoa)%field = tmpbud2d - drydepos(region)%f2dslast(:,:,2) drydepos(region)%f2dslast(:,:,2) = tmpbud2d ! SOA (NUS + AIS + ACS + COS + AII) tmpbud2d = ( collect3d(:,:,iSOAais ) + collect3d(:,:,iSOAacs ) + & collect3d(:,:,iSOAcos ) + collect3d(:,:,iSOAaii ) + & collect3d(:,:,iSOAnus )) * xmpom mixf(region)%f2d(drysoa)%field = tmpbud2d - drydepos(region)%f2dslast(:,:,8) drydepos(region)%f2dslast(:,:,8) = tmpbud2d ! BC (AIS + ACS + COS + AII) tmpbud2d = ( collect3d(:,:,iBCais ) + collect3d(:,:,iBCacs ) + & collect3d(:,:,iBCcos ) + collect3d(:,:,iBCaii ) ) * xmbc mixf(region)%f2d(drybc)%field = tmpbud2d - drydepos(region)%f2dslast(:,:,3) drydepos(region)%f2dslast(:,:,3) = tmpbud2d ! SO4 (NUS + AIS + ACS + COS) tmpbud2d = ( collect3d(:,:,iSO4nus) + collect3d(:,:,iSO4ais) + & collect3d(:,:,iSO4acs) + collect3d(:,:,iSO4cos) ) * xmso4 mixf(region)%f2d(dryso4)%field = tmpbud2d - drydepos(region)%f2dslast(:,:,4) drydepos(region)%f2dslast(:,:,4) = tmpbud2d ! Dust (ACS + COS + ACI + COI) tmpbud2d = ( collect3d(:,:,iDUacs) + collect3d(:,:,iDUcos) + & collect3d(:,:,iDUaci) + collect3d(:,:,iDUcoi) ) * xmdust mixf(region)%f2d(drydust)%field = tmpbud2d - drydepos(region)%f2dslast(:,:,5) drydepos(region)%f2dslast(:,:,5) = tmpbud2d ! DMS ! DMS is NOT deposited in TM5 --> the fields will contain zeros tmpbud2d = collect3d(:,:,iDMS) * xmdms mixf(region)%f2d(drydms)%field = tmpbud2d - drydepos(region)%f2dslast(:,:,6) drydepos(region)%f2dslast(:,:,6) = tmpbud2d ! Seasalt: (ACS + COS) tmpbud2d = ( collect3d(:,:,iSSacs) + collect3d(:,:,iSScos) ) * xmnacl mixf(region)%f2d(dryss)%field = tmpbud2d - drydepos(region)%f2dslast(:,:,7) drydepos(region)%f2dslast(:,:,7) = tmpbud2d deallocate( tmpbud2d ) deallocate( collect3d ) ! -------------- ! WET DEPOSITION ! -------------- allocate( collect4d (i1:i2,j1:j2,nbud_vg,ntracet) ); collect4d = 0.0 ! deposition is in [mole] ! convert first to kilomole per area [k-mole/m2] do j = j1, j2 collect4d(:,j,:,:) = ( buddep_dat(region)%lsp(i1:i2,j,:,:) + buddep_dat(region)%cp(i1:i2,j,:,:) ) / dxyp(j) * 1.E-3 end do allocate( tmpbud2d(i1:i2,j1:j2) ) ! on the way: convert from kilomole/m2 to kg/m2 via molar mass [g/mole] ! POM (AIS + ACS + COS + AII) tmpbud2d = ( sum(collect4d(:,:,:,iPOMais),3) + sum(collect4d(:,:,:,iPOMacs),3) + & sum(collect4d(:,:,:,iPOMcos),3) + sum(collect4d(:,:,:,iPOMaii),3)) * xmpom mixf(region)%f2d(wetoa)%field = tmpbud2d - wetdepos(region)%f2dslast(:,:,1) wetdepos(region)%f2dslast(:,:,1) = tmpbud2d ! SOA (NUS + AIS + ACS + COS + AII) tmpbud2d = ( sum(collect4d(:,:,:,iSOAais),3) + sum(collect4d(:,:,:,iSOAacs),3) + & sum(collect4d(:,:,:,iSOAcos),3) + sum(collect4d(:,:,:,iSOAaii),3) + & sum(collect4d(:,:,:,iSOAnus),3)) * xmpom mixf(region)%f2d(wetsoa)%field = tmpbud2d - wetdepos(region)%f2dslast(:,:,8) wetdepos(region)%f2dslast(:,:,8) = tmpbud2d ! BC (AIS + ACS + COS + AII) tmpbud2d = ( sum(collect4d(:,:,:,iBCais ),3) + sum(collect4d(:,:,:,iBCacs ),3) + & sum(collect4d(:,:,:,iBCcos ),3) + sum(collect4d(:,:,:,iBCaii ),3) ) * xmbc mixf(region)%f2d(wetbc)%field = tmpbud2d - wetdepos(region)%f2dslast(:,:,2) wetdepos(region)%f2dslast(:,:,2) = tmpbud2d ! SO2 tmpbud2d = sum(collect4d(:,:,:,iSO2 ),3) * xmso2 mixf(region)%f2d(wetso2)%field = tmpbud2d - wetdepos(region)%f2dslast(:,:,3) wetdepos(region)%f2dslast(:,:,3) = tmpbud2d ! SO4 (NUS + AIS + ACS + COS) tmpbud2d = ( sum(collect4d(:,:,:,iSO4nus),3) + sum(collect4d(:,:,:,iSO4ais),3) + & sum(collect4d(:,:,:,iSO4acs),3) + sum(collect4d(:,:,:,iSO4cos),3) ) * xmso4 mixf(region)%f2d(wetso4)%field = tmpbud2d - wetdepos(region)%f2dslast(:,:,4) wetdepos(region)%f2dslast(:,:,4) = tmpbud2d ! Dust (ACS + COS + ACI + COI) tmpbud2d = ( sum(collect4d(:,:,:,iDUacs ),3) + sum(collect4d(:,:,:,iDUcos ),3) + & sum(collect4d(:,:,:,iDUaci ),3) + sum(collect4d(:,:,:,iDUcoi ),3) ) * xmdust mixf(region)%f2d(wetdust)%field = tmpbud2d - wetdepos(region)%f2dslast(:,:,5) wetdepos(region)%f2dslast(:,:,5) = tmpbud2d ! DMS ! DMS is NOT deposited in TM5 --> the fields will contain zeros tmpbud2d = sum(collect4d(:,:,:,iDMS ),3) * xmdms mixf(region)%f2d(wetdms)%field = tmpbud2d - wetdepos(region)%f2dslast(:,:,6) wetdepos(region)%f2dslast(:,:,6) = tmpbud2d ! Seasalt: (ACS + COS) tmpbud2d = ( sum(collect4d(:,:,:,iSSacs ),3) + sum(collect4d(:,:,:,iSScos ),3) ) * xmnacl mixf(region)%f2d(wetss)%field = tmpbud2d - wetdepos(region)%f2dslast(:,:,7) wetdepos(region)%f2dslast(:,:,7) = tmpbud2d deallocate( tmpbud2d ) deallocate( collect4d ) nullify(dxyp) call goLabel() ; status = 0 end subroutine collect_budgets !EOC subroutine add_variable(region,itm5,varname,longname,unit,dims,n_out,status) use chem_param, only: mode_end_so4,mode_end_pom,mode_end_bc,mode_end_ss,mode_end_dust implicit none integer:: itm5 integer:: fileunit integer:: varid character(len=*)::varname character(len=*)::longname character(len=*)::unit integer:: dims integer:: region integer,intent(out)::n_out integer,intent(out)::status logical :: writeout=.true. character(len=*), parameter :: rname = mname//'/output_aerchemmip_add_variable' integer ::i1,i2,j1,j2,jmr,imr,lmr call goLabel(rname) !call GO_Timer_Start( itim_addvar, status ) !IF_NOTOK_RETURN(status=1) !write(gol,'("add_variable1")') !call goErr if( (n_2d_vars+1>ntracer_2d) .or.(n_3d_vars+1>ntracer_3d) ) then status=1 IF_ERROR_RETURN(status=1) end if call Get_DistGrid( dgrid(1), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) ! define sizes for arrays imr=i2-i1+1 jmr=j2-j1+1 lmr = levi%nlev ! 2D variables if (dims==2) then !init variable n_2d_vars=n_2d_vars+1 mixf(region)%f2d(n_2d_vars)%mf = metafields( itm5 , varname,longname,unit,'','') mixf(region)%f2d(n_2d_vars)%writeout=writeout !Later make the allocation here !allocate( mixf(region)%f2d(n_2d_vars)%field(i1:i2,j1:j2) ) mixf(region)%f2d(n_2d_vars)%field(i1:i2,j1:j2)=0.0 n_out=n_2d_vars else if (dims==3) then !init variable n_3d_vars=n_3d_vars+1 mixf(region)%f3d(n_3d_vars)%mf = metafields( itm5 , varname,longname,unit,'','') mixf(region)%f3d(n_3d_vars)%writeout=writeout !Later make the allocation here !allocate( mixf(region)%f3d(n_3d_vars)%field(i1:i2,j1:j2,lmr) ) mixf(region)%f3d(n_3d_vars)%field(i1:i2,j1:j2,lmr)=0.0 n_out=n_3d_vars end if call goLabel() status = 0 !call GO_Timer_End( itim_addvar, status ) !IF_NOTOK_RETURN(status=1) end subroutine add_variable subroutine CCN_Diag(CCN,mass_num,cmr,temper) use chem_param, only : iso4ais, iso4acs, iso4cos, ibcais, ibcacs, ibccos, ipomais, ipomacs use chem_param, only : ipomcos, issacs, isscos, iduacs, iducos, isoaais, isoaacs, isoacos use chem_param, only : iais_n, iacs_n, icos_n, sigma_lognormal use chem_param, only : ino3_a, imsa, nh4no3_factor, nh4no3_density, msa_density use chem_param, only : SurfTen, mode_tracers use chem_param, only : Kap_su, Kap_pom, Kap_soa, Kap_bc, Kap_du, Kap_ss, Kap_na2so4, Kap_no3, Kap_msa use chem_param, only : so4_density, pom_density, soa_density, carbon_density, ss_density, dust_density use mo_aero_m7, only : cmr2ram, nsol, ram2cmr use mo_aero_m7, only : dNaCl, dH2SO4, dNa2SO4, wNaCl, wH2SO4, wNa2SO4, wSO4 use binas, only : pi, Rgas, xm_h2o, rol implicit none integer :: iSat, iMode real :: CCN(:,:,:,:), mass_num(:,:,:,:), cmr(:,:,:,:), temper(:,:,:) real, allocatable :: Vol(:,:,:), Kappa(:,:,:), A_Koehler(:,:,:) real, allocatable :: N_act(:,:,:), rc(:,:,:) real, allocatable :: nNa(:,:,:), nCl(:,:,:), nNaCl(:,:,:), nH2SO4(:,:,:), nNa2SO4(:,:,:), nSO4(:,:,:) CCN(:,:,:,:) = 0.e0 cmr(:,:,:,:) = 0.e0 allocate(Vol(size(CCN,2),size(CCN,3),size(CCN,4))) allocate(Kappa(size(CCN,2),size(CCN,3),size(CCN,4))) allocate(A_Koehler(size(CCN,2),size(CCN,3),size(CCN,4))) allocate(N_act(size(CCN,2),size(CCN,3),size(CCN,4))) allocate(rc(size(CCN,2),size(CCN,3),size(CCN,4))) allocate(nNa(size(CCN,2),size(CCN,3),size(CCN,4))) allocate(nCl(size(CCN,2),size(CCN,3),size(CCN,4))) allocate(nSO4(size(CCN,2),size(CCN,3),size(CCN,4))) allocate(nNaCl(size(CCN,2),size(CCN,3),size(CCN,4))) allocate(nNa2SO4(size(CCN,2),size(CCN,3),size(CCN,4))) allocate(nH2SO4(size(CCN,2),size(CCN,3),size(CCN,4))) do iMode = 2,nsol ! Kappa SELECT CASE (iMode) CASE(2) ! no sea salt and dust in soluble aitken mode nSO4= mass_num(:,:,:,iso4ais)/(1.e-3*wSO4) nH2SO4 = nSO4 Vol = ( nH2SO4*1.e-3*wH2SO4/(dH2SO4*1.e3)*(wH2SO4/wSO4) & !mass_num(:,:,:,iso4ais)/so4_density & + mass_num(:,:,:,ipomais)/pom_density & + mass_num(:,:,:,isoaais)/soa_density & + mass_num(:,:,:,ibcais)/carbon_density) ! m3/gridbox Kappa= ( Kap_su*nH2SO4*1.e-3*wH2SO4/(dH2SO4*1.e3)*(wH2SO4/wSO4) &!mass_num(:,:,:,iso4ais)/so4_density & + Kap_pom*mass_num(:,:,:,ipomais)/pom_density & + Kap_soa*mass_num(:,:,:,isoaais)/soa_density & + Kap_bc*mass_num(:,:,:,ibcais)/carbon_density) / Vol cmr(iMode,:,:,:) = ram2cmr(iMode)*(Vol/mass_num(:,:,:,iais_n)*0.75/pi)**(1./3.) ! m CASE(3) nNa = mass_num(:,:,:,issacs)/(1.e-3*wNaCl) ! mol/gridbox nCl = nNa nSO4= mass_num(:,:,:,iso4acs)/(1.e-3*wSO4) nNa2SO4 = MIN(nNa/2.e0, nSO4) nNa = nNa - 2.e0*nNa2SO4 ! remaining nNa nNaCl = MIN(nCl, nNa) nCl = nNaCl ! overshoot nCl is supposed to evaporate nH2SO4 = nSO4 - nNa2SO4 Vol = ( nNaCl*1.e-3*wNaCl/(dNaCl*1.e3) & + nNa2SO4*1.e-3*wNa2SO4/(dNa2SO4*1.e3) & + nH2SO4*1.e-3*wH2SO4/(dH2SO4*1.e3) & + mass_num(:,:,:,ipomacs)/pom_density & + mass_num(:,:,:,isoaacs)/soa_density & + mass_num(:,:,:,ibcacs)/carbon_density & + mass_num(:,:,:,iduacs)/dust_density & + mass_num(:,:,:,ino3_a)*nh4no3_factor/nh4no3_density & + mass_num(:,:,:,imsa)/msa_density) ! m3/gridbox Kappa= ( Kap_su*nH2SO4*1.e-3*wH2SO4/(dH2SO4*1.e3) & + Kap_na2so4*nNa2SO4*1.e-3*wNa2SO4/(dNa2SO4*1.e3) & + Kap_pom*mass_num(:,:,:,ipomacs)/pom_density & + Kap_soa*mass_num(:,:,:,isoaacs)/soa_density & + Kap_bc*mass_num(:,:,:,ibcacs)/carbon_density & + Kap_ss*nNaCl*1.e-3*wNaCl/(dNaCl*1.e3) & + Kap_du*mass_num(:,:,:,iduacs)/dust_density & + Kap_no3*mass_num(:,:,:,ino3_a)*nh4no3_factor/nh4no3_density & + Kap_msa*mass_num(:,:,:,imsa)/msa_density) / Vol cmr(iMode,:,:,:) = ( Vol/mass_num(:,:,:,iacs_n)*0.75/pi)**(1./3.)*ram2cmr(iMode) ! m CASE(4) nNa = mass_num(:,:,:,isscos)/(1.e-3*wNaCl) ! mol/gridbox nCl = nNa nSO4= mass_num(:,:,:,iso4cos)/(1.e-3*wSO4) nNa2SO4 = MIN(nNa/2.e0, nSO4) nNa = nNa - 2.e0*nNa2SO4 ! remaining nNa nNaCl = MIN(nCl, nNa) nCl = nNaCl ! overshoot nCl is supposed to evaporate nH2SO4 = nSO4 - nNa2SO4 Vol = ( nNaCl*1.e-3*wNaCl/(dNaCl*1.e3) & + nNa2SO4*1.e-3*wNa2SO4/(dNa2SO4*1.e3) & + nH2SO4*1.e-3*wH2SO4/(dH2SO4*1.e3) & + mass_num(:,:,:,ipomcos)/pom_density & + mass_num(:,:,:,isoacos)/soa_density & + mass_num(:,:,:,ibccos)/carbon_density & + mass_num(:,:,:,iducos)/dust_density) ! m3/gridbox Kappa= ( Kap_su*nH2SO4*1.e-3*wH2SO4/(dH2SO4*1.e3) & + Kap_na2so4*nNa2SO4*1.e-3*wNa2SO4/(dNa2SO4*1.e3) & + Kap_pom*mass_num(:,:,:,ipomcos)/pom_density & + Kap_soa*mass_num(:,:,:,isoacos)/soa_density & + Kap_bc*mass_num(:,:,:,ibccos)/carbon_density & + Kap_ss*nNaCl*1.e-3*wNaCl/(dNaCl*1.e3) & + Kap_du*mass_num(:,:,:,iducos)/dust_density) / Vol cmr(iMode,:,:,:) = ( Vol/mass_num(:,:,:,icos_n)*0.75/pi)**(1./3.)*ram2cmr(iMode) ! m END SELECT where(Kappa < 0.04) Kappa = 0.04 A_Koehler = 2.0 * xm_h2o * SurfTen / Rgas / rol / temper do iSat=1,nSat ! critical radius per mode rc = A_Koehler/3.0 * (2.0/SuperSat(iSat)/sqrt(Kappa))**(2.0/3.0) ! number of activated particles N_act = mass_num(:,:,:,mode_tracers(0,iMode)) * 0.5*erfc( log(rc/cmr(iMode,:,:,:)) / & sqrt(2.e0)/log(sigma_lognormal(iMode)) ) CCN(iSat,:,:,:) = CCN(iSat,:,:,:) + N_act ! #/gridbox enddo enddo deallocate(Vol) deallocate(Kappa) deallocate(A_Koehler) deallocate(rc) deallocate(N_act) end subroutine CCN_Diag subroutine Kappa_diag(CCN,mass_num,cmr,temper,N_r,Kappa_r) use chem_param, only : iso4ais, iso4acs, iso4cos, ibcais, ibcacs, ibccos, ipomais, ipomacs use chem_param, only : ipomcos, issacs, isscos, iduacs, iducos, isoaais, isoaacs, isoacos use chem_param, only : iais_n, iacs_n, icos_n, sigma_lognormal use chem_param, only : ino3_a, imsa, nh4no3_factor, nh4no3_density, msa_density use chem_param, only : SurfTen, mode_tracers use chem_param, only : Kap_su, Kap_pom, Kap_soa, Kap_bc, Kap_du, Kap_ss, Kap_na2so4, Kap_no3, Kap_msa use chem_param, only : so4_density, pom_density, soa_density, carbon_density, ss_density, dust_density use mo_aero_m7, only : cmr2ram, nsol, ram2cmr use mo_aero_m7, only : dNaCl, dH2SO4, dNa2SO4, wNaCl, wH2SO4, wNa2SO4, wSO4 use binas, only : pi, Rgas, xm_h2o, rol implicit none integer :: iSat, iMode,i_rad real :: CCN(:,:,:,:), mass_num(:,:,:,:), cmr(:,:,:,:), temper(:,:,:) real :: N_r(:,:,:,:), Kappa_r(:,:,:,:) real, allocatable :: Vol(:,:,:), Kappa(:,:,:), A_Koehler(:,:,:) real, allocatable :: N_act(:,:,:), rc(:) real, allocatable :: nNa(:,:,:), nCl(:,:,:), nNaCl(:,:,:), nH2SO4(:,:,:), nNa2SO4(:,:,:), nSO4(:,:,:) real, allocatable :: Volss(:,:,:) real, allocatable :: Voldu(:,:,:) real, allocatable :: Volsu(:,:,:) real, allocatable :: Volbc(:,:,:) real, allocatable :: Volpom(:,:,:) real, allocatable :: Volsoa(:,:,:) real, allocatable :: Volna2so4(:,:,:) real, allocatable :: Volmsa(:,:,:) real, allocatable :: Volno3(:,:,:) real, allocatable :: Volsum(:,:,:) real, allocatable :: over_rad_frac(:,:,:) ! CCN(:,:,:,:) = 0.e0 cmr(:,:,:,:) = 0.e0 N_r(:,:,:,:) = 0.e0 Kappa_r(:,:,:,:) = 0.e0 allocate(Vol(size(CCN,2),size(CCN,3),size(CCN,4))) allocate(Volsu(size(CCN,2),size(CCN,3),size(CCN,4))) allocate(Volss(size(CCN,2),size(CCN,3),size(CCN,4))) allocate(Voldu(size(CCN,2),size(CCN,3),size(CCN,4))) allocate(Volbc(size(CCN,2),size(CCN,3),size(CCN,4))) allocate(Volsoa(size(CCN,2),size(CCN,3),size(CCN,4))) allocate(Volpom(size(CCN,2),size(CCN,3),size(CCN,4))) allocate(Volna2so4(size(CCN,2),size(CCN,3),size(CCN,4))) allocate(Volmsa(size(CCN,2),size(CCN,3),size(CCN,4))) allocate(Volno3(size(CCN,2),size(CCN,3),size(CCN,4))) allocate(Volsum(size(CCN,2),size(CCN,3),size(CCN,4))) allocate(over_rad_frac(size(CCN,2),size(CCN,3),size(CCN,4))) volsum(:,:,:)=0.e0 volss(:,:,:)=0.e0 voldu(:,:,:)=0.e0 volbc(:,:,:)=0.e0 volsu(:,:,:)=0.e0 volpom(:,:,:)=0.e0 volsoa(:,:,:)=0.e0 volna2so4(:,:,:)=0.e0 volmsa(:,:,:)=0.e0 volno3(:,:,:)=0.e0 over_rad_frac(:,:,:)=0.e0 ! allocate(N_r(3,size(CCN,2),size(CCN,3),size(CCN,4))) ! allocate(Kappa_r(size(CCN,2),size(CCN,3),size(CCN,4))) allocate(Kappa(size(CCN,2),size(CCN,3),size(CCN,4))) allocate(A_Koehler(size(CCN,2),size(CCN,3),size(CCN,4))) allocate(N_act(size(CCN,2),size(CCN,3),size(CCN,4))) ! allocate(rc(size(CCN,2),size(CCN,3),size(CCN,4))) allocate(rc(3)) allocate(nNa(size(CCN,2),size(CCN,3),size(CCN,4))) allocate(nCl(size(CCN,2),size(CCN,3),size(CCN,4))) allocate(nSO4(size(CCN,2),size(CCN,3),size(CCN,4))) allocate(nNaCl(size(CCN,2),size(CCN,3),size(CCN,4))) allocate(nNa2SO4(size(CCN,2),size(CCN,3),size(CCN,4))) allocate(nH2SO4(size(CCN,2),size(CCN,3),size(CCN,4))) !Diam for N !BACCHUS defs D=50,80,120nm->r=25,40,60 rc(1)= 25.0e-9 rc(2)= 40.0e-9 rc(3)= 60.0e-9 do i_rad =1,3 !neglect nucleation mode do iMode = 2,nsol ! Kappa SELECT CASE (iMode) CASE(2) ! no sea salt and dust in soluble aitken mode nSO4= mass_num(:,:,:,iso4ais)/(1.e-3*wSO4) nH2SO4 = nSO4 Vol = ( nH2SO4*1.e-3*wH2SO4/(dH2SO4*1.e3) &! mass_num(:,:,:,iso4ais)/so4_density & + mass_num(:,:,:,ipomais)/pom_density & + mass_num(:,:,:,isoaais)/soa_density & + mass_num(:,:,:,ibcais)/carbon_density) ! m3/gridbox Kappa= ( Kap_su*nH2SO4*1.e-3*wH2SO4/(dH2SO4*1.e3) &!mass_num(:,:,:,iso4ais)/so4_density & + Kap_pom*mass_num(:,:,:,ipomais)/pom_density & + Kap_soa*mass_num(:,:,:,isoaais)/soa_density & + Kap_bc*mass_num(:,:,:,ibcais)/carbon_density) / Vol cmr(iMode,:,:,:) = ram2cmr(iMode)*(Vol/mass_num(:,:,:,iais_n)*0.75/pi)**(1./3.) ! m !fraction over threshold over_rad_frac= 0.5*erfc( log(rc(i_rad)/cmr(iMode,:,:,:)) / & sqrt(2.e0)/log(sigma_lognormal(iMode)) ) volsu= volsu + nH2SO4*1.e-3*wH2SO4/(dH2SO4*1.e3) * over_rad_frac volpom = volpom + mass_num(:,:,:,ipomais)/pom_density * over_rad_frac volsoa = volsoa + mass_num(:,:,:,isoaais)/soa_density * over_rad_frac volbc = volbc + mass_num(:,:,:,ibcais)/carbon_density * over_rad_frac N_r(i_rad,:,:,:) = N_r(i_rad,:,:,:)+mass_num(:,:,:,mode_tracers(0,iMode)) * 0.5*erfc( log(rc(i_rad)/cmr(iMode,:,:,:)) / & sqrt(2.e0)/log(sigma_lognormal(iMode)) ) CASE(3) nNa = mass_num(:,:,:,issacs)/(1.e-3*wNaCl) ! mol/gridbox nCl = nNa nSO4= mass_num(:,:,:,iso4acs)/(1.e-3*wSO4) nNa2SO4 = MIN(nNa/2.e0, nSO4) nNa = nNa - 2.e0*nNa2SO4 ! remaining nNa nNaCl = MIN(nCl, nNa) nCl = nNaCl ! overshoot nCl is supposed to evaporate nH2SO4 = nSO4 - nNa2SO4 Vol = ( nNaCl*1.e-3*wNaCl/(dNaCl*1.e3) & + nNa2SO4*1.e-3*wNa2SO4/(dNa2SO4*1.e3) & + nH2SO4*1.e-3*wH2SO4/(dH2SO4*1.e3) & + mass_num(:,:,:,ipomacs)/pom_density & + mass_num(:,:,:,isoaacs)/soa_density & + mass_num(:,:,:,ibcacs)/carbon_density & + mass_num(:,:,:,iduacs)/dust_density & + mass_num(:,:,:,ino3_a)*nh4no3_factor/nh4no3_density & + mass_num(:,:,:,imsa)/msa_density) ! m3/gridbox Kappa= ( Kap_su*nH2SO4*1.e-3*wH2SO4/(dH2SO4*1.e3) & + Kap_na2so4*nNa2SO4*1.e-3*wNa2SO4/(dNa2SO4*1.e3) & + Kap_pom*mass_num(:,:,:,ipomacs)/pom_density & + Kap_soa*mass_num(:,:,:,isoaacs)/soa_density & + Kap_bc*mass_num(:,:,:,ibcacs)/carbon_density & + Kap_ss*nNaCl*1.e-3*wNaCl/(dNaCl*1.e3) & + Kap_du*mass_num(:,:,:,iduacs)/dust_density & + Kap_no3*mass_num(:,:,:,ino3_a)*nh4no3_factor/nh4no3_density & + Kap_msa*mass_num(:,:,:,imsa)/msa_density) / Vol cmr(iMode,:,:,:) = ( Vol/mass_num(:,:,:,iacs_n)*0.75/pi)**(1./3.)*ram2cmr(iMode) ! m over_rad_frac= 0.5*erfc( log(rc(i_rad)/cmr(iMode,:,:,:)) / & sqrt(2.e0)/log(sigma_lognormal(iMode)) ) volsu= volsu + nH2SO4*1.e-3*wH2SO4/(dH2SO4*1.e3) * over_rad_frac volna2so4= volna2so4 + nNa2SO4*1.e-3*wNa2SO4/(dNa2SO4*1.e3) * over_rad_frac volpom = volpom + mass_num(:,:,:,ipomacs)/pom_density * over_rad_frac volsoa = volsoa + mass_num(:,:,:,isoaacs)/soa_density * over_rad_frac volbc = volbc + mass_num(:,:,:,ibcacs)/carbon_density * over_rad_frac volss =volss + nNaCl*1.e-3*wNaCl/(dNaCl*1.e3) * over_rad_frac voldu = voldu + mass_num(:,:,:,iduacs)/dust_density * over_rad_frac volno3 = volno3 + mass_num(:,:,:,ino3_a)*nh4no3_factor/nh4no3_density * over_rad_frac volmsa = volmsa + mass_num(:,:,:,imsa)/msa_density * over_rad_frac N_r(i_rad,:,:,:) = N_r(i_rad,:,:,:)+mass_num(:,:,:,mode_tracers(0,iMode)) * 0.5*erfc( log(rc(i_rad)/cmr(iMode,:,:,:)) / & sqrt(2.e0)/log(sigma_lognormal(iMode)) ) CASE(4) nNa = mass_num(:,:,:,isscos)/(1.e-3*wNaCl) ! mol/gridbox nCl = nNa nSO4= mass_num(:,:,:,iso4cos)/(1.e-3*wSO4) nNa2SO4 = MIN(nNa/2.e0, nSO4) nNa = nNa - 2.e0*nNa2SO4 ! remaining nNa nNaCl = MIN(nCl, nNa) nCl = nNaCl ! overshoot nCl is supposed to evaporate nH2SO4 = nSO4 - nNa2SO4 Vol = ( nNaCl*1.e-3*wNaCl/(dNaCl*1.e3) & + nNa2SO4*1.e-3*wNa2SO4/(dNa2SO4*1.e3) & + nH2SO4*1.e-3*wH2SO4/(dH2SO4*1.e3) & + mass_num(:,:,:,ipomcos)/pom_density & + mass_num(:,:,:,isoacos)/soa_density & + mass_num(:,:,:,ibccos)/carbon_density & + mass_num(:,:,:,iducos)/dust_density) ! m3/gridbox cmr(iMode,:,:,:) = ( Vol/mass_num(:,:,:,icos_n)*0.75/pi)**(1./3.)*ram2cmr(iMode)! m volsu= volsu + nH2SO4*1.e-3*wH2SO4/(dH2SO4*1.e3) volna2so4= volna2so4 + nNa2SO4*1.e-3*wNa2SO4/(dNa2SO4*1.e3) volpom = volpom + mass_num(:,:,:,ipomcos)/pom_density volsoa = volsoa + mass_num(:,:,:,isoacos)/soa_density volbc = volbc + mass_num(:,:,:,ibccos)/carbon_density volss =volss + nNaCl*1.e-3*wNaCl/(dNaCl*1.e3) voldu = voldu + mass_num(:,:,:,iducos)/dust_density N_r(i_rad,:,:,:) = N_r(i_rad,:,:,:)+mass_num(:,:,:,mode_tracers(0,iMode)) * 0.5*erfc( log(rc(i_rad)/cmr(iMode,:,:,:)) / & sqrt(2.e0)/log(sigma_lognormal(iMode)) ) END SELECT volsum=volsum+vol end do Kappa_r(i_rad,:,:,:) = (& Kap_su*volsu& + Kap_na2so4*volna2so4& + Kap_pom*volpom& + Kap_soa*volsoa& + Kap_bc*volbc& + Kap_ss*volss& + Kap_no3*volno3& + Kap_du*voldu& ) / volsum !where (Kappa_r>1.0)Kappa_r=1.0 ! where(Kappa_r < 0.04) Kappa_r = 0.04 ! Kappa_r(i,:,:,:) =kappa(i,:,:,:)+ Kappa/volsum end do deallocate(Vol) deallocate(Volsu) deallocate(Volss) deallocate(Voldu) deallocate(Volbc) deallocate(Volsoa) deallocate(Volpom) deallocate(Volna2so4) deallocate(Volmsa) deallocate(Volno3) deallocate(Volsum) deallocate(over_rad_frac) deallocate(Kappa) deallocate(A_Koehler) deallocate(N_act) deallocate(rc) deallocate(nNa) deallocate(nCl) deallocate(nSO4) deallocate(nNaCl) deallocate(nNa2SO4) deallocate(nH2SO4) end subroutine Kappa_diag end module user_output_general