SICOPOLIS V3.1
 All Classes Files Functions Variables Macros
output_nc.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Subroutine : o u t p u t _ n c
4 !
5 !> @file
6 !!
7 !! Writing of time-slice files in NetCDF format.
8 !!
9 !! @section Copyright
10 !!
11 !! Copyright 2009-2013 Ralf Greve, Thomas Goelles
12 !!
13 !! @section License
14 !!
15 !! This file is part of SICOPOLIS.
16 !!
17 !! SICOPOLIS is free software: you can redistribute it and/or modify
18 !! it under the terms of the GNU General Public License as published by
19 !! the Free Software Foundation, either version 3 of the License, or
20 !! (at your option) any later version.
21 !!
22 !! SICOPOLIS is distributed in the hope that it will be useful,
23 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
24 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25 !! GNU General Public License for more details.
26 !!
27 !! You should have received a copy of the GNU General Public License
28 !! along with SICOPOLIS. If not, see <http://www.gnu.org/licenses/>.
29 !<
30 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
31 
32 !-------------------------------------------------------------------------------
33 !> Writing of time-slice files in NetCDF format.
34 !<------------------------------------------------------------------------------
35 subroutine output_nc(runname, time, delta_ts, glac_index, z_sl, &
36  flag_3d_output, ndat2d, ndat3d)
37 
38 use sico_types
40 use netcdf
41 use nc_check
42 
43 implicit none
44 
45 real(dp), intent (in) :: time, delta_ts, glac_index, z_sl
46 character (len=100), intent (in) :: runname
47 logical, intent (in) :: flag_3d_output
48 
49 integer(i4b), intent (inout) :: ndat2d, ndat3d
50 
51 integer(i4b) :: i, j, kc, kt, kr
52 integer(i4b) :: ndat, ndat_help, ndat_1000s, ndat_100s, ndat_10s, ndat_1s
53 real(dp) :: v_tot, a_grounded, a_floating
54 real(sp) :: lon0, lat0
55 character (len=256) :: filename
56 character :: ch_1000s, ch_100s, ch_10s, ch_1s
57 
58 integer(i2b), dimension(0:IMAX,0:JMAX) :: maske_conv, n_cts_conv
59 real(sp) :: time_conv, delta_ts_conv, glac_index_conv, z_sl_conv, &
60  v_tot_conv, a_grounded_conv, a_floating_conv, &
61  h_r_conv, &
62  xi_conv(0:imax), eta_conv(0:jmax), &
63  sigma_level_c_conv(0:kcmax), sigma_level_t_conv(0:ktmax), &
64  sigma_level_r_conv(0:krmax)
65 real(sp), dimension(0:IMAX,0:JMAX) :: lambda_conv, phi_conv, &
66  lon_conv, lat_conv, &
67  temp_s_conv, as_perp_conv, &
68  zs_conv, zm_conv, zb_conv, zl_conv, h_c_conv, h_t_conv, h_conv, &
69  q_bm_conv, q_tld_conv, &
70  am_perp_conv, &
71  qx_conv, qy_conv, &
72  dzs_dtau_conv, dzm_dtau_conv, dzb_dtau_conv, dzl_dtau_conv, &
73  dh_c_dtau_conv, dh_t_dtau_conv, dh_dtau_conv, &
74  vx_b_g_conv, vy_b_g_conv, vz_b_conv, vh_b_conv, &
75  vx_s_g_conv, vy_s_g_conv, vz_s_conv, vh_s_conv, &
76  temp_b_conv, temph_b_conv, &
77  p_b_w_conv, h_w_conv, q_gl_g_conv
78 real(sp), dimension(0:IMAX,0:JMAX,0:KCMAX) :: vx_c_conv, vy_c_conv, vz_c_conv, &
79  temp_c_conv, age_c_conv
80 real(sp), dimension(0:IMAX,0:JMAX,0:KTMAX) :: vx_t_conv, vy_t_conv, vz_t_conv, &
81  omega_t_conv, age_t_conv
82 real(sp), dimension(0:IMAX,0:JMAX,0:KRMAX) :: temp_r_conv
83 
84 integer(i4b) :: ncid, ncv
85 ! ncid: ID of the output file
86 ! ncv: Variable ID
87 integer(i4b) :: ncd, nc1d, nc2d(2), nc3d(3)
88 ! ncd: Dimension ID
89 ! nc1d: Dimension of a 1-d array
90 ! nc2d: Vector with the dimensions of a 2-d array
91 ! nc3d: Vector with the dimensions of a 3-d array
92 integer(i4b) :: nc3flag(3), nc4flag(4)
93 ! nc3flag: Vector with the 3 possible values of a flag variable
94 ! nc4flag: Vector with the 4 possible values of a flag variable
95 integer(i4b) :: nc1cor(1), nc2cor(2), nc3cor(3)
96 ! nc1cor(1): Corner of a 1-d array
97 ! nc2cor(2): Corner of a 2-d array
98 ! nc3cor(3): Corner of a 3-d array
99 integer(i4b) :: nc1cnt(1), nc2cnt(2), nc3cnt(3)
100 ! nc1cnt(1): Count of a 1-d array
101 ! nc2cnt(2): Count of a 2-d array
102 ! nc3cnt(3): Count of a 3-d array
103 character (len= 16) :: ch_date, ch_time, ch_zone
104 character (len=256) :: buffer
105 
106 character(len=64), parameter :: thisroutine = 'output_nc'
107 
108 !-------- Create consecutively numbered file names --------
109 
110 if (flag_3d_output) then
111  ndat = ndat3d
112 else
113  ndat = ndat2d
114 end if
115 
116 if (ndat > 9999) stop ' output_nc: Too many time-slice files!'
117 
118 ndat_help = ndat
119 ndat_1000s = ndat_help/1000
120 ndat_help = ndat_help-ndat_1000s*1000
121 ndat_100s = ndat_help/100
122 ndat_help = ndat_help-ndat_100s*100
123 ndat_10s = ndat_help/10
124 ndat_help = ndat_help-ndat_10s*10
125 ndat_1s = ndat_help
126 
127 ch_1000s = char(ndat_1000s+ichar('0'))
128 ch_100s = char(ndat_100s +ichar('0'))
129 ch_10s = char(ndat_10s +ichar('0'))
130 ch_1s = char(ndat_1s +ichar('0'))
131 
132 if (flag_3d_output) then
133  filename = trim(runname)//ch_1000s//ch_100s//ch_10s//ch_1s//'.nc'
134 else
135  filename = trim(runname)//'_2d_'//ch_1000s//ch_100s//ch_10s//ch_1s//'.nc'
136 end if
137 
138 !-------- NetCDF initialization --------
139 
140 ! ------ Open NetCDF file
141 
142 buffer = outpath//'/'//trim(filename)
143 call check( nf90_create(trim(buffer), nf90_noclobber, ncid), thisroutine )
144 
145 ! ------ Global attributes
146 
147 buffer = 'Time-slice output no. '//ch_1000s//ch_100s//ch_10s//ch_1s//' '// &
148  'of simulation '//trim(runname)
149 call check( nf90_put_att(ncid, nf90_global, 'title', trim(buffer)), &
150  thisroutine )
151 
152 buffer = 'Institute of Low Temperature Science, Hokkaido University, '// &
153  'Sapporo, Japan'
154 call check( nf90_put_att(ncid, nf90_global, 'institution', trim(buffer)), &
155  thisroutine )
156 
157 buffer = 'SICOPOLIS Version '//version
158 call check( nf90_put_att(ncid, nf90_global, 'source', trim(buffer)), &
159  thisroutine )
160 
161 call date_and_time(ch_date, ch_time, ch_zone)
162 buffer = ch_date(1:4)//'-'//ch_date(5:6)//'-'//ch_date(7:8)//' '// &
163  ch_time(1:2)//':'//ch_time(3:4)//':'//ch_time(5:6)//' '// &
164  ch_zone(1:3)//':'//ch_zone(4:5)//' - Data produced'
165 call check( nf90_put_att(ncid, nf90_global, 'history', trim(buffer)), &
166  thisroutine )
167 
168 buffer = 'http://sicopolis.greveweb.net/'
169 call check( nf90_put_att(ncid, nf90_global, 'references', trim(buffer)), &
170  thisroutine )
171 
172 ! ------ Definition of the dimensions
173 
174 call check( nf90_def_dim(ncid, 'IMAX', imax+1, ncd), thisroutine )
175 call check( nf90_def_dim(ncid, 'JMAX', jmax+1, ncd), thisroutine )
176 call check( nf90_def_dim(ncid, 'KCMAX', kcmax+1, ncd), thisroutine )
177 call check( nf90_def_dim(ncid, 'KTMAX', ktmax+1, ncd), thisroutine )
178 call check( nf90_def_dim(ncid, 'KRMAX', krmax+1, ncd), thisroutine )
179 
180 ! ------ Definition of the variables
181 
182 call check( nf90_def_var(ncid, 'crs', nf90_short, ncv), thisroutine )
183 #if ( (GRID == 0) || (GRID == 1) )
184 buffer = 'polar_stereographic'
185 call check( nf90_put_att(ncid, ncv, 'grid_mapping_name', trim(buffer)), &
186  thisroutine )
187 #elif (GRID == 2)
188 buffer = 'latitude_longitude'
189 call check( nf90_put_att(ncid, ncv, 'grid_mapping_name', trim(buffer)), &
190  thisroutine )
191 #endif
192 #if ( defined(ANT) \
193  || defined(asf) \
194  || defined(grl) \
195  || defined(scand) \
196  || defined(nhem) \
197  || defined(tibet) \
198  || defined(emtp2sge) \
199  || defined(heino) ) /* terrestrial ice sheet */
200 buffer = 'WGS84'
201 call check( nf90_put_att(ncid, ncv, 'ellipsoid', trim(buffer)), &
202  thisroutine )
203 #elif ( defined(NMARS) \
204  || defined(smars) ) /* martian ice sheet */
205 buffer = 'Mars_ellipsoid'
206 call check( nf90_put_att(ncid, ncv, 'ellipsoid', trim(buffer)), &
207  thisroutine )
208 #else
209 stop ' output_nc: No valid domain (ANT, GRL etc.) specified!'
210 #endif
211 #if ( (GRID == 0) || (GRID == 1) )
212 call check( nf90_put_att(ncid, ncv, 'false_easting', 0.0), &
213  thisroutine )
214 call check( nf90_put_att(ncid, ncv, 'false_northing', 0.0), &
215  thisroutine )
216 lon0 = lambda0 *pi_180_inv
217 lon0 = modulo(lon0+180.0_sp, 360.0_sp)-180.0_sp
218 lon0 = nint(lon0*1.0e+04_sp)*1.0e-04_sp
219 lat0 = phi0 *pi_180_inv
220 if (lat0 > 90.0_sp) lat0 = 90.0_sp
221 if (lat0 < -90.0_sp) lat0 = -90.0_sp
222 lat0 = nint(lat0*1.0e+04_sp)*1.0e-04_sp
223  ! reference longitude and standard parallel in deg rounded to 4 digits
224 if (lat0 >= 0.0_sp) then
225  call check( nf90_put_att(ncid, ncv, &
226  'latitude_of_projection_origin', 90.0), &
227  thisroutine )
228 else
229  call check( nf90_put_att(ncid, ncv, &
230  'latitude_of_projection_origin', -90.0), &
231  thisroutine )
232 end if
233 call check( nf90_put_att(ncid, ncv, &
234  'straight_vertical_longitude_from_pole', lon0), &
235  thisroutine )
236 call check( nf90_put_att(ncid, ncv, &
237  'standard_parallel', lat0), &
238  thisroutine )
239 #endif
240 
241 call check( nf90_def_var(ncid, 'time', nf90_float, ncv), &
242  thisroutine )
243 buffer = 'a'
244 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
245  thisroutine )
246 buffer = 'time'
247 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
248  thisroutine )
249 buffer = 'Time'
250 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
251  thisroutine )
252 
253 if (forcing_flag == 1) then
254 
255  call check( nf90_def_var(ncid, 'delta_ts', nf90_float, ncv), &
256  thisroutine )
257  buffer = 'degC'
258  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
259  thisroutine )
260  buffer = 'surface_temperature_anomaly'
261  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
262  thisroutine )
263  buffer = 'Surface temperature anomaly'
264  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
265  thisroutine )
266 
267 else if (forcing_flag == 2) then
268 
269  call check( nf90_def_var(ncid, 'glac_index', nf90_float, ncv), &
270  thisroutine )
271  buffer = '1'
272  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
273  thisroutine )
274  buffer = 'glacial_index'
275  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
276  thisroutine )
277  buffer = 'Glacial index'
278  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
279  thisroutine )
280 
281 else if (forcing_flag == 3) then
282 
283  call check( nf90_def_var(ncid, 'glac_index', nf90_float, ncv), &
284  thisroutine )
285  buffer = '1'
286  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
287  thisroutine )
288  buffer = 'glacial_index'
289  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
290  thisroutine )
291  buffer = 'Glacial index'
292  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
293  thisroutine )
294  buffer = 'This variable will be assigned a dummy value only!'
295  call check( nf90_put_att(ncid, ncv, 'comment', trim(buffer)), &
296  thisroutine )
297 
298 end if
299 
300 call check( nf90_def_var(ncid, 'z_sl', nf90_float, ncv), &
301  thisroutine )
302 buffer = 'm'
303 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
304  thisroutine )
305 buffer = 'global_average_sea_level_change'
306 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
307  thisroutine )
308 buffer = 'Sea level'
309 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
310  thisroutine )
311 
312 call check( nf90_inq_dimid(ncid, 'IMAX', nc1d), &
313  thisroutine )
314 call check( nf90_def_var(ncid, 'xi', nf90_float, nc1d, ncv), &
315  thisroutine )
316 buffer = 'm'
317 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
318  thisroutine )
319 buffer = 'projection_x_coordinate'
320 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
321  thisroutine )
322 buffer = 'x-coordinate of the grid point i'
323 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
324  thisroutine )
325 
326 call check( nf90_inq_dimid(ncid, 'JMAX', nc1d), &
327  thisroutine )
328 call check( nf90_def_var(ncid, 'eta', nf90_float, nc1d, ncv), &
329  thisroutine )
330 buffer = 'm'
331 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
332  thisroutine )
333 buffer = 'projection_y_coordinate'
334 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
335  thisroutine )
336 buffer = 'y-coordinate of the grid point j'
337 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
338  thisroutine )
339 
340 call check( nf90_inq_dimid(ncid, 'KCMAX', nc1d), &
341  thisroutine )
342 call check( nf90_def_var(ncid, 'sigma_level_c', nf90_float, nc1d, ncv), &
343  thisroutine )
344 buffer = 'up'
345 call check( nf90_put_att(ncid, ncv, 'positive', trim(buffer)), &
346  thisroutine )
347 buffer = 'land_ice_cold_layer_sigma_coordinate'
348 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
349  thisroutine )
350 buffer = 'sigma-coordinate of the grid point kc'
351 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
352  thisroutine )
353 
354 call check( nf90_inq_dimid(ncid, 'KTMAX', nc1d), &
355  thisroutine )
356 call check( nf90_def_var(ncid, 'sigma_level_t', nf90_float, nc1d, ncv), &
357  thisroutine )
358 buffer = 'up'
359 call check( nf90_put_att(ncid, ncv, 'positive', trim(buffer)), &
360  thisroutine )
361 buffer = 'land_ice_temperate_layer_sigma_coordinate'
362 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
363  thisroutine )
364 buffer = 'sigma-coordinate of the grid point kt'
365 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
366  thisroutine )
367 
368 call check( nf90_inq_dimid(ncid, 'KRMAX', nc1d), &
369  thisroutine )
370 call check( nf90_def_var(ncid, 'sigma_level_r', nf90_float, nc1d, ncv), &
371  thisroutine )
372 buffer = 'up'
373 call check( nf90_put_att(ncid, ncv, 'positive', trim(buffer)), &
374  thisroutine )
375 buffer = 'lithosphere_layer_sigma_coordinate'
376 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
377  thisroutine )
378 buffer = 'sigma-coordinate of the grid point kr'
379 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
380  thisroutine )
381 
382 call check( nf90_inq_dimid(ncid, 'IMAX', nc2d(1)), &
383  thisroutine )
384 call check( nf90_inq_dimid(ncid, 'JMAX', nc2d(2)), &
385  thisroutine )
386 call check( nf90_def_var(ncid, 'lon', nf90_float, nc2d, ncv), &
387  thisroutine )
388 buffer = 'degrees_E'
389 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
390  thisroutine )
391 buffer = 'longitude'
392 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
393  thisroutine )
394 buffer = 'Geographical longitude'
395 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
396  thisroutine )
397 
398 call check( nf90_inq_dimid(ncid, 'IMAX', nc2d(1)), &
399  thisroutine )
400 call check( nf90_inq_dimid(ncid, 'JMAX', nc2d(2)), &
401  thisroutine )
402 call check( nf90_def_var(ncid, 'lat', nf90_float, nc2d, ncv), &
403  thisroutine )
404 buffer = 'degrees_N'
405 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
406  thisroutine )
407 buffer = 'latitude'
408 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
409  thisroutine )
410 buffer = 'Geographical latitude'
411 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
412  thisroutine )
413 
414 call check( nf90_inq_dimid(ncid, 'IMAX', nc2d(1)), &
415  thisroutine )
416 call check( nf90_inq_dimid(ncid, 'JMAX', nc2d(2)), &
417  thisroutine )
418 call check( nf90_def_var(ncid, 'lambda', nf90_float, nc2d, ncv), &
419  thisroutine )
420 buffer = 'rad'
421 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
422  thisroutine )
423 buffer = 'longitude'
424 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
425  thisroutine )
426 buffer = 'Geographical longitude'
427 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
428  thisroutine )
429 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
430  thisroutine )
431 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
432  thisroutine )
433 
434 call check( nf90_inq_dimid(ncid, 'IMAX', nc2d(1)), &
435  thisroutine )
436 call check( nf90_inq_dimid(ncid, 'JMAX', nc2d(2)), &
437  thisroutine )
438 call check( nf90_def_var(ncid, 'phi', nf90_float, nc2d, ncv), &
439  thisroutine )
440 buffer = 'rad'
441 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
442  thisroutine )
443 buffer = 'latitude'
444 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
445  thisroutine )
446 buffer = 'Geographical latitude'
447 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
448  thisroutine )
449 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
450  thisroutine )
451 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
452  thisroutine )
453 
454 call check( nf90_inq_dimid(ncid, 'IMAX', nc2d(1)), &
455  thisroutine )
456 call check( nf90_inq_dimid(ncid, 'JMAX', nc2d(2)), &
457  thisroutine )
458 call check( nf90_def_var(ncid, 'temp_s', nf90_float, nc2d, ncv), &
459  thisroutine )
460 buffer = 'degC'
461 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
462  thisroutine )
463 buffer = 'surface_temperature'
464 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
465  thisroutine )
466 buffer = 'Temperature at the ice surface'
467 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
468  thisroutine )
469 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
470  thisroutine )
471 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
472  thisroutine )
473 
474 call check( nf90_inq_dimid(ncid, 'IMAX', nc2d(1)), &
475  thisroutine )
476 call check( nf90_inq_dimid(ncid, 'JMAX', nc2d(2)), &
477  thisroutine )
478 call check( nf90_def_var(ncid, 'as_perp', nf90_float, nc2d, ncv), &
479  thisroutine )
480 buffer = 'm a-1'
481 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
482  thisroutine )
483 buffer = 'land_ice_surface_mass_balance'
484 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
485  thisroutine )
486 buffer = 'Mass balance at the ice surface'
487 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
488  thisroutine )
489 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
490  thisroutine )
491 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
492  thisroutine )
493 
494 call check( nf90_inq_dimid(ncid, 'IMAX', nc2d(1)), &
495  thisroutine )
496 call check( nf90_inq_dimid(ncid, 'JMAX', nc2d(2)), &
497  thisroutine )
498 call check( nf90_def_var(ncid, 'maske', nf90_short, nc2d, ncv), &
499  thisroutine )
500 buffer = 'ice_land_sea_mask'
501 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
502  thisroutine )
503 buffer = 'Ice-land-sea mask'
504 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
505  thisroutine )
506 nc4flag(1) = 0
507 nc4flag(2) = 1
508 nc4flag(3) = 2
509 nc4flag(4) = 3
510 call check( nf90_put_att(ncid, ncv, 'flag_values', nc4flag), &
511  thisroutine )
512 buffer = 'glaciated_land '// &
513  'ice_free_land '// &
514  'sea '// &
515  'floating_ice'
516 call check( nf90_put_att(ncid, ncv, 'flag_meanings', trim(buffer)), &
517  thisroutine )
518 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
519  thisroutine )
520 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
521  thisroutine )
522 
523 call check( nf90_inq_dimid(ncid, 'IMAX', nc2d(1)), &
524  thisroutine )
525 call check( nf90_inq_dimid(ncid, 'JMAX', nc2d(2)), &
526  thisroutine )
527 call check( nf90_def_var(ncid, 'n_cts', nf90_short, nc2d, ncv), &
528  thisroutine )
529 buffer = 'polythermal_condition_mask'
530 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
531  thisroutine )
532 buffer = 'Mask for polythermal conditions'
533 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
534  thisroutine )
535 nc3flag(1) = -1
536 nc3flag(2) = 0
537 nc3flag(3) = 1
538 call check( nf90_put_att(ncid, ncv, 'flag_values', nc3flag), &
539  thisroutine )
540 buffer = 'cold_base '// &
541  'temperate_base_with_cold_ice_above '// &
542  'temperate_base_with_temperate_ice_above'
543 call check( nf90_put_att(ncid, ncv, 'flag_meanings', trim(buffer)), &
544  thisroutine )
545 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
546  thisroutine )
547 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
548  thisroutine )
549 
550 call check( nf90_inq_dimid(ncid, 'IMAX', nc2d(1)), &
551  thisroutine )
552 call check( nf90_inq_dimid(ncid, 'JMAX', nc2d(2)), &
553  thisroutine )
554 call check( nf90_def_var(ncid, 'zs', nf90_float, nc2d, ncv), &
555  thisroutine )
556 buffer = 'm'
557 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
558  thisroutine )
559 buffer = 'surface_altitude'
560 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
561  thisroutine )
562 buffer = 'Topography of the free surface'
563 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
564  thisroutine )
565 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
566  thisroutine )
567 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
568  thisroutine )
569 
570 call check( nf90_inq_dimid(ncid, 'IMAX', nc2d(1)), &
571  thisroutine )
572 call check( nf90_inq_dimid(ncid, 'JMAX', nc2d(2)), &
573  thisroutine )
574 call check( nf90_def_var(ncid, 'zm', nf90_float, nc2d, ncv), &
575  thisroutine )
576 buffer = 'm'
577 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
578  thisroutine )
579 buffer = 'cts_altitude'
580 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
581  thisroutine )
582 buffer = 'Topography of the CTS (cold-temperate transition surface)'
583 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
584  thisroutine )
585 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
586  thisroutine )
587 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
588  thisroutine )
589 
590 call check( nf90_inq_dimid(ncid, 'IMAX', nc2d(1)), &
591  thisroutine )
592 call check( nf90_inq_dimid(ncid, 'JMAX', nc2d(2)), &
593  thisroutine )
594 call check( nf90_def_var(ncid, 'zb', nf90_float, nc2d, ncv), &
595  thisroutine )
596 buffer = 'm'
597 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
598  thisroutine )
599 buffer = 'ice_base_altitude'
600 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
601  thisroutine )
602 buffer = 'Topography of the ice base'
603 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
604  thisroutine )
605 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
606  thisroutine )
607 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
608  thisroutine )
609 
610 call check( nf90_inq_dimid(ncid, 'IMAX', nc2d(1)), &
611  thisroutine )
612 call check( nf90_inq_dimid(ncid, 'JMAX', nc2d(2)), &
613  thisroutine )
614 call check( nf90_def_var(ncid, 'zl', nf90_float, nc2d, ncv), &
615  thisroutine )
616 buffer = 'm'
617 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
618  thisroutine )
619 buffer = 'bedrock_altitude'
620 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
621  thisroutine )
622 buffer = 'Topography of the lithosphere surface'
623 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
624  thisroutine )
625 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
626  thisroutine )
627 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
628  thisroutine )
629 
630 call check( nf90_inq_dimid(ncid, 'IMAX', nc2d(1)), &
631  thisroutine )
632 call check( nf90_inq_dimid(ncid, 'JMAX', nc2d(2)), &
633  thisroutine )
634 call check( nf90_def_var(ncid, 'H_c', nf90_float, nc2d, ncv), &
635  thisroutine )
636 buffer = 'm'
637 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
638  thisroutine )
639 buffer = 'land_ice_cold_layer_thickness'
640 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
641  thisroutine )
642 buffer = 'Thickness of the cold ice layer'
643 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
644  thisroutine )
645 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
646  thisroutine )
647 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
648  thisroutine )
649 
650 call check( nf90_inq_dimid(ncid, 'IMAX', nc2d(1)), &
651  thisroutine )
652 call check( nf90_inq_dimid(ncid, 'JMAX', nc2d(2)), &
653  thisroutine )
654 call check( nf90_def_var(ncid, 'H_t', nf90_float, nc2d, ncv), &
655  thisroutine )
656 buffer = 'm'
657 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
658  thisroutine )
659 buffer = 'land_ice_temperate_layer_thickness'
660 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
661  thisroutine )
662 buffer = 'Thickness of the temperate ice layer'
663 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
664  thisroutine )
665 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
666  thisroutine )
667 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
668  thisroutine )
669 
670 call check( nf90_inq_dimid(ncid, 'IMAX', nc2d(1)), &
671  thisroutine )
672 call check( nf90_inq_dimid(ncid, 'JMAX', nc2d(2)), &
673  thisroutine )
674 call check( nf90_def_var(ncid, 'H', nf90_float, nc2d, ncv), &
675  thisroutine )
676 buffer = 'm'
677 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
678  thisroutine )
679 buffer = 'land_ice_thickness'
680 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
681  thisroutine )
682 buffer = 'Ice thickness'
683 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
684  thisroutine )
685 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
686  thisroutine )
687 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
688  thisroutine )
689 
690 call check( nf90_def_var(ncid, 'H_R', nf90_float, ncv), &
691  thisroutine )
692 buffer = 'm'
693 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
694  thisroutine )
695 buffer = 'lithosphere_layer_thickness'
696 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
697  thisroutine )
698 buffer = 'Thickness of the lithosphere layer'
699 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
700  thisroutine )
701 
702 if (flag_3d_output) then
703 
704  call check( nf90_inq_dimid(ncid, 'IMAX', nc3d(1)), &
705  thisroutine )
706  call check( nf90_inq_dimid(ncid, 'JMAX', nc3d(2)), &
707  thisroutine )
708  call check( nf90_inq_dimid(ncid, 'KCMAX', nc3d(3)), &
709  thisroutine )
710  call check( nf90_def_var(ncid, 'vx_c', nf90_float, nc3d, ncv), &
711  thisroutine )
712  buffer = 'm a-1'
713  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
714  thisroutine )
715  buffer = 'land_ice_cold_layer_x_velocity'
716  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
717  thisroutine )
718  buffer = 'Horizontal velocity vx in the cold ice layer'
719  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
720  thisroutine )
721  buffer = 'Staggered grid variable, defined at (kc,j,i+1/2)'
722  call check( nf90_put_att(ncid, ncv, 'comment', trim(buffer)), &
723  thisroutine )
724  call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
725  thisroutine )
726  call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
727  thisroutine )
728 
729  call check( nf90_inq_dimid(ncid, 'IMAX', nc3d(1)), &
730  thisroutine )
731  call check( nf90_inq_dimid(ncid, 'JMAX', nc3d(2)), &
732  thisroutine )
733  call check( nf90_inq_dimid(ncid, 'KCMAX', nc3d(3)), &
734  thisroutine )
735  call check( nf90_def_var(ncid, 'vy_c', nf90_float, nc3d, ncv), &
736  thisroutine )
737  buffer = 'm a-1'
738  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
739  thisroutine )
740  buffer = 'land_ice_cold_layer_y_velocity'
741  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
742  thisroutine )
743  buffer = 'Horizontal velocity vy in the cold ice layer'
744  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
745  thisroutine )
746  buffer = 'Staggered grid variable, defined at (kc,j+1/2,i)'
747  call check( nf90_put_att(ncid, ncv, 'comment', trim(buffer)), &
748  thisroutine )
749  call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
750  thisroutine )
751  call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
752  thisroutine )
753 
754  call check( nf90_inq_dimid(ncid, 'IMAX', nc3d(1)), &
755  thisroutine )
756  call check( nf90_inq_dimid(ncid, 'JMAX', nc3d(2)), &
757  thisroutine )
758  call check( nf90_inq_dimid(ncid, 'KCMAX', nc3d(3)), &
759  thisroutine )
760  call check( nf90_def_var(ncid, 'vz_c', nf90_float, nc3d, ncv), &
761  thisroutine )
762  buffer = 'm a-1'
763  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
764  thisroutine )
765  buffer = 'land_ice_cold_layer_z_velocity'
766  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
767  thisroutine )
768  buffer = 'Vertical velocity vz in the cold ice layer'
769  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
770  thisroutine )
771  buffer = 'Staggered grid variable, defined at (kc+1/2,j,i)'
772  call check( nf90_put_att(ncid, ncv, 'comment', trim(buffer)), &
773  thisroutine )
774  call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
775  thisroutine )
776  call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
777  thisroutine )
778 
779  call check( nf90_inq_dimid(ncid, 'IMAX', nc3d(1)), &
780  thisroutine )
781  call check( nf90_inq_dimid(ncid, 'JMAX', nc3d(2)), &
782  thisroutine )
783  call check( nf90_inq_dimid(ncid, 'KTMAX', nc3d(3)), &
784  thisroutine )
785  call check( nf90_def_var(ncid, 'vx_t', nf90_float, nc3d, ncv), &
786  thisroutine )
787  buffer = 'm a-1'
788  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
789  thisroutine )
790  buffer = 'land_ice_temperate_layer_x_velocity'
791  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
792  thisroutine )
793  buffer = 'Horizontal velocity vx in the temperate ice layer'
794  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
795  thisroutine )
796  buffer = 'Staggered grid variable, defined at (kt,j,i+1/2)'
797  call check( nf90_put_att(ncid, ncv, 'comment', trim(buffer)), &
798  thisroutine )
799  call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
800  thisroutine )
801  call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
802  thisroutine )
803 
804  call check( nf90_inq_dimid(ncid, 'IMAX', nc3d(1)), &
805  thisroutine )
806  call check( nf90_inq_dimid(ncid, 'JMAX', nc3d(2)), &
807  thisroutine )
808  call check( nf90_inq_dimid(ncid, 'KTMAX', nc3d(3)), &
809  thisroutine )
810  call check( nf90_def_var(ncid, 'vy_t', nf90_float, nc3d, ncv), &
811  thisroutine )
812  buffer = 'm a-1'
813  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
814  thisroutine )
815  buffer = 'land_ice_temperate_layer_y_velocity'
816  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
817  thisroutine )
818  buffer = 'Horizontal velocity vy in the temperate ice layer'
819  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
820  thisroutine )
821  buffer = 'Staggered grid variable, defined at (kt,j+1/2,i)'
822  call check( nf90_put_att(ncid, ncv, 'comment', trim(buffer)), &
823  thisroutine )
824  call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
825  thisroutine )
826  call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
827  thisroutine )
828 
829  call check( nf90_inq_dimid(ncid, 'IMAX', nc3d(1)), &
830  thisroutine )
831  call check( nf90_inq_dimid(ncid, 'JMAX', nc3d(2)), &
832  thisroutine )
833  call check( nf90_inq_dimid(ncid, 'KTMAX', nc3d(3)), &
834  thisroutine )
835  call check( nf90_def_var(ncid, 'vz_t', nf90_float, nc3d, ncv), &
836  thisroutine )
837  buffer = 'm a-1'
838  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
839  thisroutine )
840  buffer = 'land_ice_temperate_layer_z_velocity'
841  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
842  thisroutine )
843  buffer = 'Vertical velocity vz in the temperate ice layer'
844  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
845  thisroutine )
846  buffer = 'Staggered grid variable, defined at (kt+1/2,j,i)'
847  call check( nf90_put_att(ncid, ncv, 'comment', trim(buffer)), &
848  thisroutine )
849  call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
850  thisroutine )
851  call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
852  thisroutine )
853 
854  call check( nf90_inq_dimid(ncid, 'IMAX', nc3d(1)), &
855  thisroutine )
856  call check( nf90_inq_dimid(ncid, 'JMAX', nc3d(2)), &
857  thisroutine )
858  call check( nf90_inq_dimid(ncid, 'KCMAX', nc3d(3)), &
859  thisroutine )
860  call check( nf90_def_var(ncid, 'temp_c', nf90_float, nc3d, ncv), &
861  thisroutine )
862  buffer = 'degC'
863  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
864  thisroutine )
865  buffer = 'land_ice_cold_layer_temperature'
866  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
867  thisroutine )
868  buffer = 'Temperature in the cold ice layer'
869  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
870  thisroutine )
871  call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
872  thisroutine )
873  call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
874  thisroutine )
875 
876  call check( nf90_inq_dimid(ncid, 'IMAX', nc3d(1)), &
877  thisroutine )
878  call check( nf90_inq_dimid(ncid, 'JMAX', nc3d(2)), &
879  thisroutine )
880  call check( nf90_inq_dimid(ncid, 'KTMAX', nc3d(3)), &
881  thisroutine )
882  call check( nf90_def_var(ncid, 'omega_t', nf90_float, nc3d, ncv), &
883  thisroutine )
884  buffer = '1'
885  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
886  thisroutine )
887  buffer = 'land_ice_temperate_layer_water_content'
888  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
889  thisroutine )
890  buffer = 'Water content in the temperate ice layer'
891  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
892  thisroutine )
893  call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
894  thisroutine )
895  call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
896  thisroutine )
897 
898  call check( nf90_inq_dimid(ncid, 'IMAX', nc3d(1)), &
899  thisroutine )
900  call check( nf90_inq_dimid(ncid, 'JMAX', nc3d(2)), &
901  thisroutine )
902  call check( nf90_inq_dimid(ncid, 'KRMAX', nc3d(3)), &
903  thisroutine )
904  call check( nf90_def_var(ncid, 'temp_r', nf90_float, nc3d, ncv), &
905  thisroutine )
906  buffer = 'degC'
907  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
908  thisroutine )
909  buffer = 'lithosphere_layer_temperature'
910  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
911  thisroutine )
912  buffer = 'Temperature in the lithosphere layer'
913  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
914  thisroutine )
915  call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
916  thisroutine )
917  call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
918  thisroutine )
919 
920 end if
921 
922 call check( nf90_inq_dimid(ncid, 'IMAX', nc2d(1)), &
923  thisroutine )
924 call check( nf90_inq_dimid(ncid, 'JMAX', nc2d(2)), &
925  thisroutine )
926 call check( nf90_def_var(ncid, 'Q_bm', nf90_float, nc2d, ncv), &
927  thisroutine )
928 buffer = 'm a-1'
929 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
930  thisroutine )
931 buffer = 'land_ice_basal_melt_rate'
932 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
933  thisroutine )
934 buffer = 'Basal melting rate'
935 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
936  thisroutine )
937 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
938  thisroutine )
939 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
940  thisroutine )
941 
942 call check( nf90_inq_dimid(ncid, 'IMAX', nc2d(1)), &
943  thisroutine )
944 call check( nf90_inq_dimid(ncid, 'JMAX', nc2d(2)), &
945  thisroutine )
946 call check( nf90_def_var(ncid, 'Q_tld', nf90_float, nc2d, ncv), &
947  thisroutine )
948 buffer = 'm a-1'
949 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
950  thisroutine )
951 buffer = 'land_ice_temperate_layer_water_drainage'
952 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
953  thisroutine )
954 buffer = 'Water drainage from the temperate layer'
955 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
956  thisroutine )
957 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
958  thisroutine )
959 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
960  thisroutine )
961 
962 call check( nf90_inq_dimid(ncid, 'IMAX', nc2d(1)), &
963  thisroutine )
964 call check( nf90_inq_dimid(ncid, 'JMAX', nc2d(2)), &
965  thisroutine )
966 call check( nf90_def_var(ncid, 'am_perp', nf90_float, nc2d, ncv), &
967  thisroutine )
968 buffer = 'm a-1'
969 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
970  thisroutine )
971 buffer = 'land_ice_volume_flux_across_cts'
972 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
973  thisroutine )
974 buffer = 'Volume flux across the CTS (cold-temperate transition surface)'
975 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
976  thisroutine )
977 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
978  thisroutine )
979 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
980  thisroutine )
981 
982 call check( nf90_inq_dimid(ncid, 'IMAX', nc2d(1)), &
983  thisroutine )
984 call check( nf90_inq_dimid(ncid, 'JMAX', nc2d(2)), &
985  thisroutine )
986 call check( nf90_def_var(ncid, 'qx', nf90_float, nc2d, ncv), &
987  thisroutine )
988 buffer = 'm2 a-1'
989 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
990  thisroutine )
991 buffer = 'land_ice_vertical_integral_x_velocity'
992 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
993  thisroutine )
994 buffer = 'Horizontal volume flux qx'
995 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
996  thisroutine )
997 buffer = 'Staggered grid variable, defined at (j,i+1/2)'
998 call check( nf90_put_att(ncid, ncv, 'comment', trim(buffer)), &
999  thisroutine )
1000 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1001  thisroutine )
1002 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1003  thisroutine )
1004 
1005 call check( nf90_inq_dimid(ncid, 'IMAX', nc2d(1)), &
1006  thisroutine )
1007 call check( nf90_inq_dimid(ncid, 'JMAX', nc2d(2)), &
1008  thisroutine )
1009 call check( nf90_def_var(ncid, 'qy', nf90_float, nc2d, ncv), &
1010  thisroutine )
1011 buffer = 'm2 a-1'
1012 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1013  thisroutine )
1014 buffer = 'land_ice_vertical_integral_y_velocity'
1015 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1016  thisroutine )
1017 buffer = 'Horizontal volume flux qy'
1018 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1019  thisroutine )
1020 buffer = 'Staggered grid variable, defined at (j+1/2,i)'
1021 call check( nf90_put_att(ncid, ncv, 'comment', trim(buffer)), &
1022  thisroutine )
1023 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1024  thisroutine )
1025 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1026  thisroutine )
1027 
1028 if (flag_3d_output) then
1029 
1030  call check( nf90_inq_dimid(ncid, 'IMAX', nc3d(1)), &
1031  thisroutine )
1032  call check( nf90_inq_dimid(ncid, 'JMAX', nc3d(2)), &
1033  thisroutine )
1034  call check( nf90_inq_dimid(ncid, 'KCMAX', nc3d(3)), &
1035  thisroutine )
1036  call check( nf90_def_var(ncid, 'age_c', nf90_float, nc3d, ncv), &
1037  thisroutine )
1038  buffer = 'a'
1039  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1040  thisroutine )
1041  buffer = 'land_ice_cold_layer_age'
1042  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1043  thisroutine )
1044  buffer = 'Age in the cold ice layer'
1045  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1046  thisroutine )
1047  call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1048  thisroutine )
1049  call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1050  thisroutine )
1051 
1052  call check( nf90_inq_dimid(ncid, 'IMAX', nc3d(1)), &
1053  thisroutine )
1054  call check( nf90_inq_dimid(ncid, 'JMAX', nc3d(2)), &
1055  thisroutine )
1056  call check( nf90_inq_dimid(ncid, 'KTMAX', nc3d(3)), &
1057  thisroutine )
1058  call check( nf90_def_var(ncid, 'age_t', nf90_float, nc3d, ncv), &
1059  thisroutine )
1060  buffer = 'a'
1061  call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1062  thisroutine )
1063  buffer = 'land_ice_temperate_layer_age'
1064  call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1065  thisroutine )
1066  buffer = 'Age in the temperate ice layer'
1067  call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1068  thisroutine )
1069  call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1070  thisroutine )
1071  call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1072  thisroutine )
1073 
1074 end if
1075 
1076 call check( nf90_inq_dimid(ncid, 'IMAX', nc2d(1)), &
1077  thisroutine )
1078 call check( nf90_inq_dimid(ncid, 'JMAX', nc2d(2)), &
1079  thisroutine )
1080 call check( nf90_def_var(ncid, 'dzs_dtau', nf90_float, nc2d, ncv), &
1081  thisroutine )
1082 buffer = 'm a-1'
1083 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1084  thisroutine )
1085 buffer = 'tendency_of_surface_altitude'
1086 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1087  thisroutine )
1088 buffer = 'Rate of change of the topography of the free surface'
1089 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1090  thisroutine )
1091 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1092  thisroutine )
1093 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1094  thisroutine )
1095 
1096 call check( nf90_inq_dimid(ncid, 'IMAX', nc2d(1)), &
1097  thisroutine )
1098 call check( nf90_inq_dimid(ncid, 'JMAX', nc2d(2)), &
1099  thisroutine )
1100 call check( nf90_def_var(ncid, 'dzm_dtau', nf90_float, nc2d, ncv), &
1101  thisroutine )
1102 buffer = 'm a-1'
1103 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1104  thisroutine )
1105 buffer = 'tendency_of_cts_altitude'
1106 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1107  thisroutine )
1108 buffer = 'Rate of change of the topography of the CTS '// &
1109  '(cold-temperate transition surface)'
1110 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1111  thisroutine )
1112 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1113  thisroutine )
1114 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1115  thisroutine )
1116 
1117 call check( nf90_inq_dimid(ncid, 'IMAX', nc2d(1)), &
1118  thisroutine )
1119 call check( nf90_inq_dimid(ncid, 'JMAX', nc2d(2)), &
1120  thisroutine )
1121 call check( nf90_def_var(ncid, 'dzb_dtau', nf90_float, nc2d, ncv), &
1122  thisroutine )
1123 buffer = 'm a-1'
1124 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1125  thisroutine )
1126 buffer = 'tendency_of_ice_base_altitude'
1127 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1128  thisroutine )
1129 buffer = 'Rate of change of the topography of the ice base'
1130 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1131  thisroutine )
1132 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1133  thisroutine )
1134 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1135  thisroutine )
1136 
1137 call check( nf90_inq_dimid(ncid, 'IMAX', nc2d(1)), &
1138  thisroutine )
1139 call check( nf90_inq_dimid(ncid, 'JMAX', nc2d(2)), &
1140  thisroutine )
1141 call check( nf90_def_var(ncid, 'dzl_dtau', nf90_float, nc2d, ncv), &
1142  thisroutine )
1143 buffer = 'm a-1'
1144 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1145  thisroutine )
1146 buffer = 'tendency_of_bedrock_altitude'
1147 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1148  thisroutine )
1149 buffer = 'Rate of change of the topography of the lithosphere surface'
1150 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1151  thisroutine )
1152 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1153  thisroutine )
1154 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1155  thisroutine )
1156 
1157 call check( nf90_inq_dimid(ncid, 'IMAX', nc2d(1)), &
1158  thisroutine )
1159 call check( nf90_inq_dimid(ncid, 'JMAX', nc2d(2)), &
1160  thisroutine )
1161 call check( nf90_def_var(ncid, 'dH_c_dtau', nf90_float, nc2d, ncv), &
1162  thisroutine )
1163 buffer = 'm a-1'
1164 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1165  thisroutine )
1166 buffer = 'tendency_of_land_ice_cold_layer_thickness'
1167 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1168  thisroutine )
1169 buffer = 'Rate of change of the thickness of the cold ice layer'
1170 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1171  thisroutine )
1172 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1173  thisroutine )
1174 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1175  thisroutine )
1176 
1177 call check( nf90_inq_dimid(ncid, 'IMAX', nc2d(1)), &
1178  thisroutine )
1179 call check( nf90_inq_dimid(ncid, 'JMAX', nc2d(2)), &
1180  thisroutine )
1181 call check( nf90_def_var(ncid, 'dH_t_dtau', nf90_float, nc2d, ncv), &
1182  thisroutine )
1183 buffer = 'm a-1'
1184 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1185  thisroutine )
1186 buffer = 'tendency_of_land_ice_temperate_layer_thickness'
1187 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1188  thisroutine )
1189 buffer = 'Rate of change of the thickness of the temperate ice layer'
1190 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1191  thisroutine )
1192 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1193  thisroutine )
1194 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1195  thisroutine )
1196 
1197 call check( nf90_inq_dimid(ncid, 'IMAX', nc2d(1)), &
1198  thisroutine )
1199 call check( nf90_inq_dimid(ncid, 'JMAX', nc2d(2)), &
1200  thisroutine )
1201 call check( nf90_def_var(ncid, 'dH_dtau', nf90_float, nc2d, ncv), &
1202  thisroutine )
1203 buffer = 'm a-1'
1204 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1205  thisroutine )
1206 buffer = 'tendency_of_land_ice_thickness'
1207 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1208  thisroutine )
1209 buffer = 'Rate of change of the ice thickness'
1210 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1211  thisroutine )
1212 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1213  thisroutine )
1214 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1215  thisroutine )
1216 
1217 call check( nf90_inq_dimid(ncid, 'IMAX', nc2d(1)), &
1218  thisroutine )
1219 call check( nf90_inq_dimid(ncid, 'JMAX', nc2d(2)), &
1220  thisroutine )
1221 call check( nf90_def_var(ncid, 'vx_b_g', nf90_float, nc2d, ncv), &
1222  thisroutine )
1223 buffer = 'm a-1'
1224 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1225  thisroutine )
1226 buffer = 'land_ice_base_x_velocity'
1227 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1228  thisroutine )
1229 buffer = 'Horizontal velocity vx at the ice base'
1230 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1231  thisroutine )
1232 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1233  thisroutine )
1234 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1235  thisroutine )
1236 
1237 call check( nf90_inq_dimid(ncid, 'IMAX', nc2d(1)), &
1238  thisroutine )
1239 call check( nf90_inq_dimid(ncid, 'JMAX', nc2d(2)), &
1240  thisroutine )
1241 call check( nf90_def_var(ncid, 'vy_b_g', nf90_float, nc2d, ncv), &
1242  thisroutine )
1243 buffer = 'm a-1'
1244 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1245  thisroutine )
1246 buffer = 'land_ice_base_y_velocity'
1247 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1248  thisroutine )
1249 buffer = 'Horizontal velocity vy at the ice base'
1250 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1251  thisroutine )
1252 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1253  thisroutine )
1254 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1255  thisroutine )
1256 
1257 call check( nf90_inq_dimid(ncid, 'IMAX', nc2d(1)), &
1258  thisroutine )
1259 call check( nf90_inq_dimid(ncid, 'JMAX', nc2d(2)), &
1260  thisroutine )
1261 call check( nf90_def_var(ncid, 'vz_b', nf90_float, nc2d, ncv), &
1262  thisroutine )
1263 buffer = 'm a-1'
1264 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1265  thisroutine )
1266 buffer = 'land_ice_base_z_velocity'
1267 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1268  thisroutine )
1269 buffer = 'Vertical velocity vz at the ice base'
1270 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1271  thisroutine )
1272 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1273  thisroutine )
1274 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1275  thisroutine )
1276 
1277 call check( nf90_inq_dimid(ncid, 'IMAX', nc2d(1)), &
1278  thisroutine )
1279 call check( nf90_inq_dimid(ncid, 'JMAX', nc2d(2)), &
1280  thisroutine )
1281 call check( nf90_def_var(ncid, 'vh_b', nf90_float, nc2d, ncv), &
1282  thisroutine )
1283 buffer = 'm a-1'
1284 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1285  thisroutine )
1286 buffer = 'land_ice_base_horizontal_velocity'
1287 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1288  thisroutine )
1289 buffer = 'Horizontal velocity vh at the ice base'
1290 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1291  thisroutine )
1292 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1293  thisroutine )
1294 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1295  thisroutine )
1296 
1297 call check( nf90_inq_dimid(ncid, 'IMAX', nc2d(1)), &
1298  thisroutine )
1299 call check( nf90_inq_dimid(ncid, 'JMAX', nc2d(2)), &
1300  thisroutine )
1301 call check( nf90_def_var(ncid, 'vx_s_g', nf90_float, nc2d, ncv), &
1302  thisroutine )
1303 buffer = 'm a-1'
1304 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1305  thisroutine )
1306 buffer = 'land_ice_surface_x_velocity'
1307 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1308  thisroutine )
1309 buffer = 'Horizontal velocity vx at the ice surface'
1310 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1311  thisroutine )
1312 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1313  thisroutine )
1314 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1315  thisroutine )
1316 
1317 call check( nf90_inq_dimid(ncid, 'IMAX', nc2d(1)), &
1318  thisroutine )
1319 call check( nf90_inq_dimid(ncid, 'JMAX', nc2d(2)), &
1320  thisroutine )
1321 call check( nf90_def_var(ncid, 'vy_s_g', nf90_float, nc2d, ncv), &
1322  thisroutine )
1323 buffer = 'm a-1'
1324 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1325  thisroutine )
1326 buffer = 'land_ice_surface_y_velocity'
1327 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1328  thisroutine )
1329 buffer = 'Horizontal velocity vy at the ice surface'
1330 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1331  thisroutine )
1332 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1333  thisroutine )
1334 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1335  thisroutine )
1336 
1337 call check( nf90_inq_dimid(ncid, 'IMAX', nc2d(1)), &
1338  thisroutine )
1339 call check( nf90_inq_dimid(ncid, 'JMAX', nc2d(2)), &
1340  thisroutine )
1341 call check( nf90_def_var(ncid, 'vz_s', nf90_float, nc2d, ncv), &
1342  thisroutine )
1343 buffer = 'm a-1'
1344 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1345  thisroutine )
1346 buffer = 'land_ice_surface_z_velocity'
1347 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1348  thisroutine )
1349 buffer = 'Vertical velocity vz at the ice surface'
1350 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1351  thisroutine )
1352 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1353  thisroutine )
1354 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1355  thisroutine )
1356 
1357 call check( nf90_inq_dimid(ncid, 'IMAX', nc2d(1)), &
1358  thisroutine )
1359 call check( nf90_inq_dimid(ncid, 'JMAX', nc2d(2)), &
1360  thisroutine )
1361 call check( nf90_def_var(ncid, 'vh_s', nf90_float, nc2d, ncv), &
1362  thisroutine )
1363 buffer = 'm a-1'
1364 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1365  thisroutine )
1366 buffer = 'land_ice_surface_horizontal_velocity'
1367 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1368  thisroutine )
1369 buffer = 'Horizontal velocity vh at the ice surface'
1370 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1371  thisroutine )
1372 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1373  thisroutine )
1374 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1375  thisroutine )
1376 
1377 call check( nf90_inq_dimid(ncid, 'IMAX', nc2d(1)), &
1378  thisroutine )
1379 call check( nf90_inq_dimid(ncid, 'JMAX', nc2d(2)), &
1380  thisroutine )
1381 call check( nf90_def_var(ncid, 'temp_b', nf90_float, nc2d, ncv), &
1382  thisroutine )
1383 buffer = 'degC'
1384 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1385  thisroutine )
1386 buffer = 'basal_temperature'
1387 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1388  thisroutine )
1389 buffer = 'Temperature at the ice base'
1390 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1391  thisroutine )
1392 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1393  thisroutine )
1394 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1395  thisroutine )
1396 
1397 call check( nf90_inq_dimid(ncid, 'IMAX', nc2d(1)), &
1398  thisroutine )
1399 call check( nf90_inq_dimid(ncid, 'JMAX', nc2d(2)), &
1400  thisroutine )
1401 call check( nf90_def_var(ncid, 'temph_b', nf90_float, nc2d, ncv), &
1402  thisroutine )
1403 buffer = 'degC'
1404 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1405  thisroutine )
1406 buffer = 'basal_temperature_rel_to_pmp'
1407 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1408  thisroutine )
1409 buffer = 'Temperature at the ice base relative to the pressure melting point'
1410 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1411  thisroutine )
1412 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1413  thisroutine )
1414 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1415  thisroutine )
1416 
1417 call check( nf90_inq_dimid(ncid, 'IMAX', nc2d(1)), &
1418  thisroutine )
1419 call check( nf90_inq_dimid(ncid, 'JMAX', nc2d(2)), &
1420  thisroutine )
1421 call check( nf90_def_var(ncid, 'p_b_w', nf90_float, nc2d, ncv), &
1422  thisroutine )
1423 buffer = 'Pa'
1424 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1425  thisroutine )
1426 buffer = 'basal_water_pressure'
1427 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1428  thisroutine )
1429 buffer = 'Basal water pressure'
1430 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1431  thisroutine )
1432 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1433  thisroutine )
1434 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1435  thisroutine )
1436 
1437 call check( nf90_inq_dimid(ncid, 'IMAX', nc2d(1)), &
1438  thisroutine )
1439 call check( nf90_inq_dimid(ncid, 'JMAX', nc2d(2)), &
1440  thisroutine )
1441 call check( nf90_def_var(ncid, 'H_w', nf90_float, nc2d, ncv), &
1442  thisroutine )
1443 buffer = 'm'
1444 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1445  thisroutine )
1446 buffer = 'subglacial_water_thickness'
1447 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1448  thisroutine )
1449 buffer = 'Effective thickness of subglacial water'
1450 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1451  thisroutine )
1452 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1453  thisroutine )
1454 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1455  thisroutine )
1456 
1457 call check( nf90_inq_dimid(ncid, 'IMAX', nc2d(1)), &
1458  thisroutine )
1459 call check( nf90_inq_dimid(ncid, 'JMAX', nc2d(2)), &
1460  thisroutine )
1461 call check( nf90_def_var(ncid, 'q_gl_g', nf90_float, nc2d, ncv), &
1462  thisroutine )
1463 buffer = 'm2 a-1'
1464 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1465  thisroutine )
1466 buffer = 'land_ice_volume_flux_across_gl'
1467 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1468  thisroutine )
1469 buffer = 'Horizontal volume flux across the grounding line'
1470 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1471  thisroutine )
1472 call check( nf90_put_att(ncid, ncv, 'grid_mapping', 'crs'), &
1473  thisroutine )
1474 call check( nf90_put_att(ncid, ncv, 'coordinates', 'lat lon'), &
1475  thisroutine )
1476 
1477 call check( nf90_def_var(ncid, 'V_tot', nf90_float, ncv), &
1478  thisroutine )
1479 buffer = 'm3'
1480 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1481  thisroutine )
1482 buffer = 'land_ice_volume'
1483 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1484  thisroutine )
1485 buffer = 'Ice volume'
1486 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1487  thisroutine )
1488 
1489 call check( nf90_def_var(ncid, 'A_grounded', nf90_float, ncv), &
1490  thisroutine )
1491 buffer = 'm2'
1492 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1493  thisroutine )
1494 buffer = 'land_ice_area_grounded'
1495 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1496  thisroutine )
1497 buffer = 'Area covered by grounded ice'
1498 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1499  thisroutine )
1500 
1501 call check( nf90_def_var(ncid, 'A_floating', nf90_float, ncv), &
1502  thisroutine )
1503 buffer = 'm2'
1504 call check( nf90_put_att(ncid, ncv, 'units', trim(buffer)), &
1505  thisroutine )
1506 buffer = 'land_ice_area_floating'
1507 call check( nf90_put_att(ncid, ncv, 'standard_name', trim(buffer)), &
1508  thisroutine )
1509 buffer = 'Area covered by floating ice'
1510 call check( nf90_put_att(ncid, ncv, 'long_name', trim(buffer)), &
1511  thisroutine )
1512 
1513 call check( nf90_enddef(ncid), thisroutine )
1514 
1515 !-------- Computation of the ice volume,
1516 ! the area covered by grounded ice
1517 ! and the area covered by floating ice --------
1518 
1519 v_tot = 0.0_dp
1520 a_grounded = 0.0_dp
1521 a_floating = 0.0_dp
1522 
1523 do i=0, imax
1524 do j=0, jmax
1525 
1526  if ( (maske(j,i)==0).or.(maske(j,i)==3) ) & ! grounded or floating ice
1527  v_tot = v_tot + (h_c(j,i)+h_t(j,i))*area(j,i)
1528 
1529  if (maske(j,i)==0) & ! grounded ice
1530  a_grounded = a_grounded + area(j,i)
1531 
1532  if (maske(j,i)==3) & ! floating ice
1533  a_floating = a_floating + area(j,i)
1534 
1535 end do
1536 end do
1537 
1538 !-------- Convert data to real*4 and seconds to years --------
1539 
1540 #if ( !defined(OUT_TIMES) || OUT_TIMES==1 )
1541 time_conv = real(time/year_sec,sp)
1542 #elif OUT_TIMES == 2
1543 time_conv = real((time+year_zero)/year_sec,sp)
1544 #else
1545 stop ' output_nc: OUT_TIMES must be either 1 or 2!'
1546 #endif
1547 
1548 delta_ts_conv = real(delta_ts,sp)
1549 glac_index_conv = real(glac_index,sp)
1550 z_sl_conv = real(z_sl,sp)
1551 v_tot_conv = real(v_tot,sp)
1552 a_grounded_conv = real(a_grounded,sp)
1553 a_floating_conv = real(a_floating,sp)
1554 h_r_conv = real(h_r,sp)
1555 
1556 do i=0, imax
1557  xi_conv(i) = real(xi(i),sp)
1558 end do
1559 
1560 do j=0, jmax
1561  eta_conv(j) = real(eta(j),sp)
1562 end do
1563 
1564 do kc=0, kcmax
1565  sigma_level_c_conv(kc) = real(eaz_c_quotient(kc),sp)
1566 end do
1567 
1568 do kt=0, ktmax
1569  sigma_level_t_conv(kt) = real(zeta_t(kt),sp)
1570 end do
1571 
1572 do kr=0, krmax
1573  sigma_level_r_conv(kr) = real(kr,sp)/real(krmax,sp)
1574 end do
1575 
1576 do i=0, imax
1577 do j=0, jmax
1578 
1579  maske_conv(i,j) = maske(j,i)
1580  n_cts_conv(i,j) = n_cts(j,i)
1581 
1582  lambda_conv(i,j) = real(lambda(j,i),sp)
1583  phi_conv(i,j) = real(phi(j,i),sp)
1584  lon_conv(i,j) = real(lambda(j,i)*pi_180_inv,sp) ! longitude in deg
1585  lon_conv(i,j) = modulo(lon_conv(i,j)+180.0_sp, 360.0_sp)-180.0_sp
1586  ! mapping to interval [-180 deg, +180 deg)
1587  lat_conv(i,j) = real(phi(j,i) *pi_180_inv,sp) ! latitute in deg
1588  if (lat_conv(i,j) > 90.0_sp) lat_conv(i,j) = 90.0_sp
1589  if (lat_conv(i,j) < -90.0_sp) lat_conv(i,j) = -90.0_sp
1590  ! constraining to interval [-90 deg, +90 deg]
1591  temp_s_conv(i,j) = real(temp_s(j,i),sp)
1592  as_perp_conv(i,j) = real(as_perp(j,i)*year_sec,sp)
1593  zs_conv(i,j) = real(zs(j,i),sp)
1594  zm_conv(i,j) = real(zm(j,i),sp)
1595  zb_conv(i,j) = real(zb(j,i),sp)
1596  zl_conv(i,j) = real(zl(j,i),sp)
1597  h_c_conv(i,j) = real(H_c(j,i),sp)
1598  h_t_conv(i,j) = real(H_t(j,i),sp)
1599  h_conv(i,j) = real(H_c(j,i)+H_t(j,i),sp)
1600  q_bm_conv(i,j) = real(q_bm(j,i)*year_sec,sp)
1601  q_tld_conv(i,j) = real(q_tld(j,i)*year_sec,sp)
1602  am_perp_conv(i,j) = real(am_perp(j,i)*year_sec,sp)
1603  qx_conv(i,j) = real(qx(j,i)*year_sec,sp)
1604  qy_conv(i,j) = real(qy(j,i)*year_sec,sp)
1605  dzs_dtau_conv(i,j) = real(dzs_dtau(j,i)*year_sec,sp)
1606  dzm_dtau_conv(i,j) = real(dzm_dtau(j,i)*year_sec,sp)
1607  dzb_dtau_conv(i,j) = real(dzb_dtau(j,i)*year_sec,sp)
1608  dzl_dtau_conv(i,j) = real(dzl_dtau(j,i)*year_sec,sp)
1609  dh_c_dtau_conv(i,j) = real(dh_c_dtau(j,i)*year_sec,sp)
1610  dh_t_dtau_conv(i,j) = real(dh_t_dtau(j,i)*year_sec,sp)
1611  dh_dtau_conv(i,j) = real((dh_c_dtau(j,i)+dh_t_dtau(j,i))*year_sec,sp)
1612  vx_b_g_conv(i,j) = real(vx_b_g(j,i)*year_sec,sp)
1613  vy_b_g_conv(i,j) = real(vy_b_g(j,i)*year_sec,sp)
1614  vz_b_conv(i,j) = real(vz_b(j,i)*year_sec,sp)
1615  vh_b_conv(i,j) = sqrt( vx_b_g_conv(i,j)**2 + vy_b_g_conv(i,j)**2 )
1616  vx_s_g_conv(i,j) = real(vx_s_g(j,i)*year_sec,sp)
1617  vy_s_g_conv(i,j) = real(vy_s_g(j,i)*year_sec,sp)
1618  vz_s_conv(i,j) = real(vz_s(j,i)*year_sec,sp)
1619  vh_s_conv(i,j) = sqrt( vx_s_g_conv(i,j)**2 + vy_s_g_conv(i,j)**2 )
1620  temp_b_conv(i,j) = real(temp_b(j,i),sp)
1621  temph_b_conv(i,j) = real(temph_b(j,i),sp)
1622  p_b_w_conv(i,j) = real(p_b_w(j,i),sp)
1623  h_w_conv(i,j) = real(H_w(j,i),sp)
1624  q_gl_g_conv(i,j) = real(q_gl_g(j,i)*year_sec,sp)
1625 
1626  do kr=0, krmax
1627  temp_r_conv(i,j,kr) = real(temp_r(kr,j,i),sp)
1628  end do
1629 
1630  do kt=0, ktmax
1631  vx_t_conv(i,j,kt) = real(vx_t(kt,j,i)*year_sec,sp)
1632  vy_t_conv(i,j,kt) = real(vy_t(kt,j,i)*year_sec,sp)
1633  vz_t_conv(i,j,kt) = real(vz_t(kt,j,i)*year_sec,sp)
1634  omega_t_conv(i,j,kt) = real(omega_t(kt,j,i),sp)
1635  age_t_conv(i,j,kt) = real(age_t(kt,j,i)/year_sec,sp)
1636  end do
1637 
1638  do kc=0, kcmax
1639  vx_c_conv(i,j,kc) = real(vx_c(kc,j,i)*year_sec,sp)
1640  vy_c_conv(i,j,kc) = real(vy_c(kc,j,i)*year_sec,sp)
1641  vz_c_conv(i,j,kc) = real(vz_c(kc,j,i)*year_sec,sp)
1642  temp_c_conv(i,j,kc) = real(temp_c(kc,j,i),sp)
1643  age_c_conv(i,j,kc) = real(age_c(kc,j,i)/year_sec,sp)
1644  end do
1645 
1646 end do
1647 end do
1648 
1649 !-------- Write data on file --------
1650 
1651 call check( nf90_inq_varid(ncid, 'crs', ncv), thisroutine )
1652 call check( nf90_put_var(ncid, ncv, 0), thisroutine )
1653 
1654 call check( nf90_inq_varid(ncid, 'time', ncv), thisroutine )
1655 call check( nf90_put_var(ncid, ncv, time_conv), thisroutine )
1656 
1657 if (forcing_flag == 1) then
1658 
1659  call check( nf90_inq_varid(ncid, 'delta_ts', ncv), thisroutine )
1660  call check( nf90_put_var(ncid, ncv, delta_ts_conv), thisroutine )
1661 
1662 else if (forcing_flag == 2) then
1663 
1664  call check( nf90_inq_varid(ncid, 'glac_index', ncv), thisroutine )
1665  call check( nf90_put_var(ncid, ncv, glac_index_conv), thisroutine )
1666 
1667 else if (forcing_flag == 3) then
1668 
1669  call check( nf90_inq_varid(ncid, 'glac_index', ncv), thisroutine )
1670  glac_index_conv = 1.11e+11 ! dummy value
1671  call check( nf90_put_var(ncid, ncv, glac_index_conv), thisroutine )
1672 
1673 end if
1674 
1675 call check( nf90_inq_varid(ncid, 'z_sl', ncv), thisroutine )
1676 call check( nf90_put_var(ncid, ncv, z_sl_conv), thisroutine )
1677 
1678 call check( nf90_inq_varid(ncid, 'xi', ncv), thisroutine )
1679 nc1cor(1) = 1
1680 nc1cnt(1) = imax + 1
1681 call check( nf90_put_var(ncid, ncv, xi_conv, start=nc1cor, count=nc1cnt), &
1682  thisroutine )
1683 
1684 call check( nf90_inq_varid(ncid, 'eta', ncv), thisroutine )
1685 nc1cor(1) = 1
1686 nc1cnt(1) = jmax + 1
1687 call check( nf90_put_var(ncid, ncv, eta_conv, start=nc1cor, count=nc1cnt), &
1688  thisroutine )
1689 
1690 call check( nf90_inq_varid(ncid, 'sigma_level_c', ncv), thisroutine )
1691 nc1cor(1) = 1
1692 nc1cnt(1) = kcmax + 1
1693 call check( nf90_put_var(ncid, ncv, sigma_level_c_conv, start=nc1cor, count=nc1cnt), &
1694  thisroutine )
1695 
1696 call check( nf90_inq_varid(ncid, 'sigma_level_t', ncv), thisroutine )
1697 nc1cor(1) = 1
1698 nc1cnt(1) = ktmax + 1
1699 call check( nf90_put_var(ncid, ncv, sigma_level_t_conv, start=nc1cor, count=nc1cnt), &
1700  thisroutine )
1701 
1702 call check( nf90_inq_varid(ncid, 'sigma_level_r', ncv), thisroutine )
1703 nc1cor(1) = 1
1704 nc1cnt(1) = krmax + 1
1705 call check( nf90_put_var(ncid, ncv, sigma_level_r_conv, start=nc1cor, count=nc1cnt), &
1706  thisroutine )
1707 
1708 call check( nf90_inq_varid(ncid, 'lon', ncv), thisroutine )
1709 nc2cor(1) = 1
1710 nc2cor(2) = 1
1711 nc2cnt(1) = imax + 1
1712 nc2cnt(2) = jmax + 1
1713 call check( nf90_put_var(ncid, ncv, lon_conv, start=nc2cor, count=nc2cnt), &
1714  thisroutine )
1715 
1716 call check( nf90_inq_varid(ncid, 'lat', ncv), thisroutine )
1717 nc2cor(1) = 1
1718 nc2cor(2) = 1
1719 nc2cnt(1) = imax + 1
1720 nc2cnt(2) = jmax + 1
1721 call check( nf90_put_var(ncid, ncv, lat_conv, start=nc2cor, count=nc2cnt), &
1722  thisroutine )
1723 
1724 call check( nf90_inq_varid(ncid, 'lambda', ncv), thisroutine )
1725 nc2cor(1) = 1
1726 nc2cor(2) = 1
1727 nc2cnt(1) = imax + 1
1728 nc2cnt(2) = jmax + 1
1729 call check( nf90_put_var(ncid, ncv, lambda_conv, start=nc2cor, count=nc2cnt), &
1730  thisroutine )
1731 
1732 call check( nf90_inq_varid(ncid, 'phi', ncv), thisroutine )
1733 nc2cor(1) = 1
1734 nc2cor(2) = 1
1735 nc2cnt(1) = imax + 1
1736 nc2cnt(2) = jmax + 1
1737 call check( nf90_put_var(ncid, ncv, phi_conv, start=nc2cor, count=nc2cnt), &
1738  thisroutine )
1739 
1740 call check( nf90_inq_varid(ncid, 'temp_s', ncv), thisroutine )
1741 nc2cor(1) = 1
1742 nc2cor(2) = 1
1743 nc2cnt(1) = imax + 1
1744 nc2cnt(2) = jmax + 1
1745 call check( nf90_put_var(ncid, ncv, temp_s_conv, start=nc2cor, count=nc2cnt), &
1746  thisroutine )
1747 
1748 call check( nf90_inq_varid(ncid, 'as_perp', ncv), thisroutine )
1749 nc2cor(1) = 1
1750 nc2cor(2) = 1
1751 nc2cnt(1) = imax + 1
1752 nc2cnt(2) = jmax + 1
1753 call check( nf90_put_var(ncid, ncv, as_perp_conv, start=nc2cor, count=nc2cnt), &
1754  thisroutine )
1755 
1756 call check( nf90_inq_varid(ncid, 'maske', ncv), thisroutine )
1757 nc2cor(1) = 1
1758 nc2cor(2) = 1
1759 nc2cnt(1) = imax + 1
1760 nc2cnt(2) = jmax + 1
1761 call check( nf90_put_var(ncid, ncv, maske_conv, start=nc2cor, count=nc2cnt), &
1762  thisroutine )
1763 
1764 call check( nf90_inq_varid(ncid, 'n_cts', ncv), thisroutine )
1765 nc2cor(1) = 1
1766 nc2cor(2) = 1
1767 nc2cnt(1) = imax + 1
1768 nc2cnt(2) = jmax + 1
1769 call check( nf90_put_var(ncid, ncv, n_cts_conv, start=nc2cor, count=nc2cnt), &
1770  thisroutine )
1771 
1772 call check( nf90_inq_varid(ncid, 'zs', ncv), thisroutine )
1773 nc2cor(1) = 1
1774 nc2cor(2) = 1
1775 nc2cnt(1) = imax + 1
1776 nc2cnt(2) = jmax + 1
1777 call check( nf90_put_var(ncid, ncv, zs_conv, start=nc2cor, count=nc2cnt), &
1778  thisroutine )
1779 
1780 call check( nf90_inq_varid(ncid, 'zm', ncv), thisroutine )
1781 nc2cor(1) = 1
1782 nc2cor(2) = 1
1783 nc2cnt(1) = imax + 1
1784 nc2cnt(2) = jmax + 1
1785 call check( nf90_put_var(ncid, ncv, zm_conv, start=nc2cor, count=nc2cnt), &
1786  thisroutine )
1787 
1788 call check( nf90_inq_varid(ncid, 'zb', ncv), thisroutine )
1789 nc2cor(1) = 1
1790 nc2cor(2) = 1
1791 nc2cnt(1) = imax + 1
1792 nc2cnt(2) = jmax + 1
1793 call check( nf90_put_var(ncid, ncv, zb_conv, start=nc2cor, count=nc2cnt), &
1794  thisroutine )
1795 
1796 call check( nf90_inq_varid(ncid, 'zl', ncv), thisroutine )
1797 nc2cor(1) = 1
1798 nc2cor(2) = 1
1799 nc2cnt(1) = imax + 1
1800 nc2cnt(2) = jmax + 1
1801 call check( nf90_put_var(ncid, ncv, zl_conv, start=nc2cor, count=nc2cnt), &
1802  thisroutine )
1803 
1804 call check( nf90_inq_varid(ncid, 'H_c', ncv), thisroutine )
1805 nc2cor(1) = 1
1806 nc2cor(2) = 1
1807 nc2cnt(1) = imax + 1
1808 nc2cnt(2) = jmax + 1
1809 call check( nf90_put_var(ncid, ncv, h_c_conv, start=nc2cor, count=nc2cnt), &
1810  thisroutine )
1811 
1812 call check( nf90_inq_varid(ncid, 'H_t', ncv), thisroutine )
1813 nc2cor(1) = 1
1814 nc2cor(2) = 1
1815 nc2cnt(1) = imax + 1
1816 nc2cnt(2) = jmax + 1
1817 call check( nf90_put_var(ncid, ncv, h_t_conv, start=nc2cor, count=nc2cnt), &
1818  thisroutine )
1819 
1820 call check( nf90_inq_varid(ncid, 'H', ncv), thisroutine )
1821 nc2cor(1) = 1
1822 nc2cor(2) = 1
1823 nc2cnt(1) = imax + 1
1824 nc2cnt(2) = jmax + 1
1825 call check( nf90_put_var(ncid, ncv, h_conv, start=nc2cor, count=nc2cnt), &
1826  thisroutine )
1827 
1828 call check( nf90_inq_varid(ncid, 'H_R', ncv), thisroutine )
1829 call check( nf90_put_var(ncid, ncv, h_r_conv), thisroutine )
1830 
1831 if (flag_3d_output) then
1832 
1833  call check( nf90_inq_varid(ncid, 'vx_c', ncv), thisroutine )
1834  nc3cor(1) = 1
1835  nc3cor(2) = 1
1836  nc3cor(3) = 1
1837  nc3cnt(1) = imax + 1
1838  nc3cnt(2) = jmax + 1
1839  nc3cnt(3) = kcmax + 1
1840  call check( nf90_put_var(ncid, ncv, vx_c_conv, start=nc3cor, count=nc3cnt), &
1841  thisroutine )
1842 
1843  call check( nf90_inq_varid(ncid, 'vy_c', ncv), thisroutine )
1844  nc3cor(1) = 1
1845  nc3cor(2) = 1
1846  nc3cor(3) = 1
1847  nc3cnt(1) = imax + 1
1848  nc3cnt(2) = jmax + 1
1849  nc3cnt(3) = kcmax + 1
1850  call check( nf90_put_var(ncid, ncv, vy_c_conv, start=nc3cor, count=nc3cnt), &
1851  thisroutine )
1852 
1853  call check( nf90_inq_varid(ncid, 'vz_c', ncv), thisroutine )
1854  nc3cor(1) = 1
1855  nc3cor(2) = 1
1856  nc3cor(3) = 1
1857  nc3cnt(1) = imax + 1
1858  nc3cnt(2) = jmax + 1
1859  nc3cnt(3) = kcmax + 1
1860  call check( nf90_put_var(ncid, ncv, vz_c_conv, start=nc3cor, count=nc3cnt), &
1861  thisroutine )
1862 
1863  call check( nf90_inq_varid(ncid, 'vx_t', ncv), thisroutine )
1864  nc3cor(1) = 1
1865  nc3cor(2) = 1
1866  nc3cor(3) = 1
1867  nc3cnt(1) = imax + 1
1868  nc3cnt(2) = jmax + 1
1869  nc3cnt(3) = ktmax + 1
1870  call check( nf90_put_var(ncid, ncv, vx_t_conv, start=nc3cor, count=nc3cnt), &
1871  thisroutine )
1872 
1873  call check( nf90_inq_varid(ncid, 'vy_t', ncv), thisroutine )
1874  nc3cor(1) = 1
1875  nc3cor(2) = 1
1876  nc3cor(3) = 1
1877  nc3cnt(1) = imax + 1
1878  nc3cnt(2) = jmax + 1
1879  nc3cnt(3) = ktmax + 1
1880  call check( nf90_put_var(ncid, ncv, vy_t_conv, start=nc3cor, count=nc3cnt), &
1881  thisroutine )
1882 
1883  call check( nf90_inq_varid(ncid, 'vz_t', ncv), thisroutine )
1884  nc3cor(1) = 1
1885  nc3cor(2) = 1
1886  nc3cor(3) = 1
1887  nc3cnt(1) = imax + 1
1888  nc3cnt(2) = jmax + 1
1889  nc3cnt(3) = ktmax + 1
1890  call check( nf90_put_var(ncid, ncv, vz_t_conv, start=nc3cor, count=nc3cnt), &
1891  thisroutine )
1892 
1893  call check( nf90_inq_varid(ncid, 'temp_c', ncv), thisroutine )
1894  nc3cor(1) = 1
1895  nc3cor(2) = 1
1896  nc3cor(3) = 1
1897  nc3cnt(1) = imax + 1
1898  nc3cnt(2) = jmax + 1
1899  nc3cnt(3) = kcmax + 1
1900  call check( nf90_put_var(ncid, ncv, temp_c_conv, start=nc3cor, count=nc3cnt), &
1901  thisroutine )
1902 
1903  call check( nf90_inq_varid(ncid, 'omega_t', ncv), thisroutine )
1904  nc3cor(1) = 1
1905  nc3cor(2) = 1
1906  nc3cor(3) = 1
1907  nc3cnt(1) = imax + 1
1908  nc3cnt(2) = jmax + 1
1909  nc3cnt(3) = ktmax + 1
1910  call check( nf90_put_var(ncid, ncv, omega_t_conv, start=nc3cor, count=nc3cnt), &
1911  thisroutine )
1912 
1913  call check( nf90_inq_varid(ncid, 'temp_r', ncv), thisroutine )
1914  nc3cor(1) = 1
1915  nc3cor(2) = 1
1916  nc3cor(3) = 1
1917  nc3cnt(1) = imax + 1
1918  nc3cnt(2) = jmax + 1
1919  nc3cnt(3) = krmax + 1
1920  call check( nf90_put_var(ncid, ncv, temp_r_conv, start=nc3cor, count=nc3cnt), &
1921  thisroutine )
1922 
1923 end if
1924 
1925 call check( nf90_inq_varid(ncid, 'Q_bm', ncv), thisroutine )
1926 nc2cor(1) = 1
1927 nc2cor(2) = 1
1928 nc2cnt(1) = imax + 1
1929 nc2cnt(2) = jmax + 1
1930 call check( nf90_put_var(ncid, ncv, q_bm_conv, start=nc2cor, count=nc2cnt), &
1931  thisroutine )
1932 
1933 call check( nf90_inq_varid(ncid, 'Q_tld', ncv), thisroutine )
1934 nc2cor(1) = 1
1935 nc2cor(2) = 1
1936 nc2cnt(1) = imax + 1
1937 nc2cnt(2) = jmax + 1
1938 call check( nf90_put_var(ncid, ncv, q_tld_conv, start=nc2cor, count=nc2cnt), &
1939  thisroutine )
1940 
1941 call check( nf90_inq_varid(ncid, 'am_perp', ncv), thisroutine )
1942 nc2cor(1) = 1
1943 nc2cor(2) = 1
1944 nc2cnt(1) = imax + 1
1945 nc2cnt(2) = jmax + 1
1946 call check( nf90_put_var(ncid, ncv, am_perp_conv, start=nc2cor, count=nc2cnt), &
1947  thisroutine )
1948 
1949 call check( nf90_inq_varid(ncid, 'qx', ncv), thisroutine )
1950 nc2cor(1) = 1
1951 nc2cor(2) = 1
1952 nc2cnt(1) = imax + 1
1953 nc2cnt(2) = jmax + 1
1954 call check( nf90_put_var(ncid, ncv, qx_conv, start=nc2cor, count=nc2cnt), &
1955  thisroutine )
1956 
1957 call check( nf90_inq_varid(ncid, 'qy', ncv), thisroutine )
1958 nc2cor(1) = 1
1959 nc2cor(2) = 1
1960 nc2cnt(1) = imax + 1
1961 nc2cnt(2) = jmax + 1
1962 call check( nf90_put_var(ncid, ncv, qy_conv, start=nc2cor, count=nc2cnt), &
1963  thisroutine )
1964 
1965 if (flag_3d_output) then
1966 
1967  call check( nf90_inq_varid(ncid, 'age_c', ncv), thisroutine )
1968  nc3cor(1) = 1
1969  nc3cor(2) = 1
1970  nc3cor(3) = 1
1971  nc3cnt(1) = imax + 1
1972  nc3cnt(2) = jmax + 1
1973  nc3cnt(3) = kcmax + 1
1974  call check( nf90_put_var(ncid, ncv, age_c_conv, start=nc3cor, count=nc3cnt), &
1975  thisroutine )
1976 
1977  call check( nf90_inq_varid(ncid, 'age_t', ncv), thisroutine )
1978  nc3cor(1) = 1
1979  nc3cor(2) = 1
1980  nc3cor(3) = 1
1981  nc3cnt(1) = imax + 1
1982  nc3cnt(2) = jmax + 1
1983  nc3cnt(3) = ktmax + 1
1984  call check( nf90_put_var(ncid, ncv, age_t_conv, start=nc3cor, count=nc3cnt), &
1985  thisroutine )
1986 
1987 end if
1988 
1989 call check( nf90_inq_varid(ncid, 'dzs_dtau', ncv), thisroutine )
1990 nc2cor(1) = 1
1991 nc2cor(2) = 1
1992 nc2cnt(1) = imax + 1
1993 nc2cnt(2) = jmax + 1
1994 call check( nf90_put_var(ncid, ncv, dzs_dtau_conv, start=nc2cor, count=nc2cnt), &
1995  thisroutine )
1996 
1997 call check( nf90_inq_varid(ncid, 'dzm_dtau', ncv), thisroutine )
1998 nc2cor(1) = 1
1999 nc2cor(2) = 1
2000 nc2cnt(1) = imax + 1
2001 nc2cnt(2) = jmax + 1
2002 call check( nf90_put_var(ncid, ncv, dzm_dtau_conv, start=nc2cor, count=nc2cnt), &
2003  thisroutine )
2004 
2005 call check( nf90_inq_varid(ncid, 'dzb_dtau', ncv), thisroutine )
2006 nc2cor(1) = 1
2007 nc2cor(2) = 1
2008 nc2cnt(1) = imax + 1
2009 nc2cnt(2) = jmax + 1
2010 call check( nf90_put_var(ncid, ncv, dzb_dtau_conv, start=nc2cor, count=nc2cnt), &
2011  thisroutine )
2012 
2013 call check( nf90_inq_varid(ncid, 'dzl_dtau', ncv), thisroutine )
2014 nc2cor(1) = 1
2015 nc2cor(2) = 1
2016 nc2cnt(1) = imax + 1
2017 nc2cnt(2) = jmax + 1
2018 call check( nf90_put_var(ncid, ncv, dzl_dtau_conv, start=nc2cor, count=nc2cnt), &
2019  thisroutine )
2020 
2021 call check( nf90_inq_varid(ncid, 'dH_c_dtau', ncv), thisroutine )
2022 nc2cor(1) = 1
2023 nc2cor(2) = 1
2024 nc2cnt(1) = imax + 1
2025 nc2cnt(2) = jmax + 1
2026 call check( nf90_put_var(ncid, ncv, dh_c_dtau_conv, start=nc2cor, count=nc2cnt), &
2027  thisroutine )
2028 
2029 call check( nf90_inq_varid(ncid, 'dH_t_dtau', ncv), thisroutine )
2030 nc2cor(1) = 1
2031 nc2cor(2) = 1
2032 nc2cnt(1) = imax + 1
2033 nc2cnt(2) = jmax + 1
2034 call check( nf90_put_var(ncid, ncv, dh_t_dtau_conv, start=nc2cor, count=nc2cnt), &
2035  thisroutine )
2036 
2037 call check( nf90_inq_varid(ncid, 'dH_dtau', ncv), thisroutine )
2038 nc2cor(1) = 1
2039 nc2cor(2) = 1
2040 nc2cnt(1) = imax + 1
2041 nc2cnt(2) = jmax + 1
2042 call check( nf90_put_var(ncid, ncv, dh_dtau_conv, start=nc2cor, count=nc2cnt), &
2043  thisroutine )
2044 
2045 call check( nf90_inq_varid(ncid, 'vx_b_g', ncv), thisroutine )
2046 nc2cor(1) = 1
2047 nc2cor(2) = 1
2048 nc2cnt(1) = imax + 1
2049 nc2cnt(2) = jmax + 1
2050 call check( nf90_put_var(ncid, ncv, vx_b_g_conv, start=nc2cor, count=nc2cnt), &
2051  thisroutine )
2052 
2053 call check( nf90_inq_varid(ncid, 'vy_b_g', ncv), thisroutine )
2054 nc2cor(1) = 1
2055 nc2cor(2) = 1
2056 nc2cnt(1) = imax + 1
2057 nc2cnt(2) = jmax + 1
2058 call check( nf90_put_var(ncid, ncv, vy_b_g_conv, start=nc2cor, count=nc2cnt), &
2059  thisroutine )
2060 
2061 call check( nf90_inq_varid(ncid, 'vz_b', ncv), thisroutine )
2062 nc2cor(1) = 1
2063 nc2cor(2) = 1
2064 nc2cnt(1) = imax + 1
2065 nc2cnt(2) = jmax + 1
2066 call check( nf90_put_var(ncid, ncv, vz_b_conv, start=nc2cor, count=nc2cnt), &
2067  thisroutine )
2068 
2069 call check( nf90_inq_varid(ncid, 'vh_b', ncv), thisroutine )
2070 nc2cor(1) = 1
2071 nc2cor(2) = 1
2072 nc2cnt(1) = imax + 1
2073 nc2cnt(2) = jmax + 1
2074 call check( nf90_put_var(ncid, ncv, vh_b_conv, start=nc2cor, count=nc2cnt), &
2075  thisroutine )
2076 
2077 call check( nf90_inq_varid(ncid, 'vx_s_g', ncv), thisroutine )
2078 nc2cor(1) = 1
2079 nc2cor(2) = 1
2080 nc2cnt(1) = imax + 1
2081 nc2cnt(2) = jmax + 1
2082 call check( nf90_put_var(ncid, ncv, vx_s_g_conv, start=nc2cor, count=nc2cnt), &
2083  thisroutine )
2084 
2085 call check( nf90_inq_varid(ncid, 'vy_s_g', ncv), thisroutine )
2086 nc2cor(1) = 1
2087 nc2cor(2) = 1
2088 nc2cnt(1) = imax + 1
2089 nc2cnt(2) = jmax + 1
2090 call check( nf90_put_var(ncid, ncv, vy_s_g_conv, start=nc2cor, count=nc2cnt), &
2091  thisroutine )
2092 
2093 call check( nf90_inq_varid(ncid, 'vz_s', ncv), thisroutine )
2094 nc2cor(1) = 1
2095 nc2cor(2) = 1
2096 nc2cnt(1) = imax + 1
2097 nc2cnt(2) = jmax + 1
2098 call check( nf90_put_var(ncid, ncv, vz_s_conv, start=nc2cor, count=nc2cnt), &
2099  thisroutine )
2100 
2101 call check( nf90_inq_varid(ncid, 'vh_s', ncv), thisroutine )
2102 nc2cor(1) = 1
2103 nc2cor(2) = 1
2104 nc2cnt(1) = imax + 1
2105 nc2cnt(2) = jmax + 1
2106 call check( nf90_put_var(ncid, ncv, vh_s_conv, start=nc2cor, count=nc2cnt), &
2107  thisroutine )
2108 
2109 call check( nf90_inq_varid(ncid, 'temp_b', ncv), thisroutine )
2110 nc2cor(1) = 1
2111 nc2cor(2) = 1
2112 nc2cnt(1) = imax + 1
2113 nc2cnt(2) = jmax + 1
2114 call check( nf90_put_var(ncid, ncv, temp_b_conv, start=nc2cor, count=nc2cnt), &
2115  thisroutine )
2116 
2117 call check( nf90_inq_varid(ncid, 'temph_b', ncv), thisroutine )
2118 nc2cor(1) = 1
2119 nc2cor(2) = 1
2120 nc2cnt(1) = imax + 1
2121 nc2cnt(2) = jmax + 1
2122 call check( nf90_put_var(ncid, ncv, temph_b_conv, start=nc2cor, count=nc2cnt), &
2123  thisroutine )
2124 
2125 call check( nf90_inq_varid(ncid, 'p_b_w', ncv), thisroutine )
2126 nc2cor(1) = 1
2127 nc2cor(2) = 1
2128 nc2cnt(1) = imax + 1
2129 nc2cnt(2) = jmax + 1
2130 call check( nf90_put_var(ncid, ncv, p_b_w_conv, start=nc2cor, count=nc2cnt), &
2131  thisroutine )
2132 
2133 call check( nf90_inq_varid(ncid, 'H_w', ncv), thisroutine )
2134 nc2cor(1) = 1
2135 nc2cor(2) = 1
2136 nc2cnt(1) = imax + 1
2137 nc2cnt(2) = jmax + 1
2138 call check( nf90_put_var(ncid, ncv, h_w_conv, start=nc2cor, count=nc2cnt), &
2139  thisroutine )
2140 
2141 call check( nf90_inq_varid(ncid, 'q_gl_g', ncv), thisroutine )
2142 nc2cor(1) = 1
2143 nc2cor(2) = 1
2144 nc2cnt(1) = imax + 1
2145 nc2cnt(2) = jmax + 1
2146 call check( nf90_put_var(ncid, ncv, q_gl_g_conv, start=nc2cor, count=nc2cnt), &
2147  thisroutine )
2148 
2149 call check( nf90_inq_varid(ncid, 'V_tot', ncv), thisroutine )
2150 call check( nf90_put_var(ncid, ncv, v_tot_conv), thisroutine )
2151 
2152 call check( nf90_inq_varid(ncid, 'A_grounded', ncv), thisroutine )
2153 call check( nf90_put_var(ncid, ncv, a_grounded_conv), thisroutine )
2154 
2155 call check( nf90_inq_varid(ncid, 'A_floating', ncv), thisroutine )
2156 call check( nf90_put_var(ncid, ncv, a_floating_conv), thisroutine )
2157 
2158 !-------- Close NetCDF file --------
2159 
2160 call check( nf90_sync(ncid), thisroutine )
2161 
2162 call check( nf90_close(ncid), thisroutine )
2163 
2164 !-------- Increase file counter --------
2165 
2166 ndat = ndat+1
2167 
2168 if (flag_3d_output) then
2169  ndat3d = ndat
2170 else
2171  ndat2d = ndat
2172 end if
2173 
2174 end subroutine output_nc
2175 !