iceprognostic

PURPOSE ^

ICEPROGNOSTIC - prognostic solution

SYNOPSIS ^

function md=iceprognostic(md)

DESCRIPTION ^

ICEPROGNOSTIC - prognostic solution

   This routine is used to compute the evolution of the geometry
   of a model

   Usage:
      md=iceprognostic(md)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function md=iceprognostic(md)
0002 %ICEPROGNOSTIC - prognostic solution
0003 %
0004 %   This routine is used to compute the evolution of the geometry
0005 %   of a model
0006 %
0007 %   Usage:
0008 %      md=iceprognostic(md)
0009 
0010 %define global variables
0011 iceglobal
0012 
0013 %determine if run is parallel
0014 if strcmpi(md.cluster,'yes'), cluster=1; else cluster=0;end;
0015 
0016 %for now, only serial support is in
0017 if cluster,
0018     error('iceprognostic error message: parallel support not implemented yet');
0019 end
0020 
0021 %setup some parameters to be passed to the core solution
0022 params.sparsity=md.sparsity;
0023 params.solver_type=md.solver_type;
0024 params.eps_rel=md.eps_rel;
0025 params.eps_abs=md.eps_abs;
0026 params.debug=md.debug;
0027 
0028 %First, build elements,grids,loads, etc ... for prognostic
0029 m=CreateFemModel(md,'prognostic');
0030 
0031 %compute depth averaged horizontal velocity
0032 u_g=zeros(gridset.gsize,1);
0033 u_g(1:6:end)=md.vx/md.yts; %from m/yr to m/s
0034 u_g(2:6:end)=md.vy/md.yts; %from m/yr to m/s
0035 u_g(3:6:end)=md.vz/md.yts; %from m/yr to m/s
0036 
0037 if strcmpi(md.type,'3d'),
0038     velocity_average=HorizontalVelocityDepthAverage(md,u_g);
0039 else
0040     %the velocity is already depthaverdged
0041     velocity_average=u_g;
0042 end
0043 
0044 %move surface and bed velocities to first lower layer. This will be needed by the prognostic model
0045 melting=zeros(gridset.gsize,1);
0046 accumulation=zeros(gridset.gsize,1);
0047 if strcmpi(md.type,'3d'),
0048     melting(1:6:end)=ShiftLayers(md,md.melting,1,1);
0049     accumulation(1:6:end)=ShiftLayers(md,md.accumulation,md.numlayers,1);
0050 else
0051     melting(1:6:end)=md.melting/md.yts;           %from m/yr to m/s
0052     accumulation(1:6:end)=md.accumulation/md.yts; %from m/yr to m/s
0053 end
0054 
0055 %prepare inputs
0056 thickness=zeros(gridset.gsize,1);thickness(1:6:end)=md.thickness;
0057 %inputs=struct('thickness',thickness,'melting',melting,'accumulation',accumulation,...
0058 %'surface_vertical_velocity',ws,'basal_vertical_velocity',wb,'dt',md.dt,'velocity_average',velocity_average);
0059 inputs=struct('thickness',thickness,'melting',melting,'accumulation',accumulation,'dt',md.dt,'velocity_average',velocity_average);
0060 
0061 %Run core solution
0062 h_g=iceprognostic_core(m,'prognostic',inputs);
0063 
0064 %project collapsed thickness onto 3d mesh
0065 indx=1:6:m.gridset.gsize; indx=indx(m.tpart);
0066 if strcmpi(md.type,'3d'),
0067     md.new_thickness=project3d(md,project2d(md,h_g(indx),1),'node');
0068 else
0069     md.new_thickness=h_g(indx);
0070 end

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