SICOPOLIS V5-dev  Revision 1173
output_m.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Module : o u t p u t _ m
4 !
5 !> @file
6 !!
7 !! Writing of output data on files.
8 !!
9 !! @section Copyright
10 !!
11 !! Copyright 2009-2017 Ralf Greve, Reinhard Calov, Thomas Goelles,
12 !! Thorben Dunse
13 !!
14 !! @section License
15 !!
16 !! This file is part of SICOPOLIS.
17 !!
18 !! SICOPOLIS is free software: you can redistribute it and/or modify
19 !! it under the terms of the GNU General Public License as published by
20 !! the Free Software Foundation, either version 3 of the License, or
21 !! (at your option) any later version.
22 !!
23 !! SICOPOLIS is distributed in the hope that it will be useful,
24 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
25 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
26 !! GNU General Public License for more details.
27 !!
28 !! You should have received a copy of the GNU General Public License
29 !! along with SICOPOLIS. If not, see <http://www.gnu.org/licenses/>.
30 !<
31 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
32 
33 !-------------------------------------------------------------------------------
34 !> Writing of output data on files.
35 !<------------------------------------------------------------------------------
36 module output_m
37 
38  use sico_types_m
40  use sico_vars_m
41 
42  implicit none
43 
44  real(dp), parameter :: sec_to_year = 1.0_dp/year_sec
45 
46  private
47  public :: output1, output2, output4, borehole
48 #if (defined(ASF))
49  public :: output5
50 #endif
51 
52 contains
53 
54 !-------------------------------------------------------------------------------
55 !> Writing of time-slice files in native binary or NetCDF format.
56 !<------------------------------------------------------------------------------
57 subroutine output1(runname, time, delta_ts, glac_index, z_sl, &
58  flag_3d_output, ndat2d, ndat3d)
59 
60 #if (CALCMOD==1 || CALCMOD==0 || CALCMOD==-1)
62 #endif
63 
64 #if (NETCDF==2) /* time-slice file in NetCDF format */
65  use netcdf
66  use nc_check_m
67 #endif
68 
69 #if (DISC>0)
71 #endif
72 
73 implicit none
74 
75 real(dp), intent(in) :: time, delta_ts, glac_index, z_sl
76 character(len=100), intent(in) :: runname
77 logical, intent(in) :: flag_3d_output
78 
79 integer(i4b), intent(inout) :: ndat2d, ndat3d
80 
81 integer(i4b) :: i, j, kc, kt, kr
82 integer(i4b) :: ndat, ndat_help, ndat_1000s, ndat_100s, ndat_10s, ndat_1s
83 real(dp), dimension(0:JMAX,0:IMAX) :: H, H_cold, H_temp, dH_dtau
84 real(dp), dimension(0:JMAX,0:IMAX) :: vx_m_g, vy_m_g
85 real(dp), dimension(0:JMAX,0:IMAX) :: tau_b_driving, tau_b_drag
86 real(dp) :: V_tot, V_grounded, V_floating, V_gr_redu, V_af
87 real(dp) :: A_grounded, A_floating
88 real(sp) :: lon0, lat0
89 real(dp) :: rhosw_rho_ratio
90 character(len=256) :: filename, filename_with_path
91 character :: ch_1000s, ch_100s, ch_10s, ch_1s
92 character(len=16) :: ch_date, ch_time, ch_zone
93 
94 integer(i2b), dimension(0:IMAX,0:JMAX) :: maske_conv, maske_old_conv, &
95  mask_run_conv, n_cts_conv, kc_cts_conv
96 integer(i2b), dimension(0:IMAX,0:JMAX) :: mask_mar_conv
97 integer(i1b), dimension(0:IMAX,0:JMAX) :: flag_shelfy_stream_x_conv, &
98  flag_shelfy_stream_y_conv, &
99  flag_shelfy_stream_conv
100 
101 real(sp) :: time_conv, delta_ts_conv, glac_index_conv, z_sl_conv, &
102  V_tot_conv, V_af_conv, A_grounded_conv, A_floating_conv, &
103  H_R_conv, &
104  xi_conv(0:imax), eta_conv(0:jmax), &
105  sigma_level_c_conv(0:kcmax), sigma_level_t_conv(0:ktmax), &
106  sigma_level_r_conv(0:krmax)
107 real(sp), dimension(0:IMAX,0:JMAX) :: lambda_conv, phi_conv, &
108  lon_conv, lat_conv, &
109  temp_s_conv, accum_conv, &
110  as_perp_conv, as_perp_apl_conv, smb_corr_conv, &
111  q_geo_conv, &
112  zs_conv, zm_conv, zb_conv, zl_conv, zl0_conv, &
113  H_cold_conv, H_temp_conv, H_conv, &
114  Q_bm_conv, Q_tld_conv, &
115  am_perp_conv, &
116  qx_conv, qy_conv, &
117  dzs_dtau_conv, dzm_dtau_conv, dzb_dtau_conv, dzl_dtau_conv, &
118  dH_c_dtau_conv, dH_t_dtau_conv, dH_dtau_conv, &
119  vx_b_g_conv, vy_b_g_conv, vz_b_conv, vh_b_conv, &
120  vx_s_g_conv, vy_s_g_conv, vz_s_conv, vh_s_conv, &
121  vx_m_g_conv, vy_m_g_conv, vh_m_conv, &
122  temp_b_conv, temph_b_conv, &
123  tau_b_driving_conv, tau_b_drag_conv, &
124  p_b_w_conv, H_w_conv, q_gl_g_conv, q_cf_g_conv, &
125  cst_dist_conv, cos_grad_tc_conv, dis_perp_conv, &
126  ratio_sl_x_conv, ratio_sl_y_conv, &
127  vis_int_g_conv
128 
129 real(sp), dimension(0:IMAX,0:JMAX,0:KCMAX) :: vx_c_conv, vy_c_conv, vz_c_conv, &
130  temp_c_conv, age_c_conv, &
131  enth_c_conv, omega_c_conv, &
132  enh_c_conv
133 real(sp), dimension(0:IMAX,0:JMAX,0:KTMAX) :: vx_t_conv, vy_t_conv, vz_t_conv, &
134  omega_t_conv, age_t_conv, &
135  enth_t_conv, &
136  enh_t_conv
137 real(sp), dimension(0:IMAX,0:JMAX,0:KRMAX) :: temp_r_conv
138 
139 #if (NETCDF==1) /* time-slice file in native binary format */
140 
141 integer(i4b) :: ios
142 character(len=256) :: ch_attr_title, ch_attr_institution, ch_attr_source, &
143  ch_attr_history, ch_attr_references
144 character(len= 16), parameter :: filename_extension = '.erg'
145 
146 #elif (NETCDF==2) /* time-slice file in NetCDF format */
147 
148 integer(i4b) :: ncid, ncv
149 ! ncid: ID of the output file
150 ! ncv: Variable ID
151 integer(i4b) :: ncd, nc1d, nc2d(2), nc3d(3)
152 ! ncd: Dimension ID
153 ! nc1d: Dimension of a 1-d array
154 ! nc2d: Vector with the dimensions of a 2-d array
155 ! nc3d: Vector with the dimensions of a 3-d array
156 integer(i4b) :: nc2flag(2), nc3flag(3), nc4flag(4)
157 ! nc2flag: Vector with the 2 possible values of a flag variable
158 ! nc3flag: Vector with the 3 possible values of a flag variable
159 ! nc4flag: Vector with the 4 possible values of a flag variable
160 integer(i4b) :: nc1cor_i(1), nc1cor_j(1), &
161  nc1cor_kc(1), nc1cor_kt(1), nc1cor_kr(1), &
162  nc2cor_ij(2), &
163  nc3cor_ijkc(3), nc3cor_ijkt(3), nc3cor_ijkr(3)
164 ! nc1cor(1): Corner of a 1-d array
165 ! nc2cor(2): Corner of a 2-d array
166 ! nc3cor(3): Corner of a 3-d array
167 integer(i4b) :: nc1cnt_i(1), nc1cnt_j(1), &
168  nc1cnt_kc(1), nc1cnt_kt(1), nc1cnt_kr(1), &
169  nc2cnt_ij(2), &
170  nc3cnt_ijkc(3), nc3cnt_ijkt(3), nc3cnt_ijkr(3)
171 ! nc1cnt(1): Count of a 1-d array
172 ! nc2cnt(2): Count of a 2-d array
173 ! nc3cnt(3): Count of a 3-d array
174 character(len=256) :: buffer
175 character(len= 16), parameter :: filename_extension = '.nc'
176 character(len= 16), allocatable :: coord_id(:)
177 
178 #else
179 stop ' >>> output1: Parameter NETCDF must be either 1 or 2!'
180 #endif
181 
182 character(len=64), parameter :: thisroutine = 'output1'
183 
184 #if (NETCDF==2) /* time-slice file in NetCDF format */
185 
186 nc1cor_i = (/ 1 /)
187 nc1cnt_i = (/ imax+1 /)
188 
189 nc1cor_j = (/ 1 /)
190 nc1cnt_j = (/ jmax+1 /)
191 
192 nc1cor_kc = (/ 1 /)
193 nc1cnt_kc = (/ kcmax+1 /)
194 
195 nc1cor_kt = (/ 1 /)
196 nc1cnt_kt = (/ ktmax+1 /)
197 
198 nc1cor_kr = (/ 1 /)
199 nc1cnt_kr = (/ krmax+1 /)
200 
201 nc2cor_ij = (/ 1, 1 /)
202 nc2cnt_ij = (/ imax+1, jmax+1 /)
203 
204 nc3cor_ijkc = (/ 1, 1, 1 /)
205 nc3cnt_ijkc = (/ imax+1, jmax+1, kcmax+1 /)
206 
207 nc3cor_ijkt = (/ 1, 1, 1 /)
208 nc3cnt_ijkt = (/ imax+1, jmax+1, ktmax+1 /)
209 
210 nc3cor_ijkr = (/ 1, 1, 1 /)
211 nc3cnt_ijkr = (/ imax+1, jmax+1, krmax+1 /)
212 
213 #endif
214 
215 !-------- Create consecutively numbered file names --------
216 
217 if (flag_3d_output) then
218  ndat = ndat3d
219 else
220  ndat = ndat2d
221 end if
222 
223 if (ndat > 9999) stop ' >>> output1: Too many time-slice files!'
224 
225 ndat_help = ndat
226 ndat_1000s = ndat_help/1000
227 ndat_help = ndat_help-ndat_1000s*1000
228 ndat_100s = ndat_help/100
229 ndat_help = ndat_help-ndat_100s*100
230 ndat_10s = ndat_help/10
231 ndat_help = ndat_help-ndat_10s*10
232 ndat_1s = ndat_help
233 
234 ch_1000s = char(ndat_1000s+ichar('0'))
235 ch_100s = char(ndat_100s +ichar('0'))
236 ch_10s = char(ndat_10s +ichar('0'))
237 ch_1s = char(ndat_1s +ichar('0'))
238 
239 if (flag_3d_output) then
240  filename = trim(runname)//ch_1000s//ch_100s//ch_10s//ch_1s &
241  //trim(filename_extension)
242 else
243  filename = trim(runname)//'_2d_'//ch_1000s//ch_100s//ch_10s//ch_1s &
244  //trim(filename_extension)
245 end if
246 
247 filename_with_path = trim(outpath)//'/'//trim(filename)
248 
249 !-------- File initialization --------
250 
251 #if (NETCDF==1) /* time-slice file in native binary format */
252 
253 ! ------ Open native binary file
254 
255 open(unit=11, iostat=ios, file=trim(filename_with_path), status='new', &
256  form='unformatted')
257 
258 if (ios /= 0) stop ' >>> output1: Error when opening an erg file!'
259 
260 ! ------ Global attributes
261 
262 ch_attr_title = 'Time-slice output no. '//ch_1000s//ch_100s//ch_10s//ch_1s// &
263  ' of simulation '//trim(runname)
264 write(unit=11) ch_attr_title
265 
266 ch_attr_institution = 'Institute of Low Temperature Science, '// &
267  'Hokkaido University, Sapporo, Japan'
268 write(unit=11) ch_attr_institution
269 
270 ch_attr_source = 'SICOPOLIS Version '//version
271 write(unit=11) ch_attr_source
272 
273 call date_and_time(ch_date, ch_time, ch_zone)
274 ch_attr_history = ch_date(1:4)//'-'//ch_date(5:6)//'-'//ch_date(7:8)//' '// &
275  ch_time(1:2)//':'//ch_time(3:4)//':'//ch_time(5:6)//' '// &
276  ch_zone(1:3)//':'//ch_zone(4:5)//' - Data produced'
277 write(unit=11) ch_attr_history
278 
279 ch_attr_references = 'http://www.sicopolis.net/'
280 write(unit=11) ch_attr_references
281 
282 #elif (NETCDF==2) /* time-slice file in NetCDF format */
283 
284 if (allocated(coord_id)) deallocate(coord_id); allocate(coord_id(5))
285 coord_id(1) = 'x'; coord_id(2) = 'y'
286 coord_id(3) = 'zeta_c'; coord_id(4) = 'zeta_t'; coord_id(5) = 'zeta_r'
287 
288 ! ------ Open NetCDF file
289 
290 call check( nf90_create(trim(filename_with_path), nf90_noclobber, ncid), &
291  thisroutine )
292 
293 ! ------ Global attributes
294 
295 buffer = 'Time-slice output no. '//ch_1000s//ch_100s//ch_10s//ch_1s// &
296  ' of simulation '//trim(runname)
297 call check( nf90_put_att(ncid, nf90_global, 'title', trim(buffer)), &
298  thisroutine )
299 
300 buffer = 'Institute of Low Temperature Science, Hokkaido University, '// &
301  'Sapporo, Japan'
302 call check( nf90_put_att(ncid, nf90_global, 'institution', trim(buffer)), &
303  thisroutine )
304 
305 buffer = 'SICOPOLIS Version '//version
306 call check( nf90_put_att(ncid, nf90_global, 'source', trim(buffer)), &
307  thisroutine )
308 
309 call date_and_time(ch_date, ch_time, ch_zone)
310 buffer = ch_date(1:4)//'-'//ch_date(5:6)//'-'//ch_date(7:8)//' '// &
311  ch_time(1:2)//':'//ch_time(3:4)//':'//ch_time(5:6)//' '// &
312  ch_zone(1:3)//':'//ch_zone(4:5)//' - Data produced'
313 call check( nf90_put_att(ncid, nf90_global, 'history', trim(buffer)), &
314  thisroutine )
315 
316 buffer = 'http://www.sicopolis.net/'
317 call check( nf90_put_att(ncid, nf90_global, 'references', trim(buffer)), &
318  thisroutine )
319 
320 ! ------ Definition of the dimensions
321 
322 call check( nf90_def_dim(ncid, trim(coord_id(1)), imax+1, ncd), thisroutine )
323 call check( nf90_def_dim(ncid, trim(coord_id(2)), jmax+1, ncd), thisroutine )
324 call check( nf90_def_dim(ncid, trim(coord_id(3)), kcmax+1, ncd), thisroutine )
325 call check( nf90_def_dim(ncid, trim(coord_id(4)), ktmax+1, ncd), thisroutine )
326 call check( nf90_def_dim(ncid, trim(coord_id(5)), krmax+1, ncd), thisroutine )
327 
328 ! ------ Definition of the variables
329 
330 ! ---- crs
331 
332 call check( nf90_def_var(ncid, 'crs', nf90_short, ncv), thisroutine )
333 #if (GRID==0 || GRID==1)
334 buffer = 'polar_stereographic'
335 call check( nf90_put_att(ncid, ncv, 'grid_mapping_name', trim(buffer)), &
336  thisroutine )
337 #elif (GRID==2)
338 buffer = 'latitude_longitude'
339 call check( nf90_put_att(ncid, ncv, 'grid_mapping_name', trim(buffer)), &
340  thisroutine )
341 #endif
342 #if (defined(ANT) \
343  || defined(asf) \
344  || defined(grl) \
345  || defined(scand) \
346  || defined(nhem) \
347  || defined(tibet) \
348  || defined(emtp2sge) \
349  || defined(xyz)) /* terrestrial ice sheet */
350 buffer = 'WGS84'
351 call check( nf90_put_att(ncid, ncv, 'ellipsoid', trim(buffer)), &
352  thisroutine )
353 #elif (defined(NMARS) || defined(SMARS)) /* Martian ice sheet */
354 buffer = 'Mars_ellipsoid'
355 call check( nf90_put_att(ncid, ncv, 'ellipsoid', trim(buffer)), &
356  thisroutine )
357 #else
358 stop ' >>> output1: No valid domain (ANT, GRL etc.) specified!'
359 #endif
360 #if (GRID==0 || GRID==1)
361 call check( nf90_put_att(ncid, ncv, 'false_easting', 0.0), &
362  thisroutine )
363 call check( nf90_put_att(ncid, ncv, 'false_northing', 0.0), &
364  thisroutine )
365 lon0 = lambda0 *pi_180_inv
366 lon0 = modulo(lon0+180.0_sp, 360.0_sp)-180.0_sp
367 lon0 = nint(lon0*1.0e+04_sp)*1.0e-04_sp
368 lat0 = phi0 *pi_180_inv
369 if (lat0 > 90.0_sp) lat0 = 90.0_sp
370 if (lat0 < -90.0_sp) lat0 = -90.0_sp
371 lat0 = nint(lat0*1.0e+04_sp)*1.0e-04_sp
372  ! reference longitude and standard parallel in deg rounded to 4 digits
373 if (lat0 >= 0.0_sp) then
374  call check( nf90_put_att(ncid, ncv, &
375  'latitude_of_projection_origin', 90.0), &
376  thisroutine )
377 else
378  call check( nf90_put_att(ncid, ncv, &
379  'latitude_of_projection_origin', -90.0), &
380  thisroutine )
381 end if
382 call check( nf90_put_att(ncid, ncv, &
383  'straight_vertical_longitude_from_pole', lon0), &
384  thisroutine )
385 call check( nf90_put_att(ncid, ncv, &
386  'standard_parallel', lat0), &
387  thisroutine )
388 #endif
389 
390 ! ---- time
391 
392 call check( nf90_def_var(ncid, 'time', nf90_float, ncv), &
393  thisroutine )
394 buffer = 'a'
395 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
396  thisroutine )
397 buffer = 'time'
398 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
399  thisroutine )
400 buffer = 'Time'
401 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
402  thisroutine )
403 
404 if (forcing_flag == 1) then
405 
406 ! ---- delta_ts
407 
408  call check( nf90_def_var(ncid, 'delta_ts', nf90_float, ncv), &
409  thisroutine )
410  buffer = 'degC'
411  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
412  thisroutine )
413  buffer = 'surface_temperature_anomaly'
414  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
415  thisroutine )
416  buffer = 'Surface temperature anomaly'
417  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
418  thisroutine )
419 
420 else if (forcing_flag == 2) then
421 
422 ! ---- glac_index
423 
424  call check( nf90_def_var(ncid, 'glac_index', nf90_float, ncv), &
425  thisroutine )
426  buffer = '1'
427  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
428  thisroutine )
429  buffer = 'glacial_index'
430  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
431  thisroutine )
432  buffer = 'Glacial index'
433  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
434  thisroutine )
435 
436 else if (forcing_flag == 3) then
437 
438 ! ---- glac_index
439 
440  call check( nf90_def_var(ncid, 'glac_index', nf90_float, ncv), &
441  thisroutine )
442  buffer = '1'
443  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
444  thisroutine )
445  buffer = 'glacial_index'
446  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
447  thisroutine )
448  buffer = 'Glacial index'
449  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
450  thisroutine )
451  buffer = 'This variable will be assigned a dummy value only!'
452  call check( nf90_put_att(ncid, ncv, 'comment', trim(buffer)), &
453  thisroutine )
454 
455 end if
456 
457 ! ---- z_sl
458 
459 call check( nf90_def_var(ncid, 'z_sl', nf90_float, ncv), &
460  thisroutine )
461 buffer = 'm'
462 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
463  thisroutine )
464 buffer = 'global_average_sea_level_change'
465 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
466  thisroutine )
467 buffer = 'Sea level'
468 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
469  thisroutine )
470 
471 ! ---- V_tot
472 
473 call check( nf90_def_var(ncid, 'V_tot', nf90_float, ncv), &
474  thisroutine )
475 buffer = 'm3'
476 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
477  thisroutine )
478 buffer = 'land_ice_volume'
479 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
480  thisroutine )
481 buffer = 'Ice volume'
482 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
483  thisroutine )
484 
485 ! ---- V_af
486 
487 call check( nf90_def_var(ncid, 'V_af', nf90_float, ncv), &
488  thisroutine )
489 buffer = 'm3'
490 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
491  thisroutine )
492 buffer = 'land_ice_volume_not_displacing_sea_water'
493 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
494  thisroutine )
495 buffer = 'Ice volume above flotation'
496 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
497  thisroutine )
498 
499 ! ---- A_grounded
500 
501 call check( nf90_def_var(ncid, 'A_grounded', nf90_float, ncv), &
502  thisroutine )
503 buffer = 'm2'
504 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
505  thisroutine )
506 buffer = 'grounded_land_ice_area'
507 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
508  thisroutine )
509 buffer = 'Area covered by grounded ice'
510 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
511  thisroutine )
512 
513 ! ---- A_floating
514 
515 call check( nf90_def_var(ncid, 'A_floating', nf90_float, ncv), &
516  thisroutine )
517 buffer = 'm2'
518 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
519  thisroutine )
520 buffer = 'floating_ice_shelf_area'
521 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
522  thisroutine )
523 buffer = 'Area covered by floating ice'
524 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
525  thisroutine )
526 
527 ! ---- x (= xi)
528 
529 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc1d), &
530  thisroutine )
531 call check( nf90_def_var(ncid, 'x', nf90_float, nc1d, ncv), &
532  thisroutine )
533 buffer = 'm'
534 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
535  thisroutine )
536 buffer = 'projection_x_coordinate'
537 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
538  thisroutine )
539 buffer = 'x-coordinate of the grid point i'
540 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
541  thisroutine )
542 call check( nf90_put_att(ncid, ncv, 'axis', 'x'), &
543  thisroutine )
544 
545 ! ---- y (= eta)
546 
547 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc1d), &
548  thisroutine )
549 call check( nf90_def_var(ncid, 'y', nf90_float, nc1d, ncv), &
550  thisroutine )
551 buffer = 'm'
552 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
553  thisroutine )
554 buffer = 'projection_y_coordinate'
555 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
556  thisroutine )
557 buffer = 'y-coordinate of the grid point j'
558 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
559  thisroutine )
560 call check( nf90_put_att(ncid, ncv, 'axis', 'y'), &
561  thisroutine )
562 
563 ! ---- sigma_level_c
564 
565 call check( nf90_inq_dimid(ncid, trim(coord_id(3)), nc1d), &
566  thisroutine )
567 call check( nf90_def_var(ncid, 'sigma_level_c', nf90_float, nc1d, ncv), &
568  thisroutine )
569 buffer = 'up'
570 call check( nf90_put_att(ncid, ncv, 'positive', trim(buffer)), &
571  thisroutine )
572 buffer = 'land_ice_kc_layer_sigma_coordinate'
573 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
574  thisroutine )
575 buffer = 'sigma-coordinate of the grid point kc'
576 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
577  thisroutine )
578 
579 ! ---- sigma_level_t
580 
581 call check( nf90_inq_dimid(ncid, trim(coord_id(4)), nc1d), &
582  thisroutine )
583 call check( nf90_def_var(ncid, 'sigma_level_t', nf90_float, nc1d, ncv), &
584  thisroutine )
585 buffer = 'up'
586 call check( nf90_put_att(ncid, ncv, 'positive', trim(buffer)), &
587  thisroutine )
588 buffer = 'land_ice_kt_layer_sigma_coordinate'
589 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
590  thisroutine )
591 buffer = 'sigma-coordinate of the grid point kt'
592 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
593  thisroutine )
594 
595 ! ---- sigma_level_r
596 
597 call check( nf90_inq_dimid(ncid, trim(coord_id(5)), nc1d), &
598  thisroutine )
599 call check( nf90_def_var(ncid, 'sigma_level_r', nf90_float, nc1d, ncv), &
600  thisroutine )
601 buffer = 'up'
602 call check( nf90_put_att(ncid, ncv, 'positive', trim(buffer)), &
603  thisroutine )
604 buffer = 'lithosphere_layer_sigma_coordinate'
605 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
606  thisroutine )
607 buffer = 'sigma-coordinate of the grid point kr'
608 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
609  thisroutine )
610 
611 ! ---- lon
612 
613 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
614  thisroutine )
615 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
616  thisroutine )
617 call check( nf90_def_var(ncid, 'lon', nf90_float, nc2d, ncv), &
618  thisroutine )
619 buffer = 'degrees_E'
620 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
621  thisroutine )
622 buffer = 'longitude'
623 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
624  thisroutine )
625 buffer = 'Geographical longitude'
626 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
627  thisroutine )
628 
629 ! ---- lat
630 
631 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
632  thisroutine )
633 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
634  thisroutine )
635 call check( nf90_def_var(ncid, 'lat', nf90_float, nc2d, ncv), &
636  thisroutine )
637 buffer = 'degrees_N'
638 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
639  thisroutine )
640 buffer = 'latitude'
641 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
642  thisroutine )
643 buffer = 'Geographical latitude'
644 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
645  thisroutine )
646 
647 ! ---- lambda
648 
649 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
650  thisroutine )
651 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
652  thisroutine )
653 call check( nf90_def_var(ncid, 'lambda', nf90_float, nc2d, ncv), &
654  thisroutine )
655 buffer = 'rad'
656 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
657  thisroutine )
658 buffer = 'longitude'
659 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
660  thisroutine )
661 buffer = 'Geographical longitude'
662 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
663  thisroutine )
664 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
665  thisroutine )
666 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
667  thisroutine )
668 
669 ! ---- phi
670 
671 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
672  thisroutine )
673 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
674  thisroutine )
675 call check( nf90_def_var(ncid, 'phi', nf90_float, nc2d, ncv), &
676  thisroutine )
677 buffer = 'rad'
678 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
679  thisroutine )
680 buffer = 'latitude'
681 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
682  thisroutine )
683 buffer = 'Geographical latitude'
684 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
685  thisroutine )
686 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
687  thisroutine )
688 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
689  thisroutine )
690 
691 ! ---- temp_s
692 
693 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
694  thisroutine )
695 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
696  thisroutine )
697 call check( nf90_def_var(ncid, 'temp_s', nf90_float, nc2d, ncv), &
698  thisroutine )
699 buffer = 'degC'
700 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
701  thisroutine )
702 buffer = 'surface_temperature'
703 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
704  thisroutine )
705 buffer = 'Temperature at the ice surface'
706 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
707  thisroutine )
708 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
709  thisroutine )
710 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
711  thisroutine )
712 
713 ! ---- accum
714 
715 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
716  thisroutine )
717 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
718  thisroutine )
719 call check( nf90_def_var(ncid, 'prec', nf90_float, nc2d, ncv), &
720  thisroutine )
721 buffer = 'm a-1'
722 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
723  thisroutine )
724 buffer = 'land_ice_precipitation'
725 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
726  thisroutine )
727 buffer = 'Annual precipitation at the ice surface'
728 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
729  thisroutine )
730 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
731  thisroutine )
732 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
733  thisroutine )
734 
735 ! ---- as_perp
736 
737 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
738  thisroutine )
739 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
740  thisroutine )
741 call check( nf90_def_var(ncid, 'as_perp', nf90_float, nc2d, ncv), &
742  thisroutine )
743 buffer = 'm a-1'
744 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
745  thisroutine )
746 buffer = 'land_ice_surface_mass_balance'
747 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
748  thisroutine )
749 buffer = 'Mass balance at the ice surface'
750 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
751  thisroutine )
752 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
753  thisroutine )
754 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
755  thisroutine )
756 
757 ! ---- as_perp_apl
758 
759 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
760  thisroutine )
761 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
762  thisroutine )
763 call check( nf90_def_var(ncid, 'as_perp_apl', nf90_float, nc2d, ncv), &
764  thisroutine )
765 buffer = 'm a-1'
766 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
767  thisroutine )
768 buffer = 'applied_land_ice_surface_mass_balance'
769 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
770  thisroutine )
771 buffer = 'Applied mass balance at the ice surface'
772 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
773  thisroutine )
774 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
775  thisroutine )
776 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
777  thisroutine )
778 
779 ! ---- smb_corr
780 
781 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
782  thisroutine )
783 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
784  thisroutine )
785 call check( nf90_def_var(ncid, 'smb_corr', nf90_float, nc2d, ncv), &
786  thisroutine )
787 buffer = 'm a-1'
788 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
789  thisroutine )
790 buffer = 'land_ice_surface_mass_balance_diagnosed_correction'
791 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
792  thisroutine )
793 buffer = 'Diagnosed correction of the mass balance at the ice surface'
794 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
795  thisroutine )
796 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
797  thisroutine )
798 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
799  thisroutine )
800 
801 #if (DISC>0) /* Ice discharge parameterisation */
802 
803 ! ---- dis_perp
804 
805 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
806  thisroutine )
807 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
808  thisroutine )
809 call check( nf90_def_var(ncid, 'dis_perp', nf90_float, nc2d, ncv), &
810  thisroutine )
811 buffer = 'm a-1'
812 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
813  thisroutine )
814 buffer = 'ice_discharge'
815 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
816  thisroutine )
817 buffer = 'Ice discharge'
818 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
819  thisroutine )
820 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
821  thisroutine )
822 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
823  thisroutine )
824 
825 ! ---- cst_dist
826 
827 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
828  thisroutine )
829 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
830  thisroutine )
831 call check( nf90_def_var(ncid, 'cst_dist', nf90_float, nc2d, ncv), &
832  thisroutine )
833 buffer = 'km'
834 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
835  thisroutine )
836 buffer = 'coastal_distance'
837 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
838  thisroutine )
839 buffer = 'Coastal distance'
840 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
841  thisroutine )
842 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
843  thisroutine )
844 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
845  thisroutine )
846 
847 ! ---- cos_grad_tc
848 
849 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
850  thisroutine )
851 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
852  thisroutine )
853 call check( nf90_def_var(ncid, 'cos_grad_tc', nf90_float, nc2d, ncv), &
854  thisroutine )
855 buffer = '1'
856 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
857  thisroutine )
858 buffer = 'cos_alpha'
859 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
860  thisroutine )
861 buffer = 'Cosine of angle between surface gradient and cst dist gradient'
862 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
863  thisroutine )
864 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
865  thisroutine )
866 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
867  thisroutine )
868 
869 ! ---- mask_mar
870 
871 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
872  thisroutine )
873 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
874  thisroutine )
875 call check( nf90_def_var(ncid, 'mask_mar', nf90_short, nc2d, ncv), &
876  thisroutine )
877 buffer = 'marginal_ring_mask'
878 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
879  thisroutine )
880 buffer = 'Marginal ring mask'
881 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
882  thisroutine )
883 nc2flag = (/ 0, 1 /)
884 call check( nf90_put_att(ncid, ncv, 'flag_values', nc2flag), &
885  thisroutine )
886 buffer = 'no_ring '// &
887  'ring'
888 call check( nf90_put_att(ncid, ncv, 'flag_meanings', trim(buffer)), &
889  thisroutine )
890 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
891  thisroutine )
892 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
893  thisroutine )
894 
895 #endif
896 
897 ! ---- q_geo
898 
899 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
900  thisroutine )
901 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
902  thisroutine )
903 call check( nf90_def_var(ncid, 'q_geo', nf90_float, nc2d, ncv), &
904  thisroutine )
905 buffer = 'W m-2'
906 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
907  thisroutine )
908 buffer = 'upward_geothermal_heat_flux_at_ground_level'
909 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
910  thisroutine )
911 buffer = 'Geothermal heat flux'
912 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
913  thisroutine )
914 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
915  thisroutine )
916 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
917  thisroutine )
918 
919 ! ---- maske
920 
921 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
922  thisroutine )
923 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
924  thisroutine )
925 call check( nf90_def_var(ncid, 'maske', nf90_short, nc2d, ncv), &
926  thisroutine )
927 buffer = 'ice_land_sea_mask'
928 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
929  thisroutine )
930 buffer = 'Ice-land-sea mask'
931 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
932  thisroutine )
933 nc4flag = (/ 0, 1, 2, 3 /)
934 call check( nf90_put_att(ncid, ncv, 'flag_values', nc4flag), &
935  thisroutine )
936 buffer = 'glaciated_land '// &
937  'ice_free_land '// &
938  'sea '// &
939  'floating_ice'
940 call check( nf90_put_att(ncid, ncv, 'flag_meanings', trim(buffer)), &
941  thisroutine )
942 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
943  thisroutine )
944 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
945  thisroutine )
946 
947 ! ---- maske_old
948 
949 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
950  thisroutine )
951 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
952  thisroutine )
953 call check( nf90_def_var(ncid, 'maske_old', nf90_short, nc2d, ncv), &
954  thisroutine )
955 buffer = 'ice_land_sea_mask_old'
956 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
957  thisroutine )
958 buffer = 'Ice-land-sea mask (old)'
959 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
960  thisroutine )
961 nc4flag = (/ 0, 1, 2, 3 /)
962 call check( nf90_put_att(ncid, ncv, 'flag_values', nc4flag), &
963  thisroutine )
964 buffer = 'glaciated_land '// &
965  'ice_free_land '// &
966  'sea '// &
967  'floating_ice'
968 call check( nf90_put_att(ncid, ncv, 'flag_meanings', trim(buffer)), &
969  thisroutine )
970 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
971  thisroutine )
972 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
973  thisroutine )
974 
975 ! ---- mask_run
976 
977 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
978  thisroutine )
979 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
980  thisroutine )
981 call check( nf90_def_var(ncid, 'mask_run', nf90_short, nc2d, ncv), &
982  thisroutine )
983 buffer = 'indicating_melt_type_mask'
984 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
985  thisroutine )
986 buffer = 'Mask indicating melt type'
987 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
988  thisroutine )
989 nc4flag = (/ -2, -1, 1, 2 /)
990 call check( nf90_put_att(ncid, ncv, 'flag_values', nc4flag), &
991  thisroutine )
992 buffer = 'hidden in ocean '// &
993  'hidden on land '// &
994  'visible on land '// &
995  'visible in ocean'
996 call check( nf90_put_att(ncid, ncv, 'flag_meanings', trim(buffer)), &
997  thisroutine )
998 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
999  thisroutine )
1000 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1001  thisroutine )
1002 
1003 ! ---- n_cts
1004 
1005 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1006  thisroutine )
1007 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1008  thisroutine )
1009 call check( nf90_def_var(ncid, 'n_cts', nf90_short, nc2d, ncv), &
1010  thisroutine )
1011 buffer = 'polythermal_condition_mask'
1012 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1013  thisroutine )
1014 buffer = 'Mask for polythermal conditions'
1015 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1016  thisroutine )
1017 nc3flag = (/ -1, 0, 1 /)
1018 call check( nf90_put_att(ncid, ncv, 'flag_values', nc3flag), &
1019  thisroutine )
1020 buffer = 'cold_base '// &
1021  'temperate_base_with_cold_ice_above '// &
1022  'temperate_base_with_temperate_ice_above'
1023 call check( nf90_put_att(ncid, ncv, 'flag_meanings', trim(buffer)), &
1024  thisroutine )
1025 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1026  thisroutine )
1027 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1028  thisroutine )
1029 
1030 ! ---- kc_cts
1031 
1032 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1033  thisroutine )
1034 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1035  thisroutine )
1036 call check( nf90_def_var(ncid, 'kc_cts', nf90_short, nc2d, ncv), &
1037  thisroutine )
1038 buffer = 'CTS_position_grid_index'
1039 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1040  thisroutine )
1041 buffer = 'Grid index of the CTS position'
1042 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1043  thisroutine )
1044 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1045  thisroutine )
1046 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1047  thisroutine )
1048 
1049 ! ---- zs
1050 
1051 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1052  thisroutine )
1053 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1054  thisroutine )
1055 call check( nf90_def_var(ncid, 'zs', nf90_float, nc2d, ncv), &
1056  thisroutine )
1057 buffer = 'm'
1058 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1059  thisroutine )
1060 buffer = 'surface_altitude'
1061 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1062  thisroutine )
1063 buffer = 'Topography of the free surface'
1064 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1065  thisroutine )
1066 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1067  thisroutine )
1068 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1069  thisroutine )
1070 
1071 ! ---- zm
1072 
1073 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1074  thisroutine )
1075 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1076  thisroutine )
1077 call check( nf90_def_var(ncid, 'zm', nf90_float, nc2d, ncv), &
1078  thisroutine )
1079 buffer = 'm'
1080 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1081  thisroutine )
1082 buffer = 'zm_interface_altitude'
1083 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1084  thisroutine )
1085 buffer = 'Topography of the z=zm interface'
1086 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1087  thisroutine )
1088 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1089  thisroutine )
1090 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1091  thisroutine )
1092 
1093 ! ---- zb
1094 
1095 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1096  thisroutine )
1097 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1098  thisroutine )
1099 call check( nf90_def_var(ncid, 'zb', nf90_float, nc2d, ncv), &
1100  thisroutine )
1101 buffer = 'm'
1102 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1103  thisroutine )
1104 buffer = 'ice_base_altitude'
1105 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1106  thisroutine )
1107 buffer = 'Topography of the ice base'
1108 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1109  thisroutine )
1110 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1111  thisroutine )
1112 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1113  thisroutine )
1114 
1115 ! ---- zl
1116 
1117 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1118  thisroutine )
1119 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1120  thisroutine )
1121 call check( nf90_def_var(ncid, 'zl', nf90_float, nc2d, ncv), &
1122  thisroutine )
1123 buffer = 'm'
1124 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1125  thisroutine )
1126 buffer = 'bedrock_altitude'
1127 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1128  thisroutine )
1129 buffer = 'Topography of the lithosphere surface'
1130 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1131  thisroutine )
1132 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1133  thisroutine )
1134 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1135  thisroutine )
1136 
1137 ! ---- zl0
1138 
1139 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1140  thisroutine )
1141 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1142  thisroutine )
1143 call check( nf90_def_var(ncid, 'zl0', nf90_float, nc2d, ncv), &
1144  thisroutine )
1145 buffer = 'm'
1146 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1147  thisroutine )
1148 buffer = 'isostatically_relaxed_bedrock_altitude'
1149 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1150  thisroutine )
1151 buffer = 'Topography of the isostatically relaxed lithosphere surface'
1152 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1153  thisroutine )
1154 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1155  thisroutine )
1156 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1157  thisroutine )
1158 
1159 ! ---- H_cold
1160 
1161 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1162  thisroutine )
1163 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1164  thisroutine )
1165 call check( nf90_def_var(ncid, 'H_cold', nf90_float, nc2d, ncv), &
1166  thisroutine )
1167 buffer = 'm'
1168 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1169  thisroutine )
1170 buffer = 'land_ice_cold_layer_thickness'
1171 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1172  thisroutine )
1173 buffer = 'Thickness of the cold ice layer'
1174 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1175  thisroutine )
1176 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1177  thisroutine )
1178 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1179  thisroutine )
1180 
1181 ! ---- H_temp
1182 
1183 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1184  thisroutine )
1185 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1186  thisroutine )
1187 call check( nf90_def_var(ncid, 'H_temp', nf90_float, nc2d, ncv), &
1188  thisroutine )
1189 buffer = 'm'
1190 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1191  thisroutine )
1192 buffer = 'land_ice_temperate_layer_thickness'
1193 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1194  thisroutine )
1195 buffer = 'Thickness of the temperate ice layer'
1196 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1197  thisroutine )
1198 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1199  thisroutine )
1200 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1201  thisroutine )
1202 
1203 ! ---- H
1204 
1205 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1206  thisroutine )
1207 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1208  thisroutine )
1209 call check( nf90_def_var(ncid, 'H', nf90_float, nc2d, ncv), &
1210  thisroutine )
1211 buffer = 'm'
1212 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1213  thisroutine )
1214 buffer = 'land_ice_thickness'
1215 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1216  thisroutine )
1217 buffer = 'Ice thickness'
1218 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1219  thisroutine )
1220 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1221  thisroutine )
1222 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1223  thisroutine )
1224 
1225 ! ---- H_R
1226 
1227 call check( nf90_def_var(ncid, 'H_R', nf90_float, ncv), &
1228  thisroutine )
1229 buffer = 'm'
1230 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1231  thisroutine )
1232 buffer = 'lithosphere_layer_thickness'
1233 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1234  thisroutine )
1235 buffer = 'Thickness of the lithosphere layer'
1236 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1237  thisroutine )
1238 
1239 ! ---- Q_bm
1240 
1241 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1242  thisroutine )
1243 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1244  thisroutine )
1245 call check( nf90_def_var(ncid, 'Q_bm', nf90_float, nc2d, ncv), &
1246  thisroutine )
1247 buffer = 'm a-1'
1248 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1249  thisroutine )
1250 buffer = 'land_ice_basal_melt_rate'
1251 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1252  thisroutine )
1253 buffer = 'Basal melting rate'
1254 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1255  thisroutine )
1256 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1257  thisroutine )
1258 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1259  thisroutine )
1260 
1261 ! ---- Q_tld
1262 
1263 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1264  thisroutine )
1265 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1266  thisroutine )
1267 call check( nf90_def_var(ncid, 'Q_tld', nf90_float, nc2d, ncv), &
1268  thisroutine )
1269 buffer = 'm a-1'
1270 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1271  thisroutine )
1272 buffer = 'land_ice_temperate_layer_water_drainage'
1273 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1274  thisroutine )
1275 buffer = 'Water drainage from the temperate layer'
1276 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1277  thisroutine )
1278 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1279  thisroutine )
1280 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1281  thisroutine )
1282 
1283 ! ---- am_perp
1284 
1285 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1286  thisroutine )
1287 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1288  thisroutine )
1289 call check( nf90_def_var(ncid, 'am_perp', nf90_float, nc2d, ncv), &
1290  thisroutine )
1291 buffer = 'm a-1'
1292 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1293  thisroutine )
1294 buffer = 'land_ice_volume_flux_across_zm_interface'
1295 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1296  thisroutine )
1297 buffer = 'Volume flux across the z=zm interface'
1298 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1299  thisroutine )
1300 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1301  thisroutine )
1302 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1303  thisroutine )
1304 
1305 ! ---- qx
1306 
1307 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1308  thisroutine )
1309 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1310  thisroutine )
1311 call check( nf90_def_var(ncid, 'qx', nf90_float, nc2d, ncv), &
1312  thisroutine )
1313 buffer = 'm2 a-1'
1314 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1315  thisroutine )
1316 buffer = 'land_ice_vertical_integral_x_velocity'
1317 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1318  thisroutine )
1319 buffer = 'Horizontal volume flux qx'
1320 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1321  thisroutine )
1322 buffer = 'Staggered grid variable, defined at (j,i+1/2)'
1323 call check( nf90_put_att(ncid, ncv, 'comment', trim(buffer)), &
1324  thisroutine )
1325 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1326  thisroutine )
1327 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1328  thisroutine )
1329 
1330 ! ---- qy
1331 
1332 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1333  thisroutine )
1334 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1335  thisroutine )
1336 call check( nf90_def_var(ncid, 'qy', nf90_float, nc2d, ncv), &
1337  thisroutine )
1338 buffer = 'm2 a-1'
1339 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1340  thisroutine )
1341 buffer = 'land_ice_vertical_integral_y_velocity'
1342 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1343  thisroutine )
1344 buffer = 'Horizontal volume flux qy'
1345 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1346  thisroutine )
1347 buffer = 'Staggered grid variable, defined at (j+1/2,i)'
1348 call check( nf90_put_att(ncid, ncv, 'comment', trim(buffer)), &
1349  thisroutine )
1350 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1351  thisroutine )
1352 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1353  thisroutine )
1354 
1355 ! ---- dzs_dt (= dzs_dtau)
1356 
1357 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1358  thisroutine )
1359 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1360  thisroutine )
1361 call check( nf90_def_var(ncid, 'dzs_dt', nf90_float, nc2d, ncv), &
1362  thisroutine )
1363 buffer = 'm a-1'
1364 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1365  thisroutine )
1366 buffer = 'tendency_of_surface_altitude'
1367 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1368  thisroutine )
1369 buffer = 'Rate of change of the topography of the free surface'
1370 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1371  thisroutine )
1372 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1373  thisroutine )
1374 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1375  thisroutine )
1376 
1377 ! ---- dzm_dt (= dzm_dtau)
1378 
1379 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1380  thisroutine )
1381 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1382  thisroutine )
1383 call check( nf90_def_var(ncid, 'dzm_dt', nf90_float, nc2d, ncv), &
1384  thisroutine )
1385 buffer = 'm a-1'
1386 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1387  thisroutine )
1388 buffer = 'tendency_of_zm_interface_altitude'
1389 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1390  thisroutine )
1391 buffer = 'Rate of change of the topography of the z=zm interface'
1392 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1393  thisroutine )
1394 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1395  thisroutine )
1396 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1397  thisroutine )
1398 
1399 ! ---- dzb_dt (= dzb_dtau)
1400 
1401 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1402  thisroutine )
1403 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1404  thisroutine )
1405 call check( nf90_def_var(ncid, 'dzb_dt', nf90_float, nc2d, ncv), &
1406  thisroutine )
1407 buffer = 'm a-1'
1408 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1409  thisroutine )
1410 buffer = 'tendency_of_ice_base_altitude'
1411 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1412  thisroutine )
1413 buffer = 'Rate of change of the topography of the ice base'
1414 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1415  thisroutine )
1416 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1417  thisroutine )
1418 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1419  thisroutine )
1420 
1421 ! ---- dzl_dt (= dzl_dtau)
1422 
1423 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1424  thisroutine )
1425 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1426  thisroutine )
1427 call check( nf90_def_var(ncid, 'dzl_dt', nf90_float, nc2d, ncv), &
1428  thisroutine )
1429 buffer = 'm a-1'
1430 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1431  thisroutine )
1432 buffer = 'tendency_of_bedrock_altitude'
1433 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1434  thisroutine )
1435 buffer = 'Rate of change of the topography of the lithosphere surface'
1436 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1437  thisroutine )
1438 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1439  thisroutine )
1440 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1441  thisroutine )
1442 
1443 ! ---- dH_c_dt (= dH_c_dtau)
1444 
1445 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1446  thisroutine )
1447 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1448  thisroutine )
1449 call check( nf90_def_var(ncid, 'dH_c_dt', nf90_float, nc2d, ncv), &
1450  thisroutine )
1451 buffer = 'm a-1'
1452 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1453  thisroutine )
1454 buffer = 'tendency_of_land_ice_kc_layer_thickness'
1455 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1456  thisroutine )
1457 buffer = 'Rate of change of the thickness of the upper (kc) ice layer'
1458 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1459  thisroutine )
1460 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1461  thisroutine )
1462 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1463  thisroutine )
1464 
1465 ! ---- dH_t_dt (= dH_t_dtau)
1466 
1467 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1468  thisroutine )
1469 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1470  thisroutine )
1471 call check( nf90_def_var(ncid, 'dH_t_dt', nf90_float, nc2d, ncv), &
1472  thisroutine )
1473 buffer = 'm a-1'
1474 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1475  thisroutine )
1476 buffer = 'tendency_of_land_ice_kt_layer_thickness'
1477 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1478  thisroutine )
1479 buffer = 'Rate of change of the thickness of the lower (kt) ice layer'
1480 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1481  thisroutine )
1482 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1483  thisroutine )
1484 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1485  thisroutine )
1486 
1487 ! ---- dH_dt (= dH_dtau)
1488 
1489 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1490  thisroutine )
1491 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1492  thisroutine )
1493 call check( nf90_def_var(ncid, 'dH_dt', nf90_float, nc2d, ncv), &
1494  thisroutine )
1495 buffer = 'm a-1'
1496 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1497  thisroutine )
1498 buffer = 'tendency_of_land_ice_thickness'
1499 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1500  thisroutine )
1501 buffer = 'Rate of change of the ice thickness'
1502 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1503  thisroutine )
1504 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1505  thisroutine )
1506 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1507  thisroutine )
1508 
1509 ! ---- vx_b_g
1510 
1511 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1512  thisroutine )
1513 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1514  thisroutine )
1515 call check( nf90_def_var(ncid, 'vx_b_g', nf90_float, nc2d, ncv), &
1516  thisroutine )
1517 buffer = 'm a-1'
1518 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1519  thisroutine )
1520 buffer = 'land_ice_base_x_velocity'
1521 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1522  thisroutine )
1523 buffer = 'Horizontal velocity vx at the ice base'
1524 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1525  thisroutine )
1526 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1527  thisroutine )
1528 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1529  thisroutine )
1530 
1531 ! ---- vy_b_g
1532 
1533 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1534  thisroutine )
1535 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1536  thisroutine )
1537 call check( nf90_def_var(ncid, 'vy_b_g', nf90_float, nc2d, ncv), &
1538  thisroutine )
1539 buffer = 'm a-1'
1540 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1541  thisroutine )
1542 buffer = 'land_ice_base_y_velocity'
1543 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1544  thisroutine )
1545 buffer = 'Horizontal velocity vy at the ice base'
1546 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1547  thisroutine )
1548 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1549  thisroutine )
1550 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1551  thisroutine )
1552 
1553 ! ---- vz_b
1554 
1555 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1556  thisroutine )
1557 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1558  thisroutine )
1559 call check( nf90_def_var(ncid, 'vz_b', nf90_float, nc2d, ncv), &
1560  thisroutine )
1561 buffer = 'm a-1'
1562 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1563  thisroutine )
1564 buffer = 'land_ice_base_z_velocity'
1565 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1566  thisroutine )
1567 buffer = 'Vertical velocity vz at the ice base'
1568 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1569  thisroutine )
1570 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1571  thisroutine )
1572 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1573  thisroutine )
1574 
1575 ! ---- vh_b
1576 
1577 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1578  thisroutine )
1579 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1580  thisroutine )
1581 call check( nf90_def_var(ncid, 'vh_b', nf90_float, nc2d, ncv), &
1582  thisroutine )
1583 buffer = 'm a-1'
1584 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1585  thisroutine )
1586 buffer = 'land_ice_base_horizontal_velocity'
1587 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1588  thisroutine )
1589 buffer = 'Horizontal velocity vh at the ice base'
1590 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1591  thisroutine )
1592 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1593  thisroutine )
1594 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1595  thisroutine )
1596 
1597 ! ---- vx_s_g
1598 
1599 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1600  thisroutine )
1601 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1602  thisroutine )
1603 call check( nf90_def_var(ncid, 'vx_s_g', nf90_float, nc2d, ncv), &
1604  thisroutine )
1605 buffer = 'm a-1'
1606 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1607  thisroutine )
1608 buffer = 'land_ice_surface_x_velocity'
1609 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1610  thisroutine )
1611 buffer = 'Horizontal velocity vx at the ice surface'
1612 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1613  thisroutine )
1614 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1615  thisroutine )
1616 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1617  thisroutine )
1618 
1619 ! ---- vy_s_g
1620 
1621 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1622  thisroutine )
1623 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1624  thisroutine )
1625 call check( nf90_def_var(ncid, 'vy_s_g', nf90_float, nc2d, ncv), &
1626  thisroutine )
1627 buffer = 'm a-1'
1628 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1629  thisroutine )
1630 buffer = 'land_ice_surface_y_velocity'
1631 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1632  thisroutine )
1633 buffer = 'Horizontal velocity vy at the ice surface'
1634 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1635  thisroutine )
1636 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1637  thisroutine )
1638 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1639  thisroutine )
1640 
1641 ! ---- vz_s
1642 
1643 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1644  thisroutine )
1645 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1646  thisroutine )
1647 call check( nf90_def_var(ncid, 'vz_s', nf90_float, nc2d, ncv), &
1648  thisroutine )
1649 buffer = 'm a-1'
1650 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1651  thisroutine )
1652 buffer = 'land_ice_surface_z_velocity'
1653 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1654  thisroutine )
1655 buffer = 'Vertical velocity vz at the ice surface'
1656 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1657  thisroutine )
1658 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1659  thisroutine )
1660 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1661  thisroutine )
1662 
1663 ! ---- vh_s
1664 
1665 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1666  thisroutine )
1667 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1668  thisroutine )
1669 call check( nf90_def_var(ncid, 'vh_s', nf90_float, nc2d, ncv), &
1670  thisroutine )
1671 buffer = 'm a-1'
1672 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1673  thisroutine )
1674 buffer = 'land_ice_surface_horizontal_velocity'
1675 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1676  thisroutine )
1677 buffer = 'Horizontal velocity vh at the ice surface'
1678 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1679  thisroutine )
1680 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1681  thisroutine )
1682 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1683  thisroutine )
1684 
1685 ! ---- vx_m_g
1686 
1687 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1688  thisroutine )
1689 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1690  thisroutine )
1691 call check( nf90_def_var(ncid, 'vx_m_g', nf90_float, nc2d, ncv), &
1692  thisroutine )
1693 buffer = 'm a-1'
1694 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1695  thisroutine )
1696 buffer = 'land_ice_vertical_mean_x_velocity'
1697 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1698  thisroutine )
1699 buffer = 'Vertical mean of horizontal velocity vx'
1700 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1701  thisroutine )
1702 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1703  thisroutine )
1704 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1705  thisroutine )
1706 
1707 ! ---- vy_m_g
1708 
1709 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1710  thisroutine )
1711 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1712  thisroutine )
1713 call check( nf90_def_var(ncid, 'vy_m_g', nf90_float, nc2d, ncv), &
1714  thisroutine )
1715 buffer = 'm a-1'
1716 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1717  thisroutine )
1718 buffer = 'land_ice_vertical_mean_y_velocity'
1719 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1720  thisroutine )
1721 buffer = 'Vertical mean of horizontal velocity vy'
1722 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1723  thisroutine )
1724 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1725  thisroutine )
1726 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1727  thisroutine )
1728 
1729 ! ---- vh_m
1730 
1731 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1732  thisroutine )
1733 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1734  thisroutine )
1735 call check( nf90_def_var(ncid, 'vh_m', nf90_float, nc2d, ncv), &
1736  thisroutine )
1737 buffer = 'm a-1'
1738 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1739  thisroutine )
1740 buffer = 'land_ice_vertical_mean_horizontal_velocity'
1741 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1742  thisroutine )
1743 buffer = 'Vertical mean of horizontal velocity vh'
1744 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1745  thisroutine )
1746 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1747  thisroutine )
1748 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1749  thisroutine )
1750 
1751 ! ---- temp_b
1752 
1753 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1754  thisroutine )
1755 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1756  thisroutine )
1757 call check( nf90_def_var(ncid, 'temp_b', nf90_float, nc2d, ncv), &
1758  thisroutine )
1759 buffer = 'degC'
1760 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1761  thisroutine )
1762 buffer = 'basal_temperature'
1763 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1764  thisroutine )
1765 buffer = 'Temperature at the ice base'
1766 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1767  thisroutine )
1768 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1769  thisroutine )
1770 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1771  thisroutine )
1772 
1773 ! ---- temph_b
1774 
1775 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1776  thisroutine )
1777 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1778  thisroutine )
1779 call check( nf90_def_var(ncid, 'temph_b', nf90_float, nc2d, ncv), &
1780  thisroutine )
1781 buffer = 'degC'
1782 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1783  thisroutine )
1784 buffer = 'basal_temperature_rel_to_pmp'
1785 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1786  thisroutine )
1787 buffer = 'Temperature at the ice base relative to the pressure melting point'
1788 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1789  thisroutine )
1790 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1791  thisroutine )
1792 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1793  thisroutine )
1794 
1795 ! ---- tau_b_driving
1796 
1797 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1798  thisroutine )
1799 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1800  thisroutine )
1801 call check( nf90_def_var(ncid, 'tau_b_driving', nf90_float, nc2d, ncv), &
1802  thisroutine )
1803 buffer = 'Pa'
1804 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1805  thisroutine )
1806 buffer = 'magnitude_of_land_ice_driving_stress'
1807 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1808  thisroutine )
1809 buffer = 'Driving stress'
1810 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1811  thisroutine )
1812 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1813  thisroutine )
1814 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1815  thisroutine )
1816 
1817 ! ---- tau_b_drag
1818 
1819 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1820  thisroutine )
1821 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1822  thisroutine )
1823 call check( nf90_def_var(ncid, 'tau_b_drag', nf90_float, nc2d, ncv), &
1824  thisroutine )
1825 buffer = 'Pa'
1826 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1827  thisroutine )
1828 buffer = 'magnitude_of_land_ice_basal_drag'
1829 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1830  thisroutine )
1831 buffer = 'Basal drag'
1832 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1833  thisroutine )
1834 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1835  thisroutine )
1836 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1837  thisroutine )
1838 
1839 ! ---- p_b_w
1840 
1841 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1842  thisroutine )
1843 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1844  thisroutine )
1845 call check( nf90_def_var(ncid, 'p_b_w', nf90_float, nc2d, ncv), &
1846  thisroutine )
1847 buffer = 'Pa'
1848 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1849  thisroutine )
1850 buffer = 'basal_water_pressure'
1851 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1852  thisroutine )
1853 buffer = 'Basal water pressure'
1854 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1855  thisroutine )
1856 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1857  thisroutine )
1858 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1859  thisroutine )
1860 
1861 ! ---- H_w
1862 
1863 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1864  thisroutine )
1865 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1866  thisroutine )
1867 call check( nf90_def_var(ncid, 'H_w', nf90_float, nc2d, ncv), &
1868  thisroutine )
1869 buffer = 'm'
1870 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1871  thisroutine )
1872 buffer = 'water_column_thickness'
1873 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1874  thisroutine )
1875 buffer = 'Thickness of the water column under the ice base'
1876 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1877  thisroutine )
1878 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1879  thisroutine )
1880 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1881  thisroutine )
1882 
1883 ! ---- q_gl_g
1884 
1885 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1886  thisroutine )
1887 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1888  thisroutine )
1889 call check( nf90_def_var(ncid, 'q_gl_g', nf90_float, nc2d, ncv), &
1890  thisroutine )
1891 buffer = 'm2 a-1'
1892 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1893  thisroutine )
1894 buffer = 'land_ice_volume_flux_across_gl'
1895 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1896  thisroutine )
1897 buffer = 'Horizontal volume flux across the grounding line'
1898 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1899  thisroutine )
1900 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1901  thisroutine )
1902 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1903  thisroutine )
1904 
1905 ! ---- q_cf_g
1906 
1907 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1908  thisroutine )
1909 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1910  thisroutine )
1911 call check( nf90_def_var(ncid, 'q_cf_g', nf90_float, nc2d, ncv), &
1912  thisroutine )
1913 buffer = 'm a-1'
1914 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1915  thisroutine )
1916 buffer = 'land_ice_volume_flux_due_to_calving'
1917 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1918  thisroutine )
1919 buffer = 'Calving flux'
1920 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1921  thisroutine )
1922 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1923  thisroutine )
1924 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1925  thisroutine )
1926 
1927 ! ---- ratio_sl_x
1928 
1929 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1930  thisroutine )
1931 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1932  thisroutine )
1933 call check( nf90_def_var(ncid, 'ratio_sl_x', nf90_float, nc2d, ncv), &
1934  thisroutine )
1935 buffer = '-'
1936 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1937  thisroutine )
1938 buffer = 'land_ice_x_slip_ratio'
1939 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1940  thisroutine )
1941 buffer = 'Ratio of basal to surface velocity (slip ratio) in x-direction, ' &
1942  // 'at (i+1/2,j)'
1943 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1944  thisroutine )
1945 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1946  thisroutine )
1947 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1948  thisroutine )
1949 
1950 ! ---- ratio_sl_y
1951 
1952 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1953  thisroutine )
1954 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1955  thisroutine )
1956 call check( nf90_def_var(ncid, 'ratio_sl_y', nf90_float, nc2d, ncv), &
1957  thisroutine )
1958 buffer = '-'
1959 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1960  thisroutine )
1961 buffer = 'land_ice_y_slip_ratio'
1962 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1963  thisroutine )
1964 buffer = 'Ratio of basal to surface velocity (slip ratio) in y-direction, ' &
1965  // 'at (i+1/2,j)'
1966 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1967  thisroutine )
1968 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1969  thisroutine )
1970 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1971  thisroutine )
1972 
1973 ! ---- flag_shelfy_stream_x
1974 
1975 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1976  thisroutine )
1977 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
1978  thisroutine )
1979 call check( nf90_def_var(ncid, 'flag_shelfy_stream_x', nf90_byte, nc2d, ncv), &
1980  thisroutine )
1981 buffer = '-'
1982 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1983  thisroutine )
1984 buffer = 'land_ice_x_shelfy_stream_flag'
1985 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1986  thisroutine )
1987 buffer = 'Shelfy stream flag in x-direction, at (i+1/2,j)'
1988 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1989  thisroutine )
1990 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1991  thisroutine )
1992 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1993  thisroutine )
1994 
1995 ! ---- flag_shelfy_stream_y
1996 
1997 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
1998  thisroutine )
1999 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
2000  thisroutine )
2001 call check( nf90_def_var(ncid, 'flag_shelfy_stream_y', nf90_byte, nc2d, ncv), &
2002  thisroutine )
2003 buffer = '-'
2004 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
2005  thisroutine )
2006 buffer = 'land_ice_y_shelfy_stream_flag'
2007 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
2008  thisroutine )
2009 buffer = 'Shelfy stream flag in y-direction, at (i,j+1/2)'
2010 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
2011  thisroutine )
2012 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
2013  thisroutine )
2014 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
2015  thisroutine )
2016 
2017 ! ---- flag_shelfy_stream
2018 
2019 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
2020  thisroutine )
2021 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
2022  thisroutine )
2023 call check( nf90_def_var(ncid, 'flag_shelfy_stream', nf90_byte, nc2d, ncv), &
2024  thisroutine )
2025 buffer = '-'
2026 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
2027  thisroutine )
2028 buffer = 'land_ice_shelfy_stream_flag'
2029 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
2030  thisroutine )
2031 buffer = 'Shelfy stream flag'
2032 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
2033  thisroutine )
2034 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
2035  thisroutine )
2036 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
2037  thisroutine )
2038 
2039 ! ---- vis_int_g
2040 
2041 call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc2d(1)), &
2042  thisroutine )
2043 call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc2d(2)), &
2044  thisroutine )
2045 call check( nf90_def_var(ncid, 'vis_int_g', nf90_float, nc2d, ncv), &
2046  thisroutine )
2047 buffer = 'Pa s m'
2048 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
2049  thisroutine )
2050 buffer = 'land_ice_vertical_integral_viscosity'
2051 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
2052  thisroutine )
2053 buffer = 'Depth-integrated viscosity'
2054 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
2055  thisroutine )
2056 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
2057  thisroutine )
2058 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
2059  thisroutine )
2060 
2061 if (flag_3d_output) then
2062 
2063 ! ---- vx_c
2064 
2065  call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc3d(1)), &
2066  thisroutine )
2067  call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc3d(2)), &
2068  thisroutine )
2069  call check( nf90_inq_dimid(ncid, trim(coord_id(3)), nc3d(3)), &
2070  thisroutine )
2071  call check( nf90_def_var(ncid, 'vx_c', nf90_float, nc3d, ncv), &
2072  thisroutine )
2073  buffer = 'm a-1'
2074  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
2075  thisroutine )
2076  buffer = 'land_ice_kc_layer_x_velocity'
2077  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
2078  thisroutine )
2079  buffer = 'Horizontal velocity vx in the upper (kc) ice layer'
2080  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
2081  thisroutine )
2082  buffer = 'Staggered grid variable, defined at (kc,j,i+1/2)'
2083  call check( nf90_put_att(ncid, ncv, 'comment', trim(buffer)), &
2084  thisroutine )
2085  call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
2086  thisroutine )
2087  call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
2088  thisroutine )
2089 
2090 ! ---- vy_c
2091 
2092  call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc3d(1)), &
2093  thisroutine )
2094  call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc3d(2)), &
2095  thisroutine )
2096  call check( nf90_inq_dimid(ncid, trim(coord_id(3)), nc3d(3)), &
2097  thisroutine )
2098  call check( nf90_def_var(ncid, 'vy_c', nf90_float, nc3d, ncv), &
2099  thisroutine )
2100  buffer = 'm a-1'
2101  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
2102  thisroutine )
2103  buffer = 'land_ice_kc_layer_y_velocity'
2104  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
2105  thisroutine )
2106  buffer = 'Horizontal velocity vy in the upper (kc) ice layer'
2107  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
2108  thisroutine )
2109  buffer = 'Staggered grid variable, defined at (kc,j+1/2,i)'
2110  call check( nf90_put_att(ncid, ncv, 'comment', trim(buffer)), &
2111  thisroutine )
2112  call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
2113  thisroutine )
2114  call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
2115  thisroutine )
2116 
2117 ! ---- vz_c
2118 
2119  call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc3d(1)), &
2120  thisroutine )
2121  call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc3d(2)), &
2122  thisroutine )
2123  call check( nf90_inq_dimid(ncid, trim(coord_id(3)), nc3d(3)), &
2124  thisroutine )
2125  call check( nf90_def_var(ncid, 'vz_c', nf90_float, nc3d, ncv), &
2126  thisroutine )
2127  buffer = 'm a-1'
2128  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
2129  thisroutine )
2130  buffer = 'land_ice_kc_layer_z_velocity'
2131  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
2132  thisroutine )
2133  buffer = 'Vertical velocity vz in the upper (kc) ice layer'
2134  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
2135  thisroutine )
2136  buffer = 'Staggered grid variable, defined at (kc+1/2,j,i)'
2137  call check( nf90_put_att(ncid, ncv, 'comment', trim(buffer)), &
2138  thisroutine )
2139  call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
2140  thisroutine )
2141  call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
2142  thisroutine )
2143 
2144 ! ---- vx_t
2145 
2146  call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc3d(1)), &
2147  thisroutine )
2148  call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc3d(2)), &
2149  thisroutine )
2150  call check( nf90_inq_dimid(ncid, trim(coord_id(4)), nc3d(3)), &
2151  thisroutine )
2152  call check( nf90_def_var(ncid, 'vx_t', nf90_float, nc3d, ncv), &
2153  thisroutine )
2154  buffer = 'm a-1'
2155  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
2156  thisroutine )
2157  buffer = 'land_ice_kt_layer_x_velocity'
2158  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
2159  thisroutine )
2160  buffer = 'Horizontal velocity vx in the lower (kt) ice layer'
2161  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
2162  thisroutine )
2163  buffer = 'Staggered grid variable, defined at (kt,j,i+1/2)'
2164  call check( nf90_put_att(ncid, ncv, 'comment', trim(buffer)), &
2165  thisroutine )
2166  call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
2167  thisroutine )
2168  call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
2169  thisroutine )
2170 
2171 ! ---- vy_t
2172 
2173  call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc3d(1)), &
2174  thisroutine )
2175  call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc3d(2)), &
2176  thisroutine )
2177  call check( nf90_inq_dimid(ncid, trim(coord_id(4)), nc3d(3)), &
2178  thisroutine )
2179  call check( nf90_def_var(ncid, 'vy_t', nf90_float, nc3d, ncv), &
2180  thisroutine )
2181  buffer = 'm a-1'
2182  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
2183  thisroutine )
2184  buffer = 'land_ice_kt_layer_y_velocity'
2185  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
2186  thisroutine )
2187  buffer = 'Horizontal velocity vy in the lower (kt) ice layer'
2188  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
2189  thisroutine )
2190  buffer = 'Staggered grid variable, defined at (kt,j+1/2,i)'
2191  call check( nf90_put_att(ncid, ncv, 'comment', trim(buffer)), &
2192  thisroutine )
2193  call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
2194  thisroutine )
2195  call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
2196  thisroutine )
2197 
2198 ! ---- vz_t
2199 
2200  call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc3d(1)), &
2201  thisroutine )
2202  call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc3d(2)), &
2203  thisroutine )
2204  call check( nf90_inq_dimid(ncid, trim(coord_id(4)), nc3d(3)), &
2205  thisroutine )
2206  call check( nf90_def_var(ncid, 'vz_t', nf90_float, nc3d, ncv), &
2207  thisroutine )
2208  buffer = 'm a-1'
2209  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
2210  thisroutine )
2211  buffer = 'land_ice_kt_layer_z_velocity'
2212  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
2213  thisroutine )
2214  buffer = 'Vertical velocity vz in the lower (kt) ice layer'
2215  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
2216  thisroutine )
2217  buffer = 'Staggered grid variable, defined at (kt+1/2,j,i)'
2218  call check( nf90_put_att(ncid, ncv, 'comment', trim(buffer)), &
2219  thisroutine )
2220  call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
2221  thisroutine )
2222  call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
2223  thisroutine )
2224 
2225 ! ---- temp_c
2226 
2227  call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc3d(1)), &
2228  thisroutine )
2229  call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc3d(2)), &
2230  thisroutine )
2231  call check( nf90_inq_dimid(ncid, trim(coord_id(3)), nc3d(3)), &
2232  thisroutine )
2233  call check( nf90_def_var(ncid, 'temp_c', nf90_float, nc3d, ncv), &
2234  thisroutine )
2235  buffer = 'degC'
2236  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
2237  thisroutine )
2238  buffer = 'land_ice_kc_layer_temperature'
2239  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
2240  thisroutine )
2241  buffer = 'Temperature in the upper (kc) ice layer'
2242  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
2243  thisroutine )
2244  call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
2245  thisroutine )
2246  call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
2247  thisroutine )
2248 
2249 ! ---- omega_t
2250 
2251  call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc3d(1)), &
2252  thisroutine )
2253  call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc3d(2)), &
2254  thisroutine )
2255  call check( nf90_inq_dimid(ncid, trim(coord_id(4)), nc3d(3)), &
2256  thisroutine )
2257  call check( nf90_def_var(ncid, 'omega_t', nf90_float, nc3d, ncv), &
2258  thisroutine )
2259  buffer = '1'
2260  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
2261  thisroutine )
2262  buffer = 'land_ice_kt_layer_water_content'
2263  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
2264  thisroutine )
2265  buffer = 'Water content in the lower (kt) ice layer'
2266  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
2267  thisroutine )
2268  call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
2269  thisroutine )
2270  call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
2271  thisroutine )
2272 
2273 ! ---- temp_r
2274 
2275  call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc3d(1)), &
2276  thisroutine )
2277  call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc3d(2)), &
2278  thisroutine )
2279  call check( nf90_inq_dimid(ncid, trim(coord_id(5)), nc3d(3)), &
2280  thisroutine )
2281  call check( nf90_def_var(ncid, 'temp_r', nf90_float, nc3d, ncv), &
2282  thisroutine )
2283  buffer = 'degC'
2284  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
2285  thisroutine )
2286  buffer = 'lithosphere_layer_temperature'
2287  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
2288  thisroutine )
2289  buffer = 'Temperature in the lithosphere layer'
2290  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
2291  thisroutine )
2292  call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
2293  thisroutine )
2294  call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
2295  thisroutine )
2296 
2297 ! ---- enth_c
2298 
2299  call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc3d(1)), &
2300  thisroutine )
2301  call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc3d(2)), &
2302  thisroutine )
2303  call check( nf90_inq_dimid(ncid, trim(coord_id(3)), nc3d(3)), &
2304  thisroutine )
2305  call check( nf90_def_var(ncid, 'enth_c', nf90_float, nc3d, ncv), &
2306  thisroutine )
2307  buffer = 'J kg-1'
2308  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
2309  thisroutine )
2310  buffer = 'land_ice_kc_layer_enthalpy'
2311  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
2312  thisroutine )
2313  buffer = 'Enthalpy in the upper (kc) ice layer'
2314  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
2315  thisroutine )
2316  call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
2317  thisroutine )
2318  call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
2319  thisroutine )
2320 
2321 ! ---- enth_t
2322 
2323  call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc3d(1)), &
2324  thisroutine )
2325  call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc3d(2)), &
2326  thisroutine )
2327  call check( nf90_inq_dimid(ncid, trim(coord_id(4)), nc3d(3)), &
2328  thisroutine )
2329  call check( nf90_def_var(ncid, 'enth_t', nf90_float, nc3d, ncv), &
2330  thisroutine )
2331  buffer = 'J kg-1'
2332  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
2333  thisroutine )
2334  buffer = 'land_ice_kt_layer_enthalpy'
2335  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
2336  thisroutine )
2337  buffer = 'Enthalpy in the lower (kt) ice layer'
2338  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
2339  thisroutine )
2340  call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
2341  thisroutine )
2342  call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
2343  thisroutine )
2344 
2345 ! ---- omega_c
2346 
2347  call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc3d(1)), &
2348  thisroutine )
2349  call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc3d(2)), &
2350  thisroutine )
2351  call check( nf90_inq_dimid(ncid, trim(coord_id(3)), nc3d(3)), &
2352  thisroutine )
2353  call check( nf90_def_var(ncid, 'omega_c', nf90_float, nc3d, ncv), &
2354  thisroutine )
2355  buffer = '1'
2356  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
2357  thisroutine )
2358  buffer = 'land_ice_kc_layer_water_content'
2359  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
2360  thisroutine )
2361  buffer = 'Water content in the upper (kc) ice layer'
2362  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
2363  thisroutine )
2364  call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
2365  thisroutine )
2366  call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
2367  thisroutine )
2368 
2369 ! ---- enh_c
2370 
2371  call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc3d(1)), &
2372  thisroutine )
2373  call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc3d(2)), &
2374  thisroutine )
2375  call check( nf90_inq_dimid(ncid, trim(coord_id(3)), nc3d(3)), &
2376  thisroutine )
2377  call check( nf90_def_var(ncid, 'enh_c', nf90_float, nc3d, ncv), &
2378  thisroutine )
2379  buffer = '1'
2380  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
2381  thisroutine )
2382  buffer = 'land_ice_kc_layer_flow_enhancement_factor'
2383  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
2384  thisroutine )
2385  buffer = 'Flow enhancement factor in the upper (kc) ice layer'
2386  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
2387  thisroutine )
2388  call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
2389  thisroutine )
2390  call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
2391  thisroutine )
2392 
2393 ! ---- enh_t
2394 
2395  call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc3d(1)), &
2396  thisroutine )
2397  call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc3d(2)), &
2398  thisroutine )
2399  call check( nf90_inq_dimid(ncid, trim(coord_id(4)), nc3d(3)), &
2400  thisroutine )
2401  call check( nf90_def_var(ncid, 'enh_t', nf90_float, nc3d, ncv), &
2402  thisroutine )
2403  buffer = '1'
2404  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
2405  thisroutine )
2406  buffer = 'land_ice_kt_layer_flow_enhancement_factor'
2407  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
2408  thisroutine )
2409  buffer = 'Flow enhancement factor in the lower (kt) ice layer'
2410  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
2411  thisroutine )
2412  call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
2413  thisroutine )
2414  call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
2415  thisroutine )
2416 
2417 ! ---- age_c
2418 
2419  call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc3d(1)), &
2420  thisroutine )
2421  call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc3d(2)), &
2422  thisroutine )
2423  call check( nf90_inq_dimid(ncid, trim(coord_id(3)), nc3d(3)), &
2424  thisroutine )
2425  call check( nf90_def_var(ncid, 'age_c', nf90_float, nc3d, ncv), &
2426  thisroutine )
2427  buffer = 'a'
2428  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
2429  thisroutine )
2430  buffer = 'land_ice_kc_layer_age'
2431  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
2432  thisroutine )
2433  buffer = 'Age in the upper (kc) ice layer'
2434  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
2435  thisroutine )
2436  call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
2437  thisroutine )
2438  call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
2439  thisroutine )
2440 
2441 ! ---- age_t
2442 
2443  call check( nf90_inq_dimid(ncid, trim(coord_id(1)), nc3d(1)), &
2444  thisroutine )
2445  call check( nf90_inq_dimid(ncid, trim(coord_id(2)), nc3d(2)), &
2446  thisroutine )
2447  call check( nf90_inq_dimid(ncid, trim(coord_id(4)), nc3d(3)), &
2448  thisroutine )
2449  call check( nf90_def_var(ncid, 'age_t', nf90_float, nc3d, ncv), &
2450  thisroutine )
2451  buffer = 'a'
2452  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
2453  thisroutine )
2454  buffer = 'land_ice_kt_layer_age'
2455  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
2456  thisroutine )
2457  buffer = 'Age in the lower (kt) ice layer'
2458  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
2459  thisroutine )
2460  call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
2461  thisroutine )
2462  call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
2463  thisroutine )
2464 
2465 end if
2466 
2467 ! ---- End of the definitions
2468 
2469 call check( nf90_enddef(ncid), thisroutine )
2470 
2471 #endif
2472 
2473 !-------- Ice thickness and time derivative --------
2474 
2475 h = h_c + h_t
2476 dh_dtau = dh_c_dtau + dh_t_dtau
2477 
2478 !-------- Thickness of the cold and temperate layers --------
2479 
2480 h_cold = 0.0_dp
2481 h_temp = 0.0_dp
2482 
2483 #if (CALCMOD==1)
2484 do i=1, imax-1
2485 do j=1, jmax-1
2486  h_temp(j,i) = h_t(j,i)
2487 end do
2488 end do
2489 #elif (CALCMOD==0 || CALCMOD==2 || CALCMOD==3 || CALCMOD==-1)
2490 do i=1, imax-1
2491 do j=1, jmax-1
2492  h_temp(j,i) = h_c(j,i)*eaz_c_quotient(kc_cts(j,i))
2493 end do
2494 end do
2495 #endif
2496 
2497 h_cold = h - h_temp
2498 
2499 !-------- Enthalpies for the non-enthalpy methods (POLY, COLD, ISOT) --------
2500 
2501 #if (CALCMOD==1)
2502 
2503 do i=0, imax
2504 do j=0, jmax
2505 
2506  do kc=0, kcmax
2507  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
2508  end do
2509 
2510  if ( (maske(j,i)==0_i2b).and.(n_cts(j,i)==1_i2b) ) then
2511  do kt=0, ktmax
2512  enth_t(kt,j,i) = enth_fct_temp_omega(temp_t_m(kt,j,i), omega_t(kt,j,i))
2513  end do
2514  else
2515  do kt=0, ktmax
2516  enth_t(kt,j,i) = enth_c(0,j,i)
2517  end do
2518  end if
2519 
2520 end do
2521 end do
2522 
2523 #elif (CALCMOD==0 || CALCMOD==-1)
2524 
2525 do i=0, imax
2526 do j=0, jmax
2527 
2528  do kc=0, kcmax
2529  enth_c(kc,j,i) = enth_fct_temp_omega(temp_c(kc,j,i), 0.0_dp)
2530  end do
2531 
2532  do kt=0, ktmax
2533  enth_t(kt,j,i) = enth_c(0,j,i)
2534  end do
2535 
2536 end do
2537 end do
2538 
2539 #endif
2540 
2541 !-------- Vertical mean of horizontal velocities --------
2542 
2543 vx_m_g = 0.0_dp
2544 vy_m_g = 0.0_dp
2545 
2546 do i=1, imax-1
2547 do j=1, jmax-1
2548 
2549  if ( (maske(j,i)==0_i2b).or.(maske(j,i)==3_i2b) ) then
2550  vx_m_g(j,i) = 0.5_dp*(vx_m(j,i)+vx_m(j,i-1))
2551  vy_m_g(j,i) = 0.5_dp*(vy_m(j,i)+vy_m(j-1,i))
2552  end if
2553 
2554 end do
2555 end do
2556 
2557 !-------- Driving stress and basal drag --------
2558 
2559 tau_b_driving = 0.0_dp
2560 tau_b_drag = 0.0_dp
2561 
2562 do i=0, imax
2563 do j=0, jmax
2564 
2565  if (maske(j,i)==0_i2b) then ! grounded ice
2566 
2567  tau_b_driving(j,i) = rho*g*h(j,i) &
2568  * sqrt( dzs_dxi_g(j,i)**2 + dzs_deta_g(j,i)**2 )
2569 
2570  if (.not.flag_shelfy_stream(j,i)) then
2571  tau_b_drag(j,i) = tau_b_driving(j,i)
2572  else
2573  tau_b_drag(j,i) = no_value_neg_2 ! dummy value
2574  end if
2575 
2576  else if (maske(j,i)==3_i2b) then ! floating ice
2577 
2578  tau_b_driving(j,i) = rho*g*h(j,i) &
2579  * sqrt( dzs_dxi_g(j,i)**2 + dzs_deta_g(j,i)**2 )
2580 
2581  tau_b_drag(j,i) = 0.0_dp
2582 
2583  end if
2584 
2585 end do
2586 end do
2587 
2588 !-------- Computation of the ice volume, the volume above flotation,
2589 ! the area covered by grounded ice, the area covered by floating ice,
2590 ! the total surface mass balance, the total basal mass balance
2591 ! and the total calving flux --------
2592 
2593 rhosw_rho_ratio = rho_sw/rho
2594 
2595 v_grounded = 0.0_dp
2596 v_floating = 0.0_dp
2597 v_gr_redu = 0.0_dp
2598 a_grounded = 0.0_dp
2599 a_floating = 0.0_dp
2600 
2601 do i=0, imax
2602 do j=0, jmax
2603 
2604  if (maske(j,i)==0_i2b) then ! grounded ice
2605 
2606  v_grounded = v_grounded + h(j,i)*area(j,i)
2607  a_grounded = a_grounded + area(j,i)
2608  v_gr_redu = v_gr_redu &
2609  + rhosw_rho_ratio*max((z_sl-zl(j,i)),0.0_dp)*area(j,i)
2610 
2611  else if (maske(j,i)==3_i2b) then ! floating ice
2612 
2613  v_floating = v_floating + h(j,i)*area(j,i)
2614  a_floating = a_floating + area(j,i)
2615 
2616  end if
2617 
2618 end do
2619 end do
2620 
2621 v_tot = v_grounded + v_floating
2622 v_af = v_grounded - v_gr_redu
2623 
2624 !-------- Convert data to real*4 and seconds to years --------
2625 
2626 #if (!defined(OUT_TIMES) || OUT_TIMES==1)
2627 time_conv = real(time*sec_to_year,sp)
2628 #elif (OUT_TIMES==2)
2629 time_conv = real((time+year_zero)*sec_to_year,sp)
2630 #else
2631 stop ' >>> output1: OUT_TIMES must be either 1 or 2!'
2632 #endif
2633 
2634 delta_ts_conv = real(delta_ts,sp)
2635 glac_index_conv = real(glac_index,sp)
2636 z_sl_conv = real(z_sl,sp)
2637 v_tot_conv = real(v_tot,sp)
2638 v_af_conv = real(v_af,sp)
2639 a_grounded_conv = real(a_grounded,sp)
2640 a_floating_conv = real(a_floating,sp)
2641 h_r_conv = real(h_r,sp)
2642 
2643 do i=0, imax
2644  xi_conv(i) = real(xi(i),sp)
2645 end do
2646 
2647 do j=0, jmax
2648  eta_conv(j) = real(eta(j),sp)
2649 end do
2650 
2651 do kc=0, kcmax
2652  sigma_level_c_conv(kc) = real(eaz_c_quotient(kc),sp)
2653 end do
2654 
2655 do kt=0, ktmax
2656  sigma_level_t_conv(kt) = real(zeta_t(kt),sp)
2657 end do
2658 
2659 do kr=0, krmax
2660  sigma_level_r_conv(kr) = real(kr,sp)/real(krmax,sp)
2661 end do
2662 
2663 do i=0, imax
2664 do j=0, jmax
2665 
2666  maske_conv(i,j) = maske(j,i)
2667  maske_old_conv(i,j) = maske_old(j,i)
2668  mask_run_conv(i,j) = mask_run(j,i)
2669  n_cts_conv(i,j) = n_cts(j,i)
2670  kc_cts_conv(i,j) = kc_cts(j,i)
2671 
2672  lambda_conv(i,j) = real(lambda(j,i),sp)
2673  phi_conv(i,j) = real(phi(j,i),sp)
2674  lon_conv(i,j) = real(lambda(j,i)*pi_180_inv,sp) ! longitude in deg
2675  lon_conv(i,j) = modulo(lon_conv(i,j)+180.0_sp, 360.0_sp)-180.0_sp
2676  ! mapping to interval [-180 deg, +180 deg)
2677  lat_conv(i,j) = real(phi(j,i) *pi_180_inv,sp) ! latitute in deg
2678  if (lat_conv(i,j) > 90.0_sp) lat_conv(i,j) = 90.0_sp
2679  if (lat_conv(i,j) < -90.0_sp) lat_conv(i,j) = -90.0_sp
2680  ! constraining to interval [-90 deg, +90 deg]
2681  temp_s_conv(i,j) = real(temp_s(j,i),sp)
2682  accum_conv(i,j) = real(accum(j,i)*year_sec,sp)
2683  as_perp_conv(i,j) = real(as_perp(j,i)*year_sec,sp)
2684  as_perp_apl_conv(i,j) = real(as_perp_apl(j,i)*year_sec,sp)
2685  smb_corr_conv(i,j) = real(smb_corr(j,i)*year_sec,sp)
2686 
2687 #if (DISC>0) /* Ice discharge parameterisation */
2688  dis_perp_conv(i,j) = real(dis_perp(j,i)*year_sec,sp)
2689  cst_dist_conv(i,j) = real(cst_dist(j,i)*0.001_dp,sp)
2690  cos_grad_tc_conv(i,j) = real(cos_grad_tc(j,i),sp)
2691  mask_mar_conv(i,j) = mask_mar(j,i)
2692 #endif
2693 
2694  q_geo_conv(i,j) = real(q_geo(j,i),sp)
2695  zs_conv(i,j) = real(zs(j,i),sp)
2696  zm_conv(i,j) = real(zm(j,i),sp)
2697  zb_conv(i,j) = real(zb(j,i),sp)
2698  zl_conv(i,j) = real(zl(j,i),sp)
2699  zl0_conv(i,j) = real(zl0(j,i),sp)
2700  h_cold_conv(i,j) = real(H_cold(j,i),sp)
2701  h_temp_conv(i,j) = real(H_temp(j,i),sp)
2702  h_conv(i,j) = real(H(j,i),sp)
2703  q_bm_conv(i,j) = real(q_bm(j,i)*year_sec,sp)
2704  q_tld_conv(i,j) = real(q_tld(j,i)*year_sec,sp)
2705  am_perp_conv(i,j) = real(am_perp(j,i)*year_sec,sp)
2706  qx_conv(i,j) = real(qx(j,i)*year_sec,sp)
2707  qy_conv(i,j) = real(qy(j,i)*year_sec,sp)
2708  dzs_dtau_conv(i,j) = real(dzs_dtau(j,i)*year_sec,sp)
2709  dzm_dtau_conv(i,j) = real(dzm_dtau(j,i)*year_sec,sp)
2710  dzb_dtau_conv(i,j) = real(dzb_dtau(j,i)*year_sec,sp)
2711  dzl_dtau_conv(i,j) = real(dzl_dtau(j,i)*year_sec,sp)
2712  dh_c_dtau_conv(i,j) = real(dh_c_dtau(j,i)*year_sec,sp)
2713  dh_t_dtau_conv(i,j) = real(dh_t_dtau(j,i)*year_sec,sp)
2714  dh_dtau_conv(i,j) = real(dh_dtau(j,i)*year_sec,sp)
2715  vx_b_g_conv(i,j) = real(vx_b_g(j,i)*year_sec,sp)
2716  vy_b_g_conv(i,j) = real(vy_b_g(j,i)*year_sec,sp)
2717  vz_b_conv(i,j) = real(vz_b(j,i)*year_sec,sp)
2718  vh_b_conv(i,j) = sqrt( vx_b_g_conv(i,j)**2 + vy_b_g_conv(i,j)**2 )
2719  vx_s_g_conv(i,j) = real(vx_s_g(j,i)*year_sec,sp)
2720  vy_s_g_conv(i,j) = real(vy_s_g(j,i)*year_sec,sp)
2721  vz_s_conv(i,j) = real(vz_s(j,i)*year_sec,sp)
2722  vh_s_conv(i,j) = sqrt( vx_s_g_conv(i,j)**2 + vy_s_g_conv(i,j)**2 )
2723  vx_m_g_conv(i,j) = real(vx_m_g(j,i)*year_sec,sp)
2724  vy_m_g_conv(i,j) = real(vy_m_g(j,i)*year_sec,sp)
2725  vh_m_conv(i,j) = sqrt( vx_m_g_conv(i,j)**2 + vy_m_g_conv(i,j)**2 )
2726  temp_b_conv(i,j) = real(temp_b(j,i),sp)
2727  temph_b_conv(i,j) = real(temph_b(j,i),sp)
2728  tau_b_driving_conv(i,j) = real(tau_b_driving(j,i),sp)
2729  tau_b_drag_conv(i,j) = real(tau_b_drag(j,i),sp)
2730  p_b_w_conv(i,j) = real(p_b_w(j,i),sp)
2731  h_w_conv(i,j) = real(H_w(j,i),sp)
2732  q_gl_g_conv(i,j) = real(q_gl_g(j,i)*year_sec,sp)
2733  q_cf_g_conv(i,j) = real(calv_grounded(j,i)*year_sec,sp)
2734  ratio_sl_x_conv(i,j) = real(ratio_sl_x(j,i),sp)
2735  ratio_sl_y_conv(i,j) = real(ratio_sl_y(j,i),sp)
2736 
2737  if (flag_shelfy_stream_x(j,i)) then
2738  flag_shelfy_stream_x_conv(i,j) = 1_i1b
2739  else
2740  flag_shelfy_stream_x_conv(i,j) = 0_i1b
2741  end if
2742 
2743  if (flag_shelfy_stream_y(j,i)) then
2744  flag_shelfy_stream_y_conv(i,j) = 1_i1b
2745  else
2746  flag_shelfy_stream_y_conv(i,j) = 0_i1b
2747  end if
2748 
2749  if (flag_shelfy_stream(j,i)) then
2750  flag_shelfy_stream_conv(i,j) = 1_i1b
2751  else
2752  flag_shelfy_stream_conv(i,j) = 0_i1b
2753  end if
2754 
2755  vis_int_g_conv(i,j) = real(vis_int_g(j,i),sp)
2756 
2757  do kr=0, krmax
2758  temp_r_conv(i,j,kr) = real(temp_r(kr,j,i),sp)
2759  end do
2760 
2761  do kt=0, ktmax
2762  vx_t_conv(i,j,kt) = real(vx_t(kt,j,i)*year_sec,sp)
2763  vy_t_conv(i,j,kt) = real(vy_t(kt,j,i)*year_sec,sp)
2764  vz_t_conv(i,j,kt) = real(vz_t(kt,j,i)*year_sec,sp)
2765  omega_t_conv(i,j,kt) = real(omega_t(kt,j,i),sp)
2766  age_t_conv(i,j,kt) = real(age_t(kt,j,i)*sec_to_year,sp)
2767  enth_t_conv(i,j,kt) = real(enth_t(kt,j,i),sp)
2768  enh_t_conv(i,j,kt) = real(enh_t(kt,j,i),sp)
2769  end do
2770 
2771  do kc=0, kcmax
2772  vx_c_conv(i,j,kc) = real(vx_c(kc,j,i)*year_sec,sp)
2773  vy_c_conv(i,j,kc) = real(vy_c(kc,j,i)*year_sec,sp)
2774  vz_c_conv(i,j,kc) = real(vz_c(kc,j,i)*year_sec,sp)
2775  temp_c_conv(i,j,kc) = real(temp_c(kc,j,i),sp)
2776  age_c_conv(i,j,kc) = real(age_c(kc,j,i)*sec_to_year,sp)
2777  enth_c_conv(i,j,kc) = real(enth_c(kc,j,i),sp)
2778  omega_c_conv(i,j,kc) = real(omega_c(kc,j,i),sp)
2779  enh_c_conv(i,j,kc) = real(enh_c(kc,j,i),sp)
2780  end do
2781 
2782 end do
2783 end do
2784 
2785 !-------- Write data on file --------
2786 
2787 #if (NETCDF==1) /* time-slice file in native binary format */
2788 
2789 write(unit=11) time_conv
2790 if (forcing_flag == 1) then
2791  write(unit=11) delta_ts_conv
2792 else if (forcing_flag == 2) then
2793  write(unit=11) glac_index_conv
2794 else if (forcing_flag == 3) then
2795  glac_index_conv = real(no_value_pos_1,sp) ! dummy value
2796  write(unit=11) glac_index_conv
2797 end if
2798 write(unit=11) z_sl_conv
2799 
2800 write(unit=11) v_tot_conv
2801 write(unit=11) v_af_conv
2802 write(unit=11) a_grounded_conv
2803 write(unit=11) a_floating_conv
2804 
2805 write(unit=11) xi_conv
2806 write(unit=11) eta_conv
2807 write(unit=11) sigma_level_c_conv
2808 write(unit=11) sigma_level_t_conv
2809 write(unit=11) sigma_level_r_conv
2810 write(unit=11) lon_conv
2811 write(unit=11) lat_conv
2812 write(unit=11) lambda_conv
2813 write(unit=11) phi_conv
2814 write(unit=11) temp_s_conv
2815 write(unit=11) accum_conv
2816 write(unit=11) as_perp_conv
2817 write(unit=11) as_perp_apl_conv
2818 write(unit=11) smb_corr_conv
2819 
2820 #if (DISC>0) /* Ice discharge parameterisation */
2821 
2822 write(unit=11) dis_perp_conv
2823 write(unit=11) cst_dist_conv
2824 write(unit=11) cos_grad_tc_conv
2825 write(unit=11) mask_mar_conv
2826 
2827 #endif
2828 
2829 write(unit=11) q_geo_conv
2830 write(unit=11) maske_conv
2831 write(unit=11) maske_old_conv
2832 write(unit=11) n_cts_conv
2833 write(unit=11) kc_cts_conv
2834 write(unit=11) zs_conv
2835 write(unit=11) zm_conv
2836 write(unit=11) zb_conv
2837 write(unit=11) zl_conv
2838 write(unit=11) zl0_conv
2839 write(unit=11) h_cold_conv
2840 write(unit=11) h_temp_conv
2841 write(unit=11) h_conv
2842 write(unit=11) h_r_conv
2843 write(unit=11) q_bm_conv
2844 write(unit=11) q_tld_conv
2845 write(unit=11) am_perp_conv
2846 write(unit=11) qx_conv
2847 write(unit=11) qy_conv
2848 write(unit=11) dzs_dtau_conv
2849 write(unit=11) dzm_dtau_conv
2850 write(unit=11) dzb_dtau_conv
2851 write(unit=11) dzl_dtau_conv
2852 write(unit=11) dh_c_dtau_conv
2853 write(unit=11) dh_t_dtau_conv
2854 write(unit=11) dh_dtau_conv
2855 write(unit=11) vx_b_g_conv
2856 write(unit=11) vy_b_g_conv
2857 write(unit=11) vz_b_conv
2858 write(unit=11) vh_b_conv
2859 write(unit=11) vx_s_g_conv
2860 write(unit=11) vy_s_g_conv
2861 write(unit=11) vz_s_conv
2862 write(unit=11) vh_s_conv
2863 write(unit=11) vx_m_g_conv
2864 write(unit=11) vy_m_g_conv
2865 write(unit=11) vh_m_conv
2866 write(unit=11) temp_b_conv
2867 write(unit=11) temph_b_conv
2868 write(unit=11) tau_b_driving_conv
2869 write(unit=11) tau_b_drag_conv
2870 write(unit=11) p_b_w_conv
2871 write(unit=11) h_w_conv
2872 write(unit=11) q_gl_g_conv
2873 write(unit=11) q_cf_g_conv
2874 
2875 if (flag_3d_output) then
2876 
2877  write(unit=11) vx_c_conv
2878  write(unit=11) vy_c_conv
2879  write(unit=11) vz_c_conv
2880  write(unit=11) vx_t_conv
2881  write(unit=11) vy_t_conv
2882  write(unit=11) vz_t_conv
2883  write(unit=11) temp_c_conv
2884  write(unit=11) omega_t_conv
2885  write(unit=11) temp_r_conv
2886  write(unit=11) enth_c_conv
2887  write(unit=11) enth_t_conv
2888  write(unit=11) omega_c_conv
2889  write(unit=11) enh_c_conv
2890  write(unit=11) enh_t_conv
2891  write(unit=11) age_c_conv
2892  write(unit=11) age_t_conv
2893 
2894 end if
2895 
2896 #elif (NETCDF==2) /* time-slice file in NetCDF format */
2897 
2898 call check( nf90_inq_varid(ncid, 'crs', ncv), thisroutine )
2899 call check( nf90_put_var(ncid, ncv, 0), thisroutine )
2900 
2901 call check( nf90_inq_varid(ncid, 'time', ncv), thisroutine )
2902 call check( nf90_put_var(ncid, ncv, time_conv), thisroutine )
2903 
2904 if (forcing_flag == 1) then
2905 
2906  call check( nf90_inq_varid(ncid, 'delta_ts', ncv), thisroutine )
2907  call check( nf90_put_var(ncid, ncv, delta_ts_conv), thisroutine )
2908 
2909 else if (forcing_flag == 2) then
2910 
2911  call check( nf90_inq_varid(ncid, 'glac_index', ncv), thisroutine )
2912  call check( nf90_put_var(ncid, ncv, glac_index_conv), thisroutine )
2913 
2914 else if (forcing_flag == 3) then
2915 
2916  call check( nf90_inq_varid(ncid, 'glac_index', ncv), thisroutine )
2917  glac_index_conv = real(no_value_pos_1,sp) ! dummy value
2918  call check( nf90_put_var(ncid, ncv, glac_index_conv), thisroutine )
2919 
2920 end if
2921 
2922 call check( nf90_inq_varid(ncid, 'z_sl', ncv), thisroutine )
2923 call check( nf90_put_var(ncid, ncv, z_sl_conv), thisroutine )
2924 
2925 call check( nf90_inq_varid(ncid, 'V_tot', ncv), thisroutine )
2926 call check( nf90_put_var(ncid, ncv, v_tot_conv), thisroutine )
2927 
2928 call check( nf90_inq_varid(ncid, 'V_af', ncv), thisroutine )
2929 call check( nf90_put_var(ncid, ncv, v_af_conv), thisroutine )
2930 
2931 call check( nf90_inq_varid(ncid, 'A_grounded', ncv), thisroutine )
2932 call check( nf90_put_var(ncid, ncv, a_grounded_conv), thisroutine )
2933 
2934 call check( nf90_inq_varid(ncid, 'A_floating', ncv), thisroutine )
2935 call check( nf90_put_var(ncid, ncv, a_floating_conv), thisroutine )
2936 
2937 call check( nf90_inq_varid(ncid, 'x', ncv), thisroutine )
2938 call check( nf90_put_var(ncid, ncv, xi_conv, &
2939  start=nc1cor_i, count=nc1cnt_i), &
2940  thisroutine )
2941 
2942 call check( nf90_inq_varid(ncid, 'y', ncv), thisroutine )
2943 call check( nf90_put_var(ncid, ncv, eta_conv, &
2944  start=nc1cor_j, count=nc1cnt_j), &
2945  thisroutine )
2946 
2947 call check( nf90_inq_varid(ncid, 'sigma_level_c', ncv), thisroutine )
2948 call check( nf90_put_var(ncid, ncv, sigma_level_c_conv, &
2949  start=nc1cor_kc, count=nc1cnt_kc), &
2950  thisroutine )
2951 
2952 call check( nf90_inq_varid(ncid, 'sigma_level_t', ncv), thisroutine )
2953 call check( nf90_put_var(ncid, ncv, sigma_level_t_conv, &
2954  start=nc1cor_kt, count=nc1cnt_kt), &
2955  thisroutine )
2956 
2957 call check( nf90_inq_varid(ncid, 'sigma_level_r', ncv), thisroutine )
2958 call check( nf90_put_var(ncid, ncv, sigma_level_r_conv, &
2959  start=nc1cor_kr, count=nc1cnt_kr), &
2960  thisroutine )
2961 
2962 call check( nf90_inq_varid(ncid, 'lon', ncv), thisroutine )
2963 call check( nf90_put_var(ncid, ncv, lon_conv, &
2964  start=nc2cor_ij, count=nc2cnt_ij), &
2965  thisroutine )
2966 
2967 call check( nf90_inq_varid(ncid, 'lat', ncv), thisroutine )
2968 call check( nf90_put_var(ncid, ncv, lat_conv, &
2969  start=nc2cor_ij, count=nc2cnt_ij), &
2970  thisroutine )
2971 
2972 call check( nf90_inq_varid(ncid, 'lambda', ncv), thisroutine )
2973 call check( nf90_put_var(ncid, ncv, lambda_conv, &
2974  start=nc2cor_ij, count=nc2cnt_ij), &
2975  thisroutine )
2976 
2977 call check( nf90_inq_varid(ncid, 'phi', ncv), thisroutine )
2978 call check( nf90_put_var(ncid, ncv, phi_conv, &
2979  start=nc2cor_ij, count=nc2cnt_ij), &
2980  thisroutine )
2981 
2982 call check( nf90_inq_varid(ncid, 'temp_s', ncv), thisroutine )
2983 call check( nf90_put_var(ncid, ncv, temp_s_conv, &
2984  start=nc2cor_ij, count=nc2cnt_ij), &
2985  thisroutine )
2986 
2987 call check( nf90_inq_varid(ncid, 'prec', ncv), thisroutine )
2988 call check( nf90_put_var(ncid, ncv, accum_conv, &
2989  start=nc2cor_ij, count=nc2cnt_ij), &
2990  thisroutine )
2991 
2992 call check( nf90_inq_varid(ncid, 'as_perp', ncv), thisroutine )
2993 call check( nf90_put_var(ncid, ncv, as_perp_conv, &
2994  start=nc2cor_ij, count=nc2cnt_ij), &
2995  thisroutine )
2996 
2997 call check( nf90_inq_varid(ncid, 'as_perp_apl', ncv), thisroutine )
2998 call check( nf90_put_var(ncid, ncv, as_perp_apl_conv, &
2999  start=nc2cor_ij, count=nc2cnt_ij), &
3000  thisroutine )
3001 
3002 call check( nf90_inq_varid(ncid, 'smb_corr', ncv), thisroutine )
3003 call check( nf90_put_var(ncid, ncv, smb_corr_conv, &
3004  start=nc2cor_ij, count=nc2cnt_ij), &
3005  thisroutine )
3006 
3007 #if (DISC>0) /* Ice discharge parameterisation */
3008 
3009 call check( nf90_inq_varid(ncid, 'dis_perp', ncv), thisroutine )
3010 call check( nf90_put_var(ncid, ncv, dis_perp_conv, &
3011  start=nc2cor_ij, count=nc2cnt_ij), &
3012  thisroutine )
3013 
3014 call check( nf90_inq_varid(ncid, 'cst_dist', ncv), thisroutine )
3015 call check( nf90_put_var(ncid, ncv, cst_dist_conv, &
3016  start=nc2cor_ij, count=nc2cnt_ij), &
3017  thisroutine )
3018 
3019 call check( nf90_inq_varid(ncid, 'cos_grad_tc', ncv), thisroutine )
3020 call check( nf90_put_var(ncid, ncv, cos_grad_tc_conv, &
3021  start=nc2cor_ij, count=nc2cnt_ij), &
3022  thisroutine )
3023 
3024 call check( nf90_inq_varid(ncid, 'mask_mar', ncv), thisroutine )
3025 call check( nf90_put_var(ncid, ncv, mask_mar_conv, &
3026  start=nc2cor_ij, count=nc2cnt_ij), &
3027  thisroutine )
3028 
3029 #endif
3030 
3031 call check( nf90_inq_varid(ncid, 'q_geo', ncv), thisroutine )
3032 call check( nf90_put_var(ncid, ncv, q_geo_conv, &
3033  start=nc2cor_ij, count=nc2cnt_ij), &
3034  thisroutine )
3035 
3036 call check( nf90_inq_varid(ncid, 'maske', ncv), thisroutine )
3037 call check( nf90_put_var(ncid, ncv, maske_conv, &
3038  start=nc2cor_ij, count=nc2cnt_ij), &
3039  thisroutine )
3040 
3041 call check( nf90_inq_varid(ncid, 'maske_old', ncv), thisroutine )
3042 call check( nf90_put_var(ncid, ncv, maske_old_conv, &
3043  start=nc2cor_ij, count=nc2cnt_ij), &
3044  thisroutine )
3045 
3046 call check( nf90_inq_varid(ncid, 'mask_run', ncv), thisroutine )
3047 call check( nf90_put_var(ncid, ncv, mask_run_conv, &
3048  start=nc2cor_ij, count=nc2cnt_ij), &
3049  thisroutine )
3050 
3051 call check( nf90_inq_varid(ncid, 'n_cts', ncv), thisroutine )
3052 call check( nf90_put_var(ncid, ncv, n_cts_conv, &
3053  start=nc2cor_ij, count=nc2cnt_ij), &
3054  thisroutine )
3055 
3056 call check( nf90_inq_varid(ncid, 'kc_cts', ncv), thisroutine )
3057 call check( nf90_put_var(ncid, ncv, kc_cts_conv, &
3058  start=nc2cor_ij, count=nc2cnt_ij), &
3059  thisroutine )
3060 
3061 call check( nf90_inq_varid(ncid, 'zs', ncv), thisroutine )
3062 call check( nf90_put_var(ncid, ncv, zs_conv, &
3063  start=nc2cor_ij, count=nc2cnt_ij), &
3064  thisroutine )
3065 
3066 call check( nf90_inq_varid(ncid, 'zm', ncv), thisroutine )
3067 call check( nf90_put_var(ncid, ncv, zm_conv, &
3068  start=nc2cor_ij, count=nc2cnt_ij), &
3069  thisroutine )
3070 
3071 call check( nf90_inq_varid(ncid, 'zb', ncv), thisroutine )
3072 call check( nf90_put_var(ncid, ncv, zb_conv, &
3073  start=nc2cor_ij, count=nc2cnt_ij), &
3074  thisroutine )
3075 
3076 call check( nf90_inq_varid(ncid, 'zl', ncv), thisroutine )
3077 call check( nf90_put_var(ncid, ncv, zl_conv, &
3078  start=nc2cor_ij, count=nc2cnt_ij), &
3079  thisroutine )
3080 
3081 call check( nf90_inq_varid(ncid, 'zl0', ncv), thisroutine )
3082 call check( nf90_put_var(ncid, ncv, zl0_conv, &
3083  start=nc2cor_ij, count=nc2cnt_ij), &
3084  thisroutine )
3085 
3086 call check( nf90_inq_varid(ncid, 'H_cold', ncv), thisroutine )
3087 call check( nf90_put_var(ncid, ncv, h_cold_conv, &
3088  start=nc2cor_ij, count=nc2cnt_ij), &
3089  thisroutine )
3090 
3091 call check( nf90_inq_varid(ncid, 'H_temp', ncv), thisroutine )
3092 call check( nf90_put_var(ncid, ncv, h_temp_conv, &
3093  start=nc2cor_ij, count=nc2cnt_ij), &
3094  thisroutine )
3095 
3096 call check( nf90_inq_varid(ncid, 'H', ncv), thisroutine )
3097 call check( nf90_put_var(ncid, ncv, h_conv, &
3098  start=nc2cor_ij, count=nc2cnt_ij), &
3099  thisroutine )
3100 
3101 call check( nf90_inq_varid(ncid, 'H_R', ncv), thisroutine )
3102 call check( nf90_put_var(ncid, ncv, h_r_conv), thisroutine )
3103 
3104 call check( nf90_inq_varid(ncid, 'Q_bm', ncv), thisroutine )
3105 call check( nf90_put_var(ncid, ncv, q_bm_conv, &
3106  start=nc2cor_ij, count=nc2cnt_ij), &
3107  thisroutine )
3108 
3109 call check( nf90_inq_varid(ncid, 'Q_tld', ncv), thisroutine )
3110 call check( nf90_put_var(ncid, ncv, q_tld_conv, &
3111  start=nc2cor_ij, count=nc2cnt_ij), &
3112  thisroutine )
3113 
3114 call check( nf90_inq_varid(ncid, 'am_perp', ncv), thisroutine )
3115 call check( nf90_put_var(ncid, ncv, am_perp_conv, &
3116  start=nc2cor_ij, count=nc2cnt_ij), &
3117  thisroutine )
3118 
3119 call check( nf90_inq_varid(ncid, 'qx', ncv), thisroutine )
3120 call check( nf90_put_var(ncid, ncv, qx_conv, &
3121  start=nc2cor_ij, count=nc2cnt_ij), &
3122  thisroutine )
3123 
3124 call check( nf90_inq_varid(ncid, 'qy', ncv), thisroutine )
3125 call check( nf90_put_var(ncid, ncv, qy_conv, &
3126  start=nc2cor_ij, count=nc2cnt_ij), &
3127  thisroutine )
3128 
3129 call check( nf90_inq_varid(ncid, 'dzs_dt', ncv), thisroutine )
3130 call check( nf90_put_var(ncid, ncv, dzs_dtau_conv, &
3131  start=nc2cor_ij, count=nc2cnt_ij), &
3132  thisroutine )
3133 
3134 call check( nf90_inq_varid(ncid, 'dzm_dt', ncv), thisroutine )
3135 call check( nf90_put_var(ncid, ncv, dzm_dtau_conv, &
3136  start=nc2cor_ij, count=nc2cnt_ij), &
3137  thisroutine )
3138 
3139 call check( nf90_inq_varid(ncid, 'dzb_dt', ncv), thisroutine )
3140 call check( nf90_put_var(ncid, ncv, dzb_dtau_conv, &
3141  start=nc2cor_ij, count=nc2cnt_ij), &
3142  thisroutine )
3143 
3144 call check( nf90_inq_varid(ncid, 'dzl_dt', ncv), thisroutine )
3145 call check( nf90_put_var(ncid, ncv, dzl_dtau_conv, &
3146  start=nc2cor_ij, count=nc2cnt_ij), &
3147  thisroutine )
3148 
3149 call check( nf90_inq_varid(ncid, 'dH_c_dt', ncv), thisroutine )
3150 call check( nf90_put_var(ncid, ncv, dh_c_dtau_conv, &
3151  start=nc2cor_ij, count=nc2cnt_ij), &
3152  thisroutine )
3153 
3154 call check( nf90_inq_varid(ncid, 'dH_t_dt', ncv), thisroutine )
3155 call check( nf90_put_var(ncid, ncv, dh_t_dtau_conv, &
3156  start=nc2cor_ij, count=nc2cnt_ij), &
3157  thisroutine )
3158 
3159 call check( nf90_inq_varid(ncid, 'dH_dt', ncv), thisroutine )
3160 call check( nf90_put_var(ncid, ncv, dh_dtau_conv, &
3161  start=nc2cor_ij, count=nc2cnt_ij), &
3162  thisroutine )
3163 
3164 call check( nf90_inq_varid(ncid, 'vx_b_g', ncv), thisroutine )
3165 call check( nf90_put_var(ncid, ncv, vx_b_g_conv, &
3166  start=nc2cor_ij, count=nc2cnt_ij), &
3167  thisroutine )
3168 
3169 call check( nf90_inq_varid(ncid, 'vy_b_g', ncv), thisroutine )
3170 call check( nf90_put_var(ncid, ncv, vy_b_g_conv, &
3171  start=nc2cor_ij, count=nc2cnt_ij), &
3172  thisroutine )
3173 
3174 call check( nf90_inq_varid(ncid, 'vz_b', ncv), thisroutine )
3175 call check( nf90_put_var(ncid, ncv, vz_b_conv, &
3176  start=nc2cor_ij, count=nc2cnt_ij), &
3177  thisroutine )
3178 
3179 call check( nf90_inq_varid(ncid, 'vh_b', ncv), thisroutine )
3180 call check( nf90_put_var(ncid, ncv, vh_b_conv, &
3181  start=nc2cor_ij, count=nc2cnt_ij), &
3182  thisroutine )
3183 
3184 call check( nf90_inq_varid(ncid, 'vx_s_g', ncv), thisroutine )
3185 call check( nf90_put_var(ncid, ncv, vx_s_g_conv, &
3186  start=nc2cor_ij, count=nc2cnt_ij), &
3187  thisroutine )
3188 
3189 call check( nf90_inq_varid(ncid, 'vy_s_g', ncv), thisroutine )
3190 call check( nf90_put_var(ncid, ncv, vy_s_g_conv, &
3191  start=nc2cor_ij, count=nc2cnt_ij), &
3192  thisroutine )
3193 
3194 call check( nf90_inq_varid(ncid, 'vz_s', ncv), thisroutine )
3195 call check( nf90_put_var(ncid, ncv, vz_s_conv, &
3196  start=nc2cor_ij, count=nc2cnt_ij), &
3197  thisroutine )
3198 
3199 call check( nf90_inq_varid(ncid, 'vh_s', ncv), thisroutine )
3200 call check( nf90_put_var(ncid, ncv, vh_s_conv, &
3201  start=nc2cor_ij, count=nc2cnt_ij), &
3202  thisroutine )
3203 
3204 call check( nf90_inq_varid(ncid, 'vx_m_g', ncv), thisroutine )
3205 call check( nf90_put_var(ncid, ncv, vx_m_g_conv, &
3206  start=nc2cor_ij, count=nc2cnt_ij), &
3207  thisroutine )
3208 
3209 call check( nf90_inq_varid(ncid, 'vy_m_g', ncv), thisroutine )
3210 call check( nf90_put_var(ncid, ncv, vy_m_g_conv, &
3211  start=nc2cor_ij, count=nc2cnt_ij), &
3212  thisroutine )
3213 
3214 call check( nf90_inq_varid(ncid, 'vh_m', ncv), thisroutine )
3215 call check( nf90_put_var(ncid, ncv, vh_m_conv, &
3216  start=nc2cor_ij, count=nc2cnt_ij), &
3217  thisroutine )
3218 
3219 call check( nf90_inq_varid(ncid, 'temp_b', ncv), thisroutine )
3220 call check( nf90_put_var(ncid, ncv, temp_b_conv, &
3221  start=nc2cor_ij, count=nc2cnt_ij), &
3222  thisroutine )
3223 
3224 call check( nf90_inq_varid(ncid, 'temph_b', ncv), thisroutine )
3225 call check( nf90_put_var(ncid, ncv, temph_b_conv, &
3226  start=nc2cor_ij, count=nc2cnt_ij), &
3227  thisroutine )
3228 
3229 call check( nf90_inq_varid(ncid, 'tau_b_driving', ncv), thisroutine )
3230 call check( nf90_put_var(ncid, ncv, tau_b_driving_conv, &
3231  start=nc2cor_ij, count=nc2cnt_ij), &
3232  thisroutine )
3233 
3234 call check( nf90_inq_varid(ncid, 'tau_b_drag', ncv), thisroutine )
3235 call check( nf90_put_var(ncid, ncv, tau_b_drag_conv, &
3236  start=nc2cor_ij, count=nc2cnt_ij), &
3237  thisroutine )
3238 
3239 call check( nf90_inq_varid(ncid, 'p_b_w', ncv), thisroutine )
3240 call check( nf90_put_var(ncid, ncv, p_b_w_conv, &
3241  start=nc2cor_ij, count=nc2cnt_ij), &
3242  thisroutine )
3243 
3244 call check( nf90_inq_varid(ncid, 'H_w', ncv), thisroutine )
3245 call check( nf90_put_var(ncid, ncv, h_w_conv, &
3246  start=nc2cor_ij, count=nc2cnt_ij), &
3247  thisroutine )
3248 
3249 call check( nf90_inq_varid(ncid, 'q_gl_g', ncv), thisroutine )
3250 call check( nf90_put_var(ncid, ncv, q_gl_g_conv, &
3251  start=nc2cor_ij, count=nc2cnt_ij), &
3252  thisroutine )
3253 
3254 call check( nf90_inq_varid(ncid, 'q_cf_g', ncv), thisroutine )
3255 call check( nf90_put_var(ncid, ncv, q_cf_g_conv, &
3256  start=nc2cor_ij, count=nc2cnt_ij), &
3257  thisroutine )
3258 
3259 call check( nf90_inq_varid(ncid, 'ratio_sl_x', ncv), thisroutine )
3260 call check( nf90_put_var(ncid, ncv, ratio_sl_x_conv, &
3261  start=nc2cor_ij, count=nc2cnt_ij), &
3262  thisroutine )
3263 
3264 call check( nf90_inq_varid(ncid, 'ratio_sl_y', ncv), thisroutine )
3265 call check( nf90_put_var(ncid, ncv, ratio_sl_y_conv, &
3266  start=nc2cor_ij, count=nc2cnt_ij), &
3267  thisroutine )
3268 
3269 call check( nf90_inq_varid(ncid, 'flag_shelfy_stream_x', ncv), thisroutine )
3270 call check( nf90_put_var(ncid, ncv, flag_shelfy_stream_x_conv, &
3271  start=nc2cor_ij, count=nc2cnt_ij), &
3272  thisroutine )
3273 
3274 call check( nf90_inq_varid(ncid, 'flag_shelfy_stream_y', ncv), thisroutine )
3275 call check( nf90_put_var(ncid, ncv, flag_shelfy_stream_y_conv, &
3276  start=nc2cor_ij, count=nc2cnt_ij), &
3277  thisroutine )
3278 
3279 call check( nf90_inq_varid(ncid, 'flag_shelfy_stream', ncv), thisroutine )
3280 call check( nf90_put_var(ncid, ncv, flag_shelfy_stream_conv, &
3281  start=nc2cor_ij, count=nc2cnt_ij), &
3282  thisroutine )
3283 
3284 call check( nf90_inq_varid(ncid, 'vis_int_g', ncv), thisroutine )
3285 call check( nf90_put_var(ncid, ncv, vis_int_g_conv, &
3286  start=nc2cor_ij, count=nc2cnt_ij), &
3287  thisroutine )
3288 
3289 if (flag_3d_output) then
3290 
3291  call check( nf90_inq_varid(ncid, 'vx_c', ncv), thisroutine )
3292  call check( nf90_put_var(ncid, ncv, vx_c_conv, &
3293  start=nc3cor_ijkc, count=nc3cnt_ijkc), &
3294  thisroutine )
3295 
3296  call check( nf90_inq_varid(ncid, 'vy_c', ncv), thisroutine )
3297  call check( nf90_put_var(ncid, ncv, vy_c_conv, &
3298  start=nc3cor_ijkc, count=nc3cnt_ijkc), &
3299  thisroutine )
3300 
3301  call check( nf90_inq_varid(ncid, 'vz_c', ncv), thisroutine )
3302  call check( nf90_put_var(ncid, ncv, vz_c_conv, &
3303  start=nc3cor_ijkc, count=nc3cnt_ijkc), &
3304  thisroutine )
3305 
3306  call check( nf90_inq_varid(ncid, 'vx_t', ncv), thisroutine )
3307  call check( nf90_put_var(ncid, ncv, vx_t_conv, &
3308  start=nc3cor_ijkt, count=nc3cnt_ijkt), &
3309  thisroutine )
3310 
3311  call check( nf90_inq_varid(ncid, 'vy_t', ncv), thisroutine )
3312  call check( nf90_put_var(ncid, ncv, vy_t_conv, &
3313  start=nc3cor_ijkt, count=nc3cnt_ijkt), &
3314  thisroutine )
3315 
3316  call check( nf90_inq_varid(ncid, 'vz_t', ncv), thisroutine )
3317  call check( nf90_put_var(ncid, ncv, vz_t_conv, &
3318  start=nc3cor_ijkt, count=nc3cnt_ijkt), &
3319  thisroutine )
3320 
3321  call check( nf90_inq_varid(ncid, 'temp_c', ncv), thisroutine )
3322  call check( nf90_put_var(ncid, ncv, temp_c_conv, &
3323  start=nc3cor_ijkc, count=nc3cnt_ijkc), &
3324  thisroutine )
3325 
3326  call check( nf90_inq_varid(ncid, 'omega_t', ncv), thisroutine )
3327  call check( nf90_put_var(ncid, ncv, omega_t_conv, &
3328  start=nc3cor_ijkt, count=nc3cnt_ijkt), &
3329  thisroutine )
3330 
3331  call check( nf90_inq_varid(ncid, 'temp_r', ncv), thisroutine )
3332  call check( nf90_put_var(ncid, ncv, temp_r_conv, &
3333  start=nc3cor_ijkr, count=nc3cnt_ijkr), &
3334  thisroutine )
3335 
3336  call check( nf90_inq_varid(ncid, 'enth_c', ncv), thisroutine )
3337  call check( nf90_put_var(ncid, ncv, enth_c_conv, &
3338  start=nc3cor_ijkc, count=nc3cnt_ijkc), &
3339  thisroutine )
3340 
3341  call check( nf90_inq_varid(ncid, 'enth_t', ncv), thisroutine )
3342  call check( nf90_put_var(ncid, ncv, enth_t_conv, &
3343  start=nc3cor_ijkt, count=nc3cnt_ijkt), &
3344  thisroutine )
3345 
3346  call check( nf90_inq_varid(ncid, 'omega_c', ncv), thisroutine )
3347  call check( nf90_put_var(ncid, ncv, omega_c_conv, &
3348  start=nc3cor_ijkc, count=nc3cnt_ijkc), &
3349  thisroutine )
3350 
3351  call check( nf90_inq_varid(ncid, 'enh_c', ncv), thisroutine )
3352  call check( nf90_put_var(ncid, ncv, enh_c_conv, &
3353  start=nc3cor_ijkc, count=nc3cnt_ijkc), &
3354  thisroutine )
3355 
3356  call check( nf90_inq_varid(ncid, 'enh_t', ncv), thisroutine )
3357  call check( nf90_put_var(ncid, ncv, enh_t_conv, &
3358  start=nc3cor_ijkt, count=nc3cnt_ijkt), &
3359  thisroutine )
3360 
3361  call check( nf90_inq_varid(ncid, 'age_c', ncv), thisroutine )
3362  call check( nf90_put_var(ncid, ncv, age_c_conv, &
3363  start=nc3cor_ijkc, count=nc3cnt_ijkc), &
3364  thisroutine )
3365 
3366  call check( nf90_inq_varid(ncid, 'age_t', ncv), thisroutine )
3367  call check( nf90_put_var(ncid, ncv, age_t_conv, &
3368  start=nc3cor_ijkt, count=nc3cnt_ijkt), &
3369  thisroutine )
3370 
3371 end if
3372 
3373 #endif
3374 
3375 !-------- Close file --------
3376 
3377 #if (NETCDF==1) /* time-slice file in native binary format */
3378 
3379 close(unit=11, status='keep')
3380 
3381 #elif (NETCDF==2) /* time-slice file in NetCDF format */
3382 
3383 call check( nf90_sync(ncid), thisroutine )
3384 call check( nf90_close(ncid), thisroutine )
3385 
3386 deallocate(coord_id)
3387 
3388 #endif
3389 
3390 !-------- Increase file counter --------
3391 
3392 ndat = ndat+1
3393 
3394 if (flag_3d_output) then
3395  ndat3d = ndat
3396 else
3397  ndat2d = ndat
3398 end if
3399 
3400 end subroutine output1
3401 
3402 !-------------------------------------------------------------------------------
3403 !> Writing of time-series data on file in ASCII format
3404 !! (and optionally in NetCDF format).
3405 !<------------------------------------------------------------------------------
3406 subroutine output2(time, dxi, deta, delta_ts, glac_index, z_sl)
3408 #if (NETCDF>1)
3409  use netcdf
3410  use nc_check_m
3411 #endif
3412 
3413 #if (DISC>0)
3414  use discharge_workers_m, only: dt_glann, dt_sub
3415 #endif
3416 
3417 implicit none
3418 
3419 real(dp), intent(in) :: time, dxi, deta, delta_ts, glac_index, z_sl
3420 
3421 integer(i4b) :: i, j
3422 integer(i4b) :: n_base, n_tempbase
3423 real(dp) :: time_val, &
3424  V_tot, V_grounded, V_floating, A_tot, A_grounded, A_floating, &
3425  V_af, V_sle, V_temp, A_temp, &
3426  Q_s, Q_b, Q_temp, &
3427  H_max, H_t_max, zs_max, vs_max, Tbh_max, &
3428  accum_tot, calv_tot, mbp, &
3429  MB, LMT, OMT, LMB, OMB, LMH, OMH, ODT, ODH, &
3430  disc_lsc, disc_ssc
3431 real(dp) :: x_pos, y_pos
3432 real(dp), dimension(0:JMAX,0:IMAX) :: H, H_cold, H_temp
3433 real(dp) :: V_gr_redu, A_surf, rhosw_rho_ratio
3434 real(dp) :: vs_help, Tbh_help
3435 real(dp) :: H_ave_sed, Tbh_ave_sed, Atb_sed
3436 real(dp) :: sum_area_sed
3437 
3438 #if (NETCDF>1)
3439 integer(i4b), save :: ncid
3440 integer(i4b) :: ncd, ncv, nc1d
3441 integer(i4b) :: nc1cor(1), nc1cnt(1)
3442 integer(i4b) :: n_sync
3443 real(dp) :: dV_dt, precip_tot, runoff_tot, bmb_tot, mb_resid
3444 real(dp), save :: time_add_offset_val
3445 character(len= 16) :: ch_date, ch_time, ch_zone
3446 character(len=256) :: filename, filename_with_path, buffer
3447 logical, save :: grads_nc_tweaks
3448 #endif
3449 
3450 integer(i4b), save :: counter = 0
3451 logical, save :: firstcall = .true.
3452 
3453 character(len=64), parameter :: thisroutine = 'output2'
3454 
3455 character(len=128), parameter :: &
3456  fmt1 = '(1pe13.6,2(1pe13.4),/,' // &
3457  '13x,6(1pe13.4),/,' // &
3458  '26x,3(1pe13.4),/,' // &
3459  '26x,5(1pe13.4),/)'
3460 
3461 character(len=128), parameter :: &
3462  fmt2 = '(1pe13.6,2(1pe13.4),/,' // &
3463  '13x,3(1pe13.4),/)'
3464 
3465 character(len=128), parameter :: &
3466  fmt3 = '(1pe13.6,3(1pe13.4),/)'
3467 
3468 counter = counter + 1
3469 
3470 !-------- Ice thickness --------
3471 
3472 h = h_c + h_t
3473 
3474 !-------- Thickness of the cold and temperate layers --------
3475 
3476 h_cold = 0.0_dp
3477 h_temp = 0.0_dp
3478 
3479 #if (CALCMOD==1)
3480 do i=1, imax-1
3481 do j=1, jmax-1
3482  h_temp(j,i) = h_t(j,i)
3483 end do
3484 end do
3485 #elif (CALCMOD==0 || CALCMOD==2 || CALCMOD==3 || CALCMOD==-1)
3486 do i=1, imax-1
3487 do j=1, jmax-1
3488  h_temp(j,i) = h_c(j,i)*eaz_c_quotient(kc_cts(j,i))
3489 end do
3490 end do
3491 #endif
3492 
3493 h_cold = h - h_temp
3494 
3495 !-------- Maximum ice elevation and thickness, ice volume,
3496 ! sea-level equivalent, volume of the temperate ice,
3497 ! area, area covered by temperate ice,
3498 ! freshwater production due to melting and calving,
3499 ! water drainage due to basal melting,
3500 ! water drainage from the temperate layer,
3501 ! maximum thickness of the temperate layer,
3502 ! maximum surface velocity,
3503 ! maximum basal temperature rel. to pmp --------
3504 
3505 #if (defined(ANT) \
3506  || defined(grl) \
3507  || defined(scand) \
3508  || defined(nhem) \
3509  || defined(tibet) \
3510  || defined(asf) \
3511  || defined(emtp2sge) \
3512  || defined(xyz)) /* terrestrial ice sheet */
3513 
3514 rhosw_rho_ratio = rho_sw/rho
3515 
3516 #elif (defined(NMARS) || defined(SMARS)) /* Martian ice sheet */
3517 
3518 rhosw_rho_ratio = 0.0_dp ! dummy value
3519 
3520 #else
3521 
3522 stop ' >>> output2: No valid domain specified!'
3523 
3524 #endif
3525 
3526 v_grounded = 0.0_dp
3527 v_gr_redu = 0.0_dp
3528 v_floating = 0.0_dp
3529 a_grounded = 0.0_dp
3530 a_floating = 0.0_dp
3531 v_temp = 0.0_dp
3532 a_temp = 0.0_dp
3533 q_s = 0.0_dp
3534 q_b = 0.0_dp
3535 q_temp = 0.0_dp
3536 h_max = 0.0_dp
3537 h_t_max = 0.0_dp
3538 zs_max = no_value_neg_2
3539 vs_max = 0.0_dp
3540 tbh_max = no_value_neg_2
3541 accum_tot = 0.0_dp
3542 
3543 do i=1, imax-1
3544 do j=1, jmax-1
3545 
3546  if (maske(j,i)==0_i2b) then ! grounded ice
3547 
3548  if (zs(j,i) > zs_max) zs_max = zs(j,i)
3549  if (h(j,i) > h_max ) h_max = h(j,i)
3550  if (h_temp(j,i) > h_t_max) h_t_max = h_temp(j,i)
3551 
3552  v_grounded = v_grounded + h(j,i) *area(j,i)
3553  v_temp = v_temp + h_temp(j,i)*area(j,i)
3554 
3555 #if (defined(ANT) \
3556  || defined(grl) \
3557  || defined(scand) \
3558  || defined(nhem) \
3559  || defined(tibet) \
3560  || defined(asf) \
3561  || defined(emtp2sge) \
3562  || defined(xyz)) /* terrestrial ice sheet */
3563 
3564  v_gr_redu = v_gr_redu &
3565  + rhosw_rho_ratio*max((z_sl-zl(j,i)),0.0_dp)*area(j,i)
3566 
3567 #endif
3568 
3569  a_grounded = a_grounded + area(j,i)
3570 
3571  if (n_cts(j,i) /= -1) a_temp = a_temp + area(j,i)
3572 
3573  vs_help = sqrt( &
3574  0.25_dp*(vx_c(kcmax,j,i)+vx_c(kcmax,j,i-1))**2 &
3575  +0.25_dp*(vy_c(kcmax,j,i)+vy_c(kcmax,j-1,i))**2 )
3576  if (vs_help > vs_max) vs_max = vs_help
3577 
3578  if (n_cts(j,i) >= 0) then ! temperate base
3579  tbh_max = 0.0_dp
3580  else ! cold base
3581  tbh_help = min((temp_c(0,j,i)-temp_c_m(0,j,i)), 0.0_dp)
3582  if (tbh_help > tbh_max) tbh_max = tbh_help
3583  end if
3584 
3585  accum_tot = accum_tot + accum(j,i)*area(j,i)
3586 
3587  else if (maske(j,i)==3_i2b) then ! floating ice
3588  ! (basal temperature assumed to be below
3589  ! the pressure melting point for pure ice)
3590 
3591  if (zs(j,i) > zs_max) zs_max = zs(j,i)
3592  if (h(j,i) > h_max) h_max = h(j,i)
3593 
3594  v_floating = v_floating + h(j,i)*area(j,i)
3595  a_floating = a_floating + area(j,i)
3596 
3597  vs_help = sqrt( &
3598  0.25_dp*(vx_c(kcmax,j,i)+vx_c(kcmax,j,i-1))**2 &
3599  +0.25_dp*(vy_c(kcmax,j,i)+vy_c(kcmax,j-1,i))**2 )
3600  if (vs_help > vs_max) vs_max = vs_help
3601 
3602  tbh_help = min((temp_c(0,j,i)-temp_c_m(0,j,i)), 0.0_dp)
3603  if (tbh_help > tbh_max) tbh_max = tbh_help
3604 
3605  accum_tot = accum_tot + accum(j,i)*area(j,i)
3606 
3607  end if
3608 
3609  if ( (maske(j,i)==0_i2b).or.(maske_old(j,i)==0_i2b) &
3610  .or.(maske(j,i)==3_i2b).or.(maske_old(j,i)==3_i2b) ) then
3611  ! grounded or floating ice before or after the time step
3612 
3613  q_s = q_s + as_perp_apl(j,i) * area(j,i)
3614  q_b = q_b + q_bm(j,i) * area(j,i)
3615 
3616  end if
3617 
3618  if ( (maske(j,i)==0_i2b).or.(maske_old(j,i)==0_i2b) ) then
3619  ! grounded ice before or after the time step
3620 
3621  q_temp = q_temp + q_tld(j,i) * area(j,i)
3622 
3623  end if
3624 
3625 end do
3626 end do
3627 
3628 !-------- Conversion --------
3629 
3630 #if (defined(ANT) \
3631  || defined(grl) \
3632  || defined(scand) \
3633  || defined(nhem) \
3634  || defined(tibet) \
3635  || defined(asf) \
3636  || defined(emtp2sge) \
3637  || defined(xyz)) /* terrestrial ice sheet */
3638 
3639 a_surf = 3.61132e+14_dp ! global ocean area, in m^2
3640 
3641 v_af = v_grounded - v_gr_redu
3642 v_sle = v_af*(rho/rho_w)/a_surf ! m^3 ice equiv./m^2 -> m water equiv.
3643 
3644 #elif (defined(NMARS) || defined(SMARS)) /* Martian ice sheet */
3645 
3646 a_surf = 1.4441e+14_dp ! surface area of Mars, in m^2
3647 v_sle = v_grounded*(rho_i/rho_w)*(1.0_dp-frac_dust)/a_surf
3648  ! m^3 (ice+dust) equiv./m^2 -> m water equiv.
3649 
3650 #endif
3651 
3652 #if (!defined(OUT_TIMES) || OUT_TIMES==1)
3653 time_val = time *sec_to_year ! s -> a
3654 #elif (OUT_TIMES==2)
3655 time_val = (time+year_zero) *sec_to_year ! s -> a
3656 #else
3657 stop ' >>> output2: OUT_TIMES must be either 1 or 2!'
3658 #endif
3659 
3660 v_tot = v_grounded + v_floating ! in m^3
3661 a_tot = a_grounded + a_floating ! in m^2
3662 
3663 vs_max = vs_max *year_sec ! m/s -> m/a
3664 
3665 #if (defined(ANT) \
3666  || defined(grl) \
3667  || defined(scand) \
3668  || defined(nhem) \
3669  || defined(tibet) \
3670  || defined(asf) \
3671  || defined(emtp2sge) \
3672  || defined(xyz)) /* terrestrial ice sheet */
3673 
3674 q_s = q_s *year_sec *(rho/rho_w)
3675  ! m^3/s ice equiv. -> m^3/a water equiv.
3676 q_b = q_b *year_sec *(rho/rho_w)
3677  ! m^3/s ice equiv. -> m^3/a water equiv.
3678 q_temp = q_temp *year_sec *(rho/rho_w)
3679  ! m^3/s ice equiv. -> m^3/a water equiv.
3680 
3681 accum_tot = accum_tot *year_sec *(rho/rho_w)
3682  ! m^3/s ice equiv. -> m^3/a water equiv.
3683 
3684 #elif (defined(NMARS) || defined(SMARS)) /* Martian ice sheet */
3685 
3686 q_s = q_s *year_sec *(rho_i/rho_w)*(1.0_dp-frac_dust)
3687  ! m^3/s (ice+dust) equiv. -> m^3/a water equiv.
3688 q_b = q_b *year_sec *(rho_i/rho_w)*(1.0_dp-frac_dust)
3689  ! m^3/s (ice+dust) equiv. -> m^3/a water equiv.
3690 q_temp = q_temp *year_sec *(rho_i/rho_w)*(1.0_dp-frac_dust)
3691  ! m^3/s (ice+dust) equiv. -> m^3/a water equiv.
3692 #endif
3693 
3694 !-------- Writing of data on time-series file --------
3695 
3696 if (forcing_flag == 1) then
3697 
3698  write(unit=12, fmt=trim(fmt1)) time_val, delta_ts, z_sl, &
3699  v_tot, v_grounded, v_floating, a_tot, a_grounded, a_floating, &
3700  v_sle, v_temp, a_temp, &
3701  h_max, h_t_max, zs_max, vs_max, tbh_max
3702 
3703 else if (forcing_flag == 2) then
3704 
3705  write(unit=12, fmt=trim(fmt1)) time_val, glac_index, z_sl, &
3706  v_tot, v_grounded, v_floating, a_tot, a_grounded, a_floating, &
3707  v_sle, v_temp, a_temp, &
3708  h_max, h_t_max, zs_max, vs_max, tbh_max
3709 
3710 else if (forcing_flag == 3) then
3711 
3712  write(unit=12, fmt=trim(fmt1)) time_val, no_value_pos_1, z_sl, &
3713  v_tot, v_grounded, v_floating, a_tot, a_grounded, a_floating, &
3714  v_sle, v_temp, a_temp, &
3715  h_max, h_t_max, zs_max, vs_max, tbh_max
3716 
3717 end if
3718 
3719 !-------- Special output --------
3720 
3721 ! ------ Time-series data for the sediment region of HEINO
3722 
3723 #if (defined(XYZ))
3724 
3725 #if (defined(HEINO))
3726 
3727 ! ---- Average ice thickness, average basal temperature rel. to pmp,
3728 ! temperate basal area
3729 
3730 h_ave_sed = 0.0_dp
3731 tbh_ave_sed = 0.0_dp
3732 atb_sed = 0.0_dp
3733 sum_area_sed = 0.0_dp
3734 
3735 do i=1, imax-1
3736 do j=1, jmax-1
3737 
3738  if (maske_sedi(j,i)==7_i2b) then ! sediment region
3739 
3740  sum_area_sed = sum_area_sed + area(j,i)
3741 
3742  h_ave_sed = h_ave_sed + area(j,i)*h(j,i)
3743 
3744  if (n_cts(j,i) /= -1) then ! temperate base
3745  tbh_help = 0.0_dp
3746  else ! cold base
3747  tbh_help = min((temp_c(0,j,i)-temp_c_m(0,j,i)), 0.0_dp)
3748  end if
3749  tbh_ave_sed = tbh_ave_sed + area(j,i)*tbh_help
3750 
3751  if (n_cts(j,i) /= -1) atb_sed = atb_sed + area(j,i)
3752 
3753  end if
3754 
3755 end do
3756 end do
3757 
3758 if (sum_area_sed > eps) then
3759  h_ave_sed = h_ave_sed / sum_area_sed
3760  tbh_ave_sed = tbh_ave_sed / sum_area_sed
3761 else
3762  stop ' >>> output2: No sediment area found!'
3763 end if
3764 
3765 ! ---- Writing of data on time-series file
3766 
3767 if (forcing_flag == 1) then
3768  write(unit=15, fmt=trim(fmt2)) time_val, delta_ts, z_sl, &
3769  h_ave_sed, tbh_ave_sed, atb_sed
3770 else if (forcing_flag == 2) then
3771  write(unit=15, fmt=trim(fmt2)) time_val, glac_index, z_sl, &
3772  h_ave_sed, tbh_ave_sed, atb_sed
3773 else if (forcing_flag == 3) then
3774  write(unit=15, fmt=trim(fmt2)) time_val, no_value_pos_1, z_sl, &
3775  h_ave_sed, tbh_ave_sed, atb_sed
3776 end if
3777 
3778 #endif
3779 
3780 #endif
3781 
3782 !-------- Extended time-series file in NetCDF format --------
3783 
3784 #if (NETCDF>1)
3785 
3786 if (firstcall) then
3787 
3788 ! ------ Open NetCDF file
3789 
3790  filename = trim(runname)//'_ser.nc'
3791  filename_with_path = trim(outpath)//'/'//trim(filename)
3792 
3793  call check( nf90_create(trim(filename_with_path), nf90_noclobber, ncid), &
3794  thisroutine )
3795 
3796  ncid_ser = ncid
3797 
3798 ! ------ Global attributes
3799 
3800  buffer = 'Time-series output of simulation '//trim(runname)
3801  call check( nf90_put_att(ncid, nf90_global, 'title', trim(buffer)), &
3802  thisroutine )
3803 
3804  buffer = 'Institute of Low Temperature Science, Hokkaido University, '// &
3805  'Sapporo, Japan'
3806  call check( nf90_put_att(ncid, nf90_global, 'institution', trim(buffer)), &
3807  thisroutine )
3808 
3809  buffer = 'SICOPOLIS Version '//version
3810  call check( nf90_put_att(ncid, nf90_global, 'source', trim(buffer)), &
3811  thisroutine )
3812 
3813  call date_and_time(ch_date, ch_time, ch_zone)
3814  buffer = ch_date(1:4)//'-'//ch_date(5:6)//'-'//ch_date(7:8)//' '// &
3815  ch_time(1:2)//':'//ch_time(3:4)//':'//ch_time(5:6)//' '// &
3816  ch_zone(1:3)//':'//ch_zone(4:5)//' - Data produced'
3817  call check( nf90_put_att(ncid, nf90_global, 'history', trim(buffer)), &
3818  thisroutine )
3819 
3820  buffer = 'http://www.sicopolis.net/'
3821  call check( nf90_put_att(ncid, nf90_global, 'references', trim(buffer)), &
3822  thisroutine )
3823 
3824 ! ------ Definition of the dimensions
3825 
3826  call set_grads_nc_tweaks(grads_nc_tweaks)
3827 
3828  if (grads_nc_tweaks) then
3829  call check( nf90_def_dim(ncid, 'x', 1, ncd), thisroutine )
3830  call check( nf90_def_dim(ncid, 'y', 1, ncd), thisroutine )
3831  end if
3832 
3833  call check( nf90_def_dim(ncid, 't', nf90_unlimited, ncd), thisroutine )
3834 
3835 ! ------ Definition of the variables
3836 
3837  if (grads_nc_tweaks) then
3838 
3839 ! ---- x
3840 
3841  call check( nf90_inq_dimid(ncid, 'x', nc1d), thisroutine )
3842  call check( nf90_def_var(ncid, 'x', nf90_float, nc1d, ncv), &
3843  thisroutine )
3844  buffer = 'm'
3845  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
3846  thisroutine )
3847  buffer = 'x'
3848  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
3849  thisroutine )
3850  buffer = 'Dummy x-coordinate for one point'
3851  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
3852  thisroutine )
3853  call check( nf90_put_att(ncid, ncv, 'axis', 'x'), thisroutine )
3854 
3855 ! ---- y
3856 
3857  call check( nf90_inq_dimid(ncid, 'y', nc1d), thisroutine )
3858  call check( nf90_def_var(ncid, 'y', nf90_float, nc1d, ncv), &
3859  thisroutine )
3860  buffer = 'm'
3861  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
3862  thisroutine )
3863  buffer = 'y'
3864  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
3865  thisroutine )
3866  buffer = 'Dummy y-coordinate for one point'
3867  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
3868  thisroutine )
3869  call check( nf90_put_att(ncid, ncv, 'axis', 'y'), thisroutine )
3870 
3871  end if
3872 
3873 ! ---- Time
3874 
3875  call check( nf90_inq_dimid(ncid, 't', nc1d), thisroutine )
3876  call check( nf90_def_var(ncid, 't', nf90_float, nc1d, ncv), &
3877  thisroutine )
3878  buffer = 'a'
3879  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
3880  thisroutine )
3881  buffer = 'time'
3882  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
3883  thisroutine )
3884  buffer = 'Time'
3885  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
3886  thisroutine )
3887  call check( nf90_put_att(ncid, ncv, 'axis', 't'), thisroutine )
3888 
3889  if (grads_nc_tweaks) then
3890 
3891 ! ---- Time offset
3892 
3893  call check( nf90_inq_dimid(ncid, 't', nc1d), thisroutine )
3894  call check( nf90_def_var(ncid, 't_add_offset', nf90_float, nc1d, ncv), &
3895  thisroutine )
3896  buffer = 'a'
3897  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
3898  thisroutine )
3899  buffer = 'time_add_offset'
3900  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
3901  thisroutine )
3902  buffer = 'Time offset'
3903  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
3904  thisroutine )
3905 
3906  time_add_offset_val = min(time_val, 0.0_dp)
3907 
3908  end if
3909 
3910  if (forcing_flag == 1) then
3911 
3912 ! ---- delta_ts
3913 
3914  call check( nf90_inq_dimid(ncid, 't', nc1d), thisroutine )
3915  call check( nf90_def_var(ncid, 'delta_ts', nf90_float, nc1d, ncv), &
3916  thisroutine )
3917  buffer = 'degC'
3918  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
3919  thisroutine )
3920  buffer = 'surface_temperature_anomaly'
3921  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
3922  thisroutine )
3923  buffer = 'Surface temperature anomaly'
3924  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
3925  thisroutine )
3926 
3927  else if ((forcing_flag == 2).or.(forcing_flag == 3)) then
3928 
3929 ! ---- glac_index
3930 
3931  call check( nf90_inq_dimid(ncid, 't', nc1d), thisroutine )
3932  call check( nf90_def_var(ncid, 'glac_index', nf90_float, nc1d, ncv), &
3933  thisroutine )
3934  buffer = '1'
3935  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
3936  thisroutine )
3937  buffer = 'glacial_index'
3938  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
3939  thisroutine )
3940  buffer = 'Glacial index'
3941  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
3942  thisroutine )
3943 
3944  end if
3945 
3946 ! ---- z_sl
3947 
3948  call check( nf90_inq_dimid(ncid, 't', nc1d), thisroutine )
3949  call check( nf90_def_var(ncid, 'z_sl', nf90_float, nc1d, ncv), &
3950  thisroutine )
3951  buffer = 'm'
3952  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
3953  thisroutine )
3954  buffer = 'global_average_sea_level_change'
3955  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
3956  thisroutine )
3957  buffer = 'Sea level'
3958  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
3959  thisroutine )
3960 
3961 ! ---- V_tot
3962 
3963  call check( nf90_inq_dimid(ncid, 't', nc1d), thisroutine )
3964  call check( nf90_def_var(ncid, 'V_tot', nf90_float, nc1d, ncv), &
3965  thisroutine )
3966  buffer = 'm3'
3967  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
3968  thisroutine )
3969  buffer = 'land_ice_volume'
3970  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
3971  thisroutine )
3972  buffer = 'Ice volume'
3973  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
3974  thisroutine )
3975 
3976 ! ---- V_grounded
3977 
3978  call check( nf90_inq_dimid(ncid, 't', nc1d), thisroutine )
3979  call check( nf90_def_var(ncid, 'V_grounded', nf90_float, nc1d, ncv), &
3980  thisroutine )
3981  buffer = 'm3'
3982  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
3983  thisroutine )
3984  buffer = 'grounded_land_ice_volume'
3985  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
3986  thisroutine )
3987  buffer = 'Volume of grounded ice'
3988  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
3989  thisroutine )
3990 
3991 ! ---- V_floating
3992 
3993  call check( nf90_inq_dimid(ncid, 't', nc1d), thisroutine )
3994  call check( nf90_def_var(ncid, 'V_floating', nf90_float, nc1d, ncv), &
3995  thisroutine )
3996  buffer = 'm3'
3997  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
3998  thisroutine )
3999  buffer = 'floating_ice_shelf_volume'
4000  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
4001  thisroutine )
4002  buffer = 'Volume of floating ice'
4003  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
4004  thisroutine )
4005 
4006 ! ---- A_tot
4007 
4008  call check( nf90_inq_dimid(ncid, 't', nc1d), thisroutine )
4009  call check( nf90_def_var(ncid, 'A_tot', nf90_float, nc1d, ncv), &
4010  thisroutine )
4011  buffer = 'm2'
4012  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
4013  thisroutine )
4014  buffer = 'land_ice_area'
4015  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
4016  thisroutine )
4017  buffer = 'Ice area'
4018  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
4019  thisroutine )
4020 
4021 ! ---- A_grounded
4022 
4023  call check( nf90_inq_dimid(ncid, 't', nc1d), thisroutine )
4024  call check( nf90_def_var(ncid, 'A_grounded', nf90_float, nc1d, ncv), &
4025  thisroutine )
4026  buffer = 'm2'
4027  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
4028  thisroutine )
4029  buffer = 'grounded_land_ice_area'
4030  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
4031  thisroutine )
4032  buffer = 'Area covered by grounded ice'
4033  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
4034  thisroutine )
4035 
4036 ! ---- A_floating
4037 
4038  call check( nf90_inq_dimid(ncid, 't', nc1d), thisroutine )
4039  call check( nf90_def_var(ncid, 'A_floating', nf90_float, nc1d, ncv), &
4040  thisroutine )
4041  buffer = 'm2'
4042  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
4043  thisroutine )
4044  buffer = 'floating_ice_shelf_area'
4045  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
4046  thisroutine )
4047  buffer = 'Area covered by floating ice'
4048  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
4049  thisroutine )
4050 
4051 ! ---- V_af
4052 
4053  call check( nf90_inq_dimid(ncid, 't', nc1d), thisroutine )
4054  call check( nf90_def_var(ncid, 'V_af', nf90_float, nc1d, ncv), &
4055  thisroutine )
4056  buffer = 'm3'
4057  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
4058  thisroutine )
4059  buffer = 'land_ice_volume_not_displacing_sea_water'
4060  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
4061  thisroutine )
4062  buffer = 'Ice volume above flotation'
4063  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
4064  thisroutine )
4065 
4066 ! ---- V_sle
4067 
4068  call check( nf90_inq_dimid(ncid, 't', nc1d), thisroutine )
4069  call check( nf90_def_var(ncid, 'V_sle', nf90_float, nc1d, ncv), &
4070  thisroutine )
4071  buffer = 'm SLE'
4072  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
4073  thisroutine )
4074  buffer = 'land_ice_volume_sle'
4075  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
4076  thisroutine )
4077  buffer = 'Ice volume in SLE'
4078  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
4079  thisroutine )
4080 
4081 ! ---- V_temp
4082 
4083  call check( nf90_inq_dimid(ncid, 't', nc1d), thisroutine )
4084  call check( nf90_def_var(ncid, 'V_temp', nf90_float, nc1d, ncv), &
4085  thisroutine )
4086  buffer = 'm3'
4087  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
4088  thisroutine )
4089  buffer = 'temperate_land_ice_volume'
4090  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
4091  thisroutine )
4092  buffer = 'Volume of temperate ice'
4093  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
4094  thisroutine )
4095 
4096 ! ---- A_temp
4097 
4098  call check( nf90_inq_dimid(ncid, 't', nc1d), thisroutine )
4099  call check( nf90_def_var(ncid, 'A_temp', nf90_float, nc1d, ncv), &
4100  thisroutine )
4101  buffer = 'm2'
4102  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
4103  thisroutine )
4104  buffer = 'temperate_land_ice_area'
4105  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
4106  thisroutine )
4107  buffer = 'Area covered by temperate ice'
4108  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
4109  thisroutine )
4110 
4111 ! ---- Q_s
4112 
4113  call check( nf90_inq_dimid(ncid, 't', nc1d), thisroutine )
4114  call check( nf90_def_var(ncid, 'Q_s', nf90_float, nc1d, ncv), &
4115  thisroutine )
4116  buffer = 'm3 ice equiv. a-1'
4117  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
4118  thisroutine )
4119  buffer = 'tendency_of_land_ice_volume_due_to_surface_mass_balance'
4120  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
4121  thisroutine )
4122  buffer = 'Total surface mass balance'
4123  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
4124  thisroutine )
4125 
4126 ! ---- precip_tot
4127 
4128  call check( nf90_inq_dimid(ncid, 't', nc1d), thisroutine )
4129  call check( nf90_def_var(ncid, 'precip_tot', nf90_float, nc1d, ncv), &
4130  thisroutine )
4131  buffer = 'm3 ice equiv. a-1'
4132  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
4133  thisroutine )
4134  buffer = 'precipitation_rate'
4135  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
4136  thisroutine )
4137  buffer = 'Precipitation rate'
4138  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
4139  thisroutine )
4140 
4141 ! ---- runoff_tot
4142 
4143  call check( nf90_inq_dimid(ncid, 't', nc1d), thisroutine )
4144  call check( nf90_def_var(ncid, 'runoff_tot', nf90_float, nc1d, ncv), &
4145  thisroutine )
4146  buffer = 'm3 ice equiv. a-1'
4147  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
4148  thisroutine )
4149  buffer = 'runoff_rate'
4150  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
4151  thisroutine )
4152  buffer = 'Runoff rate'
4153  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
4154  thisroutine )
4155 
4156 ! ---- bmb_tot
4157 
4158  call check( nf90_inq_dimid(ncid, 't', nc1d), thisroutine )
4159  call check( nf90_def_var(ncid, 'bmb_tot', nf90_float, nc1d, ncv), &
4160  thisroutine )
4161  buffer = 'm3 ice equiv. a-1'
4162  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
4163  thisroutine )
4164  buffer = 'tendency_of_land_ice_volume_due_to_basal_mass_balance'
4165  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
4166  thisroutine )
4167  buffer = 'Total basal mass balance'
4168  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
4169  thisroutine )
4170 
4171 ! ---- Q_b
4172 
4173  call check( nf90_inq_dimid(ncid, 't', nc1d), thisroutine )
4174  call check( nf90_def_var(ncid, 'Q_b', nf90_float, nc1d, ncv), &
4175  thisroutine )
4176  buffer = 'm3 ice equiv. a-1'
4177  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
4178  thisroutine )
4179  buffer = 'basal_melting_rate'
4180  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
4181  thisroutine )
4182  buffer = 'Basal melting rate'
4183  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
4184  thisroutine )
4185 
4186 ! ---- Q_temp
4187 
4188  call check( nf90_inq_dimid(ncid, 't', nc1d), thisroutine )
4189  call check( nf90_def_var(ncid, 'Q_temp', nf90_float, nc1d, ncv), &
4190  thisroutine )
4191  buffer = 'm3 ice equiv. a-1'
4192  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
4193  thisroutine )
4194  buffer = 'temperate_layer_drainage_rate'
4195  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
4196  thisroutine )
4197  buffer = 'Drainage rate from the temperate ice layer'
4198  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
4199  thisroutine )
4200 
4201 ! ---- calv_tot
4202 
4203  call check( nf90_inq_dimid(ncid, 't', nc1d), thisroutine )
4204  call check( nf90_def_var(ncid, 'calv_tot', nf90_float, nc1d, ncv), &
4205  thisroutine )
4206  buffer = 'm3 ice equiv. a-1'
4207  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
4208  thisroutine )
4209  buffer = 'tendency_of_land_ice_volume_due_to_calving'
4210  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
4211  thisroutine )
4212  buffer = 'Total calving rate'
4213  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
4214  thisroutine )
4215 
4216 #if (DISC>0)
4217 
4218 ! ---- disc_lsc
4219 
4220  call check( nf90_inq_dimid(ncid, 't', nc1d), thisroutine )
4221  call check( nf90_def_var(ncid, 'disc_lsc', nf90_float, nc1d, ncv), &
4222  thisroutine )
4223  buffer = 'm3 ice equiv. a-1'
4224  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
4225  thisroutine )
4226  buffer = 'large_scale_ice_lost_into_the_ocean'
4227  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
4228  thisroutine )
4229  buffer = 'Large scale ice lost into the ocean'
4230  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
4231  thisroutine )
4232 
4233 ! ---- disc_ssc
4234 
4235  call check( nf90_inq_dimid(ncid, 't', nc1d), thisroutine )
4236  call check( nf90_def_var(ncid, 'disc_ssc', nf90_float, nc1d, ncv), &
4237  thisroutine )
4238  buffer = 'm3 ice equiv. a-1'
4239  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
4240  thisroutine )
4241  buffer = 'small_scale_ice_lost_into_the_ocean'
4242  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
4243  thisroutine )
4244  buffer = 'Small scale ice lost into the ocean'
4245  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
4246  thisroutine )
4247 
4248 ! ---- dT_glann
4249 
4250  call check( nf90_inq_dimid(ncid, 't', nc1d), thisroutine )
4251  call check( nf90_def_var(ncid, 'dT_glann', nf90_float, nc1d, ncv), &
4252  thisroutine )
4253  buffer = 'degC'
4254  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
4255  thisroutine )
4256  buffer = 'global_annual_temperature_anomaly'
4257  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
4258  thisroutine )
4259  buffer = 'Global annual temperature anomaly'
4260  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
4261  thisroutine )
4262 
4263 ! ---- dT_sub
4264 
4265  call check( nf90_inq_dimid(ncid, 't', nc1d), thisroutine )
4266  call check( nf90_def_var(ncid, 'dT_sub', nf90_float, nc1d, ncv), &
4267  thisroutine )
4268  buffer = 'degC'
4269  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
4270  thisroutine )
4271  buffer = 'subsurface_ocean_temperature_anomaly'
4272  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
4273  thisroutine )
4274  buffer = 'Subsurface ocean temperature anomaly'
4275  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
4276  thisroutine )
4277 
4278 #endif
4279 
4280 ! ---- dV_dt
4281 
4282  call check( nf90_inq_dimid(ncid, 't', nc1d), thisroutine )
4283  call check( nf90_def_var(ncid, 'dV_dt', nf90_float, nc1d, ncv), &
4284  thisroutine )
4285  buffer = 'm3 a-1'
4286  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
4287  thisroutine )
4288  buffer = 'tendency_of_land_ice_volume'
4289  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
4290  thisroutine )
4291  buffer = 'Rate of ice volume change'
4292  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
4293  thisroutine )
4294 
4295 ! ---- mb_resid
4296 
4297  call check( nf90_inq_dimid(ncid, 't', nc1d), thisroutine )
4298  call check( nf90_def_var(ncid, 'mb_resid', nf90_float, nc1d, ncv), &
4299  thisroutine )
4300  buffer = 'm3 ice equiv. a-1'
4301  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
4302  thisroutine )
4303  buffer = 'total_mass_balance_residual'
4304  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
4305  thisroutine )
4306  buffer = 'Residual of the total mass balance'
4307  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
4308  thisroutine )
4309 
4310 ! ---- mbp
4311 
4312  call check( nf90_inq_dimid(ncid, 't', nc1d), thisroutine )
4313  call check( nf90_def_var(ncid, 'mbp', nf90_float, nc1d, ncv), &
4314  thisroutine )
4315  buffer = '1'
4316  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
4317  thisroutine )
4318  buffer = 'mass_balance_partition'
4319  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
4320  thisroutine )
4321  buffer = 'Mass balance partition'
4322  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
4323  thisroutine )
4324 
4325 ! ---- H_max
4326 
4327  call check( nf90_inq_dimid(ncid, 't', nc1d), thisroutine )
4328  call check( nf90_def_var(ncid, 'H_max', nf90_float, nc1d, ncv), &
4329  thisroutine )
4330  buffer = 'm'
4331  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
4332  thisroutine )
4333  buffer = 'maximum_ice_thickness'
4334  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
4335  thisroutine )
4336  buffer = 'Maximum ice thickness'
4337  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
4338  thisroutine )
4339 
4340 ! ---- H_t_max
4341 
4342  call check( nf90_inq_dimid(ncid, 't', nc1d), thisroutine )
4343  call check( nf90_def_var(ncid, 'H_t_max', nf90_float, nc1d, ncv), &
4344  thisroutine )
4345  buffer = 'm'
4346  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
4347  thisroutine )
4348  buffer = 'maximum_thickness_of_temperate_ice'
4349  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
4350  thisroutine )
4351  buffer = 'Maximum thickness of temperate ice'
4352  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
4353  thisroutine )
4354 
4355 ! ---- zs_max
4356 
4357  call check( nf90_inq_dimid(ncid, 't', nc1d), thisroutine )
4358  call check( nf90_def_var(ncid, 'zs_max', nf90_float, nc1d, ncv), &
4359  thisroutine )
4360  buffer = 'm'
4361  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
4362  thisroutine )
4363  buffer = 'maximum_surface_elevation'
4364  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
4365  thisroutine )
4366  buffer = 'Maximum surface elevation'
4367  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
4368  thisroutine )
4369 
4370 ! ---- vs_max
4371 
4372  call check( nf90_inq_dimid(ncid, 't', nc1d), thisroutine )
4373  call check( nf90_def_var(ncid, 'vs_max', nf90_float, nc1d, ncv), &
4374  thisroutine )
4375  buffer = 'm a-1'
4376  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
4377  thisroutine )
4378  buffer = 'maximum_surface_speed'
4379  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
4380  thisroutine )
4381  buffer = 'Maximum surface speed'
4382  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
4383  thisroutine )
4384 
4385 ! ---- Tbh_max
4386 
4387  call check( nf90_inq_dimid(ncid, 't', nc1d), thisroutine )
4388  call check( nf90_def_var(ncid, 'Tbh_max', nf90_float, nc1d, ncv), &
4389  thisroutine )
4390  buffer = 'degC'
4391  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
4392  thisroutine )
4393  buffer = 'maximum_basal_temperature_relative_to_pmp'
4394  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
4395  thisroutine )
4396  buffer = 'Maximum basal temperature relative to pmp'
4397  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
4398  thisroutine )
4399 
4400 ! ---- End of the definitions
4401 
4402  call check( nf90_enddef(ncid), thisroutine )
4403 
4404 end if
4405 
4406 ! ------ Computation of further data
4407 
4408 ! Q_s, Q_b, Q_temp from earlier computations (now in ice equiv.)
4409 q_s = q_s *(rho_w/rho)
4410  ! m^3/a water equiv. -> m^3/a ice equiv.
4411 q_b = q_b *(rho_w/rho)
4412  ! m^3/a water equiv. -> m^3/a ice equiv.
4413 q_temp = q_temp *(rho_w/rho)
4414  ! m^3/a water equiv. -> m^3/a ice equiv.
4415 
4416 ! Computation of volume fluxes
4417 
4418 dv_dt = 0.0_dp
4419 precip_tot = 0.0_dp
4420 
4421 do i=1, imax-1
4422 do j=1, jmax-1
4423 
4424  if ( (maske(j,i)==0_i2b).or.(maske_old(j,i)==0_i2b) &
4425  .or.(maske(j,i)==3_i2b).or.(maske_old(j,i)==3_i2b) ) then
4426  ! grounded or floating ice before or after the time step
4427 
4428  ! Volume change
4429  dv_dt = dv_dt + (dzs_dtau(j,i)-dzb_dtau(j,i))*area(j,i)
4430  precip_tot = precip_tot + accum(j,i)*area(j,i)
4431 
4432  end if
4433 
4434 end do
4435 end do
4436 
4437 precip_tot = precip_tot *year_sec
4438  ! m^3/s ice equiv. -> m^3/a ice equiv.
4439 dv_dt = dv_dt *year_sec
4440  ! m^3/s ice equiv. -> m^3/a ice equiv.
4441 
4442 #if (MB_ACCOUNT==0)
4443 
4444 runoff_tot = 0.0_dp
4445 calv_tot = 0.0_dp
4446 
4447 do i=1, imax-1
4448 do j=1, jmax-1
4449 
4450  if ( (maske(j,i)==0_i2b).or.(maske_old(j,i)==0_i2b) &
4451  .or.(maske(j,i)==3_i2b).or.(maske_old(j,i)==3_i2b) ) then
4452  ! grounded or floating ice before or after the time step
4453 
4454  runoff_tot = runoff_tot + runoff(j,i)*area(j,i)
4455  calv_tot = calv_tot + calv_grounded(j,i)*area(j,i)
4456 
4457  end if
4458 
4459 end do
4460 end do
4461 
4462 bmb_tot = -(q_b+q_temp) ! positive for supply, negative for loss
4463 
4464 #elif (MB_ACCOUNT==1)
4465 
4466 ! MB: total mass balance as computed in subroutine apply_smb
4467 ! LMT: total lost on land at the top of the ice sheet
4468 ! LMB: total lost on land at the base of the ice sheet
4469 ! LMH: total hidden lost on land
4470 ! OMT: total lost in the ocean at the top of the ice sheet
4471 ! OMB: total lost in the ocean at the base of the ice sheet
4472 ! OMH: total hidden lost in the ocean
4473 ! ODT: total lost through ice discharge parameterization
4474 ! ODH: total hidden lost through ice discharge parameterization
4475 
4476 mb = 0.0_dp
4477 lmt = 0.0_dp; omt = 0.0_dp
4478 lmb = 0.0_dp; omb = 0.0_dp
4479 lmh = 0.0_dp; omh = 0.0_dp
4480 odt = 0.0_dp; odh = 0.0_dp
4481 
4482 do i=1, imax-1
4483 do j=1, jmax-1
4484 
4485  if ( mask_run(j,i) .ne. 0 ) then
4486  ! glaciated land and ocean (including hidden melt points)
4487 
4488  ! Quantify what types of melt occurred
4489  select case ( mask_run(j,i) )
4490  case( 2 )
4491  omt = omt + runoff_top(j,i) * area(j,i)
4492  odt = odt + disc_top(j,i) * area(j,i)
4493  omb = omb + runoff_bot(j,i) * area(j,i)
4494  case( 1 )
4495  lmt = lmt + runoff_top(j,i) * area(j,i)
4496  odt = odt + disc_top(j,i) * area(j,i) ! ok
4497  lmb = lmb + runoff_bot(j,i) * area(j,i)
4498  case( -1 )
4499  odh = odh + disc_top(j,i) * area(j,i)
4500  lmh = lmh + runoff_top(j,i) * area(j,i)
4501 ! note: runoff_bot(j,i) not counted as hidden yet
4502  case( -2 )
4503  odh = odh + disc_top(j,i) * area(j,i)
4504  omh = omh + runoff_top(j,i) * area(j,i)
4505 ! note: runoff_bot(j,i) not counted as hidden yet
4506  end select
4507 
4508  end if
4509 
4510  ! Actual ice mass balance (from top melt, bottom melt and calving)
4511  mb = mb + mb_source_apl(j,i)*area(j,i)
4512 
4513 end do
4514 end do
4515 
4516 ! Runoff on land (excluding basal melt)
4517 runoff_tot = lmt + lmh
4518 ! Ice discharge (excluding basal melt)
4519 calv_tot = omt + omh + odt + odh
4520 ! Ice discharge from ice flow, large scale
4521 disc_lsc = omt + omh
4522 ! Ice discharge from parameterization, small scale
4523 disc_ssc = odt + odh
4524 ! Basal mass balance (here basal melt only)
4525 bmb_tot = -lmb-omb
4526 
4527 mb = mb * year_sec
4528  ! m^3/s ice equiv. -> m^3/a ice equiv.
4529 disc_lsc = disc_lsc * year_sec
4530  ! m^3/s ice equiv. -> m^3/a ice equiv.
4531 disc_ssc = disc_ssc * year_sec
4532  ! m^3/s ice equiv. -> m^3/a ice equiv.
4533 bmb_tot = bmb_tot * year_sec
4534  ! m^3/s ice equiv. -> m^3/a ice equiv.
4535 
4536 #endif /* if (MB_ACCOUNT==0) elif (MB_ACCOUNT==1) */
4537 
4538 runoff_tot = runoff_tot *year_sec
4539  ! m^3/s ice equiv. -> m^3/a ice equiv.
4540 calv_tot = calv_tot * year_sec
4541  ! m^3/s ice equiv. -> m^3/a ice equiv.
4542 
4543 if (precip_tot.ne.0.0_dp) then
4544  mbp = calv_tot/precip_tot
4545 else
4546  mbp = 0.0_dp
4547 end if
4548 
4549 #if (MB_ACCOUNT==0)
4550 mb_resid = q_s + bmb_tot - calv_tot - dv_dt
4551 #elif(MB_ACCOUNT==1)
4552 mb_resid = mb - dv_dt
4553 #endif
4554 
4555 ! ------ Writing of data on NetCDF file
4556 
4557 if (firstcall.and.grads_nc_tweaks) then
4558 
4559  nc1cor = (/ 1 /)
4560 
4561  call check( nf90_inq_varid(ncid, 'x', ncv), thisroutine )
4562  call check( nf90_put_var(ncid, ncv, 0.0_sp, start=nc1cor), thisroutine )
4563 
4564  call check( nf90_inq_varid(ncid, 'y', ncv), thisroutine )
4565  call check( nf90_put_var(ncid, ncv, 0.0_sp, start=nc1cor), thisroutine )
4566 
4567  call check( nf90_sync(ncid), thisroutine )
4568 
4569 end if
4570 
4571 nc1cor(1) = counter
4572 ! nc1cnt(1) = 1 ! (not needed, and causes troubles for whatever reason)
4573 
4574 call check( nf90_inq_varid(ncid, 't', ncv), thisroutine )
4575 
4576 if (.not.grads_nc_tweaks) then
4577  call check( nf90_put_var(ncid, ncv, real(time_val,sp), &
4578  start=nc1cor), thisroutine )
4579 else
4580  call check( nf90_put_var(ncid, ncv, &
4581  real(time_val-time_add_offset_val,sp), &
4582  start=nc1cor), thisroutine )
4583  ! this makes sure that all times are >=0
4584  ! (GrADS doesn't like negative numbers)
4585  call check( nf90_inq_varid(ncid, 't_add_offset', ncv), thisroutine )
4586  call check( nf90_put_var(ncid, ncv, &
4587  real(time_add_offset_val,sp), &
4588  start=nc1cor), thisroutine )
4589 end if
4590 
4591 if (forcing_flag == 1) then
4592 
4593  call check( nf90_inq_varid(ncid, 'delta_ts', ncv), thisroutine )
4594  call check( nf90_put_var(ncid, ncv, real(delta_ts,sp), &
4595  start=nc1cor), thisroutine )
4596 
4597 else if (forcing_flag == 2) then
4598 
4599  call check( nf90_inq_varid(ncid, 'glac_index', ncv), thisroutine )
4600  call check( nf90_put_var(ncid, ncv, real(glac_index,sp), &
4601  start=nc1cor), thisroutine )
4602 
4603 else if (forcing_flag == 3) then
4604 
4605  call check( nf90_inq_varid(ncid, 'glac_index', ncv), thisroutine )
4606  call check( nf90_put_var(ncid, ncv, real(no_value_neg_2,sp), &
4607  start=nc1cor), thisroutine )
4608 
4609 end if
4610 
4611 call check( nf90_inq_varid(ncid, 'z_sl', ncv), thisroutine )
4612 call check( nf90_put_var(ncid, ncv, real(z_sl,sp), &
4613  start=nc1cor), thisroutine )
4614 
4615 call check( nf90_inq_varid(ncid, 'V_tot', ncv), thisroutine )
4616 call check( nf90_put_var(ncid, ncv, real(V_tot,sp), &
4617  start=nc1cor), thisroutine )
4618 
4619 call check( nf90_inq_varid(ncid, 'V_grounded', ncv), thisroutine )
4620 call check( nf90_put_var(ncid, ncv, real(V_grounded,sp), &
4621  start=nc1cor), thisroutine )
4622 
4623 call check( nf90_inq_varid(ncid, 'V_floating', ncv), thisroutine )
4624 call check( nf90_put_var(ncid, ncv, real(V_floating,sp), &
4625  start=nc1cor), thisroutine )
4626 
4627 call check( nf90_inq_varid(ncid, 'A_tot', ncv), thisroutine )
4628 call check( nf90_put_var(ncid, ncv, real(A_tot,sp), &
4629  start=nc1cor), thisroutine )
4630 
4631 call check( nf90_inq_varid(ncid, 'A_grounded', ncv), thisroutine )
4632 call check( nf90_put_var(ncid, ncv, real(A_grounded,sp), &
4633  start=nc1cor), thisroutine )
4634 
4635 call check( nf90_inq_varid(ncid, 'A_floating', ncv), thisroutine )
4636 call check( nf90_put_var(ncid, ncv, real(A_floating,sp), &
4637  start=nc1cor), thisroutine )
4638 
4639 call check( nf90_inq_varid(ncid, 'V_af', ncv), thisroutine )
4640 call check( nf90_put_var(ncid, ncv, real(V_af,sp), &
4641  start=nc1cor), thisroutine )
4642 
4643 call check( nf90_inq_varid(ncid, 'V_sle', ncv), thisroutine )
4644 call check( nf90_put_var(ncid, ncv, real(V_sle,sp), &
4645  start=nc1cor), thisroutine )
4646 
4647 call check( nf90_inq_varid(ncid, 'V_temp', ncv), thisroutine )
4648 call check( nf90_put_var(ncid, ncv, real(V_temp,sp), &
4649  start=nc1cor), thisroutine )
4650 
4651 call check( nf90_inq_varid(ncid, 'A_temp', ncv), thisroutine )
4652 call check( nf90_put_var(ncid, ncv, real(A_temp,sp), &
4653  start=nc1cor), thisroutine )
4654 
4655 call check( nf90_inq_varid(ncid, 'Q_s', ncv), thisroutine )
4656 call check( nf90_put_var(ncid, ncv, real(Q_s,sp), &
4657  start=nc1cor), thisroutine )
4658 
4659 call check( nf90_inq_varid(ncid, 'precip_tot', ncv), thisroutine )
4660 call check( nf90_put_var(ncid, ncv, real(precip_tot,sp), &
4661  start=nc1cor), thisroutine )
4662 
4663 call check( nf90_inq_varid(ncid, 'runoff_tot', ncv), thisroutine )
4664 call check( nf90_put_var(ncid, ncv, real(runoff_tot,sp), &
4665  start=nc1cor), thisroutine )
4666 
4667 call check( nf90_inq_varid(ncid, 'bmb_tot', ncv), thisroutine )
4668 call check( nf90_put_var(ncid, ncv, real(bmb_tot,sp), &
4669  start=nc1cor), thisroutine )
4670 
4671 call check( nf90_inq_varid(ncid, 'Q_b', ncv), thisroutine )
4672 call check( nf90_put_var(ncid, ncv, real(Q_b,sp), &
4673  start=nc1cor), thisroutine )
4674 
4675 call check( nf90_inq_varid(ncid, 'Q_temp', ncv), thisroutine )
4676 call check( nf90_put_var(ncid, ncv, real(Q_temp,sp), &
4677  start=nc1cor), thisroutine )
4678 
4679 call check( nf90_inq_varid(ncid, 'calv_tot', ncv), thisroutine )
4680 call check( nf90_put_var(ncid, ncv, real(calv_tot,sp), &
4681  start=nc1cor), thisroutine )
4682 
4683 #if (DISC>0)
4684 call check( nf90_inq_varid(ncid, 'disc_lsc', ncv), thisroutine )
4685 call check( nf90_put_var(ncid, ncv, real(disc_lsc,sp), &
4686  start=nc1cor), thisroutine )
4687 
4688 call check( nf90_inq_varid(ncid, 'disc_ssc', ncv), thisroutine )
4689 call check( nf90_put_var(ncid, ncv, real(disc_ssc,sp), &
4690  start=nc1cor), thisroutine )
4691 
4692 call check( nf90_inq_varid(ncid, 'dT_glann', ncv), thisroutine )
4693 call check( nf90_put_var(ncid, ncv, real(dT_glann,sp), &
4694  start=nc1cor), thisroutine )
4695 
4696 call check( nf90_inq_varid(ncid, 'dT_sub', ncv), thisroutine )
4697 call check( nf90_put_var(ncid, ncv, real(dT_sub,sp), &
4698  start=nc1cor), thisroutine )
4699 #endif
4700 
4701 call check( nf90_inq_varid(ncid, 'dV_dt', ncv), thisroutine )
4702 call check( nf90_put_var(ncid, ncv, real(dV_dt,sp), &
4703  start=nc1cor), thisroutine )
4704 
4705 call check( nf90_inq_varid(ncid, 'mb_resid', ncv), thisroutine )
4706 call check( nf90_put_var(ncid, ncv, real(mb_resid,sp), &
4707  start=nc1cor), thisroutine )
4708 
4709 call check( nf90_inq_varid(ncid, 'mbp', ncv), thisroutine )
4710 call check( nf90_put_var(ncid, ncv, real(mbp,sp), &
4711  start=nc1cor), thisroutine )
4712 
4713 call check( nf90_inq_varid(ncid, 'H_max', ncv), thisroutine )
4714 call check( nf90_put_var(ncid, ncv, real(H_max,sp), &
4715  start=nc1cor), thisroutine )
4716 
4717 call check( nf90_inq_varid(ncid, 'H_t_max', ncv), thisroutine )
4718 call check( nf90_put_var(ncid, ncv, real(H_t_max,sp), &
4719  start=nc1cor), thisroutine )
4720 
4721 call check( nf90_inq_varid(ncid, 'zs_max', ncv), thisroutine )
4722 call check( nf90_put_var(ncid, ncv, real(zs_max,sp), &
4723  start=nc1cor), thisroutine )
4724 
4725 call check( nf90_inq_varid(ncid, 'vs_max', ncv), thisroutine )
4726 call check( nf90_put_var(ncid, ncv, real(vs_max,sp), &
4727  start=nc1cor), thisroutine )
4728 
4729 call check( nf90_inq_varid(ncid, 'Tbh_max', ncv), thisroutine )
4730 call check( nf90_put_var(ncid, ncv, real(Tbh_max,sp), &
4731  start=nc1cor), thisroutine )
4732 
4733 ! ------ Syncing NetCDF file (every 100th time)
4734 
4735 n_sync = 100
4736 
4737 if ( mod((counter-1), n_sync) == 0 ) &
4738  call check( nf90_sync(ncid), thisroutine )
4739 
4740 #endif
4741 
4742 if (firstcall) firstcall = .false.
4743 
4744 end subroutine output2
4745 
4746 !-------------------------------------------------------------------------------
4747 !> Writing of time-series data of the deep ice cores on file in ASCII format
4748 !! (and optionally in NetCDF format).
4749 !<------------------------------------------------------------------------------
4750 subroutine output4(time, dxi, deta, delta_ts, glac_index, z_sl)
4752 #if (NETCDF>1)
4753  use netcdf
4754  use nc_check_m
4755 #endif
4756 
4757 implicit none
4758 
4759 real(dp), intent(in) :: time, dxi, deta, delta_ts, glac_index, z_sl
4760 
4761 integer(i4b) :: i, j, n
4762 real(sp), dimension(:), allocatable :: r_n_core
4763 real(dp) :: time_val
4764 real(dp), dimension(0:JMAX,0:IMAX) :: field
4765 real(dp), dimension(:), allocatable :: H_core, temp_b_core, &
4766  vx_b_core, vy_b_core, vh_b_core, &
4767  vx_s_core, vy_s_core, vh_s_core, &
4768  Rx_b_core, Ry_b_core, R_b_core, &
4769  bmb_core
4770 
4771 #if (NETCDF>1)
4772 integer(i4b), save :: ncid
4773 integer(i4b) :: ncd, ncv, nc1d, nc2d(2)
4774 integer(i4b) :: nc1cor(1), nc1cnt(1), nc2cor(2), nc2cnt(2)
4775 integer(i4b) :: n_sync
4776 real(dp), save :: time_add_offset_val
4777 character(len= 16) :: ch_date, ch_time, ch_zone
4778 character(len=256) :: filename, filename_with_path, buffer
4779 logical, save :: grads_nc_tweaks
4780 #endif
4781 
4782 integer(i4b), save :: counter = 0
4783 logical, save :: firstcall = .true.
4784 
4785 character(len=64), parameter :: thisroutine = 'output4'
4786 
4787 counter = counter + 1
4788 
4789 if (n_core >= 1) then
4790 
4791  allocate(r_n_core(n_core), h_core(n_core), temp_b_core(n_core), &
4792  vx_b_core(n_core), vy_b_core(n_core), vh_b_core(n_core), &
4793  vx_s_core(n_core), vy_s_core(n_core), vh_s_core(n_core), &
4794  rx_b_core(n_core), ry_b_core(n_core), r_b_core(n_core), &
4795  bmb_core(n_core))
4796 
4797 !-------- Determination of ice-core data --------
4798 
4799  do n=1, n_core
4800 
4801 ! ------Ice core number
4802 
4803  r_n_core(n) = real(n,sp)
4804 
4805 ! ------ Ice thickness
4806 
4807  field = h_c + h_t
4808  call borehole(field, x_core(n), y_core(n), dxi, deta, 'grid', h_core(n))
4809 
4810 ! ------ Basal velocity
4811 
4812  do i=0, imax
4813  do j=0, jmax
4814  field(j,i) = vx_t(0,j,i)
4815  end do
4816  end do
4817 
4818  call borehole(field, x_core(n), y_core(n), dxi, deta, 'sg_x', vx_b_core(n))
4819 
4820  do i=0, imax
4821  do j=0, jmax
4822  field(j,i) = vy_t(0,j,i)
4823  end do
4824  end do
4825 
4826  call borehole(field, x_core(n), y_core(n), dxi, deta, 'sg_y', vy_b_core(n))
4827 
4828  vh_b_core(n) = sqrt(vx_b_core(n)**2+vy_b_core(n)**2)
4829 
4830 ! ------ Surface velocity
4831 
4832  do i=0, imax
4833  do j=0, jmax
4834  field(j,i) = vx_c(kcmax,j,i)
4835  end do
4836  end do
4837 
4838  call borehole(field, x_core(n), y_core(n), dxi, deta, 'sg_x', vx_s_core(n))
4839 
4840  do i=0, imax
4841  do j=0, jmax
4842  field(j,i) = vy_c(kcmax,j,i)
4843  end do
4844  end do
4845 
4846  call borehole(field, x_core(n), y_core(n), dxi, deta, 'sg_y', vy_s_core(n))
4847 
4848  vh_s_core(n) = sqrt(vx_s_core(n)**2+vy_s_core(n)**2)
4849 
4850 ! ------ Basal temperature
4851 
4852  do i=0, imax
4853  do j=0, jmax
4854  field(j,i) = temp_r(krmax,j,i)
4855  end do
4856  end do
4857 
4858  call borehole(field, x_core(n), y_core(n), dxi, deta, 'grid', temp_b_core(n))
4859 
4860 ! ------ Basal frictional heating
4861 
4862  do i=0, imax
4863  do j=0, jmax
4864  field(j,i) = vx_t(0,j,i)*txz_t(0,j,i)
4865  end do
4866  end do
4867 
4868  call borehole(field, x_core(n), y_core(n), dxi, deta, 'sg_x', rx_b_core(n))
4869 
4870  do i=0, imax
4871  do j=0, jmax
4872  field(j,i) = vy_t(0,j,i)*tyz_t(0,j,i)
4873  end do
4874  end do
4875 
4876  call borehole(field, x_core(n), y_core(n), dxi, deta, 'sg_y', ry_b_core(n))
4877 
4878  r_b_core(n) = rx_b_core(n) + ry_b_core(n)
4879 
4880 ! ------ Basal mass balance
4881 
4882  field = -(q_bm+q_tld) ! positive for supply, negative for loss
4883  call borehole(field, x_core(n), y_core(n), dxi, deta, 'grid', bmb_core(n))
4884 
4885  end do
4886 
4887 ! ------ Conversion
4888 
4889 #if (!defined(OUT_TIMES) || OUT_TIMES==1)
4890  time_val = time *sec_to_year ! s -> a
4891 #elif (OUT_TIMES==2)
4892  time_val = (time+year_zero) *sec_to_year ! s -> a
4893 #else
4894  stop ' >>> output4: OUT_TIMES must be either 1 or 2!'
4895 #endif
4896 
4897  vh_b_core = vh_b_core *year_sec ! m/s -> m/a
4898  vh_s_core = vh_s_core *year_sec ! m/s -> m/a
4899 
4900  bmb_core = bmb_core *year_sec ! m ice equiv./s -> m ice equiv./a
4901 
4902 !-------- Writing of data on file --------
4903 
4904  if (forcing_flag == 1) then
4905  write(unit=14, fmt='(1pe13.6,2(1pe13.4))') time_val, delta_ts, z_sl
4906  else if (forcing_flag == 2) then
4907  write(unit=14, fmt='(1pe13.6,2(1pe13.4))') time_val, glac_index, z_sl
4908  else if (forcing_flag == 3) then
4909  write(unit=14, fmt='(1pe13.6,2(1pe13.4))') time_val, no_value_pos_1, z_sl
4910  end if
4911 
4912  n=1
4913  write(unit=14, fmt='(13x,1pe13.4)', advance='no') h_core(n)
4914  do n=2, n_core-1
4915  write(unit=14, fmt='(1pe13.4)', advance='no') h_core(n)
4916  end do
4917  n=n_core
4918  write(unit=14, fmt='(1pe13.4)') h_core(n)
4919 
4920  n=1
4921  write(unit=14, fmt='(13x,1pe13.4)', advance='no') vh_s_core(n)
4922  do n=2, n_core-1
4923  write(unit=14, fmt='(1pe13.4)', advance='no') vh_s_core(n)
4924  end do
4925  n=n_core
4926  write(unit=14, fmt='(1pe13.4)') vh_s_core(n)
4927 
4928  n=1
4929  write(unit=14, fmt='(13x,1pe13.4)', advance='no') temp_b_core(n)
4930  do n=2, n_core-1
4931  write(unit=14, fmt='(1pe13.4)', advance='no') temp_b_core(n)
4932  end do
4933  n=n_core
4934  write(unit=14, fmt='(1pe13.4,/)') temp_b_core(n)
4935 
4936 !-------- Extended time-series file in NetCDF format --------
4937 
4938 #if (NETCDF>1)
4939 
4940  if (firstcall) then
4941 
4942 ! ------ Open NetCDF file
4943 
4944  filename = trim(runname)//'_core.nc'
4945  filename_with_path = trim(outpath)//'/'//trim(filename)
4946 
4947  call check( nf90_create(trim(filename_with_path), nf90_noclobber, ncid), &
4948  thisroutine )
4949 
4950  ncid_core = ncid
4951 
4952 ! ------ Global attributes
4953 
4954  buffer = 'Time-series output for the deep ice cores of simulation '// &
4955  trim(runname)
4956  call check( nf90_put_att(ncid, nf90_global, 'title', trim(buffer)), &
4957  thisroutine )
4958 
4959  buffer = 'Institute of Low Temperature Science, Hokkaido University, '// &
4960  'Sapporo, Japan'
4961  call check( nf90_put_att(ncid, nf90_global, 'institution', &
4962  trim(buffer)), &
4963  thisroutine )
4964 
4965  buffer = 'SICOPOLIS Version '//version
4966  call check( nf90_put_att(ncid, nf90_global, 'source', trim(buffer)), &
4967  thisroutine )
4968 
4969  call date_and_time(ch_date, ch_time, ch_zone)
4970  buffer = ch_date(1:4)//'-'//ch_date(5:6)//'-'//ch_date(7:8)//' '// &
4971  ch_time(1:2)//':'//ch_time(3:4)//':'//ch_time(5:6)//' '// &
4972  ch_zone(1:3)//':'//ch_zone(4:5)//' - Data produced'
4973  call check( nf90_put_att(ncid, nf90_global, 'history', trim(buffer)), &
4974  thisroutine )
4975 
4976  buffer = 'http://www.sicopolis.net/'
4977  call check( nf90_put_att(ncid, nf90_global, 'references', trim(buffer)), &
4978  thisroutine )
4979 
4980 ! ------ Definition of the dimensions
4981 
4982  call set_grads_nc_tweaks(grads_nc_tweaks)
4983 
4984  if (grads_nc_tweaks) then
4985  call check( nf90_def_dim(ncid, 'x', 1, ncd), thisroutine )
4986  call check( nf90_def_dim(ncid, 'y', 1, ncd), thisroutine )
4987  end if
4988 
4989  call check( nf90_def_dim(ncid, 'n', n_core, ncd), thisroutine )
4990  call check( nf90_def_dim(ncid, 't', nf90_unlimited, ncd), thisroutine )
4991 
4992 ! ------ Definition of the variables
4993 
4994  if (grads_nc_tweaks) then
4995 
4996 ! ---- x
4997 
4998  call check( nf90_inq_dimid(ncid, 'x', nc1d), thisroutine )
4999  call check( nf90_def_var(ncid, 'x', nf90_float, nc1d, ncv), &
5000  thisroutine )
5001  buffer = 'm'
5002  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
5003  thisroutine )
5004  buffer = 'x'
5005  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
5006  thisroutine )
5007  buffer = 'Dummy x-coordinate for one point'
5008  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
5009  thisroutine )
5010  call check( nf90_put_att(ncid, ncv, 'axis', 'x'), thisroutine )
5011 
5012 ! ---- y
5013 
5014  call check( nf90_inq_dimid(ncid, 'y', nc1d), thisroutine )
5015  call check( nf90_def_var(ncid, 'y', nf90_float, nc1d, ncv), &
5016  thisroutine )
5017  buffer = 'm'
5018  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
5019  thisroutine )
5020  buffer = 'y'
5021  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
5022  thisroutine )
5023  buffer = 'Dummy y-coordinate for one point'
5024  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
5025  thisroutine )
5026  call check( nf90_put_att(ncid, ncv, 'axis', 'y'), thisroutine )
5027 
5028  end if
5029 
5030 ! ---- Time
5031 
5032  call check( nf90_inq_dimid(ncid, 't', nc1d), thisroutine )
5033  call check( nf90_def_var(ncid, 't', nf90_float, nc1d, ncv), &
5034  thisroutine )
5035  buffer = 'a'
5036  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
5037  thisroutine )
5038  buffer = 'time'
5039  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
5040  thisroutine )
5041  buffer = 'Time'
5042  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
5043  thisroutine )
5044  call check( nf90_put_att(ncid, ncv, 'axis', 't'), thisroutine )
5045 
5046  if (grads_nc_tweaks) then
5047 
5048 ! ---- Time offset
5049 
5050  call check( nf90_inq_dimid(ncid, 't', nc1d), thisroutine )
5051  call check( nf90_def_var(ncid, 't_add_offset', &
5052  nf90_float, nc1d, ncv), &
5053  thisroutine )
5054  buffer = 'a'
5055  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
5056  thisroutine )
5057  buffer = 'time_add_offset'
5058  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
5059  thisroutine )
5060  buffer = 'Time offset'
5061  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
5062  thisroutine )
5063 
5064  time_add_offset_val = min(time_val, 0.0_dp)
5065 
5066  end if
5067 
5068 ! ---- Ice core number
5069 
5070  call check( nf90_inq_dimid(ncid, 'n', nc1d), thisroutine )
5071  call check( nf90_def_var(ncid, 'n', nf90_float, nc1d, ncv), &
5072  thisroutine )
5073  buffer = '1'
5074  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
5075  thisroutine )
5076  buffer = 'ice_core_number'
5077  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
5078  thisroutine )
5079  buffer = 'Ice core number'
5080  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
5081  thisroutine )
5082  buffer = trim(ch_core(1))
5083  do n=2, n_core; buffer = trim(buffer)//', '//trim(ch_core(n)); end do
5084  call check( nf90_put_att(ncid, ncv, 'ice_core_names', trim(buffer)), &
5085  thisroutine )
5086  call check( nf90_put_att(ncid, ncv, 'axis', 'n'), thisroutine )
5087 
5088  if (forcing_flag == 1) then
5089 
5090 ! ---- delta_ts
5091 
5092  call check( nf90_inq_dimid(ncid, 't', nc1d), thisroutine )
5093  call check( nf90_def_var(ncid, 'delta_ts', nf90_float, nc1d, ncv), &
5094  thisroutine )
5095  buffer = 'degC'
5096  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
5097  thisroutine )
5098  buffer = 'surface_temperature_anomaly'
5099  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
5100  thisroutine )
5101  buffer = 'Surface temperature anomaly'
5102  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
5103  thisroutine )
5104 
5105  else if ((forcing_flag == 2).or.(forcing_flag == 3)) then
5106 
5107 ! ---- glac_index
5108 
5109  call check( nf90_inq_dimid(ncid, 't', nc1d), thisroutine )
5110  call check( nf90_def_var(ncid, 'glac_index', nf90_float, nc1d, ncv), &
5111  thisroutine )
5112  buffer = '1'
5113  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
5114  thisroutine )
5115  buffer = 'glacial_index'
5116  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
5117  thisroutine )
5118  buffer = 'Glacial index'
5119  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
5120  thisroutine )
5121 
5122  end if
5123 
5124 ! ---- z_sl
5125 
5126  call check( nf90_inq_dimid(ncid, 't', nc1d), thisroutine )
5127  call check( nf90_def_var(ncid, 'z_sl', nf90_float, nc1d, ncv), &
5128  thisroutine )
5129  buffer = 'm'
5130  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
5131  thisroutine )
5132  buffer = 'global_average_sea_level_change'
5133  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
5134  thisroutine )
5135  buffer = 'Sea level'
5136  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
5137  thisroutine )
5138 
5139 ! ---- H_core
5140 
5141  call check( nf90_inq_dimid(ncid, 'n', nc2d(1)), thisroutine )
5142  call check( nf90_inq_dimid(ncid, 't', nc2d(2)), thisroutine )
5143  call check( nf90_def_var(ncid, 'H_core', nf90_float, nc2d, ncv), &
5144  thisroutine )
5145  buffer = 'm'
5146  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
5147  thisroutine )
5148  buffer = 'land_ice_thickness'
5149  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
5150  thisroutine )
5151  buffer = 'Ice thickness'
5152  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
5153  thisroutine )
5154 
5155 ! ---- vh_b_core
5156 
5157  call check( nf90_inq_dimid(ncid, 'n', nc2d(1)), thisroutine )
5158  call check( nf90_inq_dimid(ncid, 't', nc2d(2)), thisroutine )
5159  call check( nf90_def_var(ncid, 'vh_b_core', nf90_float, nc2d, ncv), &
5160  thisroutine )
5161  buffer = 'm a-1'
5162  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
5163  thisroutine )
5164  buffer = 'land_ice_basal_horizontal_velocity'
5165  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
5166  thisroutine )
5167  buffer = 'Horizontal velocity at the ice base'
5168  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
5169  thisroutine )
5170 
5171 ! ---- vh_s_core
5172 
5173  call check( nf90_inq_dimid(ncid, 'n', nc2d(1)), thisroutine )
5174  call check( nf90_inq_dimid(ncid, 't', nc2d(2)), thisroutine )
5175  call check( nf90_def_var(ncid, 'vh_s_core', nf90_float, nc2d, ncv), &
5176  thisroutine )
5177  buffer = 'm a-1'
5178  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
5179  thisroutine )
5180  buffer = 'land_ice_surface_horizontal_velocity'
5181  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
5182  thisroutine )
5183  buffer = 'Horizontal velocity at the ice surface'
5184  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
5185  thisroutine )
5186 
5187 ! ---- temp_b_core
5188 
5189  call check( nf90_inq_dimid(ncid, 'n', nc2d(1)), thisroutine )
5190  call check( nf90_inq_dimid(ncid, 't', nc2d(2)), thisroutine )
5191  call check( nf90_def_var(ncid, 'temp_b_core', nf90_float, nc2d, ncv), &
5192  thisroutine )
5193  buffer = 'degC'
5194  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
5195  thisroutine )
5196  buffer = 'basal_temperature'
5197  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
5198  thisroutine )
5199  buffer = 'Temperature at the ice base'
5200  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
5201  thisroutine )
5202 
5203 ! ---- R_b_core
5204 
5205  call check( nf90_inq_dimid(ncid, 'n', nc2d(1)), thisroutine )
5206  call check( nf90_inq_dimid(ncid, 't', nc2d(2)), thisroutine )
5207  call check( nf90_def_var(ncid, 'R_b_core', nf90_float, nc2d, ncv), &
5208  thisroutine )
5209  buffer = 'W m-2'
5210  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
5211  thisroutine )
5212  buffer = 'basal_frictional_heating'
5213  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
5214  thisroutine )
5215  buffer = 'Basal frictional heating'
5216  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
5217  thisroutine )
5218 
5219 ! ---- bmb_core
5220 
5221  call check( nf90_inq_dimid(ncid, 'n', nc2d(1)), thisroutine )
5222  call check( nf90_inq_dimid(ncid, 't', nc2d(2)), thisroutine )
5223  call check( nf90_def_var(ncid, 'bmb_core', nf90_float, nc2d, ncv), &
5224  thisroutine )
5225  buffer = 'm ice equiv. a-1'
5226  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
5227  thisroutine )
5228  buffer = 'land_ice_basal_mass_balance'
5229  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
5230  thisroutine )
5231  buffer = 'Basal mass balance'
5232  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
5233  thisroutine )
5234 
5235 ! ---- End of the definitions
5236 
5237  call check( nf90_enddef(ncid), thisroutine )
5238 
5239  end if
5240 
5241 ! ------ Writing of data on NetCDF file
5242 
5243  if (firstcall.and.grads_nc_tweaks) then
5244 
5245  nc1cor = (/ 1 /)
5246 
5247  call check( nf90_inq_varid(ncid, 'x', ncv), thisroutine )
5248  call check( nf90_put_var(ncid, ncv, 0.0_sp, start=nc1cor), thisroutine )
5249 
5250  call check( nf90_inq_varid(ncid, 'y', ncv), thisroutine )
5251  call check( nf90_put_var(ncid, ncv, 0.0_sp, start=nc1cor), thisroutine )
5252 
5253  call check( nf90_sync(ncid), thisroutine )
5254 
5255  end if
5256 
5257  if (firstcall) then
5258 
5259  nc1cor = (/ 1 /)
5260 
5261  call check( nf90_inq_varid(ncid, 'n', ncv), thisroutine )
5262  call check( nf90_put_var(ncid, ncv, r_n_core, &
5263  start=nc1cor), thisroutine )
5264 
5265  end if
5266 
5267  nc1cor(1) = counter
5268  ! nc1cnt(1) = 1 ! (not needed, and causes troubles for whatever reason)
5269 
5270  nc2cor(1) = 1
5271  nc2cor(2) = counter
5272  nc2cnt(1) = n_core
5273  nc2cnt(2) = 1
5274 
5275  call check( nf90_inq_varid(ncid, 't', ncv), thisroutine )
5276 
5277  if (.not.grads_nc_tweaks) then
5278  call check( nf90_put_var(ncid, ncv, real(time_val,sp), &
5279  start=nc1cor), thisroutine )
5280  else
5281  call check( nf90_put_var(ncid, ncv, &
5282  real(time_val-time_add_offset_val,sp), &
5283  start=nc1cor), thisroutine )
5284  ! this makes sure that all times are >=0
5285  ! (GrADS doesn't like negative numbers)
5286  call check( nf90_inq_varid(ncid, 't_add_offset', ncv), thisroutine )
5287  call check( nf90_put_var(ncid, ncv, &
5288  real(time_add_offset_val,sp), &
5289  start=nc1cor), thisroutine )
5290  end if
5291 
5292  if (forcing_flag == 1) then
5293 
5294  call check( nf90_inq_varid(ncid, 'delta_ts', ncv), thisroutine )
5295  call check( nf90_put_var(ncid, ncv, real(delta_ts,sp), &
5296  start=nc1cor), thisroutine )
5297 
5298  else if (forcing_flag == 2) then
5299 
5300  call check( nf90_inq_varid(ncid, 'glac_index', ncv), thisroutine )
5301  call check( nf90_put_var(ncid, ncv, real(glac_index,sp), &
5302  start=nc1cor), thisroutine )
5303 
5304  else if (forcing_flag == 3) then
5305 
5306  call check( nf90_inq_varid(ncid, 'glac_index', ncv), thisroutine )
5307  call check( nf90_put_var(ncid, ncv, real(no_value_neg_2,sp), &
5308  start=nc1cor), thisroutine )
5309 
5310  end if
5311 
5312  call check( nf90_inq_varid(ncid, 'z_sl', ncv), thisroutine )
5313  call check( nf90_put_var(ncid, ncv, real(z_sl,sp), &
5314  start=nc1cor), thisroutine )
5315 
5316  call check( nf90_inq_varid(ncid, 'H_core', ncv), thisroutine )
5317  call check( nf90_put_var(ncid, ncv, real(H_core,sp), &
5318  start=nc2cor, count=nc2cnt), thisroutine )
5319 
5320  call check( nf90_inq_varid(ncid, 'vh_b_core', ncv), thisroutine )
5321  call check( nf90_put_var(ncid, ncv, real(vh_b_core,sp), &
5322  start=nc2cor, count=nc2cnt), thisroutine )
5323 
5324  call check( nf90_inq_varid(ncid, 'vh_s_core', ncv), thisroutine )
5325  call check( nf90_put_var(ncid, ncv, real(vh_s_core,sp), &
5326  start=nc2cor, count=nc2cnt), thisroutine )
5327 
5328  call check( nf90_inq_varid(ncid, 'temp_b_core', ncv), thisroutine )
5329  call check( nf90_put_var(ncid, ncv, real(temp_b_core,sp), &
5330  start=nc2cor, count=nc2cnt), thisroutine )
5331 
5332  call check( nf90_inq_varid(ncid, 'R_b_core', ncv), thisroutine )
5333  call check( nf90_put_var(ncid, ncv, real(R_b_core,sp), &
5334  start=nc2cor, count=nc2cnt), thisroutine )
5335 
5336  call check( nf90_inq_varid(ncid, 'bmb_core', ncv), thisroutine )
5337  call check( nf90_put_var(ncid, ncv, real(bmb_core,sp), &
5338  start=nc2cor, count=nc2cnt), thisroutine )
5339 
5340 ! ------ Syncing NetCDF file (every 100th time)
5341 
5342  n_sync = 100
5343 
5344  if ( mod((counter-1), n_sync) == 0 ) &
5345  call check( nf90_sync(ncid), thisroutine )
5346 
5347 #endif
5348 
5349  deallocate(r_n_core, h_core, &
5350  vx_b_core, vy_b_core, vh_b_core, &
5351  vx_s_core, vy_s_core, vh_s_core, &
5352  temp_b_core, &
5353  rx_b_core, ry_b_core, r_b_core, &
5354  bmb_core)
5355 
5356 !!! else ! (n_core == 0 -> do nothing)
5357 !!!
5358 !!! continue
5359 
5360 end if
5361 
5362 if (firstcall) firstcall = .false.
5363 
5364 end subroutine output4
5365 
5366 #if (defined(ASF))
5367 
5368 !-------------------------------------------------------------------------------
5369 !> Writing of time-series data for all defined surface points on file
5370 !! in ASCII format. Modification of Tolly's output7 by Thorben Dunse.
5371 !<------------------------------------------------------------------------------
5372 subroutine output5(time, dxi, deta, delta_ts, glac_index, z_sl)
5373 
5374 implicit none
5375 
5376 real(dp), intent(in) :: time, dxi, deta, delta_ts, glac_index, z_sl
5377 
5378 integer(i4b) :: n, k
5379 real(dp) :: time_val
5380 real(dp), dimension(0:JMAX,0:IMAX) :: field
5381 real(dp), dimension(:), allocatable :: zl_surf, zs_surf, &
5382  accum_surf, as_perp_surf, &
5383  snowfall_surf, rainfall_surf, runoff_surf, &
5384  vx_surf, vy_surf, vz_surf, &
5385  vx_base, vy_base, vz_base, &
5386  temp_base_pmp
5387 
5388 allocate(zl_surf(n_surf), zs_surf(n_surf), &
5389  accum_surf(n_surf), &
5390  as_perp_surf(n_surf), snowfall_surf(n_surf), &
5391  rainfall_surf(n_surf), runoff_surf(n_surf), &
5392  vx_surf(n_surf), vy_surf(n_surf), vz_surf(n_surf), &
5393  vx_base(n_surf), vy_base(n_surf), vz_base(n_surf), &
5394  temp_base_pmp(n_surf))
5395 
5396 !-------- Determination of ice-core data --------
5397 
5398 do n=1, n_surf
5399 
5400 ! ------ Bedrock elevation
5401 
5402  field = zl
5403  call borehole(field, x_surf(n), y_surf(n), dxi, deta, &
5404  'grid', zl_surf(n))
5405 
5406 ! ------ Surface elevation
5407 
5408  field = zs
5409  call borehole(field, x_surf(n), y_surf(n), dxi, deta, &
5410  'grid', zs_surf(n))
5411 
5412 
5413 ! ------ Accumulation
5414 
5415  field = accum
5416  call borehole(field, x_surf(n), y_surf(n), dxi, deta, &
5417  'grid', accum_surf(n))
5418 
5419 ! ------ Surface mass balance
5420 
5421  field = as_perp
5422  call borehole(field, x_surf(n), y_surf(n), dxi, deta, &
5423  'grid', as_perp_surf(n))
5424 
5425 ! ------ Snowfall
5426 
5427  field = snowfall
5428  call borehole(field, x_surf(n), y_surf(n), dxi, deta, &
5429  'grid', snowfall_surf(n))
5430 
5431 ! ------ Rainfall
5432 
5433  field = rainfall
5434  call borehole(field, x_surf(n), y_surf(n), dxi, deta, &
5435  'grid', rainfall_surf(n))
5436 
5437 ! ------ Runoff
5438 
5439  field = runoff
5440  call borehole(field, x_surf(n), y_surf(n), dxi, deta, &
5441  'grid', runoff_surf(n))
5442 
5443 ! ------ Surface velocities
5444 
5445  field = vx_s_g
5446  call borehole(field, x_surf(n), y_surf(n), dxi, deta, &
5447  'grid', vx_surf(n))
5448 
5449  field = vy_s_g
5450  call borehole(field, x_surf(n), y_surf(n), dxi, deta, &
5451  'grid', vy_surf(n))
5452 
5453  field = vz_s
5454  call borehole(field, x_surf(n), y_surf(n), dxi, deta, &
5455  'grid', vz_surf(n))
5456 
5457 ! ------ Basal velocities
5458 
5459  field = vx_b_g
5460  call borehole(field, x_surf(n), y_surf(n), dxi, deta, &
5461  'grid', vx_base(n))
5462 
5463  field = vy_b_g
5464  call borehole(field, x_surf(n), y_surf(n), dxi, deta, &
5465  'grid', vy_base(n))
5466 
5467  field = vz_b
5468  call borehole(field, x_surf(n), y_surf(n), dxi, deta, &
5469  'grid', vz_base(n))
5470 
5471 ! ------ Basal temperature relative to pressure melting point
5472 
5473  field = temph_b
5474  call borehole(field, x_surf(n), y_surf(n), dxi, deta, &
5475  'grid', temp_base_pmp(n))
5476 
5477 
5478 end do
5479 
5480 ! ------ Conversion
5481 
5482 #if (!defined(OUT_TIMES) || OUT_TIMES==1)
5483 time_val = time *sec_to_year ! s -> a
5484 #elif (OUT_TIMES==2)
5485 time_val = (time+year_zero) *sec_to_year ! s -> a
5486 #else
5487 stop ' >>> output5: OUT_TIMES must be either 1 or 2!'
5488 #endif
5489 
5490 do n=1, n_surf
5491  accum_surf(n) = accum_surf(n) *year_sec ! m/s -> m/a
5492  as_perp_surf(n) = as_perp_surf(n) *year_sec ! m/s -> m/a
5493  snowfall_surf(n) = snowfall_surf(n) *year_sec ! m/s -> m/a
5494  rainfall_surf(n) = rainfall_surf(n) *year_sec ! m/s -> m/a
5495  runoff_surf(n) = runoff_surf(n) *year_sec ! m/s -> m/a
5496  vx_surf(n) = vx_surf(n) *year_sec ! m/s -> m/a
5497  vy_surf(n) = vy_surf(n) *year_sec ! m/s -> m/a
5498  vz_surf(n) = vz_surf(n) *year_sec ! m/s -> m/a
5499  vx_base(n) = vx_base(n) *year_sec ! m/s -> m/a
5500  vy_base(n) = vy_base(n) *year_sec ! m/s -> m/a
5501  vz_base(n) = vz_base(n) *year_sec ! m/s -> m/a
5502 end do
5503 
5504 !-------- Writing of data on file --------
5505 
5506 if (forcing_flag == 1) then
5507  write(41,'(1pe13.6,2(1pe13.4))') time_val, delta_ts, z_sl
5508  write(42,'(1pe13.6,2(1pe13.4))') time_val, delta_ts, z_sl
5509  write(43,'(1pe13.6,2(1pe13.4))') time_val, delta_ts, z_sl
5510  write(44,'(1pe13.6,2(1pe13.4))') time_val, delta_ts, z_sl
5511  write(45,'(1pe13.6,2(1pe13.4))') time_val, delta_ts, z_sl
5512  write(46,'(1pe13.6,2(1pe13.4))') time_val, delta_ts, z_sl
5513  write(47,'(1pe13.6,2(1pe13.4))') time_val, delta_ts, z_sl
5514  write(48,'(1pe13.6,2(1pe13.4))') time_val, delta_ts, z_sl
5515  write(49,'(1pe13.6,2(1pe13.4))') time_val, delta_ts, z_sl
5516  write(50,'(1pe13.6,2(1pe13.4))') time_val, delta_ts, z_sl
5517  write(51,'(1pe13.6,2(1pe13.4))') time_val, delta_ts, z_sl
5518  write(52,'(1pe13.6,2(1pe13.4))') time_val, delta_ts, z_sl
5519  write(53,'(1pe13.6,2(1pe13.4))') time_val, delta_ts, z_sl
5520  write(54,'(1pe13.6,2(1pe13.4))') time_val, delta_ts, z_sl
5521 else if (forcing_flag == 2) then
5522  write(41,'(1pe13.6,2(1pe13.4))') time_val, glac_index, z_sl
5523  write(42,'(1pe13.6,2(1pe13.4))') time_val, glac_index, z_sl
5524  write(43,'(1pe13.6,2(1pe13.4))') time_val, glac_index, z_sl
5525  write(44,'(1pe13.6,2(1pe13.4))') time_val, glac_index, z_sl
5526  write(45,'(1pe13.6,2(1pe13.4))') time_val, glac_index, z_sl
5527  write(46,'(1pe13.6,2(1pe13.4))') time_val, glac_index, z_sl
5528  write(47,'(1pe13.6,2(1pe13.4))') time_val, glac_index, z_sl
5529  write(48,'(1pe13.6,2(1pe13.4))') time_val, glac_index, z_sl
5530  write(49,'(1pe13.6,2(1pe13.4))') time_val, glac_index, z_sl
5531  write(50,'(1pe13.6,2(1pe13.4))') time_val, glac_index, z_sl
5532  write(51,'(1pe13.6,2(1pe13.4))') time_val, glac_index, z_sl
5533  write(52,'(1pe13.6,2(1pe13.4))') time_val, glac_index, z_sl
5534  write(53,'(1pe13.6,2(1pe13.4))') time_val, glac_index, z_sl
5535  write(54,'(1pe13.6,2(1pe13.4))') time_val, glac_index, z_sl
5536 end if
5537 
5538 do n=1, n_surf-1
5539  write(41,'(1pe13.4)',advance='no') zl_surf(n)
5540  write(42,'(1pe13.4)',advance='no') zs_surf(n)
5541  write(43,'(1pe13.4)',advance='no') accum_surf(n)
5542  write(44,'(1pe13.4)',advance='no') as_perp_surf(n)
5543  write(45,'(1pe13.4)',advance='no') snowfall_surf(n)
5544  write(46,'(1pe13.4)',advance='no') rainfall_surf(n)
5545  write(47,'(1pe13.4)',advance='no') runoff_surf(n)
5546  write(48,'(1pe13.4)',advance='no') vx_surf(n)
5547  write(49,'(1pe13.4)',advance='no') vy_surf(n)
5548  write(50,'(1pe13.4)',advance='no') vz_surf(n)
5549  write(51,'(1pe13.4)',advance='no') vx_base(n)
5550  write(52,'(1pe13.4)',advance='no') vy_base(n)
5551  write(53,'(1pe13.4)',advance='no') vz_base(n)
5552  write(54,'(1pe13.4)',advance='no') temp_base_pmp(n)
5553 end do
5554 
5555 n=n_surf
5556  write(41,'(1pe13.4)') zl_surf(n)
5557  write(42,'(1pe13.4)') zs_surf(n)
5558  write(43,'(1pe13.4)') accum_surf(n)
5559  write(44,'(1pe13.4)') as_perp_surf(n)
5560  write(45,'(1pe13.4)') snowfall_surf(n)
5561  write(46,'(1pe13.4)') rainfall_surf(n)
5562  write(47,'(1pe13.4)') runoff_surf(n)
5563  write(48,'(1pe13.4)') vx_surf(n)
5564  write(49,'(1pe13.4)') vy_surf(n)
5565  write(50,'(1pe13.4)') vz_surf(n)
5566  write(51,'(1pe13.4)') vx_base(n)
5567  write(52,'(1pe13.4)') vy_base(n)
5568  write(53,'(1pe13.4)') vz_base(n)
5569  write(54,'(1pe13.4)') temp_base_pmp(n)
5570 
5571 deallocate(zl_surf, zs_surf, accum_surf, as_perp_surf, &
5572  snowfall_surf, rainfall_surf, runoff_surf, &
5573  vx_surf, vy_surf, vz_surf, &
5574  vx_base, vy_base, vz_base,temp_base_pmp)
5575 
5576 end subroutine output5
5577 
5578 #endif /* (defined(ASF)) */
5579 
5580 !-------------------------------------------------------------------------------
5581 !> Computation of an arbitrary field quantity for a given borehole
5582 !! position x_pos, y_pos by weighed averaging of the corresponding
5583 !! gridded 2-d field.
5584 !<------------------------------------------------------------------------------
5585  subroutine borehole(field, x_pos, y_pos, dxi, deta, ch_grid, field_val)
5587  use sico_maths_m, only : bilinint
5588 
5589  implicit none
5590 
5591  real(dp), dimension(0:JMAX,0:IMAX), intent(in) :: field
5592  real(dp), intent(in) :: x_pos, y_pos, dxi, deta
5593  character (len=*), intent(in) :: ch_grid
5594 
5595  real(dp), intent(out) :: field_val
5596 
5597  integer(i4b) :: i1, i2, j1, j2
5598  real(dp) :: real_i, real_j
5599 
5600 !-------- Neighbour points --------
5601 
5602  real_i = (x_pos-xi(0)) /dxi
5603  real_j = (y_pos-eta(0))/deta
5604 
5605  if (ch_grid=='sg_x') real_i = real_i - 0.5_dp
5606  if (ch_grid=='sg_y') real_j = real_j - 0.5_dp
5607 
5608  if (real_i < 0.5_dp*real(imax,dp)) then
5609  i1 = floor(real_i)
5610  i2 = i1 + 1
5611  else
5612  i2 = ceiling(real_i)
5613  i1 = i2 - 1
5614  end if
5615 
5616  if (real_j < 0.5_dp*real(jmax,dp)) then
5617  j1 = floor(real_j)
5618  j2 = j1 + 1
5619  else
5620  j2 = ceiling(real_j)
5621  j1 = j2 - 1
5622  end if
5623 
5624  if (ch_grid=='grid') then
5625  ! field(j,i) defined on grid points
5626 
5627  if ((i1 < 0).or.(j1 < 0).or.(i2 > imax).or.(j2 > jmax)) &
5628  stop ' >>> borehole: Borehole position out of domain!'
5629 
5630  else if (ch_grid=='sg_x') then
5631  ! field(j,i) defined on staggered grid in x direction
5632 
5633  if ((i1 < 0).or.(j1 < 0).or.(i2 > imax-1).or.(j2 > jmax)) &
5634  stop ' >>> borehole: Borehole position out of domain!'
5635 
5636  else if (ch_grid=='sg_y') then
5637  ! field(j,i) defined on staggered grid in y direction
5638 
5639  if ((i1 < 0).or.(j1 < 0).or.(i2 > imax).or.(j2 > jmax-1)) &
5640  stop ' >>> borehole: Borehole position out of domain!'
5641 
5642  else
5643 
5644  stop ' >>> borehole: Parameter ch_grid has undefined value!'
5645 
5646  end if
5647 
5648 !-------- Weighing of the four adjacent grid values --------
5649 
5650  field_val = bilinint(real(i1,dp), real(i2,dp), real(j1,dp), real(j2,dp), &
5651  field(j1,i1), field(j2,i1), field(j1,i2), field(j2,i2), &
5652  real_i, real_j)
5653 
5654  end subroutine borehole
5655 
5656 !-------------------------------------------------------------------------------
5657 !> Set the value of the auxiliary variable grads_nc_tweaks.
5658 !<------------------------------------------------------------------------------
5659  subroutine set_grads_nc_tweaks(grads_nc_tweaks)
5661  implicit none
5662 
5663  logical, intent(out) :: grads_nc_tweaks
5664 
5665  character(len=16) :: ch_value
5666 
5667 #if (NETCDF>1)
5668 
5669  grads_nc_tweaks = .false.
5670 
5671 !-------- Try environment variable --------
5672 
5673  call get_environment_variable('SICO_GRADS_NC_TWEAKS', ch_value)
5674 
5675  if ( (trim(ch_value)=='true') &
5676  .or.(trim(ch_value)=='True').or.(trim(ch_value)=='TRUE') ) &
5677  grads_nc_tweaks = .true.
5678 
5679  if ( (trim(ch_value)=='yes') & ! obsolete, but still supported
5680  .or.(trim(ch_value)=='Yes').or.(trim(ch_value)=='YES') &
5681  .or.(trim(ch_value)=='y').or.(trim(ch_value)=='Y') ) &
5682  grads_nc_tweaks = .true.
5683 
5684 !-------- Try preprocessor switch --------
5685 
5686 #if (defined(GRADS_NC_TWEAKS))
5687 #if (GRADS_NC_TWEAKS==1)
5688  grads_nc_tweaks = .true.
5689 #else
5690  grads_nc_tweaks = .false.
5691 #endif
5692 #endif
5693 
5694 #else /* NetCDF not used */
5695  grads_nc_tweaks = .false. ! dummy value
5696 #endif
5697 
5698  end subroutine set_grads_nc_tweaks
5699 
5700 !-------------------------------------------------------------------------------
5701 
5702 end module output_m
5703 !
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
integer(i2b), dimension(0:jmax, 0:imax), public mask_mar
integer(i4b) ncid_ser
ncid_ser: ID of the NetCDF time-series output file
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), dimension(0:jmax, 0:imax), public cos_grad_tc
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) rho_w
RHO_W: Density of pure water.
real(dp), dimension(0:jmax, 0:imax) vx_m
vx_m(j,i): Mean (depth-averaged) velocity in x-direction, at (i+1/2,j)
real(dp), dimension(0:krmax, 0:jmax, 0:imax) temp_r
temp_r(kr,j,i): Temperature in the bedrock
real(dp), dimension(0:jmax, 0:imax) mb_source_apl
mb_source_apl(j,i): Applied mass balance source (SMB, BMB, calving)
real(dp) h_r
H_R: Thickness of the modelled lithosphere layer.
real(dp), dimension(0:jmax, 0:imax) vz_b
vz_b(j,i): Velocity in z-direction at the ice base, at (i,j)
subroutine, public output1(runname, time, delta_ts, glac_index, z_sl, flag_3d_output, ndat2d, ndat3d)
Writing of time-slice files in native binary or NetCDF format.
Definition: output_m.F90:59
real(dp), parameter eps
eps: Small number
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)
integer(i2b) forcing_flag
forcing_flag: Flag for the forcing type. 1: forcing by a spatially constant surface temperature anoma...
real(dp), dimension(0:jmax, 0:imax) calv_grounded
calv_grounded(j,i): Calving rate of grounded ice
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) enth_t
enth_t(kt,j,i): Enthalpy in the lower (kt) ice domain
Conversion from temperature (temp) and water content (omega) to enthalpy (enth) and vice versa...
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) snowfall
snowfall(j,i): Snowfall rate at the ice surface
Definition: sico_vars_m.F90:53
Several mathematical tools used by SICOPOLIS.
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) as_perp
as_perp(j,i): Accumulation-ablation function at the ice surface (SMB)
real(dp), dimension(0:jmax, 0:imax) area
area(j,i): Area of grid cell associated with grid point (i,j)
real(dp), dimension(0:jmax, 0:imax) dzs_deta_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
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
real(dp), dimension(0:jmax, 0:imax) runoff
runoff(j,i): Runoff rate at the ice surface
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) temp_t_m
temp_t_m(kt,j,i): Melting temperature in the lower (kt) ice domain
real(dp), dimension(0:jmax, 0:imax), public cst_dist
integer(i2b), dimension(0:jmax, 0:imax) maske_sedi
maske_sedi(j,i): Sediment mask. 1: hard rock, 7: soft sediment, 2: ocean.
Definition: sico_vars_m.F90:46
real(dp), dimension(:), allocatable y_surf
y_surf(n): Coordinate eta (= y) of the prescribed surface points
Definition: sico_vars_m.F90:51
real(dp), parameter pi_180_inv
pi_180_inv: 180 divided by pi (-> rad to deg)
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.
subroutine, public output4(time, dxi, deta, delta_ts, glac_index, z_sl)
Writing of time-series data of the deep ice cores on file in ASCII format (and optionally in NetCDF f...
Definition: output_m.F90:4751
real(dp), dimension(0:jmax, 0:imax) vy_m
vy_m(j,i): Mean (depth-averaged) velocity in y-direction, at (i,j+1/2)
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.
real(dp), parameter sec_to_year
Definition: output_m.F90:44
real(dp), dimension(0:jmax, 0:imax) runoff_bot
runoff_bot(j,i): Runoff rate (accounting at ice bottom)
Declarations of global variables for SICOPOLIS (for the ANT domain).
Definition: sico_vars_m.F90:35
subroutine, public borehole(field, x_pos, y_pos, dxi, deta, ch_grid, field_val)
Computation of an arbitrary field quantity for a given borehole position x_pos, y_pos by weighed aver...
Definition: output_m.F90:5586
real(dp), parameter no_value_neg_2
no_value_neg_2: Negative no-value parameter