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