SICOPOLIS V3.2
 All Classes Files Functions Variables Macros
calc_qbm.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Subroutine : c a l c _ q b m
4 !
5 !> @file
6 !!
7 !! Computation of the basal melting rate Q_bm.
8 !! Summation of Q_bm and Q_tld.
9 !!
10 !! @section Copyright
11 !!
12 !! Copyright 2009-2016 Ralf Greve, Tatsuru Sato
13 !!
14 !! @section License
15 !!
16 !! This file is part of SICOPOLIS.
17 !!
18 !! SICOPOLIS is free software: you can redistribute it and/or modify
19 !! it under the terms of the GNU General Public License as published by
20 !! the Free Software Foundation, either version 3 of the License, or
21 !! (at your option) any later version.
22 !!
23 !! SICOPOLIS is distributed in the hope that it will be useful,
24 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
25 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
26 !! GNU General Public License for more details.
27 !!
28 !! You should have received a copy of the GNU General Public License
29 !! along with SICOPOLIS. If not, see <http://www.gnu.org/licenses/>.
30 !<
31 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
32 
33 !-------------------------------------------------------------------------------
34 !> Computation of the basal melting rate Q_bm.
35 !! Summation of Q_bm and Q_tld.
36 !<------------------------------------------------------------------------------
37 subroutine calc_qbm(time, z_sl, dzeta_c, dzeta_r)
38 
39 use sico_types
40 use sico_variables
41 use sico_vars
42 use ice_material_quantities, only : kappa_val
43 
44 implicit none
45 
46 real(dp), intent(in) :: time, z_sl
47 real(dp), intent(in) :: dzeta_c, dzeta_r
48 
49 integer(i4b) :: i, j
50 integer(i4b) :: n_ocean, n_float
51 real(dp) :: aqbm1, aqbm2, aqbm3a, aqbm3b, aqbm4
52 real(dp) :: frictional_heating
53 real(dp) :: year_sec_inv, rhow_rho_ratio
54 real(dp) :: q_bm_grounded, q_bm_marine, q_bm_floating
55 real(dp) :: qbm_min, qbm_max
56 
57 !-------- Term abbreviations --------
58 
59 year_sec_inv = 1.0_dp/year_sec
60 rhow_rho_ratio = rho_w/rho
61 
62 if (flag_aa_nonzero) then
63  aqbm1 = (ea-1.0_dp)/(rho*l*aa*dzeta_c)
64 else
65  aqbm1 = 1.0_dp/(rho*l*dzeta_c)
66 end if
67 
68 aqbm2 = kappa_r/(rho*l*h_r*dzeta_r)
69 aqbm3a = g/l
70 aqbm3b = 1.0_dp/(rho*l)
71 aqbm4 = beta/(rho*l)
72 
73 !-------- Computation of Q_bm --------
74 
75 q_bm = 0.0_dp ! initialisation
76 
77 do i=1, imax-1
78 do j=1, jmax-1
79 
80  if (maske(j,i)==0) then ! grounded ice
81 
82  if (n_cts(j,i)==-1) then
83 
84  frictional_heating = 0.0_dp
85  q_bm(j,i) = 0.0_dp
86 
87  else if (n_cts(j,i)==0) then
88 
89 #if (DYNAMICS==2)
90  if (.not.flag_shelfy_stream(j,i)) then
91 #endif
92  frictional_heating &
93  = -aqbm3a*h_c(j,i) &
94  *0.5_dp*(vx_t(0,j,i)+vx_t(0,j,i-1)) &
95  *dzs_dxi_g(j,i) &
96  -aqbm3a*h_c(j,i) &
97  *0.5_dp*(vy_t(0,j,i)+vy_t(0,j-1,i)) &
98  *dzs_deta_g(j,i)
99 #if (DYNAMICS==2)
100  else ! flag_shelfy_stream(j,i) == .true.
101 
102  frictional_heating &
103  = aqbm3b &
104  * c_drag(j,i) &
105  * sqrt(vx_b_g(j,i)**2 &
106  +vy_b_g(j,i)**2) &
107  **(1.0_dp+p_weert_inv(j,i))
108  end if
109 #endif
110  q_bm(j,i) = aqbm1*(temp_c(1,j,i)-temp_c(0,j,i))/h_c(j,i) &
111  *kappa_val(temp_c(0,j,i)) &
112  - aqbm2*(temp_r(krmax,j,i)-temp_r(krmax-1,j,i)) &
113  + frictional_heating
114 
115  else ! n_cts(j,i)==1
116 
117 #if (DYNAMICS==2)
118  if (.not.flag_shelfy_stream(j,i)) then
119 #endif
120  frictional_heating &
121  = -aqbm3a*(h_c(j,i)+h_t(j,i)) &
122  *0.5_dp*(vx_t(0,j,i)+vx_t(0,j,i-1)) &
123  *dzs_dxi_g(j,i) &
124  -aqbm3a*(h_c(j,i)+h_t(j,i)) &
125  *0.5_dp*(vy_t(0,j,i)+vy_t(0,j-1,i)) &
126  *dzs_deta_g(j,i)
127 #if (DYNAMICS==2)
128  else ! flag_shelfy_stream(j,i) == .true.
129 
130  frictional_heating &
131  = aqbm3b &
132  * c_drag(j,i) &
133  * sqrt(vx_b_g(j,i)**2 &
134  +vy_b_g(j,i)**2) &
135  **(1.0_dp+p_weert_inv(j,i))
136  end if
137 #endif
138  q_bm(j,i) = aqbm4*kappa_val(temp_t_m(0,j,i)) &
139  - aqbm2*(temp_r(krmax,j,i)-temp_r(krmax-1,j,i)) &
140  + frictional_heating
141 
142  end if
143 
144 #if (MARGIN==1 || MARGIN==2)
145 
146 #if (!defined(MARINE_ICE_BASAL_MELTING) || MARINE_ICE_BASAL_MELTING==1)
147 
148  continue ! do nothing
149 
150 #elif (MARINE_ICE_BASAL_MELTING==2)
151 
152  if ( (zb(j,i) < z_sl) & ! marine ice
153  .and. &
154  ( (maske(j,i+1)==2) & ! at least one
155  .or.(maske(j,i-1)==2) & ! nearest neighbour
156  .or.(maske(j+1,i)==2) & ! is
157  .or.(maske(j-1,i)==2) & ! ocean
158  ) &
159  ) then
160 
161  q_bm(j,i) = qbm_marine *year_sec_inv*rhow_rho_ratio
162  ! m/a water equiv. -> m/s ice equiv.
163 
164  end if
165 
166 #elif (MARINE_ICE_BASAL_MELTING==3)
167 
168  if ( (zb(j,i) < z_sl) & ! marine ice
169  .and. &
170  ( (maske(j,i+1)==2) & ! at least one
171  .or.(maske(j,i-1)==2) & ! nearest neighbour
172  .or.(maske(j+1,i)==2) & ! is
173  .or.(maske(j-1,i)==2) & ! ocean
174  ) &
175  ) then
176 
177  n_ocean = 0
178  if (maske(j,i+1)==2) n_ocean = n_ocean+1
179  if (maske(j,i-1)==2) n_ocean = n_ocean+1
180  if (maske(j+1,i)==2) n_ocean = n_ocean+1
181  if (maske(j-1,i)==2) n_ocean = n_ocean+1
182 
183  if ( n_ocean > 0 ) then
184 
185  q_bm_grounded = q_bm(j,i)
186  q_bm_marine = qbm_marine *year_sec_inv*rhow_rho_ratio
187  ! m/a water equiv. -> m/s ice equiv.
188 
189  q_bm(j,i) = (1.0_dp-0.25_dp*real(n_ocean,dp)) * q_bm_grounded &
190  +0.25_dp*real(n_ocean,dp) * q_bm_marine
191  ! weighed average of grounded ice melting (computed)
192  ! and marine ice melting (prescribed)
193  else
194  write(6,'(a)') ' calc_qbm: Marine ice margin point does not have'
195  write(6,'(a)') ' an ocean neighbour!'
196  stop
197  end if
198 
199  end if
200 
201 #else
202  stop ' calc_qbm: MARINE_ICE_BASAL_MELTING must be either 1, 2 or 3!'
203 #endif
204 
205 #elif (MARGIN==3)
206 
207  if (flag_grounding_line_1(j,i)) then
208  ! grounding line (grounded-ice side)
209 
210 #if (FLOATING_ICE_BASAL_MELTING==1 || FLOATING_ICE_BASAL_MELTING==4)
211 
212  continue ! do nothing
213 
214 #elif (FLOATING_ICE_BASAL_MELTING==2)
215 
216  q_bm(j,i) = qbm_float_1 *year_sec_inv*rhow_rho_ratio
217  ! m/a water equiv. -> m/s ice equiv.
218 
219 #elif (FLOATING_ICE_BASAL_MELTING==3)
220 
221  n_float = 0
222  if (maske(j,i+1)>=2) n_float = n_float+1
223  if (maske(j,i-1)>=2) n_float = n_float+1
224  if (maske(j+1,i)>=2) n_float = n_float+1
225  if (maske(j-1,i)>=2) n_float = n_float+1
226 
227  if ( n_float > 0 ) then
228 
229  q_bm_grounded = q_bm(j,i)
230  q_bm_floating = qbm_float_1 *year_sec_inv*rhow_rho_ratio
231  ! m/a water equiv. -> m/s ice equiv.
232 
233  q_bm(j,i) = (1.0_dp-0.25_dp*real(n_float,dp)) * q_bm_grounded &
234  +0.25_dp*real(n_float,dp) * q_bm_floating
235  ! weighed average of melting for grounded and floating ice
236  else
237  write(6,'(a)') ' calc_qbm: Grounding line point does not have'
238  write(6,'(a)') ' a floating ice or ocean neighbour!'
239  stop
240  end if
241 
242 #else
243  write(6, fmt='(a)') ' calc_qbm: FLOATING_ICE_BASAL_MELTING'
244  write(6, fmt='(a)') ' must be either 1, 2, 3 or 4!'
245  stop
246 #endif
247 
248  else if ( (zb(j,i) < z_sl) & ! marine ice margin
249  .and. &
250  ( (maske(j,i+1)>=2) & ! (at least one
251  .or.(maske(j,i-1)>=2) & ! nearest neighbour
252  .or.(maske(j+1,i)>=2) & ! is
253  .or.(maske(j-1,i)>=2) & ! ocean)
254  ) &
255  ) then
256 
257 #if (MARINE_ICE_BASAL_MELTING==1)
258 
259  continue ! do nothing
260 
261 #elif (MARINE_ICE_BASAL_MELTING==2)
262 
263  q_bm(j,i) = qbm_marine *year_sec_inv*rhow_rho_ratio
264  ! m/a water equiv. -> m/s ice equiv.
265 
266 #elif (MARINE_ICE_BASAL_MELTING==3)
267 
268  n_ocean = 0
269  if (maske(j,i+1)>=2) n_ocean = n_ocean+1
270  if (maske(j,i-1)>=2) n_ocean = n_ocean+1
271  if (maske(j+1,i)>=2) n_ocean = n_ocean+1
272  if (maske(j-1,i)>=2) n_ocean = n_ocean+1
273 
274  if ( n_ocean > 0 ) then
275 
276  q_bm_grounded = q_bm(j,i)
277  q_bm_marine = qbm_marine *year_sec_inv*rhow_rho_ratio
278  ! m/a water equiv. -> m/s ice equiv.
279 
280  q_bm(j,i) = (1.0_dp-0.25_dp*real(n_ocean,dp)) * q_bm_grounded &
281  +0.25_dp*real(n_ocean,dp) * q_bm_marine
282  ! weighed average of grounded ice melting (computed)
283  ! and marine ice melting (prescribed)
284  else
285  write(6,'(a)') ' calc_qbm: Marine ice margin point does not have'
286  write(6,'(a)') ' a floating ice or ocean neighbour!'
287  stop
288  end if
289 
290 #else
291  stop ' calc_qbm: MARINE_ICE_BASAL_MELTING must be either 1, 2 or 3!'
292 #endif
293 
294  end if
295 
296  else if ( (maske(j,i)==2).or.(maske(j,i)==3) ) then ! floating ice or ocean
297 
298 #if (FLOATING_ICE_BASAL_MELTING<=3)
299 
300  if (flag_grounding_line_2(j,i)) then
301  ! grounding line (floating-ice side)
302 
303  q_bm(j,i) = qbm_float_1 *year_sec_inv*rhow_rho_ratio
304  ! m/a water equiv. -> m/s ice equiv.
305 
306  else if ( zl(j,i) > z_abyss ) then ! floating ice over continental shelf
307 
308  q_bm(j,i) = qbm_float_2 *year_sec_inv*rhow_rho_ratio
309  ! m/a water equiv. -> m/s ice equiv.
310 
311  else ! ( zl(j,i) <= Z_ABYSS ), floating ice over abyssal ocean
312 
313  q_bm(j,i) = qbm_float_3 *year_sec_inv*rhow_rho_ratio
314  ! m/a water equiv. -> m/s ice equiv.
315 
316  end if
317 
318 #elif (FLOATING_ICE_BASAL_MELTING==4)
319 
320  if ( zl(j,i) > z_abyss ) then ! floating ice over continental shelf
321  call sub_ice_shelf_melting_param(z_sl, i, j, q_bm_floating)
322  q_bm(j,i) = q_bm_floating
323  else ! ( zl(j,i) <= Z_ABYSS ), floating ice over abyssal ocean
324  q_bm(j,i) = qbm_float_3 *year_sec_inv*rhow_rho_ratio
325  ! m/a water equiv. -> m/s ice equiv.
326  end if
327 
328 #else
329  write(6, fmt='(a)') ' calc_qbm: FLOATING_ICE_BASAL_MELTING'
330  write(6, fmt='(a)') ' must be either 1, 2, 3 or 4!'
331  stop
332 #endif
333 
334 #else
335  stop ' calc_qbm: MARGIN must be either 1, 2 or 3!'
336 #endif
337 
338  end if
339 
340 end do
341 end do
342 
343 !-------- Sum of Q_bm and Q_tld --------
344 
345 q_b_tot = q_bm + q_tld
346 
347 !-------- Limitation of Q_bm, Q_tld and Q_b_tot
348 ! (only for grounded ice) --------
349 
350 #if (defined(QBM_MIN))
351  qbm_min = qbm_min *year_sec_inv*rhow_rho_ratio
352  ! m/a water equiv. -> m/s ice equiv.
353 #elif (defined(QB_MIN)) /* obsolete */
354  qbm_min = qb_min ! m/s ice equiv.
355 #else
356  stop ' calc_qbm: Limiter QBM_MIN is undefined!'
357 #endif
358 
359 #if (defined(QBM_MAX))
360  qbm_max = qbm_max *year_sec_inv*rhow_rho_ratio
361  ! m/a water equiv. -> m/s ice equiv.
362 #elif (defined(QB_MAX)) /* obsolete */
363  qbm_max = qb_max ! m/s ice equiv.
364 #else
365  stop ' calc_qbm: Limiter QBM_MAX is undefined!'
366 #endif
367 
368 #if ( defined(MARINE_ICE_BASAL_MELTING) \
369  && ( marine_ice_basal_melting==2 || marine_ice_basal_melting==3 ) )
370 
371 if (qbm_marine*year_sec_inv*rhow_rho_ratio > qbm_max) then
372  write (6,'(a)') ' calc_qbm: QBM_MARINE (basal melting rate at the marine'
373  write (6,'(a)') ' ice front) is larger than the limiter qbm_max!'
374  stop
375 end if
376 
377 #endif
378 
379 #if ( MARGIN==3 \
380  && defined(floating_ice_basal_melting) \
381  && floating_ice_basal_melting<=3 )
382 
383 if (qbm_float_1*year_sec_inv*rhow_rho_ratio > qbm_max) then
384  write (6,'(a)') ' calc_qbm: QBM_FLOAT_1 (basal melting rate in the grounding'
385  write (6,'(a)') ' line zone) is larger than the limiter qbm_max!'
386  stop
387 end if
388 
389 #endif
390 
391 do i=0, imax
392 do j=0, jmax
393  if (maske(j,i)==0) then ! grounded ice
394  if (q_bm(j,i) < qbm_min) q_bm(j,i) = 0.0_dp
395  if (q_bm(j,i) > qbm_max) q_bm(j,i) = qbm_max
396  if (q_tld(j,i) < qbm_min) q_tld(j,i) = 0.0_dp
397  if (q_tld(j,i) > qbm_max) q_tld(j,i) = qbm_max
398  if (q_b_tot(j,i) < qbm_min) q_b_tot(j,i) = 0.0_dp
399  if (q_b_tot(j,i) > qbm_max) q_b_tot(j,i) = qbm_max
400  end if
401 end do
402 end do
403 
404 end subroutine calc_qbm
405 !
sico_types
Declarations of kind types for SICOPOLIS.
Definition: sico_types.F90:35
sico_vars
Declarations of global variables for SICOPOLIS (for the ANT domain).
Definition: sico_vars.F90:35
calc_qbm
subroutine calc_qbm(time, z_sl, dzeta_c, dzeta_r)
Computation of the basal melting rate Q_bm. Summation of Q_bm and Q_tld.
Definition: calc_qbm.F90:37
sub_ice_shelf_melting_param
subroutine sub_ice_shelf_melting_param(z_sl, i, j, Q_bm_floating)
Sub-ice-shelf melting parameterisation for Antarctica.
Definition: sub_ice_shelf_melting_param.F90:35
ice_material_quantities::kappa_val
real(dp) function, public kappa_val(temp_val)
Heat conductivity of ice: Linear interpolation of tabulated values in KAPPA(.).
Definition: ice_material_quantities.F90:217
ice_material_quantities
Material quantities of ice: Rate factor, heat conductivity, specific heat (heat capacity).
Definition: ice_material_quantities.F90:37
sico_variables
Declarations of global variables for SICOPOLIS.
Definition: sico_variables.F90:35
Generated on Wed Mar 2 2016 11:26:54 for SICOPOLIS V3.2 by   doxygen 1.8.5