CISM exercise III: add a module to CISM and run prognostic test case

From Interactive System for Ice sheet Simulation
Revision as of 13:14, 16 March 2011 by Sprice (Talk | contribs)

Jump to: navigation, search

Contents

Overview

This page contains step-by-step instructions for adding a new module to CISM. In this case, the module is a first-order, upwinding advection scheme for mass transport (or, for calculating temporal changes in the ice thickness, H(x,y)) using velocities calculated from a higher-order dynamics model. The procedure, however, is fairly 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 necessary steps for adding a module to CISM in an incremental, structured way that allows you to track small errors and fix them at each step in the process.

The general outline of operations in glide.F90, the main driver subroutine for thermodynamics in CISM, is as follows:

[1] initialize the model (read in configuration file options, read in any other input files, allocate arrays, etc.)
[2] run the model
    <start time stepping loop>
        (a) calculate the temperature
        (b) calculate the velocity field and evolve the ice thickness
        (c) calculate isostatic adjustment
    <stop time stepping loop>
[3] finalize the model (write model output, de-allocate arrays, etc.)


The new module that we want to add will contain initialization, "driver", and finalization subroutines (somewhat standard coding practice, at least for the Community Earth System Model). Thus we will add the relevant subroutine calls in steps [1], [2b], and [3] (note that for now, we are assuming isothermal conditions and are ignoring isostatic adjustment so that steps [2a] and [2c] are left alone). The subroutine calls will access these subroutines through appropriate "use" statements (more on that below).

For mathematical background on "upwinded" transport, see the page on solving the equation for thickness evolution.

Step 0

We first need to download the code, configure it, and build it. Presumably we've done that already if you've been able to do the higher-order test suite exercises. If not, then follow the steps outlined here.

Step 1

In step 1, we will create a module containing "empty" subroutines and add this module to the build.

Using your favorite editor, create a module (for now, just a text file) fo_upwind_advect.F90 that contains the necessary subroutines. These subroutines will just be "stubs" for starters. We will fill them in later. Your fo_upwind_advect.F90 text file should contain the following:

module fo_upwind_advect
 
! subroutines for mass advection scheme based on 1st order upwinding
 
    contains
 
!----------------------------------------------------------------------
 
    subroutine fo_upwind_advect_init( )
 
    ! initialization for 1st-order upwinding mass advection
 
    end subroutine fo_upwind_advect_init
 
!----------------------------------------------------------------------
 
    subroutine fo_upwind_advect_final( )
 
    ! finalization for 1st-order upwinding mass advection
 
    end subroutine fo_upwind_advect_final
 
!----------------------------------------------------------------------
 
   subroutine fo_upwind_advect_driver( )
 
    ! driver for 1st-order upwind mass advection
 
   end subroutine fo_upwind_advect_driver
 
!----------------------------------------------------------------------
 
    subroutine fo_upwind_advect_main( )
 
    ! 1st-order upwinding mass advection
 
    end subroutine fo_upwind_advect_main
 
!----------------------------------------------------------------------
 
end module fo_upwind_advect

Note that the module contains all of the subroutines that we'll need, including an initialization subroutine (here, mainly for allocation of "work" arrays - those arrays not passed in or out but needed for calculations internal to the subroutine), a "driver" subroutine (called from glide.F90 - the one that will do all the work and interact with the larger body of code) holding the new subroutine that will do the thickness evolution calculation (that is, "advect" mass around the domain according to the velocity field), and a finalization routine (again, here mainly for de-allocation of work arrays).

Add this new module to the build by editing Makefile.am in libglide/ ...

libglide_la_SOURCES = glide.F90 glide_setup.F90 glide_types.F90 glide_temp.F90 \
                      glide_temp_utils.F90 glide_thck.F90 glide_velo.F90 glide_mask.F90 \
                      glide_stop.F90 glide_io.F90 glide_lithot_io.F90 glide_nc_custom.F90 \
                      isostasy.F90 isostasy_el.F90 isostasy_setup.F90 \
                      isostasy_types.F90 glide_lithot.F90 glide_lithot3d.F90 \
                      glide_lithot1d.F90 glide_profile.F90 glide_diagnostics.F90 \
                      glimmer_routing.F90 glimmer_ncio.F90 glimmer_ncparams.F90 \
                      glam.F90 glam_strs2.F90 \
                      glide_deriv.F90 glide_mask.F90 glide_nonlin.F90 \
                      glide_grids.F90 glide_stress.F90 glide_ground.F90 glide_vertint.F90\
                      glide_thckmask.F90 glide_velo_higher.F90 xls.f90 \
                      remap_advection.F90 remap_glamutils.F90 \
                      glide_mask.inc glide_nan.inc glide_bwater.F90 \
                      glimmer_restart_gcm.F90 glissade_temp.F90 \
                      glam_Basal_Proc.F90 fgmresD.f90

becomes ...

libglide_la_SOURCES = glide.F90 glide_setup.F90 glide_types.F90 glide_temp.F90 \
                      glide_temp_utils.F90 glide_thck.F90 glide_velo.F90 glide_mask.F90 \
                      glide_stop.F90 glide_io.F90 glide_lithot_io.F90 glide_nc_custom.F90 \
                      isostasy.F90 isostasy_el.F90 isostasy_setup.F90 \
                      isostasy_types.F90 glide_lithot.F90 glide_lithot3d.F90 \
                      glide_lithot1d.F90 glide_profile.F90 glide_diagnostics.F90 \
                      glimmer_routing.F90 glimmer_ncio.F90 glimmer_ncparams.F90 \
                      glam.F90 glam_strs2.F90 \
                      glide_deriv.F90 glide_mask.F90 glide_nonlin.F90 \
                      glide_grids.F90 glide_stress.F90 glide_ground.F90 glide_vertint.F90\
                      glide_thckmask.F90 glide_velo_higher.F90 xls.f90 \
                      remap_advection.F90 remap_glamutils.F90 fo_upwind_advect.F90 \
                      glide_mask.inc glide_nan.inc glide_bwater.F90 \
                      glimmer_restart_gcm.F90 glissade_temp.F90 \
                      glam_Basal_Proc.F90 fgmresD.f90

Note that the new module file fo_upwind_advect.F90 has been added to the 4th to the last line.

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.

Step 2

In step 2, we will do some minor updates to other parts of the code so that there is a run time option available to call our 1st-order upwinding scheme from a ".config" file.

In glide_types.F90, add the following to the transport scheme options around line 100 ...

  integer, parameter :: EVOL_PSEUDO_DIFF = 0
  integer, parameter :: EVOL_ADI = 1
  integer, parameter :: EVOL_DIFFUSION = 2
  integer, parameter :: EVOL_INC_REMAP = 3
  integer, parameter :: EVOL_FO_UPWIND = 4         !! ADDED THIS LINE

... similarly, update some lines around 212 for documentation (grep for "evolution" ) ...

    !*FD Thickness evolution method:
    !*FD \begin{description}
    !*FD \item[0] Pseudo-diffusion approach 
    !*FD \item[2] Diffusion approach (also calculates velocities)
    !*FD \item[3] Incremental remapping
    !*FD \item[4] 1st-order upwinding scheme       ! ADDED THIS LINE
    !*FD \end{description}

In glide_setup.F90, around line 493, change (grep for "evolution") ...

    character(len=*), dimension(0:4), parameter :: evolution = (/ &    !! CHANGED 0:3 TO 0:4
         'pseudo-diffusion                      ', &
         'ADI scheme                            ', &
         'iterated diffusion                    ', &
         'remap thickness                       ', &  
         '1st order upwind                      ' /) 			!! ADDED THIS LINE

Type make again and check for a successful build. If so, proceed to Step 3.

Step 3

Now we will add use statements, if constructs, and dummy calls to the advection subroutines from the appropriate places in the code. The use statements tell other parts of the code (e.g. glide.F90) where to find the subroutines that we call from there. The if constructs surround the various subroutine calls so that they are only activated when a flag is triggered in the .config file. The dummy calls let us test that all of the communication between the various parts of the code are working before we actually try to do anything. Thus, we are setting up all of the necessary structure to access the subroutines in fo_upwind_avect.F90 but are not passing any arguments or asking them to do anything yet.

In libglide/ glide.F90, add the following use statement to subroutine glide_initialise ...

    use glam_strs2, only : glam_velo_fordsiapstr_init
    use remap_glamutils, only : horizontal_remap_init
    use fo_upwind_advect, only : fo_upwind_advect_init         ! ADDED

Now add an if construct and dummy 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 ...

    if (model%options%whichevol== EVOL_FO_UPWIND ) then
        call fo_upwind_advect_init( )
    end if

In glide_stop.F90, add the following to the use statements ...

module glide_stop
 
  use glide_types
  use glimmer_log
  use remap_glamutils
  use fo_upwind_advect, only : fo_upwind_advect_final       ! ADDED

Add the necessary if construct to glide_stop.F90 around line 150, as we did with the initialization routine ...

    if (model%options%whichevol== EVOL_FO_UPWIND ) then
        call fo_upwind_advect_final( )
    endif

Finally, update the use statements in glide_tstep_p2 (the ice dynamics step) in glide.F90 around line 485 ...

    use glide_thck
    use glide_velo
    use glide_ground
    use glide_setup
    use glide_temp
    use glide_mask
    use isostasy
    use glam, only: inc_remap_driver
    use fo_upwind_advect, only: fo_upwind_advect_driver      ! ADDED

... and the case construct for ice sheet evolution, around line 520 (so that the 1st order upwinding subroutines can be called from here):

    case(EVOL_FO_UPWIND) ! Use first order upwind scheme for mass transport
 
       call fo_upwind_advect_driver( )

Make sure to place it BEFORE the existing "end select" statement!

Type make again and check for a successful build. If so, proceed to Step 4.

Step 4

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 with those arguments within the subroutines themselves.

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.

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 (which live in glide_velo_higher.F90) 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 the domain accordingly in order to calculate a new thickness field. In subroutine fo_upwind_advect_driver we need to add the following ...

subroutine fo_upwind_advect_driver( model )      ! ADDED ARG
 
    ! driver routine for the 1st order, upwind mass transport scheme
 
        type(glide_global_type), intent(inout) :: model        ! ADDED 
 
        call run_ho_diagnostic(model)          ! ADDED CALL AND ARG
 
        call fo_upwind_advect_main(  )          ! ADDED CALL 
 
        ...

Note that we also passed the derived type "model" (a data structure containing nearly all of the necessary fields in the model) to the driver and added a type definition for "model". We also need to add the appropriate use statements for "model" and to access the higher-order dynamics routines, which we'll do below.

Next, we will add argument definitions to the main subroutine in fo_upwind_advect.F90 ...

subroutine fo_upwind_advect_main( thck, stagthck, acab, dt, uflx, vflx, ewn, nsn, dew, dns )  ! CHANGED
 
    ! 1st-order upwinding mass advection
 
    implicit none							! ADDED
 
    real (kind = dp), intent(in) :: dt					 ! time step (sec)
    real (kind = dp), dimension(:,:), intent(inout) :: thck             ! thickness on normal grid (m)
    real (kind = dp), dimension(:,:), intent(in) :: stagthck            ! thickness on staggered grid (m)
    real (kind = sp), dimension(:,:), intent(in) :: acab                ! surf mass balance (accum/ablation) (m/sec)
    real (kind = dp), dimension(:,:), intent(in) :: uflx, vflx          ! flux in x,y directions (m^2/sec)
    real (kind = dp), intent(in) :: dew, dns                            ! grid spacing in x,y directions (m)
    integer, intent(in)  :: ewn, nsn                                    ! no. of grid points in x,y directions
 
    real (kind = dp) :: He, Hw, Hn, Hs, ue, uw, vn, vs		! ADDED
    integer :: ew, ns
 
end subroutine fo_upwind_advect_main

Because these definitions require kind "dp" we need to add a use statement 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 ...

module fo_upwind_advect
 
! subroutines for mass advection scheme based on 1st order upwinding
 
    ! ADDED
    use glimmer_paramets, only: sp, dp, len0, thk0, tim0, vel0, tim0, acc0, scyr       ! scales for swapping between dim and non-dim vars
    use glide_types
    use glide_velo_higher
 
    private
    public :: fo_upwind_advect_init, fo_upwind_advect_driver, fo_upwind_advect_final
 
    ... etc ...

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.

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 ...

    subroutine fo_upwind_advect_driver( model )
 
    ! driver routine for the 1st order, upwind mass transport scheme
 
        type(glide_global_type), intent(inout) :: model
 
        ! Compute the new geometry derivatives for this time step, since the thickness has been updated
        call geometry_derivs(model)
        call geometry_derivs_unstag(model)
 
        ! print some useful information to the screen 
        print *, ' '
        print *, '(dH/dt using first-order upwinding)'
        print *, 'time = ', model%numerics%time
 
        call run_ho_diagnostic(model)   ! get velocities and fluxes from HO dynamic subroutines
 
        call fo_upwind_advect_main( model%geometry%thck,    model%geomderv%stagthck,    &
                                    model%climate%acab,     model%numerics%dt,          &
                                    model%velocity%uflx,    model%velocity%vflx,    &
                                    model%general%ewn,      model%general%nsn,          &
                                    model%numerics%dew,     model%numerics%dns )
 
    end subroutine fo_upwind_advect_driver


Finally, we pass the derived type "model" during the call to the driver subroutine in glide.F90 ...

    case(EVOL_FO_UPWIND) ! Use first order upwind scheme for mass transport
 
       call fo_upwind_advect_driver( model )        ! ADDED ARGUMENT
 
    end select

Cross your fingers, type make again, and hope for a successful build. If so, proceed to Step 5.

Step 5

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.

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 ...

module fo_upwind_advect
 
!----------------------------------------------------------------------
 
    ! init, finalize, and driver subroutines for mass advection based on 1st order upwinding
 
    use glimmer_paramets, only: sp, dp, len0, thk0, tim0, vel0, tim0, acc0, scyr
    use glide_types
    use glide_velo_higher
    use glide_thck
 
 
    private
    public :: fo_upwind_advect_init, fo_upwind_advect_driver, fo_upwind_advect_final
 
    ! allocatable work arrays                                                    ! ADDED
    real (kind = dp), allocatable, dimension(:,:) ::    &
                       ubar, vbar,                      &
                       ubar_grid, vbar_grid,            &
                       flux_net, thck_grid,             &
                       mask, thck_old
 
    contains
 
     ... etc ...

Now change the initialization and finalization subroutines ...

subroutine fo_upwind_advect_init( ewn, nsn )                  ! ADDED ARGS HERE
 
    ! initialization for 1st-order upwinding mass advection
 
    implicit none                    ! ADDED
 
    integer, intent(in) :: ewn, nsn   ! horizontal grid dimensions                  ! ADDED TYPE DEF
 
    integer :: errstat                                                              ! ADDED FOR ERROR HANDLING
 
    ! allocate work arrays                                                          ! ADDED THESE
    allocate( ubar(ewn-1,nsn-1), stat=errstat ); ubar = 0.0_dp
    allocate( vbar(ewn-1,nsn-1), stat=errstat ); vbar = 0.0_dp
    allocate( ubar_grid(ewn+1,nsn+1), stat=errstat ); ubar_grid = 0.0_dp
    allocate( vbar_grid(ewn+1,nsn+1), stat=errstat ); vbar_grid = 0.0_dp
    allocate( thck_grid(ewn+2,nsn+2), stat=errstat ); thck_grid = 0.0_dp
    allocate( flux_net(ewn,nsn), stat=errstat ); flux_net = 0.0_dp
    allocate( mask(ewn,nsn), stat=errstat ); mask = 0.0_dp
    allocate( thck_old(ewn,nsn), stat=errstat ); thck_old = 0.0_dp
 
    if ( errstat /= 0 ) then                                                        ! ADDED FOR ERROR HANDLING
        print *, 'error: allocation in fo_upwind_advect failed!'
        stop
    end if
 
end subroutine fo_upwind_advect_init

Note that ewn, nsn are passed in above, so we need to make sure they are passed from the main code where this call sits. In glide.F90 we add these arguments to the initialization subroutine ...

     call fo_upwind_advect_init( model%general%ewn, model%general%nsn )   ! ADDED ARGS

Analogous to the allocation statements in the initialization subroutine, we add de-allocation statements for work arrays in the finalization subroutine ...

subroutine fo_upwind_advect_final( )
 
    ! finalization for 1st-order upwinding mass advection
 
    implicit none                           ! ADDED
 
    ! deallocate work arrays                    ! ADDED THESE
    if( allocated( ubar ) ) deallocate( ubar )
    if( allocated( vbar ) ) deallocate( vbar )
    if( allocated( ubar_grid ) ) deallocate( ubar_grid )
    if( allocated( vbar_grid ) ) deallocate( vbar_grid )
    if( allocated( thck_grid ) ) deallocate( thck_grid )
    if( allocated( flux_net ) ) deallocate( flux_net )
    if( allocated( mask ) ) deallocate( mask )
    if( allocated( thck_old ) ) deallocate( thck_old )
 
 
end subroutine fo_upwind_advect_final

Note that this subroutine does not require any arguments to be passed in or out.

Type make again and check for a successful build. If so, proceed to Step 6. Almost done now ...

Step 6

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 HERE.

! ----------------------------------
 
subroutine fo_upwind_advect_main( thck, stagthck, acab, dt, uflx, vflx, ewn, nsn, dew, dns )
 
    ! 1st-order upwinding mass advection that uses a finite-volume like scheme for 
    ! mass conservation. Velocities from the staggered grid (B-grid) are averaged onto the 
    ! faces of the non-staggered grid (i.e. faces of the grid where scalers like thickness live).  
    ! Thus, the averaged velocities exist on a C-grid, allowing mass transport to be treated  
    ! in a finite-volume manner; depth averaged velocities give the fluxes out of each cell 
    ! centered on a thickness point and the thickness advected is chosen according to upwinding.
    ! 
    ! Note that this works at the calving front because a non-zero staggered thickness there 
    ! defines the velocities there. These velocites can be used to define the velocity at
    ! the face of the last non-zero thickness cell (on the normal grid) which corresponds to
    ! the location of the calving front. 
 
    implicit none
 
    real (kind = dp), intent(in) :: dt
    real (kind = dp), dimension(:,:), intent(inout) :: thck
    real (kind = dp), dimension(:,:), intent(in) :: stagthck
    real (kind = sp), dimension(:,:), intent(in) :: acab
    real (kind = dp), dimension(:,:), intent(in) :: uflx, vflx
    real (kind = dp), intent(in) :: dew, dns
    integer, intent(in)  :: ewn, nsn
 
    real (kind = dp) :: He, Hw, Hn, Hs, ue, uw, vn, vs  ! upwinding variables and interface velocities
 
    integer :: ew, ns
 
    where( stagthck > 0.0_dp )  ! calculate the depth-ave velocities
        ubar = uflx / stagthck
        vbar = vflx / stagthck
    end where
 
    where( thck > 0.0_dp )      ! mask for eventually removing flux outside of the original domain
        mask = 1.0_dp           ! (i.e. stuff that moves past the calving front goes away)
    else where
        mask = 0.0_dp
    end where
 
    thck_old = thck             ! save the old thickness for debugging purposes
 
    ! fill in the interior values on the extended velocity grid (extended B-grid)
    ubar_grid(2:ewn,2:nsn) = ubar
    vbar_grid(2:ewn,2:nsn) = vbar
 
    ! fill in the interior values on the extended thickness grid
    thck_grid(2:ewn+1,2:nsn+1) = thck(:,:)
 
    ! calculate the interface velocities from the extended B-grid, then use upwinding
    ! criterion to advect thickness in or out of cells (NOTE that parts of this could
    ! probably be vectorized at some point)
    do ns = 1, nsn
        do ew = 1, ewn
 
            ! interface depth-ave velocities
            ue = ( ubar_grid(ew+1,ns+1) + ubar_grid(ew+1,ns) ) / 2.0d0
            uw = ( ubar_grid(ew,ns+1) + ubar_grid(ew,ns) ) / 2.0d0
            vn = ( vbar_grid(ew,ns+1) + vbar_grid(ew+1,ns+1) ) / 2.0d0
            vs = ( vbar_grid(ew,ns) + vbar_grid(ew+1,ns) ) / 2.0d0
 
            ! choose thickness to advect based on upwinding
            if( ue > 0.0d0 )then
                He = - thck_grid(ew+1,ns+1) ! negative signs necessary so that flux to the east
            else                            ! results in mass loss in this volume (and vice versa)
                He = - thck_grid(ew+2,ns+1)
            end if
            if( uw > 0.0d0 )then
                Hw = thck_grid(ew,ns+1)
            else
                Hw = thck_grid(ew+1,ns+1)
            end if
            if( vn > 0.0d0 )then
                Hn = - thck_grid(ew+1,ns+1) ! negative signs here as above for ue, and He
            else
                Hn = - thck_grid(ew+1,ns+2)
            end if
            if( vs > 0.0d0 )then
                Hs = thck_grid(ew+1,ns)
            else
                Hs = thck_grid(ew+1,ns+1)
            end if
 
            ! net flux into/out of each cell
            flux_net(ew,ns) = ( ue*He*dns + uw*Hw*dns + vn*Hn*dew + vs*Hs*dew )
 
        end do
    end do
 
    thck = thck_old + ( 1 / (dns * dew) * flux_net ) * dt + (acab * dt)
 
    ! debugging
    print *, ' '
    print *, 'net volume change = ', sum( (thck-thck_old)*mask )*thk0 *dew*dns*len0**2
    print *, 'net calving flux = ', sum( thck * (1.0d0-mask) )*thk0*dew*dns*len0**2
    print *, '(for the confined shelf experiment, the above two should sum to ~0)'
    print *, 'mean accum/ablat rate = ', sum( acab * mask ) / sum(mask) / (dt*tim0) * scyr
    print *, 'mean dH/dt = ', sum( (thck-thck_old)*mask )*thk0 / sum(mask) / (dt*tim0) * scyr
    print *, 'sum of flux change (should be ~0) = ', sum( flux_net*vel0*thk0*len0 )
    print *, ' '
!    pause
 
    thck = thck * mask               ! remove any mass advected outside of initial domain
 
    where( thck < 0.0_dp )           ! guard against thickness going negative
        thck = 0.0_dp
    end where
 
end subroutine fo_upwind_advect_main
 
! ----------------------------------

If you are having problems, you can copy the entire module containing the initialization, driver, and finalization subroutines into a new, blank fo_upwind.F90 file. Those can be found here:

In a dire emergency, you can untar a new version of the code that already has the first-order upwinding code integrated in it and copy all of the ".F90" files and the "Makefile.am" from the libglide/ directory into your own, followed by a "make clean; make". To get that tarball, click on the following link:

Evolution Experiments

Now the moment of truth. If you have had one last successful build, then it is time to actually try the new thickness evolution in some experiments. For this, we will return to our "confined-shelf" test case in tests/higher-order/shelf/. Move back into that directory and copy or link to your new executable file.

In confined-shelf.config, first set the option for evolution so that your new evolution module and subroutines are triggered,

[options]
evolution = 4      # was 3, now set to 4 to use new upwinding scheme

Also, for time stepping to take place, you will need to set the start and end times differently,

[time]
tstart = 0.
tend = 10.         # was previously set to 0
dt = 1.

It may take some playing around to figure out the "stable" value for dt, the time step. As discussed previously, this is an explicit scheme and so subject to the advective CFL condition.

Note also that the input file, output/confined-shelf.nc has a nonzero surface mass balance (field acab). With these changes, the ice shelf will evolve over time until its thickness and velocity are in steady-state with the specified boundary conditions and accumulation rate.

Now lets run the model forward in time for a few time steps (10 yrs as we specified in the configuration file) to see what happens.

python confined-shelf.py

Note that if you have a old executable file in this directory and you try to run the test case with it you will get an error; the old executable does not know what the .config file option evolution = 4 means. This is why in most cases it is better to set up a virtual link to the executable file rather than copy it (the link always points to the most recent version of the executable file).

Assuming the code is running, you can now sit back and watch the residuals spool by as the domain evolves through time. After the short 10 year run you can look at your output netCDF file to confirm that the thickness and velocity fields are changing over time. Since 10 years is not nearly long enough to reach a new equilibrium, you should reset the value of tend in your ".config" file to a larger value; try 100 years with a time step of 2 years. This still won't be long enough to reach a new steady-state, but it is at least long enough for you to confirm that a new equilibrium is being approached.

A Higher-order Advection Scheme

As implied by the name "first-order upwinding", the scheme we have just implemented is only first-order accurate. In fact, first-order upwinding schemes are highly diffusive; they do not do a good job of preserving "sharp" features that are simply being advected by the flow field. A good example of this is shown in the figure Comparison of upwinding schemes. In this example, we are advecting a "block" of uniform horizontal and vertical dimensions in a uniform flow field oriented at 45 degrees, up and to the right (or to the northeast). The first-order scheme is clearly diffusive, relative to the higher-order [|incremental remapping scheme] it is compared with.

You can compare the results from evolving the geometry in your model with the first-order scheme versus the higher-order incremental remapping scheme by changing the evolution flag in your .config file from 4 to 3. Save your 100 year run from the previous experiment (using the upwinding scheme) and re-run the same experiment using incremental remapping to evolve the geometry. Do you notice any significant differences at the end of the model run?

Comparison of upwinding schemes:Time series (right-hand panels) of vertical profiles showing a block of uniform horizontal and vertical dimensions being advected through a uniform flow field oriented at 45 degrees, up and to the right. Color plots on the left-hand side show elevation contours for the block after starting in the lower-left hand corner and being advected through the flow field. The upwinding scheme, on the bottom, is highly diffusive while the higher-order incremental remapping scheme does a much better job of preserving the initial shape and dimensions of the block (click for higher-resolution image).