source: palm/trunk/SOURCE/virtual_flight_mod.f90 @ 3274

!> @file virtual_flights_mod.f90
!------------------------------------------------------------------------------!
! This file is part of the PALM model system.
!
! PALM is free software: you can redistribute it and/or modify it under the
! terms of the GNU General Public License as published by the Free Software
! Foundation, either version 3 of the License, or (at your option) any later
! version.
!
! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License along with
! PALM. If not, see <http://www.gnu.org/licenses/>.
!
! Copyright 1997-2018 Leibniz Universitaet Hannover
!------------------------------------------------------------------------------!
!
! Current revisions:
! ------------------
! 
! 
! Former revisions:
! -----------------
! $Id: virtual_flight_mod.f90 3274 2018-09-24 15:42:55Z knoop $
! Modularization of all bulk cloud physics code components
! 
! 3248 2018-09-14 09:42:06Z sward
! Minor formating changes
! 
! 3246 2018-09-13 15:14:50Z sward
! Added error handling for input namelist via parin_fail_message
! 
! 3241 2018-09-12 15:02:00Z raasch
! unused variables removed
! 
! 3182 2018-07-27 13:36:03Z suehring
! IF statement revised to consider new vertical grid stretching
! 
! 3049 2018-05-29 13:52:36Z Giersch
! Error messages revised
! 
! 2932 2018-03-26 09:39:22Z Giersch
! renamed flight_par to virtual_flight_parameters
! 
! 2894 2018-03-15 09:17:58Z Giersch
! variable named found has been introduced for checking if restart data was 
! found, reading of restart strings has been moved completely to 
! read_restart_data_mod, virtual_flight_prerun flag has been removed, redundant
! skipping function has been removed, flight_read/write_restart_data 
! have been renamed to flight_r/wrd_global, flight_rrd_global is called in 
! read_restart_data_mod now, marker *** end flight *** was removed, strings and
! their respective lengths are written out and read now in case of restart runs
! to get rid of prescribed character lengths, CASE DEFAULT was added if restart
! data is read
! 
! 2872 2018-03-12 10:05:47Z Giersch
! virutal_flight_prerun flag is used to define if module 
! related parameters were outputted as restart data
! 
! 2718 2018-01-02 08:49:38Z maronga
! Corrected "Former revisions" section
! 
! 2696 2017-12-14 17:12:51Z kanani
! Change in file header (GPL part)
!
! 2576 2017-10-24 13:49:46Z Giersch
! Definition of a new function called flight_skip_var_list to skip module 
! parameters during reading restart data
! 
! 2563 2017-10-19 15:36:10Z Giersch
! flight_read_restart_data is called in flight_parin in the case of a restart 
! run. flight_skip_var_list is not necessary anymore due to marker changes in
! restart files.
! 
! 2271 2017-06-09 12:34:55Z sward
! Todo added
! 
! 2101 2017-01-05 16:42:31Z suehring
!
! 2000 2016-08-20 18:09:15Z knoop
! Forced header and separation lines into 80 columns
! 
! 1960 2016-07-12 16:34:24Z suehring
! Separate humidity and passive scalar.
! Bugfix concerning labeling of timeseries.
! 
! 1957 2016-07-07 10:43:48Z suehring
! Initial revision
!
! Description:
! ------------
!> Module for virtual flight measurements.
!> @todo Err msg PA0438: flight can be inside topography -> extra check?
!--------------------------------------------------------------------------------!
 MODULE flight_mod
 
    USE control_parameters,                                                    &
        ONLY:  fl_max, num_leg, num_var_fl, num_var_fl_user, virtual_flight
 
    USE kinds

    CHARACTER(LEN=6), DIMENSION(fl_max) ::  leg_mode = 'cyclic'  !< flight mode through the model domain, either 'cyclic' or 'return' 

    INTEGER(iwp) ::  l           !< index for flight leg
    INTEGER(iwp) ::  var_index   !< index for measured variable

    LOGICAL, DIMENSION(:), ALLOCATABLE  ::  cyclic_leg !< flag to identify fly mode

    REAL(wp) ::  flight_end = 9999999.9_wp  !< end time of virtual flight
    REAL(wp) ::  flight_begin = 0.0_wp      !< end time of virtual flight

    REAL(wp), DIMENSION(fl_max) ::  flight_angle = 45.0_wp   !< angle determining the horizontal flight direction
    REAL(wp), DIMENSION(fl_max) ::  flight_level = 100.0_wp  !< flight level 
    REAL(wp), DIMENSION(fl_max) ::  max_elev_change = 0.0_wp !< maximum elevation change for the respective flight leg 
    REAL(wp), DIMENSION(fl_max) ::  rate_of_climb = 0.0_wp   !< rate of climb or descent
    REAL(wp), DIMENSION(fl_max) ::  speed_agl = 25.0_wp      !< absolute horizontal flight speed above ground level (agl)
    REAL(wp), DIMENSION(fl_max) ::  x_start = 999999999.0_wp !< start x position 
    REAL(wp), DIMENSION(fl_max) ::  x_end   = 999999999.0_wp !< end x position
    REAL(wp), DIMENSION(fl_max) ::  y_start = 999999999.0_wp !< start y position
    REAL(wp), DIMENSION(fl_max) ::  y_end   = 999999999.0_wp !< end y position

    REAL(wp), DIMENSION(:), ALLOCATABLE ::  u_agl      !< u-component of flight speed
    REAL(wp), DIMENSION(:), ALLOCATABLE ::  v_agl      !< v-component of flight speed
    REAL(wp), DIMENSION(:), ALLOCATABLE ::  w_agl      !< w-component of flight speed
    REAL(wp), DIMENSION(:), ALLOCATABLE ::  x_pos      !< aircraft x-position 
    REAL(wp), DIMENSION(:), ALLOCATABLE ::  y_pos      !< aircraft y-position 
    REAL(wp), DIMENSION(:), ALLOCATABLE ::  z_pos      !< aircraft z-position 

    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sensor_l !< measured data on local PE
    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sensor   !< measured data

    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  var_u  !< dummy array for possibly user-defined quantities

    SAVE

    PRIVATE

    INTERFACE flight_header
       MODULE PROCEDURE flight_header
    END INTERFACE flight_header
    
    INTERFACE flight_init
       MODULE PROCEDURE flight_init
    END INTERFACE flight_init

    INTERFACE flight_init_output
       MODULE PROCEDURE flight_init_output
    END INTERFACE flight_init_output

    INTERFACE flight_check_parameters
       MODULE PROCEDURE flight_check_parameters
    END INTERFACE flight_check_parameters

    INTERFACE flight_parin
       MODULE PROCEDURE flight_parin
    END INTERFACE flight_parin

    INTERFACE interpolate_xyz
       MODULE PROCEDURE interpolate_xyz
    END INTERFACE interpolate_xyz

    INTERFACE flight_measurement
       MODULE PROCEDURE flight_measurement
    END INTERFACE flight_measurement
    
    INTERFACE flight_rrd_global 
       MODULE PROCEDURE flight_rrd_global
    END INTERFACE flight_rrd_global
    
    INTERFACE flight_wrd_global 
       MODULE PROCEDURE flight_wrd_global 
    END INTERFACE flight_wrd_global 

!
!-- Private interfaces
    PRIVATE flight_check_parameters, flight_init_output, interpolate_xyz
!
!-- Public interfaces
    PUBLIC flight_init, flight_header, flight_parin, flight_measurement,       &
           flight_wrd_global, flight_rrd_global                   
!
!-- Public variables
    PUBLIC fl_max, sensor, x_pos, y_pos, z_pos

 CONTAINS

!------------------------------------------------------------------------------!
! Description:
! ------------
!> Header output for flight module.
!------------------------------------------------------------------------------!
    SUBROUTINE flight_header ( io )

    
       IMPLICIT NONE

       INTEGER(iwp), INTENT(IN) ::  io            !< Unit of the output file

       WRITE ( io, 1  )
       WRITE ( io, 2  )
       WRITE ( io, 3  ) num_leg
       WRITE ( io, 4  ) flight_begin
       WRITE ( io, 5  ) flight_end
       
       DO l=1, num_leg
          WRITE ( io, 6   ) l
          WRITE ( io, 7   ) speed_agl(l)
          WRITE ( io, 8   ) flight_level(l)
          WRITE ( io, 9   ) max_elev_change(l)
          WRITE ( io, 10  ) rate_of_climb(l)
          WRITE ( io, 11  ) leg_mode(l)
       ENDDO

       
     1   FORMAT (' Virtual flights:'/                                           &
               ' ----------------')
     2   FORMAT ('       Output every timestep')
     3   FORMAT ('       Number of flight legs:',    I3                )
     4   FORMAT ('       Begin of measurements:',    F10.1    , ' s'   )
     5   FORMAT ('       End of measurements:',      F10.1    , ' s'   )
     6   FORMAT ('       Leg', I3/,                                             &
                '       ------' )
     7   FORMAT ('          Flight speed            : ', F5.1, ' m/s' )
     8   FORMAT ('          Flight level            : ', F5.1, ' m'   )
     9   FORMAT ('          Maximum elevation change: ', F5.1, ' m/s' )
     10  FORMAT ('          Rate of climb / descent : ', F5.1, ' m/s' )
     11  FORMAT ('          Leg mode                : ', A/           )
 
    END SUBROUTINE flight_header 
 
!------------------------------------------------------------------------------!
! Description:
! ------------
!> Reads the namelist flight_par.
!------------------------------------------------------------------------------!
    SUBROUTINE flight_parin 

       USE control_parameters,                                                 &
           ONLY:  message_string
      
       IMPLICIT NONE

       CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file

       NAMELIST /flight_par/  flight_angle, flight_end, flight_begin, leg_mode,&
                              flight_level, max_elev_change, rate_of_climb,    &
                              speed_agl, x_end, x_start, y_end, y_start
  

       NAMELIST /virtual_flight_parameters/                                    &
                              flight_angle, flight_end, flight_begin, leg_mode,&
                              flight_level, max_elev_change, rate_of_climb,    &
                              speed_agl, x_end, x_start, y_end, y_start
!
!--    Try to find the namelist flight_par
       REWIND ( 11 )
       line = ' '
       DO WHILE ( INDEX( line, '&virtual_flight_parameters' ) == 0 )
          READ ( 11, '(A)', END=12 )  line
       ENDDO
       BACKSPACE ( 11 )

!
!--    Read namelist
       READ ( 11, virtual_flight_parameters, ERR = 10 )
!
!--    Set switch that virtual flights shall be carried out
       virtual_flight = .TRUE.

       GOTO 14

 10    BACKSPACE( 11 )
       READ( 11 , '(A)') line
       CALL parin_fail_message( 'virtual_flight_parameters', line )
!
!--    Try to find the old namelist
 12    REWIND ( 11 )
       line = ' '
       DO WHILE ( INDEX( line, '&flight_par' ) == 0 )
          READ ( 11, '(A)', END=14 )  line
       ENDDO
       BACKSPACE ( 11 )

!
!--    Read namelist
       READ ( 11, flight_par, ERR = 13, END = 14 )
       
       message_string = 'namelist flight_par is deprecated and will be ' // &
                     'removed in near future.& Please use namelist ' //     &
                     'virtual_flight_parameters instead' 
       CALL message( 'flight_parin', 'PA0487', 0, 1, 0, 6, 0 )     
!
!--    Set switch that virtual flights shall be carried out
       virtual_flight = .TRUE.

       GOTO 14

 13    BACKSPACE( 11 )
       READ( 11 , '(A)') line
       CALL parin_fail_message( 'flight_par', line )

 14    CONTINUE

    END SUBROUTINE flight_parin

!------------------------------------------------------------------------------!
! Description:
! ------------
!> Inititalization of required arrays, number of legs and flags. Moreover,
!> initialize flight speed and -direction, as well as start positions.
!------------------------------------------------------------------------------!
    SUBROUTINE flight_init
 
       USE basic_constants_and_equations_mod,                                  &
           ONLY:  pi
    
       USE control_parameters,                                                 &
           ONLY:  initializing_actions 
                  
       USE indices,                                                            &
           ONLY:  nxlg, nxrg, nysg, nyng, nzb, nzt

       IMPLICIT NONE

       REAL(wp) ::  distance  !< distance between start and end position of a flight leg
!
!--    Determine the number of flight legs
       l = 1
       DO WHILE ( x_start(l) /= 999999999.0_wp  .AND.  l <= SIZE(x_start) )
          l       = l + 1
       ENDDO
       num_leg = l-1
!
!--    Check for proper parameter settings
       CALL flight_check_parameters
!
!--    Allocate and initialize logical array for flight pattern
       ALLOCATE( cyclic_leg(1:num_leg) )
!
!--    Initialize flags for cyclic/return legs
       DO l = 1, num_leg
          cyclic_leg(l) = MERGE( .TRUE., .FALSE.,                              &
                                 TRIM( leg_mode(l) ) == 'cyclic'               &
                               )
       ENDDO
!
!--    Allocate and initialize arraxs for flight position and speed. In case of 
!--    restart runs these data are read by the routine read_flight_restart_data
!--    instead.
       IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN 
       
          ALLOCATE( x_pos(1:num_leg), y_pos(1:num_leg ), z_pos(1:num_leg) )
!
!--       Inititalize x-, y-, and z-positions with initial start position          
          x_pos(1:num_leg) = x_start(1:num_leg)
          y_pos(1:num_leg) = y_start(1:num_leg)
          z_pos(1:num_leg) = flight_level(1:num_leg)
!
!--       Allocate arrays for flight-speed components
          ALLOCATE( u_agl(1:num_leg),                                          &
                    v_agl(1:num_leg),                                          &
                    w_agl(1:num_leg) )
!
!--       Inititalize u-, v- and w-component. 
          DO  l = 1, num_leg
!
!--          In case of return-legs, the flight direction, i.e. the horizontal
!--          flight-speed components, are derived from the given start/end 
!--          positions.
             IF (  .NOT.  cyclic_leg(l) )  THEN
                distance = SQRT( ( x_end(l) - x_start(l) )**2                  &
                               + ( y_end(l) - y_start(l) )**2 )
                u_agl(l) = speed_agl(l) * ( x_end(l) - x_start(l) ) / distance
                v_agl(l) = speed_agl(l) * ( y_end(l) - y_start(l) ) / distance
                w_agl(l) = rate_of_climb(l)
!
!--          In case of cyclic-legs, flight direction is directly derived from 
!--          the given flight angle.
             ELSE
                u_agl(l) = speed_agl(l) * COS( flight_angle(l) * pi / 180.0_wp )
                v_agl(l) = speed_agl(l) * SIN( flight_angle(l) * pi / 180.0_wp )
                w_agl(l) = rate_of_climb(l)
             ENDIF

          ENDDO
             
       ENDIF    
!
!--    Initialized data output
       CALL flight_init_output       
!
!--    Allocate array required for user-defined quantities if necessary.
       IF ( num_var_fl_user  > 0 )                                             &
          ALLOCATE( var_u(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
!
!--    Allocate and initialize arrays containing the measured data
       ALLOCATE( sensor_l(1:num_var_fl,1:num_leg) )
       ALLOCATE( sensor(1:num_var_fl,1:num_leg)   )
       sensor_l = 0.0_wp
       sensor   = 0.0_wp

    END SUBROUTINE flight_init
    
    
!------------------------------------------------------------------------------!
! Description:
! ------------
!> Initialization of output-variable names and units.
!------------------------------------------------------------------------------!
    SUBROUTINE flight_init_output
    
       USE control_parameters,                                                 &
          ONLY:  cloud_droplets, humidity, neutral, passive_scalar

       USE bulk_cloud_model_mod,                                               &
           ONLY:  bulk_cloud_model

       USE netcdf_interface
    
       IMPLICIT NONE

       CHARACTER(LEN=10) ::  label_leg  !< dummy argument to convert integer to string
       
       INTEGER(iwp) ::  i               !< loop variable
       INTEGER(iwp) ::  id_pt           !< identifyer for labeling
       INTEGER(iwp) ::  id_q            !< identifyer for labeling
       INTEGER(iwp) ::  id_ql           !< identifyer for labeling
       INTEGER(iwp) ::  id_s            !< identifyer for labeling        
       INTEGER(iwp) ::  id_u = 1        !< identifyer for labeling  
       INTEGER(iwp) ::  id_v = 2        !< identifyer for labeling
       INTEGER(iwp) ::  id_w = 3        !< identifyer for labeling
       INTEGER(iwp) ::  k               !< index variable
       
       LOGICAL      ::  init = .TRUE.   !< flag to distiquish calls of user_init_flight
!
!--    Define output quanities, at least three variables are measured (u,v,w)
       num_var_fl = 3
       IF ( .NOT. neutral                     )  THEN
          num_var_fl = num_var_fl + 1
          id_pt      = num_var_fl
       ENDIF
       IF ( humidity                          )  THEN
          num_var_fl = num_var_fl + 1
          id_q       = num_var_fl
       ENDIF
       IF ( bulk_cloud_model .OR. cloud_droplets )  THEN
          num_var_fl = num_var_fl + 1  
          id_ql      = num_var_fl
       ENDIF
       IF ( passive_scalar                    )  THEN
          num_var_fl = num_var_fl + 1
          id_s       = num_var_fl
       ENDIF
!
!--    Write output strings for dimensions x, y, z
       DO l=1, num_leg

          IF ( l < 10                    )  WRITE( label_leg, '(I1)')  l
          IF ( l >= 10   .AND.  l < 100  )  WRITE( label_leg, '(I2)')  l
          IF ( l >= 100  .AND.  l < 1000 )  WRITE( label_leg, '(I3)')  l

          dofl_dim_label_x(l)  = 'x_' // TRIM( label_leg )
          dofl_dim_label_y(l)  = 'y_' // TRIM( label_leg )
          dofl_dim_label_z(l)  = 'z_' // TRIM( label_leg )

       ENDDO
       
!
!--    Call user routine to initialize further variables
       CALL user_init_flight( init )
!
!--    Write output labels and units for the quanities
       k = 1
       DO l=1, num_leg
       
          IF ( l < 10                    )  WRITE( label_leg, '(I1)')  l
          IF ( l >= 10   .AND.  l < 100  )  WRITE( label_leg, '(I2)')  l
          IF ( l >= 100  .AND.  l < 1000 )  WRITE( label_leg, '(I3)')  l
          
          label_leg = 'leg_' // TRIM(label_leg) 
          DO i=1, num_var_fl

             IF ( i == id_u      )  THEN         
                dofl_label(k) = TRIM( label_leg ) // '_u'
                dofl_unit(k)  = 'm/s'
                k             = k + 1
             ELSEIF ( i == id_v  )  THEN       
                dofl_label(k) = TRIM( label_leg ) // '_v'
                dofl_unit(k)  = 'm/s'
                k             = k + 1
             ELSEIF ( i == id_w  )  THEN          
                dofl_label(k) = TRIM( label_leg ) // '_w'
                dofl_unit(k)  = 'm/s'
                k             = k + 1
             ELSEIF ( i == id_pt )  THEN       
                dofl_label(k) = TRIM( label_leg ) // '_pt'
                dofl_unit(k)  = 'K'
                k             = k + 1
             ELSEIF ( i == id_q  )  THEN       
                dofl_label(k) = TRIM( label_leg ) // '_q'
                dofl_unit(k)  = 'kg/kg'
                k             = k + 1
             ELSEIF ( i == id_ql )  THEN       
                dofl_label(k) = TRIM( label_leg ) // '_ql'
                dofl_unit(k)  = 'kg/kg'
                k             = k + 1
             ELSEIF ( i == id_s  )  THEN                          
                dofl_label(k) = TRIM( label_leg ) // '_s'
                dofl_unit(k)  = 'kg/kg'
                k             = k + 1
             ENDIF
          ENDDO
          
          DO i=1, num_var_fl_user
             CALL user_init_flight( init, k, i, label_leg )
          ENDDO
          
       ENDDO 
!
!--    Finally, set the total number of flight-output quantities.
       num_var_fl = num_var_fl + num_var_fl_user
       
    END SUBROUTINE flight_init_output

!------------------------------------------------------------------------------!
! Description:
! ------------
!> This routine calculates the current flight positions and calls the 
!> respective interpolation routine to measures the data at the current 
!> flight position.
!------------------------------------------------------------------------------!
    SUBROUTINE flight_measurement

       USE arrays_3d,                                                          &
           ONLY:  ddzu, ddzw, pt, q, ql, s, u, v, w, zu, zw

       USE control_parameters,                                                 &
           ONLY:  cloud_droplets, dt_3d, humidity, neutral,                    &
                  passive_scalar, simulated_time
                  
       USE cpulog,                                                             &
           ONLY:  cpu_log, log_point

       USE grid_variables,                                                     &
           ONLY:  ddx, ddy, dx, dy

       USE indices,                                                            &
           ONLY:  nx, nxl, nxr, ny, nys, nyn

       USE bulk_cloud_model_mod,                                               &
           ONLY:  bulk_cloud_model

       USE pegrid

       IMPLICIT NONE

       LOGICAL  ::  on_pe  !< flag to check if current flight position is on current PE

       REAL(wp) ::  x  !< distance between left edge of current grid box and flight position
       REAL(wp) ::  y  !< distance between south edge of current grid box and flight position

       INTEGER(iwp) ::  i   !< index of current grid box along x
       INTEGER(iwp) ::  j   !< index of current grid box along y
       INTEGER(iwp) ::  n   !< loop variable for number of user-defined output quantities
       
       CALL cpu_log( log_point(65), 'flight_measurement', 'start' )
!
!--    Perform flight measurement if start time is reached.
       IF ( simulated_time >= flight_begin  .AND.                              &
            simulated_time <= flight_end )  THEN

          sensor_l = 0.0_wp
          sensor   = 0.0_wp
!
!--       Loop over all flight legs
          DO l=1, num_leg
!
!--          Update location for each leg
             x_pos(l) = x_pos(l) + u_agl(l) * dt_3d 
             y_pos(l) = y_pos(l) + v_agl(l) * dt_3d 
             z_pos(l) = z_pos(l) + w_agl(l) * dt_3d
!
!--          Check if location must be modified for return legs.  
!--          Carry out horizontal reflection if required.
             IF ( .NOT. cyclic_leg(l) )  THEN
!
!--             Outward flight, i.e. from start to end
                IF ( u_agl(l) >= 0.0_wp  .AND.  x_pos(l) > x_end(l)      )  THEN
                   x_pos(l) = 2.0_wp * x_end(l)   - x_pos(l)
                   u_agl(l) = - u_agl(l)
!
!--             Return flight, i.e. from end to start
                ELSEIF ( u_agl(l) < 0.0_wp  .AND.  x_pos(l) < x_start(l) )  THEN 
                   x_pos(l) = 2.0_wp * x_start(l) - x_pos(l)
                   u_agl(l) = - u_agl(l)
                ENDIF
!
!--             Outward flight, i.e. from start to end
                IF ( v_agl(l) >= 0.0_wp  .AND.  y_pos(l) > y_end(l)      )  THEN
                   y_pos(l) = 2.0_wp * y_end(l)   - y_pos(l)
                   v_agl(l) = - v_agl(l)
!
!--             Return flight, i.e. from end to start                  
                ELSEIF ( v_agl(l) < 0.0_wp  .AND.  y_pos(l) < y_start(l) )  THEN 
                   y_pos(l) = 2.0_wp * y_start(l) - y_pos(l)
                   v_agl(l) = - v_agl(l)
                ENDIF
!
!--          Check if flight position is out of the model domain and apply
!--          cyclic conditions if required
             ELSEIF ( cyclic_leg(l) )  THEN
!
!--             Check if aircraft leaves the model domain at the right boundary
                IF ( ( flight_angle(l) >= 0.0_wp     .AND.                     &
                       flight_angle(l) <= 90.0_wp )  .OR.                      &  
                     ( flight_angle(l) >= 270.0_wp   .AND.                     &
                       flight_angle(l) <= 360.0_wp ) )  THEN
                   IF ( x_pos(l) >= ( nx + 0.5_wp ) * dx )                     &
                      x_pos(l) = x_pos(l) - ( nx + 1 ) * dx 
!
!--             Check if aircraft leaves the model domain at the left boundary
                ELSEIF ( flight_angle(l) > 90.0_wp  .AND.                      &
                         flight_angle(l) < 270.0_wp )  THEN
                   IF ( x_pos(l) < -0.5_wp * dx )                             &
                      x_pos(l) = ( nx + 1 ) * dx + x_pos(l)  
                ENDIF 
!
!--             Check if aircraft leaves the model domain at the north boundary
                IF ( flight_angle(l) >= 0.0_wp  .AND.                          &
                     flight_angle(l) <= 180.0_wp )  THEN
                   IF ( y_pos(l) >= ( ny + 0.5_wp ) * dy )                     &
                      y_pos(l) = y_pos(l) - ( ny + 1 ) * dy 
!
!--             Check if aircraft leaves the model domain at the south boundary
                ELSEIF ( flight_angle(l) > 180.0_wp  .AND.                     &
                         flight_angle(l) < 360.0_wp )  THEN
                   IF ( y_pos(l) < -0.5_wp * dy )                              &
                      y_pos(l) = ( ny + 1 ) * dy + y_pos(l) 
                ENDIF
                
             ENDIF
!
!--          Check if maximum elevation change is already reached. If required
!--          reflect vertically.
             IF ( rate_of_climb(l) /= 0.0_wp )  THEN
!
!--             First check if aircraft is too high
                IF (  w_agl(l) > 0.0_wp  .AND.                                 &
                      z_pos(l) - flight_level(l) > max_elev_change(l) )  THEN
                   z_pos(l) = 2.0_wp * ( flight_level(l) + max_elev_change(l) )&
                              - z_pos(l)
                   w_agl(l) = - w_agl(l)
!
!--             Check if aircraft is too low
                ELSEIF (  w_agl(l) < 0.0_wp  .AND.  z_pos(l) < flight_level(l) )  THEN
                   z_pos(l) = 2.0_wp * flight_level(l) - z_pos(l)
                   w_agl(l) = - w_agl(l)
                ENDIF
                
             ENDIF 
!
!--          Determine grid indices for flight position along x- and y-direction.
!--          Please note, there is a special treatment for the index 
!--          along z-direction, which is due to vertical grid stretching. 
             i = ( x_pos(l) + 0.5_wp * dx ) * ddx
             j = ( y_pos(l) + 0.5_wp * dy ) * ddy
!
!--          Check if indices are on current PE
             on_pe = ( i >= nxl  .AND.  i <= nxr  .AND.                        &
                       j >= nys  .AND.  j <= nyn )

             IF ( on_pe )  THEN

                var_index = 1
!
!--             Recalculate indices, as u is shifted by -0.5*dx.
                i =   x_pos(l) * ddx
                j = ( y_pos(l) + 0.5_wp * dy ) * ddy
!
!--             Calculate distance from left and south grid-box edges.
                x  = x_pos(l) - ( 0.5_wp - i ) * dx
                y  = y_pos(l) - j * dy
!
!--             Interpolate u-component onto current flight position.
                CALL interpolate_xyz( u, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
                var_index = var_index + 1
!
!--             Recalculate indices, as v is shifted by -0.5*dy.
                i = ( x_pos(l) + 0.5_wp * dx ) * ddx
                j =   y_pos(l) * ddy

                x  = x_pos(l) - i * dx
                y  = y_pos(l) - ( 0.5_wp - j ) * dy
                CALL interpolate_xyz( v, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
                var_index = var_index + 1
!
!--             Interpolate w and scalar quantities. Recalculate indices.
                i  = ( x_pos(l) + 0.5_wp * dx ) * ddx
                j  = ( y_pos(l) + 0.5_wp * dy ) * ddy
                x  = x_pos(l) - i * dx
                y  = y_pos(l) - j * dy
!
!--             Interpolate w-velocity component.
                CALL interpolate_xyz( w, zw, ddzw, 0.0_wp, x, y, var_index, j, i )
                var_index = var_index + 1
!
!--             Potential temerature
                IF ( .NOT. neutral )  THEN
                   CALL interpolate_xyz( pt, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
                   var_index = var_index + 1
                ENDIF
!
!--             Humidity
                IF ( humidity )  THEN
                   CALL interpolate_xyz( q, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
                   var_index = var_index + 1
                ENDIF
!
!--             Liquid water content
                IF ( bulk_cloud_model .OR. cloud_droplets )  THEN
                   CALL interpolate_xyz( ql, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
                   var_index = var_index + 1
                ENDIF
!
!--             Passive scalar
                IF ( passive_scalar )  THEN
                   CALL interpolate_xyz( s, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
                   var_index = var_index + 1
                ENDIF
!
!--             Treat user-defined variables if required
                DO n = 1, num_var_fl_user                
                   CALL user_flight( var_u, n )
                   CALL interpolate_xyz( var_u, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
                   var_index = var_index + 1
                ENDDO
             ENDIF

          ENDDO
!
!--       Write local data on global array.
#if defined( __parallel )
          CALL MPI_ALLREDUCE(sensor_l(1,1), sensor(1,1),                       &
                             num_var_fl*num_leg, MPI_REAL, MPI_SUM,               &
                             comm2d, ierr )
#else
          sensor     = sensor_l
#endif
       ENDIF
       
       CALL cpu_log( log_point(65), 'flight_measurement', 'stop' )

    END SUBROUTINE flight_measurement

!------------------------------------------------------------------------------!
! Description:
! ------------
!> This subroutine bi-linearly interpolates the respective data onto the current
!> flight position.
!------------------------------------------------------------------------------!
    SUBROUTINE interpolate_xyz( var, z_uw, ddz_uw, fac, x, y, var_ind, j, i )

       USE control_parameters,                                                 &
           ONLY:  dz, dz_stretch_level_start   
  
      USE grid_variables,                                                      &
          ONLY:  dx, dy
    
       USE indices,                                                            &
           ONLY:  nzb, nzt, nxlg, nxrg, nysg, nyng

       IMPLICIT NONE

       INTEGER(iwp) ::  i        !< index along x
       INTEGER(iwp) ::  j        !< index along y
       INTEGER(iwp) ::  k        !< index along z
       INTEGER(iwp) ::  k1       !< dummy variable
       INTEGER(iwp) ::  var_ind  !< index variable for output quantity

       REAL(wp) ::  aa        !< dummy argument for horizontal interpolation   
       REAL(wp) ::  bb        !< dummy argument for horizontal interpolation
       REAL(wp) ::  cc        !< dummy argument for horizontal interpolation
       REAL(wp) ::  dd        !< dummy argument for horizontal interpolation
       REAL(wp) ::  gg        !< dummy argument for horizontal interpolation
       REAL(wp) ::  fac       !< flag to indentify if quantity is on zu or zw level
       REAL(wp) ::  var_int   !< horizontal interpolated variable at current position
       REAL(wp) ::  var_int_l !< horizontal interpolated variable at k-level
       REAL(wp) ::  var_int_u !< horizontal interpolated variable at (k+1)-level
       REAL(wp) ::  x         !< distance between left edge of current grid box and flight position
       REAL(wp) ::  y         !< distance between south edge of current grid box and flight position

       REAL(wp), DIMENSION(1:nzt+1)   ::  ddz_uw !< inverse vertical grid spacing
       REAL(wp), DIMENSION(nzb:nzt+1) ::  z_uw   !< height level

       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var !< treted quantity
!
!--    Calculate interpolation coefficients
       aa = x**2          + y**2
       bb = ( dx - x )**2 + y**2
       cc = x**2          + ( dy - y )**2
       dd = ( dx - x )**2 + ( dy - y )**2
       gg = aa + bb + cc + dd
!
!--    Obtain vertical index. Special treatment for grid index along z-direction
!--    if flight position is above the vertical grid-stretching level.
!--    fac=1 if variable is on scalar grid level, fac=0 for w-component.
       IF ( z_pos(l) < dz_stretch_level_start(1) )  THEN
          k = ( z_pos(l) + fac * 0.5_wp * dz(1) ) / dz(1)
       ELSE
!
!--       Search for k-index
          DO k1=nzb, nzt
             IF ( z_pos(l) >= z_uw(k1) .AND. z_pos(l) < z_uw(k1+1) )  THEN
                k = k1
                EXIT
             ENDIF                   
          ENDDO
       ENDIF
!
!--    (x,y)-interpolation at k-level
       var_int_l = ( ( gg - aa ) * var(k,j,i)       +                          &
                     ( gg - bb ) * var(k,j,i+1)     +                          &
                     ( gg - cc ) * var(k,j+1,i)     +                          &
                     ( gg - dd ) * var(k,j+1,i+1)                              &
                   ) / ( 3.0_wp * gg )
!
!--    (x,y)-interpolation on (k+1)-level
       var_int_u = ( ( gg - aa ) * var(k+1,j,i)     +                          &
                     ( gg - bb ) * var(k+1,j,i+1)   +                          &
                     ( gg - cc ) * var(k+1,j+1,i)   +                          &
                     ( gg - dd ) * var(k+1,j+1,i+1)                            &
                   ) / ( 3.0_wp * gg )
!
!--    z-interpolation onto current flight postion
       var_int = var_int_l                                                     &
                           + ( z_pos(l) - z_uw(k) ) * ddz_uw(k+1)              &
                           * (var_int_u - var_int_l )
!
!--    Store on local data array
       sensor_l(var_ind,l) = var_int

    END SUBROUTINE interpolate_xyz

!------------------------------------------------------------------------------!
! Description:
! ------------
!> Perform parameter checks.
!------------------------------------------------------------------------------!
    SUBROUTINE flight_check_parameters

       USE arrays_3d,                                                          &
           ONLY:  zu
    
       USE control_parameters,                                                 &
           ONLY:  bc_lr_cyc, bc_ns_cyc, message_string

       USE grid_variables,                                                     &
           ONLY:  dx, dy   
          
       USE indices,                                                            &
           ONLY:  nx, ny, nz
           
       USE netcdf_interface,                                                   &
           ONLY:  netcdf_data_format

       IMPLICIT NONE
       
!
!--    Check if start positions are properly set.
       DO l=1, num_leg
          IF ( x_start(l) < 0.0_wp  .OR.  x_start(l) > ( nx + 1 ) * dx )  THEN
             message_string = 'Start x position is outside the model domain'
             CALL message( 'flight_check_parameters', 'PA0431', 1, 2, 0, 6, 0 )
          ENDIF 
          IF ( y_start(l) < 0.0_wp  .OR.  y_start(l) > ( ny + 1 ) * dy )  THEN
             message_string = 'Start y position is outside the model domain'
             CALL message( 'flight_check_parameters', 'PA0432', 1, 2, 0, 6, 0 )
          ENDIF 

       ENDDO
!
!--    Check for leg mode
       DO l=1, num_leg
!
!--       Check if leg mode matches the overall lateral model boundary 
!--       conditions.
          IF ( TRIM( leg_mode(l) ) == 'cyclic' )  THEN
             IF ( .NOT. bc_lr_cyc  .OR.  .NOT. bc_ns_cyc )  THEN
                message_string = 'Cyclic flight leg does not match ' //        &
                                 'lateral boundary condition'
                CALL message( 'flight_check_parameters', 'PA0433', 1, 2, 0, 6, 0 )
             ENDIF 
!
!--       Check if end-positions are inside the model domain in case of 
!..       return-legs.  
          ELSEIF ( TRIM( leg_mode(l) ) == 'return' )  THEN
             IF ( x_end(l) > ( nx + 1 ) * dx  .OR.                             &
                  y_end(l) > ( ny + 1 ) * dx )  THEN
                message_string = 'Flight leg or parts of it are outside ' //   &
                                 'the model domain'
                CALL message( 'flight_check_parameters', 'PA0434', 1, 2, 0, 6, 0 )
             ENDIF 
          ELSE
             message_string = 'Unknown flight mode'
             CALL message( 'flight_check_parameters', 'PA0435', 1, 2, 0, 6, 0 )
          ENDIF

       ENDDO          
!
!--    Check if start and end positions are properly set in case of return legs.
       DO l=1, num_leg

          IF ( x_start(l) > x_end(l) .AND. leg_mode(l) == 'return' )  THEN
             message_string = 'x_start position must be <= x_end ' //          &
                              'position for return legs'
             CALL message( 'flight_check_parameters', 'PA0436', 1, 2, 0, 6, 0 )
          ENDIF 
          IF ( y_start(l) > y_end(l) .AND. leg_mode(l) == 'return' )  THEN
             message_string = 'y_start position must be <= y_end ' //          &
                              'position for return legs'
             CALL message( 'flight_check_parameters', 'PA0437', 1, 2, 0, 6, 0 )
          ENDIF 
       ENDDO
!
!--    Check if given flight object remains inside model domain if a rate of 
!--    climb / descent is prescribed. 
       DO l=1, num_leg
          IF ( flight_level(l) + max_elev_change(l) > zu(nz) .OR.              &
               flight_level(l) + max_elev_change(l) <= 0.0_wp )  THEN
             message_string = 'Flight level is outside the model domain '
             CALL message( 'flight_check_parameters', 'PA0438', 1, 2, 0, 6, 0 )
          ENDIF 
       ENDDO       
!
!--    Check for appropriate NetCDF format. Definition of more than one 
!--    unlimited dimension is unfortunately only possible with NetCDF4/HDF5. 
       IF (  netcdf_data_format <= 2 )  THEN
          message_string = 'netcdf_data_format must be > 2'
          CALL message( 'flight_check_parameters', 'PA0439', 1, 2, 0, 6, 0 )
       ENDIF


    END SUBROUTINE flight_check_parameters
        

!------------------------------------------------------------------------------!
! Description:
! ------------
!> This routine reads the respective restart data.
!------------------------------------------------------------------------------!
    SUBROUTINE flight_rrd_global( found )  


       USE control_parameters,                                                 &
           ONLY: length, restart_string

    
       IMPLICIT NONE
       
       LOGICAL, INTENT(OUT)  ::  found 


       found = .TRUE.


       SELECT CASE ( restart_string(1:length) )
          
          CASE ( 'u_agl' )
             IF ( .NOT. ALLOCATED( u_agl ) )  ALLOCATE( u_agl(1:num_leg) )
             READ ( 13 )  u_agl   
          CASE ( 'v_agl' )
             IF ( .NOT. ALLOCATED( v_agl ) )  ALLOCATE( v_agl(1:num_leg) )
             READ ( 13 )  v_agl
          CASE ( 'w_agl' )
             IF ( .NOT. ALLOCATED( w_agl ) )  ALLOCATE( w_agl(1:num_leg) )
             READ ( 13 )  w_agl
          CASE ( 'x_pos' )
             IF ( .NOT. ALLOCATED( x_pos ) )  ALLOCATE( x_pos(1:num_leg) )
             READ ( 13 )  x_pos
          CASE ( 'y_pos' )
             IF ( .NOT. ALLOCATED( y_pos ) )  ALLOCATE( y_pos(1:num_leg) )
             READ ( 13 )  y_pos
          CASE ( 'z_pos' )
             IF ( .NOT. ALLOCATED( z_pos ) )  ALLOCATE( z_pos(1:num_leg) )
             READ ( 13 )  z_pos

          CASE DEFAULT

             found = .FALSE.
          
       END SELECT


    END SUBROUTINE flight_rrd_global  
    
!------------------------------------------------------------------------------!
! Description:
! ------------
!> This routine writes the respective restart data.
!------------------------------------------------------------------------------!
    SUBROUTINE flight_wrd_global  


       IMPLICIT NONE
 
       CALL wrd_write_string( 'u_agl' )
       WRITE ( 14 )  u_agl

       CALL wrd_write_string( 'v_agl' )

       WRITE ( 14 )  v_agl

       CALL wrd_write_string( 'w_agl' )
       WRITE ( 14 )  w_agl

       CALL wrd_write_string( 'x_pos' )
       WRITE ( 14 )  x_pos

       CALL wrd_write_string( 'y_pos' )
       WRITE ( 14 )  y_pos

       CALL wrd_write_string( 'z_pos' )
       WRITE ( 14 )  z_pos
       
    END SUBROUTINE flight_wrd_global    
    

 END MODULE flight_mod