0001 function u_g=icediagnostic3d(md,fem),
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011 global gridset
0012
0013
0014 inputs=fem.inputs;
0015
0016
0017 if fem.ishutter,
0018 m_ss=fem.m_ss;
0019 m_dhu=fem.m_dhu;
0020
0021
0022 if m_ss.params.debug
0023 disp(sprintf('\n%s',['computing surface slope (x and y derivatives)...']));
0024 end
0025 slopex=icediagnostic_core_linear(m_ss,'surface_slope_compute_x',inputs);
0026 slopey=icediagnostic_core_linear(m_ss,'surface_slope_compute_y',inputs);
0027
0028
0029 slopex=project3d(md,project2d(md,slopex(1:6:end),md.numlayers),'node');
0030 slopey=project3d(md,project2d(md,slopey(1:6:end),md.numlayers),'node');
0031
0032
0033 slopesurface=zeros(m_ss.gridset.gsize,1);
0034 slopesurface(1:6:m_ss.gridset.gsize,1)=slopex;
0035 slopesurface(2:6:m_ss.gridset.gsize,1)=slopey;
0036 inputs.slopesurface=slopesurface;
0037
0038
0039 disp(sprintf('\n%s',['computing hutter velocities...']));
0040 u_g=icediagnostic_core_linear(m_dhu,'diagnostic_hutter',inputs);
0041 end
0042
0043 if fem.ismacayealpattyn,
0044 m_dh=fem.m_dh;
0045
0046
0047 disp(sprintf('\n%s',['computing horizontal velocities...']));
0048
0049
0050 u_g=icediagnostic_core_nonlinear(m_dh,'diagnostic_horiz',inputs);
0051 end
0052
0053
0054 u_g=VelocityExtrude(md,u_g);
0055
0056
0057 gridset=fem.m_dv.gridset;
0058 velocity_average=HorizontalVelocityDepthAverage(md,u_g);
0059 inputs.velocity=u_g;
0060 inputs.velocity_average=velocity_average;
0061
0062
0063 disp(sprintf('\n%s',['computing basal vertical velocity...']));
0064 u_g_basevert=icediagnostic_core_linear(fem.m_dbv,'diagnostic_basevert',inputs);
0065
0066
0067 m_dv=fem.m_dv;
0068 gridset=m_dv.gridset;
0069 m_dv.y_g=u_g_basevert; m_dv.ys=Reducevector_g(m_dv.y_g);
0070
0071
0072 disp(sprintf('\n%s',['computing vertical velocity...']));
0073 u_g_vert=icediagnostic_core_linear(m_dv,'diagnostic_vert',inputs);
0074
0075
0076 u_g=u_g+u_g_vert;
0077
0078
0079 u_g(4:6:m_dv.gridset.gsize)=md.rho_ice*md.g*(md.surface-md.z)/md.stokesreconditioning;
0080
0081 if fem.isstokes,
0082 m_bs=fem.m_bs;
0083 m_ds=fem.m_ds;
0084
0085
0086 if m_bs.params.debug
0087 disp(sprintf('\n%s',['computing bed slope (x and y derivatives)...']));
0088 end
0089 slopex=icediagnostic_core_linear(m_bs,'bed_slope_compute_x',inputs);
0090 slopey=icediagnostic_core_linear(m_bs,'bed_slope_compute_y',inputs);
0091
0092 slopex=project3d(md,project2d(md,slopex(1:6:end),1),'node');
0093 slopey=project3d(md,project2d(md,slopey(1:6:end),1),'node');
0094
0095 slopebed(1:6:m_bs.gridset.gsize,1)=slopex;
0096 slopebed(2:6:m_bs.gridset.gsize,1)=slopey;
0097 inputs.slopebed=slopebed;
0098 inputs.velocity=u_g;
0099
0100
0101 gridset=m_ds.gridset;
0102 m_ds.y_g=u_g; m_ds.ys=Reducevector_g(m_ds.y_g);
0103
0104
0105 disp(sprintf('\n%s',['computing stokes velocities and pressure ...']));
0106 u_g=icediagnostic_core_nonlinear(m_ds,'diagnostic_stokes',inputs);
0107 end
0108
0109
0110 u_g(4:6:m_dv.gridset.gsize)=u_g(4:6:m_dv.gridset.gsize)*md.stokesreconditioning;