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.

## Step 1

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

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.

module fo_upwind_advect_mod

! subroutines for mass advection scheme based on 1st order upwinding

contains

!----------------------------------------------------------------------

! initialization for 1st-order upwinding mass advection

!----------------------------------------------------------------------

! finalization for 1st-order upwinding mass advection

!----------------------------------------------------------------------

! driver for 1st-order upwind mass advection

!----------------------------------------------------------------------

!----------------------------------------------------------------------

end module fo_upwind_advect_mod

Add this new module to the build by altering the file Makefile.am ...

libglide_a_SOURCES = glide.F90 glide_setup.F90 glide_types.F90 glide_temp.F90 \
glide_bwater.F90 glide_deriv.F90 xls.F90 ice3d_lib.F90 \
glide_velo_higher.F90 glide_thck.F90 glide_velo.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\
glam.F90 glam_strs2.F90 glam_thck_ppm.F90\
remap_advection.F90 remap_glamutils.F90 glide_ground.F90

becomes ...

libglide_a_SOURCES = glide.F90 glide_setup.F90 glide_types.F90 glide_temp.F90 \
glide_bwater.F90 glide_deriv.F90 xls.F90 ice3d_lib.F90 \
glide_velo_higher.F90 glide_thck.F90 glide_velo.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\
glam.F90 glam_strs2.F90 glam_thck_ppm.F90\
remap_advection.F90 remap_glamutils.F90 glide_ground.F90 \     ! CHANGED LAST TWO LINES
fo_upwind_advect.F90

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

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 207 for documentation ...

    !*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 475, change ...

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

Time to 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. 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).

In 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_mod, only : fo_upwind_advect_init         ! ADDED

Now add if construct and call to the initialization stub in fo_upwind_advect_mod.F90, around line 190 of glide.F90, just after the if construct for EVOL_INC_REMAP ...

    if (model%options%whichevol== EVOL_FO_UPWIND ) then
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_mod, 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
endif

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

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

... and add a call to the driver routine to the case construct (passing no args yet) ...

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

end select

Make sure to place it BEFORE the "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 w/ 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.

    subroutine fo_upwind_advect_main( thck, stagthck, acab, dt, uflx, vflx, ewn, nsn, dew, dns )  ! CHANGED

real (kind = dp), intent(in) :: dt					! ADDED THESE DEFINITIONS
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, Hu, Hd				! ADDED
integer :: ew, ns

end subroutine fo_upwind_advect_main

Because these definitions require kind "dp", we need to add a use statment at the start of the module (I've also added a few more we'll use later) ...

module fo_upwind_advect_mod

! subroutines for mass advection scheme based on 1st order upwinding

use glimmer_paramets, only: sp, dp, len0, thk0, tim0, vel0        ! ADDED

... etc ...

Then, pass these same arguements from the driver routine fo_upwind_advect_driver in fo_upwind_advect.F90. Note that this requires accessing the variables from the derived type "model", which must be added to the use statements, and also given a type definition ...

module fo_upwind_advect

! this is the driver routine for the 1st order, upwind mass transport scheme

use glide_types									! ADDED USE STATEMENT

implicit none
private

contains

type(glide_global_type), intent(inout) :: model		        ! ADDED DEFINITION

model%climate%acab,     model%numerics%dt,          &
model%velocity_hom%uflx,model%velocity_hom%vflx,    &
model%general%ewn,      model%general%nsn,          &
model%numerics%dew,     model%numerics%dns )

end module fo_upwind_advect

Finally, 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

end select

Type make again and check 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_mod.F90, which (in this case) deal only with work arrays that we need during the calculation stage. 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 the necessary variables at the start of the fo_upwind_advection_mod.F90. The beginning of the module becomes ...

module fo_upwind_advect_mod

! subroutines for mass advection scheme based on 1st order upwinding

use glimmer_paramets, only: sp, dp, len0, thk0, tim0, vel0

! allocatable work arrays								! ADDED
real (kind = dp), allocatable, dimension(:,:) ::    &
ubar_e, ubar_w,                  &
vbar_u, vbar_d,                  &
ubar_grid, vbar_grid,          &
flux_net, thck_grid

contains

... etc ...

Now change the initialization and finalization subroutines ...

   subroutine fo_upwind_advect_init( ewn, nsn )			! ADDED THESE ARGS

! initialization for 1st-order upwinding mass advection

integer, intent(in) :: ewn, nsn   ! horizontal grid dimensions	! ADDED

! allocate work arrays						! ADDED THESE
allocate( ubar_e(ewn,nsn) ); ubar_e = 0.0_dp
allocate( ubar_w(ewn,nsn) ); ubar_w = 0.0_dp
allocate( vbar_u(ewn,nsn) ); vbar_u = 0.0_dp
allocate( vbar_d(ewn,nsn) ); vbar_d = 0.0_dp
allocate( ubar_grid(ewn+1,nsn+1) ); ubar_grid = 0.0_dp
allocate( vbar_grid(ewn+1,nsn+1) ); vbar_grid = 0.0_dp
allocate( thck_grid(ewn+2,nsn+2) ); thck_grid = 0.0_dp
allocate( flux_net(ewn,nsn) ); flux_net = 0.0_dp

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

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

As with the initialization subroutine, we add deallocation statements for work arrays in the finalization subroutine ...

    subroutine fo_upwind_advect_final( )

! finalization for 1st-order upwinding mass advection

! deallocate work arrays				! ADDED
deallocate( ubar_e )
deallocate( ubar_w )
deallocate( vbar_u )
deallocate( vbar_d )
deallocate( ubar_grid )
deallocate( vbar_grid )
deallocate( thck_grid )
deallocate( flux_net )

end subroutine fo_upwind_advect_final

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

## 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 having you figure out how to do the 1st order upwinding on your own, we are going to give you an outline of the details and provide you with a chunk of code to plop in at the appropriate location. That code is here ...

! ----------------------------------

! code here

! ----------------------------------

... and it should be copy/pasted in fo_upwind_advect_main.F90 below the comment "guts of code starts here" and before the statement "end subroutine".