0001 function md=icetransient2d(md);
0002
0003
0004
0005
0006
0007 global gridset
0008
0009
0010 fem=struct();
0011
0012
0013 [fem.ishutter,fem.ismacayealpattyn,fem.isstokes]=DiagnosticSolutionType(md.elements_type);
0014
0015
0016 if fem.ishutter,
0017 fem.m_ss=CreateFemModel(md,'surface_slope_compute');
0018 fem.m_dhu=CreateFemModel(md,'diagnostic_hutter');
0019 end
0020 if fem.ismacayealpattyn,
0021 fem.m_dh=CreateFemModel(md,'diagnostic_horiz');
0022 end
0023 m_p=CreateFemModel(md,'prognostic');
0024
0025
0026 solution.u_g=zeros(m_p.gridset.gsize,1);
0027 if ~isnan(md.vx), solution.u_g(1:6:m_p.gridset.gsize)=md.vx/md.yts; end
0028 if ~isnan(md.vy), solution.u_g(2:6:m_p.gridset.gsize)=md.vy/md.yts; end
0029
0030
0031 solution.thickness=zeros(m_p.gridset.gsize,1);solution.thickness(1:6:m_p.gridset.gsize)=md.thickness;
0032 solution.surface=zeros(m_p.gridset.gsize,1);solution.surface(1:6:m_p.gridset.gsize)=md.surface;
0033 solution.bed=zeros(m_p.gridset.gsize,1);solution.bed(1:6:m_p.gridset.gsize)=md.bed;
0034 solution.time=0;
0035
0036
0037 melting=zeros(m_p.gridset.gsize,1);
0038 accumulation=zeros(m_p.gridset.gsize,1);
0039 melting(1:6:end)=md.melting/md.yts;
0040 accumulation(1:6:end)=md.accumulation/md.yts;
0041
0042
0043 fem.inputs=struct();
0044 fem.inputs.accumulation=accumulation;
0045 fem.inputs.melting=melting;
0046 fem.inputs.dt=md.dt;
0047
0048
0049 dt=md.dt;
0050 finaltime=md.ndt;
0051 time=dt;
0052 n=1;
0053
0054 while time<finaltime+dt,
0055
0056 disp(sprintf('\n%s%g%s%g%s%g\n','time [yr]: ',time/md.yts,' iteration number: ',n,'/',floor(finaltime/dt)));
0057
0058 solution(n+1).time=time;
0059
0060
0061 fem.inputs.velocity=solution(n).u_g;
0062 fem.inputs.thickness=solution(n).thickness;
0063 fem.inputs.bed=solution(n).bed;
0064 fem.inputs.surface=solution(n).surface;
0065
0066
0067
0068
0069 solution(n+1).u_g=icediagnostic2d(md,fem);
0070
0071
0072 disp(sprintf('%s',' computing new thickness...'));
0073 fem.inputs.velocity_average=solution(n+1).u_g;
0074 new_thickness=iceprognostic_core(m_p,'prognostic',fem.inputs);
0075
0076
0077 disp(sprintf('%s',' updating geometry...'));
0078 indx=1:6:m_p.gridset.gsize; indx=indx(m_p.tpart);
0079 [new_bed,new_surface,new_thickness]=UpdateGeometry(md,new_thickness(indx),solution(n).bed(indx),solution(n).surface(indx));
0080
0081
0082 solution(n+1).thickness=zeros(m_p.gridset.gsize,1);solution(n+1).thickness(1:6:m_p.gridset.gsize)=new_thickness;
0083 solution(n+1).bed=zeros(m_p.gridset.gsize,1);solution(n+1).bed(1:6:m_p.gridset.gsize)=new_bed;
0084 solution(n+1).surface=zeros(m_p.gridset.gsize,1);solution(n+1).surface(1:6:m_p.gridset.gsize)=new_surface;
0085
0086
0087 disp(sprintf('%s',' checking time stepping...'));
0088 [back,dt,time]=TimeStepping(md,solution,dt,time);
0089 if back,
0090 continue;
0091 end
0092
0093
0094 time=time+dt;
0095 n=n+1;
0096
0097 end
0098
0099
0100 indx=1:6:m_p.gridset.gsize; indx=indx(m_p.tpart);
0101 indy=2:6:m_p.gridset.gsize; indy=indy(m_p.tpart);
0102 for i=1:length(solution),
0103 solution2(i).vx=solution(i).u_g(indx)*md.yts;
0104 solution2(i).vy=solution(i).u_g(indy)*md.yts;
0105 solution2(i).vel=sqrt(solution2(i).vx.^2+solution2(i).vy.^2);
0106 solution2(i).thickness=solution(i).thickness(indx);
0107 solution2(i).surface=solution(i).surface(indx);
0108 solution2(i).bed=solution(i).bed(indx);
0109 solution2(i).time=solution(i).time/md.yts;
0110 end
0111 md.transient_results=solution2;