SICOPOLIS V5-dev  Revision 1279
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-2018 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 !<------------------------------------------------------------------------------
36 
37  use sico_types_m
39  use sico_vars_m
40  use error_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
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  errormsg = ' >>> calc_qbm: Marine ice margin point does not' &
208  // end_of_line &
209  //' have an ocean neighbour!'
210  call error(errormsg)
211  end if
212 
213  end if
214 
215 #else
216  errormsg = ' >>> calc_qbm: MARINE_ICE_BASAL_MELTING must be 1, 2 or 3!'
217  call error(errormsg)
218 #endif
219 
220 #elif (MARGIN==3)
221 
222  if (flag_grounding_line_1(j,i)) then
223  ! grounding line (grounded-ice side)
224 
225 #if (FLOATING_ICE_BASAL_MELTING==1 || FLOATING_ICE_BASAL_MELTING==4 \
226  || floating_ice_basal_melting==5)
227 
228  !!! continue
229 
230 #elif (FLOATING_ICE_BASAL_MELTING==2)
231 
232  q_bm(j,i) = qbm_float_1 *year_sec_inv*rhow_rho_ratio
233  ! m/a water equiv. -> m/s ice equiv.
234 
235 #elif (FLOATING_ICE_BASAL_MELTING==3)
236 
237  n_float = 0
238  if (maske(j,i+1)>=2) n_float = n_float+1
239  if (maske(j,i-1)>=2) n_float = n_float+1
240  if (maske(j+1,i)>=2) n_float = n_float+1
241  if (maske(j-1,i)>=2) n_float = n_float+1
242 
243  if ( n_float > 0 ) then
244 
245  q_bm_grounded = q_bm(j,i)
246  q_bm_floating = qbm_float_1 *year_sec_inv*rhow_rho_ratio
247  ! m/a water equiv. -> m/s ice equiv.
248 
249  q_bm(j,i) = (1.0_dp-0.25_dp*real(n_float,dp)) * q_bm_grounded &
250  +0.25_dp*real(n_float,dp) * Q_bm_floating
251  ! weighed average of melting for grounded and floating ice
252  else
253  errormsg = ' >>> calc_qbm: Grounding line point does not have' &
254  // end_of_line &
255  //' a floating ice or ocean neighbour!'
256  call error(errormsg)
257  end if
258 
259 #else
260  errormsg = ' >>> calc_qbm: FLOATING_ICE_BASAL_MELTING' &
261  // end_of_line &
262  //' must be 1, 2, 3, 4 or 5!'
263  call error(errormsg)
264 #endif
265 
266  else if ( (zb(j,i) < z_sl) & ! marine ice margin
267  .and. &
268  ( (maske(j,i+1)>=2) & ! (at least one
269  .or.(maske(j,i-1)>=2) & ! nearest neighbour
270  .or.(maske(j+1,i)>=2) & ! is
271  .or.(maske(j-1,i)>=2) & ! ocean)
272  ) &
273  ) then
274 
275 #if (MARINE_ICE_BASAL_MELTING==1)
276 
277  !!! continue
278 
279 #elif (MARINE_ICE_BASAL_MELTING==2)
280 
281  q_bm(j,i) = qbm_marine *year_sec_inv*rhow_rho_ratio
282  ! m/a water equiv. -> m/s ice equiv.
283 
284 #elif (MARINE_ICE_BASAL_MELTING==3)
285 
286  n_ocean = 0
287  if (maske(j,i+1)>=2) n_ocean = n_ocean+1
288  if (maske(j,i-1)>=2) n_ocean = n_ocean+1
289  if (maske(j+1,i)>=2) n_ocean = n_ocean+1
290  if (maske(j-1,i)>=2) n_ocean = n_ocean+1
291 
292  if ( n_ocean > 0 ) then
293 
294  q_bm_grounded = q_bm(j,i)
295  q_bm_marine = qbm_marine *year_sec_inv*rhow_rho_ratio
296  ! m/a water equiv. -> m/s ice equiv.
297 
298  q_bm(j,i) = (1.0_dp-0.25_dp*real(n_ocean,dp)) * q_bm_grounded &
299  +0.25_dp*real(n_ocean,dp) * Q_bm_marine
300  ! weighed average of grounded ice melting (computed)
301  ! and marine ice melting (prescribed)
302  else
303  errormsg = ' >>> calc_qbm: Marine ice margin point does not' &
304  // end_of_line &
305  //' have a floating ice or ocean neighbour!'
306  call error(errormsg)
307  end if
308 
309 #else
310  errormsg = ' >>> calc_qbm: MARINE_ICE_BASAL_MELTING must be 1, 2 or 3!'
311  call error(errormsg)
312 #endif
313 
314  end if
315 
316  else if ( (maske(j,i)==2).or.(maske(j,i)==3) ) then ! floating ice or ocean
317 
318 #if (FLOATING_ICE_BASAL_MELTING<=3)
319 
320  if (flag_grounding_line_2(j,i)) then
321  ! grounding line (floating-ice side)
322 
323  q_bm(j,i) = qbm_float_1 *year_sec_inv*rhow_rho_ratio
324  ! m/a water equiv. -> m/s ice equiv.
325 
326  else if ( zl(j,i) > z_abyss ) then ! floating ice over continental shelf
327 
328  q_bm(j,i) = qbm_float_2 *year_sec_inv*rhow_rho_ratio
329  ! m/a water equiv. -> m/s ice equiv.
330 
331  else ! ( zl(j,i) <= Z_ABYSS ), floating ice over abyssal ocean
332 
333  q_bm(j,i) = qbm_float_3 *year_sec_inv*rhow_rho_ratio
334  ! m/a water equiv. -> m/s ice equiv.
335 
336  end if
337 
338 #elif (FLOATING_ICE_BASAL_MELTING==4 || FLOATING_ICE_BASAL_MELTING==5)
339 
340  if ( zl(j,i) > z_abyss ) then ! floating ice over continental shelf
341  call sub_ice_shelf_melting_param(time, z_sl, &
342  year_sec_inv, rhow_rho_ratio, &
343  i, j, q_bm_floating)
344  q_bm(j,i) = q_bm_floating
345  else ! ( zl(j,i) <= Z_ABYSS ), floating ice over abyssal ocean
346  q_bm(j,i) = qbm_float_3 *year_sec_inv*rhow_rho_ratio
347  ! m/a water equiv. -> m/s ice equiv.
348  end if
349 
350 #else
351  errormsg = ' >>> calc_qbm: FLOATING_ICE_BASAL_MELTING' &
352  // end_of_line &
353  //' must be 1, 2, 3, 4 or 5!'
354  call error(errormsg)
355 #endif
356 
357 #else
358  errormsg = ' >>> calc_qbm: MARGIN must be 1, 2 or 3!'
359  call error(errormsg)
360 #endif
361 
362  end if
363 
364 end do
365 end do
366 
367 !-------- Sum of Q_bm and Q_tld --------
368 
369 q_b_tot = q_bm + q_tld
370 
371 !-------- Limitation of Q_bm, Q_tld and Q_b_tot
372 ! (only for grounded ice) --------
373 
374 #if (defined(QBM_MIN))
375  qbm_min = qbm_min *year_sec_inv*rhow_rho_ratio
376  ! m/a water equiv. -> m/s ice equiv.
377 #elif (defined(QB_MIN)) /* obsolete */
378  qbm_min = qb_min ! m/s ice equiv.
379 #else
380  errormsg = ' >>> calc_qbm: Limiter QBM_MIN is undefined!'
381  call error(errormsg)
382 #endif
383 
384 #if (defined(QBM_MAX))
385  qbm_max = qbm_max *year_sec_inv*rhow_rho_ratio
386  ! m/a water equiv. -> m/s ice equiv.
387 #elif (defined(QB_MAX)) /* obsolete */
388  qbm_max = qb_max ! m/s ice equiv.
389 #else
390  errormsg = ' >>> calc_qbm: Limiter QBM_MAX is undefined!'
391  call error(errormsg)
392 #endif
393 
394 #if ( defined(MARINE_ICE_BASAL_MELTING) \
395  && ( marine_ice_basal_melting==2 || marine_ice_basal_melting==3 ) )
396 
397 if (qbm_marine*year_sec_inv*rhow_rho_ratio > qbm_max) then
398  errormsg = ' >>> calc_qbm: QBM_MARINE' &
399  // end_of_line &
400  //' (basal melting rate at the ice front)' &
401  // end_of_line &
402  //' is larger than the limiter qbm_max!'
403  call error(errormsg)
404 end if
405 
406 #endif
407 
408 #if ( MARGIN==3 \
409  && defined(floating_ice_basal_melting) \
410  && floating_ice_basal_melting<=3 )
411 
412 if (qbm_float_1*year_sec_inv*rhow_rho_ratio > qbm_max) then
413  errormsg = ' >>> calc_qbm: QBM_FLOAT_1' &
414  // end_of_line &
415  //' (basal melting rate in the grounding line zone)' &
416  // end_of_line &
417  //' is larger than the limiter qbm_max!'
418  call error(errormsg)
419 end if
420 
421 #endif
422 
423 do i=0, imax
424 do j=0, jmax
425  if (maske(j,i)==0) then ! grounded ice
426  if (q_bm(j,i) < qbm_min) q_bm(j,i) = 0.0_dp
427  if (q_bm(j,i) > qbm_max) q_bm(j,i) = qbm_max
428  if (q_tld(j,i) < qbm_min) q_tld(j,i) = 0.0_dp
429  if (q_tld(j,i) > qbm_max) q_tld(j,i) = qbm_max
430  if (q_b_tot(j,i) < qbm_min) q_b_tot(j,i) = 0.0_dp
431  if (q_b_tot(j,i) > qbm_max) q_b_tot(j,i) = qbm_max
432  end if
433 end do
434 end do
435 
436 end subroutine calc_qbm
437 
438 !-------------------------------------------------------------------------------
439 !> Sub-ice-shelf melting parameterisation.
440 !<------------------------------------------------------------------------------
441 subroutine sub_ice_shelf_melting_param(time, z_sl, &
442  year_sec_inv, rhow_rho_ratio, &
443  i, j, q_bm_floating)
445 #ifdef ALLOW_OPENAD /* OpenAD */
446  use ctrl_m, only: myfloor
447 #endif /* OpenAD */
448 
449 implicit none
450 
451 #ifndef ALLOW_OPENAD /* Normal */
452 integer(i4b), intent(in) :: i, j
453 #else /* OpenAD */
454 integer(i4b), intent(inout) :: i, j
455 #endif /* Normal vs. OpenAD */
456 real(dp), intent(in) :: time, z_sl
457 real(dp), intent(in) :: year_sec_inv, rhow_rho_ratio
458 
459 real(dp), intent(out) :: Q_bm_floating
460 
461 integer(i4b) :: n
462 real(dp) :: time_in_years
463 real(dp) :: Toc, Tmb, T_forcing, Omega, alpha, draft0, draft
464 real(dp) :: lon_d, lat_d, Phi_par
465 real(dp) :: H_w_now, Q_bm_scaling_factor
466 
467 real(dp), parameter :: beta_sw = 7.61e-04_dp
468 
469 #if (FLOATING_ICE_BASAL_MELTING==5)
470 
471 real(dp), parameter :: c_sw = 3974.0_dp, &
472  g_t = 5.0e-05_dp
473 
474 #endif
475 
476 #ifdef ALLOW_OPENAD /* OpenAD */
477 integer(i4b) :: i_time_in_years
478 #endif /* OpenAD */
479 
480 time_in_years = time*year_sec_inv
481 
482 #if (FLOATING_ICE_BASAL_MELTING==4)
483 
484 toc = temp_ocean
485 omega = omega_qbm *year_sec_inv*rhow_rho_ratio
486  ! m/[a*degC^alpha] water equiv. -> m/[s*degC^alpha] ice equiv.
487 alpha = alpha_qbm
488 
489 #elif (FLOATING_ICE_BASAL_MELTING==5)
490 
491 #if (defined(ANT)) /* Antarctic ice sheet */
492 
493 !-------- Definition of the sectors --------
494 
495 lon_d = lambda(j,i) *pi_180_inv ! rad -> deg
496 lat_d = phi(j,i) *pi_180_inv ! rad -> deg
497 
498 #ifndef ALLOW_OPENAD /* Normal */
499 lon_d = modulo(lon_d+180.0_dp, 360.0_dp)-180.0_dp
500  ! mapping to interval
501  ! [-180 deg, +180 deg)
502 #else /* OpenAD */
503 lon_d = (lon_d+180.0_sp) - &
504  ((360.0_sp * int(lon_d+180.0_sp)/360.0_sp)) - &
505  180.0_sp
506 #endif /* Normal vs. OpenAD */
507 
508 if ((lon_d>=-10.0_dp).and.(lon_d<60.0_dp)) then
509  ! Western East Antarctica
510  n_sector(j,i) = 1_i2b
511 
512  toc = 0.23_dp
513  omega = 0.0076_dp
514  alpha = 1.33_dp
515  draft0 = 200.0_dp
516 
517 else if ((lon_d>=60.0_dp).and.(lon_d<80.0_dp)) then
518  ! Amery/Prydz Bay
519  n_sector(j,i) = 2_i2b
520 
521  toc = -1.61_dp
522  omega = 0.0128_dp
523  alpha = 0.94_dp
524  draft0 = 200.0_dp
525 
526 else if ((lon_d>=80.0_dp).and.(lon_d<130.0_dp)) then
527  ! Sabrina Coast/
528  ! Aurora subglacial basin
529  n_sector(j,i) = 3_i2b
530 
531  toc = -0.16_dp
532  omega = 0.0497_dp
533  alpha = 0.85_dp
534  draft0 = 200.0_dp
535 
536 else if ( ((lon_d>=130.0_dp).and.(lon_d<159.0_dp)) &
537  .or. &
538  ((lon_d>=159.0_dp).and.(lon_d<170.0_dp) &
539  .and.(lat_d>=-72.0_dp)) ) then
540  ! George V Coast/
541  ! Wilkes subglacial basin
542  n_sector(j,i) = 4_i2b
543 
544  toc = -0.22_dp
545  omega = 0.0423_dp
546  alpha = 0.11_dp
547  draft0 = 200.0_dp
548 
549 else if ( (lon_d>=159.0_dp) &
550  .or. &
551  (lon_d<-140.0_dp) &
552  .or. &
553  ((lon_d>=-140.0_dp).and.(lon_d<-120.0_dp) &
554  .and.(lat_d<-77.0_dp)) ) then
555  ! Ross Sea
556  n_sector(j,i) = 5_i2b
557 
558  toc = -1.8_dp
559  omega = 0.0181_dp
560  alpha = 0.73_dp
561  draft0 = 200.0_dp
562 
563 else if ( ((lon_d>=-140.0_dp).and.(lon_d<-120.0_dp) &
564  .and.(lat_d>=-77.0_dp)) &
565  .or. &
566  ((lon_d>=-120.0_dp).and.(lon_d<-90.0_dp)) ) then
567  ! Amundsen Sea
568  n_sector(j,i) = 6_i2b
569 
570  toc = 0.42_dp
571  omega = 0.1023_dp
572  alpha = 0.26_dp
573  draft0 = 200.0_dp
574 
575 else if ( ((lon_d>=-90.0_dp).and.(lon_d<-66.0_dp) &
576  .and.(lat_d>=-74.0_dp)) ) then
577  ! Bellingshausen Sea
578  n_sector(j,i) = 7_i2b
579 
580  toc = 0.62_dp
581  omega = 0.0644_dp
582  alpha = 0.16_dp
583  draft0 = 200.0_dp
584 
585 else
586  ! Weddell Sea
587  n_sector(j,i) = 8_i2b
588 
589  toc = -1.55_dp
590  omega = 0.0197_dp
591  alpha = 0.26_dp
592  draft0 = 200.0_dp
593 
594 endif
595 
596 #else /* not Antarctic ice sheet */
597 
598 q_bm_floating = 0.0_dp ! dummy return value
599 
600 errormsg = ' >>> sub_ice_shelf_melting_param:' &
601  // end_of_line &
602  //' Case FLOATING_ICE_BASAL_MELTING==5' &
603  // end_of_line &
604  //' only defined for Antarctica!'
605 call error(errormsg)
606 
607 #endif
608 
609 #else
610 
611 !!! continue
612 
613 #endif
614 
615 !-------- Computation of the sub-ice-shelf melting rate --------
616 
617 if (maske(j,i)==2_i2b) then ! ocean
618  draft = 0.0_dp
619  h_w_now = max((z_sl-zl(j,i)), 0.0_dp)
620 else if (maske(j,i)==3_i2b) then ! floating ice
621  draft = max((z_sl -zb(j,i)), 0.0_dp)
622  h_w_now = max((zb(j,i)-zl(j,i)), 0.0_dp)
623 else
624  errormsg = ' >>> sub_ice_shelf_melting_param:' &
625  // end_of_line &
626  //' Routine must not be called for maske(j,i) < 2!'
627  call error(errormsg)
628 end if
629 
630 tmb = -beta_sw*draft - delta_tm_sw ! temperature at the ice shelf base
631 
632 t_forcing = max((toc-tmb), 0.0_dp) ! thermal forcing
633 
634 #if (defined(H_W_0))
635 
636 if (h_w_0 > eps) then
637  q_bm_scaling_factor = tanh(h_w_now/h_w_0)
638 else
639  q_bm_scaling_factor = 1.0_dp
640 end if
641 
642 #else
643 q_bm_scaling_factor = 1.0_dp
644 #endif
645 
646 #if (FLOATING_ICE_BASAL_MELTING==4)
647 
648 q_bm_floating = q_bm_scaling_factor*omega*t_forcing**alpha
649 
650 #elif (FLOATING_ICE_BASAL_MELTING==5)
651 
652 phi_par = rho_sw*c_sw*g_t/(rho*l)
653 
654 q_bm_floating = q_bm_scaling_factor &
655  *phi_par*omega*t_forcing*(draft/draft0)**alpha
656 
657 #else
658 
659 q_bm_floating = 0.0_dp ! dummy return value
660 
661 errormsg = ' >>> sub_ice_shelf_melting_param:' &
662  // end_of_line &
663  //' FLOATING_ICE_BASAL_MELTING must be 4 or 5!'
664 call error(errormsg)
665 
666 #endif
667 
668 ! ------ Correction for ISMIP InitMIP
669 
670 #if (defined(INITMIP_BMB_ANOM_FILE))
671 
672 if ((time_in_years > 0.0_dp).and.(time_in_years <= 40.0_dp)) then
673 
674  q_bm_floating = q_bm_floating &
675 #ifndef ALLOW_OPENAD /* Normal */
676  + 0.025_dp*floor(time_in_years) * ab_anom_initmip(j,i)
677 #else /* OpenAD */
678 call myfloor(time_in_years,i_time_in_years)
679  + 0.025_dp*i_time_in_years * ab_anom_initmip(j,i)
680 #endif /* Normal vs. OpenAD */
681 
682 else if (time_in_years > 40.0_dp) then
683 
684  q_bm_floating = q_bm_floating + ab_anom_initmip(j,i)
685 
686 end if
687 
688 #endif
689 
690 ! ------ Correction for ISMIP LARMIP
691 
692 #if (defined(LARMIP_REGIONS_FILE))
693  n = int(n_regions_larmip(j,i))
694  q_bm_floating = q_bm_floating + ab_anom_larmip(n)
695 #endif
696 
697 end subroutine sub_ice_shelf_melting_param
698 
699 !-------------------------------------------------------------------------------
700 
701 end module calc_bas_melt_m
702 !
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:67
Definition: ctrl_m.F90:9
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.
Definition: error_m.F90:47
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)
Writing of error messages and stopping execution.
Definition: error_m.F90:35
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
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, 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)) ...