# Difference between revisions of "Adding dH/dt module"

Jump to: navigation, search

## Overview

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.

The general outline of operations in glide.F90, the main dynamics routine in Glimmer/CISM is as follows:

[1] initialize the model (open the door, turn on the lights, get all the toys out of the toybox)
[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 (put the toys away, turn off the lights, shut the door)


The new module that we want to add will contain initialization, "driver", and finalization subroutines (somewhat standard coding practice, at least for the Community Climate System Model). Thus we will have to 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; steps [2a] and [2c] are left alone). The subroutine calls will access these subroutines through appropriate "use" statements (more on that below).

To go directly to some exercises using the higher-order dynamics and the evolution scheme, go to Ice Sheet Evolution Experiments.

For mathematical background on upwinded transport, consult the Solving the equation for thickness evolution pages.

## Step 0

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:

./bootstrap
./configure --with-netcdf=/path/to_netcdf            !!! Note: specify your own path to netCDF libs here !!!
make


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

! 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, including calling code for the velocity calculation and the thickness evolution), the new subroutine we are creating that will do the evolution calculation (that is, "advect" mass around according to the velocity field), and a finalization routine (again, here mainly for deallocation of work arrays).

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

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_mask.F90 glide_stop.F90 glide_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 glissade.F90 glissade_remap.F90\
glissade_velo.F90 glissade_constants.F90 glide_vertint.F90\
glide_thckmask.F90 glide_nonlin.F90 glide_grids.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_mask.F90 glide_stop.F90 glide_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 glissade.F90 glissade_remap.F90\
glissade_velo.F90 glissade_constants.F90 glide_vertint.F90\
glide_thckmask.F90 glide_nonlin.F90 glide_grids.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

[Note: the comment in the 2nd to last line may cause a problem for make.]

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                      ', &
'ADI scheme                            ', &
'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, only : fo_upwind_advect_init         ! ADDED

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

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

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

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

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

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

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

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_hom%uflx,model%velocity_hom%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

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

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

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

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

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

The entire module containing the initialization, driver, and finalization subroutines can be found here:

## Ice Sheet Evolution

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: