Changeset 4591 for palm


Ignore:
Timestamp:
Jul 6, 2020 3:56:08 PM (4 months ago)
Author:
raasch
Message:

files re-formatted to follow the PALM coding standard

Location:
palm/trunk/SOURCE
Files:
8 edited

Legend:

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

    r4539 r4591  
    11!> @file restart_data_mpi_io_mod.f90
    2 !------------------------------------------------------------------------------!
     2!--------------------------------------------------------------------------------------------------!
    33! This file is part of the PALM model system.
    44!
     
    1515!
    1616! Copyright 1997-2020 Leibniz Universitaet Hannover
    17 ! -------------------------------------------------------------------------------------------------!
     17!--------------------------------------------------------------------------------------------------!
     18!
    1819!
    1920! Current revisions:
    2021! -----------------
    21 ! 
    22 ! 
     22!
     23!
    2324! Former revisions:
    2425! -----------------
    2526! $Id$
    26 ! checks added, if index limits in header are exceeded
    27 ! bugfix in rrd_mpi_io_int_2d
    28 !
     27! File re-formatted to follow the PALM coding standard
     28!
     29!
     30! 4539 2020-05-18 14:05:17Z raasch
     31! Checks added, if index limits in header are exceeded
     32! Bugfix in rrd_mpi_io_int_2d
     33!
    2934! 4536 2020-05-17 17:24:13Z raasch
    30 ! messages and debug output converted to PALM routines
    31 ! 
     35! Messages and debug output converted to PALM routines
     36!
    3237! 4534 2020-05-14 18:35:22Z raasch
    3338! I/O on reduced number of cores added (using shared memory MPI)
    34 ! 
     39!
    3540! 4500 2020-04-17 10:12:45Z suehring
    3641! Fix too long lines
    37 ! 
     42!
    3843! 4498 2020-04-15 14:26:31Z raasch
    39 ! bugfix for creation of filetypes, argument removed from rd_mpi_io_open
    40 ! 
     44! Bugfix for creation of filetypes, argument removed from rd_mpi_io_open
     45!
    4146! 4497 2020-04-15 10:20:51Z raasch
    42 ! last bugfix deactivated because of compile problems
    43 ! 
     47! Last bugfix deactivated because of compile problems
     48!
    4449! 4496 2020-04-15 08:37:26Z raasch
    45 ! problem with posix read arguments for surface data fixed
    46 ! 
     50! Problem with posix read arguments for surface data fixed
     51!
    4752! 4495 2020-04-13 20:11:20Z raasch
    4853! Initial version (K. Ketelsen), adjusted to PALM formatting standards (s. Raasch)
    4954!
    50 ! 
     55!
    5156!
    5257! Description:
     
    6469#else
    6570    USE posix_interface,                                                                           &
    66         ONLY:  posix_close, posix_lseek, posix_open, posix_read, posix_write
     71        ONLY:  posix_close,                                                                        &
     72               posix_lseek,                                                                        &
     73               posix_open,                                                                         &
     74               posix_read,                                                                         &
     75               posix_write
    6776#endif
    6877
     
    7079
    7180    USE control_parameters,                                                                        &
    72         ONLY:  debug_output, debug_string, include_total_domain_boundaries, message_string,        &
    73                restart_data_format_input, restart_data_format_output, restart_file_size
     81        ONLY:  debug_output,                                                                       &
     82               debug_string,                                                                       &
     83               include_total_domain_boundaries,                                                    &
     84               message_string,                                                                     &
     85               restart_data_format_input,                                                          &
     86               restart_data_format_output,                                                         &
     87               restart_file_size
    7488
    7589    USE exchange_horiz_mod,                                                                        &
    76         ONLY:  exchange_horiz, exchange_horiz_2d
     90        ONLY:  exchange_horiz,                                                                     &
     91               exchange_horiz_2d
    7792
    7893    USE indices,                                                                                   &
    79         ONLY:  nbgp, nnx, nny, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nz, nzb, nzt
     94        ONLY:  nbgp,                                                                               &
     95               nnx,                                                                                &
     96               nny,                                                                                &
     97               nx,                                                                                 &
     98               nxl,                                                                                &
     99               nxlg,                                                                               &
     100               nxr,                                                                                &
     101               nxrg,                                                                               &
     102               ny,                                                                                 &
     103               nyn,                                                                                &
     104               nyng,                                                                               &
     105               nys,                                                                                &
     106               nysg,                                                                               &
     107               nz,                                                                                 &
     108               nzb,                                                                                &
     109               nzt
    80110
    81111    USE kinds
    82112
    83113    USE pegrid,                                                                                    &
    84         ONLY:  comm1dx, comm1dy, comm2d, myid, myidx, myidy, npex, npey, numprocs, pdims
     114        ONLY:  comm1dx,                                                                            &
     115               comm1dy,                                                                            &
     116               comm2d,                                                                             &
     117               myid,                                                                               &
     118               myidx,                                                                              &
     119               myidy,                                                                              &
     120               npex,                                                                               &
     121               npey,                                                                               &
     122               numprocs,                                                                           &
     123               pdims
    85124
    86125    USE shared_memory_io_mod,                                                                      &
    87         ONLY:  local_boundaries, sm_class
    88 
    89 
    90     IMPLICIT NONE
    91 
    92     CHARACTER(LEN=128) :: io_file_name  !> internal variable to communicate filename between
    93                                         !> different subroutines
     126        ONLY:  local_boundaries,                                                                   &
     127               sm_class
     128
     129
     130    IMPLICIT NONE
     131
     132    CHARACTER(LEN=128) ::  io_file_name  !> internal variable to communicate filename between different subroutines
    94133
    95134#if defined( __parallel )
     
    99138#else
    100139    INTEGER(iwp), PARAMETER ::  rd_offset_kind = C_SIZE_T         !<
    101     INTEGER(iwp), PARAMETER ::  rd_status_size = 1       !< Not required in sequential mode
    102 #endif
    103 
    104     INTEGER(iwp)            ::  debug_level = 1 !< TODO: replace with standard debug output steering
    105 
    106     INTEGER(iwp)            ::  comm_io       !< Communicator for MPI-IO
    107     INTEGER(iwp)            ::  fh            !< MPI-IO file handle
    108 #if defined( __parallel )
    109     INTEGER(iwp)            ::  fhs = -1      !< MPI-IO file handle to open file with comm2d always
    110 #endif
    111     INTEGER(iwp)            ::  ft_surf = -1  !< MPI filetype surface data
    112 #if defined( __parallel )
    113     INTEGER(iwp)            ::  ft_2di_nb     !< MPI filetype 2D array INTEGER no outer boundary
    114     INTEGER(iwp)            ::  ft_2d         !< MPI filetype 2D array REAL with outer boundaries
    115     INTEGER(iwp)            ::  ft_3d         !< MPI filetype 3D array REAL with outer boundaries
    116     INTEGER(iwp)            ::  ft_3dsoil     !< MPI filetype for 3d-soil array
    117 #endif
    118     INTEGER(iwp)            ::  glo_start     !< global start index on this PE
    119 #if defined( __parallel )
    120     INTEGER(iwp)            ::  local_start   !<
    121 #endif
    122     INTEGER(iwp)            ::  nr_iope       !<
    123     INTEGER(iwp)            ::  nr_val        !< local number of values in x and y direction
    124 #if defined( __parallel )
    125     INTEGER(iwp)            ::  win_2di
    126     INTEGER(iwp)            ::  win_2dr
    127     INTEGER(iwp)            ::  win_3dr
    128     INTEGER(iwp)            ::  win_3ds
    129     INTEGER(iwp)            ::  win_surf = -1
    130 #endif
    131     INTEGER(iwp)            ::  total_number_of_surface_values    !< total number of values for one variable
     140    INTEGER(iwp), PARAMETER ::  rd_status_size = 1                !< Not required in sequential mode
     141#endif
     142
     143    INTEGER(iwp)            ::  debug_level = 1  !< TODO: replace with standard debug output steering
     144
     145    INTEGER(iwp)            ::  comm_io          !< Communicator for MPI-IO
     146    INTEGER(iwp)            ::  fh               !< MPI-IO file handle
     147#if defined( __parallel )
     148    INTEGER(iwp)            ::  fhs = -1         !< MPI-IO file handle to open file with comm2d always
     149#endif
     150    INTEGER(iwp)            ::  ft_surf = -1     !< MPI filetype surface data
     151#if defined( __parallel )
     152    INTEGER(iwp)            ::  ft_2di_nb        !< MPI filetype 2D array INTEGER no outer boundary
     153    INTEGER(iwp)            ::  ft_2d            !< MPI filetype 2D array REAL with outer boundaries
     154    INTEGER(iwp)            ::  ft_3d            !< MPI filetype 3D array REAL with outer boundaries
     155    INTEGER(iwp)            ::  ft_3dsoil        !< MPI filetype for 3d-soil array
     156#endif
     157    INTEGER(iwp)            ::  glo_start        !< global start index on this PE
     158#if defined( __parallel )
     159    INTEGER(iwp)            ::  local_start      !<
     160#endif
     161    INTEGER(iwp)            ::  nr_iope          !<
     162    INTEGER(iwp)            ::  nr_val           !< local number of values in x and y direction
     163#if defined( __parallel )
     164    INTEGER(iwp)            ::  win_2di          !<
     165    INTEGER(iwp)            ::  win_2dr          !<
     166    INTEGER(iwp)            ::  win_3dr          !<
     167    INTEGER(iwp)            ::  win_3ds          !<
     168    INTEGER(iwp)            ::  win_surf = -1    !<
     169#endif
     170    INTEGER(iwp)            ::  total_number_of_surface_values  !< total number of values for one variable
    132171
    133172    INTEGER(KIND=rd_offset_kind) ::  array_position   !<
     
    137176
    138177    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  m_end_index     !<
     178    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  m_global_start  !<
    139179    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  m_start_index   !<
    140     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  m_global_start  !<
    141 
    142     LOGICAL ::  all_pes_write                     !< all PEs have data to write
    143     LOGICAL ::  filetypes_created                 !<
    144     LOGICAL ::  io_on_limited_cores_per_node      !< switch to shared memory MPI-IO
    145     LOGICAL ::  rd_flag                           !< file is opened for read
    146     LOGICAL ::  wr_flag                           !< file is opened for write
     180
     181
     182    LOGICAL ::  all_pes_write                 !< all PEs have data to write
     183    LOGICAL ::  filetypes_created             !<
     184    LOGICAL ::  io_on_limited_cores_per_node  !< switch to shared memory MPI-IO
     185    LOGICAL ::  rd_flag                       !< file is opened for read
     186    LOGICAL ::  wr_flag                       !< file is opened for write
    147187
    148188#if defined( __parallel )
     
    160200!-- General Header (first 32 byte in restart file)
    161201    TYPE general_header
     202       INTEGER(iwp) :: endian         !< little endian (1) or big endian (2) internal format
     203       INTEGER(iwp) :: i_outer_bound  !< if 1, outer boundaries are stored in restart file
     204       INTEGER(iwp) :: nr_arrays      !< number of arrays in restart files
     205       INTEGER(iwp) :: nr_char        !< number of Text strings entries in header
    162206       INTEGER(iwp) :: nr_int         !< number of INTEGER entries in header
    163        INTEGER(iwp) :: nr_char        !< number of Text strings entries in header
    164207       INTEGER(iwp) :: nr_real        !< number of REAL entries in header
    165        INTEGER(iwp) :: nr_arrays      !< number of arrays in restart files
    166208       INTEGER(iwp) :: total_nx       !< total number of points in x-direction
    167209       INTEGER(iwp) :: total_ny       !< total number of points in y-direction
    168        INTEGER(iwp) :: i_outer_bound  !< if 1, outer boundaries are stored in restart file
    169        INTEGER(iwp) :: endian         !< little endian (1) or big endian (2) internal format
    170210    END TYPE general_header
    171211
    172     TYPE(general_header), TARGET ::  tgh
    173 
    174     TYPE(sm_class)               ::  sm_io
     212    TYPE(general_header), TARGET ::  tgh    !<
     213
     214    TYPE(sm_class)               ::  sm_io  !<
    175215
    176216!
     
    194234    CHARACTER(LEN=32), DIMENSION(max_nr_arrays) ::  array_names
    195235    INTEGER(KIND=rd_offset_kind), DIMENSION(max_nr_arrays) :: array_offset
    196 
    197236    SAVE
    198237
     
    267306    END INTERFACE wrd_mpi_io_surface
    268307
    269     PUBLIC  rd_mpi_io_check_array, rd_mpi_io_close, rd_mpi_io_open, rrd_mpi_io,                    &
    270             rrd_mpi_io_global_array, rrd_mpi_io_surface, rd_mpi_io_surface_filetypes, wrd_mpi_io,  &
    271             wrd_mpi_io_global_array, wrd_mpi_io_surface
     308    PUBLIC  rd_mpi_io_check_array,                                                                 &
     309            rd_mpi_io_close,                                                                       &
     310            rd_mpi_io_open,                                                                        &
     311            rrd_mpi_io,                                                                            &
     312            rrd_mpi_io_global_array,                                                               &
     313            rrd_mpi_io_surface,                                                                    &
     314            rd_mpi_io_surface_filetypes,                                                           &
     315            wrd_mpi_io,                                                                            &
     316            wrd_mpi_io_global_array,                                                               &
     317            wrd_mpi_io_surface
    272318
    273319
     
    284330    IMPLICIT NONE
    285331
    286     CHARACTER(LEN=*), INTENT(IN)  ::  action                           !<
    287     CHARACTER(LEN=*), INTENT(IN)  ::  file_name                        !<
    288 
    289     INTEGER(iwp)                  ::  i                                !<
    290     INTEGER(iwp)                  ::  gh_size                          !<
    291 
    292     INTEGER(KIND=rd_offset_kind)  ::  offset                           !<
    293 
    294 #if defined( __parallel )
    295     INTEGER, DIMENSION(rd_status_size) ::  status                      !<
    296 #endif
    297 
    298     LOGICAL, INTENT(IN), OPTIONAL ::  open_for_global_io_only          !<
    299     LOGICAL                       ::  set_filetype                     !<
     332    CHARACTER(LEN=*), INTENT(IN) ::  action     !<
     333    CHARACTER(LEN=*), INTENT(IN) ::  file_name  !<
     334
     335    INTEGER(iwp)                 ::  i          !<
     336    INTEGER(iwp)                 ::  gh_size    !<
     337
     338    INTEGER(KIND=rd_offset_kind) ::  offset     !<
     339
     340#if defined( __parallel )
     341    INTEGER, DIMENSION(rd_status_size) ::  status  !<
     342#endif
     343
     344    LOGICAL, INTENT(IN), OPTIONAL ::  open_for_global_io_only  !<
     345    LOGICAL                       ::  set_filetype             !<
    300346
    301347#if ! defined( __parallel )
    302     TYPE(C_PTR)                   ::  buf_ptr                          !<
     348    TYPE(C_PTR)                   ::  buf_ptr  !<
    303349#endif
    304350
     
    377423
    378424          IF ( ierr /= 0 )  THEN
    379              message_string = 'error opening restart file "' // TRIM( io_file_name ) //      &
     425             message_string = 'error opening restart file "' // TRIM( io_file_name ) //            &
    380426                              '" for reading with MPI-IO'
    381427             CALL message( 'rrd_mpi_io_open', 'PA0727', 3, 2, 0, 6, 0 )
     
    400446
    401447          IF ( ierr /= 0 )  THEN
    402              message_string = 'error opening restart file "' // TRIM( io_file_name ) //      &
     448             message_string = 'error opening restart file "' // TRIM( io_file_name ) //            &
    403449                              '" for writing with MPI-IO'
    404450             CALL message( 'rrd_mpi_io_open', 'PA0728', 3, 2, 0, 6, 0 )
     
    418464
    419465       IF ( debug_output )  THEN
    420           WRITE( debug_string, * )  'open restart file "' // TRIM( io_file_name ) //            &
     466          WRITE( debug_string, * )  'open restart file "' // TRIM( io_file_name ) //               &
    421467                                    '" for read in serial mode (posix)'
    422468          CALL debug_message( debug_string, 'start' )
     
    426472
    427473       IF ( debug_output )  THEN
    428           WRITE( debug_string, * )  'open restart file "' // TRIM( io_file_name ) //            &
     474          WRITE( debug_string, * )  'open restart file "' // TRIM( io_file_name ) //               &
    429475                                    '" for read in serial mode (posix)'
    430476          CALL debug_message( debug_string, 'end' )
     
    434480
    435481       IF ( debug_output )  THEN
    436           WRITE( debug_string, * )  'open restart file "' // TRIM( io_file_name ) //            &
     482          WRITE( debug_string, * )  'open restart file "' // TRIM( io_file_name ) //               &
    437483                                    '" for write in serial mode (posix)'
    438484          CALL debug_message( debug_string, 'start' )
     
    442488
    443489       IF ( debug_output )  THEN
    444           WRITE( debug_string, * )  'open restart file "' // TRIM( io_file_name ) //            &
     490          WRITE( debug_string, * )  'open restart file "' // TRIM( io_file_name ) //               &
    445491                                    '" for write in serial mode (posix)'
    446492          CALL debug_message( debug_string, 'end' )
     
    479525    header_int_index = header_int_index+3
    480526
    481     DO   i = 1, max_nr_arrays
     527    DO  i = 1, max_nr_arrays
    482528       array_offset(i) = 0
    483529       array_names(i)  = ' '
     
    489535       IF ( sm_io%iam_io_pe )  THEN
    490536!
    491 !--       File is open for read.
     537!--       File is open for reading
    492538#if defined( __parallel )
    493539!--       Set the default view
     
    513559
    514560!
    515 !--    File types depend on if boundaries of the total domain is included in data. This has been
     561!--    File types depend on boundaries of the total domain being included in data. This has been
    516562!--    checked with the previous statement.
    517563       IF ( set_filetype )  THEN
     
    535581          CALL MPI_FILE_SEEK( fh, header_position, MPI_SEEK_SET, ierr )
    536582          CALL MPI_FILE_READ( fh, text_lines, SIZE( text_lines ) * 128, MPI_CHAR, status, ierr )
    537           header_position = header_position+size(text_lines) * 128
     583          header_position = header_position + SIZE ( text_lines ) * 128
    538584!
    539585!--       REAL values
     
    640686    INTEGER(iwp) ::  i  !<
    641687
    642     LOGICAl      ::  found  !<
     688    LOGICAl ::  found  !<
    643689
    644690
     
    666712    IMPLICIT NONE
    667713
    668     CHARACTER(LEN=*), INTENT(IN)   :: name
    669 
    670     INTEGER(iwp)                   ::  i
    671     INTEGER(KIND=iwp), INTENT(OUT) ::  value
    672 
    673     LOGICAL                        ::  found
     714    CHARACTER(LEN=*), INTENT(IN) :: name  !<
     715
     716    INTEGER(iwp)                   ::  i      !<
     717    INTEGER(KIND=iwp), INTENT(OUT) ::  value  !<
     718
     719    LOGICAL ::  found  !<
    674720
    675721
     
    703749    IMPLICIT NONE
    704750
    705     CHARACTER(LEN=*), INTENT(IN)   ::  name
    706 
    707     INTEGER(iwp)                   ::  i
    708 
    709     LOGICAL                        ::  found
    710 
    711     REAL(KIND=wp), INTENT(OUT)     ::  value
     751    CHARACTER(LEN=*), INTENT(IN) ::  name   !<
     752
     753    INTEGER(iwp)                 ::  i      !<
     754
     755    LOGICAL                      ::  found  !<
     756
     757    REAL(KIND=wp), INTENT(OUT)   ::  value  !<
    712758
    713759
     
    741787    IMPLICIT NONE
    742788
    743     CHARACTER(LEN=*), INTENT(IN)       ::  name
    744 
    745 #if defined( __parallel )
    746     INTEGER, DIMENSION(rd_status_size) ::  status
    747 #endif
    748     INTEGER(iwp)                       ::  i
    749 
    750     LOGICAL                            ::  found
    751 
    752     REAL(wp), INTENT(INOUT), DIMENSION(nysg:nyng,nxlg:nxrg) ::  data
     789    CHARACTER(LEN=*), INTENT(IN)       ::  name    !<
     790
     791#if defined( __parallel )
     792    INTEGER, DIMENSION(rd_status_size) ::  status  !<
     793#endif
     794    INTEGER(iwp)                       ::  i       !<
     795
     796    LOGICAL                            ::  found   !<
     797
     798    REAL(wp), INTENT(INOUT), DIMENSION(nysg:nyng,nxlg:nxrg) ::  data  !<
    753799
    754800
     
    765811    IF ( found )  THEN
    766812#if defined( __parallel )
    767        CALL sm_io%sm_node_barrier()  ! has no effect if I/O on limited number of cores is inactive
     813       CALL sm_io%sm_node_barrier()  ! Has no effect if I/O on limited number of cores is inactive
    768814       IF ( sm_io%iam_io_pe )  THEN
    769           CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_REAL, ft_2d, 'native', MPI_INFO_NULL,   &
     815          CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_REAL, ft_2d, 'native', MPI_INFO_NULL,    &
    770816                                  ierr )
    771817          CALL MPI_FILE_READ_ALL( fh, array_2d, SIZE( array_2d ), MPI_REAL, status, ierr )
     
    809855    IMPLICIT NONE
    810856
    811     CHARACTER(LEN=*), INTENT(IN)        ::  name
    812 
    813     INTEGER(iwp)                        ::  i
    814     INTEGER(iwp)                        ::  j
    815 
    816 #if defined( __parallel )
    817     INTEGER, DIMENSION(rd_status_size)  ::  status
    818 #endif
    819 
    820     INTEGER(KIND=iwp), INTENT(INOUT), DIMENSION(:,:) ::  data
    821 
    822     LOGICAL                             ::  found
     857    CHARACTER(LEN=*), INTENT(IN)       ::  name    !<
     858
     859    INTEGER(iwp)                       ::  i       !<
     860    INTEGER(iwp)                       ::  j       !<
     861
     862#if defined( __parallel )
     863    INTEGER, DIMENSION(rd_status_size) ::  status  !<
     864#endif
     865
     866    INTEGER(KIND=iwp), INTENT(INOUT), DIMENSION(:,:) ::  data  !<
     867
     868    LOGICAL ::  found  !<
    823869
    824870
     
    835881    IF ( found )  THEN
    836882
    837        IF ( ( nxr - nxl + 1 + 2*nbgp ) == SIZE( data, 2 ) )  THEN
     883       IF ( ( nxr - nxl + 1 + 2 * nbgp ) == SIZE( data, 2 ) )  THEN
    838884!
    839885!--       Output array with Halos.
    840 !--       ATTENTION: INTEGER array with ghost boundaries are not implemented yet. This kind of array
    841 !--                  would be dimensioned in the caller subroutine like this:
     886!--       ATTENTION: INTEGER arrays with ghost boundaries are not implemented yet. This kind of
     887!--                  array would be dimensioned in the caller subroutine like this:
    842888!--                  INTEGER, DIMENSION(nysg:nyng,nxlg:nxrg)::  data
    843889          message_string = '2d-INTEGER array "' // TRIM( name ) // '" to be read from restart ' // &
     
    845891          CALL message( 'rrd_mpi_io_int_2d', 'PA0723', 3, 2, 0, 6, 0 )
    846892
    847        ELSEIF ( (nxr-nxl+1) == SIZE( data, 2 ) )  THEN
     893       ELSEIF ( (nxr - nxl + 1) == SIZE( data, 2 ) )  THEN
    848894!
    849895!--       INTEGER input array without Halos.
     
    852898
    853899#if defined( __parallel )
    854           CALL sm_io%sm_node_barrier() ! has no effect if I/O on limited number of cores is inactive
     900          CALL sm_io%sm_node_barrier() ! Has no effect if I/O on limited number of cores is inactive
    855901          IF ( sm_io%iam_io_pe )  THEN
    856902             CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_INTEGER, ft_2di_nb, 'native',         &
     
    898944    IMPLICIT NONE
    899945
    900     CHARACTER(LEN=*), INTENT(IN)       ::  name
    901 
    902     INTEGER(iwp)                       ::  i
    903 
    904 #if defined( __parallel )
    905     INTEGER, DIMENSION(rd_status_size) ::  status
    906 #endif
    907 
    908     LOGICAL                            ::  found
    909 
    910     REAL(wp), INTENT(INOUT), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  data
     946    CHARACTER(LEN=*), INTENT(IN)       ::  name    !<
     947
     948    INTEGER(iwp)                       ::  i       !<
     949
     950#if defined( __parallel )
     951    INTEGER, DIMENSION(rd_status_size) ::  status  !<
     952#endif
     953
     954    LOGICAL                            ::  found   !<
     955
     956    REAL(wp), INTENT(INOUT), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  data  !<
    911957
    912958
     
    923969    IF ( found )  THEN
    924970#if defined( __parallel )
    925        CALL sm_io%sm_node_barrier()  ! has no effect if I/O on limited number of cores is inactive
     971       CALL sm_io%sm_node_barrier()  ! Has no effect if I/O on limited number of cores is inactive
    926972       IF( sm_io%iam_io_pe )  THEN
    927973          CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_REAL, ft_3d, 'native', MPI_INFO_NULL,    &
     
    9681014    IMPLICIT NONE
    9691015
    970     CHARACTER(LEN=*), INTENT(IN)       ::  name
    971 
    972     INTEGER(iwp)                       ::  i
    973     INTEGER, INTENT(IN)                ::  nzb_soil
    974     INTEGER, INTENT(IN)                ::  nzt_soil
    975 
    976 #if defined( __parallel )
    977     INTEGER, DIMENSION(rd_status_size) ::  status
    978 #endif
    979 
    980     LOGICAL                            ::  found
    981 
    982     REAL(wp), INTENT(INOUT), DIMENSION(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) ::  data
     1016    CHARACTER(LEN=*), INTENT(IN)       ::  name      !<
     1017
     1018    INTEGER(iwp)                       ::  i         !<
     1019    INTEGER, INTENT(IN)                ::  nzb_soil  !<
     1020    INTEGER, INTENT(IN)                ::  nzt_soil  !<
     1021
     1022#if defined( __parallel )
     1023    INTEGER, DIMENSION(rd_status_size) ::  status    !<
     1024#endif
     1025
     1026    LOGICAL                            ::  found     !<
     1027
     1028    REAL(wp), INTENT(INOUT), DIMENSION(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) ::  data  !<
    9831029
    9841030
     
    9961042#if defined( __parallel )
    9971043       CALL rd_mpi_io_create_filetypes_3dsoil( nzb_soil, nzt_soil )
    998        CALL sm_io%sm_node_barrier()  ! has no effect if I/O on limited number of cores is inactive
     1044       CALL sm_io%sm_node_barrier()  ! Has no effect if I/O on limited number of cores is inactive
    9991045       IF ( sm_io%iam_io_pe )  THEN
    1000           CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_REAL, ft_3dsoil, 'native', MPI_INFO_NULL,&
    1001                                   ierr )
     1046          CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_REAL, ft_3dsoil, 'native',               &
     1047                                  MPI_INFO_NULL, ierr )
    10021048          CALL MPI_FILE_READ_ALL( fh, array_3d_soil, SIZE( array_3d_soil ), MPI_REAL, status, ierr )
    10031049          CALL MPI_TYPE_FREE( ft_3dsoil, ierr )
     
    10381084    IMPLICIT NONE
    10391085
    1040     CHARACTER(LEN=*), INTENT(IN)   ::  name
    1041     CHARACTER(LEN=*), INTENT(OUT)  ::  text
    1042     CHARACTER(LEN=128)             ::  line
    1043 
    1044     INTEGER(iwp)                   ::  i
    1045 
    1046     LOGICAL                        ::  found
     1086    CHARACTER(LEN=*), INTENT(IN)  ::  name   !<
     1087    CHARACTER(LEN=*), INTENT(OUT) ::  text   !<
     1088    CHARACTER(LEN=128)            ::  line   !<
     1089
     1090    INTEGER(iwp)                  ::  i      !<
     1091
     1092    LOGICAL                       ::  found  !<
    10471093
    10481094
     
    10771123    IMPLICIT NONE
    10781124
    1079     CHARACTER(LEN=*), INTENT(IN) ::  name
    1080 
    1081     INTEGER(iwp)                 ::  logical_as_integer
    1082 
    1083     LOGICAL, INTENT(OUT)         ::  value
     1125    CHARACTER(LEN=*), INTENT(IN) ::  name                !<
     1126
     1127    INTEGER(iwp)                 ::  logical_as_integer  !<
     1128
     1129    LOGICAL, INTENT(OUT)         ::  value               !<
    10841130
    10851131
     
    11001146    IMPLICIT NONE
    11011147
    1102     CHARACTER(LEN=*), INTENT(IN)  ::  name
    1103 
    1104     INTEGER(KIND=iwp), INTENT(IN) ::  value
     1148    CHARACTER(LEN=*), INTENT(IN)  ::  name   !<
     1149
     1150    INTEGER(KIND=iwp), INTENT(IN) ::  value  !<
    11051151
    11061152
     
    11161162
    11171163
    1118 
     1164!--------------------------------------------------------------------------------------------------!
     1165! Description:
     1166! ------------
     1167!> To do: Description missing!
     1168!--------------------------------------------------------------------------------------------------!
    11191169 SUBROUTINE wrd_mpi_io_real( name, value )
    11201170
    11211171    IMPLICIT NONE
    11221172
    1123     CHARACTER(LEN=*), INTENT(IN) ::  name
    1124 
    1125     REAL(wp), INTENT(IN)         ::  value
     1173    CHARACTER(LEN=*), INTENT(IN) ::  name   !<
     1174
     1175    REAL(wp), INTENT(IN)         ::  value  !<
    11261176
    11271177
     
    11471197    IMPLICIT NONE
    11481198
    1149     CHARACTER(LEN=*), INTENT(IN)       ::  name
    1150 
    1151     INTEGER(iwp)                       ::  i
    1152 
    1153 #if defined( __parallel )
    1154     INTEGER, DIMENSION(rd_status_size) ::  status
    1155 #endif
    1156 
    1157     REAL(wp), INTENT(IN), DIMENSION(nysg:nyng,nxlg:nxrg)    :: data
     1199    CHARACTER(LEN=*), INTENT(IN)       ::  name    !<
     1200
     1201    INTEGER(iwp)                       ::  i       !<
     1202
     1203#if defined( __parallel )
     1204    INTEGER, DIMENSION(rd_status_size) ::  status  !<
     1205#endif
     1206
     1207    REAL(wp), INTENT(IN), DIMENSION(nysg:nyng,nxlg:nxrg) ::  data  !<
    11581208
    11591209
     
    11681218    IF ( include_total_domain_boundaries )  THEN
    11691219!
    1170 !--    Prepare Output with outer boundaries
     1220!--    Prepare output with outer boundaries
    11711221       DO  i = lb%nxl, lb%nxr
    11721222          array_2d(i,lb%nys:lb%nyn) = data(lb%nys-nbgp:lb%nyn-nbgp,i-nbgp)
     
    11751225    ELSE
    11761226!
    1177 !--    Prepare Output without outer boundaries
     1227!--    Prepare output without outer boundaries
    11781228       DO  i = nxl,nxr
    11791229          array_2d(i,lb%nys:lb%nyn) = data(nys:nyn,i)
     
    11831233
    11841234#if defined( __parallel )
    1185     CALL sm_io%sm_node_barrier()  ! has no effect if I/O on limited number of cores is inactive
     1235    CALL sm_io%sm_node_barrier()  ! Has no effect if I/O on limited number of cores is inactive
    11861236    IF ( sm_io%iam_io_pe )  THEN
    11871237       CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_REAL, ft_2d, 'native', MPI_INFO_NULL, ierr )
     
    11941244#endif
    11951245!
    1196 !-- Type conversion required, otherwise rigth hand side brackets are calculated assuming 4 byte INT.
     1246!-- Type conversion required, otherwise right hand side brackets are calculated assuming 4 byte INT.
    11971247!-- Maybe a compiler problem.
    11981248    array_position = array_position + ( INT( lb%ny, KIND=rd_offset_kind ) + 1 ) *                  &
     
    12121262    IMPLICIT NONE
    12131263
    1214     CHARACTER(LEN=*), INTENT(IN)                  ::  name
    1215 
    1216     INTEGER(iwp)                                  ::  i
    1217     INTEGER(iwp)                                  ::  j
    1218 
    1219 #if defined( __parallel )
    1220     INTEGER, DIMENSION(rd_status_size)            ::  status
    1221 #endif
    1222     INTEGER(KIND=iwp), INTENT(IN), DIMENSION(:,:) ::  data
     1264    CHARACTER(LEN=*), INTENT(IN)                  ::  name    !<
     1265
     1266    INTEGER(iwp)                                  ::  i       !<
     1267    INTEGER(iwp)                                  ::  j       !<
     1268
     1269#if defined( __parallel )
     1270    INTEGER, DIMENSION(rd_status_size)            ::  status  !<
     1271#endif
     1272    INTEGER(KIND=iwp), INTENT(IN), DIMENSION(:,:) ::  data    !<
    12231273
    12241274
     
    12311281    header_array_index = header_array_index + 1
    12321282
    1233     IF ( ( nxr-nxl + 1 + 2 * nbgp ) == SIZE( data, 2 ) )  THEN
     1283    IF ( ( nxr - nxl + 1 + 2 * nbgp ) == SIZE( data, 2 ) )  THEN
    12341284!
    12351285!--    Integer arrays with ghost layers are not implemented yet. These kind of arrays would be
     
    12511301       ENDDO
    12521302#if defined( __parallel )
    1253        CALL sm_io%sm_node_barrier()  ! has no effect if I/O on limited number of cores is inactive
     1303       CALL sm_io%sm_node_barrier()  ! Has no effect if I/O on limited number of cores is inactive
    12541304       IF ( sm_io%iam_io_pe )  THEN
    12551305          CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_INTEGER, ft_2di_nb, 'native',            &
     
    12891339    IMPLICIT NONE
    12901340
    1291     CHARACTER(LEN=*), INTENT(IN)       ::  name
    1292 
    1293     INTEGER(iwp)                       ::  i
    1294 #if defined( __parallel )
    1295     INTEGER, DIMENSION(rd_status_size) ::  status
    1296 #endif
    1297     REAL(wp), INTENT(IN), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  ::  data
     1341    CHARACTER(LEN=*), INTENT(IN)       ::  name    !<
     1342
     1343    INTEGER(iwp)                       ::  i       !<
     1344#if defined( __parallel )
     1345    INTEGER, DIMENSION(rd_status_size) ::  status  !<
     1346#endif
     1347    REAL(wp), INTENT(IN), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  data !<
    12981348
    12991349
     
    13081358    IF ( include_total_domain_boundaries )  THEN
    13091359!
    1310 !--    Prepare output of 3d-REAL-array with ghost layers.
    1311 !--    In the virtual PE grid, the first dimension is PEs along x, and the second along y.
    1312 !--    For MPI-IO it is recommended to have the index order of the array in the same way, i.e.
    1313 !--    the first dimension should be along x and the second along y.
    1314 !--    For this reason, the original PALM data need to be swaped.
     1360!--    Prepare output of 3d-REAL-array with ghost layers. In the virtual PE grid, the first
     1361!--    dimension is PEs along x, and the second along y. For MPI-IO it is recommended to have the
     1362!--    index order of the array in the same way, i.e. the first dimension should be along x and the
     1363!--    second along y. For this reason, the original PALM data need to be swaped.
    13151364       DO  i = lb%nxl, lb%nxr
    13161365          array_3d(:,i,lb%nys:lb%nyn) = data(:,lb%nys-nbgp:lb%nyn-nbgp,i-nbgp)
     
    13261375    ENDIF
    13271376#if defined( __parallel )
    1328     CALL sm_io%sm_node_barrier()  ! has no effect if I/O on limited number of cores is inactive
     1377    CALL sm_io%sm_node_barrier()  ! Has no effect if I/O on limited number of cores is inactive
    13291378    IF ( sm_io%iam_io_pe )  THEN
    13301379       CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_REAL, ft_3d, 'native', MPI_INFO_NULL, ierr )
     
    13371386#endif
    13381387!
    1339 !-- Type conversion required, otherwise rigth hand side brackets are calculated assuming 4 byte INT.
     1388!-- Type conversion required, otherwise right hand side brackets are calculated assuming 4 byte INT.
    13401389!-- Maybe a compiler problem.
    1341     array_position = array_position + INT(    (nz+2), KIND=rd_offset_kind ) *                      &
    1342                                       INT( (lb%ny+1), KIND=rd_offset_kind ) *                      &
    1343                                       INT( (lb%nx+1), KIND=rd_offset_kind ) * wp
     1390    array_position = array_position + INT(    (nz+2), KIND = rd_offset_kind ) *                    &
     1391                                      INT( (lb%ny+1), KIND = rd_offset_kind ) *                    &
     1392                                      INT( (lb%nx+1), KIND = rd_offset_kind ) * wp
    13441393
    13451394 END SUBROUTINE wrd_mpi_io_real_3d
     
    13581407    IMPLICIT NONE
    13591408
    1360     CHARACTER(LEN=*), INTENT(IN)       ::  name
    1361 
    1362     INTEGER(iwp)                       ::  i
    1363     INTEGER, INTENT(IN)                ::  nzb_soil
    1364     INTEGER, INTENT(IN)                ::  nzt_soil
    1365 
    1366 #if defined( __parallel )
    1367     INTEGER, DIMENSION(rd_status_size) ::  status
    1368 #endif
    1369 
    1370     REAL(wp), INTENT(IN), DIMENSION(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg)  ::  data
     1409    CHARACTER(LEN=*), INTENT(IN)       ::  name      !<
     1410
     1411    INTEGER(iwp)                       ::  i         !<
     1412    INTEGER, INTENT(IN)                ::  nzb_soil  !<
     1413    INTEGER, INTENT(IN)                ::  nzt_soil  !<
     1414
     1415#if defined( __parallel )
     1416    INTEGER, DIMENSION(rd_status_size) ::  status  !<
     1417#endif
     1418
     1419    REAL(wp), INTENT(IN), DIMENSION(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) ::  data  !<
    13711420
    13721421
     
    13851434    IF ( include_total_domain_boundaries)  THEN
    13861435!
    1387 !--    Prepare output of 3d-REAL-array with ghost layers.
    1388 !--    In the virtual PE grid, the first dimension is PEs along x, and the second along y.
    1389 !--    For MPI-IO it is recommended to have the index order of the array in the same way, i.e.
    1390 !--    the first dimension should be along x and the second along y.
    1391 !--    For this reason, the original PALM data need to be swaped.
     1436!--    Prepare output of 3d-REAL-array with ghost layers. In the virtual PE grid, the first
     1437!--    dimension is PEs along x, and the second along y. For MPI-IO it is recommended to have the
     1438!--    index order of the array in the same way, i.e. the first dimension should be along x and the
     1439!--    second along y. For this reason, the original PALM data need to be swaped.
    13921440       DO  i = lb%nxl, lb%nxr
    13931441          array_3d_soil(:,i,lb%nys:lb%nyn) = data(:,lb%nys-nbgp:lb%nyn-nbgp,i-nbgp)
     
    14031451    ENDIF
    14041452#if defined( __parallel )
    1405     CALL sm_io%sm_node_barrier()  ! has no effect if I/O on limited number of cores is inactive
     1453    CALL sm_io%sm_node_barrier()  ! Has no effect if I/O on limited number of cores is inactive
    14061454    IF ( sm_io%iam_io_pe )  THEN
    14071455       CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_REAL, ft_3dsoil, 'native', MPI_INFO_NULL,   &
     
    14151463#endif
    14161464!
    1417 !-- Type conversion required, otherwise rigth hand side brackets are calculated assuming 4 byte INT.
     1465!-- Type conversion required, otherwise right hand side brackets are calculated assuming 4 byte INT.
    14181466!-- Maybe a compiler problem.
    1419     array_position = array_position + INT( (nzt_soil-nzb_soil+1), KIND=rd_offset_kind ) *          &
    1420                                       INT( (lb%ny+1),             KIND=rd_offset_kind ) *          &
    1421                                       INT( (lb%nx+1),             KIND=rd_offset_kind ) * wp
     1467    array_position = array_position + INT( (nzt_soil-nzb_soil+1), KIND = rd_offset_kind ) *        &
     1468                                      INT( (lb%ny+1),             KIND = rd_offset_kind ) *        &
     1469                                      INT( (lb%nx+1),             KIND = rd_offset_kind ) * wp
    14221470
    14231471 END SUBROUTINE wrd_mpi_io_real_3d_soil
     
    14341482    IMPLICIT NONE
    14351483
    1436     CHARACTER(LEN=128)           ::  lo_line
    1437     CHARACTER(LEN=*), INTENT(IN) ::  name
    1438     CHARACTER(LEN=*), INTENT(IN) ::  text
     1484    CHARACTER(LEN=128)           ::  lo_line  !<
     1485    CHARACTER(LEN=*), INTENT(IN) ::  name     !<
     1486    CHARACTER(LEN=*), INTENT(IN) ::  text     !<
    14391487
    14401488
     
    14611509    IMPLICIT NONE
    14621510
    1463     CHARACTER(LEN=*), INTENT(IN) ::  name
    1464 
    1465     INTEGER(iwp)                 ::  logical_as_integer
    1466 
    1467     LOGICAL, INTENT(IN)          ::  value
     1511    CHARACTER(LEN=*), INTENT(IN) ::  name                !<
     1512
     1513    INTEGER(iwp)                 ::  logical_as_integer  !<
     1514
     1515    LOGICAL, INTENT(IN)          ::  value               !<
    14681516
    14691517
     
    14901538    IMPLICIT NONE
    14911539
    1492     CHARACTER(LEN=*), INTENT(IN)       ::  name
    1493 
    1494     INTEGER(iwp)                       ::  i
    1495     INTEGER(KIND=rd_offset_kind)       ::  offset
    1496 
    1497 #if defined( __parallel )
    1498     INTEGER, DIMENSION(rd_status_size) ::  status
    1499 #endif
    1500 
    1501     LOGICAL                            ::  found
    1502 
    1503     REAL(KIND=wp), INTENT(INOUT), DIMENSION(:) ::  data
     1540    CHARACTER(LEN=*), INTENT(IN)       ::  name    !<
     1541
     1542    INTEGER(iwp)                       ::  i       !<
     1543    INTEGER(KIND=rd_offset_kind)       ::  offset  !<
     1544
     1545#if defined( __parallel )
     1546    INTEGER, DIMENSION(rd_status_size) ::  status  !<
     1547#endif
     1548
     1549    LOGICAL                            ::  found   !<
     1550
     1551    REAL(KIND=wp), INTENT(INOUT), DIMENSION(:) ::  data  !<
    15041552
    15051553
     
    15251573       ENDIF
    15261574       IF ( sm_io%is_sm_active() )  THEN
    1527           CALL MPI_BCAST( data, SIZE(data), MPI_REAL, 0, sm_io%comm_shared, ierr )
     1575          CALL MPI_BCAST( data, SIZE( data ), MPI_REAL, 0, sm_io%comm_shared, ierr )
    15281576       ENDIF
    15291577#else
     
    15541602    IMPLICIT NONE
    15551603
    1556     CHARACTER(LEN=*), INTENT(IN)                      ::  name
    1557 
    1558     INTEGER, DIMENSION(1)                             ::  bufshape
    1559 
    1560     REAL(KIND=wp), INTENT(IN), DIMENSION(:,:), TARGET ::  data
    1561     REAL(KIND=wp), POINTER, DIMENSION(:)              ::  buf
    1562 
    1563     TYPE(C_PTR)                                       ::  c_data
     1604    CHARACTER(LEN=*), INTENT(IN)                      ::  name      !<
     1605
     1606    INTEGER, DIMENSION(1)                             ::  bufshape  !<
     1607
     1608    REAL(KIND=wp), INTENT(IN), DIMENSION(:,:), TARGET ::  data      !<
     1609    REAL(KIND=wp), POINTER, DIMENSION(:)              ::  buf       !<
     1610
     1611    TYPE(C_PTR)                                       ::  c_data    !<
    15641612
    15651613
     
    15841632    IMPLICIT NONE
    15851633
    1586     CHARACTER(LEN=*), INTENT(IN)                        ::  name
    1587 
    1588     INTEGER, DIMENSION(1)                               ::  bufshape
    1589 
    1590     REAL(KIND=wp), INTENT(IN), DIMENSION(:,:,:), TARGET ::  data
    1591     REAL(KIND=wp), POINTER, DIMENSION(:)                ::  buf
    1592 
    1593     TYPE(C_PTR)                                         ::  c_data
     1634    CHARACTER(LEN=*), INTENT(IN)                        ::  name      !<
     1635
     1636    INTEGER, DIMENSION(1)                               ::  bufshape  !<
     1637
     1638    REAL(KIND=wp), INTENT(IN), DIMENSION(:,:,:), TARGET ::  data      !<
     1639    REAL(KIND=wp), POINTER, DIMENSION(:)                ::  buf       !<
     1640
     1641    TYPE(C_PTR)                                         ::  c_data    !<
    15941642
    15951643
     
    16141662    IMPLICIT NONE
    16151663
    1616     CHARACTER(LEN=*), INTENT(IN)                          ::  name
    1617 
    1618     INTEGER, DIMENSION(1)                                 ::  bufshape
    1619 
    1620     REAL(KIND=wp), INTENT(IN), DIMENSION(:,:,:,:), TARGET ::  data
    1621     REAL(KIND=wp), POINTER, DIMENSION(:)                  ::  buf
    1622 
    1623     TYPE(C_PTR)                                           ::  c_data
     1664    CHARACTER(LEN=*), INTENT(IN)                          ::  name      !<
     1665
     1666    INTEGER, DIMENSION(1)                                 ::  bufshape  !<
     1667
     1668    REAL(KIND=wp), INTENT(IN), DIMENSION(:,:,:,:), TARGET ::  data      !<
     1669    REAL(KIND=wp), POINTER, DIMENSION(:)                  ::  buf       !<
     1670
     1671    TYPE(C_PTR)                                           ::  c_data    !<
    16241672
    16251673
     
    16441692    IMPLICIT NONE
    16451693
    1646     CHARACTER(LEN=*), INTENT(IN)                   ::  name
    1647 
    1648     INTEGER(iwp)                                   ::  i
    1649     INTEGER(KIND=rd_offset_kind)                   ::  offset
    1650 
    1651 #if defined( __parallel )
    1652     INTEGER, DIMENSION(rd_status_size)             ::  status
    1653 #endif
    1654     INTEGER(KIND=iwp), INTENT(INOUT), DIMENSION(:) ::  data
    1655 
    1656     LOGICAL                                        ::  found
     1694    CHARACTER(LEN=*), INTENT(IN)                   ::  name    !<
     1695
     1696    INTEGER(iwp)                                   ::  i       !<
     1697    INTEGER(KIND=rd_offset_kind)                   ::  offset  !<
     1698
     1699#if defined( __parallel )
     1700    INTEGER, DIMENSION(rd_status_size)             ::  status  !<
     1701#endif
     1702    INTEGER(KIND=iwp), INTENT(INOUT), DIMENSION(:) ::  data    !<
     1703
     1704    LOGICAL                                        ::  found   !<
    16571705
    16581706
     
    16611709
    16621710    DO  i = 1, tgh%nr_arrays
    1663        IF ( TRIM(array_names(i)) == TRIM( name ) )  THEN
     1711       IF ( TRIM( array_names(i) ) == TRIM( name ) )  THEN
    16641712          array_position = array_offset(i)
    16651713          found = .TRUE.
     
    16781726       ENDIF
    16791727       IF ( sm_io%is_sm_active() )  THEN
    1680           CALL MPI_BCAST( data, SIZE(data), MPI_INTEGER, 0, sm_io%comm_shared, ierr )
     1728          CALL MPI_BCAST( data, SIZE( data ), MPI_INTEGER, 0, sm_io%comm_shared, ierr )
    16811729       ENDIF
    16821730#else
     
    17061754    IMPLICIT NONE
    17071755
    1708     CHARACTER(LEN=*), INTENT(IN)            ::  name
    1709 
    1710     INTEGER(KIND=rd_offset_kind)            ::  offset
    1711 
    1712 #if defined( __parallel )
    1713     INTEGER, DIMENSION(rd_status_size)      ::  status
    1714 #endif
    1715 
    1716     REAL(KIND=wp), INTENT(IN), DIMENSION(:) ::  data
     1756    CHARACTER(LEN=*), INTENT(IN)            ::  name    !<
     1757
     1758    INTEGER(KIND=rd_offset_kind)            ::  offset  !<
     1759
     1760#if defined( __parallel )
     1761    INTEGER, DIMENSION(rd_status_size)      ::  status  !<
     1762#endif
     1763
     1764    REAL(KIND=wp), INTENT(IN), DIMENSION(:) ::  data    !<
    17171765
    17181766
     
    17371785    IF ( myid == 0 )  THEN
    17381786       CALL MPI_FILE_SEEK( fh, array_position, MPI_SEEK_SET, ierr )
    1739        CALL MPI_FILE_WRITE( fh, data, SIZE( data), MPI_REAL, status, ierr )
     1787       CALL MPI_FILE_WRITE( fh, data, SIZE( data ), MPI_REAL, status, ierr )
    17401788    ENDIF
    17411789#else
     
    17591807    IMPLICIT NONE
    17601808
    1761     CHARACTER(LEN=*), INTENT(IN)                      ::  name
    1762 
    1763     INTEGER, DIMENSION(1)                             ::  bufshape
    1764 
    1765     REAL(KIND=wp), POINTER, DIMENSION(:)              ::  buf
    1766     REAL(KIND=wp), INTENT(IN), DIMENSION(:,:), TARGET ::  data
    1767 
    1768     TYPE(C_PTR)                                       ::  c_data
     1809    CHARACTER(LEN=*), INTENT(IN)                      ::  name      !<
     1810
     1811    INTEGER, DIMENSION(1)                             ::  bufshape  !<
     1812
     1813    REAL(KIND=wp), POINTER, DIMENSION(:)              ::  buf       !<
     1814    REAL(KIND=wp), INTENT(IN), DIMENSION(:,:), TARGET ::  data      !<
     1815
     1816    TYPE(C_PTR)                                       ::  c_data    !<
    17691817
    17701818
     
    17891837    IMPLICIT NONE
    17901838
    1791     CHARACTER(LEN=*), INTENT(IN)                        ::  name
    1792 
    1793     INTEGER, DIMENSION(1)                               ::  bufshape
    1794 
    1795     REAL(KIND=wp), POINTER, DIMENSION(:)                ::  buf
    1796     REAL(KIND=wp), INTENT(IN), DIMENSION(:,:,:), TARGET ::  data
    1797 
    1798     TYPE(C_PTR)                                         ::  c_data
     1839    CHARACTER(LEN=*), INTENT(IN)                        ::  name      !<
     1840
     1841    INTEGER, DIMENSION(1)                               ::  bufshape  !<
     1842
     1843    REAL(KIND=wp), POINTER, DIMENSION(:)                ::  buf       !<
     1844    REAL(KIND=wp), INTENT(IN), DIMENSION(:,:,:), TARGET ::  data      !<
     1845
     1846    TYPE(C_PTR)                                         ::  c_data    !<
    17991847
    18001848
     
    18191867    IMPLICIT NONE
    18201868
    1821     CHARACTER(LEN=*), INTENT(IN)                          ::  name
    1822 
    1823     INTEGER, DIMENSION(1)                                 ::  bufshape
    1824 
    1825     REAL(KIND=wp), POINTER, DIMENSION(:)                  ::  buf
    1826     REAL(KIND=wp), INTENT(IN), DIMENSION(:,:,:,:), TARGET ::  data
    1827 
    1828     TYPE(C_PTR)                                           ::  c_data
     1869    CHARACTER(LEN=*), INTENT(IN)                          ::  name      !<
     1870
     1871    INTEGER, DIMENSION(1)                                 ::  bufshape  !<
     1872
     1873    REAL(KIND=wp), POINTER, DIMENSION(:)                  ::  buf       !<
     1874    REAL(KIND=wp), INTENT(IN), DIMENSION(:,:,:,:), TARGET ::  data      !<
     1875
     1876    TYPE(C_PTR)                                           ::  c_data    !<
    18291877
    18301878
     
    18491897    IMPLICIT NONE
    18501898
    1851     CHARACTER(LEN=*), INTENT(IN)                ::  name
    1852 
    1853     INTEGER(KIND=rd_offset_kind)                ::  offset
    1854 
    1855     INTEGER(KIND=iwp), INTENT(IN), DIMENSION(:) ::  data
    1856 #if defined( __parallel )
    1857     INTEGER, DIMENSION(rd_status_size)          ::  status
     1899    CHARACTER(LEN=*), INTENT(IN)                ::  name    !<
     1900
     1901    INTEGER(KIND=rd_offset_kind)                ::  offset  !<
     1902
     1903    INTEGER(KIND=iwp), INTENT(IN), DIMENSION(:) ::  data    !<
     1904#if defined( __parallel )
     1905    INTEGER, DIMENSION(rd_status_size)          ::  status  !<
    18581906#endif
    18591907
     
    18751923!
    18761924!-- Only PE 0 writes replicated data
    1877     IF ( myid == 0 )  THEN                        !
     1925    IF ( myid == 0 )  THEN
    18781926       CALL MPI_FILE_SEEK( fh, array_position, MPI_SEEK_SET, ierr )
    18791927       CALL MPI_FILE_WRITE( fh, data, SIZE( data), MPI_INTEGER, status, ierr )
     
    18981946    IMPLICIT NONE
    18991947
    1900     CHARACTER(LEN=*), INTENT(IN) ::  name
    1901 
    1902     INTEGER(KIND=rd_offset_kind) ::  disp          !< displacement of actual indices
    1903     INTEGER(KIND=rd_offset_kind) ::  disp_f        !< displacement in file
    1904     INTEGER(KIND=rd_offset_kind) ::  disp_n        !< displacement of next column
    1905     INTEGER(iwp), OPTIONAL       ::  first_index
    1906 
    1907     INTEGER(iwp)                 ::  i
    1908     INTEGER(iwp)                 ::  i_f
    1909     INTEGER(iwp)                 ::  j
    1910     INTEGER(iwp)                 ::  j_f
    1911     INTEGER(iwp)                 ::  lo_first_index
    1912     INTEGER(iwp)                 ::  nr_bytes
    1913     INTEGER(iwp)                 ::  nr_bytes_f
    1914     INTEGER(iwp)                 ::  nr_words
    1915 #if defined( __parallel )
    1916     INTEGER, DIMENSION(rd_status_size)  ::  status
    1917 #endif
    1918 
    1919     LOGICAL                             ::  found
    1920 
    1921     REAL(wp), INTENT(OUT), DIMENSION(:) ::  data
     1948    CHARACTER(LEN=*), INTENT(IN) ::  name            !<
     1949
     1950    INTEGER(KIND=rd_offset_kind) ::  disp            !< displacement of actual indices
     1951    INTEGER(KIND=rd_offset_kind) ::  disp_f          !< displacement in file
     1952    INTEGER(KIND=rd_offset_kind) ::  disp_n          !< displacement of next column
     1953    INTEGER(iwp), OPTIONAL       ::  first_index     !<
     1954
     1955    INTEGER(iwp)                 ::  i               !<
     1956    INTEGER(iwp)                 ::  i_f             !<
     1957    INTEGER(iwp)                 ::  j               !<
     1958    INTEGER(iwp)                 ::  j_f             !<
     1959    INTEGER(iwp)                 ::  lo_first_index  !<
     1960    INTEGER(iwp)                 ::  nr_bytes        !<
     1961    INTEGER(iwp)                 ::  nr_bytes_f      !<
     1962    INTEGER(iwp)                 ::  nr_words        !<
     1963#if defined( __parallel )
     1964    INTEGER, DIMENSION(rd_status_size)  ::  status   !<
     1965#endif
     1966
     1967    LOGICAL                             ::  found    !<
     1968
     1969    REAL(wp), INTENT(OUT), DIMENSION(:) ::  data     !<
    19221970
    19231971
     
    19251973    lo_first_index = 1
    19261974
    1927     IF ( MAXVAL( m_global_start ) == -1 )   RETURN   ! nothing to do on this PE
     1975    IF ( MAXVAL( m_global_start ) == -1 )   RETURN   ! Nothing to do on this PE
    19281976
    19291977    IF ( PRESENT( first_index ) )   THEN
     
    19341982        IF ( TRIM( array_names(i) ) == TRIM( name ) )  THEN
    19351983           array_position = array_offset(i) + ( lo_first_index - 1 ) *                             &
    1936                                               total_number_of_surface_values * wp
     1984                            total_number_of_surface_values * wp
    19371985           found = .TRUE.
    19381986           EXIT
     
    19532001                nr_bytes = nr_words * wp
    19542002             ENDIF
    1955              IF ( disp >= 0  .AND.  disp_f == -1 )  THEN   ! first Entry
     2003             IF ( disp >= 0  .AND.  disp_f == -1 )  THEN   ! First entry
    19562004                disp_f     = disp
    19572005                nr_bytes_f = 0
     
    19592007                j_f = j
    19602008             ENDIF
    1961              IF ( j == nyn  .AND.  i == nxr )  THEN        ! last Entry
     2009             IF ( j == nyn  .AND.  i == nxr )  THEN        ! Last entry
    19622010                disp_n = -1
    19632011                IF (  nr_bytes > 0 )  THEN
    19642012                   nr_bytes_f = nr_bytes_f+nr_bytes
    19652013                ENDIF
    1966              ELSEIF ( j == nyn )  THEN                     ! next x
     2014             ELSEIF ( j == nyn )  THEN                     ! Next x
    19672015                IF ( m_global_start(nys,i+1) > 0  .AND.  disp > 0 )  THEN
    19682016                   disp_n = array_position + ( m_global_start(nys,i+1) - 1 ) * wp
     
    19792027
    19802028
    1981              IF ( disp + nr_bytes == disp_n )  THEN        ! contiguous block
     2029             IF ( disp + nr_bytes == disp_n )  THEN        ! Contiguous block
    19822030                nr_bytes_f = nr_bytes_f + nr_bytes
    1983              ELSE                                          ! read
     2031             ELSE                                          ! Read
    19842032#if defined( __parallel )
    19852033                CALL MPI_FILE_SEEK( fhs, disp_f, MPI_SEEK_SET, ierr )
    19862034                nr_words = nr_bytes_f / wp
    1987                 CALL MPI_FILE_READ( fhs, data(m_start_index(j_f,i_f)), nr_words, MPI_REAL, status, ierr )
     2035                CALL MPI_FILE_READ( fhs, data(m_start_index(j_f,i_f)), nr_words, MPI_REAL, status, &
     2036                                    ierr )
    19882037#else
    19892038                CALL posix_lseek( fh, disp_f )
     
    20272076    IMPLICIT NONE
    20282077
    2029     CHARACTER(LEN=*), INTENT(IN)          ::  name
    2030 
    2031     INTEGER(iwp)                          ::  i
    2032 
    2033     REAL(wp), INTENT(OUT), DIMENSION(:,:) ::  data
    2034     REAL(wp), DIMENSION(SIZE( data,2))    ::  tmp
    2035 
    2036 
    2037     DO  i = 1, SIZE( data,1)
     2078    CHARACTER(LEN=*), INTENT(IN)          ::  name  !<
     2079
     2080    INTEGER(iwp)                          ::  i     !<
     2081
     2082    REAL(wp), INTENT(OUT), DIMENSION(:,:) ::  data  !<
     2083    REAL(wp), DIMENSION(SIZE( data,2))    ::  tmp   !<
     2084
     2085
     2086    DO  i = 1, SIZE( data, 1 )
    20382087       CALL rrd_mpi_io_surface( name, tmp, first_index = i )
    20392088       data(i,:) = tmp
     
    20422091!
    20432092!-- Comment from Klaus Ketelsen (September 2018)
    2044 !-- The intention of the following loop was let the compiler do the copying on return.
    2045 !-- In my understanding is it standard conform to pass the second dimension to a 1d-
    2046 !-- array inside a subroutine and the compiler is responsible to generate code for
    2047 !-- copying. Acually this works fine for INTENT(IN) variables (wrd_mpi_io_surface_2d).
    2048 !-- For INTENT(OUT) like in this case the code works on pgi compiler. But both, the Intel 16
    2049 !-- and the Cray compiler show wrong answers using this loop.
    2050 !-- That is the reason why the above auxiliary array tmp was introduced.
     2093!-- The intention of the following loop was to let the compiler do the copying on return.
     2094!-- In my understanding it is standard conform to pass the second dimension to a 1d-array inside a
     2095!-- subroutine and the compiler is responsible to generate code for copying. Acually this works fine
     2096!-- for INTENT(IN) variables (wrd_mpi_io_surface_2d). For INTENT(OUT) like in this case the code
     2097!-- works on pgi compiler. But both, the Intel 16 and the Cray compiler show wrong answers using
     2098!-- this loop. That is the reason why the above auxiliary array tmp was introduced.
    20512099!    DO  i = 1, SIZE(  data,1)
    20522100!       CALL rrd_mpi_io_surface( name, data(i,:), first_index = i )
     
    20662114    IMPLICIT NONE
    20672115
    2068     CHARACTER(LEN=*), INTENT(IN)       ::  name
    2069 
    2070 #if defined( __parallel )
    2071     INTEGER(KIND=rd_offset_kind)       ::  disp
    2072 #endif
    2073     INTEGER(iwp), OPTIONAL             ::  first_index
    2074 #if defined( __parallel )
    2075     INTEGER(iwp)                       ::  i
    2076 #endif
    2077     INTEGER(iwp)                       ::  lo_first_index
    2078     INTEGER(KIND=rd_offset_kind)       ::  offset
    2079 
    2080 #if defined( __parallel )
    2081     INTEGER, DIMENSION(rd_status_size) ::  status
    2082 #endif
    2083 
    2084     REAL(wp), INTENT(IN), DIMENSION(:), TARGET ::  data
     2116    CHARACTER(LEN=*), INTENT(IN)       ::  name            !<
     2117
     2118#if defined( __parallel )
     2119    INTEGER(KIND=rd_offset_kind)       ::  disp            !<
     2120#endif
     2121    INTEGER(iwp), OPTIONAL             ::  first_index     !<
     2122#if defined( __parallel )
     2123    INTEGER(iwp)                       ::  i               !<
     2124#endif
     2125    INTEGER(iwp)                       ::  lo_first_index  !<
     2126    INTEGER(KIND=rd_offset_kind)       ::  offset          !<
     2127
     2128#if defined( __parallel )
     2129    INTEGER, DIMENSION(rd_status_size) ::  status          !<
     2130#endif
     2131
     2132    REAL(wp), INTENT(IN), DIMENSION(:), TARGET ::  data    !<
    20852133
    20862134
     
    20882136    lo_first_index = 1
    20892137
    2090     IF ( PRESENT(first_index) )  THEN
     2138    IF ( PRESENT( first_index ) )  THEN
    20912139       lo_first_index = first_index
    20922140    ENDIF
     
    21162164    ENDIF
    21172165
    2118     CALL sm_io%sm_node_barrier()  ! has no effect if I/O on limited number of cores is inactive
     2166    CALL sm_io%sm_node_barrier()  ! Has no effect if I/O on limited number of cores is inactive
    21192167    IF ( sm_io%iam_io_pe )  THEN
    21202168       IF ( all_pes_write )  THEN
     
    21542202! ------------
    21552203!> Read 2d-REAL surface data array with MPI-IO.
    2156 !> These consist of multiple 1d-REAL surface data arrays.
     2204!> This consist of multiple 1d-REAL surface data arrays.
    21572205!--------------------------------------------------------------------------------------------------!
    21582206 SUBROUTINE wrd_mpi_io_surface_2d( name, data )
     
    21602208    IMPLICIT NONE
    21612209
    2162     CHARACTER(LEN=*), INTENT(IN)         ::  name
    2163 
    2164     INTEGER(iwp)                         ::  i
    2165 
    2166     REAL(wp), INTENT(IN), DIMENSION(:,:) ::  data
    2167 
    2168 
    2169     DO  i = 1, SIZE( data,1)
     2210    CHARACTER(LEN=*), INTENT(IN)         ::  name  !<
     2211
     2212    INTEGER(iwp)                         ::  i     !<
     2213
     2214    REAL(wp), INTENT(IN), DIMENSION(:,:) ::  data  !<
     2215
     2216
     2217    DO  i = 1, SIZE( data, 1 )
    21702218       CALL wrd_mpi_io_surface( name, data(i,:), first_index = i )
    21712219    ENDDO
     
    21842232    IMPLICIT NONE
    21852233
    2186     INTEGER(iwp)                       ::  gh_size
    2187     INTEGER(KIND=rd_offset_kind)       ::  offset
    2188 #if defined( __parallel )
    2189     INTEGER, DIMENSION(rd_status_size) ::  status
     2234    INTEGER(iwp)                       ::  gh_size  !<
     2235    INTEGER(KIND=rd_offset_kind)       ::  offset   !<
     2236#if defined( __parallel )
     2237    INTEGER, DIMENSION(rd_status_size) ::  status   !<
    21902238#endif
    21912239
    21922240#if ! defined( __parallel )
    2193     TYPE(C_PTR)                        ::  buf_ptr
     2241    TYPE(C_PTR)                        ::  buf_ptr  !<
    21942242#endif
    21952243
     
    22052253       tgh%total_nx  = lb%nx + 1
    22062254       tgh%total_ny  = lb%ny + 1
    2207        IF ( include_total_domain_boundaries )  THEN   ! not sure, if LOGICAL interpretation is the same for all compilers,
    2208           tgh%i_outer_bound = 1        ! therefore store as INTEGER in general header
     2255       IF ( include_total_domain_boundaries )  THEN   ! Not sure, if LOGICAL interpretation is the same for all compilers,
     2256          tgh%i_outer_bound = 1                       ! therefore store as INTEGER in general header
    22092257       ELSE
    22102258          tgh%i_outer_bound = 0
     
    22312279!--       INTEGER values
    22322280          CALL MPI_FILE_SEEK( fh, header_position, MPI_SEEK_SET, ierr )
    2233           CALL MPI_FILE_WRITE( fh, int_names, SIZE( int_names )*32, MPI_CHAR, status, ierr )
     2281          CALL MPI_FILE_WRITE( fh, int_names, SIZE( int_names ) * 32, MPI_CHAR, status, ierr )
    22342282          header_position = header_position + SIZE( int_names ) * 32
    22352283
     
    22402288!--       Character entries
    22412289          CALL MPI_FILE_SEEK( fh, header_position, MPI_SEEK_SET, ierr )
    2242           CALL MPI_FILE_WRITE( fh, text_lines, SIZE( text_lines )*128, MPI_CHAR, status, ierr )
     2290          CALL MPI_FILE_WRITE( fh, text_lines, SIZE( text_lines ) * 128, MPI_CHAR, status, ierr )
    22432291          header_position = header_position + SIZE( text_lines ) * 128
    22442292!
    22452293!---      REAL values
    22462294          CALL MPI_FILE_SEEK( fh, header_position, MPI_SEEK_SET, ierr )
    2247           CALL MPI_FILE_WRITE( fh, real_names, SIZE( real_names )*32, MPI_CHAR, status, ierr )
     2295          CALL MPI_FILE_WRITE( fh, real_names, SIZE( real_names ) * 32, MPI_CHAR, status, ierr )
    22482296          header_position = header_position + SIZE( real_names ) * 32
    22492297
     
    22542302!--       2d- and 3d- distributed array headers, all replicated array headers
    22552303          CALL MPI_FILE_SEEK( fh, header_position, MPI_SEEK_SET, ierr )
    2256           CALL MPI_FILE_WRITE( fh, array_names, SIZE( array_names )*32, MPI_CHAR, status, ierr )
     2304          CALL MPI_FILE_WRITE( fh, array_names, SIZE( array_names ) * 32, MPI_CHAR, status, ierr )
    22572305          header_position = header_position + SIZE( array_names ) * 32
    22582306
    22592307          CALL MPI_FILE_SEEK( fh, header_position, MPI_SEEK_SET, ierr )
    2260           CALL MPI_FILE_WRITE( fh, array_offset, SIZE( array_offset )*MPI_OFFSET_KIND, MPI_BYTE,  &
     2308          CALL MPI_FILE_WRITE( fh, array_offset, SIZE( array_offset ) * MPI_OFFSET_KIND, MPI_BYTE, &
    22612309                               status, ierr )  ! There is no I*8 datatype in Fortran
    22622310          header_position = header_position + SIZE( array_offset ) * rd_offset_kind
     
    23502398    IMPLICIT NONE
    23512399
    2352     INTEGER(iwp)                          ::  i            !<  loop index
    2353     INTEGER(KIND=rd_offset_kind)          ::  offset
    2354 
    2355     INTEGER(iwp), DIMENSION(1)            ::  dims1
    2356     INTEGER(iwp), DIMENSION(1)            ::  lize1
    2357     INTEGER(iwp), DIMENSION(1)            ::  start1
    2358     INTEGER(iwp), DIMENSION(0:numprocs-1) ::  lo_nr_val    !< local number of values in x and y direction
    2359     INTEGER(iwp), DIMENSION(0:numprocs-1) ::  all_nr_val   !< number of values for all PEs
    2360 
    2361     INTEGER, INTENT(IN), DIMENSION(nys:nyn,nxl:nxr)  ::  end_index
    2362     INTEGER, INTENT(OUT), DIMENSION(nys:nyn,nxl:nxr) ::  global_start
    2363     INTEGER, INTENT(IN), DIMENSION(nys:nyn,nxl:nxr)  ::  start_index
    2364 
    2365     LOGICAL, INTENT(OUT)                  ::  data_to_write  !< returns, if surface data have to be written
     2400    INTEGER(iwp)                          ::  i           !<  loop index
     2401    INTEGER(KIND=rd_offset_kind)          ::  offset      !<
     2402
     2403    INTEGER(iwp), DIMENSION(1)            ::  dims1       !<
     2404    INTEGER(iwp), DIMENSION(1)            ::  lize1       !<
     2405    INTEGER(iwp), DIMENSION(1)            ::  start1      !<
     2406
     2407    INTEGER(iwp), DIMENSION(0:numprocs-1) ::  all_nr_val  !< number of values for all PEs
     2408    INTEGER(iwp), DIMENSION(0:numprocs-1) ::  lo_nr_val   !< local number of values in x and y direction
     2409
     2410
     2411    INTEGER, INTENT(IN), DIMENSION(nys:nyn,nxl:nxr)  ::  end_index     !<
     2412    INTEGER, INTENT(OUT), DIMENSION(nys:nyn,nxl:nxr) ::  global_start  !<
     2413    INTEGER, INTENT(IN), DIMENSION(nys:nyn,nxl:nxr)  ::  start_index   !<
     2414
     2415    LOGICAL, INTENT(OUT) ::  data_to_write  !< returns, if surface data have to be written
    23662416
    23672417
     
    23722422    CALL MPI_ALLREDUCE( lo_nr_val, all_nr_val, numprocs, MPI_INTEGER, MPI_SUM, comm2d, ierr )
    23732423    IF ( ft_surf /= -1  .AND.  sm_io%iam_io_pe )  THEN
    2374        CALL MPI_TYPE_FREE( ft_surf, ierr )    ! if set, free last surface filetype
     2424       CALL MPI_TYPE_FREE( ft_surf, ierr )    ! If set, free last surface filetype
    23752425    ENDIF
    23762426
     
    24042454
    24052455!
    2406 !-- Actions during read
     2456!-- Actions during reading
    24072457    IF ( rd_flag )  THEN
    24082458       IF ( .NOT. ALLOCATED( m_start_index )  )  ALLOCATE( m_start_index(nys:nyn,nxl:nxr)  )
     
    24222472
    24232473!
    2424 !-- Actions during write
     2474!-- Actions during writing
    24252475    IF ( wr_flag )  THEN
    24262476!
     
    25022552    IMPLICIT NONE
    25032553
    2504     INTEGER, DIMENSION(2) ::  dims2
    2505     INTEGER, DIMENSION(2) ::  lize2
    2506     INTEGER, DIMENSION(2) ::  start2
    2507 
    2508     INTEGER, DIMENSION(3) ::  dims3
    2509     INTEGER, DIMENSION(3) ::  lize3
    2510     INTEGER, DIMENSION(3) ::  start3
     2554    INTEGER, DIMENSION(2) ::  dims2   !<
     2555    INTEGER, DIMENSION(2) ::  lize2   !<
     2556    INTEGER, DIMENSION(2) ::  start2  !<
     2557
     2558    INTEGER, DIMENSION(3) ::  dims3   !<
     2559    INTEGER, DIMENSION(3) ::  lize3   !<
     2560    INTEGER, DIMENSION(3) ::  start3  !<
    25112561
    25122562    TYPE(local_boundaries) ::  save_io_grid  !< temporary variable to store grid settings
     
    26562706! Description:
    26572707! ------------
    2658 !> This subroutine creates file types to access 3d-soil arrays
    2659 !> distributed in blocks among processes to a single file that contains the global arrays.
    2660 !> It is not required for the serial mode.
     2708!> This subroutine creates file types to access 3d-soil arrays distributed in blocks among processes
     2709!> to a single file that contains the global arrays. It is not required for the serial mode.
    26612710!--------------------------------------------------------------------------------------------------!
    26622711#if defined( __parallel )
     
    26652714    IMPLICIT NONE
    26662715
    2667     INTEGER, INTENT(IN)  ::  nzb_soil
    2668     INTEGER, INTENT(IN)  ::  nzt_soil
    2669 
    2670     INTEGER, DIMENSION(3) ::  dims3
    2671     INTEGER, DIMENSION(3) ::  lize3
    2672     INTEGER, DIMENSION(3) ::  start3
     2716    INTEGER, INTENT(IN)   ::  nzb_soil  !<
     2717    INTEGER, INTENT(IN)   ::  nzt_soil  !<
     2718
     2719    INTEGER, DIMENSION(3) ::  dims3     !<
     2720    INTEGER, DIMENSION(3) ::  lize3     !<
     2721    INTEGER, DIMENSION(3) ::  start3    !<
    26732722
    26742723
     
    27642813    IMPLICIT NONE
    27652814
    2766     INTEGER(iwp) ::  i
     2815    INTEGER(iwp) ::  i  !<
    27672816
    27682817
     
    28162865    IMPLICIT NONE
    28172866
    2818     INTEGER, INTENT(out)                   ::  i_endian
    2819     INTEGER(KIND=8), TARGET                ::  int8
    2820 
    2821     INTEGER, DIMENSION(1)                  ::  bufshape
    2822     INTEGER(KIND=4), POINTER, DIMENSION(:) ::  int4
    2823 
    2824     TYPE(C_PTR)                            ::  ptr
     2867    INTEGER, INTENT(out)                   ::  i_endian  !<
     2868    INTEGER(KIND=8), TARGET                ::  int8      !<
     2869
     2870    INTEGER, DIMENSION(1)                  ::  bufshape  !<
     2871    INTEGER(KIND=4), POINTER, DIMENSION(:) ::  int4      !<
     2872
     2873    TYPE(C_PTR)                            ::  ptr       !<
    28252874
    28262875
     
    28322881
    28332882    IF ( int4(1) == 1 )  THEN
    2834        i_endian = 1    ! little endian
     2883       i_endian = 1    ! Little endian
    28352884    ELSE
    2836        i_endian = 2    ! big endian
     2885       i_endian = 2    ! Big endian
    28372886    ENDIF
    28382887
  • palm/trunk/SOURCE/run_control.f90

    r4360 r4591  
    11!> @file run_control.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:
    2121! -----------------
    22 ! 
    23 ! 
     22!
     23!
    2424! Former 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! 3655 2019-01-07 16:51:22Z knoop
    3034! Corrected "Former revisions" section
     
    3741! ------------
    3842!> Computation and output of run-control quantities
    39 !------------------------------------------------------------------------------!
     43!--------------------------------------------------------------------------------------------------!
    4044 SUBROUTINE run_control
    41  
    4245
    43     USE cpulog,                                                                &
    44         ONLY:  cpu_log, log_point
    4546
    46     USE control_parameters,                                                    &
    47         ONLY:  advected_distance_x, advected_distance_y,                       &
    48                current_timestep_number, disturbance_created, dt_3d, mgcycles,  &
    49                run_control_header, runnr, simulated_time, simulated_time_chr,  &
     47    USE cpulog,                                                                                    &
     48        ONLY:  cpu_log,                                                                            &
     49               log_point
     50
     51    USE control_parameters,                                                                        &
     52        ONLY:  advected_distance_x,                                                                &
     53               advected_distance_y,                                                                &
     54               current_timestep_number,                                                            &
     55               disturbance_created,                                                                &
     56               dt_3d,                                                                              &
     57               mgcycles,                                                                           &
     58               run_control_header,                                                                 &
     59               runnr,                                                                              &
     60               simulated_time,                                                                     &
     61               simulated_time_chr,                                                                 &
    5062               timestep_reason
    5163
    52     USE indices,                                                               &
     64    USE indices,                                                                                   &
    5365        ONLY:  nzb
    5466
     
    5769    USE pegrid
    5870
    59     USE statistics,                                                            &
    60         ONLY:  flow_statistics_called, hom, pr_palm, u_max, u_max_ijk, v_max,  &
    61                v_max_ijk, w_max, w_max_ijk
     71    USE statistics,                                                                                &
     72        ONLY:  flow_statistics_called,                                                             &
     73               hom,                                                                                &
     74               pr_palm,                                                                            &
     75               u_max,                                                                              &
     76               u_max_ijk,                                                                          &
     77               v_max,                                                                              &
     78               v_max_ijk,                                                                          &
     79               w_max,                                                                              &
     80               w_max_ijk
    6281
    6382    IMPLICIT NONE
    6483
    65     CHARACTER (LEN=1) ::  disturb_chr
     84    CHARACTER (LEN=1) ::  disturb_chr  !<
    6685
    6786!
     
    7897
    7998!
    80 !--    Check, whether file unit is already open (may have been opened in header
    81 !--    before)
     99!--    Check, whether file unit is already open (may have been opened in header before)
    82100       CALL check_open( 15 )
    83101
     
    96114          disturb_chr = ' '
    97115       ENDIF
    98        WRITE ( 15, 101 )  runnr, current_timestep_number, simulated_time_chr,  &
    99                           INT( ( simulated_time-INT( simulated_time ) ) * 100),&
    100                           dt_3d, timestep_reason, u_max, disturb_chr,          &
    101                           v_max, disturb_chr, w_max, hom(nzb,1,pr_palm,0),     &
    102                           hom(nzb+8,1,pr_palm,0), hom(nzb+3,1,pr_palm,0),      &
    103                           hom(nzb+6,1,pr_palm,0), hom(nzb+4,1,pr_palm,0),      &
    104                           hom(nzb+5,1,pr_palm,0), hom(nzb+9,1,pr_palm,0),      &
    105                           hom(nzb+10,1,pr_palm,0), u_max_ijk(1:3),             &
    106                           v_max_ijk(1:3), w_max_ijk(1:3),                      &
    107                           advected_distance_x/1000.0_wp,                       &
     116       WRITE ( 15, 101 )  runnr, current_timestep_number, simulated_time_chr,                      &
     117                          INT( ( simulated_time-INT( simulated_time ) ) * 100),                    &
     118                          dt_3d, timestep_reason, u_max, disturb_chr, v_max, disturb_chr, w_max,   &
     119                          hom(nzb,1,pr_palm,0), hom(nzb+8,1,pr_palm,0), hom(nzb+3,1,pr_palm,0),    &
     120                          hom(nzb+6,1,pr_palm,0), hom(nzb+4,1,pr_palm,0), hom(nzb+5,1,pr_palm,0),  &
     121                          hom(nzb+9,1,pr_palm,0), hom(nzb+10,1,pr_palm,0), u_max_ijk(1:3),         &
     122                          v_max_ijk(1:3), w_max_ijk(1:3), advected_distance_x/1000.0_wp,           &
    108123                          advected_distance_y/1000.0_wp, mgcycles
    109124!
     
    113128    ENDIF
    114129!
    115 !-- If required, reset disturbance flag. This has to be done outside the above 
     130!-- If required, reset disturbance flag. This has to be done outside the above
    116131!-- IF-loop, because the flag would otherwise only be reset on PE0
    117132    IF ( disturbance_created )  disturbance_created = .FALSE.
     
    131146          &'----------------------------------------------------------------', &
    132147          &'---------')
    133 101 FORMAT (I3,1X,I6,1X,A8,'.',I2.2,1X,F8.4,A1,1X,F8.4,A1,F8.4,A1,F8.4,1X,     &
    134             F6.3,1X,F5.2, &
     148101 FORMAT (I3,1X,I6,1X,A8,'.',I2.2,1X,F8.4,A1,1X,F8.4,A1,F8.4,A1,F8.4,1X, F6.3,1X,F5.2,           &
    135149            2X,E10.3,2X,F6.0,1X,4(E10.3,1X),3(3(I4),1X),F8.3,1X,F8.3,5X,I3)
    136150
  • palm/trunk/SOURCE/shared_memory_io_mod.f90

    r4536 r4591  
    11!> @file shared_memory_io_mod.f90
    2 !------------------------------------------------------------------------------!
    3 ! This file is part of PALM.
    4 !
    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/>.
     2!--------------------------------------------------------------------------------------------------!
     3! This file is part of the PALM model system.
     4!
     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:
    2121! -----------------
    22 ! 
    23 ! 
     22!
     23!
    2424! Former revisions:
    2525! $Id$
    2626!
     27! File re-formatted to follow the PALM coding standard
     28!
     29!
     30!
    2731! Initial version (Klaus Ketelsen)
    2832!
    29 ! 
    30 !
    31 ! Description:
    32 ! ------------
    33 !> handle MPI-IO or NetCDF-IO shared memory arrays.
    34 !> This module performs the organization of new communicators, adapted PE-grids
    35 !> and allocation of shared memory arrays. The IO itself is not done here.
    36 !------------------------------------------------------------------------------!
     33!
     34!
     35! Description:
     36! ------------
     37!> Handle MPI-IO or NetCDF-IO shared memory arrays.
     38!> This module performs the organization of new communicators, adapted PE-grids and allocation of
     39!> shared memory arrays. The IO itself is not done here.
     40!--------------------------------------------------------------------------------------------------!
    3741 MODULE shared_memory_io_mod
    3842
     
    4852
    4953    USE control_parameters,                                                                        &
    50         ONLY:  maximum_grid_level, mg_switch_to_pe0_level, message_string
     54        ONLY: maximum_grid_level,                                                                  &
     55              message_string,                                                                      &
     56              mg_switch_to_pe0_level
     57
    5158
    5259    USE indices,                                                                                   &
    53         ONLY: nbgp, nnx, nny, nnz, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb, nzt
     60        ONLY: nbgp,                                                                                &
     61              nnx,                                                                                 &
     62              nny,                                                                                 &
     63              nnz,                                                                                 &
     64              nx,                                                                                  &
     65              nxl,                                                                                 &
     66              nxlg,                                                                                &
     67              nxr,                                                                                 &
     68              nxrg,                                                                                &
     69              ny,                                                                                  &
     70              nyn,                                                                                 &
     71              nyng,                                                                                &
     72              nys,                                                                                 &
     73              nysg,                                                                                &
     74              nzb,                                                                                 &
     75              nzt
    5476
    5577    USE kinds,                                                                                     &
    56         ONLY: wp, iwp
     78        ONLY: iwp,                                                                                 &
     79              wp
     80
    5781
    5882    USE transpose_indices,                                                                         &
    59         ONLY:  nys_x, nyn_x, nys_z, nyn_z, nxl_z, nxr_z
     83        ONLY: nxl_z,                                                                               &
     84              nxr_z,                                                                               &
     85              nyn_x,                                                                               &
     86              nyn_z,                                                                               &
     87              nys_x,                                                                               &
     88              nys_z
     89
     90
    6091
    6192    USE pegrid,                                                                                    &
    62         ONLY: comm1dx, comm1dy, comm2d, ierr, myid, myidx, myidy, npex, npey, numprocs, pdims,     &
    63               pleft, pnorth, pright, psouth, sendrecvcount_xy
     93        ONLY: comm1dx,                                                                             &
     94              comm1dy,                                                                             &
     95              comm2d,                                                                              &
     96              ierr,                                                                                &
     97              myid,                                                                                &
     98              myidx,                                                                               &
     99              myidy,                                                                               &
     100              npex,                                                                                &
     101              npey,                                                                                &
     102              numprocs,                                                                            &
     103              pdims,                                                                               &
     104              pleft,                                                                               &
     105              pnorth,                                                                              &
     106              pright,                                                                              &
     107              psouth,                                                                              &
     108              sendrecvcount_xy
     109
    64110#if defined( __parallel )
    65111    USE pegrid,                                                                                    &
    66         ONLY: pcoord, reorder
     112        ONLY: pcoord,                                                                              &
     113              reorder
    67114#endif
    68115
     
    75122!
    76123!-- Type to store grid information
    77     TYPE, PUBLIC ::  local_boundaries
    78 
    79        INTEGER(iwp) ::  nxl
    80        INTEGER(iwp) ::  nxr
    81        INTEGER(iwp) ::  nys
    82        INTEGER(iwp) ::  nyn
    83        INTEGER(iwp) ::  nnx
    84        INTEGER(iwp) ::  nny
    85        INTEGER(iwp) ::  nx
    86        INTEGER(iwp) ::  ny
     124    TYPE, PUBLIC ::  local_boundaries  !<
     125
     126       INTEGER(iwp) ::  nnx  !<
     127       INTEGER(iwp) ::  nny  !<
     128       INTEGER(iwp) ::  nx   !<
     129       INTEGER(iwp) ::  nxl  !<
     130       INTEGER(iwp) ::  nxr  !<
     131       INTEGER(iwp) ::  ny   !<
     132       INTEGER(iwp) ::  nyn  !<
     133       INTEGER(iwp) ::  nys  !<
     134
     135
     136
    87137
    88138    END TYPE local_boundaries
     
    91141!-- Class definition for shared memory instances.
    92142!-- For every use of shared memory IO, one instance of this class is created.
    93     TYPE, PUBLIC ::  sm_class
    94 
    95        INTEGER(iwp)         :: nr_io_pe_per_node = 2   !< typical configuration, 2 sockets per node
    96        LOGICAL              :: no_shared_Memory_in_this_run
     143    TYPE, PUBLIC ::  sm_class  !<
     144
     145       INTEGER(iwp) ::  nr_io_pe_per_node = 2         !< typical configuration, 2 sockets per node
     146       LOGICAL      ::  no_shared_Memory_in_this_run  !<
    97147!
    98148!--    Variables for the shared memory communicator
    99        INTEGER(iwp), PUBLIC ::  comm_shared    !< Communicator for processes with shared array
    100        INTEGER(iwp), PUBLIC ::  sh_npes
    101        INTEGER(iwp), PUBLIC ::  sh_rank
     149       INTEGER(iwp), PUBLIC ::  comm_shared   !< Communicator for processes with shared array
     150       INTEGER(iwp), PUBLIC ::  sh_npes       !<
     151       INTEGER(iwp), PUBLIC ::  sh_rank       !<
    102152
    103153       LOGICAL, PUBLIC ::  iam_io_pe = .TRUE.  !< This PE is an IO-PE
    104154!
    105155!--    Variables for the I/O virtual grid
    106        INTEGER(iwp), PUBLIC ::  comm_io        !< Communicator for all IO processes
    107        INTEGER(iwp), PUBLIC ::  io_npes
    108        INTEGER(iwp), PUBLIC ::  io_rank
     156       INTEGER(iwp), PUBLIC ::  comm_io  !< Communicator for all IO processes
     157       INTEGER(iwp), PUBLIC ::  io_npes  !<
     158       INTEGER(iwp), PUBLIC ::  io_rank  !<
    109159
    110160       TYPE( local_boundaries ), PUBLIC ::  io_grid
     
    112162!
    113163!--    Variables for the node local communicator
    114        INTEGER(iwp)         ::  comm_node      !< Communicator for all processes of current node
    115        INTEGER(iwp)         ::  io_pe_global_rank
    116        INTEGER(iwp)         ::  n_npes
    117        INTEGER(iwp)         ::  n_rank
    118 
    119        CONTAINS
    120 
    121           PRIVATE
    122 
    123           PROCEDURE, PASS(this), PUBLIC ::  is_sm_active
    124           PROCEDURE, PASS(this), PUBLIC ::  sm_adjust_outer_boundary
    125           PROCEDURE, PASS(this), PUBLIC ::  sm_free_shared
    126           PROCEDURE, PASS(this), PUBLIC ::  sm_init_comm
    127           PROCEDURE, PASS(this), PUBLIC ::  sm_node_barrier
    128 #if defined( __parallel )
    129           PROCEDURE, PASS(this), PUBLIC ::  sm_allocate_shared_1d
    130           PROCEDURE, PASS(this), PUBLIC ::  sm_allocate_shared_2d
    131           PROCEDURE, PASS(this), PUBLIC ::  sm_allocate_shared_2di
    132           PROCEDURE, PASS(this), PUBLIC ::  sm_allocate_shared_3d
     164       INTEGER(iwp) ::  comm_node          !< Communicator for all processes of current node
     165       INTEGER(iwp) ::  io_pe_global_rank  !<
     166       INTEGER(iwp) ::  n_npes             !<
     167       INTEGER(iwp) ::  n_rank             !<
     168
     169 CONTAINS
     170
     171       PRIVATE
     172
     173          PROCEDURE, PASS(this), PUBLIC ::  is_sm_active              !<
     174          PROCEDURE, PASS(this), PUBLIC ::  sm_adjust_outer_boundary  !<
     175          PROCEDURE, PASS(this), PUBLIC ::  sm_free_shared            !<
     176          PROCEDURE, PASS(this), PUBLIC ::  sm_init_comm              !<
     177          PROCEDURE, PASS(this), PUBLIC ::  sm_node_barrier           !<
     178#if defined( __parallel )
     179          PROCEDURE, PASS(this), PUBLIC ::  sm_allocate_shared_1d   !<
     180          PROCEDURE, PASS(this), PUBLIC ::  sm_allocate_shared_2d   !<
     181          PROCEDURE, PASS(this), PUBLIC ::  sm_allocate_shared_2di  !<
     182          PROCEDURE, PASS(this), PUBLIC ::  sm_allocate_shared_3d   !<
    133183
    134184          GENERIC, PUBLIC ::  sm_allocate_shared =>  sm_allocate_shared_1d, sm_allocate_shared_2d, &
    135                                                      sm_allocate_shared_2di, sm_allocate_shared_3d
     185                                                  sm_allocate_shared_2di, sm_allocate_shared_3d  !<
    136186#endif
    137187    END TYPE sm_class
    138188
    139189
    140 CONTAINS
     190 CONTAINS
    141191
    142192
     
    154204
    155205#if defined( __parallel )
    156     INTEGER             :: color
    157     INTEGER             :: max_n_npes  !< Maximum number of PEs/node
     206    INTEGER ::  color       !<
     207    INTEGER :: max_n_npes  !< Maximum number of PEs/node
    158208#endif
    159209
     
    244294    IMPLICIT NONE
    245295
    246     INTEGER(iwp), INTENT(OUT) :: color
    247 
    248     INTEGER(iwp) ::  group_start
    249     INTEGER(iwp) ::  n
    250     INTEGER(iwp) ::  my_color
    251     INTEGER(iwp) ::  pe
    252     INTEGER(iwp) ::  sh_group_size
    253 
    254     INTEGER(iwp), DIMENSION(4,0:this%n_npes-1) ::  local_dim_s
    255     INTEGER(iwp), DIMENSION(4,0:this%n_npes-1) ::  local_dim_r
    256 
    257     TYPE(local_boundaries), DIMENSION(32) ::  node_grid
    258 
    259 !
    260 !-- Nn shared memory I/O on one node jobs
     296    INTEGER(iwp), INTENT(OUT) ::  color  !<
     297
     298    INTEGER(iwp) ::  group_start    !<
     299    INTEGER(iwp) ::  my_color       !<
     300    INTEGER(iwp) ::  n              !<
     301    INTEGER(iwp) ::  pe             !<
     302    INTEGER(iwp) ::  sh_group_size  !<
     303
     304    INTEGER(iwp), DIMENSION(4,0:this%n_npes-1) ::  local_dim_s   !<
     305    INTEGER(iwp), DIMENSION(4,0:this%n_npes-1) ::  local_dim_r   !<
     306
     307    TYPE(local_boundaries), DIMENSION(32) ::  node_grid  !<
     308
     309!
     310!-- No shared memory I/O on one node jobs
    261311    IF ( numprocs < this%n_npes )  THEN
    262312       this%no_shared_memory_in_this_run = .TRUE.
     
    337387!> Function to return if shared Memory IO is active.
    338388!--------------------------------------------------------------------------------------------------!
    339  FUNCTION is_sm_active( this) RESULT( ac )
    340 
    341     IMPLICIT NONE
    342 
    343     CLASS(sm_class), INTENT(inout) ::  this
    344 
    345     LOGICAL ::  ac
     389 FUNCTION is_sm_active( this ) RESULT( ac )
     390
     391    IMPLICIT NONE
     392
     393    CLASS(sm_class), INTENT(inout) ::  this  !<
     394
     395    LOGICAL ::  ac  !<
    346396
    347397    ac = .NOT. this%no_shared_memory_in_this_run
     
    360410    IMPLICIT NONE
    361411
    362     CLASS(sm_class), INTENT(inout)  ::  this
    363 
    364     INTEGER(iwp)                    ::  disp_unit
    365     INTEGER(iwp), INTENT(IN)        ::  d1
    366     INTEGER(iwp), INTENT(IN)        ::  d2
    367     INTEGER(iwp), SAVE              ::  pe_from = 0
    368     INTEGER(KIND=MPI_ADDRESS_KIND)  ::  rem_size
    369     INTEGER(iwp), INTENT(OUT)       ::  win
    370     INTEGER(KIND=MPI_ADDRESS_KIND)  ::  wsize
    371 
    372     INTEGER, DIMENSION(1)           ::  buf_shape
    373 
    374     REAL(wp), DIMENSION(:), POINTER ::  buf
    375     REAL(wp), DIMENSION(:), POINTER ::  p1
    376 
    377     TYPE(C_PTR), SAVE               ::  base_ptr
    378     TYPE(C_PTR), SAVE               ::  rem_ptr
     412    CLASS(sm_class), INTENT(inout) ::  this         !<
     413                                                    !<
     414    INTEGER(iwp)                   ::  disp_unit    !<
     415    INTEGER(iwp), INTENT(IN)       ::  d1           !<
     416    INTEGER(iwp), INTENT(IN)       ::  d2           !<
     417    INTEGER(iwp), SAVE             ::  pe_from = 0  !<
     418    INTEGER(KIND=MPI_ADDRESS_KIND) ::  rem_size     !<
     419    INTEGER(iwp), INTENT(OUT)      ::  win          !<
     420    INTEGER(KIND=MPI_ADDRESS_KIND) ::  wsize        !<
     421
     422    INTEGER, DIMENSION(1)           ::  buf_shape   !<
     423
     424    REAL(wp), DIMENSION(:), POINTER ::  buf         !<
     425    REAL(wp), DIMENSION(:), POINTER ::  p1          !<
     426
     427    TYPE(C_PTR), SAVE               ::  base_ptr    !<
     428    TYPE(C_PTR), SAVE               ::  rem_ptr     !<
    379429
    380430
     
    415465    IMPLICIT NONE
    416466
    417     CLASS(sm_class), INTENT(INOUT)    ::  this
    418 
    419     INTEGER(iwp)                      ::  disp_unit
    420     INTEGER(iwp), INTENT(IN)          ::  n_nxlg
    421     INTEGER(iwp), INTENT(IN)          ::  n_nxrg
    422     INTEGER(iwp), INTENT(IN)          ::  n_nyng
    423     INTEGER(iwp), INTENT(IN)          ::  n_nysg
    424     INTEGER(iwp), SAVE                ::  pe_from = 0
    425     INTEGER(KIND=MPI_ADDRESS_KIND)    ::  rem_size
    426     INTEGER(iwp), INTENT(OUT)         ::  win
    427     INTEGER(KIND=MPI_ADDRESS_KIND)    ::  wsize
    428 
    429     INTEGER(iwp), DIMENSION(2)        ::  buf_shape
    430 
    431     REAL(wp), DIMENSION(:,:), POINTER ::  buf
    432     REAL(wp), DIMENSION(:,:), POINTER ::  p2
    433 
    434     TYPE(C_PTR),SAVE                  ::  base_ptr
    435     TYPE(C_PTR),SAVE                  ::  rem_ptr
     467    CLASS(sm_class), INTENT(INOUT)    ::  this         !<
     468
     469    INTEGER(iwp)                      ::  disp_unit    !<
     470    INTEGER(iwp), INTENT(IN)          ::  n_nxlg       !<
     471    INTEGER(iwp), INTENT(IN)          ::  n_nxrg       !<
     472    INTEGER(iwp), INTENT(IN)          ::  n_nyng       !<
     473    INTEGER(iwp), INTENT(IN)          ::  n_nysg       !<
     474    INTEGER(iwp), SAVE                ::  pe_from = 0  !<
     475    INTEGER(KIND=MPI_ADDRESS_KIND)    ::  rem_size     !<
     476    INTEGER(iwp), INTENT(OUT)         ::  win          !<
     477    INTEGER(KIND=MPI_ADDRESS_KIND)    ::  wsize        !<
     478
     479    INTEGER(iwp), DIMENSION(2)        ::  buf_shape    !<
     480
     481    REAL(wp), DIMENSION(:,:), POINTER ::  buf          !<
     482    REAL(wp), DIMENSION(:,:), POINTER ::  p2           !<
     483
     484    TYPE(C_PTR),SAVE                  ::  base_ptr     !<
     485    TYPE(C_PTR),SAVE                  ::  rem_ptr      !<
    436486
    437487
     
    474524    IMPLICIT NONE
    475525
    476     CLASS(sm_class), INTENT(inout)      :: this
    477 
    478     INTEGER(iwp)                          ::  disp_unit
    479     INTEGER(iwp), INTENT(IN)              ::  n_nxlg
    480     INTEGER(iwp), INTENT(IN)              ::  n_nxrg
    481     INTEGER(iwp), INTENT(IN)              ::  n_nyng
    482     INTEGER(iwp), INTENT(IN)              ::  n_nysg
    483     INTEGER(iwp), SAVE                    ::  pe_from = 0
    484     INTEGER(kind=MPI_ADDRESS_KIND)        ::  rem_size
    485     INTEGER(iwp), INTENT(OUT)             ::  win
    486     INTEGER(kind=MPI_ADDRESS_KIND)        ::  wsize
    487 
    488     INTEGER(iwp), DIMENSION(2)            ::  buf_shape
    489 
    490     INTEGER(iwp), DIMENSION(:,:), POINTER ::  buf
    491     INTEGER(iwp), DIMENSION(:,:), POINTER ::  p2i
    492 
    493     TYPE(C_PTR),SAVE                      ::  base_ptr
    494     TYPE(C_PTR),SAVE                      ::  rem_ptr
     526    CLASS(sm_class), INTENT(inout)        ::  this         !<
     527
     528    INTEGER(iwp)                          ::  disp_unit    !<
     529    INTEGER(iwp), INTENT(IN)              ::  n_nxlg       !<
     530    INTEGER(iwp), INTENT(IN)              ::  n_nxrg       !<
     531    INTEGER(iwp), INTENT(IN)              ::  n_nyng       !<
     532    INTEGER(iwp), INTENT(IN)              ::  n_nysg       !<
     533    INTEGER(iwp), SAVE                    ::  pe_from = 0  !<
     534    INTEGER(kind=MPI_ADDRESS_KIND)        ::  rem_size     !<
     535    INTEGER(iwp), INTENT(OUT)             ::  win          !<
     536    INTEGER(kind=MPI_ADDRESS_KIND)        ::  wsize        !<
     537
     538    INTEGER(iwp), DIMENSION(2)            ::  buf_shape    !<
     539
     540    INTEGER(iwp), DIMENSION(:,:), POINTER ::  buf          !<
     541    INTEGER(iwp), DIMENSION(:,:), POINTER ::  p2i          !<
     542
     543    TYPE(C_PTR),SAVE                      ::  base_ptr     !<
     544    TYPE(C_PTR),SAVE                      ::  rem_ptr      !<
    495545
    496546
     
    533583    IMPLICIT NONE
    534584
    535     CLASS(sm_class), INTENT(inout)      ::  this
    536 
    537     INTEGER                             ::  disp_unit
    538     INTEGER, INTENT(IN)                 ::  d1e
    539     INTEGER, INTENT(IN)                 ::  d1s
    540     INTEGER, INTENT(IN)                 ::  d2e
    541     INTEGER, INTENT(IN)                 ::  d2s
    542     INTEGER, INTENT(IN)                 ::  d3e
    543     INTEGER, INTENT(IN)                 ::  d3s
    544     INTEGER, SAVE                       ::  pe_from = 0
    545     INTEGER(KIND=MPI_ADDRESS_KIND)      ::  rem_size
    546     INTEGER, INTENT(OUT)                ::  win
    547     INTEGER(KIND=MPI_ADDRESS_KIND)      ::  wsize
    548 
    549     INTEGER, DIMENSION(3)               ::  buf_shape
    550 
    551     REAL(wp), DIMENSION(:,:,:), POINTER ::  buf
    552     REAL(wp), DIMENSION(:,:,:), POINTER ::  p3
    553 
    554     TYPE(C_PTR), SAVE                   ::  base_ptr
    555     TYPE(C_PTR), SAVE                   ::  rem_ptr
     585    CLASS(sm_class), INTENT(inout)      ::  this         !<
     586
     587    INTEGER                             ::  disp_unit    !<
     588    INTEGER, INTENT(IN)                 ::  d1e          !<
     589    INTEGER, INTENT(IN)                 ::  d1s          !<
     590    INTEGER, INTENT(IN)                 ::  d2e          !<
     591    INTEGER, INTENT(IN)                 ::  d2s          !<
     592    INTEGER, INTENT(IN)                 ::  d3e          !<
     593    INTEGER, INTENT(IN)                 ::  d3s          !<
     594    INTEGER, SAVE                       ::  pe_from = 0  !<
     595    INTEGER(KIND=MPI_ADDRESS_KIND)      ::  rem_size     !<
     596    INTEGER, INTENT(OUT)                ::  win          !<
     597    INTEGER(KIND=MPI_ADDRESS_KIND)      ::  wsize        !<
     598
     599    INTEGER, DIMENSION(3)               ::  buf_shape    !<
     600
     601    REAL(wp), DIMENSION(:,:,:), POINTER ::  buf          !<
     602    REAL(wp), DIMENSION(:,:,:), POINTER ::  p3           !<
     603
     604    TYPE(C_PTR), SAVE                   ::  base_ptr     !<
     605    TYPE(C_PTR), SAVE                   ::  rem_ptr      !<
    556606
    557607
     
    596646    IMPLICIT NONE
    597647
    598     CLASS(sm_class), INTENT(inout) ::  this
     648    CLASS(sm_class), INTENT(inout) ::  this  !<
    599649
    600650
     
    640690    IMPLICIT NONE
    641691
    642     CLASS(sm_class), INTENT(inout) ::  this
    643 
    644     INTEGER(iwp), INTENT(INOUT)    ::  win
     692    CLASS(sm_class), INTENT(inout) ::  this  !<
     693
     694    INTEGER(iwp), INTENT(INOUT)    ::  win   !<
    645695
    646696    IF ( this%no_shared_memory_in_this_run  .OR.  win == -1234567890 )  RETURN
     
    662712    IMPLICIT NONE
    663713
    664     CLASS(sm_class), INTENT(inout)     :: this
     714    CLASS(sm_class), INTENT(inout) ::  this  !<
    665715
    666716
  • palm/trunk/SOURCE/singleton_mod.f90

    r4182 r4591  
    11!> @file singleton_mod.f90
    2 !------------------------------------------------------------------------------!
     2!--------------------------------------------------------------------------------------------------!
     3! This file is part of the PALM model system.
     4!
     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!
     19!
    320! Current revisions:
    421! -----------------
    5 ! 
    6 ! 
     22!
     23!
    724! Former revisions:
    825! -----------------
    926! $Id$
     27! File re-formatted to follow the PALM coding standard
     28!
     29!
     30! 4182 2019-08-22 15:20:23Z scharf
    1031! Corrected "Former revisions" section
    11 ! 
     32!
    1233! 3761 2019-02-25 15:31:42Z raasch
    13 ! statement added to prevent compiler warning about unused variables
     34! Statement added to prevent compiler warning about unused variables
    1435!
    1536! Revision 1.1  2002/05/02 18:56:59  raasch
     
    1738!
    1839!
     40!--------------------------------------------------------------------------------------------------!
    1941! Description:
    2042! ------------
    2143!> Multivariate Fast Fourier Transform
    2244!>
    23 !> Fortran 90 Implementation of Singleton's mixed-radix algorithm,
    24 !> RC Singleton, Stanford Research Institute, Sept. 1968.
    25 !>
    26 !> Adapted from fftn.c, translated from Fortran 66 to C by Mark Olesen and
    27 !> John Beale.
    28 !>
    29 !> Fourier transforms can be computed either in place, using assumed size
    30 !> arguments, or by generic function, using assumed shape arguments.
    31 !>
     45!> Fortran 90 Implementation of Singleton's mixed-radix algorithm, RC Singleton, Stanford Research
     46!> Institute, Sept. 1968.
     47!>
     48!> Adapted from fftn.c, translated from Fortran 66 to C by Mark Olesen and John Beale.
     49!>
     50!> Fourier transforms can be computed either in place, using assumed size arguments, or by generic
     51!> function, using assumed shape arguments.
    3252!>
    3353!> Public:
    3454!>
    35 !>   fftkind                              kind parameter of complex arguments
    36 !>                                        and function results.
     55!>   fftkind                              kind parameter of complex arguments and function results.
    3756!>
    3857!>   fft(array, dim, inv, stat)           generic transform function
     
    5271!> Formal Parameters:
    5372!>
    54 !>   array    The complex array to be transformed. array can be of arbitrary
    55 !>            rank (i.e. up to seven).
    56 !>
    57 !>   shape    With subroutine fftn, the shape of the array to be transformed
    58 !>            has to be passed separately, since fftradix - the internal trans-
    59 !>            formation routine - will treat array always as one dimensional.
    60 !>            The product of elements in shape must be the number of
     73!>   array    The complex array to be transformed. Array can be of arbitrary rank (i.e. up to seven).
     74!>
     75!>   shape    With subroutine fftn, the shape of the array to be transformed has to be passed
     76!>            separately, since fftradix - the internal transformation routine - will always treat
     77!>            array as one dimensional. The product of elements in shape must be the number of
    6178!>            elements in array.
    62 !>            Although passing array with assumed shape would have been nicer,
    63 !>            I prefered assumed size in order to prevent the compiler from
    64 !>            using a copy-in-copy-out mechanism. That would generally be
    65 !>            necessary with fftn passing array to fftradix and with fftn
    66 !>            being prepared for accepting non consecutive array sections.
    67 !>            Using assumed size, it's up to the user to pass an array argu-
    68 !>            ment, that can be addressed as continous one dimensional array
    69 !>            without copying. Otherwise, transformation will not really be
    70 !>            performed in place.
    71 !>            On the other hand, since the rank of array and the size of
    72 !>            shape needn't match, fftn is appropriate for handling more than
    73 !>            seven dimensions.
    74 !>            As far as function fft is concerned all this doesn't matter,
    75 !>            because the argument will be copied anyway. Thus no extra
    76 !>            shape argument is needed for fft.
     79!>            Although passing array with assumed shape would have been nicer, I prefered assumed
     80!>            size in order to prevent the compiler from using a copy-in-copy-out mechanism. That
     81!>            would generally be necessary with fftn passing array to fftradix and with fftn being
     82!>            prepared for accepting non consecutive array sections. Using assumed size, it's up to
     83!>            the user to pass an array argument, that can be addressed as continous one dimensional
     84!>            array without copying. Otherwise, transformation will not really be performed in place.
     85!>            On the other hand, since the rank of array and the size of shape needn't match, fftn
     86!>            is appropriate for handling more than seven dimensions. As far as function fft is
     87!>            concerned all this doesn't matter, because the argument will be copied anyway. Thus no
     88!>            extra shape argument is needed for fft.
    7789!>
    7890!> Optional Parameters:
    7991!>
    80 !>   dim      One dimensional integer array, containing the dimensions to be
    81 !>            transformed. Default is (/1,...,N/) with N being the rank of
    82 !>            array, i.e. complete transform. dim can restrict transformation
    83 !>            to a subset of available dimensions. Its size must not exceed the
    84 !>            rank of array or the size of shape respectivly.
    85 !>
    86 !>   inv      If .true., inverse transformation will be performed. Default is
    87 !>            .false., i.e. forward transformation.
    88 !>
    89 !>   stat     If present, a system dependent nonzero status value will be
    90 !>            returned in stat, if allocation of temporary storage failed.
     92!>   dim      One dimensional integer array, containing the dimensions to be transformed. Default
     93!>            is (/1,...,N/) with N being the rank of array, i.e. complete transform. dim can
     94!>            restrict transformation to a subset of available dimensions. Its size must not exceed
     95!>            the rank of array or the size of shape respectivly.
     96!>
     97!>   inv      If .true., inverse transformation will be performed. Default is .false., i.e. forward
     98!>            transformation.
     99!>
     100!>   stat     If present, a system dependent nonzero status value will be returned in stat, if
     101!>            allocation of temporary storage failed.
    91102!>
    92103!>
    93104!> Scaling:
    94105!>
    95 !>   Transformation results will always be scaled by the square root of the
    96 !>   product of sizes of each dimension in dim. (See examples below)
     106!>   Transformation results will always be scaled by the square root of the product of sizes of each
     107!>   dimension in dim. (See examples below)
    97108!>
    98109!>
     
    111122!>     result = fft(A, dim=(/1,3/))
    112123!>
    113 !>   will transform with respect to the first and the third dimension, scaled
    114 !>   by sqrt(L*N).
     124!>   will transform with respect to the first and the third dimension, scaled by sqrt(L*N).
    115125!>
    116126!>     result = fft(fft(A), inv=.true.)
     
    127137!>
    128138!>   Following changes have been introduced with respect to fftn.c:
    129 !>   - complex arguments and results are of type complex, rather than
    130 !>     real an imaginary part separately.
    131 !>   - increment parameter (magnitude of isign) has been dropped,
    132 !>     inc is always one, direction of transform is given by inv.     
    133 !>   - maxf and maxp have been dropped. The amount of temporary storage
    134 !>     needed is determined by the fftradix routine. Both fftn and fft
    135 !>     can handle any size of array. (Maybe they take a lot of time and
    136 !>     memory, but they will do it)
    137 !>
    138 !>   Redesigning fftradix in a way, that it handles assumed shape arrays
    139 !>   would have been desirable. However, I found it rather hard to do this
    140 !>   in an efficient way. Problems were:
    141 !>   - to prevent stride multiplications when indexing arrays. At least our
    142 !>     compiler was not clever enough to discover that in fact additions
    143 !>     would do the job as well. On the other hand, I haven't been clever
    144 !>     enough to find an implementation using array operations.
    145 !>   - fftradix is rather large and different versions would be necessaray
    146 !>     for each possible rank of array.
    147 !>   Consequently, in place transformation still needs the argument stored
    148 !>   in a consecutive bunch of memory and can't be performed on array
    149 !>   sections like A(100:199:-3, 50:1020). Calling fftn with such sections
    150 !>   will most probably imply copy-in-copy-out. However, the function fft
    151 !>   works with everything it gets and should be convenient to use.
     139!>   - Complex arguments and results are of type complex, rather than real an imaginary part
     140!>     separately.
     141!>   - Increment parameter (magnitude of isign) has been dropped, inc is always one, direction of
     142!>     transform is given by inv.
     143!>   - maxf and maxp have been dropped. The amount of temporary storage needed is determined by the
     144!>     fftradix routine. Both fftn and fft can handle any size of array. (Maybe they take a lot of
     145!>     time and memory, but they will do it)
     146!>
     147!>   Redesigning fftradix in a way, that it handles assumed shape arrays would have been desirable.
     148!>   However, I found it rather hard to do this in an efficient way. Problems were:
     149!>   - To prevent stride multiplications when indexing arrays. At least our compiler was not clever
     150!>     enough to discover that in fact additions would do the job as well. On the other hand, I
     151!>     haven't been clever enough to find an implementation using array operations.
     152!>   - fftradix is rather large and different versions would be necessaray for each possible rank of
     153!>     array.
     154!>   Consequently, in place transformation still needs the argument stored in a consecutive bunch of
     155!>   memory and can't be performed on array sections like A(100:199:-3, 50:1020). Calling fftn with
     156!>   such sections will most probably imply copy-in-copy-out. However, the function fft works with
     157!>   everything it gets and should be convenient to use.
    152158!>
    153159!> Michael Steffens, 09.12.96, <Michael.Steffens@mbox.muk.uni-hannover.de>
    154160!> Restructured fftradix for better optimization. M. Steffens, 4 June 1997
    155 !------------------------------------------------------------------------------!
     161!--------------------------------------------------------------------------------------------------!
    156162 MODULE singleton
    157  
     163
    158164
    159165    USE kinds
     
    162168
    163169    PRIVATE
    164     PUBLIC:: fft, fftn
    165 
    166     REAL(wp), PARAMETER:: sin60 = 0.86602540378443865_wp
    167     REAL(wp), PARAMETER:: cos72 = 0.30901699437494742_wp
    168     REAL(wp), PARAMETER:: sin72 = 0.95105651629515357_wp
    169     REAL(wp), PARAMETER:: pi    = 3.14159265358979323_wp
     170    PUBLIC ::  fft   !<
     171    PUBLIC ::  fftn  !<
     172
     173    REAL(wp), PARAMETER ::  cos72 = 0.30901699437494742_wp  !<
     174    REAL(wp), PARAMETER ::  pi    = 3.14159265358979323_wp  !<
     175    REAL(wp), PARAMETER ::  sin60 = 0.86602540378443865_wp  !<
     176    REAL(wp), PARAMETER ::  sin72 = 0.95105651629515357_wp  !<
    170177
    171178    INTERFACE fft
     
    183190
    184191
    185 !------------------------------------------------------------------------------!
     192!--------------------------------------------------------------------------------------------------!
    186193! Description:
    187194! ------------
    188195!> @todo Missing function description.
    189 !------------------------------------------------------------------------------!
    190  FUNCTION fft1d(array, dim, inv, stat) RESULT(ft)
     196!--------------------------------------------------------------------------------------------------!
     197 FUNCTION fft1d( array, dim, inv, stat ) RESULT( ft )
    191198!
    192199!-- Formal parameters
    193     COMPLEX(wp),  DIMENSION(:), INTENT(IN)           :: array
    194     INTEGER(iwp), DIMENSION(:), INTENT(IN),  OPTIONAL:: dim
    195     INTEGER(iwp),               INTENT(OUT), OPTIONAL:: stat
    196     LOGICAL,                    INTENT(IN),  OPTIONAL:: inv
     200    COMPLEX(wp),  DIMENSION(:), INTENT(IN)            ::  array  !<
     201    INTEGER(iwp), DIMENSION(:), INTENT(IN),  OPTIONAL ::  dim    !<
     202    INTEGER(iwp),               INTENT(OUT), OPTIONAL ::  stat   !<
     203    LOGICAL,                    INTENT(IN),  OPTIONAL ::  inv    !<
    197204!
    198205!-- Function result
    199     COMPLEX(wp), DIMENSION(SIZE(array, 1)):: ft
    200 
    201     INTEGER(iwp)::  idum
    202     INTEGER(iwp)::  ishape(1)
     206    COMPLEX(wp), DIMENSION(SIZE(array, 1)) ::  ft  !<
     207
     208    INTEGER(iwp) ::  idum       !<
     209    INTEGER(iwp) ::  ishape(1)  !<
    203210
    204211!
     
    208215    ft = array
    209216    ishape = SHAPE( array )
    210     CALL fftn(ft, ishape, inv = inv,  stat = stat)
     217    CALL fftn( ft, ishape, inv = inv,  stat = stat )
    211218!
    212219!-- Next statement to prevent compiler warning about unused variable
     
    216223
    217224
    218 !------------------------------------------------------------------------------!
     225!--------------------------------------------------------------------------------------------------!
    219226! Description:
    220227! ------------
    221228!> @todo Missing function description.
    222 !------------------------------------------------------------------------------!
    223  FUNCTION fft2d(array, dim, inv, stat) RESULT(ft)
     229!--------------------------------------------------------------------------------------------------!
     230 FUNCTION fft2d( array, dim, inv, stat ) RESULT( ft )
    224231!
    225232!-- Formal parameters
    226     COMPLEX(wp),  DIMENSION(:,:), INTENT(IN)           :: array
    227     INTEGER(iwp), DIMENSION(:),   INTENT(IN),  OPTIONAL:: dim
    228     INTEGER(iwp),                 INTENT(OUT), OPTIONAL:: stat
    229     LOGICAL,                      INTENT(IN),  OPTIONAL:: inv
     233    COMPLEX(wp),  DIMENSION(:,:), INTENT(IN)            ::  array  !<
     234    INTEGER(iwp), DIMENSION(:),   INTENT(IN),  OPTIONAL ::  dim    !<
     235    INTEGER(iwp),                 INTENT(OUT), OPTIONAL ::  stat   !<
     236    LOGICAL,                      INTENT(IN),  OPTIONAL ::  inv    !<
    230237!
    231238!-- Function result
    232     COMPLEX(wp), DIMENSION(SIZE(array, 1), SIZE(array, 2)):: ft
    233 
    234     INTEGER(iwp) ::  ishape(2)
     239    COMPLEX(wp), DIMENSION(SIZE(array, 1), SIZE(array, 2)) ::  ft  !<
     240
     241    INTEGER(iwp) ::  ishape(2)  !<
     242!
     243!-- Intrinsics used
     244    INTRINSIC SIZE, SHAPE
     245
     246    ft = array
     247    ishape = SHAPE( array )
     248    CALL fftn( ft, ishape, dim, inv, stat )
     249
     250 END FUNCTION fft2d
     251
     252
     253!--------------------------------------------------------------------------------------------------!
     254! Description:
     255! ------------
     256!> @todo Missing function description.
     257!--------------------------------------------------------------------------------------------------!
     258 FUNCTION fft3d( array, dim, inv, stat ) RESULT( ft )
     259!
     260!-- Formal parameters
     261    COMPLEX(wp),  DIMENSION(:,:,:), INTENT(IN)            ::  array  !<
     262    INTEGER(iwp), DIMENSION(:),     INTENT(IN),  OPTIONAL ::  dim    !<
     263    INTEGER(iwp),                   INTENT(OUT), OPTIONAL ::  stat   !<
     264    LOGICAL,                        INTENT(IN),  OPTIONAL ::  inv    !<
     265!
     266!-- Function result
     267    COMPLEX(wp), DIMENSION(SIZE(array, 1), SIZE(array, 2), SIZE(array, 3)) ::  ft  !<
     268
     269    INTEGER(iwp) ::  ishape(3)  !<
     270
     271!
     272!-- Intrinsics used
     273    INTRINSIC SIZE, SHAPE
     274
     275    ft = array
     276    ishape = SHAPE( array)
     277    CALL fftn(ft, ishape, dim, inv, stat)
     278
     279 END FUNCTION fft3d
     280
     281
     282!--------------------------------------------------------------------------------------------------!
     283! Description:
     284! ------------
     285!> @todo Missing function description.
     286!--------------------------------------------------------------------------------------------------!
     287 FUNCTION fft4d( array, dim, inv, stat ) RESULT( ft )
     288!
     289!-- Formal parameters
     290    COMPLEX(wp),  DIMENSION(:,:,:,:), INTENT(IN)            ::  array  !<
     291    INTEGER(iwp), DIMENSION(:),       INTENT(IN),  OPTIONAL ::  dim    !<
     292    INTEGER(iwp),                     INTENT(OUT), OPTIONAL ::  stat   !<
     293    LOGICAL,                          INTENT(IN),  OPTIONAL ::  inv    !<
     294!
     295!-- Function result
     296    COMPLEX(wp), DIMENSION(SIZE(array, 1), SIZE(array, 2), SIZE(array, 3), SIZE(array, 4)) ::  ft  !<
     297
     298    INTEGER(iwp) ::  ishape(4)  !<
    235299!
    236300!-- Intrinsics used
     
    241305    CALL fftn(ft, ishape, dim, inv, stat)
    242306
    243  END FUNCTION fft2d
    244 
    245 
    246 !------------------------------------------------------------------------------!
     307 END FUNCTION fft4d
     308
     309
     310!--------------------------------------------------------------------------------------------------!
    247311! Description:
    248312! ------------
    249313!> @todo Missing function description.
    250 !------------------------------------------------------------------------------!
    251  FUNCTION fft3d(array, dim, inv, stat) RESULT(ft)
     314!--------------------------------------------------------------------------------------------------!
     315 FUNCTION fft5d( array, dim, inv, stat ) RESULT( ft )
    252316!
    253317!-- Formal parameters
    254     COMPLEX(wp),  DIMENSION(:,:,:), INTENT(IN)           :: array
    255     INTEGER(iwp), DIMENSION(:),     INTENT(IN),  OPTIONAL:: dim
    256     INTEGER(iwp),                   INTENT(OUT), OPTIONAL:: stat
    257     LOGICAL,                        INTENT(IN),  OPTIONAL:: inv
     318    COMPLEX(wp),  DIMENSION(:,:,:,:,:), INTENT(IN)            ::  array  !<
     319    INTEGER(iwp), DIMENSION(:),         INTENT(IN),  OPTIONAL ::  dim    !<
     320    INTEGER(iwp),                       INTENT(OUT), OPTIONAL ::  stat   !<
     321    LOGICAL,                            INTENT(IN),  OPTIONAL ::  inv    !<
    258322!
    259323!-- Function result
    260     COMPLEX(wp), &
    261          DIMENSION(SIZE(array, 1), SIZE(array, 2), SIZE(array, 3)):: ft
    262 
    263     INTEGER(iwp) ::  ishape(3)
    264 
    265 !
    266 !-- Intrinsics used
    267     INTRINSIC SIZE, SHAPE
    268 
    269     ft = array
    270     ishape = SHAPE( array)
    271     CALL fftn(ft, ishape, dim, inv, stat)
    272 
    273  END FUNCTION fft3d
    274 
    275 
    276 !------------------------------------------------------------------------------!
    277 ! Description:
    278 ! ------------
    279 !> @todo Missing function description.
    280 !------------------------------------------------------------------------------!
    281  FUNCTION fft4d(array, dim, inv, stat) RESULT(ft)
    282 !
    283 !-- Formal parameters
    284     COMPLEX(wp),  DIMENSION(:,:,:,:), INTENT(IN)           :: array
    285     INTEGER(iwp), DIMENSION(:),       INTENT(IN),  OPTIONAL:: dim
    286     INTEGER(iwp),                     INTENT(OUT), OPTIONAL:: stat
    287     LOGICAL,                          INTENT(IN),  OPTIONAL:: inv
    288 !
    289 !-- Function result
    290     COMPLEX(wp), DIMENSION( &
    291          SIZE(array, 1), SIZE(array, 2), SIZE(array, 3), SIZE(array, 4)):: ft
    292 
    293     INTEGER(iwp) ::  ishape(4)
     324    COMPLEX(wp), DIMENSION(SIZE(array, 1), SIZE(array, 2), SIZE(array, 3), SIZE(array, 4), SIZE(array, 5)) ::  ft  !<
     325
     326    INTEGER(iwp) ::  ishape(5)  !<
     327
    294328!
    295329!-- Intrinsics used
     
    300334    CALL fftn(ft, ishape, dim, inv, stat)
    301335
    302  END FUNCTION fft4d
    303 
    304 
    305 !------------------------------------------------------------------------------!
     336 END FUNCTION fft5d
     337
     338
     339!--------------------------------------------------------------------------------------------------!
    306340! Description:
    307341! ------------
    308342!> @todo Missing function description.
    309 !------------------------------------------------------------------------------!
    310  FUNCTION fft5d(array, dim, inv, stat) RESULT(ft)
     343!--------------------------------------------------------------------------------------------------!
     344 FUNCTION fft6d( array, dim, inv, stat ) RESULT( ft )
    311345!
    312346!-- Formal parameters
    313     COMPLEX(wp),  DIMENSION(:,:,:,:,:), INTENT(IN)           :: array
    314     INTEGER(iwp), DIMENSION(:),         INTENT(IN),  OPTIONAL:: dim
    315     INTEGER(iwp),                       INTENT(OUT), OPTIONAL:: stat
    316     LOGICAL,                            INTENT(IN),  OPTIONAL:: inv
     347    COMPLEX(wp),  DIMENSION(:,:,:,:,:,:), INTENT(IN)            ::  array  !<
     348    INTEGER(iwp), DIMENSION(:),           INTENT(IN),  OPTIONAL ::  dim    !<
     349    INTEGER(iwp),                         INTENT(OUT), OPTIONAL ::  stat   !<
     350    LOGICAL,                              INTENT(IN),  OPTIONAL ::  inv    !<
    317351!
    318352!-- Function result
    319     COMPLEX(wp), DIMENSION( &
    320          SIZE(array, 1), SIZE(array, 2), SIZE(array, 3), SIZE(array, 4), &
    321          SIZE(array, 5)):: ft
    322 
    323     INTEGER(iwp) ::  ishape(5)
     353    COMPLEX(wp), DIMENSION( SIZE(array, 1), SIZE(array, 2), SIZE(array, 3), SIZE(array, 4),        &
     354                            SIZE(array, 5), SIZE(array, 6)) ::  ft  !<
     355
     356    INTEGER(iwp) ::  ishape(6)  !<
    324357
    325358!
     
    331364    CALL fftn(ft, ishape, dim, inv, stat)
    332365
    333  END FUNCTION fft5d
    334 
    335 
    336 !------------------------------------------------------------------------------!
     366 END FUNCTION fft6d
     367
     368
     369!--------------------------------------------------------------------------------------------------!
    337370! Description:
    338371! ------------
    339372!> @todo Missing function description.
    340 !------------------------------------------------------------------------------!
    341  FUNCTION fft6d(array, dim, inv, stat) RESULT(ft)
     373!--------------------------------------------------------------------------------------------------!
     374 FUNCTION fft7d( array, dim, inv, stat ) RESULT( ft )
    342375!
    343376!-- Formal parameters
    344     COMPLEX(wp),  DIMENSION(:,:,:,:,:,:), INTENT(IN)           :: array
    345     INTEGER(iwp), DIMENSION(:),           INTENT(IN),  OPTIONAL:: dim
    346     INTEGER(iwp),                         INTENT(OUT), OPTIONAL:: stat
    347     LOGICAL,                              INTENT(IN),  OPTIONAL:: inv
     377    COMPLEX(wp), DIMENSION(:,:,:,:,:,:,:), INTENT(IN)            ::  array  !<
     378    INTEGER(iwp),            DIMENSION(:), INTENT(IN),  OPTIONAL ::  dim    !<
     379    INTEGER(iwp),                          INTENT(OUT), OPTIONAL ::  stat   !<
     380    LOGICAL,                               INTENT(IN),  OPTIONAL ::  inv    !<
    348381!
    349382!-- Function result
    350     COMPLEX(wp), DIMENSION( &
    351          SIZE(array, 1), SIZE(array, 2), SIZE(array, 3), SIZE(array, 4), &
    352          SIZE(array, 5), SIZE(array, 6)):: ft
    353 
    354     INTEGER(iwp) ::  ishape(6)
     383    COMPLEX(wp), DIMENSION( SIZE(array, 1), SIZE(array, 2), SIZE(array, 3), SIZE(array, 4),        &
     384                            SIZE(array, 5), SIZE(array, 6), SIZE(array, 7)) ::  ft  !<
     385
     386    INTEGER(iwp) ::  ishape(7)  !<
    355387
    356388!
     
    362394    CALL fftn(ft, ishape, dim, inv, stat)
    363395
    364  END FUNCTION fft6d
    365 
    366 
    367 !------------------------------------------------------------------------------!
    368 ! Description:
    369 ! ------------
    370 !> @todo Missing function description.
    371 !------------------------------------------------------------------------------!
    372  FUNCTION fft7d(array, dim, inv, stat) RESULT(ft)
     396 END FUNCTION fft7d
     397
     398
     399!--------------------------------------------------------------------------------------------------!
     400! Description:
     401! ------------
     402!> @todo Missing subroutine description.
     403!--------------------------------------------------------------------------------------------------!
     404 SUBROUTINE fftn( array, shape, dim, inv, stat )
    373405!
    374406!-- Formal parameters
    375     COMPLEX(wp), DIMENSION(:,:,:,:,:,:,:), INTENT(IN)           :: array
    376     INTEGER(iwp),          DIMENSION(:),   INTENT(IN),  OPTIONAL:: dim
    377     INTEGER(iwp),                          INTENT(OUT), OPTIONAL:: stat
    378     LOGICAL,                               INTENT(IN),  OPTIONAL:: inv
    379 !
    380 !-- Function result
    381     COMPLEX(wp), DIMENSION( &
    382          SIZE(array, 1), SIZE(array, 2), SIZE(array, 3), SIZE(array, 4), &
    383          SIZE(array, 5), SIZE(array, 6), SIZE(array, 7)):: ft
    384 
    385     INTEGER(iwp) ::  ishape(7)
    386 
    387 !
    388 !-- Intrinsics used
    389     INTRINSIC SIZE, SHAPE
    390 
    391     ft = array
    392     ishape = SHAPE( array )
    393     CALL fftn(ft, ishape, dim, inv, stat)
    394 
    395  END FUNCTION fft7d
    396 
    397 
    398 !------------------------------------------------------------------------------!
    399 ! Description:
    400 ! ------------
    401 !> @todo Missing subroutine description.
    402 !------------------------------------------------------------------------------!
    403  SUBROUTINE fftn(array, shape, dim, inv, stat)
    404 !
    405 !-- Formal parameters
    406     COMPLEX(wp),  DIMENSION(*), INTENT(INOUT)        :: array
    407     INTEGER(iwp), DIMENSION(:), INTENT(IN)           :: shape
    408     INTEGER(iwp), DIMENSION(:), INTENT(IN),  OPTIONAL:: dim
    409     INTEGER(iwp),               INTENT(OUT), OPTIONAL:: stat
    410     LOGICAL,                    INTENT(IN),  OPTIONAL:: inv
     407    COMPLEX(wp),  DIMENSION(*), INTENT(INOUT)         ::  array  !<
     408    INTEGER(iwp), DIMENSION(:), INTENT(IN)            ::  shape  !<
     409    INTEGER(iwp), DIMENSION(:), INTENT(IN),  OPTIONAL ::  dim    !<
     410    INTEGER(iwp),               INTENT(OUT), OPTIONAL ::  stat   !<
     411    LOGICAL,                    INTENT(IN),  OPTIONAL ::  inv    !<
    411412!
    412413!-- Local arrays
    413     INTEGER(iwp), DIMENSION(SIZE(shape)):: d
     414    INTEGER(iwp), DIMENSION(SIZE(shape)) ::  d  !<
    414415!
    415416!-- Local scalars
    416     LOGICAL      :: inverse
    417     INTEGER(iwp) :: i, ndim, ntotal
    418     REAL(wp):: scale
     417    LOGICAL      ::  inverse          !<
     418    INTEGER(iwp) ::  i, ndim, ntotal  !<
     419    REAL(wp)     ::  scale            !<
    419420!
    420421!-- Intrinsics used
     
    423424!
    424425!-- Optional parameter settings
    425     IF (PRESENT(inv)) THEN
     426    IF ( PRESENT( inv ) ) THEN
    426427       inverse = inv
    427428    ELSE
    428429       inverse = .FALSE.
    429430    END IF
    430     IF (PRESENT(dim)) THEN
    431        ndim = MIN(SIZE(dim), SIZE(d))
    432        d(1:ndim) = DIM(1:ndim)
     431    IF ( PRESENT( dim ) ) THEN
     432       ndim = MIN( SIZE( dim ), SIZE( d ) )
     433       d(1:ndim) = DIM( 1:ndim )
    433434    ELSE
    434        ndim = SIZE(d)
    435        d = (/(i, i = 1, SIZE(d))/)
     435       ndim = SIZE( d )
     436        d = (/( i, i = 1, SIZE( d ) )/)
    436437    END IF
    437438
    438     ntotal = PRODUCT(shape)
    439     scale = SQRT(1.0_wp / PRODUCT(shape(d(1:ndim))))
    440     DO i = 1, ntotal
    441        array(i) = CMPLX(REAL(array(i)) * scale, AIMAG(array(i)) * scale, &
    442             KIND=wp)
     439    ntotal = PRODUCT( shape )
     440    scale = SQRT( 1.0_wp / PRODUCT( shape( d(1:ndim) ) ) )
     441    DO  i = 1, ntotal
     442       array(i) = CMPLX( REAL( array(i) ) * scale, AIMAG( array(i) ) * scale, KIND = wp )
    443443    END DO
    444444
    445     DO i = 1, ndim
    446        CALL fftradix(array, ntotal, shape(d(i)), PRODUCT(shape(1:d(i))), &
    447             inverse, stat)
    448        IF (PRESENT(stat)) THEN
    449           IF (stat /=0) RETURN
     445    DO  i = 1, ndim
     446       CALL fftradix( array, ntotal, shape( d(i) ), PRODUCT( shape( 1:d(i) ) ), inverse, stat )
     447       IF ( PRESENT( stat ) )  THEN
     448          IF ( stat /= 0 )  RETURN
    450449       END IF
    451450    END DO
     
    454453
    455454
    456 !------------------------------------------------------------------------------!
     455!--------------------------------------------------------------------------------------------------!
    457456! Description:
    458457! ------------
    459458!> @todo Missing subroutine description.
    460 !------------------------------------------------------------------------------!
    461  SUBROUTINE fftradix(array, ntotal, npass, nspan, inv, stat)
     459!--------------------------------------------------------------------------------------------------!
     460 SUBROUTINE fftradix( array, ntotal, npass, nspan, inv, stat )
    462461!
    463462!-- Formal parameters
    464     COMPLEX(wp),  DIMENSION(*), INTENT(INOUT)        :: array
    465     INTEGER(iwp),               INTENT(IN)           :: ntotal, npass, nspan
    466     INTEGER(iwp),               INTENT(OUT), OPTIONAL:: stat
    467     LOGICAL,                    INTENT(IN)           :: inv
     463    COMPLEX(wp),  DIMENSION(*), INTENT(INOUT)         ::  array                 !<
     464    INTEGER(iwp),               INTENT(IN)            ::  ntotal, npass, nspan  !<
     465    INTEGER(iwp),               INTENT(OUT), OPTIONAL ::  stat                  !<
     466    LOGICAL,                    INTENT(IN)            ::  inv                   !<
    468467!
    469468!-- Local arrays
    470     COMPLEX(wp),  DIMENSION(:), ALLOCATABLE  :: ctmp
    471     INTEGER(iwp), DIMENSION(BIT_SIZE(0))     :: factor
    472     INTEGER(iwp), DIMENSION(:), ALLOCATABLE  :: perm
    473     REAL(wp),     DIMENSION(:), ALLOCATABLE  :: sine, cosine
     469    COMPLEX(wp),  DIMENSION(:), ALLOCATABLE ::  ctmp          !<
     470    INTEGER(iwp), DIMENSION(BIT_SIZE(0))    ::  factor        !<
     471    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  perm          !<
     472    REAL(wp),     DIMENSION(:), ALLOCATABLE ::  sine, cosine  !<
    474473!
    475474!-- Local scalars
    476     INTEGER(iwp)         :: maxfactor, nfactor, nsquare, nperm
     475    INTEGER(iwp) ::  maxfactor, nfactor, nsquare, nperm  !<
    477476!
    478477!-- Intrinsics used
    479     INTRINSIC MAXVAL, MOD, PRESENT, ISHFT, BIT_SIZE, SIN, COS, &
    480          CMPLX, REAL, AIMAG
    481 
    482     IF (npass <= 1) RETURN
    483 
    484     CALL factorize(npass, factor, nfactor, nsquare)
    485 
    486     maxfactor = MAXVAL(factor(:nfactor))
    487     IF (nfactor - ISHFT(nsquare, 1) > 0) THEN
    488        nperm = MAX(nfactor + 1, PRODUCT(factor(nsquare+1: nfactor-nsquare)) - 1)
     478    INTRINSIC MAXVAL, MOD, PRESENT, ISHFT, BIT_SIZE, SIN, COS, CMPLX, REAL, AIMAG
     479
     480    IF ( npass <= 1 )  RETURN
     481
     482    CALL factorize( npass, factor, nfactor, nsquare )
     483
     484    maxfactor = MAXVAL( factor(:nfactor) )
     485    IF ( nfactor - ISHFT( nsquare, 1 ) > 0 )  THEN
     486       nperm = MAX( nfactor + 1, PRODUCT( factor(nsquare+1: nfactor-nsquare) ) - 1 )
    489487    ELSE
    490488       nperm = nfactor + 1
    491489    END IF
    492490
    493     IF (PRESENT(stat)) THEN
    494        ALLOCATE(ctmp(maxfactor), sine(maxfactor), cosine(maxfactor), STAT=stat)
    495        IF (stat /= 0) RETURN
    496        CALL transform(array, ntotal, npass, nspan, &
    497             factor, nfactor, ctmp, sine, cosine, inv)
    498        DEALLOCATE(sine, cosine, STAT=stat)
    499        IF (stat /= 0) RETURN
    500        ALLOCATE(perm(nperm), STAT=stat)
    501        IF (stat /= 0) RETURN
    502        CALL permute(array, ntotal, npass, nspan, &
    503             factor, nfactor, nsquare, maxfactor, &
    504             ctmp, perm)
    505        DEALLOCATE(perm, ctmp, STAT=stat)
    506        IF (stat /= 0) RETURN
     491    IF ( PRESENT( stat ) )  THEN
     492       ALLOCATE( ctmp(maxfactor), sine(maxfactor), cosine(maxfactor), STAT = stat )
     493       IF ( stat /= 0 )  RETURN
     494       CALL transform( array, ntotal, npass, nspan, factor, nfactor, ctmp, sine, cosine, inv )
     495       DEALLOCATE( sine, cosine, STAT = stat )
     496       IF ( stat /= 0 )  RETURN
     497       ALLOCATE( perm(nperm), STAT = stat )
     498       IF ( stat /= 0 )  RETURN
     499       CALL permute( array, ntotal, npass, nspan, factor, nfactor, nsquare, maxfactor, ctmp, perm )
     500       DEALLOCATE( perm, ctmp, STAT = stat )
     501       IF ( stat /= 0 )  RETURN
    507502    ELSE
    508        ALLOCATE(ctmp(maxfactor), sine(maxfactor), cosine(maxfactor))
    509        CALL transform(array, ntotal, npass, nspan, &
    510             factor, nfactor, ctmp, sine, cosine, inv)
    511        DEALLOCATE(sine, cosine)
    512        ALLOCATE(perm(nperm))
    513        CALL permute(array, ntotal, npass, nspan, &
    514             factor, nfactor, nsquare, maxfactor, &
    515             ctmp, perm)
    516        DEALLOCATE(perm, ctmp)
     503       ALLOCATE( ctmp(maxfactor), sine(maxfactor), cosine(maxfactor) )
     504       CALL transform( array, ntotal, npass, nspan, factor, nfactor, ctmp, sine, cosine, inv )
     505       DEALLOCATE( sine, cosine )
     506       ALLOCATE( perm(nperm) )
     507       CALL permute( array, ntotal, npass, nspan, factor, nfactor, nsquare, maxfactor, ctmp, perm )
     508       DEALLOCATE( perm, ctmp )
    517509    END IF
    518510
    519511
    520   CONTAINS
    521 
    522 
    523 !------------------------------------------------------------------------------!
     512 CONTAINS
     513
     514
     515!--------------------------------------------------------------------------------------------------!
    524516! Description:
    525517! ------------
    526518!> @todo Missing subroutine description.
    527 !------------------------------------------------------------------------------!
    528     SUBROUTINE factorize(npass, factor, nfactor, nsquare)
    529 !
    530 !--   Formal parameters
    531       INTEGER(iwp),               INTENT(IN) :: npass
    532       INTEGER(iwp), DIMENSION(*), INTENT(OUT):: factor
    533       INTEGER(iwp),               INTENT(OUT):: nfactor, nsquare
    534 !
    535 !--   Local scalars
    536       INTEGER(iwp):: j, jj, k
    537 
    538       nfactor = 0
    539       k = npass
    540       DO WHILE (MOD(k, 16) == 0)
    541          nfactor = nfactor + 1
    542          factor(nfactor) = 4
    543          k = k / 16
    544       END DO
    545       j = 3
    546       jj = 9
    547       DO
    548          DO WHILE (MOD(k, jj) == 0)
    549             nfactor = nfactor + 1
    550             factor(nfactor) = j
    551             k = k / jj
    552          END DO
    553          j = j + 2
    554          jj = j * j
    555          IF (jj > k) EXIT
    556       END DO
    557       IF (k <= 4) THEN
    558          nsquare = nfactor
    559          factor(nfactor + 1) = k
    560          IF (k /= 1) nfactor = nfactor + 1
    561       ELSE
    562          IF (k - ISHFT(k / 4, 2) == 0) THEN
    563             nfactor = nfactor + 1
    564             factor(nfactor) = 2
    565             k = k / 4
    566          END IF
    567          nsquare = nfactor
    568          j = 2
    569          DO
    570             IF (MOD(k, j) == 0) THEN
    571                nfactor = nfactor + 1
    572                factor(nfactor) = j
    573                k = k / j
    574             END IF
    575             j = ISHFT((j + 1) / 2, 1) + 1
    576             IF (j > k) EXIT
    577          END DO
    578       END IF
    579       IF (nsquare > 0) THEN
    580          j = nsquare
    581          DO
    582             nfactor = nfactor + 1
    583             factor(nfactor) = factor(j)
    584             j = j - 1
    585             IF (j==0) EXIT
    586          END DO
    587       END IF
    588 
    589     END SUBROUTINE factorize
    590 
    591 
    592 !------------------------------------------------------------------------------!
     519!--------------------------------------------------------------------------------------------------!
     520 SUBROUTINE factorize( npass, factor, nfactor, nsquare )
     521!
     522!-- Formal parameters
     523    INTEGER(iwp),               INTENT(IN)  ::  npass             !<
     524    INTEGER(iwp),               INTENT(OUT) ::  nfactor, nsquare  !<
     525    INTEGER(iwp), DIMENSION(*), INTENT(OUT) ::  factor            !<
     526!
     527!-- Local scalars
     528    INTEGER(iwp) ::  j, jj, k  !<
     529
     530    nfactor = 0
     531    k = npass
     532    DO WHILE ( MOD( k, 16 ) == 0 )
     533       nfactor = nfactor + 1
     534       factor(nfactor) = 4
     535       k = k / 16
     536    END DO
     537    j = 3
     538    jj = 9
     539    DO
     540       DO WHILE ( MOD( k, jj ) == 0 )
     541          nfactor = nfactor + 1
     542          factor(nfactor) = j
     543          k = k / jj
     544       END DO
     545       j = j + 2
     546       jj = j * j
     547       IF ( jj > k ) EXIT
     548    END DO
     549    IF ( k <= 4 ) THEN
     550       nsquare = nfactor
     551       factor(nfactor + 1) = k
     552       IF ( k /= 1 ) nfactor = nfactor + 1
     553    ELSE
     554       IF ( k - ISHFT( k / 4, 2 ) == 0 ) THEN
     555          nfactor = nfactor + 1
     556          factor(nfactor) = 2
     557          k = k / 4
     558       END IF
     559       nsquare = nfactor
     560       j = 2
     561       DO
     562          IF ( MOD(k, j) == 0 ) THEN
     563             nfactor = nfactor + 1
     564             factor(nfactor) = j
     565             k = k / j
     566          END IF
     567          j = ISHFT( (j + 1) / 2, 1 ) + 1
     568          IF ( j > k ) EXIT
     569       END DO
     570    END IF
     571    IF ( nsquare > 0 ) THEN
     572       j = nsquare
     573       DO
     574          nfactor = nfactor + 1
     575          factor(nfactor) = factor(j)
     576          j = j - 1
     577          IF ( j == 0 ) EXIT
     578       END DO
     579    END IF
     580
     581 END SUBROUTINE factorize
     582
     583
     584!--------------------------------------------------------------------------------------------------!
    593585! Description:
    594586! ------------
    595587!> @todo Missing subroutine description.
    596 !------------------------------------------------------------------------------!
    597     SUBROUTINE transform(array, ntotal, npass, nspan, &
    598          factor, nfactor, ctmp, sine, cosine, inv) !-- compute fourier transform
    599 !
    600 !--   Formal parameters
    601       COMPLEX(wp),  DIMENSION(*), INTENT(IN OUT):: array
    602       COMPLEX(wp),  DIMENSION(*), INTENT(OUT)   :: ctmp
    603       INTEGER(iwp),               INTENT(IN)    :: ntotal, npass, nspan
    604       INTEGER(iwp), DIMENSION(*), INTENT(IN)    :: factor
    605       INTEGER(iwp),               INTENT(IN)    :: nfactor
    606       LOGICAL,                    INTENT(IN)    :: inv
    607       REAL(wp),     DIMENSION(*), INTENT(OUT)   :: sine, cosine
    608 !
    609 !--   Local scalars
    610       INTEGER(iwp):: ii, ispan
    611       INTEGER(iwp):: j, jc, jf, jj
    612       INTEGER(iwp):: k, kk, kspan, k1, k2, k3, k4
    613       INTEGER(iwp):: nn, nt
    614       REAL(wp)    :: s60, c72, s72, pi2, radf
    615       REAL(wp)    :: c1, s1, c2, s2, c3, s3, cd, sd, ak
    616       COMPLEX(wp) :: cc, cj, ck, cjp, cjm, ckp, ckm
    617 
    618       c72 = cos72
    619       IF (inv) THEN
    620          s72 = sin72
    621          s60 = sin60
    622          pi2 = pi
    623       ELSE
    624          s72 = -sin72
    625          s60 = -sin60
    626          pi2 = -pi
    627       END IF
    628 
    629       nt = ntotal
    630       nn = nt - 1
    631       kspan = nspan
    632       jc = nspan / npass
    633       radf = pi2 * jc
    634       pi2 = pi2 * 2.0_wp !-- use 2 PI from here on
    635 
    636       ii = 0
    637       jf = 0
    638       DO
    639          sd = radf / kspan
    640          cd = SIN(sd)
    641          cd = 2.0_wp * cd * cd
    642          sd = SIN(sd + sd)
    643          kk = 1
    644          ii = ii + 1
    645 
    646          SELECT CASE (factor(ii))
    647          CASE (2)
    648 !
    649 !--         Transform for factor of 2 (including rotation factor)
    650             kspan = kspan / 2
    651             k1 = kspan + 2
    652             DO
    653                DO
    654                   k2 = kk + kspan
    655                   ck = array(k2)
    656                   array(k2) = array(kk)-ck
    657                   array(kk) = array(kk) + ck
    658                   kk = k2 + kspan
    659                   IF (kk > nn) EXIT
    660                END DO
    661                kk = kk - nn
    662                IF (kk > jc) EXIT
    663             END DO
    664             IF (kk > kspan) RETURN
    665             DO
    666                c1 = 1.0_wp - cd
    667                s1 = sd
    668                DO
    669                   DO
    670                      DO
    671                         k2 = kk + kspan
    672                         ck = array(kk) - array(k2)
    673                         array(kk) = array(kk) + array(k2)
    674                         array(k2) = ck * CMPLX(c1, s1, KIND=wp)
    675                         kk = k2 + kspan
    676                         IF (kk >= nt) EXIT
    677                      END DO
    678                      k2 = kk - nt
    679                      c1 = -c1
    680                      kk = k1 - k2
    681                      IF (kk <= k2) EXIT
    682                   END DO
    683                   ak = c1 - (cd * c1 + sd * s1)
    684                   s1 = sd * c1 - cd * s1 + s1
    685                   c1 = 2.0_wp - (ak * ak + s1 * s1)
    686                   s1 = s1 * c1
    687                   c1 = c1 * ak
    688                   kk = kk + jc
    689                   IF (kk >= k2) EXIT
    690                END DO
    691                k1 = k1 + 1 + 1
    692                kk = (k1 - kspan) / 2 + jc
    693                IF (kk > jc + jc) EXIT
    694             END DO
    695 
    696          CASE (4) !-- transform for factor of 4
    697             ispan = kspan
    698             kspan = kspan / 4
    699 
    700             DO
    701                c1 = 1.0_wp
    702                s1 = 0.0_wp
    703                DO
    704                   DO
    705                      k1 = kk + kspan
    706                      k2 = k1 + kspan
    707                      k3 = k2 + kspan
    708                      ckp = array(kk) + array(k2)
    709                      ckm = array(kk) - array(k2)
    710                      cjp = array(k1) + array(k3)
    711                      cjm = array(k1) - array(k3)
    712                      array(kk) = ckp + cjp
    713                      cjp = ckp - cjp
    714                      IF (inv) THEN
    715                         ckp = ckm + CMPLX(-AIMAG(cjm), REAL(cjm), KIND=wp)
    716                         ckm = ckm + CMPLX(AIMAG(cjm), -REAL(cjm), KIND=wp)
    717                      ELSE
    718                         ckp = ckm + CMPLX(AIMAG(cjm), -REAL(cjm), KIND=wp)
    719                         ckm = ckm + CMPLX(-AIMAG(cjm), REAL(cjm), KIND=wp)
    720                      END IF
    721 !
    722 !--                  Avoid useless multiplies
    723                      IF (s1 == 0.0_wp) THEN
    724                         array(k1) = ckp
    725                         array(k2) = cjp
    726                         array(k3) = ckm
    727                      ELSE
    728                         array(k1) = ckp * CMPLX(c1, s1, KIND=wp)
    729                         array(k2) = cjp * CMPLX(c2, s2, KIND=wp)
    730                         array(k3) = ckm * CMPLX(c3, s3, KIND=wp)
    731                      END IF
    732                      kk = k3 + kspan
    733                      IF (kk > nt) EXIT
    734                   END DO
    735 
    736                   c2 = c1 - (cd * c1 + sd * s1)
    737                   s1 = sd * c1 - cd * s1 + s1
    738                   c1 = 2.0_wp - (c2 * c2 + s1 * s1)
    739                   s1 = s1 * c1
    740                   c1 = c1 * c2
    741 !
    742 !--               Values of c2, c3, s2, s3 that will get used next time
    743                   c2 = c1 * c1 - s1 * s1
    744                   s2 = 2.0_wp * c1 * s1
    745                   c3 = c2 * c1 - s2 * s1
    746                   s3 = c2 * s1 + s2 * c1
    747                   kk = kk - nt + jc
    748                   IF (kk > kspan) EXIT
    749                END DO
    750                kk = kk - kspan + 1
    751                IF (kk > jc) EXIT
    752             END DO
    753             IF (kspan == jc) RETURN
    754 
    755          CASE default
    756 !
    757 !--         Transform for odd factors
    758             k = factor(ii)
    759             ispan = kspan
    760             kspan = kspan / k
    761 
    762             SELECT CASE (k)
    763             CASE (3) !-- transform for factor of 3 (optional code)
    764                DO
    765                   DO
    766                      k1 = kk + kspan
    767                      k2 = k1 + kspan
    768                      ck = array(kk)
    769                      cj = array(k1) + array(k2)
    770                      array(kk) = ck + cj
    771                      ck = ck - CMPLX( &
    772                           0.5_wp * REAL (cj), &
    773                           0.5_wp * AIMAG(cj), &
    774                           KIND=wp)
    775                      cj = CMPLX( &
    776                           (REAL (array(k1)) - REAL (array(k2))) * s60, &
    777                           (AIMAG(array(k1)) - AIMAG(array(k2))) * s60, &
    778                           KIND=wp)
    779                      array(k1) = ck + CMPLX(-AIMAG(cj), REAL(cj), KIND=wp)
    780                      array(k2) = ck + CMPLX(AIMAG(cj), -REAL(cj), KIND=wp)
    781                      kk = k2 + kspan
    782                      IF (kk >= nn) EXIT
    783                   END DO
    784                   kk = kk - nn
    785                   IF (kk > kspan) EXIT
    786                END DO
    787 
    788             CASE (5) !-- transform for factor of 5 (optional code)
    789                c2 = c72 * c72 - s72 * s72
    790                s2 = 2.0_wp * c72 * s72
    791                DO
    792                   DO
    793                      k1 = kk + kspan
    794                      k2 = k1 + kspan
    795                      k3 = k2 + kspan
    796                      k4 = k3 + kspan
    797                      ckp = array(k1) + array(k4)
    798                      ckm = array(k1) - array(k4)
    799                      cjp = array(k2) + array(k3)
    800                      cjm = array(k2) - array(k3)
    801                      cc = array(kk)
    802                      array(kk) = cc + ckp + cjp
    803                      ck = CMPLX(REAL(ckp) * c72, AIMAG(ckp) * c72, &
    804                           KIND=wp) + &
    805                           CMPLX(REAL(cjp) * c2,  AIMAG(cjp) * c2,  &
    806                           KIND=wp) + cc
    807                      cj = CMPLX(REAL(ckm) * s72, AIMAG(ckm) * s72, &
    808                           KIND=wp) + &
    809                           CMPLX(REAL(cjm) * s2,  AIMAG(cjm) * s2,  &
    810                           KIND=wp)
    811                      array(k1) = ck + CMPLX(-AIMAG(cj), REAL(cj), KIND=wp)
    812                      array(k4) = ck + CMPLX(AIMAG(cj), -REAL(cj), KIND=wp)
    813                      ck = CMPLX(REAL(ckp) * c2,  AIMAG(ckp) * c2,  &
    814                           KIND=wp) + &
    815                           CMPLX(REAL(cjp) * c72, AIMAG(cjp) * c72, &
    816                           KIND=wp) + cc
    817                      cj = CMPLX(REAL(ckm) * s2,  AIMAG(ckm) * s2,  &
    818                           KIND=wp) - &
    819                           CMPLX(REAL(cjm) * s72, AIMAG(cjm) * s72, &
    820                           KIND=wp)
    821                      array(k2) = ck + CMPLX(-AIMAG(cj), REAL(cj), KIND=wp)
    822                      array(k3) = ck + CMPLX(AIMAG(cj), -REAL(cj), KIND=wp)
    823                      kk = k4 + kspan
    824                      IF (kk >= nn) EXIT
    825                   END DO
    826                   kk = kk - nn
    827                   IF (kk > kspan) EXIT
    828                END DO
    829 
    830             CASE default
    831                IF (k /= jf) THEN
    832                   jf = k
    833                   s1 = pi2 / k
    834                   c1 = COS(s1)
    835                   s1 = SIN(s1)
    836                   cosine (jf) = 1.0_wp
    837                   sine (jf) = 0.0_wp
    838                   j = 1
    839                   DO
    840                      cosine (j) = cosine (k) * c1 + sine (k) * s1
    841                      sine (j) = cosine (k) * s1 - sine (k) * c1
    842                      k = k-1
    843                      cosine (k) = cosine (j)
    844                      sine (k) = -sine (j)
    845                      j = j + 1
    846                      IF (j >= k) EXIT
    847                   END DO
    848                END IF
    849                DO
    850                   DO
    851                      k1 = kk
    852                      k2 = kk + ispan
    853                      cc = array(kk)
    854                      ck = cc
    855                      j = 1
    856                      k1 = k1 + kspan
    857                      DO
    858                         k2 = k2 - kspan
    859                         j = j + 1
    860                         ctmp(j) = array(k1) + array(k2)
    861                         ck = ck + ctmp(j)
    862                         j = j + 1
    863                         ctmp(j) = array(k1) - array(k2)
    864                         k1 = k1 + kspan
    865                         IF (k1 >= k2) EXIT
    866                      END DO
    867                      array(kk) = ck
    868                      k1 = kk
    869                      k2 = kk + ispan
    870                      j = 1
    871                      DO
    872                         k1 = k1 + kspan
    873                         k2 = k2 - kspan
    874                         jj = j
    875                         ck = cc
    876                         cj = (0.0_wp, 0.0_wp)
    877                         k = 1
    878                         DO
    879                            k = k + 1
    880                            ck = ck + CMPLX( &
    881                                 REAL (ctmp(k)) * cosine(jj), &
    882                                 AIMAG(ctmp(k)) * cosine(jj), KIND=wp)
    883                            k = k + 1
    884                            cj = cj + CMPLX( &
    885                                 REAL (ctmp(k)) * sine(jj), &
    886                                 AIMAG(ctmp(k)) * sine(jj), KIND=wp)
    887                            jj = jj + j
    888                            IF (jj > jf) jj = jj - jf
    889                            IF (k >= jf) EXIT
    890                         END DO
    891                         k = jf - j
    892                         array(k1) = ck + CMPLX(-AIMAG(cj), REAL(cj), &
    893                              KIND=wp)
    894                         array(k2) = ck + CMPLX(AIMAG(cj), -REAL(cj), &
    895                              KIND=wp)
    896                         j = j + 1
    897                         IF (j >= k) EXIT
    898                      END DO
    899                      kk = kk + ispan
    900                      IF (kk > nn) EXIT
    901                   END DO
    902                   kk = kk - nn
    903                   IF (kk > kspan) EXIT
    904                END DO
    905 
    906             END SELECT
    907 !
    908 !--         Multiply by rotation factor (except for factors of 2 and 4)
    909             IF (ii == nfactor) RETURN
    910             kk = jc + 1
    911             DO
    912                c2 = 1.0_wp - cd
    913                s1 = sd
    914                DO
    915                   c1 = c2
    916                   s2 = s1
    917                   kk = kk + kspan
    918                   DO
    919                      DO
    920                         array(kk) = CMPLX(c2, s2, KIND=wp) * array(kk)
    921                         kk = kk + ispan
    922                         IF (kk > nt) EXIT
    923                      END DO
    924                      ak = s1 * s2
    925                      s2 = s1 * c2 + c1 * s2
    926                      c2 = c1 * c2 - ak
    927                      kk = kk - nt + kspan
    928                      IF (kk > ispan) EXIT
    929                   END DO
    930                   c2 = c1 - (cd * c1 + sd * s1)
    931                   s1 = s1 + sd * c1 - cd * s1
    932                   c1 = 2.0_wp - (c2 * c2 + s1 * s1)
    933                   s1 = s1 * c1
    934                   c2 = c2 * c1
    935                   kk = kk - ispan + jc
    936                   IF (kk > kspan) EXIT
    937                END DO
    938                kk = kk - kspan + jc + 1
    939                IF (kk > jc + jc) EXIT
    940             END DO
    941 
    942          END SELECT
    943       END DO
    944     END SUBROUTINE transform
    945 
    946 
    947 !------------------------------------------------------------------------------!
     588!--------------------------------------------------------------------------------------------------!
     589 SUBROUTINE transform( array, ntotal, npass, nspan, factor, nfactor, ctmp, sine, cosine, inv )
     590!-- Compute fourier transform
     591
     592!
     593!-- Formal parameters
     594    COMPLEX(wp),  DIMENSION(*), INTENT(IN OUT) ::  array                 !<
     595    COMPLEX(wp),  DIMENSION(*), INTENT(OUT)    ::  ctmp                  !<
     596    INTEGER(iwp),               INTENT(IN)     ::  ntotal, npass, nspan  !<
     597    INTEGER(iwp), DIMENSION(*), INTENT(IN)     ::  factor                !<
     598    INTEGER(iwp),               INTENT(IN)     ::  nfactor               !<
     599    LOGICAL,                    INTENT(IN)     ::  inv                   !<
     600    REAL(wp),     DIMENSION(*), INTENT(OUT)    ::  sine, cosine          !<
     601!
     602!-- Local scalars
     603    COMPLEX(wp)  ::  cc, cj, ck, cjp, cjm, ckp, ckm      !<
     604    INTEGER(iwp) ::  ii, ispan                           !<
     605    INTEGER(iwp) ::  j, jc, jf, jj                       !<
     606    INTEGER(iwp) ::  k, kk, kspan, k1, k2, k3, k4        !<
     607    INTEGER(iwp) ::  nn, nt                              !<
     608    REAL(wp)     ::  s60, c72, s72, pi2, radf            !<
     609    REAL(wp)     ::  c1, s1, c2, s2, c3, s3, cd, sd, ak  !<
     610
     611    c72 = cos72
     612    IF ( inv )  THEN
     613       s72 = sin72
     614       s60 = sin60
     615       pi2 = pi
     616    ELSE
     617       s72 = - sin72
     618       s60 = - sin60
     619       pi2 = - pi
     620    END IF
     621
     622    nt = ntotal
     623    nn = nt - 1
     624    kspan = nspan
     625    jc = nspan / npass
     626    radf = pi2 * jc
     627    pi2 = pi2 * 2.0_wp !-- Use 2 PI from here on
     628
     629    ii = 0
     630    jf = 0
     631    DO
     632       sd = radf / kspan
     633       cd = SIN( sd )
     634       cd = 2.0_wp * cd * cd
     635       sd = SIN( sd + sd )
     636       kk = 1
     637       ii = ii + 1
     638
     639       SELECT CASE ( factor(ii) )
     640       CASE ( 2 )
     641!
     642!--       Transform for factor of 2 (including rotation factor)
     643          kspan = kspan / 2
     644          k1 = kspan + 2
     645          DO
     646             DO
     647                k2 = kk + kspan
     648                ck = array(k2)
     649                array(k2) = array(kk) - ck
     650                array(kk) = array(kk) + ck
     651                kk = k2 + kspan
     652                IF ( kk > nn )  EXIT
     653             END DO
     654             kk = kk - nn
     655             IF ( kk > jc )  EXIT
     656          END DO
     657          IF ( kk > kspan )  RETURN
     658          DO
     659             c1 = 1.0_wp - cd
     660             s1 = sd
     661             DO
     662                DO
     663                   DO
     664                      k2 = kk + kspan
     665                      ck = array(kk) - array(k2)
     666                      array(kk) = array(kk) + array(k2)
     667                      array(k2) = ck * CMPLX( c1, s1, KIND = wp )
     668                      kk = k2 + kspan
     669                      IF ( kk >= nt )  EXIT
     670                   END DO
     671                   k2 = kk - nt
     672                   c1 = - c1
     673                   kk = k1 - k2
     674                   IF ( kk <= k2 )  EXIT
     675                END DO
     676                ak = c1 - (cd * c1 + sd * s1)
     677                s1 = sd * c1 - cd * s1 + s1
     678                c1 = 2.0_wp - ( ak * ak + s1 * s1 )
     679                s1 = s1 * c1
     680                c1 = c1 * ak
     681                kk = kk + jc
     682                IF ( kk >= k2 ) EXIT
     683             END DO
     684             k1 = k1 + 1 + 1
     685             kk = ( k1 - kspan ) / 2 + jc
     686             IF ( kk > jc + jc )  EXIT
     687          END DO
     688!
     689!--    Transform for factor of 4
     690       CASE ( 4 )
     691          ispan = kspan
     692          kspan = kspan / 4
     693
     694          DO
     695             c1 = 1.0_wp
     696             s1 = 0.0_wp
     697             DO
     698                DO
     699                   k1 = kk + kspan
     700                   k2 = k1 + kspan
     701                   k3 = k2 + kspan
     702                   ckp = array(kk) + array(k2)
     703                   ckm = array(kk) - array(k2)
     704                   cjp = array(k1) + array(k3)
     705                   cjm = array(k1) - array(k3)
     706                   array(kk) = ckp + cjp
     707                   cjp = ckp - cjp
     708                   IF ( inv )  THEN
     709                      ckp = ckm + CMPLX( - AIMAG( cjm ), REAL( cjm ), KIND = wp )
     710                      ckm = ckm + CMPLX( AIMAG( cjm ), - REAL( cjm ), KIND = wp )
     711                   ELSE
     712                      ckp = ckm + CMPLX( AIMAG( cjm ), - REAL( cjm ), KIND = wp )
     713                      ckm = ckm + CMPLX( - AIMAG( cjm ), REAL( cjm ), KIND = wp )
     714                   END IF
     715!
     716!--                Avoid useless multiplies
     717                   IF ( s1 == 0.0_wp )  THEN
     718                      array(k1) = ckp
     719                      array(k2) = cjp
     720                      array(k3) = ckm
     721                   ELSE
     722                      array(k1) = ckp * CMPLX( c1, s1, KIND = wp )
     723                      array(k2) = cjp * CMPLX( c2, s2, KIND = wp )
     724                      array(k3) = ckm * CMPLX( c3, s3, KIND = wp )
     725                   END IF
     726                   kk = k3 + kspan
     727                   IF ( kk > nt )  EXIT
     728                END DO
     729
     730                c2 = c1 - ( cd * c1 + sd * s1 )
     731                s1 = sd * c1 - cd * s1 + s1
     732                c1 = 2.0_wp - ( c2 * c2 + s1 * s1 )
     733                s1 = s1 * c1
     734                c1 = c1 * c2
     735!
     736!--             Values of c2, c3, s2, s3 that will get used next time
     737                c2 = c1 * c1 - s1 * s1
     738                s2 = 2.0_wp * c1 * s1
     739                c3 = c2 * c1 - s2 * s1
     740                s3 = c2 * s1 + s2 * c1
     741                kk = kk - nt + jc
     742                IF ( kk > kspan )  EXIT
     743             END DO
     744             kk = kk - kspan + 1
     745             IF ( kk > jc )  EXIT
     746          END DO
     747          IF ( kspan == jc )  RETURN
     748
     749       CASE default
     750!
     751!--       Transform for odd factors
     752          k = factor(ii)
     753          ispan = kspan
     754          kspan = kspan / k
     755
     756          SELECT CASE ( k )
     757!
     758!--       Transform for factor of 3 (optional code)
     759          CASE ( 3 )
     760             DO
     761                DO
     762                   k1 = kk + kspan
     763                   k2 = k1 + kspan
     764                   ck = array(kk)
     765                   cj = array(k1) + array(k2)
     766                   array(kk) = ck + cj
     767                   ck = ck - CMPLX( 0.5_wp * REAL( cj ), 0.5_wp * AIMAG( cj ), KIND = wp )
     768                   cj = CMPLX( ( REAL( array(k1) ) - REAL( array(k2) ) ) * s60,                    &
     769                               ( AIMAG( array(k1) ) - AIMAG( array(k2) ) ) * s60, KIND = wp )
     770                   array(k1) = ck + CMPLX( - AIMAG( cj ), REAL( cj ), KIND = wp )
     771                   array(k2) = ck + CMPLX( AIMAG( cj ), - REAL( cj ), KIND = wp )
     772                   kk = k2 + kspan
     773                   IF ( kk >= nn )  EXIT
     774                END DO
     775                kk = kk - nn
     776                IF ( kk > kspan )  EXIT
     777             END DO
     778!
     779!--       Transform for factor of 5 (optional code)
     780          CASE ( 5 )
     781             c2 = c72 * c72 - s72 * s72
     782             s2 = 2.0_wp * c72 * s72
     783             DO
     784                DO
     785                   k1 = kk + kspan
     786                   k2 = k1 + kspan
     787                   k3 = k2 + kspan
     788                   k4 = k3 + kspan
     789                   ckp = array(k1) + array(k4)
     790                   ckm = array(k1) - array(k4)
     791                   cjp = array(k2) + array(k3)
     792                   cjm = array(k2) - array(k3)
     793                   cc = array(kk)
     794                   array(kk) = cc + ckp + cjp
     795                   ck = CMPLX( REAL( ckp ) * c72, AIMAG( ckp ) * c72, KIND = wp ) +                &
     796                        CMPLX( REAL( cjp ) * c2,  AIMAG( cjp ) * c2,  KIND = wp ) + cc
     797                   cj = CMPLX( REAL( ckm ) * s72, AIMAG( ckm ) * s72, KIND = wp) +                 &
     798                        CMPLX( REAL( cjm ) * s2,  AIMAG( cjm ) * s2,  KIND = wp )
     799                   array(k1) = ck + CMPLX( - AIMAG( cj ), REAL( cj ), KIND = wp )
     800                   array(k4) = ck + CMPLX( AIMAG( cj ), - REAL( cj ), KIND = wp )
     801                   ck = CMPLX( REAL( ckp ) * c2,  AIMAG( ckp ) * c2,  KIND = wp ) +                &
     802                        CMPLX( REAL( cjp ) * c72, AIMAG( cjp ) * c72, KIND = wp ) + cc
     803                   cj = CMPLX( REAL( ckm ) * s2,  AIMAG( ckm ) * s2,  KIND = wp ) -                &
     804                        CMPLX( REAL( cjm ) * s72, AIMAG( cjm ) * s72, KIND = wp )
     805                   array(k2) = ck + CMPLX( -AIMAG( cj ), REAL( cj ), KIND = wp )
     806                   array(k3) = ck + CMPLX( AIMAG( cj ), - REAL( cj ), KIND = wp )
     807                   kk = k4 + kspan
     808                   IF ( kk >= nn )  EXIT
     809                END DO
     810                kk = kk - nn
     811                IF ( kk > kspan )  EXIT
     812             END DO
     813
     814          CASE default
     815             IF ( k /= jf )  THEN
     816                jf = k
     817                s1 = pi2 / k
     818                c1 = COS( s1 )
     819                s1 = SIN( s1 )
     820                cosine(jf) = 1.0_wp
     821                sine(jf) = 0.0_wp
     822                j = 1
     823                DO
     824                   cosine(j) = cosine(k) * c1 + sine(k) * s1
     825                   sine(j) = cosine(k) * s1 - sine(k) * c1
     826                   k = k - 1
     827                   cosine(k) = cosine(j)
     828                   sine(k) = - sine(j)
     829                   j = j + 1
     830                   IF ( j >= k )  EXIT
     831                END DO
     832             END IF
     833             DO
     834                DO
     835                   k1 = kk
     836                   k2 = kk + ispan
     837                   cc = array(kk)
     838                   ck = cc
     839                   j = 1
     840                   k1 = k1 + kspan
     841                   DO
     842                      k2 = k2 - kspan
     843                      j = j + 1
     844                      ctmp(j) = array(k1) + array(k2)
     845                      ck = ck + ctmp(j)
     846                      j = j + 1
     847                      ctmp(j) = array(k1) - array(k2)
     848                      k1 = k1 + kspan
     849                      IF ( k1 >= k2 )  EXIT
     850                   END DO
     851                   array(kk) = ck
     852                   k1 = kk
     853                   k2 = kk + ispan
     854                   j = 1
     855                   DO
     856                      k1 = k1 + kspan
     857                      k2 = k2 - kspan
     858                      jj = j
     859                      ck = cc
     860                      cj = ( 0.0_wp, 0.0_wp )
     861                      k = 1
     862                      DO
     863                         k = k + 1
     864                         ck = ck + CMPLX( REAL( ctmp(k) ) * cosine(jj), AIMAG( ctmp(k) ) *         &
     865                                          cosine(jj), KIND = wp )
     866                         k = k + 1
     867                         cj = cj + CMPLX( REAL( ctmp(k) ) * sine(jj), AIMAG( ctmp(k) ) * sine(jj), &
     868                                          KIND = wp )
     869                         jj = jj + j
     870                         IF ( jj > jf )  jj = jj - jf
     871                         IF ( k >= jf )  EXIT
     872                      END DO
     873                      k = jf - j
     874                      array(k1) = ck + CMPLX( - AIMAG( cj ), REAL( cj ), KIND = wp )
     875                      array(k2) = ck + CMPLX( AIMAG( cj ), -REAL( cj ), KIND = wp )
     876                      j = j + 1
     877                      IF ( j >= k )  EXIT
     878                   END DO
     879                   kk = kk + ispan
     880                   IF ( kk > nn )  EXIT
     881                END DO
     882                kk = kk - nn
     883                IF ( kk > kspan )  EXIT
     884             END DO
     885
     886          END SELECT
     887!
     888!--       Multiply by rotation factor (except for factors of 2 and 4)
     889          IF ( ii == nfactor )  RETURN
     890          kk = jc + 1
     891          DO
     892             c2 = 1.0_wp - cd
     893             s1 = sd
     894             DO
     895                c1 = c2
     896                s2 = s1
     897                kk = kk + kspan
     898                DO
     899                   DO
     900                      array(kk) = CMPLX( c2, s2, KIND = wp ) * array(kk)
     901                      kk = kk + ispan
     902                      IF ( kk > nt )  EXIT
     903                   END DO
     904                   ak = s1 * s2
     905                   s2 = s1 * c2 + c1 * s2
     906                   c2 = c1 * c2 - ak
     907                   kk = kk - nt + kspan
     908                   IF ( kk > ispan )  EXIT
     909                END DO
     910                c2 = c1 - ( cd * c1 + sd * s1 )
     911                s1 = s1 + sd * c1 - cd * s1
     912                c1 = 2.0_wp - ( c2 * c2 + s1 * s1 )
     913                s1 = s1 * c1
     914                c2 = c2 * c1
     915                kk = kk - ispan + jc
     916                IF ( kk > kspan )  EXIT
     917             END DO
     918             kk = kk - kspan + jc + 1
     919             IF ( kk > jc + jc )  EXIT
     920          END DO
     921
     922       END SELECT
     923    END DO
     924 END SUBROUTINE transform
     925
     926
     927!--------------------------------------------------------------------------------------------------!
    948928! Description:
    949929! ------------
    950930!> @todo Missing subroutine description.
    951 !------------------------------------------------------------------------------!
    952     SUBROUTINE permute(array, ntotal, npass, nspan, &
    953          factor, nfactor, nsquare, maxfactor, &
    954          ctmp, perm)
    955 !
    956 !--   Formal parameters
    957       COMPLEX(wp),  DIMENSION(*), INTENT(IN OUT):: array
    958       COMPLEX(wp),  DIMENSION(*), INTENT(OUT)   :: ctmp
    959       INTEGER(iwp),               INTENT(IN)    :: ntotal, npass, nspan
    960       INTEGER(iwp), DIMENSION(*), INTENT(IN OUT):: factor
    961       INTEGER(iwp),               INTENT(IN)    :: nfactor, nsquare
    962       INTEGER(iwp),               INTENT(IN)    :: maxfactor
    963       INTEGER(iwp), DIMENSION(*), INTENT(OUT)   :: perm
    964 !
    965 !--   Local scalars
    966       COMPLEX(wp) :: ck
    967       INTEGER(iwp):: ii, ispan
    968       INTEGER(iwp):: j, jc, jj
    969       INTEGER(iwp):: k, kk, kspan, kt, k1, k2, k3
    970       INTEGER(iwp):: nn, nt
    971 !
    972 !--   Permute the results to normal order---done in two stages
    973 !--   Permutation for square factors of n
    974 
    975       nt = ntotal
    976       nn = nt - 1
    977       kt = nsquare
    978       kspan = nspan
    979       jc = nspan / npass
    980 
    981       perm (1) = nspan
    982       IF (kt > 0) THEN
    983          k = kt + kt + 1
    984          IF (nfactor < k) k = k - 1
    985          j = 1
    986          perm (k + 1) = jc
    987          DO
    988             perm (j + 1) = perm (j) / factor(j)
    989             perm (k) = perm (k + 1) * factor(j)
    990             j = j + 1
    991             k = k - 1
    992             IF (j >= k) EXIT
    993          END DO
    994          k3 = perm (k + 1)
    995          kspan = perm (2)
    996          kk = jc + 1
    997          k2 = kspan + 1
    998          j = 1
    999 
    1000          IF (npass /= ntotal) THEN
    1001             permute_multi: DO
    1002                DO
    1003                   DO
    1004                      k = kk + jc
    1005                      DO
    1006 !
    1007 !--                     Swap array(kk) <> array(k2)
    1008                         ck = array(kk)
    1009                         array(kk) = array(k2)
    1010                         array(k2) = ck
    1011                         kk = kk + 1
    1012                         k2 = k2 + 1
    1013                         IF (kk >= k) EXIT
    1014                      END DO
    1015                      kk = kk + nspan - jc
    1016                      k2 = k2 + nspan - jc
    1017                      IF (kk >= nt) EXIT
    1018                   END DO
    1019                   kk = kk - nt + jc
    1020                   k2 = k2 - nt + kspan
    1021                   IF (k2 >= nspan) EXIT
    1022                END DO
    1023                DO
    1024                   DO
    1025                      k2 = k2 - perm (j)
    1026                      j = j + 1
    1027                      k2 = perm (j + 1) + k2
    1028                      IF (k2 <= perm (j)) EXIT
    1029                   END DO
    1030                   j = 1
    1031                   DO
    1032                      IF (kk < k2) CYCLE permute_multi
    1033                      kk = kk + jc
    1034                      k2 = k2 + kspan
    1035                      IF (k2 >= nspan) EXIT
    1036                   END DO
    1037                   IF (kk >= nspan) EXIT
    1038                END DO
    1039                EXIT
    1040             END DO permute_multi
    1041          ELSE
    1042             permute_single: DO
    1043                DO
    1044 !
    1045 !--               Swap array(kk) <> array(k2)
    1046                   ck = array(kk)
    1047                   array(kk) = array(k2)
    1048                   array(k2) = ck
    1049                   kk = kk + 1
    1050                   k2 = k2 + kspan
    1051                   IF (k2 >= nspan) EXIT
    1052                END DO
    1053                DO
    1054                   DO
    1055                      k2 = k2 - perm (j)
    1056                      j = j + 1
    1057                      k2 = perm (j + 1) + k2
    1058                      IF (k2 <= perm (j)) EXIT
    1059                   END DO
    1060                   j = 1
    1061                   DO
    1062                      IF (kk < k2) CYCLE permute_single
    1063                      kk = kk + 1
    1064                      k2 = k2 + kspan
    1065                      IF (k2 >= nspan) EXIT
    1066                   END DO
    1067                   IF (kk >= nspan) EXIT
    1068                END DO
    1069                EXIT
    1070             END DO permute_single
    1071          END IF
    1072          jc = k3
    1073       END IF
    1074 
    1075       IF (ISHFT(kt, 1) + 1 >= nfactor) RETURN
    1076 
    1077       ispan = perm (kt + 1)
    1078 !
    1079 !--   Permutation for square-free factors of n
    1080       j = nfactor - kt
    1081       factor(j + 1) = 1
    1082       DO
    1083          factor(j) = factor(j) * factor(j+1)
    1084          j = j - 1
    1085          IF (j == kt) EXIT
    1086       END DO
    1087       kt = kt + 1
    1088       nn = factor(kt) - 1
    1089       j = 0
    1090       jj = 0
    1091       DO
    1092          k = kt + 1
    1093          k2 = factor(kt)
    1094          kk = factor(k)
    1095          j = j + 1
    1096          IF (j > nn) EXIT !-- exit infinite loop
    1097          jj = jj + kk
    1098          DO WHILE (jj >= k2)
    1099             jj = jj - k2
    1100             k2 = kk
    1101             k = k + 1
    1102             kk = factor(k)
    1103             jj = jj + kk
    1104          END DO
    1105          perm (j) = jj
    1106       END DO
    1107 !
    1108 !--   Determine the permutation cycles of length greater than 1
    1109       j = 0
    1110       DO
    1111          DO
    1112             j = j + 1
    1113             kk = perm(j)
    1114             IF (kk >= 0) EXIT
    1115          END DO
    1116          IF (kk /= j) THEN
    1117             DO
    1118                k = kk
    1119                kk = perm (k)
    1120                perm (k) = -kk
    1121                IF (kk == j) EXIT
    1122             END DO
    1123             k3 = kk
    1124          ELSE
    1125             perm (j) = -j
    1126             IF (j == nn) EXIT !-- exit infinite loop
    1127          END IF
    1128       END DO
    1129 !
    1130 !--   Reorder a and b, following the permutation cycles
    1131       DO
    1132          j = k3 + 1
    1133          nt = nt - ispan
    1134          ii = nt - 1 + 1
    1135          IF (nt < 0) EXIT !-- exit infinite loop
    1136          DO
    1137             DO
    1138                j = j-1
    1139                IF (perm(j) >= 0) EXIT
    1140             END DO
    1141             jj = jc
    1142             DO
    1143                kspan = jj
    1144                IF (jj > maxfactor) kspan = maxfactor
    1145                jj = jj - kspan
    1146                k = perm(j)
    1147                kk = jc * k + ii + jj
    1148                k1 = kk + kspan
    1149                k2 = 0
    1150                DO
    1151                   k2 = k2 + 1
    1152                   ctmp(k2) = array(k1)
    1153                   k1 = k1 - 1
    1154                   IF (k1 == kk) EXIT
    1155                END DO
    1156                DO
    1157                   k1 = kk + kspan
    1158                   k2 = k1 - jc * (k + perm(k))
    1159                   k = -perm(k)
    1160                   DO
    1161                      array(k1) = array(k2)
    1162                      k1 = k1 - 1
    1163                      k2 = k2 - 1
    1164                      IF (k1 == kk) EXIT
    1165                   END DO
    1166                   kk = k2
    1167                   IF (k == j) EXIT
    1168                END DO
    1169                k1 = kk + kspan
    1170                k2 = 0
    1171                DO
    1172                   k2 = k2 + 1
    1173                   array(k1) = ctmp(k2)
    1174                   k1 = k1 - 1
    1175                   IF (k1 == kk) EXIT
    1176                END DO
    1177                IF (jj == 0) EXIT
    1178             END DO
    1179             IF (j == 1) EXIT
    1180          END DO
    1181       END DO
    1182 
    1183     END SUBROUTINE permute
     931!--------------------------------------------------------------------------------------------------!
     932 SUBROUTINE permute( array, ntotal, npass, nspan, factor, nfactor, nsquare, maxfactor, ctmp, perm )
     933!
     934!-- Formal parameters
     935    COMPLEX(wp),  DIMENSION(*), INTENT(IN OUT) ::  array                 !<
     936    COMPLEX(wp),  DIMENSION(*), INTENT(OUT)    ::  ctmp                  !<
     937    INTEGER(iwp),               INTENT(IN)     ::  ntotal, npass, nspan  !<
     938    INTEGER(iwp),               INTENT(IN)     ::  nfactor, nsquare      !<
     939    INTEGER(iwp),               INTENT(IN)     ::  maxfactor             !<
     940    INTEGER(iwp), DIMENSION(*), INTENT(IN OUT) ::  factor                !<
     941    INTEGER(iwp), DIMENSION(*), INTENT(OUT)    ::  perm                  !<
     942!
     943!-- Local scalars
     944    COMPLEX(wp)  ::  ck                            !<
     945    INTEGER(iwp) ::  ii, ispan                     !<
     946    INTEGER(iwp) ::  j, jc, jj                     !<
     947    INTEGER(iwp) ::  k, kk, kspan, kt, k1, k2, k3  !<
     948    INTEGER(iwp) ::  nn, nt                        !<
     949!
     950!-- Permute the results to normal order---done in two stages
     951!-- Permutation for square factors of n
     952
     953    nt = ntotal
     954    nn = nt - 1
     955    kt = nsquare
     956    kspan = nspan
     957    jc = nspan / npass
     958
     959    perm (1) = nspan
     960    IF ( kt > 0 )  THEN
     961       k = kt + kt + 1
     962       IF ( nfactor < k )  k = k - 1
     963       j = 1
     964       perm(k + 1) = jc
     965       DO
     966          perm(j + 1) = perm(j) / factor(j)
     967          perm(k) = perm(k + 1) * factor(j)
     968          j = j + 1
     969          k = k - 1
     970          IF ( j >= k )  EXIT
     971       END DO
     972       k3 = perm(k + 1)
     973       kspan = perm(2)
     974       kk = jc + 1
     975       k2 = kspan + 1
     976       j = 1
     977
     978       IF ( npass /= ntotal )  THEN
     979          permute_multi: DO
     980             DO
     981                DO
     982                   k = kk + jc
     983                   DO
     984!
     985!--                   Swap array(kk) <> array(k2)
     986                      ck = array(kk)
     987                      array(kk) = array(k2)
     988                      array(k2) = ck
     989                      kk = kk + 1
     990                      k2 = k2 + 1
     991                      IF ( kk >= k )  EXIT
     992                   END DO
     993                   kk = kk + nspan - jc
     994                   k2 = k2 + nspan - jc
     995                   IF ( kk >= nt )  EXIT
     996                END DO
     997                kk = kk - nt + jc
     998                k2 = k2 - nt + kspan
     999                IF ( k2 >= nspan )  EXIT
     1000             END DO
     1001             DO
     1002                DO
     1003                   k2 = k2 - perm(j)
     1004                   j = j + 1
     1005                   k2 = perm(j + 1) + k2
     1006                   IF ( k2 <= perm(j) )  EXIT
     1007                END DO
     1008                j = 1
     1009                DO
     1010                   IF ( kk < k2 )  CYCLE permute_multi
     1011                   kk = kk + jc
     1012                   k2 = k2 + kspan
     1013                   IF ( k2 >= nspan )  EXIT
     1014                END DO
     1015                IF ( kk >= nspan )  EXIT
     1016             END DO
     1017             EXIT
     1018          END DO permute_multi
     1019       ELSE
     1020          permute_single: DO
     1021             DO
     1022!
     1023!--             Swap array(kk) <> array(k2)
     1024                ck = array(kk)
     1025                array(kk) = array(k2)
     1026                array(k2) = ck
     1027                kk = kk + 1
     1028                k2 = k2 + kspan
     1029                IF ( k2 >= nspan )  EXIT
     1030             END DO
     1031             DO
     1032                DO
     1033                   k2 = k2 - perm(j)
     1034                   j = j + 1
     1035                   k2 = perm(j + 1) + k2
     1036                   IF ( k2 <= perm(j) )  EXIT
     1037                END DO
     1038                j = 1
     1039                DO
     1040                   IF ( kk < k2 )  CYCLE permute_single
     1041                   kk = kk + 1
     1042                   k2 = k2 + kspan
     1043                   IF ( k2 >= nspan )  EXIT
     1044                END DO
     1045                IF ( kk >= nspan )  EXIT
     1046             END DO
     1047             EXIT
     1048          END DO permute_single
     1049       END IF
     1050       jc = k3
     1051    END IF
     1052
     1053    IF ( ISHFT( kt, 1 ) + 1 >= nfactor )  RETURN
     1054
     1055    ispan = perm(kt + 1)
     1056!
     1057!-- Permutation for square-free factors of n
     1058    j = nfactor - kt
     1059    factor( j + 1 ) = 1
     1060    DO
     1061       factor(j) = factor(j) * factor(j+1)
     1062       j = j - 1
     1063       IF ( j == kt )  EXIT
     1064    END DO
     1065    kt = kt + 1
     1066    nn = factor( kt ) - 1
     1067    j = 0
     1068    jj = 0
     1069    DO
     1070       k = kt + 1
     1071       k2 = factor(kt)
     1072       kk = factor(k)
     1073       j = j + 1
     1074       IF ( j > nn )  EXIT !-- Exit infinite loop
     1075       jj = jj + kk
     1076       DO WHILE ( jj >= k2 )
     1077          jj = jj - k2
     1078          k2 = kk
     1079          k = k + 1
     1080          kk = factor(k)
     1081          jj = jj + kk
     1082       END DO
     1083       perm(j) = jj
     1084    END DO
     1085!
     1086!-- Determine the permutation cycles of length greater than 1
     1087    j = 0
     1088    DO
     1089       DO
     1090          j = j + 1
     1091          kk = perm(j)
     1092          IF ( kk >= 0 )  EXIT
     1093       END DO
     1094       IF ( kk /= j )  THEN
     1095          DO
     1096             k = kk
     1097             kk = perm(k)
     1098             perm(k) = - kk
     1099             IF ( kk == j )  EXIT
     1100          END DO
     1101          k3 = kk
     1102       ELSE
     1103          perm(j) = - j
     1104          IF ( j == nn )  EXIT !-- Exit infinite loop
     1105       END IF
     1106    END DO
     1107!
     1108!-- Reorder a and b, following the permutation cycles
     1109    DO
     1110       j = k3 + 1
     1111       nt = nt - ispan
     1112       ii = nt - 1 + 1
     1113       IF ( nt < 0 )  EXIT !-- Exit infinite loop
     1114       DO
     1115          DO
     1116             j = j - 1
     1117             IF ( perm(j) >= 0 )  EXIT
     1118          END DO
     1119          jj = jc
     1120          DO
     1121             kspan = jj
     1122             IF ( jj > maxfactor )  kspan = maxfactor
     1123             jj = jj - kspan
     1124             k = perm(j)
     1125             kk = jc * k + ii + jj
     1126             k1 = kk + kspan
     1127             k2 = 0
     1128             DO
     1129                k2 = k2 + 1
     1130                ctmp(k2) = array(k1)
     1131                k1 = k1 - 1
     1132                IF ( k1 == kk )  EXIT
     1133             END DO
     1134             DO
     1135                k1 = kk + kspan
     1136                k2 = k1 - jc * ( k + perm(k) )
     1137                k = - perm(k)
     1138                DO
     1139                   array(k1) = array(k2)
     1140                   k1 = k1 - 1
     1141                   k2 = k2 - 1
     1142                   IF ( k1 == kk )  EXIT
     1143                END DO
     1144                kk = k2
     1145                IF ( k == j )  EXIT
     1146             END DO
     1147             k1 = kk + kspan
     1148             k2 = 0
     1149             DO
     1150                k2 = k2 + 1
     1151                array(k1) = ctmp(k2)
     1152                k1 = k1 - 1
     1153                IF ( k1 == kk )  EXIT
     1154             END DO
     1155             IF ( jj == 0 )  EXIT
     1156          END DO
     1157          IF ( j == 1 )  EXIT
     1158       END DO
     1159    END DO
     1160
     1161 END SUBROUTINE permute
    11841162
    11851163 END SUBROUTINE fftradix
  • palm/trunk/SOURCE/sor.f90

    r4457 r4591  
    11!> @file sor.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:
    2121! -----------------
    22 ! 
    23 ! 
     22!
     23!
    2424! Former revisions:
    2525! -----------------
    2626! $Id$
    27 ! use statement for exchange horiz added
    28 !
     27! File re-formatted to follow the PALM coding standard
     28!
     29!
     30! 4457 2020-03-11 14:20:43Z raasch
     31! Use statement for exchange horiz added
     32!
    2933! 4360 2020-01-07 11:25:50Z suehring
    3034! Corrected "Former revisions" section
    31 ! 
     35!
    3236! 3655 2019-01-07 16:51:22Z knoop
    3337! Rename variables in mesoscale-offline nesting mode
     
    3640! Initial revision
    3741!
    38 !
     42!--------------------------------------------------------------------------------------------------!
    3943! Description:
    4044! ------------
    4145!> Solve the Poisson-equation with the SOR-Red/Black-scheme.
    42 !------------------------------------------------------------------------------!
     46!--------------------------------------------------------------------------------------------------!
    4347 SUBROUTINE sor( d, ddzu, ddzw, p )
    4448
    45     USE arrays_3d,                                                             &
    46         ONLY:  rho_air, rho_air_zw
    47 
    48     USE control_parameters,                                                    &
    49         ONLY:  bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r,                 &
    50                bc_dirichlet_s, bc_lr_cyc, bc_ns_cyc, bc_radiation_l,           &
    51                bc_radiation_n, bc_radiation_r, bc_radiation_s, ibc_p_b,        &
    52                ibc_p_t, n_sor, omega_sor
    53 
    54     USE exchange_horiz_mod,                                                    &
     49    USE arrays_3d,                                                                                 &
     50        ONLY:  rho_air,                                                                            &
     51               rho_air_zw
     52
     53    USE control_parameters,                                                                        &
     54        ONLY:  bc_dirichlet_l,                                                                     &
     55               bc_dirichlet_n,                                                                     &
     56               bc_dirichlet_r,                                                                     &
     57               bc_dirichlet_s,                                                                     &
     58               bc_lr_cyc,                                                                          &
     59               bc_ns_cyc,                                                                          &
     60               bc_radiation_l,                                                                     &
     61               bc_radiation_n,                                                                     &
     62               bc_radiation_r,                                                                     &
     63               bc_radiation_s,                                                                     &
     64               ibc_p_b,                                                                            &
     65               ibc_p_t,                                                                            &
     66               n_sor,                                                                              &
     67               omega_sor
     68
     69    USE exchange_horiz_mod,                                                                        &
    5570        ONLY:  exchange_horiz
    5671
    57     USE grid_variables,                                                        &
    58         ONLY:  ddx2, ddy2
    59 
    60     USE indices,                                                               &
    61         ONLY:  nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nz, nzb, nzt
     72    USE grid_variables,                                                                            &
     73        ONLY:  ddx2,                                                                               &
     74               ddy2
     75
     76    USE indices,                                                                                   &
     77        ONLY:  nbgp,                                                                               &
     78               nxl,                                                                                &
     79               nxlg,                                                                               &
     80               nxr,                                                                                &
     81               nxrg,                                                                               &
     82               nyn,                                                                                &
     83               nyng,                                                                               &
     84               nys,                                                                                &
     85               nysg,                                                                               &
     86               nz,                                                                                 &
     87               nzb,                                                                                &
     88               nzt
    6289
    6390    USE kinds
     
    6592    IMPLICIT NONE
    6693
    67     INTEGER(iwp) ::  i              !<
    68     INTEGER(iwp) ::  j              !<
    69     INTEGER(iwp) ::  k              !<
    70     INTEGER(iwp) ::  n              !<
    71     INTEGER(iwp) ::  nxl1           !<
    72     INTEGER(iwp) ::  nxl2           !<
    73     INTEGER(iwp) ::  nys1           !<
    74     INTEGER(iwp) ::  nys2           !<
    75 
    76     REAL(wp)     ::  ddzu(1:nz+1)   !<
    77     REAL(wp)     ::  ddzw(1:nzt+1)  !<
    78 
    79     REAL(wp)     ::  d(nzb+1:nzt,nys:nyn,nxl:nxr)      !<
    80     REAL(wp)     ::  p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  !<
    81 
    82     REAL(wp), DIMENSION(:), ALLOCATABLE ::  f1         !<
    83     REAL(wp), DIMENSION(:), ALLOCATABLE ::  f2         !<
    84     REAL(wp), DIMENSION(:), ALLOCATABLE ::  f3         !<
     94    INTEGER(iwp) ::  i     !<
     95    INTEGER(iwp) ::  j     !<
     96    INTEGER(iwp) ::  k     !<
     97    INTEGER(iwp) ::  n     !<
     98    INTEGER(iwp) ::  nxl1  !<
     99    INTEGER(iwp) ::  nxl2  !<
     100    INTEGER(iwp) ::  nys1  !<
     101    INTEGER(iwp) ::  nys2  !<
     102
     103    REAL(wp) ::  ddzu(1:nz+1)   !<
     104    REAL(wp) ::  ddzw(1:nzt+1)  !<
     105
     106    REAL(wp) ::  d(nzb+1:nzt,nys:nyn,nxl:nxr)      !<
     107    REAL(wp) ::  p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  !<
     108
     109    REAL(wp), DIMENSION(:), ALLOCATABLE ::  f1  !<
     110    REAL(wp), DIMENSION(:), ALLOCATABLE ::  f2  !<
     111    REAL(wp), DIMENSION(:), ALLOCATABLE ::  f3  !<
    85112
    86113    ALLOCATE( f1(1:nz), f2(1:nz), f3(1:nz) )
     
    118145          DO  j = nys2, nyn, 2
    119146             DO  k = nzb+1, nzt
    120                 p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * (            &
    121                            rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +   &
    122                            rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +   &
    123                            f2(k) * p(k+1,j,i)                              +   &
    124                            f3(k) * p(k-1,j,i)                              -   &
    125                            d(k,j,i)                                        -   &
    126                            f1(k) * p(k,j,i)           )
     147                p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * (                                        &
     148                           rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +                       &
     149                           rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +                       &
     150                           f2(k) * p(k+1,j,i)                              +                       &
     151                           f3(k) * p(k-1,j,i)                              -                       &
     152                           d(k,j,i)                                        -                       &
     153                           f1(k) * p(k,j,i)               )
    127154             ENDDO
    128155          ENDDO
     
    132159          DO  j = nys1, nyn, 2
    133160             DO  k = nzb+1, nzt
    134                 p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * (                    &
    135                            rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +   &
    136                            rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +   &
    137                            f2(k) * p(k+1,j,i)                              +   &
    138                            f3(k) * p(k-1,j,i)                              -   &
    139                            d(k,j,i)                                        -   &
    140                            f1(k) * p(k,j,i)           )
     161                p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * (                                        &
     162                           rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +                       &
     163                           rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +                       &
     164                           f2(k) * p(k+1,j,i)                              +                       &
     165                           f3(k) * p(k-1,j,i)                              -                       &
     166                           d(k,j,i)                                        -                       &
     167                           f1(k) * p(k,j,i)               )
    141168             ENDDO
    142169          ENDDO
     
    146173!--    Exchange of boundary values for p.
    147174       CALL exchange_horiz( p, nbgp )
     175
    148176
    149177!
     
    163191          DO  j = nys1, nyn, 2
    164192             DO  k = nzb+1, nzt
    165                 p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * (            &
    166                            rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +   &
    167                            rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +   &
    168                            f2(k) * p(k+1,j,i)                              +   &
    169                            f3(k) * p(k-1,j,i)                              -   &
    170                            d(k,j,i)                                        -   &
    171                            f1(k) * p(k,j,i)           )
     193                p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * (                                        &
     194                           rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +                       &
     195                           rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +                       &
     196                           f2(k) * p(k+1,j,i)                              +                       &
     197                           f3(k) * p(k-1,j,i)                              -                       &
     198                           d(k,j,i)                                        -                       &
     199                           f1(k) * p(k,j,i)               )
    172200             ENDDO
    173201          ENDDO
     
    177205          DO  j = nys2, nyn, 2
    178206             DO  k = nzb+1, nzt
    179                 p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * (            &
    180                            rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +   &
    181                            rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +   &
    182                            f2(k) * p(k+1,j,i)                              +   &
    183                            f3(k) * p(k-1,j,i)                              -   &
    184                            d(k,j,i)                                        -   &
    185                            f1(k) * p(k,j,i)           )
     207                p(k,j,i) = p(k,j,i) + omega_sor / f1(k) * (                                        &
     208                           rho_air(k) * ddx2 * ( p(k,j,i+1) + p(k,j,i-1) ) +                       &
     209                           rho_air(k) * ddy2 * ( p(k,j+1,i) + p(k,j-1,i) ) +                       &
     210                           f2(k) * p(k+1,j,i)                              +                       &
     211                           f3(k) * p(k-1,j,i)                              -                       &
     212                           d(k,j,i)                                        -                       &
     213                           f1(k) * p(k,j,i)               )
    186214             ENDDO
    187215          ENDDO
     
    195223!--    Boundary conditions top/bottom.
    196224!--    Bottom boundary
    197        IF ( ibc_p_b == 1 )  THEN       !       Neumann
     225       IF ( ibc_p_b == 1 )  THEN       ! Neumann
    198226          p(nzb,:,:) = p(nzb+1,:,:)
    199        ELSE                            !       Dirichlet
     227       ELSE                            ! Dirichlet
    200228          p(nzb,:,:) = 0.0_wp
    201229       ENDIF
     
    203231!
    204232!--    Top boundary
    205        IF ( ibc_p_t == 1 )  THEN                 ! Neumann
     233       IF ( ibc_p_t == 1 )  THEN       ! Neumann
    206234          p(nzt+1,:,:) = p(nzt,:,:)
    207        ELSE                      ! Dirichlet
     235       ELSE                            ! Dirichlet
    208236          p(nzt+1,:,:) = 0.0_wp
    209237       ENDIF
  • palm/trunk/SOURCE/spectra_mod.f90

    r4429 r4591  
    11!> @file spectra_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:
    2121! -----------------
    22 ! 
    23 ! 
     22!
     23!
    2424! Former revisions:
    2525! -----------------
    2626! $Id$
    27 ! bugfix: preprocessor directives rearranged for serial mode
    28 !
     27! File re-formatted to follow the PALM coding standard
     28!
     29! 4429 2020-02-27 15:24:30Z raasch
     30! Bugfix: preprocessor directives rearranged for serial mode
     31!
    2932! 4360 2020-01-07 11:25:50Z suehring
    3033! Corrected "Former revisions" section
    31 ! 
     34!
    3235! 3805 2019-03-20 15:26:35Z raasch
    33 ! unused variable removed
    34 ! 
     36! Unused variable removed
     37!
    3538! 3655 2019-01-07 16:51:22Z knoop
    3639! Renamed output variables
     
    3942! Initial revision
    4043!
    41 !
     44!--------------------------------------------------------------------------------------------------!
    4245! Description:
    4346! ------------
    44 !> Calculate horizontal spectra along x and y.
    45 !> ATTENTION: 1d-decomposition along y still needs improvement, because in that
    46 !>            case the gridpoint number along z still depends on the PE number
    47 !>            because transpose_xz has to be used (and possibly also
    48 !>            transpose_zyd needs modification).
    49 !------------------------------------------------------------------------------!
     47!> Calculate horizontal spectra along x and y.
     48!> ATTENTION: 1d-decomposition along y still needs improvement, because in that case the gridpoint
     49!>            number along z still depends on the PE number because transpose_xz has to be used
     50!>            (and possibly also transpose_zyd needs modification).
     51!--------------------------------------------------------------------------------------------------!
    5052 MODULE spectra_mod
    5153
     
    5456    PRIVATE
    5557
    56     CHARACTER (LEN=2),  DIMENSION(10) ::  spectra_direction = 'x'
    57     CHARACTER (LEN=10), DIMENSION(10) ::  data_output_sp  = ' '
    58 
    59     INTEGER(iwp) ::  average_count_sp = 0
    60     INTEGER(iwp) ::  dosp_time_count = 0
    61     INTEGER(iwp) ::  n_sp_x = 0, n_sp_y = 0
    62 
    63     INTEGER(iwp) ::  comp_spectra_level(100) = 999999
     58    CHARACTER (LEN=10), DIMENSION(10) ::  data_output_sp  = ' '    !<
     59    CHARACTER (LEN=2),  DIMENSION(10) ::  spectra_direction = 'x'  !<
     60
     61    INTEGER(iwp) ::  average_count_sp = 0              !<
     62    INTEGER(iwp) ::  comp_spectra_level(100) = 999999  !<
     63    INTEGER(iwp) ::  dosp_time_count = 0               !<
     64    INTEGER(iwp) ::  n_sp_x = 0, n_sp_y = 0            !<
    6465
    6566    LOGICAL ::  calculate_spectra   = .FALSE.  !< internal switch that spectra are calculated
     
    7071    REAL(wp) ::  skip_time_dosp = 9999999.9_wp         !< no output of spectra data before this interval has passed
    7172
    72     REAL(wp), DIMENSION(:), ALLOCATABLE ::  var_d
    73 
    74     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  spectrum_x, spectrum_y
     73    REAL(wp), DIMENSION(:), ALLOCATABLE ::  var_d  !<
     74
     75    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  spectrum_x  !<
     76    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  spectrum_y  !<
    7577
    7678    SAVE
     
    110112    END INTERFACE spectra_parin
    111113
    112     PUBLIC average_count_sp, averaging_interval_sp, calc_spectra,              &
    113            calculate_spectra, comp_spectra_level, data_output_sp,              &
    114            dosp_time_count, dt_dosp, n_sp_x, n_sp_y,                           &
    115            skip_time_dosp, spectra_check_parameters, spectra_direction,        &
    116            spectra_header, spectra_init, spectra_parin, spectrum_x,            &
    117            spectrum_y, var_d
     114    PUBLIC average_count_sp,                                                                       &
     115           averaging_interval_sp,                                                                  &
     116           calc_spectra,                                                                           &
     117           calculate_spectra,                                                                      &
     118           comp_spectra_level,                                                                     &
     119           data_output_sp,                                                                         &
     120           dosp_time_count,                                                                        &
     121           dt_dosp,                                                                                &
     122           n_sp_x,                                                                                 &
     123           n_sp_y,                                                                                 &
     124           skip_time_dosp,                                                                         &
     125           spectra_check_parameters,                                                               &
     126           spectra_direction,                                                                      &
     127           spectra_header,                                                                         &
     128           spectra_init,                                                                           &
     129           spectra_parin,                                                                          &
     130           spectrum_x,                                                                             &
     131           spectrum_y,                                                                             &
     132           var_d
    118133
    119134
    120135 CONTAINS
    121136
    122 !------------------------------------------------------------------------------!
     137!--------------------------------------------------------------------------------------------------!
    123138! Description:
    124139! ------------
    125140!> Parin for &spectra_par for calculating spectra
    126 !------------------------------------------------------------------------------!
    127     SUBROUTINE spectra_parin
    128 
    129        USE control_parameters,                                                 &
    130            ONLY:  dt_data_output, message_string
    131 
    132        IMPLICIT NONE
    133 
    134        CHARACTER (LEN=80) ::  line  !< dummy string that contains the current  &
    135                                     !< line of the parameter file
    136 
    137        NAMELIST /spectra_par/  averaging_interval_sp, comp_spectra_level,      &
    138                                data_output_sp, dt_dosp, skip_time_dosp,        &
    139                                spectra_direction
    140 
    141        NAMELIST /spectra_parameters/                                           &
    142                                averaging_interval_sp, comp_spectra_level,      &
    143                                data_output_sp, dt_dosp, skip_time_dosp,        &
    144                                spectra_direction
    145 !
    146 !--    Position the namelist-file at the beginning (it was already opened in
    147 !--    parin), search for the namelist-group of the package and position the
    148 !--    file at this line.
    149        line = ' '
    150 
    151 !
    152 !--    Try to find the spectra package
    153        REWIND ( 11 )
    154        line = ' '
    155        DO WHILE ( INDEX( line, '&spectra_parameters' ) == 0 )
    156           READ ( 11, '(A)', END=12 )  line
    157        ENDDO
    158        BACKSPACE ( 11 )
    159 
    160 !
    161 !--    Read namelist
    162        READ ( 11, spectra_parameters, ERR = 10 )
    163 
    164 !
    165 !--    Default setting of dt_dosp here (instead of check_parameters), because
    166 !--    its current value is needed in init_pegrid
    167        IF ( dt_dosp == 9999999.9_wp )  dt_dosp = dt_data_output
    168 
    169 !
    170 !--    Set general switch that spectra shall be calculated
    171        calculate_spectra = .TRUE.
    172 
    173        GOTO 14
    174 
    175  10    BACKSPACE( 11 )
    176        READ( 11 , '(A)') line
    177        CALL parin_fail_message( 'spectra_parameters', line )
    178 !
    179 !--    Try to find the old namelist
    180  12    REWIND ( 11 )
    181        line = ' '
    182        DO WHILE ( INDEX( line, '&spectra_par' ) == 0 )
    183           READ ( 11, '(A)', END=14 )  line
    184        ENDDO
    185        BACKSPACE ( 11 )
    186 
    187 !
    188 !--    Read namelist
    189        READ ( 11, spectra_par, ERR = 13, END = 14 )
    190 
    191        
    192        message_string = 'namelist spectra_par is deprecated and will be ' // &
    193                      'removed in near future. Please use namelist ' //       &
    194                      'spectra_parameters instead'
    195        CALL message( 'spectra_parin', 'PA0487', 0, 1, 0, 6, 0 )
    196 !
    197 !--    Default setting of dt_dosp here (instead of check_parameters), because
    198 !--    its current value is needed in init_pegrid
    199        IF ( dt_dosp == 9999999.9_wp )  dt_dosp = dt_data_output
    200 
    201 !
    202 !--    Set general switch that spectra shall be calculated
    203        calculate_spectra = .TRUE.
    204        
    205        GOTO 14
    206 
    207  13    BACKSPACE( 11 )
    208        READ( 11 , '(A)') line
    209        CALL parin_fail_message( 'spectra_par', line )
    210        
    211        
    212  14    CONTINUE
    213 
    214     END SUBROUTINE spectra_parin
    215 
    216 
    217 
    218 !------------------------------------------------------------------------------!
     141!--------------------------------------------------------------------------------------------------!
     142 SUBROUTINE spectra_parin
     143
     144    USE control_parameters,                                                                        &
     145        ONLY:  dt_data_output,                                                                     &
     146               message_string
     147
     148    IMPLICIT NONE
     149
     150    CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file
     151
     152    NAMELIST /spectra_par/  averaging_interval_sp,                                                 &
     153                            comp_spectra_level,                                                    &
     154                            data_output_sp,                                                        &
     155                            dt_dosp,                                                               &
     156                            skip_time_dosp,                                                        &
     157                            spectra_direction
     158
     159    NAMELIST /spectra_parameters/                                                                  &
     160                            averaging_interval_sp,                                                 &
     161                            comp_spectra_level,                                                    &
     162                            data_output_sp,                                                        &
     163                            dt_dosp,                                                               &
     164                            skip_time_dosp,                                                        &
     165                            spectra_direction
     166!
     167!-- Position the namelist-file at the beginning (it was already opened in parin), search for the
     168!-- namelist-group of the package and position the file at this line.
     169    line = ' '
     170
     171!
     172!-- Try to find the spectra package
     173    REWIND ( 11 )
     174    line = ' '
     175    DO WHILE ( INDEX( line, '&spectra_parameters' ) == 0 )
     176       READ ( 11, '(A)', END=12 )  line
     177    ENDDO
     178    BACKSPACE ( 11 )
     179
     180!
     181!-- Read namelist
     182    READ ( 11, spectra_parameters, ERR = 10 )
     183
     184!
     185!-- Default setting of dt_dosp here (instead of check_parameters), because its current value is
     186!-- needed in init_pegrid
     187    IF ( dt_dosp == 9999999.9_wp )  dt_dosp = dt_data_output
     188
     189!
     190!-- Set general switch that spectra shall be calculated
     191    calculate_spectra = .TRUE.
     192
     193    GOTO 14
     194
     195 10 BACKSPACE( 11 )
     196    READ( 11 , '(A)') line
     197    CALL parin_fail_message( 'spectra_parameters', line )
     198!
     199!-- Try to find the old namelist
     200 12 REWIND ( 11 )
     201    line = ' '
     202    DO WHILE ( INDEX( line, '&spectra_par' ) == 0 )
     203       READ ( 11, '(A)', END=14 )  line
     204    ENDDO
     205    BACKSPACE ( 11 )
     206
     207!
     208!-- Read namelist
     209    READ ( 11, spectra_par, ERR = 13, END = 14 )
     210
     211
     212    message_string = 'namelist spectra_par is deprecated and will be removed in near future.' //   &
     213                     ' Please use namelist spectra_parameters instead'
     214    CALL message( 'spectra_parin', 'PA0487', 0, 1, 0, 6, 0 )
     215!
     216!-- Default setting of dt_dosp here (instead of check_parameters), because its current value is
     217!-- needed in init_pegrid
     218    IF ( dt_dosp == 9999999.9_wp )  dt_dosp = dt_data_output
     219
     220!
     221!-- Set general switch that spectra shall be calculated
     222    calculate_spectra = .TRUE.
     223
     224    GOTO 14
     225
     226 13 BACKSPACE( 11 )
     227    READ( 11 , '(A)') line
     228    CALL parin_fail_message( 'spectra_par', line )
     229
     230
     231 14 CONTINUE
     232
     233 END SUBROUTINE spectra_parin
     234
     235
     236
     237!--------------------------------------------------------------------------------------------------!
    219238! Description:
    220239! ------------
    221240!> Initialization of spectra related variables
    222 !------------------------------------------------------------------------------!
    223     SUBROUTINE spectra_init
    224 
    225        USE indices,                                                            &
    226            ONLY:  nx, ny, nzb, nzt
    227 
    228        IMPLICIT NONE
    229 
    230        IF ( spectra_initialized )  RETURN
    231 
    232        IF ( dt_dosp /= 9999999.9_wp )  THEN
    233 
    234           IF ( .NOT. ALLOCATED( spectrum_x ) )  THEN
    235              ALLOCATE( spectrum_x( 1:nx/2, 1:100, 1:10 ) )
    236              spectrum_x = 0.0_wp
    237           ENDIF
    238 
    239           IF ( .NOT. ALLOCATED( spectrum_y ) )  THEN
    240              ALLOCATE( spectrum_y( 1:ny/2, 1:100, 1:10 ) )
    241              spectrum_y = 0.0_wp
    242           ENDIF
    243 
    244           ALLOCATE( var_d(nzb:nzt+1) )
    245           var_d = 0.0_wp
     241!--------------------------------------------------------------------------------------------------!
     242 SUBROUTINE spectra_init
     243
     244    USE indices,                                                                                   &
     245        ONLY:  nx,                                                                                 &
     246               ny,                                                                                 &
     247               nzb,                                                                                &
     248               nzt
     249
     250    IMPLICIT NONE
     251
     252    IF ( spectra_initialized )  RETURN
     253
     254    IF ( dt_dosp /= 9999999.9_wp )  THEN
     255
     256       IF ( .NOT. ALLOCATED( spectrum_x ) )  THEN
     257          ALLOCATE( spectrum_x( 1:nx/2, 1:100, 1:10 ) )
     258          spectrum_x = 0.0_wp
    246259       ENDIF
    247260
    248        spectra_initialized = .TRUE.
    249 
    250     END SUBROUTINE spectra_init
    251 
    252 
    253 
    254 !------------------------------------------------------------------------------!
     261       IF ( .NOT. ALLOCATED( spectrum_y ) )  THEN
     262          ALLOCATE( spectrum_y( 1:ny/2, 1:100, 1:10 ) )
     263          spectrum_y = 0.0_wp
     264       ENDIF
     265
     266       ALLOCATE( var_d(nzb:nzt+1) )
     267       var_d = 0.0_wp
     268    ENDIF
     269
     270    spectra_initialized = .TRUE.
     271
     272 END SUBROUTINE spectra_init
     273
     274
     275
     276!--------------------------------------------------------------------------------------------------!
    255277! Description:
    256278! ------------
    257279!> Check spectra related quantities
    258 !------------------------------------------------------------------------------!
    259     SUBROUTINE spectra_check_parameters
    260 
    261        USE control_parameters,                                                 &
    262            ONLY:  averaging_interval, message_string, skip_time_data_output
    263 
    264        IMPLICIT NONE
    265 
    266 !
    267 !--    Check the average interval
    268        IF ( averaging_interval_sp == 9999999.9_wp )  THEN
    269           averaging_interval_sp = averaging_interval
    270        ENDIF
    271 
    272        IF ( averaging_interval_sp > dt_dosp )  THEN
    273           WRITE( message_string, * )  'averaging_interval_sp = ',              &
    274                 averaging_interval_sp, ' must be <= dt_dosp = ', dt_dosp
    275           CALL message( 'spectra_check_parameters', 'PA0087', 1, 2, 0, 6, 0 )
    276        ENDIF
    277 
    278 !
    279 !--    Set the default skip time interval for data output, if necessary
    280        IF ( skip_time_dosp == 9999999.9_wp )                                   &
    281                                           skip_time_dosp = skip_time_data_output
    282 
    283     END SUBROUTINE spectra_check_parameters
    284 
    285 
    286 
    287 !------------------------------------------------------------------------------!
     280!--------------------------------------------------------------------------------------------------!
     281 SUBROUTINE spectra_check_parameters
     282
     283    USE control_parameters,                                                                        &
     284        ONLY:  averaging_interval,