SICOPOLIS V5-dev  Revision 1420
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-2019 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 #if (defined(ALLOW_GRDCHK) || defined(ALLOW_OPENAD))
46  public :: read_ad_data
47 #endif
48 
49 contains
50 
51 !-------------------------------------------------------------------------------
52 !> Reading of time-slice files in native binary or NetCDF format.
53 !<------------------------------------------------------------------------------
54  subroutine read_erg_nc(z_sl, filename, &
55  opt_maske, opt_n_cts, opt_kc_cts, &
56  opt_h_cold, opt_h_temp, opt_h, &
57  opt_temp_r, opt_omega_t, opt_age_t, &
58  opt_temp_c, opt_age_c, opt_omega_c, &
59  opt_flag_temp_age_only)
60 
62  use sico_vars_m
63 
64 #if (NETCDF==2) /* time-slice file in NetCDF format */
65  use netcdf
66  use nc_check_m
67 #endif
68 
69  implicit none
70 
71  character(len=100), intent(in) :: filename
72  real(dp), intent(inout) :: z_sl
73 
74  logical, optional, intent(in) :: opt_flag_temp_age_only
75  integer(i1b), optional, dimension(0:JMAX,0:IMAX), &
76  intent(inout) :: opt_maske, opt_n_cts
77  integer(i4b), optional, dimension(0:JMAX,0:IMAX), &
78  intent(inout) :: opt_kc_cts
79  real(dp), optional, dimension(0:JMAX,0:IMAX), &
80  intent(inout) :: opt_H_cold, opt_H_temp, opt_H
81  real(dp), optional, dimension(0:KRMAX,0:JMAX,0:IMAX), &
82  intent(inout) :: opt_temp_r
83  real(dp), optional, dimension(0:KTMAX,0:JMAX,0:IMAX), &
84  intent(inout) :: opt_omega_t, opt_age_t
85  real(dp), optional, dimension(0:KCMAX,0:JMAX,0:IMAX), &
86  intent(inout) :: opt_temp_c, opt_age_c, opt_omega_c
87 
88  integer(i4b) :: i, j, kc, kt, kr, n
89  character(len=256) :: filename_with_path
90  logical :: flag_temp_age_only
91 
92 #if (NETCDF==1) /* time-slice file in native binary format */
93 
94  integer(i4b) :: ios
95  character(len=256) :: ch_attr_title, ch_attr_institution, ch_attr_source, &
96  ch_attr_history, ch_attr_references
97 
98 #elif (NETCDF==2) /* time-slice file in NetCDF format */
99 
100  integer(i4b) :: ncid, ncv
101  ! ncid: ID of the output file
102  ! ncv: Variable ID
103 
104 #else
105 errormsg = ' >>> read_erg_nc: Parameter NETCDF must be either 1 or 2!'
106 call error(errormsg)
107 #endif
108 
109  integer(i1b), dimension(0:IMAX,0:JMAX) :: maske_conv, maske_old_conv, &
110  n_cts_conv, &
111  flag_shelfy_stream_x_conv, &
112  flag_shelfy_stream_y_conv, &
113  flag_shelfy_stream_conv, &
114  flag_grounding_line_1_conv, &
115  flag_grounding_line_2_conv, &
116  flag_calving_front_1_conv, &
117  flag_calving_front_2_conv, &
118  flag_grounded_front_a_1_conv, &
119  flag_grounded_front_a_2_conv, &
120  flag_grounded_front_b_1_conv, &
121  flag_grounded_front_b_2_conv
122  integer(i4b), dimension(0:IMAX,0:JMAX) :: kc_cts_conv
123  real(sp) :: time_conv, dummy_conv, z_sl_conv, &
124  V_tot_conv, A_grounded_conv, A_floating_conv, &
125  H_R_conv, &
126  xi_conv(0:imax), eta_conv(0:jmax), &
127  sigma_level_c_conv(0:kcmax), sigma_level_t_conv(0:ktmax), &
128  sigma_level_r_conv(0:krmax)
129  real(sp), dimension(0:IMAX,0:JMAX) :: lambda_conv, phi_conv, &
130  lon_conv, lat_conv, &
131  temp_s_conv, accum_conv, &
132  snowfall_conv, rainfall_conv, pdd_conv, &
133  as_perp_conv, as_perp_apl_conv, smb_corr_conv, &
134  q_geo_conv, &
135  zs_conv, zm_conv, zb_conv, zl_conv, zl0_conv, &
136  H_cold_conv, H_temp_conv, H_conv, &
137  Q_bm_conv, Q_tld_conv, &
138  am_perp_conv, &
139  qx_conv, qy_conv, &
140  vx_m_sia_conv, vy_m_sia_conv, vx_m_ssa_conv, vy_m_ssa_conv, &
141  dzs_dtau_conv, dzm_dtau_conv, dzb_dtau_conv, dzl_dtau_conv, &
142  dH_c_dtau_conv, dH_t_dtau_conv, dH_dtau_conv, &
143  vx_b_g_conv, vy_b_g_conv, vz_b_conv, vh_b_conv, &
144  vx_s_g_conv, vy_s_g_conv, vz_s_conv, vh_s_conv, &
145  vx_m_g_conv, vy_m_g_conv, vh_m_conv, &
146  temp_b_conv, temph_b_conv, &
147  tau_b_driving_conv, tau_b_drag_conv, &
148  p_b_w_conv, q_w_conv, q_w_x_conv, q_w_y_conv, H_w_conv, &
149  q_gl_g_conv, q_cf_g_conv, &
150  ratio_sl_x_conv, ratio_sl_y_conv, &
151  vis_ave_g_conv, vis_int_g_conv
152  real(sp), dimension(0:IMAX,0:JMAX,0:KCMAX) :: vx_c_conv, vy_c_conv, vz_c_conv, &
153  temp_c_conv, age_c_conv, &
154  enth_c_conv, omega_c_conv, &
155  enh_c_conv
156  real(sp), dimension(0:IMAX,0:JMAX,0:KTMAX) :: vx_t_conv, vy_t_conv, vz_t_conv, &
157  omega_t_conv, age_t_conv, &
158  enth_t_conv, &
159  enh_t_conv
160  real(sp), dimension(0:IMAX,0:JMAX,0:KRMAX) :: temp_r_conv
161 
162 #if (DISC>0) /* Ice discharge parameterisation */
163  integer(i1b), dimension(0:IMAX,0:JMAX) :: mask_mar_conv
164  real(sp), dimension(0:IMAX,0:JMAX) :: dis_perp_conv, &
165  cst_dist_conv, cos_grad_tc_conv
166 #endif
167 
168  real(dp) :: visc_min, visc_max, visc_init
169  logical :: flag_vis_ave_g
170 
171 !-------- Read data from time-slice file of previous simulation --------
172 
173  if (present(opt_flag_temp_age_only)) then
174  flag_temp_age_only = opt_flag_temp_age_only
175  else
176  flag_temp_age_only = .false.
177  end if
178 
179  filename_with_path = trim(anfdatpath)//'/'//trim(filename)
180 
181 #if (NETCDF==1) /* time-slice file in native binary format */
182 
183 #if !defined(ALLOW_OPENAD) /* Normal */
184  open(unit=11, iostat=ios, file=trim(filename_with_path), status='old', &
185  form='unformatted')
186 #else /* OpenAD */
187  open(unit=11, iostat=ios, file=trim(filename_with_path))
188 #endif /* Normal vs. OpenAD */
189 
190  if (ios /= 0) then
191  errormsg = ' >>> read_erg_nc: Error when opening the time-slice file!'
192  call error(errormsg)
193  end if
194 
195  read(unit=11) ch_attr_title
196  read(unit=11) ch_attr_institution
197  read(unit=11) ch_attr_source
198  read(unit=11) ch_attr_history
199  read(unit=11) ch_attr_references
200 
201  read(unit=11) time_conv
202  read(unit=11) dummy_conv ! this is either delta_ts or glac_index; not needed
203  read(unit=11) z_sl_conv
204 
205  read(unit=11) dummy_conv ! this is V_tot; not needed
206  read(unit=11) dummy_conv ! this is V_af; not needed
207  read(unit=11) dummy_conv ! this is A_grounded; not needed
208  read(unit=11) dummy_conv ! this is A_floating; not needed
209 
210  read(unit=11) xi_conv
211  read(unit=11) eta_conv
212  read(unit=11) sigma_level_c_conv
213  read(unit=11) sigma_level_t_conv
214  read(unit=11) sigma_level_r_conv
215  read(unit=11) lon_conv
216  read(unit=11) lat_conv
217  read(unit=11) lambda_conv
218  read(unit=11) phi_conv
219  read(unit=11) temp_s_conv
220  read(unit=11) accum_conv
221  read(unit=11) snowfall_conv
222  read(unit=11) rainfall_conv
223  read(unit=11) pdd_conv
224  read(unit=11) as_perp_conv
225  read(unit=11) as_perp_apl_conv
226  read(unit=11) smb_corr_conv
227 
228 #if (DISC>0) /* Ice discharge parameterisation */
229 
230  read(unit=11) dis_perp_conv
231  read(unit=11) cst_dist_conv
232  read(unit=11) cos_grad_tc_conv
233  read(unit=11) mask_mar_conv
234 
235 #endif
236 
237  read(unit=11) q_geo_conv
238  read(unit=11) maske_conv
239  read(unit=11) maske_old_conv
240  read(unit=11) n_cts_conv
241  read(unit=11) kc_cts_conv
242  read(unit=11) zs_conv
243  read(unit=11) zm_conv
244  read(unit=11) zb_conv
245  read(unit=11) zl_conv
246  read(unit=11) zl0_conv
247  read(unit=11) h_cold_conv
248  read(unit=11) h_temp_conv
249  read(unit=11) h_conv
250  read(unit=11) h_r_conv
251  read(unit=11) q_bm_conv
252  read(unit=11) q_tld_conv
253  read(unit=11) am_perp_conv
254  read(unit=11) qx_conv
255  read(unit=11) qy_conv
256  read(unit=11) vx_m_sia_conv
257  read(unit=11) vy_m_sia_conv
258  read(unit=11) vx_m_ssa_conv
259  read(unit=11) vy_m_ssa_conv
260  read(unit=11) dzs_dtau_conv
261  read(unit=11) dzm_dtau_conv
262  read(unit=11) dzb_dtau_conv
263  read(unit=11) dzl_dtau_conv
264  read(unit=11) dh_c_dtau_conv
265  read(unit=11) dh_t_dtau_conv
266  read(unit=11) dh_dtau_conv
267  read(unit=11) vx_b_g_conv
268  read(unit=11) vy_b_g_conv
269  read(unit=11) vz_b_conv
270  read(unit=11) vh_b_conv
271  read(unit=11) vx_s_g_conv
272  read(unit=11) vy_s_g_conv
273  read(unit=11) vz_s_conv
274  read(unit=11) vh_s_conv
275  read(unit=11) vx_m_g_conv
276  read(unit=11) vy_m_g_conv
277  read(unit=11) vh_m_conv
278  read(unit=11) temp_b_conv
279  read(unit=11) temph_b_conv
280  read(unit=11) tau_b_driving_conv
281  read(unit=11) tau_b_drag_conv
282  read(unit=11) p_b_w_conv
283  read(unit=11) q_w_conv
284  read(unit=11) q_w_x_conv
285  read(unit=11) q_w_y_conv
286  read(unit=11) h_w_conv
287  read(unit=11) q_gl_g_conv
288  read(unit=11) q_cf_g_conv
289  read(unit=11) ratio_sl_x_conv
290  read(unit=11) ratio_sl_y_conv
291  read(unit=11) flag_shelfy_stream_x_conv
292  read(unit=11) flag_shelfy_stream_y_conv
293  read(unit=11) flag_shelfy_stream_conv
294  read(unit=11) flag_grounding_line_1_conv
295  read(unit=11) flag_grounding_line_2_conv
296  read(unit=11) flag_calving_front_1_conv
297  read(unit=11) flag_calving_front_2_conv
298  read(unit=11) flag_grounded_front_a_1_conv
299  read(unit=11) flag_grounded_front_a_2_conv
300  read(unit=11) flag_grounded_front_b_1_conv
301  read(unit=11) flag_grounded_front_b_2_conv
302  read(unit=11) vis_ave_g_conv; flag_vis_ave_g = .true.
303  ! if not present in file, error message will occur
304  read(unit=11) vis_int_g_conv
305 
306  read(unit=11) vx_c_conv
307  read(unit=11) vy_c_conv
308  read(unit=11) vz_c_conv
309  read(unit=11) vx_t_conv
310  read(unit=11) vy_t_conv
311  read(unit=11) vz_t_conv
312  read(unit=11) temp_c_conv
313  read(unit=11) omega_t_conv
314  read(unit=11) temp_r_conv
315  read(unit=11) enth_c_conv
316  read(unit=11) enth_t_conv
317  read(unit=11) omega_c_conv
318  read(unit=11) enh_c_conv
319  read(unit=11) enh_t_conv
320  read(unit=11) age_c_conv
321  read(unit=11) age_t_conv
322 
323  close(unit=11,status='keep')
324 
325 #elif (NETCDF==2) /* time-slice file in NetCDF format */
326 
327  call check( nf90_open(trim(filename_with_path), nf90_nowrite, ncid) )
328 
329  call check( nf90_inq_varid(ncid, 'time', ncv) )
330  call check( nf90_get_var(ncid, ncv, time_conv) )
331 
332  if ( nf90_inq_varid(ncid, 'delta_ts', ncv) == nf90_noerr ) then
333  call check( nf90_get_var(ncid, ncv, dummy_conv) )
334  else if ( nf90_inq_varid(ncid, 'glac_index', ncv) == nf90_noerr ) then
335  call check( nf90_get_var(ncid, ncv, dummy_conv) )
336  else
337  dummy_conv = 0.0_sp
338  end if
339 
340  call check( nf90_inq_varid(ncid, 'z_sl', ncv) )
341  call check( nf90_get_var(ncid, ncv, z_sl_conv) )
342 
343  call check( nf90_inq_varid(ncid, 'x', ncv) )
344  call check( nf90_get_var(ncid, ncv, xi_conv) )
345 
346  call check( nf90_inq_varid(ncid, 'y', ncv) )
347  call check( nf90_get_var(ncid, ncv, eta_conv) )
348 
349  call check( nf90_inq_varid(ncid, 'sigma_level_c', ncv) )
350  call check( nf90_get_var(ncid, ncv, sigma_level_c_conv) )
351 
352  call check( nf90_inq_varid(ncid, 'sigma_level_t', ncv) )
353  call check( nf90_get_var(ncid, ncv, sigma_level_t_conv) )
354 
355  call check( nf90_inq_varid(ncid, 'sigma_level_r', ncv) )
356  call check( nf90_get_var(ncid, ncv, sigma_level_r_conv) )
357 
358  call check( nf90_inq_varid(ncid, 'lon', ncv) )
359  call check( nf90_get_var(ncid, ncv, lon_conv) )
360 
361  call check( nf90_inq_varid(ncid, 'lat', ncv) )
362  call check( nf90_get_var(ncid, ncv, lat_conv) )
363 
364  call check( nf90_inq_varid(ncid, 'lambda', ncv) )
365  call check( nf90_get_var(ncid, ncv, lambda_conv) )
366 
367  call check( nf90_inq_varid(ncid, 'phi', ncv) )
368  call check( nf90_get_var(ncid, ncv, phi_conv) )
369 
370  call check( nf90_inq_varid(ncid, 'temp_s', ncv) )
371  call check( nf90_get_var(ncid, ncv, temp_s_conv) )
372 
373  call check( nf90_inq_varid(ncid, 'prec', ncv) )
374  call check( nf90_get_var(ncid, ncv, accum_conv) )
375 
376  call check( nf90_inq_varid(ncid, 'snowfall', ncv) )
377  call check( nf90_get_var(ncid, ncv, snowfall_conv) )
378 
379  call check( nf90_inq_varid(ncid, 'rainfall', ncv) )
380  call check( nf90_get_var(ncid, ncv, rainfall_conv) )
381 
382  call check( nf90_inq_varid(ncid, 'pdd', ncv) )
383  call check( nf90_get_var(ncid, ncv, pdd_conv) )
384 
385  call check( nf90_inq_varid(ncid, 'as_perp', ncv) )
386  call check( nf90_get_var(ncid, ncv, as_perp_conv) )
387 
388  call check( nf90_inq_varid(ncid, 'as_perp_apl', ncv) )
389  call check( nf90_get_var(ncid, ncv, as_perp_apl_conv) )
390 
391  call check( nf90_inq_varid(ncid, 'smb_corr', ncv) )
392  call check( nf90_get_var(ncid, ncv, smb_corr_conv) )
393 
394 #if (DISC>0) /* Ice discharge parameterisation */
395 
396  call check( nf90_inq_varid(ncid, 'dis_perp', ncv) )
397  call check( nf90_get_var(ncid, ncv, dis_perp_conv) )
398 
399  call check( nf90_inq_varid(ncid, 'cst_dist', ncv) )
400  call check( nf90_get_var(ncid, ncv, cst_dist_conv) )
401 
402  call check( nf90_inq_varid(ncid, 'cos_grad_tc', ncv) )
403  call check( nf90_get_var(ncid, ncv, cos_grad_tc_conv) )
404 
405  call check( nf90_inq_varid(ncid, 'mask_mar', ncv) )
406  call check( nf90_get_var(ncid, ncv, mask_mar_conv) )
407 
408 #endif
409 
410  call check( nf90_inq_varid(ncid, 'q_geo', ncv) )
411  call check( nf90_get_var(ncid, ncv, q_geo_conv) )
412 
413  call check( nf90_inq_varid(ncid, 'maske', ncv) )
414  call check( nf90_get_var(ncid, ncv, maske_conv) )
415 
416  call check( nf90_inq_varid(ncid, 'maske_old', ncv) )
417  call check( nf90_get_var(ncid, ncv, maske_old_conv) )
418 
419  call check( nf90_inq_varid(ncid, 'n_cts', ncv) )
420  call check( nf90_get_var(ncid, ncv, n_cts_conv) )
421 
422  call check( nf90_inq_varid(ncid, 'kc_cts', ncv) )
423  call check( nf90_get_var(ncid, ncv, kc_cts_conv) )
424 
425  call check( nf90_inq_varid(ncid, 'zs', ncv) )
426  call check( nf90_get_var(ncid, ncv, zs_conv) )
427 
428  call check( nf90_inq_varid(ncid, 'zm', ncv) )
429  call check( nf90_get_var(ncid, ncv, zm_conv) )
430 
431  call check( nf90_inq_varid(ncid, 'zb', ncv) )
432  call check( nf90_get_var(ncid, ncv, zb_conv) )
433 
434  call check( nf90_inq_varid(ncid, 'zl', ncv) )
435  call check( nf90_get_var(ncid, ncv, zl_conv) )
436 
437  call check( nf90_inq_varid(ncid, 'zl0', ncv) )
438  call check( nf90_get_var(ncid, ncv, zl0_conv) )
439 
440  call check( nf90_inq_varid(ncid, 'H_cold', ncv) )
441  call check( nf90_get_var(ncid, ncv, h_cold_conv) )
442 
443  call check( nf90_inq_varid(ncid, 'H_temp', ncv) )
444  call check( nf90_get_var(ncid, ncv, h_temp_conv) )
445 
446  call check( nf90_inq_varid(ncid, 'H', ncv) )
447  call check( nf90_get_var(ncid, ncv, h_conv) )
448 
449  call check( nf90_inq_varid(ncid, 'H_R', ncv) )
450  call check( nf90_get_var(ncid, ncv, h_r_conv) )
451 
452  call check( nf90_inq_varid(ncid, 'Q_bm', ncv) )
453  call check( nf90_get_var(ncid, ncv, q_bm_conv) )
454 
455  call check( nf90_inq_varid(ncid, 'Q_tld', ncv) )
456  call check( nf90_get_var(ncid, ncv, q_tld_conv) )
457 
458  call check( nf90_inq_varid(ncid, 'am_perp', ncv) )
459  call check( nf90_get_var(ncid, ncv, am_perp_conv) )
460 
461  call check( nf90_inq_varid(ncid, 'qx', ncv) )
462  call check( nf90_get_var(ncid, ncv, qx_conv) )
463 
464  call check( nf90_inq_varid(ncid, 'qy', ncv) )
465  call check( nf90_get_var(ncid, ncv, qy_conv) )
466 
467  if ( nf90_inq_varid(ncid, 'vx_m_sia', ncv) == nf90_noerr ) then
468  call check( nf90_get_var(ncid, ncv, vx_m_sia_conv) )
469  else
470  write(6,'(/1x,a)') '>>> read_erg_nc: Variable vx_m_sia'
471  write(6, '(1x,a)') ' not available in read file *.nc.'
472  vx_m_sia_conv = 0.0_sp
473  end if
474 
475  if ( nf90_inq_varid(ncid, 'vy_m_sia', ncv) == nf90_noerr ) then
476  call check( nf90_get_var(ncid, ncv, vy_m_sia_conv) )
477  else
478  write(6,'(/1x,a)') '>>> read_erg_nc: Variable vy_m_sia'
479  write(6, '(1x,a)') ' not available in read file *.nc.'
480  vy_m_sia_conv = 0.0_sp
481  end if
482 
483  if ( nf90_inq_varid(ncid, 'vx_m_ssa', ncv) == nf90_noerr ) then
484  call check( nf90_get_var(ncid, ncv, vx_m_ssa_conv) )
485  else
486  write(6,'(/1x,a)') '>>> read_erg_nc: Variable vx_m_ssa'
487  write(6, '(1x,a)') ' not available in read file *.nc.'
488  vx_m_ssa_conv = 0.0_sp
489  end if
490 
491  if ( nf90_inq_varid(ncid, 'vy_m_ssa', ncv) == nf90_noerr ) then
492  call check( nf90_get_var(ncid, ncv, vy_m_ssa_conv) )
493  else
494  write(6,'(/1x,a)') '>>> read_erg_nc: Variable vy_m_ssa'
495  write(6, '(1x,a)') ' not available in read file *.nc.'
496  vy_m_ssa_conv = 0.0_sp
497  end if
498 
499  call check( nf90_inq_varid(ncid, 'dzs_dt', ncv) )
500  call check( nf90_get_var(ncid, ncv, dzs_dtau_conv) )
501 
502  call check( nf90_inq_varid(ncid, 'dzm_dt', ncv) )
503  call check( nf90_get_var(ncid, ncv, dzm_dtau_conv) )
504 
505  call check( nf90_inq_varid(ncid, 'dzb_dt', ncv) )
506  call check( nf90_get_var(ncid, ncv, dzb_dtau_conv) )
507 
508  call check( nf90_inq_varid(ncid, 'dzl_dt', ncv) )
509  call check( nf90_get_var(ncid, ncv, dzl_dtau_conv) )
510 
511  call check( nf90_inq_varid(ncid, 'dH_c_dt', ncv) )
512  call check( nf90_get_var(ncid, ncv, dh_c_dtau_conv) )
513 
514  call check( nf90_inq_varid(ncid, 'dH_t_dt', ncv) )
515  call check( nf90_get_var(ncid, ncv, dh_t_dtau_conv) )
516 
517  call check( nf90_inq_varid(ncid, 'dH_dt', ncv) )
518  call check( nf90_get_var(ncid, ncv, dh_dtau_conv) )
519 
520  call check( nf90_inq_varid(ncid, 'vx_b_g', ncv) )
521  call check( nf90_get_var(ncid, ncv, vx_b_g_conv) )
522 
523  call check( nf90_inq_varid(ncid, 'vy_b_g', ncv) )
524  call check( nf90_get_var(ncid, ncv, vy_b_g_conv) )
525 
526  call check( nf90_inq_varid(ncid, 'vz_b', ncv) )
527  call check( nf90_get_var(ncid, ncv, vz_b_conv) )
528 
529  call check( nf90_inq_varid(ncid, 'vh_b', ncv) )
530  call check( nf90_get_var(ncid, ncv, vh_b_conv) )
531 
532  call check( nf90_inq_varid(ncid, 'vx_s_g', ncv) )
533  call check( nf90_get_var(ncid, ncv, vx_s_g_conv) )
534 
535  call check( nf90_inq_varid(ncid, 'vy_s_g', ncv) )
536  call check( nf90_get_var(ncid, ncv, vy_s_g_conv) )
537 
538  call check( nf90_inq_varid(ncid, 'vz_s', ncv) )
539  call check( nf90_get_var(ncid, ncv, vz_s_conv) )
540 
541  call check( nf90_inq_varid(ncid, 'vh_s', ncv) )
542  call check( nf90_get_var(ncid, ncv, vh_s_conv) )
543 
544  call check( nf90_inq_varid(ncid, 'vx_m_g', ncv) )
545  call check( nf90_get_var(ncid, ncv, vx_m_g_conv) )
546 
547  call check( nf90_inq_varid(ncid, 'vy_m_g', ncv) )
548  call check( nf90_get_var(ncid, ncv, vy_m_g_conv) )
549 
550  call check( nf90_inq_varid(ncid, 'vh_m', ncv) )
551  call check( nf90_get_var(ncid, ncv, vh_m_conv) )
552 
553  call check( nf90_inq_varid(ncid, 'temp_b', ncv) )
554  call check( nf90_get_var(ncid, ncv, temp_b_conv) )
555 
556  call check( nf90_inq_varid(ncid, 'temph_b', ncv) )
557  call check( nf90_get_var(ncid, ncv, temph_b_conv) )
558 
559  call check( nf90_inq_varid(ncid, 'tau_b_driving', ncv) )
560  call check( nf90_get_var(ncid, ncv, tau_b_driving_conv) )
561 
562  call check( nf90_inq_varid(ncid, 'tau_b_drag', ncv) )
563  call check( nf90_get_var(ncid, ncv, tau_b_drag_conv) )
564 
565  call check( nf90_inq_varid(ncid, 'p_b_w', ncv) )
566  call check( nf90_get_var(ncid, ncv, p_b_w_conv) )
567 
568  call check( nf90_inq_varid(ncid, 'q_w', ncv) )
569  call check( nf90_get_var(ncid, ncv, q_w_conv) )
570 
571  call check( nf90_inq_varid(ncid, 'q_w_x', ncv) )
572  call check( nf90_get_var(ncid, ncv, q_w_x_conv) )
573 
574  call check( nf90_inq_varid(ncid, 'q_w_y', ncv) )
575  call check( nf90_get_var(ncid, ncv, q_w_y_conv) )
576 
577  call check( nf90_inq_varid(ncid, 'H_w', ncv) )
578  call check( nf90_get_var(ncid, ncv, h_w_conv) )
579 
580  call check( nf90_inq_varid(ncid, 'q_gl_g', ncv) )
581  call check( nf90_get_var(ncid, ncv, q_gl_g_conv) )
582 
583  call check( nf90_inq_varid(ncid, 'q_cf_g', ncv) )
584  call check( nf90_get_var(ncid, ncv, q_cf_g_conv) )
585 
586  call check( nf90_inq_varid(ncid, 'ratio_sl_x', ncv) )
587  call check( nf90_get_var(ncid, ncv, ratio_sl_x_conv) )
588 
589  call check( nf90_inq_varid(ncid, 'ratio_sl_y', ncv) )
590  call check( nf90_get_var(ncid, ncv, ratio_sl_y_conv) )
591 
592  call check( nf90_inq_varid(ncid, 'flag_shelfy_stream_x', ncv) )
593  call check( nf90_get_var(ncid, ncv, flag_shelfy_stream_x_conv) )
594 
595  call check( nf90_inq_varid(ncid, 'flag_shelfy_stream_y', ncv) )
596  call check( nf90_get_var(ncid, ncv, flag_shelfy_stream_y_conv) )
597 
598  call check( nf90_inq_varid(ncid, 'flag_shelfy_stream', ncv) )
599  call check( nf90_get_var(ncid, ncv, flag_shelfy_stream_conv) )
600 
601  call check( nf90_inq_varid(ncid, 'flag_grounding_line_1', ncv) )
602  call check( nf90_get_var(ncid, ncv, flag_grounding_line_1_conv) )
603 
604  call check( nf90_inq_varid(ncid, 'flag_grounding_line_2', ncv) )
605  call check( nf90_get_var(ncid, ncv, flag_grounding_line_2_conv) )
606 
607  call check( nf90_inq_varid(ncid, 'flag_calving_front_1', ncv) )
608  call check( nf90_get_var(ncid, ncv, flag_calving_front_1_conv) )
609 
610  call check( nf90_inq_varid(ncid, 'flag_calving_front_2', ncv) )
611  call check( nf90_get_var(ncid, ncv, flag_calving_front_2_conv) )
612 
613  call check( nf90_inq_varid(ncid, 'flag_grounded_front_a_1', ncv) )
614  call check( nf90_get_var(ncid, ncv, flag_grounded_front_a_1_conv) )
615 
616  call check( nf90_inq_varid(ncid, 'flag_grounded_front_a_2', ncv) )
617  call check( nf90_get_var(ncid, ncv, flag_grounded_front_a_2_conv) )
618 
619  call check( nf90_inq_varid(ncid, 'flag_grounded_front_b_1', ncv) )
620  call check( nf90_get_var(ncid, ncv, flag_grounded_front_b_1_conv) )
621 
622  call check( nf90_inq_varid(ncid, 'flag_grounded_front_b_2', ncv) )
623  call check( nf90_get_var(ncid, ncv, flag_grounded_front_b_2_conv) )
624 
625  if ( nf90_inq_varid(ncid, 'vis_ave_g', ncv) == nf90_noerr ) then
626  call check( nf90_get_var(ncid, ncv, vis_ave_g_conv) )
627  flag_vis_ave_g = .true.
628  else
629  write(6,'(/1x,a)') '>>> read_erg_nc: Variable vis_ave_g'
630  write(6, '(1x,a)') ' not available in read file *.nc.'
631  vis_ave_g_conv = 0.0_sp
632  flag_vis_ave_g = .false.
633  end if
634 
635  call check( nf90_inq_varid(ncid, 'vis_int_g', ncv) )
636  call check( nf90_get_var(ncid, ncv, vis_int_g_conv) )
637 
638  call check( nf90_inq_varid(ncid, 'vx_c', ncv) )
639  call check( nf90_get_var(ncid, ncv, vx_c_conv) )
640 
641  call check( nf90_inq_varid(ncid, 'vy_c', ncv) )
642  call check( nf90_get_var(ncid, ncv, vy_c_conv) )
643 
644  call check( nf90_inq_varid(ncid, 'vz_c', ncv) )
645  call check( nf90_get_var(ncid, ncv, vz_c_conv) )
646 
647  call check( nf90_inq_varid(ncid, 'vx_t', ncv) )
648  call check( nf90_get_var(ncid, ncv, vx_t_conv) )
649 
650  call check( nf90_inq_varid(ncid, 'vy_t', ncv) )
651  call check( nf90_get_var(ncid, ncv, vy_t_conv) )
652 
653  call check( nf90_inq_varid(ncid, 'vz_t', ncv) )
654  call check( nf90_get_var(ncid, ncv, vz_t_conv) )
655 
656  call check( nf90_inq_varid(ncid, 'temp_c', ncv) )
657  call check( nf90_get_var(ncid, ncv, temp_c_conv) )
658 
659  call check( nf90_inq_varid(ncid, 'omega_t', ncv) )
660  call check( nf90_get_var(ncid, ncv, omega_t_conv) )
661 
662  call check( nf90_inq_varid(ncid, 'temp_r', ncv) )
663  call check( nf90_get_var(ncid, ncv, temp_r_conv) )
664 
665  call check( nf90_inq_varid(ncid, 'enth_c', ncv) )
666  call check( nf90_get_var(ncid, ncv, enth_c_conv) )
667 
668  call check( nf90_inq_varid(ncid, 'enth_t', ncv) )
669  call check( nf90_get_var(ncid, ncv, enth_t_conv) )
670 
671  call check( nf90_inq_varid(ncid, 'omega_c', ncv) )
672  call check( nf90_get_var(ncid, ncv, omega_c_conv) )
673 
674  call check( nf90_inq_varid(ncid, 'enh_c', ncv) )
675  call check( nf90_get_var(ncid, ncv, enh_c_conv) )
676 
677  call check( nf90_inq_varid(ncid, 'enh_t', ncv) )
678  call check( nf90_get_var(ncid, ncv, enh_t_conv) )
679 
680  call check( nf90_inq_varid(ncid, 'age_c', ncv) )
681  call check( nf90_get_var(ncid, ncv, age_c_conv) )
682 
683  call check( nf90_inq_varid(ncid, 'age_t', ncv) )
684  call check( nf90_get_var(ncid, ncv, age_t_conv) )
685 
686  call check( nf90_close(ncid) )
687 
688 #endif
689 
690 !-------- Convert data to real*8 and years to seconds --------
691 
692  if (.not.flag_temp_age_only) then
693 
694  z_sl = real(z_sl_conv,dp)
695 
696  h_r = real(h_r_conv,dp)
697 
698  do i=0, imax
699  xi(i) = real(xi_conv(i),dp)
700  end do
701 
702  do j=0, jmax
703  eta(j) = real(eta_conv(j),dp)
704  end do
705 
706  do i=0, imax
707  do j=0, jmax
708 
709  maske(j,i) = maske_conv(i,j)
710  maske_old(j,i) = maske_old_conv(i,j)
711  n_cts(j,i) = n_cts_conv(i,j)
712  kc_cts(j,i) = kc_cts_conv(i,j)
713  zs(j,i) = real(zs_conv(i,j),dp)
714  zm(j,i) = real(zm_conv(i,j),dp)
715  zb(j,i) = real(zb_conv(i,j),dp)
716  zl(j,i) = real(zl_conv(i,j),dp)
717 #if (CALCMOD==1)
718  h_c(j,i) = real(H_cold_conv(i,j),dp)
719  h_t(j,i) = real(H_temp_conv(i,j),dp)
720 #elif (CALCMOD==0 || CALCMOD==2 || CALCMOD==3 || CALCMOD==-1)
721  h_c(j,i) = real(H_conv(i,j),dp)
722  h_t(j,i) = 0.0_dp
723 #endif
724  q_bm(j,i) = real(Q_bm_conv(i,j),dp)/year_sec
725  q_tld(j,i) = real(Q_tld_conv(i,j),dp)/year_sec
726  am_perp(j,i) = real(am_perp_conv(i,j),dp)/year_sec
727  qx(j,i) = real(qx_conv(i,j),dp)/year_sec
728  qy(j,i) = real(qy_conv(i,j),dp)/year_sec
729  vx_m_sia(j,i) = real(vx_m_sia_conv(i,j),dp)/year_sec
730  vy_m_sia(j,i) = real(vy_m_sia_conv(i,j),dp)/year_sec
731  vx_m_ssa(j,i) = real(vx_m_ssa_conv(i,j),dp)/year_sec
732  vy_m_ssa(j,i) = real(vy_m_ssa_conv(i,j),dp)/year_sec
733  dzs_dtau(j,i) = real(dzs_dtau_conv(i,j),dp)/year_sec
734  dzm_dtau(j,i) = real(dzm_dtau_conv(i,j),dp)/year_sec
735  dzb_dtau(j,i) = real(dzb_dtau_conv(i,j),dp)/year_sec
736  dzl_dtau(j,i) = real(dzl_dtau_conv(i,j),dp)/year_sec
737  dh_c_dtau(j,i) = real(dH_c_dtau_conv(i,j),dp)/year_sec
738  dh_t_dtau(j,i) = real(dH_t_dtau_conv(i,j),dp)/year_sec
739  vx_b_g(j,i) = real(vx_b_g_conv(i,j),dp)/year_sec
740  vy_b_g(j,i) = real(vy_b_g_conv(i,j),dp)/year_sec
741  vz_b(j,i) = real(vz_b_conv(i,j),dp)/year_sec
742  vx_s_g(j,i) = real(vx_s_g_conv(i,j),dp)/year_sec
743  vy_s_g(j,i) = real(vy_s_g_conv(i,j),dp)/year_sec
744  vz_s(j,i) = real(vz_s_conv(i,j),dp)/year_sec
745  temp_b(j,i) = real(temp_b_conv(i,j),dp)
746  temph_b(j,i) = real(temph_b_conv(i,j),dp)
747  p_b_w(j,i) = real(p_b_w_conv(i,j),dp)
748  q_w(j,i) = real(q_w_conv(i,j),dp)/year_sec
749  q_w_x(j,i) = real(q_w_x_conv(i,j),dp)/year_sec
750  q_w_y(j,i) = real(q_w_y_conv(i,j),dp)/year_sec
751  h_w(j,i) = real(H_w_conv(i,j),dp)
752  ratio_sl_x(j,i) = real(ratio_sl_x_conv(i,j),dp)
753  ratio_sl_y(j,i) = real(ratio_sl_y_conv(i,j),dp)
754 
755  if (flag_shelfy_stream_x_conv(i,j) == 1_i1b) then
756  flag_shelfy_stream_x(j,i) = .true.
757  else
758  flag_shelfy_stream_x(j,i) = .false.
759  end if
760 
761  if (flag_shelfy_stream_y_conv(i,j) == 1_i1b) then
762  flag_shelfy_stream_y(j,i) = .true.
763  else
764  flag_shelfy_stream_y(j,i) = .false.
765  end if
766 
767  if (flag_shelfy_stream_conv(i,j) == 1_i1b) then
768  flag_shelfy_stream(j,i) = .true.
769  else
770  flag_shelfy_stream(j,i) = .false.
771  end if
772 
773  if (flag_grounding_line_1_conv(i,j) == 1_i1b) then
774  flag_grounding_line_1(j,i) = .true.
775  else
776  flag_grounding_line_1(j,i) = .false.
777  end if
778 
779  if (flag_grounding_line_2_conv(i,j) == 1_i1b) then
780  flag_grounding_line_2(j,i) = .true.
781  else
782  flag_grounding_line_2(j,i) = .false.
783  end if
784 
785  if (flag_calving_front_1_conv(i,j) == 1_i1b) then
786  flag_calving_front_1(j,i) = .true.
787  else
788  flag_calving_front_1(j,i) = .false.
789  end if
790 
791  if (flag_calving_front_2_conv(i,j) == 1_i1b) then
792  flag_calving_front_2(j,i) = .true.
793  else
794  flag_calving_front_2(j,i) = .false.
795  end if
796 
797  if (flag_grounded_front_a_1_conv(i,j) == 1_i1b) then
798  flag_grounded_front_a_1(j,i) = .true.
799  else
800  flag_grounded_front_a_1(j,i) = .false.
801  end if
802 
803  if (flag_grounded_front_a_2_conv(i,j) == 1_i1b) then
804  flag_grounded_front_a_2(j,i) = .true.
805  else
806  flag_grounded_front_a_2(j,i) = .false.
807  end if
808 
809  if (flag_grounded_front_b_1_conv(i,j) == 1_i1b) then
810  flag_grounded_front_b_1(j,i) = .true.
811  else
812  flag_grounded_front_b_1(j,i) = .false.
813  end if
814 
815  if (flag_grounded_front_b_2_conv(i,j) == 1_i1b) then
816  flag_grounded_front_b_2(j,i) = .true.
817  else
818  flag_grounded_front_b_2(j,i) = .false.
819  end if
820 
821  vis_ave_g(j,i) = real(vis_ave_g_conv(i,j),dp)
822  vis_int_g(j,i) = real(vis_int_g_conv(i,j),dp)
823 
824  do kr=0, krmax
825  temp_r(kr,j,i) = real(temp_r_conv(i,j,kr),dp)
826  end do
827 
828  do kt=0, ktmax
829  vx_t(kt,j,i) = real(vx_t_conv(i,j,kt),dp)/year_sec
830  vy_t(kt,j,i) = real(vy_t_conv(i,j,kt),dp)/year_sec
831  vz_t(kt,j,i) = real(vz_t_conv(i,j,kt),dp)/year_sec
832  omega_t(kt,j,i) = real(omega_t_conv(i,j,kt),dp)
833  age_t(kt,j,i) = real(age_t_conv(i,j,kt),dp)*year_sec
834  enth_t(kt,j,i) = real(enth_t_conv(i,j,kt),dp)
835  enh_t(kt,j,i) = real(enh_t_conv(i,j,kt),dp)
836  end do
837 
838  do kc=0, kcmax
839  vx_c(kc,j,i) = real(vx_c_conv(i,j,kc),dp)/year_sec
840  vy_c(kc,j,i) = real(vy_c_conv(i,j,kc),dp)/year_sec
841  vz_c(kc,j,i) = real(vz_c_conv(i,j,kc),dp)/year_sec
842  temp_c(kc,j,i) = real(temp_c_conv(i,j,kc),dp)
843  age_c(kc,j,i) = real(age_c_conv(i,j,kc),dp)*year_sec
844  enth_c(kc,j,i) = real(enth_c_conv(i,j,kc),dp)
845  omega_c(kc,j,i) = real(omega_c_conv(i,j,kc),dp)
846  enh_c(kc,j,i) = real(enh_c_conv(i,j,kc),dp)
847  end do
848 
849  end do
850  end do
851 
852  if (.not.flag_vis_ave_g) then ! reconstruct vis_ave_g from vis_int_g
853 
854 #if (defined(VISC_MIN) && defined(VISC_MAX))
855  visc_min = visc_min
856  visc_max = visc_max
857 #else
858  visc_min = 1.0e+10_dp ! Pa s
859  visc_max = 1.0e+25_dp ! Pa s
860 #endif
861 
862 #if (defined(VISC_INIT_SSA))
863  visc_init = visc_init_ssa
864 #else
865  visc_init = 1.0e+15_dp ! Pa s
866 #endif
867 
868  do i=0, imax
869  do j=0, jmax
870 
871  if ((maske(j,i)==0_i1b).or.(maske(j,i)==3_i1b)) then
872  ! grounded or floating ice
873 
874  vis_ave_g(j,i) = vis_int_g(j,i)/max((h_c(j,i)+h_t(j,i)), eps_dp)
875 
876  vis_ave_g(j,i) = max(min(vis_ave_g(j,i), visc_max), visc_min)
877 
878  else ! (maske(j,i)==1_i1b).or.(maske(j,i)==2_i1b),
879  ! ice-free land or ocean
880 
881  vis_ave_g(j,i) = visc_init ! dummy value
882 
883  end if
884 
885  end do
886  end do
887 
888  end if ! (.not.flag_vis_ave_g)
889 
890  else ! flag_temp_age_only == .true.
891 
892  if ( present(opt_maske) &
893  .and.present(opt_n_cts) &
894  .and.present(opt_kc_cts) &
895  .and.present(opt_h_cold) &
896  .and.present(opt_h_temp) &
897  .and.present(opt_h) &
898  .and.present(opt_temp_r) &
899  .and.present(opt_omega_t) &
900  .and.present(opt_age_t) &
901  .and.present(opt_temp_c) &
902  .and.present(opt_age_c) &
903  .and.present(opt_omega_c) &
904  ) then
905 
906  do i=0, imax
907  do j=0, jmax
908 
909  opt_maske(j,i) = maske_conv(i,j)
910  opt_n_cts(j,i) = n_cts_conv(i,j)
911  opt_kc_cts(j,i) = kc_cts_conv(i,j)
912 
913  opt_h_cold(j,i) = real(H_cold_conv(i,j),dp)
914  opt_h_temp(j,i) = real(H_temp_conv(i,j),dp)
915  opt_h(j,i) = real(H_conv(i,j),dp)
916 
917  do kr=0, krmax
918  opt_temp_r(kr,j,i) = real(temp_r_conv(i,j,kr),dp)
919  end do
920 
921  do kt=0, ktmax
922  opt_omega_t(kt,j,i) = real(omega_t_conv(i,j,kt),dp)
923  opt_age_t(kt,j,i) = real(age_t_conv(i,j,kt),dp)*year_sec
924  end do
925 
926  do kc=0, kcmax
927  opt_temp_c(kc,j,i) = real(temp_c_conv(i,j,kc),dp)
928  opt_age_c(kc,j,i) = real(age_c_conv(i,j,kc),dp)*year_sec
929  opt_omega_c(kc,j,i) = real(omega_c_conv(i,j,kc),dp)
930  end do
931 
932  end do
933  end do
934 
935  else
936 
937  errormsg = ' >>> read_erg_nc: optional argument(s) missing!'
938  call error(errormsg)
939 
940  end if
941 
942  end if ! flag_temp_age_only == .false./.true.
943 
944  end subroutine read_erg_nc
945 
946 !-------------------------------------------------------------------------------
947 !> Reading of the target-topography file (in NetCDF format).
948 !<------------------------------------------------------------------------------
949  subroutine read_target_topo_nc(target_topo_dat_name)
951  use sico_variables_m, only : maske_target, &
954 
955 #if (NETCDF==2)
956  use netcdf
957  use nc_check_m
958 #endif
959 
960  implicit none
961 
962  character(len=100), intent(in) :: target_topo_dat_name
963 
964 ! Return variables
965 ! (defined as global variables in module sico_variables_m):
966 !
967 ! maske_target, zs_target, zb_target, zl_target, H_target
968 
969  integer(i1b), dimension(0:IMAX,0:JMAX) :: maske_conv
970  real(sp), dimension(0:IMAX,0:JMAX) :: zs_conv, zb_conv, zl_conv, H_conv
971  character(len=256) :: target_topo_dat_path
972  character(len=256) :: filename_with_path
973 
974 #if (NETCDF==2)
975  integer(i4b) :: ncid, ncv
976  ! ncid: ID of the output file
977  ! ncv: Variable ID
978 #endif
979 
980  character(len=64), parameter :: thisroutine = 'read_target_topo_nc'
981 
982 #if defined(ALLOW_OPENAD) /* OpenAD */
983  integer(i4b) :: i, j
984 #endif /* OpenAD */
985 
986 !-------- Read target-topography data
987 ! (from a time-slice file of a previous simulation) --------
988 
989 #if (NETCDF==2)
990 
991 #if (defined(TARGET_TOPO_DAT_PATH))
992  target_topo_dat_path = trim(target_topo_dat_path)
993 #else
994  target_topo_dat_path = 'dummy'
995  errormsg = ' >>> read_target_topo_nc: TARGET_TOPO_DAT_PATH must be defined!'
996  call error(errormsg)
997 #endif
998 
999  filename_with_path &
1000  = trim(target_topo_dat_path)//'/'//trim(target_topo_dat_name)
1001 
1002  call check( nf90_open(trim(filename_with_path), nf90_nowrite, ncid), &
1003  thisroutine )
1004 
1005  call check( nf90_inq_varid(ncid, 'maske', ncv), thisroutine )
1006  call check( nf90_get_var(ncid, ncv, maske_conv), thisroutine )
1007 
1008  call check( nf90_inq_varid(ncid, 'zs', ncv), thisroutine )
1009  call check( nf90_get_var(ncid, ncv, zs_conv), thisroutine )
1010 
1011  call check( nf90_inq_varid(ncid, 'zb', ncv), thisroutine )
1012  call check( nf90_get_var(ncid, ncv, zb_conv), thisroutine )
1013 
1014  call check( nf90_inq_varid(ncid, 'zl', ncv), thisroutine )
1015  call check( nf90_get_var(ncid, ncv, zl_conv), thisroutine )
1016 
1017  call check( nf90_inq_varid(ncid, 'H', ncv), thisroutine )
1018  call check( nf90_get_var(ncid, ncv, h_conv), thisroutine )
1019 
1020  call check( nf90_close(ncid), thisroutine )
1021 
1022 #else /* NETCDF != 2, not allowed */
1023 
1024  maske_conv = 0_i1b ! dummy values
1025  zs_conv = 0.0_sp ! dummy values
1026  zb_conv = 0.0_sp ! dummy values
1027  zl_conv = 0.0_sp ! dummy values
1028  h_conv = 0.0_sp ! dummy values
1029 
1030  errormsg = ' >>> read_target_topo_nc: Parameter NETCDF must be set to 2!'
1031  call error(errormsg)
1032 
1033 #endif
1034 
1035 !-------- Convert data to double precision --------
1036 
1037 #if !defined(ALLOW_OPENAD) /* Normal */
1038  maske_target = transpose(maske_conv)
1039  zs_target = real(transpose(zs_conv),dp)
1040  zb_target = real(transpose(zb_conv),dp)
1041  zl_target = real(transpose(zl_conv),dp)
1042  h_target = real(transpose(H_conv) ,dp)
1043 #else /* OpenAD */
1044  do i=0, imax
1045  do j=0, jmax
1046  maske_target(j,i) = maske_conv(i,j)
1047  zs_target(j,i) = real(zs_conv(i,j),dp)
1048  zb_target(j,i) = real(zb_conv(i,j),dp)
1049  zl_target(j,i) = real(zl_conv(i,j),dp)
1050  h_target(j,i) = real(H_conv(i,j) ,dp)
1051  end do
1052  end do
1053 #endif /* Normal vs. OpenAD */
1054 
1055  end subroutine read_target_topo_nc
1056 
1057 !-------------------------------------------------------------------------------
1058 !> Reading of the tabulated kei function.
1059 !<------------------------------------------------------------------------------
1060  subroutine read_kei()
1063 
1064  implicit none
1065 
1066  integer(i4b) :: n, n_data_kei_max
1067  integer(i4b) :: ios
1068  real(dp) :: r_val, d_dummy
1069  character :: ch_dummy
1070  character(len=256) :: filename_with_path
1071 
1072  n_data_kei_max = 10000
1073 
1074  kei = 0.0_dp
1075 
1076 !-------- Reading of tabulated values --------
1077 
1078  filename_with_path = trim(inpath)//'/general/kei.dat'
1079 
1080 #if !defined(ALLOW_OPENAD) /* Normal */
1081  open(unit=11, iostat=ios, file=trim(filename_with_path), status='old')
1082 #else /* OpenAD */
1083  open(unit=11, iostat=ios, file=trim(filename_with_path))
1084 #endif /* Normal vs. OpenAD */
1085 
1086  if (ios /= 0) then
1087  errormsg = ' >>> read_kei: Error when opening the kei file!'
1088  call error(errormsg)
1089  end if
1090 
1091  read(unit=11, fmt='(a)') ch_dummy
1092  read(unit=11, fmt='(15x,f5.0)') kei_r_max
1093  read(unit=11, fmt='(15x,f5.0)') kei_r_incr
1094  read(unit=11, fmt='(a)') ch_dummy
1095 
1097 
1098  if (n_data_kei > n_data_kei_max) then
1099  errormsg = ' >>> read_kei: Array kei too small!'
1100  call error(errormsg)
1101  end if
1102 
1103  read(unit=11, fmt='(a)') ch_dummy
1104  read(unit=11, fmt='(a)') ch_dummy
1105  read(unit=11, fmt='(a)') ch_dummy
1106 
1107  do n=-n_data_kei, n_data_kei
1108  read(unit=11, fmt=*) d_dummy, kei(n)
1109  end do
1110 
1111  close(unit=11, status='keep')
1112 
1113  end subroutine read_kei
1114 
1115 !-------------------------------------------------------------------------------
1116 !> Reading of physical parameters.
1117 !<------------------------------------------------------------------------------
1118  subroutine read_phys_para()
1120  use sico_variables_m
1121  use sico_vars_m
1122 
1123  implicit none
1124 
1125  integer(i4b), parameter :: n_unit=31
1126  integer(i4b) :: ios
1127  integer(i4b) :: n
1128  character(len=256) :: filename_with_path
1129 
1130  filename_with_path = trim(inpath)//'/'//trim(ch_domain_short)//'/'// &
1131  trim(phys_para_file)
1132 
1133 #if !defined(ALLOW_OPENAD) /* Normal */
1134  open(n_unit, iostat=ios, file=trim(filename_with_path), status='old')
1135 #else /* OpenAD */
1136  open(n_unit, iostat=ios, file=trim(filename_with_path))
1137 #endif /* Normal vs. OpenAD */
1138 
1139  if (ios /= 0) then
1140  errormsg = ' >>> read_phys_para: Error when opening the phys_para file!'
1141  call error(errormsg)
1142  end if
1143 
1144 #if (!defined(NMARS) && !defined(SMARS)) /* not Martian polar caps */
1145  call read_phys_para_value(rho, n_unit)
1146 #else /* Martian polar caps */
1147  call read_phys_para_value(rho_i, n_unit)
1148 #endif
1149  call read_phys_para_value(rho_w, n_unit)
1150  call read_phys_para_value(rho_sw, n_unit)
1151  call read_phys_para_value(l, n_unit)
1152  call read_phys_para_value(g, n_unit)
1153  call read_phys_para_value(nue, n_unit)
1154  call read_phys_para_value(beta, n_unit)
1155 #if (!defined(NMARS) && !defined(SMARS)) /* not Martian polar caps */
1156  call read_phys_para_value(delta_tm_sw, n_unit)
1157 #endif
1158  call read_phys_para_value(omega_max, n_unit)
1159 #if (defined(NMARS) || defined(SMARS)) /* Martian polar caps */
1160  call read_phys_para_value(rho_c, n_unit)
1161  call read_phys_para_value(kappa_c, n_unit)
1162  call read_phys_para_value(c_c, n_unit)
1163 #endif
1164  call read_phys_para_value(h_r, n_unit)
1165  call read_phys_para_value(rho_c_r, n_unit)
1166  call read_phys_para_value(kappa_r, n_unit)
1167  call read_phys_para_value(rho_a, n_unit)
1168  call read_phys_para_value(r_t, n_unit)
1169 
1170  call read_phys_para_value(r, n_unit)
1171  call read_phys_para_value(a, n_unit)
1172  call read_phys_para_value(b, n_unit)
1173  call read_phys_para_value(phi0, n_unit)
1174  call read_phys_para_value(lambda0, n_unit)
1175 
1176 #if (!defined(NMARS) && !defined(SMARS)) /* not Martian polar caps */
1177  call read_phys_para_value(s_stat_0, n_unit)
1178 #if (!defined(GRL)) /* not Greenland */
1179  call read_phys_para_value(beta1_0, n_unit)
1180  call read_phys_para_value(beta2_0, n_unit)
1181 #else /* Greenland */
1182  call read_phys_para_value(beta1_lt_0, n_unit)
1183  call read_phys_para_value(beta1_ht_0, n_unit)
1184  call read_phys_para_value(beta2_lt_0, n_unit)
1185  call read_phys_para_value(beta2_ht_0, n_unit)
1186  call read_phys_para_value(phi_sep_0, n_unit)
1187 #endif
1188  call read_phys_para_value(pmax_0, n_unit)
1189  call read_phys_para_value(mu_0, n_unit)
1190 #endif
1191 
1192  do n=10, -190, -1
1193  call read_phys_para_value(rf(n), n_unit)
1194  end do
1195 
1196  do n=10, -190, -1
1197  call read_phys_para_value(kappa(n), n_unit)
1198  end do
1199 
1200  do n=10, -190, -1
1201  call read_phys_para_value(c(n), n_unit)
1202  end do
1203 
1204  close(n_unit, status='keep')
1205 
1206  end subroutine read_phys_para
1207 
1208 !-------------------------------------------------------------------------------
1209 !> Reading of a value of a physical parameter from the phys_para file.
1210 !<------------------------------------------------------------------------------
1211  subroutine read_phys_para_value(para, n_unit)
1213  implicit none
1214 
1215  integer(i4b), intent(in) :: n_unit
1216  real(dp), intent(out) :: para
1217 
1218  character :: check
1219 
1220 #if !defined(ALLOW_OPENAD) /* Normal */
1221 
1222  integer :: i0
1223  character(len=1024) :: txt
1224 
1225 #else /* OpenAD */
1226 
1227  character :: ch_dummy
1228  character(len=*), parameter :: fmt1='(a)', fmt3='(15x)'
1229  logical :: first_read
1230 
1231  first_read = .true.
1232 
1233 #endif /* Normal vs. OpenAD */
1234 
1235  check = '%'
1236 
1237  do while (check == '%')
1238 
1239 #if !defined(ALLOW_OPENAD) /* Normal */
1240 
1241  read(unit=n_unit, fmt='(a)') txt
1242 
1243  check = txt(1:1)
1244 
1245  if (check /= '%') then ! no comment line
1246  i0 = index(txt, '=') + 1
1247  read(txt(i0:), *) para
1248  end if
1249 
1250 #else /* OpenAD: cannot differentiate through index function */
1251 
1252  if (first_read) then
1253  read(unit=n_unit, fmt=trim(fmt1), advance='no') check
1254  first_read=.false.
1255  else
1256  read(unit=n_unit, fmt=trim(fmt1)) ch_dummy
1257  read(unit=n_unit, fmt=trim(fmt1), advance='no') check
1258  end if
1259 
1260  if (check /= '%') then ! no comment line
1261  read(unit=n_unit, fmt=trim(fmt3), advance='no')
1262  read(unit=n_unit, fmt=*) para
1263  end if
1264 
1265 #endif /* Normal vs. OpenAD */
1266 
1267  end do
1268 
1269  end subroutine read_phys_para_value
1270 
1271 #if (defined(ALLOW_GRDCHK) || defined(ALLOW_OPENAD))
1272 !-------------------------------------------------------------------------------
1273 !> Reading in of adjoint related data (only set up for ages right now).
1274 !<------------------------------------------------------------------------------
1275  subroutine read_ad_data()
1276 
1277  use sico_variables_m
1278  use sico_vars_m
1279 
1280  implicit none
1281 
1282  integer(i4b) :: ios, i, j, k, kc, kt, KDATA
1283  logical :: debug = .true.
1284 
1285 #if (defined(AGE_COST)) /* Only designed for AGE_COST right now */
1286 ! Warning: code below is not robust against
1287 ! CALCMOD.NE.3,2
1288 #if (CALCMOD==3)
1289  kdata = kcmax
1290 #else
1291  kdata = kcmax + ktmax
1292 #endif
1293 
1294  ! Note: the last lines of these files
1295  ! correspond to the surface of the ice sheet:
1296  open(unit=90, iostat=ios, &
1297 file='subroutines/openad/gridded_age_obs.dat', status='old')
1298  open(unit=91, iostat=ios, &
1299 file='subroutines/openad/gridded_age_unc.dat', status='old')
1300  do k=0,kdata
1301  do j=0,jmax
1302  do i=0,imax
1303 !write(6, fmt='(a,i3,a,i3,a,i3,a)', advance="no") '(', k, ',', j, ',', i, ') '
1304  if (i.lt.imax) then
1305  read(unit=90, fmt='(es11.4)', advance="no") age_data(k,j,i)
1306  read(unit=91, fmt='(es11.4)', advance="no") age_unc(k,j,i)
1307  else ! do a carriage return
1308  read(unit=90, fmt='(es11.4)' ) age_data(k,j,i)
1309  read(unit=91, fmt='(es11.4)' ) age_unc(k,j,i)
1310  end if
1311  end do
1312 !write(6, fmt='(a)') ' '
1313  end do
1314  end do
1315  close(unit=90)
1316  close(unit=91)
1317 
1318  ! ages yr --> s
1319  do k=0,kdata
1320  do j=0,jmax
1321  do i=0,imax
1322  age_data(k,j,i) = age_data(k,j,i) * year_sec
1323  age_unc(k,j,i) = age_unc(k,j,i) * year_sec
1324  if (debug) then ! cheat the FDS by making initial ages = to data
1325  age_c (k,j,i) = age_data(k,j,i)
1326  age_c_neu(k,j,i) = age_data(k,j,i)
1327  end if
1328  end do
1329  end do
1330  end do
1331 
1332 #endif /* Only designed for AGE_COST right now */
1333 
1334  end subroutine read_ad_data
1335 
1336 #endif /* inclusion of OpenAD only routine */
1337 
1338 !-------------------------------------------------------------------------------
1339 
1340 end module read_m
1341 !
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.
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) 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:1061
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) 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:1119
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.
logical, dimension(0:jmax, 0:imax) flag_grounding_line_1
flag_grounding_line_1(j,i): Grounding line flag. .true.: grounding line point (grounded ice point wit...
real(dp), dimension(0:jmax, 0:imax) vis_ave_g
vis_ave_g(j,i): Depth-averaged viscosity of the SIA/SSA, at (i,j)
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)
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) p_b_w
p_b_w(j,i): Basal water pressure
integer(i1b), 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) vx_m_ssa
vx_m_ssa(j,i): Mean (depth-averaged) SSA velocity in x-direction, at (i+1/2,j)
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)
logical, dimension(0:jmax, 0:imax) flag_shelfy_stream_x
flag_shelfy_stream_x(j,i): Shelfy stream flag in x-direction, at (i+1/2,j). .true.: shelfy stream point .false.: otherwise
logical, dimension(0:jmax, 0:imax) flag_grounded_front_b_1
flag_grounded_front_b_1(j,i): Marine-terminating grounded front flag. .true.: grounded front point (g...
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(i1b), dimension(0:jmax, 0:imax) maske_target
maske_target(j,i): Target topography (ice-land-ocean mask)
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
real(dp), dimension(0:jmax, 0:imax) vis_int_g
vis_int_g(j,i): Depth-integrated viscosity of the SIA/SSA, at (i,j)
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:1212
real(dp), dimension(0:jmax, 0:imax) ratio_sl_y
ratio_sl_y(j,i): Ratio of basal to surface velocity (slip ratio) in y-direction, at (i...
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:jmax, 0:imax) vy_m_ssa
vy_m_ssa(j,i): Mean (depth-averaged) SSA velocity in y-direction, at (i,j+1/2)
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), dimension(0:kcmax, 0:jmax, 0:imax) age_c_neu
(.)_neu: New value of quantity (.) computed during an integration step
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.
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: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
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), parameter eps_dp
eps_dp: Small number to double-precision accuracy
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.
logical, dimension(0:jmax, 0:imax) flag_shelfy_stream
flag_shelfy_stream(j,i): Shelfy stream flag on the main grid. .true.: grounded ice, and at least one neighbour on the staggered grid is a shelfy stream point .false.: otherwise
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)
real(dp), dimension(0:jmax, 0:imax) vy_m_sia
vy_m_sia(j,i): Mean (depth-averaged) SIA velocity in y-direction, at (i,j+1/2)
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.
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:jmax, 0:imax) q_bm
Q_bm(j,i): Basal melting rate.
logical, dimension(0:jmax, 0:imax) flag_grounded_front_b_2
flag_grounded_front_b_2(j,i): Marine-terminating grounded front flag. .true.: grounded front point (o...
real(dp), dimension(0:jmax, 0:imax) q_w
q_w(j,i): Scalar volume flux of the basal water
real(dp), dimension(-190:10) c
C(n): Tabulated values for the specific heat of ice.
real(dp), dimension(0:jmax, 0:imax) q_w_y
q_w_y(j,i): Scalar volume flux of the basal water in y-direction
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
logical, dimension(0:jmax, 0:imax) flag_grounded_front_a_2
flag_grounded_front_a_2(j,i): Land-terminating grounded front flag. .true.: grounded front point (ice...
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.
logical, dimension(0:jmax, 0:imax) flag_calving_front_1
flag_calving_front_1(j,i): Calving front flag. .true.: calving front point (floating ice point with a...
real(dp), dimension(0:jmax, 0:imax) ratio_sl_x
ratio_sl_x(j,i): Ratio of basal to surface velocity (slip ratio) in x-direction, at (i+1/2...
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.
logical, dimension(0:jmax, 0:imax) flag_calving_front_2
flag_calving_front_2(j,i): Calving front flag. .true.: calving front point (ocean point with at least...
logical, dimension(0:jmax, 0:imax) flag_shelfy_stream_y
flag_shelfy_stream_y(j,i): Shelfy stream flag in y-direction, at (i,j+1/2). .true.: shelfy stream point .false.: otherwise
real(dp) l
L: Latent heat of ice.
integer(i4b) n_data_kei
n_data_kei: Number of tabulated values of the kei function
real(dp), dimension(0:jmax, 0:imax) q_w_x
q_w_x(j,i): Scalar volume flux of the basal water in x-direction
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
logical, dimension(0:jmax, 0:imax) flag_grounding_line_2
flag_grounding_line_2(j,i): Grounding line flag. .true.: grounding line point (floating ice point wit...
real(dp), dimension(0:jmax, 0:imax) vx_m_sia
vx_m_sia(j,i): Mean (depth-averaged) SIA velocity in x-direction, at (i+1/2,j)
real(dp), dimension(0:jmax, 0:imax) dzb_dtau
dzb_dtau(j,i): Derivative of zb by tau (time)
subroutine, public read_target_topo_nc(target_topo_dat_name)
Reading of the target-topography file (in NetCDF format).
Definition: read_m.F90:950
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.
logical, dimension(0:jmax, 0:imax) flag_grounded_front_a_1
flag_grounded_front_a_1(j,i): Land-terminating grounded front flag. .true.: grounded front point (gro...
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))