How to generate and use composite grids?

[AlexisEspinosa]

[The critical question in this post is labelled as Q.3.b. If you want to skip all the introductory explanation, you can jump directly there.] The complete set of questions are labelled Q.1.a, Q.3.a, Q.3.b, Q.4.a and Q.5.a.

Hello,

We are trying to generate a realistic nesting configuration in ROMS for the simulation of some estuaries of our interest. Our aim is to study a case similar to what is shown for the 'US East Coast' example figures in the wiki: with a big domain, a nested grid just outside the estuary, and a composite grid for the estuary (in contact with the nested grid and with the same resolution):
https://www.myroms.org/wiki/index.php/Nested_Grids


1) But, unfortunately, I have not found any documented, or partially-documented example for a case like this. The only available source of information are the dogbone test example and the wiki page mentioned above, but this information is not enough for generating a complex ROMS case like this (or at least I have not been able to find the correct source of information).
Q.1.a) Does anyone knows a good source of information for generating this type of cases?

2) As a first step I decided generate my own replicas of the 'nesting' and 'composite' branches of the dogbone test example. After confirming that both original test examples branches ran properly, I REPLICATED SUCCESSFULLY THE NESTING CASE starting from the basic "whole" case files. In the replicating process I generated the refined nested grid with "coarse2fine.m", and generated the contact-point file with "c_contact.m" and "contact.m". The initial conditions were generated with a modified version of "d_roms2roms.m". The replicated case ran properly! So far so good, but ...

3) When trying to replicate the 'composite' branch of the dogbone test case, I'm having some errors while generating the contact-point file. First, I succeeded to generate proper replicas of the 'left' and 'right' grids from the original 'whole' grid. This was done using "grid_extract.m" twice dividing the original grid in two. But THE CONTACT-POINTS FILE GENERATION IS FAILING. It even fails when using the original grids.
Q.3.a) What is wrong with the generation of the contact-points file? My matlab-function usage is the following:

Code: Select all

c_contact('dogbone_ngc_AlexisComposite.nc',false,2); %AEG: creating the contact file
contact({'dogbone_grd_left.nc','dogbone_grd_right.nc'}, ...
         'dogbone_ngc_AlexisComposite.nc', ...
         false,false,true); %AEG: Filling the contact file

and the output I get is (with errors) is:

Code: Select all

Summary of Contact Regions Processed:
  
   Grid 01: D:\AlexisEspinosa\dogbone_grd_left.nc
   Grid 02: D:\AlexisEspinosa\dogbone_grd_right.nc
  
   Contact   Donor   Receiver
    Region    Grid       Grid
  
     01        01        02
     02        02        01
  
Error using sub2ind (line 52)
Out of range subscript.

Error in contact>Hweights (line 2453)
            RindexB = sub2ind([Ir, Jr],
            ...

Error in contact (line 352)
S = Hweights(G, S, ImposeMask);

Error in refineMasterScript_RecreateComposite (line 163)
contact({'dogbone_grd_left.c','dogbone_grd_right.nc'}, ...

The sub2ind matlab function is finding indexes out of range when generating RindexB within the function Hweights within "contact.m". Here is that part of the code. The code is failing in the line for RindexB and can't move forward to the rest.

Code: Select all

    [Ir,Jr] = size(S.grid(dg).I_rho);
    RindexB = sub2ind([Ir, Jr],                                         ...
                      S.contact(cr).point.Idg_rho+1,                    ...
                      min(Jr,S.contact(cr).point.Jdg_rho));  % (Idg, Jdg-1)
    RindexO = sub2ind([Ir, Jr],                                         ...
                      S.contact(cr).point.Idg_rho+1,                    ...
                      S.contact(cr).point.Jdg_rho+1);        % (Idg, Jdg)
    RindexT = sub2ind([Ir, Jr],                                         ...
                      S.contact(cr).point.Idg_rho+1,                    ...
                      min(Jr,S.contact(cr).point.Jdg_rho+2));% (Idg, Jdg+1)
    RindexR = sub2ind([Ir, Jr],                                         ...
                      min(Ir, S.contact(cr).point.Idg_rho+2),           ...
                      S.contact(cr).point.Jdg_rho+1);        % (Idg+1, Jdg)

I debugged the code (and indeed commented RindexB line for a test) and found that for RindexO, RindexT and RindexR the indexes do not go out of range, even when they are increased by 2. No out of range happens because of the usage of the "min" matlab function. But in the case of RindexB, the min does nothing because it still allows indexes to have "0" values (ZEROS) for the "Jdg-1" case. The zeros are, obviously, out of range for matlab arrays and therefore the error.

I tried a "blind" solution for the indexes problem by modifying the code with "max" functions to impede the appearance of zeros in the generation of RindexB and UindexB:

Code: Select all

    [Ir,Jr] = size(S.grid(dg).I_rho);
% %     RindexB = sub2ind([Ir, Jr],                                         ...
% %                       S.contact(cr).point.Idg_rho+1,                    ...
% %                       min(Jr,S.contact(cr).point.Jdg_rho));  % (Idg, Jdg-1)
    RindexB = sub2ind([Ir, Jr],                                         ...
                      S.contact(cr).point.Idg_rho+1,                    ...
                      max(1,S.contact(cr).point.Jdg_rho));  % (Idg, Jdg-1)

.
.
.
.


    [Iu,Ju] = size(S.grid(dg).I_u);
% %     UindexB = sub2ind([Iu, Ju],                                         ...
% %                       S.contact(cr).point.Idg_u  ,                      ...
% %                       min(Ju, S.contact(cr).point.Jdg_u));   % (Idg, Jdg-1)
    UindexB = sub2ind([Iu, Ju],                                         ...
                      S.contact(cr).point.Idg_u  ,                      ...
                      max(1, S.contact(cr).point.Jdg_u));   % (Idg, Jdg-1)

but additional errors appeared, so "blind" corrections are clearly not the way to go. So I preferred to ask the forum for the correct usage of this function for the dogbone case and for general (realistic) composite grids for estuary simulations.

Q.3.b) I have tried different combinations for the logical, and numerical parameters of "c_contact" and "contact" (obligatory and optional) but have not succeeded yet to generate the contact-points file for the 'composite' dogbone test case. Maybe I'm still doing something wrong (read Q.3.a to see my settings). Am I doing something wrong? or Is there a bug in the "contact.m" function? I wonder how you generated the "dogbone_ngc_composite.nc" file for the dogbone composite test case in the first place?

Q.4.a) If someone has a solution for Q.3.a and Q.3.b above, can you go further and explain how can I create a new dogbone case: 'nesting+composite dogbone' (similar configuration to the realistic case): one with a nested refined grid just out of the "bone neck" together with a composite refined grid for the rest of the other side of the bone? (I'm wishing to generate a dogbone case with a grid-arrangement similar to the one for the realistic estuary case mentioned above.)

Q.5.a) If someone has an answer for the questions above, can you go further and explain what else should be done for generating a realistic case with nested-composite grids for estuary simulations? Will the process be exactly the same as for the new 'nesting+composite dogbone' case proposed in Q.4.a?

Thanks a lot for your patience and interest reading this post. Any kind of help is very welcome (Please do not forge Q.1.a very above!).

Alexis Espinosa
UWA

[arango]

Okay, I need to track this problem in contact.m. I was able to reproduce you error. I will look at this soon to see what is the solution for this out-of-range values when computing the interpolation weights.

I also have in my to do list to have a realistic application with composite grids. However, this will take some time.

[AlexisEspinosa]

Thanks a lot Prof. Arango,

We are looking forward to the solution of this issue, as we are ready to incorporate further refinement in our studies of several coastal regions with estuaries!

I have just read again your latest post about the latest major update of Nesting algorithms. I understand that floats simulations are not possible right now among nested grids. That is not a big issue for us, although it would be wonderful to have it! But what about biological models?

Q.1) Are biological models ready to be used with coarse-nested-composite grids for estuaries? We are planning to keep using NPZD model.

I'm still thinking in the realistic configuration for estuaries: coarse-nested-composite grid. In this case, the estuary (composite grid) will receive the river input of fresh water and bio-chemical variables, and the outer ocean (coarse grid) will receive the outer currents and tidal forcing. The intermediate (nested grid) will communicate both regions. And the atmospheric forcing will work over the whole domain. We have ready a coarse simulation which is giving us very encouraging results for the full hydrodynamics & biochemistry simulation, but we definitively need further refinement for the estuary. I'm hoping that if I can set this configuration for the dogbone case (as soon as contact.m is fixed), then I would be able to do it for our coastal-estuary region which is already working fine with a coarse grid.

Thank you very much for your interest and help,

Alexis Espinosa
UWA

[arango]

Although I have not tried passive tracers yet in nesting, the algorithm is designed to work on all tracers (active and passive). We operate on the full 5D tracer array. If works on temperature and salinity, it will work also on passive tracers. As designed, you need passive tracers (like biology) in all grids We need to fill the contact point with the coarser grid solution. We apply also rivers to all grids but notice that each river has it own local (i,j) location. On principle, you can put the rivers in the finer grid but you need two-way to feedback information to the coarser grid.

I will recommend you to start simple and then you add grid complexity to your application. That's what I do in my application even if I am a developer.

The more that I think about the nesting algorithm, I think that the way to build nested grids is from fine to coarse grids. That is, we generate a intermediate grid at the finest resolution for the desired application and from it we extract coarse grid, refined grids, and perhaps composite grids. The will give you the data needed for all contact points. Also, we will work the land/sea masking in this intermediate grid (it will be tedious). This is the best representation of the coastline that we have at finest resolution.

It will take several trial and error to extract the optimal grids. Recall that we offer coarse2fine.m and fine2coarse.m Matlab scripts. The extraction indices Imin, Imax, Jmin, and Jmax are crucial when extracting the nested grids from the intermediate fine resolution grid.

[AlexisEspinosa]

OK. I get it.

Then, I will start working to build a full-domain refined grid. And start from there for building the composite-nested-coarse grid sequence (from fine to coarse). This will take some time and, hopefully, in the meantime, the solution for the bug in "contact.m" interpolation-weights can be solved.

It's good to know that, in principle, the biogeochemical algorithms should not suffer with the 'estuary type' grids.

And yes, I will take care and keep moving forward in simple steps. Ha ha ha, that is why I started by trying to replicate the dogbone case. I still want to create the new 'nested-composite dogbone' case before moving forward into my realistic type application. Hopefully most of the possible problems will appear during the testing of the new dogbone case, before going into the realistic one. In this way it will be easier to find help on the forum. If everything goes fine with the new dogbone, I will then move forward and modify my realistic case, also in simple steps.

Thank you very very VERY much for your help and interest in this issue.

Alexis Espinosa
UWA

[AlexisEspinosa]

OK. I found a possible solution for the out-of-range values when computing the interpolation weights, but this still needs to be revised by ROMS developers.

As stated above, indexes are becoming 0, or exceeding the limits of arrays. The modifications I suggest are described below:

0) qr is an already defined function in matlab (Orthogonal-triangular decomposition). So, I suggest using another name for the variable. I'm using qR instead.

1) RindexB (..etc) definition needs some restriction in the indexes values for not becoming zero. So, as for Rindex0, RindexT and RindexR (where a min() 'safety' function is used), a max() safety function needs to be used in this case:

Code: Select all

    RindexB = sub2ind([Ir, Jr],                                         ...
                      S.contact(cr).point.Idg_rho+1,                    ...
                      max(1,S.contact(cr).point.Jdg_rho));  % (Idg, Jdg-1)

and

Code: Select all

    UindexB = sub2ind([Iu, Ju],                                         ...
                      S.contact(cr).point.Idg_u  ,                      ...
                      max(1, S.contact(cr).point.Jdg_u));   % (Idg, Jdg-1) 

and

Code: Select all

    VindexB = sub2ind([Iv, Jv],                                         ...
                      min(Iv, S.contact(cr).point.Idg_v+1),             ...
                      max(1, S.contact(cr).point.Jdg_v-1)); % (Idg, Jdg-1)

2) index3 in linear_weights can exceed the size of the mask array. So, again, I added a safety min() function:

Code: Select all

%AEG: This index was going out of bounds: index3  = index4 + 1;                % (Idg+1, Jdg+1)
index3  = min(length(mask),index4 + 1);                % (Idg+1, Jdg+1)

3) index1,3,4,6,7,9 in quadratic_weights may go out of range. So I'm using safety min() and max() functions:

Code: Select all

index1  = max(1,index2 - 1);               % (Idg-1, Jdg-1)
index3  = min(length(mask),index2 + 1);    % (Idg+1, Jdg-1)
index4  = max(1,index5 - 1);               % (Idg-1, Jdg  )
index6  = min(length(mask),index5 + 1);    % (Idg+1, Jdg  )
index7  = max(1,index8 - 1);               % (Idg-1, Jdg+1)
index9  = min(length(mask),index8 + 1);    % (Idg+1, Jdg+1)

4) qr (or qR) usage in linear_weights and quadratic_weights is incorrect, the variable inside the function is q so the code needs the following modification:

Code: Select all

if (any(isnan(q))),
  %AEG: the use of qR here seems wrong: q(isnan(qR)) = 0;
  q(isnan(q)) = 0;
end

With these modifications I was able to generate the composite contact points (ngc) file. In principle, everything seems fine, but the results are not quite similar to the 'whole' grid case (as I will explain in a following post). Anyway, as far as I have checked, the recently generated ncg file looks pretty similar to the original test case file (dogbone_ngc_composite.nc). The only differences I can appreciate are new variables inside the file, which seem to be neccesary for the new version of ROMS. By the way, the original set of files from the example test case do not work anymore because they come from a previous ROMS version and the new neccessary variables are not included (as I will explain in the following post).

But what concerns to this post, I hope that ROMS developers can revise and test the suggested modifications and include them into a new version of the contact.m matlab tool.

[AlexisEspinosa]

I'll show a comparison below when the figures are ready.

[arango]

I already fixed this problem in contact.m. I pretty much rewrote this logic so it is more robust. However, I am not ready to release it yet. There are other changes that we need to think carefully about and discuss between us before we can release the updated scripts. I don't have the time to check your fix. You need to be patience with this. I will probably will release all these Matlab scripts by the end on the month. I will be traveling and I have other commitments that have much higher priority.

[AlexisEspinosa]

Ok. Thank you very much Prof. Arango.

Here I'm just showing the results of the 'whole' grid dogbone, 'refined' grid dogbone and 'composite' grid dogbone obtained with my modifications to contact.m and ROMS 758. Please consider the problems shown here while tackling the issues mentioned in this thread. Both the composite and the refined case show some problems, although the composite case problems are more severe.

As I said in the post above, the refined and composite test cases do not run because the files were prepared for a previous version of ROMS and some variables are missing from the files. These are the messages for the composite attempt to run:

Code: Select all

Output/Input Files:

             Output Restart File:  dogbone_rst_right.nc
             Output History File:  dogbone_his_right.nc
            Output Averages File:  dogbone_avg_right.nc
                 Input Grid File:  ./Data/dogbone_grd_right.nc
  Nesting grid connectivity File:  ./Data/dogbone_ngc_composite.nc
    Input Nonlinear Initial File:  ./Data/dogbone_ini_right.nc

 NETCDF_GET_DIM - requested dimension: nLweights
                  not found in input file:  ./Data/dogbone_ngc_composite.nc
                  call from:

 Elapsed CPU time (seconds):


 ROMS/TOMS - Output NetCDF summary for Grid 01:

 ROMS/TOMS - Output NetCDF summary for Grid 02:

 ROMS/TOMS - Input error ............. exit_flag:   2


 ERROR: Abnormal termination: NetCDF INPUT.
 REASON: No error

So, clearly an update of the cases is needed. For that, first the contact.m function needs to be fixed. (Thank you very much for your help and interest.)




The composite and refined grids (copied from 2013 Arango's post viewtopic.php?f=13&t=2912) are:












When I use my suggested modifications, the following results are obtained.

Free surface at time: 9 hrs from initial time. (4th record in the history file)


Fig1. Free surface, WHOLE




Fig2. Free surface, COMPOSITE




Fig3. Free surface, REFINED

The three figures look very similar, so it seems that the hydrodynamics is working fine.







But for temperature and salinity results look weird and wrong.

Salinity after 57 hours from the initial condition:


Fig4. Salinity, WHOLE.




Fig5. Salinity, COMPOSITE




Fig6. Salinity, REFINED

Salinity is being affected on the east-west closed boundaries in the composite grid case. This is weird and unexpected, as I would expect problems close to the contact plane (mid vertical axis of the bone). Haha, but the problems appear on the non-contact boundaries.

Salinity is also being affected in the refined case. Although the difference with respect to the whole case is minimal, two clear non-physical spots of lower salinity appear close to the refined grid boundary (which is the neck of the bone).






Now temperature.
Temperature after 57 hours from the initial condition:


Fig7. Temperature, WHOLE




Fig8. Temperature, COMPOSITE




Fig9. Temperature, REFINED

(The initial condition for temperature has a region of low temperature close to the west side of the bone (clear blue zone)).

For the refined case, two non-phyical spots of higher temperature appear close to the refined grid (neck of the bone). This may seem minimal, but they affect the solution and the field is not symmetrical anywhere in the domain (not shown).

For the composite case, temperature is being strongly affected in the east-west boundaries too. Again, this looks very weird.




Thank you very much. I'm very thankful for your efforts in solving these issues. I'll be patient and wait for the new corrections & modifications to come. Thank you very much.

Alexis Espinosa
UWA

[TorresGarcia]

I already fixed this problem in contact.m. I pretty much rewrote this logic so it is more robust. However, I am not ready to release it yet. There are other changes that we need to think carefully about and discuss between us before we can release the updated scripts. I don't have the time to check your fix. You need to be patience with this. I will probably will release all these Matlab scripts by the end on the month. I will be traveling and I have other commitments that have much higher priority.

Dear Dr. Arando,

I have been working in a four nested refine grid and I am having problem running my model once I incorporate the four grid. When I create the contact file for three grids I had to change the pn and pm values I was obtaining and it worked fine. However, now that I ran the contact.m for the four grid is not working.

PUT_REFINE2D - unbounded contact points temporal: interpolation:
cr = 03 dg = 02 ng = 03
iic(dg) = ******* told = 1 tnew = 2
iic(ng) = 0000008 Wold = 1.83333 Wnew = -0.83333
time(ng) = ******* time(told) = ******* time(tnew) = *******

PUT_REFINE2D - unbounded contact points temporal: interpolation:
cr = 03 dg = 02 ng = 03
iic(dg) = ******* told = 1 tnew = 2
iic(ng) = 0000008 Wold = 1.83333 Wnew = -0.83333
time(ng) = ******* time(told) = ******* time(tnew) = *******

PUT_REFINE3D - unbounded contact points temporal: interpolation:
cr = 03 dg = 02 ng = 03
iic(dg) = ******* told = 1 tnew = 2
iic(ng) = 0000008 Wold = 1.83333 Wnew = -0.83333
time(ng) = ******* time(told) = ******* time(tnew) = *******

I ran it using Lmask true and it ran for few times step, howver is not taking the right contact points.

Searching on this page I found this message that you wrote about fixing the contact.m function. Is the later version the one from January 2015 (that is the one I checked at svn), if so does that takes four grids?

Thank you for your time,
Legna

[arango]

That's your problem! You need to use Lmask=false when generating the contact points with the Matlab script. Please read the instructions! Notice that the default value is false.

I have said several times in this forum that if we don't have experience in nesting applications for a particular region, we should always always always start simple. Do first two grids, then three, and so on. This recommendation applies to everybody: us developers and experienced and novice users !!!!

I always ignore posts of users going immediately to very complex nesting configurations. This stuff is difficult and requires a lot of try and error and the success of a nesting application dependents on the grids topology, bathymetry, and dynamical circulation. There is a lot to consider.

[TorresGarcia]

Thank you for replying to my question. This is an application that we have been working on for a number of months, starting from a single grid and building nested applications one at a time. The three nested grid system worked but this error occurred when adding the fourth grid and colleagues I asked from Woods Hole have not seen this error before. Strangely, the error for the four grid nest references grids 2 and 3 but the three nest system runs with no problems.

I should have mentioned that in all of our initial testing we use the default Lmask (false). Our fourth grid (centered on a coral reef) does not include land so changing the Lmask value was a list ditch attempt and allowed the model to run for a few timesteps, but clearly it is not the solution to the problem.

It seems we can get the model to run by making the Dt for grids 1 and 2 the same, otherwise we get the error previously posted.

Our primary concern was that we are not using your updated scripts and if that would solve the problem. The March 2015 post (above) indicated they would be released in April so we wanted to check if they had been released or if their was a new timeline. The svn repository appears to have codes older than March 2015. We're happy to be patient for the updated codes but just wanted to check on their release.

We will continue to debug but any insight you may have is appreciated!

Thanks!

Legna

[kate]

I am trying this with two files:

Code: Select all

grids = {'grid_Palau_1.nc', 'grid_Palau_2.nc'};

name = 'contacts_Palau.nc'
c_contact(name, true, 2)
contact(grids, name)

and this is what I get:

Code: Select all

 Summary of Contact Regions Processed:
  
   Grid 01: grid_Palau_1.nc
   Grid 02: grid_Palau_2.nc
  
   Contact   Donor   Receiver
    Region    Grid       Grid
  
  
Reference to non-existent field 'contact'.

Error in contact (line 316)
  dg = S.contact(cr).donor_grid;

Error in foo (line 5)
contact(grids, name)

This is coming about because grid_connections.m can't figure out the grid connectivity. Looking deeper, it seems coarse2fine adds stuff that grid_connections needs:

Code: Select all

		:parent_grid = "/Users/arango/ocean/repository/TestCases/dogbone/Data/dogbone_grd_hfull.nc" ;
		:parent_Imin = 14 ;
		:parent_Imax = 30 ;
		:parent_Jmin = 6 ;
		:parent_Jmax = 11 ;
		:refine_factor = 3 ;

I didn't use coarse2fine, but rather used Jamie's fine2coarse plus extract_grid, both with parent grid of grid_Palau_0.nc covering the whole region at high resolution. Has anyone gotten things to work this way?

[arango]

Yes, it needs that information from the global attributes. Otherwise, contact.m will not work. These scripts are complex and there is a lot of internal details. They need to be used as a full package. I won't guarantee that a mixture of third party scripts will work together. There is a lot of subtlety here.

The ones that I haven't released make an attend to improve volume conservation. I have been busy working on other codes. I need to come back to this and release the updated scripts. There some issues and opinions that I need to consolidate.

[jpringle]

Hernan & Kate --

Kate -- to start to get things to work, I have had to add the refine_factor attributes by hand. My efforts have been interrupted by my vacation, but will pick up in two weeks or so. I am not sure if I will end up going fine2coarse or coarse2fine, then edit the bathymetry in the resulting fine-grid. I am still struggling with how to get things to work when the refined grid rests on a closed boundary.

Hernan -- there are those of us who are working carefully from simple examples, and reading your comments closely. In them, in several places, you advocate starting with a fine grid and using fine2coarse. However, the script in the matlab repository has syntax and logical errors. I attempted to fix those and provide test cases in my script, which follows yours as closely as possible, given my limited experience. We are trying in good faith to do things properly.

Hernan -- can you help us with how you want refine_factor to be used? There are hints in the code of a coarse_factor, but it is not used in contacts.m. Right now I am just adding refine_factor as appropriate by hand to the grids, but it would be helpful, perhaps, to include an option to set refine_factor to grid_extract.m. I can do so, but unsure of the exact roll of these attributes (in particular, coarse_factor).

Thanks, Jamie

[kate]

As far as I can tell, these two are equivalent:

Code: Select all

grid_extract('grid_Palau_0.nc', 'grid_Palau_3.nc', 781,1381,652,1132)

Code: Select all

ncks -d xi_rho,781,1382 -d xi_psi,781,1381 -d xi_u,781,1381 -d
xi_v,781,1382 -d eta_rho,652,1133 -d eta_psi,652,1132 -d eta_u,652,1133 -d
eta_v,652,1132 grid_Palau_0.nc grid_Palau_2.nc

Is that true?

I did get the contact.m thing to work after:

Code: Select all

ncatted -O -a refine_factor,global,a,l,3 grid_Palau_2.nc
ncatted -O -a parent_Imin,global,a,l,261 grid_Palau_2.nc
ncatted -O -a parent_Imax,global,a,l,461 grid_Palau_2.nc
ncatted -O -a parent_Jmin,global,a,l,218 grid_Palau_2.nc
ncatted -O -a parent_Jmax,global,a,l,378 grid_Palau_2.nc
ncatted -O -a parent_grid,global,a,c,'grid_Palau_1.nc' grid_Palau_2.nc

I hope to be testing it later today, but there's more setup to do.

If one were to build tools from scratch, I've been hearing good things about Julia. Notably, it is the first language since C/Fortran with raw scientific computing speed as a stated goal. The lead developer has a Matlab background so it has some Matlabish choices, but the best thing is that it's open source. Just saying.

[TorresGarcia]

Thanks all for the discussion on this. Just for information, we used coarse2fine to create all of our nested grids. We are using the COAWST malab tools (current repository version) and are able to run the contact.m file to generate the nests. All the contact points appear correct and we have compared water level and velocity output in at overlapping points when running the two and three nested grid systems.
Legna