0001 function md=cielotransient(md)
0002
0003
0004 analysis_t='thermaltransient';
0005 analysis_m='melting';
0006 analysis_dh='diagnostic_horiz';
0007 analysis_dbv='diagnostic_basevert';
0008 analysis_dv='diagnostic_vert';
0009 analysis_p='prognostic';
0010
0011
0012 SetParameterDefaults;
0013
0014
0015 md.vx=zeros(md.numberofgrids,1); md.vy=zeros(md.numberofgrids,1); md.vz=zeros(md.numberofgrids,1);
0016
0017
0018 disp(sprintf('\n reading data from model %s...',md.name));
0019
0020
0021 m_t=CreateFEMModel(md,analysis_t);
0022 m_m=CreateFEMModel(md,analysis_m);
0023 m_dh=CreateFEMModel(md,analysis_dh);
0024 m_dbv=CreateFEMModel(md,analysis_dbv);
0025 m_dv=CreateFEMModel(md,analysis_dv);
0026 m_p=CreateFEMModel(md,analysis_p);
0027
0028
0029 tpart=m_t.tpart;
0030 part=m_t.part;
0031 indx=1:6:m_t.uset.gsize; indx=indx(tpart);
0032 indy=2:6:m_t.uset.gsize; indy=indy(tpart);
0033 indz=3:6:m_t.uset.gsize; indz=indz(tpart);
0034 indp=4:6:m_t.uset.gsize; indp=indp(tpart);
0035
0036
0037 md.dof=m_dh.uset.fsize;
0038
0039 nsteps=md.ndt/md.dt;
0040
0041
0042 solution.u_g=zeros(uset.gsize,1);
0043 if ~isempty(md.vx), solution.u_g(1:6:uset.gsize)=md.vx(part)/md.yts; end
0044 if ~isempty(md.vy), solution.u_g(2:6:uset.gsize)=md.vy(part)/md.yts; end
0045 if ~isempty(md.vz), solution.u_g(3:6:uset.gsize)=md.vz(part)/md.yts; end
0046 solution.pressure=zeros(uset.gsize,1); solution.pressure(1:6:uset.gsize)=md.rho_ice*md.g*(md.surface(part)-md.z(part));
0047
0048 solution.t_g=zeros(uset.gsize,1);
0049 if ~isempty(md.temperature),
0050 solution.t_g(1:6:uset.gsize)=md.temperature(part);
0051 else
0052 solution.t_g(1:6:uset.gsize)=md.meltingpoint;
0053 end
0054 if ~isempty(md.melting), solution.melting_g(1:6:uset.gsize)=md.melting(part);end;
0055
0056
0057 solution.thickness=zeros(uset.gsize,1);solution.thickness(1:6:uset.gsize)=md.thickness(part);
0058 solution.surface=zeros(uset.gsize,1);solution.surface(1:6:uset.gsize)=md.surface(part);
0059 solution.bed=zeros(uset.gsize,1);solution.bed(1:6:uset.gsize)=md.bed(part);
0060 solution.time=0;
0061
0062
0063 thermalparams.beta=md.beta;
0064 thermalparams.meltingpoint=md.meltingpoint;
0065 thermalparams.latentheat=md.latentheat;
0066 thermalparams.heatcapacity=md.heatcapacity;
0067 thermalparams.penalty=md.penalty_thermal;
0068
0069 for n=1:nsteps,
0070
0071 disp(sprintf('\n%s%i/%i\n','time step: ',n,md.ndt/md.dt));
0072
0073 solution(n+1).time=(n)*md.dt;
0074
0075
0076 disp(sprintf('%s',' computing temperature...'));
0077
0078
0079 inputs=struct('velocity',solution(n).u_g,'pressure',solution(n).pressure,'temperature',solution(n).t_g,'dt',md.dt,...
0080 'thickness',solution(n).thickness,'bed',solution(n).bed,'surface',solution(n).surface);
0081
0082 [solution(n+1).t_g,m_t]=cielothermal_core(m_t,m_t.params,thermalparams,inputs,analysis_t);
0083
0084
0085 disp(sprintf('%s',' computing melting...'));
0086 inputs=struct('pressure',solution(n).pressure,'temperature',solution(n+1).t_g,'dt',md.dt,...
0087 'thickness',solution(n).thickness,'bed',solution(n).bed,'surface',solution(n).surface);
0088
0089 solution(n+1).melting_g=cielomelting_core(m_m,m_m.params,thermalparams,inputs,analysis_m);
0090
0091
0092 temperature_workspace=solution(n+1).t_g(indx);
0093 temperature_average2d_workspace=DepthAverage(md,temperature_workspace);
0094 temperature_average_workspace=project3d(md,temperature_average2d_workspace,'node');
0095 temperature_average=zeros(uset.gsize,1);temperature_average(1:6:uset.gsize)=temperature_average_workspace(part);
0096
0097 if md.debug,
0098 t_g_workspace=SwitchPartitioning(solution(n+1).t_g,'workspace',part,tpart,[1]);
0099 melting_g_workspace=SwitchPartitioning(solution(n+1).melting_g,'workspace',part,tpart,[1]);
0100 plot(md,'data',t_g_workspace(1:6:uset.gsize),'title','Temperature (K)','view',3,'data',melting_g_workspace(1:6:uset.gsize)*md.yts,'title','Melting (m/a)','view',3);pause(2);
0101 end
0102
0103
0104 disp(sprintf('%s',' computing horizontal velocities...'));
0105
0106 m_dh.rifts=md.rifts;
0107 inputs=struct('thickness',solution(n).thickness,'bed',solution(n).bed,'surface',solution(n).surface,...
0108 'temperature',solution(n+1).t_g,'temperature_average',temperature_average);
0109 solution(n+1).u_g=cielodiagnostic_core_nonlinear(m_dh,m_dh.params,inputs,analysis_dh);
0110
0111
0112 u_g_workspace=SwitchPartitioning(solution(n+1).u_g,'workspace',part,tpart,[1 2]); u_g_workspace=CieloVelocityExtrude(md,u_g_workspace);
0113
0114
0115 solution(n+1).u_g=SwitchPartitioning(u_g_workspace,'cluster',part,tpart,[1 2]);
0116
0117 if md.debug,
0118 plot(md,'data',sqrt(u_g_workspace(1:6:uset.gsize).^2+u_g_workspace(2:6:uset.gsize))*md.yts,'view',3,'title','Horizontal velocity (m/a)');pause(2);
0119 end
0120
0121
0122 disp(sprintf('%s',' computing vertical base velocity...'));
0123 inputs=struct('velocity',solution(n+1).u_g,'thickness',solution(n).thickness,'bed',solution(n).bed,'surface',solution(n).surface);
0124 u_g_basevert=cielodiagnostic_core_linear(m_dbv,m_dbv.params,inputs,analysis_dbv);
0125
0126 if md.debug,
0127 u_g_basevert_workspace=SwitchPartitioning(u_g_basevert,'workspace',part,tpart,[3]);
0128 plot(md,'data',u_g_basevert_workspace(3:6:uset.gsize)*md.yts,'layer',1,'title','Basal Vertical velocity (m/a)');pause(2);
0129 end
0130
0131
0132 disp(sprintf('%s',' computing vertical velocity...'));
0133 m_dv.y_g=u_g_basevert;
0134
0135
0136 m_dv.y_s = Reducevectorg( m_dv.y_g,m_dv.uset);
0137
0138
0139 velocity_average_workspace=CieloHorizontalVelocityDepthAverage(md,u_g_workspace); velocity_average=SwitchPartitioning(velocity_average_workspace,'cluster',part,tpart,[1 2]);
0140
0141
0142 inputs=struct('velocity',solution(n+1).u_g,'velocity_average',velocity_average,'thickness',solution(n).thickness,'bed',solution(n).bed,'surface',solution(n).surface);
0143 u_g_vert=cielodiagnostic_core_linear(m_dv,m_dv.params,inputs,analysis_dv);
0144
0145
0146
0147 solution(n+1).u_g(3:6:uset.gsize)=u_g_vert(3:6:uset.gsize);
0148
0149 if md.debug,
0150 u_g_workspace=SwitchPartitioning(solution(n+1).u_g,'workspace',part,tpart,[3]);
0151 plot(md,'data',u_g_workspace(3:6:uset.gsize)*md.yts,'view',3,'title','Vertical velocity (m/a)');pause(2);
0152 end
0153
0154
0155 disp(sprintf('%s',' computing pressure...'));
0156 solution(n+1).pressure=zeros(uset.gsize,1); solution(n+1).pressure(1:6:uset.gsize)=md.rho_ice*md.g*(md.surface(part)-md.z(part));
0157
0158
0159
0160 u_g_vert_workspace=SwitchPartitioning(u_g_vert,'workspace',part,tpart,[3]);
0161 ws_workspace=ShiftLayers(md,u_g_vert_workspace(3:6:uset.gsize),md.numlayers,1);
0162 wb_workspace=ShiftLayers(md,u_g_vert_workspace(3:6:uset.gsize),1,1);
0163
0164 melting_g_workspace=SwitchPartitioning(solution(n+1).melting_g,'workspace',part,tpart,[1]);
0165 melting_workspace=ShiftLayers(md,melting_g_workspace(1:6:uset.gsize),1,1);
0166
0167 accumulation_workspace=ShiftLayers(md,md.accumulation,md.numlayers,1);
0168
0169
0170
0171 ws_g=zeros(uset.gsize,1);ws_g(1:6:uset.gsize)=ws_workspace(part);
0172 wb_g=zeros(uset.gsize,1);wb_g(1:6:uset.gsize)=wb_workspace(part);
0173 melting_g=zeros(uset.gsize,1);melting_g(1:6:uset.gsize)=melting_workspace(part);
0174 accumulation_g=zeros(uset.gsize,1);accumulation_g(1:6:uset.gsize)=accumulation_workspace(part)/md.yts;
0175
0176
0177 if md.debug,
0178 plot(md,'data',ws_workspace*md.yts,'title','Surface velocity projected on layer 1','view',3,'data',wb_workspace*md.yts,'title','Basal velocity projected on layer 1','view',3,...
0179 'data',melting_workspace*md.yts,'title','Bed melting projected on layer 1 (m/yr)','view',3,'data',accumulation_workspace,'title','Surface accumulation projected on layer 1 (m/yr)','view',3);pause(2);
0180 end
0181
0182
0183 disp(sprintf('%s',' computing new thickness...'));
0184 inputs=struct('thickness',solution(n).thickness,'melting',melting_g,'accumulation',accumulation_g,...
0185 'surface_vertical_velocity',ws_g,'basal_vertical_velocity',wb_g,'dt',md.dt,'velocity_average',velocity_average);
0186 new_thickness=cieloprognostic_core(m_p,m_p.params,inputs,analysis_p);
0187
0188
0189 new_thickness_workspace=SwitchPartitioning(new_thickness,'workspace',part,tpart,[1]);
0190 new_thickness_workspace=project3d(md,project2d(md,new_thickness_workspace(1:6:end),1),'node');
0191
0192
0193
0194 [new_bed_workspace,new_surface_workspace,]=UpdateGeometry(md,new_thickness_workspace,solution(n).bed(indx),solution(n).surface(indx));
0195
0196 if md.debug,
0197 plot(md,'data',new_thickness_workspace,'view',3,'title','New thickness (m)','data',new_surface_workspace,'view',3,'title','New surface (m)','data',new_bed_workspace,'view',3,'title','New bed (m)');pause(2);
0198 end
0199
0200
0201 solution(n+1).thickness=zeros(uset.gsize,1);solution(n+1).thickness(1:6:uset.gsize)=new_thickness_workspace(part);
0202 solution(n+1).bed=zeros(uset.gsize,1);solution(n+1).bed(1:6:uset.gsize)=new_bed_workspace(part);
0203 solution(n+1).surface=zeros(uset.gsize,1);solution(n+1).surface(1:6:uset.gsize)=new_surface_workspace(part);
0204
0205
0206 [m_t.bgpdt,m_t.bgpdtb]=CieloUpdateGridPosition(m_t.bgpdt,m_t.bgpdtb,solution(n+1).bed,solution(n+1).thickness);
0207
0208 [m_dh.bgpdt,m_dh.bgpdtb]=CieloUpdateGridPosition(m_dh.bgpdt,m_dh.bgpdtb,solution(n+1).bed,solution(n+1).thickness);
0209
0210 [m_dbv.bgpdt,m_dbv.bgpdtb]=CieloUpdateGridPosition(m_dbv.bgpdt,m_dbv.bgpdtb,solution(n+1).bed,solution(n+1).thickness);
0211
0212 [m_dv.bgpdt,m_dv.bgpdtb]=CieloUpdateGridPosition(m_dv.bgpdt,m_dv.bgpdtb,solution(n+1).bed,solution(n+1).thickness);
0213
0214 [m_p.bgpdt,m_p.bgpdtb]=CieloUpdateGridPosition(m_p.bgpdt,m_p.bgpdtb,solution(n+1).bed,solution(n+1).thickness);
0215 end
0216
0217
0218
0219
0220 for i=1:length(solution),
0221 solution2(i).vx=solution(i).u_g(indx)*md.yts;
0222 solution2(i).vy=solution(i).u_g(indy)*md.yts;
0223 solution2(i).vz=solution(i).u_g(indz)*md.yts;
0224 solution2(i).vel=sqrt(solution2(i).vx.^2+solution2(i).vy.^2+solution2(i).vz.^2);
0225 solution2(i).pressure=solution(i).pressure(indx)/10^5;
0226 solution2(i).temperature=solution(i).t_g(indx);
0227 solution2(i).melting=solution(i).melting_g(indx);
0228 solution2(i).thickness=solution(i).thickness(indx);
0229 solution2(i).surface=solution(i).surface(indx);
0230 solution2(i).bed=solution(i).bed(indx);
0231 solution2(i).time=solution(i).time/md.yts;
0232 end
0233 md.transient_results=solution2;