import os
import errno # used in silent_remove
import sys
from datetime import datetime, timedelta
import shutil
import subprocess
import glob
# Logging functions in order of increasing severity:
# debug, info, warning, error, critical
# A local logger instance is used in run(...) for development purposes, but
# it could also be initialised after the import statement. More info:
# *) https://docs.python.org/3.7/library/logging.html
# *) https://docs.python.org/3.7/howto/logging.html
# *) https://stackoverflow.com/a/48787602
import logging
# JvP 20210202: defined MODULE_NAME and module logger
MODULE_NAME = __name__[__name__.index('TM5'):] if 'TM5' in __name__ else __name__
logger = logging.getLogger(MODULE_NAME)
# JvP 20210723: added import of netCDF and numpy
import netCDF4 as nc4
import numpy as np
[docs]
def read_point_output(filename):
"""
PURPOSE
Read a TM5 point_output.nc4 file.
IN/OUT
filename = string with the filename to read.
VERSION CHANGE HISTORY
1.0 23-07-2021 by J.C.A. van Peet.
Original code
"""
# Name of the program
PROG_NAME = MODULE_NAME+".read_point_output"
# Local logger, see function run below for more information.
logger = logging.getLogger(PROG_NAME)
logger.setLevel(logging.DEBUG) # Lowest level possible
# ************************************
# * NO MORE SETTINGS BELOW THIS LINE *
# ************************************
# Start message
logger.debug("")
logger.debug(PROG_NAME+" => start")
# Open file
nc_id = nc4.Dataset( filename, 'r' )
nc_id.set_auto_maskandscale(False)
# Read datasets
mixing_ratio = nc_id['glb600x400/mixing_ratio' ][()]
mixing_ratio_sigma_spat = nc_id['glb600x400/mixing_ratio_sigma_spat'][()]
mixing_ratio_sigma_temp = nc_id['glb600x400/mixing_ratio_sigma_temp'][()]
nsamples = nc_id['glb600x400/nsamples' ][()]
sample_id = nc_id['glb600x400/id' ][()]-1 # Convert to 0-based indices
sampling_strategy = nc_id['glb600x400/sampling_strategy' ][()]
station_id = nc_id['glb600x400/station_id' ][()] # Not converted to 0-based indices...
# Close file
nc_id.close()
# Stop message
logger.debug(PROG_NAME+" => stop")
logger.debug("")
# Return data
return {'mixing_ratio' : mixing_ratio,
'mixing_ratio_sigma_spat' : mixing_ratio_sigma_spat,
'mixing_ratio_sigma_temp' : mixing_ratio_sigma_temp,
'nsamples' : nsamples,
'sample_id' : sample_id,
'sampling_strategy' : sampling_strategy,
'station_id' : station_id,
}
# end function read_point_output
# end function read_point_input
[docs]
def calc_and_write_point_dep(pod, apod, apid, filename, adj_out=None):
"""
PURPOSE
Calculate the point departure data and write it to file.
In TM5, this is calculated in the "dep" job by TMVar.py, which calls
Observations_JRC.PointObs_JRC.applyPointObs(...) as:
point.applyPointObs( trackfile=trackfile, trackfile_apri=trackfile_apri,
mismatchfile=departurefile,
Params=None, montecarlo=False )
where:
trackfile = some/path/to/point_output.nc4 (written by TM5)
trackile_apri = some/path/to/point_output.nc4 from the apri run (i.e. iter 1)
departurefile = some/path/to/point_departures.nc4 (which is read by the
adjoint run of TM5, and which is set in params.py when
the config file is written).
IN/OUT
pod = Point Output Data, the dict produced by read_point_output() above
apod = A-priori Point Output Data, the dict produced by read_point_output() above
apid = A-priori Point Input Data, the dict produced by read_point_input() above
filename = Name of file to write
KWARGS
adj_out = Forcings calculated by CIF as "R^(-1) x (H(xb) - y)". If present,
replace forcings as TM5 would have calculated them with these
values.
VERSION CHANGE HISTORY
2.0 28-10-2021 by J.C.A. van Peet.
Added the adj_out keyword, which is written to file as 'forcing' while
the original TM5 focings are retained in the 'forcing_TM5' variable for
comparison purposes.
1.0 23-07-2021 by J.C.A. van Peet.
Original code, based on Observations_JRC.PointObs_JRC.applyPointObs(...)
"""
# Name of the program
PROG_NAME = MODULE_NAME+".calc_and_write_point_dep"
# Local logger, see function run below for more information.
logger = logging.getLogger(PROG_NAME)
logger.setLevel(logging.DEBUG) # Lowest level possible
# ************************************
# * NO MORE SETTINGS BELOW THIS LINE *
# ************************************
# Start message
logger.debug("")
logger.debug(PROG_NAME+" => start")
# In TM5, you can have multiple tracers, although only CH4 is used at the
# moment. For compatibility, set i_tracer to 0.
i_tracer = 0
# Count:
nsamples = pod['nsamples'].size
# Extract observation operator variables, observations, and obs. err. std.dev.
# for the observation id's in this group:
date_components = apid['date_components' ]
lat = apid['lat' ]
lon = apid['lon' ]
alt = apid['alt' ]
obs_mixing = apid['mixing_ratio' ]
obs_error = apid['mixing_ratio_error']
# also extract non-standard variables originating from point input:
obs_stdv = apid['mixing_ratio_stdv']
tower = apid['tower' ]
bias_id = apid['bias_id' ]
# Extract mixing ratio simulated by model for selected tracer:
mod_mixing = pod['mixing_ratio'][:, i_tracer]
# Calculate model error; extract this from the apriori model run since it
# depends on the (per iteration) changing concentrations.
#
# individual spatial and temporal contributions:
mod_error_spat = apod['mixing_ratio_sigma_spat'][:, i_tracer]
mod_error_temp = apod['mixing_ratio_sigma_temp'][:, i_tracer]
#
# combined including temporal std.dev. in observations:
mod_error = np.sqrt( mod_error_spat**2 + np.maximum(mod_error_temp, obs_stdv)**2 )
# init array with bias values, set to zero:
bias = np.zeros( (nsamples), float )
# Compute and write error: R + HPH'
total_error = np.sqrt( obs_error**2 + mod_error**2 )
# Calculate mismatch : (Hx+b) - y - pert_y
mismatch = mod_mixing + bias - obs_mixing
# Forcing : (HPH'+R)**-1 ( (Hx+b) - y )
forcing = mismatch / total_error**2
# >>> Write to file
# Delete any pre-existing files, and open the output file. Usually,
# selecting "w" when opening a file should create a new file overwriting
# any existing file. But if the file is opened in hdfview (which happens
# regularly), this issues an error. Deleting the file explicitly supresses
# that error. But os.remove raises an error if the file does not exist
# (at least in some Python versions I've tried this), so wrap it in a try
# statement to suppress it.
try:
os.remove(filename)
except:
pass
# end try
# Open file
nc_id = nc4.Dataset( filename, 'w' )
# Global dimensions
nc_id.createDimension('tracer', 1)
nc_id.createDimension('cdate', 6)
# Create group
grp_id = nc_id.createGroup('glb600x400')
# Group dimensions
grp_id.createDimension( 'samples', nsamples )
# Copy observatoin id's
var = grp_id.createVariable('id', 'i', ('samples',))
var[()] = apod['sample_id'] + 1 # Convert back to 1-based indices
# Copy station id's
var = grp_id.createVariable('station_id', 'i', ('samples',))
var[()] = apod['station_id']
# Write latitudes
var = grp_id.createVariable('lat', 'd', ('samples',))
var[()] = lat
# Write longitudes
var = grp_id.createVariable('lon', 'd', ('samples',))
var[()] = lon
# Write altitudes
var = grp_id.createVariable('alt', 'd', ('samples',))
var[()] = alt
# Write tower flags
var = grp_id.createVariable( 'tower', 'i', ('samples',) )
var[()] = tower
# Write date components
var = grp_id.createVariable('date_components', 'i', ('samples', 'cdate'))
var[()] = np.int32( date_components )
# Write bias type
var = grp_id.createVariable( 'bias_id', 'i', ('samples',) )
var[()] = bias_id
# Write model bias:
var = grp_id.createVariable( 'bias', 'd', ('samples', 'tracer') )
var[:, i_tracer] = bias[:]
# Write number of samples used in temporal sampling strategy
var = grp_id.createVariable('nsamples', 'i', ('samples',))
var[()] = pod['nsamples']
# Write sampling strategy codes
var = grp_id.createVariable('sampling_strategy', 'i', ('samples',))
var[()] = pod['sampling_strategy']
# Write model error std.dev.
var = grp_id.createVariable('model_error', 'd', ('samples', 'tracer'))
var[:, i_tracer] = mod_error[:]
# Write observation error
var = grp_id.createVariable('obs_error', 'd', ('samples', 'tracer'))
var[:, i_tracer] = obs_error[:]
# Write total error: R + HPH'
var = grp_id.createVariable('error', 'd', ('samples', 'tracer'))
var[:, i_tracer] = total_error[:]
# Write mismatch : (Hx+b) - y - pert_y
var = grp_id.createVariable('mismatch', 'd', ('samples', 'tracer'))
var[:, i_tracer] = mismatch[:]
# JvP 20211028
# Replaced the TM5 forcings with the CIF adj_out (if present).
# As far as I understand the TM5 code, the netCDF file is read by the file
# called 'user_output_flask_adj.F90', and only the variables 'samples',
# 'nsamples', 'sampling_strategy', 'forcing', 'date_components', 'alt',
# 'tower', 'lat', and 'lon'.
if( adj_out is None ):
# Original TM5 forcing
# Write TM5 forcing: (HPH'+R)**-1 ( (Hx+b) - y )
var = grp_id.createVariable('forcing', 'd', ('samples', 'tracer'))
var[:, i_tracer] = forcing[:]
else:
# First write TM5 forcing: (HPH'+R)**-1 ( (Hx+b) - y )
# Note the '_TM5' suffix in the variable name. The variable 'forcing'
# will contain the CIF 'adj_out'.
var = grp_id.createVariable('forcing_TM5', 'd', ('samples', 'tracer'))
var[:, i_tracer] = forcing[:]
var.note = "Original TM5 forcing: (HPH'+R)**-1 ( (Hx+b) - y )"
# Now write the observation increments as calculated by CIF:
# R^(-1) x (H(xb) - y)
var = grp_id.createVariable('forcing', 'd', ('samples', 'tracer'))
var[:, i_tracer] = adj_out.to_numpy()[:]
var.note = "CIF observation increments: R^(-1) x (H(xb) - y)"
# For debugging...
#print(__file__)
#import code
#code.interact(local=dict(locals(), **globals()))
#raise RuntimeError
# end if
# Close file
nc_id.close()
# <<< Write to file
# Stop message
logger.debug(PROG_NAME+" => stop")
logger.debug("")
# Return data.
#return
# end function calc_and_write_point_dep
