31 views (last 30 days)

Hey!

I'm trying to remove points from my mesh grid in order to get a graph similar to this one:

Basically I want to remove all the points below the line given in the code below by h1.

clear all; close all; clc % clear variables, close plots, clear screen

%% Constants

D = 60; % Ekman layer depth (m)

rho_0 = 1025; % density (kg/m^3)

f = 0.545/(10^4); % Coriolis parameter

c = (1+i)*pi/D; % c parameter

ug = 0; % zonal geostrophic wind speed

tau_x = 0; % x-component of wind shear stress

tau_y = 0.07; % y-component of wind shear stress

tau = tau_x + i*tau_y; % complex wind shear stress

u0 = (1-i)*pi*tau/(rho_0*f*D); % surface Ekman velocity in the deep ocean limit

hDmax = 2.5; % maximum h/D value

%% Bathymetry

h1 = @(x) -x/1000; % bathymetry profile function

x_grid = linspace(-200000,0,1000); % linearly spaced grid for the horizontal coordinate (offshore distance, x)

z_grid = linspace(-100,0,1000); % linearly spaced grid for the vertical coordinate (depth, z)

h = h1(x_grid); % vector with bathymetry

tmp = find(h/D>hDmax); % reduce offshore grid

xindex = tmp(end);

h = h(xindex:end);

L = length(h);

N = length(z_grid);

zeta = pi*h/D; % vector with zeta function values

%% Horizontally-varying structure functions

alpha = (cosh(zeta).*cos(zeta)).^2 + (sinh(zeta).*sin(zeta)).^2;

S1 = cosh(zeta).*cos(zeta)./alpha;

S2 = sinh(zeta).*sin(zeta)./alpha;

T1 = cosh(zeta).*sinh(zeta)./alpha;

T2 = cos(zeta).*sin(zeta)./alpha;

%% Velocities, transport and stream function

vg = (2*pi*tau_y/(rho_0*f*D)).*(1-S1)./(T1-T2) - ug.*(T1+T2-2.*pi.*h./D)./(T1-T2); % meridional geostrophic wind speed

ugu = h*0 + ug; % adjusting the size of the vector

ugc = complex(ugu,vg); % complex total geostrophic velocity

Uek = abs(tau_y/(rho_0*f)); % Ekman transport

psi = zeros(N,L);

psi = complex(psi,psi);

for jj=1:L

for ii=1:N

psi(ii,jj) = (1/c).*(u0.*(1-(cosh(c.*(z_grid(ii)+h(jj))))./(cosh(c.*h(jj)))) +...

ugc(jj).*sinh(c.*z_grid(ii))./cosh(c.*h(jj)));

end

end

psi = real(psi); % real part of stream function

%% Plotting

X_grid = repmat(-h/D,[N 1]);

Z_grid = repmat(z_grid/D, [L 1]);

Z_grid = Z_grid'

Z = psi/Uek;

contour(X_grid,Z_grid,Z,'ShowText','on','LevelStep',0.1)

The figure produced at the moment is the following:

What can I do?

Sign in to answer this question.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.