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, read_kei

  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
integer(i4b)       :: ndata_insol
real(dp)           :: dtime0, dtime_temp0, dtime_wss0, dtime_out0, dtime_ser0
real(dp)           :: time_init0, time_end0, time_output0(100)
real(dp)           :: time_chasm_init0, time_chasm_end0
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)'

character(len=  8) :: ch_imax
character(len=128) :: fmt4

write(ch_imax, fmt='(i8)') imax
write(fmt4,    fmt='(a)')  '('//trim(adjustl(ch_imax))//'(i1),i1)'

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, rho_i, rho_c, kappa_c, c_c)

!  ------ Density of ice-dust mixture

rho = (1.0_dp-frac_dust)*rho_i + frac_dust*rho_c

rho_inv = 1.0_dp/rho

!  ------ 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 || GRID==1)

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

   if ((imax /= 90).or.(jmax /= 90)) then
      stop ' >>> sico_init: Wrong values for IMAX and JMAX!'
   end if

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

   if ((imax /= 180).or.(jmax /= 180)) then
      stop ' >>> sico_init: Wrong values for IMAX and JMAX!'
   end if

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

   if ((imax /= 360).or.(jmax /= 360)) then
      stop ' >>> sico_init: Wrong values for IMAX and JMAX!'
   end if

else
      stop ' >>> sico_init: Wrong value for DX!'
end if

#elif (GRID==2)

stop ' >>> sico_init: GRID==2 not allowed for nmars 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
#if (REBOUND==2)
dtime_wss0  = dtime_wss0
#endif

!-------- 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 (CHASM==2)
time_chasm_init0 = time_chasm_init0
time_chasm_end0  = time_chasm_end0
#endif

#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)) ' '

write(10, fmt=trim(fmt3)) 'x0      =', x0
write(10, fmt=trim(fmt3)) 'y0      =', y0
write(10, fmt=trim(fmt3)) 'dx      =', dx
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
#if (REBOUND==2)
write(10, fmt=trim(fmt3)) 'dtime_wss  =', dtime_wss0
#endif
write(10, fmt=trim(fmt1)) ' '

write(10, fmt=trim(fmt2)) 'ANF_DAT =', anf_dat
write(10, fmt=trim(fmt1)) 'zs_present file   = '//zs_present_file
#if (ANF_DAT==1)
write(10, fmt=trim(fmt1)) 'zl_present file   = '//zl_present_file
#endif
write(10, fmt=trim(fmt1)) 'zl0 file          = '//zl0_file
write(10, fmt=trim(fmt1)) 'mask_present file = '//mask_present_file
#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 (TSURFACE==1 || TSURFACE==6)
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
#endif
#if (TSURFACE==1 || TSURFACE==2 || TSURFACE==3)
write(10, fmt=trim(fmt3)) 'temp0_ma_90N    =', temp0_ma_90n
#endif
#if (TSURFACE==1 || TSURFACE==2 || TSURFACE==3 || TSURFACE==4 || TSURFACE==5)
write(10, fmt=trim(fmt3)) 'c_ma            =', c_ma
write(10, fmt=trim(fmt3)) 'gamma_ma        =', gamma_ma
#endif
#if (TSURFACE==5)
write(10, fmt=trim(fmt1)) 'Insolation file = '//insol_ma_90n_file
#endif
#if (TSURFACE==4 || TSURFACE==5 || TSURFACE==6)
write(10, fmt=trim(fmt3)) 'albedo          =', albedo
#endif

#if (ACC_PRESENT==1)
write(10, fmt=trim(fmt3)) 'acc_present_val =', acc_present_val
#endif

#if (ACCSURFACE==1 || ACCSURFACE==2)
write(10, fmt=trim(fmt3)) 'gamma_s         =', gamma_s
#endif

#if (ABLSURFACE==1 || ABLSURFACE==2)
write(10, fmt=trim(fmt3)) 'eld_0     =', eld_0
write(10, fmt=trim(fmt3)) 'g_mb_0    =', g_0
write(10, fmt=trim(fmt3)) 'gamma_eld =', gamma_eld
write(10, fmt=trim(fmt3)) 'gamma_g   =', gamma_g
#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
#elif (Q_GEO_MOD==2)
write(10, fmt=trim(fmt1)) 'q_geo file = '//q_geo_file
#endif
write(10, fmt=trim(fmt1)) ' '

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(fmt3)) 'frac_dust =', frac_dust
write(10, fmt=trim(fmt1)) ' '

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(fmt2a)) 'CHASM = ', chasm
#if (CHASM==2)
write(10, fmt=trim(fmt1)) 'mask_chasm file = '//mask_chasm_file
write(10, fmt=trim(fmt3)) 'time_chasm_init =', time_chasm_init0
write(10, fmt=trim(fmt3)) 'time_chasm_end  =', time_chasm_end0
write(10, fmt=trim(fmt3)) 'q_geo_chasm     =', q_geo_chasm
write(10, fmt=trim(fmt3)) 'erosion_chasm   =', erosion_chasm
write(10, fmt=trim(fmt1)) ' '
#endif

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
write(10, fmt=trim(fmt2)) 'ACC_UNIT   =', acc_unit
write(10, fmt=trim(fmt2)) 'ACC_PRESENT=', acc_present
write(10, fmt=trim(fmt2)) 'ACCSURFACE =', accsurface
write(10, fmt=trim(fmt2)) 'ABLSURFACE =', ablsurface
#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
#if (REBOUND==2)
dtime_wss  = dtime_wss0*year_sec    ! a --> s
#endif
#if (CHASM==2)
time_chasm_init = time_chasm_init0 *year_sec    ! a --> s
time_chasm_end  = time_chasm_end0  *year_sec    ! a --> s
#endif
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 (REBOUND==2)
if (.not.approx_integer_multiple(dtime_wss, dtime, eps_sp_dp)) &
   stop ' >>> sico_init: dtime_wss must be a multiple of dtime!'
#endif

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

!-------- Read accumulation measurements --------

#if (ACC_PRESENT==1)

do i=0, imax
do j=0, jmax
   accum_present(j,i) = acc_present_val
end do
end do

#endif

!  ------ Conversion mm/a water or (ice+dust) equivalent
!                          --> m/s (ice+dust) equivalent

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

#if (ACC_UNIT==1)

   accum_present(j,i) = accum_present(j,i) &
                           *(1.0e-03_dp/year_sec)*(rho_w/rho_i) &
                           *(1.0_dp/(1.0_dp-frac_dust))

#elif (ACC_UNIT==2)

   accum_present(j,i) = accum_present(j,i)*(1.0e-03_dp/year_sec)

#endif

end do
end do

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

#if (ACC_UNIT==1)

mean_accum = mean_accum*(1.0e-03_dp/year_sec)*(rho_w/rho_i) &
                       *(1.0_dp/(1.0_dp-frac_dust))
!                      ! mm/a water equiv. --> m/s (ice+dust) equiv.

#elif (ACC_UNIT==2)

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

#endif

!-------- Read chasm mask --------

#if (CHASM==2)

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

open(24, iostat=ios, file=trim(filename_with_path), recl=rcl2, status='old')

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

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

do j=jmax, 0, -1
   read(24, fmt=trim(fmt4)) (maske_chasm(j,i), i=0,imax)
end do

close(24, status='keep')

#endif

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

#if (SEA_LEVEL==3)

stop ' >>> sico_init: SEA_LEVEL==3 not allowed for nmars application!'

#endif

!-------- Reading of insolation data  --------

insol_ma_90 = 0.0_dp   ! Assignment of dummy values
obl_data    = 0.0_dp
ecc_data    = 0.0_dp
ave_data    = 0.0_dp
cp_data     = 0.0_dp

#if (TSURFACE==5 || TSURFACE==6)

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

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

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

read(21, fmt=*) ch_dummy, insol_time_min, insol_time_stp, insol_time_max

if (ch_dummy /= '#') then
   write(6, fmt=*) ' >>> sico_init: insol_time_min, insol_time_stp, insol_time_max'
   write(6, fmt=*) '                not defined in insolation-data file!'
   stop
end if

ndata_insol = (insol_time_max-insol_time_min)/insol_time_stp

if (ndata_insol > 100000) &
   stop ' >>> sico_init: Too many data in insolation-data file!'

do n=0, ndata_insol
   read(21, fmt=*) d_dummy, ecc_data(n), obl_data(n), cp_data(n), ave_data(n), insol_ma_90(n)
   obl_data(n) = obl_data(n) *pi_180
   ave_data(n) = ave_data(n) *pi_180
end do

close(21, status='keep')

#endif

!-------- Determination of the normal geothermal heat flux
!         (without the possible contribution from an active chasm area) --------

#if (Q_GEO_MOD==1)

!  ------ Constant value

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

#elif (Q_GEO_MOD==2)

!  ------ Read data from file

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

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

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

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

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

close(21, status='keep')

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

#endif

!-------- Reading of tabulated kei function--------

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

kei = 0.0_dp; n_data_kei = 0; kei_r_max = 0.0_dp; kei_r_incr = 0.0_dp
                                                    ! dummy values
#elif (REBOUND==2)

call read_kei()

#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)

call topography1(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)

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_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

#if (!defined(TEMP_INIT) || TEMP_INIT==1)
  call init_temp_water_age_1_1()
#elif (TEMP_INIT==2)
  call init_temp_water_age_1_2()
#elif (TEMP_INIT==3)
  call init_temp_water_age_1_3()
#elif (TEMP_INIT==4)
  call init_temp_water_age_1_4()
#else
  write(6, fmt='(a)') ' >>> sico_init:'
  write(6, fmt='(a)') '     TEMP_INIT must be set to either 1, 2, 3 or 4!'
  stop
#endif

#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

!  ------ 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(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 = 3   ! Points NP (north pole), C1, C2 (in Chasma Borealis)

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

ch_core(1)     = 'North Pole'
lambda_core(1) =    0.0_dp  ! dummy
phi_core(1)    =    0.0_dp  ! dummy
x_core(1)      =    0.0_dp *1.0e+03_dp    ! North pole,
y_core(1)      =    0.0_dp *1.0e+03_dp    ! conversion km -> m

ch_core(2)     = 'C1'
lambda_core(2) =    0.0_dp  ! dummy
phi_core(2)    =    0.0_dp  ! dummy
x_core(2)      = -150.0_dp *1.0e+03_dp    ! Point C1,
y_core(2)      = -290.0_dp *1.0e+03_dp    ! conversion km -> m

ch_core(3)     = 'C2'
lambda_core(3) =    0.0_dp  ! dummy
phi_core(3)    =    0.0_dp  ! dummy
x_core(3)      = -300.0_dp *1.0e+03_dp    ! Point C2,
y_core(3)      = -280.0_dp *1.0e+03_dp    ! 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_NP(m)      H_C1(m)      H_C2(m)',/, &
               '                 v_NP(m/a)    v_C1(m/a)    v_C2(m/a)')
   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 present-day initial topography.
!<------------------------------------------------------------------------------
subroutine topography1(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, n_dummy
real(dp)     :: d_dummy
real(dp)     :: xi0, eta0
character    :: ch_dummy

character(len=  8) :: ch_imax
character(len=128) :: fmt4
character(len=256) :: filename_with_path

write(ch_imax, fmt='(i8)') imax
write(fmt4,    fmt='(a)')  '('//trim(adjustl(ch_imax))//'(i1),i1)'

!-------- Read topography --------

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

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

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

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

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

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

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

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

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

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

open(24, iostat=ios, file=trim(filename_with_path), recl=rcl2, status='old')

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

do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
do n=1, 6; read(22, fmt='(a)') ch_dummy; end do
do n=1, 6; read(23, fmt='(a)') ch_dummy; end do
do n=1, 6; read(24, fmt='(a)') ch_dummy; end do

do j=jmax, 0, -1
   read(21, fmt=*)          (zs(j,i),    i=0,imax)
   read(22, fmt=*)          (zl(j,i),    i=0,imax)
   read(23, fmt=*)          (zl0(j,i),   i=0,imax)
   read(24, fmt=trim(fmt4)) (maske(j,i), i=0,imax)
end do

close(21, status='keep')
close(22, status='keep')
close(23, status='keep')
close(24, status='keep')

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

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

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

   if (maske(j,i) <= 1) then
      zb(j,i) = zl(j,i)
   else   ! (maske(j,i)>=2)
      stop ' >>> topography1: maske(j,i)>=2 not allowed for initial topography!'
   end if

   zs(j,i)  = zs(j,i) *1000.0_dp
   zb(j,i)  = zb(j,i) *1000.0_dp    ! km --> m
   zl(j,i)  = zl(j,i) *1000.0_dp
   zl0(j,i) = zl0(j,i)*1000.0_dp

   xi(i)  = xi0  + real(i,dp)*dxi
   eta(j) = eta0 + real(j,dp)*deta

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

   h_c(j,i) = zs(j,i)-zm(j,i)
   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 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
character    :: ch_dummy

character(len=  8) :: ch_imax
character(len=128) :: fmt4
character(len=256) :: filename_with_path

write(ch_imax, fmt='(i8)') imax
write(fmt4,    fmt='(a)')  '('//trim(adjustl(ch_imax))//'(i1),i1)'

!-------- Read topography --------

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

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

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

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

open(24, iostat=ios, file=trim(filename_with_path), recl=rcl2, status='old')

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

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

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

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

do j=jmax, 0, -1
   read(24, fmt=trim(fmt4)) (maske(j,i), i=0,imax)
end do

close(23, status='keep')
close(24, status='keep')

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

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

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)
      stop ' >>> topography2: maske(j,i)>=2 not allowed for initial topography!'
   end if

   zs(j,i)  = zs(j,i) *1000.0_dp
   zb(j,i)  = zb(j,i) *1000.0_dp   ! km --> m
   zl(j,i)  = zl(j,i) *1000.0_dp
   zl0(j,i) = zl0(j,i)*1000.0_dp

   xi(i)  = xi0  + real(i,dp)*dxi
   eta(j) = eta0 + real(j,dp)*deta

   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
integer(i4b) :: ios
character(len=256) :: filename_with_path
character :: ch_dummy

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

call read_erg_nc(z_sl, anfdatname)

!-------- Read topography of the relaxed bedrock --------

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

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

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

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

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

close(23, status='keep')

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

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

zl0 = zl0 *1000.0_dp   ! km -> m

!-------- 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
!