diff --git a/AOTnPOD_ml.f90 b/AOTnPOD_ml.f90 old mode 100644 new mode 100755 diff --git a/Advection_ml.f90 b/Advection_ml.f90 old mode 100644 new mode 100755 diff --git a/Aero_water_ml.f90 b/Aero_water_ml.f90 old mode 100644 new mode 100755 diff --git a/AirEmis_ml.f90 b/AirEmis_ml.f90 old mode 100644 new mode 100755 diff --git a/Ammonium_ml.f90 b/Ammonium_ml.f90 old mode 100644 new mode 100755 diff --git a/Aqueous_n_WetDep_ml.f90 b/Aqueous_n_WetDep_ml.f90 old mode 100644 new mode 100755 diff --git a/Biogenics_ml.f90 b/Biogenics_ml.f90 old mode 100644 new mode 100755 index dd01205..0549503 --- a/Biogenics_ml.f90 +++ b/Biogenics_ml.f90 @@ -69,6 +69,7 @@ module Biogenics_ml use ModelConstants_ml, only : NPROC, MasterProc, TINY, & USE_PFT_MAPS, NLANDUSEMAX, IOU_INST, & KT => KCHEMTOP, KG => KMAX_MID, & + EURO_SOILNOX_DEPSCALE, & DEBUG_BIO, BVOC_USED, MasterProc use NetCDF_ml, only : ReadField_CDF, printCDF use OwnDataTypes_ml, only : Deriv, TXTLEN_SHORT @@ -113,10 +114,6 @@ module Biogenics_ml AnnualNdep, & ! N-dep in mgN/m2/ SoilNOx, SoilNH3 - ! To avoid pre-runs of the model for scenario years, we assume that changes - ! can ve approximated by EU N emission changes - real, public, save :: Ndep_trends = 1.0 ! for scaling assumed soil N-dep - ! Set true if LCC read from e.g. EMEP_EuroBVOC.nc: ! (Currently for 1st four LCC, CF, DF, BF, NF) @@ -644,7 +641,7 @@ subroutine Set_SoilNOx() if( DEBUG_SOILNOX .and. debug_proc ) then write(*,*)"Biogenic_ml DEBUG_SOILNOX EURO: ",& current_date%day, current_date%hour, current_date%seconds,& - USE_EURO_SOILNOX, Ndep_trends + USE_EURO_SOILNOX, EURO_SOILNOX_DEPSCALE end if if ( .not. USE_EURO_SOILNOX ) return ! and fSW has been set to 1. at start @@ -671,7 +668,7 @@ subroutine Set_SoilNOx() ! We use a factor normalised to 1.0 at 5000 mgN/m2/a fn = AnnualNdep(i,j)/5000.0 ! scale for now - fn = fn * Ndep_trends ! For e.g. 2030, see Emissions_ml + fn = fn * EURO_SOILNOX_DEPSCALE ! See ModelConstants_ml ftn = ft * fn * hfac diff --git a/BoundaryConditions_ml.f90 b/BoundaryConditions_ml.f90 old mode 100644 new mode 100755 diff --git a/Chem_ml.f90 b/Chem_ml.f90 old mode 100644 new mode 100755 diff --git a/Country_ml.f90 b/Country_ml.f90 old mode 100644 new mode 100755 index b1e36e2..04e6717 --- a/Country_ml.f90 +++ b/Country_ml.f90 @@ -203,6 +203,7 @@ module Country_ml integer, parameter, public :: IC_TAIW = 234 ! Taiwan integer, parameter, public :: IC_THAI = 235 ! Thailand integer, parameter, public :: IC_VIET = 236 ! Vietnam +integer, parameter, public :: IC_EGYP = 237 ! Egypt ! extra subdivisions of ship emissions into shipping categories: ! Baltic Sea (30) @@ -514,6 +515,7 @@ subroutine Country_Init() Country(IC_TAIW) = cc( "TAIW", 234, F,234, -100, "Taiwan") Country(IC_THAI) = cc( "THAI", 235, F,235, -100, "Thailand") Country(IC_VIET) = cc( "VIET", 236, F,236, -100, "Vietnam") +Country(IC_EGYP) = cc( "EGYP", 237, F,237, -100, "Egypt") Country(IC_INTSHIPS ) = cc( "INTSHIPS" ,350 ,T, 350, -100 , "International ships, RCP6" ) end subroutine Country_Init diff --git a/DefPhotolysis_ml.f90 b/DefPhotolysis_ml.f90 old mode 100644 new mode 100755 diff --git a/Derived_ml.f90 b/Derived_ml.f90 old mode 100644 new mode 100755 index 30dfc88..e48cda2 --- a/Derived_ml.f90 +++ b/Derived_ml.f90 @@ -99,7 +99,7 @@ module Derived_ml ,FORECAST & ! only dayly (and hourly) output on FORECAST mode ,NTDAY & ! Number of 2D O3 to be saved each day (for SOMO) ! output types corresponding to instantaneous,year,month,day - ,IOU_INST, IOU_YEAR, IOU_MON, IOU_DAY, IOU_HOUR_PREVIOUS, IOU_HOUR, IOU_HOUR_MEAN + ,IOU_INST, IOU_YEAR, IOU_MON, IOU_DAY, IOU_YEAR_LASTHH, IOU_HOUR, IOU_HOUR_MEAN use MosaicOutputs_ml, only: nMosaic, MosaicOutput use OwnDataTypes_ml, only: Deriv,TXTLEN_DERIV,TXTLEN_SHORT ! type & length of names @@ -146,7 +146,7 @@ module Derived_ml ! fields we use daily outputs. For the big 3d fields, monthly output ! is sufficient. - integer, public, parameter :: LENOUT2D = IOU_HOUR_PREVIOUS ! Allows INST..DAY,H.PREV. for 2d fields + integer, public, parameter :: LENOUT2D = IOU_YEAR_LASTHH ! Allows INST..DAY,H.PREV. for 2d fields integer, public, parameter :: LENOUT3D = IOU_DAY ! Allows INST..DAY for 3d fields !will be used for: diff --git a/DryDep_ml.f90 b/DryDep_ml.f90 old mode 100644 new mode 100755 diff --git a/EQSAM_ml.f90 b/EQSAM_ml.f90 old mode 100644 new mode 100755 diff --git a/EcoSystem_ml.f90 b/EcoSystem_ml.f90 index 10e72e1..ff11fab 100644 --- a/EcoSystem_ml.f90 +++ b/EcoSystem_ml.f90 @@ -74,7 +74,7 @@ subroutine Init_EcoSystems() !Deriv(name, class, subc, txt, unit !Deriv index, f2d, dt_scale, scale, avg? Inst Yr Mn Day DepEcoSystem(iEco) = Deriv( & - name, "EcoFrac", "Area",DEF_ECOSYSTEMS(iEco) , unit, & + trim(name), "EcoFrac", "Area",trim(DEF_ECOSYSTEMS(iEco)) , trim(unit), & iEco, -99, F, 1.0, F, IOU_YEAR ) if(DEBUG_ECOSYSTEMS .and. MasterProc) & diff --git a/Emergency_ml.f90 b/Emergency_ml.f90 old mode 100644 new mode 100755 diff --git a/EmisDef_ml.f90 b/EmisDef_ml.f90 old mode 100644 new mode 100755 diff --git a/EmisGet_ml.f90 b/EmisGet_ml.f90 old mode 100644 new mode 100755 index a9c64c8..891c648 --- a/EmisGet_ml.f90 +++ b/EmisGet_ml.f90 @@ -201,33 +201,37 @@ subroutine EmisGetCdf(iem, fname,incl,excl) !/ We can include or exclude countries: if ( present(excl) ) then - if(MasterProc) print "(a,50a4)", "INCEXCTEST-E "// & - trim(code)!, (trim(excl(i)),i=1,size(excl)) + !if(MasterProc) print "(a,50a4)", "INCEXCTEST-E "// & + ! trim(code)!, (trim(excl(i)),i=1,size(excl)) if ( find_index( code, excl ) >0 ) then !if(MasterProc) print *, "INCEXCTEST-XX ", & ! trim(code), find_index( code, excl ) + !! if(MasterProc) write(*,"(a,a,i5)") "EmisGetCdf"//trim(fname)// & + !! "exludes:", trim(code), find_index( code, excl ) cycle - else - if(MasterProc) print *, "INCEXCTEST-XY ", & - trim(code), find_index( code, Country(:)%code ) end if end if if ( present(incl) ) then if ( find_index( code, incl ) <1 ) then - if(MasterProc) print "(a,i5,50a4)", "INCEXCTEST-I "//trim(code), & - find_index( code, incl ) !!, (trim(incl(i)),i=1,size(incl)) + !! if(MasterProc) write(*,"(a,a,i5)") "EmisGetCdf"//trim(fname)// & + !! "exludes:", trim(code), find_index( code, excl ) cycle end if end if + !if(MasterProc) write(*,"(3a,i5)") "EmisGetCdf"//trim(fname), & + ! "includes:", trim(code), find_index( code, excl ) + !if(MasterProc) print *, "!!EmisGetCdf"//trim(fname), & + ! "includes:", trim(code), find_index( code, excl ) + ic = find_index( code, Country(:)%code ) !from Country_ml - if(MasterProc) print "(a)", "INCEXCTEST-A "//trim(code) + if(MasterProc) print "(a)", "INCEXCTEST-A "//trim(fname)//":"//trim(code) if ( Country(ic)%code == "N/A" ) then if(MasterProc) print *, "CDFCYCLE ", ic, trim(Country(ic)%code) cycle ! see Country_ml end if - call CheckStop( ic < 1 , "CDFEMIS NegIC"//trim(code) ) + call CheckStop( ic < 1 , "CDFEMIS NegIC:"//trim(fname)//":"//trim(code) ) !if( DEBUG .and. debug_proc ) write( *,*) 'EmisGetCdf ', trim(fname), varid,trim(varname), & ! " ", trim(code), ic, isec @@ -293,7 +297,7 @@ subroutine EmisGetCdf(iem, fname,incl,excl) sumcdfemis(:,iem) = sumcdfemis(:,iem) + sumcdfemis_iem(:) do ic = 1, NLAND if ( sumcdfemis(ic,iem) > 1.0e-10 ) & - write(*,"(a,i3,f12.3)") "CDFSUM "//trim(fname), ic, sumcdfemis(ic,iem) + write(*,"(a,i5,f12.3)") "CDFSUM "//trim(fname), ic, sumcdfemis(ic,iem) end do end if diff --git a/Emissions_ml.f90 b/Emissions_ml.f90 old mode 100644 new mode 100755 index 5275937..7d80a06 --- a/Emissions_ml.f90 +++ b/Emissions_ml.f90 @@ -28,9 +28,7 @@ !_____________________________________________________________________________ ! >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> ! MOD MOD MOD MOD MOD MOD MOD MOD MOD MOD MOD MOD MOD MOD MOD MOD MOD MOD MOD - - module Emissions_ml - +module Emissions_ml ! MOD MOD MOD MOD MOD MOD MOD MOD MOD MOD MOD MOD MOD MOD MOD MOD MOD MOD MOD ! >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> ! Calls up emission read/set routines @@ -38,988 +36,848 @@ module Emissions_ml ! with the 3D model. !_____________________________________________________________________________ - use Biogenics_ml, only: SoilNOx, AnnualNdep, Ndep_trends - use CheckStop_ml,only : CheckStop,StopAll - use ChemSpecs_shl_ml, only: NSPEC_SHL - use ChemSpecs_tot_ml, only: NSPEC_TOT,NO2 - use ChemChemicals_ml, only: species - use Country_ml, only : NLAND,Country_Init,Country,IC_NAT,IC_FI,IC_NO,IC_SE - use Country_ml, only : EU27, EUMACC2 !CdfSnap - use EmisDef_ml, only : NSECTORS & ! No. sectors - ,NEMIS_FILE & ! No. emission files - ,EMIS_FILE & ! Names of species ("sox ",...) - ,NCMAX & ! Max. No. countries per grid - ,FNCMAX & ! Max. No. countries (with flat emissions) - ! per grid - ,ISNAP_DOM & ! snap index for domestic/resid emis - ,ISNAP_TRAF & ! snap index for road-traffic (SNAP7) - ,ISNAP_SHIP & ! snap index for ship emissions - ,ISNAP_NAT & ! snap index for nat. (dms) emissions - ,IQ_DMS & ! code for DMS emissions - ,NROAD_FILES & ! No. road dust emis potential files - ,ROAD_FILE & ! Names of road dust emission files - ,NROADDUST & ! No. road dust components - ,QROADDUST_FI & ! fine road dust emissions (PM2.5) - ,QROADDUST_CO & ! coarse road dust emis - ,ROADDUST_FINE_FRAC & ! fine (PM2.5) fraction of road dust emis - , ROADDUST_CLIMATE_FILE ! TEMPORARY! file for road dust climate factors - use EmisGet_ml, only : EmisGet, EmisSplit & - ,EmisGetCdf & ! cdfemis - ,GridEmis & ! cdfemis - ,nGridEmisCodes & ! cdfemis - ,GridEmisCodes & ! cdfemis - ,sumcdfemis & ! cdfemis - ,femis & ! Gets scaling factors -> e_fact - ,EmisHeights & ! Generates vertical distrib - ,nrcemis, nrcsplit, emisfrac & ! speciation routines and array - ,nemis_kprofile, emis_kprofile &! Vertical emissions profile - ,iqrc2itot & ! maps from split index to total index - ,emis_masscorr & ! 1/molwt for most species - ,emis_nsplit & ! No. species per emis file - ,RoadDustGet & - ,roaddust_masscorr ! 1/200. Hard coded at the moment, needs proper setting in EmisGet_ml... - use GridValues_ml, only: GRIDWIDTH_M & ! size of grid (m) +use Biogenics_ml, only: SoilNOx, AnnualNdep +use CheckStop_ml, only: CheckStop,StopAll +use ChemSpecs_shl_ml, only: NSPEC_SHL +use ChemSpecs_tot_ml, only: NSPEC_TOT,NO2 +use ChemChemicals_ml, only: species +use Country_ml, only: NLAND,Country_Init,Country,IC_NAT,IC_FI,IC_NO,IC_SE +use Country_ml, only: EU27, EUMACC2 !CdfSnap +use EmisDef_ml, only: & + NSECTORS & ! No. sectors + ,NEMIS_FILE & ! No. emission files + ,EMIS_FILE & ! Names of species ("sox ",...) + ,NCMAX & ! Max. No. countries per grid + ,FNCMAX & ! Max. No. countries (with flat emissions) per grid + ,ISNAP_DOM & ! snap index for domestic/resid emis + ,ISNAP_TRAF & ! snap index for road-traffic (SNAP7) + ,ISNAP_SHIP & ! snap index for ship emissions + ,ISNAP_NAT & ! snap index for nat. (dms) emissions + ,IQ_DMS & ! code for DMS emissions + ,NROAD_FILES & ! No. road dust emis potential files + ,ROAD_FILE & ! Names of road dust emission files + ,NROADDUST & ! No. road dust components + ,QROADDUST_FI & ! fine road dust emissions (PM2.5) + ,QROADDUST_CO & ! coarse road dust emis + ,ROADDUST_FINE_FRAC & ! fine (PM2.5) fraction of road dust emis + ,ROADDUST_CLIMATE_FILE ! TEMPORARY! file for road dust climate factors +use EmisGet_ml, only: & + EmisGet, EmisSplit & + ,EmisGetCdf,GridEmis,nGridEmisCodes,GridEmisCodes,sumcdfemis & ! cdfemis + ,femis & ! Gets scaling factors -> e_fact + ,EmisHeights & ! Generates vertical distrib + ,nrcemis, nrcsplit, emisfrac & ! speciation routines and array + ,nemis_kprofile, emis_kprofile &! Vertical emissions profile + ,iqrc2itot & ! maps from split index to total index + ,emis_masscorr & ! 1/molwt for most species + ,emis_nsplit & ! No. species per emis file + ,RoadDustGet & + ,roaddust_masscorr ! 1/200. Hard coded at the moment, needs proper setting in EmisGet_ml... +use GridValues_ml, only: GRIDWIDTH_M & ! size of grid (m) ,xm2 & ! map factor squared ,debug_proc,debug_li,debug_lj & - ,sigma_bnd, xmd, glat, glon,dA,dB,i_fdom,j_fdom - use Io_Nums_ml, only : IO_LOG, IO_DMS, IO_EMIS, IO_TMP - use Io_Progs_ml, only : ios, open_file, datewrite - use MetFields_ml, only : roa, ps, z_bnd, surface_precip ! ps in Pa, roa in kg/m3 - use MetFields_ml, only : t2_nwp - use ModelConstants_ml,only : KMAX_MID, KMAX_BND, PT ,dt_advec, & - EMIS_SOURCE, & ! emislist, CdfFractions - EMIS_TEST, & ! CdfSnap or none - EMIS_OUT, & ! output emissions in ASCII or not - IS_GLOBAL, & - MONTHLY_GRIDEMIS, & !NML - NBVOC, & ! > 0 if forest voc wanted - INERIS_SNAP2 , & ! INERIS/TFMM HDD20 method - DEBUG => DEBUG_EMISSIONS, MasterProc, & - DEBUG_SOILNOX , DEBUG_EMISTIMEFACS, & - DEBUG_ROADDUST , & - DEBUG_I,DEBUG_J, & - USE_DEGREEDAY_FACTORS, & - NPROC, IIFULLDOM,JJFULLDOM , & - SEAFIX_GEA_NEEDED, & ! see below - USE_LIGHTNING_EMIS,USE_AIRCRAFT_EMIS,USE_ROADDUST, & - USE_EURO_SOILNOX, USE_GLOBAL_SOILNOX ! one or the other - use NetCDF_ml, only : ReadField_CDF - use Par_ml, only : MAXLIMAX,MAXLJMAX,me,gi0,gi1,gj0,gj1, & - GIMAX, GJMAX, IRUNBEG, JRUNBEG, & - limax,ljmax, & - MSG_READ1,MSG_READ7 - use PhysicalConstants_ml, only : GRAV, AVOG - use Setup_1dfields_ml, only : rcemis ! ESX - use SmallUtils_ml, only : find_index - - use ReadField_ml, only : ReadField ! Reads ascii fields - use TimeDate_ml, only : nydays, nmdays, date, current_date, &! No. days per - daynumber,day_of_week ! year, date-type, weekday - use Timefactors_ml, only : & - NewDayFactors & ! subroutines - ,DegreeDayFactors & ! degree-days used for SNAP-2 - ,Gridded_SNAP2_Factors, gridfac_HDD & - ,fac_min,timefactors & ! subroutine - ,fac_ehh24x7 ,fac_emm, fac_edd, timefac ! time-factors - use Volcanos_ml - - - implicit none - private - - - ! subroutines: - - public :: Emissions ! Main emissions module - public :: newmonth - public :: EmisSet ! Sets emission rates every hour/time-step - public :: EmisOut ! Outputs emissions in ascii - - - ! The main code does not need to know about the following - private :: consistency_check ! Safety-checks - - INCLUDE 'mpif.h' - INTEGER STATUS(MPI_STATUS_SIZE),INFO - + ,xmd,dA,dB,i_fdom,j_fdom,glon +use Io_Nums_ml, only: IO_LOG, IO_DMS, IO_EMIS, IO_TMP +use Io_Progs_ml, only: ios, open_file, datewrite, PrintLog +use MetFields_ml, only: roa, ps, z_bnd, surface_precip ! ps in Pa, roa in kg/m3 +use MetFields_ml, only: t2_nwp ! DS_TEST SOILNO - was zero! +use ModelConstants_ml,only: & + KMAX_MID, KMAX_BND, PT ,dt_advec, & + EMIS_SOURCE, & ! emislist, CdfFractions + EMIS_TEST, & ! CdfSnap or none + EMIS_OUT, & ! output emissions in ASCII or not + IS_GLOBAL, MONTHLY_GRIDEMIS, & !NML + NBVOC, & ! > 0 if forest voc wanted + INERIS_SNAP2 , & ! INERIS/TFMM HDD20 method + DEBUG => DEBUG_EMISSIONS, MasterProc, & + DEBUG_SOILNOX, DEBUG_EMISTIMEFACS, DEBUG_ROADDUST, DEBUG_I,DEBUG_J, & + USE_DEGREEDAY_FACTORS, & + SEAFIX_GEA_NEEDED, & ! see below + USE_LIGHTNING_EMIS,USE_AIRCRAFT_EMIS,USE_ROADDUST, & + USE_EURO_SOILNOX, USE_GLOBAL_SOILNOX, EURO_SOILNOX_DEPSCALE! one or the other +use NetCDF_ml, only: ReadField_CDF +use Par_ml, only: MAXLIMAX,MAXLJMAX, GIMAX,GJMAX, IRUNBEG,JRUNBEG,& + me,limax,ljmax, MSG_READ1,MSG_READ7 +use PhysicalConstants_ml,only: GRAV, AVOG +use Setup_1dfields_ml,only: rcemis ! ESX +use SmallUtils_ml, only: find_index +use ReadField_ml, only: ReadField ! Reads ascii fields +use TimeDate_ml, only: nydays, nmdays, date, current_date, &! No. days per + daynumber,day_of_week ! year, date-type, weekday +use Timefactors_ml, only: & + NewDayFactors & ! subroutines + ,DegreeDayFactors & ! degree-days used for SNAP-2 + ,Gridded_SNAP2_Factors, gridfac_HDD & + ,fac_min,timefactors & ! subroutine + ,fac_ehh24x7 ,fac_emm, fac_edd, timefac ! time-factors +use Volcanos_ml + +implicit none +private + +! subroutines: +public :: Emissions ! Main emissions module +public :: newmonth +public :: EmisSet ! Sets emission rates every hour/time-step +public :: EmisOut ! Outputs emissions in ascii + +! The main code does not need to know about the following +private :: consistency_check ! Safety-checks + +INCLUDE 'mpif.h' +INTEGER STATUS(MPI_STATUS_SIZE),INFO - ! land-code information in each grid square - needed to know which country - ! is emitting. - ! nlandcode = No. countries in grid square - ! landcode = Country codes for that grid square - integer, private, save, allocatable, dimension(:,:) :: nlandcode - integer, private, save, allocatable, dimension(:,:,:) :: landcode - ! for flat emissions, i.e. no vertical extent: - integer, private, save, allocatable, dimension(:,:) :: flat_nlandcode - integer, private, save, allocatable, dimension(:,:,:):: flat_landcode - ! for road dust emission potentials: - integer, private, save, allocatable, dimension(:,:) :: road_nlandcode - integer, private, save, allocatable, dimension(:,:,:) :: road_landcode - - ! - ! The output emission matrix for the 11-SNAP data is snapemis: - ! - - real, private, allocatable, dimension(:,:,:,:,:) & - , save :: snapemis ! main emission arrays, in kg/m2/s +! land-code information in each grid square - needed to know which country +! is emitting. +! nlandcode = No. countries in grid square +! landcode = Country codes for that grid square +integer, private, save, allocatable, dimension(:,:) :: nlandcode +integer, private, save, allocatable, dimension(:,:,:) :: landcode +! for flat emissions, i.e. no vertical extent: +integer, private, save, allocatable, dimension(:,:) :: flat_nlandcode +integer, private, save, allocatable, dimension(:,:,:) :: flat_landcode +! for road dust emission potentials: +integer, private, save, allocatable, dimension(:,:) :: road_nlandcode +integer, private, save, allocatable, dimension(:,:,:) :: road_landcode - real, private, allocatable, dimension(:,:,:,:) & - , save :: snapemis_flat ! main emission arrays, in kg/m2/s - - real, private, allocatable, dimension(:,:,:,:) & ! Not sure if it is really necessary to keep the country info; gives rather messy code but consistent with the rest at least (and can do the seasonal scaling for Nordic countries in the code instead of as preprocessing) - , save :: roaddust_emis_pot ! main road dust emission potential arrays, in kg/m2/s (to be scaled!) - - ! We store the emissions for output to d_2d files and netcdf in kg/m2/s - - real, public, allocatable, dimension(:,:,:), save :: SumSnapEmis +! +! The output emission matrix for the 11-SNAP data is snapemis: +! +real, private, allocatable, dimension(:,:,:,:,:), save :: & + snapemis ! main emission arrays, in kg/m2/s - logical, save, private :: first_dms_read +real, private, allocatable, dimension(:,:,:,:), save :: & + snapemis_flat ! main emission arrays, in kg/m2/s - ! Emissions for input to chemistry routines +real, private, allocatable, dimension(:,:,:,:), save :: & +! Not sure if it is really necessary to keep the country info; gives rather messy code but consistent with the rest at least (and can do the seasonal scaling for Nordic countries in the code instead of as preprocessing) + roaddust_emis_pot ! main road dust emission potential arrays, in kg/m2/s (to be scaled!) - ! KEMISTOP added to avoid hard-coded KMAX_MID-3: +! We store the emissions for output to d_2d files and netcdf in kg/m2/s +real, public, allocatable, dimension(:,:,:), save :: SumSnapEmis - !integer, public, parameter :: KEMISTOP = KMAX_MID - NEMISLAYERS + 1 - integer, public, save :: KEMISTOP ! not defined yet= KMAX_MID - nemis_kprofile + 1 - real, public, allocatable, save, dimension(:,:,:,:) :: & - gridrcemis & ! varies every time-step (as ps changes) - ,gridrcemis0 ! varies every hour - real, public, allocatable, save, dimension(:,:,:) :: & - gridrcroadd & ! Road dust emissions - ,gridrcroadd0 ! varies every hour +logical, save, private :: first_dms_read - ! and for budgets (not yet used - not changed dimension) - real, public, save, dimension(NSPEC_SHL+1:NSPEC_TOT) :: totemadd +! Emissions for input to chemistry routines +! KEMISTOP added to avoid hard-coded KMAX_MID-3: +!integer, public, parameter :: KEMISTOP = KMAX_MID - NEMISLAYERS + 1 +integer, public, save :: KEMISTOP ! not defined yet= KMAX_MID - nemis_kprofile + 1 +real, public, allocatable, save, dimension(:,:,:,:) :: & + gridrcemis, & ! varies every time-step (as ps changes) + gridrcemis0 ! varies every hour +real, public, allocatable, save, dimension(:,:,:) :: & + gridrcroadd, & ! Road dust emissions + gridrcroadd0 ! varies every hour - integer, private, save :: iemCO ! index of CO emissions, for debug +! and for budgets (not yet used - not changed dimension) +real, public, save, dimension(NSPEC_SHL+1:NSPEC_TOT) :: totemadd +integer, private, save :: iemCO ! index of CO emissions, for debug - logical, parameter ::USE_OLDSCHEME_ROADDUST=.false. !temporary until the new scheme is validated +logical, parameter :: USE_OLDSCHEME_ROADDUST=.false. !temporary until the new scheme is validated contains - !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> - subroutine Emissions(year) - - - ! Calls main emission reading routines - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !*********************************************************************** - ! DESCRIPTION: - ! 0) Call set_molwts and set_emisconv_and_iq, followed by - ! consistency check - ! 1) Call some setups: - ! Country_Init - ! timefactors: monthly and daily factors, + time zone - ! -> fac_emm, fac_edd arrays, timezone - ! 2) Read in emission correction file femis - ! 3) Call emis_split for speciations - ! 4) Read the annual emission totals in each grid-square, converts - ! units and corrects using femis data. - ! - ! The output emission matrix for the 11-SNAP data is snapemis: - ! - ! real snapemis (NSECTORS,MAXLIMAX,MAXLJMAX,NCMAX,NEMIS_FILES) - ! - !********************************************************************** - - !--arguments - integer, intent(in) :: year ! Year ( 4-digit) - - !-- local variables - real :: conv ! Conversion factor - integer :: i, j ! Loop variables - real :: tonne_to_kgm2s ! Converts tonnes/grid to kg/m2/s - real :: ccsum ! Sum of emissions for one country - - ! arrays for whole EMEP area: - ! additional arrays on host only for landcode, nlandcode - ! BIG arrays ... will be needed only on me=0. Make allocatable - ! to reduce static memory requirements. - - real, allocatable, dimension(:,:,:,:) :: globemis - integer, allocatable, dimension(:,:) :: globnland - integer, allocatable, dimension(:,:,:) :: globland - real, allocatable, dimension(:,:,:) :: globemis_flat - integer, allocatable, dimension(:,:) :: flat_globnland - integer, allocatable, dimension(:,:,:) :: flat_globland - real, allocatable, dimension(:,:,:) :: globroad_dust_pot ! Road dust emission potentials - integer, allocatable, dimension(:,:) :: road_globnland - integer, allocatable, dimension(:,:,:) :: road_globland - real, allocatable, dimension(:,:) :: RoadDustEmis_climate_factor ! Climatic factor for scaling road dust emissions (in TNO model based on yearly average soil water) - integer :: err1, err2, err3, err4, err5, err6, err7, err8, err9 ! Error messages - integer :: fic ,insec,inland,iemis - integer :: iic,ic,n ! country codes - integer :: isec ! loop variables: emission sectors - integer :: iem,iem2 ! loop variable over pollutants (1..NEMIS_FILE) - integer :: icc,ncc ! loop variables over sources - character(len=40) :: fname ! File name - character(len=300) :: inputline - real :: tmpclimfactor - - ! Emission sums (after e_fact adjustments): - real, dimension(NEMIS_FILE) :: emsum ! Sum emis over all countries - real, dimension(NLAND,NEMIS_FILE) :: sumemis, sumemis_local ! Sum of emissions per country - real, dimension(NEMIS_FILE) :: sumEU ! Sum of emissions in EU - - ! Road dust emission potential sums (just for testing the code, the actual emissions are weather dependent!) - real, dimension(NROAD_FILES) :: roaddustsum ! Sum emission potential over all countries - real, dimension(NLAND,NROAD_FILES) :: sumroaddust ! Sum of emission potentials per country - real, dimension(NLAND,NROAD_FILES) :: sumroaddust_local ! Sum of emission potentials per country in subdomain - real :: fractions(MAXLIMAX,MAXLJMAX,NCMAX),SMI(MAXLIMAX,MAXLJMAX),SMI_roadfactor - logical ::SMI_defined=.false. - logical :: my_first_call=.true. ! Used for femis call - real, allocatable ::emis_tot(:,:) - character(len=40) :: varname - - if (MasterProc) write(6,*) "Reading emissions for year", year - - ! 0) set molwts, conversion factors (e.g. tonne NO2 -> tonne N), and - ! emission indices (IQSO2=.., ) - - allocate(nlandcode(MAXLIMAX,MAXLJMAX),landcode(MAXLIMAX,MAXLJMAX,NCMAX)) - nlandcode=0 - landcode=0 - allocate(flat_nlandcode(MAXLIMAX,MAXLJMAX),flat_landcode(MAXLIMAX,MAXLJMAX,FNCMAX)) - flat_nlandcode=0 - flat_landcode=0 - allocate(road_nlandcode(MAXLIMAX,MAXLJMAX),road_landcode(MAXLIMAX,MAXLJMAX,NCMAX)) - road_nlandcode=0 - road_landcode=0 - allocate(snapemis(NSECTORS,MAXLIMAX,MAXLJMAX,NCMAX,NEMIS_FILE)) - snapemis=0.0 - allocate(snapemis_flat(MAXLIMAX,MAXLJMAX,FNCMAX,NEMIS_FILE)) - snapemis_flat=0.0 - allocate(roaddust_emis_pot(MAXLIMAX,MAXLJMAX,NCMAX,NROAD_FILES)) - roaddust_emis_pot=0.0 - allocate(SumSnapEmis(MAXLIMAX,MAXLJMAX,NEMIS_FILE)) - SumSnapEmis=0.0 - - !========================= - call Country_Init() ! In Country_ml, => NLAND, country codes and - ! names, timezone - ! The GEA emission data, which is used for EUCAARI runs on the HIRHAM - ! domain have in several sea grid cells non-zero emissions in other sectors - ! than SNAP8 and there are also NH3 emission over sea areas. The former - ! problem makes the code crash if the sea areas are defined as - ! sea (sea=T), so we treat them as land in the EUCAARI/HIRHAM runs - ! (sea=F). This is a problem with GEA emission data only, not the - ! HIRHAM domain! When e.g. interpolated EMEP emissions are used on - ! the HIRHAM domain, this is not a problem. - - if ( SEAFIX_GEA_NEEDED ) then ! Special fix for HIRHAM/GEA - - ! Could have hard-coded 30-34, but best to avoid: - ic= find_index("BAS", Country(:)%code ) - Country( ic )%is_sea = .false. - ic= find_index("NOS", Country(:)%code ) - Country( ic )%is_sea = .false. - ic= find_index("ATL", Country(:)%code ) - Country( ic )%is_sea = .false. - ic= find_index("MED", Country(:)%code ) - Country( ic )%is_sea = .false. - ic= find_index("BLS", Country(:)%code ) - Country( ic )%is_sea = .false. - - end if ! HIRHAM/GEA fix - - !========================= - - call consistency_check() ! Below - !========================= - - ios = 0 - - if ( USE_DEGREEDAY_FACTORS ) & - call DegreeDayFactors(0) ! See if we have gridded SNAP-2 - - call EmisHeights() ! vertical emissions profile - KEMISTOP = KMAX_MID - nemis_kprofile + 1 +!*********************************************************************** +subroutine Emissions(year) +! Calls main emission reading routines +!----------------------------------------------------------------------! +! DESCRIPTION: +! 0) Call set_molwts and set_emisconv_and_iq, followed by +! consistency check +! 1) Call some setups: +! Country_Init +! timefactors: monthly and daily factors, + time zone +! -> fac_emm, fac_edd arrays, timezone +! 2) Read in emission correction file femis +! 3) Call emis_split for speciations +! 4) Read the annual emission totals in each grid-square, converts +! units and corrects using femis data. +! +! The output emission matrix for the 11-SNAP data is snapemis: +! +! real snapemis (NSECTORS,MAXLIMAX,MAXLJMAX,NCMAX,NEMIS_FILES) +! +!----------------------------------------------------------------------! + !--arguments + integer, intent(in) :: year ! Year ( 4-digit) + + !-- local variables + real :: conv ! Conversion factor + integer :: i, j ! Loop variables + real :: tonne_to_kgm2s ! Converts tonnes/grid to kg/m2/s + real :: ccsum ! Sum of emissions for one country + + ! arrays for whole EMEP area: + ! additional arrays on host only for landcode, nlandcode + ! BIG arrays ... will be needed only on me=0. Make allocatable + ! to reduce static memory requirements. + real, allocatable, dimension(:,:,:,:) :: globemis + integer, allocatable, dimension(:,:) :: globnland + integer, allocatable, dimension(:,:,:) :: globland + real, allocatable, dimension(:,:,:) :: globemis_flat + integer, allocatable, dimension(:,:) :: flat_globnland + integer, allocatable, dimension(:,:,:) :: flat_globland + real, allocatable, dimension(:,:,:) :: globroad_dust_pot ! Road dust emission potentials + integer, allocatable, dimension(:,:) :: road_globnland + integer, allocatable, dimension(:,:,:) :: road_globland + real, allocatable, dimension(:,:) :: RoadDustEmis_climate_factor ! Climatic factor for scaling road dust emissions (in TNO model based on yearly average soil water) + integer :: err1, err2, err3, err4, err5, err6, err7, err8, err9 ! Error messages + integer :: fic ,insec,inland,iemis + integer :: iic,ic,n ! country codes + integer :: isec ! loop variables: emission sectors + integer :: iem,iem2 ! loop variable over pollutants (1..NEMIS_FILE) + integer :: icc,ncc ! loop variables over sources + character(len=40) :: fname ! File name + character(len=300) :: inputline + real :: tmpclimfactor + + ! Emission sums (after e_fact adjustments): + real, dimension(NEMIS_FILE) :: emsum ! Sum emis over all countries + real, dimension(NLAND,NEMIS_FILE) :: sumemis, sumemis_local ! Sum of emissions per country + real, dimension(NEMIS_FILE) :: sumEU ! Sum of emissions in EU + + ! Road dust emission potential sums (just for testing the code, the actual emissions are weather dependent!) + real, dimension(NROAD_FILES) :: roaddustsum ! Sum emission potential over all countries + real, dimension(NLAND,NROAD_FILES) :: sumroaddust ! Sum of emission potentials per country + real, dimension(NLAND,NROAD_FILES) :: sumroaddust_local ! Sum of emission potentials per country in subdomain + real :: fractions(MAXLIMAX,MAXLJMAX,NCMAX),SMI(MAXLIMAX,MAXLJMAX),SMI_roadfactor + logical ::SMI_defined=.false. + logical :: my_first_call=.true. ! Used for femis call + real, allocatable ::emis_tot(:,:) + character(len=40) :: varname + + if (MasterProc) write(6,*) "Reading emissions for year", year + + ! 0) set molwts, conversion factors (e.g. tonne NO2 -> tonne N), and + ! emission indices (IQSO2=.., ) + + allocate(nlandcode(MAXLIMAX,MAXLJMAX),landcode(MAXLIMAX,MAXLJMAX,NCMAX)) + nlandcode=0 + landcode=0 + allocate(flat_nlandcode(MAXLIMAX,MAXLJMAX),flat_landcode(MAXLIMAX,MAXLJMAX,FNCMAX)) + flat_nlandcode=0 + flat_landcode=0 + allocate(road_nlandcode(MAXLIMAX,MAXLJMAX),road_landcode(MAXLIMAX,MAXLJMAX,NCMAX)) + road_nlandcode=0 + road_landcode=0 + allocate(snapemis(NSECTORS,MAXLIMAX,MAXLJMAX,NCMAX,NEMIS_FILE)) + snapemis=0.0 + allocate(snapemis_flat(MAXLIMAX,MAXLJMAX,FNCMAX,NEMIS_FILE)) + snapemis_flat=0.0 + allocate(roaddust_emis_pot(MAXLIMAX,MAXLJMAX,NCMAX,NROAD_FILES)) + roaddust_emis_pot=0.0 + allocate(SumSnapEmis(MAXLIMAX,MAXLJMAX,NEMIS_FILE)) + SumSnapEmis=0.0 + + !========================= + call Country_Init() ! In Country_ml, => NLAND, country codes and names, timezone + + ! The GEA emission data, which is used for EUCAARI runs on the HIRHAM + ! domain have in several sea grid cells non-zero emissions in other sectors + ! than SNAP8 and there are also NH3 emission over sea areas. The former + ! problem makes the code crash if the sea areas are defined as + ! sea (sea=T), so we treat them as land in the EUCAARI/HIRHAM runs + ! (sea=F). This is a problem with GEA emission data only, not the + ! HIRHAM domain! When e.g. interpolated EMEP emissions are used on + ! the HIRHAM domain, this is not a problem. + if(SEAFIX_GEA_NEEDED) then ! Special fix for HIRHAM/GEA + ! Could have hard-coded 30-34, but best to avoid: + do n=1,5 + select case(n) + case(1);ic=find_index("BAS",Country(:)%code) + case(2);ic=find_index("NOS",Country(:)%code) + case(3);ic=find_index("ATL",Country(:)%code) + case(4);ic=find_index("MED",Country(:)%code) + case(5);ic=find_index("BLS",Country(:)%code) + endselect + call CheckStop(ic<1,"Country_Init error in HIRHAM/GEA fix") + Country(ic)%is_sea = .false. + enddo + endif ! HIRHAM/GEA fix - if( MasterProc) then !::::::: ALL READ-INS DONE IN HOST PROCESSOR :::: + !========================= + call consistency_check() ! Below + !========================= + ios = 0 - write(6,*) "Reading monthly and daily timefactors" - !========================= - call timefactors(year) ! => fac_emm, fac_edd - !========================= + if(USE_DEGREEDAY_FACTORS) call DegreeDayFactors(0)! See if we have gridded SNAP-2 - endif + call EmisHeights() ! vertical emissions profile + KEMISTOP = KMAX_MID - nemis_kprofile + 1 + if(MasterProc) then !::::::: ALL READ-INS DONE IN HOST PROCESSOR :::: + write(*,*) "Reading monthly and daily timefactors" !========================= - call EmisSplit() ! In EmisGet_ml, => emisfrac + call timefactors(year) ! => fac_emm, fac_edd !========================= - call CheckStop(ios, "ioserror: EmisSplit") - - - ! #################################### - ! Broadcast monthly and Daily factors (and hourly factors if needed/wanted) - CALL MPI_BCAST( fac_emm ,8*NLAND*12*NSECTORS*NEMIS_FILE,MPI_BYTE, 0,MPI_COMM_WORLD,INFO) - CALL MPI_BCAST( fac_edd ,8*NLAND*7*NSECTORS*NEMIS_FILE,MPI_BYTE, 0,MPI_COMM_WORLD,INFO) - CALL MPI_BCAST( fac_ehh24x7 ,8*NSECTORS*24*7,MPI_BYTE, 0,MPI_COMM_WORLD,INFO) - - !define fac_min for all processors - do iemis = 1, NEMIS_FILE - do insec = 1, NSECTORS - do inland = 1, NLAND - fac_min(inland,insec,iemis) = minval( fac_emm(inland,:,insec,iemis) ) - enddo - enddo - enddo - if(INERIS_SNAP2 )THEN ! INERIS do not use any base-line for SNAP2 - fac_min(:,ISNAP_DOM,:) = 0. - end if - - !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: - ! c4b) Set up DMS factors here - to be used in newmonth - ! Taken from IQ_DMS=35 for SO2 nature (sector 11) - ! first_dms_read is true until first call to newmonth finished. - - - first_dms_read = .true. - - ! 4) Read emission files - - ! allocate for me=0 only: - - err1 = 0 - if ( MasterProc ) then - - if (DEBUG) write(unit=6,fmt=*) "TTT me ", me , "pre-allocate" - allocate(globnland(GIMAX,GJMAX),stat=err1) - allocate(globland(GIMAX,GJMAX,NCMAX),stat=err2) - allocate(globemis(NSECTORS,GIMAX,GJMAX,NCMAX),stat=err3) - - allocate(flat_globnland(GIMAX,GJMAX),stat=err4) - allocate(flat_globland(GIMAX,GJMAX,FNCMAX),stat=err5) - allocate(globemis_flat(GIMAX,GJMAX,FNCMAX),stat=err6) - - call CheckStop(err1, "Allocation error 1 - globland") - call CheckStop(err2, "Allocation error 2 - globland") - call CheckStop(err3, "Allocation error 3 - globland") - call CheckStop(err4, "Allocation error 4 - globland") - call CheckStop(err5, "Allocation error 5 - globland") - call CheckStop(err6, "Allocation error 6 - globland") - - if(USE_ROADDUST)then - allocate(road_globnland(GIMAX,GJMAX),stat=err7) - allocate(road_globland(GIMAX,GJMAX,NCMAX),stat=err8) - allocate(globroad_dust_pot(GIMAX,GJMAX,NCMAX),stat=err9) - allocate(RoadDustEmis_climate_factor(GIMAX,GJMAX),stat=err1) - - call CheckStop(err7, "Allocation error 7 - globroadland") - call CheckStop(err8, "Allocation error 8 - globroadland") - call CheckStop(err9, "Allocation error 9 - globroad_dust_pot") - call CheckStop(err1, "Allocation error 1 - RoadDustEmis_climate_factor") - endif ! road dust - - ! Initialise with 0 - globnland(:,:) = 0 - flat_globnland(:,:)=0 - globland(:,:,:) = 0 - globemis(:,:,:,:) = 0 - flat_globland(:,:,:)=0 - globemis_flat(:,:,:) =0 - - if(USE_ROADDUST)then - road_globnland(:,:)=0 - road_globland(:,:,:)=0 - globroad_dust_pot(:,:,:)=0. - RoadDustEmis_climate_factor(:,:)=1.0 ! default, no scaling - endif ! road dust - - end if - - ! Get emission scaling factors - !>============================ - - if ( my_first_call ) then - sumemis(:,:) = 0.0 ! initialize sums - ios = 0 - call femis() ! emission factors (femis.dat file) - if ( ios /= 0 )return - my_first_call = .false. - endif - !>============================ - - - do iem = 1, NEMIS_FILE - - !TESTING NEW SYSTEM, no effect on results so far - ! NB SLOW! - if( EMIS_TEST == "CdfSnap" ) then - !if ( iem < 3) then - - ! ************************* - fname="GriddedSnapEmis_" // trim(EMIS_FILE(iem)) // ".nc" - call EmisGetCdf(iem,fname,incl=EUMACC2) - ! ************************* - fname="GlobalSnapEmis_" // trim(EMIS_FILE(iem)) // ".nc" - call EmisGetCdf(iem,fname,excl=EUMACC2) - ! ************************* - + endif + !========================= + call EmisSplit() ! In EmisGet_ml, => emisfrac + !========================= + call CheckStop(ios, "ioserror: EmisSplit") + ! #################################### + ! Broadcast monthly and Daily factors (and hourly factors if needed/wanted) + CALL MPI_BCAST(fac_emm,8*NLAND*12*NSECTORS*NEMIS_FILE,MPI_BYTE,0,MPI_COMM_WORLD,INFO) + CALL MPI_BCAST(fac_edd,8*NLAND*7*NSECTORS*NEMIS_FILE,MPI_BYTE,0,MPI_COMM_WORLD,INFO) + CALL MPI_BCAST(fac_ehh24x7,8*NSECTORS*24*7,MPI_BYTE,0,MPI_COMM_WORLD,INFO) + + !define fac_min for all processors + forall(iemis=1:NEMIS_FILE,insec=1:NSECTORS,inland=1:NLAND) & + fac_min(inland,insec,iemis) = minval(fac_emm(inland,:,insec,iemis)) + if(INERIS_SNAP2) & ! INERIS do not use any base-line for SNAP2 + fac_min(:,ISNAP_DOM,:) = 0. + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + ! c4b) Set up DMS factors here - to be used in newmonth + ! Taken from IQ_DMS=35 for SO2 nature (sector 11) + ! first_dms_read is true until first call to newmonth finished. + first_dms_read = .true. + + ! 4) Read emission files + + ! allocate for me=0 only: + err1 = 0 + if(MasterProc) then + if(DEBUG) write(*,*) "TTT me ", me , "pre-allocate" + allocate(globnland(GIMAX,GJMAX),stat=err1) + allocate(globland(GIMAX,GJMAX,NCMAX),stat=err2) + allocate(globemis(NSECTORS,GIMAX,GJMAX,NCMAX),stat=err3) + + allocate(flat_globnland(GIMAX,GJMAX),stat=err4) + allocate(flat_globland(GIMAX,GJMAX,FNCMAX),stat=err5) + allocate(globemis_flat(GIMAX,GJMAX,FNCMAX),stat=err6) + + call CheckStop(err1, "Allocation error 1 - globland") + call CheckStop(err2, "Allocation error 2 - globland") + call CheckStop(err3, "Allocation error 3 - globland") + call CheckStop(err4, "Allocation error 4 - globland") + call CheckStop(err5, "Allocation error 5 - globland") + call CheckStop(err6, "Allocation error 6 - globland") + + if(USE_ROADDUST)then + allocate(road_globnland(GIMAX,GJMAX),stat=err7) + allocate(road_globland(GIMAX,GJMAX,NCMAX),stat=err8) + allocate(globroad_dust_pot(GIMAX,GJMAX,NCMAX),stat=err9) + allocate(RoadDustEmis_climate_factor(GIMAX,GJMAX),stat=err1) + + call CheckStop(err7, "Allocation error 7 - globroadland") + call CheckStop(err8, "Allocation error 8 - globroadland") + call CheckStop(err9, "Allocation error 9 - globroad_dust_pot") + call CheckStop(err1, "Allocation error 1 - RoadDustEmis_climate_factor") + endif ! road dust + + ! Initialise with 0 + globnland(:,:) = 0 + flat_globnland(:,:)=0 + globland(:,:,:) = 0 + globemis(:,:,:,:) = 0 + flat_globland(:,:,:)=0 + globemis_flat(:,:,:) =0 + + if(USE_ROADDUST)then + road_globnland(:,:)=0 + road_globland(:,:,:)=0 + globroad_dust_pot(:,:,:)=0. + RoadDustEmis_climate_factor(:,:)=1.0 ! default, no scaling + endif ! road dust + endif + + ! Get emission scaling factors + !>============================ + if(my_first_call) then + sumemis(:,:) = 0.0 ! initialize sums + ios = 0 + call femis() ! emission factors (femis.dat file) + if(ios/=0) return + my_first_call = .false. + endif + !>============================ + + do iem = 1, NEMIS_FILE + !TESTING NEW SYSTEM, no effect on results so far + ! NB SLOW! + if(EMIS_TEST=="CdfSnap") then + !if ( iem < 3) then + ! ************************* + fname="GriddedSnapEmis_"//trim(EMIS_FILE(iem))//".nc" + call EmisGetCdf(iem,fname,incl=EUMACC2) + ! ************************* + fname="GlobalSnapEmis_"//trim(EMIS_FILE(iem))//".nc" + call EmisGetCdf(iem,fname,excl=EUMACC2) + ! ************************* - if( debug_proc ) then - ncc = nGridEmisCodes(debug_li, debug_lj) - do isec = 1, 11 - write(*,"(a,2i4,4i3,2i4,9f12.4)")"CDF out "// & - trim(EMIS_FILE(iem)), i_fdom(debug_li), j_fdom(debug_lj), & - (GridEmisCodes(debug_li, debug_lj,icc),icc=1,4), & - ncc, isec, (GridEmis(isec,debug_li, debug_lj,i,iem), i=1,4) - end do - end if ! debug - !end if ! iem - end if ! EMIS_TEST - - ! now again test for me=0 - if(EMIS_SOURCE == "emislist" )then - if ( MasterProc ) then - - ! read in global emissions for one pollutant - ! ***************** - call EmisGet(iem,EMIS_FILE(iem),IRUNBEG,JRUNBEG,GIMAX,GJMAX, & - globemis,globnland,globland,sumemis,& - globemis_flat,flat_globnland,flat_globland) - ! ***************** - - - emsum(iem) = sum( globemis(:,:,:,:) ) + & - sum( globemis_flat(:,:,:) ) ! hf - endif ! MasterProc - - call CheckStop(ios, "ios error: EmisGet") - - ! Send data to processors - ! as e.g. snapemis (NSECTORS,MAXLIMAX,MAXLJMAX,NCMAX,iem) - ! send to nodes + if(debug_proc) then + ncc = nGridEmisCodes(debug_li, debug_lj) + do isec = 1, 11 + write(*,"(a,2i4,4i3,2i4,9f12.4)")"CDF out "// & + trim(EMIS_FILE(iem)), i_fdom(debug_li), j_fdom(debug_lj), & + (GridEmisCodes(debug_li, debug_lj,icc),icc=1,4), & + ncc, isec, (GridEmis(isec,debug_li, debug_lj,i,iem), i=1,4) + enddo + endif ! debug + !end if ! iem + endif ! EMIS_TEST + + ! now again test for me=0 + select case(EMIS_SOURCE) + case("emislist") + if(MasterProc) then + ! read in global emissions for one pollutant + ! ***************** + call EmisGet(iem,EMIS_FILE(iem),IRUNBEG,JRUNBEG,GIMAX,GJMAX, & + globemis,globnland,globland,sumemis,& + globemis_flat,flat_globnland,flat_globland) + ! ***************** + + emsum(iem) = sum(globemis(:,:,:,:)) + sum(globemis_flat(:,:,:)) ! hf + endif ! MasterProc + call CheckStop(ios, "ios error: EmisGet") - call global2local(globemis,snapemis(1,1,1,1,iem),MSG_READ1, & - NSECTORS,GIMAX,GJMAX,NCMAX,1,1) - - call global2local(globemis_flat,snapemis_flat(1,1,1,iem),MSG_READ1, & - 1,GIMAX,GJMAX,FNCMAX,1,1) - - ! real snapemis (NSECTORS,MAXLIMAX,MAXLJMAX,NCMAX,NEMIS_FILES) - ! GridEmis (NSECTORS,MAXLIMAX,MAXLJMAX,NCMAX,NEMIS_FILE) - - - if ( EMIS_TEST == "CdfSnap" .and. debug_proc ) then - i=debug_li; j=debug_lj - do isec = 1, 11 - if( snapemis(isec,i,j,1,iem)+GridEmis(isec,i,j,1,iem) > 1.0e-10) & - write(*,"(a,i3,2es10.3)") "CDFCOMP", isec, & - snapemis(isec,i,j,1,iem) ,GridEmis(isec,i,j,1,iem) - end do - end if - else if ( EMIS_SOURCE == "CdfFractions" ) then - !use grid independent netcdf emission file - !experimental so far. Needs a lot of reorganization - if(.not.allocated(emis_tot))then - allocate(emis_tot(MAXLIMAX,MAXLJMAX)) - endif -77 format(A,I2.2) - sumemis_local(:,iem)=0.0 - do isec=1,NSECTORS - - if(iem==1)varname='SOx_sec' - if(iem==2)varname='NOx_sec' - if(iem==3)varname='CO_sec' - if(iem==4)varname='NMVOC_sec' - if(iem==5)varname='NH3_sec' - if(iem==6)varname='PM25_sec' - if(iem==7)varname='PMco_sec' - write(varname,77)trim(varname),isec - call ReadField_CDF('/global/work/mifapw/emep/Data/Emis_TNO7.nc',varname,emis_tot(1,1),& - nstart=1,interpol='mass_conservative', & - fractions_out=fractions,CC_out=landcode,Ncc_out=nlandcode,needed=.true.,debug_flag=.true.,Undef=0.0) - - if(debug_proc) write(*,*) "CDFTNO ", me, iem, isec, trim(varname) - do j=1,ljmax - do i=1,limax - if(nlandcode(i,j)>NCMAX)then - write(*,*)"To many emitter countries in one gridcell: ", me,i,j,nlandcode(i,j) - call StopAll("To many countries in one gridcell ") - endif - do n=1,nlandcode(i,j) - snapemis(isec,i,j,n,iem)=snapemis(isec,i,j,n,iem)+fractions(i,j,n)*emis_tot(i,j) - ic=1 - if(landcode(i,j,n)<=NLAND)ic=landcode(i,j,n)!most country_index are equal to country_code + ! Send data to processors + ! as e.g. snapemis (NSECTORS,MAXLIMAX,MAXLJMAX,NCMAX,iem) + ! send to nodes + call global2local(globemis,snapemis(1,1,1,1,iem),MSG_READ1,& + NSECTORS,GIMAX,GJMAX,NCMAX,1,1) + call global2local(globemis_flat,snapemis_flat(1,1,1,iem),MSG_READ1,& + 1,GIMAX,GJMAX,FNCMAX,1,1) + +! real snapemis (NSECTORS,MAXLIMAX,MAXLJMAX,NCMAX,NEMIS_FILES) +! GridEmis (NSECTORS,MAXLIMAX,MAXLJMAX,NCMAX,NEMIS_FILE) + if(EMIS_TEST=="CdfSnap".and.debug_proc) then + i=debug_li;j=debug_lj + do isec = 1, 11 + if(snapemis(isec,i,j,1,iem)+GridEmis(isec,i,j,1,iem)>1.0e-10) & + write(*,"(a,i3,2es10.3)") "CDFCOMP", isec, & + snapemis(isec,i,j,1,iem) ,GridEmis(isec,i,j,1,iem) + enddo + endif + + case("CdfFractions") + !use grid independent netcdf emission file + !experimental so far. Needs a lot of reorganization + if(.not.allocated(emis_tot)) allocate(emis_tot(MAXLIMAX,MAXLJMAX)) + sumemis_local(:,iem)=0.0 + + do isec=1,NSECTORS + select case(iem) + case(1);varname='SOx_sec' + case(2);varname='NOx_sec' + case(3);varname='CO_sec' + case(4);varname='NMVOC_sec' + case(5);varname='NH3_sec' + case(6);varname='PM25_sec' + case(7);varname='PMco_sec' + endselect + write(varname,"(A,I2.2)")trim(varname),isec + call ReadField_CDF('EmisFracs_TNO7.nc',varname,emis_tot(1,1),nstart=1,& + interpol='mass_conservative',fractions_out=fractions,& + CC_out=landcode,Ncc_out=nlandcode,needed=.true.,debug_flag=.true.,& + Undef=0.0) + + if(debug_proc) write(*,*) "CDFTNO ",me,iem,isec,trim(varname) + do j=1,ljmax + do i=1,limax + if(nlandcode(i,j)>NCMAX)then + write(*,*)"To many emitter countries in one gridcell: ",& + me,i,j,nlandcode(i,j) + call StopAll("To many countries in one gridcell ") + endif + do n=1,nlandcode(i,j) + snapemis(isec,i,j,n,iem)=snapemis(isec,i,j,n,iem)& + +fractions(i,j,n)*emis_tot(i,j) + ic=1 + if(landcode(i,j,n)<=NLAND) & + ic=landcode(i,j,n) ! most country_index are equal to country_code !TESTE - if ( debug_proc .and. i==debug_li .and. j==debug_lj ) then ! .and. iem == iemCO ) then - write(*,"(a,2i3,4es10.3)") "FracCdf:" // trim(EMIS_FILE(iem)), & - isec, n, emis_tot(i,j), fractions(i,j,n), fractions(i,j,n)*emis_tot(i,j) - end if + if(debug_proc.and.i==debug_li.and.j==debug_lj) & ! .and. iem == iemCO ) then + write(*,"(a,2i3,4es10.3)") "FracCdf: "//trim(EMIS_FILE(iem)), & + isec,n,emis_tot(i,j),fractions(i,j,n),fractions(i,j,n)*emis_tot(i,j) !TESTE - - if(Country(ic)%index/=landcode(i,j,n))then - !if not equal, find which index correspond to country_code - do ic=1,NLAND - if((Country(ic)%index==landcode(i,j,n)))exit - enddo - if(ic>NLAND)then - write(*,*)"COUNTRY CODE NOT RECOGNIZED OR UNDEFINED: ", landcode(i,j,n) - call StopAll("COUNTRY CODE NOT RECOGNIZED ") - endif - endif - sumemis_local(ic,iem)= sumemis_local(ic,iem)+0.001*snapemis(isec,i,j,n,iem)!for diagnostics, mass balance - enddo + if(Country(ic)%index/=landcode(i,j,n))then + !if not equal, find which index correspond to country_code + do ic=1,NLAND + if(Country(ic)%index==landcode(i,j,n))exit enddo - enddo - - - enddo!sectors - CALL MPI_REDUCE(sumemis_local(1,iem),sumemis(1,iem),NLAND,MPI_REAL8,MPI_SUM,0,MPI_COMM_WORLD,INFO) - emsum(iem)= sum(sumemis(:,iem)) - - if(EMIS_OUT)call EmisOut ("Frac", iem, nlandcode, landcode, snapemis(:,:,:,:,iem) ) ! cdfemis - else - call CheckStop("EMIS_SOURCE not set"//trim(EMIS_SOURCE)) - endif - - - end do ! iem = 1, NEMIS_FILE-loop - - if(USE_ROADDUST) then - - if(USE_OLDSCHEME_ROADDUST)then - !use scheme with ASCII and grid dependent input data - - ! First Temporary/Test handling of climate factors. Read from file. Should be enough to do this on MasterProc - - if ( MasterProc )then - ! if(.true.)then - if (DEBUG_ROADDUST) write(unit=6,fmt=*) "Setting road dust climate scaling factors from", & - trim(roaddust_climate_file) - fname = roaddust_climate_file - call open_file(IO_EMIS,"r",fname,needed=.true.) - call CheckStop(ios,"RoadDust climate file: ios error in emission file") - - read(unit=IO_EMIS,fmt="(a200)",iostat=ios) inputline - if( inputline(1:1) .ne. "#" ) then ! Is a comment - write(*,*)'ERROR in road dust climate factor file!' - write(*,*)'First line should be a comment line, starting with #' - else - IF(DEBUG_ROADDUST) write(*,*)'I read the comment line:',inputline - endif - - READCLIMATEFACTOR: do ! ************* Loop over emislist files ******************* - - read(unit=IO_EMIS,fmt=*,iostat=ios) i,j, tmpclimfactor - - if ( ios < 0 ) exit READCLIMATEFACTOR ! End of file - call CheckStop(ios > 0,"RoadDust climate file: ios error2 in climate factor file") - - i = i-IRUNBEG+1 ! for RESTRICTED domain - j = j-JRUNBEG+1 ! for RESTRICTED domain - - if ( i <= 0 .or. i > GIMAX .or. & - j <= 0 .or. j > GJMAX )& - cycle READCLIMATEFACTOR - - RoadDustEmis_climate_factor(i,j) = tmpclimfactor - - - if( DEBUG_ROADDUST .and. i==DEBUG_i .and. j==DEBUG_j ) write(*,*) & - "DEBUG RoadDust climate factor (read from file)", RoadDustEmis_climate_factor(i,j) - - enddo READCLIMATEFACTOR - ! else - ! write(unit=6,fmt=*) "Test run set road dust climate factor to 1 everywhere!" - ! RoadDustEmis_climate_factor(:,:) = 1.0 - ! endif - - endif !MasterProc - - do iem = 1, NROAD_FILES - ! now again test for me=0 - if ( MasterProc ) then - - ! read in road dust emission potentials from one file - ! There will be two road dust files; one for highways(plus some extras) and one for non-highways - ! Each file contains spring (Mar-May) and rest-of-the-year (June-February) emission potentials - ! However, this will be changed so that only "rest-of-the-year" data are read and the spring scaling - ! for the Nordic countries is handled in the EmisSet routine! - ! Emission potentials are for PM10 so should be split into PM-fine (10%) and PM-coarse (90%) - ! There should also be some modifications to take into account temporal variations due to different - ! traffic intensities and differences due to surface wetness (handled in EmisSet) - ! and a climatological factor (handled here). - - ! ***************** - call RoadDustGet(iem,ROAD_FILE(iem),IRUNBEG,JRUNBEG,GIMAX,GJMAX, & - sumroaddust,& - globroad_dust_pot,road_globnland,road_globland) - ! ***************** - - roaddustsum(iem) = sum( globroad_dust_pot(:,:,:) ) ! - - ! ToDo-2012-0913 - - do i=1,GIMAX - do j=1,GJMAX - if(DEBUG_ROADDUST)then - WRITE(*,*)"i,j,RDECF:",i,j,RoadDustEmis_climate_factor(i,j) - endif - globroad_dust_pot(i,j,:)=RoadDustEmis_climate_factor(i,j)*globroad_dust_pot(i,j,:) - enddo - enddo - - - endif ! MasterProc - - call CheckStop(ios, "ios error: RoadDustGet") - - ! Send data to processors - ! as e.g. snapemis (NSECTORS,MAXLIMAX,MAXLJMAX,NCMAX,iem) - ! send to nodes - - call global2local(globroad_dust_pot,roaddust_emis_pot(1,1,1,iem),MSG_READ1, & - 1,GIMAX,GJMAX,NCMAX,1,1) - - end do ! iem = 1, NROAD_FILES-loop + if(ic>NLAND)then + write(*,*)"COUNTRY CODE NOT RECOGNIZED OR UNDEFINED: ",landcode(i,j,n) + call StopAll("COUNTRY CODE NOT RECOGNIZED ") + endif + endif + sumemis_local(ic,iem)=sumemis_local(ic,iem)& + +0.001*snapemis(isec,i,j,n,iem)!for diagnostics, mass balance + enddo + enddo + enddo + enddo!sectors + CALL MPI_REDUCE(sumemis_local(1,iem),sumemis(1,iem),NLAND,MPI_REAL8,& + MPI_SUM,0,MPI_COMM_WORLD,INFO) + emsum(iem)= sum(sumemis(:,iem)) + + if(EMIS_OUT)& + call EmisOut("Frac",iem,nlandcode,landcode,snapemis(:,:,:,:,iem)) !cdfemis + + case default + call CheckStop("EMIS_SOURCE not set"//trim(EMIS_SOURCE)) + endselect + enddo ! iem = 1, NEMIS_FILE-loop + + if(USE_ROADDUST) then + if(USE_OLDSCHEME_ROADDUST)then + !use scheme with ASCII and grid dependent input data + ! First Temporary/Test handling of climate factors. Read from file. Should be enough to do this on MasterProc + if(MasterProc)then +! if(.true.)then + if(DEBUG_ROADDUST) write(*,*) & + "Setting road dust climate scaling factors from",trim(roaddust_climate_file) + fname = roaddust_climate_file + call open_file(IO_EMIS,"r",fname,needed=.true.) + call CheckStop(ios,"RoadDust climate file: ios error in emission file") + + read(unit=IO_EMIS,fmt="(a200)",iostat=ios) inputline + if(inputline(1:1).ne."#") then ! Is a comment + write(*,*)'ERROR in road dust climate factor file!' + write(*,*)'First line should be a comment line, starting with #' + else + IF(DEBUG_ROADDUST) write(*,*)'I read the comment line:',inputline + endif - - call global2local_int(road_globnland,road_nlandcode,326,& - GIMAX,GJMAX,1,1,1) !extra array - call global2local_int(road_globland,road_landcode,326,& - GIMAX,GJMAX,NCMAX,1,1) - - - else - !Use grid-independent Netcdf input files - - do iem = 1, NROAD_FILES - !Read data from NetCDF file - if(iem==1)varname='HighwayRoadDustPM10_Jun-Feb' - if(iem==2)varname='nonHighwayRoadDustPM10_Jun-Feb' - call CheckStop(iem>2, "TOO MANY ROADFILES") - - roaddust_emis_pot(:,:,:,iem)=0.0 - call ReadField_CDF('RoadMap.nc',varname,roaddust_emis_pot(1,1,1,iem),& - nstart=1,interpol='mass_conservative', & - fractions_out=fractions,CC_out=road_landcode,Ncc_out=road_nlandcode,needed=.true.,debug_flag=.true.,Undef=0.0) - if(.not.SMI_defined)then - varname='SMI1' - call ReadField_CDF('AVG_SMI_2005_2010.nc',varname,SMI,nstart=1,& - interpol='conservative',needed=.true.,debug_flag=.true.) - SMI_defined=.true. - endif - - do i=1,LIMAX - do j=1,LJMAX - SMI_roadfactor=3.325-(min(1.0,max(0.5,SMI(i,j)))-0.5) *2*(3.325-1.0)!Peter: Rough estimate to get something varying between 3.325 (SMI<0.5) and 1.0 (SMI>1) - if(DEBUG_ROADDUST)then - ! WRITE(*,*)"i,j,RDECF:",i_fdom(i)-IRUNBEG+1,j_fdom(j)-JRUNBEG+1,SMI_roadfactor - endif - do iic=road_nlandcode(i,j),1,-1 - roaddust_emis_pot(i,j,iic,iem)=roaddust_emis_pot(i,j,1,iem)*fractions(i,j,iic)*SMI_roadfactor - enddo - enddo + ! ************* Loop over emislist files ******************* + READCLIMATEFACTOR: do + read(IO_EMIS,*,iostat=ios) i,j, tmpclimfactor + if(ios<0) exit READCLIMATEFACTOR ! End of file + call CheckStop(ios>0,"RoadDust climate file: ios error2 in climate factor file") + + i = i-IRUNBEG+1 ! for RESTRICTED domain + j = j-JRUNBEG+1 ! for RESTRICTED domain + if(i<=0 .or. i>GIMAX .or. & + j<=0 .or. j>GJMAX) cycle READCLIMATEFACTOR + + RoadDustEmis_climate_factor(i,j) = tmpclimfactor + if(DEBUG_ROADDUST.and.i==DEBUG_i.and.j==DEBUG_j) write(*,*) & + "DEBUG RoadDust climate factor (read from file)",& + RoadDustEmis_climate_factor(i,j) + + enddo READCLIMATEFACTOR +! else +! write(*,*) "Test run set road dust climate factor to 1 everywhere!" +! RoadDustEmis_climate_factor(:,:) = 1.0 +! endif + endif !MasterProc + + do iem=1,NROAD_FILES + ! now again test for me=0 + if(MasterProc) then +! read in road dust emission potentials from one file +! There will be two road dust files; one for highways(plus some extras) and one for non-highways +! Each file contains spring (Mar-May) and rest-of-the-year (June-February) emission potentials +! However, this will be changed so that only "rest-of-the-year" data are read and the spring scaling +! for the Nordic countries is handled in the EmisSet routine! +! Emission potentials are for PM10 so should be split into PM-fine (10%) and PM-coarse (90%) +! There should also be some modifications to take into account temporal variations due to different +! traffic intensities and differences due to surface wetness (handled in EmisSet) +! and a climatological factor (handled here). + ! ***************** + call RoadDustGet(iem,ROAD_FILE(iem),IRUNBEG,JRUNBEG,GIMAX,GJMAX, & + sumroaddust,globroad_dust_pot,road_globnland,road_globland) + ! ***************** + roaddustsum(iem) = sum(globroad_dust_pot(:,:,:)) ! + + ! ToDo-2012-0913 + do i=1,GIMAX + do j=1,GJMAX + if(DEBUG_ROADDUST)& + WRITE(*,*)"i,j,RDECF:",i,j,RoadDustEmis_climate_factor(i,j) + globroad_dust_pot(i,j,:)=RoadDustEmis_climate_factor(i,j)*globroad_dust_pot(i,j,:) enddo - sumroaddust_local(:,iem)=0.0 - do i=1,LIMAX - do j=1,LJMAX - do iic=1,road_nlandcode(i,j) - ic=1 - if(road_landcode(i,j,iic)<=NLAND)ic=road_landcode(i,j,iic)!most country_index are equal to country_code - if(Country(ic)%index/=road_landcode(i,j,iic))then - !if not equal, find which index correspond to country_code - do ic=1,NLAND - if((Country(ic)%index==road_landcode(i,j,iic)))exit - enddo - if(ic>NLAND)then - write(*,*)"COUNTRY CODE NOT RECOGNIZED OR UNDEFINED: ", road_landcode(i,j,iic) - call StopAll("COUNTRY CODE NOT RECOGNIZED ") - - endif - endif - - sumroaddust_local(ic,iem)=sumroaddust_local(ic,iem)+0.001*roaddust_emis_pot(i,j,iic,iem) - - enddo + enddo + endif ! MasterProc + call CheckStop(ios, "ios error: RoadDustGet") + + ! Send data to processors as + ! e.g. snapemis (NSECTORS,MAXLIMAX,MAXLJMAX,NCMAX,iem) send to nodes + call global2local(globroad_dust_pot,roaddust_emis_pot(1,1,1,iem),& + MSG_READ1,1,GIMAX,GJMAX,NCMAX,1,1) + enddo ! iem = 1, NROAD_FILES-loop + + call global2local_int(road_globnland,road_nlandcode,326,GIMAX,GJMAX,1,1,1) !extra array + call global2local_int(road_globland,road_landcode,326,GIMAX,GJMAX,NCMAX,1,1) + + else + !Use grid-independent Netcdf input files + call CheckStop(NROAD_FILES>2, "TOO MANY ROADFILES") + do iem = 1, NROAD_FILES + !Read data from NetCDF file + select case(iem) + case(1);varname='HighwayRoadDustPM10_Jun-Feb' + case(2);varname='nonHighwayRoadDustPM10_Jun-Feb' + endselect + roaddust_emis_pot(:,:,:,iem)=0.0 + call ReadField_CDF('RoadMap.nc',varname,roaddust_emis_pot(1,1,1,iem),& + nstart=1,interpol='mass_conservative',fractions_out=fractions,& + CC_out=road_landcode,Ncc_out=road_nlandcode,needed=.true.,& + debug_flag=.true.,Undef=0.0) + if(.not.SMI_defined)then + varname='SMI1' + call ReadField_CDF('AVG_SMI_2005_2010.nc',varname,SMI,nstart=1,& + interpol='conservative',needed=.true.,debug_flag=.true.) + SMI_defined=.true. + endif + + do i=1,LIMAX + do j=1,LJMAX + !Peter: Rough estimate to get something varying between 3.325 (SMI<0.5) and 1.0 (SMI>1) + SMI_roadfactor=3.325-(min(1.0,max(0.5,SMI(i,j)))-0.5)*2*(3.325-1.0) + !if(DEBUG_ROADDUST)& + ! WRITE(*,*)"i,j,RDECF:",i_fdom(i)-IRUNBEG+1,j_fdom(j)-JRUNBEG+1,SMI_roadfactor + do iic=road_nlandcode(i,j),1,-1 + roaddust_emis_pot(i,j,iic,iem)=roaddust_emis_pot(i,j,1,iem) & + *fractions(i,j,iic)*SMI_roadfactor + enddo + enddo + enddo + sumroaddust_local(:,iem)=0.0 + do i=1,LIMAX + do j=1,LJMAX + do iic=1,road_nlandcode(i,j) + ic=1 + if(road_landcode(i,j,iic)<=NLAND) & + ic=road_landcode(i,j,iic) ! most country_index are equal to country_code + if(Country(ic)%index/=road_landcode(i,j,iic))then + !if not equal, find which index correspond to country_code + do ic=1,NLAND + if((Country(ic)%index==road_landcode(i,j,iic)))exit enddo + if(ic>NLAND)then + write(*,*)"COUNTRY CODE NOT RECOGNIZED OR UNDEFINED: ", & + road_landcode(i,j,iic) + call StopAll("COUNTRY CODE NOT RECOGNIZED ") + endif + endif + sumroaddust_local(ic,iem)=sumroaddust_local(ic,iem)& + +0.001*roaddust_emis_pot(i,j,iic,iem) enddo - - end do ! iem = 1, NROAD_FILES-loop - - sumroaddust=0.0 - CALL MPI_REDUCE(sumroaddust_local,sumroaddust,NLAND*NROAD_FILES,MPI_REAL8,MPI_SUM,0,MPI_COMM_WORLD,INFO) - - endif - - end if !USE_ROADDUST - - if ( MasterProc ) then - write(unit=6,fmt=*) "Total emissions by countries:" - write(unit=IO_LOG,fmt=*) "Total emissions by countries:" - write(unit=6,fmt="(2a4,3x,30a12)") " N "," CC ",(EMIS_FILE(iem),iem=1,NEMIS_FILE) - write(unit=IO_LOG,fmt="(2a4,3x,30a12)") " N "," CC ",(EMIS_FILE(iem),iem=1,NEMIS_FILE) - - sumEU(:) = 0.0 - do ic = 1, NLAND - ccsum = sum( sumemis(ic,:) ) - !if ( ccsum > 0.0 ) then - if ( ccsum > 0.0 .or. sum( sumcdfemis(ic,:)) > 0.0 ) then - if ( EMIS_TEST == "None" ) then - write(unit=6,fmt="(i3,1x,a4,3x,30f12.2)") & - ic, Country(ic)%code, (sumemis(ic,i),i=1,NEMIS_FILE) - else - - write(unit=6,fmt="(a,i3,1x,a4,3x,30f12.2)") & - "ORIG:",ic, Country(ic)%code, (sumemis(ic,i),i=1,NEMIS_FILE) - write(unit=6,fmt="(a,i3,1x,a4,3x,30f12.2)") & !cdfemis - "CDFS:",ic, Country(ic)%code, (sumcdfemis(ic,i),i=1,NEMIS_FILE) - end if - write(unit=IO_LOG,fmt="(i3,1x,a4,3x,30f12.2)")& - ic, Country(ic)%code, (sumemis(ic,i),i=1,NEMIS_FILE) - if( find_index( Country(ic)%code , EU27(:) ) > 0 ) then - sumEU = sumEU + sumemis(ic,:) - end if - end if - end do - write(unit=6 ,fmt="(i3,1x,a4,3x,30f12.2)") 0, "EU", (sumEU(i),i=1,NEMIS_FILE) - write(unit=IO_LOG,fmt="(i3,1x,a4,3x,30f12.2)") 0, "EU", (sumEU(i),i=1,NEMIS_FILE) - !GP_work, TMP solution with hard-coded emissions - iem= find_index( "nox", EMIS_FILE(:) ) - iem2=find_index( "nh3", EMIS_FILE(:) ) - Ndep_trends = ( sumEU(iem) + sumEU(iem2) ) / 15167.48 - write(unit=6,fmt="(a,2i3,f12.2,f8.4)") "Ndep_trends", & - iem, iem2, sumEU(iem) + sumEU(iem2), Ndep_trends - write(unit=IO_LOG,fmt="(a,f8.4)") "Ndep_trends", Ndep_trends + enddo + enddo + enddo ! iem = 1, NROAD_FILES-loop + sumroaddust=0.0 + CALL MPI_REDUCE(sumroaddust_local,sumroaddust,NLAND*NROAD_FILES,MPI_REAL8,& + MPI_SUM,0,MPI_COMM_WORLD,INFO) + endif + endif !USE_ROADDUST + + if(MasterProc) then + call PrintLog("Total emissions by countries:") + write(* ,"(2a4,3x,30(a12,:))")" N "," CC ",EMIS_FILE(:) + write(IO_LOG,"(2a4,3x,30(a12,:))")" N "," CC ",EMIS_FILE(:) + + sumEU(:) = 0.0 + do ic = 1, NLAND + ccsum = sum( sumemis(ic,:) ) + !if ( ccsum > 0.0 ) then + if(ccsum>0.0 .or. sum(sumcdfemis(ic,:))>0.0) then + if(EMIS_TEST=="None") then + write(*,"(i3,1x,a4,3x,30(f12.2,:))")ic, Country(ic)%code, sumemis(ic,:) + else + write(*,"(a,i3,1x,a4,3x,30(f12.2,:))")"ORIG:",ic, Country(ic)%code, sumemis(ic,:) + write(*,"(a,i3,1x,a4,3x,30(f12.2,:))")"CDFS:",ic, Country(ic)%code, sumcdfemis(ic,:) + endif + write(IO_LOG,"(i3,1x,a4,3x,30(f12.2,:))")ic, Country(ic)%code, sumemis(ic,:) + if(find_index(Country(ic)%code,EU27(:))>0) sumEU = sumEU + sumemis(ic,:) + endif + enddo + write(* ,"(i3,1x,a4,3x,30(f12.2,:))") 0, "EU", sumEU(:) + write(IO_LOG,"(i3,1x,a4,3x,30(f12.2,:))") 0, "EU", sumEU(:) - if(USE_ROADDUST)THEN - - write(unit=6,fmt=*) "Total road dust emission potentials by countries (before precipitation and land corrections):" - write(unit=IO_LOG,fmt=*) "Total road dust emission potentials by countries (before precipitation and land corrections):" - write(unit=6,fmt="(2a4,11x,30a12)") " N "," CC ",(ROAD_FILE(iem),iem=1,NROAD_FILES) - write(unit=IO_LOG,fmt="(2a4,11x,30a12)") " N "," CC ",(ROAD_FILE(iem),iem=1,NROAD_FILES) - - do ic = 1, NLAND - ccsum = sum( sumroaddust(ic,:) ) - if ( ccsum > 0.0 ) then - write(unit=6,fmt="(i3,1x,a4,3x,30f12.2)") & - ic, Country(ic)%code, (sumroaddust(ic,i),i=1,NROAD_FILES) - write(unit=IO_LOG,fmt="(i3,1x,a4,3x,30f12.2)")& - ic, Country(ic)%code, (sumroaddust(ic,i),i=1,NROAD_FILES) - end if - end do - endif ! ROAD DUST - end if - - ! now all values are read, snapemis is distributed, globnland and - ! globland are ready for distribution - ! print *, "calling glob2local_int for iem", iem, " me ", me - - if( EMIS_SOURCE == "emislist" ) then - call global2local_int(globnland,nlandcode,326, GIMAX,GJMAX,1,1,1) - call global2local_int(globland, landcode ,326, GIMAX,GJMAX,NCMAX,1,1) - - call global2local_int(flat_globnland,flat_nlandcode,326,& - GIMAX,GJMAX,1,1,1) !extra array - call global2local_int(flat_globland,flat_landcode,326,& - GIMAX,GJMAX,FNCMAX,1,1) - else - !EMIS_SOURCE = CdfFractions - !emissions directly defined into - !nlandcode,landcode and snapemis - endif + if(USE_ROADDUST)THEN + call PrintLog("Total road dust emission potentials by countries & + &(before precipitation and land corrections):") + write(* ,"(2a4,11x,30(a12,:))")" N "," CC ",ROAD_FILE(:) + write(IO_LOG,"(2a4,11x,30(a12,:))")" N "," CC ",ROAD_FILE(:) + + do ic = 1, NLAND + ccsum = sum( sumroaddust(ic,:)) + if(ccsum>0.0) then + write(* ,"(i3,1x,a4,3x,30(f12.2,:))")ic, Country(ic)%code, sumroaddust(ic,:) + write(IO_LOG,"(i3,1x,a4,3x,30(f12.2,:))")ic, Country(ic)%code, sumroaddust(ic,:) + endif + enddo + endif ! ROAD DUST + endif + + ! now all values are read, snapemis is distributed, globnland and + ! globland are ready for distribution + ! print *, "calling glob2local_int for iem", iem, " me ", me + select case(EMIS_SOURCE) + case("emislist") + call global2local_int(globnland,nlandcode,326,GIMAX,GJMAX,1,1,1) + call global2local_int(globland ,landcode ,326,GIMAX,GJMAX,NCMAX,1,1) + + call global2local_int(flat_globnland,flat_nlandcode,326,& + GIMAX,GJMAX,1,1,1) !extra array + call global2local_int(flat_globland,flat_landcode,326,& + GIMAX,GJMAX,FNCMAX,1,1) + case("CdfFractions") + ! emissions directly defined into nlandcode,landcode and snapemis + endselect ! Create emislist-type files for both snap emissions and Cdf ! Useful for export to other codes, including production of ! new emislist for current NWP grid. - do iem = 1, NEMIS_FILE + do iem = 1, NEMIS_FILE + if((EMIS_SOURCE=="emislist").and.EMIS_OUT) & + call EmisOut("Snap",iem,nlandcode,landcode,snapemis(:,:,:,:,iem)) + + if(EMIS_TEST=="CdfSnap")& + write(*,"(a,i3,2i4,2es12.3)") "CALLED CDF PRE "//& + trim(EMIS_FILE(iem)),me,maxval(nlandcode),maxval(nGridEmisCodes), & + maxval(snapemis(:,:,:,:,iem)),maxval(GridEmis(:,:,:,:,iem)) + + if((EMIS_TEST=="CdfSnap").and.EMIS_OUT) & + call EmisOut("Cdf",iem,nGridEmisCodes,GridEmisCodes,GridEmis(:,:,:,:,iem)) + enddo - if( EMIS_SOURCE == "emislist" ) then - if(EMIS_OUT)call EmisOut ("Snap", iem, nlandcode, landcode, snapemis(:,:,:,:,iem) ) - end if - - if( EMIS_TEST == "CdfSnap" ) then - write(6, "(a,i3,2i4,2es12.3)") "CALLED CDF PRE "//& - trim(EMIS_FILE(iem)), me, & - maxval( nlandcode ),maxval( nGridEmisCodes ), & - maxval( snapemis(:,:,:,:,iem)), maxval( GridEmis(:,:,:,:,iem)) - if(EMIS_OUT)call EmisOut ("Cdf", iem, nGridEmisCodes, & - GridEmisCodes, GridEmis(:,:,:,:,iem) ) - end if - end do - - !------------------------------------------------------------------------- - ! Broadcast volcanoe info derived in EmisGet - - CALL MPI_BCAST(nvolc,4*1,MPI_BYTE,0,MPI_COMM_WORLD,INFO) - CALL MPI_BCAST(i_volc,4*nvolc,MPI_BYTE,0,MPI_COMM_WORLD,INFO) - CALL MPI_BCAST(j_volc,4*nvolc,MPI_BYTE,0,MPI_COMM_WORLD,INFO) - CALL MPI_BCAST(emis_volc,8*nvolc,MPI_BYTE,0,MPI_COMM_WORLD,INFO) + !------------------------------------------------------------------------- + ! Broadcast volcanoe info derived in EmisGet + CALL MPI_BCAST(nvolc,4*1,MPI_BYTE,0,MPI_COMM_WORLD,INFO) + CALL MPI_BCAST(i_volc,4*nvolc,MPI_BYTE,0,MPI_COMM_WORLD,INFO) + CALL MPI_BCAST(j_volc,4*nvolc,MPI_BYTE,0,MPI_COMM_WORLD,INFO) + CALL MPI_BCAST(emis_volc,8*nvolc,MPI_BYTE,0,MPI_COMM_WORLD,INFO) - ! And N-dep ttrends - CALL MPI_BCAST(Ndep_trends,8,MPI_BYTE,0,MPI_COMM_WORLD,INFO) - - ! Conversions - ! - ! The emission-data file are so far in units of - ! tonnes per grid-square. The conversion factor from tonnes per 50*50km2 - ! annual emission values to surface flux (kg/m2/s) is found by division - ! with (nydays*24*60*60)s and (h*h)m2 and multiply by 1.e+3. - ! The conversion factor then equals 1.27e-14 - - tonne_to_kgm2s = 1.0e3 / (nydays * 24.0 * 3600.0 * & - GRIDWIDTH_M * GRIDWIDTH_M) - - if ( DEBUG .and. MasterProc ) then - write(unit=6,fmt=*) "CONV:me, nydays, gridwidth = ",me,nydays,GRIDWIDTH_M - write(unit=6,fmt=*) "No. days in Emissions: ", nydays - write(unit=6,fmt=*) "tonne_to_kgm2s in Emissions: ", tonne_to_kgm2s - write(unit=6,fmt=*) "Emissions sums:" - do iem = 1, NEMIS_FILE - write(unit=6,fmt="(a15,f12.2)") EMIS_FILE(iem),emsum(iem) - end do - endif - + !** Conversions: + ! The emission-data file are so far in units of + ! tonnes per grid-square. The conversion factor from tonnes per 50*50km2 + ! annual emission values to surface flux (kg/m2/s) is found by division + ! with (nydays*24*60*60)s and (h*h)m2 and multiply by 1.e+3. + ! The conversion factor then equals 1.27e-14 + tonne_to_kgm2s = 1.0e3 / (nydays * 24.0 * 3600.0 * GRIDWIDTH_M * GRIDWIDTH_M) + if(DEBUG.and.MasterProc) then + write(*,*) "CONV:me, nydays, gridwidth = ",me,nydays,GRIDWIDTH_M + write(*,*) "No. days in Emissions: ", nydays + write(*,*) "tonne_to_kgm2s in Emissions: ", tonne_to_kgm2s + write(*,*) "Emissions sums:" do iem = 1, NEMIS_FILE - conv = tonne_to_kgm2s - if ( trim(EMIS_FILE(iem)) == "co" ) iemCO = iem ! save this index - - if ( DEBUG .and. debug_proc .and. iem == iemCO ) then - write(*,"(a,2es10.3)") "SnapPre:" // trim(EMIS_FILE(iem)), & - sum( snapemis (:,debug_li,debug_lj,:,iem) ) & - ,sum( snapemis_flat (debug_li,debug_lj,:,iem) ) - end if - - forall (ic=1:NCMAX, j=1:ljmax, i=1:limax, isec=1:NSECTORS) - snapemis (isec,i,j,ic,iem) = & - snapemis (isec,i,j,ic,iem) * conv * xm2(i,j) - end forall - - forall (fic=1:FNCMAX, j=1:ljmax, i=1:limax) - snapemis_flat(i,j,fic,iem) = & - snapemis_flat(i,j,fic,iem) * conv * xm2(i,j) - end forall - - - if ( DEBUG .and. debug_proc .and. iem == iemCO ) then - write(*,"(a,2es10.3)") "SnapPos:" // trim(EMIS_FILE(iem)), & - sum( snapemis (:,debug_li,debug_lj,:,iem) ) & - ,sum( snapemis_flat (debug_li,debug_lj,:,iem) ) - end if - enddo !iem - - if(USE_ROADDUST)THEN - do iem = 1, NROAD_FILES - conv = tonne_to_kgm2s - - forall (ic=1:NCMAX, j=1:ljmax, i=1:limax) - roaddust_emis_pot(i,j,ic,iem) = & - roaddust_emis_pot(i,j,ic,iem) * conv * xm2(i,j) - end forall - enddo ! iem - endif !road dust - - - ! if ( VOLCANOES ) then - - ! Read Volcanos.dat or VolcanoesLL.dat to get volcano height - ! and magnitude in the case of VolcanoesLL.dat - call VolcGet(height_volc) - - ! endif ! VOLCANOES - - err1 = 0 - if ( MasterProc ) then - deallocate(globnland,stat=err1) - deallocate(globland ,stat=err2) - deallocate(globemis ,stat=err3) - - deallocate(flat_globnland,stat=err4) - deallocate(flat_globland,stat=err5) - deallocate(globemis_flat,stat=err6) - - if(USE_ROADDUST)THEN - deallocate(road_globnland,stat=err7) - deallocate(road_globland,stat=err8) - deallocate(globroad_dust_pot,stat=err9) - endif - - call CheckStop(err1, "De-Allocation error 1 - globland") - call CheckStop(err2, "De-Allocation error 2 - globland") - call CheckStop(err3, "De-Allocation error 3 - globland") - call CheckStop(err4, "De-Allocation error 4 - globland") - call CheckStop(err5, "De-Allocation error 5 - globland") - call CheckStop(err6, "De-Allocation error 6 - globland") - if(USE_ROADDUST)THEN - call CheckStop(err7, "De-Allocation error 7 - roadglob") - call CheckStop(err8, "De-Allocation error 8 - roadglob") - call CheckStop(err9, "De-Allocation error 9 - roadglob") - endif - end if - - ! now we have nrecmis and can allocate for gridrcemis: - ! print *, "ALLOCATING GRIDRC", me, NRCEMIS - allocate(gridrcemis(NRCEMIS,KEMISTOP:KMAX_MID,MAXLIMAX,MAXLJMAX),stat=err1) - allocate(gridrcemis0(NRCEMIS,KEMISTOP:KMAX_MID,MAXLIMAX,MAXLJMAX),stat=err2) - call CheckStop(err1, "Allocation error 1 - gridrcemis") - call CheckStop(err2, "Allocation error 2 - gridrcemis0") + write(*,"(a15,f12.2)") EMIS_FILE(iem),emsum(iem) + enddo + endif + + iemCO=find_index("co",EMIS_FILE(:)) ! save this index + conv = tonne_to_kgm2s + + if(DEBUG.and.debug_proc.and.iemCO>0) & + write(*,"(a,2es10.3)") "SnapPre:" // trim(EMIS_FILE(iemCO)), & + sum(snapemis (:,debug_li,debug_lj,:,iemCO)), & + sum(snapemis_flat(debug_li,debug_lj,:,iemCO)) + + forall (ic=1:NCMAX, j=1:ljmax, i=1:limax, isec=1:NSECTORS,iem=1:NEMIS_FILE) + snapemis (isec,i,j,ic,iem) = snapemis (isec,i,j,ic,iem) * conv * xm2(i,j) + endforall + + forall (fic=1:FNCMAX, j=1:ljmax, i=1:limax,iem=1:NEMIS_FILE) + snapemis_flat(i,j,fic,iem) = snapemis_flat(i,j,fic,iem) * conv * xm2(i,j) + endforall + + if(DEBUG.and.debug_proc.and.iemCO>0) & + write(*,"(a,2es10.3)") "SnapPos:" // trim(EMIS_FILE(iem)), & + sum(snapemis (:,debug_li,debug_lj,:,iem)), & + sum(snapemis_flat(debug_li,debug_lj,:,iem)) + + if(USE_ROADDUST)THEN + conv = tonne_to_kgm2s + forall (ic=1:NCMAX, j=1:ljmax, i=1:limax, iem=1:NROAD_FILES) + roaddust_emis_pot(i,j,ic,iem) = & + roaddust_emis_pot(i,j,ic,iem) * conv * xm2(i,j) + endforall + endif !road dust + + +!if ( VOLCANOES ) then + ! Read Volcanos.dat or VolcanoesLL.dat to get volcano height + ! and magnitude in the case of VolcanoesLL.dat + call VolcGet(height_volc) +!endif ! VOLCANOES + + err1 = 0 + if(MasterProc) then + deallocate(globnland ,stat=err1) + deallocate(globland ,stat=err2) + deallocate(globemis ,stat=err3) + deallocate(flat_globnland,stat=err4) + deallocate(flat_globland ,stat=err5) + deallocate(globemis_flat ,stat=err6) if(USE_ROADDUST)THEN - allocate(gridrcroadd(NROADDUST,MAXLIMAX,MAXLJMAX),stat=err3) - allocate(gridrcroadd0(NROADDUST,MAXLIMAX,MAXLJMAX),stat=err4) - call CheckStop(err3, "Allocation error 3 - gridrcroadd") - call CheckStop(err4, "Allocation error 4 - gridrcroadd0") + deallocate(road_globnland ,stat=err7) + deallocate(road_globland ,stat=err8) + deallocate(globroad_dust_pot,stat=err9) endif - end subroutine Emissions - - !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> - subroutine consistency_check() - !------------------------------------------------------------------! - ! checks that all the values given so far are consistent ! - !------------------------------------------------------------------! - - character(len=30) :: errormsg - + call CheckStop(err1, "De-Allocation error 1 - globland") + call CheckStop(err2, "De-Allocation error 2 - globland") + call CheckStop(err3, "De-Allocation error 3 - globland") + call CheckStop(err4, "De-Allocation error 4 - globland") + call CheckStop(err5, "De-Allocation error 5 - globland") + call CheckStop(err6, "De-Allocation error 6 - globland") + if(USE_ROADDUST)THEN + call CheckStop(err7, "De-Allocation error 7 - roadglob") + call CheckStop(err8, "De-Allocation error 8 - roadglob") + call CheckStop(err9, "De-Allocation error 9 - roadglob") + endif + endif + + ! now we have nrecmis and can allocate for gridrcemis: + ! print *, "ALLOCATING GRIDRC", me, NRCEMIS + allocate(gridrcemis(NRCEMIS,KEMISTOP:KMAX_MID,MAXLIMAX,MAXLJMAX),stat=err1) + allocate(gridrcemis0(NRCEMIS,KEMISTOP:KMAX_MID,MAXLIMAX,MAXLJMAX),stat=err2) + call CheckStop(err1, "Allocation error 1 - gridrcemis") + call CheckStop(err2, "Allocation error 2 - gridrcemis0") + if(USE_ROADDUST)THEN + allocate(gridrcroadd(NROADDUST,MAXLIMAX,MAXLJMAX),stat=err3) + allocate(gridrcroadd0(NROADDUST,MAXLIMAX,MAXLJMAX),stat=err4) + call CheckStop(err3, "Allocation error 3 - gridrcroadd") + call CheckStop(err4, "Allocation error 4 - gridrcroadd0") + endif +endsubroutine Emissions +!*********************************************************************** +subroutine consistency_check() +!----------------------------------------------------------------------! +! checks that all the values given so far are consistent +!----------------------------------------------------------------------! + character(len=30) :: errormsg errormsg = "ok" - if ( size(EMIS_FILE) /= NEMIS_FILE ) errormsg = " size EMISNAME wrong " - + if(size(EMIS_FILE)/=NEMIS_FILE) errormsg = " size EMISNAME wrong " call CheckStop(errormsg,"Failed consistency check") - - end subroutine consistency_check - !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< - - - - !***************************************************************** - +endsubroutine consistency_check +!*********************************************************************** subroutine EmisSet(indate) ! emission re-set every time-step/hour - - ! - !*********************************************************************** - ! DESCRIPTION: - ! Calculates the emmision-tendencies and the local (instantaneous) dry - ! deposition in the emission squares. - ! emis set once per hour to allow for day/night variation (and voc - ! speciation) (based on local time) for each snap sector. - ! gridrcemis0 calculated every time-step to allow for ps changes. - ! inputs from Emissions in EMISSIONS_ML: - ! country and snap-specific array : - ! snapemis (NSECTORS,MAXLIMAX,MAXLJMAX,NCMAX,NEMIS_FILES) - ! - ! Units: - ! snapemis has units of kg/m2/s, SO2 as S, NO2 as N, NH3 as N. - ! Map factor (xm2) already accounted for. - ! - ! Data on how many countries contribute per grid square is stored in - ! nlandcode(MAXLIMAX,MAXLJMAX) , with the country-codes given by - ! landcode(MAXLIMAX,MAXLJMAX,NCMAX). - ! - ! Monthly and weekday factors are pre-multiplied and stored in: - ! real timefac(NLAND,NSECTORS,NEMIS_FILES) - ! And day-hour factors in fac_ehh24x7 - ! - !************************************************************************* - +!----------------------------------------------------------------------! +! DESCRIPTION: +! Calculates the emmision-tendencies and the local (instantaneous) dry +! deposition in the emission squares. +! emis set once per hour to allow for day/night variation (and voc +! speciation) (based on local time) for each snap sector. +! gridrcemis0 calculated every time-step to allow for ps changes. +! inputs from Emissions in EMISSIONS_ML: +! country and snap-specific array : +! snapemis (NSECTORS,MAXLIMAX,MAXLJMAX,NCMAX,NEMIS_FILES) +! +! Units: +! snapemis has units of kg/m2/s, SO2 as S, NO2 as N, NH3 as N. +! Map factor (xm2) already accounted for. +! +! Data on how many countries contribute per grid square is stored in +! nlandcode(MAXLIMAX,MAXLJMAX) , with the country-codes given by +! landcode(MAXLIMAX,MAXLJMAX,NCMAX). +! +! Monthly and weekday factors are pre-multiplied and stored in: +! real timefac(NLAND,NSECTORS,NEMIS_FILES) +! And day-hour factors in fac_ehh24x7 +!----------------------------------------------------------------------! implicit none - type(date), intent(in) :: indate ! Gives year..seconds - integer :: i, j, k, f ! cooridnates, loop variables - integer :: icc, ncc ! No. of countries in grid. - - integer :: ficc,fncc ! No. of countries with - integer :: iqrc ! emis indices + type(date), intent(in) :: indate ! Gives year..seconds + integer :: i, j, k, f ! cooridnates, loop variables + integer :: icc, ncc ! No. of countries in grid. + integer :: ficc,fncc ! No. of countries with + integer :: iqrc ! emis indices integer :: isec ! loop variables: emission sectors integer :: iem ! loop variable over 1..NEMIS_FILE integer :: itot ! index in xn() @@ -1035,9 +893,7 @@ subroutine EmisSet(indate) ! emission re-set every time-step/hour real :: tfac, dtgrid ! time-factor (tmp variable); dt*h*h for scaling real :: s ! source term (emis) before splitting integer :: iland, iland_timefac ! country codes, and codes for timefac - - real :: sf ! source term (emis) before splitting - ! (for flat emissions) + real :: sf ! source term (emis) before splitting (flat emissions) integer :: flat_iland ! country codes (countries with flat emissions) integer, save :: oldday = -1, oldhour = -1 @@ -1049,270 +905,200 @@ subroutine EmisSet(indate) ! emission re-set every time-step/hour integer :: daytime_longitude, daytime_iland, hour_longitude, hour_iland ! Initialize - ehlpcom0 = GRAV* 0.001*AVOG!/ (sigma_bnd(k+1) - sigma_bnd(k)) - - ! Scaling for totemadd: - dtgrid = dt_advec * GRIDWIDTH_M * GRIDWIDTH_M - - ! The emis array only needs to be updated every full hour. The - ! time-factor calculation needs to know if a local-time causes a shift - ! from day to night. In addition, we reset an overall day's time-factors - ! at midnight every day. - - hourchange = .false. - if ( indate%hour /= oldhour .or. indate%day /= oldday ) then - - hourchange = .true. - oldhour = indate%hour - - if ( indate%day /= oldday )then - - !========================== - call NewDayFactors(indate) - if ( USE_DEGREEDAY_FACTORS) & - call DegreeDayFactors(daynumber) ! => fac_emm, fac_edd - !========================== - - ! for ROADDUST - wday=day_of_week(indate%year,indate%month,indate%day) - if(wday==0)wday=7 ! Sunday -> 7 - - oldday = indate%day - endif - end if - - if( DEBUG_EMISTIMEFACS .and. MasterProc ) then - write(*,"(a,2f8.3)") " EmisSet traffic 24x7", & - fac_ehh24x7(ISNAP_TRAF,1,4),fac_ehh24x7(ISNAP_TRAF,13,4) - end if - - - !.......................................... - ! Look for day-night changes, after local time correction - ! (daytime(iland) not used if LONGITUDE_TIME=true) - - do iland = 1, NLAND - - daytime(iland) = 0 - hourloc = indate%hour + Country(iland)%timezone - - localhour(iland) = hourloc ! here from 0 to 23 - - if ( hourloc >= 7 .and. hourloc <= 18 ) daytime(iland)=1 - - end do ! iland - - - if ( hourchange ) then - - totemadd(:) = 0. - gridrcemis0(:,:,:,:) = 0.0 - SumSnapEmis(:,:,:) = 0.0 - if(USE_ROADDUST)gridrcroadd0(:,:,:) = 0.0 - - !.......................................... - ! Process each grid: - - do j = 1,ljmax - do i = 1,limax - - ncc = nlandcode(i,j) ! No. of countries in grid - debug_tfac = ( DEBUG_EMISTIMEFACS .and. debug_proc .and. & - i==DEBUG_li .and. j==DEBUG_lj ) - - ! find the approximate local time: - hourloc= mod(nint(indate%hour+24*(1+glon(i,j)/360.0)),24) - hour_longitude=hourloc - daytime_longitude=0 - if( hourloc>=7.and.hourloc<= 18) daytime_longitude=1 - - !************************************************* - ! First loop over non-flat (one sector) emissions - !************************************************* - - tmpemis(:)=0. - do icc = 1, ncc - iland = landcode(i,j,icc) ! 1=Albania, etc. - iland_timefac = Country(iland)%timefac_index - - if(Country(iland)%timezone==-100)then - daytime_iland=daytime_longitude - hour_iland=hour_longitude + 1 ! add 1 to get 1..24 - else - daytime_iland=daytime(iland) - hour_iland=localhour(iland) + 1 - endif - !if( hour_iland > 24 ) hour_iland = 1 !DSA12 - wday_loc=wday - if( hour_iland > 24 ) then - hour_iland = hour_iland - 24 - wday_loc=wday + 1 - if(wday_loc==0)wday_loc=7 ! Sunday -> 7 - if(wday_loc>7 )wday_loc=1 - end if - - call CheckStop( hour_iland < 1, & - "ERROR: HOUR Zero in EmisSet") - - if( debug_tfac ) then - write(*,"(a,i4,2i3,i5,2i4,3x,4i3)") "EmisSet DAYS times ", daynumber, & - wday, wday_loc, iland, daytime_longitude, daytime_iland,& - hour_longitude, hour_iland, hourloc, Country(iland)%timezone - call datewrite("EmisSet DAY 24x7:", & - (/ icc, iland, wday, wday_loc, hour_iland /), & - (/ fac_ehh24x7(ISNAP_TRAF,hour_iland,wday_loc) /) ) - end if - - - ! As each emission sector has a different diurnal profile - ! and possibly speciation, we loop over each sector, adding - ! the found emission rates to gridrcemis as we go. - ! ================================================== - - - do isec = 1, NSECTORS ! Loop over snap codes - - ! Calculate emission rates from snapemis, time-factors, - ! and if appropriate any speciation fraction (NEMIS_FRAC) - - iqrc = 0 ! index over emisfrac - - do iem = 1, NEMIS_FILE - - tfac = timefac(iland_timefac,isec,iem) * & - fac_ehh24x7(isec,hour_iland,wday_loc) - - if( debug_tfac .and. iem == 1 ) then - write(*,"(a,2i4,f8.3)") "EmisSet DAY TFAC:", & - isec, hour_iland, tfac - end if - - !Degree days - only SNAP-2 - if ( USE_DEGREEDAY_FACTORS .and. & - isec == ISNAP_DOM .and. Gridded_SNAP2_Factors ) then - - oldtfac = tfac - - ! If INERIS_SNAP2 set, the fac_min will be zero, otherwise - ! we make use of a baseload even for SNAP2 - - tfac = ( fac_min(iland,isec,iem) + & ! constant baseload - ( 1.0-fac_min(iland,isec,iem) )* gridfac_HDD(i,j) ) & - * fac_ehh24x7(isec,hour_iland,wday_loc) - - if ( debug_tfac .and. indate%hour == 12 .and. iem==1 ) then ! - write(*,"(a,2i3,2i4,7f8.3)") "SNAPHDD tfac ", & - isec, iland, daynumber, indate%hour, & - timefac(iland_timefac,isec,iem), t2_nwp(i,j,2)-273.15, & - fac_min(iland,isec,iem), gridfac_HDD(i,j), tfac - end if - end if ! =============== HDD - - s = tfac * snapemis(isec,i,j,icc,iem) - - ! prelim emis sum kg/m2/s - SumSnapEmis(i,j,iem) = SumSnapEmis(i,j,iem) + s - - do f = 1, emis_nsplit( iem ) - iqrc = iqrc + 1 - itot = iqrc2itot(iqrc) - tmpemis(iqrc) = s * emisfrac(iqrc,isec,iland) - ! Add up emissions in ktonne - totemadd(itot) = totemadd(itot) + & - tmpemis(iqrc) * dtgrid * xmd(i,j) - end do ! f - - end do ! iem - - - ! Assign to height levels 1-KEMISTOP - - !do k=KEMISTOP,KMAX_MID - do k=KEMISTOP,KMAX_MID - do iqrc =1, nrcemis - gridrcemis0(iqrc,k,i,j) = & - gridrcemis0(iqrc,k,i,j) + tmpemis(iqrc)* & - ehlpcom0*emis_kprofile(KMAX_BND-k,isec) & - !ehlpcom0(k)*VERTFAC(KMAX_BND-k,isec) & - * emis_masscorr(iqrc) - !if( debug_tfac.and. iqrc==1 ) then - ! write(*,"(a,2i3,2f8.3)") "KPROF ", & - ! isec, KMAX_BND-k, & - ! VERTFAC(KMAX_BND-k,isec), & - ! emis_kprofile(KMAX_BND-k,isec) - !end if - - end do ! iem - end do ! k + ehlpcom0 = GRAV* 0.001*AVOG!/(sigma_bnd(k+1)-sigma_bnd(k)) + +! Scaling for totemadd: + dtgrid = dt_advec * GRIDWIDTH_M * GRIDWIDTH_M + +! The emis array only needs to be updated every full hour. The +! time-factor calculation needs to know if a local-time causes a shift +! from day to night. In addition, we reset an overall day's time-factors +! at midnight every day. + + hourchange = (indate%hour/=oldhour).or.(indate%day/=oldday) + if(hourchange) then + oldhour = indate%hour + if(indate%day/=oldday)then + !========================== + call NewDayFactors(indate) + if(USE_DEGREEDAY_FACTORS) call DegreeDayFactors(daynumber) ! => fac_emm, fac_edd + !========================== + ! for ROADDUST + wday=day_of_week(indate%year,indate%month,indate%day) + if(wday==0)wday=7 ! Sunday -> 7 + oldday = indate%day + endif + endif + + if(DEBUG_EMISTIMEFACS.and.MasterProc) & + write(*,"(a,2f8.3)") " EmisSet traffic 24x7", & + fac_ehh24x7(ISNAP_TRAF,1,4),fac_ehh24x7(ISNAP_TRAF,13,4) +!.......................................... +! Look for day-night changes, after local time correction +! (daytime(iland) not used if LONGITUDE_TIME=true) + do iland = 1, NLAND + daytime(iland) = 0 + hourloc = indate%hour + Country(iland)%timezone + localhour(iland) = hourloc ! here from 0 to 23 + if(hourloc>=7 .and. hourloc<=18) daytime(iland)=1 + enddo ! iland + + if(hourchange) then + totemadd(:) = 0. + gridrcemis0(:,:,:,:) = 0.0 + SumSnapEmis(:,:,:) = 0.0 + if(USE_ROADDUST)gridrcroadd0(:,:,:) = 0.0 +!.......................................... +! Process each grid: + do j = 1,ljmax + do i = 1,limax + ncc = nlandcode(i,j) ! No. of countries in grid + debug_tfac=(DEBUG_EMISTIMEFACS.and.debug_proc.and.i==DEBUG_li.and.j==DEBUG_lj) + ! find the approximate local time: + hourloc= mod(nint(indate%hour+24*(1+glon(i,j)/360.0)),24) + hour_longitude=hourloc + daytime_longitude=0 + if(hourloc>=7 .and. hourloc<= 18) daytime_longitude=1 + !************************************************* + ! First loop over non-flat (one sector) emissions + !************************************************* + tmpemis(:)=0. + do icc = 1, ncc + iland = landcode(i,j,icc) ! 1=Albania, etc. + iland_timefac = Country(iland)%timefac_index + if(Country(iland)%timezone==-100)then + daytime_iland=daytime_longitude + hour_iland=hour_longitude + 1 ! add 1 to get 1..24 + else + daytime_iland=daytime(iland) + hour_iland=localhour(iland) + 1 + endif + !if( hour_iland > 24 ) hour_iland = 1 !DSA12 + wday_loc=wday + if(hour_iland>24) then + hour_iland = hour_iland - 24 + wday_loc=wday + 1 + if(wday_loc==0)wday_loc=7 ! Sunday -> 7 + if(wday_loc>7 )wday_loc=1 + endif + call CheckStop(hour_iland<1,"ERROR: HOUR Zero in EmisSet") + if(debug_tfac) then + write(*,"(a,i4,2i3,i5,2i4,3x,4i3)") "EmisSet DAYS times ", daynumber, & + wday, wday_loc, iland, daytime_longitude, daytime_iland,& + hour_longitude, hour_iland, hourloc, Country(iland)%timezone + call datewrite("EmisSet DAY 24x7:", & + (/ icc, iland, wday, wday_loc, hour_iland /), & + (/ fac_ehh24x7(ISNAP_TRAF,hour_iland,wday_loc) /) ) + endif +! As each emission sector has a different diurnal profile +! and possibly speciation, we loop over each sector, adding +! the found emission rates to gridrcemis as we go. +! ================================================== + do isec = 1, NSECTORS ! Loop over snap codes + ! Calculate emission rates from snapemis, time-factors, + ! and if appropriate any speciation fraction (NEMIS_FRAC) + iqrc = 0 ! index over emisfrac + do iem = 1, NEMIS_FILE + tfac = timefac(iland_timefac,isec,iem) & + * fac_ehh24x7(isec,hour_iland,wday_loc) - enddo ! isec + if(debug_tfac.and.iem==1) & + write(*,"(a,2i4,f8.3)")"EmisSet DAY TFAC:",isec,hour_iland,tfac + + !Degree days - only SNAP-2 + if(USE_DEGREEDAY_FACTORS .and. & + isec==ISNAP_DOM .and. Gridded_SNAP2_Factors) then + oldtfac = tfac + ! If INERIS_SNAP2 set, the fac_min will be zero, otherwise + ! we make use of a baseload even for SNAP2 + tfac = ( fac_min(iland,isec,iem) & ! constant baseload + + ( 1.0-fac_min(iland,isec,iem) )* gridfac_HDD(i,j) ) & + * fac_ehh24x7(isec,hour_iland,wday_loc) + + if(debug_tfac .and. indate%hour==12 .and. iem==1) & + write(*,"(a,2i3,2i4,7f8.3)") "SNAPHDD tfac ", & + isec, iland, daynumber, indate%hour, & + timefac(iland_timefac,isec,iem), t2_nwp(i,j,2)-273.15, & + fac_min(iland,isec,iem), gridfac_HDD(i,j), tfac + endif ! =============== HDD + + s = tfac * snapemis(isec,i,j,icc,iem) + ! prelim emis sum kg/m2/s + SumSnapEmis(i,j,iem) = SumSnapEmis(i,j,iem) + s + do f = 1,emis_nsplit(iem) + iqrc = iqrc + 1 + itot = iqrc2itot(iqrc) + tmpemis(iqrc) = s * emisfrac(iqrc,isec,iland) + ! Add up emissions in ktonne + totemadd(itot) = totemadd(itot) & + + tmpemis(iqrc) * dtgrid * xmd(i,j) + enddo ! f + enddo ! iem + + ! Assign to height levels 1-KEMISTOP + do k=KEMISTOP,KMAX_MID + do iqrc =1, nrcemis + gridrcemis0(iqrc,k,i,j) = gridrcemis0(iqrc,k,i,j) & + + tmpemis(iqrc)*ehlpcom0 & + *emis_kprofile(KMAX_BND-k,isec) & + !*ehlpcom0(k)*VERTFAC(KMAX_BND-k,isec) & + *emis_masscorr(iqrc) + !if( debug_tfac.and. iqrc==1 ) then + ! write(*,"(a,2i3,2f8.3)") "KPROF ", & + ! isec, KMAX_BND-k, & + ! VERTFAC(KMAX_BND-k,isec), & + ! emis_kprofile(KMAX_BND-k,isec) + !end if + enddo ! iem + enddo ! k + enddo ! isec ! ================================================== - - end do ! icc - - !************************************ - ! Then loop over flat emissions - !************************************ - tmpemis(:)=0. - fncc = flat_nlandcode(i,j) ! No. of countries with flat - ! emissions in grid - - do ficc = 1, fncc + enddo ! icc +!************************************ +! Then loop over flat emissions +!************************************ + tmpemis(:)=0. + fncc = flat_nlandcode(i,j) ! No. of countries with flat emissions in grid + do ficc = 1, fncc flat_iland = flat_landcode(i,j,ficc) ! 30=BAS etc. - - if ( Country(flat_iland)%is_sea ) then ! saves if statements below - isec = ISNAP_SHIP + if(Country(flat_iland)%is_sea) then ! saves if statements below + isec = ISNAP_SHIP else - isec = ISNAP_NAT - end if - + isec = ISNAP_NAT + endif ! As each emission sector has a different diurnal profile ! and possibly speciation, we loop over each sector, adding ! the found emission rates to gridrcemis as we go. ! ================================================== - - ! Calculate emission rates from snapemis, time-factors, ! and if appropriate any speciation fraction (NEMIS_FRAC) - - iqrc = 0 ! index over emis - - do iem = 1, NEMIS_FILE - - sf = snapemis_flat(i,j,ficc,iem) - + iqrc = 0 ! index over emis + do iem = 1, NEMIS_FILE + sf = snapemis_flat(i,j,ficc,iem) ! prelim emis sum kg/m2/s - SumSnapEmis(i,j,iem) = SumSnapEmis(i,j,iem) + sf - - do f = 1, emis_nsplit( iem ) - iqrc = iqrc + 1 - itot = iqrc2itot(iqrc) - tmpemis(iqrc) = sf * emisfrac(iqrc,isec,flat_iland) - - ! Add flat emissions in ktonne - totemadd(itot) = totemadd(itot) + & - tmpemis(iqrc) * dtgrid * xmd(i,j) - - end do ! f - - end do ! iem - - ! Assign flat emissions to height levels 1-4 - ! Note, no VERTFAC - - do iqrc =1, nrcemis - - gridrcemis0(iqrc,KMAX_MID,i,j) = & - gridrcemis0(iqrc,KMAX_MID,i,j) + tmpemis(iqrc)*& - ehlpcom0 * emis_masscorr(iqrc) - end do ! iem - + SumSnapEmis(i,j,iem) = SumSnapEmis(i,j,iem) + sf + do f = 1,emis_nsplit(iem) + iqrc = iqrc + 1 + itot = iqrc2itot(iqrc) + tmpemis(iqrc) = sf * emisfrac(iqrc,isec,flat_iland) + ! Add flat emissions in ktonne + totemadd(itot) = totemadd(itot) & + + tmpemis(iqrc) * dtgrid * xmd(i,j) + enddo ! f + enddo ! iem + ! Assign flat emissions to height levels 1-4. Note, no VERTFAC + do iqrc =1, nrcemis + gridrcemis0(iqrc,KMAX_MID,i,j) = gridrcemis0(iqrc,KMAX_MID,i,j) & + + tmpemis(iqrc)*ehlpcom0*emis_masscorr(iqrc) + enddo ! iem ! ================================================== - end do !ficc - - if(USE_ROADDUST)then - NO_PRECIP: if(surface_precip(i,j) < 0.1) then ! Limit as in TNO-model (but Lotos/Euros has precip in mm/3h) In the EMEP case this is in mm/h, so should be equivalent with 2.4mm per day + enddo !ficc + if(USE_ROADDUST)then + ! Limit as in TNO-model (but Lotos/Euros has precip in mm/3h) + ! In the EMEP case this is in mm/h, so should be equivalent with 2.4mm per day + NO_PRECIP: if(surface_precip(i,j) < 0.1) then + ! should use the temporal variation for road dust (SNAP_HOURFAC(HH,7)) ! and a weekday factor (initially taken from TNO model, could use country ! dependent factor in the future) @@ -1321,174 +1107,146 @@ subroutine EmisSet(indate) ! emission re-set every time-step/hour ! weekday and diurnal variation (same for all countries) ! -> Need to know day_of_week ! Relatively weak variation with day of week so use a simplified approach -! - -! if( DEBUG_ROADDUST .and. debug_proc .and. i==DEBUG_li .and. j==DEBUG_lj )THEN -! write(*,*)"DEBUG ROADDUST! Dry! ncc=", road_nlandcode(i,j) -! endif - - ncc = road_nlandcode(i,j) ! number of countries in grid point - do icc = 1, ncc - iland = road_landcode(i,j,icc) - - if(Country(iland)%timezone==-100)then - hour_iland=hour_longitude+1 - else - hour_iland=localhour(iland)+1 - endif - wday_loc = wday ! DS added here also, for fac_ehh24x7 - if( hour_iland > 24 ) then + ! if( DEBUG_ROADDUST .and. debug_proc .and. i==DEBUG_li .and. j==DEBUG_lj )THEN + ! write(*,*)"DEBUG ROADDUST! Dry! ncc=", road_nlandcode(i,j) + ! endif + + ncc = road_nlandcode(i,j) ! number of countries in grid point + do icc = 1, ncc + iland = road_landcode(i,j,icc) + if(Country(iland)%timezone==-100)then + hour_iland=hour_longitude+1 + else + hour_iland=localhour(iland)+1 + endif + + wday_loc = wday ! DS added here also, for fac_ehh24x7 + if( hour_iland > 24 ) then hour_iland = 1 if(wday_loc==0)wday_loc=7 ! Sunday -> 7 if(wday_loc>7 )wday_loc=1 - end if + endif - if(((icc.eq.IC_FI).or.(icc.eq.IC_NO).or.(icc.eq.IC_SE)).and. & ! Nordic countries - ((indate%month.eq.3).or.(indate%month.eq.4).or.(indate%month.eq.5)))then ! spring road dust - tfac = fac_ehh24x7(ISNAP_TRAF,hour_iland,wday_loc) *2.0 ! Doubling in Mar-May (as in TNO model) - else - tfac = fac_ehh24x7(ISNAP_TRAF,hour_iland,wday_loc) - endif + if(ANY(icc==(/IC_FI,IC_NO,IC_SE/)).and. & ! Nordic countries + ANY(indate%month==(/3,4,5/)))then ! spring road dust + tfac = fac_ehh24x7(ISNAP_TRAF,hour_iland,wday_loc)*2.0 ! Doubling in Mar-May (as in TNO model) + else + tfac = fac_ehh24x7(ISNAP_TRAF,hour_iland,wday_loc) + endif - do iem = 1, NROAD_FILES + do iem = 1, NROAD_FILES s = tfac * roaddust_emis_pot(i,j,icc,iem) - if( DEBUG_ROADDUST .and. debug_proc .and. i==DEBUG_li .and. j==DEBUG_lj )THEN - write(*,*)"DEBUG ROADDUST! iem,tfac,icc,roaddust_emis_pot,s", & - iem,tfac,icc,roaddust_emis_pot(i,j,icc,iem),s - endif - - gridrcroadd0(QROADDUST_FI,i,j)=gridrcroadd0(QROADDUST_FI,i,j)+ & - ROADDUST_FINE_FRAC*s - gridrcroadd0(QROADDUST_CO,i,j)=gridrcroadd0(QROADDUST_CO,i,j)+ & - (1.-ROADDUST_FINE_FRAC)*s - - if ( DEBUG_ROADDUST .AND. debug_proc .and.i==debug_li .and. j==debug_lj) then ! - write(*,*)"gridrcroadfine",gridrcroadd0(QROADDUST_FI,i,j) - write(*,*)"gridrcroadcoarse",gridrcroadd0(QROADDUST_CO,i,j) - end if - - - enddo ! nroad files - - enddo ! icc - + if(DEBUG_ROADDUST.and.debug_proc.and.i==DEBUG_li.and.j==DEBUG_lj)& + write(*,*)"DEBUG ROADDUST! iem,tfac,icc,roaddust_emis_pot,s", & + iem,tfac,icc,roaddust_emis_pot(i,j,icc,iem),s + + gridrcroadd0(QROADDUST_FI,i,j)=gridrcroadd0(QROADDUST_FI,i,j) & + +ROADDUST_FINE_FRAC*s + gridrcroadd0(QROADDUST_CO,i,j)=gridrcroadd0(QROADDUST_CO,i,j) & + +(1.-ROADDUST_FINE_FRAC)*s + + if(DEBUG_ROADDUST.AND.debug_proc.AND.& + i==debug_li.AND.j==debug_lj)then ! + write(*,*)"gridrcroadfine" ,gridrcroadd0(QROADDUST_FI,i,j) + write(*,*)"gridrcroadcoarse",gridrcroadd0(QROADDUST_CO,i,j) + endif + enddo ! nroad files + enddo ! icc ! should pick the correct emissions (spring or rest of year) ! and add the emissions from HIGHWAYplus and NONHIGHWAYS, ! using correct fine and coarse fractions. - else ! precipitation - gridrcroadd0(:,i,j)=0. - endif NO_PRECIP - - endif ! ROADDUST - - end do ! i - end do ! j - if ( DEBUG .and. debug_proc ) then ! emis sum kg/m2/s - call datewrite("SnapSum, kg/m2/s:" // trim(EMIS_FILE(iemCO)), & - (/ SumSnapEmis(debug_li,debug_lj,iemCO) /) ) - end if - - call Set_Volc !set hourly volcano emission(rcemis_volc0) - - - end if ! hourchange - + else ! precipitation + gridrcroadd0(:,i,j)=0. + endif NO_PRECIP + endif ! ROADDUST + enddo ! i + enddo ! j + if(DEBUG.and.debug_proc) & ! emis sum kg/m2/s + call datewrite("SnapSum, kg/m2/s:"//trim(EMIS_FILE(iemCO)), & + (/ SumSnapEmis(debug_li,debug_lj,iemCO) /) ) + + call Set_Volc !set hourly volcano emission(rcemis_volc0) + endif ! hourchange ! We now scale gridrcemis to get emissions in molecules/cm3/s - - do k= KEMISTOP, KMAX_MID - do j = 1,ljmax - do i = 1,limax - - ehlpcom= roa(i,j,k,1)/(dA(k)+dB(k)*ps(i,j,1)) + do k= KEMISTOP, KMAX_MID + do j = 1,ljmax + do i = 1,limax + ehlpcom= roa(i,j,k,1)/(dA(k)+dB(k)*ps(i,j,1)) !RB: This should also be done for the road dust emissions - do iqrc =1, NRCEMIS - gridrcemis(iqrc,k,i,j) = & - gridrcemis0(iqrc,k,i,j)* ehlpcom - enddo ! iqrc - end do ! i - end do ! j - end do ! k - - if(USE_ROADDUST)THEN - if( DEBUG_ROADDUST .and. debug_proc)then - write(*,*)"Before the unit scaling", & - gridrcroadd(1,DEBUG_li,DEBUG_lj), & - gridrcroadd(2,DEBUG_li,DEBUG_lj) - endif - do j = 1,ljmax - do i = 1,limax - ehlpcom= roa(i,j,KMAX_MID,1)/(ps(i,j,1)-PT) - do iqrc =1, NROADDUST - gridrcroadd(iqrc,i,j) = & - gridrcroadd0(iqrc,i,j)* ehlpcom * ehlpcom0 * roaddust_masscorr(iqrc) - enddo ! iqrc - end do ! i - end do ! j - if( DEBUG_ROADDUST .and. debug_proc)then - write(*,*)"After the unit scaling", & - gridrcroadd(1,DEBUG_li,DEBUG_lj), & - gridrcroadd(2,DEBUG_li,DEBUG_lj) - endif - endif + do iqrc =1, NRCEMIS + gridrcemis(iqrc,k,i,j) = gridrcemis0(iqrc,k,i,j)* ehlpcom + enddo ! iqrc + enddo ! i + enddo ! j + enddo ! k + + if(USE_ROADDUST)THEN + if(DEBUG_ROADDUST.and.debug_proc) & + write(*,*)"Before the unit scaling",gridrcroadd(1:2,DEBUG_li,DEBUG_lj) + do j = 1,ljmax + do i = 1,limax + ehlpcom= roa(i,j,KMAX_MID,1)/(ps(i,j,1)-PT) + do iqrc =1, NROADDUST + gridrcroadd(iqrc,i,j) = gridrcroadd0(iqrc,i,j)* ehlpcom * ehlpcom0 & + * roaddust_masscorr(iqrc) + enddo ! iqrc + enddo ! i + enddo ! j + if(DEBUG_ROADDUST.and.debug_proc) & + write(*,*)"After the unit scaling",gridrcroadd(1:2,DEBUG_li,DEBUG_lj) + endif ! Scale volc emissions to get emissions in molecules/cm3/s (rcemis_volc) - - call Scale_Volc - - end subroutine EmisSet - !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< - - subroutine newmonth - -!..................................................................... -! DESCRIPTION: + call Scale_Volc +endsubroutine EmisSet +!*********************************************************************** +subroutine newmonth +!----------------------------------------------------------------------! +! DESCRIPTION: ! Reads in natural DMS emissions at start of each month. Update ! landcode and nlandcode arrays as needed. ! ! Reads in snow cover at start of each month. ! ! April 2010: read monthly aircraft NOx emissions -! -!........................................................................... - use AirEmis_ml, only : airn - use ModelConstants_ml, only : KCHEMTOP, KMAX_MID - use NetCDF_ml, only : ReadField_CDF - - integer i, j,k, iyr - integer n, flat_ncmaxfound ! Max. no. countries w/flat emissions - real :: rdemis(MAXLIMAX,MAXLJMAX) ! Emissions read from file - character(len=20) :: fname - real ktonne_to_kgm2s, tonnemonth_to_kgm2s ! Units conversion - integer :: IQSO2 ! Index of sox in EMIS_FILE - integer errcode - real, allocatable, dimension(:,:,:,:) :: globemis - integer :: month,iem,ic,iic,isec, err3,icc - real :: duml,dumh,tmpsec(NSECTORS),conv - logical , save :: first_call=.true. - real, dimension(NSECTORS,MAXLIMAX,MAXLJMAX,NCMAX,NEMIS_FILE) & - :: snapemis_month ! monthly emissions tonne/month - logical :: needed_found - - ! For now, only the global runs use the Monthly files - - !NML logical, parameter :: MONTHLY_GRIDEMIS= IS_GLOBAL - integer :: kstart,kend,nstart,Nyears - real :: buffer(MAXLIMAX,MAXLJMAX),SumSoilNOx +!----------------------------------------------------------------------! + use AirEmis_ml, only : airn + use ModelConstants_ml, only : KCHEMTOP, KMAX_MID + use NetCDF_ml, only : ReadField_CDF + + integer i, j,k, iyr + integer n, flat_ncmaxfound ! Max. no. countries w/flat emissions + real :: rdemis(MAXLIMAX,MAXLJMAX) ! Emissions read from file + character(len=20) :: fname + real ktonne_to_kgm2s, tonnemonth_to_kgm2s ! Units conversion + integer :: IQSO2 ! Index of sox in EMIS_FILE + integer errcode + real, allocatable, dimension(:,:,:,:) :: globemis + integer :: month,iem,ic,iic,isec, err3,icc + real :: duml,dumh,tmpsec(NSECTORS),conv + logical , save :: first_call=.true. + real, dimension(NSECTORS,MAXLIMAX,MAXLJMAX,NCMAX,NEMIS_FILE) :: & + snapemis_month ! monthly emissions tonne/month + logical :: needed_found + + ! For now, only the global runs use the Monthly files + !NML logical, parameter :: MONTHLY_GRIDEMIS= IS_GLOBAL + integer :: kstart,kend,nstart,Nyears + real :: buffer(MAXLIMAX,MAXLJMAX),SumSoilNOx - if(.not.allocated(airn).and.(USE_LIGHTNING_EMIS.or.USE_AIRCRAFT_EMIS))then - allocate(airn(KCHEMTOP:KMAX_MID,MAXLIMAX,MAXLJMAX)) - endif + if(.not.allocated(airn).and.(USE_LIGHTNING_EMIS.or.USE_AIRCRAFT_EMIS))& + allocate(airn(KCHEMTOP:KMAX_MID,MAXLIMAX,MAXLJMAX)) -if( USE_AIRCRAFT_EMIS )then - airn = 0.0 - kstart=KCHEMTOP - kend=KMAX_MID + if(USE_AIRCRAFT_EMIS)then + airn = 0.0 + kstart=KCHEMTOP + kend=KMAX_MID - call ReadField_CDF('AircraftEmis_FL.nc','NOx',airn,& - nstart=current_date%month,kstart=kstart,kend=kend,& - interpol='mass_conservative', needed=.true.,debug_flag=.false.) + call ReadField_CDF('AircraftEmis_FL.nc','NOx',airn,& + nstart=current_date%month,kstart=kstart,kend=kend,& + interpol='mass_conservative', needed=.true.,debug_flag=.false.) ! convert from kg(NO2)/month into molecules/cm3/s ! from kg to molecules: 0.001*AVOG/species(NO2)%molwt @@ -1498,112 +1256,108 @@ subroutine newmonth ! dV=dz*dx*dy=dz*gridwidth**2/xm**2 *1e6 (1e6 for m3->cm3) ! from month to seconds: ndaysmonth*24*3600 -conv=0.001*AVOG/species(NO2)%molwt*GRAV/gridwidth_m**2*1.0e-6/(nmdays(current_date%month)*24*3600) - -do k=KCHEMTOP,KMAX_MID -do j=1,ljmax -do i=1,limax + conv=0.001*AVOG/species(NO2)%molwt*GRAV/gridwidth_m**2*1.0e-6& + /(nmdays(current_date%month)*24*3600) -airn(k,i,j)=airn(k,i,j)*conv*(roa(i,j,k,1))/(dA(k) + dB(k)*ps(i,j,1))*xm2(i,j) - -enddo -enddo -enddo + do k=KCHEMTOP,KMAX_MID + do j=1,ljmax + do i=1,limax + airn(k,i,j)=airn(k,i,j)*conv*(roa(i,j,k,1))& + /(dA(k)+dB(k)*ps(i,j,1))*xm2(i,j) -endif + enddo + enddo + enddo + endif -if(USE_EURO_SOILNOX)then ! European Soil NOx emissions - if (DEBUG_SOILNOX .and. debug_proc ) write(*,*) "Emissions DEBUG_SOILNOX START" + if(USE_EURO_SOILNOX)then ! European Soil NOx emissions + if(DEBUG_SOILNOX.and.debug_proc) write(*,*)"Emissions DEBUG_SOILNOX START" - ! read in map of annual N-deposition produced from pre-runs of EMEP model - ! with script mkcdo.annualNdep - ! - call ReadField_CDF('annualNdep.nc',& - 'Ndep_m2',AnnualNdep,1, interpol='zero_order',needed=.true.,debug_flag=.false.,UnDef=0.0) + ! read in map of annual N-deposition produced from pre-runs of EMEP model + ! with script mkcdo.annualNdep + call ReadField_CDF('annualNdep.nc','Ndep_m2',AnnualNdep,1,& + interpol='zero_order',needed=.true.,debug_flag=.false.,UnDef=0.0) - if (DEBUG_SOILNOX .and. debug_proc ) then + if(DEBUG_SOILNOX.and.debug_proc) then write(*,"(a,4es12.3)") "Emissions_ml: SOILNOX AnnualDEBUG ", & AnnualNdep(debug_li, debug_lj), maxval(AnnualNdep), minval(AnnualNdep) - end if - call CheckStop(USE_GLOBAL_SOILNOX, "SOILNOX - cannot use global with Euro") - ! We then calculate SoulNOx in Biogenics_ml - -else if ( USE_GLOBAL_SOILNOX ) then ! Global soil NOx - - do j=1,ljmax - do i=1,limax - SoilNOx(i,j)=0.0 - buffer(i,j)=0.0 - enddo - enddo - - nstart=(current_date%year-1996)*12 + current_date%month - if(nstart>0.and.nstart<=120)then - !the month is defined - call ReadField_CDF('nox_emission_1996-2005.nc','NOX_EMISSION',SoilNOx,nstart=nstart,& - interpol='conservative',known_projection="lon lat",needed=.true.,debug_flag=.false.) - if ( DEBUG_SOILNOX .and.debug_proc ) write(*,*) "PROPER YEAR of SOILNO ", current_date%year, nstart - else - !the year is not defined; average over all years - Nyears=10 !10 years defined - do iyr=1,Nyears - nstart=12*(iyr-1) + current_date%month - call ReadField_CDF('nox_emission_1996-2005.nc','NOX_EMISSION',buffer,nstart=nstart,& - interpol='conservative',known_projection="lon lat",needed=.true.,debug_flag=.false.,UnDef=0.0) - do j=1,ljmax - do i=1,limax - SoilNOx(i,j)=SoilNOx(i,j)+buffer(i,j) - end do - end do - if ( DEBUG_SOILNOX .and.debug_proc ) then - write(*,"(a,2i6,es10.3,a,2es10.3)") "Averaging SOILNO inputs", & - 1995+(iyr-1), nstart,SoilNOx(debug_li, debug_lj), & + endif + call CheckStop(USE_GLOBAL_SOILNOX, "SOILNOX - cannot use global with Euro") + ! We then calculate SoulNOx in Biogenics_ml + + elseif(USE_GLOBAL_SOILNOX) then ! Global soil NOx + + SoilNOx(:,:)=0.0 + buffer(:,:)=0.0 + nstart=(current_date%year-1996)*12 + current_date%month + if(nstart>0.and.nstart<=120)then + !the month is defined + call ReadField_CDF('nox_emission_1996-2005.nc','NOX_EMISSION',SoilNOx,& + !nstart=nstart,interpol='conservative',known_projection="lon lat",needed=.true.,debug_flag=.false.) + nstart=nstart,interpol='conservative',known_projection="lon lat",& + needed=.true.,debug_flag=.false.,UnDef=0.0) + if(DEBUG_SOILNOX.and.debug_proc) & + write(*,*) "PROPER YEAR of SOILNO ", current_date%year, nstart + else + !the year is not defined; average over all years + Nyears=10 !10 years defined + do iyr=1,Nyears + nstart=12*(iyr-1) + current_date%month + call ReadField_CDF('nox_emission_1996-2005.nc','NOX_EMISSION',buffer,& + nstart=nstart,interpol='conservative',known_projection="lon lat",& + needed=.true.,debug_flag=.false.,UnDef=0.0) + do j=1,ljmax + do i=1,limax + SoilNOx(i,j)=SoilNOx(i,j)+buffer(i,j) + enddo + enddo + if(DEBUG_SOILNOX.and.debug_proc) & + write(*,"(a,2i6,es10.3,a,2es10.3)") "Averaging SOILNO inputs", & + 1995+(iyr-1), nstart,SoilNOx(debug_li, debug_lj), & "max: ", maxval(buffer), maxval(SoilNOx) - end if - enddo - SoilNOx=SoilNOx/Nyears - endif ! nstart test - - if ( DEBUG_SOILNOX .and. debug_proc ) then - write(*,"(a,i3,3es10.3)") "After Global SOILNO ", me, maxval(SoilNOx), & - SoilNOx(debug_li, debug_lj) + enddo + SoilNOx=SoilNOx/Nyears + endif ! nstart test + + if(DEBUG_SOILNOX.and.debug_proc) then + write(*,"(a,i3,3es10.3)") "After Global SOILNO ",& + me,maxval(SoilNOx),SoilNOx(debug_li,debug_lj) !write(*,"(a,i3,3es10.3)") "After Global SOILNO ", me, maxval(SoilNOx), SoilNOx(3, 3) - end if -else ! no soil NO - if (DEBUG_SOILNOX .and. debug_proc ) write(*,*) "Emissions DEBUG_SOILNOX - none" -end if ! SOIL NO + endif + else ! no soil NO + if(DEBUG_SOILNOX.and.debug_proc) & + write(*,*) "Emissions DEBUG_SOILNOX - none" + endif ! SOIL NO !for testing, compute total soil NOx emissions within domain !convert from g/m2/day into kg/day -if ( USE_GLOBAL_SOILNOX ) then - SumSoilNOx=0.0 - SoilNOx = max(0.0, SoilNOx) ! Stops the NEGs! - do j=1,ljmax - do i=1,limax - SumSoilNOx=SumSoilNOx+0.001*SoilNOx(i,j)*gridwidth_m**2*xmd(i,j) - enddo - enddo - CALL MPI_ALLREDUCE(MPI_IN_PLACE, SumSoilNOx , 1, & - MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_WORLD, INFO) - if(MasterProc)write(*,*)'GLOBAL SOILNOX emissions this month within domain',& - SumSoilNOx,' kg per day' - - ! convert from g(N)/m2/day into molecules/cm3/s from g to molecules: - ! AVOG/14 14=molweight N, use roa to find dz for consistency with other - ! emissions (otherwise could have used z_bnd directly) dz=dP/(roa*GRAV) - ! dP=dA(k) + dB(k)*ps(i,j,1) dV=dz*1e6 (1e6 for m3->cm3) from month to - ! seconds: ndaysmonth*24*3600 - - conv=AVOG/14.0*GRAV*1.0e-6/(24*3600) - k=KMAX_MID!surface - do j=1,ljmax - do i=1,limax - SoilNOx(i,j)=SoilNOx(i,j)*conv*(roa(i,j,k,1))/(dA(k) + dB(k)*ps(i,j,1)) - enddo - enddo - -end if - + if(USE_GLOBAL_SOILNOX) then + SumSoilNOx=0.0 + SoilNOx = max(0.0, SoilNOx) ! Stops the NEGs! + do j=1,ljmax + do i=1,limax + SumSoilNOx=SumSoilNOx+0.001*SoilNOx(i,j)*gridwidth_m**2*xmd(i,j) + enddo + enddo + CALL MPI_ALLREDUCE(MPI_IN_PLACE,SumSoilNOx,1,MPI_DOUBLE_PRECISION, & + MPI_SUM,MPI_COMM_WORLD,INFO) + if(MasterProc)& + write(*,*)'GLOBAL SOILNOX emissions this month within domain',& + SumSoilNOx,' kg per day' + +! convert from g(N)/m2/day into molecules/cm3/s from g to molecules: +! AVOG/14 14=molweight N, use roa to find dz for consistency with other +! emissions (otherwise could have used z_bnd directly) dz=dP/(roa*GRAV) +! dP=dA(k) + dB(k)*ps(i,j,1) dV=dz*1e6 (1e6 for m3->cm3) from month to +! seconds: ndaysmonth*24*3600 + conv=AVOG/14.0*GRAV*1.0e-6/(24*3600) + k=KMAX_MID!surface + do j=1,ljmax + do i=1,limax + SoilNOx(i,j)=SoilNOx(i,j)*conv*roa(i,j,k,1)/(dA(k)+dB(k)*ps(i,j,1)) + enddo + enddo + endif ! DMS ! Units: @@ -1617,290 +1371,235 @@ subroutine newmonth ! than the sum). ! More precisely: year_emis=sum_months(emis_month*nmdays/nydays) - ktonne_to_kgm2s = 1.0e6 / & - (nydays*24.*60.*60.*GRIDWIDTH_M*GRIDWIDTH_M) - - if ( MasterProc .and. DEBUG) then - write(6,*) 'Enters newmonth, mm, ktonne_to_kgm2s = ', & - current_date%month,ktonne_to_kgm2s - write(6,*) ' first_dms_read = ', first_dms_read - end if ! me -!........................................................................... + ktonne_to_kgm2s = 1.0e6/(nydays*24.*60.*60.*GRIDWIDTH_M*GRIDWIDTH_M) + if(MasterProc.and.DEBUG) then + write(*,*) 'Enters newmonth, mm, ktonne_to_kgm2s = ', & + current_date%month,ktonne_to_kgm2s + write(*,*) ' first_dms_read = ', first_dms_read + endif !........................................................................... ! DMS Input - land 35 - SNAP sector 11 !........................................................................... - flat_ncmaxfound = 0 ! Max. no. countries(w/flat emissions) per grid - + flat_ncmaxfound = 0 ! Max. no. countries(w/flat emissions) per grid ! Natural SO2 emissions + IQSO2=find_index("sox",EMIS_FILE(:)) + if(IQSO2<1) then + write(*,*) " No SO2 emissions - need to skip DMS also" + return ! No need to read DMS fields + else + ! We have so2 emission so need DMS also + if(MasterProc) then + write(fname,'(''natso2'',i2.2,''.dat'')')current_date%month + write(*,*) 'Reading DMS emissions from ',trim(fname) + endif + needed_found=.false. + call ReadField(IO_DMS,fname,rdemis,needed_found) + if(needed_found)then + errcode = 0 + do j=1,ljmax + do i=1,limax +! Add DMS for country code IQ_DMS=35 to snap sector 11=Nature. +! First time we read we must add DMS to the "countries" +! contributing within the grid square. + ! - for flat emissions: + if(first_dms_read) then + flat_nlandcode(i,j) = flat_nlandcode(i,j) + 1 + n = flat_nlandcode(i,j) + flat_landcode(i,j,n) = IQ_DMS ! country code 35 + if(n>flat_ncmaxfound) then + flat_ncmaxfound = n + if (DEBUG) write(6,*)'DMS Increased flat_ncmaxfound to ',n + call CheckStop( n > FNCMAX, "IncreaseFNCMAX for dms") + endif + else ! We know that DMS lies last in the array, so: + n = flat_nlandcode(i,j) + call CheckStop(flat_landcode(i,j,n),IQ_DMS,"Newmonth:DMS not last!") + endif + snapemis_flat(i,j,n,IQSO2) = rdemis(i,j) * ktonne_to_kgm2s * xm2(i,j) + enddo ! i + enddo ! j + + if(first_dms_read) then + if(DEBUG) & + write(*,*)'me ',me, ' Increased flat_ncmaxfound to ',flat_ncmaxfound + first_dms_read = .false. + endif + else!no dms file found + call PrintLog("WARNING: NO DMS emissions found",MasterProc) + endif + endif ! IQSO2>0 - IQSO2 = 0 - do i = 1, NEMIS_FILE - if ( trim( EMIS_FILE(i) ) == "sox" ) IQSO2 = i - end do - - if ( IQSO2 < 1 ) then - write(*,*) " No SO2 emissions - need to skip DMS also" - return ! No need to read DMS fields - - else - ! We have so2 emission so need DMS also - - if ( MasterProc ) then - - write(fname,fmt='(''natso2'',i2.2,''.dat'')') & - current_date%month - write(6,*) 'Reading DMS emissions from ',trim(fname) - endif - - needed_found=.false. - call ReadField(IO_DMS,fname,rdemis,needed_found) - if(needed_found)then - errcode = 0 - do j=1,ljmax - do i=1,limax - - ! Add DMS for country code IQ_DMS=35 to snap sector 11=Nature. - ! First time we read we must add DMS to the "countries" - ! contributing within the grid square. - - ! - for flat emissions: - - if ( first_dms_read ) then - flat_nlandcode(i,j) = flat_nlandcode(i,j) + 1 - n = flat_nlandcode(i,j) - flat_landcode(i,j,n) = IQ_DMS ! country code 35 - if ( n > flat_ncmaxfound ) then - flat_ncmaxfound = n - if (DEBUG) write(6,*)'DMS Increased flat_ncmaxfound to ',n - call CheckStop( n > FNCMAX, "IncreaseFNCMAX for dms") - endif - else ! We know that DMS lies last in the array, so: - n = flat_nlandcode(i,j) - call CheckStop(flat_landcode(i,j,n), IQ_DMS, & - "Newmonth:DMS not last!") - endif - - snapemis_flat(i,j,n,IQSO2) = rdemis(i,j) * ktonne_to_kgm2s & - * xm2(i,j) - enddo ! i - enddo ! j - - - if ( first_dms_read ) then - if (DEBUG) write(6,*)'me ',me, ' Increased flat_ncmaxfound to ' & - ,flat_ncmaxfound - first_dms_read = .false. - end if - else!no dms file found - if ( MasterProc ) write(6,*) 'WARNING: NO DMS emissions found ' - if ( MasterProc ) write(unit=IO_LOG,fmt=*) "WARNING: NO DMS emissions found " - end if - - end if ! IQSO2>0 - - - if(MONTHLY_GRIDEMIS)then - - !Read monthly emission files - - if(first_call)then + if(MONTHLY_GRIDEMIS)then + !Read monthly emission files + if(first_call)then do j=1,ljmax - do i=1,limax - nlandcode(i,j)=nlandcode(i,j)+1 - icc=nlandcode(i,j) - landcode(i,j,icc)=67 - enddo + do i=1,limax + nlandcode(i,j)=nlandcode(i,j)+1 + icc=nlandcode(i,j) + landcode(i,j,icc)=67 + enddo enddo first_call=.false. - endif + endif month = current_date%month - tonnemonth_to_kgm2s= 1.0e3 / & - (nmdays(month)*24.*60.*60.*GRIDWIDTH_M*GRIDWIDTH_M) - - if ( MasterProc ) then - allocate(globemis(NSECTORS,GIMAX,GJMAX,NCMAX),stat=err3) - call CheckStop(err3, "Allocation error err3 - globland") - end if - do iem = 1, NEMIS_FILE + tonnemonth_to_kgm2s = 1.0e3 & + /(nmdays(month)*24.*60.*60.*GRIDWIDTH_M*GRIDWIDTH_M) + if(MasterProc)then + allocate(globemis(NSECTORS,GIMAX,GJMAX,NCMAX),stat=err3) + call CheckStop(err3, "Allocation error err3 - globland") + endif + do iem = 1,NEMIS_FILE ! if (trim(EMIS_NAME(iem)).ne.'nox' .and. trim(EMIS_NAME(iem)).ne.'co'.and.& ! trim(EMIS_NAME(iem)).ne.'pm25'.and.& ! trim(EMIS_NAME(iem)).ne.'voc'.and.trim(EMIS_NAME(iem)).ne.'nh3'.and.trim(EMIS_NAME(iem)).ne.'sox')cycle ! ! snapemis (:,:,:,:,iem) = 0.0 !NB all previous (snap)emis set to 0 - if ( MasterProc ) then - - globemis = 0.0 - - write(fname,fmt='(''grid'',A,i2.2)') & - trim(EMIS_FILE(iem))//'.',month - write(6,*) 'filename for GLOBAL emission',fname - call open_file(IO_EMIS,"r",fname,needed=.true.) - call CheckStop( ios , "ios error: emislist" // fname ) - - READEMIS: do ! Loop over emislist files - - read(unit=IO_EMIS,fmt=*,iostat=ios) iic, i, j, duml, dumh, & - (tmpsec(isec),isec=1,NSECTORS) - if ( ios < 0 ) exit READEMIS ! End of file - call CheckStop( ios , "GetEmis ios: error on " // fname ) ! exits if ios>0 - - ic=1 ! default country - - - - i = i-IRUNBEG+1 ! for RESTRICTED domain - j = j-JRUNBEG+1 ! for RESTRICTED domain - - if ( i <= 0 .or. i > GIMAX .or. & - j <= 0 .or. j > GJMAX .or. & - ic <= 0 .or. ic > NLAND .or. & - ic == IC_NAT ) & !Excludes DMS - cycle READEMIS - globemis(:,i,j,ic) = globemis(:,i,j,ic) & !Put everything in land "1" - + tmpsec(:) - - end do READEMIS - ! - close(IO_EMIS) - ios = 0 - - endif ! MasterProc + if(MasterProc) then + globemis = 0.0 + write(fname,fmt='(''grid'',A,i2.2)')trim(EMIS_FILE(iem))//'.',month + write(*,*) 'filename for GLOBAL emission',trim(fname) + call open_file(IO_EMIS,"r",fname,needed=.true.) + call CheckStop( ios , "ios error: emislist" // fname ) + + READEMIS: do ! Loop over emislist files + read(IO_EMIS,*,iostat=ios) & + iic, i, j, duml, dumh, (tmpsec(isec),isec=1,NSECTORS) + if(ios<0) exit READEMIS ! End of file + call CheckStop( ios , "GetEmis ios: error on " // fname ) ! exits if ios>0 + + ic=1 ! default country + i = i-IRUNBEG+1 ! for RESTRICTED domain + j = j-JRUNBEG+1 ! for RESTRICTED domain + if(i<=0 .or. i>GIMAX .or. j<=0 .or. j>GJMAX .or. & + ic<=0 .or. ic> NLAND .or. & + ic==IC_NAT ) & !Excludes DMS + cycle READEMIS + globemis(:,i,j,ic) = globemis(:,i,j,ic) + tmpsec(:) !Put everything in land "1" + enddo READEMIS + close(IO_EMIS) + ios = 0 + endif ! MasterProc call global2local(globemis,snapemis_month(1,1,1,1,iem),MSG_READ1, & - NSECTORS,GIMAX,GJMAX,NCMAX,1,1) - end do ! iem = 1, NEMIS-loop + NSECTORS,GIMAX,GJMAX,NCMAX,1,1) + enddo ! iem = 1, NEMIS-loop - ic=1 - do iem = 1, NEMIS_FILE + ic=1 + do iem = 1,NEMIS_FILE ! write(*,*)'iem=',iem ! if (trim(EMIS_NAME(iem)).ne.'nox' .and. trim(EMIS_NAME(iem)).ne.'co'.and.& ! trim(EMIS_NAME(iem)).ne.'pm25'.and.& ! trim(EMIS_NAME(iem)).ne.'voc'.and.trim(EMIS_NAME(iem)).ne.'nh3'.and.trim(EMIS_NAME(iem)).ne.'sox')cycle ! conv = tonnemonth_to_kgm2s do j=1,ljmax - do i=1,limax - icc=nlandcode(i,j) !67 - do isec=1,NSECTORS - snapemis (isec,i,j,icc,iem) = & - snapemis_month (isec,i,j,ic,iem) * conv * xm2(i,j) - enddo - enddo + do i=1,limax + icc=nlandcode(i,j) !67 + do isec=1,NSECTORS + snapemis (isec,i,j,icc,iem) = snapemis_month (isec,i,j,ic,iem) & + * conv * xm2(i,j) + enddo + enddo enddo - enddo !iem - if ( MasterProc ) then - deallocate(globemis) - end if - endif - - end subroutine newmonth - - !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< - subroutine EmisOut(label, iem,nsources,sources,emis) - - ! Ascii output of emissions fields (after interpolation and femis. - ! To print out emissions, we need to get fields for each country - ! in turn, needing info from GridCodes - ! e.g. 11-SNAP data is snapemis: - ! real snapemis (NSECTORS,MAXLIMAX,MAXLJMAX,NCMAX,NEMIS_FILES) - - character(len=*) :: label - integer, intent(in) :: iem - real,intent(in), dimension(:,:,:,:):: emis ! Emission values - integer, dimension(:,:), intent(in) :: nsources ! Emission values - integer, dimension(:,:,:), intent(in) :: sources ! Emission values - real, allocatable, dimension(:,:):: lemis,gemis ! 2-D emission fields - real, allocatable, dimension(:,:,:):: locemis ! Emission values - real, allocatable, dimension(:,:,:):: globemis ! Emission values - character(len=100) :: txt - real :: low, high - integer :: msg(4), i,j, ii, jj, isec, icc, ncc, iland - - - txt = trim(label)//"."//trim(EMIS_FILE(iem)) - msg(:) = 0 - - if(DEBUG) write(*,*) "CALLED "//trim(txt), me, & - maxval(emis), maxval(nsources), maxval(sources) - - if ( MasterProc ) then - allocate(globemis(GIMAX, GJMAX,NSECTORS), stat=msg(1) ) - open(IO_TMP,file="EmisOut"//trim(txt)) - end if - - allocate(locemis(MAXLIMAX, MAXLJMAX,NSECTORS), stat=msg(2) ) - allocate(lemis(MAXLIMAX, MAXLJMAX), stat=msg(3) ) - allocate(gemis(GIMAX, GJMAX), stat=msg(4) ) - - call CheckStop( any(msg(:) /= 0),"EmisOut alloc error:"//trim(label)) - - ! TEST LOCAL2GLOBAL. Found that the local array is changed - !lemis = 100 + me - !print *, "LEMIS-IN ", me, lemis(2,2) - !call local2global( lemis(:,:), gemis(:,:), msg ) - !if( MasterProc ) then - ! do j = 1, GJMAX, 3 - ! print *, "LEMIS-UT ", j, lemis(2,j), gemis(2,j) - ! end do - !end if - -EMLAND: do iland = 1, NLAND - locemis = 0.0 -! print *, trim(txt)//" iland ", me, iland, maxval(emis(:,:,:,:)) - - !/ Collect emissions locally by country code iland - do j = 1, ljmax - do i = 1, limax - ncc = nsources(i,j) - do icc = 1, ncc - if ( sources(i,j,icc) == iland ) then - locemis(i,j,: ) = emis(:, i,j,icc) - if ( DEBUG ) & - call CheckStop( any( locemis(i,j,:) < 0.0), "NEG LOCEMIS") - end if - end do - end do - end do ! j - - ! Should never happen, but... - !call CheckStop( any( lemis < 0.0 ) , "NEG LEMIS") - + enddo !iem + if(MasterProc) deallocate(globemis) + endif +endsubroutine newmonth +!*********************************************************************** +subroutine EmisOut(label, iem,nsources,sources,emis) +!----------------------------------------------------------------------! +! Ascii output of emissions fields (after interpolation and femis. +! To print out emissions, we need to get fields for each country +! in turn, needing info from GridCodes +! e.g. 11-SNAP data is snapemis: +! real snapemis (NSECTORS,MAXLIMAX,MAXLJMAX,NCMAX,NEMIS_FILES) + character(len=*) :: label + integer, intent(in) :: iem + real,intent(in), dimension(:,:,:,:):: emis ! Emission values + integer, dimension(:,:), intent(in) :: nsources ! Emission values + integer, dimension(:,:,:), intent(in) :: sources ! Emission values + real, allocatable, dimension(:,:):: lemis,gemis ! 2-D emission fields + real, allocatable, dimension(:,:,:):: locemis ! Emission values + real, allocatable, dimension(:,:,:):: globemis ! Emission values + character(len=100) :: txt + real :: low, high + integer :: msg(4), i,j, ii, jj, isec, icc, ncc, iland + + txt = trim(label)//"."//trim(EMIS_FILE(iem)) + msg(:) = 0 + + if(DEBUG) write(*,*)"CALLED "//trim(txt),me,& + maxval(emis),maxval(nsources),maxval(sources) + + if(MasterProc) then + allocate(globemis(GIMAX, GJMAX,NSECTORS), stat=msg(1) ) + open(IO_TMP,file="EmisOut"//trim(txt)) + endif + + allocate(locemis(MAXLIMAX, MAXLJMAX,NSECTORS), stat=msg(2) ) + allocate(lemis(MAXLIMAX, MAXLJMAX), stat=msg(3) ) + allocate(gemis(GIMAX, GJMAX), stat=msg(4) ) + + call CheckStop( any(msg(:) /= 0),"EmisOut alloc error:"//trim(label)) + + ! TEST LOCAL2GLOBAL. Found that the local array is changed + !lemis = 100 + me + !print *, "LEMIS-IN ", me, lemis(2,2) + !call local2global( lemis(:,:), gemis(:,:), msg ) + !if( MasterProc ) then + ! do j = 1, GJMAX, 3 + ! print *, "LEMIS-UT ", j, lemis(2,j), gemis(2,j) + ! end do + !end if + + EMLAND: do iland = 1, NLAND + locemis = 0.0 +! print *, trim(txt)//" iland ", me, iland, maxval(emis(:,:,:,:)) + + !/ Collect emissions locally by country code iland + do j = 1, ljmax + do i = 1, limax + ncc = nsources(i,j) + do icc = 1, ncc + if(sources(i,j,icc)==iland) then + locemis(i,j,: ) = emis(:, i,j,icc) + if(DEBUG) call CheckStop(any(locemis(i,j,:)< 0.0),"NEG LOCEMIS") + endif + enddo + enddo + enddo ! j - !/ Transmit to MasterProc - !/ (Need to be careful, local2global changes local arrays. Hence lemis) + ! Should never happen, but... + !call CheckStop( any( lemis < 0.0 ) , "NEG LEMIS") + - do isec = 1, NSECTORS - lemis = locemis(:,:,isec) - if( DEBUG) write(*,*) trim(txt)//" lemis ", & - me, iland, maxval(lemis(:,:)) - call local2global( lemis, gemis, msg ) - if( MasterProc ) globemis(:,:,isec) = gemis(:,:) !! for output - end do ! isec + !/ Transmit to MasterProc + !/ (Need to be careful, local2global changes local arrays. Hence lemis) + do isec = 1, NSECTORS + lemis = locemis(:,:,isec) + if(DEBUG) write(*,*) trim(txt)//" lemis ",me,iland,maxval(lemis(:,:)) + call local2global(lemis,gemis,msg) + if(MasterProc) globemis(:,:,isec) = gemis(:,:) !! for output + enddo ! isec - if( MasterProc ) then ! output emislist-type file - - do i = 1, GIMAX - do j = 1, GJMAX - ii=i+IRUNBEG-1 - jj= j+JRUNBEG-1 - if(sum(globemis(i,j,:)) > 1.0e-10) then - high = globemis(i,j,1) - low = sum(globemis(i,j,2:NSECTORS)) - write(IO_TMP,"(i3,2i4,2x,13es10.3)") iland, ii,jj, & - low, high, (globemis(i,j,isec),isec=1,NSECTORS) - end if - end do - end do - end if ! masterProc - end do EMLAND + if(MasterProc) then ! output emislist-type file + do i = 1,GIMAX + do j = 1,GJMAX + ii=i+IRUNBEG-1 + jj=j+JRUNBEG-1 + if(sum(globemis(i,j,:))>1.0e-10)then + high = globemis(i,j,1) + low = sum(globemis(i,j,2:NSECTORS)) + write(IO_TMP,"(i3,2i4,2x,13es10.3)") iland, ii,jj, & + low, high, (globemis(i,j,isec),isec=1,NSECTORS) + endif + enddo + enddo + endif ! masterProc + enddo EMLAND - deallocate(locemis) - deallocate(gemis) - deallocate(lemis) - if(MasterProc) then - close(IO_TMP) - deallocate(globemis) - end if - - end subroutine EmisOut - - -end module Emissions_ml + deallocate(locemis,gemis,lemis) + if(MasterProc) then + close(IO_TMP) + deallocate(globemis) + endif +endsubroutine EmisOut +endmodule Emissions_ml diff --git a/ExternalBICs_ml.f90 b/ExternalBICs_ml.f90 index cf4d0f6..3dc2777 100644 --- a/ExternalBICs_ml.f90 +++ b/ExternalBICs_ml.f90 @@ -29,6 +29,9 @@ module ExternalBICs_ml character(len=80),public, save :: & filename_eta = 'EMEP_IN_BC_eta.zaxis' + +character(len=*),public, parameter :: & + ICBC_FMT="(A24,'=',A24,'*',F7.2,2L2,'=',I3)" type, public :: icbc ! Inital (IC) & Boundary Conditions (BC) character(len=24) :: spcname="none",varname="none" ! Unimod,BC_file names real :: frac=1.0 ! fraction to unimod variable diff --git a/ForestFire_ml.f90 b/ForestFire_ml.f90 old mode 100644 new mode 100755 diff --git a/Functions_ml.f90 b/Functions_ml.f90 old mode 100644 new mode 100755 diff --git a/GlobalBCs_ml.f90 b/GlobalBCs_ml.f90 old mode 100644 new mode 100755 index af94b09..b0e6a26 --- a/GlobalBCs_ml.f90 +++ b/GlobalBCs_ml.f90 @@ -268,7 +268,7 @@ subroutine GetGlobalData(year,month,ibc,used, & SIAtrend%nh4 =f0*SIAtrends(i0)%nh4 + f1*SIAtrends(i1)%nh4 if (MasterProc.and.first_call) then - write(unit=txtmsg,fmt="(a,i5,4f8.3)") & + write(unit=txtmsg,fmt="(a,2i5,3f8.3)") & "BC:trends SOx,NOx,NH3 ", iyr_trend, SIAtrend call PrintLog(txtmsg) end if diff --git a/GridValues_ml.f90 b/GridValues_ml.f90 old mode 100644 new mode 100755 diff --git a/Io_ml.f90 b/Io_ml.f90 old mode 100644 new mode 100755 diff --git a/LandPFT_ml.f90 b/LandPFT_ml.f90 old mode 100644 new mode 100755 diff --git a/Landuse_ml.f90 b/Landuse_ml.f90 old mode 100644 new mode 100755 diff --git a/MARS_Aero_water_ml.f90 b/MARS_Aero_water_ml.f90 old mode 100644 new mode 100755 diff --git a/MARS_ml.f90 b/MARS_ml.f90 old mode 100644 new mode 100755 index 24f70ef..ff26ddc --- a/MARS_ml.f90 +++ b/MARS_ml.f90 @@ -41,6 +41,11 @@ module MARS_ml implicit none private +!to test without the smoothing between different TNH4/TSO4 regimes +!before rv4.2beta (Jan 2013), the code was (or should be) equivalent to MARS_RATIO_SMOOTH=.false. + logical, parameter :: MARS_RATIO_SMOOTH=.true. + + real, parameter :: FLOOR = 1.0E-30 ! minimum concentration ! -30 from RPM @@ -515,7 +520,7 @@ subroutine rpmares ( SO4, HNO3, NO3, NH3, NH4, RH, TEMP, & High_Factor=0.0 endif - IF ( RATIO >RATIO_Low .and. RATIO < RATIO_High)then + IF ( RATIO >RATIO_Low .and. RATIO < RATIO_High .and. MARS_RATIO_SMOOTH)then DO_RATIO_High_2=.true. DO_RATIO_Low_2=.true. TSO4_HighA=TSO4*Ratio/RATIO_High diff --git a/Makefile b/Makefile old mode 100644 new mode 100755 index 3db27a9..af0c368 --- a/Makefile +++ b/Makefile @@ -6,13 +6,12 @@ PROG = Unimod include Makefile.SRCS ################################################### + LIBS = -lnetcdf -lnetcdff -INCL = -I/global/apps/netcdf/4.2.1.1/intel/13.0/include -LLIB = -L/global/apps/netcdf/4.2.1.1/intel/13.0/lib -#LIBS = -lnetcdf -#INCL = -I/global/apps/netcdf/3.6.2/include -#LLIB = -L/global/apps/netcdf/3.6.2/lib +INCL = -I$(NETCDF_ROOT)/include +LLIB = -L$(NETCDF_ROOT)/lib + F90 = mpif90 diff --git a/MassBudget_ml.f90 b/MassBudget_ml.f90 old mode 100644 new mode 100755 diff --git a/Met_ml.f90 b/Met_ml.f90 old mode 100644 new mode 100755 index 113695a..5e68a38 --- a/Met_ml.f90 +++ b/Met_ml.f90 @@ -206,9 +206,23 @@ subroutine MeteoRead(numt) real :: nsec ! step in seconds real :: buff(MAXLIMAX,MAXLJMAX)!temporary metfields - integer :: i, j + integer :: i, j, isw logical :: fexist + ! Soil water has many names. Some we can deal with: + ! (and all need to end up as SMI) + character(len=*), dimension(4), parameter :: possible_soilwater_uppr = (/& + "SMI1 " & + ,"SMI " & + ,"soil_water_content " & + ,"soil_wetness_surface" & + /) + character(len=*), dimension(3), parameter :: possible_soilwater_deep = (/& + "SMI3 " & + ,"SMI " & + ,"deep_soil_water_content" & + /) + nr=2 !set to one only when the first time meteo is read call_msg = "Meteoread" @@ -525,94 +539,147 @@ subroutine MeteoRead(numt) ! ! Start with shallow + if( USE_DUST .and. .not.USE_SOILWATER ) call StopAll("Inconsistent SM, DUST") if( USE_SOILWATER ) then SoilWaterSource = "IFS"! use as default? - namefield='SMI1' - call Getmeteofield(meteoname,namefield,nrec,ndim,& + do isw = 1, size(possible_soilwater_uppr) + namefield=possible_soilwater_uppr(isw) + if(MasterProc) write(*,*) "Met_ml: soil water search ", isw, trim(namefield) + call Getmeteofield(meteoname,namefield,nrec,ndim,& unit,validity, SoilWater_uppr(:,:,nr)) - - if(validity/=field_not_found)then - if( .not.foundSMI1 ) call PrintLog("Met: found SMI1", MasterProc) - foundSMI1=.true. - foundSoilWater_uppr = .true. - else - - namefield='soil_water_content' - call Getmeteofield(meteoname,namefield,nrec,ndim,& - unit,validity, SoilWater_uppr(:,:,nr)) - if(validity==field_not_found)then - namefield='soil_wetness_surface' - call Getmeteofield(meteoname,namefield,nrec,ndim,& - unit,validity, SoilWater_uppr(:,:,nr)) - endif - if(validity==field_not_found)then - if(MasterProc.and.numt==1)write(*,*)' WARNING: SoilWater_uppr not found ' - foundSoilWater_uppr = .false. - - else if ( trim(unit) /= "m" ) then - ! PARLAM/HIRLAM has metres of water in top 7.2 cm - - if(MasterProc.and.numt==1) & - write(*,*)'WARNING: Assuming SoilWater from IFS' - foundSoilWater_uppr = .true. - SoilWaterSource = "IFS" - else - foundSoilWater_uppr = .true. - SoilWaterSource = "PARLAM" - endif - endif + if(validity/=field_not_found) then ! found + if( index( namefield, "SMI" ) > 0 ) foundSMI1=.true. + foundSoilWater_uppr = .true. + if( .not.foundSMI1 ) & ! = 1st call + call PrintLog("Met: found SMI1:" // trim( namefield), MasterProc) + exit + end if + end do + !?if( .not.foundSMI1 ) call PrintLog("Met: found SMI1", MasterProc) + +!SW if(validity/=field_not_found)then +!SW if( .not.foundSMI1 ) call PrintLog("Met: found SMI1", MasterProc) +!SW foundSMI1=.true. +!SW foundSoilWater_uppr = .true. +!SW else +!SW namefield='SMI' ! e.g. RCA data have just SMI. Use for SMI1 and SMI3 +!SW call Getmeteofield(meteoname,namefield,nrec,ndim,& +!SW unit,validity, SoilWater_uppr(:,:,nr)) +!SW foundSMI1=.true. +!SW foundSoilWater_uppr = .true. +!SW print *, "TEST RCA SMI ", me, validity +!SW +!SW if ( validity==field_not_found)then +!SW +!SW namefield='soil_water_content' +!SW call Getmeteofield(meteoname,namefield,nrec,ndim,& +!SW unit,validity, SoilWater_uppr(:,:,nr)) +!SW if(validity==field_not_found)then +!SW namefield='soil_wetness_surface' +!SW call Getmeteofield(meteoname,namefield,nrec,ndim,& +!SW unit,validity, SoilWater_uppr(:,:,nr)) +!SW endif +!SW if(validity==field_not_found)then +!SW if(MasterProc.and.numt==1)write(*,*)' WARNING: SoilWater_uppr not found ' +!SW foundSoilWater_uppr = .false. +!SW +!SW else + + if( foundSoilWater_uppr .and. trim(unit) == "m" ) SoilWaterSource = "PARLAM" + +!SW ! PARLAM/HIRLAM has metres of water in top 7.2 cm +!SW +!SW if(MasterProc.and.numt==1) & +!SW write(*,*)'WARNING: Assuming SoilWater from IFS' +!SW foundSoilWater_uppr = .true. +!SW SoilWaterSource = "IFS" +!SW else +!SW foundSoilWater_uppr = .true. +!SW SoilWaterSource = "PARLAM" +!SW endif +!SW endif end if ! USE_SOILWATER first one + if ( USE_SOILWATER ) then !just deep here + do isw = 1, size(possible_soilwater_deep) + namefield=possible_soilwater_deep(isw) + if(DomainName == "HIRHAM" ) then + if(MasterProc.and.numt==1)write(*,*) " Rename soil water in HIRHAM" + namefield='soil_water_second_layer' + end if + if(MasterProc) write(*,*) "Met_ml: deep soil water search ", isw, trim(namefield) + + call Getmeteofield(meteoname,namefield,nrec,ndim,& + unit,validity, SoilWater_deep(:,:,nr)) + + if(validity/=field_not_found) then ! found + if( index( namefield, "SMI" ) > 0 ) foundSMI3=.true. + foundSoilWater_deep = .true. + if( .not.foundSMI3 ) & ! = 1st call + call PrintLog("Met: found SMI3:" // trim( namefield), MasterProc) + exit + end if + end do !======================================== !In the long term, all meteofile should have Soil Moisture Index defined - namefield='SMI3' - call Getmeteofield(meteoname,namefield,nrec,ndim,& - unit,validity, SoilWater_deep(:,:,nr)) - - if(validity/=field_not_found)then - if( .not.foundSMI3 ) call PrintLog("Met: found SMI3", MasterProc) - foundSMI3=.true. - foundSoilWater_deep = .true. - - else - - !Search for other fields which can be used for making SMI - namefield='deep_soil_water_content' - if(DomainName == "HIRHAM" ) then - if(MasterProc.and.numt==1)write(*,*) " Rename soil water in HIRHAM" - namefield='soil_water_second_layer' - end if - call Getmeteofield(meteoname,namefield,nrec,ndim,& - unit,validity, SoilWater_deep(:,:,nr)) - if(validity==field_not_found)then - if(MasterProc.and.numt==1)write(*,*)' WARNING: ',trim(namefield),' not found ' - foundSoilWater_deep = .false. - else - !<<<<<<< process SW <<<<<<<<<<<<<<<<<<<<<<< - foundSoilWater_deep = .true. - if ( trim(unit) == "m" ) then ! PARLAM has metres of water - SoilWaterSource = "PARLAM" - else if(unit(1:5)=='m3/m3')then - !IFS has a fairly complex soil water system, with field capacity of - ! up to 0.766 for organic soils. More medium soils have ca. 0.43 - ! Typical values in January are even down to 0.2. Best to use - ! SMI.... - SoilWaterSource = "IFS" - else ! units not defined yet - if(numt==1)write(*,*)trim(unit) - call StopAll("Need units for deep soil water") - endif - - if(MasterProc.and.numt==1) write(*,*)'max Met_ml Soilwater_deep: ' // & - trim(SoilWaterSource), maxval( SoilWater_deep(:,:,nr) ) - - endif !found deep_soil_water_content - - endif !SMI3 found +!SW namefield='SMI3' +!SW call Getmeteofield(meteoname,namefield,nrec,ndim,& +!SW unit,validity, SoilWater_deep(:,:,nr)) +!SW +!SW if(validity/=field_not_found)then +!SW if( .not.foundSMI3 ) call PrintLog("Met: found SMI3", MasterProc) +!SW foundSMI3=.true. +!SW foundSoilWater_deep = .true. +!SW +!SW else +!SW namefield='SMI' ! e.g. RCA data have just SMI. Use for SMI1 and SMI3 +!SW call Getmeteofield(meteoname,namefield,nrec,ndim,& +!SW unit,validity, SoilWater_uppr(:,:,nr)) +!SW foundSMI3=.true. +!SW foundSoilWater_deep = .true. +!SW print *, "TEST RCA SMI DEEP ", me, validity +!SW +!SW if ( validity==field_not_found)then +!SW +!SW !Search for other fields which can be used for making SMI +!SW namefield='deep_soil_water_content' +!SW if(DomainName == "HIRHAM" ) then +!SW if(MasterProc.and.numt==1)write(*,*) " Rename soil water in HIRHAM" +!SW namefield='soil_water_second_layer' +!SW end if +!SW call Getmeteofield(meteoname,namefield,nrec,ndim,& +!SW unit,validity, SoilWater_deep(:,:,nr)) +!SW if(validity==field_not_found)then +!SW if(MasterProc.and.numt==1)write(*,*)' WARNING: ',trim(namefield),' not found ' +!SW foundSoilWater_deep = .false. +!SW else +!SW !<<<<<<< process SW <<<<<<<<<<<<<<<<<<<<<<< +!SW foundSoilWater_deep = .true. + if( foundSoilWater_deep ) then + if ( trim(unit) == "m" ) then ! PARLAM has metres of water + SoilWaterSource = "PARLAM" + else if(unit(1:5)=='m3/m3')then +!SW !IFS has a fairly complex soil water system, with field capacity of +!SW ! up to 0.766 for organic soils. More medium soils have ca. 0.43 +!SW ! Typical values in January are even down to 0.2. Best to use +!SW ! SMI.... + SoilWaterSource = "IFS" + end if +!SW else ! units not defined yet +!SW if(numt==1)write(*,*)trim(unit) +!SW call StopAll("Need units for deep soil water") +!SW endif +!SW +!SW if(MasterProc.and.numt==1) write(*,*)'max Met_ml Soilwater_deep: ' // & +!SW trim(SoilWaterSource), maxval( SoilWater_deep(:,:,nr) ) +!SW + endif !found deep_soil_water_content +!SW +!SW endif !SMI3 found if ( DEBUG_SOILWATER.and.debug_proc ) then @@ -691,7 +758,7 @@ subroutine metvar(numt) real, dimension(MAXLIMAX,KMAX_MID) :: vrcv ! and in y direction real divt, p1, p2 - real prhelp_sum,divk(KMAX_MID),sumdiv + real prhelp_sum,divk(KMAX_MID),sumdiv,dB_sum real inv_METSTEP integer :: i, j, k, kk, nr,info, ii,jj,ii2,jj2 @@ -1204,6 +1271,8 @@ subroutine metvar(numt) !see http://www.ecmwf.int/research/ifsdocs/DYNAMICS/Chap2_Discretization3.html#959545 !(note that u_xmj and v_xmi have already been divided by xm here) + dB_sum=1.0/(B_bnd(KMAX_MID+1)-B_bnd(1))!normalisation factor for dB (should be one if entire atmosphere is included) + do j = 1,ljmax do i = 1,limax Pmid=Ps_extended(i,j)! without "-PT" @@ -1223,7 +1292,7 @@ subroutine metvar(numt) Etadot(i,j,KMAX_MID+1,nr)=0.0 do k=KMAX_MID,1,-1 - Etadot(i,j,k,nr)=Etadot(i,j,k+1,nr)-dB(k)*sumdiv+divk(k) + Etadot(i,j,k,nr)=Etadot(i,j,k+1,nr)-dB(k)*dB_sum*sumdiv+divk(k) Etadot(i,j,k+1,nr)=Etadot(i,j,k+1,nr)*(dA(k)/Pref+dB(k))/(dA(k)+dB(k)*Pmid) ! Etadot(i,j,k+1,nr)=Etadot(i,j,k+1,nr)/(Ps_extended(i,j)-PT) gives same result as sdot for sigma coordinates enddo @@ -2953,7 +3022,7 @@ subroutine GetCDF_short(varname,fileName,var,GIMAX,IRUNBEG,GJMAX,JRUNBEG & integer :: varID,ndims integer :: ncFileID,status real :: scale,offset - character *100 :: period_read + character *100 :: period_read=' ' character *200,save :: filename_save='notsaved' integer,save :: ncFileID_save=-99 diff --git a/ModelConstants_ml.f90 b/ModelConstants_ml.f90 old mode 100644 new mode 100755 index 66d46ae..ec3f748 --- a/ModelConstants_ml.f90 +++ b/ModelConstants_ml.f90 @@ -68,45 +68,60 @@ module ModelConstants_ml ! Namelist controlled: ! Some flags for model setup -! First, those which may be reset by emep_namelist.nml file +!------------ NAMELIST VARIABLES - can be reset by emep_namelist.nml file logical, public, save :: & FORECAST = .false. &! reset in namelist ,USE_SOILWATER = .false. &! ,USE_DEGREEDAY_FACTORS = .false. &! ,USE_FOREST_FIRES = .false. &! -! Soil NOx. Choose EURO for better spatial and temp res, but for -! global runs need global monthly. Variable USE_SOILNOX set from -! these below. - ,USE_EURO_SOILNOX = .true. & ! ok, but diff for global + Euro runs - ,USE_GLOBAL_SOILNOX = .false. & ! Need to design better switch - ,USE_SOILNOX = .true. & ! DO NOT ALTER: Set after config -! ,USE_SEASALT = .true. & ! ,USE_CONVECTION = .false. & ! false works best for Euro runs, +! +! Might sometimes change for scenario runs (e.g. EnsClim): + ,USE_AIRCRAFT_EMIS = .true. & ! Needs global file, see manual + ,USE_LIGHTNING_EMIS = .true. & ! +! ! More experimental: ,DO_SAHARA = .false. & ! Turn on/off BG Saharan Dust ,USE_ROADDUST = .false. & ! TNO Road Dust routine. So far with simplified "climate-correction" factor ,USE_DUST = .false. &! Experimental -!TFMM ,INERIS_SNAP1 = .false. & !(EXP_NAME=="TFMM"), & ! Switches off decadal trend ,INERIS_SNAP2 = .false. & !(EXP_NAME=="TFMM"), & ! Allows near-zero summer values ,USE_EMERGENCY = .false. & ! Emergency: Volcanic Eruption & Nuclear Accident. Under development. ,USE_AOD = .false. & ,USE_POLLEN = .false. & ! EXPERIMENTAL. Only works if start Jan 1 ,ANALYSIS = .false. & ! EXPERIMENTAL: 3DVar data assimilation +! ! Output flags ,SELECT_LEVELS_HOURLY = .false. ! for FORECAST, 3DPROFILES +! Soil NOx. Choose EURO for better spatial and temp res, but for +! global runs need global monthly. Variable USE_SOILNOX set from +! these below. +! +! Also, is scaling needed for EURO_SOILNOX? +! The Euro soil NO emissions are based upon average Nr-deposition calculated +! for the 2000s, as given in the AnnualNdep.nc files. For future years a +! new AnnualNdep.nc could be pre-calculated. A simpler but approximate +! way is to scale with some other factor, e.g. the ratio of emissions over +! some area (EMEP, or EU) in year YYYY divided by year 2005 values. +! Remember, soil-NO emissions are *very* uncertain. + + logical, public, save :: & + USE_EURO_SOILNOX = .true. & ! ok, but diff for global + Euro runs + ,USE_GLOBAL_SOILNOX = .false. & ! Need to design better switch + ,USE_SOILNOX = .true. ! DO NOT ALTER: Set after config + real, public, save :: EURO_SOILNOX_DEPSCALE = 1.0 ! + ! Methane background. real, public, save :: BGND_CH4 = -1 ! -1 gives defaults in BoundaryConditions_ml, - ! otherwise set to positive value in NML +! +!------------ END OF NAMELIST VARIABLES ------------------------------------! ! Some flags for model setup ! will be removed when code is sufficiently tested ! (for convection use foundconv in permanent code) logical, public, parameter :: & - USE_AIRCRAFT_EMIS = .false., & ! Needs global file, see manual - USE_LIGHTNING_EMIS = .true., & ! ok NO_CROPNH3DEP = .true., & ! Stop NH3 deposition for growing crops USE_SOILNH3 = .false., & ! DUMMY VALUES, DO NOT USE! USE_ZREF = .false., & ! testing @@ -305,7 +320,7 @@ module ModelConstants_ml ,DEBUG_DO3SE = .false. & ,DEBUG_DRYRUN = .false. & ! Skips fast chemistry to save some CPU ,DEBUG_ECOSYSTEMS = .false. & - ,DEBUG_FORESTFIRE = .false. & + ,DEBUG_FORESTFIRE = .true. & ,DEBUG_Kz = .false. & ,DEBUG_MY_DERIVED = .false. & ,DEBUG_DRYDEP = .false. & @@ -317,7 +332,7 @@ module ModelConstants_ml ,DEBUG_EMISTIMEFACS = .false. & ,DEBUG_EQUIB = .false. & !MARS, EQSAM etc. ,DEBUG_GETEMIS = .false. & - ,DEBUG_GRIDVALUES = .false. & + ,DEBUG_GRIDVALUES = .true. & ,DEBUG_IOPROG = .false. & ,DEBUG_LANDDEFS = .false. & ,DEBUG_LANDUSE = .false. & @@ -471,10 +486,10 @@ module ModelConstants_ml ! refer to output variables defined in Derived_ml. ! Hourly output types: Last 2 types (hourly inst.,hourly mean), ! refer to output variables defined in My_Outputs_ml. -! IOU_HOUR_PREVIOUS: Auxiliary field for hourly accumulated Derived output +! IOU_YEAR_LASTHH: Auxiliary field for hourly accumulated Derived output integer, public, parameter :: & IOU_INST=1, IOU_YEAR=2, IOU_MON=3, IOU_DAY=4, & ! Derived output - IOU_HOUR_PREVIOUS=5, & ! Aux. field + IOU_YEAR_LASTHH=5, & ! Aux. field IOU_HOUR=6, IOU_HOUR_MEAN=7 & ! Hourly output ,IOU_MAX_MAX=7 ! Max values for of IOU (for array declarations) @@ -483,8 +498,8 @@ module ModelConstants_ml !---------------------------------------------------------------------------- contains subroutine Config_ModelConstants(iolog) - character(len=120) :: txt, errmsg - logical :: file_exists + character(len=120) :: txt + !logical :: file_exists integer, intent(in) :: iolog ! for Log file NAMELIST /ModelConstants_config/ & @@ -492,8 +507,10 @@ subroutine Config_ModelConstants(iolog) ,MY_OUTPUTS & ! e.g. EMEPSTD, FORECAST, TFMM ,USE_SOILWATER, USE_DEGREEDAY_FACTORS, USE_FOREST_FIRES & ,USE_CONVECTION & + ,USE_AIRCRAFT_EMIS,USE_LIGHTNING_EMIS & ,DO_SAHARA, USE_ROADDUST, USE_DUST & - ,USE_EURO_SOILNOX, USE_GLOBAL_SOILNOX, USE_SEASALT ,USE_POLLEN & + ,USE_EURO_SOILNOX, USE_GLOBAL_SOILNOX, EURO_SOILNOX_DEPSCALE & + ,USE_SEASALT ,USE_POLLEN & ,INERIS_SNAP1, INERIS_SNAP2 & ! Used for TFMM time-factors ,IS_GLOBAL & ! Also for EMERGENCY ,MONTHLY_GRIDEMIS & ! was set in Emissions_ml diff --git a/MosaicOutputs_ml.f90 b/MosaicOutputs_ml.f90 old mode 100644 new mode 100755 diff --git a/My_Aerosols_ml.f90 b/My_Aerosols_ml.f90 old mode 100644 new mode 100755 diff --git a/My_Derived_ml.f90 b/My_Derived_ml.f90 old mode 100644 new mode 100755 index 4a407b0..6a5e0ba --- a/My_Derived_ml.f90 +++ b/My_Derived_ml.f90 @@ -415,43 +415,43 @@ module My_Derived_ml ! character(len=TXTLEN_DERIV), public, parameter, dimension(29) :: & VEGO3_WANTED = (/ & - "POD1_IAM_DF ",& - "POD1_IAM_MF ",& - "POD1_DF ",& - "POD1_CF ",& - "POD3_TC ",& - "SPOD15_birch ",& - "SPOD10_birch ",& - "SPOD15_spruce ",& + "POD1_IAM_DF ",& + "POD1_IAM_MF ",& + "POD1_DF ",& + "POD1_CF ",& + "POD3_TC ",& + "SPOD15_birch ",& + "SPOD10_birch ",& + "SPOD15_spruce ",& "SPOD10_spruce ",& - "SPOD15_crops ",& - "SPOD25_crops ",& + "SPOD15_crops ",& + "SPOD25_crops ",& !WIMMAX: "POD1_NEUR_SPRUCE",& - "POD1_NEUR_BIRCH",& - "POD1_ACE_PINE ",& - "POD1_ACE_OAK ",& - "POD1_ACE_BEECH ",& + "POD1_NEUR_BIRCH ",& + "POD1_ACE_PINE ",& + "POD1_ACE_OAK ",& + "POD1_ACE_BEECH ",& "POD1_CCE_SPRUCE ",& - "POD1_CCE_BEECH ",& - "POD1_MED_OAK ",& - "POD1_MED_PINE ",& - "POD1_MED_BEECH ",& + "POD1_CCE_BEECH ",& + "POD1_MED_OAK ",& + "POD1_MED_PINE ",& + "POD1_MED_BEECH ",& ! "POD3_TC30d ",& ! "POD3_TC55d ",& - "POD3_IAM_CR ",& + "POD3_IAM_CR ",& ! "POD3_IAM_CR30d ",& ! "POD3_IAM_CR55d ",& ! "POD6_IAM_CR ",& ! Not recommended - not robust ! "POD6_IAM_CR30d ",& ! Not recommended - not robust ! "POD6_IAM_CR55d ",& ! Not recommended - not robust - "MMAOT40_TC ",& - "MMAOT40_IAM_DF ",& - "MMAOT40_IAM_MF ",& - "MMAOT40_IAM_CR ",& - "EUAOT40_Crops ", & - "EUAOT40_Forests", & - "MMAOT40_IAM_WH " & + "MMAOT40_TC ",& + "MMAOT40_IAM_DF ",& + "MMAOT40_IAM_MF ",& + "MMAOT40_IAM_CR ",& + "EUAOT40_Crops ", & + "EUAOT40_Forests ", & + "MMAOT40_IAM_WH " & /) !NB -last not found. Could just be skipped, but kept !to show behaviour diff --git a/My_Outputs_ml.f90 b/My_Outputs_ml.f90 old mode 100644 new mode 100755 index e9cbfa1..25b2044 --- a/My_Outputs_ml.f90 +++ b/My_Outputs_ml.f90 @@ -70,12 +70,10 @@ module My_Outputs_ml integer, public, parameter :: & NSITES_MAX = 99 & ! Max. no surface sites allowed ,FREQ_SITE = 1 & ! Interval (hrs) between outputs -! ,NADV_SITE = NSPEC_ADV & ! No. advected species (1 up to NSPEC_ADV) -! ,NSHL_SITE = NSPEC_SHL & ! No. short-lived species - ,NADV_SITE = 0 & ! No. advected species (1 up to NSPEC_ADV) - ,NSHL_SITE = 0 & ! No. short-lived species + ,NADV_SITE = NSPEC_ADV & ! No. advected species (1 up to NSPEC_ADV) + ,NSHL_SITE = NSPEC_SHL & ! No. short-lived species ,NXTRA_SITE_MISC = 2 & ! No. Misc. met. params ( e.g. T2, d_2d) - ,NXTRA_SITE_D2D = 6 ! No. Misc. met. params ( e.g. T2, d_2d) + ,NXTRA_SITE_D2D = 9 ! No. Misc. met. params ( e.g. T2, d_2d) integer, public, parameter, dimension(NADV_SITE) :: & SITE_ADV = (/ (isite, isite=1,NADV_SITE) /) ! Everything @@ -101,8 +99,17 @@ module My_Outputs_ml "HMIX ",& "PSURF ", & "ws_10m ", & - "rh2mi", "O3 ", & - "COLUMN_NO2_k20 " /) + "rh2m ", & + "Emis_mgm2_BioNatC5H8 ", & + "Emis_mgm2_BioNatAPINENE", & + "Emis_mgm2_BioNatNO ",& + "Emis_mgm2_nox ",& +! "SoilWater_deep ","EVAP_CF ","EVAP_DF ", & +! "EVAP_BF ","EVAP_NF ","WDEP_PREC ", & +! "RH_GR ","CanopyO3_GR ","VPD_GR ","FstO3_GR ", & +! "RH_IAM_DF ","CanopyO3_IAM_DF","VPD_IAM_DF ","FstO3_IAM_DF ", & +! "COLUMN_CO_k20 ","COLUMN_C2H6_k20","COLUMN_HCHO_k20","COLUMN_CH4_k20 ", + "COLUMN_NO2_k20 " /) !/*** Aircraft outputs (used in Polinat_ml) !============================================================== @@ -199,7 +206,7 @@ subroutine set_output_defs implicit none character(len=144) :: errmsg ! Local error message - integer :: i,j,ash,nuc_conc,nuc_wdep,nuc_ddep,rn222,pm25,pm10,poll,psurf,igrp,idx ! Loop & group indexes + integer :: i,j,ash,nuc_conc,nuc_wdep,nuc_ddep,rn222,pm25,pm10,poll,igrp,idx ! Loop & group indexes character(len=256) :: name ! Volcano (vent) name real, parameter :: atwC=12.0 @@ -249,7 +256,7 @@ subroutine set_output_defs select case(MY_OUTPUTS) case("EMERGENCY") nlevels_hourly = 1+18 - nhourly_out=4+2 !PM*,AOD (&Z, PS) + nhourly_out=5 !PM*,AOD & Z ash=find_index("ASH",chemgroups(:)%name) nuc_conc=find_index("NUCRACT",chemgroups(:)%name) nuc_wdep=find_index("DDEP_NUCRACT",chemgroups(:)%name) @@ -312,28 +319,34 @@ subroutine set_output_defs case("TFMM") nhourly_out=28 nlevels_hourly = 1 ! nb zero is *not* one of levels + case("average_vs_instantaneous") + nhourly_out = 6 + nlevels_hourly = 1 case default - nhourly_out=1 + nhourly_out=2 nlevels_hourly = 1 ! nb zero is *not* one of levels endselect if(allocated(hr_out)) deallocate(hr_out) if(allocated(levels_hourly))deallocate(levels_hourly) - allocate(hr_out(nhourly_out),levels_hourly(nlevels_hourly)) + allocate(hr_out(0:nhourly_out),levels_hourly(nlevels_hourly)) hr_out(:)=Asc2D("none","none",-99,-99,-99,-99,-99,-99,"none",-99.9,-99.9) +! Write out surface pressure at each time step to define vertical coordinates. +! Varaible must be named PS as defined in formula-terms for sigma. + hr_out(0)=Asc2D("PS","D2D_inst",find_index("PSURF",f_2d(:)%name),& + ix1,ix2,iy1,iy2,1,"hPa",1.0,-9999.9) select case(MY_OUTPUTS) case("EMERGENCY") ! ix1=IRUNBEG;ix2=IRUNBEG+GIMAX-1 ! iy1=JRUNBEG;iy2=JRUNBEG+GJMAX-1 levels_hourly = (/(i-1,i=1,nlevels_hourly)/) - pm25 =find_index("PMFINE",chemgroups(:)%name) !NB There is no "PM25" group - pm10 =find_index("PM10" ,chemgroups(:)%name) - psurf=find_index("PSURF" ,f_2d(:)%name) + pm25=find_index("PMFINE",chemgroups(:)%name) !NB There is no "PM25" group + pm10=find_index("PM10" ,chemgroups(:)%name) !** name type ofmt !** ispec ix1 ix2 iy1 iy2 nk sellev? unit conv max - j=6;hr_out(:j) = (/& + j=5;hr_out(1:j) = (/& Asc2D("pm25_3km" ,"BCVugXXgroup",pm25,& ix1,ix2,iy1,iy2,NLEVELS_HOURLY,"ug/m3",1.0,-999.9),& Asc2D("pm10_3km" ,"BCVugXXgroup",pm10,& @@ -343,11 +356,7 @@ subroutine set_output_defs Asc2D("AOD_550nm" ,"AOD" ,00 ,& ix1,ix2,iy1,iy2,1," ",1.0 ,-9999.9),& Asc2D("z" ,"Z_MID" ,00 ,& - ix1,ix2,iy1,iy2,NLEVELS_HOURLY,"km",1e-3,-9999.9),& -!** PSURF needs always to be named PS as defined in formula-terms for sigma -!** PSURF needed for konversion from sigma to pressure or height - Asc2D("PS" ,"D2D" ,psurf ,& - ix1,ix2,iy1,iy2,1,"hPa",1.0,-9999.9)/) + ix1,ix2,iy1,iy2,NLEVELS_HOURLY,"km",1e-3,-9999.9)/) if(ash>0)then name="none" do i=1,size(chemgroups(ash)%ptr) @@ -386,7 +395,7 @@ subroutine set_output_defs ! Asc2D(name ,"BCVugXXgroup",igrp,& ! ix1,ix2,iy1,iy2,1,"uBqh/m2",1.0,-999.9)& ! /) - hr_out(j)=Asc2D(trim(name)//"_WDEP","D2D",& + hr_out(j)=Asc2D(trim(name)//"_WDEP","D2D_accum",& find_index("WDEP_"//trim(name),f_2d(:)%name),ix1,ix2,iy1,iy2,1,"",1.0,-999.9) enddo @@ -402,7 +411,7 @@ subroutine set_output_defs ! Asc2D(name ,"BCVugXXgroup",igrp,& ! ix1,ix2,iy1,iy2,1,"uBqh/m2",1.0,-999.9)& ! /) - hr_out(j)=Asc2D(trim(name)//"_DDEP","D2D",& + hr_out(j)=Asc2D(trim(name)//"_DDEP","D2D_accum",& find_index("DDEP_"//name,f_2d(:)%name),ix1,ix2,iy1,iy2,1,"",1.0,-999.9) enddo endif @@ -417,7 +426,7 @@ subroutine set_output_defs poll =find_index("POLLEN_B",species_adv(:)%name) !** name type ofmt !** ispec ix1 ix2 iy1 iy2 nk sellev? unit conv max - j=10;hr_out(:j) = (/& + j=10;hr_out(1:j) = (/& Asc2D("o3_3km" ,"BCVugXX",IXADV_O3 ,& ix1,ix2,iy1,iy2,NLEVELS_HOURLY,"ug",to_ug_ADV(IXADV_O3) ,600.0*2.0),& Asc2D("no_3km" ,"BCVugXX",IXADV_NO ,& @@ -468,22 +477,18 @@ subroutine set_output_defs pm10 =find_index("SURF_ug_PM10_rh50" ,f_2d(:)%name) !** name type ofmt ispec !** ix1 ix2 iy1 iy2 nk sellev? unit conv max - hr_out = (/& + hr_out(1:) = (/& Asc2D("O3" ,"BCVugXX",IXADV_O3 ,& ix1,ix2,iy1,iy2,NLEVELS_HOURLY,"ug",to_ug_ADV(IXADV_O3) ,600.0*2.0),& Asc2D("NO2" ,"BCVugXX",IXADV_NO2 ,& ix1,ix2,iy1,iy2,NLEVELS_HOURLY,"ug",to_ug_ADV(IXADV_NO2),600.0*1.91),& - Asc2D("PM25","D2D",pm25,& + Asc2D("PM25","D2D_inst",pm25,& ix1,ix2,iy1,iy2,1 ,"ug",1.0 ,-999.9),& - Asc2D("PM10","D2D" ,pm10,& + Asc2D("PM10","D2D_inst",pm10,& ix1,ix2,iy1,iy2,1 ,"ug",1.0 ,-999.9)/) -! if(MasterProc)then -! write(*,*)(i,hr_out(i),i=1,nhourly_out) -! call CheckStop("AMVB") -! endif case("IMPACT2C") ! Dave's starting set. Uses Out3D to get 3m and 45m concs. levels_hourly = (/ (i, i= 0,nlevels_hourly-1) /) ! -1 will give surfac - hr_out= (/ & + hr_out(1:)= (/ & Asc2D("o3_3dppb" ,"Out3D",O3 ,& ix1,ix2,iy1,iy2,nlevels_hourly,"ppbv", PPBINV,600.0*2.0) & ,Asc2D("no_3dppb" ,"Out3D",& @@ -492,14 +497,14 @@ subroutine set_output_defs NO2 ,ix1,ix2,iy1,iy2,nlevels_hourly,"ppbv",PPBINV ,600.0*1.91) & ,Asc2D("ho2_3dppt" ,"Out3D",& !NOTE ppt HO2 ,ix1,ix2,iy1,iy2,nlevels_hourly,"pptv",PPBINV ,600.0*1.91) & - ,Asc2D("HMIX","D2D",find_index("HMIX",f_2d(:)%name), & + ,Asc2D("HMIX","D2D_inst",find_index("HMIX",f_2d(:)%name), & ix1,ix2,iy1,iy2,1,"m",1.0,10000.0) & /) case("3DPROFILES") ! nb Out3D uses totals, e.g. O3, not IXADV_O3 ! Number of definitions must match nhourly_out set above levels_hourly = (/ (i, i= 0,nlevels_hourly-1) /) ! -1 will give surfac - hr_out= (/ & + hr_out(1:)= (/ & Asc2D("o3_3dppb" ,"Out3D",O3 ,& ix1,ix2,iy1,iy2,nlevels_hourly,"ppbv", PPBINV,600.0*2.0) & ,Asc2D("no2_3dppb" ,"Out3D",& @@ -519,7 +524,7 @@ subroutine set_output_defs !!!! (As a test I tried both pmfine two ways, one as D2D and the other !!!! as ADVugXXXgroup. The results were identical.) - hr_out = (/ & + hr_out(1:) = (/ & Asc2D("o3_3m", "ADVppbv", IXADV_o3, & ix1,ix2,iy1,iy2,1, "ppbv",PPBINV,600.0) & ,Asc2D("no3_f" ,"ADVugXX",IXADV_NO3_F ,& @@ -531,38 +536,38 @@ subroutine set_output_defs ,Asc2D("so4_f" ,"ADVugXX",IXADV_SO4 ,& ix1,ix2,iy1,iy2,1,"ug/m3",to_ug_ADV(IXADV_SO4) ,-999.9) & ! Organics - ,Asc2D("OM25_3m" ,"D2D", find_index("SURF_ug_PART_OM_F",f_2d(:)%name), & + ,Asc2D("OM25_3m" ,"D2D_inst",find_index("SURF_ug_PART_OM_F",f_2d(:)%name), & ix1,ix2,iy1,iy2,1, "ug/m3",1.0,-999.9) & - ,Asc2D("EC25_3m" ,"D2D", find_index("SURF_ug_ECFINE",f_2d(:)%name), & + ,Asc2D("EC25_3m" ,"D2D_inst",find_index("SURF_ug_ECFINE",f_2d(:)%name), & ix1,ix2,iy1,iy2,1, "ug/m3",1.0,-999.9) & - ,Asc2D("OM_c_3m" ,"D2D", find_index("SURF_ug_OMCOARSE",f_2d(:)%name), & + ,Asc2D("OM_c_3m" ,"D2D_inst",find_index("SURF_ug_OMCOARSE",f_2d(:)%name), & ix1,ix2,iy1,iy2,1, "ug/m3",1.0,-999.9) & - ,Asc2D("EC_c_3m" ,"D2D", find_index("SURF_ug_ECCOARSE",f_2d(:)%name), & + ,Asc2D("EC_c_3m" ,"D2D_inst",find_index("SURF_ug_ECCOARSE",f_2d(:)%name), & ix1,ix2,iy1,iy2,1, "ug/m3",1.0,-999.9) & - ,Asc2D("NatDust_f","D2D",find_index("SURF_ug_DUST_NAT_F",f_2d(:)%name), & + ,Asc2D("NatDust_f","D2D_inst",find_index("SURF_ug_DUST_NAT_F",f_2d(:)%name), & ix1,ix2,iy1,iy2,1, "ug/m3",1.0,-999.9) & - ,Asc2D("NatDust_c","D2D",find_index("SURF_ug_DUST_NAT_C",f_2d(:)%name), & + ,Asc2D("NatDust_c","D2D_inst",find_index("SURF_ug_DUST_NAT_C",f_2d(:)%name), & ix1,ix2,iy1,iy2,1, "ug/m3",1.0,-999.9) & - ,Asc2D("RoadDust_f","D2D",find_index("SURF_ug_DUST_ROAD_F",f_2d(:)%name), & + ,Asc2D("RoadDust_f","D2D_inst",find_index("SURF_ug_DUST_ROAD_F",f_2d(:)%name), & ix1,ix2,iy1,iy2,1, "ug/m3",1.0,-999.9) & - ,Asc2D("RoadDust_c","D2D",find_index("SURF_ug_DUST_ROAD_C",f_2d(:)%name), & + ,Asc2D("RoadDust_c","D2D_inst",find_index("SURF_ug_DUST_ROAD_C",f_2d(:)%name), & ix1,ix2,iy1,iy2,1, "ug/m3",1.0,-999.9) & - ,Asc2D("SeaSalt_f","D2D",find_index("SURF_ug_SEASALT_F",f_2d(:)%name), & + ,Asc2D("SeaSalt_f","D2D_inst",find_index("SURF_ug_SEASALT_F",f_2d(:)%name), & ix1,ix2,iy1,iy2,1, "ug/m3",1.0,-999.9) & - ,Asc2D("SeaSalt_c","D2D",find_index("SURF_ug_SEASALT_C",f_2d(:)%name), & + ,Asc2D("SeaSalt_c","D2D_inst",find_index("SURF_ug_SEASALT_C",f_2d(:)%name), & ix1,ix2,iy1,iy2,1, "ug/m3",1.0,-999.9) & ! Sums - ,Asc2D("PM25_3m" ,"D2D", find_index("SURF_ug_PM25_rh50",f_2d(:)%name), & + ,Asc2D("PM25_3m" ,"D2D_inst",find_index("SURF_ug_PM25_rh50",f_2d(:)%name), & ix1,ix2,iy1,iy2,1, "ug/m3",1.0,-999.9) & - ,Asc2D("PMFINE" ,"D2D", find_index("SURF_ug_PMFINE",f_2d(:)%name), & + ,Asc2D("PMFINE" ,"D2D_inst",find_index("SURF_ug_PMFINE",f_2d(:)%name), & ix1,ix2,iy1,iy2,1, "ug/m3",1.0,-999.9) & - ,Asc2D("pm_h2o_3m","D2D",find_index("SURF_PM25water",f_2d(:)%name), & + ,Asc2D("pm_h2o_3m","D2D_inst",find_index("SURF_PM25water",f_2d(:)%name), & ix1,ix2,iy1,iy2,1,"ug",1.0,-999.9) & - ,Asc2D("PM25dry_3m","D2D",find_index("SURF_ug_PM25",f_2d(:)%name), & + ,Asc2D("PM25dry_3m","D2D_inst",find_index("SURF_ug_PM25",f_2d(:)%name), & ix1,ix2,iy1,iy2,1, "ug/m3",1.0,-999.9) & - ,Asc2D("PM10_3m" ,"D2D", find_index("SURF_ug_PM10_rh50",f_2d(:)%name), & + ,Asc2D("PM10_3m" ,"D2D_inst",find_index("SURF_ug_PM10_rh50",f_2d(:)%name), & ix1,ix2,iy1,iy2,1, "ug/m3",1.0,-999.9) & - ,Asc2D("PM10dry_3m","D2D",find_index("SURF_ug_PM10",f_2d(:)%name), & + ,Asc2D("PM10dry_3m","D2D_inst",find_index("SURF_ug_PM10",f_2d(:)%name), & ix1,ix2,iy1,iy2,1, "ug/m3",1.0,-999.9) & ,Asc2D("no_3m" ,"ADVppbv", IXADV_NO,& ix1,ix2,iy1,iy2,1,"ppbv",PPBINV ,600.0*1.91) & @@ -572,17 +577,38 @@ subroutine set_output_defs ix1,ix2,iy1,iy2,1, "degC",1.0 ,100.0) & ,Asc2D("ws_10m", "ws_10m", 00, & ix1,ix2,iy1,iy2,1, "m/s",1.0 ,100.0) & - ,Asc2D("HMIX","D2D", find_index("HMIX",f_2d(:)%name), & + ,Asc2D("HMIX","D2D_inst",find_index("HMIX",f_2d(:)%name), & ix1,ix2,iy1,iy2,1, "m",1.0,10000.0) & - ,Asc2D("USTAR_NWP","D2D", find_index("USTAR_NWP",f_2d(:)%name), & + ,Asc2D("USTAR_NWP","D2D_inst",find_index("USTAR_NWP",f_2d(:)%name), & ix1,ix2,iy1,iy2,1, "m/s",1.0,-999.9) & - ,Asc2D("Kz_m2s","D2D",find_index("Kz_m2s",f_2d(:)%name), & + ,Asc2D("Kz_m2s","D2D_inst",find_index("Kz_m2s",f_2d(:)%name), & ix1,ix2,iy1,iy2,1, "m2/s",1.0,-999.9) & /) + case("average_vs_instantaneous") + ! Example of different hourly output types - hourly average and hourly instantaneous value + ! variables of type "ADVppbv" will be outputted instantaneous and + ! variables of type "D2D_mean" will be hourly averages. + hr_out(1:) = (/ & + Asc2D("SURF_ppb_O3_inst" ,"ADVppbv", IXADV_O3, & + ix1,ix2,iy1,iy2,1,"ppbv",PPBINV,600.0) ,& + Asc2D("SURF_ppb_O3_mean" ,"D2D_mean",find_index("SURF_ppb_O3",f_2d(:)%name), & + ix1,ix2,iy1,iy2,1,"ppbv",1.0 ,600.0), & + Asc2D("SURF_ug_NO3_F_inst","ADVugXX",IXADV_NO3_F,& + ix1,ix2,iy1,iy2,1,"ug/m3",to_ug_ADV(IXADV_NO3_F),-999.9), & + Asc2D("SURF_ug_NO3_F_mean","D2D_mean",find_index("SURF_ug_NO3_F",f_2d(:)%name),& + ix1,ix2,iy1,iy2,1,"ug/m3",1.0,-999.9), & + Asc2D("SURF_ug_PMFINE_inst" ,"ADVugXX",IXADV_NO3_F ,& + ix1,ix2,iy1,iy2,1,"ug/m3",to_ug_ADV(IXADV_NO3_F) ,-999.9), & + Asc2D("SURF_ug_PMFINE_mean","D2D_mean",find_index("SURF_ug_PMFINE",f_2d(:)%name),& + ix1,ix2,iy1,iy2,1,"ug/m3",1.0,-999.9) & + /) case default - hr_out = (/ & + hr_out(1:) = (/ & Asc2D("o3_3m", "ADVppbv", IXADV_o3, & - ix1,ix2,iy1,iy2,1, "ppbv",PPBINV,600.0) & + ix1,ix2,iy1,iy2,1, "ppbv",PPBINV,600.0), & + Asc2D("no2_col" ,"COLUMN",IXADV_NO2 ,& + ix1,ix2,iy1,iy2,1,"molec/cm2",to_molec_cm2,-999.9) & +!! ix1,ix2,iy1,iy2,1,"ug",to_ug_ADV(IXADV_NO2),-999.9) & /) endselect @@ -590,9 +616,9 @@ subroutine set_output_defs ! Was the array set? R: %name/=none ! Was the D2D/Group found? R: %spec>=0 ! Is the variable name unique? R: find_index((i)%name,(:)%name)==i - do i=1,nhourly_out + do i=0,nhourly_out if(MasterProc) write(*,*) "TESTHH O3 ATEND", i, nlevels_hourly - idx = find_index(hr_out(i)%name,hr_out(:)%name) + idx = find_index(hr_out(i)%name,hr_out(:)%name)-1 if((hr_out(i)%name/="none").and.(hr_out(i)%spec>=0).and.& (idx==i)) cycle write(errmsg,"(A,4(1X,A,'=',I0),2(1X,A,':',A))")& diff --git a/My_SOA_ml.f90 b/My_SOA_ml.f90 old mode 100644 new mode 100755 index 6fb0ea3..b5adfea --- a/My_SOA_ml.f90 +++ b/My_SOA_ml.f90 @@ -1,484 +1,78 @@ +! +!*****************************************************************************! +!* +!* Copyright (C) 2007-2011 met.no +!* +!* Contact information: +!* Norwegian Meteorological Institute +!* Box 43 Blindern +!* 0313 OSLO +!* NORWAY +!* email: emep.mscw@met.no +!* http://www.emep.int +!* +!* This program is free software: you can redistribute it and/or modify +!* it under the terms of the GNU General Public License as published by +!* the Free Software Foundation, either version 3 of the License, or +!* (at your option) any later version. +!* +!* This program is distributed in the hope that it will be useful, +!* but WITHOUT ANY WARRANTY; without even the implied warranty of +!* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +!* GNU General Public License for more details. +!* +!* You should have received a copy of the GNU General Public License +!* along with this program. If not, see . +!*****************************************************************************! module OrganicAerosol_ml - ! Calculates the amount of condensible species in the gas and aerosol phases. - ! - ! References: - ! S2007: Simpson, D. et al., JGR, 2007, - ! B2012: Bergström, R. et al., ACP (in press Sept.), 2012, - ! S2012: Simpson, D. et al., ACP, 2012, - ! - ! - ! Usage: call OrganicAerosol from Runchem, after setup of column data - ! (for meteorology, etc.). The subroutine initialises itself on the - ! first call and thereafter modifies two external variables: - ! xn(k,SOA) : the concentrations of SOA (For a 1-d column array(in EMEP)) - ! Fgas(k,X) : The fraction of X which is gas and not aeorosol - ! - ! ---------------------------------------------------------------------------- - ! - ! From Gas/Particle theory, A/G = K.COA, - ! therefore, Fgas = G/(G+A) = 1/(1+K.COA) - ! - !----------------------------------------------------------------------------- - ! NB- we exclude use of gamma for now, but leave commented out code - !----------------------------------------------------------------------------- - ! - ! Dave Simpson, August 2001 -- 2012 - ! Robert Bergström 2010 -- 2012 - ! - !-------------------------------------------------------------------------- - - ! Functions + GridValues + PT only for BGNDOC - use Functions_ml, only: StandardAtmos_kPa_2_km !ds for use in Hz scaling - use ChemFields_ml, only : Fgas3d ! stores 3-d between time-steps - use ChemChemicals_ml, only : species ! for molwts - use ChemSpecs_tot_ml, S1 => FIRST_SEMIVOL , S2 => LAST_SEMIVOL - - use ChemGroups_ml, only : & - NONVOLPCM_GROUP & - ,NVABSOM_GROUP & ! nonvolatile absorbing species for OA partitioning - ,ASOA => ASOA_GROUP & - ,BSOA => BSOA_GROUP & - ,ECFINE => ECFINE_GROUP & - ,SVFFUELOA25 => SVFFUELOA25_GROUP & ! semi-volatile FFUELOA25 - ,SVWOODOA25 => SVWOODOA25_GROUP & ! for VBS semivolatile WOOD BURNING OA (primary + aged!) - ,SVFFIREOA25 => SVFFIREOA25_GROUP & ! for VBS semivolatile Wildfire OA (primary + aged!) - ,FFUELEC => FFUELEC_GROUP & - ,NVFFUELOC25 => NVFFUELOC25_GROUP & ! non-vol. FFUELOC emis. (in PM2.5) - ,NVFFUELOCCO => NVFFUELOC_COARSE_GROUP & ! non-vol. " " , (PM2.5-10 frac.) - ,NVWOODOC25 => NVWOODOC25_GROUP & ! non-vol. WOODOC emissions - ! (zero in VBS-PAX type runs) - ,NVFFIREOC25 => NVFFIREOC25_GROUP ! only non-vol. FFIREOC emissions - ! (zero in VBS-PAX type runs) - - use GridValues_ml, only: A_mid,B_mid - use ModelConstants_ml, only : PT - + !-------------------------------------------------------------------------- + ! This module is fake - for initial 2011 public-domain ozone model only, pending + ! decision as to which SOA scheme to release as default. + ! Contact David.Simpson@met.no for more information if interested in SOA + ! schemes + !-------------------------------------------------------------------------- use ModelConstants_ml, only : CHEMTMIN, CHEMTMAX, & - MasterProc, DEBUG => DEBUG_SOA, & K2 => KMAX_MID, K1 => KCHEMTOP - use Par_ml, only : LIDIM => MAXLIMAX, LJDIM => MAXLJMAX - use PhysicalConstants_ml, only : AVOG, RGAS_J - use Setup_1dfields_ml, only : itemp, xn => xn_2d, Fgas, Fpart - use TimeDate_ml, only: current_date + use PhysicalConstants_ml, only : AVOG + use Setup_1dfields_ml, only : itemp, xn => xn_2d + use ChemChemicals_ml, only : species ! for molwts + use ChemSpecs_tot_ml, A1 => FIRST_SEMIVOL , A2 => LAST_SEMIVOL implicit none - private - !/-- subroutines - - public :: Init_OrganicAerosol - public :: OrganicAerosol + public :: Init_OrganicAerosol + public :: OrganicAerosol !/-- public - logical, public, parameter :: ORGANIC_AEROSOLS = .true. - - ! We store some values in 3-D fields, to allow the next G/P partitioning - ! calculation to start off with values of COA, mw and Fgas which - ! are about right. Ensures that very few iterations are needed. - -! real,public, save, dimension(S1:S2,LIDIM,LJDIM,K1:K2) :: & -! Grid_SOA_Fgas !EXC Grid_SOA_gamma - - real,public, save, allocatable, dimension(:,:,:) :: Grid_COA - - real, private, allocatable, dimension(:), save :: & - COA & ! Org. aerosol, ug/m3 - ! (this version does not include EC as absorber) - ,BGND_OC & ! FAKE FOR NOW, 0.50 ugC/m3 at surface - ,BGND_OA ! Assumed OA/OC=2, -> 1 ug/m3 - - - !VBS real, private, dimension(S1:S2,K1:K2), save :: & - ! TMP - we assign Fpart for all species for now, since - ! it makes it easier to code for nonvol and vol + logical, public, parameter :: ORGANIC_AEROSOLS = .false. ! From Setup_1dfields now -! real, private, dimension(1:NSPEC_TOT,K1:K2), save :: & -! Fgas & ! Fraction in gas-phase -! ,Fpart! & ! Fraction in gas-phase - !VBS ,tabRTpL ! = 1.0e-6 * R.T(i)/pL(Ti) for all temps i: - !EXC ,gamma & ! activity coefficient - - real, parameter, public :: SMALLFN = 1.0e-20 ! Minimum value of ug allowed - - !/-- private - - ! ug = array for aerosol masses (ug/m3). Includes non-volatile compounds: - ! TMP??? Excluding NVOL for now? - - real, private,allocatable, dimension(:,:), save :: ug_semivol -! - use new NONVOLOC grpup to define: - integer, private, parameter :: NUM_NONVOLPCM = size(NONVOLPCM_GROUP) - integer, private, parameter :: NUM_NVABSOM = size(NVABSOM_GROUP) -! integer, private, parameter, dimension(NUM_NONVOL) :: & -! NONVOL = (/ NONVOLOC_GROUP, NONVOLEC_GROUP /) ! OC+EC in partitioning OM - real, private,allocatable, dimension(:,:), save :: ug_nonvol - !real, private, dimension(K1:K2), save :: ug_ecf ! CityZen added +! real, public, dimension(A1:A2,K1:K2), save :: Fgas ! Fraction in gas-phase - real, private, save, dimension(S1:S2,CHEMTMIN:CHEMTMAX) :: tabCiStar - - integer, private, save :: NITER = 2 ! No. iterations for Ksoa - real, private, save :: xn2molem3 ! Conversion from molec/cm3 to mole/m3 - real, private, save :: xn2ugC, ugC2xn - - - !/-- DEBUG variables - ! Usually, DEBUG_SOA gives extra outputs. debug_flag is used to allow - ! some extra outputs for a gven i,j - set in CTM model. - - character(len=20), public, save :: soa_errmsg = "ok" - character(len=*), public, parameter :: SOA_MODULE_FLAG="VBS" + character(len=*), public, parameter :: SOA_MODULE_FLAG="NotUsed" contains !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< - subroutine Init_OrganicAerosol(i,j,debug_flag) - integer, intent(in) :: i,j - logical, intent(in) :: debug_flag - integer :: is, it, k - real, parameter :: kJ = 1000.0 - real,allocatable, dimension(:), save :: p_kPa, h_km ! for standard atmosphere - logical, save :: my_first_call = .true. - - - if ( my_first_call ) then - allocate(COA(K1:K2),BGND_OC(K1:K2),BGND_OA(K1:K2)) - allocate(ug_semivol(S1:S2,K1:K2)) - allocate(Grid_COA(LIDIM,LJDIM,K1:K2)) - allocate(ug_nonvol(NUM_NVABSOM,K1:K2)) - allocate(p_kPa(K2), h_km(K2)) - - !========================================================================= - ! Set up background OM - ! Need to convert aeros to ug/m3 or ugC/m3. Since xn is in molecules/cm3 - ! we divide by AVOG to get mole/cm3, multiply by 1e6 to get mole/m3, - ! by mol. weight to get g/m3 and by 1e6 to get ug/m3 - - xn2molem3 = 1.0/AVOG * 1.0e6 - xn2ugC = xn2molem3 * 12.0 * 1.0e6 - ugC2xn = 1/xn2ugC - - ! Use Standard Atmosphere to get average heights of layers - - p_kPa(:) = 0.001*( A_mid(:)+B_mid(:)*101325.0 ) ! Pressure in kPa - h_km = StandardAtmos_kPa_2_km(p_kPa) - BGND_OC(:)= 0.5 * 1.005 ! ng/m3 !!! will give 0.5 ugC/m3 at z=0 m - - do k = K1, K2 - BGND_OC(k) = BGND_OC(k) * exp( -h_km(k)/9.1 ) - if(DEBUG .and. MasterProc ) write(*,"(a,i4,2f8.3)") & - "BGND_OC ", k, h_km(k), BGND_OC(k) - end do - BGND_OA(:) = 2*BGND_OC(:) ! Assume OA/OC = 2 for bgnd - - do k = K1, K2 - Grid_COA(:,:,k) = BGND_OA(k) ! Use OA, not OC here - end do - - !========================================================================= - ! Set up Tables for Fcond - ! Ci = 1.0e6*P0/RT - ! Now, pi(T) = Ai exp(-Hi/RT) - ! And pi(T) = Pi(Tref) * exp( H/RT * (1/Tref - 1/T) ) - ! -> Ci(T) = Ci(Tref) * Tref/T * exp(...) - - do is=S1,S2 - do it=CHEMTMIN,CHEMTMAX - - ! C*-values are given for 298K according to most(?) publications. - tabCiStar(is,it) = species(is)%CiStar * 298./it * & - exp( species(is)%DeltaH * kJ/RGAS_J * (1.0/298. - 1.0/it) ) - end do - end do - - - if ( MasterProc ) then - do is = S1, S2 - write(6,"(a,i4,a20,f7.1,i3,8es10.2)") & - " Tab SOA: MW, Carbons, C*:", is, trim(species(is)%name), & - species(is)%molwt, species(is)%carbons, & - tabCiStar(is,273), tabCiStar(is,303) - end do - end if - - - ! print *, "FEB2012 ALLOCATE ", S1,S2, K1, K2 - allocate( Fgas3d(S1:S2,LIDIM,LJDIM,K1:K2) ) - - !+ initial guess (1st time-step only) - ! Fgas3D is only defined for the semivol stuff, so no need for nonvol here - ! We need to assume something on 1st time-step though: - Fgas3d = 1.0 - - ! Initial values. Should not change except for semi-volatiles - - Fpart(:,:) = 0.0 - Fpart(NONVOLPCM_GROUP,:) = 1.0 - Fgas(:,:) = max(0.0, 1.0 - Fpart(:,:) ) - - !VBS Grid_avg_mw = 250.0 ! Da - !EXC Grid_SOA_gamma = 1.0 - !========================================================================= - my_first_call = .false. - end if ! my_first_call - - ! We need to set Fgas at start of each Runchem i,j loop, as it is - ! used for rcemis: - - Fgas(S1:S2,:) = Fgas3d(S1:S2,i,j,:) ! Semivolatiles only in 3D Fgas - Fpart(S1:S2,:) = 1-Fgas(S1:S2,:) - -! Fgas(NONVOLPCM_GROUP,:) = 0.0 ! not needed, shouldn't change - - end subroutine Init_OrganicAerosol - - !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< - subroutine OrganicAerosol(i_pos,j_pos,debug_flag) - - integer, intent(in) :: i_pos, j_pos - logical, intent(in) :: debug_flag + !+ Driver routine for Secondary Organic Aerosol module - integer :: i, k, iter, ispec ! loop variables -! real, save :: molcc2ngm3 = 1.0e9*1.0e6/AVOG !molecules/cc-> ng/m3 - real, save :: molcc2ugm3 = 1.0e12/AVOG !molecules/cc-> ng/m3 -! real, save :: ngm32molcc = AVOG/1.0e15 !molecules/cc <- ng/m3 -! real, save :: ugm32molcc = AVOG/1.0e12 !molecules/cc <- ng/m3 - real :: Ksoa - integer :: nmonth, nday, nhour, seconds - - ! Outputs: -! real :: surfOC25, surfASOA, surfBSOA, surfFFUELOC25, surfWOODOC25, surfBGNDOC, surfOFFUELOA25_OC, surfFFUELOA25_OC, surfOWOODOA25, surfWOODOA25, surfOFFIREOA25, surfFFIREOA25 - real :: surfOC25, surfASOA, surfBSOA, surfBGNDOC, surfFFUELOA25_OC, surfWOODOA25_OC, surfFFIREOA25_OC - - - nmonth = current_date%month - nday = current_date%day - nhour = current_date%hour - seconds = current_date%seconds - - if( DEBUG .and. debug_flag) write(unit=*,fmt=*) "Into SOA" - -! Note that xn(SOA) is not strictly needed with the method we have, but it -! is a good first guess of the sum of the condensed phases, enables easy output -! and saves the need for iteration. -! -! Remember also that Fgas is saved, so will initially have been set by the -! preceding call to OrganicAerosol for a different set of i,j. - - ! 1st guesses: - - COA(:) = Grid_COA(i_pos,j_pos,:) - - - ! ============ Non-volatile species first ============================ - ! NVABSOM - Only include fine OM! That is no EC and no coarse OM! - - do i = 1, NUM_NVABSOM ! OA/OC for POC about 1.333 - ispec = NVABSOM_GROUP(i) - - ug_nonvol(i,:) = molcc2ugm3 * xn(ispec,:)*species(ispec)%molwt - - end do - - ! ============ SOA species now, iteration needed =================== - - do iter = 1, NITER - - - ! Fgas = G/(G+A) = 1/(1+K.COA) - ! K = tabRTpL/(mw*gamma) - - do ispec = S1, S2 - - !VBS Fgas(ispec,:) = 1.0/(1.0+tabRTpL(ispec,:)*COA(:)/avg_mw(:) ) - !EXC *gamma(:,ispec)) ) - - Fpart(ispec,:) = COA(:)/( COA(:)+tabCiStar(ispec,itemp(:)) ) - - ug_semivol(ispec,:) = molcc2ugm3 * xn(ispec,:)*species(ispec)%molwt & - * Fpart(ispec,:) - - end do ! ispec - - !ng(:,:) = max(SMALLFN, ng(:,:)) - - ! New estimate of COA (in ug/m3) and avg_mw (g/mole): - ! (nb. xn in molecules/cm3) - - do k = K1,K2 - - COA(k) = sum( ug_semivol(:,k) ) + sum( ug_nonvol(:,k) ) + BGND_OA(k) - - !VBS Nmoles = ( sum( Fpart(VOL,k) * xn(VOL,k) ) + sum( xn(NONVOLOC,k) ) & - !VBS + BGND_OC_ng(k)*ngm32molcc/250.0 ) & ! Assumed MW - !VBS * xn2molem3 - !VBS avg_mw(k) = 1.0e-6 * COA(k)/Nmoles - - end do !k - ! ==================================================================== - - if( DEBUG .and. debug_flag ) then - - if( iter == NITER .and. seconds == 0 ) then - write(unit=6,fmt="(a,i2,a,3i3,f7.2)") "Iteration ", Niter, & - " ======== ",nmonth, nday, nhour, itemp(K2) - write(unit=6,fmt="(a3,a15,3a10,a4,4a10)") "SOA","Species", "xn", & - "Ci* ", "Ki"," ", "Fpart", "ng" - - do i = 1, NUM_NONVOLPCM - ispec = NONVOLPCM_GROUP(i) -! write(unit=6,fmt="(a4,i3,a15,es10.2,2f10.3,a4,es10.3,f13.4)")& -! "NVOL", ispec,& -! species(ispec)%name, xn(ispec,K2),-999.999, & -! -999.999, " => ", Fpart(ispec,K2), 1000.0*ug_nonvol(i, K2) - write(unit=6,fmt="(a4,i3,a15,es10.2,2f10.3)")& - "NVOL", ispec,& - species(ispec)%name, xn(ispec,K2),-999.999, & - -999.999 - end do - - do ispec = S1,S2 - ! K = tabRTpL*COA/(mw*gamma) QUERY COA===!!!! - !Ksoa = tabRTpL(ispec,K2)*COA(K2)/(avg_mw(K2)) - Ksoa = 1.0/tabCiStar(ispec,itemp(K2)) !just for printout - write(unit=6,fmt="(a4,i3,a15,3es10.2,a4,es10.3,f13.4)") "SOA ",ispec,& - species(ispec)%name, xn(ispec,K2), & - tabCiStar(ispec,itemp(K2)),& !VBStabVpsoa(ispec,298), - Ksoa, " => ", Fpart(ispec,K2), 1000.0*ug_semivol(ispec, K2) - end do ! ispec - end if - - write(unit=6,fmt="(a,i2,f12.6)") "COA: ", iter, COA(K2) - - end if ! DEBUG - - end do ! ITER - - ! The above iteration has now given new values to: - ! - ! 1) COA(1:nz) - ! 3) Fgas(FIRST_SEMIVOL:LAST_SEMIVOL,1:nz) - ! 4) Fpart(FIRST_SEMIVOL:LAST_SEMIVOL,1:nz) - ! 5) ng(FIRST_SEMIVOL:LAST_SEMIVOL,1:nz) - ! - ! Note: ng(FIRST_NONVOLOC:LAST_NONVOLOC,1:nz) - ! and xn(1:LAST_SEMIVOL,1:nz) are unaffected by the - ! iteration. - !========================================================================= - - ! Set Fgas for later chemistry, and eset 3-D fields - - Fgas(S1:S2,:) = 1.0 - Fpart(S1:S2,:) - Grid_COA(i_pos,j_pos,:) = COA(:) - Fgas3d(S1:S2,i_pos,j_pos,:) = Fgas(S1:S2,:) !FEB2012 - -! PCM_F is for output only. Has MW 1 to avoid confusion with OC -! do not use ugC outputs, just ug - - xn(PART_OM_F,:) = COA(:) * ugC2xn * 12.0 - - !Grid_SOA_Fgas(S1:S2, i_pos,j_pos,:) = Fgas(S1:S2,:) - !VBS Grid_avg_mw(i_pos,j_pos,:) = avg_mw(:) - !SOA_gamma(i_pos,j_pos,:,:) = gamma(:,:) - - ! Outputs, ugC/m3 - - ! Would like to be able to store also total OM (not only OC) - ! at least for some components. And for a "total" OM and/or OM2.5 and OM10. - ! Total TC and EC (and TC2.5, TC10, EC2.5 and EC10) would also be useful. - ! Also perhaps the names of these species should reflect - ! that they are in units of C? - do k = K1, K2 - xn(PART_ASOA_OC,k) = sum ( Fpart(ASOA,k) *xn(ASOA,k) *species(ASOA)%carbons ) - xn(GAS_ASOA_OC,k) = sum ( Fgas(ASOA,k) * xn(ASOA,k) *species(ASOA)%carbons ) - xn(PART_BSOA_OC,k) = sum ( Fpart(BSOA,k) *xn(BSOA,k) *species(BSOA)%carbons ) - xn(GAS_BSOA_OC,k) = sum ( Fgas( BSOA,k) *xn(BSOA,k) *species(BSOA)%carbons ) - xn(NONVOL_FFUELOC25,k) = sum ( xn(NVFFUELOC25,k) *species(NVFFUELOC25)%carbons) - xn(NONV_FFUELOC_COARSE,k) = sum ( xn(NVFFUELOCCO,k) *species(NVFFUELOCCO)%carbons) - xn(NONVOL_WOODOC25,k) = sum ( xn(NVWOODOC25,k) *species(NVWOODOC25)%carbons ) - xn(NONVOL_FFIREOC25,k) = sum ( xn(NVFFIREOC25,k) *species(NVFFIREOC25)%carbons ) -! want to have PART_FFUELOA25/FFIREOA/WOODOA_OC working also with nonvolatile POA emissions, test this hard coded version first - xn(PART_FFUELOA25_OC,k) = sum ( Fpart(SVFFUELOA25,k) *xn(SVFFUELOA25,k) *species(SVFFUELOA25)%carbons ) + & - xn(NONVOL_FFUELOC25,k) - xn(PART_WOODOA25_OC,k) = sum ( Fpart(SVWOODOA25,k) *xn(SVWOODOA25,k) *species(SVWOODOA25)%carbons ) + & - xn(NONVOL_WOODOC25,k) - xn(PART_FFIREOA25_OC,k) = sum ( Fpart(SVFFIREOA25,k) *xn(SVFFIREOA25,k) *species(SVFFIREOA25)%carbons ) + & - xn(NONVOL_FFIREOC25,k) -!Test for storing in ug/m3, Use with caution! - xn(PART_ASOA_OM,k) = sum ( ug_semivol(ASOA,k) ) * ugC2xn * 12.0 - xn(PART_BSOA_OM,k) = sum ( ug_semivol(BSOA,k) ) * ugC2xn * 12.0 -! want to have PART_FFUELOA25/FFIREOA/WOODOA_OM working also with nonvolatile POA emissions, test this hard coded version first - xn(PART_FFUELOA25_OM,k) = sum ( ug_semivol(SVFFUELOA25,k) ) * ugC2xn * 12.0 + xn(NONVOL_FFUELOC25,k) * 1.25 * 12.0 ! factor 12.0 from M(OC25-components)=12 and M(OM-components)=1, OM/OC=1.25 assumed for Primary FFUELOC emissions - xn(PART_WOODOA25_OM,k) = sum ( ug_semivol(SVWOODOA25,k) ) * ugC2xn * 12.0 + xn(NONVOL_WOODOC25,k) * 1.7 * 12.0 ! OM/OC=1.7 assumed for Primary WOODOC and FFIRE emissions - xn(PART_FFIREOA25_OM,k) = sum ( ug_semivol(SVFFIREOA25,k) ) * ugC2xn * 12.0 + xn(NONVOL_FFIREOC25,k) * 1.7 * 12.0 - -!HARDCODE -! xn(AER_TBSOA,k) = xn(AER_BSOA,k) ! Just in case TBSOA is wanted for kam -!............................................................................... -!Research xn(PART_TBSOA_OC,k) = & -!Research Fpart( TERP_ng100,k) *xn(TERP_ng100,k) *species(TERP_ng100)%carbons + & -!Research Fpart( TERP_ug1,k) *xn(TERP_ug1,k) *species(TERP_ug1)%carbons + & -!Research Fpart( TERP_ug10,k) *xn(TERP_ug10,k) *species(TERP_ug10)%carbons + & -!Research Fpart( TERP_ug1e2,k) *xn(TERP_ug1e2,k) *species(TERP_ug1e2)%carbons + & -!Research Fpart( TERP_ug1e3,k) *xn(TERP_ug1e3,k) *species(TERP_ug1e3)%carbons -!Research xn(PART_IBSOA_OC,k) = & -!Research Fpart( ISOP_ng100,k) *xn(ISOP_ng100,k) *species(ISOP_ng100)%carbons + & -!Research Fpart( ISOP_ug1,k) *xn(ISOP_ug1,k) *species(ISOP_ug1)%carbons + & -!Research Fpart( ISOP_ug10,k) *xn(ISOP_ug10,k) *species(ISOP_ug10)%carbons + & -!Research Fpart( ISOP_ug1e2,k) *xn(ISOP_ug1e2,k) *species(ISOP_ug1e2)%carbons + & -!Research Fpart( ISOP_ug1e3,k) *xn(ISOP_ug1e3,k) *species(ISOP_ug1e3)%carbons -! xn(PART_SBSOA,k) = & -! Fpart( SESQ_ug1,k) *xn(SESQ_ug1,k) *species(SESQ_ug1)%carbons + & -! Fpart( SESQ_ug10,k) *xn(SESQ_ug10,k) *species(SESQ_ug10)%carbons + & -! Fpart( SESQ_ug1e2,k) *xn(SESQ_ug1e2,k) *species(SESQ_ug1e2)%carbons + & -! Fpart( SESQ_ug1e3,k) *xn(SESQ_ug1e3,k) *species(SESQ_ug1e3)%carbons -!............................................................................... - end do - xn(NONVOL_BGNDOC,:) = ugC2xn * BGND_OC(:) ! FAKE FOR NOW, 0.5 ug/m3 at surface - - ! for convencience: -! WARNING! Test changing to WOODOC25, FFIREOC25 and NONVOLOC25 may cause problems here, especially if coarse components are inlcuded later! But this is rather hard coded anyway... - xn(PART_OC10,:) = xn(PART_ASOA_OC,:)+xn(PART_BSOA_OC,:)+xn(NONV_FFUELOC_COARSE,:) + & - xn(NONVOL_BGNDOC,:)+ & - xn(PART_FFUELOA25_OC,:)+xn(PART_WOODOA25_OC,:)+ & - xn(PART_FFIREOA25_OC,:) -!WARNING! The below will NOT work for NPNA (nonvolatile) type runs. These include fine and coarse OC in NONVOL_FFUELOC. So for these runs the FFUELOC-contribution has to be added separately!!! -!Test shifting to NONVOL_FFUELOC25! - xn(PART_OC25,:) = xn(PART_ASOA_OC,:)+ xn(PART_BSOA_OC,:)+ & - xn(NONVOL_BGNDOC,:)+ & - xn(PART_FFUELOA25_OC,:)+xn(PART_WOODOA25_OC,:)+ & - xn(PART_FFIREOA25_OC,:) - - surfASOA = xn2ugC* xn(PART_ASOA_OC,K2) ! sum ( Fpart(ASOA,K2) * xn(ASOA,K2)*species(ASOA)%carbons ) - surfBSOA = xn2ugC* xn(PART_BSOA_OC,K2) ! sum ( Fpart(BSOA,K2) * xn(BSOA,K2)*species(BSOA)%carbons ) -! - surfFFUELOA25_OC = xn2ugC* xn(PART_FFUELOA25_OC,K2) ! - surfWOODOA25_OC = xn2ugC* xn(PART_WOODOA25_OC,K2) ! - surfFFIREOA25_OC = xn2ugC* xn(PART_FFIREOA25_OC,K2) ! - - surfBGNDOC = BGND_OC(K2) - surfOC25 = surfASOA+surfBSOA+surfFFUELOA25_OC+surfWOODOA25_OC+surfBGNDOC+surfFFIREOA25_OC ! - - !/ Sum of Biogenics ----------------------- + subroutine Init_OrganicAerosol(i,j,debug_flag) + integer, intent(in) :: i,j + logical, intent(in) :: debug_flag ! for debugging purposes only - if ( debug_flag .and. seconds == 0 ) then + ! empty + + end subroutine Init_OrganicAerosol - k=20 - write(unit=6,fmt="(a,3i3,2f7.2,f5.2,11es9.2)")"xns ug ", & - nmonth, nday, nhour, & - xn(PART_OM_F,20)*xn2ugC, & - COA(20), surfOC25,surfBGNDOC, surfASOA, surfBSOA, surfFFUELOA25_OC, & - surfWOODOA25_OC, & - surfFFIREOA25_OC!, & -! COA(20), surfOC25,surfBGNDOC, surfASOA, surfBSOA,surfFFUELOA25_OC, & -! surfOFFUELOA25_OC,surfFFUELOC25, surfWOODOA25, surfOWOODOA25, surfWOODOC25, & -! surfFFIREOA25, surfOFFIREOA25!, & -!vbs xn2ugC*xn(PART_IBSOA,k), xn2ugC*xn(PART_TBSOA,k), xn2ugC*xn(PART_SBSOA,k) + subroutine OrganicAerosol(i,j,debug_flag) + integer, intent(in) :: i,j + logical, intent(in) :: debug_flag ! for debugging purposes only - endif + ! empty - end subroutine OrganicAerosol + end subroutine OrganicAerosol end module OrganicAerosol_ml !xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx diff --git a/Nest_ml.f90 b/Nest_ml.f90 old mode 100644 new mode 100755 index 90ed703..de70f17 --- a/Nest_ml.f90 +++ b/Nest_ml.f90 @@ -57,7 +57,7 @@ module Nest_ml ! ExternalBICs_ml should handle different for different external sources. ! Experiment specific information must be set on ExternalBICs namelists. ! So far coded for FORECAST and EnsClimRCA(?) work. -use ExternalBICs_ml, only: set_extbic, icbc, & +use ExternalBICs_ml, only: set_extbic, icbc, ICBC_FMT,& EXTERNAL_BIC_SET, EXTERNAL_BC, EXTERNAL_BIC_NAME, TOP_BC, & iw, ie, js, jn, kt, &! i West/East bnd; j North/South bnd; k Top filename_eta @@ -108,8 +108,7 @@ module Nest_ml public :: wrtxn logical, private, save :: mydebug = .false. -integer, public, parameter :: NHOURSAVE=3 !time between two saves. should be a fraction of 24 -integer, public, parameter :: NHOURREAD=1 !time between two reads. should be a fraction of 24 +integer, private, save :: NHOURSAVE,NHOURREAD ! write/read frequency !if(NHOURREAD',I3,'=',A24,'*',F7.2,2L2))") & + write(*,"((1X,A,I3,'->',"//ICBC_FMT//"))") & ('Nest: ADV_IC',n,adv_ic(n),n=1,size(adv_ic)),& ('Nest: ADV_BC',n,adv_bc(n),n=1,size(adv_bc)) endif diff --git a/NetCDF_ml.f90 b/NetCDF_ml.f90 old mode 100644 new mode 100755 index 074cf9b..784985b --- a/NetCDF_ml.f90 +++ b/NetCDF_ml.f90 @@ -55,7 +55,8 @@ module NetCDF_ml ,debug_proc, debug_li, debug_lj & ,yp_EMEP_official,fi_EMEP,GRIDWIDTH_M_EMEP& ,grid_north_pole_latitude& - ,grid_north_pole_longitude,GlobalPosition& + ,grid_north_pole_longitude,dx_rot,dx_roti,x1_rot,y1_rot& + ,GlobalPosition& ,glat_fdom,glon_fdom,ref_latitude& ,projection, sigma_mid,gb_stagg,gl_stagg,glon& ,sigma_bnd& @@ -282,7 +283,20 @@ subroutine CreatenetCDFfile(fileName,GIMAXcdf,GJMAXcdf,ISMBEGcdf,JSMBEGcdf,& call check(nf90_put_att(ncFileID, jVarID, "long_name", "latitude")) call check(nf90_put_att(ncFileID, jVarID, "units", "degrees_north")) - else !general projection + elseif(UsedProjection=='Rotated_Spherical')then + call check(nf90_def_dim(ncid = ncFileID, name = "i", len = GIMAXcdf, dimid = iDimID)) + call check(nf90_def_var(ncFileID, "i", nf90_float, dimids = iDimID, varID = iVarID) ) + call check(nf90_put_att(ncFileID, iVarID, "standard_name", "grid_longitude")) + call check(nf90_put_att(ncFileID, iVarID, "long_name", "Rotated longitude")) + call check(nf90_put_att(ncFileID, iVarID, "units", "degrees")) + call check(nf90_put_att(ncFileID, iVarID, "axis", "X")) + call check(nf90_def_dim(ncid = ncFileID, name = "j", len = GJMAXcdf, dimid = jDimID)) + call check(nf90_def_var(ncFileID, "j", nf90_float, dimids = jDimID, varID = jVarID) ) + call check(nf90_put_att(ncFileID, jVarID, "standard_name", "grid_latitude")) + call check(nf90_put_att(ncFileID, jVarID, "long_name", "Rotated latitude")) + call check(nf90_put_att(ncFileID, jVarID, "units", "degrees")) + call check(nf90_put_att(ncFileID, jVarID, "axis", "Y")) + else !general projection call check(nf90_def_dim(ncid = ncFileID, name = "i", len = GIMAXcdf, dimid = iDimID)) call check(nf90_def_dim(ncid = ncFileID, name = "j", len = GJMAXcdf, dimid = jDimID)) call check(nf90_def_var(ncFileID, "i", nf90_float, dimids = iDimID, varID = iVarID) ) @@ -398,7 +412,7 @@ subroutine CreatenetCDFfile(fileName,GIMAXcdf,GJMAXcdf,ISMBEGcdf,JSMBEGcdf,& call check(nf90_put_att(ncFileID, VarID, "latitude_of_projection_origin", 90.0)) call check(nf90_put_att(ncFileID, VarID, "scale_factor_at_projection_origin", scale_at_projection_origin)) elseif(UsedProjection=='lon lat')then - + !no additional attributes elseif(UsedProjection=='Rotated_Spherical')then call check(nf90_def_var(ncid = ncFileID, name = "Rotated_Spherical", xtype = nf90_int, varID=varID ) ) call check(nf90_put_att(ncFileID, VarID, "grid_mapping_name", "rotated_latitude_longitude")) @@ -466,6 +480,16 @@ subroutine CreatenetCDFfile(fileName,GIMAXcdf,GJMAXcdf,ISMBEGcdf,JSMBEGcdf,& glon_fdom(ISMBEGcdf:ISMBEGcdf+GIMAXcdf-1,JSMBEGcdf:JSMBEGcdf+GJMAXcdf-1))) endif + elseif(UsedProjection=='Rotated_Spherical')then + do i=1,GIMAXcdf + xcoord(i)= (i+ISMBEGcdf-1)*dx_rot+x1_rot + enddo + do j=1,GJMAXcdf + ycoord(j)= (j+JSMBEGcdf-1)*dx_rot+y1_rot + enddo + call check(nf90_put_var(ncFileID, iVarID, xcoord(1:GIMAXcdf)) ) + call check(nf90_put_var(ncFileID, jVarID, ycoord(1:GJMAXcdf)) ) + elseif(UsedProjection=='lon lat') then do i=1,GIMAXcdf xcoord(i)= glon_fdom(i+ISMBEGcdf-1,1) @@ -602,6 +626,20 @@ subroutine CreatenetCDFfile_Eta(fileName,GIMAXcdf,GJMAXcdf,ISMBEGcdf,JSMBEGcdf,& call check(nf90_put_att(ncFileID, jVarID, "long_name", "latitude")) call check(nf90_put_att(ncFileID, jVarID, "units", "degrees_north")) + elseif(UsedProjection=='Rotated_Spherical')then + call check(nf90_def_dim(ncid = ncFileID, name = "i", len = GIMAXcdf, dimid = iDimID)) + call check(nf90_def_var(ncFileID, "i", nf90_float, dimids = iDimID, varID = iVarID) ) + call check(nf90_put_att(ncFileID, iVarID, "standard_name", "grid_longitude")) + call check(nf90_put_att(ncFileID, iVarID, "long_name", "Rotated longitude")) + call check(nf90_put_att(ncFileID, iVarID, "units", "degrees")) + call check(nf90_put_att(ncFileID, iVarID, "axis", "X")) + call check(nf90_def_dim(ncid = ncFileID, name = "j", len = GJMAXcdf, dimid = jDimID)) + call check(nf90_def_var(ncFileID, "j", nf90_float, dimids = jDimID, varID = jVarID) ) + call check(nf90_put_att(ncFileID, jVarID, "standard_name", "grid_latitude")) + call check(nf90_put_att(ncFileID, jVarID, "long_name", "Rotated latitude")) + call check(nf90_put_att(ncFileID, jVarID, "units", "degrees")) + call check(nf90_put_att(ncFileID, jVarID, "axis", "Y")) + else !general projection call check(nf90_def_dim(ncid = ncFileID, name = "i", len = GIMAXcdf, dimid = iDimID)) call check(nf90_def_dim(ncid = ncFileID, name = "j", len = GJMAXcdf, dimid = jDimID)) @@ -809,6 +847,16 @@ subroutine CreatenetCDFfile_Eta(fileName,GIMAXcdf,GJMAXcdf,ISMBEGcdf,JSMBEGcdf,& glon_fdom(ISMBEGcdf:ISMBEGcdf+GIMAXcdf-1,JSMBEGcdf:JSMBEGcdf+GJMAXcdf-1))) endif + elseif(UsedProjection=='Rotated_Spherical')then + do i=1,GIMAXcdf + xcoord(i)= (i+ISMBEGcdf-1)*dx_rot+x1_rot + enddo + do j=1,GJMAXcdf + ycoord(j)= (j+JSMBEGcdf-1)*dx_rot+y1_rot + enddo + call check(nf90_put_var(ncFileID, iVarID, xcoord(1:GIMAXcdf)) ) + call check(nf90_put_var(ncFileID, jVarID, ycoord(1:GJMAXcdf)) ) + elseif(UsedProjection=='lon lat') then do i=1,GIMAXcdf xcoord(i)= glon_fdom(i+ISMBEGcdf-1,1) @@ -2252,6 +2300,7 @@ recursive subroutine ReadField_CDF(fileName,varname,Rvar,nstart,kstart,kend,inte if( present(known_projection) ) then data_projection = trim(known_projection) + if(trim(known_projection)=="longitude latitude")data_projection = "lon lat" if ( debug ) write(*,*) 'data known_projection ',trim(data_projection) else call check(nf90_get_att(ncFileID, nf90_global, "projection", data_projection ),"Proj") @@ -2259,7 +2308,8 @@ recursive subroutine ReadField_CDF(fileName,varname,Rvar,nstart,kstart,kend,inte end if if(MasterProc)write(*,*)'Interpolating ',trim(varname),' from ',trim(filename),' to present grid' - if(trim(data_projection)=="lon lat")then ! here we have simple 1-D lat, lon + if(trim(data_projection)=="lon lat")then + ! here we have simple 1-D lat, lon allocate(Rlon(dims(1)), stat=alloc_err) allocate(Rlat(dims(2)), stat=alloc_err) if ( debug ) write(*,"(a,a,i5,i5,a,i5)") 'alloc_err lon lat ',& @@ -2273,13 +2323,19 @@ recursive subroutine ReadField_CDF(fileName,varname,Rvar,nstart,kstart,kend,inte status=nf90_inq_varid(ncid = ncFileID, name = 'lon', varID = lonVarID) if(status /= nf90_noerr) then status=nf90_inq_varid(ncid = ncFileID, name = 'LON', varID = lonVarID) - call CheckStop(status /= nf90_noerr,'did not find longitude variable') + if(status /= nf90_noerr) then + status=nf90_inq_varid(ncid = ncFileID, name = 'longitude', varID = lonVarID) + call CheckStop(status /= nf90_noerr,'did not find longitude variable') + endif endif status=nf90_inq_varid(ncid = ncFileID, name = 'lat', varID = latVarID) if(status /= nf90_noerr) then status=nf90_inq_varid(ncid = ncFileID, name = 'LAT', varID = latVarID) - call CheckStop(status /= nf90_noerr,'did not find latitude variable') + if(status /= nf90_noerr) then + status=nf90_inq_varid(ncid = ncFileID, name = 'latitude', varID = latVarID) + call CheckStop(status /= nf90_noerr,'did not find latitude variable') + endif endif if(trim(data_projection)=="lon lat")then call check(nf90_get_var(ncFileID, lonVarID, Rlon), 'Getting Rlon') @@ -2310,9 +2366,9 @@ recursive subroutine ReadField_CDF(fileName,varname,Rvar,nstart,kstart,kend,inte !check that there are dimensions called lon and lat call check(nf90_inquire_dimension(ncid = ncFileID, dimID = dimids(1), name=name ),name) - call CheckStop(trim(name)/='lon',"longitude not found") + call CheckStop(trim(name)/='lon'.and.trim(name)/='longitude',"longitude not found") call check(nf90_inquire_dimension(ncid = ncFileID, dimID = dimids(2), name=name ),name) - call CheckStop(trim(name)/='lat',"latitude not found") + call CheckStop(trim(name)/='lat'.and.trim(name)/='latitude',"latitude not found") if(data3D)then call check(nf90_inquire_dimension(ncid = ncFileID, dimID = dimids(3), name=name )) @@ -2530,7 +2586,7 @@ recursive subroutine ReadField_CDF(fileName,varname,Rvar,nstart,kstart,kend,inte Ndiv=max(1,Ndiv) Ndiv2=Ndiv*Ndiv ! - if(projection/='Stereographic'.and.projection/='lon lat')then + if(projection/='Stereographic'.and.projection/='lon lat'.and.projection/='Rotated_Spherical')then !the method should be revised or used only occasionally if(me==0)write(*,*)'WARNING: interpolation method may be CPU demanding' endif @@ -2808,9 +2864,9 @@ recursive subroutine ReadField_CDF(fileName,varname,Rvar,nstart,kstart,kend,inte yp_ext_div=(yp_ext+0.5)*Ndiv-0.5 an_ext_div=an_ext*Ndiv - if(projection/='Stereographic'.and.projection/='lon lat'.and.projection=='Rotated_Spherical')then + if(projection/='Stereographic'.and.projection/='lon lat'.and.projection/='Rotated_Spherical')then !the method should be revised or used only occasionally - if(me==0)write(*,*)'WARNING: interpolation method may be CPU demanding' + if(me==0)write(*,*)'WARNING: interpolation method may be CPU demanding:',projection endif k2=1 if(data3D)k2=kend-kstart+1 @@ -2923,9 +2979,9 @@ recursive subroutine ReadField_CDF(fileName,varname,Rvar,nstart,kstart,kend,inte an_ext_div=an_ext*Ndiv if(MasterProc.and.debug)write(*,*)'zero_order interpolation ',an_ext_div,xp_ext_div,yp_ext_div,dims(1),dims(2) - if(projection/='Stereographic'.and.projection/='lon lat')then + if(projection/='Stereographic'.and.projection/='lon lat'.and.projection/='Rotated_Spherical')then !the method should be revised or used only occasionally - if(me==0)write(*,*)'WARNING: interpolation method may be CPU demanding' + if(me==0)write(*,*)'WARNING: interpolation method may be CPU demanding',projection endif diff --git a/Output_hourly.f90 b/Output_hourly.f90 old mode 100644 new mode 100755 index eb786bb..18a6c35 --- a/Output_hourly.f90 +++ b/Output_hourly.f90 @@ -32,16 +32,31 @@ subroutine hourly_out() !! spec,ofmt,ix1,ix2,iy1,iy2,unitfac) ! Calculates and ! Outputs hourly concentration (or met) values for a sub-set of the grid. ! -!** REVISION HISTORY: -! Extended to produce new file, Hourly.mmyy, every month, 10/5/01 ds -! stop_test used instead of stop_all, su, 05/01 -! Extended for variable format, met, xn_adv or xn_shl, ds, and to use -! Asc2D type 19/4/01 -! Corrected for IRUNBEG, etc., su, 4/01 -! New, ds, 5/3/99 +!** Surface output +! Onyl relevant for lowermost model level, meaningless for higher levels. +! - ADVppbv: instantaneous surface concentrations in ppb. +! - ADVugXX: instantaneous surface concentrations in ug (ug/m3, ugC/m3, ugS/m3, ugN/m3). +! For ug/m3 output, set hr_out%unitconv=to_ug_ADV(ixadv). +! For ugX/m3 output, set hr_out%unitconv=to_ug_X(ixadv). +! - D2D_mean: hourly means comparable to daily output defined on My_Derived_ml. +! - D2D_inst: instantaneous value used in daily output defined on My_Derived_ml. +! - D2D_accum: accumulated alue used in daily output defined on My_Derived_ml. +! - D2D: D2D_mean/D2D_accum according to the corresponding My_Derived_ml definition. ! -!************************************************************************* +!** Column integrated output +! - COLUMN: +! For ug/m2 output, set hr_out%unitconv=to_ug_ADV(ixadv). +! For ugX/m2 output, set hr_out%unitconv=to_ug_X(ixadv). +! For molec/cm2 output, set hr_out%unitconv=to_molec_cm2. ! +!** Multi-layer output +! NLEVELS_HOURLY (My_Outputs_ml) max levels +! hr_out%nk levels for each output +! - BCVppbv: instantaneous grid-centre concentrations in ppb. +! - BCVugXX: instantaneous grid-centre concentrations in ug (ug/m3, ugC/m3, ugS/m3, ugN/m3). +! For ug/m3 output, set hr_out%unitconv=to_ug_ADV(ixadv). +! For ugX/m3 output, set hr_out%unitconv=to_ug_X(ixadv). +!************************************************************************* use My_Outputs_ml, only: NHOURLY_OUT, & ! No. outputs NLEVELS_HOURLY, & ! No. output levels hr_out, & ! Required outputs @@ -51,7 +66,7 @@ subroutine hourly_out() !! spec,ofmt,ix1,ix2,iy1,iy2,unitfac) use CheckStop_ml, only: CheckStop use Chemfields_ml, only: xn_adv,xn_shl,cfac,PM25_water,PM25_water_rh50,AOD use ChemGroups_ml, only: chemgroups - use Derived_ml, only: num_deriv2d ! D2D houtly output type + use Derived_ml, only: num_deriv2d,nav_2d ! D2D houtly output type use DerivedFields_ml, only: f_2d,d_2d ! D2D houtly output type use OwnDataTypes_ml, only: Asc2D, Deriv use ChemSpecs_shl_ml ,only: NSPEC_SHL ! Maps indices @@ -60,7 +75,7 @@ subroutine hourly_out() !! spec,ofmt,ix1,ix2,iy1,iy2,unitfac) debug_proc, debug_li,debug_lj use Io_ml, only: IO_HOURLY use ModelConstants_ml,only: KMAX_MID, MasterProc, & - IOU_INST, IOU_HOUR, IOU_YEAR, IOU_HOUR_PREVIOUS, & + IOU_INST, IOU_HOUR, IOU_YEAR, IOU_YEAR_LASTHH, & DEBUG => DEBUG_OUT_HOUR,runlabel1,HOURLYFILE_ending,& FORECAST use ModelConstants_ml,only: SELECT_LEVELS_HOURLY !NML @@ -115,11 +130,17 @@ subroutine hourly_out() !! spec,ofmt,ix1,ix2,iy1,iy2,unitfac) integer, pointer, dimension(:) :: gspec=>null() ! group array of indexes real, pointer, dimension(:) :: gunit_conv=>null() ! & unit conv. factors + integer, allocatable, dimension(:), save :: navg ! D2D average counter + character(len=len(fileName_hour)) :: filename + logical :: file_exist=.false. logical, save :: debug_flag ! = ( MasterProc .and. DEBUG ) logical :: surf_corrected ! to get 3m values - logical :: file_exist=.false. + character(len=*), parameter :: & + SRF_TYPE(9)=(/"ADVppbv ","ADVugXX ","ADVugXXgroup",& + "COLUMN ","COLUMNgroup ","D2D ",& + "D2D_mean ","D2D_inst ","D2D_accum "/) if(NHOURLY_OUT<= 0) then if(my_first_call.and.MasterProc.and.DEBUG) & @@ -141,6 +162,9 @@ subroutine hourly_out() !! spec,ofmt,ix1,ix2,iy1,iy2,unitfac) ! hr_out(ih)%nk = min(KMAX_MID,hr_out(ih)%nk) if(debug_flag) write(*,*) "DEBUG Hourly nk ", ih, hr_out(ih)%nk enddo ! ih + + allocate(navg(NHOURLY_OUT)) ! allocate and initialize + navg(:)=0.0 ! D2D average counter endif ! first_call filename=trim(runlabel1)//date2string(HOURLYFILE_ending,current_date) @@ -161,8 +185,7 @@ subroutine hourly_out() !! spec,ofmt,ix1,ix2,iy1,iy2,unitfac) nk = hr_out(ih)%nk CDFtype=Real4 ! can be choosen as Int1,Int2,Int4,Real4 or Real8 scale=1. - if(any(hr_out(ih)%type==(/"ADVppbv ","ADVugXX ","ADVugXXgroup",& - "COLUMN ","COLUMNgroup ","D2D "/)))nk=1 + if(any(hr_out(ih)%type==SRF_TYPE))nk=1 select case(nk) case(1) ! write as 2D call Out_netCDF(IOU_HOUR,def1,2,1,hourly,scale,CDFtype,ist,jst,ien,jen,& @@ -188,12 +211,11 @@ subroutine hourly_out() !! spec,ofmt,ix1,ix2,iy1,iy2,unitfac) hourly(:,:) = 0.0 ! initialize - HLOOP: do ih = 1, NHOURLY_OUT + HLOOP: do ih = 0, NHOURLY_OUT hr_out_type=trim(hr_out(ih)%type) hr_out_nk=hr_out(ih)%nk - if(any(hr_out_type==(/"ADVppbv ","ADVugXX ","ADVugXXgroup",& - "COLUMN ","COLUMNgroup ","D2D "/)))hr_out_nk=1 + if(any(hr_out_type==SRF_TYPE))hr_out_nk=1 KVLOOP: do k = 1,hr_out_nk @@ -222,42 +244,30 @@ subroutine hourly_out() !! spec,ofmt,ix1,ix2,iy1,iy2,unitfac) if(ik==0)then ik=KMAX_MID ! surface/lowermost level if(debug_flag) write(*,*)"DEBUG LOWEST LEVELS", ik, hr_out_type - if(any(hr_out_type==(/"BCVppbv ","BCVugXX ",& - "BCVugXXgroup"/)))& -! "BCVugXXgroup","Out3D "/)))& + select case(hr_out_type) + case("BCVppbv","BCVugXX","BCVugXXgroup") hr_out_type(1:3)="ADV" ! ensure surface output - if(any(hr_out_type==(/"PMwater"/)))& + case("PMwater") hr_out_type=trim(hr_out_type)//"SRF" + case("D2D") + if(f_2d(ispec)%avg)then ! averaged variables +! hr_out_type="D2D_inst" ! output instantaneous values + hr_out_type="D2D_mean" ! output mean values + else ! accumulated variables + hr_out_type="D2D_accum" + endif + endselect else !TESTHH QUERY: ik=KMAX_MID-ik+1 ! model level to be outputed - if(any(hr_out_type==(/"ADVppbv ","ADVugXX ","ADVugXXgroup"/)))& + select case(hr_out_type) + case("ADVppbv","ADVugXX","ADVugXXgroup") ik=KMAX_MID ! all ADV* types represent surface output + endselect endif endif endif - !---------------------------------------------------------------- - ! Multi-layer output. - ! Specify NLEVELS_HOURLY here, and in hr_out defs use either: - ! ADVppbv to get surface concentrations (onyl relevant for - ! layer k=20 of course - gives meaningless number - ! for higher levels). - ! Or, - ! BCVppbv to get grid-centre concentrations (relevant for all layers). - ! - ! For ug output (ug/m3, ugC/m3, ugS/m3, ugN/m3) use - ! ADVugXX to get surface concentrations (only lowermost model level). - ! Or, - ! BCVppbv to get grid-centre concentrations (model levels). - ! For ug/m3 output, set hr_out%unitconv=to_ug_ADV(ixadv). - ! For ugX/m3 output, set hr_out%unitconv=to_ug_X(ixadv). - ! - ! For Column integrated output use COLUMN - ! For ug/m2 output, set hr_out%unitconv=to_ug_ADV(ixadv). - ! For ugX/m2 output, set hr_out%unitconv=to_ug_X(ixadv). - ! For molec/cm2 output, set hr_out%unitconv=to_molec_cm2. - !---------------------------------------------------------------- if(debug_flag) write(*,"(5a,i4)") "DEBUG Hourly MULTI ",& trim(hr_out(ih)%name), " case ", trim(hr_out_type), " k: ", ik @@ -469,34 +479,79 @@ subroutine hourly_out() !! spec,ofmt,ix1,ix2,iy1,iy2,unitfac) case("Idiffus") ! Diffuse radiation (W/m2); Skip Units conv. forall(i=1:limax,j=1:ljmax) hourly(i,j) = Idiffuse(i,j) - case ( "D2D" ) + case("D2D_inst") + ! Here ispec is the index in the f_2d arrays + call CheckStop(ispec<1.or.ispec>num_deriv2d,& + "ERROR-DEF! Hourly_out: "//trim(hr_out(ih)%name)//", wrong D2D id!") + if(hr_out(ih)%unit=="") hr_out(ih)%unit = f_2d(ispec)%unit + unit_conv = hr_out(ih)%unitconv*f_2d(ispec)%scale + + if(debug_flag) then + i=debug_li + j=debug_lj + write(*,"(2a,2i4,a,3g12.3)") "OUTHOUR "//trim(hr_out_type),& + trim(hr_out(ih)%name), ih, ispec,trim(f_2d(ispec)%name),& + d_2d(ispec,i,j,IOU_YEAR), d_2d(ispec,i,j,IOU_YEAR_LASTHH),& + unit_conv + endif + forall(i=1:limax,j=1:ljmax) + hourly(i,j) = d_2d(ispec,i,j,IOU_INST) * unit_conv + endforall + if(debug_flag) & + write(*,'(a,2i3,2es12.3)')"HHH DEBUG D2D", ispec, ih, & + hr_out(ih)%unitconv, hourly(debug_li,debug_lj) + + case("D2D_accum") ! Here ispec is the index in the f_2d arrays call CheckStop(ispec<1.or.ispec>num_deriv2d,& "ERROR-DEF! Hourly_out: "//trim(hr_out(ih)%name)//", wrong D2D id!") if(hr_out(ih)%unit=="") hr_out(ih)%unit = f_2d(ispec)%unit unit_conv = hr_out(ih)%unitconv*f_2d(ispec)%scale - if(f_2d(ispec)%avg)then ! non accumulated variables - if( debug_flag ) write(*,*) " D2Davg ",& - trim(hr_out(ih)%name), ih, ispec, trim(f_2d(ispec)%name), f_2d(ispec)%avg - forall(i=1:limax,j=1:ljmax) - hourly(i,j) = d_2d(ispec,i,j,IOU_INST) * unit_conv - endforall - else ! hourly accumulated variables - if(debug_flag) then - i=debug_li - j=debug_lj - write(*,"(2a,2i4,a,3g12.3)") "OUTHOUR D2Dpre ",& - trim(hr_out(ih)%name), ih, ispec,trim(f_2d(ispec)%name),& - d_2d(ispec,i,j,IOU_YEAR), d_2d(ispec,i,j,IOU_HOUR_PREVIOUS),& - unit_conv - endif - forall(i=1:limax,j=1:ljmax) - hourly(i,j) = (d_2d(ispec,i,j,IOU_YEAR)& - -d_2d(ispec,i,j,IOU_HOUR_PREVIOUS)) * unit_conv - d_2d(ispec,i,j,IOU_HOUR_PREVIOUS)=d_2d(ispec,i,j,IOU_YEAR) - endforall + if(debug_flag) then + i=debug_li + j=debug_lj + write(*,"(2a,2i4,a,3g12.3)") "OUTHOUR "//trim(hr_out_type),& + trim(hr_out(ih)%name), ih, ispec,trim(f_2d(ispec)%name),& + d_2d(ispec,i,j,IOU_YEAR), d_2d(ispec,i,j,IOU_YEAR_LASTHH),& + unit_conv + endif + forall(i=1:limax,j=1:ljmax) + hourly(i,j) = (d_2d(ispec,i,j,IOU_YEAR)& + -d_2d(ispec,i,j,IOU_YEAR_LASTHH)) * unit_conv + d_2d(ispec,i,j,IOU_YEAR_LASTHH)=d_2d(ispec,i,j,IOU_YEAR) + endforall + if(debug_flag) & + write(*,'(a,2i3,2es12.3)')"HHH DEBUG D2D", ispec, ih, & + hr_out(ih)%unitconv, hourly(debug_li,debug_lj) + + case("D2D_mean") + ! Here ispec is the index in the f_2d arrays + call CheckStop(ispec<1.or.ispec>num_deriv2d,& + "ERROR-DEF! Hourly_out: "//trim(hr_out(ih)%name)//", wrong D2D id!") + if(hr_out(ih)%unit=="") hr_out(ih)%unit = f_2d(ispec)%unit + unit_conv = hr_out(ih)%unitconv*f_2d(ispec)%scale & + /MAX(nav_2d(ispec,IOU_YEAR)-navg(ih),1) + navg(ih)=nav_2d(ispec,IOU_YEAR) + + if(debug_flag) then + if(f_2d(ispec)%avg)then ! averaged variables + write(*,"(a,1x,$)") "OUTHOUR D2D_mean avg" + else ! accumulated variables --> mean + write(*,"(a,1x,$)") "OUTHOUR D2D_mean acc" + endif + i=debug_li + j=debug_lj + write(*,"(a,2i4,a,3g12.3)"),& + trim(hr_out(ih)%name), ih, ispec,trim(f_2d(ispec)%name),& + d_2d(ispec,i,j,IOU_YEAR), d_2d(ispec,i,j,IOU_YEAR_LASTHH),& + unit_conv endif + forall(i=1:limax,j=1:ljmax) + hourly(i,j) = (d_2d(ispec,i,j,IOU_YEAR)& + -d_2d(ispec,i,j,IOU_YEAR_LASTHH)) * unit_conv + d_2d(ispec,i,j,IOU_YEAR_LASTHH)=d_2d(ispec,i,j,IOU_YEAR) + endforall if(debug_flag) & write(*,'(a,2i3,2es12.3)')"HHH DEBUG D2D", ispec, ih, & hr_out(ih)%unitconv, hourly(debug_li,debug_lj) @@ -556,18 +611,8 @@ subroutine hourly_out() !! spec,ofmt,ix1,ix2,iy1,iy2,unitfac) call Out_netCDF(IOU_HOUR,def1,3,1,hourly,scale,CDFtype,ist,jst,ien,jen,klevel) !case default ! no output endselect - enddo KVLOOP - + enddo KVLOOP enddo HLOOP - !write out surface pressure at each time step to define vertical coordinates - if(NHOURLY_OUT>0)then - def1%name='PS' - def1%unit='hPa' - def1%class='Surface pressure' - CDFtype=Real4 ! can be choosen as Int1,Int2,Int4,Real4 or Real8 - scale=1. - call Out_netCDF(IOU_HOUR,def1,2,1,ps(:,:,1)*0.01,scale,CDFtype,ist,jst,ien,jen) - endif !Not closing seems to give a segmentation fault when opening the daily file !Probably just a bug in the netcdf4/hdf5 library. diff --git a/PhyChem_ml.f90 b/PhyChem_ml.f90 old mode 100644 new mode 100755 diff --git a/PhysicalConstants_ml.f90 b/PhysicalConstants_ml.f90 old mode 100644 new mode 100755 diff --git a/Radiation_ml.f90 b/Radiation_ml.f90 old mode 100644 new mode 100755 diff --git a/ReadField_ml.f90 b/ReadField_ml.f90 old mode 100644 new mode 100755 diff --git a/Rsurface_ml.f90 b/Rsurface_ml.f90 old mode 100644 new mode 100755 diff --git a/Runchem_ml.f90 b/Runchem_ml.f90 old mode 100644 new mode 100755 diff --git a/SeaSalt_ml.f90 b/SeaSalt_ml.f90 old mode 100644 new mode 100755 index 224acfe..ed527ac --- a/SeaSalt_ml.f90 +++ b/SeaSalt_ml.f90 @@ -264,7 +264,7 @@ subroutine SeaSalt_flux (i,j, debug_flag) !.. Fine particles emission [molec/cm3/s] need to be scaled to get units kg/m2/s consistent with ! Emissions_ml (snapemis). Scaling factor is do ii = 1, NFIN - rcss( iSSFI) = rcss( iSSFI) + & + rcss( iSSFI) = rcss( iSSFI) & !! ESX SS_prod(QSSFI,i,j) = SS_prod(QSSFI,i,j) & + ss_flux(ii) * d3(ii) * n2m & * water_fraction(i,j) @@ -274,7 +274,7 @@ subroutine SeaSalt_flux (i,j, debug_flag) !..Coarse particles emission [molec/cm3/s] do ii = NFIN+1, NFIN+NCOA - rcss( iSSCO ) = rcss( iSSCO ) + & + rcss( iSSCO ) = rcss( iSSCO ) & !!ESX SS_prod(QSSCO,i,j) = SS_prod(QSSCO,i,j) & + ss_flux(ii) * d3(ii) * n2m & * water_fraction(i,j) diff --git a/Setup_1d_ml.f90 b/Setup_1d_ml.f90 old mode 100644 new mode 100755 diff --git a/Setup_1dfields_ml.f90 b/Setup_1dfields_ml.f90 old mode 100644 new mode 100755 diff --git a/Sites_ml.f90 b/Sites_ml.f90 old mode 100644 new mode 100755 diff --git a/SoilWater_ml.f90 b/SoilWater_ml.f90 old mode 100644 new mode 100755 diff --git a/Solver.f90 b/Solver.f90 old mode 100644 new mode 100755 diff --git a/StoFlux_ml.f90 b/StoFlux_ml.f90 old mode 100644 new mode 100755 diff --git a/SubMet_ml.f90 b/SubMet_ml.f90 old mode 100644 new mode 100755 diff --git a/Tabulations_ml.f90 b/Tabulations_ml.f90 old mode 100644 new mode 100755 diff --git a/TimeDate_ExtraUtil_ml.f90 b/TimeDate_ExtraUtil_ml.f90 old mode 100644 new mode 100755 diff --git a/TimeDate_ml.f90 b/TimeDate_ml.f90 old mode 100644 new mode 100755 diff --git a/Timefactors_ml.f90 b/Timefactors_ml.f90 old mode 100644 new mode 100755 diff --git a/Timing_ml.f90 b/Timing_ml.f90 old mode 100644 new mode 100755 diff --git a/Unimod.f90 b/Unimod.f90 old mode 100644 new mode 100755 diff --git a/Units_ml.f90 b/Units_ml.f90 index b810dd3..f3956ce 100644 --- a/Units_ml.f90 +++ b/Units_ml.f90 @@ -84,9 +84,9 @@ module Units_ml umap("mgN","mgN/m2",mgNm2),& umap("mgS","mgS/m2",mgSm2),& ! Exposure to radioactive material - umap("uBq" ,"uBq/m3" ,ugXm3 ),& ! inst/mean exposure - umap("uBqh","uBq h/m3",ugXm3*s2h),& ! accumulated exposure - umap("mBq" ,"mBq/m2" ,mgXm2 ),& ! deposition + umap("uBq" ,"uBq/m3" ,ugXm3),& ! inst/mean exposure + umap("uBqh","uBq h/m3",ugXm3),& ! accumulated exposure over 1 hour + umap("mBq" ,"mBq/m2" ,mgXm2),& ! deposition ! Aerosol optical properties umap("ext" ,"ext550nm",extX),&! ext* units need to be further multiplied... ! Coulumn output @@ -209,9 +209,11 @@ function Units_Scale(txtin,iadv,unitstxt,volunit,needroa,debug_msg) result(units case("ugSS","ugSS/m3","ugP","ugP/m3",& "mgSS","mgSS/m2","mgP","mgP/m2") txt=txt(1:2) + case("micro g/m3") + txt="ug" case("mol/mol","mole mole-1","mixratio") txt="mix_ratio" - case("ppbv") + case("ppbv","ppbV") txt="ppb" endselect i=find_index(txt,unit_map(:)%utxt) diff --git a/Volcanos_ml.f90 b/Volcanos_ml.f90 old mode 100644 new mode 100755 diff --git a/Wesely_ml.f90 b/Wesely_ml.f90 old mode 100644 new mode 100755 diff --git a/config_emep.nml b/config_emep.nml index 5180b70..42e7a48 100644 --- a/config_emep.nml +++ b/config_emep.nml @@ -10,10 +10,15 @@ USE_SOILWATER = T, ! Uses SMI from meteo data USE_CONVECTION = F, USE_DEGREEDAY_FACTORS = T, - USE_FOREST_FIRES = T, + USE_FOREST_FIRES = F, ! See the user guide dor details USE_SEASALT = T, + USE_AIRCRAFT_EMIS = F, ! Needs global file, see user guide + USE_LIGHTNING_EMIS = T, ! +! USE_EURO_SOILNOX = T, ! diff for global + Euro runs + EURO_SOILNOX_DEPSCALE = 1.0, ! Scale, default 1 (see ModelConstants_ml) USE_GLOBAL_SOILNOX = F, ! diff for global + Euro runs +! USE_POLLEN = F, ! EXPERIMENTAL. Only works if start Jan 1 DO_SAHARA = T, USE_ROADDUST = T, ! Only EECCA? @@ -33,14 +38,14 @@ USE_EMERGENCY = F, ! Used for FORECASTs usually, EMEP2010 ANALYSIS = F, ! EXPERIMENTAL: 3DVar data assimilation USE_AOD = F, ! Used for FORECASTs usually -!-- +! BGND_CH4 = -1, ! Reset CH4 (use ppb). (Use -1 for defaults) !_- !--- "fake" vegetation for ozone POD calculations FLUX_VEGS = "IAM_CR" "IAM_DF" "IAM_MF" ! do not need all array elements !FLUX_VEGS = "IAM_CR" "IAM_DF" "IAM_MF" "CCE_SPRUCE" "CCE_BEECH" "ACE_PINE" "ACE_OAK" "ACE_BEECH" "NEUR_SPRUCE" "NEUR_BIRCH" "MED_OAK", "MED_PINE" "MED_BEECH" ! Do not need to fill array :-) !-------- Sub domain (x0,x1,y0,y1) -! RUNDOMAIN = 1, 100, 1, 100, ! Orig EMEP domain in EECCA (for benchmarks) + RUNDOMAIN = 1, 100, 1, 100, ! Orig EMEP domain in EECCA (for benchmarks) ! RUNDOMAIN = 40, 210, 12, 184, ! SR TNO28 area ! RUNDOMAIN = 240, 720, 48, 736, ! TNO07 reduced (15W-45E;30N-73N) ! RUNDOMAIN = 120, 360, 24, 368, ! TNO14 reduced (15W-45E;30N-73N) @@ -61,6 +66,8 @@ FLUX_VEGS = "IAM_CR" "IAM_DF" "IAM_MF" ! do not need all array elements ! 10=write at end of run; 11=read at start; 12=read at start and write at end (BIC) !------------------------------ MODE = 0, +! NHOURSAVE = 3, ! hours between saves. Fraction of 24 +! NHOURREAD = 1, ! hours between reads. Fraction of 24 !-------- File name templates for Nest I/O template_read_3D = 'EMEP_IN.nc', ! a different path can be set here template_read_BC = 'EMEP_IN.nc', ! for each of the IO IC/BC files, diff --git a/dependencies b/dependencies index 1b34d82..acbb292 100644 --- a/dependencies +++ b/dependencies @@ -56,7 +56,7 @@ ModelConstants_ml.o : ModelConstants_ml.f90 Io_Nums_ml.o PhysicalConstants_ml.o MosaicOutputs_ml.o : MosaicOutputs_ml.f90 Wesely_ml.o Units_ml.o TimeDate_ml.o SmallUtils_ml.o OwnDataTypes_ml.o ModelConstants_ml.o MetFields_ml.o LocalVariables_ml.o Landuse_ml.o LandDefs_ml.o Io_Progs_ml.o EcoSystem_ml.o DerivedFields_ml.o CM_ChemGroups_ml.o CM_ChemSpecs_ml.o CM_ChemSpecs_ml.o CheckStop_ml.o AOTnPOD_ml.o My_Aerosols_ml.o : My_Aerosols_ml.f90 Setup_1dfields_ml.o PhysicalConstants_ml.o ModelConstants_ml.o MARS_ml.o EQSAM_ml.o Chem_ml.o CM_ChemSpecs_ml.o CM_ChemSpecs_ml.o CM_ChemSpecs_ml.o CheckStop_ml.o My_Derived_ml.o : My_Derived_ml.f90 TimeDate_ml.o SmallUtils_ml.o Par_ml.o OwnDataTypes_ml.o MosaicOutputs_ml.o ModelConstants_ml.o MetFields_ml.o LandDefs_ml.o Io_Progs_ml.o GridValues_ml.o EmisDef_ml.o CM_ChemSpecs_ml.o CM_ChemGroups_ml.o CM_ChemSpecs_ml.o CM_ChemSpecs_ml.o CM_ChemSpecs_ml.o Chem_ml.o CheckStop_ml.o AOTnPOD_ml.o -My_SOA_ml.o : My_SOA_ml.f90 TimeDate_ml.o Setup_1dfields_ml.o PhysicalConstants_ml.o Par_ml.o ModelConstants_ml.o GridValues_ml.o CM_ChemGroups_ml.o CM_ChemSpecs_ml.o CM_ChemSpecs_ml.o Chem_ml.o Functions_ml.o +My_SOA_ml.o : My_SOA_ml.f90 CM_ChemSpecs_ml.o CM_ChemSpecs_ml.o Setup_1dfields_ml.o PhysicalConstants_ml.o ModelConstants_ml.o My_Outputs_ml.o : My_Outputs_ml.f90 Units_ml.o TimeDate_ml.o SmallUtils_ml.o Par_ml.o OwnDataTypes_ml.o ModelConstants_ml.o DerivedFields_ml.o CM_ChemGroups_ml.o CM_ChemSpecs_ml.o CM_ChemSpecs_ml.o CM_ChemSpecs_ml.o CM_ChemSpecs_ml.o CheckStop_ml.o NetCDF_ml.o : NetCDF_ml.f90 SmallUtils_ml.o Functions_ml.o TimeDate_ExtraUtil_ml.o TimeDate_ml.o PhysicalConstants_ml.o Par_ml.o OwnDataTypes_ml.o ModelConstants_ml.o InterpolationRoutines_ml.o GridValues_ml.o CM_ChemSpecs_ml.o CM_ChemSpecs_ml.o CM_ChemSpecs_ml.o CM_ChemSpecs_ml.o CheckStop_ml.o Chem_ml.o My_Outputs_ml.o Nest_ml.o : Nest_ml.f90 Units_ml.o TimeDate_ExtraUtil_ml.o TimeDate_ml.o Par_ml.o OwnDataTypes_ml.o NetCDF_ml.o MetFields_ml.o ModelConstants_ml.o Io_ml.o GridValues_ml.o Functions_ml.o CM_ChemSpecs_ml.o Chem_ml.o CM_ChemSpecs_ml.o CheckStop_ml.o ExternalBICs_ml.o diff --git a/global2local.f90 b/global2local.f90 old mode 100644 new mode 100755 diff --git a/local2global.f90 b/local2global.f90 old mode 100644 new mode 100755 diff --git a/modrun.sh b/modrun.sh index e8084be..d8acba7 100755 --- a/modrun.sh +++ b/modrun.sh @@ -1,17 +1,19 @@ #!/bin/bash + # Minimalistic script for run the Unified EMEP model # Link the input data +#inputdir=/home/metno/mifasv/OpenSource201310 inputdir=. ln -s $inputdir/met/* . # Driving meteorology ln -s $inputdir/input/* . # Other input files # Define some run parameters -trendyear=2010 # emission year +trendyear=2011 # emission year runlabel1=Base # short label runlabel2=Opensource_setup # long label -startdate="2010 01 01" # start date (metdata) - enddate="2010 01 01" # end date (metdata) +startdate="2011 01 01" # start date (metdata) + enddate="2011 01 01" # end date (metdata) # Put the run parameters in a temporary file cat > INPUT.PARA << EOF @@ -23,7 +25,7 @@ $enddate EOF # Run the model -mpirun $inputdir/code/Unimod +mpiexec $inputdir/code/Unimod # Clean the links to the input data and remove INPUT.PARA ls $inputdir/met |xargs rm