%  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