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