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