source: issm/trunk/src/m/solutions/ice/transient2d.m@ 744

function md=transient2d(md);
%ICETRANSIENT2D - transient solution of 2d models
%
%   Usage:
%      md=transient2d(md)

global gridset

%Create fem structure (input of diagnostic2d)
fem=struct();

%Figure out which type of elements are present
fem.ishutter=md.ishutter;
fem.ismacayealpattyn=md.ismacayealpattyn;
fem.isstokes=md.isstokes;

%First, build elements,grids,loads, etc ... for horizontal model
if fem.ishutter,
        fem.m_ss=CreateFemModel(md,'surface_slope_compute');
        fem.m_dhu=CreateFemModel(md,'diagnostic_hutter');
end
if fem.ismacayealpattyn,
        fem.m_dh=CreateFemModel(md,'diagnostic_horiz');
end
m_p=CreateFemModel(md,'prognostic');

%initialize (velocity,pressure)
solution.u_g=zeros(m_p.gridset.gsize,1);
if ~isnan(md.vx), solution.u_g(1:6:m_p.gridset.gsize)=md.vx/md.yts; end
if ~isnan(md.vy), solution.u_g(2:6:m_p.gridset.gsize)=md.vy/md.yts; end

%initialize geometry and time
solution.thickness=zeros(m_p.gridset.gsize,1);solution.thickness(1:6:m_p.gridset.gsize)=md.thickness;
solution.surface=zeros(m_p.gridset.gsize,1);solution.surface(1:6:m_p.gridset.gsize)=md.surface;
solution.bed=zeros(m_p.gridset.gsize,1);solution.bed(1:6:m_p.gridset.gsize)=md.bed;
solution.time=0;

%melting and accumulation are taken constant for all iterations
melting=zeros(m_p.gridset.gsize,1);
accumulation=zeros(m_p.gridset.gsize,1);
melting(1:6:end)=md.melting/md.yts;           %from m/yr to m/s
accumulation(1:6:end)=md.accumulation/md.yts; %from m/yr to m/s

%build constant inputs
fem.inputs=struct();
fem.inputs.accumulation=accumulation;
fem.inputs.melting=melting;
fem.inputs.dt=md.dt;

%first time step is given by model. 
dt=md.dt;
finaltime=md.ndt;
time=dt;
n=1; %counter

while  time<finaltime+dt, %make sure we run up to finaltime.

        disp(sprintf('\n%s%g%s%g%s%g\n','time [yr]: ',time/md.yts,'    iteration number: ',n,'/',floor(finaltime/dt)));
        
        solution(n+1).time=time;

        %update inputs
        fem.inputs.velocity=solution(n).u_g;
        fem.inputs.thickness=solution(n).thickness;
        fem.inputs.bed=solution(n).bed;
        fem.inputs.surface=solution(n).surface;

        %Deal with velocities.

        %Get horizontal solution. 
        solution(n+1).u_g=diagnostic2d(md,fem);
        %remove the pressure
        solution(n+1).u_g(4:6:end)=0;

        %compute new thickness
        disp(sprintf('%s','   computing new thickness...'));
        fem.inputs.velocity_average=solution(n+1).u_g;
        new_thickness=prognostic_core(m_p,'prognostic',fem.inputs);

        %update surface and bed using the new thickness
        disp(sprintf('%s','   updating geometry...'));
        indx=1:6:m_p.gridset.gsize; indx=indx(m_p.tpart);
        [new_bed,new_surface,new_thickness]=UpdateGeometry(md,new_thickness(indx),solution(n).bed(indx),solution(n).surface(indx));
        
        %Record bed surface and thickness in the solution
        solution(n+1).thickness=zeros(m_p.gridset.gsize,1);solution(n+1).thickness(1:6:m_p.gridset.gsize)=new_thickness;
        solution(n+1).bed=zeros(m_p.gridset.gsize,1);solution(n+1).bed(1:6:m_p.gridset.gsize)=new_bed;
        solution(n+1).surface=zeros(m_p.gridset.gsize,1);solution(n+1).surface(1:6:m_p.gridset.gsize)=new_surface;

        %figure out if time stepping is good
        disp(sprintf('%s','   checking time stepping...'));
        [back,dt,time]=TimeStepping(md,solution,dt,time);
        if back,
                continue;
        end

        %update time and counter
        time=time+dt;
        n=n+1;

end

%load results onto model
indx=1:6:m_p.gridset.gsize; indx=indx(m_p.tpart);
indy=2:6:m_p.gridset.gsize; indy=indy(m_p.tpart);
for i=1:length(solution),
        solution2(i).vx=solution(i).u_g(indx)*md.yts;
        solution2(i).vy=solution(i).u_g(indy)*md.yts;
        solution2(i).vel=sqrt(solution2(i).vx.^2+solution2(i).vy.^2);
        solution2(i).thickness=solution(i).thickness(indx);
        solution2(i).surface=solution(i).surface(indx);
        solution2(i).bed=solution(i).bed(indx);
        solution2(i).time=solution(i).time/md.yts;
        solution2(i).step=i;
end
md.results.transient=solution2;