| Description |
- The Matlab script contact.m was updated to compute the interpolation weights correctly when there is land/sea masking points in the contact regions. A new local function, interp_weights, is added to compute the weights at the contact points within the donor grid cell:
index4 index3
4______________3 (Idg+1,Jdg+1)
| : |
| : p |
| 1-q <----->|
| : |
Jr|..:....x : |
| . : |
| 1-p . : |
|<------> q |
| . : |
| . : |
|___________:__|
(Idg,Jdg) 1 Ir 2
index1 index2
For linear interpolation and all water points, the weights are:
W(1,:) = (1 - p) * (1 - q)
W(2,:) = p * (1 - q)
W(3,:) = p * q
W(4,:) = (1 - p) * q
If any of the donor grid points is masked, the associated weight is
set to zero and the other weights are re-scaled accordingly.
The following figure shows a two-way nesting application in the South China Sea with a 1:3 refinement grid over the Luzon Straits. Both grids are rotated. This is a tidally forced application. Tides are imposed at the coarser, large scale grid open boundaries.The figure shows both coarser and finer grids free-surface fields after 10 days of solution. Notice the train of internal waves generated at Luzon Straits. It is very hard to see where is the finer grid in this application so I have to draw the perimeters (red outline). The two-way nesting is perfect between the coarser and finer grids... We are very happy with this configuration and solution. Many thanks to John Wilkin for his help with this application.
WARNING: You will need to regenerate the contact points NetCDF file using this updated version of contact.m if the donor grid has land/sea masking points in the contact region!!! This is the case in the above application. There is land/sea making in the southern and northern boundaries of the Luzon Straits refinement grid.
I also added few Matlab scripts:
- matlab/utility/pcolorjw.m: modified version of PCOLOR that expands the dimension the inputs X, Y and C so that the checkerboard will show all entries in C. The default PCOLOR clips the last row and column of C unless the shading('interp') is set, and the color of each square is determined by the value of C at the squares lower left corner. PCOLORJW recomputes X,Y at the mid points of the input, and pads C so that the effect is to shift each square one half space to the lower left. The perimeter squares are only half the width/height of all the others as can be seen if one sets shading('faceted'). Many thanks to John Wilkin for providing this modified function.
- matlab/utility/plot_nesting.m: plots requested ROMS nesting variable from input history NetCDF files. This function is very useful when debugging a ROMS nesting application. The plotting is not that fancy but it provides enough information for browsing ROMS variables very quickly.
- matlab/utility/plot_perimeter.m: adds a grid perimeter outline to an existing figure plotted with plot_nesting.m. It is used for browsing quickly output ROMS variables during nesting.
The above plot was done as follows:
Hnames={'scs_his_0010.nc', 'luzon_his_0010.nc'};
G=grids_structure(Hnames);
F=plot_nesting(G, Hnames, 'zeta', 8, 40, 1, [-2 4]);
shading flat; colormap(zebra);
h=plot_perimeter(G, 'r-');
print -r300 -dpng scs_ssh_day10.png
|
| Description |
nesting.F lines 850-855 has a call using DU_avg1 but it should DV_avg1:
CALL get_contact2d (dg, model, tile, &
& v2dvar, 'DV_avg1', &
& cr, Vcontact(cr)%Npoints, Vcontact, &
& LBi, UBi, LBj, UBj, &
& COUPLING(dg) % DV_avg1, &
& COMPOSITE(cr) % DV_avg1)
|
| Description |
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:
- 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.
|
