ROMS grid for conservative remapping

[c.drinkorn]

Hi all!

I am currently trying to use conservative remapping of ROMS data instead of my previous bilinear method. For conservation the cell areas are needed and therefore the cell corners. I thought those were lon_psi and lat_psi, respectively. But the remapping routine complains about missing corner points in my source grid. When I make a grid description file I realize that lon_psi and lat_psi are all zeros. How is that?

Thank you!

[kate]

How did you make your grid? Some grid generators add more points to the outside of the model domain than others. You can just extrapolate what you need from the interior. Or if all the lat_rho, lon_rho values are zero, then there's a problem in the generation phase.

[c.drinkorn]

I didn't make the grid myself, it's the one from the former metno Arctic20km model. The original grid file has no psi coordinates. How can I extrapolate them from the rho coordinates?
Also: I have experienced difficulties making a grid of my own with the North Pole in it (that's why I re-use Arctic20km). Any suggestions on how to achieve that?
Thank you!

[kate]

Rob Cermak has been working on https://github.com/ESMG/gridtools for making MOM6 grids, including handling polar grids. There's a ROMS-to-MOM6 grid tool and one could easily make the reverse.

I make my Arctic grids using old, old grid tools which still work and still produce the psi fields along with the rest.

[c.drinkorn]

Thanks Kate!
I am not much of a python expert (on my list!) but I worked my way through the installation tutorials for a local version of the grid tool. However, when I try to load A20 into the tool I get an "unloadable" error. Maybe I should also look into the "old old" grid generation tools instead since I suppose A20 was made with them, too. Do you have a link for that? I think, I came across a manual once about a grid generator of which you were one of the authors but I can't find it anymore...

[c.drinkorn]

I think I found it! Is it gridpak?

[wilkin]

I struck the same problem working from the Met Norway grids and wrote simple code to add the lon_psi etc.

Code: Select all

kiwi:sandbox wilkin$ more roms_add_psi_latlon.m
function g = roms_add_psi_latlon(ncfile,~)
% grd = roms_add_psi_latlon(grd_file,[writeit])
%    or roms_add_psi_latlon(grd_structure,[writeit])
%
% If the grid structure returned by roms_get_grid does not have 
%    lon_psi and lat_psi, this function will add them to the structure
%
% If a ROMS grid file does not have lon_psi and lat_psi, compute them and
%    add them to the output grid structure, and optionally add them to the
%    original grid file by giving a second input argument.
%
% If there is a rho mask, compute the psi mask for it.
%
% John Wilkin - 24 May 2018 (for MARTINI)

if nargin > 1
  writeit = 1;
else
  writeit = 0;
end

if ~isstruct(ncfile)
  g = roms_get_grid(ncfile);
else
  g = ncfile;
  writeit = 0;
end

addmask = 0;
if isfield(g,'mask_rho')
  addmask = 1;
end

% The psi points grid is between opposite u or v points, or at the center
% of 4 rho points. For no particular reason I choose to avewrage the
% position estimated separately from the u and v grids. 
lon_psi1 = av2(g.lon_u);
lon_psi2 = av2(g.lon_v')';
g.lon_psi = 0.5*(lon_psi1+lon_psi2);
lat_psi1 = av2(g.lat_u);
lat_psi2 = av2(g.lat_v')';
g.lat_psi = 0.5*(lat_psi1+lat_psi2);

if addmask
  [~,~,tmp] = uvp_masks(g.mask_rho');
  g.mask_psi = tmp';
end

if writeit
  I = nc_inq(ncfile);
  S.Dimensions = I.Dimensions;
  S.Variables(1) = roms_metadata('lon_psi', 1);
  S.Variables(2) = roms_metadata('lat_psi', 1);
  if addmask
    S.Variables(3) = roms_metadata('mask_psi', 1);
  end
  nc_append(ncfile, S);
  nc_write(ncfile,'lon_psi',g.lon_psi');
  nc_write(ncfile,'lat_psi',g.lat_psi');
  if addmask
    nc_write(ncfile,'mask_psi',g.mask_psi');
  end
else
  % Don't add the new lon/lat to the file, just return the modified grid
  % structure
end

[kate]

Yes, gridpak on github. What do you want in an Arctic grid? I've got two you could try as well.

[c.drinkorn]

@John: Thanks! I am trying to make use of this now but can't seem to get around what av2 is supposed to be. Am I lacking a particular set of functions here?

@Kate: We set up an Arctic/N.Atlantic sed-ocean-ice model for several sensitivity tests. E:g. our first paper: https://doi.org/10.1016/j.palaeo.2021.110391
You need the entire Arctic plus ice etc in order to reconstruct the least representation of the entire system. It still has its deficiencies but I guess that's modelling.

[kate]

I bet av2 is simply an average of the two arguments.

Your domain is similar but different from my Arctic domains. Mine cover all the Bering, but less of the Atlantic.

[wilkin]

I've attached here another approach to making a simple grid using the Matlab mapping toolbox to set a Lambert conformal projection, and create a rectangular grid in Lambert space (optionally rotated w.r.t. east-north) that is then projected back to geographic coordinates.

It takes a bit of trial and error to select the grid spacing in Lambert space.

This uses Hernan Arango's Matlab routines to create the gird file so it will have all the necessary variables and metadata to play nicely with the most recent version of the code.

(For completeness and other users av2.m attached here also, but I sent that via private message too).

[kate]

You'd have to try a polar projection instead of Lambert conic. Is that supported by your software?

[wilkin]

Good point, Kate.

The mapping toolbox does include the stereographic projection (there are 72 projections in total) but I have no idea how to drive it, and will resist the temptation to delve into some recreational Matlab today to try and figure it out.