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