Index: palm/trunk/SOURCE/netcdf_interface.f90
===================================================================
--- palm/trunk/SOURCE/netcdf.f90	(revision 1746)
+++ palm/trunk/SOURCE/netcdf_interface.f90	(revision 1783)
@@ -1,3 +1,3 @@
-!> @file netcdf.f90
+!> @file netcdf_interface.f90
 !--------------------------------------------------------------------------------!
 ! This file is part of PALM.
@@ -19,5 +19,9 @@
 ! Current revisions:
 ! ------------------
-! 
+! netcdf interface has been modularized, former file netcdf renamed to
+! netcdf_interface, creation of netcdf-dimensions and -variables moved to
+! specific new subroutines create_netcdf_dim and create_netcdf_var,
+! compression (deflation) of variables implemented,
+! ibmy special cpp directive removed
 ! 
 ! Former revisions:
@@ -120,8 +124,210 @@
 !>       cover every exception (change of dt_do, end_time in restart)
 !------------------------------------------------------------------------------!
-#if defined( __ibmy_special )
-@PROCESS NOOPTimize
+ MODULE netcdf_interface
+
+    USE control_parameters, ONLY: max_masks
+    USE kinds
+#if defined( __netcdf )
+    USE NETCDF
 #endif
- SUBROUTINE define_netcdf_header( callmode, extend, av )
+
+    PRIVATE
+
+    INTEGER(iwp), PARAMETER ::  dopr_norm_num = 7, dopts_num = 29, dots_max = 100
+
+    CHARACTER (LEN=6), DIMENSION(dopr_norm_num) ::  dopr_norm_names =          &
+         (/ 'wpt0  ', 'ws2   ', 'tsw2  ', 'ws3   ', 'ws2tsw', 'wstsw2',        &
+            'z_i   ' /)
+
+    CHARACTER (LEN=6), DIMENSION(dopr_norm_num) ::  dopr_norm_longnames =      &
+         (/ 'wpt0  ', 'w*2   ', 't*w2  ', 'w*3   ', 'w*2t*w', 'w*t*w2',        &
+            'z_i   ' /)
+
+    CHARACTER (LEN=7), DIMENSION(dopts_num) :: dopts_label =                   &
+          (/ 'tnpt   ', 'x_     ', 'y_     ', 'z_     ', 'z_abs  ', 'u      ', &
+             'v      ', 'w      ', 'u"     ', 'v"     ', 'w"     ', 'npt_up ', &
+             'w_up   ', 'w_down ', 'radius ', 'r_min  ', 'r_max  ', 'npt_max', &
+             'npt_min', 'x*2    ', 'y*2    ', 'z*2    ', 'u*2    ', 'v*2    ', &
+             'w*2    ', 'u"2    ', 'v"2    ', 'w"2    ', 'npt*2  ' /)
+
+    CHARACTER (LEN=7), DIMENSION(dopts_num) :: dopts_unit =                    &
+          (/ 'number ', 'm      ', 'm      ', 'm      ', 'm      ', 'm/s    ', &
+             'm/s    ', 'm/s    ', 'm/s    ', 'm/s    ', 'm/s    ', 'number ', &
+             'm/s    ', 'm/s    ', 'm      ', 'm      ', 'm      ', 'number ', &
+             'number ', 'm2     ', 'm2     ', 'm2     ', 'm2/s2  ', 'm2/s2  ', &
+             'm2/s2  ', 'm2/s2  ', 'm2/s2  ', 'm2/s2  ', 'number2' /)
+
+    INTEGER(iwp) ::  dots_num  = 31  !< number of timeseries defined by default
+    INTEGER(iwp) ::  dots_soil = 24  !< starting index for soil-timeseries
+    INTEGER(iwp) ::  dots_rad  = 32  !< starting index for radiation-timeseries
+
+    CHARACTER (LEN=13), DIMENSION(dots_max) :: dots_label =                    &
+          (/ 'E            ', 'E*           ', 'dt           ',                &
+             'u*           ', 'th*          ', 'umax         ',                &
+             'vmax         ', 'wmax         ', 'div_new      ',                &
+             'div_old      ', 'z_i_wpt      ', 'z_i_pt       ',                &
+             'w*           ', 'w"pt"0       ', 'w"pt"        ',                &
+             'wpt          ', 'pt(0)        ', 'pt(z_mo)     ',                &
+             'w"u"0        ', 'w"v"0        ', 'w"q"0        ',                &
+             'ol           ', 'q*           ', 'ghf_eb       ',                &
+             'shf_eb       ', 'qsws_eb      ', 'qsws_liq_eb  ',                &
+             'qsws_soil_eb ', 'qsws_veg_eb  ', 'r_a          ',                &
+             'r_s          ', 'rad_net      ', 'rad_lw_in    ',                &
+             'rad_lw_out   ', 'rad_sw_in    ', 'rad_sw_out   ',                &
+             'rrtm_aldif   ', 'rrtm_aldir   ', 'rrtm_asdif   ',                &
+             'rrtm_asdir   ',                                                  &
+             ( 'unknown      ', i9 = 1, dots_max-40 ) /)
+
+    CHARACTER (LEN=13), DIMENSION(dots_max) :: dots_unit =                     &
+          (/ 'm2/s2        ', 'm2/s2        ', 's            ',                &
+             'm/s          ', 'K            ', 'm/s          ',                &
+             'm/s          ', 'm/s          ', 's-1          ',                &
+             's-1          ', 'm            ', 'm            ',                &
+             'm/s          ', 'K m/s        ', 'K m/s        ',                &
+             'K m/s        ', 'K            ', 'K            ',                &
+             'm2/s2        ', 'm2/s2        ', 'kg m/s       ',                &
+             'm            ', 'kg/kg        ', '             ',                &
+             '             ', '             ', '             ',                &
+             '             ', 'W/m2         ', 's/m          ',                &
+             '             ', 'W/m2         ', 'W/m2         ',                &
+             'W/m2         ', 'W/m2         ', 'W/m2         ',                &
+             '             ', '             ', '             ',                &
+             '             ',                                                  &
+             ( 'unknown      ', i9 = 1, dots_max-40 ) /)
+
+    CHARACTER (LEN=9), DIMENSION(300) ::  dopr_unit = 'unknown'
+
+    CHARACTER (LEN=7), DIMENSION(0:1,100) ::  do2d_unit, do3d_unit
+
+    CHARACTER (LEN=16), DIMENSION(25) ::  prt_var_names = &
+          (/ 'pt_age          ', 'pt_dvrp_size    ', 'pt_origin_x     ', &
+             'pt_origin_y     ', 'pt_origin_z     ', 'pt_radius       ', &
+             'pt_speed_x      ', 'pt_speed_y      ', 'pt_speed_z      ', &
+             'pt_weight_factor', 'pt_x            ', 'pt_y            ', &
+             'pt_z            ', 'pt_color        ', 'pt_group        ', &
+             'pt_tailpoints   ', 'pt_tail_id      ', 'pt_density_ratio', &
+             'pt_exp_arg      ', 'pt_exp_term     ', 'not_used        ', &
+             'not_used        ', 'not_used        ', 'not_used        ', &
+             'not_used        ' /)
+
+    CHARACTER (LEN=16), DIMENSION(25) ::  prt_var_units = &
+          (/ 'seconds         ', 'meters          ', 'meters          ', &
+             'meters          ', 'meters          ', 'meters          ', &
+             'm/s             ', 'm/s             ', 'm/s             ', &
+             'factor          ', 'meters          ', 'meters          ', &
+             'meters          ', 'none            ', 'none            ', &
+             'none            ', 'none            ', 'ratio           ', &
+             'none            ', 'none            ', 'not_used        ', &
+             'not_used        ', 'not_used        ', 'not_used        ', &
+             'not_used        ' /)
+
+    CHARACTER(LEN=20), DIMENSION(11) ::  netcdf_precision = ' '
+    CHARACTER(LEN=40) ::  netcdf_data_format_string
+
+    INTEGER(iwp) ::  id_dim_prtnum, id_dim_time_pr, id_dim_time_prt, &
+                     id_dim_time_pts, id_dim_time_sp, id_dim_time_ts, id_dim_x_sp, &
+                     id_dim_y_sp, id_dim_zu_sp, id_dim_zw_sp, id_set_pr, &
+                     id_set_prt, id_set_pts, id_set_sp, id_set_ts, id_var_prtnum, &
+                     id_var_rnop_prt, id_var_time_pr, id_var_time_prt, &
+                     id_var_time_pts, id_var_time_sp, id_var_time_ts, id_var_x_sp, &
+                     id_var_y_sp, id_var_zu_sp, id_var_zw_sp, nc_stat
+
+    INTEGER(iwp), DIMENSION(0:1) ::  id_dim_time_xy, id_dim_time_xz, &
+                    id_dim_time_yz, id_dim_time_3d, id_dim_x_xy, id_dim_xu_xy, &
+                    id_dim_x_xz, id_dim_xu_xz, id_dim_x_yz, id_dim_xu_yz, &
+                    id_dim_x_3d, id_dim_xu_3d, id_dim_y_xy, id_dim_yv_xy, &
+                    id_dim_y_xz, id_dim_yv_xz, id_dim_y_yz, id_dim_yv_yz, &
+                    id_dim_y_3d, id_dim_yv_3d, id_dim_zs_xy, id_dim_zs_xz, &
+                    id_dim_zs_yz, id_dim_zs_3d, id_dim_zu_xy, id_dim_zu1_xy, &
+                    id_dim_zu_xz, id_dim_zu_yz, id_dim_zu_3d, id_dim_zw_xy, &
+                    id_dim_zw_xz, id_dim_zw_yz, id_dim_zw_3d, id_set_xy, &
+                    id_set_xz, id_set_yz, id_set_3d, id_var_ind_x_yz, &
+                    id_var_ind_y_xz, id_var_ind_z_xy, id_var_time_xy, &
+                    id_var_time_xz, id_var_time_yz, id_var_time_3d, id_var_x_xy, &
+                    id_var_xu_xy, id_var_x_xz, id_var_xu_xz, id_var_x_yz, &
+                    id_var_xu_yz, id_var_x_3d, id_var_xu_3d, id_var_y_xy, &
+                    id_var_yv_xy, id_var_y_xz, id_var_yv_xz, id_var_y_yz, &
+                    id_var_yv_yz, id_var_y_3d, id_var_yv_3d, id_var_zs_xy, &
+                    id_var_zs_xz, id_var_zs_yz, id_var_zs_3d, id_var_zusi_xy, &
+                    id_var_zusi_3d, id_var_zu_xy, id_var_zu1_xy, id_var_zu_xz, &
+                    id_var_zu_yz, id_var_zu_3d, id_var_zwwi_xy, id_var_zwwi_3d, &
+                    id_var_zw_xy, id_var_zw_xz, id_var_zw_yz, id_var_zw_3d
+
+    INTEGER ::  netcdf_data_format = 2  !< NetCDF3 64bit offset format
+    INTEGER ::  netcdf_deflate = 0      !< NetCDF compression, default: no
+                                        !< compression
+
+    INTEGER(iwp), DIMENSION(10)  ::  id_var_dospx, id_var_dospy
+    INTEGER(iwp), DIMENSION(20)  ::  id_var_prt
+    INTEGER(iwp), DIMENSION(11)  ::  nc_precision
+    INTEGER(iwp), DIMENSION(dopr_norm_num) ::  id_var_norm_dopr
+
+    INTEGER(iwp), DIMENSION(dopts_num,0:10) ::  id_var_dopts
+    INTEGER(iwp), DIMENSION(0:1,100)        ::  id_var_do2d, id_var_do3d
+    INTEGER(iwp), DIMENSION(100,0:9)        ::  id_dim_z_pr, id_var_dopr, &
+                                                id_var_z_pr
+    INTEGER(iwp), DIMENSION(dots_max,0:9)   ::  id_var_dots
+
+!
+!-- Masked output
+    CHARACTER (LEN=7), DIMENSION(max_masks,0:1,100) ::  domask_unit
+
+    LOGICAL ::  output_for_t0 = .FALSE.
+
+    INTEGER(iwp), DIMENSION(1:max_masks,0:1) ::  id_dim_time_mask, id_dim_x_mask, &
+                   id_dim_xu_mask, id_dim_y_mask, id_dim_yv_mask, id_dim_zs_mask, &
+                   id_dim_zu_mask, id_dim_zw_mask, &
+                   id_set_mask, &
+                   id_var_time_mask, id_var_x_mask, id_var_xu_mask, &
+                   id_var_y_mask, id_var_yv_mask, id_var_zs_mask, &
+                   id_var_zu_mask, id_var_zw_mask, &
+                   id_var_zusi_mask, id_var_zwwi_mask
+
+    INTEGER(iwp), DIMENSION(1:max_masks,0:1,100) ::  id_var_domask
+
+
+    PUBLIC  domask_unit, dopr_unit, dopts_num, dots_label, dots_max, dots_num, &
+            dots_rad, dots_soil, dots_unit, do2d_unit, do3d_unit, id_set_mask, &
+            id_set_pr, id_set_prt, id_set_pts, id_set_sp, id_set_ts,           &
+            id_set_xy, id_set_xz, id_set_yz, id_set_3d, id_var_domask,         &
+            id_var_dopr, id_var_dopts, id_var_dospx, id_var_dospy,             &
+            id_var_dots, id_var_do2d, id_var_do3d, id_var_norm_dopr,           &
+            id_var_time_mask, id_var_time_pr, id_var_time_pts, id_var_time_ts, &
+            id_var_time_xy, id_var_time_xz, id_var_time_yz,id_var_time_3d,     &
+            nc_stat, netcdf_data_format, netcdf_data_format_string,            &
+            netcdf_deflate, netcdf_precision, output_for_t0
+
+    SAVE
+
+    INTERFACE netcdf_create_dim
+       MODULE PROCEDURE netcdf_create_dim
+    END INTERFACE netcdf_create_dim
+
+    INTERFACE netcdf_create_file
+       MODULE PROCEDURE netcdf_create_file
+    END INTERFACE netcdf_create_file
+
+    INTERFACE netcdf_create_var
+       MODULE PROCEDURE netcdf_create_var
+    END INTERFACE netcdf_create_var
+
+    INTERFACE netcdf_define_header
+       MODULE PROCEDURE netcdf_define_header
+    END INTERFACE netcdf_define_header
+
+    INTERFACE netcdf_handle_error
+       MODULE PROCEDURE netcdf_handle_error
+    END INTERFACE netcdf_handle_error
+
+    INTERFACE netcdf_open_write_file
+       MODULE PROCEDURE netcdf_open_write_file
+    END INTERFACE netcdf_open_write_file
+
+    PUBLIC netcdf_create_file, netcdf_define_header, netcdf_handle_error,      &
+           netcdf_open_write_file
+
+ CONTAINS
+
+ SUBROUTINE netcdf_define_header( callmode, extend, av )
  
 #if defined( __netcdf )
@@ -141,6 +347,5 @@
                dt_do2d_yz, dt_do3d, mask_size, do2d_xy_time_count,             &
                do3d_time_count, domask_time_count, end_time, mask_i_global,    &
-               mask_j_global, mask_k_global, message_string, mid,              &
-               netcdf_data_format, netcdf_precision, ntdim_2d_xy,              &
+               mask_j_global, mask_k_global, message_string, mid, ntdim_2d_xy, &
                ntdim_2d_xz, ntdim_2d_yz, ntdim_3d, nz_do3d, prt_time_count,    &
                run_description_header, section, simulated_time,                &
@@ -155,13 +360,8 @@
         ONLY:  nx, ny, nz ,nzb, nzt
 
-    USE netcdf_control
-
     USE kinds
 
     USE land_surface_model_mod,                                                &
-        ONLY: land_surface, nzb_soil, nzt_soil, id_dim_zs_xy, id_dim_zs_xz,    &
-              id_dim_zs_yz, id_dim_zs_3d, id_dim_zs_mask, id_var_zs_xy,        &
-              id_var_zs_xz, id_var_zs_yz ,id_var_zs_3d, id_var_zs_mask,        &
-              nzs, zs
+        ONLY: land_surface, nzb_soil, nzt_soil, nzs, zs
 
     USE pegrid
@@ -172,9 +372,4 @@
     USE profil_parameter,                                                      &
         ONLY:  crmax, cross_profiles, dopr_index, profile_columns, profile_rows
-
-    USE radiation_model_mod,                                                   &
-        ONLY:  rad_lw_in, rad_lw_out, rad_lw_cs_hr, rad_lw_hr,                 &
-               rad_sw_in, rad_sw_out, rad_sw_cs_hr, rad_sw_hr
-
 
     USE spectrum,                                                              &
@@ -250,5 +445,5 @@
 
 !
-!-- Initializing actions (return to calling routine check_parameters afterwards)
+!-- Initializing actions
     IF ( .NOT. init_netcdf )  THEN
 !
@@ -263,5 +458,5 @@
                                          '"_"netcdf_precision(', i, ')="',   &
                                          TRIM( netcdf_precision(i) ),'"'
-             CALL message( 'define_netcdf_header', 'PA0241', 1, 2, 0, 6, 0 )
+             CALL message( 'netcdf_define_header', 'PA0241', 2, 2, 0, 6, 0 )
           ENDIF
 
@@ -277,5 +472,5 @@
                                          'netcdf_precision(', i, ')="', &
                                          TRIM( netcdf_precision(i) ),'"'
-             CALL message( 'define_netcdf_header', 'PA0242', 1, 2, 0, 6, 0 )
+             CALL message( 'netcdf_define_header', 'PA0242', 1, 2, 0, 6, 0 )
           ENDIF
 
@@ -308,5 +503,5 @@
                                   'assignment: netcdf_precision(', i, ')="',  &
                                             TRIM( netcdf_precision(i) ),'"'
-                CALL message( 'define_netcdf_header', 'PA0243', 1, 2, 0, 6, 0 ) 
+                CALL message( 'netcdf_define_header', 'PA0243', 1, 2, 0, 6, 0 )
 
           END SELECT
@@ -316,7 +511,22 @@
        ENDDO
 
+!
+!--    Check for allowed parameter range
+       IF ( netcdf_deflate < 0  .OR.  netcdf_deflate > 9 )  THEN
+          WRITE ( message_string, '(A,I3,A)' ) 'netcdf_deflate out of ' //     &
+                                      'range & given value: ', netcdf_deflate, &
+                                      ', allowed range: 0-9'
+          CALL message( 'netcdf_define_header', 'PA0355', 2, 2, 0, 6, 0 )
+       ENDIF
+!
+!--    Data compression does not work with parallel NetCDF/HDF5
+       IF ( netcdf_deflate > 0  .AND.  netcdf_data_format /= 3 )  THEN
+          message_string = 'netcdf_deflate reset to 0'
+          CALL message( 'netcdf_define_header', 'PA0356', 0, 1, 0, 6, 0 )
+
+          netcdf_deflate = 0
+       ENDIF
+
        init_netcdf = .TRUE.
-
-       RETURN
 
     ENDIF
@@ -352,6 +562,6 @@
 !--       Define some global attributes of the dataset
           nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), NF90_GLOBAL, &
-               'Conventions', 'COARDS' )
-          CALL handle_netcdf_error( 'netcdf', 464 )
+                                  'Conventions', 'COARDS' )
+          CALL netcdf_handle_error( 'netcdf_define_header', 464 )
 
           IF ( av == 0 )  THEN
@@ -364,117 +574,71 @@
                                   TRIM( run_description_header ) //    &
                                   TRIM( time_average_text ) )
-          CALL handle_netcdf_error( 'netcdf', 465 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 465 )
           IF ( av == 1 )  THEN
              WRITE ( time_average_text,'(F7.1,'' s avg'')' )  averaging_interval
              nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), NF90_GLOBAL, &
                                      'time_avg', TRIM( time_average_text ) )
-             CALL handle_netcdf_error( 'netcdf', 466 )
+             CALL netcdf_handle_error( 'netcdf_define_header', 466 )
           ENDIF
 
 !
 !--       Define time coordinate for volume data (unlimited dimension)
-          nc_stat = NF90_DEF_DIM( id_set_mask(mid,av), 'time', NF90_UNLIMITED, &
-                                  id_dim_time_mask(mid,av) )
-          CALL handle_netcdf_error( 'netcdf', 467 )
-
-          nc_stat = NF90_DEF_VAR( id_set_mask(mid,av), 'time', NF90_DOUBLE, &
-                                  id_dim_time_mask(mid,av), &
-                                  id_var_time_mask(mid,av) )
-          CALL handle_netcdf_error( 'netcdf', 468 )
-
-          nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), &
-                                  id_var_time_mask(mid,av), 'units', &
-                                  'seconds')
-          CALL handle_netcdf_error( 'netcdf', 469 )
-
+          CALL netcdf_create_dim( id_set_mask(mid,av), 'time', NF90_UNLIMITED, &
+                                  id_dim_time_mask(mid,av), 467 )
+          CALL netcdf_create_var( id_set_mask(mid,av),                         &
+                                  (/ id_dim_time_mask(mid,av) /), 'time',      &
+                                  NF90_DOUBLE, id_var_time_mask(mid,av),       &
+                                 'seconds', '', 468, 469, 000 )
 !
 !--       Define spatial dimensions and coordinates:
 !--       Define vertical coordinate grid (zu grid)
-          nc_stat = NF90_DEF_DIM( id_set_mask(mid,av), 'zu_3d', &
-                                  mask_size(mid,3), id_dim_zu_mask(mid,av) )
-          CALL handle_netcdf_error( 'netcdf', 470 )
-
-          nc_stat = NF90_DEF_VAR( id_set_mask(mid,av), 'zu_3d', NF90_DOUBLE, &
-                                  id_dim_zu_mask(mid,av), &
-                                  id_var_zu_mask(mid,av) )
-          CALL handle_netcdf_error( 'netcdf', 471 )
-
-          nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), id_var_zu_mask(mid,av), &
-                                  'units', 'meters' )
-          CALL handle_netcdf_error( 'netcdf', 472 )
-
+          CALL netcdf_create_dim( id_set_mask(mid,av), 'zu_3d',                &
+                                  mask_size(mid,3), id_dim_zu_mask(mid,av),    &
+                                  470 )
+          CALL netcdf_create_var( id_set_mask(mid,av),                         &
+                                  (/ id_dim_zu_mask(mid,av) /), 'zu_3d',       &
+                                  NF90_DOUBLE, id_var_zu_mask(mid,av),         &
+                                  'meters', '', 471, 472, 000 )
 !
 !--       Define vertical coordinate grid (zw grid)
-          nc_stat = NF90_DEF_DIM( id_set_mask(mid,av), 'zw_3d', &
-                                  mask_size(mid,3), id_dim_zw_mask(mid,av) )
-          CALL handle_netcdf_error( 'netcdf', 473 )
-
-          nc_stat = NF90_DEF_VAR( id_set_mask(mid,av), 'zw_3d', NF90_DOUBLE, &
-                                  id_dim_zw_mask(mid,av), &
-                                  id_var_zw_mask(mid,av) )
-          CALL handle_netcdf_error( 'netcdf', 474 )
-
-          nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), id_var_zw_mask(mid,av), &
-                                  'units', 'meters' )
-          CALL handle_netcdf_error( 'netcdf', 475 )
-
+          CALL netcdf_create_dim( id_set_mask(mid,av), 'zw_3d',                &
+                                  mask_size(mid,3), id_dim_zw_mask(mid,av),    &
+                                  473 )
+          CALL netcdf_create_var( id_set_mask(mid,av),                         &
+                                  (/ id_dim_zw_mask(mid,av) /), 'zw_3d',       &
+                                  NF90_DOUBLE, id_var_zw_mask(mid,av),         &
+                                 'meters', '', 474, 475, 000 )
 !
 !--       Define x-axis (for scalar position)
-          nc_stat = NF90_DEF_DIM( id_set_mask(mid,av), 'x', &
-                                  mask_size(mid,1), id_dim_x_mask(mid,av) )
-          CALL handle_netcdf_error( 'netcdf', 476 )
-
-          nc_stat = NF90_DEF_VAR( id_set_mask(mid,av), 'x', NF90_DOUBLE, &
-                                  id_dim_x_mask(mid,av), id_var_x_mask(mid,av) )
-          CALL handle_netcdf_error( 'netcdf', 477 )
-
-          nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), id_var_x_mask(mid,av), &
-                                  'units', 'meters' )
-          CALL handle_netcdf_error( 'netcdf', 478 )
-
+          CALL netcdf_create_dim( id_set_mask(mid,av), 'x', mask_size(mid,1),  &
+                                  id_dim_x_mask(mid,av), 476 )
+          CALL netcdf_create_var( id_set_mask(mid,av),                         &
+                                  (/ id_dim_x_mask(mid,av) /), 'x',            &
+                                  NF90_DOUBLE, id_var_x_mask(mid,av),          &
+                                  'meters', '', 477, 478, 000 )
 !
 !--       Define x-axis (for u position)
-          nc_stat = NF90_DEF_DIM( id_set_mask(mid,av), 'xu', &
-                                  mask_size(mid,1), id_dim_xu_mask(mid,av) )
-          CALL handle_netcdf_error( 'netcdf', 479 )
-
-          nc_stat = NF90_DEF_VAR( id_set_mask(mid,av), 'xu', NF90_DOUBLE, &
-                                  id_dim_xu_mask(mid,av), &
-                                  id_var_xu_mask(mid,av) )
-          CALL handle_netcdf_error( 'netcdf', 480 )
-
-          nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), id_var_xu_mask(mid,av), &
-                                  'units', 'meters' )
-          CALL handle_netcdf_error( 'netcdf', 481 )
-
+          CALL netcdf_create_dim( id_set_mask(mid,av), 'xu', mask_size(mid,1), &
+                                  id_dim_xu_mask(mid,av), 479 )
+          CALL netcdf_create_var( id_set_mask(mid,av),                         &
+                                  (/ id_dim_xu_mask(mid,av) /), 'xu',          &
+                                  NF90_DOUBLE, id_var_xu_mask(mid,av),         &
+                                  'meters', '', 480, 481, 000 )
 !
 !--       Define y-axis (for scalar position)
-          nc_stat = NF90_DEF_DIM( id_set_mask(mid,av), 'y', &
-                                  mask_size(mid,2), id_dim_y_mask(mid,av) )
-          CALL handle_netcdf_error( 'netcdf', 482 )
-
-          nc_stat = NF90_DEF_VAR( id_set_mask(mid,av), 'y', NF90_DOUBLE, &
-                                  id_dim_y_mask(mid,av), id_var_y_mask(mid,av) )
-          CALL handle_netcdf_error( 'netcdf', 483 )
-
-          nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), id_var_y_mask(mid,av), &
-                                  'units', 'meters' )
-          CALL handle_netcdf_error( 'netcdf', 484 )
-
+          CALL netcdf_create_dim( id_set_mask(mid,av), 'y', mask_size(mid,2),  &
+                                  id_dim_y_mask(mid,av), 482 )
+          CALL netcdf_create_var( id_set_mask(mid,av),                         &
+                                  (/ id_dim_y_mask(mid,av) /), 'y',            &
+                                  NF90_DOUBLE, id_var_y_mask(mid,av),          &
+                                  'meters', '', 483, 484, 000 )
 !
 !--       Define y-axis (for v position)
-          nc_stat = NF90_DEF_DIM( id_set_mask(mid,av), 'yv', &
-                                  mask_size(mid,2), id_dim_yv_mask(mid,av) )
-          CALL handle_netcdf_error( 'netcdf', 485 )
-
-          nc_stat = NF90_DEF_VAR( id_set_mask(mid,av), 'yv', NF90_DOUBLE, &
-                                  id_dim_yv_mask(mid,av), &
-                                  id_var_yv_mask(mid,av) )
-          CALL handle_netcdf_error( 'netcdf', 486 )
-
-          nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), id_var_yv_mask(mid,av), &
-                                  'units', 'meters' )
-          CALL handle_netcdf_error( 'netcdf', 487 )
-
+          CALL netcdf_create_dim( id_set_mask(mid,av), 'yv', mask_size(mid,2), &
+                                  id_dim_yv_mask(mid,av), 485 )
+          CALL netcdf_create_var( id_set_mask(mid,av),                         &
+                                  (/ id_dim_yv_mask(mid,av) /),                &
+                                  'yv', NF90_DOUBLE, id_var_yv_mask(mid,av),   &
+                                  'meters', '', 486, 487, 000 )
 !
 !--       In case of non-flat topography define 2d-arrays containing the height
@@ -483,38 +647,18 @@
 !
 !--          Define zusi = zu(nzb_s_inner)
-             nc_stat = NF90_DEF_VAR( id_set_mask(mid,av), 'zusi', NF90_DOUBLE, &
-                                     (/ id_dim_x_mask(mid,av),    &
-                                        id_dim_y_mask(mid,av) /), &
-                                     id_var_zusi_mask(mid,av) )
-             CALL handle_netcdf_error( 'netcdf', 488 )
-             
-             nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), &
-                                     id_var_zusi_mask(mid,av), &
-                                     'units', 'meters' )
-             CALL handle_netcdf_error( 'netcdf', 489 )
-             
-             nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), &
-                                     id_var_zusi_mask(mid,av), &
-                                     'long_name', 'zu(nzb_s_inner)' )
-             CALL handle_netcdf_error( 'netcdf', 490 )
-
+             CALL netcdf_create_var( id_set_mask(mid,av),                      &
+                                     (/ id_dim_x_mask(mid,av),                 &
+                                        id_dim_y_mask(mid,av) /), 'zusi',      &
+                                     NF90_DOUBLE, id_var_zusi_mask(mid,av),    &
+                                     'meters', 'zu(nzb_s_inner)', 488, 489,    &
+                                     490 )
 !             
 !--          Define zwwi = zw(nzb_w_inner)
-             nc_stat = NF90_DEF_VAR( id_set_mask(mid,av), 'zwwi', NF90_DOUBLE, &
-                                     (/ id_dim_x_mask(mid,av),    &
-                                        id_dim_y_mask(mid,av) /), &
-                                     id_var_zwwi_mask(mid,av) )
-             CALL handle_netcdf_error( 'netcdf', 491 )
-             
-             nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), &
-                                     id_var_zwwi_mask(mid,av), &
-                                     'units', 'meters' )
-             CALL handle_netcdf_error( 'netcdf', 492 )
-             
-             nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), &
-                                     id_var_zwwi_mask(mid,av), &
-                                     'long_name', 'zw(nzb_w_inner)' )
-             CALL handle_netcdf_error( 'netcdf', 493 )
-
+             CALL netcdf_create_var( id_set_mask(mid,av),                      &
+                                     (/ id_dim_x_mask(mid,av),                 &
+                                        id_dim_y_mask(mid,av) /), 'zwwi',      &
+                                     NF90_DOUBLE, id_var_zwwi_mask(mid,av),    &
+                                     'meters', 'zw(nzb_w_inner)', 491, 492,    &
+                                     493 )
           ENDIF             
  
@@ -522,17 +666,11 @@
 !
 !--          Define vertical coordinate grid (zw grid)
-             nc_stat = NF90_DEF_DIM( id_set_mask(mid,av), 'zs_3d', &
-                                     mask_size(mid,3), id_dim_zs_mask(mid,av) )
-             CALL handle_netcdf_error( 'netcdf', 536 )
-
-             nc_stat = NF90_DEF_VAR( id_set_mask(mid,av), 'zs_3d', NF90_DOUBLE, &
-                                     id_dim_zs_mask(mid,av), &
-                                     id_var_zs_mask(mid,av) )
-             CALL handle_netcdf_error( 'netcdf', 536 )
-
-             nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), id_var_zs_mask(mid,av), &
-                                  'units', 'meters' )
-             CALL handle_netcdf_error( 'netcdf', 537 )
-
+             CALL netcdf_create_dim( id_set_mask(mid,av), 'zs_3d',             &
+                                     mask_size(mid,3), id_dim_zs_mask(mid,av), &
+                                     536 )
+             CALL netcdf_create_var( id_set_mask(mid,av),                      &
+                                     (/ id_dim_zs_mask(mid,av) /), 'zs_3d',    &
+                                     NF90_DOUBLE, id_var_zs_mask(mid,av),      &
+                                     'meters', '', 537, 555, 000 )
           ENDIF
 
@@ -591,6 +729,13 @@
 !
 !--                Check for user-defined quantities
-                   CALL user_define_netcdf_grid( domask(mid,av,i), found, &
+                   CALL user_define_netcdf_grid( domask(mid,av,i), found,      &
                                                  grid_x, grid_y, grid_z )
+
+                   IF ( .NOT. found )  THEN
+                      WRITE ( message_string, * ) 'no grid defined for',       &
+                           ' variable ', TRIM( domask(mid,av,i) )
+                      CALL message( 'define_netcdf_header', 'PA0244', 0, 1, 0, &
+                                    6, 0 )
+                   ENDIF
 
              END SELECT
@@ -620,33 +765,12 @@
 !
 !--          Define the grid
-             nc_stat = NF90_DEF_VAR( id_set_mask(mid,av), domask(mid,av,i), &
-                                     nc_precision(11),            &
-                                     (/ id_x, id_y, id_z,                   &
-                                        id_dim_time_mask(mid,av) /),        &
-                                     id_var_domask(mid,av,i) )
-
-             IF ( .NOT. found )  THEN
-                WRITE ( message_string, * ) 'no grid defined for', &
-                     ' variable ', TRIM( domask(mid,av,i) )
-                CALL message( 'define_netcdf_header', 'PA0244', 0, 1, 0, 6, 0 )
-             ENDIF
+             CALL netcdf_create_var( id_set_mask(mid,av), (/ id_x, id_y, id_z, &
+                                     id_dim_time_mask(mid,av) /),              &
+                                     domask(mid,av,i), nc_precision(11),       &
+                                     id_var_domask(mid,av,i),                  &
+                                     TRIM( domask_unit(mid,av,i) ),            &
+                                     domask(mid,av,i), 494, 495, 496 )
 
              var_list = TRIM( var_list ) // TRIM( domask(mid,av,i) ) // ';'
-
-             CALL handle_netcdf_error( 'netcdf', 494 )
-!
-!--          Store the 'real' name of the variable (with *, for example)
-!--          in the long_name attribute. This is evaluated by Ferret,
-!--          for example.
-             nc_stat = NF90_PUT_ATT( id_set_mask(mid,av),     &
-                                     id_var_domask(mid,av,i), &
-                                     'long_name', domask(mid,av,i) )
-             CALL handle_netcdf_error( 'netcdf', 495 )
-!
-!--          Define the variable's unit
-             nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), &
-                                     id_var_domask(mid,av,i), &
-                                     'units', TRIM( domask_unit(mid,av,i) ) )
-             CALL handle_netcdf_error( 'netcdf', 496 )
 
              i = i + 1
@@ -663,10 +787,10 @@
           nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), NF90_GLOBAL, &
                                   'VAR_LIST', var_list )
-          CALL handle_netcdf_error( 'netcdf', 497 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 497 )
 
 !
 !--       Leave netCDF define mode
           nc_stat = NF90_ENDDEF( id_set_mask(mid,av) )
-          CALL handle_netcdf_error( 'netcdf', 498 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 498 )
 
 !
@@ -679,5 +803,5 @@
                                   netcdf_data, start = (/ 1 /),               &
                                   count = (/ mask_size(mid,1) /) )
-          CALL handle_netcdf_error( 'netcdf', 499 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 499 )
 
           netcdf_data = mask_i_global(mid,:mask_size(mid,1)) * dx
@@ -686,5 +810,5 @@
                                   netcdf_data, start = (/ 1 /),               &
                                   count = (/ mask_size(mid,1) /) )
-          CALL handle_netcdf_error( 'netcdf', 500 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 500 )
 
           DEALLOCATE( netcdf_data )
@@ -699,5 +823,5 @@
                                   netcdf_data, start = (/ 1 /),               &
                                   count = (/ mask_size(mid,2) /))
-          CALL handle_netcdf_error( 'netcdf', 501 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 501 )
 
           netcdf_data = mask_j_global(mid,:mask_size(mid,2)) * dy
@@ -706,5 +830,5 @@
                                   netcdf_data, start = (/ 1 /),    &
                                   count = (/ mask_size(mid,2) /))
-          CALL handle_netcdf_error( 'netcdf', 502 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 502 )
 
           DEALLOCATE( netcdf_data )
@@ -719,5 +843,5 @@
                                   netcdf_data, start = (/ 1 /), &
                                   count = (/ mask_size(mid,3) /) )
-          CALL handle_netcdf_error( 'netcdf', 503 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 503 )
 
           netcdf_data = zw( mask_k_global(mid,:mask_size(mid,3)) )
@@ -726,5 +850,5 @@
                                   netcdf_data, start = (/ 1 /), &
                                   count = (/ mask_size(mid,3) /) )
-          CALL handle_netcdf_error( 'netcdf', 504 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 504 )
 
           DEALLOCATE( netcdf_data )
@@ -744,5 +868,5 @@
                                      count = (/ mask_size(mid,1), &
                                                 mask_size(mid,2) /) )
-             CALL handle_netcdf_error( 'netcdf', 505 )
+             CALL netcdf_handle_error( 'netcdf_define_header', 505 )
 
              netcdf_data_2d = zw_w_inner( mask_i_global(mid,:mask_size(mid,1)),&
@@ -755,5 +879,5 @@
                                      count = (/ mask_size(mid,1), &
                                                 mask_size(mid,2) /) )
-             CALL handle_netcdf_error( 'netcdf', 506 )
+             CALL netcdf_handle_error( 'netcdf_define_header', 506 )
 
              DEALLOCATE( netcdf_data_2d )
@@ -772,5 +896,5 @@
                                      netcdf_data, start = (/ 1 /), &
                                      count = (/ mask_size(mid,3) /) )
-             CALL handle_netcdf_error( 'netcdf', 538 )
+             CALL netcdf_handle_error( 'netcdf_define_header', 538 )
 
              DEALLOCATE( netcdf_data )
@@ -804,5 +928,5 @@
           nc_stat = NF90_GET_ATT( id_set_mask(mid,av), NF90_GLOBAL, 'VAR_LIST',&
                                   var_list_old )
-          CALL handle_netcdf_error( 'netcdf', 507 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 507 )
 
           var_list = ';'
@@ -833,10 +957,10 @@
           nc_stat = NF90_INQ_VARID( id_set_mask(mid,av), 'zu_3d', &
                                     id_var_zu_mask(mid,av) )
-          CALL handle_netcdf_error( 'netcdf', 508 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 508 )
 
           nc_stat = NF90_INQUIRE_VARIABLE( id_set_mask(mid,av),     &
                                            id_var_zu_mask(mid,av),  &
                                            dimids = id_dim_zu_mask_old )
-          CALL handle_netcdf_error( 'netcdf', 509 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 509 )
           id_dim_zu_mask(mid,av) = id_dim_zu_mask_old(1)
 
@@ -844,5 +968,5 @@
                                             id_dim_zu_mask(mid,av), &
                                             len = nz_old )
-          CALL handle_netcdf_error( 'netcdf', 510 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 510 )
 
           IF ( mask_size(mid,3) /= nz_old )  THEN
@@ -864,10 +988,10 @@
           nc_stat = NF90_INQ_VARID( id_set_mask(mid,av), 'time', &
                                     id_var_time_mask(mid,av) )
-          CALL handle_netcdf_error( 'netcdf', 511 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 511 )
 
           nc_stat = NF90_INQUIRE_VARIABLE( id_set_mask(mid,av), &
                                            id_var_time_mask(mid,av), &
                                            dimids = id_dim_time_old )
-          CALL handle_netcdf_error( 'netcdf', 512 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 512 )
           id_dim_time_mask(mid,av) = id_dim_time_old(1)
 
@@ -875,5 +999,5 @@
                                             id_dim_time_mask(mid,av), &
                                             len = domask_time_count(mid,av) )
-          CALL handle_netcdf_error( 'netcdf', 513 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 513 )
 
           nc_stat = NF90_GET_VAR( id_set_mask(mid,av), &
@@ -882,5 +1006,5 @@
                                   start = (/ domask_time_count(mid,av) /), &
                                   count = (/ 1 /) )
-          CALL handle_netcdf_error( 'netcdf', 514 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 514 )
 
           IF ( last_time_coordinate(1) >= simulated_time )  THEN
@@ -904,5 +1028,5 @@
                                        TRIM( domask(mid,av,i) ), &
                                        id_var_domask(mid,av,i) )
-             CALL handle_netcdf_error( 'netcdf', 515 )
+             CALL netcdf_handle_error( 'netcdf_define_header', 515 )
              i = i + 1
           ENDDO
@@ -922,11 +1046,11 @@
           ENDIF
           nc_stat = NF90_REDEF( id_set_mask(mid,av) )
-          CALL handle_netcdf_error( 'netcdf', 516 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 516 )
           nc_stat = NF90_PUT_ATT( id_set_mask(mid,av), NF90_GLOBAL, 'title', &
                                   TRIM( run_description_header ) //    &
                                   TRIM( time_average_text ) )
-          CALL handle_netcdf_error( 'netcdf', 517 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 517 )
           nc_stat = NF90_ENDDEF( id_set_mask(mid,av) )
-          CALL handle_netcdf_error( 'netcdf', 518 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 518 )
           WRITE ( message_string, * ) 'netCDF file for ', TRIM( var ), &
                ' data for mask', mid, ' from previous run found.', &
@@ -944,5 +1068,5 @@
           nc_stat = NF90_PUT_ATT( id_set_3d(av), NF90_GLOBAL, 'Conventions', &
                                   'COARDS' )
-          CALL handle_netcdf_error( 'netcdf', 62 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 62 )
 
           IF ( av == 0 )  THEN
@@ -955,10 +1079,10 @@
                                   TRIM( run_description_header ) //    &
                                   TRIM( time_average_text ) )
-          CALL handle_netcdf_error( 'netcdf', 63 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 63 )
           IF ( av == 1 )  THEN
              WRITE ( time_average_text,'(F7.1,'' s avg'')' )  averaging_interval
              nc_stat = NF90_PUT_ATT( id_set_3d(av), NF90_GLOBAL, 'time_avg', &
                                      TRIM( time_average_text ) )
-             CALL handle_netcdf_error( 'netcdf', 63 )
+             CALL netcdf_handle_error( 'netcdf_define_header', 63 )
           ENDIF
 
@@ -968,102 +1092,57 @@
 !--       the performance drops significantly.
           IF ( netcdf_data_format < 5 )  THEN
-             nc_stat = NF90_DEF_DIM( id_set_3d(av), 'time', NF90_UNLIMITED, &
-                                     id_dim_time_3d(av) )
-             CALL handle_netcdf_error( 'netcdf', 64 )
+             CALL netcdf_create_dim( id_set_3d(av), 'time', NF90_UNLIMITED,    &
+                                     id_dim_time_3d(av), 64 )
           ELSE
-             nc_stat = NF90_DEF_DIM( id_set_3d(av), 'time', ntdim_3d(av), &
-                                     id_dim_time_3d(av) )
-             CALL handle_netcdf_error( 'netcdf', 523 )
-          ENDIF
-
-          nc_stat = NF90_DEF_VAR( id_set_3d(av), 'time', NF90_DOUBLE, &
-                                  id_dim_time_3d(av), id_var_time_3d(av) )
-          CALL handle_netcdf_error( 'netcdf', 65 )
-
-          nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_time_3d(av), 'units', &
-                                  'seconds')
-          CALL handle_netcdf_error( 'netcdf', 66 )
-
+             CALL netcdf_create_dim( id_set_3d(av), 'time', ntdim_3d(av),      &
+                                     id_dim_time_3d(av), 523 )
+          ENDIF
+
+          CALL netcdf_create_var( id_set_3d(av), (/ id_dim_time_3d(av) /),     &
+                                  'time', NF90_DOUBLE, id_var_time_3d(av),     &
+                                  'seconds', '', 65, 66, 00 )
 !
 !--       Define spatial dimensions and coordinates:
 !--       Define vertical coordinate grid (zu grid)
-          nc_stat = NF90_DEF_DIM( id_set_3d(av), 'zu_3d', nz_do3d-nzb+1, &
-                                  id_dim_zu_3d(av) )
-          CALL handle_netcdf_error( 'netcdf', 67 )
-
-          nc_stat = NF90_DEF_VAR( id_set_3d(av), 'zu_3d', NF90_DOUBLE, &
-                                  id_dim_zu_3d(av), id_var_zu_3d(av) )
-          CALL handle_netcdf_error( 'netcdf', 68 )
-
-          nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_zu_3d(av), 'units', &
-                                  'meters' )
-          CALL handle_netcdf_error( 'netcdf', 69 )
-
+          CALL netcdf_create_dim( id_set_3d(av), 'zu_3d', nz_do3d-nzb+1,       &
+                                  id_dim_zu_3d(av), 67 )
+          CALL netcdf_create_var( id_set_3d(av), (/ id_dim_zu_3d(av) /),       &
+                                  'zu_3d', NF90_DOUBLE, id_var_zu_3d(av),      &
+                                  'meters', '', 68, 69, 00 )
 !
 !--       Define vertical coordinate grid (zw grid)
-          nc_stat = NF90_DEF_DIM( id_set_3d(av), 'zw_3d', nz_do3d-nzb+1, &
-                                  id_dim_zw_3d(av) )
-          CALL handle_netcdf_error( 'netcdf', 70 )
-
-          nc_stat = NF90_DEF_VAR( id_set_3d(av), 'zw_3d', NF90_DOUBLE, &
-                                  id_dim_zw_3d(av), id_var_zw_3d(av) )
-          CALL handle_netcdf_error( 'netcdf', 71 )
-
-          nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_zw_3d(av), 'units', &
-                                  'meters' )
-          CALL handle_netcdf_error( 'netcdf', 72 )
-
+          CALL netcdf_create_dim( id_set_3d(av), 'zw_3d', nz_do3d-nzb+1,       &
+                                  id_dim_zw_3d(av), 70 )
+          CALL netcdf_create_var( id_set_3d(av), (/ id_dim_zw_3d(av) /),       &
+                                  'zw_3d', NF90_DOUBLE, id_var_zw_3d(av),      &
+                                  'meters', '', 71, 72, 00 )
 !
 !--       Define x-axis (for scalar position)
-          nc_stat = NF90_DEF_DIM( id_set_3d(av), 'x', nx+2, id_dim_x_3d(av) )
-          CALL handle_netcdf_error( 'netcdf', 73 )
-
-          nc_stat = NF90_DEF_VAR( id_set_3d(av), 'x', NF90_DOUBLE, &
-                                  id_dim_x_3d(av), id_var_x_3d(av) )
-          CALL handle_netcdf_error( 'netcdf', 74 )
-
-          nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_x_3d(av), 'units', &
-                                  'meters' )
-          CALL handle_netcdf_error( 'netcdf', 75 )
-
+          CALL netcdf_create_dim( id_set_3d(av), 'x', nx+2, id_dim_x_3d(av),   &
+                                  73 )
+          CALL netcdf_create_var( id_set_3d(av), (/ id_dim_x_3d(av) /), 'x',   &
+                                  NF90_DOUBLE, id_var_x_3d(av), 'meters', '',  &
+                                  74, 75, 00 )
 !
 !--       Define x-axis (for u position)
-          nc_stat = NF90_DEF_DIM( id_set_3d(av), 'xu', nx+2, id_dim_xu_3d(av) )
-          CALL handle_netcdf_error( 'netcdf', 358 )
-
-          nc_stat = NF90_DEF_VAR( id_set_3d(av), 'xu', NF90_DOUBLE, &
-                                  id_dim_xu_3d(av), id_var_xu_3d(av) )
-          CALL handle_netcdf_error( 'netcdf', 359 )
-
-          nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_xu_3d(av), 'units', &
-                                  'meters' )
-          CALL handle_netcdf_error( 'netcdf', 360 )
-
+          CALL netcdf_create_dim( id_set_3d(av), 'xu', nx+2, id_dim_xu_3d(av), &
+                                  358 )
+          CALL netcdf_create_var( id_set_3d(av), (/ id_dim_xu_3d(av) /), 'xu', &
+                                  NF90_DOUBLE, id_var_xu_3d(av), 'meters', '', &
+                                  359, 360, 000 )
 !
 !--       Define y-axis (for scalar position)
-          nc_stat = NF90_DEF_DIM( id_set_3d(av), 'y', ny+2, id_dim_y_3d(av) )
-          CALL handle_netcdf_error( 'netcdf', 76 )
-
-          nc_stat = NF90_DEF_VAR( id_set_3d(av), 'y', NF90_DOUBLE, &
-                                  id_dim_y_3d(av), id_var_y_3d(av) )
-          CALL handle_netcdf_error( 'netcdf', 77 )
-
-          nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_y_3d(av), 'units', &
-                                  'meters' )
-          CALL handle_netcdf_error( 'netcdf', 78 )
-
+          CALL netcdf_create_dim( id_set_3d(av), 'y', ny+2, id_dim_y_3d(av),   &
+                                  76 )
+          CALL netcdf_create_var( id_set_3d(av), (/ id_dim_y_3d(av) /), 'y',   &
+                                  NF90_DOUBLE, id_var_y_3d(av), 'meters', '',  &
+                                  77, 78, 00 )
 !
 !--       Define y-axis (for v position)
-          nc_stat = NF90_DEF_DIM( id_set_3d(av), 'yv', ny+2, id_dim_yv_3d(av) )
-          CALL handle_netcdf_error( 'netcdf', 361 )
-
-          nc_stat = NF90_DEF_VAR( id_set_3d(av), 'yv', NF90_DOUBLE, &
-                                  id_dim_yv_3d(av), id_var_yv_3d(av) )
-          CALL handle_netcdf_error( 'netcdf', 362 )
-
-          nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_yv_3d(av), 'units', &
-                                  'meters' )
-          CALL handle_netcdf_error( 'netcdf', 363 )
-
+          CALL netcdf_create_dim( id_set_3d(av), 'yv', ny+2, id_dim_yv_3d(av), &
+                                  361 )
+          CALL netcdf_create_var( id_set_3d(av), (/ id_dim_yv_3d(av) /), 'yv', &
+                                  NF90_DOUBLE, id_var_yv_3d(av), 'meters', '', &
+                                  362, 363, 000 )
 !
 !--       In case of non-flat topography define 2d-arrays containing the height
@@ -1072,31 +1151,14 @@
 !
 !--          Define zusi = zu(nzb_s_inner)
-             nc_stat = NF90_DEF_VAR( id_set_3d(av), 'zusi', NF90_DOUBLE,     &
-                                     (/ id_dim_x_3d(av), id_dim_y_3d(av) /), &
-                                     id_var_zusi_3d(av) )
-             CALL handle_netcdf_error( 'netcdf', 413 )
-             
-             nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_zusi_3d(av), &
-                                     'units', 'meters' )
-             CALL handle_netcdf_error( 'netcdf', 414 )
-             
-             nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_zusi_3d(av), &
-                                     'long_name', 'zu(nzb_s_inner)' )
-             CALL handle_netcdf_error( 'netcdf', 415 )
-
+             CALL netcdf_create_var( id_set_3d(av), (/ id_dim_x_3d(av),        &
+                                     id_dim_y_3d(av) /), 'zusi', NF90_DOUBLE,  &
+                                     id_var_zusi_3d(av), 'meters',             &
+                                     'zu(nzb_s_inner)', 413, 414, 415 )
 !             
 !--          Define zwwi = zw(nzb_w_inner)
-             nc_stat = NF90_DEF_VAR( id_set_3d(av), 'zwwi', NF90_DOUBLE,     &
-                                     (/ id_dim_x_3d(av), id_dim_y_3d(av) /), &
-                                     id_var_zwwi_3d(av) )
-             CALL handle_netcdf_error( 'netcdf', 416 )
-             
-             nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_zwwi_3d(av), &
-                                     'units', 'meters' )
-             CALL handle_netcdf_error( 'netcdf', 417 )
-             
-             nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_zwwi_3d(av), &
-                                     'long_name', 'zw(nzb_w_inner)' )
-             CALL handle_netcdf_error( 'netcdf', 418 )
+             CALL netcdf_create_var( id_set_3d(av), (/ id_dim_x_3d(av),        &
+                                     id_dim_y_3d(av) /), 'zwwi', NF90_DOUBLE,  &
+                                     id_var_zwwi_3d(av), 'meters',             &
+                                     'zw(nzb_w_inner)', 416, 417, 418 )
 
           ENDIF             
@@ -1105,15 +1167,9 @@
 !
 !--          Define vertical coordinate grid (zs grid)
-             nc_stat = NF90_DEF_DIM( id_set_3d(av), 'zs_3d', nzt_soil-nzb_soil+1, &
-                                     id_dim_zs_3d(av) )
-             CALL handle_netcdf_error( 'netcdf', 70 )
-
-             nc_stat = NF90_DEF_VAR( id_set_3d(av), 'zs_3d', NF90_DOUBLE, &
-                                     id_dim_zs_3d(av), id_var_zs_3d(av) )
-             CALL handle_netcdf_error( 'netcdf', 71 )
-
-             nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_zs_3d(av), 'units', &
-                                     'meters' )
-             CALL handle_netcdf_error( 'netcdf', 72 )
+             CALL netcdf_create_dim( id_set_3d(av), 'zs_3d',                   &
+                                     nzt_soil-nzb_soil+1, id_dim_zs_3d(av), 70 )
+             CALL netcdf_create_var( id_set_3d(av), (/ id_dim_zs_3d(av) /),    &
+                                     'zs_3d', NF90_DOUBLE, id_var_zs_3d(av),   &
+                                     'meters', '', 71, 72, 00 )
 
           ENDIF
@@ -1175,4 +1231,11 @@
                    CALL user_define_netcdf_grid( do3d(av,i), found, grid_x, &
                                                  grid_y, grid_z )
+
+                   IF ( .NOT. found )  THEN
+                      WRITE ( message_string, * ) 'no grid defined for varia', &
+                                                  'ble ', TRIM( do3d(av,i) )
+                      CALL message( 'define_netcdf_header', 'PA0244', 0, 1, 0, &
+                                    6, 0 )
+                   ENDIF
 
              END SELECT
@@ -1202,30 +1265,9 @@
 !
 !--          Define the grid
-             nc_stat = NF90_DEF_VAR( id_set_3d(av), do3d(av,i),                &
-                                     nc_precision(4),                          &
-                                   (/ id_x, id_y, id_z, id_dim_time_3d(av) /), &
-                                     id_var_do3d(av,i) )
-
-             IF ( .NOT. found )  THEN
-                WRITE ( message_string, * ) 'no grid defined for', &
-                                            ' variable ', TRIM( do3d(av,i) )
-                CALL message( 'define_netcdf_header', 'PA0244', 0, 1, 0, 6, 0 )
-             ENDIF
-
-             var_list = TRIM( var_list ) // TRIM( do3d(av,i) ) // ';'
-
-             CALL handle_netcdf_error( 'netcdf', 79 )
-!
-!--          Store the 'real' name of the variable (with *, for example)
-!--          in the long_name attribute. This is evaluated by Ferret,
-!--          for example.
-             nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_do3d(av,i), &
-                                     'long_name', do3d(av,i) )
-             CALL handle_netcdf_error( 'netcdf', 80 )
-!
-!--          Define the variable's unit
-             nc_stat = NF90_PUT_ATT( id_set_3d(av), id_var_do3d(av,i), &
-                                     'units', TRIM( do3d_unit(av,i) ) )
-             CALL handle_netcdf_error( 'netcdf', 357 )
+             CALL netcdf_create_var( id_set_3d(av),(/ id_x, id_y, id_z,        &
+                                     id_dim_time_3d(av) /), do3d(av,i),        &
+                                     nc_precision(4), id_var_do3d(av,i),       &
+                                     TRIM( do3d_unit(av,i) ), do3d(av,i), 79,  &
+                                     80, 357 )
 #if defined( __netcdf4_parallel )
              IF ( netcdf_data_format > 4 )  THEN
@@ -1235,5 +1277,5 @@
                                              id_var_do3d(av,i), &
                                              1, 0 )
-                CALL handle_netcdf_error( 'netcdf', 532 )
+                CALL netcdf_handle_error( 'netcdf_define_header', 532 )
 !
 !--             Set collective io operations for parallel io
@@ -1241,7 +1283,9 @@
                                                id_var_do3d(av,i), &
                                                NF90_COLLECTIVE )
-                CALL handle_netcdf_error( 'netcdf', 445 )
+                CALL netcdf_handle_error( 'netcdf_define_header', 445 )
              ENDIF
 #endif
+             var_list = TRIM( var_list ) // TRIM( do3d(av,i) ) // ';'
+
              i = i + 1
 
@@ -1257,5 +1301,5 @@
           nc_stat = NF90_PUT_ATT( id_set_3d(av), NF90_GLOBAL, 'VAR_LIST', &
                                   var_list )
-          CALL handle_netcdf_error( 'netcdf', 81 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 81 )
 
 !
@@ -1263,10 +1307,10 @@
 !--       parallel output.
           nc_stat = NF90_SET_FILL( id_set_3d(av), NF90_NOFILL, oldmode )
-          CALL handle_netcdf_error( 'netcdf', 528 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 528 )
 
 !
 !--       Leave netCDF define mode
           nc_stat = NF90_ENDDEF( id_set_3d(av) )
-          CALL handle_netcdf_error( 'netcdf', 82 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 82 )
 
 !
@@ -1285,5 +1329,5 @@
                                      netcdf_data, start = (/ 1 /),    &
                                      count = (/ nx+2 /) )
-             CALL handle_netcdf_error( 'netcdf', 83 )
+             CALL netcdf_handle_error( 'netcdf_define_header', 83 )
 
              DO  i = 0, nx+1
@@ -1294,5 +1338,5 @@
                                      netcdf_data, start = (/ 1 /),    &
                                      count = (/ nx+2 /) )
-             CALL handle_netcdf_error( 'netcdf', 385 )
+             CALL netcdf_handle_error( 'netcdf_define_header', 385 )
 
              DEALLOCATE( netcdf_data )
@@ -1309,5 +1353,5 @@
                                      netcdf_data, start = (/ 1 /),    &
                                      count = (/ ny+2 /) )
-             CALL handle_netcdf_error( 'netcdf', 84 )
+             CALL netcdf_handle_error( 'netcdf_define_header', 84 )
 
              DO  i = 0, ny+1
@@ -1318,5 +1362,5 @@
                                      netcdf_data, start = (/ 1 /),    &
                                      count = (/ ny+2 /))
-             CALL handle_netcdf_error( 'netcdf', 387 )
+             CALL netcdf_handle_error( 'netcdf_define_header', 387 )
 
              DEALLOCATE( netcdf_data )
@@ -1327,5 +1371,5 @@
                                      zu(nzb:nz_do3d), start = (/ 1 /), &
                                      count = (/ nz_do3d-nzb+1 /) )
-             CALL handle_netcdf_error( 'netcdf', 85 )
+             CALL netcdf_handle_error( 'netcdf_define_header', 85 )
 
 
@@ -1333,5 +1377,5 @@
                                      zw(nzb:nz_do3d), start = (/ 1 /), &
                                      count = (/ nz_do3d-nzb+1 /) )
-             CALL handle_netcdf_error( 'netcdf', 86 )
+             CALL netcdf_handle_error( 'netcdf_define_header', 86 )
 
 !
@@ -1343,5 +1387,5 @@
                                         start = (/ 1, 1 /), &
                                         count = (/ nx+2, ny+2 /) )
-                CALL handle_netcdf_error( 'netcdf', 419 )
+                CALL netcdf_handle_error( 'netcdf_define_header', 419 )
 
                 nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zwwi_3d(av), &
@@ -1349,5 +1393,5 @@
                                         start = (/ 1, 1 /), &
                                         count = (/ nx+2, ny+2 /) )
-                CALL handle_netcdf_error( 'netcdf', 420 )
+                CALL netcdf_handle_error( 'netcdf_define_header', 420 )
 
              ENDIF
@@ -1359,5 +1403,5 @@
                                         - zs(nzb_soil:nzt_soil), start = (/ 1 /), &
                                         count = (/ nzt_soil-nzb_soil+1 /) )
-                CALL handle_netcdf_error( 'netcdf', 86 )
+                CALL netcdf_handle_error( 'netcdf_define_header', 86 )
              ENDIF
 
@@ -1373,5 +1417,5 @@
           nc_stat = NF90_GET_ATT( id_set_3d(av), NF90_GLOBAL, 'VAR_LIST', &
                                   var_list_old )
-          CALL handle_netcdf_error( 'netcdf', 87 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 87 )
 
           var_list = ';'
@@ -1402,14 +1446,14 @@
 !--       Get and compare the number of vertical gridpoints
           nc_stat = NF90_INQ_VARID( id_set_3d(av), 'zu_3d', id_var_zu_3d(av) )
-          CALL handle_netcdf_error( 'netcdf', 88 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 88 )
 
           nc_stat = NF90_INQUIRE_VARIABLE( id_set_3d(av), id_var_zu_3d(av), &
                                            dimids = id_dim_zu_3d_old )
-          CALL handle_netcdf_error( 'netcdf', 89 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 89 )
           id_dim_zu_3d(av) = id_dim_zu_3d_old(1)
 
           nc_stat = NF90_INQUIRE_DIMENSION( id_set_3d(av), id_dim_zu_3d(av), &
                                             len = nz_old )
-          CALL handle_netcdf_error( 'netcdf', 90 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 90 )
 
           IF ( nz_do3d-nzb+1 /= nz_old )  THEN
@@ -1431,9 +1475,9 @@
 !--       on the file.
           nc_stat = NF90_INQ_VARID( id_set_3d(av), 'time', id_var_time_3d(av) )
-          CALL handle_netcdf_error( 'netcdf', 91 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 91 )
 
           nc_stat = NF90_INQUIRE_VARIABLE( id_set_3d(av), id_var_time_3d(av), &
                                            dimids = id_dim_time_old )
-          CALL handle_netcdf_error( 'netcdf', 92 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 92 )
 
           id_dim_time_3d(av) = id_dim_time_old(1)
@@ -1441,5 +1485,5 @@
           nc_stat = NF90_INQUIRE_DIMENSION( id_set_3d(av), id_dim_time_3d(av), &
                                             len = ntime_count )
-          CALL handle_netcdf_error( 'netcdf', 93 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 93 )
 
 !
@@ -1454,5 +1498,5 @@
                                   start = (/ do3d_time_count(av) /), &
                                   count = (/ 1 /) )
-          CALL handle_netcdf_error( 'netcdf', 94 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 94 )
 
           IF ( last_time_coordinate(1) >= simulated_time )  THEN
@@ -1524,5 +1568,5 @@
              nc_stat = NF90_INQ_VARID( id_set_3d(av), TRIM( do3d(av,i) ), &
                                        id_var_do3d(av,i) )
-             CALL handle_netcdf_error( 'netcdf', 95 )
+             CALL netcdf_handle_error( 'netcdf_define_header', 95 )
 #if defined( __netcdf4_parallel )
 !
@@ -1532,5 +1576,5 @@
                                                id_var_do3d(av,i), &
                                                NF90_COLLECTIVE )
-                CALL handle_netcdf_error( 'netcdf', 453 )
+                CALL netcdf_handle_error( 'netcdf_define_header', 453 )
              ENDIF
 #endif
@@ -1552,11 +1596,11 @@
           ENDIF
           nc_stat = NF90_REDEF( id_set_3d(av) )
-          CALL handle_netcdf_error( 'netcdf', 429 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 429 )
           nc_stat = NF90_PUT_ATT( id_set_3d(av), NF90_GLOBAL, 'title', &
                                   TRIM( run_description_header ) //    &
                                   TRIM( time_average_text ) )
-          CALL handle_netcdf_error( 'netcdf', 96 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 96 )
           nc_stat = NF90_ENDDEF( id_set_3d(av) )
-          CALL handle_netcdf_error( 'netcdf', 430 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 430 )
           message_string = 'netCDF file for volume data ' //             &
                            TRIM( var ) // ' from previous run found.' // &
@@ -1570,5 +1614,5 @@
           nc_stat = NF90_PUT_ATT( id_set_xy(av), NF90_GLOBAL, 'Conventions', &
                                   'COARDS' )
-          CALL handle_netcdf_error( 'netcdf', 97 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 97 )
 
           IF ( av == 0 )  THEN
@@ -1581,10 +1625,10 @@
                                   TRIM( run_description_header ) //    &
                                   TRIM( time_average_text ) )
-          CALL handle_netcdf_error( 'netcdf', 98 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 98 )
           IF ( av == 1 )  THEN
              WRITE ( time_average_text,'(F7.1,'' s avg'')' )  averaging_interval
              nc_stat = NF90_PUT_ATT( id_set_xy(av), NF90_GLOBAL, 'time_avg', &
                                      TRIM( time_average_text ) )
-             CALL handle_netcdf_error( 'netcdf', 98 )
+             CALL netcdf_handle_error( 'netcdf_define_header', 98 )
           ENDIF
 
@@ -1594,21 +1638,14 @@
 !--       the performance drops significantly.
           IF ( netcdf_data_format < 5 )  THEN
-             nc_stat = NF90_DEF_DIM( id_set_xy(av), 'time', NF90_UNLIMITED, &
-                                     id_dim_time_xy(av) )
-             CALL handle_netcdf_error( 'netcdf', 99 )
+             CALL netcdf_create_dim( id_set_xy(av), 'time', NF90_UNLIMITED,    &
+                                     id_dim_time_xy(av), 99 )
           ELSE
-             nc_stat = NF90_DEF_DIM( id_set_xy(av), 'time', ntdim_2d_xy(av), &
-                                     id_dim_time_xy(av) )
-             CALL handle_netcdf_error( 'netcdf', 524 )
-          ENDIF
-
-          nc_stat = NF90_DEF_VAR( id_set_xy(av), 'time', NF90_DOUBLE, &
-                                  id_dim_time_xy(av), id_var_time_xy(av) )
-          CALL handle_netcdf_error( 'netcdf', 100 )
-
-          nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_time_xy(av), 'units', &
-                                  'seconds')
-          CALL handle_netcdf_error( 'netcdf', 101 )
-
+             CALL netcdf_create_dim( id_set_xy(av), 'time', ntdim_2d_xy(av),   &
+                                     id_dim_time_xy(av), 524 )
+          ENDIF
+
+          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_time_xy(av) /),     &
+                                  'time', NF90_DOUBLE, id_var_time_xy(av),     &
+                                  'seconds', '', 100, 101, 000 )
 !
 !--       Define the spatial dimensions and coordinates for xy-sections.
@@ -1626,28 +1663,16 @@
 !
 !--       Define vertical coordinate grid (zu grid)
-          nc_stat = NF90_DEF_DIM( id_set_xy(av), 'zu_xy', ns, id_dim_zu_xy(av) )
-          CALL handle_netcdf_error( 'netcdf', 102 )
-
-          nc_stat = NF90_DEF_VAR( id_set_xy(av), 'zu_xy', NF90_DOUBLE, &
-                                  id_dim_zu_xy(av), id_var_zu_xy(av) )
-          CALL handle_netcdf_error( 'netcdf', 103 )
-
-          nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_zu_xy(av), 'units', &
-                                  'meters' )
-          CALL handle_netcdf_error( 'netcdf', 104 )
-
+          CALL netcdf_create_dim( id_set_xy(av), 'zu_xy', ns,                  &
+                                  id_dim_zu_xy(av), 102 )
+          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_zu_xy(av) /),       &
+                                  'zu_xy', NF90_DOUBLE, id_var_zu_xy(av),      &
+                                  'meters', '', 103, 104, 000 )
 !
 !--       Define vertical coordinate grid (zw grid)
-          nc_stat = NF90_DEF_DIM( id_set_xy(av), 'zw_xy', ns, id_dim_zw_xy(av) )
-          CALL handle_netcdf_error( 'netcdf', 105 )
-
-          nc_stat = NF90_DEF_VAR( id_set_xy(av), 'zw_xy', NF90_DOUBLE, &
-                                  id_dim_zw_xy(av), id_var_zw_xy(av) )
-          CALL handle_netcdf_error( 'netcdf', 106 )
-
-          nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_zw_xy(av), 'units', &
-                                  'meters' )
-          CALL handle_netcdf_error( 'netcdf', 107 )
-
+          CALL netcdf_create_dim( id_set_xy(av), 'zw_xy', ns,                  &
+                                  id_dim_zw_xy(av), 105 )
+          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_zw_xy(av) /),       &
+                                  'zw_xy', NF90_DOUBLE, id_var_zw_xy(av),      &
+                                  'meters', '', 106, 107, 000 )
 
           IF ( land_surface )  THEN
@@ -1659,14 +1684,9 @@
 !
 !--          Define vertical coordinate grid (zs grid)
-             nc_stat = NF90_DEF_DIM( id_set_xy(av), 'zs_xy', ns_do, id_dim_zs_xy(av) )
-             CALL handle_netcdf_error( 'netcdf', 539 )
-
-             nc_stat = NF90_DEF_VAR( id_set_xy(av), 'zs_xy', NF90_DOUBLE, &
-                                     id_dim_zs_xy(av), id_var_zs_xy(av) )
-             CALL handle_netcdf_error( 'netcdf', 540 )
-
-             nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_zs_xy(av), 'units', &
-                                  'meters' )
-             CALL handle_netcdf_error( 'netcdf', 541 )
+             CALL netcdf_create_dim( id_set_xy(av), 'zs_xy', ns_do,            &
+                                     id_dim_zs_xy(av), 539 )
+             CALL netcdf_create_var( id_set_xy(av), (/ id_dim_zs_xy(av) /),    &
+                                     'zs_xy', NF90_DOUBLE, id_var_zs_xy(av),   &
+                                     'meters', '', 540, 541, 000 )
 
           ENDIF
@@ -1675,80 +1695,44 @@
 !--       Define a pseudo vertical coordinate grid for the surface variables
 !--       u* and t* to store their height level
-          nc_stat = NF90_DEF_DIM( id_set_xy(av), 'zu1_xy', 1, &
-                                  id_dim_zu1_xy(av) )
-          CALL handle_netcdf_error( 'netcdf', 108 )
-
-          nc_stat = NF90_DEF_VAR( id_set_xy(av), 'zu1_xy', NF90_DOUBLE, &
-                                  id_dim_zu1_xy(av), id_var_zu1_xy(av) )
-          CALL handle_netcdf_error( 'netcdf', 109 )
-
-          nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_zu1_xy(av), 'units', &
-                                  'meters' )
-          CALL handle_netcdf_error( 'netcdf', 110 )
-
+          CALL netcdf_create_dim( id_set_xy(av), 'zu1_xy', 1,                  &
+                                  id_dim_zu1_xy(av), 108 )
+          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_zu1_xy(av) /),      &
+                                  'zu1_xy', NF90_DOUBLE, id_var_zu1_xy(av),    &
+                                  'meters', '', 109, 110, 000 )
 !
 !--       Define a variable to store the layer indices of the horizontal cross
 !--       sections, too
-          nc_stat = NF90_DEF_VAR( id_set_xy(av), 'ind_z_xy', NF90_DOUBLE, &
-                                  id_dim_zu_xy(av), id_var_ind_z_xy(av) )
-          CALL handle_netcdf_error( 'netcdf', 111 )
-
-          nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_ind_z_xy(av), 'units', &
-                                  'gridpoints')
-          CALL handle_netcdf_error( 'netcdf', 112 )
-
+          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_zu_xy(av) /),       &
+                                  'ind_z_xy', NF90_DOUBLE,                     &
+                                  id_var_ind_z_xy(av), 'gridpoints', '', 111,  &
+                                  112, 000 )
 !
 !--       Define x-axis (for scalar position)
-          nc_stat = NF90_DEF_DIM( id_set_xy(av), 'x', nx+2, id_dim_x_xy(av) )
-          CALL handle_netcdf_error( 'netcdf', 113 )
-
-          nc_stat = NF90_DEF_VAR( id_set_xy(av), 'x', NF90_DOUBLE, &
-                                  id_dim_x_xy(av), id_var_x_xy(av) )
-          CALL handle_netcdf_error( 'netcdf', 114 )
-
-
-          nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_x_xy(av), 'units', &
-                                  'meters' )
-          CALL handle_netcdf_error( 'netcdf', 115 )
-
+          CALL netcdf_create_dim( id_set_xy(av), 'x', nx+2, id_dim_x_xy(av),   &
+                                  113 )
+          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_x_xy(av) /), 'x',   &
+                                  NF90_DOUBLE, id_var_x_xy(av), 'meters', '',  &
+                                  114, 115, 000 )
 !
 !--       Define x-axis (for u position)
-          nc_stat = NF90_DEF_DIM( id_set_xy(av), 'xu', nx+2, id_dim_xu_xy(av) )
-          CALL handle_netcdf_error( 'netcdf', 388 )
-
-          nc_stat = NF90_DEF_VAR( id_set_xy(av), 'xu', NF90_DOUBLE, &
-                                  id_dim_xu_xy(av), id_var_xu_xy(av) )
-          CALL handle_netcdf_error( 'netcdf', 389 )
-
-          nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_xu_xy(av), 'units', &
-                                  'meters' )
-          CALL handle_netcdf_error( 'netcdf', 390 )
-
+          CALL netcdf_create_dim( id_set_xy(av), 'xu', nx+2,                   &
+                                  id_dim_xu_xy(av), 388 )
+          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_xu_xy(av) /), 'xu', &
+                                  NF90_DOUBLE, id_var_xu_xy(av), 'meters', '', &
+                                  389, 390, 000 )
 !
 !--       Define y-axis (for scalar position)
-          nc_stat = NF90_DEF_DIM( id_set_xy(av), 'y', ny+2, id_dim_y_xy(av) )
-          CALL handle_netcdf_error( 'netcdf', 116 )
-
-          nc_stat = NF90_DEF_VAR( id_set_xy(av), 'y', NF90_DOUBLE, &
-                                  id_dim_y_xy(av), id_var_y_xy(av) )
-          CALL handle_netcdf_error( 'netcdf', 117 )
-
-          nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_y_xy(av), 'units', &
-                                  'meters' )
-          CALL handle_netcdf_error( 'netcdf', 118 )
-
+          CALL netcdf_create_dim( id_set_xy(av), 'y', ny+2, id_dim_y_xy(av),   &
+                                  116 )
+          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_y_xy(av) /), 'y',   &
+                                  NF90_DOUBLE, id_var_y_xy(av), 'meters', '',  &
+                                  117, 118, 000 )
 !
 !--       Define y-axis (for scalar position)
-          nc_stat = NF90_DEF_DIM( id_set_xy(av), 'yv', ny+2, id_dim_yv_xy(av) )
-          CALL handle_netcdf_error( 'netcdf', 364 )
-
-          nc_stat = NF90_DEF_VAR( id_set_xy(av), 'yv', NF90_DOUBLE, &
-                                  id_dim_yv_xy(av), id_var_yv_xy(av) )
-          CALL handle_netcdf_error( 'netcdf', 365 )
-
-          nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_yv_xy(av), 'units', &
-                                  'meters' )
-          CALL handle_netcdf_error( 'netcdf', 366 )
-
+          CALL netcdf_create_dim( id_set_xy(av), 'yv', ny+2,                   &
+                                  id_dim_yv_xy(av), 364 )
+          CALL netcdf_create_var( id_set_xy(av), (/ id_dim_yv_xy(av) /), 'yv', &
+                                  NF90_DOUBLE, id_var_yv_xy(av), 'meters', '', &
+                                  365, 366, 000 )
 !
 !--       In case of non-flat topography define 2d-arrays containing the height
@@ -1757,31 +1741,14 @@
 !
 !--          Define zusi = zu(nzb_s_inner)
-             nc_stat = NF90_DEF_VAR( id_set_xy(av), 'zusi', NF90_DOUBLE, &
-                                     (/ id_dim_x_xy(av), id_dim_y_xy(av) /), &
-                                     id_var_zusi_xy(av) )
-             CALL handle_netcdf_error( 'netcdf', 421 )
-             
-             nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_zusi_xy(av), &
-                                     'units', 'meters' )
-             CALL handle_netcdf_error( 'netcdf', 422 )
-             
-             nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_zusi_xy(av), &
-                                     'long_name', 'zu(nzb_s_inner)' )
-             CALL handle_netcdf_error( 'netcdf', 423 )
-
+             CALL netcdf_create_var( id_set_xy(av), (/ id_dim_x_xy(av),        &
+                                     id_dim_y_xy(av) /), 'zusi', NF90_DOUBLE,  &
+                                     id_var_zusi_xy(av), 'meters',             &
+                                     'zu(nzb_s_inner)', 421, 422, 423 )
 !             
 !--          Define zwwi = zw(nzb_w_inner)
-             nc_stat = NF90_DEF_VAR( id_set_xy(av), 'zwwi', NF90_DOUBLE, &
-                                     (/ id_dim_x_xy(av), id_dim_y_xy(av) /), &
-                                     id_var_zwwi_xy(av) )
-             CALL handle_netcdf_error( 'netcdf', 424 )
-             
-             nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_zwwi_xy(av), &
-                                     'units', 'meters' )
-             CALL handle_netcdf_error( 'netcdf', 425 )
-             
-             nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_zwwi_xy(av), &
-                                     'long_name', 'zw(nzb_w_inner)' )
-             CALL handle_netcdf_error( 'netcdf', 426 )
+             CALL netcdf_create_var( id_set_xy(av), (/ id_dim_x_xy(av),        &
+                                     id_dim_y_xy(av) /), 'zwwi', NF90_DOUBLE,  &
+                                     id_var_zwwi_xy(av), 'meters',             &
+                                     'zw(nzb_w_inner)', 424, 425, 426 )
 
           ENDIF
@@ -1800,12 +1767,10 @@
                 IF ( INDEX( do2d(av,i), '*' ) /= 0 )  THEN
 
-                   nc_stat = NF90_DEF_VAR( id_set_xy(av),  do2d(av,i),         &
-                                           nc_precision(1),                    &
-                                           (/ id_dim_x_xy(av), id_dim_y_xy(av),&
-                                              id_dim_zu1_xy(av),               &
-                                              id_dim_time_xy(av) /),           &
-                                           id_var_do2d(av,i) )
-
-                   var_list = TRIM(var_list) // TRIM( do2d(av,i) ) // ';'
+                   CALL netcdf_create_var( id_set_xy(av), (/ id_dim_x_xy(av),  &
+                                           id_dim_y_xy(av), id_dim_zu1_xy(av), &
+                                           id_dim_time_xy(av) /), do2d(av,i),  &
+                                           nc_precision(1), id_var_do2d(av,i), &
+                                           TRIM( do2d_unit(av,i) ),            &
+                                           do2d(av,i), 119, 120, 354 )
 
                 ELSE
@@ -1862,4 +1827,11 @@
                                                        grid_x, grid_y, grid_z )
 
+                         IF ( .NOT. found )  THEN
+                            WRITE ( message_string, * ) 'no grid defined for', &
+                                                ' variable ', TRIM( do2d(av,i) )
+                            CALL message( 'define_netcdf_header', 'PA0244',    &
+                                          0, 1, 0, 6, 0 )
+                         ENDIF
+
                    END SELECT
 
@@ -1888,33 +1860,12 @@
 !
 !--                Define the grid
-                   nc_stat = NF90_DEF_VAR( id_set_xy(av), do2d(av,i),          &
-                                           nc_precision(1),                    &
-                                   (/ id_x, id_y, id_z, id_dim_time_xy(av) /), &
-                                           id_var_do2d(av,i) )
-
-                   IF ( .NOT. found )  THEN
-                      WRITE ( message_string, * ) 'no grid defined for', &
-                                            ' variable ', TRIM( do2d(av,i) )
-                      CALL message( 'define_netcdf_header', 'PA0244',&
-                                                                0, 1, 0, 6, 0 )
-                   ENDIF
-
-                   var_list = TRIM( var_list ) // TRIM( do2d(av,i) ) // ';'
+                   CALL netcdf_create_var( id_set_xy(av), (/ id_x, id_y, id_z, &
+                                           id_dim_time_xy(av) /), do2d(av,i),  &
+                                           nc_precision(1), id_var_do2d(av,i), &
+                                           TRIM( do2d_unit(av,i) ),            &
+                                           do2d(av,i), 119, 120, 354 )
 
                 ENDIF
 
-                CALL handle_netcdf_error( 'netcdf', 119 )
-!
-!--             Store the 'real' name of the variable (with *, for example)
-!--             in the long_name attribute. This is evaluated by Ferret,
-!--             for example.
-                nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_do2d(av,i), &
-                                        'long_name', do2d(av,i) )
-                CALL handle_netcdf_error( 'netcdf', 120 )
-!
-!--             Define the variable's unit
-                nc_stat = NF90_PUT_ATT( id_set_xy(av), id_var_do2d(av,i), &
-                                        'units', TRIM( do2d_unit(av,i) ) )
-                CALL handle_netcdf_error( 'netcdf', 354 ) 
 #if defined( __netcdf4_parallel )
                 IF ( netcdf_data_format > 4 )  THEN
@@ -1924,5 +1875,5 @@
                                                 id_var_do2d(av,i), &
                                                 1, 0 )
-                   CALL handle_netcdf_error( 'netcdf', 533 )
+                   CALL netcdf_handle_error( 'netcdf_define_header', 533 )
 !
 !--                Set collective io operations for parallel io
@@ -1930,7 +1881,9 @@
                                                   id_var_do2d(av,i), &
                                                   NF90_COLLECTIVE )
-                   CALL handle_netcdf_error( 'netcdf', 448 )
+                   CALL netcdf_handle_error( 'netcdf_define_header', 448 )
                 ENDIF
 #endif
+                var_list = TRIM( var_list) // TRIM( do2d(av,i) ) // ';'
+
              ENDIF
 
@@ -1948,5 +1901,5 @@
           nc_stat = NF90_PUT_ATT( id_set_xy(av), NF90_GLOBAL, 'VAR_LIST', &
                                   var_list )
-          CALL handle_netcdf_error( 'netcdf', 121 ) 
+          CALL netcdf_handle_error( 'netcdf_define_header', 121 )
 
 !
@@ -1954,10 +1907,10 @@
 !--       parallel output.
           nc_stat = NF90_SET_FILL( id_set_xy(av), NF90_NOFILL, oldmode )
-          CALL handle_netcdf_error( 'netcdf', 529 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 529 )
 
 !
 !--       Leave netCDF define mode
           nc_stat = NF90_ENDDEF( id_set_xy(av) )
-          CALL handle_netcdf_error( 'netcdf', 122 ) 
+          CALL netcdf_handle_error( 'netcdf_define_header', 122 )
 
 !
@@ -1982,5 +1935,5 @@
                                      netcdf_data, start = (/ 1 /),    &
                                      count = (/ ns /) )
-             CALL handle_netcdf_error( 'netcdf', 123 )
+             CALL netcdf_handle_error( 'netcdf_define_header', 123 )
 
 !
@@ -1996,5 +1949,5 @@
                                      netcdf_data, start = (/ 1 /),    &
                                      count = (/ ns /) )
-             CALL handle_netcdf_error( 'netcdf', 124 )
+             CALL netcdf_handle_error( 'netcdf_define_header', 124 )
 
 !
@@ -2015,5 +1968,5 @@
                                         netcdf_data(1:ns_do), start = (/ 1 /),    &
                                         count = (/ ns_do /) )
-                CALL handle_netcdf_error( 'netcdf', 124 )
+                CALL netcdf_handle_error( 'netcdf_define_header', 124 )
 
              ENDIF
@@ -2025,5 +1978,5 @@
                                      netcdf_data, start = (/ 1 /),       &
                                      count = (/ ns /) )
-             CALL handle_netcdf_error( 'netcdf', 125 )
+             CALL netcdf_handle_error( 'netcdf_define_header', 125 )
 
              DEALLOCATE( netcdf_data )
@@ -2034,5 +1987,5 @@
                                      (/ zu(nzb+1) /), start = (/ 1 /), &
                                      count = (/ 1 /) )
-             CALL handle_netcdf_error( 'netcdf', 126 )
+             CALL netcdf_handle_error( 'netcdf_define_header', 126 )
 
 !
@@ -2047,5 +2000,5 @@
                                      netcdf_data, start = (/ 1 /),   &
                                      count = (/ nx+2 /) )
-             CALL handle_netcdf_error( 'netcdf', 127 )
+             CALL netcdf_handle_error( 'netcdf_define_header', 127 )
 
              DO  i = 0, nx+1
@@ -2056,5 +2009,5 @@
                                      netcdf_data, start = (/ 1 /),    &
                                      count = (/ nx+2 /) )
-             CALL handle_netcdf_error( 'netcdf', 367 )
+             CALL netcdf_handle_error( 'netcdf_define_header', 367 )
 
              DEALLOCATE( netcdf_data )
@@ -2071,5 +2024,5 @@
                                      netcdf_data, start = (/ 1 /),   &
                                      count = (/ ny+2 /))
-             CALL handle_netcdf_error( 'netcdf', 128 )
+             CALL netcdf_handle_error( 'netcdf_define_header', 128 )
 
              DO  i = 0, ny+1
@@ -2080,5 +2033,5 @@
                                      netcdf_data, start = (/ 1 /),    &
                                      count = (/ ny+2 /))
-             CALL handle_netcdf_error( 'netcdf', 368 )
+             CALL netcdf_handle_error( 'netcdf_define_header', 368 )
 
              DEALLOCATE( netcdf_data )
@@ -2092,5 +2045,5 @@
                                         start = (/ 1, 1 /), &
                                         count = (/ nx+2, ny+2 /) )
-                CALL handle_netcdf_error( 'netcdf', 427 )
+                CALL netcdf_handle_error( 'netcdf_define_header', 427 )
 
                 nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_zwwi_xy(av), &
@@ -2098,5 +2051,5 @@
                                         start = (/ 1, 1 /), &
                                         count = (/ nx+2, ny+2 /) )
-                CALL handle_netcdf_error( 'netcdf', 428 )
+                CALL netcdf_handle_error( 'netcdf_define_header', 428 )
 
              ENDIF
@@ -2115,5 +2068,5 @@
           nc_stat = NF90_GET_ATT( id_set_xy(av), NF90_GLOBAL, 'VAR_LIST', &
                                   var_list_old )
-          CALL handle_netcdf_error( 'netcdf', 129 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 129 )
 
           var_list = ';'
@@ -2154,14 +2107,14 @@
 !--       Get and compare the number of horizontal cross sections
           nc_stat = NF90_INQ_VARID( id_set_xy(av), 'zu_xy', id_var_zu_xy(av) )
-          CALL handle_netcdf_error( 'netcdf', 130 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 130 )
 
           nc_stat = NF90_INQUIRE_VARIABLE( id_set_xy(av), id_var_zu_xy(av), &
                                            dimids = id_dim_zu_xy_old )
-          CALL handle_netcdf_error( 'netcdf', 131 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 131 )
           id_dim_zu_xy(av) = id_dim_zu_xy_old(1)
 
           nc_stat = NF90_INQUIRE_DIMENSION( id_set_xy(av), id_dim_zu_xy(av), &
                                             len = ns_old )
-          CALL handle_netcdf_error( 'netcdf', 132 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 132 )
 
           IF ( ns /= ns_old )  THEN
@@ -2182,5 +2135,5 @@
 
           nc_stat = NF90_GET_VAR( id_set_xy(av), id_var_zu_xy(av), netcdf_data )
-          CALL handle_netcdf_error( 'netcdf', 133 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 133 )
 
           DO  i = 1, ns
@@ -2222,14 +2175,14 @@
 !--       on the file.
           nc_stat = NF90_INQ_VARID( id_set_xy(av), 'time', id_var_time_xy(av) )
-          CALL handle_netcdf_error( 'netcdf', 134 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 134 )
 
           nc_stat = NF90_INQUIRE_VARIABLE( id_set_xy(av), id_var_time_xy(av), &
                                            dimids = id_dim_time_old )
-          CALL handle_netcdf_error( 'netcdf', 135 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 135 )
           id_dim_time_xy(av) = id_dim_time_old(1)
 
           nc_stat = NF90_INQUIRE_DIMENSION( id_set_xy(av), id_dim_time_xy(av), &
                                             len = ntime_count )
-          CALL handle_netcdf_error( 'netcdf', 136 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 136 )
 
 !
@@ -2244,5 +2197,5 @@
                                   start = (/ do2d_xy_time_count(av) /), &
                                   count = (/ 1 /) )
-          CALL handle_netcdf_error( 'netcdf', 137 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 137 )
 
           IF ( last_time_coordinate(1) >= simulated_time )  THEN
@@ -2315,5 +2268,5 @@
                 nc_stat = NF90_INQ_VARID( id_set_xy(av), do2d(av,i), &
                                           id_var_do2d(av,i) )
-                CALL handle_netcdf_error( 'netcdf', 138 )
+                CALL netcdf_handle_error( 'netcdf_define_header', 138 )
 #if defined( __netcdf4_parallel )
 !
@@ -2323,5 +2276,5 @@
                                                   id_var_do2d(av,i), &
                                                   NF90_COLLECTIVE )
-                   CALL handle_netcdf_error( 'netcdf', 454 )
+                   CALL netcdf_handle_error( 'netcdf_define_header', 454 )
                 ENDIF
 #endif
@@ -2344,11 +2297,11 @@
           ENDIF
           nc_stat = NF90_REDEF( id_set_xy(av) )
-          CALL handle_netcdf_error( 'netcdf', 431 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 431 )
           nc_stat = NF90_PUT_ATT( id_set_xy(av), NF90_GLOBAL, 'title', &
                                   TRIM( run_description_header ) //    &
                                   TRIM( time_average_text ) )
-          CALL handle_netcdf_error( 'netcdf', 139 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 139 )
           nc_stat = NF90_ENDDEF( id_set_xy(av) )
-          CALL handle_netcdf_error( 'netcdf', 432 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 432 )
           message_string = 'netCDF file for cross-sections ' //           &
                             TRIM( var ) // ' from previous run found.' // &
@@ -2363,5 +2316,5 @@
           nc_stat = NF90_PUT_ATT( id_set_xz(av), NF90_GLOBAL, 'Conventions', &
                                   'COARDS' )
-          CALL handle_netcdf_error( 'netcdf', 140 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 140 )
 
           IF ( av == 0 )  THEN
@@ -2374,10 +2327,10 @@
                                   TRIM( run_description_header )  //   &
                                   TRIM( time_average_text ) )
-          CALL handle_netcdf_error( 'netcdf', 141 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 141 )
           IF ( av == 1 )  THEN
              WRITE ( time_average_text,'(F7.1,'' s avg'')' )  averaging_interval
              nc_stat = NF90_PUT_ATT( id_set_xz(av), NF90_GLOBAL, 'time_avg', &
                                      TRIM( time_average_text ) )
-             CALL handle_netcdf_error( 'netcdf', 141 )
+             CALL netcdf_handle_error( 'netcdf_define_header', 141 )
           ENDIF
 
@@ -2387,21 +2340,14 @@
 !--       the performance drops significantly.
           IF ( netcdf_data_format < 5 )  THEN
-             nc_stat = NF90_DEF_DIM( id_set_xz(av), 'time', NF90_UNLIMITED, &
-                                     id_dim_time_xz(av) )
-             CALL handle_netcdf_error( 'netcdf', 142 )
+             CALL netcdf_create_dim( id_set_xz(av), 'time', NF90_UNLIMITED,    &
+                                     id_dim_time_xz(av), 142 )
           ELSE
-             nc_stat = NF90_DEF_DIM( id_set_xz(av), 'time', ntdim_2d_xz(av), &
-                                     id_dim_time_xz(av) )
-             CALL handle_netcdf_error( 'netcdf', 525 )
-          ENDIF
-
-          nc_stat = NF90_DEF_VAR( id_set_xz(av), 'time', NF90_DOUBLE, &
-                                  id_dim_time_xz(av), id_var_time_xz(av) )
-          CALL handle_netcdf_error( 'netcdf', 143 )
-
-          nc_stat = NF90_PUT_ATT( id_set_xz(av), id_var_time_xz(av), 'units', &
-                                  'seconds')
-          CALL handle_netcdf_error( 'netcdf', 144 )
-
+             CALL netcdf_create_dim( id_set_xz(av), 'time', ntdim_2d_xz(av),   &
+                                     id_dim_time_xz(av), 525 )
+          ENDIF
+
+          CALL netcdf_create_var( id_set_xz(av), (/ id_dim_time_xz(av) /),     &
+                                  'time', NF90_DOUBLE, id_var_time_xz(av),     &
+                                  'seconds', '', 143, 144, 000 )
 !
 !--       Define the spatial dimensions and coordinates for xz-sections.
@@ -2419,101 +2365,57 @@
 !
 !--       Define y-axis (for scalar position)
-          nc_stat = NF90_DEF_DIM( id_set_xz(av), 'y_xz', ns, id_dim_y_xz(av) )
-          CALL handle_netcdf_error( 'netcdf', 145 )
-
-          nc_stat = NF90_DEF_VAR( id_set_xz(av), 'y_xz', NF90_DOUBLE, &
-                                  id_dim_y_xz(av), id_var_y_xz(av) )
-          CALL handle_netcdf_error( 'netcdf', 146 )
-
-          nc_stat = NF90_PUT_ATT( id_set_xz(av), id_var_y_xz(av), 'units', &
-                                  'meters' )
-          CALL handle_netcdf_error( 'netcdf', 147 )
-
+          CALL netcdf_create_dim( id_set_xz(av), 'y_xz', ns, id_dim_y_xz(av),  &
+                                  145 )
+          CALL netcdf_create_var( id_set_xz(av), (/ id_dim_y_xz(av) /),        &
+                                  'y_xz', NF90_DOUBLE, id_var_y_xz(av),        &
+                                  'meters', '', 146, 147, 000 )
 !
 !--       Define y-axis (for v position)
-          nc_stat = NF90_DEF_DIM( id_set_xz(av), 'yv_xz', ns, id_dim_yv_xz(av) )
-          CALL handle_netcdf_error( 'netcdf', 369 )
-
-          nc_stat = NF90_DEF_VAR( id_set_xz(av), 'yv_xz', NF90_DOUBLE, &
-                                  id_dim_yv_xz(av), id_var_yv_xz(av) )
-          CALL handle_netcdf_error( 'netcdf', 370 )
-
-          nc_stat = NF90_PUT_ATT( id_set_xz(av), id_var_yv_xz(av), 'units', &
-                                  'meters' )
-          CALL handle_netcdf_error( 'netcdf', 371 )
-
+          CALL netcdf_create_dim( id_set_xz(av), 'yv_xz', ns,                  &
+                                  id_dim_yv_xz(av), 369 )
+          CALL netcdf_create_var( id_set_xz(av), (/ id_dim_yv_xz(av) /),       &
+                                  'yv_xz', NF90_DOUBLE, id_var_yv_xz(av),      &
+                                  'meters', '', 370, 371, 000 )
 !
 !--       Define a variable to store the layer indices of the vertical cross
 !--       sections
-          nc_stat = NF90_DEF_VAR( id_set_xz(av), 'ind_y_xz', NF90_DOUBLE, &
-                                  id_dim_y_xz(av), id_var_ind_y_xz(av) )
-          CALL handle_netcdf_error( 'netcdf', 148 )
-
-          nc_stat = NF90_PUT_ATT( id_set_xz(av), id_var_ind_y_xz(av), 'units', &
-                                  'gridpoints')
-          CALL handle_netcdf_error( 'netcdf', 149 )
-
+          CALL netcdf_create_var( id_set_xz(av), (/ id_dim_y_xz(av) /),        &
+                                  'ind_y_xz', NF90_DOUBLE,                     &
+                                  id_var_ind_y_xz(av), 'gridpoints', '', 148,  &
+                                  149, 000 )
 !
 !--       Define x-axis (for scalar position)
-          nc_stat = NF90_DEF_DIM( id_set_xz(av), 'x', nx+2, id_dim_x_xz(av) )
-          CALL handle_netcdf_error( 'netcdf', 150 )
-
-          nc_stat = NF90_DEF_VAR( id_set_xz(av), 'x', NF90_DOUBLE, &
-                                  id_dim_x_xz(av), id_var_x_xz(av) )
-          CALL handle_netcdf_error( 'netcdf', 151 )
-
-          nc_stat = NF90_PUT_ATT( id_set_xz(av), id_var_x_xz(av), 'units', &
-                                  'meters' )
-          CALL handle_netcdf_error( 'netcdf', 152 )
-
+          CALL netcdf_create_dim( id_set_xz(av), 'x', nx+2, id_dim_x_xz(av),   &
+                                  150 )
+          CALL netcdf_create_var( id_set_xz(av), (/ id_dim_x_xz(av) /), 'x',   &
+                                  NF90_DOUBLE, id_var_x_xz(av), 'meters', '',  &
+                                  151, 152, 000 )
 !
 !--       Define x-axis (for u position)
-          nc_stat = NF90_DEF_DIM( id_set_xz(av), 'xu', nx+2, id_dim_xu_xz(av) )
-          CALL handle_netcdf_error( 'netcdf', 372 )
-
-          nc_stat = NF90_DEF_VAR( id_set_xz(av), 'xu', NF90_DOUBLE, &
-                                  id_dim_xu_xz(av), id_var_xu_xz(av) )
-          CALL handle_netcdf_error( 'netcdf', 373 )
-
-          nc_stat = NF90_PUT_ATT( id_set_xz(av), id_var_xu_xz(av), 'units', &
-                                  'meters' )
-          CALL handle_netcdf_error( 'netcdf', 374 )
-
+          CALL netcdf_create_dim( id_set_xz(av), 'xu', nx+2, id_dim_xu_xz(av), &
+                                  372 )
+          CALL netcdf_create_var( id_set_xz(av), (/ id_dim_xu_xz(av) /), 'xu', &
+                                  NF90_DOUBLE, id_var_xu_xz(av), 'meters', '', &
+                                  373, 374, 000 )
 !
 !--       Define the three z-axes (zu, zw, and zs)
-          nc_stat = NF90_DEF_DIM( id_set_xz(av), 'zu', nz+2, id_dim_zu_xz(av) )
-          CALL handle_netcdf_error( 'netcdf', 153 )
-
-          nc_stat = NF90_DEF_VAR( id_set_xz(av), 'zu', NF90_DOUBLE, &
-                                  id_dim_zu_xz(av), id_var_zu_xz(av) )
-          CALL handle_netcdf_error( 'netcdf', 154 )
-
-          nc_stat = NF90_PUT_ATT( id_set_xz(av), id_var_zu_xz(av), 'units', &
-                                  'meters' )
-          CALL handle_netcdf_error( 'netcdf', 155 )
-
-          nc_stat = NF90_DEF_DIM( id_set_xz(av), 'zw', nz+2, id_dim_zw_xz(av) )
-          CALL handle_netcdf_error( 'netcdf', 156 )
-
-          nc_stat = NF90_DEF_VAR( id_set_xz(av), 'zw', NF90_DOUBLE, &
-                                  id_dim_zw_xz(av), id_var_zw_xz(av) )
-          CALL handle_netcdf_error( 'netcdf', 157 )
-
-          nc_stat = NF90_PUT_ATT( id_set_xz(av), id_var_zw_xz(av), 'units', &
-                                  'meters' )
-          CALL handle_netcdf_error( 'netcdf', 158 )
+          CALL netcdf_create_dim( id_set_xz(av), 'zu', nz+2, id_dim_zu_xz(av), &
+                                  153 )
+          CALL netcdf_create_var( id_set_xz(av), (/ id_dim_zu_xz(av) /), 'zu', &
+                                  NF90_DOUBLE, id_var_zu_xz(av), 'meters', '', &
+                                  154, 155, 000 )
+          CALL netcdf_create_dim( id_set_xz(av), 'zw', nz+2, id_dim_zw_xz(av), &
+                                  156 )
+          CALL netcdf_create_var( id_set_xz(av), (/ id_dim_zw_xz(av) /), 'zw', &
+                                  NF90_DOUBLE, id_var_zw_xz(av), 'meters', '', &
+                                  157, 158, 000 )
 
           IF ( land_surface )  THEN
 
-             nc_stat = NF90_DEF_DIM( id_set_xz(av), 'zs', nzs, id_dim_zs_xz(av) )
-             CALL handle_netcdf_error( 'netcdf', 542 )
-
-             nc_stat = NF90_DEF_VAR( id_set_xz(av), 'zs', NF90_DOUBLE, &
-                                     id_dim_zs_xz(av), id_var_zs_xz(av) )
-             CALL handle_netcdf_error( 'netcdf', 543 )
-
-             nc_stat = NF90_PUT_ATT( id_set_xz(av), id_var_zs_xz(av), 'units', &
-                                     'meters' )
-             CALL handle_netcdf_error( 'netcdf', 544 )
+             CALL netcdf_create_dim( id_set_xz(av), 'zs', nzs,                 &
+                                     id_dim_zs_xz(av), 542 )
+             CALL netcdf_create_var( id_set_xz(av), (/ id_dim_zs_xz(av) /),    &
+                                     'zs', NF90_DOUBLE, id_var_zs_xz(av),      &
+                                     'meters', '', 543, 544, 000 )
 
           ENDIF
@@ -2580,4 +2482,10 @@
                       CALL user_define_netcdf_grid( do2d(av,i), found, &
                                                     grid_x, grid_y, grid_z )
+                      IF ( .NOT. found )  THEN
+                         WRITE ( message_string, * ) 'no grid defined for', &
+                                                ' variable ', TRIM( do2d(av,i) )
+                         CALL message( 'define_netcdf_header', 'PA0244', &
+                                       0, 1, 0, 6, 0 )
+                      ENDIF
 
                 END SELECT
@@ -2607,31 +2515,10 @@
 !
 !--             Define the grid
-                nc_stat = NF90_DEF_VAR( id_set_xz(av), do2d(av,i),             &
-                                        nc_precision(2),                       &
-                                   (/ id_x, id_y, id_z, id_dim_time_xz(av) /), &
-                                        id_var_do2d(av,i) )
-
-                IF ( .NOT. found )  THEN
-                   WRITE ( message_string, * ) 'no grid defined for', &
-                                            ' variable ', TRIM( do2d(av,i) )
-                   CALL message( 'define_netcdf_header', 'PA0244',&
-                                                                 0, 1, 0, 6, 0 )
-                ENDIF
-
-                var_list = TRIM( var_list ) // TRIM( do2d(av,i) ) // ';'
-
-                CALL handle_netcdf_error( 'netcdf', 159 )
-!
-!--             Store the 'real' name of the variable (with *, for example)
-!--             in the long_name attribute. This is evaluated by Ferret,
-!--             for example.
-                nc_stat = NF90_PUT_ATT( id_set_xz(av), id_var_do2d(av,i), &
-                                        'long_name', do2d(av,i) )
-                CALL handle_netcdf_error( 'netcdf', 160 )
-!
-!--             Define the variable's unit
-                nc_stat = NF90_PUT_ATT( id_set_xz(av), id_var_do2d(av,i), &
-                                        'units', TRIM( do2d_unit(av,i) ) )
-                CALL handle_netcdf_error( 'netcdf', 355 )
+                CALL netcdf_create_var( id_set_xz(av), (/ id_x, id_y, id_z,    &
+                                        id_dim_time_xz(av) /), do2d(av,i),     &
+                                        nc_precision(2), id_var_do2d(av,i),    &
+                                        TRIM( do2d_unit(av,i) ), do2d(av,i),   &
+                                        159, 160, 355 )
+
 #if defined( __netcdf4_parallel )
 
@@ -2642,5 +2529,5 @@
                                                 id_var_do2d(av,i), &
                                                 1, 0 )
-                   CALL handle_netcdf_error( 'netcdf', 534 )
+                   CALL netcdf_handle_error( 'netcdf_define_header', 534 )
 !
 !--                Set independent io operations for parallel io. Collective io
@@ -2664,7 +2551,9 @@
 !                                                     NF90_INDEPENDENT )
                    ENDIF
-                   CALL handle_netcdf_error( 'netcdf', 449 )
+                   CALL netcdf_handle_error( 'netcdf_define_header', 449 )
                 ENDIF
 #endif
+                var_list = TRIM( var_list ) // TRIM( do2d(av,i) ) // ';'
+
              ENDIF
 
@@ -2682,5 +2571,5 @@
           nc_stat = NF90_PUT_ATT( id_set_xz(av), NF90_GLOBAL, 'VAR_LIST', &
                                   var_list )
-          CALL handle_netcdf_error( 'netcdf', 161 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 161 )
 
 !
@@ -2688,10 +2577,10 @@
 !--       parallel output.
           nc_stat = NF90_SET_FILL( id_set_xz(av), NF90_NOFILL, oldmode )
-          CALL handle_netcdf_error( 'netcdf', 530 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 530 )
 
 !
 !--       Leave netCDF define mode
           nc_stat = NF90_ENDDEF( id_set_xz(av) )
-          CALL handle_netcdf_error( 'netcdf', 162 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 162 )
 
 !
@@ -2716,5 +2605,5 @@
                                      netcdf_data, start = (/ 1 /),   &
                                      count = (/ ns /) )
-             CALL handle_netcdf_error( 'netcdf', 163 )
+             CALL netcdf_handle_error( 'netcdf_define_header', 163 )
 
 !
@@ -2730,5 +2619,5 @@
                                      netcdf_data, start = (/ 1 /),    &
                                      count = (/ ns /) )
-             CALL handle_netcdf_error( 'netcdf', 375 )
+             CALL netcdf_handle_error( 'netcdf_define_header', 375 )
 
 !
@@ -2738,5 +2627,5 @@
                                      netcdf_data, start = (/ 1 /),       &
                                      count = (/ ns /) )
-             CALL handle_netcdf_error( 'netcdf', 164 )
+             CALL netcdf_handle_error( 'netcdf_define_header', 164 )
 
 
@@ -2754,5 +2643,5 @@
                                      netcdf_data, start = (/ 1 /),   &
                                      count = (/ nx+2 /) )
-             CALL handle_netcdf_error( 'netcdf', 165 )
+             CALL netcdf_handle_error( 'netcdf_define_header', 165 )
 
              DO  i = 0, nx+1
@@ -2763,5 +2652,5 @@
                                      netcdf_data, start = (/ 1 /),    &
                                      count = (/ nx+2 /) )
-             CALL handle_netcdf_error( 'netcdf', 377 )
+             CALL netcdf_handle_error( 'netcdf_define_header', 377 )
 
              DEALLOCATE( netcdf_data )
@@ -2775,5 +2664,5 @@
                                      netcdf_data, start = (/ 1 /),    &
                                      count = (/ nz+2 /) )
-             CALL handle_netcdf_error( 'netcdf', 166 )
+             CALL netcdf_handle_error( 'netcdf_define_header', 166 )
 
              netcdf_data(0:nz+1) = zw(nzb:nzt+1)
@@ -2781,5 +2670,5 @@
                                      netcdf_data, start = (/ 1 /),    &
                                      count = (/ nz+2 /) )
-             CALL handle_netcdf_error( 'netcdf', 167 )
+             CALL netcdf_handle_error( 'netcdf_define_header', 167 )
 
 !
@@ -2790,5 +2679,5 @@
                                         netcdf_data(0:nzs), start = (/ 1 /),    &
                                         count = (/ nzt_soil-nzb_soil+1 /) )
-               CALL handle_netcdf_error( 'netcdf', 548 )
+               CALL netcdf_handle_error( 'netcdf_define_header', 548 )
              ENDIF
 
@@ -2808,5 +2697,5 @@
           nc_stat = NF90_GET_ATT( id_set_xz(av), NF90_GLOBAL, 'VAR_LIST', &
                                   var_list_old )
-          CALL handle_netcdf_error( 'netcdf', 168 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 168 )
 
           var_list = ';'
@@ -2847,14 +2736,14 @@
 !--       Get and compare the number of vertical cross sections
           nc_stat = NF90_INQ_VARID( id_set_xz(av), 'y_xz', id_var_y_xz(av) )
-          CALL handle_netcdf_error( 'netcdf', 169 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 169 )
 
           nc_stat = NF90_INQUIRE_VARIABLE( id_set_xz(av), id_var_y_xz(av), &
                                            dimids = id_dim_y_xz_old )
-          CALL handle_netcdf_error( 'netcdf', 170 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 170 )
           id_dim_y_xz(av) = id_dim_y_xz_old(1)
 
           nc_stat = NF90_INQUIRE_DIMENSION( id_set_xz(av), id_dim_y_xz(av), &
                                             len = ns_old )
-          CALL handle_netcdf_error( 'netcdf', 171 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 171 )
 
           IF ( ns /= ns_old )  THEN
@@ -2875,5 +2764,5 @@
 
           nc_stat = NF90_GET_VAR( id_set_xz(av), id_var_y_xz(av), netcdf_data )
-          CALL handle_netcdf_error( 'netcdf', 172 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 172 )
 
           DO  i = 1, ns
@@ -2915,14 +2804,14 @@
 !--       on the file.
           nc_stat = NF90_INQ_VARID( id_set_xz(av), 'time', id_var_time_xz(av) )
-          CALL handle_netcdf_error( 'netcdf', 173 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 173 )
 
           nc_stat = NF90_INQUIRE_VARIABLE( id_set_xz(av), id_var_time_xz(av), &
                                            dimids = id_dim_time_old )
-          CALL handle_netcdf_error( 'netcdf', 174 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 174 )
           id_dim_time_xz(av) = id_dim_time_old(1)
 
           nc_stat = NF90_INQUIRE_DIMENSION( id_set_xz(av), id_dim_time_xz(av), &
                                             len = ntime_count )
-          CALL handle_netcdf_error( 'netcdf', 175 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 175 )
 
 !
@@ -2937,5 +2826,5 @@
                                   start = (/ do2d_xz_time_count(av) /), &
                                   count = (/ 1 /) )
-          CALL handle_netcdf_error( 'netcdf', 176 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 176 )
 
           IF ( last_time_coordinate(1) >= simulated_time )  THEN
@@ -3008,5 +2897,5 @@
                 nc_stat = NF90_INQ_VARID( id_set_xz(av), do2d(av,i), &
                                           id_var_do2d(av,i) )
-                CALL handle_netcdf_error( 'netcdf', 177 )
+                CALL netcdf_handle_error( 'netcdf_define_header', 177 )
 #if defined( __netcdf4_parallel )
 !
@@ -3032,5 +2921,5 @@
 !                                                     NF90_INDEPENDENT )
                    ENDIF
-                   CALL handle_netcdf_error( 'netcdf', 455 )
+                   CALL netcdf_handle_error( 'netcdf_define_header', 455 )
                 ENDIF
 #endif
@@ -3053,11 +2942,11 @@
           ENDIF
           nc_stat = NF90_REDEF( id_set_xz(av) )
-          CALL handle_netcdf_error( 'netcdf', 433 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 433 )
           nc_stat = NF90_PUT_ATT( id_set_xz(av), NF90_GLOBAL, 'title', &
                                   TRIM( run_description_header ) //    &
                                   TRIM( time_average_text ) )
-          CALL handle_netcdf_error( 'netcdf', 178 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 178 )
           nc_stat = NF90_ENDDEF( id_set_xz(av) )
-          CALL handle_netcdf_error( 'netcdf', 434 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 434 )
           message_string = 'netCDF file for cross-sections ' //           &
                             TRIM( var ) // ' from previous run found.' // &
@@ -3072,5 +2961,5 @@
           nc_stat = NF90_PUT_ATT( id_set_yz(av), NF90_GLOBAL, 'Conventions', &
                                   'COARDS' )
-          CALL handle_netcdf_error( 'netcdf', 179 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 179 )
 
           IF ( av == 0 )  THEN
@@ -3083,10 +2972,10 @@
                                   TRIM( run_description_header ) //    &
                                   TRIM( time_average_text ) )
-          CALL handle_netcdf_error( 'netcdf', 180 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 180 )
           IF ( av == 1 )  THEN
              WRITE ( time_average_text,'(F7.1,'' s avg'')' )  averaging_interval
              nc_stat = NF90_PUT_ATT( id_set_yz(av), NF90_GLOBAL, 'time_avg', &
                                      TRIM( time_average_text ) )
-             CALL handle_netcdf_error( 'netcdf', 180 )
+             CALL netcdf_handle_error( 'netcdf_define_header', 180 )
           ENDIF
 
@@ -3096,21 +2985,14 @@
 !--       the performance drops significantly.
           IF ( netcdf_data_format < 5 )  THEN
-             nc_stat = NF90_DEF_DIM( id_set_yz(av), 'time', NF90_UNLIMITED, &
-                                     id_dim_time_yz(av) )
-             CALL handle_netcdf_error( 'netcdf', 181 )
+             CALL netcdf_create_dim( id_set_yz(av), 'time', NF90_UNLIMITED,    &
+                                     id_dim_time_yz(av), 181 )
           ELSE
-             nc_stat = NF90_DEF_DIM( id_set_yz(av), 'time', ntdim_2d_yz(av), &
-                                     id_dim_time_yz(av) )
-             CALL handle_netcdf_error( 'netcdf', 526 )
-          ENDIF
-
-          nc_stat = NF90_DEF_VAR( id_set_yz(av), 'time', NF90_DOUBLE, &
-                                  id_dim_time_yz(av), id_var_time_yz(av) )
-          CALL handle_netcdf_error( 'netcdf', 182 )
-
-          nc_stat = NF90_PUT_ATT( id_set_yz(av), id_var_time_yz(av), 'units', &
-                                  'seconds')
-          CALL handle_netcdf_error( 'netcdf', 183 )
-
+             CALL netcdf_create_dim( id_set_yz(av), 'time', ntdim_2d_yz(av),   &
+                                     id_dim_time_yz(av), 526 )
+          ENDIF
+
+          CALL netcdf_create_var( id_set_yz(av), (/ id_dim_time_yz(av) /),     &
+                                  'time', NF90_DOUBLE, id_var_time_yz(av),     &
+                                  'seconds', '', 182, 183, 000 )
 !
 !--       Define the spatial dimensions and coordinates for yz-sections.
@@ -3128,104 +3010,60 @@
 !
 !--       Define x axis (for scalar position)
-          nc_stat = NF90_DEF_DIM( id_set_yz(av), 'x_yz', ns, id_dim_x_yz(av) )
-          CALL handle_netcdf_error( 'netcdf', 184 )
-
-          nc_stat = NF90_DEF_VAR( id_set_yz(av), 'x_yz', NF90_DOUBLE, &
-                                  id_dim_x_yz(av), id_var_x_yz(av) )
-          CALL handle_netcdf_error( 'netcdf', 185 )
-
-          nc_stat = NF90_PUT_ATT( id_set_yz(av), id_var_x_yz(av), 'units', &
-                                  'meters' )
-          CALL handle_netcdf_error( 'netcdf', 186 )
-
+          CALL netcdf_create_dim( id_set_yz(av), 'x_yz', ns, id_dim_x_yz(av),  &
+                                  184 )
+          CALL netcdf_create_var( id_set_yz(av), (/ id_dim_x_yz(av) /),        &
+                                  'x_yz', NF90_DOUBLE, id_var_x_yz(av),        &
+                                  'meters', '', 185, 186, 000 )
 !
 !--       Define x axis (for u position)
-          nc_stat = NF90_DEF_DIM( id_set_yz(av), 'xu_yz', ns, id_dim_xu_yz(av) )
-          CALL handle_netcdf_error( 'netcdf', 377 )
-
-          nc_stat = NF90_DEF_VAR( id_set_yz(av), 'xu_yz', NF90_DOUBLE, &
-                                  id_dim_xu_yz(av), id_var_xu_yz(av) )
-          CALL handle_netcdf_error( 'netcdf', 378 )
-
-          nc_stat = NF90_PUT_ATT( id_set_yz(av), id_var_xu_yz(av), 'units', &
-                                  'meters' )
-          CALL handle_netcdf_error( 'netcdf', 379 )
-
+          CALL netcdf_create_dim( id_set_yz(av), 'xu_yz', ns,                  &
+                                  id_dim_xu_yz(av), 377 )
+          CALL netcdf_create_var( id_set_yz(av), (/ id_dim_xu_yz(av) /),       &
+                                  'xu_yz', NF90_DOUBLE, id_var_xu_yz(av),      &
+                                  'meters', '', 378, 379, 000 )
 !
 !--       Define a variable to store the layer indices of the vertical cross
 !--       sections
-          nc_stat = NF90_DEF_VAR( id_set_yz(av), 'ind_x_yz', NF90_DOUBLE, &
-                                  id_dim_x_yz(av), id_var_ind_x_yz(av) )
-          CALL handle_netcdf_error( 'netcdf', 187 )
-
-          nc_stat = NF90_PUT_ATT( id_set_yz(av), id_var_ind_x_yz(av), 'units', &
-                                  'gridpoints')
-          CALL handle_netcdf_error( 'netcdf', 188 )
-
+          CALL netcdf_create_var( id_set_yz(av), (/ id_dim_x_yz(av) /),        &
+                                  'ind_x_yz', NF90_DOUBLE,                     &
+                                  id_var_ind_x_yz(av), 'gridpoints', '', 187,  &
+                                  188, 000 )
 !
 !--       Define y-axis (for scalar position)
-          nc_stat = NF90_DEF_DIM( id_set_yz(av), 'y', ny+2, id_dim_y_yz(av) )
-          CALL handle_netcdf_error( 'netcdf', 189 )
-
-          nc_stat = NF90_DEF_VAR( id_set_yz(av), 'y', NF90_DOUBLE, &
-                                  id_dim_y_yz(av), id_var_y_yz(av) )
-          CALL handle_netcdf_error( 'netcdf', 190 )
-
-          nc_stat = NF90_PUT_ATT( id_set_yz(av), id_var_y_yz(av), 'units', &
-                                  'meters' )
-          CALL handle_netcdf_error( 'netcdf', 191 )
-
+          CALL netcdf_create_dim( id_set_yz(av), 'y', ny+2, id_dim_y_yz(av),   &
+                                  189 )
+          CALL netcdf_create_var( id_set_yz(av), (/ id_dim_y_yz(av) /), 'y',   &
+                                  NF90_DOUBLE, id_var_y_yz(av), 'meters', '',  &
+                                  190, 191, 000 )
 !
 !--       Define y-axis (for v position)
-          nc_stat = NF90_DEF_DIM( id_set_yz(av), 'yv', ny+2, id_dim_yv_yz(av) )
-          CALL handle_netcdf_error( 'netcdf', 380 )
-
-          nc_stat = NF90_DEF_VAR( id_set_yz(av), 'yv', NF90_DOUBLE, &
-                                  id_dim_yv_yz(av), id_var_yv_yz(av) )
-          CALL handle_netcdf_error( 'netcdf', 381 )
-
-          nc_stat = NF90_PUT_ATT( id_set_yz(av), id_var_yv_yz(av), 'units', &
-                                  'meters' )
-          CALL handle_netcdf_error( 'netcdf', 382 )
-
+          CALL netcdf_create_dim( id_set_yz(av), 'yv', ny+2, id_dim_yv_yz(av), &
+                                  380 )
+          CALL netcdf_create_var( id_set_yz(av), (/ id_dim_yv_yz(av) /), 'yv', &
+                                  NF90_DOUBLE, id_var_yv_yz(av), 'meters', '', &
+                                  381, 382, 000 )
 !
 !--       Define the two z-axes (zu and zw)
-          nc_stat = NF90_DEF_DIM( id_set_yz(av), 'zu', nz+2, id_dim_zu_yz(av) )
-          CALL handle_netcdf_error( 'netcdf', 192 )
-
-          nc_stat = NF90_DEF_VAR( id_set_yz(av), 'zu', NF90_DOUBLE, &
-                                  id_dim_zu_yz(av), id_var_zu_yz(av) )
-          CALL handle_netcdf_error( 'netcdf', 193 )
-
-          nc_stat = NF90_PUT_ATT( id_set_yz(av), id_var_zu_yz(av), 'units', &
-                                  'meters' )
-          CALL handle_netcdf_error( 'netcdf', 194 )
-
-          nc_stat = NF90_DEF_DIM( id_set_yz(av), 'zw', nz+2, id_dim_zw_yz(av) )
-          CALL handle_netcdf_error( 'netcdf', 195 )
-
-          nc_stat = NF90_DEF_VAR( id_set_yz(av), 'zw', NF90_DOUBLE, &
-                                  id_dim_zw_yz(av), id_var_zw_yz(av) )
-          CALL handle_netcdf_error( 'netcdf', 196 )
-
-          nc_stat = NF90_PUT_ATT( id_set_yz(av), id_var_zw_yz(av), 'units', &
-                                  'meters' )
-          CALL handle_netcdf_error( 'netcdf', 197 )
+          CALL netcdf_create_dim( id_set_yz(av), 'zu', nz+2, id_dim_zu_yz(av), &
+                                  192 )
+          CALL netcdf_create_var( id_set_yz(av), (/ id_dim_zu_yz(av) /), 'zu', &
+                                  NF90_DOUBLE, id_var_zu_yz(av), 'meters', '', &
+                                  193, 194, 000 )
+
+          CALL netcdf_create_dim( id_set_yz(av), 'zw', nz+2, id_dim_zw_yz(av), &
+                                  195 )
+          CALL netcdf_create_var( id_set_yz(av), (/ id_dim_zw_yz(av) /), 'zw', &
+                                  NF90_DOUBLE, id_var_zw_yz(av), 'meters', '', &
+                                  196, 197, 000 )
 
           IF ( land_surface )  THEN
 
-             nc_stat = NF90_DEF_DIM( id_set_yz(av), 'zs', nzs, id_dim_zs_yz(av) )
-             CALL handle_netcdf_error( 'netcdf', 545 )
-
-             nc_stat = NF90_DEF_VAR( id_set_yz(av), 'zs', NF90_DOUBLE, &
-                                     id_dim_zs_yz(av), id_var_zs_yz(av) )
-             CALL handle_netcdf_error( 'netcdf', 546 )
-
-             nc_stat = NF90_PUT_ATT( id_set_yz(av), id_var_zs_yz(av), 'units', &
-                                     'meters' )
-             CALL handle_netcdf_error( 'netcdf', 547 )
-
-          ENDIF
-
+             CALL netcdf_create_dim( id_set_yz(av), 'zs', nzs,                 &
+                                     id_dim_zs_yz(av), 545 )
+             CALL netcdf_create_var( id_set_yz(av), (/ id_dim_zs_yz(av) /),    &
+                                     'zs', NF90_DOUBLE, id_var_zs_yz(av),      &
+                                     'meters', '', 546, 547, 000 )
+
+          ENDIF
 
 !
@@ -3290,4 +3128,11 @@
                                                     grid_x, grid_y, grid_z )
 
+                      IF ( .NOT. found )  THEN
+                         WRITE ( message_string, * ) 'no grid defined for',    &
+                                                ' variable ', TRIM( do2d(av,i) )
+                         CALL message( 'define_netcdf_header', 'PA0244',       &
+                                       0, 1, 0, 6, 0 )
+                      ENDIF
+
                 END SELECT
 
@@ -3316,31 +3161,10 @@
 !
 !--             Define the grid
-                nc_stat = NF90_DEF_VAR( id_set_yz(av), do2d(av,i),             &
-                                        nc_precision(3),                       &
-                                   (/ id_x, id_y, id_z, id_dim_time_yz(av) /), &
-                                        id_var_do2d(av,i) )
-
-                IF ( .NOT. found )  THEN
-                   WRITE ( message_string, * ) 'no grid defined for', &
-                                            ' variable ', TRIM( do2d(av,i) )
-                   CALL message( 'define_netcdf_header', 'PA0244',&
-                                                                 0, 1, 0, 6, 0 )
-                ENDIF
-
-                var_list = TRIM( var_list ) // TRIM( do2d(av,i) ) // ';'
-
-                CALL handle_netcdf_error( 'netcdf', 198 )
-!
-!--             Store the 'real' name of the variable (with *, for example)
-!--             in the long_name attribute. This is evaluated by Ferret,
-!--             for example.
-                nc_stat = NF90_PUT_ATT( id_set_yz(av), id_var_do2d(av,i), &
-                                        'long_name', do2d(av,i) )
-                CALL handle_netcdf_error( 'netcdf', 199 )
-!
-!--             Define the variable's unit
-                nc_stat = NF90_PUT_ATT( id_set_yz(av), id_var_do2d(av,i), &
-                                        'units', TRIM( do2d_unit(av,i) ) )
-                CALL handle_netcdf_error( 'netcdf', 356 )
+                CALL netcdf_create_var( id_set_yz(av),  (/ id_x, id_y, id_z,   &
+                                        id_dim_time_yz(av) /), do2d(av,i),     &
+                                        nc_precision(3), id_var_do2d(av,i),    &
+                                        TRIM( do2d_unit(av,i) ), do2d(av,i),   &
+                                        198, 199, 356 )
+
 #if defined( __netcdf4_parallel )
                 IF ( netcdf_data_format > 4 )  THEN
@@ -3350,5 +3174,5 @@
                                                 id_var_do2d(av,i), &
                                                 1, 0 )
-                   CALL handle_netcdf_error( 'netcdf', 535 )
+                   CALL netcdf_handle_error( 'netcdf_define_header', 535 )
 !
 !--                Set independent io operations for parallel io. Collective io
@@ -3372,7 +3196,9 @@
 !                                                     NF90_INDEPENDENT )
                    ENDIF
-                   CALL handle_netcdf_error( 'netcdf', 450 )
+                   CALL netcdf_handle_error( 'netcdf_define_header', 450 )
                 ENDIF
 #endif
+                var_list = TRIM( var_list ) // TRIM( do2d(av,i) ) // ';'
+
              ENDIF
 
@@ -3390,5 +3216,5 @@
           nc_stat = NF90_PUT_ATT( id_set_yz(av), NF90_GLOBAL, 'VAR_LIST', &
                                   var_list )
-          CALL handle_netcdf_error( 'netcdf', 200 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 200 )
 
 !
@@ -3396,10 +3222,10 @@
 !--       parallel output.
           nc_stat = NF90_SET_FILL( id_set_yz(av), NF90_NOFILL, oldmode )
-          CALL handle_netcdf_error( 'netcdf', 531 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 531 )
 
 !
 !--       Leave netCDF define mode
           nc_stat = NF90_ENDDEF( id_set_yz(av) )
-          CALL handle_netcdf_error( 'netcdf', 201 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 201 )
 
 !
@@ -3424,5 +3250,5 @@
                                      netcdf_data, start = (/ 1 /),   &
                                      count = (/ ns /) )
-             CALL handle_netcdf_error( 'netcdf', 202 ) 
+             CALL netcdf_handle_error( 'netcdf_define_header', 202 )
 
 !
@@ -3438,5 +3264,5 @@
                                      netcdf_data, start = (/ 1 /),    &
                                      count = (/ ns /) )
-             CALL handle_netcdf_error( 'netcdf', 383 ) 
+             CALL netcdf_handle_error( 'netcdf_define_header', 383 )
 
 !
@@ -3446,5 +3272,5 @@
                                      netcdf_data, start = (/ 1 /),       &
                                      count = (/ ns /) )
-             CALL handle_netcdf_error( 'netcdf', 203 ) 
+             CALL netcdf_handle_error( 'netcdf_define_header', 203 )
 
              DEALLOCATE( netcdf_data )
@@ -3461,5 +3287,5 @@
                                      netcdf_data, start = (/ 1 /),   &
                                      count = (/ ny+2 /) )
-             CALL handle_netcdf_error( 'netcdf', 204 ) 
+             CALL netcdf_handle_error( 'netcdf_define_header', 204 )
 
              DO  j = 0, ny+1
@@ -3470,5 +3296,5 @@
                                      netcdf_data, start = (/ 1 /),    &
                                      count = (/ ny+2 /) )
-             CALL handle_netcdf_error( 'netcdf', 384 ) 
+             CALL netcdf_handle_error( 'netcdf_define_header', 384 )
 
              DEALLOCATE( netcdf_data )
@@ -3482,5 +3308,5 @@
                                      netcdf_data, start = (/ 1 /),    &
                                      count = (/ nz+2 /) )
-             CALL handle_netcdf_error( 'netcdf', 205 ) 
+             CALL netcdf_handle_error( 'netcdf_define_header', 205 )
 
              netcdf_data(0:nz+1) = zw(nzb:nzt+1)
@@ -3488,5 +3314,5 @@
                                      netcdf_data, start = (/ 1 /),    &
                                      count = (/ nz+2 /) )
-             CALL handle_netcdf_error( 'netcdf', 206 ) 
+             CALL netcdf_handle_error( 'netcdf_define_header', 206 )
 
              DEALLOCATE( netcdf_data )
@@ -3504,5 +3330,5 @@
           nc_stat = NF90_GET_ATT( id_set_yz(av), NF90_GLOBAL, 'VAR_LIST', &
                                   var_list_old )
-          CALL handle_netcdf_error( 'netcdf', 207 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 207 )
 
           var_list = ';'
@@ -3543,14 +3369,14 @@
 !--       Get and compare the number of vertical cross sections
           nc_stat = NF90_INQ_VARID( id_set_yz(av), 'x_yz', id_var_x_yz(av) )
-          CALL handle_netcdf_error( 'netcdf', 208 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 208 )
 
           nc_stat = NF90_INQUIRE_VARIABLE( id_set_yz(av), id_var_x_yz(av), &
                                            dimids = id_dim_x_yz_old )
-          CALL handle_netcdf_error( 'netcdf', 209 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 209 )
           id_dim_x_yz(av) = id_dim_x_yz_old(1)
 
           nc_stat = NF90_INQUIRE_DIMENSION( id_set_yz(av), id_dim_x_yz(av), &
                                             len = ns_old )
-          CALL handle_netcdf_error( 'netcdf', 210 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 210 )
 
           IF ( ns /= ns_old )  THEN
@@ -3571,5 +3397,5 @@
 
           nc_stat = NF90_GET_VAR( id_set_yz(av), id_var_x_yz(av), netcdf_data )
-          CALL handle_netcdf_error( 'netcdf', 211 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 211 )
 
           DO  i = 1, ns
@@ -3611,14 +3437,14 @@
 !--       on the file.
           nc_stat = NF90_INQ_VARID( id_set_yz(av), 'time', id_var_time_yz(av) )
-          CALL handle_netcdf_error( 'netcdf', 212 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 212 )
 
           nc_stat = NF90_INQUIRE_VARIABLE( id_set_yz(av), id_var_time_yz(av), &
                                            dimids = id_dim_time_old )
-          CALL handle_netcdf_error( 'netcdf', 213 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 213 )
           id_dim_time_yz(av) = id_dim_time_old(1)
 
           nc_stat = NF90_INQUIRE_DIMENSION( id_set_yz(av), id_dim_time_yz(av), &
                                             len = ntime_count )
-          CALL handle_netcdf_error( 'netcdf', 214 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 214 )
 
 !
@@ -3633,5 +3459,5 @@
                                   start = (/ do2d_yz_time_count(av) /), &
                                   count = (/ 1 /) )
-          CALL handle_netcdf_error( 'netcdf', 215 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 215 )
 
           IF ( last_time_coordinate(1) >= simulated_time )  THEN
@@ -3704,5 +3530,5 @@
                 nc_stat = NF90_INQ_VARID( id_set_yz(av), do2d(av,i), &
                                           id_var_do2d(av,i) )
-                CALL handle_netcdf_error( 'netcdf', 216 )
+                CALL netcdf_handle_error( 'netcdf_define_header', 216 )
 #if defined( __netcdf4_parallel )
 !
@@ -3728,5 +3554,5 @@
 !                                                     NF90_INDEPENDENT )
                    ENDIF
-                   CALL handle_netcdf_error( 'netcdf', 450 )
+                   CALL netcdf_handle_error( 'netcdf_define_header', 450 )
                 ENDIF
 #endif
@@ -3749,11 +3575,11 @@
           ENDIF
           nc_stat = NF90_REDEF( id_set_yz(av) )
-          CALL handle_netcdf_error( 'netcdf', 435 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 435 )
           nc_stat = NF90_PUT_ATT( id_set_yz(av), NF90_GLOBAL, 'title', &
                                   TRIM( run_description_header ) //    &
                                   TRIM( time_average_text ) )
-          CALL handle_netcdf_error( 'netcdf', 217 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 217 )
           nc_stat = NF90_ENDDEF( id_set_yz(av) )
-          CALL handle_netcdf_error( 'netcdf', 436 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 436 )
           message_string = 'netCDF file for cross-sections ' //           &
                             TRIM( var ) // ' from previous run found.' // &
@@ -3772,5 +3598,5 @@
                                      TRIM( run_description_header ) //  &
                                      TRIM( time_average_text ) )
-             CALL handle_netcdf_error( 'netcdf', 218 )
+             CALL netcdf_handle_error( 'netcdf_define_header', 218 )
 
              WRITE ( time_average_text,'(F7.1,'' s avg'')' ) averaging_interval_pr
@@ -3781,5 +3607,5 @@
                                      TRIM( run_description_header ) )
           ENDIF
-          CALL handle_netcdf_error( 'netcdf', 219 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 219 )
 
 !
@@ -3790,10 +3616,10 @@
                                   'no_rows',                                  & 
                                   profile_rows ) 
-          CALL handle_netcdf_error( 'netcdf', 519 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 519 )
 
           nc_stat = NF90_PUT_ATT( id_set_pr, NF90_GLOBAL,                     &
                                   'no_columns',                               & 
                                   profile_columns ) 
-          CALL handle_netcdf_error( 'netcdf', 520 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 520 )
 
 
@@ -3831,5 +3657,4 @@
 !--       Check if all profiles defined in cross_profiles are defined in 
 !--       data_output_pr. If not, they will be skipped.
-
           DO  i = 1, cross_profiles_count
              DO  j = 1, dopr_n
@@ -3869,5 +3694,4 @@
 !--       Check if all profiles defined in data_output_pr are defined in 
 !--       cross_profiles. If not, they will be added to cross_profiles.
-
           DO  i = 1, dopr_n
              DO  j = 1, cross_profiles_count
@@ -3929,19 +3753,13 @@
           nc_stat = NF90_PUT_ATT( id_set_pr, NF90_GLOBAL, 'cross_profiles',   &
                                   TRIM( char_cross_profiles ) )
-          CALL handle_netcdf_error( 'netcdf', 521 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 521 )
 
 !
 !--       Define time coordinate for profiles (unlimited dimension)
-          nc_stat = NF90_DEF_DIM( id_set_pr, 'time', NF90_UNLIMITED, &
-                                  id_dim_time_pr )
-          CALL handle_netcdf_error( 'netcdf', 220 )
-
-          nc_stat = NF90_DEF_VAR( id_set_pr, 'time', NF90_DOUBLE, &
-                                  id_dim_time_pr, id_var_time_pr )
-          CALL handle_netcdf_error( 'netcdf', 221 )
-
-          nc_stat = NF90_PUT_ATT( id_set_pr, id_var_time_pr, 'units', 'seconds')
-          CALL handle_netcdf_error( 'netcdf', 222 )
-
+          CALL netcdf_create_dim( id_set_pr, 'time', NF90_UNLIMITED,           &
+                                  id_dim_time_pr, 220 )
+          CALL netcdf_create_var( id_set_pr, (/ id_dim_time_pr /), 'time',     &
+                                  NF90_DOUBLE, id_var_time_pr, 'seconds', '',  &
+                                  221, 222, 000 )
 !
 !--       Define the variables
@@ -3953,32 +3771,19 @@
 !
 !--             Define the z-axes (each variable gets its own z-axis)
-                nc_stat = NF90_DEF_DIM( id_set_pr,                        &
-                                        'z' // TRIM( data_output_pr(i) ), &
-                                        nzt+2-nzb, id_dim_z_pr(i,0) )
-                CALL handle_netcdf_error( 'netcdf', 223 )
-
-                nc_stat = NF90_DEF_VAR( id_set_pr,                        &
-                                        'z' // TRIM( data_output_pr(i) ), &
-                                        NF90_DOUBLE, id_dim_z_pr(i,0),    &
-                                        id_var_z_pr(i,0) )
-                CALL handle_netcdf_error( 'netcdf', 224 )
-
-                nc_stat = NF90_PUT_ATT( id_set_pr, id_var_z_pr(i,0), 'units', &
-                                        'meters' )
-                CALL handle_netcdf_error( 'netcdf', 225 )
+                CALL netcdf_create_dim( id_set_pr,                             &
+                                        'z' // TRIM( data_output_pr(i) ),      &
+                                        nzt+2-nzb, id_dim_z_pr(i,0), 223 )
+                CALL netcdf_create_var( id_set_pr, (/ id_dim_z_pr(i,0) /),     &
+                                        'z' // TRIM( data_output_pr(i) ),      &
+                                       NF90_DOUBLE, id_var_z_pr(i,0),          &
+                                       'meters', '', 224, 225, 000 )
 !
 !--             Define the variable
-                nc_stat = NF90_DEF_VAR( id_set_pr, data_output_pr(i),         &
-                                        nc_precision(5), (/ id_dim_z_pr(i,0), &
-                                        id_dim_time_pr /), id_var_dopr(i,0) )
-                CALL handle_netcdf_error( 'netcdf', 226 )
-
-                nc_stat = NF90_PUT_ATT( id_set_pr, id_var_dopr(i,0), &
-                                        'long_name', TRIM( data_output_pr(i) ) )
-                CALL handle_netcdf_error( 'netcdf', 227 )
-
-                nc_stat = NF90_PUT_ATT( id_set_pr, id_var_dopr(i,0), &
-                                        'units', TRIM( dopr_unit(i) ) )
-                CALL handle_netcdf_error( 'netcdf', 228 )
+                CALL netcdf_create_var( id_set_pr, (/ id_dim_z_pr(i,0),        &
+                                        id_dim_time_pr /), data_output_pr(i),  &
+                                        nc_precision(5), id_var_dopr(i,0),     &
+                                        TRIM( dopr_unit(i) ),                  &
+                                        TRIM( data_output_pr(i) ), 226, 227,   &
+                                        228 )
 
                 var_list = TRIM( var_list ) // TRIM( data_output_pr(i) ) //  ';'
@@ -3993,37 +3798,21 @@
 !
 !--                Define the z-axes (each variable gets it own z-axis)
-                   nc_stat = NF90_DEF_DIM( id_set_pr,                          &
-                                         'z'//TRIM(data_output_pr(i))//suffix, &
-                                           nzt+2-nzb, id_dim_z_pr(i,j) )
-                   CALL handle_netcdf_error( 'netcdf', 229 )
-
-                   nc_stat = NF90_DEF_VAR( id_set_pr,                          &
-                                         'z'//TRIM(data_output_pr(i))//suffix, &
-                                           nc_precision(5), id_dim_z_pr(i,j),  &
-                                           id_var_z_pr(i,j) )
-                   CALL handle_netcdf_error( 'netcdf', 230 )
-
-                   nc_stat = NF90_PUT_ATT( id_set_pr, id_var_z_pr(i,j), &
-                                           'units', 'meters' )
-                   CALL handle_netcdf_error( 'netcdf', 231 )
-
+                   CALL netcdf_create_dim( id_set_pr, 'z' //                   &
+                                           TRIM(data_output_pr(i)) // suffix,  &
+                                           nzt+2-nzb, id_dim_z_pr(i,j), 229 )
+                   CALL netcdf_create_var( id_set_pr, (/ id_dim_z_pr(i,j) /),  &
+                                           'z' // TRIM(data_output_pr(i)) //   &
+                                           suffix, NF90_DOUBLE,                &
+                                           id_var_z_pr(i,j), 'meters', '',     &
+                                           230, 231, 000 )
 !
 !--                Define the variable
-                   nc_stat = NF90_DEF_VAR( id_set_pr,                         &
-                                           TRIM(data_output_pr(i)) // suffix, &
-                                           nc_precision(5),                   &
-                                           (/ id_dim_z_pr(i,j),               &
-                                           id_dim_time_pr /), id_var_dopr(i,j) )
-                   CALL handle_netcdf_error( 'netcdf', 232 )
-
-                   nc_stat = NF90_PUT_ATT( id_set_pr, id_var_dopr(i,j),        &
-                                           'long_name',                        &
-                                           TRIM( data_output_pr(i) ) // ' SR ' &
-                                           // suffix(2:2) )
-                   CALL handle_netcdf_error( 'netcdf', 233 )
-
-                   nc_stat = NF90_PUT_ATT( id_set_pr, id_var_dopr(i,j), &
-                                           'units', TRIM( dopr_unit(i) ) )
-                   CALL handle_netcdf_error( 'netcdf', 234 )
+                   CALL netcdf_create_var( id_set_pr, (/ id_dim_z_pr(i,j),     &
+                                           id_dim_time_pr /),                  &
+                                           TRIM(data_output_pr(i)) // suffix,  &
+                                           nc_precision(5), id_var_dopr(i,j),  &
+                                           TRIM( dopr_unit(i) ),               &
+                                           TRIM( data_output_pr(i) ) //        &
+                                           ' SR ', 232, 233, 234 )
 
                    var_list = TRIM( var_list ) // TRIM( data_output_pr(i) ) // &
@@ -4040,5 +3829,5 @@
 !--       restart runs)
           nc_stat = NF90_PUT_ATT( id_set_pr, NF90_GLOBAL, 'VAR_LIST', var_list )
-          CALL handle_netcdf_error( 'netcdf', 235 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 235 )
 
 !
@@ -4046,14 +3835,9 @@
           DO  i = 1, dopr_norm_num
 
-             nc_stat = NF90_DEF_VAR( id_set_pr, 'NORM_' // &
-                                     TRIM( dopr_norm_names(i) ), &
-                                     nc_precision(5), (/ id_dim_time_pr /), &
-                                     id_var_norm_dopr(i) )
-             CALL handle_netcdf_error( 'netcdf', 236 )
-
-             nc_stat = NF90_PUT_ATT( id_set_pr, id_var_norm_dopr(i), &
-                                     'long_name',                    &
-                                     TRIM( dopr_norm_longnames(i) ) )
-             CALL handle_netcdf_error( 'netcdf', 237 )
+             CALL netcdf_create_var( id_set_pr, (/ id_dim_time_pr /),          &
+                                     'NORM_' // TRIM( dopr_norm_names(i) ),    &
+                                     nc_precision(5), id_var_norm_dopr(i),     &
+                                     '', TRIM( dopr_norm_longnames(i) ), 236,  &
+                                     000, 237 )
 
           ENDDO
@@ -4062,5 +3846,5 @@
 !--       Leave netCDF define mode
           nc_stat = NF90_ENDDEF( id_set_pr )
-          CALL handle_netcdf_error( 'netcdf', 238 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 238 )
 
 !
@@ -4073,5 +3857,5 @@
                                         start = (/ 1 /),                  &
                                         count = (/ nzt-nzb+2 /) )
-                CALL handle_netcdf_error( 'netcdf', 239 )
+                CALL netcdf_handle_error( 'netcdf_define_header', 239 )
 
              ENDDO
@@ -4088,5 +3872,5 @@
           nc_stat = NF90_GET_ATT( id_set_pr, NF90_GLOBAL, 'VAR_LIST', &
                                   var_list_old )
-          CALL handle_netcdf_error( 'netcdf', 240 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 240 )
 
           var_list = ';'
@@ -4122,14 +3906,14 @@
 !--       on the file.
           nc_stat = NF90_INQ_VARID( id_set_pr, 'time', id_var_time_pr )
-          CALL handle_netcdf_error( 'netcdf', 241 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 241 )
 
           nc_stat = NF90_INQUIRE_VARIABLE( id_set_pr, id_var_time_pr, &
                                            dimids = id_dim_time_old )
-          CALL handle_netcdf_error( 'netcdf', 242 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 242 )
           id_dim_time_pr = id_dim_time_old(1)
 
           nc_stat = NF90_INQUIRE_DIMENSION( id_set_pr, id_dim_time_pr, &
                                             len = dopr_time_count )
-          CALL handle_netcdf_error( 'netcdf', 243 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 243 )
 
           nc_stat = NF90_GET_VAR( id_set_pr, id_var_time_pr,        &
@@ -4137,5 +3921,5 @@
                                   start = (/ dopr_time_count /), &
                                   count = (/ 1 /) )
-          CALL handle_netcdf_error( 'netcdf', 244 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 244 )
 
           IF ( last_time_coordinate(1) >= simulated_time )  THEN
@@ -4162,5 +3946,5 @@
                 nc_stat = NF90_INQ_VARID( id_set_pr, data_output_pr(i), &
                                           id_var_dopr(i,0) )
-                CALL handle_netcdf_error( 'netcdf', 245 )
+                CALL netcdf_handle_error( 'netcdf_define_header', 245 )
              ELSE
                 DO  j = 0, statistic_regions
@@ -4169,5 +3953,5 @@
                    nc_stat = NF90_INQ_VARID( id_set_pr, netcdf_var_name, &
                                              id_var_dopr(i,j) )
-                   CALL handle_netcdf_error( 'netcdf', 246 )
+                   CALL netcdf_handle_error( 'netcdf_define_header', 246 )
                 ENDDO
              ENDIF
@@ -4181,5 +3965,5 @@
                                        'NORM_' // TRIM( dopr_norm_names(i) ), &
                                        id_var_norm_dopr(i) )
-             CALL handle_netcdf_error( 'netcdf', 247 )
+             CALL netcdf_handle_error( 'netcdf_define_header', 247 )
           ENDDO
 
@@ -4198,12 +3982,12 @@
           ENDIF
           nc_stat = NF90_REDEF( id_set_pr )
-          CALL handle_netcdf_error( 'netcdf', 437 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 437 )
           nc_stat = NF90_PUT_ATT( id_set_pr, NF90_GLOBAL, 'title', &
                                   TRIM( run_description_header ) //    &
                                   TRIM( time_average_text ) )
-          CALL handle_netcdf_error( 'netcdf', 248 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 248 )
 
           nc_stat = NF90_ENDDEF( id_set_pr )
-          CALL handle_netcdf_error( 'netcdf', 438 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 438 )
           message_string = 'netCDF file for vertical profiles ' // &
                            'from previous run found.' //           &
@@ -4218,19 +4002,13 @@
           nc_stat = NF90_PUT_ATT( id_set_ts, NF90_GLOBAL, 'title', &
                                   TRIM( run_description_header ) )
-          CALL handle_netcdf_error( 'netcdf', 249 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 249 )
 
 !
 !--       Define time coordinate for time series (unlimited dimension)
-          nc_stat = NF90_DEF_DIM( id_set_ts, 'time', NF90_UNLIMITED, &
-                                  id_dim_time_ts )
-          CALL handle_netcdf_error( 'netcdf', 250 )
-
-          nc_stat = NF90_DEF_VAR( id_set_ts, 'time', NF90_DOUBLE, &
-                                  id_dim_time_ts, id_var_time_ts )
-          CALL handle_netcdf_error( 'netcdf', 251 )
-
-          nc_stat = NF90_PUT_ATT( id_set_ts, id_var_time_ts, 'units', 'seconds')
-          CALL handle_netcdf_error( 'netcdf', 252 )
-
+          CALL netcdf_create_dim( id_set_ts, 'time', NF90_UNLIMITED,           &
+                                  id_dim_time_ts, 250 )
+          CALL netcdf_create_var( id_set_ts, (/ id_dim_time_ts /), 'time',     &
+                                  NF90_DOUBLE, id_var_time_ts, 'seconds', '',  &
+                                  251, 252, 000 )
 !
 !--       Define the variables
@@ -4240,18 +4018,9 @@
              IF ( statistic_regions == 0 )  THEN
 
-                nc_stat = NF90_DEF_VAR( id_set_ts, dots_label(i),              &
-                                        nc_precision(6), (/ id_dim_time_ts /), &
-                                        id_var_dots(i,0) )
-                CALL handle_netcdf_error( 'netcdf', 253 )
-
-                nc_stat = NF90_PUT_ATT( id_set_ts, id_var_dots(i,0), &
-                                        'long_name', TRIM( dots_label(i) ) )
-                CALL handle_netcdf_error( 'netcdf', 254 )
-
-                nc_stat = NF90_PUT_ATT( id_set_ts, id_var_dots(i,0), &
-                                        'units', TRIM( dots_unit(i) ) )
-                CALL handle_netcdf_error( 'netcdf', 255 )
-
-                var_list = TRIM( var_list ) // TRIM( dots_label(i) ) // ';'
+                CALL netcdf_create_var( id_set_ts, (/ id_dim_time_ts /),       &
+                                        dots_label(i), nc_precision(6),        &
+                                        id_var_dots(i,0),                      &
+                                        TRIM( dots_unit(i) ),                  &
+                                        TRIM( dots_label(i) ), 253, 254, 255 )
 
              ELSE
@@ -4262,23 +4031,10 @@
                    WRITE ( suffix, '(''_'',I1)' )  j
 
-                   nc_stat = NF90_DEF_VAR( id_set_ts,                         &
-                                           TRIM( dots_label(i) ) // suffix,   &
-                                           nc_precision(6),                   &
-                                           (/ id_dim_time_ts /),              &
-                                           id_var_dots(i,j) )
-                   CALL handle_netcdf_error( 'netcdf', 256 )
-
-                   nc_stat = NF90_PUT_ATT( id_set_ts, id_var_dots(i,j),       &
-                                           'long_name',                       &
-                                           TRIM( dots_label(i) ) // ' SR ' // &
-                                           suffix(2:2) )
-                   CALL handle_netcdf_error( 'netcdf', 257 )
-
-                   nc_stat = NF90_PUT_ATT( id_set_ts, id_var_dots(i,j), &
-                                           'units', TRIM( dots_unit(i) ) )
-                   CALL handle_netcdf_error( 'netcdf', 347 )
-
-                   var_list = TRIM( var_list ) // TRIM( dots_label(i) ) // &
-                              suffix // ';'
+                   CALL netcdf_create_var( id_set_ts, (/ id_dim_time_ts /),    &
+                                           TRIM( dots_label(i) ) // suffix,    &
+                                           nc_precision(6), id_var_dots(i,j),  &
+                                           TRIM( dots_unit(i) ),               &
+                                           TRIM( dots_label(i) ) // ' SR ' //  &
+                                           suffix(2:2), 256, 257, 347)
 
                 ENDDO
@@ -4292,10 +4048,10 @@
 !--       restart runs)
           nc_stat = NF90_PUT_ATT( id_set_ts, NF90_GLOBAL, 'VAR_LIST', var_list )
-          CALL handle_netcdf_error( 'netcdf', 258 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 258 )
 
 !
 !--       Leave netCDF define mode
           nc_stat = NF90_ENDDEF( id_set_ts )
-          CALL handle_netcdf_error( 'netcdf', 259 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 259 )
 
 
@@ -4309,5 +4065,5 @@
           nc_stat = NF90_GET_ATT( id_set_ts, NF90_GLOBAL, 'VAR_LIST', &
                                   var_list_old )
-          CALL handle_netcdf_error( 'netcdf', 260 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 260 )
 
           var_list = ';'
@@ -4344,14 +4100,14 @@
 !--       on the file.
           nc_stat = NF90_INQ_VARID( id_set_ts, 'time', id_var_time_ts )
-          CALL handle_netcdf_error( 'netcdf', 261 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 261 )
 
           nc_stat = NF90_INQUIRE_VARIABLE( id_set_ts, id_var_time_ts, &
                                            dimids = id_dim_time_old )
-          CALL handle_netcdf_error( 'netcdf', 262 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 262 )
           id_dim_time_ts = id_dim_time_old(1)
 
           nc_stat = NF90_INQUIRE_DIMENSION( id_set_ts, id_dim_time_ts, &
                                             len = dots_time_count )
-          CALL handle_netcdf_error( 'netcdf', 263 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 263 )
 
           nc_stat = NF90_GET_VAR( id_set_ts, id_var_time_ts,        &
@@ -4359,5 +4115,5 @@
                                   start = (/ dots_time_count /), &
                                   count = (/ 1 /) )
-          CALL handle_netcdf_error( 'netcdf', 264 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 264 )
 
           IF ( last_time_coordinate(1) >= simulated_time )  THEN
@@ -4384,5 +4140,5 @@
                 nc_stat = NF90_INQ_VARID( id_set_ts, dots_label(i), &
                                           id_var_dots(i,0) )
-                CALL handle_netcdf_error( 'netcdf', 265 )
+                CALL netcdf_handle_error( 'netcdf_define_header', 265 )
              ELSE
                 DO  j = 0, statistic_regions
@@ -4391,5 +4147,5 @@
                    nc_stat = NF90_INQ_VARID( id_set_ts, netcdf_var_name, &
                                              id_var_dots(i,j) )
-                   CALL handle_netcdf_error( 'netcdf', 266 )
+                   CALL netcdf_handle_error( 'netcdf_define_header', 266 )
                 ENDDO
              ENDIF
@@ -4405,10 +4161,10 @@
 !--       to ensure equal attribute size in a job chain. Maybe revise later.
           nc_stat = NF90_REDEF( id_set_ts )
-          CALL handle_netcdf_error( 'netcdf', 439 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 439 )
           nc_stat = NF90_PUT_ATT( id_set_ts, NF90_GLOBAL, 'title', &
                                   TRIM( run_description_header ) )
-          CALL handle_netcdf_error( 'netcdf', 267 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 267 )
           nc_stat = NF90_ENDDEF( id_set_ts )
-          CALL handle_netcdf_error( 'netcdf', 440 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 440 )
           message_string = 'netCDF file for time series ' // &
                            'from previous run found.' //     &
@@ -4427,5 +4183,5 @@
                                      TRIM( run_description_header ) // &
                                      TRIM( time_average_text ) )
-             CALL handle_netcdf_error( 'netcdf', 268 )
+             CALL netcdf_handle_error( 'netcdf_define_header', 268 )
 
              WRITE ( time_average_text,'(F7.1,'' s avg'')' )  averaging_interval_sp
@@ -4436,19 +4192,13 @@
                                      TRIM( run_description_header ) )
           ENDIF
-          CALL handle_netcdf_error( 'netcdf', 269 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 269 )
 
 !
 !--       Define time coordinate for spectra (unlimited dimension)
-          nc_stat = NF90_DEF_DIM( id_set_sp, 'time', NF90_UNLIMITED, &
-                                  id_dim_time_sp )
-          CALL handle_netcdf_error( 'netcdf', 270 )
-
-          nc_stat = NF90_DEF_VAR( id_set_sp, 'time', NF90_DOUBLE, &
-                                  id_dim_time_sp, id_var_time_sp )
-          CALL handle_netcdf_error( 'netcdf', 271 )
-
-          nc_stat = NF90_PUT_ATT( id_set_sp, id_var_time_sp, 'units', 'seconds')
-          CALL handle_netcdf_error( 'netcdf', 272 )
-
+          CALL netcdf_create_dim( id_set_sp, 'time', NF90_UNLIMITED,           &
+                                  id_dim_time_sp, 270 )
+          CALL netcdf_create_var( id_set_sp, (/ id_dim_time_sp /), 'time',     &
+                                  NF90_DOUBLE, id_var_time_sp, 'seconds', '',  &
+                                  271, 272, 000 )
 !
 !--       Define the spatial dimensions and coordinates for spectra.
@@ -4463,50 +4213,26 @@
 !
 !--       Define vertical coordinate grid (zu grid)
-          nc_stat = NF90_DEF_DIM( id_set_sp, 'zu_sp', ns, id_dim_zu_sp )
-          CALL handle_netcdf_error( 'netcdf', 273 )
-
-          nc_stat = NF90_DEF_VAR( id_set_sp, 'zu_sp', NF90_DOUBLE, &
-                                  id_dim_zu_sp, id_var_zu_sp )
-          CALL handle_netcdf_error( 'netcdf', 274 )
-
-          nc_stat = NF90_PUT_ATT( id_set_sp, id_var_zu_sp, 'units', 'meters' )
-          CALL handle_netcdf_error( 'netcdf', 275 )
-
+          CALL netcdf_create_dim( id_set_sp, 'zu_sp', ns, id_dim_zu_sp, 273 )
+          CALL netcdf_create_var( id_set_sp, (/ id_dim_zu_sp /), 'zu_sp',      &
+                                  NF90_DOUBLE, id_var_zu_sp, 'meters', '',     &
+                                  274, 275, 000 )
 !
 !--       Define vertical coordinate grid (zw grid)
-          nc_stat = NF90_DEF_DIM( id_set_sp, 'zw_sp', ns, id_dim_zw_sp )
-          CALL handle_netcdf_error( 'netcdf', 276 )
-
-          nc_stat = NF90_DEF_VAR( id_set_sp, 'zw_sp', NF90_DOUBLE, &
-                                  id_dim_zw_sp, id_var_zw_sp )
-          CALL handle_netcdf_error( 'netcdf', 277 )
-
-          nc_stat = NF90_PUT_ATT( id_set_sp, id_var_zw_sp, 'units', 'meters' )
-          CALL handle_netcdf_error( 'netcdf', 278 )
-
+          CALL netcdf_create_dim( id_set_sp, 'zw_sp', ns, id_dim_zw_sp, 276 )
+          CALL netcdf_create_var( id_set_sp, (/ id_dim_zw_sp /), 'zw_sp',      &
+                                  NF90_DOUBLE, id_var_zw_sp, 'meters', '',     &
+                                  277, 278, 000 )
 !
 !--       Define x-axis
-          nc_stat = NF90_DEF_DIM( id_set_sp, 'k_x', nx/2, id_dim_x_sp )
-          CALL handle_netcdf_error( 'netcdf', 279 )
-
-          nc_stat = NF90_DEF_VAR( id_set_sp, 'k_x', NF90_DOUBLE, id_dim_x_sp, &
-                                  id_var_x_sp )
-          CALL handle_netcdf_error( 'netcdf', 280 )
-
-          nc_stat = NF90_PUT_ATT( id_set_sp, id_var_x_sp, 'units', 'm-1' )
-          CALL handle_netcdf_error( 'netcdf', 281 )
-
+          CALL netcdf_create_dim( id_set_sp, 'k_x', nx/2, id_dim_x_sp, 279 )
+          CALL netcdf_create_var( id_set_sp, (/ id_dim_x_sp /), 'k_x',         &
+                                  NF90_DOUBLE, id_var_x_sp, 'm-1', '', 280,    &
+                                  281, 000 )
 !
 !--       Define y-axis
-          nc_stat = NF90_DEF_DIM( id_set_sp, 'k_y', ny/2, id_dim_y_sp )
-          CALL handle_netcdf_error( 'netcdf', 282 )
-
-          nc_stat = NF90_DEF_VAR( id_set_sp, 'k_y', NF90_DOUBLE, id_dim_y_sp, &
-                                  id_var_y_sp )
-          CALL handle_netcdf_error( 'netcdf', 283 )
-
-          nc_stat = NF90_PUT_ATT( id_set_sp, id_var_y_sp, 'units', 'm-1' )
-          CALL handle_netcdf_error( 'netcdf', 284 )
-
+          CALL netcdf_create_dim(id_set_sp, 'k_y', ny/2, id_dim_y_sp, 282 )
+          CALL netcdf_create_var( id_set_sp, (/ id_dim_y_sp /), 'k_y',         &
+                                  NF90_DOUBLE, id_var_y_sp, 'm-1', '', 283,    &
+                                  284, 000 )
 !
 !--       Define the variables
@@ -4546,23 +4272,16 @@
                 netcdf_var_name = TRIM( data_output_sp(i) ) // '_x'
                 IF ( TRIM( grid_z ) == 'zw' )  THEN
-                   nc_stat = NF90_DEF_VAR( id_set_sp, netcdf_var_name,      &
-                                           nc_precision(7), (/ id_dim_x_sp, &
-                                           id_dim_zw_sp, id_dim_time_sp /), &
-                                           id_var_dospx(i) )
+                   CALL netcdf_create_var( id_set_sp, (/ id_dim_x_sp,          &
+                                           id_dim_zw_sp, id_dim_time_sp /),    &
+                                           netcdf_var_name, nc_precision(7),   &
+                                           id_var_dospx(i), 'unknown',         &
+                                           netcdf_var_name, 285, 286, 287 )
                 ELSE
-                   nc_stat = NF90_DEF_VAR( id_set_sp, netcdf_var_name,      &
-                                           nc_precision(7), (/ id_dim_x_sp, &
-                                           id_dim_zu_sp, id_dim_time_sp /), &
-                                           id_var_dospx(i) )
+                   CALL netcdf_create_var( id_set_sp, (/ id_dim_x_sp,          &
+                                           id_dim_zu_sp, id_dim_time_sp /),    &
+                                           netcdf_var_name, nc_precision(7),   &
+                                           id_var_dospx(i), 'unknown',         &
+                                           netcdf_var_name, 285, 286, 287 )
                 ENDIF
-                CALL handle_netcdf_error( 'netcdf', 285 )
-
-                nc_stat = NF90_PUT_ATT( id_set_sp, id_var_dospx(i), &
-                                        'long_name', netcdf_var_name )
-                CALL handle_netcdf_error( 'netcdf', 286 )
-
-                nc_stat = NF90_PUT_ATT( id_set_sp, id_var_dospx(i), &
-                                        'units', 'unknown' )
-                CALL handle_netcdf_error( 'netcdf', 287 )
 
                 var_list = TRIM( var_list ) // TRIM( netcdf_var_name ) // ';'
@@ -4576,23 +4295,16 @@
                 netcdf_var_name = TRIM( data_output_sp(i) ) // '_y'
                 IF ( TRIM( grid_z ) == 'zw' )  THEN
-                   nc_stat = NF90_DEF_VAR( id_set_sp, netcdf_var_name,      &
-                                           nc_precision(7), (/ id_dim_y_sp, &
-                                           id_dim_zw_sp, id_dim_time_sp /), &
-                                           id_var_dospy(i) )
+                   CALL netcdf_create_var( id_set_sp, (/ id_dim_y_sp,          &
+                                           id_dim_zw_sp, id_dim_time_sp /),    &
+                                           netcdf_var_name, nc_precision(7),   &
+                                           id_var_dospy(i), 'unknown',         &
+                                           netcdf_var_name, 288, 289, 290 )
                 ELSE
-                   nc_stat = NF90_DEF_VAR( id_set_sp, netcdf_var_name,      &
-                                           nc_precision(7), (/ id_dim_y_sp, &
-                                           id_dim_zu_sp, id_dim_time_sp /), &
-                                           id_var_dospy(i) )
+                   CALL netcdf_create_var( id_set_sp, (/ id_dim_y_sp,          &
+                                           id_dim_zu_sp, id_dim_time_sp /),    &
+                                           netcdf_var_name, nc_precision(7),   &
+                                           id_var_dospy(i), 'unknown',         &
+                                           netcdf_var_name, 288, 289, 290 )
                 ENDIF
-                CALL handle_netcdf_error( 'netcdf', 288 )
-
-                nc_stat = NF90_PUT_ATT( id_set_sp, id_var_dospy(i), &
-                                        'long_name', netcdf_var_name )
-                CALL handle_netcdf_error( 'netcdf', 289 )
-
-                nc_stat = NF90_PUT_ATT( id_set_sp, id_var_dospy(i), &
-                                        'units', 'unknown' )
-                CALL handle_netcdf_error( 'netcdf', 290 )
 
                 var_list = TRIM( var_list ) // TRIM( netcdf_var_name ) // ';'
@@ -4608,10 +4320,10 @@
 !--       restart runs)
           nc_stat = NF90_PUT_ATT( id_set_sp, NF90_GLOBAL, 'VAR_LIST', var_list )
-          CALL handle_netcdf_error( 'netcdf', 291 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 291 )
 
 !
 !--       Leave netCDF define mode
           nc_stat = NF90_ENDDEF( id_set_sp )
-          CALL handle_netcdf_error( 'netcdf', 292 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 292 )
 
 !
@@ -4624,5 +4336,5 @@
           nc_stat = NF90_PUT_VAR( id_set_sp, id_var_zu_sp, netcdf_data, &
                                   start = (/ 1 /), count = (/ ns /) )
-          CALL handle_netcdf_error( 'netcdf', 293 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 293 )
 
 !
@@ -4631,5 +4343,5 @@
           nc_stat = NF90_PUT_VAR( id_set_sp, id_var_zw_sp, netcdf_data, &
                                   start = (/ 1 /), count = (/ ns /) )
-          CALL handle_netcdf_error( 'netcdf', 294 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 294 )
 
           DEALLOCATE( netcdf_data )
@@ -4644,5 +4356,5 @@
           nc_stat = NF90_PUT_VAR( id_set_sp, id_var_x_sp, netcdf_data, &
                                   start = (/ 1 /), count = (/ nx/2 /) )
-          CALL handle_netcdf_error( 'netcdf', 295 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 295 )
 
           DEALLOCATE( netcdf_data )
@@ -4655,5 +4367,5 @@
           nc_stat = NF90_PUT_VAR( id_set_sp, id_var_y_sp, netcdf_data, &
                                   start = (/ 1 /), count = (/ ny/2 /) )
-          CALL handle_netcdf_error( 'netcdf', 296 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 296 )
 
           DEALLOCATE( netcdf_data )
@@ -4669,5 +4381,5 @@
           nc_stat = NF90_GET_ATT( id_set_sp, NF90_GLOBAL, 'VAR_LIST', &
                                   var_list_old )
-          CALL handle_netcdf_error( 'netcdf', 297 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 297 )
           var_list = ';'
           i = 1
@@ -4711,14 +4423,14 @@
 !--       Get and compare the number of vertical levels
           nc_stat = NF90_INQ_VARID( id_set_sp, 'zu_sp', id_var_zu_sp )
-          CALL handle_netcdf_error( 'netcdf', 298 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 298 )
 
           nc_stat = NF90_INQUIRE_VARIABLE( id_set_sp, id_var_zu_sp, &
                                            dimids = id_dim_zu_sp_old )
-          CALL handle_netcdf_error( 'netcdf', 299 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 299 )
           id_dim_zu_sp = id_dim_zu_sp_old(1)
 
           nc_stat = NF90_INQUIRE_DIMENSION( id_set_sp, id_dim_zu_sp, &
                                             len = ns_old )
-          CALL handle_netcdf_error( 'netcdf', 300 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 300 )
 
           IF ( ns /= ns_old )  THEN
@@ -4739,5 +4451,5 @@
 
           nc_stat = NF90_GET_VAR( id_set_sp, id_var_zu_sp, netcdf_data )
-          CALL handle_netcdf_error( 'netcdf', 301 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 301 )
 
           DO  i = 1, ns
@@ -4763,14 +4475,14 @@
 !--       on the file.
           nc_stat = NF90_INQ_VARID( id_set_sp, 'time', id_var_time_sp )
-          CALL handle_netcdf_error( 'netcdf', 302 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 302 )
 
           nc_stat = NF90_INQUIRE_VARIABLE( id_set_sp, id_var_time_sp, &
                                            dimids = id_dim_time_old )
-          CALL handle_netcdf_error( 'netcdf', 303 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 303 )
           id_dim_time_sp = id_dim_time_old(1)
 
           nc_stat = NF90_INQUIRE_DIMENSION( id_set_sp, id_dim_time_sp, &
                                             len = dosp_time_count )
-          CALL handle_netcdf_error( 'netcdf', 304 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 304 )
 
           nc_stat = NF90_GET_VAR( id_set_sp, id_var_time_sp,        &
@@ -4778,5 +4490,5 @@
                                   start = (/ dosp_time_count /), &
                                   count = (/ 1 /) )
-          CALL handle_netcdf_error( 'netcdf', 305 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 305 )
 
           IF ( last_time_coordinate(1) >= simulated_time )  THEN
@@ -4804,5 +4516,5 @@
                 nc_stat = NF90_INQ_VARID( id_set_sp, netcdf_var_name, &
                                           id_var_dospx(i) )
-                CALL handle_netcdf_error( 'netcdf', 306 )
+                CALL netcdf_handle_error( 'netcdf_define_header', 306 )
              ENDIF
 
@@ -4811,5 +4523,5 @@
                 nc_stat = NF90_INQ_VARID( id_set_sp, netcdf_var_name, &
                                           id_var_dospy(i) )
-                CALL handle_netcdf_error( 'netcdf', 307 )
+                CALL netcdf_handle_error( 'netcdf_define_header', 307 )
              ENDIF
 
@@ -4826,5 +4538,5 @@
 !--       to ensure equal attribute size in a job chain. Maybe revise later.
           nc_stat = NF90_REDEF( id_set_sp )
-          CALL handle_netcdf_error( 'netcdf', 441 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 441 )
           IF ( averaging_interval_sp /= 0.0_wp )  THEN
              WRITE (time_average_text,'('', '',F7.1,'' s average'')') &
@@ -4833,5 +4545,5 @@
                                      TRIM( run_description_header ) // &
                                      TRIM( time_average_text ) )
-             CALL handle_netcdf_error( 'netcdf', 308 )
+             CALL netcdf_handle_error( 'netcdf_define_header', 308 )
 
              WRITE ( time_average_text,'(F7.1,'' s avg'')' )  averaging_interval_sp
@@ -4842,7 +4554,7 @@
                                      TRIM( run_description_header ) )
           ENDIF
-          CALL handle_netcdf_error( 'netcdf', 309 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 309 )
           nc_stat = NF90_ENDDEF( id_set_sp )
-          CALL handle_netcdf_error( 'netcdf', 442 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 442 )
           message_string = 'netCDF file for spectra ' //     &
                            'from previous run found.' //     &
@@ -4857,68 +4569,45 @@
           nc_stat = NF90_PUT_ATT( id_set_prt, NF90_GLOBAL, 'title', &
                                   TRIM( run_description_header ) )
-          CALL handle_netcdf_error( 'netcdf', 310 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 310 )
 
 !
 !--       Define time coordinate for particles (unlimited dimension)
-          nc_stat = NF90_DEF_DIM( id_set_prt, 'time', NF90_UNLIMITED, &
-                                  id_dim_time_prt )
-          CALL handle_netcdf_error( 'netcdf', 311 )
-
-          nc_stat = NF90_DEF_VAR( id_set_prt, 'time', NF90_DOUBLE, &
-                                  id_dim_time_prt, id_var_time_prt )
-          CALL handle_netcdf_error( 'netcdf', 312 )
-
-          nc_stat = NF90_PUT_ATT( id_set_prt, id_var_time_prt, 'units', &
-                                  'seconds' )
-          CALL handle_netcdf_error( 'netcdf', 313 )
-
+          CALL netcdf_create_dim( id_set_prt, 'time', NF90_UNLIMITED,          &
+                                  id_dim_time_prt, 311 )
+          CALL netcdf_create_var( id_set_prt, (/ id_dim_time_prt /), 'time',   &
+                                  NF90_DOUBLE, id_var_time_prt, 'seconds', '', &
+                                  312, 313, 000 )
 !
 !--       Define particle coordinate (maximum particle number)
           IF ( netcdf_data_format < 3 )  THEN
-             nc_stat = NF90_DEF_DIM( id_set_prt, 'particle_number', &
-                                     maximum_number_of_particles, id_dim_prtnum)
+             CALL netcdf_create_dim( id_set_prt, 'particle_number',            &
+                                     maximum_number_of_particles,              &
+                                     id_dim_prtnum, 314 )
           ELSE
 !
 !--          netCDF4 allows more than one unlimited dimension
-             nc_stat = NF90_DEF_DIM( id_set_prt, 'particle_number', &
-                                     NF90_UNLIMITED, id_dim_prtnum)
-          ENDIF
-          CALL handle_netcdf_error( 'netcdf', 314 )
-
-          nc_stat = NF90_DEF_VAR( id_set_prt, 'particle_number', NF90_DOUBLE, &
-                                  id_dim_prtnum, id_var_prtnum )
-          CALL handle_netcdf_error( 'netcdf', 315 )
-
-          nc_stat = NF90_PUT_ATT( id_set_prt, id_var_prtnum, 'units', &
-                                  'particle number' )
-          CALL handle_netcdf_error( 'netcdf', 316 )
-
+             CALL netcdf_create_dim( id_set_prt, 'particle_number',            &
+                                     NF90_UNLIMITED, id_dim_prtnum, 314 )
+          ENDIF
+
+          CALL netcdf_create_var( id_set_prt, (/ id_dim_prtnum /),             &
+                                  'particle_number', NF90_DOUBLE,              &
+                                  id_var_prtnum, 'particle number', '', 315,   &
+                                  316, 000 )
 !
 !--       Define variable which contains the real number of particles in use
-          nc_stat = NF90_DEF_VAR( id_set_prt, 'real_num_of_prt', NF90_INT, &
-                                  id_dim_time_prt, id_var_rnop_prt )
-          CALL handle_netcdf_error( 'netcdf', 317 )
-
-          nc_stat = NF90_PUT_ATT( id_set_prt, id_var_rnop_prt, 'units', &
-                                  'particle number' )
-          CALL handle_netcdf_error( 'netcdf', 318 )
-
+          CALL netcdf_create_var( id_set_prt, (/ id_dim_time_prt /),           &
+                                  'real_num_of_prt', NF90_DOUBLE,              &
+                                  id_var_rnop_prt, 'particle number', '', 317, &
+                                  318, 000 )
 !
 !--       Define the variables
           DO  i = 1, 17
 
-             nc_stat = NF90_DEF_VAR( id_set_prt, prt_var_names(i),         &
-                                     nc_precision(8),                      &
-                                     (/ id_dim_prtnum, id_dim_time_prt /), &
-                                     id_var_prt(i) )
-             CALL handle_netcdf_error( 'netcdf', 319 )
-
-             nc_stat = NF90_PUT_ATT( id_set_prt, id_var_prt(i), &
-                                     'long_name', TRIM( prt_var_names(i) ) )
-             CALL handle_netcdf_error( 'netcdf', 320 )
-
-             nc_stat = NF90_PUT_ATT( id_set_prt, id_var_prt(i), &
-                                     'units', TRIM( prt_var_units(i) ) )
-             CALL handle_netcdf_error( 'netcdf', 321 )
+             CALL netcdf_create_var( id_set_prt, (/ id_dim_prtnum,             &
+                                     id_dim_time_prt /), prt_var_names(i),     &
+                                     nc_precision(8), id_var_prt(i),           &
+                                     TRIM( prt_var_units(i) ),                 &
+                                     TRIM( prt_var_names(i) ), 319, 320, 321 )
 
           ENDDO
@@ -4927,5 +4616,5 @@
 !--       Leave netCDF define mode
           nc_stat = NF90_ENDDEF( id_set_prt )
-          CALL handle_netcdf_error( 'netcdf', 322 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 322 )
 
 
@@ -4938,14 +4627,14 @@
 !--       on the file.
           nc_stat = NF90_INQ_VARID( id_set_prt, 'time', id_var_time_prt )
-          CALL handle_netcdf_error( 'netcdf', 323 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 323 )
 
           nc_stat = NF90_INQUIRE_VARIABLE( id_set_prt, id_var_time_prt, &
                                            dimids = id_dim_time_old )
-          CALL handle_netcdf_error( 'netcdf', 324 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 324 )
           id_dim_time_prt = id_dim_time_old(1)
 
           nc_stat = NF90_INQUIRE_DIMENSION( id_set_prt, id_dim_time_prt, &
                                             len = prt_time_count )
-          CALL handle_netcdf_error( 'netcdf', 325 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 325 )
 
           nc_stat = NF90_GET_VAR( id_set_prt, id_var_time_prt,  &
@@ -4953,5 +4642,5 @@
                                   start = (/ prt_time_count /), &
                                   count = (/ 1 /) )
-          CALL handle_netcdf_error( 'netcdf', 326 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 326 )
 
           IF ( last_time_coordinate(1) >= simulated_time )  THEN
@@ -4974,5 +4663,5 @@
           nc_stat = NF90_INQ_VARID( id_set_prt, 'real_num_of_prt', &
                                     id_var_rnop_prt )
-          CALL handle_netcdf_error( 'netcdf', 327 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 327 )
 
           DO  i = 1, 17
@@ -4980,5 +4669,5 @@
              nc_stat = NF90_INQ_VARID( id_set_prt, prt_var_names(i), &
                                        id_var_prt(i) )
-             CALL handle_netcdf_error( 'netcdf', 328 )
+             CALL netcdf_handle_error( 'netcdf_define_header', 328 )
 
           ENDDO
@@ -4997,20 +4686,13 @@
           nc_stat = NF90_PUT_ATT( id_set_pts, NF90_GLOBAL, 'title', &
                                   TRIM( run_description_header ) )
-          CALL handle_netcdf_error( 'netcdf', 396 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 396 )
 
 !
 !--       Define time coordinate for particle time series (unlimited dimension)
-          nc_stat = NF90_DEF_DIM( id_set_pts, 'time', NF90_UNLIMITED, &
-                                  id_dim_time_pts )
-          CALL handle_netcdf_error( 'netcdf', 397 )
-
-          nc_stat = NF90_DEF_VAR( id_set_pts, 'time', NF90_DOUBLE, &
-                                  id_dim_time_pts, id_var_time_pts )
-          CALL handle_netcdf_error( 'netcdf', 398 )
-
-          nc_stat = NF90_PUT_ATT( id_set_pts, id_var_time_pts, 'units', &
-                                  'seconds')
-          CALL handle_netcdf_error( 'netcdf', 399 )
-
+          CALL netcdf_create_dim( id_set_pts, 'time', NF90_UNLIMITED,          &
+                                  id_dim_time_pts, 397 )
+          CALL netcdf_create_var( id_set_pts, (/ id_dim_time_pts /), 'time',   &
+                                  NF90_DOUBLE, id_var_time_pts, 'seconds', '', &
+                                  398, 399, 000 )
 !
 !--       Define the variables. If more than one particle group is defined,
@@ -5027,26 +4709,19 @@
                 ENDIF
 
-                nc_stat = NF90_DEF_VAR( id_set_pts,                         &
-                                        TRIM( dopts_label(i) ) // suffix1,  &
-                                        nc_precision(6),                    &
-                                        (/ id_dim_time_pts /),              &
-                                        id_var_dopts(i,j) )
-                CALL handle_netcdf_error( 'netcdf', 400 )
-
                 IF ( j == 0 )  THEN
-                   nc_stat = NF90_PUT_ATT( id_set_pts, id_var_dopts(i,j), &
-                                           'long_name', &
-                                           TRIM( dopts_label(i) ) )
+                   CALL netcdf_create_var( id_set_pts, (/ id_dim_time_pts /),  &
+                                           TRIM( dopts_label(i) ) // suffix1,  &
+                                           nc_precision(6), id_var_dopts(i,j), &
+                                           TRIM( dopts_unit(i) ),              &
+                                           TRIM( dopts_label(i) ), 400, 401,   &
+                                           402 )
                 ELSE
-                   nc_stat = NF90_PUT_ATT( id_set_pts, id_var_dopts(i,j), &
-                                           'long_name', &
+                   CALL netcdf_create_var( id_set_pts, (/ id_dim_time_pts /),  &
+                                           TRIM( dopts_label(i) ) // suffix1,  &
+                                           nc_precision(6), id_var_dopts(i,j), &
+                                           TRIM( dopts_unit(i) ),              &
                                            TRIM( dopts_label(i) ) // ' PG ' // &
-                                           suffix1(2:3) )
+                                           suffix1(2:3), 400, 401, 402 )
                 ENDIF
-                CALL handle_netcdf_error( 'netcdf', 401 )
-
-                nc_stat = NF90_PUT_ATT( id_set_pts, id_var_dopts(i,j), &
-                                        'units', TRIM( dopts_unit(i) ) )
-                CALL handle_netcdf_error( 'netcdf', 402 )
 
                 var_list = TRIM( var_list ) // TRIM( dopts_label(i) ) // &
@@ -5064,5 +4739,5 @@
           nc_stat = NF90_PUT_ATT( id_set_pts, NF90_GLOBAL, 'VAR_LIST', &
                                   var_list )
-          CALL handle_netcdf_error( 'netcdf', 403 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 403 )
 
 
@@ -5070,5 +4745,5 @@
 !--       Leave netCDF define mode
           nc_stat = NF90_ENDDEF( id_set_pts )
-          CALL handle_netcdf_error( 'netcdf', 404 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 404 )
 
 
@@ -5082,5 +4757,5 @@
           nc_stat = NF90_GET_ATT( id_set_pts, NF90_GLOBAL, 'VAR_LIST', &
                                   var_list_old )
-          CALL handle_netcdf_error( 'netcdf', 405 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 405 )
 
           var_list = ';'
@@ -5122,14 +4797,14 @@
 !--       on the file.
           nc_stat = NF90_INQ_VARID( id_set_pts, 'time', id_var_time_pts )
-          CALL handle_netcdf_error( 'netcdf', 406 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 406 )
 
           nc_stat = NF90_INQUIRE_VARIABLE( id_set_pts, id_var_time_pts, &
                                            dimids = id_dim_time_old )
-          CALL handle_netcdf_error( 'netcdf', 407 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 407 )
           id_dim_time_pts = id_dim_time_old(1)
 
           nc_stat = NF90_INQUIRE_DIMENSION( id_set_pts, id_dim_time_pts, &
                                             len = dopts_time_count )
-          CALL handle_netcdf_error( 'netcdf', 408 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 408 )
 
           nc_stat = NF90_GET_VAR( id_set_pts, id_var_time_pts,    &
@@ -5137,5 +4812,5 @@
                                   start = (/ dopts_time_count /), &
                                   count = (/ 1 /) )
-          CALL handle_netcdf_error( 'netcdf', 409 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 409 )
 
           IF ( last_time_coordinate(1) >= simulated_time )  THEN
@@ -5171,5 +4846,5 @@
                 nc_stat = NF90_INQ_VARID( id_set_pts, netcdf_var_name, &
                                           id_var_dopts(i,j) )
-                CALL handle_netcdf_error( 'netcdf', 410 )
+                CALL netcdf_handle_error( 'netcdf_define_header', 410 )
 
                 IF ( number_of_particle_groups == 1 )  EXIT
@@ -5187,14 +4862,14 @@
 !--       to ensure equal attribute size in a job chain. Maybe revise later.
           nc_stat = NF90_REDEF( id_set_pts )
-          CALL handle_netcdf_error( 'netcdf', 443 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 443 )
           nc_stat = NF90_PUT_ATT( id_set_pts, NF90_GLOBAL, 'title', &
                                   TRIM( run_description_header ) )
-          CALL handle_netcdf_error( 'netcdf', 411 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 411 )
           nc_stat = NF90_ENDDEF( id_set_pts )
-          CALL handle_netcdf_error( 'netcdf', 444 )
+          CALL netcdf_handle_error( 'netcdf_define_header', 444 )
           message_string = 'netCDF file for particle time series ' // &
                            'from previous run found.' //              &
                            '&This file will be extended.'
-          CALL message( 'define_netcdf_header', 'PA0269', 0, 0, 0, 6, 0 )
+          CALL message( 'netcdf_define_header', 'PA0269', 0, 0, 0, 6, 0 )
 
 
@@ -5202,10 +4877,10 @@
 
           message_string = 'mode "' // TRIM( mode) // '" not supported'
-          CALL message( 'define_netcdf_header', 'PA0270', 0, 0, 0, 6, 0 )
+          CALL message( 'netcdf_define_header', 'PA0270', 0, 0, 0, 6, 0 )
 
     END SELECT
 
 #endif
- END SUBROUTINE define_netcdf_header
+ END SUBROUTINE netcdf_define_header
 
 
@@ -5217,11 +4892,7 @@
 !------------------------------------------------------------------------------!
  
- SUBROUTINE create_netcdf_file( filename , id, parallel, errno )
+ SUBROUTINE netcdf_create_file( filename , id, parallel, errno )
 #if defined( __netcdf )
 
-
-    USE control_parameters
-    USE netcdf
-    USE netcdf_control
     USE pegrid
 
@@ -5284,7 +4955,7 @@
     ENDIF
 
-    CALL handle_netcdf_error( 'create_netcdf_file', errno )
+    CALL netcdf_handle_error( 'netcdf_create_file', errno )
 #endif
- END SUBROUTINE create_netcdf_file
+ END SUBROUTINE netcdf_create_file
 
 
@@ -5296,14 +4967,6 @@
 !------------------------------------------------------------------------------!
  
- SUBROUTINE open_write_netcdf_file( filename, id, parallel, errno )
+ SUBROUTINE netcdf_open_write_file( filename, id, parallel, errno )
 #if defined( __netcdf )
-
-
-    USE control_parameters,                                                    &
-        ONLY:  netcdf_data_format
-
-    USE netcdf
-
-    USE netcdf_control
 
     USE pegrid
@@ -5328,7 +4991,7 @@
     ENDIF
 
-    CALL handle_netcdf_error( 'open_write_netcdf_file', errno )
+    CALL netcdf_handle_error( 'netcdf_open_write_file', errno )
 #endif
- END SUBROUTINE open_write_netcdf_file
+ END SUBROUTINE netcdf_open_write_file
 
 
@@ -5339,5 +5002,5 @@
 !------------------------------------------------------------------------------!
  
- SUBROUTINE handle_netcdf_error( routine_name, errno )
+ SUBROUTINE netcdf_handle_error( routine_name, errno )
 #if defined( __netcdf )
 
@@ -5346,10 +5009,4 @@
         ONLY:  message_string
 
-    USE netcdf
-
-    USE netcdf_control
-
-    USE pegrid
-
     IMPLICIT NONE
 
@@ -5369,3 +5026,94 @@
 
 #endif
- END SUBROUTINE handle_netcdf_error
+ END SUBROUTINE netcdf_handle_error
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Create a dimension in NetCDF file
+!------------------------------------------------------------------------------!
+
+ SUBROUTINE netcdf_create_dim(ncid, dim_name, ncdim_type, ncdim_id, error_no)
+
+#if defined( __netcdf )
+
+    USE kinds
+
+    IMPLICIT NONE
+
+    CHARACTER(LEN=*), INTENT(IN) ::  dim_name
+
+    INTEGER, INTENT(IN)  ::  error_no
+    INTEGER, INTENT(IN)  ::  ncid
+    INTEGER, INTENT(OUT) ::  ncdim_id
+    INTEGER, INTENT(IN)  ::  ncdim_type
+
+!
+!-- Define time coordinate for volume data (unlimited dimension)
+    nc_stat = NF90_DEF_DIM( ncid, dim_name, ncdim_type, ncdim_id )
+    CALL netcdf_handle_error( 'netcdf_create_dim', error_no )
+
+#endif
+
+ END SUBROUTINE netcdf_create_dim
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Create a one dimensional variable in specific units in NetCDF file
+!------------------------------------------------------------------------------!
+
+ SUBROUTINE netcdf_create_var( ncid, dim_id, var_name, var_type, var_id,       &
+                               unit_name, long_name, error_no1, error_no2,     &
+                               error_no3 )
+
+#if defined( __netcdf )
+    IMPLICIT NONE
+
+    CHARACTER(LEN=*), INTENT(IN) ::  long_name
+    CHARACTER(LEN=*), INTENT(IN) ::  unit_name
+    CHARACTER(LEN=*), INTENT(IN) ::  var_name
+
+    INTEGER, INTENT(IN)  ::  error_no1
+    INTEGER, INTENT(IN)  ::  error_no2
+    INTEGER, INTENT(IN)  ::  error_no3
+    INTEGER, INTENT(IN)  ::  ncid
+    INTEGER, INTENT(OUT) ::  var_id
+    INTEGER, INTENT(IN)  ::  var_type
+
+    INTEGER, DIMENSION(:), INTENT(IN) ::  dim_id
+
+!
+!-- Define variable
+    nc_stat = NF90_DEF_VAR( ncid, var_name, var_type, dim_id, var_id )
+    CALL netcdf_handle_error( 'netcdf_create_var', error_no1 )
+
+#if defined( __netcdf4 )
+!
+!-- Check if variable should be deflate (including shuffling)
+!-- and if it is possible (only NetCDF4 with HDF5 supports compression)
+    IF ( netcdf_data_format > 2  .AND.  netcdf_deflate > 0 )  THEN
+       nc_stat = NF90_DEF_VAR_DEFLATE( ncid, var_id, 1, 1, netcdf_deflate )
+       CALL netcdf_handle_error( 'netcdf_create_var_deflate', error_no1 )
+    ENDIF
+#endif
+!
+!-- Set unit name if set
+    IF ( unit_name /= '' )  THEN
+       nc_stat = NF90_PUT_ATT( ncid, var_id, 'units', unit_name )
+       CALL netcdf_handle_error( 'netcdf_create_var', error_no2 )
+    ENDIF
+
+!
+!-- Set long name if set
+    IF ( long_name /= '' )  THEN
+       nc_stat = NF90_PUT_ATT( ncid, var_id, 'long_name', long_name )
+       CALL netcdf_handle_error( 'netcdf_create_var', error_no3 )
+    ENDIF
+
+#endif
+ END SUBROUTINE netcdf_create_var
+
+ END MODULE netcdf_interface
