SICOPOLIS V5-dev  Revision 1420
init_temp_water_age_m.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Module : i n i t _ t e m p _ w a t e r _ a g e _ m
4 !
5 !> @file
6 !!
7 !! Initial temperature, water content and age.
8 !!
9 !! @section Copyright
10 !!
11 !! Copyright 2009-2019 Ralf Greve, Thorben Dunse
12 !!
13 !! @section License
14 !!
15 !! This file is part of SICOPOLIS.
16 !!
17 !! SICOPOLIS is free software: you can redistribute it and/or modify
18 !! it under the terms of the GNU General Public License as published by
19 !! the Free Software Foundation, either version 3 of the License, or
20 !! (at your option) any later version.
21 !!
22 !! SICOPOLIS is distributed in the hope that it will be useful,
23 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
24 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25 !! GNU General Public License for more details.
26 !!
27 !! You should have received a copy of the GNU General Public License
28 !! along with SICOPOLIS. If not, see <http://www.gnu.org/licenses/>.
29 !<
30 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
31 
32 !-------------------------------------------------------------------------------
33 !> Initial temperature, water content and age.
34 !<------------------------------------------------------------------------------
36 
37  use sico_types_m
39  use sico_vars_m
40 
41  implicit none
42 
43  private
46  public :: init_temp_water_age_1_5
47  public :: init_temp_water_age_2
48 
49 contains
50 
51 !-------------------------------------------------------------------------------
52 !> Initial temperature, water content and age
53 !! (case ANF_DAT==1, TEMP_INIT==1:
54 !! present-day initial topography, isothermal conditions).
55 !<------------------------------------------------------------------------------
56  subroutine init_temp_water_age_1_1()
57 
58  implicit none
59 
60  integer(i4b) :: i, j, kc
61 
62 !-------- Initial ice temperature --------
63 
64  do i=0, imax
65  do j=0, jmax
66 
67 #if (defined(NMARS) || defined(SMARS)) /* Polar caps of Mars */
68 
69  do kc=0, kcmax
70  temp_c(kc,j,i) = -100.0_dp
71  end do
72 
73 #else /* all other domains */
74 
75  do kc=0, kcmax
76  temp_c(kc,j,i) = -10.0_dp
77  end do
78 
79 #endif
80 
81  end do
82  end do
83 
84 !-------- Initial lithosphere temperature, water content and age --------
85 
86  call init_temp_r()
87  call init_water()
88  call init_age()
89 
90  end subroutine init_temp_water_age_1_1
91 
92 !-------------------------------------------------------------------------------
93 !> Initial temperature, water content and age
94 !! (case ANF_DAT==1, TEMP_INIT==2:
95 !! present-day initial topography,
96 !! ice temperature equal to local surface temperature).
97 !<------------------------------------------------------------------------------
98  subroutine init_temp_water_age_1_2()
99 
100  implicit none
101 
102  integer(i4b) :: i, j, kc
103 
104 !-------- Initial ice temperature --------
105 
106  do i=0, imax
107  do j=0, jmax
108 
109  do kc=0, kcmax
110  temp_c(kc,j,i) = temp_s(j,i)
111  end do
112 
113  end do
114  end do
115 
116 !-------- Initial lithosphere temperature, water content and age --------
117 
118  call init_temp_r()
119  call init_water()
120  call init_age()
121 
122  end subroutine init_temp_water_age_1_2
123 
124 !-------------------------------------------------------------------------------
125 !> Initial temperature, water content and age
126 !! (case ANF_DAT==1, TEMP_INIT==3:
127 !! present-day initial topography,
128 !! ice temperature linearly increasing with depth).
129 !<------------------------------------------------------------------------------
130  subroutine init_temp_water_age_1_3()
133 
134  implicit none
135 
136  integer(i4b) :: i, j, kc
137  real(dp) :: kappa_const_val
138  real(dp) :: temp_ice_base
139 
140 !-------- Initial ice temperature --------
141 
142 #if (defined(NMARS) || defined(SMARS)) /* Polar caps of Mars */
143  kappa_const_val = kappa_val(-100.0_dp)
144 #else /* all other domains */
145  kappa_const_val = kappa_val(-10.0_dp)
146 #endif
147 
148  do i=0, imax
149  do j=0, jmax
150 
151  if (maske(j,i)<=2_i1b) then
152 
153  do kc=0, kcmax
154 
155  temp_c(kc,j,i) = temp_s(j,i) &
156  + (q_geo(j,i)/kappa_const_val) &
157  *h_c(j,i)*(1.0_dp-eaz_c_quotient(kc))
158  ! linear temperature distribution according to the
159  ! geothermal heat flux
160  end do
161 
162  if (temp_c(0,j,i) >= -beta*h_c(j,i)) then
163 
164  temp_ice_base = -beta*h_c(j,i)
165 
166  do kc=0, kcmax
167  temp_c(kc,j,i) = temp_s(j,i) &
168  + (temp_ice_base-temp_s(j,i)) &
169  *(1.0_dp-eaz_c_quotient(kc))
170  end do
171 
172  end if
173 
174  else ! maske(j,i)==3_i1b, floating ice
175 
176  temp_ice_base = -beta*h_c(j,i) - delta_tm_sw
177 
178  do kc=0, kcmax
179  temp_c(kc,j,i) = temp_s(j,i) &
180  + (temp_ice_base-temp_s(j,i)) &
181  *(1.0_dp-eaz_c_quotient(kc))
182  end do
183 
184  end if
185 
186  end do
187  end do
188 
189 !-------- Initial lithosphere temperature, water content and age --------
190 
191  call init_temp_r()
192  call init_water()
193  call init_age()
194 
195  end subroutine init_temp_water_age_1_3
196 
197 !-------------------------------------------------------------------------------
198 !> Initial temperature, water content and age
199 !! (case ANF_DAT==1, TEMP_INIT==4:
200 !! present-day initial topography, ice temperature from Robin (1955) solution).
201 !<------------------------------------------------------------------------------
202  subroutine init_temp_water_age_1_4()
205 
206  implicit none
207 
208  integer(i4b) :: i, j, kc
209  real(dp) :: kappa_const_val, c_const_val
210  real(dp) :: as_val, H_val, qgeo_val, K, z_above_base
211  real(dp) :: erf_val_1, erf_val_2
212  real(dp) :: temp_ice_base, temp_scale_factor
213 
214 !-------- Initial ice temperature --------
215 
216 #if (defined(NMARS) || defined(SMARS)) /* Polar caps of Mars */
217  kappa_const_val = kappa_val(-100.0_dp)
218  c_const_val = c_val(-100.0_dp)
219 #else /* all other domains */
220  kappa_const_val = kappa_val(-10.0_dp)
221  c_const_val = c_val(-10.0_dp)
222 #endif
223 
224  do i=0, imax
225  do j=0, jmax
226 
227  if (maske(j,i)<=2_i1b) then
228  as_val = max(as_perp(j,i), epsi)
229  else ! maske(j,i)==3_i1b, floating ice
230  as_val = epsi ! this will produce an almost linear temperature profile
231  end if
232 
233  h_val = max(h_c(j,i) , eps)
234  qgeo_val = max(q_geo(j,i), eps)
235 
236  k = sqrt( (kappa_const_val/(rho*c_const_val)) * (h_val/as_val) )
237 
238  erf_val_1 = erf(h_c(j,i)/(sqrt(2.0_dp)*k))
239 
240  do kc=0, kcmax
241  z_above_base = h_c(j,i)*eaz_c_quotient(kc)
242  erf_val_2 = erf(z_above_base/(sqrt(2.0_dp)*k))
243  temp_c(kc,j,i) = temp_s(j,i) &
244  + (qgeo_val/kappa_const_val) &
245  * sqrt(0.5_dp*pi)*k*(erf_val_1-erf_val_2)
246  end do
247 
248  if ( (maske(j,i) <= 2_i1b).and.(temp_c(0,j,i) >= -beta*h_c(j,i)) ) then
249  temp_ice_base = -beta*h_c(j,i)
250  temp_scale_factor = (temp_ice_base-temp_s(j,i)) &
251  /(temp_c(0,j,i)-temp_s(j,i))
252  else if (maske(j,i) == 3_i1b) then
253  temp_ice_base = -beta*h_c(j,i)-delta_tm_sw
254  temp_scale_factor = (temp_ice_base-temp_s(j,i)) &
255  /(temp_c(0,j,i)-temp_s(j,i))
256  else
257  temp_scale_factor = 1.0_dp
258  end if
259 
260  do kc=0, kcmax
261  temp_c(kc,j,i) = temp_s(j,i) &
262  + temp_scale_factor*(temp_c(kc,j,i)-temp_s(j,i))
263  end do
264 
265  end do
266  end do
267 
268 !-------- Initial lithosphere temperature, water content and age --------
269 
270  call init_temp_r()
271  call init_water()
272  call init_age()
273 
274  end subroutine init_temp_water_age_1_4
275 
276 !-------------------------------------------------------------------------------
277 !> Initial temperature, water content and age
278 !! (case ANF_DAT==1, TEMP_INIT==5:
279 !! present-day initial topography,
280 !! ice temperature, water content and age from previous simulation).
281 !<------------------------------------------------------------------------------
282  subroutine init_temp_water_age_1_5(z_sl, filename)
284  use read_m, only : read_erg_nc
285 
286  implicit none
287 
288  character(len=100), intent(in) :: filename
289  real(dp), intent(inout) :: z_sl
290 
291  integer(i4b) :: i, j, kc, kt, kr
292  real(dp) :: temp_ice_base, temp_scale_factor
293 
294  integer(i1b), dimension(0:JMAX,0:IMAX) :: maske_aux, n_cts_aux
295  integer(i4b), dimension(0:JMAX,0:IMAX) :: kc_cts_aux
296  real(dp), dimension(0:JMAX,0:IMAX) :: H
297  real(dp), dimension(0:JMAX,0:IMAX) :: H_cold_aux, H_temp_aux, H_aux
298  real(dp), dimension(0:KRMAX,0:JMAX,0:IMAX) :: temp_r_aux
299  real(dp), dimension(0:KTMAX,0:JMAX,0:IMAX) :: omega_t_aux, age_t_aux
300  real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: temp_c_aux, age_c_aux
301  real(dp), dimension(0:KCMAX,0:JMAX,0:IMAX) :: omega_c_aux
302 
303  h = h_c + h_t
304 
306 
307  call read_erg_nc(z_sl, filename, &
308  opt_maske = maske_aux , &
309  opt_n_cts = n_cts_aux , &
310  opt_kc_cts = kc_cts_aux , &
311  opt_h_cold = h_cold_aux , &
312  opt_h_temp = h_temp_aux , &
313  opt_h = h_aux , &
314  opt_temp_r = temp_r_aux , &
315  opt_omega_t = omega_t_aux , &
316  opt_age_t = age_t_aux , &
317  opt_temp_c = temp_c_aux , &
318  opt_age_c = age_c_aux , &
319  opt_omega_c = omega_c_aux , &
320  opt_flag_temp_age_only = .true.)
321 
322  do i=0, imax
323  do j=0, jmax
324 
325  if ( ((maske(j,i)==0_i1b).or.(maske(j,i)==3_i1b)) &
326  .and. &
327  ((maske_aux(j,i)==0_i1b).or.(maske_aux(j,i)==3_i1b)) ) then
328 
329  n_cts(j,i) = n_cts_aux(j,i)
330  kc_cts(j,i) = kc_cts_aux(j,i)
331 
332 #if (CALCMOD==1)
333  h_t(j,i) = h(j,i) &
334  * min( &
335  max(h_temp_aux(j,i),0.0_dp)/max(h_aux(j,i),eps_dp), &
336  1.0_dp &
337  )
338  h_c(j,i) = h(j,i) - h_t(j,i)
339  zm(j,i) = zb(j,i) + h_t(j,i)
340 #endif
341 
342  temp_r(:,j,i) = temp_r_aux(:,j,i)
343 
344  omega_t(:,j,i) = omega_t_aux(:,j,i)
345  age_t(:,j,i) = age_t_aux(:,j,i)
346 
347  temp_c(:,j,i) = temp_c_aux(:,j,i)
348  age_c(:,j,i) = age_c_aux(:,j,i)
349  omega_c(:,j,i) = omega_c_aux(:,j,i)
350 
351  end if
352 
353  if ( (maske(j,i)==3_i1b).and.(maske_aux(j,i)==0_i1b) ) then
354  ! correction for ice shelves
355 
356  n_cts(j,i) = -1_i1b
357  kc_cts(j,i) = 0
358 
359  h_t(j,i) = 0.0_dp
360  h_c(j,i) = h(j,i)
361  zm(j,i) = zb(j,i)
362 
363  omega_t(:,j,i) = 0.0_dp
364  omega_c(:,j,i) = 0.0_dp
365 
366  temp_ice_base = -beta*h_c(j,i)-delta_tm_sw
367  temp_scale_factor = (temp_ice_base-temp_s(j,i)) &
368  /(temp_c(0,j,i)-temp_s(j,i))
369 
370  do kc=0, kcmax
371  temp_c(kc,j,i) = temp_s(j,i) &
372  + temp_scale_factor*(temp_c(kc,j,i)-temp_s(j,i))
373  end do
374 
375 
376  end if
377 
378  end do
379  end do
380 
381  end subroutine init_temp_water_age_1_5
382 
383 !-------------------------------------------------------------------------------
384 !> Initial temperature, water content and age
385 !! (case ANF_DAT==2: ice-free conditions with relaxed bedrock).
386 !<------------------------------------------------------------------------------
387  subroutine init_temp_water_age_2()
389  implicit none
390 
391  integer(i4b) :: i, j, kc
392 
393 !-------- Initial ice temperature --------
394 
395  do i=0, imax
396  do j=0, jmax
397 
398  do kc=0, kcmax
399  temp_c(kc,j,i) = temp_s(j,i)
400  end do
401 
402  end do
403  end do
404 
405 !-------- Initial lithosphere temperature, water content and age --------
406 
407  call init_temp_r()
408  call init_water()
409  call init_age()
410 
411  end subroutine init_temp_water_age_2
412 
413 !-------------------------------------------------------------------------------
414 !> Initial lithosphere temperature.
415 !<------------------------------------------------------------------------------
416  subroutine init_temp_r()
417 
418  implicit none
419 
420  integer(i4b) :: i, j, kr
421 
422  do i=0, imax
423  do j=0, jmax
424 
425  do kr=0, krmax
426  temp_r(kr,j,i) = temp_c(0,j,i) &
427  + (q_geo(j,i)/kappa_r) &
428  *h_r*(1.0_dp-zeta_r(kr))
429  ! linear temperature distribution according to the
430  ! geothermal heat flux
431  end do
432 
433  end do
434  end do
435 
436  end subroutine init_temp_r
437 
438 !-------------------------------------------------------------------------------
439 !> Initial water content.
440 !<------------------------------------------------------------------------------
441  subroutine init_water()
442 
443  implicit none
444 
445  omega_c = 0.0_dp ! only required for the enthalpy method
446  omega_t = 0.0_dp
447 
448  end subroutine init_water
449 
450 !-------------------------------------------------------------------------------
451 !> Initial age.
452 !<------------------------------------------------------------------------------
453  subroutine init_age()
454 
455  implicit none
456 
457 #if (defined(ASF)) /* Austfonna */
458 
459  age_c = 3500.0_dp*year_sec
460  age_t = 3500.0_dp*year_sec
461 
462 #elif (defined(NMARS) || defined(SMARS)) /* Polar caps of Mars */
463 
464  age_c = 1.0e+06_dp*year_sec
465  age_t = 1.0e+06_dp*year_sec
466 
467 #else /* all other domains */
468 
469  age_c = 15000.0_dp*year_sec
470  age_t = 15000.0_dp*year_sec
471 
472 #endif
473 
474  end subroutine init_age
475 
476 !-------------------------------------------------------------------------------
477 
478 end module init_temp_water_age_m
479 !
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) temp_c
temp_c(kc,j,i): Temperature in the upper (kc) ice domain
real(dp) kappa_r
KAPPA_R: Heat conductivity of the lithosphere.
real(dp) delta_tm_sw
DELTA_TM_SW: Melting point depression of sea water due to its average salinity.
integer(i1b), dimension(0:jmax, 0:imax) n_cts
n_cts(j,i): Mask for thermal conditions. -1: cold ice base, 0: temperate ice base with cold ice above...
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) age_c
age_c(kc,j,i): Age in the upper (kc) ice domain
real(dp), parameter epsi
epsi: Very small number
real(dp), dimension(0:jmax, 0:imax) q_geo
q_geo(j,i): Geothermal heat flux
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) age_t
age_t(kt,j,i): Age in the lower (kt) ice domain
real(dp), dimension(0:krmax, 0:jmax, 0:imax) temp_r
temp_r(kr,j,i): Temperature in the bedrock
real(dp) h_r
H_R: Thickness of the modelled lithosphere layer.
Material properties of ice: Rate factor, heat conductivity, specific heat (heat capacity), creep function, viscosity.
real(dp), parameter eps
eps: Small number
subroutine, public init_temp_water_age_2()
Initial temperature, water content and age (case ANF_DAT==2: ice-free conditions with relaxed bedrock...
real(dp) function, public kappa_val(temp_val)
Heat conductivity of ice: Linear interpolation of tabulated values in KAPPA(.).
subroutine, public read_erg_nc(z_sl, filename, opt_maske, opt_n_cts, opt_kc_cts, opt_H_cold, opt_H_temp, opt_H, opt_temp_r, opt_omega_t, opt_age_t, opt_temp_c, opt_age_c, opt_omega_c, opt_flag_temp_age_only)
Reading of time-slice files in native binary or NetCDF format.
Definition: read_m.F90:60
real(dp), dimension(0:jmax, 0:imax) as_perp
as_perp(j,i): Accumulation-ablation function at the ice surface (SMB)
subroutine, public init_temp_water_age_1_2()
Initial temperature, water content and age (case ANF_DAT==1, TEMP_INIT==2: present-day initial topogr...
real(dp), dimension(0:jmax, 0:imax) temp_s
temp_s(j,i): Ice surface temperature
Declarations of global variables for SICOPOLIS (for the ANT domain).
Definition: sico_vars_m.F90:35
Initial temperature, water content and age.
Reading of several input data.
Definition: read_m.F90:35
Declarations of kind types for SICOPOLIS.
integer(i4b), dimension(0:jmax, 0:imax) kc_cts
kc_cts(j,i): Position kc of the CTS (for COLD, ENTC, ENTM)
real(dp), dimension(0:jmax, 0:imax) h_c
H_c(j,i): Thickness of ice in the upper (kc) domain (thickness of the cold-ice layer for POLY...
real(dp), dimension(0:krmax) zeta_r
zeta_r(kr): Sigma coordinate zeta_r of grid point kr
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) omega_c
omega_c(kc,j,i): Water content in the upper (kc) ice domain
real(dp), parameter eps_dp
eps_dp: Small number to double-precision accuracy
real(dp) function, public c_val(temp_val)
Specific heat of ice: Linear interpolation of tabulated values in C(.).
real(dp), dimension(0:jmax, 0:imax) zm
zm(j,i): Coordinate z of the bottom of the upper (kc) ice domain = top of the lower (kt) ice domain (...
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 ...
real(dp), dimension(0:kcmax) eaz_c_quotient
eaz_c_quotient(kc): Abbreviation for (eaz_c(kc)-1.0)/(ea-1.0)
subroutine, public init_temp_water_age_1_3()
Initial temperature, water content and age (case ANF_DAT==1, TEMP_INIT==3: present-day initial topogr...
real(dp) beta
BETA: Clausius-Clapeyron gradient of ice.
real(dp), parameter pi
pi: Constant pi
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) omega_t
omega_t(kt,j,i): Water content in the lower (kt) ice domain
subroutine, public init_temp_water_age_1_5(z_sl, filename)
Initial temperature, water content and age (case ANF_DAT==1, TEMP_INIT==5: present-day initial topogr...
real(dp) rho
RHO: Density of ice.
subroutine, public init_temp_water_age_1_4()
Initial temperature, water content and age (case ANF_DAT==1, TEMP_INIT==4: present-day initial topogr...
subroutine, public init_temp_water_age_1_1()
Initial temperature, water content and age (case ANF_DAT==1, TEMP_INIT==1: present-day initial topogr...
real(dp), dimension(0:jmax, 0:imax) zb
zb(j,i): Coordinate z of the ice base
real(dp), dimension(0:jmax, 0:imax) h_t
H_t(j,i): Thickness of ice in the lower (kt) domain (thickness of the temperate layer for POLY...
Declarations of global variables for SICOPOLIS.