54 #if (NETCDF==2) /* time-slice file in NetCDF format */ 61 character(len=100),
intent(in) :: filename
63 real(dp),
intent(out) :: z_sl
65 integer(i4b) :: i, j, kc, kt, kr, n
66 character(len=256) :: filename_with_path
68 #if (NETCDF==1) /* time-slice file in native binary format */ 71 character(len=256) :: ch_attr_title, ch_attr_institution, ch_attr_source, &
72 ch_attr_history, ch_attr_references
74 #elif (NETCDF==2) /* time-slice file in NetCDF format */ 76 integer(i4b) :: ncid, ncv
81 stop
' >>> read_erg_nc: Parameter NETCDF must be either 1 or 2!' 84 integer(i2b),
dimension(0:IMAX,0:JMAX) :: maske_conv, maske_old_conv, &
85 n_cts_conv, kc_cts_conv
86 real(sp) :: time_conv, dummy_conv, z_sl_conv, &
87 V_tot_conv, A_grounded_conv, A_floating_conv, &
89 xi_conv(0:imax), eta_conv(0:jmax), &
90 sigma_level_c_conv(0:kcmax), sigma_level_t_conv(0:ktmax), &
91 sigma_level_r_conv(0:krmax)
92 real(sp),
dimension(0:IMAX,0:JMAX) :: lambda_conv, phi_conv, &
94 temp_s_conv, accum_conv, as_perp_conv, as_perp_apl_conv, &
96 zs_conv, zm_conv, zb_conv, zl_conv, zl0_conv, &
97 H_cold_conv, H_temp_conv, H_conv, &
98 Q_bm_conv, Q_tld_conv, &
101 dzs_dtau_conv, dzm_dtau_conv, dzb_dtau_conv, dzl_dtau_conv, &
102 dH_c_dtau_conv, dH_t_dtau_conv, dH_dtau_conv, &
103 vx_b_g_conv, vy_b_g_conv, vz_b_conv, vh_b_conv, &
104 vx_s_g_conv, vy_s_g_conv, vz_s_conv, vh_s_conv, &
105 vx_m_g_conv, vy_m_g_conv, vh_m_conv, &
106 temp_b_conv, temph_b_conv, &
107 tau_b_driving_conv, tau_b_drag_conv, &
108 p_b_w_conv, H_w_conv, q_gl_g_conv, q_cf_g_conv
109 real(sp),
dimension(0:IMAX,0:JMAX,0:KCMAX) :: vx_c_conv, vy_c_conv, vz_c_conv, &
110 temp_c_conv, age_c_conv, &
111 enth_c_conv, omega_c_conv, &
113 real(sp),
dimension(0:IMAX,0:JMAX,0:KTMAX) :: vx_t_conv, vy_t_conv, vz_t_conv, &
114 omega_t_conv, age_t_conv, &
117 real(sp),
dimension(0:IMAX,0:JMAX,0:KRMAX) :: temp_r_conv
119 #if (DISC>0) /* Ice discharge parameterisation */ 120 integer(i2b),
dimension(0:IMAX,0:JMAX) :: mask_mar_conv
121 real(sp),
dimension(0:IMAX,0:JMAX) :: dis_perp_conv, &
122 cst_dist_conv, cos_grad_tc_conv
127 filename_with_path = trim(anfdatpath)//
'/'//trim(filename)
129 #if (NETCDF==1) /* time-slice file in native binary format */ 131 open(unit=11, iostat=ios, file=trim(filename_with_path), status=
'old', &
134 if (ios /= 0) stop
' >>> read_erg_nc: Error when opening the time-slice file!' 136 read(unit=11) ch_attr_title
137 read(unit=11) ch_attr_institution
138 read(unit=11) ch_attr_source
139 read(unit=11) ch_attr_history
140 read(unit=11) ch_attr_references
142 read(unit=11) time_conv
143 read(unit=11) dummy_conv
144 read(unit=11) z_sl_conv
146 read(unit=11) dummy_conv
147 read(unit=11) dummy_conv
148 read(unit=11) dummy_conv
149 read(unit=11) dummy_conv
151 read(unit=11) xi_conv
152 read(unit=11) eta_conv
153 read(unit=11) sigma_level_c_conv
154 read(unit=11) sigma_level_t_conv
155 read(unit=11) sigma_level_r_conv
156 read(unit=11) lon_conv
157 read(unit=11) lat_conv
158 read(unit=11) lambda_conv
159 read(unit=11) phi_conv
160 read(unit=11) temp_s_conv
161 read(unit=11) accum_conv
162 read(unit=11) as_perp_conv
163 read(unit=11) as_perp_apl_conv
165 #if (DISC>0) /* Ice discharge parameterisation */ 167 read(unit=11) dis_perp_conv
168 read(unit=11) cst_dist_conv
169 read(unit=11) cos_grad_tc_conv
170 read(unit=11) mask_mar_conv
174 read(unit=11) q_geo_conv
175 read(unit=11) maske_conv
176 read(unit=11) maske_old_conv
177 read(unit=11) n_cts_conv
178 read(unit=11) kc_cts_conv
179 read(unit=11) zs_conv
180 read(unit=11) zm_conv
181 read(unit=11) zb_conv
182 read(unit=11) zl_conv
183 read(unit=11) zl0_conv
184 read(unit=11) h_cold_conv
185 read(unit=11) h_temp_conv
187 read(unit=11) h_r_conv
188 read(unit=11) q_bm_conv
189 read(unit=11) q_tld_conv
190 read(unit=11) am_perp_conv
191 read(unit=11) qx_conv
192 read(unit=11) qy_conv
193 read(unit=11) dzs_dtau_conv
194 read(unit=11) dzm_dtau_conv
195 read(unit=11) dzb_dtau_conv
196 read(unit=11) dzl_dtau_conv
197 read(unit=11) dh_c_dtau_conv
198 read(unit=11) dh_t_dtau_conv
199 read(unit=11) dh_dtau_conv
200 read(unit=11) vx_b_g_conv
201 read(unit=11) vy_b_g_conv
202 read(unit=11) vz_b_conv
203 read(unit=11) vh_b_conv
204 read(unit=11) vx_s_g_conv
205 read(unit=11) vy_s_g_conv
206 read(unit=11) vz_s_conv
207 read(unit=11) vh_s_conv
208 read(unit=11) vx_m_g_conv
209 read(unit=11) vy_m_g_conv
210 read(unit=11) vh_m_conv
211 read(unit=11) temp_b_conv
212 read(unit=11) temph_b_conv
213 read(unit=11) tau_b_driving_conv
214 read(unit=11) tau_b_drag_conv
215 read(unit=11) p_b_w_conv
216 read(unit=11) h_w_conv
217 read(unit=11) q_gl_g_conv
218 read(unit=11) q_cf_g_conv
220 read(unit=11) vx_c_conv
221 read(unit=11) vy_c_conv
222 read(unit=11) vz_c_conv
223 read(unit=11) vx_t_conv
224 read(unit=11) vy_t_conv
225 read(unit=11) vz_t_conv
226 read(unit=11) temp_c_conv
227 read(unit=11) omega_t_conv
228 read(unit=11) temp_r_conv
229 read(unit=11) enth_c_conv
230 read(unit=11) enth_t_conv
231 read(unit=11) omega_c_conv
232 read(unit=11) enh_c_conv
233 read(unit=11) enh_t_conv
234 read(unit=11) age_c_conv
235 read(unit=11) age_t_conv
237 close(unit=11,status=
'keep')
239 #elif (NETCDF==2) /* time-slice file in NetCDF format */ 241 call check( nf90_open(trim(filename_with_path), nf90_nowrite, ncid) )
243 call check( nf90_inq_varid(ncid,
'time', ncv) )
244 call check( nf90_get_var(ncid, ncv, time_conv) )
246 if ( nf90_inq_varid(ncid,
'delta_ts', ncv) == nf90_noerr )
then 247 call check( nf90_get_var(ncid, ncv, dummy_conv) )
248 else if ( nf90_inq_varid(ncid,
'glac_index', ncv) == nf90_noerr )
then 249 call check( nf90_get_var(ncid, ncv, dummy_conv) )
254 call check( nf90_inq_varid(ncid,
'z_sl', ncv) )
255 call check( nf90_get_var(ncid, ncv, z_sl_conv) )
257 call check( nf90_inq_varid(ncid,
'x', ncv) )
258 call check( nf90_get_var(ncid, ncv, xi_conv) )
260 call check( nf90_inq_varid(ncid,
'y', ncv) )
261 call check( nf90_get_var(ncid, ncv, eta_conv) )
263 call check( nf90_inq_varid(ncid,
'sigma_level_c', ncv) )
264 call check( nf90_get_var(ncid, ncv, sigma_level_c_conv) )
266 call check( nf90_inq_varid(ncid,
'sigma_level_t', ncv) )
267 call check( nf90_get_var(ncid, ncv, sigma_level_t_conv) )
269 call check( nf90_inq_varid(ncid,
'sigma_level_r', ncv) )
270 call check( nf90_get_var(ncid, ncv, sigma_level_r_conv) )
272 call check( nf90_inq_varid(ncid,
'lon', ncv) )
273 call check( nf90_get_var(ncid, ncv, lon_conv) )
275 call check( nf90_inq_varid(ncid,
'lat', ncv) )
276 call check( nf90_get_var(ncid, ncv, lat_conv) )
278 call check( nf90_inq_varid(ncid,
'lambda', ncv) )
279 call check( nf90_get_var(ncid, ncv, lambda_conv) )
281 call check( nf90_inq_varid(ncid,
'phi', ncv) )
282 call check( nf90_get_var(ncid, ncv, phi_conv) )
284 call check( nf90_inq_varid(ncid,
'temp_s', ncv) )
285 call check( nf90_get_var(ncid, ncv, temp_s_conv) )
287 if ( nf90_inq_varid(ncid,
'prec', ncv) == nf90_noerr )
then 288 call check( nf90_get_var(ncid, ncv, accum_conv) )
290 write(6,
'(/1x,a)')
'>>> read_erg_nc: Variable prec' 291 write(6,
'(1x,a)')
' not available in read file *.nc.' 295 call check( nf90_inq_varid(ncid,
'as_perp', ncv) )
296 call check( nf90_get_var(ncid, ncv, as_perp_conv) )
298 if ( nf90_inq_varid(ncid,
'as_perp_apl', ncv) == nf90_noerr )
then 299 call check( nf90_get_var(ncid, ncv, as_perp_apl_conv) )
301 write(6,
'(/1x,a)')
'>>> read_erg_nc: Variable as_perp_apl' 302 write(6,
'(1x,a)')
' not available in read file *.nc.' 303 as_perp_apl_conv = 0.0_sp
306 #if (DISC>0) /* Ice discharge parameterisation */ 308 if ( nf90_inq_varid(ncid,
'dis_perp', ncv) == nf90_noerr )
then 309 call check( nf90_get_var(ncid, ncv, dis_perp_conv) )
311 write(6,
'(/1x,a)')
'>>> read_erg_nc: Variable dis_perp' 312 write(6,
'(1x,a)')
' not available in read file *.nc.' 313 dis_perp_conv = 0.0_sp
316 if ( nf90_inq_varid(ncid,
'cst_dist', ncv) == nf90_noerr )
then 317 call check( nf90_get_var(ncid, ncv, cst_dist_conv) )
319 write(6,
'(/1x,a)')
'>>> read_erg_nc: Variable cst_dist' 320 write(6,
'(1x,a)')
' not available in read file *.nc.' 324 if ( nf90_inq_varid(ncid,
'cos_grad_tc', ncv) == nf90_noerr )
then 325 call check( nf90_get_var(ncid, ncv, cos_grad_tc_conv) )
327 write(6,
'(/1x,a)')
'>>> read_erg_nc: Variable cos_grad_tc' 328 write(6,
'(1x,a)')
' not available in read file *.nc.' 329 cos_grad_tc_conv = 1.0_sp
332 if ( nf90_inq_varid(ncid,
'mask_mar', ncv) == nf90_noerr )
then 333 call check( nf90_get_var(ncid, ncv, mask_mar_conv) )
335 write(6,
'(/1x,a)')
'>>> read_erg_nc: Variable mask_mar' 336 write(6,
'(1x,a)')
' not available in read file *.nc.' 337 mask_mar_conv = 0_i2b
342 if ( nf90_inq_varid(ncid,
'q_geo', ncv) == nf90_noerr )
then 343 call check( nf90_get_var(ncid, ncv, q_geo_conv) )
345 write(6,
'(/1x,a)')
'>>> read_erg_nc: Variable q_geo' 346 write(6,
'(1x,a)')
' not available in read file *.nc.' 350 call check( nf90_inq_varid(ncid,
'maske', ncv) )
351 call check( nf90_get_var(ncid, ncv, maske_conv) )
353 if ( nf90_inq_varid(ncid,
'maske_old', ncv) == nf90_noerr )
then 354 call check( nf90_get_var(ncid, ncv, maske_old_conv) )
356 write(6,
'(/1x,a)')
'>>> read_erg_nc: Variable maske_old' 357 write(6,
'(1x,a)')
' not available in read file *.nc.' 358 maske_old_conv = maske_conv
361 call check( nf90_inq_varid(ncid,
'n_cts', ncv) )
362 call check( nf90_get_var(ncid, ncv, n_cts_conv) )
364 call check( nf90_inq_varid(ncid,
'kc_cts', ncv) )
365 call check( nf90_get_var(ncid, ncv, kc_cts_conv) )
367 call check( nf90_inq_varid(ncid,
'zs', ncv) )
368 call check( nf90_get_var(ncid, ncv, zs_conv) )
370 call check( nf90_inq_varid(ncid,
'zm', ncv) )
371 call check( nf90_get_var(ncid, ncv, zm_conv) )
373 call check( nf90_inq_varid(ncid,
'zb', ncv) )
374 call check( nf90_get_var(ncid, ncv, zb_conv) )
376 call check( nf90_inq_varid(ncid,
'zl', ncv) )
377 call check( nf90_get_var(ncid, ncv, zl_conv) )
379 if ( nf90_inq_varid(ncid,
'zl0', ncv) == nf90_noerr )
then 380 call check( nf90_get_var(ncid, ncv, zl0_conv) )
382 write(6,
'(/1x,a)')
'>>> read_erg_nc: Variable zl0' 383 write(6,
'(1x,a)')
' not available in read file *.nc.' 387 call check( nf90_inq_varid(ncid,
'H_cold', ncv) )
388 call check( nf90_get_var(ncid, ncv, h_cold_conv) )
390 call check( nf90_inq_varid(ncid,
'H_temp', ncv) )
391 call check( nf90_get_var(ncid, ncv, h_temp_conv) )
393 call check( nf90_inq_varid(ncid,
'H', ncv) )
394 call check( nf90_get_var(ncid, ncv, h_conv) )
396 call check( nf90_inq_varid(ncid,
'H_R', ncv) )
397 call check( nf90_get_var(ncid, ncv, h_r_conv) )
399 call check( nf90_inq_varid(ncid,
'Q_bm', ncv) )
400 call check( nf90_get_var(ncid, ncv, q_bm_conv) )
402 call check( nf90_inq_varid(ncid,
'Q_tld', ncv) )
403 call check( nf90_get_var(ncid, ncv, q_tld_conv) )
405 call check( nf90_inq_varid(ncid,
'am_perp', ncv) )
406 call check( nf90_get_var(ncid, ncv, am_perp_conv) )
408 call check( nf90_inq_varid(ncid,
'qx', ncv) )
409 call check( nf90_get_var(ncid, ncv, qx_conv) )
411 call check( nf90_inq_varid(ncid,
'qy', ncv) )
412 call check( nf90_get_var(ncid, ncv, qy_conv) )
414 call check( nf90_inq_varid(ncid,
'dzs_dt', ncv) )
415 call check( nf90_get_var(ncid, ncv, dzs_dtau_conv) )
417 call check( nf90_inq_varid(ncid,
'dzm_dt', ncv) )
418 call check( nf90_get_var(ncid, ncv, dzm_dtau_conv) )
420 call check( nf90_inq_varid(ncid,
'dzb_dt', ncv) )
421 call check( nf90_get_var(ncid, ncv, dzb_dtau_conv) )
423 call check( nf90_inq_varid(ncid,
'dzl_dt', ncv) )
424 call check( nf90_get_var(ncid, ncv, dzl_dtau_conv) )
426 call check( nf90_inq_varid(ncid,
'dH_c_dt', ncv) )
427 call check( nf90_get_var(ncid, ncv, dh_c_dtau_conv) )
429 call check( nf90_inq_varid(ncid,
'dH_t_dt', ncv) )
430 call check( nf90_get_var(ncid, ncv, dh_t_dtau_conv) )
432 call check( nf90_inq_varid(ncid,
'dH_dt', ncv) )
433 call check( nf90_get_var(ncid, ncv, dh_dtau_conv) )
435 call check( nf90_inq_varid(ncid,
'vx_b_g', ncv) )
436 call check( nf90_get_var(ncid, ncv, vx_b_g_conv) )
438 call check( nf90_inq_varid(ncid,
'vy_b_g', ncv) )
439 call check( nf90_get_var(ncid, ncv, vy_b_g_conv) )
441 call check( nf90_inq_varid(ncid,
'vz_b', ncv) )
442 call check( nf90_get_var(ncid, ncv, vz_b_conv) )
444 call check( nf90_inq_varid(ncid,
'vh_b', ncv) )
445 call check( nf90_get_var(ncid, ncv, vh_b_conv) )
447 call check( nf90_inq_varid(ncid,
'vx_s_g', ncv) )
448 call check( nf90_get_var(ncid, ncv, vx_s_g_conv) )
450 call check( nf90_inq_varid(ncid,
'vy_s_g', ncv) )
451 call check( nf90_get_var(ncid, ncv, vy_s_g_conv) )
453 call check( nf90_inq_varid(ncid,
'vz_s', ncv) )
454 call check( nf90_get_var(ncid, ncv, vz_s_conv) )
456 call check( nf90_inq_varid(ncid,
'vh_s', ncv) )
457 call check( nf90_get_var(ncid, ncv, vh_s_conv) )
459 if ( nf90_inq_varid(ncid,
'vx_m_g', ncv) == nf90_noerr )
then 460 call check( nf90_get_var(ncid, ncv, vx_m_g_conv) )
462 write(6,
'(/1x,a)')
'>>> read_erg_nc: Variable vx_m_g' 463 write(6,
'(1x,a)')
' not available in read file *.nc.' 467 if ( nf90_inq_varid(ncid,
'vy_m_g', ncv) == nf90_noerr )
then 468 call check( nf90_get_var(ncid, ncv, vy_m_g_conv) )
470 write(6,
'(/1x,a)')
'>>> read_erg_nc: Variable vy_m_g' 471 write(6,
'(1x,a)')
' not available in read file *.nc.' 475 if ( nf90_inq_varid(ncid,
'vh_m', ncv) == nf90_noerr )
then 476 call check( nf90_get_var(ncid, ncv, vh_m_conv) )
478 write(6,
'(/1x,a)')
'>>> read_erg_nc: Variable vh_m' 479 write(6,
'(1x,a)')
' not available in read file *.nc.' 483 call check( nf90_inq_varid(ncid,
'temp_b', ncv) )
484 call check( nf90_get_var(ncid, ncv, temp_b_conv) )
486 call check( nf90_inq_varid(ncid,
'temph_b', ncv) )
487 call check( nf90_get_var(ncid, ncv, temph_b_conv) )
489 if ( nf90_inq_varid(ncid,
'tau_b_driving', ncv) == nf90_noerr )
then 490 call check( nf90_get_var(ncid, ncv, tau_b_driving_conv) )
492 write(6,
'(/1x,a)')
'>>> read_erg_nc: Variable tau_b_driving' 493 write(6,
'(1x,a)')
' not available in read file *.nc.' 494 tau_b_driving_conv = 0.0_sp
497 if ( nf90_inq_varid(ncid,
'tau_b_drag', ncv) == nf90_noerr )
then 498 call check( nf90_get_var(ncid, ncv, tau_b_drag_conv) )
500 write(6,
'(/1x,a)')
'>>> read_erg_nc: Variable tau_b_drag' 501 write(6,
'(1x,a)')
' not available in read file *.nc.' 502 tau_b_drag_conv = 0.0_sp
505 call check( nf90_inq_varid(ncid,
'p_b_w', ncv) )
506 call check( nf90_get_var(ncid, ncv, p_b_w_conv) )
508 call check( nf90_inq_varid(ncid,
'H_w', ncv) )
509 call check( nf90_get_var(ncid, ncv, h_w_conv) )
511 call check( nf90_inq_varid(ncid,
'q_gl_g', ncv) )
512 call check( nf90_get_var(ncid, ncv, q_gl_g_conv) )
514 if ( nf90_inq_varid(ncid,
'q_cf_g', ncv) == nf90_noerr )
then 515 call check( nf90_get_var(ncid, ncv, q_cf_g_conv) )
517 write(6,
'(/1x,a)')
'>>> read_erg_nc: Variable q_cf_g' 518 write(6,
'(1x,a)')
' not available in read file *.nc.' 522 call check( nf90_inq_varid(ncid,
'vx_c', ncv) )
523 call check( nf90_get_var(ncid, ncv, vx_c_conv) )
525 call check( nf90_inq_varid(ncid,
'vy_c', ncv) )
526 call check( nf90_get_var(ncid, ncv, vy_c_conv) )
528 call check( nf90_inq_varid(ncid,
'vz_c', ncv) )
529 call check( nf90_get_var(ncid, ncv, vz_c_conv) )
531 call check( nf90_inq_varid(ncid,
'vx_t', ncv) )
532 call check( nf90_get_var(ncid, ncv, vx_t_conv) )
534 call check( nf90_inq_varid(ncid,
'vy_t', ncv) )
535 call check( nf90_get_var(ncid, ncv, vy_t_conv) )
537 call check( nf90_inq_varid(ncid,
'vz_t', ncv) )
538 call check( nf90_get_var(ncid, ncv, vz_t_conv) )
540 call check( nf90_inq_varid(ncid,
'temp_c', ncv) )
541 call check( nf90_get_var(ncid, ncv, temp_c_conv) )
543 call check( nf90_inq_varid(ncid,
'omega_t', ncv) )
544 call check( nf90_get_var(ncid, ncv, omega_t_conv) )
546 call check( nf90_inq_varid(ncid,
'temp_r', ncv) )
547 call check( nf90_get_var(ncid, ncv, temp_r_conv) )
549 call check( nf90_inq_varid(ncid,
'enth_c', ncv) )
550 call check( nf90_get_var(ncid, ncv, enth_c_conv) )
552 call check( nf90_inq_varid(ncid,
'enth_t', ncv) )
553 call check( nf90_get_var(ncid, ncv, enth_t_conv) )
555 call check( nf90_inq_varid(ncid,
'omega_c', ncv) )
556 call check( nf90_get_var(ncid, ncv, omega_c_conv) )
558 if ( nf90_inq_varid(ncid,
'enh_c', ncv) == nf90_noerr )
then 559 call check( nf90_get_var(ncid, ncv, enh_c_conv) )
561 write(6,
'(/1x,a)')
'>>> read_erg_nc: Variable enh_c' 562 write(6,
'(1x,a)')
' not available in read file *.nc.' 566 if ( nf90_inq_varid(ncid,
'enh_t', ncv) == nf90_noerr )
then 567 call check( nf90_get_var(ncid, ncv, enh_t_conv) )
569 write(6,
'(/1x,a)')
'>>> read_erg_nc: Variable enh_t' 570 write(6,
'(1x,a)')
' not available in read file *.nc.' 574 call check( nf90_inq_varid(ncid,
'age_c', ncv) )
575 call check( nf90_get_var(ncid, ncv, age_c_conv) )
577 call check( nf90_inq_varid(ncid,
'age_t', ncv) )
578 call check( nf90_get_var(ncid, ncv, age_t_conv) )
580 call check( nf90_close(ncid) )
586 z_sl =
real(z_sl_conv,
dp)
588 h_r =
real(h_r_conv,
dp)
591 xi(i) =
real(xi_conv(i),dp)
595 eta(j) =
real(eta_conv(j),dp)
601 maske(j,i) = maske_conv(i,j)
603 n_cts(j,i) = n_cts_conv(i,j)
604 kc_cts(j,i) = kc_cts_conv(i,j)
605 zs(j,i) =
real(zs_conv(i,j),dp)
606 zm(j,i) =
real(zm_conv(i,j),dp)
607 zb(j,i) =
real(zb_conv(i,j),dp)
608 zl(j,i) =
real(zl_conv(i,j),dp)
610 h_c(j,i) =
real(H_cold_conv(i,j),dp)
611 h_t(j,i) =
real(H_temp_conv(i,j),dp)
612 #elif (CALCMOD==0 || CALCMOD==2 || CALCMOD==3 || CALCMOD==-1) 613 h_c(j,i) =
real(H_conv(i,j),dp)
616 q_bm(j,i) =
real(Q_bm_conv(i,j),dp)/year_sec
617 q_tld(j,i) =
real(Q_tld_conv(i,j),dp)/year_sec
618 am_perp(j,i) =
real(am_perp_conv(i,j),dp)/year_sec
619 qx(j,i) =
real(qx_conv(i,j),dp)/year_sec
620 qy(j,i) =
real(qy_conv(i,j),dp)/year_sec
621 dzs_dtau(j,i) =
real(dzs_dtau_conv(i,j),dp)/year_sec
622 dzm_dtau(j,i) =
real(dzm_dtau_conv(i,j),dp)/year_sec
623 dzb_dtau(j,i) =
real(dzb_dtau_conv(i,j),dp)/year_sec
624 dzl_dtau(j,i) =
real(dzl_dtau_conv(i,j),dp)/year_sec
625 dh_c_dtau(j,i) =
real(dH_c_dtau_conv(i,j),dp)/year_sec
626 dh_t_dtau(j,i) =
real(dH_t_dtau_conv(i,j),dp)/year_sec
627 vx_b_g(j,i) =
real(vx_b_g_conv(i,j),dp)/year_sec
628 vy_b_g(j,i) =
real(vy_b_g_conv(i,j),dp)/year_sec
629 vz_b(j,i) =
real(vz_b_conv(i,j),dp)/year_sec
630 vx_s_g(j,i) =
real(vx_s_g_conv(i,j),dp)/year_sec
631 vy_s_g(j,i) =
real(vy_s_g_conv(i,j),dp)/year_sec
632 vz_s(j,i) =
real(vz_s_conv(i,j),dp)/year_sec
633 temp_b(j,i) =
real(temp_b_conv(i,j),dp)
634 temph_b(j,i) =
real(temph_b_conv(i,j),dp)
635 p_b_w(j,i) =
real(p_b_w_conv(i,j),dp)
636 h_w(j,i) =
real(H_w_conv(i,j),dp)
639 temp_r(kr,j,i) =
real(temp_r_conv(i,j,kr),dp)
643 vx_t(kt,j,i) =
real(vx_t_conv(i,j,kt),dp)/year_sec
644 vy_t(kt,j,i) =
real(vy_t_conv(i,j,kt),dp)/year_sec
645 vz_t(kt,j,i) =
real(vz_t_conv(i,j,kt),dp)/year_sec
646 omega_t(kt,j,i) =
real(omega_t_conv(i,j,kt),dp)
647 age_t(kt,j,i) =
real(age_t_conv(i,j,kt),dp)*year_sec
648 enth_t(kt,j,i) =
real(enth_t_conv(i,j,kt),dp)
649 enh_t(kt,j,i) =
real(enh_t_conv(i,j,kt),dp)
653 vx_c(kc,j,i) =
real(vx_c_conv(i,j,kc),dp)/year_sec
654 vy_c(kc,j,i) =
real(vy_c_conv(i,j,kc),dp)/year_sec
655 vz_c(kc,j,i) =
real(vz_c_conv(i,j,kc),dp)/year_sec
656 temp_c(kc,j,i) =
real(temp_c_conv(i,j,kc),dp)
657 age_c(kc,j,i) =
real(age_c_conv(i,j,kc),dp)*year_sec
658 enth_c(kc,j,i) =
real(enth_c_conv(i,j,kc),dp)
659 omega_c(kc,j,i) =
real(omega_c_conv(i,j,kc),dp)
660 enh_c(kc,j,i) =
real(enh_c_conv(i,j,kc),dp)
684 character(len=100),
intent(in) :: target_topo_dat_name
691 integer(i2b),
dimension(0:IMAX,0:JMAX) :: maske_conv
692 real(sp),
dimension(0:IMAX,0:JMAX) :: zs_conv, zb_conv, zl_conv, H_conv
693 character(len=256) :: target_topo_dat_path
694 character(len=256) :: filename_with_path
697 integer(i4b) :: ncid, ncv
702 character(len=64),
parameter :: thisroutine =
'read_target_topo_nc' 709 #if (defined(TARGET_TOPO_DAT_PATH)) 710 target_topo_dat_path = trim(target_topo_dat_path)
712 target_topo_dat_path =
'dummy' 713 stop
' >>> read_target_topo_nc: TARGET_TOPO_DAT_PATH must be defined!' 717 = trim(target_topo_dat_path)//
'/'//trim(target_topo_dat_name)
719 call check( nf90_open(trim(filename_with_path), nf90_nowrite, ncid), &
722 call check( nf90_inq_varid(ncid,
'maske', ncv), thisroutine )
723 call check( nf90_get_var(ncid, ncv, maske_conv), thisroutine )
725 call check( nf90_inq_varid(ncid,
'zs', ncv), thisroutine )
726 call check( nf90_get_var(ncid, ncv, zs_conv), thisroutine )
728 call check( nf90_inq_varid(ncid,
'zb', ncv), thisroutine )
729 call check( nf90_get_var(ncid, ncv, zb_conv), thisroutine )
731 call check( nf90_inq_varid(ncid,
'zl', ncv), thisroutine )
732 call check( nf90_get_var(ncid, ncv, zl_conv), thisroutine )
734 call check( nf90_inq_varid(ncid,
'H', ncv), thisroutine )
735 call check( nf90_get_var(ncid, ncv, h_conv), thisroutine )
737 call check( nf90_close(ncid), thisroutine )
739 #else /* NETCDF != 2, not allowed */ 747 stop
' >>> read_target_topo_nc: Parameter NETCDF must be set to 2!' 757 h_target =
real(transpose(H_conv) ,dp)
770 integer(i4b) :: n, n_data_kei_max
772 real(dp) :: r_val, d_dummy
773 character :: ch_dummy
774 character(len=256) :: filename_with_path
776 n_data_kei_max = 10000
782 filename_with_path = trim(inpath)//
'/general/kei.dat' 784 open(unit=11, iostat=ios, file=trim(filename_with_path), status=
'old')
786 if (ios /= 0) stop
' >>> read_kei: Error when opening the kei file!' 788 read(unit=11, fmt=
'(a)') ch_dummy
789 read(unit=11, fmt=
'(15x,f5.0)')
kei_r_max 791 read(unit=11, fmt=
'(a)') ch_dummy
795 if (
n_data_kei > n_data_kei_max) stop
' >>> read_kei: Array kei too small!' 797 read(unit=11, fmt=
'(a)') ch_dummy
798 read(unit=11, fmt=
'(a)') ch_dummy
799 read(unit=11, fmt=
'(a)') ch_dummy
802 read(unit=11, fmt=*) d_dummy,
kei(n)
805 close(unit=11, status=
'keep')
819 integer(i4b),
parameter :: n_unit=31
822 character(len=256) :: filename_with_path
827 open(n_unit, iostat=ios, file=trim(filename_with_path), status=
'old')
830 stop
' >>> read_phys_para: Error when opening the phys_para file!' 832 #if (!defined(NMARS) && !defined(SMARS)) /* not Martian polar caps */ 834 #else /* Martian polar caps */ 843 #if (!defined(NMARS) && !defined(SMARS)) /* not Martian polar caps */ 847 #if (defined(NMARS) || defined(SMARS)) /* Martian polar caps */ 864 #if (!defined(NMARS) && !defined(SMARS)) /* not Martian polar caps */ 866 #if (!defined(GRL)) /* not Greenland */ 869 #else /* Greenland */ 892 close(n_unit, status=
'keep')
903 integer(i4b),
intent(in) :: n_unit
904 real(dp),
intent(out) :: para
906 character :: check, ch_dummy
907 character(len=*),
parameter :: fmt1=
'(a)', fmt3=
'(15x)' 908 logical :: first_read
914 read(unit=n_unit, fmt=trim(fmt1), advance=
'no') check
917 read(unit=n_unit, fmt=trim(fmt1)) ch_dummy
918 read(unit=n_unit, fmt=trim(fmt1), advance=
'no') check
921 if (check /=
'%')
then 922 read(unit=n_unit, fmt=trim(fmt3), advance=
'no')
923 read(unit=n_unit, fmt=*) para
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) vx_c
vx_c(kc,j,i): Velocity in x-direction in the upper (kc) ice domain (at (i+1/2,j,kc)) ...
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) temp_c
temp_c(kc,j,i): Temperature in the upper (kc) ice domain
real(dp) kappa_r
KAPPA_R: Heat conductivity of the lithosphere.
real(dp) delta_tm_sw
DELTA_TM_SW: Melting point depression of sea water due to its average salinity.
subroutine, public read_erg_nc(z_sl, filename)
Reading of time-slice files in native binary or NetCDF format.
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.
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) age_c
age_c(kc,j,i): Age in the upper (kc) ice domain
integer(i2b), dimension(0:jmax, 0:imax) n_cts
n_cts(j,i): Mask for thermal conditions. -1: cold ice base, 0: temperate ice base with cold ice above...
real(dp) beta2_ht_0
BETA2_HT_0: Degree-day factor for ice at high summer temperatures.
real(dp) beta1_ht_0
BETA1_HT_0: Degree-day factor for snow at high summer temperatures.
real(dp), dimension(0:jmax, 0:imax) vx_b_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) age_t
age_t(kt,j,i): Age in the lower (kt) ice domain
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) enh_c
enh_c(kc,j,i): Flow enhancement factor in the upper (kc) ice domain
real(dp) rho_w
RHO_W: Density of pure water.
real(dp), dimension(0:krmax, 0:jmax, 0:imax) temp_r
temp_r(kr,j,i): Temperature in the bedrock
real(dp) h_r
H_R: Thickness of the modelled lithosphere layer.
real(dp) a
A: Semi-major axis of the planet.
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.
real(dp), dimension(0:jmax, 0:imax) vx_s_g
vx_s_g(j,i): Velocity in x-direction at the ice surface, at (i,j)
real(dp) omega_max
OMEGA_MAX: Threshold value for the water content of temperate ice.
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) enth_t
enth_t(kt,j,i): Enthalpy in the lower (kt) ice domain
real(dp), dimension(0:jmax, 0:imax) vz_s
vz_s(j,i): Velocity in z-direction at the ice surface, at (i,j)
real(dp), dimension(0:jmax, 0:imax) dh_t_dtau
dH_t_dtau(j,i): Derivative of H_t by tau (time)
real(dp), dimension(0:jmax, 0:imax) dh_c_dtau
dH_c_dtau(j,i): Derivative of H_c by tau (time)
integer(i2b), dimension(0:jmax, 0:imax) maske
maske(j,i): Ice-land-ocean mask. 0: grounded ice, 1: ice-free land, 2: ocean, 3: floating ice ...
real(dp), dimension(0:jmax, 0:imax) p_b_w
p_b_w(j,i): Basal water pressure
real(dp) s_stat_0
S_STAT_0: Standard deviation of the air termperature for the degree-day model.
real(dp), dimension(0:jmax, 0:imax) zl_target
zl_target(j,i): Target topography (lithosphere surface)
real(dp) beta2_lt_0
BETA2_LT_0: Degree-day factor for ice at low summer temperatures.
real(dp), dimension(0:jmax, 0:imax) zs_target
zs_target(j,i): Target topography (ice surface)
real(dp) b
B: Semi-minor axis of the planet.
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) vz_t
vz_t(kt,j,i): Velocity in z-direction in the lower (kt) ice domain (at (i,j,kt+1/2)) ...
real(dp), dimension(0:jmax, 0:imax) zl
zl(j,i): Coordinate z of the lithosphere surface
real(dp) g
G: Acceleration due to gravity.
real(dp) r
R: Radius of the planet.
integer(i2b), dimension(0:jmax, 0:imax) kc_cts
kc_cts(j,i): Position kc of the CTS (for COLD, ENTC, ENTM)
real(dp), dimension(0:jmax, 0:imax) q_tld
Q_tld(j,i): Water drainage rate from the temperate layer.
Declarations of global variables for SICOPOLIS (for the ANT domain).
subroutine read_phys_para_value(para, n_unit)
Reading of a value of a physical parameter from the phys_para file.
real(dp), dimension(0:jmax, 0:imax) temp_b
temp_b(j,i): Basal temperature
real(dp) beta1_0
BETA1_0: Degree-day factor for snow.
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) enth_c
enth_c(kc,j,i): Enthalpy in the upper (kc) ice domain
real(dp) rho_c_r
RHO_C_R: Density times specific heat of the lithosphere.
Reading of several input data.
Declarations of kind types for SICOPOLIS.
real(dp), dimension(0:jmax, 0:imax) zs
zs(j,i): Coordinate z of the surface topography
real(dp) pmax_0
PMAX_0: Saturation factor for the formation of superimposed ice.
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) enh_t
enh_t(kt,j,i): Flow enhancement factor in the lower (kt) ice domain
real(dp) rho_c
RHO_C: Density of crustal material (dust)
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.
real(dp), dimension(0:jmax, 0:imax) dzm_dtau
dzm_dtau(j,i): Derivative of zm by tau (time)
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) omega_c
omega_c(kc,j,i): Water content in the upper (kc) ice domain
integer(i2b), dimension(0:jmax, 0:imax) maske_old
maske_old(j,i): Old value of maske (at the previous time step)
real(dp), dimension(0:jmax, 0:imax) qx
qx(j,i): Volume flux in x-direction (depth-integrated vx, at (i+1/2,j))
real(dp) kappa_c
KAPPA_C: Heat conductivity of crustal material (dust)
real(dp) rho_i
RHO_I: Density of ice.
real(dp) nue
NUE: Water diffusivity in ice.
real(dp) c_c
C_C: Specific heat of crustal material (dust)
real(dp) beta1_lt_0
BETA1_LT_0: Degree-day factor for snow at low summer temperatures.
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) vz_c
vz_c(kc,j,i): Velocity in z-direction in the upper (kc) ice domain (at (i,j,kc+1/2)) ...
real(dp) phi0
PHI0: Standard parallel of the stereographic projection.
real(dp), dimension(0:jmax, 0:imax) vy_s_g
vy_s_g(j,i): Velocity in x-direction at the ice surface, at (i,j)
real(dp), dimension(0:jmax, 0:imax) zm
zm(j,i): Coordinate z of the bottom of the upper (kc) ice domain = top of the lower (kt) ice domain (...
real(dp), dimension(-190:10) rf
RF(n): Tabulated values for the rate factor of cold ice.
real(dp), dimension(0:jmax, 0:imax) q_bm
Q_bm(j,i): Basal melting rate.
real(dp), dimension(-190:10) c
C(n): Tabulated values for the specific heat of ice.
real(dp), dimension(-10000:10000) kei
kei(n): Tabulated values of the kei function (Kelvin function of zero order)
subroutine check(status, ch_calling_routine)
NetCDF error capturing.
real(dp), dimension(0:jmax, 0:imax) temph_b
temph_b(j,i): Basal temperature relative to the pressure melting point
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) vy_t
vy_t(kt,j,i): Velocity in y-direction in the lower (kt) ice domain (at (i,j+1/2,kt)) ...
real(dp), dimension(0:jmax, 0:imax) am_perp
am_perp(j,i): Ice volume flux across the z=zm interface
real(dp) beta
BETA: Clausius-Clapeyron gradient of ice.
real(dp), parameter no_value_pos_1
no_value_pos_1: Positive no-value parameter
real(dp), dimension(0:jmax) eta
eta(j): Coordinate eta (= y) of grid point j
real(dp), dimension(0:imax) xi
xi(i): Coordinate xi (= x) of grid point i
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) vy_c
vy_c(kc,j,i): Velocity in y-direction in the upper (kc) ice domain (at (i,j+1/2,kc)) ...
character(len=16) ch_domain_short
ch_domain_short: Short name of the computational domain
real(dp) r_t
R_T: Coefficient of the water-content dependence in the rate factor for temperate ice...
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) omega_t
omega_t(kt,j,i): Water content in the lower (kt) ice domain
real(dp) kei_r_incr
kei_r_incr: Increment of the argument r of the tabulated kei function
real(dp) rho_sw
RHO_SW: Density of sea water.
real(dp) l
L: Latent heat of ice.
integer, parameter sp
Single-precision reals.
integer(i4b) n_data_kei
n_data_kei: Number of tabulated values of the kei function
real(dp) rho
RHO: Density of ice.
real(dp) rho_a
RHO_A: Density of the asthenosphere.
real(dp), dimension(0:jmax, 0:imax) h_target
H_target(j,i): Target topography (ice thickness)
integer, parameter dp
Double-precision reals.
real(dp), dimension(0:jmax, 0:imax) zb
zb(j,i): Coordinate z of the ice base
real(dp), dimension(0:jmax, 0:imax) dzb_dtau
dzb_dtau(j,i): Derivative of zb by tau (time)
integer(i2b), dimension(0:jmax, 0:imax) maske_target
maske_target(j,i): Target topography (ice-land-ocean mask)
subroutine, public read_target_topo_nc(target_topo_dat_name)
Reading of the target-topography file (in NetCDF format).
real(dp) lambda0
LAMBDA0: Reference longitude (central meridian) of the stereographic projection.
real(dp), dimension(0:jmax, 0:imax) h_t
H_t(j,i): Thickness of ice in the lower (kt) domain (thickness of the temperate layer for POLY...
real(dp), dimension(0:jmax, 0:imax) h_w
H_w(j,i): Thickness of the water column under the ice base.
Declarations of global variables for SICOPOLIS.
real(dp), dimension(0:jmax, 0:imax) dzs_dtau
dzs_dtau(j,i): Derivative of zs by tau (time)
real(dp), dimension(0:jmax, 0:imax) dzl_dtau
dzl_dtau(j,i): Derivative of zl by tau (time)
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) vx_t
vx_t(kt,j,i): Velocity in x-direction in the lower (kt) ice domain (at (i+1/2,j,kt)) ...
real(dp), dimension(0:jmax, 0:imax) qy
qy(j,i): Volume flux in y-direction (depth-integrated vy, at (i,j+1/2))