SICOPOLIS V5-dev  Revision 1288
read_m.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Module : r e a d _ m
4 !
5 !> @file
6 !!
7 !! Reading of several input data.
8 !!
9 !! @section Copyright
10 !!
11 !! Copyright 2009-2018 Ralf Greve, Fuyuki Saito
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 !> Reading of several input data.
34 !<------------------------------------------------------------------------------
35 module read_m
36 
37  use sico_types_m
38  use error_m
39 
40  implicit none
41 
42  private
44 
45 contains
46 
47 !-------------------------------------------------------------------------------
48 !> Reading of time-slice files in native binary or NetCDF format.
49 !<------------------------------------------------------------------------------
50  subroutine read_erg_nc(z_sl, filename)
51 
53  use sico_vars_m
54 
55 #if (NETCDF==2) /* time-slice file in NetCDF format */
56  use netcdf
57  use nc_check_m
58 #endif
59 
60  implicit none
61 
62  character(len=100), intent(in) :: filename
63 
64  real(dp), intent(out) :: z_sl
65 
66  integer(i4b) :: i, j, kc, kt, kr, n
67  character(len=256) :: filename_with_path
68 
69 #if (NETCDF==1) /* time-slice file in native binary format */
70 
71  integer(i4b) :: ios
72  character(len=256) :: ch_attr_title, ch_attr_institution, ch_attr_source, &
73  ch_attr_history, ch_attr_references
74 
75 #elif (NETCDF==2) /* time-slice file in NetCDF format */
76 
77  integer(i4b) :: ncid, ncv
78  ! ncid: ID of the output file
79  ! ncv: Variable ID
80 
81 #else
82 errormsg = ' >>> read_erg_nc: Parameter NETCDF must be either 1 or 2!'
83 call error(errormsg)
84 #endif
85 
86  integer(i2b), dimension(0:IMAX,0:JMAX) :: maske_conv, maske_old_conv, &
87  n_cts_conv, kc_cts_conv
88  real(sp) :: time_conv, dummy_conv, z_sl_conv, &
89  V_tot_conv, A_grounded_conv, A_floating_conv, &
90  H_R_conv, &
91  xi_conv(0:imax), eta_conv(0:jmax), &
92  sigma_level_c_conv(0:kcmax), sigma_level_t_conv(0:ktmax), &
93  sigma_level_r_conv(0:krmax)
94  real(sp), dimension(0:IMAX,0:JMAX) :: lambda_conv, phi_conv, &
95  lon_conv, lat_conv, &
96  temp_s_conv, accum_conv, &
97  as_perp_conv, as_perp_apl_conv, smb_corr_conv, &
98  q_geo_conv, &
99  zs_conv, zm_conv, zb_conv, zl_conv, zl0_conv, &
100  H_cold_conv, H_temp_conv, H_conv, &
101  Q_bm_conv, Q_tld_conv, &
102  am_perp_conv, &
103  qx_conv, qy_conv, &
104  dzs_dtau_conv, dzm_dtau_conv, dzb_dtau_conv, dzl_dtau_conv, &
105  dH_c_dtau_conv, dH_t_dtau_conv, dH_dtau_conv, &
106  vx_b_g_conv, vy_b_g_conv, vz_b_conv, vh_b_conv, &
107  vx_s_g_conv, vy_s_g_conv, vz_s_conv, vh_s_conv, &
108  vx_m_g_conv, vy_m_g_conv, vh_m_conv, &
109  temp_b_conv, temph_b_conv, &
110  tau_b_driving_conv, tau_b_drag_conv, &
111  p_b_w_conv, H_w_conv, q_gl_g_conv, q_cf_g_conv
112  real(sp), dimension(0:IMAX,0:JMAX,0:KCMAX) :: vx_c_conv, vy_c_conv, vz_c_conv, &
113  temp_c_conv, age_c_conv, &
114  enth_c_conv, omega_c_conv, &
115  enh_c_conv
116  real(sp), dimension(0:IMAX,0:JMAX,0:KTMAX) :: vx_t_conv, vy_t_conv, vz_t_conv, &
117  omega_t_conv, age_t_conv, &
118  enth_t_conv, &
119  enh_t_conv
120  real(sp), dimension(0:IMAX,0:JMAX,0:KRMAX) :: temp_r_conv
121 
122 #if (DISC>0) /* Ice discharge parameterisation */
123  integer(i2b), dimension(0:IMAX,0:JMAX) :: mask_mar_conv
124  real(sp), dimension(0:IMAX,0:JMAX) :: dis_perp_conv, &
125  cst_dist_conv, cos_grad_tc_conv
126 #endif
127 
128 !-------- Read data from time-slice file of previous simulation --------
129 
130  filename_with_path = trim(anfdatpath)//'/'//trim(filename)
131 
132 #if (NETCDF==1) /* time-slice file in native binary format */
133 
134 #ifndef ALLOW_OPENAD /* Normal */
135  open(unit=11, iostat=ios, file=trim(filename_with_path), status='old', &
136  form='unformatted')
137 #else /* OpenAD */
138  open(unit=11, iostat=ios, file=trim(filename_with_path))
139 #endif /* Normal vs. OpenAD */
140 
141  if (ios /= 0) then
142  errormsg = ' >>> read_erg_nc: Error when opening the time-slice file!'
143  call error(errormsg)
144  end if
145 
146  read(unit=11) ch_attr_title
147  read(unit=11) ch_attr_institution
148  read(unit=11) ch_attr_source
149  read(unit=11) ch_attr_history
150  read(unit=11) ch_attr_references
151 
152  read(unit=11) time_conv
153  read(unit=11) dummy_conv ! this is either delta_ts or glac_index; not needed
154  read(unit=11) z_sl_conv
155 
156  read(unit=11) dummy_conv ! this is V_tot; not needed
157  read(unit=11) dummy_conv ! this is V_af; not needed
158  read(unit=11) dummy_conv ! this is A_grounded; not needed
159  read(unit=11) dummy_conv ! this is A_floating; not needed
160 
161  read(unit=11) xi_conv
162  read(unit=11) eta_conv
163  read(unit=11) sigma_level_c_conv
164  read(unit=11) sigma_level_t_conv
165  read(unit=11) sigma_level_r_conv
166  read(unit=11) lon_conv
167  read(unit=11) lat_conv
168  read(unit=11) lambda_conv
169  read(unit=11) phi_conv
170  read(unit=11) temp_s_conv
171  read(unit=11) accum_conv
172  read(unit=11) as_perp_conv
173  read(unit=11) as_perp_apl_conv
174  read(unit=11) smb_corr_conv
175 
176 #if (DISC>0) /* Ice discharge parameterisation */
177 
178  read(unit=11) dis_perp_conv
179  read(unit=11) cst_dist_conv
180  read(unit=11) cos_grad_tc_conv
181  read(unit=11) mask_mar_conv
182 
183 #endif
184 
185  read(unit=11) q_geo_conv
186  read(unit=11) maske_conv
187  read(unit=11) maske_old_conv
188  read(unit=11) n_cts_conv
189  read(unit=11) kc_cts_conv
190  read(unit=11) zs_conv
191  read(unit=11) zm_conv
192  read(unit=11) zb_conv
193  read(unit=11) zl_conv
194  read(unit=11) zl0_conv
195  read(unit=11) h_cold_conv
196  read(unit=11) h_temp_conv
197  read(unit=11) h_conv
198  read(unit=11) h_r_conv
199  read(unit=11) q_bm_conv
200  read(unit=11) q_tld_conv
201  read(unit=11) am_perp_conv
202  read(unit=11) qx_conv
203  read(unit=11) qy_conv
204  read(unit=11) dzs_dtau_conv
205  read(unit=11) dzm_dtau_conv
206  read(unit=11) dzb_dtau_conv
207  read(unit=11) dzl_dtau_conv
208  read(unit=11) dh_c_dtau_conv
209  read(unit=11) dh_t_dtau_conv
210  read(unit=11) dh_dtau_conv
211  read(unit=11) vx_b_g_conv
212  read(unit=11) vy_b_g_conv
213  read(unit=11) vz_b_conv
214  read(unit=11) vh_b_conv
215  read(unit=11) vx_s_g_conv
216  read(unit=11) vy_s_g_conv
217  read(unit=11) vz_s_conv
218  read(unit=11) vh_s_conv
219  read(unit=11) vx_m_g_conv
220  read(unit=11) vy_m_g_conv
221  read(unit=11) vh_m_conv
222  read(unit=11) temp_b_conv
223  read(unit=11) temph_b_conv
224  read(unit=11) tau_b_driving_conv
225  read(unit=11) tau_b_drag_conv
226  read(unit=11) p_b_w_conv
227  read(unit=11) h_w_conv
228  read(unit=11) q_gl_g_conv
229  read(unit=11) q_cf_g_conv
230 
231  read(unit=11) vx_c_conv
232  read(unit=11) vy_c_conv
233  read(unit=11) vz_c_conv
234  read(unit=11) vx_t_conv
235  read(unit=11) vy_t_conv
236  read(unit=11) vz_t_conv
237  read(unit=11) temp_c_conv
238  read(unit=11) omega_t_conv
239  read(unit=11) temp_r_conv
240  read(unit=11) enth_c_conv
241  read(unit=11) enth_t_conv
242  read(unit=11) omega_c_conv
243  read(unit=11) enh_c_conv
244  read(unit=11) enh_t_conv
245  read(unit=11) age_c_conv
246  read(unit=11) age_t_conv
247 
248  close(unit=11,status='keep')
249 
250 #elif (NETCDF==2) /* time-slice file in NetCDF format */
251 
252  call check( nf90_open(trim(filename_with_path), nf90_nowrite, ncid) )
253 
254  call check( nf90_inq_varid(ncid, 'time', ncv) )
255  call check( nf90_get_var(ncid, ncv, time_conv) )
256 
257  if ( nf90_inq_varid(ncid, 'delta_ts', ncv) == nf90_noerr ) then
258  call check( nf90_get_var(ncid, ncv, dummy_conv) )
259  else if ( nf90_inq_varid(ncid, 'glac_index', ncv) == nf90_noerr ) then
260  call check( nf90_get_var(ncid, ncv, dummy_conv) )
261  else
262  dummy_conv = 0.0_sp
263  end if
264 
265  call check( nf90_inq_varid(ncid, 'z_sl', ncv) )
266  call check( nf90_get_var(ncid, ncv, z_sl_conv) )
267 
268  call check( nf90_inq_varid(ncid, 'x', ncv) )
269  call check( nf90_get_var(ncid, ncv, xi_conv) )
270 
271  call check( nf90_inq_varid(ncid, 'y', ncv) )
272  call check( nf90_get_var(ncid, ncv, eta_conv) )
273 
274  call check( nf90_inq_varid(ncid, 'sigma_level_c', ncv) )
275  call check( nf90_get_var(ncid, ncv, sigma_level_c_conv) )
276 
277  call check( nf90_inq_varid(ncid, 'sigma_level_t', ncv) )
278  call check( nf90_get_var(ncid, ncv, sigma_level_t_conv) )
279 
280  call check( nf90_inq_varid(ncid, 'sigma_level_r', ncv) )
281  call check( nf90_get_var(ncid, ncv, sigma_level_r_conv) )
282 
283  call check( nf90_inq_varid(ncid, 'lon', ncv) )
284  call check( nf90_get_var(ncid, ncv, lon_conv) )
285 
286  call check( nf90_inq_varid(ncid, 'lat', ncv) )
287  call check( nf90_get_var(ncid, ncv, lat_conv) )
288 
289  call check( nf90_inq_varid(ncid, 'lambda', ncv) )
290  call check( nf90_get_var(ncid, ncv, lambda_conv) )
291 
292  call check( nf90_inq_varid(ncid, 'phi', ncv) )
293  call check( nf90_get_var(ncid, ncv, phi_conv) )
294 
295  call check( nf90_inq_varid(ncid, 'temp_s', ncv) )
296  call check( nf90_get_var(ncid, ncv, temp_s_conv) )
297 
298  if ( nf90_inq_varid(ncid, 'prec', ncv) == nf90_noerr ) then
299  call check( nf90_get_var(ncid, ncv, accum_conv) )
300  else
301  write(6,'(/1x,a)') '>>> read_erg_nc: Variable prec'
302  write(6, '(1x,a)') ' not available in read file *.nc.'
303  accum_conv = 0.0_sp
304  end if
305 
306  call check( nf90_inq_varid(ncid, 'as_perp', ncv) )
307  call check( nf90_get_var(ncid, ncv, as_perp_conv) )
308 
309  if ( nf90_inq_varid(ncid, 'as_perp_apl', ncv) == nf90_noerr ) then
310  call check( nf90_get_var(ncid, ncv, as_perp_apl_conv) )
311  else
312  write(6,'(/1x,a)') '>>> read_erg_nc: Variable as_perp_apl'
313  write(6, '(1x,a)') ' not available in read file *.nc.'
314  as_perp_apl_conv = 0.0_sp
315  end if
316 
317  if ( nf90_inq_varid(ncid, 'smb_corr', ncv) == nf90_noerr ) then
318  call check( nf90_get_var(ncid, ncv, smb_corr_conv) )
319  else
320  write(6,'(/1x,a)') '>>> read_erg_nc: Variable smb_corr'
321  write(6, '(1x,a)') ' not available in read file *.nc.'
322  smb_corr_conv = 0.0_sp
323  end if
324 
325 #if (DISC>0) /* Ice discharge parameterisation */
326 
327  if ( nf90_inq_varid(ncid, 'dis_perp', ncv) == nf90_noerr ) then
328  call check( nf90_get_var(ncid, ncv, dis_perp_conv) )
329  else
330  write(6,'(/1x,a)') '>>> read_erg_nc: Variable dis_perp'
331  write(6, '(1x,a)') ' not available in read file *.nc.'
332  dis_perp_conv = 0.0_sp
333  end if
334 
335  if ( nf90_inq_varid(ncid, 'cst_dist', ncv) == nf90_noerr ) then
336  call check( nf90_get_var(ncid, ncv, cst_dist_conv) )
337  else
338  write(6,'(/1x,a)') '>>> read_erg_nc: Variable cst_dist'
339  write(6, '(1x,a)') ' not available in read file *.nc.'
340  cst_dist_conv = real(no_value_pos_1,sp)
341  end if
342 
343  if ( nf90_inq_varid(ncid, 'cos_grad_tc', ncv) == nf90_noerr ) then
344  call check( nf90_get_var(ncid, ncv, cos_grad_tc_conv) )
345  else
346  write(6,'(/1x,a)') '>>> read_erg_nc: Variable cos_grad_tc'
347  write(6, '(1x,a)') ' not available in read file *.nc.'
348  cos_grad_tc_conv = 1.0_sp
349  end if
350 
351  if ( nf90_inq_varid(ncid, 'mask_mar', ncv) == nf90_noerr ) then
352  call check( nf90_get_var(ncid, ncv, mask_mar_conv) )
353  else
354  write(6,'(/1x,a)') '>>> read_erg_nc: Variable mask_mar'
355  write(6, '(1x,a)') ' not available in read file *.nc.'
356  mask_mar_conv = 0_i2b
357  end if
358 
359 #endif
360 
361  if ( nf90_inq_varid(ncid, 'q_geo', ncv) == nf90_noerr ) then
362  call check( nf90_get_var(ncid, ncv, q_geo_conv) )
363  else
364  write(6,'(/1x,a)') '>>> read_erg_nc: Variable q_geo'
365  write(6, '(1x,a)') ' not available in read file *.nc.'
366  q_geo_conv = 0.0_sp
367  end if
368 
369  call check( nf90_inq_varid(ncid, 'maske', ncv) )
370  call check( nf90_get_var(ncid, ncv, maske_conv) )
371 
372  if ( nf90_inq_varid(ncid, 'maske_old', ncv) == nf90_noerr ) then
373  call check( nf90_get_var(ncid, ncv, maske_old_conv) )
374  else
375  write(6,'(/1x,a)') '>>> read_erg_nc: Variable maske_old'
376  write(6, '(1x,a)') ' not available in read file *.nc.'
377  maske_old_conv = maske_conv
378  end if
379 
380  call check( nf90_inq_varid(ncid, 'n_cts', ncv) )
381  call check( nf90_get_var(ncid, ncv, n_cts_conv) )
382 
383  call check( nf90_inq_varid(ncid, 'kc_cts', ncv) )
384  call check( nf90_get_var(ncid, ncv, kc_cts_conv) )
385 
386  call check( nf90_inq_varid(ncid, 'zs', ncv) )
387  call check( nf90_get_var(ncid, ncv, zs_conv) )
388 
389  call check( nf90_inq_varid(ncid, 'zm', ncv) )
390  call check( nf90_get_var(ncid, ncv, zm_conv) )
391 
392  call check( nf90_inq_varid(ncid, 'zb', ncv) )
393  call check( nf90_get_var(ncid, ncv, zb_conv) )
394 
395  call check( nf90_inq_varid(ncid, 'zl', ncv) )
396  call check( nf90_get_var(ncid, ncv, zl_conv) )
397 
398  if ( nf90_inq_varid(ncid, 'zl0', ncv) == nf90_noerr ) then
399  call check( nf90_get_var(ncid, ncv, zl0_conv) )
400  else
401  write(6,'(/1x,a)') '>>> read_erg_nc: Variable zl0'
402  write(6, '(1x,a)') ' not available in read file *.nc.'
403  zl0_conv = 0.0_sp
404  end if
405 
406  call check( nf90_inq_varid(ncid, 'H_cold', ncv) )
407  call check( nf90_get_var(ncid, ncv, h_cold_conv) )
408 
409  call check( nf90_inq_varid(ncid, 'H_temp', ncv) )
410  call check( nf90_get_var(ncid, ncv, h_temp_conv) )
411 
412  call check( nf90_inq_varid(ncid, 'H', ncv) )
413  call check( nf90_get_var(ncid, ncv, h_conv) )
414 
415  call check( nf90_inq_varid(ncid, 'H_R', ncv) )
416  call check( nf90_get_var(ncid, ncv, h_r_conv) )
417 
418  call check( nf90_inq_varid(ncid, 'Q_bm', ncv) )
419  call check( nf90_get_var(ncid, ncv, q_bm_conv) )
420 
421  call check( nf90_inq_varid(ncid, 'Q_tld', ncv) )
422  call check( nf90_get_var(ncid, ncv, q_tld_conv) )
423 
424  call check( nf90_inq_varid(ncid, 'am_perp', ncv) )
425  call check( nf90_get_var(ncid, ncv, am_perp_conv) )
426 
427  call check( nf90_inq_varid(ncid, 'qx', ncv) )
428  call check( nf90_get_var(ncid, ncv, qx_conv) )
429 
430  call check( nf90_inq_varid(ncid, 'qy', ncv) )
431  call check( nf90_get_var(ncid, ncv, qy_conv) )
432 
433  call check( nf90_inq_varid(ncid, 'dzs_dt', ncv) )
434  call check( nf90_get_var(ncid, ncv, dzs_dtau_conv) )
435 
436  call check( nf90_inq_varid(ncid, 'dzm_dt', ncv) )
437  call check( nf90_get_var(ncid, ncv, dzm_dtau_conv) )
438 
439  call check( nf90_inq_varid(ncid, 'dzb_dt', ncv) )
440  call check( nf90_get_var(ncid, ncv, dzb_dtau_conv) )
441 
442  call check( nf90_inq_varid(ncid, 'dzl_dt', ncv) )
443  call check( nf90_get_var(ncid, ncv, dzl_dtau_conv) )
444 
445  call check( nf90_inq_varid(ncid, 'dH_c_dt', ncv) )
446  call check( nf90_get_var(ncid, ncv, dh_c_dtau_conv) )
447 
448  call check( nf90_inq_varid(ncid, 'dH_t_dt', ncv) )
449  call check( nf90_get_var(ncid, ncv, dh_t_dtau_conv) )
450 
451  call check( nf90_inq_varid(ncid, 'dH_dt', ncv) )
452  call check( nf90_get_var(ncid, ncv, dh_dtau_conv) )
453 
454  call check( nf90_inq_varid(ncid, 'vx_b_g', ncv) )
455  call check( nf90_get_var(ncid, ncv, vx_b_g_conv) )
456 
457  call check( nf90_inq_varid(ncid, 'vy_b_g', ncv) )
458  call check( nf90_get_var(ncid, ncv, vy_b_g_conv) )
459 
460  call check( nf90_inq_varid(ncid, 'vz_b', ncv) )
461  call check( nf90_get_var(ncid, ncv, vz_b_conv) )
462 
463  call check( nf90_inq_varid(ncid, 'vh_b', ncv) )
464  call check( nf90_get_var(ncid, ncv, vh_b_conv) )
465 
466  call check( nf90_inq_varid(ncid, 'vx_s_g', ncv) )
467  call check( nf90_get_var(ncid, ncv, vx_s_g_conv) )
468 
469  call check( nf90_inq_varid(ncid, 'vy_s_g', ncv) )
470  call check( nf90_get_var(ncid, ncv, vy_s_g_conv) )
471 
472  call check( nf90_inq_varid(ncid, 'vz_s', ncv) )
473  call check( nf90_get_var(ncid, ncv, vz_s_conv) )
474 
475  call check( nf90_inq_varid(ncid, 'vh_s', ncv) )
476  call check( nf90_get_var(ncid, ncv, vh_s_conv) )
477 
478  if ( nf90_inq_varid(ncid, 'vx_m_g', ncv) == nf90_noerr ) then
479  call check( nf90_get_var(ncid, ncv, vx_m_g_conv) )
480  else
481  write(6,'(/1x,a)') '>>> read_erg_nc: Variable vx_m_g'
482  write(6, '(1x,a)') ' not available in read file *.nc.'
483  vx_m_g_conv = 0.0_sp
484  end if
485 
486  if ( nf90_inq_varid(ncid, 'vy_m_g', ncv) == nf90_noerr ) then
487  call check( nf90_get_var(ncid, ncv, vy_m_g_conv) )
488  else
489  write(6,'(/1x,a)') '>>> read_erg_nc: Variable vy_m_g'
490  write(6, '(1x,a)') ' not available in read file *.nc.'
491  vy_m_g_conv = 0.0_sp
492  end if
493 
494  if ( nf90_inq_varid(ncid, 'vh_m', ncv) == nf90_noerr ) then
495  call check( nf90_get_var(ncid, ncv, vh_m_conv) )
496  else
497  write(6,'(/1x,a)') '>>> read_erg_nc: Variable vh_m'
498  write(6, '(1x,a)') ' not available in read file *.nc.'
499  vh_m_conv = 0.0_sp
500  end if
501 
502  call check( nf90_inq_varid(ncid, 'temp_b', ncv) )
503  call check( nf90_get_var(ncid, ncv, temp_b_conv) )
504 
505  call check( nf90_inq_varid(ncid, 'temph_b', ncv) )
506  call check( nf90_get_var(ncid, ncv, temph_b_conv) )
507 
508  if ( nf90_inq_varid(ncid, 'tau_b_driving', ncv) == nf90_noerr ) then
509  call check( nf90_get_var(ncid, ncv, tau_b_driving_conv) )
510  else
511  write(6,'(/1x,a)') '>>> read_erg_nc: Variable tau_b_driving'
512  write(6, '(1x,a)') ' not available in read file *.nc.'
513  tau_b_driving_conv = 0.0_sp
514  end if
515 
516  if ( nf90_inq_varid(ncid, 'tau_b_drag', ncv) == nf90_noerr ) then
517  call check( nf90_get_var(ncid, ncv, tau_b_drag_conv) )
518  else
519  write(6,'(/1x,a)') '>>> read_erg_nc: Variable tau_b_drag'
520  write(6, '(1x,a)') ' not available in read file *.nc.'
521  tau_b_drag_conv = 0.0_sp
522  end if
523 
524  call check( nf90_inq_varid(ncid, 'p_b_w', ncv) )
525  call check( nf90_get_var(ncid, ncv, p_b_w_conv) )
526 
527  call check( nf90_inq_varid(ncid, 'H_w', ncv) )
528  call check( nf90_get_var(ncid, ncv, h_w_conv) )
529 
530  call check( nf90_inq_varid(ncid, 'q_gl_g', ncv) )
531  call check( nf90_get_var(ncid, ncv, q_gl_g_conv) )
532 
533  if ( nf90_inq_varid(ncid, 'q_cf_g', ncv) == nf90_noerr ) then
534  call check( nf90_get_var(ncid, ncv, q_cf_g_conv) )
535  else
536  write(6,'(/1x,a)') '>>> read_erg_nc: Variable q_cf_g'
537  write(6, '(1x,a)') ' not available in read file *.nc.'
538  q_cf_g_conv = 0.0_sp
539  end if
540 
541  call check( nf90_inq_varid(ncid, 'vx_c', ncv) )
542  call check( nf90_get_var(ncid, ncv, vx_c_conv) )
543 
544  call check( nf90_inq_varid(ncid, 'vy_c', ncv) )
545  call check( nf90_get_var(ncid, ncv, vy_c_conv) )
546 
547  call check( nf90_inq_varid(ncid, 'vz_c', ncv) )
548  call check( nf90_get_var(ncid, ncv, vz_c_conv) )
549 
550  call check( nf90_inq_varid(ncid, 'vx_t', ncv) )
551  call check( nf90_get_var(ncid, ncv, vx_t_conv) )
552 
553  call check( nf90_inq_varid(ncid, 'vy_t', ncv) )
554  call check( nf90_get_var(ncid, ncv, vy_t_conv) )
555 
556  call check( nf90_inq_varid(ncid, 'vz_t', ncv) )
557  call check( nf90_get_var(ncid, ncv, vz_t_conv) )
558 
559  call check( nf90_inq_varid(ncid, 'temp_c', ncv) )
560  call check( nf90_get_var(ncid, ncv, temp_c_conv) )
561 
562  call check( nf90_inq_varid(ncid, 'omega_t', ncv) )
563  call check( nf90_get_var(ncid, ncv, omega_t_conv) )
564 
565  call check( nf90_inq_varid(ncid, 'temp_r', ncv) )
566  call check( nf90_get_var(ncid, ncv, temp_r_conv) )
567 
568  call check( nf90_inq_varid(ncid, 'enth_c', ncv) )
569  call check( nf90_get_var(ncid, ncv, enth_c_conv) )
570 
571  call check( nf90_inq_varid(ncid, 'enth_t', ncv) )
572  call check( nf90_get_var(ncid, ncv, enth_t_conv) )
573 
574  call check( nf90_inq_varid(ncid, 'omega_c', ncv) )
575  call check( nf90_get_var(ncid, ncv, omega_c_conv) )
576 
577  if ( nf90_inq_varid(ncid, 'enh_c', ncv) == nf90_noerr ) then
578  call check( nf90_get_var(ncid, ncv, enh_c_conv) )
579  else
580  write(6,'(/1x,a)') '>>> read_erg_nc: Variable enh_c'
581  write(6, '(1x,a)') ' not available in read file *.nc.'
582  enh_c_conv = 1.0_sp
583  end if
584 
585  if ( nf90_inq_varid(ncid, 'enh_t', ncv) == nf90_noerr ) then
586  call check( nf90_get_var(ncid, ncv, enh_t_conv) )
587  else
588  write(6,'(/1x,a)') '>>> read_erg_nc: Variable enh_t'
589  write(6, '(1x,a)') ' not available in read file *.nc.'
590  enh_t_conv = 1.0_sp
591  end if
592 
593  call check( nf90_inq_varid(ncid, 'age_c', ncv) )
594  call check( nf90_get_var(ncid, ncv, age_c_conv) )
595 
596  call check( nf90_inq_varid(ncid, 'age_t', ncv) )
597  call check( nf90_get_var(ncid, ncv, age_t_conv) )
598 
599  call check( nf90_close(ncid) )
600 
601 #endif
602 
603 !-------- Convert data to real*8 and years to seconds --------
604 
605  z_sl = real(z_sl_conv,dp)
606 
607  h_r = real(h_r_conv,dp)
608 
609  do i=0, imax
610  xi(i) = real(xi_conv(i),dp)
611  end do
612 
613  do j=0, jmax
614  eta(j) = real(eta_conv(j),dp)
615  end do
616 
617  do i=0, imax
618  do j=0, jmax
619 
620  maske(j,i) = maske_conv(i,j)
621  maske_old(j,i) = maske_old_conv(i,j)
622  n_cts(j,i) = n_cts_conv(i,j)
623  kc_cts(j,i) = kc_cts_conv(i,j)
624  zs(j,i) = real(zs_conv(i,j),dp)
625  zm(j,i) = real(zm_conv(i,j),dp)
626  zb(j,i) = real(zb_conv(i,j),dp)
627  zl(j,i) = real(zl_conv(i,j),dp)
628 #if (CALCMOD==1)
629  h_c(j,i) = real(H_cold_conv(i,j),dp)
630  h_t(j,i) = real(H_temp_conv(i,j),dp)
631 #elif (CALCMOD==0 || CALCMOD==2 || CALCMOD==3 || CALCMOD==-1)
632  h_c(j,i) = real(H_conv(i,j),dp)
633  h_t(j,i) = 0.0_dp
634 #endif
635  q_bm(j,i) = real(Q_bm_conv(i,j),dp)/year_sec
636  q_tld(j,i) = real(Q_tld_conv(i,j),dp)/year_sec
637  am_perp(j,i) = real(am_perp_conv(i,j),dp)/year_sec
638  qx(j,i) = real(qx_conv(i,j),dp)/year_sec
639  qy(j,i) = real(qy_conv(i,j),dp)/year_sec
640  dzs_dtau(j,i) = real(dzs_dtau_conv(i,j),dp)/year_sec
641  dzm_dtau(j,i) = real(dzm_dtau_conv(i,j),dp)/year_sec
642  dzb_dtau(j,i) = real(dzb_dtau_conv(i,j),dp)/year_sec
643  dzl_dtau(j,i) = real(dzl_dtau_conv(i,j),dp)/year_sec
644  dh_c_dtau(j,i) = real(dH_c_dtau_conv(i,j),dp)/year_sec
645  dh_t_dtau(j,i) = real(dH_t_dtau_conv(i,j),dp)/year_sec
646  vx_b_g(j,i) = real(vx_b_g_conv(i,j),dp)/year_sec
647  vy_b_g(j,i) = real(vy_b_g_conv(i,j),dp)/year_sec
648  vz_b(j,i) = real(vz_b_conv(i,j),dp)/year_sec
649  vx_s_g(j,i) = real(vx_s_g_conv(i,j),dp)/year_sec
650  vy_s_g(j,i) = real(vy_s_g_conv(i,j),dp)/year_sec
651  vz_s(j,i) = real(vz_s_conv(i,j),dp)/year_sec
652  temp_b(j,i) = real(temp_b_conv(i,j),dp)
653  temph_b(j,i) = real(temph_b_conv(i,j),dp)
654  p_b_w(j,i) = real(p_b_w_conv(i,j),dp)
655  h_w(j,i) = real(H_w_conv(i,j),dp)
656 
657  do kr=0, krmax
658  temp_r(kr,j,i) = real(temp_r_conv(i,j,kr),dp)
659  end do
660 
661  do kt=0, ktmax
662  vx_t(kt,j,i) = real(vx_t_conv(i,j,kt),dp)/year_sec
663  vy_t(kt,j,i) = real(vy_t_conv(i,j,kt),dp)/year_sec
664  vz_t(kt,j,i) = real(vz_t_conv(i,j,kt),dp)/year_sec
665  omega_t(kt,j,i) = real(omega_t_conv(i,j,kt),dp)
666  age_t(kt,j,i) = real(age_t_conv(i,j,kt),dp)*year_sec
667  enth_t(kt,j,i) = real(enth_t_conv(i,j,kt),dp)
668  enh_t(kt,j,i) = real(enh_t_conv(i,j,kt),dp)
669  end do
670 
671  do kc=0, kcmax
672  vx_c(kc,j,i) = real(vx_c_conv(i,j,kc),dp)/year_sec
673  vy_c(kc,j,i) = real(vy_c_conv(i,j,kc),dp)/year_sec
674  vz_c(kc,j,i) = real(vz_c_conv(i,j,kc),dp)/year_sec
675  temp_c(kc,j,i) = real(temp_c_conv(i,j,kc),dp)
676  age_c(kc,j,i) = real(age_c_conv(i,j,kc),dp)*year_sec
677  enth_c(kc,j,i) = real(enth_c_conv(i,j,kc),dp)
678  omega_c(kc,j,i) = real(omega_c_conv(i,j,kc),dp)
679  enh_c(kc,j,i) = real(enh_c_conv(i,j,kc),dp)
680  end do
681 
682  end do
683  end do
684 
685  end subroutine read_erg_nc
686 
687 !-------------------------------------------------------------------------------
688 !> Reading of the target-topography file (in NetCDF format).
689 !<------------------------------------------------------------------------------
690  subroutine read_target_topo_nc(target_topo_dat_name)
692  use sico_variables_m, only : maske_target, &
695 
696 #if (NETCDF==2)
697  use netcdf
698  use nc_check_m
699 #endif
700 
701  implicit none
702 
703  character(len=100), intent(in) :: target_topo_dat_name
704 
705 ! Return variables
706 ! (defined as global variables in module sico_variables_m):
707 !
708 ! maske_target, zs_target, zb_target, zl_target, H_target
709 
710  integer(i2b), dimension(0:IMAX,0:JMAX) :: maske_conv
711  real(sp), dimension(0:IMAX,0:JMAX) :: zs_conv, zb_conv, zl_conv, H_conv
712  character(len=256) :: target_topo_dat_path
713  character(len=256) :: filename_with_path
714 
715 #if (NETCDF==2)
716  integer(i4b) :: ncid, ncv
717  ! ncid: ID of the output file
718  ! ncv: Variable ID
719 #endif
720 
721  character(len=64), parameter :: thisroutine = 'read_target_topo_nc'
722 
723 !-------- Read target-topography data
724 ! (from a time-slice file of a previous simulation) --------
725 
726 #if (NETCDF==2)
727 
728 #if (defined(TARGET_TOPO_DAT_PATH))
729  target_topo_dat_path = trim(target_topo_dat_path)
730 #else
731  target_topo_dat_path = 'dummy'
732  errormsg = ' >>> read_target_topo_nc: TARGET_TOPO_DAT_PATH must be defined!'
733  call error(errormsg)
734 #endif
735 
736  filename_with_path &
737  = trim(target_topo_dat_path)//'/'//trim(target_topo_dat_name)
738 
739  call check( nf90_open(trim(filename_with_path), nf90_nowrite, ncid), &
740  thisroutine )
741 
742  call check( nf90_inq_varid(ncid, 'maske', ncv), thisroutine )
743  call check( nf90_get_var(ncid, ncv, maske_conv), thisroutine )
744 
745  call check( nf90_inq_varid(ncid, 'zs', ncv), thisroutine )
746  call check( nf90_get_var(ncid, ncv, zs_conv), thisroutine )
747 
748  call check( nf90_inq_varid(ncid, 'zb', ncv), thisroutine )
749  call check( nf90_get_var(ncid, ncv, zb_conv), thisroutine )
750 
751  call check( nf90_inq_varid(ncid, 'zl', ncv), thisroutine )
752  call check( nf90_get_var(ncid, ncv, zl_conv), thisroutine )
753 
754  call check( nf90_inq_varid(ncid, 'H', ncv), thisroutine )
755  call check( nf90_get_var(ncid, ncv, h_conv), thisroutine )
756 
757  call check( nf90_close(ncid), thisroutine )
758 
759 #else /* NETCDF != 2, not allowed */
760 
761  maske_conv = 0_i2b ! dummy values
762  zs_conv = 0.0_sp ! dummy values
763  zb_conv = 0.0_sp ! dummy values
764  zl_conv = 0.0_sp ! dummy values
765  h_conv = 0.0_sp ! dummy values
766 
767  errormsg = ' >>> read_target_topo_nc: Parameter NETCDF must be set to 2!'
768  call error(errormsg)
769 
770 #endif
771 
772 !-------- Convert data to double precision --------
773 
774 #ifndef ALLOW_OPENAD /* Normal */
775  maske_target = transpose(maske_conv)
776  zs_target = real(transpose(zs_conv),dp)
777  zb_target = real(transpose(zb_conv),dp)
778  zl_target = real(transpose(zl_conv),dp)
779  h_target = real(transpose(H_conv) ,dp)
780 #else /* OpenAD */
781  do i=0, imax
782  do j=0, jmax
783  maske_target(j,i) = maske_conv(i,j)
784  zs_target(j,i) = real(zs_conv(i,j),dp)
785  zb_target(j,i) = real(zb_conv(i,j),dp)
786  zl_target(j,i) = real(zl_conv(i,j),dp)
787  h_target(j,i) = real(H_conv(i,j) ,dp)
788  end do
789  end do
790 #endif /* Normal vs. OpenAD */
791 
792  end subroutine read_target_topo_nc
793 
794 !-------------------------------------------------------------------------------
795 !> Reading of the tabulated kei function.
796 !<------------------------------------------------------------------------------
797  subroutine read_kei()
800 
801  implicit none
802 
803  integer(i4b) :: n, n_data_kei_max
804  integer(i4b) :: ios
805  real(dp) :: r_val, d_dummy
806  character :: ch_dummy
807  character(len=256) :: filename_with_path
808 
809  n_data_kei_max = 10000
810 
811  kei = 0.0_dp
812 
813 !-------- Reading of tabulated values --------
814 
815  filename_with_path = trim(inpath)//'/general/kei.dat'
816 
817 #ifndef ALLOW_OPENAD /* Normal */
818  open(unit=11, iostat=ios, file=trim(filename_with_path), status='old')
819 #else /* OpenAD */
820  open(unit=11, iostat=ios, file=trim(filename_with_path))
821 #endif /* Normal vs. OpenAD */
822 
823  if (ios /= 0) then
824  errormsg = ' >>> read_kei: Error when opening the kei file!'
825  call error(errormsg)
826  end if
827 
828  read(unit=11, fmt='(a)') ch_dummy
829  read(unit=11, fmt='(15x,f5.0)') kei_r_max
830  read(unit=11, fmt='(15x,f5.0)') kei_r_incr
831  read(unit=11, fmt='(a)') ch_dummy
832 
834 
835  if (n_data_kei > n_data_kei_max) then
836  errormsg = ' >>> read_kei: Array kei too small!'
837  call error(errormsg)
838  end if
839 
840  read(unit=11, fmt='(a)') ch_dummy
841  read(unit=11, fmt='(a)') ch_dummy
842  read(unit=11, fmt='(a)') ch_dummy
843 
844  do n=-n_data_kei, n_data_kei
845  read(unit=11, fmt=*) d_dummy, kei(n)
846  end do
847 
848  close(unit=11, status='keep')
849 
850  end subroutine read_kei
851 
852 !-------------------------------------------------------------------------------
853 !> Reading of physical parameters.
854 !<------------------------------------------------------------------------------
855  subroutine read_phys_para()
857  use sico_variables_m
858  use sico_vars_m
859 
860  implicit none
861 
862  integer(i4b), parameter :: n_unit=31
863  integer(i4b) :: ios
864  integer(i4b) :: n
865  character(len=256) :: filename_with_path
866 
867  filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
868  trim(phys_para_file)
869 
870 #ifndef ALLOW_OPENAD /* Normal */
871  open(n_unit, iostat=ios, file=trim(filename_with_path), status='old')
872 #else /* OpenAD */
873  open(n_unit, iostat=ios, file=trim(filename_with_path))
874 #endif /* Normal vs. OpenAD */
875 
876  if (ios /= 0) then
877  errormsg = ' >>> read_phys_para: Error when opening the phys_para file!'
878  call error(errormsg)
879  end if
880 
881 #if (!defined(NMARS) && !defined(SMARS)) /* not Martian polar caps */
882  call read_phys_para_value(rho, n_unit)
883 #else /* Martian polar caps */
884  call read_phys_para_value(rho_i, n_unit)
885 #endif
886  call read_phys_para_value(rho_w, n_unit)
887  call read_phys_para_value(rho_sw, n_unit)
888  call read_phys_para_value(l, n_unit)
889  call read_phys_para_value(g, n_unit)
890  call read_phys_para_value(nue, n_unit)
891  call read_phys_para_value(beta, n_unit)
892 #if (!defined(NMARS) && !defined(SMARS)) /* not Martian polar caps */
893  call read_phys_para_value(delta_tm_sw, n_unit)
894 #endif
895  call read_phys_para_value(omega_max, n_unit)
896 #if (defined(NMARS) || defined(SMARS)) /* Martian polar caps */
897  call read_phys_para_value(rho_c, n_unit)
898  call read_phys_para_value(kappa_c, n_unit)
899  call read_phys_para_value(c_c, n_unit)
900 #endif
901  call read_phys_para_value(h_r, n_unit)
902  call read_phys_para_value(rho_c_r, n_unit)
903  call read_phys_para_value(kappa_r, n_unit)
904  call read_phys_para_value(rho_a, n_unit)
905  call read_phys_para_value(r_t, n_unit)
906 
907  call read_phys_para_value(r, n_unit)
908  call read_phys_para_value(a, n_unit)
909  call read_phys_para_value(b, n_unit)
910  call read_phys_para_value(phi0, n_unit)
911  call read_phys_para_value(lambda0, n_unit)
912 
913 #if (!defined(NMARS) && !defined(SMARS)) /* not Martian polar caps */
914  call read_phys_para_value(s_stat_0, n_unit)
915 #if (!defined(GRL)) /* not Greenland */
916  call read_phys_para_value(beta1_0, n_unit)
917  call read_phys_para_value(beta2_0, n_unit)
918 #else /* Greenland */
919  call read_phys_para_value(beta1_lt_0, n_unit)
920  call read_phys_para_value(beta1_ht_0, n_unit)
921  call read_phys_para_value(beta2_lt_0, n_unit)
922  call read_phys_para_value(beta2_ht_0, n_unit)
923  call read_phys_para_value(phi_sep_0, n_unit)
924 #endif
925  call read_phys_para_value(pmax_0, n_unit)
926  call read_phys_para_value(mu_0, n_unit)
927 #endif
928 
929  do n=10, -190, -1
930  call read_phys_para_value(rf(n), n_unit)
931  end do
932 
933  do n=10, -190, -1
934  call read_phys_para_value(kappa(n), n_unit)
935  end do
936 
937  do n=10, -190, -1
938  call read_phys_para_value(c(n), n_unit)
939  end do
940 
941  close(n_unit, status='keep')
942 
943  end subroutine read_phys_para
944 
945 !-------------------------------------------------------------------------------
946 !> Reading of a value of a physical parameter from the phys_para file.
947 !<------------------------------------------------------------------------------
948  subroutine read_phys_para_value(para, n_unit)
950  implicit none
951 
952  integer(i4b), intent(in) :: n_unit
953  real(dp), intent(out) :: para
954 
955  character :: check
956 
957 #ifndef ALLOW_OPENAD /* Normal */
958 
959  integer :: i0
960  character(len=1024) :: txt
961 
962 #else /* OpenAD */
963 
964  character :: ch_dummy
965  character(len=*), parameter :: fmt1='(a)', fmt3='(15x)'
966  logical :: first_read
967 
968  first_read = .true.
969 
970 #endif /* Normal vs. OpenAD */
971 
972  check = '%'
973  do while (check == '%')
974 
975 #ifndef ALLOW_OPENAD /* Normal */
976 
977  read(unit=n_unit, fmt='(a)') txt
978 
979  check = txt(1:1)
980 
981  if (check /= '%') then ! no comment line
982  i0 = index(txt, '=') + 1
983  read(txt(i0:), *) para
984  end if
985 
986 #else /* OpenAD */
987 
988  if (first_read) then
989  read(unit=n_unit, fmt=trim(fmt1), advance='no') check
990  first_read=.false.
991  else
992  read(unit=n_unit, fmt=trim(fmt1)) ch_dummy
993  read(unit=n_unit, fmt=trim(fmt1), advance='no') check
994  end if
995 
996  if (check /= '%') then ! no comment line
997  read(unit=n_unit, fmt=trim(fmt3), advance='no')
998  read(unit=n_unit, fmt=*) para
999  end if
1000 
1001 #endif /* Normal vs. OpenAD */
1002 
1003  end do
1004 
1005  end subroutine read_phys_para_value
1006 
1007 !-------------------------------------------------------------------------------
1008 
1009 end module read_m
1010 !
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) vx_c
vx_c(kc,j,i): Velocity in x-direction in the upper (kc) ice domain (at (i+1/2,j,kc)) ...
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.
subroutine, public read_erg_nc(z_sl, filename)
Reading of time-slice files in native binary or NetCDF format.
Definition: read_m.F90:51
real(dp) phi_sep_0
PHI_SEP_0: Separation latitude for the computation of the degree-day factors beta1 and beta2: Equator...
subroutine, public read_kei()
Reading of the tabulated kei function.
Definition: read_m.F90:798
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) age_c
age_c(kc,j,i): Age in the upper (kc) ice domain
integer(i2b), 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) beta2_ht_0
BETA2_HT_0: Degree-day factor for ice at high summer temperatures.
real(dp) beta1_ht_0
BETA1_HT_0: Degree-day factor for snow at high summer temperatures.
real(dp), dimension(0:jmax, 0:imax) vx_b_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
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:kcmax, 0:jmax, 0:imax) enh_c
enh_c(kc,j,i): Flow enhancement factor in the upper (kc) ice domain
real(dp) rho_w
RHO_W: Density of pure water.
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.
real(dp) a
A: Semi-major axis of the planet.
subroutine error(error_message)
Main routine of error_m: Writing of error messages and stopping execution.
Definition: error_m.F90:47
real(dp), dimension(0:jmax, 0:imax) vz_b
vz_b(j,i): Velocity in z-direction at the ice base, at (i,j)
real(dp) mu_0
MU_0: Firn-warming correction.
subroutine, public read_phys_para()
Reading of physical parameters.
Definition: read_m.F90:856
real(dp), dimension(0:jmax, 0:imax) vx_s_g
vx_s_g(j,i): Velocity in x-direction at the ice surface, at (i,j)
real(dp) omega_max
OMEGA_MAX: Threshold value for the water content of temperate ice.
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) enth_t
enth_t(kt,j,i): Enthalpy in the lower (kt) ice domain
real(dp), dimension(0:jmax, 0:imax) vz_s
vz_s(j,i): Velocity in z-direction at the ice surface, at (i,j)
real(dp), dimension(0:jmax, 0:imax) dh_t_dtau
dH_t_dtau(j,i): Derivative of H_t by tau (time)
real(dp), dimension(0:jmax, 0:imax) dh_c_dtau
dH_c_dtau(j,i): Derivative of H_c by tau (time)
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 ...
real(dp), dimension(0:jmax, 0:imax) p_b_w
p_b_w(j,i): Basal water pressure
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) zl_target
zl_target(j,i): Target topography (lithosphere surface)
real(dp) beta2_lt_0
BETA2_LT_0: Degree-day factor for ice at low summer temperatures.
real(dp), dimension(0:jmax, 0:imax) zs_target
zs_target(j,i): Target topography (ice surface)
real(dp) b
B: Semi-minor axis of the planet.
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) vz_t
vz_t(kt,j,i): Velocity in z-direction in the lower (kt) ice domain (at (i,j,kt+1/2)) ...
real(dp), dimension(0:jmax, 0:imax) zl
zl(j,i): Coordinate z of the lithosphere surface
real(dp) g
G: Acceleration due to gravity.
real(dp) r
R: Radius of the planet.
integer(i2b), 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) q_tld
Q_tld(j,i): Water drainage rate from the temperate layer.
Declarations of global variables for SICOPOLIS (for the ANT domain).
Definition: sico_vars_m.F90:35
subroutine read_phys_para_value(para, n_unit)
Reading of a value of a physical parameter from the phys_para file.
Definition: read_m.F90:949
real(dp), dimension(0:jmax, 0:imax) temp_b
temp_b(j,i): Basal temperature
real(dp) beta1_0
BETA1_0: Degree-day factor for snow.
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) enth_c
enth_c(kc,j,i): Enthalpy in the upper (kc) ice domain
real(dp) rho_c_r
RHO_C_R: Density times specific heat of the lithosphere.
Reading of several input data.
Definition: read_m.F90:35
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:ktmax, 0:jmax, 0:imax) enh_t
enh_t(kt,j,i): Flow enhancement factor in the lower (kt) ice domain
real(dp) rho_c
RHO_C: Density of crustal material (dust)
Definition: sico_vars_m.F90:79
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:jmax, 0:imax) zb_target
zb_target(j,i): Target topography (ice base)
real(dp) kei_r_max
kei_r_max: Maximum value of the argument r of the tabulated kei function
real(dp), dimension(0:jmax, 0:imax) vy_b_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
real(dp) beta2_0
BETA2_0: Degree-day factor for ice.
real(dp), dimension(-190:10) kappa
KAPPA(n): Tabulated values for the heat conductivity of ice.
Writing of error messages and stopping execution.
Definition: error_m.F90:35
real(dp), dimension(0:jmax, 0:imax) dzm_dtau
dzm_dtau(j,i): Derivative of zm by tau (time)
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) omega_c
omega_c(kc,j,i): Water content in the upper (kc) ice domain
integer(i2b), dimension(0:jmax, 0:imax) maske_old
maske_old(j,i): Old value of maske (at the previous time step)
real(dp), dimension(0:jmax, 0:imax) qx
qx(j,i): Volume flux in x-direction (depth-integrated vx, at (i+1/2,j))
real(dp) kappa_c
KAPPA_C: Heat conductivity of crustal material (dust)
Definition: sico_vars_m.F90:81
real(dp) rho_i
RHO_I: Density of ice.
Definition: sico_vars_m.F90:77
real(dp) nue
NUE: Water diffusivity in ice.
real(dp) c_c
C_C: Specific heat of crustal material (dust)
Definition: sico_vars_m.F90:83
real(dp) beta1_lt_0
BETA1_LT_0: Degree-day factor for snow at low summer temperatures.
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) vz_c
vz_c(kc,j,i): Velocity in z-direction in the upper (kc) ice domain (at (i,j,kc+1/2)) ...
real(dp) phi0
PHI0: Standard parallel of the stereographic projection.
real(dp), dimension(0:jmax, 0:imax) vy_s_g
vy_s_g(j,i): Velocity in x-direction at the ice surface, at (i,j)
NetCDF error capturing.
Definition: nc_check_m.F90:35
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 (...
real(dp), dimension(-190:10) rf
RF(n): Tabulated values for the rate factor of cold ice.
real(dp), dimension(0:jmax, 0:imax) q_bm
Q_bm(j,i): Basal melting rate.
real(dp), dimension(-190:10) c
C(n): Tabulated values for the specific heat of ice.
real(dp), dimension(-10000:10000) kei
kei(n): Tabulated values of the kei function (Kelvin function of zero order)
subroutine check(status, ch_calling_routine)
NetCDF error capturing.
Definition: nc_check_m.F90:48
real(dp), dimension(0:jmax, 0:imax) temph_b
temph_b(j,i): Basal temperature relative to the pressure melting point
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) vy_t
vy_t(kt,j,i): Velocity in y-direction in the lower (kt) ice domain (at (i,j+1/2,kt)) ...
real(dp), dimension(0:jmax, 0:imax) am_perp
am_perp(j,i): Ice volume flux across the z=zm interface
real(dp) beta
BETA: Clausius-Clapeyron gradient of ice.
real(dp), parameter no_value_pos_1
no_value_pos_1: Positive no-value parameter
real(dp), dimension(0:jmax) eta
eta(j): Coordinate eta (= y) of grid point j
real(dp), dimension(0:imax) xi
xi(i): Coordinate xi (= x) of grid point i
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) vy_c
vy_c(kc,j,i): Velocity in y-direction in the upper (kc) ice domain (at (i,j+1/2,kc)) ...
character(len=16) ch_domain_short
ch_domain_short: Short name of the computational domain
real(dp) r_t
R_T: Coefficient of the water-content dependence in the rate factor for temperate ice...
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) omega_t
omega_t(kt,j,i): Water content in the lower (kt) ice domain
real(dp) kei_r_incr
kei_r_incr: Increment of the argument r of the tabulated kei function
real(dp) rho_sw
RHO_SW: Density of sea water.
real(dp) l
L: Latent heat of ice.
integer, parameter sp
Single-precision reals.
integer(i4b) n_data_kei
n_data_kei: Number of tabulated values of the kei function
character(len=256) errormsg
errormsg: Error-message string
real(dp) rho
RHO: Density of ice.
real(dp) rho_a
RHO_A: Density of the asthenosphere.
real(dp), dimension(0:jmax, 0:imax) h_target
H_target(j,i): Target topography (ice thickness)
integer, parameter dp
Double-precision reals.
real(dp), dimension(0:jmax, 0:imax) zb
zb(j,i): Coordinate z of the ice base
real(dp), dimension(0:jmax, 0:imax) dzb_dtau
dzb_dtau(j,i): Derivative of zb by tau (time)
integer(i2b), dimension(0:jmax, 0:imax) maske_target
maske_target(j,i): Target topography (ice-land-ocean mask)
subroutine, public read_target_topo_nc(target_topo_dat_name)
Reading of the target-topography file (in NetCDF format).
Definition: read_m.F90:691
real(dp) lambda0
LAMBDA0: Reference longitude (central meridian) of the stereographic projection.
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...
real(dp), dimension(0:jmax, 0:imax) h_w
H_w(j,i): Thickness of the water column under the ice base.
Declarations of global variables for SICOPOLIS.
real(dp), dimension(0:jmax, 0:imax) dzs_dtau
dzs_dtau(j,i): Derivative of zs by tau (time)
real(dp), dimension(0:jmax, 0:imax) dzl_dtau
dzl_dtau(j,i): Derivative of zl by tau (time)
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) vx_t
vx_t(kt,j,i): Velocity in x-direction in the lower (kt) ice domain (at (i+1/2,j,kt)) ...
real(dp), dimension(0:jmax, 0:imax) qy
qy(j,i): Volume flux in y-direction (depth-integrated vy, at (i,j+1/2))