EnviFlux Vn1.0: Technical and User Guide

EnviFlux Vn1.0: Technical and User Guide November 14, 2025 Ross Bannister, National Centre for Earth Observation, University of Reading, Reading, UK [email protected] 1 Introduction EnviFlux (Environmental Flux inversion tool) is a global 4DVar system to infer surface uxes of trace gases in the Earth's atmosphere from measurements of the gases within a time window. This document describes general, technical, and operational aspects of the system including:  How the problem is represented. This includes the structure of the state vector, the underlying grid, and the processes represented.  The formulation of the forward model (a transport/source problem). This includes the way that the tracer transport equations are solved.  The inverse problem (details of the 4DVar algorithm). This includes the structure of the cost function, the derivation of its gradient with respect to the state vector, and the changes of variables introduced to make this ecient.  How to run the various master programs. This includes:  Generating wind elds at the suitable resolution which will drive the dynamics.  Generating suitable initial conditions and surface uxes for the purposes of testing the forward model and generating data for a nature run.  Running the forward model.  Running back-trajectories.  Generating a-priori state vectors.  Generating synthetic observations.  Running the full 4DVar system to infer surface uxes. The code is designed to be relatively light and ecient, and can be run in a Linux environment (from a laptop to a supercomputer). No parallelisation is exploited. 1 2 How the forward problem is represented  There are four kinds of time-step in EnviFlux:  Minor time-step: this is δt (or dt in the code), the time-step used by the Runge-Kutta solver in the back-trajectory part of the semi-Lagrangian advection scheme.  Wind le time-step (or major time-step): this is ∆t (or Dt in the code), the time separation between wind les. It is also the semi-Lagrangian time step (how often the semi-Lagrangian scheme is run). ∆t must be a multiple of δt .  Diusion time-step: this is τκ (or kappa _ dt in the code), the time-step of the diusion scheme. ∆t must be a multiple of τκ .  Source time-step: this is τρ . This represents how frequently the source eld is updated in the inverse problem.  Let γ= ∆t/δt and is the number of minor steps per major step.  Let there be γ(T+ 1) minor time levels and T+ 1 major time levels.  Let the augmented state vector be x , comprising χ(0) (the initial conditions ( t= 0 ) for the tracer) and ρ(t) (the surface ux at major time t ): x=χ(0) ρ, where ρ=     ρ(0) ρ(τρ) . . . ρ(Tτρ)      . (1) In EnviFlux, χ(0) has units of ppb, and ρ has units of µgm−2s−1 .  The model evolution between two major time steps is done in three stages.  The rst stage is the advection stage (performed by semi-Lagrangian ( SL ) over γ minor timesteps, labelled with k ). χt+1 (k∆t) = SLk→k+1χt(k∆t), where the argument represents the major time-step and the super/subscript represent the minor time-steps.  The second stage is the diusion stage. This solves the following diusion equation over pseudo time τ : ∂χt+1(τ) ∂τ =κh∇2 hχt+1(τ) + κv ∂2χt+1(τ) ∂z2, (2) where κh is the horizontal diusion coecient and κv is the vertical diusion coecient.  The third stage is the source/sink stage. χ([k+ 1]∆t) = Mk→k+1χk(k∆t) + Sρ(I(k∆t/τρ)) ∆t, (3) where I(x) is the integer operator, used here to select the relevant ux eld, S is the `surface operator', which converts the 2D source eld to the surface level of a 3D eld and divides by ∆z0d0 at the surface layer (layer thickness ∆z0 , surface air density d0 ), and Mk→k+1 is the combined eect of the sequence of SLk→k+1 operations and diusion operations over a complete major time-step. For compactness, the full evolution is written as χ([k+ 1]∆t) = Mx k→k+1χ(k∆t), (4) where the x in the superscript indicates that the state vector itself (namely the ux part) contributes to the tracer evolution. 2 The combined advection/diusion/ux equation may be written as ∂χ(t) ∂t =−u(t)· ∇χ(t) + κh∇2 hχ(t) + kv ∂2χ(t) ∂z2+SρI(t/τρ), (5) where κ is the diusion coecient (actually in the developments below, this has dierent values for horizontal and vertical diusions). 3 The inverse problem 3.1 The cost function EnviFlux works by minimising the following cost function: Jδx(δx) = 1 2δx⊺B−1δx+ 1 2 NT X k=0 yk−HkMxb+δx 0→khχb(0) + δχi⊺ R−1 kyk−HkMxb+δx 0→khχb(0) + δχi, (6) where δx is the proposed analysis increment with tracer component δχ (the actual analysis increment is the one that minimises Jδx(δx) ), B is the background error covariance matrix, k is the major timestep index, Hk is the observation operator for time-step k , Mxb+δx 0→k is the compact form of the evolution operator (see Eq. (4)), and Rk is the observation error covariance matrix for time-step k . The observation error covariance matrix is set-up for uncorrelated observation errors at present. δx is an increment and is related to the full state by x=xb+δx , where xb is the background state. The analysis, xa , corresponds to the special value of δx that minimises Jδx(δx) (call this special value δxa ). The analysis is then xa=xb+δxa . Even though the notation in Eq. (6) implies that observations are at major timesteps only, EnviFlux allows observations to be made at intermediate times, which are treated with linear interpolation (see Sect. 3.2). The gradient of the cost function is ∇δxJ(δx) = B−1δx+ NT X k=0 Mxb+δx 0→k⊺ R−1 kyk−HkMxb+δx 0→khχb(0) + δχi, (7) which is used by EnviFlux's conjugate gradient-based descent algorithm ([ ? , ? ]) to minimise the cost function. A version of the gradient calculation that is more ecient to compute than (7) is used in the code, which is documented in Sect. 3.2. 3.2 Intermediate observation times and an ecient gradient computation The cost function, Eq. (6), and the gradient, Eq. (7) are written above in their simplest forms. In the EnviFlux code there are two dierences that are not evident in these equations. The rst is that observations may be made at arbitrary times (and not just at the major time-steps). The second is that the gradient in Eq. (7) is not ecient as some factorisation can be done so that a separate adjoint of Mxb+δx 0→k does not need to be done for each k . 3.2.1 Intermediate observation times To allow intermediate observation times, dene a new operator, Hk that acts on the state between timesteps k−1 and k (in order to do an interpolations in time, and to allow for observations of the ux itself). 3 Figure 1: Working grid for the system. Panel (a) is for the horizontal grid. The labels for the u -latitudes are shown on the right of the panel and for the v -latitudes on the left. The labels for the v -longitudes are shown at the bottom and for the u -longitudes at the top. Panel (b) is for the vertical grid. The labels for the u, v levels (half-levels) are shown on the left and for the w -levels (full-levels) on the right. The tracer χ shares its longitude positions on the grid with v , its latitude positions with u and its level positions with u and v . The grid staggering is `Arakawa C' in the horizontal and `CharneyPhillips' grid in the vertical. 4 To do this we introduce the state labeled by a capital X(k∆t) , which contains tracer pairs ( χ([k−1]∆t) and χ(k∆t) ), and the uxes: X(k∆t) =   χ([k−1]∆t) χ(k∆t) ρ , (8) with X(0) =   0 χ(0) ρ . (9) This X(0) is dierent from x in (1), but X(0) can be constructed from x as follows: X(0) = Cx =  0 0 I 0 0 I x. (10) The observation operator acts as follows: HkXk=Hk  χ([k−1]∆t) χ(k∆t) ρ . (11) The time-stepping of the state vector is dened as X(k∆t) = M0→kX(0) = Mk−1→k· · · M1→2M0→1X(0) . (12) This (calligraphic Mk→k+1 ) propagates the tracer pairs, but keeps the uxes in the state vector. In terms of the propagation of the tracer χ(k∆t) to χ([k+ 1]∆t) , the time-stepping is (Roman Mk→k+1 ) δχ([k+ 1]∆t) = Mk→k+1δχ(k∆t) + SρI(k∆t/γ)∆t, (13) where S is dened after Eq. (3) and the function I(•) integerises its argument, so time k∆t is associated with ux time I(k∆t/γ) . Re-indexing (13) gives δχ(k∆t) = Mk−1→kδχ([k−1]∆t) + SρI([k−1]∆t/γ)∆t. (14) Putting this together gives the relationship between Mk−1→k and Mk−1→k : X(k∆t) = Mk−1→kX([k−1]∆t)   χ([k−1]∆t) χ(k∆t) ρ =  0 I 0 0 Mk−1→kS 0 0 I   χ([k−2]∆t) χ([k−1]∆t) ρ , (15) where S now also picks-out the relevant ux time. Alternatively, the above can be written X([k+ 1]∆t) = Mk→k+1X(k∆t)   χ(k∆t) χ([k+ 1]∆t) ρ =  0 I 0 0 Mk→k+1 S 0 0 I   χ([k−1]∆t) χ(k∆t) ρ . (16) The adjoint of this step is ˆ X([k−1]∆t):=M⊺ k−1→kˆ X(k∆t)   ˆ χ([k−2]∆t) ˆ χ([k−1]∆t) ˆ ρ :=  000 I M⊺ k−1→k0 0 S⊺I   ˆ χ([k−1]∆t) ˆ χ(k∆t) ˆ ρ  (17) 5 or ˆ X(k∆t):=M⊺ k→k+1 ˆ X([k+ 1]∆t)   ˆ χ([k−1]∆t) ˆ χ(k∆t) ˆ ρ :=  000 I M⊺ k→k+1 0 0 S⊺I   ˆ χ(k∆t) χ([k+ 1]∆t) ˆ ρ . (18) In terms of M , the cost function (6) is Jδx(δx) = 1 2δx⊺B−1δx | {z } Jb + 1 2 NT X k=0 (yk−HkM0→kC[xb+δx])⊺R−1 k(yk−HkM0→kC[xb+δx]) | {z } Jo , (19) where yt and Rk now relate to observations between time-steps k−1 and k (or just at 0 for k= 0 ). 3.2.2 Ecient gradient calculation The gradient of Jb is trivial: ∂Jb ∂δx=B−1δx. (20) The gradient of Jo is non-trivial and is found as follows. ∂Jo ∂δx=∂ ∂δx"1 2 NT X k=0 (yk−HkM0→kC[xb+δx])⊺R−1 k(yk−HkM0→kC[xb+δx])# =−C⊺ NT X k=0 M⊺ 0→kH⊺ kR−1 k(yk−HkM0→kC[xb+δx]) =C⊺ NT X k=0 M⊺ 0→kH⊺ k ∂Jo ∂ym k , (21) where ∂Jo/∂ym k=−R−1 k(yk−ym k) and ym t=HkM0→kC[xb+δx] . ∂Jo ∂δx=C⊺"M⊺ 0→1M⊺ 1→2· · · M⊺ k−1→k· · · M⊺ T−2→T−1M⊺ T−2→T−1M⊺ T−1→TH⊺ NT ∂Jo ∂ym NT + M⊺ 0→1M⊺ 1→2· · · M⊺ k−1→k· · · M⊺ T−2→T−1M⊺ T−2→T−1H⊺ NT−1∂Jo ∂ym NT−1+ M⊺ 0→1M⊺ 1→2· · · M⊺ k−1→kM⊺ T−2→T−1H⊺ NT−2∂Jo ∂ym NT−2+ . . .. . ..... . . M⊺ 0→1M⊺ 1→2· · · M⊺ k−1→kH⊺ k∂Jo ∂ym k . . .. . . M⊺ 0→1M⊺ 1→2H⊺ 2∂Jo ∂ym 2+ M⊺ 0→1H⊺ 1∂Jo ∂ym 1 H⊺ 0∂Jo ∂ym 0# 6 4 Building the master programs 4.1 Software required This software requires the following free software and libraries to be installed on the host system.  C++ compiler (gnu): this is produces most of the executables as EnviFlux is written in C++.  NetCDF library : this is needed to handle the input and output of elds.  FFTW library : this is needed to do the fast Fourier transforms.  SHTools library : this is needed to produce a data le of Gaussian co-latitudes, Gaussian weights, associated Legendre polynomials (for a spectral transform that is used by EnviFlux), and other data. This library is called only once per grid conguration and produces a le, cvt_ass_legendre_poly.dat , used by EnviFlux.  Fortran compiler (gnu): this is needed to produce the CVT les (standing for Control Variable Transform, which is the way that the B-matrix is encoded).  Python: this is needed for post-processing and to visualise the outputs of EnviFlux. 4.2 Contents of the EnviFlux download package  Calibration This directory contains ForTran code to produce the le that contains all the information to represent the B-matrix of the problem (the CVT le  see Sect. 4.1). The code is legacy ForTran code used to generate covariance les for INVICAT.  data This directory contains data needed to run the ux inversion system. This includes calibration les, wind data les, and the associated Legendre polynomial data le needed to perform spectral transforms.  docs This directory contains documentation, including this document.  graphics This directory contains python scripts to visualise output from EnviFlux.  src This directory contains the C++ source code of EnviFlux. Most of the system is written in C++. The only exception is the calibration code, which is written in Fortran-90 (this is the same code as used for InviCat, but EnviFlux does not use all the options that can be set), and the visualisation code, which is written in python. 4.3 Build commands To build one of the programs (each program is documented below), use the following command: make Master_program_name This will compile and link the required program. The running instructions are provided in each subsection. 5 Running the master test programs There is a suite of test programs to test coding of (i) read/write routines, (ii) adjoint routines, (iii) implied covariances, and (iv) the gradient of the cost function. In these programs, the input/output/specications are hardwired into the code. Here is the list of master test programs that do these tasks. 7  Master_TestReadWrite.cpp Tests reading and writing of state and wind data.  Master_AdjointTests_CVT.cpp Tests the adjoints of the following routines: halos, spherical, cvt_h, cvt_v, cvt_t, cvt_total.  Master_AdjointTestsObObs.cpp Tests the adjoint of the observation operator.  Master_AdjointTests_SemiLagrangian.cpp Tests the adjoint of the semi-Lagrangian advection scheme.  Master_AdjointTests_SL+ObsOp.cpp Tests the adjoint of the combined semi-Lagrangian advection scheme and the observation operator.  Master_ImpliedCov.cpp Tests the CVT to give some implied covariances.  Master_GradTest.cpp Tests the gradient produced by the pen-and-grad routine in EnviFlux. The remaining suite of programs concern the running of the system.  Master_MakeWinds.cpp Makes suitable wind les from ERA-5 winds to the format and resolution of EnviFlux. These drive the transport of the chosen tracer.  Master_MakeFields.cpp Generates initial conditions and surface ux elds for a run of the forward model.  Master_Invicat2EnviFlux.cpp Converts data from an INVICAT le to an EnviFlux le.  Master_ReplaceFirstField.cpp makes a copy of an EnviFlux le, but replacing the initial conditions from another le.  Master_ForwardTraj.cpp computes the trajectory of a particle advected by the prescribed winds.  Master_GenerateObs.cpp outpus model observations from a given EnviFlux le.  Master_MakeBG.cpp makes a background state by perturbing a specied truth with a random perturbation.  Master_Assimilate.cpp runs 4DVar using a covariance le, a bakckground le, and some observations to produce an analysis.  Master_CalcError.cpp calculates the dierences between two states. 6 Running the master core routines 6.1 Master_MakeWinds: Generating wind elds at the suitable resolution to drive the dynamics This program takes ECMWF ERA-5 wind elds and generates les of lower resolution winds. To build the code issue the command: make Master_MakeWinds . out in the directory containing the source code and makele. This will create the executable le Master_MakeWinds.out . To run this code, the user rst needs to prepare a number of things. The raw u , v , w winds need to be downloaded from the ECWMF by the user and placed in a specic directory, a wind le list needs to 8 be created, an output directory needs to be made, and a CVT le needs to be available (see below). The run command is ./Master_MakeWinds . out \ <CVT f i l e name ( for meta data)> \ <filename_containing_full = resolution_ECMWF_winds> \ <output_directory_for_reduced_resolution_winds>  The CVT le (or control variable transform le) is a netCDF le containing information that describes the background error covariance matrix used in the inverse problem. This le is specied here, even though no inverse problem is being solved by this module, as it contains the specication of the longitude/latitude/height grid that the ECMWF winds will be interpolated to (we call this meta data ). The grid staggering is `Arakawa C' in the horizontal and `CharneyPhillips' grid in the vertical. Fields u and v are both placed on so-called `half-levels', but at staggered locations in the horizontal, and w is placed on `full-levels' at yet another horizontally staggered position. See Fig. 1.  The lename containing the specications of the full-resolution ECMWF wind netcdf les (the wind le list) is a text le of the following strict format. Wind le list Comment line 1 Comment line 2 Comment line 3 Comment line 4 Comment line 5 Comment line 6 <Directory containing the ECMWF wind files > <Filename of 1 st set of winds> <No. of time s l i c e s in this f il e > <Filename of 2nd set of winds> <No. of time s li c es in this file > <Filename of 3rd set of winds> <No. of time s l i c e s in this f i le > ...  The directory name to output the reduced resolution winds. The full-resolution input ECMWF les are assumed to have a structure such that each le contains a dierent number of time slices. For example, if the ECMWF wind les have a time slice every 12 hours, and each input le contains a month of winds then the rst le (January) is expected to contain 62 wind slices, the second le (February) is expected to contain 56 wind slices, etc. The output (reduced resolution) winds, however, are output as one time slice per le. These will have lename Windsxxxx.nc , so Winds0000.nc to Winds0061.nc will be the wind les extracted for January, Winds0062.nc to Winds0117.nc will be the wind les for February, etc. The user will need to extract as many wind les as required for the length of the model run or assimilation window. 6.2 Master_MakeFields: Generate initial conditions and surface ux elds for a run of the forward model This program prepares a eld of initial conditions (3D) and surface ux (2D+time)  together called a state . The output state can be directly fed into (e.g.) the program that runs the forward model. To build the code, issue the command: 9 ===== If "s" ================================= <orbit inclination (deg)> \ <orbit period (min)> \ <t=0 long (deg)> \ <height (m) ( tracer only)> \ <separation in time (min)> \ <number in time> \ <error standard deviation> *** See special note below  Many of the program's arguments are the same as for Master_SemiLagrangian, see Sect. 6.5.  *** Warning: there may be an inconsistency between the entries marked above. If a departure points le is specied (departure points lename) to be read ((r)ead or (w)rite departure points le), the dynamical settings include vertical transport? and w factor will be ignored. In this case the dynamical settings that were used in the run that produced the departure points lename will be used. The output le of the observations has the following overall structure. General structure of observations le (example with 56 vertical levels) ===== Observation f i l e ===== nlevs : 56 ===== Mass profile ========= 001 29.470579 002 36.358128 ... 056 11.223949 ===== Observation number 1 ===== OBSERVATION ENTRY ===== Observation number 2 ===== OBSERVATION ENTRY ... The mass prole is related to the density of each level (needed for total columns observations), and is set by the code. No special ag is needed to the end list of observations (when the observations are read-in by other routines, the end of le is used to terminate the list). Each OBSERVATION ENTRY comprises a number of lines with a structure that indicates the obersation type. There are three observation types ( tracer , total column tracer , and ux ), which have structures as follows. 16 Observation entry for a tracer observation ob_of : t ob_type : i=individually specified obs , g=part of a grid of obs , s=part of a set of s a t e l l i t e obs ( for info only ) time : time_index day hour min secs t_alpha t_beta lon : lon_index longitude lon_alpha lon_beta lat : lat_index latitude lat_alpha lat_beta lev : lev_index level lev_alpha lev_beta ob : ob_value ob_err_stddev model_ob : model_ob_value innov : ob_value = model_ob_value grad : gradient Notes:  The time_index is not used for tracer observations, so is set to zero.  The t_alpha and t_beta are temporal interpolation coecients (these are set by the code).  The lon_index is the longitude index associated with this observation.  The lon_alpha, lon_beta are longitudinal interpolation coecients (these are set by the code).  The lat_index is the latitude index associated with this observation.  The lat_alpha, lat_beta are latitudinal interpolation coecients (these are set by the code).  The lev_index is the level index associated with this observation.  The lev_alpha, lev_beta are level interpolation coecients (these are set by the code).  The model_ob_value is the model's version of the observation. For the Master_GenerateObs routine, this contains the model observation corresponding to the true state (this is set by the code, as is the innovation).  The innov is the observation − model observation(-9999.0 means that it has not been set).  The gradient is the value of ∂Jo/∂ model_ob (this is set by the code, -9999.0 means that it has not been set). Observation entry for a total column tracer observation ob_of : t ob_type : i=individually specified obs , g=part of a grid of obs , s=part of a set of s a t e l l i t e obs ( for info only ) time : time_index day hour min secs t_alpha t_beta lon : lon_index longitude lon_alpha lon_beta lat : lat_index latitude lat_alpha lat_beta ob : ob_value ob_err_stddev model_ob : model_ob_value innov : ob_value = model_ob_value grad : gradient 17 The above notes apply (there is no level information needed for a total column observation). Observation entry for a ux observation ob_of : t ob_type : i=individually specified obs , g=part of a grid of obs , s=part of a set of s a t e l l i t e obs ( for info only ) time : time_index day hour min secs t_alpha t_beta lon : lon_index longitude lon_alpha lon_beta lat : lat_index latitude lat_alpha lat_beta lev : lev_index level lev_alpha lev_beta ob : ob_value ob_err_stddev model_ob : model_ob_value innov : ob_value = model_ob_value grad : gradient The above notes apply (there is no level information needed for a ux observation). The exception is the following:  The time_index is the time index associated with this observation (which indicates which ux eld is associated with the observation). 6.8 Master_MakeBG: Make a background state by incrementing a truth with a random perturbation This program generates a synthetic background state, xb , by adding a random perturbation, δxb , (with statistics consistent with background errors) on to a truth state, xt : xb=xt+Uχ, where δxb=Uχ , U being the control variable transform related to the background error covariance matrix B via, B=UU⊺ , and χ is a unit Gaussian random vector, χ∼N(0,I) . To build the code, issue the command: make Master_MakeBG. out in the directory containing the source code and makele. This will create the executable le Master_MakeBG.out . The run command is ./Master_MakeBG. out \ <truth filename ( i n i t i a l and boundary conds)> \ <output filename of background state> \ <output filename of background perturbation> \ <CVT file > \ <multiplication factor for tracer err std> \ <multiplication factor for flux err std> \ <multiplication factor for tracer pert> \ <multiplication factor for flux pert> 18  The CVT le is produced by the Calibration.f90 program. It contains the details that describe the background error covariance matrix.  The error standard deviations of the initial tracer and ux from the CVT le are multiplied by the respective std factors.  The generated background error perturbations of the initial tracer and ux are multiplied by the respective perturbation factors. 6.9 Master_Assimilate: Run a 4D-Var assimilation This program combines observations of the tracer (and/or ux if required) with a background state to produce a 4D-variational analysis. To build the code, issue the command: make Master_Assimilate . out in the directory containing the source code and makele. This will create the executable le Master_Assimilate.out . The run command is ./ Master_Assimilate . out \ <filename of background state> \ <wind directory> \ <obs filename> \ <departure points filename > *** \ <CVT file > \ <multipication factor for tracer err std> \ <multipication factor for flux err std> \ <obs output filename , inc model obs at bg> \ <output filename of analysis state> \ <output filename of analysis increment> \ <output filename of analysis diagnostics> \ <obs output filename , inc model obs at anal>\ <major time = step ( between winds , s)> \ <minor time = step ( integration , s)> \ <diffusion time = step ( s ) \ <horiz diffusion coefficient > \ <vert diffusion coefficient > \ <include vertical transport?> *** \ 0 or 1 <w factor > *** \ <interpolation type ( l or c)> \ linear or cubic <descent algorithm type (c or q)> \ c=conjugate grad <convergence criterion > \ <max number of var iterations > *** See special note below  The background le contains the background initial and boundary (surface ux) conditions.  The wind directory contains the sequence of driving winds. 19  The observation le contains the observations to be assimilated.  The departure points lename is a le of semi-Lagrangian departure points output for every major time-step. If this le does not exist at the start of the run it is created from the rst Var iteration.  The CVT le is produced by the Calibration.f90 program. It contains the details that describe the background error covariance matrix.  The multiplication factors respectively multiply the initial tracer and ux standard deviations read in from the CVT le.  The next observation le is for output; it contains the same information as the input le, but additionally with the model observations for the background, with innovation and gradient information. Specily nil to not output.  The analysis state le contains the 4D-Var analysis.  The analysis increment le contains the 4D-Var analysis increment.  The analysis diagnostics le contains the following information: iteration, Jb , Jo , J , ∥∇χJ∥2 , ∥∇χJ∥ .  The last observation le is for output; it contains the same information as the input le, but additionally with the model observations for the analysis, with residual and gradient information (the residual is placed in the eld labelled innov). Specily nil to not output.  The major time-step is the time-step between wind les (seconds).  The minor time-step is the time-step of the backward trajectories in the semi-Lagrangian scheme (seconds).  The diusion time-step is the time-step of the diusion process (seconds, applied after the semiLagrangian.  The horizontal diusion coecient is for the horizontal diusion scheme (0.0 for no horizontal diusion).  The vertical diusion coecient is for the vertical diusion scheme (0.0 for no vertical diusion).  Include the vertical dimension? (0 means use one level, 1 means use all vertical levels).  The w factor multiplies the vertical wind component by this number.  The interpolation type is either l for linear interpolation or c for cubic interpolation.  The interpolation type can be selected as either c for conjugate gradient or q for quasi-Newton (latter not yet implemented).  The convergence criterion sets the critical value of r=∥∇χJ∥2( iteration n)/∥∇χJ∥2( iteration 0 ) for termination of the iterations. For example if the criterion is set to 0.01 then the iterations will terminate when r < 0.01 .  The maximum number of iterations performs as specied.  *** Warning: there may be an inconsistency between the entries marked above. If a departure points le is specied (departure points lename), the dynamical settings include vertical transport? and w factor will be ignored. In this case the dynamical settings that were used in the run that produced the departure points lename will be used. 20 6.10 Master_CalcError: Calculate norms of the dierences between two states This program reads in two states and computes some norms of the dierence between them. To build the code, issue the command: make Master_CalcError . out in the directory containing the source code and makele. This will create the executable le Master_CalcError.out . The run command is ./ Master_CalcError . out \ <filename of state 1> \ <filename of state 2> <output filename> The following dierence norms are computed for the two states, x1 and x2 : norm1 = X p (x1(p)−x2(p)) norm2 = sX p (x1(p)−x2(p))2, where p represents the `location' in the eld. These norms are computed separately for the initial tracer and for the uxes. For the initial tracer, p represents longitude, latitude, and height; for the uxes, p represents longitude, latitude, and time. When state 2 is the `truth', the norms represent the error in state 1. 7 Running procedure Running and testing a system like EnviFlux is complex. In order to do an assimilation run, many les need to be produced rst. 1. Generate some wind les . How to do this is detailed in Sect. 6.1. You need to generate wind elds that cover the period of interest. 2. Generate a covariance le for the system . EnviFlux needs a CVT le, which is generated with the same FORTRAN program used for the INVICAT system. Some special options are needed for use with the toy system. a) The calibration program is at Calibration/Calibration.f90 . Here are the settings that need to be made to this le. If temporal correlations (for the source/sink eld) are required, set CVT_TemporalCors % temporal_covs ( f i el d ) = 1 otherwise CVT_TemporalCors % temporal_covs ( f i el d ) = 0 21 If temporal correlations are switched on (as above), set the shape of the temporal correlations (1=Lorentzian, 2=Gaussian, 3=SOAR, 4=exponential), e.g. CVT_TemporalCors % temporal_cor_tpe ( fie ld ) = 3 If temporal correlations are switched on (as above), set the timescale of the correlation (months), e.g. CVT_TemporalCors % timescale_flux ( fi eld ) = 1.0 Set the standard deviation of the initial eld, e.g. CVT_std % std_tracer ( lat , lon , lev ) = 15.0 Set the standard deviation of the source/sink eld, e.g. CVT_std % std_flux (1: ylat , 1: ylon , 1:nmonth , fie ld ) = 0.00001 Set the lengthscale of the initial eld (m), e.g. CVT_HorizCors % lengthscale_tracer (1: ylev ) = 600000.0 Set the shape of the horizontal correlations for the source/sink eld (1=Lorentzian, 2=Gaussian, 3=SOAR, 4=exponential), e.g. CVT_HorizCors % horizflux_cor_tpe ( fi e ld ) = 3 Set the lengthscale of the source/sink eld (m), e.g. CVT_HorizCors % lengthscale_flux (time1 , fi eld ) = 400000.0 b) The next le to edit is Calibration/main_data.f90 . Here are the settings that need to be made to this le. Make sure that the following are set (e.g. for L= 32 ) ylon=64, ylat =32, ylev=56 nmonth=37 22 c) The next le to edit is Calibration/onedvar_data.f90 . Here are the settings that need to be made to this le. Set the shape of the horizontal correlations for the initial eld (1=Lorentzian, 2=Gaussian, 3=SOAR, 4=exponential), e.g. horiztracer_cor_tpe = 3 Make sure that the following is set (to use global eigenvalues in the vertical covariances) vert_covs = 3 d) Go back to the directory containing the Calibration.90 le and run make Calibration . out ./ Calibration . out The le CVT_calib.nc is output, which can be moved elsewhere. 3. Generate a truth state . A prescription for doing this from scratch is detailed here. a) Optional: User the Master_WorldFlux utility to convert a given Invicat ux eld to the grid and units used by EnviFlux. This creates a state le containing the Invicat uxes, and with zero for the initial eld. b) Use the Master_MakeFields utility to generate a state eld with a constant initial eld (achieved by speciying a single blob with extremely large horizontal and vertical lengthscales), and localised source blobs. See Sect. 6.2. c) Run this eld forward in time using the Master_SemiLagrangian utility. See Sect. 6.5. d) Use the Master_ReplaceFirstField utility to take the le produced in step 3b above and replace the initial eld with the output from step 3c. See Sect. 6.4. e) Alternatively, use the Master_WorldFlux utility to generate elds of ux from an existing Invicat le. See Sect. 6.3. 4. Generate a background state . A prescription for producing a background state as a perturbed version of the initial eld of the true state, and zero background for the ux, use the following prescription. a) Use the Master_MakeFields utility to generate a state eld with zero values for the initial eld and ux. See Sect. 6.2. b) Use the Master_ReplaceFirstField utility to take the le produced in step 4a above and replace the initial eld with the output from step 3c when generating the truth. See Sect. 6.4. c) Use the Master_MakeBG utility to perturb the output of step 4b with random background perturbations specied in the covariance le. The standard deviations of the initial and ux elds can be scaled (separately). In order to have zero perturbations made to the ux eld, the scale for the ux can be set to zero. See Sect. 6.8. 5. Use the utility to Master_GenerateObs utility to generate some synthetic observations. Observations are generated by integrating the truth forward in time, computing sets of synthetic observations, and adding noise to produce observations that can later be assimilated. Observations can 23 be individual point measurements, regular grids of measurements, or patterns of observations that represent the nadir of a satellite orbit. Observations can be direct measurement of the tracer or the ux, or total column measurements of the tracer. See Sect. 6.7. 24 Index 4D-Var, 19 4DVar, 1 advection/diusion/ux equation, 3 analysis, 3 Arakawa C grid, 4, 9 associated Legendre polynomial, 7 augmented state vector, 2 back-trajectory, 2 background error covariance matrix, 3, 7, 9 blobs, 10 build commands, 7 building the master programs, 7 C++, 7 C++ compiler, 7 calculate norms between two states, 21 calibration les, 7 CharneyPhillips grid, 4, 9 compile and link, 7 convergence criterion, 20 cost function, 3, 6 CVT, 18 CVT le, 7, 911, 19, 20 departure points le, 14, 20 diusion stage, 2 diusion time-step, 2, 13, 20 documentation, 7 ECMWF ERA-5 wind elds, 8, 9 ecient gradient calculation, 6 FFTW library, 7 ux observation, 18 ForTran compiler, 7 forward model, 12 forward trajectory code, 14 full-levels, 4, 9 generate initial conditions and surface ux elds, 9 generate synthetic observations, 14 generating wind elds, 8 gradient, 6 gradient of the cost function, 3 graphics, 7 half-levels, 4, 9 horizontal diusion coecient, 2, 13, 20 horizontal grid, 4 intermediate observation times, 3 interpolate INVICAT le to EnviFlux, 11 INVICAT, 7 libraries, 7 major time-step, 2, 13, 14, 20 make, 8, 11, 12, 14, 15, 18, 19, 21 make a background state, 18 make a new state le, 12 mass prole, 16 master test programs, 7 Master_Assimilate, 19 Master_CalcError, 21 Master_ForwardTraj, 14 Master_GenerateObs, 14 Master_Invicat2EnviFlux, 11 Master_MakeBG, 18 Master_MakeFields, 9 Master_ReplaceFirstField, 12 Master_SemiLagrangian, 12 maximum number of iterations, 20 minor time-step, 2, 13, 14, 20 netCDF library, 7 observation error covariance matrix, 3 observation operator, 3 observations le, 16 post-processing, 7 python, 7 run 4D-Var assimilation, 19 running procedure, 21 semi-Lagrangian, 2 semi-Lagrangian advection, 2, 12, 14 SHTools library, 7 software required, 7 source time-step, 2 source/sink stage, 2 state vector, 2 synthetic background state, 18 synthetic observations, 14 testing, 21 time-steps, 2 total column tracer observation, 17 25