sico_init.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Subroutine :  s i c o _ i n i t
!
!> @file
!!
!! Initialisations for SICOPOLIS.
!!
!! @section Copyright
!!
!! Copyright 2009-2013 Ralf Greve
!!
!! @section License
!!
!! This file is part of SICOPOLIS.
!!
!! SICOPOLIS is free software: you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! SICOPOLIS is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with SICOPOLIS.  If not, see <http://www.gnu.org/licenses/>.
!<
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!-------------------------------------------------------------------------------
!> Initialisations for SICOPOLIS.
!<------------------------------------------------------------------------------
subroutine sico_init(delta_ts, glac_index, &
               sum_accum, mean_accum, mean_accum_inv, &
               dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, &
               time, time_init, time_end, time_output, &
               dxi, deta, dzeta_c, dzeta_t, dzeta_r, &
               z_sl, dzsl_dtau, z_mar, &
               ii, jj, nn, &
               ndat2d, ndat3d, n_output, &
               runname)

use sico_types
use sico_variables

implicit none
integer(i4b) :: i, j, kc, kt, kr, m, n, ir, jr
integer(i4b) :: ios, ios1, ios2, ios3, ios4
integer(i4b) :: ndat2d, ndat3d
integer(i4b) :: n_output
integer(i4b) :: ierr
integer(i4b) :: ii((imax+1)*(jmax+1)), jj((imax+1)*(jmax+1)), nn(0:jmax,0:imax)
real(dp) :: delta_ts, glac_index
real(dp) :: sum_accum, mean_accum, mean_accum_inv
real(dp) :: dtime0, dtime_temp0, dtime_wss0, dtime_out0, dtime_ser0
real(dp) :: dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser
real(dp) :: time_init0, time_end0, time_output0(100)
real(dp) :: time, time_init, time_end, time_output(100)
real(dp) :: dxi, deta, dzeta_c, dzeta_t, dzeta_r
real(dp) :: z_sl, dzsl_dtau, z_mar
real(dp) :: d_dummy
character (len=100) :: runname, anfdatname
character (len=256) :: shell_command
character :: ch_dummy
logical :: flag_3d_output

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

#if (CALCZS==4 || MARGIN==3)

#include "lisf.h"
call lis_initialize(ierr)

#endif

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

call phys_para

temp_min  = temp_min
s_t       = s_t       *1.0e-09_dp            ! K/km3 -> K/m3
x_hat     = x_hat     *1.0e+03_dp            ! km -> m
y_hat     = y_hat     *1.0e+03_dp            ! km -> m
rad       = rad       *1.0e+03_dp            ! km -> m
b_min     = b_min     /year_sec              ! m/a -> m/s
b_max     = b_max     /year_sec              ! m/a -> m/s

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

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

#if GRID==0

if (abs(dx-50.0_dp) < eps) then

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

else if (abs(dx-25.0_dp) < eps) then

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

else
      stop ' sico_init: DX wrong!'
end if

#elif GRID==1

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

#elif GRID==2

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

#endif

!-------- Compatibility check of 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

#if ((ADV_HOR == 3) && (ADV_VERT != 3))
stop ' sico_init: Options ADV_HOR==3 and ADV_VERT==3 must be used together!'
#endif

#if ((ADV_HOR != 3) && (ADV_VERT == 3))
stop ' sico_init: Options ADV_HOR==3 and ADV_VERT==3 must be used together!'
#endif

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

forcing_flag = 1   ! forcing by delta_ts

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

dtime0      = dtime0
dtime_temp0 = dtime_temp0

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

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

ndat2d = 1
ndat3d = 1

!-------- 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.le.imax).and.(j.le.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.

open(10, iostat=ios, &
     file=outpath//'/'//trim(runname)//'.log', &
     status='new')
if (ios /= 0) stop ' sico_init: Error when opening the log file!'

write(10,1000) 'Computational domain:'
#if defined(ANT)
write(10,1000) 'Antarctica'
#elif defined(GRL)
write(10,1000) 'Greenland'
#elif defined(SCAND)
write(10,1000) 'Scandinavia and Eurasia'
#elif defined(NHEM)
write(10,1000) 'Northern hemisphere'
#elif defined(EMTP2SGE)
write(10,1000) 'EISMINT Phase 2 Simplified Geometry Experiment'
#elif defined(HEINO)
write(10,1000) 'ISMIP HEINO'
#elif defined(NMARS)
write(10,1000) 'North polar cap of Mars'
#elif defined(SMARS)
write(10,1000) 'South polar cap of Mars'
#else
stop ' sico_init: Subroutine sico_init: No valid domain specified!'
#endif
write(10,1000) ' '

write(10,1001) 'imax  =', imax
write(10,1001) 'jmax  =', jmax
write(10,1001) 'kcmax =', kcmax
write(10,1001) 'ktmax =', ktmax
write(10,1001) 'krmax =', krmax
write(10,1000) ' '

write(10,1002) 'a =', deform
write(10,1000) ' '

#if (GRID==0 || GRID==1)
write(10,1002) 'x0      =', x0
write(10,1002) 'y0      =', y0
write(10,1002) 'dx      =', dx
#elif GRID==2
write(10,1002) 'dlambda =', dlambda
write(10,1002) 'dphi    =', dphi
#endif
write(10,1000) ' '

write(10,1002) 'year_zero  =', year_zero
write(10,1002) 'time_init  =', time_init0
write(10,1002) 'time_end   =', time_end0
write(10,1002) 'dtime      =', dtime0
write(10,1002) 'dtime_temp =', dtime_temp0
#if ( OUTPUT==1 || OUTPUT==3 )
write(10,1002) 'dtime_out  =', dtime_out0
#endif
write(10,1002) 'dtime_ser  =', dtime_ser0
write(10,1000) ' '

write(10,1000) 'mask_present file = '//mask_present_file
write(10,1000) ' '

write(10,1000) 'Physical-parameter file = '//phys_para_file
write(10,1000) ' '

#if (CALCZS==3 || CALCZS==4)
write(10,1002) 'ovi_weight =', ovi_weight
#if CALCZS==3
write(10,1002) 'omega_sor  =', omega_sor
#endif
write(10,1000) ' '
#endif

write(10,1002) 'temp_min =', temp_min
write(10,1002) 's_t      =', s_t
write(10,1002) 'x_hat    =', x_hat
write(10,1002) 'y_hat    =', y_hat
write(10,1002) 'rad      =', rad
write(10,1002) 'b_min    =', b_min
write(10,1002) 'b_max    =', b_max
write(10,1000) ' '

#if TSURFACE==1
write(10,1002) 'delta_ts0      =', delta_ts0
#elif TSURFACE==3
write(10,1002) 'sine_amplit    =', sine_amplit
write(10,1002) 'sine_period    =', sine_period
#elif TSURFACE==4
write(10,1000) 'GRIP file      = '//grip_temp_file
write(10,1002) 'grip_temp_fact =', grip_temp_fact
#endif

#if SEA_LEVEL==1
write(10,1002) 'z_sl0          =', z_sl0
#elif SEA_LEVEL==3
write(10,1000) 'sea-level file = '//sea_level_file
#endif
write(10,1000) ' '

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

write(10,1000) 'Sediment-mask file = '//mask_sedi_file
write(10,1000) ' '

write(10,1002) 'c_slide          =', c_slide
write(10,1002) 'gamma_slide      =', gamma_slide
write(10,1001) 'p_weert          =', p_weert
write(10,1001) 'q_weert          =', q_weert
#if SLIDE_LAW==2
write(10,1002) 'red_pres_limit_fact =', red_pres_limit_fact
#endif
write(10,1000) ' '

write(10,1002) 'c_slide_sedi     =', c_slide_sedi
write(10,1002) 'gamma_slide_sedi =', gamma_slide_sedi
write(10,1001) 'p_weert_sedi     =', p_weert_sedi
write(10,1001) 'q_weert_sedi     =', q_weert_sedi
write(10,1000) ' '

#if Q_GEO_MOD==1
write(10,1002) 'q_geo     =', q_geo
write(10,1000) ' '
#endif

#if REBOUND==1
write(10,1002) 'frac_llra =', frac_llra
write(10,1000) ' '
#endif

#if FLOW_LAW==2
write(10,1002) 'gr_size   =', gr_size
write(10,1000) ' '
#endif
#if FIN_VISC==2
write(10,1002) 'sigma_res =', sigma_res
write(10,1000) ' '
#endif

write(10,1001) 'ENHMOD =', enhmod
write(10,1002) 'enh_fact    =', enh_fact
#if ( ENHMOD==2 || ENHMOD==3 )
write(10,1002) 'enh_intg    =', enh_intg
#endif
#if ENHMOD==2
write(10,1002) 'age_trans   =', age_trans_0
#endif
#if ENHMOD==3
write(10,1002) 'date_trans1 =', date_trans1_0
write(10,1002) 'date_trans2 =', date_trans2_0
write(10,1002) 'date_trans3 =', date_trans3_0
#endif
#if MARGIN==3
write(10,1002) 'enh_shelf   =', enh_shelf
#endif
write(10,1000) ' '

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

write(10,1001) 'GRID       =', grid
write(10,1001) 'CALCMOD    =', calcmod
write(10,1001) 'FLOW_LAW   =', flow_law
write(10,1001) 'FIN_VISC   =', fin_visc
write(10,1001) 'MARGIN     =', margin
#if MARGIN==2
write(10,1001) 'MARINE_ICE_FORMATION =', marine_ice_formation
write(10,1001) 'MARINE_ICE_CALVING   =', marine_ice_calving
#elif MARGIN==3
write(10,1001) 'ICE_SHELF_CALVING =', ice_shelf_calving
#endif
write(10,1001) 'ANF_DAT    =', anf_dat
write(10,1001) 'REBOUND    =', rebound
write(10,1001) 'Q_LITHO    =', q_litho
write(10,1001) 'ZS_EVOL    =', zs_evol
write(10,1001) 'CALCZS     =', calczs
write(10,1001) 'ADV_HOR    =', adv_hor
write(10,1001) 'ADV_VERT   =', adv_vert
write(10,1001) 'TOPOGRAD   =', topograd
write(10,1001) 'TSURFACE   =', tsurface
write(10,1001) 'SEA_LEVEL  =', sea_level
write(10,1001) 'SLIDE_LAW  =', slide_law
write(10,1001) 'Q_GEO_MOD  =', q_geo_mod
write(10,1000) ' '

#if defined(OUT_TIMES)
write(10,1001) 'OUT_TIMES  =', out_times
#endif
write(10,1001) 'OUTPUT     =', output
write(10,1001) 'OUTSER     =', outser
#if defined(OUTSER_V_A)
write(10,1001) 'OUTSER_V_A =', outser_v_a
#endif
#if ( OUTPUT==1 || OUTPUT==2 )
write(10,1001) 'ERGDAT     =', ergdat
#endif
write(10,1001) 'NETCDF     =', netcdf
write(10,1000) ' '
#if ( OUTPUT==2 || OUTPUT==3 )
write(10,1001) 'n_output =', n_output
do n=1, n_output
   write(10,1002) 'time_output =', time_output0(n)
end do
write(10,1000) ' '
#endif

write(10,1000) 'Program version and date: '//version//' / '//date

 1000 format(a)
 1001 format(a,i4)
 1002 format(a,es12.4)

close(10, status='keep')

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

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

if (abs(dtime_temp/dtime-nint(dtime_temp/dtime)) > eps) &
   stop ' sico_init: dtime_temp must be a multiple of dtime!'

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

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

mean_accum_inv = 1.0_dp/mean_accum

!-------- Read present topography mask --------

open(24, iostat=ios, &
     file=inpath//'/heino/'//mask_present_file, &
     recl=8192, status='old')

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

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

do j=jmax, 0, -1
   read(24,2300) (maske_help(j,i), i=0,imax)
end do

close(24, status='keep')

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

open(24, iostat=ios, &
     file=inpath//'/heino/'//mask_sedi_file, &
     recl=8192, status='old')

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

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

do j=jmax, 0, -1
   read(24,2300) (maske_sedi(j,i), i=0,imax)
end do

close(24, status='keep')

2300 format(imax(i1),i1)

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

#if TSURFACE==4

open(21, iostat=ios, &
     file=inpath//'/general/'//grip_temp_file, &
     status='old')
if (ios /= 0) stop ' sico_init: Error when opening the data file for delta_ts!'

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

if (ch_dummy /= '#') then
   write(6,*) ' sico_init: grip_time_min, grip_time_stp, grip_time_max'
   write(6,*) '            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,*) d_dummy, griptemp(n)
end do

close(21, status='keep')

#endif

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

#if SEA_LEVEL==3

open(21, iostat=ios, &
     file=inpath//'/general/'//sea_level_file, &
     status='old')
if (ios /= 0) stop ' sico_init: Error when opening the data file for z_sl!'

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

if (ch_dummy /= '#') then
   write(6,*) ' sico_init: specmap_time_min, specmap_time_stp, specmap_time_max'
   write(6,*) '            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,*) d_dummy, specmap_zsl(n)
end do

close(21, status='keep')

#endif

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

#if Q_GEO_MOD==1

!  ------ Constant value

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

#elif Q_GEO_MOD==2

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

#endif

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

!  ------ Initial topography with a thin ice layer everywhere

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

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

   q_bm(j,i)    = 0.0_dp
   q_tld(j,i)   = 0.0_dp
   q_b_tot(j,i) = 0.0_dp

   p_b_w(j,i)   = 0.0_dp
   h_w(j,i)     = 0.0_dp

   do kc=0, kcmax
      temp_c(kc,j,i) = temp_s(j,i)
      age_c(kc,j,i)  = 15000.0_dp*year_sec
   end do

   do kt=0, ktmax
      omega_t(kt,j,i) = 0.0_dp
      age_t(kt,j,i)   = 15000.0_dp*year_sec
   end do

   do kr=0, krmax
      temp_r(kr,j,i) = temp_s(j,i)+q_geo(j,i)*h_r/kappa_r &
                       *(1.0_dp-real(kr,dp)/real(krmax,dp))
!           (linear temperature distribution according to the
!            geothermal heat flux)
   end do

end do
end do

call calc_enhance(time_init)

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

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

   q_bm(j,i)    = 0.0_dp
   q_tld(j,i)   = 0.0_dp
   q_b_tot(j,i) = 0.0_dp

   p_b_w(j,i)   = 0.0_dp
   h_w(j,i)     = 0.0_dp

   do kc=0, kcmax
      temp_c(kc,j,i) = temp_s(j,i)
      age_c(kc,j,i)  = 15000.0_dp*year_sec
   end do

   do kt=0, ktmax
      omega_t(kt,j,i) = 0.0_dp
      age_t(kt,j,i)   = 15000.0_dp*year_sec
   end do

   do kr=0, krmax
      temp_r(kr,j,i) = temp_s(j,i)+q_geo(j,i)*h_r/kappa_r &
                       *(1.0_dp-real(kr,dp)/real(krmax,dp))
!           (linear temperature distribution according to the
!            geothermal heat flux)
   end do

end do
end do

call calc_enhance(time_init)

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

#elif ANF_DAT==3

#if NETCDF==1
call topography3(dxi, deta, z_sl, anfdatname)
#elif NETCDF==2
call topography3_nc(dxi, deta, z_sl, anfdatname)
#else
stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
#endif

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

q_b_tot = q_bm + q_tld

call calc_enhance(time_init)

#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

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

ea = exp(deform)

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

do kt=0, ktmax
   zeta_t(kt) = kt*dzeta_t
end do

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

call calc_temp_melt

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

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

call calc_vz_sia(dxi, deta, dzeta_c, dzeta_t)

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

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

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

open(12, iostat=ios, &
     file=outpath//'/'//trim(runname)//'.ser', &
     status='new')
if (ios /= 0) stop ' sico_init: Error when opening the ser file!'

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

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

   1102 format('         t(a)  D_Ts(deg C)      z_sl(m)',/, &
               '                  H_max(m)    zs_max(m)     V_g(m^3)', &
               '     V_t(m^3)  V_fw(m^3/a)     V_sle(m)',/, &
               '                  A_g(m^2)     A_t(m^2)  V_bm(m^3/a)', &
               ' V_tld(m^3/a)   H_t_max(m)  vs_max(m/a)')
   1103 format('----------------------------------------------------', &
               '---------------------------------------')

#elif OUTSER_V_A==2

   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)',/, &
               '                               H_max(m)    zs_max(m)', &
               '     V_t(m^3)  V_fw(m^3/a)     V_sle(m)',/, &
               '                               A_t(m^2)  V_bm(m^3/a)', &
               ' V_tld(m^3/a)   H_t_max(m)  vs_max(m/a)')
   1103 format('----------------------------------------------------', &
               '---------------------------------------')

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

!  ------ Time-series file for the sediment area ("Hudson Bay, Hudson Strait")

open(15, iostat=ios, &
     file=outpath//'/'//trim(runname)//'.sed', &
     status='new')
if (ios /= 0) stop ' sico_init: Error when opening the sed file!'

write(15,1108)
write(15,1109)

1108 format('         t(a)  D_Ts(deg C)      z_sl(m)',/, &
            '                  H_ave(m)   Tbh_ave(C)     Atb(m^2)')
1109 format('----------------------------------------------------')

!  ------ Time-series file for the ice-thickness field

#if OUTSER==2

open(13, iostat=ios, &
     file=outpath//'/'//trim(runname)//'.thk', &
     status='new')
if (ios /= 0) stop ' sico_init: Error when opening the thk file!'

write(13,1104)
write(13,1105)

1104 format(' t(a)  D_Ts(deg C)  z_sl(m)',/, &
            ' H(m): JMAX+1 records [j = JMAX (-1) 0] with',/, &
            '       IMAX+1 values [i = 0 (1) IMAX] in each record')
1105 format('----------------------------------------------------')

#endif

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

#if OUTSER==3

n_core = 7   ! Points P1 - P7

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

lambda_core(1) =    0.0_dp  ! dummy
phi_core(1)    =    0.0_dp  ! dummy
x_core(1)      = 3900.0_dp *1.0e+03_dp    ! Point P1,
y_core(1)      = 2000.0_dp *1.0e+03_dp    ! conversion km -> m

lambda_core(2) =    0.0_dp  ! dummy
phi_core(2)    =    0.0_dp  ! dummy
x_core(2)      = 3800.0_dp *1.0e+03_dp    ! Point P2,
y_core(2)      = 2000.0_dp *1.0e+03_dp    ! conversion km -> m

lambda_core(3) =    0.0_dp  ! dummy
phi_core(3)    =    0.0_dp  ! dummy
x_core(3)      = 3700.0_dp *1.0e+03_dp    ! Point P3,
y_core(3)      = 2000.0_dp *1.0e+03_dp    ! conversion km -> m

lambda_core(4) =    0.0_dp  ! dummy
phi_core(4)    =    0.0_dp  ! dummy
x_core(4)      = 3500.0_dp *1.0e+03_dp    ! Point P4,
y_core(4)      = 2000.0_dp *1.0e+03_dp    ! conversion km -> m

lambda_core(5) =    0.0_dp  ! dummy
phi_core(5)    =    0.0_dp  ! dummy
x_core(5)      = 3200.0_dp *1.0e+03_dp    ! Point P5,
y_core(5)      = 2000.0_dp *1.0e+03_dp    ! conversion km -> m

lambda_core(6) =    0.0_dp  ! dummy
phi_core(6)    =    0.0_dp  ! dummy
x_core(6)      = 2900.0_dp *1.0e+03_dp    ! Point P6,
y_core(6)      = 2000.0_dp *1.0e+03_dp    ! conversion km -> m

lambda_core(7) =    0.0_dp  ! dummy
phi_core(7)    =    0.0_dp  ! dummy
x_core(7)      = 2600.0_dp *1.0e+03_dp    ! Point P7,
y_core(7)      = 2000.0_dp *1.0e+03_dp    ! conversion km -> m

open(14, iostat=ios, &
     file=outpath//'/'//trim(runname)//'.core', &
     status='new')
if (ios /= 0) stop ' sico_init: Error when opening the core file!'

if (forcing_flag == 1) then

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

   1106 format('         t(a)      D_Ts(C)      z_sl(m)',/, &
               '                   H_P1(m)      H_P2(m)      H_P3(m)      H_P4(m)', &
               '      H_P5(m)      H_P6(m)      H_P7(m)',/, &
               '                 v_P1(m/a)    v_P2(m/a)    v_P3(m/a)    v_P4(m/a)', &
               '    v_P5(m/a)    v_P6(m/a)    v_P7(m/a)',/, &
               '                   T_P1(C)      T_P2(C)      T_P3(C)      T_P4(C)', &
               '      T_P5(C)      T_P6(C)      T_P7(C)',/, &
               '               Rb_P1(W/m2)  Rb_P2(W/m2)  Rb_P3(W/m2)  Rb_P4(W/m2)', &
               '  Rb_P5(W/m2)  Rb_P6(W/m2)  Rb_P7(W/m2)')
   1107 format('-----------------------------------------------------------------', &
               '---------------------------------------')

end if

#endif

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

#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 NETCDF==1
   call output1(runname, time_init, delta_ts, glac_index, z_sl, &
                flag_3d_output, ndat2d, ndat3d)
#elif NETCDF==2
   call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
                  flag_3d_output, ndat2d, ndat3d)
#else
   stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
#endif

#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

#if NETCDF==1
   call output1(runname, time_init, delta_ts, glac_index, z_sl, &
                flag_3d_output, ndat2d, ndat3d)
#elif NETCDF==2
   call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
                  flag_3d_output, ndat2d, ndat3d)
#else
   stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
#endif

end if

#elif OUTPUT==3

   flag_3d_output = .false.

#if NETCDF==1
   call output1(runname, time_init, delta_ts, glac_index, z_sl, &
                flag_3d_output, ndat2d, ndat3d)
#elif NETCDF==2
   call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
                  flag_3d_output, ndat2d, ndat3d)
#else
   stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
#endif

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

   flag_3d_output = .true.

#if NETCDF==1
   call output1(runname, time_init, delta_ts, glac_index, z_sl, &
                flag_3d_output, ndat2d, ndat3d)
#elif NETCDF==2
   call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
                  flag_3d_output, ndat2d, ndat3d)
#else
   stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
#endif

end if

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

call output2(time_init, dxi, deta, delta_ts, glac_index, z_sl)
#if OUTSER==2
call output3(time_init, delta_ts, glac_index, z_sl)
#endif
#if OUTSER==3
call output4(time_init, dxi, deta, delta_ts, glac_index, z_sl)
#endif

end subroutine sico_init
!