SICOPOLIS V5-dev  Revision 1368
sico_init_m.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Module : s i c o _ i n i t _ m
4 !
5 !> @file
6 !!
7 !! Initialisations for SICOPOLIS.
8 !!
9 !! @section Copyright
10 !!
11 !! Copyright 2009-2018 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 module sico_init_m
36 
37  use sico_types_m
39  use sico_vars_m
40  use error_m
41 
42  implicit none
43 
44  public
45 
46 contains
47 
48 !-------------------------------------------------------------------------------
49 !> Main routine of sico_init_m: Initialisations for SICOPOLIS.
50 !<------------------------------------------------------------------------------
51 subroutine sico_init(delta_ts, glac_index, &
52  mean_accum, &
53  dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, &
54  time, time_init, time_end, time_output, &
55  dxi, deta, dzeta_c, dzeta_t, dzeta_r, &
56  z_sl, dzsl_dtau, z_mar, &
57  ii, jj, nn, &
58  ndat2d, ndat3d, n_output, &
59  runname)
60 
61  use compare_float_m
65 
66  use read_m, only : read_phys_para, read_kei
67 
68  use boundary_m
70  use calc_enhance_m
72  use calc_vxy_m
73  use calc_vz_m
74  use calc_dxyz_m
76 
77  use output_m
78 
79 implicit none
80 
81 integer(i4b), intent(out) :: ii((imax+1)*(jmax+1)), &
82  jj((imax+1)*(jmax+1)), &
83  nn(0:jmax,0:imax)
84 integer(i4b), intent(out) :: ndat2d, ndat3d
85 integer(i4b), intent(out) :: n_output
86 real(dp), intent(out) :: delta_ts, glac_index
87 real(dp), intent(out) :: mean_accum
88 real(dp), intent(out) :: dtime, dtime_temp, dtime_wss, &
89  dtime_out, dtime_ser
90 real(dp), intent(out) :: time, time_init, time_end, time_output(100)
91 real(dp), intent(out) :: dxi, deta, dzeta_c, dzeta_t, dzeta_r
92 real(dp), intent(out) :: z_sl, dzsl_dtau, z_mar
93 character(len=100), intent(out) :: runname
94 
95 integer(i4b) :: i, j, kc, kt, kr, m, n, ir, jr
96 integer(i4b) :: ios
97 integer(i4b) :: ierr
98 integer(i1b), dimension(0:JMAX,0:IMAX) :: maske_ref
99 real(dp) :: dtime0, dtime_temp0, dtime_wss0, dtime_out0, dtime_ser0
100 real(dp) :: time_init0, time_end0, time_output0(100)
101 real(dp) :: d_dummy
102 character(len=100) :: anfdatname
103 character(len=256) :: filename_with_path
104 character(len=256) :: shell_command
105 character :: ch_dummy
106 logical :: flag_init_output, flag_3d_output
107 
108 character(len=64), parameter :: fmt1 = '(a)', &
109  fmt2 = '(a,i4)', &
110  fmt2a = '(a,i0)', &
111  fmt3 = '(a,es12.4)'
112 
113 character(len= 8) :: ch_imax
114 character(len=128) :: fmt4
115 
116 write(ch_imax, fmt='(i8)') imax
117 write(fmt4, fmt='(a)') '('//trim(adjustl(ch_imax))//'(i1),i1)'
118 
119 write(unit=6, fmt='(a)') ' '
120 write(unit=6, fmt='(a)') ' -------- sico_init --------'
121 write(unit=6, fmt='(a)') ' '
122 
123 !-------- Name of the computational domain --------
124 
125 #if (defined(ANT))
126 ch_domain_long = 'Antarctica'
127 ch_domain_short = 'ant'
128 
129 #elif (defined(ASF))
130 ch_domain_long = 'Austfonna'
131 ch_domain_short = 'asf'
132 
133 #elif (defined(EMTP2SGE))
134 ch_domain_long = 'EISMINT Phase 2 Simplified Geometry Experiment'
135 ch_domain_short = 'emtp2sge'
136 
137 #elif (defined(GRL))
138 ch_domain_long = 'Greenland'
139 ch_domain_short = 'grl'
140 
141 #elif (defined(NHEM))
142 ch_domain_long = 'Northern hemisphere'
143 ch_domain_short = 'nhem'
144 
145 #elif (defined(SCAND))
146 ch_domain_long = 'Scandinavia and Eurasia'
147 ch_domain_short = 'scand'
148 
149 #elif (defined(TIBET))
150 ch_domain_long = 'Tibet'
151 ch_domain_short = 'tibet'
152 
153 #elif (defined(NMARS))
154 ch_domain_long = 'North polar cap of Mars'
155 ch_domain_short = 'nmars'
156 
157 #elif (defined(SMARS))
158 ch_domain_long = 'South polar cap of Mars'
159 ch_domain_short = 'smars'
160 
161 #elif (defined(XYZ))
162 ch_domain_long = 'XYZ'
163 ch_domain_short = 'xyz'
164 #if (defined(HEINO))
165 ch_domain_long = trim(ch_domain_long)//'/ISMIP HEINO'
166 #endif
167 
168 #else
169 
170 errormsg = ' >>> sico_init: No valid domain specified!'
171 call error(errormsg)
172 
173 #endif
174 
175 !-------- Some initial values --------
176 
177 n_output = 0
178 
179 dtime = 0.0_dp
180 dtime_temp = 0.0_dp
181 dtime_wss = 0.0_dp
182 dtime_out = 0.0_dp
183 dtime_ser = 0.0_dp
184 
185 time = 0.0_dp
186 time_init = 0.0_dp
187 time_end = 0.0_dp
188 time_output = 0.0_dp
189 
190 !-------- Initialisation of the Library of Iterative Solvers Lis,
191 ! if required --------
192 
193 #if (CALCTHK==3 || CALCTHK==6 || MARGIN==3 || DYNAMICS==2)
194  call lis_initialize(ierr)
195 #endif
196 
197 !-------- Read physical parameters --------
198 
199 call read_phys_para()
200 
201 call ice_mat_eqs_pars(rf, r_t, kappa, c, -190, 10)
202 
203 ! ------ Some auxiliary quantities required for the enthalpy method
204 
205 call calc_c_int_table(c, -190, 10, l)
207 
208 !-------- Compatibility check of the SICOPOLIS version with the header file
209 
210 if ( trim(version) /= trim(sico_version) ) then
211  errormsg = ' >>> sico_init: ' &
212  //'SICOPOLIS version not compatible with header file!'
213  call error(errormsg)
214 end if
215 
216 !-------- Check whether the dynamics and thermodynamics modes are defined
217 
218 #if (!defined(DYNAMICS))
219 errormsg = ' >>> sico_init: DYNAMICS not defined in the header file!'
220 call error(errormsg)
221 #endif
222 
223 #if (!defined(CALCMOD))
224 errormsg = ' >>> sico_init: CALCMOD not defined in the header file!'
225 call error(errormsg)
226 #endif
227 
228 #if (defined(ENTHMOD))
229 errormsg = ' >>> sico_init: ENTHMOD must not be defined any more.' &
230  // end_of_line &
231  //' Please update your header file!'
232 call error(errormsg)
233 #endif
234 
235 !-------- Compatibility check of the horizontal resolution with the
236 ! number of grid points --------
237 
238 #if (GRID==0 || GRID==1)
239 
240  errormsg = ' >>> sico_init: ' &
241  //'GRID==0 or 1 not allowed for the Tibet application!'
242  call error(errormsg)
243 
244 #elif (GRID==2)
245 
246 
247 if ( (approx_equal(dlambda, 1.0_dp/3.0_dp, eps_sp_dp)) &
248  .and.(approx_equal(dphi , 1.0_dp/3.0_dp, eps_sp_dp)) ) then
249 
250  if ((imax /= 135).or.(jmax /= 51)) then
251  errormsg = ' >>> sico_init: IMAX and/or JMAX wrong!'
252  call error(errormsg)
253  end if
254 
255 else if ( (approx_equal(dlambda, 1.0_dp/6.0_dp, eps_sp_dp)) &
256  .and.(approx_equal(dphi , 1.0_dp/6.0_dp, eps_sp_dp)) ) then
257 
258  if ((imax /= 270).or.(jmax /= 102)) then
259  errormsg = ' >>> sico_init: IMAX and/or JMAX wrong!'
260  call error(errormsg)
261  end if
262 
263 else if ( (approx_equal(dlambda, 1.0_dp/12.0_dp, eps_sp_dp)) &
264  .and.(approx_equal(dphi , 1.0_dp/12.0_dp, eps_sp_dp)) ) then
265 
266  if ((imax /= 540).or.(jmax /= 204)) then
267  errormsg = ' >>> sico_init: IMAX and/or JMAX wrong!'
268  call error(errormsg)
269  end if
270 
271 else
272 
273  errormsg = ' >>> sico_init: DLAMBDA / DPHI wrong!'
274  call error(errormsg)
275 
276 end if
277 
278 #endif
279 
280 !-------- Compatibility check of the thermodynamics mode
281 ! (cold vs. polythermal vs. enthalpy method)
282 ! and the number of grid points in the lower (kt) ice domain --------
283 
284 #if (CALCMOD==0 || CALCMOD==2 || CALCMOD==3 || CALCMOD==-1)
285 
286 if (ktmax > 2) then
287  write(6, fmt='(a)') ' >>> sico_init: For options CALCMOD==0, 2, 3 or -1,'
288  write(6, fmt='(a)') ' the separate kt domain is redundant.'
289  write(6, fmt='(a)') ' Therefore, consider setting KTMAX to 2.'
290  write(6, fmt='(a)') ' '
291 end if
292 
293 #endif
294 
295 !-------- Compatibility check of surface-temperature and precipitation
296 ! determination by interpolation between present and LGM values
297 ! with a glacial index --------
298 
299 #if (TSURFACE == 5 && ACCSURFACE != 5)
300 errormsg = ' >>> sico_init: ' &
301  //'Options TSURFACE==5 and ACCSURFACE==5 must be used together!'
302 call error(errormsg)
303 #endif
304 
305 #if (TSURFACE != 5 && ACCSURFACE == 5)
306 errormsg = ' >>> sico_init: ' &
307  //'Options TSURFACE==5 and ACCSURFACE==5 must be used together!'
308 call error(errormsg)
309 #endif
310 
311 !-------- Compatibility check of discretization schemes for the horizontal and
312 ! vertical advection terms in the temperature and age equations --------
313 
314 #if (ADV_HOR==1)
315 errormsg = ' >>> sico_init: ' &
316  //'Option ADV_HOR==1 (central differences) not defined!'
317 call error(errormsg)
318 #endif
319 
320 !-------- Check whether for the shallow shelf
321 ! or shelfy stream approximation
322 ! the chosen grid is Cartesian coordinates
323 ! without distortion correction (GRID==0) --------
324 
325 #if ((MARGIN==3 && DYNAMICS==1) || DYNAMICS==2) /* requires SSA or SStA */
326 #if (GRID != 0)
327 errormsg = ' >>> sico_init: Distortion correction not implemented' &
328  // end_of_line &
329  //' for the shallow shelf approximation (SSA)' &
330  // end_of_line &
331  //' or the shelfy stream approximation (SStA)' &
332  // end_of_line &
333  //' -> GRID==0 required!'
334 call error(errormsg)
335 #endif
336 #endif
337 
338 !-------- Setting of forcing flag --------
339 
340 #if (TSURFACE <= 4)
341 
342 forcing_flag = 1 ! forcing by delta_ts
343 
344 #elif (TSURFACE == 5)
345 
346 forcing_flag = 2 ! forcing by glac_index
347 
348 #endif
349 
350 !-------- Initialization of numerical time steps --------
351 
352 dtime0 = dtime0
353 dtime_temp0 = dtime_temp0
354 #if (REBOUND==2)
355 dtime_wss0 = dtime_wss0
356 #endif
357 
358 !-------- Further initializations --------
359 
360 dzeta_c = 1.0_dp/real(kcmax,dp)
361 dzeta_t = 1.0_dp/real(ktmax,dp)
362 dzeta_r = 1.0_dp/real(krmax,dp)
363 
364 ndat2d = 1
365 ndat3d = 1
366 
367 !-------- General abbreviations --------
368 
369 ! ------ kc domain
370 
371 if (deform >= eps) then
372 
373  flag_aa_nonzero = .true. ! non-equidistant grid
374 
375  aa = deform
376  ea = exp(aa)
377 
378  kc=0
379  zeta_c(kc) = 0.0_dp
380  eaz_c(kc) = 1.0_dp
381  eaz_c_quotient(kc) = 0.0_dp
382 
383  do kc=1, kcmax-1
384  zeta_c(kc) = kc*dzeta_c
385  eaz_c(kc) = exp(aa*zeta_c(kc))
386  eaz_c_quotient(kc) = (eaz_c(kc)-1.0_dp)/(ea-1.0_dp)
387  end do
388 
389  kc=kcmax
390  zeta_c(kc) = 1.0_dp
391  eaz_c(kc) = exp(aa)
392  eaz_c_quotient(kc) = 1.0_dp
393 
394 else
395 
396  flag_aa_nonzero = .false. ! equidistant grid
397 
398  aa = 0.0_dp
399  ea = 1.0_dp
400 
401  kc=0
402  zeta_c(kc) = 0.0_dp
403  eaz_c(kc) = 1.0_dp
404  eaz_c_quotient(kc) = 0.0_dp
405 
406  do kc=1, kcmax-1
407  zeta_c(kc) = kc*dzeta_c
408  eaz_c(kc) = 1.0_dp
409  eaz_c_quotient(kc) = zeta_c(kc)
410  end do
411 
412  kc=kcmax
413  zeta_c(kc) = 1.0_dp
414  eaz_c(kc) = 1.0_dp
415  eaz_c_quotient(kc) = 1.0_dp
416 
417 end if
418 
419 ! ------ kt domain
420 
421 kt=0
422 zeta_t(kt) = 0.0_dp
423 
424 do kt=1, ktmax-1
425  zeta_t(kt) = kt*dzeta_t
426 end do
427 
428 kt=ktmax
429 zeta_t(kt) = 1.0_dp
430 
431 ! ------ kr domain
432 
433 kr=0
434 zeta_r(kr) = 0.0_dp
435 
436 do kr=1, krmax-1
437  zeta_r(kr) = kr*dzeta_r
438 end do
439 
440 kr=krmax
441 zeta_r(kr) = 1.0_dp
442 
443 !-------- Construction of a vector (with index n) from a 2-d array
444 ! (with indices i, j) by diagonal numbering --------
445 
446 n=1
447 
448 do m=0, imax+jmax
449  do i=m, 0, -1
450  j = m-i
451  if ((i <= imax).and.(j <= jmax)) then
452  ii(n) = i
453  jj(n) = j
454  nn(j,i) = n
455  n=n+1
456  end if
457  end do
458 end do
459 
460 !-------- Specification of current simulation --------
461 
462 runname = runname
463 anfdatname = anfdatname
464 
465 #if (defined(YEAR_ZERO))
467 #else
468 year_zero = 2000.0_dp ! default value 2000 CE
469 #endif
470 
471 time_init0 = time_init0
472 time_end0 = time_end0
473 dtime_ser0 = dtime_ser0
474 
475 #if (OUTPUT==1 || OUTPUT==3)
476 dtime_out0 = dtime_out0
477 #endif
478 
479 #if (OUTPUT==2 || OUTPUT==3)
480 n_output = n_output
481 time_output0( 1) = time_out0_01
482 time_output0( 2) = time_out0_02
483 time_output0( 3) = time_out0_03
484 time_output0( 4) = time_out0_04
485 time_output0( 5) = time_out0_05
486 time_output0( 6) = time_out0_06
487 time_output0( 7) = time_out0_07
488 time_output0( 8) = time_out0_08
489 time_output0( 9) = time_out0_09
490 time_output0(10) = time_out0_10
491 time_output0(11) = time_out0_11
492 time_output0(12) = time_out0_12
493 time_output0(13) = time_out0_13
494 time_output0(14) = time_out0_14
495 time_output0(15) = time_out0_15
496 time_output0(16) = time_out0_16
497 time_output0(17) = time_out0_17
498 time_output0(18) = time_out0_18
499 time_output0(19) = time_out0_19
500 time_output0(20) = time_out0_20
501 #endif
502 
503 !-------- Write log file --------
504 
505 shell_command = 'if [ ! -d'
506 shell_command = trim(shell_command)//' '//outpath
507 shell_command = trim(shell_command)//' '//'] ; then mkdir'
508 shell_command = trim(shell_command)//' '//outpath
509 shell_command = trim(shell_command)//' '//'; fi'
510 call system(trim(shell_command))
511  ! Check whether directory OUTPATH exists. If not, it is created.
512 
513 filename_with_path = trim(outpath)//'/'//trim(runname)//'.log'
514 
515 open(10, iostat=ios, file=trim(filename_with_path), status='new')
516 
517 if (ios /= 0) then
518  errormsg = ' >>> sico_init: Error when opening the log file!'
519  call error(errormsg)
520 end if
521 
522 write(10, fmt=trim(fmt1)) 'Computational domain:'
523 write(10, fmt=trim(fmt1)) trim(ch_domain_long)
524 write(10, fmt=trim(fmt1)) ' '
525 
526 write(10, fmt=trim(fmt2a)) 'GRID = ', grid
527 write(10, fmt=trim(fmt1)) ' '
528 
529 write(10, fmt=trim(fmt2)) 'imax =', imax
530 write(10, fmt=trim(fmt2)) 'jmax =', jmax
531 write(10, fmt=trim(fmt2)) 'kcmax =', kcmax
532 write(10, fmt=trim(fmt2)) 'ktmax =', ktmax
533 write(10, fmt=trim(fmt2)) 'krmax =', krmax
534 write(10, fmt=trim(fmt1)) ' '
535 
536 write(10, fmt=trim(fmt3)) 'a =', aa
537 write(10, fmt=trim(fmt1)) ' '
538 
539 #if (GRID==0 || GRID==1)
540 errormsg = ' >>> sico_init: GRID==0 or 1 not allowed for this application!'
541 call error(errormsg)
542 #elif (GRID==2)
543 write(10, fmt=trim(fmt3)) 'lambda0 =', lambda_0
544 write(10, fmt=trim(fmt3)) 'phi0 =', phi_0
545 write(10, fmt=trim(fmt3)) 'dlambda =', dlambda
546 write(10, fmt=trim(fmt3)) 'dphi =', dphi
547 #endif
548 write(10, fmt=trim(fmt1)) ' '
549 
550 write(10, fmt=trim(fmt3)) 'year_zero =', year_zero
551 write(10, fmt=trim(fmt3)) 'time_init =', time_init0
552 write(10, fmt=trim(fmt3)) 'time_end =', time_end0
553 write(10, fmt=trim(fmt3)) 'dtime =', dtime0
554 write(10, fmt=trim(fmt3)) 'dtime_temp =', dtime_temp0
555 #if (REBOUND==2)
556 write(10, fmt=trim(fmt3)) 'dtime_wss =', dtime_wss0
557 #endif
558 write(10, fmt=trim(fmt1)) ' '
559 
560 write(10, fmt=trim(fmt2)) 'ANF_DAT =', anf_dat
561 write(10, fmt=trim(fmt1)) 'zs_present file = '//zs_present_file
562 #if (ANF_DAT==1)
563 #if (defined(ZB_PRESENT_FILE))
564 write(10, fmt=trim(fmt1)) 'zb_present file = '//zb_present_file
565 #endif
566 write(10, fmt=trim(fmt1)) 'zl_present file = '//zl_present_file
567 #endif
568 write(10, fmt=trim(fmt1)) 'zl0 file = '//zl0_file
569 write(10, fmt=trim(fmt1)) 'mask_present file = '//mask_present_file
570 #if (ANF_DAT==1 && defined(TEMP_INIT))
571 write(10, fmt=trim(fmt2)) 'TEMP_INIT =', temp_init
572 #endif
573 #if (ANF_DAT==3)
574 write(10, fmt=trim(fmt1)) 'Initial-value file = '//anfdatname
575 write(10, fmt=trim(fmt1)) 'Path to initial-value file = '//anfdatpath
576 #endif
577 write(10, fmt=trim(fmt1)) ' '
578 
579 write(10, fmt=trim(fmt1)) 'Physical-parameter file = '//phys_para_file
580 write(10, fmt=trim(fmt1)) ' '
581 
582 #if (CALCTHK==2 || CALCTHK==3 || CALCTHK==5 || CALCTHK==6)
583 write(10, fmt=trim(fmt3)) 'ovi_weight =', ovi_weight
584 #if (CALCTHK==2 || CALCTHK==5)
585 write(10, fmt=trim(fmt3)) 'omega_sor =', omega_sor
586 #if (ITER_MAX_SOR>0)
587 write(10, fmt=trim(fmt2a)) 'iter_max_sor = ', iter_max_sor
588 #endif
589 #endif
590 write(10, fmt=trim(fmt1)) ' '
591 #endif
592 
593 write(10, fmt=trim(fmt1)) 'temp_mm_present file = '//temp_mm_present_file
594 #if (TSURFACE==1)
595 write(10, fmt=trim(fmt3)) 'delta_ts0 =', delta_ts0
596 #elif (TSURFACE==3)
597 write(10, fmt=trim(fmt3)) 'sine_amplit =', sine_amplit
598 write(10, fmt=trim(fmt3)) 'sine_period =', sine_period
599 #elif (TSURFACE==4)
600 write(10, fmt=trim(fmt1)) 'GRIP file = '//grip_temp_file
601 write(10, fmt=trim(fmt3)) 'grip_temp_fact =', grip_temp_fact
602 #elif (TSURFACE==5)
603 write(10, fmt=trim(fmt1)) 'Glacial-index file = '//glac_ind_file
604 write(10, fmt=trim(fmt1)) 'temp_mm_anom file = '//temp_mm_anom_file
605 write(10, fmt=trim(fmt3)) 'temp_mm_anom fact = ', temp_mm_anom_fact
606 #endif
607 
608 write(10, fmt=trim(fmt1)) 'precip_mm_present file = '//precip_mm_present_file
609 #if (ACCSURFACE==1)
610 write(10, fmt=trim(fmt3)) 'accfact =', accfact
611 #elif (ACCSURFACE==2 || ACCSURFACE==3)
612 write(10, fmt=trim(fmt3)) 'gamma_s =', gamma_s
613 #elif (ACCSURFACE==5)
614 write(10, fmt=trim(fmt1)) 'precip_mm_anom file = '//precip_mm_anom_file
615 write(10, fmt=trim(fmt3)) 'precip_mm_anom fact = ', precip_mm_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 write(10, fmt=trim(fmt2a)) 'SEA_LEVEL = ', sea_level
631 #if (SEA_LEVEL==1)
632 write(10, fmt=trim(fmt3)) 'z_sl0 =', z_sl0
633 #elif (SEA_LEVEL==3)
634 write(10, fmt=trim(fmt1)) 'sea-level file = '//sea_level_file
635 #endif
636 write(10, fmt=trim(fmt1)) ' '
637 
638 #if (MARGIN==2)
639 #if (MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3)
640 write(10, fmt=trim(fmt3)) 'z_mar =', z_mar
641 write(10, fmt=trim(fmt1)) ' '
642 #elif (MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5 \
643  || marine_ice_calving==6 || marine_ice_calving==7)
644 write(10, fmt=trim(fmt3)) 'fact_z_mar =', fact_z_mar
645 write(10, fmt=trim(fmt1)) ' '
646 #elif (MARINE_ICE_FORMATION==2 && MARINE_ICE_CALVING==9)
647 write(10, fmt=trim(fmt3)) 'calv_uw_coeff =', calv_uw_coeff
648 write(10, fmt=trim(fmt3)) 'r1_calv_uw =', r1_calv_uw
649 write(10, fmt=trim(fmt3)) 'r2_calv_uw =', r2_calv_uw
650 write(10, fmt=trim(fmt1)) ' '
651 #endif
652 #elif (MARGIN==3)
653 #if (ICE_SHELF_CALVING==2)
654 write(10, fmt=trim(fmt3)) 'H_calv =', h_calv
655 write(10, fmt=trim(fmt1)) ' '
656 #endif
657 #endif
658 
659 #if (defined(BASAL_HYDROLOGY))
660 write(10, fmt=trim(fmt2a)) 'BASAL_HYDROLOGY = ', basal_hydrology
661 #endif
662 write(10, fmt=trim(fmt2a)) 'SLIDE_LAW = ', slide_law
663 write(10, fmt=trim(fmt3)) 'c_slide =', c_slide
664 write(10, fmt=trim(fmt3)) 'gamma_slide =', gamma_slide
665 write(10, fmt=trim(fmt2a)) 'p_weert = ', p_weert
666 write(10, fmt=trim(fmt2a)) 'q_weert = ', q_weert
667 #if (defined(TIME_RAMP_UP_SLIDE))
668 write(10, fmt=trim(fmt3)) 'time_ramp_up_slide =', time_ramp_up_slide
669 #endif
670 #if (SLIDE_LAW==2 || SLIDE_LAW==3)
671 write(10, fmt=trim(fmt3)) 'red_pres_limit_fact =', red_pres_limit_fact
672 #endif
673 #if (BASAL_HYDROLOGY==1 \
674  && defined(hydro_slide_sat_fct) \
675  && defined(c_hw_slide) && defined(hw0_slide))
676 write(10, fmt=trim(fmt2a)) 'HYDRO_SLIDE_SAT_FCT = ', hydro_slide_sat_fct
677 write(10, fmt=trim(fmt3)) 'c_Hw_slide =', c_hw_slide
678 write(10, fmt=trim(fmt3)) 'Hw0_slide =', hw0_slide
679 #endif
680 write(10, fmt=trim(fmt1)) ' '
681 
682 write(10, fmt=trim(fmt2a)) 'Q_GEO_MOD = ', q_geo_mod
683 #if (Q_GEO_MOD==1)
684 write(10, fmt=trim(fmt3)) 'q_geo =', q_geo
685 #elif (Q_GEO_MOD==2)
686 write(10, fmt=trim(fmt1)) 'q_geo file = '//q_geo_file
687 #endif
688 write(10, fmt=trim(fmt1)) ' '
689 
690 #if (defined(MARINE_ICE_BASAL_MELTING))
691 write(10, fmt=trim(fmt2)) 'MARINE_ICE_BASAL_MELTING =', marine_ice_basal_melting
692 #if (MARINE_ICE_BASAL_MELTING==2 || MARINE_ICE_BASAL_MELTING==3)
693 write(10, fmt=trim(fmt3)) 'qbm_marine =', qbm_marine
694 #endif
695 write(10, fmt=trim(fmt1)) ' '
696 #endif
697 
698 write(10, fmt=trim(fmt2)) 'REBOUND =', rebound
699 #if (REBOUND==1)
700 write(10, fmt=trim(fmt3)) 'frac_llra =', frac_llra
701 #endif
702 #if (REBOUND==1 || REBOUND==2)
703 write(10, fmt=trim(fmt2)) 'TIME_LAG_MOD =', time_lag_mod
704 #if (TIME_LAG_MOD==1)
705 write(10, fmt=trim(fmt3)) 'time_lag =', time_lag
706 #elif (TIME_LAG_MOD==2)
707 write(10, fmt=trim(fmt1)) 'time_lag_file = '//time_lag_file
708 #else
709 errormsg = ' >>> sico_init: TIME_LAG_MOD must be either 1 or 2!'
710 call error(errormsg)
711 #endif
712 #endif
713 #if (REBOUND==2)
714 write(10, fmt=trim(fmt2)) 'FLEX_RIG_MOD =', flex_rig_mod
715 #if (FLEX_RIG_MOD==1)
716 write(10, fmt=trim(fmt3)) 'flex_rig =', flex_rig
717 #elif (FLEX_RIG_MOD==2)
718 write(10, fmt=trim(fmt1)) 'flex_rig_file = '//flex_rig_file
719 #else
720 errormsg = ' >>> sico_init: FLEX_RIG_MOD must be either 1 or 2!'
721 call error(errormsg)
722 #endif
723 #endif
724 write(10, fmt=trim(fmt2)) 'Q_LITHO =', q_litho
725 write(10, fmt=trim(fmt1)) ' '
726 
727 #if (FLOW_LAW==2)
728 write(10, fmt=trim(fmt3)) 'gr_size =', gr_size
729 write(10, fmt=trim(fmt1)) ' '
730 #endif
731 #if (FIN_VISC==2)
732 write(10, fmt=trim(fmt3)) 'sigma_res =', sigma_res
733 write(10, fmt=trim(fmt1)) ' '
734 #endif
735 
736 write(10, fmt=trim(fmt2)) 'ENHMOD =', enhmod
737 #if (ENHMOD==1 || ENHMOD==2 || ENHMOD==3)
738 write(10, fmt=trim(fmt3)) 'enh_fact =', enh_fact
739 #endif
740 #if (ENHMOD==2 || ENHMOD==3)
741 write(10, fmt=trim(fmt3)) 'enh_intg =', enh_intg
742 #endif
743 #if (ENHMOD==2)
744 write(10, fmt=trim(fmt3)) 'age_trans =', age_trans_0
745 #endif
746 #if (ENHMOD==3)
747 write(10, fmt=trim(fmt3)) 'date_trans1 =', date_trans1_0
748 write(10, fmt=trim(fmt3)) 'date_trans2 =', date_trans2_0
749 write(10, fmt=trim(fmt3)) 'date_trans3 =', date_trans3_0
750 #endif
751 #if (ENHMOD==4 || ENHMOD==5)
752 write(10, fmt=trim(fmt3)) 'enh_compr =', enh_compr
753 write(10, fmt=trim(fmt3)) 'enh_shear =', enh_shear
754 #endif
755 #if ((ENHMOD==1 || ENHMOD==2 || ENHMOD==3 || ENHMOD==4) && MARGIN==3)
756 write(10, fmt=trim(fmt3)) 'enh_shelf =', enh_shelf
757 #endif
758 write(10, fmt=trim(fmt1)) ' '
759 
760 write(10, fmt=trim(fmt3)) 'numdiff_H_t =', numdiff_h_t
761 write(10, fmt=trim(fmt3)) 'tau_cts =', tau_cts
762 write(10, fmt=trim(fmt3)) 'vh_max =', vh_max
763 write(10, fmt=trim(fmt3)) 'hd_min =', hd_min
764 write(10, fmt=trim(fmt3)) 'hd_max =', hd_max
765 #if (defined(QBM_MIN))
766 write(10, fmt=trim(fmt3)) 'qbm_min =', qbm_min
767 #elif (defined(QB_MIN)) /* obsolete */
768 write(10, fmt=trim(fmt3)) 'qb_min =', qb_min
769 #endif
770 #if (defined(QBM_MAX))
771 write(10, fmt=trim(fmt3)) 'qbm_max =', qbm_max
772 #elif (defined(QB_MAX)) /* obsolete */
773 write(10, fmt=trim(fmt3)) 'qb_max =', qb_max
774 #endif
775 write(10, fmt=trim(fmt3)) 'age_min =', age_min
776 write(10, fmt=trim(fmt3)) 'age_max =', age_max
777 write(10, fmt=trim(fmt3)) 'mean_accum =', mean_accum
778 #if (ADV_VERT==1)
779 write(10, fmt=trim(fmt3)) 'age_diff =', agediff
780 #endif
781 write(10, fmt=trim(fmt1)) ' '
782 
783 write(10, fmt=trim(fmt2a)) 'DYNAMICS = ', dynamics
784 #if ((DYNAMICS==1 && MARGIN==3) || DYNAMICS==2)
785 #if (defined(LIS_OPTS))
786 write(10, fmt=trim(fmt1)) 'lis_opts = '//lis_opts
787 #endif
788 #if (defined(N_ITER_SSA))
789 write(10, fmt=trim(fmt2a)) 'n_iter_ssa = ', n_iter_ssa
790 #endif
791 #if (defined(RELAX_FACT_SSA))
792 write(10, fmt=trim(fmt3)) 'relax_fact_ssa =', relax_fact_ssa
793 #endif
794 #endif
795 #if (DYNAMICS==2 && defined(RATIO_SL_THRESH))
796 write(10, fmt=trim(fmt3)) 'ratio_sl_thresh =', ratio_sl_thresh
797 #endif
798 #if (DYNAMICS==2 && defined(SSTA_SIA_WEIGH_FCT))
799 write(10, fmt=trim(fmt2a)) 'SSTA_SIA_WEIGH_FCT = ', ssta_sia_weigh_fct
800 #endif
801 write(10, fmt=trim(fmt1)) ' '
802 
803 write(10, fmt=trim(fmt2)) 'CALCMOD =', calcmod
804 #if (CALCMOD==-1 && defined(TEMP_CONST))
805 write(10, fmt=trim(fmt3)) 'TEMP_CONST =', temp_const
806 #endif
807 #if (CALCMOD==-1 && defined(AGE_CONST))
808 write(10, fmt=trim(fmt3)) 'AGE_CONST =', age_const
809 #endif
810 #if (CALCMOD==1 && defined(CTS_MELTING_FREEZING))
811 write(10, fmt=trim(fmt2)) 'CTS_MELTING_FREEZING =', cts_melting_freezing
812 #endif
813 write(10, fmt=trim(fmt2)) 'FLOW_LAW =', flow_law
814 write(10, fmt=trim(fmt2)) 'FIN_VISC =', fin_visc
815 write(10, fmt=trim(fmt2)) 'MARGIN =', margin
816 #if (MARGIN==2)
817 write(10, fmt=trim(fmt2)) 'MARINE_ICE_FORMATION =', marine_ice_formation
818 write(10, fmt=trim(fmt2)) 'MARINE_ICE_CALVING =', marine_ice_calving
819 #elif (MARGIN==3)
820 write(10, fmt=trim(fmt2)) 'ICE_SHELF_CALVING =', ice_shelf_calving
821 #endif
822 #if (defined(THK_EVOL))
823 write(10, fmt=trim(fmt2)) 'THK_EVOL =', thk_evol
824 #else
825 errormsg = ' >>> sico_init: Define THK_EVOL in header file!'
826 call error(errormsg)
827 #endif
828 #if (defined(CALCTHK))
829 write(10, fmt=trim(fmt2)) 'CALCTHK =', calcthk
830 #else
831 errormsg = ' >>> sico_init: Define CALCTHK in header file!'
832 call error(errormsg)
833 #endif
834 write(10, fmt=trim(fmt2)) 'ADV_HOR =', adv_hor
835 write(10, fmt=trim(fmt2)) 'ADV_VERT =', adv_vert
836 write(10, fmt=trim(fmt2)) 'TOPOGRAD =', topograd
837 write(10, fmt=trim(fmt2)) 'TSURFACE =', tsurface
838 write(10, fmt=trim(fmt2)) 'ACCSURFACE =', accsurface
839 #if (ACCSURFACE==5)
840 write(10, fmt=trim(fmt2)) 'PRECIP_ANOM_INTERPOL =', precip_anom_interpol
841 #endif
842 write(10, fmt=trim(fmt2)) 'SOLID_PRECIP =', solid_precip
843 write(10, fmt=trim(fmt2)) 'ABLSURFACE =', ablsurface
844 #if (defined(MB_ACCOUNT))
845 write(10, fmt=trim(fmt2a)) 'MB_ACCOUNT = ', mb_account
846 #endif
847 write(10, fmt=trim(fmt1)) ' '
848 
849 #if (defined(OUT_TIMES))
850 write(10, fmt=trim(fmt2a)) 'OUT_TIMES = ', out_times
851 #endif
852 #if (defined(OUTPUT_INIT))
853 write(10, fmt=trim(fmt2a)) 'OUTPUT_INIT = ', output_init
854 #endif
855 write(10, fmt=trim(fmt2a)) 'OUTPUT = ', output
856 #if (OUTPUT==1 || OUTPUT==3)
857 write(10, fmt=trim(fmt3)) 'dtime_out =' , dtime_out0
858 #endif
859 write(10, fmt=trim(fmt3)) 'dtime_ser =' , dtime_ser0
860 #if (OUTPUT==1 || OUTPUT==2)
861 write(10, fmt=trim(fmt2a)) 'ERGDAT = ', ergdat
862 #endif
863 write(10, fmt=trim(fmt2a)) 'NETCDF = ', netcdf
864 #if (OUTPUT==2 || OUTPUT==3)
865 write(10, fmt=trim(fmt2a)) 'n_output = ', n_output
866 do n=1, n_output
867  write(10, fmt=trim(fmt3)) 'time_output =' , time_output0(n)
868 end do
869 #endif
870 write(10, fmt=trim(fmt1)) ' '
871 
872 write(10, fmt=trim(fmt1)) 'Program version and date: '//version//' / '//date
873 
874 close(10, status='keep')
875 
876 !-------- Conversion of time quantities --------
877 
878 year_zero = year_zero*year_sec ! a --> s
879 time_init = time_init0*year_sec ! a --> s
880 time_end = time_end0*year_sec ! a --> s
881 dtime = dtime0*year_sec ! a --> s
882 dtime_temp = dtime_temp0*year_sec ! a --> s
883 #if (REBOUND==2)
884 dtime_wss = dtime_wss0*year_sec ! a --> s
885 #endif
886 dtime_ser = dtime_ser0*year_sec ! a --> s
887 #if (OUTPUT==1 || OUTPUT==3)
888 dtime_out = dtime_out0*year_sec ! a --> s
889 #endif
890 #if (OUTPUT==2 || OUTPUT==3)
891 do n=1, n_output
892  time_output(n) = time_output0(n)*year_sec ! a --> s
893 end do
894 #endif
895 
896 if (.not.approx_integer_multiple(dtime_temp, dtime, eps_sp_dp)) then
897  errormsg = ' >>> sico_init: dtime_temp must be a multiple of dtime!'
898  call error(errormsg)
899 end if
900 
901 #if (REBOUND==2)
902 if (.not.approx_integer_multiple(dtime_wss, dtime, eps_sp_dp)) then
903  errormsg = ' >>> sico_init: dtime_wss must be a multiple of dtime!'
904  call error(errormsg)
905 end if
906 #endif
907 
908 if (.not.approx_integer_multiple(dtime_ser, dtime, eps_sp_dp)) then
909  errormsg = ' >>> sico_init: dtime_ser must be a multiple of dtime!'
910  call error(errormsg)
911 end if
912 
913 #if (OUTPUT==1 || OUTPUT==3)
914 if (.not.approx_integer_multiple(dtime_out, dtime, eps_sp_dp)) then
915  errormsg = ' >>> sico_init: dtime_out must be a multiple of dtime!'
916  call error(errormsg)
917 end if
918 #endif
919 
920 time = time_init
921 
922 !-------- Reading of measurements for present monthly-mean precipitation --------
923 
924 #if (GRID==0 || GRID==1)
925 
926 errormsg = ' >>> sico_init: GRID==0 or 1 not allowed for the Tibet application!'
927 call error(errormsg)
928 
929 #elif (GRID==2)
930 
931 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
932  trim(precip_mm_present_file)
933 
934 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
935 
936 #endif
937 
938 if (ios /= 0) then
939  errormsg = ' >>> sico_init: Error when opening the precipitation file!'
940  call error(errormsg)
941 end if
942 
943 do m=1, 6; read(21, fmt='(a)') ch_dummy; end do
944 
945 do n=1, 12 ! month counter
946  do m=1, 3; read(21, fmt='(a)') ch_dummy; end do
947  do j=jmax, 0, -1
948  read(21, fmt=*) (precip_present(j,i,n), i=0,imax)
949  end do
950 end do
951 
952 close(21, status='keep')
953 
954 ! ------ Conversion mm/a water equivalent --> m/s ice equivalent
955 
956 precip_present = precip_present *(1.0e-03_dp/year_sec)*(rho_w/rho)
957  ! mm/a water equiv. -> m/s ice equiv.
958 
959 !-------- Reading of LGM monthly-mean precipitation-rate anomalies --------
960 
961 #if (ACCSURFACE==5)
962 
963 #if (GRID==2)
964 
965 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
966  trim(precip_mm_anom_file)
967 
968 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
969 
970 #endif
971 
972 if (ios /= 0) then
973  errormsg = ' >>> sico_init: Error when opening the precip anomaly file!'
974  call error(errormsg)
975 end if
976 
977 do m=1, 6; read(21, fmt='(a)') ch_dummy; end do
978 
979 do n=1, 12 ! month counter
980  do m=1, 3; read(21, fmt='(a)') ch_dummy; end do
981  do j=jmax, 0, -1
982  read(21, fmt=*) (precip_lgm_anom(j,i,n), i=0,imax)
983  end do
984 end do
985 
986 close(21, status='keep')
987 
988 precip_lgm_anom = precip_lgm_anom * precip_mm_anom_fact
989 
990 do i=0, imax
991 do j=0, jmax
992 
993 #if (PRECIP_ANOM_INTERPOL==1)
994  do n=1, 12 ! month counter
995  gamma_precip_lgm_anom(j,i,n) = 0.0_dp ! dummy values
996  end do
997 #elif (PRECIP_ANOM_INTERPOL==2)
998  do n=1, 12 ! month counter
999  gamma_precip_lgm_anom(j,i,n) = -log(precip_lgm_anom(j,i,n))
1000  end do
1001 #else
1002  errormsg = ' >>> sico_init: Wrong value of switch PRECIP_ANOM_INTERPOL!'
1003  call error(errormsg)
1004 #endif
1005 
1006 end do
1007 end do
1008 
1009 #endif
1010 
1011 !-------- Mean accumulation --------
1012 
1013 mean_accum = mean_accum*(1.0e-03_dp/year_sec)*(rho_w/rho)
1014  ! mm/a water equiv. -> m/s ice equiv.
1015 
1016 !-------- Reading of present topography mask --------
1017 
1018 #if (GRID==0 || GRID==1)
1019 
1020 errormsg = ' >>> sico_init: GRID==0 or 1 not allowed for the Tibet application!'
1021 call error(errormsg)
1022 
1023 #elif (GRID==2)
1024 
1025 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1026  trim(mask_present_file)
1027 
1028 open(24, iostat=ios, file=trim(filename_with_path), recl=rcl2, status='old')
1029 
1030 #endif
1031 
1032 if (ios /= 0) then
1033  errormsg = ' >>> sico_init: Error when opening the mask file!'
1034  call error(errormsg)
1035 end if
1036 
1037 do m=1, 6; read(24, fmt='(a)') ch_dummy; end do
1038 
1039 do j=jmax, 0, -1
1040  read(24, fmt=trim(fmt4)) (maske_ref(j,i), i=0,imax)
1041 end do
1042 
1043 close(24, status='keep')
1044 
1045 !-------- Reading of data for present monthly-mean surface temperature --------
1046 
1047 #if (GRID==2)
1048 
1049 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1050  trim(temp_mm_present_file)
1051 
1052 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1053 
1054 #endif
1055 
1056 if (ios /= 0) then
1057  errormsg = ' >>> sico_init: Error when opening the temperature file!'
1058  call error(errormsg)
1059 end if
1060 
1061 do m=1, 6; read(21, fmt='(a)') ch_dummy; end do
1062 
1063 do n=1, 12 ! month counter
1064  do m=1, 3; read(21, fmt='(a)') ch_dummy; end do
1065  do j=jmax, 0, -1
1066  read(21, fmt=*) (temp_mm_present(j,i,n), i=0,imax)
1067  end do
1068 end do
1069 
1070 close(21, status='keep')
1071 
1072 !-------- Reading of LGM monthly-mean surface-temperature anomalies --------
1073 
1074 #if (TSURFACE==5)
1075 
1076 #if (GRID==2)
1077 
1078 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1079  trim(temp_mm_anom_file)
1080 
1081 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1082 
1083 #endif
1084 
1085 if (ios /= 0) then
1086  errormsg = ' >>> sico_init: Error when opening the temperature anomaly file!'
1087  call error(errormsg)
1088 end if
1089 
1090 do m=1, 6; read(21, fmt='(a)') ch_dummy; end do
1091 
1092 do n=1, 12 ! month counter
1093  do m=1, 3; read(21, fmt='(a)') ch_dummy; end do
1094  do j=jmax, 0, -1
1095  read(21, fmt=*) (temp_mm_lgm_anom(j,i,n), i=0,imax)
1096  end do
1097 end do
1098 
1099 close(21, status='keep')
1100 
1101 temp_mm_lgm_anom = temp_mm_lgm_anom * temp_mm_anom_fact
1102 
1103 #endif
1104 
1105 !-------- Present reference elevation
1106 ! (for precipitation and surface-temperature data) --------
1107 
1108 #if (GRID==0 || GRID==1)
1109 
1110 errormsg = ' >>> sico_init: GRID==0 or 1 not allowed for the Tibet application!'
1111 call error(errormsg)
1112 
1113 #elif (GRID==2)
1114 
1115 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1116  trim(zs_present_file)
1117 
1118 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1119 
1120 #endif
1121 
1122 if (ios /= 0) then
1123  errormsg = ' >>> sico_init: Error when opening the zs file!'
1124  call error(errormsg)
1125 end if
1126 
1127 do m=1, 6; read(21, fmt='(a)') ch_dummy; end do
1128 
1129 do j=jmax, 0, -1
1130  read(21, fmt=*) (zs_ref(j,i), i=0,imax)
1131 end do
1132 
1133 close(21, status='keep')
1134 
1135 ! ------ Reset bathymetry data (sea floor elevation) to the
1136 ! sea surface
1137 
1138 do i=0, imax
1139 do j=0, jmax
1140  if (maske_ref(j,i) >= 2_i1b) zs_ref(j,i) = 0.0_dp
1141 end do
1142 end do
1143 
1144 !-------- Read data for delta_ts --------
1145 
1146 #if (TSURFACE==4)
1147 
1148 filename_with_path = trim(inpath)//'/general/'//trim(grip_temp_file)
1149 
1150 open(21, iostat=ios, file=trim(filename_with_path), status='old')
1151 
1152 if (ios /= 0) then
1153  errormsg = ' >>> sico_init: Error when opening the data file for delta_ts!'
1154  call error(errormsg)
1155 end if
1156 
1157 read(21, fmt=*) ch_dummy, grip_time_min, grip_time_stp, grip_time_max
1158 
1159 if (ch_dummy /= '#') then
1160  errormsg = ' >>> sico_init: grip_time_min, grip_time_stp, grip_time_max' &
1161  // end_of_line &
1162  //' not defined in data file for delta_ts!'
1163  call error(errormsg)
1164 end if
1165 
1167 
1168 allocate(griptemp(0:ndata_grip))
1169 
1170 do n=0, ndata_grip
1171  read(21, fmt=*) d_dummy, griptemp(n)
1172 end do
1173 
1174 close(21, status='keep')
1175 
1176 #endif
1177 
1178 !-------- Read data for the glacial index --------
1179 
1180 #if (TSURFACE==5)
1181 
1182 filename_with_path = trim(inpath)//'/general/'//trim(glac_ind_file)
1183 
1184 open(21, iostat=ios, file=trim(filename_with_path), status='old')
1185 
1186 if (ios /= 0) then
1187  errormsg = ' >>> sico_init: Error when opening the glacial-index file!'
1188  call error(errormsg)
1189 end if
1190 
1191 read(21, fmt=*) ch_dummy, gi_time_min, gi_time_stp, gi_time_max
1192 
1193 if (ch_dummy /= '#') then
1194  errormsg = ' >>> sico_init: gi_time_min, gi_time_stp, gi_time_max' &
1195  // end_of_line &
1196  //' not defined in glacial-index file!'
1197  call error(errormsg)
1198 end if
1199 
1201 
1202 allocate(glacial_index(0:ndata_gi))
1203 
1204 do n=0, ndata_gi
1205  read(21, fmt=*) d_dummy, glacial_index(n)
1206 end do
1207 
1208 close(21, status='keep')
1209 
1210 #endif
1211 
1212 !-------- Read data for z_sl --------
1213 
1214 #if (SEA_LEVEL==3)
1215 
1216 filename_with_path = trim(inpath)//'/general/'//trim(sea_level_file)
1217 
1218 open(21, iostat=ios, file=trim(filename_with_path), status='old')
1219 
1220 if (ios /= 0) then
1221  errormsg = ' >>> sico_init: Error when opening the data file for z_sl!'
1222  call error(errormsg)
1223 end if
1224 
1225 read(21, fmt=*) ch_dummy, specmap_time_min, specmap_time_stp, specmap_time_max
1226 
1227 if (ch_dummy /= '#') then
1228  errormsg = ' >>> sico_init:' &
1229  // end_of_line &
1230  //' specmap_time_min, specmap_time_stp, specmap_time_max' &
1231  // end_of_line &
1232  //' not defined in data file for z_sl!'
1233  call error(errormsg)
1234 end if
1235 
1237 
1238 allocate(specmap_zsl(0:ndata_specmap))
1239 
1240 do n=0, ndata_specmap
1241  read(21, fmt=*) d_dummy, specmap_zsl(n)
1242 end do
1243 
1244 close(21, status='keep')
1245 
1246 #endif
1247 
1248 !-------- Determination of the geothermal heat flux --------
1249 
1250 #if (Q_GEO_MOD==1)
1251 
1252 ! ------ Constant value
1253 
1254 do i=0, imax
1255 do j=0, jmax
1256  q_geo(j,i) = q_geo *1.0e-03_dp ! mW/m2 -> W/m2
1257 end do
1258 end do
1259 
1260 #elif (Q_GEO_MOD==2)
1261 
1262 ! ------ Read data from file
1263 
1264 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1265  trim(q_geo_file)
1266 
1267 open(21, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1268 
1269 if (ios /= 0) then
1270  errormsg = ' >>> sico_init: Error when opening the qgeo file!'
1271  call error(errormsg)
1272 end if
1273 
1274 do m=1, 6; read(21, fmt='(a)') ch_dummy; end do
1275 
1276 do j=jmax, 0, -1
1277  read(21, fmt=*) (q_geo(j,i), i=0,imax)
1278 end do
1279 
1280 close(21, status='keep')
1281 
1282 do i=0, imax
1283 do j=0, jmax
1284  q_geo(j,i) = q_geo(j,i) *1.0e-03_dp ! mW/m2 -> W/m2
1285 end do
1286 end do
1287 
1288 #endif
1289 
1290 !-------- Reading of tabulated kei function--------
1291 
1292 #if (REBOUND==0 || REBOUND==1)
1293 
1294 kei = 0.0_dp; n_data_kei = 0; kei_r_max = 0.0_dp; kei_r_incr = 0.0_dp
1295  ! dummy values
1296 #elif (REBOUND==2)
1297 
1298 call read_kei()
1299 
1300 #endif
1301 
1302 !-------- Determination of the time lag
1303 ! of the relaxing asthenosphere --------
1304 
1305 #if (REBOUND==1 || REBOUND==2)
1306 
1307 #if (TIME_LAG_MOD==1)
1308 
1309 time_lag_asth = time_lag*year_sec ! a -> s
1310 
1311 #elif (TIME_LAG_MOD==2)
1312 
1313 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1314  trim(time_lag_file)
1315 
1316 open(29, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1317 
1318 if (ios /= 0) then
1319  errormsg = ' >>> sico_init: Error when opening the time-lag file!'
1320  call error(errormsg)
1321 end if
1322 
1323 do n=1, 6; read(29, fmt='(a)') ch_dummy; end do
1324 
1325 do j=jmax, 0, -1
1326  read(29, fmt=*) (time_lag_asth(j,i), i=0,imax)
1327 end do
1328 
1329 close(29, status='keep')
1330 
1331 time_lag_asth = time_lag_asth*year_sec ! a -> s
1332 
1333 #endif
1334 
1335 #elif (REBOUND==0)
1336 
1337 time_lag_asth = 0.0_dp ! dummy values
1338 
1339 #endif
1340 
1341 !-------- Determination of the flexural rigidity of the lithosphere --------
1342 
1343 #if (REBOUND==2)
1344 
1345 #if (FLEX_RIG_MOD==1)
1346 
1347 flex_rig_lith = flex_rig
1348 
1349 #elif (FLEX_RIG_MOD==2)
1350 
1351 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1352  trim(flex_rig_file)
1353 
1354 open(29, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1355 
1356 if (ios /= 0) then
1357  errormsg = ' >>> sico_init: Error when opening the flex-rig file!'
1358  call error(errormsg)
1359 end if
1360 
1361 do n=1, 6; read(29, fmt='(a)') ch_dummy; end do
1362 
1363 do j=jmax, 0, -1
1364  read(29, fmt=*) (flex_rig_lith(j,i), i=0,imax)
1365 end do
1366 
1367 close(29, status='keep')
1368 
1369 #endif
1370 
1371 #elif (REBOUND==0 || REBOUND==1)
1372 
1373 flex_rig_lith = 0.0_dp ! dummy values
1374 
1375 #endif
1376 
1377 !-------- Definition of initial values --------
1378 
1379 ! ------ Present topography
1380 
1381 #if (ANF_DAT==1)
1382 
1383 errormsg = ' >>> sico_init: ANF_DAT==1 not allowed for Tibet application!'
1384 call error(errormsg)
1385 
1386 ! ------ Ice-free, relaxed bedrock
1387 
1388 #elif (ANF_DAT==2)
1389 
1390 call topography2(dxi, deta)
1391 
1392 z_sl = -1.11e+11_dp ! dummy value for initial call of subroutine boundary
1393 
1394 call boundary(time_init, dtime, dxi, deta, &
1395  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1396 
1397 as_perp_apl = 0.0_dp
1398 
1399 smb_corr = 0.0_dp
1400 
1401 q_bm = 0.0_dp
1402 q_tld = 0.0_dp
1403 q_b_tot = 0.0_dp
1404 
1405 p_b_w = 0.0_dp
1406 h_w = 0.0_dp
1407 
1408 call init_temp_water_age_2()
1409 
1410 #if (ENHMOD==1)
1411  call calc_enhance_1()
1412 #elif (ENHMOD==2)
1413  call calc_enhance_2()
1414 #elif (ENHMOD==3)
1415  call calc_enhance_3(time_init)
1416 #elif (ENHMOD==4)
1417  call calc_enhance_4()
1418 #elif (ENHMOD==5)
1419  call calc_enhance_5()
1420 #else
1421  errormsg = ' >>> sico_init: Parameter ENHMOD must be between 1 and 5!'
1422  call error(errormsg)
1423 #endif
1424 
1425 ! ------ Read initial state from output of previous simulation
1426 
1427 #elif (ANF_DAT==3)
1428 
1429 call topography3(dxi, deta, z_sl, anfdatname)
1430 
1431 call boundary(time_init, dtime, dxi, deta, &
1432  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1433 
1434 where ((maske==0_i1b).or.(maske==3_i1b))
1435  ! grounded or floating ice
1437 elsewhere ! maske==1_i1b or 2_i1b, ice-free land or sea
1438  as_perp_apl = 0.0_dp
1439 end where
1440 
1441 smb_corr = 0.0_dp
1442 
1443 q_b_tot = q_bm + q_tld
1444 
1445 #endif
1446 
1447 !-------- Inner-point flag --------
1448 
1449 flag_inner_point = .true.
1450 
1451 flag_inner_point(0,:) = .false.
1452 flag_inner_point(jmax,:) = .false.
1453 
1454 flag_inner_point(:,0) = .false.
1455 flag_inner_point(:,imax) = .false.
1456 
1457 !-------- Distance between grid points with delta_i=ir, delta_j=jr --------
1458 
1459 #if (GRID==2)
1460 
1461 do ir=-imax, imax
1462 do jr=-jmax, jmax
1463 
1464  dist_dxdy(jr,ir) = sqrt( (sq_g11_g(jmax/2,imax/2)*real(ir,dp)*dxi)**2 &
1465  + (sq_g22_g(jmax/2,imax/2)*real(jr,dp)*deta)**2 )
1466 
1467  ! This uses the metric tensor in the center of the Tibet
1468  ! domain for the entire domain; DIRTY TRICK !!!
1469 
1470 end do
1471 end do
1472 
1473 #endif
1474 
1475 !-------- Initial velocities --------
1476 
1477 call calc_temp_melt()
1478 call flag_update_gf_gl_cf()
1479 
1480 #if (DYNAMICS==1 || DYNAMICS==2)
1481 
1482 call calc_vxy_b_sia(time, z_sl)
1483 call calc_vxy_sia(dzeta_c, dzeta_t)
1484 
1485 #if (MARGIN==3 || DYNAMICS==2)
1486 call calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t, ii, jj, nn)
1487 #endif
1488 
1489 call calc_vz_grounded(dxi, deta, dzeta_c, dzeta_t)
1490 
1491 #if (MARGIN==3)
1492 call calc_vz_floating(z_sl, dxi, deta, dzeta_c)
1493 #endif
1494 
1495 #elif (DYNAMICS==0)
1496 
1497 call calc_vxy_static()
1498 call calc_vz_static()
1499 
1500 #else
1501 
1502 errormsg = ' >>> sico_init: DYNAMICS must be either 0, 1 or 2!'
1503 call error(errormsg)
1504 
1505 #endif
1506 
1507 call calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
1508 
1509 !-------- Initial enthalpies --------
1510 
1511 #if (CALCMOD==0 || CALCMOD==-1)
1512 
1513 do i=0, imax
1514 do j=0, jmax
1515 
1516  do kc=0, kcmax
1517  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
1518  end do
1519 
1520  do kt=0, ktmax
1521  enth_t(kt,j,i) = enth_c(0,j,i)
1522  end do
1523 
1524 end do
1525 end do
1526 
1527 #elif (CALCMOD==1)
1528 
1529 do i=0, imax
1530 do j=0, jmax
1531 
1532  do kc=0, kcmax
1533  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
1534  end do
1535 
1536  if ( (maske(j,i) == 0_i1b).and.(n_cts(j,i) == 1_i1b) ) then
1537  do kt=0, ktmax
1538  enth_t(kt,j,i) = enth_fct_temp_omega(temp_t_m(kt,j,i), omega_t(kt,j,i))
1539  end do
1540  else
1541  do kt=0, ktmax
1542  enth_t(kt,j,i) = enth_c(0,j,i)
1543  end do
1544  end if
1545 
1546 end do
1547 end do
1548 
1549 #elif (CALCMOD==2 || CALCMOD==3)
1550 
1551 do i=0, imax
1552 do j=0, jmax
1553 
1554  do kc=0, kcmax
1555  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), omega_c(kc,j,i))
1556  end do
1557 
1558  do kt=0, ktmax
1559  enth_t(kt,j,i) = enth_c(0,j,i)
1560  end do
1561 
1562 end do
1563 end do
1564 
1565 #else
1566 
1567 errormsg = ' >>> sico_init: Parameter CALCMOD must be either -1, 0, 1, 2 or 3!'
1568 call error(errormsg)
1569 
1570 #endif
1571 
1572 !-------- Initialize time-series files --------
1573 
1574 ! ------ Time-series file for the ice sheet on the whole
1575 
1576 filename_with_path = trim(outpath)//'/'//trim(runname)//'.ser'
1577 
1578 open(12, iostat=ios, file=trim(filename_with_path), status='new')
1579 
1580 if (ios /= 0) then
1581  errormsg = ' >>> sico_init: Error when opening the ser file!'
1582  call error(errormsg)
1583 end if
1584 
1585 if (forcing_flag == 1) then
1586 
1587  write(12,1102)
1588  write(12,1103)
1589 
1590  1102 format(' t(a) D_Ts(deg C) z_sl(m)',/, &
1591  ' V(m^3) V_g(m^3) V_f(m^3)', &
1592  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1593  ' V_sle(m) V_t(m^3)', &
1594  ' A_t(m^2)',/, &
1595  ' H_max(m) H_t_max(m)', &
1596  ' zs_max(m) vs_max(m/a) Tbh_max(C)')
1597  1103 format('----------------------------------------------------', &
1598  '---------------------------------------')
1599 
1600 else if (forcing_flag == 2) then
1601 
1602  write(12,1112)
1603  write(12,1113)
1604 
1605  1112 format(' t(a) glac_ind(1) z_sl(m)',/, &
1606  ' V(m^3) V_g(m^3) V_f(m^3)', &
1607  ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1608  ' V_sle(m) V_t(m^3)', &
1609  ' A_t(m^2)',/, &
1610  ' H_max(m) H_t_max(m)', &
1611  ' zs_max(m) vs_max(m/a) Tbh_max(C)')
1612  1113 format('----------------------------------------------------', &
1613  '---------------------------------------')
1614 
1615 end if
1616 
1617 ! ------ Time-series file for deep boreholes
1618 
1619 n_core = 0 ! No boreholes defined
1620 
1621 filename_with_path = trim(outpath)//'/'//trim(runname)//'.core'
1622 
1623 open(14, iostat=ios, file=trim(filename_with_path), status='new')
1624 
1625 write(14,'(1x,a)') '---------------------'
1626 write(14,'(1x,a)') 'No boreholes defined.'
1627 write(14,'(1x,a)') '---------------------'
1628 
1629 !-------- Output of the initial state --------
1630 
1631 #if (defined(OUTPUT_INIT))
1632 
1633 #if (OUTPUT_INIT==0)
1634  flag_init_output = .false.
1635 #elif (OUTPUT_INIT==1)
1636  flag_init_output = .true.
1637 #else
1638  errormsg = ' >>> sico_init: OUTPUT_INIT must be either 0 or 1!'
1639  call error(errormsg)
1640 #endif
1641 
1642 #else
1643  flag_init_output = .true. ! default for undefined parameter OUTPUT_INIT
1644 #endif
1645 
1646 #if (OUTPUT==1)
1647 
1648 #if (ERGDAT==0)
1649  flag_3d_output = .false.
1650 #elif (ERGDAT==1)
1651  flag_3d_output = .true.
1652 #else
1653  errormsg = ' >>> sico_init: ERGDAT must be either 0 or 1!'
1654  call error(errormsg)
1655 #endif
1656 
1657  if (flag_init_output) &
1658  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1659  flag_3d_output, ndat2d, ndat3d)
1660 
1661 #elif (OUTPUT==2)
1662 
1663 if (time_output(1) <= time_init+eps) then
1664 
1665 #if (ERGDAT==0)
1666  flag_3d_output = .false.
1667 #elif (ERGDAT==1)
1668  flag_3d_output = .true.
1669 #else
1670  errormsg = ' >>> sico_init: ERGDAT must be either 0 or 1!'
1671  call error(errormsg)
1672 #endif
1673 
1674  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1675  flag_3d_output, ndat2d, ndat3d)
1676 
1677 end if
1678 
1679 #elif (OUTPUT==3)
1680 
1681  flag_3d_output = .false.
1682 
1683  if (flag_init_output) &
1684  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1685  flag_3d_output, ndat2d, ndat3d)
1686 
1687 if (time_output(1) <= time_init+eps) then
1688 
1689  flag_3d_output = .true.
1690 
1691  call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1692  flag_3d_output, ndat2d, ndat3d)
1693 
1694 end if
1695 
1696 #else
1697 
1698  errormsg = ' >>> sico_init: OUTPUT must be either 1, 2 or 3!'
1699  call error(errormsg)
1700 
1701 #endif
1702 
1703 if (flag_init_output) then
1704  call output2(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1705  call output4(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1706 end if
1707 
1708 end subroutine sico_init
1709 
1710 !-------------------------------------------------------------------------------
1711 !> Definition of the initial surface and bedrock topography
1712 !! (including gradients) and of the horizontal grid spacings dxi, deta.
1713 !! For ice-free initial topography with relaxed lithosphere
1714 !! (not defined for this domain, thus routine stops execution of SICOPOLIS).
1715 !<------------------------------------------------------------------------------
1716 subroutine topography1(dxi, deta)
1717 
1718 implicit none
1719 
1720 real(dp), intent(out) :: dxi, deta
1721 
1722 ! integer(i4b) :: i, j
1723 ! integer(i4b) :: ios
1724 ! real(dp) :: xi0, eta0
1725 
1726 dxi=0.0_dp; deta=0.0_dp ! dummy values
1727 
1728 errormsg = ' >>> topography1: ANF_DAT==1 not defined for this domain!'
1729 call error(errormsg)
1730 
1731 end subroutine topography1
1732 
1733 !-------------------------------------------------------------------------------
1734 !> Definition of the initial surface and bedrock topography
1735 !! (including gradients) and of the horizontal grid spacings dxi, deta.
1736 !! For ice-free initial topography with relaxed lithosphere.
1737 !<------------------------------------------------------------------------------
1738 subroutine topography2(dxi, deta)
1739 
1740 #if (GRID==0 || GRID==1)
1742 #endif
1743 
1744  use metric_m
1745  use topograd_m
1746 
1747 implicit none
1748 
1749 real(dp), intent(out) :: dxi, deta
1750 
1751 integer(i4b) :: i, j, n
1752 integer(i4b) :: ios
1753 real(dp) :: xi0, eta0
1754 character :: ch_dummy
1755 
1756 character(len= 8) :: ch_imax
1757 character(len=128) :: fmt4
1758 character(len=256) :: filename_with_path
1759 
1760 write(ch_imax, fmt='(i8)') imax
1761 write(fmt4, fmt='(a)') '('//trim(adjustl(ch_imax))//'(i1),i1)'
1762 
1763 !-------- Read topography --------
1764 
1765 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1766  trim(zl0_file)
1767 
1768 open(23, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1769 
1770 if (ios /= 0) then
1771  errormsg = ' >>> topography2: Error when opening the zl0 file!'
1772  call error(errormsg)
1773 end if
1774 
1775 do n=1, 6; read(23, fmt='(a)') ch_dummy; end do
1776 
1777 do j=jmax, 0, -1
1778  read(23, fmt=*) (zl0(j,i), i=0,imax)
1779 end do
1780 
1781 close(23, status='keep')
1782 
1783 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1784  trim(mask_present_file)
1785 
1786 open(24, iostat=ios, file=trim(filename_with_path), recl=rcl2, status='old')
1787 
1788 if (ios /= 0) then
1789  errormsg = ' >>> topography2: Error when opening the mask file!'
1790  call error(errormsg)
1791 end if
1792 
1793 do n=1, 6; read(24, fmt='(a)') ch_dummy; end do
1794 
1795 do j=jmax, 0, -1
1796  read(24, fmt=trim(fmt4)) (maske(j,i), i=0,imax)
1797 end do
1798 
1799 close(24, status='keep')
1800 
1801 !-------- Further stuff --------
1802 
1803 dxi = dlambda *pi_180 ! deg -> rad
1804 deta = dphi *pi_180 ! deg -> rad
1805 
1806 xi0 = lambda_0 *pi_180 ! deg -> rad
1807 eta0 = phi_0 *pi_180 ! deg -> rad
1808 
1809 do i=0, imax
1810 do j=0, jmax
1811 
1812  if (maske(j,i) <= 1_i1b) then
1813  maske(j,i) = 1_i1b
1814  zs(j,i) = zl0(j,i)
1815  zb(j,i) = zl0(j,i)
1816  zl(j,i) = zl0(j,i)
1817  else ! (maske(j,i) >= 2_i1b)
1818  maske(j,i) = 2_i1b
1819 #if (MARGIN==1 || MARGIN==2)
1820  zs(j,i) = zl0(j,i)
1821  zb(j,i) = zl0(j,i)
1822 #elif (MARGIN==3)
1823  zs(j,i) = 0.0_dp ! present-day
1824  zb(j,i) = 0.0_dp ! sea level
1825 #endif
1826  zl(j,i) = zl0(j,i)
1827  end if
1828 
1829  xi(i) = xi0 + real(i,dp)*dxi
1830  eta(j) = eta0 + real(j,dp)*deta
1831 
1832  zm(j,i) = zb(j,i)
1833  n_cts(j,i) = -1_i1b
1834  kc_cts(j,i) = 0
1835 
1836  h_c(j,i) = 0.0_dp
1837  h_t(j,i) = 0.0_dp
1838 
1839  dzs_dtau(j,i) = 0.0_dp
1840  dzm_dtau(j,i) = 0.0_dp
1841  dzb_dtau(j,i) = 0.0_dp
1842  dzl_dtau(j,i) = 0.0_dp
1843  dh_c_dtau(j,i) = 0.0_dp
1844  dh_t_dtau(j,i) = 0.0_dp
1845 
1846 end do
1847 end do
1848 
1849 maske_old = maske
1850 
1851 !-------- Geographic coordinates, metric tensor,
1852 ! gradients of the topography --------
1853 
1854 do i=0, imax
1855 do j=0, jmax
1856 
1857 #if (GRID==0 || GRID==1) /* Stereographic projection */
1858 
1859  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
1860  lambda0, phi0, lambda(j,i), phi(j,i))
1861 
1862 #elif (GRID==2) /* Geographic coordinates */
1863 
1864  lambda(j,i) = xi(i)
1865  phi(j,i) = eta(j)
1866 
1867 #endif
1868 
1869 end do
1870 end do
1871 
1872 call metric()
1873 
1874 #if (TOPOGRAD==0)
1875 call topograd_1(dxi, deta, 1)
1876 #elif (TOPOGRAD==1)
1877 call topograd_2(dxi, deta, 1)
1878 #endif
1879 
1880 !-------- Corresponding area of grid points --------
1881 
1882 do i=0, imax
1883 do j=0, jmax
1884  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
1885 end do
1886 end do
1887 
1888 end subroutine topography2
1889 
1890 !-------------------------------------------------------------------------------
1891 !> Definition of the initial surface and bedrock topography
1892 !! (including gradients) and of the horizontal grid spacings dxi, deta.
1893 !! For initial topography from previous simulation.
1894 !<------------------------------------------------------------------------------
1895 subroutine topography3(dxi, deta, z_sl, anfdatname)
1896 
1897  use read_m, only : read_erg_nc
1898 
1899 #if (GRID==0 || GRID==1)
1901 #endif
1902 
1903  use metric_m
1904  use topograd_m
1905 
1906 implicit none
1907 
1908 character(len=100), intent(in) :: anfdatname
1909 
1910 real(dp), intent(out) :: dxi, deta, z_sl
1911 
1912 integer(i4b) :: i, j, n
1913 integer(i4b) :: ios
1914 character(len=256) :: filename_with_path
1915 character :: ch_dummy
1916 
1917 !-------- Read data from time-slice file of previous simulation --------
1918 
1919 call read_erg_nc(z_sl, anfdatname)
1920 
1921 !-------- Read topography of the relaxed bedrock --------
1922 
1923 filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1924  trim(zl0_file)
1925 
1926 open(23, iostat=ios, file=trim(filename_with_path), recl=rcl1, status='old')
1927 
1928 if (ios /= 0) then
1929  errormsg = ' >>> topography3: Error when opening the zl0 file!'
1930  call error(errormsg)
1931 end if
1932 
1933 do n=1, 6; read(23, fmt='(a)') ch_dummy; end do
1934 
1935 do j=jmax, 0, -1
1936  read(23, fmt=*) (zl0(j,i), i=0,imax)
1937 end do
1938 
1939 close(23, status='keep')
1940 
1941 !-------- Further stuff --------
1942 
1943 dxi = dlambda *pi_180 ! deg -> rad
1944 deta = dphi *pi_180 ! deg -> rad
1945 
1946 !-------- Geographic coordinates, metric tensor,
1947 ! gradients of the topography --------
1948 
1949 do i=0, imax
1950 do j=0, jmax
1951 
1952 #if (GRID==0 || GRID==1) /* Stereographic projection */
1953 
1954  call stereo_inv_ellipsoid(xi(i), eta(j), a, b, &
1955  lambda0, phi0, lambda(j,i), phi(j,i))
1956 
1957 #elif (GRID==2) /* Geographic coordinates */
1958 
1959  lambda(j,i) = xi(i)
1960  phi(j,i) = eta(j)
1961 
1962 #endif
1963 
1964 end do
1965 end do
1966 
1967 call metric()
1968 
1969 #if (TOPOGRAD==0)
1970 call topograd_1(dxi, deta, 1)
1971 #elif (TOPOGRAD==1)
1972 call topograd_2(dxi, deta, 1)
1973 #endif
1974 
1975 !-------- Corresponding area of grid points --------
1976 
1977 do i=0, imax
1978 do j=0, jmax
1979  area(j,i) = sq_g11_g(j,i)*sq_g22_g(j,i)*dxi*deta
1980 end do
1981 end do
1982 
1983 end subroutine topography3
1984 
1985 !-------------------------------------------------------------------------------
1986 
1987 end module sico_init_m
1988 !
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) temp_c
temp_c(kc,j,i): Temperature in the upper (kc) ice domain
real(dp), dimension(0:jmax, 0:imax) c_slide
c_slide(j,i): Basal sliding coefficient
integer(i1b), dimension(0:jmax, 0:imax) n_cts
n_cts(j,i): Mask for thermal conditions. -1: cold ice base, 0: temperate ice base with cold ice above...
Update of the flags for the land-terminating grounded front, marine-terminating grounded front...
real(dp), dimension(0:jmax, 0:imax) phi
phi(j,i): Geographic latitude of grid point (i,j)
subroutine, public read_erg_nc(z_sl, filename)
Reading of time-slice files in native binary or NetCDF format.
Definition: read_m.F90:55
subroutine, public read_kei()
Reading of the tabulated kei function.
Definition: read_m.F90:845
logical, dimension(0:jmax, 0:imax) flag_inner_point
flag_inner_point(j,i): Inner-point flag. .true.: inner point, .false.: margin point ...
integer(i4b) grip_time_stp
grip_time_stp: Time step of the data values for the surface temperature anomaly
integer(i4b) ndata_specmap
ndata_specmap: Number of data values for the sea level
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_...
real(dp), dimension(0:kcmax) eaz_c
eaz_c(kc): Abbreviation for exp(aa*zeta(kc))
real(dp), dimension(0:jmax, 0:imax) q_geo
q_geo(j,i): Geothermal heat flux
real(dp), dimension(0:jmax, 0:imax, 12) temp_mm_present
temp_mm_present(j,i,n): Present-day mean monthly surface temperature
subroutine, public flag_update_gf_gl_cf()
Main subroutine of flag_update_gf_gl_cf_m: Update of the flags for the land-terminating grounded fron...
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)
Main routine of sico_init_m: Initialisations for SICOPOLIS.
Definition: sico_init_m.F90:60
real(dp) rho_w
RHO_W: Density of pure water.
subroutine, public calc_enhance_2()
Computation of the flow enhancement factor. Case ENHMOD==2: two different values depending on age for...
real(dp) a
A: Semi-major axis of the planet.
subroutine, public calc_vz_floating(z_sl, dxi, deta, dzeta_c)
Computation of the vertical velocity vz for floating ice.
Definition: calc_vz_m.F90:259
subroutine error(error_message)
Main routine of error_m: Writing of error messages and stopping execution.
Definition: error_m.F90:47
subroutine, public output1(runname, time, delta_ts, glac_index, z_sl, flag_3d_output, ndat2d, ndat3d)
Writing of time-slice files in native binary or NetCDF format.
Definition: output_m.F90:60
integer(i4b) specmap_time_max
specmap_time_max: Maximum time of the data values for the sea level
Material properties of ice: Rate factor, heat conductivity, specific heat (heat capacity), creep function, viscosity.
real(dp), parameter eps
eps: Small number
subroutine, public read_phys_para()
Reading of physical parameters.
Definition: read_m.F90:903
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) enth_t
enth_t(kt,j,i): Enthalpy in the lower (kt) ice domain
Conversion from temperature (temp) and water content (omega) to enthalpy (enth) and vice versa...
subroutine, public init_temp_water_age_2()
Initial temperature, water content and age (case ANF_DAT==2: ice-free conditions with relaxed bedrock...
real(dp), dimension(0:jmax, 0:imax) dh_t_dtau
dH_t_dtau(j,i): Derivative of H_t by tau (time)
real(dp), dimension(0:jmax, 0:imax) lambda
lambda(j,i): Geographic longitude of grid point (i,j)
real(dp), dimension(0:jmax, 0:imax) dh_c_dtau
dH_c_dtau(j,i): Derivative of H_c by tau (time)
real(dp), dimension(0:jmax, 0:imax) p_b_w
p_b_w(j,i): Basal water pressure
integer(i4b) specmap_time_min
specmap_time_min: Minimum time of the data values for the sea level
integer(i1b), dimension(0:jmax, 0:imax) maske_old
maske_old(j,i): Old value of maske (at the previous time step)
real(dp), dimension(0:jmax, 0:imax) as_perp
as_perp(j,i): Accumulation-ablation function at the ice surface (SMB)
real(dp), dimension(0:jmax, 0:imax) area
area(j,i): Area of grid cell associated with grid point (i,j)
real(dp), dimension(:), allocatable glacial_index
glacial_index(n): Data values for the glacial index
Computation of all components of the strain-rate tensor, the full effective strain rate and the shear...
Definition: calc_dxyz_m.F90:37
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) temp_t_m
temp_t_m(kt,j,i): Melting temperature in the lower (kt) ice domain
subroutine, public stereo_inv_ellipsoid(x_val, y_val, A, B, lambda0, phi0, lambda_val, phi_val)
Inverse stereographic projection for an ellipsoidal planet.
subroutine topograd_1(dxi, deta, n_switch)
Calculation of topography gradients on the staggered grid and on the grid points (the latter by secon...
Definition: topograd_m.F90:57
real(dp), dimension(0:jmax, 0:imax) zs_ref
zs_ref(j,i): Reference elevation for precip_present, temp_ma_present and temp_mj_present ...
subroutine, public calc_enhance_5()
Computation of the flow enhancement factor. Case ENHMOD==5: minimal anisotropic enhancement factor fo...
real(dp) b
B: Semi-minor axis of the planet.
real(dp), dimension(0:jmax, 0:imax) zl
zl(j,i): Coordinate z of the lithosphere surface
subroutine, public calc_vxy_static()
Computation of the horizontal velocity vx, vy, the horizontal volume flux qx, qy etc. for static ice.
subroutine, public 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 (and optionally in NetCDF f...
Definition: output_m.F90:5214
integer(i4b) grip_time_max
grip_time_max: Maximum time of the data values for the surface temperature anomaly ...
real(dp), dimension(0:jmax, 0:imax) q_tld
Q_tld(j,i): Water drainage rate from the temperate layer.
subroutine, public 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 ...
Declarations of global variables for SICOPOLIS (for the ANT domain).
Definition: sico_vars_m.F90:35
integer(i4b), parameter rcl1
rcl1: Maximum length of record for input files (with factor 3 safety margin)
Initialisations for SICOPOLIS.
Definition: sico_init_m.F90:35
integer(i4b) ndata_gi
ndata_gi: Number of data values for the glacial index
Computation of the flow enhancement factor.
integer(i4b) ndata_grip
ndata_grip: Number of data values for the surface temperature anomaly
Computation of the horizontal velocity vx, vy.
Definition: calc_vxy_m.F90:36
integer(i4b), dimension(0:jmax, 0:imax) q_weert
q_weert(j,i): Weertman exponent for the basal pressure
subroutine, public calc_temp_melt()
Computation of the melting temperatures.
real(dp) ea
ea: Abbreviation for exp(aa)
real(dp), parameter pi_180
pi_180: pi divided by 180 (-> deg to rad)
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) enth_c
enth_c(kc,j,i): Enthalpy in the upper (kc) ice domain
logical flag_aa_nonzero
flag_aa_nonzero: Flag for the exponential stretch parameter aa. .true.: aa greater than zero (non-equ...
Initial temperature, water content and age.
subroutine, public calc_vz_grounded(dxi, deta, dzeta_c, dzeta_t)
Computation of the vertical velocity vz for grounded ice.
Definition: calc_vz_m.F90:52
real(dp) aa
aa: Exponential stretch parameter of the non-equidistant vertical grid in the upper (kc) ice domain ...
Reading of several input data.
Definition: read_m.F90:35
Declarations of kind types for SICOPOLIS.
real(dp), dimension(0:jmax, 0:imax) zs
zs(j,i): Coordinate z of the surface topography
integer(i4b), dimension(0:jmax, 0:imax) kc_cts
kc_cts(j,i): Position kc of the CTS (for COLD, ENTC, ENTM)
integer(i4b) gi_time_min
gi_time_min: Minimum time of the data values for the glacial index
real(dp), dimension(0:jmax, 0:imax) h_c
H_c(j,i): Thickness of ice in the upper (kc) domain (thickness of the cold-ice layer for POLY...
real(dp), dimension(0:jmax, 0:imax) smb_corr
smb_corr(j,i): Diagnosed SMB correction
real(dp) kei_r_max
kei_r_max: Maximum value of the argument r of the tabulated kei function
real(dp) year_zero
year_zero: SICOPOLIS year zero in astronomical year numbering [ = signed year CE (AD) ] ...
real(dp), dimension(:), allocatable specmap_zsl
specmap_zsl(n): Data values for the sea level
Computation of the melting and basal temperatures.
subroutine topograd_2(dxi, deta, n_switch)
Calculation of topography gradients on the staggered grid and on the grid points (the latter by fourt...
Definition: topograd_m.F90:203
real(dp), dimension(-190:10) kappa
KAPPA(n): Tabulated values for the heat conductivity of ice.
Computation of the vertical velocity vz.
Definition: calc_vz_m.F90:35
integer(i4b) gi_time_stp
gi_time_stp: Time step of the data values for the glacial index
Writing of error messages and stopping execution.
Definition: error_m.F90:35
real(dp), dimension(0:jmax, 0:imax) dzm_dtau
dzm_dtau(j,i): Derivative of zm by tau (time)
real(dp), dimension(0:krmax) zeta_r
zeta_r(kr): Sigma coordinate zeta_r of grid point kr
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) omega_c
omega_c(kc,j,i): Water content in the upper (kc) ice domain
character(len=64) ch_domain_long
ch_domain_long: Long name of the computational domain
Definition of the components g11 and g22 of the metric tensor of the applied numerical coordinates...
Definition: metric_m.F90:37
Computation of the forward or inverse stereographic projection, alternatively for a spherical or an e...
real(dp), dimension(0:jmax, 0:imax) zl0
zl0(j,i): zl for isostatically relaxed ice-free conditions
subroutine, public calc_vz_static()
Computation of the vertical velocity vz for static ice.
Definition: calc_vz_m.F90:374
integer(i4b) gi_time_max
gi_time_max: Maximum time of the data values for the glacial index
real(dp) phi0
PHI0: Standard parallel of the stereographic projection.
character, parameter end_of_line
end_of_line: End-of-line string
subroutine, public ice_mat_eqs_pars(RF_table, R_T_val, KAPPA_table, C_table, n_tmp_min, n_tmp_max, RHO_I_val, RHO_C_val, KAPPA_C_val, C_C_val)
Setting of required physical parameters.
real(dp), dimension(0:jmax, 0:imax, 12) gamma_precip_lgm_anom
gamma_precip_lgm_anom(j,i,n): negative natural logarithm of precip_lgm_anom(j,i,n) ...
Comparison of single- or double-precision floating-point numbers.
subroutine, public calc_enhance_1()
Computation of the flow enhancement factor. Case ENHMOD==1: constant for grounded ice...
real(dp), dimension(0:jmax, 0:imax) zm
zm(j,i): Coordinate z of the bottom of the upper (kc) ice domain = top of the lower (kt) ice domain (...
real(dp), dimension(-190:10) rf
RF(n): Tabulated values for the rate factor of cold ice.
integer(i1b), dimension(0:jmax, 0:imax) maske
maske(j,i): Ice-land-ocean mask. 0: grounded ice, 1: ice-free land, 2: ocean, 3: floating ice ...
real(dp), dimension(0:jmax, 0:imax) q_b_tot
Q_b_tot(j,i): Sum of Q_bm and Q_tld.
real(dp), dimension(0:jmax, 0:imax) q_bm
Q_bm(j,i): Basal melting rate.
real(dp) function, public enth_fct_temp_omega(temp_val, omega_val)
Enthalpy as a function of temperature and water content.
subroutine, public 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_m.F90:561
real(dp), dimension(-190:10) c
C(n): Tabulated values for the specific heat of ice.
real(dp), dimension(0:kcmax) eaz_c_quotient
eaz_c_quotient(kc): Abbreviation for (eaz_c(kc)-1.0)/(ea-1.0)
real(dp), dimension(-10000:10000) kei
kei(n): Tabulated values of the kei function (Kelvin function of zero order)
subroutine boundary(time, dtime, dxi, deta, delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
Main routine of boundary_m: Computation of the surface temperature (must be less than 0 deg C!) and o...
Definition: boundary_m.F90:61
integer(i4b) specmap_time_stp
specmap_time_stp: Time step of the data values for the sea level
integer(i4b), parameter rcl2
rcl2: Maximum length of record for input mask files (with factor 3 safety margin) ...
subroutine, public output2(time, dxi, deta, delta_ts, glac_index, z_sl)
Writing of time-series data on file in ASCII format (and optionally in NetCDF format).
Definition: output_m.F90:3869
integer(i1b) forcing_flag
forcing_flag: Flag for the forcing type. 1: forcing by a spatially constant surface temperature anoma...
real(dp), dimension(-jmax:jmax,-imax:imax) dist_dxdy
dist_dxdy(jr,ir): Distance between grid points with delta_i=ir, delta_j=jr
Calculation of topography gradients on the staggered grid and on the grid points (including length re...
Definition: topograd_m.F90:39
real(dp), dimension(0:jmax, 0:imax) sq_g11_g
sq_g11_g(j,i): Square root of the coefficient g11 of the metric tensor on grid point (i...
real(dp), dimension(0:jmax, 0:imax) as_perp_apl
as_perp_apl(j,i): Applied accumulation-ablation function (SMB)
subroutine, public metric()
Main routine of module metric_m: Definition of the components g11 and g22 of the metric tensor of the...
Definition: metric_m.F90:54
subroutine, public calc_enhance_3(time)
Computation of the flow enhancement factor. Case ENHMOD==3: two different values depending on time of...
integer(i4b) grip_time_min
grip_time_min: Minimum time of the data values for the surface temperature anomaly ...
real(dp), dimension(0:jmax, 0:imax) time_lag_asth
time_lag_asth(j,i): Time lag of the relaxing asthenosphere
real(dp), dimension(0:jmax) eta
eta(j): Coordinate eta (= y) of grid point j
real(dp), dimension(0:imax) xi
xi(i): Coordinate xi (= x) of grid point i
character(len=16) ch_domain_short
ch_domain_short: Short name of the computational domain
real(dp) r_t
R_T: Coefficient of the water-content dependence in the rate factor for temperate ice...
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) omega_t
omega_t(kt,j,i): Water content in the lower (kt) ice domain
subroutine, public calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
Main subroutine of calc_dxyz_m: Computation of all components of the strain-rate tensor, the full effective strain rate and the shear fraction.
Definition: calc_dxyz_m.F90:57
subroutine, public calc_enhance_4()
Computation of the flow enhancement factor. Case ENHMOD==4: minimal anisotropic enhancement factor fo...
real(dp) kei_r_incr
kei_r_incr: Increment of the argument r of the tabulated kei function
real(dp), dimension(0:ktmax) zeta_t
zeta_t(kt): Sigma coordinate zeta_t of grid point kt
Computation of the surface temperature (must be less than 0 deg C!) and of the accumulation-ablation ...
Definition: boundary_m.F90:37
real(dp) l
L: Latent heat of ice.
real(dp), parameter eps_sp_dp
eps_sp_dp: Small number to single-precision accuracy in double precision
integer(i4b) n_data_kei
n_data_kei: Number of tabulated values of the kei function
character(len=256) errormsg
errormsg: Error-message string
real(dp) rho
RHO: Density of ice.
subroutine topography3(dxi, deta, z_sl, anfdatname)
Definition of the initial surface and bedrock topography (including gradients) and of the horizontal ...
integer(i4b), dimension(0:jmax, 0:imax) p_weert
p_weert(j,i): Weertman exponent for the basal shear stress
integer(i4b) n_core
n_core: Number of positions to be considered in the time-series file for deep boreholes ...
integer, parameter dp
Double-precision reals.
subroutine topography2(dxi, deta)
Definition of the initial surface and bedrock topography (including gradients) and of the horizontal ...
real(dp), dimension(0:jmax, 0:imax) zb
zb(j,i): Coordinate z of the ice base
real(dp), dimension(0:jmax, 0:imax, 12) precip_present
precip_present(j,i,n): Present-day mean monthly precipitation rate at the ice surface ...
Writing of output data on files.
Definition: output_m.F90:36
real(dp), dimension(0:jmax, 0:imax) dzb_dtau
dzb_dtau(j,i): Derivative of zb by tau (time)
real(dp), dimension(0:kcmax) zeta_c
zeta_c(kc): Sigma coordinate zeta_c of grid point kc
real(dp) lambda0
LAMBDA0: Reference longitude (central meridian) of the stereographic projection.
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.
real(dp), dimension(0:jmax, 0:imax) h_t
H_t(j,i): Thickness of ice in the lower (kt) domain (thickness of the temperate layer for POLY...
real(dp), dimension(0:jmax, 0:imax, 12) precip_lgm_anom
precip_lgm_anom(j,i,n): LGM anomaly (ratio LGM/present) of the mean monthly precipitation rate at the...
subroutine topography1(dxi, deta)
Definition of the initial surface and bedrock topography (including gradients) and of the horizontal ...
real(dp), dimension(0:jmax, 0:imax, 12) temp_mm_lgm_anom
temp_mm_lgm_anom(j,i,n): LGM anomaly (difference LGM - present) of the mean monthly surface temperatu...
real(dp), dimension(0:jmax, 0:imax) h_w
H_w(j,i): Thickness of the water column under the ice base.
subroutine, public 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_m.F90:55
Declarations of global variables for SICOPOLIS.
real(dp), dimension(0:jmax, 0:imax) dzs_dtau
dzs_dtau(j,i): Derivative of zs by tau (time)
real(dp), dimension(0:jmax, 0:imax) flex_rig_lith
flex_rig_lith(j,i): Flexural rigidity of the lithosphere
real(dp), dimension(:), allocatable griptemp
griptemp(n): Data values for the surface temperature anomaly
real(dp), dimension(0:jmax, 0:imax) dzl_dtau
dzl_dtau(j,i): Derivative of zl by tau (time)
real(dp), dimension(0:jmax, 0:imax) sq_g22_g
sq_g22_g(j,i): Square root of the coefficient g22 of the metric tensor on grid point (i...