Commit 72f67d46 authored by Espen Sollum's avatar Espen Sollum
Browse files

Added background plugin for FLEXPART

parent dbcb9c1b
......@@ -33,22 +33,8 @@ def create_domain(domain,
lon = lonc + dlon/2.
lat = latc + dlat/2.
# Meshgrids
zlon, zlat = np.meshgrid(lon, lat)
lonc = np.append(lonc, lonc[-1] + dlon)
latc = np.append(latc, latc[-1] + dlat)
zlonc, zlatc = np.meshgrid(lonc, latc)
domain.zlon = zlon
domain.zlat = zlat
domain.zlonc = zlonc
domain.zlatc = zlatc
domain.lon = lon
domain.lat = lat
# Adding global grid definition
# Global grid definition
xmin_glob = domain.xmin_glob
ymin_glob = domain.ymin_glob
nlon_glob = domain.nlon_glob
......@@ -64,6 +50,31 @@ def create_domain(domain,
lon_glob = lonc_glob + dlon_glob/2.
lat_glob = latc_glob + dlat_glob/2.
# Meshgrids
zlon, zlat = np.meshgrid(lon, lat)
lonc = np.append(lonc, lonc[-1] + dlon)
latc = np.append(latc, latc[-1] + dlat)
zlonc, zlatc = np.meshgrid(lonc, latc)
zlon_glob, zlat_glob = np.meshgrid(lon_glob, lat_glob)
domain.zlon = zlon
domain.zlat = zlat
domain.zlonc = zlonc
domain.zlatc = zlatc
domain.lon = lon
domain.lat = lat
domain.lon_glob = lon_glob
domain.lat_glob = lat_glob
domain.lonc_glob = lonc_glob
domain.latc_glob = latc_glob
domain.zlon_glob = zlon_glob
domain.zlat_glob = zlat_glob
# Indices for inversion domain into global domain (relative to lower left corner)
ix1 = np.argmin(np.abs(lonc_glob - domain.lon[0] - 1))
iy1 = np.argmin(np.abs(latc_glob - domain.lat[0] - 1))
......
......@@ -42,7 +42,6 @@ def read(self, name, tracdir, tracfic, dates,
trcr_flx = []
times = []
# for dd, fic_flx in zip(dates, list_fic_flx):
trcr_flx = []
for file_flx in list(list_fic_flx):
data, lat, lon, time_jd = readnc(os.path.join(tracdir, file_flx),
[self.varname_flx, self.latname_flx,
......
......@@ -10,8 +10,9 @@ from .outputs2native import outputs2native
from .run import run
requirements = {'domain': {'name': 'FLEXPART', 'version': 'std', 'empty': False},
'fluxes': {'name': 'FLEXPART', 'version': 'nc',
'empty': True, 'any': False}}
'fluxes': {'name': 'FLEXPART', 'version': 'nc', 'empty': True, 'any': False}
# 'background': {'name': 'FLEXPART', 'version': 'nc', 'empty': True, 'any': False}
}
def ini_data(plugin, **kwargs):
"""Initializes FLEXPART
......
......@@ -9,8 +9,13 @@ def ini_periods(self, **kwargs):
datef = self.datef
# List of sub-simulation windows
#self.subsimu_dates = \
# date_range(datei, datef, period=self.periods)
# TODO: implement sub-periods? For now just keep the full period
self.subsimu_dates = \
date_range(datei, datef, period=self.periods)
date_range(datei, datef, period='')
# Time steps per day
freq = '{:.0f}s'.format(3600)
......
......@@ -8,7 +8,6 @@
! Compile this by
! f2py -c mod_flexpart.f90 -m mod_flexpart
!
! TODO: add makefile/automatic compile on setup
!
!====================================================================
......@@ -570,7 +569,7 @@ end subroutine read_factor
!!
! --------------------------------------------------
subroutine read_init(filename, gridinit, nxgrid,nygrid,nzgrid)
subroutine read_init(gridinit, filename, nxgrid, nygrid, nzgrid)
! use mod_var
......@@ -578,7 +577,7 @@ subroutine read_init(filename, gridinit, nxgrid,nygrid,nzgrid)
character(len=256), intent(in) :: filename
integer, intent(in) :: nxgrid, nygrid, nzgrid
real, dimension(nxgrid,nygrid,nzgrid), intent(in out) :: gridinit
real, dimension(nxgrid,nygrid,nzgrid), intent(out) :: gridinit
real, parameter :: scaleconc=1.e12
real, parameter :: smallnum=1.e-38
......
......@@ -87,17 +87,35 @@ def read_flexpart_grid(subdir, file_name, fp_header, **kwargs):
# Convert grid times to datetime format
gtime_dt = []
for i in range(len(gtime)):
# gtime_dt[i] = flexpart_header.Flexpartheader.j2d(gtime[i])
if gtime[i] == 0.:
break
# gtime_dt.append(fp_header.j2d(gtime[i]))
gtime_dt.append(j2d(gtime[i]))
return grid_fp, gtime_dt, ngrid
def read_init(subdir, file_name, nx, ny, nz, **kwargs):
""" Reads FLEXPART grid_initial file.
convert from ppt to fraction of one
gridini = gridini*1.e-12
Inputs:
subdir - directory where files are read
file_name - full path name to file name
Returns:
grid_init: Array (nlon, nlat, nheight)
"""
grid_init = mod_flexpart.read_init(os.path.join(subdir, file_name), nx, ny, nz)*1.e-12
return grid_init
def get_spec(subdir, **kwargs):
""" Get species name from simulation output
......
......@@ -19,7 +19,7 @@ def obsoper(self, inputs, mode,
This function maps information from the control space to the observation
space
This version of 'obsoper' was developed for use with FLEXPART
This version of 'obsoper' was developed for use with FLEXPART
backward runs.
Gets concentrations/mixing ratios from FLEXPART flux sensitivities.
......@@ -67,6 +67,7 @@ def obsoper(self, inputs, mode,
controlvect = self.controlvect
obsvect = self.obsvect
subsimu_dates = model.subsimu_dates
dump_debug = obsvect.dump_debug
if mode == 'fwd':
......@@ -89,7 +90,7 @@ def obsoper(self, inputs, mode,
os.path.join(model.run_dir_nest, obsvect.datastore.head(1)['station'][0].upper(),
subdir, 'header_nest'))
# Get the domain definition
species = getattr(controlvect.components, 'fluxes').parameters.attributes[0]
tracer = getattr(getattr(controlvect.components, 'fluxes').parameters, species)
......@@ -109,15 +110,14 @@ def obsoper(self, inputs, mode,
iy2 = model.domain.iy2
else:
raise Exception("For FLEXPART, only hresol:regions are implemented in controlvect")
# Loop through model periods and read model output
self.missingperiod = False
for di, df in zip(subsimu_dates[:-1], subsimu_dates[1:]):
subdir = di.strftime("%Y%m")
# Save to datastore for debugging purposes
obs_ghg = np.empty(len(obsvect.datastore)); obs_ghg[:]=np.nan
obs_bkg = np.empty(len(obsvect.datastore)); obs_bkg[:]=np.nan
......@@ -131,6 +131,9 @@ def obsoper(self, inputs, mode,
print("di, df", di, df)
for obs_i, row in enumerate(obsvect.datastore.itertuples()):
# Subdirectory for FLEXPART footprints
subdir = row.Index.strftime("%Y%m")
# For debugging
obs_check[obs_i] = obs_i
......@@ -166,12 +169,13 @@ def obsoper(self, inputs, mode,
grid_glob *= model.coeff*model.mmair/model.molarmass
# Array for global transport backhground
# Background contribution from fluxes outside domain
hbkg = np.sum(grid_glob[:, :, 0:ngrid-1], axis=2)
hbkg[ix1:ix2, iy1:iy2] = 0.0
# Index to state vector
ind = np.argmin(np.abs(tracer.dates[0::-1] - gtime_glob[0]))
obs_bkg[obs_i] = np.sum(hbkg[:,:].T*controlvect.flx_bkg[ind,:,:])
obs_bkgerr[obs_i] = obs_bkg[obs_i]*tracer.err
......@@ -251,16 +255,17 @@ def obsoper(self, inputs, mode,
obsvect.dx = grad_o
# Add the different components to datastore
obsvect.datastore['obs_bkgerr'] = obs_bkgerr
obsvect.datastore['sim'] = obs_sim
obsvect.datastore['obs_ghg'] = obs_ghg
obsvect.datastore['obs_bkg'] = obs_bkg
obsvect.datastore['obs_model'] = obs_model
obsvect.datastore['obs_sim'] = obs_sim
obsvect.datastore['obs_check'] = obs_check
obsvect.datastore['obs_bkgerr'] = obs_bkgerr
obsvect.datastore['obs_err'] = np.sqrt(obsvect.datastore['obs_bkgerr']**2 +
obsvect.datastore['obserror']**2)
if dump_debug:
obsvect.datastore['obs_ghg'] = obs_ghg
obsvect.datastore['obs_bkg'] = obs_bkg
obsvect.datastore['obs_model'] = obs_model
obsvect.datastore['obs_sim'] = obs_sim
obsvect.datastore['obs_check'] = obs_check
# Save grad_o for inspection
rundir = "{}/obsoperator".format(workdir)
......@@ -320,11 +325,15 @@ def obsoper(self, inputs, mode,
rundir = "{}/obsoperator/{}".format(workdir, mode)
path.init_dir(rundir)
dump_type = obsvect.dump_type
dump_datastore(obsvect.datastore,
file_monit='{}/monitor_{}_.{}'.format(rundir, run_id, dump_type),
mode='w', dump_type=dump_type,
col2dump=['obs_ghg', 'obs_bkg', 'obs_model', 'obs_sim', 'obs_check',
'obs_bkgerr', 'obs_err', 'obs_hx'])
import pdb; pdb.set_trace()
if dump_debug:
dump_datastore(obsvect.datastore,
file_monit='{}/monitor_{}_.{}'.format(rundir, run_id, dump_type),
mode='w', dump_type=dump_type,
col2dump=['obs_ghg', 'obs_bkg', 'obs_model', 'obs_sim', 'obs_check',
'obs_bkgerr', 'obs_err', 'obs_hx'])
# Returning obsvect to the simulator
if mode is 'fwd':
......
......@@ -4,6 +4,7 @@ from .init_rinvprod import init_rinvprod
from .native2obsvect import native2obsvect
from .obsvect2native import obsvect2native
from .rinvprod import rinvprod
from .init_background import init_background
from pycif.plugins.obsvects.standard.utils.check_monitor import check_monitor
from pycif.utils import path
......@@ -15,6 +16,7 @@ import numpy as np
# to initialize the observation vector
requirements = {'measurements': {'any': True, 'empty': True},
'model': {'any': True, 'empty': False}}
# 'background': {'any': True, 'empty': False}}
def ini_data(plugin, **kwargs):
......@@ -57,6 +59,12 @@ def ini_data(plugin, **kwargs):
.format(plugin.workdir, os.path.basename(obs_file))
path.link(obs_file, target)
plugin.has_satellites = False
if np.where((plugin.datastore['level'] < 0))[0].size > 0:
plugin.has_satellites = True
# Initialize background if present
if hasattr(plugin, 'background'):
init_background(plugin, **kwargs)
......@@ -5,4 +5,7 @@ def rinvprod(obsvect, dy):
# At the moment, this is valid for diagonal observation matrices only
# TODO: Include these lines into rinvprod.py, with option for
# non-diagonal matrices, eventually
# ESO TODO: for FLEXPART we would add bkgerr to obserror
return old_div(dy, obsvect.datastore['obserror'] ** 2)
......@@ -61,7 +61,8 @@ def simul(self, chi, grad=True, run_id=-1, **kwargs):
# TODO: Include these lines into rinvprod.py, with option for
# non-diagonal matrices, eventually
departures = obsvect.datastore['sim'] - obsvect.datastore['obs']
j_o = 0.5 * (departures ** 2 / (obsvect.datastore['obserror']**2 + obsvect.datastore['obs_bkgerr']** 2) ).sum()
# j_o = 0.5 * (departures ** 2 / (obsvect.datastore['obserror']**2 + obsvect.datastore['obs_bkgerr']** 2) ).sum()
j_o = 0.5 * (departures ** 2 / obsvect.datastore['obs_err'] ** 2).sum()
zcost = j_b + j_o
......
from types import MethodType
from .baseclass import Plugin
from pycif.utils.check.errclass import PluginError
class Background(Plugin):
@classmethod
def register_plugin(cls, name, version, module, **kwargs):
"""Register a module for a plugin and version with possibly options
Args:
name (str): name of the plugin
version (str): version of the plugin
module (types.ModuleType): module defining the interface
between pyCIF and the plugin
plugin_type (str): type of plugin
**kwargs (dictionary): default options for module
"""
super(Background, cls).register_plugin(name, version, module,
plugin_type='background')
def read_conc(self, name, dir_initconc, file_initconc):
"""Get concentrations from global concentration fields and load them into pyCIF
variables
Args:
self: the model Plugin
name: the name of the component
dir_initconc, file_initconc: flux directory and file format
interpol_conc (bool): if True, interpolates concentrations at time t from
values of surrounding available files
"""
raise PluginError('The function read was not defined')
def calc_init(self, name, dir_init, dates):
""" Calculate initial concentration sensitivities
Args:
self: the model Plugin
name: the name of the component
dir_init: directory for init files
dates: list of dates to extract
"""
raise PluginError('The function read was not defined')
def initiate_template(self):
module = super(Background, self).initiate(plg_type='background')
# Replacing auxiliary functions:
if hasattr(module, 'write'):
self.write = MethodType(module.write, self)
if hasattr(module, 'read_conc'):
self.read_conc = MethodType(module.read_conc, self)
if hasattr(module, 'calc_init'):
self.calc_init = MethodType(module.calc_init, self)
......@@ -19,6 +19,7 @@ class Plugin(object):
plugin_types = {'mode': ['.modes', 'Mode'],
'meteo': ['.meteos', 'Meteo'],
'model': ['.models', 'Model'],
'background' : ['.background', 'Background'],
'fields': ['.fields', 'Fields'],
'fluxes': ['.fluxes', 'Fluxes'],
'domain': ['.domains', 'Domain'],
......
......@@ -45,7 +45,9 @@ def date_range(datei, datef, period='1000000H', subperiod=''):
# If the period is empty, keeps the full period as a sub-period
if period == '':
list_dates = [datei]
# list_dates = [datei]
# ESO: return list of start and end date
list_dates = [datei, datef]
return np.array(list_dates)
# Otherwise, split simulation window along base periods
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment