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