#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_AERCHEMMIP ! ! !DESCRIPTION: ! isoanus,isoaais,isoaacs,isoacos,isoaaii ! This module provides the code needed to produce the CMIP6 AERCHEMMIP ! output from TM5. Code is based on the user_output_aerocom. ! ! output_aerchemmip_init: ! - initialise the list of variables (allvars) ! - initialise the data holder within each variable (2Ddata/3Ddata,...) ! - initialise the output netcdf files, one for eacht variable ! output_aerchemmip_accumulate: ! - do accumulation for all variables and save data to either ! 2Ddata or 3Ddata holder of the variable (excluding optical vars) ! output_aerchemmip_write ! - write the monthly variable data to netcdf-file ! - initialise the dataholders for accumulation of new month ! output_aerchemmip_write_hourly ! - write the hourly variable data to netcdf-file ! - initialise the dataholders for accumulation of new hour ! output_aerchemmip_write_daily ! - write the daily variable data to netcdf-file ! - initialise the dataholders for accumulation of new day ! write_var ! ! calculate_optics ! - calculate the abss,aods and extinctions ! - directly accumulate the data containers of these variables ! ! mode_fraction ! - calculate the fraction of a M7 mode that is under a size limit !\\ !\\ ! !INTERFACE: ! MODULE USER_OUTPUT_AERCHEMMIP ! ! !USES: ! use go, only : gol, goErr, goPr, goLabel use GO, ONLY : GO_Timer_Def, GO_Timer_End, GO_Timer_Start use dims, only : nregions !=1, support for zooming with larger values not available for AERCHEMMIP use optics, only : wavelendep use meteodata, only : global_lli, levi !use chem_param, only : nmod,xmc2h6,xmc3h6,xmc3h8,xmch4,xmco,xmdms,xmh2o,xmhno3,xmisop,xmno,xmno2,xmo3,xmoh,xmpan,xmso2 use MDF use TM5_DISTGRID, only : dgrid, Get_DistGrid, update_halo, update_halo_iband use chem_param use m7_data, only : h2o_mode implicit none private ! ! !PUBLIC MEMBER FUNCTIONS: ! public :: output_aerchemmip_init !public :: output_aerchemmip_step public :: output_aerchemmip_write public :: output_aerchemmip_write_hourly public :: output_aerchemmip_write_6hourly public :: output_aerchemmip_write_daily public :: output_aerchemmip_done public :: accumulate_data public :: wdep_out character(len=*), parameter :: mname = 'user_output_aerchemmip' ! max indices, at maximum 7, one for each mode integer,parameter :: n_indices=11 type varfile integer :: itm5 ! character(len=16) :: vname ! character(len=64) :: lname ! character(len=11) :: unit ! character(len=10) :: positive ! character(len=130) :: standard_name ! ! real,dimension (:,:),pointer :: dataZonal ! real,dimension (:,:),pointer :: data2D ! real,dimension (:,:,:),pointer :: data3D ! real,dimension (:,:,:),pointer :: budgetdata ! integer :: varid ! integer :: fileunit ! file unit number character(len=200) :: filename ! integer :: dimensions ! integer :: lon_varid ! integer :: lat_varid ! integer :: lev_varid ! integer :: time_varid integer :: hyam_varid integer :: hybm_varid integer :: hyai_varid integer :: hybi_varid integer :: bounds_varid integer :: dims character(len=10) :: freq character(len=9) :: class ! which class of variable :emi, ddep, wdep,conc,aod,met,crescendo integer,dimension(n_indices) :: indices real :: xmgas character(len=20) :: table_id end type varfile type dimdata integer :: nlon ! size of x dimension of requested field integer :: nlat ! size of y dimension of requested field integer :: nlev ! size of 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 dimdata type(dimdata)::dimension_data type budgetstore real, dimension(:,:,:), allocatable :: f2dslast integer :: lasttime end type budgetstore type(budgetstore), dimension(nregions), save :: drydepos, wetdepos, emission ! type experimentdata ! character(len=16) :: ! end type experimentdata ! wavelength information type(wavelendep), dimension(:), allocatable :: wdep_out !!!! integer::test_fileunit !!!! integer :: n_vars=0 type(varfile), dimension(:), allocatable :: allvars type(varfile), dimension(:), allocatable :: fixedvars integer :: n_var_max=300 integer :: n_fixed=3 integer, public :: n_days_in_month character(len=20), public :: aerchemmip_exper ! AeroCom experiment name integer, save :: od550aer, & ec550aer,& ec550aermon,& od550aerdaily, & od550so4, & od550bc, & od550oa,& od550soa,& od550ss,& od550dust,& od440dustday,& od550dust2dday,& od550dust3dday,& od550no3,& od550aerh2o,& od550lt1aer,& abs550aer,& od440aer,& od350aer,& od870aer,& od440aerhr,& od440aermonthly,& od440aerdaily,& od550aerhr,& areacella,& sftlf,& orog integer :: fid ! file id for IF_NOTOK_MDF macro integer :: access_mode integer :: accumulation_mon,accumulation_day,accumulation_hr,accumulation_6hr integer :: timeidx_mon,timeidx_hr,timeidx_day,timeidx_6hr logical,public::crescendo_out=.false. integer,parameter::iemiunit=1 integer,parameter::iddepunit=1 !same dimensions as emi integer,parameter::iwdepunit=1 !same dimensions as emi integer,parameter::iprod3dunit=2 integer,parameter::immrunit=3 integer,parameter::idimensionlessunit=4 integer,parameter::iheightunit=5 integer,parameter::itempunit=6 integer,parameter::io3unit=7 integer,parameter::ipresunit=8 integer,parameter::ivmrunit=9 integer,parameter::irateunit=10 integer,parameter::iloadunit=11 integer,parameter::iextunit=12 integer,parameter::iccunit=13 integer,parameter::im3unit=14 integer,parameter::imassunit=15 character(len=11),dimension(15),parameter::units=(/'kg m-2 s-1','kg m-3 s-1','kg kg-1','1','m','K','DU','Pa','mole mole-1',& 's-1','kg m-2','m-1','cm-3','m-3','kg'/) !CRESCENDO !3D Character(len=11),dimension(40),parameter :: crescendo3D= (/'nd1', 'nd2', 'nd3', 'nd4', 'nd5', 'nd6', 'nd7', & 'mmrsu1', 'mmrsu2', 'mmrsu3', 'mmrsu4', 'mmrsoa1', 'mmrsoa2', 'mmrsoa3', 'mmrsoa4', 'mmrsoa5', 'mmroa2', & 'mmroa3', 'mmroa4', 'mmroa5', 'mmrbc2', 'mmrbc3', 'mmrbc4', 'mmrbc5', 'mmrss3', 'mmrss4', 'mmrdu3', 'mmrdu4', & 'mmrdu6', 'mmrdu7', 'mmraerh2o_1', 'mmraerh2o_2', 'mmraerh2o_3', 'mmraerh2o_4', 'mono', 'nh3', 'ndtot', & 'ec550aer', 'chepsoa3d','emilnox'/) Character(len=11),dimension(40),parameter :: crescendo3Dunit=(/units(im3unit), units(im3unit), units(im3unit), & units(im3unit), units(im3unit), units(im3unit), units(im3unit), units(immrunit), units(immrunit), units(immrunit), & units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit),& units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit), & units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit), & units(immrunit), units(immrunit), units(immrunit), units(ivmrunit), units(ivmrunit), units(im3unit), units(iextunit), & units(iemiunit),units(iemiunit)/) Character(len=11),dimension(2),parameter :: crescendo3Dday=(/'co', 'od5503ddust'/) Character(len=11),dimension(2),parameter :: crescendo3Ddayunit=(/units(ivmrunit),units(idimensionlessunit)/) !2D !mon Character(len=8),dimension(14),parameter :: crescendo2Dmon=(/'drydms', 'wetdms', 'wetno3', 'dryhno3', 'wethno3', & 'dryno2', 'dryno', 'drypan', 'emimono', 'dmsos', 'seddust', 'uas', 'vas', 'sfcWind'/) Character(len=11),dimension(14),parameter :: crescendo2Dmonunit=(/units(iddepunit), units(iddepunit), units(iddepunit), & units(iddepunit), units(iddepunit), units(iddepunit), units(iddepunit), units(iddepunit), units(iddepunit), & 'kg m-3',units(iddepunit),'m s-1', 'm s-1', 'm s-1'/) ! 2d ! 6hr Character(len=16),dimension(19),parameter :: crescendo2D6hr=(/'sfmmrso4', 'sfmmrss', 'sfmmroa', 'sfmmrsoa', 'sfmmrbc', & 'sfmmrdustabv1', 'sfmmrdustabv10', 'sfmmrdustbel1', 'sfdms', 'sfisop', 'sfmono', 'emidms', 'emiss', 'emioa', & 'emiisop', 'emimono', 'sfndtot', 'sfnd100', 'chepsoa2d'/) Character(len=11),dimension(19),parameter :: crescendo2D6hrunit=(/units(immrunit), units(immrunit), units(immrunit),& units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(ivmrunit), & units(ivmrunit), units(ivmrunit), units(iemiunit), units(iemiunit), units(iemiunit), units(iemiunit), & units(iemiunit), units(im3unit), units(im3unit),'kgm-2s-1'/) !2d !day Character(len=16),dimension(52),parameter :: crescendo2Dday=(/'sfnd1', 'sfnd2', 'sfnd3', 'sfnd4', 'sfnd5', 'sfnd6', & 'sfnd7', 'sfmmrsu1', 'sfmmrsu2', 'sfmmrsu3', 'sfmmrsu4', 'sfmmrsoa1', 'sfmmrsoa2', 'sfmmrsoa3', 'sfmmrsoa4', & 'sfmmrsoa5', 'sfmmroa2', 'sfmmroa3', 'sfmmroa4', 'sfmmroa5', 'sfmmrbc2', 'sfmmrbc3', 'sfmmrbc4', 'sfmmrbc5', & 'sfmmrss3', 'sfmmrss4', 'sfmmrdu3', 'sfmmrdu4', 'sfmmrdu6', 'sfmmrdu7', 'sfmmraerh2o_1', 'sfmmraerh2o_2', & 'sfmmraerh2o_3', 'sfmmraerh2o_4', 'sfmmrno3', 'sfmmrnh4', 'sfmmrnh3', 'sfsh', 'od440aer', 'od440dust', & 'od550dust', 'depdustbel1', 'depdustabv1', 'depdustabv10', 'sfmmrdust', 'drynh4', 'wetnh4', 'dryno3', & 'wetno3', 'dryhno3','drydust','wetdust'/) Character(len=11),dimension(52),parameter :: crescendo2Ddayunit=(/units(im3unit), units(im3unit), units(im3unit), & units(im3unit), units(im3unit), units(im3unit), units(im3unit), units(immrunit), units(immrunit), units(immrunit), & units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit),& units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit), & units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit),& units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit),& units(idimensionlessunit), units(idimensionlessunit), units(idimensionlessunit), units(idimensionlessunit), units(iddepunit),& units(iddepunit), units(iddepunit), units(immrunit), units(iddepunit), units(iddepunit), units(iddepunit), units(iddepunit), & units(iddepunit), units(iddepunit), units(iddepunit)/) !2d !hr Character(len=9), dimension(5), parameter :: crescendo2Dhr=(/'od550aer', 'od440aer', 'sfno', 'sfnh3', 'sfhno3'/) Character(len=11), dimension(5), parameter :: crescendo2Dhrunit=(/units(idimensionlessunit), units(idimensionlessunit), & units(ivmrunit), units(ivmrunit), units(ivmrunit)/) !Character(len=11),dimension(6,2),parameter :: crescendo2Dtest=reshape(& ! (/'od550aer', 'od440aer', 'sfno' ,'sfnh3' ,'sfhno3' ,'chepsoa2d' ,& ! '1' ,'1' ,units(ivmrunit), units(ivmrunit), units(ivmrunit),'kgm-2s-1'/),(/6,2/)) !METEO !3D Character(len=7),dimension(11),parameter :: meteo3D=(/'ta', 'pfull', 'phalf', 'hus', 'zg', 'airmass', 'emilnox', & 'jno2', 'ua', 'va', 'wa'/) Character(len=11),dimension(11),parameter :: meteo3Dunit=(/units(itempunit), units(ipresunit), units(ipresunit), & units(idimensionlessunit), units(iheightunit), units(iloadunit), 'mol s-1', units(irateunit),'ms-1', 'ms-1', 'ms-1'/) !2D Character(len=7),dimension(9),parameter :: meteo2D=(/'ps', 'ptp', 'tatp', 'ztp', 'bldep', 'pr', 'tropoz', 'toz', 'albsrfc'/) Character(len=11),dimension(9),parameter :: meteo2Dunit=(/units(ipresunit), units(ipresunit), units(itempunit), & units(iheightunit), units(iheightunit), units(iemiunit), units(io3unit), units(io3unit), units(idimensionlessunit)/) !OPTICS Character(len=11),dimension(13),parameter :: opticscomp=(/'od550aer', 'od550so4', 'od550bc', 'od550oa', 'od550soa',& 'od550ss', 'od550dust', 'od550no3', 'od550aerh2o', 'od550lt1aer', 'od440aer', 'od870aer', 'abs550aer'/)! !AEROSOL compounds Character(len=3),dimension(6),parameter :: comp=(/'BC', 'POM', 'SO4', 'DU', 'SS', 'SOA'/)!SOA !CHEMICAL Character(len=6),dimension(2),parameter :: checomp=(/'aqpso4', 'gpso4'/) Character(len=6),dimension(1),parameter :: chepcomp=(/'soa'/) Character(len=7),dimension(4),parameter :: o3chepcomp=(/'o3loss', 'o3prod','lossch4','lossco'/) !Emon Character(len=9),dimension(8),parameter :: emonout=(/'flashrate', 'depdust','seddustCI','md','loaddust','concdust','conccn','sconcss'/) Character(len=11),dimension(8),parameter :: emonoutunit=(/'km-2 s-1', units(iddepunit),units(iddepunit),'m','kg m-2','kg m-3','m-3','kg m-3'/) !BUDGET (EMI+REMOVAL) Character(len=4),dimension(14),parameter :: emicomp=(/'bvoc', 'co', 'dms', 'isop', 'nox', 'nh3', 'oa', 'so2', 'bc', & 'so4', 'dust', 'ss', 'voc','poa'/) Character(len=4),dimension(12),parameter :: ddepcomp=(/'nh3', 'noy', 'o3', 'oa', 'so2', 'bc', 'so4', 'dust', 'ss', & 'soa', 'nh4', 'no3'/) Character(len=4),dimension(10),parameter :: wdepcomp=(/'bc', 'dust', 'nh3', 'nh4', 'noy', 'oa', 'so2', 'so4', 'ss', 'soa'/) !MIXING ratios Character(len=6),dimension(13),parameter :: aerommrcomp=(/'bc', 'dust', 'nh3', 'nh4', 'no3', 'oa', 'so4', 'ss', 'pm1', & 'pm2p5', 'pm10', 'aerh2o', 'soa'/) Character(len=8),dimension(20),parameter :: gascomp=(/'c2h6', 'c3h6', 'c3h8', 'ch3coch3', 'ch4', 'co', 'co2', 'dms', & 'h2o', 'hno3', 'isop', 'no', 'no2', 'o3', 'oh', 'pan', 'so2', 'voc', 'hcho', 'o3ste'/) !Molecular weights real,dimension(20),parameter :: xmgascomp=(/xmc2h6, xmc3h6, xmc3h8, xmacet, xmch4, xmco, -1.0, xmdms, xmh2o, xmhno3, & xmisop, xmno, xmno2, xmo3, xmoh, xmpan, xmso2, 1.0, xmch2o,xmo3/)!VOC=1, !hcho=ch2o in TM5, but output for hcho is needed. !SPECIAL !HOURLY character(len=8),dimension(1),parameter:: ps6hr=(/'ps'/) character(len=11),dimension(1),parameter:: ps6hrunit=(/units(ipresunit)/) character(len=8),dimension(6),parameter:: hourly_var=(/'ps', 'sfno2', 'sfo3', 'sfpm25', 'tas', 'ec550aer'/) character(len=11),dimension(6),parameter:: hourly_varunit=(/units(ipresunit), units(ivmrunit), units(ivmrunit), & units(immrunit), units(itempunit), units(iextunit)/) !DAILY character(len=8),dimension(10),parameter:: daily_var=(/'od550aer', 'toz', 'maxpblz', 'minpblz', 'tasmin', 'tasmax', & 'pr', 'sfo3max', 'tas','ps'/) character(len=11),dimension(10),parameter:: daily_varunit=(/units(idimensionlessunit), units(io3unit), & units(iheightunit), units(iheightunit), units(itempunit), units(itempunit), units(iemiunit), units(ivmrunit), units(itempunit),units(ipresunit)/) !ZONAL character(len=6),dimension(8),parameter:: zonal_var=(/'ch4', 'hno3', 'ho2', 'noy', 'ta', 'zg', 'o3', 'oh'/) character(len=11),dimension(8),parameter:: zonal_varunit=(/units(ivmrunit), units(ivmrunit), units(ivmrunit), & units(ivmrunit), units(itempunit), units(iheightunit), units(ivmrunit), units(ivmrunit)/) integer,dimension(8),parameter:: zonal_idx=(/ich4,ihno3,iho2,-1,-1,-1,io3,ioh/) !AERCHEMMIP global attributes that might change with run or something else character(len=3),parameter::grid='3x2'!'250 km' character(len=3),parameter::grid_label='gn'!'gnz' for zonal means character(len=300),parameter::source='EC-Earth3-AerChem (2017): atmosphere: IFS cy36r4 (TL255, linearly & &reduced Gaussian grid equivalent to 512 x 256, 91 levels, top level: 0.01 hPa);atmospheric_chemistry: & &TM5 (3 deg. (long.) x 2 deg. (lat.), 34 levels, top level: 0.1 hPa; aerosol: TM5' character(len=17),parameter::source_id='EC-Earth3-AerChem' character(len=20),public ::source_type!='AOGCM CHEM AER' !or 'AGCM CHEM AER' for amip-runs character(len=60),public ::realm character(len=60),public::experiment_id character(len=60),public::experiment character(len=1),public::realization_i='1' character(len=1),public::physics_i='1' character(len=1),public::forcing_i='1' character(len=1),public::initialization_i='1' integer, public:: aerchemmip_dhour ! Timers integer :: itim_init, itim_addvar, itim_write, itim_accu, itim_attr, itim_accu_opt, itim_write_hour, itim_write_day, & itim_write_mon, itim_write_gather contains subroutine output_aerchemmip_init(status) ! Open files ! allocate dataholders use dims, only : newsrun,itau,mlen use global_data, only : outdir use datetime, only : tau2date, date2tau use partools, only : MPI_INFO_NULL, localComm use optics, only : Optics_Init use sedimentation, only : ised,nsed use partools , only : isRoot,myid use global_data, only : region_dat use tm5_distgrid, only : gather use meteodata , only : lsmask_dat,oro_dat use Binas , only : grav implicit none !OUTPUT parameters integer, intent(out) :: status !LOCAL parameters integer :: region !iterator for regions integer :: nlon_region integer :: nlat_region integer :: nlev_region ! also global integer :: j_var !integer :: nlev_region ! also global !integer :: nlev_region ! also global integer :: itrac integer :: i_sed integer :: i,i1,i2,j1,j2,k,j,imr,jmr character(len=*), parameter :: rname = mname//'/output_aerchemmip_init' !FIXED VARS real, dimension(:),pointer :: dxyp real, allocatable :: arr2d(:,:) real ::xmcomp call goLabel(rname) ! define timers: call GO_Timer_Def( itim_init, 'output aerchemmip init', status ) IF_NOTOK_RETURN(status=1) call GO_Timer_Def( itim_write, 'output aerchemmip write', status ) IF_NOTOK_RETURN(status=1) call GO_Timer_Def( itim_write_gather, 'output aerchemmip write gather', status ) IF_NOTOK_RETURN(status=1) call GO_Timer_Def( itim_write_day, 'output aerchemmip write day', status ) IF_NOTOK_RETURN(status=1) call GO_Timer_Def( itim_write_hour, 'output aerchemmip write hour', status ) IF_NOTOK_RETURN(status=1) call GO_Timer_Def( itim_write_mon, 'output aerchemmip write mon', status ) IF_NOTOK_RETURN(status=1) call GO_Timer_Def( itim_accu, 'output aerchemmip accu', status ) IF_NOTOK_RETURN(status=1) call GO_Timer_Def( itim_accu_opt, 'output aerchemmip accu _optics', status ) IF_NOTOK_RETURN(status=1) call GO_Timer_Def( itim_attr, 'output aerchemmip attr', status ) IF_NOTOK_RETURN(status=1) call GO_Timer_Def( itim_addvar, 'output aerchemmip addvar', status ) IF_NOTOK_RETURN(status=1) call Get_DistGrid( dgrid(1), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2) if (newsrun) then !optics? ! 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) !TB Have to keep insitu part, since optics-module usest it for comparisons. allocate( wdep_out( 3 ) ) wdep_out(1)%wl = 0.550 ; wdep_out(1)%split = .true. ; wdep_out(1)%insitu = .false. wdep_out(2)%wl = 0.440 ; wdep_out(2)%split = .true. ; wdep_out(2)%insitu = .false. wdep_out(3)%wl = 0.870 ; wdep_out(3)%split = .false. ; wdep_out(3)%insitu = .false. !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) accumulation_mon=0.0 accumulation_hr=0.0 accumulation_6hr=0.0 accumulation_day=0.0 region=1 ! intermediate structures for budgets allocate(drydepos(region)%f2dslast(i1:i2,j1:j2,18)) allocate(wetdepos(region)%f2dslast(i1:i2,j1:j2,13)) allocate(emission(region)%f2dslast(i1:i2,j1:j2,13)) !imr = i2-i1+1 !jmr = j2-j1+1 ! these here are the initial budgets (monthly): 0.0 drydepos(region)%f2dslast = 0.0 wetdepos(region)%f2dslast = 0.0 emission(region)%f2dslast = 0.0 imr = global_lli(1)%nlon jmr = global_lli(1)%nlat ! for areacella,orog,sftlf if (isRoot) then allocate( arr2d(imr,jmr) ) else allocate( arr2d(1,1) ) endif arr2d(:,:)=0.0 ! for monthly output ! initialise with 31 for january n_days_in_month=31 end if call GO_Timer_Start( itim_init, status ) IF_NOTOK_RETURN(status=1) ! AERCHEMMIP only available for global-> region=1 region=1 !Initialise grid definitions nlon_region = global_lli(region)%nlon nlat_region = global_lli(region)%nlat nlev_region = levi%nlev dimension_data%nlon= nlon_region dimension_data%nlat= nlat_region dimension_data%nlev= nlev_region allocate(dimension_data%lon(nlon_region)) allocate(dimension_data%lat(nlat_region)) allocate(dimension_data%lev(nlev_region)) dimension_data%lon=global_lli(region)%lon_deg dimension_data%lat=global_lli(region)%lat_deg ! initialise output timeidx used for keeping track which output steps is written timeidx_mon=1 timeidx_day=1 timeidx_hr=1 timeidx_6hr=1 ! allocate room for variables allocate(allvars(n_var_max)) allocate(fixedvars(n_fixed)) if (crescendo_out)then do i=1,size(crescendo3D) if (trim(crescendo3D(i))=='mono')then xmcomp=xmterp else if (trim(crescendo3D(i))=='nh3')then xmcomp=xmnh3 else write(gol,*) 'CRESCENDO 3D monthly with negative molar mass' xmcomp=-1.0 end if call add_variable(-1,trim(crescendo3D(i)),trim(crescendo3D(i)),crescendo3Dunit(i),3,status,'crescendo','AERmon',xmcomp) end do do i=1,size(crescendo3Dday) if (trim(crescendo3Dday(i))=='co')then xmcomp=xmco else write(gol,*) 'CRESCENDO 3D daily with negative molar mass' xmcomp=-1.0 end if call add_variable(-1,trim(crescendo3Dday(i)),trim(crescendo3Dday(i)),crescendo3Ddayunit(i),3,status,'crescendo','AERday',xmcomp) end do do i=1,size(crescendo2D6hr) if (trim(crescendo2D6hr(i))=='sfdms')then xmcomp=xmdms else if (trim(crescendo2D6hr(i))=='sfisop')then xmcomp=xmisop else if (trim(crescendo2D6hr(i))=='sfmono')then xmcomp=xmterp else write(gol,*) 'CRESCENDO 2D 6hr with negative molar mass' write(gol,*) 'gascomp with negative molar mass' xmcomp=-1.0 end if call add_variable(-1,trim(crescendo2D6hr(i)),trim(crescendo2D6hr(i)),crescendo2D6hrunit(i),2,status,'crescendo','AER6hr',xmcomp) end do do i=1,size(crescendo2Dmon) call add_variable(-1,trim(crescendo2Dmon(i)),trim(crescendo2Dmon(i)),crescendo2Dmonunit(i),2,status,'crescendo','AERmon',-1.0) end do do i=1,size(crescendo2Dhr) if (trim(crescendo2Dhr(i))=='sfno')then xmcomp=xmno else if (trim(crescendo2Dhr(i))=='sfnh3')then xmcomp=xmnh3 else if (trim(crescendo2Dhr(i))=='sfhno3')then xmcomp=xmhno3 else ! -1 so that missing molar mass will be noticed easily write(gol,*) 'CRESCENDO 2D hr with negative molar mass' xmcomp=-1.0 end if call add_variable(-1,trim(crescendo2Dhr(i)),trim(crescendo2Dhr(i)),crescendo2Dhrunit(i),2,status,'crescendo','AERhr',xmcomp) end do do i=1,size(crescendo2Dday) call add_variable(-1,trim(crescendo2Dday(i)),trim(crescendo2Dday(i)),crescendo2Ddayunit(i),2,status,'crescendo','AERday',-1.0) !call add_variable(-1,trim(crescendo2Dday_new(i,1)),trim(crescendo2Dday_new(i,1)),crescendo2Dday_new(i,2),2,status,'crescendo','AERday',-1.0) end do end if do i=1,size(ps6hr) call add_variable(-1,trim(ps6hr(i)),trim(ps6hr(i)),ps6hrunit(i),2,status,'ps6h','AER6hr',-1.0) end do ! add deposition variables do i=1,size(emicomp) call add_variable(-1,'emi'//trim(emicomp(i)),'emission '//trim(emicomp(i)), units(iemiunit),2,status,'emi','AERmon',-1.0) end do ! add 3D chemical production fields do i=1,size(checomp) call add_variable(-1,'che'//trim(checomp(i)),'chemical production of '//trim(checomp(i)), units(iemiunit),3,status,'chep','AERmon',-1.0) end do ! add 3D ozone prod loss do i=1,size(o3chepcomp) call add_variable(-1,trim(o3chepcomp(i)),'tendency_'//trim(o3chepcomp(i)),'mol m-3 s-1',3,status,'chep','AERmon',-1.0) end do ! add 2D chemical production fields do i=1,size(chepcomp) call add_variable(-1,'chep'//trim(chepcomp(i)),'chemical production of '//trim(chepcomp(i)), units(iemiunit),2,status,'chep','AERmon',-1.0) end do ! add Emon fields do i=1,size(emonout) select case (trim(emonout(i))) case( 'md','concdust','conccn') call add_variable(-1,trim(emonout(i)),trim(emonout(i)), emonoutunit(i),3,status,'emon','Emon',-1.0) case('flashrate','depdust','seddustCI','loaddust','sconcss') call add_variable(-1,trim(emonout(i)),trim(emonout(i)), emonoutunit(i),2,status,'emon','Emon',-1.0) end select end do ! add dry deposition fields do i=1,size(ddepcomp) call add_variable(-1,'dry'//trim(ddepcomp(i)),'dry deposition '//trim(ddepcomp(i)), units(iddepunit),2,status,'dry','AERmon',-1.0) end do ! add wetdep variables do i=1,size(wdepcomp) call add_variable(-1,'wet'//trim(wdepcomp(i)),'wet deposition '//trim(wdepcomp(i)), units(iwdepunit),2,status,'wet','AERmon',-1.0) end do ! add optics fields do i=1,size(opticscomp) !if (trim(opticscomp(i))=='ec550aer') then ! call add_variable(-1,trim(opticscomp(i)),'optics '//trim(opticscomp(i)), units(iextunit),3,status,'optics','AER6hr',-1.0) !else call add_variable(-1,trim(opticscomp(i)),'optics '//trim(opticscomp(i)), units(idimensionlessunit),2,status,'optics','AERmon',-1.0) end do ! Aerosol species mass mixing ratios do i=1,size(aerommrcomp) call add_variable(-1,'mmr'//trim(aerommrcomp(i)),'mass mixing ratio of '//trim(aerommrcomp(i)), units(immrunit),3,status,'mmr','AERmon',-1.0) end do ! Gas-phase species volume mixingratios do i=1,size(gascomp) if (xmgascomp(i)>0.0) then call add_variable(-1,trim(gascomp(i)),'volume mixing ratio of '//trim(gascomp(i)), units(ivmrunit),3,status,'gas','AERmon',xmgascomp(i)) else write(gol,*) 'gascomp with negative molar mass' end if end do ! add meterorological fields do i=1,size(meteo3D) call add_variable(-1,trim(meteo3D(i)),trim(meteo3D(i)),meteo3Dunit(i),3,status,'meteo3D','AERmon',-1.0) end do ! surface meteorological fields do i=1,size(meteo2D) call add_variable(-1,trim(meteo2D(i)),trim(meteo2D(i)),meteo2Dunit(i),2,status,'meteo2D','AERmon',-1.0) end do ! add hourly output do i=1,size(hourly_var) if (trim(hourly_var(i))=='ec550aer')then call add_variable(-1,trim(hourly_var(i)),'optics '//trim(hourly_var(i)), units(iextunit),3,status,'optics','AER6hr',-1.0) else call add_variable(-1,trim(hourly_var(i)),trim(hourly_var(i)),hourly_varunit(i),2,status,'sf1h','AERhr',-1.0) end if end do ! add daily fields do i=1,size(daily_var) call add_variable(-1,trim(daily_var(i)),trim(daily_var(i)),daily_varunit(i),2,status,'sf1d','AERday',-1.0) end do ! add zonal fields, some fields are repeated (3d/zonal) do i=1,size(zonal_var) call add_variable(zonal_idx(i),trim(zonal_var(i)),trim(zonal_var(i)),zonal_varunit(i),3,status,'zonal','AERmonZ',-1.0) end do call add_variable(-1,'areacella','cell area','m2',2,status,'fixed','AERfx',-1.0) call add_variable(-1,'orog','surface_altitude','m',2,status,'fixed','AERfx',-1.0) call add_variable(-1,'sftlf','land_area_fraction','1',2,status,'fixed','AERfx',-1.0) ! do j_var = 1, n_vars ! initialise a single file for each variables as per AERCHEMMIP specification ! overwrite existing files (clobber) if (isroot)call MDF_Create( allvars(j_var)%filename, MDF_NETCDF4, MDF_REPLACE, allvars(j_var)%fileunit, status ) IF_NOTOK_RETURN(status=1) !For each file ! write grid dimension attributes if (isroot) call write_dimensions(status, j_var) IF_NOTOK_RETURN(status=1) ! write global attributes if (isroot) call write_attributes(status, j_var) IF_NOTOK_RETURN(status=1) !write dimension variables if (isroot) call write_var(status,j_var) IF_NOTOK_RETURN(status=1) if (allvars(j_var)%table_id=='AERfx')then if (trim(allvars(j_var)%vname)=='areacella')then ! Gridbox area dxyp => region_dat(1)%dxyp do j=j1,j2 allvars(j_var)%data2D(i1:i2,j)=dxyp(j) end do call gather( dgrid(1), allvars(j_var)%data2D , arr2d(:,:), 0,status) if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%varid, arr2d, status, start=(/1,1/), count=(/imr,jmr/) ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) else if (trim(allvars(j_var)%vname)=='orog')then ! Gridbox area allvars(j_var)%data2D(i1:i2,j1:j2) =oro_dat(region)%data(i1:i2,j1:j2,1)/grav call gather( dgrid(1), allvars(j_var)%data2D , arr2d(:,:), 0,status) if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%varid, arr2d, status, start=(/1,1/), count=(/imr,jmr/) ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) else if (trim(allvars(j_var)%vname)=='sftlf')then ! Gridbox area allvars(j_var)%data2D(i1:i2,j1:j2)=lsmask_dat(1)%data(i1:i2,j1:j2,1)/100. call gather( dgrid(1), allvars(j_var)%data2D , arr2d(:,:), 0,status) if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%varid, arr2d, status, start=(/1,1/), count=(/imr,jmr/) ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) end if end if end do deallocate(arr2d) call GO_Timer_End( itim_init, status ) IF_NOTOK_RETURN(status=1) call goLabel() status = 0 end subroutine output_aerchemmip_init subroutine output_aerchemmip_write(region,newhour,status) use MeteoData, only : temper_dat,sst_dat,albedo_dat,ci_dat,lsp_dat,cp_dat,ssr_dat,str_dat,& blh_dat,gph_dat,lwc_dat,iwc_dat,cco_dat,cc_dat,humid_dat, m_dat,phlb_dat,sp_dat ! use global_data, only : conv_dat use GO, only : TDate, NewDate use go_date,only: days_in_month! use datetime, only : tau2date,date2tau,julday ! use dims, only : itau,iyear0 !current time use ebischeme, only : AC_diag_prod,AC_O3_lp,AC_loss use tm5_distgrid, only : dgrid, Get_DistGrid ,gather use partools , only : isRoot,myid ! use domain_decomp, only: gather implicit none logical,intent(in) ::newhour integer,intent(out)::status integer::region integer:: j_var integer:: imr,jmr,i,i1,i2,j1,j2,lmr character(len=*), parameter :: rname = mname//'/output_aerchemmip_write' real, allocatable :: arr3d(:,:,:),arr3dh(:,:,:),arr2d(:,:) integer, dimension(6) :: curdate ! reference time: integer, parameter :: time_reftime6(6) = (/1750,01,01,00,00,00/) integer(kind=8) :: itauref ! reftime in seconds real :: reftime ! seconds from reference time real :: rangemon type(Tdate)::curdate_tdate call goLabel(rname) call GO_Timer_Start( itim_write_mon, status ) IF_NOTOK_RETURN(status=1) if (region > 1) then write(gol,*) 'output_aerhemmip_write: region >1, only available in global mode!' call goErr status=1 return end if call Get_DistGrid( dgrid(1), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2) ! entire region grid size imr = global_lli(1)%nlon jmr = global_lli(1)%nlat lmr = levi%nlev ! define the reference time in seconds (itauref) ! (for now in days since 1850-01-01 00:00, local variable) ! returns the difference to current time, beginning of new step ! call date2tau( time_reftime6, itauref ) ! calculate reftime as fractional days from itauref, hence division by 86400 ! call date2tau( idater, itaucur ) ! itau is the 1st day of new month at 00:00 so we need to fix the reftime back half a month (15th day) ! ((cursecond - reftimesecond) / seconds_in_day) - days_in_last_month + 15days !reftime = (itau - itauref -n_days_in_month*24*3600 + 15*24*3600) / 86400. reftime = (itau - itauref ) / 86400. !get current date in integers call tau2date(itau, curdate) ! create a TDATE date variable of the previous month (curdate(3)-1) curdate_tdate=newdate(curdate(1),curdate(2),curdate(3)-1,curdate(4),curdate(5),curdate(6),calender='greg') ! get days in month and save for next step n_days_in_month=days_in_month(curdate_tdate) ! change reftime to beginning of last month (the month data is from) reftime=reftime-n_days_in_month !length of the month-1s(in days) for the time_bounds rangemon=n_days_in_month !-(1.0/86400.0) ! allocate 3D and 4D global arrays for gathering data ! only root needs the full array, but it must be allocated in all tasks if (isRoot) then allocate( arr3d(imr,jmr,lmr) ) allocate( arr3dh(imr,jmr,lmr+1) ) allocate( arr2d(imr,jmr) ) else allocate( arr3d(1,1,1) ) allocate( arr3dh(1,1,1) ) allocate( arr2d(1,1) ) endif arr2d(:,:)=0.0 arr3d(:,:,:)=0.0 arr3dh(:,:,:)=0.0 do j_var =1,n_vars ! hourly and daily variables are handled separately if (allvars(j_var)%table_id=='AERhr'.or.allvars(j_var)%table_id=='AER6hr'.or.& allvars(j_var)%table_id=='AERday'.or.allvars(j_var)%table_id=='AERfx')then cycle end if if (allvars(j_var)%dims==2)then !2D if (trim(allvars(j_var)%vname)/='minpblz'.and.trim(allvars(j_var)%vname)/='tasmin'.and. & trim(allvars(j_var)%vname)/='maxpblz'.and.trim(allvars(j_var)%vname)/='tasmax')then allvars(j_var)%data2D(i1:i2,j1:j2)=allvars(j_var)%data2D(i1:i2,j1:j2)/real(accumulation_mon) end if call GO_Timer_Start( itim_write_gather, status ) IF_NOTOK_RETURN(status=1) call gather( dgrid(1), allvars(j_var)%data2D , arr2d(:,:), 0,status) call GO_Timer_End( itim_write_gather, status ) IF_NOTOK_RETURN(status=1) IF_NOTOK_RETURN(status=1) if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%varid, arr2d, status, start=(/1,1,timeidx_mon/), & count=(/imr,jmr,1/) ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) ! Zero out the accumulated and written data for the next interval if (trim(allvars(j_var)%vname)=='minpblz' .or. trim(allvars(j_var)%vname)=='tasmin') then ! put high number so simple comparison is needed for finding minimum allvars(j_var)%data2D(i1:i2,j1:j2)=1e10 else allvars(j_var)%data2D(i1:i2,j1:j2)=0.0 end if else !3D if (trim( allvars(j_var)%vname)=='phalf') then !lmr+1 allvars(j_var)%data3D(i1:i2,j1:j2,1:lmr+1)=allvars(j_var)%data3D(i1:i2,j1:j2,1:lmr+1)/real(accumulation_mon) call GO_Timer_Start( itim_write_gather, status ) IF_NOTOK_RETURN(status=1) call gather( dgrid(1), allvars(j_var)%data3D , arr3dh(:,:,:), 0, status) IF_NOTOK_RETURN(status=1) call GO_Timer_End( itim_write_gather, status ) IF_NOTOK_RETURN(status=1) if (isroot) call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%varid, arr3dh , status, start=(/1,1,1,timeidx_mon/), & count=(/imr,jmr,lmr+1,1/) ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) else allvars(j_var)%data3D(i1:i2,j1:j2,1:lmr)=allvars(j_var)%data3D(i1:i2,j1:j2,1:lmr)/real(accumulation_mon) call GO_Timer_Start( itim_write_gather, status ) IF_NOTOK_RETURN(status=1) call gather( dgrid(1), allvars(j_var)%data3D , arr3d(:,:,:), 0, status) IF_NOTOK_RETURN(status=1) call GO_Timer_End( itim_write_gather, status ) IF_NOTOK_RETURN(status=1) if (isroot) call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%varid, arr3d , status, start=(/1,1,1,timeidx_mon/), & count=(/imr,jmr,lmr,1/) ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) end if ! Zero out the accumulated and written data for the next interval allvars(j_var)%data3D(i1:i2,j1:j2,:)=0.0 end if !end if ! write the date for timestep if (isroot) call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%time_varid,(/ reftime+real(rangemon/2)/) , status, start=(/timeidx_mon/), count=(/1/) ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) if (isroot) call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%bounds_varid,(/ reftime,reftime+rangemon/) , status, & start=(/1,timeidx_mon/), count=(/2,1/) ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) end do deallocate( arr3d,arr3dh,arr2d) ! change time index timeidx_mon= timeidx_mon + 1 ! accululated time to zero accumulation_mon=0 ! zero out the chemical production (for mongthly output) !AC_diag_prod(region)%prod(i1:i2,j1:j2,:,1:3)=0.0 ! zero out the chemical production !AC_O3_lp(region)%prod(i1:i2,j1:j2,:,1:2)=0.0 ! zero out the chemical production !AC_loss(region)%prod(i1:i2,j1:j2,:,1:2)=0.0 !status=1 !return call GO_Timer_End( itim_write_mon, status ) IF_NOTOK_RETURN(status=1) call goLabel() status = 0 end subroutine output_aerchemmip_write subroutine output_aerchemmip_write_daily(region,newday,status) use MeteoData, only : temper_dat, sst_dat, albedo_dat, ci_dat, lsp_dat, cp_dat, ssr_dat, str_dat, & blh_dat, gph_dat, lwc_dat, iwc_dat, cco_dat, cc_dat, humid_dat, m_dat, phlb_dat, sp_dat ! use meteodata , only : global_lli, levi use partools , only : isRoot,myid use GO, only : TDate, NewDate! use datetime, only : tau2date,date2tau,julday ! use dims, only : itau,iyear0 !current time use tm5_distgrid, only : dgrid, Get_DistGrid ,gather implicit none logical,intent(in) ::newday integer,intent(out)::status integer::region integer:: imr,jmr,i,i1,i2,j1,j2,lmr character(len=*), parameter :: rname = mname//'/output_aerchemmip_write' integer:: j_var ! reference time: integer, parameter :: time_reftime6(6) = (/1750,01,01,00,00,00/) integer(kind=8) :: itauref ! reftime in seconds real :: reftime ! seconds from reference time real :: rangeday ! for bounds ! root writes netcdf arrays real, allocatable :: arr3d(:,:,:), arr2d(:,:) integer:: imr_f,jmr_f,lmr_f call goLabel(rname) call GO_Timer_Start( itim_write_day, status ) IF_NOTOK_RETURN(status=1) if (region > 1) then write(gol,*) 'output_aerhemmip_write: region >1, only available in global mode!' call goErr status=1 return end if call Get_DistGrid( dgrid(1), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2) ! entire region grid size imr = global_lli(1)%nlon jmr = global_lli(1)%nlat lmr = levi%nlev ! allocate 3D and 4D global arrays for gathering data if (isRoot) then allocate( arr3d(imr,jmr,lmr) ) allocate( arr2d(imr,jmr) ) else allocate( arr3d(1,1,1) ) allocate( arr2d(1,1) ) endif arr2d(:,:)=0.0 arr3d(:,:,:)=0.0 ! define the reference time in seconds (itauref) ! (for now in days since 1850-01-01 00:00, local variable) call date2tau( time_reftime6, itauref ) ! calculate reftime as fractional days from itauref, hence division by 86400 ! call date2tau( idater, itaucur ) reftime = (itau - itauref) / 86400. - 1.0 !23h59 as days rangeday=1.0!(23.0*3600.0+59.0*60.0+59.0)/86400.0 do j_var =1,n_vars if (allvars(j_var)%table_id/='AERday')then cycle end if if (allvars(j_var)%dims==2)then if ( trim(allvars(j_var)%vname)/='minpblz' .and. trim(allvars(j_var)%vname)/='tasmin'.and.trim(allvars(j_var)%vname)/='maxpblz'.and. trim(allvars(j_var)%vname)/='tasmax')then allvars(j_var)%data2D(i1:i2,j1:j2)=allvars(j_var)%data2D(i1:i2,j1:j2)/real(accumulation_day) end if call GO_Timer_Start( itim_write_gather, status ) IF_NOTOK_RETURN(status=1) call gather( dgrid(1), allvars(j_var)%data2D , arr2d(:,:), 0, status) IF_NOTOK_RETURN(status=1) call GO_Timer_End( itim_write_gather, status ) IF_NOTOK_RETURN(status=1) if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%varid, arr2d, status, start=(/1,1,timeidx_day/), count=(/imr,jmr,1/) ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) if (trim(allvars(j_var)%vname)=='minpblz' .or. trim(allvars(j_var)%vname)=='tasmin') then ! put high number so simple comparison is needed for finding minimum allvars(j_var)%data2D(i1:i2,j1:j2)=1e10 else ! Zero out the accumulated and written data for the next interval allvars(j_var)%data2D(i1:i2,j1:j2)=0.0 end if else allvars(j_var)%data3D(i1:i2,j1:j2,1:lmr)=allvars(j_var)%data3D(i1:i2,j1:j2,1:lmr)/real(accumulation_day) !end if call GO_Timer_Start( itim_write_gather, status ) IF_NOTOK_RETURN(status=1) call gather( dgrid(1), allvars(j_var)%data3D , arr3d(:,:,:), 0, status) call GO_Timer_End( itim_write_gather, status ) IF_NOTOK_RETURN(status=1) !if (trim(allvars(j_var)%vname)=='od5503ddust')then IF_NOTOK_RETURN(status=1) if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%varid, arr3d, status, start=(/1,1,1,timeidx_day/), count=(/imr,jmr,lmr,1/) ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) allvars(j_var)%data3D(i1:i2,j1:j2,1:lmr)=0.0 end if ! write the date for timestep if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%time_varid,(/ reftime+0.5/) , status, start=(/timeidx_day/), count=(/1/) ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) if (isroot) call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%bounds_varid,(/ reftime,reftime+ rangeday/) , status, start=(/1,timeidx_day/), count=(/2,1/) ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) end do deallocate(arr3d, arr2d) ! Timeindex to next day timeidx_day= timeidx_day + 1 ! daily accumulated time to zero accumulation_day=0.0 !status=1 !return call GO_Timer_End( itim_write_day, status ) IF_NOTOK_RETURN(status=1) call goLabel() status = 0 end subroutine output_aerchemmip_write_daily subroutine output_aerchemmip_write_hourly(region,newhour,status) use MeteoData, only : temper_dat,sst_dat,albedo_dat,ci_dat,lsp_dat,cp_dat,ssr_dat,str_dat,blh_dat,gph_dat,lwc_dat,iwc_dat,cco_dat,cc_dat,humid_dat, m_dat,phlb_dat,sp_dat ! use GO, only : TDate, NewDate! use datetime, only : tau2date,date2tau,julday ! use dims, only : itau,iyear0 !current time use tm5_distgrid, only : dgrid, Get_DistGrid ,gather use partools , only : isRoot,myid implicit none logical,intent(in) ::newhour integer,intent(out)::status integer:: j_var integer::region integer:: imr,jmr,i,i1,i2,j1,j2,lmr character(len=*), parameter :: rname = mname//'/output_aerchemmip_write' real :: rangehr,range6hr ! hour in days for bounds in NC-file ! reference time: integer, parameter :: time_reftime6(6) = (/1750,01,01,00,00,00/) integer(kind=8) :: itauref ! reftime in seconds real :: reftime ! seconds from reference time ! root writes netcdf arrays real, allocatable :: arr3d(:,:,:) , arr2d(:,:) call goLabel(rname) call GO_Timer_Start( itim_write_hour, status ) IF_NOTOK_RETURN(status=1) if (region > 1) then write(gol,*) 'output_aerhemmip_write: region >1, only available in global mode!' call goErr status=1 return end if call Get_DistGrid( dgrid(1), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2) ! entire region grid size imr = global_lli(1)%nlon jmr = global_lli(1)%nlat lmr = levi%nlev ! allocate 3D and 4D global arrays for gathering data if (isRoot) then allocate( arr3d(imr,jmr,lmr) ) allocate( arr2d(imr,jmr) ) else ! other than root need the variable, but no space allocate( arr3d(1,1,1) ) allocate( arr2d(1,1) ) endif arr2d(:,:)=0.0 arr3d(:,:,:)=0.0 ! define the reference time in seconds (itauref) ! (for now in days since 1850-01-01 00:00, local variable) call date2tau( time_reftime6, itauref ) ! call date2tau( idater, itaucur ) rangehr=1.0/24.0!(3600)/86400.0 do j_var =1,n_vars if (allvars(j_var)%table_id/='AERhr' .and. allvars(j_var)%table_id/='AER6hr' )then cycle end if if (allvars(j_var)%dims==2)then if ( trim(allvars(j_var)%table_id)=='AERhr') then reftime = (itau - itauref) / 86400. - (1./24) allvars(j_var)%data2D(i1:i2,j1:j2)=allvars(j_var)%data2D(i1:i2,j1:j2)/real(accumulation_hr) call GO_Timer_Start( itim_write_gather, status ) IF_NOTOK_RETURN(status=1) call gather( dgrid(1), allvars(j_var)%data2D , arr3d(:,:,1), 0, status) IF_NOTOK_RETURN(status=1) call GO_Timer_End( itim_write_gather, status ) IF_NOTOK_RETURN(status=1) if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%varid,arr3d(:,:,1), status, start=(/1,1,timeidx_hr/), count=(/imr,jmr,1/) ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) ! write the date for timestep if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%time_varid,(/ reftime+(0.5/24.0)/) , status, start=(/timeidx_hr/), count=(/1/) ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) if (isroot) call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%bounds_varid,(/ reftime,reftime+rangehr/) , status, start=(/1,timeidx_hr/), count=(/2,1/) ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) ! Zero out the accumulated and written data for the next interval allvars(j_var)%data2D(i1:i2,j1:j2)=0.0 end if else if ( trim(allvars(j_var)%table_id)=='AERhr') then reftime = (itau - itauref) / 86400. - (1./24) allvars(j_var)%data3D(i1:i2,j1:j2,1:lmr)=allvars(j_var)%data3D(i1:i2,j1:j2,1:lmr)/real(accumulation_hr) call GO_Timer_Start( itim_write_gather, status ) IF_NOTOK_RETURN(status=1) call gather( dgrid(1), allvars(j_var)%data3D , arr3d(:,:,:), 0, status) call GO_Timer_End( itim_write_gather, status ) IF_NOTOK_RETURN(status=1) if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%varid,arr3d, status, start=(/1,1,1,timeidx_hr/), count=(/imr,jmr,lmr,1/) ) ! write the date for timestep if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%time_varid,(/ reftime+(0.5/24.0)/) , status, start=(/timeidx_hr/), count=(/1/) ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) if (isroot) call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%bounds_varid,(/ reftime,reftime+(1./24.)/) , status, start=(/1,timeidx_hr/), count=(/2,1/) ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) ! Zero out the accumulated and written data for the next interval allvars(j_var)%data3D(i1:i2,j1:j2,1:lmr)=0.0 end if end if !end if end do deallocate(arr3d, arr2d) ! changed index to next hour timeidx_hr= timeidx_hr + 1 ! accumulated hours to zero, actually this will always be 1h accumulation_hr=0.0 !status=1 !return call GO_Timer_End( itim_write_hour, status ) IF_NOTOK_RETURN(status=1) call goLabel() status = 0 end subroutine output_aerchemmip_write_hourly subroutine output_aerchemmip_write_6hourly(region,newhour,status) use MeteoData, only : temper_dat,sst_dat,albedo_dat,ci_dat,lsp_dat,cp_dat,ssr_dat,str_dat,blh_dat,gph_dat,lwc_dat,iwc_dat,cco_dat,cc_dat,humid_dat, m_dat,phlb_dat,sp_dat ! use GO, only : TDate, NewDate! use datetime, only : tau2date,date2tau,julday ! use dims, only : itau,iyear0 !current time use tm5_distgrid, only : dgrid, Get_DistGrid ,gather use partools , only : isRoot,myid use ebischeme, only : AC_diag_prod,iprod_soa2dhour implicit none logical,intent(in) ::newhour integer,intent(out)::status integer::region integer:: j_var integer:: imr,jmr,i,i1,i2,j1,j2,lmr character(len=*), parameter :: rname = mname//'/output_aerchemmip_write' ! reference time: integer, parameter :: time_reftime6(6) = (/1750,01,01,00,00,00/) integer(kind=8) :: itauref ! reftime in seconds real :: reftime ! seconds from reference time ! root writes netcdf arrays real, allocatable :: arr3d(:,:,:) , arr2d(:,:) call goLabel(rname) call GO_Timer_Start( itim_write_hour, status ) IF_NOTOK_RETURN(status=1) if (region > 1) then write(gol,*) 'output_aerhemmip_write: region >1, only available in global mode!' call goErr status=1 return end if call Get_DistGrid( dgrid(1), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2) ! entire region grid size imr = global_lli(1)%nlon jmr = global_lli(1)%nlat lmr = levi%nlev ! allocate 3D and 4D global arrays for gathering data if (isRoot) then allocate( arr3d(imr,jmr,lmr) ) allocate( arr2d(imr,jmr) ) else ! other than root need the variable, but no space allocate( arr3d(1,1,1) ) allocate( arr2d(1,1) ) endif arr2d(:,:)=0.0 arr3d(:,:,:)=0.0 ! define the reference time in seconds (itauref) ! (for now in days since 1850-01-01 00:00, local variable) call date2tau( time_reftime6, itauref ) ! call date2tau( idater, itaucur ) reftime = (itau - itauref) / 86400. do j_var =1,n_vars if ( allvars(j_var)%table_id/='AER6hr' )then cycle end if if (allvars(j_var)%dims==2)then !allvars(j_var)%data2D(i1:i2,j1:j2)=allvars(j_var)%data2D(i1:i2,j1:j2)/real(accumulation_6hr) call GO_Timer_Start( itim_write_gather, status ) IF_NOTOK_RETURN(status=1) call gather( dgrid(1), allvars(j_var)%data2D , arr2d(:,:), 0, status) IF_NOTOK_RETURN(status=1) call GO_Timer_End( itim_write_gather, status ) IF_NOTOK_RETURN(status=1) if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%varid,arr2d, status, start=(/1,1,timeidx_6hr/), count=(/imr,jmr,1/) ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) ! write the date for timestep if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%time_varid,(/ reftime/) , status, start=(/timeidx_6hr/), count=(/1/) ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) ! Zero out the accumulated and written data for the next interval allvars(j_var)%data2D(i1:i2,j1:j2)=0.0 else !allvars(j_var)%data3D(i1:i2,j1:j2,1:lmr)=allvars(j_var)%data3D(i1:i2,j1:j2,1:lmr)/real(accumulation_6hr) call GO_Timer_Start( itim_write_gather, status ) IF_NOTOK_RETURN(status=1) call gather( dgrid(1), allvars(j_var)%data3D , arr3d(:,:,:), 0, status) call GO_Timer_End( itim_write_gather, status ) IF_NOTOK_RETURN(status=1) if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%varid,arr3d, status, start=(/1,1,1,timeidx_6hr/), count=(/imr,jmr,lmr,1/) ) ! write the date for timestep if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%time_varid,(/ reftime/) , status, start=(/timeidx_6hr/), count=(/1/) ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) !start=(/i1,j1,1,timeidx_mon/), count=(/imr,jmr,lmr,1/) ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) ! Zero out the accumulated and written data for the next interval allvars(j_var)%data3D(i1:i2,j1:j2,1:lmr)=0.0 end if end do deallocate(arr3d,arr2d) ! changed index to next 6hour timeidx_6hr= timeidx_6hr + 1 ! exception for one 6hr field, no need to do another subroutine for it accumulation_6hr=0.0 ! zero out the chemical production (SOA for hourly output) ! now used fro 6h output !AC_diag_prod(region)%prod(i1:i2,j1:j2,:,iprod_soa2dhour)=0.0 call GO_Timer_End( itim_write_hour, status ) IF_NOTOK_RETURN(status=1) call goLabel() status = 0 end subroutine output_aerchemmip_write_6hourly subroutine accumulate_data(dhour,l_do_ec550aer_only,status) use GO , only : TDate, NewDate, rTotal, operator(-) use Grid , only : FPressure,HPressure use binas, only : rgas, rol,xmair,Dobs,Avog USE toolbox, only : ltropo_ifs, lvlpress !use datetime, only : tau2date use MeteoData, only : temper_dat, sst_dat, albedo_dat, ci_dat, lsp_dat, cp_dat, ssr_dat, str_dat, blh_dat, & gph_dat, lwc_dat, iwc_dat, cco_dat, cc_dat, humid_dat, m_dat, phlb_dat, sp_dat, pu_dat, pv_dat,pw_dat use photolysis_data,only:phot_dat ! use global_data, only : mass_dat, region_dat,conv_dat use tracer_data, only : chem_dat use dims, only : lm,sec_month use chem_param, only : iso4acs, iso4cos, iduacs, iso4ais, ibccos, ibcaii, xmair, xmo3,nmhc,xmcb5,ncb5, o3p, o3l,ino3_a,xmn,ra !use chem_param, only : iso4nus, isscos, ino3_a, issacs, iducos, iduaci, nmod !use chem_param, only : iducoi, ibcacs, ipomcos, ipomaii, ibcais, ipomacs, ipomais !use chem_param, only : imsa, inh4 !use chem_param, only : ntrace,names,mode_end_so4 !use mo_aero_m7, only : nmod!,nsol !use optics, only : optics_aop_get use m7_data, only : h2o_mode, rw_mode, rwd_mode use sedimentation, only : nsed,ised ,sindex use sedimentation, only : deposition => buddep_m7_dat !(i,j,lev,nsed) use wet_deposition,only : wetdep=>buddep_dat !(i,j,lev,ntrace) use emission_data, only : budemi_dat use ebischeme, only : buddry_dat => buddep_dat use chem_param, only : xmso2, xmso4, xmdms, xmpom, xmbc, xmdust, xmnacl,xmo3,xmnh4,ino,ino2,ino3,ihno3,ino3_a,ihno4,ihono,ich3o2no2,ipan,iorgntr,ra use calc_pm, only : PMx_Integrate_3d,PMx_integrate_0d use emission_nox, only : eminox_lightning use ebischeme, only : diag_prod,AC_diag_prod,AC_O3_lp,AC_loss,iloss_ch4,iloss_co,iprod_gasso4,iprod_aqso4,iprod_soa3dmon,iprod_soa2dhour use partools , only : isRoot,myid use dims, only: gtor, dx, dy, ybeg, xref, yref,ndyn use binas, only: ae use emission_nox, only: flashrate_lightning implicit none character(len=*), parameter :: rname = mname//'/output_aerchemmip_accumulate_data' ! integer :: indices(7) integer :: itrac,gasind integer :: i_sed integer :: i_emi integer :: i, j, k, imr, jmr, lmr, lwl, dtime,index,imode,m integer :: i1, i2, j1, j2 integer :: ie, je ! indices for subdomain extended with halo cells integer,intent(in) :: dhour integer :: status integer :: j_var,region,i_var,i_wdep,sedindex,icomp real :: dens real :: vol real :: tempbud,xmcomp,temp real, dimension(:,:,:,:), pointer :: tracers ! transported tracers real, dimension(:,:,:,:), pointer :: tracers_c ! non-transported real, dimension(:), pointer :: dxyp integer, dimension(n_indices) :: indices real::xmgas real, dimension(:,:,:), allocatable :: pres3d real, dimension(:,:,:), allocatable :: temp_pm real, dimension(:,:,:), allocatable :: pres3dh real :: pm_sizelimit integer::emi_index,wet_index,dry_index ! tropopause calculations real, dimension(:), allocatable :: gphx, tx real, dimension(:,:,:), pointer :: gph ! height (incl. oro) real, dimension(:,:,:), pointer :: t ! temperature (K) integer :: itropo real::xres,yres,dxx,dyy,uwind,vwind,lat,wwind,dz,meanwind integer, dimension(4) :: modes_dust=(/3,4,6,7/) ! !EC550AER logical,intent(in) :: l_do_ec550aer_only call goLabel(rname) call GO_Timer_Start( itim_accu, status ) IF_NOTOK_RETURN(status=1) region=1 ! 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 ! for tropopause variables allocate(gphx(0:lm(region))) ! note now from 0-->lm allocate(tx(lm(region))) t => temper_dat (region)%data gph => gph_dat (region)%data allocate( temp_pm(i1:i2,j1:j2,lmr) ) temp_pm=0.0 allocate( pres3d(i1:i2,j1:j2,lmr) ) allocate( pres3dh(i1:i2,j1:j2,0:lmr) ) ! fill mid level pressure call FPressure( levi, sp_dat(region)%data(i1:i2,j1:j2,1), pres3d, status ) IF_NOTOK_RETURN(status=1) ! fill interface pressure call HPressure( levi, sp_dat(region)%data(i1:i2,j1:j2,1), pres3dh, status ) IF_NOTOK_RETURN(status=1) accumulation_6hr=0.0!accumulation_6hr+dtime ! Gridbox area dxyp => region_dat(region)%dxyp ! mass_dat and chem_data keep data in kg/gridbox ! tracers => mass_dat(region)%rm tracers_c => chem_dat(region)%rm if (.not. l_do_ec550aer_only) then dtime = dhour*3600 accumulation_mon=accumulation_mon+dtime accumulation_hr=accumulation_hr+dtime accumulation_day=accumulation_day+dtime do j_var = 1, n_vars indices(:)=allvars(j_var)%indices(:) !Here we use the tm5-indices to collect data for output ! if (trim(allvars(j_var)%class)=='meteo2D') then do j = j1,j2 do i=i1,i2 tx(:)=t(i,j,:) gphx(0:lm(region))=gph(i,j,1:lm(region)+1) !note the bounds !ibase = lubottom(i,j) !itop = lutop(i,j) itropo = ltropo_ifs(region,i,j,tx,lm(region)) ! 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) if (trim(allvars(j_var)%vname)=='ps')then allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*sp_dat(1)%data(i,j,1) else if (trim(allvars(j_var)%vname)=='ptp')then allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+pres3d(i,j,itropo)*dtime else if (trim(allvars(j_var)%vname)=='tropoz')then ! conversion (in order of execution): !kg->g (1e3), g->molec (xmo3), m2->cm2(1e-4) , molec-> molec/m2(dxyp) , molec/cm2->dobson (DOBS) allvars(j_var)%data2D(i,j)=allvars(j_var)%data2D(i,j)+dtime*sum(tracers(i,j,1:itropo,io3)*1e3/xmo3*Avog*1e-4/dxyp(j)/Dobs) else if (trim(allvars(j_var)%vname)=='toz')then ! conversion (in order of execution): !kg->g (1e3), g->molec (xmo3), m2->cm2(1e-4) , molec-> molec/m2(dxyp) , molec/cm2->dobson (DOBS) allvars(j_var)%data2D(i,j)=allvars(j_var)%data2D(i,j)+dtime*sum(tracers(i,j,:,io3)*1e3/xmo3*Avog*1e-4/dxyp(j)/Dobs) else if (trim(allvars(j_var)%vname)=='tos')then allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)!+dtime*sst_dat(1)%data(i,j,1) else if (trim(allvars(j_var)%vname)=='sic')then allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*ci_dat(1)%data(i,j,1) else if (trim(allvars(j_var)%vname)=='tatp')then ! t at the gridpoint center ->mean at the interface value allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*t(i,j,itropo) else if (trim(allvars(j_var)%vname)=='ztp')then !gph at interface allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(gph(i,j,itropo+1)+gph(i,j,itropo))/2 else if (trim(allvars(j_var)%vname)=='bldep')then allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*conv_dat(1)%blh(i,j) else if (trim(allvars(j_var)%vname)=='pr')then allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(lsp_dat(1)%data(i,j,1)+ cp_dat(1)%data(i,j,1)) else if (trim(allvars(j_var)%vname)=='albsrfc')then allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*albedo_dat(1)%data(i,j,1) end if end do end do !end if else if (trim(allvars(j_var)%class)=='meteo3D') then do j = j1,j2 do i=i1,i2 if (trim(allvars(j_var)%vname)=='phalf')then do k=1,lmr+1 allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*pres3dh(i,j,k-1) end do else do k =1,lmr if(trim(allvars(j_var)%vname)=='ta')then allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*temper_dat(1)%data(i,j,k) else if (trim(allvars(j_var)%vname)=='emilnox')then allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime* eminox_lightning(region)%d3(i,j,k) / dxyp(j) else if (trim(allvars(j_var)%vname)=='jno2')then allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*phot_dat(region)%jno2(i,j,k) else if (trim(allvars(j_var)%vname)=='airmass')then allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime* m_dat(region)%data(i,j,k) / dxyp(j) else if (trim(allvars(j_var)%vname)=='pfull')then allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*pres3d(i,j,k) ! else if (trim(allvars(j_var)%vname)=='phalf')then ! allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*pres3dh(i,j,k) else if (trim(allvars(j_var)%vname)=='hus')then allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*humid_dat(1)%data(i,j,k) else if (trim(allvars(j_var)%vname)=='zg')then allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*gph_dat(1)%data(i,j,k) else if (trim(allvars(j_var)%vname)=='ua')then ! distance of single gridbox yres = dy/yref(1) xres = dx/xref(1) lat = ybeg(1) + 0.5 * yres + yres * (j-1) dxx = ae * xres * gtor * cos(lat*gtor) uwind=dxx*(pu_dat(1)%data(i,j,k) + pu_dat(1)%data(i-1,j,k))*0.5 / m_dat(1)%data(i,j,k) allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*uwind else if (trim(allvars(j_var)%vname)=='va')then yres = dy/yref(1) dyy = ae * yres * gtor vwind= dyy *(pv_dat(1)%data(i,j,k) + pv_dat(1)%data(i,j+1,k))*0.5 / m_dat(1)%data(i,j,k) allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*vwind else if (trim(allvars(j_var)%vname)=='wa')then dz = gph_dat(1)%data(i,j,k+1) - gph_dat(1)%data(i,j,k) wwind= dz*(pw_dat(1)%data1(i,j,k-1) + pw_dat(1)%data1(i,j,k))*0.5 / m_dat(1)%data(i,j,k) allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*wwind end if end do end if end do end do !end if else if (trim(allvars(j_var)%class)=='emon') then !Sedimented indices in deposition%sed !do i_sed=1,nsed !gridpoints do j = j1,j2 do i=i1,i2 select case ( trim(allvars(j_var)%vname)) case('loaddust') do m=1,n_indices index=allvars(j_var)%indices(m) if (index==0) cycle do k=1,lmr allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tracers(i,j,k,index))/ dxyp(j) end do end do case('concdust') do m=1,n_indices index=allvars(j_var)%indices(m) if (index==0) cycle do k=1,lmr vol = (gph_dat(region)%data(i,j,k+1)-gph_dat(region)%data(i,j,k)) * dxyp(j) allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*(tracers(i,j,k,index))/ vol end do end do case('seddustCI') sedindex=sindex(iducoi) xmcomp=xmdust allvars(j_var)%data2D(i,j)=allvars(j_var)%data2D(i,j)+ & deposition(region)%sed(i,j,1,sedindex) / dxyp(j) * 1.E-03*xmcomp*-100000 case('depdust') tempbud=0.0 xmcomp=xmdust do m=1,n_indices !exit if indices ==0 !since indices after 0 will all be also 0 if (indices(m)>0) then if (allvars(j_var)%indices(m)>69)then sedindex=-1 else sedindex=sindex(allvars(j_var)%indices(m)) end if if (sedindex>0)then tempbud=tempbud+( & sum(deposition(region)%sed(i,j,:,sedindex))) / dxyp(j) * 1.E-03*xmcomp end if end if allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tempbud) end do case('md') do k=1,lmr ! tendency of atm_mass content of dust dry aerosol particles due to emission ! coarse mode only 3d hmmm allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)-1!dtime*rw_mode(1,7)%d3(i,j,k) end do case('flashrate') allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*flashrate_lightning(1)%d23(i,j)/(dxyp(j)*1e-4) !km-2s-1 (flashrate_lightning in [s-1] case('conccn') do k=1,lmr temp=0.0 dens = pres3d(i,j,k) / temper_dat(region)%data(i,j,k) * xmair * 1.E-3 / (m_dat(region)%data(i,j,k) * rgas) do m=1,n_indices index=allvars(j_var)%indices(m) if (index==0) cycle temp=tracers(i,j,k,index)*dens end do allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*temp end do case('sconcss') temp=0.0 k=1 !Surface dens = pres3d(i,j,k) / temper_dat(region)%data(i,j,k) * xmair * 1.E-3 / (m_dat(region)%data(i,j,k) * rgas) do m=1,n_indices index=allvars(j_var)%indices(m) if (index==0) cycle temp=tracers(i,j,k,index)*dens end do allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*temp case DEFAULT !('default') write(gol,'("emon accumulate missing: ",a)') trim(allvars(j_var)%vname); call goPr end select end do end do !end if else if (trim(allvars(j_var)%class)=='dry') then !Sedimented indices in deposition%sed !do i_sed=1,nsed !gridpoints do j = j1,j2 do i=i1,i2 select case ( trim(allvars(j_var)%vname)) case('dryoa') dry_index=1 xmcomp=xmpom case('drybc') dry_index=2 xmcomp=xmbc case ('dryso2') dry_index=3 xmcomp=xmso2 case('dryso4') dry_index=4 xmcomp=xmso4 case('drydust') dry_index=5 xmcomp=xmdust case('drydms') dry_index=6 xmcomp=xmdms case('dryss') dry_index=7 xmcomp=xmnacl case('drysoa') dry_index=8 xmcomp=xmpom!soa case('drynh3') dry_index=9 xmcomp=xmnh3 case('drynh4') dry_index=10 xmcomp=xmnh4 case('drynoy') dry_index=11 xmcomp=1 case('dryo3') dry_index=12 xmcomp=xmo3 case('dryno3') xmcomp=xmno3 case('dryno2') xmcomp=xmno2 case('dryno') xmcomp=xmno case DEFAULT !('default') write(gol,'("dry xm-var missing: ",a)') trim(allvars(j_var)%vname); call goPr end select tempbud=0.0 !if (xx='dryoa') then !end if do m=1,n_indices !exit if indices ==0 !since indices after 0 will all be also 0 if (indices(m)<=0) then cycle else if (indices(m)>69) then sedindex=-1 else sedindex=sindex(allvars(j_var)%indices(m)) end if !only M7 aerosol tracers (and msa/nh4/no3) are sedimented and deposited if ( trim(allvars(j_var)%vname)=='drynoy')then select case(allvars(j_var)%indices(m)) case(ino,ino2,ino3,ihno3,ino3_a,ihno4,ihono,ich3o2no2,ipan,iorgntr) xmcomp=xmn case(in2o5) xmcomp=xmn*2.0 end select end if if (sedindex>0)then tempbud=tempbud+(buddry_dat(region)%dry(i,j,indices(m)) +& sum(deposition(region)%sed(i,j,:,sedindex))) / dxyp(j) * 1.E-03*xmcomp !kg m-2s-1 (s-1 at the time of writing) !for others only deposition applies else tempbud=tempbud+buddry_dat(region)%dry(i,j,allvars(j_var)%indices(m)) / dxyp(j) * 1.E-03*xmcomp !kg m-2s-1 (s-1 at the time of writing) end if end if end do !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tempbud-drydepos(region)%f2dslast(i,j,dry_index)) !drydepos(region)%f2dslast(i,j,dry_index)=tempbud allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+(tempbud-allvars(j_var)%budgetdata(i,j,1)) allvars(j_var)%budgetdata(i,j,1)=tempbud end do end do !end if else if (trim(allvars(j_var)%class)=='wet') then !gridpoints do j = j1,j2 do i=i1,i2 select case ( trim(allvars(j_var)%vname)) case('wetoa') wet_index=1 xmcomp=xmpom case('wetbc') wet_index=2 xmcomp=xmbc case ('wetso2') wet_index=3 xmcomp=xmso2 case('wetso4') wet_index=4 xmcomp=xmso4 case('wetdust') wet_index=5 xmcomp=xmdust case('wetnoy') !wet_index=6 !DMS wet_index=13 xmcomp=1 case('wetss') wet_index=7 xmcomp=xmnacl case('wetsoa') wet_index=8 xmcomp=xmpom!soa case('wetnh3') wet_index=9 xmcomp=xmnh3 case('wetnh4') wet_index=10 xmcomp=xmnh4 case DEFAULT !('default') write(gol,'("wetdep xm-var missing: ",a)') trim(allvars(j_var)%vname); call goPr end select tempbud=0.0 if ( trim(allvars(j_var)%vname)=='wetnoy')then do m=1,n_indices select case(allvars(j_var)%indices(m)) case(ino,ino2,ino3,ihno3,ino3_a,ihno4,ihono,ich3o2no2,ipan,iorgntr) xmcomp=xmn case(in2o5) xmcomp=xmn*2.0 end select if(allvars(j_var)%indices(m)>0 .and. allvars(j_var)%indices(m)<=ntracet)then tempbud=tempbud+(sum(wetdep(region)%lsp(i,j,:,indices(m))) + sum(wetdep(region)%cp(i,j,:,indices(m)))) / dxyp(j) * 1.E-03*xmcomp end if end do !write(gol,'("wet: ",a,i3,2E30.7e5)') trim(allvars(j_var)%vname),indices(m),sum(wetdep(region)%lsp(i,j,:,indices(m)))/ dxyp(j) * 1.E-03*xmcomp, sum(wetdep(region)%cp(i,j,:,indices(m))) / dxyp(j) * 1.E-03*xmcomp; call goPr allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+(tempbud-allvars(j_var)%budgetdata(i,j,1)) !kg m-2s-1 (s-1 at the time of writing) allvars(j_var)%budgetdata(i,j,1)=tempbud else do m=1,n_indices !exit if indices ==0 !since indices after 0 will all be also 0, and wetdep not defined for non-transported vars if (indices(m)<=0.or.indices(m)>ntracet) then cycle else if ( trim(allvars(j_var)%vname)=='wetnoy')then select case(indices(m)) case(ino,ino2,ino3,ihno3,ino3_a,ihno4,ihono,ich3o2no2,ipan,iorgntr) xmcomp=xmn case(in2o5) xmcomp=xmn*2.0 end select ! go through gridpoints end if tempbud=tempbud+(sum(wetdep(region)%lsp(i,j,:,indices(m))) + sum(wetdep(region)%cp(i,j,:,indices(m)))) / dxyp(j) * 1.E-03*xmcomp end if end do allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+(tempbud-allvars(j_var)%budgetdata(i,j,1)) !kg m-2s-1 (s-1 at the time of writing) allvars(j_var)%budgetdata(i,j,1)=tempbud end if end do end do !end if else if (trim(allvars(j_var)%class)=='che') then do m=1,n_indices !exit if indices ==0 !since indices after 0 will all be also 0 if (indices(m)<=0) then cycle else ! go through gridpoints do j = j1,j2 do i=i1,i2 do k=1,lmr ! add emitted mass from different size ranges (i_emi) allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime! / dxyp(j)) * 1.E-03 end do end do end do end if end do !end if else if (trim(allvars(j_var)%class)=='chep') then ! go through gridpoints do j = j1,j2 do i=i1,i2 do k=1,lmr ! density for conversion of aerosol mass densities ( ---> 1/m3 ) ! (conversion factor 1.E-03 is for g --> kg) !volume(i1:ie,j,1:lmr) = (gph_dat(region)%data(i1:ie,j,2:lmr+1)-gph_dat(region)%data(i1:ie,j,1:lmr)) * dxyp(j) vol = (gph_dat(region)%data(i,j,k+1)-gph_dat(region)%data(i,j,k)) * dxyp(j) ! dtime not needed as time is already included in calculation of prod-field ! gas-phase so4 production if (trim(allvars(j_var)%vname)=='chegpso4')then allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+AC_diag_prod(region)%prod(i,j,k,iprod_gasso4)/ dxyp(j) !liquid so4 production else if (trim(allvars(j_var)%vname)=='cheaqpso4')then allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+AC_diag_prod(region)%prod(i,j,k,iprod_aqso4)/ dxyp(j) else if (trim(allvars(j_var)%vname)=='lossch4')then allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+AC_loss(region)%prod(i,j,k,iloss_ch4)/vol else if (trim(allvars(j_var)%vname)=='lossco')then allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+AC_loss(region)%prod(i,j,k,iloss_co)/vol else if (trim(allvars(j_var)%vname)=='o3loss')then allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+AC_O3_lp(region)%prod(i,j,k,1) / vol else if (trim(allvars(j_var)%vname)=='o3prod')then allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+AC_O3_lp(region)%prod(i,j,k,2) /vol end if end do !if (trim(allvars(j_var)%vname)=='chepsoa')then ! !All soa (svoc+elvoc) ! allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+sum(AC_diag_prod(region)%prod(i,j,:,iprod_soa3dmon))/ dxyp(j) !end if end do end do if (trim(allvars(j_var)%vname)=='chepsoa')then do j =j1,j2 do i =i1,i2 !All soa (svoc+elvoc) allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+sum(AC_diag_prod(region)%prod(i,j,:,iprod_soa3dmon))/ dxyp(j) end do end do end if !end if else if (trim(allvars(j_var)%class)=='emi') then !Data can be found in emis3d in emission_pom.f90, only when declaring emissions ! which requires to write from there. ! or in pom_emis_2/3d(region,sectors), available all the time. !Could also use mode based index holder: nmode_end_XXX(nmod), where XXX=[so4,bc,pom,ss,dust] !Sedimented indices in deposition%sed do j = j1,j2 do i=i1,i2 ! add emitted mass from different size ranges (i_emi) select case ( trim(allvars(j_var)%vname)) case('emioa') emi_index=1 xmcomp=xmpom case('emibc') emi_index=2 xmcomp=xmbc case ('emiso2') emi_index=3 xmcomp=xmso2 case('emiso4') emi_index=4 xmcomp=xmso4 case('emidust') emi_index=5 xmcomp=xmdust case('emidms') emi_index=6 xmcomp=xmdms case('emiss') emi_index=7 xmcomp=xmnacl case('emiisop') emi_index=8 xmcomp=xmisop case('emivoc') emi_index=9 xmcomp=1!voc case('eminh3') emi_index=10 xmcomp=xmnh3 case('eminh4') emi_index=11 xmcomp=xmnh4 case('emico') emi_index=12 xmcomp=xmco case('emibvoc') emi_index=13 xmcomp=xmterp case('eminox') !emi_index=13 xmcomp=xmn case('emipoa') emi_index=1 xmcomp=xmpom case DEFAULT write(gol,'("emi xm-var missing: ",a)') trim(allvars(j_var)%vname); call goPr end select tempbud=0.0 if (trim(allvars(j_var)%vname)=='emivoc')then do m=1,ncb5 select case(nmhc(m)) case(ipar,ich2o,ich3oh,ihcooh) xmcomp=xmc case(ieth,iole,iald2,imcooh,ic2h6,iethoh) xmcomp=xmc*2.0 case(imgly,ic3h8,ic3h6,iacet) xmcomp=xmc*3.0 end select tempbud=tempbud+ sum(budemi_dat(region)%emi(i,j,:,nmhc(m))) / dxyp(j) * 1.E-03*xmcomp!kg(C)m-2s-1 end do do m=1,2 if (m==1)then gasind=iterp xmcomp=xmc*5*2 else gasind=iisop xmcomp=xmc*5.0 end if tempbud=tempbud+ sum(budemi_dat(region)%emi(i,j,:,gasind)) / dxyp(j) * 1.E-03*xmcomp!kg(C)m-2s-1 end do !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tempbud-emission(region)%f2dslast(i,j,emi_index)) !emission(region)%f2dslast(i,j,emi_index)=tempbud allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+(tempbud-allvars(j_var)%budgetdata(i,j,1)) allvars(j_var)%budgetdata(i,j,1)=tempbud else if (trim(allvars(j_var)%vname)=='emibvoc')then tempbud=0.0 do m=1,2!ncb5 if (m==1)then gasind=iterp xmcomp=2*5*xmc!xmterp=2*cxmisop and as C: 2*5*xmc else gasind=iisop xmcomp=5*xmc!xmisop as mass oc C end if tempbud=tempbud+ sum(budemi_dat(region)%emi(i,j,:,gasind)) / dxyp(j) * 1.E-03*xmcomp end do !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tempbud-emission(region)%f2dslast(i,j,emi_index)) !emission(region)%f2dslast(i,j,emi_index)=tempbud allvars(j_var)%data2D(i,j)= 0.0!allvars(j_var)%data2D(i,j)+(tempbud-allvars(j_var)%budgetdata(i,j,1)) allvars(j_var)%budgetdata(i,j,1)=tempbud else if (trim(allvars(j_var)%vname)=='eminox') then tempbud=0.0 indices=allvars(j_var)%indices do m=1,n_indices if (indices(m)<=0)then cycle else xmcomp=xmn !write(gol,'("eminox: ",voca,i3,E20.7e1)')trim(allvars(j_var)%vname), indices(m),sum(budemi_dat(region)%emi(i,j,:,indices(m))) / dxyp(j) * 1.E-03*xmcomp; call goPr tempbud=tempbud+ sum(budemi_dat(region)%emi(i,j,:,indices(m))) / dxyp(j) * 1.E-03*xmcomp !kg(N) m-2s-1 (s-1 at the time of writing) end if end do allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+(tempbud-allvars(j_var)%budgetdata(i,j,1)) allvars(j_var)%budgetdata(i,j,1)=tempbud else tempbud=0.0 ! add soa for emioa (emission + production) do m=1,n_indices !exit if indices ==0 !since indices after 0 will all be also 0 if (indices(m)<=0) then cycle else ! go through gridpoints tempbud = tempbud + sum(budemi_dat(region)%emi(i,j,:,indices(m))) / dxyp(j) * 1.E-03*xmcomp allvars(j_var)%data2D(i,j) = allvars(j_var)%data2D(i,j)+(tempbud-allvars(j_var)%budgetdata(i,j,1)) allvars(j_var)%budgetdata(i,j,1)=tempbud end if end do end if ! add soa production to oa emission after poa is calculated. if (allvars(j_var)%vname == 'emioa' .and. trim(allvars(j_var)%freq)=='mon') then allvars(j_var)%data2D(i,j) = allvars(j_var)%data2D(i,j)+sum(AC_diag_prod(region)%prod(i,j,:,3))/ dxyp(j) end if end do end do !end if ! CRESCENDO extensions ! else if (trim(allvars(j_var)%class)=='crescendo')then if (trim(allvars(j_var)%vname)=='od5503ddust')then cycle else if (trim(allvars(j_var)%vname)=='od550dust')then cycle else if (trim(allvars(j_var)%vname)=='od550aer')then cycle else if (trim(allvars(j_var)%vname)=='ec550aer')then cycle else if (trim(allvars(j_var)%vname)=='od440aer')then cycle else if (trim(allvars(j_var)%vname)=='od440dust')then cycle else if (trim(allvars(j_var)%vname)=='od870aer')then cycle else if (trim(allvars(j_var)%vname)=='dmsos')then cycle END if index=allvars(j_var)%indices(1) indices=allvars(j_var)%indices(:) if (allvars(j_var)%dims==3)then do k=1,lmr do j=j1,j2 do i=i1,i2 if (trim(allvars(j_var)%vname)=='ndtot')then dens = pres3d(i,j,k) / temper_dat(region)%data(i,j,k) * xmair * 1.E-3 / (m_dat(region)%data(i,j,k) * rgas) do m=1,7 allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*(tracers(i,j,k,indices(m)))*dens end do else if (trim(allvars(j_var)%vname)=='emilnox')then ! eminoc_lightning is /month so dtime/sec_month gives us per this step allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime/sec_month* eminox_lightning(region)%d3(i,j,k) / dxyp(j) else if (trim(allvars(j_var)%vname)=='mmraerh2o_1')then allvars(j_var)%data3D(i,j,k)=allvars(j_var)%data3D(i,j,k)+dtime*(h2o_mode(region,1)%d3(i,j,k))/m_dat(region)%data(i,j,k) else if (trim(allvars(j_var)%vname)=='mmraerh2o_2')then allvars(j_var)%data3D(i,j,k)=allvars(j_var)%data3D(i,j,k)+dtime*(h2o_mode(region,2)%d3(i,j,k))/m_dat(region)%data(i,j,k) else if (trim(allvars(j_var)%vname)=='mmraerh2o_3')then allvars(j_var)%data3D(i,j,k)=allvars(j_var)%data3D(i,j,k)+dtime*(h2o_mode(region,3)%d3(i,j,k))/m_dat(region)%data(i,j,k) else if (trim(allvars(j_var)%vname)=='mmraerh2o_4')then allvars(j_var)%data3D(i,j,k)=allvars(j_var)%data3D(i,j,k)+dtime*(h2o_mode(region,4)%d3(i,j,k))/m_dat(region)%data(i,j,k) else if (trim(allvars(j_var)%vname)=='chepsoa3d')then !All soa (svoc+elvoc) allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+AC_diag_prod(region)%prod(i,j,k,iprod_soa3dmon)/ dxyp(j) !number concentrations (ndtot, sfndtot and sfnd100 handled separately) else if (trim(allvars(j_var)%unit)==units(im3unit).and. .not. (trim(allvars(j_var)%vname)=='sfndtot' .or. trim(allvars(j_var)%vname)=='sfnd100' .or. trim(allvars(j_var)%vname)=='ndtot' ) ) then dens = pres3d(i,j,k) / temper_dat(region)%data(i,j,k) * xmair * 1.E-3 / (m_dat(region)%data(i,j,k) * rgas) allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*tracers(i,j,k,index)*dens else if (trim(allvars(j_var)%unit)==units(im3unit)) then !nd? output processed here dens = pres3d(i,j,k) / temper_dat(region)%data(i,j,k) * xmair * 1.E-3 / (m_dat(region)%data(i,j,k) * rgas) !1e-6 for m3->cm-3 allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*tracers(i,j,k,index)*dens else if (index>0 .and. index<=ntracet) then if (trim(allvars(j_var)%unit)==units(ivmrunit)) then xmcomp=allvars(j_var)%xmgas allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*tracers(i,j,k,index)*xmair/xmcomp/m_dat(region)%data(i,j,k) !kmole kmole-1 else if (trim(allvars(j_var)%unit)==units(immrunit)) then allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*tracers(i,j,k,index)/m_dat(region)%data(i,j,k)!kg kg-1 end if else if (index>ntracet)then ! than non transported (tracers_c=> chem_dat(region)%rm) if (trim(allvars(j_var)%unit)==units(ivmrunit)) then xmcomp=allvars(j_var)%xmgas allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*tracers_c(i,j,k,index)*xmair/xmcomp/m_dat(region)%data(i,j,k) else if (trim(allvars(j_var)%unit)==units(immrunit)) then allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*tracers_c(i,j,k,index)/m_dat(region)%data(i,j,k)!kg kg-1 else write(gol,'("user_output_aerchemmip: no case for output ",a)') trim(allvars(j_var)%vname) ; call goErr TRACEBACK status=1; return end if else write(gol,*) 'user_output_aerchemmip: missing accumulation for CRESCENDO 2D variable ',allvars(j_var)%vname ; call goErr TRACEBACK status=1; return end if end do end do end do else !2D vars do j=j1,j2 do i=i1,i2 tempbud=0.0 if (trim(allvars(j_var)%vname)=='sfmmrss')then tempbud=0.0 do m=1,size(allvars(j_var)%indices) index=allvars(j_var)%indices(m) if (index>0)then tempbud=tempbud+tracers(i,j,1,index)/m_dat(region)%data(i,j,1)!kg kg-1 end if end do !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tracers(i,j,1,index))/m_dat(region)%data(i,j,1)!*den allvars(j_var)%data2D(i,j)= tempbud else if (trim(allvars(j_var)%vname)=='co2mass')then !NOTAVAILABLE do k=1,lmr !mass of CO2, sum over levels !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tracers_c(i,j,k,iaco2))!/m_dat(region)%data(i,j,k) end do else if (trim(allvars(j_var)%vname)=='chepsoa2d')then !All soa (svoc+elvoc) !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+sum(AC_diag_prod(region)%prod(i,j,:,iprod_soa2dhour))/ dxyp(j) allvars(j_var)%data2D(i,j)= sum(AC_diag_prod(region)%prod(i,j,:,iprod_soa2dhour))/ dxyp(j)/dtime ! !!! FOLLOWING CRESCENDO VARS: Instantaneous, no accumulation else if (trim(allvars(j_var)%vname)=='sfmmrso4')then !!! Instantaneous, no accumulation tempbud=0.0 do m=1,size(allvars(j_var)%indices) index=allvars(j_var)%indices(m) if (index>0)then !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tracers(i,j,1,index))/m_dat(region)%data(i,j,1)!*den tempbud=tempbud+(tracers(i,j,1,index))/m_dat(region)%data(i,j,1)!kg kg-1 end if end do !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tracers(i,j,1,index))/m_dat(region)%data(i,j,1)!*den allvars(j_var)%data2D(i,j)= tempbud else if (trim(allvars(j_var)%vname)=='sfmmrbc')then !!! Instantaneous, no accumulation tempbud=0.0 do m=1,size(allvars(j_var)%indices) index=allvars(j_var)%indices(m) if (index>0)then tempbud=tempbud+(tracers(i,j,1,index))/m_dat(region)%data(i,j,1)!kg kg-1 end if end do allvars(j_var)%data2D(i,j)=tempbud else if (trim(allvars(j_var)%vname)=='sfmmrsoa')then !!! Instantaneous, no accumulation tempbud=0.0 do m=1,size(allvars(j_var)%indices) index=allvars(j_var)%indices(m) if (index>0)then !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tracers(i,j,1,index))/m_dat(region)%data(i,j,1)!*den tempbud=tempbud+(tracers(i,j,1,index))/m_dat(region)%data(i,j,1)!kg kg-1 end if end do allvars(j_var)%data2D(i,j)= tempbud else if (trim(allvars(j_var)%vname)=='sfmmroa')then !!! Instantaneous, no accumulation tempbud=0.0 do m=1,size(allvars(j_var)%indices) index=allvars(j_var)%indices(m) if (index>0)then !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tracers(i,j,1,index))/m_dat(region)%data(i,j,1)!*den tempbud=tempbud+(tracers(i,j,1,index))/m_dat(region)%data(i,j,1)!kg kg-1 end if end do allvars(j_var)%data2D(i,j)= tempbud else if (trim(allvars(j_var)%vname)=='sfisop')then !!! Instantaneous, no accumulation !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*tracers(i,j,1,iisop)*xmair/xmisop/m_dat(region)%data(i,j,1) allvars(j_var)%data2D(i,j)=tracers(i,j,1,iisop)*xmair/xmisop/m_dat(region)%data(i,j,1) !kilomole kilomole-1 else if (trim(allvars(j_var)%vname)=='sfdms')then !!! Instantaneous, no accumulation !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*tracers(i,j,1,idms)*xmair/xmdms/m_dat(region)%data(i,j,1) allvars(j_var)%data2D(i,j)= tracers(i,j,1,idms)*xmair/xmdms/m_dat(region)%data(i,j,1) !kmole kmole-1 else if (trim(allvars(j_var)%vname)=='sfmono')then !!! Instantaneous, no accumulation !tracers_c(i,j,k,index)*xmair/xmgas/m_dat(region)%data(i,j,k) !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*tracers(i,j,1,iterp)*xmair/xmterp/m_dat(region)%data(i,j,1) allvars(j_var)%data2D(i,j)= tracers(i,j,1,iterp)*xmair/xmterp/m_dat(region)%data(i,j,1) !kmole kmole-1 !!! no accumulation else if (trim(allvars(j_var)%vname)=='sfmmrdustabv1')then !!! Instantaneous, no accumulation allvars(j_var)%data2D(i,j)= (& (1.0-mode_fraction(rw_mode(1,3)%d3(i,j,1),1000e-9,3))*tracers(i,j,1,iduacs)+& (1.0-mode_fraction(rw_mode(1,4)%d3(i,j,1),1000e-9,4))*tracers(i,j,1,iducos)+& (1.0-mode_fraction(rw_mode(1,6)%d3(i,j,1),1000e-9,6))*tracers(i,j,1,iduaci)+& (1.0-mode_fraction(rw_mode(1,7)%d3(i,j,1),1000e-9,7))*tracers(i,j,1,iducoi))& /m_dat(region)%data(i,j,1)!kg kg-1 else if (trim(allvars(j_var)%vname)=='sfmmrdustabv10')then !!! Instantaneous, no accumulation allvars(j_var)%data2D(i,j)=(& (1.0-mode_fraction(rw_mode(1,3)%d3(i,j,1),10000e-9,3))*tracers(i,j,1,iduacs)+& (1.0-mode_fraction(rw_mode(1,4)%d3(i,j,1),10000e-9,4))*tracers(i,j,1,iducos)+& (1.0-mode_fraction(rw_mode(1,6)%d3(i,j,1),10000e-9,6))*tracers(i,j,1,iduaci)+& (1.0-mode_fraction(rw_mode(1,7)%d3(i,j,1),10000e-9,7))*tracers(i,j,1,iducoi))& /m_dat(region)%data(i,j,1)!kg kg-1 else if (trim(allvars(j_var)%vname)=='sfmmrdustbel1')then !!! Instantaneous, no accumulation allvars(j_var)%data2D(i,j)= (& (mode_fraction(rw_mode(1,3)%d3(i,j,1),1000e-9,3))*tracers(i,j,1,iduacs)+& (mode_fraction(rw_mode(1,4)%d3(i,j,1),1000e-9,4))*tracers(i,j,1,iducos)+& (mode_fraction(rw_mode(1,6)%d3(i,j,1),1000e-9,6))*tracers(i,j,1,iduaci)+& (mode_fraction(rw_mode(1,7)%d3(i,j,1),1000e-9,7))*tracers(i,j,1,iducoi))& /m_dat(region)%data(i,j,1)!kg kg-1 else if (trim(allvars(j_var)%vname)=='uas')then ! distance of single gridbox yres = dy/yref(1) xres = dx/xref(1) lat = ybeg(1) + 0.5 * yres + yres * (j-1) dxx = ae * xres * gtor * cos(lat*gtor) uwind=dxx*(pu_dat(1)%data(i,j,1) + pu_dat(1)%data(i-1,j,1))*0.5 / m_dat(1)%data(i,j,1) allvars(j_var)%data2D(i,j)= uwind else if (trim(allvars(j_var)%vname)=='vas')then yres = dy/yref(1) dyy = ae * yres * gtor vwind= dyy *(pv_dat(1)%data(i,j,1) + pv_dat(1)%data(i,j+1,1))*0.5 / m_dat(1)%data(i,j,1) allvars(j_var)%data2D(i,j)= +vwind !!!6hr output ending here else if (trim(allvars(j_var)%vname)=='ps' .and. trim(allvars(j_var)%freq) =='1hr')then allvars(j_var)%data2D(i,j)=sp_dat(1)%data(i,j,1)!Pa else if (trim(allvars(j_var)%vname)=='sfno'.and. trim(allvars(j_var)%freq) =='1hr')then allvars(j_var)%data2D(i,j)=tracers_c(i,j,1,ino)*xmair/xmno/m_dat(region)%data(i,j,1)!mole mole-1 else if (trim(allvars(j_var)%vname)=='sfnh3'.and. trim(allvars(j_var)%freq) =='1hr')then allvars(j_var)%data2D(i,j)= tracers(i,j,1,inh3)*xmair/xmnh3/m_dat(region)%data(i,j,1)!mole mole-1 else if (trim(allvars(j_var)%vname)=='sfhno3'.and. trim(allvars(j_var)%freq) =='1hr')then allvars(j_var)%data2D(i,j)= tracers(i,j,1,ihno3)*xmair/xmhno3/m_dat(region)%data(i,j,1)!mole mole-1 else if (trim(allvars(j_var)%vname)=='tas'.and. trim(allvars(j_var)%freq) =='1hr')then allvars(j_var)%data2D(i,j)= temper_dat(1)%data(i,j,1)!K else if (trim(allvars(j_var)%vname)=='sfo3'.and. trim(allvars(j_var)%freq) =='1hr')then allvars(j_var)%data2D(i,j)= tracers(i,j,1,io3)*xmair/xmo3/m_dat(region)%data(i,j,1)!mole mole-1 !!$ else if (trim(allvars(j_var)%vname)=='sfpm25'.and. trim(allvars(j_var)%freq) =='1hr')then !!$ pm_sizelimit=2.5!micrometers !!$ call PMx_Integrate_0d(region,i,j,1,pm_sizelimit,temp,status)!kgm-3 !!$ allvars(j_var)%data2D(i,j)= temp/dens/m_dat(region)%data(i,j,1)!kg kg-1 else if (trim(allvars(j_var)%vname)=='sfno2'.and. trim(allvars(j_var)%freq) =='1hr')then allvars(j_var)%data2D(i,j)=tracers_c(i,j,1,ino2)*xmair/xmno2/m_dat(region)%data(i,j,1) else if (trim(allvars(j_var)%vname)=='sfndtot'.and. trim(allvars(j_var)%freq) =='6hr')then tempbud=0.0 dens = pres3d(i,j,1) / temper_dat(region)%data(i,j,1) * xmair * 1.E-3 / (m_dat(region)%data(i,j,1) * rgas) do m=1,7 !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tracers(i,j,1,indices(m)))*dens tempbud=tempbud+(tracers(i,j,1,indices(m)))*dens end do !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tracers(i,j,1,indices(m)))*dens allvars(j_var)%data2D(i,j)=tempbud else if (trim(allvars(j_var)%vname)=='sfnd100'.and. trim(allvars(j_var)%freq) =='6hr')then tempbud=0.0 dens = pres3d(i,j,1) / temper_dat(region)%data(i,j,1) * xmair * 1.E-3 / (m_dat(region)%data(i,j,1) * rgas) do m=1,7 tempbud=tempbud+(1.0-mode_fraction(rw_mode(1,m)%d3(i,j,1),100e-9,m))*(tracers(i,j,1,indices(m)))*dens end do !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(temp) allvars(j_var)%data2D(i,j)= tempbud !! daily ouptu else if (trim(allvars(j_var)%vname)=='sfmmrdust')then tempbud=0.0 do m=1,size(allvars(j_var)%indices) index=allvars(j_var)%indices(m) if (index>0)then allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tracers(i,j,1,index))/m_dat(region)%data(i,j,1)!kg kg-1 !tempbud=tempbud+(tracers(i,j,1,index))/m_dat(region)%data(i,j,1)!*den end if end do else if (trim(allvars(j_var)%vname)=='depdustabv1')then tempbud=0.0 xmcomp=xmdust !modes_dust=(/3,4,6,7/) do m=1,size(allvars(j_var)%indices) !exit if indices ==0 !since indices after 0 will all be also 0 if (allvars(j_var)%indices(m)>0) then if (allvars(j_var)%indices(m)>69)then sedindex=-1 else sedindex=sindex(allvars(j_var)%indices(m)) end if if (sedindex>0)then do k=1,lmr tempbud=tempbud+( & (deposition(region)%sed(i,j,k,sedindex)*(1.0-mode_fraction(rw_mode(1,modes_dust(m))%d3(i,j,k),1000e-9,modes_dust(m))))) / dxyp(j) * 1.E-03*xmcomp end do end if end if !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tempbud) allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+(tempbud-allvars(j_var)%budgetdata(i,j,1)) allvars(j_var)%budgetdata(i,j,1)=tempbud end do else if (trim(allvars(j_var)%vname)=='depdustabv10')then tempbud=0.0 xmcomp=xmdust !modes_dust=(/3,4,6,7/) do m=1,size(allvars(j_var)%indices) !exit if indices ==0 !since indices after 0 will all be also 0 if (allvars(j_var)%indices(m)>0) then if (allvars(j_var)%indices(m)>69)then sedindex=-1 else sedindex=sindex(allvars(j_var)%indices(m)) end if if (sedindex>0)then do k=1,lmr tempbud=tempbud+( & (deposition(region)%sed(i,j,k,sedindex)*(1.0-mode_fraction(rw_mode(1,modes_dust(m))%d3(i,j,k),10000e-9,modes_dust(m))))) / dxyp(j) * 1.E-03*xmcomp end do end if end if !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tempbud) allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+(tempbud-allvars(j_var)%budgetdata(i,j,1)) allvars(j_var)%budgetdata(i,j,1)=tempbud end do else if (trim(allvars(j_var)%vname)=='depdustbel1')then tempbud=0.0 xmcomp=xmdust modes_dust=(/3,4,6,7/) do m=1,size(allvars(j_var)%indices) !exit if indices ==0 !since indices after 0 will all be also 0 if (allvars(j_var)%indices(m)>0) then if (allvars(j_var)%indices(m)>69)then sedindex=-1 else sedindex=sindex(allvars(j_var)%indices(m)) end if if (sedindex>0)then do k=1,lmr tempbud=tempbud+( & (deposition(region)%sed(i,j,k,sedindex)*(mode_fraction(rw_mode(1,modes_dust(m))%d3(i,j,k),1000e-9,modes_dust(m))))) / dxyp(j) * 1.E-03*xmcomp end do end if end if !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tempbud) allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+(tempbud-allvars(j_var)%budgetdata(i,j,1)) allvars(j_var)%budgetdata(i,j,1)=tempbud end do else if (trim(allvars(j_var)%vname)=='sfcWind')then yres = dy/yref(1) xres = dx/xref(1) lat = ybeg(1) + 0.5 * yres + yres * (j-1) dxx = ae * xres * gtor * cos(lat*gtor) dyy = ae * yres * gtor vwind= dyy *(pv_dat(1)%data(i,j,1) + pv_dat(1)%data(i,j+1,1))*0.5 / m_dat(1)%data(i,j,1) uwind=dxx*(pu_dat(1)%data(i,j,1) + pu_dat(1)%data(i-1,j,1))*0.5 / m_dat(1)%data(i,j,1) meanwind=(uwind**2+vwind**2)**(1./2.) allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*meanwind else if (trim(allvars(j_var)%vname)=='sfndtot')then tempbud=0.0 dens = pres3d(i,j,1) / temper_dat(region)%data(i,j,1) * xmair * 1.E-3 / (m_dat(region)%data(i,j,1) * rgas) do m=1,7 !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tracers(i,j,1,indices(m)))*dens tempbud=tempbud+(tracers(i,j,1,indices(m)))*dens end do !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tracers(i,j,1,indices(m)))*dens allvars(j_var)%data2D(i,j)=tempbud else if (trim(allvars(j_var)%vname)=='sfnd100')then tempbud=0.0 dens = pres3d(i,j,1) / temper_dat(region)%data(i,j,1) * xmair * 1.E-3 / (m_dat(region)%data(i,j,1) * rgas) do m=1,7 tempbud=tempbud+(1.0-mode_fraction(rw_mode(1,m)%d3(i,j,1),100e-9,m))*(tracers(i,j,1,indices(m)))*dens end do !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(temp) allvars(j_var)%data2D(i,j)= tempbud else if (trim(allvars(j_var)%vname)=='sfsh')then allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*humid_dat(1)%data(i,j,1) else if (trim(allvars(j_var)%vname)=='sfmmraerh2o_1')then allvars(j_var)%data2D(i,j)=allvars(j_var)%data2D(i,j)+dtime*(h2o_mode(region,1)%d3(i,j,1))/m_dat(region)%data(i,j,1)!kg kg-1 else if (trim(allvars(j_var)%vname)=='sfmmraerh2o_2')then allvars(j_var)%data2D(i,j)=allvars(j_var)%data2D(i,j)+dtime*(h2o_mode(region,2)%d3(i,j,1))/m_dat(region)%data(i,j,1)!kg kg-1 else if (trim(allvars(j_var)%vname)=='sfmmraerh2o_3')then allvars(j_var)%data2D(i,j)=allvars(j_var)%data2D(i,j)+dtime*(h2o_mode(region,3)%d3(i,j,1))/m_dat(region)%data(i,j,1)!kg kg-1 else if (trim(allvars(j_var)%vname)=='sfmmraerh2o_4')then allvars(j_var)%data2D(i,j)=allvars(j_var)%data2D(i,j)+dtime*(h2o_mode(region,4)%d3(i,j,1))/m_dat(region)%data(i,j,1)!kg kg-1 else if (trim(allvars(j_var)%vname)=='tas')then allvars(j_var)%data2D(i,j)=allvars(j_var)%data2D(i,j)+dtime*temper_dat(1)%data(i,j,1)!K else if (trim(allvars(j_var)%vname)=='sfo3')then allvars(j_var)%data2D(i,j)=allvars(j_var)%data2D(i,j)+dtime*tracers(i,j,1,io3)*xmair/xmo3/m_dat(region)%data(i,j,1)!mole mole-1 !!$ else if (trim(allvars(j_var)%vname)=='sfpm25')then !!$ pm_sizelimit=2.5!micrometers !!$ call PMx_Integrate_0d(region,i,j,1,pm_sizelimit,temp,status)!kg !!$ allvars(j_var)%data2D(i,j)=allvars(j_var)%data2D(i,j) + dtime*temp/m_dat(region)%data(i,j,1)/dens!kg kg-1 else if (trim(allvars(j_var)%vname)=='sfno2')then allvars(j_var)%data2D(i,j)=allvars(j_var)%data2D(i,j) + dtime*tracers_c(i,j,1,ino2)*xmair/xmno2/m_dat(region)%data(i,j,1) else if (trim(allvars(j_var)%vname)=='wetdms' .or.& trim(allvars(j_var)%vname)=='wetno3' .or.& trim(allvars(j_var)%vname)=='wetnh4' .or.& trim(allvars(j_var)%vname)=='wethno3' .or.& trim(allvars(j_var)%vname)=='wetdust' ) then ! Wet index not used for CRESCENDO vars ! Changed to own budget variable to fix budget calculation on the ! variable holder in allvars(jvar). This way you can use the same var ! in more than one time means (hr, daily, monthly) select case ( trim(allvars(j_var)%vname)) case('wetdms') !wet_index=6 xmcomp=xmdms case('wetno3') !wet_index=11 xmcomp=xmno3 case('wethno3') !wet_index=12 xmcomp=xmhno3 case('wetnh4') !wet_index=10 xmcomp=xmnh4 case('wetdust') !wet_index=5 xmcomp=xmdust case DEFAULT write(gol,'("CRESCENDO wetdep xm-var missing: ",a)') trim(allvars(j_var)%vname); call goPr end select tempbud=0.0 do m=1,n_indices !exit if indices ==0 !since indices after 0 will all be also 0 if (allvars(j_var)%indices(m)<70)then if (allvars(j_var)%indices(m)<=0) then cycle else ! go through gridpoints tempbud=tempbud+(sum(wetdep(region)%lsp(i,j,:,indices(m))) + sum(wetdep(region)%cp(i,j,:,indices(m)))) / dxyp(j) * 1.E-03*xmcomp !kg m-2s-1 (s-1 at the time of writing) end if end if end do !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tempbud-wetdepos(region)%f2dslast(i,j,wet_index)) !wetdepos(region)%f2dslast(i,j,wet_index)=tempbud allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+(tempbud-allvars(j_var)%budgetdata(i,j,1)) allvars(j_var)%budgetdata(i,j,1)=tempbud else if (trim(allvars(j_var)%vname)=='seddust')then tempbud=0.0 xmcomp=xmdust do m=1,n_indices !exit if indices ==0 !since indices after 0 will all be also 0 if (indices(m)<=0) then cycle else if (allvars(j_var)%indices(m)>69)then sedindex=-1 else sedindex=sindex(allvars(j_var)%indices(m)) end if if (sedindex>0)then tempbud=tempbud+( & sum(deposition(region)%sed(i,j,:,sedindex))) / dxyp(j) * 1.E-03*xmcomp !kg m-2s-1 (s-1 at the time of writing) end if end if !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tempbud) allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+(tempbud-allvars(j_var)%budgetdata(i,j,1)) allvars(j_var)%budgetdata(i,j,1)=tempbud end do else if (trim(allvars(j_var)%vname)=='drydms' .or.& trim(allvars(j_var)%vname)=='dryno3' .or.& trim(allvars(j_var)%vname)=='drynh4' .or.& trim(allvars(j_var)%vname)=='dryhno3' .or.& trim(allvars(j_var)%vname)=='dryno2' .or.& trim(allvars(j_var)%vname)=='dryno' .or.& trim(allvars(j_var)%vname)=='drypan' .or.& trim(allvars(j_var)%vname)=='drydust' & ) then ! Dry index not used for CRESCENDO vars ! Changed to own budget variable to fix budget calculation on the ! variable holder in allvars(jvar). This way you can use the same var ! in more than one time means (hr, daily, monthly) select case ( trim(allvars(j_var)%vname)) case('drynh4') !dry_index=10 xmcomp=xmnh4 case('dryno3') !dry_index=13 xmcomp=xmno3 case('drydms') !dry_index=6 xmcomp=xmdms case('dryhno3') !dry_index=14 xmcomp=xmhno3 case('dryno2') !dry_index=15 xmcomp=xmno2 case('dryno') !dry_index=16 xmcomp=xmno case('drypan') !dry_index=17 xmcomp=xmpan case('drydust') !dry_index=5 xmcomp=xmdust case DEFAULT write(gol,'("CRESCENDO drydep xm-var missing: ",a)') trim(allvars(j_var)%vname); call goPr end select tempbud=0.0 do m=1,n_indices !exit if indices ==0 !since indices after 0 will all be also 0 if (indices(m)<=0) then cycle else if (allvars(j_var)%indices(m)>69)then sedindex=-1 else sedindex=sindex(allvars(j_var)%indices(m)) end if !only M7 aerosol tracers (and msa/nh4/no3) are sedimented and deposited if (sedindex>0)then tempbud=tempbud+(buddry_dat(region)%dry(i,j,indices(m)) +& sum(deposition(region)%sed(i,j,:,sedindex))) / dxyp(j) * 1.E-03*xmcomp !kg m-2s-1 (s-1 at the time of writing) else tempbud=tempbud+buddry_dat(region)%dry(i,j,indices(m)) / dxyp(j) * 1.E-03*xmcomp !kg m-2s-1 (s-1 at the time of writing) end if end if end do !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tempbud-drydepos(region)%f2dslast(i,j,dry_index)) !drydepos(region)%f2dslast(i,j,dry_index)=tempbud allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+(tempbud-allvars(j_var)%budgetdata(i,j,1)) allvars(j_var)%budgetdata(i,j,1)=tempbud else if (trim(allvars(j_var)%vname)=='emioa')then emi_index=1 xmcomp=xmpom tempbud=0.0 do m=1,n_indices !exit if indices ==0 !since indices after 0 will all be also 0 if (indices(m)<=0) then cycle else tempbud=tempbud+ sum(budemi_dat(region)%emi(i,j,:,indices(m))) / dxyp(j) * 1.E-03*xmcomp !kg m-2s-1 (s-1 at the time of writing) !OLDWAY which is not in use for CRESCENDO vars: !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tempbud-emission(region)%f2dslast(i,j,emi_index)) !emission(region)%f2dslast(i,j,emi_index)=tempbud end if end do ! emioa is requested as emi+chepsoa so add soa production !if (trim(allvars(j_var)%freq)=='mon')then ! tempbud= tempbud + AC_diag_prod(region)%prod(i,j,k,3)/ dxyp(j) !kg m-2s-1 (s-1 at the time of writing) !else if (trim(allvars(j_var)%freq)=='6hr')then ! tempbud= tempbud + AC_diag_prod(region)%prod(i,j,k,4)/ dxyp(j) !kg m-2s-1 (s-1 at the time of writing) !else ! write(gol,'("Uknown output frequency for production of soa in calculation of : ",a)')trim(allvars(j_var)%freq); call goErr ! TRACEBACK ! status=1; return !end if if (trim(allvars(j_var)%freq)=='6hr'.or.trim(allvars(j_var)%freq)=='1hr')then allvars(j_var)%data2D(i,j)= (tempbud-allvars(j_var)%budgetdata(i,j,1))/dtime!timerange for instantaneous data is half of dynamical timestep /dtime allvars(j_var)%budgetdata(i,j,1)=tempbud else if (trim(allvars(j_var)%freq)=='mon'.or.trim(allvars(j_var)%freq)=='day')then allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+tempbud-allvars(j_var)%budgetdata(i,j,1) allvars(j_var)%budgetdata(i,j,1)=tempbud else write(gol,'("Uknown output frequency: ",a)')trim(allvars(j_var)%vname); call goErr TRACEBACK status=1; return end if else if (trim(allvars(j_var)%vname)=='emiisop')then tempbud=0.0 gasind=iisop xmcomp=xmisop emi_index=8 tempbud=tempbud+ sum(budemi_dat(region)%emi(i,j,:,gasind)) / dxyp(j) * 1.E-03*xmcomp !kg m-2s-1 (s-1 at the time of writing) if (trim(allvars(j_var)%freq)=='6hr'.or.trim(allvars(j_var)%freq)=='1hr')then allvars(j_var)%data2D(i,j)= (tempbud-allvars(j_var)%budgetdata(i,j,1))/dtime!timerange for instantaneous data is half of dynamical timestep /dtime! instantaneous values -> /dtime !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tempbud-allvars(j_var)%budgetdata(i,j,1)) allvars(j_var)%budgetdata(i,j,1)=tempbud !kg m-2s-1 (s-1 at the time of writing) else if (trim(allvars(j_var)%freq)=='mon'.or.trim(allvars(j_var)%freq)=='day')then allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j) + (tempbud-allvars(j_var)%budgetdata(i,j,1)) !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tempbud-allvars(j_var)%budgetdata(i,j,1)) allvars(j_var)%budgetdata(i,j,1)=tempbud !kg m-2s-1 (s-1 at the time of writing) else write(gol,'("Uknown output frequency: ",a)')trim(allvars(j_var)%vname); call goErr TRACEBACK status=1; return end if !allvars(j_var)%data2D(i,j)= tempbud-allvars(j_var)%budgetdata(i,j,1) !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tempbud-allvars(j_var)%budgetdata(i,j,1)) !allvars(j_var)%budgetdata(i,j,1)=tempbud else if (trim(allvars(j_var)%vname)=='emimono')then tempbud=0.0 gasind=iterp xmcomp=xmterp emi_index=13 tempbud=tempbud+ sum(budemi_dat(region)%emi(i,j,:,gasind)) / dxyp(j) * 1.E-03*xmcomp if (trim(allvars(j_var)%freq)=='6hr'.or.trim(allvars(j_var)%freq)=='1hr')then !Instantaneous values requested allvars(j_var)%data2D(i,j)= (tempbud-allvars(j_var)%budgetdata(i,j,1))/dtime!timerange for instantaneous data is half of dynamical timestep /dtime! instantaneous values -> /dtime : kgm-2s-1 allvars(j_var)%budgetdata(i,j,1)=tempbud !kg m-2s-1 (s-1 at the time of writing) else if (trim(allvars(j_var)%freq)=='mon'.or.trim(allvars(j_var)%freq)=='day')then !mean values allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j) + (tempbud-allvars(j_var)%budgetdata(i,j,1)) allvars(j_var)%budgetdata(i,j,1)=tempbud !kg m-2s-1 (s-1 at the time of writing) else write(gol,'("Uknown output frequency: ",a)')trim(allvars(j_var)%vname); call goErr TRACEBACK status=1; return end if else if (trim(allvars(j_var)%vname)=='emidms')then tempbud=0.0 gasind=idms emi_index=6 xmcomp=xmdms tempbud=tempbud+ sum(budemi_dat(region)%emi(i,j,:,gasind)) / dxyp(j) * 1.E-03*xmcomp if (trim(allvars(j_var)%freq)=='6hr')then allvars(j_var)%data2D(i,j)= (tempbud-allvars(j_var)%budgetdata(i,j,1))/dtime!timerange for instantaneous data is half of dynamical timestep /dtime allvars(j_var)%budgetdata(i,j,1)=tempbud !kg m-2s-1 (s-1 at the time of writing) else write(gol,'("CRESCENDO output frequency not implemented: ",a)')trim(allvars(j_var)%vname); call goErr TRACEBACK status=1; return end if else if (trim(allvars(j_var)%vname)=='emiss')then emi_index=7 xmcomp=xmnacl do m=1,n_indices !exit if indices ==0 !since indices after 0 will all be also 0 if (indices(m)<=0) then cycle else ! go through gridpoints tempbud=tempbud+ sum(budemi_dat(region)%emi(i,j,:,indices(m))) / dxyp(j) * 1.E-03*xmcomp end if end do if (trim(allvars(j_var)%freq)=='6hr')then allvars(j_var)%data2D(i,j)= (tempbud-allvars(j_var)%budgetdata(i,j,1))/dtime!timerange for instantaneous data is half of dynamical timestep /dtime !kg m-2s-1 (s-1 at the time of writing) allvars(j_var)%budgetdata(i,j,1)=tempbud else if (trim(allvars(j_var)%freq)=='mon'.or.trim(allvars(j_var)%freq)=='day')then allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j) + (tempbud-allvars(j_var)%budgetdata(i,j,1))*dtime allvars(j_var)%budgetdata(i,j,1)=tempbud end if else if (trim(allvars(j_var)%unit)=='cm-3') then dens = pres3d(i,j,1) / temper_dat(region)%data(i,j,1) * xmair * 1.E-3 / (m_dat(region)%data(i,j,1) * rgas) !1e-6 for m3->cm-3 allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*tracers(i,j,1,index)*dens*1e-6 else if (index>0 .and. index<=ntracet) then allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*tracers(i,j,1,index)/m_dat(region)%data(i,j,1)!kg kg-1 else if (index>ntracet) then ! than non transported (tracers_c=> chem_dat(region)%rm) allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*tracers_c(i,j,1,index)/m_dat(region)%data(i,j,1)!kg kg-1 else write(gol,*) 'user_output_aerchemmip: missing accumulation for CRESCENDO 3D variable ',allvars(j_var)%vname ; call goErr TRACEBACK status=1; return end if end do end do end if !end if else if (trim(allvars(j_var)%class)=='mmr')then 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) !do all mode-specific tracers for given compound dens = pres3d(i,j,k) / temper_dat(region)%data(i,j,k) * xmair * 1.E-3 / (m_dat(region)%data(i,j,k) * rgas) do imode=1,size(allvars(j_var)%indices) index=allvars(j_var)%indices(imode) ! first variables that that are not tracers (transported or non-transported) if (index<=0) then if (trim(allvars(j_var)%vname)=='mmrpm2p5')then pm_sizelimit=2.5!micrometers !call PMx_Integrate_3d(region,pm_sizelimit,temp_pm,status) CALL pmx_integrate_0d( region, i, j, k, pm_sizelimit, temp, status ) !returns kg/m3 allvars(j_var)%data3D(i,j,k)=allvars(j_var)%data3D(i,j,k)+dtime*temp/dens/m_dat(region)%data(i,j,k) !kg(pm)/m3 /(1/m3) /kg(air)->kg(pm)/kg(air) else if (trim(allvars(j_var)%vname)=='mmrpm1')then pm_sizelimit=1.0!micrometers !call PMx_Integrate_3d(region,pm_sizelimit,temp_pm,status) CALL pmx_integrate_0d( region, i, j, k, pm_sizelimit, temp, status ) !returns kg/m3 allvars(j_var)%data3D(i,j,k)=allvars(j_var)%data3D(i,j,k)+dtime*temp/dens/m_dat(region)%data(i,j,k)!kg kg-1 else if (trim(allvars(j_var)%vname)=='mmrpm10')then pm_sizelimit=10.0!micrometers !call PMx_Integrate_3d(region,pm_sizelimit,temp_pm,status) CALL pmx_integrate_0d( region, i, j, k, pm_sizelimit, temp, status ) !returns kg/m3 !kg/m3/m-3/kg allvars(j_var)%data3D(i,j,k)=allvars(j_var)%data3D(i,j,k)+dtime*temp/dens/m_dat(region)%data(i,j,k)!kg kg-1 else if (trim(allvars(j_var)%vname)=='mmraerh2o')then allvars(j_var)%data3D(i,j,k)=allvars(j_var)%data3D(i,j,k)+dtime*(h2o_mode(region,1)%d3(i,j,k)+h2o_mode(region,2)%d3(i,j,k)+h2o_mode(region,3)%d3(i,j,k)+h2o_mode(region,4)%d3(i,j,k))/m_dat(region)%data(i,j,k)!kg kg-1 end if !exit else !Transported variables are stored in different array (tracers=> mass_dat(region)%rm) if (index<=ntracet) then allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*tracers(i,j,k,index)/m_dat(region)%data(i,j,k)!kg kg-1 else ! than non transported (tracers_c=> chem_dat(region)%rm) allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*tracers_c(i,j,k,index)/m_dat(region)%data(i,j,k)!kg kg-1 end if end if end do end do end do end do !end if else if (trim(allvars(j_var)%class)=='gas')then do k=1,lmr do j=j1,j2 do i=i1,i2 index=allvars(j_var)%indices(1) if (index<=0) then cycle else if (trim(allvars(j_var)%vname)=='h2o') then allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*humid_dat(1)%data(i,j,k)/xmh2o*xmair !mole mole-1 !tracers_c(i,j,k,index)*xmair/xmgas/m_dat(region)%data(i,j,k) else if (trim(allvars(j_var)%vname)=='voc')then ! ACTUALLY not requested do m=1,ncb5 gasind=nmhc(m) xmgas=xmcb5(m) select case(nmhc(m)) case(ipar,ich2o,ich3oh,ihcooh) xmcomp=xmc case(ieth,iole,iald2,imcooh,ic2h6,iethoh) xmcomp=xmc*2.0 case(imgly,ic3h8,ic3h6,iacet) xmcomp=xmc*3.0 end select allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*tracers(i,j,k,gasind)*xmair/xmcomp/m_dat(region)%data(i,j,k) !kmole kmole-1 end do do m=1,2 if (m==1)then gasind=iterp xmcomp=xmc*5*2 else gasind=iisop xmcomp=xmc*5.0 end if allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*tracers(i,j,k,gasind)*xmair/xmcomp/m_dat(region)%data(i,j,k) !kmole kmole-1 end do else xmgas=allvars(j_var)%xmgas if (index >= ntracet) then allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*tracers_c(i,j,k,index)*xmair/xmgas/m_dat(region)%data(i,j,k) !kmole kmole-1 else allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*tracers(i,j,k,index)*xmair/xmgas/m_dat(region)%data(i,j,k) !kmole kmole-1 end if end if end do end do end do !end if else if (trim(allvars(j_var)%class)=='ps6h')then do i=i1,i2 do j=j1,j2 if (trim(allvars(j_var)%vname)=='ps')then allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*sp_dat(1)%data(i,j,1)!Pa end if end do end do !end if ! 1 hourly surface variables else if (trim(allvars(j_var)%class)=='sf1h')then do i=i1,i2 do j=j1,j2 if (trim(allvars(j_var)%vname)=='ps')then allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*sp_dat(1)%data(i,j,1)!Pa else if (trim(allvars(j_var)%vname)=='sfno2')then allvars(j_var)%data2D(i,j)=allvars(j_var)%data2D(i,j)+dtime*tracers_c(i,j,1,ino2)*xmair/xmno2/m_dat(region)%data(i,j,1)*dtime !kmole kmole-1 else if (trim(allvars(j_var)%vname)=='sfo3')then allvars(j_var)%data2D(i,j)=allvars(j_var)%data2D(i,j)+dtime*tracers(i,j,1,io3)*xmair/xmo3/m_dat(region)%data(i,j,1) !kmole kmole-1 else if (trim(allvars(j_var)%vname)=='sfpm25')then pm_sizelimit=2.5!micrometers call PMx_Integrate_0d(region,i,j,1,pm_sizelimit,temp,status) !allvars(j_var)%data2D(i,j)=allvars(j_var)%data2D(i,j)+dtime*temp/m_dat(region)%data(i,j,1)/dens!kg kg-1 allvars(j_var)%data2D(i,j)=temp/m_dat(region)%data(i,j,1)/dens*dtime!kg kg-1 else if (trim(allvars(j_var)%vname)=='tas')then allvars(j_var)%data2D(i,j)=allvars(j_var)%data2D(i,j)+dtime*temper_dat(1)%data(i,j,1)!K end if end do end do !end if ! surface daily variables else if (trim(allvars(j_var)%class)=='sf1d')then do i=i1,i2 do j=j1,j2 if (trim(allvars(j_var)%vname)=='ps')then allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*sp_dat(1)%data(i,j,1)!Pa else if (trim(allvars(j_var)%vname)=='toz')then ! conversion (in order of execution): !kg->g (1e3), g->molec (xmo3), m2->cm2(1e-4) , molec-> molec/m2(dxyp) , molec/cm2->dobson (DOBS) allvars(j_var)%data2D(i,j)=allvars(j_var)%data2D(i,j)+dtime*sum(tracers(i,j,:,io3)*1e3/xmo3*Avog*1e-4/dxyp(j)/Dobs) !GU else if (trim(allvars(j_var)%vname)=='sfo3max')then if (tracers(i,j,1,io3)> allvars(j_var)%data2D(i,j)) then allvars(j_var)%data2D(i,j)=tracers(i,j,1,io3)*xmair/xmo3/m_dat(region)%data(i,j,1) !kmole kmole-1 end if else if (trim(allvars(j_var)%vname)=='maxpblz')then if (conv_dat(1)%blh(i,j)> allvars(j_var)%data2D(i,j)) then allvars(j_var)%data2D(i,j)= conv_dat(1)%blh(i,j)!m end if else if (trim(allvars(j_var)%vname)=='minpblz')then if (conv_dat(1)%blh(i,j)< allvars(j_var)%data2D(i,j)) then allvars(j_var)%data2D(i,j)= conv_dat(1)%blh(i,j)!m end if else if (trim(allvars(j_var)%vname)=='tas')then allvars(j_var)%data2D(i,j)=allvars(j_var)%data2D(i,j)+dtime*temper_dat(1)%data(i,j,1)!K else if (trim(allvars(j_var)%vname)=='tasmin')then if (temper_dat(1)%data(i,j,1)< allvars(j_var)%data2D(i,j)) then allvars(j_var)%data2D(i,j)= temper_dat(1)%data(i,j,1)!K end if else if (trim(allvars(j_var)%vname)=='tasmax')then if (temper_dat(1)%data(i,j,1)> allvars(j_var)%data2D(i,j)) then allvars(j_var)%data2D(i,j)= temper_dat(1)%data(i,j,1)!K end if else if (trim(allvars(j_var)%vname)=='pr')then allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(lsp_dat(1)%data(i,j,1)+ cp_dat(1)%data(i,j,1))/rol!kgm-2s-1 end if end do end do !end if else if (trim(allvars(j_var)%class)=='zonal')then ! zonal mean needs to be done afterwards... ! So just save the variables as 3D do i=i1,i2 do j=j1,j2 do k=1,lmr if (trim(allvars(j_var)%vname)=='ch4')then allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*tracers(i,j,k,ich4)/m_dat(region)%data(i,j,k)*xmair/xmch4!mol mol-1 else if (trim(allvars(j_var)%vname)=='hno3')then allvars(j_var)%data3D(i,j,k)=allvars(j_var)%data3D(i,j,k)+dtime*tracers(i,j,k,ihno3)/m_dat(region)%data(i,j,k)*xmair/xmhno3!mol mol-1 else if (trim(allvars(j_var)%vname)=='o3')then allvars(j_var)%data3D(i,j,k)=allvars(j_var)%data3D(i,j,k)+dtime*tracers(i,j,k,io3)/m_dat(region)%data(i,j,k)*xmair/xmo3!mol mol-1 else if (trim(allvars(j_var)%vname)=='ta')then allvars(j_var)%data3D(i,j,k)=allvars(j_var)%data3D(i,j,k)+dtime*temper_dat(1)%data(i,j,k)!K else if (trim(allvars(j_var)%vname)=='h2o')then allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+0.0!dtime*humid_dat(1)%data(i,j,k)/xmh2o*xmair!mol mol-1 else if (trim(allvars(j_var)%vname)=='zg')then allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*gph_dat(1)%data(i,j,k)!m else if (trim(allvars(j_var)%vname)=='ho2')then allvars(j_var)%data3D(i,j,k)=allvars(j_var)%data3D(i,j,k)+dtime*tracers_c(i,j,k,iho2)/m_dat(region)%data(i,j,k)*xmair/xmho2!mol mol-1 else if (trim(allvars(j_var)%vname)=='oh')then allvars(j_var)%data3D(i,j,k)=allvars(j_var)%data3D(i,j,k)+dtime*tracers_c(i,j,k,ioh)/m_dat(region)%data(i,j,k)*xmair/xmoh!mol mol-1 else if (trim(allvars(j_var)%vname)=='noy')then do icomp=1,size( allvars(j_var)%indices(:)) index= allvars(j_var)%indices(icomp) xmcomp=ra(index) if (allvars(j_var)%indices(icomp)>0) then !write(gol,'(": ",a,i3,E20.7e2)') trim(allvars(j_var)%vname),index,tracers(i,j,k,index)/m_dat(region)%data(i,j,k)*xmair/xmcomp; call goPr if (allvars(j_var)%indices(icomp)<70) then allvars(j_var)%data3D(i,j,k)=allvars(j_var)%data3D(i,j,k)+dtime*tracers(i,j,k,index)/m_dat(region)%data(i,j,k)*xmair/xmcomp!??? else allvars(j_var)%data3D(i,j,k)=allvars(j_var)%data3D(i,j,k)+dtime*tracers_c(i,j,k,index)/m_dat(region)%data(i,j,k)*xmair/xmcomp!??? end if end if end do end if end do end do end do else ! optics and fixed are not accumulated here ! optics in optics_calc ! fixed not at all if (trim(allvars(j_var)%class)/='optics' .and. trim(allvars(j_var)%class)/='fixed') then write(gol,*) 'output_aerhemmip_accumulate: output class not found!!!',trim(allvars(j_var)%vname),trim(allvars(j_var)%class) !call goPr call goErr status=1 return end if end if end do end if ! zero accumulated budget variables for the amount between output steps AC_diag_prod(region)%prod(i1:i2,j1:j2,:,1)=0.0 AC_diag_prod(region)%prod(i1:i2,j1:j2,:,2)=0.0 AC_diag_prod(region)%prod(i1:i2,j1:j2,:,3)=0.0 AC_diag_prod(region)%prod(i1:i2,j1:j2,:,4)=0.0 AC_loss(region)%prod(i1:i2,j1:j2,:,1:2)=0.0 AC_O3_lp(region)%prod(i1:i2,j1:j2,:,1:2)=0.0 deallocate(gphx) deallocate(tx) deallocate(temp_pm) deallocate(pres3d) deallocate(pres3dh) call GO_Timer_Start( itim_accu_opt, status ) call calculate_optics(1,dhour,l_do_ec550aer_only,status) call GO_Timer_End( itim_accu_opt, status ) call GO_Timer_End( itim_accu, status ) ! IF_NOTOK_RETURN(status=1) call goLabel() !status = 1 end subroutine accumulate_data subroutine output_aerchemmip_done(status) use partools, only: isRoot,myid implicit none integer :: status character(len=*), parameter :: rname = mname//'/output_aerchemmip_done' integer :: j_var, region call goLabel(rname) region = 1 if (isroot) then do j_var=1,n_vars call MDF_Close( allvars(j_var)%fileunit, status ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) end do end if deallocate(wdep_out ) deallocate(drydepos(region)%f2dslast) deallocate(wetdepos(region)%f2dslast) deallocate(emission(region)%f2dslast) deallocate(dimension_data%lon) deallocate(dimension_data%lat) deallocate(dimension_data%lev) do j_var=1,n_vars deallocate(allvars(j_var)%data2D) deallocate(allvars(j_var)%data3D) end do deallocate(allvars) deallocate(fixedvars) call goLabel() status = 0 end subroutine output_aerchemmip_done !subroutine add_variable(itm5,shortname,longname,unit,positive,standardname,varid,fileunit,filename,data_dims,status,class,table,pxmgas) subroutine add_variable(itm5,shortname,longname,unit,data_dims,status,class,table,pxmgas) use chem_param, only: mode_end_so4,mode_end_pom,mode_end_bc,mode_end_ss,mode_end_dust use global_data, only: outdir use datetime, only: tau2date, date2tau use dims, only: itau,itaue,itaut implicit none integer:: itm5 character(len=*),intent(in)::shortname character(len=*),intent(in)::longname character(len=*)::unit character(len=30)::standardname character(len=*)::table character(len=*),optional::class real,optional::pxmgas integer:: data_dims integer,intent(out)::status !LOCAL character(len=4)::positive='' integer:: fileunit=-1 !defined only when creating a file integer:: varid=-1! defined when opening ncfile !character(len=120)::filename character(len=30)::table_id !character(len=30)::source_id !character(len=30)::experiment_id character(len=30)::member_id !character(len=30)::grid_label character(len=30)::time_range character(len=200)::filename1 character(len=10)::freq real,dimension(:,:),pointer::data2D ! real,dimension(:,:),pointer::dataZonal real,dimension(:,:,:),pointer::data3D real,dimension(:,:,:),pointer::budget character(len=*), parameter :: rname = mname//'/output_aerchemmip_add_variable' integer ::i1,i2,j1,j2,jmr,imr,lmr integer, dimension(6) :: idater, idateend, idatett integer :: endmonth, endday character(len=30) :: idates call tau2date(itau,idater) ! define frequency from table if (trim(table)=='AERhr')then !table id depends on variable table_id=table freq='1hr' else if (trim(table)=='AER6hr')then !table id depends on variable table_id=table freq='6hr' else if( trim(table)=='AERmon'.or.trim(table)=='AERmonZ'.or.trim(table)=='Emon')then !table id depends on variable table_id=table freq='mon' else if(trim(table)=='AERday')then !table id depends on variable table_id=table freq='day' else if(trim(table)=='AERfx')then !table id depends on variable table_id=table freq='na' else freq='freq-nondefined' table_id='table-nondefined' end if ! CREATE date string for output ! ! ATM assumed that the output is initilised at the beginninh of year endmonth=12 endday=31 ! if (freq=='mon')then ! By default AERCHEMMIP runs are done by 1-year chunks -> idater(2) - idater(2)+11 write(idates, '(i4,i2.2,a,i4,i2.2)') idater(1), idater(2),'-', idater(1),endmonth else if(freq=='day')then !time range form Jan-1 ->Dec-31x write(idates, '(i4,2i2.2,a,i4,2i2.2)') idater(1), idater(2), idater(3),'-', idater(1), endmonth, endday else if(freq=='1hr')then write(idates, '(i4,2i2.2,2a2,a,i4,2i2.2,2a2)') idater(1), idater(2), idater(3),'00','00','-', idater(1), endmonth, endday, '23', '59' else if (freq=='6hr')then write(idates, '(i4,2i2.2,2a2,a,i4,2i2.2,2a2)') idater(1), idater(2), idater(3),'00','00','-', idater(1), endmonth, endday,'18','00' end if !statf(region)%fname = trim(outdir)//'/'//& ! trim(f_dataid) //'_'//& ! trim(f_model) //'_'//& ! trim(aerocom_exper) //'_'//& ! trim(f_dimstat)//'_'//& ! trim(idates) //'_'//& ! trim(aerocom_freq) //'.nc' call goLabel(rname) call GO_Timer_Start( itim_addvar, status ) IF_NOTOK_RETURN(status=1) !outdir='output' ! temporary standardname=longname ! source_id constant !source_id='EC-EARTH-AerChem' ! experiment depends on run !experiment_id='exp_dynamic' member_id='r'//trim(realization_i)//'i'//trim(initialization_i)//'p'//trim(physics_i)//'f'//trim(forcing_i) !grid_label='3x2_degrees' ! time range has divider in place since it can be omitted alltogether with non-time dependendent variables if (trim(table)=='AERfx')then time_range='' else time_range='_'//trim(idates) end if ! for fixed variables time range should not be written n_vars=n_vars+1 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 if (trim(shortname)=='phalf') then allocate(budget(i1:i2,j1:j2,1:lmr+1)) else allocate(budget(i1:i2,j1:j2,1:lmr)) end if budget(:,:,:)=0.0 ! 2D variables if (data_dims==2) then !Allocate holders for data allocate(allvars(n_vars)%data2D(i1:i2,j1:j2)) allocate(allvars(n_vars)%data3D(1,1,1)) ! allocate local variables allocate(data2D(i1:i2,j1:j2)) allocate(data3D(1,1,1)) ! zero local data holders data2D(:,:)=0.0 ! dataZonal(:,:)=0.0 data3D(:,:,:)=0.0 !init variable ! set default for minimum variables to high value if (shortname=='minpblz' .or. shortname=='tasmin')then data2D(:,:)=1000000.0 end if !create filename if (trim(class)=='crescendo')then ! filename1= trim(outdir)//'/'//trim(shortname)//'_'//trim(table_id)//'_'//trim(source_id)//'_'//trim(experiment_id)//'_'//trim(member_id)//'_'//trim(grid_label)//'_'//trim(time_range)//trim('.nc') filename1= trim(outdir)//'/'//trim(shortname)//'_'//trim(class)//'_'//trim(table_id)//'_'//trim(source_id)//'_'//trim(experiment_id)//'_'//trim(member_id)//'_'//trim(grid_label)//trim(time_range)//trim('.nc') else filename1= trim(outdir)//'/'//trim(shortname)//'_'//trim(table_id)//'_'//trim(source_id)//'_'//trim(experiment_id)//'_'//trim(member_id)//'_'//trim(grid_label)//trim(time_range)//trim('.nc') end if allvars(n_vars)=varfile(itm5,shortname,longname,unit,positive,standardname,data2D,data3D,budget,-1,-1,& filename1,2,-1,-1,-1,-1,-1,-1,-1,-1,-1,data_dims,freq,class,(/0,0,0,0,0,0,0,0,0,0,0/),pxmgas,table_id) !! LEFT HERE on purpose to see what variables go where in above statement !!$ allvars(n_vars)%itm5=itm5 !!$ allvars(n_vars)%vname=shortname !!$ allvars(n_vars)%lname=longname !!$ allvars(n_vars)%unit=unit !!$ allvars(n_vars)%positive=positive !!$ allvars(n_vars)%standard_name=standardname !!$ allvars(n_vars)%data2D=data2D !!$ allvars(n_vars)%data3D=data3D !!$ allvars(n_vars)%budgetdata=data3D !!$ allvars(n_vars)varid=-1 !!$ allvars(n_vars)%filenunit=-1 !!$ allvars(n_vars)%filename=filename1 !!$ allvars(n_vars)%dimensions=2 !!$ allvars(n_vars)%lon_varid=-1 !!$ allvars(n_vars)%lat_varid=-1 !!$ allvars(n_vars)%lev_varid=-1 !!$ allvars(n_vars)%time_varid=-1 !!$ allvars(n_vars)%bounds_varid=-1 !!$ allvars(n_vars)%dims=dims !!$ allvars(n_vars)%class=class !!$ allvars(n_vars)%indices=(/0,0,0,0,0,0,0/)) !!$ allvars(n_vars)%xmgas=molarweight !!$ allvars(n_vars)%table_id= ! 3D variables else if (data_dims==3) then ! allocate holders for data allocate(allvars(n_vars)%data2D(1,1)) if (trim(shortname)=='phalf') then allocate(allvars(n_vars)%data3D(i1:i2,j1:j2,1:lmr+1)) allocate(data3D(i1:i2,j1:j2,1:lmr+1)) else allocate(allvars(n_vars)%data3D(i1:i2,j1:j2,1:lmr)) allocate(data3D(i1:i2,j1:j2,1:lmr)) end if ! allocate local variables ! maybe remove these allocate(data2D(1,1)) !allocate(data3D(i1:i2,j1:j2,1:lmr)) ! zero local data holders data2D(:,:)=0.0 data3D(:,:,:)=0.0 !init variable if (trim(class)=='crescendo')then ! filename1= trim(outdir)//'/'//trim(shortname)//'_'//trim(table_id)//'_'//trim(source_id)//'_'//trim(experiment_id)//'_'//trim(member_id)//'_'//trim(grid_label)//'_'//trim(time_range)//trim('.nc') filename1= trim(outdir)//'/'//trim(shortname)//'_'//trim(class)//'_'//trim(table_id)//'_'//trim(source_id)//'_'//trim(experiment_id)//'_'//trim(member_id)//'_'//trim(grid_label)//trim(time_range)//trim('.nc') else filename1= trim(outdir)//'/'//trim(shortname)//'_'//trim(table_id)//'_'//trim(source_id)//'_'//trim(experiment_id)//'_'//trim(member_id)//'_'//trim(grid_label)//trim(time_range)//trim('.nc') end if allvars(n_vars)=varfile(itm5,shortname,longname,unit,positive,standardname,data2D,data3D,budget,-1,-1,& filename1,3,-1,-1,-1,-1,-1,-1,-1,-1,-1,data_dims,freq,class,(/0,0,0,0,0,0,0,0,0,0,0/),pxmgas,table) end if ! add chemical info also: (vars beginning with emi,wet,dry) select case (trim(shortname(4:))) case ('so4') allvars(n_vars)%indices(1:7)=(/iso4nus,iso4ais,iso4acs,iso4cos,0,0,0/)!mode_end_so4 case ('so2') allvars(n_vars)%indices(1)=iso2 case ('oa','emioa') !allvars(n_vars)%indices(1:7)=(/ipomais,ipomacs,ipomcos,ipomaii,0,0,0/)!mode_end_pom allvars(n_vars)%indices(1:9)=(/ipomais,ipomacs,ipomcos,ipomaii,isoanus,isoaais,isoaacs,isoacos,isoaaii/)!mode_end_pom case ('poa','emipoa') allvars(n_vars)%indices(1:4)=(/ipomais,ipomacs,ipomcos,ipomaii/)!mode_end_pom case ('soa') allvars(n_vars)%indices(1:7)=(/isoanus,isoaais,isoaacs,isoacos,isoaaii,0,0/)!mode_end_pom case ('bc') allvars(n_vars)%indices(1:7)=(/ibcaii,ibcais,ibcacs,ibccos,0,0,0/)!mode_end_bc case ('ss','emiss') allvars(n_vars)%indices(1:7)=(/issacs,isscos,0,0,0,0,0/)!mode_end_ss case ('dust') allvars(n_vars)%indices(1:7)=(/iduacs,iducos,iduaci,iducoi,0,0,0/)!mode_end_dust case('nox') !allvars(n_vars)%indices(1:2)=(/ino,ino2/) allvars(n_vars)%indices(1)=inox case('voc') allvars(n_vars)%indices(1)=1!(/ipar,ieth,iole,iald2,imgly,0,0/) case('bvoc') allvars(n_vars)%indices(1:2)=(/iterp,iisop/)!ibvoc case('nh3','sfnh3') allvars(n_vars)%indices(1)=inh3 case('nh4') allvars(n_vars)%indices(1)=inh4 case('noy') allvars(n_vars)%indices(1:11)=(/ ino, iNO2, ino3, iHNO3, iNO3_a, ihno4, in2o5, iPAN, iOrgNtr, ihono, ich3o2no2/) !allvars(n_vars)%indices(1)=ino2!(/ ino, iNO2, ino3, iHNO3, iNO3_a, ihno4, in2o5, iPAN, iOrgNtr, ihono, ich3o2no2/) case('no3') allvars(n_vars)%indices(1)=ino3_a case('pm1') allvars(n_vars)%indices(1)=-1 case('pm2p5','sfpm25') allvars(n_vars)%indices(1)=-1 case('pm10') allvars(n_vars)%indices(1)=-1 case('o3') allvars(n_vars)%indices(1)=io3 case('co') allvars(n_vars)%indices(1)=ico case('dms') allvars(n_vars)%indices(1)=idms case('isop') allvars(n_vars)%indices(1)=iisop !case('terp') ! allvars(n_vars)%indices(1)=iterp end select select case (trim(shortname)) case('noy') allvars(n_vars)%indices(1:11)=(/ iNO, ino2, ino3, iHNO3, iNO3_a, ihno4, in2o5, iPAN, iOrgNtr, ihono, ich3o2no2/) !allvars(n_vars)%indices(1)=ino2!(/ iNO, ino3, iHNO3, iNO3_a, ihno4, in2o5, iPAN, iOrgNtr, ihono, ich3o2no2/) case('areacella') allvars(n_vars)%indices(:)=0 areacella=n_vars case('c2h2') allvars(n_vars)%indices(1)=-1 case('c2h6') allvars(n_vars)%indices(1)=ic2h6 case('c3h6') allvars(n_vars)%indices(1)=ic3h6 case('c3h8') allvars(n_vars)%indices(1)=ic3h8 case('ch3coch3') allvars(n_vars)%indices(1)=iacet!ich3coch3 case('ch4') allvars(n_vars)%indices(1)=ich4 case('co') allvars(n_vars)%indices(1)=ico case('so2') allvars(n_vars)%indices(1)=iso2 case('co2') allvars(n_vars)%indices(1)=-1 case('dms') allvars(n_vars)%indices(1)=idms case('h2o') allvars(n_vars)%indices(1)=-1!ih2o case('hcho') allvars(n_vars)%indices(1)=ich2o case('hcl') allvars(n_vars)%indices(1)=-1!ihcl case('hno3','sfmmrhno3') allvars(n_vars)%indices(1)=ihno3 case('isop') allvars(n_vars)%indices(1)=iisop case('n2o') allvars(n_vars)%indices(1)=-1!in2o case('no', 'sfno') allvars(n_vars)%indices(1)=ino case('no2','sfno2') allvars(n_vars)%indices(1)=ino2 case('o3','sfo3') allvars(n_vars)%indices(1)=io3 case('o3ste') allvars(n_vars)%indices(1)=io3s case('oh') allvars(n_vars)%indices(1)=ioh case('pan') allvars(n_vars)%indices(1)=ipan !CRESCENDO variables case('nh3','sfmmrnh3','sfnh3') allvars(n_vars)%indices(1)=inh3 case('nh4','sfmmrnh4') allvars(n_vars)%indices(1)=inh4 case('sfmmrno3') allvars(n_vars)%indices(1)=ino3_a case('pm1') allvars(n_vars)%indices(1)=-1 case('pm2p5','sfpm25') allvars(n_vars)%indices(1)=-1 case ('mmraerh2o_1','mmraerh2o_2','mmraerh2o_3','mmraerh2o_4','mmraerh2o') allvars(n_vars)%indices(1)=-1 case ('sfmmrso4') allvars(n_vars)%indices(1:4)=(/iso4nus,iso4ais,iso4acs,iso4cos/) case ('sfmmrsoa') allvars(n_vars)%indices(1:5)=(/isoanus,isoaais,isoaacs,isoacos,isoaaii/) case ('sfmmroa') allvars(n_vars)%indices(1:4)=(/ipomais,ipomacs,ipomcos,ipomaii/) case ('sfmmrbc') allvars(n_vars)%indices(1:4)=(/ibcais,ibcacs,ibccos,ibcaii/) case ('sfmmrdustbel1','sfmmrdustabv1','sfmmrdustabv10','seddust','sfmmrdust') allvars(n_vars)%indices(1:4)=(/iduacs,iducos,iduaci,iducoi/) case ('sfmmrss') allvars(n_vars)%indices(1:2)=(/issacs,isscos/) case ('mmrsu1','sfmmrsu1') allvars(n_vars)%indices(1)=iso4nus case ('mmrsu2','sfmmrsu2') allvars(n_vars)%indices(1)=iso4ais case ('mmrsu3','sfmmrsu3') allvars(n_vars)%indices(1)=iso4acs case ('mmrsu4','sfmmrsu4') allvars(n_vars)%indices(1)=iso4cos case ('mmrsoa1','sfmmrsoa1') allvars(n_vars)%indices(1)=isoanus case ('mmrsoa2','sfmmrsoa2') allvars(n_vars)%indices(1)=isoaais case ('mmrsoa3','sfmmrsoa3') allvars(n_vars)%indices(1)=isoaacs case ('mmrsoa4','sfmmrsoa4') allvars(n_vars)%indices(1)=isoacos case ('mmrsoa5','sfmmrsoa5') allvars(n_vars)%indices(1)=isoaaii case ('mmroa2','sfmmroa2') allvars(n_vars)%indices(1)=ipomais case ('mmroa3','sfmmroa3') allvars(n_vars)%indices(1)=ipomacs case ('mmroa4','sfmmroa4') allvars(n_vars)%indices(1)=ipomcos case ('mmroa5','sfmmroa5') allvars(n_vars)%indices(1)=ipomaii case ('mmrbc2','sfmmrbc2') allvars(n_vars)%indices(1)=ibcais case ('mmrbc3','sfmmrbc3') allvars(n_vars)%indices(1)=ibcacs case ('mmrbc4','sfmmrbc4') allvars(n_vars)%indices(1)=ibccos case ('mmrbc5','sfmmrbc5') allvars(n_vars)%indices(1)=ibcaii case ('mmrss3','sfmmrss3') allvars(n_vars)%indices(1)=issacs case ('mmrss4','sfmmrss4') allvars(n_vars)%indices(1)=isscos case ('mmrdu3','sfmmrdu3') allvars(n_vars)%indices(1)=iduacs case ('mmrdu4','sfmmrdu4') allvars(n_vars)%indices(1)=iducos case ('mmrdu6','sfmmrdu6') allvars(n_vars)%indices(1)=iduaci case ('mmrdu7','sfmmrdu7') allvars(n_vars)%indices(1)=iducoi case ('sfndtot','sfnd100','ndtot') allvars(n_vars)%indices(1:7)=(/inus_n,iais_n,iacs_n,icos_n,iaii_n,iaci_n,icoi_n/) case ('nd1','sfnd1','sfnd1_tst') allvars(n_vars)%indices(1)=inus_n case ('nd2','sfnd2') allvars(n_vars)%indices(1)=iais_n case ('nd3','sfnd3') allvars(n_vars)%indices(1)=iacs_n case ('nd4','sfnd4') allvars(n_vars)%indices(1)=icos_n case ('nd5','sfnd5') allvars(n_vars)%indices(1)=iaii_n case ('nd6','sfnd6') allvars(n_vars)%indices(1)=iaci_n case ('nd7','sfnd7') allvars(n_vars)%indices(1)=icoi_n ! case('sfmmrnh3') ! allvars(n_vars)%indices(1)=inh3 ! case('sfmmrnh4') ! allvars(n_vars)%indices(1)=inh4 case('drydms') allvars(n_vars)%indices(1)=idms case('wetdms') allvars(n_vars)%indices(1)=idms case('dryno3') allvars(n_vars)%indices(1)=ino3_a case('wetno3') allvars(n_vars)%indices(1)=ino3_a case('dryhno3') allvars(n_vars)%indices(1)=ihno3 case('wethno3') allvars(n_vars)%indices(1)=ihno3 case('dryno2') allvars(n_vars)%indices(1)=ino2 case('dryno') allvars(n_vars)%indices(1)=ino case('drypan') allvars(n_vars)%indices(1)=ipan case('sfo3max') allvars(n_vars)%indices(1)=io3 case('mono') allvars(n_vars)%indices(1)=iterp case('co2mass') allvars(n_vars)%indices(1)=-1! not available case('chepsoa','chepsoa3d','chepsoa2d') allvars(n_vars)%indices(1)=-1 case('flashrate') allvars(n_vars)%indices(1)=-1 case('md') allvars(n_vars)%indices(1)=-1 case('sedustCI') allvars(n_vars)%indices(1)=iducoi case('depdust') allvars(n_vars)%indices(1:4)=(/iduacs,iducos,iduaci,iducoi/) case('concdust') allvars(n_vars)%indices(1:4)=(/iduacs,iducos,iduaci,iducoi/) case('conccn') allvars(n_vars)%indices(1:7)=(/inus_n,iais_n,iacs_n,icos_n,iaii_n,iaci_n,icoi_n/) case('sconcss') allvars(n_vars)%indices(1:2)=(/issacs,isscos/) case('loaddust') allvars(n_vars)%indices(1:4)=(/iduacs,iducos,iduaci,iducoi/) case('lossch4') allvars(n_vars)%indices(1)=ich4 case('lossco') allvars(n_vars)%indices(1)=ico end select if (class=='crescendo') then select case (trim(shortname)) case('od440dust') allvars(n_vars)%indices(1)=-1 if (freq=='day' .and. data_dims==2)then od440dustday=n_vars end if case('od440aer') allvars(n_vars)%indices(1)=-1 if (freq=='1hr') then od440aerhr=n_vars else if (freq=='day') then od440aerdaily=n_vars else if (freq=='mon') then od440aermonthly=n_vars end if case('od550aer') allvars(n_vars)%indices(1)=-1 if (freq=='1hr')then od550aerhr=n_vars end if case('od550dust') if (freq=='day' .and. data_dims==2)then od550dust2dday=n_vars end if case('od5503ddust') if (freq=='day' .and. data_dims==3)then od550dust3dday=n_vars end if case('ec550aer') allvars(n_vars)%indices(1)=-1 ec550aermon=n_vars end select end if if (class=='sf1d') then select case (trim(shortname)) case('od550aer') allvars(n_vars)%indices(1)=-1 od550aerdaily=n_vars case('sfo3max') allvars(n_vars)%indices(1)=io3 end select end if if(class=='optics') then select case (trim(shortname)) case('ec550aer') allvars(n_vars)%indices(1)=-1 ec550aer=n_vars case('od550aer') allvars(n_vars)%indices(1)=-1 od550aer=n_vars case('od550so4') allvars(n_vars)%indices(1)=-1 od550so4=n_vars case('od550bc') allvars(n_vars)%indices(1)=-1 od550bc=n_vars case('od550oa') allvars(n_vars)%indices(1)=-1 od550oa=n_vars case('od550soa') allvars(n_vars)%indices(1)=-1 od550soa=n_vars case('od550ss') allvars(n_vars)%indices(1)=-1 od550ss=n_vars case('od550dust') allvars(n_vars)%indices(1)=-1 od550dust=n_vars case('od550no3') allvars(n_vars)%indices(1)=-1 od550no3=n_vars case('od550aerh2o') allvars(n_vars)%indices(1)=-1 od550aerh2o=n_vars case('od550lt1aer') allvars(n_vars)%indices(1)=-1 od550lt1aer=n_vars case('abs550aer') allvars(n_vars)%indices(1)=-1 abs550aer=n_vars case('od440aer') allvars(n_vars)%indices(1)=-1 if (freq=='mon') then od440aer=n_vars !else if (freq=='day') then ! od440aerdaily=n_vars end if case('od350aer') allvars(n_vars)%indices(1)=-1 od350aer=n_vars case('od870aer') allvars(n_vars)%indices(1)=-1 od870aer=n_vars end select end if call goLabel() status = 0 call GO_Timer_End( itim_addvar, status ) IF_NOTOK_RETURN(status=1) end subroutine add_variable subroutine write_attributes(status,j_var) implicit none integer :: j_var integer,intent(out)::status character(len=*), parameter :: rname = mname//'/output_aerchemmip_writeattr' call goLabel(rname) call GO_Timer_Start( itim_attr, status ) IF_NOTOK_RETURN(status=1) call MDF_Put_Att( allvars(j_var)%fileunit, MDF_GLOBAL, 'title', 'Model output for Aerchemmip', status ) IF_NOTOK_MDF(fid= allvars(j_var)%fileunit) call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lon_varid , 'units', 'degrees_east', status) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lon_varid , 'axis', 'X', status) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lon_varid , 'long_name', 'longitude', status) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lon_varid , 'standard_name', 'longitude', status) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lat_varid , 'units', 'degrees_north', status) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lat_varid , 'axis', 'Y', status) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lat_varid , 'long_name', 'latitude', status) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lat_varid , 'standard_name', 'latitude', status) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) ! allvars(j_var)%lev_varid if (allvars(j_var)%dims==3) then call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lev_varid , 'units', 'level', status) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lev_varid , 'axis', 'Z', status) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lev_varid , 'long_name', 'hybrid model level at layer midpoints', status) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lev_varid , 'standard_name', 'hybrid_model_level', status) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lev_varid , 'formula', 'a+b*ps', status) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lev_varid , 'positive', 'up', status) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%hyam_varid , 'long_name', 'hybrid A coefficient at layer midpoints', status) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%hyam_varid , 'units', 'Pa', status) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%hybm_varid , 'long_name', 'hybrid B coefficient at layer midpoints', status) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%hybm_varid , 'units', '1', status) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%hyai_varid , 'long_name', 'hybrid A coefficient at layer interfaces', status) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%hyai_varid , 'units', 'Pa', status) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%hybi_varid , 'long_name', 'hybrid B coefficient at layer interfaces', status) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%hybi_varid , 'units', '1', status) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) END if if (allvars(j_var)%table_id/='AERfx')then call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%time_varid , 'units', 'days since 1750-01-01 00:00:00', status) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) ! call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%time_varid , 'axis', 'X', status) ! IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%time_varid , 'calendar', 'gregorian', status) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%time_varid , 'long_name', 'time', status) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%time_varid , 'axis', 'T', status) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) !time bounds call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%time_varid , 'bounds', 'time_bounds', status) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) end if !experiment= !CMIP6 global attributes: call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL , 'Conventions', 'CF-1.7 CMIP-6.0 UGRID-0.9', status) call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL , 'activity_id', 'AerChemMIP', status) call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL , 'branch_method','', status) call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL , 'creation_date','', status) call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL , 'data_specs_version','1.0.0', status) call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL , 'experiment',trim(experiment), status) call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL , 'experiment_id',trim(experiment_id), status) call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL , 'forcing_index',trim(forcing_i), status) call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL , 'frequency',trim(allvars(j_var)%freq), status) call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL , 'further_info_url','MISSING', status) call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL , 'grid','native '//trim(grid)//' degree grid', status)!module variables call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL , 'grid_label',trim(grid_label), status)!module variables call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL , 'initialization_index',trim(initialization_i), status) call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL , 'institution','KNMI, The Netherlands; SMHI, Sweden; DMI, Denmark; & &AEMET, Spain; Met Eireann, Ireland; CNR-ISAC, Italy; Instituto de Meteorologia, Portugal; FMI, Finland; BSC, Spain; & &Centro de Geofisica, University of Lisbon, Portugal; ENEA, Italy; Geomar, Germany; Geophysical Institute, University of Bergen, Norway; & &ICHEC, Ireland; ICTP, Italy; IMAU, The Netherlands; IRV, Sweden; Lund University, Sweden; & &Meteorologiska Institutionen, Stockholms University, Sweden; Niels Bohr Institute, University of Copenhagen, Denmark; & &NTNU, Norway; SARA, The Netherlands; Unite ASTR, Belgium; Universiteit Utrecht, The Netherlands; & &Universiteit Wageningen, The Netherlands; University College Dublin, Ireland; Vrije Universiteit Amsterdam, the Netherlands; & &University of Helsinki, Finland; KIT, Karlsruhe, Germany; USC, University of Santiago de Compostela, Spain; & &Uppsala Universitet, Sweden; NLeSC, Netherlands eScience Center, The Netherlands', status) call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL , 'institution_id','EC-Earth-Consortium', status) call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL , 'license','NEEDS DEFINING', status) call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL , 'mip_era','CMIP6', status) call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL , 'nominal_resolution','250 km', status)!dmax !dmax=r*phi/2*(1+((phi**2+lamb**2)/(phi*lamb))*np.arctan(lamb/phi))=348 r=6371, phi=2(lat), lamb=3(long) !CMIP6 global attributes: 100 < dmax < 350 -> nominal resolution 250 km call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL , 'physics_index',trim(physics_i), status) call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL , 'product','output', status)!only choice call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL , 'realization_index',trim(realization_i), status)!1 for primary or single realization call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL , 'realm',trim(realm), status)! depends on run, some are AGCM call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL , 'source',trim(source), status)! call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL , 'source_id',trim(source_id), status) call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL , 'source_type',trim(source_type), status) call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL , 'sub_experiment','', status) call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL , 'sub_experiment_id','', status) call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL , 'table_id',trim(allvars(j_var)%table_id), status) call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL , 'tracking_id','', status) call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL , 'variable_id',trim(allvars(j_var)%vname), status) call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL , 'variant_label','', status) call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL , 'variant_label','', status) call GO_Timer_End( itim_attr, status ) IF_NOTOK_RETURN(status=1) call goLabel() status = 0 end subroutine write_attributes subroutine write_dimensions(status,j_var) use dims, only : iyear0 !current year use go_date, only : days_in_year,newDate use partools , only : isRoot,myid implicit none integer :: j_var integer,intent(out)::status integer :: i1,i2,j1,j2,imr,jmr,lmr integer :: lon_varid,lonid,lon_dimid integer :: lat_varid,latid,lat_dimid integer :: lev_varid,levid,lev_dimid integer :: hym_varid,hym_dimid integer :: hyi_varid,hyi_dimid integer :: time_varid,timeid,time_dimid,bounds_dimid,bounds_varid,boudid ! most of data is monthly mean, but change to dynamic number of output steps needed integer :: nout_steps!=12 integer :: nhym integer :: nhyi character(len=*), parameter :: rname = mname//'/output_aerchemmip_write_dim' call goLabel(rname) ! define dimensions !call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) imr=dimension_data%nlon jmr=dimension_data%nlat lmr=dimension_data%nlev nhym=lmr nhyi=lmr+1 ! With parallel netcdf whole netcdf must be reserved at the time of initialisation ! therefore we need to know the number of output steps per file. ! Define number of output steps in a file depending on the output frequency ! use newdate to create TDate structure, and use that in days_in_year if (allvars(j_var)%table_id=='AERhr')then nout_steps=24*days_in_year(newdate(iyear0)) else if (allvars(j_var)%table_id=='AER6hr')then nout_steps=4*days_in_year(newdate(iyear0)) else if (allvars(j_var)%table_id=='AERday')then nout_steps=days_in_year(newdate(iyear0)) else if (allvars(j_var)%table_id=='AERmon'.or. (allvars(j_var)%table_id=='AERmonZ'))then nout_steps=12 end if if (isroot) then !DEFINE DIMENSIONS call MDF_Def_Dim( allvars(j_var)%fileunit, 'lon', imr,lon_dimid, status ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) call MDF_Def_Dim( allvars(j_var)%fileunit, 'lat', jmr, lat_dimid, status ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) if (allvars(j_var)%dims==3) then if (trim(allvars(j_var)%vname)=='phalf') then call MDF_Def_Dim( allvars(j_var)%fileunit, 'lev', lmr+1, lev_dimid, status ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) else call MDF_Def_Dim( allvars(j_var)%fileunit, 'lev', lmr, lev_dimid, status ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) end if end if if (allvars(j_var)%table_id/='AERfx')then !call MDF_Def_Dim( allvars(j_var)%fileunit, 'time', nout_steps, time_dimid, status ) call MDF_Def_Dim( allvars(j_var)%fileunit, 'time', MDF_UNLIMITED, time_dimid, status ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) call MDF_Def_Dim( allvars(j_var)%fileunit, 'bounds', 2, bounds_dimid, status ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) end if if (allvars(j_var)%dims==3) then call MDF_Def_Dim( allvars(j_var)%fileunit, 'nhym', nhym, hym_dimid, status ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) call MDF_Def_Dim( allvars(j_var)%fileunit, 'nhyi', nhyi, hyi_dimid, status ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) end if ! define dimension variables ! longitude call MDF_Def_Var( allvars(j_var)%fileunit, 'lon', MDF_DOUBLE, & (/ lon_dimid/), allvars(j_var)%lon_varid, status ) ! define latitude variable call MDF_Def_Var( allvars(j_var)%fileunit, 'lat', MDF_DOUBLE, & (/ lat_dimid/), allvars(j_var)%lat_varid, status ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) ! level if (allvars(j_var)%dims==3) then call MDF_Def_Var( allvars(j_var)%fileunit, 'lev', MDF_DOUBLE, & (/ lev_dimid/), allvars(j_var)%lev_varid, status ) end if if (allvars(j_var)%table_id/='AERfx')then call MDF_Def_Var( allvars(j_var)%fileunit, 'time', MDF_DOUBLE, & (/ time_dimid/), allvars(j_var)%time_varid, status ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) call MDF_Def_Var( allvars(j_var)%fileunit, 'time_bounds', MDF_DOUBLE, & (/ bounds_dimid,time_dimid/), allvars(j_var)%bounds_varid, status ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) end if if (allvars(j_var)%dims==3) then ! define hybm variable call MDF_Def_Var( allvars(j_var)%fileunit, 'hybm', MDF_DOUBLE, & (/ hym_dimid/), allvars(j_var)%hybm_varid, status ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) ! define hyam variable call MDF_Def_Var( allvars(j_var)%fileunit, 'hyam', MDF_DOUBLE, & (/ hym_dimid/), allvars(j_var)%hyam_varid, status ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) ! define hybi variable call MDF_Def_Var( allvars(j_var)%fileunit, 'hybi', MDF_DOUBLE, & (/ hyi_dimid/), allvars(j_var)%hybi_varid, status ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) ! define hyai variable call MDF_Def_Var( allvars(j_var)%fileunit, 'hyai', MDF_DOUBLE, & (/ hyi_dimid/), allvars(j_var)%hyai_varid, status ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) end if ! * Write out dimension variable values ! Write out hybm if (allvars(j_var)%dims==3) then ! midpoints if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit, allvars(j_var)%hybm_varid,levi%fb, status) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) ! Write out hyam if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit, allvars(j_var)%hyam_varid,levi%fa, status) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) ! interfaces if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit, allvars(j_var)%hybi_varid,levi%b, status) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) ! Write out hyai if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit, allvars(j_var)%hyai_varid,levi%a, status) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) end if ! Write out the longitudes if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit, allvars(j_var)%lon_varid, dimension_data%lon, status) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) !write out the latitude to variable if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit, allvars(j_var)%lat_varid, dimension_data%lat, status) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) if (allvars(j_var)%dims==3) then if (trim(allvars(j_var)%vname)=='phalf') then ! Write out the levels if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit, allvars(j_var)%lev_varid, (/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,35/), status) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) else if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit, allvars(j_var)%lev_varid, (/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) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) end if end if ! Time will be written during output write -steps end if call goLabel() status = 0 end subroutine write_dimensions subroutine write_var(status,j_var) implicit none integer :: j_var integer,intent(out)::status integer :: i1,i2,j1,j2,imr,jmr,lmr integer :: lon_varid,lonid integer :: lat_varid,latid integer :: lev_varid,levid integer :: time_varid,timeid character(len=*), parameter :: rname = mname//'/output_aerchemmip_write_dim' call goLabel(rname) ! define dimensions !call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) imr=dimension_data%nlon jmr=dimension_data%nlat lmr=dimension_data%nlev ! define dimension variables ! longitude if (allvars(j_var)%dims==2.and.allvars(j_var)%table_id/='AERfx') then call MDF_Def_Var( allvars(j_var)%fileunit, allvars(j_var)%vname, MDF_DOUBLE, & (/allvars(j_var)%lon_varid, allvars(j_var)%lat_varid, allvars(j_var)%time_varid/), allvars(j_var)%varid, status ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) else if (allvars(j_var)%dims==2.and.allvars(j_var)%table_id=='AERfx') then call MDF_Def_Var( allvars(j_var)%fileunit, allvars(j_var)%vname, MDF_DOUBLE, & (/allvars(j_var)%lon_varid, allvars(j_var)%lat_varid/), allvars(j_var)%varid, status ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) else !(allvars(j_var)%dims==3) call MDF_Def_Var( allvars(j_var)%fileunit, allvars(j_var)%vname, MDF_DOUBLE, & (/allvars(j_var)%lon_varid, allvars(j_var)%lat_varid, allvars(j_var)%lev_varid,allvars(j_var)%time_varid/), allvars(j_var)%varid, status ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) end if ! Write out the longitudes call MDF_Put_Att( allvars(j_var)%fileunit, allvars(j_var)%varid, 'long_name', trim(allvars(j_var)%lname), status ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) call MDF_Put_Att(allvars(j_var)%fileunit,allvars(j_var)%varid, 'standard_name', trim(allvars(j_var)%standard_name), status ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) call MDF_Put_Att(allvars(j_var)%fileunit , allvars(j_var)%varid, 'units', trim(allvars(j_var)%unit), status ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) call MDF_EndDef( allvars(j_var)%fileunit, status ) IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) call goLabel() status = 0 end subroutine write_var subroutine calculate_optics( region, dhour, l_do_ec550aer_only,status) use optics, only : optics_aop_get use MeteoData, only : gph_dat,temper_dat USE toolbox, only : ltropo_ifs, lvlpress use dims, only : lm implicit none integer :: i, j, k, imr, jmr, lmr, lwl, dtime integer :: i1, i2, j1, j2,ivar integer :: ie, je ! indices for subdomain extended with halo cells integer, parameter :: l_halo=1 logical :: polar integer :: istat, imode real :: dens, load_tmp 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 :: tempdust2dday real, dimension(:,:), allocatable :: tempdust2d440day integer :: ncontr, lvec, lct, l, indoffset, nwl integer :: nadd, nadd_max, indoffset_stat, indabs integer :: iadd integer :: region,dhour,status real :: no3fac, wght, dwght, tabs real,dimension(:),allocatable :: tx integer :: itropo real, dimension(:,:,:), pointer :: gph ! height (incl. oro) real,dimension(:,:,:), pointer :: t ! temperature (K) logical :: l_do_ec550aer_only character(len=*), parameter :: rname = mname//'/output_aerchemmip_optics' 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 allocate(tx(lm(region))) t => temper_dat (region)%data gph => gph_dat (region)%data ! --------------------- ! 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 AEROCOM 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 ! Total(water), Total(nowater) !ncontr = 10 ncontr = 12 ! Total, SO4, BC, OC,SOA, SS, DU, NO3, Water, lt1aer dtime=dhour*3600 !TB Following additional fields are nnecessary but call to optics needs it... ! 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, .false., 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/) ) 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_aerocom: 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 (.not. l_do_ec550aer_only)then ! if split and if wavelength=550nm if( wdep_out(lwl)%split ) then if (wdep_out(lwl)%wl == 0.550) then ! for 550nm: ! 1) AODs ! 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)) allocate(tempdust2dday(i1:i2,j1:j2)) do lct = 1, ncontr temp2d = 0.0 tempdust2dday=0.0 do j = j1, j2 do i = i1, i2 ! multiply with height elements and sum up tabs = 0.0 ! AOD output will only be for troposphere in data request tx(:)=t(i,j,:) itropo = ltropo_ifs(region,i,j,tx,lm(region)) do l = 1, itropo!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 ! TOTAL AOD ! 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)) else if (lct==4) then ! OAAOD ! add SOA aod (5)to POM aod (4) (AerChemMIP expects total oa in aod550oa) temp2d(i,j) = temp2d(i,j) + opt_ext(i,j,l,lwl,5) * & (gph_dat(region)%data(i,j,l+1) - & gph_dat(region)%data(i,j,l )) else if (lct==7) then ! DUSTAOD if ( wdep_out(lwl)%wl == 0.550) then if (crescendo_out) allvars(od550dust3dday)%data3D(i,j,l)= opt_ext(i,j,l,lwl,lct) * & (gph_dat(region)%data(i,j,l+1) - & gph_dat(region)%data(i,j,l )) tempdust2dday(i,j)=tempdust2dday(i,j)+ opt_ext(i,j,l,lwl,lct) * & (gph_dat(region)%data(i,j,l+1) - & gph_dat(region)%data(i,j,l )) else tempdust2dday(i,j)=tempdust2dday(i,j)+ opt_ext(i,j,l,lwl,lct) * & (gph_dat(region)%data(i,j,l+1) - & gph_dat(region)%data(i,j,l )) end if end if end do ! Absorption: do this only once for the totals if( lct == 1 ) then allvars(abs550aer)%data2D(i,j)=allvars(abs550aer)%data2D(i,j) + tabs*dtime !extra output for total od550aer (ncontr==1) allvars(od550aerdaily)%data2D(i,j)=allvars(od550aerdaily)%data2D(i,j) + temp2d(i,j)*dtime if (crescendo_out)then allvars(od550aerhr)%data2D(i,j)= dtime*temp2d(i,j) end if allvars(od550aer)%data2D(i,j)=allvars(od550aer)%data2D(i,j) + temp2d(i,j)*dtime else if (lct<11)Then !AOD by compounds allvars(indoffset+lct-1)%data2D(i,j)=allvars(indoffset+lct-1)%data2D(i,j) + temp2d(i,j)*dtime if (crescendo_out .and. lct==7 .and. wdep_out(lwl)%wl == 0.550) then !DUSTAOD for 550nm allvars(od550dust2dday)%data2D(i,j)=allvars(od550dust2dday)%data2D(i,j) + tempdust2dday(i,j)*dtime end if end if end do end do end do deallocate( temp2d ) deallocate( tempdust2dday) !if 440 has splits do this else if (wdep_out(lwl)%wl == 0.440 ) then indoffset = od440aer !not needed for AERCHEMMIP indabs = -1 !abs440aer ! fill AODs allocate(tempdust2d440day(i1:i2,j1:j2)) allocate(temp2d(i1:i2,j1:j2)) tempdust2d440day=0.0 temp2d = 0.0 tempdust2d440day=0.0 do j = j1, j2 do i = i1, i2 ! AOD output will only be for troposphere in data request tx(:)=t(i,j,:) itropo = ltropo_ifs(region,i,j,tx,lm(region)) ! multiply with height elements and sum up tabs = 0.0 do l = 1, itropo!lmr ! od440aer lct=1 ! total aod 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 )) 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)) !OD440DUST lct=7 tempdust2d440day(i,j)=tempdust2d440day(i,j)+ opt_ext(i,j,l,lwl,lct) * & (gph_dat(region)%data(i,j,l+1) - & gph_dat(region)%data(i,j,l )) end do if (indabs>0) then allvars(indabs)%data2D(i,j)=allvars(indabs)%data2D(i,j) + tabs*dtime end if end do end do allvars(od440aer)%data2D=allvars(od440aer)%data2D + temp2d*dtime if (crescendo_out)then allvars(od440dustday)%data2D=allvars(od440dustday)%data2D + tempdust2d440day*dtime allvars(od440aerdaily)%data2D=allvars(od440aerdaily)%data2D + temp2d*dtime allvars(od440aerhr)%data2D=temp2d*dtime end if deallocate( tempdust2d440day) deallocate( temp2d ) end if else !NON SPLITS ! 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.870 ) then indoffset = od870aer !not needed for AERCHEMMIP indabs = -1 !abs870aer elseif ( wdep_out(lwl)%wl == 0.440 ) then indoffset = -1 !od440aer !not needed for AERCHEMMIP indabs = -1 !abs350aer elseif ( wdep_out(lwl)%wl == 0.350 ) then indoffset = od350aer !not needed for AERCHEMMIP indabs = -1 !abs350aer else write(gol,*) 'user_output_aerchemmip: wrong wavelength given for aerocom diagnostics' ; call goErr TRACEBACK status=1; return end if ! fill AODs, total only 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 if (indabs>0) then allvars(indabs)%data2D(i,j)=allvars(indabs)%data2D(i,j) + tabs*dtime end if end do end do if ((indoffset)==od870aer) then allvars(od870aer)%data2D=allvars(od870aer)%data2D + temp2d*dtime end if deallocate(temp2d) endif end if ! Extinction in 3D if ( wdep_out(lwl)%wl == 0.550 ) then allvars(ec550aer)%data3D= opt_ext(:,:,:,lwl,1) if (crescendo_out) allvars(ec550aermon)%data3D=allvars(ec550aermon)%data3D + opt_ext(:,:,:,lwl,1)*dtime endif end do deallocate( opt_ext, opt_a, opt_g, opt_add ) deallocate( tx ) call goLabel() ; status = 0 end subroutine calculate_optics real function mode_fraction(r,limit,imode) result(modfrac) Use mo_aero_m7, only : nmod, cmedr2mmedr, sigma implicit none !real :: var real :: r real :: z real :: limit ! real :: sigma real :: hr2=0.5*sqrt(2.0) integer::imode z=0.0 if( r<100*tiny(1.0))then modfrac=1.0 else z=(log(limit)-log(r*cmedr2mmedr(imode)))/log(sigma(imode)) modfrac=0.5+0.5*erf(z*hr2) end if end function mode_fraction end MODULE USER_OUTPUT_AERCHEMMIP