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