ROMS/TOMS Developers

Made few changes to allow model-observation comparison statistics after the 4DVAR data assimilation. Many thanks to Andy and Brian for their suggestions.

• Modified is4dvar_ocean.h so the nonlinear model is run at the end after data assimilation, when using IS4DVAR. The nonlinear model is initialized with the estimated initial conditions from data assimilation. The nonlinear model is then interpolated at the observation locations for comparison statistics.
• Added a new routine Utility/stats_4dvar.F to compute model-observation comparisons. The following measurements are computed for each state variable:

  1) Observations mean (obs_mean) and standard deviation (obs_std).
  2) Model mean (model_mean) and standard deviation (model_std).
  3) Model bias (model_bias).
  4) Model-Observations standard deviation error (SDE).
  5) Model-Observations cross-correlation (CC).
  6) Model-Observations mean squared error (MSE).

  All these statistical quantities are written into 4DVAR output NetCDF file (MODname) and standard output file.

• Added stats_4dvar calls to all 4DVAR data assimilation drivers: is4dvar_ocean.h, s4dvar_ocean.h, w4dpsas_ocean.h, and w4dvar_ocean.h.
• Fixed an out-range index bound in arrayd A2d and A3d in the exact background error normalization coefficients routine normalization.F. This is a benign bug.
• Made changes to inp_par.F to insure that there is only one record in the restart file after the nonlinear model is run at the end of an assimilation cycle. The parameter nRST is set to ntimes. The restart file can be used as the first guess for the next assimilation file in sequential data assimilation.

Here are today’s changes:

• I updated the codes to write out initial nonlinear model at observation locations and initial and final misfit (innovation) vectors. These vectors are written into output 4DVAR NetCDF file (MODname) in variables NLmodel_initial, misfit_initial and misfit_final, respectively. The misfit vectors at output are equal to the model minus the observations and then divided by the observation’s standard deviation (SQRT(ObsError(:))). They are nondimensional. See obs_write.F for more details. The initial nonlinear model vector is used to compute several model-observations comparison statistics. These misfit vectors can be plotted to evaluate the minimization.
• I changed the following files to compute the misfit vectors described above: mod_ncparam.F, def_mod.F, obs_write.F, s4dvar_ocean.h, is4dvar_ocean.h, w4dpsas_ocean.h, and w4dvar_ocean.h. Also, I changed varinfo.dat. Warning: you need to use the new varinfo.dat, otherwise you can get a segmentation fault when creating the output 4DVAR NetCDF file.
• Following Brian’s suggestion, I deactivated the bulk flux computations in the tangent linear, representer, and adjoint models. The model is blowing-up in all the 4DVAR drivers. The bulk fluxes are still computed in the basic state when calling the nonlinear model. For now, the BULK_FLUX_NOT_YET was added to tl_main3d.F, rp_main3d, ad_main3d.F, tl_bulk_flux.F, rp_bulk_flux.F, ad_bulk_flux.F, tl_rho_eos.F, rp_rho_eos.F, and ad_rho_eos.F
• Made a small change in get_state.F so the time variable is ignored in the tangent linear, representer, and adjoint models when different than dstart. Recall, that usually these models are initialized from a generic zero-field NetCDF file in mostly all our adjoint applications. If this is not the case, the user needs to make sure that ocean_time and dstart are the same.

Corrected a couple of small bugs that affect the Generalized Stability Theory (GST) algorithms for CPP options AFT_EIGENMODES, FORCING_SV, FT_EIGENMODES, OPT_PERTURBATION.

• Corrected the declaration of input/output arrays A and Aout in routine mp_gather_state (in distribute.F) to A(Nstr:Nend) and Aout(Asize), respectively. The declarion was reversed between these two arrays.
• Corrected a small bug in PARPACK routine pdnaupd.f in line 482. The

  if ((ido .eq. 0) .or. (ido .eq. -2)) then

  statement had an .eq. instead of .or.; I think that I introduced this bug when implementing checkpointing.

• Changed get_state.F so the initial time is set to dstart. This allows to use the zero-initial NetCDF with any starting time. Recall that in several adjoint applications zero-initial fields are used for the tangent linear and adjoint models.

• Corrected a bug in congrad.F so the solution to the representer coefficients, cg_x, is loaded into vector ADmodVal at the end of the inner loop, even if the convergence criteria (CGeps) has not been satisfied. That is, the final inner loop solution is loaded if either the algorithm has converged or after Ninner iterations.
• Changed drivers w4dpsas_ocean.h and w4dvar_ocean.h to add an additional argument (Ninner) to the congrad call.

The W4DPSAS algorithm seems to be working well now. It has better convergence behaviour than the W4DVAR algorithm. The stability of the finite amplitude tangent linear model (representer tangent linear model) is affecting the convergence of the W4DVAR algorithm. Some improvement can be achieved by activating the diffusive relaxation in the Picard iterations (RPM_RELAXATION).

The remaining issue in both W4DPSAS and W4DPSAS algorithms is the conjugate gradient preconditioning. The convergence is affected as the number of observation increases.

Fixed a bug in routine nf_fread2d.F when interpolating coarser data (usually 2D forcing fields) into the model grid. Many thanks to Kate for finding this bug. In the past, I was assuming that the coarser data at input was smaller in size than the finner grid to interpolate. Now, the scratch work array wrk is allocated to the needed length.