Matlab nesting pre-processing program error.

[d.kobashi]

Hi,

I am trying to create nested grid and contact netcdf file using ROMS Matlab routine.

I made a child grid, but when I tried to create contact file, I had a couple of errors.
By the way, I've created a grid file with spherical coordinate.

First error is near line 720 in contact.m

Code: Select all

Undefined function or variable 'FCr'.
Error in contact>refine_coordinates (line 724)
    FSr.Values = G(dg).angle(:);    R.angle = FCr(XrF, YrF);

The error is related to FCr array.

Originally,

Code: Select all

  elseif (~isempty(G(dg).lon_rho) && ~isempty(G(dg).lat_rho))

    FSr.Values = G(dg).angle(:);    R.angle = FCr(XrF, YrF); 
    FSr.Values = G(dg).f(:);        R.f     = FCr(XrF, YrF); 
    FSr.Values = G(dg).h(:);        R.h     = FCr(XrF, YrF); 
    

The above lines should be

Code: Select all

  elseif (~isempty(G(dg).lon_rho) && ~isempty(G(dg).lat_rho))

    FSr.Values = G(dg).angle(:);    R.angle = FSr(XrF, YrF); 
    FSr.Values = G(dg).f(:);        R.f     = FSr(XrF, YrF); 
    FSr.Values = G(dg).h(:);        R.h     = FSr(XrF, YrF); 

as I guess FCr is for cartesian coordinate whereas, FSr is for spherical coordinate. Please correct me if I am wrong.

Then, after I modified the above lines and ran the script, I had another error for contact.m.

Code: Select all

Subscripted assignment between dissimilar structures.
Error in contact (line 329)
    S.refined(cr) = R;

The cell array, R is as follow.

Code: Select all

    spherical: 1
      uniform: 0
       xi_rho: [1117×877 double]
      eta_rho: [1117×877 double]
      lon_rho: [1117×877 double]
      lat_rho: [1117×877 double]
      Irg_rho: [1117×877 double]
      Jrg_rho: [1117×877 double]
       xi_psi: [1116×876 double]
      eta_psi: [1116×876 double]
      lon_psi: [1116×876 double]
      lat_psi: [1116×876 double]
      Irg_psi: [1116×876 double]
      Jrg_psi: [1116×876 double]
         xi_u: [1116×877 double]
        eta_u: [1116×877 double]
        lon_u: [1116×877 double]
        lat_u: [1116×877 double]
        Irg_u: [1116×877 double]
        Jrg_u: [1116×877 double]
         xi_v: [1117×876 double]
        eta_v: [1117×876 double]
        lon_v: [1117×876 double]
        lat_v: [1117×876 double]
        Irg_v: [1117×876 double]
        Jrg_v: [1117×876 double]
        angle: [1117×877 double]
            f: [1117×877 double]
            h: [1117×877 double]
           pm: [1117×877 double]
           pn: [1117×877 double]
         dndx: [1117×877 double]
         dmde: [1117×877 double]
     mask_rho: [1117×877 double]
     mask_psi: [1116×876 double]
       mask_u: [1116×877 double]
       mask_v: [1117×876 double]
    spherical: []
      uniform: []
       xi_rho: []
      eta_rho: []
        x_rho: []
        y_rho: []
      lon_rho: []
      lat_rho: []
      Irg_rho: []
      Jrg_rho: []
       xi_psi: []
      eta_psi: []
        x_psi: []
        y_psi: []
      lon_psi: []
      lat_psi: []
      Irg_psi: []
      Jrg_psi: []
         xi_u: []
        eta_u: []
          x_u: []
          y_u: []
        lon_u: []
        lat_u: []
        Irg_u: []
        Jrg_u: []
         xi_v: []
        eta_v: []
          x_v: []
          y_v: []
        lon_v: []
        lat_v: []
        Irg_v: []
        Jrg_v: []
            h: []
            f: []
        angle: []
           pm: []
           pn: []
         dndx: []
         dmde: []
     mask_rho: []
     mask_psi: []
       mask_u: []
       mask_v: []

Not sure why there are empty cell arrays for cr=2 (Ncontact=2). That may be the problem. Does anyone have any idea how to resolve the problem?

I downloaded the Matlab routine today from
https://www.myroms.org/svn/src/matlab/grid/

Thanks in advance,

-DJ

[arango]

Yes, your first correction is fine. Thank you I updated the repository. We didn't catch that one before because all our grids always have the (XI,ETA) coordinates even if the applications grid is spherical. It turns out that in curvilinear grids is better to do the interpolation in (XI,ETA) space instead of (lon,lat) to avoid interpolant errors due to the curvature of the grid cell. The (XI,ETA) space is rectangular and Cartesian. We need to think hard about this issue, but it is the best way to proceed. That's the reason that conditional if-statement is ordered in the following way:

Code: Select all

 if (~isempty(G(dg).x_rho) && ~isempty(G(dg).y_rho))
    if (UseGriddedInterpolant)
      ...
    else
      ...
    end
  end

  if (spherical && ~isempty(G(dg).lon_rho) && ~isempty(G(dg).lat_rho))
    if (UseGriddedInterpolant)
      ...
    else
      ...
    end
  end 

So the Cartesian interpolant function is always available and it advantageous to use it for particular grid variables. The error that you get implies that your original grid does not have x_rho, y_rho, x_u, y_u, x_v, y_v. It is better to compute the inverse grid spacing variables pm and pn from x_rho and y_rho.

Now, your second error implies that something is wrong with your coarse grid structure that is used to generate the refined grid. The S and R structures are dissimilar and do not include the same fields. You need to check S and R in the debugger to compare the field components. Notice that S is an array of structures.

[d.kobashi]

Hernan,

Thanks for your input.

Now, your second error implies that something is wrong with your coarse grid structure that is used to generate the refined grid. The S and R structures are dissimilar and do not include the same fields. You need to check S and R in the debugger to compare the field components. Notice that S is an array of structures.

The error is gone after I added x_/y_ (rho, u, v, psi) to the parent grid file.
To be honest, this is a strange error and should not happen because if we run ROMS with spherical coordinate, including x_/y_ (rho, u, v, psi) is an option, not a requirement. Should that be a bug? or if we do nesting, adding x_/y_ is required? Anyway, I will start testing the nesting model with my refined grids and contact file soon.

For non-related matter, you mentioned that it is better to compute pm, pn from x_rho, y_rho. However, I usually run ROMS with spherical coordinate and I need to convert lon/lat to x/y after I create ROMS grids. And that depends on map projection. What map projection gives the most accurate on distance? I thought about UTM, but usually my grids cover multiple UTM zones so I cannot use UTM. Do you have any script to convert lon/lat to x/y? I use Basemap or Cartopy to convert with the Mercator projection, but I guess Mercator is not the best option. What would ROMS users usually do on this matter?

Also, when creating refined grids, I noticed that the value for "spherical" in a grid file is 0 or 1 (integer), not 'T' or ' ' (char). Have you changed it recently or was it that way from the beginning? I use spherical = 'T' in all of my grid files and have not had any issues.

Thanks in advance.

-DJ