sicopolis.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Program :  s i c o p o l i s . F 9 0
!             (SImulation COde for POLythermal Ice Sheets)
#define       VERSION '3.1'
#define       DATE    '2013-09-10'
!
!> @mainpage
!!
!! @section Description
!!
!! SICOPOLIS (SImulation COde for POLythermal Ice Sheets) is a 3-d
!! dynamic/thermodynamic model which simulates the evolution of large ice
!! sheets. It was originally created as a part of the doctoral thesis by
!! Greve (1995) in a version for the Greenland ice sheet. Since then,
!! SICOPOLIS has been developed continuously and applied to problems of
!! past, present and future glaciation of Greenland, Antarctica, the
!! entire northern hemisphere and also the polar ice caps of the planet
!! Mars.
!!
!! The model is based on the shallow ice approximation for grounded ice
!! and the shallow shelf approximation for floating ice (e.g., Greve and
!! Blatter 2009). It is coded in Fortran 90 and uses finite-difference
!! discretization on a staggered grid, the velocity components being
!! taken between grid points. Its particularity is the detailed treatment
!! of basal temperate layers (that is, regions with a temperature at the
!! pressure melting point), which are positioned by fulfilling a
!! Stefan-type jump condition at the interface to the cold ice
!! regions. Within the temperate layers, the water content is computed,
!! and its influence on the ice viscosity is taken into account.
!!
!! Required model forcing:
!! @li Surface mass balance (precipitation, evaporation, runoff).
!! @li Mean annual air temperature above the ice.
!! @li Eustatic sea level.
!! @li Geothermal heat flux.
!!      
!! Output (as functions of position and time):
!! @li Extent and thickness of the ice sheet.
!! @li Velocity field.
!! @li Temperature field.
!! @li Water-content field (temperate regions).
!! @li Age of the ice.
!! @li Isostatic displacement and temperature of the lithosphere.\n
!!      
!! References:
!! @li Greve, R. 1995.\n
!!     Thermomechanisches Verhalten polythermer Eisschilde - Theorie, Analytik,
!!     Numerik.\n
!!     Doctoral thesis, Department of Mechanics, Darmstadt University of
!!     Technology, Germany, 226 pp.\n
!!     ISBN: 3-8265-0999-4.
!! @li Greve, R. and H. Blatter. 2009.\n
!!     Dynamics of Ice Sheets and Glaciers.\n
!!     Springer, Berlin, Germany etc., 287 pp. ISBN: 978-3-642-03414-5.
!! @li SICOPOLIS website: http://sicopolis.greveweb.net/
!!
!! @section Copyright
!!
!! Copyright 2009-2013 Ralf Greve\n
!! (with contributions by Reinhard Calov, Thorben Dunse, Thomas Goelles,
!! Philipp Hancke, Nina Kirchner, Sascha Knell, Alex Robinson, Tatsuru Sato,
!! Malte Thoma)
!!
!! @section License
!!
!! 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/>.
!!
!! @file
!!
!! Main program file of 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/>.
!!
!<
!              *****************************************************************
!              Version history:
!              *****************************************************************
!              V1.0: Modular version of SICOPOLIS, applicable to any
!                    ice sheet.
!                    -----------------------------------------------------------
!                    Based on sicopolis_nhem.F95 V4.0
!                    (northern-hemisphere version of SICOPOLIS).
!              -----------------------------------------------------------------
!              V1.1: Optional new boundary conditions for the ice
!                    surface at the margin of the computational domain
!                    (parameter BC_MARG).
!              -----------------------------------------------------------------
!              V1.2: Optional first-order-upstream discretization
!                    of vertical advection terms without artificial
!                    diffusion in the age equation
!                    (parameter ADV_VERT).
!                    -----------------------------------------------------------
!                    Optional exponential coupling between accumulation
!                    and delta_ts (parameter ACCSURFACE==3).
!              -----------------------------------------------------------------
!              V1.3: New module for the southern hemisphere included
!                    (parameter DOMAIN==4).
!                    -----------------------------------------------------------
!                    Improved reading of ice-core-temperature and
!                    sea-level data, no longer restricted to the
!                    interval from 250 kyr BP to today.
!                    -----------------------------------------------------------
!                    New model domain for Scandinavia included
!                    (parameter DOMAIN==2).
!              -----------------------------------------------------------------
!              V1.4: New module for Antarctica included.
!                    -----------------------------------------------------------
!                    Bugfix in subroutines sico_init
!                    (for grl, shem, nmars):
!                    Parameter list in calls of subroutine boundary
!                    for parameter ANF_DAT==3 corrected.
!                    -----------------------------------------------------------
!                    Mars module: Degree-day model replaced by simpler
!                    accumulation-ablation parameterization.
!              -----------------------------------------------------------------
!              V2.0: All file extensions changed from '.F95' to '.F90'.
!                    -----------------------------------------------------------
!                    COMMON blocks in input file 'sico.common' replaced
!                    by module 'sico_variables'.
!                    -----------------------------------------------------------
!                    Character variable 'runname' (name of current run)
!                    can have an arbitrary length of less than
!                    100 characters.
!                    -----------------------------------------------------------
!                    Variables d_help_c, d_help_t, h_diff,
!                    maske, maske_help, n_cts, n_cts_neu, kc_cts,
!                    as_perp, temp_s
!                    now declared globally in module 'sico_variables'.
!                    -----------------------------------------------------------
!                    The header file 'sicoheader.F95' does no longer
!                    exist; however, 'sicopolis.F90' is compiled
!                    directly.
!                    -----------------------------------------------------------
!                    Header file with simulation specifications must
!                    be called 'sico_specs.h'.
!                    -----------------------------------------------------------
!                    Explicit definition of paths INPATH (for input
!                    files), OUTPATH (for output files), ANFDATPATH
!                    (for initial-value file) and Z_SLE_PATH (for
!                    previously computed time-series file with
!                    sea-level equivalents) in specifications file.
!                    -----------------------------------------------------------
!                    Parameter DOMAIN removed. Instead, the
!                    computational domain is now defined by a macro
!                    ANT, GRL, SCAND, NHEM, SHEM, EMTP2SGE or NMARS.
!                    -----------------------------------------------------------
!                    Parameter BC_MARG (see version 1.1) removed, back
!                    to zero-thickness boundary condition at the margin
!                    of the numerical domain as the only possibility.
!                    -----------------------------------------------------------
!                    Geothermal heat flux (Q_GEO), flow-enhancement factor
!                    (ENH) and minimum ice thickness interpreted as glaciation
!                    (H_MIN) now defined in header file sico_specs.h.
!                    -----------------------------------------------------------
!                    Physical parameters now read from a file by the
!                    new subroutine 'phys_para'.
!                    -----------------------------------------------------------
!                    Indexing of arrays RF(T), KAPPA(T) and C(T)
!                    changed such that the index T now corresponds to
!                    the temperature in deg C.
!                    -----------------------------------------------------------
!                    Computation of the horizontal volume flux qx, qy
!                    moved from subroutine 'calc_top' to subroutine
!                    'calc_vxy'.
!                    -----------------------------------------------------------
!                    Forcing quantities delta_ts and z_sl added to
!                    the output for the time-slice files produced by
!                    subroutine 'output1'.
!                    -----------------------------------------------------------
!                    Bugfix in subroutine 'calc_top':
!                    In the diffusive smoothing scheme for zs_neu, the
!                    denominator of the correction factor 'korrfakt'
!                    can no longer become zero due to the insertion
!                    of a volume offset eps_vol.
!              -----------------------------------------------------------------
!              V2.0.1: Maximum surface velocity now written to the time-series
!                      files.
!              -----------------------------------------------------------------
!              V2.0.2: Module for the southern hemisphere deleted.
!                      ---------------------------------------------------------
!                      New Greenland topography data based on Letreguilly and
!                      ETOPO5 with slightly smoothed relaxed bedrock used.
!              -----------------------------------------------------------------
!              V2.0.3: Limitation of the computed horizontal velocities to an
!                      interval [-VH_MAX, VH_MAX] to be prescribed in the
!                      header file.
!              -----------------------------------------------------------------
!              V2.0.4: Implementation of alternative flow laws
!                      (Goldsby-Kohlstedt, Durham) in addition to the usual
!                      Glen's flow law (parameter FLOW_LAW).
!                      ---------------------------------------------------------
!                      For the north polar cap of Mars, the average dust content
!                      of the ice can be specified (parameter FRAC_DUST).
!              -----------------------------------------------------------------
!              V2.1: Spatially variable geothermal heat flux implemented, can
!                    be chosen by new parameter Q_GEO_MOD.
!                    -----------------------------------------------------------
!                    Output of time-series file for ice-core locations
!                    prescribed in subroutine sico_init implemented
!                    (parameter OUTSER==3).
!                    -----------------------------------------------------------
!                    New module for the south-polar cap of Mars included
!                    (macro SMARS).
!                    -----------------------------------------------------------
!                    For the Martian ice caps (NMARS, SMARS), the mass balance
!                    is now converted from water equiv. to (ice+dust) equiv.
!                    (previously the dust was ignored).
!              -----------------------------------------------------------------
!              V2.2: Sub-melt basal sliding implemented (parameter GAMMA_SLIDE
!                    in phys_para file; 0: no sub-melt sliding).
!                    -----------------------------------------------------------
!                    Greenland:
!                    - New Bamber topography implemented, horizontal resolution
!                      40, 20 or 10 km.
!                    - Surface temperatures (mean annual, mean July)
!                      parameterised by Ritz et al. (1997).
!                    - Precipitation instead of solid accumulation used as
!                      input.
!                    - Reeh's degree day model with statistical fluctuations,
!                      supplemented by explicit consideration of rainfall
!                      and different ice/snow degree-day factors for warm
!                      and cold conditions (parameters BETA1_LT_0, BETA1_HT_0,
!                      BETA2_LT_0, BETA2_HT_0 in phys_para file).
!                    - Net surface mass balance computed as precipitation minus
!                      runoff (previously solid accumulation minus melting).
!                    - Optional climate forcing by glacial index and prescribed
!                      LGM anomalies for surface temperature and precipitation
!                      (parameters TSURFACE==5 and ACCSURFACE==5).
!                    -----------------------------------------------------------
!                    Martian ice caps (NMARS, SMARS):
!                    - Accumulation input either in water equiv.
!                      or in (ice+dust) equiv. (parameter ACC_UNIT).
!                    - Lower limit of -125/-128 deg C (sublimation temperature
!                      of CO2 over the NPC/SPC) introduced for the ice surface
!                      temperature.
!              -----------------------------------------------------------------
!              V2.3: Implementation of Bjoern Grieger's surface-temperature
!                    parameterisation for the Martian ice caps (parameter
!                    TSURFACE==6).
!                    -----------------------------------------------------------
!                    Greenland: Optional linear or exponential interpolation
!                    function for glacial-index-based precipitation forcing
!                    (parameter PRECIP_ANOM_INTERPOL).
!              -----------------------------------------------------------------
!              V2.4: Implementation of an elastic lithosphere (ELRA model),
!                    based on Sascha Knell's diploma thesis
!                    (parameter REBOUND==2)
!                    -----------------------------------------------------------
!                    Implementation of Bernd Muegge's conservative first-order
!                    upstream scheme for the age equation
!                    (parameters ADV_HOR==3 and ADV_VERT==3).
!                    -----------------------------------------------------------
!                    Antarctica, Scandinavia (Greenland see V2.2):
!                    - Precipitation instead of solid accumulation used as
!                      input.
!                    - Reeh's degree day model with statistical fluctuations,
!                      supplemented by explicit consideration of rainfall.
!                    - Net surface mass balance computed as precipitation minus
!                      runoff (previously solid accumulation minus melting).
!                    - Optional climate forcing by glacial index and prescribed
!                      LGM anomalies for surface temperature and precipitation
!                      (parameters TSURFACE==5 and ACCSURFACE==5).
!                    -----------------------------------------------------------
!                    Northern-hemisphere module no longer up-to-date, therefore
!                    deleted for the time being.
!                    -----------------------------------------------------------
!                    Compatibility check between version number of SICOPOLIS
!                    and version number in header file.
!              -----------------------------------------------------------------
!              V2.4.1: Antarctica, Scandinavia, Greenland:
!                      - Optional computation of surface melting (ablation)
!                        by the positive-degree-day (PDD) method
!                        or the linear-temperature-index (LTI) method
!                        (parameter ABLSURFACE).
!                      - Computation of positive degree days (subroutine pdd)
!                        changed from annual sinus cycle of the surface
!                        temperature to direct use of monthly-mean temperatures.
!                      ---------------------------------------------------------
!                      Scandinavia:
!                      - Climate input (surface temperature, precipitation)
!                        changed from seasonal to monthly-mean values.
!                      - Optional horizontal resolutions 40, 20 or 10 km.
!              -----------------------------------------------------------------
!              V2.5: Basal sliding now computed in new subroutine calc_vxy_b.
!                    Parameters of sliding law are set in the header file
!                    (previously in phys_para file). New options: Weertman
!                    exponents can be chosen freely, alternative use of
!                    full or reduced pressure in sliding law.
!                    -----------------------------------------------------------
!                    Antarctica, Greenland, Scandinavia:
!                    Two variants of PDD, either with instantaneous runoff of
!                    rainfall (ABLSURFACE==1), or with contribution of rainfall
!                    to formation of superimposed ice (ABLSURFACE==2).
!                    Linear-temperature-index (LTI) method is now chosen by
!                    ABLSURFACE==3.
!                    -----------------------------------------------------------
!                    North- and south-polar cap of Mars:
!                    In addition to the density and the heat conductivity,
!                    now the specific heat is also computed as a
!                    volume-fraction-weighed mean of the ice and dust values,
!                    respectively. Function c_val changed accordingly.
!                    -----------------------------------------------------------
!                    EISMINT Phase 2 Simplified Geometry Experiment:
!                    Input file surf_para_emtp2sge.dat and subroutine surf_para
!                    deleted. Parameters now defined in header file.
!                    -----------------------------------------------------------
!                    Minimum bedrock elevation for extent of marine ice
!                    defined by new variable z_mar. Can be chosen as a constant
!                    (MARGIN==2, old option), as proportional to the sea-level
!                    stand (MARGIN==3, new option) or as a piecewise linear
!                    function of sea level (MARGIN==4, new option; as proposed
!                    by Zweck and Huybrechts, 2005).
!              -----------------------------------------------------------------
!              V2.6: New module for ISMIP HEINO included (defined by macro
!                    HEINO).
!                    -----------------------------------------------------------
!                    Additional output of the basal frictional heating in the
!                    time-series file for the ice-core locations.
!                    -----------------------------------------------------------
!                    All constants throughout the program defined consistently
!                    as double-precision constants.
!                    -----------------------------------------------------------
!                    New global parameters pi, pi_inv, pi_180, pi_180_inv, eps.
!              -----------------------------------------------------------------
!              V2.7: Special treatment for the chasms of the Martian NPC (Chasma
!                    Borealis) and SPC (Chasma Australis) implemented (optional,
!                    parameter CHASM==2). Chasm regions defined in additional
!                    mask file -> larger geothermal heat flux and wind erosion
!                    rate active within prescribed time interval (specified in
!                    header file).
!                    -----------------------------------------------------------
!                    ISMIP HEINO: Time-series output for the sediment region
!                    ("Hudson Bay, Hudson Strait") in new file runname.sed.
!                    -----------------------------------------------------------
!                    Improved restart: Time derivatives dzs_dtau, dzm_dtau,
!                    dzb_dtau, dH_c_dtau and dH_t_dtau written to time-slice
!                    files and read by restart routine topography3. Further,
!                    Q_b_tot is computed as the sum of Q_bm and Q_tld in 
!                    initialization routine sico_init.
!                    -----------------------------------------------------------
!                    In addition to Q_b_tot (formerly called Q_mean), now also
!                    limitations of Q_bm and Q_tld with the same lower and upper
!                    limits, QB_MIN and QB_MAX, respectively.
!                    -----------------------------------------------------------
!                    In subroutine output2 Reinhard Calov's method for computing
!                    the ice equivalent of the water volume below the sea level
!                    (V_ges_redu) implemented.
!              -----------------------------------------------------------------
!              V2.7.1: Martian ice caps (NMARS, SMARS):
!                      New topography inputs based on the MOLA MEGDR data with
!                      64 pix/deg resolution. Optional horizontal resolutions
!                      20, 10 or 5 km.
!              -----------------------------------------------------------------
!              V2.7.2: x and y coordinates of the origin point (i,j) = (0,0)
!                      now defined in header file as X0 and Y0.
!                      ---------------------------------------------------------
!                      Optional computation of the evolution of the free
!                      surface as usual (parameter ZS_EVOL==1),
!                      or keeping the free surface fixed on the initial
!                      topography (parameter ZS_EVOL==0).
!                      ---------------------------------------------------------
!                      Greenland:
!                      Update of the ice-core positions according to Dahl-Jensen
!                      (pers. comm. 2006) and inclusion of the planned NEEM
!                      core.
!                      ---------------------------------------------------------
!                      Martian ice caps (NMARS, SMARS):
!                      New version of module instemp for the LIT scheme
!                      implemented.
!              -----------------------------------------------------------------
!              V2.8: Greenland:
!                    Optional special sliding law for the North-East Greenland
!                    Ice Stream NEGIS (parameter ICE_STREAM==2).
!                    Surface meltwater effect on basal sliding implemented.
!                    -----------------------------------------------------------
!                    ISMIP HEINO:
!                    Optional 25-km resolution implemented (in addition to the
!                    standard resolution of 50 km).
!                    -----------------------------------------------------------
!                    Antarctica:
!                    Re-formatting and re-naming of the input data (now same
!                    style and units like for Greenland). Topography input newly
!                    created from RAMPDEM V2 and BEDMAP data-sets, with an
!                    additional 20-km resolution implemented.
!              -----------------------------------------------------------------
!              V2.8.1: Northern-hemisphere module reintroduced, topography
!                      files newly created for horizontal resolutions of
!                      80, 40 and 20 km.
!                      ---------------------------------------------------------
!                      Explicit introduction of the evaporation (variable evap).
!                      In the modules for the terrestrial ice sheets, evap is
!                      set to zero, whereas in the modules for the Martian
!                      ice caps, evap accounts for the entire ablation.
!                      ---------------------------------------------------------
!                      ISMIP HEINO:
!                      Introduction of the option ANF_DAT==1, which defines the
!                      initial topography as a thin ice layer (with a thickness
!                      of 2*H_MIN) everywhere on the land area.
!              -----------------------------------------------------------------
!              V2.9: New module for Tibet included (defined by macro TIBET).
!                    -----------------------------------------------------------
!                    Greenland:
!                    Parameterization for the surface meltwater effect on basal
!                    sliding changed; additional dependence on ice thickness
!                    and possible nonlinear dependence on runoff.
!                    -----------------------------------------------------------
!                    Topography files (for zs, zb, zb0, mask) can now be
!                    specified in the header file.
!              -----------------------------------------------------------------
!              V3.0: Handling of ice shelf and marine ice dynamics.
!                    -----------------------------------------------------------
!                    SeaRISE set-up implemented.
!                    -----------------------------------------------------------
!                    Changing water load now considered for the LLRA model
!                    (REBOUND==1) in the same way as for the ELRA model.
!                    -----------------------------------------------------------
!                    Optional output of time-slice files in NetCDF format
!                    (controlled by parameter NETCDF).
!                    -----------------------------------------------------------
!                    Name of the physical-parameter file is now specified in
!                    the header file.
!                    -----------------------------------------------------------
!                    New module 'sico_types' for declarations of kind types
!                    introduced.
!                    -----------------------------------------------------------
!                    Subroutine 'calc_top':
!                    Diffusive smoothing scheme for zs_neu deleted.
!                    -----------------------------------------------------------
!                    Introduction of module Austfonna (ASF) available at 1, 2
!                    and 4 km horizontal resolution - Thorben Dunse.
!                    Output5.F90 writes time-series data at specific surface 
!                    points.   
!                    -----------------------------------------------------------
!                    Stereographic projection changed from a spherical planet
!                    to an ellipsoidal planet.   
!              -----------------------------------------------------------------
!              V3.1: Discretisation of the horizontal and vertical advection
!                    terms in the 3-d temperature and age equations:
!                    Improved upwind scheme implemented that uses interpolated
!                    velocities on the main grid (ADV_HOR==3, ADV_VERT==3).
!                    -----------------------------------------------------------
!                    Slightly improved boundary and transition conditions for
!                    the computation of the age.
!                    -----------------------------------------------------------
!                    Index order of 2-d and 3-d arrays in the time-slice output
!                    files *.nc (netCDF) and *.erg (native binary) changed:
!                    (j,i) -> (i,j), (kc/t/r,j,i) -> (i,j,kc/t/r).
!                    Output time units changed from seconds to years.
!                    -----------------------------------------------------------
!                    Routines that use Lis changed such that they are compatible
!                    with Lis version 1.4.13 and newer.
!                    -----------------------------------------------------------
!                    Solvers for systems of linear equations assembled in the
!                    new module sico_sle_solvers.
!                    -----------------------------------------------------------
!                    All subroutines and functions are now included in the main
!                    program sicopolis following a 'contains' statement.
!                    -----------------------------------------------------------
!                    Intent attributes inserted in all subroutines and
!                    functions.
!                    -----------------------------------------------------------
!                    New topography input data for Antarctica (BEDMAP2) and
!                    Greenland (Bamber Version 2).
!                    -----------------------------------------------------------
!                    New input data for the geothermal heat flux for Antarctica
!                    and Greenland [Michael Purucker (pers. comm. 2012), based
!                    on Fox Maule et al., (2005)].
!              *****************************************************************
!
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!-------- Include specification header --------

#include "sico_specs.h"

!-------- Include kind-type and global-variable modules --------

#include "subroutines/general/sico_types.F90"
#include "subroutines/general/sico_variables.F90"

!-------- Include module that contains the solvers for systems
!                                          of linear equations --------

#include "subroutines/general/sico_sle_solvers.F90"

!-------- Include NetCDF error capturing module --------

#if (NETCDF > 1)
#include "subroutines/general/nc_check.F90"
#endif

!-------- Include Bjoern Grieger's module for the surface-temperature
!         parameterisation for the Martian ice caps --------

#if defined(NMARS) 
#include "subroutines/nmars/instemp.f90"
#elif defined(SMARS)
#include "subroutines/smars/instemp.f90"
#endif

!-------------------------------------------------------------------------------
!> Main program of SICOPOLIS.
!<------------------------------------------------------------------------------
program sicopolis

!-------- Description of variables declared locally --------

!  i                 : Index for coordinate xi
!  j                 : Index for coordinate eta
!  kc                : Index for coordinate zeta_c
!  kt                : Index for coordinate zeta_t
!  kr                : Index for coordinate zeta_r
!                      (in the bedrock)
!  ios               : IOSTAT variable for files
!  itercount         : Counter for the time steps
!  ndat2d/3d         : Counters for the time-slice files
!  n_output          : For OUTPUT==2 number of time-slice files
!                      to be produced
!  dH_t_smooth(j,i)  : Amount of smoothing of H_t
!  delta_ts          : Time-dependent surface-temperature variation
!  glac_index        : Time-dependent glacial index
!  forcing_flag      : 1 - forcing by delta_ts, 2 - forcing by glac_index
!  z_sl              : Sea level
!  dzsl_dtau         : Derivative of z_sl by tau (time)
!  precip_mam_present(j,i) : Measured present spring precipitation
!  precip_jja_present(j,i) : Measured present summer precipitation
!  precip_son_present(j,i) : Measured present autumn precipitation
!  precip_djf_present(j,i) : Measured present winter precipitation
!  temp_mam_present(j,i)   : Present spring surface temperature
!                            from data
!  temp_jja_present(j,i)   : Present summer surface temperature
!                            from data
!  temp_son_present(j,i)   : Present autumn surface temperature
!                            from data
!  temp_djf_present(j,i)   : Present winter surface temperature
!                            from data
!  mean_accum        : Mean present accumulation over land
!  mean_accum_inv    : Inverse of mean_accum
!  time              : Current time of simulation
!  time_init         : Initial time of simulation
!  time_end          : Final time of simulation
!  dtime             : Time step for computation of velocity and
!                      topography
!  dtime_temp        : Time step for computation of temperature,
!                      water content and age of the ice
!  dtime_wss         : Time step for computation of isostatic steady-state
!                      displacement of the lithosphere (ELRA model)
!  dtime_out         : Time step for output of time-slice files
!  dtime_ser         : Time step for writing of data in time-series
!                      file
!  time_output(n)    : For OUTPUT==2 specified times for output of
!                      time-slice files
!  (.)0              : Quantity (.) before its conversion to MKS units
!  dxi               : Grid spacing in x-direction
!  deta              : Grid spacing in y-direction
!  dzeta_c           : Grid spacing in z-direction in cold ice
!                      (in sigma-coordinate zeta_c)
!  dzeta_t           : Grid spacing in z-direction in temperate ice
!                      (in sigma-coordinate zeta_t)
!  dzeta_r           : Grid spacing in z-direction in the bedrock
!                      (in sigma-coordinate zeta_r)
!  runname           : Name of simulation
!  anfdatname        : Name of initial-value file
!                      (only for ANF_DAT==3)
!  datname           : Auxiliary string for generation of file names

!-------- Declaration of variables --------

use sico_types
use sico_variables

implicit none

integer(i4b) :: ndat2d, ndat3d
integer(i4b) :: n_output
integer(i4b), dimension((IMAX+1)*(JMAX+1)) :: ii, jj
integer(i4b), dimension(0:JMAX,0:IMAX)     :: nn
real(dp) :: delta_ts, glac_index
real(dp) :: mean_accum, mean_accum_inv
real(dp) :: dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser
real(dp) :: time, time_init, time_end
real(dp), dimension(100) :: time_output
real(dp) :: dxi, deta, dzeta_c, dzeta_t, dzeta_r
real(dp) :: z_sl, dzsl_dtau, z_mar
character(len=100) :: runname

!-------- Initialisations --------

call sico_init(delta_ts, glac_index, &
     mean_accum, mean_accum_inv, &
     dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, &
     time, time_init, time_end, time_output, &
     dxi, deta, dzeta_c, dzeta_t, dzeta_r, &
     z_sl, dzsl_dtau, z_mar, &
     ii, jj, nn, &
     ndat2d, ndat3d, n_output, &
     runname)

!-------- Main loop --------

call sico_main_loop(delta_ts, glac_index, &
     mean_accum, mean_accum_inv, &
     dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, &
     time, time_init, time_end, time_output, &
     dxi, deta, dzeta_c, dzeta_t, dzeta_r, &
     z_sl, dzsl_dtau, z_mar, &
     ii, jj, nn, &
     ndat2d, ndat3d, n_output, &
     runname)

!-------- Endings --------

call sico_end()

!-------- Include subroutines --------

contains

#if defined(NHEM)
#include "subroutines/nhem/sico_init.F90"
#elif defined(SCAND)
#include "subroutines/scand/sico_init.F90"
#elif defined(GRL)
#include "subroutines/grl/sico_init.F90"
#elif defined(NMARS)
#include "subroutines/nmars/sico_init.F90"
#elif defined(SMARS)
#include "subroutines/smars/sico_init.F90"
#elif defined(EMTP2SGE)
#include "subroutines/emtp2sge/sico_init.F90"
#elif defined(HEINO)
#include "subroutines/heino/sico_init.F90"
#elif defined(ANT)
#include "subroutines/ant/sico_init.F90"
#elif defined(TIBET)
#include "subroutines/tibet/sico_init.F90"
#elif defined(ASF)
#include "subroutines/asf/sico_init.F90"
#endif

#include "subroutines/general/sico_main_loop.F90"
#include "subroutines/general/sico_end.F90"

#if defined(NHEM)
#include "subroutines/nhem/phys_para.F90"
#elif defined(SCAND)
#include "subroutines/scand/phys_para.F90"
#elif defined(GRL)
#include "subroutines/grl/phys_para.F90"
#elif defined(NMARS)
#include "subroutines/nmars/phys_para.F90"
#elif defined(SMARS)
#include "subroutines/smars/phys_para.F90"
#elif defined(EMTP2SGE)
#include "subroutines/emtp2sge/phys_para.F90"
#elif defined(HEINO)
#include "subroutines/heino/phys_para.F90"
#elif defined(ANT)
#include "subroutines/ant/phys_para.F90"
#elif defined(TIBET)
#include "subroutines/tibet/phys_para.F90"
#elif defined(ASF)
#include "subroutines/asf/phys_para.F90"
#endif

#include "subroutines/general/read_phys_para_value.F90"
#include "subroutines/general/read_kei.F90"

#if defined(NHEM)
#include "subroutines/nhem/topography2.F90"
#include "subroutines/nhem/topography3.F90"
#if (NETCDF > 1)
#include "subroutines/nhem/topography3_nc.F90"
#endif
#elif defined(SCAND)
#include "subroutines/scand/topography2.F90"
#include "subroutines/scand/topography3.F90"
#if (NETCDF > 1)
#include "subroutines/scand/topography3_nc.F90"
#endif
#elif defined(GRL)
#include "subroutines/grl/topography1.F90"
#include "subroutines/grl/topography2.F90"
#include "subroutines/grl/topography3.F90"
#if (NETCDF > 1)
#include "subroutines/grl/topography3_nc.F90"
#if (ZS_EVOL==2)
#include "subroutines/grl/read_target_topo_nc.F90"
#endif
#endif
#elif defined(NMARS)
#include "subroutines/nmars/topography1.F90"
#include "subroutines/nmars/topography2.F90"
#include "subroutines/nmars/topography3.F90"
#if (NETCDF > 1)
#include "subroutines/nmars/topography3_nc.F90"
#endif
#elif defined(SMARS)
#include "subroutines/smars/topography1.F90"
#include "subroutines/smars/topography2.F90"
#include "subroutines/smars/topography3.F90"
#if (NETCDF > 1)
#include "subroutines/smars/topography3_nc.F90"
#endif
#elif defined(EMTP2SGE)
#include "subroutines/emtp2sge/topography2.F90"
#include "subroutines/emtp2sge/topography3.F90"
#if (NETCDF > 1)
#include "subroutines/emtp2sge/topography3_nc.F90"
#endif
#elif defined(HEINO)
#include "subroutines/heino/topography1.F90"
#include "subroutines/heino/topography2.F90"
#include "subroutines/heino/topography3.F90"
#if (NETCDF > 1)
#include "subroutines/heino/topography3_nc.F90"
#endif
#elif defined(ANT)
#include "subroutines/ant/topography1.F90"
#include "subroutines/ant/topography2.F90"
#include "subroutines/ant/topography3.F90"
#if (NETCDF > 1)
#include "subroutines/ant/topography3_nc.F90"
#if (ZS_EVOL==2)
#include "subroutines/ant/read_target_topo_nc.F90"
#endif
#endif
#elif defined(TIBET)
#include "subroutines/tibet/topography2.F90"
#include "subroutines/tibet/topography3.F90"
#if (NETCDF > 1)
#include "subroutines/tibet/topography3_nc.F90"
#endif
#elif defined(ASF)
#include "subroutines/asf/topography1.F90"
#include "subroutines/asf/topography2.F90"
#include "subroutines/asf/topography3.F90"
#if (NETCDF > 1)
#include "subroutines/asf/topography3_nc.F90"
#endif
#endif

#if defined(NHEM)
#include "subroutines/nhem/boundary.F90"
#elif defined(SCAND)
#include "subroutines/scand/boundary.F90"
#elif defined(GRL)
#include "subroutines/grl/boundary.F90"
#elif defined(NMARS)
#include "subroutines/nmars/boundary.F90"
#elif defined(SMARS)
#include "subroutines/smars/boundary.F90"
#elif defined(EMTP2SGE)
#include "subroutines/emtp2sge/boundary.F90"
#elif defined(HEINO)
#include "subroutines/heino/boundary.F90"
#elif defined(ANT)
#include "subroutines/ant/boundary.F90"
#elif defined(TIBET)
#include "subroutines/tibet/boundary.F90"
#elif defined(ASF)
#include "subroutines/asf/boundary.F90"
#endif

#include "subroutines/general/mask_update.F90"
#include "subroutines/general/pdd.F90"
#include "subroutines/general/erfcc.F90"
#include "subroutines/general/calc_temp_poly.F90"
#include "subroutines/general/calc_temp_cold.F90"
#include "subroutines/general/calc_temp1.F90"
#include "subroutines/general/calc_temp2.F90"
#include "subroutines/general/calc_temp3.F90"
#include "subroutines/general/calc_temp_r.F90"
#include "subroutines/general/shift_cts_upward.F90"
#include "subroutines/general/shift_cts_downward.F90"
#include "subroutines/general/calc_enhance.F90"
#include "subroutines/general/calc_vxy_b_sia.F90"
#include "subroutines/general/calc_vxy_sia.F90"
#include "subroutines/general/calc_vz_sia.F90"
#if MARGIN==3
#include "subroutines/general/calc_temp_ssa.F90"
#include "subroutines/general/calc_vxy_ssa.F90"
#include "subroutines/general/calc_vxy_ssa_matrix.F90"
#include "subroutines/general/calc_vis_ssa.F90"
#include "subroutines/general/calc_vz_ssa.F90"
#include "subroutines/general/viscosity.F90"
#endif
#include "subroutines/general/calc_qbm.F90"
#include "subroutines/general/calc_water_bas.F90"
#include "subroutines/general/calc_top.F90"
#include "subroutines/general/calc_elra.F90"
#include "subroutines/general/topograd_1.F90"
#include "subroutines/general/topograd_2.F90"
#include "subroutines/general/calc_temp_melt.F90"
#include "subroutines/general/calc_temp_bas.F90"
#include "subroutines/general/creep.F90"
#include "subroutines/general/ratefac.F90"
#include "subroutines/general/ratefac_t.F90"
#include "subroutines/general/kappa_val.F90"
#include "subroutines/general/c_val.F90"
#include "subroutines/general/metric.F90"
#include "subroutines/general/geo_coord.F90"
#include "subroutines/general/num_coord.F90"
#include "subroutines/general/output1.F90"
#include "subroutines/general/output2.F90"
#include "subroutines/general/output3.F90"
#include "subroutines/general/output4.F90"

#if defined(ASF)
#include "subroutines/asf/output5.F90"
#endif

#if (NETCDF > 1)
#include "subroutines/general/output_nc.F90"
#endif

#include "subroutines/general/borehole.F90"

!-------- End of program --------

end program sicopolis

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!                       End of sicopolis.F90
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!