ROMS/TOMS Developers

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.

Corrected a couple of parallel and observation-processing bugs in the weak contraint algorithms:

• Corrected a parallel bug in the convolution of the adjoint solution in drivers: convolution.h, symmetry.h, w4dpsas_ocean.h, and w4dvar_ocean.h. Now, the adjoint solution is read in the tangent linear state arrays with a call to routine get_state, using iTLM instead of iADM in the calling arguments. Then, the routines load_TLtoAD and load_ADtoTL are used to load the solution to convolve into the appropriated state arrays before applying the adjoint and tangent linear squared-root diffusion operators, repectively. Only the interior solution is loaded (recall that get_state also puts the appropriate data in the ghost points). This is critial in parallel adjoint applications. Many thanks to Andy for finding this ellusive bug.
• Changed the order of the squared-root adjoint and tangent linear diffusion operators in the weak constraint drivers to be consistent with the computation of the error covariance normalization coefficients and other 4DVAR related algorithms. The squared-root adjoint operator is applied first for half of the diffusion steps. The resulting solution is then filtered with the tangent linear operator for the other half of diffusion steps. Recall, that both squared-root operators are used to impose symmetry.
• Corrected a parallel bug in the processing of the tangent linear vector TLmodVal in the weak constraint algorithms. This tangent linear vector is no longer read in routine obs_read to avoid accumulating their values with mp_collect during parallel exchanges. In weak constraint, the entire observation operator is maintained in memory. It includes all the observations for different time surveys. This requires additional logic in routines obs_write.F and ad_htobs.F. In addition, another parallel bug was fixed for more than one observation time survey applications. Only the relevant section of the observation vector is proccessed, for a particular time, using the datum indices Mstr and Mend.
• Corrected the processing of the tangent linear model forcing terms with the adjoint vector ADmodVal in ad_htobs.F. Previously, it was processing all time surveys at once in the f_* arrays. Now, only the relevant observation time survey is processed and applied as impulse forcing.

If you get a compilation error about missing tl_balance.F and/or ad_balance.F, just do a make depend to solve the problem.

Fixed a couple of bugs:

• Modified obs_initial.F, obs_read.F, and obs_write.F, so the observation arrays NLmodVal, TLmodVal and ObsScale are exchanged correctly in MPI applications. These arrays are now zero-out in obs_initial to allow exchanging data correctly between all tiles using mp_collect. Recall that the strategy here is to use a parallel reduction to process all these nontiled arrays. Otherwise, we will add nonzero and incorrect values to these processing observation arrays. Many thanks to Brian and Manu for reporting this problem.
• Fixed a bug in output.F and ad_output.F when creating multiple average or diagnostic terms NetCDF files. We were ignoring the first time-averaging cycle so the file was empty because another file was created before writing the first time-averaged cycle. Now, the creation of the new NetCDF files is delayed by one time step. This bug only applies to time-averaged NetCDF. Many thanks to Andy for reporting this bug.
• Added submit_is4dvar.bash which is required to run in some clusters. Therefore, we have the choice between csh and batch scripts.