source: palm/trunk/SOURCE/netcdf_interface_mod.f90 @ 4281

Last change on this file since 4281 was 4227, checked in by gronemeier, 21 months ago

implement new palm_date_time_mod; replaced namelist parameters time_utc_init and day_of_year_init by origin_date_time

  • Property svn:keywords set to Id
File size: 317.3 KB
RevLine 
[1850]1!> @file netcdf_interface_mod.f90
[2000]2!------------------------------------------------------------------------------!
[2696]3! This file is part of the PALM model system.
[1036]4!
[2000]5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
[1036]9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
[3655]17! Copyright 1997-2019 Leibniz Universitaet Hannover
[2000]18!------------------------------------------------------------------------------!
[1036]19!
[1]20! Current revisions:
21! ------------------
[3464]22!
[3995]23!
[3464]24! Former revisions:
25! -----------------
26! $Id: netcdf_interface_mod.f90 4227 2019-09-10 18:04:34Z schwenkel $
[4227]27! Replace function date_time_string by call to get_date_time
28!
29! 4223 2019-09-10 09:20:47Z gronemeier
[4196]30! replaced rotation angle from input-netCDF file
31! by namelist parameter 'rotation_angle'
32!
33! 4182 2019-08-22 15:20:23Z scharf
[4182]34! Corrected "Former revisions" section
35!
36! 4127 2019-07-30 14:47:10Z suehring
[4127]37! -Introduce new vertical dimension for plant-canopy output.
38! -Temporarlily disable masked output for soil (merge from branch resler)
39!
40! 4069 2019-07-01 14:05:51Z Giersch
[4069]41! Masked output running index mid has been introduced as a local variable to
42! avoid runtime error (Loop variable has been modified) in time_integration
43!
44! 4046 2019-06-21 17:32:04Z knoop
[4046]45! removal of special treatment for usm_define_netcdf_grid call
46!
47! 4039 2019-06-18 10:32:41Z suehring
[4039]48! Rename subroutines in module for diagnostic quantities
[4046]49!
[4039]50! 4029 2019-06-14 14:04:35Z raasch
[4029]51! netcdf variable NF90_NOFILL is used as argument instead of "1" in calls to NF90_DEF_VAR_FILL
[4046]52!
[4029]53! 3995 2019-05-22 18:59:54Z suehring
[3995]54! output of turbulence intensity added
[4046]55!
[3995]56! 3994 2019-05-22 18:08:09Z suehring
[3966]57! remove origin time from time unit, compose origin_time_string within
58! subroutine netcdf_create_global_atts
[4046]59!
[3966]60! 3954 2019-05-06 12:49:42Z gronemeier
[3954]61! bugfix: corrected format for date_time_string
[4046]62!
[3954]63! 3953 2019-05-06 12:11:55Z gronemeier
[3953]64! bugfix: set origin_time and starting point of time coordinate according to
65!         day_of_year_init and time_utc_init
[4046]66!
[3953]67! 3942 2019-04-30 13:08:30Z kanani
[3942]68! Add specifier to netcdf_handle_error to simplify identification of attribute
69! causing the error
[4046]70!
[3942]71! 3766 2019-02-26 16:23:41Z raasch
[3766]72! bugfix in im_define_netcdf_grid argument list
[4046]73!
[3766]74! 3745 2019-02-15 18:57:56Z suehring
[3745]75! Add indoor model
[4046]76!
[3745]77! 3744 2019-02-15 18:38:58Z suehring
[3729]78! Bugfix: - initialize return values to ensure they are set before returning
79!           (routine define_geo_coordinates)
80!         - change order of dimensions for some variables
[4046]81!
[3729]82! 3727 2019-02-08 14:52:10Z gronemeier
[3727]83! make several routines publicly available
[4046]84!
[3727]85! 3701 2019-01-26 18:57:21Z knoop
[3665]86! Statement added to prevent compiler warning about unused variable
[4046]87!
[3665]88! 3655 2019-01-07 16:51:22Z knoop
[4046]89! Move the control parameter "salsa" from salsa_mod to control_parameters
[3589]90! (M. Kurppa)
[4046]91!
[4182]92! Revision 1.1  2005/05/18 15:37:16  raasch
93! Initial revision
[4046]94!
[4182]95!
[1]96! Description:
97! ------------
[1682]98!> In case of extend = .FALSE.:
99!> Define all necessary dimensions, axes and variables for the different
100!> netCDF datasets. This subroutine is called from check_open after a new
101!> dataset is created. It leaves the open netCDF files ready to write.
102!>
103!> In case of extend = .TRUE.:
104!> Find out if dimensions and variables of an existing file match the values
105!> of the actual run. If so, get all necessary information (ids, etc.) from
106!> this file.
107!>
108!> Parameter av can assume values 0 (non-averaged data) and 1 (time averaged
109!> data)
[1745]110!>
111!> @todo calculation of output time levels for parallel NetCDF still does not
112!>       cover every exception (change of dt_do, end_time in restart)
[1976]113!> @todo timeseries and profile output still needs to be rewritten to allow
114!>       modularization
[3459]115!> @todo output 2d UTM coordinates without global arrays
116!> @todo output longitude/latitude also with non-parallel output (3d and xy)
[1]117!------------------------------------------------------------------------------!
[1783]118 MODULE netcdf_interface
119
[2232]120    USE control_parameters,                                                    &
[3448]121        ONLY:  biometeorology, fl_max,                                         &
122               max_masks, multi_agent_system_end,                              &
[4196]123               multi_agent_system_start,                                       &
124               rotation_angle,                                                 &
125               var_fl_max, varnamelength
[1783]126    USE kinds
127#if defined( __netcdf )
128    USE NETCDF
[1682]129#endif
[3235]130    USE mas_global_attributes,                                                 &
131        ONLY:  dim_size_agtnum
[1783]132
[3459]133    USE netcdf_data_input_mod,                                                 &
[3485]134        ONLY: coord_ref_sys, init_model
[3459]135
[1783]136    PRIVATE
137
[3159]138    CHARACTER (LEN=16), DIMENSION(13) ::  agt_var_names =                      &
[3235]139          (/ 'ag_id           ', 'ag_x            ', 'ag_y            ',       &
140             'ag_wind         ', 'ag_temp         ', 'ag_group        ',       &
[3448]141             'PM10            ', 'PM25            ', 'ag_iPT          ',       &
[3165]142             'ag_uv           ', 'not_used        ', 'not_used        ',       &
[3159]143             'not_used        ' /)
144
145    CHARACTER (LEN=16), DIMENSION(13) ::  agt_var_units = &
[3235]146          (/ 'dim_less        ', 'meters          ', 'meters          ',       &
147             'm/s             ', 'K               ', 'dim_less        ',       &
[3159]148             'tbd             ', 'tbd             ', 'tbd             ',       &
[3165]149             'tbd             ', 'not_used        ', 'not_used        ',       &
[3159]150             'not_used        ' /)
151
[1783]152    INTEGER(iwp), PARAMETER ::  dopr_norm_num = 7, dopts_num = 29, dots_max = 100
153
[3421]154    CHARACTER (LEN=7), DIMENSION(dopr_norm_num) ::  dopr_norm_names =          &
155         (/ 'wtheta0', 'ws2    ', 'tsw2   ', 'ws3    ', 'ws2tsw ', 'wstsw2 ',  &
156            'z_i    ' /)
[1783]157
[3421]158    CHARACTER (LEN=7), DIMENSION(dopr_norm_num) ::  dopr_norm_longnames =      &
159         (/ 'wtheta0', 'w*2    ', 't*w2   ', 'w*3    ', 'w*2t*w ', 'w*t*w2 ',  &
160            'z_i    ' /)
[1783]161
162    CHARACTER (LEN=7), DIMENSION(dopts_num) :: dopts_label =                   &
163          (/ 'tnpt   ', 'x_     ', 'y_     ', 'z_     ', 'z_abs  ', 'u      ', &
164             'v      ', 'w      ', 'u"     ', 'v"     ', 'w"     ', 'npt_up ', &
165             'w_up   ', 'w_down ', 'radius ', 'r_min  ', 'r_max  ', 'npt_max', &
166             'npt_min', 'x*2    ', 'y*2    ', 'z*2    ', 'u*2    ', 'v*2    ', &
167             'w*2    ', 'u"2    ', 'v"2    ', 'w"2    ', 'npt*2  ' /)
168
169    CHARACTER (LEN=7), DIMENSION(dopts_num) :: dopts_unit =                    &
170          (/ 'number ', 'm      ', 'm      ', 'm      ', 'm      ', 'm/s    ', &
171             'm/s    ', 'm/s    ', 'm/s    ', 'm/s    ', 'm/s    ', 'number ', &
172             'm/s    ', 'm/s    ', 'm      ', 'm      ', 'm      ', 'number ', &
173             'number ', 'm2     ', 'm2     ', 'm2     ', 'm2/s2  ', 'm2/s2  ', &
174             'm2/s2  ', 'm2/s2  ', 'm2/s2  ', 'm2/s2  ', 'number2' /)
175
[2270]176    INTEGER(iwp) ::  dots_num  = 25  !< number of timeseries defined by default
[1960]177    INTEGER(iwp) ::  dots_soil = 26  !< starting index for soil-timeseries
[2270]178    INTEGER(iwp) ::  dots_rad  = 32  !< starting index for radiation-timeseries
[1783]179
180    CHARACTER (LEN=13), DIMENSION(dots_max) :: dots_label =                    &
181          (/ 'E            ', 'E*           ', 'dt           ',                &
[3421]182             'us*          ', 'th*          ', 'umax         ',                &
[1783]183             'vmax         ', 'wmax         ', 'div_new      ',                &
[3421]184             'div_old      ', 'zi_wtheta    ', 'zi_theta     ',                &
185             'w*           ', 'w"theta"0    ', 'w"theta"     ',                &
186             'wtheta       ', 'theta(0)     ', 'theta(z_mo)  ',                &
[1783]187             'w"u"0        ', 'w"v"0        ', 'w"q"0        ',                &
[1960]188             'ol           ', 'q*           ', 'w"s"         ',                &
[2270]189             's*           ', 'ghf          ', 'qsws_liq     ',                &
190             'qsws_soil    ', 'qsws_veg     ', 'r_a          ',                &
191             'r_s          ',                                                  &
[1960]192             'rad_net      ', 'rad_lw_in    ', 'rad_lw_out   ',                &
193             'rad_sw_in    ', 'rad_sw_out   ', 'rrtm_aldif   ',                &
[4046]194             'rrtm_aldir   ', 'rrtm_asdif   ', 'rrtm_asdir   ',                &
[2270]195             ( 'unknown      ', i9 = 1, dots_max-40 ) /)
[1783]196
197    CHARACTER (LEN=13), DIMENSION(dots_max) :: dots_unit =                     &
198          (/ 'm2/s2        ', 'm2/s2        ', 's            ',                &
199             'm/s          ', 'K            ', 'm/s          ',                &
200             'm/s          ', 'm/s          ', 's-1          ',                &
201             's-1          ', 'm            ', 'm            ',                &
202             'm/s          ', 'K m/s        ', 'K m/s        ',                &
203             'K m/s        ', 'K            ', 'K            ',                &
204             'm2/s2        ', 'm2/s2        ', 'kg m/s       ',                &
[1960]205             'm            ', 'kg/kg        ', 'kg m/(kg s)  ',                &
[1972]206             'kg/kg        ', 'W/m2         ', 'W/m2         ',                &
[2270]207             'W/m2         ', 'W/m2         ', 's/m          ',                &
208             's/m          ',                                                  &
[1783]209             'W/m2         ', 'W/m2         ', 'W/m2         ',                &
[1960]210             'W/m2         ', 'W/m2         ', '             ',                &
[1783]211             '             ', '             ', '             ',                &
[2270]212             ( 'unknown      ', i9 = 1, dots_max-40 ) /)
[1783]213
[2037]214    CHARACTER (LEN=16) :: heatflux_output_unit     !< unit for heatflux output
215    CHARACTER (LEN=16) :: waterflux_output_unit    !< unit for waterflux output
216    CHARACTER (LEN=16) :: momentumflux_output_unit !< unit for momentumflux output
217
[1783]218    CHARACTER (LEN=9), DIMENSION(300) ::  dopr_unit = 'unknown'
219
[2007]220    CHARACTER (LEN=7), DIMENSION(0:1,500) ::  do2d_unit, do3d_unit
[1783]221
[3241]222!    CHARACTER (LEN=16), DIMENSION(25) ::  prt_var_names = &
223!          (/ 'pt_age          ', 'pt_dvrp_size    ', 'pt_origin_x     ', &
224!             'pt_origin_y     ', 'pt_origin_z     ', 'pt_radius       ', &
225!             'pt_speed_x      ', 'pt_speed_y      ', 'pt_speed_z      ', &
226!             'pt_weight_factor', 'pt_x            ', 'pt_y            ', &
227!             'pt_z            ', 'pt_color        ', 'pt_group        ', &
228!             'pt_tailpoints   ', 'pt_tail_id      ', 'pt_density_ratio', &
229!             'pt_exp_arg      ', 'pt_exp_term     ', 'not_used        ', &
230!             'not_used        ', 'not_used        ', 'not_used        ', &
231!             'not_used        ' /)
[1783]232
[3241]233!    CHARACTER (LEN=16), DIMENSION(25) ::  prt_var_units = &
234!          (/ 'seconds         ', 'meters          ', 'meters          ', &
235!             'meters          ', 'meters          ', 'meters          ', &
236!             'm/s             ', 'm/s             ', 'm/s             ', &
237!             'factor          ', 'meters          ', 'meters          ', &
238!             'meters          ', 'none            ', 'none            ', &
239!             'none            ', 'none            ', 'ratio           ', &
240!             'none            ', 'none            ', 'not_used        ', &
241!             'not_used        ', 'not_used        ', 'not_used        ', &
242!             'not_used        ' /)
[1783]243
244    CHARACTER(LEN=20), DIMENSION(11) ::  netcdf_precision = ' '
245    CHARACTER(LEN=40) ::  netcdf_data_format_string
246
[3241]247    INTEGER(iwp) ::  id_dim_agtnum, id_dim_time_agt,                           &
248                     id_dim_time_fl, id_dim_time_pr,                           &
[3159]249                     id_dim_time_pts, id_dim_time_sp, id_dim_time_ts,          &
250                     id_dim_x_sp, id_dim_y_sp, id_dim_zu_sp, id_dim_zw_sp,     &
251                     id_set_agt, id_set_fl, id_set_pr, id_set_prt, id_set_pts, &
252                     id_set_sp, id_set_ts, id_var_agtnum, id_var_time_agt,     &
[3241]253                     id_var_time_fl, id_var_rnoa_agt, id_var_time_pr,          &
[3159]254                     id_var_time_pts, id_var_time_sp, id_var_time_ts,          &
255                     id_var_x_sp, id_var_y_sp, id_var_zu_sp, id_var_zw_sp,     &
256                     nc_stat
[1783]257
[1957]258
[1783]259    INTEGER(iwp), DIMENSION(0:1) ::  id_dim_time_xy, id_dim_time_xz, &
260                    id_dim_time_yz, id_dim_time_3d, id_dim_x_xy, id_dim_xu_xy, &
261                    id_dim_x_xz, id_dim_xu_xz, id_dim_x_yz, id_dim_xu_yz, &
262                    id_dim_x_3d, id_dim_xu_3d, id_dim_y_xy, id_dim_yv_xy, &
263                    id_dim_y_xz, id_dim_yv_xz, id_dim_y_yz, id_dim_yv_yz, &
264                    id_dim_y_3d, id_dim_yv_3d, id_dim_zs_xy, id_dim_zs_xz, &
[4127]265                    id_dim_zs_yz, id_dim_zs_3d, id_dim_zpc_3d, &
266                    id_dim_zu_xy, id_dim_zu1_xy, &
[1783]267                    id_dim_zu_xz, id_dim_zu_yz, id_dim_zu_3d, id_dim_zw_xy, &
268                    id_dim_zw_xz, id_dim_zw_yz, id_dim_zw_3d, id_set_xy, &
269                    id_set_xz, id_set_yz, id_set_3d, id_var_ind_x_yz, &
270                    id_var_ind_y_xz, id_var_ind_z_xy, id_var_time_xy, &
271                    id_var_time_xz, id_var_time_yz, id_var_time_3d, id_var_x_xy, &
272                    id_var_xu_xy, id_var_x_xz, id_var_xu_xz, id_var_x_yz, &
273                    id_var_xu_yz, id_var_x_3d, id_var_xu_3d, id_var_y_xy, &
274                    id_var_yv_xy, id_var_y_xz, id_var_yv_xz, id_var_y_yz, &
275                    id_var_yv_yz, id_var_y_3d, id_var_yv_3d, id_var_zs_xy, &
[4127]276                    id_var_zs_xz, id_var_zs_yz, id_var_zs_3d, id_var_zpc_3d, &
277                    id_var_zusi_xy, id_var_zusi_3d, id_var_zu_xy, id_var_zu1_xy, id_var_zu_xz, &
[1783]278                    id_var_zu_yz, id_var_zu_3d, id_var_zwwi_xy, id_var_zwwi_3d, &
279                    id_var_zw_xy, id_var_zw_xz, id_var_zw_yz, id_var_zw_3d
280
[3729]281    INTEGER(iwp), DIMENSION(0:2,0:1) ::  id_var_eutm_3d, id_var_nutm_3d, &
[3421]282                                         id_var_eutm_xy, id_var_nutm_xy, &
283                                         id_var_eutm_xz, id_var_nutm_xz, &
[3459]284                                         id_var_eutm_yz, id_var_nutm_yz
[3421]285
[3729]286    INTEGER(iwp), DIMENSION(0:2,0:1) ::  id_var_lat_3d, id_var_lon_3d, &
[3459]287                                         id_var_lat_xy, id_var_lon_xy, &
288                                         id_var_lat_xz, id_var_lon_xz, &
289                                         id_var_lat_yz, id_var_lon_yz
290
[1783]291    INTEGER ::  netcdf_data_format = 2  !< NetCDF3 64bit offset format
292    INTEGER ::  netcdf_deflate = 0      !< NetCDF compression, default: no
293                                        !< compression
294
[1957]295    INTEGER(iwp)                 ::  dofl_time_count
[1783]296    INTEGER(iwp), DIMENSION(10)  ::  id_var_dospx, id_var_dospy
[3159]297    INTEGER(iwp), DIMENSION(20)  ::  id_var_agt
[3241]298!    INTEGER(iwp), DIMENSION(20)  ::  id_var_prt
[1783]299    INTEGER(iwp), DIMENSION(11)  ::  nc_precision
300    INTEGER(iwp), DIMENSION(dopr_norm_num) ::  id_var_norm_dopr
[4046]301
[1957]302    INTEGER(iwp), DIMENSION(fl_max) ::  id_dim_x_fl, id_dim_y_fl, id_dim_z_fl
303    INTEGER(iwp), DIMENSION(fl_max) ::  id_var_x_fl, id_var_y_fl, id_var_z_fl
[4046]304
[1957]305    CHARACTER (LEN=20), DIMENSION(fl_max*var_fl_max) :: dofl_label
[4046]306    CHARACTER (LEN=20), DIMENSION(fl_max*var_fl_max) :: dofl_unit
[1957]307    CHARACTER (LEN=20), DIMENSION(fl_max) :: dofl_dim_label_x
308    CHARACTER (LEN=20), DIMENSION(fl_max) :: dofl_dim_label_y
309    CHARACTER (LEN=20), DIMENSION(fl_max) :: dofl_dim_label_z
[1783]310
[4046]311    INTEGER(iwp), DIMENSION(fl_max*var_fl_max) :: id_var_dofl
[1957]312
[1783]313    INTEGER(iwp), DIMENSION(dopts_num,0:10) ::  id_var_dopts
[2007]314    INTEGER(iwp), DIMENSION(0:1,500)        ::  id_var_do2d, id_var_do3d
[2039]315    INTEGER(iwp), DIMENSION(100,0:99)       ::  id_dim_z_pr, id_var_dopr, &
[1783]316                                                id_var_z_pr
[2039]317    INTEGER(iwp), DIMENSION(dots_max,0:99)  ::  id_var_dots
[1783]318
319!
320!-- Masked output
321    CHARACTER (LEN=7), DIMENSION(max_masks,0:1,100) ::  domask_unit
322
323    LOGICAL ::  output_for_t0 = .FALSE.
324
325    INTEGER(iwp), DIMENSION(1:max_masks,0:1) ::  id_dim_time_mask, id_dim_x_mask, &
326                   id_dim_xu_mask, id_dim_y_mask, id_dim_yv_mask, id_dim_zs_mask, &
327                   id_dim_zu_mask, id_dim_zw_mask, &
328                   id_set_mask, &
329                   id_var_time_mask, id_var_x_mask, id_var_xu_mask, &
330                   id_var_y_mask, id_var_yv_mask, id_var_zs_mask, &
331                   id_var_zu_mask, id_var_zw_mask, &
332                   id_var_zusi_mask, id_var_zwwi_mask
[3448]333
[3729]334    INTEGER(iwp), DIMENSION(0:2,1:max_masks,0:1) ::  id_var_eutm_mask, &
[3421]335                                                     id_var_nutm_mask
[1783]336
[3729]337    INTEGER(iwp), DIMENSION(0:2,1:max_masks,0:1) ::  id_var_lat_mask, &
[3459]338                                                     id_var_lon_mask
339
[1783]340    INTEGER(iwp), DIMENSION(1:max_masks,0:1,100) ::  id_var_domask
341
[3529]342    REAL(wp) ::  fill_value = -9999.0_wp    !< value for the _FillValue attribute
[1783]343
[2696]344
[1957]345    PUBLIC  dofl_dim_label_x, dofl_dim_label_y, dofl_dim_label_z, dofl_label,  &
346            dofl_time_count, dofl_unit, domask_unit, dopr_unit, dopts_num,     &
347            dots_label, dots_max, dots_num, dots_rad, dots_soil, dots_unit,    &
[3159]348            do2d_unit, do3d_unit, fill_value, id_set_agt, id_set_fl,           &
349            id_set_mask, id_set_pr, id_set_prt, id_set_pts, id_set_sp,         &
350            id_set_ts, id_set_xy, id_set_xz, id_set_yz, id_set_3d, id_var_agt, &
351            id_var_domask, id_var_dofl, id_var_dopr, id_var_dopts,             &
352            id_var_dospx, id_var_dospy, id_var_dots, id_var_do2d, id_var_do3d, &
353            id_var_norm_dopr, id_var_time_agt, id_var_time_fl,                 &
354            id_var_time_mask, id_var_time_pr, id_var_rnoa_agt, id_var_time_pts,&
355            id_var_time_sp, id_var_time_ts,                                    &
[1957]356            id_var_time_xy, id_var_time_xz, id_var_time_yz, id_var_time_3d,    &
357            id_var_x_fl, id_var_y_fl, id_var_z_fl,  nc_stat,                   &
358            netcdf_data_format, netcdf_data_format_string, netcdf_deflate,     &
[2037]359            netcdf_precision, output_for_t0, heatflux_output_unit,             &
360            waterflux_output_unit, momentumflux_output_unit
361
[1783]362    SAVE
363
364    INTERFACE netcdf_create_dim
365       MODULE PROCEDURE netcdf_create_dim
366    END INTERFACE netcdf_create_dim
367
368    INTERFACE netcdf_create_file
369       MODULE PROCEDURE netcdf_create_file
370    END INTERFACE netcdf_create_file
371
[3727]372    INTERFACE netcdf_create_global_atts
373       MODULE PROCEDURE netcdf_create_global_atts
374    END INTERFACE netcdf_create_global_atts
375
[1783]376    INTERFACE netcdf_create_var
377       MODULE PROCEDURE netcdf_create_var
378    END INTERFACE netcdf_create_var
379
[3529]380    INTERFACE netcdf_create_att
381       MODULE PROCEDURE netcdf_create_att_string
382    END INTERFACE netcdf_create_att
383
[1783]384    INTERFACE netcdf_define_header
385       MODULE PROCEDURE netcdf_define_header
386    END INTERFACE netcdf_define_header
387
388    INTERFACE netcdf_handle_error
389       MODULE PROCEDURE netcdf_handle_error
390    END INTERFACE netcdf_handle_error
391
392    INTERFACE netcdf_open_write_file
393       MODULE PROCEDURE netcdf_open_write_file
394    END INTERFACE netcdf_open_write_file
395
[3727]396    PUBLIC netcdf_create_att, netcdf_create_dim, netcdf_create_file,           &
397           netcdf_create_global_atts, netcdf_create_var, netcdf_define_header, &
[2696]398           netcdf_handle_error, netcdf_open_write_file
[1783]399
400 CONTAINS
401
402 SUBROUTINE netcdf_define_header( callmode, extend, av )
[4046]403
[1]404#if defined( __netcdf )
405
[1691]406    USE arrays_3d,                                                             &
[1320]407        ONLY:  zu, zw
408
[3448]409    USE biometeorology_mod,                                                    &
[3525]410        ONLY:  bio_define_netcdf_grid
[3448]411
[2696]412    USE chemistry_model_mod,                                                   &
[4046]413        ONLY:  chem_define_netcdf_grid
[2696]414
[3274]415    USE basic_constants_and_equations_mod,                                     &
[1320]416        ONLY:  pi
417
[1691]418    USE control_parameters,                                                    &
[3159]419        ONLY:  agent_time_unlimited, air_chemistry, averaging_interval,        &
420               averaging_interval_pr, data_output_pr, domask, dopr_n,          &
[1691]421               dopr_time_count, dopts_time_count, dots_time_count,             &
[2769]422               do2d, do2d_at_begin, do2d_xz_time_count, do3d, do3d_at_begin,   &
[1745]423               do2d_yz_time_count, dt_data_output_av, dt_do2d_xy, dt_do2d_xz,  &
[3159]424               dt_do2d_yz, dt_do3d, dt_write_agent_data, mask_size,            &
[3744]425               do2d_xy_time_count, do3d_time_count, domask_time_count,         &
426               end_time, indoor_model, land_surface,                           &
[2696]427               mask_size_l, mask_i, mask_i_global, mask_j, mask_j_global,      &
[3435]428               mask_k_global, mask_surface,                                    &
[4069]429               message_string, ntdim_2d_xy, ntdim_2d_xz,                       &
[3294]430               ntdim_2d_yz, ntdim_3d, nz_do3d, ocean_mode, plant_canopy,       &
[3582]431               run_description_header, salsa, section, simulated_time,         &
[1745]432               simulated_time_at_begin, skip_time_data_output_av,              &
433               skip_time_do2d_xy, skip_time_do2d_xz, skip_time_do2d_yz,        &
[2011]434               skip_time_do3d, topography, num_leg, num_var_fl,                &
[3569]435               urban_surface
[1320]436
[4046]437    USE diagnostic_output_quantities_mod,                                      &
438        ONLY:  doq_define_netcdf_grid
[4039]439
[1691]440    USE grid_variables,                                                        &
[1320]441        ONLY:  dx, dy, zu_s_inner, zw_w_inner
442
[2817]443    USE gust_mod,                                                              &
444        ONLY: gust_define_netcdf_grid, gust_module_enabled
445
[1691]446    USE indices,                                                               &
[2232]447        ONLY:  nx, nxl, nxr, ny, nys, nyn, nz ,nzb, nzt
[1320]448
449    USE kinds
450
[3744]451    USE indoor_model_mod,                                                      &
452        ONLY: im_define_netcdf_grid
[4046]453
[1551]454    USE land_surface_model_mod,                                                &
[2232]455        ONLY: lsm_define_netcdf_grid, nzb_soil, nzt_soil, nzs, zs
[1551]456
[3294]457    USE ocean_mod,                                                             &
458        ONLY:  ocean_define_netcdf_grid
459
[1]460    USE pegrid
461
[1691]462    USE particle_attributes,                                                   &
[2265]463        ONLY:  number_of_particle_groups
[1]464
[2209]465    USE plant_canopy_model_mod,                                                &
[4127]466        ONLY:  pch_index, pcm_define_netcdf_grid
[2209]467
[1691]468    USE profil_parameter,                                                      &
469        ONLY:  crmax, cross_profiles, dopr_index, profile_columns, profile_rows
[1320]470
[1976]471    USE radiation_model_mod,                                                   &
[4046]472        ONLY: radiation, radiation_define_netcdf_grid
473
[3467]474    USE salsa_mod,                                                             &
[4046]475        ONLY:  salsa_define_netcdf_grid
[1976]476
[1833]477    USE spectra_mod,                                                           &
478        ONLY:  averaging_interval_sp, comp_spectra_level, data_output_sp, dosp_time_count, spectra_direction
[1320]479
[1691]480    USE statistics,                                                            &
[1320]481        ONLY:  hom, statistic_regions
482
[2696]483    USE turbulence_closure_mod,                                                &
484        ONLY:  tcm_define_netcdf_grid
485
[2007]486    USE urban_surface_mod,                                                     &
[2011]487        ONLY:  usm_define_netcdf_grid
[1320]488
[3701]489    USE user,                                                                  &
490        ONLY:  user_module_enabled, user_define_netcdf_grid
[2007]491
[2696]492
[3701]493
[1]494    IMPLICIT NONE
495
[2039]496    CHARACTER (LEN=3)              ::  suffix                !<
[1682]497    CHARACTER (LEN=2), INTENT (IN) ::  callmode              !<
498    CHARACTER (LEN=4)              ::  grid_x                !<
499    CHARACTER (LEN=4)              ::  grid_y                !<
500    CHARACTER (LEN=4)              ::  grid_z                !<
501    CHARACTER (LEN=6)              ::  mode                  !<
502    CHARACTER (LEN=10)             ::  precision             !<
503    CHARACTER (LEN=10)             ::  var                   !<
[2109]504    CHARACTER (LEN=20)             ::  netcdf_var_name       !<
[2011]505    CHARACTER (LEN=varnamelength)  ::  trimvar               !< TRIM of output-variable string
[1682]506    CHARACTER (LEN=80)             ::  time_average_text     !<
[2007]507    CHARACTER (LEN=4000)           ::  char_cross_profiles   !<
508    CHARACTER (LEN=4000)           ::  var_list              !<
509    CHARACTER (LEN=4000)           ::  var_list_old          !<
[1]510
[1682]511    CHARACTER (LEN=100), DIMENSION(1:crmax) ::  cross_profiles_adj   !<
512    CHARACTER (LEN=100), DIMENSION(1:crmax) ::  cross_profiles_char  !<
[1]513
[1682]514    INTEGER(iwp) ::  av                                      !<
515    INTEGER(iwp) ::  cross_profiles_count                    !<
516    INTEGER(iwp) ::  cross_profiles_maxi                     !<
517    INTEGER(iwp) ::  delim                                   !<
518    INTEGER(iwp) ::  delim_old                               !<
519    INTEGER(iwp) ::  file_id                                 !<
520    INTEGER(iwp) ::  i                                       !<
521    INTEGER(iwp) ::  id_last                                 !<
522    INTEGER(iwp) ::  id_x                                    !<
523    INTEGER(iwp) ::  id_y                                    !<
524    INTEGER(iwp) ::  id_z                                    !<
525    INTEGER(iwp) ::  j                                       !<
526    INTEGER(iwp) ::  k                                       !<
527    INTEGER(iwp) ::  kk                                      !<
[4069]528    INTEGER(iwp) ::  mid                                     !< masked output running index
[1682]529    INTEGER(iwp) ::  ns                                      !<
530    INTEGER(iwp) ::  ns_do                                   !< actual value of ns for soil model data
531    INTEGER(iwp) ::  ns_old                                  !<
[1745]532    INTEGER(iwp) ::  ntime_count                             !< number of time levels found in file
[1682]533    INTEGER(iwp) ::  nz_old                                  !<
[1957]534    INTEGER(iwp) ::  l                                       !<
[951]535
[1682]536    INTEGER(iwp), SAVE ::  oldmode                           !<
[1308]537
[1682]538    INTEGER(iwp), DIMENSION(1) ::  id_dim_time_old           !<
539    INTEGER(iwp), DIMENSION(1) ::  id_dim_x_yz_old           !<
540    INTEGER(iwp), DIMENSION(1) ::  id_dim_y_xz_old           !<
541    INTEGER(iwp), DIMENSION(1) ::  id_dim_zu_sp_old          !<
542    INTEGER(iwp), DIMENSION(1) ::  id_dim_zu_xy_old          !<
543    INTEGER(iwp), DIMENSION(1) ::  id_dim_zu_3d_old          !<
544    INTEGER(iwp), DIMENSION(1) ::  id_dim_zu_mask_old        !<
[1]545
546
[1682]547    INTEGER(iwp), DIMENSION(1:crmax) ::  cross_profiles_numb !<
[951]548
[1682]549    LOGICAL ::  found                                        !<
[1]550
[1682]551    LOGICAL, INTENT (INOUT) ::  extend                       !<
[1]552
[1682]553    LOGICAL, SAVE ::  init_netcdf = .FALSE.                  !<
[1]554
[4196]555    REAL(wp) ::  cos_rot_angle                               !< cosine of rotation_angle
[3459]556    REAL(wp) ::  eutm                                        !< easting (UTM)
557    REAL(wp) ::  nutm                                        !< northing (UTM)
[3421]558    REAL(wp) ::  shift_x                                     !< shift of x coordinate
559    REAL(wp) ::  shift_y                                     !< shift of y coordinate
[4196]560    REAL(wp) ::  sin_rot_angle                               !< sine of rotation_angle
[3421]561
[1745]562    REAL(wp), DIMENSION(1) ::  last_time_coordinate          !< last time value in file
[3459]563    REAL(wp), DIMENSION(8) ::  crs_list                      !< list of coord_ref_sys values
[1]564
[1682]565    REAL(wp), DIMENSION(:), ALLOCATABLE   ::  netcdf_data    !<
566    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  netcdf_data_2d !<
[3459]567    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lat            !< latitude
568    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lon            !< longitude
[1320]569
[2007]570
[1]571!
[1783]572!-- Initializing actions
[1]573    IF ( .NOT. init_netcdf )  THEN
574!
[1031]575!--    Check and set accuracy for netCDF output. First set default value
[1]576       nc_precision = NF90_REAL4
577
578       i = 1
579       DO  WHILE ( netcdf_precision(i) /= ' ' )
580          j = INDEX( netcdf_precision(i), '_' )
581          IF ( j == 0 )  THEN
[274]582             WRITE ( message_string, * ) 'netcdf_precision must contain a ', &
583                                         '"_"netcdf_precision(', i, ')="',   &
[257]584                                         TRIM( netcdf_precision(i) ),'"'
[1783]585             CALL message( 'netcdf_define_header', 'PA0241', 2, 2, 0, 6, 0 )
[1]586          ENDIF
587
588          var       = netcdf_precision(i)(1:j-1)
589          precision = netcdf_precision(i)(j+1:)
590
591          IF ( precision == 'NF90_REAL4' )  THEN
592             j = NF90_REAL4
593          ELSEIF ( precision == 'NF90_REAL8' )  THEN
594             j = NF90_REAL8
595          ELSE
[257]596             WRITE ( message_string, * ) 'illegal netcdf precision: ',  &
597                                         'netcdf_precision(', i, ')="', &
598                                         TRIM( netcdf_precision(i) ),'"'
[1783]599             CALL message( 'netcdf_define_header', 'PA0242', 1, 2, 0, 6, 0 )
[1]600          ENDIF
601
602          SELECT CASE ( var )
603             CASE ( 'xy' )
604                nc_precision(1) = j
605             CASE ( 'xz' )
606                nc_precision(2) = j
607             CASE ( 'yz' )
608                nc_precision(3) = j
609             CASE ( '2d' )
610                nc_precision(1:3) = j
611             CASE ( '3d' )
612                nc_precision(4) = j
613             CASE ( 'pr' )
614                nc_precision(5) = j
615             CASE ( 'ts' )
616                nc_precision(6) = j
617             CASE ( 'sp' )
618                nc_precision(7) = j
619             CASE ( 'prt' )
620                nc_precision(8) = j
[410]621             CASE ( 'masks' )
622                nc_precision(11) = j
[1957]623             CASE ( 'fl' )
624                nc_precision(9) = j
[1]625             CASE ( 'all' )
626                nc_precision    = j
627
628             CASE DEFAULT
[2932]629                WRITE ( message_string, * ) 'unknown variable in ' //          &
[4046]630                                  'initialization_parameters ',                &
[2932]631                                  'assignment: netcdf_precision(', i, ')="',   &
[257]632                                            TRIM( netcdf_precision(i) ),'"'
[1783]633                CALL message( 'netcdf_define_header', 'PA0243', 1, 2, 0, 6, 0 )
[1]634
635          END SELECT
636
637          i = i + 1
[410]638          IF ( i > 50 )  EXIT
[1]639       ENDDO
640
[1783]641!
642!--    Check for allowed parameter range
643       IF ( netcdf_deflate < 0  .OR.  netcdf_deflate > 9 )  THEN
644          WRITE ( message_string, '(A,I3,A)' ) 'netcdf_deflate out of ' //     &
[3046]645                                      'range & given value: ', netcdf_deflate, &
[1783]646                                      ', allowed range: 0-9'
647          CALL message( 'netcdf_define_header', 'PA0355', 2, 2, 0, 6, 0 )
648       ENDIF
649!
650!--    Data compression does not work with parallel NetCDF/HDF5
651       IF ( netcdf_deflate > 0  .AND.  netcdf_data_format /= 3 )  THEN
652          message_string = 'netcdf_deflate reset to 0'
653          CALL message( 'netcdf_define_header', 'PA0356', 0, 1, 0, 6, 0 )
654
655          netcdf_deflate = 0
656       ENDIF
657
[1]658       init_netcdf = .TRUE.
659
660    ENDIF
661!
[3459]662!-- Convert coord_ref_sys into vector (used for lat/lon calculation)
663    crs_list = (/ coord_ref_sys%semi_major_axis,                  &
664                  coord_ref_sys%inverse_flattening,               &
665                  coord_ref_sys%longitude_of_prime_meridian,      &
666                  coord_ref_sys%longitude_of_central_meridian,    &
667                  coord_ref_sys%scale_factor_at_central_meridian, &
668                  coord_ref_sys%latitude_of_projection_origin,    &
669                  coord_ref_sys%false_easting,                    &
670                  coord_ref_sys%false_northing /)
671
672!
[1]673!-- Determine the mode to be processed
674    IF ( extend )  THEN
675       mode = callmode // '_ext'
676    ELSE
677       mode = callmode // '_new'
678    ENDIF
679
680!
[4046]681!-- Select the mode to be processed. Possibilities are 3d, ma (mask), xy, xz,
682!-- yz, pr (profiles), ps (particle timeseries), fl (flight data), ts
[3045]683!-- (timeseries) or sp (spectra)
[1]684    SELECT CASE ( mode )
685
[410]686       CASE ( 'ma_new' )
687
688!
[4046]689!--       decompose actual parameter file_id (=formal parameter av) into
[410]690!--       mid and av
691          file_id = av
[564]692          IF ( file_id <= 200+max_masks )  THEN
693             mid = file_id - 200
[410]694             av = 0
695          ELSE
[564]696             mid = file_id - (200+max_masks)
[410]697             av = 1
698          ENDIF
699
700!
701!--       Define some global attributes of the dataset
702          IF ( av == 0 )  THEN
[3529]703             CALL netcdf_create_global_atts( id_set_mask(mid,av), 'podsmasked', TRIM( run_description_header ), 464 )
[410]704             time_average_text = ' '
705          ELSE
[3529]706             CALL netcdf_create_global_atts( id_set_mask(mid,av), 'podsmasked', TRIM( run_description_header ), 464 )
[410]707             WRITE ( time_average_text,'(F7.1,'' s avg'')' )  averaging_interval
[3529]708             nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), NF90_GLOBAL, 'time_avg',   &
709                                     TRIM( time_average_text ) )
[1783]710             CALL netcdf_handle_error( 'netcdf_define_header', 466 )
[410]711          ENDIF
712
713!
714!--       Define time coordinate for volume data (unlimited dimension)
[1783]715          CALL netcdf_create_dim( id_set_mask(mid,av), 'time', NF90_UNLIMITED, &
716                                  id_dim_time_mask(mid,av), 467 )
717          CALL netcdf_create_var( id_set_mask(mid,av),                         &
718                                  (/ id_dim_time_mask(mid,av) /), 'time',      &
719                                  NF90_DOUBLE, id_var_time_mask(mid,av),       &
[3966]720                                 'seconds', 'time', 468, 469, 000 )
[3529]721          CALL netcdf_create_att( id_set_mask(mid,av), id_var_time_mask(mid,av), 'standard_name', 'time', 000)
722          CALL netcdf_create_att( id_set_mask(mid,av), id_var_time_mask(mid,av), 'axis', 'T', 000)
723
[410]724!
725!--       Define spatial dimensions and coordinates:
[3435]726          IF ( mask_surface(mid) )  THEN
[410]727!
[3435]728!--          In case of terrain-following output, the vertical dimensions are
729!--          indices, not meters
730             CALL netcdf_create_dim( id_set_mask(mid,av), 'ku_above_surf',     &
731                                     mask_size(mid,3), id_dim_zu_mask(mid,av), &
732                                     470 )
733             CALL netcdf_create_var( id_set_mask(mid,av),                      &
734                                     (/ id_dim_zu_mask(mid,av) /),             &
735                                     'ku_above_surf',                          &
736                                     NF90_DOUBLE, id_var_zu_mask(mid,av),      &
737                                     '1', 'grid point above terrain',          &
738                                     471, 472, 000 )
739             CALL netcdf_create_dim( id_set_mask(mid,av), 'kw_above_surf',     &
740                                     mask_size(mid,3), id_dim_zw_mask(mid,av), &
741                                     473 )
742             CALL netcdf_create_var( id_set_mask(mid,av),                      &
743                                     (/ id_dim_zw_mask(mid,av) /),             &
744                                     'kw_above_surf',                          &
745                                     NF90_DOUBLE, id_var_zw_mask(mid,av),      &
746                                    '1', 'grid point above terrain',           &
747                                    474, 475, 000 )
748          ELSE
[410]749!
[3435]750!--          Define vertical coordinate grid (zu grid)
751             CALL netcdf_create_dim( id_set_mask(mid,av), 'zu_3d',             &
752                                     mask_size(mid,3), id_dim_zu_mask(mid,av), &
753                                     470 )
754             CALL netcdf_create_var( id_set_mask(mid,av),                      &
755                                     (/ id_dim_zu_mask(mid,av) /), 'zu_3d',    &
756                                     NF90_DOUBLE, id_var_zu_mask(mid,av),      &
757                                     'meters', '', 471, 472, 000 )
758!
759!--          Define vertical coordinate grid (zw grid)
760             CALL netcdf_create_dim( id_set_mask(mid,av), 'zw_3d',             &
761                                     mask_size(mid,3), id_dim_zw_mask(mid,av), &
762                                     473 )
763             CALL netcdf_create_var( id_set_mask(mid,av),                      &
764                                     (/ id_dim_zw_mask(mid,av) /), 'zw_3d',    &
765                                     NF90_DOUBLE, id_var_zw_mask(mid,av),      &
766                                    'meters', '', 474, 475, 000 )
767          ENDIF
768!
[410]769!--       Define x-axis (for scalar position)
[1783]770          CALL netcdf_create_dim( id_set_mask(mid,av), 'x', mask_size(mid,1),  &
771                                  id_dim_x_mask(mid,av), 476 )
772          CALL netcdf_create_var( id_set_mask(mid,av),                         &
773                                  (/ id_dim_x_mask(mid,av) /), 'x',            &
774                                  NF90_DOUBLE, id_var_x_mask(mid,av),          &
775                                  'meters', '', 477, 478, 000 )
[410]776!
777!--       Define x-axis (for u position)
[1783]778          CALL netcdf_create_dim( id_set_mask(mid,av), 'xu', mask_size(mid,1), &
779                                  id_dim_xu_mask(mid,av), 479 )
780          CALL netcdf_create_var( id_set_mask(mid,av),                         &
781                                  (/ id_dim_xu_mask(mid,av) /), 'xu',          &
782                                  NF90_DOUBLE, id_var_xu_mask(mid,av),         &
783                                  'meters', '', 480, 481, 000 )
[410]784!
785!--       Define y-axis (for scalar position)
[1783]786          CALL netcdf_create_dim( id_set_mask(mid,av), 'y', mask_size(mid,2),  &
787                                  id_dim_y_mask(mid,av), 482 )
788          CALL netcdf_create_var( id_set_mask(mid,av),                         &
789                                  (/ id_dim_y_mask(mid,av) /), 'y',            &
790                                  NF90_DOUBLE, id_var_y_mask(mid,av),          &
791                                  'meters', '', 483, 484, 000 )
[410]792!
793!--       Define y-axis (for v position)
[1783]794          CALL netcdf_create_dim( id_set_mask(mid,av), 'yv', mask_size(mid,2), &
795                                  id_dim_yv_mask(mid,av), 485 )
796          CALL netcdf_create_var( id_set_mask(mid,av),                         &
797                                  (/ id_dim_yv_mask(mid,av) /),                &
798                                  'yv', NF90_DOUBLE, id_var_yv_mask(mid,av),   &
799                                  'meters', '', 486, 487, 000 )
[410]800!
[3529]801!--       Define UTM and geographic coordinates
802          CALL define_geo_coordinates( id_set_mask(mid,av),               &
803                  (/ id_dim_x_mask(mid,av), id_dim_xu_mask(mid,av) /),    &
804                  (/ id_dim_y_mask(mid,av), id_dim_yv_mask(mid,av) /),    &
[3729]805                  id_var_eutm_mask(:,mid,av), id_var_nutm_mask(:,mid,av), &
806                  id_var_lat_mask(:,mid,av), id_var_lon_mask(:,mid,av)    )
[3421]807!
[3459]808!--       Define coordinate-reference system
809          CALL netcdf_create_crs( id_set_mask(mid,av), 000 )
810!
[410]811!--       In case of non-flat topography define 2d-arrays containing the height
[2232]812!--       information. Only for parallel netcdf output.
813          IF ( TRIM( topography ) /= 'flat'  .AND.                             &
814               netcdf_data_format > 4 )  THEN
[410]815!
816!--          Define zusi = zu(nzb_s_inner)
[1783]817             CALL netcdf_create_var( id_set_mask(mid,av),                      &
818                                     (/ id_dim_x_mask(mid,av),                 &
819                                        id_dim_y_mask(mid,av) /), 'zusi',      &
820                                     NF90_DOUBLE, id_var_zusi_mask(mid,av),    &
821                                     'meters', 'zu(nzb_s_inner)', 488, 489,    &
822                                     490 )
[4046]823!
[410]824!--          Define zwwi = zw(nzb_w_inner)
[1783]825             CALL netcdf_create_var( id_set_mask(mid,av),                      &
826                                     (/ id_dim_x_mask(mid,av),                 &
827                                        id_dim_y_mask(mid,av) /), 'zwwi',      &
828                                     NF90_DOUBLE, id_var_zwwi_mask(mid,av),    &
829                                     'meters', 'zw(nzb_w_inner)', 491, 492,    &
830                                     493 )
[4046]831          ENDIF
832
[1551]833          IF ( land_surface )  THEN
834!
835!--          Define vertical coordinate grid (zw grid)
[1783]836             CALL netcdf_create_dim( id_set_mask(mid,av), 'zs_3d',             &
837                                     mask_size(mid,3), id_dim_zs_mask(mid,av), &
838                                     536 )
839             CALL netcdf_create_var( id_set_mask(mid,av),                      &
840                                     (/ id_dim_zs_mask(mid,av) /), 'zs_3d',    &
841                                     NF90_DOUBLE, id_var_zs_mask(mid,av),      &
842                                     'meters', '', 537, 555, 000 )
[1551]843          ENDIF
844
[410]845!
846!--       Define the variables
847          var_list = ';'
848          i = 1
849
850          DO WHILE ( domask(mid,av,i)(1:1) /= ' ' )
[4046]851
[2007]852             trimvar = TRIM( domask(mid,av,i) )
853!
[410]854!--          Check for the grid
[1980]855             found = .FALSE.
[2007]856             SELECT CASE ( trimvar )
[410]857!
858!--             Most variables are defined on the scalar grid
[3421]859                CASE ( 'e', 'nc', 'nr', 'p', 'pc', 'pr', 'prr',                &
[2292]860                       'q', 'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',   &
[4039]861                       's', 'theta', 'thetal', 'thetav' )
[410]862
863                   grid_x = 'x'
864                   grid_y = 'y'
865                   grid_z = 'zu'
866!
867!--             u grid
868                CASE ( 'u' )
869
870                   grid_x = 'xu'
871                   grid_y = 'y'
872                   grid_z = 'zu'
873!
874!--             v grid
875                CASE ( 'v' )
876
877                   grid_x = 'x'
878                   grid_y = 'yv'
879                   grid_z = 'zu'
880!
881!--             w grid
[1976]882                CASE ( 'w' )
[410]883
884                   grid_x = 'x'
885                   grid_y = 'y'
886                   grid_z = 'zw'
887
[1551]888
[410]889                CASE DEFAULT
[3294]890!
[4046]891!--                Check for quantities defined in other modules
892                   CALL tcm_define_netcdf_grid( trimvar, found, grid_x, grid_y, grid_z )
[3435]893
[3294]894                   IF ( .NOT. found  .AND.  air_chemistry )  THEN
[4046]895                      CALL chem_define_netcdf_grid( trimvar, found, grid_x, grid_y, grid_z )
[3294]896                   ENDIF
[1972]897
[4039]898                   IF ( .NOT. found )                                          &
[4046]899                      CALL doq_define_netcdf_grid( trimvar, found, grid_x, grid_y, grid_z )
900
[3294]901                   IF ( .NOT. found  .AND.  gust_module_enabled )  THEN
[4046]902                      CALL gust_define_netcdf_grid( trimvar, found, grid_x, grid_y, grid_z )
[3294]903                   ENDIF
[2696]904
[4046]905                   IF ( .NOT. found  .AND. land_surface )  THEN
906                      CALL lsm_define_netcdf_grid( trimvar, found, grid_x, grid_y, grid_z )
[1972]907                   ENDIF
[3294]908
909                   IF ( .NOT. found  .AND.  ocean_mode )  THEN
[4046]910                      CALL ocean_define_netcdf_grid( trimvar, found, grid_x, grid_y, grid_z )
[3294]911                   ENDIF
912
[2696]913                   IF ( .NOT. found  .AND.  plant_canopy )  THEN
[4046]914                      CALL pcm_define_netcdf_grid( trimvar, found, grid_x, grid_y, grid_z )
[2209]915                   ENDIF
[1972]916
[1976]917                   IF ( .NOT. found  .AND.  radiation )  THEN
[4046]918                      CALL radiation_define_netcdf_grid( trimvar, found, grid_x, grid_y, grid_z )
[1976]919                   ENDIF
[3337]920!
[3467]921!--                Check for SALSA quantities
922                   IF ( .NOT. found  .AND.  salsa )  THEN
[4046]923                      CALL salsa_define_netcdf_grid( trimvar, found, grid_x, grid_y, grid_z )
924                   ENDIF
925
926                   IF ( .NOT. found  .AND.  urban_surface )  THEN
927                      CALL usm_define_netcdf_grid( trimvar, found, grid_x, grid_y, grid_z )
928                   ENDIF
[3467]929!
[3294]930!--                Now check for user-defined quantities
[3701]931                   IF ( .NOT. found  .AND.  user_module_enabled )  THEN
[4046]932                      CALL user_define_netcdf_grid( trimvar, found, grid_x, grid_y, grid_z )
[1972]933                   ENDIF
[4046]934
[1783]935                   IF ( .NOT. found )  THEN
[4046]936                      WRITE ( message_string, * ) 'no grid defined for variable ', TRIM( trimvar )
937                      CALL message( 'define_netcdf_header', 'PA0244', 0, 1, 0, 6, 0 )
[1783]938                   ENDIF
939
[410]940             END SELECT
941
942!
943!--          Select the respective dimension ids
944             IF ( grid_x == 'x' )  THEN
945                id_x = id_dim_x_mask(mid,av)
946             ELSEIF ( grid_x == 'xu' )  THEN
947                id_x = id_dim_xu_mask(mid,av)
948             ENDIF
949
950             IF ( grid_y == 'y' )  THEN
951                id_y = id_dim_y_mask(mid,av)
952             ELSEIF ( grid_y == 'yv' )  THEN
953                id_y = id_dim_yv_mask(mid,av)
954             ENDIF
955
956             IF ( grid_z == 'zu' )  THEN
957                id_z = id_dim_zu_mask(mid,av)
958             ELSEIF ( grid_z == 'zw' )  THEN
959                id_z = id_dim_zw_mask(mid,av)
[1551]960             ELSEIF ( grid_z == "zs" )  THEN
961                id_z = id_dim_zs_mask(mid,av)
[410]962             ENDIF
963
964!
965!--          Define the grid
[1783]966             CALL netcdf_create_var( id_set_mask(mid,av), (/ id_x, id_y, id_z, &
967                                     id_dim_time_mask(mid,av) /),              &
968                                     domask(mid,av,i), nc_precision(11),       &
969                                     id_var_domask(mid,av,i),                  &
970                                     TRIM( domask_unit(mid,av,i) ),            &
[2696]971                                     domask(mid,av,i), 494, 495, 496, .TRUE. )
[410]972
973             var_list = TRIM( var_list ) // TRIM( domask(mid,av,i) ) // ';'
974
975             i = i + 1
976
977          ENDDO
978
979!
980!--       No arrays to output
981          IF ( i == 1 )  RETURN
982
983!
984!--       Write the list of variables as global attribute (this is used by
985!--       restart runs and by combine_plot_fields)
986          nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), NF90_GLOBAL, &
987                                  'VAR_LIST', var_list )
[1783]988          CALL netcdf_handle_error( 'netcdf_define_header', 497 )
[410]989
990!
[1031]991!--       Leave netCDF define mode
[410]992          nc_stat = NF90_ENDDEF( id_set_mask(mid,av) )
[1783]993          CALL netcdf_handle_error( 'netcdf_define_header', 498 )
[410]994
995!
996!--       Write data for x (shifted by +dx/2) and xu axis
997          ALLOCATE( netcdf_data(mask_size(mid,1)) )
998
[1353]999          netcdf_data = ( mask_i_global(mid,:mask_size(mid,1)) + 0.5_wp ) * dx
[410]1000
1001          nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_x_mask(mid,av), &
1002                                  netcdf_data, start = (/ 1 /),               &
1003                                  count = (/ mask_size(mid,1) /) )
[1783]1004          CALL netcdf_handle_error( 'netcdf_define_header', 499 )
[410]1005
1006          netcdf_data = mask_i_global(mid,:mask_size(mid,1)) * dx
1007
1008          nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_xu_mask(mid,av),&
1009                                  netcdf_data, start = (/ 1 /),               &
1010                                  count = (/ mask_size(mid,1) /) )
[1783]1011          CALL netcdf_handle_error( 'netcdf_define_header', 500 )
[410]1012
1013          DEALLOCATE( netcdf_data )
1014
1015!
1016!--       Write data for y (shifted by +dy/2) and yv axis
1017          ALLOCATE( netcdf_data(mask_size(mid,2)) )
1018
[1353]1019          netcdf_data = ( mask_j_global(mid,:mask_size(mid,2)) + 0.5_wp ) * dy
[410]1020
1021          nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_y_mask(mid,av), &
1022                                  netcdf_data, start = (/ 1 /),               &
1023                                  count = (/ mask_size(mid,2) /))
[1783]1024          CALL netcdf_handle_error( 'netcdf_define_header', 501 )
[410]1025
1026          netcdf_data = mask_j_global(mid,:mask_size(mid,2)) * dy
1027
1028          nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_yv_mask(mid,av), &
1029                                  netcdf_data, start = (/ 1 /),    &
1030                                  count = (/ mask_size(mid,2) /))
[1783]1031          CALL netcdf_handle_error( 'netcdf_define_header', 502 )
[410]1032
1033          DEALLOCATE( netcdf_data )
1034
1035!
[3421]1036!--       Write UTM coordinates
[4196]1037          IF ( rotation_angle == 0.0_wp )  THEN
[3421]1038!
1039!--          1D in case of no rotation
[4196]1040             cos_rot_angle = COS( rotation_angle * pi / 180.0_wp )
[3421]1041!
1042!--          x coordinates
1043             ALLOCATE( netcdf_data(mask_size(mid,1)) )
1044             DO  k = 0, 2
[4046]1045!
[3421]1046!--             Scalar grid points
1047                IF ( k == 0 )  THEN
1048                   shift_x = 0.5
[4046]1049!
[3421]1050!--             u grid points
1051                ELSEIF ( k == 1 )  THEN
1052                   shift_x = 0.0
[4046]1053!
[3421]1054!--             v grid points
1055                ELSEIF ( k == 2 )  THEN
1056                   shift_x = 0.5
1057                ENDIF
1058
[4196]1059                netcdf_data = init_model%origin_x + cos_rot_angle              &
[3421]1060                       * ( mask_i_global(mid,:mask_size(mid,1)) + shift_x ) * dx
1061
1062                nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), &
[3729]1063                                        id_var_eutm_mask(k,mid,av), &
[3421]1064                                        netcdf_data, start = (/ 1 /), &
1065                                        count = (/ mask_size(mid,1) /) )
1066                CALL netcdf_handle_error( 'netcdf_define_header', 555 )
1067
1068             ENDDO
1069             DEALLOCATE( netcdf_data )
1070!
1071!--          y coordinates
1072             ALLOCATE( netcdf_data(mask_size(mid,2)) )
1073             DO  k = 0, 2
1074!
1075!--             Scalar grid points
1076                IF ( k == 0 )  THEN
1077                   shift_y = 0.5
1078!
1079!--             u grid points
1080                ELSEIF ( k == 1 )  THEN
1081                   shift_y = 0.5
1082!
1083!--             v grid points
1084                ELSEIF ( k == 2 )  THEN
1085                   shift_y = 0.0
1086                ENDIF
1087
[4196]1088                netcdf_data = init_model%origin_y + cos_rot_angle              &
[3421]1089                       * ( mask_j_global(mid,:mask_size(mid,2)) + shift_y ) * dy
1090
1091                nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), &
[3729]1092                                        id_var_nutm_mask(k,mid,av), &
[3421]1093                                        netcdf_data, start = (/ 1 /), &
1094                                        count = (/ mask_size(mid,2) /) )
1095                CALL netcdf_handle_error( 'netcdf_define_header', 556 )
1096
1097             ENDDO
1098             DEALLOCATE( netcdf_data )
1099
1100          ELSE
1101!
1102!--          2D in case of rotation
1103             ALLOCATE( netcdf_data_2d(mask_size(mid,1),mask_size(mid,2)) )
[4196]1104             cos_rot_angle = COS( rotation_angle * pi / 180.0_wp )
1105             sin_rot_angle = SIN( rotation_angle * pi / 180.0_wp )
[3421]1106
1107             DO  k = 0, 2
[4046]1108!
[3421]1109!--            Scalar grid points
1110               IF ( k == 0 )  THEN
1111                  shift_x = 0.5 ; shift_y = 0.5
[4046]1112!
[3421]1113!--            u grid points
1114               ELSEIF ( k == 1 )  THEN
1115                  shift_x = 0.0 ; shift_y = 0.5
[4046]1116!
[3421]1117!--            v grid points
1118               ELSEIF ( k == 2 )  THEN
1119                  shift_x = 0.5 ; shift_y = 0.0
1120               ENDIF
1121
1122               DO  j = 1, mask_size(mid,2)
1123                  DO  i = 1, mask_size(mid,1)
[4196]1124                     netcdf_data_2d(i,j) = init_model%origin_x                        &
1125                           + cos_rot_angle * ( mask_i_global(mid,i) + shift_x ) * dx  &
1126                           + sin_rot_angle * ( mask_j_global(mid,j) + shift_y ) * dy
[3421]1127                  ENDDO
1128               ENDDO
1129
1130               nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), &
[3729]1131                                       id_var_eutm_mask(k,mid,av), &
[3421]1132                                       netcdf_data_2d, start = (/ 1, 1 /), &
1133                                       count = (/ mask_size(mid,1), &
1134                                                  mask_size(mid,2) /) )
1135               CALL netcdf_handle_error( 'netcdf_define_header', 555 )
1136
1137               DO  j = 1, mask_size(mid,2)
1138                  DO  i = 1, mask_size(mid,1)
[4196]1139                     netcdf_data_2d(i,j) = init_model%origin_y                        &
1140                           - sin_rot_angle * ( mask_i_global(mid,i) + shift_x ) * dx  &
1141                           + cos_rot_angle * ( mask_j_global(mid,j) + shift_y ) * dy
[3421]1142                  ENDDO
1143               ENDDO
[4046]1144
[3421]1145               nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), &
[3729]1146                                       id_var_nutm_mask(k,mid,av), &
[3421]1147                                       netcdf_data_2d, start = (/ 1, 1 /), &
1148                                       count = (/ mask_size(mid,1), &
1149                                                  mask_size(mid,2) /) )
1150               CALL netcdf_handle_error( 'netcdf_define_header', 556 )
[4046]1151
[3421]1152             ENDDO
1153             DEALLOCATE( netcdf_data_2d )
1154          ENDIF
[3459]1155!
1156!--       Write lon and lat data.
1157          ALLOCATE( lat(mask_size(mid,1),mask_size(mid,2)) )
1158          ALLOCATE( lon(mask_size(mid,1),mask_size(mid,2)) )
[4196]1159          cos_rot_angle = COS( rotation_angle * pi / 180.0_wp )
1160          sin_rot_angle = SIN( rotation_angle * pi / 180.0_wp )
[3421]1161
[3459]1162          DO  k = 0, 2
[4046]1163!
[3459]1164!--          Scalar grid points
1165             IF ( k == 0 )  THEN
1166                shift_x = 0.5 ; shift_y = 0.5
[4046]1167!
[3459]1168!--          u grid points
1169             ELSEIF ( k == 1 )  THEN
1170                shift_x = 0.0 ; shift_y = 0.5
[4046]1171!
[3459]1172!--          v grid points
1173             ELSEIF ( k == 2 )  THEN
1174                shift_x = 0.5 ; shift_y = 0.0
1175             ENDIF
1176
1177             DO  j = 1, mask_size(mid,2)
1178                DO  i = 1, mask_size(mid,1)
[4196]1179                   eutm = init_model%origin_x                                      &
1180                        + cos_rot_angle * ( mask_i_global(mid,i) + shift_x ) * dx  &
1181                        + sin_rot_angle * ( mask_j_global(mid,j) + shift_y ) * dy
1182                   nutm = init_model%origin_y                                      &
1183                        - sin_rot_angle * ( mask_i_global(mid,i) + shift_x ) * dx  &
1184                        + cos_rot_angle * ( mask_j_global(mid,j) + shift_y ) * dy
[3459]1185
1186                   CALL  convert_utm_to_geographic( crs_list,          &
1187                                                    eutm, nutm,        &
1188                                                    lon(i,j), lat(i,j) )
1189                ENDDO
1190             ENDDO
1191
1192             nc_stat = NF90_PUT_VAR( id_set_mask(mid,av),           &
[3729]1193                                     id_var_lon_mask(k,mid,av),     &
[3459]1194                                     lon, start = (/ 1, 1 /),       &
1195                                     count = (/ mask_size(mid,1),   &
1196                                                mask_size(mid,2) /) )
1197             CALL netcdf_handle_error( 'netcdf_define_header', 556 )
1198
1199             nc_stat = NF90_PUT_VAR( id_set_mask(mid,av),           &
[3729]1200                                     id_var_lat_mask(k,mid,av),     &
[3459]1201                                     lat, start = (/ 1, 1 /),       &
1202                                     count = (/ mask_size(mid,1),   &
1203                                                mask_size(mid,2) /) )
1204             CALL netcdf_handle_error( 'netcdf_define_header', 556 )
1205          ENDDO
1206
1207          DEALLOCATE( lat )
1208          DEALLOCATE( lon )
[3421]1209!
[410]1210!--       Write zu and zw data (vertical axes)
1211          ALLOCATE( netcdf_data(mask_size(mid,3)) )
1212
[3435]1213          IF ( mask_surface(mid) )  THEN
[410]1214
[3435]1215             netcdf_data = mask_k_global(mid,:mask_size(mid,3))
[410]1216
[3435]1217             nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_zu_mask(mid,av), &
1218                                     netcdf_data, start = (/ 1 /), &
1219                                     count = (/ mask_size(mid,3) /) )
1220             CALL netcdf_handle_error( 'netcdf_define_header', 503 )
[410]1221
[3435]1222             netcdf_data = mask_k_global(mid,:mask_size(mid,3))
[410]1223
[3435]1224             nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_zw_mask(mid,av), &
1225                                     netcdf_data, start = (/ 1 /), &
1226                                     count = (/ mask_size(mid,3) /) )
1227             CALL netcdf_handle_error( 'netcdf_define_header', 504 )
1228
1229          ELSE
1230
1231             netcdf_data = zu( mask_k_global(mid,:mask_size(mid,3)) )
1232
1233             nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_zu_mask(mid,av), &
1234                                     netcdf_data, start = (/ 1 /), &
1235                                     count = (/ mask_size(mid,3) /) )
1236             CALL netcdf_handle_error( 'netcdf_define_header', 503 )
1237
1238             netcdf_data = zw( mask_k_global(mid,:mask_size(mid,3)) )
1239
1240             nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_zw_mask(mid,av), &
1241                                     netcdf_data, start = (/ 1 /), &
1242                                     count = (/ mask_size(mid,3) /) )
1243             CALL netcdf_handle_error( 'netcdf_define_header', 504 )
1244
1245          ENDIF
1246
[410]1247          DEALLOCATE( netcdf_data )
1248
1249!
1250!--       In case of non-flat topography write height information
[2232]1251          IF ( TRIM( topography ) /= 'flat'  .AND.                             &
1252               netcdf_data_format > 4 )  THEN
[410]1253
[2232]1254             ALLOCATE( netcdf_data_2d(mask_size_l(mid,1),mask_size_l(mid,2)) )
1255             netcdf_data_2d = zu_s_inner( mask_i(mid,:mask_size_l(mid,1)),     &
1256                                          mask_j(mid,:mask_size_l(mid,2)) )
[410]1257
[2232]1258             nc_stat = NF90_PUT_VAR( id_set_mask(mid,av),                      &
1259                                     id_var_zusi_mask(mid,av),                 &
1260                                     netcdf_data_2d,                           &
1261                                     start = (/ 1, 1 /),                       &
1262                                     count = (/ mask_size_l(mid,1),            &
1263                                                mask_size_l(mid,2) /) )
[1783]1264             CALL netcdf_handle_error( 'netcdf_define_header', 505 )
[410]1265
[2232]1266             netcdf_data_2d = zw_w_inner( mask_i(mid,:mask_size_l(mid,1)),     &
1267                                          mask_j(mid,:mask_size_l(mid,2)) )
[410]1268
[2232]1269             nc_stat = NF90_PUT_VAR( id_set_mask(mid,av),                      &
1270                                     id_var_zwwi_mask(mid,av),                 &
1271                                     netcdf_data_2d,                           &
1272                                     start = (/ 1, 1 /),                       &
1273                                     count = (/ mask_size_l(mid,1),            &
1274                                                mask_size_l(mid,2) /) )
[1783]1275             CALL netcdf_handle_error( 'netcdf_define_header', 506 )
[410]1276
1277             DEALLOCATE( netcdf_data_2d )
1278
1279          ENDIF
1280!
[4127]1281!--       soil is not in masked output for now - disable temporary this block
1282!          IF ( land_surface )  THEN
1283!
[1551]1284!--          Write zs data (vertical axes for soil model), use negative values
1285!--          to indicate soil depth
[4127]1286!             ALLOCATE( netcdf_data(mask_size(mid,3)) )
1287!
1288!             netcdf_data = zs( mask_k_global(mid,:mask_size(mid,3)) )
1289!
1290!             nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_zs_mask(mid,av), &
1291!                                     netcdf_data, start = (/ 1 /), &
1292!                                     count = (/ mask_size(mid,3) /) )
1293!             CALL netcdf_handle_error( 'netcdf_define_header', 538 )
1294!
1295!             DEALLOCATE( netcdf_data )
1296!
1297!          ENDIF
[1551]1298
1299!
[410]1300!--       restore original parameter file_id (=formal parameter av) into av
1301          av = file_id
1302
1303
1304       CASE ( 'ma_ext' )
1305
1306!
[4046]1307!--       decompose actual parameter file_id (=formal parameter av) into
[410]1308!--       mid and av
1309          file_id = av
[564]1310          IF ( file_id <= 200+max_masks )  THEN
1311             mid = file_id - 200
[410]1312             av = 0
1313          ELSE
[564]1314             mid = file_id - (200+max_masks)
[410]1315             av = 1
1316          ENDIF
1317
1318!
1319!--       Get the list of variables and compare with the actual run.
1320!--       First var_list_old has to be reset, since GET_ATT does not assign
1321!--       trailing blanks.
1322          var_list_old = ' '
1323          nc_stat = NF90_GET_ATT( id_set_mask(mid,av), NF90_GLOBAL, 'VAR_LIST',&
1324                                  var_list_old )
[1783]1325          CALL netcdf_handle_error( 'netcdf_define_header', 507 )
[410]1326
1327          var_list = ';'
1328          i = 1
1329          DO WHILE ( domask(mid,av,i)(1:1) /= ' ' )
1330             var_list = TRIM(var_list) // TRIM( domask(mid,av,i) ) // ';'
1331             i = i + 1
1332          ENDDO
1333
1334          IF ( av == 0 )  THEN
1335             var = '(mask)'
1336          ELSE
1337             var = '(mask_av)'
1338          ENDIF
1339
1340          IF ( TRIM( var_list ) /= TRIM( var_list_old ) )  THEN
[3045]1341             WRITE ( message_string, * ) 'netCDF file for ', TRIM( var ),       &
1342                  ' data for mask', mid, ' from previous run found,',           &
[3046]1343                  '&but this file cannot be extended due to variable ',         &
1344                  'mismatch.&New file is created instead.'
[564]1345             CALL message( 'define_netcdf_header', 'PA0335', 0, 1, 0, 6, 0 )
[410]1346             extend = .FALSE.
1347             RETURN
1348          ENDIF
1349
1350!
1351!--       Get and compare the number of vertical gridpoints
[1596]1352          nc_stat = NF90_INQ_VARID( id_set_mask(mid,av), 'zu_3d', &
[410]1353                                    id_var_zu_mask(mid,av) )
[1783]1354          CALL netcdf_handle_error( 'netcdf_define_header', 508 )
[410]1355
1356          nc_stat = NF90_INQUIRE_VARIABLE( id_set_mask(mid,av),     &
1357                                           id_var_zu_mask(mid,av),  &
1358                                           dimids = id_dim_zu_mask_old )
[1783]1359          CALL netcdf_handle_error( 'netcdf_define_header', 509 )
[410]1360          id_dim_zu_mask(mid,av) = id_dim_zu_mask_old(1)
1361
[3045]1362          nc_stat = NF90_INQUIRE_DIMENSION( id_set_mask(mid,av),               &
1363                                            id_dim_zu_mask(mid,av),            &
[410]1364                                            len = nz_old )
[1783]1365          CALL netcdf_handle_error( 'netcdf_define_header', 510 )
[410]1366
1367          IF ( mask_size(mid,3) /= nz_old )  THEN
[3045]1368             WRITE ( message_string, * ) 'netCDF file for ', TRIM( var ),      &
[3046]1369                  '&data for mask', mid, ' from previous run found,',          &
[3045]1370                  ' but this file cannot be extended due to mismatch in ',     &
1371                  ' number of vertical grid points.',                          &
[3046]1372                  '&New file is created instead.'
[564]1373             CALL message( 'define_netcdf_header', 'PA0336', 0, 1, 0, 6, 0 )
[410]1374             extend = .FALSE.
1375             RETURN
1376          ENDIF
1377
1378!
1379!--       Get the id of the time coordinate (unlimited coordinate) and its
1380!--       last index on the file. The next time level is plmask..count+1.
1381!--       The current time must be larger than the last output time
1382!--       on the file.
[3045]1383          nc_stat = NF90_INQ_VARID( id_set_mask(mid,av), 'time',               &
[410]1384                                    id_var_time_mask(mid,av) )
[1783]1385          CALL netcdf_handle_error( 'netcdf_define_header', 511 )
[410]1386
[3045]1387          nc_stat = NF90_INQUIRE_VARIABLE( id_set_mask(mid,av),                &
1388                                           id_var_time_mask(mid,av),           &
[410]1389                                           dimids = id_dim_time_old )
[1783]1390          CALL netcdf_handle_error( 'netcdf_define_header', 512 )
[410]1391          id_dim_time_mask(mid,av) = id_dim_time_old(1)
1392
[3045]1393          nc_stat = NF90_INQUIRE_DIMENSION( id_set_mask(mid,av),               &
1394                                            id_dim_time_mask(mid,av),          &
[410]1395                                            len = domask_time_count(mid,av) )
[1783]1396          CALL netcdf_handle_error( 'netcdf_define_header', 513 )
[410]1397
[3045]1398          nc_stat = NF90_GET_VAR( id_set_mask(mid,av),                         &
1399                                  id_var_time_mask(mid,av),                    &
1400                                  last_time_coordinate,                        &
1401                                  start = (/ domask_time_count(mid,av) /),     &
[410]1402                                  count = (/ 1 /) )
[1783]1403          CALL netcdf_handle_error( 'netcdf_define_header', 514 )
[410]1404
1405          IF ( last_time_coordinate(1) >= simulated_time )  THEN
[3045]1406             WRITE ( message_string, * ) 'netCDF file for ', TRIM( var ),      &
1407                  ' data for mask', mid, ' from previous run found,',          &
[3046]1408                  '&but this file cannot be extended because the current ',    &
[3045]1409                  'output time is less or equal than the last output time ',   &
[3046]1410                  'on this file.&New file is created instead.'
[564]1411             CALL message( 'define_netcdf_header', 'PA0337', 0, 1, 0, 6, 0 )
[410]1412             domask_time_count(mid,av) = 0
1413             extend = .FALSE.
1414             RETURN
1415          ENDIF
1416
1417!
1418!--       Dataset seems to be extendable.
1419!--       Now get the variable ids.
1420          i = 1
1421          DO WHILE ( domask(mid,av,i)(1:1) /= ' ' )
1422             nc_stat = NF90_INQ_VARID( id_set_mask(mid,av), &
1423                                       TRIM( domask(mid,av,i) ), &
1424                                       id_var_domask(mid,av,i) )
[1783]1425             CALL netcdf_handle_error( 'netcdf_define_header', 515 )
[410]1426             i = i + 1
1427          ENDDO
1428
1429!
1430!--       Update the title attribute on file
[4046]1431!--       In order to avoid 'data mode' errors if updated attributes are larger
1432!--       than their original size, NF90_PUT_ATT is called in 'define mode'
1433!--       enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a possible
1434!--       performance loss due to data copying; an alternative strategy would be
[410]1435!--       to ensure equal attribute size in a job chain. Maybe revise later.
1436          IF ( av == 0 )  THEN
1437             time_average_text = ' '
1438          ELSE
[3045]1439             WRITE (time_average_text, '('', '',F7.1,'' s average'')')         &
[410]1440                                                            averaging_interval
1441          ENDIF
1442          nc_stat = NF90_REDEF( id_set_mask(mid,av) )
[1783]1443          CALL netcdf_handle_error( 'netcdf_define_header', 516 )
[3045]1444          nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), NF90_GLOBAL, 'title',   &
1445                                  TRIM( run_description_header ) //            &
[410]1446                                  TRIM( time_average_text ) )
[1783]1447          CALL netcdf_handle_error( 'netcdf_define_header', 517 )
[410]1448          nc_stat = NF90_ENDDEF( id_set_mask(mid,av) )
[1783]1449          CALL netcdf_handle_error( 'netcdf_define_header', 518 )
[3045]1450          WRITE ( message_string, * ) 'netCDF file for ', TRIM( var ),         &
1451               ' data for mask', mid, ' from previous run found.',             &
[3046]1452               ' &This file will be extended.'
[564]1453          CALL message( 'define_netcdf_header', 'PA0338', 0, 0, 0, 6, 0 )
[410]1454!
1455!--       restore original parameter file_id (=formal parameter av) into av
1456          av = file_id
1457
1458
[1]1459       CASE ( '3d_new' )
1460
1461!
1462!--       Define some global attributes of the dataset
1463          IF ( av == 0 )  THEN
[3529]1464             CALL netcdf_create_global_atts( id_set_3d(av), '3d', TRIM( run_description_header ), 62 )
[1]1465             time_average_text = ' '
1466          ELSE
[3529]1467             CALL netcdf_create_global_atts( id_set_3d(av), '3d_av', TRIM( run_description_header ), 62 )
[1]1468             WRITE ( time_average_text,'(F7.1,'' s avg'')' )  averaging_interval
[3529]1469             nc_stat = NF90_PUT_ATT( id_set_3d(av), NF90_GLOBAL, 'time_avg',   &
[1]1470                                     TRIM( time_average_text ) )
[1783]1471             CALL netcdf_handle_error( 'netcdf_define_header', 63 )
[1]1472          ENDIF
1473
1474!
[1308]1475!--       Define time coordinate for volume data.
1476!--       For parallel output the time dimensions has to be limited, otherwise
1477!--       the performance drops significantly.
1478          IF ( netcdf_data_format < 5 )  THEN
[1783]1479             CALL netcdf_create_dim( id_set_3d(av), 'time', NF90_UNLIMITED,    &
1480                                     id_dim_time_3d(av), 64 )
[1308]1481          ELSE
[1783]1482             CALL netcdf_create_dim( id_set_3d(av), 'time', ntdim_3d(av),      &
1483                                     id_dim_time_3d(av), 523 )
[1308]1484          ENDIF
[1]1485
[1783]1486          CALL netcdf_create_var( id_set_3d(av), (/ id_dim_time_3d(av) /),     &
1487                                  'time', NF90_DOUBLE, id_var_time_3d(av),     &
[3966]1488                                  'seconds', 'time', 65, 66, 00 )
[3529]1489          CALL netcdf_create_att( id_set_3d(av), id_var_time_3d(av), 'standard_name', 'time', 000)
1490          CALL netcdf_create_att( id_set_3d(av), id_var_time_3d(av), 'axis', 'T', 000)
[1]1491!
1492!--       Define spatial dimensions and coordinates:
1493!--       Define vertical coordinate grid (zu grid)
[1783]1494          CALL netcdf_create_dim( id_set_3d(av), 'zu_3d', nz_do3d-nzb+1,       &
1495                                  id_dim_zu_3d(av), 67 )
1496          CALL netcdf_create_var( id_set_3d(av), (/ id_dim_zu_3d(av) /),       &
1497                                  'zu_3d', NF90_DOUBLE, id_var_zu_3d(av),      &
1498                                  'meters', '', 68, 69, 00 )
[1]1499!
1500!--       Define vertical coordinate grid (zw grid)
[1783]1501          CALL netcdf_create_dim( id_set_3d(av), 'zw_3d', nz_do3d-nzb+1,       &
1502                                  id_dim_zw_3d(av), 70 )
1503          CALL netcdf_create_var( id_set_3d(av), (/ id_dim_zw_3d(av) /),       &
1504                                  'zw_3d', NF90_DOUBLE, id_var_zw_3d(av),      &
1505                                  'meters', '', 71, 72, 00 )
[1]1506!
1507!--       Define x-axis (for scalar position)
[2512]1508          CALL netcdf_create_dim( id_set_3d(av), 'x', nx+1, id_dim_x_3d(av),   &
[1783]1509                                  73 )
1510          CALL netcdf_create_var( id_set_3d(av), (/ id_dim_x_3d(av) /), 'x',   &
1511                                  NF90_DOUBLE, id_var_x_3d(av), 'meters', '',  &
1512                                  74, 75, 00 )
[1]1513!
1514!--       Define x-axis (for u position)
[2512]1515          CALL netcdf_create_dim( id_set_3d(av), 'xu', nx+1, id_dim_xu_3d(av), &
[1783]1516                                  358 )
1517          CALL netcdf_create_var( id_set_3d(av), (/ id_dim_xu_3d(av) /), 'xu', &
1518                                  NF90_DOUBLE, id_var_xu_3d(av), 'meters', '', &
1519                                  359, 360, 000 )
[1]1520!
1521!--       Define y-axis (for scalar position)
[2512]1522          CALL netcdf_create_dim( id_set_3d(av), 'y', ny+1, id_dim_y_3d(av),   &
[1783]1523                                  76 )
1524          CALL netcdf_create_var( id_set_3d(av), (/ id_dim_y_3d(av) /), 'y',   &
1525                                  NF90_DOUBLE, id_var_y_3d(av), 'meters', '',  &
1526                                  77, 78, 00 )
[1]1527!
1528!--       Define y-axis (for v position)
[2512]1529          CALL netcdf_create_dim( id_set_3d(av), 'yv', ny+1, id_dim_yv_3d(av), &
[1783]1530                                  361 )
1531          CALL netcdf_create_var( id_set_3d(av), (/ id_dim_yv_3d(av) /), 'yv', &
1532                                  NF90_DOUBLE, id_var_yv_3d(av), 'meters', '', &
1533                                  362, 363, 000 )
[1]1534!
[3529]1535!--       Define UTM and geographic coordinates
1536          CALL define_geo_coordinates( id_set_3d(av),         &
1537                  (/ id_dim_x_3d(av), id_dim_xu_3d(av) /),    &
1538                  (/ id_dim_y_3d(av), id_dim_yv_3d(av) /),    &
[3729]1539                  id_var_eutm_3d(:,av), id_var_nutm_3d(:,av), &
1540                  id_var_lat_3d(:,av), id_var_lon_3d(:,av)    )
[3421]1541!
[3459]1542!--       Define coordinate-reference system
1543          CALL netcdf_create_crs( id_set_3d(av), 000 )
1544!
[48]1545!--       In case of non-flat topography define 2d-arrays containing the height
[2232]1546!--       information. Only output 2d topography information in case of parallel
1547!--       output.
1548          IF ( TRIM( topography ) /= 'flat'  .AND.                             &
1549               netcdf_data_format > 4 )  THEN
[48]1550!
1551!--          Define zusi = zu(nzb_s_inner)
[1783]1552             CALL netcdf_create_var( id_set_3d(av), (/ id_dim_x_3d(av),        &
1553                                     id_dim_y_3d(av) /), 'zusi', NF90_DOUBLE,  &
1554                                     id_var_zusi_3d(av), 'meters',             &
1555                                     'zu(nzb_s_inner)', 413, 414, 415 )
[4046]1556!
[48]1557!--          Define zwwi = zw(nzb_w_inner)
[1783]1558             CALL netcdf_create_var( id_set_3d(av), (/ id_dim_x_3d(av),        &
1559                                     id_dim_y_3d(av) /), 'zwwi', NF90_DOUBLE,  &
1560                                     id_var_zwwi_3d(av), 'meters',             &
1561                                     'zw(nzb_w_inner)', 416, 417, 418 )
[48]1562
[4046]1563          ENDIF
[48]1564
[1551]1565          IF ( land_surface )  THEN
1566!
1567!--          Define vertical coordinate grid (zs grid)
[1783]1568             CALL netcdf_create_dim( id_set_3d(av), 'zs_3d',                   &
1569                                     nzt_soil-nzb_soil+1, id_dim_zs_3d(av), 70 )
1570             CALL netcdf_create_var( id_set_3d(av), (/ id_dim_zs_3d(av) /),    &
1571                                     'zs_3d', NF90_DOUBLE, id_var_zs_3d(av),   &
1572                                     'meters', '', 71, 72, 00 )
[48]1573
[1551]1574          ENDIF
1575
[4127]1576          IF ( plant_canopy )  THEN
[48]1577!
[4127]1578!--          Define vertical coordinate grid (zpc grid)
1579             CALL netcdf_create_dim( id_set_3d(av), 'zpc_3d',                  &
1580                                     pch_index+1, id_dim_zpc_3d(av), 70 )
1581             !netcdf_create_dim(ncid, dim_name, ncdim_type, ncdim_id, error_no)
1582             CALL netcdf_create_var( id_set_3d(av), (/ id_dim_zpc_3d(av) /),    &
1583                                     'zpc_3d', NF90_DOUBLE, id_var_zpc_3d(av),   &
1584                                     'meters', '', 71, 72, 00 )
1585
1586          ENDIF
1587
1588!
[1]1589!--       Define the variables
1590          var_list = ';'
1591          i = 1
1592
1593          DO WHILE ( do3d(av,i)(1:1) /= ' ' )
1594!
[4046]1595!--          Temporary solution to account for data output within the new urban
[2007]1596!--          surface model (urban_surface_mod.f90), see also SELECT CASE ( trimvar )
1597             trimvar = TRIM( do3d(av,i) )
[2011]1598             IF ( urban_surface  .AND.  trimvar(1:4) == 'usm_' )  THEN
[2007]1599                trimvar = 'usm_output'
1600             ENDIF
1601!
[1]1602!--          Check for the grid
[1980]1603             found = .FALSE.
[2007]1604             SELECT CASE ( trimvar )
[1]1605!
1606!--             Most variables are defined on the scalar grid
[3421]1607                CASE ( 'e', 'nc', 'nr', 'p', 'pc', 'pr', 'prr',   &
[2292]1608                       'q', 'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',   &
[4039]1609                       's', 'theta', 'thetal', 'thetav' )
[1]1610
1611                   grid_x = 'x'
1612                   grid_y = 'y'
1613                   grid_z = 'zu'
1614!
1615!--             u grid
1616                CASE ( 'u' )
1617
1618                   grid_x = 'xu'
1619                   grid_y = 'y'
1620                   grid_z = 'zu'
1621!
1622!--             v grid
1623                CASE ( 'v' )
1624
1625                   grid_x = 'x'
1626                   grid_y = 'yv'
1627                   grid_z = 'zu'
1628!
1629!--             w grid
[1976]1630                CASE ( 'w' )
[1]1631
1632                   grid_x = 'x'
1633                   grid_y = 'y'
1634                   grid_z = 'zw'
1635
[4046]1636!
1637!--             Block of urban surface model outputs
[2007]1638                CASE ( 'usm_output' )
1639                   CALL usm_define_netcdf_grid( do3d(av,i), found, &
1640                                                   grid_x, grid_y, grid_z )
[1551]1641
[1]1642                CASE DEFAULT
[1972]1643
[2696]1644                   CALL tcm_define_netcdf_grid( do3d(av,i), found, &
1645                                                   grid_x, grid_y, grid_z )
1646
[1]1647!
[1972]1648!--                Check for land surface quantities
[3004]1649                   IF ( .NOT. found .AND. land_surface )  THEN
[1972]1650                      CALL lsm_define_netcdf_grid( do3d(av,i), found, grid_x,  &
[1976]1651                                                   grid_y, grid_z )
[1972]1652                   ENDIF
[3421]1653!
1654!--                Check for ocean quantities
1655                   IF ( .NOT. found  .AND.  ocean_mode )  THEN
1656                      CALL ocean_define_netcdf_grid( do3d(av,i), found,  &
1657                                                     grid_x, grid_y, grid_z )
1658                   ENDIF
[1976]1659
1660!
[2209]1661!--                Check for plant canopy quantities
[2696]1662                   IF ( .NOT. found  .AND.  plant_canopy )  THEN
[2209]1663                      CALL pcm_define_netcdf_grid( do3d(av,i), found, grid_x,  &
1664                                                   grid_y, grid_z )
1665                   ENDIF
[3337]1666
[2209]1667!
[1976]1668!--                Check for radiation quantities
1669                   IF ( .NOT. found  .AND.  radiation )  THEN
1670                      CALL radiation_define_netcdf_grid( do3d(av,i), found,    &
1671                                                         grid_x, grid_y,       &
1672                                                         grid_z )
1673                   ENDIF
[2696]1674
[2817]1675!--                Check for gust module quantities
1676                   IF ( .NOT. found  .AND.  gust_module_enabled )  THEN
1677                      CALL gust_define_netcdf_grid( do3d(av,i), found, grid_x, &
1678                                                    grid_y, grid_z )
1679                   ENDIF
[3744]1680!
[4046]1681!--                Check for indoor model quantities
1682                   IF ( .NOT. found .AND. indoor_model ) THEN
[3744]1683                      CALL im_define_netcdf_grid( do3d(av,i), found,           &
[3766]1684                                                  grid_x, grid_y, grid_z )
[3744]1685                   ENDIF
[2817]1686
[2696]1687!
[3448]1688!--                Check for biometeorology quantities
1689                   IF ( .NOT. found  .AND.  biometeorology )  THEN
[3525]1690                      CALL bio_define_netcdf_grid( do3d(av,i), found,          &
1691                                                   grid_x, grid_y, grid_z )
[3448]1692                   ENDIF
1693
1694!
[4046]1695!--                Check for chemistry quantities
[2696]1696                   IF ( .NOT. found  .AND.  air_chemistry )  THEN
1697                      CALL chem_define_netcdf_grid( do3d(av,i), found,         &
1698                                                    grid_x, grid_y, grid_z )
1699                   ENDIF
1700
[3467]1701!
1702!--                Check for SALSA quantities
1703                   IF ( .NOT. found  .AND.  salsa )  THEN
1704                      CALL salsa_define_netcdf_grid( do3d(av,i), found, grid_x,&
1705                                                     grid_y, grid_z )
[4046]1706                   ENDIF
1707!
[1]1708!--                Check for user-defined quantities
[3701]1709                   IF ( .NOT. found  .AND.  user_module_enabled )  THEN
[1972]1710                      CALL user_define_netcdf_grid( do3d(av,i), found, grid_x, &
1711                                                    grid_y, grid_z )
1712                   ENDIF
[4046]1713
[4039]1714                   IF ( .NOT. found )                                          &
1715                      CALL doq_define_netcdf_grid( do3d(av,i), found, grid_x,  &
1716                                                   grid_y, grid_z        )
[4046]1717
[1783]1718                   IF ( .NOT. found )  THEN
1719                      WRITE ( message_string, * ) 'no grid defined for varia', &
1720                                                  'ble ', TRIM( do3d(av,i) )
1721                      CALL message( 'define_netcdf_header', 'PA0244', 0, 1, 0, &
1722                                    6, 0 )
1723                   ENDIF
1724
[1]1725             END SELECT
1726
1727!
1728!--          Select the respective dimension ids
1729             IF ( grid_x == 'x' )  THEN
1730                id_x = id_dim_x_3d(av)
1731             ELSEIF ( grid_x == 'xu' )  THEN
1732                id_x = id_dim_xu_3d(av)
1733             ENDIF
1734
1735             IF ( grid_y == 'y' )  THEN
1736                id_y = id_dim_y_3d(av)
1737             ELSEIF ( grid_y == 'yv' )  THEN
1738                id_y = id_dim_yv_3d(av)
1739             ENDIF
1740
1741             IF ( grid_z == 'zu' )  THEN
1742                id_z = id_dim_zu_3d(av)
1743             ELSEIF ( grid_z == 'zw' )  THEN
1744                id_z = id_dim_zw_3d(av)
[1551]1745             ELSEIF ( grid_z == 'zs' )  THEN
1746                id_z = id_dim_zs_3d(av)
[4127]1747             ELSEIF ( grid_z == 'zpc' )  THEN
1748                id_z = id_dim_zpc_3d(av)
[1]1749             ENDIF
1750
1751!
1752!--          Define the grid
[1783]1753             CALL netcdf_create_var( id_set_3d(av),(/ id_x, id_y, id_z,        &
1754                                     id_dim_time_3d(av) /), do3d(av,i),        &
1755                                     nc_precision(4), id_var_do3d(av,i),       &
1756                                     TRIM( do3d_unit(av,i) ), do3d(av,i), 79,  &
[2696]1757                                     80, 357, .TRUE. )
[1031]1758#if defined( __netcdf4_parallel )
[1308]1759             IF ( netcdf_data_format > 4 )  THEN
[493]1760!
[1308]1761!--             Set no fill for every variable to increase performance.
1762                nc_stat = NF90_DEF_VAR_FILL( id_set_3d(av),     &
1763                                             id_var_do3d(av,i), &
[4029]1764                                             NF90_NOFILL, 0 )
[1783]1765                CALL netcdf_handle_error( 'netcdf_define_header', 532 )
[1308]1766!
1767!--             Set collective io operations for parallel io
[493]1768                nc_stat = NF90_VAR_PAR_ACCESS( id_set_3d(av),     &
1769                                               id_var_do3d(av,i), &
1770                                               NF90_COLLECTIVE )
[1783]1771                CALL netcdf_handle_error( 'netcdf_define_header', 445 )
[493]1772             ENDIF
1773#endif
[1783]1774             var_list = TRIM( var_list ) // TRIM( do3d(av,i) ) // ';'
1775
[1]1776             i = i + 1
1777
1778          ENDDO
1779
1780!
1781!--       No arrays to output
1782          IF ( i == 1 )  RETURN
1783
1784!
1785!--       Write the list of variables as global attribute (this is used by
1786!--       restart runs and by combine_plot_fields)
1787          nc_stat = NF90_PUT_ATT( id_set_3d(av), NF90_GLOBAL, 'VAR_LIST', &
1788                                  var_list )
[1783]1789          CALL netcdf_handle_error( 'netcdf_define_header', 81 )
[1]1790
1791!
[1308]1792!--       Set general no fill, otherwise the performance drops significantly for
1793!--       parallel output.
1794          nc_stat = NF90_SET_FILL( id_set_3d(av), NF90_NOFILL, oldmode )
[1783]1795          CALL netcdf_handle_error( 'netcdf_define_header', 528 )
[1308]1796
1797!
[1031]1798!--       Leave netCDF define mode
[1]1799          nc_stat = NF90_ENDDEF( id_set_3d(av) )
[1783]1800          CALL netcdf_handle_error( 'netcdf_define_header', 82 )
[1]1801
1802!
[1308]1803!--       These data are only written by PE0 for parallel output to increase
1804!--       the performance.
1805          IF ( myid == 0  .OR.  netcdf_data_format < 5 )  THEN
1806!
1807!--          Write data for x (shifted by +dx/2) and xu axis
[2512]1808             ALLOCATE( netcdf_data(0:nx) )
[1]1809
[2512]1810             DO  i = 0, nx
[1308]1811                netcdf_data(i) = ( i + 0.5 ) * dx
1812             ENDDO
[1]1813
[1308]1814             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_x_3d(av),  &
1815                                     netcdf_data, start = (/ 1 /),    &
[2512]1816                                     count = (/ nx+1 /) )
[1783]1817             CALL netcdf_handle_error( 'netcdf_define_header', 83 )
[1]1818
[2512]1819             DO  i = 0, nx
[1308]1820                netcdf_data(i) = i * dx
1821             ENDDO
[1]1822
[1308]1823             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_xu_3d(av), &
1824                                     netcdf_data, start = (/ 1 /),    &
[2512]1825                                     count = (/ nx+1 /) )
[1783]1826             CALL netcdf_handle_error( 'netcdf_define_header', 385 )
[1]1827
[1308]1828             DEALLOCATE( netcdf_data )
[1]1829
1830!
[1308]1831!--          Write data for y (shifted by +dy/2) and yv axis
[2512]1832             ALLOCATE( netcdf_data(0:ny) )
[1]1833
[2512]1834             DO  i = 0, ny
[1353]1835                netcdf_data(i) = ( i + 0.5_wp ) * dy
[1308]1836             ENDDO
[1]1837
[1308]1838             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_y_3d(av),  &
1839                                     netcdf_data, start = (/ 1 /),    &
[2512]1840                                     count = (/ ny+1 /) )
[1783]1841             CALL netcdf_handle_error( 'netcdf_define_header', 84 )
[1]1842
[2512]1843             DO  i = 0, ny
[1308]1844                netcdf_data(i) = i * dy
1845             ENDDO
[1]1846
[1308]1847             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_yv_3d(av), &
1848                                     netcdf_data, start = (/ 1 /),    &
[2512]1849                                     count = (/ ny+1 /))
[1783]1850             CALL netcdf_handle_error( 'netcdf_define_header', 387 )
[1]1851
[1308]1852             DEALLOCATE( netcdf_data )
[1]1853
1854!
[3421]1855!--          Write UTM coordinates
[4196]1856             IF ( rotation_angle == 0.0_wp )  THEN
[3421]1857!
1858!--             1D in case of no rotation
[4196]1859                cos_rot_angle = COS( rotation_angle * pi / 180.0_wp )
[3421]1860!
1861!--             x coordinates
1862                ALLOCATE( netcdf_data(0:nx) )
1863                DO  k = 0, 2
[4046]1864!
[3421]1865!--                Scalar grid points
1866                   IF ( k == 0 )  THEN
1867                      shift_x = 0.5
[4046]1868!
[3421]1869!--                u grid points
1870                   ELSEIF ( k == 1 )  THEN
1871                      shift_x = 0.0
[4046]1872!
[3421]1873!--                v grid points
1874                   ELSEIF ( k == 2 )  THEN
1875                      shift_x = 0.5
1876                   ENDIF
[4046]1877
[3421]1878                   DO  i = 0, nx
[4196]1879                     netcdf_data(i) = init_model%origin_x                      &
1880                                    + cos_rot_angle * ( i + shift_x ) * dx
[3421]1881                   ENDDO
[4046]1882
[3729]1883                   nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_eutm_3d(k,av),&
[4196]1884                                           netcdf_data, start = (/ 1 /),       &
[3421]1885                                           count = (/ nx+1 /) )
1886                   CALL netcdf_handle_error( 'netcdf_define_header', 555 )
1887
1888                ENDDO
1889                DEALLOCATE( netcdf_data )
1890!
1891!--             y coordinates
1892                ALLOCATE( netcdf_data(0:ny) )
1893                DO  k = 0, 2
1894!
1895!--                Scalar grid points
1896                   IF ( k == 0 )  THEN
1897                      shift_y = 0.5
1898!
1899!--                u grid points
1900                   ELSEIF ( k == 1 )  THEN
1901                      shift_y = 0.5
1902!
1903!--                v grid points
1904                   ELSEIF ( k == 2 )  THEN
1905                      shift_y = 0.0
1906                   ENDIF
1907
1908                   DO  j = 0, ny
[4196]1909                      netcdf_data(j) = init_model%origin_y                     &
1910                                     + cos_rot_angle * ( j + shift_y ) * dy
[3421]1911                   ENDDO
1912
[3729]1913                   nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_nutm_3d(k,av),&
[4196]1914                                           netcdf_data, start = (/ 1 /),       &
[3421]1915                                           count = (/ ny+1 /) )
1916                   CALL netcdf_handle_error( 'netcdf_define_header', 556 )
1917
1918                ENDDO
1919                DEALLOCATE( netcdf_data )
1920
1921             ELSE
1922!
1923!--             2D in case of rotation
1924                ALLOCATE( netcdf_data_2d(0:nx,0:ny) )
[4196]1925                cos_rot_angle = COS( rotation_angle * pi / 180.0_wp )
1926                sin_rot_angle = SIN( rotation_angle * pi / 180.0_wp )
[4046]1927
[3421]1928                DO  k = 0, 2
[4046]1929!
[3421]1930!--               Scalar grid points
1931                  IF ( k == 0 )  THEN
1932                     shift_x = 0.5 ; shift_y = 0.5
[4046]1933!
[3421]1934!--               u grid points
1935                  ELSEIF ( k == 1 )  THEN
1936                     shift_x = 0.0 ; shift_y = 0.5
[4046]1937!
[3421]1938!--               v grid points
1939                  ELSEIF ( k == 2 )  THEN
1940                     shift_x = 0.5 ; shift_y = 0.0
1941                  ENDIF
[4046]1942
[3421]1943                  DO  j = 0, ny
1944                     DO  i = 0, nx
[4196]1945                        netcdf_data_2d(i,j) = init_model%origin_x                   &
1946                                            + cos_rot_angle * ( i + shift_x ) * dx  &
1947                                            + sin_rot_angle * ( j + shift_y ) * dy
[3421]1948                     ENDDO
1949                  ENDDO
[4046]1950
[3729]1951                  nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_eutm_3d(k,av),  &
[3421]1952                                          netcdf_data_2d, start = (/ 1, 1 /),   &
1953                                          count = (/ nx+1, ny+1 /) )
1954                  CALL netcdf_handle_error( 'netcdf_define_header', 555 )
[4046]1955
[3421]1956                  DO  j = 0, ny
1957                     DO  i = 0, nx
[4196]1958                        netcdf_data_2d(i,j) = init_model%origin_y                   &
1959                                            - sin_rot_angle * ( i + shift_x ) * dx  &
1960                                            + cos_rot_angle * ( j + shift_y ) * dy
[3421]1961                     ENDDO
1962                  ENDDO
[4046]1963
[3729]1964                  nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_nutm_3d(k,av),  &
[3421]1965                                          netcdf_data_2d, start = (/ 1, 1 /),   &
1966                                          count = (/ nx+1, ny+1 /) )
1967                  CALL netcdf_handle_error( 'netcdf_define_header', 556 )
[4046]1968
[3421]1969                ENDDO
1970                DEALLOCATE( netcdf_data_2d )
1971             ENDIF
1972!
[1308]1973!--          Write zu and zw data (vertical axes)
1974             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zu_3d(av),  &
1975                                     zu(nzb:nz_do3d), start = (/ 1 /), &
1976                                     count = (/ nz_do3d-nzb+1 /) )
[1783]1977             CALL netcdf_handle_error( 'netcdf_define_header', 85 )
[1]1978
[263]1979
[1308]1980             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zw_3d(av),  &
1981                                     zw(nzb:nz_do3d), start = (/ 1 /), &
1982                                     count = (/ nz_do3d-nzb+1 /) )
[1783]1983             CALL netcdf_handle_error( 'netcdf_define_header', 86 )
[1]1984
[1551]1985             IF ( land_surface )  THEN
1986!
1987!--             Write zs grid
1988                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zs_3d(av),  &
1989                                        - zs(nzb_soil:nzt_soil), start = (/ 1 /), &
1990                                        count = (/ nzt_soil-nzb_soil+1 /) )
[1783]1991                CALL netcdf_handle_error( 'netcdf_define_header', 86 )
[1551]1992             ENDIF
1993
[4127]1994             IF ( plant_canopy )  THEN
1995!
1996!--             Write zpc grid
1997                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zpc_3d(av),  &
1998                                        zu(nzb:nzb+pch_index), start = (/ 1 /), &
1999                                        count = (/ pch_index+1 /) )
2000                CALL netcdf_handle_error( 'netcdf_define_header', 86 )
2001             ENDIF
2002
[48]2003          ENDIF
[2232]2004!
[3459]2005!--       Write lon and lat data. Only for parallel output.
2006          IF ( netcdf_data_format > 4 )  THEN
2007
2008             ALLOCATE( lat(nxl:nxr,nys:nyn) )
2009             ALLOCATE( lon(nxl:nxr,nys:nyn) )
[4196]2010             cos_rot_angle = COS( rotation_angle * pi / 180.0_wp )
2011             sin_rot_angle = SIN( rotation_angle * pi / 180.0_wp )
[3459]2012
2013             DO  k = 0, 2
[4046]2014!
[3459]2015!--             Scalar grid points
2016                IF ( k == 0 )  THEN
2017                   shift_x = 0.5 ; shift_y = 0.5
[4046]2018!
[3459]2019!--             u grid points
2020                ELSEIF ( k == 1 )  THEN
2021                   shift_x = 0.0 ; shift_y = 0.5
[4046]2022!
[3459]2023!--             v grid points
2024                ELSEIF ( k == 2 )  THEN
2025                   shift_x = 0.5 ; shift_y = 0.0
2026                ENDIF
2027
2028                DO  j = nys, nyn
2029                   DO  i = nxl, nxr
[4196]2030                      eutm = init_model%origin_x                   &
2031                           + cos_rot_angle * ( i + shift_x ) * dx  &
2032                           + sin_rot_angle * ( j + shift_y ) * dy
2033                      nutm = init_model%origin_y                   &
2034                           - sin_rot_angle * ( i + shift_x ) * dx  &
2035                           + cos_rot_angle * ( j + shift_y ) * dy
[3459]2036
2037                      CALL  convert_utm_to_geographic( crs_list,          &
2038                                                       eutm, nutm,        &
2039                                                       lon(i,j), lat(i,j) )
2040                   ENDDO
2041                ENDDO
2042
[3729]2043                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_lon_3d(k,av), &
[3459]2044                                     lon, start = (/ nxl+1, nys+1 /),       &
2045                                     count = (/ nxr-nxl+1, nyn-nys+1 /) )
2046                CALL netcdf_handle_error( 'netcdf_define_header', 556 )
2047
[3729]2048                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_lat_3d(k,av), &
[3459]2049                                     lat, start = (/ nxl+1, nys+1 /),       &
2050                                     count = (/ nxr-nxl+1, nyn-nys+1 /) )
2051                CALL netcdf_handle_error( 'netcdf_define_header', 556 )
2052             ENDDO
2053
2054             DEALLOCATE( lat )
2055             DEALLOCATE( lon )
2056
2057          ENDIF
2058!
[2232]2059!--       In case of non-flat topography write height information. Only for
2060!--       parallel netcdf output.
2061          IF ( TRIM( topography ) /= 'flat'  .AND.                             &
2062               netcdf_data_format > 4 )  THEN
[48]2063
[2512]2064!             IF ( nxr == nx  .AND.  nyn /= ny )  THEN
2065!                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zusi_3d(av),     &
2066!                                        zu_s_inner(nxl:nxr+1,nys:nyn),         &
2067!                                        start = (/ nxl+1, nys+1 /),            &
2068!                                        count = (/ nxr-nxl+2, nyn-nys+1 /) )
2069!             ELSEIF ( nxr /= nx  .AND.  nyn == ny )  THEN
2070!                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zusi_3d(av),     &
2071!                                        zu_s_inner(nxl:nxr,nys:nyn+1),         &
2072!                                        start = (/ nxl+1, nys+1 /),            &
2073!                                        count = (/ nxr-nxl+1, nyn-nys+2 /) )
2074!             ELSEIF ( nxr == nx  .AND.  nyn == ny )  THEN
2075!                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zusi_3d(av),     &
2076!                                        zu_s_inner(nxl:nxr+1,nys:nyn+1),       &
2077!                                        start = (/ nxl+1, nys+1 /),            &
2078!                                        count = (/ nxr-nxl+2, nyn-nys+2 /) )
2079!             ELSE
[2232]2080                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zusi_3d(av),     &
2081                                        zu_s_inner(nxl:nxr,nys:nyn),           &
2082                                        start = (/ nxl+1, nys+1 /),            &
2083                                        count = (/ nxr-nxl+1, nyn-nys+1 /) )
[2512]2084!             ENDIF
[2232]2085             CALL netcdf_handle_error( 'netcdf_define_header', 419 )
2086
[2512]2087!             IF ( nxr == nx  .AND.  nyn /= ny )  THEN
2088!                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zwwi_3d(av),     &
2089!                                        zw_w_inner(nxl:nxr+1,nys:nyn),         &
2090!                                        start = (/ nxl+1, nys+1 /),            &
2091!                                        count = (/ nxr-nxl+2, nyn-nys+1 /) )
2092!             ELSEIF ( nxr /= nx  .AND.  nyn == ny )  THEN
2093!                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zwwi_3d(av),     &
2094!                                        zw_w_inner(nxl:nxr,nys:nyn+1),         &
2095!                                        start = (/ nxl+1, nys+1 /),            &
2096!                                        count = (/ nxr-nxl+1, nyn-nys+2 /) )
2097!             ELSEIF ( nxr == nx  .AND.  nyn == ny )  THEN
2098!                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zwwi_3d(av),     &
2099!                                        zw_w_inner(nxl:nxr+1,nys:nyn+1),       &
2100!                                        start = (/ nxl+1, nys+1 /),            &
2101!                                        count = (/ nxr-nxl+2, nyn-nys+2 /) )
2102!             ELSE
[2232]2103                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zwwi_3d(av),     &
2104                                        zw_w_inner(nxl:nxr,nys:nyn),           &
2105                                        start = (/ nxl+1, nys+1 /),            &
2106                                        count = (/ nxr-nxl+1, nyn-nys+1 /) )
[2512]2107!             ENDIF
[2232]2108             CALL netcdf_handle_error( 'netcdf_define_header', 420 )
2109
2110          ENDIF
2111
[1]2112       CASE ( '3d_ext' )
2113
2114!
2115!--       Get the list of variables and compare with the actual run.
2116!--       First var_list_old has to be reset, since GET_ATT does not assign
2117!--       trailing blanks.
2118          var_list_old = ' '
2119          nc_stat = NF90_GET_ATT( id_set_3d(av), NF90_GLOBAL, 'VAR_LIST', &
2120                                  var_list_old )
[1783]2121          CALL netcdf_handle_error( 'netcdf_define_header', 87 )
[1]2122
2123          var_list = ';'
2124          i = 1
2125          DO WHILE ( do3d(av,i)(1:1) /= ' ' )
2126             var_list = TRIM(var_list) // TRIM( do3d(av,i) ) // ';'
2127             i = i + 1
2128          ENDDO
2129
2130          IF ( av == 0 )  THEN
2131             var = '(3d)'
2132          ELSE
2133             var = '(3d_av)'
2134          ENDIF
2135
2136          IF ( TRIM( var_list ) /= TRIM( var_list_old ) )  THEN
[1031]2137             message_string = 'netCDF file for volume data ' //             &
[292]2138                              TRIM( var ) // ' from previous run found,' // &
[3046]2139                              '&but this file cannot be extended due to' // &
[257]2140                              ' variable mismatch.' //                      &
[3046]2141                              '&New file is created instead.'
[257]2142             CALL message( 'define_netcdf_header', 'PA0245', 0, 1, 0, 6, 0 )
[1]2143             extend = .FALSE.
2144             RETURN
2145          ENDIF
2146
2147!
2148!--       Get and compare the number of vertical gridpoints
2149          nc_stat = NF90_INQ_VARID( id_set_3d(av), 'zu_3d', id_var_zu_3d(av) )
[1783]2150          CALL netcdf_handle_error( 'netcdf_define_header', 88 )
[1]2151
2152          nc_stat = NF90_INQUIRE_VARIABLE( id_set_3d(av), id_var_zu_3d(av), &
2153                                           dimids = id_dim_zu_3d_old )
[1783]2154          CALL netcdf_handle_error( 'netcdf_define_header', 89 )
[1]2155          id_dim_zu_3d(av) = id_dim_zu_3d_old(1)
2156
2157          nc_stat = NF90_INQUIRE_DIMENSION( id_set_3d(av), id_dim_zu_3d(av), &
2158                                            len = nz_old )
[1783]2159          CALL netcdf_handle_error( 'netcdf_define_header', 90 )
[1]2160
2161          IF ( nz_do3d-nzb+1 /= nz_old )  THEN
[1031]2162              message_string = 'netCDF file for volume data ' //             &
[292]2163                               TRIM( var ) // ' from previous run found,' // &
[3046]2164                               '&but this file cannot be extended due to' // &
[257]2165                               ' mismatch in number of' //                   &
[3045]2166                               ' vertical grid points (nz_do3d).' //         &
[3046]2167                               '&New file is created instead.'
[257]2168             CALL message( 'define_netcdf_header', 'PA0246', 0, 1, 0, 6, 0 )
[1]2169             extend = .FALSE.
2170             RETURN
2171          ENDIF
2172
2173!
2174!--       Get the id of the time coordinate (unlimited coordinate) and its
2175!--       last index on the file. The next time level is pl3d..count+1.
2176!--       The current time must be larger than the last output time
2177!--       on the file.
2178          nc_stat = NF90_INQ_VARID( id_set_3d(av), 'time', id_var_time_3d(av) )
[1783]2179          CALL netcdf_handle_error( 'netcdf_define_header', 91 )
[1]2180
2181          nc_stat = NF90_INQUIRE_VARIABLE( id_set_3d(av), id_var_time_3d(av), &
2182                                           dimids = id_dim_time_old )
[1783]2183          CALL netcdf_handle_error( 'netcdf_define_header', 92 )
[263]2184
[1]2185          id_dim_time_3d(av) = id_dim_time_old(1)
2186
2187          nc_stat = NF90_INQUIRE_DIMENSION( id_set_3d(av), id_dim_time_3d(av), &
[1308]2188                                            len = ntime_count )
[1783]2189          CALL netcdf_handle_error( 'netcdf_define_header', 93 )
[1]2190
[1308]2191!
2192!--       For non-parallel output use the last output time level of the netcdf
2193!--       file because the time dimension is unlimited. In case of parallel
2194!--       output the variable ntime_count could get the value of 9*10E36 because
2195!--       the time dimension is limited.
2196          IF ( netcdf_data_format < 5 ) do3d_time_count(av) = ntime_count
2197
[1]2198          nc_stat = NF90_GET_VAR( id_set_3d(av), id_var_time_3d(av), &
2199                                  last_time_coordinate,              &
2200                                  start = (/ do3d_time_count(av) /), &
2201                                  count = (/ 1 /) )
[1783]2202          CALL netcdf_handle_error( 'netcdf_define_header', 94 )
[1]2203
2204          IF ( last_time_coordinate(1) >= simulated_time )  THEN
[3045]2205             message_string = 'netCDF file for volume data ' //             &
[292]2206                              TRIM( var ) // ' from previous run found,' // &
[3046]2207                              '&but this file cannot be extended becaus' // &
[257]2208                              'e the current output time' //                &
[3046]2209                              '&is less or equal than the last output t' // &
[257]2210                              'ime on this file.' //                        &
[3046]2211                              '&New file is created instead.'
[257]2212             CALL message( 'define_netcdf_header', 'PA0247', 0, 1, 0, 6, 0 )
[1]2213             do3d_time_count(av) = 0
2214             extend = .FALSE.
2215             RETURN
2216          ENDIF
2217
[1308]2218          IF ( netcdf_data_format > 4 )  THEN
[1745]2219!
2220!--          Check if the needed number of output time levels is increased
2221!--          compared to the number of time levels in the existing file.
[1308]2222             IF ( ntdim_3d(av) > ntime_count )  THEN
2223                message_string = 'netCDF file for volume data ' // &
2224                                 TRIM( var ) // ' from previous run found,' // &
[3046]2225                                 '&but this file cannot be extended becaus' // &
[3045]2226                                 'e the number of output time levels has b' // &
[1308]2227                                 'een increased compared to the previous s' // &
[1745]2228                                 'imulation.' //                               &
[3046]2229                                 '&New file is created instead.'
[1308]2230                CALL message( 'define_netcdf_header', 'PA0388', 0, 1, 0, 6, 0 )
2231                do3d_time_count(av) = 0
2232                extend = .FALSE.
[1745]2233!
2234!--             Recalculate the needed time levels for the new file.
2235                IF ( av == 0 )  THEN
2236                   ntdim_3d(0) = CEILING(                               &
2237                           ( end_time - MAX( skip_time_do3d,            &
2238                                             simulated_time_at_begin )  &
2239                           ) / dt_do3d )
[2769]2240                   IF ( do3d_at_begin )  ntdim_3d(0) = ntdim_3d(0) + 1
[1745]2241                ELSE
2242                   ntdim_3d(1) = CEILING(                               &
2243                           ( end_time - MAX( skip_time_data_output_av,  &
2244                                             simulated_time_at_begin )  &
2245                           ) / dt_data_output_av )
2246                ENDIF
[1308]2247                RETURN
2248             ENDIF
2249          ENDIF
2250
[1]2251!
2252!--       Dataset seems to be extendable.
2253!--       Now get the variable ids.
2254          i = 1
2255          DO WHILE ( do3d(av,i)(1:1) /= ' ' )
2256             nc_stat = NF90_INQ_VARID( id_set_3d(av), TRIM( do3d(av,i) ), &
2257                                       id_var_do3d(av,i) )
[1783]2258             CALL netcdf_handle_error( 'netcdf_define_header', 95 )
[1031]2259#if defined( __netcdf4_parallel )
[493]2260!
2261!--          Set collective io operations for parallel io
[1031]2262             IF ( netcdf_data_format > 4 )  THEN
[493]2263                nc_stat = NF90_VAR_PAR_ACCESS( id_set_3d(av),     &
2264                                               id_var_do3d(av,i), &
2265                                               NF90_COLLECTIVE )
[1783]2266                CALL netcdf_handle_error( 'netcdf_define_header', 453 )
[493]2267             ENDIF
2268#endif
[1]2269             i = i + 1
2270          ENDDO
2271
2272!
[359]2273!--       Update the title attribute on file
[4046]2274!--       In order to avoid 'data mode' errors if updated attributes are larger
2275!--       than their original size, NF90_PUT_ATT is called in 'define mode'
2276!--       enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a possible
2277!--       performance loss due to data copying; an alternative strategy would be
[359]2278!--       to ensure equal attribute size. Maybe revise later.
2279          IF ( av == 0 )  THEN
2280             time_average_text = ' '
2281          ELSE
2282             WRITE (time_average_text, '('', '',F7.1,'' s average'')') &
2283                                                            averaging_interval
2284          ENDIF
2285          nc_stat = NF90_REDEF( id_set_3d(av) )
[1783]2286          CALL netcdf_handle_error( 'netcdf_define_header', 429 )
[1]2287          nc_stat = NF90_PUT_ATT( id_set_3d(av), NF90_GLOBAL, 'title', &
[359]2288                                  TRIM( run_description_header ) //    &
2289                                  TRIM( time_average_text ) )
[1783]2290          CALL netcdf_handle_error( 'netcdf_define_header', 96 )
[359]2291          nc_stat = NF90_ENDDEF( id_set_3d(av) )
[1783]2292          CALL netcdf_handle_error( 'netcdf_define_header', 430 )
[1031]2293          message_string = 'netCDF file for volume data ' //             &
[257]2294                           TRIM( var ) // ' from previous run found.' // &
[3046]2295                           '&This file will be extended.'
[257]2296          CALL message( 'define_netcdf_header', 'PA0248', 0, 0, 0, 6, 0 )
[1]2297
[3159]2298
2299       CASE ( 'ag_new' )
2300
2301!
2302!--       Define some global attributes of the dataset
2303          nc_stat = NF90_PUT_ATT( id_set_agt, NF90_GLOBAL, 'title', &
2304                                  TRIM( run_description_header ) )
2305          CALL netcdf_handle_error( 'netcdf_define_header', 330 )
2306!
2307!--       Switch for unlimited time dimension
2308          IF ( agent_time_unlimited ) THEN
2309             CALL netcdf_create_dim( id_set_agt, 'time', NF90_UNLIMITED,       &
2310                                     id_dim_time_agt, 331 )
2311          ELSE
2312             CALL netcdf_create_dim( id_set_agt, 'time',                       &
[3198]2313                                     INT( ( MIN( multi_agent_system_end,       &
2314                                                 end_time ) -                  &
2315                                            multi_agent_system_start ) /       &
2316                                            dt_write_agent_data * 1.1 ),       &
[3159]2317                                     id_dim_time_agt, 331 )
2318          ENDIF
2319
2320          CALL netcdf_create_var( id_set_agt, (/ id_dim_time_agt /), 'time',   &
[3966]2321                                  NF90_REAL4, id_var_time_agt, 'seconds', 'time',  &
[3159]2322                                  332, 333, 000 )
[3529]2323          CALL netcdf_create_att( id_set_agt, id_var_time_agt, 'standard_name', 'time', 000)
2324          CALL netcdf_create_att( id_set_agt, id_var_time_agt, 'axis', 'T', 000)
[3159]2325
[3235]2326          CALL netcdf_create_dim( id_set_agt, 'agent_number',                  &
2327                                  dim_size_agtnum, id_dim_agtnum, 334 )
2328
[3159]2329          CALL netcdf_create_var( id_set_agt, (/ id_dim_agtnum /),             &
[3235]2330                                  'agent_number', NF90_REAL4,                  &
2331                                  id_var_agtnum, 'agent number', '', 335,      &
[3159]2332                                  336, 000 )
2333!
2334!--       Define variable which contains the real number of agents in use
2335          CALL netcdf_create_var( id_set_agt, (/ id_dim_time_agt /),           &
[3235]2336                                  'real_num_of_agt', NF90_REAL4,               &
2337                                  id_var_rnoa_agt, 'agent number', '', 337,    &
[3159]2338                                  338, 000 )
[3235]2339          i = 1
2340          CALL netcdf_create_var( id_set_agt, (/ id_dim_agtnum,                &
2341                                  id_dim_time_agt /), agt_var_names(i),        &
2342                                  NF90_DOUBLE, id_var_agt(i),                  &
2343                                  TRIM( agt_var_units(i) ),                    &
2344                                  TRIM( agt_var_names(i) ), 339, 340, 341 )
[3159]2345!
2346!--       Define the variables
[3235]2347          DO  i = 2, 6
[3159]2348             CALL netcdf_create_var( id_set_agt, (/ id_dim_agtnum,             &
2349                                     id_dim_time_agt /), agt_var_names(i),     &
[3235]2350                                     NF90_REAL4, id_var_agt(i),                &
[3159]2351                                     TRIM( agt_var_units(i) ),                 &
2352                                     TRIM( agt_var_names(i) ), 339, 340, 341 )
2353
2354          ENDDO
2355!
[3448]2356!--       Define vars for biometeorology
2357          CALL netcdf_create_var( id_set_agt, (/ id_dim_agtnum,                &
2358                                  id_dim_time_agt /), agt_var_names(9),        &
2359                                  nc_precision(8), id_var_agt(9),              &
2360                                  TRIM( agt_var_units(9) ),                    &
2361                                  TRIM( agt_var_names(9) ), 339, 340, 341 )
2362
2363!
[3159]2364!--       Leave netCDF define mode
2365          nc_stat = NF90_ENDDEF( id_set_agt )
2366          CALL netcdf_handle_error( 'netcdf_define_header', 342 )
2367
2368
2369!        CASE ( 'ag_ext' )
2370! !+?agent extend output for restart runs has to be adapted
[4046]2371!
[3159]2372! !
2373! !--       Get the id of the time coordinate (unlimited coordinate) and its
2374! !--       last index on the file. The next time level is prt..count+1.
2375! !--       The current time must be larger than the last output time
2376! !--       on the file.
2377!           nc_stat = NF90_INQ_VARID( id_set_agt, 'time', id_var_time_agt )
2378!           CALL netcdf_handle_error( 'netcdf_define_header', 343 )
[4046]2379!
[3159]2380!           nc_stat = NF90_INQUIRE_VARIABLE( id_set_agt, id_var_time_agt, &
2381!                                            dimids = id_dim_time_old )
2382!           CALL netcdf_handle_error( 'netcdf_define_header', 344 )
2383!           id_dim_time_agt = id_dim_time_old(1)
[4046]2384!
[3159]2385!           nc_stat = NF90_INQUIRE_DIMENSION( id_set_agt, id_dim_time_agt, &
2386!                                             len = agt_time_count )
2387!           CALL netcdf_handle_error( 'netcdf_define_header', 345 )
[4046]2388!
[3159]2389!           nc_stat = NF90_GET_VAR( id_set_agt, id_var_time_agt,  &
2390!                                   last_time_coordinate,         &
2391!                                   start = (/ agt_time_count /), &
2392!                                   count = (/ 1 /) )
2393!           CALL netcdf_handle_error( 'netcdf_define_header', 346 )
[4046]2394!
[3159]2395!           IF ( last_time_coordinate(1) >= simulated_time )  THEN
2396!              message_string = 'netCDF file for agents ' //                  &
2397!                               'from previous run found,' //                 &
2398!                               '&but this file cannot be extended becaus' // &
2399!                               'e the current output time' //                &
2400!                               '&is less or equal than the last output t' // &
2401!                               'ime on this file.' //                        &
2402!                               '&New file is created instead.'
2403!              CALL message( 'define_netcdf_header', 'PA0265', 0, 1, 0, 6, 0 )
2404!              agt_time_count = 0
2405!              extend = .FALSE.
2406!              RETURN
2407!           ENDIF
[4046]2408!
[3159]2409! !
2410! !--       Dataset seems to be extendable.
2411! !--       Now get the variable ids.
2412!           nc_stat = NF90_INQ_VARID( id_set_agt, 'real_num_of_agt', &
2413!                                     id_var_rnoa_agt )
2414!           CALL netcdf_handle_error( 'netcdf_define_header', 347 )
[4046]2415!
[3159]2416!           DO  i = 1, 17
[4046]2417!
[3159]2418!              nc_stat = NF90_INQ_VARID( id_set_agt, agt_var_names(i), &
2419!                                        id_var_prt(i) )
2420!              CALL netcdf_handle_error( 'netcdf_define_header', 348 )
[4046]2421!
[3159]2422!           ENDDO
[4046]2423!
[3159]2424!           message_string = 'netCDF file for particles ' // &
2425!                            'from previous run found.' //   &
2426!                            '&This file will be extended.'
2427!           CALL message( 'define_netcdf_header', 'PA0266', 0, 0, 0, 6, 0 )
2428
[4046]2429
[1]2430       CASE ( 'xy_new' )
2431
2432!
2433!--       Define some global attributes of the dataset
2434          IF ( av == 0 )  THEN
[3529]2435             CALL netcdf_create_global_atts( id_set_xy(av), 'xy', TRIM( run_description_header ), 97 )
[1]2436             time_average_text = ' '
2437          ELSE
[3529]2438             CALL netcdf_create_global_atts( id_set_xy(av), 'xy_av', TRIM( run_description_header ), 97 )
[1]2439             WRITE ( time_average_text,'(F7.1,'' s avg'')' )  averaging_interval
[3529]2440             nc_stat = NF90_PUT_ATT( id_set_xy(av), NF90_GLOBAL, 'time_avg',   &
[1]2441                                     TRIM( time_average_text ) )
[1783]2442             CALL netcdf_handle_error( 'netcdf_define_header', 98 )
[1]2443          ENDIF
2444
2445!
[1308]2446!--       Define time coordinate for xy sections.
2447!--       For parallel output the time dimensions has to be limited, otherwise
2448!--       the performance drops significantly.
2449          IF ( netcdf_data_format < 5 )  THEN
[1783]2450             CALL netcdf_create_dim( id_set_xy(av), 'time', NF90_UNLIMITED,    &
2451                                     id_dim_time_xy(av), 99 )
[1308]2452          ELSE
[1783]2453             CALL netcdf_create_dim( id_set_xy(av), 'time', ntdim_2d_xy(av),   &
2454                                     id_dim_time_xy(av), 524 )
[1308]2455          ENDIF
[1]2456
[1783]2457          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_time_xy(av) /),     &
2458                                  'time', NF90_DOUBLE, id_var_time_xy(av),     &
[3966]2459                                  'seconds', 'time', 100, 101, 000 )
[3529]2460          CALL netcdf_create_att( id_set_xy(av), id_var_time_xy(av), 'standard_name', 'time', 000)
2461          CALL netcdf_create_att( id_set_xy(av), id_var_time_xy(av), 'axis', 'T', 000)
[1]2462!
2463!--       Define the spatial dimensions and coordinates for xy-sections.
2464!--       First, determine the number of horizontal sections to be written.
2465          IF ( section(1,1) == -9999 )  THEN
2466             RETURN
2467          ELSE
2468             ns = 1
2469             DO WHILE ( section(ns,1) /= -9999  .AND.  ns <= 100 )
2470                ns = ns + 1
2471             ENDDO
2472             ns = ns - 1
2473          ENDIF
2474
2475!
2476!--       Define vertical coordinate grid (zu grid)
[1783]2477          CALL netcdf_create_dim( id_set_xy(av), 'zu_xy', ns,                  &
2478                                  id_dim_zu_xy(av), 102 )
2479          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_zu_xy(av) /),       &
2480                                  'zu_xy', NF90_DOUBLE, id_var_zu_xy(av),      &
2481                                  'meters', '', 103, 104, 000 )
[1]2482!
2483!--       Define vertical coordinate grid (zw grid)
[1783]2484          CALL netcdf_create_dim( id_set_xy(av), 'zw_xy', ns,                  &
2485                                  id_dim_zw_xy(av), 105 )
2486          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_zw_xy(av) /),       &
2487                                  'zw_xy', NF90_DOUBLE, id_var_zw_xy(av),      &
2488                                  'meters', '', 106, 107, 000 )
[1]2489
[1551]2490          IF ( land_surface )  THEN
2491
[2232]2492             ns_do = 1
[2239]2493             DO WHILE ( section(ns_do,1) /= -9999  .AND.  ns_do < nzs )
[1551]2494                ns_do = ns_do + 1
2495             ENDDO
[1]2496!
[1551]2497!--          Define vertical coordinate grid (zs grid)
[1783]2498             CALL netcdf_create_dim( id_set_xy(av), 'zs_xy', ns_do,            &
2499                                     id_dim_zs_xy(av), 539 )
2500             CALL netcdf_create_var( id_set_xy(av), (/ id_dim_zs_xy(av) /),    &
2501                                     'zs_xy', NF90_DOUBLE, id_var_zs_xy(av),   &
2502                                     'meters', '', 540, 541, 000 )
[1551]2503
2504          ENDIF
2505
2506!
[1]2507!--       Define a pseudo vertical coordinate grid for the surface variables
2508!--       u* and t* to store their height level
[1783]2509          CALL netcdf_create_dim( id_set_xy(av), 'zu1_xy', 1,                  &
2510                                  id_dim_zu1_xy(av), 108 )
2511          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_zu1_xy(av) /),      &
2512                                  'zu1_xy', NF90_DOUBLE, id_var_zu1_xy(av),    &
2513                                  'meters', '', 109, 110, 000 )
[1]2514!
2515!--       Define a variable to store the layer indices of the horizontal cross
2516!--       sections, too
[1783]2517          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_zu_xy(av) /),       &
2518                                  'ind_z_xy', NF90_DOUBLE,                     &
2519                                  id_var_ind_z_xy(av), 'gridpoints', '', 111,  &
2520                                  112, 000 )
[1]2521!
2522!--       Define x-axis (for scalar position)
[2512]2523          CALL netcdf_create_dim( id_set_xy(av), 'x', nx+1, id_dim_x_xy(av),   &
[1783]2524                                  113 )
2525          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_x_xy(av) /), 'x',   &
2526                                  NF90_DOUBLE, id_var_x_xy(av), 'meters', '',  &
2527                                  114, 115, 000 )
[1]2528!
2529!--       Define x-axis (for u position)
[2512]2530          CALL netcdf_create_dim( id_set_xy(av), 'xu', nx+1,                   &
[1783]2531                                  id_dim_xu_xy(av), 388 )
2532          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_xu_xy(av) /), 'xu', &
2533                                  NF90_DOUBLE, id_var_xu_xy(av), 'meters', '', &
2534                                  389, 390, 000 )
[1]2535!
2536!--       Define y-axis (for scalar position)
[2512]2537          CALL netcdf_create_dim( id_set_xy(av), 'y', ny+1, id_dim_y_xy(av),   &
[1783]2538                                  116 )
2539          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_y_xy(av) /), 'y',   &
2540                                  NF90_DOUBLE, id_var_y_xy(av), 'meters', '',  &
2541                                  117, 118, 000 )
[1]2542!
2543!--       Define y-axis (for scalar position)
[2512]2544          CALL netcdf_create_dim( id_set_xy(av), 'yv', ny+1,                   &
[1783]2545                                  id_dim_yv_xy(av), 364 )
2546          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_yv_xy(av) /), 'yv', &
2547                                  NF90_DOUBLE, id_var_yv_xy(av), 'meters', '', &
2548                                  365, 366, 000 )
[1]2549!
[3529]2550!--       Define UTM and geographic coordinates
2551          CALL define_geo_coordinates( id_set_xy(av),         &
2552                  (/ id_dim_x_xy(av), id_dim_xu_xy(av) /),    &
2553                  (/ id_dim_y_xy(av), id_dim_yv_xy(av) /),    &
[3729]2554                  id_var_eutm_xy(:,av), id_var_nutm_xy(:,av), &
2555                  id_var_lat_xy(:,av), id_var_lon_xy(:,av)    )
[3421]2556!
[3459]2557!--       Define coordinate-reference system
2558          CALL netcdf_create_crs( id_set_xy(av), 000 )
2559!
[48]2560!--       In case of non-flat topography define 2d-arrays containing the height
[2232]2561!--       information. Only for parallel netcdf output.
2562          IF ( TRIM( topography ) /= 'flat'  .AND.                             &
2563               netcdf_data_format > 4  )  THEN
[48]2564!
2565!--          Define zusi = zu(nzb_s_inner)
[1783]2566             CALL netcdf_create_var( id_set_xy(av), (/ id_dim_x_xy(av),        &
2567                                     id_dim_y_xy(av) /), 'zusi', NF90_DOUBLE,  &
2568                                     id_var_zusi_xy(av), 'meters',             &
2569                                     'zu(nzb_s_inner)', 421, 422, 423 )
[4046]2570!
[48]2571!--          Define zwwi = zw(nzb_w_inner)
[1783]2572             CALL netcdf_create_var( id_set_xy(av), (/ id_dim_x_xy(av),        &
2573                                     id_dim_y_xy(av) /), 'zwwi', NF90_DOUBLE,  &
2574                                     id_var_zwwi_xy(av), 'meters',             &
2575                                     'zw(nzb_w_inner)', 424, 425, 426 )
[48]2576
2577          ENDIF
2578
2579!
[1]2580!--       Define the variables
2581          var_list = ';'
2582          i = 1
2583
2584          DO WHILE ( do2d(av,i)(1:1) /= ' ' )
2585
2586             IF ( INDEX( do2d(av,i), 'xy' ) /= 0 )  THEN
2587!
2588!--             If there is a star in the variable name (u* or t*), it is a
2589!--             surface variable. Define it with id_dim_zu1_xy.
2590                IF ( INDEX( do2d(av,i), '*' ) /= 0 )  THEN
2591
[1783]2592                   CALL netcdf_create_var( id_set_xy(av), (/ id_dim_x_xy(av),  &
2593                                           id_dim_y_xy(av), id_dim_zu1_xy(av), &
2594                                           id_dim_time_xy(av) /), do2d(av,i),  &
2595                                           nc_precision(1), id_var_do2d(av,i), &
2596                                           TRIM( do2d_unit(av,i) ),            &
[2696]2597                                           do2d(av,i), 119, 120, 354, .TRUE. )
[1]2598
2599                ELSE
2600
2601!
2602!--                Check for the grid
[1980]2603                   found = .FALSE.
[1]2604                   SELECT CASE ( do2d(av,i) )
2605!
2606!--                   Most variables are defined on the zu grid
[3421]2607                      CASE ( 'e_xy', 'nc_xy', 'nr_xy', 'p_xy',                 &
2608                             'pc_xy', 'pr_xy', 'prr_xy', 'q_xy',               &
[2292]2609                             'qc_xy', 'ql_xy', 'ql_c_xy', 'ql_v_xy',           &
[3421]2610                             'ql_vp_xy', 'qr_xy', 'qv_xy',                     &
2611                             's_xy',                                           &
[4039]2612                             'theta_xy', 'thetal_xy', 'thetav_xy' )
[1]2613
2614                         grid_x = 'x'
2615                         grid_y = 'y'
2616                         grid_z = 'zu'
2617!
2618!--                   u grid
2619                      CASE ( 'u_xy' )
2620
2621                         grid_x = 'xu'
2622                         grid_y = 'y'
2623                         grid_z = 'zu'
2624!
2625!--                   v grid
2626                      CASE ( 'v_xy' )
2627
2628                         grid_x = 'x'
2629                         grid_y = 'yv'
2630                         grid_z = 'zu'
2631!
2632!--                   w grid
[1976]2633                      CASE ( 'w_xy' )
[1]2634
2635                         grid_x = 'x'
2636                         grid_y = 'y'
2637                         grid_z = 'zw'
2638
[1976]2639
[1]2640                      CASE DEFAULT
2641!
[1976]2642!--                      Check for land surface quantities
2643                         IF ( land_surface )  THEN
2644                            CALL lsm_define_netcdf_grid( do2d(av,i), found,    &
2645                                                   grid_x, grid_y, grid_z )
2646                         ENDIF
2647
[2696]2648                         IF ( .NOT. found )  THEN
2649                            CALL tcm_define_netcdf_grid( do2d(av,i), found,    &
2650                                                         grid_x, grid_y,       &
2651                                                         grid_z )
2652                         ENDIF
2653
[1976]2654!
[3421]2655!--                      Check for ocean quantities
2656                         IF ( .NOT. found  .AND.  ocean_mode )  THEN
2657                            CALL ocean_define_netcdf_grid( do2d(av,i), found,  &
2658                                                           grid_x, grid_y,     &
2659                                                           grid_z )
2660                         ENDIF
2661!
[1976]2662!--                      Check for radiation quantities
2663                         IF ( .NOT. found  .AND.  radiation )  THEN
2664                            CALL radiation_define_netcdf_grid( do2d(av,i),     &
2665                                                         found, grid_x, grid_y,&
2666                                                         grid_z )
2667                         ENDIF
[3582]2668
[3467]2669!
2670!--                      Check for SALSA quantities
2671                         IF ( .NOT. found  .AND.  salsa )  THEN
2672                            CALL salsa_define_netcdf_grid( do2d(av,i), found,  &
2673                                                           grid_x, grid_y,     &
2674                                                           grid_z )
[4046]2675                         ENDIF
[1976]2676
2677!
[2817]2678!--                      Check for gust module quantities
2679                         IF ( .NOT. found  .AND.  gust_module_enabled )  THEN
2680                            CALL gust_define_netcdf_grid( do2d(av,i), found,   &
2681                                                          grid_x, grid_y,      &
2682                                                          grid_z )
2683                         ENDIF
2684!
[3569]2685!--                      Check for biometeorology quantities
[3448]2686                         IF ( .NOT. found  .AND.  biometeorology )  THEN
[3525]2687                            CALL bio_define_netcdf_grid( do2d( av, i), found,  &
2688                                                         grid_x, grid_y,       &
2689                                                         grid_z )
[3448]2690                         ENDIF
2691!
[2696]2692!--                      Check for chemistry quantities
2693                         IF ( .NOT. found  .AND.  air_chemistry )  THEN
2694                            CALL chem_define_netcdf_grid( do2d(av,i), found,   &
2695                                                          grid_x, grid_y,      &
2696                                                          grid_z )
2697                         ENDIF
[4046]2698
[4039]2699                         IF ( .NOT. found )                                    &
2700                            CALL doq_define_netcdf_grid(                       &
2701                                                    do2d(av,i), found, grid_x, &
2702                                                    grid_y, grid_z              )
[2696]2703!
[1]2704!--                      Check for user-defined quantities
[3701]2705                         IF ( .NOT. found  .AND.  user_module_enabled )  THEN
[1976]2706                            CALL user_define_netcdf_grid( do2d(av,i), found,   &
2707                                                          grid_x, grid_y,      &
2708                                                          grid_z )
2709                         ENDIF
[1]2710
[1783]2711                         IF ( .NOT. found )  THEN
2712                            WRITE ( message_string, * ) 'no grid defined for', &
2713                                                ' variable ', TRIM( do2d(av,i) )
2714                            CALL message( 'define_netcdf_header', 'PA0244',    &
2715                                          0, 1, 0, 6, 0 )
2716                         ENDIF
2717
[1]2718                   END SELECT
2719
2720!
2721!--                Select the respective dimension ids
2722                   IF ( grid_x == 'x' )  THEN
2723                      id_x = id_dim_x_xy(av)
2724                   ELSEIF ( grid_x == 'xu' )  THEN
2725                      id_x = id_dim_xu_xy(av)
2726                   ENDIF
2727
2728                   IF ( grid_y == 'y' )  THEN
2729                      id_y = id_dim_y_xy(av)
2730                   ELSEIF ( grid_y == 'yv' )  THEN
2731                      id_y = id_dim_yv_xy(av)
2732                   ENDIF
2733
2734                   IF ( grid_z == 'zu' )  THEN
2735                      id_z = id_dim_zu_xy(av)
2736                   ELSEIF ( grid_z == 'zw' )  THEN
2737                      id_z = id_dim_zw_xy(av)
[1551]2738                   ELSEIF ( grid_z == 'zs' )  THEN
2739                      id_z = id_dim_zs_xy(av)
[3448]2740                   ELSEIF ( grid_z == 'zu1' )  THEN
2741                      id_z = id_dim_zu1_xy(av)
[1]2742                   ENDIF
2743
2744!
2745!--                Define the grid
[1783]2746                   CALL netcdf_create_var( id_set_xy(av), (/ id_x, id_y, id_z, &
2747                                           id_dim_time_xy(av) /), do2d(av,i),  &
2748                                           nc_precision(1), id_var_do2d(av,i), &
2749                                           TRIM( do2d_unit(av,i) ),            &
[2696]2750                                           do2d(av,i), 119, 120, 354, .TRUE. )
[1]2751
2752                ENDIF
2753
[1031]2754#if defined( __netcdf4_parallel )
[1308]2755                IF ( netcdf_data_format > 4 )  THEN
[493]2756!
[1308]2757!--                Set no fill for every variable to increase performance.
2758                   nc_stat = NF90_DEF_VAR_FILL( id_set_xy(av),     &
2759                                                id_var_do2d(av,i), &
[4029]2760                                                NF90_NOFILL, 0 )
[1783]2761                   CALL netcdf_handle_error( 'netcdf_define_header', 533 )
[1308]2762!
2763!--                Set collective io operations for parallel io
[493]2764                   nc_stat = NF90_VAR_PAR_ACCESS( id_set_xy(av),     &
2765                                                  id_var_do2d(av,i), &
2766                                                  NF90_COLLECTIVE )
[1783]2767                   CALL netcdf_handle_error( 'netcdf_define_header', 448 )
[493]2768                ENDIF
2769#endif
[1783]2770                var_list = TRIM( var_list) // TRIM( do2d(av,i) ) // ';'
2771
[1]2772             ENDIF
2773
2774             i = i + 1
2775
2776          ENDDO
2777
2778!
2779!--       No arrays to output. Close the netcdf file and return.
2780          IF ( i == 1 )  RETURN
2781
2782!
2783!--       Write the list of variables as global attribute (this is used by
2784!--       restart runs and by combine_plot_fields)
2785          nc_stat = NF90_PUT_ATT( id_set_xy(av), NF90_GLOBAL, 'VAR_LIST', &
2786                                  var_list )
[1783]2787          CALL netcdf_handle_error( 'netcdf_define_header', 121 )
[1]2788
2789!
[1308]2790!--       Set general no fill, otherwise the performance drops significantly for
2791!--       parallel output.
2792          nc_stat = NF90_SET_FILL( id_set_xy(av), NF90_NOFILL, oldmode )
[1783]2793          CALL netcdf_handle_error( 'netcdf_define_header', 529 )
[1308]2794
2795!
[1031]2796!--       Leave netCDF define mode
[1]2797          nc_stat = NF90_ENDDEF( id_set_xy(av) )
[1783]2798          CALL netcdf_handle_error( 'netcdf_define_header', 122 )
[1]2799
2800!
[1308]2801!--       These data are only written by PE0 for parallel output to increase
2802!--       the performance.
2803          IF ( myid == 0  .OR.  netcdf_data_format < 5 )  THEN
[1]2804
2805!
[1308]2806!--          Write axis data: z_xy, x, y
2807             ALLOCATE( netcdf_data(1:ns) )
[1]2808
2809!
[1308]2810!--          Write zu data
2811             DO  i = 1, ns
2812                IF( section(i,1) == -1 )  THEN
[1353]2813                   netcdf_data(i) = -1.0_wp  ! section averaged along z
[1308]2814                ELSE
2815                   netcdf_data(i) = zu( section(i,1) )
2816                ENDIF
2817             ENDDO
2818             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zu_xy(av), &
2819                                     netcdf_data, start = (/ 1 /),    &
2820                                     count = (/ ns /) )
[1783]2821             CALL netcdf_handle_error( 'netcdf_define_header', 123 )
[1]2822
2823!
[1308]2824!--          Write zw data
2825             DO  i = 1, ns
2826                IF( section(i,1) == -1 )  THEN
[1353]2827                   netcdf_data(i) = -1.0_wp  ! section averaged along z
[1308]2828                ELSE
2829                   netcdf_data(i) = zw( section(i,1) )
2830                ENDIF
2831             ENDDO
2832             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zw_xy(av), &
2833                                     netcdf_data, start = (/ 1 /),    &
2834                                     count = (/ ns /) )
[1783]2835             CALL netcdf_handle_error( 'netcdf_define_header', 124 )
[1]2836
[1308]2837!
[1972]2838!--          Write zs data
[1551]2839             IF ( land_surface )  THEN
2840                ns_do = 0
2841                DO  i = 1, ns
2842                   IF( section(i,1) == -1 )  THEN
2843                      netcdf_data(i) = 1.0_wp  ! section averaged along z
2844                      ns_do = ns_do + 1
2845                   ELSEIF ( section(i,1) < nzs )  THEN
2846                      netcdf_data(i) = - zs( section(i,1) )
2847                      ns_do = ns_do + 1
2848                   ENDIF
2849                ENDDO
2850
2851                nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zs_xy(av), &
2852                                        netcdf_data(1:ns_do), start = (/ 1 /),    &
2853                                        count = (/ ns_do /) )
[1783]2854                CALL netcdf_handle_error( 'netcdf_define_header', 124 )
[1551]2855
2856             ENDIF
2857
2858!
[1308]2859!--          Write gridpoint number data
2860             netcdf_data(1:ns) = section(1:ns,1)
2861             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_ind_z_xy(av), &
2862                                     netcdf_data, start = (/ 1 /),       &
2863                                     count = (/ ns /) )
[1783]2864             CALL netcdf_handle_error( 'netcdf_define_header', 125 )
[1]2865
[1308]2866             DEALLOCATE( netcdf_data )
2867
[1]2868!
[1308]2869!--          Write the cross section height u*, t*
2870             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zu1_xy(av), &
2871                                     (/ zu(nzb+1) /), start = (/ 1 /), &
2872                                     count = (/ 1 /) )
[1783]2873             CALL netcdf_handle_error( 'netcdf_define_header', 126 )
[1]2874
2875!
[1308]2876!--          Write data for x (shifted by +dx/2) and xu axis
[2512]2877             ALLOCATE( netcdf_data(0:nx) )
[1]2878
[2512]2879             DO  i = 0, nx
[1353]2880                netcdf_data(i) = ( i + 0.5_wp ) * dx
[1308]2881             ENDDO
[1]2882
[1308]2883             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_x_xy(av), &
2884                                     netcdf_data, start = (/ 1 /),   &
[2512]2885                                     count = (/ nx+1 /) )
[1783]2886             CALL netcdf_handle_error( 'netcdf_define_header', 127 )
[1]2887
[2512]2888             DO  i = 0, nx
[1308]2889                netcdf_data(i) = i * dx
2890             ENDDO
[1]2891
[1308]2892             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_xu_xy(av), &
2893                                     netcdf_data, start = (/ 1 /),    &
[2512]2894                                     count = (/ nx+1 /) )
[1783]2895             CALL netcdf_handle_error( 'netcdf_define_header', 367 )
[1]2896
[1308]2897             DEALLOCATE( netcdf_data )
[1]2898
2899!
[1308]2900!--          Write data for y (shifted by +dy/2) and yv axis
2901             ALLOCATE( netcdf_data(0:ny+1) )
[1]2902
[2512]2903             DO  i = 0, ny
[1353]2904                netcdf_data(i) = ( i + 0.5_wp ) * dy
[1308]2905             ENDDO
[1]2906
[1308]2907             nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_y_xy(av), &
2908                                     netcdf_data, start = (/ 1 /),   &
[2512]2909                                     count = (/ ny+1 /))
[1783]2910             CALL netcdf_handle_error( 'netcdf_define_header', 128 )
[1]2911
[2512]2912