http://websrv.cs.umt.edu/isis/api.php?action=feedcontributions&user=Meierbtw&feedformat=atomInteractive System for Ice sheet Simulation - User contributions [en]2020-02-26T07:52:10ZUser contributionsMediaWiki 1.21.1http://websrv.cs.umt.edu/isis/index.php/Adding_dH/dt_moduleAdding dH/dt module2009-08-13T22:57:00Z<p>Meierbtw: /* Step 6 */</p>
<hr />
<div>==Overview==<br />
<br />
This page contains step-by-step instructions for adding a new module to Glimmer-CISM. In this case, the module is a first-order, upwinding advection scheme for mass transport (dH/dt) using velocities calculated from a higher-order dynamics model. The procedure, however, is generic and could apply to adding almost any module. The goal here is not to get lost in the details of the upwinding scheme (we'll provide the code for that later) but rather to understand the steps for adding a module in an incremental, structured way that allows you to track small errors and fix them at each step. <br />
<br />
*To go directly to some exercises using the higher-order dynamics and the evolution scheme, go to [[Ice Sheet Evolution Experiments]].<br />
*For mathematical background on upwinded transport, consult the [[Solving the equation for thickness evolution]] pages.<br />
<br />
<br />
== '''Step 0''' ==<br />
<br />
We first need to download the code, configure it, and check that we can get a successful initial build (presumably we've done that already if you've been able to do the higher-order test suite exercises). If not, then type the following from within the highest level directory containing the source code:<br />
<br />
./bootstrap<br />
./configure --with-netcdf=/path/to_netcdf !!! Note: specify your own path to netCDF libs here !!!<br />
make<br />
<br />
== '''Step 1''' ==<br />
<br />
In step 1, we will create a module containing "empty" subroutines and add this module to the build.<br />
<br />
Create a module '''fo_upwind_advect.F90''' that contains the necessary subroutines. For now, these will just be "stubs", which we will fill in later.<br />
<br />
<source lang=fortran>module fo_upwind_advect<br />
<br />
! subroutines for mass advection scheme based on 1st order upwinding<br />
<br />
contains<br />
<br />
!----------------------------------------------------------------------<br />
<br />
subroutine fo_upwind_advect_init( )<br />
<br />
! initialization for 1st-order upwinding mass advection<br />
<br />
end subroutine fo_upwind_advect_init<br />
<br />
!----------------------------------------------------------------------<br />
<br />
subroutine fo_upwind_advect_final( )<br />
<br />
! finalization for 1st-order upwinding mass advection<br />
<br />
end subroutine fo_upwind_advect_final<br />
<br />
!----------------------------------------------------------------------<br />
<br />
subroutine fo_upwind_advect_driver( )<br />
<br />
! driver for 1st-order upwind mass advection<br />
<br />
end subroutine fo_upwind_advect_driver<br />
<br />
!----------------------------------------------------------------------<br />
<br />
subroutine fo_upwind_advect_main( )<br />
<br />
! 1st-order upwinding mass advection<br />
<br />
end subroutine fo_upwind_advect_main<br />
<br />
!----------------------------------------------------------------------<br />
<br />
end module fo_upwind_advect</source><br />
<br />
<br />
Add this new module to the build by editing '''Makefile.am''' in '''src/fortran/''' ...<br />
<br />
<source lang=fortran>libglide_a_SOURCES = glide.F90 glide_setup.F90 glide_types.F90 glide_temp.F90 \<br />
glide_bwater.F90 glide_deriv.F90 xls.F90 ice3d_lib.F90 \<br />
glide_velo_higher.F90 glide_thck.F90 glide_velo.F90 \<br />
glide_mask.F90 glide_stop.F90 glide_io.F90 \<br />
glide_nc_custom.F90 isostasy.F90 isostasy_el.F90\<br />
isostasy_setup.F90 isostasy_types.F90 glide_lithot.F90\<br />
glide_lithot3d.F90 glide_lithot1d.F90 glide_profile.F90\<br />
glide_diagnostics.F90 glissade.F90 glissade_remap.F90\<br />
glissade_velo.F90 glissade_constants.F90 glide_vertint.F90\<br />
glide_thckmask.F90 glide_nonlin.F90 glide_grids.F90\<br />
glam.F90 glam_strs2.F90 glam_thck_ppm.F90\<br />
remap_advection.F90 remap_glamutils.F90 glide_ground.F90</source><br />
<br />
becomes ...<br />
<br />
<source lang=fortran>libglide_a_SOURCES = glide.F90 glide_setup.F90 glide_types.F90 glide_temp.F90 \<br />
glide_bwater.F90 glide_deriv.F90 xls.F90 ice3d_lib.F90 \<br />
glide_velo_higher.F90 glide_thck.F90 glide_velo.F90 \<br />
glide_mask.F90 glide_stop.F90 glide_io.F90 \<br />
glide_nc_custom.F90 isostasy.F90 isostasy_el.F90\<br />
isostasy_setup.F90 isostasy_types.F90 glide_lithot.F90\<br />
glide_lithot3d.F90 glide_lithot1d.F90 glide_profile.F90\<br />
glide_diagnostics.F90 glissade.F90 glissade_remap.F90\<br />
glissade_velo.F90 glissade_constants.F90 glide_vertint.F90\<br />
glide_thckmask.F90 glide_nonlin.F90 glide_grids.F90\<br />
glam.F90 glam_strs2.F90 glam_thck_ppm.F90\<br />
remap_advection.F90 remap_glamutils.F90 glide_ground.F90 \ ! CHANGED LAST TWO LINES<br />
fo_upwind_advect.F90</source><br />
<br />
[Note: the comment in the 2nd to last line may cause a problem for make.]<br />
<br />
Now we stop adding/editing things for a minute and type '''make''' to rebuild and make sure we haven't broken the code yet. If you get a successful build, then proceed to '''Step 2'''.<br />
<br />
== '''Step 2''' ==<br />
<br />
In step 2, we will do some minor updates to other parts of the code so that there is an option to call our 1st order upwinding scheme when we want to evolve the ice thickness (that is, rather than using some other scheme).<br />
<br />
In '''glide_types.F90''', add the following to the transport scheme options around line 100 ...<br />
<br />
<source lang=fortran> integer, parameter :: EVOL_PSEUDO_DIFF = 0<br />
integer, parameter :: EVOL_ADI = 1<br />
integer, parameter :: EVOL_DIFFUSION = 2<br />
integer, parameter :: EVOL_INC_REMAP = 3<br />
integer, parameter :: EVOL_FO_UPWIND = 4 ! ADDED THIS LINE</source><br />
<br />
... similarly, update some lines around 207 for documentation ...<br />
<br />
<source lang=fortran> !*FD Thickness evolution method:<br />
!*FD \begin{description}<br />
!*FD \item[0] Pseudo-diffusion approach <br />
!*FD \item[2] Diffusion approach (also calculates velocities)<br />
!*FD \item[3] Incremental remapping<br />
!*FD \item[4] 1st-order upwinding scheme ! ADDED THIS LINE<br />
!*FD \end{description}</source><br />
<br />
In '''glide_setup.F90''', around line 475, change ...<br />
<br />
<source lang=fortran> character(len=*), dimension(0:4), parameter :: evolution = (/ & ! CHANGED 0:3 TO 0:4<br />
'pseudo-diffusion ', &<br />
'ADI scheme ', &<br />
'iterated diffusion ', &<br />
'remap thickness ', & <br />
'1st order upwind ' /) ! ADDED THIS LINE</source><br />
<br />
Time to type '''make''' again and check for a successful build. If so, proceed to '''Step 3'''.<br />
<br />
<br />
== '''Step 3''' ==<br />
<br />
Now we will add use statements, if constructs, and dummy calls to the advection subroutines from the appropriate places in the code. Hence, we are setting up all the necessary structure to call the subroutines but are not passing any arguments to them yet (nor are they actually doing anything yet).<br />
<br />
In '''glide.F90''', add the following use statement to subroutine '''glide_initialise''' ...<br />
<br />
<source lang=fortran> use glam_strs2, only : glam_velo_fordsiapstr_init<br />
use remap_glamutils, only : horizontal_remap_init<br />
use fo_upwind_advect, only : fo_upwind_advect_init ! ADDED<br />
</source><br />
<br />
Now add if construct and call to the initialization stub in '''fo_upwind_advect_mod.F90''', around line 217 of '''glide.F90''', just after the if construct for EVOL_INC_REMAP ...<br />
<br />
<source lang=fortran> if (model%options%whichevol== EVOL_FO_UPWIND ) then<br />
call fo_upwind_advect_init( )<br />
end if</source><br />
<br />
In '''glide_stop.F90''', add the following to the use statements ...<br />
<br />
<source lang=fortran>module glide_stop<br />
<br />
use glide_types<br />
use glimmer_log<br />
use remap_glamutils<br />
use fo_upwind_advect, only : fo_upwind_advect_final ! ADDED</source><br />
<br />
Add the necessary if construct to '''glide_stop.F90''' around line 150, as we did with the initialization routine ...<br />
<br />
<source lang=fortran> if (model%options%whichevol== EVOL_FO_UPWIND ) then<br />
call fo_upwind_advect_final( )<br />
endif</source><br />
<br />
Finally, update the use statements in '''glide.F90''' around line 340, and the case construct for ice sheet evolution, around line 380, so that the 1st order upwinding subroutines can be called ...<br />
<br />
<source lang=fortran> use glide_thck<br />
use glide_velo<br />
use glide_ground<br />
use glide_setup<br />
use glide_temp<br />
use glide_mask<br />
use isostasy<br />
use glam, only: inc_remap_driver<br />
use fo_upwind_advect, only: fo_upwind_advect_driver ! ADDED </source><br />
<br />
... and add a call to the driver routine to the case construct (passing no args yet) ...<br />
<br />
<source lang=fortran> case(EVOL_FO_UPWIND) ! Use first order upwind scheme for mass transport<br />
<br />
call fo_upwind_advect_driver( )</source><br />
<br />
Make sure to place it BEFORE the existing "end select" statement!<br />
<br />
Type '''make''' again and check for a successful build. If so, proceed to '''Step 4'''.<br />
<br />
== '''Step 4''' ==<br />
<br />
At this point we've done nothing but build the necessary structures in the code. Now we will start filling in the subroutines with the necessary variable definitions. We'll follow this by passing arguments and, if everything still works, we can actually start to do something w/ those arguments within the subroutines themselves.<br />
<br />
Note that, in practice, one doesn't always know ahead of time exactly what variables are needed within a particular subroutine, and there is some trial and error while figuring it out. Here, we'll assume we've thought this out really well ahead of time and we know exactly what we need to do before we ever started writing the code (like punch-card programmers of old), in which case we know exactly which arguments are needed when we start.<br />
<br />
<br />
First, we'll fill out the rest of the "driver" subroutine, which has two parts, (1) calling the higher-order dynamics subroutines '''run_ho_diagnostic''' to get the appropriate velocity fields and, (2) calling '''fo_upwind_advect_main''', which uses those velocity fields, the thickness field, and then moves mass around to calculate a new thickness field. In subroutine '''fo_upwind_advect_driver''' we add the following ...<br />
<br />
<source lang=fortran> subroutine fo_upwind_advect_driver( model ) ! ADDED ARG<br />
<br />
! driver routine for the 1st order, upwind mass transport scheme<br />
<br />
type(glide_global_type), intent(inout) :: model ! ADDED <br />
<br />
call run_ho_diagnostic(model) ! ADDED CALL AND ARG<br />
<br />
call fo_upwind_advect_main( ) ! ADDED CALL <br />
<br />
... </source><br />
<br />
Note that we also passed the model in to the driver and gave "model" a type definition. We also need to add the appropriate use statements for "model" and to access the higher-order dynamics routines, which we'll do below. <br />
<br />
Next, add argument definitions to the main subroutine in '''fo_upwind_advect.F90''' ...<br />
<br />
<source lang=fortran> subroutine fo_upwind_advect_main( thck, stagthck, acab, dt, uflx, vflx, ewn, nsn, dew, dns ) ! CHANGED<br />
<br />
! 1st-order upwinding mass advection<br />
<br />
implicit none ! ADDED<br />
<br />
real (kind = dp), intent(in) :: dt ! time step (sec)<br />
real (kind = dp), dimension(:,:), intent(inout) :: thck ! thickness on normal grid (m)<br />
real (kind = dp), dimension(:,:), intent(in) :: stagthck ! thickness on staggered grid (m)<br />
real (kind = sp), dimension(:,:), intent(in) :: acab ! surf mass balance (accum/ablation) (m/sec)<br />
real (kind = dp), dimension(:,:), intent(in) :: uflx, vflx ! flux in x,y directions (m^2/sec)<br />
real (kind = dp), intent(in) :: dew, dns ! grid spacing in x,y directions (m)<br />
integer, intent(in) :: ewn, nsn ! no. of grid points in x,y directions<br />
<br />
real (kind = dp) :: He, Hw, Hn, Hs, ue, uw, vn, vs ! ADDED<br />
integer :: ew, ns<br />
<br />
end subroutine fo_upwind_advect_main</source><br />
<br />
<br />
Because these definitions require kind "dp", we need to add a use statment at the start of the module. We've also added the necessary use statement for accessing "model" and the higher-order dynamics subroutines (from above) and a few more we'll use later on ...<br />
<br />
<source lang=fortran>module fo_upwind_advect<br />
<br />
! subroutines for mass advection scheme based on 1st order upwinding<br />
<br />
! ADDED<br />
use glimmer_paramets, only: sp, dp, len0, thk0, tim0, vel0, tim0, acc0, scyr ! scales for swapping between dim and non-dim vars<br />
use glide_types<br />
use glide_velo_higher<br />
<br />
private<br />
public :: fo_upwind_advect_init, fo_upwind_advect_driver, fo_upwind_advect_final<br />
<br />
... etc ... </source><br />
<br />
To avoid someone mucking around with variables in the module that we don't want them to touch, we've also added the "private" and "public" statements; the only parts of the module that can be accessed from outside are through calls to the three "public" subroutines. <br />
<br />
<br />
The arguments that we defined above in '''fo_upwind_advect_main''' need to be passed in from the driver routine '''fo_upwind_advect_driver''' to the subroutine that does all of the work, '''fo_upwind_advect_main'''. To do this, we must access them from the derived type "model". The entire driver subroutine then becomes ...<br />
<br />
<source lang=fortran> subroutine fo_upwind_advect_driver( model )<br />
<br />
! driver routine for the 1st order, upwind mass transport scheme<br />
<br />
type(glide_global_type), intent(inout) :: model<br />
<br />
call run_ho_diagnostic(model) ! get velocities and fluxes from HO dynamic subroutines<br />
<br />
call fo_upwind_advect_main( model%geometry%thck, model%geomderv%stagthck, &<br />
model%climate%acab, model%numerics%dt, &<br />
model%velocity_hom%uflx,model%velocity_hom%vflx, &<br />
model%general%ewn, model%general%nsn, &<br />
model%numerics%dew, model%numerics%dns )<br />
<br />
end subroutine fo_upwind_advect_driver</source><br />
<br />
<br />
Finally, we pass the derived type "model" during the call to the driver subroutine in '''glide.F90''' ...<br />
<br />
<source lang=fortran> case(EVOL_FO_UPWIND) ! Use first order upwind scheme for mass transport<br />
<br />
call fo_upwind_advect_driver( model ) ! ADDED ARGUMENT<br />
<br />
end select</source><br />
<br />
Type '''make''' again and check for a successful build. If so, proceed to '''Step 5'''.<br />
<br />
== '''Step 5''' ==<br />
<br />
Now let's fill in the init and finalization subroutines in '''fo_upwind_advect.F90'''. This takes some thought ahead of time and might require a bit of trial and error. Again, for now we'll assume that we've thought this out really well ahead of time and we know exactly what we need.<br />
<br />
First, declare any other necessary variables at the start of the module '''fo_upwind_advection.F90'''. Here, these are the allocatable work arrays that we'll use in '''fo_upwind_advect_main'''. The beginning of the module becomes ...<br />
<br />
<source lang=fortran>module fo_upwind_advect<br />
<br />
!----------------------------------------------------------------------<br />
<br />
! init, finalize, and driver subroutines for mass advection based on 1st order upwinding<br />
<br />
use glimmer_paramets, only: sp, dp, len0, thk0, tim0, vel0, tim0, acc0, scyr<br />
use glide_types<br />
use glide_velo_higher<br />
<br />
private<br />
public :: fo_upwind_advect_init, fo_upwind_advect_driver, fo_upwind_advect_final<br />
<br />
! allocatable work arrays ! ADDED<br />
real (kind = dp), allocatable, dimension(:,:) :: &<br />
ubar, vbar, &<br />
ubar_grid, vbar_grid, &<br />
flux_net, thck_grid, &<br />
mask, thck_old<br />
<br />
contains<br />
<br />
... etc ...</source><br />
<br />
Now change the initialization and finalization subroutines ...<br />
<br />
<source lang=fortran> subroutine fo_upwind_advect_init( ewn, nsn ) ! ADDED ARGS HERE<br />
<br />
! initialization for 1st-order upwinding mass advection<br />
<br />
implicit none ! ADDED<br />
<br />
integer, intent(in) :: ewn, nsn ! horizontal grid dimensions ! ADDED TYPE DEF<br />
<br />
integer :: errstat ! ADDED FOR ERROR HANDLING<br />
<br />
! allocate work arrays ! ADDED THESE<br />
allocate( ubar(ewn-1,nsn-1), stat=errstat ); ubar = 0.0_dp<br />
allocate( vbar(ewn-1,nsn-1), stat=errstat ); vbar = 0.0_dp<br />
allocate( ubar_grid(ewn+1,nsn+1), stat=errstat ); ubar_grid = 0.0_dp<br />
allocate( vbar_grid(ewn+1,nsn+1), stat=errstat ); vbar_grid = 0.0_dp<br />
allocate( thck_grid(ewn+2,nsn+2), stat=errstat ); thck_grid = 0.0_dp<br />
allocate( flux_net(ewn,nsn), stat=errstat ); flux_net = 0.0_dp<br />
allocate( mask(ewn,nsn), stat=errstat ); mask = 0.0_dp<br />
allocate( thck_old(ewn,nsn), stat=errstat ); thck_old = 0.0_dp<br />
<br />
if ( errstat /= 0 ) then ! ADDED FOR ERROR HANDLING<br />
print *, 'error: allocation in fo_upwind_advect failed!'<br />
stop<br />
end if<br />
<br />
end subroutine fo_upwind_advect_init</source><br />
<br />
Note that ewn, nsn are passed in above, so we need to make sure they are passed from the main code <br />
where this call sits. In '''glide.F90''' we have ...<br />
<br />
<source lang=fortran> call fo_upwind_advect_init( model%general%ewn, model%general%nsn ) ! ADDED ARGS</source><br />
<br />
As with the initialization subroutine, we add deallocation statements for work arrays in the finalization subroutine ...<br />
<br />
<source lang=fortran> subroutine fo_upwind_advect_final( )<br />
<br />
! finalization for 1st-order upwinding mass advection<br />
<br />
implicit none ! ADDED<br />
<br />
! deallocate work arrays ! ADDED THESE<br />
if( allocated( ubar ) ) deallocate( ubar )<br />
if( allocated( vbar ) ) deallocate( vbar )<br />
if( allocated( ubar_grid ) ) deallocate( ubar_grid )<br />
if( allocated( vbar_grid ) ) deallocate( vbar_grid )<br />
if( allocated( thck_grid ) ) deallocate( thck_grid )<br />
if( allocated( flux_net ) ) deallocate( flux_net )<br />
if( allocated( mask ) ) deallocate( mask )<br />
if( allocated( thck_old ) ) deallocate( thck_old )<br />
<br />
<br />
end subroutine fo_upwind_advect_final</source><br />
<br />
Type '''make''' again and check for a successful build. If so, proceed to '''Step 6'''.<br />
<br />
== '''Step 6''' ==<br />
<br />
Now we are actually ready to do something ... that is, fill in the guts of the new subroutine, where the thickness evolution calculation takes place. Rather than have you figure out how to code up the 1st-order upwinding scheme on your own, we'll provide you with the chunk of code to do that (below), noting that the details of the calculation scheme are discussed [[Solving the equation for thickness evolution|HERE]]. <br />
<br />
<source lang=fortran>! ----------------------------------<br />
<br />
subroutine fo_upwind_advect_main( thck, stagthck, acab, dt, uflx, vflx, ewn, nsn, dew, dns )<br />
<br />
! 1st-order upwinding mass advection that uses a finite-volume like scheme for <br />
! mass conservation. Velocities from the staggered grid (B-grid) are averaged onto the <br />
! faces of the non-staggered grid (i.e. faces of the grid where scalers like thickness live). <br />
! Thus, the averaged velocities exist on a C-grid, allowing mass transport to be treated <br />
! in a finite-volume manner; depth averaged velocities give the fluxes out of each cell <br />
! centered on a thickness point and the thickness advected is chosen according to upwinding.<br />
! <br />
! Note that this works at the calving front because a non-zero staggered thickness there <br />
! defines the velocities there. These velocites can be used to define the velocity at<br />
! the face of the last non-zero thickness cell (on the normal grid) which corresponds to<br />
! the location of the calving front. <br />
<br />
implicit none<br />
<br />
real (kind = dp), intent(in) :: dt<br />
real (kind = dp), dimension(:,:), intent(inout) :: thck<br />
real (kind = dp), dimension(:,:), intent(in) :: stagthck<br />
real (kind = sp), dimension(:,:), intent(in) :: acab<br />
real (kind = dp), dimension(:,:), intent(in) :: uflx, vflx<br />
real (kind = dp), intent(in) :: dew, dns<br />
integer, intent(in) :: ewn, nsn<br />
<br />
real (kind = dp) :: He, Hw, Hn, Hs, ue, uw, vn, vs ! upwinding variables and interface velocities<br />
<br />
integer :: ew, ns<br />
<br />
where( stagthck > 0.0_dp ) ! calculate the depth-ave velocities<br />
ubar = uflx / stagthck<br />
vbar = vflx / stagthck<br />
end where<br />
<br />
where( thck > 0.0_dp ) ! mask for eventually removing flux outside of the original domain<br />
mask = 1.0_dp ! (i.e. stuff that moves past the calving front goes away)<br />
else where<br />
mask = 0.0_dp<br />
end where<br />
<br />
thck_old = thck ! save the old thickness for debugging purposes<br />
<br />
! fill in the interior values on the extended velocity grid (extended B-grid)<br />
ubar_grid(2:ewn,2:nsn) = ubar<br />
vbar_grid(2:ewn,2:nsn) = vbar<br />
<br />
! fill in the interior values on the extended thickness grid<br />
thck_grid(2:ewn+1,2:nsn+1) = thck(:,:)<br />
<br />
! calculate the interface velocities from the extended B-grid, then use upwinding<br />
! criterion to advect thickness in or out of cells (NOTE that parts of this could<br />
! probably be vectorized at some point)<br />
do ns = 1, nsn<br />
do ew = 1, ewn<br />
<br />
! interface depth-ave velocities<br />
ue = ( ubar_grid(ew+1,ns+1) + ubar_grid(ew+1,ns) ) / 2.0d0<br />
uw = ( ubar_grid(ew,ns+1) + ubar_grid(ew,ns) ) / 2.0d0<br />
vn = ( vbar_grid(ew,ns+1) + vbar_grid(ew+1,ns+1) ) / 2.0d0<br />
vs = ( vbar_grid(ew,ns) + vbar_grid(ew+1,ns) ) / 2.0d0<br />
<br />
! choose thickness to advect based on upwinding<br />
if( ue > 0.0d0 )then<br />
He = - thck_grid(ew+1,ns+1) ! negative signs necessary so that flux to the east<br />
else ! results in mass loss in this volume (and vice versa)<br />
He = - thck_grid(ew+2,ns+1)<br />
end if<br />
if( uw > 0.0d0 )then<br />
Hw = thck_grid(ew,ns+1)<br />
else<br />
Hw = thck_grid(ew+1,ns+1)<br />
end if<br />
if( vn > 0.0d0 )then<br />
Hn = - thck_grid(ew+1,ns+1) ! negative signs here as above for ue, and He<br />
else<br />
Hn = - thck_grid(ew+1,ns+2)<br />
end if<br />
if( vs > 0.0d0 )then<br />
Hs = thck_grid(ew+1,ns)<br />
else<br />
Hs = thck_grid(ew+1,ns+1)<br />
end if<br />
<br />
! net flux into/out of each cell<br />
flux_net(ew,ns) = ( ue*He*dns + uw*Hw*dns + vn*Hn*dew + vs*Hs*dew )<br />
<br />
end do<br />
end do<br />
<br />
thck = thck_old + ( 1 / (dns * dew) * flux_net ) * dt + (acab * dt)<br />
<br />
! debugging<br />
print *, ' '<br />
print *, 'net volume change = ', sum( (thck-thck_old)*mask )*thk0 *dew*dns*len0**2<br />
print *, 'net calving flux = ', sum( thck * (1.0d0-mask) )*thk0*dew*dns*len0**2<br />
print *, '(for the confined shelf experiment, the above two should sum to ~0)'<br />
print *, 'mean accum/ablat rate = ', sum( acab * mask ) / sum(mask) / (dt*tim0) * scyr<br />
print *, 'mean dH/dt = ', sum( (thck-thck_old)*mask )*thk0 / sum(mask) / (dt*tim0) * scyr<br />
print *, 'sum of flux change (should be ~0) = ', sum( flux_net*vel0*thk0*len0 )<br />
print *, ' '<br />
! pause<br />
<br />
thck = thck * mask ! remove any mass advected outside of initial domain<br />
<br />
where( thck < 0.0_dp ) ! guard against thickness going negative<br />
thck = 0.0_dp<br />
end where<br />
<br />
end subroutine fo_upwind_advect_main<br />
<br />
! ----------------------------------</source><br />
<br />
==Ice Sheet Evolution==<br />
<br />
At this point, if you have one final successful build, we should be ready to use the code to actually evolve the ice sheet thickness. To try some simple test cases with the higher-order model and ice sheet evolution using the 1st-order scheme, go here:<br />
<br />
* [[Ice Sheet Evolution Experiments]]</div>Meierbtwhttp://websrv.cs.umt.edu/isis/index.php/Adding_dH/dt_moduleAdding dH/dt module2009-08-13T22:04:36Z<p>Meierbtw: Changed line number after EVOL_INC_REMAP</p>
<hr />
<div>==Overview==<br />
<br />
This page contains step-by-step instructions for adding a new module to Glimmer-CISM. In this case, the module is a first-order, upwinding advection scheme for mass transport (dH/dt) using velocities calculated from a higher-order dynamics model. The procedure, however, is generic and could apply to adding almost any module. The goal here is not to get lost in the details of the upwinding scheme (we'll provide the code for that later) but rather to understand the steps for adding a module in an incremental, structured way that allows you to track small errors and fix them at each step. <br />
<br />
*To go directly to some exercises using the higher-order dynamics and the evolution scheme, go to [[Ice Sheet Evolution Experiments]].<br />
*For mathematical background on upwinded transport, consult the [[Solving the equation for thickness evolution]] pages.<br />
<br />
<br />
== '''Step 0''' ==<br />
<br />
We first need to download the code, configure it, and check that we can get a successful initial build (presumably we've done that already if you've been able to do the higher-order test suite exercises). If not, then type the following from within the highest level directory containing the source code:<br />
<br />
./bootstrap<br />
./configure --with-netcdf=/path/to_netcdf !!! Note: specify your own path to netCDF libs here !!!<br />
make<br />
<br />
== '''Step 1''' ==<br />
<br />
In step 1, we will create a module containing "empty" subroutines and add this module to the build.<br />
<br />
Create a module '''fo_upwind_advect.F90''' that contains the necessary subroutines. For now, these will just be "stubs", which we will fill in later.<br />
<br />
<source lang=fortran>module fo_upwind_advect<br />
<br />
! subroutines for mass advection scheme based on 1st order upwinding<br />
<br />
contains<br />
<br />
!----------------------------------------------------------------------<br />
<br />
subroutine fo_upwind_advect_init( )<br />
<br />
! initialization for 1st-order upwinding mass advection<br />
<br />
end subroutine fo_upwind_advect_init<br />
<br />
!----------------------------------------------------------------------<br />
<br />
subroutine fo_upwind_advect_final( )<br />
<br />
! finalization for 1st-order upwinding mass advection<br />
<br />
end subroutine fo_upwind_advect_final<br />
<br />
!----------------------------------------------------------------------<br />
<br />
subroutine fo_upwind_advect_driver( )<br />
<br />
! driver for 1st-order upwind mass advection<br />
<br />
end subroutine fo_upwind_advect_driver<br />
<br />
!----------------------------------------------------------------------<br />
<br />
subroutine fo_upwind_advect_main( )<br />
<br />
! 1st-order upwinding mass advection<br />
<br />
end subroutine fo_upwind_advect_main<br />
<br />
!----------------------------------------------------------------------<br />
<br />
end module fo_upwind_advect</source><br />
<br />
<br />
Add this new module to the build by editing '''Makefile.am''' in '''src/fortran/''' ...<br />
<br />
<source lang=fortran>libglide_a_SOURCES = glide.F90 glide_setup.F90 glide_types.F90 glide_temp.F90 \<br />
glide_bwater.F90 glide_deriv.F90 xls.F90 ice3d_lib.F90 \<br />
glide_velo_higher.F90 glide_thck.F90 glide_velo.F90 \<br />
glide_mask.F90 glide_stop.F90 glide_io.F90 \<br />
glide_nc_custom.F90 isostasy.F90 isostasy_el.F90\<br />
isostasy_setup.F90 isostasy_types.F90 glide_lithot.F90\<br />
glide_lithot3d.F90 glide_lithot1d.F90 glide_profile.F90\<br />
glide_diagnostics.F90 glissade.F90 glissade_remap.F90\<br />
glissade_velo.F90 glissade_constants.F90 glide_vertint.F90\<br />
glide_thckmask.F90 glide_nonlin.F90 glide_grids.F90\<br />
glam.F90 glam_strs2.F90 glam_thck_ppm.F90\<br />
remap_advection.F90 remap_glamutils.F90 glide_ground.F90</source><br />
<br />
becomes ...<br />
<br />
<source lang=fortran>libglide_a_SOURCES = glide.F90 glide_setup.F90 glide_types.F90 glide_temp.F90 \<br />
glide_bwater.F90 glide_deriv.F90 xls.F90 ice3d_lib.F90 \<br />
glide_velo_higher.F90 glide_thck.F90 glide_velo.F90 \<br />
glide_mask.F90 glide_stop.F90 glide_io.F90 \<br />
glide_nc_custom.F90 isostasy.F90 isostasy_el.F90\<br />
isostasy_setup.F90 isostasy_types.F90 glide_lithot.F90\<br />
glide_lithot3d.F90 glide_lithot1d.F90 glide_profile.F90\<br />
glide_diagnostics.F90 glissade.F90 glissade_remap.F90\<br />
glissade_velo.F90 glissade_constants.F90 glide_vertint.F90\<br />
glide_thckmask.F90 glide_nonlin.F90 glide_grids.F90\<br />
glam.F90 glam_strs2.F90 glam_thck_ppm.F90\<br />
remap_advection.F90 remap_glamutils.F90 glide_ground.F90 \ ! CHANGED LAST TWO LINES<br />
fo_upwind_advect.F90</source><br />
<br />
[Note: the comment in the 2nd to last line may cause a problem for make.]<br />
<br />
Now we stop adding/editing things for a minute and type '''make''' to rebuild and make sure we haven't broken the code yet. If you get a successful build, then proceed to '''Step 2'''.<br />
<br />
== '''Step 2''' ==<br />
<br />
In step 2, we will do some minor updates to other parts of the code so that there is an option to call our 1st order upwinding scheme when we want to evolve the ice thickness (that is, rather than using some other scheme).<br />
<br />
In '''glide_types.F90''', add the following to the transport scheme options around line 100 ...<br />
<br />
<source lang=fortran> integer, parameter :: EVOL_PSEUDO_DIFF = 0<br />
integer, parameter :: EVOL_ADI = 1<br />
integer, parameter :: EVOL_DIFFUSION = 2<br />
integer, parameter :: EVOL_INC_REMAP = 3<br />
integer, parameter :: EVOL_FO_UPWIND = 4 ! ADDED THIS LINE</source><br />
<br />
... similarly, update some lines around 207 for documentation ...<br />
<br />
<source lang=fortran> !*FD Thickness evolution method:<br />
!*FD \begin{description}<br />
!*FD \item[0] Pseudo-diffusion approach <br />
!*FD \item[2] Diffusion approach (also calculates velocities)<br />
!*FD \item[3] Incremental remapping<br />
!*FD \item[4] 1st-order upwinding scheme ! ADDED THIS LINE<br />
!*FD \end{description}</source><br />
<br />
In '''glide_setup.F90''', around line 475, change ...<br />
<br />
<source lang=fortran> character(len=*), dimension(0:4), parameter :: evolution = (/ & ! CHANGED 0:3 TO 0:4<br />
'pseudo-diffusion ', &<br />
'ADI scheme ', &<br />
'iterated diffusion ', &<br />
'remap thickness ', & <br />
'1st order upwind ' /) ! ADDED THIS LINE</source><br />
<br />
Time to type '''make''' again and check for a successful build. If so, proceed to '''Step 3'''.<br />
<br />
<br />
== '''Step 3''' ==<br />
<br />
Now we will add use statements, if constructs, and dummy calls to the advection subroutines from the appropriate places in the code. Hence, we are setting up all the necessary structure to call the subroutines but are not passing any arguments to them yet (nor are they actually doing anything yet).<br />
<br />
In '''glide.F90''', add the following use statement to subroutine '''glide_initialise''' ...<br />
<br />
<source lang=fortran> use glam_strs2, only : glam_velo_fordsiapstr_init<br />
use remap_glamutils, only : horizontal_remap_init<br />
use fo_upwind_advect, only : fo_upwind_advect_init ! ADDED<br />
</source><br />
<br />
Now add if construct and call to the initialization stub in '''fo_upwind_advect_mod.F90''', around line 217 of '''glide.F90''', just after the if construct for EVOL_INC_REMAP ...<br />
<br />
<source lang=fortran> if (model%options%whichevol== EVOL_FO_UPWIND ) then<br />
call fo_upwind_advect_init( )<br />
end if</source><br />
<br />
In '''glide_stop.F90''', add the following to the use statements ...<br />
<br />
<source lang=fortran>module glide_stop<br />
<br />
use glide_types<br />
use glimmer_log<br />
use remap_glamutils<br />
use fo_upwind_advect, only : fo_upwind_advect_final ! ADDED</source><br />
<br />
Add the necessary if construct to '''glide_stop.F90''' around line 150, as we did with the initialization routine ...<br />
<br />
<source lang=fortran> if (model%options%whichevol== EVOL_FO_UPWIND ) then<br />
call fo_upwind_advect_final( )<br />
endif</source><br />
<br />
Finally, update the use statements in '''glide.F90''' around line 340, and the case construct for ice sheet evolution, around line 380, so that the 1st order upwinding subroutines can be called ...<br />
<br />
<source lang=fortran> use glide_thck<br />
use glide_velo<br />
use glide_ground<br />
use glide_setup<br />
use glide_temp<br />
use glide_mask<br />
use isostasy<br />
use glam, only: inc_remap_driver<br />
use fo_upwind_advect, only: fo_upwind_advect_driver ! ADDED </source><br />
<br />
... and add a call to the driver routine to the case construct (passing no args yet) ...<br />
<br />
<source lang=fortran> case(EVOL_FO_UPWIND) ! Use first order upwind scheme for mass transport<br />
<br />
call fo_upwind_advect_driver( )</source><br />
<br />
Make sure to place it BEFORE the existing "end select" statement!<br />
<br />
Type '''make''' again and check for a successful build. If so, proceed to '''Step 4'''.<br />
<br />
== '''Step 4''' ==<br />
<br />
At this point we've done nothing but build the necessary structures in the code. Now we will start filling in the subroutines with the necessary variable definitions. We'll follow this by passing arguments and, if everything still works, we can actually start to do something w/ those arguments within the subroutines themselves.<br />
<br />
Note that, in practice, one doesn't always know ahead of time exactly what variables are needed within a particular subroutine, and there is some trial and error while figuring it out. Here, we'll assume we've thought this out really well ahead of time and we know exactly what we need to do before we ever started writing the code (like punch-card programmers of old), in which case we know exactly which arguments are needed when we start.<br />
<br />
<br />
First, we'll fill out the rest of the "driver" subroutine, which has two parts, (1) calling the higher-order dynamics subroutines '''run_ho_diagnostic''' to get the appropriate velocity fields and, (2) calling '''fo_upwind_advect_main''', which uses those velocity fields, the thickness field, and then moves mass around to calculate a new thickness field. In subroutine '''fo_upwind_advect_driver''' we add the following ...<br />
<br />
<source lang=fortran> subroutine fo_upwind_advect_driver( model ) ! ADDED ARG<br />
<br />
! driver routine for the 1st order, upwind mass transport scheme<br />
<br />
type(glide_global_type), intent(inout) :: model ! ADDED <br />
<br />
call run_ho_diagnostic(model) ! ADDED CALL AND ARG<br />
<br />
call fo_upwind_advect_main( )<br />
<br />
... </source><br />
<br />
Note that we also passed the model in to the driver and gave "model" a type definition. We also need to add the appropriate use statements for "model" and to access the higher-order dynamics routines, which we'll do below. <br />
<br />
Next, add argument definitions to the main subroutine in '''fo_upwind_advect.F90''' ...<br />
<br />
<source lang=fortran> subroutine fo_upwind_advect_main( thck, stagthck, acab, dt, uflx, vflx, ewn, nsn, dew, dns ) ! CHANGED<br />
<br />
! 1st-order upwinding mass advection<br />
<br />
implicit none ! ADDED<br />
<br />
real (kind = dp), intent(in) :: dt ! time step (sec)<br />
real (kind = dp), dimension(:,:), intent(inout) :: thck ! thickness on normal grid (m)<br />
real (kind = dp), dimension(:,:), intent(in) :: stagthck ! thickness on staggered grid (m)<br />
real (kind = sp), dimension(:,:), intent(in) :: acab ! surf mass balance (accum/ablation) (m/sec)<br />
real (kind = dp), dimension(:,:), intent(in) :: uflx, vflx ! flux in x,y directions (m^2/sec)<br />
real (kind = dp), intent(in) :: dew, dns ! grid spacing in x,y directions (m)<br />
integer, intent(in) :: ewn, nsn ! no. of grid points in x,y directions<br />
<br />
real (kind = dp) :: He, Hw, Hn, Hs, ue, uw, vn, vs ! ADDED<br />
integer :: ew, ns<br />
<br />
end subroutine fo_upwind_advect_main</source><br />
<br />
<br />
Because these definitions require kind "dp", we need to add a use statment at the start of the module. We've also added the necessary use statement for accessing "model" and the higher-order dynamics subroutines (from above) and a few more we'll use later on ...<br />
<br />
<source lang=fortran>module fo_upwind_advect<br />
<br />
! subroutines for mass advection scheme based on 1st order upwinding<br />
<br />
! ADDED<br />
use glimmer_paramets, only: sp, dp, len0, thk0, tim0, vel0, tim0, acc0, scyr ! scales for swapping between dim and non-dim vars<br />
use glide_types<br />
use glide_velo_higher<br />
<br />
private<br />
public :: fo_upwind_advect_init, fo_upwind_advect_driver, fo_upwind_advect_final<br />
<br />
... etc ... </source><br />
<br />
To avoid someone mucking around with variables in the module that we don't want them to touch, we've also added the "private" and "public" statements; the only parts of the module that can be accessed from outside are through calls to the three "public" subroutines. <br />
<br />
<br />
The arguments that we defined above in '''fo_upwind_advect_main''' need to be passed in from the driver routine '''fo_upwind_advect_driver''' to the subroutine that does all of the work, '''fo_upwind_advect_main'''. To do this, we must access them from the derived type "model". The entire driver subroutine then becomes ...<br />
<br />
<source lang=fortran> subroutine fo_upwind_advect_driver( model )<br />
<br />
! driver routine for the 1st order, upwind mass transport scheme<br />
<br />
type(glide_global_type), intent(inout) :: model<br />
<br />
call run_ho_diagnostic(model) ! get velocities and fluxes from HO dynamic subroutines<br />
<br />
call fo_upwind_advect_main( model%geometry%thck, model%geomderv%stagthck, &<br />
model%climate%acab, model%numerics%dt, &<br />
model%velocity_hom%uflx,model%velocity_hom%vflx, &<br />
model%general%ewn, model%general%nsn, &<br />
model%numerics%dew, model%numerics%dns )<br />
<br />
end subroutine fo_upwind_advect_driver</source><br />
<br />
<br />
Finally, we pass the derived type "model" during the call to the driver subroutine in '''glide.F90''' ...<br />
<br />
<source lang=fortran> case(EVOL_FO_UPWIND) ! Use first order upwind scheme for mass transport<br />
<br />
call fo_upwind_advect_driver( model ) ! ADDED ARGUMENT<br />
<br />
end select</source><br />
<br />
Type '''make''' again and check for a successful build. If so, proceed to '''Step 5'''.<br />
<br />
== '''Step 5''' ==<br />
<br />
Now let's fill in the init and finalization subroutines in '''fo_upwind_advect.F90'''. This takes some thought ahead of time and might require a bit of trial and error. Again, for now we'll assume that we've thought this out really well ahead of time and we know exactly what we need.<br />
<br />
First, declare any other necessary variables at the start of the module '''fo_upwind_advection.F90'''. Here, these are the allocatable work arrays that we'll use in '''fo_upwind_advect_main'''. The beginning of the module becomes ...<br />
<br />
<source lang=fortran>module fo_upwind_advect<br />
<br />
!----------------------------------------------------------------------<br />
<br />
! init, finalize, and driver subroutines for mass advection based on 1st order upwinding<br />
<br />
use glimmer_paramets, only: sp, dp, len0, thk0, tim0, vel0, tim0, acc0, scyr<br />
use glide_types<br />
use glide_velo_higher<br />
<br />
private<br />
public :: fo_upwind_advect_init, fo_upwind_advect_driver, fo_upwind_advect_final<br />
<br />
! allocatable work arrays ! ADDED<br />
real (kind = dp), allocatable, dimension(:,:) :: &<br />
ubar, vbar, &<br />
ubar_grid, vbar_grid, &<br />
flux_net, thck_grid, &<br />
mask, thck_old<br />
<br />
contains<br />
<br />
... etc ...</source><br />
<br />
Now change the initialization and finalization subroutines ...<br />
<br />
<source lang=fortran> subroutine fo_upwind_advect_init( ewn, nsn ) ! ADDED ARGS HERE<br />
<br />
! initialization for 1st-order upwinding mass advection<br />
<br />
implicit none ! ADDED<br />
<br />
integer, intent(in) :: ewn, nsn ! horizontal grid dimensions ! ADDED TYPE DEF<br />
<br />
integer :: errstat ! ADDED FOR ERROR HANDLING<br />
<br />
! allocate work arrays ! ADDED THESE<br />
allocate( ubar(ewn-1,nsn-1), stat=errstat ); ubar = 0.0_dp<br />
allocate( vbar(ewn-1,nsn-1), stat=errstat ); vbar = 0.0_dp<br />
allocate( ubar_grid(ewn+1,nsn+1), stat=errstat ); ubar_grid = 0.0_dp<br />
allocate( vbar_grid(ewn+1,nsn+1), stat=errstat ); vbar_grid = 0.0_dp<br />
allocate( thck_grid(ewn+2,nsn+2), stat=errstat ); thck_grid = 0.0_dp<br />
allocate( flux_net(ewn,nsn), stat=errstat ); flux_net = 0.0_dp<br />
allocate( mask(ewn,nsn), stat=errstat ); mask = 0.0_dp<br />
allocate( thck_old(ewn,nsn), stat=errstat ); thck_old = 0.0_dp<br />
<br />
if ( errstat /= 0 ) then ! ADDED FOR ERROR HANDLING<br />
print *, 'error: allocation in fo_upwind_advect failed!'<br />
stop<br />
end if<br />
<br />
end subroutine fo_upwind_advect_init</source><br />
<br />
Note that ewn, nsn are passed in above, so we need to make sure they are passed from the main code <br />
where this call sits. In '''glide.F90''' we have ...<br />
<br />
<source lang=fortran> call fo_upwind_advect_init( model%general%ewn, model%general%nsn ) ! ADDED ARGS</source><br />
<br />
As with the initialization subroutine, we add deallocation statements for work arrays in the finalization subroutine ...<br />
<br />
<source lang=fortran> subroutine fo_upwind_advect_final( )<br />
<br />
! finalization for 1st-order upwinding mass advection<br />
<br />
implicit none ! ADDED<br />
<br />
! deallocate work arrays ! ADDED THESE<br />
if( allocated( ubar ) ) deallocate( ubar )<br />
if( allocated( vbar ) ) deallocate( vbar )<br />
if( allocated( ubar_grid ) ) deallocate( ubar_grid )<br />
if( allocated( vbar_grid ) ) deallocate( vbar_grid )<br />
if( allocated( thck_grid ) ) deallocate( thck_grid )<br />
if( allocated( flux_net ) ) deallocate( flux_net )<br />
if( allocated( mask ) ) deallocate( mask )<br />
if( allocated( thck_old ) ) deallocate( thck_old )<br />
<br />
<br />
end subroutine fo_upwind_advect_final</source><br />
<br />
Type '''make''' again and check for a successful build. If so, proceed to '''Step 6'''.<br />
<br />
== '''Step 6''' ==<br />
<br />
Now we are actually ready to do something ... that is, fill in the guts of the new subroutine, where the thickness evolution calculation takes place. Rather than have you figure out how to code up the 1st-order upwinding scheme on your own, we'll provide you with the chunk of code to do that (below), noting that the details of the calculation scheme are discussed [[Solving the equation for thickness evolution|HERE]]. <br />
<br />
<source lang=fortran>! ----------------------------------<br />
<br />
subroutine fo_upwind_advect_main( thck, stagthck, acab, dt, uflx, vflx, ewn, nsn, dew, dns )<br />
<br />
! 1st-order upwinding mass advection that uses a finite-volume like scheme for <br />
! mass conservation. Velocities from the staggered grid (B-grid) are averaged onto the <br />
! faces of the non-staggered grid (i.e. faces of the grid where scalers like thickness live). <br />
! Thus, the averaged velocities exist on a C-grid, allowing mass transport to be treated <br />
! in a finite-volume manner; depth averaged velocities give the fluxes out of each cell <br />
! centered on a thickness point and the thickness advected is chosen according to upwinding.<br />
! <br />
! Note that this works at the calving front because a non-zero staggered thickness there <br />
! defines the velocities there. These velocites can be used to define the velocity at<br />
! the face of the last non-zero thickness cell (on the normal grid) which corresponds to<br />
! the location of the calving front. <br />
<br />
implicit none<br />
<br />
real (kind = dp), intent(in) :: dt<br />
real (kind = dp), dimension(:,:), intent(inout) :: thck<br />
real (kind = dp), dimension(:,:), intent(in) :: stagthck<br />
real (kind = sp), dimension(:,:), intent(in) :: acab<br />
real (kind = dp), dimension(:,:), intent(in) :: uflx, vflx<br />
real (kind = dp), intent(in) :: dew, dns<br />
integer, intent(in) :: ewn, nsn<br />
<br />
real (kind = dp) :: He, Hw, Hn, Hs, ue, uw, vn, vs ! upwinding variables and interface velocities<br />
<br />
integer :: ew, ns<br />
<br />
where( stagthck > 0.0_dp ) ! calculate the depth-ave velocities<br />
ubar = uflx / stagthck<br />
vbar = vflx / stagthck<br />
end where<br />
<br />
where( thck > 0.0_dp ) ! mask for eventually removing flux outside of the original domain<br />
mask = 1.0_dp ! (i.e. stuff that moves past the calving front goes away)<br />
else where<br />
mask = 0.0_dp<br />
end where<br />
<br />
thck_old = thck ! save the old thickness for debugging purposes<br />
<br />
! fill in the interior values on the extended velocity grid (extended B-grid)<br />
ubar_grid(2:ewn,2:nsn) = ubar<br />
vbar_grid(2:ewn,2:nsn) = vbar<br />
<br />
! fill in the interior values on the extended thickness grid<br />
thck_grid(2:ewn+1,2:nsn+1) = thck(:,:)<br />
<br />
! calculate the interface velocities from the extended B-grid, then use upwinding<br />
! criterion to advect thickness in or out of cells (NOTE that parts of this could<br />
! probably be vectorized at some point)<br />
do ns = 1, nsn<br />
do ew = 1, ewn<br />
<br />
! interface depth-ave velocities<br />
ue = ( ubar_grid(ew+1,ns+1) + ubar_grid(ew+1,ns) ) / 2.0d0<br />
uw = ( ubar_grid(ew,ns+1) + ubar_grid(ew,ns) ) / 2.0d0<br />
vn = ( vbar_grid(ew,ns+1) + vbar_grid(ew+1,ns+1) ) / 2.0d0<br />
vs = ( vbar_grid(ew,ns) + vbar_grid(ew+1,ns) ) / 2.0d0<br />
<br />
! choose thickness to advect based on upwinding<br />
if( ue > 0.0d0 )then<br />
He = - thck_grid(ew+1,ns+1) ! negative signs necessary so that flux to the east<br />
else ! results in mass loss in this volume (and vice versa)<br />
He = - thck_grid(ew+2,ns+1)<br />
end if<br />
if( uw > 0.0d0 )then<br />
Hw = thck_grid(ew,ns+1)<br />
else<br />
Hw = thck_grid(ew+1,ns+1)<br />
end if<br />
if( vn > 0.0d0 )then<br />
Hn = - thck_grid(ew+1,ns+1) ! negative signs here as above for ue, and He<br />
else<br />
Hn = - thck_grid(ew+1,ns+2)<br />
end if<br />
if( vs > 0.0d0 )then<br />
Hs = thck_grid(ew+1,ns)<br />
else<br />
Hs = thck_grid(ew+1,ns+1)<br />
end if<br />
<br />
! net flux into/out of each cell<br />
flux_net(ew,ns) = ( ue*He*dns + uw*Hw*dns + vn*Hn*dew + vs*Hs*dew )<br />
<br />
end do<br />
end do<br />
<br />
thck = thck_old + ( 1 / (dns * dew) * flux_net ) * dt + (acab * dt)<br />
<br />
! debugging<br />
print *, ' '<br />
print *, 'net volume change = ', sum( (thck-thck_old)*mask )*thk0 *dew*dns*len0**2<br />
print *, 'net calving flux = ', sum( thck * (1.0d0-mask) )*thk0*dew*dns*len0**2<br />
print *, '(for the confined shelf experiment, the above two should sum to ~0)'<br />
print *, 'mean accum/ablat rate = ', sum( acab * mask ) / sum(mask) / (dt*tim0) * scyr<br />
print *, 'mean dH/dt = ', sum( (thck-thck_old)*mask )*thk0 / sum(mask) / (dt*tim0) * scyr<br />
print *, 'sum of flux change (should be ~0) = ', sum( flux_net*vel0*thk0*len0 )<br />
print *, ' '<br />
! pause<br />
<br />
thck = thck * mask ! remove any mass advected outside of initial domain<br />
<br />
where( thck < 0.0_dp ) ! gaurd against thickness going negative<br />
thck = 0.0_dp<br />
end where<br />
<br />
end subroutine fo_upwind_advect_main<br />
<br />
! ----------------------------------</source><br />
<br />
==Ice Sheet Evolution==<br />
<br />
At this point, if you have one final successful build, we should be ready to use the code to actually evolve the ice sheet thickness. To try some simple test cases with the higher-order model and ice sheet evolution using the 1st-order scheme, go here:<br />
<br />
* [[Ice Sheet Evolution Experiments]]</div>Meierbtwhttp://websrv.cs.umt.edu/isis/index.php/PDX_afterhoursPDX afterhours2009-08-13T15:23:40Z<p>Meierbtw: /* Thursday Aug 13, 2009 */</p>
<hr />
<div>[[Image:Erin_adam.JPG|thumb|right|300px|The PDX Afterhours king & queen in Tube Bar, home of Wednesday night $1 Miller High Life.]]<br />
<br />
==tuesday Aug 4, 2009==<br />
* meet at [http://paccinirestaurant.com/ Paccini's] pub at 7:30<br />
<br />
==wednesday Aug 5, 2009==<br />
* after FORTRAN session, leave PSU to go to [http://maps.google.com/maps?hl=en&client=firefox-a&rls=com.ubuntu:en-US:unofficial&hs=zuU&um=1&ie=UTF-8&q=sushi+ichiban+portland&fb=1&split=1&gl=us&view=text&latlng=899833397715679289 Sushi Ichiban]. Adam Campbell will guide you.<br />
* go to [http://www.groundkontrol.com Ground Kontrol] a retro arcade with beer<br />
* go to [http://www.voodoodoughnut.com Voodoo Doughnut], please someone buy Ian Rutt the $5 doughnut.<br />
<br />
<br />
==thursday Aug 6, 2009==<br />
[http://amontobin.com/field/ Amon Tobin] and [http://www.pitchblack.co.nz/?s1=index Pitch Black], along with two opening bands, are playing at the Roseland Theatre (8 NW 6th Ave) Thursday night (starting at 9:00 pm or thereabouts). Tickets are $26 (available online at [http://ticketswest.rdln.com/Venue.aspx?ven=ROS TicketsWest]). The music is best described as sampled electronica (Amon Tobin) and Kiwi-style dub (Pitch Black). I'd expect a late night of electronic music: an evening nap may be in order! Jeremy's already got his ticket and can fill you in with more, including a sample of the music.<br />
<br />
*other local music recommendations from Adam<br />
<br />
'''Boy Eats Drum Machine, French Miami, Southern Belle and Electric Opera Company''' - Indie Rock -<br />
Thu., Aug. 6, 9 p.m.<br />
$6-8<br />
Berbati's Pan<br />
10 SW 3rd Ave.<br />
Downtown<br />
<br />
'''Nurses, Inside Voices and Slaves''' - Indie Rock -<br />
Thu., Aug. 6, 8:30 p.m.<br />
$7<br />
Holocene <br />
1001 SE Morrison<br />
Southeast<br />
<br />
==friday Aug 7, 2009==<br />
[http://www.biteoforegon.com/ The Bite of Oregon] is a food festival that takes place at Tom McCall Waterfront Park. Featuring food, wine, beer and entertainment from Oregon. Entry is $8, food and beverages are extra.<br />
<br />
[http://www.pioneercourthousesquare.org/calendar_august.htm Flicks on the Bricks] will be showing Jurassic Park at dusk outside at Pioneer Courthouse Square, FREE (including popcorn). (10 minute walk)<br />
<br />
[http://www.portlandtwilight.com/ Portland Downtown Twilight Bike Criterium]: Professional (by U.S. standards) bike race through downtown Portland. Complete with beer garden, food, and an expected 15,000 spectators. Pro race starts at 7:30.<br />
<br />
==saturday Aug 8, 2009==<br />
[http://www.biteoforegon.com/ The Bite of Oregon] is a food festival that takes place at Tom McCall Waterfront Park. Featuring food, wine, beer and entertainment from Oregon. Entry is $8, food and beverages are extra.<br />
<br />
==sunday Aug 9, 2009==<br />
The [http://providence.org/bridgepedal/ Portland Bridge Pedal] is a fun event where you can go on 14, 24, or 37 mi. bike ride over Portland's Bridges. Adam is trying to assemble a group to go on Sunday morning. Please speak with him if you are interested in going to this. I tentatively have 4 bikes I can get ahold of.<br />
<br />
RIDERS: Register [http://providence.org/bridgepedal/ here]. Make sure to click the team rider option and the 11 bridge option! Ok so here is the deal I tried to register a team:<br />
Team Name The Icepocalypse<br />
Team password BPT943<br />
<br />
HEY It works now!<br />
<br />
[http://www.biteoforegon.com/ The Bite of Oregon] is a food festival that takes place at Tom McCall Waterfront Park. Featuring food, wine, beer and entertainment from Oregon. Entry is $8, food and beverages are extra.<br />
<br />
==tuesday Aug 11, 2009==<br />
[http://www.rootsorganicbrewing.com/ Roots Brewery] has a $2.50 imperial pint night. The beer is great! It's on the eastside but it's not very far. There are carts and more beer nearby.<br />
<br />
==wednesday Aug 12, 2009==<br />
[[Image:420.png|thumb|right|300px|The Glimmer model output gives suggestions for illicit afterhours activities in Portland. I suggest looking in the North park blocks.]]<br />
<br />
==Thursday Aug 13, 2009==<br />
We attempted to make an expedition to the movie theater last Tuesday but I forgot my ID (:D). Tonight, we will make another try... So let's meet at the hotel at 6:30 pm. Don't forget to bring your ID!!!!!<br />
<br />
[http://www.freedomridersthemovie.com/ Freedom Riders] is a mountain bike documentary about the evolution of freeriding near Jackson Hole, WY (I think). The show is at the [http://www.clintonsttheater.com/ Clinton St. Theater]. It starts at 7 pm.</div>Meierbtwhttp://websrv.cs.umt.edu/isis/index.php/Student_Presentation_DevelopmentStudent Presentation Development2009-08-07T21:31:26Z<p>Meierbtw: /* Roundtable 1: Study regions */</p>
<hr />
<div>Students: Please email [[User:Mankoff|Ken Mankoff]] the Lat/Long of your institution and study area by the end of the day. (You can click 'email this user' in the lower left toolbox from his user page.<br />
<br />
==Roundtable 1: Study regions==<br />
* Antarctica. 6 students: [[User:Hoffman|Matt]], [[User:Adamc,|Adam]], [[User:mankoff|Ken Mankoff]]...<br />
* Greenland. 6 students: [[User:Kpoinar|Kristin]], [[User:meierbtw|Toby]]...<br />
* Mountain glaciers. 4 students: <br />
* Other/Global. 2 students:<br />
Possible discussion questions:<br />
* What questions are the climate change community pressuring us to answer?<br />
* What do we know now that would have been a big surprise 10 years ago?<br />
* How important is field data to your research? If you could collect any field data/observations to progress your work, what would it/they be? <br />
* <br />
*<br />
...<br />
<br />
==Roundtable 2: Background==<br />
* Geology. 5 students: [[User:meierbtw,|Toby]],...<br />
* Physics. 5 students: [[User:Adamc,|Adam]], [[User:Kpoinar|Kristin]]...<br />
* Math(s). 4 students:<br />
* Engineering/CS/Other. 4 students: [[User:Hoffman|Matt]], [[User:mankoff|Ken Mankoff]], ...<br />
Possible discussion questions:<br />
*<br />
* Given your background in XX, what have you done/would you do outside of your area of expertise to make yourself a better cryosphere scientist?<br />
* Why have we come to glaciology, over the other careers / topics available to people with our background?<br />
* Glaciology is 'hot' now... how hot will it stay? What are possible exit ramps once it cools down?<br />
* What should we do immediately post-PhD? travel, stay at home institution, ship off to postdoc abroad, ...<br />
* How can I continue to deal with people who say, "Oh, you study glaciers, huh? Better hurry up, hohnk hohnk"<br />
...<br />
<br />
==Other Ways to Sort Ourselves==<br />
Seneca - notes from yesterday?<br />
<br />
<br />
==Lists of valuable skills/knowledge learned, connections made, etc.?==</div>Meierbtwhttp://websrv.cs.umt.edu/isis/index.php/Student_Presentation_DevelopmentStudent Presentation Development2009-08-07T21:24:05Z<p>Meierbtw: /* Roundtable 2: Background */</p>
<hr />
<div>Students: Please email [[User:Mankoff|Ken Mankoff]] the Lat/Long of your institution and study area by the end of the day. (You can click 'email this user' in the lower left toolbox from his user page.<br />
<br />
==Roundtable 1: Study regions==<br />
* Antarctica. 6 students: [[User:Hoffman|Matt]], , [[User:Adamc,|Adam]],...<br />
* Greenland. 6 students: [[User:Kpoinar|Kristin]], [[User:meierbtw|Toby]]...<br />
* Mountain glaciers. 4 students: <br />
* Other/Global. 2 students:<br />
Possible discussion questions:<br />
* What questions are the climate change community pressuring us to answer?<br />
* What do we know now that would have been a big surprise 10 years ago?<br />
* <br />
* <br />
*<br />
...<br />
<br />
==Roundtable 2: Background==<br />
* Geology. 5 students: [[User:meierbtw,|Toby]],...<br />
* Physics. 5 students: [[User:Adamc,|Adam]],...<br />
* Math(s). 4 students:<br />
* Engineering/CS/Other. 4 students: [[User:Hoffman|Matt]], ...<br />
Possible discussion questions:<br />
*<br />
* Given your background in XX, what have you done/would you do outside of your area of expertise to make yourself a better cryosphere scientist?<br />
...<br />
<br />
==Other Ways to Sort Ourselves==<br />
Seneca - notes from yesterday?<br />
<br />
<br />
==Lists of valuable skills/knowledge learned, connections made, etc.?==</div>Meierbtwhttp://websrv.cs.umt.edu/isis/index.php/Student_Presentation_DevelopmentStudent Presentation Development2009-08-07T21:23:38Z<p>Meierbtw: /* Roundtable 1: Study regions */</p>
<hr />
<div>Students: Please email [[User:Mankoff|Ken Mankoff]] the Lat/Long of your institution and study area by the end of the day. (You can click 'email this user' in the lower left toolbox from his user page.<br />
<br />
==Roundtable 1: Study regions==<br />
* Antarctica. 6 students: [[User:Hoffman|Matt]], , [[User:Adamc,|Adam]],...<br />
* Greenland. 6 students: [[User:Kpoinar|Kristin]], [[User:meierbtw|Toby]]...<br />
* Mountain glaciers. 4 students: <br />
* Other/Global. 2 students:<br />
Possible discussion questions:<br />
* What questions are the climate change community pressuring us to answer?<br />
* What do we know now that would have been a big surprise 10 years ago?<br />
* <br />
* <br />
*<br />
...<br />
<br />
==Roundtable 2: Background==<br />
* Geology. 5 students:<br />
* Physics. 5 students: [[User:Adamc,|Adam]],...<br />
* Math(s). 4 students:<br />
* Engineering/CS/Other. 4 students: [[User:Hoffman|Matt]], ...<br />
Possible discussion questions:<br />
*<br />
* Given your background in XX, what have you done/would you do outside of your area of expertise to make yourself a better cryosphere scientist?<br />
...<br />
<br />
==Other Ways to Sort Ourselves==<br />
Seneca - notes from yesterday?<br />
<br />
<br />
==Lists of valuable skills/knowledge learned, connections made, etc.?==</div>Meierbtwhttp://websrv.cs.umt.edu/isis/index.php/Student_Presentation_DevelopmentStudent Presentation Development2009-08-07T21:22:58Z<p>Meierbtw: /* Roundtable 2: Background */</p>
<hr />
<div>Students: Please email [[User:Mankoff|Ken Mankoff]] the Lat/Long of your institution and study area by the end of the day. (You can click 'email this user' in the lower left toolbox from his user page.<br />
<br />
==Roundtable 1: Study regions==<br />
* Antarctica. 6 students: [[User:Hoffman|Matt]], , [[User:Adamc,|Adam]],...<br />
* Greenland. 6 students: [[User:Kpoinar|Kristin]], ...<br />
* Mountain glaciers. 4 students: <br />
* Other/Global. 2 students:<br />
Possible discussion questions:<br />
* What questions are the climate change community pressuring us to answer?<br />
* What do we know now that would have been a big surprise 10 years ago?<br />
* <br />
* <br />
*<br />
...<br />
<br />
==Roundtable 2: Background==<br />
* Geology. 5 students:<br />
* Physics. 5 students: [[User:Adamc,|Adam]],...<br />
* Math(s). 4 students:<br />
* Engineering/CS/Other. 4 students: [[User:Hoffman|Matt]], ...<br />
Possible discussion questions:<br />
*<br />
* Given your background in XX, what have you done/would you do outside of your area of expertise to make yourself a better cryosphere scientist?<br />
...<br />
<br />
==Other Ways to Sort Ourselves==<br />
Seneca - notes from yesterday?<br />
<br />
<br />
==Lists of valuable skills/knowledge learned, connections made, etc.?==</div>Meierbtwhttp://websrv.cs.umt.edu/isis/index.php/PDX_afterhoursPDX afterhours2009-08-07T01:21:45Z<p>Meierbtw: /* friday Aug 7, 2009 */</p>
<hr />
<div>[[Image:Erin_adam.JPG|thumb|right|300px|The PDX Afterhours king & queen in Tube Bar, home of Wednesday night $1 Miller High Life.]]<br />
<br />
==tuesday Aug 4, 2009==<br />
* meet at [http://paccinirestaurant.com/ Paccini's] pub at 7:30<br />
<br />
==wednesday Aug 5, 2009==<br />
* after FORTRAN session, leave PSU to go to [http://maps.google.com/maps?hl=en&client=firefox-a&rls=com.ubuntu:en-US:unofficial&hs=zuU&um=1&ie=UTF-8&q=sushi+ichiban+portland&fb=1&split=1&gl=us&view=text&latlng=899833397715679289 Sushi Ichiban]. Adam Campbell will guide you.<br />
* go to [http://www.groundkontrol.com Ground Kontrol] a retro arcade with beer<br />
* go to [http://www.voodoodoughnut.com Voodoo Doughnut], please someone buy Ian Rutt the $5 doughnut.<br />
<br />
<br />
==thursday Aug 6, 2009==<br />
[http://amontobin.com/field/ Amon Tobin] and [http://www.pitchblack.co.nz/?s1=index Pitch Black], along with two opening bands, are playing at the Roseland Theatre (8 NW 6th Ave) Thursday night (starting at 9:00 pm or thereabouts). Tickets are $26 (available online at [http://ticketswest.rdln.com/Venue.aspx?ven=ROS TicketsWest]). The music is best described as sampled electronica (Amon Tobin) and Kiwi-style dub (Pitch Black). I'd expect a late night of electronic music: an evening nap may be in order! Jeremy's already got his ticket and can fill you in with more, including a sample of the music.<br />
<br />
*other local music recommendations from Adam<br />
<br />
'''Boy Eats Drum Machine, French Miami, Southern Belle and Electric Opera Company''' - Indie Rock -<br />
Thu., Aug. 6, 9 p.m.<br />
$6-8<br />
Berbati's Pan<br />
10 SW 3rd Ave.<br />
Downtown<br />
<br />
'''Nurses, Inside Voices and Slaves''' - Indie Rock -<br />
Thu., Aug. 6, 8:30 p.m.<br />
$7<br />
Holocene <br />
1001 SE Morrison<br />
Southeast<br />
<br />
==friday Aug 7, 2009==<br />
[http://www.biteoforegon.com/ The Bite of Oregon] is a food festival that takes place at Tom McCall Waterfront Park. Featuring food, wine, beer and entertainment from Oregon. Entry is $8, food and beverages are extra.<br />
<br />
[http://www.pioneercourthousesquare.org/calendar_august.htm Flicks on the Bricks] will be showing Jurassic Park at dusk outside at Pioneer Courthouse Square, FREE (including popcorn). (10 minute walk)<br />
<br />
[http://www.portlandtwilight.com/ Portland Downtown Twilight Bike Criterium]: Professional (by U.S. standards) bike race through downtown Portland. Complete with beer garden, food, and an expected 15,000 spectators. Pro race starts at 7:30.<br />
<br />
==saturday Aug 8, 2009==<br />
[http://www.biteoforegon.com/ The Bite of Oregon] is a food festival that takes place at Tom McCall Waterfront Park. Featuring food, wine, beer and entertainment from Oregon. Entry is $8, food and beverages are extra.<br />
<br />
==sunday Aug 9, 2009==<br />
The [http://providence.org/bridgepedal/ Portland Bridge Pedal] is a fun event where you can go on 14, 24, or 37 mi. bike ride over Portland's Bridges. Adam is trying to assemble a group to go on Sunday morning. Please speak with him if you are interested in going to this. I tentatively have 4 bikes I can get ahold of.<br />
<br />
[http://www.biteoforegon.com/ The Bite of Oregon] is a food festival that takes place at Tom McCall Waterfront Park. Featuring food, wine, beer and entertainment from Oregon. Entry is $8, food and beverages are extra.</div>Meierbtwhttp://websrv.cs.umt.edu/isis/index.php/PDX_afterhoursPDX afterhours2009-08-07T01:20:55Z<p>Meierbtw: /* friday Aug 7, 2009 */</p>
<hr />
<div>[[Image:Erin_adam.JPG|thumb|right|300px|The PDX Afterhours king & queen in Tube Bar, home of Wednesday night $1 Miller High Life.]]<br />
<br />
==tuesday Aug 4, 2009==<br />
* meet at [http://paccinirestaurant.com/ Paccini's] pub at 7:30<br />
<br />
==wednesday Aug 5, 2009==<br />
* after FORTRAN session, leave PSU to go to [http://maps.google.com/maps?hl=en&client=firefox-a&rls=com.ubuntu:en-US:unofficial&hs=zuU&um=1&ie=UTF-8&q=sushi+ichiban+portland&fb=1&split=1&gl=us&view=text&latlng=899833397715679289 Sushi Ichiban]. Adam Campbell will guide you.<br />
* go to [http://www.groundkontrol.com Ground Kontrol] a retro arcade with beer<br />
* go to [http://www.voodoodoughnut.com Voodoo Doughnut], please someone buy Ian Rutt the $5 doughnut.<br />
<br />
<br />
==thursday Aug 6, 2009==<br />
[http://amontobin.com/field/ Amon Tobin] and [http://www.pitchblack.co.nz/?s1=index Pitch Black], along with two opening bands, are playing at the Roseland Theatre (8 NW 6th Ave) Thursday night (starting at 9:00 pm or thereabouts). Tickets are $26 (available online at [http://ticketswest.rdln.com/Venue.aspx?ven=ROS TicketsWest]). The music is best described as sampled electronica (Amon Tobin) and Kiwi-style dub (Pitch Black). I'd expect a late night of electronic music: an evening nap may be in order! Jeremy's already got his ticket and can fill you in with more, including a sample of the music.<br />
<br />
*other local music recommendations from Adam<br />
<br />
'''Boy Eats Drum Machine, French Miami, Southern Belle and Electric Opera Company''' - Indie Rock -<br />
Thu., Aug. 6, 9 p.m.<br />
$6-8<br />
Berbati's Pan<br />
10 SW 3rd Ave.<br />
Downtown<br />
<br />
'''Nurses, Inside Voices and Slaves''' - Indie Rock -<br />
Thu., Aug. 6, 8:30 p.m.<br />
$7<br />
Holocene <br />
1001 SE Morrison<br />
Southeast<br />
<br />
==friday Aug 7, 2009==<br />
[http://www.biteoforegon.com/ The Bite of Oregon] is a food festival that takes place at Tom McCall Waterfront Park. Featuring food, wine, beer and entertainment from Oregon. Entry is $8, food and beverages are extra.<br />
<br />
[http://www.pioneercourthousesquare.org/calendar_august.htm Flicks on the Bricks] will be showing Jurassic Park at dusk outside at Pioneer Courthouse Square, FREE (including popcorn). (10 minute walk)<br />
<br />
[http://www.portlandtwilight.com/ Portland Downtown Twilight Bike Criterium]: Professional (by U.S. standards) bike race through downtown Portland. Complete with beer garden, food, and an expected 15,000 spectators.<br />
<br />
==saturday Aug 8, 2009==<br />
[http://www.biteoforegon.com/ The Bite of Oregon] is a food festival that takes place at Tom McCall Waterfront Park. Featuring food, wine, beer and entertainment from Oregon. Entry is $8, food and beverages are extra.<br />
<br />
==sunday Aug 9, 2009==<br />
The [http://providence.org/bridgepedal/ Portland Bridge Pedal] is a fun event where you can go on 14, 24, or 37 mi. bike ride over Portland's Bridges. Adam is trying to assemble a group to go on Sunday morning. Please speak with him if you are interested in going to this. I tentatively have 4 bikes I can get ahold of.<br />
<br />
[http://www.biteoforegon.com/ The Bite of Oregon] is a food festival that takes place at Tom McCall Waterfront Park. Featuring food, wine, beer and entertainment from Oregon. Entry is $8, food and beverages are extra.</div>Meierbtwhttp://websrv.cs.umt.edu/isis/index.php/Finite_differencing_IFinite differencing I2009-08-05T15:19:43Z<p>Meierbtw: /* Final program and exercise */</p>
<hr />
<div>==Overview==<br />
<br />
In one dimension, the general form of the [http://en.wikipedia.org/wiki/Convection_diffusion_equation convection diffusion equation] is<br />
<br />
:<math>\frac{\partial u(x,t) }{\partial t} - \frac{\partial}{\partial x} D(x) \frac{\partial}{\partial x} u(x,t) - C(x)\frac{\partial}{\partial x} u(x,t) = S(x,t), </math><br />
<br />
<math>u</math> is a general variable, <math>D</math> is a spatially-varying diffusivity, <math>C</math> is a spatially-varying convection rate, and <math>S</math> is a source term. The second term on the left represents diffusion of a solute or other material property, the third term represent convection.<br />
<br />
This equation can be used to model a wide range of phenomena, including the distribution of temperatures (or energy conservation) in an ice sheet. It also bears similarity to the equations expressing conservation of momentum, and analysis of the numerical solutions to this equation are representative of the analysis of numerous numerical treatments in computational fluid dynamics. For these reasons, we will make convection diffusion our first "model problem", or problem to solve in order to strengthen intuition.<br />
<br />
We will take a step-wise approach to solving this equation, first solving for the diffusion, or parabolic equation, then solving for the convection portion. Finally, we will solve the complete equation. Through this process, we will be looking at the stability of the numerical methods used to solve the equation. <br />
<br />
You should think of this as a starting point for both your learning to program, as well as your learning to solve PDEs with programs.<br />
<br />
==Diffusion and explicit solution==<br />
First, we will solve a simplified version of the equation ''explicitly''. Explicit here refers to the way (or what) the differentiation operators are applied to. In this situation they are directly applied to the solution at the present time step in order to determine the next time. <br />
<br />
[[Image:Explicit_method-stencil.png|left|thumb|500 px|The Stencil for the most common explicit method for the parabolic equation.]]<br />
<br />
To better understand, apply the idea to what is called the parabolic, diffusion, or sometimes heat equation. In terms of convection diffusion this is <math>D(x,t) = 1</math>, <math>C(x,t) = 0</math> and <math>S(x,t) = 0</math>,<br />
<br />
:<math> \frac{\partial u(x,t) }{\partial t} = \frac{\partial ^2 u(x,t)}{\partial x^2},</math><br />
<br />
The finite difference approximation of the equation is<br />
<br />
:<math> \frac{u(x,t+\Delta t) - u(x,t)}{\Delta t} = \frac{u(x+\Delta x,t) - 2u(x,t) + u(x-\Delta x,t)}{\Delta x^2}.</math><br />
<br />
Where both derivative approximations are known from the previous lesson. One is called the 'forward Euler' approximation of the time derivative, and the other is the second order accurate, centered second derivative. <br />
<br />
The equation is then algebraically solved for <math> u(x,t)</math> <br />
<br />
:<math>u(x,t + \Delta t) = u(x,t) + \Delta t \frac{u(x+\Delta x,t) - 2u(x,t) +<br />
u(x-\Delta x,t)}{\Delta x^2}</math><br />
<br />
There you have it, a way to compute the future, using the present. Your first task will be to change this into an algorithm.<br />
<br />
One final note, the stencil at the left makes great sense, and understanding it will help make other algorithms clear. However, to understand it we must modify the notation a little. Make super-scripts <math>n</math> refer to time and subscripts <math>j</math> refer to space. The previous equation becomes<br />
<br />
: <math>u_j^{n+1} = u_j^n + \frac{\Delta t}{\Delta x ^2} \left( u_{j-1}^{n} - 2 u_j^{n} + u_{j+1}^{n} \right).</math><br />
<br />
Make sure you recognize how this corresponds with the diagram of the explicit stencil.<br />
<br />
==Numerical solution==<br />
One can see from the above equation that one way to numerically solve the parabolic equation is to use stencils or operators for computing second derivatives. Once the derivative is computed, finding the solution corresponding to the next time step, <math>u(x,t+\Delta t)</math> is just a matter of multiplying the derivative by <math>\Delta t</math> and adding <math> u(x,t)</math>. <br />
<br />
In psuedocode, the solution looks something like<br />
<br />
<source lang=text><br />
<br />
Initialize variables<br />
<br />
loop t over time:<br />
loop i over space:<br />
u(t,i) = u(t-1,i) + delta_t * (u(t-1,i-1) - 2*u(t-1,i) + u(t-1,i+1) ) / delta_x**2<br />
store solution as needed<br />
end loop over space<br />
end loop over time<br />
<br />
</source><br />
<br />
Note that in psuedo-code the final result is an T by N array (T is time steps, N space points). This is typical, the data structures used to store output is often as complex as the algorithm. Sometimes more. <br />
<br />
You'll need to do this in fortran 90, and you'll find Gethin's [[Pragmatic Programming]] very helpful. Plotting is also a key to understanding simulation output. Again, Gethin's got what you need, although you'll also find the [http://matplotlib.sourceforge.net/ matplotlib] documentation helpful (should be familiar if you've used Matlab).<br />
<br />
===Plotting results===<br />
The above line <br />
<source lang=text><br />
store solution as needed<br />
</source><br />
<br />
is perhaps, frustratingly vague. Gethin's [http://source.ggy.bris.ac.uk/wiki/Fortran1 Fortran1: Fortran for beginners] discusses input/output, IO, so the writing should not be a problem. Even<br />
<br />
<source lang=fortran><br />
write(*,*) u<br />
</source><br />
<br />
will suffice, if you can re-direct program output to a file. Assuming your program is called '''cd_prg''', this is done with<br />
<br />
<source lang=bash><br />
./cd_prg > data<br />
</source><br />
<br />
where the '''>''' operator will replace whatever was previously in the file '''data'''.<br />
<br />
As for reading it into Python and plotting it, the main challenges are reading in the data and animating the time series data. Here is some code to do that<br />
<br />
<source lang=python><br />
#!/usr/bin/env python<br />
<br />
# Import only what is needed<br />
from numpy import loadtxt,shape,linspace<br />
from pylab import plot,show,clf,show,ion<br />
<br />
# Import data file, called 'data'<br />
d=loadtxt('data')<br />
<br />
#Determine how much data came in<br />
dims = shape(d)<br />
<br />
clf() # Clears the screen<br />
ion() # Interactive plot mode, critical for animation<br />
<br />
# x data, note that this must correspond to program's domain<br />
x = linspace(0,1,dims[1]) <br />
<br />
# Initial plot, very Matlab(ish), note return of plot handle that allows plot to<br />
# be altered elsewhere in code.<br />
ph,=plot(x,d[0,:],'k') <br />
ph.figure.show() # matplot lib requires show to be called<br />
<br />
# Loop to plot each time step<br />
for i in range(1,dims[0]):<br />
ph.set_ydata(d[i,:]) # Only update y data (faster than replot)<br />
ph.figure.show()<br />
</source><br />
<br />
==Exercise==<br />
#Using the algorithm for the 'explicit' method, find a numerical solution to this heat conduction problem:<br />
<br />
:<math> \frac{\partial u(x,t) }{\partial t} = \frac{\partial ^2 u(x,t)}{\partial x^2},</math><br />
:<math> u(x,0) = sin(\pi x)</math><br />
:<math> u(0,t) = u(1,t) = 0</math><br />
<br />
Use <math>\Delta x</math> = 0.1, <math>\Delta t</math> = 0.005125, and <math>t_{end}</math> = 1.025. Compare the computed solution to the exact solution <math> u(x,t) = \exp(-\pi^2 t) \sin(\pi x)</math>. Repeat the experiment with <math>\Delta t</math> = 0.006 and <math>t_{end}</math> = 1.026.<br />
* [[Group one, parabolic, explicit]]<br />
* [[Group two, parabolic, explicit]]<br />
* [[Group three, parabolic, explicit]]<br />
* [[Group four, parabolic, explicit]]<br />
* [[Group five, parabolic, explicit]]<br />
* [[Group six, parabolic, explicit]]<br />
<br />
==Convection and numerical stability==<br />
For this unit consider the first-order hyperbolic PDE<br />
:<math> \frac{\partial u}{\partial t} = - v \frac{\partial u}{\partial x}.</math><br />
<br />
Mathematically, this statement is saying that a quantity <math>u(x,t)</math> exists on some grid, and is being carried along by a wind with velocity <math>v</math>. Before applying finite difference operators, clean up the notation so that super-scripts (<math>u^n</math>) refer to time and subscripts (<math>u_j</math>) refer to space.<br />
<br />
Continuing to work with explicit schemes, the machinery of descritization allows us to quickly move to the form<br />
:<math> \frac{u_j^{n+1} - u_j^n}{\Delta t} = -v \left( \frac{u^n_{j+1} - u_{j-1}^n}{2\Delta x} \right ) </math><br />
<br />
and solve to give a recurrence relation<br />
<br />
:<math> u_j^{n+1} =u_j^n - \frac{v \Delta t}{2\Delta x} \left(u^n_{j+1} - u_{j-1}^n \right ) </math><br />
<br />
===von Neumann Stability Analysis===<br />
Before implementing this, consider the stability of the solutions by assuming a very generic form of solution<br />
:<math> u_j^n = A(k)^n e^{ijk}</math>.<br />
<br />
Maybe you can recall a course in differential equations where you spent the better part of a semester making similar substitutions into equations to find solutions? This [[Wikipedia:Euler's formula | complex exponential]] is the Swiss Army knife of functions, and satisfies many equations.<br />
<br />
In our assumed solutions the amplitude is <math>A^n(k)</math> (exponentiated to higher powers with time) and the ''wave number'' is <math>k = \frac{2 \pi}{\lambda}</math>. Said in words, we assume that the solution will be oscillatory (recall <math> e^{ikx} = cos(kx) + i sin(kx)</math>) and that the solution's amplitude will depend on the frequency, or <math>k</math>. In our discrete case <math>j</math> serves as a proxy for space, <math>x</math>.<br />
<br />
:<math> A^{n+1} e^{ijk} = A^{n} e^{ijk} - \frac{v \Delta t}{2\Delta x} \left( A^{n} e^{i(j+1)k} - A^{n} e^{i(j-1)k} \right ) </math><br />
<br />
divide through by <math> A^n e^{ijk}</math><br />
<br />
:<math> A = 1 - \frac{v \Delta t}{2\Delta x} (e^{ik} - e^{-ik}) = 1 - i\frac{v \Delta t}{\Delta x} sin k</math><br />
<br />
''What does that mean?'' If <math>|A|^2>1</math>, the solution grows without bound in time, because each time step applies an higher exponent to <math> A(k)^n</math>. So, this solution scheme is unstable for all time steps, and space steps. Bummer.<br />
<br />
===Getting stability===<br />
Now, let's try that again, and when discretizing do what will become a favorite trick, average or smear the values of the function. A new discretization will be <br />
<br />
:<math> u_j^{n+1} = \frac{1}{2} (u_{j+1}^n + u_{j-1}^n) - \frac{v \Delta t}{2\Delta x} \left(u^n_{j+1} - u_{j-1}^n \right ), </math><br />
<br />
this is called the ''Lax method''. Now consider stability in the same way. Omitting some algebra<br />
:<math> A = cos k - i \frac{v \Delta t}{2\Delta x} sin k, </math><br />
is the amplitude. Requiring <br />
:<math> |A|^2 = cos^2 k + \left(\frac{v \Delta t}{2\Delta x}\right)^2 sin^2 k \leq 1, </math><br />
to avoid unbound growth, yields<br />
:<math> \frac{|v|\Delta t}{\Delta x} \leq 1.</math><br />
<br />
This is called the ''[[Wikipedia:Courantâ€“Friedrichsâ€“Lewy condition |Courant-Friedrichs-Levy]]'' stability criterion. It states that the information on a grid has a velocity of <math> \frac{\Delta x}{\Delta t} </math> and that the velocity in the system can not be exceeded by it (causing the ratio to exceed one). For such a thing to happen would be completely unphysical. Consider what happens when an object exceeds the velocity of waves in the media that carries it, a sonic boom. This is a "numerical boom".<br />
<br />
==Exercises==<br />
#Implement the Lax method for a linear system. Is this method explicit or implicit? Use a 10 unit domain and begin with a height 1.0 square wave between 4.5 <math>\leq x \leq</math> 5.5. Fix the ends at 0. let v be 1.0. Also track the sum of the solution before and after the simulation. End the simulation after 4 seconds. Report the behavior with and without the CFL being satisfied. If the CFL is very small, do things improve. What about when it's just under 1.0. What's going on here. Try subtracting <math>u^j_n</math> from both sides of the discretization and inspect for differences between the original discretization and Lax. See an extra term?<br />
#Try the ''leapfrog method'' for descretization<br />
:<math>u^{n+1}_j = u^{n-1}_j - \frac{v\Delta t}{\Delta x} (u^n_{j+1} - u^n_{j-1}).</math><br />
:Does this improve the numerical diffusion?<br />
<br />
* [[Group2 Convection Exercise]]<br />
<br />
==Final program and exercise==<br />
Bring the descretization schemes for diffusion and convection to solve the convection-diffusion equation explicitly, with finite differences. Consider a non-dimensional form of the equation.<br />
<br />
:<math>\frac{\partial \phi}{\partial t} + u \frac{\partial \phi}{\partial x} - \frac{\partial }{\partial x}\mathrm{Pe}^{-1} \frac{\partial }{\partial x} \phi = q </math><br />
<br />
The <math>\mathrm{Pe}^{-1} </math>, the inverse of the [[Wikipedia:Peclet number|Peclet number]], the ratio of the velocity scale <math>U</math> times the length scale <math>L</math> to the diffusivity <math>D</math>, <br />
:<math> \mathrm{Pe} = \frac{UL}{D}.</math><br />
*On a unit domain, specify <math>\phi_(t,0)</math>=a= 0, <math>\phi_(t,1)</math>=b=1, and Pe = 10.<br />
*Compare the solution to <math>c_a(x)</math>, the analytic solution of the equation, is<br />
:<math> c_a(x) = a + (b-a)*\frac{\exp((x-1)\mathrm{Pe}) - \exp(\mathrm{-Pe})}{1-\exp(\mathrm{-Pe})} </math><br />
* Experiment with the Peclet number and mesh resolution to determine how stable your numerical scheme is.<br />
<br />
* [[Group2 Convection + Diffusion Exercise]]</div>Meierbtwhttp://websrv.cs.umt.edu/isis/index.php/Student_BiosStudent Bios2009-07-31T21:05:40Z<p>Meierbtw: </p>
<hr />
<div>*[[User:Mankoff|Ken Mankoff]] will begin his PhD. this fall at UCSC and as such does not have a very well defined research topic. He will likely work on projects involving subglacial lakes and grounding lines. He is currently analyzing data from the terminal face of Pine Island Glacier, and oceanographic and sea ice data from the larger Amundsen Sea area.<br />
<br />
*[http://www.victoria.ac.nz/antarctic/people/jeremy-fyke/index.aspx Jeremy Fyke] is working on a PhD with the Antarctic Research Centre in Wellington, New Zealand. My project involves coupling an ice sheet model to an Earth System model 'of intermediate complexity' (the University of Victoria Earth System Climate Model) in order to have a go at simulating coupled climate/ice sheet interactions over millennial time scales.<br />
<br />
*[http://flo-colleoni.ifrance.com/ Florence Colleoni] will defend her Ph.D. in paleoclimate modeling at [http://www-lgge.obs.ujf-grenoble.fr/ LGGE] (Grenoble, Fr) in early September. She will then start a post-doctorate at the [http://www.cmcc.it/welcome-at-cmccs-web-site?set_language=en Centro Euro-Mediterraneo per i Cambiamenti Climatici] in Bologna (Italy) to couple the CISM Glimmer to the Earth System model composed of the AGCM of NCAR and of the OGCM NEMO. The final aim is to carry out transient paleoclimate simulations to understand and reproduce the interglacial/glacial transition mechanisms. This will be done in collaboration with NCAR. - My entire Ph.D. thesis is available [ftp://ftp-lgge.obs.ujf-grenoble.fr/pub/depot/florence/ here]-<br />
<br />
* [http://homepages.ucalgary.ca/~adhikars/ Surendra Adhikari] is currently in his second year of PhD at the University of Calgary, Canada. He is trying to develop a 3-D higher-order numerical ice-flow model applied to valley glaciers and alpine ice-fields. This HO-model will then be coupled to the traditional SIA-model to simulate the large ice sheets such as Greenland Ice Sheet.<br />
<br />
*[http://bigice.apl.washington.edu/people_poinar.html Kristin Poinar] is a second-year Ph.D. student at the University of Washington who is working on two "learning curve" ice sheet modelling projects. One is writing a thermal model to apply to the Greenland ice sheet, where surface lake drainages make basal thermodynamics interesting; the second is your standard model-perturbations-at-the-terminus study, on Petermann Glacier in NW Greenland.<br />
<br />
*[[User:adamc|Adam Campbell]] is entering a PhD program at the University of Washington in Fall 2009. I have just completed a Masters Degree in Geology at Portland State University where I examined the physics of the reaction of Crane Glacier to the disintegration of the Larsen B Ice Shelf using a steady state 2-D flow model with a basal sliding law. I am presently investigating structures on the Kamb Ice Shelf to determine if they were developed by a pinch and swell mechanism. I am also uncomfortable writing about myself in the third person.<br />
<br />
*[[User:papplega|Patrick Applegate]]: I am a glacial geomorphologist and geochronologist with a taste for modeling. My Ph. D. work involves the use of geomorphic process modeling to parse out the real meaning of cosmogenic exposure dates from moraines. I am asymptotically approaching the completion of my Ph. D. at Penn State. I'm attending the Summer School because I anticipate taking a new direction for my research in the near future.<br />
<br />
*[[User:hoffman|Matt Hoffman]] is in his fifth and final (?) year of a PhD at Portland State University. I am developing a spatially-distributed energy balance model for the glaciers of the McMurdo Dry Valleys, Antarctica. The glaciers of the Dry Valleys are near the threshold of melt during summer, such that sublimation and melt are of similar magnitude. I anticipate the Summer School will develop my skills as a modeler and help me think about the relationships between surface mass balance and ice dynamics.<br />
<br />
*[[Toby Meierbachtol]]: I am beginning my PhD at The University of Montana this fall. While still in the beginning stages, my research will likely be focused on the subglacial hydrology of the Greenland ice sheet and controls on sliding through direct borehole observations. Additionally, I anticipate a modeling component to my research that could include incorporating field findings to constrain model results, or investigating uncertainties in boundary conditions. The Summer School is a great way for me to jump in with both feet.</div>Meierbtw