Data assimilation algorithms

3D-Var

The default upper-air data assimilation algorithm is three-dimensional variational assimilation (3D-Var). To use 3D-Var no changes to ecf/config_exp.h should be required assuming structure function data files are available for your domain. Structure function input is defined in scr/include.ass.

Settings

List of 3D-Var settings that the user should be concerned about.

4D-Var

With 43h2.2 four-dimensional variational assimilation (4D-Var) is available as a non-default option. In order to setup an experiment to use 4D-Var one should issue the following commands:

mkdir -p /path/to/home/hm_home/MY_EXP
cd /path/to/home/hm_home/MY_EXP
/path/to/Harmonie/config-sh/Harmonie setup -r /path/to/Harmonie -c AROME_4DVAR

To use 4D-Var no changes to ecf/config_exp.h should be required assuming suitable structure function data files are available for your domain (see also ILRES in the settings section). Structure function input is defined in scr/include.ass.

Settings

List of 4D-Var settings that the user should be concerned about.

With the following settings the working of the 4D-Var can be changed. Defaults values are given

NOUTERLOOP=2                            # number of 4DVAR outer loops
ILRES=6,3                               # Resolution (in parts of full) of outer loops as compared to the forecast resolution. The domain NLATxNLON should have the property that for the settings of ILRES  both NLON/ILRES and NLAT/ILRES are of the form  5^c^ 3^d^ 2^e^, where c, d and e are integers >= 0 and e>=1.
TSTEP4D=300,150                         # Timestep length (seconds) of outer loops TL+AD
TSTEPTRAJ=300,300                       # How often the model state is saved for linearization
NITER4D=10,15                           # Maximum number of inner loop iterations in 4D-Var outer loops
NSIMU4D=15,20                           # Maximum number of cost function simulations in 4D-Var outer loops
CH_RES_SPEC=yes                         # yes => change of resolution of the increment spectrally; no => by FULLPOS
FORCE1=no                               # yes => tendency increment; no => analysis increment in loop 1
FORCE2=no                               # yes => tendency increment; no => analysis increment in loop 2

Schematic work flow of 4D-Var

• After the screening 4DVscreen for the default 3h observation window (2h before and 1h after the start of the main forecast),
• 4DVprolog prepares the initial conditions at the appropriate resolution (ILRES settings) for the forecasts used in minimization. The input here is the background BGHR (ICMSHHARM+0001 fields of the former cycle supplemented with some surface fields).
• Subsequently 4DVminim produces initially the low resolution at the beginning of the observation window (an_lowres_$loop), which is then (CH_RES_SPEC=yes in config_exp.h) transformed to the field at forecast resolution (an_hr_begwin).
• This field is complemented in Blendhr with necessary surface fields and the resulting field (anb_hr_begwin) acts as the initial condition for the trajectory run 4DVtraj.
• From the 2h forecast of 4DVtraj (ICMSHTRAJ+0002 fields) the main forecast is started.

Cloudingest - Cloud Initialization

The cloud initialization modifies temperature and humidity fields with help of MSGcloudYYYYMMDDHH.grib file, which contains 2-d fields of cloudtop-temperature [K], cloudmask [0-1] and cloudbase [m]. Pre cy46h the Cloudingest happened within MASTERODB in src/arpifs specifically in src/arpifs/phys_dmn/msginit.F90 routine. Since cy46h the Cloudingest uses pysurfex and gl to do the job.

Settings

Enable Cloudingest in ecf/config_exp.h:

USE_MSG=yes                              # Use MSG data for adjustment of inital profiles, EXPERIMENTAL! (no|yes), expects MSGcloudYYYYMMDDHH.grib in $OBDIR

And in src/include.ass:

export CLOUD_BASE=1                # 1 and USE_MSG=true (ecf/config_exp.h) => msginit

to replace the cloudbase field in MSGcloudYYYYMMDDHH.grib with an OI interpolated field of SYNOP observations of CloudBaseHeights and postprocessed field of cloudbases of the first-guess file.

Ecflow and schematic work flow of Cloudingest

In case of 3DVar assimilation, the ecflow AnUA family should look like this:

And a schematic work-flow of tasks and files involved:

Incremental Analysis Updates - IAU

This is a modified and experimental version of the default implemented Forward-Facing Incremental Analysis Update, that is described elsewhere (and mainly controlled via namelist settings). This algoritm is using the same code but is changed in a few fundamental ways on the logistical plane. Instead of introducing the innovations gradually in the forecast from the analysis point (as in default 3Dvar) the analysis is done at "t=0" and then the forecast is started at an earlier point so that the center of the introductions of the innovations are located at the above mentioned analysis point.

This is done through shifting the forecast start to 1h before "t=0" and activating a namelist change in the forecast namelist:

NAMIAU=>{
  'LIAU' => '.TRUE.,',
  'TSTARTIAU' => '1800.0,',
  'TSTOPIAU' => '5400.0,',
  'ALPHAIAU' => '1.0,',

Hence the forecast is running freely for 30min and then starts adding the increments during 1h centered around the "t=0" point, during 1h. If these times are subject of change (for example if a larger window is desired, mind that you also need to have change what files are linked to various process such as the forecast so that the correct start files are linked to that process.

Settings

If you want to try this method, or want to test it in order to develop it further there is a few thing to keep in mind. Before doing anything else you have to modify the write up times in ecf/config_exp.h so that the model have access to the startfiles that it needs inorder to start.

From the DEFAULT:

    LL_LIST="12,3"                          # Forecast lengths for the cycles [h], replaces LL, LLMAIN
    HWRITUPTIMES="00-21:3,24-60:6"          # History file output times
    FULLFAFTIMES=$HWRITUPTIMES              # History FA file IO server gather times
    PWRITUPTIMES="00-60:3"                  # Postprocessing times
    PFFULLWFTIMES=$PWRITUPTIMES             # Postprocessing FA file IO server gathering times
    VERITIMES="00-60:1"                     # Verification output times, changes PWRITUPTIMES/SFXSELTIMES
    SFXSELTIMES=$PWRITUPTIMES               # Surfex select file output times

to:

    LL_LIST="12,6"                          # Forecast lengths for the cycles [h], replaces LL, LLMAIN
    HWRITUPTIMES="00-06:1,06-21:3,24-60:6"  # History file output times
    FULLFAFTIMES=$HWRITUPTIMES              # History FA file IO server gather times
    PWRITUPTIMES="00-06:1,06-60:3"          # Postprocessing times
    PFFULLWFTIMES=$PWRITUPTIMES             # Postprocessing FA file IO server gathering times
    VERITIMES="00-60:1"                     # Verification output times, changes PWRITUPTIMES/SFXSELTIMES
    SFXSELTIMES=$PWRITUPTIMES               # Surfex select file output times
                                            # Only meaningful if SURFEX_LSELECT=yes
    SFXSWFTIMES=$SFXSELTIMES                # SURFEX select FA file IO server gathering times
    SWRITUPTIMES="00-06:1"                  # Surfex model state output times
    SFXWFTIMES=$SWRITUPTIMES                # SURFEX history FA file IO server gathering times

• Then run one cycle for at least 6h so that all the prerequesite files are generated before turning on IAUVAR.

• Second thing, when the previous run has completed, is to set the IAUVAR ecf/config_exp.h to yes and keep running. It should be automatic to continue as usual. Make sure that you do not manually set LL to be shorter then 6h.

The surface assimilation is moved to the start of the forecast and hence it is only the upper air assimilation that is involved in the IAU.

Flow diagram of IAU (Magnus will help)

Variational Constraints (VC)

Balances among analysed meteorological variables in the data assimilation process can be effectively introduced by means of the well-known method of variational constraints (VC). In the case of the ALADIN-NH dynamics, its semi-implicit linear system for the non-hydrostatic fully compressible Euler equations (SI) appears as a convenient way of giving a precise definition to these constraints. An interesting feature of this method is the integration in the analysis of the vertical velocity field, which clearly must be important in convection permitting NWP. In this on-line doc resource only practical aspects on how to use this option will be mentioned, comprehensive information on this method can be found , for instance, in Variational Constraints for Data Assimilation in ALADIN-NH Dynamics.

The VC algorithm is activated by setting TRUE the LVARCONS switch in ecf/config_exp.h

# **** Assimilation ****
LOOPS="no"
ANAATMO=3DVAR                           # Atmospheric analysis (3DVAR|4DVAR|blending|none)
LVARCONS='.TRUE.'                       # apply VC to analysis increments for a better balanced initial state
IAUVAR=no

Obviously this will have any effect on the analysis only when 3DVAR is active and other initialization schemes (e.g. DFI) are off. The VC algorithm is applied after the 3DVAR minimization to the analysis increments. When any of the other options for ANAATMO are active, VC will not be applied. Application of DFI or other filters on the analysis will be redundant as VC is supposed to remove much of the noise in principle present in the t=0 dynamical fields: horizontal and vertical divergence, temperature, pressure departure and surface pressure. Note that vorticity is not part of the SI system, and therefore is not processed by VC.

The whole VC run configuration is controlled by means of a single the namelist NEMVARCONS to be found in nam/harmonie_namelists.pm . The default values are:

%minimization=(
...

NEMVARCONS=>{
'VCVERBOSELEV' => '0,',
'LVCTESTID' => '.FALSE.,',
'lVCPLOTGF' => '.FALSE.,',
'LVCFHD' => '.TRUE.,',
'LVCFT'  => '.TRUE.,',
'LVCFPS' => '.TRUE.,',
'LVCFGW' => '.FALSE.,',
'VCWEIGHT' => '10.0,',
'VCWEIGHTHD' => '-1.50,',
'VCWEIGHTT' => '-1.50,',
'VCWEIGHTPS' => '-1.50,',
'VCWEIGHTGW' => '-1.50,',
'LVCFGWSCALE' => '.TRUE.,',
'VCLOGPS' => '0,',
'VCTSTEP' => '60.0,',
'VCPSREF' => '90000.0,',
'VCTREF' => '350.0,',
'VCTREFA' => '70.0,',
},
)

Modification of these values should be done in nam/harmonie_namelists.pm itself. To use effectively VC, the user has to be aware of only a subset of these parameters :

'VCVERBOSELEV' => '0,', # for more detailed printouts regarding VC , raise this value ( up to 4 )

'VCWEIGHT' => '10.0,',  # Critical parameter of the algorithm. Gives the relative weight between increments to be
                        # filtered and dynamical constraints. Values should be between 0.5 (zero not allowed !) and 10.
                        # The higher the value, the less constrained is the analysis. This default value of 10
                        # guarantees that final analysis increments are very close to 3DVAR increments.
                        # In the many different tests done during validation and verification experiemnts,
                        # a value of about 1.0 has given best results.

'VCTSTEP' => '60.0,',     # These parameters correspond to the SI parameters and to the time step. Although the  
'VCPSREF' => '90000.0,',  # SI in forecast mode is different from the SI used by VC ( in particular the numerical   
'VCTREF' => '350.0,',     # method is different ) these values should agree with those used in forecast mode for  
'VCTREFA' => '70.0,',     # optimal performance
},

Maybe the only other information that is required to effectively use this VC algorithm, concerns the parameters :

'VCWEIGHTHD' => '-1.50,',
'VCWEIGHTT' => '-1.50,',
'VCWEIGHTPS' => '-1.50,',

They enable the flexibility of considering different weights for different 3DVAR analysis increments. When VC operates with a low VCWEIGHT value (strongly constrained mode), it can remove some overfitting to wind, temperature and/or surface pressure observations, and this may produce an apparent degradation in operational verification curves close to t=0. These three parameters permit adjust individually each variable. Note that negative values ( as in the default ) automatically reverts to the VCWEIGHT value.