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