SICOPOLIS V3.2
 All Classes Files Functions Variables Macros
sico_init.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Subroutine : s i c o _ i n i t
4 !
5 !> @file
6 !!
7 !! Initialisations for SICOPOLIS.
8 !!
9 !! @section Copyright
10 !!
11 !! Copyright 2009-2016 Ralf Greve
12 !!
13 !! @section License
14 !!
15 !! This file is part of SICOPOLIS.
16 !!
17 !! SICOPOLIS is free software: you can redistribute it and/or modify
18 !! it under the terms of the GNU General Public License as published by
19 !! the Free Software Foundation, either version 3 of the License, or
20 !! (at your option) any later version.
21 !!
22 !! SICOPOLIS is distributed in the hope that it will be useful,
23 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
24 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25 !! GNU General Public License for more details.
26 !!
27 !! You should have received a copy of the GNU General Public License
28 !! along with SICOPOLIS. If not, see <http://www.gnu.org/licenses/>.
29 !<
30 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
31 
32 !-------------------------------------------------------------------------------
33 !> Initialisations for SICOPOLIS.
34 !<------------------------------------------------------------------------------
35 subroutine sico_init(delta_ts, glac_index, &
36  mean_accum, &
37  dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, &
38  time, time_init, time_end, time_output, &
39  dxi, deta, dzeta_c, dzeta_t, dzeta_r, &
40  z_sl, dzsl_dtau, z_mar, &
41  ii, jj, nn, &
42  ndat2d, ndat3d, n_output, &
43  runname)
44 
45 use sico_types
46 use sico_variables
47 use sico_vars
48 use ice_material_quantities, only : ice_mat_eqs_pars
49 use enth_temp_omega, only : calc_c_int_table, calc_c_int_inv_table, &
50  enth_fct_temp_omega
51 
52 #if (NETCDF > 1)
53 use netcdf
54 use nc_check
55 #endif
56 
57 implicit none
58 
59 integer(i4b), intent(out) :: ii((imax+1)*(jmax+1)), &
60  jj((imax+1)*(jmax+1)), &
61  nn(0:jmax,0:imax)
62 integer(i4b), intent(out) :: ndat2d, ndat3d
63 integer(i4b), intent(out) :: n_output
64 real(dp), intent(out) :: delta_ts, glac_index
65 real(dp), intent(out) :: mean_accum
66 real(dp), intent(out) :: dtime, dtime_temp, dtime_wss, &
67  dtime_out, dtime_ser
68 real(dp), intent(out) :: time, time_init, time_end, time_output(100)
69 real(dp), intent(out) :: dxi, deta, dzeta_c, dzeta_t, dzeta_r
70 real(dp), intent(out) :: z_sl, dzsl_dtau, z_mar
71 character(len=100), intent(out) :: runname
72 
73 integer(i4b) :: i, j, kc, kt, kr, m, n, ir, jr
74 integer(i4b) :: ios, ios1, ios2, ios3, ios4
75 integer(i4b) :: ierr
76 real(dp) :: dtime0, dtime_temp0, dtime_wss0, dtime_out0, dtime_ser0
77 real(dp) :: time_init0, time_end0, time_output0(100)
78 real(dp), dimension(0:JMAX,0:IMAX) :: precip_factor
79 real(dp) :: d_dummy
80 real(dp) :: transition_length_sliding, quasi_time, d_quasi_time
81 real(dp), dimension(0:JMAX,0:IMAX) :: r_mask_sedi_new
82 character(len=100) :: anfdatname, target_topo_dat_name
83 character(len=256) :: filename_with_path
84 character(len=256) :: shell_command
85 character :: ch_dummy
86 logical :: flag_3d_output
87 
88 #if (NETCDF > 1)
89 integer(i4b) :: dimid
90 ! dimid: Dimension ID
91 integer(i4b) :: ncv
92 ! ncv: Variable ID
93 #endif
94 
95 character(len=64), parameter :: thisroutine = 'sico_init'
96 
97 character(len=64), parameter :: fmt1 = '(a)', &
98  fmt2 = '(a,i4)', &
99  fmt3 = '(a,es12.4)'
100 
101 character(len= 8) :: ch_imax
102 character(len=128) :: fmt4
103 
104 write(ch_imax, fmt='(i8)') imax
105 write(fmt4, fmt='(a)') '('//trim(adjustl(ch_imax))//'(i1),i1)'
106 
107 !-------- Name of the computational domain --------
108 
109 #if defined(ANT)
110 ch_domain_long = 'Antarctica'
111 ch_domain_short = 'ant'
112 
113 #elif defined(ASF)
114 ch_domain_long = 'Austfonna'
115 ch_domain_short = 'asf'
116 
117 #elif defined(EMTP2SGE)
118 ch_domain_long = 'EISMINT Phase 2 Simplified Geometry Experiment'
119 ch_domain_short = 'emtp2sge'
120 
121 #elif defined(GRL)
122 ch_domain_long = 'Greenland'
123 ch_domain_short = 'grl'
124 
125 #elif defined(NHEM)
126 ch_domain_long = 'Northern hemisphere'
127 ch_domain_short = 'nhem'
128 
129 #elif defined(SCAND)
130 ch_domain_long = 'Scandinavia and Eurasia'
131 ch_domain_short = 'scand'
132 
133 #elif defined(TIBET)
134 ch_domain_long = 'Tibet'
135 ch_domain_short = 'tibet'
136 
137 #elif defined(NMARS)
138 ch_domain_long = 'North polar cap of Mars'
139 ch_domain_short = 'nmars'
140 
141 #elif defined(SMARS)
142 ch_domain_long = 'South polar cap of Mars'
143 ch_domain_short = 'smars'
144 
145 #elif (defined(XYZ))
146 ch_domain_long = 'XYZ'
147 ch_domain_short = 'xyz'
148 #if (defined(HEINO))
149 ch_domain_long = trim(ch_domain_long)//'/ISMIP HEINO'
150 #endif
151 
152 #else
153 stop ' sico_init: No valid domain specified!'
154 #endif
155 
156 !-------- Some initial values --------
157 
158 n_output = 0
159 
160 dtime = 0.0_dp
161 dtime_temp = 0.0_dp
162 dtime_wss = 0.0_dp
163 dtime_out = 0.0_dp
164 dtime_ser = 0.0_dp
165 
166 time = 0.0_dp
167 time_init = 0.0_dp
168 time_end = 0.0_dp
169 time_output = 0.0_dp
170 
171 !-------- Initialisation of the Library of Iterative Solvers Lis,
172 ! if required --------
173 
174 #if (CALCZS==4 || MARGIN==3 || DYNAMICS==2)
175 
176 #include "lisf.h"
177 call lis_initialize(ierr)
178 
179 #endif
180 
181 !-------- Read physical parameters --------
182 
183 call phys_para()
184 
185 call ice_mat_eqs_pars(rf, r_t, kappa, c, -190, 10)
186 
187 ! ------ Some auxiliary quantities required for the enthalpy method
188 
189 call calc_c_int_table(c, -190, 10, l)
190 call calc_c_int_inv_table()
191 
192 !-------- Compatibility check of the SICOPOLIS version with the header file
193 
194 if ( trim(version) /= trim(sico_version) ) &
195  stop ' sico_init: SICOPOLIS version not compatible with header file!'
196 
197 !-------- Check whether the dynamics and thermodynamics modes are defined
198 
199 #if ( !defined(DYNAMICS) )
200 stop ' sico_init: DYNAMICS not defined in the header file!'
201 #endif
202 
203 #if ( !defined(CALCMOD) )
204 stop ' sico_init: CALCMOD not defined in the header file!'
205 #endif
206 
207 #if ( defined(ENTHMOD) )
208 write(6, fmt='(a)') ' sico_init: ENTHMOD must not be defined any more.'
209 write(6, fmt='(a)') ' Please update your header file!'
210 stop
211 #endif
212 
213 !-------- Compatibility check of the horizontal resolution with the
214 ! number of grid points --------
215 
216 #if (GRID==0 || GRID==1)
217 
218 if (abs(dx-40.0_dp) < eps) then
219 
220  if ( (imax /= 150).or.(jmax /= 150) ) then
221  stop ' sico_init: IMAX and/or JMAX wrong!'
222  end if
223 
224 else if (abs(dx-20.0_dp) < eps) then
225 
226  if ( (imax /= 300).or.(jmax /= 300) ) then
227  stop ' sico_init: IMAX and/or JMAX wrong!'
228  end if
229 
230 else if (abs(dx-10.0_dp) < eps) then
231 
232  if ( (imax /= 600).or.(jmax /= 600) ) then
233  stop ' sico_init: IMAX and/or JMAX wrong!'
234  end if
235 
236 else
237  stop ' sico_init: DX wrong!'
238 end if
239 
240 #elif GRID==2
241 
242 stop ' sico_init: GRID==2 not allowed for this application!'
243 
244 #endif
245 
246 !-------- Compatibility check of the thermodynamics mode
247 ! (cold vs. polythermal vs. enthalpy method)
248 ! and the number of grid points in the lower (kt) ice domain --------
249 
250 #if (CALCMOD==0 || CALCMOD==2 || CALCMOD==3 || CALCMOD==-1)
251 
252 if (ktmax /= 1) then
253  write(6, fmt='(a)') ' sico_init: For options CALCMOD==0, 2, 3 or -1,'
254  write(6, fmt='(a)') ' the separate kt domain is redundant.'
255  write(6, fmt='(a)') ' Therefore, consider setting KTMAX to 1.'
256  write(6, fmt='(a)') ' '
257 end if
258 
259 #endif
260 
261 !-------- Compatibility check of surface-temperature
262 ! and precipitation switches --------
263 
264 #if ((TSURFACE == 5) && (ACCSURFACE != 5))
265 stop ' sico_init: Options TSURFACE==5 and ACCSURFACE==5 must be used together!'
266 #endif
267 
268 #if ((TSURFACE != 5) && (ACCSURFACE == 5))
269 stop ' sico_init: Options TSURFACE==5 and ACCSURFACE==5 must be used together!'
270 #endif
271 
272 #if (TSURFACE == 6)
273 #if (ACCSURFACE != 6)
274 write(6, fmt='(a)') ' sico_init: Options TSURFACE==6, ACCSURFACE==6'
275 write(6, fmt='(a)') ' and NETCDF==2 must be used together!'
276 stop
277 #elif (NETCDF != 2)
278 write(6, fmt='(a)') ' sico_init: Options TSURFACE==6, ACCSURFACE==6'
279 write(6, fmt='(a)') ' and NETCDF==2 must be used together!'
280 stop
281 #endif
282 #endif
283 
284 #if (ACCSURFACE == 6)
285 #if (TSURFACE != 6)
286 write(6, fmt='(a)') ' sico_init: Options TSURFACE==6, ACCSURFACE==6'
287 write(6, fmt='(a)') ' and NETCDF==2 must be used together!'
288 stop
289 #elif (NETCDF != 2)
290 write(6, fmt='(a)') ' sico_init: Options TSURFACE==6, ACCSURFACE==6'
291 write(6, fmt='(a)') ' and NETCDF==2 must be used together!'
292 stop
293 #endif
294 #endif
295 
296 !-------- Compatibility check of discretization schemes for the horizontal and
297 ! vertical advection terms in the temperature and age equations --------
298 
299 #if (ADV_HOR==1)
300 stop ' sico_init: Option ADV_HOR==1 (central differences) not defined!'
301 #endif
302 
303 !-------- Check whether for the shallow shelf
304 ! or shelfy stream approximation
305 ! the chosen grid is Cartesian coordinates
306 ! without distortion correction (GRID==0) --------
307 
308 #if ((MARGIN==3 && DYNAMICS==1) || DYNAMICS==2) /* requires SSA or SStA */
309 #if (GRID != 0)
310 write(6, fmt='(a)') ' sico_init: Distortion correction not implemented'
311 write(6, fmt='(a)') ' for the shallow shelf approximation (SSA)'
312 write(6, fmt='(a)') ' or the shelfy stream approximation (SStA)'
313 write(6, fmt='(a)') ' -> GRID==0 required!'
314 stop
315 #endif
316 #endif
317 
318 !-------- Setting of forcing flag --------
319 
320 #if (TSURFACE <= 4)
321 
322 forcing_flag = 1 ! forcing by delta_ts
323 
324 #elif (TSURFACE == 5)
325 
326 forcing_flag = 2 ! forcing by glac_index
327 
328 #elif (TSURFACE == 6)
329 
330 forcing_flag = 3 ! forcing by time-dependent surface temperature
331  ! and precipitation data
332 
333 #endif
334 
335 !-------- Initialization of numerical time steps --------
336 
337 dtime0 = dtime0
338 dtime_temp0 = dtime_temp0
339 #if (REBOUND==2)
340 dtime_wss0 = dtime_wss0
341 #endif
342 
343 !-------- Further initializations --------
344 
345 dzeta_c = 1.0_dp/real(kcmax,dp)
346 dzeta_t = 1.0_dp/real(ktmax,dp)
347 dzeta_r = 1.0_dp/real(krmax,dp)
348 
349 ndat2d = 1
350 ndat3d = 1
351 
352 !-------- General abbreviations --------
353 
354 ! ------ kc domain
355 
356 if (deform >= eps) then
357 
358  flag_aa_nonzero = .true. ! non-equidistant grid
359 
360  aa = deform
361  ea = exp(aa)
362 
363  kc=0
364  zeta_c(kc) = 0.0_dp
365  eaz_c(kc) = 1.0_dp
366  eaz_c_quotient(kc) = 0.0_dp
367 
368  do kc=1, kcmax-1
369  zeta_c(kc) = kc*dzeta_c
370  eaz_c(kc) = exp(aa*zeta_c(kc))
371  eaz_c_quotient(kc) = (eaz_c(kc)-1.0_dp)/(ea-1.0_dp)
372  end do
373 
374  kc=kcmax
375  zeta_c(kc) = 1.0_dp
376  eaz_c(kc) = exp(aa)
377  eaz_c_quotient(kc) = 1.0_dp
378 
379 else
380 
381  flag_aa_nonzero = .false. ! equidistant grid
382 
383  aa = 0.0_dp
384  ea = 1.0_dp
385 
386  kc=0
387  zeta_c(kc) = 0.0_dp
388  eaz_c(kc) = 1.0_dp
389  eaz_c_quotient(kc) = 0.0_dp
390 
391  do kc=1, kcmax-1
392  zeta_c(kc) = kc*dzeta_c
393  eaz_c(kc) = 1.0_dp
394  eaz_c_quotient(kc) = zeta_c(kc)
395  end do
396 
397  kc=kcmax
398  zeta_c(kc) = 1.0_dp
399  eaz_c(kc) = 1.0_dp
400  eaz_c_quotient(kc) = 1.0_dp
401 
402 end if
403 
404 ! ------ kt domain
405 
406 kt=0
407 zeta_t(kt) = 0.0_dp
408 
409 do kt=1, ktmax-1
410  zeta_t(kt) = kt*dzeta_t
411 end do
412 
413 kt=ktmax
414 zeta_t(kt) = 1.0_dp
415 
416 ! ------ kr domain
417 
418 kr=0
419 zeta_r(kr) = 0.0_dp
420 
421 do kr=1, krmax-1
422  zeta_r(kr) = kr*dzeta_r
423 end do
424 
425 kr=krmax
426 zeta_r(kr) = 1.0_dp
427 
428 !-------- Construction of a vector (with index n) from a 2-d array
429 ! (with indices i, j) by diagonal numbering --------
430 
431 n=1
432 
433 do m=0, imax+jmax
434  do i=m, 0, -1
435  j = m-i
436  if ((i <= imax).and.(j <= jmax)) then
437  ii(n) = i
438  jj(n) = j
439  nn(j,i) = n
440  n=n+1
441  end if
442  end do
443 end do
444 
445 !-------- Specification of current simulation --------
446 
447 runname = runname
448 anfdatname = anfdatname
449 
450 #if defined(YEAR_ZERO)
451 year_zero = year_zero
452 #else
453 year_zero = 2000.0_dp ! default value 2000 CE
454 #endif
455 
456 time_init0 = time_init0
457 time_end0 = time_end0
458 dtime_ser0 = dtime_ser0
459 
460 #if ( OUTPUT==1 || OUTPUT==3 )
461 dtime_out0 = dtime_out0
462 #endif
463 
464 #if ( OUTPUT==2 || OUTPUT==3 )
465 n_output = n_output
466 time_output0( 1) = time_out0_01
467 time_output0( 2) = time_out0_02
468 time_output0( 3) = time_out0_03
469 time_output0( 4) = time_out0_04
470 time_output0( 5) = time_out0_05
471 time_output0( 6) = time_out0_06
472 time_output0( 7) = time_out0_07
473 time_output0( 8) = time_out0_08
474 time_output0( 9) = time_out0_09
475 time_output0(10) = time_out0_10
476 time_output0(11) = time_out0_11
477 time_output0(12) = time_out0_12
478 time_output0(13) = time_out0_13
479 time_output0(14) = time_out0_14
480 time_output0(15) = time_out0_15
481 time_output0(16) = time_out0_16
482 time_output0(17) = time_out0_17
483 time_output0(18) = time_out0_18
484 time_output0(19) = time_out0_19
485 time_output0(20) = time_out0_20
486 #endif
487 
488 !-------- Write log file --------
489 
490 shell_command = 'if [ ! -d'
491 shell_command = trim(shell_command)//' '//outpath
492 shell_command = trim(shell_command)//' '//'] ; then mkdir'
493 shell_command = trim(shell_command)//' '//outpath
494 shell_command = trim(shell_command)//' '//'; fi'
495 call system(trim(shell_command))
496  ! Check whether directory OUTPATH exists. If not, it is created.
497 
498 open(10, iostat=ios, &
499  file=outpath//'/'//trim(runname)//'.log', &
500  status='new')
501 if (ios /= 0) &
502  stop ' sico_init: Error when opening the log file!'
503 
504 write(10, fmt=trim(fmt1)) 'Computational domain:'
505 write(10, fmt=trim(fmt1)) trim(ch_domain_long)
506 write(10, fmt=trim(fmt1)) ' '
507 
508 write(10, fmt=trim(fmt2)) 'imax =', imax
509 write(10, fmt=trim(fmt2)) 'jmax =', jmax
510 write(10, fmt=trim(fmt2)) 'kcmax =', kcmax
511 write(10, fmt=trim(fmt2)) 'ktmax =', ktmax
512 write(10, fmt=trim(fmt2)) 'krmax =', krmax
513 write(10, fmt=trim(fmt1)) ' '
514 
515 write(10, fmt=trim(fmt3)) 'a =', aa
516 write(10, fmt=trim(fmt1)) ' '
517 
518 #if (GRID==0 || GRID==1)
519 write(10, fmt=trim(fmt3)) 'x0 =', x0
520 write(10, fmt=trim(fmt3)) 'y0 =', y0
521 write(10, fmt=trim(fmt3)) 'dx =', dx
522 #elif GRID==2
523 stop ' sico_init: GRID==2 not allowed for this application!'
524 #endif
525 write(10, fmt=trim(fmt1)) ' '
526 
527 write(10, fmt=trim(fmt3)) 'year_zero =', year_zero
528 write(10, fmt=trim(fmt3)) 'time_init =', time_init0
529 write(10, fmt=trim(fmt3)) 'time_end =', time_end0
530 write(10, fmt=trim(fmt3)) 'dtime =', dtime0
531 write(10, fmt=trim(fmt3)) 'dtime_temp =', dtime_temp0
532 #if (REBOUND==2)
533 write(10, fmt=trim(fmt3)) 'dtime_wss =', dtime_wss0
534 #endif
535 #if ( OUTPUT==1 || OUTPUT==3 )
536 write(10, fmt=trim(fmt3)) 'dtime_out =', dtime_out0
537 #endif
538 write(10, fmt=trim(fmt3)) 'dtime_ser =', dtime_ser0
539 write(10, fmt=trim(fmt1)) ' '
540 
541 write(10, fmt=trim(fmt2)) 'ANF_DAT =', anf_dat
542 write(10, fmt=trim(fmt1)) 'zs_present file = '//zs_present_file
543 #if (ANF_DAT==1)
544 write(10, fmt=trim(fmt1)) 'zl_present file = '//zl_present_file
545 #endif
546 write(10, fmt=trim(fmt1)) 'zl0 file = '//zl0_file
547 write(10, fmt=trim(fmt1)) 'mask_present file = '//mask_present_file
548 #if (ANF_DAT==1 && defined(TEMP_INIT))
549 write(10, fmt=trim(fmt2)) 'TEMP_INIT =', temp_init
550 #endif
551 #if (ANF_DAT==3)
552 write(10, fmt=trim(fmt1)) 'Initial-value file = '//anfdatname
553 write(10, fmt=trim(fmt1)) 'Path to initial-value file = '//anfdatpath
554 #endif
555 write(10, fmt=trim(fmt1)) ' '
556 
557 #if ZS_EVOL==2
558 write(10, fmt=trim(fmt3)) 'time_target_topo_init =', time_target_topo_init0
559 write(10, fmt=trim(fmt3)) 'time_target_topo_final =', time_target_topo_final0
560 write(10, fmt=trim(fmt1)) 'Target-topography file = '//target_topo_dat_name
561 write(10, fmt=trim(fmt1)) 'Path to target-topography file = '//target_topo_dat_path
562 write(10, fmt=trim(fmt1)) ' '
563 #endif
564 
565 #if ZS_EVOL==3
566 write(10, fmt=trim(fmt1)) 'Maximum ice extent mask file = '//mask_maxextent_file
567 write(10, fmt=trim(fmt1)) ' '
568 #endif
569 
570 write(10, fmt=trim(fmt1)) 'Physical-parameter file = '//phys_para_file
571 write(10, fmt=trim(fmt1)) ' '
572 
573 #if (CALCZS==3 || CALCZS==4)
574 write(10, fmt=trim(fmt3)) 'ovi_weight =', ovi_weight
575 #if (CALCZS==3)
576 write(10, fmt=trim(fmt3)) 'omega_sor =', omega_sor
577 #endif
578 write(10, fmt=trim(fmt1)) ' '
579 #endif
580 
581 #if (TSURFACE==1)
582 write(10, fmt=trim(fmt3)) 'delta_ts0 =', delta_ts0
583 #elif (TSURFACE==3)
584 write(10, fmt=trim(fmt3)) 'sine_amplit =', sine_amplit
585 write(10, fmt=trim(fmt3)) 'sine_period =', sine_period
586 #elif (TSURFACE==4)
587 write(10, fmt=trim(fmt1)) 'GRIP file = '//grip_temp_file
588 write(10, fmt=trim(fmt3)) 'grip_temp_fact =', grip_temp_fact
589 #elif (TSURFACE==5)
590 write(10, fmt=trim(fmt1)) 'Glacial-index file = '//glac_ind_file
591 write(10, fmt=trim(fmt1)) 'temp_ma_anom file = '//temp_ma_anom_file
592 write(10, fmt=trim(fmt3)) 'temp_ma_anom fact = ', temp_ma_anom_fact
593 write(10, fmt=trim(fmt1)) 'temp_mj_anom file = '//temp_mj_anom_file
594 write(10, fmt=trim(fmt3)) 'temp_mj_anom fact = ', temp_mj_anom_fact
595 #elif (TSURFACE==6)
596 write(10, fmt=trim(fmt1)) 'temp_precip_anom file = '//temp_precip_anom_file
597 write(10, fmt=trim(fmt3)) 'temp_ma_anom fact = ', temp_ma_anom_fact
598 write(10, fmt=trim(fmt3)) 'temp_mj_anom fact = ', temp_mj_anom_fact
599 #endif
600 
601 write(10, fmt=trim(fmt1)) 'precip_present_file = '//precip_present_file
602 #if (defined(PRECIP_FACTOR_FILE))
603 if ( trim(adjustl(precip_factor_file)) /= 'none' ) &
604  write(10, fmt=trim(fmt1)) 'precip_factor_file = '//precip_factor_file
605 #endif
606 #if (ACCSURFACE==1)
607 write(10, fmt=trim(fmt3)) 'accfact =', accfact
608 #elif (ACCSURFACE==2 || ACCSURFACE==3)
609 write(10, fmt=trim(fmt3)) 'gamma_s =', gamma_s
610 #elif (ACCSURFACE==5)
611 write(10, fmt=trim(fmt1)) 'precip_anom file = '//precip_anom_file
612 write(10, fmt=trim(fmt3)) 'precip_anom fact = ', precip_anom_fact
613 #elif (ACCSURFACE==6)
614 write(10, fmt=trim(fmt1)) 'temp_precip_anom file = '//temp_precip_anom_file
615 write(10, fmt=trim(fmt3)) 'precip_anom fact = ', precip_anom_fact
616 #endif
617 #if (ACCSURFACE <= 3)
618 write(10, fmt=trim(fmt2)) 'ELEV_DESERT =', elev_desert
619 #if (ELEV_DESERT == 1)
620 write(10, fmt=trim(fmt3)) 'gamma_p =', gamma_p
621 write(10, fmt=trim(fmt3)) 'zs_thresh =', zs_thresh
622 #endif
623 #endif
624 
625 #if (ABLSURFACE==3)
626 write(10, fmt=trim(fmt3)) 'lambda_lti =', lambda_lti
627 write(10, fmt=trim(fmt3)) 'temp_lti =', temp_lti
628 #endif
629 
630 #if (SEA_LEVEL==1)
631 write(10, fmt=trim(fmt3)) 'z_sl0 =', z_sl0
632 #elif (SEA_LEVEL==3)
633 write(10, fmt=trim(fmt1)) 'sea-level file = '//sea_level_file
634 #endif
635 write(10, fmt=trim(fmt1)) ' '
636 
637 #if (MARGIN==2)
638 #if ( MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3 )
639 write(10, fmt=trim(fmt3)) 'z_mar =', z_mar
640 write(10, fmt=trim(fmt1)) ' '
641 #elif ( MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5 \
642  || marine_ice_calving==6 || marine_ice_calving==7 )
643 write(10, fmt=trim(fmt3)) 'fact_z_mar =', fact_z_mar
644 write(10, fmt=trim(fmt1)) ' '
645 #elif ( ( MARINE_ICE_FORMATION==2 ) && ( MARINE_ICE_CALVING==9 ) )
646 write(10, fmt=trim(fmt3)) 'calv_uw_coeff =', calv_uw_coeff
647 write(10, fmt=trim(fmt3)) 'r1_calv_uw =', r1_calv_uw
648 write(10, fmt=trim(fmt3)) 'r2_calv_uw =', r2_calv_uw
649 write(10, fmt=trim(fmt1)) ' '
650 #endif
651 #elif (MARGIN==3)
652 #if ICE_SHELF_CALVING==2
653 write(10, fmt=trim(fmt3)) 'H_calv =', h_calv
654 write(10, fmt=trim(fmt1)) ' '
655 #endif
656 #endif
657 
658 write(10, fmt=trim(fmt3)) 'c_slide =', c_slide
659 write(10, fmt=trim(fmt3)) 'gamma_slide =', gamma_slide
660 write(10, fmt=trim(fmt2)) 'p_weert =', p_weert
661 write(10, fmt=trim(fmt2)) 'q_weert =', q_weert
662 #if SLIDE_LAW==2
663 write(10, fmt=trim(fmt3)) 'red_pres_limit_fact =', red_pres_limit_fact
664 #endif
665 write(10, fmt=trim(fmt1)) ' '
666 
667 #if ( defined(SEDI_SLIDE) && SEDI_SLIDE==2 )
668 write(10, fmt=trim(fmt1)) 'Sediment-mask file = '//mask_sedi_file
669 write(10, fmt=trim(fmt1)) ' '
670 
671 write(10, fmt=trim(fmt3)) 'c_slide_sedi =', c_slide_sedi
672 write(10, fmt=trim(fmt3)) 'gamma_slide_sedi =', gamma_slide_sedi
673 write(10, fmt=trim(fmt2)) 'p_weert_sedi =', p_weert_sedi
674 write(10, fmt=trim(fmt2)) 'q_weert_sedi =', q_weert_sedi
675 #if ( defined(TRANS_LENGTH_SL) )
676 write(10, fmt=trim(fmt3)) 'trans_length_sl =', trans_length_sl
677 #endif
678 write(10, fmt=trim(fmt1)) ' '
679 #endif
680 
681 #if Q_GEO_MOD==1
682 write(10, fmt=trim(fmt3)) 'q_geo =', q_geo
683 #elif Q_GEO_MOD==2
684 write(10, fmt=trim(fmt1)) 'q_geo file = '//q_geo_file
685 #endif
686 write(10, fmt=trim(fmt1)) ' '
687 
688 #if (defined(MARINE_ICE_BASAL_MELTING))
689 write(10, fmt=trim(fmt2)) 'MARINE_ICE_BASAL_MELTING =', marine_ice_basal_melting
690 #if (MARINE_ICE_BASAL_MELTING==2 || MARINE_ICE_BASAL_MELTING==3)
691 write(10, fmt=trim(fmt3)) 'qbm_marine =', qbm_marine
692 #endif
693 write(10, fmt=trim(fmt1)) ' '
694 #endif
695 
696 #if (MARGIN==3)
697 write(10, fmt=trim(fmt2)) 'FLOATING_ICE_BASAL_MELTING =', floating_ice_basal_melting
698 #if (FLOATING_ICE_BASAL_MELTING<=3)
699 write(10, fmt=trim(fmt3)) 'qbm_float_1 =', qbm_float_1
700 write(10, fmt=trim(fmt3)) 'qbm_float_2 =', qbm_float_2
701 #endif
702 write(10, fmt=trim(fmt3)) 'qbm_float_3 =', qbm_float_3
703 write(10, fmt=trim(fmt3)) 'z_abyss =', z_abyss
704 #if (FLOATING_ICE_BASAL_MELTING==4)
705 #if (defined(H_W_0))
706 write(10, fmt=trim(fmt3)) 'H_w_0 =', h_w_0
707 #endif
708 #endif
709 write(10, fmt=trim(fmt1)) ' '
710 #endif
711 
712 write(10, fmt=trim(fmt2)) 'REBOUND =', rebound
713 #if (REBOUND==1)
714 write(10, fmt=trim(fmt3)) 'frac_llra =', frac_llra
715 #endif
716 #if (REBOUND==1 || REBOUND==2)
717 write(10, fmt=trim(fmt2)) 'TIME_LAG_MOD =', time_lag_mod
718 #if (TIME_LAG_MOD==1)
719 write(10, fmt=trim(fmt3)) 'time_lag =', time_lag
720 #elif (TIME_LAG_MOD==2)
721 write(10, fmt=trim(fmt1)) 'time_lag_file = '//time_lag_file
722 #else
723 stop ' sico_init: TIME_LAG_MOD must be either 1 or 2!'
724 #endif
725 #endif
726 #if (REBOUND==2)
727 write(10, fmt=trim(fmt2)) 'FLEX_RIG_MOD =', flex_rig_mod
728 #if (FLEX_RIG_MOD==1)
729 write(10, fmt=trim(fmt3)) 'flex_rig =', flex_rig
730 #elif (FLEX_RIG_MOD==2)
731 write(10, fmt=trim(fmt1)) 'flex_rig_file = '//flex_rig_file
732 #else
733 stop ' sico_init: FLEX_RIG_MOD must be either 1 or 2!'
734 #endif
735 #endif
736 write(10, fmt=trim(fmt2)) 'Q_LITHO =', q_litho
737 write(10, fmt=trim(fmt1)) ' '
738 
739 #if FLOW_LAW==2
740 write(10, fmt=trim(fmt3)) 'gr_size =', gr_size
741 write(10, fmt=trim(fmt1)) ' '
742 #endif
743 #if FIN_VISC==2
744 write(10, fmt=trim(fmt3)) 'sigma_res =', sigma_res
745 write(10, fmt=trim(fmt1)) ' '
746 #endif
747 
748 write(10, fmt=trim(fmt2)) 'ENHMOD =', enhmod
749 #if (ENHMOD==1 || ENHMOD==2 || ENHMOD==3)
750 write(10, fmt=trim(fmt3)) 'enh_fact =', enh_fact
751 #endif
752 #if (ENHMOD==2 || ENHMOD==3)
753 write(10, fmt=trim(fmt3)) 'enh_intg =', enh_intg
754 #endif
755 #if (ENHMOD==2)
756 write(10, fmt=trim(fmt3)) 'age_trans =', age_trans_0
757 #endif
758 #if (ENHMOD==3)
759 write(10, fmt=trim(fmt3)) 'date_trans1 =', date_trans1_0
760 write(10, fmt=trim(fmt3)) 'date_trans2 =', date_trans2_0
761 write(10, fmt=trim(fmt3)) 'date_trans3 =', date_trans3_0
762 #endif
763 #if (ENHMOD==4 || ENHMOD==5)
764 write(10, fmt=trim(fmt3)) 'enh_compr =', enh_compr
765 write(10, fmt=trim(fmt3)) 'enh_shear =', enh_shear
766 #endif
767 #if ((ENHMOD==1 || ENHMOD==2 || ENHMOD==3 || ENHMOD==4) && MARGIN==3)
768 write(10, fmt=trim(fmt3)) 'enh_shelf =', enh_shelf
769 #endif
770 write(10, fmt=trim(fmt1)) ' '
771 
772 write(10, fmt=trim(fmt3)) 'numdiff_H_t =', numdiff_h_t
773 write(10, fmt=trim(fmt3)) 'tau_cts =', tau_cts
774 write(10, fmt=trim(fmt3)) 'H_min =', h_min
775 write(10, fmt=trim(fmt3)) 'vh_max =', vh_max
776 write(10, fmt=trim(fmt3)) 'hd_min =', hd_min
777 write(10, fmt=trim(fmt3)) 'hd_max =', hd_max
778 #if defined(QBM_MIN)
779 write(10, fmt=trim(fmt3)) 'qbm_min =', qbm_min
780 #elif defined(QB_MIN) /* obsolete */
781 write(10, fmt=trim(fmt3)) 'qb_min =', qb_min
782 #endif
783 #if defined(QBM_MAX)
784 write(10, fmt=trim(fmt3)) 'qbm_max =', qbm_max
785 #elif defined(QB_MAX) /* obsolete */
786 write(10, fmt=trim(fmt3)) 'qb_max =', qb_max
787 #endif
788 write(10, fmt=trim(fmt3)) 'age_min =', age_min
789 write(10, fmt=trim(fmt3)) 'age_max =', age_max
790 write(10, fmt=trim(fmt3)) 'mean_accum =', mean_accum
791 #if ADV_VERT==1
792 write(10, fmt=trim(fmt3)) 'age_diff =', agediff
793 #endif
794 write(10, fmt=trim(fmt1)) ' '
795 
796 write(10, fmt=trim(fmt2)) 'GRID =', grid
797 write(10, fmt=trim(fmt2)) 'DYNAMICS =', dynamics
798 #if (DYNAMICS==2 && defined(RATIO_SL_THRESH))
799 write(10, fmt=trim(fmt3)) 'ratio_sl_thresh =', ratio_sl_thresh
800 #endif
801 write(10, fmt=trim(fmt2)) 'CALCMOD =', calcmod
802 #if ( CALCMOD==-1 && defined(TEMP_CONST) )
803 write(10, fmt=trim(fmt3)) 'TEMP_CONST =', temp_const
804 #endif
805 #if ( CALCMOD==-1 && defined(AGE_CONST) )
806 write(10, fmt=trim(fmt3)) 'AGE_CONST =', age_const
807 #endif
808 #if ( CALCMOD==1 && defined(CTS_MELTING_FREEZING) )
809 write(10, fmt=trim(fmt2)) 'CTS_MELTING_FREEZING =', cts_melting_freezing
810 #endif
811 write(10, fmt=trim(fmt2)) 'FLOW_LAW =', flow_law
812 write(10, fmt=trim(fmt2)) 'FIN_VISC =', fin_visc
813 write(10, fmt=trim(fmt2)) 'MARGIN =', margin
814 #if (MARGIN==2)
815 write(10, fmt=trim(fmt2)) 'MARINE_ICE_FORMATION =', marine_ice_formation
816 write(10, fmt=trim(fmt2)) 'MARINE_ICE_CALVING =', marine_ice_calving
817 #elif (MARGIN==3)
818 write(10, fmt=trim(fmt2)) 'ICE_SHELF_CALVING =', ice_shelf_calving
819 #endif
820 write(10, fmt=trim(fmt2)) 'ZS_EVOL =', zs_evol
821 write(10, fmt=trim(fmt2)) 'CALCZS =', calczs
822 write(10, fmt=trim(fmt2)) 'ADV_HOR =', adv_hor
823 write(10, fmt=trim(fmt2)) 'ADV_VERT =', adv_vert
824 write(10, fmt=trim(fmt2)) 'TOPOGRAD =', topograd
825 write(10, fmt=trim(fmt2)) 'TEMP_PRESENT_PARA =', temp_present_para
826 write(10, fmt=trim(fmt2)) 'TSURFACE =', tsurface
827 write(10, fmt=trim(fmt2)) 'ACCSURFACE =', accsurface
828 #if ACCSURFACE==5
829 write(10, fmt=trim(fmt2)) 'PRECIP_ANOM_INTERPOL =', precip_anom_interpol
830 #endif
831 write(10, fmt=trim(fmt2)) 'SOLID_PRECIP =', solid_precip
832 write(10, fmt=trim(fmt2)) 'ABLSURFACE =', ablsurface
833 write(10, fmt=trim(fmt2)) 'SEA_LEVEL =', sea_level
834 write(10, fmt=trim(fmt2)) 'SLIDE_LAW =', slide_law
835 #if defined(SEDI_SLIDE)
836 write(10, fmt=trim(fmt2)) 'SEDI_SLIDE =', sedi_slide
837 #endif
838 write(10, fmt=trim(fmt2)) 'Q_GEO_MOD =', q_geo_mod
839 write(10, fmt=trim(fmt1)) ' '
840 
841 #if defined(OUT_TIMES)
842 write(10, fmt=trim(fmt2)) 'OUT_TIMES =', out_times
843 #endif
844 write(10, fmt=trim(fmt2)) 'OUTPUT =', output
845 write(10, fmt=trim(fmt2)) 'OUTSER =', outser
846 #if defined(OUTSER_V_A)
847 write(10, fmt=trim(fmt2)) 'OUTSER_V_A =', outser_v_a
848 #endif
849 #if ( OUTPUT==1 || OUTPUT==2 )
850 write(10, fmt=trim(fmt2)) 'ERGDAT =', ergdat
851 #endif
852 write(10, fmt=trim(fmt2)) 'NETCDF =', netcdf
853 write(10, fmt=trim(fmt1)) ' '
854 #if ( OUTPUT==2 || OUTPUT==3 )
855 write(10, fmt=trim(fmt2)) 'n_output =', n_output
856 do n=1, n_output
857  write(10, fmt=trim(fmt3)) 'time_output =', time_output0(n)
858 end do
859 write(10, fmt=trim(fmt1)) ' '
860 #endif
861 
862 write(10, fmt=trim(fmt1)) 'Program version and date: '//version//' / '//date
863 
864 close(10, status='keep')
865 
866 !-------- Conversion of time quantities --------
867 
868 year_zero = year_zero*year_sec ! a --> s
869 time_init = time_init0*year_sec ! a --> s
870 time_end = time_end0*year_sec ! a --> s
871 dtime = dtime0*year_sec ! a --> s
872 dtime_temp = dtime_temp0*year_sec ! a --> s
873 #if (REBOUND==2)
874 dtime_wss = dtime_wss0*year_sec ! a --> s
875 #endif
876 dtime_ser = dtime_ser0*year_sec ! a --> s
877 #if ( OUTPUT==1 || OUTPUT==3 )
878 dtime_out = dtime_out0*year_sec ! a --> s
879 #endif
880 #if ( OUTPUT==2 || OUTPUT==3 )
881 do n=1, n_output
882  time_output(n) = time_output0(n)*year_sec ! a --> s
883 end do
884 #endif
885 
886 if (abs(dtime_temp/dtime-nint(dtime_temp/dtime)) > eps) &
887  stop ' sico_init: dtime_temp must be a multiple of dtime!'
888 
889 #if (REBOUND==2)
890 if (abs(dtime_wss/dtime-nint(dtime_wss/dtime)) > eps) &
891  stop ' sico_init: dtime_wss must be a multiple of dtime!'
892 #endif
893 
894 #if ZS_EVOL==2
895 time_target_topo_init = time_target_topo_init0 *year_sec ! a --> s
896 time_target_topo_final = time_target_topo_final0*year_sec ! a --> s
897 #endif
898 
899 time = time_init
900 
901 !-------- Reading of present mean-annual precipitation rate --------
902 
903 #if (GRID==0 || GRID==1)
904 
905 open(21, iostat=ios, &
906  file=inpath//'/'//trim(ch_domain_short)//'/'//precip_present_file, &
907  recl=8192, status='old')
908 
909 #elif (GRID==2)
910 
911 stop ' sico_init: GRID==2 not allowed for this application!'
912 
913 #endif
914 
915 if (ios /= 0) stop ' sico_init: Error when opening the precip file!'
916 
917 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
918 
919 do j=jmax, 0, -1
920  read(21, fmt=*) (precip_ma_present(j,i), i=0,imax)
921 end do
922 
923 close(21, status='keep')
924 
925 precip_ma_present = precip_ma_present *(1.0e-03_dp/year_sec)*(rho_w/rho)
926  ! mm/a water equiv. -> m/s ice equiv.
927 
928 precip_factor = 1.0_dp
929 
930 #if (defined(PRECIP_FACTOR_FILE))
931 
932 if ( trim(adjustl(precip_factor_file)) /= 'none' ) then
933 
934  open(21, iostat=ios, &
935  file=inpath//'/'//trim(ch_domain_short)//'/'//precip_factor_file, &
936  recl=8192, status='old')
937 
938  if (ios /= 0) stop ' sico_init: Error when opening the precip factor file!'
939 
940  do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
941 
942  do j=jmax, 0, -1
943  read(21, fmt=*) (precip_factor(j,i), i=0,imax)
944  end do
945 
946  close(21, status='keep')
947 
948 end if
949 
950 #endif
951 
952 precip_ma_present = precip_ma_present * precip_factor
953 
954 !-------- Present monthly precipitation rates
955 ! (assumed to be equal to the mean annual precipitation rate) --------
956 
957 do n=1, 12 ! month counter
958 do i=0, imax
959 do j=0, jmax
960  precip_present(j,i,n) = precip_ma_present(j,i)
961 end do
962 end do
963 end do
964 
965 !-------- Reading of LGM mean-annual precipitation-rate anomaly --------
966 
967 #if ACCSURFACE==5
968 
969 #if (GRID==0 || GRID==1)
970 
971 open(21, iostat=ios, &
972  file=inpath//'/'//trim(ch_domain_short)//'/'//precip_anom_file, &
973  recl=8192, status='old')
974 
975 #endif
976 
977 if (ios /= 0) stop ' sico_init: Error when opening the precip anomaly file!'
978 
979 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
980 
981 do j=jmax, 0, -1
982  read(21, fmt=*) (precip_ma_lgm_anom(j,i), i=0,imax)
983 end do
984 
985 close(21, status='keep')
986 
987 precip_ma_lgm_anom = precip_ma_lgm_anom * precip_anom_fact
988 
989 !-------- LGM monthly precipitation-rate anomalies (assumed to be
990 ! equal to the mean annual precipitation-rate anomaly) --------
991 
992 do n=1, 12 ! month counter
993 do i=0, imax
994 do j=0, jmax
995 
996  precip_lgm_anom(j,i,n) = max(precip_ma_lgm_anom(j,i), eps)
997  ! positive values ensured
998 
999 #if (PRECIP_ANOM_INTERPOL==1)
1000  gamma_precip_lgm_anom(j,i,n) = 0.0_dp ! dummy values
1001 #elif (PRECIP_ANOM_INTERPOL==2)
1002  gamma_precip_lgm_anom(j,i,n) = -log(precip_lgm_anom(j,i,n))
1003 #else
1004  stop ' sico_init: Wrong value of switch PRECIP_ANOM_INTERPOL!'
1005 #endif
1006 
1007 end do
1008 end do
1009 end do
1010 
1011 #endif
1012 
1013 !-------- Mean accumulation --------
1014 
1015 mean_accum = mean_accum*(1.0e-03_dp/year_sec)*(rho_w/rho)
1016  ! mm/a water equiv. -> m/s ice equiv.
1017 
1018 !-------- Read present topography mask --------
1019 
1020 #if (GRID==0 || GRID==1)
1021 
1022 open(24, iostat=ios, &
1023  file=inpath//'/'//trim(ch_domain_short)//'/'//mask_present_file, &
1024  recl=2048, status='old')
1025 
1026 #elif GRID==2
1027 
1028 stop ' sico_init: GRID==2 not allowed for this application!'
1029 
1030 #endif
1031 
1032 if (ios /= 0) stop ' sico_init: Error when opening the mask file!'
1033 
1034 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
1035 
1036 do j=jmax, 0, -1
1037  read(24, fmt=trim(fmt4)) (maske_help(j,i), i=0,imax)
1038 end do
1039 
1040 close(24, status='keep')
1041 
1042 !-------- Read rock/sediment mask --------
1043 
1044 #if ( defined(SEDI_SLIDE) && SEDI_SLIDE==2 )
1045 
1046 open(24, iostat=ios, &
1047  file=inpath//'/'//trim(ch_domain_short)//'/'//mask_sedi_file, &
1048  recl=2048, status='old')
1049 
1050 if (ios /= 0) stop ' sico_init: Error when opening the rock/sediment mask file!'
1051 
1052 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
1053 
1054 do j=jmax, 0, -1
1055  read(24, fmt=trim(fmt4)) (maske_sedi(j,i), i=0,imax)
1056 end do
1057 
1058 close(24, status='keep')
1059 
1060 do i=0, imax
1061 do j=0, jmax
1062  if (maske_sedi(j,i) == 2) maske_sedi(j,i) = 7
1063  ! ocean to be treated like soft sediment
1064 end do
1065 end do
1066 
1067 #endif
1068 
1069 !-------- Initialisation of the marker for the different sectors for
1070 ! ice shelf basal melting --------
1071 
1072 n_sector = 0_i2b
1073 
1074 !-------- Present reference elevation (for precipitation data) --------
1075 
1076 #if (GRID==0 || GRID==1)
1077 
1078 open(21, iostat=ios, &
1079  file=inpath//'/'//trim(ch_domain_short)//'/'//zs_present_file, &
1080  recl=8192, status='old')
1081 
1082 #elif GRID==2
1083 
1084 stop ' sico_init: GRID==2 not allowed for this application!'
1085 
1086 #endif
1087 
1088 if (ios /= 0) stop ' sico_init: Error when opening the zs file!'
1089 
1090 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
1091 
1092 do j=jmax, 0, -1
1093  read(21, fmt=*) (zs_ref(j,i), i=0,imax)
1094 end do
1095 
1096 close(21, status='keep')
1097 
1098 ! ------ Reset bathymetry data (sea floor elevation) to the
1099 ! sea surface
1100 
1101 do i=0, imax
1102 do j=0, jmax
1103  if (maske_help(j,i) >= 2) zs_ref(j,i) = 0.0_dp
1104 end do
1105 end do
1106 
1107 !-------- Reading of the prescribed target topography --------
1108 
1109 #if ZS_EVOL==2
1110 
1111 target_topo_dat_name = trim(target_topo_dat_name)
1112 
1113 #if NETCDF==1
1114 write(6, fmt='(a)') ' sico_init: Reading of target topography'
1115 write(6, fmt='(a)') ' only implemented for NetCDF files (NETCDF==2)!'
1116 stop
1117 #elif NETCDF==2
1118 call read_target_topo_nc(target_topo_dat_name)
1119 #else
1120 stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
1121 #endif
1122 
1123 #endif
1124 
1125 !-------- Reading of the maximum ice extent mask --------
1126 
1127 #if ZS_EVOL==3
1128 
1129 #if (GRID==0 || GRID==1)
1130 
1131 open(24, iostat=ios, &
1132  file=inpath//'/'//trim(ch_domain_short)//'/'//mask_maxextent_file, &
1133  recl=2048, status='old')
1134 
1135 #elif GRID==2
1136 
1137 stop ' sico_init: GRID==2 not allowed for this application!'
1138 
1139 #endif
1140 
1141 if (ios /= 0) &
1142  stop ' sico_init: Error when opening the maximum ice extent mask file!'
1143 
1144 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
1145 
1146 do j=jmax, 0, -1
1147  read(24, fmt=trim(fmt4)) (maske_maxextent(j,i), i=0,imax)
1148 end do
1149 
1150 close(24, status='keep')
1151 
1152 #endif
1153 
1154 !-------- Reading of LGM mean-annual and mean-January (summer)
1155 ! surface-temperature anomalies --------
1156 
1157 #if TSURFACE==5
1158 
1159 #if (GRID==0 || GRID==1)
1160 
1161 open(21, iostat=ios, &
1162  file=inpath//'/'//trim(ch_domain_short)//'/'//temp_ma_anom_file, &
1163  recl=8192, status='old')
1164 
1165 #endif
1166 
1167 if (ios /= 0) stop ' sico_init: Error when opening the temp_ma anomaly file!'
1168 
1169 #if (GRID==0 || GRID==1)
1170 
1171 open(22, iostat=ios, &
1172  file=inpath//'/'//trim(ch_domain_short)//'/'//temp_mj_anom_file, &
1173  recl=8192, status='old')
1174 
1175 #endif
1176 
1177 if (ios /= 0) stop ' sico_init: Error when opening the temp_mj anomaly file!'
1178 
1179 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
1180 do n=1, 6; read(22, fmt='(a)') ch_dummy; end do
1181 
1182 do j=jmax, 0, -1
1183  read(21, fmt=*) (temp_ma_lgm_anom(j,i), i=0,imax)
1184  read(22, fmt=*) (temp_mj_lgm_anom(j,i), i=0,imax)
1185 end do
1186 
1187 close(21, status='keep')
1188 close(22, status='keep')
1189 
1190 temp_ma_lgm_anom = temp_ma_lgm_anom * temp_ma_anom_fact
1191 temp_mj_lgm_anom = temp_mj_lgm_anom * temp_mj_anom_fact
1192 
1193 #endif
1194 
1195 !-------- Read data for delta_ts --------
1196 
1197 #if TSURFACE==4
1198 
1199 open(21, iostat=ios, &
1200  file=inpath//'/general/'//grip_temp_file, &
1201  status='old')
1202 if (ios /= 0) &
1203  stop ' sico_init: Error when opening the data file for delta_ts!'
1204 
1205 read(21, fmt=*) ch_dummy, grip_time_min, grip_time_stp, grip_time_max
1206 
1207 if (ch_dummy /= '#') then
1208  write(6, fmt=*) ' sico_init: grip_time_min, grip_time_stp, grip_time_max'
1209  write(6, fmt=*) ' not defined in data file for delta_ts!'
1210  stop
1211 end if
1212 
1213 ndata_grip = (grip_time_max-grip_time_min)/grip_time_stp
1214 
1215 allocate(griptemp(0:ndata_grip))
1216 
1217 do n=0, ndata_grip
1218  read(21, fmt=*) d_dummy, griptemp(n)
1219 end do
1220 
1221 close(21, status='keep')
1222 
1223 #endif
1224 
1225 !-------- Read data for the glacial index --------
1226 
1227 #if TSURFACE==5
1228 
1229 open(21, iostat=ios, &
1230  file=inpath//'/general/'//glac_ind_file, status='old')
1231 if (ios /= 0) stop ' sico_init: Error when opening the glacial-index file!'
1232 
1233 read(21, fmt=*) ch_dummy, gi_time_min, gi_time_stp, gi_time_max
1234 
1235 if (ch_dummy /= '#') then
1236  write(6, fmt=*) ' sico_init: gi_time_min, gi_time_stp, gi_time_max'
1237  write(6, fmt=*) ' not defined in glacial-index file!'
1238  stop
1239 end if
1240 
1241 ndata_gi = (gi_time_max-gi_time_min)/gi_time_stp
1242 
1243 allocate(glacial_index(0:ndata_gi))
1244 
1245 do n=0, ndata_gi
1246  read(21, fmt=*) d_dummy, glacial_index(n)
1247 end do
1248 
1249 close(21, status='keep')
1250 
1251 #endif
1252 
1253 !-------- Reading of time-dependent surface-temperature
1254 ! and precipitation anomaly data --------
1255 
1256 #if ( (TSURFACE==6) && (ACCSURFACE==6) )
1257 
1258 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)// &
1259  '/'//trim(temp_precip_anom_file)
1260 
1261 call check( nf90_open(trim(filename_with_path), nf90_nowrite, &
1262  ncid_temp_precip), thisroutine )
1263 
1264 call check( nf90_inquire(ncid_temp_precip, unlimiteddimid = dimid), &
1265  thisroutine )
1266 
1267 call check( nf90_inquire_dimension(ncid_temp_precip, dimid, &
1268  len = ndata_temp_precip), thisroutine )
1269 
1270 ndata_temp_precip = ndata_temp_precip-1
1271 
1272 allocate(temp_precip_time(0:ndata_temp_precip))
1273 
1274 call check( nf90_inq_varid(ncid_temp_precip, 'time', ncv), thisroutine )
1275 call check( nf90_get_var(ncid_temp_precip, ncv, temp_precip_time), thisroutine )
1276  ! in a
1277 
1278 temp_precip_time_min = temp_precip_time(0) ! in a
1279 temp_precip_time_stp = temp_precip_time(1)-temp_precip_time(0)
1280  ! in a; constant time step assumed
1281 temp_precip_time_max = temp_precip_time(ndata_temp_precip) ! in a
1282 
1283 #endif
1284 
1285 !-------- Read data for z_sl --------
1286 
1287 #if (SEA_LEVEL==3)
1288 
1289 open(21, iostat=ios, &
1290  file=inpath//'/general/'//sea_level_file, &
1291  status='old')
1292 if (ios /= 0) &
1293  stop ' sico_init: Error when opening the data file for z_sl!'
1294 
1295 read(21, fmt=*) ch_dummy, specmap_time_min, specmap_time_stp, specmap_time_max
1296 
1297 if (ch_dummy /= '#') then
1298  write(6, fmt=*) ' sico_init: specmap_time_min, specmap_time_stp, specmap_time_max'
1299  write(6, fmt=*) ' not defined in data file for z_sl!'
1300  stop
1301 end if
1302 
1303 ndata_specmap = (specmap_time_max-specmap_time_min)/specmap_time_stp
1304 
1305 allocate(specmap_zsl(0:ndata_specmap))
1306 
1307 do n=0, ndata_specmap
1308  read(21, fmt=*) d_dummy, specmap_zsl(n)
1309 end do
1310 
1311 close(21, status='keep')
1312 
1313 #endif
1314 
1315 !-------- Determination of the geothermal heat flux --------
1316 
1317 #if Q_GEO_MOD==1
1318 
1319 ! ------ Constant value
1320 
1321 do i=0, imax
1322 do j=0, jmax
1323  q_geo(j,i) = q_geo *1.0e-03_dp ! mW/m2 -> W/m2
1324 end do
1325 end do
1326 
1327 #elif Q_GEO_MOD==2
1328 
1329 ! ------ Read data from file
1330 
1331 open(21, iostat=ios, &
1332  file=inpath//'/'//trim(ch_domain_short)//'/'//q_geo_file, &
1333  recl=8192, status='old')
1334 
1335 if (ios /= 0) stop ' sico_init: Error when opening the qgeo file!'
1336 
1337 do n=1, 6; read(21, fmt='(a)') ch_dummy; end do
1338 
1339 do j=jmax, 0, -1
1340  read(21, fmt=*) (q_geo(j,i), i=0,imax)
1341 end do
1342 
1343 close(21, status='keep')
1344 
1345 do i=0, imax
1346 do j=0, jmax
1347  q_geo(j,i) = q_geo(j,i) *1.0e-03_dp ! mW/m2 -> W/m2
1348 end do
1349 end do
1350 
1351 #endif
1352 
1353 !-------- Reading of tabulated kei function--------
1354 
1355 #if (REBOUND==0 || REBOUND==1)
1356 
1357 kei = 0.0_dp; n_data_kei = 0; kei_r_max = 0.0_dp; kei_r_incr = 0.0_dp
1358  ! dummy values
1359 #elif (REBOUND==2)
1360 
1361 call read_kei()
1362 
1363 #endif
1364 
1365 !-------- Determination of the time lag
1366 ! of the relaxing asthenosphere --------
1367 
1368 #if (REBOUND==1 || REBOUND==2)
1369 
1370 #if (TIME_LAG_MOD==1)
1371 
1372 time_lag_asth = time_lag*year_sec ! a -> s
1373 
1374 #elif (TIME_LAG_MOD==2)
1375 
1376 open(29, iostat=ios, &
1377  file=inpath//'/'//trim(ch_domain_short)//'/'//time_lag_file, &
1378  recl=8192, status='old')
1379 
1380 if (ios /= 0) stop ' sico_init: Error when opening the time-lag file!'
1381 
1382 do n=1, 6; read(29, fmt='(a)') ch_dummy; end do
1383 
1384 do j=jmax, 0, -1
1385  read(29, fmt=*) (time_lag_asth(j,i), i=0,imax)
1386 end do
1387 
1388 close(29, status='keep')
1389 
1390 time_lag_asth = time_lag_asth*year_sec ! a -> s
1391 
1392 #endif
1393 
1394 #elif (REBOUND==0)
1395 
1396 time_lag_asth = 0.0_dp ! dummy values
1397 
1398 #endif
1399 
1400 !-------- Determination of the flexural rigidity of the lithosphere --------
1401 
1402 #if (REBOUND==2)
1403 
1404 #if (FLEX_RIG_MOD==1)
1405 
1406 flex_rig_lith = flex_rig
1407 
1408 #elif (FLEX_RIG_MOD==2)
1409 
1410 open(29, iostat=ios, &
1411  file=inpath//'/'//trim(ch_domain_short)//'/'//flex_rig_file, &
1412  recl=8192, status='old')
1413 
1414 if (ios /= 0) stop ' sico_init: Error when opening the flex-rig file!'
1415 
1416 do n=1, 6; read(29, fmt='(a)') ch_dummy; end do
1417 
1418 do j=jmax, 0, -1
1419  read(29, fmt=*) (flex_rig_lith(j,i), i=0,imax)
1420 end do
1421 
1422 close(29, status='keep')
1423 
1424 #endif
1425 
1426 #elif (REBOUND==0 || REBOUND==1)
1427 
1428 flex_rig_lith = 0.0_dp ! dummy values
1429 
1430 #endif
1431 
1432 !-------- Definition of initial values --------
1433 
1434 ! ------ Present topography
1435 
1436 #if (ANF_DAT==1)
1437 
1438 call topography1(dxi, deta)
1439 
1440 z_sl = -1.11e+11_dp ! dummy value for initial call of subroutine boundary
1441 
1442 call boundary(time_init, dtime, dxi, deta, &
1443  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1444 
1445 q_bm = 0.0_dp
1446 q_tld = 0.0_dp
1447 q_b_tot = 0.0_dp
1448 
1449 p_b_w = 0.0_dp
1450 h_w = 0.0_dp
1451 
1452 call init_temp_water_age()
1453 
1454 call calc_enhance(time_init)
1455 
1456 ! ------ Ice-free, relaxed bedrock
1457 
1458 #elif (ANF_DAT==2)
1459 
1460 call topography2(dxi, deta)
1461 
1462 z_sl = -1.11e+11_dp ! dummy value for initial call of subroutine boundary
1463 
1464 call boundary(time_init, dtime, dxi, deta, &
1465  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1466 
1467 q_bm = 0.0_dp
1468 q_tld = 0.0_dp
1469 q_b_tot = 0.0_dp
1470 
1471 p_b_w = 0.0_dp
1472 h_w = 0.0_dp
1473 
1474 call init_temp_water_age()
1475 
1476 call calc_enhance(time_init)
1477 
1478 ! ------ Read initial state from output of previous simulation
1479 
1480 #elif (ANF_DAT==3)
1481 
1482 #if (NETCDF==1)
1483 call topography3(dxi, deta, z_sl, anfdatname)
1484 #elif (NETCDF==2)
1485 call topography3_nc(dxi, deta, z_sl, anfdatname)
1486 #else
1487 stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
1488 #endif
1489 
1490 call boundary(time_init, dtime, dxi, deta, &
1491  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1492 
1493 q_b_tot = q_bm + q_tld
1494 
1495 #endif
1496 
1497 !-------- Grounding line and calving front flags --------
1498 
1499 flag_grounding_line_1 = .false.
1500 flag_grounding_line_2 = .false.
1501 
1502 flag_calving_front_1 = .false.
1503 flag_calving_front_2 = .false.
1504 
1505 #if (MARGIN==1 || MARGIN==2)
1506 
1507 continue
1508 
1509 #elif (MARGIN==3)
1510 
1511 do i=1, imax-1
1512 do j=1, jmax-1
1513 
1514  if ( (maske(j,i)==0_i2b) & ! grounded ice
1515  .and. &
1516  ( (maske(j,i+1)==3_i2b) & ! with
1517  .or.(maske(j,i-1)==3_i2b) & ! one
1518  .or.(maske(j+1,i)==3_i2b) & ! neighbouring
1519  .or.(maske(j-1,i)==3_i2b) ) & ! floating ice point
1520  ) &
1521  flag_grounding_line_1(j,i) = .true.
1522 
1523  if ( (maske(j,i)==3_i2b) & ! floating ice
1524  .and. &
1525  ( (maske(j,i+1)==0_i2b) & ! with
1526  .or.(maske(j,i-1)==0_i2b) & ! one
1527  .or.(maske(j+1,i)==0_i2b) & ! neighbouring
1528  .or.(maske(j-1,i)==0_i2b) ) & ! grounded ice point
1529  ) &
1530  flag_grounding_line_2(j,i) = .true.
1531 
1532  if ( (maske(j,i)==3_i2b) & ! floating ice
1533  .and. &
1534  ( (maske(j,i+1)==2_i2b) & ! with
1535  .or.(maske(j,i-1)==2_i2b) & ! one
1536  .or.(maske(j+1,i)==2_i2b) & ! neighbouring
1537  .or.(maske(j-1,i)==2_i2b) ) & ! ocean point
1538  ) &
1539  flag_calving_front_1(j,i) = .true.
1540 
1541  if ( (maske(j,i)==2_i2b) & ! ocean
1542  .and. &
1543  ( (maske(j,i+1)==3_i2b) & ! with
1544  .or.(maske(j,i-1)==3_i2b) & ! one
1545  .or.(maske(j+1,i)==3_i2b) & ! neighbouring
1546  .or.(maske(j-1,i)==3_i2b) ) & ! floating ice point
1547  ) &
1548  flag_calving_front_2(j,i) = .true.
1549 
1550 end do
1551 end do
1552 
1553 do i=1, imax-1
1554 
1555  j=0
1556  if ( (maske(j,i)==2_i2b) & ! ocean
1557  .and. (maske(j+1,i)==3_i2b) & ! with one neighbouring
1558  ) & ! floating ice point
1559  flag_calving_front_2(j,i) = .true.
1560 
1561  j=jmax
1562  if ( (maske(j,i)==2_i2b) & ! ocean
1563  .and. (maske(j-1,i)==3_i2b) & ! with one neighbouring
1564  ) & ! floating ice point
1565  flag_calving_front_2(j,i) = .true.
1566 
1567 end do
1568 
1569 do j=1, jmax-1
1570 
1571  i=0
1572  if ( (maske(j,i)==2_i2b) & ! ocean
1573  .and. (maske(j,i+1)==3_i2b) & ! with one neighbouring
1574  ) & ! floating ice point
1575  flag_calving_front_2(j,i) = .true.
1576 
1577  i=imax
1578  if ( (maske(j,i)==2_i2b) & ! ocean
1579  .and. (maske(j,i-1)==3_i2b) & ! with one neighbouring
1580  ) & ! floating ice point
1581  flag_calving_front_2(j,i) = .true.
1582 
1583 end do
1584 
1585 #else
1586 stop ' sico_init: MARGIN must be either 1, 2 or 3!'
1587 #endif
1588 
1589 !-------- Rock/sediment mask: smoothing over the transition length --------
1590 
1591 #if ( defined(SEDI_SLIDE) && SEDI_SLIDE==2 && defined(TRANS_LENGTH_SL) )
1592 
1593 do i=0, imax
1594 do j=0, jmax
1595  if (maske_sedi(j,i) /= 7) then ! no soft sediment
1596  r_mask_sedi(j,i) = 0.0_dp
1597  else ! maske_sedi(j,i) == 7, soft sediment
1598  r_mask_sedi(j,i) = 1.0_dp
1599  end if
1600 end do
1601 end do
1602 
1603 if (trans_length_sl > eps) then
1604 
1605  if ((p_weert /= p_weert_sedi).or.(q_weert /= q_weert_sedi)) then
1606  write(6, fmt='(a)') ' sico_init: Specifying a transition length'
1607  write(6, fmt='(a)') ' TRANS_LENGTH_SL > 0 for the basal sliding'
1608  write(6, fmt='(a)') ' regimes requires P_WEERT==P_WEERT_SEDI'
1609  write(6, fmt='(a)') ' and Q_WEERT==Q_WEERT_SEDI!'
1610  stop
1611  end if
1612 
1613  transition_length_sliding = trans_length_sl *1.0e+03_dp ! km -> m
1614 
1615  quasi_time = 0.25_dp*transition_length_sliding**2
1616  ! Time scale resulting from the fundamental solution
1617  ! of the diffusion equation for diffusitivity = 1
1618  d_quasi_time = 0.1_dp*quasi_time
1619  ! First guess for the time step
1620 
1621  do
1622  if ( d_quasi_time < 0.1_dp*min(dxi,deta)**2 ) exit
1623  ! CFL condition fulfilled (with some leeway)
1624  d_quasi_time = 0.5_dp*d_quasi_time
1625  end do
1626 
1627  m = nint(quasi_time/d_quasi_time)
1628 
1629  do n=1, m
1630 
1631  r_mask_sedi_new = r_mask_sedi
1632 
1633  do i=1, imax-1
1634  do j=1, jmax-1
1635  r_mask_sedi_new(j,i) = r_mask_sedi(j,i) &
1636  + d_quasi_time/dxi**2 &
1637  * (r_mask_sedi(j,i+1)-2.0_dp*r_mask_sedi(j,i)+r_mask_sedi(j,i-1)) &
1638  + d_quasi_time/deta**2 &
1639  * (r_mask_sedi(j+1,i)-2.0_dp*r_mask_sedi(j,i)+r_mask_sedi(j-1,i))
1640  ! Solving the diffusion equation for diffusitivity = 1
1641  end do
1642  end do
1643 
1644  r_mask_sedi = r_mask_sedi_new
1645 
1646  end do
1647 
1648 end if
1649 
1650 #endif
1651 
1652 !-------- Distance between grid points with delta_i=ir, delta_j=jr --------
1653 
1654 #if (GRID==0 || GRID==1)
1655 
1656 do ir=-imax, imax
1657 do jr=-jmax, jmax
1658  dist_dxdy(jr,ir) = sqrt( (real(ir,dp)*dxi)**2 + (real(jr,dp)*deta)**2 )
1659  ! distortion due to stereographic projection not accounted for
1660 end do
1661 end do
1662 
1663 #endif
1664 
1665 !-------- Initial velocities --------
1666 
1667 call calc_temp_melt()
1668 
1669 #if (DYNAMICS==1 || DYNAMICS==2)
1670 
1671 call calc_vxy_b_sia(time, z_sl)
1672 call calc_vxy_sia(dzeta_c, dzeta_t)
1673 
1674 #if (MARGIN==3 || DYNAMICS==2)
1675 call calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t, ii, jj, nn)
1676 #endif
1677 
1678 call calc_vz_sia(dxi, deta, dzeta_c, dzeta_t)
1679 
1680 #if (MARGIN==3)
1681 call calc_vz_ssa(z_sl, dxi, deta, dzeta_c)
1682 #endif
1683 
1684 #elif (DYNAMICS==0)
1685 
1686 call calc_vxy_static()
1687 call calc_vz_static()
1688 
1689 #else
1690  stop ' sico_init: DYNAMICS must be either 0, 1 or 2!'
1691 #endif
1692 
1693 call calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
1694 
1695 !-------- Initial enthalpies --------
1696 
1697 #if (CALCMOD==0 || CALCMOD==-1)
1698 
1699 do i=0, imax
1700 do j=0, jmax
1701 
1702  do kc=0, kcmax
1703  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
1704  end do
1705 
1706  do kt=0, ktmax
1707  enth_t(kt,j,i) = enth_c(0,j,i)
1708  end do
1709 
1710 end do
1711 end do
1712 
1713 #elif (CALCMOD==1)
1714 
1715 do i=0, imax
1716 do j=0, jmax
1717 
1718  do kc=0, kcmax
1719  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
1720  end do
1721 
1722  if ( (maske(j,i) == 0_i2b).and.(n_cts(j,i) == 1_i2b) ) then
1723  do kt=0, ktmax
1724  enth_t(kt,j,i) = enth_fct_temp_omega(temp_t_m(kt,j,i), omega_t(kt,j,i))
1725  end do
1726  else
1727  do kt=0, ktmax
1728  enth_t(kt,j,i) = enth_c(0,j,i)
1729  end do
1730  end if
1731 
1732 end do
1733 end do
1734 
1735 #elif (CALCMOD==2 || CALCMOD==3)
1736 
1737 do i=0, imax
1738 do j=0, jmax
1739 
1740  do kc=0, kcmax
1741  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), omega_c(kc,j,i))
1742  end do
1743 
1744  do kt=0, ktmax
1745  enth_t(kt,j,i) = enth_c(0,j,i)
1746  end do
1747 
1748 end do
1749 end do
1750 
1751 #else
1752 
1753 stop ' sico_init: Parameter CALCMOD must be either -1, 0, 1, 2 or 3!'
1754 
1755 #endif
1756 
1757 !-------- Initialize time-series files --------
1758 
1759 ! ------ Time-series file for the ice sheet on the whole
1760 
1761 open(12, iostat=ios, &
1762  file=outpath//'/'//trim(runname)//'.ser', &
1763  status='new')
1764 if (ios /= 0) &
1765  stop ' sico_init: Error when opening the ser file!'
1766 
1767 if (forcing_flag == 1) then
1768 
1769  write(12,1102)
1770  write(12,1103)
1771 
1772 #if ( !defined(OUTSER_V_A) || OUTSER_V_A==1 )
1773 
1774  1102 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
1775  ' H_max(m) zs_max(m) V_g(m^3)', &
1776  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1777  ' A_g(m^2) A_t(m^2) V_bm(m^3/a)', &
1778  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1779  1103 format('----------------------------------------------------', &
1780  '---------------------------------------')
1781 
1782 #elif OUTSER_V_A==2
1783 
1784  1102 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
1785  ' V(m^3) V_g(m^3) V_f(m^3)', &
1786  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1787  ' H_max(m) zs_max(m)', &
1788  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1789  ' A_t(m^2) V_bm(m^3/a)', &
1790  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1791  1103 format('----------------------------------------------------', &
1792  '---------------------------------------')
1793 
1794 #else
1795  stop ' sico_init: OUTSER_V_A must be either 1 or 2!'
1796 #endif
1797 
1798 else if (forcing_flag == 2) then
1799 
1800  write(12,1112)
1801  write(12,1113)
1802 
1803 #if ( !defined(OUTSER_V_A) || OUTSER_V_A==1 )
1804 
1805  1112 format(' t(a) glac_ind(1) z_sl(m)',/, &
1806  ' H_max(m) zs_max(m) V_g(m^3)', &
1807  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1808  ' A_g(m^2) A_t(m^2) V_bm(m^3/a)', &
1809  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1810  1113 format('----------------------------------------------------', &
1811  '---------------------------------------')
1812 
1813 #elif OUTSER_V_A==2
1814 
1815  1112 format(' t(a) glac_ind(1) z_sl(m)',/, &
1816  ' V(m^3) V_g(m^3) V_f(m^3)', &
1817  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1818  ' H_max(m) zs_max(m)', &
1819  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1820  ' A_t(m^2) V_bm(m^3/a)', &
1821  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1822  1113 format('----------------------------------------------------', &
1823  '---------------------------------------')
1824 
1825 #else
1826  stop ' sico_init: OUTSER_V_A must be either 1 or 2!'
1827 #endif
1828 
1829 else if (forcing_flag == 3) then
1830 
1831  write(12,1122)
1832  write(12,1123)
1833 
1834 #if ( !defined(OUTSER_V_A) || OUTSER_V_A==1 )
1835 
1836  1122 format(' t(a) (Dummy)(1) z_sl(m)',/, &
1837  ' H_max(m) zs_max(m) V_g(m^3)', &
1838  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1839  ' A_g(m^2) A_t(m^2) V_bm(m^3/a)', &
1840  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1841  1123 format('----------------------------------------------------', &
1842  '---------------------------------------')
1843 
1844 #elif OUTSER_V_A==2
1845 
1846  1122 format(' t(a) (Dummy)(1) z_sl(m)',/, &
1847  ' V(m^3) V_g(m^3) V_f(m^3)', &
1848  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1849  ' H_max(m) zs_max(m)', &
1850  ' V_t(m^3) V_fw(m^3/a) V_sle(m)',/, &
1851  ' A_t(m^2) V_bm(m^3/a)', &
1852  ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1853  1123 format('----------------------------------------------------', &
1854  '---------------------------------------')
1855 
1856 #else
1857  stop ' sico_init: OUTSER_V_A must be either 1 or 2!'
1858 #endif
1859 
1860 end if
1861 
1862 ! ------ Time-series file for the ice-thickness field
1863 
1864 #if OUTSER==2
1865 
1866 open(13, iostat=ios, &
1867  file=outpath//'/'//trim(runname)//'.thk', &
1868  status='new')
1869 if (ios /= 0) stop ' sico_init: Error when opening the thk file!'
1870 
1871 if (forcing_flag == 1) then
1872 
1873  write(13,1104)
1874  write(13,1105)
1875 
1876  1104 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
1877  ' H(m): JMAX+1 records [j = JMAX (-1) 0] with',/, &
1878  ' IMAX+1 values [i = 0 (1) IMAX] in each record')
1879  1105 format('----------------------------------------------------')
1880 
1881 else if (forcing_flag == 2) then
1882 
1883  write(13,1114)
1884  write(13,1115)
1885 
1886  1114 format(' t(a) glac_ind(1) z_sl(m)',/, &
1887  ' H(m): JMAX+1 records [j = JMAX (-1) 0] with',/, &
1888  ' IMAX+1 values [i = 0 (1) IMAX] in each record')
1889  1115 format('----------------------------------------------------')
1890 
1891 else if (forcing_flag == 3) then
1892 
1893  write(13,1124)
1894  write(13,1125)
1895 
1896  1124 format(' t(a) (Dummy)(1) z_sl(m)',/, &
1897  ' H(m): JMAX+1 records [j = JMAX (-1) 0] with',/, &
1898  ' IMAX+1 values [i = 0 (1) IMAX] in each record')
1899  1125 format('----------------------------------------------------')
1900 
1901 end if
1902 
1903 #endif
1904 
1905 ! ------ Time-series file for deep boreholes
1906 
1907 #if OUTSER==3
1908 
1909 n_core = 6 ! Vostok, Dome A, Dome C, Dome F, Kohnen, Byrd
1910 
1911 allocate(lambda_core(n_core), phi_core(n_core), &
1912  x_core(n_core), y_core(n_core))
1913 
1914 phi_core(1) = -78.467_dp *pi_180 ! Geographical position of Vostok,
1915 lambda_core(1) = 106.8_dp *pi_180 ! conversion deg -> rad
1916 call num_coord(lambda_core(1), phi_core(1), x_core(1), y_core(1))
1917 
1918 phi_core(2) = -80.367_dp *pi_180 ! Geographical position of Dome A,
1919 lambda_core(2) = 77.35_dp *pi_180 ! conversion deg -> rad
1920 call num_coord(lambda_core(2), phi_core(2), x_core(2), y_core(2))
1921 
1922 phi_core(3) = -75.1_dp *pi_180 ! Geographical position of Dome C,
1923 lambda_core(3) = 123.4_dp *pi_180 ! conversion deg -> rad
1924 call num_coord(lambda_core(3), phi_core(3), x_core(3), y_core(3))
1925 
1926 phi_core(4) = -77.317_dp *pi_180 ! Geographical position of Dome F,
1927 lambda_core(4) = 39.7_dp *pi_180 ! conversion deg -> rad
1928 call num_coord(lambda_core(4), phi_core(4), x_core(4), y_core(4))
1929 
1930 phi_core(5) = -75.0_dp *pi_180 ! Geographical position of Kohnen,
1931 lambda_core(5) = 0.067_dp *pi_180 ! conversion deg -> rad
1932 call num_coord(lambda_core(5), phi_core(5), x_core(5), y_core(5))
1933 
1934 phi_core(6) = -80.017_dp *pi_180 ! Geographical position of Byrd,
1935 lambda_core(6) = -119.517_dp *pi_180 ! conversion deg -> rad
1936 call num_coord(lambda_core(6), phi_core(6), x_core(6), y_core(6))
1937 
1938 open(14, iostat=ios, &
1939  file=outpath//'/'//trim(runname)//'.core', &
1940  status='new')
1941 if (ios /= 0) stop ' sico_init: Error when opening the core file!'
1942 
1943 if (forcing_flag == 1) then
1944 
1945  write(14,1106)
1946  write(14,1107)
1947 
1948  1106 format(' t(a) D_Ts(C) z_sl(m)',/, &
1949  ' H_Vo(m) H_DA(m) H_DC(m)', &
1950  ' H_DF(m) H_Ko(m) H_By(m)',/, &
1951  ' v_Vo(m/a) v_DA(m/a) v_DC(m/a)', &
1952  ' v_DF(m/a) v_Ko(m/a) v_By(m/a)',/, &
1953  ' T_Vo(C) T_DA(C) T_DC(C)', &
1954  ' T_DF(C) T_Ko(C) T_By(C)',/, &
1955  ' Rb_Vo(W/m2) Rb_DA(W/m2) Rb_DC(W/m2)', &
1956  ' Rb_DF(W/m2) Rb_Ko(W/m2) Rb_By(W/m2)')
1957  1107 format('----------------------------------------------------', &
1958  '---------------------------------------')
1959 
1960 else if (forcing_flag == 2) then
1961 
1962  write(14,1116)
1963  write(14,1117)
1964 
1965  1116 format(' t(a) glac_ind(1) z_sl(m)',/, &
1966  ' H_Vo(m) H_DA(m) H_DC(m)', &
1967  ' H_DF(m) H_Ko(m) H_By(m)',/, &
1968  ' v_Vo(m/a) v_DA(m/a) v_DC(m/a)', &
1969  ' v_DF(m/a) v_Ko(m/a) v_By(m/a)',/, &
1970  ' T_Vo(C) T_DA(C) T_DC(C)', &
1971  ' T_DF(C) T_Ko(C) T_By(C)',/, &
1972  ' Rb_Vo(W/m2) Rb_DA(W/m2) Rb_DC(W/m2)', &
1973  ' Rb_DF(W/m2) Rb_Ko(W/m2) Rb_By(W/m2)')
1974  1117 format('----------------------------------------------------', &
1975  '---------------------------------------')
1976 
1977 else if (forcing_flag == 3) then
1978 
1979  write(14,1126)
1980  write(14,1127)
1981 
1982  1126 format(' t(a) (Dummy)(1) z_sl(m)',/, &
1983  ' H_Vo(m) H_DA(m) H_DC(m)', &
1984  ' H_DF(m) H_Ko(m) H_By(m)',/, &
1985  ' v_Vo(m/a) v_DA(m/a) v_DC(m/a)', &
1986  ' v_DF(m/a) v_Ko(m/a) v_By(m/a)',/, &
1987  ' T_Vo(C) T_DA(C) T_DC(C)', &
1988  ' T_DF(C) T_Ko(C) T_By(C)',/, &
1989  ' Rb_Vo(W/m2) Rb_DA(W/m2) Rb_DC(W/m2)', &
1990  ' Rb_DF(W/m2) Rb_Ko(W/m2) Rb_By(W/m2)')
1991  1127 format('----------------------------------------------------', &
1992  '---------------------------------------')
1993 
1994 end if
1995 
1996 #endif
1997 
1998 !-------- Output of the initial state --------
1999 
2000 #if OUTPUT==1
2001 
2002 #if ERGDAT==0
2003  flag_3d_output = .false.
2004 #elif ERGDAT==1
2005  flag_3d_output = .true.
2006 #else
2007  stop ' sico_init: ERGDAT must be either 0 or 1!'
2008 #endif
2009 
2010 #if NETCDF==1
2011  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
2012  flag_3d_output, ndat2d, ndat3d)
2013 #elif NETCDF==2
2014  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
2015  flag_3d_output, ndat2d, ndat3d)
2016 #else
2017  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
2018 #endif
2019 
2020 #elif OUTPUT==2
2021 
2022 if (time_output(1) <= time_init+eps) then
2023 
2024 #if ERGDAT==0
2025  flag_3d_output = .false.
2026 #elif ERGDAT==1
2027  flag_3d_output = .true.
2028 #else
2029  stop ' sico_init: ERGDAT must be either 0 or 1!'
2030 #endif
2031 
2032 #if NETCDF==1
2033  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
2034  flag_3d_output, ndat2d, ndat3d)
2035 #elif NETCDF==2
2036  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
2037  flag_3d_output, ndat2d, ndat3d)
2038 #else
2039  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
2040 #endif
2041 
2042 end if
2043 
2044 #elif OUTPUT==3
2045 
2046  flag_3d_output = .false.
2047 
2048 #if NETCDF==1
2049  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
2050  flag_3d_output, ndat2d, ndat3d)
2051 #elif NETCDF==2
2052  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
2053  flag_3d_output, ndat2d, ndat3d)
2054 #else
2055  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
2056 #endif
2057 
2058 if (time_output(1) <= time_init+eps) then
2059 
2060  flag_3d_output = .true.
2061 
2062 #if NETCDF==1
2063  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
2064  flag_3d_output, ndat2d, ndat3d)
2065 #elif NETCDF==2
2066  call output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
2067  flag_3d_output, ndat2d, ndat3d)
2068 #else
2069  stop ' sico_init: Parameter NETCDF must be either 1 or 2!'
2070 #endif
2071 
2072 end if
2073 
2074 #else
2075  stop ' sico_init: OUTPUT must be either 1, 2 or 3!'
2076 #endif
2077 
2078 call output2(time_init, dxi, deta, delta_ts, glac_index, z_sl)
2079 #if OUTSER==2
2080 call output3(time_init, delta_ts, glac_index, z_sl)
2081 #endif
2082 #if OUTSER==3
2083 call output4(time_init, dxi, deta, delta_ts, glac_index, z_sl)
2084 #endif
2085 
2086 end subroutine sico_init
2087 !
nc_check
NetCDF error capturing.
Definition: nc_check.F90:35
calc_vz_ssa
subroutine calc_vz_ssa(z_sl, dxi, deta, dzeta_c)
Computation of the vertical velocity vz in the shallow shelf approximation.
Definition: calc_vz_ssa.F90:35
calc_vxy_static
subroutine calc_vxy_static()
Computation of the horizontal velocity vx, vy, the horizontal volume flux qx, qy etc. for static ice.
Definition: calc_vxy_static.F90:37
calc_vxy_sia
subroutine calc_vxy_sia(dzeta_c, dzeta_t)
Computation of the shear stresses txz, tyz, the effective shear stress sigma, the depth-averaged flui...
Definition: calc_vxy_sia.F90:41
topography1
subroutine topography1(dxi, deta)
Definition of the initial surface and bedrock topography (including gradients) and of the horizontal ...
Definition: topography1.F90:39
ice_material_quantities::ice_mat_eqs_pars
subroutine, public ice_mat_eqs_pars(RF_table, R_T_val, KAPPA_table, C_table, n_temp_min, n_temp_max, RHO_I_val, RHO_C_val, KAPPA_C_val, C_C_val)
Setting of required physical parameters.
Definition: ice_material_quantities.F90:74
output2
subroutine output2(time, dxi, deta, delta_ts, glac_index, z_sl)
Writing of time-series data on file in ASCII format.
Definition: output2.F90:35
phys_para
subroutine phys_para()
Reading of physical parameters.
Definition: phys_para.F90:35
sico_types
Declarations of kind types for SICOPOLIS.
Definition: sico_types.F90:35
topography3
subroutine topography3(dxi, deta, z_sl, anfdatname)
Definition of the initial surface and bedrock topography (including gradients) and of the horizontal ...
Definition: topography3.F90:39
output3
subroutine output3(time, delta_ts, glac_index, z_sl)
Writing of time-series data of the ice thickness field on file in ASCII format.
Definition: output3.F90:38
num_coord
subroutine num_coord(lambda_val, phi_val, x_val, y_val)
Computation of position (x,y) in the numerical domain for longitude lambda and latitude phi (in rad)...
Definition: num_coord.F90:37
topography3_nc
subroutine topography3_nc(dxi, deta, z_sl, anfdatname)
Definition of the initial surface and bedrock topography (including gradients) and of the horizontal ...
Definition: topography3_nc.F90:39
calc_temp_melt
subroutine calc_temp_melt()
Computation of the melting temperatures.
Definition: calc_temp_melt.F90:35
output_nc
subroutine output_nc(runname, time, delta_ts, glac_index, z_sl, flag_3d_output, ndat2d, ndat3d)
Writing of time-slice files in NetCDF format.
Definition: output_nc.F90:35
calc_vxy_ssa
subroutine calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t, ii, jj, nn)
Computation of the horizontal velocity vx, vy, the horizontal volume flux qx, qy and the flux across ...
Definition: calc_vxy_ssa.F90:39
nc_check::check
subroutine check(status, ch_calling_routine)
NetCDF error capturing.
Definition: nc_check.F90:45
sico_vars
Declarations of global variables for SICOPOLIS (for the ANT domain).
Definition: sico_vars.F90:35
topography2
subroutine topography2(dxi, deta)
Definition of the initial surface and bedrock topography (including gradients) and of the horizontal ...
Definition: topography2.F90:39
calc_vz_static
subroutine calc_vz_static()
Computation of the vertical velocity vz for static ice.
Definition: calc_vz_static.F90:35
calc_dxyz
subroutine calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
Computation of all components of the strain-rate tensor, the full effective strain rate and the shear...
Definition: calc_dxyz.F90:37
calc_enhance
subroutine calc_enhance(time)
Computation of the flow enhancement factor.
Definition: calc_enhance.F90:35
enth_temp_omega::calc_c_int_table
subroutine, public calc_c_int_table(c_table, n_tmp_min, n_tmp_max, L_val)
Computation of the temperature integral of the specific heat of ice as a table (c_int_table). Further, definition of the latent heat of ice.
Definition: enth_temp_omega.F90:84
calc_vxy_b_sia
subroutine calc_vxy_b_sia(time, z_sl)
Computation of the basal horizontal velocity vx_b, vy_b in the shallow ice approximation.
Definition: calc_vxy_b_sia.F90:37
sico_init
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)
Initialisations for SICOPOLIS.
Definition: sico_init.F90:35
boundary
subroutine boundary(time, dtime, dxi, deta, delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
Computation of the surface temperature (must be less than 0 deg C!) and of the accumulation-ablation ...
Definition: boundary.F90:37
read_kei
subroutine read_kei()
Reading of the tabulated kei function.
Definition: read_kei.F90:35
init_temp_water_age
subroutine init_temp_water_age()
Initial temperature, water content and age.
Definition: init_temp_water_age.F90:35
ice_material_quantities
Material quantities of ice: Rate factor, heat conductivity, specific heat (heat capacity).
Definition: ice_material_quantities.F90:37
output4
subroutine output4(time, dxi, deta, delta_ts, glac_index, z_sl)
Writing of time-series data of the deep ice cores on file in ASCII format.
Definition: output4.F90:35
read_target_topo_nc
subroutine read_target_topo_nc(target_topo_dat_name)
Reading of the target-topography file (in NetCDF format).
Definition: read_target_topo_nc.F90:35
sico_variables
Declarations of global variables for SICOPOLIS.
Definition: sico_variables.F90:35
enth_temp_omega::calc_c_int_inv_table
subroutine, public calc_c_int_inv_table()
Computation of the inverse of the temperature integral of the specific heat of ice as a table (c_int_...
Definition: enth_temp_omega.F90:141
enth_temp_omega
Conversion from temperature (temp) and water content (omega) to enthalpy (enth) and vice versa...
Definition: enth_temp_omega.F90:37
calc_vz_sia
subroutine calc_vz_sia(dxi, deta, dzeta_c, dzeta_t)
Computation of the vertical velocity vz in the shallow ice approximation.
Definition: calc_vz_sia.F90:35
output1
subroutine output1(runname, time, delta_ts, glac_index, z_sl, flag_3d_output, ndat2d, ndat3d)
Writing of time-slice files in binary format.
Definition: output1.F90:35
enth_temp_omega::enth_fct_temp_omega
real(dp) function, public enth_fct_temp_omega(temp_val, omega_val)
Enthalpy as a function of temperature and water content.
Definition: enth_temp_omega.F90:259
Generated on Wed Mar 2 2016 11:26:49 for SICOPOLIS V3.2 by   doxygen 1.8.5