SICOPOLIS V5-dev  Revision 1420
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-2019 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)*0.5_dp &
107  * ( (vx_t(0,j,i)+vx_t(0,j,i-1))*dzs_dxi_g(j,i) &
108  +(vy_t(0,j,i)+vy_t(0,j-1,i))*dzs_deta_g(j,i) )
109 #if (DYNAMICS==2)
110  else ! flag_shelfy_stream(j,i) == .true.
111 
112  frictional_heating &
113  = aqbm3b &
114  * c_drag(j,i) &
115  * sqrt(vx_b_g(j,i)**2 &
116  +vy_b_g(j,i)**2) &
117  **(1.0_dp+p_weert_inv(j,i))
118  end if
119 #endif
120  q_bm(j,i) = aqbm1*(temp_c(1,j,i)-temp_c(0,j,i))/h_c(j,i) &
121  *kappa_val(temp_c(0,j,i)) &
122  - aqbm2*(temp_r(krmax,j,i)-temp_r(krmax-1,j,i)) &
123  + frictional_heating
124 
125  else ! n_cts(j,i)==1_i1b
126 
127 #if (DYNAMICS==2)
128  if (.not.flag_shelfy_stream(j,i)) then
129 #endif
130  frictional_heating &
131  = -aqbm3a*(h_c(j,i)+h_t(j,i))*0.5_dp &
132  * ( (vx_t(0,j,i)+vx_t(0,j,i-1))*dzs_dxi_g(j,i) &
133  +(vy_t(0,j,i)+vy_t(0,j-1,i))*dzs_deta_g(j,i) )
134 #if (DYNAMICS==2)
135  else ! flag_shelfy_stream(j,i) == .true.
136 
137  frictional_heating &
138  = aqbm3b &
139  * c_drag(j,i) &
140  * sqrt(vx_b_g(j,i)**2 &
141  +vy_b_g(j,i)**2) &
142  **(1.0_dp+p_weert_inv(j,i))
143  end if
144 #endif
145  q_bm(j,i) = aqbm4*kappa_val(temp_t_m(0,j,i)) &
146  - aqbm2*(temp_r(krmax,j,i)-temp_r(krmax-1,j,i)) &
147  + frictional_heating
148 
149  end if
150 
151 #if (MARGIN==1 || MARGIN==2)
152 
153 #if (!defined(MARINE_ICE_BASAL_MELTING) || MARINE_ICE_BASAL_MELTING==1)
154 
155  !!! continue
156 
157 #elif (MARINE_ICE_BASAL_MELTING==2)
158 
159  if ( (zb(j,i) < z_sl) & ! marine ice
160  .and. &
161  ( (maske(j,i+1)==2_i1b) & ! at least one
162  .or.(maske(j,i-1)==2_i1b) & ! nearest neighbour
163  .or.(maske(j+1,i)==2_i1b) & ! is
164  .or.(maske(j-1,i)==2_i1b) & ! ocean
165  ) &
166  ) then
167 
168  q_bm(j,i) = qbm_marine *year_sec_inv*rhow_rho_ratio
169  ! m/a water equiv. -> m/s ice equiv.
170 
171  end if
172 
173 #elif (MARINE_ICE_BASAL_MELTING==3)
174 
175  if ( (zb(j,i) < z_sl) & ! marine ice
176  .and. &
177  ( (maske(j,i+1)==2_i1b) & ! at least one
178  .or.(maske(j,i-1)==2_i1b) & ! nearest neighbour
179  .or.(maske(j+1,i)==2_i1b) & ! is
180  .or.(maske(j-1,i)==2_i1b) & ! ocean
181  ) &
182  ) then
183 
184  n_ocean = 0
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
189 
190  if ( n_ocean > 0 ) then
191 
192  q_bm_grounded = q_bm(j,i)
193  q_bm_marine = qbm_marine *year_sec_inv*rhow_rho_ratio
194  ! m/a water equiv. -> m/s ice equiv.
195 
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
198  ! weighed average of grounded ice melting (computed)
199  ! and marine ice melting (prescribed)
200  else
201  errormsg = ' >>> calc_qbm: Marine ice margin point does not' &
202  // end_of_line &
203  //' have an ocean neighbour!'
204  call error(errormsg)
205  end if
206 
207  end if
208 
209 #else
210  errormsg = ' >>> calc_qbm: MARINE_ICE_BASAL_MELTING must be 1, 2 or 3!'
211  call error(errormsg)
212 #endif
213 
214 #elif (MARGIN==3)
215 
216  if (flag_grounding_line_1(j,i)) then
217  ! grounding line (grounded-ice side)
218 
219 #if (FLOATING_ICE_BASAL_MELTING==1 || FLOATING_ICE_BASAL_MELTING==4 \
220  || floating_ice_basal_melting==5)
221 
222  !!! continue
223 
224 #elif (FLOATING_ICE_BASAL_MELTING==2)
225 
226  q_bm(j,i) = qbm_float_1 *year_sec_inv*rhow_rho_ratio
227  ! m/a water equiv. -> m/s ice equiv.
228 
229 #elif (FLOATING_ICE_BASAL_MELTING==3)
230 
231  n_float = 0
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
236 
237  if ( n_float > 0 ) then
238 
239  q_bm_grounded = q_bm(j,i)
240  q_bm_floating = qbm_float_1 *year_sec_inv*rhow_rho_ratio
241  ! m/a water equiv. -> m/s ice equiv.
242 
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
245  ! weighed average of melting for grounded and floating ice
246  else
247  errormsg = ' >>> calc_qbm: Grounding line point does not have' &
248  // end_of_line &
249  //' a floating ice or ocean neighbour!'
250  call error(errormsg)
251  end if
252 
253 #else
254  errormsg = ' >>> calc_qbm: FLOATING_ICE_BASAL_MELTING' &
255  // end_of_line &
256  //' must be 1, 2, 3, 4 or 5!'
257  call error(errormsg)
258 #endif
259 
260  else if ( (zb(j,i) < z_sl) & ! marine ice margin
261  .and. &
262  ( (maske(j,i+1)>=2_i1b) & ! (at least one
263  .or.(maske(j,i-1)>=2_i1b) & ! nearest neighbour
264  .or.(maske(j+1,i)>=2_i1b) & ! is
265  .or.(maske(j-1,i)>=2_i1b) & ! ocean)
266  ) &
267  ) then
268 
269 #if (MARINE_ICE_BASAL_MELTING==1)
270 
271  !!! continue
272 
273 #elif (MARINE_ICE_BASAL_MELTING==2)
274 
275  q_bm(j,i) = qbm_marine *year_sec_inv*rhow_rho_ratio
276  ! m/a water equiv. -> m/s ice equiv.
277 
278 #elif (MARINE_ICE_BASAL_MELTING==3)
279 
280  n_ocean = 0
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
285 
286  if ( n_ocean > 0 ) then
287 
288  q_bm_grounded = q_bm(j,i)
289  q_bm_marine = qbm_marine *year_sec_inv*rhow_rho_ratio
290  ! m/a water equiv. -> m/s ice equiv.
291 
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
294  ! weighed average of grounded ice melting (computed)
295  ! and marine ice melting (prescribed)
296  else
297  errormsg = ' >>> calc_qbm: Marine ice margin point does not' &
298  // end_of_line &
299  //' have a floating ice or ocean neighbour!'
300  call error(errormsg)
301  end if
302 
303 #else
304  errormsg = ' >>> calc_qbm: MARINE_ICE_BASAL_MELTING must be 1, 2 or 3!'
305  call error(errormsg)
306 #endif
307 
308  end if
309 
310  else if ( (maske(j,i)==2_i1b).or.(maske(j,i)==3_i1b) ) then
311  ! floating ice or ocean
312 
313 #if (FLOATING_ICE_BASAL_MELTING<=3)
314 
315  if (flag_grounding_line_2(j,i)) then
316  ! grounding line (floating-ice side)
317 
318  q_bm(j,i) = qbm_float_1 *year_sec_inv*rhow_rho_ratio
319  ! m/a water equiv. -> m/s ice equiv.
320 
321  else if ( zl(j,i) > z_abyss ) then ! floating ice over continental shelf
322 
323  q_bm(j,i) = qbm_float_2 *year_sec_inv*rhow_rho_ratio
324  ! m/a water equiv. -> m/s ice equiv.
325 
326  else ! ( zl(j,i) <= Z_ABYSS ), floating ice over abyssal ocean
327 
328  q_bm(j,i) = qbm_float_3 *year_sec_inv*rhow_rho_ratio
329  ! m/a water equiv. -> m/s ice equiv.
330 
331  end if
332 
333 #elif (FLOATING_ICE_BASAL_MELTING==4 || FLOATING_ICE_BASAL_MELTING==5)
334 
335  if ( zl(j,i) > z_abyss ) then ! floating ice over continental shelf
336  call sub_ice_shelf_melting_param(time, z_sl, &
337  year_sec_inv, rhow_rho_ratio, &
338  i, j, q_bm_floating)
339  q_bm(j,i) = q_bm_floating
340  else ! ( zl(j,i) <= Z_ABYSS ), floating ice over abyssal ocean
341  q_bm(j,i) = qbm_float_3 *year_sec_inv*rhow_rho_ratio
342  ! m/a water equiv. -> m/s ice equiv.
343  end if
344 
345 #else
346  errormsg = ' >>> calc_qbm: FLOATING_ICE_BASAL_MELTING' &
347  // end_of_line &
348  //' must be 1, 2, 3, 4 or 5!'
349  call error(errormsg)
350 #endif
351 
352 #else
353  errormsg = ' >>> calc_qbm: MARGIN must be 1, 2 or 3!'
354  call error(errormsg)
355 #endif
356 
357  end if
358 
359 end do
360 end do
361 
362 !-------- Sum of Q_bm and Q_tld --------
363 
364 q_b_tot = q_bm + q_tld
365 
366 !-------- Limitation of Q_bm, Q_tld and Q_b_tot
367 ! (only for grounded ice) --------
368 
369 #if (defined(QBM_MIN))
370  qbm_min = qbm_min *year_sec_inv*rhow_rho_ratio
371  ! m/a water equiv. -> m/s ice equiv.
372 #elif (defined(QB_MIN)) /* obsolete */
373  qbm_min = qb_min ! m/s ice equiv.
374 #else
375  errormsg = ' >>> calc_qbm: Limiter QBM_MIN is undefined!'
376  call error(errormsg)
377 #endif
378 
379 #if (defined(QBM_MAX))
380  qbm_max = qbm_max *year_sec_inv*rhow_rho_ratio
381  ! m/a water equiv. -> m/s ice equiv.
382 #elif (defined(QB_MAX)) /* obsolete */
383  qbm_max = qb_max ! m/s ice equiv.
384 #else
385  errormsg = ' >>> calc_qbm: Limiter QBM_MAX is undefined!'
386  call error(errormsg)
387 #endif
388 
389 #if ( defined(MARINE_ICE_BASAL_MELTING) \
390  && ( marine_ice_basal_melting==2 || marine_ice_basal_melting==3 ) )
391 
392 if (qbm_marine*year_sec_inv*rhow_rho_ratio > qbm_max) then
393  errormsg = ' >>> calc_qbm: QBM_MARINE' &
394  // end_of_line &
395  //' (basal melting rate at the ice front)' &
396  // end_of_line &
397  //' is larger than the limiter qbm_max!'
398  call error(errormsg)
399 end if
400 
401 #endif
402 
403 #if ( MARGIN==3 \
404  && defined(floating_ice_basal_melting) \
405  && floating_ice_basal_melting<=3 )
406 
407 if (qbm_float_1*year_sec_inv*rhow_rho_ratio > qbm_max) then
408  errormsg = ' >>> calc_qbm: QBM_FLOAT_1' &
409  // end_of_line &
410  //' (basal melting rate in the grounding line zone)' &
411  // end_of_line &
412  //' is larger than the limiter qbm_max!'
413  call error(errormsg)
414 end if
415 
416 #endif
417 
418 do i=0, imax
419 do j=0, jmax
420  if (maske(j,i)==0_i1b) then ! grounded ice
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
423  if (q_tld(j,i) < qbm_min) q_tld(j,i) = 0.0_dp
424  if (q_tld(j,i) > qbm_max) q_tld(j,i) = qbm_max
425  if (q_b_tot(j,i) < qbm_min) q_b_tot(j,i) = 0.0_dp
426  if (q_b_tot(j,i) > qbm_max) q_b_tot(j,i) = qbm_max
427  end if
428 end do
429 end do
430 
431 end subroutine calc_qbm
432 
433 !-------------------------------------------------------------------------------
434 !> Sub-ice-shelf melting parameterisation.
435 !<------------------------------------------------------------------------------
436 subroutine sub_ice_shelf_melting_param(time, z_sl, &
437  year_sec_inv, rhow_rho_ratio, &
438  i, j, q_bm_floating)
440 #if defined(ALLOW_OPENAD) /* OpenAD */
441  use ctrl_m, only: myfloor
442 #endif /* OpenAD */
443 
444 implicit none
445 
446 #if !defined(ALLOW_OPENAD) /* Normal */
447 integer(i4b), intent(in) :: i, j
448 #else /* OpenAD */
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
453 
454 real(dp), intent(out) :: Q_bm_floating
455 
456 integer(i4b) :: n
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
461 
462 real(dp), parameter :: beta_sw = 7.61e-04_dp
463 
464 #if (FLOATING_ICE_BASAL_MELTING==5)
465 
466 real(dp), parameter :: c_sw = 3974.0_dp, &
467  g_t = 5.0e-05_dp
468 
469 #endif
470 
471 #if defined(ALLOW_OPENAD) /* OpenAD */
472 integer(i4b) :: i_time_in_years
473 #endif /* OpenAD */
474 
475 time_in_years = time*year_sec_inv
476 
477 #if (FLOATING_ICE_BASAL_MELTING==4)
478 
479 toc = temp_ocean
480 omega = omega_qbm *year_sec_inv*rhow_rho_ratio
481  ! m/[a*degC^alpha] water equiv. -> m/[s*degC^alpha] ice equiv.
482 alpha = alpha_qbm
483 
484 #elif (FLOATING_ICE_BASAL_MELTING==5)
485 
486 #if (defined(ANT)) /* Antarctic ice sheet */
487 
488 !-------- Definition of the sectors --------
489 
490 lon_d = lambda(j,i) *pi_180_inv ! rad -> deg
491 lat_d = phi(j,i) *pi_180_inv ! rad -> deg
492 
493 #if !defined(ALLOW_OPENAD) /* Normal */
494 lon_d = modulo(lon_d+180.0_dp, 360.0_dp)-180.0_dp
495  ! mapping to interval
496  ! [-180 deg, +180 deg)
497 #else /* OpenAD */
498 lon_d = (lon_d+180.0_sp) - &
499  ((360.0_sp * int(lon_d+180.0_sp)/360.0_sp)) - &
500  180.0_sp
501 #endif /* Normal vs. OpenAD */
502 
503 if ((lon_d>=-10.0_dp).and.(lon_d<60.0_dp)) then
504  ! Western East Antarctica
505  n_sector(j,i) = 1_i1b
506 
507  toc = 0.23_dp
508  omega = 0.0076_dp
509  alpha = 1.33_dp
510  draft0 = 200.0_dp
511 
512 else if ((lon_d>=60.0_dp).and.(lon_d<80.0_dp)) then
513  ! Amery/Prydz Bay
514  n_sector(j,i) = 2_i1b
515 
516  toc = -1.61_dp
517  omega = 0.0128_dp
518  alpha = 0.94_dp
519  draft0 = 200.0_dp
520 
521 else if ((lon_d>=80.0_dp).and.(lon_d<130.0_dp)) then
522  ! Sabrina Coast/
523  ! Aurora subglacial basin
524  n_sector(j,i) = 3_i1b
525 
526  toc = -0.16_dp
527  omega = 0.0497_dp
528  alpha = 0.85_dp
529  draft0 = 200.0_dp
530 
531 else if ( ((lon_d>=130.0_dp).and.(lon_d<159.0_dp)) &
532  .or. &
533  ((lon_d>=159.0_dp).and.(lon_d<170.0_dp) &
534  .and.(lat_d>=-72.0_dp)) ) then
535  ! George V Coast/
536  ! Wilkes subglacial basin
537  n_sector(j,i) = 4_i1b
538 
539  toc = -0.22_dp
540  omega = 0.0423_dp
541  alpha = 0.11_dp
542  draft0 = 200.0_dp
543 
544 else if ( (lon_d>=159.0_dp) &
545  .or. &
546  (lon_d<-140.0_dp) &
547  .or. &
548  ((lon_d>=-140.0_dp).and.(lon_d<-120.0_dp) &
549  .and.(lat_d<-77.0_dp)) ) then
550  ! Ross Sea
551  n_sector(j,i) = 5_i1b
552 
553  toc = -1.8_dp
554  omega = 0.0181_dp
555  alpha = 0.73_dp
556  draft0 = 200.0_dp
557 
558 else if ( ((lon_d>=-140.0_dp).and.(lon_d<-120.0_dp) &
559  .and.(lat_d>=-77.0_dp)) &
560  .or. &
561  ((lon_d>=-120.0_dp).and.(lon_d<-90.0_dp)) ) then
562  ! Amundsen Sea
563  n_sector(j,i) = 6_i1b
564 
565  toc = 0.42_dp
566  omega = 0.1023_dp
567  alpha = 0.26_dp
568  draft0 = 200.0_dp
569 
570 else if ( ((lon_d>=-90.0_dp).and.(lon_d<-66.0_dp) &
571  .and.(lat_d>=-74.0_dp)) ) then
572  ! Bellingshausen Sea
573  n_sector(j,i) = 7_i1b
574 
575  toc = 0.62_dp
576  omega = 0.0644_dp
577  alpha = 0.16_dp
578  draft0 = 200.0_dp
579 
580 else
581  ! Weddell Sea
582  n_sector(j,i) = 8_i1b
583 
584  toc = -1.55_dp
585  omega = 0.0197_dp
586  alpha = 0.26_dp
587  draft0 = 200.0_dp
588 
589 endif
590 
591 #else /* not Antarctic ice sheet */
592 
593 q_bm_floating = 0.0_dp ! dummy return value
594 
595 errormsg = ' >>> sub_ice_shelf_melting_param:' &
596  // end_of_line &
597  //' Case FLOATING_ICE_BASAL_MELTING==5' &
598  // end_of_line &
599  //' only defined for Antarctica!'
600 call error(errormsg)
601 
602 #endif
603 
604 #else
605 
606 !!! continue
607 
608 #endif
609 
610 !-------- Computation of the sub-ice-shelf melting rate --------
611 
612 if (maske(j,i)==2_i1b) then ! ocean
613  draft = 0.0_dp
614  h_w_now = max((z_sl-zl(j,i)), 0.0_dp)
615 else if (maske(j,i)==3_i1b) then ! floating ice
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)
618 else
619  errormsg = ' >>> sub_ice_shelf_melting_param:' &
620  // end_of_line &
621  //' Routine must not be called for maske(j,i) < 2!'
622  call error(errormsg)
623 end if
624 
625 tmb = -beta_sw*draft - delta_tm_sw ! temperature at the ice shelf base
626 
627 t_forcing = max((toc-tmb), 0.0_dp) ! thermal forcing
628 
629 #if (defined(H_W_0))
630 
631 if (h_w_0 > eps) then
632  q_bm_scaling_factor = tanh(h_w_now/h_w_0)
633 else
634  q_bm_scaling_factor = 1.0_dp
635 end if
636 
637 #else
638 q_bm_scaling_factor = 1.0_dp
639 #endif
640 
641 #if (FLOATING_ICE_BASAL_MELTING==4)
642 
643 q_bm_floating = q_bm_scaling_factor*omega*t_forcing**alpha
644 
645 #elif (FLOATING_ICE_BASAL_MELTING==5)
646 
647 phi_par = rho_sw*c_sw*g_t/(rho*l)
648 
649 #if !defined(ALLOW_OPENAD)
650 q_bm_floating = q_bm_scaling_factor &
651  *phi_par*omega*t_forcing*(draft/draft0)**alpha
652 #else
653 if (draft.eq.0) then
654 q_bm_floating = 0.0_dp
655 else
656 q_bm_floating = q_bm_scaling_factor &
657  *phi_par*omega*t_forcing*(draft/draft0)**alpha
658 end if
659 #endif
660 #else
661 
662 q_bm_floating = 0.0_dp ! dummy return value
663 
664 errormsg = ' >>> sub_ice_shelf_melting_param:' &
665  // end_of_line &
666  //' FLOATING_ICE_BASAL_MELTING must be 4 or 5!'
667 call error(errormsg)
668 
669 #endif
670 
671 ! ------ Correction for ISMIP InitMIP
672 
673 #if (defined(INITMIP_BMB_ANOM_FILE))
674 
675 if ((time_in_years > 0.0_dp).and.(time_in_years <= 40.0_dp)) then
676 
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)
680 #else /* OpenAD */
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 */
684 
685 else if (time_in_years > 40.0_dp) then
686 
687  q_bm_floating = q_bm_floating + ab_anom_initmip(j,i)
688 
689 end if
690 
691 #endif
692 
693 ! ------ Correction for ISMIP LARMIP
694 
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)
698 #endif
699 
700 end subroutine sub_ice_shelf_melting_param
701 
702 !-------------------------------------------------------------------------------
703 
704 end module calc_bas_melt_m
705 !
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)) ...