Index: palm/trunk/SOURCE/netcdf_interface_mod.f90
===================================================================
--- palm/trunk/SOURCE/netcdf_interface_mod.f90	(revision 3337)
+++ palm/trunk/SOURCE/netcdf_interface_mod.f90	(revision 3421)
@@ -25,4 +25,9 @@
 ! -----------------
 ! $Id$
+! Bugfix: move ocean output variables to ocean_mod
+! Renamed output variables
+! Add UTM coordinates to mask, 3d, xy, xz, yz output
+! 
+! 3337 2018-10-12 15:17:09Z kanani
 ! (from branch resler)
 ! Add biometeorology
@@ -330,11 +335,11 @@
     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(dopr_norm_num) ::  dopr_norm_names =          &
+         (/ 'wtheta0', 'ws2    ', 'tsw2   ', 'ws3    ', 'ws2tsw ', 'wstsw2 ',  &
+            'z_i    ' /)
+
+    CHARACTER (LEN=7), DIMENSION(dopr_norm_num) ::  dopr_norm_longnames =      &
+         (/ 'wtheta0', 'w*2    ', 't*w2   ', 'w*3    ', 'w*2t*w ', 'w*t*w2 ',  &
+            'z_i    ' /)
 
     CHARACTER (LEN=7), DIMENSION(dopts_num) :: dopts_label =                   &
@@ -358,9 +363,9 @@
     CHARACTER (LEN=13), DIMENSION(dots_max) :: dots_label =                    &
           (/ 'E            ', 'E*           ', 'dt           ',                &
-             'u*           ', 'th*          ', 'umax         ',                &
+             'us*          ', '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)     ',                &
+             'div_old      ', 'zi_wtheta    ', 'zi_theta     ',                &
+             'w*           ', 'w"theta"0    ', 'w"theta"     ',                &
+             'wtheta       ', 'theta(0)     ', 'theta(z_mo)  ',                &
              'w"u"0        ', 'w"v"0        ', 'w"q"0        ',                &
              'ol           ', 'q*           ', 'w"s"         ',                &
@@ -456,4 +461,10 @@
                     id_var_zw_xy, id_var_zw_xz, id_var_zw_yz, id_var_zw_3d
 
+    INTEGER(iwp), DIMENSION(0:1,0:2) ::  id_var_eutm_3d, id_var_nutm_3d, &
+                                         id_var_eutm_xy, id_var_nutm_xy, &
+                                         id_var_eutm_xz, id_var_nutm_xz, &
+                                         id_var_eutm_yz, id_var_nutm_yz, &
+                                         id_var_eutm_ma, id_var_nutm_ma
+
     INTEGER ::  netcdf_data_format = 2  !< NetCDF3 64bit offset format
     INTEGER ::  netcdf_deflate = 0      !< NetCDF compression, default: no
@@ -498,4 +509,7 @@
                    id_var_zu_mask, id_var_zw_mask, &
                    id_var_zusi_mask, id_var_zwwi_mask
+                   
+    INTEGER(iwp), DIMENSION(1:max_masks,0:1,0:2) ::  id_var_eutm_mask, &
+                                                     id_var_nutm_mask
 
     INTEGER(iwp), DIMENSION(1:max_masks,0:1,100) ::  id_var_domask
@@ -596,4 +610,7 @@
     USE land_surface_model_mod,                                                &
         ONLY: lsm_define_netcdf_grid, nzb_soil, nzt_soil, nzs, zs
+
+    USE netcdf_data_input_mod,                                                 &
+        ONLY: init_model
 
     USE ocean_mod,                                                             &
@@ -692,4 +709,9 @@
 
     LOGICAL, SAVE ::  init_netcdf = .FALSE.                  !<
+
+    REAL(wp) ::  cos_ra                                      !< cosine of rotation_angle
+    REAL(wp) ::  shift_x                                     !< shift of x coordinate
+    REAL(wp) ::  shift_y                                     !< shift of y coordinate
+    REAL(wp) ::  sin_ra                                      !< sine of rotation_angle
 
     REAL(wp), DIMENSION(1) ::  last_time_coordinate          !< last time value in file
@@ -901,4 +923,57 @@
                                   'meters', '', 486, 487, 000 )
 !
+!--       Define UTM coordinates
+          IF ( init_model%rotation_angle == 0.0_wp )  THEN
+             CALL netcdf_create_var( id_set_mask(mid,av), &
+                                 (/ id_dim_x_mask(mid,av) /),      &
+                                 'E_UTM', NF90_DOUBLE, id_var_eutm_mask(mid,av,0),  &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_mask(mid,av), &
+                                 (/ id_dim_y_mask(mid,av) /),      &
+                                 'N_UTM', NF90_DOUBLE, id_var_nutm_mask(mid,av,0),  &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_mask(mid,av), &
+                                 (/ id_dim_xu_mask(mid,av) /),     &
+                                 'Eu_UTM', NF90_DOUBLE, id_var_eutm_mask(mid,av,1), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_mask(mid,av), &
+                                 (/ id_dim_y_mask(mid,av) /),     &
+                                 'Nu_UTM', NF90_DOUBLE, id_var_nutm_mask(mid,av,1), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_mask(mid,av), &
+                                 (/ id_dim_x_mask(mid,av) /),     &
+                                 'Ev_UTM', NF90_DOUBLE, id_var_eutm_mask(mid,av,2), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_mask(mid,av), &
+                                 (/ id_dim_yv_mask(mid,av) /),     &
+                                 'Nv_UTM', NF90_DOUBLE, id_var_nutm_mask(mid,av,2), &
+                                 'm', '', 000, 000, 000 )
+          ELSE
+             CALL netcdf_create_var( id_set_mask(mid,av), &
+                                 (/ id_dim_x_mask(mid,av), id_dim_y_mask(mid,av) /),      &
+                                 'E_UTM', NF90_DOUBLE, id_var_eutm_mask(mid,av,0),  &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_mask(mid,av), &
+                                 (/ id_dim_x_mask(mid,av), id_dim_y_mask(mid,av) /),      &
+                                 'N_UTM', NF90_DOUBLE, id_var_nutm_mask(mid,av,0),  &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_mask(mid,av), &
+                                 (/ id_dim_xu_mask(mid,av), id_dim_y_mask(mid,av) /),     &
+                                 'Eu_UTM', NF90_DOUBLE, id_var_eutm_mask(mid,av,1), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_mask(mid,av), &
+                                 (/ id_dim_xu_mask(mid,av), id_dim_y_mask(mid,av) /),     &
+                                 'Nu_UTM', NF90_DOUBLE, id_var_nutm_mask(mid,av,1), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_mask(mid,av), &
+                                 (/ id_dim_x_mask(mid,av), id_dim_yv_mask(mid,av) /),     &
+                                 'Ev_UTM', NF90_DOUBLE, id_var_eutm_mask(mid,av,2), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_mask(mid,av), &
+                                 (/ id_dim_x_mask(mid,av), id_dim_yv_mask(mid,av) /),     &
+                                 'Nv_UTM', NF90_DOUBLE, id_var_nutm_mask(mid,av,2), &
+                                 'm', '', 000, 000, 000 )
+          ENDIF
+!
 !--       In case of non-flat topography define 2d-arrays containing the height
 !--       information. Only for parallel netcdf output.
@@ -954,7 +1029,7 @@
 !
 !--             Most variables are defined on the scalar grid
-                CASE ( 'e', 'lpt', 'nc', 'nr', 'p', 'pc', 'pr', 'prr', 'pt',   &
+                CASE ( 'e', 'nc', 'nr', 'p', 'pc', 'pr', 'prr',                &
                        'q', 'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',   &
-                       's', 'vpt'  )
+                       's', 'theta', 'thetal', 'thetav' )
 
                    grid_x = 'x'
@@ -1144,4 +1219,125 @@
 
 !
+!--       Write UTM coordinates
+          IF ( init_model%rotation_angle == 0.0_wp )  THEN
+!
+!--          1D in case of no rotation
+             cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
+!
+!--          x coordinates
+             ALLOCATE( netcdf_data(mask_size(mid,1)) )
+             DO  k = 0, 2
+!            
+!--             Scalar grid points
+                IF ( k == 0 )  THEN
+                   shift_x = 0.5
+!            
+!--             u grid points
+                ELSEIF ( k == 1 )  THEN
+                   shift_x = 0.0
+!            
+!--             v grid points
+                ELSEIF ( k == 2 )  THEN
+                   shift_x = 0.5
+                ENDIF
+
+                netcdf_data = init_model%origin_x + cos_ra                     &
+                       * ( mask_i_global(mid,:mask_size(mid,1)) + shift_x ) * dx
+
+                nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), &
+                                        id_var_eutm_mask(mid,av,k), &
+                                        netcdf_data, start = (/ 1 /), &
+                                        count = (/ mask_size(mid,1) /) )
+                CALL netcdf_handle_error( 'netcdf_define_header', 555 )
+
+             ENDDO
+             DEALLOCATE( netcdf_data )
+!
+!--          y coordinates
+             ALLOCATE( netcdf_data(mask_size(mid,2)) )
+             DO  k = 0, 2
+!
+!--             Scalar grid points
+                IF ( k == 0 )  THEN
+                   shift_y = 0.5
+!
+!--             u grid points
+                ELSEIF ( k == 1 )  THEN
+                   shift_y = 0.5
+!
+!--             v grid points
+                ELSEIF ( k == 2 )  THEN
+                   shift_y = 0.0
+                ENDIF
+
+                netcdf_data = init_model%origin_y + cos_ra                     &
+                       * ( mask_j_global(mid,:mask_size(mid,2)) + shift_y ) * dy
+
+                nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), &
+                                        id_var_nutm_mask(mid,av,k), &
+                                        netcdf_data, start = (/ 1 /), &
+                                        count = (/ mask_size(mid,2) /) )
+                CALL netcdf_handle_error( 'netcdf_define_header', 556 )
+
+             ENDDO
+             DEALLOCATE( netcdf_data )
+
+          ELSE
+!
+!--          2D in case of rotation
+             ALLOCATE( netcdf_data_2d(mask_size(mid,1),mask_size(mid,2)) )
+             cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
+             sin_ra = SIN( init_model%rotation_angle * pi / 180.0_wp )
+
+             DO  k = 0, 2
+!            
+!--            Scalar grid points
+               IF ( k == 0 )  THEN
+                  shift_x = 0.5 ; shift_y = 0.5
+!            
+!--            u grid points
+               ELSEIF ( k == 1 )  THEN
+                  shift_x = 0.0 ; shift_y = 0.5
+!            
+!--            v grid points
+               ELSEIF ( k == 2 )  THEN
+                  shift_x = 0.5 ; shift_y = 0.0
+               ENDIF
+
+               DO  j = 1, mask_size(mid,2)
+                  DO  i = 1, mask_size(mid,1)
+                     netcdf_data_2d(i,j) = init_model%origin_x                 &
+                           + cos_ra * ( mask_i_global(mid,i) + shift_x ) * dx  &
+                           + sin_ra * ( mask_j_global(mid,j) + shift_y ) * dy
+                  ENDDO
+               ENDDO
+
+               nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), &
+                                       id_var_eutm_mask(mid,av,k), &
+                                       netcdf_data_2d, start = (/ 1, 1 /), &
+                                       count = (/ mask_size(mid,1), &
+                                                  mask_size(mid,2) /) )
+               CALL netcdf_handle_error( 'netcdf_define_header', 555 )
+
+               DO  j = 1, mask_size(mid,2)
+                  DO  i = 1, mask_size(mid,1)
+                     netcdf_data_2d(i,j) = init_model%origin_y                 &
+                           - sin_ra * ( mask_i_global(mid,i) + shift_x ) * dx  &
+                           + cos_ra * ( mask_j_global(mid,j) + shift_y ) * dy
+                  ENDDO
+               ENDDO
+             
+               nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), &
+                                       id_var_nutm_mask(mid,av,k), &
+                                       netcdf_data_2d, start = (/ 1, 1 /), &
+                                       count = (/ mask_size(mid,1), &
+                                                  mask_size(mid,2) /) )
+               CALL netcdf_handle_error( 'netcdf_define_header', 556 )
+             
+             ENDDO
+             DEALLOCATE( netcdf_data_2d )
+          ENDIF
+
+!
 !--       Write zu and zw data (vertical axes)
           ALLOCATE( netcdf_data(mask_size(mid,3)) )
@@ -1456,4 +1652,57 @@
                                   362, 363, 000 )
 !
+!--       Define UTM coordinates
+          IF ( init_model%rotation_angle == 0.0_wp )  THEN
+             CALL netcdf_create_var( id_set_3d(av), &
+                                 (/ id_dim_x_3d(av) /),      &
+                                 'E_UTM', NF90_DOUBLE, id_var_eutm_3d(av,0),  &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_3d(av), &
+                                 (/ id_dim_y_3d(av) /),      &
+                                 'N_UTM', NF90_DOUBLE, id_var_nutm_3d(av,0),  &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_3d(av), &
+                                 (/ id_dim_xu_3d(av) /),     &
+                                 'Eu_UTM', NF90_DOUBLE, id_var_eutm_3d(av,1), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_3d(av), &
+                                 (/ id_dim_y_3d(av) /),     &
+                                 'Nu_UTM', NF90_DOUBLE, id_var_nutm_3d(av,1), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_3d(av), &
+                                 (/ id_dim_x_3d(av) /),     &
+                                 'Ev_UTM', NF90_DOUBLE, id_var_eutm_3d(av,2), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_3d(av), &
+                                 (/ id_dim_yv_3d(av) /),     &
+                                 'Nv_UTM', NF90_DOUBLE, id_var_nutm_3d(av,2), &
+                                 'm', '', 000, 000, 000 )
+          ELSE
+             CALL netcdf_create_var( id_set_3d(av), &
+                                 (/ id_dim_x_3d(av), id_dim_y_3d(av) /),      &
+                                 'E_UTM', NF90_DOUBLE, id_var_eutm_3d(av,0),  &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_3d(av), &
+                                 (/ id_dim_x_3d(av), id_dim_y_3d(av) /),      &
+                                 'N_UTM', NF90_DOUBLE, id_var_nutm_3d(av,0),  &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_3d(av), &
+                                 (/ id_dim_xu_3d(av), id_dim_y_3d(av) /),     &
+                                 'Eu_UTM', NF90_DOUBLE, id_var_eutm_3d(av,1), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_3d(av), &
+                                 (/ id_dim_xu_3d(av), id_dim_y_3d(av) /),     &
+                                 'Nu_UTM', NF90_DOUBLE, id_var_nutm_3d(av,1), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_3d(av), &
+                                 (/ id_dim_x_3d(av), id_dim_yv_3d(av) /),     &
+                                 'Ev_UTM', NF90_DOUBLE, id_var_eutm_3d(av,2), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_3d(av), &
+                                 (/ id_dim_x_3d(av), id_dim_yv_3d(av) /),     &
+                                 'Nv_UTM', NF90_DOUBLE, id_var_nutm_3d(av,2), &
+                                 'm', '', 000, 000, 000 )
+          ENDIF
+!
 !--       In case of non-flat topography define 2d-arrays containing the height
 !--       information. Only output 2d topography information in case of parallel
@@ -1506,7 +1755,7 @@
 !
 !--             Most variables are defined on the scalar grid
-                CASE ( 'e', 'lpt', 'nc', 'nr', 'p', 'pc', 'pr', 'prr', 'pt',   &
+                CASE ( 'e', 'nc', 'nr', 'p', 'pc', 'pr', 'prr',   &
                        'q', 'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',   &
-                       'rho_ocean', 's', 'sa', 'vpt' )
+                       's', 'theta', 'thetal', 'thetav' )
 
                    grid_x = 'x'
@@ -1551,4 +1800,10 @@
                       CALL lsm_define_netcdf_grid( do3d(av,i), found, grid_x,  &
                                                    grid_y, grid_z )
+                   ENDIF
+!
+!--                Check for ocean quantities
+                   IF ( .NOT. found  .AND.  ocean_mode )  THEN
+                      CALL ocean_define_netcdf_grid( do3d(av,i), found,  &
+                                                     grid_x, grid_y, grid_z )
                    ENDIF
 
@@ -1729,4 +1984,122 @@
              DEALLOCATE( netcdf_data )
 
+!
+!--          Write UTM coordinates
+             IF ( init_model%rotation_angle == 0.0_wp )  THEN
+!
+!--             1D in case of no rotation
+                cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
+!
+!--             x coordinates
+                ALLOCATE( netcdf_data(0:nx) )
+                DO  k = 0, 2
+!               
+!--                Scalar grid points
+                   IF ( k == 0 )  THEN
+                      shift_x = 0.5
+!               
+!--                u grid points
+                   ELSEIF ( k == 1 )  THEN
+                      shift_x = 0.0
+!               
+!--                v grid points
+                   ELSEIF ( k == 2 )  THEN
+                      shift_x = 0.5
+                   ENDIF
+                
+                   DO  i = 0, nx
+                     netcdf_data(i) = init_model%origin_x            &
+                                    + cos_ra * ( i + shift_x ) * dx
+                   ENDDO
+                
+                   nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_eutm_3d(av,k),&
+                                           netcdf_data, start = (/ 1 /),   &
+                                           count = (/ nx+1 /) )
+                   CALL netcdf_handle_error( 'netcdf_define_header', 555 )
+
+                ENDDO
+                DEALLOCATE( netcdf_data )
+!
+!--             y coordinates
+                ALLOCATE( netcdf_data(0:ny) )
+                DO  k = 0, 2
+!
+!--                Scalar grid points
+                   IF ( k == 0 )  THEN
+                      shift_y = 0.5
+!
+!--                u grid points
+                   ELSEIF ( k == 1 )  THEN
+                      shift_y = 0.5
+!
+!--                v grid points
+                   ELSEIF ( k == 2 )  THEN
+                      shift_y = 0.0
+                   ENDIF
+
+                   DO  j = 0, ny
+                      netcdf_data(j) = init_model%origin_y            &
+                                     + cos_ra * ( j + shift_y ) * dy
+                   ENDDO
+
+                   nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_nutm_3d(av,k),&
+                                           netcdf_data, start = (/ 1 /),   &
+                                           count = (/ ny+1 /) )
+                   CALL netcdf_handle_error( 'netcdf_define_header', 556 )
+
+                ENDDO
+                DEALLOCATE( netcdf_data )
+
+             ELSE
+!
+!--             2D in case of rotation
+                ALLOCATE( netcdf_data_2d(0:nx,0:ny) )
+                cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
+                sin_ra = SIN( init_model%rotation_angle * pi / 180.0_wp )
+                
+                DO  k = 0, 2
+!               
+!--               Scalar grid points
+                  IF ( k == 0 )  THEN
+                     shift_x = 0.5 ; shift_y = 0.5
+!               
+!--               u grid points
+                  ELSEIF ( k == 1 )  THEN
+                     shift_x = 0.0 ; shift_y = 0.5
+!               
+!--               v grid points
+                  ELSEIF ( k == 2 )  THEN
+                     shift_x = 0.5 ; shift_y = 0.0
+                  ENDIF
+                
+                  DO  j = 0, ny
+                     DO  i = 0, nx
+                        netcdf_data_2d(i,j) = init_model%origin_x            &
+                                            + cos_ra * ( i + shift_x ) * dx  &
+                                            + sin_ra * ( j + shift_y ) * dy
+                     ENDDO
+                  ENDDO
+                
+                  nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_eutm_3d(av,k),  &
+                                          netcdf_data_2d, start = (/ 1, 1 /),   &
+                                          count = (/ nx+1, ny+1 /) )
+                  CALL netcdf_handle_error( 'netcdf_define_header', 555 )
+                
+                  DO  j = 0, ny
+                     DO  i = 0, nx
+                        netcdf_data_2d(i,j) = init_model%origin_y            &
+                                            - sin_ra * ( i + shift_x ) * dx  &
+                                            + cos_ra * ( j + shift_y ) * dy
+                     ENDDO
+                  ENDDO
+                
+                  nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_nutm_3d(av,k),  &
+                                          netcdf_data_2d, start = (/ 1, 1 /),   &
+                                          count = (/ nx+1, ny+1 /) )
+                  CALL netcdf_handle_error( 'netcdf_define_header', 556 )
+                
+                ENDDO
+                DEALLOCATE( netcdf_data_2d )
+             ENDIF
 !
 !--          Write zu and zw data (vertical axes)
@@ -2242,4 +2615,57 @@
                                   365, 366, 000 )
 !
+!--       Define UTM coordinates
+          IF ( init_model%rotation_angle == 0.0_wp )  THEN
+             CALL netcdf_create_var( id_set_xy(av), &
+                                 (/ id_dim_x_xy(av) /),      &
+                                 'E_UTM', NF90_DOUBLE, id_var_eutm_xy(av,0),  &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_xy(av), &
+                                 (/ id_dim_y_xy(av) /),      &
+                                 'N_UTM', NF90_DOUBLE, id_var_nutm_xy(av,0),  &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_xy(av), &
+                                 (/ id_dim_xu_xy(av) /),     &
+                                 'Eu_UTM', NF90_DOUBLE, id_var_eutm_xy(av,1), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_xy(av), &
+                                 (/ id_dim_y_xy(av) /),     &
+                                 'Nu_UTM', NF90_DOUBLE, id_var_nutm_xy(av,1), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_xy(av), &
+                                 (/ id_dim_x_xy(av) /),     &
+                                 'Ev_UTM', NF90_DOUBLE, id_var_eutm_xy(av,2), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_xy(av), &
+                                 (/ id_dim_yv_xy(av) /),     &
+                                 'Nv_UTM', NF90_DOUBLE, id_var_nutm_xy(av,2), &
+                                 'm', '', 000, 000, 000 )
+          ELSE
+             CALL netcdf_create_var( id_set_xy(av), &
+                                 (/ id_dim_x_xy(av), id_dim_y_xy(av) /),      &
+                                 'E_UTM', NF90_DOUBLE, id_var_eutm_xy(av,0),  &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_xy(av), &
+                                 (/ id_dim_x_xy(av), id_dim_y_xy(av) /),      &
+                                 'N_UTM', NF90_DOUBLE, id_var_nutm_xy(av,0),  &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_xy(av), &
+                                 (/ id_dim_xu_xy(av), id_dim_y_xy(av) /),     &
+                                 'Eu_UTM', NF90_DOUBLE, id_var_eutm_xy(av,1), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_xy(av), &
+                                 (/ id_dim_xu_xy(av), id_dim_y_xy(av) /),     &
+                                 'Nu_UTM', NF90_DOUBLE, id_var_nutm_xy(av,1), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_xy(av), &
+                                 (/ id_dim_x_xy(av), id_dim_yv_xy(av) /),     &
+                                 'Ev_UTM', NF90_DOUBLE, id_var_eutm_xy(av,2), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_xy(av), &
+                                 (/ id_dim_x_xy(av), id_dim_yv_xy(av) /),     &
+                                 'Nv_UTM', NF90_DOUBLE, id_var_nutm_xy(av,2), &
+                                 'm', '', 000, 000, 000 )
+          ENDIF
+!
 !--       In case of non-flat topography define 2d-arrays containing the height
 !--       information. Only for parallel netcdf output.
@@ -2289,9 +2715,10 @@
 !
 !--                   Most variables are defined on the zu grid
-                      CASE ( 'e_xy', 'lpt_xy', 'nc_xy', 'nr_xy', 'p_xy',       &
-                             'pc_xy', 'pr_xy', 'prr_xy', 'pt_xy', 'q_xy',      &
+                      CASE ( 'e_xy', 'nc_xy', 'nr_xy', 'p_xy',                 &
+                             'pc_xy', 'pr_xy', 'prr_xy', 'q_xy',               &
                              'qc_xy', 'ql_xy', 'ql_c_xy', 'ql_v_xy',           &
-                             'ql_vp_xy', 'qr_xy', 'qv_xy', 'rho_ocean_xy',     &
-                             's_xy', 'sa_xy', 'vpt_xy' )
+                             'ql_vp_xy', 'qr_xy', 'qv_xy',                     &
+                             's_xy',                                           &
+                             'theta_xy', 'thetal_xy', 'thetav_xy' )
 
                          grid_x = 'x'
@@ -2335,4 +2762,11 @@
                          ENDIF
 
+!
+!--                      Check for ocean quantities
+                         IF ( .NOT. found  .AND.  ocean_mode )  THEN
+                            CALL ocean_define_netcdf_grid( do2d(av,i), found,  &
+                                                           grid_x, grid_y,     &
+                                                           grid_z )
+                         ENDIF
 !
 !--                      Check for radiation quantities
@@ -2591,7 +3025,124 @@
 
              DEALLOCATE( netcdf_data )
-
-          ENDIF
-
+!
+!--          Write UTM coordinates
+             IF ( init_model%rotation_angle == 0.0_wp )  THEN
+!
+!--             1D in case of no rotation
+                cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
+!
+!--             x coordinates
+                ALLOCATE( netcdf_data(0:nx) )
+                DO  k = 0, 2
+!               
+!--                Scalar grid points
+                   IF ( k == 0 )  THEN
+                      shift_x = 0.5
+!               
+!--                u grid points
+                   ELSEIF ( k == 1 )  THEN
+                      shift_x = 0.0
+!               
+!--                v grid points
+                   ELSEIF ( k == 2 )  THEN
+                      shift_x = 0.5
+                   ENDIF
+                
+                   DO  i = 0, nx
+                     netcdf_data(i) = init_model%origin_x            &
+                                    + cos_ra * ( i + shift_x ) * dx
+                   ENDDO
+                
+                   nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_eutm_xy(av,k),&
+                                           netcdf_data, start = (/ 1 /),   &
+                                           count = (/ nx+1 /) )
+                   CALL netcdf_handle_error( 'netcdf_define_header', 555 )
+
+                ENDDO
+                DEALLOCATE( netcdf_data )
+!
+!--             y coordinates
+                ALLOCATE( netcdf_data(0:ny) )
+                DO  k = 0, 2
+!
+!--                Scalar grid points
+                   IF ( k == 0 )  THEN
+                      shift_y = 0.5
+!
+!--                u grid points
+                   ELSEIF ( k == 1 )  THEN
+                      shift_y = 0.5
+!
+!--                v grid points
+                   ELSEIF ( k == 2 )  THEN
+                      shift_y = 0.0
+                   ENDIF
+
+                   DO  j = 0, ny
+                      netcdf_data(j) = init_model%origin_y            &
+                                     + cos_ra * ( j + shift_y ) * dy
+                   ENDDO
+
+                   nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_nutm_xy(av,k),&
+                                           netcdf_data, start = (/ 1 /),   &
+                                           count = (/ ny+1 /) )
+                   CALL netcdf_handle_error( 'netcdf_define_header', 556 )
+
+                ENDDO
+                DEALLOCATE( netcdf_data )
+
+             ELSE
+!
+!--             2D in case of rotation
+                ALLOCATE( netcdf_data_2d(0:nx,0:ny) )
+                cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
+                sin_ra = SIN( init_model%rotation_angle * pi / 180.0_wp )
+                
+                DO  k = 0, 2
+!               
+!--               Scalar grid points
+                  IF ( k == 0 )  THEN
+                     shift_x = 0.5 ; shift_y = 0.5
+!               
+!--               u grid points
+                  ELSEIF ( k == 1 )  THEN
+                     shift_x = 0.0 ; shift_y = 0.5
+!               
+!--               v grid points
+                  ELSEIF ( k == 2 )  THEN
+                     shift_x = 0.5 ; shift_y = 0.0
+                  ENDIF
+                
+                  DO  j = 0, ny
+                     DO  i = 0, nx
+                        netcdf_data_2d(i,j) = init_model%origin_x            &
+                                            + cos_ra * ( i + shift_x ) * dx  &
+                                            + sin_ra * ( j + shift_y ) * dy
+                     ENDDO
+                  ENDDO
+                
+                  nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_eutm_xy(av,k),  &
+                                          netcdf_data_2d, start = (/ 1, 1 /),   &
+                                          count = (/ nx+1, ny+1 /) )
+                  CALL netcdf_handle_error( 'netcdf_define_header', 555 )
+                
+                  DO  j = 0, ny
+                     DO  i = 0, nx
+                        netcdf_data_2d(i,j) = init_model%origin_y            &
+                                            - sin_ra * ( i + shift_x ) * dx  &
+                                            + cos_ra * ( j + shift_y ) * dy
+                     ENDDO
+                  ENDDO
+                
+                  nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_nutm_xy(av,k),  &
+                                          netcdf_data_2d, start = (/ 1, 1 /),   &
+                                          count = (/ nx+1, ny+1 /) )
+                  CALL netcdf_handle_error( 'netcdf_define_header', 556 )
+                
+                ENDDO
+                DEALLOCATE( netcdf_data_2d )
+             ENDIF
+
+          ENDIF
 !
 !--       In case of non-flat topography write height information. Only for 
@@ -2985,4 +3536,57 @@
                                   NF90_DOUBLE, id_var_zw_xz(av), 'meters', '', &
                                   157, 158, 000 )
+!
+!--       Define UTM coordinates
+          IF ( init_model%rotation_angle == 0.0_wp )  THEN
+             CALL netcdf_create_var( id_set_xz(av), &
+                                 (/ id_dim_x_xz(av) /),      &
+                                 'E_UTM', NF90_DOUBLE, id_var_eutm_xz(av,0),  &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_xz(av), &
+                                 (/ id_dim_y_xz(av) /),      &
+                                 'N_UTM', NF90_DOUBLE, id_var_nutm_xz(av,0),  &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_xz(av), &
+                                 (/ id_dim_xu_xz(av) /),     &
+                                 'Eu_UTM', NF90_DOUBLE, id_var_eutm_xz(av,1), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_xz(av), &
+                                 (/ id_dim_y_xz(av) /),     &
+                                 'Nu_UTM', NF90_DOUBLE, id_var_nutm_xz(av,1), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_xz(av), &
+                                 (/ id_dim_x_xz(av) /),     &
+                                 'Ev_UTM', NF90_DOUBLE, id_var_eutm_xz(av,2), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_xz(av), &
+                                 (/ id_dim_yv_xz(av) /),     &
+                                 'Nv_UTM', NF90_DOUBLE, id_var_nutm_xz(av,2), &
+                                 'm', '', 000, 000, 000 )
+          ELSE
+             CALL netcdf_create_var( id_set_xz(av), &
+                                 (/ id_dim_x_xz(av), id_dim_y_xz(av) /),      &
+                                 'E_UTM', NF90_DOUBLE, id_var_eutm_xz(av,0),  &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_xz(av), &
+                                 (/ id_dim_x_xz(av), id_dim_y_xz(av) /),      &
+                                 'N_UTM', NF90_DOUBLE, id_var_nutm_xz(av,0),  &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_xz(av), &
+                                 (/ id_dim_xu_xz(av), id_dim_y_xz(av) /),     &
+                                 'Eu_UTM', NF90_DOUBLE, id_var_eutm_xz(av,1), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_xz(av), &
+                                 (/ id_dim_xu_xz(av), id_dim_y_xz(av) /),     &
+                                 'Nu_UTM', NF90_DOUBLE, id_var_nutm_xz(av,1), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_xz(av), &
+                                 (/ id_dim_x_xz(av), id_dim_yv_xz(av) /),     &
+                                 'Ev_UTM', NF90_DOUBLE, id_var_eutm_xz(av,2), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_xz(av), &
+                                 (/ id_dim_x_xz(av), id_dim_yv_xz(av) /),     &
+                                 'Nv_UTM', NF90_DOUBLE, id_var_nutm_xz(av,2), &
+                                 'm', '', 000, 000, 000 )
+          ENDIF
 
           IF ( land_surface )  THEN
@@ -3011,8 +3615,9 @@
 !
 !--                Most variables are defined on the zu grid
-                   CASE ( 'e_xz', 'lpt_xz', 'nc_xz', 'nr_xz', 'p_xz', 'pc_xz', &
-                          'pr_xz', 'prr_xz', 'pt_xz', 'q_xz', 'qc_xz',         &
+                   CASE ( 'e_xz', 'nc_xz', 'nr_xz', 'p_xz', 'pc_xz',           &
+                          'pr_xz', 'prr_xz', 'q_xz', 'qc_xz',                  &
                           'ql_xz', 'ql_c_xz', 'ql_v_xz', 'ql_vp_xz', 'qr_xz',  &
-                          'qv_xz', 'rho_ocean_xz', 's_xz', 'sa_xz', 'vpt_xz' )
+                          'qv_xz', 's_xz',                                     &
+                          'theta_xz', 'thetal_xz', 'thetav_xz'                 )
 
                       grid_x = 'x'
@@ -3055,4 +3660,10 @@
                       ENDIF
 
+!
+!--                   Check for ocean quantities
+                      IF ( .NOT. found  .AND.  ocean_mode )  THEN
+                         CALL ocean_define_netcdf_grid( do2d(av,i), found,  &
+                                                        grid_x, grid_y, grid_z )
+                      ENDIF
 !
 !--                   Check for radiation quantities
@@ -3294,6 +3905,135 @@
              ENDIF
 
-
              DEALLOCATE( netcdf_data )
+!
+!--          Write UTM coordinates
+             IF ( init_model%rotation_angle == 0.0_wp )  THEN
+!
+!--             1D in case of no rotation
+                cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
+!
+!--             x coordinates
+                ALLOCATE( netcdf_data(0:nx) )
+                DO  k = 0, 2
+!               
+!--                Scalar grid points
+                   IF ( k == 0 )  THEN
+                      shift_x = 0.5
+!               
+!--                u grid points
+                   ELSEIF ( k == 1 )  THEN
+                      shift_x = 0.0
+!               
+!--                v grid points
+                   ELSEIF ( k == 2 )  THEN
+                      shift_x = 0.5
+                   ENDIF
+                
+                   DO  i = 0, nx
+                     netcdf_data(i) = init_model%origin_x            &
+                                    + cos_ra * ( i + shift_x ) * dx
+                   ENDDO
+                
+                   nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_eutm_xz(av,k),&
+                                           netcdf_data, start = (/ 1 /),   &
+                                           count = (/ nx+1 /) )
+                   CALL netcdf_handle_error( 'netcdf_define_header', 555 )
+
+                ENDDO
+                DEALLOCATE( netcdf_data )
+!
+!--             y coordinates
+                ALLOCATE( netcdf_data(1:ns) )
+                DO  k = 0, 2
+!
+!--                Scalar grid points
+                   IF ( k == 0 )  THEN
+                      shift_y = 0.5
+!
+!--                u grid points
+                   ELSEIF ( k == 1 )  THEN
+                      shift_y = 0.5
+!
+!--                v grid points
+                   ELSEIF ( k == 2 )  THEN
+                      shift_y = 0.0
+                   ENDIF
+
+                   DO  i = 1, ns
+                      IF( section(i,2) == -1 )  THEN
+                         netcdf_data(i) = -1.0_wp  ! section averaged along y
+                      ELSE
+                         netcdf_data(i) = init_model%origin_y &
+                                     + cos_ra * ( section(i,2) + shift_y ) * dy
+                      ENDIF
+                   ENDDO
+
+                   nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_nutm_xz(av,k),&
+                                           netcdf_data, start = (/ 1 /),   &
+                                           count = (/ ns /) )
+                   CALL netcdf_handle_error( 'netcdf_define_header', 556 )
+
+                ENDDO
+                DEALLOCATE( netcdf_data )
+
+             ELSE
+!
+!--             2D in case of rotation
+                ALLOCATE( netcdf_data_2d(0:nx,1:ns) )
+                cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
+                sin_ra = SIN( init_model%rotation_angle * pi / 180.0_wp )
+                
+                DO  k = 0, 2
+!               
+!--                Scalar grid points
+                   IF ( k == 0 )  THEN
+                      shift_x = 0.5 ; shift_y = 0.5
+!                  
+!--                u grid points
+                   ELSEIF ( k == 1 )  THEN
+                      shift_x = 0.0 ; shift_y = 0.5
+!                  
+!--                v grid points
+                   ELSEIF ( k == 2 )  THEN
+                      shift_x = 0.5 ; shift_y = 0.0
+                   ENDIF
+
+                   DO  j = 1, ns
+                      IF( section(j,2) == -1 )  THEN
+                         netcdf_data_2d(:,j) = -1.0_wp  ! section averaged along y
+                      ELSE
+                         DO  i = 0, nx
+                            netcdf_data_2d(i,j) = init_model%origin_x          &
+                                    + cos_ra * ( i + shift_x ) * dx            &
+                                    + sin_ra * ( section(j,2) + shift_y ) * dy
+                         ENDDO
+                      ENDIF
+                   ENDDO
+                   
+                   nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_eutm_xz(av,k),  &
+                                           netcdf_data_2d, start = (/ 1, 1 /),   &
+                                           count = (/ nx+1, ns /) )
+                   CALL netcdf_handle_error( 'netcdf_define_header', 555 )
+                   
+                   DO  j = 1, ns
+                      IF( section(j,2) == -1 )  THEN
+                         netcdf_data_2d(:,j) = -1.0_wp  ! section averaged along y
+                      ELSE
+                         DO  i = 0, nx
+                            netcdf_data_2d(i,j) = init_model%origin_y          &
+                                    - sin_ra * ( i + shift_x ) * dx            &
+                                    + cos_ra * ( section(j,2) + shift_y ) * dy
+                         ENDDO
+                      ENDIF
+                   ENDDO
+                   
+                   nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_nutm_xz(av,k),  &
+                                           netcdf_data_2d, start = (/ 1, 1 /),   &
+                                           count = (/ nx+1, ns /) )
+                   CALL netcdf_handle_error( 'netcdf_define_header', 556 )
+                
+                ENDDO
+                DEALLOCATE( netcdf_data_2d )
+             ENDIF
 
           ENDIF
@@ -3654,4 +4394,57 @@
                                   NF90_DOUBLE, id_var_zw_yz(av), 'meters', '', &
                                   196, 197, 000 )
+!
+!--       Define UTM coordinates
+          IF ( init_model%rotation_angle == 0.0_wp )  THEN
+             CALL netcdf_create_var( id_set_yz(av), &
+                                 (/ id_dim_x_yz(av) /),      &
+                                 'E_UTM', NF90_DOUBLE, id_var_eutm_yz(av,0),  &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_yz(av), &
+                                 (/ id_dim_y_yz(av) /),      &
+                                 'N_UTM', NF90_DOUBLE, id_var_nutm_yz(av,0),  &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_yz(av), &
+                                 (/ id_dim_xu_yz(av) /),     &
+                                 'Eu_UTM', NF90_DOUBLE, id_var_eutm_yz(av,1), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_yz(av), &
+                                 (/ id_dim_y_yz(av) /),     &
+                                 'Nu_UTM', NF90_DOUBLE, id_var_nutm_yz(av,1), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_yz(av), &
+                                 (/ id_dim_x_yz(av) /),     &
+                                 'Ev_UTM', NF90_DOUBLE, id_var_eutm_yz(av,2), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_yz(av), &
+                                 (/ id_dim_yv_yz(av) /),     &
+                                 'Nv_UTM', NF90_DOUBLE, id_var_nutm_yz(av,2), &
+                                 'm', '', 000, 000, 000 )
+          ELSE
+             CALL netcdf_create_var( id_set_yz(av), &
+                                 (/ id_dim_x_yz(av), id_dim_y_yz(av) /),      &
+                                 'E_UTM', NF90_DOUBLE, id_var_eutm_yz(av,0),  &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_yz(av), &
+                                 (/ id_dim_x_yz(av), id_dim_y_yz(av) /),      &
+                                 'N_UTM', NF90_DOUBLE, id_var_nutm_yz(av,0),  &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_yz(av), &
+                                 (/ id_dim_xu_yz(av), id_dim_y_yz(av) /),     &
+                                 'Eu_UTM', NF90_DOUBLE, id_var_eutm_yz(av,1), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_yz(av), &
+                                 (/ id_dim_xu_yz(av), id_dim_y_yz(av) /),     &
+                                 'Nu_UTM', NF90_DOUBLE, id_var_nutm_yz(av,1), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_yz(av), &
+                                 (/ id_dim_x_yz(av), id_dim_yv_yz(av) /),     &
+                                 'Ev_UTM', NF90_DOUBLE, id_var_eutm_yz(av,2), &
+                                 'm', '', 000, 000, 000 )
+             CALL netcdf_create_var( id_set_yz(av), &
+                                 (/ id_dim_x_yz(av), id_dim_yv_yz(av) /),     &
+                                 'Nv_UTM', NF90_DOUBLE, id_var_nutm_yz(av,2), &
+                                 'm', '', 000, 000, 000 )
+          ENDIF
 
           IF ( land_surface )  THEN
@@ -3680,8 +4473,9 @@
 !
 !--                Most variables are defined on the zu grid
-                   CASE ( 'e_yz', 'lpt_yz', 'nc_yz', 'nr_yz', 'p_yz', 'pc_yz', &
-                          'pr_yz','prr_yz', 'pt_yz', 'q_yz', 'qc_yz', 'ql_yz', &
+                   CASE ( 'e_yz', 'nc_yz', 'nr_yz', 'p_yz', 'pc_yz',           &
+                          'pr_yz','prr_yz', 'q_yz', 'qc_yz', 'ql_yz',          &
                           'ql_c_yz', 'ql_v_yz', 'ql_vp_yz', 'qr_yz', 'qv_yz',  &
-                          'rho_ocean_yz', 's_yz', 'sa_yz', 'vpt_yz' )
+                          's_yz',                                              &
+                          'theta_yz', 'thetal_yz', 'thetav_yz' )
 
                       grid_x = 'x'
@@ -3724,4 +4518,10 @@
                       ENDIF
 
+!
+!--                   Check for ocean quantities
+                      IF ( .NOT. found  .AND.  ocean_mode )  THEN
+                         CALL ocean_define_netcdf_grid( do2d(av,i), found,     &
+                                                       grid_x, grid_y, grid_z )
+                      ENDIF
 !
 !--                   Check for radiation quantities
@@ -3951,4 +4751,134 @@
 
              DEALLOCATE( netcdf_data )
+!
+!--          Write UTM coordinates
+             IF ( init_model%rotation_angle == 0.0_wp )  THEN
+!
+!--             1D in case of no rotation
+                cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
+!
+!--             x coordinates
+                ALLOCATE( netcdf_data(1:ns) )
+                DO  k = 0, 2
+!               
+!--                Scalar grid points
+                   IF ( k == 0 )  THEN
+                      shift_x = 0.5
+!               
+!--                u grid points
+                   ELSEIF ( k == 1 )  THEN
+                      shift_x = 0.0
+!               
+!--                v grid points
+                   ELSEIF ( k == 2 )  THEN
+                      shift_x = 0.5
+                   ENDIF
+                
+                   DO  i = 1, ns
+                      IF( section(i,3) == -1 )  THEN
+                         netcdf_data(i) = -1.0_wp  ! section averaged along x
+                      ELSE
+                         netcdf_data(i) = init_model%origin_x &
+                                     + cos_ra * ( section(i,3) + shift_x ) * dx
+                      ENDIF
+                   ENDDO
+                
+                   nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_eutm_yz(av,k),&
+                                           netcdf_data, start = (/ 1 /),   &
+                                           count = (/ ns /) )
+                   CALL netcdf_handle_error( 'netcdf_define_header', 555 )
+
+                ENDDO
+                DEALLOCATE( netcdf_data )
+!
+!--             y coordinates
+                ALLOCATE( netcdf_data(0:ny) )
+                DO  k = 0, 2
+!
+!--                Scalar grid points
+                   IF ( k == 0 )  THEN
+                      shift_y = 0.5
+!
+!--                u grid points
+                   ELSEIF ( k == 1 )  THEN
+                      shift_y = 0.5
+!
+!--                v grid points
+                   ELSEIF ( k == 2 )  THEN
+                      shift_y = 0.0
+                   ENDIF
+
+                   DO  i = 0, ny
+                     netcdf_data(i) = init_model%origin_y            &
+                                    + cos_ra * ( i + shift_y ) * dy
+                   ENDDO
+
+                   nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_nutm_yz(av,k),&
+                                           netcdf_data, start = (/ 1 /),   &
+                                           count = (/ ny+1 /) )
+                   CALL netcdf_handle_error( 'netcdf_define_header', 556 )
+
+                ENDDO
+                DEALLOCATE( netcdf_data )
+
+             ELSE
+!
+!--             2D in case of rotation
+                ALLOCATE( netcdf_data_2d(1:ns,0:ny) )
+                cos_ra = COS( init_model%rotation_angle * pi / 180.0_wp )
+                sin_ra = SIN( init_model%rotation_angle * pi / 180.0_wp )
+
+                DO  k = 0, 2
+!               
+!--                Scalar grid points
+                   IF ( k == 0 )  THEN
+                      shift_x = 0.5 ; shift_y = 0.5
+!                  
+!--                u grid points
+                   ELSEIF ( k == 1 )  THEN
+                      shift_x = 0.0 ; shift_y = 0.5
+!                  
+!--                v grid points
+                   ELSEIF ( k == 2 )  THEN
+                      shift_x = 0.5 ; shift_y = 0.0
+                   ENDIF
+
+                   DO  j = 0, ny
+                      DO  i = 1, ns
+                         IF( section(i,3) == -1 )  THEN
+                            netcdf_data_2d(i,:) = -1.0_wp !section averaged along x
+                         ELSE
+                            netcdf_data_2d(i,j) = init_model%origin_x          &
+                                    + cos_ra * ( section(i,3) + shift_x ) * dx &
+                                    + sin_ra * ( j + shift_y ) * dy
+                         ENDIF
+                      ENDDO
+                   ENDDO
+                   
+                   nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_eutm_yz(av,k),  &
+                                           netcdf_data_2d, start = (/ 1, 1 /),   &
+                                           count = (/ ns, ny+1 /) )
+                   CALL netcdf_handle_error( 'netcdf_define_header', 555 )
+                   
+                   DO  j = 0, ny
+                      DO  i = 1, ns
+                         IF( section(i,3) == -1 )  THEN
+                            netcdf_data_2d(i,:) = -1.0_wp !section averaged along x
+                         ELSE
+                            netcdf_data_2d(i,j) = init_model%origin_y          &
+                                    - sin_ra * ( section(i,3) + shift_x ) * dx &
+                                    + cos_ra * ( j + shift_y ) * dy
+                         ENDIF
+                      ENDDO
+                   ENDDO
+
+                   nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_nutm_yz(av,k),  &
+                                           netcdf_data_2d, start = (/ 1, 1 /),   &
+                                           count = (/ ns, ny+1 /) )
+                   CALL netcdf_handle_error( 'netcdf_define_header', 556 )
+                
+                ENDDO
+                DEALLOCATE( netcdf_data_2d )
+             ENDIF
 
           ENDIF
@@ -4871,7 +5801,8 @@
 !
 !--             Most variables are defined on the zu levels
-                CASE ( 'e', 'lpt', 'nc', 'nr', 'p', 'pc', 'pr', 'prr', 'pt',   &
+                CASE ( 'e', 'nc', 'nr', 'p', 'pc', 'pr', 'prr',   &
                        'q', 'qc', 'ql', 'ql_c', 'ql_v', 'ql_vp', 'qr', 'qv',   &
-                       'rho_ocean', 's', 'sa', 'u', 'v', 'vpt' )
+                       'rho_sea_water', 's', 'sa', &
+                       'theta', 'thetal', 'thetav', 'u', 'v' )
 
                    grid_z = 'zu'
