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