#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_AEROCOM ! ! !DESCRIPTION: ! ! This module provides objects and methods for the AEROCOM-2 related ! diagnostics package. It is strongly related to the setup of the ! optics module. ! A list of parameters are defined for which daily averages are output ! (so far), result is one file per day containing 50 2d-fields and ! 4 3d-fields. ! ! Life cycle(s): ! 1) output_aerocom_init: ! - if( newsrun ): - allocation and initialisation of output container ! (mixf%..., drydepos/wetdepos/emission, wdep_out) ! - initialisation of optics module ! - lookuptable ! - wavelengths ! - open output file (for today), initialise it and write attributes ! 2) output_aerocom_step: ! - gather tracer ! - add current state of model to containers ! (weighted by time since last call, --> period of validity) ! - do fluxes, surface concentrations and loads ! - do optics ! - aerocom_aopio_init: initialise and fill input parameters ! - optics::optics_calculate_aop: ! calculate optical properties (per wavelength) ! (take into account additional splitting, facilitated in ! optics::optics_calculate_aop) ! - fill containers for optical properties ! - free optics input parameters ! - fill aerocom output containers ! 3) output_aerocom_write: ! - collect_budgets: ! - use current budget containers (budemi, buddep, buddry) ! - compare with saved values and infer by that the increment ! - add to temporary output containers (wetdepos,drydepos,emission) ! - write output containers to file ! - close file ! 4) output_aerocom_done: ! - free allocated containers ! (mixf%..., drydepos/wetdepos/emission, wdep_out) ! !\\ !\\ ! !INTERFACE: ! MODULE USER_OUTPUT_AEROCOM ! ! !USES: ! use go, only : gol, goErr, goPr, goLabel use dims, only : nregions use optics, only : wavelendep use meteodata, only : global_lli, levi use MDF use TM5_DISTGRID, only : dgrid, Get_DistGrid, update_halo, update_halo_iband implicit none private ! ! !PUBLIC MEMBER FUNCTIONS: ! public :: output_aerocom_init public :: output_aerocom_step public :: output_aerocom_write public :: output_aerocom_done public :: wdep_out character(len=20), public :: aerocom_exper ! AeroCom experiment name character(len=20), public :: aerocom_freq ! AeroCom output frequency logical, parameter :: aai_output=.false. ! ! !PRIVATE DATA MEMBERS: ! character(len=*), parameter :: mname = 'user_output_aerocom' ! Changed naming convections to AeroCom Control 2015 experiment character(len=20), parameter :: f_dataid='aerocom3', f_model='TM5' character(len=20), parameter :: f_dimmix='global', f_dimstat='stations' ! ! !PRIVATE TYPES: ! type metafields integer :: itm5 character(len=16) :: vname character(len=64) :: lname character(len=10) :: unit character(len=10) :: positive character(len=130) :: standard_name end type metafields type metafields2 integer :: itm5 character(len=17) :: vname character(len=64) :: lname character(len=10) :: unit character(len=130) :: standard_name end type metafields2 type field3d type( metafields ) :: mf real, dimension(:,:,:), pointer :: field integer :: varid end type field3d type field2d type( metafields ) :: mf real, dimension(:,:), pointer :: field integer :: varid end type field2d type field1d type( metafields2 ) :: mf real, dimension(:,:), pointer :: field integer :: varid end type field1d type field0d type( metafields2 ) :: mf real, dimension(:), pointer :: field integer :: varid end type field0d type mixfile type( field3d ), dimension(:), pointer :: f3d type( field2d ), dimension(:), pointer :: f2d character(len=200) :: fname integer :: acct ! accumulation time integer :: funit ! unit number integer :: nlon ! x dimension of requested field integer :: nlat ! y dimension of requested field integer :: nlev ! z dimension of requested field end type mixfile type(mixfile), dimension(nregions), save :: mixf type stationfile type( field1d ), dimension(:), pointer :: f1d type( field0d ), dimension(:), pointer :: f0d character(len=200) :: fname integer :: acct ! accumulation time integer :: funit ! unit number integer :: nstat ! dimension of requested field integer :: nlev ! z dimension of requested field end type stationfile type(stationfile), dimension(nregions), save :: statf integer, dimension(5) :: statvarid type budgetstore real, dimension(:,:,:), allocatable :: f2dslast integer :: lasttime end type budgetstore type(budgetstore), dimension(nregions), save :: drydepos, wetdepos, emission ! wavelength information type(wavelendep), dimension(:), allocatable :: wdep_out integer, parameter :: statnamelen = 22 integer :: snamelendim type station integer :: statnb ! station number character(len=statnamelen) :: statname ! station name real :: statlat ! geographic coordinate real :: statlon ! geographic coordinate real :: statalt ! geographic coordinate end type station type(station), dimension(:), allocatable :: statlist character(len=statnamelen), dimension(:), allocatable :: tmpstatnames type statintp integer,dimension(:),pointer :: ixr, ixl ! coordinates in current grid integer,dimension(:),pointer :: jyr, jyl ! coordinates in current grid integer,dimension(:),pointer :: lzr, lzl ! coordinates in current grid real, dimension(:),pointer :: wixl, wixr ! interpolation coefficients real, dimension(:),pointer :: wjyl, wjyr ! interpolation coefficients real, dimension(:),pointer :: wlzl, wlzr ! interpolation coefficients real, dimension(:),pointer :: tght ! terrain height at station end type statintp type(statintp), dimension(nregions) :: sintp logical, dimension(:), allocatable :: stat_active ! true if station is in processor domain integer :: stat_halo_loc ! 1 if halo cells are needed for the local processor, else 0 integer :: stat_halo_glob ! 1 if halo cells are needed on any processor, else 0 ! the region dimension has already been removed for the above three variables integer :: fid ! file id for IF_NOTOK_MDF macro integer :: access_mode ! netcdf-4 parallel access mode for distributed data ! generate output if true: logical :: okdebug = .true. ! reference time: integer, parameter :: time_reftime6(6) = (/2001,01,01,00,00,00/) character(len=*), parameter :: time_units = 'days since 2001-01-01 00:00' ! sizes of arrays integer, save :: ntracer_3d, ntracer_2d integer, save :: ntracer_1dstat, ntracer_0dstat, nstations integer, save :: nextra_1dstat ! index pointers for 1d/2d/3d fields in mixf integer, save :: temp, hus, airmass, pressure integer, save :: ec5503Daer, abs5503Daer, ec3503Daer, abs3503Daer integer, save :: ps , precip , sconcoa , sconcbc , sconcso4 , sconcdust integer, save :: sconcss , sconcno3 , loadoa , loadbc , loadso4 , loaddust integer, save :: loadss , loadno3 , emioa , emibc , emiso2 , emiso4 integer, save :: emidust , emidms , emiss , dryso2 , drybc , dryso4 integer, save :: drydust , drydms , dryss , wetoa , wetbc , wetso2 integer, save :: wetso4 , wetdust , wetdms , wetss , od550aer , od550so4 integer, save :: od550bc , od550oa , od550ss , od550dust , od550no3 , od550aerh2o integer, save :: od550lt1aer , od550lt1dust, od550lt1ss , abs550aer , ec550aer , asyaer integer, save :: ec550dryaer , abs550dryaer, asydryaer , ec550drylt1aer, abs550drylt1aer integer, save :: od440aer , od870aer , sconcmsa , dryoa , sconcnh4 integer, save :: abs440aer , ec440dryaer , abs440dryaer integer, save :: abs870aer , ec870dryaer , abs870dryaer integer, save :: od350aer , abs350aer integer, save :: ec550dryaer_stat, abs550dryaer_stat, asydryaer_stat integer, save :: ec550drylt1aer_stat, abs550drylt1aer_stat integer, save :: ec440dryaer_stat, abs440dryaer_stat integer, save :: ec870dryaer_stat, abs870dryaer_stat integer, save :: bso4nus, bso4ais, bso4acs, bso4cos, bbcais , bbcacs integer, save :: bbccos , bbcaii , bpomais, bpomacs, bpomcos, bpomaii integer, save :: bssacs , bsscos , bduacs , bducos , bduaci , bducoi, bno3_a integer, save :: bnus_n , bais_n , bacs_n , bcos_n , baii_n , baci_n integer, save :: bcoi_n , bh2onus, bh2oais, bh2oacs, bh2ocos ! ! !REVISION HISTORY: ! 29 Nov 2010 - Achim Strunk - Creation ! 10 Oct 2014 - Twan van Noije - Adapted for TM5-mp ! ! !REMARKS: ! (1) compiled only if with_m7 is used. ! !EOP !------------------------------------------------------------------------ contains !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: OUTPUT_AEROCOM_INIT ! ! !DESCRIPTION: Initialise various parameters, eg., station lists, ! wavelengths for optics and output containers. ! This routine should be called once per day, since fields ! are daily averages and one file per day is written. ! Additional tasks are done for newsrun==.true. . !\\ !\\ ! !INTERFACE: ! subroutine output_aerocom_init(stat_output, status, iregion) ! ! !USES: ! !use MeteoData, only : sp_dat, set use chem_param, only : ntracet use chem_param, only : inus_n, iso4nus, nmod use chem_param, only : iais_n, iso4ais, ibcais, ipomais use chem_param, only : iacs_n, iso4acs, ibcacs, ipomacs, issacs, iduacs use chem_param, only : icos_n, iso4cos, ibccos, ipomcos, isscos, iducos use chem_param, only : iaii_n, ibcaii, ipomaii use chem_param, only : iaci_n, icoi_n, iduaci, iducoi use chem_param, only : ino3_a, inh4, imsa use dims, only : nregions, newsrun, idatee, idatei use dims, only : region_name, parent, itau use dims, only : xbeg, xend, ybeg, yend, dx, dy, xref, yref use dims, only : zbeg, zend, dz, zref use global_data, only : region_dat, outdir use datetime, only : tau2date, date2tau use budget_global, only : nbud_vg use partools, only : MPI_INFO_NULL, localComm use optics, only : Optics_Init use datetime, only : date2tau use dims, only : mlen ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !INPUT PARAMETERS: ! logical, intent(in) :: stat_output ! true if stations output is requested integer, intent(in), optional :: iregion ! ! !REVISION HISTORY: ! 29 Nov 2010 - Achim Strunk - ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/output_aerocom_init' ! --- local ------------------------------ integer :: imr, jmr, lmr, access_mode_sta integer :: region, varid integer :: io, istat, i, j, n, sc integer :: i1, i2, j1, j2 integer, dimension(6) :: idater character(len=10) :: idates character(len=16) :: lidates character(len=3) :: cwavel real :: reftime integer :: itaucur, itauref real :: dlat, dlon integer :: mlength ! --- begin ------------------------------- call goLabel(rname) ! access mode for distributed data (2d and 3d) #ifdef MPI #ifdef with_netcdf4_par access_mode = MDF_COLLECTIVE #else write(gol,'("AeroCom aerosol output requires netcdf4 with parallel access enabled")') ; call goErr TRACEBACK status=1; return #endif #else access_mode = MDF_INDEPENDENT #endif ! for station access_mode_sta = MDF_INDEPENDENT ! initialize (only once) if( newsrun ) then ! ensure that required meteo is loaded: ! call Set( sp_dat(region), status, used=.true. ) ! set wavelength information ! wl: wavelength in microns ! split: whether to split into contributions from ! M7 constituents (incl. water) allocate( wdep_out( 4 ) ) wdep_out(1)%wl = 0.550 ; wdep_out(1)%split = .true. ; wdep_out(1)%insitu = .true. wdep_out(2)%wl = 0.440 ; wdep_out(2)%split = .false. ; wdep_out(2)%insitu = .true. wdep_out(3)%wl = 0.870 ; wdep_out(3)%split = .false. ; wdep_out(3)%insitu = .true. wdep_out(4)%wl = 0.350 ; wdep_out(4)%split = .false. ; wdep_out(4)%insitu = .false. ! get the optics code prepared call Optics_Init(size(wdep_out), wdep_out, status ) IF_NOTOK_RETURN(status=1) ! ----------------------- ! parameters needed to reference the different 1d/2d/3d-fields ! (in order to avoid errors in referencing) ! finally, this list here simply determines the order in the output files ntracer_3d = 4 if (aai_output) then ntracer_3d = ntracer_3d + 4 endif ! ----------------------- temp = 1 ! 3d hus = 2 ! 3d airmass = 3 ! 3d pressure = 4 ! 3d ec5503Daer = 5 ! 3d abs5503Daer = 6 ! 3d ec3503Daer = 7 ! 3d abs3503Daer = 8 ! 3d ntracer_2d = 66 ps = 1 ! 2d precip = 2 ! 2d sconcoa = 3 ! 2d sconcbc = 4 ! 2d sconcso4 = 5 ! 2d sconcdust = 6 ! 2d sconcss = 7 ! 2d sconcno3 = 8 ! 2d sconcnh4 = 9 ! 2d sconcmsa = 10 ! 2d loadoa = 11 ! 2d loadbc = 12 ! 2d loadso4 = 13 ! 2d loaddust = 14 ! 2d loadss = 15 ! 2d loadno3 = 16 ! 2d emioa = 17 ! 2d emibc = 18 ! 2d emiso2 = 19 ! 2d emiso4 = 20 ! 2d emidust = 21 ! 2d emidms = 22 ! 2d emiss = 23 ! 2d dryso2 = 24 ! 2d dryoa = 25 ! 2d drybc = 26 ! 2d dryso4 = 27 ! 2d drydust = 28 ! 2d drydms = 29 ! 2d dryss = 30 ! 2d wetoa = 31 ! 2d wetbc = 32 ! 2d wetso2 = 33 ! 2d wetso4 = 34 ! 2d wetdust = 35 ! 2d wetdms = 36 ! 2d wetss = 37 ! 2d ! --- from here onwards keep consistent with order in optics (optics.F90) ! --- begin split order od550aer = 38 ! 2d od550so4 = 39 ! 2d od550bc = 40 ! 2d od550oa = 41 ! 2d od550ss = 42 ! 2d od550dust = 43 ! 2d od550no3 = 44 ! 2d od550aerh2o = 45 ! 2d od550lt1aer = 46 ! 2d od550lt1dust = 47 ! 2d od550lt1ss = 48 ! 2d ! --- end split order abs550aer = 49 ! 2d asyaer = 50 ! 2d ec550aer = 51 ! 2d ! --- begin in-situ data order ec550dryaer = 52 ! 2d abs550dryaer = 53 ! 2d asydryaer = 54 ! 2d ec550drylt1aer = 55 ! 2d abs550drylt1aer = 56 ! 2d ! --- end in-situ data order ! od440aer = 57 ! 2d abs440aer = 58 ! 2d ec440dryaer = 59 ! 2d abs440dryaer = 60 ! 2d ! od870aer = 61 ! 2d abs870aer = 62 ! 2d ec870dryaer = 63 ! 2d abs870dryaer = 64 ! 2d ! od350aer = 65 ! 2d abs350aer = 66 ! 2d if (stat_output) then ! ----------------------- ! set station information ! ----------------------- ! this here is a kind of redundant coding, but it enables to easily check for correct ! information (instead of assigning dedicated arrays (name,lat,lon etc.) side by side) ec550dryaer_stat=33 abs550dryaer_stat=34 asydryaer_stat=35 ec550drylt1aer_stat=36 abs550drylt1aer_stat=37 ec440dryaer_stat=38 abs440dryaer_stat=39 ec870dryaer_stat=40 abs870dryaer_stat=41 nextra_1dstat=3 ! temperature, specific humidity and pressure at model levels ntracer_1dstat=41+nextra_1dstat ntracer_0dstat=ntracer_1dstat+2 ! plus surface pressure and precipitation ! Adapted station list for AeroCom Phase III in-situ projects nstations = 90 allocate( statlist( nstations ) ) sc = 1 statlist(sc) = station( sc, 'Alert ', 82.499, -62.342, 210.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Anmyeon-do ', 36.540, 126.330, 94.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Annaberg-Buchholz ', 50.570, 13.000, 545.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Appalachian State U. ', 36.213, -81.692, 1076.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Aspvreten ', 58.806, 17.388, 20.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Auchencorth Moss ', 55.793, -3.245, 260.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Backgarden ', 23.540, 113.060, -999.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Barrow ', 71.323, -156.611, 11.0 ); sc = sc + 1 statlist(sc) = station( sc, 'BEO Moussala ', 42.179, 23.586, 2925.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Birkenes ', 58.388, 8.252, 190.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Bondville ', 40.050, -88.367, 213.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Bukit Kototabang ', -0.202, 100.318, 864.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Bosel ', 52.998, 7.943, 40.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Cabauw ', 51.971, 4.927, 1.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Cape Cod ', 40.070, -70.200, 43.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Cape Grim ', -40.682, 144.688, 94.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Cape Point ', -34.353, 18.490, 230.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Cape San Juan ', 18.381, -65.618, 65.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Chacaltaya ', -16.200, -68.100, 5320.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Danum Valley ', 4.981, 117.844, 426.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Demokritos ', 37.995, 23.816, 270.0 ); sc = sc + 1 statlist(sc) = station( sc, 'East Trout Lake ', 54.350, -104.983, 500.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Egbert ', 44.230, -79.783, 255.0 ); sc = sc + 1 statlist(sc) = station( sc, 'El Arenosillo ', 37.104, -6.734, 41.0 ); sc = sc + 1 statlist(sc) = station( sc, 'El Tololo ', -30.173, -70.799, 2220.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Elandsfontein ', -25.533, 25.750, 1424.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Finokalia ', 35.333, 25.670, 250.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Giordan Lighthouse ', 36.073, 14.219, 167.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Gosan ', 33.280, 126.170, 72.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Graciosa ', 39.080, -28.030, 30.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Granada ', 37.164, -3.605, 680.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Hesselbach ', 48.540, 8.400, 51.4 ); sc = sc + 1 statlist(sc) = station( sc, 'Harwell ', 51.573, -1.317, 137.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Hohenpeissenberg ', 47.802, 11.010, 985.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Hyytiala ', 61.847, 24.295, 181.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Illmitz ', 47.767, 16.767, 117.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Ispra ', 45.803, 8.627, 209.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Izana ', 28.309, -16.499, 2373.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Jungfraujoch ', 46.547, 7.985, 3580.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Kosetice ', 49.583, 15.083, 535.0 ); sc = sc + 1 statlist(sc) = station( sc, 'K-Puszta ', 46.967, 19.583, 125.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Leibzig ', 51.353, 12.435, 126.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Leibzig-West ', 51.318, 12.298, 122.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Lille Valby ', 55.683, 12.117, 20.0 ); sc = sc + 1 statlist(sc) = station( sc, 'London - N. Kensington', 51.521, -0.214, 27.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Lulin ', 23.470, 120.870, 2862.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Mace Head ', 53.326, -9.899, 5.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Madrid ', 40.450, -3.720, 669.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Manacapuro ', -3.210, -60.590, 50.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Manaus ', -2.595, -60.209, 45.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Mauna Loa ', 19.536, -155.576, 3397.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Melpitz ', 51.530, 12.930, 86.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Montseny ', 41.767, 2.350, 700.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Monte Cimone ', 44.167, 10.683, 2165.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Mount Waliguan ', 36.283, 100.900, 3810.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Nepal Climate Obs. ', 27.958, 86.815, 5079.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Nainital ', 29.360, 79.460, 1936.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Neuglobsov ', 53.143, 13.033, 62.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Neumayer ', -70.666, -8.266, 42.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Niamey ', 13.470, 2.170, 223.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Pallas ', 67.974, 24.116, 560.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Payerne ', 46.817, 6.950, 490.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Prague-Suchdol ', 50.130, 14.380, 270.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Preila ', 55.350, 21.067, 5.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Pt. Reyes ', 38.091, -122.957, 8.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Puy de Dome ', 45.772, 2.966, 1465.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Resolute Bay ', 74.717, -94.983, 64.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Sable Island ', 43.933, -60.017, 2.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Samoa ', -14.232, -170.563, 77.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Schauinsland ', 47.900, 7.917, 1205.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Shang Dianzi ', 40.390, 117.070, 294.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Shouxian ', 32.558, 116.782, 23.0 ); sc = sc + 1 statlist(sc) = station( sc, 'SIRTA ', 48.709, 2.159, 160.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Sonnblick ', 47.054, 12.959, 3106.0 ); sc = sc + 1 statlist(sc) = station( sc, 'South Pole ', -89.997, -24.800, 2841.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Southern Great Plains ', 36.617, -97.500, 318.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Storm Peak ', 40.455, -106.744, 3220.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Summit ', 72.580, -38.480, 3238.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Tahkuse ', 58.520, 24.940, 23.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Tiksi ', 71.586, 128.919, 8.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Trinidad Head ', 41.054, -124.151, 107.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Trollhaugen ', -72.012, 2.535, 1553.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Uto ', 59.783, 21.383, 7.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Varrio ', 67.767, 29.583, 400.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Vavihill ', 56.017, 13.150, 172.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Vielsalm ', 50.304, 6.001, 496.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Waldhof ', 52.802, 10.759, 74.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Whistler Mountain ', 50.059, -122.958, 2182.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Zeppelin ', 78.907, 11.889, 475.0 ); sc = sc + 1 statlist(sc) = station( sc, 'Zugspitze ', 47.417, 10.980, 2650.0 ); sc = sc + 1 ! this copy is needed for the netcdf output of the 2d-character-array [statnamelen,nstations] allocate( tmpstatnames( nstations ) ) tmpstatnames = statlist(:)%statname allocate(stat_active(nstations)) stat_active(:)=.false. stat_halo_loc=0 endif end if regionloop: do region = 1, nregions ! if region given, cycle if other region! if (present(iregion)) then if(iregion /= region) cycle regionloop endif imr = global_lli(region)%nlon jmr = global_lli(region)%nlat lmr = levi%nlev if( nbud_vg /= lmr ) then write(gol,*)'output_aerocom_init: nbud_vg /= lmr'; call goErr write(gol,*)'output_aerocom_init: YOU MUST ADD THE PROJ/BUDGET10 TO SRC CODE !!!!!'; call goErr TRACEBACK status=1; return end if ! initialize the output: if( newsrun ) then call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) ! target structures for output allocate(mixf (region)%f3d(ntracer_3d)) allocate(mixf (region)%f2d(ntracer_2d)) ! allocate fields do i=1,ntracer_3d allocate( mixf(region)%f3d(i)%field(i1:i2,j1:j2,lmr) ) end do do i=1,ntracer_2d allocate( mixf(region)%f2d(i)%field(i1:i2,j1:j2) ) end do if (stat_output) then ! --------- ! find coordinates in grid and get interpolation coefficients ! --------- allocate( sintp(region)%ixl(nstations), sintp(region)%wixl(nstations), & sintp(region)%ixr(nstations), sintp(region)%wixr(nstations), & sintp(region)%jyl(nstations), sintp(region)%wjyl(nstations), & sintp(region)%jyr(nstations), sintp(region)%wjyr(nstations), & sintp(region)%lzl(nstations), sintp(region)%wlzl(nstations), & sintp(region)%lzr(nstations), sintp(region)%wlzr(nstations), & sintp(region)%tght(nstations) ) call aerocom_get_statintp( region, status ) IF_NOTOK_RETURN(status=1) allocate(statf(region)%f1d(ntracer_1dstat)) allocate(statf(region)%f0d(ntracer_0dstat)) do i=1,ntracer_1dstat allocate( statf(region)%f1d(i)%field(nstations,lmr) ) end do do i=1,ntracer_0dstat allocate( statf(region)%f0d(i)%field(nstations) ) end do endif ! 3d data mixf(region)%f3d(temp )%mf = metafields( -1 , 'temp ', 'Temperature ' , 'K ', '', 'air_temperature' ) mixf(region)%f3d(hus )%mf = metafields( -1 , 'hus ', 'Specific Humidity ' , '1 ', '', 'specific_humidity' ) mixf(region)%f3d(airmass )%mf = metafields( -1 , 'airmass ', 'Air Mass ' , 'kg m-2 ', '', 'atmosphere_mass_of_air_per_unit_area' ) mixf(region)%f3d(pressure )%mf = metafields( -1 , 'pressure ', 'Pressure ' , 'Pa ', '', 'air_pressure' ) if (aai_output) then mixf(region)%f3d(ec5503Daer )%mf = metafields( -1 , 'ec5503Daer ', 'Aerosol Extinction @550nm ' , 'm-1 ', '', 'atmosphere_extinction_due_to_ambient_aerosol' ) mixf(region)%f3d(abs5503Daer )%mf = metafields( -1 , 'abs5503Daer ', 'Aerosol Absorption @550nm ' , 'm-1 ', '', 'atmosphere_absorption_due_to_ambient_aerosol' ) mixf(region)%f3d(ec3503Daer )%mf = metafields( -1 , 'ec3503Daer ', 'Aerosol Extinction @350nm ' , 'm-1 ', '', 'atmosphere_extinction_due_to_ambient_aerosol' ) mixf(region)%f3d(abs3503Daer )%mf = metafields( -1 , 'abs3503Daer ', 'Aerosol Absorption @350nm ' , 'm-1 ', '', 'atmosphere_absorption_due_to_ambient_aerosol' ) endif ! 2d data mixf(region)%f2d(ps )%mf = metafields( -1 , 'ps ', 'Surface Air Pressure ' , 'Pa ', '', 'surface_air_pressure' ) mixf(region)%f2d(precip )%mf = metafields( -1 , 'precip ', 'Precipitation ' , 'kg m-2 s-1', '', 'precipitation_flux' ) mixf(region)%f2d(sconcoa )%mf = metafields( -1 , 'sconcoa ', 'Surface Concentration POM ' , 'kg m-3 ', '', 'mass_concentration_of_particulate_organic_matter_dry_aerosol_in_air' ) mixf(region)%f2d(sconcbc )%mf = metafields( -1 , 'sconcbc ', 'Surface Concentration BC ' , 'kg m-3 ', '', 'mass_concentration_of_black_carbon_dry_aerosol_in_air' ) mixf(region)%f2d(sconcso4 )%mf = metafields( -1 , 'sconcso4 ', 'Surface Concentration SO4 ' , 'kg m-3 ', '', 'mass_concentration_of_sulfate_dry_aerosol_in_air' ) mixf(region)%f2d(sconcdust )%mf = metafields( -1 , 'sconcdust ', 'Surface Concentration Dust ' , 'kg m-3 ', '', 'mass_concentration_of_dust_dry_aerosol_in_air' ) mixf(region)%f2d(sconcss )%mf = metafields( -1 , 'sconcss ', 'Surface Concentration SS ' , 'kg m-3 ', '', 'mass_concentration_of_seasalt_dry_aerosol_in_air' ) mixf(region)%f2d(sconcno3 )%mf = metafields( -1 , 'sconcno3 ', 'Surface Concentration NO3 ' , 'kg m-3 ', '', 'mass_concentration_of_nitrate_dry_aerosol_in_air' ) mixf(region)%f2d(sconcnh4 )%mf = metafields( -1 , 'sconcnh4 ', 'Surface Concentration NH4 ' , 'kg m-3 ', '', 'mass_concentration_of_ammonium_dry_aerosol_in_air' ) mixf(region)%f2d(sconcmsa )%mf = metafields( -1 , 'sconcmsa ', 'Surface Concentration MSA ' , 'kg m-3 ', '', 'mass_concentration_of_methane_sulfonic_acid_in_air' ) mixf(region)%f2d(loadoa )%mf = metafields( -1 , 'loadoa ', 'Load of POM ' , 'kg m-2 ', '', 'atmosphere_mass_content_of_particulate_organic_matter_dry_aerosol' ) mixf(region)%f2d(loadbc )%mf = metafields( -1 , 'loadbc ', 'Load of BC ' , 'kg m-2 ', '', 'atmosphere_mass_content_of_black_carbon_dry_aerosol' ) mixf(region)%f2d(loadso4 )%mf = metafields( -1 , 'loadso4 ', 'Load of SO4 ' , 'kg m-2 ', '', 'atmosphere_mass_content_of_sulfate_dry_aerosol' ) mixf(region)%f2d(loaddust )%mf = metafields( -1 , 'loaddust ', 'Load of Dust ' , 'kg m-2 ', '', 'atmosphere_mass_content_of_dust_dry_aerosol' ) mixf(region)%f2d(loadss )%mf = metafields( -1 , 'loadss ', 'Load of SS ' , 'kg m-2 ', '', 'atmosphere_mass_content_of_seasalt_dry_aerosol' ) mixf(region)%f2d(loadno3 )%mf = metafields( -1 , 'loadno3 ', 'Load of NO3 ' , 'kg m-2 ', '', 'atmosphere_mass_content_of_nitrate_dry_aerosol' ) mixf(region)%f2d(emioa )%mf = metafields( -1 , 'emioa ', 'Total Emission of POM ' , 'kg m-2 s-1', 'up', 'tendency_of_atmosphere_mass_content_of_primary_particulate_organic_matter_dry_aerosol_due_to_emission' ) mixf(region)%f2d(emibc )%mf = metafields( -1 , 'emibc ', 'Total Emission of BC ' , 'kg m-2 s-1', 'up', 'tendency_of_atmosphere_mass_content_of_black_carbon_dry_aerosol_due_to_emission' ) mixf(region)%f2d(emiso2 )%mf = metafields( -1 , 'emiso2 ', 'Total Emission of SO2 ' , 'kg m-2 s-1', 'up', 'tendency_of_atmosphere_mass_content_of_sulfur_dioxide_due_to_emission' ) mixf(region)%f2d(emiso4 )%mf = metafields( -1 , 'emiso4 ', 'Total Direct Emission of SO4 ' , 'kg m-2 s-1', 'up', 'tendency_of_atmosphere_mass_content_of_sulfate_dry_aerosol_due_to_emission' ) mixf(region)%f2d(emidust )%mf = metafields( -1 , 'emidust ', 'Total Emission of Dust ' , 'kg m-2 s-1', 'up', 'tendency_of_atmosphere_mass_content_of_dust_dry_aerosol_due_to_emission' ) mixf(region)%f2d(emidms )%mf = metafields( -1 , 'emidms ', 'Total Emission of DMS ' , 'kg m-2 s-1', 'up', 'tendency_of_atmosphere_mass_content_of_dimethyl_sulfide_due_to_emission' ) mixf(region)%f2d(emiss )%mf = metafields( -1 , 'emiss ', 'Total Emission of SeaSalt ' , 'kg m-2 s-1', 'up', 'tendency_of_atmosphere_mass_content_of_seasalt_dry_aerosol_due_to_emission' ) mixf(region)%f2d(dryso2 )%mf = metafields( -1 , 'dryso2 ', 'Dry Deposition of SO2 ' , 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_sulfur_dioxide_due_to_dry_deposition' ) mixf(region)%f2d(dryoa )%mf = metafields( -1 , 'dryoa ', 'Dry Deposition of POM ' , 'kg m-2 s-1', 'down', & 'tendency_of_atmosphere_mass_content_of_primary_particulate_organic_matter_dry_aerosol_due_to_dry_deposition' ) mixf(region)%f2d(drybc )%mf = metafields( -1 , 'drybc ', 'Dry Deposition of BC ' , 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_black_carbon_dry_aerosol_due_to_dry_deposition' ) mixf(region)%f2d(dryso4 )%mf = metafields( -1 , 'dryso4 ', 'Dry Deposition of SO4 ' , 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_sulfate_dry_aerosol_due_to_dry_deposition' ) mixf(region)%f2d(drydust )%mf = metafields( -1 , 'drydust ', 'Dry Deposition of Dust ' , 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_dust_dry_aerosol_due_to_dry_deposition' ) mixf(region)%f2d(drydms )%mf = metafields( -1 , 'drydms ', 'Dry Deposition of DMS ' , 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_dimethyl_sulfide_due_to_dry_deposition' ) mixf(region)%f2d(dryss )%mf = metafields( -1 , 'dryss ', 'Dry Deposition of SeaSalt ' , 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_seasalt_dry_aerosol_due_to_dry_deposition' ) mixf(region)%f2d(wetoa )%mf = metafields( -1 , 'wetoa ', 'Wet Deposition of POM ' , 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_particulate_organic_matter_dry_aerosol_due_to_wet_deposition' ) mixf(region)%f2d(wetbc )%mf = metafields( -1 , 'wetbc ', 'Wet Deposition of BC ' , 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_black_carbon_dry_aerosol_due_to_wet_deposition' ) mixf(region)%f2d(wetso2 )%mf = metafields( -1 , 'wetso2 ', 'Wet Deposition of SO2 ' , 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_sulfur_dioxide_due_to_wet_deposition' ) mixf(region)%f2d(wetso4 )%mf = metafields( -1 , 'wetso4 ', 'Wet Deposition of SO4 ' , 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_sulfate_dry_aerosol_due_to_wet_deposition' ) mixf(region)%f2d(wetdust )%mf = metafields( -1 , 'wetdust ', 'Wet Deposition of Dust ' , 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_dust_dry_aerosol_due_to_wet_deposition' ) mixf(region)%f2d(wetdms )%mf = metafields( -1 , 'wetdms ', 'Wet Deposition of DMS ' , 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_dimethyl_sulfide_due_to_wet_deposition' ) mixf(region)%f2d(wetss )%mf = metafields( -1 , 'wetss ', 'Wet Deposition of Seasalt ' , 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_seasalt_dry_aerosol_due_to_wet_deposition' ) mixf(region)%f2d(od550aer )%mf = metafields( -1 , 'od550aer ', 'AOD@550nm ' , '1 ', '', 'atmosphere_optical_thickness_due_to_ambient_aerosol' ) mixf(region)%f2d(od550so4 )%mf = metafields( -1 , 'od550so4 ', 'SO4 AOD@550nm ' , '1 ', '', 'atmosphere_optical_thickness_due_to_sulfate_ambient_aerosol' ) mixf(region)%f2d(od550bc )%mf = metafields( -1 , 'od550bc ', 'Black Carbon AOD@550nm ' , '1 ', '', 'atmosphere_optical_thickness_due_to_black_carbon_ambient_aerosol' ) mixf(region)%f2d(od550oa )%mf = metafields( -1 , 'od550oa ', 'POM AOD@550nm ' , '1 ', '', 'atmosphere_optical_thickness_due_to_particulate_organic_matter_ambient_aerosol' ) mixf(region)%f2d(od550ss )%mf = metafields( -1 , 'od550ss ', 'SeaSalt AOD@550nm ' , '1 ', '', 'atmosphere_optical_thickness_due_to_seasalt_ambient_aerosol' ) mixf(region)%f2d(od550dust )%mf = metafields( -1 , 'od550dust ', 'Dust AOD@550nm ' , '1 ', '', 'atmosphere_optical_thickness_due_to_dust_ambient_aerosol' ) mixf(region)%f2d(od550no3 )%mf = metafields( -1 , 'od550no3 ', 'Nitrate AOD@550nm ' , '1 ', '', 'atmosphere_optical_thickness_due_to_nitrate_ambient_aerosol' ) mixf(region)%f2d(od550aerh2o )%mf = metafields( -1 , 'od550aerh2o ', 'Aerosol Water AOD@550nm ' , '1 ', '', 'atmosphere_optical_thickness_due_to_water_in_ambient_aerosol' ) mixf(region)%f2d(od550lt1aer )%mf = metafields( -1 , 'od550lt1aer ', 'Fine Mode AOD@550nm ' , '1 ', '', 'atmosphere_optical_thickness_due_to_pm1_ambient_aerosol' ) mixf(region)%f2d(od550lt1dust)%mf = metafields( -1 , 'od550lt1dust', 'Fine Mode Dust AOD@550nm ' , '1 ', '', 'atmosphere_optical_thickness_due_to_dust_pm1_ambient_aerosol' ) mixf(region)%f2d(od550lt1ss)%mf = metafields( -1 , 'od550lt1ss ', 'Fine Mode SeaSalt AOD@550nm ' , '1 ', '', 'atmosphere_optical_thickness_due_to_seasalt_pm1_ambient_aerosol' ) mixf(region)%f2d(abs550aer )%mf = metafields( -1 , 'abs550aer ', 'Absorption AOD@550nm ' , '1 ', '', 'atmosphere_absorption_optical_thickness_due_to_ambient_aerosol' ) mixf(region)%f2d(asyaer )%mf = metafields( -1 , 'asyaer ', 'Asymmetry Parameter ' , '1 ', '', 'atmosphere_asymmetry_parameter_ambient_aerosol' ) mixf(region)%f2d(ec550aer )%mf = metafields( -1 , 'ec550aer ', 'Surface Ambient Aerosol Extinction@550nm' , 'm-1 ', '', 'surface_extinction_due_to_ambient_aerosol' ) mixf(region)%f2d(ec550dryaer )%mf = metafields( -1 , 'ec550dryaer ', 'Surface Dry Aerosol Extinction@550 nm','m-1 ', '', 'surface_extinction_due_to_dry_aerosol' ) mixf(region)%f2d(abs550dryaer)%mf = metafields( -1 , 'abs550dryaer', 'Surface Dry Aerosol Absorption@550 nm','m-1 ', '', 'surface_absorption_due_to_dry_aerosol' ) mixf(region)%f2d(asydryaer )%mf = metafields( -1 , 'asydryaer ', 'Surface Dry Aersol Asymmetry Parameter' , '1 ', '', 'surface_asymmetry_parameter_dry_aerosol' ) mixf(region)%f2d(ec550drylt1aer )%mf = metafields( -1 , 'ec550drylt1aer ', 'Surface Fine Mode Dry Aerosol Extinction@550 nm' , 'm-1', '', 'surface_extinction_due_to_pm1_dry_aerosol' ) mixf(region)%f2d(abs550drylt1aer)%mf = metafields( -1 , 'abs550drylt1aer', 'Surface Fine Mode Dry Aerosol Absorption@550 nm' , 'm-1', '', 'surface_absorption_due_to_pm1_dry_aerosol' ) mixf(region)%f2d(od440aer )%mf = metafields( -1 , 'od440aer ', 'AOD@440nm ' , '1 ', '', 'atmosphere_optical_thickness_due_to_ambient_aerosol' ) mixf(region)%f2d(abs440aer )%mf = metafields( -1 , 'abs440aer ', 'Absorption AOD@440nm ' , '1 ', '', 'atmosphere_absorption_optical_thickness_due_to_ambient_aerosol' ) mixf(region)%f2d(ec440dryaer )%mf = metafields( -1 , 'ec440dryaer ', 'Surface Dry Aerosol Extinction@440nm' , 'm-1 ', '', 'surface_extinction_due_to_dry_aerosol' ) mixf(region)%f2d(abs440dryaer)%mf = metafields( -1 , 'abs440dryaer', 'Surface Dry Aerosol Absorption@440nm' , 'm-1 ', '', 'surface_absorption_due_to_dry_aerosol' ) mixf(region)%f2d(od870aer )%mf = metafields( -1 , 'od870aer ', 'AOD@870nm ' , '1 ', '', 'atmosphere_optical_thickness_due_to_ambient_aerosol' ) mixf(region)%f2d(abs870aer )%mf = metafields( -1 , 'abs870aer ', 'Absorption AOD@870nm ' , '1 ', '', 'atmosphere_absorption_optical_thickness_due_to_ambient_aerosol' ) mixf(region)%f2d(ec870dryaer )%mf = metafields( -1 , 'ec870dryaer ', 'Surface Dry Aerosol Extinction@870nm' , 'm-1 ', '', 'surface_extinction_due_to_dry_aerosol' ) mixf(region)%f2d(abs870dryaer)%mf = metafields( -1 , 'abs870dryaer', 'Surface Dry Aerosol Absorption@870nm' , 'm-1 ', '', 'surface_absorption_due_to_dry_aerosol' ) mixf(region)%f2d(od350aer )%mf = metafields( -1 , 'od350aer ', 'AOD@350nm ' , '1 ', '', 'atmosphere_optical_thickness_due_to_ambient_aerosol' ) mixf(region)%f2d(abs350aer )%mf = metafields( -1 , 'abs350aer ', 'Absorption AOD@350nm ' , '1 ', '', 'atmosphere_absorption_optical_thickness_due_to_ambient_aerosol' ) !!$ mixf(region)%f2d(52)%mf = metafields( 'od550soa ', 'AOD@550nm SOA ' , '1 ' ) !!$ mixf(region)%f2d(53)%mf = metafields( 'od550bb ', 'AOD@550nm BB ' , '1 ' ) !!$ mixf(region)%f2d(54)%mf = metafields( 'gf90aer ', 'GF @ 90 % RH ' , '1 ' ) if (stat_output) then ! 1d fields (stations) statf(region)%f1d( 1)%mf = metafields2( iso4nus, 'mmr1Dtr01', 'mmr of tracer 01' , 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air' ) statf(region)%f1d( 2)%mf = metafields2( iso4ais, 'mmr1Dtr02', 'mmr of tracer 02' , 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air' ) statf(region)%f1d( 3)%mf = metafields2( iso4acs, 'mmr1Dtr03', 'mmr of tracer 03' , 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air' ) statf(region)%f1d( 4)%mf = metafields2( iso4cos, 'mmr1Dtr04', 'mmr of tracer 04' , 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air' ) statf(region)%f1d( 5)%mf = metafields2( ibcais , 'mmr1Dtr05', 'mmr of tracer 05' , 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air' ) statf(region)%f1d( 6)%mf = metafields2( ibcacs , 'mmr1Dtr06', 'mmr of tracer 06' , 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air' ) statf(region)%f1d( 7)%mf = metafields2( ibccos , 'mmr1Dtr07', 'mmr of tracer 07' , 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air' ) statf(region)%f1d( 8)%mf = metafields2( ibcaii , 'mmr1Dtr08', 'mmr of tracer 08' , 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air' ) statf(region)%f1d( 9)%mf = metafields2( ipomais, 'mmr1Dtr09', 'mmr of tracer 09' , 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air' ) statf(region)%f1d(10)%mf = metafields2( ipomacs, 'mmr1Dtr10', 'mmr of tracer 10' , 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air' ) statf(region)%f1d(11)%mf = metafields2( ipomcos, 'mmr1Dtr11', 'mmr of tracer 11' , 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air' ) statf(region)%f1d(12)%mf = metafields2( ipomaii, 'mmr1Dtr12', 'mmr of tracer 12' , 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air' ) statf(region)%f1d(13)%mf = metafields2( issacs , 'mmr1Dtr13', 'mmr of tracer 13' , 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air' ) statf(region)%f1d(14)%mf = metafields2( isscos , 'mmr1Dtr14', 'mmr of tracer 14' , 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air' ) statf(region)%f1d(15)%mf = metafields2( iduacs , 'mmr1Dtr15', 'mmr of tracer 15' , 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air' ) statf(region)%f1d(16)%mf = metafields2( iducos , 'mmr1Dtr16', 'mmr of tracer 16' , 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air' ) statf(region)%f1d(17)%mf = metafields2( iduaci , 'mmr1Dtr17', 'mmr of tracer 17' , 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air' ) statf(region)%f1d(18)%mf = metafields2( iducoi , 'mmr1Dtr18', 'mmr of tracer 18' , 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air' ) statf(region)%f1d(19)%mf = metafields2( inus_n , 'conccn1Dmode01', 'number concentration of mode 01' , 'm-3' , 'number_concentration_of_ambient_aerosol_in_air' ) statf(region)%f1d(20)%mf = metafields2( iais_n , 'conccn1Dmode02', 'number concentration of mode 02' , 'm-3' , 'number_concentration_of_ambient_aerosol_in_air' ) statf(region)%f1d(21)%mf = metafields2( iacs_n , 'conccn1Dmode03', 'number concentration of mode 03' , 'm-3' , 'number_concentration_of_ambient_aerosol_in_air' ) statf(region)%f1d(22)%mf = metafields2( icos_n , 'conccn1Dmode04', 'number concentration of mode 04' , 'm-3' , 'number_concentration_of_ambient_aerosol_in_air' ) statf(region)%f1d(23)%mf = metafields2( iaii_n , 'conccn1Dmode05', 'number concentration of mode 05' , 'm-3' , 'number_concentration_of_ambient_aerosol_in_air' ) statf(region)%f1d(24)%mf = metafields2( iaci_n , 'conccn1Dmode06', 'number concentration of mode 06' , 'm-3' , 'number_concentration_of_ambient_aerosol_in_air' ) statf(region)%f1d(25)%mf = metafields2( icoi_n , 'conccn1Dmode07', 'number concentration of mode 07' , 'm-3' , 'number_concentration_of_ambient_aerosol_in_air' ) statf(region)%f1d(26)%mf = metafields2( 1 , 'mmr1Daerh2o01 ', 'mmr of aerosol water of mode 01' , 'kg kg-1', 'mass_fraction_of_water_in_ambient_aerosol_in_air' ) statf(region)%f1d(27)%mf = metafields2( 2 , 'mmr1Daerh2o02 ', 'mmr of aerosol water of mode 02' , 'kg kg-1', 'mass_fraction_of_water_in_ambient_aerosol_in_air' ) statf(region)%f1d(28)%mf = metafields2( 3 , 'mmr1Daerh2o03 ', 'mmr of aerosol water of mode 03' , 'kg kg-1', 'mass_fraction_of_water_in_ambient_aerosol_in_air' ) statf(region)%f1d(29)%mf = metafields2( 4 , 'mmr1Daerh2o04 ', 'mmr of aerosol water of mode 04' , 'kg kg-1', 'mass_fraction_of_water_in_ambient_aerosol_in_air' ) statf(region)%f1d(30)%mf = metafields2( ino3_a , 'mmr1Dno3_a ', 'mmr of nitrate aerosol' , 'kg kg-1', 'mass_fraction_of_nitrate_dry_aerosol_in_air' ) statf(region)%f1d(31)%mf = metafields2( inh4 , 'mmr1Dnh4 ', 'mmr of ammonium' , 'kg kg-1', 'mass_fraction_of_ammonium_dry_aerosol_in_air' ) statf(region)%f1d(32)%mf = metafields2( imsa , 'mmr1Dmsa ', 'mmr of methane sulfonic acid' , 'kg kg-1', 'mass_fraction_of_methane_sulfonic_acid_in_air' ) statf(region)%f1d(33)%mf = metafields2( -1 , 'ec550dry1Daer ', 'dry aerosol extinction at 550 nm' , 'm-1' , 'atmosphere_extinction_due_to_dry_aerosol' ) statf(region)%f1d(34)%mf = metafields2( -1 , 'abs550dry1Daer', 'dry aerosol absorption at 550 nm' , 'm-1' , 'atmosphere_absorption_due_to_dry_aerosol' ) statf(region)%f1d(35)%mf = metafields2( -1 , 'asydry1Daer', 'dry aerosol asymmetry parameter' , '1' , 'atmosphere_absorption_due_to_dry_aerosol' ) statf(region)%f1d(36)%mf = metafields2( -1 , 'ec550dry1Dlt1aer', 'fine mode dry aerosol extinction at 550 nm' , 'm-1' , 'atmosphere_absorption_due_to_dry_aerosol' ) statf(region)%f1d(37)%mf = metafields2( -1 , 'abs550dry1Dlt1aer', 'fine mode dry aerosol absorption at 550 nm' , 'm-1' , 'atmosphere_absorption_due_to_dry_aerosol' ) statf(region)%f1d(38)%mf = metafields2( -1 , 'ec440dry1Daer', 'dry aerosol extinction at 440 nm' , 'm-1' , 'atmosphere_absorption_due_to_dry_aerosol' ) statf(region)%f1d(39)%mf = metafields2( -1 , 'abs440dry1Daer', 'dry aerosol absorption at 440 nm' , 'm-1' , 'atmosphere_absorption_due_to_dry_aerosol' ) statf(region)%f1d(40)%mf = metafields2( -1 , 'ec870dry1Daer', 'dry aerosol extinction at 870 nm' , 'm-1' , 'atmosphere_absorption_due_to_dry_aerosol' ) statf(region)%f1d(41)%mf = metafields2( -1 , 'abs870dry1Daer', 'dry aerosol absorption at 870 nm' , 'm-1' , 'atmosphere_absorption_due_to_dry_aerosol' ) statf(region)%f1d(42)%mf = metafields2( -1 , 'temp1D ', 'Temperature ' , 'K' , 'air_temperature' ) statf(region)%f1d(43)%mf = metafields2( -1 , 'hus1D ', 'Specific Humidity ' , '1' , 'specific_humidity' ) statf(region)%f1d(44)%mf = metafields2( -1 , 'pressure1D ', 'Pressure ' , 'Pa' , 'air_pressure' ) ! 0d fields (stations) statf(region)%f0d( 1)%mf = metafields2( iso4nus, 'mmrtr01', 'mmr of tracer 01' , 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air' ) statf(region)%f0d( 2)%mf = metafields2( iso4ais, 'mmrtr02', 'mmr of tracer 02' , 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air' ) statf(region)%f0d( 3)%mf = metafields2( iso4acs, 'mmrtr03', 'mmr of tracer 03' , 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air' ) statf(region)%f0d( 4)%mf = metafields2( iso4cos, 'mmrtr04', 'mmr of tracer 04' , 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air' ) statf(region)%f0d( 5)%mf = metafields2( ibcais , 'mmrtr05', 'mmr of tracer 05' , 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air' ) statf(region)%f0d( 6)%mf = metafields2( ibcacs , 'mmrtr06', 'mmr of tracer 06' , 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air' ) statf(region)%f0d( 7)%mf = metafields2( ibccos , 'mmrtr07', 'mmr of tracer 07' , 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air' ) statf(region)%f0d( 8)%mf = metafields2( ibcaii , 'mmrtr08', 'mmr of tracer 08' , 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air' ) statf(region)%f0d( 9)%mf = metafields2( ipomais, 'mmrtr09', 'mmr of tracer 09' , 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air' ) statf(region)%f0d(10)%mf = metafields2( ipomacs, 'mmrtr10', 'mmr of tracer 10' , 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air' ) statf(region)%f0d(11)%mf = metafields2( ipomcos, 'mmrtr11', 'mmr of tracer 11' , 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air' ) statf(region)%f0d(12)%mf = metafields2( ipomaii, 'mmrtr12', 'mmr of tracer 12' , 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air' ) statf(region)%f0d(13)%mf = metafields2( issacs , 'mmrtr13', 'mmr of tracer 13' , 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air' ) statf(region)%f0d(14)%mf = metafields2( isscos , 'mmrtr14', 'mmr of tracer 14' , 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air' ) statf(region)%f0d(15)%mf = metafields2( iduacs , 'mmrtr15', 'mmr of tracer 15' , 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air' ) statf(region)%f0d(16)%mf = metafields2( iducos , 'mmrtr16', 'mmr of tracer 16' , 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air' ) statf(region)%f0d(17)%mf = metafields2( iduaci , 'mmrtr17', 'mmr of tracer 17' , 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air' ) statf(region)%f0d(18)%mf = metafields2( iducoi , 'mmrtr18', 'mmr of tracer 18' , 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air' ) statf(region)%f0d(19)%mf = metafields2( inus_n , 'conccnmode01', 'number concentration of mode 01' , 'm-3' , 'number_concentration_of_ambient_aerosol_in_air' ) statf(region)%f0d(20)%mf = metafields2( iais_n , 'conccnmode02', 'number concentration of mode 02' , 'm-3' , 'number_concentration_of_ambient_aerosol_in_air' ) statf(region)%f0d(21)%mf = metafields2( iacs_n , 'conccnmode03', 'number concentration of mode 03' , 'm-3' , 'number_concentration_of_ambient_aerosol_in_air' ) statf(region)%f0d(22)%mf = metafields2( icos_n , 'conccnmode04', 'number concentration of mode 04' , 'm-3' , 'number_concentration_of_ambient_aerosol_in_air' ) statf(region)%f0d(23)%mf = metafields2( iaii_n , 'conccnmode05', 'number concentration of mode 05' , 'm-3' , 'number_concentration_of_ambient_aerosol_in_air' ) statf(region)%f0d(24)%mf = metafields2( iaci_n , 'conccnmode06', 'number concentration of mode 06' , 'm-3' , 'number_concentration_of_ambient_aerosol_in_air' ) statf(region)%f0d(25)%mf = metafields2( icoi_n , 'conccnmode07', 'number concentration of mode 07' , 'm-3' , 'number_concentration_of_ambient_aerosol_in_air' ) statf(region)%f0d(26)%mf = metafields2( 1 , 'mmraerh2o01 ', 'mmr of aerosol water of mode 01' , 'kg kg-1', 'mass_fraction_of_water_in_ambient_aerosol_in_air' ) statf(region)%f0d(27)%mf = metafields2( 2 , 'mmraerh2o02 ', 'mmr of aerosol water of mode 02' , 'kg kg-1', 'mass_fraction_of_water_in_ambient_aerosol_in_air' ) statf(region)%f0d(28)%mf = metafields2( 3 , 'mmraerh2o03 ', 'mmr of aerosol water of mode 03' , 'kg kg-1', 'mass_fraction_of_water_in_ambient_aerosol_in_air' ) statf(region)%f0d(29)%mf = metafields2( 4 , 'mmraerh2o04 ', 'mmr of aerosol water of mode 04' , 'kg kg-1', 'mass_fraction_of_water_in_ambient_aerosol_in_air' ) statf(region)%f0d(30)%mf = metafields2( ino3_a , 'mmrno3_a ', 'mmr of nitrate aerosol' , 'kg kg-1', 'mass_fraction_of_nitrate_dry_aerosol_in_air' ) statf(region)%f0d(31)%mf = metafields2( inh4 , 'mmrnh4 ', 'mmr of ammonium' , 'kg kg-1', 'mass_fraction_of_ammonium_dry_aerosol_in_air' ) statf(region)%f0d(32)%mf = metafields2( imsa , 'mmrmsa ', 'mmr of methane sulfonic acid' , 'kg kg-1', 'mass_fraction_of_methane_sulfonic_acid_in_air' ) statf(region)%f0d(33)%mf = metafields2( -1 , 'ec550dryaer ', 'dry aerosol extinction at 550 nm' , 'm-1' , 'atmosphere_extinction_due_to_dry_aerosol' ) statf(region)%f0d(34)%mf = metafields2( -1 , 'abs550dryaer', 'dry aerosol absorption at 550 nm' , 'm-1' , 'atmosphere_absorption_due_to_dry_aerosol' ) statf(region)%f0d(35)%mf = metafields2( -1 , 'asydryaer', 'dry aerosol asymmetry parameter' , '1' , 'atmosphere_absorption_due_to_dry_aerosol' ) statf(region)%f0d(36)%mf = metafields2( -1 , 'ec550drylt1aer', 'fine mode dry aerosol extinction at 550 nm' , 'm-1' , 'atmosphere_absorption_due_to_dry_aerosol' ) statf(region)%f0d(37)%mf = metafields2( -1 , 'abs550drylt1aer', 'fine mode dry aerosol absorption at 550 nm' , 'm-1' , 'atmosphere_absorption_due_to_dry_aerosol' ) statf(region)%f0d(38)%mf = metafields2( -1 , 'ec440dryaer', 'dry aerosol extinction at 440 nm' , 'm-1' , 'atmosphere_absorption_due_to_dry_aerosol' ) statf(region)%f0d(39)%mf = metafields2( -1 , 'abs440dryaer', 'dry aerosol absorption at 440 nm' , 'm-1' , 'atmosphere_absorption_due_to_dry_aerosol' ) statf(region)%f0d(40)%mf = metafields2( -1 , 'ec870dryaer', 'dry aerosol extinction at 870 nm' , 'm-1' , 'atmosphere_absorption_due_to_dry_aerosol' ) statf(region)%f0d(41)%mf = metafields2( -1 , 'abs870dryaer', 'dry aerosol absorption at 870 nm' , 'm-1' , 'atmosphere_absorption_due_to_dry_aerosol' ) statf(region)%f0d(42)%mf = metafields2( -1 , 'temp ', 'Temperature ' , 'K' , 'air_temperature' ) statf(region)%f0d(43)%mf = metafields2( -1 , 'hus ', 'Specific Humidity ' , '1' , 'specific_humidity' ) statf(region)%f0d(44)%mf = metafields2( -1 , 'pressure ', 'Pressure ' , 'Pa' , 'air_pressure' ) statf(region)%f0d(45)%mf = metafields2( -1 , 'ps ', 'Surface Air Pressure ' , 'Pa ' , 'surface_air_pressure' ) statf(region)%f0d(46)%mf = metafields2( -1 , 'precip ', 'Precipitation ' , 'kg m-2 s-1', 'precipitation_flux' ) endif ! set global dimensions of fields for netcdf definitions mixf (region)%nlon = imr mixf (region)%nlat = jmr mixf (region)%nlev = lmr if (stat_output) then statf(region)%nstat = nstations statf(region)%nlev = lmr endif ! intermediate structures for budgets allocate(drydepos(region)%f2dslast(i1:i2,j1:j2,7)) allocate(wetdepos(region)%f2dslast(i1:i2,j1:j2,7)) allocate(emission(region)%f2dslast(i1:i2,j1:j2,7)) ! these here are the initial budgets (monthly): 0.0 drydepos(region)%f2dslast = 0.0 wetdepos(region)%f2dslast = 0.0 emission(region)%f2dslast = 0.0 endif ! newsrun ! reset counters and accumulators mixf (region)%acct = 0 do i = 1, ntracer_3d mixf(region)%f3d(i)%field = 0.0 end do do i = 1, ntracer_2d mixf(region)%f2d(i)%field = 0.0 end do if (stat_output) then ! stations statf(region)%acct = 0 do i = 1, ntracer_1dstat statf(region)%f1d(i)%field = 0.0 end do do i = 1, ntracer_0dstat statf(region)%f0d(i)%field = 0.0 end do endif ! ---------------- ! open NetCDF file (mixf) ! ---------------- call tau2date(itau,idater) ! time (in days since 2001-01-01 00:00) call date2tau( time_reftime6, itauref ) select case (trim(aerocom_freq)) case ('hourly') write(idates, '(i4,3i2.2)') idater(1), idater(2), idater(3), idater(4)+1 write(lidates, '(i4,"-",i2.2,"-",i2.2," ",i2.2,":30")') idater(1), idater(2), idater(3), idater(4) idater(5:6) = (/30,0/) case ('daily') write(idates, '(i4,2i2.2)') idater(1), idater(2), idater(3) write(lidates, '(i4,"-",i2.2,"-",i2.2," 12:00")') idater(1), idater(2), idater(3) ! set noon (recommendations) idater(4:6) = (/12,0,0/) case ('monthly') ! for monthly files, set time to middle of the month write(idates, '(i4,i2.2)') idater(1), idater(2) mlength=mlen(idater(2)) if ( mod(mlength,2) == 0 ) then idater(3:6) = (/mlength/2 + 1,0,0,0/) write(lidates, '(i4,"-",i2.2,"-",i2.2," 00:00")') idater(1), idater(2), idater(3) else idater(3:6) = (/(mlength+1)/2,12,0,0/) write(lidates, '(i4,"-",i2.2,"-",i2.2," 12:00")') idater(1), idater(2), idater(3) endif case default write (gol,'("Unknown AeroCom output frequency;")'); call goErr end select call date2tau( idater, itaucur ) reftime = (itaucur - itauref) / 86400. ! Changed file name convention to AeroCom Control 2015 mixf(region)%fname = trim(outdir)//'/'//& trim(f_dataid)//'_'//& trim(f_model) //'_'//& trim(aerocom_exper)//'_'//& trim(f_dimmix)//'_'//& trim(idates) //'_'//& trim(aerocom_freq) //'.nc' #ifdef MPI ! overwrite existing files (clobber), provide MPI stuff: call MDF_Create( mixf(region)%fname, MDF_NETCDF4, MDF_REPLACE, mixf(region)%funit, status, & mpi_comm=localComm, mpi_info=MPI_INFO_NULL ) if (status/=0) then write (gol,'("from creating NetCDF4 file for writing in parallel;")'); call goErr write (gol,'("MDF module not compiled with netcdf4_par support ?")'); call goErr TRACEBACK; status=1; return end if #else ! overwrite existing files (clobber) call MDF_Create( mixf(region)%fname, MDF_NETCDF4, MDF_REPLACE, mixf(region)%funit, status ) IF_NOTOK_RETURN(status=1) #endif if(okdebug) then write(gol,*) 'output_aerocom_init: File ', trim(mixf(region)%fname), ' opened ' ; call goPr endif ! global attributes call MDF_Put_Att( mixf(region)%funit, MDF_GLOBAL, 'title', 'Model output for AeroCom', status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit, MDF_GLOBAL, 'source', 'TM5-mp', status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit, MDF_GLOBAL, 'institution', 'Royal Netherlands Meteorological Institute (KNMI), De Bilt, The Netherlands)', status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit, MDF_GLOBAL, 'contact' , 'Twan van Noije; noije@knmi.nl', status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit, MDF_GLOBAL, 'project_id', 'AeroCom Phase 3', status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit, MDF_GLOBAL, 'conventions', 'CF-1.0 - HTAP', status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit, MDF_GLOBAL, 'date', lidates, status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit, MDF_GLOBAL, 'time_units', time_units, status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit, MDF_GLOBAL, 'references', 'http://tm5.sourceforge.net/', status ) IF_NOTOK_MDF(fid=mixf(region)%funit) ! define dimensions call MDF_Def_Dim( mixf(region)%funit, 'longitude', imr, mixf(region)%nlon, status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Def_Dim( mixf(region)%funit, 'latitude', jmr, mixf(region)%nlat, status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Def_Dim( mixf(region)%funit, 'alevel', lmr, mixf(region)%nlev, status ) IF_NOTOK_MDF(fid=mixf(region)%funit) ! define variables do i = 1, ntracer_3d call MDF_Def_Var( mixf(region)%funit, trim(mixf(region)%f3d(i)%mf%vname), MDF_FLOAT, & (/mixf(region)%nlon, mixf(region)%nlat, mixf(region)%nlev/), varid, status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Var_Par_Access( mixf(region)%funit, varid, access_mode, status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit, varid, 'long_name', trim(mixf(region)%f3d(i)%mf%lname), status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit, varid, 'standard_name', trim(mixf(region)%f3d(i)%mf%standard_name), status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit, varid, 'units', trim(mixf(region)%f3d(i)%mf%unit), status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit, varid, 'time', reftime, status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit, varid, 'cell_methods', 'time: mean', status ) IF_NOTOK_MDF(fid=mixf(region)%funit) if( mixf(region)%f3d(i)%mf%positive /= '' ) then call MDF_Put_Att( mixf(region)%funit, varid, 'positive', trim(mixf(region)%f3d(i)%mf%positive), status ) IF_NOTOK_MDF(fid=mixf(region)%funit) end if mixf(region)%f3d(i)%varid = varid end do do i = 1, ntracer_2d call MDF_Def_Var( mixf(region)%funit, trim(mixf(region)%f2d(i)%mf%vname), MDF_FLOAT, & (/mixf(region)%nlon, mixf(region)%nlat/), varid, status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Var_Par_Access( mixf(region)%funit, varid, access_mode, status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit, varid, 'long_name', trim(mixf(region)%f2d(i)%mf%lname), status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit, varid, 'standard_name', trim(mixf(region)%f2d(i)%mf%standard_name), status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit, varid, 'units', trim(mixf(region)%f2d(i)%mf%unit), status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit, varid, 'time', reftime, status ) IF_NOTOK_MDF(fid=mixf(region)%funit) call MDF_Put_Att( mixf(region)%funit, varid, 'cell_methods', 'time: mean', status ) IF_NOTOK_MDF(fid=mixf(region)%funit) if( mixf(region)%f2d(i)%mf%positive /= '' ) then call MDF_Put_Att( mixf(region)%funit, varid, 'positive', trim(mixf(region)%f2d(i)%mf%positive), status ) IF_NOTOK_MDF(fid=mixf(region)%funit) end if mixf(region)%f2d(i)%varid = varid end do call MDF_EndDef( mixf(region)%funit, status ) IF_NOTOK_MDF(fid=mixf(region)%funit) if (stat_output) then ! ------------ ! station NetCDF file ! ------------ statf(region)%fname = trim(outdir)//'/'//& trim(f_dataid) //'_'//& trim(f_model) //'_'//& trim(aerocom_exper) //'_'//& trim(f_dimstat)//'_'//& trim(idates) //'_'//& trim(aerocom_freq) //'.nc' #ifdef MPI ! overwrite existing files (clobber), provide MPI stuff: call MDF_Create( statf(region)%fname, MDF_NETCDF4, MDF_REPLACE, statf(region)%funit, status, & mpi_comm=localComm, mpi_info=MPI_INFO_NULL ) if (status/=0) then write (gol,'("from creating NetCDF4 file for writing in parallel;")'); call goErr write (gol,'("MDF module not compiled with netcdf4_par support ?")'); call goErr TRACEBACK; status=1; return end if #else ! overwrite existing files (clobber) call MDF_Create( statf(region)%fname, MDF_NETCDF4, MDF_REPLACE, statf(region)%funit, status ) IF_NOTOK_RETURN(status=1) #endif if(okdebug) then write(gol,*) ' output_aerocom_init: File ', trim(statf(region)%fname), ' opened '; call goPr endif ! global attributes call MDF_Put_Att( statf(region)%funit, MDF_GLOBAL, 'title', 'Model output for AeroCom', status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Put_Att( statf(region)%funit, MDF_GLOBAL, 'source', 'TM5-mp', status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Put_Att( statf(region)%funit, MDF_GLOBAL, 'institution', 'Royal Netherlands Meteorological Institute (KNMI), De Bilt, The Netherlands)', status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Put_Att( statf(region)%funit, MDF_GLOBAL, 'contact', 'Twan van Noije; noije@knmi.nl', status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Put_Att( statf(region)%funit, MDF_GLOBAL, 'project_id', 'AeroCom Phase 3', status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Put_Att( statf(region)%funit, MDF_GLOBAL, 'conventions', 'CF-1.0 - HTAP', status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Put_Att( statf(region)%funit, MDF_GLOBAL, 'date', lidates, status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Put_Att( statf(region)%funit, MDF_GLOBAL, 'time_units', time_units, status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Put_Att( statf(region)%funit, MDF_GLOBAL, 'references', 'http://tm5.sourceforge.net/', status ) IF_NOTOK_MDF(fid=statf(region)%funit) ! define dimensions call MDF_Def_Dim( statf(region)%funit, 'station', nstations, statf(region)%nstat, status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Def_Dim( statf(region)%funit, 'stationname_len', statnamelen, snamelendim, status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Def_Dim( statf(region)%funit, 'alevel', lmr, statf(region)%nlev, status ) IF_NOTOK_MDF(fid=statf(region)%funit) ! define variables do i = 1, ntracer_1dstat call MDF_Def_Var( statf(region)%funit, trim(statf(region)%f1d(i)%mf%vname), MDF_FLOAT, & (/statf(region)%nstat,statf(region)%nlev/), varid, status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Var_Par_Access( statf(region)%funit, varid, access_mode_sta, status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Put_Att( statf(region)%funit, varid, 'long_name', trim(statf(region)%f1d(i)%mf%lname), status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Put_Att( statf(region)%funit, varid, 'standard_name', trim(statf(region)%f1d(i)%mf%standard_name), status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Put_Att( statf(region)%funit, varid, 'units', trim(statf(region)%f1d(i)%mf%unit), status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Put_Att( statf(region)%funit, varid, 'time', reftime, status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Put_Att( statf(region)%funit, varid, 'cell_methods', 'time: mean', status ) IF_NOTOK_MDF(fid=statf(region)%funit) statf(region)%f1d(i)%varid = varid end do do i = 1, ntracer_0dstat call MDF_Def_Var( statf(region)%funit, trim(statf(region)%f0d(i)%mf%vname), MDF_FLOAT, & (/statf(region)%nstat/), varid, status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Var_Par_Access( statf(region)%funit, varid, access_mode_sta, status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Put_Att( statf(region)%funit, varid, 'long_name', trim(statf(region)%f0d(i)%mf%lname), status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Put_Att( statf(region)%funit, varid, 'standard_name', trim(statf(region)%f0d(i)%mf%standard_name), status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Put_Att( statf(region)%funit, varid, 'units', trim(statf(region)%f0d(i)%mf%unit), status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Put_Att( statf(region)%funit, varid, 'time', reftime, status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Put_Att( statf(region)%funit, varid, 'cell_methods', 'time: mean', status ) IF_NOTOK_MDF(fid=statf(region)%funit) statf(region)%f0d(i)%varid = varid end do ! station information (name, nb, lat, lon, alt) call MDF_Def_Var( statf(region)%funit, 'stationname', MDF_CHAR, (/snamelendim,statf(region)%nstat/), varid, status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Var_Par_Access( statf(region)%funit, varid, access_mode_sta, status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Put_Att( statf(region)%funit, varid, 'long_name', 'station name', status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Put_Att( statf(region)%funit, varid, 'standard_name', 'station_name', status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Put_Att( statf(region)%funit, varid, 'units', '1', status ) IF_NOTOK_MDF(fid=statf(region)%funit) statvarid(1) = varid call MDF_Def_Var( statf(region)%funit, 'stationnb', MDF_INT, (/statf(region)%nstat/), varid, status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Var_Par_Access( statf(region)%funit, varid, access_mode_sta, status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Put_Att( statf(region)%funit, varid, 'long_name', 'station number', status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Put_Att( statf(region)%funit, varid, 'standard_name', 'station_nb', status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Put_Att( statf(region)%funit, varid, 'units', '1', status ) IF_NOTOK_MDF(fid=statf(region)%funit) statvarid(2) = varid call MDF_Def_Var( statf(region)%funit, 'stationlat', MDF_FLOAT, (/statf(region)%nstat/), varid, status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Var_Par_Access( statf(region)%funit, varid, access_mode_sta, status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Put_Att( statf(region)%funit, varid, 'long_name', 'location of station: station latitude', status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Put_Att( statf(region)%funit, varid, 'standard_name', 'stationlat', status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Put_Att( statf(region)%funit, varid, 'units', '1', status ) IF_NOTOK_MDF(fid=statf(region)%funit) statvarid(3) = varid call MDF_Def_Var( statf(region)%funit, 'stationlon', MDF_FLOAT, (/statf(region)%nstat/), varid, status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Var_Par_Access( statf(region)%funit, varid, access_mode_sta, status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Put_Att( statf(region)%funit, varid, 'long_name', 'location of station: station longitude', status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Put_Att( statf(region)%funit, varid, 'standard_name', 'stationlon', status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Put_Att( statf(region)%funit, varid, 'units', '1', status ) IF_NOTOK_MDF(fid=statf(region)%funit) statvarid(4) = varid call MDF_Def_Var( statf(region)%funit, 'stationalt', MDF_FLOAT, (/statf(region)%nstat/), varid, status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Var_Par_Access( statf(region)%funit, varid, access_mode_sta, status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Put_Att( statf(region)%funit, varid, 'long_name', 'location of station: station altitude', status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Put_Att( statf(region)%funit, varid, 'standard_name', 'stationalt', status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Put_Att( statf(region)%funit, varid, 'units', '1', status ) IF_NOTOK_MDF(fid=statf(region)%funit) statvarid(5) = varid call MDF_EndDef( statf(region)%funit, status ) IF_NOTOK_MDF(fid=statf(region)%funit) ! put station information (note: behavior here is like if all procs write same data) call MDF_Put_Var( statf(region)%funit, statvarid(1), tmpstatnames, status, start=(/ 1, 1/), count=(/ statnamelen,nstations /) ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Put_Var( statf(region)%funit, statvarid(2), statlist(:)%statnb, status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Put_Var( statf(region)%funit, statvarid(3), statlist(:)%statlat, status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Put_Var( statf(region)%funit, statvarid(4), statlist(:)%statlon, status ) IF_NOTOK_MDF(fid=statf(region)%funit) call MDF_Put_Var( statf(region)%funit, statvarid(5), statlist(:)%statalt, status ) IF_NOTOK_MDF(fid=statf(region)%funit) endif end do regionloop ! nregions call goLabel() ; status = 0 end subroutine output_aerocom_init !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: OUTPUT_AEROCOM_DONE ! ! !DESCRIPTION: Free parameters. !\\ !\\ ! !INTERFACE: ! subroutine output_aerocom_done(stat_output, status, iregion) ! ! !USES: ! ! !INPUT PARAMETERS: ! logical, intent(in) :: stat_output ! true if stations output is requested integer, intent(in), optional :: iregion ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 29 Nov 2010 - Achim Strunk - ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/output_aerocom_done' integer :: i, region ! --- begin ------------------------------------- call goLabel(rname) deallocate( wdep_out ) regionloop: do region = 1, nregions ! if region given, cycle if other region! if (present(iregion)) then if(iregion /= region) cycle regionloop endif do i=1,ntracer_3d deallocate( mixf(region)%f3d(i)%field ) end do do i=1,ntracer_2d deallocate( mixf(region)%f2d(i)%field ) end do deallocate( mixf (region)%f3d ) deallocate( mixf (region)%f2d ) deallocate( drydepos(region)%f2dslast ) deallocate( wetdepos(region)%f2dslast ) deallocate( emission(region)%f2dslast ) if (stat_output) then ! stations do i=1,ntracer_1dstat deallocate( statf(region)%f1d(i)%field ) end do do i=1,ntracer_0dstat deallocate( statf(region)%f0d(i)%field ) end do deallocate( statf(region)%f1d ) deallocate( statf(region)%f0d ) deallocate( sintp(region)%ixl, sintp(region)%wixl, sintp(region)%ixr, sintp(region)%wixr, & sintp(region)%jyl, sintp(region)%wjyl, sintp(region)%jyr, sintp(region)%wjyr, & sintp(region)%lzl, sintp(region)%wlzl, sintp(region)%lzr, sintp(region)%wlzr, & sintp(region)%tght ) endif end do regionloop if (stat_output) then deallocate( statlist ) deallocate( tmpstatnames ) deallocate( stat_active ) endif call goLabel() ; status = 0 end subroutine output_aerocom_done !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: output_aerocom_write ! ! !DESCRIPTION: This routines builds the average by dividing accumulated ! data by stack counter. The results are written to file. !\\ !\\ ! !INTERFACE: ! subroutine output_aerocom_write(region, stat_output, status) ! ! !USES: ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region logical, intent(in) :: stat_output ! true if stations output is requested ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 29 Nov 2010 - Achim Strunk - ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/output_aerocom_write' integer :: i, imr, jmr, lmr integer :: i1, i2, j1, j2, ilev integer :: istat ! --- begin ------------------------------------- call goLabel(rname) ! grid size call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) imr = i2-i1+1 jmr = j2-j1+1 lmr = levi%nlev ! this here is already accumulated over the time period (day) call collect_budgets( region, status ) IF_NOTOK_RETURN(status=1) ! --------------------- ! divide by accumulator ! --------------------- do i = 1, ntracer_3d mixf(region)%f3d(i)%field = mixf(region)%f3d(i)%field / real( mixf(region)%acct ) end do do i = 1, ntracer_2d mixf(region)%f2d(i)%field = mixf(region)%f2d(i)%field / real( mixf(region)%acct ) end do if (stat_output) then ! stations do istat = 1, nstations if (stat_active(istat)) then do i = 1, ntracer_1dstat statf(region)%f1d(i)%field(istat,:) = statf(region)%f1d(i)%field(istat,:) / real( statf(region)%acct ) end do do i = 1, ntracer_0dstat statf(region)%f0d(i)%field(istat) = statf(region)%f0d(i)%field(istat) / real( statf(region)%acct ) end do endif end do endif select case (trim(aerocom_freq)) case ('hourly') write(gol,'("---> AEROCOM diagnostics: write file for previous hour")'); call goPr case ('daily') write(gol,'("---> AEROCOM diagnostics: write file for previous day")'); call goPr case ('monthly') write(gol,'("---> AEROCOM diagnostics: write file for previous month")'); call goPr end select ! ------------- ! write to file ! ------------- do i=1,ntracer_2d call MDF_Put_Var( mixf(region)%funit, mixf(region)%f2d(i)%varid, mixf(region)%f2d(i)%field, status, start=(/i1,j1/), count=(/imr,jmr/) ) IF_NOTOK_MDF(fid=mixf(region)%funit) end do do i=1,ntracer_3d call MDF_Put_Var( mixf(region)%funit, mixf(region)%f3d(i)%varid, mixf(region)%f3d(i)%field, status, start=(/i1,j1,1/), count=(/imr,jmr,lmr/) ) IF_NOTOK_MDF(fid=mixf(region)%funit) end do call MDF_Close( mixf(region)%funit, status ) IF_NOTOK_MDF(fid=mixf(region)%funit) if (stat_output) then ! stations do istat = 1, nstations if (stat_active(istat)) then do i=1,ntracer_1dstat call MDF_Put_Var( statf(region)%funit, statf(region)%f1d(i)%varid, statf(region)%f1d(i)%field(istat,:), & status, start=(/istat,1/), count=(/1,lmr/) ) IF_NOTOK_MDF(fid=statf(region)%funit) end do do i=1,ntracer_0dstat call MDF_Put_Var( statf(region)%funit, statf(region)%f0d(i)%varid, (/statf(region)%f0d(i)%field(istat)/), & status, start=(/istat/), count=(/1/) ) IF_NOTOK_MDF(fid=statf(region)%funit) end do endif end do call MDF_Close( statf(region)%funit, status ) IF_NOTOK_MDF(fid=statf(region)%funit) endif call goLabel() ; status = 0 end subroutine output_aerocom_write !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: OUTPUT_AEROCOM_STEP ! ! !DESCRIPTION: This is the accumulation routine. It is called following ! the user specification aerocom.dhour in rc-file. It ! evaluates the various diagnostics and does summing. !\\ !\\ ! !INTERFACE: ! subroutine output_aerocom_step( region, dhour, stat_output, status ) ! ! !USES: ! use GO , only : TDate, NewDate, rTotal, operator(-) use Grid , only : FPressure use binas, only : rgas, rol use datetime, only : tau2date use MeteoData, only : sp_dat, lsp_dat, cp_dat, m_dat use MeteoData, only : temper_dat, humid_dat, gph_dat use global_data, only : mass_dat, region_dat use tracer_data, only : chem_dat use dims, only : itaur use chem_param, only : iso4acs, iso4cos, iduacs, iso4ais, ibccos, ibcaii, xmair use chem_param, only : iso4nus, isscos, ino3_a, issacs, iducos, iduaci, nmod use chem_param, only : iducoi, ibcacs, ipomcos, ipomaii, ibcais, ipomacs, ipomais use chem_param, only : imsa, inh4 use chem_param, only : ntracet use mo_aero_m7, only : nsol use optics, only : optics_aop_get use m7_data, only : h2o_mode, rw_mode, rwd_mode ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region integer, intent(in) :: dhour ! this is aerocom.dhour [hours] logical, intent(in) :: stat_output ! true if stations output is requested ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 29 Nov 2010 - Achim Strunk - Creation ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/output_aerocom_step' ! MPI arrays to gather fields. real, dimension(:,:,:,:), pointer :: rm, rmc real, dimension(:,:,:), allocatable, target :: mg real, dimension(:), pointer :: dxyp real, dimension(:,:,:), allocatable :: pres3d integer :: i, j, k, imr, jmr, lmr, lwl, dtime integer :: i1, i2, j1, j2 integer :: ie, je ! indices for subdomain extended with halo cells integer, parameter :: l_halo=1 logical :: polar integer :: istat, imode real :: dens, load_tmp 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 integer :: ncontr, lvec, lct, l, indoffset, nwl integer :: nadd, nadd_max, indoffset_stat, indabs integer :: iadd real :: no3fac, wght, dwght, tabs integer,dimension(6) :: idate_f type(TDate) :: t, t0 real :: time ! --- begin ------------------------------- call goLabel(rname) ! grid size call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2, hasNorthPole=polar ) imr = i2-i1+1 jmr = j2-j1+1 lmr = levi%nlev nullify(rmc) rmc => chem_dat(region)%rm dxyp => region_dat(region)%dxyp rm => mass_dat(region)%rm ! get helping pressure field in 3d ie=i2; je=j2 if (stat_output) then if (stat_halo_glob.eq.1) then CALL UPDATE_HALO( dgrid(region), sp_dat(region)%data(:,:,1), sp_dat(region)%halo, status) IF_NOTOK_RETURN(status=1) ! in the interpolation routine, ! halo cells are only needed at the eastern or northern boundary of the subdomain ie=i2+1 if (.not.polar) je=j2+1 endif endif allocate( pres3d(i1:ie,j1:je,lmr) ) ! fill mid level pressure call FPressure( levi, sp_dat(region)%data(i1:ie,j1:je,1), pres3d, status ) IF_NOTOK_RETURN(status=1) ! dtime is seconds per step ! dtime will be taken as weight for summing up ! (makes it fit for varying step lengths) dtime = dhour * 3600. ! sum up for later averaging mixf (region)%acct = mixf (region)%acct + dtime if (stat_output) then statf(region)%acct = statf(region)%acct + dtime endif ! ---------------------- ! 3D meteo fields ! ---------------------- ! temperature mixf(region)%f3d(temp)%field = mixf(region)%f3d(temp)%field + dtime * temper_dat(region)%data(i1:i2,j1:j2,:) ! specific humidity mixf(region)%f3d(hus)%field = mixf(region)%f3d(hus)%field + dtime * humid_dat(region)%data(i1:i2,j1:j2,:) ! air mass (kg -> kg/m2) do j = j1, j2 mixf(region)%f3d(airmass)%field(:,j,:) = mixf(region)%f3d(airmass)%field(:,j,:) + & dtime * m_dat(region)%data(i1:i2,j,:) / dxyp(j) end do ! pressure (already filled above) mixf(region)%f3d(pressure)%field = mixf(region)%f3d(pressure)%field + dtime * pres3d ! ---------------------- ! cycle over horizontal domain ! ---------------------- do j = j1, j2 do i = i1, i2 ! ----------------------- ! SURFACE FIELDS ! ----------------------- k = 1 ! ---------------------- ! Physical Parameters ! ---------------------- ! surface pressure [Pa] mixf(region)%f2d(ps)%field(i,j) = mixf(region)%f2d(ps)%field(i,j) + dtime * sp_dat(region)%data(i,j,k) ! precipitation ([m s-1] --> [kg m-2 s-1] with density of water `rol`) mixf(region)%f2d(precip)%field(i,j) = mixf(region)%f2d(precip)%field(i,j) + dtime * (cp_dat(region)%data(i,j,k) + lsp_dat(region)%data(i,j,k)) * rol ! density for conversion of aerosol mass densities ( ---> 1/m3 ) ! (conversion factor 1.E-03 is for g --> kg) dens = pres3d(i,j,k) / temper_dat(region)%data(i,j,k) * xmair * 1.E-3 / (m_dat(region)%data(i,j,k) * Rgas) ! ---------------------- ! Surface Concentrations in [kg m-3] ! ---------------------- ! POM: (AIS + ACS + COS + AII) mixf(region)%f2d(sconcoa)%field(i,j) = mixf(region)%f2d(sconcoa)%field(i,j) + dtime * & dens * (rm(i,j,k,iPOMais) + rm(i,j,k,iPOMacs) + rm(i,j,k,iPOMcos) + rm(i,j,k,iPOMaii)) ! BC: (AIS + ACS + COS + AII) mixf(region)%f2d(sconcbc)%field(i,j) = mixf(region)%f2d(sconcbc)%field(i,j) + dtime * & dens * (rm(i,j,k,iBCais) + rm(i,j,k,iBCacs) + rm(i,j,k,iBCcos) + rm(i,j,k,iBCaii)) ! SO4: (NUS + AIS + ACS + COS) mixf(region)%f2d(sconcso4)%field(i,j) = mixf(region)%f2d(sconcso4)%field(i,j) + dtime * & dens * (rm(i,j,k,iSO4nus) + rm(i,j,k,iSO4ais) + rm(i,j,k,iSO4acs) + rm(i,j,k,iSO4cos)) ! Dust: mixf(region)%f2d(sconcdust)%field(i,j) = mixf(region)%f2d(sconcdust)%field(i,j) + dtime * & dens * (rm(i,j,k,iDUacs) + rm(i,j,k,iDUcos) + rm(i,j,k,iDUaci) + rm(i,j,k,iDUcoi)) ! Seasalt: mixf(region)%f2d(sconcss)%field(i,j) = mixf(region)%f2d(sconcss)%field(i,j) + dtime * & dens * (rm(i,j,k,iSSacs) + rm(i,j,k,iSScos)) ! NO3: mixf(region)%f2d(sconcno3)%field(i,j) = mixf(region)%f2d(sconcno3)%field(i,j) + dtime * & dens * rm(i,j,k,iNO3_a) ! NH4: mixf(region)%f2d(sconcnh4)%field(i,j) = mixf(region)%f2d(sconcnh4)%field(i,j) + dtime * & dens * rm(i,j,k,iNH4) ! MSA: mixf(region)%f2d(sconcmsa)%field(i,j) = mixf(region)%f2d(sconcmsa)%field(i,j) + dtime * & dens * rm(i,j,k,iMSA) ! ---------------------- ! Load in [kg m-2] ! ---------------------- do k = 1, lmr ! POM: (AIS + ACS + COS + AII) load_tmp = (rm(i,j,k,iPOMais) + rm(i,j,k,iPOMacs) + rm(i,j,k,iPOMcos) + rm(i,j,k,iPOMaii)) mixf(region)%f2d(loadoa)%field(i,j) = mixf(region)%f2d(loadoa)%field(i,j) + load_tmp * dtime / dxyp(j) ! BC: (AIS + ACS + COS + AII) load_tmp = (rm(i,j,k,iBCais) + rm(i,j,k,iBCacs) + rm(i,j,k,iBCcos) + rm(i,j,k,iBCaii)) mixf(region)%f2d(loadbc)%field(i,j) = mixf(region)%f2d(loadbc)%field(i,j) + load_tmp * dtime / dxyp(j) ! SO4: (NUS + AIS + ACS + COS) load_tmp = (rm(i,j,k,iSO4nus) + rm(i,j,k,iSO4ais) + rm(i,j,k,iSO4acs) + rm(i,j,k,iSO4cos)) mixf(region)%f2d(loadso4)%field(i,j) = mixf(region)%f2d(loadso4)%field(i,j) + load_tmp * dtime / dxyp(j) ! Dust: (ACS + COS + ACI + COI) load_tmp = (rm(i,j,k,iDUacs) + rm(i,j,k,iDUcos) + rm(i,j,k,iDUaci) + rm(i,j,k,iDUcoi)) mixf(region)%f2d(loaddust)%field(i,j) = mixf(region)%f2d(loaddust)%field(i,j) + load_tmp * dtime / dxyp(j) ! Seasalt: (ACS + COS) load_tmp = (rm(i,j,k,iSSacs) + rm(i,j,k,iSScos)) mixf(region)%f2d(loadss)%field(i,j) = mixf(region)%f2d(loadss)%field(i,j) + load_tmp * dtime / dxyp(j) ! NO3: load_tmp = rm(i,j,k,iNO3_a) mixf(region)%f2d(loadno3)%field(i,j) = mixf(region)%f2d(loadno3)%field(i,j) + load_tmp * dtime / dxyp(j) end do ! k end do ! i end do ! j if (stat_output) then ! -------------------------------------- ! stations ! -------------------------------------- if (stat_halo_glob.eq.1) then CALL UPDATE_HALO( dgrid(region), rm, mass_dat(region)%halo, status) IF_NOTOK_RETURN(status=1) CALL UPDATE_HALO( dgrid(region), m_dat(region)%data, m_dat(region)%halo, status) IF_NOTOK_RETURN(status=1) CALL UPDATE_HALO( dgrid(region), gph_dat(region)%data, gph_dat(region)%halo, status) IF_NOTOK_RETURN(status=1) do imode = 1, nsol CALL UPDATE_HALO( dgrid(region), h2o_mode(region,imode)%d3, h2o_mode(region,imode)%halo, status) IF_NOTOK_RETURN(status=1) end do CALL UPDATE_HALO( dgrid(region), temper_dat(region)%data(:,:,:), temper_dat(region)%halo, status) IF_NOTOK_RETURN(status=1) CALL UPDATE_HALO( dgrid(region), humid_dat(region)%data(:,:,:), humid_dat(region)%halo, status) IF_NOTOK_RETURN(status=1) ! sp_dat has been updated earlier CALL UPDATE_HALO( dgrid(region), cp_dat(region)%data(:,:,1), cp_dat(region)%halo, status) IF_NOTOK_RETURN(status=1) CALL UPDATE_HALO( dgrid(region), lsp_dat(region)%data(:,:,1), lsp_dat(region)%halo, status) IF_NOTOK_RETURN(status=1) CALL UPDATE_HALO_IBAND ( dgrid(region), dxyp(:), region_dat(region)%halo, status) endif allocate(volume(i1:ie,j1:je,1:lmr)) do j = j1, je 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) end do do istat = 1, nstations if (stat_active(istat)) then ! masses do i = 1, 18 statf(region)%f1d(i)%field(istat,:) = statf(region)%f1d(i)%field(istat,:) + dtime * & aerocom_do_interp_1d( rm(i1:ie,j1:je,1:lmr,statf(region)%f1d(i)%mf%itm5) / & m_dat(region)%data(i1:ie,j1:je,1:lmr), region, istat ) statf(region)%f0d(i)%field(istat) = statf(region)%f0d(i)%field(istat) + dtime * & aerocom_do_interp_0d( rm(i1:ie,j1:je,1:lmr,statf(region)%f1d(i)%mf%itm5) / & m_dat(region)%data(i1:ie,j1:je,1:lmr), region, istat ) end do ! number densities do i = 19, 19+nmod-1 statf(region)%f1d(i)%field(istat,:) = statf(region)%f1d(i)%field(istat,:) + dtime * & aerocom_do_interp_1d( rm(i1:ie,j1:je,1:lmr,statf(region)%f1d(i)%mf%itm5) / & volume(i1:ie,j1:je,1:lmr), region, istat ) statf(region)%f0d(i)%field(istat) = statf(region)%f0d(i)%field(istat) + dtime * & aerocom_do_interp_0d( rm(i1:ie,j1:je,1:lmr,statf(region)%f1d(i)%mf%itm5) / & volume(i1:ie,j1:je,1:lmr), region, istat ) end do ! aerosol water of the wet modes do i = 19+nmod, 19+nmod+nsol-1 statf(region)%f1d(i)%field(istat,:) = statf(region)%f1d(i)%field(istat,:) + dtime * & aerocom_do_interp_1d( h2o_mode(region,statf(region)%f1d(i)%mf%itm5)%d3(i1:ie,j1:je,1:lmr) / & m_dat(region)%data(i1:ie,j1:je,1:lmr), region, istat ) statf(region)%f0d(i)%field(istat) = statf(region)%f0d(i)%field(istat) + dtime * & aerocom_do_interp_0d( h2o_mode(region,statf(region)%f1d(i)%mf%itm5)%d3(i1:ie,j1:je,1:lmr) / & m_dat(region)%data(i1:ie,j1:je,1:lmr), region, istat ) end do ! nitrate, ammonium and MSA do i= 19+nmod+nsol,19+nmod+nsol+3-1 statf(region)%f1d(i)%field(istat,:) = statf(region)%f1d(i)%field(istat,:) + dtime * & aerocom_do_interp_1d( rm(i1:ie,j1:je,1:lmr,statf(region)%f1d(i)%mf%itm5) / & m_dat(region)%data(i1:ie,j1:je,1:lmr), region, istat ) statf(region)%f0d(i)%field(istat) = statf(region)%f0d(i)%field(istat) + dtime * & aerocom_do_interp_0d( rm(i1:ie,j1:je,1:lmr,statf(region)%f1d(i)%mf%itm5) / & m_dat(region)%data(i1:ie,j1:je,1:lmr), region, istat ) end do ! temperature i=ntracer_1dstat-nextra_1dstat+1 statf(region)%f1d(i)%field(istat,:) = statf(region)%f1d(i)%field(istat,:) + dtime * & aerocom_do_interp_1d( temper_dat(region)%data(i1:ie,j1:je,:), region, istat ) statf(region)%f0d(i)%field(istat) = statf(region)%f0d(i)%field(istat) + dtime * & aerocom_do_interp_0d( temper_dat(region)%data(i1:ie,j1:je,:), region, istat ) ! specific humidity i=i+1 statf(region)%f1d(i)%field(istat,:) = statf(region)%f1d(i)%field(istat,:) + dtime * & aerocom_do_interp_1d( humid_dat(region)%data(i1:ie,j1:je,:), region, istat ) statf(region)%f0d(i)%field(istat) = statf(region)%f0d(i)%field(istat) + dtime * & aerocom_do_interp_0d( humid_dat(region)%data(i1:ie,j1:je,:), region, istat ) ! pressure i=i+1 statf(region)%f1d(i)%field(istat,:) = statf(region)%f1d(i)%field(istat,:) + dtime * & aerocom_do_interp_1d( pres3d(i1:ie,j1:je,:), region, istat ) statf(region)%f0d(i)%field(istat) = statf(region)%f0d(i)%field(istat) + dtime * & aerocom_do_interp_0d( pres3d(i1:ie,j1:je,:), region, istat ) ! surface pressure i=i+1 statf(region)%f0d(i)%field(istat) = statf(region)%f0d(i)%field(istat) + dtime * & aerocom_do_interp_surf( sp_dat(region)%data(i1:ie,j1:je,1), region, istat ) ! total precipitation i=i+1 statf(region)%f0d(i)%field(istat) = statf(region)%f0d(i)%field(istat) + dtime * & aerocom_do_interp_surf( ( cp_dat(region)%data(i1:ie,j1:je,1) + & lsp_dat(region)%data(i1:ie,j1:je,1) ) *rol, region, istat ) endif end do deallocate(volume) endif ! --------------------- ! 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 !ncontr = 10 ncontr = 11 ! Additional fields for surface dry aerosol ! to be used for validation against in-situ data: ! extinction, absorption, asymmetry factor ! fine-mode contribution to extinction, absorption nadd = 5 lvec = imr*jmr*lmr ! allocate output fields for optics_calculate_aop Allocate( aop_output_ext(lvec, size(wdep_out), ncontr ) ) ! extinction Allocate( aop_output_a (lvec, size(wdep_out) ) ) ! single scattering albedo Allocate( aop_output_g (lvec, size(wdep_out) ) ) ! asymmetry factor Allocate( aop_output_add(lvec, size(wdep_out), nadd ) ) ! extinction/absorption due to dry aerosol at surface call optics_aop_get(lvec, region, size(wdep_out), wdep_out, ncontr, aop_output_ext, aop_output_a, aop_output_g, & status, aop_output_add ) IF_NOTOK_RETURN(status=1) ! allocate here result arrays for ext, abs, ssa, asymmetry parameter, additional parameters (PM10) ! ncontr is number of contributors for 'split' allocate( opt_ext( i1:i2, j1:j2, lmr, size(wdep_out), ncontr ) ) ; opt_ext = 0.0 allocate( opt_a ( i1:i2, j1:j2, lmr, size(wdep_out) ) ) ; opt_a = 0.0 allocate( opt_g ( i1:i2, j1:j2, lmr, size(wdep_out) ) ) ; opt_g = 0.0 if (stat_output) then ! halo cells are needed for interpolation of opt_add allocate( opt_add( i1-l_halo:i2+l_halo, j1-l_halo:j2+l_halo, lmr, size(wdep_out), nadd ) ) ; opt_add = 0.0 else allocate( opt_add( i1:i2, j1:j2, lmr, size(wdep_out), nadd ) ) ; opt_add = 0.0 endif ! --------------------------------- ! unpack results from calculate_aop ! --------------------------------- do lwl = 1, size(wdep_out) if( wdep_out(lwl)%split ) then ! fill the array for the extinction coefficients of contributors do lct = 1, ncontr opt_ext(:,:,:,lwl,lct) = reshape( aop_output_ext(:,lwl,lct), (/imr,jmr,lmr/) ) end do else ! fill only total extinction coefficients opt_ext(:,:,:,lwl,1) = reshape( aop_output_ext(:,lwl,1), (/imr,jmr,lmr/) ) endif opt_a(:,:,:,lwl) = reshape( aop_output_a(:,lwl),(/imr,jmr,lmr/) ) opt_g(:,:,:,lwl) = reshape( aop_output_g(:,lwl),(/imr,jmr,lmr/) ) if ( wdep_out(lwl)%insitu ) then ! additional fields for split do lct = 1, nadd opt_add(i1:i2,j1:j2,:,lwl,lct) = reshape( aop_output_add(:,lwl,lct),(/imr,jmr,lmr/) ) end do endif end do ! lwl ! free temporary arrays for results from calculate_aop deallocate( aop_output_ext ) deallocate( aop_output_a ) deallocate( aop_output_g ) deallocate( aop_output_add ) ! the following sections assume that for 550nm there is splitted information ! available and that there is *NO* split for other wavelengths! if( count( (wdep_out(:)%wl == 0.550) .and. wdep_out(:)%split ) /= 1 ) then write(gol,*) 'user_output_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( wdep_out(lwl)%split ) then ! for 550nm: ! 1) AODs (35 - 44) ! fill for contributors (total, SO4, BC, POM, SS, DU, NO3, H2O, fine, fine dust, fine SS) ! 2) Absorption for 550nm (45) ! Extinction is here the sum of scattering and absorption and ! the scattering part is given by the albedo (SSA). Thus the ! remainder is absorption: Abs = Ext * (1 - SSA) indoffset = od550aer allocate(temp2d(i1:i2,j1:j2)) do lct = 1, ncontr temp2d = 0.0 do j = j1, j2 do i = i1, i2 ! multiply with height elements and sum up tabs = 0.0 do l = 1, lmr temp2d(i,j) = temp2d(i,j) + opt_ext(i,j,l,lwl,lct) * & (gph_dat(region)%data(i,j,l+1) - & gph_dat(region)%data(i,j,l )) if( lct == 1 ) then ! Absorption: do this only once for the totals tabs = tabs + opt_ext(i,j,l,lwl,lct) * (1.0 - opt_a(i,j,l,lwl)) * & (gph_dat(region)%data(i,j,l+1) - gph_dat(region)%data(i,j,l)) end if end do ! Absorption: do this only once for the totals if( lct == 1 ) then mixf(region)%f2d(abs550aer)%field(i,j) = & mixf(region)%f2d(abs550aer)%field(i,j) + tabs * dtime end if end do end do ! this here is AODs for different contributors, splitted by volumes mixf(region)%f2d(indoffset+lct-1)%field = & mixf(region)%f2d(indoffset+lct-1)%field + temp2d * dtime end do deallocate( temp2d ) ! Asymmetry Parameter, weighted by AOD and single scattering albedo at each layer allocate(temp2d(i1:i2,j1:j2)) ; temp2d = 0.0 do j = j1, j2 do i = i1, i2 wght = 0.0 do l = 1, lmr dwght = opt_ext(i,j,l,lwl,1) * (gph_dat(region)%data(i,j,l+1) - gph_dat(region)%data(i,j,l )) * opt_a(i,j,l,lwl) temp2d(i,j) = temp2d(i,j) + opt_g(i,j,l,lwl) * dwght wght = wght + dwght end do temp2d(i,j) = temp2d(i,j) / wght end do end do mixf(region)%f2d(asyaer)%field = mixf(region)%f2d(asyaer)%field + temp2d * dtime deallocate(temp2d) ! Surface Ambient Aerosol Extinction mixf(region)%f2d(ec550aer)%field = & mixf(region)%f2d(ec550aer)%field + opt_ext(:,:,1,lwl,1) * dtime else ! for 440/870/350 nm: ! 1) fill total AOD and AAOD and write to containers ! 2) dump AOD fields if ( wdep_out(lwl)%wl == 0.440 ) then indoffset = od440aer indabs = abs440aer elseif ( wdep_out(lwl)%wl == 0.550 ) then indoffset = od550aer indabs = abs550aer elseif ( wdep_out(lwl)%wl == 0.870 ) then indoffset = od870aer indabs = abs870aer elseif ( wdep_out(lwl)%wl == 0.350 ) then indoffset = od350aer indabs = abs350aer else write(gol,*) 'user_output_aerocom: wrong wavelength given for aerocom diagnostics' ; call goErr TRACEBACK status=1; return end if ! fill AODs allocate(temp2d(i1:i2,j1:j2)) temp2d = 0.0 do j = j1, j2 do i = i1, i2 ! multiply with height elements and sum up tabs = 0.0 do l = 1, lmr temp2d(i,j) = temp2d(i,j) + opt_ext(i,j,l,lwl,1) * & (gph_dat(region)%data(i,j,l+1) - & gph_dat(region)%data(i,j,l )) tabs = tabs + opt_ext(i,j,l,lwl,1) * (1.0 - opt_a(i,j,l,lwl)) * & (gph_dat(region)%data(i,j,l+1) - gph_dat(region)%data(i,j,l)) end do mixf(region)%f2d(indabs)%field(i,j) = & mixf(region)%f2d(indabs)%field(i,j) + tabs * dtime end do end do ! write to container mixf(region)%f2d(indoffset)%field = & mixf(region)%f2d(indoffset)%field + temp2d * dtime deallocate(temp2d) endif ! 3-D extinction and absorption (m-1) if ( aai_output ) then if ( wdep_out(lwl)%wl == 0.550 .or. wdep_out(lwl)%wl == 0.350 ) then if ( wdep_out(lwl)%wl == 0.550 ) then indoffset = ec5503Daer elseif ( wdep_out(lwl)%wl == 0.350 ) then indoffset = ec3503Daer endif mixf(region)%f3d(indoffset)%field = & mixf(region)%f3d(indoffset)%field + opt_ext(:,:,:,lwl,1) * dtime mixf(region)%f3d(indoffset+1)%field = & mixf(region)%f3d(indoffset+1)%field + opt_ext(:,:,:,lwl,1) * (1.0 - opt_a(:,:,:,lwl)) * dtime endif endif if ( wdep_out(lwl)%insitu ) then if ( wdep_out(lwl)%wl == 0.440 ) then indoffset = ec440dryaer indoffset_stat = ec440dryaer_stat nadd_max=2 elseif ( wdep_out(lwl)%wl == 0.550 ) then indoffset = ec550dryaer indoffset_stat = ec550dryaer_stat nadd_max=5 elseif ( wdep_out(lwl)%wl == 0.870 ) then indoffset = ec870dryaer indoffset_stat = ec870dryaer_stat nadd_max=2 else write(gol,*) 'user_output_aerocom: wrong wavelength given for aerocom diagnostics' ; call goErr TRACEBACK status=1; return end if ! Surface dry aerosol extinction & absorption do lct = 1, nadd_max mixf(region)%f2d(indoffset+lct-1)%field = mixf(region)%f2d(indoffset+lct-1)%field + & opt_add(i1:i2,j1:j2,1,lwl,lct) * dtime enddo if (stat_output) then ! -------------------------------------- ! stations ! -------------------------------------- if (stat_halo_glob.eq.1) then do lct = 1, nadd_max CALL UPDATE_HALO( dgrid(region), opt_add(:,:,:,lwl,lct), l_halo, status) IF_NOTOK_RETURN(status=1) enddo endif do istat = 1, nstations if (stat_active(istat)) then do lct = 1, nadd_max statf(region)%f1d(indoffset_stat+lct-1)%field(istat,:) = & statf(region)%f1d(indoffset_stat+lct-1)%field(istat,:) + & aerocom_do_interp_1d( opt_add(i1:ie,j1:je,:,lwl,lct), region, istat ) * dtime statf(region)%f0d(indoffset_stat+lct-1)%field(istat) = & statf(region)%f0d(indoffset_stat+lct-1)%field(istat) + & aerocom_do_interp_0d( opt_add(i1:ie,j1:je,:,lwl,lct), region, istat ) * dtime enddo endif end do endif endif end do ! done deallocate( opt_ext, opt_a, opt_g, opt_add ) nullify(rm) nullify(dxyp) deallocate( pres3d ) call goLabel() ; status = 0 end subroutine output_aerocom_step !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: COLLECT_BUDGETS ! ! !DESCRIPTION: This routine does book keeping of the budget values in ! order to extract daily contributions to ! emissions / dry deposition / wet deposition. !\\ !\\ ! !INTERFACE: ! subroutine collect_budgets(region, status) ! ! !USES: ! #ifndef without_chemistry use ebischeme, only : buddry_dat => buddep_dat #endif use wet_deposition, only : buddep_dat use emission_data, only : budemi_dat use budget_global, only : nbud_vg use global_data, only : region_dat use chem_param, only : iducoi, iduacs, iducos, iduaci, isscos, issacs use chem_param, only : ibccos, ibcaii, ibcacs, ibcais, ino3_a, xmair use chem_param, only : iso4nus, iso4acs, iso4cos, iso4ais use chem_param, only : ipomcos, ipomaii, ipomacs, ipomais use chem_param, only : idms, iso2, ntracet, ntrace use chem_param, only : xmso2, xmso4, xmdms, xmpom, xmbc, xmdust, xmnacl #ifndef without_sedimentation use sedimentation, only : buddep_m7_dat use sedimentation, only : nsed, ised #endif ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 29 Nov 2010 - Achim Strunk - ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/collect_budgets' real, dimension(:,:,:,:), allocatable :: collect4d real, dimension(:,:,:), allocatable :: collect3d real, dimension(:,:), allocatable :: tmpbud2d real, dimension(:), pointer :: dxyp integer :: imr, jmr, lmr, j integer :: i1, i2, j1, j2 !>>> TvN integer :: l, n !<<< TvN ! --- begin ------------------------------- call goLabel(rname) ! grid size call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) imr=i2-i1+1 jmr=j2-j1+1 lmr = levi%nlev ! this is area element dxyp => region_dat(region)%dxyp ! -------------- ! EMISSIONS ! -------------- ! collect emission budget allocate( collect4d(i1:i2,j1:j2,nbud_vg,ntracet) ); collect4d = 0.0 ! emissions are in [mole] ! convert first to kilomole per area [k-mole/m2] #ifndef without_emission do j = j1, j2 collect4d(:,j,:,:) = budemi_dat(region)%emi(i1:i2,j,:,:) / dxyp(j) * 1.E-03 end do #endif allocate( tmpbud2d(i1:i2,j1:j2) ) ! on the way: convert from kilomole/m2 to kg/m2 via molar mass [g/mole] ! POM (AIS + ACS + COS + AII) tmpbud2d = ( sum(collect4d(:,:,:,iPOMais),3) + sum(collect4d(:,:,:,iPOMacs),3) + & sum(collect4d(:,:,:,iPOMcos),3) + sum(collect4d(:,:,:,iPOMaii),3) ) * xmpom mixf(region)%f2d(emioa)%field = tmpbud2d - emission(region)%f2dslast(:,:,1) emission(region)%f2dslast(:,:,1) = tmpbud2d ! BC (AIS + ACS + COS + AII) tmpbud2d = ( sum(collect4d(:,:,:,iBCais ),3) + sum(collect4d(:,:,:,iBCacs ),3) + & sum(collect4d(:,:,:,iBCcos ),3) + sum(collect4d(:,:,:,iBCaii ),3) ) * xmbc mixf(region)%f2d(emibc)%field = tmpbud2d - emission(region)%f2dslast(:,:,2) emission(region)%f2dslast(:,:,2) = tmpbud2d ! SO2 tmpbud2d = sum(collect4d(:,:,:,iSO2 ),3) * xmso2 mixf(region)%f2d(emiso2)%field = tmpbud2d - emission(region)%f2dslast(:,:,3) emission(region)%f2dslast(:,:,3) = tmpbud2d ! SO4 (NUS + AIS + ACS + COS) tmpbud2d = ( sum(collect4d(:,:,:,iSO4nus),3) + sum(collect4d(:,:,:,iSO4ais),3) + & sum(collect4d(:,:,:,iSO4acs),3) + sum(collect4d(:,:,:,iSO4cos),3) ) * xmso4 mixf(region)%f2d(emiso4)%field = tmpbud2d - emission(region)%f2dslast(:,:,4) emission(region)%f2dslast(:,:,4) = tmpbud2d ! Dust (ACS + COS + ACI + COI) tmpbud2d = ( sum(collect4d(:,:,:,iDUacs ),3) + sum(collect4d(:,:,:,iDUcos ),3) + & sum(collect4d(:,:,:,iDUaci ),3) + sum(collect4d(:,:,:,iDUcoi ),3) ) * xmdust mixf(region)%f2d(emidust)%field = tmpbud2d - emission(region)%f2dslast(:,:,5) emission(region)%f2dslast(:,:,5) = tmpbud2d ! DMS tmpbud2d = sum(collect4d(:,:,:,iDMS ),3) * xmdms mixf(region)%f2d(emidms)%field = tmpbud2d - emission(region)%f2dslast(:,:,6) emission(region)%f2dslast(:,:,6) = tmpbud2d ! Seasalt: (ACS + COS) tmpbud2d = ( sum(collect4d(:,:,:,iSSacs ),3) + sum(collect4d(:,:,:,iSScos ),3) ) * xmnacl mixf(region)%f2d(emiss)%field = tmpbud2d - emission(region)%f2dslast(:,:,7) emission(region)%f2dslast(:,:,7) = tmpbud2d deallocate( tmpbud2d ) deallocate( collect4d ) ! -------------- ! DRY DEPOSITION ! -------------- allocate( collect3d(i1:i2,j1:j2,ntrace) ); collect3d = 0.0 ! deposition is in [mole] ! convert first to kilomole per area [k-mole/m2] do j = j1, j2 #ifndef without_chemistry collect3d(:,j,:) = buddry_dat(region)%dry(i1:i2,j,:) / dxyp(j) * 1.E-3 #endif ! Add sedimentation at the surface to dry depostion ! Sedimentation budgets have to be summed in the vertical ! to get the surface contribution. #ifndef without_sedimentation do l = 1, nbud_vg do n = 1, nsed collect3d(:,j,ised(n)) = collect3d(:,j,ised(n)) + & buddep_m7_dat(region)%sed(i1:i2,j,l,n) / dxyp(j) * 1.E-3 end do end do #endif end do allocate( tmpbud2d(i1:i2,j1:j2) ) ! on the way: convert from kilomole/m2 to kg/m2 via molar mass [g/mole] ! SO2 tmpbud2d = collect3d(:,:,iSO2) * xmso2 mixf(region)%f2d(dryso2)%field = tmpbud2d - drydepos(region)%f2dslast(:,:,1) drydepos(region)%f2dslast(:,:,1) = tmpbud2d ! POM (AIS + ACS + COS + AII) tmpbud2d = ( collect3d(:,:,iPOMais ) + collect3d(:,:,iPOMacs ) + & collect3d(:,:,iPOMcos ) + collect3d(:,:,iPOMaii ) ) * xmpom mixf(region)%f2d(dryoa)%field = tmpbud2d - drydepos(region)%f2dslast(:,:,2) drydepos(region)%f2dslast(:,:,2) = tmpbud2d ! BC (AIS + ACS + COS + AII) tmpbud2d = ( collect3d(:,:,iBCais ) + collect3d(:,:,iBCacs ) + & collect3d(:,:,iBCcos ) + collect3d(:,:,iBCaii ) ) * xmbc mixf(region)%f2d(drybc)%field = tmpbud2d - drydepos(region)%f2dslast(:,:,3) drydepos(region)%f2dslast(:,:,3) = tmpbud2d ! SO4 (NUS + AIS + ACS + COS) tmpbud2d = ( collect3d(:,:,iSO4nus) + collect3d(:,:,iSO4ais) + & collect3d(:,:,iSO4acs) + collect3d(:,:,iSO4cos) ) * xmso4 mixf(region)%f2d(dryso4)%field = tmpbud2d - drydepos(region)%f2dslast(:,:,4) drydepos(region)%f2dslast(:,:,4) = tmpbud2d ! Dust (ACS + COS + ACI + COI) tmpbud2d = ( collect3d(:,:,iDUacs) + collect3d(:,:,iDUcos) + & collect3d(:,:,iDUaci) + collect3d(:,:,iDUcoi) ) * xmdust mixf(region)%f2d(drydust)%field = tmpbud2d - drydepos(region)%f2dslast(:,:,5) drydepos(region)%f2dslast(:,:,5) = tmpbud2d ! DMS ! DMS is NOT deposited in TM5 --> the fields will contain zeros tmpbud2d = collect3d(:,:,iDMS) * xmdms mixf(region)%f2d(drydms)%field = tmpbud2d - drydepos(region)%f2dslast(:,:,6) drydepos(region)%f2dslast(:,:,6) = tmpbud2d ! Seasalt: (ACS + COS) tmpbud2d = ( collect3d(:,:,iSSacs) + collect3d(:,:,iSScos) ) * xmnacl mixf(region)%f2d(dryss)%field = tmpbud2d - drydepos(region)%f2dslast(:,:,7) drydepos(region)%f2dslast(:,:,7) = tmpbud2d deallocate( tmpbud2d ) deallocate( collect3d ) ! -------------- ! WET DEPOSITION ! -------------- allocate( collect4d (i1:i2,j1:j2,nbud_vg,ntracet) ); collect4d = 0.0 ! deposition is in [mole] ! convert first to kilomole per area [k-mole/m2] do j = j1, j2 collect4d(:,j,:,:) = ( buddep_dat(region)%lsp(i1:i2,j,:,:) + buddep_dat(region)%cp(i1:i2,j,:,:) ) / dxyp(j) * 1.E-3 end do allocate( tmpbud2d(i1:i2,j1:j2) ) ! on the way: convert from kilomole/m2 to kg/m2 via molar mass [g/mole] ! POM (AIS + ACS + COS + AII) tmpbud2d = ( sum(collect4d(:,:,:,iPOMais),3) + sum(collect4d(:,:,:,iPOMacs),3) + & sum(collect4d(:,:,:,iPOMcos),3) + sum(collect4d(:,:,:,iPOMaii),3) ) * xmpom mixf(region)%f2d(wetoa)%field = tmpbud2d - wetdepos(region)%f2dslast(:,:,1) wetdepos(region)%f2dslast(:,:,1) = tmpbud2d ! BC (AIS + ACS + COS + AII) tmpbud2d = ( sum(collect4d(:,:,:,iBCais ),3) + sum(collect4d(:,:,:,iBCacs ),3) + & sum(collect4d(:,:,:,iBCcos ),3) + sum(collect4d(:,:,:,iBCaii ),3) ) * xmbc mixf(region)%f2d(wetbc)%field = tmpbud2d - wetdepos(region)%f2dslast(:,:,2) wetdepos(region)%f2dslast(:,:,2) = tmpbud2d ! SO2 tmpbud2d = sum(collect4d(:,:,:,iSO2 ),3) * xmso2 mixf(region)%f2d(wetso2)%field = tmpbud2d - wetdepos(region)%f2dslast(:,:,3) wetdepos(region)%f2dslast(:,:,3) = tmpbud2d ! SO4 (NUS + AIS + ACS + COS) tmpbud2d = ( sum(collect4d(:,:,:,iSO4nus),3) + sum(collect4d(:,:,:,iSO4ais),3) + & sum(collect4d(:,:,:,iSO4acs),3) + sum(collect4d(:,:,:,iSO4cos),3) ) * xmso4 mixf(region)%f2d(wetso4)%field = tmpbud2d - wetdepos(region)%f2dslast(:,:,4) wetdepos(region)%f2dslast(:,:,4) = tmpbud2d ! Dust (ACS + COS + ACI + COI) tmpbud2d = ( sum(collect4d(:,:,:,iDUacs ),3) + sum(collect4d(:,:,:,iDUcos ),3) + & sum(collect4d(:,:,:,iDUaci ),3) + sum(collect4d(:,:,:,iDUcoi ),3) ) * xmdust mixf(region)%f2d(wetdust)%field = tmpbud2d - wetdepos(region)%f2dslast(:,:,5) wetdepos(region)%f2dslast(:,:,5) = tmpbud2d ! DMS ! DMS is NOT deposited in TM5 --> the fields will contain zeros tmpbud2d = sum(collect4d(:,:,:,iDMS ),3) * xmdms mixf(region)%f2d(wetdms)%field = tmpbud2d - wetdepos(region)%f2dslast(:,:,6) wetdepos(region)%f2dslast(:,:,6) = tmpbud2d ! Seasalt: (ACS + COS) tmpbud2d = ( sum(collect4d(:,:,:,iSSacs ),3) + sum(collect4d(:,:,:,iSScos ),3) ) * xmnacl mixf(region)%f2d(wetss)%field = tmpbud2d - wetdepos(region)%f2dslast(:,:,7) wetdepos(region)%f2dslast(:,:,7) = tmpbud2d deallocate( tmpbud2d ) deallocate( collect4d ) nullify(dxyp) call goLabel() ; status = 0 end subroutine collect_budgets !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: AEROCOM_GET_STATINTP ! ! !DESCRIPTION: This routine provides interpolation coefficients to ! (yet another) interpolation routine for the stations. ! The interpolation routine is below ( aerocom_do_interp ). !\\ ! !INTERFACE: ! subroutine aerocom_get_statintp( region, status ) ! ! !USES: ! use dims, only : dy, dx, yref, xref, xbeg, xend, ybeg, yend, xcyc use dims, only : im, jm, lm use Meteodata, only : gph_dat use partools, only : par_reduce ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 29 Nov 2010 - Achim Strunk - Creation ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/aerocom_get_statintp' real :: rt, sheight, dxr, dyr real, dimension(0:lm(region)) :: gheight integer :: i, l integer :: i1, i2, j1, j2 ! --- begin ------------------------------- call goLabel(rname) call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) dyr = dy/yref(region) dxr = dx/xref(region) do i = 1, nstations ! i-coordinates and weights rt = (statlist(i)%statlon - float(xbeg(region)))/dxr + 0.5 sintp(region)%ixl(i) = int(rt) sintp(region)%wixr(i) = rt - sintp(region)%ixl(i) sintp(region)%wixl(i) = 1.0 - sintp(region)%wixr(i) if ( sintp(region)%ixl(i) == 0 ) then sintp(region)%ixl(i) = im(region) endif sintp(region)%ixr(i) = sintp(region)%ixl(i) + 1 ! j-coordinates and weights rt = (statlist(i)%statlat - float(ybeg(region)))/dyr + 0.5 sintp(region)%jyl(i) = int(rt) sintp(region)%jyr(i) = sintp(region)%jyl(i) + 1 sintp(region)%wjyr(i) = rt - sintp(region)%jyl(i) sintp(region)%wjyl(i) = 1.0 - sintp(region)%wjyr(i) if ( sintp(region)%jyl(i) == 0 ) sintp(region)%jyl(i) = 1 if ( sintp(region)%jyr(i) == jm(region)+1 ) sintp(region)%jyr(i) = jm(region) ! check if left/lower cell of interpolation calculation is in the processor domain ! ensures that each station is handled by a single processor if ( (sintp(region)%ixl(i).ge.i1).and.(sintp(region)%ixl(i).le.i2).and. & (sintp(region)%jyl(i).ge.j1).and.(sintp(region)%jyl(i).le.j2) ) then stat_active(i)=.true. if ( (sintp(region)%ixr(i).eq.i2+1).or.(sintp(region)%jyr(i).eq.j2+1) ) then stat_halo_loc=1 endif endif enddo call PAR_REDUCE( stat_halo_loc, 'max', stat_halo_glob, status, all=.true.) if (stat_halo_glob.eq.1) then CALL UPDATE_HALO( dgrid(region), gph_dat(region)%data, gph_dat(region)%halo, status) IF_NOTOK_RETURN(status=1) endif do i = 1, nstations if ( stat_active(i) ) then ! -------- ! interpolate gph to station, use half level (center of box) & surface ! gph is from 1 -> lm(region)+1 ! -------- gheight(0) = sintp(region)%wixl(i) * sintp(region)%wjyl(i) * gph_dat(region)%data(sintp(region)%ixl(i),sintp(region)%jyl(i),1) + & sintp(region)%wixr(i) * sintp(region)%wjyl(i) * gph_dat(region)%data(sintp(region)%ixr(i),sintp(region)%jyl(i),1) + & sintp(region)%wixr(i) * sintp(region)%wjyr(i) * gph_dat(region)%data(sintp(region)%ixr(i),sintp(region)%jyr(i),1) + & sintp(region)%wixl(i) * sintp(region)%wjyr(i) * gph_dat(region)%data(sintp(region)%ixl(i),sintp(region)%jyr(i),1) do l = 1, lm(region) gheight(l) = sintp(region)%wixl(i) * sintp(region)%wjyl(i) * sum(gph_dat(region)%data(sintp(region)%ixl(i),sintp(region)%jyl(i),l:l+1))/2. + & sintp(region)%wixr(i) * sintp(region)%wjyl(i) * sum(gph_dat(region)%data(sintp(region)%ixr(i),sintp(region)%jyl(i),l:l+1))/2. + & sintp(region)%wixr(i) * sintp(region)%wjyr(i) * sum(gph_dat(region)%data(sintp(region)%ixr(i),sintp(region)%jyr(i),l:l+1))/2. + & sintp(region)%wixl(i) * sintp(region)%wjyr(i) * sum(gph_dat(region)%data(sintp(region)%ixl(i),sintp(region)%jyr(i),l:l+1))/2. end do sintp(region)%tght(i) = gheight(0) ! -------- ! determine layer of station ! -------- ! copy temporary the height of station sheight = statlist(i)%statalt ! do not allow sampleheight below model surface ! -> force station height to model surface sheight = max(sheight, gheight(0)) do l = 1, lm(region) if( gheight(l) > sheight ) exit end do ! success?? if( l > lm(region) ) then ! not successful, station higher than model (strange!!!) write(gol,'("Error in retrieving vertical layer of station")'); call goErr TRACEBACK status=1 return endif ! restrict lower index to surface sintp(region)%lzl(i) = max(l-1,1) sintp(region)%lzr(i) = l sintp(region)%wlzl(i) = (gheight(l) - sheight) / (gheight(l) - gheight(l-1)) sintp(region)%wlzr(i) = 1.0 - sintp(region)%wlzl(i) endif end do call goLabel() ; status = 0 end subroutine aerocom_get_statintp !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !FUNCTION: AEROCOM_DO_INTERP_1D ! ! !DESCRIPTION: Interpolation of a given 3d field, using module wide ! interpolation coeffients and indices for stations (sintp). !\\ !\\ ! !INTERFACE: ! function aerocom_do_interp_1d( f3d, region, istat ) ! ! !INPUT PARAMETERS: ! real, dimension(:,:,:), intent(in) :: f3d integer, intent(in) :: region integer, intent(in) :: istat ! ! !RETURN VALUE: ! real,dimension (levi%nlev) :: aerocom_do_interp_1d ! ! !REVISION HISTORY: ! 6 Feb 2011 - Achim Strunk - ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/aerocom_do_interp_1d' integer :: i integer :: i1, i2, j1, j2 integer :: il,ir,jl,jr ! --- begin ------------------------------- call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) i = istat ! shift index in the input array il=sintp(region)%ixl(i) - i1 + 1 ir=sintp(region)%ixr(i) - i1 + 1 jl=sintp(region)%jyl(i) - j1 + 1 jr=sintp(region)%jyr(i) - j1 + 1 ! ordering ! i = llllrrrr, j = llrrllrr aerocom_do_interp_1d(:) = & sintp(region)%wixl(i) * sintp(region)%wjyl(i) * f3d(il,jl,:) + & sintp(region)%wixl(i) * sintp(region)%wjyr(i) * f3d(il,jr,:) + & sintp(region)%wixr(i) * sintp(region)%wjyl(i) * f3d(ir,jl,:) + & sintp(region)%wixr(i) * sintp(region)%wjyr(i) * f3d(ir,jr,:) return end function aerocom_do_interp_1d !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !FUNCTION: AEROCOM_DO_INTERP_0D ! ! !DESCRIPTION: Interpolation of a given 3d field, using module wide ! interpolation coeffients and indices for stations (sintp). !\\ !\\ ! !INTERFACE: ! function aerocom_do_interp_0d( f3d, region, istat ) ! ! !INPUT PARAMETERS: ! real, dimension(:,:,:), intent(in) :: f3d integer, intent(in) :: region integer, intent(in) :: istat ! ! !RETURN VALUE: ! real :: aerocom_do_interp_0d ! ! !REVISION HISTORY: ! 6 Feb 2011 - Achim Strunk - ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/aerocom_do_interp_0d' integer :: i integer :: i1, i2, j1, j2 integer :: il,ir,jl,jr ! --- begin ------------------------------- call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) i = istat ! shift index in the input array il=sintp(region)%ixl(i) - i1 + 1 ir=sintp(region)%ixr(i) - i1 + 1 jl=sintp(region)%jyl(i) - j1 + 1 jr=sintp(region)%jyr(i) - j1 + 1 ! ordering ! i = llllrrrr, j = llrrllrr, l = lrlrlrlr aerocom_do_interp_0d = & sintp(region)%wixl(i) * sintp(region)%wjyl(i) * sintp(region)%wlzl(i) * f3d(il,jl,sintp(region)%lzl(i)) + & sintp(region)%wixl(i) * sintp(region)%wjyl(i) * sintp(region)%wlzr(i) * f3d(il,jl,sintp(region)%lzr(i)) + & sintp(region)%wixl(i) * sintp(region)%wjyr(i) * sintp(region)%wlzl(i) * f3d(il,jr,sintp(region)%lzl(i)) + & sintp(region)%wixl(i) * sintp(region)%wjyr(i) * sintp(region)%wlzr(i) * f3d(il,jr,sintp(region)%lzr(i)) + & sintp(region)%wixr(i) * sintp(region)%wjyl(i) * sintp(region)%wlzl(i) * f3d(ir,jl,sintp(region)%lzl(i)) + & sintp(region)%wixr(i) * sintp(region)%wjyl(i) * sintp(region)%wlzr(i) * f3d(ir,jl,sintp(region)%lzr(i)) + & sintp(region)%wixr(i) * sintp(region)%wjyr(i) * sintp(region)%wlzl(i) * f3d(ir,jr,sintp(region)%lzl(i)) + & sintp(region)%wixr(i) * sintp(region)%wjyr(i) * sintp(region)%wlzr(i) * f3d(ir,jr,sintp(region)%lzr(i)) return end function aerocom_do_interp_0d !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !FUNCTION: AEROCOM_DO_INTERP_SURF ! ! !DESCRIPTION: Interpolation of a given surface field, using module wide ! interpolation coeffients and indices for stations (sintp). !\\ !\\ ! !INTERFACE: ! function aerocom_do_interp_surf( f2d, region, istat ) ! ! !INPUT PARAMETERS: ! real, dimension(:,:), intent(in) :: f2d integer, intent(in) :: region integer, intent(in) :: istat ! ! !RETURN VALUE: ! real :: aerocom_do_interp_surf ! ! !REVISION HISTORY: ! 6 Feb 2011 - Achim Strunk - ! ! !REMARKS: ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/aerocom_do_interp_surf' integer :: i integer :: i1, i2, j1, j2 integer :: il,ir,jl,jr ! --- begin ------------------------------- call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) i = istat ! shift index in the input array il=sintp(region)%ixl(i) - i1 + 1 ir=sintp(region)%ixr(i) - i1 + 1 jl=sintp(region)%jyl(i) - j1 + 1 jr=sintp(region)%jyr(i) - j1 + 1 ! ordering ! i = llllrrrr, j = llrrllrr aerocom_do_interp_surf = & sintp(region)%wixl(i) * sintp(region)%wjyl(i) * f2d(il,jl) + & sintp(region)%wixl(i) * sintp(region)%wjyr(i) * f2d(il,jr) + & sintp(region)%wixr(i) * sintp(region)%wjyl(i) * f2d(ir,jl) + & sintp(region)%wixr(i) * sintp(region)%wjyr(i) * f2d(ir,jr) return end function aerocom_do_interp_surf !EOC end module user_output_aerocom