glide.F90

! WJS (1-30-12): The following (turning optimization off) is needed as a workaround for an
! xlf compiler bug, at least in IBM XL Fortran for AIX, V12.1 on bluefire
#ifdef CPRIBM
@PROCESS OPT(0)
#endif

!CLEANUP - glide.F90
! Moved higher-order computations to a new module, glissade.F90.
! Simplified glide.F90 to include only SIA computations.
!
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!                                                             
!   glide.F90 - part of the Community Ice Sheet Model (CISM)  
!                                                              
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!   Copyright (C) 2005-2018
!   CISM contributors - see AUTHORS file for list of contributors
!
!   This file is part of CISM.
!
!   CISM is free software: you can redistribute it and/or modify it
!   under the terms of the Lesser GNU General Public License as published
!   by the Free Software Foundation, either version 3 of the License, or
!   (at your option) any later version.
!
!   CISM is distributed in the hope that it will be useful,
!   but WITHOUT ANY WARRANTY; without even the implied warranty of
!   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!   Lesser GNU General Public License for more details.
!
!   You should have received a copy of the Lesser GNU General Public License
!   along with CISM. If not, see <http://www.gnu.org/licenses/>.
!
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

#ifdef HAVE_CONFIG_H
#include "config.inc"
#endif

#include "glide_mask.inc"

!=======================================================================

module glide

  ! Driver for Glide (serial, SIA) dynamical core
  
  use glide_types
  use glide_stop
  use glide_io
  use glide_lithot
  use glide_profile
  use glimmer_config
  use glimmer_global, only: dp

  use glimmer_paramets, only: oldglide

  implicit none

  integer, private, parameter :: dummyunit=99

!WHL - debug
  logical, parameter :: verbose_glide = .false.

contains

!=======================================================================

  subroutine glide_config(model,config,fileunit)

    ! Read glide configuration from file and print it to the log

    use glide_setup
    use isostasy
    use glimmer_ncparams
    use glimmer_config
    use glimmer_map_init
    use glimmer_filenames

    implicit none

    type(glide_global_type), intent(inout) :: model  ! model instance
    type(ConfigSection), pointer  :: config          ! structure holding sections of configuration file
    integer, intent(in), optional :: fileunit        ! fileunit for reading config file 

    type(ConfigSection), pointer :: ncconfig
    integer :: unit

    unit = 99
    if (present(fileunit)) then
       unit = fileunit
    endif

    ! read configuration file
    call glide_readconfig (model,config)
    call glide_printconfig(model)

    ! read sigma levels from config file, if present
    call glide_read_sigma(model,config)

    !WHL - Moved isostasy configuration to glide_setup
!    call isos_readconfig(model%isos,config)
!    call isos_printconfig(model%isos)

    ! read mapping from config file
    ! **** Use of dew and dns here is an ugly fudge that
    ! **** allows the use of old [GLINT projection] config section
    ! **** for backwards compatibility. It will be deleted soon.
    ! **** (You have been warned!)
    ! **** N.B. Here, dew and dns are unscaled - i.e. real distances in m

    call glimmap_readconfig(model%projection,   config,   &
                            model%numerics%dew, model%numerics%dns)

    ! netCDF I/O
    if (trim(model%funits%ncfile) == '') then
       ncconfig => config
    else
       call ConfigRead(process_path(model%funits%ncfile), ncconfig, unit)
    end if

    call glimmer_nc_readparams(model, ncconfig)

  end subroutine glide_config

!=======================================================================

  subroutine glide_initialise(model)

    ! Initialise Glide model instance

    use glide_setup
    use glimmer_ncio
    use glide_velo, only: init_velo
    use glide_thck
    use glide_temp
    use glimmer_log
    use glimmer_scales
    use glide_mask
    use isostasy
    use glimmer_map_init
    use glimmer_coordinates, only: coordsystem_new
    use glide_diagnostics, only: glide_init_diag
    use glide_bwater
    use glimmer_paramets, only: len0

    use parallel, only: distributed_grid

    type(glide_global_type), intent(inout) :: model     ! model instance

!TODO - Is glimmer_version_char still used?
!       Old Glide does not include this variable.
    character(len=100), external :: glimmer_version_char
    character(len=100) :: message
    real(dp) :: smb_maxval

    integer, parameter :: my_nhalo = 0   ! no halo layers for Glide dycore

!!!Old Glide has this:
!!!    call write_log(glimmer_version)  

    call write_log(trim(glimmer_version_char()))

    ! initialise scales
    call glimmer_init_scales

    ! scale parameters (some conversions to SI units)
    call glide_scale_params(model)

    ! set up coordinate systems

    ! Note: nhalo = 0 is included in call to distributed_grid to set other halo
    !  variables (lhalo, uhalo, etc.) to 0 instead of default values

!WHL - distributed_grid is not in old glide
      
    call distributed_grid(model%general%ewn, model%general%nsn,  &
                          nhalo_in=my_nhalo)

    model%general%ice_grid = coordsystem_new(0.d0,               0.d0, &
                                             model%numerics%dew, model%numerics%dns, &
                                             model%general%ewn,  model%general%nsn)

    model%general%velo_grid = coordsystem_new(model%numerics%dew/2.d0, model%numerics%dns/2.d0, &
                                              model%numerics%dew,      model%numerics%dns,      &
                                              model%general%ewn-1,     model%general%nsn-1)

    ! allocate arrays
    call glide_allocarr(model)

!TODO - Eliminate the bed softness parameter and set btrc to model%velowo%btrac_const in glide_velo?
    ! initialise bed softness to uniform parameter
    model%velocity%bed_softness = model%velowk%btrac_const

    ! set uniform basal heat flux (positive down)
    !NOTE: This value will be overridden if we read bheatflx from an input file 
    !      (model%options%gthf = 1) or compute it (model%options%gthf = 2)
    model%temper%bheatflx = model%paramets%geot

    ! compute sigma levels or load from external file
    ! (if not already read from config file)
    call glide_load_sigma(model,dummyunit)

    ! initialize the time step counter
    ! For restart, tstep_count will be overwritten from the restart file.
    ! Alternatively, we could initialize tstep_count as follows:
    !    model%numerics%tstep_count = nint(model%numerics%time/model%numerics%tinc)
    ! But reading from a restart file should prevent roundoff issues.
    model%numerics%tstep_count = 0

    ! open all input files and forcing files
    call openall_in(model)

    ! read first time slice
    call glide_io_readall(model,model)

    ! Compute area scale factors for stereographic map projection.
    ! This should be done after reading the input file, in case the input file contains mapping info.
    ! Note: Not yet enabled for other map projections.
    ! TODO - Tested only for Greenland (N. Hem.; projection origin offset from N. Pole). Test for other grids.

    if (associated(model%projection%stere)) then

       call glimmap_stere_area_factor(model%projection%stere,  &
                                      model%general%ewn,       &
                                      model%general%nsn,       &
                                      model%numerics%dew*len0, &
                                      model%numerics%dns*len0)

    endif

    ! write projection info to log
    call glimmap_printproj(model%projection)
   
    ! If SMB input units are mm/yr w.e., then convert to units of acab.
    ! Note: In this case the input field should be called 'smb', not 'acab'.

    if (model%options%smb_input == SMB_INPUT_MMYR_WE) then

       ! make sure a nonzero SMB was read in
       smb_maxval = maxval(abs(model%climate%smb))
       if (smb_maxval < 1.0d-11) then
          write(message,*) 'Error: Failed to read in a nonzero SMB field with smb_input =', SMB_INPUT_MMYR_WE
          call write_log(trim(message), GM_FATAL)
       endif

       ! Convert units from mm/yr w.e. to m/yr ice
       model%climate%acab(:,:) = model%climate%smb(:,:) * (rhow/rhoi) / 1000.d0

       ! Convert acab from m/yr ice to model units
       model%climate%acab(:,:) = model%climate%acab(:,:) / scale_acab

    else
       ! assume acab was read in with units of m/yr ice; do nothing
    endif

    !WHL - Should have been read from glide_io_readall
    ! read lithot if required
!!    if (model%options%gthf > 0) then
!    if (model%options%gthf == GTHF_COMPUTE) then
!       call glide_lithot_io_readall(model,model)
!    end if

    ! handle relaxed/equilibrium topo
    ! Initialise isostasy first
    call init_isostasy(model)

    select case(model%isostasy%whichrelaxed)

    case(RELAXED_TOPO_INPUT)   ! Supplied input topography is relaxed

       model%isostasy%relx = model%geometry%topg

    case(RELAXED_TOPO_COMPUTE) ! Supplied topography is in equilibrium
                               !TODO - test case RELAXED_TOPO_COMPUTE
       call isos_relaxed(model)

    end select

    ! open all output files
    call openall_out(model)

    ! create glide variables
    call glide_io_createall(model, model)

!WHL - debug
!    print*, ' '
!    print*, 'Created Glide variables'
!    print*, 'max, min bheatflx (W/m2)=', maxval(model%temper%bheatflx), minval(model%temper%bheatflx)

    ! Compute the cell areas of the grid
    model%geometry%cell_area = model%numerics%dew*model%numerics%dns

    ! If a 2D bheatflx field is present in the input file, it will have been written 
    !  to model%temper%bheatflx.  For the case model%options%gthf = 0, we want to use
    !  a uniform heat flux instead.
    ! If no bheatflx field is present in the input file, then we default to the 
    !  prescribed uniform value, model%paramets%geot.

    if (model%options%gthf == GTHF_UNIFORM) then

       ! Check to see if this flux was present in the input file
       ! (by checking whether the flux is nonuniform over the domain)
       if (abs(maxval(model%temper%bheatflx) - minval(model%temper%bheatflx)) > 1.d-6) then  
          call write_log('Setting uniform prescribed geothermal flux')
          call write_log('(Set gthf = 1 to read geothermal flux field from input file)')
       endif

       ! set uniform basal heat flux (positive down)
       model%temper%bheatflx = model%paramets%geot

    endif
 
    ! Make sure the basal heat flux follows the positive-down sign convention                                                                                
    if (maxval(model%temper%bheatflx) > 0.0d0) then
       write(message,*) 'Error, Input basal heat flux has positive values: '
       call write_log(trim(message))
       write(message,*) 'maxval =', maxval(model%temper%bheatflx)
       call write_log(trim(message))
       write(message,*) 'Basal heat flux is defined as positive down, so should be <= 0 on input'
       call write_log(trim(message), GM_FATAL)
    endif

    !TODO - Change subroutine names to glide_init_velo, glide_init_thck

    ! initialise velocity calc
    call init_velo(model)

!WHL - old glide has a call to init_temp, which is similar to glide_init_temp
!      but does not set the temperature or compute flwa until later call to timeevoltemp

    ! Initialize temperature field - this needs to happen after input file is
    !  read so we can assign artm (which could possibly be read in) if temp has not been input.
    !
    ! Note: If the temperature field has not been read already from an input or restart file, 
    !        then temperature is initialized by this subroutine based on model%options%temp_init.  
    !       If the temperature has been read already, this subroutine will *not* overwrite it.
  
    call glide_init_temp(model)

    ! Initialize basal hydrology model, if enabled
    call bwater_init(model)

    ! initialise thickness evolution calc
    call init_thck(model)

    if (model%options%gthf == GTHF_COMPUTE) then
!!       call glide_lithot_io_createall(model)  !WHL - Variables should have been created by glide_io_createall
       call init_lithot(model)
    end if

!WHL - This call will set the ice column temperature to artm as in old glide,
!       regardless of the value of model%options%temp_init
!      Commented out at least for now.  To reproduce results of old_glide, make sure
!       model%options%temp_init = TEMP_INIT_ARTM.
!!  if (oldglide) then
!!    if (model%options%is_restart.ne.1) then
!!       ! initialise Glen's flow parameter A using an isothermal temperature distribution
!!       call glide_temp_driver(model,0)
!!    endif
!!  endif  ! oldglide
    
!WHL - This option is disabled for now.
    ! *mb* added; initialization of basal proc. module
!!    if (model%options%which_bproc == BAS_PROC_FULLCALC .or. &
!!        model%options%which_bproc == BAS_PROC_FASTCALC) then
!!        call Basal_Proc_init (model%general%ewn, model%general%nsn,model%basalproc,     &
!!                              model%numerics%dttem)
!!    end if      

    call glide_set_mask(model%numerics,                                   &
                        model%geometry%thck,  model%geometry%topg,     &
                        model%general%ewn,    model%general%nsn,       &
                        model%climate%eus,    model%geometry%thkmask,  &
                        model%geometry%iarea, model%geometry%ivol)
 
    ! calculate lower and upper ice surface
    call glide_calclsrf(model%geometry%thck, model%geometry%topg, model%climate%eus,model%geometry%lsrf)

    model%geometry%usrf = max(0.d0, model%geometry%thck + model%geometry%lsrf)

    ! initialise thckwk variables; used in timeders subroutine
    model%thckwk%olds(:,:,1) = model%geometry%thck(:,:)
    model%thckwk%olds(:,:,2) = model%geometry%usrf(:,:)

    ! initialise standard glide profiling
    call glide_prof_init(model)

    !TODO - Unclear on how subroutine register_model is used - Is it needed for serial code?
    ! register the newly created model so that it can be finalised in the case
    ! of an error without needing to pass the whole thing around to every
    ! function that might cause an error

    call register_model(model)

    ! initialise model diagnostics

    call glide_init_diag(model)

!WHL - debug
!    print*, 'After glide_initialise:'
!    print*, 'max, min thck (m)=', maxval(model%geometry%thck)*thk0, minval(model%geometry%thck)*thk0
!    print*, 'max, min usrf (m)=', maxval(model%geometry%usrf)*thk0, minval(model%geometry%usrf)*thk0
!    print*, 'max, min artm =', maxval(model%climate%artm), minval(model%climate%artm)
!    print*, 'max, min temp =', maxval(model%temper%temp), minval(model%temper%temp)
!    print*, 'max, min flwa =', maxval(model%temper%flwa), minval(model%temper%flwa)

!    print*, ' '
!    print*, 'thck:'
!    do j = model%general%nsn, 1, -1
!       write(6,'(30f5.0)') thk0 * model%geometry%thck(:,j)
!    enddo
!    print*, ' '
!    print*, 'temp, k = 2:'
!    do j = model%general%nsn+1, 0, -1
!       write(6,'(32f5.0)') model%temper%temp(2,:,j)
!    enddo
!    print*, 'basal temp:'
!    do j = model%general%nsn+1, 0, -1
!       write(6,'(32f5.0)') model%temper%temp(model%general%upn,:,j)
!    enddo

  end subroutine glide_initialise

!=======================================================================

  subroutine glide_init_state_diagnostic(model, evolve_ice)

    ! Calculate diagnostic variables for the initial model state
    ! This provides calculation of output fields at time 0
    ! This is analagous to glissade_diagnostic_variable_solve but is only 
    ! called from init.  The glide tstep routines take care of these calculations
    ! during time stepping.  
    ! Note that none of this is needed on a restart - this code ensures a complete 
    ! set of diagnostic output fields for the initial state.

    use glide_thck
    use glide_velo
    use glide_mask
    use glimmer_paramets, only: tim0
    use glimmer_physcon, only: scyr
    use glide_ground, only: glide_calve_ice
    use glide_bwater, only: calcbwat
    use glide_temp, only: glide_calcbmlt, glide_calcbpmp
    use glide_grid_operators

    type(glide_global_type), intent(inout) :: model     ! model instance
    logical, intent(in), optional :: evolve_ice         ! whether ice evolution is turned on (if not present, assumed true)

    integer :: i, j
    logical :: l_evolve_ice  ! local version of evolve_ice

    if (present(evolve_ice)) then
       l_evolve_ice = evolve_ice
    else
       l_evolve_ice = .true.
    end if

    if (model%options%is_restart == RESTART_TRUE) then
       ! On a restart, just assign the basal velocity from uvel/vvel (which are restart variables)
       ! to ubas/vbas which are used by the temperature solver to calculate basal heating.
       ! During time stepping ubas/vbas are calculated by slipvelo during thickness evolution or below on a cold start.
       model%velocity%ubas = model%velocity%uvel(model%general%upn,:,:)
       model%velocity%vbas = model%velocity%vvel(model%general%upn,:,:)

    else   ! Only make the calculations on a cold start

    ! ------------------------------------------------------------------------ 
    ! ***Part 1: Make geometry consistent with calving law, if desired
    ! ------------------------------------------------------------------------       

       if (l_evolve_ice .and. model%options%calving_init == CALVING_INIT_ON) then

          !WHL - debug
          print*, 'Calving at initialization, whichcalving =', model%options%whichcalving

          ! ------------------------------------------------------------------------ 
          ! Remove ice which should calve, depending on the value of whichcalving
          ! ------------------------------------------------------------------------ 
       
          !  Note: On a cold start, glide_calve_ice needs the mask to be calculated, but a call to 
          !  glide_set_mask occurs in glide_initialise and does not need to be repeated here.
          
          call glide_calve_ice(model%options%whichcalving,     &
                               model%geometry%thck,            &
                               model%isostasy%relx,            &
                               model%geometry%topg,            &
                               model%geometry%thkmask,         &
                               model%calving%marine_limit,     &
                               model%calving%calving_fraction, &
                               model%climate%eus,              &
                               model%calving%calving_thck)
          
          ! We now need to recalculate the mask because glide_calve_ice may have modified the geometry.
          call glide_set_mask(model%numerics,                                &
                              model%geometry%thck,  model%geometry%topg,     &
                              model%general%ewn,    model%general%nsn,       &
                              model%climate%eus,    model%geometry%thkmask,  &
                              model%geometry%iarea, model%geometry%ivol)
          
       endif  ! calving_init

       ! Compute total areas of grounded and floating ice
       call calc_iareaf_iareag(model%numerics%dew,    model%numerics%dns,     &
                               model%geometry%thkmask,                        &
                               model%geometry%iareaf, model%geometry%iareag)

    ! ------------------------------------------------------------------------ 
    ! ***Part 2: Calculate geometry related fields
    ! ------------------------------------------------------------------------    

    ! ------------------------------------------------------------------------
    ! calculate upper and lower ice surface
    ! ------------------------------------------------------------------------

       call glide_calclsrf(model%geometry%thck, model%geometry%topg, &
                        model%climate%eus,   model%geometry%lsrf)

       model%geometry%usrf = max(0.d0,model%geometry%thck + model%geometry%lsrf)

    ! ------------------------------------------------------------------------ 
    ! Calculate various derivatives
    !
    ! This call is needed here to make sure stagthck is calculated 
    ! the same way as in thck_lin_evolve/thck_nonlin_evolve
    ! ------------------------------------------------------------------------     

       call glide_prof_start(model,model%glide_prof%geomderv)

       call glide_geometry_derivs(model)   ! stagvarb, geomders as in old Glide

       call glide_prof_stop(model,model%glide_prof%geomderv)

       call glide_prof_start(model,model%glide_prof%ice_mask1)

    !TREY This sets local values of dom, mask, totpts, and empty
    !EIB! call veries between lanl and gc2, this is lanl version
    !magi a hack, someone explain what whichthck=5 does

!WHL - Modified this subroutine so that ice can accumulate in regions with
!      a small positive mass balance.

       call glide_thck_index(model%geometry% thck,      &
                          model%climate%  acab,      &
                          model%geometry% thck_index,  &
                          model%geometry% totpts,    &
                          .true.,                    &
                          model%geometry% empty)

       call glide_prof_stop(model,model%glide_prof%ice_mask1)


    ! ------------------------------------------------------------------------ 
    ! Part 3: Solve velocity
    ! ------------------------------------------------------------------------    

       ! initial value for flwa should already be calculated as part of glide_init_temp()
       ! calculate the part of the vertically averaged velocity field which solely depends on the temperature

        call velo_integrate_flwa(model%velowk,             &
                                 model%geomderv%stagthck,  &
                                 model%temper%flwa)

        ! Calculate diffusivity

       call velo_calc_diffu(model%velowk,              &
                            model%geomderv%stagthck,   &
                            model%geomderv%dusrfdew,   &
                            model%geomderv%dusrfdns,   &
                            model%velocity%diffu)

       ! If necessary, compute staggered variables required for basal traction calculation

       if (model%options%whichbtrc == BTRC_CONSTANT_BWAT) then

          !TODO - I think the next two calls are not needed, given that bwat should be in restart file for this option.

          ! Calculate basal melt rate --------------------------------------------------
          ! Note: For the initial state, we won't have values for ubas/vbas (unless they were 
          !       supplied in the input file) to get an initial guess of sliding heating.
          !       We could iterate on this, but for simplicity that is not done.

          call glide_calcbmlt(model, &
                              model%temper%temp, &
                              model%geometry%thck, &
                              model%geomderv%stagthck, &
                              model%geomderv%dusrfdew, &
                              model%geomderv%dusrfdns, &
                              model%velocity%ubas, &
                              model%velocity%vbas, &
                              model%basal_melt%bmlt, &
                              GLIDE_IS_FLOAT(model%geometry%thkmask))

          ! Note: calcbwat computes stagbwat
          call calcbwat(model, &
                        model%options%whichbwat, &
                        model%basal_melt%bmlt, &
                        model%temper%bwat, &
                        model%temper%bwatflx, &
                        model%geometry%thck, &
                        model%geometry%topg, &
                        model%temper%temp(model%general%upn,:,:), &
                        GLIDE_IS_FLOAT(model%geometry%thkmask), &
                        model%tempwk%wphi)


          ! This call is redundant for now, but is needed if the call to calcbwat is removed
          call stagvarb(model%temper%bwat, &
                        model%temper%stagbwat ,&
                        model%general%ewn, &
                        model%general%nsn)

       elseif (model%options%whichbtrc == BTRC_CONSTANT_BPMP) then

          call stagvarb(model%temper%temp(model%general%upn,1:model%general%ewn,1:model%general%nsn), &
                        model%temper%stagbtemp ,&
                        model%general%  ewn, &
                        model%general%  nsn)
       
          call glide_calcbpmp(model,               &
                              model%geometry%thck, &
                              model%temper%bpmp)
 
          call stagvarb(model%temper%bpmp, &
                        model%temper%stagbpmp ,&
                        model%general%  ewn, &
                        model%general%  nsn)

       endif   ! whichbtrc

       !------------------------------------------------------------------------ 
       ! Calculate basal traction factor
       !------------------------------------------------------------------------ 

       do j = 1, model%general%nsn-1
          do i = 1, model%general%ewn-1
             if (model%geomderv%stagthck(i,j)*thk0 < 1000.d0) then
                model%temper%stagbtemp(i,j) = model%temper%stagbpmp(i,j)
             else
                model%temper%stagbtemp(i,j) = -20.d0
             endif
          enddo
       enddo

       call calc_btrc(model,                    &
                      model%options%whichbtrc,  &
                      model%velocity%btrc)

       call slipvelo(model,                     &
                     0,                         &
                     model%velocity%btrc,       &
                     model%velocity%ubas,       &
                     model%velocity%vbas)
       

       ! Calculate velocity
       call velo_calc_velo(model%velowk,            model%geomderv%stagthck,  &
                           model%geomderv%dusrfdew, model%geomderv%dusrfdns,  &
                           model%temper%flwa,       model%velocity%diffu,     &
                           model%velocity%ubas,     model%velocity%vbas,      &
                           model%velocity%uvel,     model%velocity%vvel,      &
                           model%velocity%uflx,     model%velocity%vflx,      &
                           model%velocity%velnorm)    

    endif  ! if a restart


    ! MJH: I have left these calls outside of the restart if-construct so that there will
    ! always be a velnorm field calculated, which can be helpful for debugging.

    ! ------------------------------------------------------------------------ 
    ! Part 4: Calculate other diagnostic fields that depend on velocity
    ! ------------------------------------------------------------------------    

    ! ------------------------------------------------------------------------
    ! basal shear stress calculation
    ! ------------------------------------------------------------------------

    call calc_basal_shear(model%geomderv%stagthck,                          &
                          model%geomderv%dusrfdew, model%geomderv%dusrfdns, &
                          model%velocity%tau_x,    model%velocity%tau_y)

    ! velocity norm
    model%velocity%velnorm = sqrt(model%velocity%uvel**2 + model%velocity%vvel**2)

!WHL - debug
       if (verbose_glide) then

          print*, ' '
          print*, 'stagthck:'
          do i = 1, model%general%ewn-1
             write(6,'(i7)',advance='no') i
          enddo
          print*, ' '
          do j = model%general%nsn-1, 1, -1
             write(6,'(i3)',advance='no') j
             do i = 1, model%general%ewn-1
                write(6,'(f7.2)',advance='no') model%geomderv%stagthck(i,j)*thk0
             enddo 
             print*, ' '
          enddo

          print*, ' '
          print*, 'diffu (m^2/yr):'
          do i = 1, model%general%ewn-1
             write(6,'(i8)',advance='no') i
          enddo
          print*, ' '
          do j = model%general%nsn-1, 1, -1
             write(6,'(i3)',advance='no') j
             do i = 1, model%general%ewn-1
                write(6,'(f8.0)',advance='no') -model%velocity%diffu(i,j) * vel0*len0*scyr
             enddo
             print*, ' '
          enddo

          print*, ' '
          print*, 'ubas:'
          do i = 1, model%general%ewn-1
             write(6,'(i7)',advance='no') i
          enddo
          print*, ' '
          do j = model%general%nsn-1, 1, -1
             write(6,'(i4)',advance='no') j
             do i = 1, model%general%ewn-1
                write(6,'(f7.2)',advance='no') model%velocity%uvel(model%general%upn,i,j)*(vel0*scyr)
             enddo
             print*, ' '
          enddo

          print*, ' '
          print*, 'vbas:'
          do i = 1, model%general%ewn-1
             write(6,'(i7)',advance='no') i
          enddo
          print*, ' '
          do j = model%general%nsn-1, 1, -1
             write(6,'(i4)',advance='no') j
             do i = 1, model%general%ewn-1
                write(6,'(f7.2)',advance='no') model%velocity%vvel(model%general%upn,i,j)*(vel0*scyr)
             enddo
             print*, ' '
          enddo
          
          print*, ' '
          print*, 'uvel, k = 1:'
          do i = 1, model%general%ewn-1
             write(6,'(i8)',advance='no') i
          enddo
          print*, ' '
          do j = model%general%nsn-1, 1, -1
             write(6,'(i4)',advance='no') j
             do i = 1, model%general%ewn-1
                write(6,'(f8.2)',advance='no') model%velocity%uvel(1,i,j) * (vel0*scyr)
             enddo 
             print*, ' '
          enddo

          print*, ' '
          print*, 'u=vvel, k = 1:'
          do i = 1, model%general%ewn-1
             write(6,'(i8)',advance='no') i
          enddo
          print*, ' '
          do j = model%general%nsn-1, 1, -1
             write(6,'(i4)',advance='no') j
             do i = 1, model%general%ewn-1
                write(6,'(f8.2)',advance='no') model%velocity%vvel(1,i,j) * (vel0*scyr)
             enddo 
             print*, ' '
          enddo

       endif  ! verbose_glide

  end subroutine glide_init_state_diagnostic

!=======================================================================

  subroutine glide_tstep_p1(model,time)

    ! Perform first part of time-step of an ice model instance:
    ! temperature advection, vertical conduction, and internal dissipation.

    use glide_thck
    use glide_velo
    use glide_temp
    use glide_mask
    use glimmer_paramets, only: tim0
    use glimmer_physcon, only: scyr
    use glide_grid_operators
    use glide_bwater

    type(glide_global_type), intent(inout) :: model     ! model instance
    real(dp),  intent(in)   :: time                     ! current time in years

    ! Update internal clock
    model%numerics%time = time  
    model%numerics%tstep_count = model%numerics%tstep_count + 1
    model%temper%newtemps = .false.

    model%thckwk%oldtime = model%numerics%time - (model%numerics%dt * tim0/scyr)

    call glide_prof_start(model,model%glide_prof%geomderv)

    ! Update geometric quantities: stagthck, dusrfdew/dns, dthckdew/dns

    call glide_geometry_derivs(model)  ! compute stagthck, dusrfdew/dns, dthckdew/dns

    call glide_prof_stop(model,model%glide_prof%geomderv)

    call glide_prof_start(model,model%glide_prof%ice_mask1)

    !WHL - Modified this subroutine so that ice can accumulate in regions with
    !      a small positive mass balance.

    call glide_thck_index(model%geometry% thck,        &
                          model%climate%  acab,        &
                          model%geometry% thck_index,  &
                          model%geometry% totpts,      &
                          .true.,                      &
                          model%geometry% empty)

    call glide_prof_stop(model,model%glide_prof%ice_mask1)

    ! ------------------------------------------------------------------------ 
    ! calculate geothermal heat flux
    ! ------------------------------------------------------------------------ 

    if (model%options%gthf == GTHF_COMPUTE) then
       call calc_lithot(model)
    end if

    ! ------------------------------------------------------------------------ 
    ! Calculate temperature evolution and Glen's A, if necessary
    ! ------------------------------------------------------------------------ 

    ! Note: These times have units of years.
    !       dttem has scaled units, so multiply by tim0/scyr to convert to years

    if ( model%numerics%tinc >  mod(model%numerics%time,model%numerics%dttem*tim0/scyr)) then

       call glide_prof_start(model,model%glide_prof%temperature)

       if (oldglide) then   ! compute vertical velocity in glide_tstep_p1 
                          ! In new glide, this is called in glide_tstep_p3
         
          call glide_velo_vertical(model)

       endif   ! oldglide = T

       ! temperature advection, vertical conduction, and internal dissipation

       call glide_temp_driver(model, model%options%whichtemp)

       model%temper%newtemps = .true.

       call glide_prof_stop(model,model%glide_prof%temperature)

       ! Update hydrology, if needed ------------------------------------------------
       call calcbwat(model, &
                     model%options%whichbwat, &
                     model%basal_melt%bmlt, &
                     model%temper%bwat, &
                     model%temper%bwatflx, &
                     model%geometry%thck, &
                     model%geometry%topg, &
                     model%temper%temp(model%general%upn,:,:), &
                     GLIDE_IS_FLOAT(model%geometry%thkmask), &
                     model%tempwk%wphi)

    end if

    ! ------------------------------------------------------------------------ 
    ! Calculate basal traction factor
    ! ------------------------------------------------------------------------ 

    call calc_btrc(model,                    &
                   model%options%whichbtrc,  &
                   model%velocity%btrc)

!WHL - debug
!    print*, ' '
!    print*, 'After glide_tstep_p1:'
!    print*, 'max, min temp =', maxval(model%temper%temp), minval(model%temper%temp)
!    print*, 'max, min flwa =', maxval(model%temper%flwa), minval(model%temper%flwa)

!    print*, ' '
!    print*, 'temp, k = 2:'
!    do j = model%general%nsn+1, 0, -1
!       write(6,'(14f12.7)') model%temper%temp(2,3:16,j)
!    enddo
!    print*, 'basal temp:'
!    do j = model%general%nsn+1, 0, -1
!       write(6,'(14f12.7)') model%temper%temp(model%general%upn,3:16,j)
!    enddo

  end subroutine glide_tstep_p1

!=======================================================================

  subroutine glide_tstep_p2(model)

    ! Perform second part of time-step of an ice model instance:
    ! thickness evolution by one of several methods.

    use glide_thck
    use glide_velo
    use glide_temp
    use glide_mask
    use isostasy
    use glide_ground, only: glide_calve_ice

    type(glide_global_type), intent(inout) :: model    ! model instance

    !debug
    integer :: j

    ! ------------------------------------------------------------------------ 
    ! Calculate flow evolution by various different methods
    ! ------------------------------------------------------------------------ 

    call glide_prof_start(model,model%glide_prof%ice_evo)

    select case(model%options%whichevol)

    case(EVOL_PSEUDO_DIFF) ! Use precalculated uflx, vflx -----------------------------------

       call thck_lin_evolve(model,model%temper%newtemps)

    case(EVOL_ADI) ! Use explicit leap frog method with uflx,vflx -------------------

       call stagleapthck(model,model%temper%newtemps)

    case(EVOL_DIFFUSION) ! Use non-linear calculation that incorporates velocity calc -----

       call thck_nonlin_evolve(model,model%temper%newtemps)

    end select

    call glide_prof_stop(model,model%glide_prof%ice_evo)

    ! ------------------------------------------------------------------------ 
    ! get new mask
    ! Note: A call to glide_set_mask is needed before glide_calve_ice.
    ! ------------------------------------------------------------------------ 

    call glide_prof_start(model,model%glide_prof%ice_mask2)

    !TODO - Calculate area and vol separately from glide_set_mask?

    call glide_set_mask(model%numerics,                                &
                        model%geometry%thck,  model%geometry%topg,     &
                        model%general%ewn,    model%general%nsn,       &
                        model%climate%eus,    model%geometry%thkmask,  &
                        model%geometry%iarea, model%geometry%ivol)

    call glide_prof_stop(model,model%glide_prof%ice_mask2)

    ! ------------------------------------------------------------------------ 
    ! Remove ice which is either floating, or is present below prescribed
    ! depth, depending on value of whichcalving
    ! ------------------------------------------------------------------------ 

    !TODO - Some arguments for glide_calve_ice may not be needed.
    !       Old glide includes only arguments through model%climate%calving.

    call glide_calve_ice(model%options%whichcalving,     &
                         model%geometry%thck,            &
                         model%isostasy%relx,            &
                         model%geometry%topg,            &
                         model%geometry%thkmask,         &
                         model%calving%marine_limit,     &
                         model%calving%calving_fraction, &
                         model%climate%eus,              &
                         model%calving%calving_thck)

    ! Recalculate the mask following calving
    ! Note - This call to glide_set_mask is not in old Glide, but should have been.

 if (.not. oldglide) then   ! recalculate the thickness mask after calving
    call glide_set_mask(model%numerics,                                &
                        model%geometry%thck,  model%geometry%topg,     &
                        model%general%ewn,    model%general%nsn,       &
                        model%climate%eus,    model%geometry%thkmask,  &
                        model%geometry%iarea, model%geometry%ivol)
 endif   ! oldglide = F

 if (.not. oldglide) then   ! calculate area of floating and grounded ice
    call calc_iareaf_iareag(model%numerics%dew,    model%numerics%dns,     &
                            model%geometry%thkmask,                        &
                            model%geometry%iareaf, model%geometry%iareag)
 endif   ! oldglide = F

    ! ------------------------------------------------------------------------
    ! update ice/water load if necessary
    ! Note: Suppose the update period is 100 years, and the time step is 1 year.
    !       Then the update will be done on the first time step of the simulation,
    !        (model%numerics%tstep_count = 1) and again on step 101, 201, etc.
    !       The update will not be done before writing output at t = 100, when
    !        model%numerics%tstep_count = 100.
    !       Thus the output file will contain the load that was applied during the
    !        preceding years, not the new load.
    !       In older code versions, the new load would have been computed on step 100.
    ! ------------------------------------------------------------------------

    call glide_prof_start(model,model%glide_prof%isos_water)

    if (model%options%isostasy == ISOSTASY_COMPUTE) then

       if (model%isostasy%nlith > 0) then
          if (mod(model%numerics%tstep_count-1, model%isostasy%nlith) == 0) then
             call isos_icewaterload(model)
             model%isostasy%new_load = .true.
          end if
       endif  ! nlith > 0

    end if

    call glide_prof_stop(model,model%glide_prof%isos_water)
    
    ! ------------------------------------------------------------------------
    ! basal shear stress calculation
    ! ------------------------------------------------------------------------

! Old glide just passes 'model'

    call calc_basal_shear(model%geomderv%stagthck,                          &
                          model%geomderv%dusrfdew, model%geomderv%dusrfdns, &
                          model%velocity%tau_x,    model%velocity%tau_y)

! not in old glide, but this is a useful diagnostic

    ! velocity norm
    model%velocity%velnorm = sqrt(model%velocity%uvel**2 + model%velocity%vvel**2)

    !WHL - debug
!    print*, ' '
!    print*, 'After tstep_p2:'
!    print*, 'max, min thck (m)=', maxval(model%geometry%thck)*thk0, minval(model%geometry%thck)*thk0
!    print*, 'max, min usrf (m)=', maxval(model%geometry%usrf)*thk0, minval(model%geometry%usrf)*thk0
!    print*, 'max uvel, vvel =', maxval(model%velocity%uvel), maxval(model%velocity%vvel)

!    print*, ' '
!    print*, 'thck:'
!    do j = model%general%nsn, 1, -1
!       write(6,'(14f12.7)') thk0 * model%geometry%thck(3:16,j)
!!       write(6,'(14e12.5)') thk0 * model%geometry%thck(3:16,j)
!    enddo
!    print*, ' '
!    print*, 'btemp:'
!    do j = model%general%nsn, 1, -1
!       write(6,'(14f12.7)') model%temper%btemp(3:16,j)
!    enddo
!    print*, 'sfc uvel:'
!    do j = model%general%nsn-1, 1, -1
!       write(6,'(14f12.7)') model%velocity%uvel(1,3:16,j)
!    enddo
!    print*, 'sfc vvel:'
!    do j = model%general%nsn-1, 1, -1
!       write(6,'(14f12.7)') model%velocity%vvel(1,3:16,j)
!    enddo

  end subroutine glide_tstep_p2

!=======================================================================

  subroutine glide_tstep_p3(model)

    ! Perform third part of time-step of an ice model instance:
    ! calculate isostatic adjustment and upper and lower ice surface

    use isostasy
    use glimmer_scales, only: scale_acab
    use glide_setup
    use glide_velo, only: glide_velo_vertical
    use glide_thck, only: glide_calclsrf
    implicit none

    type(glide_global_type), intent(inout) :: model     ! model instance
    
    ! ------------------------------------------------------------------------ 
    ! Calculate isostasy
    ! ------------------------------------------------------------------------ 

    call glide_prof_start(model,model%glide_prof%isos)

    if (model%options%isostasy == ISOSTASY_COMPUTE) then
       call isos_compute(model)
    end if

    call glide_prof_stop(model,model%glide_prof%isos)

    ! ------------------------------------------------------------------------
    ! calculate upper and lower ice surface
    ! ------------------------------------------------------------------------

    call glide_calclsrf(model%geometry%thck, model%geometry%topg, &
                        model%climate%eus,   model%geometry%lsrf)

    model%geometry%usrf = max(0.d0,model%geometry%thck + model%geometry%lsrf)

    ! surface mass balance in units of mm/yr w.e.
    ! (model%climate%acab * scale_acab) has units of m/yr of ice
    ! Note: This is not necessary (and can destroy exact restart) if the SMB was already input in units of mm/yr
    if (model%options%smb_input /= SMB_INPUT_MMYR_WE) then
       model%climate%smb(:,:) = (model%climate%acab(:,:) * scale_acab) * (1000.d0 * rhoi/rhow)
    endif

    !Note: The time step counter used to be updated here; now it is updated at the start
    !of glide_tstep_p1.

    !TODO - Combine these timeders and vert velo calls into a subroutine?

    ! For exact restart, compute wgrd here and write it to the restart file.
    ! (This is easier than writing thckwk quantities to the restart file.)

 if (.not. oldglide) then  ! compute vertical velocity in glide_tstep_p3

    ! compute vertical velocity
         
    call t_startf('vertical_velo')

    call glide_velo_vertical(model)

    call t_stopf('vertical_velo')

 endif  ! oldglide = F

 !WHL - Moved netCDF output to simple_glide
 !!       call glide_io_writeall(model,model)

  end subroutine glide_tstep_p3

!=======================================================================

end module glide