diff --git a/blueprints/highres_chemistry/VERSION b/blueprints/highres_chemistry/VERSION new file mode 100644 index 0000000000000000000000000000000000000000..8fdd954df9831dfd29ceec0d74829b02f3f5d8c3 --- /dev/null +++ b/blueprints/highres_chemistry/VERSION @@ -0,0 +1 @@ +22 \ No newline at end of file diff --git a/blueprints/highres_chemistry/WRF_DOMAIN_PLOT_COAWST.png b/blueprints/highres_chemistry/WRF_DOMAIN_PLOT_COAWST.png new file mode 100644 index 0000000000000000000000000000000000000000..9c845e78a3e46eb964a687bb796fb0db111bcfbf Binary files /dev/null and b/blueprints/highres_chemistry/WRF_DOMAIN_PLOT_COAWST.png differ diff --git a/blueprints/highres_chemistry/config.bash b/blueprints/highres_chemistry/config.bash new file mode 100644 index 0000000000000000000000000000000000000000..536e8c74c4a8afbe30067f2550c00ee9aa964fe4 --- /dev/null +++ b/blueprints/highres_chemistry/config.bash @@ -0,0 +1,99 @@ +#!/bin/bash +# ------------------------------------------------------------------------------ +# WRFOTRON v 2.0b +# Christoph Knote (LMU Munich, Germany) +# 06/2016 +# christoph.knote@lmu.de +# ------------------------------------------------------------------------------ + +# path to the WRFotron installation +chainDir=${HOME}/wrfotron + +# --- Executable locations --- + +# WPS installation directory +WPSDir=${WPS_SRC_PATH} +# WRF installation directory +WRFDir=${WRF_SRC_PATH} + +# --- Input data settings --- + +# path to geogrid input data +geogDir=/alcc/gpfs2/home/mbees/data/geog/ + +# meteo input +# Vtable for the chosen meteo input +metVtableFile=${WPS_SRC_PATH}/ungrib/Variable_Tables/Vtable.GFS +# time increment in hours +metInc=1 +# full path to a met input file - you can use any "%<>" abbreviations known +# to the "date" command +metFilePattern="/alcc/gpfs2/home/mbees/data/meteo/GFS_OPR/GF%Y%m%d%H" +# example: +# "/glade/p/rda/data/ds083.2/grib2/%Y/%Y.%m/fnl_%Y%m%d_%H_00.grib2" + +# --- Pre/Postprocessing settings --- + +# prepararation script +preScriptPath=NONEXISTENT.bash + +# postprocessing scripts (arbitrary) +postScriptPath=NONEXISTENT.bash +# postprocessing scripts (actions for each wrfout file) +postPerFileScriptPath=NONEXISTENT.bash + +# --- Working directories --- + +# where the WRF will be run (some fast, large disk like "scratch" or similar) +workDir=${SCRATCH}/WRF/work/ + +# where the unprocessed WRF output will be stored +stagingRootDir=${SCRATCH}/WRF/staging/ +# where the WRF output will be stored +archiveRootDir=${SCRATCH}/archive/WRF/ +# where the WRF restart files will be stored +restartRootDir=${SCRATCH}/WRF/restart/ + +# remove run directory after run is finished? +removeRunDir=false + +# --- Chemistry --- + +withChemistry=true + +# WRF-Chem installation directory +WRFChemDir=${WRF_CHEM_SRC_PATH} + +# path to utility executables +meganBioEmissBin=megan_bio_emiss +mozbcBin=mozbc +weselyBin=wesely +exo_coldensBin=exo_coldens +anthro_emisBin=anthro_emis +fire_emisBin=fire_emis + +# path to Wesely and Exo_Coldens input data +wesColdensDataPath=${WESELY_EXO_COLDENS_DATA} + +# path to MEGAN input data +MEGANdir=/alcc/gpfs2/home/mbees/data/emissions/biogenic/MEGAN + +# use anthro_emiss or predefined files? +emissUseAnthroEmiss=false +# raw emission input - the files you read in with anthro_emiss +emissDir=/alcc/gpfs2/home/mbees/data/emissions/anthropogenic/EDGARv5/MOZART_MOSAIC +# emission conversion script for anthro_emis - must match emissions in emissDir +emissInpFile=emis_edgarv5_mozmos.inp +# year the emissions are valid for (for offset calculation) +emissYear=2015 + +# FINN fires +fireFilePattern="/alcc/gpfs2/home/mbees/data/emissions/fires/FINN/GLOB_MOZ4_%Y%j.txt" +fireInpFile=finn_fires.inp + +# boundary condition input +chembcFilePattern="/alcc/gpfs2/home/mbees/data/chembc/WACCM/WACCM%Y%m%d" +chembcInpFile=waccm.inp + +# TUV photolysis option 4 data file +TUVDataPath="/alcc/gpfs2/home/mbees/data/tuv/TUV.phot.bz2" diff --git a/blueprints/highres_chemistry/emis_edgarv5_mozmos.inp b/blueprints/highres_chemistry/emis_edgarv5_mozmos.inp new file mode 100644 index 0000000000000000000000000000000000000000..5bd229eb59b94700be346daa82ac08be472d4ad3 --- /dev/null +++ b/blueprints/highres_chemistry/emis_edgarv5_mozmos.inp @@ -0,0 +1,29 @@ +&CONTROL + domains = __domains__ + anthro_dir = '__emissDir__' + src_file_prefix = 'edgar_v5_2015_' + src_file_suffix = '_0.1x0.1.nc' + src_names = 'CO(28)','NOx(30)','SO2(64)','NH3(17)','BC(12)','OC(12)','PM2.5(1)','PM10(1)','BENZENE(78)','BIGALK(72)','BIGENE(56)','C2H2(26)','C2H4(28)', + 'C2H5OH(46)','C2H6(30)','C3H6(42)','C3H8(44)','CH2O(30)', + 'CH3CHO(44)','CH3COCH3(58)','CH3OH(32)','CH3COOH(60)','HCOOH(46)', + 'MEK(72)','TOLUENE(92)','XYLENES(106)' + sub_categories = 'TOTAL', + serial_output = .true., + start_output_time = '__startDate__', + stop_output_time = '__startDate__', + output_interval = 3600, + data_yrs_offset = __emissYearOffset__, + emissions_zdim_stag = 1, + emis_map = 'CO->CO','NO->0.8*NOx','NO2->0.2*NOx','SO2->SO2','NH3->NH3','BIGALK->BIGALK','BIGENE->BIGENE', + 'C2H4->C2H4','C2H5OH->C2H5OH','C2H6->C2H6','CH2O->CH2O','CH3CHO->CH3CHO', + 'CH3COCH3->CH3COCH3','CH3OH->CH3OH','MEK->MEK','TOLUENE->TOLUENE', + 'C3H6->C3H6','C3H8->C3H8','ISOP->0.0*CO','C10H16->0.0*CO', + 'SULF->0.0*SO2','C2H2->0.00561790*CO','BENZENE->BENZENE','XYLENE->XYLENES', + 'GLY->0.0*CO','MACR->0.0*CO','MGLY->0.0*CO','MVK->0.0*CO', + 'HCOOH->0.0*CO','HONO->0.0*CO','VOCA->0.04*CO','VOCBB->0.0*CO', + 'ECI(a)->0.1*BC','ECJ(a)->0.9*BC','ORGI(a)->0.1*OC','ORGJ(a)->0.9*OC','PM25I(a)->0.1*PM2.5', + 'PM25J(a)->0.9*PM2.5','PM_10(a)->PM10 + -1.0*PM2.5','SO4I(a)->0.0*PM10','SO4J(a)->0.0*PM10','NO3I(a)->0.0*PM10', + 'NO3J(a)->0.0*PM10','NH4I(a)->0.0*PM10','NH4J(a)->0.0*PM10','NAI(a)->0.0*PM10','NAJ(a)->0.0*PM10', + 'CLI(a)->0.0*PM10','CLJ(a)->0.0*PM10','CO_A->CO','CO_BB->0.0*CO','ORGI_A(a)->0.0*PM10', + 'ORGI_BB(a)->0.0*PM10','ORGJ_A(a)->0.0*PM10','ORGJ_BB(a)->0.0*PM10' +/ diff --git a/blueprints/highres_chemistry/exo_coldens.inp b/blueprints/highres_chemistry/exo_coldens.inp new file mode 100644 index 0000000000000000000000000000000000000000..19febe530505bd8c7a3660f8b35bddbd85f1341b --- /dev/null +++ b/blueprints/highres_chemistry/exo_coldens.inp @@ -0,0 +1,6 @@ +&control + +wrf_dir = './' +domains = __domains__, + +/ diff --git a/blueprints/highres_chemistry/finn_fires.inp b/blueprints/highres_chemistry/finn_fires.inp new file mode 100644 index 0000000000000000000000000000000000000000..96ff076731d5f1a10e78f908a19614d3b5ad6b7d --- /dev/null +++ b/blueprints/highres_chemistry/finn_fires.inp @@ -0,0 +1,30 @@ +&control +fire_directory = './' +fire_filename(1) = 'FINN_FIRES.txt' +start_date = '__startYear__-__startMonth__-__startDay__', +end_date = '__endYear__-__endMonth__-__endDay__', +Model = 'WRF', +wrf_directory = './', +domains = __domain__, +defaultUnits = 'molecules/cm^2/s', +FinnVers = '1.5', + +wrf2fire_map = 'co -> CO', 'no -> NO', 'no2 -> NO2', 'so2 -> SO2', + 'nh3 -> NH3', 'open -> BIGALD', 'bigalk -> BIGALK', + 'bigene -> BIGENE', 'apin -> C10H16', 'c2h4 -> C2H4', + 'c2h5oh -> C2H5OH', 'c2h6 -> C2H6', 'c3h6 -> C3H6', + 'c3h8 -> C3H8', 'ch2o -> CH2O', 'ald -> CH3CHO', + 'ch3cn -> CH3CN', 'ACET -> CH3COCH3', 'mgly -> CH3COCHO', + 'ch3cooh -> CH3COOH', 'ch3oh -> CH3OH', 'cres -> CRESOL', + 'glyald -> GLYALD', 'hcn -> HCN', 'acetol -> HYAC', + 'isop -> ISOP', 'macr -> MACR', 'mek -> MEK', 'mvk -> MVK', + 'toluene -> TOLUENE', 'hcooh -> HCOOH', 'c2h2 -> C2H2', + 'oc_a01 -> 0.01*OC;aerosol', 'oc_a02 -> 0.09*OC;aerosol', + 'oc_a03 -> 0.70*OC;aerosol', 'oc_a04 -> 0.20*OC;aerosol', + 'bc_a01 -> 0.01*BC;aerosol', 'bc_a02 -> 0.09*BC;aerosol', + 'bc_a03 -> 0.70*BC;aerosol', 'bc_a04 -> 0.20*BC;aerosol', + 'oin_a04 -> 1.00*PM10+-1.00*PM25+-0.20*OC+-0.20*BC;aerosol', + 'oin_a01 -> 0.01*PM25+-0.01*OC+-0.01*BC;aerosol', + 'oin_a02 -> 0.09*PM25+-0.09*OC+-0.09*BC;aerosol', + 'oin_a03 -> 0.90*PM25+-0.90*OC+-0.90*BC;aerosol' +/ diff --git a/blueprints/highres_chemistry/iofields.chem b/blueprints/highres_chemistry/iofields.chem new file mode 100644 index 0000000000000000000000000000000000000000..a4590190f58881c0d348293ecd02cc79f7f73bff --- /dev/null +++ b/blueprints/highres_chemistry/iofields.chem @@ -0,0 +1,20 @@ ++:h:0:alt ++:h:0:gsw ++:h:0:extaer1 ++:h:0:extaer2 ++:h:0:extaer3 ++:h:0:extaer4 ++:h:0:tauaer1 ++:h:0:tauaer2 ++:h:0:tauaer3 ++:h:0:tauaer4 ++:h:0:gaer1 ++:h:0:gaer2 ++:h:0:gaer3 ++:h:0:gaer4 ++:h:0:waer1 ++:h:0:waer2 ++:h:0:waer3 ++:h:0:waer4 ++:h:0:brnch_rto ++:h:0:ccn1 diff --git a/blueprints/highres_chemistry/iofields.met b/blueprints/highres_chemistry/iofields.met new file mode 100644 index 0000000000000000000000000000000000000000..10455b4f9351e774714eb0564e92e27393fc1a06 --- /dev/null +++ b/blueprints/highres_chemistry/iofields.met @@ -0,0 +1,4 @@ ++:h:0:alt ++:h:0:gsw ++:h:0:tke ++:h:0:SWDDIF \ No newline at end of file diff --git a/blueprints/highres_chemistry/megan_bio_emiss.inp b/blueprints/highres_chemistry/megan_bio_emiss.inp new file mode 100644 index 0000000000000000000000000000000000000000..e54891ebcb643132915bc5eee3fc4d3fef056dea --- /dev/null +++ b/blueprints/highres_chemistry/megan_bio_emiss.inp @@ -0,0 +1,9 @@ +&control + +domains = __domains__, +start_lai_mnth = 1, +end_lai_mnth = 12, +wrf_dir = '.', +megan_dir = '__MEGANdir__' + +/ diff --git a/blueprints/highres_chemistry/mozart4.inp b/blueprints/highres_chemistry/mozart4.inp new file mode 100644 index 0000000000000000000000000000000000000000..2d33d08d2a12e0a4fefc85b253389310938d989e --- /dev/null +++ b/blueprints/highres_chemistry/mozart4.inp @@ -0,0 +1,86 @@ +&control + +do_ic = __do_ic__, +do_bc = __do_bc__, + +domain = __domain__, + +dir_wrf = './', +dir_moz = './', +fn_moz = 'moz0000.nc', +moz_var_suffix='', + +spc_map = +'o3 -> O3', +'no -> NO', +'no2 -> NO2', +'hno3 -> HNO3', +'ho -> OH', +'ho2 -> HO2', +'h2o2 -> H2O2', +'so2 -> SO2', +'co -> CO', +'hcho -> CH2O', +'c2h4 -> C2H4', +'c2h6 -> C2H6', +'c3h6 -> C3H6', +'c3h8 -> C3H8', +'ald -> CH3CHO', +'acet -> CH3COCH3', +'ch3cooh -> CH3COOH', +'ch3oh -> CH3OH', +'ch3ooh -> CH3OOH', +'pan -> PAN', +'macr -> MACR', +'mvk -> MVK', +'c2h5oh -> C2H5OH', +!'etooh -> C2H5OOH', +'bigene -> BIGENE', +'bigalk -> BIGALK', +'tol -> TOLUENE', +'benzene -> BENZENE', +'xylenes -> XYLENES', +! http://www.atmos-chem-phys.net/12/7215/2012/acp-12-7215-2012.pdf Fig 6 +'apin -> 0.75*MTERP', +'bpin -> 0.25*MTERP', +'isopr -> ISOP', +! +! aerosols +! +'oc_a01->0.1216*pom_a1+0.9886*soa1_a2+0.1216*soa1_a1+0.9886*soa2_a2+ 0.1216*soa2_a1+0.9886*soa3_a2+0.1216*soa3_a1+0.9886*soa4_a2+0.1216*soa4_a1+ 0.9886*soa5_a2+ 0.1216*soa5_a1;1.e9', +'oc_a02->0.7618*pom_a1+0.0114*soa1_a2+0.7618*soa1_a1+0.0114*soa2_a2+ 0.7618*soa2_a1+0.0114*soa3_a2+0.7618*soa3_a1+0.0114*soa4_a2+0.7618*soa4_a1+ 0.0114*soa5_a2+ 0.7618*soa5_a1;1.e9', +'oc_a03->0.1164*pom_a1+0.0000*soa1_a2+0.1164*soa1_a1+0.0000*soa2_a2+ 0.1164*soa2_a1+0.0000*soa3_a2+0.1164*soa3_a1+0.0000*soa4_a2+0.1164*soa4_a1+ 0.0000*soa5_a2+ 0.1164*soa5_a1;1.e9', +'oc_a04->0.0002*pom_a1+0.0000*soa1_a2+0.0002*soa1_a1+0.0000*soa2_a2+ 0.0002*soa2_a1+0.0000*soa3_a2+0.0002*soa3_a1+0.0000*soa4_a2+0.0002*soa4_a1+ 0.0000*soa5_a2+ 0.0002*soa5_a1;1.e9', +'bc_a01->0.1216*bc_a1+0.1216*bc_a4;1.e9', +'bc_a02->0.7618*bc_a1+0.7618*bc_a4;1.e9', +'bc_a03->0.1164*bc_a1+0.1164*bc_a4;1.e9', +'bc_a04->0.0002*bc_a1+0.0002*bc_a4;1.e9', +'so4_a01->0.9886*so4_a2+0.1216*so4_a1+0.0000*so4_a3;1.e9', +'so4_a02->0.0114*so4_a2+0.7618*so4_a1+0.0000*so4_a3;1.e9', +'so4_a03->0.0000*so4_a2+0.1164*so4_a1+0.0995*so4_a3;1.e9', +'so4_a04->0.0000*so4_a2+0.0002*so4_a1+0.9003*so4_a3;1.e9', +'nh4_a01->0.1856*so4_a2+0.0050*so4_a1+0.0000*so4_a3;1.e9', +'nh4_a02->0.0021*so4_a2+0.0930*so4_a1+0.0000*so4_a3;1.e9', +'nh4_a03->0.0000*so4_a2+0.0203*so4_a1+0.0186*so4_a3;1.e9', +'nh4_a04->0.0000*so4_a2+0.0000*so4_a1+0.1690*so4_a3;1.e9', +'no3_a01->0.0000*so4_a2+0.0000*so4_a1+0.0000*so4_a3;1.e9', +'no3_a02->0.0000*so4_a2+0.0000*so4_a1+0.0000*so4_a3;1.e9', +'no3_a03->0.0000*so4_a2+0.0000*so4_a1+0.0000*so4_a3;1.e9', +'no3_a04->0.0000*so4_a2+0.0000*so4_a1+0.0000*so4_a3;1.e9', +'na_a01->0.3889*ncl_a2+0.0479*ncl_a1+0.0000*ncl_a3;1.e9', +'na_a02->0.0045*ncl_a2+0.2997*ncl_a1+0.0000*ncl_a3;1.e9', +'na_a03->0.0000*ncl_a2+0.0458*ncl_a1+0.0391*ncl_a3;1.e9', +'na_a04->0.0000*ncl_a2+0.0000*ncl_a1+0.3542*ncl_a3;1.e9', +'cl_a01->0.5996*ncl_a2+0.0737*ncl_a1+0.0000*ncl_a3;1.e9', +'cl_a02->0.0068*ncl_a2+0.4621*ncl_a1+0.0000*ncl_a3;1.e9', +'cl_a03->0.0000*ncl_a2+0.0709*ncl_a1+0.0604*ncl_a3;1.e9', +'cl_a04->0.0000*ncl_a2+0.0001*ncl_a1+0.5462*ncl_a3;1.e9', +'oin_a01->0.9886*dst_a2+0.1216*dst_a1+0.0000*dst_a3;1.e9', +'oin_a02->0.0114*dst_a2+0.7618*dst_a1+0.0002*dst_a3;1.e9', +'oin_a03->0.0000*dst_a2+0.1164*dst_a1+0.0995*dst_a3;1.e9', +'oin_a04->0.0000*dst_a2+0.0002*dst_a1+0.9003*dst_a3;1.e9', +'num_a01->0.9996*num_a2+0.7135*num_a1+0.0000*num_a3;1.0', +'num_a02->0.0004*num_a2+0.2470*num_a1+0.2118*num_a3;1.0', +'num_a03->0.0000*num_a2+0.0016*num_a1+0.6258*num_a3;1.0', +'num_a04->0.0000*num_a2+0.0000*num_a1+0.3501*num_a3;1.0', +/ diff --git a/blueprints/highres_chemistry/mozbc.inp b/blueprints/highres_chemistry/mozbc.inp new file mode 100644 index 0000000000000000000000000000000000000000..97097238df53fcde7ceb9676d3057782dcb00271 --- /dev/null +++ b/blueprints/highres_chemistry/mozbc.inp @@ -0,0 +1,60 @@ +&control + +do_ic = __do_ic__, +do_bc = __do_bc__, + +domain = __domain__, + +dir_wrf = './', +dir_moz = './', +fn_moz = 'moz0000.nc', +moz_var_suffix='', + +def_missing_var = .true. +spc_map = 'o3 -> O3', 'n2o -> N2O', 'no -> NO', +'no2 -> NO2', 'nh3 -> NH3', 'hno3 -> HNO3', 'hno4 -> HO2NO2', 'n2o5 -> N2O5', 'h2o2 -> H2O2', +'ch4 -> CH4', 'co -> CO', 'ch3ooh -> CH3OOH', +'hcho -> CH2O', 'ch3oh -> CH3OH', 'c2h4 -> C2H4', +'ald -> CH3CHO', 'acet -> CH3COCH3', 'mgly -> CH3COCHO', +'pan -> PAN', 'mpan -> MPAN', 'macr -> MACR', +'mvk -> MVK', 'c2h6 -> C2H6', 'c3h6 -> C3H6', 'c3h8 -> C3H8', 'c2h5oh -> C2H5OH', 'c10h16 -> MTERP', +'isopr -> ISOP','acetol -> HYAC', 'mek -> MEK', +'bigene -> BIGENE', 'bigalk -> BIGALK', +'tol -> TOLUENE+BENZENE+XYLENES', 'cres -> CRESOL', 'dms -> DMS', 'so2 -> SO2', +'oc_a01->0.12*pom_a1+0.99*soa1_a2+0.12*soa1_a1+0.99*soa2_a2+0.12*soa2_a1+0.99*soa3_a2+0.12*soa3_a1+0.99*soa4_a2+0.12*soa4_a1+0.99*soa5_a2+0.12*soa5_a1;1.e9', +'oc_a02->0.76*pom_a1+0.01*soa1_a2+0.76*soa1_a1+0.01*soa2_a2+0.76*soa2_a1+0.01*soa3_a2+0.76*soa3_a1+0.01*soa4_a2+0.76*soa4_a1+0.01*soa5_a2+0.76*soa5_a1;1.e9', +'oc_a03->0.12*pom_a1+0.00*soa1_a2+0.12*soa1_a1+0.00*soa2_a2+0.12*soa2_a1+0.00*soa3_a2+0.12*soa3_a1+0.00*soa4_a2+0.12*soa4_a1+0.00*soa5_a2+0.12*soa5_a1;1.e9', +'oc_a04->0.00*pom_a1+0.00*soa1_a2+0.00*soa1_a1+0.00*soa2_a2+0.00*soa2_a1+0.00*soa3_a2+0.00*soa3_a1+0.00*soa4_a2+0.00*soa4_a1+0.00*soa5_a2+0.00*soa5_a1;1.e9', +'bc_a01->0.1216*bc_a1+0.1216*bc_a4;1.e9', +'bc_a02->0.7618*bc_a1+0.7618*bc_a4;1.e9', +'bc_a03->0.1164*bc_a1+0.1164*bc_a4;1.e9', +'bc_a04->0.0002*bc_a1+0.0002*bc_a4;1.e9', +'so4_a01->0.9886*so4_a2+0.1216*so4_a1+0.0000*so4_a3;1.e9', +'so4_a02->0.0114*so4_a2+0.7618*so4_a1+0.0002*so4_a3;1.e9', +'so4_a03->0.0000*so4_a2+0.1216*so4_a1+0.0995*so4_a3;1.e9', +'so4_a04->0.0000*so4_a2+0.0002*so4_a1+0.9003*so4_a3;1.e9', +'nh4_a01->0.1856*so4_a2+0.0050*so4_a1+0.0000*so4_a3;1.e9', +'nh4_a02->0.0021*so4_a2+0.0930*so4_a1+0.0000*so4_a3;1.e9', +'nh4_a03->0.0000*so4_a2+0.0203*so4_a1+0.0186*so4_a3;1.e9', +'nh4_a04->0.0000*so4_a2+0.0000*so4_a1+0.1690*so4_a3;1.e9', +'no3_a01->0.0000*so4_a2+0.0000*so4_a1+0.0000*so4_a3;1.e9', +'no3_a02->0.0000*so4_a2+0.0000*so4_a1+0.0000*so4_a3;1.e9', +'no3_a03->0.0000*so4_a2+0.0000*so4_a1+0.0000*so4_a3;1.e9', +'no3_a04->0.0000*so4_a2+0.0000*so4_a1+0.0000*so4_a3;1.e9', +'na_a01->0.3889*ncl_a2+0.0479*ncl_a1+0.0000*ncl_a3;1.e9', +'na_a02->0.0045*ncl_a2+0.2997*ncl_a1+0.0000*ncl_a3;1.e9', +'na_a03->0.0000*ncl_a2+0.0458*ncl_a1+0.0391*ncl_a3;1.e9', +'na_a04->0.0000*ncl_a2+0.0000*ncl_a1+0.3542*ncl_a3;1.e9', +'cl_a01->0.5995*ncl_a2+0.0737*ncl_a1+0.0000*ncl_a3;1.e9', +'cl_a02->0.0068*ncl_a2+0.4621*ncl_a1+0.0000*ncl_a3;1.e9', +'cl_a03->0.0000*ncl_a2+0.0709*ncl_a1+0.0604*ncl_a3;1.e9', +'cl_a04->0.0000*ncl_a2+0.0001*ncl_a1+0.5462*ncl_a3;1.e9', +'oin_a01->0.9886*dst_a2+0.1216*dst_a1+0.0000*dst_a3;1.e9', +'oin_a02->0.0114*dst_a2+0.7618*dst_a1+0.0002*dst_a3;1.e9', +'oin_a03->0.0000*dst_a2+0.1216*dst_a1+0.0995*dst_a3;1.e9', +'oin_a04->0.0000*dst_a2+0.0002*dst_a1+0.9003*dst_a3;1.e9', +'num_a01->0.9996*num_a2+0.7135*num_a1+0.0000*num_a3;1.0', +'num_a02->0.0004*num_a2+0.2847*num_a1+0.0239*num_a3;1.0', +'num_a03->0.0000*num_a2+0.0016*num_a1+0.6258*num_a3;1.0', +'num_a04->0.0000*num_a2+0.0000*num_a1+0.3501*num_a3;1.0', +/ diff --git a/blueprints/highres_chemistry/namelist.chem b/blueprints/highres_chemistry/namelist.chem new file mode 100644 index 0000000000000000000000000000000000000000..4ce3023c5d852fab468a4b7548be2c3cd1eec6dc --- /dev/null +++ b/blueprints/highres_chemistry/namelist.chem @@ -0,0 +1,44 @@ + &chem + kemit = 1, + chem_opt = 202, 202, 202, + bioemdt = 12., 4., 4., + photdt = 12., 4., 4., + chemdt = 12., 4., 4., + io_style_emissions = 2, + emiss_inpt_opt = 102, 102, 102, + emiss_opt = 10, 10, 10, + chem_in_opt = 1, 2, 2, + phot_opt = 4, 4, 4, + gas_drydep_opt = 1, 1, 1, + aer_drydep_opt = 1, 1, 1, + bio_emiss_opt = 3, 3, 3, + gas_bc_opt = 1, 1, 1, + gas_ic_opt = 1, 1, 1, + aer_bc_opt = 1, 1, 1, + aer_ic_opt = 1, 1, 1, + gaschem_onoff = 1, 1, 1, + aerchem_onoff = 1, 1, 1, + wetscav_onoff = 1, 1, 1, + cldchem_onoff = 0, 0, 0, + vertmix_onoff = 1, 1, 1, + chem_conv_tr = 1, 1, 1, + conv_tr_wetscav = 1, 1, 1, + conv_tr_aqchem = 1, 1, 1, + seas_opt = 2, + dust_opt = 3, + dmsemis_opt = 1, + aer_ra_feedback = 1, 1, 1, + aer_op_opt = 1, 1, 1, + ne_area = 500, + opt_pars_out = 1, + have_bcs_chem = .true., + have_bcs_upper = .false., + biomass_burn_opt = 5, 5, 5, + plumerisefire_frq = 30, 30, 30, + scale_fire_emiss = .true., .true., .true., + n2o5_hetchem = 1, + lnox_opt = 1, 1, 1, + diagnostic_chem = 0, 0, 0, + diagnostic_dep = 0, 0, 0, + / + diff --git a/blueprints/highres_chemistry/namelist.wps b/blueprints/highres_chemistry/namelist.wps new file mode 100644 index 0000000000000000000000000000000000000000..7499aec8602462364db34a2c9f4c1e6c2ddb7990 --- /dev/null +++ b/blueprints/highres_chemistry/namelist.wps @@ -0,0 +1,39 @@ +&share + wrf_core = 'ARW', + max_dom = 3, + start_date = '__startDate__','__startDate__','__startDate__' + end_date = '__endDate__','__startDate__','__startDate__', + interval_seconds = __metIncSec__, + io_form_geogrid = 2, +/ + +&geogrid + parent_id = 1, 1, 2, + parent_grid_ratio = 1, 9, 5, + i_parent_start = 1, 96, 69, + j_parent_start = 1, 75, 74, + e_we = 185, 172, 151, + e_sn = 186, 163, 156, +! geog_data_res = 'modis_fpar+modis_lai+modis_lakes+modis_30s+modis_15s+30s', 'modis_fpar+modis_lai+modis_lakes+modis_30s+modis_15s+30s', +! 'modis_fpar+modis_lai+modis_lakes+modis_30s+modis_15s+30s', + dx = 20000, + dy = 20000, + map_proj = 'lambert', + ref_lat = 50.0, + ref_lon = 7.5, + truelat1 = 30.0, + truelat2 = 60.0, + stand_lon = 7.5, + geog_data_path = '__geogDir__' +/ + +&ungrib + out_format = 'WPS', + prefix = 'FILE', +/ + +&metgrid + fg_name = 'FILE', + constants_name = 'TAVGSFC', + io_form_metgrid = 2, +/ diff --git a/blueprints/highres_chemistry/namelist.wrf b/blueprints/highres_chemistry/namelist.wrf new file mode 100644 index 0000000000000000000000000000000000000000..e1072e9e4a1059c0a74db538360df9e8378be088 --- /dev/null +++ b/blueprints/highres_chemistry/namelist.wrf @@ -0,0 +1,160 @@ + &time_control + start_year = __startYear__, __startYear__, __startYear__, + start_month = __startMonth__, __startMonth__, __startMonth__, + start_day = __startDay__, __startDay__, __startDay__, + start_hour = __startHour__, __startHour__, __startHour__, + start_minute = 00, 00, 00, + start_second = 00, 00, 00, + end_year = __endYear__, __endYear__, __endYear__, + end_month = __endMonth__, __endMonth__, __endMonth__, + end_day = __endDay__, __endDay__, __endDay__, + end_hour = __endHour__, __endHour__, __endHour__, + end_minute = 00, 00, 00, + end_second = 00, 00, 00, + interval_seconds = __metIncSec__, + input_from_file = .true., .true., .true., + history_interval = 60, 60, 60, + frames_per_outfile = 1, 1, 1, + restart = __isRestart__, + restart_interval = __restartInterval__, + override_restart_timers = .true., + io_form_history = 2 + io_form_restart = 2 + io_form_input = 2 + io_form_boundary = 2 + debug_level = 0 + iofields_filename = "iofields" ,"iofields", "iofields", + force_use_old_data = .true., + ! CHEM + / + + &domains + use_adaptive_time_step = .false., + step_to_output_time = .true., + target_cfl = 1.0, 1.0, 1.0, + max_step_increase_pct = 5, 51, 51, + starting_time_step = -1, -1, -1, + max_time_step = 200, 40, 3, + min_time_step = -1, -1, -1, + adaptation_domain = 3, + time_step = 120, + time_step_fract_num = 0, + time_step_fract_den = 1, + max_dom = 3, + e_we = 185, 172, 151, + e_sn = 186, 163, 156, + e_vert = 33, 33, 33, + num_metgrid_levels = 34, + num_metgrid_soil_levels = 4, + dx = 20000, 2222, 400, + dy = 20000, 2222, 400, + grid_id = 1, 2, 3, + parent_id = 0, 1, 2, + i_parent_start = 1, 96, 69, + j_parent_start = 1, 75, 74, + parent_grid_ratio = 1, 9, 5, + parent_time_step_ratio = 1, 9, 5, + feedback = 0, + / + + &physics + mp_physics = 10, 10, 10, + progn = 1, 1, 1, + prec_acc_dt = 60., 60., 60., + ra_lw_physics = 4, 4, 4, + ra_sw_physics = 4, 4, 4, + radt = 16, 4, 4, + aer_opt = 1, + sf_sfclay_physics = 5, 5, 5, ! Ravans suggestion + sf_surface_physics = 2, 2, 2, ! NOAH + sf_surface_mosaic = 0, + sf_lake_physics = 1, 1, 1, + sf_urban_physics = 1, 1, 1, + bl_pbl_physics = 5, 5, 5, ! MYNN 2.5 + bl_mynn_cloudpdf = 1, + bl_mynn_cloudmix = 1, + scalar_pblmix = 1, + tracer_pblmix = 1, +! grav_settling = 0, 2, 2, ! settling of fog / cloud droplets + bldt = 0, 0, 0, + cu_physics = 5, 5, 0, ! Grell 3D + cudt = 1, + cugd_avedx = 1, ! set to 3 for 4km run, 1 for 36km + ishallow = 0, + iz0tlnd = 1, ! thermal roughness length over land determined by vegetation type + isfflx = 1, + ifsnow = 1, + icloud = 1, + icloud_bl = 1, + num_soil_layers = 4, + cu_rad_feedback = .true., .true.,.false., + cu_diag = 1, 1, 0, + slope_rad = 0, 0, 1, + topo_shading = 0, 0, 0, + num_land_cat = 21, + / + + &fdda + grid_fdda = 1, 0, 0, + gfdda_inname = "wrffdda_d", + gfdda_end_h = 312, 0, 0, + gfdda_interval_m = 60, 0, 0, + if_no_pbl_nudging_uv = 1, 1, 1, + if_no_pbl_nudging_t = 1, 1, 1, + if_no_pbl_nudging_q = 1, 1, 1, + if_zfac_uv = 0, 0, 0, + k_zfac_uv = 2, + if_zfac_t = 0, 0, 0, + k_zfac_t = 2, + if_zfac_q = 0, 0, 0, + k_zfac_q = 2, + guv = 0.0006, 0.0006, 0.0006, + gt = 0.0006, 0.0006, 0.0006, + gq = 0.0006, 0.0006, 0.0006, + if_ramping = 0, + dtramp_min = 360, + io_form_gfdda = 2, + / + + &dynamics + rk_ord = 3, + time_step_sound = 4, 4, 4, + w_damping = 1, + diff_opt = 1, 1, 2, + km_opt = 4, 4, 4, + diff_6th_opt = 0, 0, 0, + diff_6th_factor = 0.12, 0.12, 0.75, + epssm = 0.1, 0.1, 0.5, + sfs_opt = 0, 0, 0, + base_temp = 290. + damp_opt = 3, + zdamp = 5000., 5000., 5000., + dampcoef = 0.15, 0.15, 0.15, + khdif = 0, 0, 0, + kvdif = 0, 0, 0, + mix_isotropic = 0, 0, 1, + non_hydrostatic = .true., + moist_adv_opt = 2, 2, 4, + scalar_adv_opt = 2, 2, 4, + tke_adv_opt = 2, 2, 4, + do_avgflx_em = 1, 1, 1, + hybrid_opt = 0, + use_theta_m = 1, + / + + &bdy_control + spec_bdy_width = 5, + spec_zone = 1, + relax_zone = 4, + specified = .true., .false., .false., + nested = .false., .true., .true., + / + + &grib2 + / + + &namelist_quilt + nio_tasks_per_group = 0, + nio_groups = 1 + / + diff --git a/blueprints/highres_chemistry/namelist_timecontrol_chem_patch.wrf b/blueprints/highres_chemistry/namelist_timecontrol_chem_patch.wrf new file mode 100644 index 0000000000000000000000000000000000000000..70e38516aa7c7de0ad3140ae6ee7a563323fef37 --- /dev/null +++ b/blueprints/highres_chemistry/namelist_timecontrol_chem_patch.wrf @@ -0,0 +1,12 @@ + io_form_auxinput5 = 2, + auxinput5_inname = 'wrfchemi_d_', + auxinput5_interval_h = 1, 1, + frames_per_auxinput5 = 1, 1, + io_form_auxinput6 = 2, + auxinput6_inname = 'wrfbiochemi_d' + auxinput6_interval_h = 2400, 2400, + io_form_auxinput7 = 2, + auxinput7_inname = 'wrffirechemi_d_', + auxinput7_interval_m = 60, 60, + frames_per_auxinput7 = 1, 1, + \ No newline at end of file diff --git a/blueprints/highres_chemistry/pp/plotting/crontab b/blueprints/highres_chemistry/pp/plotting/crontab new file mode 100644 index 0000000000000000000000000000000000000000..e6d46be2184b9d4a1599d4d01e6b96682adb0083 --- /dev/null +++ b/blueprints/highres_chemistry/pp/plotting/crontab @@ -0,0 +1,34 @@ +# Edit this file to introduce tasks to be run by cron. +# +# Each task to run has to be defined through a single line +# indicating with different fields when the task will be run +# and what command to run for the task +# +# To define the time you can provide concrete values for +# minute (m), hour (h), day of month (dom), month (mon), +# and day of week (dow) or use '*' in these fields (for 'any').# +# Notice that tasks will be started based on the cron's system +# daemon's notion of time and timezones. +# +# Output of the crontab jobs (including errors) is sent through +# email to the user the crontab file belongs to (unless redirected). +# +# For example, you can run a backup of all your user accounts +# at 5 a.m every week with: +# 0 5 * * 1 tar -zcf /var/backups/home.tgz /home/ +# +# For more information see the manual pages of crontab(5) and cron(8) +# +# m h dom mon dow command +HOME=/home/m/met-wrf-chem +MAILTO=christoph.knote@physik.uni-muenchen.de + +# wrf meteogram +52 * * * * /usr/bin/sbatch /alcc/gpfs2/home/u/knotechr/wrfotron/blueprints/operational_chemistry/pp/plotting/meteograms/run.sh + +# wrf maps +00 07 * * * /usr/bin/sbatch /alcc/gpfs2/home/u/knotechr/wrfotron/blueprints/operational_chemistry/pp/plotting/maps/run.sh + +# wrf geotiffs +17 6,12,18 * * * /usr/bin/sbatch /alcc/gpfs2/home/u/knotechr/wrfotron/blueprints/operational_chemistry/pp/plotting/geotiffs/run.sh + diff --git a/blueprints/highres_chemistry/pp/plotting/geotiffs/run.sh b/blueprints/highres_chemistry/pp/plotting/geotiffs/run.sh new file mode 100755 index 0000000000000000000000000000000000000000..dbffc99908251525e4069bc210799e6757712032 --- /dev/null +++ b/blueprints/highres_chemistry/pp/plotting/geotiffs/run.sh @@ -0,0 +1,37 @@ +#!/bin/bash -l +#SBATCH --partition=alcc1 +#SBATCH -o /alcc/gpfs2/scratch/mbees/knotechr/WRF/postprocessing/geotiffs.%j.%N.out +#SBATCH -J geotiffs +#SBATCH --ntasks=1 +#SBATCH --mem=5G +#SBATCH --mail-type=FAIL +#SBATCH --mail-user=christoph.knote@med.uni-augsburg.de +#SBATCH --time=00:30:00 + +module load gnu gdal geos proj + +scriptPath=/alcc/gpfs2/home/u/knotechr/wrfotron/tools/nc_2_geotiff.py +wrfDataPathPattern="/alcc/gpfs2/scratch/mbees/knotechr/archive/WRF/operational_chemistry/wrfout___domain___%Y-%m-%d_%H:%M:%S" +firstDateToProcess=`date --date="2 days ago" "+%Y-%m-%d"` +intervalHours=1 +rangeHours=120 +plotDirPattern="/alcc/gpfs2/scratch/mbees/knotechr/plots/WRF/operational_chemistry/__domain__/geotiffs/%Y/%m/%d" + +for domain in d01 d02 +do + plotPathPattern="${plotDirPattern/__domain__/${domain}}/{species:s}_%Y-%m-%d_%H-%M.tif" + + for anHour in $(seq 0 ${intervalHours} ${rangeHours}) + do + aWRFDataPath=$(date -u --date="${firstDateToProcess} 0 UTC +${anHour} hours" "+${wrfDataPathPattern/__domain__/${domain}}") + + aPlotDir=$(date -u --date="${firstDateToProcess} 0 UTC +${anHour} hours" "+${plotDirPattern/__domain__/${domain}}") + aPlotPath=$(date -u --date="${firstDateToProcess} 0 UTC +${anHour} hours" "+${plotPathPattern}") + + if [ -f ${aWRFDataPath} ] + then + mkdir -p $aPlotDir + python3 $scriptPath -s "o3,O3" -s "no2,NO2" -s "co,CO" -s "PM2_5_DRY,PM2.5" -s "PM10,PM10" "${aWRFDataPath}" "${aPlotPath}" + fi + done +done diff --git a/blueprints/highres_chemistry/pp/plotting/maps/plots.py b/blueprints/highres_chemistry/pp/plotting/maps/plots.py new file mode 100644 index 0000000000000000000000000000000000000000..b44e3a1b5685f90abbf1644ac228271eeb429cea --- /dev/null +++ b/blueprints/highres_chemistry/pp/plotting/maps/plots.py @@ -0,0 +1,197 @@ + +from wrf import to_np, getvar, smooth2d, latlon_coords, interplevel, get_cartopy +import numpy as np +import matplotlib +from matplotlib import pyplot as plt +import cartopy.feature as cfeature +import cartopy.crs as ccrs +from matplotlib.colors import LinearSegmentedColormap, ListedColormap + + +from tools import total_prec, cartopy_lim_buffered, make_fig_and_ax, add_title_and_description + +def slp_map(t, nc, prec_prev, output_path, plevels=np.arange(950, 1050, 5)): +# + slp = getvar(nc, "slp") + mid_clouds = getvar(nc, "mid_cloudfrac") + prec = total_prec(nc, prec_prev) +# + # Smooth the sea level pressure since it tends to be noisy near the + # mountains + smooth_slp = smooth2d(slp, 5, cenweight=4) + # Get the latitude and longitude points + lats, lons = latlon_coords(slp) +# + # Get the WRF projection (CRS) object + crs = get_cartopy(slp) +# + # Create a figure + fig, ax, ax_cbarv, ax_cbarh = make_fig_and_ax(plt, crs, add_horiz_cbar=True) +# + # Draw the contours and filled contours + cmap = LinearSegmentedColormap.from_list("", ["white", "silver"]) + cmap2 = matplotlib.cm.viridis_r +# + con = ax.contourf(to_np(lons), to_np(lats), to_np(mid_clouds), + levels=[1/10, 1/4, 1/2, 3/4, 1], cmap=cmap, zorder=2, + transform=ccrs.PlateCarree()) + con2 = ax.contourf(to_np(lons), to_np(lats), to_np(prec), + levels=[0.5, 1, 2, 5, 10, 25, 50], cmap=cmap2, zorder=2, extend='max', + transform=ccrs.PlateCarree()) +# + linewidths = np.ones(plevels.shape) + linewidths[np.where(plevels == 1015)[0]] = 3 + cs = ax.contour(to_np(lons), to_np(lats), to_np(smooth_slp), levels=plevels, + linewidths=linewidths, colors="black", zorder=4, + transform=ccrs.PlateCarree()) + plt.clabel(cs, fmt='%1.0f', fontsize=8) +# + ax.coastlines(linewidth=0.3, zorder=2) + countries = cfeature.NaturalEarthFeature(category='cultural', name='admin_0_countries', scale='50m') + ax.add_feature(countries, facecolor='none', edgecolor='black', linewidth=0.3, zorder=3) +# + cbar = plt.colorbar(con, cax=ax_cbarv, ticks=[1/10, 1/4, 1/2, 3/4, 1]) + cbar2 = plt.colorbar(con2, cax=ax_cbarh, orientation="horizontal") +# + try: + cbar.ax.set_yticklabels(["1/10", "1/4", "1/2", "3/4", "1"], fontsize=10) + except: + print("Plotting tick labels failed") +# + ax.set_xlim(cartopy_lim_buffered(slp, "x")) + ax.set_ylim(cartopy_lim_buffered(slp, "y")) + + add_title_and_description(nc, "Sea Level Pressure [hPa], Mid-Level Cloud Cover, Total prec. in the last hour [mm] \n" + "WRF {init:s} +{fcsthour:d}", ax) +# + plt.savefig(output_path, bbox_inches='tight', pad_inches=0.1) + plt.close() + + return prec + + +def plevel_map(t, nc, output_path, plevel, phi_contour_intervals=np.arange(400, 800, 8)): +# + def interp_var(nc, variablename, p, plevel): + var = getvar(nc, variablename) + return interplevel(var, p, plevel) +# + p = getvar(nc, "p", units='hPa') +# + geop = interp_var(nc, "geopt", p, plevel) + geop = smooth2d(geop, 5, cenweight=4) +# + temp = interp_var(nc, "tc", p, plevel) +# + u, v = interp_var(nc, "uvmet", p, plevel) +# + # Get the latitude and longitude points + lats, lons = latlon_coords(temp) +# + # Get the WRF projection (CRS) object + crs = get_cartopy(temp) +# + # Create a figure + fig, ax, ax_cbarv, ax_cbarh = make_fig_and_ax(plt, crs) +# + # Temperature + lab = np.array([-46, -41, -36, -31, -26, -21, -16, -11, -6, -1, 4]) + con = ax.contourf(to_np(lons), to_np(lats), temp, levels=lab, + transform=ccrs.PlateCarree(), + cmap='viridis', zorder=0, extend='both') +# + # Geopotential height + linewidths = np.ones(phi_contour_intervals.shape) + linewidths[np.where(phi_contour_intervals == 552)[0]] = 3 + cs = ax.contour(to_np(lons), to_np(lats), geop/98.0, levels=phi_contour_intervals, + linewidths=linewidths, colors="black", + zorder=1, + transform=ccrs.PlateCarree()) + plt.clabel(cs, fmt='%1.0f', fontsize=10) +# +# + ax.coastlines(linewidth=0.3, zorder=2) + countries = cfeature.NaturalEarthFeature(category='cultural', name='admin_0_countries', scale='50m') + ax.add_feature(countries, facecolor='none', edgecolor='black', linewidth=0.3, zorder=3) +# + nb = 20 + brsb = ax.barbs(to_np(lons[::nb,::nb]), to_np(lats[::nb,::nb]), + to_np(u[::nb,::nb]), to_np(v[::nb,::nb]), + transform=ccrs.PlateCarree(), + length=4.5, zorder=4, linewidth=0.3) +# + ax.set_xlim(cartopy_lim_buffered(temp, "x")) + ax.set_ylim(cartopy_lim_buffered(temp, "y")) +# + cbar = plt.colorbar(con, cax=ax_cbarv) + add_title_and_description(nc, "{:d}".format(plevel) + " hPa level Geopotential [gpdm], Temperature [°C] and Wind \n" + "WRF {init:s} +{fcsthour:d}", ax) +# + plt.savefig(output_path, bbox_inches='tight', pad_inches=0.1) + plt.close() + +# https://stackoverflow.com/questions/37327308/add-alpha-to-an-existing-matplotlib-colormap +def alphafy_colormap(name, mina=0.0, maxa=0.8): + # Choose colormap + cmap = matplotlib.cm.get_cmap(name) + + # Get the colormap colors + my_cmap = cmap(np.arange(cmap.N)) + + # Set alpha + my_cmap[:,-1] = np.linspace(mina, maxa, cmap.N) + + # Create new colormap + return ListedColormap(my_cmap) + + +def sfc_map(t, nc, get_var_fun, levels, ticks, label, bg_function, output_path, cmap="YlOrRd", extend='both'): +# + var = get_var_fun(nc) + u, v = getvar(nc, 'uvmet10') +# + # Get the latitude and longitude points + lats, lons = latlon_coords(u) +# + # Get the WRF projection (CRS) object + crs = get_cartopy(u) +# + # Create a figure + fig, ax, ax_cbarv, ax_cbarh = make_fig_and_ax(plt, crs) +# + # Draw the contours and filled contours + lab = np.array(levels) + con = ax.contourf(to_np(lons), to_np(lats), var, levels=lab, + transform=ccrs.PlateCarree(), + cmap=alphafy_colormap(cmap), zorder=1, extend=extend) +# + linewidths = np.ones(lab.shape)/2 + linewidths[np.where(lab == 0)[0]] = 1 + cs = ax.contour(to_np(lons), to_np(lats), var, levels=lab[::2], + linewidths=linewidths, colors="black", alpha=0.5, + zorder=2, + transform=ccrs.PlateCarree()) + plt.clabel(cs, fmt='%1.0f', fontsize=8) +# + # Add background (needs to happen after map is already created...) + bg_function(ax, crs) +# + nb = 10 + brsb = ax.barbs(to_np(lons[::nb,::nb]), to_np(lats[::nb,::nb]), + to_np(u[::nb,::nb]), to_np(v[::nb,::nb]), + transform=ccrs.PlateCarree(), + length=4.5, zorder=4, linewidth=0.5) + + xlim = cartopy_lim_buffered(u, "x") + ylim = cartopy_lim_buffered(u, "y") +# + ax.set_xlim(xlim) + ax.set_ylim(ylim) +# + cbar = plt.colorbar(con, cax=ax_cbarv) + cbar.set_ticks(ticks) +# + add_title_and_description(nc, label + " and Wind (10m) \n" + "WRF {init:s} +{fcsthour:d}", ax) +# + plt.savefig(output_path, bbox_inches='tight', pad_inches=0.1) + plt.close() diff --git a/blueprints/highres_chemistry/pp/plotting/maps/pois.py b/blueprints/highres_chemistry/pp/plotting/maps/pois.py new file mode 100644 index 0000000000000000000000000000000000000000..68244521e48e31c252db0ee3bdb74106d6c2ff1b --- /dev/null +++ b/blueprints/highres_chemistry/pp/plotting/maps/pois.py @@ -0,0 +1,268 @@ +southern_germany = [ { 'name':'Munich', 'lon':11.582, 'lat':48.135}, + { 'name':'Salzburg', 'lon':13.0550, 'lat':47.8095}, + { 'name':'Augsburg', 'lon':10.8978, 'lat':48.3705}, + { 'name':'Passau', 'lon':13.4319, 'lat':48.5667}, + { 'name':'Rosenheim', 'lon':12.1181, 'lat':47.8571}, + { 'name':'Ingolstadt', 'lon':11.4258, 'lat':48.7665}, + { 'name':'Garmisch-Partenkirchen', 'lon':11.0948, 'lat':47.4919}, + { 'name':'Ulm', 'lon':9.9876, 'lat':48.4011}, + { 'name':'Stuttgart', 'lon':9.1829, 'lat':48.7758}, + { 'name':'Konstanz', 'lon':9.1737, 'lat':47.6780}, + { 'name':'Zürich', 'lon':8.5417, 'lat':47.3769}, + { 'name':'Nürnberg', 'lon':11.0767, 'lat':49.4521}, + { 'name':'Regensburg', 'lon':12.101624, 'lat':49.013432}, + { 'name':'Innsbruck', 'lon':11.388397, 'lat':47.258320} ] + + +capitals = [ { 'name':'Sukhumi', 'lat':43.001525, 'lon':41.023415}, + { 'name':'Kabul', 'lat':34.575503, 'lon':69.240073}, + { 'name':'Mariehamn', 'lat':60.1, 'lon':19.933333}, + { 'name':'Tirana', 'lat':41.327546, 'lon':19.818698}, + { 'name':'Algiers', 'lat':36.752887, 'lon':3.042048}, + { 'name':'Pago Pago', 'lat':-14.275632, 'lon':-170.702036}, + { 'name':'Andorra la Vella', 'lat':42.506317, 'lon':1.521835}, + { 'name':'Luanda', 'lat':-8.839988, 'lon':13.289437}, + { 'name':'The Valley', 'lat':18.214813, 'lon':-63.057441}, + { 'name':'South Pole', 'lat':-90, 'lon':0}, + { 'name':'St. John\'s', 'lat':17.12741, 'lon':-61.846772}, + { 'name':'Buenos Aires', 'lat':-34.603684, 'lon':-58.381559}, + { 'name':'Yerevan', 'lat':40.179186, 'lon':44.499103}, + { 'name':'Oranjestad', 'lat':12.509204, 'lon':-70.008631}, + { 'name':'Canberra', 'lat':-35.282, 'lon':149.128684}, + { 'name':'Vienna', 'lat':48.208174, 'lon':16.373819}, + { 'name':'Baku', 'lat':40.409262, 'lon':49.867092}, + { 'name':'Nassau', 'lat':25.047984, 'lon':-77.355413}, + { 'name':'Manama', 'lat':26.228516, 'lon':50.58605}, + { 'name':'Dhaka', 'lat':23.810332, 'lon':90.412518}, + { 'name':'Bridgetown', 'lat':13.113222, 'lon':-59.598809}, + { 'name':'Minsk', 'lat':53.90454, 'lon':27.561524}, + { 'name':'Brussels', 'lat':50.85034, 'lon':4.35171}, + { 'name':'Belmopan', 'lat':17.251011, 'lon':-88.75902}, + { 'name':'Porto-Novo', 'lat':6.496857, 'lon':2.628852}, + { 'name':'Hamilton', 'lat':32.294816, 'lon':-64.781375}, + { 'name':'Thimphu', 'lat':27.472792, 'lon':89.639286}, + { 'name':'La Paz', 'lat':-16.489689, 'lon':-68.119294}, + { 'name':'Sarajevo', 'lat':43.856259, 'lon':18.413076}, + { 'name':'Gaborone', 'lat':-24.628208, 'lon':25.923147}, + { 'name':'Bouvet Island', 'lat':-54.43, 'lon':3.38}, + { 'name':'Brasília', 'lat':-15.794229, 'lon':-47.882166}, + { 'name':'Camp Justice', 'lat':21.3419, 'lon':55.4778}, + { 'name':'Road Town', 'lat':18.428612, 'lon':-64.618466}, + { 'name':'Bandar Seri Begawan', 'lat':4.903052, 'lon':114.939821}, + { 'name':'Sofia', 'lat':42.697708, 'lon':23.321868}, + { 'name':'Ouagadougou', 'lat':12.371428, 'lon':-1.51966}, + { 'name':'Bujumbura', 'lat':-3.361378, 'lon':29.359878}, + { 'name':'Phnom Penh', 'lat':11.544873, 'lon':104.892167}, + { 'name':'Yaoundé', 'lat':3.848033, 'lon':11.502075}, + { 'name':'Ottawa', 'lat':45.42153, 'lon':-75.697193}, + { 'name':'Praia', 'lat':14.93305, 'lon':-23.513327}, + { 'name':'George Town', 'lat':19.286932, 'lon':-81.367439}, + { 'name':'Bangui', 'lat':4.394674, 'lon':18.55819}, + { 'name':'N\'Djamena', 'lat':12.134846, 'lon':15.055742}, + { 'name':'Santiago', 'lat':-33.44889, 'lon':-70.669265}, + { 'name':'Beijing', 'lat':39.904211, 'lon':116.407395}, + { 'name':'Flying Fish Cove', 'lat':-10.420686, 'lon':105.679379}, + { 'name':'West Island', 'lat':-12.188834, 'lon':96.829316}, + { 'name':'Bogotá', 'lat':4.710989, 'lon':-74.072092}, + { 'name':'Moroni', 'lat':-11.717216, 'lon':43.247315}, + { 'name':'Kinshasa', 'lat':-4.441931, 'lon':15.266293}, + { 'name':'Brazzaville', 'lat':-4.26336, 'lon':15.242885}, + { 'name':'Avarua', 'lat':-21.212901, 'lon':-159.782306}, + { 'name':'San José', 'lat':9.928069, 'lon':-84.090725}, + { 'name':'Yamoussoukro', 'lat':6.827623, 'lon':-5.289343}, + { 'name':'Zagreb', 'lat':45.815011, 'lon':15.981919}, + { 'name':'Havana', 'lat':23.05407, 'lon':-82.345189}, + { 'name':'Willemstad', 'lat':12.122422, 'lon':-68.882423}, + { 'name':'Nicosia', 'lat':35.185566, 'lon':33.382276}, + { 'name':'Prague', 'lat':50.075538, 'lon':14.4378}, + { 'name':'Copenhagen', 'lat':55.676097, 'lon':12.568337}, + { 'name':'Djibouti', 'lat':11.572077, 'lon':43.145647}, + { 'name':'Roseau', 'lat':15.309168, 'lon':-61.379355}, + { 'name':'Santo Domingo', 'lat':18.486058, 'lon':-69.931212}, + { 'name':'Quito', 'lat':-0.180653, 'lon':-78.467838}, + { 'name':'Cairo', 'lat':30.04442, 'lon':31.235712}, + { 'name':'San Salvador', 'lat':13.69294, 'lon':-89.218191}, + { 'name':'Malabo', 'lat':3.750412, 'lon':8.737104}, + { 'name':'Asmara', 'lat':15.322877, 'lon':38.925052}, + { 'name':'Tallinn', 'lat':59.436961, 'lon':24.753575}, + { 'name':'Addis Ababa', 'lat':8.980603, 'lon':38.757761}, + { 'name':'Stanley', 'lat':-51.697713, 'lon':-57.851663}, + { 'name':'Tórshavn', 'lat':62.007864, 'lon':-6.790982}, + { 'name':'Suva', 'lat':-18.124809, 'lon':178.450079}, + { 'name':'Helsinki', 'lat':60.173324, 'lon':24.941025}, + { 'name':'Paris', 'lat':48.856614, 'lon':2.352222}, + { 'name':'Cayenne', 'lat':4.92242, 'lon':-52.313453}, + { 'name':'Papeete', 'lat':-17.551625, 'lon':-149.558476}, + { 'name':'Saint-Pierre ', 'lat':-21.3419, 'lon':55.4778}, + { 'name':'Libreville', 'lat':0.416198, 'lon':9.467268}, + { 'name':'Banjul', 'lat':13.454876, 'lon':-16.579032}, + { 'name':'Tbilisi', 'lat':41.715138, 'lon':44.827096}, + { 'name':'Berlin', 'lat':52.520007, 'lon':13.404954}, + { 'name':'Accra', 'lat':5.603717, 'lon':-0.186964}, + { 'name':'Gibraltar', 'lat':36.140773, 'lon':-5.353599}, + { 'name':'Athens', 'lat':37.983917, 'lon':23.72936}, + { 'name':'Nuuk', 'lat':64.18141, 'lon':-51.694138}, + { 'name':'St. George\'s', 'lat':12.056098, 'lon':-61.7488}, + { 'name':'Basse-Terre', 'lat':16.014453, 'lon':-61.706411}, + { 'name':'Hagåtña', 'lat':13.470891, 'lon':144.751278}, + { 'name':'Guatemala City', 'lat':14.634915, 'lon':-90.506882}, + { 'name':'St. Peter Port', 'lat':49.455443, 'lon':-2.536871}, + { 'name':'Conakry', 'lat':9.641185, 'lon':-13.578401}, + { 'name':'Bissau', 'lat':11.881655, 'lon':-15.617794}, + { 'name':'Georgetown', 'lat':6.801279, 'lon':-58.155125}, + { 'name':'Port-au-Prince', 'lat':18.594395, 'lon':-72.307433}, + { 'name':'Tegucigalpa', 'lat':14.072275, 'lon':-87.192136}, + { 'name':'Hong Kong', 'lat':22.396428, 'lon':114.109497}, + { 'name':'Budapest', 'lat':47.497912, 'lon':19.040235}, + { 'name':'Reykjavík', 'lat':64.126521, 'lon':-21.817439}, + { 'name':'New Delhi', 'lat':28.613939, 'lon':77.209021}, + { 'name':'Jakarta', 'lat':-6.208763, 'lon':106.845599}, + { 'name':'Tehran', 'lat':35.689198, 'lon':51.388974}, + { 'name':'Baghdad', 'lat':33.312806, 'lon':44.361488}, + { 'name':'Dublin', 'lat':53.349805, 'lon':-6.26031}, + { 'name':'Douglas', 'lat':54.152337, 'lon':-4.486123}, + { 'name':'Tel Aviv', 'lat':32.0853, 'lon':34.781768}, + { 'name':'Rome', 'lat':41.902784, 'lon':12.496366}, + { 'name':'Kingston', 'lat':18.042327, 'lon':-76.802893}, + { 'name':'Tokyo', 'lat':35.709026, 'lon':139.731992}, + { 'name':'St. Helier', 'lat':49.186823, 'lon':-2.106568}, + { 'name':'Amman', 'lat':31.956578, 'lon':35.945695}, + { 'name':'Astana', 'lat':51.160523, 'lon':71.470356}, + { 'name':'Nairobi', 'lat':-1.292066, 'lon':36.821946}, + { 'name':'Tarawa Atoll', 'lat':1.451817, 'lon':172.971662}, + { 'name':'Pristina', 'lat':42.662914, 'lon':21.165503}, + { 'name':'Kuwait City', 'lat':29.375859, 'lon':47.977405}, + { 'name':'Bishkek', 'lat':42.874621, 'lon':74.569762}, + { 'name':'Vientiane', 'lat':17.975706, 'lon':102.633104}, + { 'name':'Riga', 'lat':56.949649, 'lon':24.105186}, + { 'name':'Beirut', 'lat':33.888629, 'lon':35.495479}, + { 'name':'Maseru', 'lat':-29.363219, 'lon':27.51436}, + { 'name':'Monrovia', 'lat':6.290743, 'lon':-10.760524}, + { 'name':'Tripoli', 'lat':32.887209, 'lon':13.191338}, + { 'name':'Vaduz', 'lat':47.14103, 'lon':9.520928}, + { 'name':'Vilnius', 'lat':54.687156, 'lon':25.279651}, + { 'name':'Luxembourg', 'lat':49.611621, 'lon':6.131935}, + { 'name':'Macau', 'lat':22.166667, 'lon':113.55}, + { 'name':'Skopje', 'lat':41.997346, 'lon':21.427996}, + { 'name':'Antananarivo', 'lat':-18.87919, 'lon':47.507905}, + { 'name':'Lilongwe', 'lat':-13.962612, 'lon':33.774119}, + { 'name':'Kuala Lumpur', 'lat':3.139003, 'lon':101.686855}, + { 'name':'Malé', 'lat':4.175496, 'lon':73.509347}, + { 'name':'Bamako', 'lat':12.639232, 'lon':-8.002889}, + { 'name':'Valletta', 'lat':35.898909, 'lon':14.514553}, + { 'name':'Majuro', 'lat':7.116421, 'lon':171.185774}, + { 'name':'Fort-de-France', 'lat':14.616065, 'lon':-61.05878}, + { 'name':'Nouakchott', 'lat':18.07353, 'lon':-15.958237}, + { 'name':'Port Louis', 'lat':-20.166896, 'lon':57.502332}, + { 'name':'Mamoudzou', 'lat':-12.780949, 'lon':45.227872}, + { 'name':'Mexico City', 'lat':19.432608, 'lon':-99.133208}, + { 'name':'Palikir', 'lat':6.914712, 'lon':158.161027}, + { 'name':'Chisinau', 'lat':47.010453, 'lon':28.86381}, + { 'name':'Monaco', 'lat':43.737411, 'lon':7.420816}, + { 'name':'Ulaanbaatar', 'lat':47.886399, 'lon':106.905744}, + { 'name':'Podgorica', 'lat':42.43042, 'lon':19.259364}, + { 'name':'Plymouth', 'lat':16.706523, 'lon':-62.215738}, + { 'name':'Rabat', 'lat':33.97159, 'lon':-6.849813}, + { 'name':'Maputo', 'lat':-25.891968, 'lon':32.605135}, + { 'name':'Naypyidaw', 'lat':19.763306, 'lon':96.07851}, + { 'name':'Stepanakert', 'lat':39.826385, 'lon':46.763595}, + { 'name':'Windhoek', 'lat':-22.560881, 'lon':17.065755}, + { 'name':'Yaren', 'lat':-0.546686, 'lon':166.921091}, + { 'name':'Kathmandu', 'lat':27.717245, 'lon':85.323961}, + { 'name':'Amsterdam', 'lat':52.370216, 'lon':4.895168}, + { 'name':'Willemstad ', 'lat':12.1091242, 'lon':-68.9316546}, + { 'name':'Nouméa', 'lat':-22.255823, 'lon':166.450524}, + { 'name':'Wellington', 'lat':-41.28646, 'lon':174.776236}, + { 'name':'Managua', 'lat':12.114993, 'lon':-86.236174}, + { 'name':'Niamey', 'lat':13.511596, 'lon':2.125385}, + { 'name':'Abuja', 'lat':9.076479, 'lon':7.398574}, + { 'name':'Alofi', 'lat':-19.055371, 'lon':-169.917871}, + { 'name':'Kingston', 'lat':-29.056394, 'lon':167.959588}, + { 'name':'Pyongyang', 'lat':39.039219, 'lon':125.762524}, + { 'name':'Nicosia', 'lat':35.185566, 'lon':33.382276}, + { 'name':'Saipan', 'lat':15.177801, 'lon':145.750967}, + { 'name':'Oslo', 'lat':59.913869, 'lon':10.752245}, + { 'name':'Muscat', 'lat':23.58589, 'lon':58.405923}, + { 'name':'Islamabad', 'lat':33.729388, 'lon':73.093146}, + { 'name':'Ngerulmud', 'lat':7.500384, 'lon':134.624289}, + { 'name':'Ramallah', 'lat':31.9073509, 'lon':35.5354719}, + { 'name':'Panama City', 'lat':9.101179, 'lon':-79.402864}, + { 'name':'Port Moresby', 'lat':-9.4438, 'lon':147.180267}, + { 'name':'Asuncion', 'lat':-25.26374, 'lon':-57.575926}, + { 'name':'Lima', 'lat':-12.046374, 'lon':-77.042793}, + { 'name':'Manila', 'lat':14.599512, 'lon':120.98422}, + { 'name':'Adamstown', 'lat':-25.06629, 'lon':-130.100464}, + { 'name':'Warsaw', 'lat':52.229676, 'lon':21.012229}, + { 'name':'Lisbon', 'lat':38.722252, 'lon':-9.139337}, + { 'name':'San Juan', 'lat':18.466334, 'lon':-66.105722}, + { 'name':'Doha', 'lat':25.285447, 'lon':51.53104}, + { 'name':'Saint-Denis', 'lat':-20.882057, 'lon':55.450675}, + { 'name':'Bucharest', 'lat':44.426767, 'lon':26.102538}, + { 'name':'Moscow', 'lat':55.755826, 'lon':37.6173}, + { 'name':'Kigali', 'lat':-1.957875, 'lon':30.112735}, + { 'name':'St. Pierre', 'lat':46.775846, 'lon':-56.180636}, + { 'name':'Kingstown', 'lat':13.160025, 'lon':-61.224816}, + { 'name':'Apia', 'lat':-13.850696, 'lon':-171.751355}, + { 'name':'San Marino', 'lat':43.935591, 'lon':12.447281}, + { 'name':'São Tomé', 'lat':0.330192, 'lon':6.733343}, + { 'name':'Riyadh', 'lat':24.749403, 'lon':46.902838}, + { 'name':'Dakar', 'lat':14.764504, 'lon':-17.366029}, + { 'name':'Belgrade', 'lat':44.786568, 'lon':20.448922}, + { 'name':'Victoria', 'lat':-4.619143, 'lon':55.451315}, + { 'name':'Freetown', 'lat':8.465677, 'lon':-13.231722}, + { 'name':'Singapore', 'lat':1.280095, 'lon':103.850949}, + { 'name':'Bratislava', 'lat':48.145892, 'lon':17.107137}, + { 'name':'Ljubljana', 'lat':46.056947, 'lon':14.505751}, + { 'name':'Honiara', 'lat':-9.445638, 'lon':159.9729}, + { 'name':'Mogadishu', 'lat':2.046934, 'lon':45.318162}, + { 'name':'Pretoria', 'lat':-25.747868, 'lon':28.229271}, + { 'name':'King Edward Point', 'lat':-54.28325, 'lon':-36.493735}, + { 'name':'Seoul', 'lat':37.566535, 'lon':126.977969}, + { 'name':'Tskhinvali', 'lat':42.22146, 'lon':43.964405}, + { 'name':'Juba', 'lat':4.859363, 'lon':31.57125}, + { 'name':'Madrid', 'lat':40.416775, 'lon':-3.70379}, + { 'name':'Sri Jayawardenepura Kotte', 'lat':6.89407, 'lon':79.902478}, + { 'name':'Gustavia', 'lat':17.896435, 'lon':-62.852201}, + { 'name':'Basseterre', 'lat':17.302606, 'lon':-62.717692}, + { 'name':'Castries', 'lat':14.010109, 'lon':-60.987469}, + { 'name':'Marigot', 'lat':18.067519, 'lon':-63.082466}, + { 'name':'Khartoum', 'lat':15.500654, 'lon':32.559899}, + { 'name':'Paramaribo', 'lat':5.852036, 'lon':-55.203828}, + { 'name':'Longyearbyen ', 'lat':78.062, 'lon':22.055}, + { 'name':'Mbabane', 'lat':-26.305448, 'lon':31.136672}, + { 'name':'Stockholm', 'lat':59.329323, 'lon':18.068581}, + { 'name':'Bern', 'lat':46.947974, 'lon':7.447447}, + { 'name':'Damascus', 'lat':33.513807, 'lon':36.276528}, + { 'name':'Taipei', 'lat':25.032969, 'lon':121.565418}, + { 'name':'Dushanbe', 'lat':38.559772, 'lon':68.787038}, + { 'name':'Dodoma', 'lat':-6.162959, 'lon':35.751607}, + { 'name':'Bangkok', 'lat':13.756331, 'lon':100.501765}, + { 'name':'Dili', 'lat':-8.556856, 'lon':125.560314}, + { 'name':'Lomé', 'lat':6.172497, 'lon':1.231362}, + { 'name':'Nukunonu', 'lat':-9.2005, 'lon':-171.848}, + { 'name':'Nuku\'alofa', 'lat':-21.139342, 'lon':-175.204947}, + { 'name':'Tiraspol', 'lat':46.848185, 'lon':29.596805}, + { 'name':'Port of Spain', 'lat':10.654901, 'lon':-61.501926}, + { 'name':'Edinburgh of the Seven Seas', 'lat':-37.068042, 'lon':-12.311315}, + { 'name':'Tunis', 'lat':36.806495, 'lon':10.181532}, + { 'name':'Ankara', 'lat':39.933364, 'lon':32.859742}, + { 'name':'Ashgabat', 'lat':37.960077, 'lon':58.326063}, + { 'name':'Cockburn Town', 'lat':21.467458, 'lon':-71.13891}, + { 'name':'Funafuti', 'lat':-8.520066, 'lon':179.198128}, + { 'name':'Charlotte Amalie', 'lat':18.3419, 'lon':-64.930701}, + { 'name':'Kampala', 'lat':0.347596, 'lon':32.58252}, + { 'name':'Kiev', 'lat':50.4501, 'lon':30.5234}, + { 'name':'Abu Dhabi', 'lat':24.299174, 'lon':54.697277}, + { 'name':'London', 'lat':51.507351, 'lon':-0.127758}, + { 'name':'Washington', 'lat':38.907192, 'lon':-77.036871}, + { 'name':'Montevideo', 'lat':-34.901113, 'lon':-56.164531}, + { 'name':'Tashkent', 'lat':41.299496, 'lon':69.240073}, + { 'name':'Port Vila', 'lat':-17.733251, 'lon':168.327325}, + { 'name':'Vatican City', 'lat':41.902179, 'lon':12.453601}, + { 'name':'Caracas', 'lat':10.480594, 'lon':-66.903606}, + { 'name':'Hanoi', 'lat':21.027764, 'lon':105.83416}, + { 'name':'Mata-Utu', 'lat':-13.282509, 'lon':-176.176447}, + { 'name':'El Aaiún', 'lat':27.125287, 'lon':-13.1625}, + { 'name':'Sana\'a', 'lat':15.369445, 'lon':44.191007}, + { 'name':'Lusaka', 'lat':-15.387526, 'lon':28.322817}, + { 'name':'Harare', 'lat':-17.825166, 'lon':31.03351} ] \ No newline at end of file diff --git a/blueprints/highres_chemistry/pp/plotting/maps/run.py b/blueprints/highres_chemistry/pp/plotting/maps/run.py new file mode 100644 index 0000000000000000000000000000000000000000..c9820b99942b44bdd379d2617eaf091122fbc107 --- /dev/null +++ b/blueprints/highres_chemistry/pp/plotting/maps/run.py @@ -0,0 +1,110 @@ +import datetime + +from argparse import ArgumentParser, ArgumentTypeError + +def valid_date(s): + try: + return datetime.datetime.strptime(s, "%Y-%m-%d") + except ValueError: + raise ArgumentTypeError("Not a valid date: '{0}'.".format(s)) + +parser = ArgumentParser() +parser.add_argument('wrf_data_fpath_pattern', type=str, help='WRF input data path pattern (including strptime pattern for date and {domain:02d} for domain (01, 02, ...))') +parser.add_argument('date', type=valid_date, help='Date to process (YYYY-mm-dd)') +parser.add_argument('interval_hours', type=int, help='Time interval in hours for plotting') +parser.add_argument('fcst_time_hours', type=int, help='Forecast time in hours') +parser.add_argument('plot_fpath_pattern', type=str, help='Plot output path pattern (including strptime pattern for date, {plot:s} for plot identifier (e.g., 500hPa...) and {domain:02d} for domain (01, 02)') + +args = parser.parse_args() + +from plots import plevel_map, sfc_map, slp_map +import pois +import netCDF4 +from wrf import getvar +import numpy as np +import cartopy.feature as cfeature +import contextily +import rasterio + +times = [ args.date + datetime.timedelta(hours=t) for t in range(0, args.fcst_time_hours, args.interval_hours) ] + +# for calculating delta precip +prec_prev = None + +def bg_function(ax, wrfcrs): + wrfcrs_proj4 = wrfcrs.proj4_init + wrfcrs_proj4 = wrfcrs_proj4.replace('+nadgrids=@null', '') + rasteriocrs = rasterio.crs.CRS.from_proj4(wrfcrs_proj4) + contextily.add_basemap(ax, source=contextily.providers.Stamen.TonerLite, crs=rasteriocrs, zorder=0) + +for t in times: + input_path_d01 = t.strftime(args.wrf_data_fpath_pattern).format(domain=1) + input_path_d02 = t.strftime(args.wrf_data_fpath_pattern).format(domain=2) + + try: + nc_d01 = netCDF4.Dataset(input_path_d01) + nc_d02 = netCDF4.Dataset(input_path_d02) + except Exception as e: + print("Not all data available for step {:s}: {:s}".format(t.strftime("%Y-%m-%d %H"), str(e))) + continue + + print(t) + + # 500 hPa Europe domain +# output_path = t.strftime(args.plot_fpath_pattern).format(plot='500hPa', domain=1) +# plevel_map(t, nc_d01, output_path, 500) + + # Surface Europe domain + # temperature +# sfc_map(t, nc_d01, lambda nc: getvar(nc, 'T2') - 273.15, +# levels=[-25, -20, -15, -12.5, -10, -7.5, -5, -2.5, 0, 2.5, 5, 7.5, 10, 12.5, 15, 17.5, 20, 22.5, 25, 27.5, 30, 32.5, 35], +# ticks=[-25, -15, -10, -5, 0, 5, 10, 15, 20, 25, 30, 35], +# label="Temperature (2m) [°C]", pois=pois.capitals, +# output_path=t.strftime(args.plot_fpath_pattern).format(plot='T', domain=1)) + + # Ozone + for domain, nc_d, poilist in zip([1,2],[nc_d01, nc_d02], [pois.capitals, pois.southern_germany]): + sfc_map(t, nc_d, lambda nc: getvar(nc, 'o3')[0,:,:] * 1e3, + levels=range(0, 160, 10), + ticks=range(0, 150, 20), + label="Ozone (ground level) [ppbv]", bg_function=bg_function, + output_path=t.strftime(args.plot_fpath_pattern).format(plot='O3', domain=domain), + extend='max') + # NO2 + sfc_map(t, nc_d, lambda nc: getvar(nc, 'no2')[0,:,:] * 1e3, + levels=range(0, 110, 10), + ticks=range(0, 110, 20), + label="NO$_2$ (ground level) [ppbv]", bg_function=bg_function, + output_path=t.strftime(args.plot_fpath_pattern).format(plot='NO2', domain=domain), + extend='max') + # PM + def get_pm(nc, bins=[ 1, 2, 3 ]) : + pmspecs = [ "so4", "no3", "smpa", "smpbb", "glysoa_sfc", "biog1_c", "biog1_o", "cl", "co3", "nh4", "na", "cl", "oin", "oc", "bc", "water" ] # ug kg-1 + alt = getvar(nc, "ALT")[0,:,:] # inverse density m3 kg-1 + all = [ getvar(nc, spec + "_a{:02d}".format(bin))[0,:,:] / alt for spec in pmspecs for bin in bins ] + return np.sum(all, axis=0) + +# sfc_map(t, nc_d, lambda nc: get_pm(nc, bins=[1]), +# levels=range(0, 130, 10), +# ticks=range(0, 130, 20), +# label="UFP (ground level) [$\mu$g m$^{{-3}}$]", bg_function=bg_function, +# output_path=t.strftime(args.plot_fpath_pattern).format(plot='UFP', domain=domain), +# extend='max') + sfc_map(t, nc_d, lambda nc: get_pm(nc, bins=[1,2,3]), + levels=range(0, 130, 10), + ticks=range(0, 130, 20), + label="PM$_{{2.5}}$ (ground level) [$\mu$g m$^{{-3}}$]", bg_function=bg_function, + output_path=t.strftime(args.plot_fpath_pattern).format(plot='PM2_5', domain=domain), + extend='max') + sfc_map(t, nc_d, lambda nc: get_pm(nc, bins=[1,2,3,4]), + levels=range(0, 130, 10), + ticks=range(0, 130, 20), + label="PM$_{{10}}$ (ground level) [$\mu$g m$^{{-3}}$]", bg_function=bg_function, + output_path=t.strftime(args.plot_fpath_pattern).format(plot='PM10', domain=domain), + extend='max') + + # SLP Europe +# output_path = t.strftime(args.plot_fpath_pattern).format(plot='SLP', domain=1) +# prec_prev = slp_map(t, nc_d01, prec_prev, output_path) + + nc_d01.close() diff --git a/blueprints/highres_chemistry/pp/plotting/maps/run.sh b/blueprints/highres_chemistry/pp/plotting/maps/run.sh new file mode 100644 index 0000000000000000000000000000000000000000..248dc5b5b142ce64455a10d6f2baef35fed36e73 --- /dev/null +++ b/blueprints/highres_chemistry/pp/plotting/maps/run.sh @@ -0,0 +1,50 @@ +#!/bin/bash -l +#SBATCH --partition=alcc1,alcc3,epyc +#SBATCH -o /alcc/gpfs2/scratch/mbees/knotechr/WRF/postprocessing/maps/%j.%N.out +#SBATCH -J maps +#SBATCH --ntasks=1 +#SBATCH --mem=10G +#SBATCH --mail-type=FAIL +#SBATCH --mail-user=christoph.knote@med.uni-augsburg.de +#SBATCH --time=01:00:00 + +module load gnu geos ffmpeg + +export MPLBACKEND="agg" + +scriptPath=/alcc/gpfs2/home/u/knotechr/wrfotron/blueprints/operational_chemistry/pp/plotting/maps/run.py +wrfDataPathPattern="/alcc/gpfs2/scratch/mbees/knotechr/archive/WRF/operational_chemistry/wrfout_d{domain:02d}_%Y-%m-%d_%H:%M:%S" +tmpWorkDir="/alcc/gpfs2/scratch/mbees/knotechr/WRF/postprocessing/maps/movie" +plotPathPattern="${tmpWorkDir}/{plot:s}_d{domain:02d}_%Y%m%d%H.png" +webDir="/alcc/gpfs2/scratch/mbees/knotechr/plots/WRF/operational_chemistry/__domain__/maps/%Y/%m" +dateToProcess="today" +fcstIntervalHours=1 +fcstRangeHours=72 +nSecondsWait=2 + +mkdir -p ${tmpWorkDir} +rm -f ${tmpWorkDir}/* + +python3 $scriptPath ${wrfDataPathPattern} $(date -u --date="${dateToProcess}" "+%Y-%m-%d") ${fcstIntervalHours} ${fcstRangeHours} ${plotPathPattern} + +frameRate=3 + +for domain in d01 d02 +do + # get all plot type names ("500hpa_d01, o3_d02", ...) + typs=$(ls -1 ${tmpWorkDir}/*${domain}*.png | sed -e "s|${tmpWorkDir}/||" | sed -E -e "s/_[0-9]+.png//g" | uniq) + + for typ in ${typs} + do + echo "" + echo $typ $domain + outFilePath=$(date --date="${dateToProcess} 0 UTC" "+${webDir/__domain__/${domain}}/${typ/_${domain}/}_%Y%m%d.mp4") + mkdir -p $(dirname $outFilePath) + rm -f ${outFilePath} + # see here for why we need vf argument https://stackoverflow.com/questions/20847674/ffmpeg-libx264-height-not-divisible-by-2 + cat ${tmpWorkDir}/${typ}_*.png | ffmpeg -y -r ${frameRate} -f image2pipe -i - -vcodec libx264 -crf 18 -preset slow -pix_fmt yuv420p -vf "pad=ceil(iw/2)*2:ceil(ih/2)*2, tpad=stop_mode=clone:stop_duration=2, fps=fps=${frameRate}" ${outFilePath} + done +done + +rm -f ${tmpWorkDir}/* + diff --git a/blueprints/highres_chemistry/pp/plotting/maps/tools.py b/blueprints/highres_chemistry/pp/plotting/maps/tools.py new file mode 100644 index 0000000000000000000000000000000000000000..b2a7fdf244453f5afd16089d3842b9c68641624d --- /dev/null +++ b/blueprints/highres_chemistry/pp/plotting/maps/tools.py @@ -0,0 +1,68 @@ +from wrf import cartopy_xlim, cartopy_ylim +import matplotlib.gridspec as gridspec + +from matplotlib.offsetbox import AnchoredText + +import datetime +import numpy as np + +t2datetime = lambda x: datetime.datetime.strptime(x, '%Y-%m-%d_%H:%M:%S') + +def wrf_time(nc): + first = ''.join([x.decode('utf-8') for x in nc['Times'][0]]) + return t2datetime(first) + +def total_prec(nc, prec_p): + prec_vars = ['RAINC', 'RAINSH', 'RAINNC', 'SNOWNC', 'GRAUPELNC', 'HAILNC'] + start = t2datetime(nc.START_DATE) + prec = np.sum([nc[var][0, :, :] for var in prec_vars], axis=0) +# + if type(prec_p) == type(None): + prec_p = prec +# + if start == wrf_time(nc): + return prec + else: + return prec-prec_p + +# buffer 0.1 -> 10% border around the edges +def cartopy_lim_buffered(field, direction, buffer=0.1): + if direction == "x": + lims = cartopy_xlim(field) + elif direction == "y": + lims = cartopy_ylim(field) + else: + raise Exception("Direction not understood!") + dlim = lims[1] - lims[0] + lims[0] += 0.5 * buffer * dlim + lims[1] -= 0.5 * buffer * dlim + return lims + +def make_fig_and_ax(plt, crs, add_horiz_cbar=False): + fig = plt.figure(figsize=(7,7)) +# + gs = gridspec.GridSpec(2, 2, width_ratios=[20,1], height_ratios=[20,1], wspace=0.05, hspace=0.05) +# + ax = fig.add_subplot(gs[0,0], projection=crs) + ax_cbarv = fig.add_subplot(gs[0,1]) + ax_cbarh = None + if add_horiz_cbar: + ax_cbarh = fig.add_subplot(gs[1,0]) + return fig, ax, ax_cbarv, ax_cbarh + +def add_title_and_description(nc, title, ax): + ''' + title: including {init:s} for timestamp of run initialization, and {fcsthour:d} for forecast hour + ''' + fcstHour = int((wrf_time(nc) - t2datetime(nc.START_DATE)).total_seconds()/3600.) +# + ax.set_title(title.format(init=t2datetime(nc.START_DATE).strftime("%d.%m.%Y %H:%M UTC"), + fcsthour=fcstHour), fontdict={ 'fontsize':8, 'horizontalalignment': 'left'}, loc='left') +# + texttime = AnchoredText(wrf_time(nc).strftime("%H UTC %a"), + loc=2, pad=0.2, prop={'size': 10}, frameon=True) + ax.add_artist(texttime) +# + textcopyright = AnchoredText("$\copyright$ MBEES, Faculty of Medicine, University of Augsburg", + loc=4, pad=0.1, prop={'size': 6, 'color':'white'}, frameon=False) + ax.add_artist(textcopyright) diff --git a/blueprints/highres_chemistry/pp/plotting/meteograms/plot.py b/blueprints/highres_chemistry/pp/plotting/meteograms/plot.py new file mode 100755 index 0000000000000000000000000000000000000000..278d373363974298fcbb75be02b199cc2c143f57 --- /dev/null +++ b/blueprints/highres_chemistry/pp/plotting/meteograms/plot.py @@ -0,0 +1,757 @@ +#!/usr/bin/env python + +from wrf import getvar, ll_to_xy, to_np, extract_times +import numpy as np + + +def fT2m(D, xmin, xmax, ymin, ymax): + return D["T2"][0, xmin:xmax, ymin:ymax] - 273.15 + + +def fTd2m(D, xmin, xmax, ymin, ymax): + return getvar(D, "td2")[xmin:xmax, ymin:ymax] + + +def wind10m(D, lat, lon): + """Returns the wind speed and the wind direction at 10 m""" + try: + x, y = ll_to_xy(D, lat, lon) + w_speed, w_dir = getvar(D, "uvmet10_wspd_wdir") + return w_speed.values[int(x), int(y)], w_dir.values[int(x), int(y)] + except: + return np.NaN, np.NaN + + +def wind_uv(D, lat, lon): + """Returns the u and v wind components at 10 m""" + try: + x, y = ll_to_xy(D, lat, lon) + w_u, w_v = getvar(D, "uvmet10") + return w_u.values[int(x), int(y)], w_v.values[int(x), int(y)] + except: + return np.NaN, np.NaN + + +def Mslp(D, lat, lon): + try: + """Returns the mean sea level pressure""" + x, y = ll_to_xy(D, lat, lon) + slp = getvar(D, "slp") + return slp.values[int(x), int(y)] + except: + return np.NaN + + +def cin_cape(D, lat, lon): + try: + """Returns the 2D max cin and cape""" + x, y = ll_to_xy(D, lat, lon) + var = getvar(D, "cape_2d") + return (var[1].values[int(x), int(y)], var[0].values[int(x), int(y)]) + except: + return np.NaN, np.NaN + + +def vert_wind(D, lat, lon): + try: + from wrf import vertcross, CoordPair + + x, y = ll_to_xy(D, lat, lon) + xy = CoordPair(x=x, y=y) + u, v = getvar(D, "uvmet") + p = getvar(D, "p", units="hPa") + levels = np.arange(400, 1060, 60) + u_vert = vertcross( + u, + p, + levels=levels, + start_point=xy, + end_point=CoordPair(x=x + 1, y=y + 1), + meta=False, + missing=np.NaN, + ) + v_vert = vertcross( + v, + p, + levels=levels, + start_point=xy, + end_point=CoordPair(x=x + 1, y=y + 1), + meta=False, + missing=np.NaN, + ) + return list(u_vert[:, 0]), list(v_vert[:, 0]) + except: + levels = np.arange(400, 1060, 60) + return [np.NaN for x in range(len(levels))], [ + np.NaN for x in range(len(levels)) + ] + + +def vert_T(D, lat, lon): + try: + from wrf import vertcross, CoordPair + + x, y = ll_to_xy(D, lat, lon) + xy = CoordPair(x=x, y=y) + T = getvar(D, "tc") + p = getvar(D, "p", units="hPa") + levels = np.arange(400, 1060, 60) + T_vert = to_np( + vertcross( + T, + p, + levels=levels, + start_point=xy, + end_point=CoordPair(x=x + 1, y=y + 1), + meta=False, + missing=np.NaN, + ) + ) + return T_vert[:, 0][:] + except: + levels = np.arange(400, 1060, 60) + return [np.NaN for x in range(len(levels))], [ + np.NaN for x in range(len(levels)) + ] + + +def fcfrac(D, lat, lon): + from wrf import interplevel, CoordPair + + levels = np.arange(400, 1030, 30) + try: + x, y = ll_to_xy(D, lat, lon) + p = getvar(D, "p", units="hPa") + cf = interplevel(D["CLDFRA"], p, levels) + return cf.values[:, int(x), int(y)] + + except: + return [np.NaN for x in range(len(levels))] + + +def ftot_accum_rain_2d(D, xmin, xmax, ymin, ymax): + try: + return np.sum( + [ + D[var][0, xmin:xmax, ymin:ymax] + for var in [ + "RAINC", + "RAINSH", + "RAINNC", + "SNOWNC", + "GRAUPELNC", + "HAILNC", + ] + ], + axis=0, + ) + except: + return np.mgrid[xmin:xmax, ymin:ymax][0] * np.NaN + + +def fdrain(DS, xmin, xmax, ymin, ymax): + tot_accum_rain = np.array( + [ftot_accum_rain_2d(D, xmin, xmax, ymin, ymax) for D in DS] + ) + + # weird construct to remove negative rain rates + # because timeseries have restarts in them + # i.e. we search for spots where rainrate accumulation starts + + starts = np.array([t2datetime(D.START_DATE) for D in DS]) + dt = np.array([x.total_seconds() for x in (starts[1:] - starts[:-1])]) + drain = np.zeros_like(tot_accum_rain) + drain[1:] = tot_accum_rain[1:] - tot_accum_rain[:-1] + drain[1:] = np.where( + dt > 0, tot_accum_rain[1:].T, (tot_accum_rain[1:] - tot_accum_rain[:-1]).T + ).T + return drain + + +def ftot_accum_snow_2d(D, xmin, xmax, ymin, ymax): + try: + return np.sum( + [D[var][0, xmin:xmax, ymin:ymax] for var in ["SNOWNC", "GRAUPELNC"]], axis=0 + ) + except: + return np.mgrid[xmin:xmax, ymin:ymax][0] * np.NaN + + +def fdsnow(DS, xmin, xmax, ymin, ymax): + tot_accum_snow = np.array( + [ftot_accum_snow_2d(D, xmin, xmax, ymin, ymax) for D in DS] + ) + + starts = np.array([t2datetime(D.START_DATE) for D in DS]) + dt = np.array([x.total_seconds() for x in (starts[1:] - starts[:-1])]) + drain = np.zeros_like(tot_accum_snow) + drain[1:] = tot_accum_snow[1:] - tot_accum_snow[:-1] + drain[1:] = np.where( + dt > 0, tot_accum_snow[1:].T, (tot_accum_snow[1:] - tot_accum_snow[:-1]).T + ).T + return drain + + +def t2datetime(x): + import datetime + + return datetime.datetime.strptime(x, "%Y-%m-%d_%H:%M:%S") + + +def ftimes(D): + """Convert first WRF timestamp description into datetimes string""" + return t2datetime("".join([x.decode("utf-8") for x in D["Times"][0]])) + + +def corresponding_pixels(D, location_lon, location_lat, deg_round_loc): + import numpy as np + + xlon = D["XLONG"][0] + xlat = D["XLAT"][0] + + around_loc = (np.abs(xlat - location_lat) < deg_round_loc) & ( + np.abs(xlon - location_lon) < deg_round_loc + ) + + xmin, xmax = np.nonzero(around_loc)[0][[0, -1]] + ymin, ymax = np.nonzero(around_loc)[1][[0, -1]] + + xlon = xlon[xmin:xmax, ymin:ymax] + xlat = xlat[xmin:xmax, ymin:ymax] + + # print("Using {} pixels for plot".format(np.size(xlon))) + return (xmin, xmax), (ymin, ymax), xlon, xlat + + +def haversine(lat1, lon1, lat2, lon2): + """ + Compute distance between lat lon points + """ + import numpy as np + + R = 6378137 # Earth radius + dLat = np.deg2rad(lat2) - np.deg2rad(lat1) + dLon = np.deg2rad(lon2) - np.deg2rad(lon1) + a = np.sin(dLat / 2) * np.sin(dLat / 2) + np.cos(np.deg2rad(lat1)) * np.cos( + np.deg2rad(lat2) + ) * np.sin(dLon / 2.0) * np.sin(dLon / 2.0) + c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1.0 - a)) + return R * c + + +def compute_average_resolution(D, xmin, xmax, ymin, ymax): + import numpy as np + + lon = D["XLONG"][0] + lat = D["XLAT"][0] + + dy_avg = np.average( + haversine( + lat[xmin:xmax, ymin : ymax - 1], + lon[xmin:xmax, ymin : ymax - 1], + lat[xmin:xmax, ymin + 1 : ymax], + lon[xmin:xmax, ymin + 1 : ymax], + ) + ) + dx_avg = np.average( + haversine( + lat[xmin : xmax - 1, ymin:ymax], + lon[xmin : xmax - 1, ymin:ymax], + lat[xmin + 1 : xmax, ymin:ymax], + lon[xmin + 1 : xmax, ymin:ymax], + ) + ) + + return (dy_avg + dx_avg) / 2 + + +def datetime_to_float(d): + import datetime + + return d.timestamp() + + +# -------------------------------------------------------------------------------- + + +def plot_var( + filenames, outfile, deg_round_loc=0.2, location_lon=10.8978, location_lat=48.3705 +): + """ + Plot meteogram for given location (location_lon,location_lat) + Script uses all pixels around location in a radius of ''deg_round_loc'' + to attain model spread. + """ + import matplotlib as mpl + + mpl.use("Agg") + import matplotlib.pyplot as plt + import numpy as np + import netCDF4 + import datetime + + from matplotlib.colors import LinearSegmentedColormap + from matplotlib.font_manager import FontProperties + + def draw_vlines(ax, labels=False): + for t in time[1:]: + if t.hour == 0.0 and t.minute == 0.0: + ax.vlines( + datetime_to_float(t), + *ax.get_ylim(), + color="black", + linestyle=":", + lw=1 + ) + for t in time: + if t.hour == 1.0 and t.minute == 0.0: + if labels: + ycoo = ax.get_ylim() + ax.text( + datetime_to_float(t + datetime.timedelta(hours=11)), + -ycoo[1] / 3, + t.strftime("%A %-d %b"), + ha="center", + ) + + # -----------------reading input files:------------------------- + + def get_only_valid_file(f): + try: + nc = netCDF4.Dataset(f) + ts = extract_times(nc, None) + if len(ts) == 0: + nc = None + except: + nc = None + return nc + + DS = [get_only_valid_file(f) for f in sorted(filenames)] + DS = [x for x in DS if not x is None] + + (xmin, xmax), (ymin, ymax), xlon, xlat = corresponding_pixels( + DS[0], location_lon, location_lat, deg_round_loc + ) + + time = np.array([ftimes(D) for D in DS]) + + extent = (xmin, xmax, ymin, ymax) + + # -----------------reading data from input files----------------------- + T2m = [fT2m(D, *extent) for D in DS] + Tavg, Tmin, Tmax = [x(T2m, axis=(1, 2)) for x in [np.average, np.min, np.max]] + + Td2m = [fTd2m(D, *extent) for D in DS] + Tdavg, Tdmin, Tdmax = [x(Td2m, axis=(1, 2)) for x in [np.average, np.min, np.max]] + + wind_speed = [wind10m(D, location_lat, location_lon)[0] for D in DS] + wind_u = [wind_uv(D, location_lat, location_lon)[0] for D in DS] + wind_v = [wind_uv(D, location_lat, location_lon)[1] for D in DS] + mslp = [Mslp(D, location_lat, location_lon) for D in DS] + + drain = fdrain(DS, *extent) + ravg = np.average(drain, axis=(1, 2)) + + dsnow = fdsnow(DS, *extent) + snow = np.average(dsnow, axis=(1, 2)) + + conv = [cin_cape(D, location_lat, location_lon) for D in DS] + cape = [x[1] for x in conv] + cin = [x[0] for x in conv] + + w_vert = [vert_wind(D, location_lat, location_lon) for D in DS] + u_vert = [x[0] for x in w_vert] + v_vert = [x[1] for x in w_vert] + + T_vert = [vert_T(D, location_lat, location_lon) for D in DS] + + cf_vert = [list(fcfrac(D, location_lat, location_lon)) for D in DS] + + # figure settings: + fig = plt.figure(figsize=(8.27 * 1.25, 11.69 * 1.25), num=1) + plt.clf() + + mplotkw = {"ls": "-", "marker": "."} + + ax = plt.subplot2grid((8, 20), (0, 19), rowspan=3) + ax0 = plt.subplot2grid((8, 20), (0, 0), colspan=19, rowspan=3) + ax1 = plt.subplot2grid((8, 20), (3, 0), colspan=19, rowspan=2) + ax2 = plt.subplot2grid((8, 20), (5, 0), colspan=19) + ax3 = plt.subplot2grid((8, 20), (6, 0), colspan=19) + ax4 = plt.subplot2grid((8, 20), (7, 0), colspan=19) + + # for the xticks: + hours = np.array(list(map(lambda x: x.hour, time))) + ind = list(np.where(hours % 6 == 0)[0]) + + # 0.--------------vertical plot------------------------ + + ax0.grid(axis="y") + + # (vertical) cloud fraction + levels = np.arange(400, 1030, 30) + cf_vert = np.matrix(cf_vert).transpose() + cmap = LinearSegmentedColormap.from_list("", ["steelblue", "white"]) + time_mod = list(map(lambda x: datetime_to_float(x), time)) + cs = ax0.contourf( + time_mod, levels, cf_vert, levels=[x / 8 for x in range(10)], cmap=cmap + ) + ax0.contourf(time_mod, levels, cf_vert, levels=[x / 8 for x in range(9)], cmap=cmap) + cbar = fig.colorbar(cs, cax=ax, orientation="vertical") + cbar.set_ticks([x / 8 for x in range(0, 9)]) + try: + cbar.ax.set_yticklabels(["{}/8".format(x) for x in range(0, 9)]) + except: + print("Plotting tick labels failed") + ax.set_ylabel("cloud cover") + + # (vertical) temperature contours + levels = np.arange(400, 1060, 60) + T_vert = np.matrix(T_vert).transpose() + cs = ax0.contour(time_mod, levels, T_vert, cmap="plasma") + ax0.clabel(cs, fmt="%1.0f", fontsize=10) + + # (vertical) wind barbs: + for i, j in enumerate(time_mod[:]): + if i % (len(time) // 24) == 0: + for k, l in enumerate(levels): + ax0.barbs( + j, + l, + u_vert[i][k], + v_vert[i][k], + length=5, + pivot="middle", + linewidth=0.5, + ) + + ax0.set_yticks(np.arange(400, 1050, 50)) + ax0.set_ylim(980, 400) + ax0.set_xlim( + datetime_to_float(time[0] - datetime.timedelta(hours=3 * len(time) / (24 * 3))), + datetime_to_float( + time[-1] + datetime.timedelta(hours=3 * len(time) / (24 * 3)) + ), + ) + ax0.set_ylabel("pressure [hPa]") + try: + ax0.set_xticklabels([], minor=False) + ax0.set_xticklabels([], minor=True) + except: + print("Plotting tick labels failed") + ax0.set_xticks([time_mod[i] for i in ind]) + draw_vlines(ax0) + + # 1.--------------temperature plot------------------------ + + ax1.fill_between(time_mod, Tdmin, Tdmax, alpha=0.3, color="lightblue") + ax1.fill_between(time_mod, Tmin, Tmax, alpha=0.3, color="grey") + ax1.grid(axis="y") + ax1.plot(time_mod, Tdavg, color="blue", label="Td(2m)", **mplotkw) + ax1.plot(time_mod, Tavg, color="black", label="T(2m)", **mplotkw) + + actual_day = time[0].day + actual_month = time[0].month + day = 0 + + # min and max temperature for the current day + for i, j in enumerate(time): + + if time[i - 1].month != actual_month and (i != 0): + day += 1 + actual_month = time[i - 1].month + + elif ( + (day != time[i - 1].day - actual_day) + and (i != 0) + and (time[i - 1].day - actual_day >= 0) + ): + day += 1 + + try: + if Tavg[i] == min(Tavg[day * 24 : (day + 1) * 24]): + ax1.text( + datetime_to_float(time[i]), + Tavg[i] + 2, + str(np.round(Tavg[i], 1)), + color="blue", + ha="center", + ) + elif Tavg[i] == max(Tavg[day * 24 : (day + 1) * 24]): + ax1.text( + datetime_to_float(time[i]), + Tavg[i] + 1.5, + str(np.round(Tavg[i], 1)), + color="red", + ha="center", + ) + except: + pass + + ax1.set_ylabel(r"temp. [$^\circ$C]") + ax1.legend(loc="best", fontsize="small", ncol=3) + ax1.set_ylim(min(min(Tmin), min(Tdmin)) - 3, max(max(Tmax), max(Tdmax)) + 7.5) + ax1.set_xlim( + datetime_to_float(time[0] - datetime.timedelta(hours=3 * len(time) / (24 * 3))), + datetime_to_float( + time[-1] + datetime.timedelta(hours=3 * len(time) / (24 * 3)) + ), + ) + try: + ax1.xaxis.set_ticklabels([], minor=True) + ax1.xaxis.set_ticklabels([], minor=False) + except: + print("Plotting tick labels failed") + ax1.set_xticks([time_mod[i] for i in ind]) + draw_vlines(ax1) + + # 2.--------------10m wind plot------------------------ + + ax2.grid(axis="y") + ax2.plot(time_mod, wind_speed, color="black", **mplotkw) + ax2.set_xlim( + datetime_to_float(time[0] - datetime.timedelta(hours=3 * len(time) / (24 * 3))), + datetime_to_float( + time[-1] + datetime.timedelta(hours=3 * len(time) / (24 * 3)) + ), + ) + for i, j in enumerate(time_mod): + if i % (len(time) // 24) == 0: + ax2.text( + j, + wind_speed[i] - 4.9, + str(int(round(wind_speed[i]))), + color="black", + ha="center", + ) + ax2.barbs( + j, + max(15, max(wind_speed) + 5), + wind_u[i], + wind_v[i], + length=5, + pivot="middle", + linewidth=0.5, + ) + ax2.set_ylabel(r"10 m ws [ms$^{-1}$]") + ax2.set_ylim(-4.9, max(20, max(wind_speed) + 8)) + try: + ax2.xaxis.set_ticklabels([], minor=True) + ax2.xaxis.set_ticklabels([], minor=False) + except: + print("Plotting tick labels failed") + ax2.set_xticks([time_mod[i] for i in ind]) + draw_vlines(ax2) + + # 3.--------------mslp plot------------------------ + + ax3.grid(axis="y") + ax3.plot(time_mod, mslp, color="red", **mplotkw) + ax3.set_ylabel(r"MSLP [hPa]") + try: + ax3.xaxis.set_ticklabels([], minor=True) + ax3.xaxis.set_ticklabels([], minor=False) + except: + print("Plotting tick labels failed") + ax3.set_xlim( + datetime_to_float(time[0] - datetime.timedelta(hours=3 * len(time) / (24 * 3))), + datetime_to_float( + time[-1] + datetime.timedelta(hours=3 * len(time) / (24 * 3)) + ), + ) + ax3.set_ylim(min(mslp) - 2, max(mslp) + 2) + ax3.ticklabel_format(useOffset=False, axis="y") + ax3.set_xticks([time_mod[i] for i in ind]) + draw_vlines(ax3) + + # 4.--------------prec. plot------------------------ + # rain: + + ax4.grid(axis="y") + ax4.bar(time_mod, ravg, color="c", label="tot. prec.", width=3600) + ax4.set_ylabel(r"total prec. [mm/h]") + ax4.set_ylim(0, max(max(ravg), 3) + 1) + ax4.set_xlim( + datetime_to_float(time[0] - datetime.timedelta(hours=3 * len(time) / (24 * 3))), + datetime_to_float( + time[-1] + datetime.timedelta(hours=3 * len(time) / (24 * 3)) + ), + ) + ax4.set_xticks([time_mod[i] for i in ind], minor=True) + ax4.set_yticks( + [x for x in range(0, max(int(max(ravg)), 3) + 2, max(int(max(ravg)), 3) // 3)] + ) + try: + ax4.xaxis.set_ticklabels( + ["00", "06", "12", "18"] * (len(time) // 24), minor=True + ) + except: + print("Plotting tick labels failed") + draw_vlines(ax4, labels=True) + + # print a red (magenta) thunderstorm symbol if cape > 500 (> 1500) Jkg-1 and cin < 100 Jkg-1 + prop = FontProperties() + prop.set_file("/usr/share/fonts/truetype/freefont/FreeMono.ttf") + for i, j in enumerate(time_mod): + if (cin[i] < 100 or np.isnan(cin[i])) and (cape[i] >= 1500): + ax4.text( + datetime_to_float(time[i]), + ax4.get_ylim()[1] * 0.62, + "\u2608", + fontproperties=prop, + color="magenta", + ha="center", + fontsize="large", + ) + elif (cin[i] < 100 or np.isnan(cin[i])) and (cape[i] >= 500): + ax4.text( + datetime_to_float(time[i]), + ax4.get_ylim()[1] * 0.62, + "\u2608", + fontproperties=prop, + color="red", + ha="center", + fontsize="large", + ) + + # snow: + ax5 = ax4.twinx() + ax5.bar(time_mod, snow, color="magenta", label="snow", width=3600) + ax5.set_ylabel(r"snow [cm/h]") + ax5.set_ylim(0, max(max(ravg), 3) + 1) + ax5.set_xlim( + datetime_to_float(time[0] - datetime.timedelta(hours=3 * len(time) / (24 * 3))), + datetime_to_float( + time[-1] + datetime.timedelta(hours=3 * len(time) / (24 * 3)) + ), + ) + ax5.set_yticks( + [ + x + for x in range(max(int(max(ravg)), 3) + 2) + if x % (max(int(max(ravg)), 3) / 3) == 0 + ] + ) + ax5.set_xticks([time_mod[i] for i in ind]) + try: + ax5.yaxis.set_ticklabels([]) + ax5.xaxis.set_ticklabels([]) + except: + print("Plotting tick labels failed") + lines, labels = ax4.get_legend_handles_labels() + lines2, labels2 = ax5.get_legend_handles_labels() + ax4.legend(lines + lines2, labels + labels2, loc="best", fontsize="small", ncol=2) + + # ---------------additional----------------- + + date = datetime.datetime.utcnow().strftime("%H:%M, %B %d, %Y, UTC") + ax0.text( + 1, + 1.01, + "{0:} \n experimental WRF setup \n $\Delta x = {1:.3f}$ km".format( + date, compute_average_resolution(DS[0], *extent) / 1e3 + ), + fontsize="x-small", + transform=ax0.transAxes, + horizontalalignment="right", + ) + + plt.figtext( + 0.483, + 0.948, + r"Christoph Knote (MBEES, Uni Augsburg), Fabian Jakub, Matjaz Puh (Meteorological Institute, LMU Munich)", + fontsize="x-small", + alpha=0.4, + ) + ax0.text( + 0.05, + 1.15, + "MIM Forecast", + transform=ax0.transAxes, + fontsize="xx-large", + horizontalalignment="left", + ) + ax0.text( + 0.05, + 1.05, + " Meteogram @{0:.5f},{1:.5f} \n using {2:} pts for output statistics".format( + location_lon, location_lat, np.size(xlon) + ), + transform=ax0.transAxes, + fontsize="small", + horizontalalignment="left", + ) + + # plt.ticklabel_format(useoffset=False, axis='y') + plt.savefig(outfile, bbox_inches="tight") + + +if __name__ == "__main__": + import datetime + from argparse import ArgumentParser, ArgumentTypeError + + def valid_date(s): + try: + return datetime.datetime.strptime(s, "%Y-%m-%d") + except ValueError: + raise ArgumentTypeError("Not a valid date: '{0}'.".format(s)) + + parser = ArgumentParser() + parser.add_argument( + "wrf_data_fpath_pattern", + type=str, + help="WRF input data path pattern (including strptime pattern for date and {domain:s} for domain (01, 02, ...))", + ) + parser.add_argument( + "start_date", type=valid_date, help="First date to process (YYYY-mm-dd)" + ) + parser.add_argument( + "end_date", type=valid_date, help="Last date to process (YYYY-mm-dd)" + ) + parser.add_argument( + "interval_hours", type=int, help="Time interval in hours for plotting" + ) + parser.add_argument("domain", type=int, help="Which domain to process") + parser.add_argument( + "plot_fpath_pattern", + type=str, + help="Plot output path pattern (including strptime pattern for date and {plot:s} for plot identifier (e.g., 500hPa_europe...)", + ) + parser.add_argument( + "--spatial_patch_deg", + type=float, + default=0.1, + help="Size of the patch (in degrees) to use for uncertainty calculations", + ) + + args = parser.parse_args() + # args = parser.parse_args( [ "/alcc/gpfs2/scratch/mbees/knotechr/archive/WRF/operational_chemistry/wrfout_d{domain:02d}_%Y-%m-%d_%H:%M:%S", "2021-03-11", "2021-03-14", "3", "1", "{plot:s}_%Y%m%d%H.png" ]) + + import numpy as np + + filenames = [ + d.strftime(args.wrf_data_fpath_pattern).format(domain=args.domain) + for d in np.arange( + args.start_date, + args.end_date, + datetime.timedelta(hours=args.interval_hours), + dtype=datetime.datetime, + ) + ] + + try: + plot_var( + filenames, + args.start_date.strftime(args.plot_fpath_pattern.format(plot="met")), + deg_round_loc=args.spatial_patch_deg, + ) + except: + import matplotlib.font_manager as fm + + print("Rebuilding font cache") + fm._rebuild() + plot_var( + filenames, + args.start_date.strftime(args.plot_fpath_pattern.format(plot="met")), + deg_round_loc=args.spatial_patch_deg, + ) diff --git a/blueprints/highres_chemistry/pp/plotting/meteograms/run.sh b/blueprints/highres_chemistry/pp/plotting/meteograms/run.sh new file mode 100755 index 0000000000000000000000000000000000000000..bac917862bea4ad109826ea65f60682b3db868b6 --- /dev/null +++ b/blueprints/highres_chemistry/pp/plotting/meteograms/run.sh @@ -0,0 +1,26 @@ +#!/bin/bash -l +#SBATCH --partition=alcc1 +#SBATCH -o /alcc/gpfs2/scratch/mbees/knotechr/WRF/postprocessing/meteograms.%j.%N.out +#SBATCH -J meteograms +#SBATCH --ntasks=1 +#SBATCH --mem=10G +#SBATCH --mail-type=FAIL +#SBATCH --mail-user=christoph.knote@med.uni-augsburg.de +#SBATCH --time=00:30:00 + +module load gnu geos + +export MPLBACKEND="agg" + +scriptPath=/alcc/gpfs2/home/u/knotechr/wrfotron/blueprints/operational_chemistry/pp/plotting/meteograms/plot.py +wrfDataPathPattern="/alcc/gpfs2/scratch/mbees/knotechr/archive/WRF/operational_chemistry/wrfout_d{domain:02d}_%Y-%m-%d_%H:%M:%S" +firstDateToProcess=`date --date="2 days ago" "+%Y-%m-%d"` +lastDateToProcess=`date --date="3 days" "+%Y-%m-%d"` +intervalHours=1 +domain=2 +plotDir="/alcc/gpfs2/scratch/mbees/knotechr/plots/WRF/operational_chemistry/meteograms" +plotPathPattern="${plotDir}/{plot:s}_%Y%m%d%H.png" + +mkdir -p $plotDir + +python3 $scriptPath "$wrfDataPathPattern" "$firstDateToProcess" "$lastDateToProcess" "$intervalHours" "$domain" "$plotPathPattern" diff --git a/blueprints/highres_chemistry/waccm.inp b/blueprints/highres_chemistry/waccm.inp new file mode 100644 index 0000000000000000000000000000000000000000..4491518846c7d80bf5aaf82466771fa780ef30db --- /dev/null +++ b/blueprints/highres_chemistry/waccm.inp @@ -0,0 +1,67 @@ +&control + +do_ic = __do_ic__, +do_bc = __do_bc__, + +domain = __domain__, + +dir_wrf = './', +dir_moz = './', +fn_moz = 'moz0000.nc', +moz_var_suffix='', + +def_missing_var = .true. + +spc_map = 'o3 -> O3', 'n2o -> N2O', 'no -> NO', +'no2 -> NO2', 'nh3 -> NH3', 'hno3 -> HNO3', 'hno4 -> HO2NO2', +'n2o5 -> N2O5', 'h2o2 -> H2O2', +'ch4 -> CH4', 'co -> CO', 'ch3ooh -> CH3OOH', +'hcho -> CH2O', 'ch3oh -> CH3OH', 'c2h4 -> C2H4', +'ald -> CH3CHO', 'acet ->CH3COCH3', 'mgly -> CH3COCHO', +'pan -> PAN', 'mpan ->MPAN', 'macr -> MACR', +'mvk -> MVK', 'c2h6 -> C2H6', 'c3h6 -> C3H6', 'c3h8 -> C3H8', +'c2h5oh -> C2H5OH', 'apin -> 0.75*MTERP', 'bpin -> 0.25*MTERP', +'isopr -> ISOP','acetol -> HYAC', 'mek -> MEK', +'bigene -> BIGENE', 'bigalk -> BIGALK', +'tol -> TOLUENE', 'benzene -> BENZENE','xylenes ->XYLENES', +'cres -> CRESOL', 'dms -> DMS', 'so2 -> SO2' +! +! aerosols +! +'oc_a01->0.1216*pom_a1+0.9886*soa1_a2+0.1216*soa1_a1+0.9886*soa2_a2+ 0.1216*soa2_a1+0.9886*soa3_a2+0.1216*soa3_a1+0.9886*soa4_a2+0.1216*soa4_a1+ 0.9886*soa5_a2+ 0.1216*soa5_a1;1.e9', +'oc_a02->0.7618*pom_a1+0.0114*soa1_a2+0.7618*soa1_a1+0.0114*soa2_a2+ 0.7618*soa2_a1+0.0114*soa3_a2+0.7618*soa3_a1+0.0114*soa4_a2+0.7618*soa4_a1+ 0.0114*soa5_a2+ 0.7618*soa5_a1;1.e9', +'oc_a03->0.1164*pom_a1+0.0000*soa1_a2+0.1164*soa1_a1+0.0000*soa2_a2+ 0.1164*soa2_a1+0.0000*soa3_a2+0.1164*soa3_a1+0.0000*soa4_a2+0.1164*soa4_a1+ 0.0000*soa5_a2+ 0.1164*soa5_a1;1.e9', +'oc_a04->0.0002*pom_a1+0.0000*soa1_a2+0.0002*soa1_a1+0.0000*soa2_a2+ 0.0002*soa2_a1+0.0000*soa3_a2+0.0002*soa3_a1+0.0000*soa4_a2+0.0002*soa4_a1+ 0.0000*soa5_a2+ 0.0002*soa5_a1;1.e9', +'bc_a01->0.1216*bc_a1+0.1216*bc_a4;1.e9', +'bc_a02->0.7618*bc_a1+0.7618*bc_a4;1.e9', +'bc_a03->0.1164*bc_a1+0.1164*bc_a4;1.e9', +'bc_a04->0.0002*bc_a1+0.0002*bc_a4;1.e9', +'so4_a01->0.9886*so4_a2+0.1216*so4_a1+0.0000*so4_a3;1.e9', +'so4_a02->0.0114*so4_a2+0.7618*so4_a1+0.0000*so4_a3;1.e9', +'so4_a03->0.0000*so4_a2+0.1164*so4_a1+0.0995*so4_a3;1.e9', +'so4_a04->0.0000*so4_a2+0.0002*so4_a1+0.9003*so4_a3;1.e9', +'nh4_a01->0.1856*so4_a2+0.0050*so4_a1+0.0000*so4_a3;1.e9', +'nh4_a02->0.0021*so4_a2+0.0930*so4_a1+0.0000*so4_a3;1.e9', +'nh4_a03->0.0000*so4_a2+0.0203*so4_a1+0.0186*so4_a3;1.e9', +'nh4_a04->0.0000*so4_a2+0.0000*so4_a1+0.1690*so4_a3;1.e9', +'no3_a01->0.0000*so4_a2+0.0000*so4_a1+0.0000*so4_a3;1.e9', +'no3_a02->0.0000*so4_a2+0.0000*so4_a1+0.0000*so4_a3;1.e9', +'no3_a03->0.0000*so4_a2+0.0000*so4_a1+0.0000*so4_a3;1.e9', +'no3_a04->0.0000*so4_a2+0.0000*so4_a1+0.0000*so4_a3;1.e9', +'na_a01->0.3889*ncl_a2+0.0479*ncl_a1+0.0000*ncl_a3;1.e9', +'na_a02->0.0045*ncl_a2+0.2997*ncl_a1+0.0000*ncl_a3;1.e9', +'na_a03->0.0000*ncl_a2+0.0458*ncl_a1+0.0391*ncl_a3;1.e9', +'na_a04->0.0000*ncl_a2+0.0000*ncl_a1+0.3542*ncl_a3;1.e9', +'cl_a01->0.5996*ncl_a2+0.0737*ncl_a1+0.0000*ncl_a3;1.e9', +'cl_a02->0.0068*ncl_a2+0.4621*ncl_a1+0.0000*ncl_a3;1.e9', +'cl_a03->0.0000*ncl_a2+0.0709*ncl_a1+0.0604*ncl_a3;1.e9', +'cl_a04->0.0000*ncl_a2+0.0001*ncl_a1+0.5462*ncl_a3;1.e9', +'oin_a01->0.9886*dst_a2+0.1216*dst_a1+0.0000*dst_a3;1.e9', +'oin_a02->0.0114*dst_a2+0.7618*dst_a1+0.0002*dst_a3;1.e9', +'oin_a03->0.0000*dst_a2+0.1164*dst_a1+0.0995*dst_a3;1.e9', +'oin_a04->0.0000*dst_a2+0.0002*dst_a1+0.9003*dst_a3;1.e9', +'num_a01->0.9996*num_a2+0.7135*num_a1+0.0000*num_a3;1.0', +'num_a02->0.0004*num_a2+0.2470*num_a1+0.2118*num_a3;1.0', +'num_a03->0.0000*num_a2+0.0016*num_a1+0.6258*num_a3;1.0', +'num_a04->0.0000*num_a2+0.0000*num_a1+0.3501*num_a3;1.0', +/ diff --git a/blueprints/highres_chemistry/wesely.inp b/blueprints/highres_chemistry/wesely.inp new file mode 100644 index 0000000000000000000000000000000000000000..19febe530505bd8c7a3660f8b35bddbd85f1341b --- /dev/null +++ b/blueprints/highres_chemistry/wesely.inp @@ -0,0 +1,6 @@ +&control + +wrf_dir = './' +domains = __domains__, + +/ diff --git a/tools/plot_domains.py b/tools/plot_domains.py index 0b8847174e6e69df0d6d527d8f5432c92c66305f..da2df98fbdcb8deac9e1a1974bb6503a105a60a9 100644 --- a/tools/plot_domains.py +++ b/tools/plot_domains.py @@ -844,11 +844,12 @@ if __name__ == '__main__': print( "Try e_sn = {:d}".format( ceil(wps_info.e_sn_original[i] / wps_info.parent_grid_ratio[i]) * wps_info.parent_grid_ratio[i] + 1) ) print( "" ) - #ll_x, ll_y = wpsproj.transform_point(5.5, 47., latlonproj) - #ur_x, ur_y = wpsproj.transform_point(15., 55., latlonproj) + ll_x, ll_y = wpsproj.transform_point(5.5, 47., latlonproj) + ur_x, ur_y = wpsproj.transform_point(15., 55., latlonproj) # - #ax1.set_xlim([ll_x, ur_x]) - #ax1.set_ylim([ll_y, ur_y]) + ax1.set_xlim([ll_x, ur_x]) + ax1.set_ylim([ll_y, ur_y]) + # Augsburg lon, lat = wpsproj.transform_point(10.898333, 48.371667, latlonproj) ax1.scatter(lon, lat, marker='o', color='r', label='point')