% evaluate the SEL4 model % by Jim Price, April 98 clear; format compact % a matlab program to set default graphics things set(0,'DefaultTextFontSize',16) set(0,'DefaultAxesFontSize',14) set(0,'DefaultAxesLineWidth',1.4) set(0,'DefaultLineLineWidth',1.2) % for Lotus3 you could use lat = 36., Q = 600., tau = 0.08, H = 50. lat = input(' Enter the latitude (degrees) ') f = 2.*7.292e-5*sin(lat*pi/180.); Q = input(' Enter the noon heat flux (W/m*m) ') Qstar = 9.8*0.3*Q/(1025.*1025.*4100.); Pq = 12.*3.6e3; Ptau = (1/f)*sqrt(2 - 2*cos(f*Pq/2)); tau = input(' Enter the wind stress (Pa) ') % assumed northward ustar = sqrt(tau/1028.); H = input(' Enter the depth of semi-permanent stratification (m) ') Dq = ustar*ustar*Ptau/sqrt(Qstar*Pq/2); % the trapping depth P24 = 8.64e4; psir = (f*Pq - sin(f*Pq))/(f*P24); psii = (1 - cos(f*Pq))/(f*P24); psimag = sqrt(psir^2 + psii^2); Un = ustar^2/(f*H); % the neutral velocity scale Uh = psimag*sqrt(Qstar*Pq/2)/(f*Ptau); % the heating scale ang = atan2(psii, psir); V1 = Un*(1 + (H/Dq - 1)*psimag*exp(i*ang)); V2 = Un*(1 - psimag*exp(i*ang)); V1mag = abs(V1); Echeck = Dq*V1 + (H - Dq)*V2; % the transport Etran = ustar^2/f; % ratio of transports T1 = real(V1)*Dq; T2 = real(V2)*(H-Dq); Trat = T1/T2; T1frac = T1/(T1 + T2) % the Ekman transport figure(1) clf reset axis([-0.02 0.08 -0.04 0.06]) xlabel('acrosswind, m/sec') ylabel('downwind, m/sec') grid arrow6([0 0], [real(V1) imag(V1)], 0.05, 1.9); % layer 1 velocity arrow6([0 0], [real(V2) imag(V2)], 0.03, 1.9); % layer 2 velocity arrow6([0 0], [Un 0], 0.02, 1.9, pi/6, 1); % the neutral velocity text(real(V1), imag(V1)+0.002, 'V1') text(real(V2), imag(V2)-0.004, 'V2') text(Un+0.002, 0., 'Un') % now compute the profile by interpolating z = 0:1:H; z = z'; a = 0*z; az = [0 0.5*Dq 1.5*Dq H]; z = 0:1:H; z = z'; a = 0*z; az = [0 0.5*Dq 1.5*Dq H]; ao = [1 1 0 0]; a = interp1(az, ao, z, 'linear'); % interpolate to find a at all depths Vz = a*V1 + (1-a)*V2; % interpolate to find V at all depths figure(2) clf reset subplot(1,2,1) plot(real(Vz), -z); ylabel('depth, m') xlabel('crosswind, m/sec') axis([-0.02 0.08 -H 0]) grid subplot(1,2,2) plot(imag(Vz), -z); xlabel('downwind, m/sec') axis([-0.02 0.08 -H 0]) grid