50 subroutine sico_init(delta_ts, glac_index, &
52 dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, &
53 time, time_init, time_end, time_output, &
54 dxi, deta, dzeta_c, dzeta_t, dzeta_r, &
55 z_sl, dzsl_dtau, z_mar, &
57 ndat2d, ndat3d, n_output, &
79 integer(i4b),
intent(out) :: ii((imax+1)*(jmax+1)), &
80 jj((imax+1)*(jmax+1)), &
82 integer(i4b),
intent(out) :: ndat2d, ndat3d
83 integer(i4b),
intent(out) :: n_output
84 real(dp),
intent(out) :: delta_ts, glac_index
85 real(dp),
intent(out) :: mean_accum
86 real(dp),
intent(out) :: dtime, dtime_temp, dtime_wss, &
88 real(dp),
intent(out) :: time, time_init, time_end, time_output(100)
89 real(dp),
intent(out) :: dxi, deta, dzeta_c, dzeta_t, dzeta_r
90 real(dp),
intent(out) :: z_sl, dzsl_dtau, z_mar
91 character(len=100),
intent(out) :: runname
93 integer(i4b) :: i, j, kc, kt, kr, m, n, ir, jr
96 integer(i2b),
dimension(0:JMAX,0:IMAX) :: maske_ref
97 real(dp) :: dtime0, dtime_temp0, dtime_wss0, dtime_out0, dtime_ser0
98 real(dp) :: time_init0, time_end0, time_output0(100)
100 character(len=100) :: anfdatname
101 character(len=256) :: filename_with_path
102 character(len=256) :: shell_command
103 character :: ch_dummy
104 logical :: flag_init_output, flag_3d_output
106 character(len=64),
parameter :: fmt1 =
'(a)', &
111 character(len= 8) :: ch_imax
112 character(len=128) :: fmt4
114 write(ch_imax, fmt=
'(i8)') imax
115 write(fmt4, fmt=
'(a)')
'('//trim(adjustl(ch_imax))//
'(i1),i1)' 117 write(unit=6, fmt=
'(a)')
' ' 118 write(unit=6, fmt=
'(a)')
' -------- sico_init --------' 119 write(unit=6, fmt=
'(a)')
' ' 131 #elif (defined(EMTP2SGE)) 139 #elif (defined(NHEM)) 143 #elif (defined(SCAND)) 147 #elif (defined(TIBET)) 151 #elif (defined(NMARS)) 155 #elif (defined(SMARS)) 167 stop
' >>> sico_init: No valid domain specified!' 188 #if (CALCTHK==3 || CALCTHK==6 || MARGIN==3 || DYNAMICS==2) 189 call lis_initialize(ierr)
205 if ( trim(version) /= trim(sico_version) ) &
206 stop
' >>> sico_init: SICOPOLIS version not compatible with header file!' 210 #if (!defined(DYNAMICS)) 211 stop
' >>> sico_init: DYNAMICS not defined in the header file!' 214 #if (!defined(CALCMOD)) 215 stop
' >>> sico_init: CALCMOD not defined in the header file!' 218 #if (defined(ENTHMOD)) 219 write(6, fmt=
'(a)')
' >>> sico_init: ENTHMOD must not be defined any more.' 220 write(6, fmt=
'(a)')
' Please update your header file!' 227 #if (GRID==0 || GRID==1) 229 if (approx_equal(dx, 40.0_dp,
eps_sp_dp))
then 231 if ((imax /= 150).or.(jmax /= 70))
then 232 stop
' >>> sico_init: IMAX and/or JMAX wrong!' 235 else if (approx_equal(dx, 20.0_dp,
eps_sp_dp))
then 237 if ((imax /= 300).or.(jmax /= 140))
then 238 stop
' >>> sico_init: IMAX and/or JMAX wrong!' 241 else if (approx_equal(dx, 10.0_dp,
eps_sp_dp))
then 243 if ((imax /= 600).or.(jmax /= 280))
then 244 stop
' >>> sico_init: IMAX and/or JMAX wrong!' 248 stop
' >>> sico_init: DX wrong!' 253 stop
' >>> sico_init: GRID==2 not allowed for the Scandinavia application!' 261 #if (CALCMOD==0 || CALCMOD==2 || CALCMOD==3 || CALCMOD==-1) 264 write(6, fmt=
'(a)')
' >>> sico_init: For options CALCMOD==0, 2, 3 or -1,' 265 write(6, fmt=
'(a)')
' the separate kt domain is redundant.' 266 write(6, fmt=
'(a)')
' Therefore, consider setting KTMAX to 2.' 267 write(6, fmt=
'(a)')
' ' 276 #if (TSURFACE == 5 && ACCSURFACE != 5) 277 stop
' >>> sico_init: Options TSURFACE==5 and ACCSURFACE==5 must be used together!' 280 #if (TSURFACE != 5 && ACCSURFACE == 5) 281 stop
' >>> sico_init: Options TSURFACE==5 and ACCSURFACE==5 must be used together!' 288 stop
' >>> sico_init: Option ADV_HOR==1 (central differences) not defined!' 296 #if ((MARGIN==3 && DYNAMICS==1) || DYNAMICS==2) /* requires SSA or SStA */ 298 write(6, fmt=
'(a)')
' >>> sico_init: Distortion correction not implemented' 299 write(6, fmt=
'(a)')
' for the shallow shelf approximation (SSA)' 300 write(6, fmt=
'(a)')
' or the shelfy stream approximation (SStA)' 301 write(6, fmt=
'(a)')
' -> GRID==0 required!' 312 #elif (TSURFACE == 5) 321 dtime_temp0 = dtime_temp0
323 dtime_wss0 = dtime_wss0
328 dzeta_c = 1.0_dp/
real(kcmax,dp)
329 dzeta_t = 1.0_dp/
real(ktmax,dp)
330 dzeta_r = 1.0_dp/
real(krmax,dp)
339 if (deform >=
eps)
then 419 if ((i <= imax).and.(j <= jmax))
then 431 anfdatname = anfdatname
433 #if (defined(YEAR_ZERO)) 439 time_init0 = time_init0
440 time_end0 = time_end0
441 dtime_ser0 = dtime_ser0
443 #if (OUTPUT==1 || OUTPUT==3) 444 dtime_out0 = dtime_out0
447 #if (OUTPUT==2 || OUTPUT==3) 449 time_output0( 1) = time_out0_01
450 time_output0( 2) = time_out0_02
451 time_output0( 3) = time_out0_03
452 time_output0( 4) = time_out0_04
453 time_output0( 5) = time_out0_05
454 time_output0( 6) = time_out0_06
455 time_output0( 7) = time_out0_07
456 time_output0( 8) = time_out0_08
457 time_output0( 9) = time_out0_09
458 time_output0(10) = time_out0_10
459 time_output0(11) = time_out0_11
460 time_output0(12) = time_out0_12
461 time_output0(13) = time_out0_13
462 time_output0(14) = time_out0_14
463 time_output0(15) = time_out0_15
464 time_output0(16) = time_out0_16
465 time_output0(17) = time_out0_17
466 time_output0(18) = time_out0_18
467 time_output0(19) = time_out0_19
468 time_output0(20) = time_out0_20
473 shell_command =
'if [ ! -d' 474 shell_command = trim(shell_command)//
' '//outpath
475 shell_command = trim(shell_command)//
' '//
'] ; then mkdir' 476 shell_command = trim(shell_command)//
' '//outpath
477 shell_command = trim(shell_command)//
' '//
'; fi' 478 call system(trim(shell_command))
481 filename_with_path = trim(outpath)//
'/'//trim(runname)//
'.log' 483 open(10, iostat=ios, file=trim(filename_with_path), status=
'new')
485 if (ios /= 0) stop
' >>> Error when opening the log file!' 487 write(10, fmt=trim(fmt1))
'Computational domain:' 489 write(10, fmt=trim(fmt1))
' ' 491 write(10, fmt=trim(fmt2a))
'GRID = ', grid
492 write(10, fmt=trim(fmt1))
' ' 494 write(10, fmt=trim(fmt2))
'imax =', imax
495 write(10, fmt=trim(fmt2))
'jmax =', jmax
496 write(10, fmt=trim(fmt2))
'kcmax =', kcmax
497 write(10, fmt=trim(fmt2))
'ktmax =', ktmax
498 write(10, fmt=trim(fmt2))
'krmax =', krmax
499 write(10, fmt=trim(fmt1))
' ' 501 write(10, fmt=trim(fmt3))
'a =',
aa 502 write(10, fmt=trim(fmt1))
' ' 504 #if (GRID==0 || GRID==1) 505 write(10, fmt=trim(fmt3))
'x0 =', x0
506 write(10, fmt=trim(fmt3))
'y0 =', y0
507 write(10, fmt=trim(fmt3))
'dx =', dx
509 stop
' >>> sico_init: GRID==2 not allowed for this application!' 511 write(10, fmt=trim(fmt1))
' ' 513 write(10, fmt=trim(fmt3))
'year_zero =',
year_zero 514 write(10, fmt=trim(fmt3))
'time_init =', time_init0
515 write(10, fmt=trim(fmt3))
'time_end =', time_end0
516 write(10, fmt=trim(fmt3))
'dtime =', dtime0
517 write(10, fmt=trim(fmt3))
'dtime_temp =', dtime_temp0
519 write(10, fmt=trim(fmt3))
'dtime_wss =', dtime_wss0
521 write(10, fmt=trim(fmt1))
' ' 523 write(10, fmt=trim(fmt2))
'ANF_DAT =', anf_dat
524 write(10, fmt=trim(fmt1))
'zs_present file = '//zs_present_file
526 #if (defined(ZB_PRESENT_FILE)) 527 write(10, fmt=trim(fmt1))
'zb_present file = '//zb_present_file
529 write(10, fmt=trim(fmt1))
'zl_present file = '//zl_present_file
531 write(10, fmt=trim(fmt1))
'zl0 file = '//zl0_file
532 write(10, fmt=trim(fmt1))
'mask_present file = '//mask_present_file
533 #if (ANF_DAT==1 && defined(TEMP_INIT)) 534 write(10, fmt=trim(fmt2))
'TEMP_INIT =', temp_init
537 write(10, fmt=trim(fmt1))
'Initial-value file = '//anfdatname
538 write(10, fmt=trim(fmt1))
'Path to initial-value file = '//anfdatpath
540 write(10, fmt=trim(fmt1))
' ' 542 write(10, fmt=trim(fmt1))
'Physical-parameter file = '//phys_para_file
543 write(10, fmt=trim(fmt1))
' ' 545 #if (CALCTHK==2 || CALCTHK==3 || CALCTHK==5 || CALCTHK==6) 546 write(10, fmt=trim(fmt3))
'ovi_weight =', ovi_weight
547 #if (CALCTHK==2 || CALCTHK==5) 548 write(10, fmt=trim(fmt3))
'omega_sor =', omega_sor
550 write(10, fmt=trim(fmt2a))
'iter_max_sor = ', iter_max_sor
553 write(10, fmt=trim(fmt1))
' ' 556 write(10, fmt=trim(fmt1))
'temp_mm_present file = '//temp_mm_present_file
558 write(10, fmt=trim(fmt3))
'delta_ts0 =', delta_ts0
560 write(10, fmt=trim(fmt3))
'sine_amplit =', sine_amplit
561 write(10, fmt=trim(fmt3))
'sine_period =', sine_period
563 write(10, fmt=trim(fmt1))
'GRIP file = '//grip_temp_file
564 write(10, fmt=trim(fmt3))
'grip_temp_fact =', grip_temp_fact
566 write(10, fmt=trim(fmt1))
'Glacial-index file = '//glac_ind_file
567 write(10, fmt=trim(fmt1))
'temp_mm_anom file = '//temp_mm_anom_file
568 write(10, fmt=trim(fmt3))
'temp_mm_anom fact = ', temp_mm_anom_fact
571 write(10, fmt=trim(fmt1))
'precip_mm_present file = '//precip_mm_present_file
573 write(10, fmt=trim(fmt3))
'accfact =', accfact
574 #elif (ACCSURFACE==2 || ACCSURFACE==3) 575 write(10, fmt=trim(fmt3))
'gamma_s =', gamma_s
576 #elif (ACCSURFACE==5) 577 write(10, fmt=trim(fmt1))
'precip_mm_anom file = '//precip_mm_anom_file
578 write(10, fmt=trim(fmt3))
'precip_mm_anom fact = ', precip_mm_anom_fact
580 #if (ACCSURFACE <= 3) 581 write(10, fmt=trim(fmt2))
'ELEV_DESERT =', elev_desert
582 #if (ELEV_DESERT == 1) 583 write(10, fmt=trim(fmt3))
'gamma_p =', gamma_p
584 write(10, fmt=trim(fmt3))
'zs_thresh =', zs_thresh
589 write(10, fmt=trim(fmt3))
'lambda_lti =', lambda_lti
590 write(10, fmt=trim(fmt3))
'temp_lti =', temp_lti
593 write(10, fmt=trim(fmt2a))
'SEA_LEVEL = ', sea_level
595 write(10, fmt=trim(fmt3))
'z_sl0 =', z_sl0
597 write(10, fmt=trim(fmt1))
'sea-level file = '//sea_level_file
599 write(10, fmt=trim(fmt1))
' ' 602 #if (MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3) 603 write(10, fmt=trim(fmt3))
'z_mar =', z_mar
604 write(10, fmt=trim(fmt1))
' ' 605 #elif (MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5 \ 606 || marine_ice_calving==6 || marine_ice_calving==7)
607 write(10, fmt=trim(fmt3))
'fact_z_mar =', fact_z_mar
608 write(10, fmt=trim(fmt1))
' ' 609 #elif (MARINE_ICE_FORMATION==2 && MARINE_ICE_CALVING==9) 610 write(10, fmt=trim(fmt3))
'calv_uw_coeff =', calv_uw_coeff
611 write(10, fmt=trim(fmt3))
'r1_calv_uw =', r1_calv_uw
612 write(10, fmt=trim(fmt3))
'r2_calv_uw =', r2_calv_uw
613 write(10, fmt=trim(fmt1))
' ' 616 #if (ICE_SHELF_CALVING==2) 617 write(10, fmt=trim(fmt3))
'H_calv =', h_calv
618 write(10, fmt=trim(fmt1))
' ' 622 #if (defined(BASAL_HYDROLOGY)) 623 write(10, fmt=trim(fmt2a))
'BASAL_HYDROLOGY = ', basal_hydrology
625 write(10, fmt=trim(fmt2a))
'SLIDE_LAW = ', slide_law
626 write(10, fmt=trim(fmt3))
'c_slide =',
c_slide 627 write(10, fmt=trim(fmt3))
'gamma_slide =', gamma_slide
628 write(10, fmt=trim(fmt2a))
'p_weert = ',
p_weert 629 write(10, fmt=trim(fmt2a))
'q_weert = ',
q_weert 630 #if (defined(TIME_RAMP_UP_SLIDE)) 631 write(10, fmt=trim(fmt3))
'time_ramp_up_slide =', time_ramp_up_slide
634 write(10, fmt=trim(fmt3))
'red_pres_limit_fact =', red_pres_limit_fact
636 #if (BASAL_HYDROLOGY==1) 637 write(10, fmt=trim(fmt3))
'c_Hw_slide =', c_hw_slide
638 write(10, fmt=trim(fmt3))
'Hw0_slide =', hw0_slide
640 write(10, fmt=trim(fmt1))
' ' 642 write(10, fmt=trim(fmt2a))
'Q_GEO_MOD = ', q_geo_mod
644 write(10, fmt=trim(fmt3))
'q_geo =',
q_geo 646 write(10, fmt=trim(fmt1))
'q_geo file = '//q_geo_file
648 write(10, fmt=trim(fmt1))
' ' 650 #if (defined(MARINE_ICE_BASAL_MELTING)) 651 write(10, fmt=trim(fmt2))
'MARINE_ICE_BASAL_MELTING =', marine_ice_basal_melting
652 #if (MARINE_ICE_BASAL_MELTING==2 || MARINE_ICE_BASAL_MELTING==3) 653 write(10, fmt=trim(fmt3))
'qbm_marine =', qbm_marine
655 write(10, fmt=trim(fmt1))
' ' 659 write(10, fmt=trim(fmt2))
'FLOATING_ICE_BASAL_MELTING =', floating_ice_basal_melting
660 #if (FLOATING_ICE_BASAL_MELTING<=3) 661 write(10, fmt=trim(fmt3))
'qbm_float_1 =', qbm_float_1
662 write(10, fmt=trim(fmt3))
'qbm_float_2 =', qbm_float_2
664 write(10, fmt=trim(fmt3))
'qbm_float_3 =', qbm_float_3
665 write(10, fmt=trim(fmt3))
'z_abyss =', z_abyss
666 #if (FLOATING_ICE_BASAL_MELTING==4) 667 write(10, fmt=trim(fmt3))
'temp_ocean =', temp_ocean
668 write(10, fmt=trim(fmt3))
'Omega_qbm =', omega_qbm
669 write(10, fmt=trim(fmt3))
'alpha_qbm =', alpha_qbm
671 #if (FLOATING_ICE_BASAL_MELTING==4 || FLOATING_ICE_BASAL_MELTING==5) 672 write(10, fmt=trim(fmt3))
'H_w_0 =', h_w_0
674 write(10, fmt=trim(fmt1))
' ' 677 write(10, fmt=trim(fmt2))
'REBOUND =', rebound
679 write(10, fmt=trim(fmt3))
'frac_llra =', frac_llra
681 #if (REBOUND==1 || REBOUND==2) 682 write(10, fmt=trim(fmt2))
'TIME_LAG_MOD =', time_lag_mod
683 #if (TIME_LAG_MOD==1) 684 write(10, fmt=trim(fmt3))
'time_lag =', time_lag
685 #elif (TIME_LAG_MOD==2) 686 write(10, fmt=trim(fmt1))
'time_lag_file = '//time_lag_file
688 stop
' >>> sico_init: TIME_LAG_MOD must be either 1 or 2!' 692 write(10, fmt=trim(fmt2))
'FLEX_RIG_MOD =', flex_rig_mod
693 #if (FLEX_RIG_MOD==1) 694 write(10, fmt=trim(fmt3))
'flex_rig =', flex_rig
695 #elif (FLEX_RIG_MOD==2) 696 write(10, fmt=trim(fmt1))
'flex_rig_file = '//flex_rig_file
698 stop
' >>> sico_init: FLEX_RIG_MOD must be either 1 or 2!' 701 write(10, fmt=trim(fmt2))
'Q_LITHO =', q_litho
702 write(10, fmt=trim(fmt1))
' ' 705 write(10, fmt=trim(fmt3))
'gr_size =', gr_size
706 write(10, fmt=trim(fmt1))
' ' 709 write(10, fmt=trim(fmt3))
'sigma_res =', sigma_res
710 write(10, fmt=trim(fmt1))
' ' 713 write(10, fmt=trim(fmt2))
'ENHMOD =', enhmod
714 #if (ENHMOD==1 || ENHMOD==2 || ENHMOD==3) 715 write(10, fmt=trim(fmt3))
'enh_fact =', enh_fact
717 #if (ENHMOD==2 || ENHMOD==3) 718 write(10, fmt=trim(fmt3))
'enh_intg =', enh_intg
721 write(10, fmt=trim(fmt3))
'age_trans =', age_trans_0
724 write(10, fmt=trim(fmt3))
'date_trans1 =', date_trans1_0
725 write(10, fmt=trim(fmt3))
'date_trans2 =', date_trans2_0
726 write(10, fmt=trim(fmt3))
'date_trans3 =', date_trans3_0
728 #if (ENHMOD==4 || ENHMOD==5) 729 write(10, fmt=trim(fmt3))
'enh_compr =', enh_compr
730 write(10, fmt=trim(fmt3))
'enh_shear =', enh_shear
732 #if ((ENHMOD==1 || ENHMOD==2 || ENHMOD==3 || ENHMOD==4) && MARGIN==3) 733 write(10, fmt=trim(fmt3))
'enh_shelf =', enh_shelf
735 write(10, fmt=trim(fmt1))
' ' 737 write(10, fmt=trim(fmt3))
'numdiff_H_t =', numdiff_h_t
738 write(10, fmt=trim(fmt3))
'tau_cts =', tau_cts
739 write(10, fmt=trim(fmt3))
'vh_max =', vh_max
740 write(10, fmt=trim(fmt3))
'hd_min =', hd_min
741 write(10, fmt=trim(fmt3))
'hd_max =', hd_max
742 #if (defined(QBM_MIN)) 743 write(10, fmt=trim(fmt3))
'qbm_min =', qbm_min
744 #elif (defined(QB_MIN)) /* obsolete */ 745 write(10, fmt=trim(fmt3))
'qb_min =', qb_min
747 #if (defined(QBM_MAX)) 748 write(10, fmt=trim(fmt3))
'qbm_max =', qbm_max
749 #elif (defined(QB_MAX)) /* obsolete */ 750 write(10, fmt=trim(fmt3))
'qb_max =', qb_max
752 write(10, fmt=trim(fmt3))
'age_min =', age_min
753 write(10, fmt=trim(fmt3))
'age_max =', age_max
754 write(10, fmt=trim(fmt3))
'mean_accum =', mean_accum
756 write(10, fmt=trim(fmt3))
'age_diff =', agediff
758 write(10, fmt=trim(fmt1))
' ' 760 write(10, fmt=trim(fmt2a))
'DYNAMICS = ', dynamics
761 #if ((DYNAMICS==1 && MARGIN==3) || DYNAMICS==2) 762 #if (defined(LIS_OPTS)) 763 write(10, fmt=trim(fmt1))
'lis_opts = '//lis_opts
765 #if (defined(N_ITER_SSA)) 766 write(10, fmt=trim(fmt2a))
'n_iter_ssa = ', n_iter_ssa
768 #if (defined(RELAX_FACT_SSA)) 769 write(10, fmt=trim(fmt3))
'relax_fact_ssa =', relax_fact_ssa
772 #if (DYNAMICS==2 && defined(RATIO_SL_THRESH)) 773 write(10, fmt=trim(fmt3))
'ratio_sl_thresh =', ratio_sl_thresh
775 #if (DYNAMICS==2 && defined(SSTA_SIA_WEIGH_FCT)) 776 write(10, fmt=trim(fmt2a))
'SSTA_SIA_WEIGH_FCT = ', ssta_sia_weigh_fct
778 write(10, fmt=trim(fmt1))
' ' 780 write(10, fmt=trim(fmt2))
'CALCMOD =', calcmod
781 #if (CALCMOD==-1 && defined(TEMP_CONST)) 782 write(10, fmt=trim(fmt3))
'TEMP_CONST =', temp_const
784 #if (CALCMOD==-1 && defined(AGE_CONST)) 785 write(10, fmt=trim(fmt3))
'AGE_CONST =', age_const
787 #if (CALCMOD==1 && defined(CTS_MELTING_FREEZING)) 788 write(10, fmt=trim(fmt2))
'CTS_MELTING_FREEZING =', cts_melting_freezing
790 write(10, fmt=trim(fmt2))
'FLOW_LAW =', flow_law
791 write(10, fmt=trim(fmt2))
'FIN_VISC =', fin_visc
792 write(10, fmt=trim(fmt2))
'MARGIN =', margin
794 write(10, fmt=trim(fmt2))
'MARINE_ICE_FORMATION =', marine_ice_formation
795 write(10, fmt=trim(fmt2))
'MARINE_ICE_CALVING =', marine_ice_calving
797 write(10, fmt=trim(fmt2))
'ICE_SHELF_CALVING =', ice_shelf_calving
799 #if (defined(THK_EVOL)) 800 write(10, fmt=trim(fmt2))
'THK_EVOL =', thk_evol
802 stop
' >>> sico_init: Define THK_EVOL in header file!' 804 #if (defined(CALCTHK)) 805 write(10, fmt=trim(fmt2))
'CALCTHK =', calcthk
807 stop
' >>> sico_init: Define CALCTHK in header file!' 809 write(10, fmt=trim(fmt2))
'ADV_HOR =', adv_hor
810 write(10, fmt=trim(fmt2))
'ADV_VERT =', adv_vert
811 write(10, fmt=trim(fmt2))
'TOPOGRAD =', topograd
812 write(10, fmt=trim(fmt2))
'TSURFACE =', tsurface
813 write(10, fmt=trim(fmt2))
'ACCSURFACE =', accsurface
815 write(10, fmt=trim(fmt2))
'PRECIP_ANOM_INTERPOL =', precip_anom_interpol
817 write(10, fmt=trim(fmt2))
'SOLID_PRECIP =', solid_precip
818 write(10, fmt=trim(fmt2))
'ABLSURFACE =', ablsurface
819 #if (defined(MB_ACCOUNT)) 820 write(10, fmt=trim(fmt2a))
'MB_ACCOUNT = ', mb_account
822 write(10, fmt=trim(fmt1))
' ' 824 #if (defined(OUT_TIMES)) 825 write(10, fmt=trim(fmt2a))
'OUT_TIMES = ', out_times
827 #if (defined(OUTPUT_INIT)) 828 write(10, fmt=trim(fmt2a))
'OUTPUT_INIT = ', output_init
830 write(10, fmt=trim(fmt2a))
'OUTPUT = ', output
831 #if (OUTPUT==1 || OUTPUT==3) 832 write(10, fmt=trim(fmt3))
'dtime_out =' , dtime_out0
834 write(10, fmt=trim(fmt3))
'dtime_ser =' , dtime_ser0
835 #if (OUTPUT==1 || OUTPUT==2) 836 write(10, fmt=trim(fmt2a))
'ERGDAT = ', ergdat
838 write(10, fmt=trim(fmt2a))
'NETCDF = ', netcdf
839 #if (OUTPUT==2 || OUTPUT==3) 840 write(10, fmt=trim(fmt2a))
'n_output = ', n_output
842 write(10, fmt=trim(fmt3))
'time_output =' , time_output0(n)
845 write(10, fmt=trim(fmt1))
' ' 847 write(10, fmt=trim(fmt1))
'Program version and date: '//version//
' / '//date
849 close(10, status=
'keep')
854 time_init = time_init0*year_sec
855 time_end = time_end0*year_sec
856 dtime = dtime0*year_sec
857 dtime_temp = dtime_temp0*year_sec
859 dtime_wss = dtime_wss0*year_sec
861 dtime_ser = dtime_ser0*year_sec
862 #if (OUTPUT==1 || OUTPUT==3) 863 dtime_out = dtime_out0*year_sec
865 #if (OUTPUT==2 || OUTPUT==3) 867 time_output(n) = time_output0(n)*year_sec
871 if (.not.approx_integer_multiple(dtime_temp, dtime,
eps_sp_dp)) &
872 stop
' >>> sico_init: dtime_temp must be a multiple of dtime!' 875 if (.not.approx_integer_multiple(dtime_wss, dtime,
eps_sp_dp)) &
876 stop
' >>> sico_init: dtime_wss must be a multiple of dtime!' 879 if (.not.approx_integer_multiple(dtime_ser, dtime,
eps_sp_dp)) &
880 stop
' >>> sico_init: dtime_ser must be a multiple of dtime!' 882 #if (OUTPUT==1 || OUTPUT==3) 883 if (.not.approx_integer_multiple(dtime_out, dtime,
eps_sp_dp)) &
884 stop
' >>> sico_init: dtime_out must be a multiple of dtime!' 891 #if (GRID==0 || GRID==1) 894 trim(precip_mm_present_file)
896 open(21, iostat=ios, file=trim(filename_with_path), recl=
rcl1, status=
'old')
900 stop
' >>> sico_init: GRID==2 not allowed for the Scandinavia application!' 904 if (ios /= 0) stop
' >>> sico_init: Error when opening the precipitation file!' 906 do m=1, 6;
read(21, fmt=
'(a)') ch_dummy;
end do 909 do m=1, 3;
read(21, fmt=
'(a)') ch_dummy;
end do 915 close(21, status=
'keep')
926 #if (GRID==0 || GRID==1) 929 trim(precip_mm_anom_file)
931 open(21, iostat=ios, file=trim(filename_with_path), recl=
rcl1, status=
'old')
935 if (ios /= 0) stop
' >>> sico_init: Error when opening the precip anomaly file!' 937 do m=1, 6;
read(21, fmt=
'(a)') ch_dummy;
end do 940 do m=1, 3;
read(21, fmt=
'(a)') ch_dummy;
end do 946 close(21, status=
'keep')
953 #if (PRECIP_ANOM_INTERPOL==1) 957 #elif (PRECIP_ANOM_INTERPOL==2) 962 stop
' >>> sico_init: Wrong value of switch PRECIP_ANOM_INTERPOL!' 972 mean_accum = mean_accum*(1.0e-03_dp/year_sec)*(
rho_w/
rho)
977 #if (GRID==0 || GRID==1) 980 trim(mask_present_file)
982 open(24, iostat=ios, file=trim(filename_with_path), recl=
rcl2, status=
'old')
986 stop
' >>> sico_init: GRID==2 not allowed for the Scandinavia application!' 990 if (ios /= 0) stop
' >>> sico_init: Error when opening the mask file!' 992 do m=1, 6;
read(24, fmt=
'(a)') ch_dummy;
end do 995 read(24, fmt=trim(fmt4)) (maske_ref(j,i), i=0,imax)
998 close(24, status=
'keep')
1002 #if (GRID==0 || GRID==1) 1004 filename_with_path = trim(inpath)//
'/'//trim(
ch_domain_short)//
'/'// &
1005 trim(temp_mm_present_file)
1007 open(21, iostat=ios, file=trim(filename_with_path), recl=
rcl1, status=
'old')
1011 stop
' >>> sico_init: GRID==2 not allowed for the Scandinavia application!' 1015 if (ios /= 0) stop
' >>> sico_init: Error when opening the temperature file!' 1017 do m=1, 6;
read(21, fmt=
'(a)') ch_dummy;
end do 1020 do m=1, 3;
read(21, fmt=
'(a)') ch_dummy;
end do 1026 close(21, status=
'keep')
1032 #if (GRID==0 || GRID==1) 1034 filename_with_path = trim(inpath)//
'/'//trim(
ch_domain_short)//
'/'// &
1035 trim(temp_mm_anom_file)
1037 open(21, iostat=ios, file=trim(filename_with_path), recl=
rcl1, status=
'old')
1041 if (ios /= 0) stop
' >>> sico_init: Error when opening the temperature anomaly file!' 1043 do m=1, 6;
read(21, fmt=
'(a)') ch_dummy;
end do 1046 do m=1, 3;
read(21, fmt=
'(a)') ch_dummy;
end do 1052 close(21, status=
'keep')
1061 #if (GRID==0 || GRID==1) 1063 filename_with_path = trim(inpath)//
'/'//trim(
ch_domain_short)//
'/'// &
1064 trim(zs_present_file)
1066 open(21, iostat=ios, file=trim(filename_with_path), recl=
rcl1, status=
'old')
1070 stop
' >>> sico_init: GRID==2 not allowed for the Scandinavia application!' 1074 if (ios /= 0) stop
' >>> sico_init: Error when opening the zs file!' 1076 do m=1, 6;
read(21, fmt=
'(a)') ch_dummy;
end do 1079 read(21, fmt=*) (
zs_ref(j,i), i=0,imax)
1082 close(21, status=
'keep')
1089 if (maske_ref(j,i) >= 2)
zs_ref(j,i) = 0.0_dp
1097 filename_with_path = trim(inpath)//
'/general/'//trim(grip_temp_file)
1099 open(21, iostat=ios, file=trim(filename_with_path), status=
'old')
1101 if (ios /= 0) stop
' >>> sico_init: Error when opening the data file for delta_ts!' 1105 if (ch_dummy /=
'#')
then 1106 write(6, fmt=*)
' >>> sico_init: grip_time_min, grip_time_stp, grip_time_max' 1107 write(6, fmt=*)
' not defined in data file for delta_ts!' 1116 read(21, fmt=*) d_dummy,
griptemp(n)
1119 close(21, status=
'keep')
1127 filename_with_path = trim(inpath)//
'/general/'//trim(glac_ind_file)
1129 open(21, iostat=ios, file=trim(filename_with_path), status=
'old')
1131 if (ios /= 0) stop
' >>> sico_init: Error when opening the glacial-index file!' 1135 if (ch_dummy /=
'#')
then 1136 write(6, fmt=*)
' >>> sico_init: gi_time_min, gi_time_stp, gi_time_max' 1137 write(6, fmt=*)
' not defined in glacial-index file!' 1149 close(21, status=
'keep')
1157 filename_with_path = trim(inpath)//
'/general/'//trim(sea_level_file)
1159 open(21, iostat=ios, file=trim(filename_with_path), status=
'old')
1161 if (ios /= 0) stop
' >>> sico_init: Error when opening the data file for z_sl!' 1165 if (ch_dummy /=
'#')
then 1166 write(6, fmt=*)
' >>> sico_init: specmap_time_min, specmap_time_stp, specmap_time_max' 1167 write(6, fmt=*)
' not defined in data file for z_sl!' 1179 close(21, status=
'keep')
1195 #elif (Q_GEO_MOD==2) 1199 filename_with_path = trim(inpath)//
'/'//trim(
ch_domain_short)//
'/'// &
1202 open(21, iostat=ios, file=trim(filename_with_path), recl=
rcl1, status=
'old')
1204 if (ios /= 0) stop
' >>> sico_init: Error when opening the qgeo file!' 1206 do m=1, 6;
read(21, fmt=
'(a)') ch_dummy;
end do 1209 read(21, fmt=*) (
q_geo(j,i), i=0,imax)
1212 close(21, status=
'keep')
1224 #if (REBOUND==0 || REBOUND==1) 1237 #if (REBOUND==1 || REBOUND==2) 1239 #if (TIME_LAG_MOD==1) 1243 #elif (TIME_LAG_MOD==2) 1245 filename_with_path = trim(inpath)//
'/'//trim(
ch_domain_short)//
'/'// &
1248 open(29, iostat=ios, file=trim(filename_with_path), recl=
rcl1, status=
'old')
1250 if (ios /= 0) stop
' >>> sico_init: Error when opening the time-lag file!' 1252 do n=1, 6;
read(29, fmt=
'(a)') ch_dummy;
end do 1258 close(29, status=
'keep')
1274 #if (FLEX_RIG_MOD==1) 1278 #elif (FLEX_RIG_MOD==2) 1280 filename_with_path = trim(inpath)//
'/'//trim(
ch_domain_short)//
'/'// &
1283 open(29, iostat=ios, file=trim(filename_with_path), recl=
rcl1, status=
'old')
1285 if (ios /= 0) stop
' >>> sico_init: Error when opening the flex-rig file!' 1287 do n=1, 6;
read(29, fmt=
'(a)') ch_dummy;
end do 1293 close(29, status=
'keep')
1297 #elif (REBOUND==0 || REBOUND==1) 1309 stop
' >>> sico_init: ANF_DAT==1 not allowed for scand application!' 1319 call boundary(time_init, dtime, dxi, deta, &
1320 delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1344 stop
' >>> sico_init: Parameter ENHMOD must be between 1 and 5!' 1353 call boundary(time_init, dtime, dxi, deta, &
1354 delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
1385 #if (MARGIN==1 || MARGIN==2) 1394 if ( (
maske(j,i)==0_i2b) &
1396 ( (
maske(j,i+1)==3_i2b) &
1397 .or.(
maske(j,i-1)==3_i2b) &
1398 .or.(
maske(j+1,i)==3_i2b) &
1399 .or.(
maske(j-1,i)==3_i2b) ) &
1403 if ( (
maske(j,i)==3_i2b) &
1405 ( (
maske(j,i+1)==0_i2b) &
1406 .or.(
maske(j,i-1)==0_i2b) &
1407 .or.(
maske(j+1,i)==0_i2b) &
1408 .or.(
maske(j-1,i)==0_i2b) ) &
1412 if ( (
maske(j,i)==3_i2b) &
1414 ( (
maske(j,i+1)==2_i2b) &
1415 .or.(
maske(j,i-1)==2_i2b) &
1416 .or.(
maske(j+1,i)==2_i2b) &
1417 .or.(
maske(j-1,i)==2_i2b) ) &
1421 if ( (
maske(j,i)==2_i2b) &
1423 ( (
maske(j,i+1)==3_i2b) &
1424 .or.(
maske(j,i-1)==3_i2b) &
1425 .or.(
maske(j+1,i)==3_i2b) &
1426 .or.(
maske(j-1,i)==3_i2b) ) &
1436 if ( (
maske(j,i)==2_i2b) &
1437 .and. (
maske(j+1,i)==3_i2b) &
1442 if ( (
maske(j,i)==2_i2b) &
1443 .and. (
maske(j-1,i)==3_i2b) &
1452 if ( (
maske(j,i)==2_i2b) &
1453 .and. (
maske(j,i+1)==3_i2b) &
1458 if ( (
maske(j,i)==2_i2b) &
1459 .and. (
maske(j,i-1)==3_i2b) &
1466 stop
' >>> sico_init: MARGIN must be either 1, 2 or 3!' 1471 #if (GRID==0 || GRID==1) 1475 dist_dxdy(jr,ir) = sqrt( (
real(ir,dp)*dxi)**2 + (
real(jr,dp)*deta)**2 )
1486 #if (DYNAMICS==1 || DYNAMICS==2) 1491 #if (MARGIN==3 || DYNAMICS==2) 1492 call calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t, ii, jj, nn)
1507 stop
' >>> sico_init: DYNAMICS must be either 0, 1 or 2!' 1510 call calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
1514 #if (CALCMOD==0 || CALCMOD==-1) 1539 if ( (
maske(j,i) == 0_i2b).and.(
n_cts(j,i) == 1_i2b) )
then 1552 #elif (CALCMOD==2 || CALCMOD==3) 1570 stop
' >>> sico_init: Parameter CALCMOD must be either -1, 0, 1, 2 or 3!' 1578 filename_with_path = trim(outpath)//
'/'//trim(runname)//
'.ser' 1580 open(12, iostat=ios, file=trim(filename_with_path), status=
'new')
1582 if (ios /= 0) stop
' >>> sico_init: Error when opening the ser file!' 1589 1102
format(
' t(a) D_Ts(deg C) z_sl(m)',/, &
1590 ' V(m^3) V_g(m^3) V_f(m^3)', &
1591 ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1592 ' V_sle(m) V_t(m^3)', &
1594 ' H_max(m) H_t_max(m)', &
1595 ' zs_max(m) vs_max(m/a) Tbh_max(C)')
1596 1103
format(
'----------------------------------------------------', &
1597 '---------------------------------------')
1604 1112
format(
' t(a) glac_ind(1) z_sl(m)',/, &
1605 ' V(m^3) V_g(m^3) V_f(m^3)', &
1606 ' A(m^2) A_g(m^2) A_f(m^2)',/, &
1607 ' V_sle(m) V_t(m^3)', &
1609 ' H_max(m) H_t_max(m)', &
1610 ' zs_max(m) vs_max(m/a) Tbh_max(C)')
1611 1113
format(
'----------------------------------------------------', &
1612 '---------------------------------------')
1620 filename_with_path = trim(outpath)//
'/'//trim(runname)//
'.core' 1622 open(14, iostat=ios, file=trim(filename_with_path), status=
'new')
1624 write(14,
'(1x,a)')
'---------------------' 1625 write(14,
'(1x,a)')
'No boreholes defined.' 1626 write(14,
'(1x,a)')
'---------------------' 1630 #if (defined(OUTPUT_INIT)) 1632 #if (OUTPUT_INIT==0) 1633 flag_init_output = .false.
1634 #elif (OUTPUT_INIT==1) 1635 flag_init_output = .true.
1637 stop
' >>> sico_init: OUTPUT_INIT must be either 0 or 1!' 1641 flag_init_output = .true.
1647 flag_3d_output = .false.
1649 flag_3d_output = .true.
1651 stop
' >>> sico_init: ERGDAT must be either 0 or 1!' 1654 if (flag_init_output) &
1655 call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1656 flag_3d_output, ndat2d, ndat3d)
1660 if (time_output(1) <= time_init+
eps)
then 1663 flag_3d_output = .false.
1665 flag_3d_output = .true.
1667 stop
' >>> sico_init: ERGDAT must be either 0 or 1!' 1670 call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1671 flag_3d_output, ndat2d, ndat3d)
1677 flag_3d_output = .false.
1679 if (flag_init_output) &
1680 call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1681 flag_3d_output, ndat2d, ndat3d)
1683 if (time_output(1) <= time_init+
eps)
then 1685 flag_3d_output = .true.
1687 call output1(runname, time_init, delta_ts, glac_index, z_sl, &
1688 flag_3d_output, ndat2d, ndat3d)
1693 stop
' >>> sico_init: OUTPUT must be either 1, 2 or 3!' 1696 if (flag_init_output)
then 1697 call output2(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1698 call output4(time_init, dxi, deta, delta_ts, glac_index, z_sl)
1713 real(dp),
intent(out) :: dxi, deta
1719 dxi=0.0_dp; deta=0.0_dp
1721 write(6,
'(/1x,a)')
'>>> topography1: ANF_DAT==1 not defined for this domain!' 1733 #if (GRID==0 || GRID==1) 1742 real(dp),
intent(out) :: dxi, deta
1744 integer(i4b) :: i, j, n
1746 real(dp) :: xi0, eta0
1747 character :: ch_dummy
1749 character(len= 8) :: ch_imax
1750 character(len=128) :: fmt4
1751 character(len=256) :: filename_with_path
1753 write(ch_imax, fmt=
'(i8)') imax
1754 write(fmt4, fmt=
'(a)')
'('//trim(adjustl(ch_imax))//
'(i1),i1)' 1758 filename_with_path = trim(inpath)//
'/'//trim(
ch_domain_short)//
'/'// &
1761 open(23, iostat=ios, file=trim(filename_with_path), recl=
rcl1, status=
'old')
1763 if (ios /= 0) stop
' >>> topography2: Error when opening the zl0 file!' 1765 filename_with_path = trim(inpath)//
'/'//trim(
ch_domain_short)//
'/'// &
1766 trim(mask_present_file)
1768 open(24, iostat=ios, file=trim(filename_with_path), recl=
rcl2, status=
'old')
1770 if (ios /= 0) stop
' >>> topography2: Error when opening the mask file!' 1772 do n=1, 6;
read(23, fmt=
'(a)') ch_dummy;
end do 1775 read(23, fmt=*) (
zl0(j,i), i=0,imax)
1778 do n=1, 6;
read(24, fmt=
'(a)') ch_dummy;
end do 1781 read(24, fmt=trim(fmt4)) (
maske(j,i), i=0,imax)
1784 close(23, status=
'keep')
1785 close(24, status=
'keep')
1790 deta = dx *1000.0_dp
1793 eta0 = y0 *1000.0_dp
1798 if (
maske(j,i) <= 1)
then 1805 #if (MARGIN==1 || MARGIN==2) 1815 xi(i) = xi0 +
real(i,
dp)*dxi
1816 eta(j) = eta0 +
real(j,
dp)*deta
1843 #if (GRID==0 || GRID==1) /* Stereographic projection */ 1848 #elif (GRID==2) /* Geographic coordinates */ 1881 subroutine topography3(dxi, deta, z_sl, anfdatname)
1885 #if (GRID==0 || GRID==1) 1894 character(len=100),
intent(in) :: anfdatname
1896 real(dp),
intent(out) :: dxi, deta, z_sl
1898 integer(i4b) :: i, j, n
1900 character(len=256) :: filename_with_path
1901 character :: ch_dummy
1909 filename_with_path = trim(inpath)//
'/'//trim(
ch_domain_short)//
'/'// &
1912 open(23, iostat=ios, file=trim(filename_with_path), recl=
rcl1, status=
'old')
1914 if (ios /= 0) stop
' >>> topography3: Error when opening the zl0 file!' 1916 do n=1, 6;
read(23, fmt=
'(a)') ch_dummy;
end do 1919 read(23, fmt=*) (
zl0(j,i), i=0,imax)
1922 close(23, status=
'keep')
1927 deta = dx *1000.0_dp
1935 #if (GRID==0 || GRID==1) /* Stereographic projection */ 1940 #elif (GRID==2) /* Geographic coordinates */ real(dp), dimension(0:kcmax, 0:jmax, 0:imax) temp_c
temp_c(kc,j,i): Temperature in the upper (kc) ice domain
real(dp), dimension(0:jmax, 0:imax) c_slide
c_slide(j,i): Basal sliding coefficient
real(dp), dimension(0:jmax, 0:imax) phi
phi(j,i): Geographic latitude of grid point (i,j)
subroutine, public read_erg_nc(z_sl, filename)
Reading of time-slice files in native binary or NetCDF format.
subroutine, public read_kei()
Reading of the tabulated kei function.
logical, dimension(0:jmax, 0:imax) flag_inner_point
flag_inner_point(j,i): Inner-point flag. .true.: inner point, .false.: margin point ...
integer(i4b) grip_time_stp
grip_time_stp: Time step of the data values for the surface temperature anomaly
integer(i2b), dimension(0:jmax, 0:imax) n_cts
n_cts(j,i): Mask for thermal conditions. -1: cold ice base, 0: temperate ice base with cold ice above...
integer(i4b) ndata_specmap
ndata_specmap: Number of data values for the sea level
subroutine, public calc_c_int_inv_table()
Computation of the inverse of the temperature integral of the specific heat of ice as a table (c_int_...
real(dp), dimension(0:kcmax) eaz_c
eaz_c(kc): Abbreviation for exp(aa*zeta(kc))
real(dp), dimension(0:jmax, 0:imax) q_geo
q_geo(j,i): Geothermal heat flux
real(dp), dimension(0:jmax, 0:imax, 12) temp_mm_present
temp_mm_present(j,i,n): Present-day mean monthly surface temperature
subroutine sico_init(delta_ts, glac_index, mean_accum, dtime, dtime_temp, dtime_wss, dtime_out, dtime_ser, time, time_init, time_end, time_output, dxi, deta, dzeta_c, dzeta_t, dzeta_r, z_sl, dzsl_dtau, z_mar, ii, jj, nn, ndat2d, ndat3d, n_output, runname)
Main routine of sico_init_m: Initialisations for SICOPOLIS.
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.
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.
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.
integer(i2b) forcing_flag
forcing_flag: Flag for the forcing type. 1: forcing by a spatially constant surface temperature anoma...
logical, dimension(0:jmax, 0:imax) flag_grounding_line_1
flag_grounding_line_1(j,i): Grounding line flag. .true.: grounding line point (grounded ice point wit...
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) enth_t
enth_t(kt,j,i): Enthalpy in the lower (kt) ice domain
Conversion from temperature (temp) and water content (omega) to enthalpy (enth) and vice versa...
subroutine, public init_temp_water_age_2()
Initial temperature, water content and age (case ANF_DAT==2: ice-free conditions with relaxed bedrock...
real(dp), dimension(0:jmax, 0:imax) dh_t_dtau
dH_t_dtau(j,i): Derivative of H_t by tau (time)
real(dp), dimension(0:jmax, 0:imax) lambda
lambda(j,i): Geographic longitude of grid point (i,j)
real(dp), dimension(0:jmax, 0:imax) dh_c_dtau
dH_c_dtau(j,i): Derivative of H_c by tau (time)
integer(i2b), dimension(0:jmax, 0:imax) maske
maske(j,i): Ice-land-ocean mask. 0: grounded ice, 1: ice-free land, 2: ocean, 3: floating ice ...
real(dp), dimension(0:jmax, 0:imax) p_b_w
p_b_w(j,i): Basal water pressure
integer(i4b) specmap_time_min
specmap_time_min: Minimum time of the data values for the sea level
real(dp), dimension(0:jmax, 0:imax) as_perp
as_perp(j,i): Accumulation-ablation function at the ice surface (SMB)
real(dp), dimension(0:jmax, 0:imax) area
area(j,i): Area of grid cell associated with grid point (i,j)
real(dp), dimension(:), allocatable glacial_index
glacial_index(n): Data values for the glacial index
Computation of all components of the strain-rate tensor, the full effective strain rate and the shear...
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...
real(dp), dimension(0:jmax, 0:imax) zs_ref
zs_ref(j,i): Reference elevation for precip_present, temp_ma_present and temp_mj_present ...
subroutine, public calc_enhance_5()
Computation of the flow enhancement factor. Case ENHMOD==5: minimal anisotropic enhancement factor fo...
real(dp) b
B: Semi-minor axis of the planet.
real(dp), dimension(0:jmax, 0:imax) zl
zl(j,i): Coordinate z of the lithosphere surface
subroutine, public calc_vxy_static()
Computation of the horizontal velocity vx, vy, the horizontal volume flux qx, qy etc. for static ice.
subroutine, public output4(time, dxi, deta, delta_ts, glac_index, z_sl)
Writing of time-series data of the deep ice cores on file in ASCII format (and optionally in NetCDF f...
integer(i2b), dimension(0:jmax, 0:imax) kc_cts
kc_cts(j,i): Position kc of the CTS (for COLD, ENTC, ENTM)
integer(i4b) grip_time_max
grip_time_max: Maximum time of the data values for the surface temperature anomaly ...
real(dp), dimension(0:jmax, 0:imax) q_tld
Q_tld(j,i): Water drainage rate from the temperate layer.
subroutine, public calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t, ii, jj, nn)
Computation of the horizontal velocity vx, vy, the horizontal volume flux qx, qy and the flux across ...
Declarations of global variables for SICOPOLIS (for the ANT domain).
integer(i4b), parameter rcl1
rcl1: Maximum length of record for input files (with factor 3 safety margin)
Initialisations for SICOPOLIS.
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.
integer(i4b), dimension(0:jmax, 0:imax) q_weert
q_weert(j,i): Weertman exponent for the basal pressure
subroutine, public calc_temp_melt()
Computation of the melting temperatures.
real(dp) ea
ea: Abbreviation for exp(aa)
real(dp), 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.
real(dp) aa
aa: Exponential stretch parameter of the non-equidistant vertical grid in the upper (kc) ice domain ...
Reading of several input data.
Declarations of kind types for SICOPOLIS.
real(dp), dimension(0:jmax, 0:imax) zs
zs(j,i): Coordinate z of the surface topography
integer(i4b) gi_time_min
gi_time_min: Minimum time of the data values for the glacial index
real(dp), dimension(0:jmax, 0:imax) h_c
H_c(j,i): Thickness of ice in the upper (kc) domain (thickness of the cold-ice layer for POLY...
real(dp) 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...
real(dp), dimension(-190:10) kappa
KAPPA(n): Tabulated values for the heat conductivity of ice.
Computation of the vertical velocity vz.
integer(i4b) gi_time_stp
gi_time_stp: Time step of the data values for the glacial index
real(dp), dimension(0:jmax, 0:imax) dzm_dtau
dzm_dtau(j,i): Derivative of zm by tau (time)
real(dp), dimension(0:krmax) zeta_r
zeta_r(kr): Sigma coordinate zeta_r of grid point kr
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) omega_c
omega_c(kc,j,i): Water content in the upper (kc) ice domain
integer(i2b), dimension(0:jmax, 0:imax) maske_old
maske_old(j,i): Old value of maske (at the previous time step)
character(len=64) ch_domain_long
ch_domain_long: Long name of the computational domain
Definition of the components g11 and g22 of the metric tensor of the applied numerical coordinates...
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.
integer(i4b) gi_time_max
gi_time_max: Maximum time of the data values for the glacial index
real(dp) phi0
PHI0: Standard parallel of the stereographic projection.
subroutine, public ice_mat_eqs_pars(RF_table, R_T_val, KAPPA_table, C_table, n_tmp_min, n_tmp_max, RHO_I_val, RHO_C_val, KAPPA_C_val, C_C_val)
Setting of required physical parameters.
real(dp), dimension(0:jmax, 0:imax, 12) gamma_precip_lgm_anom
gamma_precip_lgm_anom(j,i,n): negative natural logarithm of precip_lgm_anom(j,i,n) ...
Comparison of single- or double-precision floating-point numbers.
subroutine, public calc_enhance_1()
Computation of the flow enhancement factor. Case ENHMOD==1: constant for grounded ice...
real(dp), dimension(0:jmax, 0:imax) zm
zm(j,i): Coordinate z of the bottom of the upper (kc) ice domain = top of the lower (kt) ice domain (...
real(dp), dimension(-190:10) rf
RF(n): Tabulated values for the rate factor of cold ice.
real(dp), dimension(0:jmax, 0:imax) q_b_tot
Q_b_tot(j,i): Sum of Q_bm and Q_tld.
real(dp), dimension(0:jmax, 0:imax) q_bm
Q_bm(j,i): Basal melting rate.
real(dp) function, public enth_fct_temp_omega(temp_val, omega_val)
Enthalpy as a function of temperature and water content.
subroutine, public calc_vxy_sia(dzeta_c, dzeta_t)
Computation of the shear stresses txz, tyz, the effective shear stress sigma, the depth-averaged flui...
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...
integer(i4b) specmap_time_stp
specmap_time_stp: Time step of the data values for the sea level
integer(i4b), parameter rcl2
rcl2: Maximum length of record for input mask files (with factor 3 safety margin) ...
subroutine, public output2(time, dxi, deta, delta_ts, glac_index, z_sl)
Writing of time-series data on file in ASCII format (and optionally in NetCDF format).
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...
real(dp), dimension(0:jmax, 0:imax) sq_g11_g
sq_g11_g(j,i): Square root of the coefficient g11 of the metric tensor on grid point (i...
real(dp), dimension(0:jmax, 0:imax) as_perp_apl
as_perp_apl(j,i): Applied accumulation-ablation function (SMB)
logical, dimension(0:jmax, 0:imax) flag_calving_front_1
flag_calving_front_1(j,i): Calving front flag. .true.: calving front point (floating ice point with a...
subroutine, public metric()
Main routine of module metric_m: Definition of the components g11 and g22 of the metric tensor of the...
subroutine, public calc_enhance_3(time)
Computation of the flow enhancement factor. Case ENHMOD==3: two different values depending on time of...
integer(i4b) grip_time_min
grip_time_min: Minimum time of the data values for the surface temperature anomaly ...
real(dp), dimension(0:jmax, 0:imax) time_lag_asth
time_lag_asth(j,i): Time lag of the relaxing asthenosphere
real(dp), dimension(0:jmax) eta
eta(j): Coordinate eta (= y) of grid point j
real(dp), dimension(0:imax) xi
xi(i): Coordinate xi (= x) of grid point i
character(len=16) ch_domain_short
ch_domain_short: Short name of the computational domain
real(dp) r_t
R_T: Coefficient of the water-content dependence in the rate factor for temperate ice...
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) omega_t
omega_t(kt,j,i): Water content in the lower (kt) ice domain
subroutine, public calc_dxyz(dxi, deta, dzeta_c, dzeta_t)
Main subroutine of calc_dxyz_m: Computation of all components of the strain-rate tensor, the full effective strain rate and the shear fraction.
subroutine, public calc_enhance_4()
Computation of the flow enhancement factor. Case ENHMOD==4: minimal anisotropic enhancement factor fo...
real(dp) kei_r_incr
kei_r_incr: Increment of the argument r of the tabulated kei function
real(dp), dimension(0:ktmax) zeta_t
zeta_t(kt): Sigma coordinate zeta_t of grid point kt
logical, dimension(0:jmax, 0:imax) flag_calving_front_2
flag_calving_front_2(j,i): Calving front flag. .true.: calving front point (ocean point with at least...
Computation of the surface temperature (must be less than 0 deg C!) and of the accumulation-ablation ...
real(dp) l
L: Latent heat of ice.
real(dp), parameter eps_sp_dp
eps_sp_dp: Small number to single-precision accuracy in double precision
integer(i4b) n_data_kei
n_data_kei: Number of tabulated values of the kei function
real(dp) rho
RHO: Density of ice.
subroutine topography3(dxi, deta, z_sl, anfdatname)
Definition of the initial surface and bedrock topography (including gradients) and of the horizontal ...
integer(i4b), dimension(0:jmax, 0:imax) p_weert
p_weert(j,i): Weertman exponent for the basal shear stress
integer(i4b) n_core
n_core: Number of positions to be considered in the time-series file for deep boreholes ...
integer, parameter dp
Double-precision reals.
subroutine topography2(dxi, deta)
Definition of the initial surface and bedrock topography (including gradients) and of the horizontal ...
real(dp), dimension(0:jmax, 0:imax) zb
zb(j,i): Coordinate z of the ice base
logical, dimension(0:jmax, 0:imax) flag_grounding_line_2
flag_grounding_line_2(j,i): Grounding line flag. .true.: grounding line point (floating ice point wit...
real(dp), dimension(0:jmax, 0:imax, 12) precip_present
precip_present(j,i,n): Present-day mean monthly precipitation rate at the ice surface ...
Writing of output data on files.
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.
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...