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