sico_variables.F90

Go to the documentation of this file.

!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Module :  s i c o _ v a r i a b l e s
!
!> @file
!!
!! Declarations of global variables for SICOPOLIS.
!!
!! @section Copyright
!!
!! Copyright 2009-2013 Ralf Greve
!!
!! @section License
!!
!! This file is part of SICOPOLIS.
!!
!! SICOPOLIS is free software: you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! SICOPOLIS 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
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with SICOPOLIS.  If not, see <http://www.gnu.org/licenses/>.
!<
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!-------------------------------------------------------------------------------
!> Declarations of global variables for SICOPOLIS.
!<------------------------------------------------------------------------------
module sico_variables

use sico_types

!-------- Field quantities --------

!> maske(j,i): Ice-land-ocean mask.
!>             0: grounded ice,
!>             1: ice-free land,
!>             2: ocean,
!>             3: floating ice
   integer(i2b), dimension(0:JMAX,0:IMAX) :: maske
!> maske_help(j,i): Present-day ice-land-ocean mask
   integer(i2b), dimension(0:JMAX,0:IMAX) :: maske_help
!> (.)_neu: New value of quantity (.) computed during an integration step
   integer(i2b), dimension(0:JMAX,0:IMAX) :: maske_neu
!> n_cts(j,i): Mask for thermal conditions.
!>             -1: cold ice base,
!>              0: temperate ice base with cold ice above,
!>              1: temperate ice base with temperate ice layer above
!>                 (only in polythermal mode)
   integer(i2b), dimension(0:JMAX,0:IMAX) :: n_cts
!> (.)_neu: New value of quantity (.) computed during an integration step
   integer(i2b), dimension(0:JMAX,0:IMAX) :: n_cts_neu
!> kc_cts(j,i): Position kc of the CTS (only in cold-ice mode)
   integer(i4b), dimension(0:JMAX,0:IMAX) :: kc_cts
!> flag_grounding_line(j,i): Grounding line flag.
!>                            .true.: grounding line point
!>                                    (grounded ice point with at least
!>                                    one floating ice neighbour),
!>                           .false.: otherwise
   logical, dimension(0:JMAX,0:IMAX) :: flag_grounding_line
!> flag_calving_front(j,i): Calving front flag.
!>                            .true.: calving front point
!>                                    (floating ice point with at least
!>                                    one ocean neighbour),
!>                           .false.: otherwise
   logical, dimension(0:JMAX,0:IMAX) :: flag_calving_front
!> flag_sf(j,i): Shelfy stream flag.
!>                            .true.: shelfy stream is applied
!>                           .false.: shelfy stream is not applied
   logical, dimension(0:JMAX,0:IMAX) :: flag_sf
!> xi(i): Coordinate xi (= x) of grid point i
   real(dp), dimension(0:IMAX) :: xi
!> eta(j): Coordinate eta (= y) of grid point j
   real(dp), dimension(0:JMAX) :: eta
!> zeta_c(kc): Sigma coordinate zeta_c of grid point kc
   real(dp), dimension(0:KCMAX) :: zeta_c
!> eaz_c(kc): Abbreviation for exp(DEFORM*zeta(kc))
   real(dp), dimension(0:KCMAX) :: eaz_c
!> eaz_c_quotient(kc): Abbreviation for (eaz_c(kc)-1.0)/(ea-1.0)
   real(dp), dimension(0:KCMAX) :: eaz_c_quotient
!> zeta_t(kt): Sigma coordinate zeta_t of grid point kt
   real(dp), dimension(0:KTMAX) :: zeta_t
!> ea: Abbreviation for exp(DEFORM)
   real(dp) :: ea

!> lambda(j,i): Geographic longitude of grid point (i,j)
   real(dp), dimension(0:JMAX,0:IMAX) :: lambda
!> phi(j,i): Geographic latitude of grid point (i,j)
   real(dp), dimension(0:JMAX,0:IMAX) :: phi
!> area(j,i): Area of grid cell associated with grid point (i,j)
   real(dp), dimension(0:JMAX,0:IMAX) :: area
!> sq_g11_g(j,i): Square root of the coefficient g11 of the metric tensor
!>                on grid point (i,j)
   real(dp), dimension(0:JMAX,0:IMAX) :: sq_g11_g
!> sq_g22_g(j,i): Square root of the coefficient g22 of the metric tensor
!>                on grid point (i,j)
   real(dp), dimension(0:JMAX,0:IMAX) :: sq_g22_g
!> insq_g11_g(j,i): Inverse square root of g11 on grid point (i,j)
   real(dp), dimension(0:JMAX,0:IMAX) :: insq_g11_g
!> insq_g22_g(j,i): Inverse square root of g22 on grid point (i,j)
   real(dp), dimension(0:JMAX,0:IMAX) :: insq_g22_g
!> sq_g11_sgx(j,i): Square root of g11, at (i+1/2,j)
   real(dp), dimension(0:JMAX,0:IMAX) :: sq_g11_sgx
!> sq_g11_sgy(j,i): Square root of g11, at (i,j+1/2)
   real(dp), dimension(0:JMAX,0:IMAX) :: sq_g11_sgy
!> sq_g22_sgx(j,i): Square root of g22, at (i+1/2,j)
   real(dp), dimension(0:JMAX,0:IMAX) :: sq_g22_sgx
!> sq_g22_sgy(j,i): Square root of g22, at (i,j+1/2)
   real(dp), dimension(0:JMAX,0:IMAX) :: sq_g22_sgy
!> insq_g11_sgx(j,i): Inverse square root of g11, at (i+1/2,j)
   real(dp), dimension(0:JMAX,0:IMAX) :: insq_g11_sgx
!> insq_g22_sgy(j,i): Inverse square root of g22, at (i,j+1/2)
   real(dp), dimension(0:JMAX,0:IMAX) :: insq_g22_sgy
!> zs(j,i): Coordinate z of the surface topography
   real(dp), dimension(0:JMAX,0:IMAX) :: zs
!> zm(j,i): Coordinate z of the CTS
   real(dp), dimension(0:JMAX,0:IMAX) :: zm
!> zb(j,i): Coordinate z of the ice base
   real(dp), dimension(0:JMAX,0:IMAX) :: zb
!> zl(j,i): Coordinate z of the lithosphere surface
   real(dp), dimension(0:JMAX,0:IMAX) :: zl
!> zl0(j,i): zl for isostatically relaxed ice-free conditions
   real(dp), dimension(0:JMAX,0:IMAX) :: zl0
!> wss(j,i): Isostatic steady-state displacement of the lithosphere
   real(dp), dimension(0:JMAX,0:IMAX) :: wss
!> H_c(j,i): Thickness of cold ice
   real(dp), dimension(0:JMAX,0:IMAX) :: H_c
!> H_t(j,i): Thickness of temperate ice
   real(dp), dimension(0:JMAX,0:IMAX) :: H_t
!> dzs_dxi(j,i): Derivative of zs by xi (at i+1/2,j)
   real(dp), dimension(0:JMAX,0:IMAX) :: dzs_dxi
!> dzm_dxi(j,i): Derivative of zm by xi (at i+1/2,j)
   real(dp), dimension(0:JMAX,0:IMAX) :: dzm_dxi
!> dzb_dxi(j,i): Derivative of zb by xi (at i+1/2,j)
   real(dp), dimension(0:JMAX,0:IMAX) :: dzb_dxi
!> dH_c_dxi(j,i): Derivative of H_c by xi (at i+1/2,j)
   real(dp), dimension(0:JMAX,0:IMAX) :: dH_c_dxi
!> dH_t_dxi(j,i): Derivative of H_t by xi (at i+1/2,j)
   real(dp), dimension(0:JMAX,0:IMAX) :: dH_t_dxi
!> dzs_deta(j,i): Derivative of zs by eta (at i,j+1/2)
   real(dp), dimension(0:JMAX,0:IMAX) :: dzs_deta
!> dzm_deta(j,i): Derivative of zm by eta (at i,j+1/2)
   real(dp), dimension(0:JMAX,0:IMAX) :: dzm_deta
!> dzb_deta(j,i): Derivative of zb by eta (at i,j+1/2)
   real(dp), dimension(0:JMAX,0:IMAX) :: dzb_deta
!> dH_c_deta(j,i): Derivative of H_c by eta (at i,j+1/2)
   real(dp), dimension(0:JMAX,0:IMAX) :: dH_c_deta
!> dH_t_deta(j,i): Derivative of H_t by eta (at i,j+1/2)
   real(dp), dimension(0:JMAX,0:IMAX) :: dH_t_deta
!> (.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
   real(dp), dimension(0:JMAX,0:IMAX) :: dzs_dxi_g
!> (.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
   real(dp), dimension(0:JMAX,0:IMAX) :: dzm_dxi_g
!> (.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
   real(dp), dimension(0:JMAX,0:IMAX) :: dzb_dxi_g
!> (.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
   real(dp), dimension(0:JMAX,0:IMAX) :: dH_c_dxi_g
!> (.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
   real(dp), dimension(0:JMAX,0:IMAX) :: dH_t_dxi_g
!> (.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
   real(dp), dimension(0:JMAX,0:IMAX) :: dzs_deta_g
!> (.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
   real(dp), dimension(0:JMAX,0:IMAX) :: dzm_deta_g
!> (.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
   real(dp), dimension(0:JMAX,0:IMAX) :: dzb_deta_g
!> (.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
   real(dp), dimension(0:JMAX,0:IMAX) :: dH_c_deta_g
!> (.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
   real(dp), dimension(0:JMAX,0:IMAX) :: dH_t_deta_g
!> dzs_dtau(j,i): Derivative of zs by tau (time)
   real(dp), dimension(0:JMAX,0:IMAX) :: dzs_dtau
!> dzm_dtau(j,i): Derivative of zm by tau (time)
   real(dp), dimension(0:JMAX,0:IMAX) :: dzm_dtau
!> dzb_dtau(j,i): Derivative of zb by tau (time)
   real(dp), dimension(0:JMAX,0:IMAX) :: dzb_dtau
!> dH_c_dtau(j,i): Derivative of zl by tau (time)
   real(dp), dimension(0:JMAX,0:IMAX) :: dzl_dtau
!> dH_c_dtau(j,i): Derivative of H_c by tau (time)
   real(dp), dimension(0:JMAX,0:IMAX) :: dH_c_dtau
!> dH_t_dtau(j,i): Derivative of H_t by tau (time)
   real(dp), dimension(0:JMAX,0:IMAX) :: dH_t_dtau
!> d_help_b(j,i): Auxiliary quantity for the computation of vx_b and vy_b
   real(dp), dimension(0:JMAX,0:IMAX) :: d_help_b
!> p_b_w(j,i): Basal water pressure
   real(dp), dimension(0:JMAX,0:IMAX) :: p_b_w
!> vx_b(j,i): Velocity in x-direction at the ice base, at (i+1/2,j)
   real(dp), dimension(0:JMAX,0:IMAX) :: vx_b
!> vy_b(j,i): Velocity in y-direction at the ice base, at (i,j+1/2)
   real(dp), dimension(0:JMAX,0:IMAX) :: vy_b
!> vx_m(j,i): Mean (depth-averaged) velocity in x-direction, at (i+1/2,j)
   real(dp), dimension(0:JMAX,0:IMAX) :: vx_m
!> vy_m(j,i): Mean (depth-averaged) velocity in y-direction, at (i,j+1/2)
   real(dp), dimension(0:JMAX,0:IMAX) :: vy_m
!> (.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
   real(dp), dimension(0:JMAX,0:IMAX) :: vx_b_g
!> (.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
   real(dp), dimension(0:JMAX,0:IMAX) :: vy_b_g
!> vz_b(j,i): Velocity in z-direction at the ice base, at (i,j)
   real(dp), dimension(0:JMAX,0:IMAX) :: vz_b
!> vx_s_g(j,i): Velocity in x-direction at the ice surface, at (i,j)
   real(dp), dimension(0:JMAX,0:IMAX) :: vx_s_g
!> vy_s_g(j,i): Velocity in x-direction at the ice surface, at (i,j)
   real(dp), dimension(0:JMAX,0:IMAX) :: vy_s_g
!> vz_s(j,i): Velocity in z-direction at the ice surface, at (i,j)
   real(dp), dimension(0:JMAX,0:IMAX) :: vz_s
!> flui_ave_sia(j,i): Depth-averaged fluidity of the SIA
   real(dp), dimension(0:JMAX,0:IMAX) :: flui_ave_sia
!> h_diff(j,i): Diffusivity of the SIA evolution equation of the ice surface
   real(dp), dimension(0:JMAX,0:IMAX) :: h_diff
!> qx(j,i): Volume flux in x-direction (depth-integrated vx, at (i+1/2,j))
   real(dp), dimension(0:JMAX,0:IMAX) :: qx
!> qy(j,i): Volume flux in y-direction (depth-integrated vy, at (i,j+1/2))
   real(dp), dimension(0:JMAX,0:IMAX) :: qy
!> q_gl_g(j,i): Volume flux across the grounding line, at (i,j)
   real(dp), dimension(0:JMAX,0:IMAX) :: q_gl_g
!> q_geo(j,i): Geothermal heat flux
   real(dp), dimension(0:JMAX,0:IMAX) :: q_geo
!> temp_b(j,i): Basal temperature
   real(dp), dimension(0:JMAX,0:IMAX) :: temp_b
!> temph_b(j,i): Basal temperature relative to the pressure melting point
   real(dp), dimension(0:JMAX,0:IMAX) :: temph_b
!> Q_bm(j,i): Basal melting rate
   real(dp), dimension(0:JMAX,0:IMAX) :: Q_bm
!> Q_tld(j,i): Water drainage rate from the temperate layer
   real(dp), dimension(0:JMAX,0:IMAX) :: Q_tld
!> Q_b_tot(j,i): Sum of Q_bm and Q_tld
   real(dp), dimension(0:JMAX,0:IMAX) :: Q_b_tot
!> H_w(j,i): Effective thickness of subglacial water
   real(dp), dimension(0:JMAX,0:IMAX) :: H_w
!> accum(j,i): Accumulation rate at the ice surface
!>             (includes liquid precipitation = rainfall!)
   real(dp), dimension(0:JMAX,0:IMAX) :: accum
!> evap(j,i): Evaporation rate at the ice surface
   real(dp), dimension(0:JMAX,0:IMAX) :: evap
!> runoff(j,i): Runoff rate at the ice surface
   real(dp), dimension(0:JMAX,0:IMAX) :: runoff
!> as_perp(j,i): Accumulation-ablation function at the ice surface
   real(dp), dimension(0:JMAX,0:IMAX) :: as_perp
!> temp_s(j,i): Ice surface temperature
   real(dp), dimension(0:JMAX,0:IMAX) :: temp_s
!> am_perp(j,i): Ice volume flux across the CTS
   real(dp), dimension(0:JMAX,0:IMAX) :: am_perp
!> am_perp_st(j,i): Steady-state part of am_perp
!>                  (without contribution of dzm_dtau)
   real(dp), dimension(0:JMAX,0:IMAX) :: am_perp_st
!> (.)_neu: New value of quantity (.) computed during an integration step
   real(dp), dimension(0:JMAX,0:IMAX) :: zs_neu
!> (.)_neu: New value of quantity (.) computed during an integration step
   real(dp), dimension(0:JMAX,0:IMAX) :: zm_neu
!> (.)_neu: New value of quantity (.) computed during an integration step
   real(dp), dimension(0:JMAX,0:IMAX) :: zb_neu
!> (.)_neu: New value of quantity (.) computed during an integration step
   real(dp), dimension(0:JMAX,0:IMAX) :: zl_neu
!> (.)_neu: New value of quantity (.) computed during an integration step
   real(dp), dimension(0:JMAX,0:IMAX) :: H_c_neu
!> (.)_neu: New value of quantity (.) computed during an integration step
   real(dp), dimension(0:JMAX,0:IMAX) :: H_t_neu

!> zs_ref(j,i): Reference elevation for precip_present, temp_ma_present
!>              and temp_mj_present
   real(dp), dimension(0:JMAX,0:IMAX) :: zs_ref

!> accum_present(j,i): Present-day accumulation rate at the ice surface
!>                     (for EISMINT, ISMIP HEINO and the north and south
!>                     polar caps of Mars)
   real(dp), dimension(0:JMAX,0:IMAX) :: accum_present
!> precip_ma_present(j,i): Present-day mean annual precipitation rate
!>                         at the ice surface
   real(dp), dimension(0:JMAX,0:IMAX) :: precip_ma_present
!> precip_ma_lgm_anom(j,i): LGM anomaly (ratio LGM/present) of the mean annual
!>                          precipitation rate at the ice surface
   real(dp), dimension(0:JMAX,0:IMAX) :: precip_ma_lgm_anom
!> temp_ma_present(j,i): Present-day mean annual surface temperature
   real(dp), dimension(0:JMAX,0:IMAX) :: temp_ma_present
!> temp_mj_present(j,i): Present-day mean summer (northern hemisphere: July,
!>                       southern hemisphere: January) surface temperature
   real(dp), dimension(0:JMAX,0:IMAX) :: temp_mj_present
!> temp_ma_lgm_anom(j,i): LGM anomaly (difference LGM - present) of the mean
!>                        annual surface temperature
   real(dp), dimension(0:JMAX,0:IMAX) :: temp_ma_lgm_anom
!> temp_mj_lgm_anom(j,i): LGM anomaly (difference LGM - present) of the mean
!>                        summer (northern hemisphere: July, southern
!>                        hemisphere: January) surface temperature
   real(dp), dimension(0:JMAX,0:IMAX) :: temp_mj_lgm_anom

!> dist_dxdy(jr,ir): Distance between grid points with delta_i=ir, delta_j=jr
   real(dp), dimension(-JMAX:JMAX,-IMAX:IMAX) :: dist_dxdy

!> precip_present(j,i,n): Present-day mean monthly precipitation rate
!>                        at the ice surface
   real(dp), dimension(0:JMAX,0:IMAX,12) :: precip_present
!> precip_lgm_anom(j,i,n): LGM anomaly (ratio LGM/present) of the mean monthly
!>                         precipitation rate at the ice surface
   real(dp), dimension(0:JMAX,0:IMAX,12) :: precip_lgm_anom
!> gamma_precip_lgm_anom(j,i,n): negative natural logarithm of
!>                               precip_lgm_anom(j,i,n)
   real(dp), dimension(0:JMAX,0:IMAX,12) :: gamma_precip_lgm_anom
!> temp_mm_present(j,i,n): Present-day mean monthly surface temperature
   real(dp), dimension(0:JMAX,0:IMAX,12) :: temp_mm_present
!> temp_mm_lgm_anom(j,i,n): LGM anomaly (difference LGM - present) of the mean
!>                          monthly surface temperature
   real(dp), dimension(0:JMAX,0:IMAX,12) :: temp_mm_lgm_anom

!> d_help_c(kc,j,i): Auxiliary quantity for the computation of vx, vy und zs
   real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: d_help_c
!> vx_c(kc,j,i): Velocity in x-direction in cold ice (at (i+1/2,j,kc))
   real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: vx_c
!> vy_c(kc,j,i): Velocity in y-direction in cold ice (at (i,j+1/2,kc))
   real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: vy_c
!> vz_c(kc,j,i): Velocity in z-direction in cold ice (at (i,j,kc+1/2))
   real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: vz_c
!> temp_c(kc,j,i): Temperature in cold ice
   real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: temp_c
!> (.)_neu: New value of quantity (.) computed during an integration step
   real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: temp_c_neu
!> temp_c_m(kc,j,i): Melting temperature in cold ice
   real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: temp_c_m
!> age_c(kc,j,i): Age in cold ice
   real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: age_c
!> (.)_neu: New value of quantity (.) computed during an integration step
   real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: age_c_neu
!> txz_c(kc,j,i): Shear stress txz in cold ice (at (i+1/2,j,kc))
   real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: txz_c
!> tyz_c(kc,j,i): Shear stress tyz in cold ice (at (i,j+1/2,kc))
   real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: tyz_c
!> sigma_c(kc,j,i): Effective stress in cold ice
   real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: sigma_c
!> enh_c(kc,j,i): Flow enhancement factor in cold ice
   real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: enh_c

!> de_ssa(j,i): Effective strain rate of the SSA, at (i,j)
   real(dp), dimension(0:JMAX,0:IMAX) :: de_ssa
!> vis_int_g(j,i): Depth-integrated viscosity of the SSA, at (i,j)
   real(dp), dimension(0:JMAX,0:IMAX) :: vis_int_g
!> vx_g(j,i): Velocity in x-direction of the SSA, at (i,j)
   real(dp), dimension(0:JMAX,0:IMAX) :: vx_g
!> vy_g(j,i): Velocity in y-direction of the SSA, at (i,j)
   real(dp), dimension(0:JMAX,0:IMAX) :: vy_g

!> d_help_t(kt,j,i): Auxiliary quantity for the computation of vx, vy und zs
   real(dp), dimension(0:KTMAX,0:JMAX,0:IMAX) :: d_help_t
!> vx_t(kt,j,i): Velocity in x-direction in temperate ice (at (i+1/2,j,kt))
   real(dp), dimension(0:KTMAX,0:JMAX,0:IMAX) :: vx_t
!> vy_t(kt,j,i): Velocity in y-direction in temperate ice (at (i,j+1/2,kt))
   real(dp), dimension(0:KTMAX,0:JMAX,0:IMAX) :: vy_t
!> vz_t(kt,j,i): Velocity in z-direction in temperate ice (at (i,j,kt+1/2))
   real(dp), dimension(0:KTMAX,0:JMAX,0:IMAX) :: vz_t
!> omega_t(kt,j,i): Water content in temperate ice
   real(dp), dimension(0:KTMAX,0:JMAX,0:IMAX) :: omega_t
!> (.)_neu: New value of quantity (.) computed during an integration step
   real(dp), dimension(0:KTMAX,0:JMAX,0:IMAX) :: omega_t_neu
!> temp_t_m(kt,j,i): Melting temperature in temperate ice (= actual temperature)
   real(dp), dimension(0:KTMAX,0:JMAX,0:IMAX) :: temp_t_m
!> age_t(kt,j,i): Age in temperate ice
   real(dp), dimension(0:KTMAX,0:JMAX,0:IMAX) :: age_t
!> (.)_neu: New value of quantity (.) computed during an integration step
   real(dp), dimension(0:KTMAX,0:JMAX,0:IMAX) :: age_t_neu
!> txz_t(kt,j,i): Shear stress txz in temperate ice (at (i+1/2,j,kt))
   real(dp), dimension(0:KTMAX,0:JMAX,0:IMAX) :: txz_t
!> tyz_t(kt,j,i): Shear stress tyz in temperate ice (at (i,j+1/2,kt))
   real(dp), dimension(0:KTMAX,0:JMAX,0:IMAX) :: tyz_t
!> sigma_t(kt,j,i): Effective stress in temperate ice
   real(dp), dimension(0:KTMAX,0:JMAX,0:IMAX) :: sigma_t
!> enh_t(kt,j,i): Flow enhancement factor in temperate ice
   real(dp), dimension(0:KTMAX,0:JMAX,0:IMAX) :: enh_t

!> temp_r(kr,j,i): Temperature in the bedrock
   real(dp), dimension(0:KRMAX,0:JMAX,0:IMAX) :: temp_r
!> (.)_neu: New value of quantity (.) computed during an integration step
   real(dp), dimension(0:KRMAX,0:JMAX,0:IMAX) :: temp_r_neu

!-------- Physical parameters --------

!> RHO: Density of ice
   real(dp) :: RHO
!> RHO_W: Density of pure water
   real(dp) :: RHO_W
!> RHO_SW: Density of sea water
   real(dp) :: RHO_SW
!> L: Latent heat of ice
   real(dp) :: L
!> G: Acceleration due to gravity
   real(dp) :: G
!> NUE: Water diffusivity in ice
   real(dp) :: NUE
!> BETA: Clausius-Clapeyron gradient of ice
   real(dp) :: BETA
!> DELTA_TM_SW: Melting point depression of sea water due to its
!>              average salinity
   real(dp) :: DELTA_TM_SW
!> OMEGA_MAX: Threshold value for the water content of temperate ice
   real(dp) :: OMEGA_MAX
!> H_R: Thickness of the modelled lithosphere layer
   real(dp) :: H_R
!> RHO_C_R: Density times specific heat of the lithosphere
   real(dp) :: RHO_C_R
!> KAPPA_R: Heat conductivity of the lithosphere
   real(dp) :: KAPPA_R
!> FLEX_RIG_L: Flexural rigidity of the lithosphere
   real(dp) :: FLEX_RIG_L
!> RHO_A: Density of the asthenosphere
   real(dp) :: RHO_A
!> TIME_LAG: Time lag of the relaxing asthenosphere
   real(dp) :: TIME_LAG
!> R_T: Coefficient of the water-content dependence in the rate factor
!>      for temperate ice
   real(dp) :: R_T
!> R: Radius of the planet
   real(dp) :: R
!> A: Semi-major axis of the planet
   real(dp) :: A
!> B: Semi-minor axis of the planet
   real(dp) :: B
!> LAMBDA0: Reference longitude (central meridian) of the stereographic
!>          projection
   real(dp) :: LAMBDA0
!> PHI0: Standard parallel of the stereographic projection
   real(dp) :: PHI0
!> S_STAT_0: Standard deviation of the air termperature for the
!>           degree-day model
   real(dp) :: S_STAT_0
!> BETA1_0: Degree-day factor for snow
   real(dp) :: BETA1_0
!> BETA1_LT_0: Degree-day factor for snow at low summer temperatures
   real(dp) :: BETA1_LT_0
!> BETA1_HT_0: Degree-day factor for snow at high summer temperatures
   real(dp) :: BETA1_HT_0
!> BETA2_0: Degree-day factor for ice
   real(dp) :: BETA2_0
!> BETA2_LT_0: Degree-day factor for ice at low summer temperatures
   real(dp) :: BETA2_LT_0
!> BETA2_HT_0: Degree-day factor for ice at high summer temperatures
   real(dp) :: BETA2_HT_0
!> PHI_SEP_0: Separation latitude for the computation of the degree-day
!>            factors beta1 and beta2: Equatorward of phi_sep, only the
!>            high-temperature values are used, whereas poleward of phi_sep,
!>            beta1 and beta2 are temperature-dependent
   real(dp) :: PHI_SEP_0
!> PMAX_0: Saturation factor for the formation of superimposed ice
   real(dp) :: PMAX_0
!> MU_0: Firn-warming correction
   real(dp) :: MU_0

#if ( defined(NMARS) \
      || defined(smars) )   /* martian ice sheet */
!> RHO_I: Density of ice
   real(dp) :: RHO_I
!> RHO_C: Density of crustal material (dust)
   real(dp) :: RHO_C
!> KAPPA_C: Heat conductivity of crustal material (dust)
   real(dp) :: KAPPA_C
!> C_C: Specific heat of crustal material (dust)
   real(dp) :: C_C
!> rho_inv: Inverse of the density of ice-dust mixture
   real(dp) :: rho_inv
#endif

!> RF(n): Tabulated values for the rate factor of cold ice
   real(dp), dimension(-190:10) :: RF
!> KAPPA(n): Tabulated values for the heat conductivity of ice
   real(dp), dimension(-190:10) :: KAPPA
!> C(n): Tabulated values for the specific heat of ice
   real(dp), dimension(-190:10) :: C

!-------- Further quantities -------- 

!> year_zero: SICOPOLIS year zero in astronomical year numbering
!>            [ = signed year CE (AD) ]
   real(dp) :: year_zero

#if ( defined(NMARS) \
      || defined(smars) )   /* martian ice sheet */
!> insol_time_min: Minimum time of the data values for the insolation etc.
   integer(i4b) :: insol_time_min
!> insol_time_stp: Time step of the data values for the insolation etc.
   integer(i4b) :: insol_time_stp
!> insol_time_max: Maximum time of the data values for the insolation etc.
   integer(i4b) :: insol_time_max
!> insol_ma_90(n): Data values for the mean-annual north- or south-polar
!>                 insolation
   real(dp), dimension(0:100000) :: insol_ma_90
!> obl_data(n): Data values for the obliquity
   real(dp), dimension(0:100000) :: obl_data
!> ecc_data(n): Data values for the eccentricity
   real(dp), dimension(0:100000) :: ecc_data
!> ave_data(n): Data values for the anomaly of vernal equinox
!>              (= 360 deg - Ls of perihelion )
   real(dp), dimension(0:100000) :: ave_data
!> cp_data(n): Data values for Laskar's climate parameter
!>             = eccentricity
!>               *sin(Laskar's longitude of perihelion from moving equinox),
!>               ( where Laskar's longitude of perihelion from moving equinox
!>                       = Ls of perihelion - 180 deg )
   real(dp), dimension(0:100000) :: cp_data
!> maske_chasm(j,i): Chasma mask.
!>                   0: grounded ice,
!>                   1: ice-free land (normal area),
!>                   7: chasma area
   integer(i2b), dimension(0:JMAX,0:IMAX) :: maske_chasm
!> time_chasm_init: Initial time for active chasma area
   real(dp) :: time_chasm_init
!> time_chasm_end: Final time for active chasma area
   real(dp) :: time_chasm_end
!> q_geo_normal(j,i): Geothermal heat flux for normal (non-chasma) areas
   real(dp), dimension(0:JMAX,0:IMAX) :: q_geo_normal
#endif

#if defined(EMTP2SGE)   /* EISMINT Phase 2 Simplified Geometry Experiment */
!> temp_min: Minimum surface temperature
   real(dp) :: temp_min
!> s_t: Gradient of surface temperature change with horizontal distance
   real(dp) :: s_t
!> x_hat: Coordinate xi (= x) of the centre of the model domain
   real(dp) :: x_hat
!> y_hat: Coordinate eta (= y) of the centre of the model domain
   real(dp) :: y_hat
!> b_max: Maximum accumulation rate
   real(dp) :: b_max
!> s_b: Gradient of accumulation rate change with horizontal distance
   real(dp) :: s_b
!> eld: Equilibrium line distance
   real(dp) :: eld
#elif defined(HEINO)    /* ISMIP HEINO */
!> maske_sedi(j,i): Sediment mask.
!>                  1: hard rock,
!>                  7: soft sediment,
!>                  2: ocean.
   integer(i2b), dimension(0:JMAX,0:IMAX) :: maske_sedi
!> temp_min: Minimum surface temperature
   real(dp) :: temp_min
!> s_t: Gradient of surface temperature change with horizontal distance
   real(dp) :: s_t
!> x_hat: Coordinate xi (= x) of the centre of the model domain
   real(dp) :: x_hat
!> y_hat: Coordinate eta (= y) of the centre of the model domain
   real(dp) :: y_hat
!> rad: Radius of the model domain
   real(dp) :: rad
!> b_min: Minimum accumulation rate
   real(dp) :: b_min
!> b_max: Maximum accumulation rate
   real(dp) :: b_max
#elif ( defined(GRL) && ICE_STREAM==2 )   /* Greenland with ice streams */
!> maske_sedi(j,i): Sediment mask.
!>                  1: hard rock,
!>                  7: soft sediment,
!>                  2: ocean.
   integer(i2b), dimension(0:JMAX,0:IMAX) :: maske_sedi
#elif ( defined(ANT) \
        && defined(sedi_slide) \
        && sedi_slide==2 )   /* antarctica with sediment sliding */
!> maske_sedi(j,i): Sediment mask.
!>                  1: hard rock,
!>                  7: soft sediment,
!>                  2: ocean.
   integer(i2b), dimension(0:JMAX,0:IMAX) :: maske_sedi
#endif

#if defined(ASF)  /* Austfonna */
!> n_surf: Number of surface points for which time-series data are written 
   integer(i4b) :: n_surf
!> lambda_surf(n): Geographical longitude of the prescribed surface points
   real(dp), dimension(:), allocatable, save :: lambda_surf
!> phi_surf(n): Geographical latitude of the prescribed surface points
   real(dp), dimension(:), allocatable, save :: phi_surf
!> x_surf(n): Coordinate xi (= x) of the prescribed surface points
   real(dp), dimension(:), allocatable, save :: x_surf
!> y_surf(n): Coordinate eta (= y) of the prescribed surface points
   real(dp), dimension(:), allocatable, save :: y_surf
!> snowfall(j,i): Snowfall rate at the ice surface
   real(dp), dimension(0:JMAX,0:IMAX) :: snowfall
!> rainfall(j,i): Rainfall rate at the ice surface
   real(dp), dimension(0:JMAX,0:IMAX) :: rainfall
#if ( ICE_STREAM==2 )   /* Austfonna with ice streams */
!> maske_sedi(j,i): Sediment mask.
!>                  1: hard rock,
!>                  7: soft sediment,
!>                  2: ocean.
   integer(i2b), dimension(0:JMAX,0:IMAX) :: maske_sedi
#endif
#endif

!> forcing_flag: Flag for the forcing type.
!>               1: forcing by a spatially constant surface temperature
!>                  anomaly (delta_ts),
!>               2: forcing by a glacial index (glac_index),
!>               3: forcing by time-dependent surface temperature
!>                  and precipitation data.
   integer(i2b) :: forcing_flag

!> n_core: Number of positions to be considered in the time-series file
!>         for deep boreholes
   integer(i4b) :: n_core
!> lambda_core(n): Geographical longitude of the prescribed borehole positions
   real(dp), dimension(:), allocatable, save :: lambda_core
!> phi_core(n): Geographical latitude of the prescribed borehole positions
   real(dp), dimension(:), allocatable, save :: phi_core
!> x_core(n): Coordinate xi (= x) of the prescribed borehole positions
   real(dp), dimension(:), allocatable, save :: x_core
!> y_core(n): Coordinate eta (= y) of the prescribed borehole positions
   real(dp), dimension(:), allocatable, save :: y_core

!> grip_time_min: Minimum time of the data values for the
!>                surface temperature anomaly
   integer(i4b) :: grip_time_min
!> grip_time_stp: Time step of the data values for the
!>                surface temperature anomaly
   integer(i4b) :: grip_time_stp
!> grip_time_max: Maximum time of the data values for the
!>                surface temperature anomaly
   integer(i4b) :: grip_time_max
!> ndata_grip: Number of data values for the surface temperature anomaly
   integer(i4b) :: ndata_grip
!> griptemp(n): Data values for the surface temperature anomaly
   real(dp), dimension(:), allocatable, save :: griptemp

!> gi_time_min: Minimum time of the data values for the glacial index
   integer(i4b) :: gi_time_min
!> gi_time_stp: Time step of the data values for the glacial index
   integer(i4b) :: gi_time_stp
!> gi_time_max: Maximum time of the data values for the glacial index
   integer(i4b) :: gi_time_max
!> ndata_gi: Number of data values for the glacial index
   integer(i4b) :: ndata_gi
!> glacial_index(n): Data values for the glacial index
   real(dp), dimension(:), allocatable, save :: glacial_index

!> specmap_time_min: Minimum time of the data values for the sea level
   integer(i4b) :: specmap_time_min
!> specmap_time_stp: Time step of the data values for the sea level
   integer(i4b) :: specmap_time_stp
!> specmap_time_max: Maximum time of the data values for the sea level
   integer(i4b) :: specmap_time_max
!> ndata_specmap: Number of data values for the sea level
   integer(i4b) :: ndata_specmap
!> specmap_zsl(n): Data values for the sea level
   real(dp), dimension(:), allocatable, save :: specmap_zsl

!> time_target_topo_init: Initial time for target-topography adjustment
   real(dp) :: time_target_topo_init
!> time_target_topo_final: Final time for target-topography adjustment
   real(dp) :: time_target_topo_final
!> maske_target(j,i): Target topography (ice-land-ocean mask)
   integer(i2b), dimension(0:JMAX,0:IMAX) :: maske_target
!> zs_target(j,i): Target topography (ice surface)
   real(dp), dimension(0:JMAX,0:IMAX) :: zs_target
!> zb_target(j,i): Target topography (ice base)
   real(dp), dimension(0:JMAX,0:IMAX) :: zb_target
!> zl_target(j,i): Target topography (lithosphere surface)
   real(dp), dimension(0:JMAX,0:IMAX) :: zl_target
!> H_target(j,i): Target topography (ice thickness)
   real(dp), dimension(0:JMAX,0:IMAX) :: H_target

!> maske_maxextent(j,i): Maximum ice extent mask.
!>                       0: not allowed to glaciate,
!>                       1: allowed to glaciate.
   integer(i2b), dimension(0:JMAX,0:IMAX) :: maske_maxextent

!> ncid_temp_precip: ID of the NetCDF file containing the
!>                   surface-temperature and precipitation data
!>                   as functions of time
   integer(i4b) :: ncid_temp_precip
!> ndata_temp_precip: Number of surface-temperature
!>                    and precipitation data
   integer(i4b) :: ndata_temp_precip
!> temp_precip_time_min: Minimum time of the surface-temperature
!>                       and precipitation data
   real(dp) :: temp_precip_time_min
!> temp_precip_time_stp: Time step of the surface-temperature
!>                       and precipitation data
   real(dp) :: temp_precip_time_stp
!> temp_precip_time_max: Maximum time of the surface-temperature
!>                       and precipitation data
   real(dp) :: temp_precip_time_max
!> temp_precip_time(n): Times of the surface-temperature
!>                      and precipitation data
   real(dp), dimension(:), allocatable, save :: temp_precip_time

!> kei(n): Tabulated values of the kei function (Kelvin function of zero order)
   real(dp), dimension(-10000:10000) :: kei
!> n_data_kei: Number of tabulated values of the kei function
   integer(i4b):: n_data_kei
!> kei_r_max: Maximum value of the argument r of the tabulated kei function
   real(dp) :: kei_r_max
!> kei_r_incr: Increment of the argument r of the tabulated kei function
   real(dp) :: kei_r_incr

!> pi: Mathematical constant
   real(dp), parameter :: pi = 3.141592653589793_dp
!> pi_inv: Inverse of pi
   real(dp), parameter :: pi_inv = 1.0_dp/pi
!> pi_180: pi divided by 180 (-> deg to rad)
   real(dp), parameter :: pi_180 = pi/180.0_dp
!> pi_180_inv: 180 divided by pi (-> rad to deg)
   real(dp), parameter :: pi_180_inv = 180.0_dp/pi
!> eps: Small number
   real(dp), parameter :: eps = 1.0e-05_dp

end module sico_variables
!