tutorial_global_oce_latlon/diags_matlab/mit_oceanmasks.m

function grid = mit_oceanmasks(grid,plotit)
% function grid = mit_oceanmasks(grid)
% create masks for the 4x4 degree run for the individual oceans
% hfacc/s/w have to be available for this
  
% $Header:  $
% $Name:  $

  if grid.nx ~= 90 | grid.ny ~= 40 
    error(sprintf( ...
        '%s\n%s\n%s\n%s', ...
        'Your horizontal grid has dimensions different from 90x40.', ...
        'mit_oceanmasks creates masks for the individual ocean basins,', ...
        'therefore it needs to be customized if you change the', ...
        'horizontal dimensions.'))
  end
    
  % masks for the zonal velocities:
  % hfacw
  
  % Atlantic mask
  grid.atlantic_hfacw = grid.hfacw;
  ixsow = 6:73; % southern ocean and most of pacific
  iysow = 1:24; 
  grid.atlantic_hfacw(ixsow,iysow,:) = NaN;
  ixnpw = 15:67; % northern pacific
  iynpw = 25:37;
  grid.atlantic_hfacw(ixnpw,iynpw,:) = NaN;
  % Pacific mask
  grid.pacific_hfacw = grid.hfacw;
  ixaw = 74:90; % most of the atlantic
  iyaw = 1:40;
  grid.pacific_hfacw(ixaw,iyaw,:) = NaN;
  ixnaw = 68:73; % north atlantic and gulf of Mexico
  iynaw = 27:40;
  grid.pacific_hfacw(ixnaw,iynaw,:) = NaN;
  ixcaw = 70:73; % carribean
  iycaw = 24:26;
  grid.pacific_hfacw(ixcaw,iycaw,:) = NaN;
  ixiw = 1:27; % most of the indic
  iyiw = 1:40;
  grid.pacific_hfacw(ixiw,iyiw,:) = NaN;
  ixauw = 28:35; % water south and west of Australia
  iyauw = 1:20;
  grid.pacific_hfacw(ixauw,iyauw,:) = NaN;
  
  
  % masks for the meridional overturning stream functions:
  % hfacs
  
  % Atlantic mask
  grid.atlantic_hfacs = grid.hfacs;
  ixso = 5:73; % southern ocean
  iyso = 1:24;
  grid.atlantic_hfacs(ixso,iyso,:) = NaN;
  ixnp = 14:65; % northern pacific
  iynp = 25:36; 
  grid.atlantic_hfacs(ixnp,iynp,:) = NaN;
  ixsoc = 1:90; % complete southern ocean
  iysoc = 1:12; % south of latg(12) = -34 degN;
  grid.atlantic_hfacs(ixsoc,iysoc,:) = NaN;
  
  % Pacific mask
  grid.pacific_hfacs = grid.hfacs;
  ixatl = [1:26 74:90]; % most of altantic and indian ocean
  iyatl = 1:grid.ny;
  grid.pacific_hfacs(ixatl,iyatl,:) = NaN;
  ixatlf = 73; % small fix
  iyatlf = 31;
  grid.pacific_hfacs(ixatlf,iyatlf,:) = NaN;
  ixgm = 66:73; % Gulf of Mexico
  iygm = 25:30;
  grid.pacific_hfacs(ixgm,iygm,:) = NaN;
  ixhb = 68:73; % Hudson Bay and Baffin Bay
  iyhb = 36:40;
  grid.pacific_hfacs(ixhb,iyhb,:) = NaN;
  ixsop = 1:grid.nx; % southern ocean
  iysop = 1:12;
  grid.pacific_hfacs(ixsop,iysop,:) = NaN;
  ixio = 27:32; % rest of indian ocean
  iyio = 13:19;
  grid.pacific_hfacs(ixio,iyio,:) = NaN;
  %pacific_hfacs(33:34,12,:) = NaN; % two points south of Australia
  
  % Indian ocean
  grid.indic_hfacs = grid.hfacs;
  ixsoi = 1:90; % southern ocean
  iysoi = 1:12;
  grid.indic_hfacs(ixsoi,iysoi,:) = NaN;
  ixap = [1:4 33:90]; % atlantic and pacific
  iyap = 1:40;
  grid.indic_hfacs(ixap,iyap,:) = NaN;
  iynps = 35:40; % north polar sea
  grid.indic_hfacs(:,iynps,:) = NaN;
  ixmeds = 1:10;  % mediterraenan
  iymeds = 29:31;
  grid.indic_hfacs(ixmeds,iymeds,:) = NaN;
  ixscs = 27:32; % south china sea
  iyscs = 21:29;
  grid.indic_hfacs(ixscs,iyscs,:) = NaN;
  grid.indic_hfacs(5,39,:) = NaN; % singular point near Murmansk
  
  % Southern Ocean
  grid.so_hfacs = change( ...
      change(grid.atlantic_hfacs,'==',NaN,-1) ...
      + change(grid.pacific_hfacs,'==',NaN,-1) ...
      + change(grid.indic_hfacs,'==',NaN,-1) ...
      + grid.hfacs-grid.hfacs,'==',-1,NaN);
  iso = find(~isnan(grid.so_hfacs));
  grid.so_hfacs(iso)= grid.hfacs(iso);

  if exist('plotit','var')
    figure;
    spy(isnan(grid.hfacs(:,:,1)'),'kx');
    hold on; 
    spy(~isnan(grid.atlantic_hfacs(:,:,1)')); 
    spy(~isnan(grid.indic_hfacs(:,:,1)'),'g.'); 
    spy(~isnan(grid.so_hfacs(:,:,1)'),'r.'); 
    spy(~isnan(grid.pacific_hfacs(:,:,1)'),'y.');
    title('hfacs')
    axis xy
  end
  
  % masks for zonal averages of C-properties:
  % hfacc
  
  % Atlantic mask
  grid.atlantic_hfacc = grid.hfacc;
  ixe = 6:65; % indian ocean
  iye = 1:28;
  grid.atlantic_hfacc(ixe,iye,:) = NaN;
  ixep = 26:65; % eastern pacific
  iyep = 1:grid.ny; 
  grid.atlantic_hfacc(ixep,iyep,:) = NaN;
  ixw = 6:73; % western pacific
  iyw = 1:23;
  grid.atlantic_hfacc(ixw,iyw,:) = NaN;
  ixgl = 66:68; % galapagos
  iygl = 24;
  grid.atlantic_hfacc(ixgl,iygl,:) = NaN;
  ixch = 73; % coast of Chile
  iych = 13:15;
  grid.atlantic_hfacc(ixch,iych,:) = NaN;
  
  % Pacific mask
  grid.pacific_hfacc = grid.hfacc;
  ixatl = [1:25 74:90]; % most of altantic and indian ocean
  iyatl = 1:grid.ny;
  grid.pacific_hfacc(ixatl,iyatl,:) = NaN;
  ixatlf = 73; % small fix
  iyatlf = 31;
  grid.pacific_hfacc(ixatlf,iyatlf,:) = NaN;
  ixgme = 69:73; % Gulf of Mexico
  iygme = 24:30;
  ixgmw = 66:69;
  iygmw = 26:28;
  grid.pacific_hfacc(ixgme,iygme,:) = NaN;
  grid.pacific_hfacc(ixgmw,iygmw,:) = NaN;
  ixhbb = 68:73; % Hudson and Baffin Bay
  iyhbb = 35:40;
  grid.pacific_hfacc(ixhbb,iyhbb,:) = NaN;
  ixio = 26:34; % rest of indian ocean
  iyio =  1:19;
  grid.pacific_hfacc(ixio,iyio,:) = NaN;
%  grid.pacific_hfacc(33:34,12,:) = NaN; % two points south of Australia
  
  % Indian ocean
  grid.indic_hfacc = ...
      change(change(grid.atlantic_hfacc,'==',NaN,0),'>',0,NaN) ...
      + change(change(grid.pacific_hfacc,'==',NaN,0),'>',0,NaN) ...
      + grid.hfacc;
  
  if exist('plotit','var')
    figure;
    spy(isnan(grid.hfacc(:,:,1)'),'kx'); axis xy
    hold on; 
    spy(~isnan(grid.atlantic_hfacc(:,:,1)')); 
    spy(~isnan(grid.pacific_hfacc(:,:,1)'),'r.');
    spy(~isnan(grid.indic_hfacc(:,:,1)'),'g.'); 
    axis xy
    title('hfacc');

  end
  
  return
1	jmc	1.1	function grid = mit_oceanmasks(grid,plotit)
2			% function grid = mit_oceanmasks(grid)
3			% create masks for the 4x4 degree run for the individual oceans
4			% hfacc/s/w have to be available for this
5
6			% $Header: $
7			% $Name: $
8
9			if grid.nx ~= 90 \| grid.ny ~= 40
10			error(sprintf( ...
11			'%s\n%s\n%s\n%s', ...
12			'Your horizontal grid has dimensions different from 90x40.', ...
13			'mit_oceanmasks creates masks for the individual ocean basins,', ...
14			'therefore it needs to be customized if you change the', ...
15			'horizontal dimensions.'))
16			end
17
18			% masks for the zonal velocities:
19			% hfacw
20
21			% Atlantic mask
22			grid.atlantic_hfacw = grid.hfacw;
23			ixsow = 6:73; % southern ocean and most of pacific
24			iysow = 1:24;
25			grid.atlantic_hfacw(ixsow,iysow,:) = NaN;
26			ixnpw = 15:67; % northern pacific
27			iynpw = 25:37;
28			grid.atlantic_hfacw(ixnpw,iynpw,:) = NaN;
29			% Pacific mask
30			grid.pacific_hfacw = grid.hfacw;
31			ixaw = 74:90; % most of the atlantic
32			iyaw = 1:40;
33			grid.pacific_hfacw(ixaw,iyaw,:) = NaN;
34			ixnaw = 68:73; % north atlantic and gulf of Mexico
35			iynaw = 27:40;
36			grid.pacific_hfacw(ixnaw,iynaw,:) = NaN;
37			ixcaw = 70:73; % carribean
38			iycaw = 24:26;
39			grid.pacific_hfacw(ixcaw,iycaw,:) = NaN;
40			ixiw = 1:27; % most of the indic
41			iyiw = 1:40;
42			grid.pacific_hfacw(ixiw,iyiw,:) = NaN;
43			ixauw = 28:35; % water south and west of Australia
44			iyauw = 1:20;
45			grid.pacific_hfacw(ixauw,iyauw,:) = NaN;
46
47
48			% masks for the meridional overturning stream functions:
49			% hfacs
50
51			% Atlantic mask
52			grid.atlantic_hfacs = grid.hfacs;
53			ixso = 5:73; % southern ocean
54			iyso = 1:24;
55			grid.atlantic_hfacs(ixso,iyso,:) = NaN;
56			ixnp = 14:65; % northern pacific
57			iynp = 25:36;
58			grid.atlantic_hfacs(ixnp,iynp,:) = NaN;
59			ixsoc = 1:90; % complete southern ocean
60			iysoc = 1:12; % south of latg(12) = -34 degN;
61			grid.atlantic_hfacs(ixsoc,iysoc,:) = NaN;
62
63			% Pacific mask
64			grid.pacific_hfacs = grid.hfacs;
65			ixatl = [1:26 74:90]; % most of altantic and indian ocean
66			iyatl = 1:grid.ny;
67			grid.pacific_hfacs(ixatl,iyatl,:) = NaN;
68			ixatlf = 73; % small fix
69			iyatlf = 31;
70			grid.pacific_hfacs(ixatlf,iyatlf,:) = NaN;
71			ixgm = 66:73; % Gulf of Mexico
72			iygm = 25:30;
73			grid.pacific_hfacs(ixgm,iygm,:) = NaN;
74			ixhb = 68:73; % Hudson Bay and Baffin Bay
75			iyhb = 36:40;
76			grid.pacific_hfacs(ixhb,iyhb,:) = NaN;
77			ixsop = 1:grid.nx; % southern ocean
78			iysop = 1:12;
79			grid.pacific_hfacs(ixsop,iysop,:) = NaN;
80			ixio = 27:32; % rest of indian ocean
81			iyio = 13:19;
82			grid.pacific_hfacs(ixio,iyio,:) = NaN;
83			%pacific_hfacs(33:34,12,:) = NaN; % two points south of Australia
84
85			% Indian ocean
86			grid.indic_hfacs = grid.hfacs;
87			ixsoi = 1:90; % southern ocean
88			iysoi = 1:12;
89			grid.indic_hfacs(ixsoi,iysoi,:) = NaN;
90			ixap = [1:4 33:90]; % atlantic and pacific
91			iyap = 1:40;
92			grid.indic_hfacs(ixap,iyap,:) = NaN;
93			iynps = 35:40; % north polar sea
94			grid.indic_hfacs(:,iynps,:) = NaN;
95			ixmeds = 1:10; % mediterraenan
96			iymeds = 29:31;
97			grid.indic_hfacs(ixmeds,iymeds,:) = NaN;
98			ixscs = 27:32; % south china sea
99			iyscs = 21:29;
100			grid.indic_hfacs(ixscs,iyscs,:) = NaN;
101			grid.indic_hfacs(5,39,:) = NaN; % singular point near Murmansk
102
103			% Southern Ocean
104			grid.so_hfacs = change( ...
105			change(grid.atlantic_hfacs,'==',NaN,-1) ...
106			+ change(grid.pacific_hfacs,'==',NaN,-1) ...
107			+ change(grid.indic_hfacs,'==',NaN,-1) ...
108			+ grid.hfacs-grid.hfacs,'==',-1,NaN);
109			iso = find(~isnan(grid.so_hfacs));
110			grid.so_hfacs(iso)= grid.hfacs(iso);
111
112			if exist('plotit','var')
113			figure;
114			spy(isnan(grid.hfacs(:,:,1)'),'kx');
115			hold on;
116			spy(~isnan(grid.atlantic_hfacs(:,:,1)'));
117			spy(~isnan(grid.indic_hfacs(:,:,1)'),'g.');
118			spy(~isnan(grid.so_hfacs(:,:,1)'),'r.');
119			spy(~isnan(grid.pacific_hfacs(:,:,1)'),'y.');
120			title('hfacs')
121			axis xy
122			end
123
124			% masks for zonal averages of C-properties:
125			% hfacc
126
127			% Atlantic mask
128			grid.atlantic_hfacc = grid.hfacc;
129			ixe = 6:65; % indian ocean
130			iye = 1:28;
131			grid.atlantic_hfacc(ixe,iye,:) = NaN;
132			ixep = 26:65; % eastern pacific
133			iyep = 1:grid.ny;
134			grid.atlantic_hfacc(ixep,iyep,:) = NaN;
135			ixw = 6:73; % western pacific
136			iyw = 1:23;
137			grid.atlantic_hfacc(ixw,iyw,:) = NaN;
138			ixgl = 66:68; % galapagos
139			iygl = 24;
140			grid.atlantic_hfacc(ixgl,iygl,:) = NaN;
141			ixch = 73; % coast of Chile
142			iych = 13:15;
143			grid.atlantic_hfacc(ixch,iych,:) = NaN;
144
145			% Pacific mask
146			grid.pacific_hfacc = grid.hfacc;
147			ixatl = [1:25 74:90]; % most of altantic and indian ocean
148			iyatl = 1:grid.ny;
149			grid.pacific_hfacc(ixatl,iyatl,:) = NaN;
150			ixatlf = 73; % small fix
151			iyatlf = 31;
152			grid.pacific_hfacc(ixatlf,iyatlf,:) = NaN;
153			ixgme = 69:73; % Gulf of Mexico
154			iygme = 24:30;
155			ixgmw = 66:69;
156			iygmw = 26:28;
157			grid.pacific_hfacc(ixgme,iygme,:) = NaN;
158			grid.pacific_hfacc(ixgmw,iygmw,:) = NaN;
159			ixhbb = 68:73; % Hudson and Baffin Bay
160			iyhbb = 35:40;
161			grid.pacific_hfacc(ixhbb,iyhbb,:) = NaN;
162			ixio = 26:34; % rest of indian ocean
163			iyio = 1:19;
164			grid.pacific_hfacc(ixio,iyio,:) = NaN;
165			% grid.pacific_hfacc(33:34,12,:) = NaN; % two points south of Australia
166
167			% Indian ocean
168			grid.indic_hfacc = ...
169			change(change(grid.atlantic_hfacc,'==',NaN,0),'>',0,NaN) ...
170			+ change(change(grid.pacific_hfacc,'==',NaN,0),'>',0,NaN) ...
171			+ grid.hfacc;
172
173			if exist('plotit','var')
174			figure;
175			spy(isnan(grid.hfacc(:,:,1)'),'kx'); axis xy
176			hold on;
177			spy(~isnan(grid.atlantic_hfacc(:,:,1)'));
178			spy(~isnan(grid.pacific_hfacc(:,:,1)'),'r.');
179			spy(~isnan(grid.indic_hfacc(:,:,1)'),'g.');
180			axis xy
181			title('hfacc');
182
183			end
184
185			return