0001 function md=icetransient3d(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.ismacayealpattyn,
0017 fem.m_dh=CreateFemModel(md,'diagnostic_horiz');
0018 end
0019 if fem.ishutter,
0020 fem.m_ss=CreateFemModel(md,'surface_slope_compute');
0021 fem.m_dhu=CreateFemModel(md,'diagnostic_hutter');
0022 end
0023 if fem.isstokes,
0024 fem.m_bs=CreateFemModel(md,'bed_slope_compute');
0025 fem.m_ds=CreateFemModel(md,'diagnostic_stokes');
0026 end
0027 fem.m_dbv=CreateFemModel(md,'diagnostic_basevert');
0028 fem.m_dv=CreateFemModel(md,'diagnostic_vert');
0029
0030
0031 fem.m_t=CreateFemModel(md,'thermaltransient');
0032 fem.m_m=CreateFemModel(md,'meltingtransient');
0033 fem.m_p=CreateFemModel(md,'prognostic');
0034
0035
0036 solution.u_g=zeros(fem.m_dhu.gridset.gsize,1);
0037 if ~isnan(md.vx), solution.u_g(1:6:fem.m_dhu.gridset.gsize)=md.vx/md.yts; end
0038 if ~isnan(md.vy), solution.u_g(2:6:fem.m_dhu.gridset.gsize)=md.vy/md.yts; end
0039 if ~isnan(md.vz), solution.u_g(3:6:fem.m_dhu.gridset.gsize)=md.vz/md.yts; end
0040 solution.pressure=zeros(fem.m_dhu.gridset.gsize,1); solution.pressure(1:6:fem.m_dhu.gridset.gsize)=md.rho_ice*md.g*(md.surface-md.z);
0041 solution.t_g=zeros(fem.m_t.gridset.gsize,1);
0042 if ~isempty(md.temperature),
0043 solution.t_g(1:6:fem.m_t.gridset.gsize)=md.temperature;
0044 else
0045 solution.t_g(1:6:fem.m_t.gridset.gsize)=md.meltingpoint;
0046 end
0047 pos=find(isnan(solution.t_g));
0048 solution.t_g(pos)=md.meltingpoint;
0049 if ~isempty(md.melting), solution.melting_g(1:6:fem.m_m.gridset.gsize)=md.melting;end;
0050
0051
0052 solution.mesh=struct('elements',md.elements,'grids',[md.x md.y md.z]);
0053 solution.thickness=zeros(fem.m_dhu.gridset.gsize,1);solution.thickness(1:6:fem.m_dhu.gridset.gsize)=md.thickness;
0054 solution.surface=zeros(fem.m_dhu.gridset.gsize,1);solution.surface(1:6:fem.m_dhu.gridset.gsize)=md.surface;
0055 solution.bed=zeros(fem.m_dhu.gridset.gsize,1);solution.bed(1:6:fem.m_dhu.gridset.gsize)=md.bed;
0056 solution.time=0;
0057
0058
0059 fem.inputs=struct();
0060
0061 fem.inputs.melting=zeros(fem.m_m.gridset.gsize,1); melting=ShiftLayers(md,md.melting,1,1); fem.inputs.melting(1:6:fem.m_m.gridset.gsize)=melting;
0062 fem.inputs.accumulation=zeros(fem.m_m.gridset.gsize,1); accumulation=ShiftLayers(md,md.accumulation,md.numlayers,1); fem.inputs.accumulation(1:6:fem.m_m.gridset.gsize)=accumulation;
0063 fem.inputs.dt=md.dt;
0064
0065
0066 dt=md.dt;
0067 finaltime=md.ndt;
0068 time=dt;
0069 n=1;
0070
0071 while time<finaltime+dt,
0072
0073 disp(sprintf('\n%s%g%s%g%s%g\n','time [yr]: ',time/md.yts,' iteration number: ',n,'/',floor(finaltime/dt)));
0074
0075 solution(n+1).time=time;
0076
0077
0078 fem.inputs.velocity=solution(n).u_g;
0079 fem.inputs.pressure=solution(n).pressure;
0080 fem.inputs.temperature=full(solution(n).t_g);
0081 fem.inputs.thickness=solution(n).thickness;
0082 fem.inputs.bed=solution(n).bed;
0083 fem.inputs.surface=solution(n).surface;
0084
0085
0086
0087
0088 solution(n+1).t_g=solution(1).t_g;
0089
0090
0091
0092
0093 solution(n+1).melting_g=solution(1).melting_g;
0094
0095
0096
0097
0098
0099
0100
0101
0102 disp(sprintf('%s',' computing velocities...'));
0103 solution(n+1).u_g=icediagnostic3d(md,fem);
0104 solution(n+1).pressure=zeros(fem.m_dv.gridset.gsize,1); solution(n+1).pressure(1:6:end)=solution(n+1).u_g(4:6:end);
0105
0106
0107 disp(sprintf('\n%s',' computing new thickness...'));
0108
0109 velocity_average=HorizontalVelocityDepthAverage(md,solution(n+1).u_g);
0110 fem.inputs.velocity_average=velocity_average;
0111 new_thickness=iceprognostic_core(fem.m_p,'prognostic',fem.inputs);
0112
0113
0114 indx=1:6:fem.m_p.gridset.gsize; indx=indx(fem.m_p.tpart);
0115 new_thickness=project3d(md,project2d(md,new_thickness(indx),1),'node');
0116
0117
0118 disp(sprintf('%s',' updating geometry...'));
0119 [new_bed,new_surface,new_thickness]=UpdateGeometry(md,new_thickness,solution(n).bed(indx),solution(n).surface(indx));
0120
0121
0122 [solution(n+1).mesh solution(n+1).u_g solution(n+1).t_g solution(n+1).pressure]=UpdateMesh(md,fem.m_dv.grids,solution(n).mesh,solution(n+1).u_g,solution(n+1).t_g,solution(n+1).pressure,new_bed,new_thickness);
0123 solution(n+1).thickness=zeros(fem.m_p.gridset.gsize,1);solution(n+1).thickness(1:6:fem.m_p.gridset.gsize)=new_thickness;
0124 solution(n+1).bed=zeros(fem.m_p.gridset.gsize,1);solution(n+1).bed(1:6:fem.m_p.gridset.gsize)=new_bed;
0125 solution(n+1).surface=zeros(fem.m_p.gridset.gsize,1);solution(n+1).surface(1:6:fem.m_p.gridset.gsize)=new_surface;
0126
0127
0128 disp(sprintf('%s',' checking time stepping...'));
0129 [back,dt,time]=TimeStepping(md,solution,dt,time);
0130 if back,
0131 continue;
0132 end
0133
0134
0135 disp(sprintf('%s',' updating grid positions...'));
0136 grids_t=UpdateGridPosition(md,fem.m_t.grids,new_bed,new_thickness);
0137 grids_dhu=UpdateGridPosition(md,fem.m_dhu.grids,new_bed,new_thickness);
0138 grids_dbv=UpdateGridPosition(md,fem.m_dbv.grids,new_bed,new_thickness);
0139 grids_dv=UpdateGridPosition(md,fem.m_dv.grids,new_bed,new_thickness);
0140 grids_p=UpdateGridPosition(md,fem.m_p.grids,new_bed,new_thickness);
0141
0142
0143 time=time+dt;
0144 n=n+1;
0145
0146 end
0147
0148
0149 indx=1:6:fem.m_dv.gridset.gsize; indx=indx(fem.m_dv.tpart);
0150 indy=2:6:fem.m_dv.gridset.gsize; indy=indy(fem.m_dv.tpart);
0151 indz=3:6:fem.m_dv.gridset.gsize; indz=indz(fem.m_dv.tpart);
0152 for i=1:length(solution),
0153 solution2(i).vx=solution(i).u_g(indx)*md.yts;
0154 solution2(i).vy=solution(i).u_g(indy)*md.yts;
0155 solution2(i).vz=solution(i).u_g(indz)*md.yts;
0156 solution2(i).vel=sqrt(solution2(i).vx.^2+solution2(i).vy.^2+solution2(i).vz.^2);
0157 solution2(i).pressure=solution(i).pressure(indx)/10^5;
0158 solution2(i).temperature=solution(i).t_g(indx);
0159 solution2(i).melting=solution(i).melting_g(indx);
0160 solution2(i).thickness=solution(i).thickness(indx);
0161 solution2(i).surface=solution(i).surface(indx);
0162 solution2(i).bed=solution(i).bed(indx);
0163 solution2(i).time=solution(i).time/md.yts;
0164 end
0165
0166 md.transient_results=solution2;