37 dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, &
38 time, time_init, time_end, time_output, &
39 dxi, deta, dzeta_c, dzeta_t, dzeta_r, &
40 z_sl, dzsl_dtau, z_mar, &
42 ndat2d, ndat3d, n_output, &
54 integer(i4b),
intent(out) :: ii((imax+1)*(jmax+1)), &
55 jj((imax+1)*(jmax+1)), &
57 integer(i4b),
intent(out) :: ndat2d, ndat3d
58 integer(i4b),
intent(out) :: n_output
59 real(dp),
intent(out) :: delta_ts, glac_index
60 real(dp),
intent(out) :: mean_accum
61 real(dp),
intent(out) :: dtime, dtime_temp, dtime_wss, &
63 real(dp),
intent(out) :: time, time_init, time_end, time_output(100)
64 real(dp),
intent(out) :: dxi, deta, dzeta_c, dzeta_t, dzeta_r
65 real(dp),
intent(out) :: z_sl, dzsl_dtau, z_mar
66 character(len=100),
intent(out) :: runname
68 integer(i4b) :: i, j, kc, kt, kr, m, n, ir, jr
69 integer(i4b) :: ios, ios1, ios2, ios3, ios4
71 integer(i4b) :: ndata_insol
72 real(dp) :: dtime0, dtime_temp0, dtime_wss0, dtime_out0, dtime_ser0
73 real(dp) :: time_init0, time_end0, time_output0(100)
74 real(dp) :: time_chasm_init0, time_chasm_end0
76 character (len=100) :: anfdatname
77 character (len=256) :: shell_command
79 logical :: flag_3d_output
81 character(len=64),
parameter :: fmt1 =
'(a)', &
85 character(len= 8) :: ch_imax
86 character(len=128) :: fmt4
88 write(ch_imax, fmt=
'(i8)') imax
89 write(fmt4, fmt=
'(a)')
'('//trim(adjustl(ch_imax))//
'(i1),i1)'
94 ch_domain_long =
'Antarctica'
95 ch_domain_short =
'ant'
98 ch_domain_long =
'Austfonna'
99 ch_domain_short =
'asf'
101 #elif defined(EMTP2SGE)
102 ch_domain_long =
'EISMINT Phase 2 Simplified Geometry Experiment'
103 ch_domain_short =
'emtp2sge'
106 ch_domain_long =
'Greenland'
107 ch_domain_short =
'grl'
110 ch_domain_long =
'Northern hemisphere'
111 ch_domain_short =
'nhem'
114 ch_domain_long =
'Scandinavia and Eurasia'
115 ch_domain_short =
'scand'
118 ch_domain_long =
'Tibet'
119 ch_domain_short =
'tibet'
122 ch_domain_long =
'North polar cap of Mars'
123 ch_domain_short =
'nmars'
126 ch_domain_long =
'South polar cap of Mars'
127 ch_domain_short =
'smars'
130 ch_domain_long =
'XYZ'
131 ch_domain_short =
'xyz'
133 ch_domain_long = trim(ch_domain_long)//
'/ISMIP HEINO'
137 stop
' sico_init: No valid domain specified!'
158 #if (CALCZS==4 || MARGIN==3 || DYNAMICS==2)
161 call lis_initialize(ierr)
169 call
ice_mat_eqs_pars(rf, r_t, kappa, c, -190, 10, rho_i, rho_c, kappa_c, c_c)
173 rho = (1.0_dp-frac_dust)*rho_i + frac_dust*rho_c
184 if ( trim(version) /= trim(sico_version) ) &
185 stop
' sico_init: SICOPOLIS version not compatible with header file!'
189 #if ( !defined(DYNAMICS) )
190 stop
' sico_init: DYNAMICS not defined in the header file!'
193 #if ( !defined(CALCMOD) )
194 stop
' sico_init: CALCMOD not defined in the header file!'
197 #if ( defined(ENTHMOD) )
198 write(6, fmt=
'(a)')
' sico_init: ENTHMOD must not be defined any more.'
199 write(6, fmt=
'(a)')
' Please update your header file!'
206 #if (GRID==0 || GRID==1)
208 if (abs(dx-20.0_dp) < eps)
then
210 if ((imax /= 90).or.(jmax /= 90))
then
211 stop
' sico_init: Wrong values for IMAX and JMAX!'
214 else if (abs(dx-10.0_dp) < eps)
then
216 if ((imax /= 180).or.(jmax /= 180))
then
217 stop
' sico_init: Wrong values for IMAX and JMAX!'
220 else if (abs(dx-5.0_dp) < eps)
then
222 if ((imax /= 360).or.(jmax /= 360))
then
223 stop
' sico_init: Wrong values for IMAX and JMAX!'
227 stop
' sico_init: Wrong value for DX!'
232 stop
' sico_init: GRID==2 not allowed for nmars application!'
240 #if (CALCMOD==0 || CALCMOD==2 || CALCMOD==3 || CALCMOD==-1)
243 write(6, fmt=
'(a)')
' sico_init: For options CALCMOD==0, 2, 3 or -1,'
244 write(6, fmt=
'(a)')
' the separate kt domain is redundant.'
245 write(6, fmt=
'(a)')
' Therefore, consider setting KTMAX to 1.'
246 write(6, fmt=
'(a)')
' '
255 stop
' sico_init: Option ADV_HOR==1 (central differences) not defined!'
263 #if ((MARGIN==3 && DYNAMICS==1) || DYNAMICS==2) /* requires SSA or SStA */
265 write(6, fmt=
'(a)')
' sico_init: Distortion correction not implemented'
266 write(6, fmt=
'(a)')
' for the shallow shelf approximation (SSA)'
267 write(6, fmt=
'(a)')
' or the shelfy stream approximation (SStA)'
268 write(6, fmt=
'(a)')
' -> GRID==0 required!'
280 dtime_temp0 = dtime_temp0
282 dtime_wss0 = dtime_wss0
287 dzeta_c = 1.0_dp/
real(kcmax,dp)
288 dzeta_t = 1.0_dp/
real(ktmax,dp)
289 dzeta_r = 1.0_dp/
real(krmax,dp)
298 if (deform >= eps)
then
300 flag_aa_nonzero = .true.
308 eaz_c_quotient(kc) = 0.0_dp
311 zeta_c(kc) = kc*dzeta_c
312 eaz_c(kc) = exp(aa*zeta_c(kc))
313 eaz_c_quotient(kc) = (eaz_c(kc)-1.0_dp)/(ea-1.0_dp)
319 eaz_c_quotient(kc) = 1.0_dp
323 flag_aa_nonzero = .false.
331 eaz_c_quotient(kc) = 0.0_dp
334 zeta_c(kc) = kc*dzeta_c
336 eaz_c_quotient(kc) = zeta_c(kc)
342 eaz_c_quotient(kc) = 1.0_dp
352 zeta_t(kt) = kt*dzeta_t
364 zeta_r(kr) = kr*dzeta_r
378 if ((i <= imax).and.(j <= jmax))
then
390 anfdatname = anfdatname
392 #if defined(YEAR_ZERO)
393 year_zero = year_zero
395 year_zero = 2000.0_dp
398 time_init0 = time_init0
399 time_end0 = time_end0
400 dtime_ser0 = dtime_ser0
403 time_chasm_init0 = time_chasm_init0
404 time_chasm_end0 = time_chasm_end0
407 #if ( OUTPUT==1 || OUTPUT==3 )
408 dtime_out0 = dtime_out0
411 #if ( OUTPUT==2 || OUTPUT==3 )
413 time_output0( 1) = time_out0_01
414 time_output0( 2) = time_out0_02
415 time_output0( 3) = time_out0_03
416 time_output0( 4) = time_out0_04
417 time_output0( 5) = time_out0_05
418 time_output0( 6) = time_out0_06
419 time_output0( 7) = time_out0_07
420 time_output0( 8) = time_out0_08
421 time_output0( 9) = time_out0_09
422 time_output0(10) = time_out0_10
423 time_output0(11) = time_out0_11
424 time_output0(12) = time_out0_12
425 time_output0(13) = time_out0_13
426 time_output0(14) = time_out0_14
427 time_output0(15) = time_out0_15
428 time_output0(16) = time_out0_16
429 time_output0(17) = time_out0_17
430 time_output0(18) = time_out0_18
431 time_output0(19) = time_out0_19
432 time_output0(20) = time_out0_20
437 shell_command =
'if [ ! -d'
438 shell_command = trim(shell_command)//
' '//outpath
439 shell_command = trim(shell_command)//
' '//
'] ; then mkdir'
440 shell_command = trim(shell_command)//
' '//outpath
441 shell_command = trim(shell_command)//
' '//
'; fi'
442 call system(trim(shell_command))
445 open(10, iostat=ios, &
446 file=outpath//
'/'//trim(runname)//
'.log', &
449 stop
' sico_init: Error when opening the log file!'
451 write(10, fmt=trim(fmt1))
'Computational domain:'
452 write(10, fmt=trim(fmt1)) trim(ch_domain_long)
453 write(10, fmt=trim(fmt1))
' '
455 write(10, fmt=trim(fmt2))
'imax =', imax
456 write(10, fmt=trim(fmt2))
'jmax =', jmax
457 write(10, fmt=trim(fmt2))
'kcmax =', kcmax
458 write(10, fmt=trim(fmt2))
'ktmax =', ktmax
459 write(10, fmt=trim(fmt2))
'krmax =', krmax
460 write(10, fmt=trim(fmt1))
' '
462 write(10, fmt=trim(fmt3))
'a =', aa
463 write(10, fmt=trim(fmt1))
' '
465 write(10, fmt=trim(fmt3))
'x0 =', x0
466 write(10, fmt=trim(fmt3))
'y0 =', y0
467 write(10, fmt=trim(fmt3))
'dx =', dx
468 write(10, fmt=trim(fmt1))
' '
470 write(10, fmt=trim(fmt3))
'year_zero =', year_zero
471 write(10, fmt=trim(fmt3))
'time_init =', time_init0
472 write(10, fmt=trim(fmt3))
'time_end =', time_end0
473 write(10, fmt=trim(fmt3))
'dtime =', dtime0
474 write(10, fmt=trim(fmt3))
'dtime_temp =', dtime_temp0
476 write(10, fmt=trim(fmt3))
'dtime_wss =', dtime_wss0
478 #if ( OUTPUT==1 || OUTPUT==3 )
479 write(10, fmt=trim(fmt3))
'dtime_out =', dtime_out0
481 write(10, fmt=trim(fmt3))
'dtime_ser =', dtime_ser0
482 write(10, fmt=trim(fmt1))
' '
484 write(10, fmt=trim(fmt2))
'ANF_DAT =', anf_dat
485 write(10, fmt=trim(fmt1))
'zs_present file = '//zs_present_file
487 write(10, fmt=trim(fmt1))
'zl_present file = '//zl_present_file
489 write(10, fmt=trim(fmt1))
'zl0 file = '//zl0_file
490 write(10, fmt=trim(fmt1))
'mask_present file = '//mask_present_file
491 #if (ANF_DAT==1 && defined(TEMP_INIT))
492 write(10, fmt=trim(fmt2))
'TEMP_INIT =', temp_init
495 write(10, fmt=trim(fmt1))
'Initial-value file = '//anfdatname
496 write(10, fmt=trim(fmt1))
'Path to initial-value file = '//anfdatpath
498 write(10, fmt=trim(fmt1))
' '
500 write(10, fmt=trim(fmt1))
'Physical-parameter file = '//phys_para_file
501 write(10, fmt=trim(fmt1))
' '
503 #if (CALCZS==3 || CALCZS==4)
504 write(10, fmt=trim(fmt3))
'ovi_weight =', ovi_weight
506 write(10, fmt=trim(fmt3))
'omega_sor =', omega_sor
508 write(10, fmt=trim(fmt1))
' '
511 #if (TSURFACE==1 || TSURFACE==6)
512 write(10, fmt=trim(fmt3))
'delta_ts0 =', delta_ts0
514 write(10, fmt=trim(fmt3))
'sine_amplit =', sine_amplit
515 write(10, fmt=trim(fmt3))
'sine_period =', sine_period
517 #if (TSURFACE==1 || TSURFACE==2 || TSURFACE==3)
518 write(10, fmt=trim(fmt3))
'temp0_ma_90N =', temp0_ma_90n
520 #if (TSURFACE==1 || TSURFACE==2 || TSURFACE==3 || TSURFACE==4 || TSURFACE==5)
521 write(10, fmt=trim(fmt3))
'c_ma =', c_ma
522 write(10, fmt=trim(fmt3))
'gamma_ma =', gamma_ma
525 write(10, fmt=trim(fmt1))
'Insolation file = '//insol_ma_90n_file
527 #if (TSURFACE==4 || TSURFACE==5 || TSURFACE==6)
528 write(10, fmt=trim(fmt3))
'albedo =', albedo
532 write(10, fmt=trim(fmt3))
'acc_present_val =', acc_present_val
535 #if ( ACCSURFACE==1 || ACCSURFACE==2 )
536 write(10, fmt=trim(fmt3))
'gamma_s =', gamma_s
539 #if ( ABLSURFACE==1 || ABLSURFACE==2 )
540 write(10, fmt=trim(fmt3))
'eld_0 =', eld_0
541 write(10, fmt=trim(fmt3))
'g_mb_0 =', g_0
542 write(10, fmt=trim(fmt3))
'gamma_eld =', gamma_eld
543 write(10, fmt=trim(fmt3))
'gamma_g =', gamma_g
547 write(10, fmt=trim(fmt3))
'z_sl0 =', z_sl0
549 write(10, fmt=trim(fmt1))
'sea-level file = '//sea_level_file
551 write(10, fmt=trim(fmt1))
' '
554 #if ( MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3 )
555 write(10, fmt=trim(fmt3))
'z_mar =', z_mar
556 write(10, fmt=trim(fmt1))
' '
557 #elif ( MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5 \
558 || marine_ice_calving==6 || marine_ice_calving==7 )
559 write(10, fmt=trim(fmt3))
'fact_z_mar =', fact_z_mar
560 write(10, fmt=trim(fmt1))
' '
561 #elif ( ( MARINE_ICE_FORMATION==2 ) && ( MARINE_ICE_CALVING==9 ) )
562 write(10, fmt=trim(fmt3))
'calv_uw_coeff =', calv_uw_coeff
563 write(10, fmt=trim(fmt3))
'r1_calv_uw =', r1_calv_uw
564 write(10, fmt=trim(fmt3))
'r2_calv_uw =', r2_calv_uw
565 write(10, fmt=trim(fmt1))
' '
568 #if ICE_SHELF_CALVING==2
569 write(10, fmt=trim(fmt3))
'H_calv =', h_calv
570 write(10, fmt=trim(fmt1))
' '
574 write(10, fmt=trim(fmt3))
'c_slide =', c_slide
575 write(10, fmt=trim(fmt3))
'gamma_slide =', gamma_slide
576 write(10, fmt=trim(fmt2))
'p_weert =', p_weert
577 write(10, fmt=trim(fmt2))
'q_weert =', q_weert
579 write(10, fmt=trim(fmt3))
'red_pres_limit_fact =', red_pres_limit_fact
581 write(10, fmt=trim(fmt1))
' '
584 write(10, fmt=trim(fmt3))
'q_geo =', q_geo
586 write(10, fmt=trim(fmt1))
'q_geo file = '//q_geo_file
588 write(10, fmt=trim(fmt1))
' '
590 write(10, fmt=trim(fmt2))
'REBOUND =', rebound
592 write(10, fmt=trim(fmt3))
'frac_llra =', frac_llra
594 #if (REBOUND==1 || REBOUND==2)
595 write(10, fmt=trim(fmt2))
'TIME_LAG_MOD =', time_lag_mod
596 #if (TIME_LAG_MOD==1)
597 write(10, fmt=trim(fmt3))
'time_lag =', time_lag
598 #elif (TIME_LAG_MOD==2)
599 write(10, fmt=trim(fmt1))
'time_lag_file = '//time_lag_file
601 stop
' sico_init: TIME_LAG_MOD must be either 1 or 2!'
605 write(10, fmt=trim(fmt2))
'FLEX_RIG_MOD =', flex_rig_mod
606 #if (FLEX_RIG_MOD==1)
607 write(10, fmt=trim(fmt3))
'flex_rig =', flex_rig
608 #elif (FLEX_RIG_MOD==2)
609 write(10, fmt=trim(fmt1))
'flex_rig_file = '//flex_rig_file
611 stop
' sico_init: FLEX_RIG_MOD must be either 1 or 2!'
614 write(10, fmt=trim(fmt2))
'Q_LITHO =', q_litho
615 write(10, fmt=trim(fmt1))
' '
618 write(10, fmt=trim(fmt3))
'gr_size =', gr_size
619 write(10, fmt=trim(fmt1))
' '
622 write(10, fmt=trim(fmt3))
'sigma_res =', sigma_res
623 write(10, fmt=trim(fmt1))
' '
625 write(10, fmt=trim(fmt3))
'frac_dust =', frac_dust
626 write(10, fmt=trim(fmt1))
' '
628 write(10, fmt=trim(fmt2))
'ENHMOD =', enhmod
629 #if (ENHMOD==1 || ENHMOD==2 || ENHMOD==3)
630 write(10, fmt=trim(fmt3))
'enh_fact =', enh_fact
632 #if (ENHMOD==2 || ENHMOD==3)
633 write(10, fmt=trim(fmt3))
'enh_intg =', enh_intg
636 write(10, fmt=trim(fmt3))
'age_trans =', age_trans_0
639 write(10, fmt=trim(fmt3))
'date_trans1 =', date_trans1_0
640 write(10, fmt=trim(fmt3))
'date_trans2 =', date_trans2_0
641 write(10, fmt=trim(fmt3))
'date_trans3 =', date_trans3_0
643 #if (ENHMOD==4 || ENHMOD==5)
644 write(10, fmt=trim(fmt3))
'enh_compr =', enh_compr
645 write(10, fmt=trim(fmt3))
'enh_shear =', enh_shear
647 #if ((ENHMOD==1 || ENHMOD==2 || ENHMOD==3 || ENHMOD==4) && MARGIN==3)
648 write(10, fmt=trim(fmt3))
'enh_shelf =', enh_shelf
650 write(10, fmt=trim(fmt1))
' '
653 write(10, fmt=trim(fmt1))
'mask_chasm file = '//mask_chasm_file
654 write(10, fmt=trim(fmt3))
'time_chasm_init =', time_chasm_init0
655 write(10, fmt=trim(fmt3))
'time_chasm_end =', time_chasm_end0
656 write(10, fmt=trim(fmt3))
'q_geo_chasm =', q_geo_chasm
657 write(10, fmt=trim(fmt3))
'erosion_chasm =', erosion_chasm
658 write(10, fmt=trim(fmt1))
' '
661 write(10, fmt=trim(fmt3))
'numdiff_H_t =', numdiff_h_t
662 write(10, fmt=trim(fmt3))
'tau_cts =', tau_cts
663 write(10, fmt=trim(fmt3))
'H_min =', h_min
664 write(10, fmt=trim(fmt3))
'vh_max =', vh_max
665 write(10, fmt=trim(fmt3))
'hd_min =', hd_min
666 write(10, fmt=trim(fmt3))
'hd_max =', hd_max
668 write(10, fmt=trim(fmt3))
'qbm_min =', qbm_min
669 #elif defined(QB_MIN) /* obsolete */
670 write(10, fmt=trim(fmt3))
'qb_min =', qb_min
673 write(10, fmt=trim(fmt3))
'qbm_max =', qbm_max
674 #elif defined(QB_MAX) /* obsolete */
675 write(10, fmt=trim(fmt3))
'qb_max =', qb_max
677 write(10, fmt=trim(fmt3))
'age_min =', age_min
678 write(10, fmt=trim(fmt3))
'age_max =', age_max
679 write(10, fmt=trim(fmt3))
'mean_accum =', mean_accum
681 write(10, fmt=trim(fmt3))
'age_diff =', agediff
683 write(10, fmt=trim(fmt1))
' '
685 write(10, fmt=trim(fmt2))
'GRID =', grid
686 write(10, fmt=trim(fmt2))
'DYNAMICS =', dynamics
687 #if (DYNAMICS==2 && defined(RATIO_SL_THRESH))
688 write(10, fmt=trim(fmt3))
'ratio_sl_thresh =', ratio_sl_thresh
690 write(10, fmt=trim(fmt2))
'CALCMOD =', calcmod
691 #if ( CALCMOD==-1 && defined(TEMP_CONST) )
692 write(10, fmt=trim(fmt3))
'TEMP_CONST =', temp_const
694 #if ( CALCMOD==-1 && defined(AGE_CONST) )
695 write(10, fmt=trim(fmt3))
'AGE_CONST =', age_const
697 #if ( CALCMOD==1 && defined(CTS_MELTING_FREEZING) )
698 write(10, fmt=trim(fmt2))
'CTS_MELTING_FREEZING =', cts_melting_freezing
700 write(10, fmt=trim(fmt2))
'FLOW_LAW =', flow_law
701 write(10, fmt=trim(fmt2))
'FIN_VISC =', fin_visc
702 write(10, fmt=trim(fmt2))
'MARGIN =', margin
704 write(10, fmt=trim(fmt2))
'MARINE_ICE_FORMATION =', marine_ice_formation
705 write(10, fmt=trim(fmt2))
'MARINE_ICE_CALVING =', marine_ice_calving
707 write(10, fmt=trim(fmt2))
'ICE_SHELF_CALVING =', ice_shelf_calving
709 write(10, fmt=trim(fmt2))
'ZS_EVOL =', zs_evol
710 write(10, fmt=trim(fmt2))
'CALCZS =', calczs
711 write(10, fmt=trim(fmt2))
'ADV_HOR =', adv_hor
712 write(10, fmt=trim(fmt2))
'ADV_VERT =', adv_vert
713 write(10, fmt=trim(fmt2))
'TOPOGRAD =', topograd
714 write(10, fmt=trim(fmt2))
'TSURFACE =', tsurface
715 write(10, fmt=trim(fmt2))
'ACC_UNIT =', acc_unit
716 write(10, fmt=trim(fmt2))
'ACC_PRESENT=', acc_present
717 write(10, fmt=trim(fmt2))
'ACCSURFACE =', accsurface
718 write(10, fmt=trim(fmt2))
'ABLSURFACE =', ablsurface
719 write(10, fmt=trim(fmt2))
'SEA_LEVEL =', sea_level
720 write(10, fmt=trim(fmt2))
'SLIDE_LAW =', slide_law
721 write(10, fmt=trim(fmt2))
'Q_GEO_MOD =', q_geo_mod
722 write(10, fmt=trim(fmt2))
'CHASM =', chasm
723 write(10, fmt=trim(fmt1))
' '
725 #if defined(OUT_TIMES)
726 write(10, fmt=trim(fmt2))
'OUT_TIMES =', out_times
728 write(10, fmt=trim(fmt2))
'OUTPUT =', output
729 write(10, fmt=trim(fmt2))
'OUTSER =', outser
730 #if ( OUTPUT==1 || OUTPUT==2 )
731 write(10, fmt=trim(fmt2))
'ERGDAT =', ergdat
733 write(10, fmt=trim(fmt2))
'NETCDF =', netcdf
734 write(10, fmt=trim(fmt1))
' '
735 #if ( OUTPUT==2 || OUTPUT==3 )
736 write(10, fmt=trim(fmt2))
'n_output =', n_output
738 write(10, fmt=trim(fmt3))
'time_output =', time_output0(n)
740 write(10, fmt=trim(fmt1))
' '
743 write(10, fmt=trim(fmt1))
'Program version and date: '//version//
' / '//date
745 close(10, status=
'keep')
749 year_zero = year_zero*year_sec
750 time_init = time_init0*year_sec
751 time_end = time_end0*year_sec
752 dtime = dtime0*year_sec
753 dtime_temp = dtime_temp0*year_sec
755 dtime_wss = dtime_wss0*year_sec
758 time_chasm_init = time_chasm_init0 *year_sec
759 time_chasm_end = time_chasm_end0 *year_sec
761 dtime_ser = dtime_ser0*year_sec
762 #if ( OUTPUT==1 || OUTPUT==3 )
763 dtime_out = dtime_out0*year_sec
765 #if ( OUTPUT==2 || OUTPUT==3 )
767 time_output(n) = time_output0(n)*year_sec
771 if (abs(dtime_temp/dtime-nint(dtime_temp/dtime)) > eps) &
772 stop
' sico_init: dtime_temp must be a multiple of dtime!'
775 if (abs(dtime_wss/dtime-nint(dtime_wss/dtime)) > eps) &
776 stop
' sico_init: dtime_wss must be a multiple of dtime!'
787 accum_present(j,i) = acc_present_val
801 accum_present(j,i) = accum_present(j,i) &
802 *(1.0e-03_dp/year_sec)*(rho_w/rho_i) &
803 *(1.0_dp/(1.0_dp-frac_dust))
807 accum_present(j,i) = accum_present(j,i)*(1.0e-03_dp/year_sec)
818 mean_accum = mean_accum*(1.0e-03_dp/year_sec)*(rho_w/rho_i) &
819 *(1.0_dp/(1.0_dp-frac_dust))
824 mean_accum = mean_accum*(1.0e-03_dp/year_sec)
831 open(24, iostat=ios, &
832 file=inpath//
'/'//trim(ch_domain_short)//
'/'//mask_present_file, &
833 recl=2048, status=
'old')
836 stop
' sico_init: Error when opening the mask file!'
838 do n=1, 6;
read(24, fmt=
'(a)') ch_dummy;
end do
841 read(24, fmt=trim(fmt4)) (maske_help(j,i), i=0,imax)
844 close(24, status=
'keep')
850 open(24, iostat=ios, &
851 file=inpath//
'/'//trim(ch_domain_short)//
'/'//mask_chasm_file, &
852 recl=2048, status=
'old')
854 if (ios /= 0) stop
' sico_init: Error when opening the chasm mask file!'
856 do n=1, 6;
read(24, fmt=
'(a)') ch_dummy;
end do
859 read(24, fmt=trim(fmt4)) (maske_chasm(j,i), i=0,imax)
862 close(24, status=
'keep')
870 stop
' sico_init: SEA_LEVEL==3 not allowed for nmars application!'
882 #if ( TSURFACE==5 || TSURFACE==6 )
884 open(21, iostat=ios, &
885 file=inpath//
'/'//trim(ch_domain_short)//
'/'//insol_ma_90n_file, &
887 if (ios /= 0) stop
' sico_init: Error when opening the insolation-data file!'
889 read(21, fmt=*) ch_dummy, insol_time_min, insol_time_stp, insol_time_max
891 if (ch_dummy /=
'#')
then
892 write(6, fmt=*)
' sico_init: insol_time_min, insol_time_stp, insol_time_max'
893 write(6, fmt=*)
' not defined in insolation-data file!'
897 ndata_insol = (insol_time_max-insol_time_min)/insol_time_stp
899 if (ndata_insol > 100000) &
900 stop
' sico_init: Too many data in insolation-data file!'
903 read(21, fmt=*) d_dummy, ecc_data(n), obl_data(n), cp_data(n), ave_data(n), insol_ma_90(n)
904 obl_data(n) = obl_data(n) *pi_180
905 ave_data(n) = ave_data(n) *pi_180
908 close(21, status=
'keep')
921 q_geo_normal(j,i) = q_geo *1.0e-03_dp
929 open(21, iostat=ios, &
930 file=inpath//
'/'//trim(ch_domain_short)//
'/'//q_geo_file, &
931 recl=8192, status=
'old')
933 if (ios /= 0) stop
' sico_init: Error when opening the qgeo file!'
935 do n=1, 6;
read(21, fmt=
'(a)') ch_dummy;
end do
938 read(21, fmt=*) (q_geo_normal(j,i), i=0,imax)
941 close(21, status=
'keep')
945 q_geo_normal(j,i) = q_geo_normal(j,i) *1.0e-03_dp
953 #if (REBOUND==0 || REBOUND==1)
955 kei = 0.0_dp; n_data_kei = 0; kei_r_max = 0.0_dp; kei_r_incr = 0.0_dp
966 #if (REBOUND==1 || REBOUND==2)
968 #if (TIME_LAG_MOD==1)
970 time_lag_asth = time_lag*year_sec
972 #elif (TIME_LAG_MOD==2)
974 open(29, iostat=ios, &
975 file=inpath//
'/'//trim(ch_domain_short)//
'/'//time_lag_file, &
976 recl=8192, status=
'old')
978 if (ios /= 0) stop
' sico_init: Error when opening the time-lag file!'
980 do n=1, 6;
read(29, fmt=
'(a)') ch_dummy;
end do
983 read(29, fmt=*) (time_lag_asth(j,i), i=0,imax)
986 close(29, status=
'keep')
988 time_lag_asth = time_lag_asth*year_sec
994 time_lag_asth = 0.0_dp
1002 #if (FLEX_RIG_MOD==1)
1004 flex_rig_lith = flex_rig
1006 #elif (FLEX_RIG_MOD==2)
1008 open(29, iostat=ios, &
1009 file=inpath//
'/'//trim(ch_domain_short)//
'/'//flex_rig_file, &
1010 recl=8192, status=
'old')
1012 if (ios /= 0) stop
' sico_init: Error when opening the flex-rig file!'
1014 do n=1, 6;
read(29, fmt=
'(a)') ch_dummy;
end do
1017 read(29, fmt=*) (flex_rig_lith(j,i), i=0,imax)
1020 close(29, status=
'keep')
1024 #elif (REBOUND==0 || REBOUND==1)
1026 flex_rig_lith = 0.0_dp
1040 call
boundary(time_init, dtime, dxi, deta, &
1041 delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1062 call
boundary(time_init, dtime, dxi, deta, &
1063 delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1085 stop
' sico_init: Parameter NETCDF must be either 1 or 2!'
1088 call
boundary(time_init, dtime, dxi, deta, &
1089 delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1091 q_b_tot = q_bm + q_tld
1097 flag_grounding_line_1 = .false.
1098 flag_grounding_line_2 = .false.
1100 flag_calving_front_1 = .false.
1101 flag_calving_front_2 = .false.
1103 #if (MARGIN==1 || MARGIN==2)
1112 if ( (maske(j,i)==0_i2b) &
1114 ( (maske(j,i+1)==3_i2b) &
1115 .or.(maske(j,i-1)==3_i2b) &
1116 .or.(maske(j+1,i)==3_i2b) &
1117 .or.(maske(j-1,i)==3_i2b) ) &
1119 flag_grounding_line_1(j,i) = .true.
1121 if ( (maske(j,i)==3_i2b) &
1123 ( (maske(j,i+1)==0_i2b) &
1124 .or.(maske(j,i-1)==0_i2b) &
1125 .or.(maske(j+1,i)==0_i2b) &
1126 .or.(maske(j-1,i)==0_i2b) ) &
1128 flag_grounding_line_2(j,i) = .true.
1130 if ( (maske(j,i)==3_i2b) &
1132 ( (maske(j,i+1)==2_i2b) &
1133 .or.(maske(j,i-1)==2_i2b) &
1134 .or.(maske(j+1,i)==2_i2b) &
1135 .or.(maske(j-1,i)==2_i2b) ) &
1137 flag_calving_front_1(j,i) = .true.
1139 if ( (maske(j,i)==2_i2b) &
1141 ( (maske(j,i+1)==3_i2b) &
1142 .or.(maske(j,i-1)==3_i2b) &
1143 .or.(maske(j+1,i)==3_i2b) &
1144 .or.(maske(j-1,i)==3_i2b) ) &
1146 flag_calving_front_2(j,i) = .true.
1154 if ( (maske(j,i)==2_i2b) &
1155 .and. (maske(j+1,i)==3_i2b) &
1157 flag_calving_front_2(j,i) = .true.
1160 if ( (maske(j,i)==2_i2b) &
1161 .and. (maske(j-1,i)==3_i2b) &
1163 flag_calving_front_2(j,i) = .true.
1170 if ( (maske(j,i)==2_i2b) &
1171 .and. (maske(j,i+1)==3_i2b) &
1173 flag_calving_front_2(j,i) = .true.
1176 if ( (maske(j,i)==2_i2b) &
1177 .and. (maske(j,i-1)==3_i2b) &
1179 flag_calving_front_2(j,i) = .true.
1184 stop
' sico_init: MARGIN must be either 1, 2 or 3!'
1189 #if (GRID==0 || GRID==1)
1193 dist_dxdy(jr,ir) = sqrt( (
real(ir,dp)*dxi)**2 + (
real(jr,dp)*deta)**2 )
1204 #if (DYNAMICS==1 || DYNAMICS==2)
1209 #if (MARGIN==3 || DYNAMICS==2)
1210 call
calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t, ii, jj, nn)
1225 stop
' sico_init: DYNAMICS must be either 0, 1 or 2!'
1228 call
calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
1232 #if (CALCMOD==0 || CALCMOD==-1)
1242 enth_t(kt,j,i) = enth_c(0,j,i)
1257 if ( (maske(j,i) == 0_i2b).and.(n_cts(j,i) == 1_i2b) )
then
1263 enth_t(kt,j,i) = enth_c(0,j,i)
1270 #elif (CALCMOD==2 || CALCMOD==3)
1280 enth_t(kt,j,i) = enth_c(0,j,i)
1288 stop
' sico_init: Parameter CALCMOD must be either -1, 0, 1, 2 or 3!'
1296 open(12, iostat=ios, &
1297 file=outpath//
'/'//trim(runname)//
'.ser', &
1300 stop
' sico_init: Error when opening the ser file!'
1305 1102
format(
' t(a) D_Ts(C) z_sl(m)',/, &
1306 ' H_max(m) zs_max(m) V(m^3)', &
1307 ' V_t(m^3) V_fw(m^3/a) Tbh_max(C)',/, &
1308 ' A(m^2) A_t(m^2) V_bm(m^3/a)', &
1309 ' V_tld(m^3/a) H_t_max(m) vs_max(m/a)')
1310 1103
format(
'----------------------------------------------------', &
1311 '---------------------------------------')
1317 open(13, iostat=ios, &
1318 file=outpath//
'/'//trim(runname)//
'.thk', &
1320 if (ios /= 0) stop
' sico_init: Error when opening the thk file!'
1325 1104
format(
' t(a) D_Ts(deg C) z_sl(m)',/, &
1326 ' H(m): JMAX+1 records [j = JMAX (-1) 0] with',/, &
1327 ' IMAX+1 values [i = 0 (1) IMAX] in each record')
1328 1105
format(
'----------------------------------------------------')
1338 allocate(lambda_core(n_core), phi_core(n_core), &
1339 x_core(n_core), y_core(n_core))
1341 lambda_core(1) = 0.0_dp
1342 phi_core(1) = 0.0_dp
1343 x_core(1) = 0.0_dp *1.0e+03_dp
1344 y_core(1) = 0.0_dp *1.0e+03_dp
1346 lambda_core(2) = 0.0_dp
1347 phi_core(2) = 0.0_dp
1348 x_core(2) = -150.0_dp *1.0e+03_dp
1349 y_core(2) = -290.0_dp *1.0e+03_dp
1351 lambda_core(3) = 0.0_dp
1352 phi_core(3) = 0.0_dp
1353 x_core(3) = -300.0_dp *1.0e+03_dp
1354 y_core(3) = -280.0_dp *1.0e+03_dp
1356 open(14, iostat=ios, &
1357 file=outpath//
'/'//trim(runname)//
'.core', &
1359 if (ios /= 0) stop
' sico_init: Error when opening the core file!'
1361 if (forcing_flag == 1)
then
1366 1106
format(
' t(a) D_Ts(C) z_sl(m)',/, &
1367 ' H_NP(m) H_C1(m) H_C2(m)',/, &
1368 ' v_NP(m/a) v_C1(m/a) v_C2(m/a)',/, &
1369 ' T_NP(C) T_C1(C) T_C2(C)',/, &
1370 ' Rb_NP(W/m2) Rb_C1(W/m2) Rb_C2(W/m2)')
1371 1107
format(
'----------------------------------------------------')
1382 flag_3d_output = .false.
1384 flag_3d_output = .true.
1386 stop
' sico_init: ERGDAT must be either 0 or 1!'
1390 call
output1(runname, time_init, delta_ts, glac_index, z_sl, &
1391 flag_3d_output, ndat2d, ndat3d)
1393 call
output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1394 flag_3d_output, ndat2d, ndat3d)
1396 stop
' sico_init: Parameter NETCDF must be either 1 or 2!'
1401 if (time_output(1) <= time_init+eps)
then
1404 flag_3d_output = .false.
1406 flag_3d_output = .true.
1408 stop
' sico_init: ERGDAT must be either 0 or 1!'
1412 call
output1(runname, time_init, delta_ts, glac_index, z_sl, &
1413 flag_3d_output, ndat2d, ndat3d)
1415 call
output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1416 flag_3d_output, ndat2d, ndat3d)
1418 stop
' sico_init: Parameter NETCDF must be either 1 or 2!'
1425 flag_3d_output = .false.
1428 call
output1(runname, time_init, delta_ts, glac_index, z_sl, &
1429 flag_3d_output, ndat2d, ndat3d)
1431 call
output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1432 flag_3d_output, ndat2d, ndat3d)
1434 stop
' sico_init: Parameter NETCDF must be either 1 or 2!'
1437 if (time_output(1) <= time_init+eps)
then
1439 flag_3d_output = .true.
1442 call
output1(runname, time_init, delta_ts, glac_index, z_sl, &
1443 flag_3d_output, ndat2d, ndat3d)
1445 call
output_nc(runname, time_init, delta_ts, glac_index, z_sl, &
1446 flag_3d_output, ndat2d, ndat3d)
1448 stop
' sico_init: Parameter NETCDF must be either 1 or 2!'
1454 stop
' sico_init: OUTPUT must be either 1, 2 or 3!'
1457 call
output2(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1459 call
output3(time_init, delta_ts, glac_index, z_sl)
1462 call
output4(time_init, dxi, deta, delta_ts, glac_index, z_sl)
subroutine calc_vz_ssa(z_sl, dxi, deta, dzeta_c)
Computation of the vertical velocity vz in the shallow shelf approximation.
subroutine calc_vxy_static()
Computation of the horizontal velocity vx, vy, the horizontal volume flux qx, qy etc. for static ice.
subroutine calc_vxy_sia(dzeta_c, dzeta_t)
Computation of the shear stresses txz, tyz, the effective shear stress sigma, the depth-averaged flui...
subroutine topography1(dxi, deta)
Definition of the initial surface and bedrock topography (including gradients) and of the horizontal ...
subroutine, public ice_mat_eqs_pars(RF_table, R_T_val, KAPPA_table, C_table, n_temp_min, n_temp_max, RHO_I_val, RHO_C_val, KAPPA_C_val, C_C_val)
Setting of required physical parameters.
subroutine output2(time, dxi, deta, delta_ts, glac_index, z_sl)
Writing of time-series data on file in ASCII format.
subroutine phys_para()
Reading of physical parameters.
Declarations of kind types for SICOPOLIS.
subroutine topography3(dxi, deta, z_sl, anfdatname)
Definition of the initial surface and bedrock topography (including gradients) and of the horizontal ...
subroutine output3(time, delta_ts, glac_index, z_sl)
Writing of time-series data of the ice thickness field on file in ASCII format.
subroutine topography3_nc(dxi, deta, z_sl, anfdatname)
Definition of the initial surface and bedrock topography (including gradients) and of the horizontal ...
subroutine calc_temp_melt()
Computation of the melting temperatures.
subroutine output_nc(runname, time, delta_ts, glac_index, z_sl, flag_3d_output, ndat2d, ndat3d)
Writing of time-slice files in NetCDF format.
subroutine calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t, ii, jj, nn)
Computation of the horizontal velocity vx, vy, the horizontal volume flux qx, qy and the flux across ...
Declarations of global variables for SICOPOLIS (for the ANT domain).
subroutine topography2(dxi, deta)
Definition of the initial surface and bedrock topography (including gradients) and of the horizontal ...
subroutine calc_vz_static()
Computation of the vertical velocity vz for static ice.
subroutine calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
Computation of all components of the strain-rate tensor, the full effective strain rate and the shear...
subroutine calc_enhance(time)
Computation of the flow enhancement factor.
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.
subroutine calc_vxy_b_sia(time, z_sl)
Computation of the basal horizontal velocity vx_b, vy_b in the shallow ice approximation.
subroutine sico_init(delta_ts, glac_index, mean_accum, dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, time, time_init, time_end, time_output, dxi, deta, dzeta_c, dzeta_t, dzeta_r, z_sl, dzsl_dtau, z_mar, ii, jj, nn, ndat2d, ndat3d, n_output, runname)
Initialisations for SICOPOLIS.
subroutine boundary(time, dtime, dxi, deta, delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
Computation of the surface temperature (must be less than 0 deg C!) and of the accumulation-ablation ...
subroutine read_kei()
Reading of the tabulated kei function.
subroutine init_temp_water_age()
Initial temperature, water content and age.
Material quantities of ice: Rate factor, heat conductivity, specific heat (heat capacity).
subroutine output4(time, dxi, deta, delta_ts, glac_index, z_sl)
Writing of time-series data of the deep ice cores on file in ASCII format.
Declarations of global variables for SICOPOLIS.
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_...
Conversion from temperature (temp) and water content (omega) to enthalpy (enth) and vice versa...
subroutine calc_vz_sia(dxi, deta, dzeta_c, dzeta_t)
Computation of the vertical velocity vz in the shallow ice approximation.
subroutine output1(runname, time, delta_ts, glac_index, z_sl, flag_3d_output, ndat2d, ndat3d)
Writing of time-slice files in binary format.
real(dp) function, public enth_fct_temp_omega(temp_val, omega_val)
Enthalpy as a function of temperature and water content.