Index: /palm/trunk/SOURCE/netcdf_data_input_mod.f90
===================================================================
--- /palm/trunk/SOURCE/netcdf_data_input_mod.f90	(revision 3458)
+++ /palm/trunk/SOURCE/netcdf_data_input_mod.f90	(revision 3459)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! -----------------
-! 
+! Define coordinate reference system (crs) and read from input dataset
+! Revise default values for reference coordinates
 ! 
 ! Former revisions:
@@ -348,9 +349,9 @@
        REAL(wp) ::  fill_v                  !< fill value for v
        REAL(wp) ::  fill_w                  !< fill value for w
-       REAL(wp) ::  latitude                !< latitude of the southern model boundary
-       REAL(wp) ::  longitude               !< longitude of the western model boundary
-       REAL(wp) ::  origin_x = 0.0_wp       !< x position of the western model boundary (default: Hannover)
-       REAL(wp) ::  origin_y = 0.0_wp       !< y position of the northern model boundary (default: Hannover)
-       REAL(wp) ::  origin_z = 0.0_wp       !< reference height of input data (default: Hannover)
+       REAL(wp) ::  latitude = 0.0_wp       !< latitude of the lower left corner
+       REAL(wp) ::  longitude = -3.0_wp     !< longitude of the lower left corner
+       REAL(wp) ::  origin_x = 500000.0_wp  !< UTM easting of the lower left corner
+       REAL(wp) ::  origin_y = 0.0_wp       !< UTM northing of the lower left corner
+       REAL(wp) ::  origin_z = 0.0_wp       !< reference height of input data
        REAL(wp) ::  rotation_angle = 0.0_wp !< rotation angle of input data
 
@@ -585,6 +586,26 @@
     END TYPE global_atts_type
 !
+!-- Define type for coordinate reference system (crs)
+    TYPE crs_type
+       CHARACTER(LEN=200) ::  epsg_code = 'EPSG:25831'                   !< EPSG code
+       CHARACTER(LEN=200) ::  grid_mapping_name = 'transverse_mercator'  !< name of grid mapping
+       CHARACTER(LEN=200) ::  long_name = 'coordinate reference system'  !< name of variable crs
+       CHARACTER(LEN=200) ::  units = 'm'                                !< unit of crs
+
+       REAL(wp) ::  false_easting = 500000.0_wp                  !< false easting
+       REAL(wp) ::  false_northing = 0.0_wp                      !< false northing
+       REAL(wp) ::  inverse_flattening = 298.257223563_wp        !< 1/f (default for WGS84)
+       REAL(wp) ::  latitude_of_projection_origin = 0.0_wp       !< latitude of projection origin
+       REAL(wp) ::  longitude_of_central_meridian = 3.0_wp       !< longitude of central meridian of UTM zone (default: zone 31)
+       REAL(wp) ::  longitude_of_prime_meridian = 0.0_wp         !< longitude of prime meridian
+       REAL(wp) ::  scale_factor_at_central_meridian = 0.9996_wp !< scale factor of UTM coordinates
+       REAL(wp) ::  semi_major_axis = 6378137.0_wp               !< length of semi major axis (default for WGS84)
+    END TYPE crs_type
+
+!
 !-- Define variables
-    TYPE(dims_xy)    ::  dim_static  !< data structure for x, y-dimension in static input file
+    TYPE(crs_type)   ::  coord_ref_sys  !< coordinate reference system
+
+    TYPE(dims_xy)    ::  dim_static     !< data structure for x, y-dimension in static input file
 
     TYPE(nest_offl_type) ::  nest_offl  !< data structure for data input at lateral and top boundaries (provided by Inifor)  
@@ -740,4 +761,5 @@
            building_id_f, building_pars_f, building_type_f,                    &
            chem_emis, chem_emis_att, chem_emis_att_type, chem_emis_val_type,   &
+           coord_ref_sys,                                                      &
            init_3d, init_model, input_file_static, input_pids_static,          &
            input_pids_dynamic, leaf_area_density_f, nest_offl,                 &
@@ -796,5 +818,6 @@
 ! Description:
 ! ------------
-!> Reads global attributes required for initialization of the model.
+!> Reads global attributes and coordinate reference system required for 
+!> initialization of the model.
 !------------------------------------------------------------------------------!
     SUBROUTINE netcdf_data_input_init
@@ -802,5 +825,6 @@
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  id_mod   !< NetCDF id of input file
+       INTEGER(iwp) ::  id_mod     !< NetCDF id of input file
+       INTEGER(iwp) ::  var_id_crs !< NetCDF id of variable crs
 
        IF ( .NOT. input_pids_static )  RETURN
@@ -833,5 +857,45 @@
        CALL get_attribute( id_mod, input_file_atts%rotation_angle_char,        &
                            input_file_atts%rotation_angle, .TRUE. )
-
+!
+!--    Read coordinate reference system if available
+       nc_stat = NF90_INQ_VARID( id_mod, 'crs', var_id_crs )
+       IF ( nc_stat == NF90_NOERR )  THEN
+          CALL get_attribute( id_mod, 'epsg_code',                             &
+                              coord_ref_sys%epsg_code,                         &
+                              .FALSE., 'crs' )
+          CALL get_attribute( id_mod, 'false_easting',                         &
+                              coord_ref_sys%false_easting,                     &
+                              .FALSE., 'crs' )
+          CALL get_attribute( id_mod, 'false_northing',                        &
+                              coord_ref_sys%false_northing,                    &
+                              .FALSE., 'crs' )
+          CALL get_attribute( id_mod, 'grid_mapping_name',                     &
+                              coord_ref_sys%grid_mapping_name,                 &
+                              .FALSE., 'crs' )
+          CALL get_attribute( id_mod, 'inverse_flattening',                    &
+                              coord_ref_sys%inverse_flattening,                &
+                              .FALSE., 'crs' )
+          CALL get_attribute( id_mod, 'latitude_of_projection_origin',         &
+                              coord_ref_sys%latitude_of_projection_origin,     &
+                              .FALSE., 'crs' )
+          CALL get_attribute( id_mod, 'long_name',                             &
+                              coord_ref_sys%long_name,                         &
+                              .FALSE., 'crs' )
+          CALL get_attribute( id_mod, 'longitude_of_central_meridian',         &
+                              coord_ref_sys%longitude_of_central_meridian,     &
+                              .FALSE., 'crs' )
+          CALL get_attribute( id_mod, 'longitude_of_prime_meridian',           &
+                              coord_ref_sys%longitude_of_prime_meridian,       &
+                              .FALSE., 'crs' )
+          CALL get_attribute( id_mod, 'scale_factor_at_central_meridian',      &
+                              coord_ref_sys%scale_factor_at_central_meridian,  &
+                              .FALSE., 'crs' )
+          CALL get_attribute( id_mod, 'semi_major_axis',                       &
+                              coord_ref_sys%semi_major_axis,                   &
+                              .FALSE., 'crs' )
+          CALL get_attribute( id_mod, 'units',                                 &
+                              coord_ref_sys%units,                             &
+                              .FALSE., 'crs' )
+       ENDIF
 !
 !--    Finally, close input file
Index: /palm/trunk/SOURCE/netcdf_interface_mod.f90
===================================================================
--- /palm/trunk/SOURCE/netcdf_interface_mod.f90	(revision 3458)
+++ /palm/trunk/SOURCE/netcdf_interface_mod.f90	(revision 3459)
@@ -20,5 +20,8 @@
 ! Current revisions:
 ! ------------------
-! 
+! - Add variable crs to output files
+! - Add long_name to UTM coordinates
+! - Add latitude/longitude coordinates. For 3d and xy output, lon and lat are
+!   only written if parallel output is used.
 ! 
 ! Former revisions:
@@ -311,4 +314,6 @@
 !> @todo timeseries and profile output still needs to be rewritten to allow
 !>       modularization
+!> @todo output 2d UTM coordinates without global arrays
+!> @todo output longitude/latitude also with non-parallel output (3d and xy)
 !------------------------------------------------------------------------------!
  MODULE netcdf_interface
@@ -324,4 +329,7 @@
     USE mas_global_attributes,                                                 &
         ONLY:  dim_size_agtnum
+
+    USE netcdf_data_input_mod,                                                 &
+        ONLY: coord_ref_sys
 
     PRIVATE
@@ -472,6 +480,10 @@
                                          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
+                                         id_var_eutm_yz, id_var_nutm_yz
+
+    INTEGER(iwp), DIMENSION(0:1,0:2) ::  id_var_lat_3d, id_var_lon_3d, &
+                                         id_var_lat_xy, id_var_lon_xy, &
+                                         id_var_lat_xz, id_var_lon_xz, &
+                                         id_var_lat_yz, id_var_lon_yz
 
     INTEGER ::  netcdf_data_format = 2  !< NetCDF3 64bit offset format
@@ -520,4 +532,7 @@
     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,0:2) ::  id_var_lat_mask, &
+                                                     id_var_lon_mask
 
     INTEGER(iwp), DIMENSION(1:max_masks,0:1,100) ::  id_var_domask
@@ -624,5 +639,5 @@
 
     USE netcdf_data_input_mod,                                                 &
-        ONLY: init_model
+        ONLY:  init_model
 
     USE ocean_mod,                                                             &
@@ -720,4 +735,6 @@
 
     REAL(wp) ::  cos_ra                                      !< cosine of rotation_angle
+    REAL(wp) ::  eutm                                        !< easting (UTM)
+    REAL(wp) ::  nutm                                        !< northing (UTM)
     REAL(wp) ::  shift_x                                     !< shift of x coordinate
     REAL(wp) ::  shift_y                                     !< shift of y coordinate
@@ -725,7 +742,10 @@
 
     REAL(wp), DIMENSION(1) ::  last_time_coordinate          !< last time value in file
+    REAL(wp), DIMENSION(8) ::  crs_list                      !< list of coord_ref_sys values
 
     REAL(wp), DIMENSION(:), ALLOCATABLE   ::  netcdf_data    !<
     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  netcdf_data_2d !<
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lat            !< latitude
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  lon            !< longitude
 
 
@@ -822,4 +842,15 @@
 
 !
+!-- Convert coord_ref_sys into vector (used for lat/lon calculation)
+    crs_list = (/ coord_ref_sys%semi_major_axis,                  &
+                  coord_ref_sys%inverse_flattening,               &
+                  coord_ref_sys%longitude_of_prime_meridian,      &
+                  coord_ref_sys%longitude_of_central_meridian,    &
+                  coord_ref_sys%scale_factor_at_central_meridian, &
+                  coord_ref_sys%latitude_of_projection_origin,    &
+                  coord_ref_sys%false_easting,                    &
+                  coord_ref_sys%false_northing /)
+
+!
 !-- Determine the mode to be processed
     IF ( extend )  THEN
@@ -962,51 +993,80 @@
                                  (/ id_dim_x_mask(mid,av) /),      &
                                  'E_UTM', NF90_DOUBLE, id_var_eutm_mask(mid,av,0),  &
-                                 'm', '', 000, 000, 000 )
+                                 'm', 'projection_x_coordinate', 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 )
+                                 'm', 'projection_y_coordinate', 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 )
+                                 'm', 'projection_x_coordinate', 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 )
+                                 'm', 'projection_y_coordinate', 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 )
+                                 'm', 'projection_x_coordinate', 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 )
+                                 'm', 'projection_y_coordinate', 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 )
+                                 (/ id_dim_x_mask(mid,av), id_dim_y_mask(mid,av) /), &
+                                 'E_UTM', NF90_DOUBLE, id_var_eutm_mask(mid,av,0),   &
+                                 'm', 'projection_x_coordinate', 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 )
+                                 (/ id_dim_x_mask(mid,av), id_dim_y_mask(mid,av) /), &
+                                 'N_UTM', NF90_DOUBLE, id_var_nutm_mask(mid,av,0),   &
+                                 'm', 'projection_y_coordinate', 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 )
+                                 (/ id_dim_xu_mask(mid,av), id_dim_y_mask(mid,av) /),&
+                                 'Eu_UTM', NF90_DOUBLE, id_var_eutm_mask(mid,av,1),  &
+                                 'm', 'projection_x_coordinate', 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 )
+                                 (/ id_dim_xu_mask(mid,av), id_dim_y_mask(mid,av) /),&
+                                 'Nu_UTM', NF90_DOUBLE, id_var_nutm_mask(mid,av,1),  &
+                                 'm', 'projection_y_coordinate', 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 )
+                                 (/ id_dim_x_mask(mid,av), id_dim_yv_mask(mid,av) /),&
+                                 'Ev_UTM', NF90_DOUBLE, id_var_eutm_mask(mid,av,2),  &
+                                 'm', 'projection_x_coordinate', 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
+                                 (/ id_dim_x_mask(mid,av), id_dim_yv_mask(mid,av) /),&
+                                 'Nv_UTM', NF90_DOUBLE, id_var_nutm_mask(mid,av,2),  &
+                                 'm', 'projection_y_coordinate', 000, 000, 000 )
+          ENDIF
+!
+!--       Define geographic coordinates
+          CALL netcdf_create_var( id_set_mask(mid,av), &
+                              (/ id_dim_x_mask(mid,av), id_dim_y_mask(mid,av) /), &
+                              'lon', NF90_DOUBLE, id_var_lon_mask(mid,av,0),      &
+                              'degrees_east', 'longitude', 000, 000, 000 )
+          CALL netcdf_create_var( id_set_mask(mid,av), &
+                              (/ id_dim_x_mask(mid,av), id_dim_y_mask(mid,av) /), &
+                              'lat', NF90_DOUBLE, id_var_lat_mask(mid,av,0),      &
+                              'degrees_north', 'latitude', 000, 000, 000 )
+          CALL netcdf_create_var( id_set_mask(mid,av), &
+                              (/ id_dim_xu_mask(mid,av), id_dim_y_mask(mid,av) /),&
+                              'lonu', NF90_DOUBLE, id_var_lon_mask(mid,av,1),     &
+                              'degrees_east', 'longitude', 000, 000, 000 )
+          CALL netcdf_create_var( id_set_mask(mid,av), &
+                              (/ id_dim_xu_mask(mid,av), id_dim_y_mask(mid,av) /),&
+                              'latu', NF90_DOUBLE, id_var_lat_mask(mid,av,1),     &
+                              'degrees_north', 'latitude', 000, 000, 000 )
+          CALL netcdf_create_var( id_set_mask(mid,av), &
+                              (/ id_dim_x_mask(mid,av), id_dim_yv_mask(mid,av) /),&
+                              'lonv', NF90_DOUBLE, id_var_lon_mask(mid,av,2),     &
+                              'degrees_east', 'longitude', 000, 000, 000 )
+          CALL netcdf_create_var( id_set_mask(mid,av), &
+                              (/ id_dim_x_mask(mid,av), id_dim_yv_mask(mid,av) /),&
+                              'latv', NF90_DOUBLE, id_var_lat_mask(mid,av,2),     &
+                              'degrees_north', 'latitude', 000, 000, 000 )
+!
+!--       Define coordinate-reference system
+          CALL netcdf_create_crs( id_set_mask(mid,av), 000 )
 !
 !--       In case of non-flat topography define 2d-arrays containing the height
@@ -1364,5 +1424,58 @@
              DEALLOCATE( netcdf_data_2d )
           ENDIF
-
+!
+!--       Write lon and lat data.
+          ALLOCATE( lat(mask_size(mid,1),mask_size(mid,2)) )
+          ALLOCATE( lon(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)
+                   eutm = init_model%origin_x                               &
+                        + cos_ra * ( mask_i_global(mid,i) + shift_x ) * dx  &
+                        + sin_ra * ( mask_j_global(mid,j) + shift_y ) * dy
+                   nutm = init_model%origin_y                               &
+                        - sin_ra * ( mask_i_global(mid,i) + shift_x ) * dx  &
+                        + cos_ra * ( mask_j_global(mid,j) + shift_y ) * dy
+
+                   CALL  convert_utm_to_geographic( crs_list,          &
+                                                    eutm, nutm,        &
+                                                    lon(i,j), lat(i,j) )
+                ENDDO
+             ENDDO
+
+             nc_stat = NF90_PUT_VAR( id_set_mask(mid,av),           &
+                                     id_var_lon_mask(mid,av,k),     &
+                                     lon, start = (/ 1, 1 /),       &
+                                     count = (/ mask_size(mid,1),   &
+                                                mask_size(mid,2) /) )
+             CALL netcdf_handle_error( 'netcdf_define_header', 556 )
+
+             nc_stat = NF90_PUT_VAR( id_set_mask(mid,av),           &
+                                     id_var_lat_mask(mid,av,k),     &
+                                     lat, start = (/ 1, 1 /),       &
+                                     count = (/ mask_size(mid,1),   &
+                                                mask_size(mid,2) /) )
+             CALL netcdf_handle_error( 'netcdf_define_header', 556 )
+          ENDDO
+
+          DEALLOCATE( lat )
+          DEALLOCATE( lon )
 !
 !--       Write zu and zw data (vertical axes)
@@ -1703,51 +1816,80 @@
                                  (/ id_dim_x_3d(av) /),      &
                                  'E_UTM', NF90_DOUBLE, id_var_eutm_3d(av,0),  &
-                                 'm', '', 000, 000, 000 )
+                                 'm', 'projection_x_coordinate', 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 )
+                                 'm', 'projection_y_coordinate', 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 )
+                                 'm', 'projection_x_coordinate', 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 )
+                                 'm', 'projection_y_coordinate', 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 )
+                                 'm', 'projection_x_coordinate', 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 )
+                                 'm', 'projection_y_coordinate', 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 )
+                                 'm', 'projection_x_coordinate', 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 )
+                                 'm', 'projection_y_coordinate', 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 )
+                                 'm', 'projection_x_coordinate', 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 )
+                                 'm', 'projection_y_coordinate', 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 )
+                                 'm', 'projection_x_coordinate', 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
+                                 'm', 'projection_y_coordinate', 000, 000, 000 )
+          ENDIF
+!
+!--       Define geographic coordinates
+          CALL netcdf_create_var( id_set_3d(av), &
+                              (/ id_dim_x_3d(av), id_dim_y_3d(av) /),      &
+                              'lon', NF90_DOUBLE, id_var_lon_3d(av,0),  &
+                              'degrees_east', 'longitude', 000, 000, 000 )
+          CALL netcdf_create_var( id_set_3d(av), &
+                              (/ id_dim_x_3d(av), id_dim_y_3d(av) /),      &
+                              'lat', NF90_DOUBLE, id_var_lat_3d(av,0),  &
+                              'degrees_north', 'latitude', 000, 000, 000 )
+          CALL netcdf_create_var( id_set_3d(av), &
+                              (/ id_dim_xu_3d(av), id_dim_y_3d(av) /),     &
+                              'lonu', NF90_DOUBLE, id_var_lon_3d(av,1), &
+                              'degrees_east', 'longitude', 000, 000, 000 )
+          CALL netcdf_create_var( id_set_3d(av), &
+                              (/ id_dim_xu_3d(av), id_dim_y_3d(av) /),     &
+                              'latu', NF90_DOUBLE, id_var_lat_3d(av,1), &
+                              'degrees_north', 'latitude', 000, 000, 000 )
+          CALL netcdf_create_var( id_set_3d(av), &
+                              (/ id_dim_x_3d(av), id_dim_yv_3d(av) /),     &
+                              'lonv', NF90_DOUBLE, id_var_lon_3d(av,2), &
+                              'degrees_east', 'longitude', 000, 000, 000 )
+          CALL netcdf_create_var( id_set_3d(av), &
+                              (/ id_dim_x_3d(av), id_dim_yv_3d(av) /),     &
+                              'latv', NF90_DOUBLE, id_var_lat_3d(av,2), &
+                              'degrees_north', 'latitude', 000, 000, 000 )
+!
+!--       Define coordinate-reference system
+          CALL netcdf_create_crs( id_set_3d(av), 000 )
 !
 !--       In case of non-flat topography define 2d-arrays containing the height
@@ -2168,4 +2310,58 @@
                 CALL netcdf_handle_error( 'netcdf_define_header', 86 )
              ENDIF
+
+          ENDIF
+!
+!--       Write lon and lat data. Only for parallel output.
+          IF ( netcdf_data_format > 4 )  THEN
+
+             ALLOCATE( lat(nxl:nxr,nys:nyn) )
+             ALLOCATE( lon(nxl:nxr,nys:nyn) )
+             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 = nys, nyn
+                   DO  i = nxl, nxr
+                      eutm = init_model%origin_x            &
+                           + cos_ra * ( i + shift_x ) * dx  &
+                           + sin_ra * ( j + shift_y ) * dy
+                      nutm = init_model%origin_y            &
+                           - sin_ra * ( i + shift_x ) * dx  &
+                           + cos_ra * ( j + shift_y ) * dy
+
+                      CALL  convert_utm_to_geographic( crs_list,          &
+                                                       eutm, nutm,        &
+                                                       lon(i,j), lat(i,j) )
+                   ENDDO
+                ENDDO
+
+                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_lon_3d(av,k), &
+                                     lon, start = (/ nxl+1, nys+1 /),       &
+                                     count = (/ nxr-nxl+1, nyn-nys+1 /) )
+                CALL netcdf_handle_error( 'netcdf_define_header', 556 )
+
+                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_lat_3d(av,k), &
+                                     lat, start = (/ nxl+1, nys+1 /),       &
+                                     count = (/ nxr-nxl+1, nyn-nys+1 /) )
+                CALL netcdf_handle_error( 'netcdf_define_header', 556 )
+             ENDDO
+
+             DEALLOCATE( lat )
+             DEALLOCATE( lon )
 
           ENDIF
@@ -2673,51 +2869,80 @@
                                  (/ id_dim_x_xy(av) /),      &
                                  'E_UTM', NF90_DOUBLE, id_var_eutm_xy(av,0),  &
-                                 'm', '', 000, 000, 000 )
+                                 'm', 'projection_x_coordinate', 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 )
+                                 'm', 'projection_y_coordinate', 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 )
+                                 'm', 'projection_x_coordinate', 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 )
+                                 'm', 'projection_y_coordinate', 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 )
+                                 'm', 'projection_x_coordinate', 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 )
+                                 'm', 'projection_y_coordinate', 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 )
+                                 'm', 'projection_x_coordinate', 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 )
+                                 'm', 'projection_y_coordinate', 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 )
+                                 'm', 'projection_x_coordinate', 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 )
+                                 'm', 'projection_y_coordinate', 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 )
+                                 'm', 'projection_x_coordinate', 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
+                                 'm', 'projection_y_coordinate', 000, 000, 000 )
+          ENDIF
+!
+!--       Define geographic coordinates
+          CALL netcdf_create_var( id_set_xy(av), &
+                              (/ id_dim_x_xy(av), id_dim_y_xy(av) /),      &
+                              'lon', NF90_DOUBLE, id_var_lon_xy(av,0),  &
+                              'degrees_east', 'longitude', 000, 000, 000 )
+          CALL netcdf_create_var( id_set_xy(av), &
+                              (/ id_dim_x_xy(av), id_dim_y_xy(av) /),      &
+                              'lat', NF90_DOUBLE, id_var_lat_xy(av,0),  &
+                              'degrees_north', 'latitude', 000, 000, 000 )
+          CALL netcdf_create_var( id_set_xy(av), &
+                              (/ id_dim_xu_xy(av), id_dim_y_xy(av) /),     &
+                              'lonu', NF90_DOUBLE, id_var_lon_xy(av,1), &
+                              'degrees_east', 'longitude', 000, 000, 000 )
+          CALL netcdf_create_var( id_set_xy(av), &
+                              (/ id_dim_xu_xy(av), id_dim_y_xy(av) /),     &
+                              'latu', NF90_DOUBLE, id_var_lat_xy(av,1), &
+                              'degrees_north', 'latitude', 000, 000, 000 )
+          CALL netcdf_create_var( id_set_xy(av), &
+                              (/ id_dim_x_xy(av), id_dim_yv_xy(av) /),     &
+                              'lonv', NF90_DOUBLE, id_var_lon_xy(av,2), &
+                              'degrees_east', 'longitude', 000, 000, 000 )
+          CALL netcdf_create_var( id_set_xy(av), &
+                              (/ id_dim_x_xy(av), id_dim_yv_xy(av) /),     &
+                              'latv', NF90_DOUBLE, id_var_lat_xy(av,2), &
+                              'degrees_north', 'latitude', 000, 000, 000 )
+!
+!--       Define coordinate-reference system
+          CALL netcdf_create_crs( id_set_xy(av), 000 )
 !
 !--       In case of non-flat topography define 2d-arrays containing the height
@@ -3199,4 +3424,58 @@
           ENDIF
 !
+!--       Write lon and lat data. Only for parallel output.
+          IF ( netcdf_data_format > 4 )  THEN
+
+             ALLOCATE( lat(nxl:nxr,nys:nyn) )
+             ALLOCATE( lon(nxl:nxr,nys:nyn) )
+             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 = nys, nyn
+                   DO  i = nxl, nxr
+                      eutm = init_model%origin_x            &
+                           + cos_ra * ( i + shift_x ) * dx  &
+                           + sin_ra * ( j + shift_y ) * dy
+                      nutm = init_model%origin_y            &
+                           - sin_ra * ( i + shift_x ) * dx  &
+                           + cos_ra * ( j + shift_y ) * dy
+
+                      CALL  convert_utm_to_geographic( crs_list,          &
+                                                       eutm, nutm,        &
+                                                       lon(i,j), lat(i,j) )
+                   ENDDO
+                ENDDO
+
+                nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_lon_xy(av,k), &
+                                     lon, start = (/ nxl+1, nys+1 /),       &
+                                     count = (/ nxr-nxl+1, nyn-nys+1 /) )
+                CALL netcdf_handle_error( 'netcdf_define_header', 556 )
+
+                nc_stat = NF90_PUT_VAR( id_set_xy(av), id_var_lat_xy(av,k), &
+                                     lat, start = (/ nxl+1, nys+1 /),       &
+                                     count = (/ nxr-nxl+1, nyn-nys+1 /) )
+                CALL netcdf_handle_error( 'netcdf_define_header', 556 )
+             ENDDO
+
+             DEALLOCATE( lat )
+             DEALLOCATE( lon )
+
+          ENDIF
+!
 !--       In case of non-flat topography write height information. Only for 
 !--       parallel netcdf output.
@@ -3595,51 +3874,80 @@
                                  (/ id_dim_x_xz(av) /),      &
                                  'E_UTM', NF90_DOUBLE, id_var_eutm_xz(av,0),  &
-                                 'm', '', 000, 000, 000 )
+                                 'm', 'projection_x_coordinate', 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 )
+                                 'm', 'projection_y_coordinate', 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 )
+                                 'm', 'projection_x_coordinate', 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 )
+                                 'm', 'projection_y_coordinate', 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 )
+                                 'm', 'projection_x_coordinate', 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 )
+                                 'm', 'projection_y_coordinate', 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 )
+                                 'm', 'projection_x_coordinate', 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 )
+                                 'm', 'projection_y_coordinate', 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 )
+                                 'm', 'projection_x_coordinate', 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 )
+                                 'm', 'projection_y_coordinate', 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 )
+                                 'm', 'projection_x_coordinate', 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
+                                 'm', 'projection_y_coordinate', 000, 000, 000 )
+          ENDIF
+!
+!--       Define geographic coordinates
+          CALL netcdf_create_var( id_set_xz(av), &
+                              (/ id_dim_x_xz(av), id_dim_y_xz(av) /),      &
+                              'lon', NF90_DOUBLE, id_var_lon_xz(av,0),  &
+                              'degrees_east', 'longitude', 000, 000, 000 )
+          CALL netcdf_create_var( id_set_xz(av), &
+                              (/ id_dim_x_xz(av), id_dim_y_xz(av) /),      &
+                              'lat', NF90_DOUBLE, id_var_lat_xz(av,0),  &
+                              'degrees_north', 'latitude', 000, 000, 000 )
+          CALL netcdf_create_var( id_set_xz(av), &
+                              (/ id_dim_xu_xz(av), id_dim_y_xz(av) /),     &
+                              'lonu', NF90_DOUBLE, id_var_lon_xz(av,1), &
+                              'degrees_east', 'longitude', 000, 000, 000 )
+          CALL netcdf_create_var( id_set_xz(av), &
+                              (/ id_dim_xu_xz(av), id_dim_y_xz(av) /),     &
+                              'latu', NF90_DOUBLE, id_var_lat_xz(av,1), &
+                              'degrees_north', 'latitude', 000, 000, 000 )
+          CALL netcdf_create_var( id_set_xz(av), &
+                              (/ id_dim_x_xz(av), id_dim_yv_xz(av) /),     &
+                              'lonv', NF90_DOUBLE, id_var_lon_xz(av,2), &
+                              'degrees_east', 'longitude', 000, 000, 000 )
+          CALL netcdf_create_var( id_set_xz(av), &
+                              (/ id_dim_x_xz(av), id_dim_yv_xz(av) /),     &
+                              'latv', NF90_DOUBLE, id_var_lat_xz(av,2), &
+                              'degrees_north', 'latitude', 000, 000, 000 )
+!
+!--       Define coordinate-reference system
+          CALL netcdf_create_crs( id_set_xz(av), 000 )
 
           IF ( land_surface )  THEN
@@ -4080,4 +4388,59 @@
                 DEALLOCATE( netcdf_data_2d )
              ENDIF
+!
+!--          Write lon and lat data
+             ALLOCATE( lat(0:nx,1:ns) )
+             ALLOCATE( lon(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
+                      lat(:,j) = -90.0_wp  ! section averaged along y
+                      lon(:,j) = -180.0_wp  ! section averaged along y
+                   ELSE
+                      DO  i = 0, nx
+                         eutm = init_model%origin_x            &
+                              + cos_ra * ( i + shift_x ) * dx  &
+                              + sin_ra * ( section(j,2) + shift_y ) * dy
+                         nutm = init_model%origin_y            &
+                              - sin_ra * ( i + shift_x ) * dx  &
+                              + cos_ra * ( section(j,2) + shift_y ) * dy
+
+                         CALL  convert_utm_to_geographic( crs_list,          &
+                                                          eutm, nutm,        &
+                                                          lon(i,j), lat(i,j) )
+                      ENDDO
+                   ENDIF
+                ENDDO
+
+                nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_lon_xz(av,k), &
+                                     lon, start = (/ 1, 1 /),       &
+                                     count = (/ nx+1, ns /) )
+                CALL netcdf_handle_error( 'netcdf_define_header', 556 )
+
+                nc_stat = NF90_PUT_VAR( id_set_xz(av), id_var_lat_xz(av,k), &
+                                     lat, start = (/ 1, 1 /),       &
+                                     count = (/ nx+1, ns /) )
+                CALL netcdf_handle_error( 'netcdf_define_header', 556 )
+             ENDDO
+
+             DEALLOCATE( lat )
+             DEALLOCATE( lon )
 
           ENDIF
@@ -4444,51 +4807,80 @@
                                  (/ id_dim_x_yz(av) /),      &
                                  'E_UTM', NF90_DOUBLE, id_var_eutm_yz(av,0),  &
-                                 'm', '', 000, 000, 000 )
+                                 'm', 'projection_x_coordinate', 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 )
+                                 'm', 'projection_y_coordinate', 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 )
+                                 'm', 'projection_x_coordinate', 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 )
+                                 'm', 'projection_y_coordinate', 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 )
+                                 'm', 'projection_x_coordinate', 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 )
+                                 'm', 'projection_y_coordinate', 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 )
+                                 'm', 'projection_x_coordinate', 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 )
+                                 'm', 'projection_y_coordinate', 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 )
+                                 'm', 'projection_x_coordinate', 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 )
+                                 'm', 'projection_y_coordinate', 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 )
+                                 'm', 'projection_x_coordinate', 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
+                                 'm', 'projection_y_coordinate', 000, 000, 000 )
+          ENDIF
+!
+!--       Define geographic coordinates
+          CALL netcdf_create_var( id_set_yz(av), &
+                              (/ id_dim_x_yz(av), id_dim_y_yz(av) /),      &
+                              'lon', NF90_DOUBLE, id_var_lon_yz(av,0),  &
+                              'degrees_east', 'longitude', 000, 000, 000 )
+          CALL netcdf_create_var( id_set_yz(av), &
+                              (/ id_dim_x_yz(av), id_dim_y_yz(av) /),      &
+                              'lat', NF90_DOUBLE, id_var_lat_yz(av,0),  &
+                              'degrees_north', 'latitude', 000, 000, 000 )
+          CALL netcdf_create_var( id_set_yz(av), &
+                              (/ id_dim_xu_yz(av), id_dim_y_yz(av) /),     &
+                              'lonu', NF90_DOUBLE, id_var_lon_yz(av,1), &
+                              'degrees_east', 'longitude', 000, 000, 000 )
+          CALL netcdf_create_var( id_set_yz(av), &
+                              (/ id_dim_xu_yz(av), id_dim_y_yz(av) /),     &
+                              'latu', NF90_DOUBLE, id_var_lat_yz(av,1), &
+                              'degrees_north', 'latitude', 000, 000, 000 )
+          CALL netcdf_create_var( id_set_yz(av), &
+                              (/ id_dim_x_yz(av), id_dim_yv_yz(av) /),     &
+                              'lonv', NF90_DOUBLE, id_var_lon_yz(av,2), &
+                              'degrees_east', 'longitude', 000, 000, 000 )
+          CALL netcdf_create_var( id_set_yz(av), &
+                              (/ id_dim_x_yz(av), id_dim_yv_yz(av) /),     &
+                              'latv', NF90_DOUBLE, id_var_lat_yz(av,2), &
+                              'degrees_north', 'latitude', 000, 000, 000 )
+!
+!--       Define coordinate-reference system
+          CALL netcdf_create_crs( id_set_yz(av), 000 )
 
           IF ( land_surface )  THEN
@@ -4916,4 +5308,59 @@
                 DEALLOCATE( netcdf_data_2d )
              ENDIF
+!
+!--          Write lon and lat data
+             ALLOCATE( lat(1:ns,0:ny) )
+             ALLOCATE( lon(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
+                         lat(i,:) = -90.0_wp   ! section averaged along x
+                         lon(i,:) = -180.0_wp  ! section averaged along x
+                      ELSE
+                         eutm = init_model%origin_x            &
+                              + cos_ra * ( section(i,3) + shift_x ) * dx  &
+                              + sin_ra * ( j + shift_y ) * dy
+                         nutm = init_model%origin_y            &
+                              - sin_ra * ( section(i,3) + shift_x ) * dx  &
+                              + cos_ra * ( j + shift_y ) * dy
+
+                         CALL  convert_utm_to_geographic( crs_list,          &
+                                                          eutm, nutm,        &
+                                                          lon(i,j), lat(i,j) )
+                      ENDIF
+                   ENDDO
+                ENDDO
+
+                nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_lon_yz(av,k), &
+                                     lon, start = (/ 1, 1 /),       &
+                                     count = (/ ns, ny+1 /) )
+                CALL netcdf_handle_error( 'netcdf_define_header', 556 )
+
+                nc_stat = NF90_PUT_VAR( id_set_yz(av), id_var_lat_yz(av,k), &
+                                     lat, start = (/ 1, 1 /),       &
+                                     count = (/ ns, ny+1 /) )
+                CALL netcdf_handle_error( 'netcdf_define_header', 556 )
+             ENDDO
+
+             DEALLOCATE( lat )
+             DEALLOCATE( lon )
 
           ENDIF
@@ -6896,3 +7343,216 @@
  END SUBROUTINE netcdf_create_var
 
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Create a variable holding the coordinate-reference-system information.
+!------------------------------------------------------------------------------!
+ SUBROUTINE netcdf_create_crs( ncid, error_no )
+
+#if defined( __netcdf )
+    IMPLICIT NONE
+
+    INTEGER, INTENT(IN)  ::  error_no  !< error number
+    INTEGER, INTENT(IN)  ::  ncid      !< file id
+    INTEGER              ::  var_id    !< variable id
+
+!
+!-- Define variable
+    nc_stat = NF90_DEF_VAR( ncid, 'crs', NF90_INT, VARID = var_id )
+    CALL netcdf_handle_error( 'netcdf_create_crs', error_no )
+!
+!-- Set attributes
+    nc_stat = NF90_PUT_ATT( ncid, var_id, 'epsg_code', &
+                            coord_ref_sys%epsg_code )
+    CALL netcdf_handle_error( 'netcdf_create_crs', error_no )
+
+    nc_stat = NF90_PUT_ATT( ncid, var_id, 'false_easting', &
+                            coord_ref_sys%false_easting )
+    CALL netcdf_handle_error( 'netcdf_create_crs', error_no )
+
+    nc_stat = NF90_PUT_ATT( ncid, var_id, 'false_northing', &
+                            coord_ref_sys%false_northing )
+    CALL netcdf_handle_error( 'netcdf_create_crs', error_no )
+
+    nc_stat = NF90_PUT_ATT( ncid, var_id, 'grid_mapping_name', &
+                            coord_ref_sys%grid_mapping_name )
+    CALL netcdf_handle_error( 'netcdf_create_crs', error_no )
+
+    nc_stat = NF90_PUT_ATT( ncid, var_id, 'inverse_flattening', &
+                            coord_ref_sys%inverse_flattening )
+    CALL netcdf_handle_error( 'netcdf_create_crs', error_no )
+
+    nc_stat = NF90_PUT_ATT( ncid, var_id, 'latitude_of_projection_origin', &
+                            coord_ref_sys%latitude_of_projection_origin )
+    CALL netcdf_handle_error( 'netcdf_create_crs', error_no )
+
+    nc_stat = NF90_PUT_ATT( ncid, var_id, 'long_name', &
+                            coord_ref_sys%long_name )
+    CALL netcdf_handle_error( 'netcdf_create_crs', error_no )
+
+    nc_stat = NF90_PUT_ATT( ncid, var_id, 'longitude_of_central_meridian', &
+                            coord_ref_sys%longitude_of_central_meridian )
+    CALL netcdf_handle_error( 'netcdf_create_crs', error_no )
+
+    nc_stat = NF90_PUT_ATT( ncid, var_id, 'longitude_of_prime_meridian', &
+                            coord_ref_sys%longitude_of_prime_meridian )
+    CALL netcdf_handle_error( 'netcdf_create_crs', error_no )
+
+    nc_stat = NF90_PUT_ATT( ncid, var_id, 'scale_factor_at_central_meridian', &
+                            coord_ref_sys%scale_factor_at_central_meridian )
+    CALL netcdf_handle_error( 'netcdf_create_crs', error_no )
+
+    nc_stat = NF90_PUT_ATT( ncid, var_id, 'semi_major_axis', &
+                            coord_ref_sys%semi_major_axis )
+    CALL netcdf_handle_error( 'netcdf_create_crs', error_no )
+
+    nc_stat = NF90_PUT_ATT( ncid, var_id, 'units', &
+                            coord_ref_sys%units )
+    CALL netcdf_handle_error( 'netcdf_create_crs', error_no )
+
+#endif
+END SUBROUTINE netcdf_create_crs
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Convert UTM coordinates into geographic latitude and longitude. Conversion
+!> is based on the work of KrÃŒger (1912) DOI: 10.2312/GFZ.b103-krueger28
+!> and Karney (2013) DOI: 10.1007/s00190-012-0578-z
+!> Based on a JavaScript of the geodesy function library written by chrisveness
+!> https://github.com/chrisveness/geodesy
+!------------------------------------------------------------------------------!
+ SUBROUTINE convert_utm_to_geographic( crs, eutm, nutm, lon, lat )
+
+    USE basic_constants_and_equations_mod,                                     &
+        ONLY:  pi
+
+    IMPLICIT NONE
+
+    REAL(wp), INTENT(in)  ::  eutm !< easting (UTM)
+    REAL(wp), INTENT(out) ::  lat  !< geographic latitude in degree
+    REAL(wp), INTENT(out) ::  lon  !< geographic longitude in degree
+    REAL(wp), INTENT(in)  ::  nutm !< northing (UTM)
+
+    REAL(wp) ::  a           !< 2*pi*a is the circumference of a meridian
+    REAL(wp) ::  cos_eta_s   !< cos(eta_s)
+    REAL(wp) ::  delta_i     !<
+    REAL(wp) ::  delta_tau_i !<
+    REAL(wp) ::  e           !< eccentricity
+    REAL(wp) ::  eta         !<
+    REAL(wp) ::  eta_s       !<
+    REAL(wp) ::  j           !< loop index
+    REAL(wp) ::  n           !< 3rd flattening
+    REAL(wp) ::  n2          !< n^2
+    REAL(wp) ::  n3          !< n^3
+    REAL(wp) ::  n4          !< n^4
+    REAL(wp) ::  n5          !< n^5
+    REAL(wp) ::  n6          !< n^6
+    REAL(wp) ::  nu          !<
+    REAL(wp) ::  nu_s        !<
+    REAL(wp) ::  sin_eta_s   !< sin(eta_s)
+    REAL(wp) ::  sinh_nu_s   !< sinush(nu_s)
+    REAL(wp) ::  tau_i       !<
+    REAL(wp) ::  tau_i_s     !<
+    REAL(wp) ::  tau_s       !<
+    REAL(wp) ::  x           !< adjusted easting
+    REAL(wp) ::  y           !< adjusted northing
+
+    REAL(wp), DIMENSION(6) ::  beta !< 6th order KrÃŒger expressions
+
+    REAL(wp), DIMENSION(8), INTENT(in) ::  crs !< coordinate reference system, consists of
+                                               !< (/semi_major_axis,
+                                               !<   inverse_flattening,
+                                               !<   longitude_of_prime_meridian,
+                                               !<   longitude_of_central_meridian,
+                                               !<   scale_factor_at_central_meridian,
+                                               !<   latitude_of_projection_origin,
+                                               !<   false_easting,
+                                               !<   false_northing /)
+
+    x = eutm - crs(7)  ! remove false easting
+    y = nutm - crs(8)  ! remove false northing
+
+!-- from Karney 2011 Eq 15-22, 36:
+    e = SQRT( 1.0_wp / crs(2) * ( 2.0_wp - 1.0_wp / crs(2) ) )
+    n = 1.0_wp / crs(2) / ( 2.0_wp - 1.0_wp / crs(2) )
+    n2 = n * n
+    n3 = n * n2
+    n4 = n * n3
+    n5 = n * n4
+    n6 = n * n5
+
+    a = crs(1) / ( 1.0_wp + n ) * ( 1.0_wp + 0.25_wp * n2       &
+                                           + 0.015625_wp * n4   &
+                                           + 3.90625E-3_wp * n6 )
+
+    nu  = x / ( crs(5) * a )
+    eta = y / ( crs(5) * a )
+
+!-- According to KrÃŒger (1912), eq. 26*
+    beta(1) =        0.5_wp                  * n  &
+            -        2.0_wp /         3.0_wp * n2 &
+            +       37.0_wp /        96.0_wp * n3 &
+            -        1.0_wp /       360.0_wp * n4 &
+            -       81.0_wp /       512.0_wp * n5 &
+            +    96199.0_wp /    604800.0_wp * n6
+
+    beta(2) =        1.0_wp /        48.0_wp * n2 &
+            +        1.0_wp /        15.0_wp * n3 &
+            -      437.0_wp /      1440.0_wp * n4 &
+            +       46.0_wp /       105.0_wp * n5 &
+            -  1118711.0_wp /   3870720.0_wp * n6
+
+    beta(3) =       17.0_wp /       480.0_wp * n3 &
+            -       37.0_wp /       840.0_wp * n4 &
+            -      209.0_wp /      4480.0_wp * n5 &
+            +     5569.0_wp /     90720.0_wp * n6
+
+    beta(4) =     4397.0_wp /    161280.0_wp * n4 &
+            -       11.0_wp /       504.0_wp * n5 &
+            -   830251.0_wp /   7257600.0_wp * n6
+
+    beta(5) =     4583.0_wp /    161280.0_wp * n5 &
+            -   108847.0_wp /   3991680.0_wp * n6
+
+    beta(6) = 20648693.0_wp / 638668800.0_wp * n6
+
+    eta_s = eta
+    nu_s  = nu
+    DO  j = 1, 6
+      eta_s = eta_s - beta(j) * SIN(2.0_wp * j * eta) * COSH(2.0_wp * j * nu)
+      nu_s  = nu_s  - beta(j) * COS(2.0_wp * j * eta) * SINH(2.0_wp * j * nu)
+    ENDDO
+
+    sinh_nu_s = SINH( nu_s )
+    sin_eta_s = SIN( eta_s )
+    cos_eta_s = COS( eta_s )
+
+    tau_s = sin_eta_s / SQRT( sinh_nu_s**2 + cos_eta_s**2 )
+
+    tau_i = tau_s
+    delta_tau_i = 1.0_wp
+
+    DO WHILE ( ABS( delta_tau_i ) > 1.0E-12_wp )
+
+      delta_i = SINH( e * ATANH( e * tau_i / SQRT( 1.0_wp + tau_i**2 ) ) )
+
+      tau_i_s = tau_i   * SQRT( 1.0_wp + delta_i**2 )  &
+               - delta_i * SQRT( 1.0_wp + tau_i**2 )
+
+      delta_tau_i = ( tau_s - tau_i_s ) / SQRT( 1.0_wp + tau_i_s**2 )  &
+                   * ( 1.0_wp + ( 1.0_wp - e**2 ) * tau_i**2 )          &
+                   / ( ( 1.0_wp - e**2 ) * SQRT( 1.0_wp + tau_i**2 ) )
+
+      tau_i = tau_i + delta_tau_i
+
+    ENDDO
+
+    lat = ATAN( tau_i ) / pi * 180.0_wp
+    lon = ATAN2( sinh_nu_s, cos_eta_s ) / pi * 180.0_wp + crs(4)
+
+ END SUBROUTINE convert_utm_to_geographic
+
  END MODULE netcdf_interface
Index: /palm/trunk/TESTS/cases/example_cbl/MONITORING/example_cbl_rc
===================================================================
--- /palm/trunk/TESTS/cases/example_cbl/MONITORING/example_cbl_rc	(revision 3458)
+++ /palm/trunk/TESTS/cases/example_cbl/MONITORING/example_cbl_rc	(revision 3459)
@@ -130,5 +130,5 @@
 
 
-       Profile: wpt, w"pt", w*pt*, w*2, pt, pt*2,
+       Profile: wtheta, w"theta", w*theta*, w*2, theta, theta*2,
        Output every               900.00 s
        Time averaged over         600.00 s
@@ -139,5 +139,5 @@
 
 
-       XY-cross-section  Arrays:  w, pt,
+       XY-cross-section  Arrays:  w, theta,
        Output every               900.00 s
        Cross sections at k = /2/10/
@@ -148,5 +148,5 @@
 
 
-       XZ-cross-section  Arrays:  w, pt,
+       XZ-cross-section  Arrays:  w, theta,
        Output every               900.00 s
        Cross sections at j = /20/
Index: /palm/trunk/TESTS/cases/example_cbl_short/MONITORING/example_cbl_short_rc
===================================================================
--- /palm/trunk/TESTS/cases/example_cbl_short/MONITORING/example_cbl_short_rc	(revision 3458)
+++ /palm/trunk/TESTS/cases/example_cbl_short/MONITORING/example_cbl_short_rc	(revision 3459)
@@ -131,5 +131,5 @@
 
 
-       Profile: wpt, w"pt", w*pt*, w*2, pt, pt*2,                            
+       Profile: wtheta, w"theta", w*theta*, w*2, theta, theta*2,                            
        Output every               900.00 s
        Time averaged over         600.00 s
@@ -140,5 +140,5 @@
 
 
-       XY-cross-section  Arrays:  w, pt,                                                               
+       XY-cross-section  Arrays:  w, theta,                                                               
        Output every               900.00 s  
        Cross sections at k = /2/10/
@@ -149,5 +149,5 @@
 
 
-       XZ-cross-section  Arrays:  w, pt,                                                               
+       XZ-cross-section  Arrays:  w, theta,                                                               
        Output every               900.00 s  
        Cross sections at j = /20/
Index: /palm/trunk/TESTS/cases/example_cbl_user_code/MONITORING/example_cbl_user_code_rc
===================================================================
--- /palm/trunk/TESTS/cases/example_cbl_user_code/MONITORING/example_cbl_user_code_rc	(revision 3458)
+++ /palm/trunk/TESTS/cases/example_cbl_user_code/MONITORING/example_cbl_user_code_rc	(revision 3459)
@@ -130,5 +130,5 @@
 
 
-       Profile: wpt, w"pt", w*pt*, w*2, pt, pt*2,
+       Profile: wtheta, w"theta", w*theta*, w*2, theta, theta*2,
        Output every               900.00 s
        Time averaged over         600.00 s
@@ -139,5 +139,5 @@
 
 
-       XY-cross-section  Arrays:  w, pt,
+       XY-cross-section  Arrays:  w, theta,
        Output every               900.00 s
        Cross sections at k = /2/10/
@@ -148,5 +148,5 @@
 
 
-       XZ-cross-section  Arrays:  w, pt,
+       XZ-cross-section  Arrays:  w, theta,
        Output every               900.00 s
        Cross sections at j = /20/
Index: /palm/trunk/TESTS/cases/test_oceanml/MONITORING/test_oceanml_rc
===================================================================
--- /palm/trunk/TESTS/cases/test_oceanml/MONITORING/test_oceanml_rc	(revision 3458)
+++ /palm/trunk/TESTS/cases/test_oceanml/MONITORING/test_oceanml_rc	(revision 3459)
@@ -130,5 +130,5 @@
 
 
-       Profile: e, e*, pt,
+       Profile: e, e*, theta,
        Output every             ******** s
        Time averaged over           0.00 s
