53 subroutine calc_qbm(time, z_sl, dzeta_c, dzeta_r)
59 real(dp),
intent(in) :: time, z_sl
60 real(dp),
intent(in) :: dzeta_c, dzeta_r
63 integer(i4b) :: n_ocean, n_float
64 real(dp) :: aqbm1, aqbm2, aqbm3a, aqbm3b, aqbm4
65 real(dp) :: frictional_heating
66 real(dp) :: year_sec_inv, rhow_rho_ratio
67 real(dp) :: Q_bm_grounded, Q_bm_marine, Q_bm_floating
68 real(dp) :: qbm_min, qbm_max
72 year_sec_inv = 1.0_dp/year_sec
76 aqbm1 = (
ea-1.0_dp)/(
rho*
l*
aa*dzeta_c)
78 aqbm1 = 1.0_dp/(
rho*
l*dzeta_c)
83 aqbm3b = 1.0_dp/(
rho*
l)
93 if (
maske(j,i)==0_i1b)
then 95 if (
n_cts(j,i)==-1_i1b)
then 97 frictional_heating = 0.0_dp
100 else if (
n_cts(j,i)==0_i1b)
then 106 = -aqbm3a*
h_c(j,i)*0.5_dp &
131 = -aqbm3a*(
h_c(j,i)+
h_t(j,i))*0.5_dp &
151 #if (MARGIN==1 || MARGIN==2) 153 #if (!defined(MARINE_ICE_BASAL_MELTING) || MARINE_ICE_BASAL_MELTING==1) 157 #elif (MARINE_ICE_BASAL_MELTING==2) 159 if ( (
zb(j,i) < z_sl) &
161 ( (
maske(j,i+1)==2_i1b) &
162 .or.(
maske(j,i-1)==2_i1b) &
163 .or.(
maske(j+1,i)==2_i1b) &
164 .or.(
maske(j-1,i)==2_i1b) &
168 q_bm(j,i) = qbm_marine *year_sec_inv*rhow_rho_ratio
173 #elif (MARINE_ICE_BASAL_MELTING==3) 175 if ( (
zb(j,i) < z_sl) &
177 ( (
maske(j,i+1)==2_i1b) &
178 .or.(
maske(j,i-1)==2_i1b) &
179 .or.(
maske(j+1,i)==2_i1b) &
180 .or.(
maske(j-1,i)==2_i1b) &
185 if (
maske(j,i+1)==2_i1b) n_ocean = n_ocean+1
186 if (
maske(j,i-1)==2_i1b) n_ocean = n_ocean+1
187 if (
maske(j+1,i)==2_i1b) n_ocean = n_ocean+1
188 if (
maske(j-1,i)==2_i1b) n_ocean = n_ocean+1
190 if ( n_ocean > 0 )
then 192 q_bm_grounded =
q_bm(j,i)
193 q_bm_marine = qbm_marine *year_sec_inv*rhow_rho_ratio
196 q_bm(j,i) = (1.0_dp-0.25_dp*
real(n_ocean,
dp)) * q_bm_grounded &
197 +0.25_dp*
real(n_ocean,dp) * Q_bm_marine
201 errormsg =
' >>> calc_qbm: Marine ice margin point does not' &
203 //
' have an ocean neighbour!' 210 errormsg =
' >>> calc_qbm: MARINE_ICE_BASAL_MELTING must be 1, 2 or 3!' 219 #if (FLOATING_ICE_BASAL_MELTING==1 || FLOATING_ICE_BASAL_MELTING==4 \ 220 || floating_ice_basal_melting==5)
224 #elif (FLOATING_ICE_BASAL_MELTING==2) 226 q_bm(j,i) = qbm_float_1 *year_sec_inv*rhow_rho_ratio
229 #elif (FLOATING_ICE_BASAL_MELTING==3) 232 if (
maske(j,i+1)>=2_i1b) n_float = n_float+1
233 if (
maske(j,i-1)>=2_i1b) n_float = n_float+1
234 if (
maske(j+1,i)>=2_i1b) n_float = n_float+1
235 if (
maske(j-1,i)>=2_i1b) n_float = n_float+1
237 if ( n_float > 0 )
then 239 q_bm_grounded =
q_bm(j,i)
240 q_bm_floating = qbm_float_1 *year_sec_inv*rhow_rho_ratio
243 q_bm(j,i) = (1.0_dp-0.25_dp*
real(n_float,
dp)) * q_bm_grounded &
244 +0.25_dp*
real(n_float,dp) * Q_bm_floating
247 errormsg =
' >>> calc_qbm: Grounding line point does not have' &
249 //
' a floating ice or ocean neighbour!' 254 errormsg =
' >>> calc_qbm: FLOATING_ICE_BASAL_MELTING' &
256 //
' must be 1, 2, 3, 4 or 5!' 260 else if ( (
zb(j,i) < z_sl) &
262 ( (
maske(j,i+1)>=2_i1b) &
263 .or.(
maske(j,i-1)>=2_i1b) &
264 .or.(
maske(j+1,i)>=2_i1b) &
265 .or.(
maske(j-1,i)>=2_i1b) &
269 #if (MARINE_ICE_BASAL_MELTING==1) 273 #elif (MARINE_ICE_BASAL_MELTING==2) 275 q_bm(j,i) = qbm_marine *year_sec_inv*rhow_rho_ratio
278 #elif (MARINE_ICE_BASAL_MELTING==3) 281 if (
maske(j,i+1)>=2_i1b) n_ocean = n_ocean+1
282 if (
maske(j,i-1)>=2_i1b) n_ocean = n_ocean+1
283 if (
maske(j+1,i)>=2_i1b) n_ocean = n_ocean+1
284 if (
maske(j-1,i)>=2_i1b) n_ocean = n_ocean+1
286 if ( n_ocean > 0 )
then 288 q_bm_grounded =
q_bm(j,i)
289 q_bm_marine = qbm_marine *year_sec_inv*rhow_rho_ratio
292 q_bm(j,i) = (1.0_dp-0.25_dp*
real(n_ocean,
dp)) * q_bm_grounded &
293 +0.25_dp*
real(n_ocean,dp) * Q_bm_marine
297 errormsg =
' >>> calc_qbm: Marine ice margin point does not' &
299 //
' have a floating ice or ocean neighbour!' 304 errormsg =
' >>> calc_qbm: MARINE_ICE_BASAL_MELTING must be 1, 2 or 3!' 310 else if ( (
maske(j,i)==2_i1b).or.(
maske(j,i)==3_i1b) )
then 313 #if (FLOATING_ICE_BASAL_MELTING<=3) 318 q_bm(j,i) = qbm_float_1 *year_sec_inv*rhow_rho_ratio
321 else if (
zl(j,i) > z_abyss )
then 323 q_bm(j,i) = qbm_float_2 *year_sec_inv*rhow_rho_ratio
328 q_bm(j,i) = qbm_float_3 *year_sec_inv*rhow_rho_ratio
333 #elif (FLOATING_ICE_BASAL_MELTING==4 || FLOATING_ICE_BASAL_MELTING==5) 335 if (
zl(j,i) > z_abyss )
then 337 year_sec_inv, rhow_rho_ratio, &
339 q_bm(j,i) = q_bm_floating
341 q_bm(j,i) = qbm_float_3 *year_sec_inv*rhow_rho_ratio
346 errormsg =
' >>> calc_qbm: FLOATING_ICE_BASAL_MELTING' &
348 //
' must be 1, 2, 3, 4 or 5!' 353 errormsg =
' >>> calc_qbm: MARGIN must be 1, 2 or 3!' 369 #if (defined(QBM_MIN)) 370 qbm_min = qbm_min *year_sec_inv*rhow_rho_ratio
372 #elif (defined(QB_MIN)) /* obsolete */ 375 errormsg =
' >>> calc_qbm: Limiter QBM_MIN is undefined!' 379 #if (defined(QBM_MAX)) 380 qbm_max = qbm_max *year_sec_inv*rhow_rho_ratio
382 #elif (defined(QB_MAX)) /* obsolete */ 385 errormsg =
' >>> calc_qbm: Limiter QBM_MAX is undefined!' 389 #if ( defined(MARINE_ICE_BASAL_MELTING) \ 390 && ( marine_ice_basal_melting==2 || marine_ice_basal_melting==3 ) )
392 if (qbm_marine*year_sec_inv*rhow_rho_ratio > qbm_max)
then 393 errormsg =
' >>> calc_qbm: QBM_MARINE' &
395 //
' (basal melting rate at the ice front)' &
397 //
' is larger than the limiter qbm_max!' 404 && defined(floating_ice_basal_melting) \
405 && floating_ice_basal_melting<=3 )
407 if (qbm_float_1*year_sec_inv*rhow_rho_ratio > qbm_max)
then 408 errormsg =
' >>> calc_qbm: QBM_FLOAT_1' &
410 //
' (basal melting rate in the grounding line zone)' &
412 //
' is larger than the limiter qbm_max!' 420 if (
maske(j,i)==0_i1b)
then 421 if (
q_bm(j,i) < qbm_min)
q_bm(j,i) = 0.0_dp
422 if (
q_bm(j,i) > qbm_max)
q_bm(j,i) = qbm_max
424 if (
q_tld(j,i) > qbm_max)
q_tld(j,i) = qbm_max
437 year_sec_inv, rhow_rho_ratio, &
440 #if defined(ALLOW_OPENAD) /* OpenAD */ 446 #if !defined(ALLOW_OPENAD) /* Normal */ 447 integer(i4b),
intent(in) :: i, j
449 integer(i4b),
intent(inout) :: i, j
450 #endif /* Normal vs. OpenAD */ 451 real(dp),
intent(in) :: time, z_sl
452 real(dp),
intent(in) :: year_sec_inv, rhow_rho_ratio
454 real(dp),
intent(out) :: Q_bm_floating
457 real(dp) :: time_in_years
458 real(dp) :: Toc, Tmb, T_forcing, Omega, alpha, draft0, draft
459 real(dp) :: lon_d, lat_d, Phi_par
460 real(dp) :: H_w_now, Q_bm_scaling_factor
462 real(dp),
parameter :: beta_sw = 7.61e-04_dp
464 #if (FLOATING_ICE_BASAL_MELTING==5) 466 real(dp),
parameter :: c_sw = 3974.0_dp, &
471 #if defined(ALLOW_OPENAD) /* OpenAD */ 472 integer(i4b) :: i_time_in_years
475 time_in_years = time*year_sec_inv
477 #if (FLOATING_ICE_BASAL_MELTING==4) 480 omega = omega_qbm *year_sec_inv*rhow_rho_ratio
484 #elif (FLOATING_ICE_BASAL_MELTING==5) 486 #if (defined(ANT)) /* Antarctic ice sheet */ 493 #if !defined(ALLOW_OPENAD) /* Normal */ 494 lon_d = modulo(lon_d+180.0_dp, 360.0_dp)-180.0_dp
498 lon_d = (lon_d+180.0_sp) - &
499 ((360.0_sp * int(lon_d+180.0_sp)/360.0_sp)) - &
501 #endif /* Normal vs. OpenAD */ 503 if ((lon_d>=-10.0_dp).and.(lon_d<60.0_dp))
then 512 else if ((lon_d>=60.0_dp).and.(lon_d<80.0_dp))
then 521 else if ((lon_d>=80.0_dp).and.(lon_d<130.0_dp))
then 531 else if ( ((lon_d>=130.0_dp).and.(lon_d<159.0_dp)) &
533 ((lon_d>=159.0_dp).and.(lon_d<170.0_dp) &
534 .and.(lat_d>=-72.0_dp)) )
then 544 else if ( (lon_d>=159.0_dp) &
548 ((lon_d>=-140.0_dp).and.(lon_d<-120.0_dp) &
549 .and.(lat_d<-77.0_dp)) )
then 558 else if ( ((lon_d>=-140.0_dp).and.(lon_d<-120.0_dp) &
559 .and.(lat_d>=-77.0_dp)) &
561 ((lon_d>=-120.0_dp).and.(lon_d<-90.0_dp)) )
then 570 else if ( ((lon_d>=-90.0_dp).and.(lon_d<-66.0_dp) &
571 .and.(lat_d>=-74.0_dp)) )
then 591 #else /* not Antarctic ice sheet */ 593 q_bm_floating = 0.0_dp
595 errormsg =
' >>> sub_ice_shelf_melting_param:' &
597 //
' Case FLOATING_ICE_BASAL_MELTING==5' &
599 //
' only defined for Antarctica!' 612 if (
maske(j,i)==2_i1b)
then 614 h_w_now = max((z_sl-
zl(j,i)), 0.0_dp)
615 else if (
maske(j,i)==3_i1b)
then 616 draft = max((z_sl -
zb(j,i)), 0.0_dp)
617 h_w_now = max((
zb(j,i)-
zl(j,i)), 0.0_dp)
619 errormsg =
' >>> sub_ice_shelf_melting_param:' &
621 //
' Routine must not be called for maske(j,i) < 2!' 627 t_forcing = max((toc-tmb), 0.0_dp)
631 if (h_w_0 >
eps)
then 632 q_bm_scaling_factor = tanh(h_w_now/h_w_0)
634 q_bm_scaling_factor = 1.0_dp
638 q_bm_scaling_factor = 1.0_dp
641 #if (FLOATING_ICE_BASAL_MELTING==4) 643 q_bm_floating = q_bm_scaling_factor*omega*t_forcing**alpha
645 #elif (FLOATING_ICE_BASAL_MELTING==5) 649 #if !defined(ALLOW_OPENAD) 650 q_bm_floating = q_bm_scaling_factor &
651 *phi_par*omega*t_forcing*(draft/draft0)**alpha
654 q_bm_floating = 0.0_dp
656 q_bm_floating = q_bm_scaling_factor &
657 *phi_par*omega*t_forcing*(draft/draft0)**alpha
662 q_bm_floating = 0.0_dp
664 errormsg =
' >>> sub_ice_shelf_melting_param:' &
666 //
' FLOATING_ICE_BASAL_MELTING must be 4 or 5!' 673 #if (defined(INITMIP_BMB_ANOM_FILE)) 675 if ((time_in_years > 0.0_dp).and.(time_in_years <= 40.0_dp))
then 677 q_bm_floating = q_bm_floating &
678 #if !defined(ALLOW_OPENAD) /* Normal */ 679 + 0.025_dp*floor(time_in_years) * ab_anom_initmip(j,i)
681 call myfloor(time_in_years,i_time_in_years)
682 + 0.025_dp*i_time_in_years * ab_anom_initmip(j,i)
683 #endif /* Normal vs. OpenAD */ 685 else if (time_in_years > 40.0_dp)
then 687 q_bm_floating = q_bm_floating + ab_anom_initmip(j,i)
695 #if (defined(LARMIP_REGIONS_FILE)) 696 n = int(n_regions_larmip(j,i))
697 q_bm_floating = q_bm_floating + ab_anom_larmip(n)
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) kappa_r
KAPPA_R: Heat conductivity of the lithosphere.
real(dp) delta_tm_sw
DELTA_TM_SW: Melting point depression of sea water due to its average salinity.
integer(i1b), dimension(0:jmax, 0:imax) n_cts
n_cts(j,i): Mask for thermal conditions. -1: cold ice base, 0: temperate ice base with cold ice above...
real(dp), dimension(0:jmax, 0:imax) phi
phi(j,i): Geographic latitude of grid point (i,j)
subroutine, public calc_qbm(time, z_sl, dzeta_c, dzeta_r)
Computation of the basal melting rate Q_bm. Summation of Q_bm and Q_tld (water drainage rate from the...
real(dp), dimension(0:jmax, 0:imax) vx_b_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
real(dp) rho_w
RHO_W: Density of pure water.
real(dp), dimension(0:krmax, 0:jmax, 0:imax) temp_r
temp_r(kr,j,i): Temperature in the bedrock
real(dp) h_r
H_R: Thickness of the modelled lithosphere layer.
subroutine error(error_message)
Main routine of error_m: Writing of error messages and stopping execution.
Material properties of ice: Rate factor, heat conductivity, specific heat (heat capacity), creep function, viscosity.
real(dp), parameter eps
eps: Small number
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:jmax, 0:imax) lambda
lambda(j,i): Geographic longitude of grid point (i,j)
real(dp) function, public kappa_val(temp_val)
Heat conductivity of ice: Linear interpolation of tabulated values in KAPPA(.).
real(dp), dimension(0:jmax, 0:imax) dzs_deta_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
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), parameter pi_180_inv
pi_180_inv: 180 divided by pi (-> rad to deg)
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.
real(dp), dimension(0:jmax, 0:imax) q_tld
Q_tld(j,i): Water drainage rate from the temperate layer.
Declarations of global variables for SICOPOLIS (for the ANT domain).
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) 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) 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)
Writing of error messages and stopping execution.
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
character, parameter end_of_line
end_of_line: End-of-line string
integer(i1b), dimension(0:jmax, 0:imax) maske
maske(j,i): Ice-land-ocean mask. 0: grounded ice, 1: ice-free land, 2: ocean, 3: floating ice ...
real(dp), dimension(0:jmax, 0:imax) q_b_tot
Q_b_tot(j,i): Sum of Q_bm and Q_tld.
real(dp), dimension(0:jmax, 0:imax) q_bm
Q_bm(j,i): Basal melting rate.
real(dp), 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)) ...
real(dp) beta
BETA: Clausius-Clapeyron gradient of ice.
real(dp), dimension(0:jmax, 0:imax) c_drag
c_drag(j,i): Auxiliary quantity for the computation of the basal drag
Computation of the basal melting rate.
real(dp) rho_sw
RHO_SW: Density of sea water.
real(dp) l
L: Latent heat of ice.
subroutine sub_ice_shelf_melting_param(time, z_sl, year_sec_inv, rhow_rho_ratio, i, j, Q_bm_floating)
Sub-ice-shelf melting parameterisation.
character(len=256) errormsg
errormsg: Error-message string
real(dp) rho
RHO: Density of ice.
integer(i1b), dimension(0:jmax, 0:imax) n_sector
n_sector(j,i): Marker for the different sectors for ice shelf basal melting.
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...
Declarations of global variables for SICOPOLIS.
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)) ...