function pwpgo(qi,qo,emp,tx,ty,m)

% The following 3 paragraphs are comments from the original PWP Fortran version.
% The entrainment parameterization mentioned in the third paragraph
% is not implemented in this or the URI PWP v1.4 versions of the model. 
% Similarly, the parameterization of inertial-internal wave dispersion included 
% in the code (lines 106-112 below) is not implemented. (Fernando Santiago-Mandujano, 6/2016)
% 
%  This subroutine is an implementation of the Price, Weller,
%  Pinkel upper ocean model described in JGR 91, C7 8411-8427
%  July 15, 1986 (PWP). This version was coded and documented by
%  Jim Price in April 1989.

%  Edited on 20 September 1993 by JFP to allow for a critical
%  gradient Richardson number other than 1/4, and to implement a
%  different and a priori better means of achieving convergence
%  of gradient ri mixing (see subroutine stir for details).
%  The major difference is that the revised scheme gives a more
%  smoothly varying mixed layer depth over a diurnal cycle.
%  Edited on 14 December 1998 by JFP to clarify the bulk Ri
%  relaxation.

%  This model also implements an energy budget form of an
%  entrainment parameterization that is very similar to that
%  described in Price, Mooers, Van Leer, JPO 8, 4, 582-599 (and
%  references therein). This part of the model should be
%  treated as developmental only, as there are several features
%  that are arbitrary to this model (i.e., the depth of the ml
%  during times of heating can be the grid interval, dz). To
%  use this parameterization set the bulk Richardson number to
%  zero, and set em1, or em2, or em3 to non-zero. The gradient
%  Richardson number can be zero or not.

global dt dz days depth dt_save lat g cpw rb rg beta1 beta2 udrag Method InstabilityMatrix richard_grad bulk
global nz z T s d 
global u v absrb
global f ucon
global ml_depth
global kcount;   % gradient Richardson number iteration counter (FS-M, 4/6/09)

%  Apply heat and fresh water fluxes to the top most grid cell.

T(1) = T(1)+(qi*absrb(1)-qo)*dt./(dz*d(1)*cpw);
s(1) = s(1)/(1-emp*dt/dz);

%  Absorb solar radiation at depth.

T(2:nz) = T(2:nz)+qi*absrb(2:nz)*dt./(dz*d(2:nz)*cpw);

%  Compute the density, and relieve static instability, if it occurs.

d = sw_dens0(s,T);

ml_index2=find(diff(d)<0);    % index of the static instability for diagnostic
InstabilityMatrix(ml_index2,m)=1;

% Remove static instabilities

remove_si;

%  At this point the density profile should be statically stable.

%  Find the index of the surface mixed-layer right after the heat/salt fluxes.

if Method==1;
        ml_index = find(diff(d)>0.0001*dz,1,'first');% density gradient difference 
                                                 % method :1*10^-4 (kg.m-3).m-1
elseif Method==2
        ml_index = 5/dz + find(d(round(5/dz):end)-d(round(5/dz))>0.01,1,'first'); 
            % Delta criteria density method 5m reference : delta = 0.01 kg.m-3
elseif Method==3
        ml_index = find(abs(diff(T))>0.02*dz,1,'first'); % temperature gradient 
                                        % difference method : 0.02 degreeC.m-1
elseif Method==4
        ml_index = find(-T+T(1)>0.8,1,'first'); % Delta criteria temperature method 
                                      %surface reference : delta = 0.8 degree C
end
    
if isempty(ml_index)
	error_text = ['Error reached in pwpgo.m: the mixed layer is =too deep!' sprintf('\n') ...
	'Possible reasons:' sprintf('\n') ...
	'1) tempearature inversion at surface is too strong' sprintf('\n') ...
	'2) wind is too strong' sprintf('\n') ...
	'3) evaporation is too high'];
	error([ sprintf('\n') error_text])
end

%  Get the depth of the surface mixed-layer.
ml_depth = z(ml_index+1);

%  Time step the momentum equation.

%  Rotate the current throughout the water column through an
%  angle equal to inertial rotation for half a time step.

ang = -f*dt/2;
rot(ang);

%  Apply the wind stress to the mixed layer as it now exists.
du = (tx/(ml_depth*d(1)))*dt;
dv = (ty/(ml_depth*d(1)))*dt;
u(1:ml_index) = u(1:ml_index)+du;
v(1:ml_index) = v(1:ml_index)+dv;


%  Apply drag to the current (this is a horrible parameterization of
%  inertial-internal wave dispersion).
ucon = 0; % (This parameterization is not implemented. FS-M, 6/2016)
if ucon > 1E-10
	u = u*(1-dt*ucon);
	v = v*(1-dt*ucon);
end

%  Rotate another half time step.
rot(ang);

%  Finished with the momentum equation for this time step.

%  Do the bulk Richardson number instability form of mixing (as in PWP).
if rb > 1E-5
	bulk_mix(ml_index,m)
end

%  Do the gradient Richardson number instability form of mixing.
if rg > 0
	grad_mix(m);
    disp(['m = ',num2str(m),'  kcount = ',num2str(kcount)])
end
