wim.m 6.68 KB
function [T,cice,Hsig,Ei,E,S_win,S_wcp,S_ice,Dmax,Dave,om_do] = wim(C,h,D,U10,Tm,Hs)
% WIM - Models the propagation of a wave spectrum along a 1D 
% ice-covered transect.
%
% Syntax: [T,cice,Hsig,Ei,E,S_win,S_wcp,S_ice] = wim_final(C,h,D,U10,Tm,Hs)
%
% Inputs:
%    C    : ice concentration. For a homogeneous concentration, C is scalar
%           in tenths. Otherwise, C is a vector of the size of the spatial
%           grid, in fraction of 1 (0 < C < 1). In this version of WIM, if
%           C is homogeneous, the transect is 5km long with spatial
%           resolution of 500m.
%    h    : ice thickness in meters. h is a scalar.
%    D    : average floe diameter in meters. D is a scalar.
%    U10  : wind speed at 10m elevation in m/sec. U10 is a scalar.
%    Tm   : mean wave period to build the spectrum, in seconds. Tm is a scalar.
%    Hs   : significant wave height, in meters. Hs is a scalar.
%
% Outputs:
%    T    : wave period in seconds. T is a vector the size of the frequency range.
%    cice : ice concentration vector in fraction of 1. cice has the size of
%           the spatial grid.If C is a vector, then cice = C.
%    Hsig : final wave heigh in each cell in meters. Hsig has the size of
%           the spatial grid.
%    Ei   : initial spectrum in m2/hz.
%    E    : final wave spectrum in each cell [m2/Hz].
%    S_win: final wind source term in each cell in m2/Hz/sec.
%    S_wcp: final white-capping source term in each cell in m2/Hz/sec.
%    S_ice: final ice source term in each cell in m2/Hz/sec.

% Parameters
g  = 9.81;  % gravitational acceleration
dx = 5000;    % spatial resolution

% Ice conditions
if length(C) == 1    % Homogeneous ice concentration case
    nx   = 10;                   % 10 cells of 500m -> 5km long transect
    cice = zeros(1,nx)+C;     % homogeneous ice concentration vector
    hice = zeros(1,nx)+h;        % homogeneous ice thickness vector
    Dmax = zeros(1,nx)+D;        % homogeneous floe diameter vector
    Dave = zeros(1,nx)+D;  
else                 % Variable ice concentration case
    nx            = length(C);   % number of cells
    cice          = C;           % ice concentration vector
    hice          = zeros(1,nx); % ice thickness vector size allocation
    hice(find(C)) = h;           % ice thickness vector
    Dmax          = zeros(1,nx); % floes diameter vector size allocation
    Dmax(find(C)) = D;           % floes diameter vector
end

% Waves
fmin  = 1/20;                   % minimum wave frequency
fmax  = 1/2.5;                  % maximum wave frequency
om1   = 2*pi*fmin;              % minimum wave radial frequency
om2   = 2*pi*fmax;              % maximum wave radial fequency
nw    = 61;                     % number of frequency bins
dw    = (om2-om1)/(nw-1);       % integral interval for wave radial frequencies
om    = om1+(0:nw-1)'*dw;       % wave radial frequencies vector
T     = 2*pi./om;               % wave periods vector
Ei    = jonswap(Tm,Hs,om);      % JONSWAP spectrum
wlng  = g.*T.^2./(2.*pi);       % wavelength as a function only of wave period
cp    = sqrt(g.*wlng./(2.*pi)); % phase speed
cg    = cp./2;                  % group speed
cgmax = max(cg);                % group speed maximum
cg(:) = cgmax;                  % no dispersion = all group speed are the same (maximum)

% Wind
U10 = repmat(U10,nx,1); % wind speed at 10m elevation vector

% Temporal grid
dt     = dx/cgmax;       % time interval (temporal resolution)
nsteps = nx+1;           % number of time steps (=nx+1 because advection is 
                         % done from one cell at each time step)
time   = 0:dt:nsteps*dt; % time vector
nt     = length(time);   % time vector size

% Memory preallocation
E        = zeros(nx,nw); % wave spectrum
Swin     = zeros(nx,nw); % wind source term
Swcp     = zeros(nx,nw); % white-capping source term
Sice     = zeros(nx,nw); % ice dissipation source term
S_win    = zeros(nx,nw); % wind source term weighted by open water fraction (for outputs)
S_ice    = zeros(nx,nw); % ice dissipation source terme weighted by ice fraction
S_wcp    = zeros(nx,nw); % white-capping source term weighted by open water fraction
Hsig     = zeros(1,nx);  % significant wave height

% Action density limiter
Slim = action_density_limiter(om,cp); % SEE action_density_limiter routine 
                                     

for n=1:nt     % Time loop
    % Advection.
    % The advection is done by solving Dt(S) = 0 using the
    % Lax-Wendroff scheme with Superbee flux limiting and a Neumann
    % boundary condition. The advection is performed over the whole domain
    % in one step on an unattenuated intermediate spectrum.
    for w=1:nw % Advection loop over each frequency
        E(:,w)   = advection(E(:,w),cg(w,:),dx,dt); % SEE advection routine
    end
    % Incident wave spectrum
    E(1,:) = Ei; % The initial wave spectrum is forced in the first cell
    % Processes for each spatial cell with ice (generation by wind, dissipation 
    % by white-capping and attenuation by ice)
    for i=2:nx % Spatial loop
%         % Generation by wind
%         Swin(i,:) = wind_gen(U10(i),E(i,:),om,cp); % SEE wind_gen routine
%         Swin(i,:) = min(Swin(i,:),Slim');           % action density limiter
%         E(i,:) = E(i,:) + Swin(i,:)*dt*(1-cice(i)); % the wave spectrum is updated (explicit scheme)
%         S_win(i,:) = Swin(i,:)*(1-cice(i));         % effective wind source term for outputs
%         % Dissipation by white-capping
%         Swcp(i,:) = white_cap(E(i,:),om,cp); % SEE white-cap routine
%         Swcp(i,:) = max(Swcp(i,:),-Slim');          % action density limiter
%         E(i,:) = E(i,:) + Swcp(i,:)*dt*(1-cice(i)); % the wave spectrum is updated (semi-implicit scheme)
%         S_wcp(i,:) = Swcp(i,:)*(1-cice(i));         % effective white-capping source term for outputs
      
        % Attenuation by ice
        if cice(i)>0
            Sice(i,:) = ice_att(E(i,:),T,cg,cice(i),hice(i),Dave(i),dt*cice(i)); % SEE ice_att routine
            E(i,:) = E(i,:) + Sice(i,:)*dt; % the wave spectrum is updated (implicit scheme)
        else
            Sice(i,:) = 0;                          % if there is no ice, there is no attenuation Sherlock !
        end
        S_ice(i,:) = Sice(i,:)*cice(i);             % effective ice dissipation source term for outputs
        [Dmax(i),om_do(i)]=floe_breaking(E(i,:),h,om,Dmax(i),dw);
           Dave(i)=FSD(Dmax(i),D);
        % Wave spectrum statistics
        E(E<0) = 0;                                 % avoiding negative energy
        m0 = trapz(om,E(i,:));                      % total energy
        Hsig(i) = 4*sqrt(m0);                       % significant wave height
    end
end

end