IMPORTANT: Updated several Matlab scripts for nesting

[arango]

This is an important update to several Matlab scripts since several scripts associated with nesting were update. Users need to update their nesting connectivity NetCDF file.

The following Matlab scripts were updated for new capabilities or to correct bugs:

• 4dvar/super_obs.m: corrected statement at line 310 to be

    Sout.origin_flag_values = Sinp.origin_flag_values;
  
  instead of
  
    Sout.origin_flag_values = Sinp.state_flag_values;
  

• boundary/obc_roms2roms.m: corrected conditional at line 296 to be

     if (isvec && R.vector_rotation),
  
  instead of
  
     if (isvec && D.vector_rotation),
  

• grid/c_contact.m: renamed linear interpolation weights from Hweight to Lweight. Added quadratic interpolation weights Qweight, which are used for conservative 9-points interpolation. The nesting connectivity NetCDF file has now the following variables:

  dimensions:
          ...
          nLweights = 4 ;
          nQweights = 9 ;
          ...
  variables:
          ...
          double Lweight(datum, nLweights) ;
                  Lweight:long_name = "horizontal linear interpolation weights" ;
                  Lweight:coordinates = "Xrg Yrg" ;
          double Qweight(datum, nQweights) ;
                  Qweight:long_name = "horizontal quadratic interpolation weights" ;
                  Qweight:coordinates = "Xrg Yrg" ;
          ...
  

• grid/contact.m: added the computation of quadratic interpolation weights, Qweight. They are computed in local function quadratic_weights.

The nine-points, quadratic, conservative interpolation following Clark and Farley (1984). In two-way nesting with refinement grids, the interpolation formula used to derive the contact point data from the donor coarse grid for use in the receiver finer grid (coarse to fine) needs to be consistent with the operators used to average the finer grid data down to the coarser grid resolution in its conjugate contact region (fine to coarse). This is referred as the conservation or reversibility condition. Clark and Farlay (1984) showed that the conservation condition reduced significantly the high frequency noise in the contact regions between coarser and finer grids when compared with nonreversible interpolation formulas. This reversibility condition is fundamental when coding the tangent linear and adjoint version of ROMS nesting routines.

The stencil diagram is:

                index7           index8           index9
 
                  Rm               Ro               Rp
 
             Jd+1 7________________8________________9     Sp
                  |                |                |
                  |                |                |
                  |                |<---p--->       |
                  |                |                |
                  |             Jr |........x   :   |
                  |                |       .    :   |
                  |                |       .    q   |
                  |                |       .    :   |
             Jd   |________________|____________:___|
                  4                5       Ir       6     So
                  |                |                |
                index4           index5           index6
                  |                |                |
                  |                |                |
                  |                |                |
                  |                |                |
                  |                |                |
             Jd-1 |________________|________________|
                  1                2                3     Sm
 
                 Id-1              Id              Id+1  
 
                index1           index2           index3

For quadratic interpolation the coefficients are:

      Rm = 0.5 * p * (p - 1) + alpha
      Ro = (1 - p^2)         - 2 * alpha
      Rp = 0.5 * p * (p + 1) + alpha

      Sm = 0.5 * q * (q - 1) + alpha
      So = (1 - q^2)         - 2 * alpha
      Sp = 0.5 * q * (q + 1) + alpha
  
where      alpha = [(1/rfactor)^2 - 1] / 24       for rfactor > 0
                                                  (conservative)

           alpha = 0                              for rfactor = 0

The quadratic interpolation weights for all water points are:

      W(1,:) = Rm * Sm
      W(2,:) = Ro * Sm
      W(3,:) = Rp * Sm
      W(4,:) = Rm * So
      W(5,:) = Ro * So
      W(6,:) = Rp * So
      W(7,:) = Rm * Sp
      W(8,:) = Ro * Sp
      W(9,:) = Rp * Sp

In ROMS, the quadratic interpolation is carried out as:

      F(Irg, Jrg) = W(1,:) * C(Idg-1, Jdg-1) +
                    W(2,:) * C(Idg  , Jdg-1) +
                    W(3,:) * C(Idg+1, Jdg-1) +
                    W(4,:) * C(Idg-1, Jdg  ) +
                    W(5,:) * C(Idg  , Jdg  ) +
                    W(6,:) * C(Idg+1, Jdg  ) +
                    W(7,:) * C(Idg-1, Jdg+1) +
                    W(8,:) * C(Idg  , Jdg+1) +
                    W(9,:) * C(Idg+1, Jdg+1)

• grid/grid_connections.m: update connectivity logic for more than two layers of refinement grids. It computes the refinement grid averaged area to help selection the appropriate nested grid contact region.

• grid/grid_perimeter.m: updated nested grid connectivity structure to include new linear and quadratic interpolation weights parameters.

• grid/grids_structure.m: added logic for processing history file which is used to compute vertical depth arrays.

• grid/plot_contact.m: update corner logical switch names from nested grid connectivity structure.

• grid/read_contact.m: updated reading of linear and quadratic interpolations weights from nested grid connectivity NetCDF file.

• grid/sponge.m: given a grid application NetCDF file, this function computes enhanced viscosity and diffusion scaling variables (visc_factor and diff_factor) that can be added input ROMS Grid NetCDF file. These scales are used in an application to set sponge areas with larger horizontal mixing coefficients for the damping of high frequency noise coming from open boundary conditions or nesting.

• grid/write_contact.m: updated writing of linear and quadratic interpolations weights into nested grid connectivity NetCDF file.

• initial/d_mercator2roms.m: corrected writing of spherical switch into NetCDF file.

• landmask/read_mask.m: updated reading of spherical switch from grid NetCDF file.

• seagrid/seagrid2roms.m: corrected typo when clipping bathymetry to maximum depth.

• utility/hslope.m: corrected reading of spherical switch from NetCDF file.

• utility/plot_nesting.m: updated computation of nested grids structure array.

• utility/sample_grid.m: corrected bug when computing the Jstr parameter. Many thanks to John Wilkin for bringing this to my attention. We need to have instead:

   Jstr = min(J(:))-offset;
   if (isnan(Jstr) || Jstr < 1),
     Jstr = 1;
   end