Difference between revisions of "RBL4DVar Analysis Observation Impact Tutorial"
m (Text replacement  "W4DPSAS_SENSITIVITY" to "RBL4DVAR_ANA_SENSITIVITY") (change visibility) 
m (Robertson moved page PSAS Analysis Observation Impact Tutorial to RBL4DVar Analysis Observation Impact Tutorial: PSAS renamed to RBL4DVar) (change visibility) 
(No difference)

Revision as of 20:36, 23 July 2020
Tutorial Menu 

1. Introduction 
2. Error Covariance Normalization 
3. I4DVar 
4. RBL4DVar 
5. RBL4DVar Analysis Observation Impact 
6. RBL4DVar Analysis Observation Sensitivity 
7. Array Modes 
8. RBL4DVar Forecast Observation Impact 
9. RBL4DVar Forecast Observation Sensitivity 
Contents
Introduction
During this exercise you will apply the strong/weak constraint, dual form 4Dimensional Variational (4DVar) data assimilation observation impact based on the Physicalspace Statistical Analysis System (PSAS) algorithm to ROMS configured for the U.S. west coast and the California Current System (CCS). This configuration, referred to as WC13, has 30 km horizontal resolution, and 30 levels in the vertical. While 30 km resolution is inadequate for capturing much of the energetic mesoscale circulation associated with the CCS, WC13 captures the broadscale features of the circulation quite well, and serves as a very useful and efficient illustrative example of 4DPSAS observation impact.
Model Setup
The WC13 model domain is shown in Fig. 1 and has open boundaries along the northern, western, and southern edges of the model domain.
In the tutorial, you will perform a 4DVar data assimilation cycle that spans the period 36 January, 2004. The 4DVar control vector δz is comprised of increments to the initial conditions, δx(t_{0}), surface forcing, δf(t), and open boundary conditions, δb(t). The prior initial conditions, x_{b}(t_{0}), are taken from the sequence of 4DVar experiments described by Moore et al. (2011b) in which data were assimilated every 7 days during the period July 2002 December 2004. The prior surface forcing, f_{b}(t), takes the form of surface wind stress, heat flux, and a freshwater flux computed using the ROMS bulk flux formulation, and using near surface air data from COAMPS (Doyle et al., 2009). Clamped open boundary conditions are imposed on (u,v) and tracers, and the prior boundary conditions, b_{b}(t), are taken from the global ECCO product (Wunsch and Heimbach, 2007). The freesurface height and vertically integrated velocity components are subject to the usual Chapman and Flather radiation conditions at the open boundaries. The prior surface forcing and open boundary conditions are provided daily and linearly interpolated in time. Similarly, the increments δf(t) and δb(t) are also computed daily and linearly interpolated in time.
The observations assimilated into the model are satellite SST, satellite SSH in the form of a gridded product from Aviso, and hydrographic observations of temperature and salinity collected from Argo floats and during the GLOBEC/LTOP and CalCOFI cruises off the coast of Oregon and southern California, respectively. The observation locations are illustrated in Fig. 2.
Running 4DPSAS Observation Impact
To run this exercise, go first to the directory WC13/PSAS_impact. Instructions for compiling and running the model are provided below or can be found in the Readme file. The recommended configuration for this exercise is one outerloop and 50 innerloops, and roms_wc13.in is configured for this default case. The number of innerloops is controlled by the parameter Ninner in roms_wc13.in.
Important CPP Options
The following Cpreprocessing options are activated in the build script:
ANA_SPONGE Analytical enhanced viscosity/diffusion sponge
AD_IMPULSE Force ADM with intermittent impulses
BGQC Backgound quality control of observations
IMPACT_INNER Writing observations impacts for each inner loop
MINRES Minimal Residual Method for 4DVar minimization
OBS_IMPACT Compute observation impact
OBS_IMPACT_SPLIT separate impact due to IC, forcing, and OBC
RPCG Restricted Bpreconditioned Lanczos minimization
TIME_CONV Weakconstraint 4DVar time convolution
WC13 Application CPP option
Input NetCDF Files
WC13 requires the following input NetCDF files:
Nonlinear Initial File: wc13_ini.nc
Forcing File 01: ../Data/coamps_wc13_lwrad_down.nc
Forcing File 02: ../Data/coamps_wc13_Pair.nc
Forcing File 03: ../Data/coamps_wc13_Qair.nc
Forcing File 04: ../Data/coamps_wc13_rain.nc
Forcing File 05: ../Data/coamps_wc13_swrad.nc
Forcing File 06: ../Data/coamps_wc13_Tair.nc
Forcing File 07: ../Data/coamps_wc13_wind.nc
Boundary File: ../Data/wc13_ecco_bry.nc
Adjoint Sensitivity File: wc13_ads.nc
Initial Conditions STD File: ../Data/wc13_std_i.nc
Model STD File: ../Data/wc13_std_m.nc
Boundary Conditions STD File: ../Data/wc13_std_b.nc
Surface Forcing STD File: ../Data/wc13_std_f.nc
Initial Conditions Norm File: ../Data/wc13_nrm_i.nc
Model Norm File: ../Data/wc13_nrm_m.nc
Boundary Conditions Norm File: ../Data/wc13_nrm_b.nc
Surface Forcing Norm File: ../Data/wc13_nrm_f.nc
Observations File: wc13_obs.nc
Lanczos Vectors File: wc13_lcz.nc
Various Scripts and Include Files
The following files will be found in WC13/PSAS_impact directory after downloading from ROMS test cases SVN repository:
build_roms.bash bash shell script to compile application
build_roms.sh csh Unix script to compile application
job_psas_impact.sh job configuration script
roms_wc13.in ROMS standard input script for WC13
s4dvar.in 4DVar standard input script template
wc13.h WC13 header with CPP options
Important parameters in standard input roms_wc13.in script
 Notice that this driver uses the following adjoint sensitivity parameters (see input script for details):
 DstrS == 0.0d0 ! starting day
DendS == 0.0d0 ! ending day
KstrS == 1 ! starting level
KendS == 30 ! ending level
Lstate(isFsur) == T ! freesurface
Lstate(isUbar) == T ! 2D Umomentum
Lstate(isVbar) == T ! 2D Vmomentum
Lstate(isUvel) == T ! 3D Umomentum
Lstate(isVvel) == T ! 3D Vmomentum
Lstate(isWvel) == F ! 3D Wmomentum
Lstate(isTvar) == T T ! tracers
 Both FWDNAME and HISNAME must be the same:
Instructions
To run this application you need to take the following steps:
 We need to run the model application for a period that is long enough to compute meaningful circulation statistics, like mean and standard deviations for all prognostic state variables (zeta, u, v, T, and S). The standard deviations are written to NetCDF files and are read by the 4DVar algorithm to convert modeled error correlations to error covariances. The error covariance matrix, D=diag(B_{x}, B_{b}, B_{f}, Q), is very large and not well known. B is modeled as the solution of a diffusion equation as in Weaver and Courtier (2001). Each covariance matrix is factorized as B = K Σ C Σ^{T} K^{T}, where C is a univariate correlation matrix, Σ is a diagonal matrix of error standard deviations, and K is a multivariate balance operator.../Data/wc13_std_i.nc initial conditions
../Data/wc13_std_m.nc model error (if weak constraint)
../Data/wc13_std_b.nc open boundary conditions
../Data/wc13_std_f.nc surface forcing (wind stress and net heat flux)  Since we are modeling the error covariance matrix, D, we need to compute the normalization coefficients to ensure that the diagonal elements of the associated correlation matrix C are equal to unity. There are two methods to compute normalization coefficients: exact and randomization (an approximation).Nmethod == 0 ! normalization method: 0=Exact (expensive) or 1=Approximated (randomization)These normalization coefficients have already been computed for you (../Normalization) using the exact method since this application has a small grid (54x53x30):
Nrandom == 5000 ! randomization iterations
LdefNRM == T T T T ! Create a new normalization files
LwrtNRM == T T T T ! Compute and write normalization
CnormM(isFsur) = T ! model error covariance, 2D variable at RHOpoints
CnormM(isUbar) = T ! model error covariance, 2D variable at Upoints
CnormM(isVbar) = T ! model error covariance, 2D variable at Vpoints
CnormM(isUvel) = T ! model error covariance, 3D variable at Upoints
CnormM(isVvel) = T ! model error covariance, 3D variable at Vpoints
CnormM(isTvar) = T T ! model error covariance, NT tracers
CnormI(isFsur) = T ! IC error covariance, 2D variable at RHOpoints
CnormI(isUbar) = T ! IC error covariance, 2D variable at Upoints
CnormI(isVbar) = T ! IC error covariance, 2D variable at Vpoints
CnormI(isUvel) = T ! IC error covariance, 3D variable at Upoints
CnormI(isVvel) = T ! IC error covariance, 3D variable at Vpoints
CnormI(isTvar) = T T ! IC error covariance, NT tracers
CnormB(isFsur) = T ! BC error covariance, 2D variable at RHOpoints
CnormB(isUbar) = T ! BC error covariance, 2D variable at Upoints
CnormB(isVbar) = T ! BC error covariance, 2D variable at Vpoints
CnormB(isUvel) = T ! BC error covariance, 3D variable at Upoints
CnormB(isVvel) = T ! BC error covariance, 3D variable at Vpoints
CnormB(isTvar) = T T ! BC error covariance, NT tracers
CnormF(isUstr) = T ! surface forcing error covariance, Umomentum stress
CnormF(isVstr) = T ! surface forcing error covariance, Vmomentum stress
CnormF(isTsur) = T T ! surface forcing error covariance, NT tracers fluxes../Data/wc13_nrm_i.nc initial conditionsNotice that the switches LdefNRM and LwrtNRM are all false (F) since we already computed these coefficients.
../Data/wc13_nrm_m.nc model error (if weak constraint)
../Data/wc13_nrm_b.nc open boundary conditions
../Data/wc13_nrm_f.nc surface forcing (wind stress and
net heat flux)  Before you run this application, you need to run the standard 4DPSAS (../PSAS directory) since we need the Lanczos vectors. Notice that in job_psas_impact.sh we have the following operation:cp p ${Dir}/PSAS/wc13_mod.nc wc13_lcz.ncIn 4DVar (observartion space minimization), the Lanczos vectors are stored in the output 4DVar NetCDF file wc13_mod.nc.
 In addition, to run this application you need an adjoint sensitivity functional. This is computed by the following Matlab script:../Data/adsen_37N_transport.mwhich creates the NetCDF file wc13_ads.nc. This file has already been created for you.
 Customize your preferred build script and provide the appropriate values for:
 Root directory, MY_ROOT_DIR
 ROMS source code, MY_ROMS_SRC
 Fortran compiler, FORT
 MPI flags, USE_MPI and USE_MPIF90
 Path of MPI, NetCDF, and ARPACK libraries according to the compiler are set in my_build_paths.sh. Notice that you need to provide the correct places of these libraries for your computer. If you want to ignore this section, set USE_MY_LIBS value to no.
 Notice that the most important CPP options for this application are specified in the build script instead of wc13.h:setenv MY_CPP_FLAGS "${MY_CPP_FLAGS} DRBL4DVAR_ANA_SENSITIVITY"This is to allow flexibility with different CPP options.
setenv MY_CPP_FLAGS "${MY_CPP_FLAGS} DANA_SPONGE"
setenv MY_CPP_FLAGS "${MY_CPP_FLAGS} DAD_IMPULSE"
setenv MY_CPP_FLAGS "${MY_CPP_FLAGS} DOBS_IMPACT"
setenv MY_CPP_FLAGS "${MY_CPP_FLAGS} DOBS_IMPACT_SPLIT"
setenv MY_CPP_FLAGS "${MY_CPP_FLAGS} DRPCG"  You MUST use the build script to compile.
 Customize the ROMS input script roms_wc13.in and specify the appropriate values for the distributedmemory partition. It is set by default to:Notice that the adjointbased algorithms can only be run in parallel using MPI. This is because of the way that the adjoint model is constructed.
 Customize the configuration script job_psas_impact.sh and provide the appropriate place for the substitute Perl script:set SUBSTITUTE=${ROMS_ROOT}/ROMS/Bin/substituteThis script is distributed with ROMS and it is found in the ROMS/Bin subdirectory. Alternatively, you can define ROMS_ROOT environmental variable in your .cshrc login script. For example, I have:setenv ROMS_ROOT /home/arango/ocean/toms/repository/trunk
 Execute the configuration job_psas_impact.sh before running the model. It copies the required files and creates psas.in input script from template s4dvar.in. This has to be done every time that you run this application. We need a clean and fresh copy of the initial conditions and observation files since they are modified by ROMS during execution.
 Run ROMS with data assimilation:mpirun np 8 romsM roms_wc13.in > & log &
 We recommend creating a new subdirectory EX5, and saving the solution in it for analysis and plotting to avoid overwriting solutions when playing with different parameters. For examplemkdir EX5where log is the ROMS standard output specified in the previous step.
mv Build_roms psas.in *.nc log EX5
cp p romsM roms_wc13.in EX5
Results
The WC13/plotting/plot_psas_analysis_impact.m Matlab script will allow you to plot the 4DPSAS analysis observation impacts: