sicopolis.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Program :  s i c o p o l i s
!             (SImulation COde for POLythermal Ice Sheets)
!
#define       MODEL_SICOPOLIS
#define       VERSION '3.3'
#define       DATE    '2017-08-03'
!
!> @mainpage
!!
!! @section Description
!!
!! SICOPOLIS (SImulation COde for POLythermal Ice Sheets) is a 3-d
!! dynamic/thermodynamic model that simulates the evolution of large ice
!! sheets and ice caps. It was originally created by Greve (1997a,b) 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, the polar ice caps of the planet Mars and others.
!!
!! The model is based on the shallow ice approximation for grounded ice, the
!! shallow shelf approximation for floating ice (e.g., Greve and Blatter 2009)
!! and, optionally, hybrid shallow-ice--shelfy-stream dynamics for ice streams
!! (Bernales et al. 2017). It is coded in Fortran and uses finite difference
!! discretisation on a staggered (Arakawa C) grid, the velocity
!! components being taken between grid points. A variety of different
!! thermodynamics solvers are available, namely the polythermal two-layer
!! method, two versions of the one-layer enthalpy method, the cold-ice
!! method and the isothermal method (Greve and Blatter 2016).
!!
!! The coding is based on a consequent low-tech philosophy. All
!! structures are kept as simple as possible, and advanced coding
!! techniques are only employed where it is deemed appropriate. The use
!! of external libraries is kept at an absolute minimum. In fact,
!! SICOPOLIS can be run without external libraries at all, which makes
!! the installation very easy and fast.
!!
!! 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.
!!      
!! Main 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.
!!      
!! References:
!! @li Bernales, J., I. Rogozhina, R. Greve and M. Thomas. 2017.\n
!!     Comparison of hybrid schemes for the combination of
!!     shallow approximations in numerical simulations of the
!!     Antarctic Ice Sheet.\n
!!     Cryosphere 11 (1), 247-265.
!! @li Greve, R. 1997a.\n
!!     A continuum-mechanical formulation for shallow polythermal ice sheets.\n
!!     Phil. Trans. R. Soc. A 355 (1726), 921-974.
!! @li Greve, R. 1997b.\n
!!     Application of a polythermal three-dimensional ice sheet model to the
!!     Greenland ice sheet: Response to steady-state and transient climate
!!     scenarios.\n
!!     J. Climate 10 (5), 901-918.
!! @li Greve, R. and H. Blatter. 2009.\n
!!     Dynamics of Ice Sheets and Glaciers.\n
!!     Springer, Berlin, Germany etc., 287 pp.
!! @li Greve, R. and H. Blatter. 2016.\n
!!     Comparison of thermodynamics solvers in the polythermal ice sheet model
!!     SICOPOLIS.\n
!!     Polar Sci. 10 (1), 11-23.
!! @li SICOPOLIS website: <http://www.sicopolis.net/>
!! @li Changelog: <http://tinyurl.com/commit-logs-sico-trunk/>
!!
!! @section Copyright
!!
!! Copyright 2009-2017 Ralf Greve\n
!! (with contributions by Jorge Bernales, Heinz Blatter, Reinhard Calov,
!! Thorben Dunse, Ben Galton-Fenzi, Thomas Goelles, Philipp Hancke,
!! Nina Kirchner, Sascha Knell, Alex Robinson, Tatsuru Sato, Malte Thoma,
!! Roland Warner)
!!
!! @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-2017 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/>.
!<
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

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

#include "sico_specs.h"

!-------- Include header for the Library of Iterative Solvers Lis
!                                               (only if required) --------

#if (CALCTHK==3 || CALCTHK==6 || MARGIN==3 || DYNAMICS==2)
#include "lisf.h"
#endif

!-------- Include modules --------

#include "subroutines/general/sico_types_m.F90"
#include "subroutines/general/sico_variables_m.F90"

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

#include "subroutines/general/ice_material_properties_m.F90"
#include "subroutines/general/stereo_proj_m.F90"
#include "subroutines/general/metric_m.F90"
#include "subroutines/general/sico_maths_m.F90"
#include "subroutines/general/compare_float_m.F90"

#if (NETCDF==2)
#include "subroutines/general/nc_check_m.F90"
#endif

#include "subroutines/general/read_m.F90"
#include "subroutines/general/init_temp_water_age_m.F90"

#include "subroutines/general/mask_update_sea_level_m.F90"
#include "subroutines/general/pdd_m.F90"

#if ((MARGIN==2) \
      && (marine_ice_formation==2) \
      && (marine_ice_calving==9))
#include "subroutines/general/calving_underwater_ice_m.F90"
#endif

#if (defined(GRL))
#if (DISC>0)
#include "subroutines/grl/discharge_workers_m.F90"
#endif
#endif

#include "subroutines/general/calc_vxy_m.F90"
#include "subroutines/general/calc_vz_m.F90"
#include "subroutines/general/calc_dxyz_m.F90"

#include "subroutines/general/calc_gia_m.F90"
#include "subroutines/general/topograd_m.F90"
#include "subroutines/general/calc_thk_m.F90"

#include "subroutines/general/enth_temp_omega_m.F90"

#if (CALCMOD==0 || CALCMOD==1 || CALCMOD==-1)
#include "subroutines/general/calc_temp_m.F90"
#elif (CALCMOD==2 || CALCMOD==3)
#include "subroutines/general/calc_temp_enth_m.F90"
#endif

#include "subroutines/general/calc_enhance_m.F90"

#if (BASAL_HYDROLOGY==1)
#include "subroutines/non-foss/hydro_m.F90"
#endif

#if (defined(NMARS) || defined(SMARS))
#include "subroutines/general/mars_instemp_m.f90"
#endif

#include "subroutines/general/calc_temp_melt_bas_m.F90"
#include "subroutines/general/calc_bas_melt_m.F90"
#include "subroutines/general/calc_thk_water_bas_m.F90"

#include "subroutines/general/output_m.F90"

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

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

#include "subroutines/general/sico_main_loop_m.F90"
#include "subroutines/general/sico_end_m.F90"

!-------------------------------------------------------------------------------
!> 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
!  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 the upper (kc) ice domain
!                      (in sigma-coordinate zeta_c)
!  dzeta_t           : Grid spacing in z-direction in the lower (kt) ice domain
!                      (in sigma-coordinate zeta_t)
!  dzeta_r           : Grid spacing in z-direction in the bedrock (kr) domain
!                      (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_m
use sico_variables_m
use sico_vars_m

use sico_init_m
use sico_main_loop_m
use sico_end_m

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
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, &
     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, &
     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()

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

end program sicopolis

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