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

#if (NETCDF > 1)
  use netcdf
  use nc_check_m
#endif

#if (DISC>0 && !defined(EXEC_MAKE_C_DIS_0))
  use discharge_workers_m, only: disc_param, disc_fields
#endif
#if (DISC>0 && defined(EXEC_MAKE_C_DIS_0))
  use discharge_workers_m, only: disc_param, disc_fields, calc_c_dis_0
#endif

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

  use read_m, only : read_phys_para, read_target_topo_nc, 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
real(dp)           :: dtime0, dtime_temp0, dtime_wss0, dtime_out0, dtime_ser0
real(dp)           :: time_init0, time_end0, time_output0(100)
real(dp), dimension(0:JMAX,0:IMAX) :: precip_factor
real(dp)           :: d_dummy
character(len=100) :: anfdatname, target_topo_dat_name
character(len=256) :: filename_with_path
character(len=256) :: shell_command
character          :: ch_dummy
logical            :: flag_precip_monthly_mean
logical            :: flag_init_output, flag_3d_output

#if (defined(INITMIP_SMB_ANOM_FILE))
real(dp), dimension(0:IMAX,0:JMAX) :: as_anom_initmip_conv
#endif

#if (NETCDF > 1)
integer(i4b) :: dimid, ncid, ncv
!   dimid:       Dimension ID
!    ncid:       File ID
!     ncv:       Variable ID
#endif

real(dp), parameter :: inv_twelve = 1.0_dp/12.0_dp

character(len=64), parameter :: thisroutine = 'sico_init'

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)

!  ------ 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 /= 75).or.(jmax /= 140) ) 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 /= 280) ) then
      stop ' >>> sico_init: IMAX and/or JMAX wrong!'
   end if

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

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

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

#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 surface-temperature
!                         and precipitation switches --------

#if (TSURFACE == 5 && ACCSURFACE != 5)
stop ' >>> sico_init: Options TSURFACE==5 and ACCSURFACE==5 must be used together!'
#endif

#if (TSURFACE != 5 && ACCSURFACE == 5)
stop ' >>> sico_init: Options TSURFACE==5 and ACCSURFACE==5 must be used together!'
#endif

#if (TSURFACE == 6)
#if (ACCSURFACE != 6)
write(6, fmt='(a)') ' >>> sico_init: Options TSURFACE==6, ACCSURFACE==6'
write(6, fmt='(a)') '                and NETCDF==2 must be used together!'
stop
#elif (NETCDF != 2)
write(6, fmt='(a)') ' >>> sico_init: Options TSURFACE==6, ACCSURFACE==6'
write(6, fmt='(a)') '                and NETCDF==2 must be used together!'
stop
#endif
#endif

#if (ACCSURFACE == 6)
#if (TSURFACE != 6)
write(6, fmt='(a)') ' >>> sico_init: Options TSURFACE==6, ACCSURFACE==6'
write(6, fmt='(a)') '                and NETCDF==2 must be used together!'
stop
#elif (NETCDF != 2)
write(6, fmt='(a)') ' >>> sico_init: Options TSURFACE==6, ACCSURFACE==6'
write(6, fmt='(a)') '                and NETCDF==2 must be used together!'
stop
#endif
#endif

#if (defined(INITMIP_SMB_ANOM_FILE) && NETCDF != 2)
if ( trim(adjustl(initmip_smb_anom_file)) /= 'none' ) then
   write(6, fmt='(a)') ' >>> sico_init: Defining INITMIP_SMB_ANOM_FILE'
   write(6, fmt='(a)') '                must be used together with NETCDF==2!'
   stop
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


!-------- Check whether for the hidden ablation scheme
!         can be applied --------

#if (MARGIN==3 && MB_ACCOUNT==1)
write(6, fmt='(a)') &
   ' >>> sico_init: MB_ACCOUNT==1 does not yet work for ice shelves!'
stop
#endif

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

#if (TSURFACE <= 4)

forcing_flag = 1   ! forcing by delta_ts

#elif (TSURFACE == 5)

forcing_flag = 2   ! forcing by glac_index

#elif (TSURFACE == 6)

forcing_flag = 3   ! forcing by time-dependent surface temperature
                   ! and precipitation data

#endif

!-------- 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 (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
stop ' >>> sico_init: GRID==2 not allowed for this application!'
#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
#if (REBOUND==2)
write(10, fmt=trim(fmt3)) 'dtime_wss  =', dtime_wss0
#endif
#if (DISC>0)
write(10, fmt=trim(fmt3)) 'dtime_mar_coa =', dtime_mar_coa0
#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)
#if (defined(ZB_PRESENT_FILE))
write(10, fmt=trim(fmt1)) 'zb_present file   = '//zb_present_file
#endif
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)) ' '


#if (THK_EVOL==2)
write(10, fmt=trim(fmt3)) 'time_target_topo_init  =', time_target_topo_init0
write(10, fmt=trim(fmt3)) 'time_target_topo_final =', time_target_topo_final0
write(10, fmt=trim(fmt1)) 'Target-topography file = '//target_topo_dat_name
write(10, fmt=trim(fmt1)) 'Path to target-topography file = '//target_topo_dat_path
write(10, fmt=trim(fmt1)) ' '
#endif

#if (THK_EVOL==3)
write(10, fmt=trim(fmt3)) 'target_topo_tau  =', target_topo_tau0
write(10, fmt=trim(fmt1)) 'Target-topography file = '//target_topo_dat_name
write(10, fmt=trim(fmt1)) 'Path to target-topography file = '//target_topo_dat_path
write(10, fmt=trim(fmt1)) ' '
#endif

#if (THK_EVOL==4)
write(10, fmt=trim(fmt1)) 'Maximum ice extent mask file = '//mask_maxextent_file
write(10, fmt=trim(fmt1)) ' '
#endif

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)
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
#elif (TSURFACE==5)
write(10, fmt=trim(fmt1)) 'Glacial-index file = '//glac_ind_file
write(10, fmt=trim(fmt1)) 'temp_ma_anom file  = '//temp_ma_anom_file
write(10, fmt=trim(fmt3)) 'temp_ma_anom fact  = ', temp_ma_anom_fact
write(10, fmt=trim(fmt1)) 'temp_mj_anom file  = '//temp_mj_anom_file
write(10, fmt=trim(fmt3)) 'temp_mj_anom fact  = ', temp_mj_anom_fact
#elif (TSURFACE==6)
write(10, fmt=trim(fmt1)) 'temp_precip_anom file = '//temp_precip_anom_file
write(10, fmt=trim(fmt3)) 'temp_ma_anom fact     = ', temp_ma_anom_fact
write(10, fmt=trim(fmt3)) 'temp_mj_anom fact     = ', temp_mj_anom_fact
#endif

#if (!defined(PRECIP_PRESENT_FILE))
stop ' >>> sico_init: PRECIP_PRESENT_FILE not defined in the header file!'
#endif
#if (!defined(PRECIP_MA_PRESENT_FILE))
stop ' >>> sico_init: PRECIP_MA_PRESENT_FILE not defined in the header file!'
#endif
#if (defined(PRECIP_PRESENT_FILE) && defined(PRECIP_MA_PRESENT_FILE))
if ( trim(adjustl(precip_present_file)) /= 'none' ) then
   write(10, fmt=trim(fmt1)) 'precip_present_file    = '//precip_present_file
   flag_precip_monthly_mean = .true.
else if ( trim(adjustl(precip_ma_present_file)) /= 'none' ) then
   write(10, fmt=trim(fmt1)) 'precip_ma_present_file = '//precip_ma_present_file
   flag_precip_monthly_mean = .false.
else
   write(6, fmt='(a)') ' >>> sico_init: Neither PRECIP_PRESENT_FILE'
   write(6, fmt='(a)') '                nor PRECIP_MA_PRESENT_FILE'
   write(6, fmt='(a)') '                specified in the header file!'
   stop
end if
#endif
#if (!defined(PRECIP_ZS_REF_FILE))
stop ' >>> sico_init: PRECIP_ZS_REF_FILE not defined in the header file!'
#else
write(10, fmt=trim(fmt1)) 'precip_zs_ref_file     = '//precip_zs_ref_file
#endif
#if (defined(PRECIP_FACTOR_FILE))
if ( trim(adjustl(precip_factor_file)) /= 'none' ) &
   write(10, fmt=trim(fmt1)) 'precip_factor_file     = '//precip_factor_file
#endif
#if (ACCSURFACE==1)
write(10, fmt=trim(fmt3)) 'accfact        =', accfact
#elif (ACCSURFACE==2 || ACCSURFACE==3)
write(10, fmt=trim(fmt3)) 'gamma_s        =', gamma_s
#elif (ACCSURFACE==5)
write(10, fmt=trim(fmt1)) 'precip_anom file    = '//precip_anom_file
write(10, fmt=trim(fmt3)) 'precip_anom fact    = ', precip_anom_fact
#elif (ACCSURFACE==6)
write(10, fmt=trim(fmt1)) 'temp_precip_anom file = '//temp_precip_anom_file
write(10, fmt=trim(fmt3)) 'precip_anom fact      = ', precip_anom_fact
#endif
#if (ACCSURFACE <= 3)
write(10, fmt=trim(fmt2)) 'ELEV_DESERT =', elev_desert
#if (ELEV_DESERT == 1)
write(10, fmt=trim(fmt3)) 'gamma_p     =', gamma_p
write(10, fmt=trim(fmt3)) 'zs_thresh   =', zs_thresh
#endif
#endif

#if (ABLSURFACE==3)
write(10, fmt=trim(fmt3)) 'lambda_lti     =', lambda_lti
write(10, fmt=trim(fmt3)) 'temp_lti       =', temp_lti
#endif

write(10, fmt=trim(fmt1)) ' '

#if (defined(INITMIP_CONST_SMB))
write(10, fmt=trim(fmt2a)) 'INITMIP_CONST_SMB = ', initmip_const_smb
#endif
#if (defined(INITMIP_SMB_ANOM_FILE))
if ( trim(adjustl(initmip_smb_anom_file)) /= 'none' ) then
   write(10, fmt=trim(fmt1)) 'initmip_smb_anom file = '//initmip_smb_anom_file
end if
#endif
#if (defined(INITMIP_CONST_SMB) || defined(INITMIP_SMB_ANOM_FILE))
write(10, fmt=trim(fmt1)) ' '
#endif

write(10, fmt=trim(fmt2a)) 'DISC = ', disc
#if (DISC>0)
write(10, fmt=trim(fmt3)) 'c_dis_0   =', c_dis_0
write(10, fmt=trim(fmt3)) 'c_dis_fac =', c_dis_fac
write(10, fmt=trim(fmt3)) 'm_H       =', m_h
write(10, fmt=trim(fmt3)) 'm_D       =', m_d
write(10, fmt=trim(fmt3)) 'r_mar_eff =', r_mar_eff
#if (defined(S_DIS))
write(10, fmt=trim(fmt3)) 's_dis     =', s_dis
#endif
#if (defined(ALPHA_SUB))
write(10, fmt=trim(fmt3)) 'alpha_sub =', alpha_sub
#endif
#if (defined(ALPHA_O))
write(10, fmt=trim(fmt3)) 'alpha_o   =', alpha_o
#endif
#endif
write(10, fmt=trim(fmt1)) ' '

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)) 'smw_coeff   =', smw_coeff
write(10, fmt=trim(fmt2a)) 'r_smw       = ', r_smw
write(10, fmt=trim(fmt2a)) 's_smw       = ', s_smw
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(fmt2a)) 'ICE_STREAM = ', ice_stream
#if (ICE_STREAM==2)
write(10, fmt=trim(fmt1)) 'Sediment-mask file = '//mask_sedi_file
write(10, fmt=trim(fmt3)) 'c_slide_sedi     =', c_slide_sedi
write(10, fmt=trim(fmt3)) 'smw_coeff_sedi   =', smw_coeff_sedi
write(10, fmt=trim(fmt2a)) 'r_smw_sedi       = ', r_smw_sedi
write(10, fmt=trim(fmt2a)) 's_smw_sedi       = ', s_smw_sedi
write(10, fmt=trim(fmt3)) 'gamma_slide_sedi =', gamma_slide_sedi
write(10, fmt=trim(fmt2a)) 'p_weert_sedi     = ', p_weert_sedi
write(10, fmt=trim(fmt2a)) 'q_weert_sedi     = ', q_weert_sedi
#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)) ' '

#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(fmt1)) ' '
write(10, fmt=trim(fmt2)) 'TEMP_PRESENT_PARA =', temp_present_para
write(10, fmt=trim(fmt2)) 'TSURFACE   =', tsurface
write(10, fmt=trim(fmt2)) 'ACCSURFACE =', accsurface
#if (ACCSURFACE==5)
write(10, fmt=trim(fmt2)) 'PRECIP_ANOM_INTERPOL =', precip_anom_interpol
#endif
write(10, fmt=trim(fmt2)) 'SOLID_PRECIP =', solid_precip
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
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

#if (THK_EVOL==2)
time_target_topo_init  = time_target_topo_init0 *year_sec   ! a --> s
time_target_topo_final = time_target_topo_final0*year_sec   ! a --> s
#endif

#if (THK_EVOL==3)
target_topo_tau  = target_topo_tau0 *year_sec   ! a --> s
#endif

time = time_init

!-------- Reading of present-day
!         monthly mean or mean annual precipitation rate --------

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

#if (defined(PRECIP_PRESENT_FILE) && defined(PRECIP_MA_PRESENT_FILE))

if (flag_precip_monthly_mean) then
   filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
                        trim(precip_present_file)
else
   filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
                        trim(precip_ma_present_file)
end if

#endif

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

#elif (GRID==2)

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

#endif

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

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

if (flag_precip_monthly_mean) then

   do n=1, 12   ! month counter
      do m=1, 3; read(21, fmt='(a)') ch_dummy; end do
      do j=jmax, 0, -1
         read(21, fmt=*) (precip_present(j,i,n), i=0,imax)
      end do
   end do

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

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

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

close(21, status='keep')

!  ------ Spatially dependent modification factor

precip_factor = 1.0_dp

#if (defined(PRECIP_FACTOR_FILE))

if ( trim(adjustl(precip_factor_file)) /= 'none' ) then

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

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

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

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

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

   close(21, status='keep')

end if

#endif

if (flag_precip_monthly_mean) then

   do n=1, 12   ! month counter
      do i=0, imax
      do j=0, jmax
         precip_present(j,i,n) = precip_present(j,i,n) * precip_factor(j,i)
      end do
      end do
   end do

else

   precip_ma_present = precip_ma_present * precip_factor

end if

!  ------ Computation of the still undefined present-day
!         mean annual or monthly mean precipitation rate

if (flag_precip_monthly_mean) then

   precip_ma_present = 0.0_dp   ! initialisation

   do n=1, 12   ! month counter
      do i=0, imax
      do j=0, jmax
         precip_ma_present(j,i) = precip_ma_present(j,i) &
                                  + precip_present(j,i,n)*inv_twelve
      end do
      end do
   end do

else

   do n=1, 12   ! month counter
      do i=0, imax
      do j=0, jmax
         precip_present(j,i,n) = precip_ma_present(j,i)
                ! monthly means assumed to be equal
                ! to the mean annual precipitation rate
      end do
      end do
   end do

end if

!  ------ Reference topography for present-day precipitation rate

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

#if (defined(PRECIP_ZS_REF_FILE))

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

#endif

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

#elif (GRID==2)

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

#endif

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

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

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

close(21, status='keep')

do i=0, imax
do j=0, jmax
   zs_ref(j,i) = max(zs_ref(j,i), 0.0_dp)
                 ! resetting negative elevations (bathymetry data) 
                 ! to the present-day sea surface
end do
end do

!-------- Reading of LGM mean-annual precipitation-rate anomaly --------

#if (ACCSURFACE==5)

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

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

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

#endif

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

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

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

close(21, status='keep')

precip_ma_lgm_anom = precip_ma_lgm_anom * precip_anom_fact

!-------- LGM monthly precipitation-rate anomalies (assumed to be
!         equal to the mean annual precipitation-rate anomaly) --------

do n=1, 12   ! month counter
do i=0, imax
do j=0, jmax

   precip_lgm_anom(j,i,n) = max(precip_ma_lgm_anom(j,i), eps)
                            ! positive values ensured

#if (PRECIP_ANOM_INTERPOL==1)
   gamma_precip_lgm_anom(j,i,n) = 0.0_dp   ! dummy values
#elif (PRECIP_ANOM_INTERPOL==2)
   gamma_precip_lgm_anom(j,i,n) = -log(precip_lgm_anom(j,i,n))
#else
   stop ' >>> sico_init: Wrong value of switch PRECIP_ANOM_INTERPOL!'
#endif

end do
end do
end do

#endif

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

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

!-------- Read rock/sediment mask --------

#if (ICE_STREAM==2)

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

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

if (ios /= 0) stop ' >>> sico_init: Error when opening the rock/sediment 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_sedi(j,i), i=0,imax)
end do

close(24, status='keep')

#endif

!-------- Reading of the prescribed target topography --------

#if ((THK_EVOL==2) || (THK_EVOL==3))

target_topo_dat_name = trim(target_topo_dat_name)

#if (NETCDF==1)
write(6, fmt='(a)') ' >>> sico_init: Reading of target topography'
write(6, fmt='(a)') '                only implemented for NetCDF files (NETCDF==2)!'
stop
#elif (NETCDF==2)
call read_target_topo_nc(target_topo_dat_name)
#else
stop ' >>> sico_init: Parameter NETCDF must be either 1 or 2!'
#endif

#endif

!-------- Reading of the maximum ice extent mask --------

#if (THK_EVOL==4)

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

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

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

#elif (GRID==2)

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

#endif

if (ios /= 0) &
   stop ' >>> sico_init: Error when opening the maximum ice extent 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_maxextent(j,i), i=0,imax)
end do

close(24, status='keep')

#endif

!-------- Reading of LGM mean-annual and mean-July surface-temperature 
!         anomalies --------

#if (TSURFACE==5)

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

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

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

#endif

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

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

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

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

#endif

if (ios /= 0) stop ' >>> sico_init: Error when opening the temp_mj anomaly 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 j=jmax, 0, -1
   read(21, fmt=*) (temp_ma_lgm_anom(j,i), i=0,imax)
   read(22, fmt=*) (temp_mj_lgm_anom(j,i), i=0,imax)
end do

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

temp_ma_lgm_anom = temp_ma_lgm_anom * temp_ma_anom_fact
temp_mj_lgm_anom = temp_mj_lgm_anom * temp_mj_anom_fact

#endif

!-------- 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 the glacial index --------

#if (TSURFACE==5)

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

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

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

read(21, fmt=*) ch_dummy, gi_time_min, gi_time_stp, gi_time_max

if (ch_dummy /= '#') then
   write(6, fmt=*) ' >>> sico_init: gi_time_min, gi_time_stp, gi_time_max'
   write(6, fmt=*) '                not defined in glacial-index file!'
   stop
end if

ndata_gi = (gi_time_max-gi_time_min)/gi_time_stp

allocate(glacial_index(0:ndata_gi))

do n=0, ndata_gi
   read(21, fmt=*) d_dummy, glacial_index(n)
end do

close(21, status='keep')

#endif

!-------- Reading of time-dependent surface-temperature
!                        and precipitation anomaly data --------

#if (TSURFACE==6 && ACCSURFACE==6)

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

call check( nf90_open(trim(filename_with_path), nf90_nowrite, &
                      ncid_temp_precip), thisroutine )

call check( nf90_inquire(ncid_temp_precip, unlimiteddimid = dimid), &
            thisroutine )

call check( nf90_inquire_dimension(ncid_temp_precip, dimid, &
                                   len = ndata_temp_precip), thisroutine )

ndata_temp_precip = ndata_temp_precip-1

allocate(temp_precip_time(0:ndata_temp_precip))

call check( nf90_inq_varid(ncid_temp_precip, 'time', ncv), thisroutine )
call check( nf90_get_var(ncid_temp_precip, ncv, temp_precip_time), thisroutine )
                                                ! in a

temp_precip_time_min = temp_precip_time(0)   ! in a
temp_precip_time_stp = temp_precip_time(1)-temp_precip_time(0)
                       ! in a; constant time step assumed
temp_precip_time_max = temp_precip_time(ndata_temp_precip)   ! in a

#endif

!-------- Reading of ISMIP6 InitMIP SMB anomaly data --------

#if (defined(INITMIP_SMB_ANOM_FILE))

if ( trim(adjustl(initmip_smb_anom_file)) /= 'none' ) then

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

   call check( nf90_open(trim(filename_with_path), nf90_nowrite, ncid), &
               thisroutine )

   call check( nf90_inq_varid(ncid, 'DSMB', ncv) )
   call check( nf90_get_var(ncid, ncv, as_anom_initmip_conv) )

   call check( nf90_close(ncid) )

   do i=0, imax
   do j=0, jmax
      as_anom_initmip(j,i) = as_anom_initmip_conv(i,j) /year_sec
                                   ! m/a ice equiv. -> m/s ice equiv.
   end do
   end do

else
   as_anom_initmip = 0.0_dp
end if

#endif

!-------- Reading of the global annual temperature anomaly
!         (for parameterising the sub-ocean temperature anomaly
!         for the ice discharge parameterisation)

#if (DISC==2)

filename_with_path = trim(inpath)//'/general/dTg_paleo.dat'

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

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

read(21, fmt=*) ch_dummy, glann_time_min, glann_time_stp, glann_time_max

if (ch_dummy /= '#') then
   write(6, fmt=*) ' >>> sico_init: glann_time_min, glann_time_stp, glann_time_max'
   write(6, fmt=*) '                not defined in data file for dT_glann_CLIMBER!'
   stop
end if

ndata_glann = (glann_time_max-glann_time_min)/glann_time_stp

allocate(dt_glann_climber(0:ndata_glann))

do n=0, ndata_glann
   read(21, fmt=*) d_dummy, dt_glann_climber(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)

!  ------ 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(j,i), i=0,imax)
end do

close(21, status='keep')

do i=0, imax
do j=0, jmax
   q_geo(j,i) = q_geo(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)

#if (DISC>0)   /* Ice discharge parameterization */
call disc_param(dtime)
call disc_fields()
#endif

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)

#if (DISC>0)   /* Ice discharge parameterization */
call disc_param(dtime)
call disc_fields()
#endif

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)

#if (DISC>0)   /* Ice discharge parameterization */
call disc_param(dtime)
call disc_fields()
#endif

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

if (forcing_flag == 1) then

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

else if (forcing_flag == 2) then

   write(12,1112)
   write(12,1113)

   1112 format('         t(a)  glac_ind(1)      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)')
   1113 format('----------------------------------------------------', &
               '---------------------------------------')

else if (forcing_flag == 3) then

   write(12,1122)
   write(12,1123)

   1122 format('         t(a)   (Dummy)(1)      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)')
   1123 format('----------------------------------------------------', &
               '---------------------------------------')

end if

!  ------ Time-series file for deep boreholes

n_core = 7   ! GRIP, GISP2, Dye3, Camp Century (CC),
             ! NorthGRIP (NGRIP), NEEM, EastGRIP (EGRIP)

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

ch_core(1)     = 'GRIP'
phi_core(1)    =  72.58722_dp *pi_180   ! Geographical position of GRIP,
lambda_core(1) = -37.64222_dp *pi_180   ! conversion deg -> rad
 
ch_core(2)     = 'GISP2'
phi_core(2)    =  72.58833_dp *pi_180   ! Geographical position of GISP2
lambda_core(2) = -38.45750_dp *pi_180   ! conversion deg -> rad

ch_core(3)     = 'Dye3'
phi_core(3)    =  65.15139_dp *pi_180   ! Geographical position of Dye3,
lambda_core(3) = -43.81722_dp *pi_180   ! conversion deg -> rad

ch_core(4)     = 'Camp Century'
phi_core(4)    =  77.17970_dp *pi_180   ! Geographical position of CC,
lambda_core(4) = -61.10975_dp *pi_180   ! conversion deg -> rad

ch_core(5)     = 'NGRIP'
phi_core(5)    =  75.09694_dp *pi_180   ! Geographical position of NGRIP,
lambda_core(5) = -42.31956_dp *pi_180   ! conversion deg -> rad

ch_core(6)     = 'NEEM'
phi_core(6)    =  77.5_dp     *pi_180   ! Geographical position of NEEM,
lambda_core(6) = -50.9_dp     *pi_180   ! conversion deg -> rad

ch_core(7)     = 'EGRIP'
phi_core(7)    =  75.6299_dp  *pi_180   ! Geographical position of EGRIP,
lambda_core(7) = -35.9867_dp  *pi_180   ! conversion deg -> rad

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

do n=1, n_core
   call stereo_forw_ellipsoid(lambda_core(n), phi_core(n), &
                              a, b, lambda0, phi0, x_core(n), y_core(n))
end do

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

x_core = lambda_core
y_core = phi_core

#endif

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_GR(m)      H_G2(m)      H_D3(m)', &
               '      H_CC(m)      H_NG(m)      H_NE(m)      H_EG(m)',/, &
               '                 v_GR(m/a)    v_G2(m/a)    v_D3(m/a)', &
               '    v_CC(m/a)    v_NG(m/a)    v_NE(m/a)    v_EG(m/a)',/, &
               '                   T_GR(C)      T_G2(C)      T_D3(C)', &
               '      T_CC(C)      T_NG(C)      T_NE(C)      T_EG(C)')
   1107 format('----------------------------------------------------', &
               '----------------------------------------------------')

else if (forcing_flag == 2) then

   write(14,1116)
   write(14,1117)

   1116 format('         t(a)  glac_ind(1)      z_sl(m)',/, &
               '                   H_GR(m)      H_G2(m)      H_D3(m)', &
               '      H_CC(m)      H_NG(m)      H_NE(m)      H_EG(m)',/, &
               '                 v_GR(m/a)    v_G2(m/a)    v_D3(m/a)', &
               '    v_CC(m/a)    v_NG(m/a)    v_NE(m/a)    v_EG(m/a)',/, &
               '                   T_GR(C)      T_G2(C)      T_D3(C)', &
               '      T_CC(C)      T_NG(C)      T_NE(C)      T_EG(C)')
   1117 format('----------------------------------------------------', &
               '----------------------------------------------------')

else if (forcing_flag == 3) then

   write(14,1126)
   write(14,1127)

   1126 format('         t(a)   (Dummy)(1)      z_sl(m)',/, &
               '                   H_GR(m)      H_G2(m)      H_D3(m)', &
               '      H_CC(m)      H_NG(m)      H_NE(m)      H_EG(m)',/, &
               '                 v_GR(m/a)    v_G2(m/a)    v_D3(m/a)', &
               '    v_CC(m/a)    v_NG(m/a)    v_NE(m/a)    v_EG(m/a)',/, &
               '                   T_GR(C)      T_G2(C)      T_D3(C)', &
               '      T_CC(C)      T_NG(C)      T_NE(C)      T_EG(C)')
   1127 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

#if (defined(EXEC_MAKE_C_DIS_0))

#if (DISC>0)

call calc_c_dis_0(dxi, deta)

write(6, fmt='(a)') ' '
write(6, fmt='(a)') ' * * * * * * * * * * * * * * * * * * * * * * * * * * * * '
write(6, fmt='(a)') ' '
write(6, fmt='(a)') ' sico_init: Routine calc_c_dis_0 successfully completed, '
write(6, fmt='(a)') '            c_dis_0 written on file out_runname.dat      '
write(6, fmt='(a)') '            (in directory specified by OUTPATH).         '
write(6, fmt='(a)') '            Execution of SICOPOLIS stopped.              '
write(6, fmt='(a)') ' '
write(6, fmt='(a)') ' * * * * * * * * * * * * * * * * * * * * * * * * * * * * '
write(6, fmt='(a)') ' '

stop

#else
  stop ' >>> sico_init: EXEC_MAKE_C_DIS_0 requires DISC>0!'
#endif

#endif

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
real(dp)     :: h_ice, freeboard_ratio
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')

#if (defined(ZB_PRESENT_FILE))

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

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

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

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

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

close(25, status='keep')

#else

write(6, fmt='(a)') ' >>> topography1: ZB_PRESENT_FILE not defined,'
write(6, fmt='(a)') '                  thus zb = zl assumed.'

zb = zl

#endif

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

freeboard_ratio = (rho_sw-rho)/rho_sw

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

   if (maske(j,i) <= 1) then

      zb(j,i) = zl(j,i)   ! ensure consistency

   else if (maske(j,i) == 2) then

#if (MARGIN==1 || MARGIN==2)
      zs(j,i) = zl(j,i)   ! ensure
      zb(j,i) = zl(j,i)   ! consistency
#elif (MARGIN==3)
      zs(j,i) = 0.0_dp    ! present-day
      zb(j,i) = 0.0_dp    ! sea level
#endif

   else if (maske(j,i) == 3) then

#if (MARGIN==1 || (MARGIN==2 && MARINE_ICE_FORMATION==1))
      maske(j,i) = 2     ! floating ice cut off
      zs(j,i) = zl(j,i)
      zb(j,i) = zl(j,i)
#elif (MARGIN==2 && MARINE_ICE_FORMATION==2)
      maske(j,i) = 0     ! floating ice becomes "underwater ice"
      h_ice   = zs(j,i)-zb(j,i)   ! ice thickness
      zs(j,i) = zl(j,i)+h_ice
      zb(j,i) = zl(j,i)
#elif (MARGIN==3)
      h_ice = zs(j,i)-zb(j,i)   ! ice thickness
      zs(j,i) = freeboard_ratio*h_ice   ! ensure properly
      zb(j,i) = zs(j,i)-h_ice           ! floating ice
#endif

   end if

   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 n=1, 6; read(24, fmt='(a)') ch_dummy; end do

do j=jmax, 0, -1
   read(23, fmt=*)          (zl0(j,i),   i=0,imax)
   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)
      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

   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

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