sico_init_m.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Module :  s i c o _ i n i t _ m
!
!> @file
!!
!! Initialisations for 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/>.
!<
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!-------------------------------------------------------------------------------
!> Initialisations for SICOPOLIS.
!<------------------------------------------------------------------------------
module sico_init_m

  use sico_types_m
  use sico_variables_m
  use sico_vars_m

  implicit none

  public

contains

!-------------------------------------------------------------------------------
!> Main routine of sico_init_m: Initialisations for SICOPOLIS.
!<------------------------------------------------------------------------------
subroutine 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)

  use compare_float_m
  use ice_material_properties_m, only : ice_mat_eqs_pars
  use enth_temp_omega_m, only : calc_c_int_table, calc_c_int_inv_table, &
                                enth_fct_temp_omega

  use read_m, only : read_phys_para

  use boundary_m
  use init_temp_water_age_m
  use calc_enhance_m
  use calc_vxy_m
  use calc_vz_m
  use calc_dxyz_m
  use calc_temp_melt_bas_m

  use output_m

implicit none

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

integer(i4b)       :: i, j, kc, kt, kr, m, n, ir, jr
integer(i4b)       :: ios, ios1, ios2, ios3, ios4
integer(i4b)       :: ierr
real(dp)           :: dtime0, dtime_temp0, dtime_wss0, dtime_out0, dtime_ser0
real(dp)           :: time_init0, time_end0, time_output0(100)
real(dp)           :: d_dummy
character(len=100) :: anfdatname
character(len=256) :: filename_with_path
character(len=256) :: shell_command
character          :: ch_dummy
logical            :: flag_init_output, flag_3d_output

character(len=64), parameter :: fmt1  = '(a)', &
                                fmt2  = '(a,i4)', &
                                fmt2a = '(a,i0)', &
                                fmt3  = '(a,es12.4)'

write(unit=6, fmt='(a)') ' '
write(unit=6, fmt='(a)') ' -------- sico_init --------'
write(unit=6, fmt='(a)') ' '

!-------- Name of the computational domain --------

#if (defined(ANT))
ch_domain_long  = 'Antarctica'
ch_domain_short = 'ant'

#elif (defined(ASF))
ch_domain_long  = 'Austfonna'
ch_domain_short = 'asf'

#elif (defined(EMTP2SGE))
ch_domain_long  = 'EISMINT Phase 2 Simplified Geometry Experiment'
ch_domain_short = 'emtp2sge'

#elif (defined(GRL))
ch_domain_long  = 'Greenland'
ch_domain_short = 'grl'

#elif (defined(NHEM))
ch_domain_long  = 'Northern hemisphere'
ch_domain_short = 'nhem'

#elif (defined(SCAND))
ch_domain_long  = 'Scandinavia and Eurasia'
ch_domain_short = 'scand'

#elif (defined(TIBET))
ch_domain_long  = 'Tibet'
ch_domain_short = 'tibet'

#elif (defined(NMARS))
ch_domain_long  = 'North polar cap of Mars'
ch_domain_short = 'nmars'

#elif (defined(SMARS))
ch_domain_long  = 'South polar cap of Mars'
ch_domain_short = 'smars'

#elif (defined(XYZ))
ch_domain_long  = 'XYZ'
ch_domain_short = 'xyz'
#if (defined(HEINO))
ch_domain_long  = trim(ch_domain_long)//'/ISMIP HEINO'
#endif

#else
stop ' >>> sico_init: No valid domain specified!'
#endif

!-------- Some initial values --------

n_output = 0

dtime       = 0.0_dp
dtime_temp  = 0.0_dp
dtime_wss   = 0.0_dp
dtime_out   = 0.0_dp
dtime_ser   = 0.0_dp

time        = 0.0_dp
time_init   = 0.0_dp
time_end    = 0.0_dp
time_output = 0.0_dp

!-------- Initialisation of the Library of Iterative Solvers Lis,
!                                                     if required --------

#if (CALCTHK==3 || CALCTHK==6 || MARGIN==3 || DYNAMICS==2)
  call lis_initialize(ierr)
#endif

!-------- Read physical parameters --------

call read_phys_para()

call ice_mat_eqs_pars(rf, r_t, kappa, c, -190, 10)

#if (!defined(SURFACE_FORCING) || SURFACE_FORCING==1)

temp_min  = temp_min                         ! deg C
s_t       = s_t       *1.0e-03_dp            ! K/km -> K/m
x_hat     = x_hat     *1.0e+03_dp            ! km -> m
y_hat     = y_hat     *1.0e+03_dp            ! km -> m
b_max     = b_max     /year_sec              ! m/a -> m/s
s_b       = s_b       *1.0e-03_dp/year_sec   ! m/(a*km) -> 1/s
eld       = eld       *1.0e+03_dp            ! km -> m

#elif (SURFACE_FORCING==2)

temp_0    = temp_0                           ! deg C
gamma_t   = gamma_t   *1.0e-03_dp            ! K/km -> K/m
s_0       = s_0       /year_sec              ! m/a -> m/s
m_0       = m_0       *1.0e-03_dp/year_sec   ! m/(a*km) -> 1/s
ela       = ela       *1.0e+03_dp            ! km -> m

#else
stop ' >>> sico_init: SURFACE_FORCING must be either 1 or 2!'
#endif

!  ------ Some auxiliary quantities required for the enthalpy method

call calc_c_int_table(c, -190, 10, l)
call calc_c_int_inv_table()

!-------- Compatibility check of the SICOPOLIS version with the header file

if ( trim(version) /= trim(sico_version) ) &
   stop ' >>> sico_init: SICOPOLIS version not compatible with header file!'

!-------- Check whether the dynamics and thermodynamics modes are defined

#if (!defined(DYNAMICS))
stop ' >>> sico_init: DYNAMICS not defined in the header file!'
#endif

#if (!defined(CALCMOD))
stop ' >>> sico_init: CALCMOD not defined in the header file!'
#endif

#if (defined(ENTHMOD))
write(6, fmt='(a)') ' >>> sico_init: ENTHMOD must not be defined any more.'
write(6, fmt='(a)') '                Please update your header file!'
stop
#endif

!-------- Compatibility check of the horizontal resolution with the
!         number of grid points --------

#if (GRID==0)

if (approx_equal(dx, 5.0_dp, eps_sp_dp)) then

   if ((imax /= 300).or.(jmax /= 300)) then
      stop ' >>> sico_init: IMAX and/or JMAX wrong!'
   end if

else if (approx_equal(dx, 10.0_dp, eps_sp_dp)) then

   if ((imax /= 150).or.(jmax /= 150)) then
      stop ' >>> sico_init: IMAX and/or JMAX wrong!'
   end if

else if (approx_equal(dx, 25.0_dp, eps_sp_dp)) then

   if ((imax /= 60).or.(jmax /= 60)) then
      stop ' >>> sico_init: IMAX and/or JMAX wrong!'
   end if

else if (approx_equal(dx, 75.0_dp, eps_sp_dp)) then

   if ((imax /= 20).or.(jmax /= 20)) then
      stop ' >>> sico_init: IMAX and/or JMAX wrong!'
   end if

else
      stop ' >>> sico_init: DX wrong!'
end if

#elif (GRID==1)

stop ' >>> sico_init: GRID==1 not allowed for this application!'

#elif (GRID==2)

stop ' >>> sico_init: GRID==2 not allowed for this application!'

#endif

!-------- Compatibility check of the thermodynamics mode
!         (cold vs. polythermal vs. enthalpy method)
!         and the number of grid points in the lower (kt) ice domain --------

#if (CALCMOD==0 || CALCMOD==2 || CALCMOD==3 || CALCMOD==-1)

if (ktmax > 2) then
   write(6, fmt='(a)') ' >>> sico_init: For options CALCMOD==0, 2, 3 or -1,'
   write(6, fmt='(a)') '                the separate kt domain is redundant.'
   write(6, fmt='(a)') '                Therefore, consider setting KTMAX to 2.'
   write(6, fmt='(a)') ' '
end if

#endif

!-------- Compatibility check of discretization schemes for the horizontal and
!         vertical advection terms in the temperature and age equations --------

#if (ADV_HOR==1)
stop ' >>> sico_init: Option ADV_HOR==1 (central differences) not defined!'
#endif

!-------- Check whether for the shallow shelf
!               or shelfy stream approximation
!                  the chosen grid is Cartesian coordinates
!                             without distortion correction (GRID==0) --------

#if ((MARGIN==3 && DYNAMICS==1) || DYNAMICS==2)   /* requires SSA or SStA */
#if (GRID != 0)
write(6, fmt='(a)') ' >>> sico_init: Distortion correction not implemented'
write(6, fmt='(a)') '                for the shallow shelf approximation (SSA)'
write(6, fmt='(a)') '                or the shelfy stream approximation (SStA)'
write(6, fmt='(a)') '                -> GRID==0 required!'
stop
#endif
#endif

!-------- Setting of forcing flag --------

forcing_flag = 1   ! forcing by delta_ts

!-------- Initialization of numerical time steps --------

dtime0      = dtime0
dtime_temp0 = dtime_temp0

!-------- Further initializations --------

dzeta_c = 1.0_dp/real(kcmax,dp)
dzeta_t = 1.0_dp/real(ktmax,dp)
dzeta_r = 1.0_dp/real(krmax,dp)

ndat2d = 1
ndat3d = 1

!-------- General abbreviations --------

!  ------ kc domain

if (deform >= eps) then

   flag_aa_nonzero = .true.   ! non-equidistant grid

   aa = deform
   ea = exp(aa)

   kc=0
   zeta_c(kc)         = 0.0_dp
   eaz_c(kc)          = 1.0_dp
   eaz_c_quotient(kc) = 0.0_dp

   do kc=1, kcmax-1
      zeta_c(kc) = kc*dzeta_c
      eaz_c(kc)  = exp(aa*zeta_c(kc))
      eaz_c_quotient(kc) = (eaz_c(kc)-1.0_dp)/(ea-1.0_dp)
   end do

   kc=kcmax
   zeta_c(kc)         = 1.0_dp
   eaz_c(kc)          = exp(aa)
   eaz_c_quotient(kc) = 1.0_dp

else

   flag_aa_nonzero = .false.   ! equidistant grid

   aa = 0.0_dp
   ea = 1.0_dp

   kc=0
   zeta_c(kc)         = 0.0_dp
   eaz_c(kc)          = 1.0_dp
   eaz_c_quotient(kc) = 0.0_dp

   do kc=1, kcmax-1
      zeta_c(kc) = kc*dzeta_c
      eaz_c(kc)  = 1.0_dp
      eaz_c_quotient(kc) = zeta_c(kc)
   end do

   kc=kcmax
   zeta_c(kc)         = 1.0_dp
   eaz_c(kc)          = 1.0_dp
   eaz_c_quotient(kc) = 1.0_dp

end if

!  ------ kt domain

kt=0
zeta_t(kt) = 0.0_dp

do kt=1, ktmax-1
   zeta_t(kt) = kt*dzeta_t
end do

kt=ktmax
zeta_t(kt) = 1.0_dp

!  ------ kr domain

kr=0
zeta_r(kr) = 0.0_dp

do kr=1, krmax-1
   zeta_r(kr) = kr*dzeta_r
end do

kr=krmax
zeta_r(kr) = 1.0_dp

!-------- Construction of a vector (with index n) from a 2-d array
!         (with indices i, j) by diagonal numbering --------

n=1

do m=0, imax+jmax
   do i=m, 0, -1
      j = m-i
      if ((i <= imax).and.(j <= jmax)) then
         ii(n)   = i
         jj(n)   = j
         nn(j,i) = n
         n=n+1
      end if
   end do
end do

!-------- Specification of current simulation --------

runname    = runname
anfdatname = anfdatname

#if (defined(YEAR_ZERO))
year_zero  = year_zero
#else
year_zero  = 2000.0_dp   ! default value 2000 CE
#endif

time_init0 = time_init0
time_end0  = time_end0
dtime_ser0 = dtime_ser0

#if (OUTPUT==1 || OUTPUT==3)
dtime_out0 = dtime_out0
#endif

#if (OUTPUT==2 || OUTPUT==3)
n_output         = n_output
time_output0( 1) = time_out0_01
time_output0( 2) = time_out0_02
time_output0( 3) = time_out0_03
time_output0( 4) = time_out0_04
time_output0( 5) = time_out0_05
time_output0( 6) = time_out0_06
time_output0( 7) = time_out0_07
time_output0( 8) = time_out0_08
time_output0( 9) = time_out0_09
time_output0(10) = time_out0_10
time_output0(11) = time_out0_11
time_output0(12) = time_out0_12
time_output0(13) = time_out0_13
time_output0(14) = time_out0_14
time_output0(15) = time_out0_15
time_output0(16) = time_out0_16
time_output0(17) = time_out0_17
time_output0(18) = time_out0_18
time_output0(19) = time_out0_19
time_output0(20) = time_out0_20
#endif

!-------- Write log file --------

shell_command = 'if [ ! -d'
shell_command = trim(shell_command)//' '//outpath
shell_command = trim(shell_command)//' '//'] ; then mkdir'
shell_command = trim(shell_command)//' '//outpath
shell_command = trim(shell_command)//' '//'; fi'
call system(trim(shell_command))
     ! Check whether directory OUTPATH exists. If not, it is created.

filename_with_path = trim(outpath)//'/'//trim(runname)//'.log'

open(10, iostat=ios, file=trim(filename_with_path), status='new')

if (ios /= 0) stop ' >>> sico_init: Error when opening the log file!'

write(10, fmt=trim(fmt1)) 'Computational domain:'
write(10, fmt=trim(fmt1)) trim(ch_domain_long)
write(10, fmt=trim(fmt1)) ' '

write(10, fmt=trim(fmt2a)) 'GRID = ', grid
write(10, fmt=trim(fmt1)) ' '

write(10, fmt=trim(fmt2)) 'imax  =', imax
write(10, fmt=trim(fmt2)) 'jmax  =', jmax
write(10, fmt=trim(fmt2)) 'kcmax =', kcmax
write(10, fmt=trim(fmt2)) 'ktmax =', ktmax
write(10, fmt=trim(fmt2)) 'krmax =', krmax
write(10, fmt=trim(fmt1)) ' '

write(10, fmt=trim(fmt3)) 'a =', aa
write(10, fmt=trim(fmt1)) ' '

#if (GRID==0 || GRID==1)
write(10, fmt=trim(fmt3)) 'x0      =', x0
write(10, fmt=trim(fmt3)) 'y0      =', y0
write(10, fmt=trim(fmt3)) 'dx      =', dx
#elif (GRID==2)
write(10, fmt=trim(fmt3)) 'dlambda =', dlambda
write(10, fmt=trim(fmt3)) 'dphi    =', dphi
#endif
write(10, fmt=trim(fmt1)) ' '

write(10, fmt=trim(fmt3)) 'year_zero  =', year_zero
write(10, fmt=trim(fmt3)) 'time_init  =', time_init0
write(10, fmt=trim(fmt3)) 'time_end   =', time_end0
write(10, fmt=trim(fmt3)) 'dtime      =', dtime0
write(10, fmt=trim(fmt3)) 'dtime_temp =', dtime_temp0
write(10, fmt=trim(fmt1)) ' '

write(10, fmt=trim(fmt2)) 'ANF_DAT =', anf_dat
#if (defined(LAND_OCEAN_TRANSITION_WIDTH))
write(10, fmt=trim(fmt3)) 'land_ocean_transition_width =', &
                           land_ocean_transition_width
#endif
#if (ANF_DAT==1 && defined(TEMP_INIT))
write(10, fmt=trim(fmt2)) 'TEMP_INIT =', temp_init
#endif
#if (ANF_DAT==3)
write(10, fmt=trim(fmt1)) 'Initial-value file = '//anfdatname
write(10, fmt=trim(fmt1)) 'Path to initial-value file = '//anfdatpath
#endif
write(10, fmt=trim(fmt1)) ' '

write(10, fmt=trim(fmt1)) 'Physical-parameter file = '//phys_para_file
write(10, fmt=trim(fmt1)) ' '

#if (CALCTHK==2 || CALCTHK==3 || CALCTHK==5 || CALCTHK==6)
write(10, fmt=trim(fmt3))  'ovi_weight   =', ovi_weight
#if (CALCTHK==2 || CALCTHK==5)
write(10, fmt=trim(fmt3))  'omega_sor    =', omega_sor
#if (ITER_MAX_SOR>0)
write(10, fmt=trim(fmt2a)) 'iter_max_sor = ', iter_max_sor
#endif
#endif
write(10, fmt=trim(fmt1)) ' '
#endif

#if (defined(SURFACE_FORCING))
write(10, fmt=trim(fmt2)) 'SURFACE_FORCING =', surface_forcing
#endif
#if (!defined(SURFACE_FORCING) || SURFACE_FORCING==1)
write(10, fmt=trim(fmt3)) 'temp_min =', temp_min
write(10, fmt=trim(fmt3)) 's_t      =', s_t
write(10, fmt=trim(fmt3)) 'x_hat    =', x_hat
write(10, fmt=trim(fmt3)) 'y_hat    =', y_hat
write(10, fmt=trim(fmt3)) 'b_max    =', b_max
write(10, fmt=trim(fmt3)) 's_b      =', s_b
write(10, fmt=trim(fmt3)) 'eld      =', eld
#elif (SURFACE_FORCING==2)
write(10, fmt=trim(fmt3)) 'temp_0   =', temp_0
write(10, fmt=trim(fmt3)) 'gamma_t  =', gamma_t
write(10, fmt=trim(fmt3)) 's_0      =', s_0
write(10, fmt=trim(fmt3)) 'm_0      =', m_0
write(10, fmt=trim(fmt3)) 'ela      =', ela
#endif
write(10, fmt=trim(fmt1)) ' '

#if (TSURFACE==1)
write(10, fmt=trim(fmt3)) 'delta_ts0      =', delta_ts0
#elif (TSURFACE==3)
write(10, fmt=trim(fmt3)) 'sine_amplit    =', sine_amplit
write(10, fmt=trim(fmt3)) 'sine_period    =', sine_period
#elif (TSURFACE==4)
write(10, fmt=trim(fmt1)) 'GRIP file      = '//grip_temp_file
write(10, fmt=trim(fmt3)) 'grip_temp_fact =', grip_temp_fact
#endif

write(10, fmt=trim(fmt2a)) 'SEA_LEVEL  = ', sea_level
#if (SEA_LEVEL==1)
write(10, fmt=trim(fmt3)) 'z_sl0          =', z_sl0
#elif (SEA_LEVEL==3)
write(10, fmt=trim(fmt1)) 'sea-level file = '//sea_level_file
#endif
write(10, fmt=trim(fmt1)) ' '

#if (MARGIN==2)
#if (MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3)
write(10, fmt=trim(fmt3)) 'z_mar          =', z_mar
write(10, fmt=trim(fmt1)) ' '
#elif (MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5 \
       || marine_ice_calving==6 || marine_ice_calving==7)
write(10, fmt=trim(fmt3)) 'fact_z_mar     =', fact_z_mar
write(10, fmt=trim(fmt1)) ' '
#elif (MARINE_ICE_FORMATION==2 && MARINE_ICE_CALVING==9)
write(10, fmt=trim(fmt3)) 'calv_uw_coeff  =', calv_uw_coeff
write(10, fmt=trim(fmt3)) 'r1_calv_uw     =', r1_calv_uw
write(10, fmt=trim(fmt3)) 'r2_calv_uw     =', r2_calv_uw
write(10, fmt=trim(fmt1)) ' '
#endif
#elif (MARGIN==3)
#if (ICE_SHELF_CALVING==2)
write(10, fmt=trim(fmt3)) 'H_calv          =', h_calv
write(10, fmt=trim(fmt1)) ' '
#endif
#endif

#if (defined(BASAL_HYDROLOGY))
write(10, fmt=trim(fmt2a)) 'BASAL_HYDROLOGY = ', basal_hydrology
#endif
write(10, fmt=trim(fmt2a)) 'SLIDE_LAW  = ', slide_law
write(10, fmt=trim(fmt3)) 'c_slide     =', c_slide
write(10, fmt=trim(fmt3)) 'gamma_slide =', gamma_slide
write(10, fmt=trim(fmt2a)) 'p_weert     = ', p_weert
write(10, fmt=trim(fmt2a)) 'q_weert     = ', q_weert
#if (defined(TIME_RAMP_UP_SLIDE))
write(10, fmt=trim(fmt3)) 'time_ramp_up_slide =', time_ramp_up_slide
#endif
#if (SLIDE_LAW==2)
write(10, fmt=trim(fmt3)) 'red_pres_limit_fact =', red_pres_limit_fact
#endif
#if (BASAL_HYDROLOGY==1)
write(10, fmt=trim(fmt3)) 'c_Hw_slide  =', c_hw_slide
write(10, fmt=trim(fmt3)) 'Hw0_slide   =', hw0_slide
#endif
write(10, fmt=trim(fmt1)) ' '

write(10, fmt=trim(fmt2a)) 'Q_GEO_MOD = ', q_geo_mod
#if (Q_GEO_MOD==1)
write(10, fmt=trim(fmt3)) 'q_geo =', q_geo
#endif
write(10, fmt=trim(fmt1)) ' '

#if (defined(MARINE_ICE_BASAL_MELTING))
write(10, fmt=trim(fmt2)) 'MARINE_ICE_BASAL_MELTING =', marine_ice_basal_melting
#if (MARINE_ICE_BASAL_MELTING==2 || MARINE_ICE_BASAL_MELTING==3)
write(10, fmt=trim(fmt3)) 'qbm_marine               =', qbm_marine
#endif
write(10, fmt=trim(fmt1)) ' '
#endif

#if (MARGIN==3)
write(10, fmt=trim(fmt2)) 'FLOATING_ICE_BASAL_MELTING =', floating_ice_basal_melting
#if (FLOATING_ICE_BASAL_MELTING<=3)
write(10, fmt=trim(fmt3)) 'qbm_float_1 =', qbm_float_1
write(10, fmt=trim(fmt3)) 'qbm_float_2 =', qbm_float_2
#endif
write(10, fmt=trim(fmt3)) 'qbm_float_3 =', qbm_float_3
write(10, fmt=trim(fmt3)) 'z_abyss     =', z_abyss
#if (FLOATING_ICE_BASAL_MELTING==4)
write(10, fmt=trim(fmt3)) 'temp_ocean  =', temp_ocean
write(10, fmt=trim(fmt3)) 'Omega_qbm   =', omega_qbm
write(10, fmt=trim(fmt3)) 'alpha_qbm   =', alpha_qbm
#endif
#if (FLOATING_ICE_BASAL_MELTING==4 || FLOATING_ICE_BASAL_MELTING==5)
write(10, fmt=trim(fmt3)) 'H_w_0       =', h_w_0
#endif
write(10, fmt=trim(fmt1)) ' '
#endif

write(10, fmt=trim(fmt2)) 'REBOUND       =', rebound
#if (REBOUND==1)
write(10, fmt=trim(fmt3)) 'frac_llra     =', frac_llra
#endif
#if (REBOUND==1 || REBOUND==2)
write(10, fmt=trim(fmt2)) 'TIME_LAG_MOD  =', time_lag_mod
#if (TIME_LAG_MOD==1)
write(10, fmt=trim(fmt3)) 'time_lag      =', time_lag
#elif (TIME_LAG_MOD==2)
write(10, fmt=trim(fmt1)) 'time_lag_file = '//time_lag_file
#else
stop ' >>> sico_init: TIME_LAG_MOD must be either 1 or 2!'
#endif
#endif
#if (REBOUND==2)
write(10, fmt=trim(fmt2)) 'FLEX_RIG_MOD  =', flex_rig_mod
#if (FLEX_RIG_MOD==1)
write(10, fmt=trim(fmt3)) 'flex_rig      =', flex_rig
#elif (FLEX_RIG_MOD==2)
write(10, fmt=trim(fmt1)) 'flex_rig_file = '//flex_rig_file
#else
stop ' >>> sico_init: FLEX_RIG_MOD must be either 1 or 2!'
#endif
#endif
write(10, fmt=trim(fmt2)) 'Q_LITHO       =', q_litho
write(10, fmt=trim(fmt1)) ' '

#if (FLOW_LAW==2)
write(10, fmt=trim(fmt3)) 'gr_size   =', gr_size
write(10, fmt=trim(fmt1)) ' '
#endif
#if (FIN_VISC==2)
write(10, fmt=trim(fmt3)) 'sigma_res =', sigma_res
write(10, fmt=trim(fmt1)) ' '
#endif

write(10, fmt=trim(fmt2)) 'ENHMOD =', enhmod
#if (ENHMOD==1 || ENHMOD==2 || ENHMOD==3)
write(10, fmt=trim(fmt3)) 'enh_fact    =', enh_fact
#endif
#if (ENHMOD==2 || ENHMOD==3)
write(10, fmt=trim(fmt3)) 'enh_intg    =', enh_intg
#endif
#if (ENHMOD==2)
write(10, fmt=trim(fmt3)) 'age_trans   =', age_trans_0
#endif
#if (ENHMOD==3)
write(10, fmt=trim(fmt3)) 'date_trans1 =', date_trans1_0
write(10, fmt=trim(fmt3)) 'date_trans2 =', date_trans2_0
write(10, fmt=trim(fmt3)) 'date_trans3 =', date_trans3_0
#endif
#if (ENHMOD==4 || ENHMOD==5)
write(10, fmt=trim(fmt3)) 'enh_compr   =', enh_compr
write(10, fmt=trim(fmt3)) 'enh_shear   =', enh_shear
#endif
#if ((ENHMOD==1 || ENHMOD==2 || ENHMOD==3 || ENHMOD==4) && MARGIN==3)
write(10, fmt=trim(fmt3)) 'enh_shelf   =', enh_shelf
#endif
write(10, fmt=trim(fmt1)) ' '

write(10, fmt=trim(fmt3)) 'numdiff_H_t =', numdiff_h_t
write(10, fmt=trim(fmt3)) 'tau_cts     =', tau_cts
write(10, fmt=trim(fmt3)) 'vh_max      =', vh_max
write(10, fmt=trim(fmt3)) 'hd_min      =', hd_min
write(10, fmt=trim(fmt3)) 'hd_max      =', hd_max
#if (defined(QBM_MIN))
write(10, fmt=trim(fmt3)) 'qbm_min     =', qbm_min
#elif (defined(QB_MIN)) /* obsolete */
write(10, fmt=trim(fmt3)) 'qb_min      =', qb_min
#endif
#if (defined(QBM_MAX))
write(10, fmt=trim(fmt3)) 'qbm_max     =', qbm_max
#elif (defined(QB_MAX)) /* obsolete */
write(10, fmt=trim(fmt3)) 'qb_max      =', qb_max
#endif
write(10, fmt=trim(fmt3)) 'age_min     =', age_min
write(10, fmt=trim(fmt3)) 'age_max     =', age_max
write(10, fmt=trim(fmt3)) 'mean_accum  =', mean_accum
#if (ADV_VERT==1)
write(10, fmt=trim(fmt3)) 'age_diff    =', agediff
#endif
write(10, fmt=trim(fmt1)) ' '

write(10, fmt=trim(fmt2a)) 'DYNAMICS = ', dynamics
#if ((DYNAMICS==1 && MARGIN==3) || DYNAMICS==2)
#if (defined(LIS_OPTS))
write(10, fmt=trim(fmt1)) 'lis_opts = '//lis_opts
#endif
#if (defined(N_ITER_SSA))
write(10, fmt=trim(fmt2a)) 'n_iter_ssa = ', n_iter_ssa
#endif
#if (defined(RELAX_FACT_SSA))
write(10, fmt=trim(fmt3)) 'relax_fact_ssa =', relax_fact_ssa
#endif
#endif
#if (DYNAMICS==2 && defined(RATIO_SL_THRESH))
write(10, fmt=trim(fmt3)) 'ratio_sl_thresh =', ratio_sl_thresh
#endif
#if (DYNAMICS==2 && defined(SSTA_SIA_WEIGH_FCT))
write(10, fmt=trim(fmt2a)) 'SSTA_SIA_WEIGH_FCT = ', ssta_sia_weigh_fct
#endif
write(10, fmt=trim(fmt1)) ' '

write(10, fmt=trim(fmt2)) 'CALCMOD    =', calcmod
#if (CALCMOD==-1 && defined(TEMP_CONST))
write(10, fmt=trim(fmt3)) 'TEMP_CONST =', temp_const
#endif
#if (CALCMOD==-1 && defined(AGE_CONST))
write(10, fmt=trim(fmt3)) 'AGE_CONST  =', age_const
#endif
#if (CALCMOD==1 && defined(CTS_MELTING_FREEZING))
write(10, fmt=trim(fmt2)) 'CTS_MELTING_FREEZING =', cts_melting_freezing
#endif
write(10, fmt=trim(fmt2)) 'FLOW_LAW   =', flow_law
write(10, fmt=trim(fmt2)) 'FIN_VISC   =', fin_visc
write(10, fmt=trim(fmt2)) 'MARGIN     =', margin
#if (MARGIN==2)
write(10, fmt=trim(fmt2)) 'MARINE_ICE_FORMATION =', marine_ice_formation
write(10, fmt=trim(fmt2)) 'MARINE_ICE_CALVING   =', marine_ice_calving
#elif (MARGIN==3)
write(10, fmt=trim(fmt2)) 'ICE_SHELF_CALVING =', ice_shelf_calving
#endif
#if (defined(THK_EVOL))
write(10, fmt=trim(fmt2)) 'THK_EVOL   =', thk_evol
#else
stop ' >>> sico_init: Define THK_EVOL in header file!'
#endif
#if (defined(CALCTHK))
write(10, fmt=trim(fmt2)) 'CALCTHK    =', calcthk
#else
stop ' >>> sico_init: Define CALCTHK in header file!'
#endif
write(10, fmt=trim(fmt2)) 'ADV_HOR    =', adv_hor
write(10, fmt=trim(fmt2)) 'ADV_VERT   =', adv_vert
write(10, fmt=trim(fmt2)) 'TOPOGRAD   =', topograd
write(10, fmt=trim(fmt2)) 'TSURFACE   =', tsurface
#if (defined(MB_ACCOUNT))
write(10, fmt=trim(fmt2a)) 'MB_ACCOUNT = ', mb_account
#endif
write(10, fmt=trim(fmt1)) ' '

#if (defined(OUT_TIMES))
write(10, fmt=trim(fmt2a)) 'OUT_TIMES   = ', out_times
#endif
#if (defined(OUTPUT_INIT))
write(10, fmt=trim(fmt2a)) 'OUTPUT_INIT = ', output_init
#endif
write(10, fmt=trim(fmt2a)) 'OUTPUT      = ', output
#if (OUTPUT==1 || OUTPUT==3)
write(10, fmt=trim(fmt3))  'dtime_out   =' , dtime_out0
#endif
write(10, fmt=trim(fmt3))  'dtime_ser   =' , dtime_ser0
#if (OUTPUT==1 || OUTPUT==2)
write(10, fmt=trim(fmt2a)) 'ERGDAT      = ', ergdat
#endif
write(10, fmt=trim(fmt2a)) 'NETCDF      = ', netcdf
#if (OUTPUT==2 || OUTPUT==3)
write(10, fmt=trim(fmt2a)) 'n_output    = ', n_output
do n=1, n_output
   write(10, fmt=trim(fmt3))  'time_output =' , time_output0(n)
end do
#endif
write(10, fmt=trim(fmt1)) ' '

write(10, fmt=trim(fmt1)) 'Program version and date: '//version//' / '//date

close(10, status='keep')

!-------- Conversion of time quantities --------

year_zero  = year_zero*year_sec     ! a --> s
time_init  = time_init0*year_sec    ! a --> s
time_end   = time_end0*year_sec     ! a --> s
dtime      = dtime0*year_sec        ! a --> s
dtime_temp = dtime_temp0*year_sec   ! a --> s
dtime_ser  = dtime_ser0*year_sec    ! a --> s
#if (OUTPUT==1 || OUTPUT==3)
dtime_out  = dtime_out0*year_sec    ! a --> s
#endif
#if (OUTPUT==2 || OUTPUT==3)
do n=1, n_output
   time_output(n) = time_output0(n)*year_sec  ! a --> s
end do
#endif

if (.not.approx_integer_multiple(dtime_temp, dtime, eps_sp_dp)) &
   stop ' >>> sico_init: dtime_temp must be a multiple of dtime!'

if (.not.approx_integer_multiple(dtime_ser, dtime, eps_sp_dp)) &
   stop ' >>> sico_init: dtime_ser must be a multiple of dtime!'

#if (OUTPUT==1 || OUTPUT==3)
if (.not.approx_integer_multiple(dtime_out, dtime, eps_sp_dp)) &
   stop ' >>> sico_init: dtime_out must be a multiple of dtime!'
#endif

time = time_init

!-------- Mean accumulation --------

mean_accum = mean_accum*(1.0e-03_dp/year_sec)*(rho_w/rho)
!                      ! mm/a water equiv. --> m/s ice equiv.

!-------- Read data for delta_ts --------

#if (TSURFACE==4)

filename_with_path = trim(inpath)//'/general/'//trim(grip_temp_file)

open(21, iostat=ios, file=trim(filename_with_path), status='old')

if (ios /= 0) stop ' >>> sico_init: Error when opening the data file for delta_ts!'

read(21, fmt=*) ch_dummy, grip_time_min, grip_time_stp, grip_time_max

if (ch_dummy /= '#') then
   write(6, fmt=*) ' >>> sico_init: grip_time_min, grip_time_stp, grip_time_max'
   write(6, fmt=*) '                not defined in data file for delta_ts!'
   stop
end if

ndata_grip = (grip_time_max-grip_time_min)/grip_time_stp

allocate(griptemp(0:ndata_grip))

do n=0, ndata_grip
   read(21, fmt=*) d_dummy, griptemp(n)
end do

close(21, status='keep')

#endif

!-------- Read data for z_sl --------

#if (SEA_LEVEL==3)

filename_with_path = trim(inpath)//'/general/'//trim(sea_level_file)

open(21, iostat=ios, file=trim(filename_with_path), status='old')

if (ios /= 0) stop ' >>> sico_init: Error when opening the data file for z_sl!'

read(21, fmt=*) ch_dummy, specmap_time_min, specmap_time_stp, specmap_time_max

if (ch_dummy /= '#') then
   write(6, fmt=*) ' >>> sico_init: specmap_time_min, specmap_time_stp, specmap_time_max'
   write(6, fmt=*) '                not defined in data file for z_sl!'
   stop
end if

ndata_specmap = (specmap_time_max-specmap_time_min)/specmap_time_stp

allocate(specmap_zsl(0:ndata_specmap))

do n=0, ndata_specmap
   read(21, fmt=*) d_dummy, specmap_zsl(n)
end do

close(21, status='keep')

#endif

!-------- Determination of the geothermal heat flux --------

#if (Q_GEO_MOD==1)

!  ------ Constant value

do i=0, imax
do j=0, jmax
   q_geo(j,i) = q_geo *1.0e-03_dp   ! mW/m2 -> W/m2
end do
end do

#elif (Q_GEO_MOD==2)

stop ' >>> sico_init: Option Q_GEO_MOD==2 not available for this application!'

#endif

!-------- Determination of the time lag
!                              of the relaxing asthenosphere --------

#if (REBOUND==1 || REBOUND==2)

#if (TIME_LAG_MOD==1)

time_lag_asth = time_lag*year_sec   ! a -> s

#elif (TIME_LAG_MOD==2)

filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
                     trim(time_lag_file)

open(29, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')

if (ios /= 0) stop ' >>> sico_init: Error when opening the time-lag file!'

do n=1, 6; read(29, fmt='(a)') ch_dummy; end do

do j=jmax, 0, -1
   read(29, fmt=*) (time_lag_asth(j,i), i=0,imax)
end do

close(29, status='keep')

time_lag_asth = time_lag_asth*year_sec   ! a -> s

#endif

#elif (REBOUND==0)

time_lag_asth = 0.0_dp   ! dummy values

#endif

!-------- Determination of the flexural rigidity of the lithosphere --------

#if (REBOUND==2)

#if (FLEX_RIG_MOD==1)

flex_rig_lith = flex_rig

#elif (FLEX_RIG_MOD==2)

filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
                     trim(flex_rig_file)

open(29, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')

if (ios /= 0) stop ' >>> sico_init: Error when opening the flex-rig file!'

do n=1, 6; read(29, fmt='(a)') ch_dummy; end do

do j=jmax, 0, -1
   read(29, fmt=*) (flex_rig_lith(j,i), i=0,imax)
end do

close(29, status='keep')

#endif

#elif (REBOUND==0 || REBOUND==1)

flex_rig_lith = 0.0_dp   ! dummy values

#endif

!-------- Definition of initial values --------

!  ------ Present topography

#if (ANF_DAT==1)

stop ' >>> sico_init: ANF_DAT==1 not allowed for this application!'

!  ------ Ice-free, relaxed bedrock

#elif (ANF_DAT==2)

call topography2(dxi, deta)

z_sl = -1.11e+11_dp   ! dummy value for initial call of subroutine boundary

call boundary(time_init, dtime, dxi, deta, &
              delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)

as_perp_apl = 0.0_dp

q_bm    = 0.0_dp
q_tld   = 0.0_dp
q_b_tot = 0.0_dp 

p_b_w   = 0.0_dp
h_w     = 0.0_dp

call init_temp_water_age_2()

#if (ENHMOD==1)
   call calc_enhance_1()
#elif (ENHMOD==2)
   call calc_enhance_2()
#elif (ENHMOD==3)
   call calc_enhance_3(time_init)
#elif (ENHMOD==4)
   call calc_enhance_4()
#elif (ENHMOD==5)
   call calc_enhance_5()
#else
   stop ' >>> sico_init: Parameter ENHMOD must be between 1 and 5!'
#endif

!  ------ Read initial state from output of previous simulation

#elif (ANF_DAT==3)

call topography3(dxi, deta, z_sl, anfdatname)

call boundary(time_init, dtime, dxi, deta, &
              delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)

where ((maske==0_i2b).or.(maske==3_i2b))
                 ! grounded or floating ice
   as_perp_apl = as_perp
elsewhere        ! maske==1_i2b or 2_i2b, ice-free land or sea
   as_perp_apl = 0.0_dp
end where

q_b_tot = q_bm + q_tld

#endif

!-------- Inner-point flag --------

flag_inner_point = .true.

flag_inner_point(0,:)    = .false.
flag_inner_point(jmax,:) = .false.

flag_inner_point(:,0)    = .false.
flag_inner_point(:,imax) = .false.

!-------- Grounding line and calving front flags --------

flag_grounding_line_1 = .false.
flag_grounding_line_2 = .false.

flag_calving_front_1  = .false.
flag_calving_front_2  = .false.

#if (MARGIN==1 || MARGIN==2)

continue

#elif (MARGIN==3)

do i=1, imax-1
do j=1, jmax-1

   if ( (maske(j,i)==0_i2b) &   ! grounded ice
        .and. &
          (    (maske(j,i+1)==3_i2b)   &   ! with
           .or.(maske(j,i-1)==3_i2b)   &   ! one
           .or.(maske(j+1,i)==3_i2b)   &   ! neighbouring
           .or.(maske(j-1,i)==3_i2b) ) &   ! floating ice point
      ) &
   flag_grounding_line_1(j,i) = .true.

   if ( (maske(j,i)==3_i2b) &   ! floating ice
        .and. &
          (    (maske(j,i+1)==0_i2b)   &   ! with
           .or.(maske(j,i-1)==0_i2b)   &   ! one
           .or.(maske(j+1,i)==0_i2b)   &   ! neighbouring
           .or.(maske(j-1,i)==0_i2b) ) &   ! grounded ice point
      ) &
   flag_grounding_line_2(j,i) = .true.

   if ( (maske(j,i)==3_i2b) &   ! floating ice
        .and. &
          (    (maske(j,i+1)==2_i2b)   &   ! with
           .or.(maske(j,i-1)==2_i2b)   &   ! one
           .or.(maske(j+1,i)==2_i2b)   &   ! neighbouring
           .or.(maske(j-1,i)==2_i2b) ) &   ! ocean point
      ) &
   flag_calving_front_1(j,i) = .true.

   if ( (maske(j,i)==2_i2b) &   ! ocean
        .and. &
          (    (maske(j,i+1)==3_i2b)   &   ! with
           .or.(maske(j,i-1)==3_i2b)   &   ! one
           .or.(maske(j+1,i)==3_i2b)   &   ! neighbouring
           .or.(maske(j-1,i)==3_i2b) ) &   ! floating ice point
      ) &
   flag_calving_front_2(j,i) = .true.

end do
end do

do i=1, imax-1

   j=0
   if ( (maske(j,i)==2_i2b) &   ! ocean
        .and. (maske(j+1,i)==3_i2b) &   ! with one neighbouring
      ) &                               ! floating ice point
   flag_calving_front_2(j,i) = .true.

   j=jmax
   if ( (maske(j,i)==2_i2b) &   ! ocean
        .and. (maske(j-1,i)==3_i2b) &   ! with one neighbouring
      ) &                               ! floating ice point
   flag_calving_front_2(j,i) = .true.

end do

do j=1, jmax-1

   i=0
   if ( (maske(j,i)==2_i2b) &   ! ocean
        .and. (maske(j,i+1)==3_i2b) &   ! with one neighbouring
      ) &                               ! floating ice point
   flag_calving_front_2(j,i) = .true.

   i=imax
   if ( (maske(j,i)==2_i2b) &   ! ocean
        .and. (maske(j,i-1)==3_i2b) &   ! with one neighbouring
      ) &                               ! floating ice point
   flag_calving_front_2(j,i) = .true.

end do

#else
stop ' >>> sico_init: MARGIN must be either 1, 2 or 3!'
#endif

!-------- Distance between grid points with delta_i=ir, delta_j=jr --------

#if (GRID==0 || GRID==1)

do ir=-imax, imax
do jr=-jmax, jmax
   dist_dxdy(jr,ir) = sqrt( (real(ir,dp)*dxi)**2 + (real(jr,dp)*deta)**2 )
                  ! distortion due to stereographic projection not accounted for
end do
end do

#endif

!-------- Initial velocities --------

call calc_temp_melt()

#if (DYNAMICS==1 || DYNAMICS==2)

call calc_vxy_b_sia(time, z_sl)
call calc_vxy_sia(dzeta_c, dzeta_t)

#if (MARGIN==3 || DYNAMICS==2)
call calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t, ii, jj, nn)
#endif

call calc_vz_grounded(dxi, deta, dzeta_c, dzeta_t)

#if (MARGIN==3)
call calc_vz_floating(z_sl, dxi, deta, dzeta_c)
#endif

#elif (DYNAMICS==0)

call calc_vxy_static()
call calc_vz_static()

#else
   stop ' >>> sico_init: DYNAMICS must be either 0, 1 or 2!'
#endif

call calc_dxyz(dxi, deta, dzeta_c, dzeta_t)

!-------- Initial enthalpies --------

#if (CALCMOD==0 || CALCMOD==-1)

do i=0, imax
do j=0, jmax

   do kc=0, kcmax
      enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
   end do

   do kt=0, ktmax
      enth_t(kt,j,i) = enth_c(0,j,i)
   end do

end do
end do

#elif (CALCMOD==1)

do i=0, imax
do j=0, jmax

   do kc=0, kcmax
      enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
   end do

   if ( (maske(j,i) == 0_i2b).and.(n_cts(j,i) == 1_i2b) ) then
      do kt=0, ktmax
         enth_t(kt,j,i) = enth_fct_temp_omega(temp_t_m(kt,j,i), omega_t(kt,j,i))
      end do
   else
      do kt=0, ktmax
         enth_t(kt,j,i) = enth_c(0,j,i)
      end do
   end if

end do
end do

#elif (CALCMOD==2 || CALCMOD==3)

do i=0, imax
do j=0, jmax

   do kc=0, kcmax
      enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), omega_c(kc,j,i))
   end do

   do kt=0, ktmax
      enth_t(kt,j,i) = enth_c(0,j,i)
   end do

end do
end do

#else

stop ' >>> sico_init: Parameter CALCMOD must be either -1, 0, 1, 2 or 3!'

#endif

!-------- Initialize time-series files --------

!  ------ Time-series file for the ice sheet on the whole

filename_with_path = trim(outpath)//'/'//trim(runname)//'.ser'

open(12, iostat=ios, file=trim(filename_with_path), status='new')

if (ios /= 0) stop ' >>> sico_init: Error when opening the ser file!'

write(12,1102)
write(12,1103)

   1102 format('         t(a)  D_Ts(deg C)      z_sl(m)',/, &
               '                    V(m^3)     V_g(m^3)     V_f(m^3)', &
               '       A(m^2)     A_g(m^2)     A_f(m^2)',/, &
               '                               V_sle(m)     V_t(m^3)', &
               '     A_t(m^2)',/, &
               '                               H_max(m)   H_t_max(m)', &
               '    zs_max(m)  vs_max(m/a)   Tbh_max(C)')
   1103 format('----------------------------------------------------', &
               '---------------------------------------')

!  ------ Time-series file for selected positions ("deep boreholes")

n_core = 2   ! central dome, position halfway to coast

allocate(lambda_core(n_core), phi_core(n_core), &
         x_core(n_core), y_core(n_core), ch_core(n_core))

ch_core(1)     = 'P1'
lambda_core(1) =    0.0_dp  ! dummy
phi_core(1)    =    0.0_dp  ! dummy
x_core(1)      =  750.0_dp *1.0e+03_dp    ! Position of the central dome
y_core(1)      =  750.0_dp *1.0e+03_dp    ! (750 km, 750 km),
                                          ! conversion km -> m

ch_core(2)     = 'P2'
lambda_core(2) =    0.0_dp  ! dummy
phi_core(2)    =    0.0_dp  ! dummy
x_core(2)      =  750.0_dp *1.0e+03_dp    ! Position halfway to the coast
y_core(2)      = 1125.0_dp *1.0e+03_dp    ! (750 km, 1125 km),
                                          ! conversion km -> m

filename_with_path = trim(outpath)//'/'//trim(runname)//'.core'

open(14, iostat=ios, file=trim(filename_with_path), status='new')

if (ios /= 0) stop ' >>> sico_init: Error when opening the core file!'

if (forcing_flag == 1) then

   write(14,1106)
   write(14,1107)

   1106 format('         t(a)      D_Ts(C)      z_sl(m)',/, &
               '                   H_P1(m)      H_P2(m)',/, &
               '                 v_P1(m/a)    v_P2(m/a)',/, &
               '                   T_P1(C)      T_P2(C)')
   1107 format('---------------------------------------')

end if

!-------- Output of the initial state --------

#if (defined(OUTPUT_INIT))

#if (OUTPUT_INIT==0)
   flag_init_output = .false.
#elif (OUTPUT_INIT==1)
   flag_init_output = .true.
#else
   stop ' >>> sico_init: OUTPUT_INIT must be either 0 or 1!'
#endif

#else
   flag_init_output = .true.   ! default for undefined parameter OUTPUT_INIT
#endif

#if (OUTPUT==1)

#if (ERGDAT==0)
   flag_3d_output = .false.
#elif (ERGDAT==1)
   flag_3d_output = .true.
#else
   stop ' >>> sico_init: ERGDAT must be either 0 or 1!'
#endif

   if (flag_init_output) &
      call output1(runname, time_init, delta_ts, glac_index, z_sl, &
                   flag_3d_output, ndat2d, ndat3d)

#elif (OUTPUT==2)

if (time_output(1) <= time_init+eps) then

#if (ERGDAT==0)
   flag_3d_output = .false.
#elif (ERGDAT==1)
   flag_3d_output = .true.
#else
   stop ' >>> sico_init: ERGDAT must be either 0 or 1!'
#endif

   call output1(runname, time_init, delta_ts, glac_index, z_sl, &
                flag_3d_output, ndat2d, ndat3d)

end if

#elif (OUTPUT==3)

   flag_3d_output = .false.

   if (flag_init_output) &
      call output1(runname, time_init, delta_ts, glac_index, z_sl, &
                   flag_3d_output, ndat2d, ndat3d)

if (time_output(1) <= time_init+eps) then

   flag_3d_output = .true.

   call output1(runname, time_init, delta_ts, glac_index, z_sl, &
                flag_3d_output, ndat2d, ndat3d)

end if

#else
   stop ' >>> sico_init: OUTPUT must be either 1, 2 or 3!'
#endif

if (flag_init_output) then
   call output2(time_init, dxi, deta, delta_ts, glac_index, z_sl)
   call output4(time_init, dxi, deta, delta_ts, glac_index, z_sl)
end if

end subroutine sico_init

!-------------------------------------------------------------------------------
!> Definition of the initial surface and bedrock topography
!! (including gradients) and of the horizontal grid spacings dxi, deta.
!! For ice-free initial topography with relaxed lithosphere
!! (not defined for this domain, thus routine stops execution of SICOPOLIS).
!<------------------------------------------------------------------------------
subroutine topography1(dxi, deta)

implicit none

real(dp), intent(out) :: dxi, deta

! integer(i4b) :: i, j
! integer(i4b) :: ios
! real(dp)     :: xi0, eta0

dxi=0.0_dp; deta=0.0_dp   ! dummy values

write(6,'(/1x,a)') '>>> topography1: ANF_DAT==1 not defined for this domain!'
stop

end subroutine topography1

!-------------------------------------------------------------------------------
!> Definition of the initial surface and bedrock topography
!! (including gradients) and of the horizontal grid spacings dxi, deta.
!! For ice-free initial topography with relaxed lithosphere.
!<------------------------------------------------------------------------------
subroutine topography2(dxi, deta)

#if (GRID==0 || GRID==1)
  use stereo_proj_m, only : stereo_inv_ellipsoid
#endif

  use metric_m
  use topograd_m

implicit none

real(dp), intent(out) :: dxi, deta

integer(i4b) :: i, j, n
integer(i4b) :: ios
real(dp)     :: xi0, eta0

real(dp)                           :: half_width
real(dp), dimension(0:JMAX,0:IMAX) :: zl0_update

real(dp), parameter :: zl0_ocean = -2000.0_dp
real(dp), parameter :: epss      =     0.01_dp

!-------- Set topography --------

#if (!defined(LAND_OCEAN_TRANSITION_WIDTH))

write(6, fmt='(a)') ' >>> topography2: LAND_OCEAN_TRANSITION_WIDTH not defined.'
write(6, fmt='(a)') '                  Value 0 assumed by default!'
write(6, fmt='(a)') ' '

half_width = 0.0_dp

#else

half_width = 0.5_dp*abs(land_ocean_transition_width) *1000.0_dp   ! km -> m

#endif

zl0 = 10.0_dp

zl0(0,:)    = zl0_ocean
zl0(jmax,:) = zl0_ocean

zl0(:,0)    = zl0_ocean
zl0(:,imax) = zl0_ocean

maske = 1

maske(0,:)    = 2
maske(jmax,:) = 2

maske(:,0)    = 2
maske(:,imax) = 2

!-------- Further stuff --------

#if (GRID==0 || GRID==1)

dxi  = dx *1000.0_dp   ! km -> m
deta = dx *1000.0_dp   ! km -> m

xi0  = x0 *1000.0_dp   ! km -> m
eta0 = y0 *1000.0_dp   ! km -> m

#elif (GRID==2)

stop ' >>> topography2: GRID==2 not allowed for this application!'

#endif

do i=0, imax
   xi(i)  = xi0  + real(i,dp)*dxi
end do

do j=0, jmax
   eta(j) = eta0 + real(j,dp)*deta
end do

do i=1, imax-1
do j=1, jmax-1

   if ( (xi(i) >= (1150.0e+03_dp+half_width-epss)) &
        .and.(eta(j) >= (625.0e+03_dp+half_width-epss)) &
        .and.(eta(j) <= (875.0e+03_dp-half_width+epss)) ) then

      maske(j,i) = 2
      zl0(j,i)   = zl0_ocean

   end if

end do
end do

if (half_width > epss) then

   zl0_update = zl0

   do n=1, 1000

      do i=1, imax-1
      do j=1, jmax-1

         if ( ( .not.( (xi(i) >= (1150.0e+03_dp+half_width-epss)) &
                       .and.(eta(j) >= (625.0e+03_dp+half_width-epss)) &
                       .and.(eta(j) <= (875.0e+03_dp-half_width+epss)) ) ) &
              .and. &
              ( (xi(i) >= (1150.0e+03_dp-half_width-epss)) &
                .and.(eta(j) >= (625.0e+03_dp-half_width-epss)) &
                .and.(eta(j) <= (875.0e+03_dp+half_width+epss)) ) ) then

            maske(j,i) = 2

            zl0_update(j,i) = zl0(j,i) &
              + 0.1_dp &
                * (zl0(j,i+1)-2.0_dp*zl0(j,i)+zl0(j,i-1)) &
              + 0.1_dp &
                * (zl0(j+1,i)-2.0_dp*zl0(j,i)+zl0(j-1,i))
              ! Solving the diffusion equation

         end if

      end do
      end do

      zl0 = zl0_update

   end do

end if

do i=0, imax
do j=0, jmax

   if (maske(j,i) <= 1) then
      maske(j,i) = 1
      zs(j,i) = zl0(j,i)
      zb(j,i) = zl0(j,i)
      zl(j,i) = zl0(j,i)
   else   ! (maske(j,i) >= 2)
      maske(j,i) = 2
#if (MARGIN==1 || MARGIN==2)
      zs(j,i) = zl0(j,i)
      zb(j,i) = zl0(j,i)
#elif (MARGIN==3)
      zs(j,i) = 0.0_dp   ! present-day
      zb(j,i) = 0.0_dp   ! sea level
#endif
      zl(j,i) = zl0(j,i)
   end if

   zm(j,i) = zb(j,i)
   n_cts(j,i) = -1
   kc_cts(j,i) = 0

   h_c(j,i) = 0.0_dp
   h_t(j,i) = 0.0_dp

   dzs_dtau(j,i)  = 0.0_dp
   dzm_dtau(j,i)  = 0.0_dp
   dzb_dtau(j,i)  = 0.0_dp
   dzl_dtau(j,i)  = 0.0_dp
   dh_c_dtau(j,i) = 0.0_dp
   dh_t_dtau(j,i) = 0.0_dp

end do
end do

maske_old = maske

!-------- Geographic coordinates, metric tensor,
!                                 gradients of the topography --------

do i=0, imax
do j=0, jmax

#if (GRID==0 || GRID==1)   /* Stereographic projection */

   call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
                             lambda0, phi0, lambda(j,i), phi(j,i))

#elif (GRID==2)   /* Geographic coordinates */

   lambda(j,i) = xi(i)
   phi(j,i)    = eta(j)

#endif

end do
end do

call metric()

#if (TOPOGRAD==0)
call topograd_1(dxi, deta, 1)
#elif (TOPOGRAD==1)
call topograd_2(dxi, deta, 1)
#endif

!-------- Corresponding area of grid points --------

do i=0, imax
do j=0, jmax
   area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
end do
end do

end subroutine topography2

!-------------------------------------------------------------------------------
!> Definition of the initial surface and bedrock topography
!! (including gradients) and of the horizontal grid spacings dxi, deta.
!! For initial topography from previous simulation.
!<------------------------------------------------------------------------------
subroutine topography3(dxi, deta, z_sl, anfdatname)

  use read_m, only : read_erg_nc

#if (GRID==0 || GRID==1)
  use stereo_proj_m, only : stereo_inv_ellipsoid
#endif

  use metric_m
  use topograd_m

implicit none

character(len=100), intent(in) :: anfdatname

real(dp),          intent(out) :: dxi, deta, z_sl

integer(i4b) :: i, j, n

real(dp)                           :: half_width
real(dp), dimension(0:JMAX,0:IMAX) :: zl0_update

real(dp), parameter :: zl0_ocean = -2000.0_dp
real(dp), parameter :: epss      =     0.01_dp

!-------- Read data from time-slice file of previous simulation --------

call read_erg_nc(z_sl, anfdatname)

!-------- Set topography of the relaxed bedrock --------

#if (!defined(LAND_OCEAN_TRANSITION_WIDTH))

write(6, fmt='(a)') ' >>> topography3: LAND_OCEAN_TRANSITION_WIDTH not defined.'
write(6, fmt='(a)') '                  Value 0 assumed by default!'
write(6, fmt='(a)') ' '

half_width = 0.0_dp

#else

half_width = 0.5_dp*abs(land_ocean_transition_width) *1000.0_dp   ! km -> m

#endif

zl0 = 10.0_dp

zl0(0,:)    = zl0_ocean
zl0(jmax,:) = zl0_ocean

zl0(:,0)    = zl0_ocean
zl0(:,imax) = zl0_ocean

!-------- Further stuff --------

#if (GRID==0 || GRID==1)

dxi  = dx *1000.0_dp   ! km -> m
deta = dx *1000.0_dp   ! km -> m

#elif (GRID==2)

stop ' >>> topography3: GRID==2 not allowed for this application!'

#endif

do i=1, imax-1
do j=1, jmax-1

   if ( (xi(i) >= (1150.0e+03_dp+half_width-epss)) &
        .and.(eta(j) >= (625.0e+03_dp+half_width-epss)) &
        .and.(eta(j) <= (875.0e+03_dp-half_width+epss)) ) then

      zl0(j,i)   = zl0_ocean

   end if

end do
end do

if (half_width > epss) then

   zl0_update = zl0

   do n=1, 1000

      do i=1, imax-1
      do j=1, jmax-1

         if ( ( .not.( (xi(i) >= (1150.0e+03_dp+half_width-epss)) &
                       .and.(eta(j) >= (625.0e+03_dp+half_width-epss)) &
                       .and.(eta(j) <= (875.0e+03_dp-half_width+epss)) ) ) &
              .and. &
              ( (xi(i) >= (1150.0e+03_dp-half_width-epss)) &
                .and.(eta(j) >= (625.0e+03_dp-half_width-epss)) &
                .and.(eta(j) <= (875.0e+03_dp+half_width+epss)) ) ) then

            zl0_update(j,i) = zl0(j,i) &
              + 0.1_dp &
                * (zl0(j,i+1)-2.0_dp*zl0(j,i)+zl0(j,i-1)) &
              + 0.1_dp &
                * (zl0(j+1,i)-2.0_dp*zl0(j,i)+zl0(j-1,i))
              ! Solving the diffusion equation

         end if

      end do
      end do

      zl0 = zl0_update

   end do

end if

!-------- Geographic coordinates, metric tensor,
!                                 gradients of the topography --------

do i=0, imax
do j=0, jmax

#if (GRID==0 || GRID==1)   /* Stereographic projection */

   call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
                             lambda0, phi0, lambda(j,i), phi(j,i))

#elif (GRID==2)   /* Geographic coordinates */

   lambda(j,i) = xi(i)
   phi(j,i)    = eta(j)

#endif

end do
end do

call metric()

#if (TOPOGRAD==0)
call topograd_1(dxi, deta, 1)
#elif (TOPOGRAD==1)
call topograd_2(dxi, deta, 1)
#endif

!-------- Corresponding area of grid points --------

do i=0, imax
do j=0, jmax
   area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
end do
end do

end subroutine topography3

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

end module sico_init_m
!