Changeset 4497 for palm


Ignore:
Timestamp:
Apr 15, 2020 10:20:51 AM (7 months ago)
Author:
raasch
Message:

last bugfix deactivated because of compile problems, files re-formatted to follow the PALM coding standard

Location:
palm/trunk/SOURCE
Files:
12 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/restart_data_mpi_io_mod.f90

    r4496 r4497  
    2424! -----------------
    2525! $Id$
     26! last bugfix deactivated because of compile problems
     27!
     28! 4496 2020-04-15 08:37:26Z raasch
    2629! problem with posix read arguments for surface data fixed
    2730!
     
    17731776#else
    17741777                CALL posix_lseek( fh, disp_f )
    1775                 CALL posix_read( fh, data(m_start_index(j_f:,i_f:)), nr_bytes_f )
     1778!                CALL posix_read( fh, data(m_start_index(j_f:,i_f:)), nr_bytes_f )
    17761779#endif
    17771780                disp_f     = disp
  • palm/trunk/SOURCE/user_init_land_surface.f90

    r4360 r4497  
    11!> @file user_init_land_surface.f90
    2 !------------------------------------------------------------------------------!
     2!--------------------------------------------------------------------------------------------------!
    33! This file is part of the PALM model system.
    44!
    5 ! PALM is free software: you can redistribute it and/or modify it under the
    6 ! terms of the GNU General Public License as published by the Free Software
    7 ! Foundation, either version 3 of the License, or (at your option) any later
    8 ! version.
     5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
     6! Public License as published by the Free Software Foundation, either version 3 of the License, or
     7! (at your option) any later version.
    98!
    10 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
    11 ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
    12 ! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
     9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
     10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
     11! Public License for more details.
    1312!
    14 ! You should have received a copy of the GNU General Public License along with
    15 ! PALM. If not, see <http://www.gnu.org/licenses/>.
     13! You should have received a copy of the GNU General Public License along with PALM. If not, see
     14! <http://www.gnu.org/licenses/>.
    1615!
    1716! Copyright 1997-2020 Leibniz Universitaet Hannover
    18 !------------------------------------------------------------------------------!
     17!--------------------------------------------------------------------------------------------------!
     18!
    1919!
    2020! Current revisions:
     
    2525! -----------------
    2626! $Id$
     27! file re-formatted to follow the PALM coding standard
     28!
     29! Current revisions:
     30! -----------------
     31!
     32!
     33! 4360 2020-01-07 11:25:50Z suehring
    2734! Corrected "Former revisions" section
    28 ! 
     35!
    2936! 3766 2019-02-26 16:23:41Z raasch
    3037! unused variables commented out
    31 ! 
     38!
    3239! 3700 2019-01-26 17:03:42Z knoop
    3340! Corrected "Former revisions" section
    34 ! 
     41!
    3542! 1496 2014-12-02 17:25:50Z maronga
    3643! Initial revision
    37 ! 
     44!
    3845! Description:
    3946! ------------
    4047!> Execution of user-defined actions to initiate the land surface model
    41 !------------------------------------------------------------------------------!
     48!--------------------------------------------------------------------------------------------------!
    4249 SUBROUTINE user_init_land_surface
    43  
     50
    4451
    4552    USE kinds
  • palm/trunk/SOURCE/user_init_plant_canopy.f90

    r4360 r4497  
    11!> @file user_init_plant_canopy.f90
    2 !------------------------------------------------------------------------------!
     2!--------------------------------------------------------------------------------------------------!
    33! This file is part of the PALM model system.
    44!
    5 ! PALM is free software: you can redistribute it and/or modify it under the
    6 ! terms of the GNU General Public License as published by the Free Software
    7 ! Foundation, either version 3 of the License, or (at your option) any later
    8 ! version.
     5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
     6! Public License as published by the Free Software Foundation, either version 3 of the License, or
     7! (at your option) any later version.
    98!
    10 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
    11 ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
    12 ! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
     9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
     10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
     11! Public License for more details.
    1312!
    14 ! You should have received a copy of the GNU General Public License along with
    15 ! PALM. If not, see <http://www.gnu.org/licenses/>.
     13! You should have received a copy of the GNU General Public License along with PALM. If not, see
     14! <http://www.gnu.org/licenses/>.
    1615!
    1716! Copyright 1997-2020 Leibniz Universitaet Hannover
    18 !------------------------------------------------------------------------------!
     17!--------------------------------------------------------------------------------------------------!
     18!
    1919!
    2020! Current revisions:
     
    2525! -----------------
    2626! $Id$
     27! file re-formatted to follow the PALM coding standard
     28!
     29! Current revisions:
     30! -----------------
     31!
     32!
     33! 4360 2020-01-07 11:25:50Z suehring
    2734! Renamed canopy_mode 'block' to 'homogeneous'
    28 ! 
     35!
    2936! 4182 2019-08-22 15:20:23Z scharf
    3037! Corrected "Former revisions" section
    31 ! 
     38!
    3239! 3768 2019-02-27 14:35:58Z raasch
    3340! unused variables commented out to avoid compiler warnings
    34 ! 
     41!
    3542! 3655 2019-01-07 16:51:22Z knoop
    3643! Corrected "Former revisions" section
     
    4148! Description:
    4249! ------------
    43 !> Initialisation of the leaf area density array (for scalar grid points) and
    44 !> the array of the canopy drag coefficient, if the user has not chosen any
    45 !> of the default cases
    46 !------------------------------------------------------------------------------!
     50!> Initialisation of the leaf area density array (for scalar grid points) and the array of the
     51!> canopy drag coefficient, if the user has not chosen any of the default cases
     52!--------------------------------------------------------------------------------------------------!
    4753 SUBROUTINE user_init_plant_canopy
    48  
     54
    4955
    5056    USE arrays_3d
    51    
     57
    5258    USE control_parameters
    53    
     59
    5460    USE indices
    55    
     61
    5662    USE kinds
    5763
    5864    USE plant_canopy_model_mod
    59    
     65
    6066    USE user
    6167
    6268    IMPLICIT NONE
    6369
    64 !    INTEGER(iwp) :: i   !< running index
    65 !    INTEGER(iwp) :: j   !< running index
     70!    INTEGER(iwp) :: i  !< running index
     71!    INTEGER(iwp) :: j  !< running index
    6672
    6773!
     
    7884       CASE ( 'user_defined_canopy_1' )
    7985!
    80 !--       Here the user can define his own forest topography.
    81 !--       The following lines contain an example, where the plant canopy extends
    82 !--       only over the second half of the model domain along x.
    83 !--       Attention: DO-loops have to include the ghost points (nxlg-nxrg,
    84 !--       nysg-nyng), because no exchange of ghost point information is intended,
     86!--       Here the user can define his own forest topography.
     87!--       The following lines contain an example, where the plant canopy extends only over the
     88!--       second half of the model domain along x. Attention: DO-loops have to include the ghost
     89!--       points (nxlg-nxrg, nysg-nyng), because no exchange of ghost point information is intended,
    8590!--       in order to minimize communication between CPUs.
    8691!          DO  i = nxlg, nxrg
    87 !             IF ( i >= INT( ( nx+1 ) / 2 ) ) THEN
     92!             IF ( i >= INT( ( nx+1 ) / 2 ) )  THEN
    8893!                DO  j = nysg, nyng
    8994!                   lad_s(:,j,i) = lad(:)
     
    9297!                lad_s(:,:,i) = 0.0_wp
    9398!             ENDIF
    94 !          ENDDO 
     99!          ENDDO
    95100!
    96 !--       After definition, please
    97 !--       remove the following three lines!
    98           message_string = 'canopy_mode "' // canopy_mode // &
    99                            '" not available yet'
     101!--       After definition, please remove the following three lines!
     102          message_string = 'canopy_mode "' // canopy_mode //  '" not available yet'
    100103          CALL message( 'user_init_plant_canopy', 'UI0007', 0, 1, 0, 6, 0 )
    101          
     104
    102105       CASE DEFAULT
    103106!
    104 !--       The DEFAULT case is reached if the parameter canopy_mode contains a
    105 !--       wrong character string that is neither recognized in init_3d_model nor
    106 !--       here in user_init_plant_canopy.
     107!--       The DEFAULT case is reached if the parameter canopy_mode contains a wrong character string
     108!--       that is neither recognized in init_3d_model nor here in user_init_plant_canopy.
    107109          message_string = 'unknown canopy_mode "' // canopy_mode // '"'
    108110          CALL message( 'user_init_plant_canopy', 'UI0008', 1, 2, 0, 6, 0 )
  • palm/trunk/SOURCE/user_init_radiation.f90

    r4360 r4497  
    11!> @file user_init_radiation.f90
    2 !------------------------------------------------------------------------------!
     2!--------------------------------------------------------------------------------------------------!
    33! This file is part of the PALM model system.
    44!
    5 ! PALM is free software: you can redistribute it and/or modify it under the
    6 ! terms of the GNU General Public License as published by the Free Software
    7 ! Foundation, either version 3 of the License, or (at your option) any later
    8 ! version.
     5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
     6! Public License as published by the Free Software Foundation, either version 3 of the License, or
     7! (at your option) any later version.
    98!
    10 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
    11 ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
    12 ! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
     9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
     10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
     11! Public License for more details.
    1312!
    14 ! You should have received a copy of the GNU General Public License along with
    15 ! PALM. If not, see <http://www.gnu.org/licenses/>.
     13! You should have received a copy of the GNU General Public License along with PALM. If not, see
     14! <http://www.gnu.org/licenses/>.
    1615!
    1716! Copyright 1997-2020 Leibniz Universitaet Hannover
    18 !------------------------------------------------------------------------------!
     17!--------------------------------------------------------------------------------------------------!
     18!
    1919!
    2020! Current revisions:
     
    2525! -----------------
    2626! $Id$
     27! file re-formatted to follow the PALM coding standard
     28!
     29!
     30! 4360 2020-01-07 11:25:50Z suehring
    2731! Corrected "Former revisions" section
    28 ! 
     32!
    2933! 3768 2019-02-27 14:35:58Z raasch
    3034! unused variables commented out to avoid compiler warnings
    31 ! 
     35!
    3236! 3655 2019-01-07 16:51:22Z knoop
    3337! Corrected "Former revisions" section
    34 ! 
     38!
    3539! 1585 2015-04-30 07:05:52Z maronga
    3640! Initial revision
    37 ! 
     41!
    3842! Description:
    3943! ------------
    4044!> Execution of user-defined actions to initiate the radiation model
    41 !------------------------------------------------------------------------------!
     45!--------------------------------------------------------------------------------------------------!
    4246 SUBROUTINE user_init_radiation
    43  
     47
    4448
    4549    USE arrays_3d
    46    
     50
    4751    USE control_parameters
    48    
     52
    4953    USE indices
    50    
     54
    5155    USE kinds
    5256
    5357    USE radiation_model_mod
    54    
     58
    5559    USE user
    5660
    5761    IMPLICIT NONE
    5862
    59 !    INTEGER(iwp) :: i   !< running index
    60 !    INTEGER(iwp) :: j   !< running index
     63!    INTEGER(iwp) :: i  !< running index
     64!    INTEGER(iwp) :: j  !< running index
    6165
    6266!
  • palm/trunk/SOURCE/user_init_urban_surface.f90

    r4481 r4497  
    11!> @file user_init_urban_surface.f90
    2 !------------------------------------------------------------------------------!
     2!--------------------------------------------------------------------------------------------------!
    33! This file is part of the PALM model system.
    44!
    5 ! PALM is free software: you can redistribute it and/or modify it under the
    6 ! terms of the GNU General Public License as published by the Free Software
    7 ! Foundation, either version 3 of the License, or (at your option) any later
    8 ! version.
     5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
     6! Public License as published by the Free Software Foundation, either version 3 of the License, or
     7! (at your option) any later version.
    98!
    10 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
    11 ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
    12 ! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
     9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
     10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
     11! Public License for more details.
    1312!
    14 ! You should have received a copy of the GNU General Public License along with
    15 ! PALM. If not, see <http://www.gnu.org/licenses/>.
     13! You should have received a copy of the GNU General Public License along with PALM. If not, see
     14! <http://www.gnu.org/licenses/>.
    1615!
    17 ! Copyright 2015-2020  Technical University in Prague
    1816! Copyright 1997-2020 Leibniz Universitaet Hannover
    19 !--------------------------------------------------------------------------------!
     17!--------------------------------------------------------------------------------------------------!
     18!
    2019!
    2120! Current revisions:
     
    2625! -----------------
    2726! $Id$
     27! file re-formatted to follow the PALM coding standard
     28!
     29!
     30! 4481 2020-03-31 18:55:54Z maronga
    2831! Corrected "Former revisions" section
    29 ! 
     32!
    3033! 3768 2019-02-27 14:35:58Z raasch
    3134! unused variables commented out to avoid compiler warnings
    32 ! 
     35!
    3336! 3655 2019-01-07 16:51:22Z knoop
    3437! Corrected "Former revisions" section
    35 ! 
     38!
    3639! 2007 2016-08-24 15:47:17Z kanani
    3740! Initial revision
     
    4144! ------------
    4245!> Execution of user-defined actions to initiate the urban surface model
    43 !------------------------------------------------------------------------------!
     46!--------------------------------------------------------------------------------------------------!
    4447 SUBROUTINE user_init_urban_surface
    4548
    4649    USE arrays_3d
    47    
    48 !    USE control_parameters,                                                    &
     50
     51!    USE control_parameters,                                                                        &
    4952!        ONLY:  urban_surface
    50    
     53
    5154    USE indices
    52    
     55
    5356    USE kinds
    5457
    5558    USE urban_surface_mod
    5659
    57     USE surface_mod   
    58    
     60    USE surface_mod
     61
    5962    USE user
    6063
  • palm/trunk/SOURCE/user_lpm_advec.f90

    r4360 r4497  
    11!> @file user_lpm_advec.f90
    2 !------------------------------------------------------------------------------!
     2!--------------------------------------------------------------------------------------------------!
    33! This file is part of the PALM model system.
    44!
    5 ! PALM is free software: you can redistribute it and/or modify it under the
    6 ! terms of the GNU General Public License as published by the Free Software
    7 ! Foundation, either version 3 of the License, or (at your option) any later
    8 ! version.
     5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
     6! Public License as published by the Free Software Foundation, either version 3 of the License, or
     7! (at your option) any later version.
    98!
    10 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
    11 ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
    12 ! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
     9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
     10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
     11! Public License for more details.
    1312!
    14 ! You should have received a copy of the GNU General Public License along with
    15 ! PALM. If not, see <http://www.gnu.org/licenses/>.
     13! You should have received a copy of the GNU General Public License along with PALM. If not, see
     14! <http://www.gnu.org/licenses/>.
    1615!
    1716! Copyright 1997-2020 Leibniz Universitaet Hannover
    18 !------------------------------------------------------------------------------!
     17!--------------------------------------------------------------------------------------------------!
     18!
    1919!
    2020! Current revisions:
     
    2525! -----------------
    2626! $Id$
     27! file re-formatted to follow the PALM coding standard
     28!
     29!
     30! 4360 2020-01-07 11:25:50Z suehring
    2731! Corrected "Former revisions" section
    28 ! 
     32!
    2933! 3768 2019-02-27 14:35:58Z raasch
    3034! variables commented + statement added to avoid compiler warnings about unused variables
    31 ! 
     35!
    3236! 3655 2019-01-07 16:51:22Z knoop
    3337! Corrected "Former revisions" section
     
    3943! ------------
    4044!> Modification of initial particles by the user.
    41 !------------------------------------------------------------------------------!
     45!--------------------------------------------------------------------------------------------------!
    4246 SUBROUTINE user_lpm_advec( ip, jp, kp )
    43  
     47
    4448
    4549    USE kinds
    46    
     50
    4751    USE particle_attributes
    48    
     52
    4953    USE user
    5054
    5155    IMPLICIT NONE
    5256
    53     INTEGER(iwp) ::  ip   !< index of particle grid box, x-direction
    54     INTEGER(iwp) ::  jp   !< index of particle grid box, y-direction
    55     INTEGER(iwp) ::  kp   !< index of particle grid box, z-direction
    56 !    INTEGER(iwp) ::  n    !< particle index
    57 !    INTEGER(iwp) ::  nb   !< index of sub-box particles are sorted in
     57    INTEGER(iwp) ::  ip  !< index of particle grid box, x-direction
     58    INTEGER(iwp) ::  jp  !< index of particle grid box, y-direction
     59    INTEGER(iwp) ::  kp  !< index of particle grid box, z-direction
     60!    INTEGER(iwp) ::  n   !< particle index
     61!    INTEGER(iwp) ::  nb  !< index of sub-box particles are sorted in
    5862
    5963!    INTEGER(iwp), DIMENSION(0:7)  ::  start_index !< start particle index for current sub-box
     
    7680!   DO  nb = 0, 7
    7781!      DO  n = start_index(nb), end_index(nb)
    78 !         particles(n)%xxx = 
     82!         particles(n)%xxx =
    7983!      ENDDO
    8084!   ENDDO
  • palm/trunk/SOURCE/user_lpm_init.f90

    r4360 r4497  
    11!> @file user_lpm_init.f90
    2 !------------------------------------------------------------------------------!
     2!--------------------------------------------------------------------------------------------------!
    33! This file is part of the PALM model system.
    44!
    5 ! PALM is free software: you can redistribute it and/or modify it under the
    6 ! terms of the GNU General Public License as published by the Free Software
    7 ! Foundation, either version 3 of the License, or (at your option) any later
    8 ! version.
     5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
     6! Public License as published by the Free Software Foundation, either version 3 of the License, or
     7! (at your option) any later version.
    98!
    10 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
    11 ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
    12 ! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
     9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
     10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
     11! Public License for more details.
    1312!
    14 ! You should have received a copy of the GNU General Public License along with
    15 ! PALM. If not, see <http://www.gnu.org/licenses/>.
     13! You should have received a copy of the GNU General Public License along with PALM. If not, see
     14! <http://www.gnu.org/licenses/>.
    1615!
    1716! Copyright 1997-2020 Leibniz Universitaet Hannover
    18 !------------------------------------------------------------------------------!
     17!--------------------------------------------------------------------------------------------------!
     18!
    1919!
    2020! Current revisions:
     
    2525! -----------------
    2626! $Id$
     27! file re-formatted to follow the PALM coding standard
     28!
     29!
     30! 4360 2020-01-07 11:25:50Z suehring
    2731! Corrected "Former revisions" section
    28 ! 
     32!
    2933! 3768 2019-02-27 14:35:58Z raasch
    3034! unused variables commented out to avoid compiler warnings
    31 ! 
     35!
    3236! 3655 2019-01-07 16:51:22Z knoop
    3337! Corrected "Former revisions" section
     
    3943! ------------
    4044!> Modification of initial particles by the user.
    41 !------------------------------------------------------------------------------!
     45!--------------------------------------------------------------------------------------------------!
    4246 SUBROUTINE user_lpm_init
    43  
     47
    4448    USE indices
    4549
    4650    USE kinds
    47    
     51
    4852    USE particle_attributes
    49    
     53
    5054    USE user
    5155
    5256    IMPLICIT NONE
    5357
    54 !    INTEGER(iwp) ::  ip   !<
    55 !    INTEGER(iwp) ::  jp   !<
    56 !    INTEGER(iwp) ::  kp   !<
    57 !    INTEGER(iwp) ::  n    !<
     58!    INTEGER(iwp) ::  ip  !<
     59!    INTEGER(iwp) ::  jp  !<
     60!    INTEGER(iwp) ::  kp  !<
     61!    INTEGER(iwp) ::  n   !<
    5862
    5963!
  • palm/trunk/SOURCE/user_module.f90

    r4495 r4497  
    11!> @file user_module.f90
    2 !------------------------------------------------------------------------------!
     2!--------------------------------------------------------------------------------------------------!
    33! This file is part of the PALM model system.
    44!
    5 ! PALM is free software: you can redistribute it and/or modify it under the
    6 ! terms of the GNU General Public License as published by the Free Software
    7 ! Foundation, either version 3 of the License, or (at your option) any later
    8 ! version.
    9 !
    10 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
    11 ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
    12 ! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
    13 !
    14 ! You should have received a copy of the GNU General Public License along with
    15 ! PALM. If not, see <http://www.gnu.org/licenses/>.
     5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
     6! Public License as published by the Free Software Foundation, either version 3 of the License, or
     7! (at your option) any later version.
     8!
     9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
     10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
     11! Public License for more details.
     12!
     13! You should have received a copy of the GNU General Public License along with PALM. If not, see
     14! <http://www.gnu.org/licenses/>.
    1615!
    1716! Copyright 1997-2020 Leibniz Universitaet Hannover
    18 !------------------------------------------------------------------------------!
     17!--------------------------------------------------------------------------------------------------!
     18!
    1919!
    2020! Current revisions:
     
    2525! -----------------
    2626! $Id$
     27! file re-formatted to follow the PALM coding standard
     28!
     29!
     30! 4495 2020-04-13 20:11:20Z raasch
    2731! restart data handling with MPI-IO added
    28 ! 
     32!
    2933! 4360 2020-01-07 11:25:50Z suehring
    30 ! Introduction of wall_flags_total_0, which currently sets bits based on static
    31 ! topography information used in wall_flags_static_0
    32 ! 
     34! Introduction of wall_flags_total_0, which currently sets bits based on static topography
     35! information used in wall_flags_static_0
     36!
    3337! 4329 2019-12-10 15:46:36Z motisi
    3438! Renamed wall_flags_0 to wall_flags_static_0
    35 ! 
     39!
    3640! 4287 2019-11-01 14:50:20Z raasch
    3741! reading of namelist file and actions in case of namelist errors revised so that statement labels
    38 ! and goto statements are not required any more; this revision also removes a previous bug
    39 ! which appeared when the namelist has been commented out in the namelist file
    40 ! 
     42! and goto statements are not required any more; this revision also removes a previous bug which
     43! appeared when the namelist has been commented out in the namelist file
     44!
    4145! 4182 2019-08-22 15:20:23Z scharf
    4246! Corrected "Former revisions" section
    43 ! 
     47!
    4448! 3986 2019-05-20 14:08:14Z Giersch
    4549! Redundant integration of control parameters in user_rrd_global removed
    46 ! 
     50!
    4751! 3911 2019-04-17 12:26:19Z knoop
    4852! Bugfix: added before_prognostic_equations case in user_actions
    49 ! 
     53!
    5054! 3768 2019-02-27 14:35:58Z raasch
    5155! variables commented + statements added to avoid compiler warnings about unused variables
     
    5357! 3767 2019-02-27 08:18:02Z raasch
    5458! unused variable for file index removed from rrd-subroutines parameter list
    55 ! 
     59!
    5660! 3747 2019-02-16 15:15:23Z gronemeier
    5761! Add routine user_init_arrays
    58 ! 
     62!
    5963! 3703 2019-01-29 16:43:53Z knoop
    6064! An example for a user defined global variable has been added (Giersch)
     
    6670! Description:
    6771! ------------
    68 !> Declaration of user-defined variables. This module may only be used
    69 !> in the user-defined routines (contained in user_interface.f90).
    70 !------------------------------------------------------------------------------!
     72!> Declaration of user-defined variables. This module may only be used in the user-defined routines
     73!> (contained in user_interface.f90).
     74!--------------------------------------------------------------------------------------------------!
    7175 MODULE user
    7276
    73 
    7477    USE arrays_3d
    7578
     
    9093    IMPLICIT NONE
    9194
    92     INTEGER(iwp) ::  dots_num_palm   !<
    93     INTEGER(iwp) ::  dots_num_user = 0  !< 
    94     INTEGER(iwp) ::  user_idummy     !<
    95    
    96     LOGICAL ::  user_module_enabled = .FALSE.   !<
    97    
    98     REAL(wp) ::  user_rdummy   !<
     95    INTEGER(iwp) ::  dots_num_palm      !<
     96    INTEGER(iwp) ::  dots_num_user = 0  !<
     97    INTEGER(iwp) ::  user_idummy        !<
     98
     99    LOGICAL ::  user_module_enabled = .FALSE.  !<
     100
     101    REAL(wp) ::  user_rdummy  !<
    99102
    100103!
    101104!-- Sample for user-defined output
    102 !    REAL(wp) :: global_parameter !< user defined global parameter
    103 !
    104 !    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  u2       !< user defined array
    105 !    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  u2_av    !< user defined array
    106 !    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  ustvst   !< user defined array
     105!    REAL(wp) :: global_parameter  !< user defined global parameter
     106!
     107!    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  u2      !< user defined array
     108!    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  u2_av   !< user defined array
     109!    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  ustvst  !< user defined array
    107110
    108111    SAVE
     
    112115!
    113116!- Public functions
    114     PUBLIC &
    115        user_parin, &
    116        user_check_parameters, &
    117        user_check_data_output_ts, &
    118        user_check_data_output_pr, &
    119        user_check_data_output, &
    120        user_define_netcdf_grid, &
    121        user_init, &
    122        user_init_arrays, &
    123        user_header, &
    124        user_actions, &
    125        user_3d_data_averaging, &
    126        user_data_output_2d, &
    127        user_data_output_3d, &
    128        user_statistics, &
    129        user_rrd_global, &
    130        user_rrd_local, &
    131        user_wrd_global, &
    132        user_wrd_local, &
    133        user_last_actions
     117    PUBLIC                                                                                         &
     118       user_actions,                                                                               &
     119       user_check_data_output,                                                                     &
     120       user_check_data_output_pr,                                                                  &
     121       user_check_data_output_ts,                                                                  &
     122       user_check_parameters,                                                                      &
     123       user_data_output_2d,                                                                        &
     124       user_data_output_3d,                                                                        &
     125       user_define_netcdf_grid,                                                                    &
     126       user_header,                                                                                &
     127       user_init,                                                                                  &
     128       user_init_arrays,                                                                           &
     129       user_last_actions,                                                                          &
     130       user_parin,                                                                                 &
     131       user_rrd_global,                                                                            &
     132       user_rrd_local,                                                                             &
     133       user_statistics,                                                                            &
     134       user_3d_data_averaging,                                                                     &
     135       user_wrd_global,                                                                            &
     136       user_wrd_local
     137
    134138
    135139!
    136140!- Public parameters, constants and initial values
    137    PUBLIC &
     141   PUBLIC                                                                                          &
    138142      user_module_enabled
    139143
     
    220224
    221225
    222 !------------------------------------------------------------------------------!
     226!--------------------------------------------------------------------------------------------------!
    223227! Description:
    224228! ------------
    225229!> Parin for &user_parameters for user module
    226 !------------------------------------------------------------------------------!
     230!--------------------------------------------------------------------------------------------------!
    227231 SUBROUTINE user_parin
    228232
    229     CHARACTER (LEN=80) ::  line        !< string containing the last line read from namelist file
    230 
    231     INTEGER(iwp) ::  i                 !<
    232     INTEGER(iwp) ::  io_status         !< status after reading the namelist file
    233     INTEGER(iwp) ::  j                 !<
     233    CHARACTER (LEN=80) ::  line  !< string containing the last line read from namelist file
     234
     235    INTEGER(iwp) ::  i          !<
     236    INTEGER(iwp) ::  io_status  !< status after reading the namelist file
     237    INTEGER(iwp) ::  j          !<
    234238
    235239
     
    247251
    248252!
    249 !-- Set revision number of this default interface version. It will be checked within
    250 !-- the main program (palm). Please change the revision number in case that the
    251 !-- current revision does not match with previous revisions (e.g. if routines
    252 !-- have been added/deleted or if parameter lists in subroutines have been changed).
     253!-- Set revision number of this default interface version. It will be checked within the main
     254!-- program (palm). Please change the revision number in case that the current revision does not
     255!-- match with previous revisions (e.g. if routines have been added/deleted or if parameter lists
     256!-- in subroutines have been changed).
    253257    user_interface_current_revision = 'r4495'
    254258
    255259!
    256 !-- Position the namelist-file at the beginning (it was already opened in
    257 !-- parin), and try to read (find) a namelist named "user_parameters".
     260!-- Position the namelist-file at the beginning (it has already been opened in parin), and try to
     261!-- read (find) a namelist named "user_parameters".
    258262    REWIND ( 11 )
    259263    READ( 11, user_parameters, IOSTAT=io_status )
     
    278282
    279283!
    280 !-- Determine the number of user-defined profiles and append them to the
    281 !-- standard data output (data_output_pr)
     284!-- Determine the number of user-defined profiles and append them to the standard data output
     285!-- (data_output_pr)
    282286    IF ( user_module_enabled )  THEN
    283287       IF ( data_output_pr_user(1) /= ' ' )  THEN
     
    300304
    301305
    302 !------------------------------------------------------------------------------!
     306!--------------------------------------------------------------------------------------------------!
    303307! Description:
    304308! ------------
    305309!> Check &userpar control parameters and deduce further quantities.
    306 !------------------------------------------------------------------------------!
     310!--------------------------------------------------------------------------------------------------!
    307311 SUBROUTINE user_check_parameters
    308312
    309 
    310 !-- Here the user may add code to check the validity of further &userpar
    311 !-- control parameters or deduce further quantities.
     313!
     314!-- Here the user may add code to check the validity of further &userpar control parameters or
     315!-- deduce further quantities.
    312316
    313317
     
    315319
    316320
    317 !------------------------------------------------------------------------------!
     321!--------------------------------------------------------------------------------------------------!
    318322! Description:
    319323! ------------
    320324!> Set module-specific timeseries units and labels
    321 !------------------------------------------------------------------------------!
     325!--------------------------------------------------------------------------------------------------!
    322326 SUBROUTINE user_check_data_output_ts( dots_max, dots_num, dots_label, dots_unit )
    323327
    324 
    325     INTEGER(iwp),      INTENT(IN)     ::  dots_max
    326     INTEGER(iwp),      INTENT(INOUT)  ::  dots_num
    327     CHARACTER (LEN=*), DIMENSION(dots_max), INTENT(INOUT)  :: dots_label
    328     CHARACTER (LEN=*), DIMENSION(dots_max), INTENT(INOUT)  :: dots_unit
     328    INTEGER(iwp),      INTENT(IN)     ::  dots_max  !<
     329    INTEGER(iwp),      INTENT(INOUT)  ::  dots_num  !<
     330
     331    CHARACTER(LEN=*), DIMENSION(dots_max), INTENT(INOUT)  ::  dots_label  !<
     332    CHARACTER(LEN=*), DIMENSION(dots_max), INTENT(INOUT)  ::  dots_unit   !<
    329333
    330334!
     
    333337
    334338!
    335 !-- Sample for user-defined time series
    336 !-- For each time series quantity you have to give a label and a unit,
    337 !-- which will be used for the NetCDF file. They must not contain more than
    338 !-- seven characters. The value of dots_num has to be increased by the
    339 !-- number of new time series quantities. Its old value has to be store in
    340 !-- dots_num_palm. See routine user_statistics on how to output calculate
    341 !-- and output these quantities.
     339!-- Sample for user-defined time series:
     340!-- For each time series quantity you have to give a label and a unit, which will be used for the
     341!-- NetCDF file. They must not contain more than seven characters. The value of dots_num has to be
     342!-- increased by the number of new time series quantities. Its old value has to be stored in
     343!-- dots_num_palm. See routine user_statistics on how to calculate and output these quantities.
    342344
    343345!    dots_num_palm = dots_num
     
    357359
    358360
    359 !------------------------------------------------------------------------------!
    360 ! Description:
    361 ! ------------
    362 !> Set the unit of user defined profile output quantities. For those variables
    363 !> not recognized by the user, the parameter unit is set to "illegal", which
    364 !> tells the calling routine that the output variable is not defined and leads
    365 !> to a program abort.
    366 !------------------------------------------------------------------------------!
     361!--------------------------------------------------------------------------------------------------!
     362! Description:
     363! ------------
     364!> Set the unit of user defined profile output quantities. For those variables not recognized by the
     365!> user, the parameter unit is set to "illegal", which tells the calling routine that the
     366!> output variable is not defined and leads to a program abort.
     367!--------------------------------------------------------------------------------------------------!
    367368 SUBROUTINE user_check_data_output_pr( variable, var_count, unit, dopr_unit )
    368369
     
    371372
    372373
    373     CHARACTER (LEN=*) ::  unit     !<
    374     CHARACTER (LEN=*) ::  variable !<
     374    CHARACTER (LEN=*) ::  unit      !<
     375    CHARACTER (LEN=*) ::  variable  !<
    375376    CHARACTER (LEN=*) ::  dopr_unit !< local value of dopr_unit
    376377
    377 !    INTEGER(iwp) ::  user_pr_index !<
    378     INTEGER(iwp) ::  var_count     !<
     378!    INTEGER(iwp) ::  user_pr_index  !<
     379    INTEGER(iwp) ::  var_count      !<
    379380
    380381!
     
    386387!
    387388!--    Uncomment and extend the following lines, if necessary.
    388 !--    Add additional CASE statements depending on the number of quantities
    389 !--    for which profiles are to be calculated. The respective calculations
    390 !--    to be performed have to be added in routine user_statistics.
    391 !--    The quantities are (internally) identified by a user-profile-number
    392 !--    (see variable "user_pr_index" below). The first user-profile must be assigned
    393 !--    the number "pr_palm+1", the second one "pr_palm+2", etc. The respective
    394 !--    user-profile-numbers have also to be used in routine user_statistics!
    395 !       CASE ( 'u*v*' )                      ! quantity string as given in
    396 !                                            ! data_output_pr_user
     389!--    Add additional CASE statements depending on the number of quantities for which profiles are
     390!--    to be calculated. The respective calculations to be performed have to be added in routine
     391!--    user_statistics. The quantities are (internally) identified by a user-profile-number
     392!--    (see variable "user_pr_index" below). The first user-profile must be assigned the number
     393!--    "pr_palm+1", the second one "pr_palm+2", etc. The respective user-profile-numbers have also
     394!--    to be used in routine user_statistics!
     395!       CASE ( 'u*v*' )                      ! quantity string as given in data_output_pr_user
    397396!          user_pr_index = pr_palm + 1
    398397!          dopr_index(var_count)  = user_pr_index    ! quantities' user-profile-number
    399398!          dopr_unit = 'm2/s2'  ! quantity unit
    400399!          unit = dopr_unit
    401 !          hom(:,2,user_pr_index,:)       = SPREAD( zu, 2, statistic_regions+1 )
    402 !                                            ! grid on which the quantity is
    403 !                                            ! defined (use zu or zw)
     400!          hom(:,2,user_pr_index,:) = SPREAD( zu, 2, statistic_regions+1 )
     401!                                            ! grid on which the quantity is defined (use zu or zw)
     402!
    404403
    405404       CASE DEFAULT
     
    412411
    413412
    414 !------------------------------------------------------------------------------!
    415 ! Description:
    416 ! ------------
    417 !> Set the unit of user defined output quantities. For those variables
    418 !> not recognized by the user, the parameter unit is set to "illegal", which
    419 !> tells the calling routine that the output variable is not defined and leads
    420 !> to a program abort.
    421 !------------------------------------------------------------------------------!
     413!--------------------------------------------------------------------------------------------------!
     414! Description:
     415! ------------
     416!> Set the unit of user defined output quantities. For those variables not recognized by the user,
     417!> the parameter unit is set to "illegal", which tells the calling routine that the output variable
     418!> is not defined and leads to a program abort.
     419!--------------------------------------------------------------------------------------------------!
    422420 SUBROUTINE user_check_data_output( variable, unit )
    423421
    424422
    425     CHARACTER (LEN=*) ::  unit     !<
    426     CHARACTER (LEN=*) ::  variable !<
     423    CHARACTER (LEN=*) ::  unit      !<
     424    CHARACTER (LEN=*) ::  variable  !<
    427425
    428426
     
    446444
    447445
    448 !------------------------------------------------------------------------------!
     446!--------------------------------------------------------------------------------------------------!
    449447! Description:
    450448! ------------
    451449!> Initialize user-defined arrays
    452 !------------------------------------------------------------------------------!
     450!--------------------------------------------------------------------------------------------------!
    453451 SUBROUTINE user_init_arrays
    454452
     
    468466!     IF ( statistic_regions >= 1 )  THEN
    469467!        region = 1
    470 ! 
     468!
    471469!        rmask(:,:,region) = 0.0_wp
    472470!        DO  i = nxl, nxr
     
    479477!           ENDIF
    480478!        ENDDO
    481 ! 
     479!
    482480!     ENDIF
    483481
     
    485483
    486484
    487 !------------------------------------------------------------------------------!
     485!--------------------------------------------------------------------------------------------------!
    488486! Description:
    489487! ------------
    490488!> Execution of user-defined initializing actions
    491 !------------------------------------------------------------------------------!
     489!--------------------------------------------------------------------------------------------------!
    492490 SUBROUTINE user_init
    493491
    494492
    495 !    CHARACTER (LEN=20) :: field_char   !<
     493!    CHARACTER(LEN=20) :: field_char  !<
    496494!
    497495!-- Here the user-defined initializing actions follow:
     
    503501
    504502
    505 !------------------------------------------------------------------------------!
    506 ! Description:
    507 ! ------------
    508 !> Set the grids on which user-defined output quantities are defined.
    509 !> Allowed values for grid_x are "x" and "xu", for grid_y "y" and "yv", and
    510 !> for grid_z "zu" and "zw".
    511 !------------------------------------------------------------------------------!
     503!--------------------------------------------------------------------------------------------------!
     504! Description:
     505! ------------
     506!> Set the grids on which user-defined output quantities are defined. Allowed values for grid_x are
     507!> "x" and "xu", for grid_y "y" and "yv", and for grid_z "zu" and "zw".
     508!--------------------------------------------------------------------------------------------------!
    512509 SUBROUTINE user_define_netcdf_grid( variable, found, grid_x, grid_y, grid_z )
    513510
     
    551548
    552549
    553 !------------------------------------------------------------------------------!
     550!--------------------------------------------------------------------------------------------------!
    554551! Description:
    555552! ------------
    556553!> Print a header with user-defined information.
    557 !------------------------------------------------------------------------------!
     554!--------------------------------------------------------------------------------------------------!
    558555 SUBROUTINE user_header( io )
    559556
    560557
    561     INTEGER(iwp) ::  i    !<
    562     INTEGER(iwp) ::  io   !<
    563 
    564 !
    565 !-- If no user-defined variables are read from the namelist-file, no
    566 !-- information will be printed.
     558    INTEGER(iwp) ::  i   !<
     559    INTEGER(iwp) ::  io  !<
     560
     561!
     562!-- If no user-defined variables are read from the namelist-file, no information will be printed.
    567563    IF ( .NOT. user_module_enabled )  THEN
    568564       WRITE ( io, 100 )
     
    584580!-- Format-descriptors
    585581100 FORMAT (//' *** no user-defined variables found'/)
    586 110 FORMAT (//1X,78('#')                                                       &
    587             //' User-defined variables and actions:'/                          &
    588               ' -----------------------------------'//)
     582110 FORMAT (//1X,78('#') // ' User-defined variables and actions:' /                               &
     583            ' -----------------------------------'//)
    589584200 FORMAT (' Output of profiles and time series for following regions:' /)
    590585201 FORMAT (4X,'Region ',I1,':   ',A)
     
    594589
    595590
    596 !------------------------------------------------------------------------------!
     591!--------------------------------------------------------------------------------------------------!
    597592! Description:
    598593! ------------
    599594!> Call for all grid points
    600 !------------------------------------------------------------------------------!
     595!--------------------------------------------------------------------------------------------------!
    601596 SUBROUTINE user_actions( location )
    602597
    603598
    604     CHARACTER (LEN=*) ::  location !<
    605 
    606 !    INTEGER(iwp) ::  i !<
    607 !    INTEGER(iwp) ::  j !<
    608 !    INTEGER(iwp) ::  k !<
     599    CHARACTER(LEN=*) ::  location  !<
     600
     601!    INTEGER(iwp) ::  i  !<
     602!    INTEGER(iwp) ::  j  !<
     603!    INTEGER(iwp) ::  k  !<
    609604
    610605    CALL cpu_log( log_point(24), 'user_actions', 'start' )
    611606
    612607!
    613 !-- Here the user-defined actions follow
    614 !-- No calls for single grid points are allowed at locations before and
    615 !-- after the timestep, since these calls are not within an i,j-loop
     608!-- Here the user-defined actions follow. No calls for single grid points are allowed at locations
     609!-- before and after the timestep, since these calls are not within an i,j-loop
    616610    SELECT CASE ( location )
    617611
     
    626620       CASE ( 'after_integration' )
    627621!
    628 !--       Enter actions to be done after every time integration (before
    629 !--       data output)
     622!--       Enter actions to be done after every time integration (before data output)
    630623!--       Sample for user-defined output:
    631624!          DO  i = nxlg, nxrg
     
    640633!                DO  k = nzb, nzt+1
    641634!                   ustvst(k,j,i) =  &
    642 !                      ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) - hom(k,1,1,0) ) * &
     635!                      ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) - hom(k,1,1,0) ) *                      &
    643636!                      ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) - hom(k,1,2,0) )
    644637!                ENDDO
     
    688681
    689682
    690 !------------------------------------------------------------------------------!
     683!--------------------------------------------------------------------------------------------------!
    691684! Description:
    692685! ------------
    693686!> Call for grid point i,j
    694 !------------------------------------------------------------------------------!
     687!--------------------------------------------------------------------------------------------------!
    695688 SUBROUTINE user_actions_ij( i, j, location )
    696689
    697690
    698     CHARACTER (LEN=*) ::  location
    699 
    700     INTEGER(iwp) ::  i
    701     INTEGER(iwp) ::  j
     691    CHARACTER(LEN=*) ::  location  !<
     692
     693    INTEGER(iwp) ::  i  !<
     694    INTEGER(iwp) ::  j  !<
    702695
    703696!
     
    744737
    745738
    746 !------------------------------------------------------------------------------!
    747 ! Description:
    748 ! ------------
    749 !> Sum up and time-average user-defined output quantities as well as allocate
    750 !> the array necessary for storing the average.
    751 !------------------------------------------------------------------------------!
     739!--------------------------------------------------------------------------------------------------!
     740! Description:
     741! ------------
     742!> Sum up and time-average user-defined output quantities as well as allocate the array necessary
     743!> for storing the average.
     744!--------------------------------------------------------------------------------------------------!
    752745 SUBROUTINE user_3d_data_averaging( mode, variable )
    753746
    754747
    755     CHARACTER (LEN=*) ::  mode    !<
    756     CHARACTER (LEN=*) :: variable !<
    757 
    758 !    INTEGER(iwp) ::  i !<
    759 !    INTEGER(iwp) ::  j !<
    760 !    INTEGER(iwp) ::  k !<
     748    CHARACTER(LEN=*) ::  mode      !<
     749    CHARACTER(LEN=*) ::  variable  !<
     750
     751!    INTEGER(iwp) ::  i  !<
     752!    INTEGER(iwp) ::  j  !<
     753!    INTEGER(iwp) ::  k  !<
    761754
    762755    IF ( mode == 'allocate' )  THEN
     
    766759!
    767760!--       Uncomment and extend the following lines, if necessary.
    768 !--       The arrays for storing the user defined quantities (here u2_av) have
    769 !--       to be declared and defined by the user!
     761!--       The arrays for storing the user defined quantities (here u2_av) have to be declared and
     762!--       defined by the user!
    770763!--       Sample for user-defined output:
    771764!          CASE ( 'u2' )
     
    786779!
    787780!--       Uncomment and extend the following lines, if necessary.
    788 !--       The arrays for storing the user defined quantities (here u2 and
    789 !--       u2_av) have to be declared and defined by the user!
     781!--       The arrays for storing the user defined quantities (here u2 and u2_av) have to be declared
     782!--       and defined by the user!
    790783!--       Sample for user-defined output:
    791784!          CASE ( 'u2' )
    792 !             IF ( ALLOCATED( u2_av ) ) THEN
     785!             IF ( ALLOCATED( u2_av ) )  THEN
    793786!                DO  i = nxlg, nxrg
    794787!                   DO  j = nysg, nyng
     
    811804!
    812805!--       Uncomment and extend the following lines, if necessary.
    813 !--       The arrays for storing the user defined quantities (here u2_av) have
    814 !--       to be declared and defined by the user!
     806!--       The arrays for storing the user defined quantities (here u2_av) have to be declared and
     807!--       defined by the user!
    815808!--       Sample for user-defined output:
    816809!          CASE ( 'u2' )
    817 !             IF ( ALLOCATED( u2_av ) ) THEN
     810!             IF ( ALLOCATED( u2_av ) )  THEN
    818811!                DO  i = nxlg, nxrg
    819812!                   DO  j = nysg, nyng
     
    833826
    834827
    835 !------------------------------------------------------------------------------!
    836 ! Description:
    837 ! ------------
    838 !> Resorts the user-defined output quantity with indices (k,j,i) to a
    839 !> temporary array with indices (i,j,k) and sets the grid on which it is defined.
    840 !> Allowed values for grid are "zu" and "zw".
    841 !------------------------------------------------------------------------------!
     828!--------------------------------------------------------------------------------------------------!
     829! Description:
     830! ------------
     831!> Resorts the user-defined output quantity with indices (k,j,i) to a temporary array with indices
     832!> (i,j,k) and sets the grid on which it is defined. Allowed values for grid are "zu" and "zw".
     833!--------------------------------------------------------------------------------------------------!
    842834 SUBROUTINE user_data_output_2d( av, variable, found, grid, local_pf, two_d, nzb_do, nzt_do )
    843835
    844836
    845     CHARACTER (LEN=*) ::  grid     !<
    846     CHARACTER (LEN=*) ::  variable !<
    847 
    848     INTEGER(iwp) ::  av     !< flag to control data output of instantaneous or time-averaged data
    849 !    INTEGER(iwp) ::  i      !< grid index along x-direction
    850 !    INTEGER(iwp) ::  j      !< grid index along y-direction
    851 !    INTEGER(iwp) ::  k      !< grid index along z-direction
    852 !    INTEGER(iwp) ::  m      !< running index surface elements
    853     INTEGER(iwp) ::  nzb_do !< lower limit of the domain (usually nzb)
    854     INTEGER(iwp) ::  nzt_do !< upper limit of the domain (usually nzt+1)
    855 
    856     LOGICAL      ::  found !<
    857     LOGICAL      ::  two_d !< flag parameter that indicates 2D variables (horizontal cross sections)
     837    CHARACTER(LEN=*) ::  grid      !<
     838    CHARACTER(LEN=*) ::  variable  !<
     839
     840    INTEGER(iwp) ::  av      !< flag to control data output of instantaneous or time-averaged data
     841!    INTEGER(iwp) ::  i       !< grid index along x-direction
     842!    INTEGER(iwp) ::  j       !< grid index along y-direction
     843!    INTEGER(iwp) ::  k       !< grid index along z-direction
     844!    INTEGER(iwp) ::  m       !< running index surface elements
     845    INTEGER(iwp) ::  nzb_do  !< lower limit of the domain (usually nzb)
     846    INTEGER(iwp) ::  nzt_do  !< upper limit of the domain (usually nzt+1)
     847
     848    LOGICAL      ::  found  !<
     849    LOGICAL      ::  two_d  !< flag parameter that indicates 2D variables (horizontal cross sections)
    858850
    859851!    REAL(wp) ::  fill_value = -999.0_wp    !< value for the _FillValue attribute
    860852
    861     REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !<
     853    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf  !<
    862854
    863855!
     
    872864!
    873865!--    Uncomment and extend the following lines, if necessary.
    874 !--    The arrays for storing the user defined quantities (here u2 and u2_av)
    875 !--    have to be declared and defined by the user!
     866!--    The arrays for storing the user defined quantities (here u2 and u2_av) have to be declared
     867!--    and defined by the user!
    876868!--    Sample for user-defined output:
    877869!       CASE ( 'u2_xy', 'u2_xz', 'u2_yz' )
     
    885877!             ENDDO
    886878!          ELSE
    887 !             IF ( .NOT. ALLOCATED( u2_av ) ) THEN
     879!             IF ( .NOT. ALLOCATED( u2_av ) )  THEN
    888880!                ALLOCATE( u2_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    889881!                u2_av = REAL( fill_value, KIND = wp )
     
    900892!          grid = 'zu'
    901893!
    902 !--    In case two-dimensional surface variables are output, the user
    903 !--    has to access related surface-type. Uncomment and extend following lines
    904 !--    appropriately (example output of vertical surface momentum flux of u-
    905 !--    component). Please note, surface elements can be distributed over
    906 !--    several data type, depending on their respective surface properties.
     894!--    In case two-dimensional surface variables are output, the user has to access related
     895!--    surface-type. Uncomment and extend following lines appropriately (example output of vertical
     896!--    surface momentum flux of u-component). Please note, surface elements can be distributed over
     897!--    several data types, depending on their respective surface properties.
    907898!       CASE ( 'usws_xy' )
    908899!          IF ( av == 0 )  THEN
     
    931922!
    932923!          grid = 'zu'
    933 !--       
     924!--
    934925
    935926
     
    944935
    945936
    946 !------------------------------------------------------------------------------!
    947 ! Description:
    948 ! ------------
    949 !> Resorts the user-defined output quantity with indices (k,j,i) to a
    950 !> temporary array with indices (i,j,k).
    951 !------------------------------------------------------------------------------!
     937!--------------------------------------------------------------------------------------------------!
     938! Description:
     939! ------------
     940!> Resorts the user-defined output quantity with indices (k,j,i) to a temporary array with indices
     941!> (i,j,k).
     942!--------------------------------------------------------------------------------------------------!
    952943 SUBROUTINE user_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
    953944
    954945
    955     CHARACTER (LEN=*) ::  variable !<
    956 
    957     INTEGER(iwp) ::  av    !<
    958 !    INTEGER(iwp) ::  i     !<
    959 !    INTEGER(iwp) ::  j     !<
    960 !    INTEGER(iwp) ::  k     !<
     946    CHARACTER(LEN=*) ::  variable  !<
     947
     948    INTEGER(iwp) ::  av     !<
     949!    INTEGER(iwp) ::  i      !<
     950!    INTEGER(iwp) ::  j      !<
     951!    INTEGER(iwp) ::  k      !<
    961952    INTEGER(iwp) ::  nzb_do !< lower limit of the data output (usually 0)
    962953    INTEGER(iwp) ::  nzt_do !< vertical upper limit of the data output (usually nz_do3d)
    963954
    964     LOGICAL      ::  found !<
    965 
    966 !    REAL(wp) ::  fill_value = -999.0_wp    !< value for the _FillValue attribute
    967 
    968     REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !<
     955    LOGICAL      ::  found  !<
     956
     957!    REAL(wp) ::  fill_value = -999.0_wp  !< value for the _FillValue attribute
     958
     959    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf  !<
    969960
    970961!
     
    979970!
    980971!--    Uncomment and extend the following lines, if necessary.
    981 !--    The arrays for storing the user defined quantities (here u2 and u2_av)
    982 !--    have to be declared and defined by the user!
     972!--    The arrays for storing the user defined quantities (here u2 and u2_av) have to be declared
     973!--    and defined by the user!
    983974!--    Sample for user-defined output:
    984975!       CASE ( 'u2' )
     
    992983!             ENDDO
    993984!          ELSE
    994 !             IF ( .NOT. ALLOCATED( u2_av ) ) THEN
     985!             IF ( .NOT. ALLOCATED( u2_av ) )  THEN
    995986!                ALLOCATE( u2_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    996987!                u2_av = REAL( fill_value, KIND = wp )
     
    10151006
    10161007
    1017 !------------------------------------------------------------------------------!
    1018 ! Description:
    1019 ! ------------
    1020 !> Calculation of user-defined statistics, i.e. horizontally averaged profiles
    1021 !> and time series.
    1022 !> This routine is called for every statistic region sr defined by the user,
    1023 !> but at least for the region "total domain" (sr=0).
    1024 !> See section 3.5.4 on how to define, calculate, and output user defined
    1025 !> quantities.
    1026 !------------------------------------------------------------------------------!
     1008!--------------------------------------------------------------------------------------------------!
     1009! Description:
     1010! ------------
     1011!> Calculation of user-defined statistics, i.e. horizontally averaged profiles and time series.
     1012!> This routine is called for every statistic region sr defined by the user, but at least for the
     1013!> region "total domain" (sr=0). See section 3.5.4 on how to define, calculate, and output user
     1014!> defined quantities.
     1015!--------------------------------------------------------------------------------------------------!
    10271016 SUBROUTINE user_statistics( mode, sr, tn )
    10281017
    10291018
    1030     CHARACTER (LEN=*) ::  mode   !<
    1031 !    INTEGER(iwp) ::  i    !<
    1032 !    INTEGER(iwp) ::  j    !<
    1033 !    INTEGER(iwp) ::  k    !<
    1034     INTEGER(iwp) ::  sr   !<
    1035     INTEGER(iwp) ::  tn   !<
    1036 
    1037 !    REAL(wp), DIMENSION(:), ALLOCATABLE ::  ts_value_l   !<
     1019    CHARACTER(LEN=*) ::  mode  !<
     1020!    INTEGER(iwp) ::  i   !<
     1021!    INTEGER(iwp) ::  j   !<
     1022!    INTEGER(iwp) ::  k   !<
     1023    INTEGER(iwp) ::  sr  !<
     1024    INTEGER(iwp) ::  tn  !<
     1025
     1026!    REAL(wp), DIMENSION(:), ALLOCATABLE ::  ts_value_l  !<
    10381027
    10391028!
     
    10441033
    10451034!
    1046 !--    Sample on how to calculate horizontally averaged profiles of user-
    1047 !--    defined quantities. Each quantity is identified by the index
    1048 !--    "pr_palm+#" where "#" is an integer starting from 1. These
    1049 !--    user-profile-numbers must also be assigned to the respective strings
    1050 !--    given by data_output_pr_user in routine user_check_data_output_pr.
     1035!--    Sample on how to calculate horizontally averaged profiles of user-defined quantities. Each
     1036!--    quantity is identified by the index "pr_palm+#" where "#" is an integer starting from 1.
     1037!--    These user-profile-numbers must also be assigned to the respective strings given by
     1038!--    data_output_pr_user in routine user_check_data_output_pr.
    10511039!       !$OMP DO
    10521040!       DO  i = nxl, nxr
     
    10541042!             DO  k = nzb+1, nzt
    10551043!!
    1056 !!--             Sample on how to calculate the profile of the resolved-scale
    1057 !!--             horizontal momentum flux u*v*
    1058 !                sums_l(k,pr_palm+1,tn) = sums_l(k,pr_palm+1,tn) +             &
    1059 !                      ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) - hom(k,1,1,sr) ) *&
    1060 !                      ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) - hom(k,1,2,sr) )  &
    1061 !                                     * rmask(j,i,sr)                          &
    1062 !                                     * MERGE( 1.0_wp, 0.0_wp,                 &
    1063 !                                              BTEST( wall_flags_total_0(k,j,i), 0 ) )
     1044!!--             Sample on how to calculate the profile of the resolved-scale horizontal momentum
     1045!!--             flux u*v*
     1046!                sums_l(k,pr_palm+1,tn) = sums_l(k,pr_palm+1,tn) +                                  &
     1047!                                         ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) - hom(k,1,1,sr) ) *  &
     1048!                                         ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) - hom(k,1,2,sr) ) *  &
     1049!                                         rmask(j,i,sr) * MERGE( 1.0_wp, 0.0_wp,                    &
     1050!                                         BTEST( wall_flags_total_0(k,j,i), 0 ) )
    10641051!!
    1065 !!--             Further profiles can be defined and calculated by increasing
    1066 !!--             the second index of array sums_l (replace ... appropriately)
    1067 !                sums_l(k,pr_palm+2,tn) = sums_l(k,pr_palm+2,tn) + ...           &
    1068 !                                         * rmask(j,i,sr)
     1052!!--             Further profiles can be defined and calculated by increasing the second index of
     1053!!--             array sums_l (replace ... appropriately)
     1054!                sums_l(k,pr_palm+2,tn) = sums_l(k,pr_palm+2,tn) + ...   * rmask(j,i,sr)
    10691055!             ENDDO
    10701056!          ENDDO
     
    10771063!
    10781064!--    Sample on how to add values for the user-defined time series quantities.
    1079 !--    These have to be defined before in routine user_init. This sample
    1080 !--    creates two time series for the absolut values of the horizontal
    1081 !--    velocities u and v.
     1065!--    These have to be defined before in routine user_init. This sample creates two time series for
     1066!--    the absolut values of the horizontal velocities u and v.
    10821067!       ts_value_l = 0.0_wp
    10831068!       ts_value_l(1) = ABS( u_max )
     
    10851070!
    10861071!--     Collect / send values to PE0, because only PE0 outputs the time series.
    1087 !--     CAUTION: Collection is done by taking the sum over all processors.
    1088 !--              You may have to normalize this sum, depending on the quantity
    1089 !--              that you like to calculate. For serial runs, nothing has to be
    1090 !--              done.
    1091 !--     HINT: If the time series value that you are calculating has the same
    1092 !--           value on all PEs, you can omit the MPI_ALLREDUCE call and
    1093 !--           assign ts_value(dots_num_palm+1:,sr) = ts_value_l directly.
     1072!--     CAUTION: Collection is done by taking the sum over all processors. You may have to normalize
     1073!--              this sum, depending on the quantity that you like to calculate. For serial runs,
     1074!--              nothing has to be done.
     1075!--     HINT: If the time series value that you are calculating has the same value on all PEs, you
     1076!--           can omit the MPI_ALLREDUCE call and assign ts_value(dots_num_palm+1:,sr) = ts_value_l directly.
    10941077!#if defined( __parallel )
    10951078!       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    1096 !       CALL MPI_ALLREDUCE( ts_value_l(1),                         &
    1097 !                           ts_value(dots_num_palm+1,sr),                        &
    1098 !                           dots_num_user, MPI_REAL, MPI_MAX, comm2d,   &
    1099 !                           ierr )
     1079!       CALL MPI_ALLREDUCE( ts_value_l(1), ts_value(dots_num_palm+1,sr), dots_num_user, MPI_REAL,   &
     1080!                           MPI_MAX, comm2d, ierr )
    11001081!#else
    11011082!       ts_value(dots_num_palm+1:dots_num_palm+dots_num_user,sr) = ts_value_l
     
    11071088
    11081089
    1109 !------------------------------------------------------------------------------!
     1090!--------------------------------------------------------------------------------------------------!
    11101091! Description:
    11111092! ------------
    11121093!> Read module-specific global restart data (Fortran binary format).
    1113 !------------------------------------------------------------------------------!
     1094!--------------------------------------------------------------------------------------------------!
    11141095 SUBROUTINE user_rrd_global_ftn( found )
    11151096
    11161097
    1117     LOGICAL, INTENT(OUT)  ::  found
     1098    LOGICAL, INTENT(OUT)  ::  found  !<
    11181099
    11191100
     
    11271108
    11281109       CASE DEFAULT
    1129  
     1110
    11301111          found = .FALSE.
    11311112
     
    11361117
    11371118
    1138 !------------------------------------------------------------------------------!
     1119!--------------------------------------------------------------------------------------------------!
    11391120! Description:
    11401121! ------------
    11411122!> Read module-specific global restart data (MPI-IO).
    1142 !------------------------------------------------------------------------------!
     1123!--------------------------------------------------------------------------------------------------!
    11431124 SUBROUTINE user_rrd_global_mpi
    11441125
     
    11491130
    11501131
    1151 !------------------------------------------------------------------------------!
    1152 ! Description:
    1153 ! ------------
    1154 !> Reading processor specific restart data from file(s) that has been defined
    1155 !> by the user.
    1156 !> Subdomain index limits on file are given by nxl_on_file, etc.
    1157 !> Indices nxlc, etc. indicate the range of gridpoints to be mapped from the
    1158 !> subdomain on file (f) to the subdomain of the current PE (c). They have been
    1159 !> calculated in routine rrd_local.
    1160 !------------------------------------------------------------------------------!
    1161  SUBROUTINE user_rrd_local( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc,         &
    1162                             nxr_on_file, nynf, nync, nyn_on_file, nysf,     &
    1163                             nysc, nys_on_file, tmp_3d, found )
     1132!--------------------------------------------------------------------------------------------------!
     1133! Description:
     1134! ------------
     1135!> Reading processor specific restart data from file(s) that has been defined by the user. Subdomain
     1136!> index limits on file are given by nxl_on_file, etc. Indices nxlc, etc. indicate the range of
     1137!> gridpoints to be mapped from the subdomain on file (f) to the subdomain of the current PE (c).
     1138!> They have been calculated in routine rrd_local.
     1139!--------------------------------------------------------------------------------------------------!
     1140 SUBROUTINE user_rrd_local( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, nxr_on_file, nynf, nync,       &
     1141                            nyn_on_file, nysf, nysc, nys_on_file, tmp_3d, found )
    11641142
    11651143
     
    11791157    INTEGER(iwp) ::  nys_on_file     !<
    11801158
    1181     LOGICAL, INTENT(OUT)  ::  found
    1182 
    1183     REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d   !<
     1159    LOGICAL, INTENT(OUT)  ::  found  !<
     1160
     1161    REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d  !<
    11841162
    11851163!
     
    11971175
    11981176       CASE ( 'u2_av' )
    1199 !          IF ( .NOT. ALLOCATED( u2_av ) ) THEN
     1177!          IF ( .NOT. ALLOCATED( u2_av ) )  THEN
    12001178!               ALLOCATE( u2_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    12011179!          ENDIF
    12021180!          IF ( k == 1 )  READ ( 13 )  tmp_3d
    1203 !             u2_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =         &
    1204 !                tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     1181!             u2_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =                                    &
     1182!             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    12051183!
    12061184       CASE DEFAULT
     
    12131191
    12141192
    1215 !------------------------------------------------------------------------------!
    1216 ! Description:
    1217 ! ------------
    1218 !> Writes global and user-defined restart data into binary file(s) for restart
    1219 !> runs.
    1220 !------------------------------------------------------------------------------!
     1193!--------------------------------------------------------------------------------------------------!
     1194! Description:
     1195! ------------
     1196!> Writes global and user-defined restart data into binary file(s) for restart runs.
     1197!--------------------------------------------------------------------------------------------------!
    12211198 SUBROUTINE user_wrd_global
    12221199
     
    12351212
    12361213
    1237 !------------------------------------------------------------------------------!
    1238 ! Description:
    1239 ! ------------
    1240 !> Writes processor specific and user-defined restart data into binary file(s)
    1241 !> for restart runs.
    1242 !------------------------------------------------------------------------------!
     1214!--------------------------------------------------------------------------------------------------!
     1215! Description:
     1216! ------------
     1217!> Writes processor specific and user-defined restart data into binary file(s) for restart runs.
     1218!--------------------------------------------------------------------------------------------------!
    12431219 SUBROUTINE user_wrd_local
    12441220
     
    12621238
    12631239
    1264 !------------------------------------------------------------------------------!
     1240!--------------------------------------------------------------------------------------------------!
    12651241! Description:
    12661242! ------------
    12671243!> Execution of user-defined actions at the end of a job.
    1268 !------------------------------------------------------------------------------!
     1244!--------------------------------------------------------------------------------------------------!
    12691245 SUBROUTINE user_last_actions
    12701246
  • palm/trunk/SOURCE/user_spectra.f90

    r4360 r4497  
    11!> @file user_spectra.f90
    2 !------------------------------------------------------------------------------!
     2!--------------------------------------------------------------------------------------------------!
    33! This file is part of the PALM model system.
    44!
    5 ! PALM is free software: you can redistribute it and/or modify it under the
    6 ! terms of the GNU General Public License as published by the Free Software
    7 ! Foundation, either version 3 of the License, or (at your option) any later
    8 ! version.
     5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
     6! Public License as published by the Free Software Foundation, either version 3 of the License, or
     7! (at your option) any later version.
    98!
    10 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
    11 ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
    12 ! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
     9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
     10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
     11! Public License for more details.
    1312!
    14 ! You should have received a copy of the GNU General Public License along with
    15 ! PALM. If not, see <http://www.gnu.org/licenses/>.
     13! You should have received a copy of the GNU General Public License along with PALM. If not, see
     14! <http://www.gnu.org/licenses/>.
    1615!
    1716! Copyright 1997-2020 Leibniz Universitaet Hannover
    18 !------------------------------------------------------------------------------!
     17!--------------------------------------------------------------------------------------------------!
     18!
    1919!
    2020! Current revisions:
     
    2525! -----------------
    2626! $Id$
     27! file re-formatted to follow the PALM coding standard
     28!
     29! 4360 2020-01-07 11:25:50Z suehring
    2730! Corrected "Former revisions" section
    28 ! 
     31!
    2932! 3768 2019-02-27 14:35:58Z raasch
    3033! variables removed + statement added to avoid compiler warnings about unused variables
    31 ! 
     34!
    3235! 3655 2019-01-07 16:51:22Z knoop
    3336! Renamed output variables
     
    3942! ------------
    4043!> Calculation of user-defined spectra.
    41 !> See section 3.5.4 on how to define, calculate, and output user defined
    42 !> quantities.
    43 !------------------------------------------------------------------------------!
     44!> See section 3.5.4 on how to define, calculate, and output user defined quantities.
     45!--------------------------------------------------------------------------------------------------!
    4446 SUBROUTINE user_spectra( mode, m, pr )
    45  
     47
    4648
    4749    USE arrays_3d
    48    
     50
    4951    USE control_parameters
    50    
     52
    5153    USE indices
    52    
     54
    5355    USE kinds
    54    
     56
    5557    USE spectra_mod
    56    
     58
    5759    USE statistics
    58    
     60
    5961    USE user
    6062
    6163    IMPLICIT NONE
    6264
    63     CHARACTER (LEN=*) ::  mode
     65    CHARACTER(LEN=*) ::  mode  !<
    6466
    65     INTEGER(iwp) ::  m    !<
    66     INTEGER(iwp) ::  pr   !<
     67    INTEGER(iwp) ::  m   !<
     68    INTEGER(iwp) ::  pr  !<
    6769
    6870
     
    7274
    7375!
    74 !-- Sample on how to calculate spectra of user-defined quantities.
    75 !-- Each quantity is identified by the corresponding user profile index
    76 !-- "pr_palm+#" where "#" is an integer starting from 1. These
    77 !-- user-profile-numbers must also be assigned to the respective strings
    78 !-- given by data_output_pr_user in routine user_check_data_output_pr.
     76!-- Sample on how to calculate spectra of user-defined quantities. Each quantity is identified by
     77!-- the corresponding user profile index "pr_palm+#" where "#" is an integer starting from 1. These
     78!-- user-profile-numbers must also be assigned to the respective strings given by
     79!-- data_output_pr_user in routine user_check_data_output_pr.
    7980    IF ( mode == 'preprocess' )  THEN
    8081
    8182       SELECT CASE ( TRIM( data_output_sp(m) ) )
    82          
     83
    8384          CASE ( 'u', 'v', 'w', 'theta', 'q', 's' )
    84 !--          Not allowed here since these are the standard quantities used in
    85 !--          preprocess_spectra.
    86        
     85!--          Not allowed here since these are the standard quantities used in preprocess_spectra.
     86
    8787!          CASE ( 'u*v*' )
    8888!             pr = pr_palm+1
    8989!             d(nzb+1:nzt,nys:nyn,nxl:nxr) = ustvst(nzb+1:nzt,nys:nyn,nxl:nxr)
    90        
     90
    9191          CASE DEFAULT
    92              message_string = 'Spectra of ' //                                 &
    93                          TRIM( data_output_sp(m) ) // ' can not be calculated'
     92             message_string = 'Spectra of ' // TRIM( data_output_sp(m) ) // ' can not be calculated'
    9493             CALL message( 'user_spectra', 'UI0010', 0, 1, 0, 6, 0 )
    95            
     94
    9695       END SELECT
    9796
     
    101100
    102101          CASE ( 'u', 'v', 'w', 'theta', 'q', 's' )
    103 !--          Not allowed here since these are the standard quantities used in
    104 !--          data_output_spectra.
     102!--          Not allowed here since these are the standard quantities used in data_output_spectra.
    105103
    106104!          CASE ( 'u*v*' )
     
    108106
    109107          CASE DEFAULT
    110              message_string = 'Spectra of ' //                                 &
    111                               TRIM( data_output_sp(m) ) // ' are not defined'
     108             message_string = 'Spectra of ' // TRIM( data_output_sp(m) ) // ' are not defined'
    112109             CALL message( 'user_spectra', 'UI0011', 0, 0, 0, 6, 0 )
    113              
     110
    114111          END SELECT
    115112
  • palm/trunk/SOURCE/vdi_internal_controls.f90

    r4481 r4497  
    11!> @file vdi_internal_controls.f90
    2 !--------------------------------------------------------------------------------!
     2!--------------------------------------------------------------------------------------------------!
    33! This file is part of the PALM model system.
    44!
    5 ! PALM is free software: you can redistribute it and/or modify it under the
    6 ! terms of the GNU General Public License as published by the Free Software
    7 ! Foundation, either version 3 of the License, or (at your option) any later
    8 ! version.
    9 !
    10 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
    11 ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
    12 ! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
    13 !
    14 ! You should have received a copy of the GNU General Public License along with
    15 ! PALM. If not, see <http://www.gnu.org/licenses/>.
    16 !
    17 ! Copyright 2019-2020 Leibniz Universitaet Hannover
    18 !--------------------------------------------------------------------------------!
     5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
     6! Public License as published by the Free Software Foundation, either version 3 of the License, or
     7! (at your option) any later version.
     8!
     9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
     10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
     11! Public License for more details.
     12!
     13! You should have received a copy of the GNU General Public License along with PALM. If not, see
     14! <http://www.gnu.org/licenses/>.
     15!
     16! Copyright 1997-2020 Leibniz Universitaet Hannover
     17!--------------------------------------------------------------------------------------------------!
     18!
    1919!
    2020! Current revisions:
     
    2525! -----------------
    2626! $Id$
     27! file re-formatted to follow the PALM coding standard
     28!
     29!
     30! 4481 2020-03-31 18:55:54Z maronga
    2731! missing preprocessor directive added
    28 ! 
     32!
    2933! 4346 2019-12-18 11:55:56Z motisi
    30 ! Introduction of wall_flags_total_0, which currently sets bits based on static
    31 ! topography information used in wall_flags_static_0
    32 ! 
     34! Introduction of wall_flags_total_0, which currently sets bits based on static topography
     35! information used in wall_flags_static_0
     36!
    3337! 4329 2019-12-10 15:46:36Z motisi
    3438! Renamed wall_flags_0 to wall_flags_static_0
    35 ! 
     39!
    3640! 4182 2019-08-22 15:20:23Z scharf
    3741! added "Authors" section
    38 ! 
     42!
    3943! 4175 2019-08-20 13:19:16Z gronemeier
    4044! bugfix: removed unused variables
     
    5054! Description:
    5155! ------------
    52 !> According to VDI Guideline 3783 Part 9, internal assessment have to be
    53 !> carried out within the program for the model to be considered as evaluated.
    54 !------------------------------------------------------------------------------!
     56!> According to VDI Guideline 3783 Part 9, internal assessment has to be carried out within the
     57!> program for the model to be considered as evaluated.
     58!--------------------------------------------------------------------------------------------------!
    5559 MODULE vdi_internal_controls
    5660
    57     USE arrays_3d,                          &
    58         ONLY:  dzw,                         &
    59                pt,                          &
    60                q,                           &
    61                u,                           &
    62                u_p,                         &
    63                v,                           &
     61    USE arrays_3d,                                                                                 &
     62        ONLY:  dzw,                                                                                &
     63               pt,                                                                                 &
     64               q,                                                                                  &
     65               u,                                                                                  &
     66               u_p,                                                                                &
     67               v,                                                                                  &
    6468               w
    6569
    66     USE control_parameters,                 &
    67         ONLY:  bc_dirichlet_l,              &
    68                bc_dirichlet_n,              &
    69                bc_dirichlet_r,              &
    70                bc_dirichlet_s,              &
    71                bc_lr_cyc,                   &
    72                bc_ns_cyc,                   &
    73                humidity,                    &
    74                end_time,                    &
    75                message_string,              &
    76                neutral,                     &
     70    USE control_parameters,                                                                        &
     71        ONLY:  bc_dirichlet_l,                                                                     &
     72               bc_dirichlet_n,                                                                     &
     73               bc_dirichlet_r,                                                                     &
     74               bc_dirichlet_s,                                                                     &
     75               bc_lr_cyc,                                                                          &
     76               bc_ns_cyc,                                                                          &
     77               end_time,                                                                           &
     78               humidity,                                                                           &
     79               message_string,                                                                     &
     80               neutral,                                                                            &
    7781               time_since_reference_point
    7882
    79     USE indices,                            &
    80         ONLY:  nx,                          &
    81                nxl,                         &
    82                nxlg,                        &
    83                nxr,                         &
    84                nxrg,                        &
    85                ny,                          &
    86                nyn,                         &
    87                nyng,                        &
    88                nys,                         &
    89                nysg,                        &
    90                nzb,                         &
    91                nzt,                         &
     83    USE indices,                                                                                   &
     84        ONLY:  nx,                                                                                 &
     85               nxl,                                                                                &
     86               nxlg,                                                                               &
     87               nxr,                                                                                &
     88               nxrg,                                                                               &
     89               ny,                                                                                 &
     90               nyn,                                                                                &
     91               nyng,                                                                               &
     92               nys,                                                                                &
     93               nysg,                                                                               &
     94               nzb,                                                                                &
     95               nzt,                                                                                &
    9296               wall_flags_total_0
    9397
     
    9599
    96100#if defined( __parallel )
    97     USE pegrid,                             &
    98         ONLY:  collective_wait,             &
    99                comm2d,                      &
    100                ierr,                        &
    101                MPI_DOUBLE_PRECISION,        &
    102                MPI_INTEGER,                 &
    103                MPI_MAX,                     &
    104                MPI_SUM,                     &
     101    USE pegrid,                                                                                    &
     102        ONLY:  collective_wait,                                                                    &
     103               comm2d,                                                                             &
     104               ierr,                                                                               &
     105               MPI_DOUBLE_PRECISION,                                                               &
     106               MPI_INTEGER,                                                                        &
     107               MPI_MAX,                                                                            &
     108               MPI_SUM,                                                                            &
    105109               myid
    106110#else
    107     USE pegrid,                             &
     111    USE pegrid,                                                                                    &
    108112        ONLY:  myid
    109113#endif
    110114
    111115
    112     USE grid_variables,                     &
    113         ONLY:  dx,                          &
     116    USE grid_variables,                                                                            &
     117        ONLY:  dx,                                                                                 &
    114118               dy
    115119
    116     USE pmc_interface,                      &
     120    USE pmc_interface,                                                                             &
    117121        ONLY: nested_run
    118122
    119123    IMPLICIT NONE
    120    
     124
    121125    INTEGER(iwp) ::  internal_count = 0  !< counts calls to this module
    122126
     
    151155!
    152156!-- Public functions
    153     PUBLIC          &
     157    PUBLIC                                                                                         &
    154158       vdi_actions
    155159
     
    157161 CONTAINS
    158162
    159 !------------------------------------------------------------------------------!
     163!--------------------------------------------------------------------------------------------------!
    160164! Description:
    161165! ------------
    162166!> Call for all grid points
    163167!> @todo Add proper description
    164 !------------------------------------------------------------------------------!
     168!--------------------------------------------------------------------------------------------------!
    165169 SUBROUTINE vdi_actions( location )
    166170
    167     CHARACTER (LEN=*), INTENT(IN) ::  location  !< call location string
     171    CHARACTER(LEN=*), INTENT(IN) ::  location  !< call location string
    168172
    169173
     
    173177
    174178          internal_count = internal_count + 1
    175        
     179
    176180          CALL vdi_2_deltat_wave
    177181
     
    191195
    192196 END SUBROUTINE vdi_actions
    193 !------------------------------------------------------------------------------!
     197!--------------------------------------------------------------------------------------------------!
    194198! Description:
    195199! ------------
    196 !> At a control grid point in the interior of the model domain,
    197 !> 2 * deltat waves are not to be generated with increasing simulation time.
    198 !------------------------------------------------------------------------------!
     200!> At a control grid point in the interior of the model domain, 2 * delta t waves are not to be
     201!> generated with increasing simulation time.
     202!--------------------------------------------------------------------------------------------------!
    199203 SUBROUTINE vdi_2_deltat_wave
    200204
    201205    INTEGER(iwp) ::  count_wave = 0  !< counts the number of consecutive waves
    202     INTEGER(iwp) ::  count_time = 0  !< counter, so that the waves follow each other without gaps
    203     INTEGER(iwp) ::  cgp_i = 0       !< x coordinate of the control grid point for testing 2deltat waves
    204     INTEGER(iwp) ::  cgp_j = 0       !< y coordinate of the control grid point for testing 2deltat waves
    205     INTEGER(iwp) ::  cgp_k = 0       !< z coordinate of the control grid point for testing 2deltat waves
    206 
    207     INTEGER(iwp), DIMENSION(4) ::  sig_arr = (/ 0, 0, 0, 0/)  !< indicates an increase(1) or a decrease (0) of u in the last four time steps
     206    INTEGER(iwp) ::  count_time = 0  !< counter, so that the waves follow one another without gaps
     207    INTEGER(iwp) ::  cgp_i      = 0  !< x coordinate of the control grid point for testing 2 delta t waves
     208    INTEGER(iwp) ::  cgp_j      = 0  !< y coordinate of the control grid point for testing 2 delta t waves
     209    INTEGER(iwp) ::  cgp_k      = 0  !< z coordinate of the control grid point for testing 2 delta t waves
     210
     211    INTEGER(iwp), DIMENSION(4) ::  sig_arr = (/ 0, 0, 0, 0 /)  !< indicates an increase(1) or a decrease (0)
     212                                                               !< of u in the last four time steps
    208213
    209214    REAL(wp) ::  random  !< random number
     
    222227!
    223228!--       If there is topography in the entire grid column, a new x coordinate is chosen
    224           IF ( cgp_k >= nzt-1 )  THEN
     229          IF ( cgp_k >= nzt -1 )  THEN
    225230             CALL RANDOM_NUMBER( random )
    226231             cgp_i = nxl + FLOOR( ( nxr + 1 - nxl ) * random )
     
    230235    ENDIF
    231236
    232     CALL testing_2_deltat_wave( u_p(cgp_k,cgp_j,cgp_i), u(cgp_k,cgp_j,cgp_i), &
     237    CALL testing_2_deltat_wave( u_p(cgp_k,cgp_j,cgp_i), u(cgp_k,cgp_j,cgp_i),                      &
    233238                                sig_arr, count_wave, count_time )
    234239
     
    236241
    237242
    238 !------------------------------------------------------------------------------!
     243!--------------------------------------------------------------------------------------------------!
    239244! Description:
    240245! ------------
    241 !> In this subroutine a quantity quant is tested for 2 delta t waves.
    242 !> For this, the size must have a wave-shaped course over 4*4 time steps
    243 !> and the amplitude of the wave has to be greater than the change of quant with
    244 !> increasing time.
    245 !------------------------------------------------------------------------------!
    246 SUBROUTINE testing_2_deltat_wave( quant_p_r, quant_r, sig_arr, count_wave, count_time )
     246!> In this subroutine the quantity quant is tested for 2 delta t waves. For this, the size must have
     247!> a wave-shaped course over 4*4 time steps and the amplitude of the wave has to be greater than the
     248!> change of quant with increasing time.
     249!--------------------------------------------------------------------------------------------------!
     250 SUBROUTINE testing_2_deltat_wave( quant_p_r, quant_r, sig_arr, count_wave, count_time )
    247251
    248252    INTEGER(iwp), INTENT(INOUT) ::  count_wave        !< counts the number of consecutive waves
    249     INTEGER(iwp), INTENT(INOUT) ::  count_time        !< counter, so that the waves follow each other without gaps
     253    INTEGER(iwp), INTENT(INOUT) ::  count_time        !< counter, so that the waves follow one another without gaps
    250254    INTEGER(iwp), PARAMETER     ::  number_wave = 10  !< number of consecutive waves that are not allowed
    251255
    252     REAL(wp), INTENT(IN) ::  quant_p_r                !< quantity from the previous time step as a real
    253     REAL(wp), INTENT(IN) ::  quant_r                  !< quantity as a real
    254     REAL(wp)             ::  quant_rel = 0.0_wp       !< rel. change of the quantity to the previous time step
    255 
    256     INTEGER(iwp), DIMENSION(4), INTENT(INOUT) ::  sig_arr !< indicates an increase(1) or a decrease (0) of
    257                                                           !> quantity quant in the last four time steps
     256    INTEGER(iwp), DIMENSION(4), INTENT(INOUT) ::  sig_arr  !< indicates an increase (1) or a decrease (0) of
     257                                                           !> quantity quant in the last four time steps
     258
     259    REAL(wp), INTENT(IN) ::  quant_p_r           !< quantity from the previous time step as a real
     260    REAL(wp), INTENT(IN) ::  quant_r             !< quantity as a real
     261    REAL(wp)             ::  quant_rel = 0.0_wp  !< rel. change of the quantity to the previous time step
     262
     263
    258264
    259265
     
    267273
    268274!
    269 !-- With this criterion 2 delta t waves are detected if the amplitude of
    270 !-- the wave is greater than the change of quant with increasing time
     275!-- With this criterion 2 delta t waves are detected if the amplitude of the wave is greater than
     276!-- the change of quant with increasing time
    271277    IF ( ALL( sig_arr(1:4) == (/ 1, 0, 1, 0 /) )  .AND.  quant_rel > 0.01 )  THEN
    272278
     
    297303
    298304
    299 !------------------------------------------------------------------------------!
     305!--------------------------------------------------------------------------------------------------!
    300306! Description:
    301307! ------------
    302 !> In this internal assessment the maxima of standarddifferences of the
    303 !> meteorological variables, computed layer by layer will be checked.
    304 !> The maxima should not to remain at the open edges of the model or
    305 !> travel from there into the interior of the domain with increasing
    306 !> simulation time.
     308!> In this internal assessment the maxima of standard differences of the meteorological variables,
     309!> computed layer by layer will be checked. The maxima should not remain at the open edges of the
     310!> model or travel from there into the interior of the domain with increasing simulation time.
    307311!> @todo try to reduce repeating code.
    308 !------------------------------------------------------------------------------!
    309 SUBROUTINE vdi_standard_differences
    310 
    311     INTEGER(iwp) ::  position_u_deviation = 0     !< position of the maximum of the standard deviation of u
    312     INTEGER(iwp) ::  position_u_deviation_p = 0   !< position of the maximum of the standard deviation of u to the previous time step
    313     INTEGER(iwp) ::  position_u_deviation_pp = 0  !< position of the maximum of the standard deviation of u two time steps ago
    314     INTEGER(iwp) ::  position_v_deviation = 0     !< position of the maximum of the standard deviation of v
    315     INTEGER(iwp) ::  position_v_deviation_p = 0   !< position of the maximum of the standard deviation of v to the previous time step
    316     INTEGER(iwp) ::  position_v_deviation_pp = 0  !< position of the maximum of the standard deviation of v two time steps ago
    317     INTEGER(iwp) ::  position_w_deviation = 0     !< position of the maximum of the standard deviation of w
    318     INTEGER(iwp) ::  position_w_deviation_p = 0   !< position of the maximum of the standard deviation of w to the previous time step
    319     INTEGER(iwp) ::  position_w_deviation_pp = 0  !< position of the maximum of the standard deviation of w two time steps ago
    320     INTEGER(iwp) ::  position_pt_deviation = 0    !< position of the maximum of the standard deviation of pt
    321     INTEGER(iwp) ::  position_pt_deviation_p = 0  !< position of the maximum of the standard deviation of pt to the previous time step
    322     INTEGER(iwp) ::  position_pt_deviation_pp = 0 !< position of the maximum of the standard deviation of pt two time steps ago
    323     INTEGER(iwp) ::  position_q_deviation = 0     !< position of the maximum of the standard deviation of q
    324     INTEGER(iwp) ::  position_q_deviation_p = 0   !< position of the maximum of the standard deviation of q to the previous time step
    325     INTEGER(iwp) ::  position_q_deviation_pp = 0  !< position of the maximum of the standard deviation of q two time steps ago
    326 
    327     REAL(wp), DIMENSION(nzb:nzt+1) ::  u_deviation  !< standard deviation of u depending on k
    328     REAL(wp), DIMENSION(nzb:nzt+1) ::  v_deviation  !< standard deviation of v depending on k
    329     REAL(wp), DIMENSION(nzb:nzt+1) ::  w_deviation  !< standard deviation of w depending on k
    330     REAL(wp), DIMENSION(nzb:nzt+1) ::  pt_deviation !< standard deviation of pt depending on k
    331     REAL(wp), DIMENSION(nzb:nzt+1) ::  q_deviation  !< standard deviation of q depending on k
     312!--------------------------------------------------------------------------------------------------!
     313 SUBROUTINE vdi_standard_differences
     314
     315    INTEGER(iwp) ::  position_pt_deviation    = 0  !< position of the maximum of the standard deviation of pt
     316    INTEGER(iwp) ::  position_pt_deviation_p  = 0  !< position of the maximum of the standard deviation of pt
     317                                                   !< to the previous time step
     318    INTEGER(iwp) ::  position_pt_deviation_pp = 0  !< position of the maximum of the standard deviation of pt two time steps ago
     319    INTEGER(iwp) ::  position_q_deviation     = 0  !< position of the maximum of the standard deviation of q
     320    INTEGER(iwp) ::  position_q_deviation_p   = 0  !< position of the maximum of the standard deviation of q to
     321                                                   !< the previous time step
     322    INTEGER(iwp) ::  position_q_deviation_pp  = 0  !< position of the maximum of the standard deviation of q two time steps ago
     323    INTEGER(iwp) ::  position_u_deviation     = 0  !< position of the maximum of the standard deviation of u
     324    INTEGER(iwp) ::  position_u_deviation_p   = 0  !< position of the maximum of the standard deviation of u to
     325                                                   !< the previous time step
     326    INTEGER(iwp) ::  position_u_deviation_pp  = 0  !< position of the maximum of the standard deviation of u two time steps ago
     327    INTEGER(iwp) ::  position_v_deviation     = 0  !< position of the maximum of the standard deviation of v
     328    INTEGER(iwp) ::  position_v_deviation_p   = 0  !< position of the maximum of the standard deviation of v
     329                                                   !< to the previous time step
     330    INTEGER(iwp) ::  position_v_deviation_pp  = 0  !< position of the maximum of the standard deviation of v two time steps ago
     331    INTEGER(iwp) ::  position_w_deviation     = 0  !< position of the maximum of the standard deviation of w
     332    INTEGER(iwp) ::  position_w_deviation_p   = 0  !< position of the maximum of the standard deviation of w
     333                                                   !< to the previous time step
     334    INTEGER(iwp) ::  position_w_deviation_pp  = 0  !< position of the maximum of the standard deviation of w two time steps ago
     335
     336    REAL(wp), DIMENSION(nzb:nzt+1) ::  pt_deviation  !< standard deviation of pt depending on k
     337    REAL(wp), DIMENSION(nzb:nzt+1) ::  q_deviation   !< standard deviation of q depending on k
     338    REAL(wp), DIMENSION(nzb:nzt+1) ::  u_deviation   !< standard deviation of u depending on k
     339    REAL(wp), DIMENSION(nzb:nzt+1) ::  v_deviation   !< standard deviation of v depending on k
     340    REAL(wp), DIMENSION(nzb:nzt+1) ::  w_deviation   !< standard deviation of w depending on k
    332341
    333342!
     
    337346!
    338347!-- Determination of the position of the maximum
    339     position_u_deviation = MAXLOC( u_deviation, DIM=1 )
     348    position_u_deviation = MAXLOC( u_deviation, DIM = 1 )
    340349
    341350!
     
    354363!
    355364!-- Determination of the position of the maximum
    356     position_v_deviation = MAXLOC( v_deviation, DIM=1 )
     365    position_v_deviation = MAXLOC( v_deviation, DIM = 1 )
    357366
    358367!
     
    371380!
    372381!-- Determination of the position of the maximum
    373     position_w_deviation = MAXLOC( w_deviation, DIM=1 )
     382    position_w_deviation = MAXLOC( w_deviation, DIM = 1 )
    374383
    375384!
     
    389398!
    390399!--    Determination of the position of the maximum
    391        position_pt_deviation = MAXLOC( pt_deviation, DIM=1 )
     400       position_pt_deviation = MAXLOC( pt_deviation, DIM = 1 )
    392401
    393402!
    394403!--    Check the position of the maximum of the standard deviation of pt
    395404       IF ( internal_count > 2 )  THEN
    396           CALL check_position( position_pt_deviation,   &
    397                                position_pt_deviation_p, &
     405          CALL check_position( position_pt_deviation,                                              &
     406                               position_pt_deviation_p,                                            &
    398407                               position_pt_deviation_pp )
    399408       ENDIF
     
    411420!
    412421!--    Determination of the position of the maximum
    413        position_q_deviation = MAXLOC( q_deviation, DIM=1 )
     422       position_q_deviation = MAXLOC( q_deviation, DIM = 1 )
    414423
    415424!
    416425!--    Check the position of the maximum of the standard deviation of q
    417426       IF ( internal_count > 2 )  THEN
    418           CALL check_position( position_q_deviation,   &
    419                                position_q_deviation_p, &
     427          CALL check_position( position_q_deviation,                                               &
     428                               position_q_deviation_p,                                             &
    420429                               position_q_deviation_pp )
    421430       ENDIF
     
    426435    ENDIF
    427436
    428 END SUBROUTINE vdi_standard_differences
    429 
    430 
    431 !------------------------------------------------------------------------------!
     437 END SUBROUTINE vdi_standard_differences
     438
     439
     440!--------------------------------------------------------------------------------------------------!
    432441! Description:
    433442! ------------
    434443!> Calculation of the standard deviation
    435 !------------------------------------------------------------------------------!
    436 SUBROUTINE calc_standard_deviation( quant, std_deviation, quant_type )
     444!--------------------------------------------------------------------------------------------------!
     445 SUBROUTINE calc_standard_deviation( quant, std_deviation, quant_type )
    437446
    438447    INTEGER(iwp)             ::  i           !< loop index
     
    468477             flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), quant_type ) )
    469478             quant_av_k_l(k) = quant_av_k_l(k) + quant(k,j,i) * flag
    470              count_2d_l(k)   = count_2d_l(k) + INT( flag, KIND=iwp )
     479             count_2d_l(k)   = count_2d_l(k) + INT( flag, KIND = iwp )
    471480          ENDDO
    472481       ENDDO
     
    474483
    475484#if defined( __parallel )
    476     CALL MPI_ALLREDUCE( quant_av_k_l, quant_av_k, nzt+1-nzb+1, &
    477                         MPI_REAL, MPI_SUM, comm2d, ierr )
    478 
    479     CALL MPI_ALLREDUCE( count_2d_l, count_2d, nzt+1-nzb+1,     &
    480                         MPI_INTEGER, MPI_SUM, comm2d, ierr )
     485    CALL MPI_ALLREDUCE( quant_av_k_l, quant_av_k, nzt+1 - nzb+1, MPI_REAL, MPI_SUM, comm2d, ierr )
     486
     487    CALL MPI_ALLREDUCE( count_2d_l, count_2d, nzt+1 - nzb+1, MPI_INTEGER, MPI_SUM, comm2d, ierr )
    481488#else
    482489    quant_av_k = quant_av_k_l
     
    485492
    486493    DO  k = nzb+1, nzt+1
    487        quant_av_k(k) = quant_av_k(k) / REAL( count_2d(k), KIND=wp )
     494       quant_av_k(k) = quant_av_k(k) / REAL( count_2d(k), KIND = wp )
    488495    ENDDO
    489496
     
    491498       DO  j = nys, nyn
    492499          DO  k = nzb+1, nzt+1
    493              std_deviation_l(k) = std_deviation_l(k)                  &
    494                                 + ( quant(k,j,i) - quant_av_k(k) )**2 &
    495                                 * MERGE( 1.0_wp, 0.0_wp,              &
    496                                          BTEST( wall_flags_total_0(k,j,i), quant_type ) )
     500             std_deviation_l(k) = std_deviation_l(k)                                               &
     501                                + ( quant(k,j,i) - quant_av_k(k) )**2                              &
     502                                * MERGE( 1.0_wp, 0.0_wp,                                           &
     503                                  BTEST( wall_flags_total_0(k,j,i), quant_type ) )
    497504          ENDDO
    498505       ENDDO
     
    501508
    502509#if defined( __parallel )
    503     CALL MPI_ALLREDUCE( std_deviation_l, std_deviation, nzt+1-nzb+1, &
    504                         MPI_REAL, MPI_SUM, comm2d, ierr )
     510    CALL MPI_ALLREDUCE( std_deviation_l, std_deviation, nzt+1-nzb+1, MPI_REAL, MPI_SUM, comm2d, ierr )
    505511#else
    506512    std_deviation = std_deviation_l
     
    508514
    509515    DO  k = nzb+1, nzt+1
    510        std_deviation(k) = SQRT( std_deviation(k) / REAL( count_2d(k), KIND=wp ) )
     516       std_deviation(k) = SQRT( std_deviation(k) / REAL( count_2d(k), KIND = wp ) )
    511517    ENDDO
    512518
    513 END SUBROUTINE calc_standard_deviation
    514 
    515 
    516 !------------------------------------------------------------------------------!
     519 END SUBROUTINE calc_standard_deviation
     520
     521
     522!--------------------------------------------------------------------------------------------------!
    517523! Description:
    518524! ------------
    519 !> Tests for the position of the maxima of the standard deviation.
    520 !> If the maxima remain at the open edges of the model or travel from
    521 !> the open edges into the interior of the domain with increasing
     525!> Tests for the position of the maxima of the standard deviation. If the maxima remain at the open
     526!> edges of the model or travel from the open edges into the interior of the domain with increasing
    522527!> simulation time, the simulation should be aborted.
    523 !------------------------------------------------------------------------------!
     528!--------------------------------------------------------------------------------------------------!
    524529 SUBROUTINE check_position( position_std_deviation, position_std_deviation_p, &
    525530                            position_std_deviation_pp )
     
    530535
    531536
    532     IF ( position_std_deviation == nzt    .AND.  &
    533          position_std_deviation_p == nzt  .AND.  &
    534          position_std_deviation_pp == nzt        )  THEN
    535        message_string = 'The maxima of the standard deviation remain ' // &
    536                         'at the open edges of the model.'
     537    IF ( position_std_deviation == nzt    .AND.                                                    &
     538         position_std_deviation_p == nzt  .AND.                                                    &
     539         position_std_deviation_pp == nzt       )  THEN
     540       message_string = 'The maxima of the standard deviation' //                                  &
     541                        'remain at the open edges of the model.'
    537542       CALL message( 'vdi_standard_differences', 'PA0663', 1, 2, 0, 6, 0 )
    538543    ENDIF
    539544
    540     IF ( position_std_deviation == nzt-2    .AND. &
    541          position_std_deviation_p == nzt-1  .AND. &
     545    IF ( position_std_deviation == nzt-2    .AND.                                                  &
     546         position_std_deviation_p == nzt-1  .AND.                                                  &
    542547         position_std_deviation_pp == nzt         )  THEN
    543        message_string = 'The maxima of the standard deviation travel ' // &
    544                         'from the open edges into the interior ' //       &
     548       message_string = 'The maxima of the standard deviation travel ' //                          &
     549                        'from the open edges into the interior ' //                                &
    545550                        'of the domain with increasing simulation time.'
    546551       CALL message( 'vdi_standard_differences', 'PA0664', 1, 2, 0, 6, 0 )
    547552    ENDIF
    548553
    549 END SUBROUTINE check_position
    550 
    551 
    552 !------------------------------------------------------------------------------!
     554 END SUBROUTINE check_position
     555
     556
     557!--------------------------------------------------------------------------------------------------!
    553558! Description:
    554559! ------------
    555 !> In this control it will be checked, if the means of the meteorological
    556 !> variables over the model grid are not to exhibit 2 deltat waves or
    557 !> monotonic increase or decrease with increasing simulation time.
    558 !------------------------------------------------------------------------------!
    559 SUBROUTINE vdi_domain_averages
    560 
    561    INTEGER(iwp) ::  mono_count_u = 0   !< counter for monotonic decrease or increase of u
    562    INTEGER(iwp) ::  mono_count_v = 0   !< counter for monotonic decrease or increase of v
    563    INTEGER(iwp) ::  mono_count_w = 0   !< counter for monotonic decrease or increase of w
    564    INTEGER(iwp) ::  mono_count_q = 0   !< counter for monotonic decrease or increase of q
    565    INTEGER(iwp) ::  mono_count_pt = 0  !< counter for monotonic decrease or increase of pt
    566    INTEGER(iwp) ::  count_time_u = 0   !< counter, so that the waves of u follow each other without gaps
    567    INTEGER(iwp) ::  count_time_v = 0   !< counter, so that the waves of v follow each other without gaps
    568    INTEGER(iwp) ::  count_time_w = 0   !< counter, so that the waves of w follow each other without gaps
    569    INTEGER(iwp) ::  count_time_q = 0   !< counter, so that the waves of q follow each other without gaps
    570    INTEGER(iwp) ::  count_time_pt = 0  !< counter, so that the waves of pt follow each other without gaps
    571    INTEGER(iwp) ::  count_wave_u = 0   !< counts the number of consecutive waves of u
    572    INTEGER(iwp) ::  count_wave_v = 0   !< counts the number of consecutive waves of v
    573    INTEGER(iwp) ::  count_wave_w = 0   !< counts the number of consecutive waves of w
    574    INTEGER(iwp) ::  count_wave_q = 0   !< counts the number of consecutive waves of q
    575    INTEGER(iwp) ::  count_wave_pt = 0  !< counts the number of consecutive waves of pt
    576 
    577    INTEGER(iwp), DIMENSION(4) ::  sig_u_arr = (/ 0, 0, 0, 0/)   !< indicates an increase(1) or a decrease (0) of u in the last four time steps
    578    INTEGER(iwp), DIMENSION(4) ::  sig_v_arr = (/ 0, 0, 0, 0/)   !< indicates an increase(1) or a decrease (0) of v in the last four time steps
    579    INTEGER(iwp), DIMENSION(4) ::  sig_w_arr = (/ 0, 0, 0, 0/)   !< indicates an increase(1) or a decrease (0) of w in the last four time steps
    580    INTEGER(iwp), DIMENSION(4) ::  sig_q_arr = (/ 0, 0, 0, 0/)   !< indicates an increase(1) or a decrease (0) of q in the last four time steps
    581    INTEGER(iwp), DIMENSION(4) ::  sig_pt_arr = (/ 0, 0, 0, 0/)  !< indicates an increase(1) or a decrease (0) of pt in the last four time steps
    582 
    583    REAL(wp) ::  u_av = 0.0_wp     !< Mean of u
    584    REAL(wp) ::  u_av_p = 0.0_wp   !< Mean of u at the previous time step
    585    REAL(wp) ::  v_av = 0.0_wp     !< Mean of v
    586    REAL(wp) ::  v_av_p = 0.0_wp   !< Mean of v at the previous time step
    587    REAL(wp) ::  w_av = 0.0_wp     !< Mean of w
    588    REAL(wp) ::  w_av_p = 0.0_wp   !< Mean of w at the previous time step
    589    REAL(wp) ::  q_av = 0.0_wp     !< Mean of q
    590    REAL(wp) ::  q_av_p = 0.0_wp   !< Mean of q at the previous time step
    591    REAL(wp) ::  pt_av = 0.0_wp    !< Mean of pt
    592    REAL(wp) ::  pt_av_p = 0.0_wp  !< Mean of pt at the previous time step
     560!> In this control it will be checked, if the means of the meteorological variables over the model
     561!> grid are not to exhibit 2 delta t waves or monotonic increase or decrease with increasing
     562!> simulation time.
     563!--------------------------------------------------------------------------------------------------!
     564 SUBROUTINE vdi_domain_averages
     565
     566    INTEGER(iwp) ::  count_time_u  = 0  !< counter, so that the waves of u follow each other without gaps
     567    INTEGER(iwp) ::  count_time_v  = 0  !< counter, so that the waves of v follow each other without gaps
     568    INTEGER(iwp) ::  count_time_w  = 0  !< counter, so that the waves of w follow each other without gaps
     569    INTEGER(iwp) ::  count_time_q  = 0  !< counter, so that the waves of q follow each other without gaps
     570    INTEGER(iwp) ::  count_time_pt = 0  !< counter, so that the waves of pt follow each other without gaps
     571    INTEGER(iwp) ::  count_wave_u  = 0  !< counts the number of consecutive waves of u
     572    INTEGER(iwp) ::  count_wave_v  = 0  !< counts the number of consecutive waves of v
     573    INTEGER(iwp) ::  count_wave_w  = 0  !< counts the number of consecutive waves of w
     574    INTEGER(iwp) ::  count_wave_q  = 0  !< counts the number of consecutive waves of q
     575    INTEGER(iwp) ::  count_wave_pt = 0  !< counts the number of consecutive waves of pt
     576    INTEGER(iwp) ::  mono_count_u  = 0  !< counter for monotonic decrease or increase of u
     577    INTEGER(iwp) ::  mono_count_v  = 0  !< counter for monotonic decrease or increase of v
     578    INTEGER(iwp) ::  mono_count_w  = 0  !< counter for monotonic decrease or increase of w
     579    INTEGER(iwp) ::  mono_count_q  = 0  !< counter for monotonic decrease or increase of q
     580    INTEGER(iwp) ::  mono_count_pt = 0  !< counter for monotonic decrease or increase of pt
     581
     582    INTEGER(iwp), DIMENSION(4) ::  sig_u_arr = (/ 0, 0, 0, 0/)   !< indicates an increase(1) or a decrease (0)
     583                                                                 !< of u in the last four time steps
     584    INTEGER(iwp), DIMENSION(4) ::  sig_v_arr = (/ 0, 0, 0, 0/)   !< indicates an increase(1) or a decrease (0)
     585                                                                 !< of v in the last four time steps
     586    INTEGER(iwp), DIMENSION(4) ::  sig_w_arr = (/ 0, 0, 0, 0/)   !< indicates an increase(1) or a decrease (0)
     587                                                                 !< of w in the last four time steps
     588    INTEGER(iwp), DIMENSION(4) ::  sig_q_arr = (/ 0, 0, 0, 0/)   !< indicates an increase(1) or a decrease (0)
     589                                                                 !< of q in the last four time steps
     590    INTEGER(iwp), DIMENSION(4) ::  sig_pt_arr = (/ 0, 0, 0, 0/)  !< indicates an increase(1) or a decrease (0)
     591                                                                 !< of pt in the last four time steps
     592
     593    REAL(wp) ::  pt_av = 0.0_wp    !< Mean of pt
     594    REAL(wp) ::  pt_av_p = 0.0_wp  !< Mean of pt at the previous time step
     595    REAL(wp) ::  q_av = 0.0_wp     !< Mean of q
     596    REAL(wp) ::  q_av_p = 0.0_wp   !< Mean of q at the previous time step
     597    REAL(wp) ::  u_av = 0.0_wp     !< Mean of u
     598    REAL(wp) ::  u_av_p = 0.0_wp   !< Mean of u at the previous time step
     599    REAL(wp) ::  v_av = 0.0_wp     !< Mean of v
     600    REAL(wp) ::  v_av_p = 0.0_wp   !< Mean of v at the previous time step
     601    REAL(wp) ::  w_av = 0.0_wp     !< Mean of w
     602    REAL(wp) ::  w_av_p = 0.0_wp   !< Mean of w at the previous time step
    593603
    594604!
     
    626636    ENDIF
    627637
    628     IF ( time_since_reference_point >= end_time  .AND.   &
     638    IF ( time_since_reference_point >= end_time  .AND.                                             &
    629639         mono_count_u > 0.9_wp * internal_count )  THEN
    630640
    631        message_string = 'Monotonic decrease or increase with ' // &
    632                         'increasing simulation time for u'
     641       message_string = 'Monotonic decrease or increase with increasing simulation time for u'
    633642       CALL message( 'vdi_domain_averages', 'PA0665', 0, 1, 0, 6, 0 )
    634643    ENDIF
     
    640649    ENDIF
    641650
    642     IF ( time_since_reference_point >= end_time  .AND.   &
     651    IF ( time_since_reference_point >= end_time  .AND.                                             &
    643652         mono_count_v > 0.9_wp * internal_count )  THEN
    644        message_string = 'Monotonic decrease or increase with ' // &
    645                         'increasing simulation time for v'
     653       message_string = 'Monotonic decrease or increase with increasing simulation time for v'
    646654       CALL message( 'vdi_domain_averages', 'PA0665', 0, 1, 0, 6, 0 )
    647655    ENDIF
     
    653661    ENDIF
    654662
    655     IF ( time_since_reference_point >= end_time  .AND.   &
     663    IF ( time_since_reference_point >= end_time  .AND.                                             &
    656664         mono_count_w > 0.9_wp * internal_count )  THEN
    657        message_string = 'Monotonic decrease or increase with ' // &
    658                         'increasing simulation time for w'
     665       message_string = 'Monotonic decrease or increase with increasing simulation time for w'
    659666       CALL message( 'vdi_domain_averages', 'PA0665', 0, 1, 0, 6, 0 )
    660667    ENDIF
     
    667674       ENDIF
    668675
    669        IF ( time_since_reference_point >= end_time  .AND.    &
     676       IF ( time_since_reference_point >= end_time  .AND.                                          &
    670677            mono_count_pt > 0.9_wp * internal_count )  THEN
    671           message_string = 'Monotonic decrease or increase with ' // &
    672                            'increasing simulation time for pt'
     678          message_string = 'Monotonic decrease or increase with increasing simulation time for pt'
    673679          CALL message( 'vdi_domain_averages', 'PA0665', 0, 1, 0, 6, 0 )
    674680       ENDIF
     
    682688       ENDIF
    683689
    684        IF ( time_since_reference_point >= end_time  .AND.   &
     690       IF ( time_since_reference_point >= end_time  .AND.                                          &
    685691            mono_count_q > 0.9_wp * internal_count )  THEN
    686           message_string = 'Monotonic decrease or increase with ' // &
    687                            'increasing simulation time for q'
     692          message_string = 'Monotonic decrease or increase with increasing simulation time for q'
    688693          CALL message( 'vdi_domain_averages', 'PA0665', 0, 1, 0, 6, 0 )
    689694       ENDIF
     
    707712
    708713
    709 !------------------------------------------------------------------------------!
     714!--------------------------------------------------------------------------------------------------!
    710715! Description:
    711716! ------------
    712717!> Calculate the average of a quantity 'quant'.
    713 !------------------------------------------------------------------------------!
     718!--------------------------------------------------------------------------------------------------!
    714719 SUBROUTINE calc_average( quant, quant_av, quant_type )
    715720
    716     INTEGER(iwp) ::  average_count = 0    !< counter for averaging
     721    INTEGER(iwp) ::  average_count   = 0  !< counter for averaging
    717722    INTEGER(iwp) ::  average_count_l = 0  !< counter for averaging (local)
    718723    INTEGER      ::  i                    !< loop index
     
    721726    INTEGER(iwp) ::  quant_type           !< bit position (1 for u, 2 for v, 3 for w and 0 for scalar)
    722727
    723     REAL(wp) ::  flag                     !< flag indicating atmosphere (1) or wall (0) grid point
    724     REAL(wp) ::  quant_av                 !< average of the quantity quant
    725     REAL(wp) ::  quant_av_l = 0.0_wp      !< average of the quantity quant (local)
    726 
    727     REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  quant
     728    REAL(wp) ::  flag                 !< flag indicating atmosphere (1) or wall (0) grid point
     729    REAL(wp) ::  quant_av             !< average of the quantity quant
     730    REAL(wp) ::  quant_av_l = 0.0_wp  !< average of the quantity quant (local)
     731
     732    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  quant  !<
    728733
    729734!
     
    736741            flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), quant_type ) )
    737742            quant_av_l = quant_av_l + quant(k,j,i) * flag
    738             average_count_l = average_count_l + INT( flag, KIND=iwp )
     743            average_count_l = average_count_l + INT( flag, KIND = iwp )
    739744         ENDDO
    740745      ENDDO
     
    742747
    743748#if defined( __parallel )
    744     CALL MPI_ALLREDUCE( quant_av_l, quant_av, 1,        &
    745                         MPI_REAL, MPI_SUM, comm2d, ierr )
    746     CALL MPI_ALLREDUCE( average_count_l, average_count, 1, &
    747                         MPI_INTEGER, MPI_SUM, comm2d, ierr )
     749    CALL MPI_ALLREDUCE( quant_av_l, quant_av, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
     750    CALL MPI_ALLREDUCE( average_count_l, average_count, 1, MPI_INTEGER, MPI_SUM, comm2d, ierr )
    748751#else
    749752    quant_av = quant_av_l
     
    756759
    757760
    758 !------------------------------------------------------------------------------!
     761!--------------------------------------------------------------------------------------------------!
    759762! Description:
    760763! ------------
    761764!> Testing for conservation of mass.
    762 !------------------------------------------------------------------------------!
     765!--------------------------------------------------------------------------------------------------!
    763766 SUBROUTINE vdi_conservation_of_mass
    764767
    765     INTEGER(iwp) ::  i              !< loop index
    766     INTEGER(iwp) ::  j              !< loop index
    767     INTEGER(iwp) ::  k              !< loop index
     768    INTEGER(iwp) ::  i  !< loop index
     769    INTEGER(iwp) ::  j  !< loop index
     770    INTEGER(iwp) ::  k  !< loop index
    768771
    769772    REAL(wp)     ::  sum_mass_flux  !< sum of the mass flow
    770773
     774    REAL(wp), DIMENSION(1:3) ::  volume_flow    !< volume flow
    771775    REAL(wp), DIMENSION(1:3) ::  volume_flow_l  !< volume flow (local)
    772     REAL(wp), DIMENSION(1:3) ::  volume_flow    !< volume flow
    773776
    774777
     
    783786       DO  j = nys, nyn
    784787          DO  k = nzb+1, nzt
    785              volume_flow_l(1) = volume_flow_l(1)                        &
    786                               + u(k,j,i) * dzw(k) * dy                  &
    787                               * MERGE( 1.0_wp, 0.0_wp,                  &
    788                                  BTEST( wall_flags_total_0(k,j,i), 1 )  &
    789                                      )
     788             volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k) * dy                          &
     789                                * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
    790790          ENDDO
    791791       ENDDO
    792792    ENDIF
    793 ! 
     793!
    794794!-- Sum up the volume flow through the right boundary
    795795    IF ( nxr == nx )  THEN
     
    797797       DO  j = nys, nyn
    798798          DO  k = nzb+1, nzt
    799              volume_flow_l(1) = volume_flow_l(1)                        &
    800                               - u(k,j,i) * dzw(k) * dy                  &
    801                               * MERGE( 1.0_wp, 0.0_wp,                  &
    802                                  BTEST( wall_flags_total_0(k,j,i), 1 )  &
    803                                      )
     799             volume_flow_l(1) = volume_flow_l(1) - u(k,j,i) * dzw(k) * dy                          &
     800                                * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
    804801          ENDDO
    805802       ENDDO
     
    812809       DO  i = nxl, nxr
    813810          DO  k = nzb+1, nzt
    814              volume_flow_l(2) = volume_flow_l(2)                        &
    815                               + v(k,j,i) * dzw(k) * dx                  &
    816                               * MERGE( 1.0_wp, 0.0_wp,                  &
    817                                  BTEST( wall_flags_total_0(k,j,i), 2 )  &
    818                                      )
     811             volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k) * dx                          &
     812                                * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
    819813          ENDDO
    820814       ENDDO
     
    826820       DO  i = nxl, nxr
    827821          DO  k = nzb+1, nzt
    828              volume_flow_l(2) = volume_flow_l(2)                        &
    829                               - v(k,j,i) * dzw(k) * dx                  &
    830                               * MERGE( 1.0_wp, 0.0_wp,                  &
    831                                  BTEST( wall_flags_total_0(k,j,i), 2 )  &
    832                                      )
     822             volume_flow_l(2) = volume_flow_l(2) - v(k,j,i) * dzw(k) * dx                          &
     823                                * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
    833824          ENDDO
    834825       ENDDO
     
    857848   ENDIF
    858849
    859 END SUBROUTINE vdi_conservation_of_mass
    860 
    861 
    862 !------------------------------------------------------------------------------!
     850 END SUBROUTINE vdi_conservation_of_mass
     851
     852
     853!--------------------------------------------------------------------------------------------------!
    863854! Description:
    864855! ------------
    865 !> The results will be checked for exceedance the specified limits.
    866 !> The controls are performed at every time step and at every grid point.
    867 !> No wind component is allowed to have a magnitude greater than ten times
    868 !> the maximum wind velocity at the approach flow profile (Vdi 3783 part 9).
    869 !> Note, that the supersaturation can not be higher than 10%. Therefore, no
    870 !> test is required.
    871 !------------------------------------------------------------------------------!
    872 SUBROUTINE vdi_plausible_values
    873 
    874     INTEGER(iwp) ::  i          !< loop index
    875     INTEGER(iwp) ::  j          !< loop index
    876     INTEGER(iwp) ::  k          !< loop index
    877 
     856!> The results will be checked for exceedance of the specified limits. The controls are performed at
     857!> every time step and at every grid point. No wind component is allowed to have a magnitude greater
     858!> than ten times the maximum wind velocity at the approach flow profile (Vdi 3783 part 9).
     859!> Note, that the supersaturation can not be higher than 10%. Therefore, no test is required.
     860!--------------------------------------------------------------------------------------------------!
     861 SUBROUTINE vdi_plausible_values
     862
     863    INTEGER(iwp) ::  i  !< loop index
     864    INTEGER(iwp) ::  j  !< loop index
     865    INTEGER(iwp) ::  k  !< loop index
     866
     867    REAL(wp)     :: max_uv      !< maximum speed of all edges
    878868    REAL(wp)     :: max_uv_l_l  !< maximum speed at the left edge (local)
    879869    REAL(wp)     :: max_uv_l    !< maximum speed at the left edge
     870    REAL(wp)     :: max_uv_n_l  !< maximum speed at the north edge (local)
     871    REAL(wp)     :: max_uv_n    !< maximum speed at the north edge
    880872    REAL(wp)     :: max_uv_r_l  !< maximum speed at the right edge (local)
    881873    REAL(wp)     :: max_uv_r    !< maximum speed at the right edge
    882874    REAL(wp)     :: max_uv_s_l  !< maximum speed at the south edge (local)
    883875    REAL(wp)     :: max_uv_s    !< maximum speed at the south edge
    884     REAL(wp)     :: max_uv_n_l  !< maximum speed at the north edge (local)
    885     REAL(wp)     :: max_uv_n    !< maximum speed at the north edge
    886     REAL(wp)     :: max_uv      !< maximum speed of all edges
    887 
    888     REAL(wp), DIMENSION(4)                 ::  max_arr    !<
    889     REAL(wp), DIMENSION(:), ALLOCATABLE    ::  uv         !< wind velocity at the approach flow
    890     REAL(wp), DIMENSION(:), ALLOCATABLE    ::  uv_l       !< wind velocity at the approach flow (local)
     876
     877    REAL(wp), DIMENSION(4)                 ::  max_arr  !<
     878    REAL(wp), DIMENSION(:), ALLOCATABLE    ::  uv       !< wind velocity at the approach flow
     879    REAL(wp), DIMENSION(:), ALLOCATABLE    ::  uv_l     !< wind velocity at the approach flow (local)
    891880
    892881    REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn) ::  uv_l_nest  !< wind profile at the left edge (nesting)
     
    916905       IF ( nxl == 0 )  THEN
    917906          i = nxl
    918           DO j = nys, nyn
    919              DO k = nzb, nzt+1
    920                 uv_l_nest(k,j) = SQRT( ( 0.5_wp * ( u(k,j,i-1) + u(k,j,i) ) )**2  &
    921                                      + ( 0.5_wp * ( v(k,j-1,i) + v(k,j,i) ) )**2  )
     907          DO  j = nys, nyn
     908             DO  k = nzb, nzt+1
     909                uv_l_nest(k,j) = SQRT( ( 0.5_wp * ( u(k,j,i-1) + u(k,j,i) ) )**2                   &
     910                                     + ( 0.5_wp * ( v(k,j-1,i) + v(k,j,i) ) )**2 )
    922911             ENDDO
    923912          ENDDO
     
    926915!
    927916!--    Right boundary
    928        IF( nxr == nx )  THEN
     917       IF ( nxr == nx )  THEN
    929918          i = nxr
    930           DO j = nys, nyn
    931              DO k = nzb, nzt+1
    932                 uv_r_nest(k,j) = SQRT( ( 0.5_wp * ( u(k,j,i-1) + u(k,j,i) ) )**2  &
    933                                      + ( 0.5_wp * ( v(k,j-1,i) + v(k,j,i) ) )**2  )
     919          DO  j = nys, nyn
     920             DO  k = nzb, nzt+1
     921                uv_r_nest(k,j) = SQRT( ( 0.5_wp * ( u(k,j,i-1) + u(k,j,i) ) )**2                   &
     922                                     + ( 0.5_wp * ( v(k,j-1,i) + v(k,j,i) ) )**2 )
    934923
    935924             ENDDO
     
    941930       IF ( nys == 0 )  THEN
    942931          j = nys
    943           DO i = nxl, nxr
    944              DO k = nzb, nzt+1
    945                 uv_s_nest(k,i) = SQRT( ( 0.5_wp * ( u(k,j,i-1) + u(k,j,i) ) )**2  &
    946                                      + ( 0.5_wp * ( v(k,j-1,i) + v(k,j,i) ) )**2  )
     932          DO  i = nxl, nxr
     933             DO  k = nzb, nzt+1
     934                uv_s_nest(k,i) = SQRT( ( 0.5_wp * ( u(k,j,i-1) + u(k,j,i) ) )**2                   &
     935                                     + ( 0.5_wp * ( v(k,j-1,i) + v(k,j,i) ) )**2 )
    947936             ENDDO
    948937          ENDDO
     
    953942       IF ( nyn == ny )  THEN
    954943          j = nyn
    955           DO i = nxl, nxr
    956              DO k = nzb, nzt+1
    957                 uv_n_nest(k,i) = SQRT( ( 0.5_wp * ( u(k,j,i-1) + u(k,j,i) ) )**2  &
    958                                      + ( 0.5_wp * ( v(k,j-1,i) + v(k,j,i) ) )**2  )
     944          DO  i = nxl, nxr
     945             DO  k = nzb, nzt+1
     946                uv_n_nest(k,i) = SQRT( ( 0.5_wp * ( u(k,j,i-1) + u(k,j,i) ) )**2                   &
     947                                     + ( 0.5_wp * ( v(k,j-1,i) + v(k,j,i) ) )**2 )
    959948
    960949             ENDDO
     
    983972          IF ( nxl == 0  .AND.  nys == 0 )  THEN
    984973             DO  k = nzb, nzt+1
    985                 uv_l(k) = SQRT( ( 0.5_wp * ( u(k,0,-1) + u(k,0,0) ) )**2  &
    986                               + ( 0.5_wp * ( v(k,-1,0) + v(k,0,0) ) )**2  )
     974                uv_l(k) = SQRT( ( 0.5_wp * ( u(k,0,-1) + u(k,0,0) ) )**2                           &
     975                              + ( 0.5_wp * ( v(k,-1,0) + v(k,0,0) ) )**2 )
    987976             ENDDO
    988977          ENDIF
     
    993982          IF ( nxl == 0 .AND. nys == 0 )  THEN
    994983             DO  k = nzb, nzt+1
    995                 uv_l(k) = SQRT( ( 0.5_wp * ( u(k,0,-1) + u(k,0,0) ) )**2  &
    996                               + ( 0.5_wp * ( v(k,-1,0) + v(k,0,0) ) )**2  )
     984                uv_l(k) = SQRT( ( 0.5_wp * ( u(k,0,-1) + u(k,0,0) ) )**2                           &
     985                              + ( 0.5_wp * ( v(k,-1,0) + v(k,0,0) ) )**2 )
    997986             ENDDO
    998987          ENDIF
     
    1001990          IF ( nxr == nx .AND. nys == 0 )  THEN
    1002991             DO  k = nzb, nzt+1
    1003                 uv_l(k) = SQRT( ( 0.5_wp * ( u(k,0,nxr) + u(k,0,nxr+1) ) )**2  &
    1004                               + ( 0.5_wp * ( v(k,-1,nxr) + v(k,0,nxr) ) )**2   )
     992                uv_l(k) = SQRT( ( 0.5_wp * ( u(k,0,nxr) + u(k,0,nxr+1) ) )**2                      &
     993                              + ( 0.5_wp * ( v(k,-1,nxr) + v(k,0,nxr) ) )**2 )
    1005994             ENDDO
    1006995          ENDIF
     
    1010999          IF ( nxl == 0 .AND. nyn == ny )  THEN
    10111000             DO  k = nzb, nzt+1
    1012                 uv_l(k) = SQRT( ( 0.5_wp * ( u(k,nyn,-1) + u(k,nyn,0) ) )**2  &
     1001                uv_l(k) = SQRT( ( 0.5_wp * ( u(k,nyn,-1) + u(k,nyn,0) ) )**2                       &
    10131002                              + ( 0.5_wp * ( v(k,nyn+1,0) + v(k,nyn,0) ) )**2 )
    10141003             ENDDO
     
    10181007          IF ( nxl == 0 .AND. nys == 0 )  THEN
    10191008             DO  k = nzb, nzt+1
    1020                 uv_l(k) = SQRT( ( 0.5_wp * ( u(k,0,-1) + u(k,0,0) ) )**2  &
    1021                               + ( 0.5_wp * ( v(k,-1,0) + v(k,0,0) ) )**2  )
     1009                uv_l(k) = SQRT( ( 0.5_wp * ( u(k,0,-1) + u(k,0,0) ) )**2                           &
     1010                              + ( 0.5_wp * ( v(k,-1,0) + v(k,0,0) ) )**2 )
    10221011             ENDDO
    10231012          ENDIF
     
    10351024
    10361025!
    1037 !-- Test for exceedance the specified limits
    1038     message_string = 'A wind component have a magnitude greater ' //  &
    1039                      'than ten times the maximum wind velocity ' //   &
    1040                      'at the approach flow profile.'
     1026!-- Test for exceedance of the specified limits
     1027    message_string = 'A wind component have a magnitude greater than ten times the maximum' //     &
     1028                     'wind velocity at the approach flow profile.'
    10411029
    10421030    IF ( MAXVAL( ABS( u ) ) > 10.0_wp * max_uv )  THEN
     
    10551043!-- Test if the potential temperature lies between 220 K and 330 K
    10561044    IF ( MAXVAL( pt ) > 330.0_wp .OR. MAXVAL( pt ) < 220.0_wp )  THEN
    1057        message_string = 'The potential temperature does not lie ' //  &
    1058                         'between 220 K and 330 K.'
     1045       message_string = 'The potential temperature does not lie between 220 K and 330 K.'
    10591046       CALL message( 'vdi_plausible_values', 'PA0668', 2, 2, myid, 6, 0 )
    10601047    ENDIF
  • palm/trunk/SOURCE/virtual_flight_mod.f90

    r4495 r4497  
    11!> @file virtual_flights_mod.f90
    2 !------------------------------------------------------------------------------!
     2!--------------------------------------------------------------------------------------------------!
    33! This file is part of the PALM model system.
    44!
    5 ! PALM is free software: you can redistribute it and/or modify it under the
    6 ! terms of the GNU General Public License as published by the Free Software
    7 ! Foundation, either version 3 of the License, or (at your option) any later
    8 ! version.
    9 !
    10 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
    11 ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
    12 ! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
    13 !
    14 ! You should have received a copy of the GNU General Public License along with
    15 ! PALM. If not, see <http://www.gnu.org/licenses/>.
     5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
     6! Public License as published by the Free Software Foundation, either version 3 of the License, or
     7! (at your option) any later version.
     8!
     9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
     10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
     11! Public License for more details.
     12!
     13! You should have received a copy of the GNU General Public License along with PALM. If not, see
     14! <http://www.gnu.org/licenses/>.
    1615!
    1716! Copyright 1997-2020 Leibniz Universitaet Hannover
    18 !------------------------------------------------------------------------------!
     17!--------------------------------------------------------------------------------------------------!
     18!
    1919!
    2020! Current revisions:
    21 ! ------------------
     21! -----------------
    2222!
    2323!
     
    2525! -----------------
    2626! $Id$
     27! file re-formatted to follow the PALM coding standard
     28!
     29!
     30! 4495 2020-04-13 20:11:20Z raasch
    2731! restart data handling with MPI-IO added
    2832!
    2933! 4360 2020-01-07 11:25:50Z suehring
    3034! Corrected "Former revisions" section
    31 ! 
     35!
    3236! 4004 2019-05-24 11:32:38Z suehring
    3337! Allow variable start- and end locations also in return mode
    34 ! 
     38!
    3539! 3885 2019-04-11 11:29:34Z kanani
    36 ! Changes related to global restructuring of location messages and introduction
    37 ! of additional debug messages
    38 ! 
     40! Changes related to global restructuring of location messages and introduction of additional
     41! debug messages
     42!
    3943! 3655 2019-01-07 16:51:22Z knoop
    4044! variables documented
    41 ! 
     45!
    4246! 1957 2016-07-07 10:43:48Z suehring
    4347! Initial revision
     
    4751!> Module for virtual flight measurements.
    4852!> @todo Err msg PA0438: flight can be inside topography -> extra check?
    49 !--------------------------------------------------------------------------------!
     53!--------------------------------------------------------------------------------------------------!
    5054 MODULE flight_mod
    51  
     55
    5256    USE control_parameters,                                                                        &
    53         ONLY:  debug_output, fl_max, num_leg, num_var_fl, num_var_fl_user,                         &
    54                restart_data_format_output, virtual_flight
    55  
     57        ONLY:  debug_output,                                                                       &
     58               fl_max, num_leg,                                                                    &
     59               num_var_fl,                                                                         &
     60               num_var_fl_user,                                                                    &
     61               restart_data_format_output,                                                         &
     62               virtual_flight
     63
    5664    USE kinds
    5765
     
    6068
    6169
    62     CHARACTER(LEN=6), DIMENSION(fl_max) ::  leg_mode = 'cyclic'  !< flight mode through the model domain, either 'cyclic' or 'return'
     70    CHARACTER(LEN=6), DIMENSION(fl_max) ::  leg_mode = 'cyclic'  !< flight mode through the model domain, either 'cyclic' or
     71                                                                 !<'return'
    6372
    6473    INTEGER(iwp) ::  l           !< index for flight leg
    6574    INTEGER(iwp) ::  var_index   !< index for measured variable
    6675
    67     LOGICAL, DIMENSION(:), ALLOCATABLE  ::  cyclic_leg !< flag to identify fly mode
     76    LOGICAL, DIMENSION(:), ALLOCATABLE  ::  cyclic_leg  !< flag to identify fly mode
    6877
    6978    REAL(wp) ::  flight_end = 9999999.9_wp  !< end time of virtual flight
    7079    REAL(wp) ::  flight_begin = 0.0_wp      !< end time of virtual flight
    7180
    72     REAL(wp), DIMENSION(fl_max) ::  flight_angle = 45.0_wp   !< angle determining the horizontal flight direction
    73     REAL(wp), DIMENSION(fl_max) ::  flight_level = 100.0_wp  !< flight level
    74     REAL(wp), DIMENSION(fl_max) ::  max_elev_change = 0.0_wp !< maximum elevation change for the respective flight leg
    75     REAL(wp), DIMENSION(fl_max) ::  rate_of_climb = 0.0_wp   !< rate of climb or descent
    76     REAL(wp), DIMENSION(fl_max) ::  speed_agl = 25.0_wp      !< absolute horizontal flight speed above ground level (agl)
    77     REAL(wp), DIMENSION(fl_max) ::  x_start = 999999999.0_wp !< start x position
    78     REAL(wp), DIMENSION(fl_max) ::  x_end   = 999999999.0_wp !< end x position
    79     REAL(wp), DIMENSION(fl_max) ::  y_start = 999999999.0_wp !< start y position
    80     REAL(wp), DIMENSION(fl_max) ::  y_end   = 999999999.0_wp !< end y position
    81 
    82     REAL(wp), DIMENSION(:), ALLOCATABLE ::  u_agl      !< u-component of flight speed
    83     REAL(wp), DIMENSION(:), ALLOCATABLE ::  v_agl      !< v-component of flight speed
    84     REAL(wp), DIMENSION(:), ALLOCATABLE ::  w_agl      !< w-component of flight speed
    85     REAL(wp), DIMENSION(:), ALLOCATABLE ::  x_pos      !< aircraft x-position
    86     REAL(wp), DIMENSION(:), ALLOCATABLE ::  y_pos      !< aircraft y-position
    87     REAL(wp), DIMENSION(:), ALLOCATABLE ::  z_pos      !< aircraft z-position
    88 
    89     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sensor_l !< measured data on local PE
    90     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sensor   !< measured data
     81    REAL(wp), DIMENSION(fl_max) ::  flight_angle = 45.0_wp    !< angle determining the horizontal flight direction
     82    REAL(wp), DIMENSION(fl_max) ::  flight_level = 100.0_wp   !< flight level
     83    REAL(wp), DIMENSION(fl_max) ::  max_elev_change = 0.0_wp  !< maximum elevation change for the respective flight leg
     84    REAL(wp), DIMENSION(fl_max) ::  rate_of_climb = 0.0_wp    !< rate of climb or descent
     85    REAL(wp), DIMENSION(fl_max) ::  speed_agl = 25.0_wp       !< absolute horizontal flight speed above ground level (agl)
     86    REAL(wp), DIMENSION(fl_max) ::  x_start = 999999999.0_wp  !< start x position
     87    REAL(wp), DIMENSION(fl_max) ::  x_end   = 999999999.0_wp  !< end x position
     88    REAL(wp), DIMENSION(fl_max) ::  y_start = 999999999.0_wp  !< start y position
     89    REAL(wp), DIMENSION(fl_max) ::  y_end   = 999999999.0_wp  !< end y position
     90
     91    REAL(wp), DIMENSION(:), ALLOCATABLE ::  u_agl  !< u-component of flight speed
     92    REAL(wp), DIMENSION(:), ALLOCATABLE ::  v_agl  !< v-component of flight speed
     93    REAL(wp), DIMENSION(:), ALLOCATABLE ::  w_agl  !< w-component of flight speed
     94    REAL(wp), DIMENSION(:), ALLOCATABLE ::  x_pos  !< aircraft x-position
     95    REAL(wp), DIMENSION(:), ALLOCATABLE ::  y_pos  !< aircraft y-position
     96    REAL(wp), DIMENSION(:), ALLOCATABLE ::  z_pos  !< aircraft z-position
     97
     98    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sensor_l  !< measured data on local PE
     99    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  sensor    !< measured data
    91100
    92101    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  var_u  !< dummy array for possibly user-defined quantities
     
    99108       MODULE PROCEDURE flight_header
    100109    END INTERFACE flight_header
    101    
     110
    102111    INTERFACE flight_init
    103112       MODULE PROCEDURE flight_init
     
    123132       MODULE PROCEDURE flight_measurement
    124133    END INTERFACE flight_measurement
    125    
    126     INTERFACE flight_rrd_global 
     134
     135    INTERFACE flight_rrd_global
    127136       MODULE PROCEDURE flight_rrd_global_ftn
    128137       MODULE PROCEDURE flight_rrd_global_mpi
    129138    END INTERFACE flight_rrd_global
    130    
    131     INTERFACE flight_wrd_global 
    132        MODULE PROCEDURE flight_wrd_global 
    133     END INTERFACE flight_wrd_global 
     139
     140    INTERFACE flight_wrd_global
     141       MODULE PROCEDURE flight_wrd_global
     142    END INTERFACE flight_wrd_global
    134143
    135144!
    136145!-- Private interfaces
    137     PRIVATE flight_check_parameters, flight_init_output, interpolate_xyz
     146    PRIVATE flight_check_parameters,                                                               &
     147            flight_init_output,                                                                    &
     148            interpolate_xyz
    138149!
    139150!-- Public interfaces
    140     PUBLIC flight_init, flight_header, flight_parin, flight_measurement,       &
    141            flight_wrd_global, flight_rrd_global                   
     151    PUBLIC flight_init,                                                                            &
     152           flight_header,                                                                          &
     153           flight_parin,                                                                           &
     154           flight_measurement,                                                                     &
     155           flight_wrd_global,                                                                      &
     156           flight_rrd_global
    142157!
    143158!-- Public variables
    144     PUBLIC fl_max, sensor, x_pos, y_pos, z_pos
     159    PUBLIC fl_max,                                                                                 &
     160           sensor,                                                                                 &
     161           x_pos,                                                                                  &
     162           y_pos,                                                                                  &
     163           z_pos
    145164
    146165 CONTAINS
    147166
    148 !------------------------------------------------------------------------------!
     167!--------------------------------------------------------------------------------------------------!
    149168! Description:
    150169! ------------
    151170!> Header output for flight module.
    152 !------------------------------------------------------------------------------!
    153     SUBROUTINE flight_header ( io )
    154 
    155    
    156        IMPLICIT NONE
    157 
    158        INTEGER(iwp), INTENT(IN) ::  io            !< Unit of the output file
    159 
    160        WRITE ( io, 1  )
    161        WRITE ( io, 2  )
    162        WRITE ( io, 3  ) num_leg
    163        WRITE ( io, 4  ) flight_begin
    164        WRITE ( io, 5  ) flight_end
    165        
    166        DO l=1, num_leg
    167           WRITE ( io, 6   ) l
    168           WRITE ( io, 7   ) speed_agl(l)
    169           WRITE ( io, 8   ) flight_level(l)
    170           WRITE ( io, 9   ) max_elev_change(l)
    171           WRITE ( io, 10  ) rate_of_climb(l)
    172           WRITE ( io, 11  ) leg_mode(l)
    173        ENDDO
    174 
    175        
    176      1   FORMAT (' Virtual flights:'/                                           &
    177                ' ----------------')
    178      2   FORMAT ('       Output every timestep')
    179      3   FORMAT ('       Number of flight legs:',    I3                )
    180      4   FORMAT ('       Begin of measurements:',    F10.1    , ' s'   )
    181      5   FORMAT ('       End of measurements:',      F10.1    , ' s'   )
    182      6   FORMAT ('       Leg', I3/,                                             &
    183                 '       ------' )
    184      7   FORMAT ('          Flight speed            : ', F5.1, ' m/s' )
    185      8   FORMAT ('          Flight level            : ', F5.1, ' m'   )
    186      9   FORMAT ('          Maximum elevation change: ', F5.1, ' m/s' )
    187      10  FORMAT ('          Rate of climb / descent : ', F5.1, ' m/s' )
    188      11  FORMAT ('          Leg mode                : ', A/           )
    189  
    190     END SUBROUTINE flight_header
    191  
    192 !------------------------------------------------------------------------------!
     171!--------------------------------------------------------------------------------------------------!
     172 SUBROUTINE flight_header ( io )
     173
     174
     175    IMPLICIT NONE
     176
     177    INTEGER(iwp), INTENT(IN) ::  io  !< Unit of the output file
     178
     179    WRITE ( io, 1  )
     180    WRITE ( io, 2  )
     181    WRITE ( io, 3  ) num_leg
     182    WRITE ( io, 4  ) flight_begin
     183    WRITE ( io, 5  ) flight_end
     184
     185    DO  l=1, num_leg
     186       WRITE ( io, 6   ) l
     187       WRITE ( io, 7   ) speed_agl(l)
     188       WRITE ( io, 8   ) flight_level(l)
     189       WRITE ( io, 9   ) max_elev_change(l)
     190       WRITE ( io, 10  ) rate_of_climb(l)
     191       WRITE ( io, 11  ) leg_mode(l)
     192    ENDDO
     193
     194
     195 1   FORMAT (' Virtual flights: ----------------' )
     196 2   FORMAT ('       Output every timestep' )
     197 3   FORMAT ('       Number of flight legs:',    I3 )
     198 4   FORMAT ('       Begin of measurements:',    F10.1    , ' s' )
     199 5   FORMAT ('       End of measurements:',      F10.1    , ' s' )
     200 6   FORMAT ('       Leg', I3/, '       ------' )
     201 7   FORMAT ('          Flight speed            : ', F5.1, ' m/s' )
     202 8   FORMAT ('          Flight level            : ', F5.1, ' m' )
     203 9   FORMAT ('          Maximum elevation change: ', F5.1, ' m/s' )
     204 10  FORMAT ('          Rate of climb / descent : ', F5.1, ' m/s' )
     205 11  FORMAT ('          Leg mode                : ', A/ )
     206
     207 END SUBROUTINE flight_header
     208
     209!--------------------------------------------------------------------------------------------------!
    193210! Description:
    194211! ------------
    195212!> Reads the namelist flight_par.
    196 !------------------------------------------------------------------------------!
    197     SUBROUTINE flight_parin
    198 
    199        USE control_parameters,                                                 &
    200            ONLY:  message_string
    201      
    202        IMPLICIT NONE
    203 
    204        CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file
    205 
    206        NAMELIST /flight_par/  flight_angle, flight_end, flight_begin, leg_mode,&
    207                               flight_level, max_elev_change, rate_of_climb,    &
    208                               speed_agl, x_end, x_start, y_end, y_start
    209  
    210 
    211        NAMELIST /virtual_flight_parameters/                                    &
    212                               flight_angle, flight_end, flight_begin, leg_mode,&
    213                               flight_level, max_elev_change, rate_of_climb,    &
    214                               speed_agl, x_end, x_start, y_end, y_start
    215 !
    216 !--    Try to find the namelist flight_par
    217        REWIND ( 11 )
    218        line = ' '
    219        DO WHILE ( INDEX( line, '&virtual_flight_parameters' ) == 0 )
    220           READ ( 11, '(A)', END=12 )  line
    221        ENDDO
    222        BACKSPACE ( 11 )
    223 
    224 !
    225 !--    Read namelist
    226        READ ( 11, virtual_flight_parameters, ERR = 10 )
    227 !
    228 !--    Set switch that virtual flights shall be carried out
    229        virtual_flight = .TRUE.
    230 
    231        GOTO 14
     213!--------------------------------------------------------------------------------------------------!
     214 SUBROUTINE flight_parin
     215
     216    USE control_parameters,                                                                        &
     217        ONLY:  message_string
     218
     219    IMPLICIT NONE
     220
     221    CHARACTER(LEN=80) ::  line  !< dummy string that contains the current line of the parameter file
     222
     223    NAMELIST /flight_par/ flight_angle,                                                            &
     224                          flight_begin,                                                            &
     225                          flight_end,                                                              &
     226                          flight_level,                                                            &
     227                          leg_mode,                                                                &
     228                          max_elev_change,                                                         &
     229                          rate_of_climb,                                                           &
     230                          speed_agl,                                                               &
     231                          x_end,                                                                   &
     232                          x_start,                                                                 &
     233                          y_end,                                                                   &
     234                          y_start
     235
     236
     237    NAMELIST /virtual_flight_parameters/ flight_angle,                                             &
     238                                         flight_begin,                                             &
     239                                         flight_end,                                               &
     240                                         flight_level,                                             &
     241                                         leg_mode,                                                 &
     242                                         max_elev_change,                                          &
     243                                         rate_of_climb,                                            &
     244                                         speed_agl,                                                &
     245                                         x_end,                                                    &
     246                                         x_start,                                                  &
     247                                         y_end,                                                    &
     248                                         y_start
     249!
     250!-- Try to find the namelist flight_par
     251    REWIND ( 11 )
     252    line = ' '
     253    DO  WHILE ( INDEX( line, '&virtual_flight_parameters' ) == 0 )
     254       READ ( 11, '(A)', END = 12 )  line
     255    ENDDO
     256    BACKSPACE ( 11 )
     257
     258!
     259!-- Read namelist
     260    READ ( 11, virtual_flight_parameters, ERR = 10 )
     261!
     262!-- Set switch so that virtual flights shall be carried out
     263    virtual_flight = .TRUE.
     264
     265    GOTO 14
    232266
    233267 10    BACKSPACE( 11 )
     
    238272 12    REWIND ( 11 )
    239273       line = ' '
    240        DO WHILE ( INDEX( line, '&flight_par' ) == 0 )
    241           READ ( 11, '(A)', END=14 )  line
     274       DO  WHILE ( INDEX( line, '&flight_par' ) == 0 )
     275          READ ( 11, '(A)', END = 14 )  line
    242276       ENDDO
    243277       BACKSPACE ( 11 )
    244278
    245279!
    246 !--    Read namelist
    247        READ ( 11, flight_par, ERR = 13, END = 14 )
    248        
    249        message_string = 'namelist flight_par is deprecated and will be ' // &
    250                      'removed in near future.& Please use namelist ' //     &
    251                      'virtual_flight_parameters instead' 
    252        CALL message( 'flight_parin', 'PA0487', 0, 1, 0, 6, 0 )     
    253 !
    254 !--    Set switch that virtual flights shall be carried out
    255        virtual_flight = .TRUE.
    256 
    257        GOTO 14
     280!-- Read namelist
     281    READ ( 11, flight_par, ERR = 13, END = 14 )
     282
     283    message_string = 'namelist flight_par is deprecated and will be ' //                          &
     284                     'removed in near future.& Please use namelist ' //                            &
     285                     'virtual_flight_parameters instead'
     286    CALL message( 'flight_parin', 'PA0487', 0, 1, 0, 6, 0 )
     287!
     288!-- Set switch so that virtual flights shall be carried out
     289    virtual_flight = .TRUE.
     290
     291    GOTO 14
    258292
    259293 13    BACKSPACE( 11 )
     
    263297 14    CONTINUE
    264298
    265     END SUBROUTINE flight_parin
    266 
    267 !------------------------------------------------------------------------------!
     299 END SUBROUTINE flight_parin
     300
     301!--------------------------------------------------------------------------------------------------!
    268302! Description:
    269303! ------------
    270 !> Inititalization of required arrays, number of legs and flags. Moreover,
    271 !> initialize flight speed and -direction, as well as start positions.
    272 !------------------------------------------------------------------------------!
    273     SUBROUTINE flight_init
    274  
    275        USE basic_constants_and_equations_mod,                                  &
    276            ONLY:  pi
    277    
    278        USE control_parameters,                                                 &
    279            ONLY:  initializing_actions
    280                  
    281        USE indices,                                                            &
    282            ONLY:  nxlg, nxrg, nysg, nyng, nzb, nzt
    283 
    284        IMPLICIT NONE
    285 
    286        REAL(wp) ::  distance  !< distance between start and end position of a flight leg
    287 
    288 
    289        IF ( debug_output )  CALL debug_message( 'flight_init', 'start' )
    290 !
    291 !--    Determine the number of flight legs
    292        l = 1
    293        DO WHILE ( x_start(l) /= 999999999.0_wp  .AND.  l <= SIZE(x_start) )
    294           l       = l + 1
     304!> Inititalization of required arrays, number of legs and flags. Moreover, initialize flight speed
     305!> and -direction, as well as start positions.
     306!--------------------------------------------------------------------------------------------------!
     307 SUBROUTINE flight_init
     308
     309    USE basic_constants_and_equations_mod,                                                         &
     310        ONLY:  pi
     311
     312    USE control_parameters,                                                                        &
     313        ONLY:  initializing_actions
     314
     315    USE indices,                                                                                   &
     316        ONLY:  nxlg,                                                                               &
     317               nxrg,                                                                               &
     318               nysg,                                                                               &
     319               nyng,                                                                               &
     320               nzb,                                                                                &
     321               nzt
     322
     323    IMPLICIT NONE
     324
     325    REAL(wp) ::  distance  !< distance between start and end position of a flight leg
     326
     327
     328    IF ( debug_output )  CALL debug_message( 'flight_init', 'start' )
     329!
     330!-- Determine the number of flight legs
     331    l = 1
     332    DO  WHILE ( x_start(l) /= 999999999.0_wp  .AND.  l <= SIZE(x_start) )
     333       l       = l + 1
     334    ENDDO
     335    num_leg = l-1
     336!
     337!-- Check for proper parameter settings
     338    CALL flight_check_parameters
     339!
     340!-- Allocate and initialize logical array for flight pattern
     341    ALLOCATE( cyclic_leg(1:num_leg) )
     342!
     343!-- Initialize flags for cyclic/return legs
     344    DO  l = 1, num_leg
     345       cyclic_leg(l) = MERGE( .TRUE., .FALSE., TRIM( leg_mode(l) ) == 'cyclic' )
     346    ENDDO
     347!
     348!-- Allocate and initialize arraxs for flight position and speed. In case of restart runs these data
     349!-- are read by the routine read_flight_restart_data instead.
     350    IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
     351
     352       ALLOCATE( x_pos(1:num_leg), y_pos(1:num_leg ), z_pos(1:num_leg) )
     353!
     354!--    Inititalize x-, y-, and z-positions with initial start position
     355       x_pos(1:num_leg) = x_start(1:num_leg)
     356       y_pos(1:num_leg) = y_start(1:num_leg)
     357       z_pos(1:num_leg) = flight_level(1:num_leg)
     358!
     359!--    Allocate arrays for flight-speed components
     360       ALLOCATE( u_agl(1:num_leg),                                                                 &
     361                 v_agl(1:num_leg),                                                                 &
     362                 w_agl(1:num_leg) )
     363!
     364!--    Inititalize u-, v- and w-component.
     365       DO  l = 1, num_leg
     366!
     367!--       In case of return-legs, the flight direction, i.e. the horizontal flight-speed components,
     368!--       are derived from the given start/end positions.
     369          IF (  .NOT.  cyclic_leg(l) )  THEN
     370             distance = SQRT( ( x_end(l) - x_start(l) )**2 + ( y_end(l) - y_start(l) )**2 )
     371             u_agl(l) = speed_agl(l) * ( x_end(l) - x_start(l) ) / distance
     372             v_agl(l) = speed_agl(l) * ( y_end(l) - y_start(l) ) / distance
     373             w_agl(l) = rate_of_climb(l)
     374!
     375!--       In case of cyclic-legs, flight direction is directly derived from the given flight angle.
     376          ELSE
     377             u_agl(l) = speed_agl(l) * COS( flight_angle(l) * pi / 180.0_wp )
     378             v_agl(l) = speed_agl(l) * SIN( flight_angle(l) * pi / 180.0_wp )
     379             w_agl(l) = rate_of_climb(l)
     380          ENDIF
     381
    295382       ENDDO
    296        num_leg = l-1
    297 !
    298 !--    Check for proper parameter settings
    299        CALL flight_check_parameters
    300 !
    301 !--    Allocate and initialize logical array for flight pattern
    302        ALLOCATE( cyclic_leg(1:num_leg) )
    303 !
    304 !--    Initialize flags for cyclic/return legs
    305        DO l = 1, num_leg
    306           cyclic_leg(l) = MERGE( .TRUE., .FALSE.,                              &
    307                                  TRIM( leg_mode(l) ) == 'cyclic'               &
    308                                )
    309        ENDDO
    310 !
    311 !--    Allocate and initialize arraxs for flight position and speed. In case of
    312 !--    restart runs these data are read by the routine read_flight_restart_data
    313 !--    instead.
    314        IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
    315        
    316           ALLOCATE( x_pos(1:num_leg), y_pos(1:num_leg ), z_pos(1:num_leg) )
    317 !
    318 !--       Inititalize x-, y-, and z-positions with initial start position         
    319           x_pos(1:num_leg) = x_start(1:num_leg)
    320           y_pos(1:num_leg) = y_start(1:num_leg)
    321           z_pos(1:num_leg) = flight_level(1:num_leg)
    322 !
    323 !--       Allocate arrays for flight-speed components
    324           ALLOCATE( u_agl(1:num_leg),                                          &
    325                     v_agl(1:num_leg),                                          &
    326                     w_agl(1:num_leg) )
    327 !
    328 !--       Inititalize u-, v- and w-component.
    329           DO  l = 1, num_leg
    330 !
    331 !--          In case of return-legs, the flight direction, i.e. the horizontal
    332 !--          flight-speed components, are derived from the given start/end
    333 !--          positions.
    334              IF (  .NOT.  cyclic_leg(l) )  THEN
    335                 distance = SQRT( ( x_end(l) - x_start(l) )**2                  &
    336                                + ( y_end(l) - y_start(l) )**2 )
    337                 u_agl(l) = speed_agl(l) * ( x_end(l) - x_start(l) ) / distance
    338                 v_agl(l) = speed_agl(l) * ( y_end(l) - y_start(l) ) / distance
    339                 w_agl(l) = rate_of_climb(l)
    340 !
    341 !--          In case of cyclic-legs, flight direction is directly derived from
    342 !--          the given flight angle.
    343              ELSE
    344                 u_agl(l) = speed_agl(l) * COS( flight_angle(l) * pi / 180.0_wp )
    345                 v_agl(l) = speed_agl(l) * SIN( flight_angle(l) * pi / 180.0_wp )
    346                 w_agl(l) = rate_of_climb(l)
    347              ENDIF
    348 
    349           ENDDO
    350              
    351        ENDIF   
    352 !
    353 !--    Initialized data output
    354        CALL flight_init_output       
    355 !
    356 !--    Allocate array required for user-defined quantities if necessary.
    357        IF ( num_var_fl_user  > 0 )                                             &
    358           ALLOCATE( var_u(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    359 !
    360 !--    Allocate and initialize arrays containing the measured data
    361        ALLOCATE( sensor_l(1:num_var_fl,1:num_leg) )
    362        ALLOCATE( sensor(1:num_var_fl,1:num_leg)   )
    363        sensor_l = 0.0_wp
    364        sensor   = 0.0_wp
    365 
    366        IF ( debug_output )  CALL debug_message( 'flight_init', 'end' )
    367 
    368     END SUBROUTINE flight_init
    369    
    370    
    371 !------------------------------------------------------------------------------!
     383
     384    ENDIF
     385!
     386!-- Initialized data output
     387    CALL flight_init_output
     388!
     389!-- Allocate array required for user-defined quantities if necessary.
     390    IF ( num_var_fl_user  > 0 )  ALLOCATE( var_u(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     391!
     392!-- Allocate and initialize arrays containing the measured data
     393    ALLOCATE( sensor_l(1:num_var_fl,1:num_leg) )
     394    ALLOCATE( sensor(1:num_var_fl,1:num_leg)   )
     395    sensor_l = 0.0_wp
     396    sensor   = 0.0_wp
     397
     398    IF ( debug_output )  CALL debug_message( 'flight_init', 'end' )
     399
     400 END SUBROUTINE flight_init
     401
     402
     403!--------------------------------------------------------------------------------------------------!
    372404! Description:
    373405! ------------
    374406!> Initialization of output-variable names and units.
    375 !------------------------------------------------------------------------------!
    376     SUBROUTINE flight_init_output
    377    
    378        USE control_parameters,                                                 &
    379           ONLY:  cloud_droplets, humidity, neutral, passive_scalar
    380 
    381        USE bulk_cloud_model_mod,                                               &
    382            ONLY:  bulk_cloud_model
    383 
    384        USE netcdf_interface
    385    
    386        IMPLICIT NONE
    387 
    388        CHARACTER(LEN=10) ::  label_leg  !< dummy argument to convert integer to string
    389        
    390        INTEGER(iwp) ::  i               !< loop variable
    391        INTEGER(iwp) ::  id_pt           !< identifyer for labeling
    392        INTEGER(iwp) ::  id_q            !< identifyer for labeling
    393        INTEGER(iwp) ::  id_ql           !< identifyer for labeling
    394        INTEGER(iwp) ::  id_s            !< identifyer for labeling       
    395        INTEGER(iwp) ::  id_u = 1        !< identifyer for labeling 
    396        INTEGER(iwp) ::  id_v = 2        !< identifyer for labeling
    397        INTEGER(iwp) ::  id_w = 3        !< identifyer for labeling
    398        INTEGER(iwp) ::  k               !< index variable
    399        
    400        LOGICAL      ::  init = .TRUE.   !< flag to distiquish calls of user_init_flight
    401 !
    402 !--    Define output quanities, at least three variables are measured (u,v,w)
    403        num_var_fl = 3
    404        IF ( .NOT. neutral                     )  THEN
    405           num_var_fl = num_var_fl + 1
    406           id_pt      = num_var_fl
    407        ENDIF
    408        IF ( humidity                          )  THEN
    409           num_var_fl = num_var_fl + 1
    410           id_q       = num_var_fl
    411        ENDIF
    412        IF ( bulk_cloud_model .OR. cloud_droplets )  THEN
    413           num_var_fl = num_var_fl + 1 
    414           id_ql      = num_var_fl
    415        ENDIF
    416        IF ( passive_scalar                    )  THEN
    417           num_var_fl = num_var_fl + 1
    418           id_s       = num_var_fl
    419        ENDIF
    420 !
    421 !--    Write output strings for dimensions x, y, z
    422        DO l=1, num_leg
    423 
    424           IF ( l < 10                    )  WRITE( label_leg, '(I1)')  l
    425           IF ( l >= 10   .AND.  l < 100  )  WRITE( label_leg, '(I2)')  l
    426           IF ( l >= 100  .AND.  l < 1000 )  WRITE( label_leg, '(I3)')  l
    427 
    428           dofl_dim_label_x(l)  = 'x_' // TRIM( label_leg )
    429           dofl_dim_label_y(l)  = 'y_' // TRIM( label_leg )
    430           dofl_dim_label_z(l)  = 'z_' // TRIM( label_leg )
    431 
     407!--------------------------------------------------------------------------------------------------!
     408 SUBROUTINE flight_init_output
     409
     410    USE control_parameters,                                                                        &
     411       ONLY:  cloud_droplets,                                                                      &
     412              humidity,                                                                            &
     413              neutral,                                                                             &
     414              passive_scalar
     415
     416    USE bulk_cloud_model_mod,                                                                      &
     417        ONLY:  bulk_cloud_model
     418
     419    USE netcdf_interface
     420
     421    IMPLICIT NONE
     422
     423    CHARACTER(LEN=10) ::  label_leg  !< dummy argument to convert integer to string
     424
     425    INTEGER(iwp) ::  i         !< loop variable
     426    INTEGER(iwp) ::  id_pt     !< identifyer for labeling
     427    INTEGER(iwp) ::  id_q      !< identifyer for labeling
     428    INTEGER(iwp) ::  id_ql     !< identifyer for labeling
     429    INTEGER(iwp) ::  id_s      !< identifyer for labeling
     430    INTEGER(iwp) ::  id_u = 1  !< identifyer for labeling
     431    INTEGER(iwp) ::  id_v = 2  !< identifyer for labeling
     432    INTEGER(iwp) ::  id_w = 3  !< identifyer for labeling
     433    INTEGER(iwp) ::  k         !< index variable
     434
     435    LOGICAL      ::  init = .TRUE.  !< flag to distiquish calls of user_init_flight
     436!
     437!-- Define output quanities, at least three variables are measured (u,v,w)
     438    num_var_fl = 3
     439    IF ( .NOT. neutral                     )  THEN
     440       num_var_fl = num_var_fl + 1
     441       id_pt      = num_var_fl
     442    ENDIF
     443    IF ( humidity                          )  THEN
     444       num_var_fl = num_var_fl + 1
     445       id_q       = num_var_fl
     446    ENDIF
     447    IF ( bulk_cloud_model .OR. cloud_droplets )  THEN
     448       num_var_fl = num_var_fl + 1
     449       id_ql      = num_var_fl
     450    ENDIF
     451    IF ( passive_scalar                    )  THEN
     452       num_var_fl = num_var_fl + 1
     453       id_s       = num_var_fl
     454    ENDIF
     455!
     456!-- Write output strings for dimensions x, y, z
     457    DO  l=1, num_leg
     458
     459       IF ( l < 10                    )  WRITE( label_leg, '(I1)' )  l
     460       IF ( l >= 10   .AND.  l < 100  )  WRITE( label_leg, '(I2)' )  l
     461       IF ( l >= 100  .AND.  l < 1000 )  WRITE( label_leg, '(I3)' )  l
     462
     463       dofl_dim_label_x(l)  = 'x_' // TRIM( label_leg )
     464       dofl_dim_label_y(l)  = 'y_' // TRIM( label_leg )
     465       dofl_dim_label_z(l)  = 'z_' // TRIM( label_leg )
     466
     467    ENDDO
     468
     469!
     470!-- Call user routine to initialize further variables
     471    CALL user_init_flight( init )
     472!
     473!-- Write output labels and units for the quanities
     474    k = 1
     475    DO  l=1, num_leg
     476
     477       IF ( l < 10                    )  WRITE( label_leg, '(I1)' )  l
     478       IF ( l >= 10   .AND.  l < 100  )  WRITE( label_leg, '(I2)' )  l
     479       IF ( l >= 100  .AND.  l < 1000 )  WRITE( label_leg, '(I3)' )  l
     480
     481       label_leg = 'leg_' // TRIM(label_leg)
     482       DO  i=1, num_var_fl
     483
     484          IF ( i == id_u      )  THEN
     485             dofl_label(k) = TRIM( label_leg ) // '_u'
     486             dofl_unit(k)  = 'm/s'
     487             k             = k + 1
     488          ELSEIF ( i == id_v  )  THEN
     489             dofl_label(k) = TRIM( label_leg ) // '_v'
     490             dofl_unit(k)  = 'm/s'
     491             k             = k + 1
     492          ELSEIF ( i == id_w  )  THEN
     493             dofl_label(k) = TRIM( label_leg ) // '_w'
     494             dofl_unit(k)  = 'm/s'
     495             k             = k + 1
     496          ELSEIF ( i == id_pt )  THEN
     497             dofl_label(k) = TRIM( label_leg ) // '_theta'
     498             dofl_unit(k)  = 'K'
     499             k             = k + 1
     500          ELSEIF ( i == id_q  )  THEN
     501             dofl_label(k) = TRIM( label_leg ) // '_q'
     502             dofl_unit(k)  = 'kg/kg'
     503             k             = k + 1
     504          ELSEIF ( i == id_ql )  THEN
     505             dofl_label(k) = TRIM( label_leg ) // '_ql'
     506             dofl_unit(k)  = 'kg/kg'
     507             k             = k + 1
     508          ELSEIF ( i == id_s  )  THEN
     509             dofl_label(k) = TRIM( label_leg ) // '_s'
     510             dofl_unit(k)  = 'kg/kg'
     511             k             = k + 1
     512          ENDIF
    432513       ENDDO
    433        
    434 !
    435 !--    Call user routine to initialize further variables
    436        CALL user_init_flight( init )
    437 !
    438 !--    Write output labels and units for the quanities
    439        k = 1
    440        DO l=1, num_leg
    441        
    442           IF ( l < 10                    )  WRITE( label_leg, '(I1)')  l
    443           IF ( l >= 10   .AND.  l < 100  )  WRITE( label_leg, '(I2)')  l
    444           IF ( l >= 100  .AND.  l < 1000 )  WRITE( label_leg, '(I3)')  l
    445          
    446           label_leg = 'leg_' // TRIM(label_leg)
    447           DO i=1, num_var_fl
    448 
    449              IF ( i == id_u      )  THEN         
    450                 dofl_label(k) = TRIM( label_leg ) // '_u'
    451                 dofl_unit(k)  = 'm/s'
    452                 k             = k + 1
    453              ELSEIF ( i == id_v  )  THEN       
    454                 dofl_label(k) = TRIM( label_leg ) // '_v'
    455                 dofl_unit(k)  = 'm/s'
    456                 k             = k + 1
    457              ELSEIF ( i == id_w  )  THEN         
    458                 dofl_label(k) = TRIM( label_leg ) // '_w'
    459                 dofl_unit(k)  = 'm/s'
    460                 k             = k + 1
    461              ELSEIF ( i == id_pt )  THEN       
    462                 dofl_label(k) = TRIM( label_leg ) // '_theta'
    463                 dofl_unit(k)  = 'K'
    464                 k             = k + 1
    465              ELSEIF ( i == id_q  )  THEN       
    466                 dofl_label(k) = TRIM( label_leg ) // '_q'
    467                 dofl_unit(k)  = 'kg/kg'
    468                 k             = k + 1
    469              ELSEIF ( i == id_ql )  THEN       
    470                 dofl_label(k) = TRIM( label_leg ) // '_ql'
    471                 dofl_unit(k)  = 'kg/kg'
    472                 k             = k + 1
    473              ELSEIF ( i == id_s  )  THEN                         
    474                 dofl_label(k) = TRIM( label_leg ) // '_s'
    475                 dofl_unit(k)  = 'kg/kg'
    476                 k             = k + 1
    477              ENDIF
    478           ENDDO
    479          
    480           DO i=1, num_var_fl_user
    481              CALL user_init_flight( init, k, i, label_leg )
    482           ENDDO
    483          
    484        ENDDO
    485 !
    486 !--    Finally, set the total number of flight-output quantities.
    487        num_var_fl = num_var_fl + num_var_fl_user
    488        
    489     END SUBROUTINE flight_init_output
    490 
    491 !------------------------------------------------------------------------------!
     514
     515       DO i=1, num_var_fl_user
     516          CALL user_init_flight( init, k, i, label_leg )
     517       ENDDO
     518
     519    ENDDO
     520!
     521!-- Finally, set the total number of flight-output quantities.
     522    num_var_fl = num_var_fl + num_var_fl_user
     523
     524 END SUBROUTINE flight_init_output
     525
     526!--------------------------------------------------------------------------------------------------!
    492527! Description:
    493528! ------------
    494 !> This routine calculates the current flight positions and calls the
    495 !> respective interpolation routine to measures the data at the current
    496 !> flight position.
    497 !------------------------------------------------------------------------------!
    498     SUBROUTINE flight_measurement
    499 
    500        USE arrays_3d,                                                          &
    501            ONLY:  ddzu, ddzw, pt, q, ql, s, u, v, w, zu, zw
    502 
    503        USE control_parameters,                                                 &
    504            ONLY:  cloud_droplets, dt_3d, humidity, neutral,                    &
    505                   passive_scalar, simulated_time
    506                  
    507        USE cpulog,                                                             &
    508            ONLY:  cpu_log, log_point
    509 
    510        USE grid_variables,                                                     &
    511            ONLY:  ddx, ddy, dx, dy
    512 
    513        USE indices,                                                            &
    514            ONLY:  nx, nxl, nxr, ny, nys, nyn
    515 
    516        USE bulk_cloud_model_mod,                                               &
    517            ONLY:  bulk_cloud_model
    518 
    519        USE pegrid
    520 
    521        IMPLICIT NONE
    522 
    523        LOGICAL  ::  on_pe  !< flag to check if current flight position is on current PE
    524 
    525        REAL(wp) ::  x  !< distance between left edge of current grid box and flight position
    526        REAL(wp) ::  y  !< distance between south edge of current grid box and flight position
    527 
    528        INTEGER(iwp) ::  i   !< index of current grid box along x
    529        INTEGER(iwp) ::  j   !< index of current grid box along y
    530        INTEGER(iwp) ::  n   !< loop variable for number of user-defined output quantities
    531        
    532        CALL cpu_log( log_point(65), 'flight_measurement', 'start' )
    533 !
    534 !--    Perform flight measurement if start time is reached.
    535        IF ( simulated_time >= flight_begin  .AND.                              &
    536             simulated_time <= flight_end )  THEN
    537 
    538           sensor_l = 0.0_wp
    539           sensor   = 0.0_wp
    540 !
    541 !--       Loop over all flight legs
    542           DO l=1, num_leg
    543 !
    544 !--          Update location for each leg
    545              x_pos(l) = x_pos(l) + u_agl(l) * dt_3d
    546              y_pos(l) = y_pos(l) + v_agl(l) * dt_3d
    547              z_pos(l) = z_pos(l) + w_agl(l) * dt_3d
    548 !
    549 !--          Check if location must be modified for return legs. 
    550 !--          Carry out horizontal reflection if required.
    551              IF ( .NOT. cyclic_leg(l) )  THEN
    552 
    553                 IF ( x_start(l) <= x_end(l) )  THEN
    554 !
    555 !--                Outward flight, i.e. from start to end
    556                    IF ( u_agl(l) >= 0.0_wp  .AND.  x_pos(l) > x_end(l)      )  THEN
    557                       x_pos(l) = 2.0_wp * x_end(l)   - x_pos(l)
    558                       u_agl(l) = - u_agl(l)
    559 !                 
    560 !--                Return flight, i.e. from end to start
    561                    ELSEIF ( u_agl(l) < 0.0_wp  .AND.  x_pos(l) < x_start(l) )  THEN
    562                       x_pos(l) = 2.0_wp * x_start(l) - x_pos(l)
    563                       u_agl(l) = - u_agl(l)
    564                    ENDIF
    565                 ELSE
    566 !
    567 !--                Outward flight, i.e. from start to end
    568                    IF ( u_agl(l) < 0.0_wp  .AND.  x_pos(l) < x_end(l)      )  THEN
    569                       x_pos(l) = 2.0_wp * x_end(l)   - x_pos(l)
    570                       u_agl(l) = - u_agl(l)
    571 !                 
    572 !--                Return flight, i.e. from end to start
    573                    ELSEIF ( u_agl(l) >= 0.0_wp  .AND.  x_pos(l) > x_start(l) )  THEN
    574                       x_pos(l) = 2.0_wp * x_start(l) - x_pos(l)
    575                       u_agl(l) = - u_agl(l)
    576                    ENDIF
     529!> This routine calculates the current flight positions and calls the respective interpolation
     530!> routine to measure the data at the current flight position.
     531!--------------------------------------------------------------------------------------------------!
     532 SUBROUTINE flight_measurement
     533
     534    USE arrays_3d,                                                                                 &
     535        ONLY:  ddzu,                                                                               &
     536               ddzw,                                                                               &
     537               pt,                                                                                 &
     538               q,                                                                                  &
     539               ql,                                                                                 &
     540               s,                                                                                  &
     541               u,                                                                                  &
     542               v,                                                                                  &
     543               w,                                                                                  &
     544               zu,                                                                                 &
     545               zw
     546
     547    USE control_parameters,                                                                        &
     548        ONLY:  cloud_droplets,                                                                     &
     549               dt_3d,                                                                              &
     550               humidity,                                                                           &
     551               neutral,                                                                            &
     552               passive_scalar,                                                                     &
     553               simulated_time
     554
     555    USE cpulog,                                                                                    &
     556        ONLY:  cpu_log,                                                                            &
     557               log_point
     558
     559    USE grid_variables,                                                                            &
     560        ONLY:  ddx,                                                                                &
     561               ddy,                                                                                &
     562               dx,                                                                                 &
     563               dy
     564
     565    USE indices,                                                                                   &
     566        ONLY:  nx,                                                                                 &
     567               nxl,                                                                                &
     568               nxr,                                                                                &
     569               ny,                                                                                 &
     570               nys,                                                                                &
     571               nyn
     572
     573    USE bulk_cloud_model_mod,                                                                      &
     574        ONLY:  bulk_cloud_model
     575
     576    USE pegrid
     577
     578    IMPLICIT NONE
     579
     580    INTEGER(iwp) ::  i  !< index of current grid box along x
     581    INTEGER(iwp) ::  j  !< index of current grid box along y
     582    INTEGER(iwp) ::  n  !< loop variable for number of user-defined output quantities
     583
     584    LOGICAL  ::  on_pe  !< flag to check if current flight position is on current PE
     585
     586    REAL(wp) ::  x  !< distance between left edge of current grid box and flight position
     587    REAL(wp) ::  y  !< distance between south edge of current grid box and flight position
     588
     589    CALL cpu_log( log_point(65), 'flight_measurement', 'start' )
     590!
     591!-- Perform flight measurement if start time is reached.
     592    IF ( simulated_time >= flight_begin  .AND.  simulated_time <= flight_end )  THEN
     593
     594       sensor_l = 0.0_wp
     595       sensor   = 0.0_wp
     596!
     597!--    Loop over all flight legs
     598       DO  l = 1, num_leg
     599!
     600!--       Update location for each leg
     601          x_pos(l) = x_pos(l) + u_agl(l) * dt_3d
     602          y_pos(l) = y_pos(l) + v_agl(l) * dt_3d
     603          z_pos(l) = z_pos(l) + w_agl(l) * dt_3d
     604!
     605!--       Check if location must be modified for return legs.
     606!--       Carry out horizontal reflection if required.
     607          IF ( .NOT. cyclic_leg(l) )  THEN
     608
     609             IF ( x_start(l) <= x_end(l) )  THEN
     610!
     611!--             Outward flight, i.e. from start to end
     612                IF ( u_agl(l) >= 0.0_wp  .AND.  x_pos(l) > x_end(l)      )  THEN
     613                   x_pos(l) = 2.0_wp * x_end(l)   - x_pos(l)
     614                   u_agl(l) = - u_agl(l)
     615!
     616!--             Return flight, i.e. from end to start
     617                ELSEIF ( u_agl(l) < 0.0_wp  .AND.  x_pos(l) < x_start(l) )  THEN
     618                   x_pos(l) = 2.0_wp * x_start(l) - x_pos(l)
     619                   u_agl(l) = - u_agl(l)
    577620                ENDIF
    578                
    579                 IF ( y_start(l) <= y_end(l) )  THEN
    580 !
    581 !--                Outward flight, i.e. from start to end
    582                    IF ( v_agl(l) >= 0.0_wp  .AND.  y_pos(l) > y_end(l)      )  THEN
    583                       y_pos(l) = 2.0_wp * y_end(l)   - y_pos(l)
    584                       v_agl(l) = - v_agl(l)
    585 !                 
    586 !--                Return flight, i.e. from end to start                 
    587                    ELSEIF ( v_agl(l) < 0.0_wp  .AND.  y_pos(l) < y_start(l) )  THEN
    588                       y_pos(l) = 2.0_wp * y_start(l) - y_pos(l)
    589                       v_agl(l) = - v_agl(l)
    590                    ENDIF
    591                 ELSE
    592 !
    593 !--                Outward flight, i.e. from start to end
    594                    IF ( v_agl(l) < 0.0_wp  .AND.  y_pos(l) < y_end(l)      )  THEN
    595                       y_pos(l) = 2.0_wp * y_end(l)   - y_pos(l)
    596                       v_agl(l) = - v_agl(l)
    597 !                 
    598 !--                Return flight, i.e. from end to start                 
    599                    ELSEIF ( v_agl(l) >= 0.0_wp  .AND.  y_pos(l) > y_start(l) )  THEN
    600                       y_pos(l) = 2.0_wp * y_start(l) - y_pos(l)
    601                       v_agl(l) = - v_agl(l)
    602                    ENDIF
     621             ELSE
     622!
     623!--             Outward flight, i.e. from start to end
     624                IF ( u_agl(l) < 0.0_wp  .AND.  x_pos(l) < x_end(l)      )  THEN
     625                   x_pos(l) = 2.0_wp * x_end(l)   - x_pos(l)
     626                   u_agl(l) = - u_agl(l)
     627!
     628!--             Return flight, i.e. from end to start
     629                ELSEIF ( u_agl(l) >= 0.0_wp  .AND.  x_pos(l) > x_start(l) )  THEN
     630                   x_pos(l) = 2.0_wp * x_start(l) - x_pos(l)
     631                   u_agl(l) = - u_agl(l)
    603632                ENDIF
    604 !
    605 !--          Check if flight position is out of the model domain and apply
    606 !--          cyclic conditions if required
    607              ELSEIF ( cyclic_leg(l) )  THEN
    608 !
    609 !--             Check if aircraft leaves the model domain at the right boundary
    610                 IF ( ( flight_angle(l) >= 0.0_wp     .AND.                     &
    611                        flight_angle(l) <= 90.0_wp )  .OR.                      & 
    612                      ( flight_angle(l) >= 270.0_wp   .AND.                     &
    613                        flight_angle(l) <= 360.0_wp ) )  THEN
    614                    IF ( x_pos(l) >= ( nx + 0.5_wp ) * dx )                     &
    615                       x_pos(l) = x_pos(l) - ( nx + 1 ) * dx
    616 !
    617 !--             Check if aircraft leaves the model domain at the left boundary
    618                 ELSEIF ( flight_angle(l) > 90.0_wp  .AND.                      &
    619                          flight_angle(l) < 270.0_wp )  THEN
    620                    IF ( x_pos(l) < -0.5_wp * dx )                             &
    621                       x_pos(l) = ( nx + 1 ) * dx + x_pos(l) 
    622                 ENDIF
    623 !
    624 !--             Check if aircraft leaves the model domain at the north boundary
    625                 IF ( flight_angle(l) >= 0.0_wp  .AND.                          &
    626                      flight_angle(l) <= 180.0_wp )  THEN
    627                    IF ( y_pos(l) >= ( ny + 0.5_wp ) * dy )                     &
    628                       y_pos(l) = y_pos(l) - ( ny + 1 ) * dy
    629 !
    630 !--             Check if aircraft leaves the model domain at the south boundary
    631                 ELSEIF ( flight_angle(l) > 180.0_wp  .AND.                     &
    632                          flight_angle(l) < 360.0_wp )  THEN
    633                    IF ( y_pos(l) < -0.5_wp * dy )                              &
    634                       y_pos(l) = ( ny + 1 ) * dy + y_pos(l)
     633             ENDIF
     634
     635             IF ( y_start(l) <= y_end(l) )  THEN
     636!
     637!--             Outward flight, i.e. from start to end
     638                IF ( v_agl(l) >= 0.0_wp  .AND.  y_pos(l) > y_end(l)      )  THEN
     639                   y_pos(l) = 2.0_wp * y_end(l)   - y_pos(l)
     640                   v_agl(l) = - v_agl(l)
     641!
     642!--             Return flight, i.e. from end to start
     643                ELSEIF ( v_agl(l) < 0.0_wp  .AND.  y_pos(l) < y_start(l) )  THEN
     644                   y_pos(l) = 2.0_wp * y_start(l) - y_pos(l)
     645                   v_agl(l) = - v_agl(l)
    635646                ENDIF
    636                
     647             ELSE
     648!
     649!--             Outward flight, i.e. from start to end
     650                IF ( v_agl(l) < 0.0_wp  .AND.  y_pos(l) < y_end(l)      )  THEN
     651                   y_pos(l) = 2.0_wp * y_end(l)   - y_pos(l)
     652                   v_agl(l) = - v_agl(l)
     653!
     654!--             Return flight, i.e. from end to start
     655                ELSEIF ( v_agl(l) >= 0.0_wp  .AND.  y_pos(l) > y_start(l) )  THEN
     656                   y_pos(l) = 2.0_wp * y_start(l) - y_pos(l)
     657                   v_agl(l) = - v_agl(l)
     658                ENDIF
    637659             ENDIF
    638660!
    639 !--          Check if maximum elevation change is already reached. If required
    640 !--          reflect vertically.
    641              IF ( rate_of_climb(l) /= 0.0_wp )  THEN
    642 !
    643 !--             First check if aircraft is too high
    644                 IF (  w_agl(l) > 0.0_wp  .AND.                                 &
    645                       z_pos(l) - flight_level(l) > max_elev_change(l) )  THEN
    646                    z_pos(l) = 2.0_wp * ( flight_level(l) + max_elev_change(l) )&
    647                               - z_pos(l)
    648                    w_agl(l) = - w_agl(l)
    649 !
    650 !--             Check if aircraft is too low
    651                 ELSEIF (  w_agl(l) < 0.0_wp  .AND.  z_pos(l) < flight_level(l) )  THEN
    652                    z_pos(l) = 2.0_wp * flight_level(l) - z_pos(l)
    653                    w_agl(l) = - w_agl(l)
    654                 ENDIF
    655                
    656              ENDIF
    657 !
    658 !--          Determine grid indices for flight position along x- and y-direction.
    659 !--          Please note, there is a special treatment for the index
    660 !--          along z-direction, which is due to vertical grid stretching.
     661!--       Check if flight position is outside the model domain and apply cyclic conditions if required
     662          ELSEIF ( cyclic_leg(l) )  THEN
     663!
     664!--          Check if aircraft leaves the model domain at the right boundary
     665             IF ( ( flight_angle(l) >= 0.0_wp     .AND.                                            &
     666                    flight_angle(l) <= 90.0_wp )  .OR.                                             &
     667                  ( flight_angle(l) >= 270.0_wp   .AND.                                            &
     668                    flight_angle(l) <= 360.0_wp ) )  THEN
     669                IF ( x_pos(l) >= ( nx + 0.5_wp ) * dx )                                            &
     670                     x_pos(l) = x_pos(l) - ( nx + 1 ) * dx
     671!
     672!--          Check if aircraft leaves the model domain at the left boundary
     673             ELSEIF ( flight_angle(l) > 90.0_wp  .AND.  flight_angle(l) < 270.0_wp )  THEN
     674                IF ( x_pos(l) < -0.5_wp * dx )                                                     &
     675                     x_pos(l) = ( nx + 1 ) * dx + x_pos(l)
     676             ENDIF
     677!
     678!--          Check if aircraft leaves the model domain at the north boundary
     679             IF ( flight_angle(l) >= 0.0_wp  .AND.  flight_angle(l) <= 180.0_wp )  THEN
     680                IF ( y_pos(l) >= ( ny + 0.5_wp ) * dy )                                            &
     681                     y_pos(l) = y_pos(l) - ( ny + 1 ) * dy
     682!
     683!--          Check if aircraft leaves the model domain at the south boundary
     684             ELSEIF ( flight_angle(l) > 180.0_wp  .AND.  flight_angle(l) < 360.0_wp )  THEN
     685                IF ( y_pos(l) < -0.5_wp * dy )                                                     &
     686                     y_pos(l) = ( ny + 1 ) * dy + y_pos(l)
     687             ENDIF
     688
     689          ENDIF
     690!
     691!--       Check if maximum elevation change is already reached. If required reflect vertically.
     692          IF ( rate_of_climb(l) /= 0.0_wp )  THEN
     693!
     694!--          First check if aircraft is too high
     695             IF (  w_agl(l) > 0.0_wp  .AND.  z_pos(l) - flight_level(l) > max_elev_change(l) )  THEN
     696                z_pos(l) = 2.0_wp * ( flight_level(l) + max_elev_change(l) ) - z_pos(l)
     697                w_agl(l) = - w_agl(l)
     698!
     699!--          Check if aircraft is too low
     700             ELSEIF (  w_agl(l) < 0.0_wp  .AND.  z_pos(l) < flight_level(l) )  THEN
     701                z_pos(l) = 2.0_wp * flight_level(l) - z_pos(l)
     702                w_agl(l) = - w_agl(l)
     703             ENDIF
     704
     705          ENDIF
     706!
     707!--       Determine grid indices for flight position along x- and y-direction. Please note, there is
     708!--       a special treatment for the index along z-direction, which is due to vertical grid stretching.
     709          i = ( x_pos(l) + 0.5_wp * dx ) * ddx
     710          j = ( y_pos(l) + 0.5_wp * dy ) * ddy
     711!
     712!--       Check if indices are on current PE
     713          on_pe = ( i >= nxl  .AND.  i <= nxr  .AND.  j >= nys  .AND.  j <= nyn )
     714
     715          IF ( on_pe )  THEN
     716
     717             var_index = 1
     718!
     719!--          Recalculate indices, as u is shifted by -0.5*dx.
     720             i =   x_pos(l) * ddx
     721             j = ( y_pos(l) + 0.5_wp * dy ) * ddy
     722!
     723!--          Calculate distance from left and south grid-box edges.
     724             x  = x_pos(l) - ( 0.5_wp - i ) * dx
     725             y  = y_pos(l) - j * dy
     726!
     727!--          Interpolate u-component onto current flight position.
     728             CALL interpolate_xyz( u, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
     729             var_index = var_index + 1
     730!
     731!--          Recalculate indices, as v is shifted by -0.5*dy.
    661732             i = ( x_pos(l) + 0.5_wp * dx ) * ddx
    662              j = ( y_pos(l) + 0.5_wp * dy ) * ddy
    663 !
    664 !--          Check if indices are on current PE
    665              on_pe = ( i >= nxl  .AND.  i <= nxr  .AND.                        &
    666                        j >= nys  .AND.  j <= nyn )
    667 
    668              IF ( on_pe )  THEN
    669 
    670                 var_index = 1
    671 !
    672 !--             Recalculate indices, as u is shifted by -0.5*dx.
    673                 i =   x_pos(l) * ddx
    674                 j = ( y_pos(l) + 0.5_wp * dy ) * ddy
    675 !
    676 !--             Calculate distance from left and south grid-box edges.
    677                 x  = x_pos(l) - ( 0.5_wp - i ) * dx
    678                 y  = y_pos(l) - j * dy
    679 !
    680 !--             Interpolate u-component onto current flight position.
    681                 CALL interpolate_xyz( u, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
     733             j =   y_pos(l) * ddy
     734
     735             x  = x_pos(l) - i * dx
     736             y  = y_pos(l) - ( 0.5_wp - j ) * dy
     737             CALL interpolate_xyz( v, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
     738             var_index = var_index + 1
     739!
     740!--          Interpolate w and scalar quantities. Recalculate indices.
     741             i  = ( x_pos(l) + 0.5_wp * dx ) * ddx
     742             j  = ( y_pos(l) + 0.5_wp * dy ) * ddy
     743             x  = x_pos(l) - i * dx
     744             y  = y_pos(l) - j * dy
     745!
     746!--          Interpolate w-velocity component.
     747             CALL interpolate_xyz( w, zw, ddzw, 0.0_wp, x, y, var_index, j, i )
     748             var_index = var_index + 1
     749!
     750!--          Potential temerature
     751             IF ( .NOT. neutral )  THEN
     752                CALL interpolate_xyz( pt, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
    682753                var_index = var_index + 1
    683 !
    684 !--             Recalculate indices, as v is shifted by -0.5*dy.
    685                 i = ( x_pos(l) + 0.5_wp * dx ) * ddx
    686                 j =   y_pos(l) * ddy
    687 
    688                 x  = x_pos(l) - i * dx
    689                 y  = y_pos(l) - ( 0.5_wp - j ) * dy
    690                 CALL interpolate_xyz( v, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
     754             ENDIF
     755!
     756!--          Humidity
     757             IF ( humidity )  THEN
     758                CALL interpolate_xyz( q, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
    691759                var_index = var_index + 1
    692 !
    693 !--             Interpolate w and scalar quantities. Recalculate indices.
    694                 i  = ( x_pos(l) + 0.5_wp * dx ) * ddx
    695                 j  = ( y_pos(l) + 0.5_wp * dy ) * ddy
    696                 x  = x_pos(l) - i * dx
    697                 y  = y_pos(l) - j * dy
    698 !
    699 !--             Interpolate w-velocity component.
    700                 CALL interpolate_xyz( w, zw, ddzw, 0.0_wp, x, y, var_index, j, i )
     760             ENDIF
     761!
     762!--          Liquid water content
     763             IF ( bulk_cloud_model .OR. cloud_droplets )  THEN
     764                CALL interpolate_xyz( ql, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
    701765                var_index = var_index + 1
    702 !
    703 !--             Potential temerature
    704                 IF ( .NOT. neutral )  THEN
    705                    CALL interpolate_xyz( pt, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
    706                    var_index = var_index + 1
    707                 ENDIF
    708 !
    709 !--             Humidity
    710                 IF ( humidity )  THEN
    711                    CALL interpolate_xyz( q, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
    712                    var_index = var_index + 1
    713                 ENDIF
    714 !
    715 !--             Liquid water content
    716                 IF ( bulk_cloud_model .OR. cloud_droplets )  THEN
    717                    CALL interpolate_xyz( ql, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
    718                    var_index = var_index + 1
    719                 ENDIF
    720 !
    721 !--             Passive scalar
    722                 IF ( passive_scalar )  THEN
    723                    CALL interpolate_xyz( s, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
    724                    var_index = var_index + 1
    725                 ENDIF
    726 !
    727 !--             Treat user-defined variables if required
    728                 DO n = 1, num_var_fl_user               
    729                    CALL user_flight( var_u, n )
    730                    CALL interpolate_xyz( var_u, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
    731                    var_index = var_index + 1
    732                 ENDDO
    733766             ENDIF
    734 
    735           ENDDO
    736 !
    737 !--       Write local data on global array.
     767!
     768!--          Passive scalar
     769             IF ( passive_scalar )  THEN
     770                CALL interpolate_xyz( s, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
     771                var_index = var_index + 1
     772             ENDIF
     773!
     774!--          Treat user-defined variables if required
     775             DO  n = 1, num_var_fl_user
     776                CALL user_flight( var_u, n )
     777                CALL interpolate_xyz( var_u, zu, ddzu, 1.0_wp, x, y, var_index, j, i )
     778                var_index = var_index + 1
     779             ENDDO
     780          ENDIF
     781
     782       ENDDO
     783!
     784!--    Write local data on global array.
    738785#if defined( __parallel )
    739           CALL MPI_ALLREDUCE(sensor_l(1,1), sensor(1,1),                       &
    740                              num_var_fl*num_leg, MPI_REAL, MPI_SUM,               &
    741                              comm2d, ierr )
     786       CALL MPI_ALLREDUCE( sensor_l(1,1), sensor(1,1), num_var_fl * num_leg, MPI_REAL, MPI_SUM,    &
     787                           comm2d, ierr )
    742788#else
    743           sensor     = sensor_l
     789       sensor     = sensor_l
    744790#endif
    745        ENDIF
    746        
    747        CALL cpu_log( log_point(65), 'flight_measurement', 'stop' )
    748 
    749     END SUBROUTINE flight_measurement
    750 
    751 !------------------------------------------------------------------------------!
     791    ENDIF
     792
     793    CALL cpu_log( log_point(65), 'flight_measurement', 'stop' )
     794
     795 END SUBROUTINE flight_measurement
     796
     797!--------------------------------------------------------------------------------------------------!
    752798! Description:
    753799! ------------
    754 !> This subroutine bi-linearly interpolates the respective data onto the current
    755 !> flight position.
    756 !------------------------------------------------------------------------------!
    757     SUBROUTINE interpolate_xyz( var, z_uw, ddz_uw, fac, x, y, var_ind, j, i )
    758 
    759        USE control_parameters,                                                 &
    760            ONLY:  dz, dz_stretch_level_start   
    761  
    762       USE grid_variables,                                                      &
    763           ONLY:  dx, dy
    764    
    765        USE indices,                                                            &
    766            ONLY:  nzb, nzt, nxlg, nxrg, nysg, nyng
    767 
    768        IMPLICIT NONE
    769 
    770        INTEGER(iwp) ::  i        !< index along x
    771        INTEGER(iwp) ::  j        !< index along y
    772        INTEGER(iwp) ::  k        !< index along z
    773        INTEGER(iwp) ::  k1       !< dummy variable
    774        INTEGER(iwp) ::  var_ind  !< index variable for output quantity
    775 
    776        REAL(wp) ::  aa        !< dummy argument for horizontal interpolation   
    777        REAL(wp) ::  bb        !< dummy argument for horizontal interpolation
    778        REAL(wp) ::  cc        !< dummy argument for horizontal interpolation
    779        REAL(wp) ::  dd        !< dummy argument for horizontal interpolation
    780        REAL(wp) ::  gg        !< dummy argument for horizontal interpolation
    781        REAL(wp) ::  fac       !< flag to indentify if quantity is on zu or zw level
    782        REAL(wp) ::  var_int   !< horizontal interpolated variable at current position
    783        REAL(wp) ::  var_int_l !< horizontal interpolated variable at k-level
    784        REAL(wp) ::  var_int_u !< horizontal interpolated variable at (k+1)-level
    785        REAL(wp) ::  x         !< distance between left edge of current grid box and flight position
    786        REAL(wp) ::  y         !< distance between south edge of current grid box and flight position
    787 
    788        REAL(wp), DIMENSION(1:nzt+1)   ::  ddz_uw !< inverse vertical grid spacing
    789        REAL(wp), DIMENSION(nzb:nzt+1) ::  z_uw   !< height level
    790 
    791        REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var !< treted quantity
    792 !
    793 !--    Calculate interpolation coefficients
    794        aa = x**2          + y**2
    795        bb = ( dx - x )**2 + y**2
    796        cc = x**2          + ( dy - y )**2
    797        dd = ( dx - x )**2 + ( dy - y )**2
    798        gg = aa + bb + cc + dd
    799 !
    800 !--    Obtain vertical index. Special treatment for grid index along z-direction
    801 !--    if flight position is above the vertical grid-stretching level.
    802 !--    fac=1 if variable is on scalar grid level, fac=0 for w-component.
    803        IF ( z_pos(l) < dz_stretch_level_start(1) )  THEN
    804           k = ( z_pos(l) + fac * 0.5_wp * dz(1) ) / dz(1)
    805        ELSE
    806 !
    807 !--       Search for k-index
    808           DO k1=nzb, nzt
    809              IF ( z_pos(l) >= z_uw(k1) .AND. z_pos(l) < z_uw(k1+1) )  THEN
    810                 k = k1
    811                 EXIT
    812              ENDIF                   
    813           ENDDO
    814        ENDIF
    815 !
    816 !--    (x,y)-interpolation at k-level
    817        var_int_l = ( ( gg - aa ) * var(k,j,i)       +                          &
    818                      ( gg - bb ) * var(k,j,i+1)     +                          &
    819                      ( gg - cc ) * var(k,j+1,i)     +                          &
    820                      ( gg - dd ) * var(k,j+1,i+1)                              &
    821                    ) / ( 3.0_wp * gg )
    822 !
    823 !--    (x,y)-interpolation on (k+1)-level
    824        var_int_u = ( ( gg - aa ) * var(k+1,j,i)     +                          &
    825                      ( gg - bb ) * var(k+1,j,i+1)   +                          &
    826                      ( gg - cc ) * var(k+1,j+1,i)   +                          &
    827                      ( gg - dd ) * var(k+1,j+1,i+1)                            &
    828                    ) / ( 3.0_wp * gg )
    829 !
    830 !--    z-interpolation onto current flight postion
    831        var_int = var_int_l                                                     &
    832                            + ( z_pos(l) - z_uw(k) ) * ddz_uw(k+1)              &
    833                            * (var_int_u - var_int_l )
    834 !
    835 !--    Store on local data array
    836        sensor_l(var_ind,l) = var_int
    837 
    838     END SUBROUTINE interpolate_xyz
    839 
    840 !------------------------------------------------------------------------------!
     800!> This subroutine bi-linearly interpolates the respective data onto the current flight position.
     801!--------------------------------------------------------------------------------------------------!
     802 SUBROUTINE interpolate_xyz( var, z_uw, ddz_uw, fac, x, y, var_ind, j, i )
     803
     804    USE control_parameters,                                                                        &
     805        ONLY:  dz,                                                                                 &
     806               dz_stretch_level_start
     807
     808    USE grid_variables,                                                                            &
     809       ONLY:  dx,                                                                                  &
     810              dy
     811
     812    USE indices,                                                                                   &
     813        ONLY:  nzb,                                                                                &
     814               nzt,                                                                                &
     815               nxlg,                                                                               &
     816               nxrg,                                                                               &
     817               nysg,                                                                               &
     818               nyng
     819
     820    IMPLICIT NONE
     821
     822    INTEGER(iwp) ::  i        !< index along x
     823    INTEGER(iwp) ::  j        !< index along y
     824    INTEGER(iwp) ::  k        !< index along z
     825    INTEGER(iwp) ::  k1       !< dummy variable
     826    INTEGER(iwp) ::  var_ind  !< index variable for output quantity
     827
     828    REAL(wp) ::  aa         !< dummy argument for horizontal interpolation
     829    REAL(wp) ::  bb         !< dummy argument for horizontal interpolation
     830    REAL(wp) ::  cc         !< dummy argument for horizontal interpolation
     831    REAL(wp) ::  dd         !< dummy argument for horizontal interpolation
     832    REAL(wp) ::  gg         !< dummy argument for horizontal interpolation
     833    REAL(wp) ::  fac        !< flag to indentify if quantity is on zu or zw level
     834    REAL(wp) ::  var_int    !< horizontal interpolated variable at current position
     835    REAL(wp) ::  var_int_l  !< horizontal interpolated variable at k-level
     836    REAL(wp) ::  var_int_u  !< horizontal interpolated variable at (k+1)-level
     837    REAL(wp) ::  x          !< distance between left edge of current grid box and flight position
     838    REAL(wp) ::  y          !< distance between south edge of current grid box and flight position
     839
     840    REAL(wp), DIMENSION(1:nzt+1)   ::  ddz_uw  !< inverse vertical grid spacing
     841    REAL(wp), DIMENSION(nzb:nzt+1) ::  z_uw    !< height level
     842
     843    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var  !< treated quantity
     844!
     845!-- Calculate interpolation coefficients
     846    aa = x**2          + y**2
     847    bb = ( dx - x )**2 + y**2
     848    cc = x**2          + ( dy - y )**2
     849    dd = ( dx - x )**2 + ( dy - y )**2
     850    gg = aa + bb + cc + dd
     851!
     852!-- Obtain vertical index. Special treatment for grid index along z-direction if flight position is
     853!-- above the vertical grid-stretching level. fac=1 if variable is on scalar grid level, fac=0 for
     854!-- w-component.
     855    IF ( z_pos(l) < dz_stretch_level_start(1) )  THEN
     856       k = ( z_pos(l) + fac * 0.5_wp * dz(1) ) / dz(1)
     857    ELSE
     858!
     859!--    Search for k-index
     860       DO  k1 = nzb, nzt
     861          IF ( z_pos(l) >= z_uw(k1) .AND. z_pos(l) < z_uw(k1+1) )  THEN
     862             k = k1
     863             EXIT
     864          ENDIF
     865       ENDDO
     866    ENDIF
     867!
     868!-- (x,y)-interpolation at k-level
     869    var_int_l = ( ( gg - aa ) * var(k,j,i)       +                                                 &
     870                  ( gg - bb ) * var(k,j,i+1)     +                                                 &
     871                  ( gg - cc ) * var(k,j+1,i)     +                                                 &
     872                  ( gg - dd ) * var(k,j+1,i+1)                                                     &
     873                ) / ( 3.0_wp * gg )
     874!
     875!-- (x,y)-interpolation on (k+1)-level
     876    var_int_u = ( ( gg - aa ) * var(k+1,j,i)     +                                                 &
     877                  ( gg - bb ) * var(k+1,j,i+1)   +                                                 &
     878                  ( gg - cc ) * var(k+1,j+1,i)   +                                                 &
     879                  ( gg - dd ) * var(k+1,j+1,i+1)                                                   &
     880                ) / ( 3.0_wp * gg )
     881!
     882!-- z-interpolation onto current flight postion
     883    var_int = var_int_l + ( z_pos(l) - z_uw(k) ) * ddz_uw(k+1) * (var_int_u - var_int_l )
     884!
     885!-- Store on local data array
     886    sensor_l(var_ind,l) = var_int
     887
     888 END SUBROUTINE interpolate_xyz
     889
     890!--------------------------------------------------------------------------------------------------!
    841891! Description:
    842892! ------------
    843893!> Perform parameter checks.
    844 !------------------------------------------------------------------------------!
    845     SUBROUTINE flight_check_parameters
    846 
    847        USE arrays_3d,                                                          &
    848            ONLY:  zu
    849    
    850        USE control_parameters,                                                 &
    851            ONLY:  bc_lr_cyc, bc_ns_cyc, message_string
    852 
    853        USE grid_variables,                                                     &
    854            ONLY:  dx, dy   
    855          
    856        USE indices,                                                            &
    857            ONLY:  nx, ny, nz
    858            
    859        USE netcdf_interface,                                                   &
    860            ONLY:  netcdf_data_format
    861 
    862        IMPLICIT NONE
    863        
    864 !
    865 !--    Check if start positions are properly set.
    866        DO l=1, num_leg
    867           IF ( x_start(l) < 0.0_wp  .OR.  x_start(l) > ( nx + 1 ) * dx )  THEN
    868              message_string = 'Start x position is outside the model domain'
    869              CALL message( 'flight_check_parameters', 'PA0431', 1, 2, 0, 6, 0 )
    870           ENDIF
    871           IF ( y_start(l) < 0.0_wp  .OR.  y_start(l) > ( ny + 1 ) * dy )  THEN
    872              message_string = 'Start y position is outside the model domain'
    873              CALL message( 'flight_check_parameters', 'PA0432', 1, 2, 0, 6, 0 )
    874           ENDIF
    875 
    876        ENDDO
    877 !
    878 !--    Check for leg mode
    879        DO l=1, num_leg
    880 !
    881 !--       Check if leg mode matches the overall lateral model boundary
    882 !--       conditions.
    883           IF ( TRIM( leg_mode(l) ) == 'cyclic' )  THEN
    884              IF ( .NOT. bc_lr_cyc  .OR.  .NOT. bc_ns_cyc )  THEN
    885                 message_string = 'Cyclic flight leg does not match ' //        &
    886                                  'lateral boundary condition'
    887                 CALL message( 'flight_check_parameters', 'PA0433', 1, 2, 0, 6, 0 )
    888              ENDIF
    889 !
    890 !--       Check if end-positions are inside the model domain in case of
    891 !..       return-legs. 
    892           ELSEIF ( TRIM( leg_mode(l) ) == 'return' )  THEN
    893              IF ( x_end(l) > ( nx + 1 ) * dx  .OR.                             &
    894                   y_end(l) > ( ny + 1 ) * dx )  THEN
    895                 message_string = 'Flight leg or parts of it are outside ' //   &
    896                                  'the model domain'
    897                 CALL message( 'flight_check_parameters', 'PA0434', 1, 2, 0, 6, 0 )
    898              ENDIF
    899           ELSE
    900              message_string = 'Unknown flight mode'
    901              CALL message( 'flight_check_parameters', 'PA0435', 1, 2, 0, 6, 0 )
     894!--------------------------------------------------------------------------------------------------!
     895 SUBROUTINE flight_check_parameters
     896
     897    USE arrays_3d,                                                                                 &
     898        ONLY:  zu
     899
     900    USE control_parameters,                                                                        &
     901        ONLY:  bc_lr_cyc,                                                                          &
     902               bc_ns_cyc,                                                                          &
     903               message_string
     904
     905    USE grid_variables,                                                                            &
     906        ONLY:  dx,                                                                                 &
     907               dy
     908
     909    USE indices,                                                                                   &
     910        ONLY:  nx,                                                                                 &
     911               ny,                                                                                 &
     912               nz
     913
     914    USE netcdf_interface,                                                                          &
     915        ONLY:  netcdf_data_format
     916
     917    IMPLICIT NONE
     918
     919!
     920!-- Check if start positions are properly set.
     921    DO  l = 1, num_leg
     922       IF ( x_start(l) < 0.0_wp  .OR.  x_start(l) > ( nx + 1 ) * dx )  THEN
     923          message_string = 'Start x position is outside the model domain'
     924          CALL message( 'flight_check_parameters', 'PA0431', 1, 2, 0, 6, 0 )
     925       ENDIF
     926       IF ( y_start(l) < 0.0_wp  .OR.  y_start(l) > ( ny + 1 ) * dy )  THEN
     927          message_string = 'Start y position is outside the model domain'
     928          CALL message( 'flight_check_parameters', 'PA0432', 1, 2, 0, 6, 0 )
     929       ENDIF
     930
     931    ENDDO
     932!
     933!-- Check for leg mode
     934    DO  l = 1, num_leg
     935!
     936!--    Check if leg mode matches the overall lateral model boundary conditions.
     937       IF ( TRIM( leg_mode(l) ) == 'cyclic' )  THEN
     938          IF ( .NOT. bc_lr_cyc  .OR.  .NOT. bc_ns_cyc )  THEN
     939             message_string = 'Cyclic flight leg does not match lateral boundary condition'
     940             CALL message( 'flight_check_parameters', 'PA0433', 1, 2, 0, 6, 0 )
    902941          ENDIF
    903 
    904        ENDDO         
    905 !
    906 !--    Check if given flight object remains inside model domain if a rate of
    907 !--    climb / descent is prescribed.
    908        DO l=1, num_leg
    909           IF ( flight_level(l) + max_elev_change(l) > zu(nz) .OR.              &
    910                flight_level(l) + max_elev_change(l) <= 0.0_wp )  THEN
    911              message_string = 'Flight level is outside the model domain '
    912              CALL message( 'flight_check_parameters', 'PA0438', 1, 2, 0, 6, 0 )
    913           ENDIF
    914        ENDDO       
    915 !
    916 !--    Check for appropriate NetCDF format. Definition of more than one
    917 !--    unlimited dimension is unfortunately only possible with NetCDF4/HDF5.
    918        IF (  netcdf_data_format <= 2 )  THEN
    919           message_string = 'netcdf_data_format must be > 2'
    920           CALL message( 'flight_check_parameters', 'PA0439', 1, 2, 0, 6, 0 )
     942!
     943!--    Check if end-positions are inside the model domain in case of return-legs.
     944       ELSEIF ( TRIM( leg_mode(l) ) == 'return' )  THEN
     945          IF ( x_end(l) > ( nx + 1 ) * dx  .OR.  y_end(l) > ( ny + 1 ) * dx )  THEN
     946             message_string = 'Flight leg or parts of it are outside the model domain'
     947             CALL message( 'flight_check_parameters', 'PA0434', 1, 2, 0, 6, 0 )
     948          ENDIF
     949       ELSE
     950          message_string = 'Unknown flight mode'
     951          CALL message( 'flight_check_parameters', 'PA0435', 1, 2, 0, 6, 0 )
    921952       ENDIF
    922953
    923 
    924     END SUBROUTINE flight_check_parameters
    925        
    926 
    927 !------------------------------------------------------------------------------!
     954    ENDDO
     955!
     956!-- Check if given flight object remains inside model domain if a rate of climb / descent is
     957!-- prescribed.
     958    DO  l = 1, num_leg
     959       IF ( flight_level(l) + max_elev_change(l) > zu(nz) .OR.                                     &
     960            flight_level(l) + max_elev_change(l) <= 0.0_wp )  THEN
     961          message_string = 'Flight level is outside the model domain '
     962          CALL message( 'flight_check_parameters', 'PA0438', 1, 2, 0, 6, 0 )
     963       ENDIF
     964    ENDDO
     965!
     966!-- Check for appropriate NetCDF format. Definition of more than one unlimited dimension is
     967!-- unfortunately only possible with NetCDF4/HDF5.
     968    IF (  netcdf_data_format <= 2 )  THEN
     969       message_string = 'netcdf_data_format must be > 2'
     970       CALL message( 'flight_check_parameters', 'PA0439', 1, 2, 0, 6, 0 )
     971    ENDIF
     972
     973
     974 END SUBROUTINE flight_check_parameters
     975
     976
     977!--------------------------------------------------------------------------------------------------!
    928978! Description:
    929979! ------------
    930980!> Read module-specific global restart data (Fortran binary format).
    931 !------------------------------------------------------------------------------!
    932     SUBROUTINE flight_rrd_global_ftn( found )
    933 
    934 
    935        USE control_parameters,                                                 &
    936            ONLY: length, restart_string
    937 
    938    
    939        IMPLICIT NONE
    940        
    941        LOGICAL, INTENT(OUT)  ::  found !< flag indicating if a variable string is found
    942 
    943 
    944        found = .TRUE.
    945 
    946 
    947        SELECT CASE ( restart_string(1:length) )
    948          
    949           CASE ( 'u_agl' )
    950              IF ( .NOT. ALLOCATED( u_agl ) )  ALLOCATE( u_agl(1:num_leg) )
    951              READ ( 13 )  u_agl   
    952           CASE ( 'v_agl' )
    953              IF ( .NOT. ALLOCATED( v_agl ) )  ALLOCATE( v_agl(1:num_leg) )
    954              READ ( 13 )  v_agl
    955           CASE ( 'w_agl' )
    956              IF ( .NOT. ALLOCATED( w_agl ) )  ALLOCATE( w_agl(1:num_leg) )
    957              READ ( 13 )  w_agl
    958           CASE ( 'x_pos' )
    959              IF ( .NOT. ALLOCATED( x_pos ) )  ALLOCATE( x_pos(1:num_leg) )
    960              READ ( 13 )  x_pos
    961           CASE ( 'y_pos' )
    962              IF ( .NOT. ALLOCATED( y_pos ) )  ALLOCATE( y_pos(1:num_leg) )
    963              READ ( 13 )  y_pos
    964           CASE ( 'z_pos' )
    965              IF ( .NOT. ALLOCATED( z_pos ) )  ALLOCATE( z_pos(1:num_leg) )
    966              READ ( 13 )  z_pos
    967 
    968           CASE DEFAULT
    969 
    970              found = .FALSE.
    971          
    972        END SELECT
    973 
    974 
    975     END SUBROUTINE flight_rrd_global_ftn
    976 
    977 
     981!--------------------------------------------------------------------------------------------------!
     982 SUBROUTINE flight_rrd_global_ftn( found )
     983
     984
     985    USE control_parameters,                                                                        &
     986        ONLY: length,                                                                              &
     987              restart_string
     988
     989
     990    IMPLICIT NONE
     991
     992    LOGICAL, INTENT(OUT)  ::  found  !< flag indicating if a variable string is found
     993
     994
     995    found = .TRUE.
     996
     997
     998    SELECT CASE ( restart_string(1:length) )
     999
     1000       CASE ( 'u_agl' )
     1001          IF ( .NOT. ALLOCATED( u_agl ) )  ALLOCATE( u_agl(1:num_leg) )
     1002          READ ( 13 )  u_agl
     1003       CASE ( 'v_agl' )
     1004          IF ( .NOT. ALLOCATED( v_agl ) )  ALLOCATE( v_agl(1:num_leg) )
     1005          READ ( 13 )  v_agl
     1006       CASE ( 'w_agl' )
     1007          IF ( .NOT. ALLOCATED( w_agl ) )  ALLOCATE( w_agl(1:num_leg) )
     1008          READ ( 13 )  w_agl
     1009       CASE ( 'x_pos' )
     1010          IF ( .NOT. ALLOCATED( x_pos ) )  ALLOCATE( x_pos(1:num_leg) )
     1011          READ ( 13 )  x_pos
     1012       CASE ( 'y_pos' )
     1013          IF ( .NOT. ALLOCATED( y_pos ) )  ALLOCATE( y_pos(1:num_leg) )
     1014          READ ( 13 )  y_pos
     1015       CASE ( 'z_pos' )
     1016          IF ( .NOT. ALLOCATED( z_pos ) )  ALLOCATE( z_pos(1:num_leg) )
     1017          READ ( 13 )  z_pos
     1018
     1019       CASE DEFAULT
     1020
     1021          found = .FALSE.
     1022
     1023    END SELECT
     1024
     1025
     1026 END SUBROUTINE flight_rrd_global_ftn
     1027
     1028
     1029 
    9781030!------------------------------------------------------------------------------!
    9791031! Description:
     
    9811033!> Read module-specific global restart data (MPI-IO).
    9821034!------------------------------------------------------------------------------!
    983     SUBROUTINE flight_rrd_global_mpi
    984 
    985 
    986        IMPLICIT NONE
    987 
    988        LOGICAL  ::  array_found  !< flag indicating if respective array is found in restart file
    989 
    990 
    991        CALL rd_mpi_io_check_array( 'u_agl', found = array_found )
    992        IF ( array_found)  THEN
    993           IF ( .NOT. ALLOCATED( u_agl ) )  ALLOCATE( u_agl(1:num_leg) )
    994           CALL rrd_mpi_io_global_array( 'u_agl', u_agl )
    995        ENDIF
    996        CALL rd_mpi_io_check_array( 'v_agl', found = array_found )
    997        IF ( array_found)  THEN
    998           IF ( .NOT. ALLOCATED( v_agl ) )  ALLOCATE( v_agl(1:num_leg) )
    999           CALL rrd_mpi_io_global_array( 'v_agl', v_agl )
    1000        ENDIF
    1001        CALL rd_mpi_io_check_array( 'w_agl', found = array_found )
    1002        IF ( array_found)  THEN
    1003           IF ( .NOT. ALLOCATED( w_agl ) )  ALLOCATE( w_agl(1:num_leg) )
    1004           CALL rrd_mpi_io_global_array( 'w_agl', w_agl )
    1005        ENDIF
    1006        CALL rd_mpi_io_check_array( 'x_pos', found = array_found )
    1007        IF ( array_found)  THEN
    1008           IF ( .NOT. ALLOCATED( x_pos ) )  ALLOCATE( x_pos(1:num_leg) )
    1009           CALL rrd_mpi_io_global_array( 'x_pos', x_pos )
    1010        ENDIF
    1011        CALL rd_mpi_io_check_array( 'y_pos', found = array_found )
    1012        IF ( array_found)  THEN
    1013           IF ( .NOT. ALLOCATED( y_pos ) )  ALLOCATE( y_pos(1:num_leg) )
    1014           CALL rrd_mpi_io_global_array( 'y_pos', y_pos )
    1015        ENDIF
    1016        CALL rd_mpi_io_check_array( 'z_pos', found = array_found )
    1017        IF ( array_found)  THEN
    1018           IF ( .NOT. ALLOCATED( z_pos ) )  ALLOCATE( z_pos(1:num_leg) )
    1019           CALL rrd_mpi_io_global_array( 'z_pos', z_pos )
    1020        ENDIF
    1021 
    1022     END SUBROUTINE flight_rrd_global_mpi
     1035 SUBROUTINE flight_rrd_global_mpi
     1036
     1037
     1038    IMPLICIT NONE
     1039
     1040    LOGICAL  ::  array_found  !< flag indicating if respective array is found in restart file
     1041
     1042
     1043    CALL rd_mpi_io_check_array( 'u_agl', found = array_found )
     1044    IF ( array_found)  THEN
     1045       IF ( .NOT. ALLOCATED( u_agl ) )  ALLOCATE( u_agl(1:num_leg) )
     1046       CALL rrd_mpi_io_global_array( 'u_agl', u_agl )
     1047    ENDIF
     1048    CALL rd_mpi_io_check_array( 'v_agl', found = array_found )
     1049    IF ( array_found)  THEN
     1050       IF ( .NOT. ALLOCATED( v_agl ) )  ALLOCATE( v_agl(1:num_leg) )
     1051       CALL rrd_mpi_io_global_array( 'v_agl', v_agl )
     1052    ENDIF
     1053    CALL rd_mpi_io_check_array( 'w_agl', found = array_found )
     1054    IF ( array_found)  THEN
     1055       IF ( .NOT. ALLOCATED( w_agl ) )  ALLOCATE( w_agl(1:num_leg) )
     1056       CALL rrd_mpi_io_global_array( 'w_agl', w_agl )
     1057    ENDIF
     1058    CALL rd_mpi_io_check_array( 'x_pos', found = array_found )
     1059    IF ( array_found)  THEN
     1060       IF ( .NOT. ALLOCATED( x_pos ) )  ALLOCATE( x_pos(1:num_leg) )
     1061       CALL rrd_mpi_io_global_array( 'x_pos', x_pos )
     1062    ENDIF
     1063    CALL rd_mpi_io_check_array( 'y_pos', found = array_found )
     1064    IF ( array_found)  THEN
     1065       IF ( .NOT. ALLOCATED( y_pos ) )  ALLOCATE( y_pos(1:num_leg) )
     1066       CALL rrd_mpi_io_global_array( 'y_pos', y_pos )
     1067    ENDIF
     1068    CALL rd_mpi_io_check_array( 'z_pos', found = array_found )
     1069    IF ( array_found)  THEN
     1070       IF ( .NOT. ALLOCATED( z_pos ) )  ALLOCATE( z_pos(1:num_leg) )
     1071       CALL rrd_mpi_io_global_array( 'z_pos', z_pos )
     1072    ENDIF
     1073
     1074 END SUBROUTINE flight_rrd_global_mpi
    10231075
    10241076   
    1025 !------------------------------------------------------------------------------!
     1077    !--------------------------------------------------------------------------------------------------!
    10261078! Description:
    10271079! ------------
    10281080!> This routine writes the respective restart data.
    1029 !------------------------------------------------------------------------------!
    1030     SUBROUTINE flight_wrd_global 
    1031 
    1032 
    1033        IMPLICIT NONE
    1034  
    1035 
    1036        IF ( TRIM( restart_data_format_output ) == 'fortran_binary' )  THEN
    1037 
    1038           CALL wrd_write_string( 'u_agl' )
    1039           WRITE ( 14 )  u_agl
    1040 
    1041           CALL wrd_write_string( 'v_agl' )
    1042           WRITE ( 14 )  v_agl
    1043 
    1044           CALL wrd_write_string( 'w_agl' )
    1045           WRITE ( 14 )  w_agl
    1046 
    1047           CALL wrd_write_string( 'x_pos' )
    1048           WRITE ( 14 )  x_pos
    1049 
    1050           CALL wrd_write_string( 'y_pos' )
    1051           WRITE ( 14 )  y_pos
    1052 
    1053           CALL wrd_write_string( 'z_pos' )
    1054           WRITE ( 14 )  z_pos
    1055 
    1056        ELSEIF ( TRIM( restart_data_format_output ) == 'mpi' )  THEN
    1057 
    1058           CALL wrd_mpi_io_global_array( 'u_agl', u_agl )
    1059           CALL wrd_mpi_io_global_array( 'v_agl', v_agl )
    1060           CALL wrd_mpi_io_global_array( 'w_agl', w_agl )
    1061           CALL wrd_mpi_io_global_array( 'x_pos', x_pos )
    1062           CALL wrd_mpi_io_global_array( 'y_pos', y_pos )
    1063           CALL wrd_mpi_io_global_array( 'z_pos', z_pos )
    1064 
    1065        ENDIF
    1066 
    1067     END SUBROUTINE flight_wrd_global   
    1068    
     1081!--------------------------------------------------------------------------------------------------!
     1082 SUBROUTINE flight_wrd_global
     1083
     1084
     1085    IMPLICIT NONE
     1086
     1087
     1088    IF ( TRIM( restart_data_format_output ) == 'fortran_binary' )  THEN
     1089
     1090       CALL wrd_write_string( 'u_agl' )
     1091       WRITE ( 14 )  u_agl
     1092
     1093       CALL wrd_write_string( 'v_agl' )
     1094       WRITE ( 14 )  v_agl
     1095
     1096       CALL wrd_write_string( 'w_agl' )
     1097       WRITE ( 14 )  w_agl
     1098
     1099       CALL wrd_write_string( 'x_pos' )
     1100       WRITE ( 14 )  x_pos
     1101
     1102       CALL wrd_write_string( 'y_pos' )
     1103       WRITE ( 14 )  y_pos
     1104
     1105       CALL wrd_write_string( 'z_pos' )
     1106       WRITE ( 14 )  z_pos
     1107
     1108    ELSEIF ( TRIM( restart_data_format_output ) == 'mpi' )  THEN
     1109
     1110       CALL wrd_mpi_io_global_array( 'u_agl', u_agl )
     1111       CALL wrd_mpi_io_global_array( 'v_agl', v_agl )
     1112       CALL wrd_mpi_io_global_array( 'w_agl', w_agl )
     1113       CALL wrd_mpi_io_global_array( 'x_pos', x_pos )
     1114       CALL wrd_mpi_io_global_array( 'y_pos', y_pos )
     1115       CALL wrd_mpi_io_global_array( 'z_pos', z_pos )
     1116
     1117    ENDIF
     1118
     1119 END SUBROUTINE flight_wrd_global
     1120
    10691121
    10701122 END MODULE flight_mod
  • palm/trunk/SOURCE/wind_turbine_model_mod.f90

    r4495 r4497  
    11!> @file wind_turbine_model_mod.f90
    2 !------------------------------------------------------------------------------!
     2!--------------------------------------------------------------------------------------------------!
    33! This file is part of the PALM model system.
    44!
    5 ! PALM is free software: you can redistribute it and/or modify it under the
    6 ! terms of the GNU General Public License as published by the Free Software
    7 ! Foundation, either version 3 of the License, or (at your option) any later
    8 ! version.
    9 !
    10 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANYr
    11 ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
    12 ! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
    13 !
    14 ! You should have received a copy of the GNU General Public License along with
    15 ! PALM. If not, see <http://www.gnu.org/licenses/>.
     5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
     6! Public License as published by the Free Software Foundation, either version 3 of the License, or
     7! (at your option) any later version.
     8!
     9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
     10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
     11! Public License for more details.
     12!
     13! You should have received a copy of the GNU General Public License along with PALM. If not, see
     14! <http://www.gnu.org/licenses/>.
    1615!
    1716! Copyright 2009-2020 Carl von Ossietzky Universitaet Oldenburg
    1817! Copyright 1997-2020 Leibniz Universitaet Hannover
    19 !------------------------------------------------------------------------------!
     18!--------------------------------------------------------------------------------------------------!
     19!
    2020!
    2121! Current revisions:
     
    2626! -----------------
    2727! $Id$
     28! file re-formatted to follow the PALM coding standard
     29!
     30!
     31! 4495 2020-04-13 20:11:20Z raasch
    2832! restart data handling with MPI-IO added
    2933!
    3034! 4481 2020-03-31 18:55:54Z maronga
    3135! ASCII output cleanup
    32 ! 
     36!
    3337! 4465 2020-03-20 11:35:48Z maronga
    34 ! Removed old ASCII outputm, some syntax layout adjustments, added output for
    35 ! rotor and tower diameters. Added warning message in case of NetCDF 3 (no
    36 ! WTM output file will be produced).
    37 !
     38! Removed old ASCII outputm, some syntax layout adjustments, added output for rotor and tower
     39! diameters. Added warning message in case of NetCDF 3 (no WTM output file will be produced).
     40!
    3841! 4460 2020-03-12 16:47:30Z oliver.maas
    3942! allow for simulating up to 10 000 wind turbines
    40 ! 
     43!
    4144! 4459 2020-03-12 09:35:23Z oliver.maas
    4245! avoid division by zero in tip loss correction factor calculation
    43 ! 
     46!
    4447! 4457 2020-03-11 14:20:43Z raasch
    4548! use statement for exchange horiz added
    46 ! 
     49!
    4750! 4447 2020-03-06 11:05:30Z oliver.maas
    4851! renamed wind_turbine_parameters namelist variables
    49 ! 
     52!
    5053! 4438 2020-03-03 20:49:28Z suehring
    5154! Bugfix: shifted netcdf preprocessor directive to correct position
    52 ! 
     55!
    5356! 4436 2020-03-03 12:47:02Z oliver.maas
    5457! added optional netcdf data input for wtm array input parameters
    55 ! 
     58!
    5659! 4426 2020-02-27 10:02:19Z oliver.maas
    57 ! define time as unlimited dimension so that no maximum number
    58 ! of time steps has to be given for wtm_data_output
    59 ! 
     60! define time as unlimited dimension so that no maximum number of time steps has to be given for
     61! wtm_data_output
     62!
    6063! 4423 2020-02-25 07:17:11Z maronga
    6164! Switched to serial output as data is aggerated before anyway.
    62 ! 
     65!
    6366! 4420 2020-02-24 14:13:56Z maronga
    6467! Added output control for wind turbine model
    65 ! 
     68!
    6669! 4412 2020-02-19 14:53:13Z maronga
    6770! Bugfix: corrected character length in dimension_names
    68 ! 
     71!
    6972! 4411 2020-02-18 14:28:02Z maronga
    7073! Added output in NetCDF format using DOM (only netcdf4-parallel is supported).
    7174! Old ASCII output is still available at the moment.
    72 ! 
     75!
    7376! 4360 2020-01-07 11:25:50Z suehring
    74 ! Introduction of wall_flags_total_0, which currently sets bits based on static
    75 ! topography information used in wall_flags_static_0
    76 ! 
     77! Introduction of wall_flags_total_0, which currently sets bits based on static topography
     78! information used in wall_flags_static_0
     79!
    7780! 4343 2019-12-17 12:26:12Z oliver.maas
    78 ! replaced <= by < in line 1464 to ensure that ialpha will not be
    79 ! greater than dlen
    80 !
     81! replaced <= by < in line 1464 to ensure that ialpha will not be greater than dlen
     82!
    8183! 4329 2019-12-10 15:46:36Z motisi
    8284! Renamed wall_flags_0 to wall_flags_static_0
    83 ! 
     85!
    8486! 4326 2019-12-06 14:16:14Z oliver.maas
    85 ! changed format of turbine control output to allow for higher
    86 ! torque and power values
    87 !
     87! changed format of turbine control output to allow for higher torque and power values
     88!
    8889! 4182 2019-08-22 15:20:23Z scharf
    8990! Corrected "Former revisions" section
    90 ! 
     91!
    9192! 4144 2019-08-06 09:11:47Z raasch
    9293! relational operators .EQ., .NE., etc. replaced by ==, /=, etc.
    93 ! 
     94!
    9495! 4056 2019-06-27 13:53:16Z Giersch
    95 ! CASE DEFAULT action in wtm_actions needs to be CONTINUE. Otherwise an abort  
    96 ! will happen for location values that are not implemented as CASE statements
    97 ! but are already realized in the code (e.g. pt-tendency)
    98 ! 
     96! CASE DEFAULT action in wtm_actions needs to be CONTINUE. Otherwise an abort will happen for
     97! location values that are not implemented as CASE statements but are already realized in the code
     98! (e.g. pt-tendency)
     99!
    99100! 3885 2019-04-11 11:29:34Z kanani
    100 ! Changes related to global restructuring of location messages and introduction
    101 ! of additional debug messages
    102 ! 
     101! Changes related to global restructuring of location messages and introduction of additional debug
     102! messages
     103!
    103104! 3875 2019-04-08 17:35:12Z knoop
    104105! Addaped wtm_tendency to fit the module actions interface
    105 ! 
     106!
    106107! 3832 2019-03-28 13:16:58Z raasch
    107108! instrumented with openmp directives
    108 ! 
     109!
    109110! 3725 2019-02-07 10:11:02Z raasch
    110111! unused variables removed
    111 ! 
     112!
    112113! 3685 2019-01-21 01:02:11Z knoop
    113114! Some interface calls moved to module_interface + cleanup
    114 ! 
     115!
    115116! 3655 2019-01-07 16:51:22Z knoop
    116117! Replace degree symbol by 'degrees'
    117 ! 
     118!
    118119! 1914 2016-05-26 14:44:07Z witha
    119120! Initial revision
    120 ! 
     121!
    121122!
    122123! Description:
    123124! ------------
    124 !> This module calculates the effect of wind turbines on the flow fields. The
    125 !> initial version contains only the advanced actuator disk with rotation method
    126 !> (ADM-R).
     125!> This module calculates the effect of wind turbines on the flow fields. The initial version
     126!> contains only the advanced actuator disk with rotation method (ADM-R).
    127127!> The wind turbines include the tower effect, can be yawed and tilted.
    128 !> The wind turbine model includes controllers for rotational speed, pitch and
    129 !> yaw.
    130 !> Currently some specifications of the NREL 5 MW reference turbine
    131 !> are hardcoded whereas most input data comes from separate files (currently