Difference between revisions of "Nested Grids"

From WikiROMS
Jump to navigationJump to search
m (Text replacement - "ocean.in" to "roms.in")   (change visibility)
(38 intermediate revisions by 2 users not shown)
Line 8: Line 8:
==Introduction==
==Introduction==


ROMS has a nested design structure for the temporal coupling between a hierarchy of donor and receiver grids.  All the state model variables are dynamically allocated and passed as arguments to the computational routines via Fortran-90 dereferenced pointer structures. All private arrays are automatic; their size is determined when the procedure is entered.  This code structure facilitates computations over nested grids of different topologies. There are three types of nesting capabilities in ROMS:
ROMS has a nested design structure for the temporal coupling between a hierarchy of '''donor''' and '''receiver''' grids.  All the state model variables are dynamically allocated and passed as arguments to the computational routines via Fortran-90 dereferenced pointer structures. All private arrays are automatic; their size is determined when the procedure is entered.  This code structure facilitates computations over nested grids of different topologies. There are three types of nesting capabilities in ROMS:


* '''Mosaics Grids''' which connect several grids along their edges,
* '''Mosaics Grids''' which connect several grids along their edges,
* '''Composite Grids''' which allow overlap regions of aligned and non-aligned grids, and
* '''Composite Grids''' which allow overlap regions of aligned and non-aligned grids, and
* '''Refinement Grids''' which provide increased resolution (3:1, 5:1, or 7:1) in a specific region.
* '''Refinement Grids''' which provide increased resolution ('''3:1''', '''5:1''', or '''7:1''') in a specific region.


The mosaic and composite grid code infrastructures are identical.  The differences are geometrical and primary based on the alignment between adjacent grids.  All the mosaic grids are exactly aligned with the adjacent grid.  In general, the mosaic grids are special case of the composite grids.
The mosaic and composite grid code infrastructures are identical.  The differences are geometric and primarily based on the alignment between adjacent grids.  All the mosaic grids are exactly aligned with the adjacent grid.  In general, the mosaic grids are a special case of the composite grids.


The nested grid implementation was divided in three sequential phases:
The nested grid implementation was divided in three sequential phases:


# Phase '''I''', released on April 2011 as ROMS version '''3.5''', included major updates to the nonlinear (NLM), perturbation tangent linear (TLM), finite-amplitude tangent linear (RPM), and adjoint (ADM) models numerical kernels to allow different horizontal '''I'''- and '''J'''-ranges in the '''DO'''-loops to permit operations on various nested grid classes (refinement, mosaics, and composite) and nesting layers (refinement and composite grid combinations).  This facilitates the computation of any horizontal operator (advection, diffusion, gradient, etc.) in the nesting contact regions and avoids the need for cumbersome lateral boundary conditions on the model variables and their associated flux/gradient values.  The advantage of this approach is that it is generic to any discrete horizontal operator. The contact region is an extended section of the grid that overlays an adjacent grid. The strategy is to compute the full horizontal operator at the contact points between nested grids instead of specifying boundary conditions.
# Phase '''I''', released in April 2011 as ROMS version '''3.5''', included major updates to the nonlinear (NLM), perturbation tangent linear (TLM), finite-amplitude tangent linear (RPM), and adjoint (ADM) model numerical kernels to allow different horizontal '''I'''- and '''J'''-ranges in the '''DO'''-loops to permit operations on various nested grid classes (refinement, mosaics, and composite) and nesting layers (refinement and composite grid combinations).  This facilitates the computation of any horizontal operator (advection, diffusion, gradient, etc.) in the nesting '''contact regions''' and avoids the need for cumbersome lateral boundary conditions on the model variables and their associated flux/gradient values.  The advantage of this approach is that it is generic to any discrete horizontal operator. The '''contact region''' is an extended section of the grid that overlays an adjacent grid. The strategy is to compute the full horizontal operator at the '''contact points''' between nested grids instead of specifying boundary conditions.
# Phase '''II''', released on September 2011 as ROMS version '''3.6''', included a major overhaul of lateral boundary conditions.  All their associated C-preprocessing options were eliminated and replaced with logical switches that depend on the nested grid, if any.  This facilitates, in a generic way, the processing or not of lateral boundary conditions in applications with nested grids.  In nesting applications, the values at the lateral boundary points are computed directly in the contact region by the numerical kernel.  Partial lateral boundary conditions are allowed by combining specified conditions (like closed, gradient, and others) with nested contact points. In addition, the logical switches allow different lateral boundary conditions types between active (temperature and salinity) and passive (biology, sediment, inert, etc.) tracers. The lateral boundary condition switches for each state variable and boundary edge are now specified in ROMS input script file, ocean.in.
# Phase '''II''', released in September 2011 as ROMS version '''3.6''', included a major overhaul of lateral boundary conditions.  All their associated C-preprocessing options were eliminated and replaced with logical switches that depend on the nested grid, if any.  This facilitates, in a generic way, the processing or not of lateral boundary conditions in applications with nested grids.  In nesting applications, the values at the lateral boundary points are computed directly in the '''contact region''' by the numerical kernel.  Partial lateral boundary conditions are allowed by combining specified conditions (like closed, gradient, and others) with nested '''contact points'''. In addition, the logical switches allow different lateral boundary conditions types between active (temperature and salinity) and passive (biology, sediment, inert, etc.) tracers. The lateral boundary condition switches for each state variable and boundary edge are now specified in ROMS input script file, roms.in.
# Phase '''III''', to be released soon as ROMS version '''3.7''', includes the nesting calls to ROMS main time-stepping routines, '''main2d''' and '''main3d'''.  The routine '''main2d''' is only used in shallow-water 2D barotropic applications whereas '''main3d''' is used in full 3D baroclinic applications.  Several new routines will be added to process the information that it is required in the overlap region, what information needs to be exchanged from/to another grid, and when to exchange it.  The exchange between grids can be one-way or two-way.
# Phase '''III''', to be released soon as ROMS version '''3.7''', includes the nesting calls to ROMS main time-stepping routines, '''main2d''' and '''main3d'''.  The routine '''main2d''' is only used in shallow-water 2D barotropic applications whereas '''main3d''' is used in full 3D baroclinic applications.  Several new routines will be added to process the information that is required in the '''contact region''', what information needs to be exchanged from/to another grid, and when to exchange it.  <br \> <br \>In grid '''refinement''' applications, the information is exchanged at the beginning of the full time-step (top of '''main2d''' or '''main3d''').  Contrarily, in '''composite''' grid applications, the information is exchanged between each sub-time step call in '''main2d''' or '''main3d'''.  That is, the composite '''donor''' and '''receiver''' grids need to sub-time step the 2D momentum equations before any of them start solving and coupling the 3D momentum equations.  Therefore, the governing equations are solved and nested in a synchronous fashion.  The concept of '''nesting layers''' is introduced to allow applications with both '''composite''' and '''refinement''' grid combinations, as shown in the diagram below for the Refinement and Partial Boundary Composite Sub-Class.  Here, two nesting layers are required.  In the first nesting layer, the composite grids '''1''' and '''2''' are synchronously sub-time stepped for each governing equation term.  Then, in the second and last nesting layer, the refinement grid '''3''' is time stepped and the information between '''donor''' ('''1''') and '''receiver''' refinement ('''3''') grids are exchanged at the beginning of the time step for both '''contact regions''': coarse-to-fine and fine-to-coarse.  The exchange between grids is two-way.


The current plan is to move to ROMS version '''4.0''' when all complex new developments are fully stable and tested.
The current plan is to move to ROMS version '''4.0''' when all complex new developments are fully stable and tested.  The design and implementation is fully adjointable and we have big plans in the future to support the adjoint-based algorithms, like 4D-Var, over selected or all nested grids.


==Nested Grids Classes==
==Nested Grids Classes==
<wikitex>
 
The ROMS nested grid design includes three Super-Classes and several Sub-Classes:
The ROMS nested grid design includes three Super-Classes and several Sub-Classes:


Line 38: Line 38:
# Composite and Refinement Combination Super-Class:
# Composite and Refinement Combination Super-Class:
#* Refinement and Partial Boundary Composite Sub-Class
#* Refinement and Partial Boundary Composite Sub-Class
#* Complex Estuary Refinement Grids Sub-Class
#* Complex Estuary Refinement-Composite Sub-Class
 
Hence, there are several possibilities and combinations. The design is flexible enough to allow complex nested grid configurations in coastal applications.  These capabilities are better illustrated with diagrams of the different Sub-Classes. The diagrams are all in terms of  computational coordinates (<math>\xi</math>, <math>\eta</math>) or fractional (<math>i</math>, <math>j</math>) coordinates.  Therefore we have squares and rectangles that may map to physical curvilinear coordinates. Here, '''Ngrids''' is the number of nested grids, '''NestLayers''' is the number of nested layers, '''GridsInLayer''' is a vector with the number of grids in each nested layer, and '''Ncontact''' is the number of nesting contact regions, ('''Ngrids'''-1)*2.


Hence, there are several possibilities and combinations. The design is flexible enough to allow complex nested grid configurations in coastal applications.  These capabilities are better illustrated with diagrams of the different Sub-Classes. The diagrams are all in terms of  computational coordinates ($\xi$, $\eta$) or fractional ($i$, $j$) coordinates.  Therefore we have squares and rectangles that may map to physical curvilinear coordinates.


<gallery widths=375px heights=200px perrow=2 caption="Composite Grids">
<gallery widths=375px heights=200px perrow=2 caption="Composite Grids">
Image:mosaic_nd.png|<b>Mosaic Sub-Class:</b><br />Ngrids = 2<br />NestLayers = 1<br />GridsInLayer = 1<br />Ncontact = 2
Image:mosaic_nd.png|<b>Mosaic Sub-Class:</b><br />Ngrids = 2<br />NestLayers = 1<br />GridsInLayer = 2<br />Ncontact = 2
Image:composite_overlap_nd.png|<b>Composite Overlap Sub-Class:</b><br />Ngrids = 2<br />NestLayers = 1<br />GridsInLayer = 1<br />Ncontact = 2
Image:composite_overlap_nd.png|<b>Composite Overlap Sub-Class:</b><br />Ngrids = 2<br />NestLayers = 1<br />GridsInLayer = 2<br />Ncontact = 2
Image:composite_nd.png|<b>Partial Boundary Composite Sub-Class:</b><br />Ngrids = 2<br />NestLayers = 1<br />GridsInLayer = 1<br />Ncontact = 2
Image:composite_nd.png|<b>Partial Boundary Composite Sub-Class:</b><br />Ngrids = 2<br />NestLayers = 1<br />GridsInLayer = 2<br />Ncontact = 2
Image:composite_estuary_nd.png|<b>Complex Estuary Composite Sub-Class:</b><br />Ngrids = 6<br />NestLayers = 1<br />GridsInLayer = 6<br />Ncontact = 10
Image:composite_estuary_nd.png|<b>Complex Estuary Composite Sub-Class:</b><br />Ngrids = 6<br />NestLayers = 1<br />GridsInLayer = 6<br />Ncontact = 10
</gallery>
</gallery>
Line 58: Line 59:
<gallery widths=375px heights=200px perrow=2 caption="Composite and Refinement Combination Grids">
<gallery widths=375px heights=200px perrow=2 caption="Composite and Refinement Combination Grids">
Image:composite-refinement_nd.png|<b>Refinement and Partial Boundary Composite Sub-Class:</b><br />Ngrids = 3<br />NestLayers = 2<br />GridsInLayer = 2 &nbsp; 1<br />Ncontact = 4
Image:composite-refinement_nd.png|<b>Refinement and Partial Boundary Composite Sub-Class:</b><br />Ngrids = 3<br />NestLayers = 2<br />GridsInLayer = 2 &nbsp; 1<br />Ncontact = 4
Image:complex_estuary_refinement_nd.png|<b>Complex Estuary Refinement Sub-Class:</b><br />Ngrids = 3<br />NestLayers = 2<br />GridsInLayer = 1 &nbsp; 2<br />Ncontact = 4
Image:complex_estuary_refinement_nd.png|<b>Complex Estuary Refinement-Composite Sub-Class:</b><br />Ngrids = 3<br />NestLayers = 2<br />GridsInLayer = 1 &nbsp; 2<br />Ncontact = 4
</gallery>
</gallery>


</wikitex>
 
{{warning}} '''Warning 1:'''  The order of nested grids in ROMS is '''extremely important'''.  We need to follow such order when specifying parameters in all ROMS standard input scripts: [[roms.in]], [[biology.in]], [[sediment.in]], and so on.  The numeric order is shown above for all nested grid Sub-Classes. It depends on a particular configuration topology, number of nested layers, and on the physical circulation.  Generally, the larger/coarser grid is first in the order list. Therefore, users need to order their grids wisely...
 
 
{{warning}} '''Warning 2:'''  In '''composite''' grids, it is assumed that the '''donor''' and the '''receiver''' grids have the same time-step size.  This is done to avoid the time interpolation between '''donor''' and '''receiver''' grids.  Only spatial interpolations are possible and carried out. This is due to the fact that '''composite''' grids need to solve the governing equations in a synchronous fashion, as explained above.  In '''refinement''' grids, it is assumed that the '''donor''' and '''receiver''' grids are an integer factor ('''RefineScale''') of the grid-size and time-step size. Therefore, time interpolation and time integration is possible between '''donor''' and '''receiver''' grids. 
 
 
{{warning}} '''Warning 3:''' In '''refinement''' grid applications, we do not recommend using refinement factors exceeding the '''1:7''' ratio. Most users should consider either the '''1:3''' or '''1:5''' ratio between '''donor''' and '''receiver''' refinement grid applications.  Notice that the refinement factor ('''RefineScale''') in C-grids needs to have an integer value given by '''2n + 1''' for '''n=1,2,3,...''' to always guarantee a '''&rho;'''-point at the center of the grid cell.
 
 
{{note}} '''Usually''', finer resolution is required to resolve estuary regions or other coastal features in a particular application.  In this case, we should consider a '''refinement''' grid with Land/Sea masking in such areas. This case is illustrated in the Complex Estuary Refinement-Composite Sub-Class diagram above.  Notice that the refinement grid '''2''' has an adjacent composite grid '''3''' of similar spatial resolution, so all the grids in the second nested layer have the same time step.  As you can see, ROMS nesting capabilities are powerful but require knowledge and expertise.  We highly recommend that users acquire this knowledge before configuring their own nested grid applications.  This knowledge is difficult to explain in words, but easier to show with the diagrams in this page.  We need to gain experience with this capability.


==Contact Regions and Contact Points==
==Contact Regions and Contact Points==
<wikitex>
 
A '''contact region''' is an extended section of the grid that overlays or is adjacent to a nested grid. A '''contact point''' is a grid cell inside the '''contact region'''. The '''contact points''' are located according to the staggered Arakawa C-grid horizontal positions.  Therefore, we have '''contact points''' at the &rho;-, &psi;-, u-, and v-points. However, the &psi;-points are only used to define the physical grid perimeters within a '''contact region'''. ROMS does not have state variables located at &psi;-points; vorticity is not explicitly specified in the governing equations. In a nested application, there are ([[Variables#Ngrids|Ngrids]]-1)*2 '''contact regions'''.  Each '''contact region''' has a '''receiver''' grid and '''donor''' grid.  The '''contact points''' of the '''receiver''' grid are processed using the '''donor''' grid cell containing the '''contact point'''.  The '''contact points''' are coincident if the '''receiver''' and '''donor''' grids occupy the same position, '''p = q = 0''' in diagram below. If coincident grids, the '''receiver''' grid '''contact point''' data is just filled using the '''donor''' data.  Otherwise, if not coincident, the '''receiver''' grid data is linearly interpolated from the '''donor''' grid cell containing the '''contact point'''.
A '''contact region''' is an extended section of the grid that overlays or is adjacent to a nested grid. A '''contact point''' is a grid cell inside the '''contact region'''. The '''contact points''' are located according to the staggered Arakawa C-grid horizontal positions.  Therefore, we have '''contact points''' at the '''&rho;'''-, '''&psi;'''-, '''u'''-, and '''v'''-points. However, the '''&psi;'''-points are only used to define the physical grid perimeters within a '''contact region'''. ROMS does not have state variables located at '''&psi;'''-points; vorticity is not explicitly specified in the governing equations. In a nested application, there are ([[Variables#Ngrids|Ngrids]]-1)*2 '''contact regions'''.  Each '''contact region''' has a '''receiver''' grid and '''donor''' grid.  The '''contact points''' of the '''receiver''' grid are processed using the '''donor''' grid cell containing the '''contact point'''.  The '''contact points''' are coincident if the '''receiver''' and '''donor''' grids occupy the same position, '''p = q = 0''' in diagram below. If coincident grids, the '''receiver''' grid '''contact point''' data is just filled using the '''donor''' data.  Otherwise (if not coincident) the '''receiver''' grid data is linearly interpolated from the '''donor''' grid cell containing the '''contact point'''.




Line 77: Line 88:
::::<span class="red">Hweight(4,:)</span> '''= (1.0 - p) * q'''
::::<span class="red">Hweight(4,:)</span> '''= (1.0 - p) * q'''


These horizontal interpolation weights are computed outside ROMS and read from '''contact points''' NetCDF file. Then, the contact point for a 2D state variable in the '''receiver''' grid is computed from the '''donor''' grid as:
These horizontal interpolation weights are computed outside ROMS and read from the '''contact points''' NetCDF file. Then, the contact point for a 2D state variable in the '''receiver''' grid is computed from the '''donor''' grid as:


::::<span class="red">V(Irg, Jrg)</span>  '''=''' <span class="red">Hweight(1,:)</span> '''*''' <span class="blue">F2D(Idg, Jdg)</span> +
::::<span class="red">V(Irg, Jrg)</span>  '''=''' <span class="red">Hweight(1,:)</span> '''*''' <span class="blue">F2D(Idg, Jdg)</span> +
Line 96: Line 107:




The contact points are processed outside of ROMS and all the connectivity within nested grids is read from an input NetCDF file.  This facilitates the configuration of various grid Sub-Classes. It simplifies tremendously the processing of such points in parallel computations.  It also give the user a full editing control of the contact points.  The contact points just to be set-up once for a particular application.
Notice the '''duality''' of a ROMS nested grid: data '''donor''' in one '''contact region''' and data '''receiver''' in its conjugate '''contact region'''. The exchange of information is always two-way.  This explains why we prefer to use the '''donor''' and '''receiver''' categories instead of the usual '''parent''' and '''child''' descriptions found in the literature. 
 
The contact points are processed outside of ROMS and all the connectivity within nested grids is read from an input NetCDF file.  This facilitates the configuration of various grid Sub-Classes. It tremendously simplifies the processing of such points in parallel computations.  It also gives the user full editing control of the contact points.  The contact points need to be set-up once for a particular application.


In ROMS, the contact region and contact point are managed with the following Fortran-90 structure in [[Modules#mod_nesting.F|mod_nesting.F]]:
In ROMS, the contact regions and contact points are managed with the following Fortran-90 structure in [[Modules#mod_nesting.F|mod_nesting.F]]:


   <span class="twilightBlue">    integer :: Ncontact                            ! total number of contact regions <br />
   <span class="twilightBlue">    integer :: Ncontact                            ! total number of contact regions <br />
Line 119: Line 132:
       TYPE (T_NGC), allocatable :: Vcontact(:)        ! contact points at V-points,  [Ncontact]</span>
       TYPE (T_NGC), allocatable :: Vcontact(:)        ! contact points at V-points,  [Ncontact]</span>


The diagrams below show the '''contact regions''' and '''contact points''' for each nesting Sub-Class.  Each diagram shows the '''contact region''' and its '''contact points''' colored '''purple''' as a transition between the '''blue''' and '''red''' grids.  Color is used to illustrate the nesting role between '''donor''' and '''receiver''' grids.  It also shows the number of nesting layers ('''NestLayer''') and number grids in each nested layer ('''GridsInLayer''').  Notice that there is an extra '''contact point''' to the '''left''' and '''bottom''' of the nested grid.  This is due to ROMS index enumeration system for staggered, Arakawa C-grids.  A '''left-bottom''' enumeration system is used: the '''i'''-index for '''u'''-points to the '''left''' of the '''&rho;'''-point have the same value; the '''j'''-index for '''v'''-points to the '''bottom''' of the '''&rho;'''-point have the same value.
The diagrams below show the '''contact regions''' and '''contact points''' for each nesting Sub-Class.  Each diagram shows the '''contact region''' and its '''contact points''' colored '''purple''' as a transition between the '''blue''' and '''red''' grids.  Color is used to illustrate the nesting role between '''donor''' and '''receiver''' grids.  You may click on the figure and the <span class="cyan">Full resolution</span> link to see a larger and better resolved diagram.  It also shows the number of nesting layers ('''NestLayer''') and number grids in each nested layer ('''GridsInLayer''').  Notice that there is an extra '''contact point''' to the '''left''' and '''bottom''' of the nested grid.  This is due to ROMS index enumeration system for staggered, Arakawa C-grids.  A '''left-bottom''' enumeration system is used: the '''i'''-index for '''u'''-points to the '''left''' of the '''&rho;'''-point have the same value; the '''j'''-index for '''v'''-points to the '''bottom''' of the '''&rho;'''-point have the same value.




Line 141: Line 154:


<gallery widths=375px heights=265px perrow=2 caption="Complex Estuary Composite Sub-Class">
<gallery widths=375px heights=265px perrow=2 caption="Complex Estuary Composite Sub-Class">
Image:composite_estuary_b_donor.png|Ngrids = 2<br />NestLayers = 1<br />GridsInLayer = 2<br />Ncontact = 2<br />Donor Grid = <span class="blue">blue</span><br />Receiver Grid = <span class="red">red</span>
Image:composite_estuary_b_donor.png|Ngrids = 6<br />NestLayers = 1<br />GridsInLayer = 6<br />Ncontact = 10<br />Donor Grid = <span class="blue">blue</span><br />Receiver Grid = <span class="red">red</span>
Image:composite_estuary_r_donor.png|Ngrids = 2<br />NestLayers = 1<br />GridsInLayer = 2<br />Ncontact = 2<br />Donor Grid = <span class="red">red</span><br />Receiver Grid = <span class="blue">blue</span>
Image:composite_estuary_r_donor.png|Ngrids = 6<br />NestLayers = 1<br />GridsInLayer = 6<br />Ncontact = 10<br />Donor Grid = <span class="red">red</span><br />Receiver Grid = <span class="blue">blue</span>
</gallery>
</gallery>


Line 165: Line 178:
</gallery>
</gallery>


<!--gallery widths=375px heights=285px perrow=2 caption="Complex Estuary Refinement Sub-Class (Coming Soon)">
Image:composite-refinement_b_donor.png|Ngrids = 3<br />NestLayers = 2<br />GridsInLayer = 2 &nbsp; 1<br />Ncontact = 4<br />Donor Grid = <span class="blue">blue</span><br />Receiver Grid = <span class="red">red</span>
Image:composite-refinement_r_donor.png|Ngrids = 3<br />NestLayers = 2<br />GridsInLayer = 2 &nbsp; 1<br />Ncontact = 4<br />Donor Grid = <span class="red">red</span><br />Receiver Grid = <span class="blue">blue</span>
</gallery-->
</wikitex>


==Metadata==
<gallery widths=375px heights=285px perrow=2 caption="Complex Estuary Refinement-Composite Sub-Class">
<wikitex>
Image:composite-refinement_estuary_b_donor.png|Ngrids = 3<br />NestLayers = 2<br />GridsInLayer = 1 &nbsp; 2<br />Ncontact = 4<br />Donor Grid = <span class="blue">blue</span><br />Receiver Grid = <span class="red">red</span>
The contact regions and contact points are processed outside ROMS and read in from a NetCDF file with the following Metadata structure:
Image:composite-refinement_estuary_r_donor.png|Ngrids = 3<br />NestLayers = 2<br />GridsInLayer = 1 &nbsp; 2<br />Ncontact = 4<br />Donor Grid = <span class="red">red</span><br />Receiver Grid = <span class="blue">blue</span>
Image:composite-refinement_estuary_comp_r_donor.png|Ngrids = 3<br />NestLayers = 2<br />GridsInLayer = 1 &nbsp; 2<br />Ncontact = 4<br />Donor Grid = <span class="red">red</span><br />Receiver Grid = <span class="green">green</span>
Image:composite-refinement_estuary_comp_g_donor.png|Ngrids = 3<br />NestLayers = 2<br />GridsInLayer = 1 &nbsp; 2<br />Ncontact = 4<br />Donor Grid = <span class="green">green</span><br />Receiver Grid = <span class="red">red</span>
</gallery>
 
 
==Nesting Test Cases==
 
<gallery widths=375px heights=275px perrow=2 caption="RHO-MASK">
Image:dogbone_rho-mask_whole_grid.png|<center>RHO-MASK Whole Grid</center>
Image:dogbone_rho-mask_comp_grids.png|<center>RHO-MASK Composite Grids</center>
</gallery>
 
 
<gallery widths=375px heights=275px perrow=2 caption="Contact Points, Composite">
Image:dogbone_comp_r_donor_noLmask.png|Donor Grid = <span class="blue">blue</span><br />Receiver Grid = <span class="red">red</span>
Image:dogbone_comp_b_donor_noLmask.png|Donor Grid = <span class="red">red</span><br />Receiver Grid = <span class="blue">blue</span>
</gallery>
 
 
<gallery widths=375px heights=275px perrow=2 caption="Contact Points, RHO-MASK, Composite">
Image:dogbone_comp_r_donor_Lmask.png|Donor Grid = <span class="blue">blue</span><br />Receiver Grid = <span class="red">red</span>
Image:dogbone_comp_b_donor_Lmask.png|Donor Grid = <span class="red">red</span><br />Receiver Grid = <span class="blue">blue</span>
</gallery>
 
 
<gallery widths=375px heights=275px perrow=3 caption="RHO-MASK">
Image:dogbone_rho-mask_whole_grid.png|<center>RHO-MASK Whole Grid</center>
Image:dogbone_rho-mask_ref_grid.png|<center>RHO-MASK Refinement Grid</center>
</gallery>
 
 
<gallery widths=375px heights=275px perrow=2 caption="Contact Points, Refinement">
Image:dogbone_ref_r_donor_noLmask.png|Donor Grid = <span class="red">red</span> (coarse)<br />Receiver Grid = <span class="blue">blue</span> (fine)
Image:dogbone_ref_b_donor_Lmask.png|Donor Grid = <span class="blue">blue</span> (fine)<br />Receiver Grid = <span class="red">red</span> (coarse)
</gallery>
 
 
<gallery widths=375px heights=275px perrow=2 caption="Contact Points, RHO-MASK, Refinement">
Image:dogbone_ref_r_donor_Lmask.png|Donor Grid = <span class="red">red</span> (coarse)<br />Receiver Grid = <span class="blue">blue</span> (fine)
Image:dogbone_ref_b_donor_Lmask.png|Donor Grid = <span class="blue">blue</span> (fine)<br />Receiver Grid = <span class="red">red</span> (coarse)
</gallery>
 
 
<gallery widths=375px heights=275px perrow=2 caption="Contact Points, Refinement">
Image:dogbone_ref_b_donor_u_avg.png|Donor Grid = <span class="blue">blue</span> (fine)<br />Receiver Grid = <span class="red">red</span> (coarse)
Image:dogbone_ref_b_donor_v_avg.png|Donor Grid = <span class="blue">blue</span> (fine)<br />Receiver Grid = <span class="red">red</span> (coarse)
</gallery>
 
 
<gallery widths=375px heights=275px perrow=2 caption="Contact Points, With (i,j) Donor Coordinates, Refinement">
Image:dogbone_ref_r_donor_noLmask_indicies.png|Donor Grid = <span class="blue">blue</span> (fine)<br />Receiver Grid = <span class="red">red</span> (coarse)
Image:dogbone_ref_r_donor_Lmask_indicies.png|Donor Grid = <span class="blue">blue</span> (fine)<br />Receiver Grid = <span class="red">red</span> (coarse)
</gallery>
 
 
 
 
==Realistic Nesting Configuration: US East Coast==
 
Several plots are shown below for a realistic nesting configuration in a US East Coast application using the Complex Estuary Refinement-Composite Sub-Class.  It illustrates ROMS nesting capabilities with complex coastlines and estuaries.  The coarser grid, '''ESPRESSO''' (130 x 82), has an average resolution of '''dx'''=7.5km, '''dy'''=5.8km.  The nested grids ('''&rho;'''-points mesh) are color coded for convenience to show the strategy used to better resolve the Delaware and Chesapeake Estuary Systems. The <span class="red">red</span> and <span class="green">green</span> are refinement grids whereas <span class="blue">blue</span> and <span class="violet">magenta</span> are composite grids. The refinement ratio is '''1:7'''.  An intermediate '''1:7''' refinement grid is created using Matlab script [[Grid_Processing_Scripts#coarse2fine|coarse2fine.m]] that included both the Delaware and Chesapeake Estuary Systems. Then,  the Matlab script [[Grid_Processing_Scripts#grid_extract|grid_extract.m]] is used to extract the Delaware Bay refinement grid (<span class="green">58 x 142</span>) and Delaware River composite grid (<span class="violet">42 x 55</span>).  Similarly, [[Grid_Processing_Scripts#grid_extract|grid_extract.m]] is used to extract the Chesapeake Bay outer refinement grid (<span class="red">135 x 142</span>) and Chesapeake Bay inner composite grid (<span class="blue">233 x 212</span>).
 
Similar strategy that combines both refinement and composite grids can be used to achieve the desired resolution in complex coastal regions. You need click on each figure and the <span class="cyan">Full resolution</span> link to see a larger and better resolved plot.
 
 
<gallery widths=375px heights=300px perrow=2 caption="ESPRESSO With Nested Delaware and Chesapeake Estuaries">
Image:ESPRESSO_nesting.png|<center>'''ESPRESSO''' domain with nested estuaries</center>
Image:Delaware_Chesapeake_zoom.png|<center>Closer look at Delaware and Chesapeake Estuary grids</center>
Image:Delaware_contact_zoom.png|<center>Contact region between the Delaware <span class="violet">River</span> and <span class="green">Estuary</span></center>
Image:Chesapeake_contact_zoom.png|<center>Contact region between <span class="blue">inner</span> and <span class="red">outer</span> Chesapeake Bay</center>
</gallery>
 
 
==Contact Points Metadata==
 
The '''contact regions''' and '''contact points''' are processed outside ROMS and read in from a NetCDF file with the following Metadata structure:


   <span class="twilightBlue"> dimensions: <br />
   <span class="twilightBlue"> dimensions: <br />
Line 179: Line 262:
         Ncontact = 2 ;
         Ncontact = 2 ;
         Nweights = 4 ;
         Nweights = 4 ;
         datum = 98 ; <br />
         datum = 350 ; <br />
   variables: <br />
   variables: <br />
         int spherical ;
         int spherical ;
Line 232: Line 315:
                 on_boundary:long_name = "contact point on receiver grid physical boundary" ;
                 on_boundary:long_name = "contact point on receiver grid physical boundary" ;
                 on_boundary:flag_values = 0, 1, 2, 3, 4 ;
                 on_boundary:flag_values = 0, 1, 2, 3, 4 ;
                 on_boundary:flag_meanings = "false western southern eastern northern" ;
                 on_boundary:flag_meanings = "other western southern eastern northern" ;
         int Idg(datum) ;
         int Idg(datum) ;
                 Idg:long_name = "I-left index of donor cell containing contact point" ;
                 Idg:long_name = "I-left index of donor cell containing contact point" ;
Line 299: Line 382:
                 :history = "Contact points created by contact.m and written by write_contact.m" ; </span>
                 :history = "Contact points created by contact.m and written by write_contact.m" ; </span>


Notice that the variable dimensions (<span class="red">Hweight</span>) are in '''raw-major''' order since this is the output from '''ncdump'''. However, the data is stored in '''column-major''' order because this file will be processed by ROMS in Fortran-90. The strategy here is to '''pack''' the contact points information for each contact region in a long vector.  Several indices are available to '''unpack''' the data into the Fortran-90 contact structures for &rho;-, u-, and v-points.  Although NetCDF-4 supports derived type structures via the HDF5 library, a simpler Metadata model is chosen for interoperability. Notice that the contact points NetCDF file contains also the values for the various grid metrics needed in ROMS numerical kernel.  Only the horizontal interpolation weights are provided here since they are always the same for a particular configuration.  The vertical interpolation weights are computed internally in ROMS at each time step because of the time dependency of the vertical coordinates to the evolution of the free-surface.
Notice that the variable dimensions (like <span class="red">Hweight</span>) are in '''row-major''' order since this is the output from '''ncdump'''. However, the data is stored in '''column-major''' order because this file will be processed by ROMS in Fortran-90. The strategy here is to '''pack''' the '''contact points''' information for each '''contact region''' in a long vector.  Several indices are available to '''unpack''' the data into the Fortran-90 '''contact point''' structures for '''&rho;'''-, '''u'''-, and '''v'''-points.  Although '''NetCDF-4''' supports derived type structures via the '''HDF5''' library, a simpler Metadata model is chosen for interoperability. Notice that the '''contact points''' NetCDF file also contains the values for the various grid metrics needed in ROMS numerical kernel.  Only the horizontal interpolation weights are provided here since they are always the same (time independent) for a particular configuration.  The vertical interpolation weights are computed internally in ROMS at each time step because of the time dependency of the vertical coordinates to the evolution of the free-surface.
 
</wikitex>

Revision as of 15:13, 17 July 2019

Nested Grids
Nested Grids Menu
  1. Introduction
  2. Matlab Processing Scripts

Introduction

ROMS has a nested design structure for the temporal coupling between a hierarchy of donor and receiver grids. All the state model variables are dynamically allocated and passed as arguments to the computational routines via Fortran-90 dereferenced pointer structures. All private arrays are automatic; their size is determined when the procedure is entered. This code structure facilitates computations over nested grids of different topologies. There are three types of nesting capabilities in ROMS:

  • Mosaics Grids which connect several grids along their edges,
  • Composite Grids which allow overlap regions of aligned and non-aligned grids, and
  • Refinement Grids which provide increased resolution (3:1, 5:1, or 7:1) in a specific region.

The mosaic and composite grid code infrastructures are identical. The differences are geometric and primarily based on the alignment between adjacent grids. All the mosaic grids are exactly aligned with the adjacent grid. In general, the mosaic grids are a special case of the composite grids.

The nested grid implementation was divided in three sequential phases:

  1. Phase I, released in April 2011 as ROMS version 3.5, included major updates to the nonlinear (NLM), perturbation tangent linear (TLM), finite-amplitude tangent linear (RPM), and adjoint (ADM) model numerical kernels to allow different horizontal I- and J-ranges in the DO-loops to permit operations on various nested grid classes (refinement, mosaics, and composite) and nesting layers (refinement and composite grid combinations). This facilitates the computation of any horizontal operator (advection, diffusion, gradient, etc.) in the nesting contact regions and avoids the need for cumbersome lateral boundary conditions on the model variables and their associated flux/gradient values. The advantage of this approach is that it is generic to any discrete horizontal operator. The contact region is an extended section of the grid that overlays an adjacent grid. The strategy is to compute the full horizontal operator at the contact points between nested grids instead of specifying boundary conditions.
  2. Phase II, released in September 2011 as ROMS version 3.6, included a major overhaul of lateral boundary conditions. All their associated C-preprocessing options were eliminated and replaced with logical switches that depend on the nested grid, if any. This facilitates, in a generic way, the processing or not of lateral boundary conditions in applications with nested grids. In nesting applications, the values at the lateral boundary points are computed directly in the contact region by the numerical kernel. Partial lateral boundary conditions are allowed by combining specified conditions (like closed, gradient, and others) with nested contact points. In addition, the logical switches allow different lateral boundary conditions types between active (temperature and salinity) and passive (biology, sediment, inert, etc.) tracers. The lateral boundary condition switches for each state variable and boundary edge are now specified in ROMS input script file, roms.in.
  3. Phase III, to be released soon as ROMS version 3.7, includes the nesting calls to ROMS main time-stepping routines, main2d and main3d. The routine main2d is only used in shallow-water 2D barotropic applications whereas main3d is used in full 3D baroclinic applications. Several new routines will be added to process the information that is required in the contact region, what information needs to be exchanged from/to another grid, and when to exchange it.

    In grid refinement applications, the information is exchanged at the beginning of the full time-step (top of main2d or main3d). Contrarily, in composite grid applications, the information is exchanged between each sub-time step call in main2d or main3d. That is, the composite donor and receiver grids need to sub-time step the 2D momentum equations before any of them start solving and coupling the 3D momentum equations. Therefore, the governing equations are solved and nested in a synchronous fashion. The concept of nesting layers is introduced to allow applications with both composite and refinement grid combinations, as shown in the diagram below for the Refinement and Partial Boundary Composite Sub-Class. Here, two nesting layers are required. In the first nesting layer, the composite grids 1 and 2 are synchronously sub-time stepped for each governing equation term. Then, in the second and last nesting layer, the refinement grid 3 is time stepped and the information between donor (1) and receiver refinement (3) grids are exchanged at the beginning of the time step for both contact regions: coarse-to-fine and fine-to-coarse. The exchange between grids is two-way.

The current plan is to move to ROMS version 4.0 when all complex new developments are fully stable and tested. The design and implementation is fully adjointable and we have big plans in the future to support the adjoint-based algorithms, like 4D-Var, over selected or all nested grids.

Nested Grids Classes

The ROMS nested grid design includes three Super-Classes and several Sub-Classes:

  1. Composite Grids Super-Class:
    • Mosaic Grids Sub-Class
    • Composite Overlap Grids Sub-Class
    • Partial Boundary Composite Grids Sub-Class
    • Complex Estuary Composite Grids Sub-Class
  2. Refinement Grids Super-Class:
    • Single Refinement Sub-Class
    • Multiple Refinement Sub-Class
  3. Composite and Refinement Combination Super-Class:
    • Refinement and Partial Boundary Composite Sub-Class
    • Complex Estuary Refinement-Composite Sub-Class

Hence, there are several possibilities and combinations. The design is flexible enough to allow complex nested grid configurations in coastal applications. These capabilities are better illustrated with diagrams of the different Sub-Classes. The diagrams are all in terms of computational coordinates (, ) or fractional (, ) coordinates. Therefore we have squares and rectangles that may map to physical curvilinear coordinates. Here, Ngrids is the number of nested grids, NestLayers is the number of nested layers, GridsInLayer is a vector with the number of grids in each nested layer, and Ncontact is the number of nesting contact regions, (Ngrids-1)*2.





Warning Warning 1: The order of nested grids in ROMS is extremely important. We need to follow such order when specifying parameters in all ROMS standard input scripts: roms.in, biology.in, sediment.in, and so on. The numeric order is shown above for all nested grid Sub-Classes. It depends on a particular configuration topology, number of nested layers, and on the physical circulation. Generally, the larger/coarser grid is first in the order list. Therefore, users need to order their grids wisely...


Warning Warning 2: In composite grids, it is assumed that the donor and the receiver grids have the same time-step size. This is done to avoid the time interpolation between donor and receiver grids. Only spatial interpolations are possible and carried out. This is due to the fact that composite grids need to solve the governing equations in a synchronous fashion, as explained above. In refinement grids, it is assumed that the donor and receiver grids are an integer factor (RefineScale) of the grid-size and time-step size. Therefore, time interpolation and time integration is possible between donor and receiver grids.


Warning Warning 3: In refinement grid applications, we do not recommend using refinement factors exceeding the 1:7 ratio. Most users should consider either the 1:3 or 1:5 ratio between donor and receiver refinement grid applications. Notice that the refinement factor (RefineScale) in C-grids needs to have an integer value given by 2n + 1 for n=1,2,3,... to always guarantee a ρ-point at the center of the grid cell.


Note Usually, finer resolution is required to resolve estuary regions or other coastal features in a particular application. In this case, we should consider a refinement grid with Land/Sea masking in such areas. This case is illustrated in the Complex Estuary Refinement-Composite Sub-Class diagram above. Notice that the refinement grid 2 has an adjacent composite grid 3 of similar spatial resolution, so all the grids in the second nested layer have the same time step. As you can see, ROMS nesting capabilities are powerful but require knowledge and expertise. We highly recommend that users acquire this knowledge before configuring their own nested grid applications. This knowledge is difficult to explain in words, but easier to show with the diagrams in this page. We need to gain experience with this capability.

Contact Regions and Contact Points

A contact region is an extended section of the grid that overlays or is adjacent to a nested grid. A contact point is a grid cell inside the contact region. The contact points are located according to the staggered Arakawa C-grid horizontal positions. Therefore, we have contact points at the ρ-, ψ-, u-, and v-points. However, the ψ-points are only used to define the physical grid perimeters within a contact region. ROMS does not have state variables located at ψ-points; vorticity is not explicitly specified in the governing equations. In a nested application, there are (Ngrids-1)*2 contact regions. Each contact region has a receiver grid and donor grid. The contact points of the receiver grid are processed using the donor grid cell containing the contact point. The contact points are coincident if the receiver and donor grids occupy the same position, p = q = 0 in diagram below. If coincident grids, the receiver grid contact point data is just filled using the donor data. Otherwise (if not coincident) the receiver grid data is linearly interpolated from the donor grid cell containing the contact point.


contact points interp flat.png   contact points interp cube.png

Here, the suffixes dg and rg are for the donor grid and receiver grid, respectively. The horizontal interpolation weights are as follows:

Hweight(1,:) = (1.0 - p) * (1.0 - q)
Hweight(2,:) = p * (1.0 - q)
Hweight(3,:) = p * q
Hweight(4,:) = (1.0 - p) * q

These horizontal interpolation weights are computed outside ROMS and read from the contact points NetCDF file. Then, the contact point for a 2D state variable in the receiver grid is computed from the donor grid as:

V(Irg, Jrg) = Hweight(1,:) * F2D(Idg, Jdg) +
Hweight(2,:) * F2D(Idg+1, Jdg) +
Hweight(3,:) * F2D(Idg+1, Jdg+1) +
Hweight(4,:) * F2D(Idg, Jdg+1)

If coincident contact points between data donor and data receiver grids, p = q = 0,

Hweight(1,:) = 1.0
Hweight(2,:) = 0.0
Hweight(3,:) = 0.0
Hweight(4,:) = 0.0

and

V(Irg, Jrg) = Hweight(1,:) * F2D(Idg, Jdg)


Notice the duality of a ROMS nested grid: data donor in one contact region and data receiver in its conjugate contact region. The exchange of information is always two-way. This explains why we prefer to use the donor and receiver categories instead of the usual parent and child descriptions found in the literature.

The contact points are processed outside of ROMS and all the connectivity within nested grids is read from an input NetCDF file. This facilitates the configuration of various grid Sub-Classes. It tremendously simplifies the processing of such points in parallel computations. It also gives the user full editing control of the contact points. The contact points need to be set-up once for a particular application.

In ROMS, the contact regions and contact points are managed with the following Fortran-90 structure in mod_nesting.F:

     integer :: Ncontact                             ! total number of contact regions 
TYPE T_NGC
logical :: coincident  ! coincident donor/receiver grid switch, p=q=0 logical :: interpolate  ! switch to perform vertical interpolation
integer :: donor_grid  ! data donor grid number integer :: receiver_grid  ! data receiver grid number
integer :: Npoints  ! number of points in contact region
integer, pointer :: Irg(:)  ! receiver grid, I-contact point integer, pointer :: Jrg(:)  ! receiver grid, J-contact point
integer, pointer :: Idg(:)  ! donor grid, contact cell I-left index integer, pointer :: Jdg(:)  ! donor grid, contact cell J-bottom index integer, pointer :: Kdg(:,:)  ! donor grid, contact cell K-index
real(r8), pointer :: Hweight(:,:)  ! horizontal interpolation weights real(r8), pointer :: Vweight(:,:,:)  ! vertical interpolation weights
END TYPE T_NGC
TYPE (T_NGC), allocatable :: Rcontact(:)  ! contact points at RHO-points, [Ncontact] TYPE (T_NGC), allocatable :: Ucontact(:)  ! contact points at U-points, [Ncontact] TYPE (T_NGC), allocatable :: Vcontact(:)  ! contact points at V-points, [Ncontact]

The diagrams below show the contact regions and contact points for each nesting Sub-Class. Each diagram shows the contact region and its contact points colored purple as a transition between the blue and red grids. Color is used to illustrate the nesting role between donor and receiver grids. You may click on the figure and the Full resolution link to see a larger and better resolved diagram. It also shows the number of nesting layers (NestLayer) and number grids in each nested layer (GridsInLayer). Notice that there is an extra contact point to the left and bottom of the nested grid. This is due to ROMS index enumeration system for staggered, Arakawa C-grids. A left-bottom enumeration system is used: the i-index for u-points to the left of the ρ-point have the same value; the j-index for v-points to the bottom of the ρ-point have the same value.










Nesting Test Cases










Realistic Nesting Configuration: US East Coast

Several plots are shown below for a realistic nesting configuration in a US East Coast application using the Complex Estuary Refinement-Composite Sub-Class. It illustrates ROMS nesting capabilities with complex coastlines and estuaries. The coarser grid, ESPRESSO (130 x 82), has an average resolution of dx=7.5km, dy=5.8km. The nested grids (ρ-points mesh) are color coded for convenience to show the strategy used to better resolve the Delaware and Chesapeake Estuary Systems. The red and green are refinement grids whereas blue and magenta are composite grids. The refinement ratio is 1:7. An intermediate 1:7 refinement grid is created using Matlab script coarse2fine.m that included both the Delaware and Chesapeake Estuary Systems. Then, the Matlab script grid_extract.m is used to extract the Delaware Bay refinement grid (58 x 142) and Delaware River composite grid (42 x 55). Similarly, grid_extract.m is used to extract the Chesapeake Bay outer refinement grid (135 x 142) and Chesapeake Bay inner composite grid (233 x 212).

Similar strategy that combines both refinement and composite grids can be used to achieve the desired resolution in complex coastal regions. You need click on each figure and the Full resolution link to see a larger and better resolved plot.



Contact Points Metadata

The contact regions and contact points are processed outside ROMS and read in from a NetCDF file with the following Metadata structure:

  dimensions: 
Ngrids = 2 ; Ncontact = 2 ; Nweights = 4 ; datum = 350 ;
variables:
int spherical ; spherical:long_name = "grid type logical switch" ; spherical:flag_values = 0, 1 ; spherical:flag_meanings = "Cartesian spherical" ; int Lm(Ngrids) ; Lm:long_name = "number of interior RHO-points in the I-direction" ; int Mm(Ngrids) ; Mm:long_name = "number of interior RHO-points in the J-direction" ; int coincident(Ngrids) ; coincident:long_name = "coincident donor and receiver grids logical switch" ; coincident:flag_values = 0, 1 ; coincident:flag_meanings = "false true" ; int composite(Ngrids) ; composite:long_name = "composite grid type logical switch" ; composite:flag_values = 0, 1 ; composite:flag_meanings = "false true" ; int mosaic(Ngrids) ; mosaic:long_name = "mosaic grid type logical switch" ; mosaic:flag_values = 0, 1 ; mosaic:flag_meanings = "false true" ; int refinement(Ngrids) ; refinement:long_name = "refinement grid type logical switch" ; refinement:flag_values = 0, 1 ; refinement:flag_meanings = "false true" ; int refine_factor(Ngrids) ; refine_factor:long_name = "refinement factor from donor grid" ; int interpolate(Ncontact) ; interpolate:long_name = "vertical interpolation at contact points logical switch" ; interpolate:flag_values = 0, 1 ; interpolate:flag_meanings = "false true" ; int donor_grid(Ncontact) ; donor_grid:long_name = "data donor grid number" ; int receiver_grid(Ncontact) ; receiver_grid:long_name = "data receiver grid number" ; int NstrR(Ncontact) ; NstrR:long_name = "starting contact RHO-point index in data vector" ; int NendR(Ncontact) ; NendR:long_name = "ending contact RHO-point index in data vector" ; int NstrU(Ncontact) ; NstrU:long_name = "starting contact U-point index in data vector" ; int NendU(Ncontact) ; NendU:long_name = "ending contact U-point index in data vector" ; int NstrV(Ncontact) ; NstrV:long_name = "starting contact V-point index in data vector" ; int NendV(Ncontact) ; NendV:long_name = "ending contact V-point index in data vector" ; int contact_region(datum) ; contact_region:long_name = "contact region number" ; int on_boundary(datum) ; on_boundary:long_name = "contact point on receiver grid physical boundary" ; on_boundary:flag_values = 0, 1, 2, 3, 4 ; on_boundary:flag_meanings = "other western southern eastern northern" ; int Idg(datum) ; Idg:long_name = "I-left index of donor cell containing contact point" ; int Jdg(datum) ; Jdg:long_name = "J-bottom index of donor cell containing contact point" ; int Irg(datum) ; Irg:long_name = "receiver grid I-index of contact point" ; Irg:coordinates = "Xrg Yrg" ; int Jrg(datum) ; Jrg:long_name = "receiver grid J-index of contact point" ; Jrg:coordinates = "Xrg Yrg" ; double Xrg(datum) ; Xrg:long_name = "X-location of receiver grid contact points" ; Xrg:units = "meter" ; double Yrg(datum) ; Yrg:long_name = "Y-location of receiver grid contact points" ; Yrg:units = "meter" ; double Hweight(datum, Nweights) ; Hweight:long_name = "horizontal interpolation weights" ; Hweight:coordinates = "Xrg Yrg" ; double h(datum) ; h:long_name = "bathymetry at RHO-points" ; h:units = "meter" ; h:coordinates = "Xrg Yrg" ; h:_FillValue = 1.e+37 ; double f(datum) ; f:long_name = "Coriolis parameter at RHO-points" ; f:units = "second-1" ; f:coordinates = "Xrg Yrg" ; f:_FillValue = 1.e+37 ; double pm(datum) ; pm:long_name = "curvilinear coordinate metric in XI" ; pm:units = "meter-1" ; pm:coordinates = "Xrg Yrg" ; pm:_FillValue = 1.e+37 ; double pn(datum) ; pn:long_name = "curvilinear coordinate metric in ETA" ; pn:units = "meter-1" ; pn:coordinates = "Xrg Yrg" ; pn:_FillValue = 1.e+37 ; double dndx(datum) ; dndx:long_name = "XI-derivative of inverse metric factor pn" ; dndx:units = "meter" ; dndx:coordinates = "Xrg Yrg" ; dndx:_FillValue = 1.e+37 ; double dmde(datum) ; dmde:long_name = "ETA-derivative of inverse metric factor pm" ; dmde:units = "meter" ; dmde:coordinates = "Xrg Yrg" ; dmde:_FillValue = 1.e+37 ; double angle(datum) ; angle:long_name = "angle between XI-axis and EAST" ; angle:units = "radians" ; angle:coordinates = "Xrg Yrg" ; angle:_FillValue = 1.e+37 ; double mask(datum) ; mask:long_name = "land-sea mask of contact points" ; mask:flag_values = 0., 1. ; mask:flag_meanings = "land water" ; mask:coordinates = "Xrg Yrg" ;
global attributes:
:type = "ROMS Nesting Contact Regions Data" ; :grid_files = " \n", "/Users/arango/ocean/repository/TestCases/dogbone/Data/dogbone_grd_left.nc \n", "/Users/arango/ocean/repository/TestCases/dogbone/Data/dogbone_grd_right.nc" ; :history = "Contact points created by contact.m and written by write_contact.m" ;

Notice that the variable dimensions (like Hweight) are in row-major order since this is the output from ncdump. However, the data is stored in column-major order because this file will be processed by ROMS in Fortran-90. The strategy here is to pack the contact points information for each contact region in a long vector. Several indices are available to unpack the data into the Fortran-90 contact point structures for ρ-, u-, and v-points. Although NetCDF-4 supports derived type structures via the HDF5 library, a simpler Metadata model is chosen for interoperability. Notice that the contact points NetCDF file also contains the values for the various grid metrics needed in ROMS numerical kernel. Only the horizontal interpolation weights are provided here since they are always the same (time independent) for a particular configuration. The vertical interpolation weights are computed internally in ROMS at each time step because of the time dependency of the vertical coordinates to the evolution of the free-surface.