Index: palm/trunk/SCRIPTS/.palm.iofiles
===================================================================
--- palm/trunk/SCRIPTS/.palm.iofiles	(revision 4464)
+++ palm/trunk/SCRIPTS/.palm.iofiles	(revision 4465)
@@ -86,4 +86,3 @@
 PARTICLE_DATA*             out:lnpe *         $base_data/$run_identifier/OUTPUT         _prt_dat
 #
-#WTM_OUTPUT_DATA*           out:tr   *         $base_data/$run_identifier/MONITORING     _wtm
-DATA_1D_TS_WTM_NETCDF*     out:tr   *         $base_data/$run_identifier/OUTPUT   _wtm       nc
+DATA_1D_TS_WTM_NETCDF*     out:tr   *         $base_data/$run_identifier/OUTPUT          _wtm       nc
Index: palm/trunk/SOURCE/wind_turbine_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/wind_turbine_model_mod.f90	(revision 4464)
+++ palm/trunk/SOURCE/wind_turbine_model_mod.f90	(revision 4465)
@@ -26,4 +26,9 @@
 ! -----------------
 ! $Id$
+! Removed old ASCII outputm, some syntax layout adjustments, added output for
+! rotor and tower diameters. Added warning message in case of NetCDF 3 (no
+! WTM output file will be produced).
+! 
+! 4460 2020-03-12 16:47:30Z oliver.maas
 ! allow for simulating up to 10 000 wind turbines
 ! 
@@ -258,12 +263,4 @@
     REAL(wp) ::  yaw_misalignment_min = 0.008726_wp    !< minimum yaw missalignment for which the actuator stops [rad]
 
-!
-!-- Set flag for output files TURBINE_PARAMETERS
-    TYPE file_status
-       LOGICAL ::  opened, opened_before
-    END TYPE file_status
-    
-    TYPE(file_status), DIMENSION(500) :: openfile_turb_mod =                   &
-                                         file_status(.FALSE.,.FALSE.)
 
 !
@@ -542,5 +539,5 @@
 !--    Set flag that indicates that the wind turbine model is switched on
        wind_turbine = .TRUE.
-       
+
        GOTO 12
 
@@ -1391,7 +1388,13 @@
 !
 !-- Create NetCDF output file
-    nc_filename = 'DATA_1D_TS_WTM_NETCDF' // TRIM( coupling_char )      
+#if defined( __netcdf4 )
+    nc_filename = 'DATA_1D_TS_WTM_NETCDF' // TRIM( coupling_char )
     return_value = dom_def_file( nc_filename, 'netcdf4-serial' )
-
+#else
+    message_string = 'Wind turbine model output requires NetCDF version 4. ' //&
+                     'No output file will be created.'
+                   CALL message( 'wtm_init_output', 'PA0672',             &
+                                                                 0, 1, 0, 6, 0 )
+#endif
     IF ( myid == 0 )  THEN
 !
@@ -1408,6 +1411,4 @@
        DEALLOCATE( ndim )
 
-!
-!--    time
        return_value = dom_def_dim( nc_filename,                                &
                                    dimension_name = 'time',                    &
@@ -1416,15 +1417,12 @@
                                    values_realwp = (/0.0_wp/) )
  
-!
-!--    x
        variable_name = 'x'
-       return_value = dom_def_var( nc_filename,                        &
+       return_value = dom_def_var( nc_filename,                                &
                                    variable_name = variable_name,              &
                                    dimension_names = (/'turbine'/),            &
                                    output_type = 'real32' )
-!
-!--    y
+
        variable_name = 'y'
-       return_value = dom_def_var( nc_filename,                        &
+       return_value = dom_def_var( nc_filename,                                &
                                    variable_name = variable_name,              &
                                    dimension_names = (/'turbine'/),            &
@@ -1432,185 +1430,214 @@
 
        variable_name = 'z'
-       return_value = dom_def_var( nc_filename,                        &
+       return_value = dom_def_var( nc_filename,                                &
                                    variable_name = variable_name,              &
                                    dimension_names = (/'turbine'/),            &
                                    output_type = 'real32' )
-                                   
-
-       return_value = dom_def_att( nc_filename,                        &
+
+
+       variable_name = 'rotor_diameter'
+       return_value = dom_def_var( nc_filename,                                &
+                                   variable_name = variable_name,              &
+                                   dimension_names = (/'turbine'/),            &
+                                   output_type = 'real32' )
+
+       variable_name = 'tower_diameter'
+       return_value = dom_def_var( nc_filename,                                &
+                                   variable_name = variable_name,              &
+                                   dimension_names = (/'turbine'/),            &
+                                   output_type = 'real32' )
+
+       return_value = dom_def_att( nc_filename,                                &
                                    variable_name = 'time',                     &
-                                   attribute_name = 'units',                 &
+                                   attribute_name = 'units',                   &
                                    value = 'seconds since ' // origin_date_time )
-                                   
-       return_value = dom_def_att( nc_filename,                        &
+
+       return_value = dom_def_att( nc_filename,                                &
                                    variable_name = 'x',                        &
-                                   attribute_name = 'units',                 &
+                                   attribute_name = 'units',                   &
                                    value = 'm' )
 
-       return_value = dom_def_att( nc_filename,                        &
+       return_value = dom_def_att( nc_filename,                                &
                                    variable_name = 'y',                        &
-                                   attribute_name = 'units',                 &
+                                   attribute_name = 'units',                   &
                                    value = 'm' )      
 
-       return_value = dom_def_att( nc_filename,                        &
+       return_value = dom_def_att( nc_filename,                                &
                                    variable_name = 'z',                        &
-                                   attribute_name = 'units',                 &
-                                   value = 'm' )                                   
-
-       return_value = dom_def_att( nc_filename,                        &
+                                   attribute_name = 'units',                   &
+                                   value = 'm' )
+
+       return_value = dom_def_att( nc_filename,                                &
+                                   variable_name = 'rotor_diameter',           &
+                                   attribute_name = 'units',                   &
+                                   value = 'm' )
+
+       return_value = dom_def_att( nc_filename,                                &
+                                   variable_name = 'tower_diameter',           &
+                                   attribute_name = 'units',                   &
+                                   value = 'm' )
+
+       return_value = dom_def_att( nc_filename,                                &
                                    variable_name = 'x',                        &
-                                   attribute_name = 'long_name',                 &
+                                   attribute_name = 'long_name',               &
                                    value = 'x location of rotor center' )
 
-       return_value = dom_def_att( nc_filename,                        &
+       return_value = dom_def_att( nc_filename,                                &
                                    variable_name = 'y',                        &
-                                   attribute_name = 'long_name',                 &
+                                   attribute_name = 'long_name',               &
                                    value = 'y location of rotor center' )      
 
-       return_value = dom_def_att( nc_filename,                        &
+       return_value = dom_def_att( nc_filename,                                &
                                    variable_name = 'z',                        &
-                                   attribute_name = 'long_name',                 &
+                                   attribute_name = 'long_name',               &
                                    value = 'z location of rotor center' )     
-                                   
-                                   
-       return_value = dom_def_att( nc_filename,                        &
+
+       return_value = dom_def_att( nc_filename,                                &
                                    variable_name = 'turbine_name',             &
-                                   attribute_name = 'long_name',                 &
-                                   value = 'turbine name')   
- 
-       return_value = dom_def_att( nc_filename,                        &
+                                   attribute_name = 'long_name',               &
+                                   value = 'turbine name')
+
+       return_value = dom_def_att( nc_filename,                                &
+                                   variable_name = 'rotor_diameter',           &
+                                   attribute_name = 'long_name',               &
+                                   value = 'rotor diameter')
+
+       return_value = dom_def_att( nc_filename,                                &
+                                   variable_name = 'tower_diameter',           &
+                                   attribute_name = 'long_name',               &
+                                   value = 'tower diameter')
+
+       return_value = dom_def_att( nc_filename,                               &
                                    variable_name = 'time',                     &
-                                   attribute_name = 'standard_name',             &
+                                   attribute_name = 'standard_name',           &
                                    value = 'time')
 
-       return_value = dom_def_att( nc_filename,                        &
+       return_value = dom_def_att( nc_filename,                                &
                                    variable_name = 'time',                     &
                                    attribute_name = 'axis',                    &
                                    value = 'T')
 
-       return_value = dom_def_att( nc_filename,                        &
+       return_value = dom_def_att( nc_filename,                                &
                                    variable_name = 'x',                        &
                                    attribute_name = 'axis',                    &
                                    value = 'X' )
 
-       return_value = dom_def_att( nc_filename,                        &
+       return_value = dom_def_att( nc_filename,                                &
                                    variable_name = 'y',                        &
                                    attribute_name = 'axis',                    &
-                                   value = 'Y' )                                  
-
-
-       variable_name = 'rotor_speed'
-       return_value = dom_def_var( nc_filename,                        &
-                                   variable_name = variable_name,      &
-                                   dimension_names = (/ 'turbine', 'time   ' /),     &
+                                   value = 'Y' )
+
+                                   
+       variable_name = 'generator_power'
+       return_value = dom_def_var( nc_filename,                                  &
+                                   variable_name = variable_name,                &
+                                   dimension_names = (/ 'turbine', 'time   ' /), &
+                                   output_type = 'real32' )
+
+       return_value = dom_def_att( nc_filename,                                &
+                                   variable_name = variable_name,              &
+                                   attribute_name = 'units',                   &
+                                   value = 'W' )  
+
+       variable_name = 'generator_speed'
+       return_value = dom_def_var( nc_filename,                                  &
+                                   variable_name = variable_name,                &
+                                   dimension_names = (/ 'turbine', 'time   ' /), &
+                                   output_type = 'real32' )
+
+       return_value = dom_def_att( nc_filename,                                &
+                                   variable_name = variable_name,              &
+                                   attribute_name = 'units',                   &
+                                   value = 'rad/s' )
+
+       variable_name = 'generator_torque'
+       return_value = dom_def_var( nc_filename,                                  &
+                                   variable_name = variable_name,                &
+                                   dimension_names = (/ 'turbine', 'time   ' /), &
                                    output_type = 'real32' )
   
-       return_value = dom_def_att( nc_filename,                  &
-                                   variable_name = variable_name,        &
-                                   attribute_name = 'units',           &
-                                   value = 'rad/s' )
-  
-       variable_name = 'generator_speed'
-       return_value = dom_def_var( nc_filename,                        &
-                                   variable_name = variable_name,      &
-                                   dimension_names = (/ 'turbine', 'time   ' /),     &
+       return_value = dom_def_att( nc_filename,                                &
+                                   variable_name = variable_name,              &
+                                   attribute_name = 'units',                   &
+                                   value = 'Nm' )
+
+       variable_name = 'pitch_angle'
+       return_value = dom_def_var( nc_filename,                                  &
+                                   variable_name = variable_name,                &
+                                   dimension_names = (/ 'turbine', 'time   ' /), &
+                                   output_type = 'real32' )
+
+       return_value = dom_def_att( nc_filename,                                &
+                                   variable_name = variable_name,              &
+                                   attribute_name = 'units',                   &
+                                   value = 'degrees' )
+
+       variable_name = 'rotor_power'
+       return_value = dom_def_var( nc_filename,                                  &
+                                   variable_name = variable_name,                &
+                                   dimension_names = (/ 'turbine', 'time   ' /), &
                                    output_type = 'real32' )
   
-       return_value = dom_def_att( nc_filename,                  &
-                                   variable_name = variable_name,        &
-                                   attribute_name = 'units',           &
-                                   value = 'rad/s' )
-      
-      
-       variable_name = 'generator_torque'
-       return_value = dom_def_var( nc_filename,                        &
-                                   variable_name = variable_name,      &
-                                   dimension_names = (/ 'turbine', 'time   ' /),     &
+       return_value = dom_def_att( nc_filename,                                &
+                                   variable_name = variable_name,              &
+                                   attribute_name = 'units',                   &
+                                   value = 'W' )  
+
+       variable_name = 'rotor_speed'
+       return_value = dom_def_var( nc_filename,                                  &
+                                   variable_name = variable_name,                &
+                                   dimension_names = (/ 'turbine', 'time   ' /), &
                                    output_type = 'real32' )
   
-       return_value = dom_def_att( nc_filename,                  &
-                                   variable_name = variable_name,        &
-                                   attribute_name = 'units',           &
-                                   value = 'Nm' )  
-        
-       variable_name = 'rotor_torque'
-       return_value = dom_def_var( nc_filename,                        &
-                                   variable_name = variable_name,      &
-                                   dimension_names = (/ 'turbine', 'time   ' /),     &
+       return_value = dom_def_att( nc_filename,                                &
+                                   variable_name = variable_name,              &
+                                   attribute_name = 'units',                   &
+                                   value = 'rad/s' )
+
+       variable_name = 'rotor_thrust'
+       return_value = dom_def_var( nc_filename,                                  &
+                                   variable_name = variable_name,                &
+                                   dimension_names = (/ 'turbine', 'time   ' /), &
                                    output_type = 'real32' )
   
-       return_value = dom_def_att( nc_filename,                  &
-                                   variable_name = variable_name,        &
-                                   attribute_name = 'units',           &
-                                   value = 'Nm' )
-
-       variable_name = 'pitch_angle'
-       return_value = dom_def_var( nc_filename,                        &
-                                   variable_name = variable_name,      &
-                                   dimension_names = (/ 'turbine', 'time   ' /),     &
+       return_value = dom_def_att( nc_filename,                                &
+                                   variable_name = variable_name,              &
+                                   attribute_name = 'units',                   &
+                                   value = 'N' )    
+
+       variable_name = 'rotor_torque'
+       return_value = dom_def_var( nc_filename,                                  &
+                                   variable_name = variable_name,                &
+                                   dimension_names = (/ 'turbine', 'time   ' /), &
                                    output_type = 'real32' )
   
-       return_value = dom_def_att( nc_filename,                  &
-                                   variable_name = variable_name,        &
-                                   attribute_name = 'units',           &
-                                   value = 'degrees' )     
-                                   
-       variable_name = 'generator_power'
-       return_value = dom_def_var( nc_filename,                        &
-                                   variable_name = variable_name,      &
-                                   dimension_names = (/ 'turbine', 'time   ' /),     &
+       return_value = dom_def_att( nc_filename,                                &
+                                   variable_name = variable_name,              &
+                                   attribute_name = 'units',                   &
+                                   value = 'Nm' )
+
+       variable_name = 'wind_direction'
+       return_value = dom_def_var( nc_filename,                                  &
+                                   variable_name = variable_name,                &
+                                   dimension_names = (/ 'turbine', 'time   ' /), &
                                    output_type = 'real32' )
   
-       return_value = dom_def_att( nc_filename,                  &
-                                   variable_name = variable_name,        &
-                                   attribute_name = 'units',           &
-                                   value = 'W' )  
-                                   
-       variable_name = 'rotor_power'
-       return_value = dom_def_var( nc_filename,                        &
-                                   variable_name = variable_name,      &
-                                   dimension_names = (/ 'turbine', 'time   ' /),     &
+       return_value = dom_def_att( nc_filename,                                &
+                                   variable_name = variable_name,              &
+                                   attribute_name = 'units',                   &
+                                   value = 'degrees' )  
+
+       variable_name = 'yaw_angle'
+       return_value = dom_def_var( nc_filename,                                  &
+                                   variable_name = variable_name,                &
+                                   dimension_names = (/ 'turbine', 'time   ' /), &
                                    output_type = 'real32' )
   
-       return_value = dom_def_att( nc_filename,                  &
-                                   variable_name = variable_name,        &
-                                   attribute_name = 'units',           &
-                                   value = 'W' )  
-                                   
-        variable_name = 'rotor_thrust'
-       return_value = dom_def_var( nc_filename,                        &
-                                   variable_name = variable_name,      &
-                                   dimension_names = (/ 'turbine', 'time   ' /),     &
-                                   output_type = 'real32' )
-  
-       return_value = dom_def_att( nc_filename,                  &
-                                   variable_name = variable_name,        &
-                                   attribute_name = 'units',           &
-                                   value = 'N' )    
-                                   
-                                   
-       variable_name = 'wind_direction'
-       return_value = dom_def_var( nc_filename,                        &
-                                   variable_name = variable_name,      &
-                                   dimension_names = (/ 'turbine', 'time   ' /),     &
-                                   output_type = 'real32' )
-  
-       return_value = dom_def_att( nc_filename,                  &
-                                   variable_name = variable_name,        &
-                                   attribute_name = 'units',           &
+       return_value = dom_def_att( nc_filename,                                &
+                                   variable_name = variable_name,              &
+                                   attribute_name = 'units',                   &
                                    value = 'degrees' )  
-                                   
-       variable_name = 'yaw_angle'
-       return_value = dom_def_var( nc_filename,                        &
-                                   variable_name = variable_name,      &
-                                   dimension_names = (/ 'turbine', 'time   ' /),     &
-                                   output_type = 'real32' )
-  
-       return_value = dom_def_att( nc_filename,                  &
-                                   variable_name = variable_name,        &
-                                   attribute_name = 'units',           &
-                                   value = 'degrees' )  
-
-   ENDIF    
+
+   ENDIF
 END SUBROUTINE
     
@@ -1897,27 +1924,27 @@
 !--    the rotor. If tilt = 0, rot_coord_trans is a unit matrix.
 
-       rot_coord_trans(inot,1,1) = rot_eigen_rad(1)**2                   *     &
+       rot_coord_trans(inot,1,1) = rot_eigen_rad(1)**2                   *           &
                                    ( 1.0_wp - COS( tilt_angle ) ) + COS( tilt_angle ) 
-       rot_coord_trans(inot,1,2) = rot_eigen_rad(1) * rot_eigen_rad(2)   *     &
+       rot_coord_trans(inot,1,2) = rot_eigen_rad(1) * rot_eigen_rad(2)   *           &
                                    ( 1.0_wp - COS( tilt_angle ) )              -     &
                                    rot_eigen_rad(3) * SIN( tilt_angle )
-       rot_coord_trans(inot,1,3) = rot_eigen_rad(1) * rot_eigen_rad(3)   *     &
+       rot_coord_trans(inot,1,3) = rot_eigen_rad(1) * rot_eigen_rad(3)   *           &
                                    ( 1.0_wp - COS( tilt_angle ) )              +     &
                                    rot_eigen_rad(2) * SIN( tilt_angle )
-       rot_coord_trans(inot,2,1) = rot_eigen_rad(2) * rot_eigen_rad(1)   *     &
+       rot_coord_trans(inot,2,1) = rot_eigen_rad(2) * rot_eigen_rad(1)   *           &
                                    ( 1.0_wp - COS( tilt_angle ) )              +     &
                                    rot_eigen_rad(3) * SIN( tilt_angle )
-       rot_coord_trans(inot,2,2) = rot_eigen_rad(2)**2                   *     &
+       rot_coord_trans(inot,2,2) = rot_eigen_rad(2)**2                   *           &
                                    ( 1.0_wp - COS( tilt_angle ) ) + COS( tilt_angle ) 
-       rot_coord_trans(inot,2,3) = rot_eigen_rad(2) * rot_eigen_rad(3)   *     &
+       rot_coord_trans(inot,2,3) = rot_eigen_rad(2) * rot_eigen_rad(3)   *           &
                                    ( 1.0_wp - COS( tilt_angle ) )              -     &
                                    rot_eigen_rad(1) * SIN( tilt_angle )
-       rot_coord_trans(inot,3,1) = rot_eigen_rad(3) * rot_eigen_rad(1)   *     &
+       rot_coord_trans(inot,3,1) = rot_eigen_rad(3) * rot_eigen_rad(1)   *           &
                                    ( 1.0_wp - COS( tilt_angle ) )              -     &
                                    rot_eigen_rad(2) * SIN( tilt_angle )
-       rot_coord_trans(inot,3,2) = rot_eigen_rad(3) * rot_eigen_rad(2)   *     &
+       rot_coord_trans(inot,3,2) = rot_eigen_rad(3) * rot_eigen_rad(2)   *           &
                                    ( 1.0_wp - COS( tilt_angle ) )              +     &
                                    rot_eigen_rad(1) * SIN( tilt_angle )
-       rot_coord_trans(inot,3,3) = rot_eigen_rad(3)**2                   *     &
+       rot_coord_trans(inot,3,3) = rot_eigen_rad(3)**2                   *           &
                                    ( 1.0_wp - COS( tilt_angle ) ) + COS( tilt_angle )
 
@@ -2409,21 +2436,21 @@
 
                    IF ( tip_loss_correction )  THEN
-                   
+!
 !--                  Tip loss correction following Schito
 !--                  Schito applies the tip loss correction only to the lift force
 !--                  Therefore, the tip loss correction is only applied to the lift 
 !--                  coefficient and not to the drag coefficient in our case 
-!--                  
+!--
                      IF ( phi_rel(rseg) == 0.0_wp ) THEN
                         tl_factor = 1.0_wp
                      ELSE
-                        tl_factor = ( 2.0 / pi ) *                                &
-                             ACOS( EXP( -1.0 * ( 3.0 * ( rotor_radius(inot) - cur_r ) /     &
+                        tl_factor = ( 2.0 / pi ) *                                      &
+                             ACOS( EXP( -1.0 * ( 3.0 * ( rotor_radius(inot) - cur_r ) / &
                              ( 2.0 * cur_r * abs( sin( phi_rel(rseg) ) ) ) ) ) )
                      ENDIF
-                     
-                     turb_cl(rseg)  = tl_factor * turb_cl(rseg)                                  
-                                 
-                   ENDIF               
+
+                     turb_cl(rseg)  = tl_factor * turb_cl(rseg)
+
+                   ENDIF
 !
 !--                !-----------------------------------------------------!
@@ -2464,6 +2491,6 @@
                    !$OMP CRITICAL
                    thrust_rotor(inot) = thrust_rotor(inot) +                   &
-                                        thrust_seg(rseg)                   
-                   
+                                        thrust_seg(rseg)
+
 
                    torque_total(inot) = torque_total(inot) + (torque_seg * cur_r)
@@ -2486,10 +2513,10 @@
              CALL cpu_log( log_point_s(62), 'wtm_controller', 'start' )
 
-             
+
              IF ( speed_control )  THEN
 !
 !--             Calculation of the current generator speed for rotor speed control
-             
-!                                     
+
+! 
 !--             The acceleration of the rotor speed is calculated from 
 !--             the force balance of the accelerating torque
@@ -2504,7 +2531,7 @@
 !--             and rotor speed
                 generator_speed(inot) = gear_ratio * ( rotor_speed(inot) + om_rate )     
-             
+
              ENDIF
-             
+
              IF ( pitch_control )  THEN
 
@@ -2523,5 +2550,5 @@
                    CYCLE pit_loop
                 ENDIF
-                
+
 !
 !--             The current pitch is saved for the next time step
@@ -2535,5 +2562,5 @@
 !
 !--          Call the rotor speed controller
-             
+
              IF ( speed_control )  THEN
 !
@@ -2546,5 +2573,5 @@
                    ENDIF
                 ENDIF
-                               
+
              ENDIF
 
@@ -2592,5 +2619,5 @@
 !--                         Predetermine flag to mask topography
                             flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
-      
+
 !
 !--                         Determine the square of the distance between the
@@ -2636,5 +2663,5 @@
                          ENDDO  ! End of loop over rseg
                       ENDDO     ! End of loop over ring
-              
+
 !
 !--                   Rotation of force components:
@@ -2644,5 +2671,5 @@
                                       torque_z(k,j,i)*rotz(inot,1)             &
                                                               ) * flag
-                                
+
                       rot_tend_y(k,j,i) = rot_tend_y(k,j,i) + (                &
                                       thrust(k,j,i)*rotx(inot,2) +             &
@@ -2650,5 +2677,5 @@
                                       torque_z(k,j,i)*rotz(inot,2)             &
                                                               ) * flag
-                                
+
                       rot_tend_z(k,j,i) = rot_tend_z(k,j,i) + (                &
                                       thrust(k,j,i)*rotx(inot,3) +             &
@@ -2663,8 +2690,8 @@
 
              CALL cpu_log( log_point_s(63), 'wtm_smearing', 'stop' )          
-                   
+
           ENDDO                  !-- end of loop over turbines
 
-                
+
           IF ( yaw_control )  THEN
 !
@@ -2676,5 +2703,5 @@
                 wd30 = 999.0_wp                  ! Set to dummy value
                 ALLOCATE( wd30_l(1:WDLON) )
-                
+
                 WDSHO = MAX( 1 , NINT( 2.0_wp / dt_3d ) )   ! 2s running mean array
                 ALLOCATE( wd2(1:n_turbines,1:WDSHO) )
@@ -2682,5 +2709,5 @@
                 ALLOCATE( wd2_l(1:WDSHO) )
                 start_up = .FALSE.
-             ENDIF          
+             ENDIF
 
 !
@@ -2707,8 +2734,8 @@
 
                       yaw_angle_l(inot) = yaw_angle(inot)
-                      
+
                    ENDIF
                 ENDIF
-                   
+
              ENDDO                                 !-- end of loop over turbines
 
@@ -2726,16 +2753,15 @@
              wdir     = wdir_l
              yaw_angle  = yaw_angle_l
-             
-             
+
 #endif
              DO inot = 1, n_turbines
-!             
-!--             Update rotor orientation                
+!
+!--             Update rotor orientation
                 CALL wtm_rotate_rotor( inot )
 
              ENDDO ! End of loop over turbines
-                           
+
           ENDIF  ! end of yaw control 
-          
+
           IF ( speed_control )  THEN
 !
@@ -2755,70 +2781,15 @@
              generator_speed_f_old = generator_speed_f
 #endif
-            
+
           ENDIF
 
-          
-          
-          
-          DO inot = 1, n_turbines
-
-
-          
-             IF ( myid == 0 ) THEN
-                IF ( openfile_turb_mod(400+inot)%opened )  THEN
-                   WRITE ( 400+inot, 106 ) time_since_reference_point, rotor_speed(inot),          &
-                             generator_speed(inot), torque_gen_old(inot),                          &
-                             torque_total(inot), pitch_angle(inot),                                &
-                             torque_gen_old(inot)*generator_speed(inot)*generator_efficiency,      &
-                             torque_total(inot)*rotor_speed(inot)*air_density,                     &
-                             thrust_rotor(inot),                                                   & 
-                             wdir(inot)*180.0_wp/pi,                                               &
-                             (yaw_angle(inot))*180.0_wp/pi                   
-                             
-                ELSE
-
-                   WRITE ( turbine_id,'(A2,I2.2)')  '_T', inot
-                   OPEN ( 400+inot, FILE=( 'WTM_OUTPUT_DATA' //                &
-                                           TRIM( coupling_char ) //            &
-                                           turbine_id ), FORM='FORMATTED' )
-                   WRITE ( 400+inot, 105 ) inot
-                   WRITE ( 400+inot, 106 ) time_since_reference_point, rotor_speed(inot),          &
-                             generator_speed(inot), torque_gen_old(inot),                          &
-                             torque_total(inot), pitch_angle(inot),                                &
-                             torque_gen_old(inot)*generator_speed(inot)*generator_efficiency,      &
-                             torque_total(inot)*rotor_speed(inot)*air_density,                     &
-                             thrust_rotor(inot),                                                   &
-                             wdir(inot)*180.0_wp/pi,                                               &
-                             (yaw_angle(inot))*180.0_wp/pi
-                ENDIF
-             ENDIF
-
-!--          Set open flag
-             openfile_turb_mod(400+inot)%opened = .TRUE.
-          ENDDO                                    !-- end of loop over turbines
-
-          
-          
-       ENDIF   
-
-       
+       ENDIF
+
        CALL cpu_log( log_point_s(61), 'wtm_forces', 'stop' )
-       
-!
-!--    Formats
-       105 FORMAT ('Turbine control data for turbine ',I2,1X,':'/ &
-              &'----------------------------------------'/ &
-              &'     Time   RSpeed   GSpeed     GenTorque     AeroTorque  ',   &
-               'Pitch    Power(Gen)    Power(Rot)    RotThrust WDirection  ',  &
-               'YawOrient')
-
-       106 FORMAT (F9.2,2X,F7.3,2X,F7.2,2X,F12.1,3X,F12.1,1X,F6.2,2X,F12.1,2X, &
-                   F12.1,1X,F12.1,4X,F7.2,4X,F7.2)
-
-           
-                   
+
+
     END SUBROUTINE wtm_forces
 
-    
+
 !------------------------------------------------------------------------------!
 ! Description:
@@ -2827,9 +2798,9 @@
 !------------------------------------------------------------------------------!
     SUBROUTINE wtm_yawcontrol( inot )
-    
+
        USE kinds
-                
+
        IMPLICIT NONE
-      
+
        INTEGER(iwp)             :: inot
        INTEGER(iwp)             :: i_wd_30
@@ -2876,8 +2847,8 @@
              ENDIF
           ENDIF
-         
+
           wd30(inot,:) = wd30_l
 
-!      
+!
 !--    If turbine is already yawing:
 !--    Initialize 2 s running mean and yaw until the missalignment is smaller
@@ -2890,5 +2861,5 @@
           wd2_l = CSHIFT( wd2_l, SHIFT = -1 )
           wd2_l(1) = wdir(inot)
-!      
+!
 !--       Check if array is full ( no more dummies )
           IF ( .NOT. ANY( wd2_l == 999.0_wp ) ) THEN
@@ -2915,9 +2886,9 @@
              yaw_angle(inot) = yaw_angle(inot) + yawdir(inot) * yaw_speed * dt_3d
           ENDIF
-      
+
           wd2(inot,:) = wd2_l
-            
+
        ENDIF
-     
+
     END SUBROUTINE wtm_yawcontrol 
 
@@ -2976,6 +2947,5 @@
 
        INTEGER(iwp)             :: inot
-       
-         
+
 
 !
@@ -2992,13 +2962,13 @@
 
        IF ( generator_speed_f(inot) <= region_15_min )  THEN
-!                        
+!
 !--       Region 1: Generator torque is set to zero to accelerate the rotor:
           torque_gen(inot) = 0
-       
+
        ELSEIF ( generator_speed_f(inot) <= region_2_min )  THEN
-!                        
+!
 !--       Region 1.5: Generator torque is increasing linearly with rotor speed:
           torque_gen(inot) = slope15 * ( generator_speed_f(inot) - region_15_min )
-                         
+
        ELSEIF ( generator_speed_f(inot) <= region_25_min )  THEN
 !
@@ -3006,31 +2976,31 @@
 !--                 speed to keep the TSR optimal:
           torque_gen(inot) = region_2_slope * generator_speed_f(inot) * generator_speed_f(inot)
-       
+
        ELSEIF ( generator_speed_f(inot) < generator_speed_rated )  THEN
-!                        
+!
 !--       Region 2.5: Slipage region between 2 and 3:
           torque_gen(inot) = region_25_slope * ( generator_speed_f(inot) - vs_sysp )
-       
+
        ELSE
-!                        
+!
 !--       Region 3: Generator torque is antiproportional to the rotor speed to
 !--                 keep the power constant:
           torque_gen(inot) = generator_power_rated / generator_speed_f(inot)
-       
+
        ENDIF
-!                        
+!
 !--    Calculate torque rate and confine with a max
        trq_rate = ( torque_gen(inot) - torque_gen_old(inot) ) / dt_3d
        trq_rate = MIN( MAX( trq_rate, -1.0_wp * generator_torque_rate_max ), &
                   generator_torque_rate_max )
-!                        
-!--    Calculate new gen torque and confine with max torque                         
+!
+!--    Calculate new gen torque and confine with max torque
        torque_gen(inot) = torque_gen_old(inot) + trq_rate * dt_3d
-       torque_gen(inot) = MIN( torque_gen(inot), generator_torque_max )                                             
-!
-!--    Overwrite values for next timestep                        
+       torque_gen(inot) = MIN( torque_gen(inot), generator_torque_max )
+!
+!--    Overwrite values for next timestep
        rotor_speed_l(inot) = generator_speed(inot) / gear_ratio
 
-    
+
     END SUBROUTINE wtm_speed_control    
 
@@ -3090,7 +3060,7 @@
                    DO  k = nzb+1, MAXVAL(k_hub) + MAXVAL(k_smear)
                       tend_nac_y = 0.5_wp * nac_cd_surf(k,j,i) *               &
-                                         SIGN( v(k,j,i)**2 , v(k,j,i) )     
+                                         SIGN( v(k,j,i)**2 , v(k,j,i) )
                       tend_tow_y   = 0.5_wp * tow_cd_surf(k,j,i) *             &
-                                         SIGN( v(k,j,i)**2 , v(k,j,i) )                      
+                                         SIGN( v(k,j,i)**2 , v(k,j,i) )
                       tend(k,j,i) = tend(k,j,i) + ( - rot_tend_y(k,j,i)        &
                                   - tend_nac_y - tend_tow_y )                  &
@@ -3156,7 +3126,7 @@
 !--             Calculate the thrust generated by the nacelle and the tower 
                 tend_nac_x = 0.5_wp * nac_cd_surf(k,j,i) *                     &
-                                   SIGN( u(k,j,i)**2 , u(k,j,i) )     
+                                   SIGN( u(k,j,i)**2 , u(k,j,i) )
                 tend_tow_x   = 0.5_wp * tow_cd_surf(k,j,i) *                   &
-                                   SIGN( u(k,j,i)**2 , u(k,j,i) ) 
+                                   SIGN( u(k,j,i)**2 , u(k,j,i) )
                 tend(k,j,i) = tend(k,j,i) + ( - rot_tend_x(k,j,i)              &
                             - tend_nac_x - tend_tow_x )                        &
@@ -3172,7 +3142,7 @@
              DO  k = nzb+1,  MAXVAL(k_hub) + MAXVAL(k_smear)
                 tend_nac_y = 0.5_wp * nac_cd_surf(k,j,i) *                     &
-                                   SIGN( v(k,j,i)**2 , v(k,j,i) )     
+                                   SIGN( v(k,j,i)**2 , v(k,j,i) )
                 tend_tow_y   = 0.5_wp * tow_cd_surf(k,j,i) *                   &
-                                   SIGN( v(k,j,i)**2 , v(k,j,i) )                      
+                                   SIGN( v(k,j,i)**2 , v(k,j,i) )
                 tend(k,j,i) = tend(k,j,i) + ( - rot_tend_y(k,j,i)              &
                             - tend_nac_y - tend_tow_y )                        &
@@ -3217,5 +3187,5 @@
              ALLOCATE ( output_values_1d_target(1:n_turbines) )
              output_values_1d_target = hub_x(1:n_turbines)
-             output_values_1d_pointer => output_values_1d_target      
+             output_values_1d_pointer => output_values_1d_target
              return_value = dom_write_var( nc_filename,                              &
                                         'x',                                         &
@@ -3225,5 +3195,5 @@
 
              output_values_1d_target = hub_y(1:n_turbines)
-             output_values_1d_pointer => output_values_1d_target      
+             output_values_1d_pointer => output_values_1d_target
              return_value = dom_write_var( nc_filename,                              &
                                         'y',                                         &
@@ -3233,15 +3203,31 @@
 
              output_values_1d_target = hub_z(1:n_turbines)
-             output_values_1d_pointer => output_values_1d_target      
+             output_values_1d_pointer => output_values_1d_target
              return_value = dom_write_var( nc_filename,                              &
                                         'z',                                         &
                                         values_realwp_1d = output_values_1d_pointer, &
                                         bounds_start = (/1/),                        &
-                                        bounds_end   = (/n_turbines/) )                                       
-                                        
+                                        bounds_end   = (/n_turbines/) )
+
+             output_values_1d_target = rotor_radius(1:n_turbines) * 2.0_wp
+             output_values_1d_pointer => output_values_1d_target
+             return_value = dom_write_var( nc_filename,                              &
+                                        'rotor_diameter',                            &
+                                        values_realwp_1d = output_values_1d_pointer, &
+                                        bounds_start = (/1/),                        &
+                                        bounds_end   = (/n_turbines/) )
+
+             output_values_1d_target = tower_diameter(1:n_turbines)
+             output_values_1d_pointer => output_values_1d_target
+             return_value = dom_write_var( nc_filename,                              &
+                                        'tower_diameter',                            &
+                                        values_realwp_1d = output_values_1d_pointer, &
+                                        bounds_start = (/1/),                        &
+                                        bounds_end   = (/n_turbines/) )
+
              initial_write_coordinates = .TRUE.
              DEALLOCATE ( output_values_1d_target )
           ENDIF
-          
+
           t_ind = t_ind + 1
           
