SICOPOLIS V5-dev  Revision 1105
calc_bas_melt_m.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Module : c a l c _ b a s _ m e l t _ m
4 !
5 !> @file
6 !!
7 !! Computation of the basal melting rate.
8 !!
9 !! @section Copyright
10 !!
11 !! Copyright 2009-2017 Ralf Greve, Ben Galton-Fenzi, Tatsuru Sato
12 !!
13 !! @section License
14 !!
15 !! This file is part of SICOPOLIS.
16 !!
17 !! SICOPOLIS is free software: you can redistribute it and/or modify
18 !! it under the terms of the GNU General Public License as published by
19 !! the Free Software Foundation, either version 3 of the License, or
20 !! (at your option) any later version.
21 !!
22 !! SICOPOLIS is distributed in the hope that it will be useful,
23 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
24 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25 !! GNU General Public License for more details.
26 !!
27 !! You should have received a copy of the GNU General Public License
28 !! along with SICOPOLIS. If not, see <http://www.gnu.org/licenses/>.
29 !<
30 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
31 
32 !-------------------------------------------------------------------------------
33 !> Computation of the basal melting rate.
34 !<------------------------------------------------------------------------------
35 
37 
38  use sico_types_m
40  use sico_vars_m
41 
42  implicit none
43 
44  private
45  public :: calc_qbm
46 
47 contains
48 
49 !-------------------------------------------------------------------------------
50 !> Computation of the basal melting rate Q_bm.
51 !! Summation of Q_bm and Q_tld (water drainage rate from the temperate layer).
52 !<------------------------------------------------------------------------------
53 subroutine calc_qbm(time, z_sl, dzeta_c, dzeta_r)
54 
56 
57 implicit none
58 
59 real(dp), intent(in) :: time, z_sl
60 real(dp), intent(in) :: dzeta_c, dzeta_r
61 
62 integer(i4b) :: i, j
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
69 
70 !-------- Term abbreviations --------
71 
72 year_sec_inv = 1.0_dp/year_sec
73 rhow_rho_ratio = rho_w/rho
74 
75 if (flag_aa_nonzero) then
76  aqbm1 = (ea-1.0_dp)/(rho*l*aa*dzeta_c)
77 else
78  aqbm1 = 1.0_dp/(rho*l*dzeta_c)
79 end if
80 
81 aqbm2 = kappa_r/(rho*l*h_r*dzeta_r)
82 aqbm3a = g/l
83 aqbm3b = 1.0_dp/(rho*l)
84 aqbm4 = beta/(rho*l)
85 
86 !-------- Computation of Q_bm --------
87 
88 q_bm = 0.0_dp ! initialisation
89 
90 do i=1, imax-1
91 do j=1, jmax-1
92 
93  if (maske(j,i)==0) then ! grounded ice
94 
95  if (n_cts(j,i)==-1) then
96 
97  frictional_heating = 0.0_dp
98  q_bm(j,i) = 0.0_dp
99 
100  else if (n_cts(j,i)==0) then
101 
102 #if (DYNAMICS==2)
103  if (.not.flag_shelfy_stream(j,i)) then
104 #endif
105  frictional_heating &
106  = -aqbm3a*h_c(j,i) &
107  *0.5_dp*(vx_t(0,j,i)+vx_t(0,j,i-1)) &
108  *dzs_dxi_g(j,i) &
109  -aqbm3a*h_c(j,i) &
110  *0.5_dp*(vy_t(0,j,i)+vy_t(0,j-1,i)) &
111  *dzs_deta_g(j,i)
112 #if (DYNAMICS==2)
113  else ! flag_shelfy_stream(j,i) == .true.
114 
115  frictional_heating &
116  = aqbm3b &
117  * c_drag(j,i) &
118  * sqrt(vx_b_g(j,i)**2 &
119  +vy_b_g(j,i)**2) &
120  **(1.0_dp+p_weert_inv(j,i))
121  end if
122 #endif
123  q_bm(j,i) = aqbm1*(temp_c(1,j,i)-temp_c(0,j,i))/h_c(j,i) &
124  *kappa_val(temp_c(0,j,i)) &
125  - aqbm2*(temp_r(krmax,j,i)-temp_r(krmax-1,j,i)) &
126  + frictional_heating
127 
128  else ! n_cts(j,i)==1
129 
130 #if (DYNAMICS==2)
131  if (.not.flag_shelfy_stream(j,i)) then
132 #endif
133  frictional_heating &
134  = -aqbm3a*(h_c(j,i)+h_t(j,i)) &
135  *0.5_dp*(vx_t(0,j,i)+vx_t(0,j,i-1)) &
136  *dzs_dxi_g(j,i) &
137  -aqbm3a*(h_c(j,i)+h_t(j,i)) &
138  *0.5_dp*(vy_t(0,j,i)+vy_t(0,j-1,i)) &
139  *dzs_deta_g(j,i)
140 #if (DYNAMICS==2)
141  else ! flag_shelfy_stream(j,i) == .true.
142 
143  frictional_heating &
144  = aqbm3b &
145  * c_drag(j,i) &
146  * sqrt(vx_b_g(j,i)**2 &
147  +vy_b_g(j,i)**2) &
148  **(1.0_dp+p_weert_inv(j,i))
149  end if
150 #endif
151  q_bm(j,i) = aqbm4*kappa_val(temp_t_m(0,j,i)) &
152  - aqbm2*(temp_r(krmax,j,i)-temp_r(krmax-1,j,i)) &
153  + frictional_heating
154 
155  end if
156 
157 #if (MARGIN==1 || MARGIN==2)
158 
159 #if (!defined(MARINE_ICE_BASAL_MELTING) || MARINE_ICE_BASAL_MELTING==1)
160 
161  continue ! do nothing
162 
163 #elif (MARINE_ICE_BASAL_MELTING==2)
164 
165  if ( (zb(j,i) < z_sl) & ! marine ice
166  .and. &
167  ( (maske(j,i+1)==2) & ! at least one
168  .or.(maske(j,i-1)==2) & ! nearest neighbour
169  .or.(maske(j+1,i)==2) & ! is
170  .or.(maske(j-1,i)==2) & ! ocean
171  ) &
172  ) then
173 
174  q_bm(j,i) = qbm_marine *year_sec_inv*rhow_rho_ratio
175  ! m/a water equiv. -> m/s ice equiv.
176 
177  end if
178 
179 #elif (MARINE_ICE_BASAL_MELTING==3)
180 
181  if ( (zb(j,i) < z_sl) & ! marine ice
182  .and. &
183  ( (maske(j,i+1)==2) & ! at least one
184  .or.(maske(j,i-1)==2) & ! nearest neighbour
185  .or.(maske(j+1,i)==2) & ! is
186  .or.(maske(j-1,i)==2) & ! ocean
187  ) &
188  ) then
189 
190  n_ocean = 0
191  if (maske(j,i+1)==2) n_ocean = n_ocean+1
192  if (maske(j,i-1)==2) n_ocean = n_ocean+1
193  if (maske(j+1,i)==2) n_ocean = n_ocean+1
194  if (maske(j-1,i)==2) n_ocean = n_ocean+1
195 
196  if ( n_ocean > 0 ) then
197 
198  q_bm_grounded = q_bm(j,i)
199  q_bm_marine = qbm_marine *year_sec_inv*rhow_rho_ratio
200  ! m/a water equiv. -> m/s ice equiv.
201 
202  q_bm(j,i) = (1.0_dp-0.25_dp*real(n_ocean,dp)) * q_bm_grounded &
203  +0.25_dp*real(n_ocean,dp) * Q_bm_marine
204  ! weighed average of grounded ice melting (computed)
205  ! and marine ice melting (prescribed)
206  else
207  write(6,'(a)') ' >>> calc_qbm: Marine ice margin point does not'
208  write(6,'(a)') ' have an ocean neighbour!'
209  stop
210  end if
211 
212  end if
213 
214 #else
215  stop ' >>> calc_qbm: MARINE_ICE_BASAL_MELTING must be 1, 2 or 3!'
216 #endif
217 
218 #elif (MARGIN==3)
219 
220  if (flag_grounding_line_1(j,i)) then
221  ! grounding line (grounded-ice side)
222 
223 #if (FLOATING_ICE_BASAL_MELTING==1 || FLOATING_ICE_BASAL_MELTING==4 \
224  || floating_ice_basal_melting==5)
225 
226  continue ! do nothing
227 
228 #elif (FLOATING_ICE_BASAL_MELTING==2)
229 
230  q_bm(j,i) = qbm_float_1 *year_sec_inv*rhow_rho_ratio
231  ! m/a water equiv. -> m/s ice equiv.
232 
233 #elif (FLOATING_ICE_BASAL_MELTING==3)
234 
235  n_float = 0
236  if (maske(j,i+1)>=2) n_float = n_float+1
237  if (maske(j,i-1)>=2) n_float = n_float+1
238  if (maske(j+1,i)>=2) n_float = n_float+1
239  if (maske(j-1,i)>=2) n_float = n_float+1
240 
241  if ( n_float > 0 ) then
242 
243  q_bm_grounded = q_bm(j,i)
244  q_bm_floating = qbm_float_1 *year_sec_inv*rhow_rho_ratio
245  ! m/a water equiv. -> m/s ice equiv.
246 
247  q_bm(j,i) = (1.0_dp-0.25_dp*real(n_float,dp)) * q_bm_grounded &
248  +0.25_dp*real(n_float,dp) * Q_bm_floating
249  ! weighed average of melting for grounded and floating ice
250  else
251  write(6,'(a)') ' >>> calc_qbm: Grounding line point does not have'
252  write(6,'(a)') ' a floating ice or ocean neighbour!'
253  stop
254  end if
255 
256 #else
257  write(6, fmt='(a)') ' >>> calc_qbm: FLOATING_ICE_BASAL_MELTING'
258  write(6, fmt='(a)') ' must be 1, 2, 3, 4 or 5!'
259  stop
260 #endif
261 
262  else if ( (zb(j,i) < z_sl) & ! marine ice margin
263  .and. &
264  ( (maske(j,i+1)>=2) & ! (at least one
265  .or.(maske(j,i-1)>=2) & ! nearest neighbour
266  .or.(maske(j+1,i)>=2) & ! is
267  .or.(maske(j-1,i)>=2) & ! ocean)
268  ) &
269  ) then
270 
271 #if (MARINE_ICE_BASAL_MELTING==1)
272 
273  continue ! do nothing
274 
275 #elif (MARINE_ICE_BASAL_MELTING==2)
276 
277  q_bm(j,i) = qbm_marine *year_sec_inv*rhow_rho_ratio
278  ! m/a water equiv. -> m/s ice equiv.
279 
280 #elif (MARINE_ICE_BASAL_MELTING==3)
281 
282  n_ocean = 0
283  if (maske(j,i+1)>=2) n_ocean = n_ocean+1
284  if (maske(j,i-1)>=2) n_ocean = n_ocean+1
285  if (maske(j+1,i)>=2) n_ocean = n_ocean+1
286  if (maske(j-1,i)>=2) n_ocean = n_ocean+1
287 
288  if ( n_ocean > 0 ) then
289 
290  q_bm_grounded = q_bm(j,i)
291  q_bm_marine = qbm_marine *year_sec_inv*rhow_rho_ratio
292  ! m/a water equiv. -> m/s ice equiv.
293 
294  q_bm(j,i) = (1.0_dp-0.25_dp*real(n_ocean,dp)) * q_bm_grounded &
295  +0.25_dp*real(n_ocean,dp) * Q_bm_marine
296  ! weighed average of grounded ice melting (computed)
297  ! and marine ice melting (prescribed)
298  else
299  write(6,'(a)') ' >>> calc_qbm: Marine ice margin point does not'
300  write(6,'(a)') ' have a floating ice or ocean neighbour!'
301  stop
302  end if
303 
304 #else
305  stop ' >>> calc_qbm: MARINE_ICE_BASAL_MELTING must be 1, 2 or 3!'
306 #endif
307 
308  end if
309 
310  else if ( (maske(j,i)==2).or.(maske(j,i)==3) ) then ! floating ice or ocean
311 
312 #if (FLOATING_ICE_BASAL_MELTING<=3)
313 
314  if (flag_grounding_line_2(j,i)) then
315  ! grounding line (floating-ice side)
316 
317  q_bm(j,i) = qbm_float_1 *year_sec_inv*rhow_rho_ratio
318  ! m/a water equiv. -> m/s ice equiv.
319 
320  else if ( zl(j,i) > z_abyss ) then ! floating ice over continental shelf
321 
322  q_bm(j,i) = qbm_float_2 *year_sec_inv*rhow_rho_ratio
323  ! m/a water equiv. -> m/s ice equiv.
324 
325  else ! ( zl(j,i) <= Z_ABYSS ), floating ice over abyssal ocean
326 
327  q_bm(j,i) = qbm_float_3 *year_sec_inv*rhow_rho_ratio
328  ! m/a water equiv. -> m/s ice equiv.
329 
330  end if
331 
332 #elif (FLOATING_ICE_BASAL_MELTING==4 || FLOATING_ICE_BASAL_MELTING==5)
333 
334  if ( zl(j,i) > z_abyss ) then ! floating ice over continental shelf
335  call sub_ice_shelf_melting_param(time, z_sl, &
336  year_sec_inv, rhow_rho_ratio, &
337  i, j, q_bm_floating)
338  q_bm(j,i) = q_bm_floating
339  else ! ( zl(j,i) <= Z_ABYSS ), floating ice over abyssal ocean
340  q_bm(j,i) = qbm_float_3 *year_sec_inv*rhow_rho_ratio
341  ! m/a water equiv. -> m/s ice equiv.
342  end if
343 
344 #else
345  write(6, fmt='(a)') ' >>> calc_qbm: FLOATING_ICE_BASAL_MELTING'
346  write(6, fmt='(a)') ' must be 1, 2, 3, 4 or 5!'
347  stop
348 #endif
349 
350 #else
351  stop ' >>> calc_qbm: MARGIN must be 1, 2 or 3!'
352 #endif
353 
354  end if
355 
356 end do
357 end do
358 
359 !-------- Sum of Q_bm and Q_tld --------
360 
361 q_b_tot = q_bm + q_tld
362 
363 !-------- Limitation of Q_bm, Q_tld and Q_b_tot
364 ! (only for grounded ice) --------
365 
366 #if (defined(QBM_MIN))
367  qbm_min = qbm_min *year_sec_inv*rhow_rho_ratio
368  ! m/a water equiv. -> m/s ice equiv.
369 #elif (defined(QB_MIN)) /* obsolete */
370  qbm_min = qb_min ! m/s ice equiv.
371 #else
372  stop ' >>> calc_qbm: Limiter QBM_MIN is undefined!'
373 #endif
374 
375 #if (defined(QBM_MAX))
376  qbm_max = qbm_max *year_sec_inv*rhow_rho_ratio
377  ! m/a water equiv. -> m/s ice equiv.
378 #elif (defined(QB_MAX)) /* obsolete */
379  qbm_max = qb_max ! m/s ice equiv.
380 #else
381  stop ' >>> calc_qbm: Limiter QBM_MAX is undefined!'
382 #endif
383 
384 #if ( defined(MARINE_ICE_BASAL_MELTING) \
385  && ( marine_ice_basal_melting==2 || marine_ice_basal_melting==3 ) )
386 
387 if (qbm_marine*year_sec_inv*rhow_rho_ratio > qbm_max) then
388  write (6,'(a)') ' >>> calc_qbm: QBM_MARINE'
389  write (6,'(a)') ' (basal melting rate at the ice front)'
390  write (6,'(a)') ' is larger than the limiter qbm_max!'
391  stop
392 end if
393 
394 #endif
395 
396 #if ( MARGIN==3 \
397  && defined(floating_ice_basal_melting) \
398  && floating_ice_basal_melting<=3 )
399 
400 if (qbm_float_1*year_sec_inv*rhow_rho_ratio > qbm_max) then
401  write (6,'(a)') ' >>> calc_qbm: QBM_FLOAT_1'
402  write (6,'(a)') ' (basal melting rate in the grounding line zone)'
403  write (6,'(a)') ' is larger than the limiter qbm_max!'
404  stop
405 end if
406 
407 #endif
408 
409 do i=0, imax
410 do j=0, jmax
411  if (maske(j,i)==0) then ! grounded ice
412  if (q_bm(j,i) < qbm_min) q_bm(j,i) = 0.0_dp
413  if (q_bm(j,i) > qbm_max) q_bm(j,i) = qbm_max
414  if (q_tld(j,i) < qbm_min) q_tld(j,i) = 0.0_dp
415  if (q_tld(j,i) > qbm_max) q_tld(j,i) = qbm_max
416  if (q_b_tot(j,i) < qbm_min) q_b_tot(j,i) = 0.0_dp
417  if (q_b_tot(j,i) > qbm_max) q_b_tot(j,i) = qbm_max
418  end if
419 end do
420 end do
421 
422 end subroutine calc_qbm
423 
424 !-------------------------------------------------------------------------------
425 !> Sub-ice-shelf melting parameterisation.
426 !<------------------------------------------------------------------------------
427 subroutine sub_ice_shelf_melting_param(time, z_sl, &
428  year_sec_inv, rhow_rho_ratio, &
429  i, j, q_bm_floating)
431 implicit none
432 
433 integer(i4b), intent(in) :: i, j
434 real(dp), intent(in) :: time, z_sl
435 real(dp), intent(in) :: year_sec_inv, rhow_rho_ratio
436 
437 real(dp), intent(out) :: Q_bm_floating
438 
439 real(dp) :: time_in_years
440 real(dp) :: Toc, Tmb, T_forcing, Omega, alpha, draft0, draft
441 real(dp) :: lon_d, lat_d, Phi_par
442 real(dp) :: H_w_now, Q_bm_scaling_factor
443 
444 real(dp), parameter :: beta_sw = 7.61e-04_dp
445 
446 #if (FLOATING_ICE_BASAL_MELTING==5)
447 
448 real(dp), parameter :: c_sw = 3974.0_dp, &
449  g_t = 5.0e-05_dp
450 
451 #endif
452 
453 time_in_years = time*year_sec_inv
454 
455 #if (FLOATING_ICE_BASAL_MELTING==4)
456 
457 toc = temp_ocean
458 omega = omega_qbm *year_sec_inv*rhow_rho_ratio
459  ! m/[a*degC^alpha] water equiv. -> m/[s*degC^alpha] ice equiv.
460 alpha = alpha_qbm
461 
462 #elif (FLOATING_ICE_BASAL_MELTING==5)
463 
464 #if (defined(ANT)) /* Antarctic ice sheet */
465 
466 !-------- Definition of the sectors --------
467 
468 lon_d = lambda(j,i) *pi_180_inv ! rad -> deg
469 lat_d = phi(j,i) *pi_180_inv ! rad -> deg
470 
471 lon_d = modulo(lon_d+180.0_dp, 360.0_dp)-180.0_dp
472  ! mapping to interval
473  ! [-180 deg, +180 deg)
474 
475 if ((lon_d>=-10.0_dp).and.(lon_d<60.0_dp)) then
476  ! Western East Antarctica
477  n_sector(j,i) = 1_i2b
478 
479  toc = 0.23_dp
480  omega = 0.0076_dp
481  alpha = 1.33_dp
482  draft0 = 200.0_dp
483 
484 else if ((lon_d>=60.0_dp).and.(lon_d<80.0_dp)) then
485  ! Amery/Prydz Bay
486  n_sector(j,i) = 2_i2b
487 
488  toc = -1.61_dp
489  omega = 0.0128_dp
490  alpha = 0.94_dp
491  draft0 = 200.0_dp
492 
493 else if ((lon_d>=80.0_dp).and.(lon_d<130.0_dp)) then
494  ! Sabrina Coast/
495  ! Aurora subglacial basin
496  n_sector(j,i) = 3_i2b
497 
498  toc = -0.16_dp
499  omega = 0.0497_dp
500  alpha = 0.85_dp
501  draft0 = 200.0_dp
502 
503 else if ( ((lon_d>=130.0_dp).and.(lon_d<159.0_dp)) &
504  .or. &
505  ((lon_d>=159.0_dp).and.(lon_d<170.0_dp) &
506  .and.(lat_d>=-72.0_dp)) ) then
507  ! George V Coast/
508  ! Wilkes subglacial basin
509  n_sector(j,i) = 4_i2b
510 
511  toc = -0.22_dp
512  omega = 0.0423_dp
513  alpha = 0.11_dp
514  draft0 = 200.0_dp
515 
516 else if ( (lon_d>=159.0_dp) &
517  .or. &
518  (lon_d<-140.0_dp) &
519  .or. &
520  ((lon_d>=-140.0_dp).and.(lon_d<-120.0_dp) &
521  .and.(lat_d<-77.0_dp)) ) then
522  ! Ross Sea
523  n_sector(j,i) = 5_i2b
524 
525  toc = -1.8_dp
526  omega = 0.0181_dp
527  alpha = 0.73_dp
528  draft0 = 200.0_dp
529 
530 else if ( ((lon_d>=-140.0_dp).and.(lon_d<-120.0_dp) &
531  .and.(lat_d>=-77.0_dp)) &
532  .or. &
533  ((lon_d>=-120.0_dp).and.(lon_d<-90.0_dp)) ) then
534  ! Amundsen Sea
535  n_sector(j,i) = 6_i2b
536 
537  toc = 0.42_dp
538  omega = 0.1023_dp
539  alpha = 0.26_dp
540  draft0 = 200.0_dp
541 
542 else if ( ((lon_d>=-90.0_dp).and.(lon_d<-66.0_dp) &
543  .and.(lat_d>=-74.0_dp)) ) then
544  ! Bellingshausen Sea
545  n_sector(j,i) = 7_i2b
546 
547  toc = 0.62_dp
548  omega = 0.0644_dp
549  alpha = 0.16_dp
550  draft0 = 200.0_dp
551 
552 else
553  ! Weddell Sea
554  n_sector(j,i) = 8_i2b
555 
556  toc = -1.55_dp
557  omega = 0.0197_dp
558  alpha = 0.26_dp
559  draft0 = 200.0_dp
560 
561 endif
562 
563 #else /* not Antarctic ice sheet */
564 
565 q_bm_floating = 0.0_dp ! dummy return value
566 
567 write(6, fmt='(a)') ' >>> sub_ice_shelf_melting_param:'
568 write(6, fmt='(a)') ' Case FLOATING_ICE_BASAL_MELTING==5'
569 write(6, fmt='(a)') ' only defined for Antarctica!'
570 stop
571 
572 #endif
573 
574 #else
575 
576 continue ! do nothing
577 
578 #endif
579 
580 !-------- Computation of the sub-ice-shelf melting rate --------
581 
582 if (maske(j,i)==2_i2b) then ! ocean
583  draft = 0.0_dp
584  h_w_now = max((z_sl-zl(j,i)), 0.0_dp)
585 else if (maske(j,i)==3_i2b) then ! floating ice
586  draft = max((z_sl -zb(j,i)), 0.0_dp)
587  h_w_now = max((zb(j,i)-zl(j,i)), 0.0_dp)
588 else
589  write(6, fmt='(a)') ' >>> sub_ice_shelf_melting_param:'
590  write(6, fmt='(a)') ' Routine must not be called for maske(j,i) < 2!'
591  stop
592 end if
593 
594 tmb = -beta_sw*draft - delta_tm_sw ! temperature at the ice shelf base
595 
596 t_forcing = max((toc-tmb), 0.0_dp) ! thermal forcing
597 
598 #if (defined(H_W_0))
599 
600 if (h_w_0 > eps) then
601  q_bm_scaling_factor = tanh(h_w_now/h_w_0)
602 else
603  q_bm_scaling_factor = 1.0_dp
604 end if
605 
606 #else
607 q_bm_scaling_factor = 1.0_dp
608 #endif
609 
610 #if (FLOATING_ICE_BASAL_MELTING==4)
611 
612 q_bm_floating = q_bm_scaling_factor*omega*t_forcing**alpha
613 
614 #elif (FLOATING_ICE_BASAL_MELTING==5)
615 
616 phi_par = rho_sw*c_sw*g_t/(rho*l)
617 
618 q_bm_floating = q_bm_scaling_factor &
619  *phi_par*omega*t_forcing*(draft/draft0)**alpha
620 
621 #else
622 
623 q_bm_floating = 0.0_dp ! dummy return value
624 
625 write(6, fmt='(a)') ' >>> sub_ice_shelf_melting_param:'
626 write(6, fmt='(a)') ' FLOATING_ICE_BASAL_MELTING must be 4 or 5!'
627 stop
628 
629 #endif
630 
631 #if (defined(INITMIP_BMB_ANOM_FILE))
632 
633 ! ------ Correction for ISMIP InitMIP
634 
635 if ((time_in_years > 0.0_dp).and.(time_in_years <= 40.0_dp)) then
636  q_bm_floating = q_bm_floating &
637  + 0.025_dp*floor(time_in_years) * ab_anom_initmip(j,i)
638 else if (time_in_years > 40.0_dp) then
639  q_bm_floating = q_bm_floating + ab_anom_initmip(j,i)
640 end if
641 
642 #endif
643 
644 end subroutine sub_ice_shelf_melting_param
645 
646 !-------------------------------------------------------------------------------
647 
648 end module calc_bas_melt_m
649 !
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.
real(dp), dimension(0:jmax, 0:imax) phi
phi(j,i): Geographic latitude of grid point (i,j)
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...
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...
integer(i2b), dimension(0:jmax, 0:imax) n_sector
n_sector(j,i): Marker for the different sectors for ice shelf basal melting.
Definition: sico_vars_m.F90:55
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.
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(.).
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) 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).
Definition: sico_vars_m.F90:35
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)
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) 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.
real(dp) rho
RHO: Density of ice.
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)) ...