56 real(dp),
intent(in) :: time, z_sl
59 integer(i4b) :: r_smw_1, s_smw_1, r_smw_2, s_smw_2
60 real(dp),
dimension(0:JMAX,0:IMAX) :: gamma_slide_inv
61 real(dp) :: gamma_slide_inv_1, gamma_slide_inv_2
62 real(dp) :: smw_coeff_1, smw_coeff_2
63 real(dp),
dimension(0:JMAX,0:IMAX) :: p_b, p_b_red, tau_b
64 real(dp) :: cvxy1, cvxy1a, cvxy1b, ctau1, ctau1a, ctau1b
66 real(dp) :: c_Hw_slide, Hw0_slide, Hw0_slide_inv
68 real(dp) :: year_sec_inv, time_in_years
69 real(dp) :: ramp_up_factor
70 logical,
dimension(0:JMAX,0:IMAX) :: sub_melt_flag
72 year_sec_inv = 1.0_dp/year_sec
73 time_in_years = time*year_sec_inv
77 #if (defined(ANT) && (!defined(SEDI_SLIDE) || SEDI_SLIDE==1)) 81 gamma_slide_inv_1 = 1.0_dp/max(gamma_slide,
eps)
88 gamma_slide_inv(j,i) = gamma_slide_inv_1
89 sub_melt_flag(j,i) = (gamma_slide >=
eps)
93 #elif (defined(ANT) && SEDI_SLIDE==2) 97 gamma_slide_inv_1 = 1.0_dp/max(gamma_slide,
eps)
98 gamma_slide_inv_2 = 1.0_dp/max(gamma_slide_sedi,
eps)
105 #if (defined(TRANS_LENGTH_SL)) 107 +r_mask_sedi(j,i) *c_slide_sedi ) &
115 gamma_slide_inv(j,i) = gamma_slide_inv_1
116 sub_melt_flag(j,i) = (gamma_slide >=
eps)
120 #if (defined(TRANS_LENGTH_SL)) 122 +r_mask_sedi(j,i) *c_slide_sedi ) &
125 c_slide(j,i) = c_slide_sedi*year_sec_inv
130 gamma_slide_inv(j,i) = gamma_slide_inv_2
131 sub_melt_flag(j,i) = (gamma_slide_sedi >=
eps)
138 #elif (defined(GRL) && ICE_STREAM==1 \ 139 || defined(asf) && ice_stream==1)
144 gamma_slide_inv_1 = 1.0_dp/max(gamma_slide,
eps)
148 smw_coeff_1 = smw_coeff*year_sec**s_smw_1
153 *(1.0_dp+smw_coeff_1*
runoff(j,i)**s_smw_1 &
157 gamma_slide_inv(j,i) = gamma_slide_inv_1
158 sub_melt_flag(j,i) = (gamma_slide >=
eps)
162 #elif (defined(GRL) && ICE_STREAM==2 \ 163 || defined(asf) && ice_stream==2)
167 gamma_slide_inv_1 = 1.0_dp/max(gamma_slide,
eps)
168 gamma_slide_inv_2 = 1.0_dp/max(gamma_slide_sedi,
eps)
172 smw_coeff_1 = smw_coeff *year_sec**s_smw_1
175 smw_coeff_2 = smw_coeff_sedi*year_sec**s_smw_2
182 *(1.0_dp+smw_coeff_1*
runoff(j,i)**s_smw_1 &
186 gamma_slide_inv(j,i) = gamma_slide_inv_1
187 sub_melt_flag(j,i) = (gamma_slide >=
eps)
189 c_slide(j,i) = c_slide_sedi*year_sec_inv &
190 *(1.0_dp+smw_coeff_2*
runoff(j,i)**s_smw_2 &
194 gamma_slide_inv(j,i) = gamma_slide_inv_2
195 sub_melt_flag(j,i) = (gamma_slide_sedi >=
eps)
201 #elif (defined(XYZ) && defined(HEINO)) 205 gamma_slide_inv_1 = 1.0_dp/max(gamma_slide,
eps)
206 gamma_slide_inv_2 = 1.0_dp/max(gamma_slide_sedi,
eps)
215 gamma_slide_inv(j,i) = gamma_slide_inv_1
216 sub_melt_flag(j,i) = (gamma_slide >=
eps)
218 c_slide(j,i) = c_slide_sedi*year_sec_inv
221 gamma_slide_inv(j,i) = gamma_slide_inv_2
222 sub_melt_flag(j,i) = (gamma_slide_sedi >=
eps)
228 #elif (defined(SCAND) \ 233 || defined(emtp2sge) \
237 gamma_slide_inv_1 = 1.0_dp/max(gamma_slide,
eps)
244 gamma_slide_inv(j,i) = gamma_slide_inv_1
245 sub_melt_flag(j,i) = (gamma_slide >=
eps)
251 stop
' >>> calc_vxy_b_sia: No valid domain specified!' 257 p_weert_inv(j,i) = 1.0_dp/max(
real(p_weert(j,i),dp), eps)
263 ramp_up_factor = 1.0_dp
265 #if (defined(TIME_RAMP_UP_SLIDE)) 267 if ( (time_ramp_up_slide >
eps_dp) &
269 ((time_in_years-(time_init0)) < (time_ramp_up_slide)) )
then 271 ramp_up_factor = (time_in_years-(time_init0))/(time_ramp_up_slide)
273 ramp_up_factor = max(min(ramp_up_factor, 1.0_dp), 0.0_dp)
276 ramp_up_factor = ramp_up_factor*ramp_up_factor*ramp_up_factor &
277 *(10.0_dp + ramp_up_factor &
278 *(-15.0_dp+6.0_dp*ramp_up_factor))
289 #if (BASAL_HYDROLOGY==1) 291 #if (defined(C_HW_SLIDE)) 292 c_hw_slide = c_hw_slide
297 #if (defined(HW0_SLIDE)) 298 hw0_slide = hw0_slide
303 hw0_slide_inv = 1.0_dp/max(hw0_slide,
eps_dp)
308 hw0_slide_inv = 1.0_dp
324 p_b_red(j,i) = p_b(j,i)-
p_b_w(j,i)
325 p_b_red(j,i) = max(p_b_red(j,i), red_pres_limit_fact*p_b(j,i))
333 p_b_red(j,i) = 0.0_dp
376 if (
n_cts(j,i) == -1)
then 378 if (sub_melt_flag(j,i))
then 380 cvxy1a = exp(-gamma_slide_inv(j,i)*temp_diff)
388 c_drag(j,i) = ctau1*ctau1a
390 else if (
n_cts(j,i) == 0)
then 404 #if (BASAL_HYDROLOGY==1) 407 + c_hw_slide*(1.0_dp-exp(-max(
h_w(j,i),0.0_dp)*hw0_slide_inv))
454 vh_max = vh_max/year_sec
483 real(dp),
intent(in) :: dzeta_c, dzeta_t
485 integer(i4b) :: i, j, kc, kt
486 real(dp) :: avxy3(0:kcmax), aqxy1(0:kcmax)
487 real(dp) :: ctxyz1(0:kcmax,0:jmax,0:imax), &
488 ctxyz2(0:ktmax,0:jmax,0:imax)
489 real(dp) :: flui_t(0:ktmax), flui_c(0:kcmax)
490 real(dp) :: cflui0(0:ktmax), cflui1(0:kcmax)
491 real(dp) :: cvxy2(0:ktmax), cvxy3(0:kcmax)
492 real(dp) :: cqxy0(0:ktmax), cqxy1(0:kcmax)
494 real(dp) :: ratio_sl_threshold
500 avxy3(kc) =
aa*
eaz_c(kc)/(
ea-1.0_dp)*dzeta_c
501 aqxy1(kc) =
aa/(
ea-1.0_dp)*
eaz_c(kc)*dzeta_c
522 if (
n_cts(j,i) == 1)
then 531 ctxyz2(kt,j,i) = 0.0_dp
539 ctxyz1(kc,j,i) = 0.0_dp
543 ctxyz2(kt,j,i) = 0.0_dp
557 txz_c(kc,j,i) = -0.5_dp*(ctxyz1(kc,j,i)+ctxyz1(kc,j,i+1)) &
563 -0.5_dp*(ctxyz2(kt,j,i)+ctxyz2(kt,j,i+1)) &
576 tyz_c(kc,j,i) = -0.5_dp*(ctxyz1(kc,j,i)+ctxyz1(kc,j+1,i)) &
582 -0.5_dp*(ctxyz2(kt,j,i)+ctxyz2(kt,j+1,i)) &
595 sigma_c(kc,j,i) = ctxyz1(kc,j,i) &
620 flui_t(kt) = 2.0_dp &
624 cflui0(kt) =
h_t(j,i)*flui_t(kt)*dzeta_t
628 flui_c(kc) = 2.0_dp &
630 #if (CALCMOD==0 || CALCMOD==1 || CALCMOD==-1) 632 #elif (CALCMOD==2 || CALCMOD==3) 637 cflui1(kc) = aqxy1(kc)*
h_c(j,i)*flui_c(kc)
644 if (
n_cts(j,i) == 1)
then 679 cvxy2(kt) = 2.0_dp*
h_t(j,i) &
683 *(ctxyz1(0,j,i)+ctxyz2(kt,j,i)) &
688 cvxy3(kc) = 2.0_dp*avxy3(kc)*
h_c(j,i) &
690 #if (CALCMOD==0 || CALCMOD==1 || CALCMOD==-1) 692 #elif (CALCMOD==2 || CALCMOD==3) 702 if (
n_cts(j,i) == -1)
then 712 +0.5_dp*(cvxy3(kc+1)+cvxy3(kc))
715 else if (
n_cts(j,i) == 0)
then 725 +0.5_dp*(cvxy3(kc+1)+cvxy3(kc))
734 +0.5_dp*(cvxy2(kt+1)+cvxy2(kt))
741 +0.5_dp*(cvxy3(kc+1)+cvxy3(kc))
810 vh_max = vh_max/year_sec
819 vx_t(kt,j,i) = max(
vx_t(kt,j,i), -vh_max)
820 vx_t(kt,j,i) = min(
vx_t(kt,j,i), vh_max)
821 vy_t(kt,j,i) = max(
vy_t(kt,j,i), -vh_max)
822 vy_t(kt,j,i) = min(
vy_t(kt,j,i), vh_max)
825 vx_c(kc,j,i) = max(
vx_c(kc,j,i), -vh_max)
826 vx_c(kc,j,i) = min(
vx_c(kc,j,i), vh_max)
827 vy_c(kc,j,i) = max(
vy_c(kc,j,i), -vh_max)
828 vy_c(kc,j,i) = min(
vy_c(kc,j,i), vh_max)
856 if (
n_cts(j,i) == 1)
then 859 h_diff(j,i) =
h_diff(j,i)+0.5_dp*(cqxy0(kt+1)+cqxy0(kt))
865 h_diff(j,i) =
h_diff(j,i)+0.5_dp*(cqxy1(kc+1)+cqxy1(kc))
890 if ( (
maske(j,i)==0).or.(
maske(j,i+1)==0) )
then 916 if ( (
maske(j,i)==0).or.(
maske(j+1,i)==0) )
then 943 #if (DYNAMICS==0 || DYNAMICS==1) 945 ratio_sl_threshold = 1.11e+11_dp
949 #if ( defined(RATIO_SL_THRESH) ) 950 ratio_sl_threshold = ratio_sl_thresh
952 ratio_sl_threshold = 0.5_dp
970 if ( (
maske(j,i) == 0) &
986 stop
' >>> calc_vxy_sia: DYNAMICS must be 0, 1 or 2!' 1068 subroutine calc_vxy_ssa(z_sl, dxi, deta, dzeta_c, dzeta_t, ii, jj, nn)
1072 integer(i4b),
intent(in) :: ii((imax+1)*(jmax+1)), jj((imax+1)*(jmax+1))
1073 integer(i4b),
intent(in) :: nn(0:jmax,0:imax)
1074 real(dp),
intent(in) :: z_sl, dxi, deta, dzeta_c, dzeta_t
1076 integer(i4b) :: i, j, kc, kt, m
1077 integer(i4b) :: iter_ssa
1078 real(dp),
dimension(0:JMAX,0:IMAX) :: vx_m_sia, vy_m_sia
1079 real(dp),
dimension(0:JMAX,0:IMAX) :: vx_m_prev, vy_m_prev
1081 real(dp) :: dxi_inv, deta_inv
1083 real(dp) :: ratio_sl_threshold, ratio_help
1084 real(dp),
dimension(0:JMAX,0:IMAX) :: weigh_ssta_sia_x, weigh_ssta_sia_y
1085 real(dp) :: qx_gl_g, qy_gl_g
1087 #if (MARGIN==3 || DYNAMICS==2) 1091 #if (defined(N_ITER_SSA)) 1092 iter_ssa = max(n_iter_ssa, 1)
1097 #if (defined(RELAX_FACT_SSA)) 1098 rel_ssa = relax_fact_ssa
1103 write(6,
'(10x,a,i0,a,f6.3)')
'calc_vxy_ssa: iter_ssa = ', iter_ssa, &
1104 ', rel_ssa =' , rel_ssa
1117 if ( (
maske(j,i)==0_i2b) &
1119 ( (
maske(j,i+1)==3_i2b) &
1120 .or.(
maske(j,i-1)==3_i2b) &
1121 .or.(
maske(j+1,i)==3_i2b) &
1122 .or.(
maske(j-1,i)==3_i2b) ) &
1126 if ( (
maske(j,i)==3_i2b) &
1128 ( (
maske(j,i+1)==0_i2b) &
1129 .or.(
maske(j,i-1)==0_i2b) &
1130 .or.(
maske(j+1,i)==0_i2b) &
1131 .or.(
maske(j-1,i)==0_i2b) ) &
1135 if ( (
maske(j,i)==3_i2b) &
1137 ( (
maske(j,i+1)==2_i2b) &
1138 .or.(
maske(j,i-1)==2_i2b) &
1139 .or.(
maske(j+1,i)==2_i2b) &
1140 .or.(
maske(j-1,i)==2_i2b) ) &
1144 if ( (
maske(j,i)==2_i2b) &
1146 ( (
maske(j,i+1)==3_i2b) &
1147 .or.(
maske(j,i-1)==3_i2b) &
1148 .or.(
maske(j+1,i)==3_i2b) &
1149 .or.(
maske(j-1,i)==3_i2b) ) &
1159 if ( (
maske(j,i)==2_i2b) &
1160 .and. (
maske(j+1,i)==3_i2b) &
1165 if ( (
maske(j,i)==2_i2b) &
1166 .and. (
maske(j-1,i)==3_i2b) &
1175 if ( (
maske(j,i)==2_i2b) &
1176 .and. (
maske(j,i+1)==3_i2b) &
1181 if ( (
maske(j,i)==2_i2b) &
1182 .and. (
maske(j,i-1)==3_i2b) &
1204 call calc_vxy_ssa_matrix(z_sl, dxi, deta, ii, jj, nn, m, iter_ssa)
1217 call calc_vis_ssa(dxi, deta, dzeta_c, dzeta_t)
1221 call calc_vxy_ssa_matrix(z_sl, dxi, deta, ii, jj, nn, m, iter_ssa)
1225 vx_m = rel_ssa*
vx_m + (1.0_dp-rel_ssa)*vx_m_prev
1226 vy_m = rel_ssa*
vy_m + (1.0_dp-rel_ssa)*vy_m_prev
1232 call calc_vis_ssa(dxi, deta, dzeta_c, dzeta_t)
1237 vh_max = vh_max/year_sec
1246 #if (DYNAMICS==0 || DYNAMICS==1) 1248 ratio_sl_threshold = 1.11e+11_dp
1253 #if ( defined(RATIO_SL_THRESH) ) 1254 ratio_sl_threshold = ratio_sl_thresh
1256 ratio_sl_threshold = 0.5_dp
1259 ratio_help = 1.0_dp/(1.0_dp-ratio_sl_threshold)
1262 stop
' >>> calc_vxy_ssa: DYNAMICS must be 0, 1 or 2!' 1265 weigh_ssta_sia_x = 0.0_dp
1266 weigh_ssta_sia_y = 0.0_dp
1276 weigh_ssta_sia_x(j,i) = (
ratio_sl_x(j,i)-ratio_sl_threshold)*ratio_help
1278 weigh_ssta_sia_x(j,i) = max(min(weigh_ssta_sia_x(j,i), 1.0_dp), 0.0_dp)
1281 #if (SSTA_SIA_WEIGH_FCT==0) 1285 #elif (SSTA_SIA_WEIGH_FCT==1) 1287 weigh_ssta_sia_x(j,i) = weigh_ssta_sia_x(j,i)*weigh_ssta_sia_x(j,i) &
1288 *(3.0_dp-2.0_dp*weigh_ssta_sia_x(j,i))
1291 #elif (SSTA_SIA_WEIGH_FCT==2) 1293 weigh_ssta_sia_x(j,i) = weigh_ssta_sia_x(j,i)*weigh_ssta_sia_x(j,i) &
1294 *weigh_ssta_sia_x(j,i) &
1295 *(10.0_dp + weigh_ssta_sia_x(j,i) &
1296 *(-15.0_dp+6.0_dp*weigh_ssta_sia_x(j,i)))
1300 stop
' >>> calc_vxy_ssa: SSTA_SIA_WEIGH_FCT must be 0, 1 or 2!' 1304 vx_t(kt,j,i) = weigh_ssta_sia_x(j,i)*
vx_m(j,i) &
1305 + (1.0_dp-weigh_ssta_sia_x(j,i))*
vx_t(kt,j,i)
1309 vx_c(kc,j,i) = weigh_ssta_sia_x(j,i)*
vx_m(j,i) &
1310 + (1.0_dp-weigh_ssta_sia_x(j,i))*
vx_c(kc,j,i)
1315 vx_m(j,i) = weigh_ssta_sia_x(j,i)*
vx_m(j,i) &
1316 + (1.0_dp-weigh_ssta_sia_x(j,i))*vx_m_sia(j,i)
1320 else if ( (
maske(j,i)==3_i2b).or.(
maske(j,i+1)==3_i2b) )
then 1347 weigh_ssta_sia_y(j,i) = (
ratio_sl_y(j,i)-ratio_sl_threshold)*ratio_help
1349 weigh_ssta_sia_y(j,i) = max(min(weigh_ssta_sia_y(j,i), 1.0_dp), 0.0_dp)
1352 #if (SSTA_SIA_WEIGH_FCT==0) 1356 #elif (SSTA_SIA_WEIGH_FCT==1) 1358 weigh_ssta_sia_y(j,i) = weigh_ssta_sia_y(j,i)*weigh_ssta_sia_y(j,i) &
1359 *(3.0_dp-2.0_dp*weigh_ssta_sia_y(j,i))
1362 #elif (SSTA_SIA_WEIGH_FCT==2) 1364 weigh_ssta_sia_y(j,i) = weigh_ssta_sia_y(j,i)*weigh_ssta_sia_y(j,i) &
1365 *weigh_ssta_sia_y(j,i) &
1366 *(10.0_dp + weigh_ssta_sia_y(j,i) &
1367 *(-15.0_dp+6.0_dp*weigh_ssta_sia_y(j,i)))
1371 stop
' >>> calc_vxy_ssa: SSTA_SIA_WEIGH_FCT must be 0, 1 or 2!' 1375 vy_t(kt,j,i) = weigh_ssta_sia_y(j,i)*
vy_m(j,i) &
1376 + (1.0_dp-weigh_ssta_sia_y(j,i))*
vy_t(kt,j,i)
1380 vy_c(kc,j,i) = weigh_ssta_sia_y(j,i)*
vy_m(j,i) &
1381 + (1.0_dp-weigh_ssta_sia_y(j,i))*
vy_c(kc,j,i)
1386 vy_m(j,i) = weigh_ssta_sia_y(j,i)*
vy_m(j,i) &
1387 + (1.0_dp-weigh_ssta_sia_y(j,i))*vy_m_sia(j,i)
1391 else if ( (
maske(j,i)==3_i2b).or.(
maske(j+1,i)==3_i2b) )
then 1424 else if (
maske(j,i)==3_i2b)
then 1444 qx_gl_g = 0.5_dp*(
qx(j,i)+
qx(j,i-1))
1445 qy_gl_g = 0.5_dp*(
qy(j,i)+
qy(j-1,i))
1447 q_gl_g(j,i) = sqrt(qx_gl_g*qx_gl_g+qy_gl_g*qy_gl_g)
1456 stop
' >>> calc_vxy_ssa: Only to be called for MARGIN==3 or DYNAMICS==2!' 1466 subroutine calc_vxy_ssa_matrix(z_sl,dxi, deta, ii, jj, nn, m, iter_ssa)
1470 integer(i4b),
intent(in) :: ii((imax+1)*(jmax+1)), jj((imax+1)*(jmax+1))
1471 integer(i4b),
intent(in) :: nn(0:jmax,0:imax)
1472 integer(i4b),
intent(in) :: m, iter_ssa
1473 real(dp),
intent(in) :: dxi, deta, z_sl
1475 integer(i4b) :: i, j, k, n
1476 integer(i4b) :: i1, j1
1477 real(dp) :: inv_dxi, inv_deta, inv_dxi_deta, inv_dxi2, inv_deta2
1478 real(dp) :: factor_rhs_1, factor_rhs_2, rhosw_rho_ratio
1479 real(dp) :: H_mid, zl_mid
1480 real(dp),
dimension(0:JMAX,0:IMAX) :: vis_int_sgxy, beta_drag
1481 character(len=256) :: ch_solver_set_option
1483 #if (MARGIN==3 || DYNAMICS==2) 1486 lis_integer :: nc, nr
1489 lis_vector :: lgs_b, lgs_x
1490 lis_solver :: solver
1492 lis_integer,
parameter :: nmax = 2*(imax+1)*(jmax+1)
1493 lis_integer,
parameter :: n_sprs = 20*(imax+1)*(jmax+1)
1494 lis_integer,
allocatable,
dimension(:) :: lgs_a_ptr, lgs_a_index
1495 lis_scalar,
allocatable,
dimension(:) :: lgs_a_value, lgs_b_value, lgs_x_value
1499 inv_dxi = 1.0_dp/dxi
1500 inv_deta = 1.0_dp/deta
1501 inv_dxi_deta = 1.0_dp/(dxi*deta)
1502 inv_dxi2 = 1.0_dp/(dxi*dxi)
1503 inv_deta2 = 1.0_dp/(deta*deta)
1506 factor_rhs_1 = 0.5_dp*
rho*
g 1512 vis_int_sgxy = 0.0_dp
1519 if ((
maske(j,i)==0_i2b).or.(
maske(j,i)==3_i2b))
then 1521 vis_int_sgxy(j,i) = vis_int_sgxy(j,i) +
vis_int_g(j,i)
1524 if ((
maske(j,i+1)==0_i2b).or.(
maske(j,i+1)==3_i2b))
then 1526 vis_int_sgxy(j,i) = vis_int_sgxy(j,i) +
vis_int_g(j,i+1)
1529 if ((
maske(j+1,i)==0_i2b).or.(
maske(j+1,i)==3_i2b))
then 1531 vis_int_sgxy(j,i) = vis_int_sgxy(j,i) +
vis_int_g(j+1,i)
1534 if ((
maske(j+1,i+1)==0_i2b).or.(
maske(j+1,i+1)==3_i2b))
then 1536 vis_int_sgxy(j,i) = vis_int_sgxy(j,i) +
vis_int_g(j+1,i+1)
1539 if (k>0) vis_int_sgxy(j,i) = vis_int_sgxy(j,i)/
real(k,
dp)
1555 beta_drag(j,i) =
c_drag(j,i) &
1562 beta_drag(j,i) =
c_drag(j,i) &
1563 / sqrt((0.5_dp*(
vx_m(j,i)+
vx_m(j,i-1)))**2 &
1564 +(0.5_dp*(
vy_m(j,i)+
vy_m(j-1,i)))**2) &
1578 allocate(lgs_a_value(n_sprs), lgs_a_index(n_sprs), lgs_a_ptr(nmax+1))
1579 allocate(lgs_b_value(nmax), lgs_x_value(nmax))
1581 lgs_a_value = 0.0_dp
1585 lgs_b_value = 0.0_dp
1586 lgs_x_value = 0.0_dp
1601 if ( (i /= imax).and.(j /= 0).and.(j /= jmax) )
then 1604 h_mid = 0.50_dp*(
h_c(j,i)+
h_t(j,i)+
h_c(j,i+1)+
h_t(j,i+1))
1605 zl_mid = 0.50_dp*(
zl(j,i)+
zl(j,i+1))
1608 ( (
maske(j,i)==3_i2b).and.(
maske(j,i+1)==3_i2b) ) &
1611 .and.(h_mid < (z_sl-zl_mid)*rhosw_rho_ratio) ) &
1614 .and.(h_mid < (z_sl-zl_mid)*rhosw_rho_ratio) ) &
1628 lgs_a_value(k) = inv_deta2*vis_int_sgxy(j-1,i)
1634 lgs_a_value(k) = inv_dxi_deta &
1635 *(2.0_dp*
vis_int_g(j,i)+vis_int_sgxy(j-1,i))
1641 lgs_a_value(k) = 4.0_dp*inv_dxi2*
vis_int_g(j,i)
1647 lgs_a_value(k) = -inv_dxi_deta &
1648 *(2.0_dp*
vis_int_g(j,i+1)+vis_int_sgxy(j-1,i))
1653 stop
' >>> calc_vxy_ssa_matrix: Check for diagonal element failed!' 1656 lgs_a_value(k) = -4.0_dp*inv_dxi2 &
1659 *(vis_int_sgxy(j,i)+vis_int_sgxy(j-1,i))
1665 lgs_a_value(k) = -inv_dxi_deta &
1666 *(2.0_dp*
vis_int_g(j,i)+vis_int_sgxy(j,i))
1672 lgs_a_value(k) = 4.0_dp*inv_dxi2*
vis_int_g(j,i+1)
1678 lgs_a_value(k) = inv_dxi_deta &
1679 *(2.0_dp*
vis_int_g(j,i+1)+vis_int_sgxy(j,i))
1685 lgs_a_value(k) = inv_deta2*vis_int_sgxy(j,i)
1688 lgs_b_value(nr) = factor_rhs_1 &
1692 lgs_x_value(nr) =
vx_m(j,i)
1700 lgs_a_value(k) = inv_deta2*vis_int_sgxy(j-1,i)
1706 lgs_a_value(k) = inv_dxi_deta &
1707 *(2.0_dp*
vis_int_g(j,i)+vis_int_sgxy(j-1,i))
1713 lgs_a_value(k) = 4.0_dp*inv_dxi2*
vis_int_g(j,i)
1719 lgs_a_value(k) = -inv_dxi_deta &
1720 *(2.0_dp*
vis_int_g(j,i+1)+vis_int_sgxy(j-1,i))
1725 stop
' >>> calc_vxy_ssa_matrix: Check for diagonal element failed!' 1728 lgs_a_value(k) = -4.0_dp*inv_dxi2 &
1731 *(vis_int_sgxy(j,i)+vis_int_sgxy(j-1,i)) &
1732 -0.5_dp*(beta_drag(j,i+1)+beta_drag(j,i))
1738 lgs_a_value(k) = -inv_dxi_deta &
1739 *(2.0_dp*
vis_int_g(j,i)+vis_int_sgxy(j,i))
1745 lgs_a_value(k) = 4.0_dp*inv_dxi2*
vis_int_g(j,i+1)
1751 lgs_a_value(k) = inv_dxi_deta &
1752 *(2.0_dp*
vis_int_g(j,i+1)+vis_int_sgxy(j,i))
1758 lgs_a_value(k) = inv_deta2*vis_int_sgxy(j,i)
1761 lgs_b_value(nr) = factor_rhs_1 &
1765 lgs_x_value(nr) =
vx_m(j,i)
1769 .and.(h_mid >= (z_sl-zl_mid)*rhosw_rho_ratio) ) &
1772 .and.(h_mid >= (z_sl-zl_mid)*rhosw_rho_ratio) ) &
1780 lgs_a_value(k) = 1.0_dp
1783 lgs_b_value(nr) =
vx_m(j,i)
1784 lgs_x_value(nr) =
vx_m(j,i)
1787 ( (
maske(j,i)==3_i2b).and.(
maske(j,i+1)==1_i2b) ) &
1789 ( (
maske(j,i)==1_i2b).and.(
maske(j,i+1)==3_i2b) ) &
1796 lgs_a_value(k) = 1.0_dp
1799 lgs_b_value(nr) = 0.0_dp
1801 lgs_x_value(nr) = 0.0_dp
1811 if (
maske(j,i)==3_i2b)
then 1818 ( (
maske(j,i1-1)==3_i2b).or.(
maske(j,i1+1)==3_i2b) ) &
1820 ( (
maske(j,i1-1)==0_i2b).or.(
maske(j,i1+1)==0_i2b) ) &
1822 ( (
maske(j,i1-1)==1_i2b).or.(
maske(j,i1+1)==1_i2b) ) &
1829 lgs_a_value(k) = -2.0_dp*inv_deta*
vis_int_g(j,i1)
1835 lgs_a_value(k) = -4.0_dp*inv_dxi*
vis_int_g(j,i1)
1841 lgs_a_value(k) = 4.0_dp*inv_dxi*
vis_int_g(j,i1)
1847 lgs_a_value(k) = 2.0_dp*inv_deta*
vis_int_g(j,i1)
1850 lgs_b_value(nr) = factor_rhs_2 &
1853 lgs_x_value(nr) =
vx_m(j,i)
1860 lgs_a_value(k) = 1.0_dp
1863 lgs_b_value(nr) = 0.0_dp
1865 lgs_x_value(nr) = 0.0_dp
1869 else if ( (
maske(j,i)==0_i2b).or.(
maske(j,i+1)==0_i2b) )
then 1875 lgs_a_value(k) = 1.0_dp
1878 lgs_b_value(nr) =
vx_m(j,i)
1880 lgs_x_value(nr) =
vx_m(j,i)
1887 lgs_a_value(k) = 1.0_dp
1890 lgs_b_value(nr) = 0.0_dp
1892 lgs_x_value(nr) = 0.0_dp
1900 lgs_a_value(k) = 1.0_dp
1903 lgs_b_value(nr) = 0.0_dp
1905 lgs_x_value(nr) = 0.0_dp
1909 lgs_a_ptr(nr+1) = k+1
1915 if ( (j /= jmax).and.(i /= 0).and.(i /= imax) )
then 1918 h_mid = 0.5_dp*(
h_c(j+1,i)+
h_t(j+1,i)+
h_c(j,i)+
h_t(j,i))
1919 zl_mid = 0.5_dp*(
zl(j+1,i)+
zl(j,i))
1922 ( (
maske(j,i)==3_i2b).and.(
maske(j+1,i)==3_i2b) ) &
1925 .and.(h_mid < (z_sl-zl_mid)*rhosw_rho_ratio) ) &
1928 .and.(h_mid < (z_sl-zl_mid)*rhosw_rho_ratio) ) &
1942 lgs_a_value(k) = 4.0_dp*inv_deta2*
vis_int_g(j,i)
1948 lgs_a_value(k) = inv_dxi_deta &
1949 *(2.0_dp*
vis_int_g(j,i)+vis_int_sgxy(j,i-1))
1955 lgs_a_value(k) = inv_dxi2*vis_int_sgxy(j,i-1)
1961 lgs_a_value(k) = -inv_dxi_deta &
1962 *(2.0_dp*
vis_int_g(j,i)+vis_int_sgxy(j,i))
1967 stop
' >>> calc_vxy_ssa_matrix: Check for diagonal element failed!' 1970 lgs_a_value(k) = -4.0_dp*inv_deta2 &
1973 *(vis_int_sgxy(j,i)+vis_int_sgxy(j,i-1))
1976 nc = 2*nn(j+1,i-1)-1
1979 lgs_a_value(k) = -inv_dxi_deta &
1980 *(2.0_dp*
vis_int_g(j+1,i)+vis_int_sgxy(j,i-1))
1986 lgs_a_value(k) = inv_dxi2*vis_int_sgxy(j,i)
1992 lgs_a_value(k) = inv_dxi_deta &
1993 *(2.0_dp*
vis_int_g(j+1,i)+vis_int_sgxy(j,i))
1999 lgs_a_value(k) = 4.0_dp*inv_deta2*
vis_int_g(j+1,i)
2002 lgs_b_value(nr) = factor_rhs_1 &
2006 lgs_x_value(nr) =
vy_m(j,i)
2014 lgs_a_value(k) = 4.0_dp*inv_deta2*
vis_int_g(j,i)
2020 lgs_a_value(k) = inv_dxi_deta &
2021 *(2.0_dp*
vis_int_g(j,i)+vis_int_sgxy(j,i-1))
2027 lgs_a_value(k) = inv_dxi2*vis_int_sgxy(j,i-1)
2033 lgs_a_value(k) = -inv_dxi_deta &
2034 *(2.0_dp*
vis_int_g(j,i)+vis_int_sgxy(j,i))
2039 stop
' >>> calc_vxy_ssa_matrix: Check for diagonal element failed!' 2042 lgs_a_value(k) = -4.0_dp*inv_deta2 &
2045 *(vis_int_sgxy(j,i)+vis_int_sgxy(j,i-1)) &
2046 -0.5_dp*(beta_drag(j+1,i)+beta_drag(j,i))
2049 nc = 2*nn(j+1,i-1)-1
2052 lgs_a_value(k) = -inv_dxi_deta &
2053 *(2.0_dp*
vis_int_g(j+1,i)+vis_int_sgxy(j,i-1))
2059 lgs_a_value(k) = inv_dxi2*vis_int_sgxy(j,i)
2065 lgs_a_value(k) = inv_dxi_deta &
2066 *(2.0_dp*
vis_int_g(j+1,i)+vis_int_sgxy(j,i))
2072 lgs_a_value(k) = 4.0_dp*inv_deta2*
vis_int_g(j+1,i)
2075 lgs_b_value(nr) = factor_rhs_1 &
2079 lgs_x_value(nr) =
vy_m(j,i)
2083 .and.(h_mid >= (z_sl-zl_mid)*rhosw_rho_ratio) ) &
2086 .and.(h_mid >= (z_sl-zl_mid)*rhosw_rho_ratio) ) &
2094 lgs_a_value(k) = 1.0_dp
2097 lgs_b_value(nr) =
vy_m(j,i)
2098 lgs_x_value(nr) =
vy_m(j,i)
2101 ( (
maske(j,i)==3_i2b).and.(
maske(j+1,i)==1_i2b) ) &
2103 ( (
maske(j,i)==1_i2b).and.(
maske(j+1,i)==3_i2b) ) &
2110 lgs_a_value(k) = 1.0_dp
2113 lgs_b_value(nr) = 0.0_dp
2115 lgs_x_value(nr) = 0.0_dp
2125 if (
maske(j,i)==3_i2b)
then 2132 ( (
maske(j1-1,i)==3_i2b).or.(
maske(j1+1,i)==3_i2b) ) &
2134 ( (
maske(j1-1,i)==0_i2b).or.(
maske(j1+1,i)==0_i2b) ) &
2136 ( (
maske(j1-1,i)==1_i2b).or.(
maske(j1+1,i)==1_i2b) ) &
2143 lgs_a_value(k) = -4.0_dp*inv_deta*
vis_int_g(j1,i)
2149 lgs_a_value(k) = -2.0_dp*inv_dxi*
vis_int_g(j1,i)
2155 lgs_a_value(k) = 2.0_dp*inv_dxi*
vis_int_g(j1,i)
2161 lgs_a_value(k) = 4.0_dp*inv_deta*
vis_int_g(j1,i)
2164 lgs_b_value(nr) = factor_rhs_2 &
2167 lgs_x_value(nr) =
vy_m(j,i)
2174 lgs_a_value(k) = 1.0_dp
2177 lgs_b_value(nr) = 0.0_dp
2179 lgs_x_value(nr) = 0.0_dp
2183 else if ( (
maske(j,i)==0_i2b).or.(
maske(j+1,i)==0_i2b) )
then 2189 lgs_a_value(k) = 1.0_dp
2192 lgs_b_value(nr) =
vy_m(j,i)
2194 lgs_x_value(nr) =
vy_m(j,i)
2201 lgs_a_value(k) = 1.0_dp
2204 lgs_b_value(nr) = 0.0_dp
2205 lgs_x_value(nr) = 0.0_dp
2213 lgs_a_value(k) = 1.0_dp
2216 lgs_b_value(nr) = 0.0_dp
2217 lgs_x_value(nr) = 0.0_dp
2221 lgs_a_ptr(nr+1) = k+1
2227 call lis_matrix_create(lis_comm_world, lgs_a, ierr)
2228 call lis_vector_create(lis_comm_world, lgs_b, ierr)
2229 call lis_vector_create(lis_comm_world, lgs_x, ierr)
2231 call lis_matrix_set_size(lgs_a, 0, nmax, ierr)
2232 call lis_vector_set_size(lgs_b, 0, nmax, ierr)
2233 call lis_vector_set_size(lgs_x, 0, nmax, ierr)
2237 do nc=lgs_a_ptr(nr), lgs_a_ptr(nr+1)-1
2238 call lis_matrix_set_value(lis_ins_value, nr, lgs_a_index(nc), &
2239 lgs_a_value(nc), lgs_a, ierr)
2242 call lis_vector_set_value(lis_ins_value, nr, lgs_b_value(nr), lgs_b, ierr)
2243 call lis_vector_set_value(lis_ins_value, nr, lgs_x_value(nr), lgs_x, ierr)
2247 call lis_matrix_set_type(lgs_a, lis_matrix_csr, ierr)
2248 call lis_matrix_assemble(lgs_a, ierr)
2252 call lis_solver_create(solver, ierr)
2255 #if (defined(LIS_OPTS)) 2256 ch_solver_set_option = trim(lis_opts)
2258 ch_solver_set_option =
'-i bicgsafe -p jacobi '// &
2259 '-maxiter 100 -tol 1.0e-4 -initx_zeros false' 2262 call lis_solver_set_option(trim(ch_solver_set_option), solver, ierr)
2265 call lis_solve(lgs_a, lgs_b, lgs_x, solver, ierr)
2268 call lis_solver_get_iter(solver, iter, ierr)
2270 if (iter_ssa == 1)
then 2271 write(6,
'(10x,a,i0)')
'calc_vxy_ssa_matrix: LIS iter = ', iter
2272 else if (m == 1)
then 2273 write(6,
'(10x,a,i0)', advance=
'no')
'calc_vxy_ssa_matrix: LIS iter = ', iter
2274 else if (m == iter_ssa)
then 2275 write(6,
'(a,i0)')
', ', iter
2277 write(6,
'(a,i0)', advance=
'no')
', ', iter
2283 lgs_x_value = 0.0_dp
2284 call lis_vector_gather(lgs_x, lgs_x_value, ierr)
2285 call lis_matrix_destroy(lgs_a, ierr)
2286 call lis_vector_destroy(lgs_b, ierr)
2287 call lis_vector_destroy(lgs_x, ierr)
2288 call lis_solver_destroy(solver, ierr)
2296 vx_m(j,i) = lgs_x_value(nr)
2299 vy_m(j,i) = lgs_x_value(nr)
2303 deallocate(lgs_a_value, lgs_a_index, lgs_a_ptr)
2304 deallocate(lgs_b_value, lgs_x_value)
2308 stop
' >>> calc_vxy_ssa_matrix: Only to be called for MARGIN==3 or DYNAMICS==2!' 2312 end subroutine calc_vxy_ssa_matrix
2318 subroutine calc_vis_ssa(dxi, deta, dzeta_c, dzeta_t)
2324 real(dp),
intent(in) :: dxi, deta, dzeta_c, dzeta_t
2326 integer(i4b) :: i, j, kc, kt
2327 real(dp) :: dxi_inv, deta_inv
2328 real(dp) :: dvx_dxi, dvx_deta, dvy_dxi, dvy_deta
2329 real(dp) :: aqxy1(0:kcmax)
2330 real(dp) :: cvis0(0:ktmax), cvis1(0:kcmax)
2332 #if (MARGIN==3 || DYNAMICS==2) 2336 dxi_inv = 1.0_dp/dxi
2337 deta_inv = 1.0_dp/deta
2341 aqxy1(kc) =
aa/(
ea-1.0_dp)*
eaz_c(kc)*dzeta_c
2359 else if ((
maske(j,i)==1_i2b).or.(
maske(j,i)==2_i2b))
then 2371 dvx_dxi = (
vx_m(j,i)-
vx_m(j,i-1))*dxi_inv
2372 dvy_deta = (
vy_m(j,i)-
vy_m(j-1,i))*deta_inv
2374 dvx_deta = 0.25_dp*deta_inv &
2376 dvy_dxi = 0.25_dp*dxi_inv &
2379 de_ssa(j,i) = sqrt( dvx_dxi*dvx_dxi &
2380 + dvy_deta*dvy_deta &
2381 + dvx_dxi*dvy_deta &
2382 + 0.25_dp*(dvx_deta+dvy_dxi)*(dvx_deta+dvy_dxi) )
2391 cvis1(kc) = aqxy1(kc)*
h_c(j,i) &
2394 enh_c(kc,j,i), 0_i2b)
2401 #if (CALCMOD==-1 || CALCMOD==0) 2404 cvis1(kc) = aqxy1(kc)*
h_c(j,i) &
2407 enh_c(kc,j,i), 0_i2b)
2413 cvis0(kt) = dzeta_t*
h_t(j,i) &
2416 enh_t(kt,j,i), 1_i2b)
2420 cvis1(kc) = aqxy1(kc)*
h_c(j,i) &
2423 enh_c(kc,j,i), 0_i2b)
2426 #elif (CALCMOD==2 || CALCMOD==3) 2429 cvis1(kc) = aqxy1(kc)*
h_c(j,i) &
2432 enh_c(kc,j,i), 2_i2b)
2436 stop
' >>> calc_vis_ssa: CALCMOD must be -1, 0, 1, 2 or 3!' 2463 stop
' >>> calc_vis_ssa: Only to be called for MARGIN==3 or DYNAMICS==2!' 2467 end subroutine calc_vis_ssa
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) vx_c
vx_c(kc,j,i): Velocity in x-direction in the upper (kc) ice domain (at (i+1/2,j,kc)) ...
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(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...
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) d_help_c
d_help_c(kc,j,i): Auxiliary quantity for the computation of vx, vy und zs
real(dp), parameter epsi
epsi: Very small number
real(dp), dimension(0:kcmax) eaz_c
eaz_c(kc): Abbreviation for exp(aa*zeta(kc))
real(dp), dimension(0:jmax, 0:imax) vx_b_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) enh_c
enh_c(kc,j,i): Flow enhancement factor in the upper (kc) ice domain
real(dp), dimension(0:jmax, 0:imax) vx_m
vx_m(j,i): Mean (depth-averaged) velocity in x-direction, at (i+1/2,j)
Material properties of ice: Rate factor, heat conductivity, specific heat (heat capacity), creep function, viscosity.
real(dp), parameter eps
eps: Small number
real(dp), dimension(0:jmax, 0:imax) vx_s_g
vx_s_g(j,i): Velocity in x-direction at the ice surface, at (i,j)
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...
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
real(dp), dimension(0:jmax, 0:imax) vy_b
vy_b(j,i): Velocity in y-direction at the ice base, at (i,j+1/2)
real(dp), dimension(0:jmax, 0:imax) dzs_deta_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
logical, dimension(0:jmax, 0:imax) flag_shelfy_stream_x
flag_shelfy_stream_x(j,i): Shelfy stream flag in x-direction, at (i+1/2,j). .true.: shelfy stream point .false.: otherwise
real(dp), dimension(0:jmax, 0:imax) runoff
runoff(j,i): Runoff rate at the ice surface
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
real(dp), dimension(0:jmax, 0:imax) de_ssa
de_ssa(j,i): Effective strain rate of the SSA, at (i,j)
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) txz_c
txz_c(kc,j,i): Shear stress txz in the upper (kc) ice domain (at (i+1/2,j,kc))
integer(i2b), dimension(0:jmax, 0:imax) maske_sedi
maske_sedi(j,i): Sediment mask. 1: hard rock, 7: soft sediment, 2: ocean.
real(dp), dimension(0:jmax, 0:imax) dzs_dxi
dzs_dxi(j,i): Derivative of zs by xi (at i+1/2,j)
real(dp) function, public creep(sigma_val)
Creep response function for ice.
real(dp), dimension(0:jmax, 0:imax) zl
zl(j,i): Coordinate z of the lithosphere surface
real(dp) g
G: Acceleration due to gravity.
subroutine, public calc_vxy_static()
Computation of the horizontal velocity vx, vy, the horizontal volume flux qx, qy etc. for static ice.
real(dp), dimension(0:jmax, 0:imax) vy_m
vy_m(j,i): Mean (depth-averaged) velocity in y-direction, at (i,j+1/2)
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).
real(dp), dimension(0:jmax, 0:imax) vis_int_g
vis_int_g(j,i): Depth-integrated viscosity of the SSA, at (i,j)
real(dp), dimension(0:jmax, 0:imax) ratio_sl_y
ratio_sl_y(j,i): Ratio of basal to surface velocity (slip ratio) in y-direction, at (i...
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
real(dp), dimension(0:jmax, 0:imax) d_help_b
d_help_b(j,i): Auxiliary quantity for the computation of vx_b and vy_b
real(dp) function, public ratefac_c(temp_val, temp_m_val)
Rate factor for cold ice: Linear interpolation of tabulated values in RF(.).
real(dp) ea
ea: Abbreviation for exp(aa)
logical flag_aa_nonzero
flag_aa_nonzero: Flag for the exponential stretch parameter aa. .true.: aa greater than zero (non-equ...
real(dp), dimension(0:jmax, 0:imax) h_diff
h_diff(j,i): Diffusivity of the SIA evolution equation of the ice surface
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) d_help_t
d_help_t(kt,j,i): Auxiliary quantity for the computation of vx, vy und zs
real(dp) aa
aa: Exponential stretch parameter of the non-equidistant vertical grid in the upper (kc) ice domain ...
Declarations of kind types for SICOPOLIS.
real(dp), dimension(0:jmax, 0:imax) flui_ave_sia
flui_ave_sia(j,i): Depth-averaged fluidity of the SIA
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) txz_t
txz_t(kt,j,i): Shear stress txz in the lower (kt) ice domain (at (i+1/2,j,kt))
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) enh_t
enh_t(kt,j,i): Flow enhancement factor in the lower (kt) ice domain
real(dp) function, public ratefac_t(omega_val)
Rate factor for temperate ice.
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) vy_b_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) omega_c
omega_c(kc,j,i): Water content in the upper (kc) ice domain
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) temp_c_m
temp_c_m(kc,j,i): Melting temperature in the upper (kc) ice domain
real(dp), dimension(0:jmax, 0:imax) qx
qx(j,i): Volume flux in x-direction (depth-integrated vx, at (i+1/2,j))
real(dp), parameter eps_dp
eps_dp: Small number to double-precision accuracy
real(dp), dimension(0:jmax, 0:imax) p_weert_inv
p_weert_inv(j,i): Inverse of p_weert
logical, dimension(0:jmax, 0:imax) flag_shelfy_stream
flag_shelfy_stream(j,i): Shelfy stream flag on the main grid. .true.: grounded ice, and at least one neighbour on the staggered grid is a shelfy stream point .false.: otherwise
real(dp), dimension(0:jmax, 0:imax) dzs_deta
dzs_deta(j,i): Derivative of zs by eta (at i,j+1/2)
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) tyz_c
tyz_c(kc,j,i): Shear stress tyz in the upper (kc) ice domain (at (i,j+1/2,kc))
real(dp), dimension(0:jmax, 0:imax) q_gl_g
q_gl_g(j,i): Volume flux across the grounding line, at (i,j)
real(dp), dimension(0:jmax, 0:imax) vy_s_g
vy_s_g(j,i): Velocity in x-direction at the ice surface, at (i,j)
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(0:kcmax) eaz_c_quotient
eaz_c_quotient(kc): Abbreviation for (eaz_c(kc)-1.0)/(ea-1.0)
real(dp), dimension(0:jmax, 0:imax) dzs_dxi_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) vy_t
vy_t(kt,j,i): Velocity in y-direction in the lower (kt) ice domain (at (i,j+1/2,kt)) ...
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...
real(dp), dimension(0:jmax, 0:imax) ratio_sl_x
ratio_sl_x(j,i): Ratio of basal to surface velocity (slip ratio) in x-direction, at (i+1/2...
real(dp), dimension(0:jmax, 0:imax) c_drag
c_drag(j,i): Auxiliary quantity for the computation of the basal drag
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) vy_c
vy_c(kc,j,i): Velocity in y-direction in the upper (kc) ice domain (at (i,j+1/2,kc)) ...
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) omega_t
omega_t(kt,j,i): Water content in the lower (kt) ice domain
real(dp), dimension(0:ktmax) zeta_t
zeta_t(kt): Sigma coordinate zeta_t of grid point kt
real(dp) rho_sw
RHO_SW: Density of sea water.
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...
real(dp) function, public viscosity(de_val, temp_val, temp_m_val, omega_val, enh_val, i_flag_cold_temp)
Ice viscosity as a function of the effective strain rate and the temperature (in cold ice) or the wat...
logical, dimension(0:jmax, 0:imax) flag_shelfy_stream_y
flag_shelfy_stream_y(j,i): Shelfy stream flag in y-direction, at (i,j+1/2). .true.: shelfy stream point .false.: otherwise
real(dp) function, public ratefac_c_t(temp_val, omega_val, temp_m_val)
Rate factor for cold and temperate ice: Combination of ratefac_c and ratefac_t (only for the enthalpy...
real(dp) rho
RHO: Density of ice.
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) sigma_c
sigma_c(kc,j,i): Effective stress in the upper (kc) ice domain
real(dp), dimension(0:jmax, 0:imax) vx_b
vx_b(j,i): Velocity in x-direction at the ice base, at (i+1/2,j)
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) sigma_t
sigma_t(kt,j,i): Effective stress in the lower (kt) ice domain
integer(i4b), dimension(0:jmax, 0:imax) p_weert
p_weert(j,i): Weertman exponent for the basal shear stress
integer, parameter dp
Double-precision reals.
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) 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) 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:ktmax, 0:jmax, 0:imax) tyz_t
tyz_t(kt,j,i): Shear stress tyz in the lower (kt) ice domain (at (i,j+1/2,kt))
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) vx_t
vx_t(kt,j,i): Velocity in x-direction in the lower (kt) ice domain (at (i+1/2,j,kt)) ...
real(dp), dimension(0:jmax, 0:imax) qy
qy(j,i): Volume flux in y-direction (depth-integrated vy, at (i,j+1/2))