Circumpolar grid

[vmeccia]

Hi!
I'm a new Roms user. I'm planning to implement the model in the Southern Ocean. So I would like to construct a circumpolar grid with variable horizontal resolution. I'm dealing with Seagrid but I'm not sure that it is the best tool for constructing such a grid since it doesn't use an azimuthal projection (or I didn't find it?).
Does anybody know about an adequate tool for doing that? Does anybody have a suggestion?
Please accept my apologies if it is not the correct place to post my doubt, but your opinion will be very valuable for me.
Cheers.

[kate]

How do you want your resolution to vary? What I would do is create an analytic polar coordinate grid then use editmask to mask out the land. There is a LAB_CANYON example in ROMS to do something similar, though you would have to make an external grid file to use editmask.

[vmeccia]

Many thanks Kate to reply!!! I finally could reach the grid I wanted (I think....)
Cheers,
Virna

[susonic]

Hi,

Is there a toolbox to generate attached figure like domain?

[Pysh]

Hi

May be this script? Sometimes it works normally

(Gurvan Madec, Maurice Imbard; A global ocean mesh to overcome the North Pole singularity; Climate Dynamics (1996) 12:381-388)

Cheers
Boris

Code: Select all

function glob_grid

global Jeq JM dd

% ---- net parameters ------------

IM = 280;   % number of longitudes (without additional overlappimg point for periodic variant)
JM = 275;   % number of latitudes from pole to pole
Jeq = 115;  % South-North latitude number 
Lan = 40;   % latitude of north pole
Lon = 80;   % longitude of north pole
Lae = 0;    % South-North bourder latitude
a2b = 2.5;  % maximum axis relation in the north pole
a_fr_fg = 0;  % 1 - minimum axis calculate from f and g, else maximum axis
Jel = 200;  % Start number of ellipse

%  ----- parameters for correction longitudes on 180 meridian (smoothing)
dnl = 6;  % --- numbers of longidudes from 180 meridian --------
dns = 3;   % ---- width od smoothing window ------------ 

%  ----- parameters of grid for latitudes limits -------------
nlan = 25;   % number of the last latitude in the grid from north pole
nlas = 9;    % number of first latitude in the grid from south pole
batmax = 0;  % maximum of depth in bathymetry data

% --------- steps for grid picture ---------------
lastp = 1;
lostp = 1;

%  ------- constraction of f and g functions for cicles on stereographic projection --------------------
% ((x/a)^2 + ((y-b-f)/b)^2 = 1)

for j = 1:Jeq
 af(j) = -pi + (pi / 2 + Lae * pi / 180) / (Jeq - 1) * (j - 1);   
 ag(j) = pi + (-Lae * pi / 180 - pi / 2 )/ (Jeq - 1) * (j - 1);  
end

alf_to = pi / 2 - Lan * pi / 180;  
alff_fr = -pi / 2 + Lae * pi / 180;  
alfg_fr = pi / 2 - Lae * pi / 180;  

delaf = alf_to - alff_fr;
df = (alff_fr - (-pi)) / (Jeq - 1);
dd = delaf / df;
cef = fzero(@cc, 1 / dd + 0.001);
bef = delaf / (exp(cef * JM) - exp(cef * Jeq));
aef = alff_fr - bef * exp(cef * Jeq);

delag = alf_to - alfg_fr;
dg = (alfg_fr - pi) / (Jeq - 1);
dd = delag / dg;
ceg = fzero(@cc, -1 / dd + 0.001);
beg = delag / (exp(ceg * JM) - exp(ceg * Jeq));
aeg = alfg_fr - beg * exp(ceg * Jeq);

for j = Jeq+1:JM
 af(j) = aef + bef * exp(cef * j);   
 ag(j) = aeg + beg * exp(ceg * j);   
end

k(1:Jel) = 1;
for j = Jel+1:JM
 k(j) = 1 + ((a2b - 1) / (JM - Jel) * (j - Jel))^2; 
end

f = sin(af);
g = sin(ag);

if a_fr_fg  == 0
 b = (g - f) / 2;
 a = b .* k;
else
 a = (g - f) / 2;
 b = a ./ k;
end    

% ---- points on equatorial cicle --------------
dx = 2 * pi / IM;
for i = 1:IM
 x(i,Jeq) = cos(-pi / 2 + (i - 1) * dx);
 y(i,Jeq) = sin(-pi / 2 + (i - 1) * dx);
end

% ------ points on initial longitude line ----------
for j = Jeq+1:JM-1
 x(1,j) = 0;
 y(1,j) = f(j);
 x(IM/2+1,j) = 0;
 y(IM/2+1,j) = g(j);
end

% ---- points coordinates constraction from condition of  
% ( {x4-x1,y4-y1} . {x-x4,y-y4}) = 0   (scalar production of the
% perpendicular vectors) and equation for ellipse{x,y} x^2/a^2 + (y^2 - (f +
% b))^2/b^2 = 1

for j = Jeq+1:JM-1
 for i = 1:IM/2
  x1 = x(i,j-1);
  x4 = x(i+1,j-1);
  y1 = y(i,j-1);
  y4 = y(i+1,j-1);
  dx = x1 - x4;
  dy = y1 - y4;
  a1 = y4 + x4 * dx / dy;
  b1 = -dx / dy;
  c1 = b(j)^2 + a(j)^2 * b1^2;
  c2 = 2 * a(j)^2 * b1 * (a1 - f(j) - b(j));
  c3 = a(j)^2 * (a1 - f(j) - b(j))^2 - (a(j)* b(j))^2;
  det = max(0, c2 * c2 - 4 * c1 * c3); 
  x(i+1,j) = (-c2 + det^0.5) / 2 / c1;
  y(i+1,j) = a1 + b1 * x(i+1,j);
 end    
end

% --------- simmetrics for longitudes -----------
for j = Jeq+1:JM-1
 for i = IM/2+2:IM
  ifr = IM - i + 2;  
  x(i,j) = -x(ifr,j); 
  y(i,j) = y(ifr,j); 
 end    
end    

% ----------- south hemisphere (not disturbed) -----------------
dla2e = (Lae + 90) / (Jeq - 1);
for j = 1:Jeq
 for i = 1:IM   
  lon(i,j) = (i - 1) * 360 / IM; 
  lat(i,j) = -90 + (j - 1) * dla2e;
 end
end

for j = Jeq+1:JM-1
 for i = 1:IM
  if x(i,j) == 0 && y(i,j) < 0    
   lon(i,j) = 0; 
  elseif x(i,j) == 0 && y(i,j) > 0    
   lon(i,j) = 180; 
  elseif x(i,j) > 0     
   lon(i,j) = 90 + 180 / pi .* atan(y(i,j)/x(i,j)); 
  elseif x(i,j) < 0    
   lon(i,j) = 270 + 180 / pi .* atan(y(i,j)/x(i,j)); 
  end
 end
end

% ---------- projection from stereograpthic plate onto sphere -----------
lat(:,Jeq+1:JM-1) = 90 - 180 / pi .* asin((x(:,Jeq+1:JM-1).^2 +y(:,Jeq+1:JM-1).^2).^0.5); 

for j = Jeq+1:JM-1
 if lon(1,j) == 0
  lon(2,j) = lon(3,j) / 2;
  lon(IM,j) = (360 + lon(IM-1,j)) / 2;
 else
  lon(2,j) = (180 + lon(3,j)) / 2;
  lon(IM,j) = (180 + lon(IM-1,j)) / 2;
 end
 lat(2,j) = correct_lat(lon(2,j), a(j), b(j), f(j), y(2,j));
 lat(IM,j) = correct_lat(lon(IM,j), a(j), b(j), f(j), y(IM,j));
end


% ----- additional point for periodic variant --------------
lon(IM+1:IM+2,:) = lon(1:2,:);
lat(IM+1:IM+2,:) = lat(1:2,:);

% ------- smoothing longitudes near 180 meridian to correct algorithm
% accumulated error
lon1 = lon;
lons = lon;
sdns = 2 * dns + 1;
for j = Jeq:JM-1
 for i = IM/2-dnl:IM/2+dnl
  lon1(i,j) = sum(lons(i-dns:i+dns,j)) / sdns;
  lat1(i,j) = correct_lat(lon1(i,j), a(j), b(j), f(j), y(i,j));
 end
end 

% -------- shift longitudes in order to initial longidute will be 0 ----
lon1 = mod(lon1+Lon+180,360);
lat1 = lat;

%nl = floor(Lon / 360 * IM);
nl = ceil(Lon / 360 * IM);
dn = IM / 2 - nl;
for i = 1:IM-dn
 lon2(i,:) = mod(lon1(i+dn,:),360);   
 lat2(i,:) = lat1(i+dn,:);    
end
for i = IM-dn+1:IM+1
 lon2(i,:) = mod(lon1(i-IM+dn,:),360);    
 lat2(i,:) = lat1(i-IM+dn,:);    
end    

lon3 = lon2(:,nlas:JM-nlan);
lat3 = lat2(:,nlas:JM-nlan);

samplefactor = 1;
[bat, refvec] = etopo('ETOPO5.dat', samplefactor);

dloe = 360 / 4320 / samplefactor;
dlae = 180 / 2160 / samplefactor;

[m, n] = size(lat3);

for i = 1:m
 for j = 1:n
  ii = max(floor(lon3(i,j) / dloe),1);
  jj = max(floor((lat3(i,j) + 90) / dlae),1);
  bat3(i,j) = bat(jj,ii);
 end
end 

for i = 1:m
 for j = 1:n
  if isnan(bat3(i,j))
   i1 = i + 1;
   if i1 == m + 1; i1 = 1; end 
   i2 = i - 1;
   if i2 == 0; i2 = m; end;
   j1 = min(n,j+1);
   j2 = max(1,j-1);
   si = 0; sb = 0;
   if ~isnan(bat3(i1,j1)); si = si + 1; sb = sb + bat3(i1,j1); end 
   if ~isnan(bat3(i2,j2)); si = si + 1; sb = sb + bat3(i2,j2); end 
   if ~isnan(bat3(i1,j2)); si = si + 1; sb = sb + bat3(i1,j2); end 
   if ~isnan(bat3(i2,j1)); si = si + 1; sb = sb + bat3(i2,j1); end 
   if si > 0
    bat3(i,j) = sb / si;
   else
    bat3(i,j) = NaN;   
   end    
   astop = 1;   
  end    
 end
end 

bat3(bat3 > batmax) = 10000; 

% ------- saving grid to file -------------

rho.lon = lon3;
rho.lat = lat3;
rho.depth = bat3;
save toroms_glob.mat rho

%  ------- pictures of grid with bathymetry -------------
figure
axesm ('MaPprojection','glob','Meridianlabel','on','Parallellabel','on');
hold on
surfm(lat3,lon3,bat3);

for j = 1:lastp:n
 plotm(lat3(:,j), lon3(:,j),'w');
 astop = 1;
end
for i = 1:lostp:m
 plotm(lat3(i,:), lon3(i,:),'w');
 astop = 1;
end
astop = 1;

function lat = correct_lat(lon, a, b, f, yo)
 b1 = -1 / tan(lon * pi / 180);
 c1 = b^2 + a^2 * b1^2;
 c2 = -2 * a^2 * b1 * (f + b);
 c3 = a^2 * (f^2 + 2 * f * b);
 det = max(0, c2 * c2 - 4 * c1 * c3); 
 dy1 = yo - b1 * (-c2 + det^0.5) / 2 / c1;
 dy2 = yo - b1 * (-c2 - det^0.5) / 2 / c1;
 if abs(dy1) > abs(dy2)
  x = (-c2 - det^0.5) / 2 / c1;
 else
  x = (-c2 + det^0.5) / 2 / c1;
 end    
 y = b1 * x;
 lat = 90 - 180 / pi * asin((x^2 +y^2)^0.5); 

function cf = cc(x)
global Jeq JM dd 
cf = x - log(1 + x * dd) / (JM - Jeq); 

[susonic]

I sincerely appreciate your file.
It will be very helpful.
Thank you so much.

Best,
-JH