% DERIV.M % function d_dt=deriv(t,posvel) % VARIABLES: % d_dt = a vector of rates of change of x, y, u, v for time, t % = [dx/dt dy/dt du/dt dv/dt] % NOTE: Time, t, is in seconds but is evaluated hourly. % The 5th & 6th elements of d_dt are constants. % t = time in seconds % posvel = a vector containing initial particle position (x0,y0) and % velocity (u0,v0), and the topography amplitude (h0) & % reference latitude (lat) (optional) % = [x0 y0 u0 v0 h0 lat] % % UNITS: % m = meters % km = kilometers % s = seconds % hr = hours % rad = radians % deg = degrees function d_dt=deriv(t,posvel) %---------------------------------------------------------------------------- % Define parameters of the problem % Define amplitude & sign of topography and reference latitude if(length(posvel)==4) h0=.05; % set default bump (if not specified in input) lat=30; % set default latitude (if not specified in input) else h0=posvel(5); % user input or h0 from main program lat=posvel(6); % user input or lat from main program end grav = 9.8; % gravity in [m/s^2] Omega = (2*pi)/(24*3600); % rotation rate of the earth in [rad/s] f = 2*Omega*sin(lat*(pi/180)); % Coriolis parameter in [1/s] rfact=0.05; r=max(5e-6,rfact*abs(f)); % scale friction by f to have units of 1/s %---------------------------------------------------------------------------- % Separate input vector into individual variables for clarity x0=posvel(1); y0=posvel(2); u0=posvel(3); v0=posvel(4); %---------------------------------------------------------------------------- % Define the topography of the sea surface to be a symmetric Gaussian "bump" % All quantities have units of [m] lx=30*10^3; % x length scale of bump ly=30*10^3; % y length scale of bump xm=400*10^3; % x location of the maximum height of bump ym=400*10^3; % y location of the maximum height of bump rfact=.02; % set default Rayleigh friction factor % Compute the height and slope of the topography at the particle position h = h0*exp(-((x0-xm)/lx)^2 -((y0-ym)/ly)^2); hx = -2*h*(x0-xm)/(lx^2); % x-derivative, or slope of topography (dh/dx) hy = -2*h*(y0-ym)/(ly^2); % y-derivative, or slope of topography (dh/dy) %---------------------------------------------------------------------------- % Define the (wind) forcing tfact=3600; % converts hours to seconds tmax=72*tfact; % length of time of forcing (3 days) if(t>0 & t