ROMS/TOMS Developers

I’m trying this option for the first time and came across a little bug: the horizontal grid for the global CCSM CORE files have more points than my little ROMS grid. In nf_fread2d.F, wrk is dimensioned:

real(r8), dimension(2+(Lm(ng)+2)*(Mm(ng)+2)) :: wrk

This isn’t big enough in my case. Why not make it allocatable?

Also, I would dearly love to see the winds handled in the same way as all the other fields. I know it would get tricky if we wanted to say use WRF on it’s rotated grid to drive ROMS on another rotated grid, but CCSM should be simple enough. I’ll look into it.

The main goal of strong contraint 4DVar is to estimate optimal initial conditions by melding the first guess (background) state and available data for the selected assimilation time window using both the statistics of the background and observations. In the incremental algorithm, IS4DVAR, the estimated initial conditions is saved in NetCDF file INIname. This NetCDF file will have two time records at the end of the assimilation run. The first record is the new estimated initial conditions from data assimilation whereas the second record is the starting first guess. Noticed that in the driver is4dvar_ocean.h these records are, respectively, processed with index variables Lini=1 and Lbck=2. In the IS4DVAR algorithm both records are used extensively to update the background state with the estimated data assimilation increments.

I have been asked why I do store these two records in this order. Obviously, it makes more sense to store the background first (record 1) and the new estimated initial conditions on record 2 of the unlimited time-dimension. Well, it is not done in this way because this file will be used latter to initialize the nonlinear model in forecast mode. It is easier to start from record 1 than 2 without special manipulations of the standard input parameters. The default in ROMS is to read initial conditions from the first record of the NetCDF file, if not specified otherwise. This strategy is extremely important in sequential data assimilation. See previous post about the sequential IS4DVAR running script, submit_is4dvar.sh.

I made a couple of changes to the code to allow real-time sequential assimilation with the incremental 4DVAR algorithm (IS4DVAR):

• Added a generic csh script, Bin/submit_is4dvar.sh, to carry out sequential 4DVAR experiments. That is, the state estimation is split in several assimilation cycles (time windows) as new data becomes available. This script runs both the IS4DVAR algorithm and nonlinear model. It initializes the nonlinear model with the estimated optimal initial conditions from the assimilation cycle. The nonlinear model is run to end of the assimilation cycle and then the restart file becomes the first guess (background state) for the next assimilation cycle. The user needs to have the following directory structure:

       $MYROOT/                    root directory
              /Data                Application configuration Input NetCDF files
              /Forward             Nonlinear model working directory
              /IS4DVAR             Data assimilation working directory
              /OBS                 Observations NetCDF file
              /Storage             Storage directory for estimation output NetCDF files
    

  See above script for more details. This script is generic and well documented.

• Removed the writing of the cost function information to the nonlinear and tangent linear initial conditions NetCDF files. This allows to use the restart file in sequential data assimilation. The diagnostic cost function information is still written to the MODname NetCDF file. This also cleaned the initial conditions files used in the 4DVAR algorithms. Therefore, the following files were changed def_ini.F, tl_def_ini.F, wrt_ini.F, and tl_write_ini.F.
• Modified inp_par.F to write horizontal and vertical convolution stability paramenters Hgamma and Vgamma to standard output.

I made few changes:

• Introduced changes in get_state.F to allow reading vertical viscosity and diffusion coefficients, if available, during initialization when GLS_MIXING, LMD_MIXING, or MY25_MIXING are activated. Otherwise, during restart or 4DVAR these coefficients will have the background value during the predictor corrector which require sometimes to reduce the time-step to avoid the model to blow-up.
• Added new CPP option RAMP_TIDES to activate ramping of tidal forcing for one day during initialization. This is usually a good idea to avoid a sharp shock to the system. This was hardwired in routine set_tides.F.
• Added open boundary sponge for applications SW06_COARSE and SW06_FINE in routine horz_mix.F.

Corrected few routines to allow assimilating observations from a master NetCDF file containing data outside the assimilation time window:

• Fixed a bug in obs_initial.F to allow processing of the observations from a NetCDF file containing surveys outside of the assimilation time window. The model now search for the available observations and ignores outbound time records. This is nice because there is not need to write the observation file to specific time records used in the assimilation window. Only one observation file is needed. Since the datum dimension is unlimited, this observation file can grow in real-time applications.
• Modified obs_write.F, ad_misfit.F, and ad_htobs.F to report better the observations as they are processed. I was not reporting the observation processing information during the adjoint execution. Now will see something like:

     Number of State Observations Processed: ObsTime = 194.8962
  
     Variable IstrObs IendObs Count Rejected
  
     temp 1047 1084 38 0
     salt 1085 1122 38 0
  
     Total 76 0
     Obs Tally 986 0
  
     Wrote background state at observation locations, datum = 0001047 - 0001122
  

  Notice that now the processing of each state variable is in terms of datum record in the observation NetCDF file. Also, I am now reporting the tally count of all observations processed.

• Corrected a segmentation fault in nf_fread2d.F and nf_fread3d.F in some new cluster systems. The problem was with the inconsistency in the use of the intent(inout) for gtype and all the subroutine calls.