icediagnostic3d

PURPOSE ^

ICEDIAGNOSTIC3D - diagnostic of a 3d model

SYNOPSIS ^

function u_g=icediagnostic3d(md,fem),

DESCRIPTION ^

ICEDIAGNOSTIC3D - diagnostic of a 3d model

   this routine computes the velocity field of a glacier in 3d

   Usage:
      u_g=icediagnostic3d(md,fem)

   See also: ICEDIAGNOSTIC, ICEDIAGNOSTIC2D, ICEDIAGNOSTIC_CORE_NONLINEAR ICEDIAGNOSTIC_CORE_LINEAR

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function u_g=icediagnostic3d(md,fem),
0002 %ICEDIAGNOSTIC3D - diagnostic of a 3d model
0003 %
0004 %   this routine computes the velocity field of a glacier in 3d
0005 %
0006 %   Usage:
0007 %      u_g=icediagnostic3d(md,fem)
0008 %
0009 %   See also: ICEDIAGNOSTIC, ICEDIAGNOSTIC2D, ICEDIAGNOSTIC_CORE_NONLINEAR ICEDIAGNOSTIC_CORE_LINEAR
0010 
0011 global gridset
0012 
0013 %recover inputs
0014 inputs=fem.inputs;
0015 
0016 %Compute surface slope if there are hutter elements
0017 if fem.ishutter,
0018     m_ss=fem.m_ss;
0019     m_dhu=fem.m_dhu;
0020 
0021     %Compute slope of the surface
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     %Project slopes computed at surface, onto the 3d volume
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     %Feed hutter solution with the slope input
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     %Compute Hutter solution
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     %Compute horizontal solution.
0047     disp(sprintf('\n%s',['computing horizontal velocities...']));
0048 
0049     %Run core solution
0050     u_g=icediagnostic_core_nonlinear(m_dh,'diagnostic_horiz',inputs);
0051 end
0052 
0053 %Extrude velocities for collapsed penta elements
0054 u_g=VelocityExtrude(md,u_g);
0055 
0056 %Compute depth averaged velocity and add it to inputs
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 %Compute basal vertical velocity, in 3d
0063 disp(sprintf('\n%s',['computing basal vertical velocity...']));
0064 u_g_basevert=icediagnostic_core_linear(fem.m_dbv,'diagnostic_basevert',inputs);
0065 
0066 %update dirichlet boundary conditions with these new base vertical velocities
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 %compute vertical velocities
0072 disp(sprintf('\n%s',['computing vertical velocity...']));
0073 u_g_vert=icediagnostic_core_linear(m_dv,'diagnostic_vert',inputs);
0074 
0075 %add contribution to vertical velocity
0076 u_g=u_g+u_g_vert;
0077 
0078 %Computation of pressure with Pattyn's assumptions (P=rho_ice*g*(s-z) in Pa)
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     %Compute slope of the bed
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     %update dirichlet boundary conditions with the velocities computed previously
0101     gridset=m_ds.gridset;
0102     m_ds.y_g=u_g; m_ds.ys=Reducevector_g(m_ds.y_g);
0103 
0104     %compute stokes velocities and pressure
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 %Reconditioning of the pressure
0110 u_g(4:6:m_dv.gridset.gsize)=u_g(4:6:m_dv.gridset.gsize)*md.stokesreconditioning;

Generated on Sun 29-Mar-2009 20:22:55 by m2html © 2003