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.
Circumpolar grid
Re: Circumpolar grid
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.
-
- Posts: 2
- Joined: Fri Nov 13, 2009 6:21 pm
- Location: Instituto Oceanográfico - University of Sao Paulo
Re: Circumpolar grid
Many thanks Kate to reply!!! I finally could reach the grid I wanted (I think....)
Cheers,
Virna
Cheers,
Virna
Re: Circumpolar grid
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
May be this script? Sometimes it works normally
![Smile :)](./images/smilies/icon_smile.gif)
(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);
- Attachments
-
glob_grid.m
- (7.21 KiB) Downloaded 295 times
Re: Circumpolar grid
I sincerely appreciate your file.
It will be very helpful.
Thank you so much.
Best,
-JH
It will be very helpful.
Thank you so much.
Best,
-JH
Joonho Lee