SICOPOLIS V5-dev  Revision 1368
boundary_m.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Module : b o u n d a r y _ m
4 !
5 !> @file
6 !!
7 !! Computation of the surface temperature (must be less than 0 deg C!)
8 !! and of the accumulation-ablation function.
9 !!
10 !! @section Copyright
11 !!
12 !! Copyright 2009-2018 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 surface temperature (must be less than 0 deg C!)
35 !! and of the accumulation-ablation function.
36 !<------------------------------------------------------------------------------
37 module boundary_m
38 
39  use sico_types_m
41  use sico_vars_m
42  use error_m
43 
44  implicit none
45 
46 #if (INITMIP_CONST_SMB==1) /* SMB held constant at its initial distribution */
47  real(dp), dimension(0:JMAX,0:IMAX), save :: accum_save, evap_save, runoff_save
48 #endif
49 
50  public
51 
52 contains
53 
54 !-------------------------------------------------------------------------------
55 !> Main routine of boundary_m:
56 !! Computation of the surface temperature (must be less than 0 deg C!)
57 !! and of the accumulation-ablation function.
58 !<------------------------------------------------------------------------------
59 subroutine boundary(time, dtime, dxi, deta, &
60  delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
61 
62 #if (NETCDF > 1)
63  use netcdf
64  use nc_check_m
65 #endif
66 
67 #if ((MARGIN==2) \
68  && (marine_ice_formation==2) \
69  && (marine_ice_calving==9))
71 #endif
72 
74  use pdd_m
75 
76 #if (defined(ALLOW_GRDCHK) || defined(ALLOW_OPENAD)) /* OpenAD */
77  use ctrl_m, only: myfloor, myceiling
78  use sico_maths_m
79 #endif /* OpenAD */
80 
81 implicit none
82 
83 real(dp), intent(in) :: time, dtime, dxi, deta
84 
85 real(dp), intent(out) :: delta_ts, glac_index, dzsl_dtau, z_mar
86 real(dp), intent(inout) :: z_sl
87 
88 ! Further return variables
89 ! (defined as global variables in module sico_variables_m):
90 !
91 ! accum(j,i), evap(j,i), runoff(j,i), as_perp(j,i),
92 ! calv_grounded(j,i), temp_s(j,i)
93 
94 integer(i4b) :: i, j, n
95 integer(i4b) :: i_gr, i_kl
96 integer(i4b) :: nrec_temp_precip
97 integer(i4b), save :: nrec_temp_precip_save = -1
98 real(dp) :: z_sl_old
99 real(dp) :: z_sl_min, t1, t2, t3, t4, t5, t6
100 real(dp) :: year_sec_inv, time_in_years
101 real(dp) :: time_gr, time_kl
102 real(dp) :: z_sle_present, z_sle_help
103 real(dp), dimension(0:JMAX,0:IMAX,0:12) :: precip
104 real(dp), dimension(0:JMAX,0:IMAX,12) :: temp_mm
105 real(dp), dimension(0:JMAX,0:IMAX) :: temp_ma, temp_ampl
106 real(dp), dimension(0:JMAX,0:IMAX) :: temp_ma_anom, temp_mj_anom, precip_ma_anom
107 real(dp), dimension(0:IMAX,0:JMAX), save :: temp_ma_anom_tra, temp_mj_anom_tra, &
108  precip_ma_anom_tra
109  ! These arrays are transposed with respect
110  ! to the usual SICOPOLIS convention;
111  ! they have i as first and j as second index
112 real(dp), dimension(12) :: temp_mm_help
113 real(dp) :: temp_jja_help
114 real(dp) :: theta_ma, theta_ma_offset, c_ma, gamma_ma, &
115  theta_ma_1, c_ma_1, gamma_ma_1, &
116  theta_ma_2, c_ma_2, gamma_ma_2, &
117  theta_ma_3, c_ma_3, gamma_ma_3, &
118  zs_sep_1, zs_sep_2, &
119  theta_mj, theta_mj_offset, c_mj, gamma_mj
120 real(dp) :: sine_factor
121 real(dp) :: gamma_p, zs_thresh, &
122  alpha_p, beta_p, temp_0, alpha_t, beta_t, &
123  temp_inv, temp_inv_present, &
124  temp_rain, temp_snow, &
125  inv_delta_temp_rain_snow, coeff(0:5), inv_sqrt2_s_stat, &
126  precip_fact, frac_solid
127 real(dp) :: s_stat, beta1, beta2, Pmax, mu, lambda_lti, temp_lti
128 logical, dimension(0:JMAX,0:IMAX) :: check_point
129 logical, save :: firstcall = .true.
130 
131 #if (NETCDF > 1)
132 integer(i4b) :: ncv
133 ! ncv: Variable ID
134 integer(i4b) :: nc3cor(3)
135 ! nc3cor(3): Corner of a 3-d array
136 integer(i4b) :: nc3cnt(3)
137 ! nc3cnt(3): Count of a 3-d array
138 #endif
139 
140 real(dp), parameter :: &
141  inv_twelve = 1.0_dp/12.0_dp, one_third = 1.0_dp/3.0_dp
142 
143 character(len=64), parameter :: thisroutine = 'boundary'
144 
145 #if (defined(ALLOW_GRDCHK) || defined(ALLOW_OPENAD)) /* OpenAD */
146 integer(i4b) :: i_time_in_years
147 real(dp) :: temp_val
148 #endif /* OpenAD */
149 
150 year_sec_inv = 1.0_dp/year_sec
151 time_in_years = time*year_sec_inv
152 
153 !-------- Initialization of intent(out) variables --------
154 
155 z_sl_old = z_sl
156 
157 delta_ts = 0.0_dp
158 glac_index = 0.0_dp
159 z_sl = 0.0_dp
160 dzsl_dtau = 0.0_dp
161 z_mar = 0.0_dp
162 
163 !-------- Surface-temperature deviation from present values --------
164 
165 #if (TSURFACE==1)
166 delta_ts = delta_ts0
167 ! ! Steady state with prescribed constant
168 ! ! air-temperature deviation
169 #elif (TSURFACE==3)
170 delta_ts = sine_amplit &
171  *cos(2.0_dp*pi*time/(sine_period*year_sec)) &
172  -sine_amplit
173 ! ! Sinusoidal air-temperature forcing
174 #elif (TSURFACE==4)
175 
176 ! ------ delta_ts from ice-core record
177 
178 if (time/year_sec < real(grip_time_min,dp)) then
179  delta_ts = griptemp(0)
180 else if (time/year_sec < real(grip_time_max,dp)) then
181 
182 #if (!defined(ALLOW_GRDCHK) && !defined(ALLOW_OPENAD)) /* Normal */
183  i_kl = floor(((time/year_sec) &
184  -real(grip_time_min,dp))/real(grip_time_stp,dp))
185  i_kl = max(i_kl, 0)
186 
187  i_gr = ceiling(((time/year_sec) &
188  -real(grip_time_min,dp))/real(grip_time_stp,dp))
189  i_gr = min(i_gr, ndata_grip)
190 #else /* OpenAD */
191  call myfloor((((time/year_sec) &
192  -real(grip_time_min,dp))/real(grip_time_stp,dp)), i_kl)
193  i_kl = max(i_kl, 0)
194 
195  call myceiling((((time/year_sec) &
196  -real(grip_time_min,dp))/real(grip_time_stp,dp)), i_gr)
197  i_gr = min(i_gr, ndata_grip)
198 #endif /* Normal vs. OpenAD */
199 
200  if (i_kl == i_gr) then
201 
202  delta_ts = griptemp(i_kl)
203 
204  else
205 
206  time_kl = (grip_time_min + i_kl*grip_time_stp) *year_sec
207  time_gr = (grip_time_min + i_gr*grip_time_stp) *year_sec
208 
209  delta_ts = griptemp(i_kl) &
210  +(griptemp(i_gr)-griptemp(i_kl)) &
211  *(time-time_kl)/(time_gr-time_kl)
212  ! linear interpolation of the ice-core data
213 
214  end if
215 
216 else
217  delta_ts = griptemp(ndata_grip)
218 end if
219 
220 delta_ts = delta_ts * grip_temp_fact
221 ! ! modification by constant factor
222 
223 !-------- Glacial index --------
224 
225 #elif (TSURFACE==5)
226 
227 if (time/year_sec < real(gi_time_min,dp)) then
228  glac_index = glacial_index(0)
229 else if (time/year_sec < real(gi_time_max,dp)) then
230 
231 #if (!defined(ALLOW_GRDCHK) && !defined(ALLOW_OPENAD)) /* Normal */
232  i_kl = floor(((time/year_sec) &
233  -real(gi_time_min,dp))/real(gi_time_stp,dp))
234  i_kl = max(i_kl, 0)
235 
236  i_gr = ceiling(((time/year_sec) &
237  -real(gi_time_min,dp))/real(gi_time_stp,dp))
238  i_gr = min(i_gr, ndata_gi)
239 #else /* OpenAD */
240  call myfloor((((time/year_sec) &
241  -real(gi_time_min,dp))/real(gi_time_stp,dp)),i_kl)
242  i_kl = max(i_kl, 0)
243 
244  call myceiling((((time/year_sec) &
245  -real(gi_time_min,dp))/real(gi_time_stp,dp)),i_gr)
246  i_gr = min(i_gr, ndata_gi)
247 #endif /* Normal vs. OpenAD */
248 
249  if (i_kl == i_gr) then
250 
251  glac_index = glacial_index(i_kl)
252 
253  else
254 
255  time_kl = (gi_time_min + i_kl*gi_time_stp) *year_sec
256  time_gr = (gi_time_min + i_gr*gi_time_stp) *year_sec
257 
258  glac_index = glacial_index(i_kl) &
259  +(glacial_index(i_gr)-glacial_index(i_kl)) &
260  *(time-time_kl)/(time_gr-time_kl)
261  ! linear interpolation of the glacial-index data
262 
263  end if
264 
265 else
266  glac_index = glacial_index(ndata_gi)
267 end if
268 
269 !-------- Reading of surface temperature and precipitation
270 ! directly from data --------
271 
272 #elif (TSURFACE==6 && ACCSURFACE==6)
273 
274 if (time/year_sec <= temp_precip_time_min) then
275  nrec_temp_precip = 0
276 else if (time/year_sec < temp_precip_time_max) then
277  nrec_temp_precip = nint( ((time/year_sec)-temp_precip_time_min) &
279 else
280  nrec_temp_precip = ndata_temp_precip
281 end if
282 
283 if ( nrec_temp_precip < 0 ) then
284  errormsg = ' >>> boundary: nrec_temp_precip < 0, not allowed!'
285  call error(errormsg)
286 else if ( nrec_temp_precip > ndata_temp_precip ) then
287  errormsg = ' >>> boundary: ' &
288  //'nrec_temp_precip > ndata_temp_precip, not allowed!'
289  call error(errormsg)
290 end if
291 
292 if ( nrec_temp_precip /= nrec_temp_precip_save ) then
293 
294  nc3cor(1) = 1
295  nc3cor(2) = 1
296  nc3cor(3) = nrec_temp_precip + 1
297  nc3cnt(1) = imax + 1
298  nc3cnt(2) = jmax + 1
299  nc3cnt(3) = 1
300 
301  call check( nf90_inq_varid(ncid_temp_precip, 'annualtemp_anom', ncv), &
302  thisroutine )
303  call check( nf90_get_var(ncid_temp_precip, ncv, temp_ma_anom_tra, &
304  start=nc3cor, count=nc3cnt), thisroutine )
305 
306  call check( nf90_inq_varid(ncid_temp_precip, 'januarytemp_anom', ncv), &
307  thisroutine )
308  call check( nf90_get_var(ncid_temp_precip, ncv, temp_mj_anom_tra, &
309  start=nc3cor, count=nc3cnt), thisroutine )
310 
311  call check( nf90_inq_varid(ncid_temp_precip, 'precipitation_anom', ncv), &
312  thisroutine )
313  call check( nf90_get_var(ncid_temp_precip, ncv, precip_ma_anom_tra, &
314  start=nc3cor, count=nc3cnt), thisroutine )
315 
316 end if
317 
318 temp_ma_anom = transpose(temp_ma_anom_tra)
319 temp_mj_anom = transpose(temp_mj_anom_tra)
320 precip_ma_anom = transpose(precip_ma_anom_tra) *(1.0e-03_dp/year_sec)*(rho_w/rho)
321  ! mm/a water equiv. -> m/s ice equiv.
322 
323 nrec_temp_precip_save = nrec_temp_precip
324 
325 #endif
326 
327 !-------- Sea level --------
328 
329 #if (SEA_LEVEL==1)
330 ! ------ constant sea level
331 z_sl = z_sl0
332 
333 #elif (SEA_LEVEL==2)
334 ! ------ saw-tooth-shaped palaeoclimatic sea-level forcing
335 
336 z_sl_min = -130.0_dp
337 
338 t1 = -250000.0_dp *year_sec
339 t2 = -140000.0_dp *year_sec
340 t3 = -125000.0_dp *year_sec
341 t4 = -21000.0_dp *year_sec
342 t5 = -8000.0_dp *year_sec
343 t6 = 0.0_dp *year_sec
344 
345 if (time < t1) then
346  z_sl = 0.0_dp
347 else if (time < t2) then
348  z_sl = z_sl_min*(time-t1)/(t2-t1)
349 else if (time < t3) then
350  z_sl = -z_sl_min*(time-t3)/(t3-t2)
351 else if (time < t4) then
352  z_sl = z_sl_min*(time-t3)/(t4-t3)
353 else if (time < t5) then
354  z_sl = -z_sl_min*(time-t5)/(t5-t4)
355 else if (time < t6) then
356  z_sl = 0.0_dp
357 else
358  z_sl = 0.0_dp
359 end if
360 
361 #elif (SEA_LEVEL==3)
362 
363 ! ------ z_sl from the SPECMAP record
364 
365 if (time/year_sec < real(specmap_time_min,dp)) then
366  z_sl = specmap_zsl(0)
367 else if (time/year_sec < real(specmap_time_max,dp)) then
368 
369 #if (!defined(ALLOW_GRDCHK) && !defined(ALLOW_OPENAD)) /* Normal */
370  i_kl = floor(((time/year_sec) &
371  -real(specmap_time_min,dp))/real(specmap_time_stp,dp))
372  i_kl = max(i_kl, 0)
373 
374  i_gr = ceiling(((time/year_sec) &
375  -real(specmap_time_min,dp))/real(specmap_time_stp,dp))
376  i_gr = min(i_gr, ndata_specmap)
377 #else /* OpenAD */
378  call myfloor((((time/year_sec) &
379  -real(specmap_time_min,dp))/real(specmap_time_stp,dp)),i_kl)
380  i_kl = max(i_kl, 0)
381 
382  call myceiling((((time/year_sec) &
383  -real(specmap_time_min,dp))/real(specmap_time_stp,dp)),i_gr)
384  i_gr = min(i_gr, ndata_specmap)
385 #endif /* Normal vs. OpenAD */
386 
387  if (i_kl == i_gr) then
388 
389  z_sl = specmap_zsl(i_kl)
390 
391  else
392 
393  time_kl = (specmap_time_min + i_kl*specmap_time_stp) *year_sec
394  time_gr = (specmap_time_min + i_gr*specmap_time_stp) *year_sec
395 
396  z_sl = specmap_zsl(i_kl) &
397  +(specmap_zsl(i_gr)-specmap_zsl(i_kl)) &
398  *(time-time_kl)/(time_gr-time_kl)
399  ! linear interpolation of the sea-level data
400 
401  end if
402 
403 else
404  z_sl = specmap_zsl(ndata_specmap)
405 end if
406 
407 #endif
408 
409 ! ------ Time derivative of the sea level
410 
411 if ( z_sl_old > -999999.9_dp ) then
412  dzsl_dtau = (z_sl-z_sl_old)/dtime
413 else ! only dummy value for z_sl_old available
414  dzsl_dtau = 0.0_dp
415 end if
416 
417 ! ------ Minimum bedrock elevation for extent of marine ice
418 
419 #if (MARGIN==2)
420 
421 #if (MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3)
422 z_mar = z_mar
423 #elif (MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5)
424 z_mar = fact_z_mar*z_sl
425 #elif (MARINE_ICE_CALVING==6 || MARINE_ICE_CALVING==7)
426 if (z_sl >= -80.0_dp) then
427  z_mar = 2.5_dp*z_sl
428 else
429  z_mar = 10.25_dp*(z_sl+80.0_dp)-200.0_dp
430 end if
431 z_mar = fact_z_mar*z_mar
432 #endif
433 
434 #endif
435 
436 ! ------ Update of the mask according to the sea level
437 
438 ! ---- Check all sea and floating-ice points and their direct
439 ! neighbours
440 
441 do i=0, imax
442 do j=0, jmax
443  check_point(j,i) = .false.
444 end do
445 end do
446 
447 do i=1, imax-1
448 do j=1, jmax-1
449  if (maske(j,i) >= 2_i1b) then
450  check_point(j ,i ) = .true.
451  check_point(j ,i+1) = .true.
452  check_point(j ,i-1) = .true.
453  check_point(j+1,i ) = .true.
454  check_point(j-1,i ) = .true.
455  end if
456 end do
457 end do
458 
459 do i=1, imax-1
460 do j=1, jmax-1
461  if (check_point(j,i)) then
462  maske_neu(j,i) = mask_update_sea_level(z_sl, i, j)
463  end if
464 end do
465 end do
466 
467 ! ---- Assign new values of the mask
468 
469 do i=1, imax-1
470 do j=1, jmax-1
471  if (check_point(j,i)) then
472  maske(j,i) = maske_neu(j,i)
473  end if
474 end do
475 end do
476 
477 !-------- Surface air temperatures --------
478 
479 #if (TEMP_PRESENT_PARA==1) /* Parameterisation by Fortuin and Oerlemans */
480  ! (1990) for the whole ice sheet
481 
482 #if (defined(TEMP_PRESENT_OFFSET))
483 theta_ma_offset = temp_present_offset
484 theta_mj_offset = temp_present_offset
485 #else
486 theta_ma_offset = 0.0_dp
487 theta_mj_offset = 0.0_dp
488 #endif
489 
490 theta_ma = 34.461_dp + theta_ma_offset
491 gamma_ma = -9.14e-03_dp
492 c_ma = -0.688_dp
493 
494 theta_mj = 16.81_dp + theta_mj_offset
495 gamma_mj = -6.92e-03_dp
496 c_mj = -0.27973_dp
497 
498 #elif (TEMP_PRESENT_PARA==2) /* Parameterisation by Fortuin and Oerlemans */
499  ! (1990), separately for three different
500  ! elevation ranges
501 
502 #if (defined(TEMP_PRESENT_OFFSET))
503 theta_ma_offset = temp_present_offset
504 theta_mj_offset = temp_present_offset
505 #else
506 theta_ma_offset = 0.0_dp
507 theta_mj_offset = 0.0_dp
508 #endif
509 
510 zs_sep_1 = 200.0_dp
511 zs_sep_2 = 1500.0_dp
512 
513 theta_ma_1 = 49.642_dp + theta_ma_offset
514 gamma_ma_1 = 0.0_dp
515 c_ma_1 = -0.943_dp
516 
517 theta_ma_2 = 36.689_dp + theta_ma_offset
518 gamma_ma_2 = -5.102e-03_dp
519 c_ma_2 = -0.725_dp
520 
521 theta_ma_3 = 7.405_dp + theta_ma_offset
522 gamma_ma_3 = -14.285e-03_dp
523 c_ma_3 = -0.180_dp
524 
525 theta_mj = 16.81_dp + theta_mj_offset
526 gamma_mj = -6.92e-03_dp
527 c_mj = -0.27937_dp
528 
529 #else
530 
531 errormsg = ' >>> boundary: Parameter TEMP_PRESENT_PARA must be either 1 or 2!'
532 call error(errormsg)
533 
534 #endif
535 
536 do i=0, imax
537 do j=0, jmax
538 
539 ! ------ Present-day mean-annual air temperature
540 
541 #if (TEMP_PRESENT_PARA==1)
542  temp_ma_present(j,i) = theta_ma + gamma_ma*zs(j,i) &
543  + c_ma*abs(phi(j,i))*pi_180_inv
544 #elif (TEMP_PRESENT_PARA==2)
545  if ( zs(j,i) <= zs_sep_1 ) then
546  temp_ma_present(j,i) = theta_ma_1 + gamma_ma_1*zs(j,i) &
547  + c_ma_1*abs(phi(j,i))*pi_180_inv
548  else if ( zs(j,i) <= zs_sep_2 ) then
549  temp_ma_present(j,i) = theta_ma_2 + gamma_ma_2*zs(j,i) &
550  + c_ma_2*abs(phi(j,i))*pi_180_inv
551  else
552  temp_ma_present(j,i) = theta_ma_3 + gamma_ma_3*zs(j,i) &
553  + c_ma_3*abs(phi(j,i))*pi_180_inv
554  end if
555 #endif
556 
557 ! ------ Present-day mean-January (summer) air temperature
558 
559  temp_mj_present(j,i) = theta_mj + gamma_mj*zs(j,i) &
560  + c_mj*abs(phi(j,i))*pi_180_inv
561 
562 #if (TSURFACE<=4)
563 
564 ! ------ Correction with deviation delta_ts
565 
566  temp_ma(j,i) = temp_ma_present(j,i) + delta_ts
567  temp_mm(j,i,7) = temp_mj_present(j,i) + delta_ts
568 
569 #elif (TSURFACE==5)
570 
571 ! ------ Correction with LGM anomaly and glacial index
572 
573  temp_ma(j,i) = temp_ma_present(j,i) + glac_index*temp_ma_lgm_anom(j,i)
574  temp_mm(j,i,7) = temp_mj_present(j,i) + glac_index*temp_mj_lgm_anom(j,i)
575 
576 #elif (TSURFACE==6)
577 
578 ! ------ Mean-annual and mean-January (summer) surface temperatures
579 ! from data read above --------
580 
581  temp_ma(j,i) = temp_ma_present(j,i) + temp_ma_anom_fact*temp_ma_anom(j,i)
582  temp_mm(j,i,7) = temp_mj_present(j,i) + temp_mj_anom_fact*temp_mj_anom(j,i)
583 
584 #endif
585 
586 ! ------ Amplitude of the annual cycle
587 
588  temp_ampl(j,i) = temp_mm(j,i,7) - temp_ma(j,i)
589 
590  if (temp_ampl(j,i) < eps) then
591  temp_ampl(j,i) = eps ! Correction of amplitude, if required
592  end if
593 
594 end do
595 end do
596 
597 ! ------ Monthly temperatures
598 
599 do n=1, 12 ! month counter
600 
601  sine_factor = sin((real(n,dp)-4.0_dp)*pi/6.0_dp)
602 
603  do i=0, imax
604  do j=0, jmax
605  temp_mm(j,i,n) = temp_ma(j,i) + sine_factor*temp_ampl(j,i)
606  end do
607  end do
608 
609 end do
610 
611 !-------- Accumulation-ablation function as_perp --------
612 
613 #if (INITMIP_CONST_SMB==1) /* SMB held constant at its initial distribution */
614 if (firstcall) then
615 #endif
616 
617 #if (ACCSURFACE<=3)
618 
619 #if (ELEV_DESERT==1)
620 
621 gamma_p = gamma_p*1.0e-03_dp ! Precipitation lapse rate
622  ! for elevation desertification, in m^(-1)
623 zs_thresh = zs_thresh ! Elevation threshold, in m
624 
625 #endif
626 
627 #elif (ACCSURFACE==4)
628 
629 alpha_p = 22.47_dp
630 beta_p = 0.046_dp
631 temp_0 = 273.15_dp
632 alpha_t = 0.67_dp
633 beta_t = 88.9_dp
634 
635 #endif
636 
637 #if (SOLID_PRECIP==1) /* Marsiat (1994) */
638 
639 temp_rain = 7.0_dp ! Threshold monthly mean temperature for
640  ! precipitation = 100% rain, in deg C
641 temp_snow = -10.0_dp ! Threshold monthly mean temperature for &
642  ! precipitation = 100% snow, in deg C
643 
644 inv_delta_temp_rain_snow = 1.0_dp/(temp_rain-temp_snow)
645 
646 #elif (SOLID_PRECIP==2) /* Bales et al. (2009) */
647 
648 temp_rain = 7.2_dp ! Threshold monthly mean temperature for &
649  ! precipitation = 100% rain, in deg C
650 temp_snow = -11.6_dp ! Threshold monthly mean temperature for &
651  ! precipitation = 100% snow, in deg C
652 
653 coeff(0) = 5.4714e-01_dp ! Coefficients
654 coeff(1) = -9.1603e-02_dp ! of
655 coeff(2) = -3.314e-03_dp ! the
656 coeff(3) = 4.66e-04_dp ! fifth-order
657 coeff(4) = 3.8e-05_dp ! polynomial
658 coeff(5) = 6.0e-07_dp ! fit
659 
660 #elif (SOLID_PRECIP==3) /* Huybrechts and de Wolde (1999) */
661 
662 temp_rain = 2.0_dp ! Threshold instantaneous temperature for &
663  ! precipitation = 100% rain, in deg C
664 temp_snow = temp_rain ! Threshold instantaneous temperature for &
665  ! precipitation = 100% snow, in deg C
666 
667 s_stat = s_stat_0 ! Standard deviation of the air termperature
668  ! (same parameter as in the PDD model)
669 
670 inv_sqrt2_s_stat = 1.0_dp/(sqrt(2.0_dp)*s_stat)
671 
672 #endif
673 
674 #if (ABLSURFACE==1 || ABLSURFACE==2)
675 
676 s_stat = s_stat_0
677 beta1 = beta1_0 *(0.001_dp/86400.0_dp)*(rho_w/rho)
678  ! (mm WE)/(d*deg C) --> (m IE)/(s*deg C)
679 beta2 = beta2_0 *(0.001_dp/86400.0_dp)*(rho_w/rho)
680  ! (mm WE)/(d*deg C) --> (m IE)/(s*deg C)
681 pmax = pmax_0
682 mu = mu_0 *(1000.0_dp*86400.0_dp)*(rho/rho_w)
683  ! (d*deg C)/(mm WE) --> (s*deg C)/(m IE)
684 
685 #elif (ABLSURFACE==3)
686 
687 lambda_lti = lambda_lti *(0.001_dp/year_sec)*(rho_w/rho)
688  ! (mm WE)/(a*deg C) --> (m IE)/(s*deg C)
689 temp_lti = temp_lti
690 mu = 0.0_dp ! no superimposed ice considered
691 
692 #endif
693 
694 do i=0, imax
695 do j=0, jmax
696 
697 ! ------ Accumulation
698 
699 #if (ACCSURFACE<=3)
700 
701 ! ---- Elevation desertification of precipitation
702 
703 #if (ELEV_DESERT==0)
704 
705  precip_fact = 1.0_dp ! no elevation desertification
706 
707 #elif (ELEV_DESERT==1)
708 
709  if (zs_ref(j,i) < zs_thresh) then
710  precip_fact &
711  = exp(gamma_p*(max(zs(j,i),zs_thresh)-zs_thresh))
712  else
713  precip_fact &
714  = exp(gamma_p*(max(zs(j,i),zs_thresh)-zs_ref(j,i)))
715  end if
716 
717 #else
718  errormsg = ' >>> boundary: Parameter ELEV_DESERT must be either 0 or 1!'
719  call error(errormsg)
720 #endif
721 
722  do n=1, 12 ! month counter
723  precip(j,i,n) = precip_present(j,i,n)*precip_fact
724  end do
725 
726 #endif
727 
728 ! ---- Precipitation change related to changing climate
729 
730 #if (ACCSURFACE==1)
731  precip_fact = accfact
732 #elif (ACCSURFACE==2)
733  precip_fact = 1.0_dp + gamma_s*delta_ts
734 #elif (ACCSURFACE==3)
735  precip_fact = exp(gamma_s*delta_ts)
736 #endif
737 
738 #if (ACCSURFACE<=3)
739 
740  precip(j,i,0) = 0.0_dp ! initialisation value for mean annual precip
741 
742  do n=1, 12 ! month counter
743  precip(j,i,n) = precip(j,i,n)*precip_fact ! monthly precip
744  precip(j,i,0) = precip(j,i,0) + precip(j,i,n)*inv_twelve
745  ! mean annual precip
746  end do
747 
748 #elif (ACCSURFACE==4)
749 
750  precip(j,i,0) = 0.0_dp ! initialisation value for mean annual precip
751 
752  temp_inv = alpha_t * (temp_ma(j,i)+temp_0) + beta_t ! in K
753  temp_inv_present = alpha_t * (temp_ma_present(j,i)+temp_0) + beta_t ! in K
754 
755  precip_fact = exp(alpha_p*(temp_0/temp_inv_present-temp_0/temp_inv)) &
756  *(temp_inv_present/temp_inv)**2 &
757  *(1.0_dp+beta_p*(temp_inv-temp_inv_present))
758 
759  do n=1, 12 ! month counter
760  precip(j,i,n) = precip_present(j,i,n)*precip_fact ! monthly precip
761  precip(j,i,0) = precip(j,i,0) + precip(j,i,n)*inv_twelve
762  ! mean annual precip
763  end do
764 
765 #elif (ACCSURFACE==5)
766 
767  precip(j,i,0) = 0.0_dp ! initialisation value for mean annual precip
768 
769  do n=1, 12 ! month counter
770 
771 #if (PRECIP_ANOM_INTERPOL==1)
772  precip_fact = 1.0_dp-glac_index+glac_index*precip_lgm_anom(j,i,n)
773  ! interpolation with a linear function
774 #elif (PRECIP_ANOM_INTERPOL==2)
775  precip_fact = exp(-glac_index*gamma_precip_lgm_anom(j,i,n))
776  ! interpolation with an exponential function
777 #endif
778 
779  precip(j,i,n) = precip_present(j,i,n)*precip_fact ! monthly precip
780  precip(j,i,0) = precip(j,i,0) + precip(j,i,n)*inv_twelve
781  ! mean annual precip
782  end do
783 
784 #elif (ACCSURFACE==6)
785 
786 ! -- Mean-annual precipitation from data already read above
787 
788  precip(j,i,0) = precip_ma_present(j,i) + precip_anom_fact*precip_ma_anom(j,i)
789  ! mean annual precip
790 
791  do n=1, 12 ! month counter
792  precip(j,i,n) = precip(j,i,0) ! monthly precip, assumed to be equal
793  ! to the mean annual precip
794  end do
795 
796 #endif
797 
798 ! ---- Annual accumulation, snowfall and rainfall rates
799 
800  accum(j,i) = precip(j,i,0)
801 
802  snowfall(j,i) = 0.0_dp ! initialisation value
803 
804  do n=1, 12 ! month counter
805 
806 #if (SOLID_PRECIP==1) /* Marsiat (1994) */
807 
808  if (temp_mm(j,i,n) >= temp_rain) then
809  frac_solid = 0.0_dp
810  else if (temp_mm(j,i,n) <= temp_snow) then
811  frac_solid = 1.0_dp
812  else
813  frac_solid = (temp_rain-temp_mm(j,i,n))*inv_delta_temp_rain_snow
814  end if
815 
816 #elif (SOLID_PRECIP==2) /* Bales et al. (2009) */
817 
818  if (temp_mm(j,i,n) >= temp_rain) then
819  frac_solid = 0.0_dp
820  else if (temp_mm(j,i,n) <= temp_snow) then
821  frac_solid = 1.0_dp
822  else
823  frac_solid = coeff(0) + temp_mm(j,i,n) * ( coeff(1) &
824  + temp_mm(j,i,n) * ( coeff(2) &
825  + temp_mm(j,i,n) * ( coeff(3) &
826  + temp_mm(j,i,n) * ( coeff(4) &
827  + temp_mm(j,i,n) * coeff(5) ) ) ) )
828  ! evaluation of 5th-order polynomial by Horner scheme
829  end if
830 
831 #elif (SOLID_PRECIP==3) /* Huybrechts and de Wolde (1999) */
832 
833 #if (!defined(ALLOW_GRDCHK) && !defined(ALLOW_OPENAD)) /* Normal */
834  frac_solid = 1.0_dp &
835  - 0.5_dp*erfc((temp_rain-temp_mm(j,i,n))*inv_sqrt2_s_stat)
836 #else /* OpenAD */
837  call my_erfc(((temp_rain-temp_mm(j,i,n))*inv_sqrt2_s_stat), temp_val)
838  frac_solid = 1.0_dp - 0.5_dp*temp_val
839 #endif /* Normal vs. OpenAD */
840 
841 #endif
842 
843  snowfall(j,i) = snowfall(j,i) + precip(j,i,n)*frac_solid*inv_twelve
844 
845  end do
846 
847  rainfall(j,i) = precip(j,i,0) - snowfall(j,i)
848 
849  if (snowfall(j,i) < 0.0_dp) snowfall(j,i) = 0.0_dp ! correction of
850  if (rainfall(j,i) < 0.0_dp) rainfall(j,i) = 0.0_dp ! negative values
851 
852 ! ------ Ablation
853 
854 ! ---- Runoff
855 
856 #if (ABLSURFACE==1 || ABLSURFACE==2)
857 
858 ! -- Temperature excess ET
859 
860  do n=1, 12 ! month counter
861  temp_mm_help(n) = temp_mm(j,i,n)
862  end do
863 
864  call pdd(temp_mm_help, s_stat, et(j,i))
865 
866 ! -- Formation rate of superimposed ice (melt_star), melt rate (melt)
867 ! and runoff rate (runoff)
868 
869 #if (ABLSURFACE==1)
870 
871  if ((beta1*et(j,i)) <= (pmax*snowfall(j,i))) then
872  melt_star(j,i) = beta1*et(j,i)
873  melt(j,i) = 0.0_dp
874  runoff(j,i) = melt(j,i)+rainfall(j,i)
875  else
876  melt_star(j,i) = pmax*snowfall(j,i)
877  melt(j,i) = beta2*(et(j,i)-melt_star(j,i)/beta1)
878  runoff(j,i) = melt(j,i)+rainfall(j,i)
879  end if
880 
881 #elif (ABLSURFACE==2)
882 
883  if ( rainfall(j,i) <= (pmax*snowfall(j,i)) ) then
884 
885  if ( (rainfall(j,i)+beta1*et(j,i)) <= (pmax*snowfall(j,i)) ) then
886  melt_star(j,i) = rainfall(j,i)+beta1*et(j,i)
887  melt(j,i) = 0.0_dp
888  runoff(j,i) = melt(j,i)
889  else
890  melt_star(j,i) = pmax*snowfall(j,i)
891  melt(j,i) = beta2 &
892  *(et(j,i)-(melt_star(j,i)-rainfall(j,i))/beta1)
893  runoff(j,i) = melt(j,i)
894  end if
895 
896  else
897 
898  melt_star(j,i) = pmax*snowfall(j,i)
899  melt(j,i) = beta2*et(j,i)
900  runoff(j,i) = melt(j,i) + rainfall(j,i)-pmax*snowfall(j,i)
901 
902  end if
903 
904 #endif
905 
906 #elif (ABLSURFACE==3)
907 
908  temp_jja_help = one_third*(temp_mm(j,i,6)+temp_mm(j,i,7)+temp_mm(j,i,8))
909 
910  melt_star(j,i) = 0.0_dp ! no superimposed ice considered
911  melt(j,i) = lambda_lti*max((temp_jja_help-temp_lti), 0.0_dp)
912  runoff(j,i) = melt(j,i) + rainfall(j,i)
913 
914 #endif
915 
916 ! ---- Evaporation
917 
918  evap(j,i) = 0.0_dp
919 
920 end do
921 end do
922 
923 #if (INITMIP_CONST_SMB==1) /* SMB held constant at its initial distribution */
924 
925 accum_save = accum
926 evap_save = evap
927 runoff_save = runoff
928 
929 else ! firstcall == .false.
930 
931 accum = accum_save
932 evap = evap_save
933 runoff = runoff_save
934 
935 end if ! firstcall == .true./.false.
936 
937 #endif
938 
939 ! ------ SMB = accumulation minus ablation
940 
941 as_perp = accum - evap - runoff
942 
943 ! ---- Prescribed SMB correction
944 
946 
947 #if (defined(INITMIP_SMB_ANOM_FILE)) /* Correction for ISMIP InitMIP */
948 
949 #if (defined(ALLOW_GRDCHK) || defined(ALLOW_OPENAD)) /* OpenAD */
950 call myfloor(time_in_years, i_time_in_years)
951 #endif /* OpenAD */
952 
953 if ((time_in_years > 0.0_dp).and.(time_in_years <= 40.0_dp)) then
954 #if (!defined(ALLOW_GRDCHK) && !defined(ALLOW_OPENAD)) /* Normal */
956  + 0.025_dp*floor(time_in_years) * smb_anom_initmip
957 #else /* OpenAD */
959  + 0.025_dp*(i_time_in_years) * smb_anom_initmip
960 #endif /* Normal vs. OpenAD */
961 else if (time_in_years > 40.0_dp) then
963  + smb_anom_initmip
964 end if
965 
966 #endif
967 
969 
970 ! ------ Ice-surface temperature (10-m firn temperature) temp_s,
971 ! including empirical firn-warming correction due to
972 ! refreezing meltwater when superimposed ice is formed
973 
974 #if (!defined(ALLOW_OPENAD) && !defined(ALLOW_GRDCHK)) /* Normal */
975 
976 where (melt_star >= melt)
977  temp_s = temp_ma + mu*(melt_star-melt)
978 elsewhere
979  temp_s = temp_ma
980 end where
981 
982 where (temp_s > -0.001_dp) temp_s = -0.001_dp
983  ! Cut-off of positive air temperatures
984 
985 #else /* OpenAD */
986 
987 do i=0, imax
988 do j=0, jmax
989  if (melt_star(j,i) >= melt(j,i)) then
990  temp_s(j,i) = temp_ma(j,i) + mu*(melt_star(j,i)-melt(j,i))
991  else
992  temp_s(j,i) = temp_ma(j,i)
993  end if
994 end do
995 end do
996 
997 do i=0, imax
998 do j=0, jmax
999  if (temp_s(j,i) > -0.001_dp) then
1000  temp_s(j,i) = -0.001_dp ! Cut-off of positive air temperatures
1001  end if
1002 end do
1003 end do
1004 
1005 #endif /* Normal vs. OpenAD */
1006 
1007 !-------- Calving rate of grounded ice --------
1008 
1009 calv_grounded = 0.0_dp
1010 
1011 #if ((MARGIN==2) \
1012  && (marine_ice_formation==2) \
1013  && (marine_ice_calving==9))
1014 
1015 call calving_underwater_ice(z_sl)
1017 
1018 #endif
1019 
1020 if (firstcall) firstcall = .false.
1021 
1022 end subroutine boundary
1023 
1024 !-------------------------------------------------------------------------------
1025 
1026 end module boundary_m
1027 !
real(dp), dimension(0:jmax, 0:imax) phi
phi(j,i): Geographic latitude of grid point (i,j)
integer(i4b) ndata_specmap
ndata_specmap: Number of data values for the sea level
Definition: ctrl_m.F90:9
real(dp) rho_w
RHO_W: Density of pure water.
subroutine error(error_message)
Main routine of error_m: Writing of error messages and stopping execution.
Definition: error_m.F90:47
integer(i4b) specmap_time_max
specmap_time_max: Maximum time of the data values for the sea level
real(dp) temp_precip_time_stp
temp_precip_time_stp: Time step of the surface-temperature and precipitation data ...
real(dp) mu_0
MU_0: Firn-warming correction.
real(dp), parameter eps
eps: Small number
real(dp) temp_precip_time_max
temp_precip_time_max: Maximum time of the surface-temperature and precipitation data ...
real(dp), dimension(0:jmax, 0:imax) calv_grounded
calv_grounded(j,i): Calving rate of grounded ice
real(dp), dimension(0:jmax, 0:imax) rainfall
rainfall(j,i): Rainfall rate at the ice surface
Several mathematical tools used by SICOPOLIS.
real(dp), dimension(0:jmax, 0:imax) temp_mj_present
temp_mj_present(j,i): Present-day mean summer (northern hemisphere: July, southern hemisphere: Januar...
integer(i4b) specmap_time_min
specmap_time_min: Minimum time of the data values for the sea level
Computation of the positive degree days (PDD) with statistical temperature fluctuations; based on sem...
Definition: pdd_m.F90:37
integer(i1b), dimension(0:jmax, 0:imax) maske_neu
maske_neu(j,i): New value of maske computed during an integration step
real(dp), dimension(0:jmax, 0:imax) as_perp
as_perp(j,i): Accumulation-ablation function at the ice surface (SMB)
real(dp), dimension(:), allocatable glacial_index
glacial_index(n): Data values for the glacial index
real(dp) s_stat_0
S_STAT_0: Standard deviation of the air termperature for the degree-day model.
real(dp), dimension(0:jmax, 0:imax) runoff
runoff(j,i): Runoff rate at the ice surface
real(dp), dimension(0:jmax, 0:imax) melt_star
melt_star(j,i): Superimposed ice formation rate at the ice surface
real(dp), dimension(0:jmax, 0:imax) smb_corr_in
smb_corr_in(j,i): Prescribed SMB correction
Definition: sico_vars_m.F90:43
real(dp), dimension(0:jmax, 0:imax) zs_ref
zs_ref(j,i): Reference elevation for precip_present, temp_ma_present and temp_mj_present ...
real(dp), dimension(0:jmax, 0:imax) snowfall
snowfall(j,i): Snowfall rate at the ice surface
real(dp), parameter pi_180_inv
pi_180_inv: 180 divided by pi (-> rad to deg)
real(dp), dimension(0:jmax, 0:imax) temp_s
temp_s(j,i): Ice surface temperature
real(dp), dimension(0:jmax, 0:imax) melt
melt(j,i): Melting rate at the ice surface
integer(i4b) grip_time_max
grip_time_max: Maximum time of the data values for the surface temperature anomaly ...
Declarations of global variables for SICOPOLIS (for the ANT domain).
Definition: sico_vars_m.F90:35
integer(i4b) ndata_gi
ndata_gi: Number of data values for the glacial index
integer(i4b) ndata_grip
ndata_grip: Number of data values for the surface temperature anomaly
real(dp) beta1_0
BETA1_0: Degree-day factor for snow.
real(dp), dimension(0:jmax, 0:imax) et
ET(j,i): Temperature excess at the ice surface (positive degree days divided by time) ...
real(dp), dimension(0:jmax, 0:imax) calv_uw_ice
Calving of "underwater ice".
Declarations of kind types for SICOPOLIS.
real(dp), dimension(0:jmax, 0:imax) zs
zs(j,i): Coordinate z of the surface topography
real(dp) pmax_0
PMAX_0: Saturation factor for the formation of superimposed ice.
real(dp), dimension(0:jmax, 0:imax) evap
evap(j,i): Evaporation rate at the ice surface
real(dp), dimension(0:jmax, 0:imax) accum
accum(j,i): Accumulation rate at the ice surface (includes liquid precipitation = rainfall!) ...
real(dp), dimension(0:jmax, 0:imax) precip_ma_present
precip_ma_present(j,i): Present-day mean annual precipitation rate at the ice surface ...
integer(i1b) function, public mask_update_sea_level(z_sl, i, j)
Main function of mask_update_m: Update of the ice-land-ocean mask due to changes of the sea level...
integer(i4b) gi_time_min
gi_time_min: Minimum time of the data values for the glacial index
real(dp) temp_precip_time_min
temp_precip_time_min: Minimum time of the surface-temperature and precipitation data ...
real(dp), dimension(:), allocatable specmap_zsl
specmap_zsl(n): Data values for the sea level
real(dp) beta2_0
BETA2_0: Degree-day factor for ice.
Writing of error messages and stopping execution.
Definition: error_m.F90:35
integer(i4b) ndata_temp_precip
ndata_temp_precip: Number of surface-temperature and precipitation data
real(dp), dimension(0:jmax, 0:imax) temp_ma_present
temp_ma_present(j,i): Present-day mean annual surface temperature
integer(i4b) gi_time_max
gi_time_max: Maximum time of the data values for the glacial index
real(dp), dimension(0:jmax, 0:imax, 12) gamma_precip_lgm_anom
gamma_precip_lgm_anom(j,i,n): negative natural logarithm of precip_lgm_anom(j,i,n) ...
NetCDF error capturing.
Definition: nc_check_m.F90:35
real(dp), dimension(0:jmax, 0:imax) temp_mj_lgm_anom
temp_mj_lgm_anom(j,i): LGM anomaly (difference LGM - present) of the mean summer (northern hemisphere...
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 ...
subroutine boundary(time, dtime, dxi, deta, delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
Main routine of boundary_m: Computation of the surface temperature (must be less than 0 deg C!) and o...
Definition: boundary_m.F90:61
real(dp), dimension(0:jmax, 0:imax) temp_ma_lgm_anom
temp_ma_lgm_anom(j,i): LGM anomaly (difference LGM - present) of the mean annual surface temperature ...
subroutine check(status, ch_calling_routine)
NetCDF error capturing.
Definition: nc_check_m.F90:48
real(dp), parameter pi
pi: Constant pi
integer(i4b) ncid_temp_precip
ncid_temp_precip: ID of the NetCDF file containing the surface-temperature and precipitation data as ...
integer(i4b) grip_time_min
grip_time_min: Minimum time of the data values for the surface temperature anomaly ...
subroutine calving_underwater_ice(z_sl)
Main routine: Calving of "underwater ice".
Computation of the surface temperature (must be less than 0 deg C!) and of the accumulation-ablation ...
Definition: boundary_m.F90:37
subroutine, public pdd(temp_mm, s_stat, ET)
Main subroutine of pdd_m: Computation of the positive degree days (PDD) with statistical temperature ...
Definition: pdd_m.F90:57
character(len=256) errormsg
errormsg: Error-message string
real(dp) rho
RHO: Density of ice.
real(dp), dimension(0:jmax, 0:imax) smb_corr_prescribed
smb_corr_prescribed(j,i): Prescribed SMB correction
Update of the ice-land-ocean mask due to changes of the sea level.
integer, parameter dp
Double-precision reals.
real(dp), dimension(0:jmax, 0:imax, 12) precip_present
precip_present(j,i,n): Present-day mean monthly precipitation rate at the ice surface ...
real(dp), dimension(0:jmax, 0:imax, 12) precip_lgm_anom
precip_lgm_anom(j,i,n): LGM anomaly (ratio LGM/present) of the mean monthly precipitation rate at the...
Declarations of global variables for SICOPOLIS.
real(dp), dimension(:), allocatable griptemp
griptemp(n): Data values for the surface temperature anomaly