0001 function md=cielodiagnostic(md)
0002
0003
0004
0005 t1=clock;
0006
0007
0008 infile='./cielo.dat';
0009
0010 analysis_dh='diagnostic_horiz';
0011 analysis_dbv='diagnostic_basevert';
0012 analysis_dv='diagnostic_vert';
0013
0014
0015 SetParameterDefaults;
0016
0017
0018
0019 m_dh=CreateFEMModel(md,analysis_dh);
0020 if strcmpi(md.type,'3d'),
0021 m_dbv=CreateFEMModel(md,analysis_dbv);
0022 m_dv=CreateFEMModel(md,analysis_dv);
0023 end
0024
0025
0026
0027 part=m_dh.part;
0028 tpart=m_dh.tpart';
0029 indx=1:6:m_dh.uset.gsize; indx=indx(tpart);
0030 indy=2:6:m_dh.uset.gsize; indy=indy(tpart);
0031 indz=3:6:m_dh.uset.gsize; indz=indz(tpart);
0032
0033
0034 md.dof=m_dh.uset.fsize;
0035
0036
0037 m_dh.rifts=md.rifts;
0038
0039
0040 disp(sprintf('\n%s',['computing horizontal velocities...']));
0041
0042 u_g=cielodiagnostic_core_nonlinear(m_dh,m_dh.params,{},analysis_dh);
0043
0044 if strcmpi(md.type,'3d'),
0045
0046
0047 u_g=SwitchPartitioning(u_g,'workspace',part,tpart,[1 2]);
0048
0049
0050 u_g=CieloVelocityExtrude(md,u_g);
0051
0052
0053 velocity_average=CieloHorizontalVelocityDepthAverage(md,u_g);
0054
0055
0056 u_g=SwitchPartitioning(u_g,'cluster',part,tpart,[1 2]);
0057 velocity_average=SwitchPartitioning(velocity_average,'cluster',part,tpart,[1 2]);
0058
0059
0060 disp(sprintf('\n%s',['computing basal vertical velocities...']));
0061 SetUset(m_dbv);
0062 u_g_basevert=cielodiagnostic_core_linear(m_dbv,m_dbv.params,struct('velocity',u_g,'velocity_average',velocity_average),analysis_dbv);
0063
0064
0065 SetUset(m_dv);
0066 m_dv.y_g=u_g_basevert;
0067
0068
0069 m_dv.y_s = Reducevectorg( m_dv.y_g);
0070
0071
0072 disp(sprintf('\n%s',['computing vertical velocities...']));
0073 u_g_vert=cielodiagnostic_core_linear(m_dv,m_dv.params,struct('velocity',u_g),analysis_dv);
0074
0075
0076 md.vx=u_g(indx)*md.yts;
0077 md.vy=u_g(indy)*md.yts;
0078 md.vz=u_g_vert(indz)*md.yts;
0079 md.vel=sqrt(md.vx.^2+md.vy.^2+md.vz.^2);
0080
0081
0082 md.pressure=md.rho_ice*md.g*(md.surface-md.z);
0083 else
0084
0085 md.vx=u_g(indx)*md.yts;
0086 md.vy=u_g(indy)*md.yts;
0087 md.vz=zeros(md.numberofgrids,1);
0088 md.vel=sqrt(md.vx.^2+md.vy.^2+md.vz.^2);
0089
0090
0091 md.pressure=md.rho_ice*md.g*(md.thickness);
0092 end
0093
0094
0095 t2=clock;
0096 disp(sprintf('\n%s\n',[' solution converged in ' num2str(etime(t2,t1)) ' seconds']));