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-2016 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, &
               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 sico_types
use sico_variables
use sico_vars
use ice_material_quantities, only : ice_mat_eqs_pars
use enth_temp_omega, only : calc_c_int_table, calc_c_int_inv_table, &
                            enth_fct_temp_omega

implicit none

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

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

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

!-------- 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 (CALCZS==4 || MARGIN==3 || DYNAMICS==2)

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

#endif

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

call phys_para()

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

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

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

#elif (SURFACE_FORCING==2)

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

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

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

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

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

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

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

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

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

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

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

#if GRID==0

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

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

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

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

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

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

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

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

else
      stop ' sico_init: DX wrong!'
end if

#elif GRID==1

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

#elif GRID==2

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

#endif

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

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

if (ktmax /= 1) 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 1.'
   write(6, fmt='(a)') ' '
end if

#endif

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

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

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

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

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

forcing_flag = 1   ! forcing by delta_ts

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

dtime0      = dtime0
dtime_temp0 = dtime_temp0

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

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

ndat2d = 1
ndat3d = 1

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

!  ------ kc domain

if (deform >= eps) then

   flag_aa_nonzero = .true.   ! non-equidistant grid

   aa = deform
   ea = exp(aa)

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

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

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

else

   flag_aa_nonzero = .false.   ! equidistant grid

   aa = 0.0_dp
   ea = 1.0_dp

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

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

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

end if

!  ------ kt domain

kt=0
zeta_t(kt) = 0.0_dp

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

kt=ktmax
zeta_t(kt) = 1.0_dp

!  ------ kr domain

kr=0
zeta_r(kr) = 0.0_dp

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

kr=krmax
zeta_r(kr) = 1.0_dp

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

n=1

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

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

runname    = runname
anfdatname = anfdatname

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

time_init0 = time_init0
time_end0  = time_end0
dtime_ser0 = dtime_ser0

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

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

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

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

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, fmt=trim(fmt1)) 'Computational domain:'
write(10, fmt=trim(fmt1)) trim(ch_domain_long)
write(10, fmt=trim(fmt1)) ' '

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

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

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

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

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

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

#if (CALCZS==3 || CALCZS==4)
write(10, fmt=trim(fmt3)) 'ovi_weight =', ovi_weight
#if CALCZS==3
write(10, fmt=trim(fmt3)) 'omega_sor  =', omega_sor
#endif
write(10, fmt=trim(fmt1)) ' '
#endif

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

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

#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

write(10, fmt=trim(fmt3)) 'c_slide     =', c_slide
write(10, fmt=trim(fmt3)) 'gamma_slide =', gamma_slide
write(10, fmt=trim(fmt2)) 'p_weert     =', p_weert
write(10, fmt=trim(fmt2)) 'q_weert     =', q_weert
#if SLIDE_LAW==2
write(10, fmt=trim(fmt3)) 'red_pres_limit_fact =', red_pres_limit_fact
#endif
write(10, fmt=trim(fmt1)) ' '

#if Q_GEO_MOD==1
write(10, fmt=trim(fmt3)) 'q_geo     =', q_geo
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)) 'H_min       =', h_min
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(fmt2)) 'GRID       =', grid
write(10, fmt=trim(fmt2)) 'DYNAMICS   =', dynamics
#if (DYNAMICS==2 && defined(RATIO_SL_THRESH))
write(10, fmt=trim(fmt3)) 'ratio_sl_thresh =', ratio_sl_thresh
#endif
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
write(10, fmt=trim(fmt2)) 'ZS_EVOL    =', zs_evol
write(10, fmt=trim(fmt2)) 'CALCZS     =', calczs
write(10, fmt=trim(fmt2)) 'ADV_HOR    =', adv_hor
write(10, fmt=trim(fmt2)) 'ADV_VERT   =', adv_vert
write(10, fmt=trim(fmt2)) 'TOPOGRAD   =', topograd
write(10, fmt=trim(fmt2)) 'TSURFACE   =', tsurface
write(10, fmt=trim(fmt2)) 'SEA_LEVEL  =', sea_level
write(10, fmt=trim(fmt2)) 'SLIDE_LAW  =', slide_law
write(10, fmt=trim(fmt2)) 'Q_GEO_MOD  =', q_geo_mod
write(10, fmt=trim(fmt1)) ' '

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

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

close(10, status='keep')

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

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

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

time = time_init

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

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

!-------- Set present topography mask --------

maske_help = 1

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

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

!-------- 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, fmt=*) ch_dummy, grip_time_min, grip_time_stp, grip_time_max

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

ndata_grip = (grip_time_max-grip_time_min)/grip_time_stp

allocate(griptemp(0:ndata_grip))

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

close(21, status='keep')

#endif

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

#if (SEA_LEVEL==3)

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, fmt=*) ch_dummy, specmap_time_min, specmap_time_stp, specmap_time_max

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

ndata_specmap = (specmap_time_max-specmap_time_min)/specmap_time_stp

allocate(specmap_zsl(0:ndata_specmap))

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

close(21, status='keep')

#endif

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

#if Q_GEO_MOD==1

!  ------ Constant value

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

#elif Q_GEO_MOD==2

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

#endif

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

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

#if (TIME_LAG_MOD==1)

time_lag_asth = time_lag*year_sec   ! a -> s

#elif (TIME_LAG_MOD==2)

open(29, iostat=ios, &
     file=inpath//'/'//trim(ch_domain_short)//'/'//time_lag_file, &
     recl=8192, 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)

open(29, iostat=ios, &
     file=inpath//'/'//trim(ch_domain_short)//'/'//flex_rig_file, &
     recl=8192, status='old')

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

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

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

close(29, status='keep')

#endif

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

flex_rig_lith = 0.0_dp   ! dummy values

#endif

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

!  ------ Present topography

#if (ANF_DAT==1)

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

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

#elif (ANF_DAT==2)

call topography2(dxi, deta)

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

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

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

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

#endif

!-------- 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_sia(dxi, deta, dzeta_c, dzeta_t)

#if (MARGIN==3)
call calc_vz_ssa(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

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 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 = 2   ! central dome, position halfway to coast

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)      =  750.0_dp *1.0e+03_dp    ! Position of the central dome
y_core(1)      =  750.0_dp *1.0e+03_dp    ! (750 km, 750 km),
                                          ! conversion km -> m

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

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)',/, &
               '                 v_P1(m/a)    v_P2(m/a)',/, &
               '                   T_P1(C)      T_P2(C)',/, &
               '               Rb_P1(W/m2)  Rb_P2(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
!