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