SICOPOLIS V3.0
 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-2013 Ralf Greve
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 
42 implicit none
43 
44 real(dp), intent(in) :: time, z_sl
45 real(dp), intent(in) :: dzeta_c, dzeta_r
46 
47 integer(i4b) :: i, j
48 integer(i4b) :: n_ocean, n_float
49 real(dp) :: kappa_val
50 real(dp) :: aqbm1, aqbm2, aqbm3, aqbm4
51 real(dp) :: year_sec_inv, rhow_rho_ratio
52 real(dp) :: q_bm_grounded, q_bm_marine, q_bm_floating
53 real(dp) :: qbm_min, qbm_max
54 
55 !-------- Term abbreviations --------
56 
57 year_sec_inv = 1.0_dp/year_sec
58 rhow_rho_ratio = rho_w/rho
59 
60 aqbm1 = (ea-1.0_dp)/(rho*l*deform*dzeta_c)
61 aqbm2 = kappa_r/(rho*l*h_r*dzeta_r)
62 aqbm3 = g/l
63 aqbm4 = beta/(rho*l)
64 
65 !-------- Computation of Q_bm --------
66 
67 q_bm = 0.0_dp ! initialisation
68 
69 do i=1, imax-1
70 do j=1, jmax-1
71 
72  if (maske(j,i)==0) then ! grounded ice
73 
74  if (n_cts(j,i)==-1) then
75  q_bm(j,i) = 0.0_dp
76  else if (n_cts(j,i)==0) then
77  q_bm(j,i) = aqbm1*(temp_c(1,j,i)-temp_c(0,j,i))/h_c(j,i) &
78  *kappa_val(temp_c(0,j,i)) &
79  -aqbm2*(temp_r(krmax,j,i)-temp_r(krmax-1,j,i)) &
80  -aqbm3*h_c(j,i)*0.5_dp*(vx_t(0,j,i)+vx_t(0,j,i-1)) &
81  *dzs_dxi_g(j,i) &
82  -aqbm3*h_c(j,i)*0.5_dp*(vy_t(0,j,i)+vy_t(0,j-1,i)) &
83  *dzs_deta_g(j,i)
84  else ! n_cts(j,i)==1
85  q_bm(j,i) = aqbm4 &
86  *kappa_val(temp_t_m(0,j,i)) &
87  -aqbm2*(temp_r(krmax,j,i)-temp_r(krmax-1,j,i)) &
88  -aqbm3*(h_c(j,i)+h_t(j,i)) &
89  *0.5_dp*(vx_t(0,j,i)+vx_t(0,j,i-1)) &
90  *dzs_dxi_g(j,i) &
91  -aqbm3*(h_c(j,i)+h_t(j,i)) &
92  *0.5_dp*(vy_t(0,j,i)+vy_t(0,j-1,i)) &
93  *dzs_deta_g(j,i)
94  end if
95 
96 #if ( MARGIN==1 || MARGIN==2 )
97 
98 #if ( !defined(MARINE_ICE_BASAL_MELTING) || MARINE_ICE_BASAL_MELTING==1 )
99 
100  continue ! do nothing
101 
102 #elif MARINE_ICE_BASAL_MELTING==2
103 
104  if ( (zb(j,i) < z_sl) & ! marine ice
105  .and. &
106  ( (maske(j,i+1)==2) & ! at least one
107  .or.(maske(j,i-1)==2) & ! nearest neighbour
108  .or.(maske(j+1,i)==2) & ! is
109  .or.(maske(j-1,i)==2) & ! ocean
110  ) &
111  ) then
112 
113  q_bm(j,i) = qbm_marine *year_sec_inv*rhow_rho_ratio
114  ! m/a water equiv. -> m/s ice equiv.
115 
116  end if
117 
118 #elif MARINE_ICE_BASAL_MELTING==3
119 
120  if ( (zb(j,i) < z_sl) & ! marine ice
121  .and. &
122  ( (maske(j,i+1)==2) & ! at least one
123  .or.(maske(j,i-1)==2) & ! nearest neighbour
124  .or.(maske(j+1,i)==2) & ! is
125  .or.(maske(j-1,i)==2) & ! ocean
126  ) &
127  ) then
128 
129  n_ocean = 0
130  if (maske(j,i+1)==2) n_ocean = n_ocean+1
131  if (maske(j,i-1)==2) n_ocean = n_ocean+1
132  if (maske(j+1,i)==2) n_ocean = n_ocean+1
133  if (maske(j-1,i)==2) n_ocean = n_ocean+1
134 
135  if ( n_ocean > 0 ) then
136 
137  q_bm_grounded = q_bm(j,i)
138  q_bm_marine = qbm_marine *year_sec_inv*rhow_rho_ratio
139  ! m/a water equiv. -> m/s ice equiv.
140 
141  q_bm(j,i) = (1.0_dp-0.25_dp*real(n_ocean,dp)) * q_bm_grounded &
142  +0.25_dp*real(n_ocean,dp) * q_bm_marine
143  ! weighed average of grounded ice melting (computed)
144  ! and marine ice melting (prescribed)
145  else
146  write(6,'(a)') ' calc_qbm: Marine ice margin point does not have'
147  write(6,'(a)') ' an ocean neighbour!'
148  stop
149  end if
150 
151  end if
152 
153 #else
154  stop ' calc_qbm: MARINE_ICE_BASAL_MELTING must be either 1, 2 or 3!'
155 #endif
156 
157 #elif MARGIN==3
158 
159  if ( flag_grounding_line(j,i) ) then ! grounding line
160 
161 #if FLOATING_ICE_BASAL_MELTING==1
162 
163  continue ! do nothing
164 
165 #elif FLOATING_ICE_BASAL_MELTING==2
166 
167  q_bm(j,i) = qbm_float_1 *year_sec_inv*rhow_rho_ratio
168  ! m/a water equiv. -> m/s ice equiv.
169 
170 #elif FLOATING_ICE_BASAL_MELTING==3
171 
172  n_float = 0
173  if (maske(j,i+1)>=2) n_float = n_float+1
174  if (maske(j,i-1)>=2) n_float = n_float+1
175  if (maske(j+1,i)>=2) n_float = n_float+1
176  if (maske(j-1,i)>=2) n_float = n_float+1
177 
178  if ( n_float > 0 ) then
179 
180  q_bm_grounded = q_bm(j,i)
181  q_bm_floating = qbm_float_1 *year_sec_inv*rhow_rho_ratio
182  ! m/a water equiv. -> m/s ice equiv.
183 
184  q_bm(j,i) = (1.0_dp-0.25_dp*real(n_float,dp)) * q_bm_grounded &
185  +0.25_dp*real(n_float,dp) * q_bm_floating
186  ! weighed average of melting for grounded and floating ice
187  else
188  write(6,'(a)') ' calc_qbm: Grounding line point does not have'
189  write(6,'(a)') ' a floating ice or ocean neighbour!'
190  stop
191  end if
192 
193 #else
194  stop ' calc_qbm: FLOATING_ICE_BASAL_MELTING must be either 1, 2 or 3!'
195 #endif
196 
197  else if ( (zb(j,i) < z_sl) & ! marine ice margin
198  .and. &
199  ( (maske(j,i+1)>=2) & ! (at least one
200  .or.(maske(j,i-1)>=2) & ! nearest neighbour
201  .or.(maske(j+1,i)>=2) & ! is
202  .or.(maske(j-1,i)>=2) & ! ocean)
203  ) &
204  ) then
205 
206 #if MARINE_ICE_BASAL_MELTING==1
207 
208  continue ! do nothing
209 
210 #elif MARINE_ICE_BASAL_MELTING==2
211 
212  q_bm(j,i) = qbm_marine *year_sec_inv*rhow_rho_ratio
213  ! m/a water equiv. -> m/s ice equiv.
214 
215 #elif MARINE_ICE_BASAL_MELTING==3
216 
217  n_ocean = 0
218  if (maske(j,i+1)>=2) n_ocean = n_ocean+1
219  if (maske(j,i-1)>=2) n_ocean = n_ocean+1
220  if (maske(j+1,i)>=2) n_ocean = n_ocean+1
221  if (maske(j-1,i)>=2) n_ocean = n_ocean+1
222 
223  if ( n_ocean > 0 ) then
224 
225  q_bm_grounded = q_bm(j,i)
226  q_bm_marine = qbm_marine *year_sec_inv*rhow_rho_ratio
227  ! m/a water equiv. -> m/s ice equiv.
228 
229  q_bm(j,i) = (1.0_dp-0.25_dp*real(n_ocean,dp)) * q_bm_grounded &
230  +0.25_dp*real(n_ocean,dp) * q_bm_marine
231  ! weighed average of grounded ice melting (computed)
232  ! and marine ice melting (prescribed)
233  else
234  write(6,'(a)') ' calc_qbm: Marine ice margin point does not have'
235  write(6,'(a)') ' a floating ice or ocean neighbour!'
236  stop
237  end if
238 
239 #else
240  stop ' calc_qbm: MARINE_ICE_BASAL_MELTING must be either 1, 2 or 3!'
241 #endif
242 
243  end if
244 
245  else if ( (maske(j,i)==2).or.(maske(j,i)==3) ) then ! floating ice or ocean
246 
247  if ( flag_grounding_line(j,i+1) & ! at least
248  .or.flag_grounding_line(j,i-1) & ! one neighboring point
249  .or.flag_grounding_line(j+1,i) & ! is grounding line
250  .or.flag_grounding_line(j-1,i) ) then
251 
252  q_bm(j,i) = qbm_float_1 *year_sec_inv*rhow_rho_ratio
253  ! m/a water equiv. -> m/s ice equiv.
254 
255  else if ( zl(j,i) > z_abyss ) then ! floating ice over continental shelf
256 
257  q_bm(j,i) = qbm_float_2 *year_sec_inv*rhow_rho_ratio
258  ! m/a water equiv. -> m/s ice equiv.
259 
260  else ! ( zl(j,i) <= Z_ABYSS ), floating ice over abyssal ocean
261 
262  q_bm(j,i) = qbm_float_3 *year_sec_inv*rhow_rho_ratio
263  ! m/a water equiv. -> m/s ice equiv.
264 
265  end if
266 
267 #else
268  stop ' calc_qbm: MARGIN must be either 1, 2 or 3!'
269 #endif
270 
271  end if
272 
273 end do
274 end do
275 
276 !-------- Sum of Q_bm and Q_tld --------
277 
278 q_b_tot = q_bm + q_tld
279 
280 !-------- Limitation of Q_bm, Q_tld and Q_b_tot
281 ! (only for grounded ice) --------
282 
283 #if defined(QBM_MIN)
284  qbm_min = qbm_min *year_sec_inv*rhow_rho_ratio
285  ! m/a water equiv. -> m/s ice equiv.
286 #elif defined(QB_MIN) /* obsolete */
287  qbm_min = qb_min ! m/s ice equiv.
288 #else
289  stop ' calc_qbm: Limiter QBM_MIN is undefined!'
290 #endif
291 
292 #if defined(QBM_MAX)
293  qbm_max = qbm_max *year_sec_inv*rhow_rho_ratio
294  ! m/a water equiv. -> m/s ice equiv.
295 #elif defined(QB_MAX) /* obsolete */
296  qbm_max = qb_max ! m/s ice equiv.
297 #else
298  stop ' calc_qbm: Limiter QBM_MAX is undefined!'
299 #endif
300 
301 #if ( defined(MARINE_ICE_BASAL_MELTING) \
302  && ( marine_ice_basal_melting==2 || marine_ice_basal_melting==3 ) )
303 
304 if (qbm_marine*year_sec_inv*rhow_rho_ratio > qbm_max) then
305  write (6,'(a)') ' calc_qbm: QBM_MARINE (basal melting rate at the marine'
306  write (6,'(a)') ' ice front) is larger than the limiter qbm_max!'
307  stop
308 end if
309 
310 #endif
311 
312 #if ( MARGIN==3 \
313  && defined(floating_ice_basal_melting) \
314  && floating_ice_basal_melting<=3 )
315 
316 if (qbm_float_1*year_sec_inv*rhow_rho_ratio > qbm_max) then
317  write (6,'(a)') ' calc_qbm: QBM_FLOAT_1 (basal melting rate in the grounding'
318  write (6,'(a)') ' line zone) is larger than the limiter qbm_max!'
319  stop
320 end if
321 
322 #endif
323 
324 do i=0, imax
325 do j=0, jmax
326  if (maske(j,i)==0) then ! grounded ice
327  if (q_bm(j,i) < qbm_min) q_bm(j,i) = 0.0_dp
328  if (q_bm(j,i) > qbm_max) q_bm(j,i) = qbm_max
329  if (q_tld(j,i) < qbm_min) q_tld(j,i) = 0.0_dp
330  if (q_tld(j,i) > qbm_max) q_tld(j,i) = qbm_max
331  if (q_b_tot(j,i) < qbm_min) q_b_tot(j,i) = 0.0_dp
332  if (q_b_tot(j,i) > qbm_max) q_b_tot(j,i) = qbm_max
333  end if
334 end do
335 end do
336 
337 end subroutine calc_qbm
338 !
Generated on Mon May 20 2013 10:18:40 for SICOPOLIS V3.0 by   doxygen 1.8.1