Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

Changeset 3866 for palm


Ignore:
Timestamp:
Apr 5, 2019 2:25:01 PM (6 years ago)
Author:
eckhard
Message:

inifor: Use PALM's working precision; improved error handling, coding style, and comments

Location:
palm/trunk/UTIL
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • palm/trunk/UTIL/Makefile_utilities

    r3795 r3866  
    2323# -----------------
    2424# $Id$
     25# Use PALM's kinds module in inifor
     26#
     27#
     28# 3795 2019-03-15 09:40:05Z eckhard
    2529# Upated inifor build dependencies
    2630#
    … …  
    141145        inifor_defs.o \
    142146        inifor_util.o
     147inifor_defs.o: \
     148        mod_kinds.o
    143149inifor_grid.o: \
    144150        inifor_control.o \

  • palm/trunk/UTIL/inifor/src/inifor.f90

    r3785 r3866  
    2626! -----------------
    2727! $Id$
     28! Use PALM's working precision
     29! Show error message if compiled without netCDF support
     30! Renamed run_control -> log_runtime
     31! Improved coding style added comments
     32!
     33!
     34! 3785 2019-03-06 10:41:14Z eckhard
    2835! Average geostrophic wind components on coarse COSMO levels instead of fine PALM levels
    2936! Remove --debug netCDF output of internal pressure profiles
    … …  
    8188!------------------------------------------------------------------------------!
    8289 PROGRAM inifor
     90
    8391#if defined ( __netcdf )
    8492
    … …  
    8694    USE inifor_defs
    8795    USE inifor_grid,                                                           &
    88         ONLY:  setup_parameters, setup_grids, setup_variable_tables,           &
    89                setup_io_groups, fini_grids, fini_variables, fini_io_groups,    &
    90                fini_file_lists, preprocess, origin_lon, origin_lat,            &
    91                output_file, io_group_list, output_var_table,                   &
    92                cosmo_grid, palm_grid, nx, ny, nz, p0, cfg, f3,                 &
    93                averaging_width_ns, averaging_width_ew, phi_n, lambda_n,        &
    94                lam_centre, phi_centre
     96        ONLY:  averaging_width_ns,                                             &
     97               averaging_width_ew,                                             &
     98               cfg,                                                            &   
     99               cosmo_grid,                                                     &
     100               f3,                                                             &
     101               fini_grids,                                                     &
     102               fini_io_groups,                                                 &
     103               fini_variables,                                                 &
     104               fini_file_lists,                                                &
     105               io_group_list,                                                  &
     106               lam_centre,                                                     &
     107               lambda_n,                                                       &
     108               nx, ny, nz,                                                     &
     109               origin_lat,                                                     &
     110               origin_lon,                                                     &
     111               output_file,                                                    &
     112               output_var_table,                                               &
     113               p0,                                                             &
     114               phi_centre,                                                     &
     115               phi_n,                                                          &
     116               preprocess,                                                     &
     117               palm_grid,                                                      &
     118               setup_grids,                                                    &
     119               setup_parameters,                                               &
     120               setup_variable_tables,                                          &
     121               setup_io_groups
    95122    USE inifor_io
    96123    USE inifor_transform,                                                      &
    97         ONLY:  average_pressure_perturbation, average_profile, interpolate_1d, &
    98                interpolate_1d_arr, interpolate_2d, interpolate_3d,             &
    99                interp_average_profile, geostrophic_winds, extrapolate_density, &
    100                extrapolate_pressure, get_surface_pressure
     124        ONLY:  average_pressure_perturbation,                                  &
     125               average_profile,                                                &
     126               extrapolate_density,                                            &
     127               extrapolate_pressure,                                           &
     128               geostrophic_winds,                                              &
     129               get_surface_pressure,                                           &
     130               interp_average_profile,                                         &
     131               interpolate_1d,                                                 &
     132               interpolate_1d_arr,                                             &
     133               interpolate_2d,                                                 &
     134               interpolate_3d
    101135    USE inifor_types
    102136   
    … …  
    107141    INTEGER ::  iter   !< loop index for time steps
    108142
    109     REAL(dp), ALLOCATABLE, DIMENSION(:,:,:)     ::  output_arr !< array buffer for interpolated quantities
    110     REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_centre !< density profile of the centre averaging domain
    111     REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  ug_cosmo   !< profile of the geostrophic wind in x direction on COSMO levels
    112     REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  vg_cosmo   !< profile of the geostrophic wind in y direction on COSMO levels
    113     REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  ug_palm    !< profile of the geostrophic wind in x direction interpolated onto PALM levels
    114     REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  vg_palm    !< profile of the geostrophic wind in y direction interpolated onto PALM levels
    115     REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_north  !< density profile of the northern averaging domain
    116     REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_south  !< density profile of the southern averaging domain
    117     REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_east   !< density profile of the eastern averaging domain
    118     REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_west   !< density profile of the western averaging domain
    119     REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_north    !< pressure profile of the northern averaging domain
    120     REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_south    !< pressure profile of the southern averaging domain
    121     REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_east     !< pressure profile of the eastern averaging domain
    122     REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_west     !< pressure profile of the western averaging domain
    123 
    124     REAL(dp), POINTER, DIMENSION(:) ::  internal_arr !< pointer to the currently processed internal array (density, pressure)
    125     REAL(dp), POINTER, DIMENSION(:) ::  ug_vg_palm   !< pointer to the currently processed geostrophic wind component
     143    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)     ::  output_arr !< array buffer for interpolated quantities
     144    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_centre !< density profile of the centre averaging domain
     145    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  ug_cosmo   !< profile of the geostrophic wind in x direction on COSMO levels
     146    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  vg_cosmo   !< profile of the geostrophic wind in y direction on COSMO levels
     147    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  ug_palm    !< profile of the geostrophic wind in x direction interpolated onto PALM levels
     148    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  vg_palm    !< profile of the geostrophic wind in y direction interpolated onto PALM levels
     149    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_north  !< density profile of the northern averaging domain
     150    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_south  !< density profile of the southern averaging domain
     151    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_east   !< density profile of the eastern averaging domain
     152    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_west   !< density profile of the western averaging domain
     153    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_north    !< pressure profile of the northern averaging domain
     154    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_south    !< pressure profile of the southern averaging domain
     155    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_east     !< pressure profile of the eastern averaging domain
     156    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_west     !< pressure profile of the western averaging domain
     157
     158    REAL(wp), POINTER, DIMENSION(:) ::  internal_arr !< pointer to the currently processed internal array (density, pressure)
     159    REAL(wp), POINTER, DIMENSION(:) ::  ug_vg_palm   !< pointer to the currently processed geostrophic wind component
    126160
    127161    TYPE(nc_var), POINTER        ::  output_var      !< pointer to the currently processed output variable
    … …  
    138172!- Section 1: Initialization
    139173!------------------------------------------------------------------------------
    140  CALL run_control('init', 'void')
     174    CALL log_runtime( 'init', 'void' )
    141175
    142176!
    143177!-- Initialize INIFOR's parameters from command-line interface and namelists
    144     CALL setup_parameters()
     178    CALL setup_parameters
    145179
    146180!
    147181!-- Initialize all grids, including interpolation neighbours and weights
    148     CALL setup_grids()
    149  CALL run_control('time', 'init')
     182    CALL setup_grids
     183    CALL log_runtime( 'time', 'init' )
    150184
    151185!
    152186!-- Initialize the netCDF output file and define dimensions
    153     CALL setup_netcdf_dimensions(output_file, palm_grid, cfg % start_date,    &
     187    CALL setup_netcdf_dimensions(output_file, palm_grid, cfg%start_date,    &
    154188                                 origin_lon, origin_lat)
    155  CALL run_control('time', 'write')
     189    CALL log_runtime( 'time', 'write' )
    156190
    157191!
    158192!-- Set up the tables containing the input and output variables and set
    159193!-- the corresponding netCDF dimensions for each output variable
    160     CALL setup_variable_tables(cfg % ic_mode)
    161  CALL run_control('time', 'write')
     194    CALL setup_variable_tables( cfg%ic_mode )
     195    CALL log_runtime( 'time', 'write' )
    162196
    163197!
    164198!-- Add the output variables to the netCDF output file
    165     CALL setup_netcdf_variables(output_file % name, output_var_table)
    166 
    167     CALL setup_io_groups()
    168  CALL run_control('time', 'init')
    169 
    170 !------------------------------------------------------------------------------
    171 !- Section 2: Main loop
    172 !------------------------------------------------------------------------------
    173     DO igroup = 1, SIZE(io_group_list)
     199    CALL setup_netcdf_variables( output_file%name, output_var_table )
     200
     201    CALL setup_io_groups
     202    CALL log_runtime( 'time', 'init' )
     203
     204!------------------------------------------------------------------------------
     205!-- Section 2: Main loop
     206!------------------------------------------------------------------------------
     207!
     208!-- Input and output variables are organized into IO groups. For instance, the
     209!-- 'thermodynamics' IO group bundles the input variaebls T, P, QV and the
     210!-- output variables p, theta, rho, and qv.
     211!-- An IO group bunldes variables that are physically dependent on each other.
     212!-- In case of the 'thermodynamics' group, theta = f(P,T), rho = f(P,T,QV).
     213    DO  igroup = 1, SIZE( io_group_list )
    174214
    175215       group => io_group_list(igroup)
    176        IF ( group % to_be_processed )  THEN
     216       IF ( group%to_be_processed )  THEN
    177217         
    178           DO iter = 1, group % nt
    179 
    180 !------------------------------------------------------------------------------
    181 !- Section 2.1: Read and preprocess input data
    182 !------------------------------------------------------------------------------
    183              CALL read_input_variables(group, iter, input_buffer)
    184  CALL run_control('time', 'read')
    185 
    186              CALL preprocess(group, input_buffer, cosmo_grid, iter)
    187  CALL run_control('time', 'comp')
     218!--       Loop over all output time steps for the current group.
     219          DO  iter = 1, group%nt
     220
     221!------------------------------------------------------------------------------
     222!-- Section 2.1: Read and preprocess input data
     223!------------------------------------------------------------------------------
     224             CALL read_input_variables( group, iter, input_buffer )
     225             CALL log_runtime( 'time', 'read' )
     226
     227!--          Carry out all required physical conversion of the input variables
     228!--          of the current IO group on the input (COSMO) grid. For instance,
     229!--          horizontal velocities are rotated to the PALM coordinate system and
     230!--          potential temperature is computed from the absolute temperature and
     231!--          pressure.
     232             CALL preprocess( group, input_buffer, cosmo_grid, iter )
     233             CALL log_runtime( 'time', 'comp' )
    188234
    189235             !TODO: move this assertion into 'preprocess'.
    190              IF ( .NOT. ALL(input_buffer(:) % is_preprocessed .AND. .TRUE.) )  THEN
    191                 message = "Input buffers for group '" // TRIM(group % kind) // &
    192                    "' could not be preprocessed sucessfully."
    193                 CALL inifor_abort('main loop', message)
     236             IF ( .NOT. ALL(input_buffer(:)%is_preprocessed .AND. .TRUE.) )  THEN
     237                message = "Input buffers for group '" // TRIM( group%kind ) // &
     238                          "' could not be preprocessed sucessfully."
     239                CALL inifor_abort( 'main loop', message )
    194240             ENDIF
    195241
    196242!------------------------------------------------------------------------------
    197 !- Section 2.2: Interpolate each output variable of the group
    198 !------------------------------------------------------------------------------
    199              DO ivar = 1, group % nv
    200 
    201                 output_var => group % out_vars( ivar )
    202 
    203                 IF ( output_var % to_be_processed .AND.                        &
    204                      iter .LE. output_var % nt )  THEN
    205 
    206                    message = "Processing '" // TRIM(output_var % name) //      &
    207                              "' (" // TRIM(output_var % kind) //               &
    208                              "), iteration " // TRIM(str(iter)) //" of " //    &
    209                              TRIM(str(output_var % nt))
    210                    CALL report('main loop', message)
    211 
    212                    SELECT CASE( TRIM(output_var % task) )
    213 
    214                    CASE( 'interpolate_2d' )
    215                    
    216                       SELECT CASE( TRIM(output_var % kind) )
    217                        
    218                       CASE( 'init soil' )
    219 
    220                          ALLOCATE( output_arr( 0:output_var % grid % nx,       &
    221                                                0:output_var % grid % ny,       &
    222                                                SIZE(output_var % grid % depths) ) )
    223 
    224                       CASE ( 'surface forcing' )
    225 
    226                          ALLOCATE( output_arr( 0:output_var % grid % nx,       &
    227                                                0:output_var % grid % ny, 1 ) )
    228 
    229                       CASE DEFAULT
    230 
    231                           message = "'" // TRIM(output_var % kind) // "' is not a soil variable"
    232                           CALL inifor_abort("main loop", message)
    233 
    234                       END SELECT
    235  CALL run_control('time', 'alloc')
    236 
    237                       CALL interpolate_2d(input_buffer(output_var % input_id) % array(:,:,:), &
    238                               output_arr(:,:,:), output_var % intermediate_grid, output_var)
    239  CALL run_control('time', 'comp')
    240 
    241 
    242                    CASE ( 'interpolate_3d' )
    243 
    244                       ALLOCATE( output_arr( 0:output_var % grid % nx,          &
    245                                             0:output_var % grid % ny,          &
    246                                             1:output_var % grid % nz ) )
    247 
    248  CALL run_control('time', 'alloc')
    249                       CALL interpolate_3d(                                     &
    250                          input_buffer(output_var % input_id) % array(:,:,:),   &
    251                          output_arr(:,:,:),                                    &
    252                          output_var % intermediate_grid,                       &
    253                          output_var % grid)
    254  CALL run_control('time', 'comp')
    255 
    256                    CASE ( 'average profile' )
    257 
    258                       ALLOCATE( output_arr( 0:output_var % grid % nx,          &
    259                                             0:output_var % grid % ny,          &
    260                                             1:output_var % grid % nz ) )
    261  CALL run_control('time', 'alloc')
     243!-- Section 2.2: Interpolate each output variable of the group
     244!------------------------------------------------------------------------------
     245             DO  ivar = 1, group%nv
     246
     247                output_var => group%out_vars(ivar)
     248
     249                IF ( output_var%to_be_processed .AND.                          &
     250                     iter .LE. output_var%nt )  THEN
     251
     252                   message = "Processing '" // TRIM( output_var%name ) //      &
     253                             "' (" // TRIM( output_var%kind ) //               &
     254                             "), iteration " // TRIM( str( iter ) ) //" of " //&
     255                             TRIM( str( output_var%nt ) )
     256                   CALL report( 'main loop', message )
     257
     258                   SELECT CASE( TRIM( output_var%task ) )
     259
     260!--                   2D horizontal interpolation
     261                      CASE( 'interpolate_2d' )
    262262                     
    263                       CALL interp_average_profile(                             &
    264                          input_buffer(output_var % input_id) % array(:,:,:),   &
    265                          output_arr(0,0,:),                                    &
    266                          output_var % averaging_grid)
    267 
    268                       IF ( TRIM(output_var % name) ==                          &
    269                            'surface_forcing_surface_pressure' )  THEN
    270 
    271                          IF ( cfg % p0_is_set )  THEN
    272                             output_arr(0,0,1) = p0
    273                          ELSE
    274                             CALL get_surface_pressure(                         &
    275                                output_arr(0,0,:), rho_centre,                  &
    276                                output_var % averaging_grid)
     263                         SELECT CASE( TRIM( output_var%kind ) )
     264                         
     265                         CASE( 'init soil' )
     266   
     267                            ALLOCATE( output_arr(0:output_var%grid%nx,            &
     268                                                 0:output_var%grid%ny,            &
     269                                                 SIZE( output_var%grid%depths )) )
     270   
     271                         CASE ( 'surface forcing' )
     272   
     273                            ALLOCATE( output_arr(0:output_var%grid%nx,            &
     274                                                 0:output_var%grid%ny, 1) )
     275   
     276                         CASE DEFAULT
     277   
     278                             message = "'" // TRIM( output_var%kind ) // "' is not a soil variable"
     279                             CALL inifor_abort( "main loop", message )
     280   
     281                         END SELECT
     282                         CALL log_runtime( 'time', 'alloc' )
     283   
     284                         CALL interpolate_2d( input_buffer(output_var%input_id)%array(:,:,:), &
     285                                 output_arr(:,:,:), output_var%intermediate_grid, output_var )
     286                         CALL log_runtime( 'time', 'comp' )
     287   
     288   
     289!--                   Interpolation in 3D, used for atmospheric initial and
     290!--                   boundary conditions.
     291                      CASE ( 'interpolate_3d' )
     292   
     293                         ALLOCATE( output_arr(0:output_var%grid % nx,           &
     294                                              0:output_var%grid % ny,           &
     295                                              1:output_var%grid % nz) )
     296   
     297                         CALL log_runtime( 'time', 'alloc' )
     298                         CALL interpolate_3d(                                     &
     299                            input_buffer(output_var%input_id)%array(:,:,:),       &
     300                            output_arr(:,:,:),                                    &
     301                            output_var%intermediate_grid,                         &
     302                            output_var%grid)
     303                         CALL log_runtime( 'time', 'comp' )
     304   
     305!--                   Compute initial avaerage profiles (if --init-mode profile
     306!--                   is used)
     307                      CASE ( 'average profile' )
     308   
     309                         ALLOCATE( output_arr(0:output_var%grid%nx,               &
     310                                              0:output_var%grid%ny,               &
     311                                              1:output_var%grid%nz) )
     312                         CALL log_runtime( 'time', 'alloc' )
     313                         
     314                         CALL interp_average_profile(                             &
     315                            input_buffer(output_var%input_id)%array(:,:,:),     &
     316                            output_arr(0,0,:),                                    &
     317                            output_var%averaging_grid )
     318   
     319                         IF ( TRIM( output_var%name ) ==                          &
     320                              'surface_forcing_surface_pressure' )  THEN
     321   
     322                            IF ( cfg%p0_is_set )  THEN
     323                               output_arr(0,0,1) = p0
     324                            ELSE
     325                               CALL get_surface_pressure(                         &
     326                                  output_arr(0,0,:), rho_centre,                  &
     327                                  output_var%averaging_grid )
     328                            ENDIF
     329   
    277330                         ENDIF
    278 
    279                       ENDIF
    280  CALL run_control('time', 'comp')
    281 
    282                    CASE ( 'internal profile' )
    283 
    284                       message = "Averaging of internal profile for variable '" //&
    285                          TRIM(output_var % name) // "' is not supported."
    286 
    287                       SELECT CASE (TRIM(output_var % name))
    288 
    289                       CASE('internal_density_centre')
    290                          ALLOCATE( rho_centre( 1:cosmo_grid % nz) )
    291                          internal_arr => rho_centre
    292 
    293                       CASE('internal_density_north')
    294                          ALLOCATE( rho_north( 1:cosmo_grid % nz) )
    295                          internal_arr => rho_north
    296 
    297                       CASE('internal_density_south')
    298                          ALLOCATE( rho_south( 1:cosmo_grid % nz) )
    299                          internal_arr => rho_south
    300 
    301                       CASE('internal_density_east')
    302                          ALLOCATE( rho_east( 1:cosmo_grid % nz) )
    303                          internal_arr => rho_east
    304 
    305                       CASE('internal_density_west')
    306                          ALLOCATE( rho_west( 1:cosmo_grid % nz) )
    307                          internal_arr => rho_west
    308 
    309                       CASE('internal_pressure_north')
    310                          ALLOCATE( p_north( 1:cosmo_grid % nz) )
    311                          internal_arr => p_north
    312 
    313                       CASE('internal_pressure_south')
    314                          ALLOCATE( p_south( 1:cosmo_grid % nz) )
    315                          internal_arr => p_south
    316 
    317                       CASE('internal_pressure_east')
    318                          ALLOCATE( p_east( 1:cosmo_grid % nz) )
    319                          internal_arr => p_east
    320 
    321                       CASE('internal_pressure_west')
    322                          ALLOCATE( p_west( 1:cosmo_grid % nz) )
    323                          internal_arr => p_west
    324 
    325                       CASE DEFAULT
    326                          CALL inifor_abort('main loop', message)
    327 
    328                       END SELECT
    329  CALL run_control('time', 'alloc')
    330 
    331 
    332                       SELECT CASE( TRIM( output_var % name ) )
    333 
    334                       CASE( 'internal_pressure_north',                         &
    335                             'internal_pressure_south',                         &
    336                             'internal_pressure_east',                          &
    337                             'internal_pressure_west' )
    338 
    339                          CALL average_pressure_perturbation(                   &
    340                             input_buffer(output_var % input_id) % array(:,:,:),&
    341                             internal_arr(:),                                   &
    342                             cosmo_grid, output_var % averaging_grid            &
    343                          )
    344 
    345                       CASE DEFAULT
    346 
    347                          CALL average_profile(                                 &
    348                             input_buffer(output_var % input_id) % array(:,:,:),&
    349                             internal_arr(:),                                   &
    350                             output_var % averaging_grid                        &
    351                          )
     331                         CALL log_runtime( 'time', 'comp' )
     332   
     333!--                   Compute internal profiles, required for differentiation of
     334!--                   geostrophic wind
     335                      CASE ( 'internal profile' )
     336   
     337                         message = "Averaging of internal profile for variable '" //&
     338                            TRIM( output_var%name ) // "' is not supported."
     339   
     340                         SELECT CASE ( TRIM( output_var%name ) )
     341   
     342                         CASE( 'internal_density_centre' )
     343                            ALLOCATE( rho_centre(1:cosmo_grid%nz) )
     344                            internal_arr => rho_centre
     345   
     346                         CASE( 'internal_density_north' )
     347                            ALLOCATE( rho_north(1:cosmo_grid%nz) )
     348                            internal_arr => rho_north
     349   
     350                         CASE( 'internal_density_south' )
     351                            ALLOCATE( rho_south(1:cosmo_grid%nz) )
     352                            internal_arr => rho_south
     353   
     354                         CASE( 'internal_density_east' )
     355                            ALLOCATE( rho_east(1:cosmo_grid%nz) )
     356                            internal_arr => rho_east
     357   
     358                         CASE( 'internal_density_west' )
     359                            ALLOCATE( rho_west(1:cosmo_grid%nz) )
     360                            internal_arr => rho_west
     361   
     362                         CASE( 'internal_pressure_north' )
     363                            ALLOCATE( p_north(1:cosmo_grid%nz) )
     364                            internal_arr => p_north
     365   
     366                         CASE( 'internal_pressure_south' )
     367                            ALLOCATE( p_south(1:cosmo_grid%nz) )
     368                            internal_arr => p_south
     369   
     370                         CASE( 'internal_pressure_east' )
     371                            ALLOCATE( p_east(1:cosmo_grid%nz) )
     372                            internal_arr => p_east
     373   
     374                         CASE( 'internal_pressure_west' )
     375                            ALLOCATE( p_west(1:cosmo_grid%nz) )
     376                            internal_arr => p_west
     377   
     378                         CASE DEFAULT
     379                            CALL inifor_abort( 'main loop', message )
     380   
     381                         END SELECT
     382                         CALL log_runtime( 'time', 'alloc' )
     383   
     384   
     385                         SELECT CASE( TRIM( output_var%name ) )
     386   
     387                         CASE( 'internal_pressure_north',                         &
     388                               'internal_pressure_south',                         &
     389                               'internal_pressure_east',                          &
     390                               'internal_pressure_west' )
     391   
     392                            CALL average_pressure_perturbation(                   &
     393                               input_buffer(output_var%input_id) % array(:,:,:),&
     394                               internal_arr(:),                                   &
     395                               cosmo_grid, output_var%averaging_grid            &
     396                            )
     397   
     398                         CASE DEFAULT
     399   
     400                            CALL average_profile(                                 &
     401                               input_buffer(output_var%input_id) % array(:,:,:),&
     402                               internal_arr(:),                                   &
     403                               output_var%averaging_grid                        &
     404                            )
    352405
    353406                      END SELECT
    … …  
    360413!--                   (requires definiton of COSMO levels in netCDF output.)
    361414                      !IF (.TRUE.)  THEN
    362                       !   ALLOCATE( output_arr(1,1,1:output_var % grid % nz) )
     415                      !   ALLOCATE( output_arr(1,1,1:output_var%grid % nz) )
    363416                      !   output_arr(1,1,:) = internal_arr(:)
    364417                      !ENDIF
    365  CALL run_control('time', 'comp')
     418                      CALL log_runtime( 'time', 'comp' )
    366419
    367420!
    … …  
    372425
    373426                      IF (.NOT. ug_vg_have_been_computed )  THEN
    374                          ALLOCATE( ug_palm(output_var % grid % nz) )
    375                          ALLOCATE( vg_palm(output_var % grid % nz) )
    376                          ALLOCATE( ug_cosmo(cosmo_grid % nz) )
    377                          ALLOCATE( vg_cosmo(cosmo_grid % nz) )
    378 
    379                          IF ( cfg % ug_defined_by_user )  THEN
    380                             ug_palm = cfg % ug
    381                             vg_palm = cfg % vg
     427                         ALLOCATE( ug_palm(output_var%grid%nz) )
     428                         ALLOCATE( vg_palm(output_var%grid%nz) )
     429                         ALLOCATE( ug_cosmo(cosmo_grid%nz) )
     430                         ALLOCATE( vg_cosmo(cosmo_grid%nz) )
     431
     432                         IF ( cfg%ug_defined_by_user )  THEN
     433                            ug_palm = cfg%ug
     434                            vg_palm = cfg%vg
    382435                         ELSE
    383436                            CALL geostrophic_winds( p_north, p_south, p_east,  &
    … …  
    390443
    391444                            CALL interpolate_1d( ug_cosmo, ug_palm,             &
    392                                                  output_var % grid )
     445                                                 output_var%grid )
    393446
    394447                            CALL interpolate_1d( vg_cosmo, vg_palm,             &
    395                                                  output_var % grid )
     448                                                 output_var%grid )
    396449                         ENDIF
    397450
    … …  
    402455!
    403456!--                   Select output array of current geostrophic wind component
    404                       SELECT CASE(TRIM(output_var % name))
    405                       CASE ('ls_forcing_ug')
     457                      SELECT CASE( TRIM( output_var%name ) )
     458                      CASE ( 'ls_forcing_ug' )
    406459                         ug_vg_palm => ug_palm
    407                       CASE ('ls_forcing_vg')
     460                      CASE ( 'ls_forcing_vg' )
    408461                         ug_vg_palm => vg_palm
    409462                      END SELECT
    410463
    411                       ALLOCATE( output_arr(1,1,output_var % grid % nz) )
     464                      ALLOCATE( output_arr(1,1,output_var%grid%nz) )
    412465                      output_arr(1,1,:) = ug_vg_palm(:)
    413466
    414                    CASE ( 'average scalar' )
    415 
    416                       ALLOCATE( output_arr(1,1,1) )
    417  CALL run_control('time', 'alloc')
    418                       output_arr(1,1,1) = p0
    419  CALL run_control('time', 'comp')
    420 
     467!--                User defined constant profiles
    421468                   CASE ( 'set profile' )
    422469                     
    423                       ALLOCATE( output_arr( 1, 1, 1:nz ) )
    424  CALL run_control('time', 'alloc')
    425 
    426                       SELECT CASE (TRIM(output_var % name))
    427 
    428                       CASE('nudging_tau')
     470                      ALLOCATE( output_arr(1,1,1:nz) )
     471                      CALL log_runtime( 'time', 'alloc' )
     472
     473                      SELECT CASE ( TRIM( output_var%name ) )
     474
     475                      CASE ( 'nudging_tau' )
    429476                          output_arr(1, 1, :) = NUDGING_TAU
    430477
    431478                      CASE DEFAULT
    432                           message = "'" // TRIM(output_var % name) //          &
    433                              "' is not a valid '" // TRIM(output_var % kind) //&
     479                          message = "'" // TRIM( output_var%name ) //          &
     480                             "' is not a valid '" // TRIM( output_var%kind ) //&
    434481                             "' variable kind."
    435                           CALL inifor_abort('main loop', message)
     482                          CALL inifor_abort( 'main loop', message )
    436483                      END SELECT
    437  CALL run_control('time', 'comp')
    438 
    439                    CASE('average large-scale profile')
    440                       message = "Averaging of large-scale forcing profiles " //&
    441                                 "has not been implemented, yet."
    442                       CALL inifor_abort('main loop', message)
     484                      CALL log_runtime( 'time', 'comp' )
    443485
    444486                   CASE DEFAULT
    445                       message = "Processing task '" // TRIM(output_var % task) //&
     487                      message = "Processing task '" // TRIM( output_var%task ) //&
    446488                               "' not recognized."
    447                       CALL inifor_abort('', message)
     489                      CALL inifor_abort( '', message )
    448490
    449491                   END SELECT
    450  CALL run_control('time', 'comp')
     492                   CALL log_runtime( 'time', 'comp' )
    451493
    452494!------------------------------------------------------------------------------
    … …  
    458500!--                defined on averaged COSMO levels instead of PALM levels
    459501!--                (requires definiton of COSMO levels in netCDF output.)
    460                    !IF (.NOT. output_var % is_internal .OR. debugging_output)  THEN
    461 
    462                    IF (.NOT. output_var % is_internal)  THEN
    463                       message = "Writing variable '" // TRIM(output_var%name) // "'."
    464                       CALL report('main loop', message)
    465                       CALL update_output(output_var, output_arr, iter,         &
    466                                          output_file, cfg)
    467  CALL run_control('time', 'write')
     502                   !IF (.NOT. output_var%is_internal .OR. debugging_output)  THEN
     503
     504                   IF ( .NOT. output_var%is_internal )  THEN
     505                      message = "Writing variable '" // TRIM( output_var%name ) // "'."
     506                      CALL report( 'main loop', message )
     507                      CALL update_output( output_var, output_arr, iter,        &
     508                                          output_file, cfg )
     509                      CALL log_runtime( 'time', 'write' )
    468510                   ENDIF
    469511
    470                    IF (ALLOCATED(output_arr))  DEALLOCATE(output_arr)
    471  CALL run_control('time', 'alloc')
     512                   IF ( ALLOCATED( output_arr ) )  DEALLOCATE( output_arr )
     513                   CALL log_runtime( 'time', 'alloc' )
    472514
    473515                ENDIF
    … …  
    478520
    479521             ug_vg_have_been_computed = .FALSE.
    480              IF ( group % kind == 'thermodynamics' )  THEN
     522             IF ( group%kind == 'thermodynamics' )  THEN
    481523                DEALLOCATE( rho_centre )
    482524                DEALLOCATE( ug_palm )
    … …  
    484526                DEALLOCATE( ug_cosmo )
    485527                DEALLOCATE( vg_cosmo )
    486                 IF ( .NOT. cfg % ug_defined_by_user )  THEN
     528                IF ( .NOT. cfg%ug_defined_by_user )  THEN
    487529                   DEALLOCATE( rho_north )
    488530                   DEALLOCATE( rho_south )
    … …  
    499541!--          Keep input buffer around for averaged (radiation) and
    500542!--          accumulated COSMO quantities (precipitation).
    501              IF ( group % kind == 'running average' .OR. &
    502                   group % kind == 'accumulated' )  THEN
     543             IF ( group%kind == 'running average' .OR. &
     544                  group%kind == 'accumulated' )  THEN
    503545             ELSE
    504                 CALL report('main loop', 'Deallocating input buffer', cfg % debug)
    505                 DEALLOCATE(input_buffer)
     546                CALL report( 'main loop', 'Deallocating input buffer', cfg%debug )
     547                DEALLOCATE( input_buffer )
    506548             ENDIF
    507  CALL run_control('time', 'alloc')
     549             CALL log_runtime( 'time', 'alloc' )
    508550
    509551!
    … …  
    511553          ENDDO
    512554
    513           IF (ALLOCATED(input_buffer))  THEN
    514              CALL report('main loop', 'Deallocating input buffer', cfg % debug)
    515              DEALLOCATE(input_buffer)
     555          IF ( ALLOCATED( input_buffer ) )  THEN
     556             CALL report( 'main loop', 'Deallocating input buffer', cfg%debug )
     557             DEALLOCATE( input_buffer )
    516558          ENDIF
    517  CALL run_control('time', 'alloc')
     559          CALL log_runtime( 'time', 'alloc' )
    518560
    519561       ELSE
    520562
    521           message = "Skipping IO group " // TRIM(str(igroup)) // " '" // TRIM(group % kind) // "'"
    522           IF ( ALLOCATED(group % in_var_list) )  THEN
    523               message = TRIM(message) // " with input variable '" //           &
    524               TRIM(group % in_var_list(1) % name) // "'."
     563          message = "Skipping IO group " // TRIM( str( igroup ) ) // " '" // TRIM( group%kind ) // "'"
     564          IF ( ALLOCATED( group%in_var_list ) )  THEN
     565              message = TRIM( message ) // " with input variable '" //         &
     566              TRIM( group%in_var_list(1)%name ) // "'."
    525567          ENDIF
    526568
    527           CALL report('main loop', message, cfg % debug)
    528 
    529 !
    530 !--    IO group % to_be_processed conditional
     569          CALL report( 'main loop', message, cfg%debug )
     570
     571!
     572!--    IO group%to_be_processed conditional
    531573       ENDIF
    532574
    … …  
    538580!- Section 3: Clean up.
    539581!------------------------------------------------------------------------------
    540     CALL fini_file_lists()
    541     CALL fini_io_groups()
    542     CALL fini_variables()
    543     !CALL fini_grids()
    544  CALL run_control('time', 'alloc')
    545  CALL run_control('report', 'void')
    546 
    547     message = "Finished writing dynamic driver '" // TRIM(output_file % name) // &
     582    CALL fini_file_lists
     583    CALL fini_io_groups
     584    CALL fini_variables
     585    !CALL fini_grids
     586    CALL log_runtime( 'time', 'alloc' )
     587    CALL log_runtime( 'report', 'void' )
     588
     589    message = "Finished writing dynamic driver '" // TRIM( output_file%name ) // &
    548590              "' successfully."
    549     CALL report('main loop', message)
    550 
    551 
     591    CALL report( 'main loop', message )
     592    CALL close_log
     593
     594#else
     595
     596    USE inifor_control
     597    IMPLICIT NONE
     598   
     599    message = "INIFOR was compiled without netCDF support, which is required for it to run. "  //     &
     600              "To use INIFOR, recompile PALM with netCDF support by adding the -D__netcdf " //        &
     601              "precompiler flag to your .palm.config file."
     602    CALL inifor_abort( 'main loop', message )
     603 
    552604#endif
     605
    553606 END PROGRAM inifor

  • palm/trunk/UTIL/inifor/src/inifor_control.f90

    r3785 r3866  
    2626! -----------------
    2727! $Id$
     28! Use PALM's working precision
     29! Renamed run_control -> log_runtime
     30! Open log file only once
     31! Improved coding style
     32!
     33!
     34! 3785 2019-03-06 10:41:14Z eckhard
    2835! Added message buffer for displaying tips to rectify encountered errors
    2936!
    … …  
    6673!> feedback to the terminal and a log file.
    6774!------------------------------------------------------------------------------!
    68 #if defined ( __netcdf )
    6975 MODULE inifor_control
    7076
    7177    USE inifor_defs,                                                           &
    72         ONLY:  LNAME, dp, VERSION, COPYRIGHT
    73 
     78        ONLY:  COPYRIGHT, LNAME, LOG_FILE_NAME, VERSION, wp
    7479    USE inifor_util,                                                           &
    7580        ONLY:  real_to_str, real_to_str_f
    … …  
    7984    CHARACTER (LEN=5000) ::  message = '' !< log message buffer
    8085    CHARACTER (LEN=5000) ::  tip     = '' !< optional log message buffer for tips on how to rectify encountered errors
     86    INTEGER, SAVE        ::  u            !< Fortran file unit for the log file
    8187
    8288 CONTAINS
    … …  
    94100!> to it.
    95101!------------------------------------------------------------------------------!
    96     SUBROUTINE report(routine, message, debug)
    97 
    98        CHARACTER(LEN=*), INTENT(IN)  ::  routine !< name of calling subroutine of function
    99        CHARACTER(LEN=*), INTENT(IN)  ::  message !< log message
    100        LOGICAL, OPTIONAL, INTENT(IN) ::  debug   !< flag the current message as debugging message
    101 
    102        INTEGER                       ::  u                     !< Fortran file unit for the log file
    103        LOGICAL, SAVE                 ::  is_first_run = .TRUE. !< control flag for file opening mode
    104        LOGICAL                       ::  suppress_message      !< control falg for additional debugging log
    105 
    106 
    107        IF ( is_first_run )  THEN
    108           OPEN( NEWUNIT=u, FILE='inifor.log', STATUS='replace' )
    109           is_first_run = .FALSE.
    110        ELSE
    111           OPEN( NEWUNIT=u, FILE='inifor.log', POSITION='append', STATUS='old' )
    112        ENDIF
    113          
    114 
    115        suppress_message = .FALSE.
    116        IF ( PRESENT(debug) )  THEN
    117           IF ( .NOT. debug )  suppress_message = .TRUE.
    118        ENDIF
    119 
    120        IF ( .NOT. suppress_message )  THEN
    121           PRINT *, "inifor: " // TRIM(message) // "  [ " // TRIM(routine) // " ]"
    122           WRITE(u, *)  TRIM(message) // "  [ " // TRIM(routine) // " ]"
    123        ENDIF
    124 
    125        CLOSE(u)
    126 
    127     END SUBROUTINE report
     102 SUBROUTINE report(routine, message, debug)
     103
     104    CHARACTER(LEN=*), INTENT(IN)  ::  routine !< name of calling subroutine of function
     105    CHARACTER(LEN=*), INTENT(IN)  ::  message !< log message
     106    LOGICAL, OPTIONAL, INTENT(IN) ::  debug   !< flag the current message as debugging message
     107
     108    LOGICAL, SAVE                 ::  is_first_run = .TRUE. !< control flag for file opening mode
     109    LOGICAL                       ::  suppress_message      !< control falg for additional debugging log
     110
     111    IF ( is_first_run )  THEN
     112       OPEN( NEWUNIT=u, FILE=LOG_FILE_NAME, STATUS='replace' )
     113       is_first_run = .FALSE.
     114    ENDIF
     115       
     116
     117    suppress_message = .FALSE.
     118    IF ( PRESENT(debug) )  THEN
     119       IF ( .NOT. debug )  suppress_message = .TRUE.
     120    ENDIF
     121
     122    IF ( .NOT. suppress_message )  THEN
     123       PRINT *, "inifor: " // TRIM(message) // "  [ " // TRIM(routine) // " ]"
     124       WRITE(u, *)  TRIM(message) // "  [ " // TRIM(routine) // " ]"
     125    ENDIF
     126
     127 END SUBROUTINE report
    128128
    129129
    … …  
    138138!> continue.
    139139!------------------------------------------------------------------------------!
    140     SUBROUTINE warn(routine, message)
    141 
    142        CHARACTER(LEN=*), INTENT(IN) ::  routine !< name of calling subroutine or function
    143        CHARACTER(LEN=*), INTENT(IN) ::  message !< log message
    144 
    145        CALL report(routine, "WARNING: " // TRIM(message))
    146 
    147     END SUBROUTINE warn
     140 SUBROUTINE warn(routine, message)
     141
     142    CHARACTER(LEN=*), INTENT(IN) ::  routine !< name of calling subroutine or function
     143    CHARACTER(LEN=*), INTENT(IN) ::  message !< log message
     144
     145    CALL report(routine, "WARNING: " // TRIM(message))
     146
     147 END SUBROUTINE warn
    148148
    149149
    … …  
    158158!> INIFOR from continueing.
    159159!------------------------------------------------------------------------------!
    160     SUBROUTINE inifor_abort(routine, message)
    161 
    162        CHARACTER(LEN=*), INTENT(IN) ::  routine !< name of calling subroutine or function
    163        CHARACTER(LEN=*), INTENT(IN) ::  message !< log message
    164 
    165        CALL report(routine, "ERROR: " // TRIM(message) // " Stopping.")
    166        STOP
    167 
    168     END SUBROUTINE inifor_abort
     160 SUBROUTINE inifor_abort(routine, message)
     161
     162    CHARACTER(LEN=*), INTENT(IN) ::  routine !< name of calling subroutine or function
     163    CHARACTER(LEN=*), INTENT(IN) ::  message !< log message
     164
     165    CALL report(routine, "ERROR: " // TRIM(message) // " Stopping.")
     166    CALL close_log
     167    STOP
     168
     169 END SUBROUTINE inifor_abort
     170
     171
     172 SUBROUTINE close_log()
     173
     174    CLOSE(u)
     175
     176 END SUBROUTINE close_log
    169177
    170178
    … …  
    175183!> print_version() prints the INIFOR version number and copyright notice.
    176184!------------------------------------------------------------------------------!
    177     SUBROUTINE print_version()
    178        PRINT *, "INIFOR " // VERSION
    179        PRINT *, COPYRIGHT
    180     END SUBROUTINE print_version
    181 
    182 
    183 !------------------------------------------------------------------------------!
    184 ! Description:
    185 ! ------------
    186 !>
    187 !> run_control() measures the run times of various parts of INIFOR and
     185 SUBROUTINE print_version()
     186    PRINT *, "INIFOR " // VERSION
     187    PRINT *, COPYRIGHT
     188 END SUBROUTINE print_version
     189
     190
     191!------------------------------------------------------------------------------!
     192! Description:
     193! ------------
     194!>
     195!> log_runtime() measures the run times of various parts of INIFOR and
    188196!> accumulates them in timing budgets.
    189197!------------------------------------------------------------------------------!
    190     SUBROUTINE run_control(mode, budget)
    191 
    192        CHARACTER(LEN=*), INTENT(IN) ::  mode   !< name of the calling mode
    193        CHARACTER(LEN=*), INTENT(IN) ::  budget !< name of the timing budget
    194 
    195        REAL(dp), SAVE ::  t0               !< begin of timing interval
    196        REAL(dp), SAVE ::  t1               !< end of timing interval
    197        REAL(dp), SAVE ::  t_comp  = 0.0_dp !< computation timing budget
    198        REAL(dp), SAVE ::  t_alloc = 0.0_dp !< allocation timing budget
    199        REAL(dp), SAVE ::  t_init  = 0.0_dp !< initialization timing budget
    200        REAL(dp), SAVE ::  t_read  = 0.0_dp !< reading timing budget
    201        REAL(dp), SAVE ::  t_total = 0.0_dp !< total time
    202        REAL(dp), SAVE ::  t_write = 0.0_dp !< writing timing budget
    203 
    204        CHARACTER(LEN=*), PARAMETER  ::  fmt='(F6.2)' !< floating-point output format
    205 
    206 
    207        SELECT CASE(TRIM(mode))
    208 
    209        CASE('init')
    210           CALL CPU_TIME(t0)
    211 
    212        CASE('time')
    213 
    214           CALL CPU_TIME(t1)
    215 
    216           SELECT CASE(TRIM(budget))
    217 
    218              CASE('alloc')
    219                 t_alloc = t_alloc + t1 - t0
    220 
    221              CASE('init')
    222                 t_init = t_init + t1 - t0
    223 
    224              CASE('read')
    225                 t_read = t_read + t1 - t0
    226 
    227              CASE('write')
    228                 t_write = t_write + t1 - t0
    229 
    230              CASE('comp')
    231                 t_comp = t_comp + t1 - t0
    232 
    233              CASE DEFAULT
    234                 CALL inifor_abort('run_control', "Time Budget '" // TRIM(mode) // "' is not supported.")
    235 
    236           END SELECT
    237 
    238           t0 = t1
    239 
    240        CASE('report')
    241            t_total = t_init + t_read + t_write + t_comp
    242 
    243            CALL report('run_control', " *** CPU time ***")
    244 
    245            CALL report('run_control', "Initialization: " // real_to_str(t_init)  // &
    246                        " s (" // TRIM(real_to_str(100*t_init/t_total, fmt))      // " %)")
    247 
    248            CALL report('run_control', "(De-)Allocation:" // real_to_str(t_alloc)  // &
    249                        " s (" // TRIM(real_to_str(100*t_alloc/t_total, fmt))      // " %)")
    250 
    251            CALL report('run_control', "Reading data:   " // real_to_str(t_read)  // &
    252                        " s (" // TRIM(real_to_str(100*t_read/t_total, fmt))      // " %)")
    253 
    254            CALL report('run_control', "Writing data:   " // real_to_str(t_write) // &
    255                        " s (" // TRIM(real_to_str(100*t_write/t_total, fmt))     // " %)")
    256 
    257            CALL report('run_control', "Computation:    " // real_to_str(t_comp)  // &
    258                        " s (" // TRIM(real_to_str(100*t_comp/t_total, fmt))      // " %)")
    259 
    260            CALL report('run_control', "Total:          " // real_to_str(t_total) // &
    261                        " s (" // TRIM(real_to_str(100*t_total/t_total, fmt))     // " %)")
    262 
    263        CASE DEFAULT
    264           CALL inifor_abort('run_control', "Mode '" // TRIM(mode) // "' is not supported.")
     198 SUBROUTINE log_runtime(mode, budget)
     199
     200    CHARACTER(LEN=*), INTENT(IN) ::  mode   !< name of the calling mode
     201    CHARACTER(LEN=*), INTENT(IN) ::  budget !< name of the timing budget
     202
     203    REAL(wp), SAVE ::  t0               !< begin of timing interval
     204    REAL(wp), SAVE ::  t1               !< end of timing interval
     205    REAL(wp), SAVE ::  t_comp  = 0.0_wp !< computation timing budget
     206    REAL(wp), SAVE ::  t_alloc = 0.0_wp !< allocation timing budget
     207    REAL(wp), SAVE ::  t_init  = 0.0_wp !< initialization timing budget
     208    REAL(wp), SAVE ::  t_read  = 0.0_wp !< reading timing budget
     209    REAL(wp), SAVE ::  t_total = 0.0_wp !< total time
     210    REAL(wp), SAVE ::  t_write = 0.0_wp !< writing timing budget
     211
     212    CHARACTER(LEN=*), PARAMETER  ::  fmt='(F6.2)' !< floating-point output format
     213
     214
     215    SELECT CASE(TRIM(mode))
     216
     217    CASE('init')
     218       CALL CPU_TIME(t0)
     219
     220    CASE('time')
     221
     222       CALL CPU_TIME(t1)
     223
     224       SELECT CASE(TRIM(budget))
     225
     226          CASE('alloc')
     227             t_alloc = t_alloc + t1 - t0
     228
     229          CASE('init')
     230             t_init = t_init + t1 - t0
     231
     232          CASE('read')
     233             t_read = t_read + t1 - t0
     234
     235          CASE('write')
     236             t_write = t_write + t1 - t0
     237
     238          CASE('comp')
     239             t_comp = t_comp + t1 - t0
     240
     241          CASE DEFAULT
     242             CALL inifor_abort('log_runtime', "Time Budget '" // TRIM(mode) // "' is not supported.")
    265243
    266244       END SELECT
    267245
    268     END SUBROUTINE run_control
     246       t0 = t1
     247
     248    CASE('report')
     249        t_total = t_init + t_read + t_write + t_comp
     250
     251        CALL report('log_runtime', " *** CPU time ***")
     252
     253        CALL report('log_runtime', "Initialization:  " // TRIM( real_to_str( t_init ) ) // &
     254                    " s  (" // TRIM( real_to_str( 100 * t_init / t_total, fmt ) ) // " %)" )
     255
     256        CALL report('log_runtime', "(De-)Allocation: " // TRIM( real_to_str( t_alloc ) ) // &
     257                    " s  (" // TRIM( real_to_str( 100 * t_alloc / t_total, fmt ) ) // " %)" )
     258
     259        CALL report('log_runtime', "Reading data:    " // TRIM( real_to_str( t_read ) )  // &
     260                    " s  (" // TRIM( real_to_str( 100 * t_read / t_total, fmt ) ) // " %)" )
     261
     262        CALL report('log_runtime', "Writing data:    " // TRIM( real_to_str( t_write ) ) // &
     263                    " s  (" // TRIM( real_to_str( 100 * t_write / t_total, fmt ) ) // " %)" )
     264
     265        CALL report('log_runtime', "Computation:     " // TRIM( real_to_str( t_comp ) )  // &
     266                    " s  (" // TRIM( real_to_str( 100 * t_comp / t_total, fmt) ) // " %)" )
     267
     268        CALL report('log_runtime', "Total:           " // TRIM( real_to_str( t_total ) ) // &
     269                    " s  (" // TRIM( real_to_str( 100 * t_total / t_total, fmt ) ) // " %)")
     270
     271    CASE DEFAULT
     272       CALL inifor_abort('log_runtime', "Mode '" // TRIM(mode) // "' is not supported.")
     273
     274    END SELECT
     275
     276 END SUBROUTINE log_runtime
    269277
    270278 END MODULE inifor_control
    271 #endif
    272 

  • palm/trunk/UTIL/inifor/src/inifor_defs.f90

    r3801 r3866  
    2626! -----------------
    2727! $Id$
     28! Added parameter for INIFOR's log file name
     29! Use PALM's working precision
     30!
     31!
     32! 3801 2019-03-15 17:14:25Z eckhard
    2833! Defined netCDF variable names for COSMO grid
    2934! Bumped version number
    … …  
    9499!> The defs module provides global constants used in INIFOR.
    95100!------------------------------------------------------------------------------!
    96 #if defined ( __netcdf )
    97101 MODULE inifor_defs
    98102 
     103 USE kinds,                                                                    &
     104     ONLY :  wp, iwp
     105
    99106 IMPLICIT NONE
    100107
    101108!
    102109!-- Parameters for type definitions
    103  INTEGER, PARAMETER  ::  dp    = 8   !< double precision (8 bytes = 64 bits)
    104  INTEGER, PARAMETER  ::  sp    = 4   !< single precision (4 bytes = 32 bits)
    105  INTEGER, PARAMETER  ::  hp    = 2   !< half precision (2 bytes = 16 bits)
    106110 INTEGER, PARAMETER  ::  PATH  = 140 !< length of file path strings
    107111 INTEGER, PARAMETER  ::  LNAME = 150 !< length of long name strings
    … …  
    111115!
    112116!-- Trigonomentry
    113  REAL(dp), PARAMETER ::  PI = 3.14159265358979323846264338_dp !< Ratio of a circle's circumference to its diamter [-]
    114  REAL(dp), PARAMETER ::  TO_RADIANS = PI / 180.0_dp           !< Conversion factor from degrees to radiant [-]
    115  REAL(dp), PARAMETER ::  TO_DEGREES = 180.0_dp / PI           !< Conversion factor from radians to degrees [-]
     117 REAL(wp), PARAMETER ::  PI = 3.14159265358979323846264338_wp !< Ratio of a circle's circumference to its diamter [-]
     118 REAL(wp), PARAMETER ::  TO_RADIANS = PI / 180.0_wp           !< Conversion factor from degrees to radiant [-]
     119 REAL(wp), PARAMETER ::  TO_DEGREES = 180.0_wp / PI           !< Conversion factor from radians to degrees [-]
    116120
    117121!
    118122!-- COSMO parameters
    119123 INTEGER, PARAMETER  ::  WATER_ID = 9                !< Integer corresponding to the water soil type in COSMO-DE [-]
    120  REAL(dp), PARAMETER ::  EARTH_RADIUS = 6371229.0_dp !< Earth radius used in COSMO-DE [m]
    121  REAL(dp), PARAMETER ::  P_SL = 1e5_dp               !< Reference pressure for computation of COSMO-DE's basic state pressure [Pa]
    122  REAL(dp), PARAMETER ::  T_SL = 288.15_dp            !< Reference temperature for computation of COSMO-DE's basic state pressure [K]
    123  REAL(dp), PARAMETER ::  BETA = 42.0_dp              !< logarithmic lapse rate, dT / d ln(p), for computation of COSMO-DE's basic
     124 REAL(wp), PARAMETER ::  EARTH_RADIUS = 6371229.0_wp !< Earth radius used in COSMO-DE [m]
     125 REAL(wp), PARAMETER ::  P_SL = 1e5_wp               !< Reference pressure for computation of COSMO-DE's basic state pressure [Pa]
     126 REAL(wp), PARAMETER ::  T_SL = 288.15_wp            !< Reference temperature for computation of COSMO-DE's basic state pressure [K]
     127 REAL(wp), PARAMETER ::  BETA = 42.0_wp              !< logarithmic lapse rate, dT / d ln(p), for computation of COSMO-DE's basic
    124128                                                     !< state pressure [K]
    125  REAL(dp), PARAMETER ::  RD   = 287.05_dp            !< specific gas constant of dry air, used in computation of COSMO-DE's basic
     129 REAL(wp), PARAMETER ::  RD   = 287.05_wp            !< specific gas constant of dry air, used in computation of COSMO-DE's basic
    126130                                                     !< state [J/kg/K]
    127  REAL(dp), PARAMETER ::  RV   = 461.51_dp            !< specific gas constant of water vapor [J/kg/K]
    128  REAL(dp), PARAMETER ::  G    = 9.80665_dp           !< acceleration of Earth's gravity, used in computation of COSMO-DE's basic
     131 REAL(wp), PARAMETER ::  RV   = 461.51_wp            !< specific gas constant of water vapor [J/kg/K]
     132 REAL(wp), PARAMETER ::  G    = 9.80665_wp           !< acceleration of Earth's gravity, used in computation of COSMO-DE's basic
    129133                                                     !< state [m/s/s]
    130  REAL(dp), PARAMETER ::  RHO_L = 1e3_dp              !< density of liquid water, used to convert W_SO from [kg/m^2] to [m^3/m^3],
     134 REAL(wp), PARAMETER ::  RHO_L = 1e3_wp              !< density of liquid water, used to convert W_SO from [kg/m^2] to [m^3/m^3],
    131135                                                     !< in [kg/m^3]
    132  REAL(dp), PARAMETER ::  HECTO = 100_dp              !< unit conversion factor from hPa to Pa
     136 REAL(wp), PARAMETER ::  HECTO = 100_wp              !< unit conversion factor from hPa to Pa
    133137
    134138!
    135139!-- PALM-4U parameters
    136  REAL(dp), PARAMETER ::  OMEGA   = 7.29e-5_dp !< angular velocity of Earth's rotation [s^-1]
    137  REAL(dp), PARAMETER ::  P_REF   = 1e5_dp     !< Reference pressure for potential temperature [Pa]
    138  REAL(dp), PARAMETER ::  RD_PALM = 287.0_dp   !< specific gas constant of dry air, used in computation of PALM-4U's potential temperature [J/kg/K]
    139  REAL(dp), PARAMETER ::  CP_PALM = 1005.0_dp  !< heat capacity of dry air at constant pressure, used in computation of PALM-4U's potential temperature [J/kg/K]
     140 REAL(wp), PARAMETER ::  OMEGA   = 7.29e-5_wp !< angular velocity of Earth's rotation [s^-1]
     141 REAL(wp), PARAMETER ::  P_REF   = 1e5_wp     !< Reference pressure for potential temperature [Pa]
     142 REAL(wp), PARAMETER ::  RD_PALM = 287.0_wp   !< specific gas constant of dry air, used in computation of PALM-4U's potential temperature [J/kg/K]
     143 REAL(wp), PARAMETER ::  CP_PALM = 1005.0_wp  !< heat capacity of dry air at constant pressure, used in computation of PALM-4U's potential temperature [J/kg/K]
    140144
    141145!
    … …  
    154158                                                              !< water cells [-]
    155159 INTEGER, PARAMETER          ::  FORCING_STEP = 1             !< Number of hours between forcing time steps [h]
    156  REAL(dp), PARAMETER         ::  NUDGING_TAU = 21600.0_dp     !< Nudging relaxation time scale [s]
    157  CHARACTER(LEN=*), PARAMETER ::  VERSION = '1.4.8'            !< INIFOR version number
     160 REAL(wp), PARAMETER         ::  NUDGING_TAU = 21600.0_wp     !< Nudging relaxation time scale [s]
    158161 CHARACTER(LEN=*), PARAMETER ::  COPYRIGHT = 'Copyright 2017-2019 Leibniz Universitaet Hannover' // &
    159      ACHAR( 10 ) // ' Copyright 2017-2019 Deutscher Wetterdienst Offenbach' !< Copyright notice
    160 
     162    ACHAR( 10 ) // ' Copyright 2017-2019 Deutscher Wetterdienst Offenbach' !< Copyright notice
     163 CHARACTER(LEN=*), PARAMETER ::  LOG_FILE_NAME = 'inifor.log' !< Name of INIFOR's log file
     164 CHARACTER(LEN=*), PARAMETER ::  VERSION = '1.4.9rc'          !< INIFOR version number
     165 
    161166 END MODULE inifor_defs
    162 #endif

  • palm/trunk/UTIL/inifor/src/inifor_grid.f90

    r3802 r3866  
    2626! -----------------
    2727! $Id$
     28! Use PALM's working precision
     29! Catch errors while reading namelists
     30! Improved coding style
     31!
     32!
     33! 3802 2019-03-17 13:33:42Z raasch
    2834! unused variable removed
    2935!
    … …  
    131137    USE inifor_control
    132138    USE inifor_defs,                                                           &
    133         ONLY:  DATE, EARTH_RADIUS, TO_RADIANS, TO_DEGREES, PI, dp, hp, sp,     &
     139        ONLY:  DATE, EARTH_RADIUS, TO_RADIANS, TO_DEGREES, PI,                 &
    134140               SNAME, LNAME, PATH, FORCING_STEP, WATER_ID, FILL_ITERATIONS,    &
    135141               BETA, P_SL, T_SL, BETA, RD, RV, G, P_REF, RD_PALM, CP_PALM,     &
    136                RHO_L, OMEGA, HECTO
     142               RHO_L, OMEGA, HECTO, wp, iwp
    137143    USE inifor_io,                                                             &
    138144        ONLY:  get_cosmo_grid, get_netcdf_attribute, get_netcdf_dim_vector,    &
    … …  
    140146               get_input_file_list, validate_config
    141147    USE inifor_transform,                                                      &
    142         ONLY:  average_2d, rotate_to_cosmo, find_horizontal_neighbours,&
     148        ONLY:  average_2d, rotate_to_cosmo, find_horizontal_neighbours,        &
    143149               compute_horizontal_interp_weights,                              &
    144150               find_vertical_neighbours_and_weights_interp,                    &
    … …  
    156162    SAVE
    157163   
    158     REAL(dp) ::  averaging_angle   = 0.0_dp       !< latitudal and longitudal width of averaging regions [rad]
    159     REAL(dp) ::  averaging_width_ns = 0.0_dp       !< longitudal width of averaging regions [m]
    160     REAL(dp) ::  averaging_width_ew = 0.0_dp       !< latitudal width of averaging regions [m]
    161     REAL(dp) ::  phi_equat         = 0.0_dp       !< latitude of rotated equator of COSMO-DE grid [rad]
    162     REAL(dp) ::  phi_n             = 0.0_dp       !< latitude of rotated pole of COSMO-DE grid [rad]
    163     REAL(dp) ::  lambda_n          = 0.0_dp       !< longitude of rotaded pole of COSMO-DE grid [rad]
    164     REAL(dp) ::  phi_c             = 0.0_dp       !< rotated-grid latitude of the center of the PALM domain [rad]
    165     REAL(dp) ::  lambda_c          = 0.0_dp       !< rotated-grid longitude of the centre of the PALM domain [rad]
    166     REAL(dp) ::  phi_cn            = 0.0_dp       !< latitude of the rotated pole relative to the COSMO-DE grid [rad]
    167     REAL(dp) ::  lambda_cn         = 0.0_dp       !< longitude of the rotated pole relative to the COSMO-DE grid [rad]
    168     REAL(dp) ::  lam_centre        = 0.0_dp       !< longitude of the PLAM domain centre in the source (COSMO rotated-pole) system [rad]
    169     REAL(dp) ::  phi_centre        = 0.0_dp       !< latitude of the PLAM domain centre in the source (COSMO rotated-pole) system [rad]
    170     REAL(dp) ::  lam_east          = 0.0_dp       !< longitude of the east central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad]
    171     REAL(dp) ::  lam_west          = 0.0_dp       !< longitude of the west central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad]
    172     REAL(dp) ::  phi_north         = 0.0_dp       !< latitude of the north central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad]
    173     REAL(dp) ::  phi_south         = 0.0_dp       !< latitude of the south central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad]
    174     REAL(dp) ::  gam               = 0.0_dp       !< angle for working around phirot2phi/rlarot2rla bug
    175     REAL(dp) ::  dx                = 0.0_dp       !< PALM-4U grid spacing in x direction [m]
    176     REAL(dp) ::  dy                = 0.0_dp       !< PALM-4U grid spacing in y direction [m]
    177     REAL(dp) ::  dz(10)            = -1.0_dp      !< PALM-4U grid spacing in z direction [m]
    178     REAL(dp) ::  dz_max            = 1000.0_dp    !< maximum vertical grid spacing [m]
    179     REAL(dp) ::  dz_stretch_factor = 1.08_dp      !< factor for vertical grid stretching [m]
    180     REAL(dp) ::  dz_stretch_level  = -9999999.9_dp!< height above which the vertical grid will be stretched [m]
    181     REAL(dp) ::  dz_stretch_level_start(9) = -9999999.9_dp !< namelist parameter
    182     REAL(dp) ::  dz_stretch_level_end(9) = 9999999.9_dp !< namelist parameter
    183     REAL(dp) ::  dz_stretch_factor_array(9) = 1.08_dp !< namelist parameter
    184     REAL(dp) ::  dxi               = 0.0_dp       !< inverse PALM-4U grid spacing in x direction [m^-1]
    185     REAL(dp) ::  dyi               = 0.0_dp       !< inverse PALM-4U grid spacing in y direction [m^-1]
    186     REAL(dp) ::  dzi               = 0.0_dp       !< inverse PALM-4U grid spacing in z direction [m^-1]
    187     REAL(dp) ::  f3                = 0.0_dp       !< Coriolis parameter
    188     REAL(dp) ::  lx                = 0.0_dp       !< PALM-4U domain size in x direction [m]
    189     REAL(dp) ::  ly                = 0.0_dp       !< PALM-4U domain size in y direction [m]
    190     REAL(dp) ::  p0                = 0.0_dp       !< PALM-4U surface pressure, at z0 [Pa]
    191     REAL(dp) ::  x0                = 0.0_dp       !< x coordinate of PALM-4U Earth tangent [m]
    192     REAL(dp) ::  y0                = 0.0_dp       !< y coordinate of PALM-4U Earth tangent [m]
    193     REAL(dp) ::  z0                = 0.0_dp       !< Elevation of the PALM-4U domain above sea level [m]
    194     REAL(dp) ::  z_top             = 0.0_dp       !< height of the scalar top boundary [m]
    195     REAL(dp) ::  zw_top            = 0.0_dp       !< height of the vertical velocity top boundary [m]
    196     REAL(dp) ::  lonmin_cosmo      = 0.0_dp       !< Minimunm longitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
    197     REAL(dp) ::  lonmax_cosmo      = 0.0_dp       !< Maximum longitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
    198     REAL(dp) ::  latmin_cosmo      = 0.0_dp       !< Minimunm latitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
    199     REAL(dp) ::  latmax_cosmo      = 0.0_dp       !< Maximum latitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
    200     REAL(dp) ::  lonmin_palm       = 0.0_dp       !< Minimunm longitude of PALM grid [COSMO rotated-pole rad]
    201     REAL(dp) ::  lonmax_palm       = 0.0_dp       !< Maximum longitude of PALM grid [COSMO rotated-pole rad]
    202     REAL(dp) ::  latmin_palm       = 0.0_dp       !< Minimunm latitude of PALM grid [COSMO rotated-pole rad]
    203     REAL(dp) ::  latmax_palm       = 0.0_dp       !< Maximum latitude of PALM grid [COSMO rotated-pole rad]
    204     REAL(dp) ::  lonmin_tot        = 0.0_dp       !< Minimunm longitude of required COSMO data [COSMO rotated-pole rad]
    205     REAL(dp) ::  lonmax_tot        = 0.0_dp       !< Maximum longitude of required COSMO data [COSMO rotated-pole rad]
    206     REAL(dp) ::  latmin_tot        = 0.0_dp       !< Minimunm latitude of required COSMO data [COSMO rotated-pole rad]
    207     REAL(dp) ::  latmax_tot        = 0.0_dp       !< Maximum latitude of required COSMO data [COSMO rotated-pole rad]
    208     REAL(dp) ::  latitude          = 0.0_dp       !< geographical latitude of the PALM-4U origin, from inipar namelist [deg]
    209     REAL(dp) ::  longitude         = 0.0_dp       !< geographical longitude of the PALM-4U origin, from inipar namelist [deg]
    210     REAL(dp) ::  origin_lat        = 0.0_dp       !< geographical latitude of the PALM-4U origin, from static driver netCDF file [deg]
    211     REAL(dp) ::  origin_lon        = 0.0_dp       !< geographical longitude of the PALM-4U origin, from static driver netCDF file [deg]
    212     REAL(dp) ::  rotation_angle    = 0.0_dp       !< clockwise angle the PALM-4U north is rotated away from geographical north [deg]
    213     REAL(dp) ::  end_time          = 0.0_dp       !< PALM-4U simulation time [s]
    214 
    215     REAL(dp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  hhl             !< heights of half layers (cell faces) above sea level in COSMO-DE, read in from external file
    216     REAL(dp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  hfl             !< heights of full layers (cell centres) above sea level in COSMO-DE, computed as arithmetic average of hhl
    217     REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET     ::  depths          !< COSMO-DE's TERRA-ML soil layer depths
    218     REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET     ::  d_depth         !< COSMO-DE's TERRA-ML soil layer thicknesses
    219     REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET     ::  d_depth_rho_inv !< inverted soil water mass
    220     REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET     ::  rlon            !< longitudes of COSMO-DE's rotated-pole grid
    221     REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET     ::  rlat            !< latitudes of COSMO-DE's rotated-pole grid
    222     REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET     ::  time            !< output times
    223     REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET     ::  x               !< base palm grid x coordinate vector pointed to by grid_definitions
    224     REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET     ::  xu              !< base palm grid xu coordinate vector pointed to by grid_definitions
    225     REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET     ::  y               !< base palm grid y coordinate vector pointed to by grid_definitions
    226     REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET     ::  yv              !< base palm grid yv coordinate vector pointed to by grid_definitions
    227     REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET     ::  z_column        !< base palm grid z coordinate vector including the top boundary coordinate (entire column)
    228     REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET     ::  zw_column       !< base palm grid zw coordinate vector including the top boundary coordinate (entire column)
    229     REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET     ::  z               !< base palm grid z coordinate vector pointed to by grid_definitions
    230     REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET     ::  zw              !< base palm grid zw coordinate vector pointed to by grid_definitions
    231 
    232     INTEGER(hp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  soiltyp      !< COSMO-DE soil type map
     164    REAL(wp) ::  averaging_angle   = 0.0_wp       !< latitudal and longitudal width of averaging regions [rad]
     165    REAL(wp) ::  averaging_width_ns = 0.0_wp       !< longitudal width of averaging regions [m]
     166    REAL(wp) ::  averaging_width_ew = 0.0_wp       !< latitudal width of averaging regions [m]
     167    REAL(wp) ::  phi_equat         = 0.0_wp       !< latitude of rotated equator of COSMO-DE grid [rad]
     168    REAL(wp) ::  phi_n             = 0.0_wp       !< latitude of rotated pole of COSMO-DE grid [rad]
     169    REAL(wp) ::  lambda_n          = 0.0_wp       !< longitude of rotaded pole of COSMO-DE grid [rad]
     170    REAL(wp) ::  phi_c             = 0.0_wp       !< rotated-grid latitude of the center of the PALM domain [rad]
     171    REAL(wp) ::  lambda_c          = 0.0_wp       !< rotated-grid longitude of the centre of the PALM domain [rad]
     172    REAL(wp) ::  phi_cn            = 0.0_wp       !< latitude of the rotated pole relative to the COSMO-DE grid [rad]
     173    REAL(wp) ::  lambda_cn         = 0.0_wp       !< longitude of the rotated pole relative to the COSMO-DE grid [rad]
     174    REAL(wp) ::  lam_centre        = 0.0_wp       !< longitude of the PLAM domain centre in the source (COSMO rotated-pole) system [rad]
     175    REAL(wp) ::  phi_centre        = 0.0_wp       !< latitude of the PLAM domain centre in the source (COSMO rotated-pole) system [rad]
     176    REAL(wp) ::  lam_east          = 0.0_wp       !< longitude of the east central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad]
     177    REAL(wp) ::  lam_west          = 0.0_wp       !< longitude of the west central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad]
     178    REAL(wp) ::  phi_north         = 0.0_wp       !< latitude of the north central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad]
     179    REAL(wp) ::  phi_south         = 0.0_wp       !< latitude of the south central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad]
     180    REAL(wp) ::  gam               = 0.0_wp       !< angle for working around phirot2phi/rlarot2rla bug
     181    REAL(wp) ::  dx                = 0.0_wp       !< PALM-4U grid spacing in x direction [m]
     182    REAL(wp) ::  dy                = 0.0_wp       !< PALM-4U grid spacing in y direction [m]
     183    REAL(wp) ::  dz(10)            = -1.0_wp      !< PALM-4U grid spacing in z direction [m]
     184    REAL(wp) ::  dz_max            = 1000.0_wp    !< maximum vertical grid spacing [m]
     185    REAL(wp) ::  dz_stretch_factor = 1.08_wp      !< factor for vertical grid stretching [m]
     186    REAL(wp) ::  dz_stretch_level  = -9999999.9_wp!< height above which the vertical grid will be stretched [m]
     187    REAL(wp) ::  dz_stretch_level_start(9) = -9999999.9_wp !< namelist parameter
     188    REAL(wp) ::  dz_stretch_level_end(9) = 9999999.9_wp !< namelist parameter
     189    REAL(wp) ::  dz_stretch_factor_array(9) = 1.08_wp !< namelist parameter
     190    REAL(wp) ::  dxi               = 0.0_wp       !< inverse PALM-4U grid spacing in x direction [m^-1]
     191    REAL(wp) ::  dyi               = 0.0_wp       !< inverse PALM-4U grid spacing in y direction [m^-1]
     192    REAL(wp) ::  dzi               = 0.0_wp       !< inverse PALM-4U grid spacing in z direction [m^-1]
     193    REAL(wp) ::  f3                = 0.0_wp       !< Coriolis parameter
     194    REAL(wp) ::  lx                = 0.0_wp       !< PALM-4U domain size in x direction [m]
     195    REAL(wp) ::  ly                = 0.0_wp       !< PALM-4U domain size in y direction [m]
     196    REAL(wp) ::  p0                = 0.0_wp       !< PALM-4U surface pressure, at z0 [Pa]
     197    REAL(wp) ::  x0                = 0.0_wp       !< x coordinate of PALM-4U Earth tangent [m]
     198    REAL(wp) ::  y0                = 0.0_wp       !< y coordinate of PALM-4U Earth tangent [m]
     199    REAL(wp) ::  z0                = 0.0_wp       !< Elevation of the PALM-4U domain above sea level [m]
     200    REAL(wp) ::  z_top             = 0.0_wp       !< height of the scalar top boundary [m]
     201    REAL(wp) ::  zw_top            = 0.0_wp       !< height of the vertical velocity top boundary [m]
     202    REAL(wp) ::  lonmin_cosmo      = 0.0_wp       !< Minimunm longitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
     203    REAL(wp) ::  lonmax_cosmo      = 0.0_wp       !< Maximum longitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
     204    REAL(wp) ::  latmin_cosmo      = 0.0_wp       !< Minimunm latitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
     205    REAL(wp) ::  latmax_cosmo      = 0.0_wp       !< Maximum latitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
     206    REAL(wp) ::  lonmin_palm       = 0.0_wp       !< Minimunm longitude of PALM grid [COSMO rotated-pole rad]
     207    REAL(wp) ::  lonmax_palm       = 0.0_wp       !< Maximum longitude of PALM grid [COSMO rotated-pole rad]
     208    REAL(wp) ::  latmin_palm       = 0.0_wp       !< Minimunm latitude of PALM grid [COSMO rotated-pole rad]
     209    REAL(wp) ::  latmax_palm       = 0.0_wp       !< Maximum latitude of PALM grid [COSMO rotated-pole rad]
     210    REAL(wp) ::  lonmin_tot        = 0.0_wp       !< Minimunm longitude of required COSMO data [COSMO rotated-pole rad]
     211    REAL(wp) ::  lonmax_tot        = 0.0_wp       !< Maximum longitude of required COSMO data [COSMO rotated-pole rad]
     212    REAL(wp) ::  latmin_tot        = 0.0_wp       !< Minimunm latitude of required COSMO data [COSMO rotated-pole rad]
     213    REAL(wp) ::  latmax_tot        = 0.0_wp       !< Maximum latitude of required COSMO data [COSMO rotated-pole rad]
     214    REAL(wp) ::  latitude          = 0.0_wp       !< geographical latitude of the PALM-4U origin, from inipar namelist [deg]
     215    REAL(wp) ::  longitude         = 0.0_wp       !< geographical longitude of the PALM-4U origin, from inipar namelist [deg]
     216    REAL(wp) ::  origin_lat        = 0.0_wp       !< geographical latitude of the PALM-4U origin, from static driver netCDF file [deg]
     217    REAL(wp) ::  origin_lon        = 0.0_wp       !< geographical longitude of the PALM-4U origin, from static driver netCDF file [deg]
     218    REAL(wp) ::  rotation_angle    = 0.0_wp       !< clockwise angle the PALM-4U north is rotated away from geographical north [deg]
     219    REAL(wp) ::  end_time          = 0.0_wp       !< PALM-4U simulation time [s]
     220
     221    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  hhl             !< heights of half layers (cell faces) above sea level in COSMO-DE, read in from external file
     222    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  hfl             !< heights of full layers (cell centres) above sea level in COSMO-DE, computed as arithmetic average of hhl
     223    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET     ::  depths          !< COSMO-DE's TERRA-ML soil layer depths
     224    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET     ::  d_depth         !< COSMO-DE's TERRA-ML soil layer thicknesses
     225    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET     ::  d_depth_rho_inv !< inverted soil water mass
     226    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET     ::  rlon            !< longitudes of COSMO-DE's rotated-pole grid
     227    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET     ::  rlat            !< latitudes of COSMO-DE's rotated-pole grid
     228    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET     ::  time            !< output times
     229    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET     ::  x               !< base palm grid x coordinate vector pointed to by grid_definitions
     230    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET     ::  xu              !< base palm grid xu coordinate vector pointed to by grid_definitions
     231    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET     ::  y               !< base palm grid y coordinate vector pointed to by grid_definitions
     232    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET     ::  yv              !< base palm grid yv coordinate vector pointed to by grid_definitions
     233    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET     ::  z_column        !< base palm grid z coordinate vector including the top boundary coordinate (entire column)
     234    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET     ::  zw_column       !< base palm grid zw coordinate vector including the top boundary coordinate (entire column)
     235    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET     ::  z               !< base palm grid z coordinate vector pointed to by grid_definitions
     236    REAL(wp), DIMENSION(:), ALLOCATABLE, TARGET     ::  zw              !< base palm grid zw coordinate vector pointed to by grid_definitions
     237
     238    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  soiltyp     !< COSMO-DE soil type map
    233239    INTEGER ::  dz_stretch_level_end_index(9)               !< vertical grid level index until which the vertical grid spacing is stretched
    234240    INTEGER ::  dz_stretch_level_start_index(9)             !< vertical grid level index above which the vertical grid spacing is stretched
     241    INTEGER ::  iostat !< return status of READ statement
    235242    INTEGER ::  nt    !< number of output time steps
    236243    INTEGER ::  nx    !< number of PALM-4U grid points in x direction
    … …  
    355362!> interplation grids later in setup_grids().
    356363!------------------------------------------------------------------------------!
    357     SUBROUTINE setup_parameters()
     364 SUBROUTINE setup_parameters()
    358365
    359366!
    … …  
    361368! Section 1: Define default parameters
    362369!------------------------------------------------------------------------------
    363        cfg % start_date = '2013072100'
    364        end_hour = 2
    365        start_hour_soil = -2
    366        start_hour_soilmoisture = - (4 * 7 * 24) - 2
    367 
    368 !
    369 !--    Defaultmain centre (_c) of the PALM-4U grid in the geographical system (_g)
    370        origin_lat = 52.325079_dp * TO_RADIANS ! south-west of Berlin, origin used for the Dec 2017 showcase simulation
    371        origin_lon = 13.082744_dp * TO_RADIANS
    372        cfg % z0 = 35.0_dp
    373 
    374 !
    375 !--    Default atmospheric parameters
    376        cfg % ug = 0.0_dp
    377        cfg % vg = 0.0_dp
    378        cfg % p0 = P_SL
    379 
    380 !
    381 !--    Parameters for file names
    382        start_hour_flow = 0
    383        start_hour_soil = 0
    384        start_hour_radiation = 0
    385        start_hour_soilmoisture = start_hour_flow - 2
    386        end_hour_soilmoisture = start_hour_flow
    387        step_hour = FORCING_STEP
    388 
    389        input_prefix = 'laf'
    390        cfg % flow_prefix = input_prefix
    391        cfg % input_prefix = input_prefix
    392        cfg % soil_prefix = input_prefix
    393        cfg % radiation_prefix = input_prefix
    394        cfg % soilmoisture_prefix  = input_prefix
    395 
    396        flow_suffix = '-flow'
    397        soil_suffix = '-soil'
    398        radiation_suffix = '-rad'
    399        soilmoisture_suffix = '-soilmoisture'
    400 
    401        cfg % debug = .FALSE.
    402        cfg % averaging_angle = 2.0_dp
     370    cfg%start_date = '2013072100'
     371    end_hour = 2
     372    start_hour_soil = -2
     373    start_hour_soilmoisture = - (4 * 7 * 24) - 2
     374
     375!
     376!-- Defaultmain centre (_c) of the PALM-4U grid in the geographical system (_g)
     377    origin_lat = 52.325079_wp * TO_RADIANS ! south-west of Berlin, origin used for the Dec 2017 showcase simulation
     378    origin_lon = 13.082744_wp * TO_RADIANS
     379    cfg%z0 = 35.0_wp
     380
     381!
     382!-- Default atmospheric parameters
     383    cfg%ug = 0.0_wp
     384    cfg%vg = 0.0_wp
     385    cfg%p0 = P_SL
     386
     387!
     388!-- Parameters for file names
     389    start_hour_flow = 0
     390    start_hour_soil = 0
     391    start_hour_radiation = 0
     392    start_hour_soilmoisture = start_hour_flow - 2
     393    end_hour_soilmoisture = start_hour_flow
     394    step_hour = FORCING_STEP
     395
     396    input_prefix = 'laf'
     397    cfg%flow_prefix = input_prefix
     398    cfg%input_prefix = input_prefix
     399    cfg%soil_prefix = input_prefix
     400    cfg%radiation_prefix = input_prefix
     401    cfg%soilmoisture_prefix  = input_prefix
     402
     403    flow_suffix = '-flow'
     404    soil_suffix = '-soil'
     405    radiation_suffix = '-rad'
     406    soilmoisture_suffix = '-soilmoisture'
     407
     408    cfg%debug = .FALSE.
     409    cfg%averaging_angle = 2.0_wp
    403410!
    404411!------------------------------------------------------------------------------
    … …  
    407414
    408415!
    409 !--    Set default paths and modes
    410        cfg % input_path         = './'
    411        cfg % hhl_file           = ''
    412        cfg % soiltyp_file       = ''
    413        cfg % namelist_file      = './namelist'
    414        cfg % static_driver_file = ''
    415        cfg % output_file = './palm-4u-input.nc'
    416        cfg % ic_mode = 'volume'
    417        cfg % bc_mode = 'real'
    418        cfg % averaging_mode = 'level'
    419 
    420 !
    421 !--    Overwrite defaults with user configuration
    422        CALL parse_command_line_arguments( cfg )
    423        CALL report('main_loop', 'Running INIFOR version ' // VERSION)
    424 
    425        flow_prefix = TRIM(cfg % input_prefix)
    426        radiation_prefix = TRIM(cfg % input_prefix)
    427        soil_prefix = TRIM(cfg % input_prefix)
    428        soilmoisture_prefix = TRIM(cfg % input_prefix)
    429        IF (cfg % flow_prefix_is_set)  flow_prefix = TRIM(cfg % flow_prefix)
    430        IF (cfg % radiation_prefix_is_set)  radiation_prefix = TRIM(cfg % radiation_prefix)
    431        IF (cfg % soil_prefix_is_set)  soil_prefix = TRIM(cfg % soil_prefix)
    432        IF (cfg % soilmoisture_prefix_is_set)  soilmoisture_prefix = TRIM(cfg % soilmoisture_prefix)
    433 
    434        output_file % name = cfg % output_file
    435        z0 = cfg % z0
    436        p0 = cfg % p0
    437 
    438        init_variables_required = .TRUE.
    439        boundary_variables_required = TRIM( cfg % bc_mode ) == 'real'
    440        ls_forcing_variables_required = TRIM( cfg % bc_mode ) == 'ideal'
    441        surface_forcing_required = .FALSE.
    442 
    443        IF ( ls_forcing_variables_required )  THEN
    444           message = "Averaging of large-scale forcing profiles " //            &
    445                     "has not been implemented, yet."
    446           CALL inifor_abort('setup_parameters', message)
    447        ENDIF
    448 
    449 !
    450 !--    Set default file paths, if not specified by user.
    451        CALL normalize_path(cfg % input_path)
    452        IF (TRIM(cfg % hhl_file) == '')  cfg % hhl_file = TRIM(cfg % input_path) // 'hhl.nc'
    453        IF (TRIM(cfg % soiltyp_file) == '')  cfg % soiltyp_file = TRIM(cfg % input_path) // 'soil.nc'
    454 
    455        CALL validate_config( cfg )
    456 
    457        CALL report('setup_parameters', "initialization mode: " // TRIM(cfg % ic_mode))
    458        CALL report('setup_parameters', "       forcing mode: " // TRIM(cfg % bc_mode))
    459        CALL report('setup_parameters', "     averaging mode: " // TRIM(cfg % averaging_mode))
    460        CALL report('setup_parameters', "    averaging angle: " // real_to_str(cfg % averaging_angle))
    461        CALL report('setup_parameters', "          data path: " // TRIM(cfg % input_path))
    462        CALL report('setup_parameters', "           hhl file: " // TRIM(cfg % hhl_file))
    463        CALL report('setup_parameters', "       soiltyp file: " // TRIM(cfg % soiltyp_file))
    464        CALL report('setup_parameters', "      namelist file: " // TRIM(cfg % namelist_file))
    465        CALL report('setup_parameters', "   output data file: " // TRIM(output_file % name))
    466        IF (cfg % debug )  CALL report('setup_parameters', "     debugging mode: enabled")
    467 
    468  CALL run_control('time', 'init')
    469  !
    470  !--   Read in namelist parameters
    471        OPEN(10, FILE=cfg % namelist_file)
    472        READ(10, NML=inipar) ! nx, ny, nz, dx, dy, dz
    473        READ(10, NML=d3par)  ! end_time
     416!-- Set default paths and modes
     417    cfg%input_path         = './'
     418    cfg%hhl_file           = ''
     419    cfg%soiltyp_file       = ''
     420    cfg%namelist_file      = './namelist'
     421    cfg%static_driver_file = ''
     422    cfg%output_file = './palm-4u-input.nc'
     423    cfg%ic_mode = 'volume'
     424    cfg%bc_mode = 'real'
     425    cfg%averaging_mode = 'level'
     426
     427!
     428!-- Overwrite defaults with user configuration
     429    CALL parse_command_line_arguments( cfg )
     430    CALL report('main_loop', 'Running INIFOR version ' // VERSION)
     431
     432    flow_prefix = TRIM(cfg%input_prefix)
     433    radiation_prefix = TRIM(cfg%input_prefix)
     434    soil_prefix = TRIM(cfg%input_prefix)
     435    soilmoisture_prefix = TRIM(cfg%input_prefix)
     436    IF (cfg%flow_prefix_is_set)  flow_prefix = TRIM(cfg%flow_prefix)
     437    IF (cfg%radiation_prefix_is_set)  radiation_prefix = TRIM(cfg%radiation_prefix)
     438    IF (cfg%soil_prefix_is_set)  soil_prefix = TRIM(cfg%soil_prefix)
     439    IF (cfg%soilmoisture_prefix_is_set)  soilmoisture_prefix = TRIM(cfg%soilmoisture_prefix)
     440
     441    output_file%name = cfg%output_file
     442    z0 = cfg%z0
     443    p0 = cfg%p0
     444
     445    init_variables_required = .TRUE.
     446    boundary_variables_required = TRIM( cfg%bc_mode ) == 'real'
     447    ls_forcing_variables_required = TRIM( cfg%bc_mode ) == 'ideal'
     448    surface_forcing_required = .FALSE.
     449
     450    IF ( ls_forcing_variables_required )  THEN
     451       message = "Averaging of large-scale forcing profiles " //            &
     452                 "has not been implemented, yet."
     453       CALL inifor_abort('setup_parameters', message)
     454    ENDIF
     455
     456!
     457!-- Set default file paths, if not specified by user.
     458    CALL normalize_path(cfg%input_path)
     459    IF (TRIM(cfg%hhl_file) == '')  cfg%hhl_file = TRIM(cfg%input_path) // 'hhl.nc'
     460    IF (TRIM(cfg%soiltyp_file) == '')  cfg%soiltyp_file = TRIM(cfg%input_path) // 'soil.nc'
     461
     462    CALL validate_config( cfg )
     463
     464    CALL report('setup_parameters', "initialization mode: " // TRIM(cfg%ic_mode))
     465    CALL report('setup_parameters', "       forcing mode: " // TRIM(cfg%bc_mode))
     466    CALL report('setup_parameters', "     averaging mode: " // TRIM(cfg%averaging_mode))
     467    CALL report('setup_parameters', "    averaging angle: " // real_to_str(cfg%averaging_angle))
     468    CALL report('setup_parameters', "          data path: " // TRIM(cfg%input_path))
     469    CALL report('setup_parameters', "           hhl file: " // TRIM(cfg%hhl_file))
     470    CALL report('setup_parameters', "       soiltyp file: " // TRIM(cfg%soiltyp_file))
     471    CALL report('setup_parameters', "      namelist file: " // TRIM(cfg%namelist_file))
     472    CALL report('setup_parameters', "   output data file: " // TRIM(output_file%name))
     473    IF (cfg%debug )  CALL report('setup_parameters', "     debugging mode: enabled")
     474
     475    CALL log_runtime('time', 'init')
     476!
     477!-- Read in namelist parameters
     478    OPEN(10, FILE=cfg%namelist_file, STATUS='old')
     479    READ(10, NML=inipar, IOSTAT=iostat) ! nx, ny, nz, dx, dy, dz   
     480    IF ( iostat > 0 )  THEN     
     481       message = "Failed to read namelist 'inipar' from file '" //             &
     482                 TRIM( cfg%namelist_file ) // "'. "
     483       CALL inifor_abort( 'setup_parameters', message )
    474484       CLOSE(10)
    475  CALL run_control('time', 'read')
    476 
    477        end_hour = CEILING( end_time / 3600.0 * step_hour )
    478 
    479 !
    480 !--    Generate input file lists
    481        CALL get_input_file_list(                                               &
    482           cfg % start_date, start_hour_flow, end_hour, step_hour,              &
    483           cfg % input_path, flow_prefix, flow_suffix, flow_files)
    484        CALL get_input_file_list(                                               &
    485           cfg % start_date, start_hour_soil, end_hour, step_hour,              &
    486           cfg % input_path, soil_prefix, soil_suffix, soil_files)
    487        CALL get_input_file_list(                                               &
    488           cfg % start_date, start_hour_radiation, end_hour, step_hour,         &
    489           cfg % input_path, radiation_prefix, radiation_suffix, radiation_files, nocheck=.TRUE.)
    490        CALL get_input_file_list(                                               &
    491           cfg % start_date, start_hour_soilmoisture, end_hour_soilmoisture, step_hour, &
    492           cfg % input_path, soilmoisture_prefix, soilmoisture_suffix, soil_moisture_files, nocheck=.TRUE.)
     485    ENDIF
     486
     487    READ(10, NML=d3par, IOSTAT=iostat)  ! end_time
     488    IF ( iostat > 0 )  THEN
     489       message = "Failed to read namelist 'd3par' from file '" //              &
     490                 TRIM( cfg%namelist_file ) // "'. "
     491       CALL inifor_abort( 'setup_parameters', message )
     492       CLOSE(10)
     493    ENDIF
     494    CLOSE(10)
     495   
     496    CALL log_runtime('time', 'read')
     497
     498    end_hour = CEILING( end_time / 3600.0 * step_hour )
     499
     500!
     501!-- Generate input file lists
     502    CALL get_input_file_list(                                               &
     503       cfg%start_date, start_hour_flow, end_hour, step_hour,              &
     504       cfg%input_path, flow_prefix, flow_suffix, flow_files)
     505    CALL get_input_file_list(                                               &
     506       cfg%start_date, start_hour_soil, end_hour, step_hour,              &
     507       cfg%input_path, soil_prefix, soil_suffix, soil_files)
     508    CALL get_input_file_list(                                               &
     509       cfg%start_date, start_hour_radiation, end_hour, step_hour,         &
     510       cfg%input_path, radiation_prefix, radiation_suffix, radiation_files, nocheck=.TRUE.)
     511    CALL get_input_file_list(                                               &
     512       cfg%start_date, start_hour_soilmoisture, end_hour_soilmoisture, step_hour, &
     513       cfg%input_path, soilmoisture_prefix, soilmoisture_suffix, soil_moisture_files, nocheck=.TRUE.)
    493514
    494515!
    … …  
    506527
    507528
    508  CALL run_control('time', 'init')
    509 !
    510 !--    Read COSMO soil type map
    511        cosmo_var % name = 'SOILTYP'
    512        CALL get_netcdf_variable(cfg % soiltyp_file, cosmo_var, soiltyp)
    513 
    514        message = 'Reading PALM-4U origin from'
    515        IF (TRIM(cfg % static_driver_file) .NE. '')  THEN
    516 
    517           origin_lon = get_netcdf_attribute(cfg % static_driver_file, 'origin_lon')
    518           origin_lat = get_netcdf_attribute(cfg % static_driver_file, 'origin_lat')
    519 
    520           message = TRIM(message) // " static driver file '"                   &
    521                                   // TRIM(cfg % static_driver_file) // "'"
    522 
    523 
    524        ELSE
    525 
    526           origin_lon = longitude
    527           origin_lat = latitude
    528 
    529           message = TRIM(message) // " namlist file '"                         &
    530                                   // TRIM(cfg % namelist_file) // "'"
    531 
    532        ENDIF
    533        origin_lon = origin_lon * TO_RADIANS
    534        origin_lat = origin_lat * TO_RADIANS
    535 
    536        CALL report('setup_parameters', message)
    537 
    538  CALL run_control('time', 'read')
    539 
    540        CALL get_cosmo_grid( cfg, soil_files(1), rlon, rlat, hhl, hfl, depths,  &
    541                             d_depth, d_depth_rho_inv, phi_n, lambda_n,         &
    542                             phi_equat,                                         &
    543                             lonmin_cosmo, lonmax_cosmo,                        &
    544                             latmin_cosmo, latmax_cosmo,                        &
    545                             nlon, nlat, nlev, ndepths )
     529    CALL log_runtime('time', 'init')
     530!
     531!-- Read COSMO soil type map
     532    cosmo_var%name = 'SOILTYP'
     533    CALL get_netcdf_variable(cfg%soiltyp_file, cosmo_var, soiltyp)
     534
     535    message = 'Reading PALM-4U origin from'
     536    IF (TRIM(cfg%static_driver_file) .NE. '')  THEN
     537
     538       origin_lon = get_netcdf_attribute(cfg%static_driver_file, 'origin_lon')
     539       origin_lat = get_netcdf_attribute(cfg%static_driver_file, 'origin_lat')
     540
     541       message = TRIM(message) // " static driver file '"                   &
     542                               // TRIM(cfg%static_driver_file) // "'"
     543
     544
     545    ELSE
     546
     547       origin_lon = longitude
     548       origin_lat = latitude
     549
     550       message = TRIM(message) // " namlist file '"                         &
     551                               // TRIM(cfg%namelist_file) // "'"
     552
     553    ENDIF
     554    origin_lon = origin_lon * TO_RADIANS
     555    origin_lat = origin_lat * TO_RADIANS
     556
     557    CALL report('setup_parameters', message)
     558
     559    CALL log_runtime('time', 'read')
     560
     561    CALL get_cosmo_grid( cfg, soil_files(1), rlon, rlat, hhl, hfl, depths,  &
     562                         d_depth, d_depth_rho_inv, phi_n, lambda_n,         &
     563                         phi_equat,                                         &
     564                         lonmin_cosmo, lonmax_cosmo,                        &
     565                         latmin_cosmo, latmax_cosmo,                        &
     566                         nlon, nlat, nlev, ndepths )
    546567
    547568
    … …  
    550571!------------------------------------------------------------------------------
    551572!
    552 !--    PALM-4U domain extents
    553        lx = (nx+1) * dx
    554        ly = (ny+1) * dy
    555        
    556 !
    557 !--    PALM-4U point of Earth tangency
    558        x0 = 0.0_dp
    559        y0 = 0.0_dp
    560 
    561 !
    562 !--    time vector
    563        nt = CEILING(end_time / (step_hour * 3600.0_dp)) + 1
    564        ALLOCATE( time(nt) )
    565        CALL linspace(0.0_dp, 3600.0_dp * (nt-1), time)
    566        output_file % time => time
    567  CALL run_control('time', 'init')
    568 
    569 !
    570 !--    Convert the PALM-4U origin coordinates to COSMO's rotated-pole grid
    571        phi_c    = TO_RADIANS *                                                 &
    572                   phi2phirot( origin_lat * TO_DEGREES, origin_lon * TO_DEGREES,&
    573                               phi_n * TO_DEGREES, lambda_n * TO_DEGREES )
    574        lambda_c = TO_RADIANS *                                                 &
    575                   rla2rlarot( origin_lat * TO_DEGREES, origin_lon * TO_DEGREES,&
    576                               phi_n * TO_DEGREES, lambda_n * TO_DEGREES,     &
    577                               0.0_dp )
    578 
    579 !
    580 !--    Set gamma according to whether PALM domain is in the northern or southern
    581 !--    hemisphere of the COSMO rotated-pole system. Gamma assumes either the
    582 !--    value 0 or PI and is needed to work around around a bug in the
    583 !--    rotated-pole coordinate transformations.
    584        gam = gamma_from_hemisphere(origin_lat, phi_equat)
    585 
    586 !
    587 !--    Compute the north pole of the rotated-pole grid centred at the PALM-4U
    588 !--    domain centre. The resulting (phi_cn, lambda_cn) are coordinates in
    589 !--    COSMO-DE's rotated-pole grid.
    590        phi_cn    = phic_to_phin(phi_c)
    591        lambda_cn = lamc_to_lamn(phi_c, lambda_c)
    592 
    593        message =   "PALM-4U origin:" // NEW_LINE('') // &
    594           "           lon (lambda) = " // &
    595           TRIM(real_to_str_f(origin_lon * TO_DEGREES)) // " deg"// NEW_LINE(' ') //&
    596           "           lat (phi   ) = " // &
    597           TRIM(real_to_str_f(origin_lat * TO_DEGREES)) // " deg (geographical)" // NEW_LINE(' ') //&
    598           "           lon (lambda) = " // &
    599           TRIM(real_to_str_f(lambda_c * TO_DEGREES)) // " deg" // NEW_LINE(' ') // &
    600           "           lat (phi   ) = " // &
    601           TRIM(real_to_str_f(phi_c * TO_DEGREES)) // " deg (COSMO-DE rotated-pole)"
    602       CALL report ('setup_parameters', message)
    603 
    604        message = "North pole of the rotated COSMO-DE system:" // NEW_LINE(' ') // &
    605           "           lon (lambda) = " // &
    606           TRIM(real_to_str_f(lambda_n * TO_DEGREES)) // " deg" // NEW_LINE(' ') //&
    607           "           lat (phi   ) = " // &
    608           TRIM(real_to_str_f(phi_n * TO_DEGREES)) // " deg (geographical)"
    609        CALL report ('setup_parameters', message)
    610           
    611        message = "North pole of the rotated palm system:" // NEW_LINE(' ') // &
    612           "           lon (lambda) = " // &
    613           TRIM(real_to_str_f(lambda_cn * TO_DEGREES)) // " deg" // NEW_LINE(' ') // &
    614           "           lat (phi   ) = " // &
    615           TRIM(real_to_str_f(phi_cn * TO_DEGREES)) // " deg (COSMO-DE rotated-pole)"
    616        CALL report ('setup_parameters', message)
    617 
    618  CALL run_control('time', 'comp')
     573!-- PALM-4U domain extents
     574    lx = (nx+1) * dx
     575    ly = (ny+1) * dy
     576   
     577!
     578!-- PALM-4U point of Earth tangency
     579    x0 = 0.0_wp
     580    y0 = 0.0_wp
     581
     582!
     583!-- time vector
     584    nt = CEILING(end_time / (step_hour * 3600.0_wp)) + 1
     585    ALLOCATE( time(nt) )
     586    CALL linspace(0.0_wp, 3600.0_wp * (nt-1), time)
     587    output_file%time => time
     588    CALL log_runtime('time', 'init')
     589
     590!
     591!-- Convert the PALM-4U origin coordinates to COSMO's rotated-pole grid
     592    phi_c    = TO_RADIANS *                                                 &
     593               phi2phirot( origin_lat * TO_DEGREES, origin_lon * TO_DEGREES,&
     594                           phi_n * TO_DEGREES, lambda_n * TO_DEGREES )
     595    lambda_c = TO_RADIANS *                                                 &
     596               rla2rlarot( origin_lat * TO_DEGREES, origin_lon * TO_DEGREES,&
     597                           phi_n * TO_DEGREES, lambda_n * TO_DEGREES,     &
     598                           0.0_wp )
     599
     600!
     601!-- Set gamma according to whether PALM domain is in the northern or southern
     602!-- hemisphere of the COSMO rotated-pole system. Gamma assumes either the
     603!-- value 0 or PI and is needed to work around around a bug in the
     604!-- rotated-pole coordinate transformations.
     605    gam = gamma_from_hemisphere(origin_lat, phi_equat)
     606
     607!
     608!-- Compute the north pole of the rotated-pole grid centred at the PALM-4U
     609!-- domain centre. The resulting (phi_cn, lambda_cn) are coordinates in
     610!-- COSMO-DE's rotated-pole grid.
     611    phi_cn    = phic_to_phin(phi_c)
     612    lambda_cn = lamc_to_lamn(phi_c, lambda_c)
     613
     614    message =   "PALM-4U origin:" // NEW_LINE('') // &
     615       "           lon (lambda) = " // &
     616       TRIM(real_to_str_f(origin_lon * TO_DEGREES)) // " deg"// NEW_LINE(' ') //&
     617       "           lat (phi   ) = " // &
     618       TRIM(real_to_str_f(origin_lat * TO_DEGREES)) // " deg (geographical)" // NEW_LINE(' ') //&
     619       "           lon (lambda) = " // &
     620       TRIM(real_to_str_f(lambda_c * TO_DEGREES)) // " deg" // NEW_LINE(' ') // &
     621       "           lat (phi   ) = " // &
     622       TRIM(real_to_str_f(phi_c * TO_DEGREES)) // " deg (COSMO-DE rotated-pole)"
     623    CALL report ('setup_parameters', message)
     624
     625    message = "North pole of the rotated COSMO-DE system:" // NEW_LINE(' ') // &
     626       "           lon (lambda) = " // &
     627       TRIM(real_to_str_f(lambda_n * TO_DEGREES)) // " deg" // NEW_LINE(' ') //&
     628       "           lat (phi   ) = " // &
     629       TRIM(real_to_str_f(phi_n * TO_DEGREES)) // " deg (geographical)"
     630    CALL report ('setup_parameters', message)
     631       
     632    message = "North pole of the rotated palm system:" // NEW_LINE(' ') // &
     633       "           lon (lambda) = " // &
     634       TRIM(real_to_str_f(lambda_cn * TO_DEGREES)) // " deg" // NEW_LINE(' ') // &
     635       "           lat (phi   ) = " // &
     636       TRIM(real_to_str_f(phi_cn * TO_DEGREES)) // " deg (COSMO-DE rotated-pole)"
     637    CALL report ('setup_parameters', message)
     638
     639    CALL log_runtime('time', 'comp')
    619640
    620641!------------------------------------------------------------------------------
    … …  
    625646!-- Compute coordiantes of the PALM centre in the source (COSMO) system
    626647    phi_centre = phirot2phi(                                                   &
    627        phirot = project(0.5_dp*ly, y0, EARTH_RADIUS) * TO_DEGREES,             &
    628        rlarot = project(0.5_dp*lx, x0, EARTH_RADIUS) * TO_DEGREES,             &
     648       phirot = project(0.5_wp*ly, y0, EARTH_RADIUS) * TO_DEGREES,             &
     649       rlarot = project(0.5_wp*lx, x0, EARTH_RADIUS) * TO_DEGREES,             &
    629650       polphi = phi_cn * TO_DEGREES,                                           &
    630651       polgam = gam * TO_DEGREES                                               &
    … …  
    632653
    633654    lam_centre = rlarot2rla(                                                   &
    634        phirot = project(0.5_dp*ly, y0, EARTH_RADIUS) * TO_DEGREES,             &
    635        rlarot = project(0.5_dp*lx, x0, EARTH_RADIUS) * TO_DEGREES,             &
     655       phirot = project(0.5_wp*ly, y0, EARTH_RADIUS) * TO_DEGREES,             &
     656       rlarot = project(0.5_wp*lx, x0, EARTH_RADIUS) * TO_DEGREES,             &
    636657       polphi = phi_cn * TO_DEGREES, pollam = lambda_cn * TO_DEGREES,          &
    637658       polgam = gam * TO_DEGREES                                               &
    … …  
    647668!
    648669!-- Compute boundaries of the central averaging box
    649     averaging_angle = cfg % averaging_angle * TO_RADIANS
    650     lam_east = lam_centre + 0.5_dp * averaging_angle
    651     lam_west = lam_centre - 0.5_dp * averaging_angle
    652     phi_north = phi_centre + 0.5_dp * averaging_angle
    653     phi_south = phi_centre - 0.5_dp * averaging_angle
     670    averaging_angle = cfg%averaging_angle * TO_RADIANS
     671    lam_east = lam_centre + 0.5_wp * averaging_angle
     672    lam_west = lam_centre - 0.5_wp * averaging_angle
     673    phi_north = phi_centre + 0.5_wp * averaging_angle
     674    phi_south = phi_centre - 0.5_wp * averaging_angle
    654675    averaging_width_ew = averaging_angle * COS(phi_centre) * EARTH_RADIUS
    655676    averaging_width_ns = averaging_angle * EARTH_RADIUS
    … …  
    674695!
    675696!-- Coriolis parameter
    676     f3 = 2.0_dp * OMEGA * SIN(                                                 &
     697    f3 = 2.0_wp * OMEGA * SIN(                                                 &
    677698       TO_RADIANS*phirot2phi( phi_centre * TO_DEGREES, lam_centre * TO_DEGREES,&
    678699                              phi_n * TO_DEGREES,                              &
    … …  
    680701    )
    681702
    682     END SUBROUTINE setup_parameters
     703 END SUBROUTINE setup_parameters
    683704
    684705
    … …  
    689710!> coordinates and interpolation weights
    690711!------------------------------------------------------------------------------!
    691     SUBROUTINE setup_grids()
    692        CHARACTER ::  interp_mode
     712 SUBROUTINE setup_grids()
     713    CHARACTER ::  interp_mode
    693714
    694715!------------------------------------------------------------------------------
    … …  
    696717!------------------------------------------------------------------------------
    697718!
    698 !--    palm x y z, we allocate the column to nz+1 in order to include the top
    699 !--    scalar boundary. The interpolation grids will be associated with
    700 !--    a shorter column that omits the top element.
    701        ALLOCATE( x(0:nx), y(0:ny), z(1:nz), z_column(1:nz+1) )
    702        CALL linspace(0.5_dp * dx, lx - 0.5_dp * dx, x)
    703        CALL linspace(0.5_dp * dy, ly - 0.5_dp * dy, y)
    704        CALL stretched_z(z_column, dz, dz_max=dz_max,                           &
    705                         dz_stretch_factor=dz_stretch_factor,                   &
    706                         dz_stretch_level=dz_stretch_level,                     &
    707                         dz_stretch_level_start=dz_stretch_level_start,         &
    708                         dz_stretch_level_end=dz_stretch_level_end,             &
    709                         dz_stretch_factor_array=dz_stretch_factor_array)
    710        z(1:nz) = z_column(1:nz)
    711        z_top = z_column(nz+1)
    712 
    713 !
    714 !--    palm xu yv zw, compared to the scalar grid, velocity coordinates
    715 !--    contain one element less.
    716        ALLOCATE( xu(1:nx),  yv(1:ny), zw(1:nz-1), zw_column(1:nz))
    717        CALL linspace(dx, lx - dx, xu)
    718        CALL linspace(dy, ly - dy, yv)
    719        CALL midpoints(z_column, zw_column)
    720        zw(1:nz-1) = zw_column(1:nz-1)
    721        zw_top     = zw_column(nz)
     719!-- palm x y z, we allocate the column to nz+1 in order to include the top
     720!-- scalar boundary. The interpolation grids will be associated with
     721!-- a shorter column that omits the top element.
     722    ALLOCATE( x(0:nx), y(0:ny), z(1:nz), z_column(1:nz+1) )
     723    CALL linspace(0.5_wp * dx, lx - 0.5_wp * dx, x)
     724    CALL linspace(0.5_wp * dy, ly - 0.5_wp * dy, y)
     725    CALL stretched_z(z_column, dz, dz_max=dz_max,                           &
     726                     dz_stretch_factor=dz_stretch_factor,                   &
     727                     dz_stretch_level=dz_stretch_level,                     &
     728                     dz_stretch_level_start=dz_stretch_level_start,         &
     729                     dz_stretch_level_end=dz_stretch_level_end,             &
     730                     dz_stretch_factor_array=dz_stretch_factor_array)
     731    z(1:nz) = z_column(1:nz)
     732    z_top = z_column(nz+1)
     733
     734!
     735!-- palm xu yv zw, compared to the scalar grid, velocity coordinates
     736!-- contain one element less.
     737    ALLOCATE( xu(1:nx),  yv(1:ny), zw(1:nz-1), zw_column(1:nz))
     738    CALL linspace(dx, lx - dx, xu)
     739    CALL linspace(dy, ly - dy, yv)
     740    CALL midpoints(z_column, zw_column)
     741    zw(1:nz-1) = zw_column(1:nz-1)
     742    zw_top     = zw_column(nz)
    722743
    723744
    … …  
    725746! Section 1: Define initialization and boundary grids
    726747!------------------------------------------------------------------------------
    727        CALL init_grid_definition('palm', grid=palm_grid,                       &
    728                xmin=0.0_dp, xmax=lx,                                           &
    729                ymin=0.0_dp, ymax=ly,                                           &
    730                x0=x0, y0=y0, z0=z0,                                            &
    731                nx=nx, ny=ny, nz=nz, z=z, zw=zw, ic_mode=cfg % ic_mode)
    732 
    733 !
    734 !--    Subtracting 1 because arrays will be allocated with nlon + 1 elements.
    735        CALL init_grid_definition('cosmo-de', grid=cosmo_grid,                  &
    736                xmin=lonmin_cosmo, xmax=lonmax_cosmo,                           &
    737                ymin=latmin_cosmo, ymax=latmax_cosmo,                           &
    738                x0=x0, y0=y0, z0=0.0_dp,                                        &
    739                nx=nlon-1, ny=nlat-1, nz=nlev-1)
    740 
    741 !
    742 !--    Define intermediate grid. This is the same as palm_grid except with a
    743 !--    much coarser vertical grid. The vertical levels are interpolated in each
    744 !--    PALM column from COSMO's secondary levels. The main levels are then
    745 !--    computed as the averages of the bounding secondary levels.
    746        CALL init_grid_definition('palm intermediate', grid=palm_intermediate,  &
    747                xmin=0.0_dp, xmax=lx,                                           &
    748                ymin=0.0_dp, ymax=ly,                                           &
    749                x0=x0, y0=y0, z0=z0,                                            &
    750                nx=nx, ny=ny, nz=nlev-2)
    751 
    752        CALL init_grid_definition('boundary', grid=u_initial_grid,              &
     748    CALL init_grid_definition('palm', grid=palm_grid,                       &
     749            xmin=0.0_wp, xmax=lx,                                           &
     750            ymin=0.0_wp, ymax=ly,                                           &
     751            x0=x0, y0=y0, z0=z0,                                            &
     752            nx=nx, ny=ny, nz=nz, z=z, zw=zw, ic_mode=cfg%ic_mode)
     753
     754!
     755!-- Subtracting 1 because arrays will be allocated with nlon + 1 elements.
     756    CALL init_grid_definition('cosmo-de', grid=cosmo_grid,                  &
     757            xmin=lonmin_cosmo, xmax=lonmax_cosmo,                           &
     758            ymin=latmin_cosmo, ymax=latmax_cosmo,                           &
     759            x0=x0, y0=y0, z0=0.0_wp,                                        &
     760            nx=nlon-1, ny=nlat-1, nz=nlev-1)
     761
     762!
     763!-- Define intermediate grid. This is the same as palm_grid except with a
     764!-- much coarser vertical grid. The vertical levels are interpolated in each
     765!-- PALM column from COSMO's secondary levels. The main levels are then
     766!-- computed as the averages of the bounding secondary levels.
     767    CALL init_grid_definition('palm intermediate', grid=palm_intermediate,  &
     768            xmin=0.0_wp, xmax=lx,                                           &
     769            ymin=0.0_wp, ymax=ly,                                           &
     770            x0=x0, y0=y0, z0=z0,                                            &
     771            nx=nx, ny=ny, nz=nlev-2)
     772
     773    CALL init_grid_definition('boundary', grid=u_initial_grid,              &
     774            xmin = dx, xmax = lx - dx,                                      &
     775            ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy,                    &
     776            x0=x0, y0=y0, z0 = z0,                                          &
     777            nx = nx-1, ny = ny, nz = nz,                                    &
     778            z=z, ic_mode=cfg%ic_mode)
     779
     780    CALL init_grid_definition('boundary', grid=v_initial_grid,              &
     781            xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx,                    &
     782            ymin = dy, ymax = ly - dy,                                      &
     783            x0=x0, y0=y0, z0 = z0,                                          &
     784            nx = nx, ny = ny-1, nz = nz,                                    &
     785            z=z, ic_mode=cfg%ic_mode)
     786
     787    CALL init_grid_definition('boundary', grid=w_initial_grid,              &
     788            xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx,                    &
     789            ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy,                    &
     790            x0=x0, y0=y0, z0 = z0,                                          &
     791            nx = nx, ny = ny, nz = nz-1,                                    &
     792            z=zw, ic_mode=cfg%ic_mode)
     793
     794    CALL init_grid_definition('boundary intermediate', grid=u_initial_intermediate,      &
     795            xmin = dx, xmax = lx - dx,                                      &
     796            ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy,                    &
     797            x0=x0, y0=y0, z0 = z0,                                          &
     798            nx = nx-1, ny = ny, nz = nlev - 2)
     799
     800    CALL init_grid_definition('boundary intermediate', grid=v_initial_intermediate,      &
     801            xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx,                    &
     802            ymin = dy, ymax = ly - dy,                                      &
     803            x0=x0, y0=y0, z0 = z0,                                          &
     804            nx = nx, ny = ny-1, nz = nlev - 2)
     805
     806    CALL init_grid_definition('boundary intermediate', grid=w_initial_intermediate,      &
     807            xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx,                    &
     808            ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy,                    &
     809            x0=x0, y0=y0, z0 = z0,                                          &
     810            nx = nx, ny = ny, nz = nlev - 1)
     811
     812    IF (boundary_variables_required)  THEN
     813!
     814!------------------------------------------------------------------------------
     815! Section 2: Define PALM-4U boundary grids
     816!------------------------------------------------------------------------------
     817       CALL init_grid_definition('boundary', grid=scalars_east_grid,           &
     818               xmin = lx + 0.5_wp * dx, xmax = lx + 0.5_wp * dx,               &
     819               ymin =  0.5_wp * dy, ymax = ly - 0.5_wp * dy,                   &
     820               x0=x0, y0=y0, z0 = z0,                                          &
     821               nx = 0, ny = ny, nz = nz, z=z)
     822
     823       CALL init_grid_definition('boundary', grid=scalars_west_grid,           &
     824               xmin = -0.5_wp * dx, xmax = -0.5_wp * dx,                       &
     825               ymin =  0.5_wp * dy, ymax = ly - 0.5_wp * dy,                   &
     826               x0=x0, y0=y0, z0 = z0,                                          &
     827               nx = 0, ny = ny, nz = nz, z=z)
     828
     829       CALL init_grid_definition('boundary', grid=scalars_north_grid,          &
     830               xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx,                    &
     831               ymin = ly + 0.5_wp * dy, ymax = ly + 0.5_wp * dy,               &
     832               x0=x0, y0=y0, z0 = z0,                                          &
     833               nx = nx, ny = 0, nz = nz, z=z)
     834
     835       CALL init_grid_definition('boundary', grid=scalars_south_grid,          &
     836               xmin =  0.5_wp * dx, xmax = lx - 0.5_wp * dx,                   &
     837               ymin = -0.5_wp * dy, ymax = -0.5_wp * dy,                       &
     838               x0=x0, y0=y0, z0 = z0,                                          &
     839               nx = nx, ny = 0, nz = nz, z=z)
     840
     841       CALL init_grid_definition('boundary', grid=scalars_top_grid,            &
     842               xmin =  0.5_wp * dx, xmax = lx - 0.5_wp * dx,                   &
     843               ymin =  0.5_wp * dy, ymax = ly - 0.5_wp * dy,                   &
     844               x0=x0, y0=y0, z0 = z0,                                          &
     845               nx = nx, ny = ny, nz = 1, z=(/z_top/))
     846
     847       CALL init_grid_definition('boundary', grid=u_east_grid,                 &
     848               xmin = lx, xmax = lx,                                           &
     849               ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy,                    &
     850               x0=x0, y0=y0, z0 = z0,                                          &
     851               nx = 0, ny = ny, nz = nz, z=z)
     852
     853       CALL init_grid_definition('boundary', grid=u_west_grid,                 &
     854               xmin = 0.0_wp, xmax = 0.0_wp,                                   &
     855               ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy,                    &
     856               x0=x0, y0=y0, z0 = z0,                                          &
     857               nx = 0, ny = ny, nz = nz, z=z)
     858
     859       CALL init_grid_definition('boundary', grid=u_north_grid,                &
    753860               xmin = dx, xmax = lx - dx,                                      &
    754                ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy,                    &
     861               ymin = ly + 0.5_wp * dy, ymax = ly + 0.5_wp * dy,               &
    755862               x0=x0, y0=y0, z0 = z0,                                          &
    756                nx = nx-1, ny = ny, nz = nz,                                    &
    757                z=z, ic_mode=cfg % ic_mode)
    758 
    759        CALL init_grid_definition('boundary', grid=v_initial_grid,              &
    760                xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
     863               nx = nx-1, ny = 0, nz = nz, z=z)
     864   
     865       CALL init_grid_definition('boundary', grid=u_south_grid,                &
     866               xmin = dx, xmax = lx - dx,                                      &
     867               ymin = -0.5_wp * dy, ymax = -0.5_wp * dy,                       &
     868               x0=x0, y0=y0, z0 = z0,                                          &
     869               nx = nx-1, ny = 0, nz = nz, z=z)
     870
     871       CALL init_grid_definition('boundary', grid=u_top_grid,                  &
     872               xmin = dx, xmax = lx - dx,                                      &
     873               ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy,                    &
     874               x0=x0, y0=y0, z0 = z0,                                          &
     875               nx = nx-1, ny = ny, nz = 1, z=(/z_top/))
     876
     877       CALL init_grid_definition('boundary', grid=v_east_grid,                 &
     878               xmin = lx + 0.5_wp * dx, xmax = lx + 0.5_wp * dx,               &
    761879               ymin = dy, ymax = ly - dy,                                      &
    762880               x0=x0, y0=y0, z0 = z0,                                          &
    763                nx = nx, ny = ny-1, nz = nz,                                    &
    764                z=z, ic_mode=cfg % ic_mode)
    765 
    766        CALL init_grid_definition('boundary', grid=w_initial_grid,              &
    767                xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
    768                ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy,                    &
     881               nx = 0, ny = ny-1, nz = nz, z=z)
     882
     883       CALL init_grid_definition('boundary', grid=v_west_grid,                 &
     884               xmin = -0.5_wp * dx, xmax = -0.5_wp * dx,                       &
     885               ymin = dy, ymax = ly - dy,                                      &
    769886               x0=x0, y0=y0, z0 = z0,                                          &
    770                nx = nx, ny = ny, nz = nz-1,                                    &
    771                z=zw, ic_mode=cfg % ic_mode)
    772 
    773        CALL init_grid_definition('boundary intermediate', grid=u_initial_intermediate,      &
     887               nx = 0, ny = ny-1, nz = nz, z=z)
     888
     889       CALL init_grid_definition('boundary', grid=v_north_grid,                &
     890               xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx,                    &
     891               ymin = ly, ymax = ly,                                           &
     892               x0=x0, y0=y0, z0 = z0,                                          &
     893               nx = nx, ny = 0, nz = nz, z=z)
     894
     895       CALL init_grid_definition('boundary', grid=v_south_grid,                &
     896               xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx,                    &
     897               ymin = 0.0_wp, ymax = 0.0_wp,                                   &
     898               x0=x0, y0=y0, z0 = z0,                                          &
     899               nx = nx, ny = 0, nz = nz, z=z)
     900
     901       CALL init_grid_definition('boundary', grid=v_top_grid,                  &
     902               xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx,                    &
     903               ymin = dy, ymax = ly - dy,                                      &
     904               x0=x0, y0=y0, z0 = z0,                                          &
     905               nx = nx, ny = ny-1, nz = 1, z=(/z_top/))
     906
     907       CALL init_grid_definition('boundary', grid=w_east_grid,                 &
     908               xmin = lx + 0.5_wp * dx, xmax = lx + 0.5_wp * dx,               &
     909               ymin =  0.5_wp * dy, ymax = ly - 0.5_wp * dy,                   &
     910               x0=x0, y0=y0, z0 = z0,                                          &
     911               nx = 0, ny = ny, nz = nz - 1, z=zw)
     912
     913       CALL init_grid_definition('boundary', grid=w_west_grid,                 &
     914               xmin = -0.5_wp * dx, xmax = -0.5_wp * dx,                       &
     915               ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy,                    &
     916               x0=x0, y0=y0, z0 = z0,                                          &
     917               nx = 0, ny = ny, nz = nz - 1, z=zw)
     918
     919       CALL init_grid_definition('boundary', grid=w_north_grid,                &
     920               xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx,                    &
     921               ymin = ly + 0.5_wp * dy, ymax = ly + 0.5_wp * dy,               &
     922               x0=x0, y0=y0, z0 = z0,                                          &
     923               nx = nx, ny = 0, nz = nz - 1, z=zw)
     924
     925       CALL init_grid_definition('boundary', grid=w_south_grid,                &
     926               xmin =  0.5_wp * dx, xmax = lx - 0.5_wp * dx,                   &
     927               ymin = -0.5_wp * dy, ymax = -0.5_wp * dy,                       &
     928               x0=x0, y0=y0, z0 = z0,                                          &
     929               nx = nx, ny = 0, nz = nz - 1, z=zw)
     930
     931       CALL init_grid_definition('boundary', grid=w_top_grid,                  &
     932               xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx,                    &
     933               ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy,                    &
     934               x0=x0, y0=y0, z0 = z0,                                          &
     935               nx = nx, ny = ny, nz = 1, z=(/zw_top/))
     936
     937       CALL init_grid_definition('boundary intermediate', grid=scalars_east_intermediate,   &
     938               xmin = lx + 0.5_wp * dx, xmax = lx + 0.5_wp * dx,               &
     939               ymin =  0.5_wp * dy, ymax = ly - 0.5_wp * dy,                   &
     940               x0=x0, y0=y0, z0 = z0,                                          &
     941               nx = 0, ny = ny, nz = nlev - 2)
     942
     943       CALL init_grid_definition('boundary intermediate', grid=scalars_west_intermediate,   &
     944               xmin = -0.5_wp * dx, xmax = -0.5_wp * dx,                       &
     945               ymin =  0.5_wp * dy, ymax = ly - 0.5_wp * dy,                   &
     946               x0=x0, y0=y0, z0 = z0,                                          &
     947               nx = 0, ny = ny, nz = nlev - 2)
     948
     949       CALL init_grid_definition('boundary intermediate', grid=scalars_north_intermediate,  &
     950               xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx,                    &
     951               ymin = ly + 0.5_wp * dy, ymax = ly + 0.5_wp * dy,               &
     952               x0=x0, y0=y0, z0 = z0,                                          &
     953               nx = nx, ny = 0, nz = nlev - 2)
     954
     955       CALL init_grid_definition('boundary intermediate', grid=scalars_south_intermediate,  &
     956               xmin =  0.5_wp * dx, xmax = lx - 0.5_wp * dx,                   &
     957               ymin = -0.5_wp * dy, ymax = -0.5_wp * dy,                       &
     958               x0=x0, y0=y0, z0 = z0,                                          &
     959               nx = nx, ny = 0, nz = nlev - 2)
     960
     961       CALL init_grid_definition('boundary intermediate', grid=scalars_top_intermediate,    &
     962               xmin =  0.5_wp * dx, xmax = lx - 0.5_wp * dx,                   &
     963               ymin =  0.5_wp * dy, ymax = ly - 0.5_wp * dy,                   &
     964               x0=x0, y0=y0, z0 = z0,                                          &
     965               nx = nx, ny = ny, nz = nlev - 2)
     966
     967       CALL init_grid_definition('boundary intermediate', grid=u_east_intermediate,         &
     968               xmin = lx, xmax = lx,                                           &
     969               ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy,                    &
     970               x0=x0, y0=y0, z0 = z0,                                          &
     971               nx = 0, ny = ny, nz = nlev - 2)
     972
     973       CALL init_grid_definition('boundary intermediate', grid=u_west_intermediate,         &
     974               xmin = 0.0_wp, xmax = 0.0_wp,                                   &
     975               ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy,                    &
     976               x0=x0, y0=y0, z0 = z0,                                          &
     977               nx = 0, ny = ny, nz = nlev - 2)
     978
     979       CALL init_grid_definition('boundary intermediate', grid=u_north_intermediate,        &
    774980               xmin = dx, xmax = lx - dx,                                      &
    775                ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy,                    &
     981               ymin = ly + 0.5_wp * dy, ymax = ly + 0.5_wp * dy,               &
     982               x0=x0, y0=y0, z0 = z0,                                          &
     983               nx = nx-1, ny = 0, nz = nlev - 2)
     984
     985       CALL init_grid_definition('boundary intermediate', grid=u_south_intermediate,        &
     986               xmin = dx, xmax = lx - dx,                                      &
     987               ymin = -0.5_wp * dy, ymax = -0.5_wp * dy,                       &
     988               x0=x0, y0=y0, z0 = z0,                                          &
     989               nx = nx-1, ny = 0, nz = nlev - 2)
     990
     991       CALL init_grid_definition('boundary intermediate', grid=u_top_intermediate,          &
     992               xmin = dx, xmax = lx - dx,                                      &
     993               ymin = 0.5_wp * dy, ymax = ly - 0.5_wp * dy,                    &
    776994               x0=x0, y0=y0, z0 = z0,                                          &
    777995               nx = nx-1, ny = ny, nz = nlev - 2)
    778996
    779        CALL init_grid_definition('boundary intermediate', grid=v_initial_intermediate,      &
    780                xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
     997       CALL init_grid_definition('boundary intermediate', grid=v_east_intermediate,         &
     998               xmin = lx + 0.5_wp * dx, xmax = lx + 0.5_wp * dx,               &
     999               ymin = dy, ymax = ly - dy,                                      &
     1000               x0=x0, y0=y0, z0 = z0,                                          &
     1001               nx = 0, ny = ny-1, nz = nlev - 2)
     1002
     1003       CALL init_grid_definition('boundary intermediate', grid=v_west_intermediate,         &
     1004               xmin = -0.5_wp * dx, xmax = -0.5_wp * dx,                       &
     1005               ymin = dy, ymax = ly - dy,                                      &
     1006               x0=x0, y0=y0, z0 = z0,                                          &
     1007               nx = 0, ny = ny-1, nz = nlev - 2)
     1008
     1009       CALL init_grid_definition('boundary intermediate', grid=v_north_intermediate,        &
     1010               xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx,                    &
     1011               ymin = ly, ymax = ly,                                           &
     1012               x0=x0, y0=y0, z0 = z0,                                          &
     1013               nx = nx, ny = 0, nz = nlev - 2)
     1014
     1015       CALL init_grid_definition('boundary intermediate', grid=v_south_intermediate,        &
     1016               xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx,                    &
     1017               ymin = 0.0_wp, ymax = 0.0_wp,                                   &
     1018               x0=x0, y0=y0, z0 = z0,                                          &
     1019               nx = nx, ny = 0, nz = nlev - 2)
     1020
     1021       CALL init_grid_definition('boundary intermediate', grid=v_top_intermediate,          &
     1022               xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx,                    &
    7811023               ymin = dy, ymax = ly - dy,                                      &
    7821024               x0=x0, y0=y0, z0 = z0,                                          &
    7831025               nx = nx, ny = ny-1, nz = nlev - 2)
    7841026
    785        CALL init_grid_definition('boundary intermediate', grid=w_initial_intermediate,      &
    786                xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
    787                ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy,                    &
     1027       CALL init_grid_definition('boundary intermediate', grid=w_east_intermediate,         &
     1028               xmin = lx + 0.5_wp * dx, xmax = lx + 0.5_wp * dx,               &
     1029               ymin =  0.5_wp * dy, ymax = ly - 0.5_wp * dy,                   &
     1030               x0=x0, y0=y0, z0 = z0,                                          &
     1031               nx = 0, ny = ny, nz = nlev - 1)
     1032
     1033       CALL init_grid_definition('boundary intermediate', grid=w_west_intermediate,         &
     1034               xmin = -0.5_wp * dx, xmax = -0.5_wp * dx,                       &
     1035               ymin =  0.5_wp * dy, ymax = ly - 0.5_wp * dy,                   &
     1036               x0=x0, y0=y0, z0 = z0,                                          &
     1037               nx = 0, ny = ny, nz = nlev - 1)
     1038
     1039       CALL init_grid_definition('boundary intermediate', grid=w_north_intermediate,        &
     1040               xmin = 0.5_wp * dx, xmax = lx - 0.5_wp * dx,                    &
     1041               ymin = ly + 0.5_wp * dy, ymax = ly + 0.5_wp * dy,               &
     1042               x0=x0, y0=y0, z0 = z0,                                          &
     1043               nx = nx, ny = 0, nz = nlev - 1)
     1044
     1045       CALL init_grid_definition('boundary intermediate', grid=w_south_intermediate,        &
     1046               xmin =  0.5_wp * dx, xmax = lx - 0.5_wp * dx,                   &
     1047               ymin = -0.5_wp * dy, ymax = -0.5_wp * dy,                       &
     1048               x0=x0, y0=y0, z0 = z0,                                          &
     1049               nx = nx, ny = 0, nz = nlev - 1)
     1050
     1051       CALL init_grid_definition('boundary intermediate', grid=w_top_intermediate,          &
     1052               xmin =  0.5_wp * dx, xmax = lx - 0.5_wp * dx,                   &
     1053               ymin =  0.5_wp * dy, ymax = ly - 0.5_wp * dy,                   &
    7881054               x0=x0, y0=y0, z0 = z0,                                          &
    7891055               nx = nx, ny = ny, nz = nlev - 1)
    790 
    791       IF (boundary_variables_required)  THEN
    792 !
    793 !------------------------------------------------------------------------------
    794 ! Section 2: Define PALM-4U boundary grids
    795 !------------------------------------------------------------------------------
    796           CALL init_grid_definition('boundary', grid=scalars_east_grid,           &
    797                   xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx,               &
    798                   ymin =  0.5_dp * dy, ymax = ly - 0.5_dp * dy,                   &
    799                   x0=x0, y0=y0, z0 = z0,                                          &
    800                   nx = 0, ny = ny, nz = nz, z=z)
    801 
    802           CALL init_grid_definition('boundary', grid=scalars_west_grid,           &
    803                   xmin = -0.5_dp * dx, xmax = -0.5_dp * dx,                       &
    804                   ymin =  0.5_dp * dy, ymax = ly - 0.5_dp * dy,                   &
    805                   x0=x0, y0=y0, z0 = z0,                                          &
    806                   nx = 0, ny = ny, nz = nz, z=z)
    807 
    808           CALL init_grid_definition('boundary', grid=scalars_north_grid,          &
    809                   xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
    810                   ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy,               &
    811                   x0=x0, y0=y0, z0 = z0,                                          &
    812                   nx = nx, ny = 0, nz = nz, z=z)
    813 
    814           CALL init_grid_definition('boundary', grid=scalars_south_grid,          &
    815                   xmin =  0.5_dp * dx, xmax = lx - 0.5_dp * dx,                   &
    816                   ymin = -0.5_dp * dy, ymax = -0.5_dp * dy,                       &
    817                   x0=x0, y0=y0, z0 = z0,                                          &
    818                   nx = nx, ny = 0, nz = nz, z=z)
    819 
    820           CALL init_grid_definition('boundary', grid=scalars_top_grid,            &
    821                   xmin =  0.5_dp * dx, xmax = lx - 0.5_dp * dx,                   &
    822                   ymin =  0.5_dp * dy, ymax = ly - 0.5_dp * dy,                   &
    823                   x0=x0, y0=y0, z0 = z0,                                          &
    824                   nx = nx, ny = ny, nz = 1, z=(/z_top/))
    825 
    826           CALL init_grid_definition('boundary', grid=u_east_grid,                 &
    827                   xmin = lx, xmax = lx,                                           &
    828                   ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy,                    &
    829                   x0=x0, y0=y0, z0 = z0,                                          &
    830                   nx = 0, ny = ny, nz = nz, z=z)
    831 
    832           CALL init_grid_definition('boundary', grid=u_west_grid,                 &
    833                   xmin = 0.0_dp, xmax = 0.0_dp,                                   &
    834                   ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy,                    &
    835                   x0=x0, y0=y0, z0 = z0,                                          &
    836                   nx = 0, ny = ny, nz = nz, z=z)
    837 
    838           CALL init_grid_definition('boundary', grid=u_north_grid,                &
    839                   xmin = dx, xmax = lx - dx,                                      &
    840                   ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy,               &
    841                   x0=x0, y0=y0, z0 = z0,                                          &
    842                   nx = nx-1, ny = 0, nz = nz, z=z)
    843    
    844           CALL init_grid_definition('boundary', grid=u_south_grid,                &
    845                   xmin = dx, xmax = lx - dx,                                      &
    846                   ymin = -0.5_dp * dy, ymax = -0.5_dp * dy,                       &
    847                   x0=x0, y0=y0, z0 = z0,                                          &
    848                   nx = nx-1, ny = 0, nz = nz, z=z)
    849 
    850           CALL init_grid_definition('boundary', grid=u_top_grid,                  &
    851                   xmin = dx, xmax = lx - dx,                                      &
    852                   ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy,                    &
    853                   x0=x0, y0=y0, z0 = z0,                                          &
    854                   nx = nx-1, ny = ny, nz = 1, z=(/z_top/))
    855 
    856           CALL init_grid_definition('boundary', grid=v_east_grid,                 &
    857                   xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx,               &
    858                   ymin = dy, ymax = ly - dy,                                      &
    859                   x0=x0, y0=y0, z0 = z0,                                          &
    860                   nx = 0, ny = ny-1, nz = nz, z=z)
    861 
    862           CALL init_grid_definition('boundary', grid=v_west_grid,                 &
    863                   xmin = -0.5_dp * dx, xmax = -0.5_dp * dx,                       &
    864                   ymin = dy, ymax = ly - dy,                                      &
    865                   x0=x0, y0=y0, z0 = z0,                                          &
    866                   nx = 0, ny = ny-1, nz = nz, z=z)
    867 
    868           CALL init_grid_definition('boundary', grid=v_north_grid,                &
    869                   xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
    870                   ymin = ly, ymax = ly,                                           &
    871                   x0=x0, y0=y0, z0 = z0,                                          &
    872                   nx = nx, ny = 0, nz = nz, z=z)
    873 
    874           CALL init_grid_definition('boundary', grid=v_south_grid,                &
    875                   xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
    876                   ymin = 0.0_dp, ymax = 0.0_dp,                                   &
    877                   x0=x0, y0=y0, z0 = z0,                                          &
    878                   nx = nx, ny = 0, nz = nz, z=z)
    879 
    880           CALL init_grid_definition('boundary', grid=v_top_grid,                  &
    881                   xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
    882                   ymin = dy, ymax = ly - dy,                                      &
    883                   x0=x0, y0=y0, z0 = z0,                                          &
    884                   nx = nx, ny = ny-1, nz = 1, z=(/z_top/))
    885 
    886           CALL init_grid_definition('boundary', grid=w_east_grid,                 &
    887                   xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx,               &
    888                   ymin =  0.5_dp * dy, ymax = ly - 0.5_dp * dy,                   &
    889                   x0=x0, y0=y0, z0 = z0,                                          &
    890                   nx = 0, ny = ny, nz = nz - 1, z=zw)
    891 
    892           CALL init_grid_definition('boundary', grid=w_west_grid,                 &
    893                   xmin = -0.5_dp * dx, xmax = -0.5_dp * dx,                       &
    894                   ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy,                    &
    895                   x0=x0, y0=y0, z0 = z0,                                          &
    896                   nx = 0, ny = ny, nz = nz - 1, z=zw)
    897 
    898           CALL init_grid_definition('boundary', grid=w_north_grid,                &
    899                   xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
    900                   ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy,               &
    901                   x0=x0, y0=y0, z0 = z0,                                          &
    902                   nx = nx, ny = 0, nz = nz - 1, z=zw)
    903 
    904           CALL init_grid_definition('boundary', grid=w_south_grid,                &
    905                   xmin =  0.5_dp * dx, xmax = lx - 0.5_dp * dx,                   &
    906                   ymin = -0.5_dp * dy, ymax = -0.5_dp * dy,                       &
    907                   x0=x0, y0=y0, z0 = z0,                                          &
    908                   nx = nx, ny = 0, nz = nz - 1, z=zw)
    909 
    910           CALL init_grid_definition('boundary', grid=w_top_grid,                  &
    911                   xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
    912                   ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy,                    &
    913                   x0=x0, y0=y0, z0 = z0,                                          &
    914                   nx = nx, ny = ny, nz = 1, z=(/zw_top/))
    915 
    916           CALL init_grid_definition('boundary intermediate', grid=scalars_east_intermediate,   &
    917                   xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx,               &
    918                   ymin =  0.5_dp * dy, ymax = ly - 0.5_dp * dy,                   &
    919                   x0=x0, y0=y0, z0 = z0,                                          &
    920                   nx = 0, ny = ny, nz = nlev - 2)
    921 
    922           CALL init_grid_definition('boundary intermediate', grid=scalars_west_intermediate,   &
    923                   xmin = -0.5_dp * dx, xmax = -0.5_dp * dx,                       &
    924                   ymin =  0.5_dp * dy, ymax = ly - 0.5_dp * dy,                   &
    925                   x0=x0, y0=y0, z0 = z0,                                          &
    926                   nx = 0, ny = ny, nz = nlev - 2)
    927 
    928           CALL init_grid_definition('boundary intermediate', grid=scalars_north_intermediate,  &
    929                   xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
    930                   ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy,               &
    931                   x0=x0, y0=y0, z0 = z0,                                          &
    932                   nx = nx, ny = 0, nz = nlev - 2)
    933 
    934           CALL init_grid_definition('boundary intermediate', grid=scalars_south_intermediate,  &
    935                   xmin =  0.5_dp * dx, xmax = lx - 0.5_dp * dx,                   &
    936                   ymin = -0.5_dp * dy, ymax = -0.5_dp * dy,                       &
    937                   x0=x0, y0=y0, z0 = z0,                                          &
    938                   nx = nx, ny = 0, nz = nlev - 2)
    939 
    940           CALL init_grid_definition('boundary intermediate', grid=scalars_top_intermediate,    &
    941                   xmin =  0.5_dp * dx, xmax = lx - 0.5_dp * dx,                   &
    942                   ymin =  0.5_dp * dy, ymax = ly - 0.5_dp * dy,                   &
    943                   x0=x0, y0=y0, z0 = z0,                                          &
    944                   nx = nx, ny = ny, nz = nlev - 2)
    945 
    946           CALL init_grid_definition('boundary intermediate', grid=u_east_intermediate,         &
    947                   xmin = lx, xmax = lx,                                           &
    948                   ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy,                    &
    949                   x0=x0, y0=y0, z0 = z0,                                          &
    950                   nx = 0, ny = ny, nz = nlev - 2)
    951 
    952           CALL init_grid_definition('boundary intermediate', grid=u_west_intermediate,         &
    953                   xmin = 0.0_dp, xmax = 0.0_dp,                                   &
    954                   ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy,                    &
    955                   x0=x0, y0=y0, z0 = z0,                                          &
    956                   nx = 0, ny = ny, nz = nlev - 2)
    957 
    958           CALL init_grid_definition('boundary intermediate', grid=u_north_intermediate,        &
    959                   xmin = dx, xmax = lx - dx,                                      &
    960                   ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy,               &
    961                   x0=x0, y0=y0, z0 = z0,                                          &
    962                   nx = nx-1, ny = 0, nz = nlev - 2)
    963 
    964           CALL init_grid_definition('boundary intermediate', grid=u_south_intermediate,        &
    965                   xmin = dx, xmax = lx - dx,                                      &
    966                   ymin = -0.5_dp * dy, ymax = -0.5_dp * dy,                       &
    967                   x0=x0, y0=y0, z0 = z0,                                          &
    968                   nx = nx-1, ny = 0, nz = nlev - 2)
    969 
    970           CALL init_grid_definition('boundary intermediate', grid=u_top_intermediate,          &
    971                   xmin = dx, xmax = lx - dx,                                      &
    972                   ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy,                    &
    973                   x0=x0, y0=y0, z0 = z0,                                          &
    974                   nx = nx-1, ny = ny, nz = nlev - 2)
    975 
    976           CALL init_grid_definition('boundary intermediate', grid=v_east_intermediate,         &
    977                   xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx,               &
    978                   ymin = dy, ymax = ly - dy,                                      &
    979                   x0=x0, y0=y0, z0 = z0,                                          &
    980                   nx = 0, ny = ny-1, nz = nlev - 2)
    981 
    982           CALL init_grid_definition('boundary intermediate', grid=v_west_intermediate,         &
    983                   xmin = -0.5_dp * dx, xmax = -0.5_dp * dx,                       &
    984                   ymin = dy, ymax = ly - dy,                                      &
    985                   x0=x0, y0=y0, z0 = z0,                                          &
    986                   nx = 0, ny = ny-1, nz = nlev - 2)
    987 
    988           CALL init_grid_definition('boundary intermediate', grid=v_north_intermediate,        &
    989                   xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
    990                   ymin = ly, ymax = ly,                                           &
    991                   x0=x0, y0=y0, z0 = z0,                                          &
    992                   nx = nx, ny = 0, nz = nlev - 2)
    993 
    994           CALL init_grid_definition('boundary intermediate', grid=v_south_intermediate,        &
    995                   xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
    996                   ymin = 0.0_dp, ymax = 0.0_dp,                                   &
    997                   x0=x0, y0=y0, z0 = z0,                                          &
    998                   nx = nx, ny = 0, nz = nlev - 2)
    999 
    1000           CALL init_grid_definition('boundary intermediate', grid=v_top_intermediate,          &
    1001                   xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
    1002                   ymin = dy, ymax = ly - dy,                                      &
    1003                   x0=x0, y0=y0, z0 = z0,                                          &
    1004                   nx = nx, ny = ny-1, nz = nlev - 2)
    1005 
    1006           CALL init_grid_definition('boundary intermediate', grid=w_east_intermediate,         &
    1007                   xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx,               &
    1008                   ymin =  0.5_dp * dy, ymax = ly - 0.5_dp * dy,                   &
    1009                   x0=x0, y0=y0, z0 = z0,                                          &
    1010                   nx = 0, ny = ny, nz = nlev - 1)
    1011 
    1012           CALL init_grid_definition('boundary intermediate', grid=w_west_intermediate,         &
    1013                   xmin = -0.5_dp * dx, xmax = -0.5_dp * dx,                       &
    1014                   ymin =  0.5_dp * dy, ymax = ly - 0.5_dp * dy,                   &
    1015                   x0=x0, y0=y0, z0 = z0,                                          &
    1016                   nx = 0, ny = ny, nz = nlev - 1)
    1017 
    1018           CALL init_grid_definition('boundary intermediate', grid=w_north_intermediate,        &
    1019                   xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx,                    &
    1020                   ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy,               &
    1021                   x0=x0, y0=y0, z0 = z0,                                          &
    1022                   nx = nx, ny = 0, nz = nlev - 1)
    1023 
    1024           CALL init_grid_definition('boundary intermediate', grid=w_south_intermediate,        &
    1025                   xmin =  0.5_dp * dx, xmax = lx - 0.5_dp * dx,                   &
    1026                   ymin = -0.5_dp * dy, ymax = -0.5_dp * dy,                       &
    1027                   x0=x0, y0=y0, z0 = z0,                                          &
    1028                   nx = nx, ny = 0, nz = nlev - 1)
    1029 
    1030           CALL init_grid_definition('boundary intermediate', grid=w_top_intermediate,          &
    1031                   xmin =  0.5_dp * dx, xmax = lx - 0.5_dp * dx,                   &
    1032                   ymin =  0.5_dp * dy, ymax = ly - 0.5_dp * dy,                   &
    1033                   x0=x0, y0=y0, z0 = z0,                                          &
    1034                   nx = nx, ny = ny, nz = nlev - 1)
    1035        ENDIF
     1056    ENDIF
    10361057
    10371058!                                                                             
    … …  
    10401061!------------------------------------------------------------------------------
    10411062
    1042        lonmin_palm = MINVAL(palm_intermediate % clon)
    1043        lonmax_palm = MAXVAL(palm_intermediate % clon)
    1044        latmin_palm = MINVAL(palm_intermediate % clat)
    1045        latmax_palm = MAXVAL(palm_intermediate % clat)
    1046 
    1047        CALL init_averaging_grid(averaged_initial_scalar_profile, cosmo_grid,   &
    1048                x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0,               &
    1049                lonmin = lonmin_palm, lonmax = lonmax_palm,                     &
    1050                latmin = latmin_palm, latmax = latmax_palm,                     &
    1051                kind='scalar', name='averaged initial scalar')
    1052 
    1053        CALL init_averaging_grid(averaged_initial_w_profile, cosmo_grid,        &
    1054                x = 0.5_dp * lx, y = 0.5_dp * ly, z = zw, z0 = z0,              &
    1055                lonmin = lonmin_palm, lonmax = lonmax_palm,                     &
    1056                latmin = latmin_palm, latmax = latmax_palm,                     &
    1057                kind='w', name='averaged initial w')
    1058 
    1059        CALL init_averaging_grid(averaged_scalar_profile, cosmo_grid,           &
    1060                x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0,               &
    1061                lonmin = lam_west, lonmax = lam_east,                           &
    1062                latmin = phi_south, latmax = phi_north,                         &
    1063                kind='scalar', name='centre geostrophic scalar')
    1064 
    1065        CALL init_averaging_grid(averaged_w_profile, cosmo_grid,                &
    1066                x = 0.5_dp * lx, y = 0.5_dp * ly, z = zw, z0 = z0,              &
    1067                lonmin = lam_west, lonmax = lam_east,                           &
    1068                latmin = phi_south, latmax = phi_north,                         &
    1069                kind='w', name='centre geostrophic w')
    1070 
    1071        CALL init_averaging_grid(south_averaged_scalar_profile, cosmo_grid,     &
    1072                x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0,               &
    1073                lonmin = lam_west, lonmax = lam_east,                           &
    1074                latmin = phi_centre - averaging_angle,                          &
    1075                latmax = phi_centre,                                            &
    1076                kind='scalar', name='south geostrophic scalar')
    1077 
    1078        CALL init_averaging_grid(north_averaged_scalar_profile, cosmo_grid,     &
    1079                x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0,               &
    1080                lonmin = lam_west, lonmax = lam_east,                           &
    1081                latmin = phi_centre,                                            &
    1082                latmax = phi_centre + averaging_angle,                          &
    1083                kind='scalar', name='north geostrophic scalar')
    1084 
    1085        CALL init_averaging_grid(west_averaged_scalar_profile, cosmo_grid,      &
    1086                x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0,               &
    1087                lonmin = lam_centre - averaging_angle,                          &
    1088                lonmax = lam_centre,                                            &
    1089                latmin = phi_south, latmax = phi_north,                         &
    1090                kind='scalar', name='west geostrophic scalar')
    1091 
    1092        CALL init_averaging_grid(east_averaged_scalar_profile, cosmo_grid,      &
    1093                x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0,               &
    1094                lonmin = lam_centre,                                            &
    1095                lonmax = lam_centre + averaging_angle,                          &
    1096                latmin = phi_south, latmax = phi_north,                         &
    1097                kind='scalar', name='east geostrophic scalar')
     1063    lonmin_palm = MINVAL(palm_intermediate%clon)
     1064    lonmax_palm = MAXVAL(palm_intermediate%clon)
     1065    latmin_palm = MINVAL(palm_intermediate%clat)
     1066    latmax_palm = MAXVAL(palm_intermediate%clat)
     1067
     1068    CALL init_averaging_grid(averaged_initial_scalar_profile, cosmo_grid,   &
     1069            x = 0.5_wp * lx, y = 0.5_wp * ly, z = z, z0 = z0,               &
     1070            lonmin = lonmin_palm, lonmax = lonmax_palm,                     &
     1071            latmin = latmin_palm, latmax = latmax_palm,                     &
     1072            kind='scalar', name='averaged initial scalar')
     1073
     1074    CALL init_averaging_grid(averaged_initial_w_profile, cosmo_grid,        &
     1075            x = 0.5_wp * lx, y = 0.5_wp * ly, z = zw, z0 = z0,              &
     1076            lonmin = lonmin_palm, lonmax = lonmax_palm,                     &
     1077            latmin = latmin_palm, latmax = latmax_palm,                     &
     1078            kind='w', name='averaged initial w')
     1079
     1080    CALL init_averaging_grid(averaged_scalar_profile, cosmo_grid,           &
     1081            x = 0.5_wp * lx, y = 0.5_wp * ly, z = z, z0 = z0,               &
     1082            lonmin = lam_west, lonmax = lam_east,                           &
     1083            latmin = phi_south, latmax = phi_north,                         &
     1084            kind='scalar', name='centre geostrophic scalar')
     1085
     1086    CALL init_averaging_grid(averaged_w_profile, cosmo_grid,                &
     1087            x = 0.5_wp * lx, y = 0.5_wp * ly, z = zw, z0 = z0,              &
     1088            lonmin = lam_west, lonmax = lam_east,                           &
     1089            latmin = phi_south, latmax = phi_north,                         &
     1090            kind='w', name='centre geostrophic w')
     1091
     1092    CALL init_averaging_grid(south_averaged_scalar_profile, cosmo_grid,     &
     1093            x = 0.5_wp * lx, y = 0.5_wp * ly, z = z, z0 = z0,               &
     1094            lonmin = lam_west, lonmax = lam_east,                           &
     1095            latmin = phi_centre - averaging_angle,                          &
     1096            latmax = phi_centre,                                            &
     1097            kind='scalar', name='south geostrophic scalar')
     1098
     1099    CALL init_averaging_grid(north_averaged_scalar_profile, cosmo_grid,     &
     1100            x = 0.5_wp * lx, y = 0.5_wp * ly, z = z, z0 = z0,               &
     1101            lonmin = lam_west, lonmax = lam_east,                           &
     1102            latmin = phi_centre,                                            &
     1103            latmax = phi_centre + averaging_angle,                          &
     1104            kind='scalar', name='north geostrophic scalar')
     1105
     1106    CALL init_averaging_grid(west_averaged_scalar_profile, cosmo_grid,      &
     1107            x = 0.5_wp * lx, y = 0.5_wp * ly, z = z, z0 = z0,               &
     1108            lonmin = lam_centre - averaging_angle,                          &
     1109            lonmax = lam_centre,                                            &
     1110            latmin = phi_south, latmax = phi_north,                         &
     1111            kind='scalar', name='west geostrophic scalar')
     1112
     1113    CALL init_averaging_grid(east_averaged_scalar_profile, cosmo_grid,      &
     1114            x = 0.5_wp * lx, y = 0.5_wp * ly, z = z, z0 = z0,               &
     1115            lonmin = lam_centre,                                            &
     1116            lonmax = lam_centre + averaging_angle,                          &
     1117            latmin = phi_south, latmax = phi_north,                         &
     1118            kind='scalar', name='east geostrophic scalar')
    10981119
    10991120
    … …  
    11021123! Section 4: Precompute neighbours and weights for interpolation             
    11031124!------------------------------------------------------------------------------
    1104        interp_mode = 's'
    1105        CALL setup_interpolation(cosmo_grid, palm_grid, palm_intermediate, interp_mode, ic_mode=cfg % ic_mode)
    1106        IF (boundary_variables_required)  THEN
    1107           CALL setup_interpolation(cosmo_grid, scalars_east_grid, scalars_east_intermediate, interp_mode)
    1108           CALL setup_interpolation(cosmo_grid, scalars_west_grid, scalars_west_intermediate, interp_mode)
    1109           CALL setup_interpolation(cosmo_grid, scalars_north_grid, scalars_north_intermediate, interp_mode)
    1110           CALL setup_interpolation(cosmo_grid, scalars_south_grid, scalars_south_intermediate, interp_mode)
    1111           CALL setup_interpolation(cosmo_grid, scalars_top_grid, scalars_top_intermediate, interp_mode)
    1112        ENDIF
    1113 
    1114        interp_mode = 'u'
    1115        CALL setup_interpolation(cosmo_grid, u_initial_grid, u_initial_intermediate, interp_mode, ic_mode=cfg % ic_mode)
    1116        IF (boundary_variables_required)  THEN
    1117           CALL setup_interpolation(cosmo_grid, u_east_grid, u_east_intermediate, interp_mode)
    1118           CALL setup_interpolation(cosmo_grid, u_west_grid, u_west_intermediate, interp_mode)
    1119           CALL setup_interpolation(cosmo_grid, u_north_grid, u_north_intermediate, interp_mode)
    1120           CALL setup_interpolation(cosmo_grid, u_south_grid, u_south_intermediate, interp_mode)
    1121           CALL setup_interpolation(cosmo_grid, u_top_grid, u_top_intermediate, interp_mode)
    1122        ENDIF
    1123 
    1124        interp_mode = 'v'
    1125        CALL setup_interpolation(cosmo_grid, v_initial_grid, v_initial_intermediate, interp_mode, ic_mode=cfg % ic_mode)
    1126        IF (boundary_variables_required)  THEN
    1127           CALL setup_interpolation(cosmo_grid, v_east_grid, v_east_intermediate, interp_mode)
    1128           CALL setup_interpolation(cosmo_grid, v_west_grid, v_west_intermediate, interp_mode)
    1129           CALL setup_interpolation(cosmo_grid, v_north_grid, v_north_intermediate, interp_mode)
    1130           CALL setup_interpolation(cosmo_grid, v_south_grid, v_south_intermediate, interp_mode)
    1131           CALL setup_interpolation(cosmo_grid, v_top_grid, v_top_intermediate, interp_mode)
    1132        ENDIF
    1133 
    1134        interp_mode = 'w'
    1135        CALL setup_interpolation(cosmo_grid, w_initial_grid, w_initial_intermediate, interp_mode, ic_mode=cfg % ic_mode)
    1136        IF (boundary_variables_required)  THEN
    1137           CALL setup_interpolation(cosmo_grid, w_east_grid, w_east_intermediate, interp_mode)
    1138           CALL setup_interpolation(cosmo_grid, w_west_grid, w_west_intermediate, interp_mode)
    1139           CALL setup_interpolation(cosmo_grid, w_north_grid, w_north_intermediate, interp_mode)
    1140           CALL setup_interpolation(cosmo_grid, w_south_grid, w_south_intermediate, interp_mode)
    1141           CALL setup_interpolation(cosmo_grid, w_top_grid, w_top_intermediate, interp_mode)
    1142        ENDIF
    1143 
    1144        IF (TRIM(cfg % ic_mode) == 'profile')  THEN
    1145            !TODO: remove this conditional if not needed.
    1146        ENDIF
    1147        
    1148 
    1149     END SUBROUTINE setup_grids
     1125    interp_mode = 's'
     1126    CALL setup_interpolation(cosmo_grid, palm_grid, palm_intermediate, interp_mode, ic_mode=cfg%ic_mode)
     1127    IF (boundary_variables_required)  THEN
     1128       CALL setup_interpolation(cosmo_grid, scalars_east_grid, scalars_east_intermediate, interp_mode)
     1129       CALL setup_interpolation(cosmo_grid, scalars_west_grid, scalars_west_intermediate, interp_mode)
     1130       CALL setup_interpolation(cosmo_grid, scalars_north_grid, scalars_north_intermediate, interp_mode)
     1131       CALL setup_interpolation(cosmo_grid, scalars_south_grid, scalars_south_intermediate, interp_mode)
     1132       CALL setup_interpolation(cosmo_grid, scalars_top_grid, scalars_top_intermediate, interp_mode)
     1133    ENDIF
     1134
     1135    interp_mode = 'u'
     1136    CALL setup_interpolation(cosmo_grid, u_initial_grid, u_initial_intermediate, interp_mode, ic_mode=cfg%ic_mode)
     1137    IF (boundary_variables_required)  THEN
     1138       CALL setup_interpolation(cosmo_grid, u_east_grid, u_east_intermediate, interp_mode)
     1139       CALL setup_interpolation(cosmo_grid, u_west_grid, u_west_intermediate, interp_mode)
     1140       CALL setup_interpolation(cosmo_grid, u_north_grid, u_north_intermediate, interp_mode)
     1141       CALL setup_interpolation(cosmo_grid, u_south_grid, u_south_intermediate, interp_mode)
     1142       CALL setup_interpolation(cosmo_grid, u_top_grid, u_top_intermediate, interp_mode)
     1143    ENDIF
     1144
     1145    interp_mode = 'v'
     1146    CALL setup_interpolation(cosmo_grid, v_initial_grid, v_initial_intermediate, interp_mode, ic_mode=cfg%ic_mode)
     1147    IF (boundary_variables_required)  THEN
     1148       CALL setup_interpolation(cosmo_grid, v_east_grid, v_east_intermediate, interp_mode)
     1149       CALL setup_interpolation(cosmo_grid, v_west_grid, v_west_intermediate, interp_mode)
     1150       CALL setup_interpolation(cosmo_grid, v_north_grid, v_north_intermediate, interp_mode)
     1151       CALL setup_interpolation(cosmo_grid, v_south_grid, v_south_intermediate, interp_mode)
     1152       CALL setup_interpolation(cosmo_grid, v_top_grid, v_top_intermediate, interp_mode)
     1153    ENDIF
     1154
     1155    interp_mode = 'w'
     1156    CALL setup_interpolation(cosmo_grid, w_initial_grid, w_initial_intermediate, interp_mode, ic_mode=cfg%ic_mode)
     1157    IF (boundary_variables_required)  THEN
     1158       CALL setup_interpolation(cosmo_grid, w_east_grid, w_east_intermediate, interp_mode)
     1159       CALL setup_interpolation(cosmo_grid, w_west_grid, w_west_intermediate, interp_mode)
     1160       CALL setup_interpolation(cosmo_grid, w_north_grid, w_north_intermediate, interp_mode)
     1161       CALL setup_interpolation(cosmo_grid, w_south_grid, w_south_intermediate, interp_mode)
     1162       CALL setup_interpolation(cosmo_grid, w_top_grid, w_top_intermediate, interp_mode)
     1163    ENDIF
     1164
     1165    IF (TRIM(cfg%ic_mode) == 'profile')  THEN
     1166        !TODO: remove this conditional if not needed.
     1167    ENDIF
     1168
     1169 END SUBROUTINE setup_grids
    11501170
    11511171
    … …  
    11561176!> vertical interpolation.
    11571177!------------------------------------------------------------------------------!
    1158     SUBROUTINE setup_interpolation(cosmo_grid, grid, intermediate_grid, kind, ic_mode)
    1159 
    1160        TYPE(grid_definition), INTENT(IN), TARGET    ::  cosmo_grid
    1161        TYPE(grid_definition), INTENT(INOUT), TARGET ::  grid, intermediate_grid
    1162        CHARACTER, INTENT(IN)                        ::  kind
    1163        CHARACTER(LEN=*), INTENT(IN), OPTIONAL       ::  ic_mode
    1164 
    1165        REAL(dp), DIMENSION(:), POINTER     ::  cosmo_lat, cosmo_lon
    1166        REAL(dp), DIMENSION(:,:,:), POINTER ::  cosmo_h
    1167 
    1168        LOGICAL :: setup_volumetric
     1178 SUBROUTINE setup_interpolation(cosmo_grid, grid, intermediate_grid, kind, ic_mode)
     1179
     1180    TYPE(grid_definition), INTENT(IN), TARGET    ::  cosmo_grid
     1181    TYPE(grid_definition), INTENT(INOUT), TARGET ::  grid, intermediate_grid
     1182    CHARACTER, INTENT(IN)                        ::  kind
     1183    CHARACTER(LEN=*), INTENT(IN), OPTIONAL       ::  ic_mode
     1184
     1185    REAL(wp), DIMENSION(:), POINTER     ::  cosmo_lat, cosmo_lon
     1186    REAL(wp), DIMENSION(:,:,:), POINTER ::  cosmo_h
     1187
     1188    LOGICAL :: setup_volumetric
    11691189
    11701190!------------------------------------------------------------------------------
    … …  
    11721192!------------------------------------------------------------------------------
    11731193!
    1174 !--    Select horizontal coordinates according to kind of points (s/w, u, v)
    1175        SELECT CASE(kind)
    1176 
    1177 !
    1178 !--    scalars
    1179        CASE('s')
    1180 
    1181           cosmo_lat => cosmo_grid % lat
    1182           cosmo_lon => cosmo_grid % lon
    1183           cosmo_h   => cosmo_grid % hfl
    1184 !
    1185 !--    vertical velocity
    1186        CASE('w')
    1187 
    1188           cosmo_lat => cosmo_grid % lat
    1189           cosmo_lon => cosmo_grid % lon
    1190           cosmo_h   => cosmo_grid % hhl
    1191 !
    1192 !--    x velocity
    1193        CASE('u')
    1194 
    1195           cosmo_lat => cosmo_grid % lat
    1196           cosmo_lon => cosmo_grid % lonu
    1197           cosmo_h   => cosmo_grid % hfl
    1198 
    1199 !
    1200 !--    y velocity
    1201        CASE('v')
    1202 
    1203           cosmo_lat => cosmo_grid % latv
    1204           cosmo_lon => cosmo_grid % lon
    1205           cosmo_h   => cosmo_grid % hfl
    1206 
    1207        CASE DEFAULT
    1208 
    1209           message = "Interpolation quantity '" // kind // "' is not supported."
    1210           CALL inifor_abort('setup_interpolation', message)
    1211 
    1212        END SELECT
    1213 
    1214        CALL find_horizontal_neighbours(cosmo_lat, cosmo_lon,                   &
    1215           intermediate_grid % clat, intermediate_grid % clon,                  &
    1216           intermediate_grid % ii, intermediate_grid % jj)
    1217 
    1218        CALL compute_horizontal_interp_weights(cosmo_lat, cosmo_lon,            &
    1219           intermediate_grid % clat, intermediate_grid % clon,                  &
    1220           intermediate_grid % ii, intermediate_grid % jj,                      &
    1221           intermediate_grid % w_horiz)
     1194!-- Select horizontal coordinates according to kind of points (s/w, u, v)
     1195    SELECT CASE(kind)
     1196
     1197!
     1198!-- scalars
     1199    CASE('s')
     1200
     1201       cosmo_lat => cosmo_grid%lat
     1202       cosmo_lon => cosmo_grid%lon
     1203       cosmo_h   => cosmo_grid%hfl
     1204!
     1205!-- vertical velocity
     1206    CASE('w')
     1207
     1208       cosmo_lat => cosmo_grid%lat
     1209       cosmo_lon => cosmo_grid%lon
     1210       cosmo_h   => cosmo_grid%hhl
     1211!
     1212!-- x velocity
     1213    CASE('u')
     1214
     1215       cosmo_lat => cosmo_grid%lat
     1216       cosmo_lon => cosmo_grid%lonu
     1217       cosmo_h   => cosmo_grid%hfl
     1218
     1219!
     1220!-- y velocity
     1221    CASE('v')
     1222
     1223       cosmo_lat => cosmo_grid%latv
     1224       cosmo_lon => cosmo_grid%lon
     1225       cosmo_h   => cosmo_grid%hfl
     1226
     1227    CASE DEFAULT
     1228
     1229       message = "Interpolation quantity '" // kind // "' is not supported."
     1230       CALL inifor_abort('setup_interpolation', message)
     1231
     1232    END SELECT
     1233
     1234    CALL find_horizontal_neighbours(cosmo_lat, cosmo_lon,                   &
     1235       intermediate_grid%clat, intermediate_grid%clon,                  &
     1236       intermediate_grid%ii, intermediate_grid%jj)
     1237
     1238    CALL compute_horizontal_interp_weights(cosmo_lat, cosmo_lon,            &
     1239       intermediate_grid%clat, intermediate_grid%clon,                  &
     1240       intermediate_grid%ii, intermediate_grid%jj,                      &
     1241       intermediate_grid%w_horiz)
    12221242
    12231243!------------------------------------------------------------------------------
    … …  
    12261246
    12271247!
    1228 !--    If profile initialization is chosen, we--somewhat counterintuitively--
    1229 !--    don't need to compute vertical interpolation weights. At least, we
    1230 !--    don't need them on the intermediate grid, which fills the entire PALM
    1231 !--    domain volume. Instead we need vertical weights for the intermediate
    1232 !--    profile grids, which get computed in setup_averaging().
    1233        setup_volumetric = .TRUE.
    1234        IF (PRESENT(ic_mode))  THEN
    1235           IF (TRIM(ic_mode) == 'profile')  setup_volumetric = .FALSE.
    1236        ENDIF
    1237 
    1238        IF (setup_volumetric)  THEN
    1239           ALLOCATE( intermediate_grid % h(0:intermediate_grid % nx,            &
    1240                                           0:intermediate_grid % ny,            &
    1241                                           0:intermediate_grid % nz) )
    1242           intermediate_grid % h(:,:,:) = - EARTH_RADIUS
    1243 
    1244 !
    1245 !--       For w points, use hhl, for scalars use hfl
    1246 !--       compute the full heights for the intermediate grids
    1247           CALL interpolate_2d(cosmo_h, intermediate_grid % h, intermediate_grid)
    1248           CALL find_vertical_neighbours_and_weights_interp(grid, intermediate_grid)
    1249        ENDIF
    1250        
    1251     END SUBROUTINE setup_interpolation
     1248!-- If profile initialization is chosen, we--somewhat counterintuitively--
     1249!-- don't need to compute vertical interpolation weights. At least, we
     1250!-- don't need them on the intermediate grid, which fills the entire PALM
     1251!-- domain volume. Instead we need vertical weights for the intermediate
     1252!-- profile grids, which get computed in setup_averaging().
     1253    setup_volumetric = .TRUE.
     1254    IF (PRESENT(ic_mode))  THEN
     1255       IF (TRIM(ic_mode) == 'profile')  setup_volumetric = .FALSE.
     1256    ENDIF
     1257
     1258    IF (setup_volumetric)  THEN
     1259       ALLOCATE( intermediate_grid%h(0:intermediate_grid%nx,            &
     1260                                     0:intermediate_grid%ny,            &
     1261                                     0:intermediate_grid%nz) )
     1262       intermediate_grid%h(:,:,:) = - EARTH_RADIUS
     1263
     1264!
     1265!--    For w points, use hhl, for scalars use hfl
     1266!--    compute the full heights for the intermediate grids
     1267       CALL interpolate_2d(cosmo_h, intermediate_grid%h, intermediate_grid)
     1268       CALL find_vertical_neighbours_and_weights_interp(grid, intermediate_grid)
     1269    ENDIF
     1270       
     1271 END SUBROUTINE setup_interpolation
    12521272
    12531273!------------------------------------------------------------------------------!
    … …  
    12811301!> grid : Grid variable to be initialized.
    12821302!------------------------------------------------------------------------------!
    1283     SUBROUTINE init_grid_definition(kind, xmin, xmax, ymin, ymax,              &
    1284                                     x0, y0, z0, nx, ny, nz, z, zw, grid, ic_mode)
    1285         CHARACTER(LEN=*), INTENT(IN)           ::  kind
    1286         CHARACTER(LEN=*), INTENT(IN), OPTIONAL ::  ic_mode
    1287         INTEGER, INTENT(IN)                    ::  nx, ny, nz
    1288         REAL(dp), INTENT(IN)                   ::  xmin, xmax, ymin, ymax
    1289         REAL(dp), INTENT(IN)                   ::  x0, y0, z0
    1290         REAL(dp), INTENT(IN), TARGET, OPTIONAL ::  z(:)
    1291         REAL(dp), INTENT(IN), TARGET, OPTIONAL ::  zw(:)
    1292         TYPE(grid_definition), INTENT(INOUT)   ::  grid
    1293 
    1294         grid % nx = nx
    1295         grid % ny = ny
    1296         grid % nz = nz
    1297 
    1298         grid % lx = xmax - xmin
    1299         grid % ly = ymax - ymin
    1300 
    1301         grid % x0 = x0
    1302         grid % y0 = y0
    1303         grid % z0 = z0
    1304 
    1305         SELECT CASE( TRIM(kind) )
    1306 
    1307         CASE('boundary')
    1308 
    1309            IF (.NOT.PRESENT(z))  THEN
    1310               message = "z has not been passed but is required for 'boundary' grids"
    1311               CALL inifor_abort('init_grid_definition', message)
    1312            ENDIF
    1313 
    1314            ALLOCATE( grid % x(0:nx) )
    1315            CALL linspace(xmin, xmax, grid % x)
    1316 
    1317            ALLOCATE( grid % y(0:ny) )
    1318            CALL linspace(ymin, ymax, grid % y)
    1319 
    1320            grid % z => z
    1321 
    1322 !
    1323 !--        Allocate neighbour indices and weights
    1324            IF (TRIM(ic_mode) .NE. 'profile')  THEN
    1325               ALLOCATE( grid % kk(0:nx, 0:ny, 1:nz, 2) )
    1326               grid % kk(:,:,:,:) = -1
    1327 
    1328               ALLOCATE( grid % w_verti(0:nx, 0:ny, 1:nz, 2) )
    1329               grid % w_verti(:,:,:,:) = 0.0_dp
    1330            ENDIF
    1331         
    1332         CASE('boundary intermediate')
    1333 
    1334            ALLOCATE( grid % x(0:nx) )
    1335            CALL linspace(xmin, xmax, grid % x)
    1336 
    1337            ALLOCATE( grid % y(0:ny) )
    1338            CALL linspace(ymin, ymax, grid % y)
    1339 
    1340            ALLOCATE( grid % clon(0:nx, 0:ny), grid % clat(0:nx, 0:ny)  )
    1341 
    1342            CALL rotate_to_cosmo(                                               &
    1343               phir = project( grid % y, y0, EARTH_RADIUS ) ,                   & ! = plate-carree latitude
    1344               lamr = project( grid % x, x0, EARTH_RADIUS ) ,                   & ! = plate-carree longitude
    1345               phip = phi_cn, lamp = lambda_cn,                                 &
    1346               phi  = grid % clat,                                              &
    1347               lam  = grid % clon,                                              &
    1348               gam  = gam                                                       &
    1349            )
    1350 
    1351 !
    1352 !--        Allocate neighbour indices and weights
    1353            ALLOCATE( grid % ii(0:nx, 0:ny, 4),                                 &
    1354                      grid % jj(0:nx, 0:ny, 4) )
    1355            grid % ii(:,:,:)   = -1
    1356            grid % jj(:,:,:)   = -1
    1357 
    1358            ALLOCATE( grid % w_horiz(0:nx, 0:ny, 4) )
    1359            grid % w_horiz(:,:,:)   = 0.0_dp
    1360         
    1361 !
    1362 !--     This mode initializes a Cartesian PALM-4U grid and adds the
    1363 !--     corresponding latitudes and longitudes of the rotated pole grid.
    1364         CASE('palm')
    1365 
    1366            IF (.NOT.PRESENT(z))  THEN
    1367               message = "z has not been passed but is required for 'palm' grids"
    1368               CALL inifor_abort('init_grid_definition', message)
    1369            ENDIF
    1370 
    1371            IF (.NOT.PRESENT(zw))  THEN
    1372               message = "zw has not been passed but is required for 'palm' grids"
    1373               CALL inifor_abort('init_grid_definition', message)
    1374            ENDIF
    1375 
    1376            grid % name(1) = 'x and lon'
    1377            grid % name(2) = 'y and lat'
    1378            grid % name(3) = 'z'
    1379 
    1380 !
    1381 !--        TODO: Remove use of global dx, dy, dz variables. Consider
    1382 !--        TODO: associating global x,y, and z arrays.
    1383            ALLOCATE( grid % x(0:nx),   grid % y(0:ny) )
    1384            ALLOCATE( grid % xu(1:nx),  grid % yv(1:ny) )
    1385            CALL linspace(xmin + 0.5_dp* dx, xmax - 0.5_dp* dx, grid % x)
    1386            CALL linspace(ymin + 0.5_dp* dy, ymax - 0.5_dp* dy, grid % y)
    1387            grid % z => z
    1388            CALL linspace(xmin +  dx, xmax -  dx, grid % xu)
    1389            CALL linspace(ymin +  dy, ymax -  dy, grid % yv)
    1390            grid % zw => zw
    1391 
    1392            grid % depths => depths
    1393 
    1394 !
    1395 !--        Allocate neighbour indices and weights
    1396            IF (TRIM(ic_mode) .NE. 'profile')  THEN
    1397               ALLOCATE( grid % kk(0:nx, 0:ny, 1:nz, 2) )
    1398               grid % kk(:,:,:,:) = -1
    1399 
    1400               ALLOCATE( grid % w_verti(0:nx, 0:ny, 1:nz, 2) )
    1401               grid % w_verti(:,:,:,:) = 0.0_dp
    1402            ENDIF
    1403 
    1404         CASE('palm intermediate')
    1405 
    1406            grid % name(1) = 'x and lon'
    1407            grid % name(2) = 'y and lat'
    1408            grid % name(3) = 'interpolated hhl or hfl'
    1409 
    1410 !
    1411 !--        TODO: Remove use of global dx, dy, dz variables. Consider
    1412 !--        TODO: associating global x,y, and z arrays.
    1413            ALLOCATE( grid % x(0:nx),   grid % y(0:ny) )
    1414            ALLOCATE( grid % xu(1:nx),  grid % yv(1:ny) )
    1415            CALL linspace(xmin + 0.5_dp*dx, xmax - 0.5_dp*dx, grid % x)
    1416            CALL linspace(ymin + 0.5_dp*dy, ymax - 0.5_dp*dy, grid % y)
    1417            CALL linspace(xmin + dx, xmax - dx, grid % xu)
    1418            CALL linspace(ymin + dy, ymax - dy, grid % yv)
    1419 
    1420            grid % depths => depths
    1421 
    1422 !
    1423 !--        Allocate rotated-pole coordinates, clon is for (c)osmo-de (lon)gitude
    1424            ALLOCATE( grid % clon(0:nx, 0:ny),   grid % clat(0:nx, 0:ny)  )
    1425            ALLOCATE( grid % clonu(1:nx, 0:ny),  grid % clatu(1:nx, 0:ny) )
    1426            ALLOCATE( grid % clonv(0:nx, 1:ny),  grid % clatv(0:nx, 1:ny) )
    1427 
    1428 !
    1429 !--        Compute rotated-pole coordinates of...
    1430 !--        ... PALM-4U centres
    1431            CALL rotate_to_cosmo(                                               &
    1432               phir = project( grid % y, y0, EARTH_RADIUS ) , & ! = plate-carree latitude
    1433               lamr = project( grid % x, x0, EARTH_RADIUS ) , & ! = plate-carree longitude
    1434               phip = phi_cn, lamp = lambda_cn,                                 &
    1435               phi  = grid % clat,                                              &
    1436               lam  = grid % clon,                                              &
    1437               gam  = gam                                                       &
    1438            )
    1439 
    1440 !
    1441 !--        ... PALM-4U u winds
    1442            CALL rotate_to_cosmo(                                               &
    1443               phir = project( grid % y,  y0, EARTH_RADIUS ), & ! = plate-carree latitude
    1444               lamr = project( grid % xu, x0, EARTH_RADIUS ), & ! = plate-carree longitude
    1445               phip = phi_cn, lamp = lambda_cn,                                 &
    1446               phi  = grid % clatu,                                             &
    1447               lam  = grid % clonu,                                             &
    1448               gam  = gam                                                       &
    1449            )
    1450 
    1451 !
    1452 !--        ... PALM-4U v winds
    1453            CALL rotate_to_cosmo(                                               &
    1454               phir = project( grid % yv, y0, EARTH_RADIUS ), & ! = plate-carree latitude
    1455               lamr = project( grid % x,  x0, EARTH_RADIUS ), & ! = plate-carree longitude
    1456               phip = phi_cn, lamp = lambda_cn,                                 &
    1457               phi  = grid % clatv,                                             &
    1458               lam  = grid % clonv,                                             &
    1459               gam  = gam                                                       &
    1460            )
    1461 
    1462 !
    1463 !--        Allocate neighbour indices and weights
    1464            ALLOCATE( grid % ii(0:nx, 0:ny, 4),                                 &
    1465                      grid % jj(0:nx, 0:ny, 4) )
    1466            grid % ii(:,:,:)   = -1
    1467            grid % jj(:,:,:)   = -1
    1468 
    1469            ALLOCATE( grid % w_horiz(0:nx, 0:ny, 4) )
    1470            grid % w_horiz(:,:,:)   = 0.0_dp
    1471 
    1472         CASE('cosmo-de')
    1473            grid % name(1) = 'rlon'         ! of COMSO-DE cell centres (scalars)
    1474            grid % name(2) = 'rlat'         ! of COMSO-DE cell centres (scalars)
    1475            grid % name(3) = 'height'
    1476 
    1477            ALLOCATE( grid % lon(0:nx),   grid % lat(0:ny)  )
    1478            ALLOCATE( grid % lonu(0:nx),  grid % latv(0:ny) )
    1479 
    1480            CALL linspace(xmin, xmax, grid % lon)
    1481            CALL linspace(ymin, ymax, grid % lat)
    1482            grid % lonu(:) = grid % lon + 0.5_dp * (grid % lx / grid % nx)
    1483            grid % latv(:) = grid % lat + 0.5_dp * (grid % ly / grid % ny)
    1484 
    1485 !
    1486 !--        Point to heights of half levels (hhl) and compute heights of full
    1487 !--        levels (hfl) as arithmetic averages
    1488            grid % hhl => hhl
    1489            grid % hfl => hfl
    1490            grid % depths => depths
    1491 
    1492         CASE DEFAULT
    1493             message = "Grid kind '" // TRIM(kind) // "' is not recognized."
    1494             CALL inifor_abort('init_grid_definition', message)
    1495 
    1496         END SELECT
    1497 
    1498     END SUBROUTINE init_grid_definition
     1303 SUBROUTINE init_grid_definition(kind, xmin, xmax, ymin, ymax,              &
     1304                                 x0, y0, z0, nx, ny, nz, z, zw, grid, ic_mode)
     1305    CHARACTER(LEN=*), INTENT(IN)           ::  kind
     1306    CHARACTER(LEN=*), INTENT(IN), OPTIONAL ::  ic_mode
     1307    INTEGER, INTENT(IN)                    ::  nx, ny, nz
     1308    REAL(wp), INTENT(IN)                   ::  xmin, xmax, ymin, ymax
     1309    REAL(wp), INTENT(IN)                   ::  x0, y0, z0
     1310    REAL(wp), INTENT(IN), TARGET, OPTIONAL ::  z(:)
     1311    REAL(wp), INTENT(IN), TARGET, OPTIONAL ::  zw(:)
     1312    TYPE(grid_definition), INTENT(INOUT)   ::  grid
     1313
     1314    grid%nx = nx
     1315    grid%ny = ny
     1316    grid%nz = nz
     1317
     1318    grid%lx = xmax - xmin
     1319    grid%ly = ymax - ymin
     1320
     1321    grid%x0 = x0
     1322    grid%y0 = y0
     1323    grid%z0 = z0
     1324
     1325    SELECT CASE( TRIM(kind) )
     1326
     1327       CASE('boundary')
     1328       
     1329          IF (.NOT.PRESENT(z))  THEN
     1330             message = "z has not been passed but is required for 'boundary' grids"
     1331             CALL inifor_abort('init_grid_definition', message)
     1332          ENDIF
     1333       
     1334          ALLOCATE( grid%x(0:nx) )
     1335          CALL linspace(xmin, xmax, grid%x)
     1336       
     1337          ALLOCATE( grid%y(0:ny) )
     1338          CALL linspace(ymin, ymax, grid%y)
     1339       
     1340          grid%z => z
     1341       
     1342!     
     1343!--       Allocate neighbour indices and weights
     1344          IF (TRIM(ic_mode) .NE. 'profile')  THEN
     1345             ALLOCATE( grid%kk(0:nx, 0:ny, 1:nz, 2) )
     1346             grid%kk(:,:,:,:) = -1
     1347       
     1348             ALLOCATE( grid%w_verti(0:nx, 0:ny, 1:nz, 2) )
     1349             grid%w_verti(:,:,:,:) = 0.0_wp
     1350          ENDIF
     1351       
     1352       CASE('boundary intermediate')
     1353       
     1354          ALLOCATE( grid%x(0:nx) )
     1355          CALL linspace(xmin, xmax, grid%x)
     1356       
     1357          ALLOCATE( grid%y(0:ny) )
     1358          CALL linspace(ymin, ymax, grid%y)
     1359       
     1360          ALLOCATE( grid%clon(0:nx, 0:ny), grid%clat(0:nx, 0:ny)  )
     1361       
     1362          CALL rotate_to_cosmo(                                               &
     1363             phir = project( grid%y, y0, EARTH_RADIUS ) ,                   & ! = plate-carree latitude
     1364             lamr = project( grid%x, x0, EARTH_RADIUS ) ,                   & ! = plate-carree longitude
     1365             phip = phi_cn, lamp = lambda_cn,                                 &
     1366             phi  = grid%clat,                                              &
     1367             lam  = grid%clon,                                              &
     1368             gam  = gam                                                       &
     1369          )
     1370       
     1371!     
     1372!--       Allocate neighbour indices and weights
     1373          ALLOCATE( grid%ii(0:nx, 0:ny, 4),                                 &
     1374                    grid%jj(0:nx, 0:ny, 4) )
     1375          grid%ii(:,:,:)   = -1
     1376          grid%jj(:,:,:)   = -1
     1377       
     1378          ALLOCATE( grid%w_horiz(0:nx, 0:ny, 4) )
     1379          grid%w_horiz(:,:,:)   = 0.0_wp
     1380       
     1381!     
     1382!--    This mode initializes a Cartesian PALM-4U grid and adds the
     1383!--    corresponding latitudes and longitudes of the rotated pole grid.
     1384       CASE('palm')
     1385       
     1386          IF (.NOT.PRESENT(z))  THEN
     1387             message = "z has not been passed but is required for 'palm' grids"
     1388             CALL inifor_abort('init_grid_definition', message)
     1389          ENDIF
     1390       
     1391          IF (.NOT.PRESENT(zw))  THEN
     1392             message = "zw has not been passed but is required for 'palm' grids"
     1393             CALL inifor_abort('init_grid_definition', message)
     1394          ENDIF
     1395       
     1396          grid%name(1) = 'x and lon'
     1397          grid%name(2) = 'y and lat'
     1398          grid%name(3) = 'z'
     1399       
     1400!     
     1401!--       TODO: Remove use of global dx, dy, dz variables. Consider
     1402!--       TODO: associating global x,y, and z arrays.
     1403          ALLOCATE( grid%x(0:nx),   grid%y(0:ny) )
     1404          ALLOCATE( grid%xu(1:nx),  grid%yv(1:ny) )
     1405          CALL linspace(xmin + 0.5_wp* dx, xmax - 0.5_wp* dx, grid%x)
     1406          CALL linspace(ymin + 0.5_wp* dy, ymax - 0.5_wp* dy, grid%y)
     1407          grid%z => z
     1408          CALL linspace(xmin +  dx, xmax -  dx, grid%xu)
     1409          CALL linspace(ymin +  dy, ymax -  dy, grid%yv)
     1410          grid%zw => zw
     1411       
     1412          grid%depths => depths
     1413       
     1414!     
     1415!--       Allocate neighbour indices and weights
     1416          IF (TRIM(ic_mode) .NE. 'profile')  THEN
     1417             ALLOCATE( grid%kk(0:nx, 0:ny, 1:nz, 2) )
     1418             grid%kk(:,:,:,:) = -1
     1419       
     1420             ALLOCATE( grid%w_verti(0:nx, 0:ny, 1:nz, 2) )
     1421             grid%w_verti(:,:,:,:) = 0.0_wp
     1422          ENDIF
     1423       
     1424       CASE('palm intermediate')
     1425       
     1426          grid%name(1) = 'x and lon'
     1427          grid%name(2) = 'y and lat'
     1428          grid%name(3) = 'interpolated hhl or hfl'
     1429       
     1430!     
     1431!--       TODO: Remove use of global dx, dy, dz variables. Consider
     1432!--       TODO: associating global x,y, and z arrays.
     1433          ALLOCATE( grid%x(0:nx),   grid%y(0:ny) )
     1434          ALLOCATE( grid%xu(1:nx),  grid%yv(1:ny) )
     1435          CALL linspace(xmin + 0.5_wp*dx, xmax - 0.5_wp*dx, grid%x)
     1436          CALL linspace(ymin + 0.5_wp*dy, ymax - 0.5_wp*dy, grid%y)
     1437          CALL linspace(xmin + dx, xmax - dx, grid%xu)
     1438          CALL linspace(ymin + dy, ymax - dy, grid%yv)
     1439       
     1440          grid%depths => depths
     1441       
     1442!     
     1443!--       Allocate rotated-pole coordinates, clon is for (c)osmo-de (lon)gitude
     1444          ALLOCATE( grid%clon(0:nx, 0:ny),   grid%clat(0:nx, 0:ny)  )
     1445          ALLOCATE( grid%clonu(1:nx, 0:ny),  grid%clatu(1:nx, 0:ny) )
     1446          ALLOCATE( grid%clonv(0:nx, 1:ny),  grid%clatv(0:nx, 1:ny) )
     1447       
     1448!     
     1449!--       Compute rotated-pole coordinates of...
     1450!--       ... PALM-4U centres
     1451          CALL rotate_to_cosmo(                                               &
     1452             phir = project( grid%y, y0, EARTH_RADIUS ) , & ! = plate-carree latitude
     1453             lamr = project( grid%x, x0, EARTH_RADIUS ) , & ! = plate-carree longitude
     1454             phip = phi_cn, lamp = lambda_cn,                                 &
     1455             phi  = grid%clat,                                              &
     1456             lam  = grid%clon,                                              &
     1457             gam  = gam                                                       &
     1458          )
     1459       
     1460!     
     1461!--       ... PALM-4U u winds
     1462          CALL rotate_to_cosmo(                                               &
     1463             phir = project( grid%y,  y0, EARTH_RADIUS ), & ! = plate-carree latitude
     1464             lamr = project( grid%xu, x0, EARTH_RADIUS ), & ! = plate-carree longitude
     1465             phip = phi_cn, lamp = lambda_cn,                                 &
     1466             phi  = grid%clatu,                                             &
     1467             lam  = grid%clonu,                                             &
     1468             gam  = gam                                                       &
     1469          )
     1470       
     1471!     
     1472!--       ... PALM-4U v winds
     1473          CALL rotate_to_cosmo(                                               &
     1474             phir = project( grid%yv, y0, EARTH_RADIUS ), & ! = plate-carree latitude
     1475             lamr = project( grid%x,  x0, EARTH_RADIUS ), & ! = plate-carree longitude
     1476             phip = phi_cn, lamp = lambda_cn,                                 &
     1477             phi  = grid%clatv,                                             &
     1478             lam  = grid%clonv,                                             &
     1479             gam  = gam                                                       &
     1480          )
     1481       
     1482!     
     1483!--       Allocate neighbour indices and weights
     1484          ALLOCATE( grid%ii(0:nx, 0:ny, 4),                                 &
     1485                    grid%jj(0:nx, 0:ny, 4) )
     1486          grid%ii(:,:,:)   = -1
     1487          grid%jj(:,:,:)   = -1
     1488       
     1489          ALLOCATE( grid%w_horiz(0:nx, 0:ny, 4) )
     1490          grid%w_horiz(:,:,:)   = 0.0_wp
     1491       
     1492       CASE('cosmo-de')
     1493          grid%name(1) = 'rlon'         ! of COMSO-DE cell centres (scalars)
     1494          grid%name(2) = 'rlat'         ! of COMSO-DE cell centres (scalars)
     1495          grid%name(3) = 'height'
     1496       
     1497          ALLOCATE( grid%lon(0:nx),   grid%lat(0:ny)  )
     1498          ALLOCATE( grid%lonu(0:nx),  grid%latv(0:ny) )
     1499       
     1500          CALL linspace(xmin, xmax, grid%lon)
     1501          CALL linspace(ymin, ymax, grid%lat)
     1502          grid%lonu(:) = grid%lon + 0.5_wp * (grid%lx / grid%nx)
     1503          grid%latv(:) = grid%lat + 0.5_wp * (grid%ly / grid%ny)
     1504       
     1505!     
     1506!--       Point to heights of half levels (hhl) and compute heights of full
     1507!--       levels (hfl) as arithmetic averages
     1508          grid%hhl => hhl
     1509          grid%hfl => hfl
     1510          grid%depths => depths
     1511       
     1512       CASE DEFAULT
     1513           message = "Grid kind '" // TRIM(kind) // "' is not recognized."
     1514           CALL inifor_abort('init_grid_definition', message)
     1515
     1516    END SELECT
     1517
     1518 END SUBROUTINE init_grid_definition
    14991519
    15001520
    … …  
    15271547!> avg_grid : averagin grid to be initialized
    15281548!------------------------------------------------------------------------------!
    1529     SUBROUTINE init_averaging_grid(avg_grid, cosmo_grid, x, y, z, z0,          &
    1530        lonmin, lonmax, latmin, latmax, kind, name)
    1531 
    1532        TYPE(grid_definition), INTENT(INOUT) ::  avg_grid
    1533        TYPE(grid_definition), INTENT(IN)    ::  cosmo_grid
    1534        REAL(dp), INTENT(IN)                 ::  x, y, z0
    1535        REAL(dp), INTENT(IN), TARGET         ::  z(:)
    1536        REAL(dp), INTENT(IN)                 ::  lonmin !< lower longitude bound of the averaging grid region [COSMO rotated-pole rad]
    1537        REAL(dp), INTENT(IN)                 ::  lonmax !< upper longitude bound of the averaging grid region [COSMO rotated-pole rad]
    1538        REAL(dp), INTENT(IN)                 ::  latmin !< lower latitude bound of the averaging grid region [COSMO rotated-pole rad]
    1539        REAL(dp), INTENT(IN)                 ::  latmax !< lower latitude bound of the averaging grid region [COSMO rotated-pole rad]
    1540 
    1541        CHARACTER(LEN=*), INTENT(IN)         ::  kind
    1542        CHARACTER(LEN=*), INTENT(IN)         ::  name
    1543 
    1544        LOGICAL                              ::  level_based_averaging
    1545 
    1546        ALLOCATE( avg_grid % x(1) )
    1547        ALLOCATE( avg_grid % y(1) )
    1548        avg_grid % x(1) = x
    1549        avg_grid % y(1) = y
    1550        avg_grid % z => z
    1551        avg_grid % z0 = z0
    1552 
    1553        avg_grid % nz = SIZE(z, 1)
    1554 
    1555        ALLOCATE( avg_grid % lon(2) )
    1556        ALLOCATE( avg_grid % lat(2) )
    1557        avg_grid % lon(1:2) = (/lonmin, lonmax/)
    1558        avg_grid % lat(1:2) = (/latmin, latmax/)
    1559 
    1560        avg_grid % kind = TRIM(kind)
    1561        avg_grid % name(1) = TRIM(name)
    1562 
    1563 !
    1564 !--    Find and store COSMO columns that fall into the coordinate range
    1565 !--    given by avg_grid % clon, % clat
    1566        CALL get_cosmo_averaging_region(avg_grid, cosmo_grid)
    1567 
    1568        ALLOCATE (avg_grid % kkk(avg_grid % n_columns, avg_grid % nz, 2) )
    1569        ALLOCATE (avg_grid % w(avg_grid % n_columns, avg_grid % nz, 2) )
    1570 !
    1571 !--    Compute average COSMO levels in the averaging region
    1572        SELECT CASE(avg_grid % kind)
     1549 SUBROUTINE init_averaging_grid(avg_grid, cosmo_grid, x, y, z, z0,          &
     1550    lonmin, lonmax, latmin, latmax, kind, name)
     1551
     1552    TYPE(grid_definition), INTENT(INOUT) ::  avg_grid
     1553    TYPE(grid_definition), INTENT(IN)    ::  cosmo_grid
     1554    REAL(wp), INTENT(IN)                 ::  x, y, z0
     1555    REAL(wp), INTENT(IN), TARGET         ::  z(:)
     1556    REAL(wp), INTENT(IN)                 ::  lonmin !< lower longitude bound of the averaging grid region [COSMO rotated-pole rad]
     1557    REAL(wp), INTENT(IN)                 ::  lonmax !< upper longitude bound of the averaging grid region [COSMO rotated-pole rad]
     1558    REAL(wp), INTENT(IN)                 ::  latmin !< lower latitude bound of the averaging grid region [COSMO rotated-pole rad]
     1559    REAL(wp), INTENT(IN)                 ::  latmax !< lower latitude bound of the averaging grid region [COSMO rotated-pole rad]
     1560
     1561    CHARACTER(LEN=*), INTENT(IN)         ::  kind
     1562    CHARACTER(LEN=*), INTENT(IN)         ::  name
     1563
     1564    LOGICAL                              ::  level_based_averaging
     1565
     1566    ALLOCATE( avg_grid%x(1) )
     1567    ALLOCATE( avg_grid%y(1) )
     1568    avg_grid%x(1) = x
     1569    avg_grid%y(1) = y
     1570    avg_grid%z => z
     1571    avg_grid%z0 = z0
     1572
     1573    avg_grid%nz = SIZE(z, 1)
     1574
     1575    ALLOCATE( avg_grid%lon(2) )
     1576    ALLOCATE( avg_grid%lat(2) )
     1577    avg_grid%lon(1:2) = (/lonmin, lonmax/)
     1578    avg_grid%lat(1:2) = (/latmin, latmax/)
     1579
     1580    avg_grid%kind = TRIM(kind)
     1581    avg_grid%name(1) = TRIM(name)
     1582
     1583!
     1584!-- Find and store COSMO columns that fall into the coordinate range
     1585!-- given by avg_grid%clon, %clat
     1586    CALL get_cosmo_averaging_region(avg_grid, cosmo_grid)
     1587
     1588    ALLOCATE (avg_grid%kkk(avg_grid%n_columns, avg_grid%nz, 2) )
     1589    ALLOCATE (avg_grid%w(avg_grid%n_columns, avg_grid%nz, 2) )
     1590!
     1591!-- Compute average COSMO levels in the averaging region
     1592    SELECT CASE(avg_grid%kind)
    15731593
    15741594       CASE('scalar', 'u', 'v')
    1575           avg_grid % cosmo_h => cosmo_grid % hfl
     1595          avg_grid%cosmo_h => cosmo_grid%hfl
    15761596
    15771597       CASE('w')
    1578           avg_grid % cosmo_h => cosmo_grid % hhl
     1598          avg_grid%cosmo_h => cosmo_grid%hhl
    15791599
    15801600       CASE DEFAULT
    1581           message = "Averaging grid kind '" // TRIM(avg_grid % kind) // &
     1601          message = "Averaging grid kind '" // TRIM(avg_grid%kind) // &
    15821602                    "' is not supported. Use 'scalar', 'u', or 'v'."
    15831603          CALL inifor_abort('get_cosmo_averaging_region', message)
    15841604
    1585        END SELECT
    1586 
    1587 !
    1588 !--    For level-besed averaging, compute average heights
    1589        !level_based_averaging = ( TRIM(mode) == 'level' )
    1590        level_based_averaging = ( TRIM(cfg % averaging_mode) == 'level' )
    1591        IF (level_based_averaging)  THEN
    1592           ALLOCATE(avg_grid % h(1,1,SIZE(avg_grid % cosmo_h, 3)) )
     1605    END SELECT
     1606
     1607!
     1608!-- For level-besed averaging, compute average heights
     1609    level_based_averaging = ( TRIM(cfg%averaging_mode) == 'level' )
     1610    IF (level_based_averaging)  THEN
     1611       ALLOCATE(avg_grid%h(1,1,SIZE(avg_grid%cosmo_h, 3)) )
    15931612 
    1594           CALL average_2d(avg_grid % cosmo_h, avg_grid % h(1,1,:),             &
    1595                           avg_grid % iii, avg_grid % jjj)
    1596 
    1597        ENDIF
    1598 
    1599 !
    1600 !--    Compute vertical weights and neighbours
    1601        CALL find_vertical_neighbours_and_weights_average(                      &
    1602           avg_grid, level_based_averaging                                      &
    1603        )
    1604 
    1605     END SUBROUTINE init_averaging_grid
    1606 
    1607 
    1608     SUBROUTINE get_cosmo_averaging_region(avg_grid, cosmo_grid)
    1609        TYPE(grid_definition), INTENT(INOUT)         ::  avg_grid
    1610        TYPE(grid_definition), TARGET, INTENT(IN)    ::  cosmo_grid
    1611 
    1612        REAL(dp), DIMENSION(:), POINTER              ::  cosmo_lon, cosmo_lat
    1613        REAL(dp)                                     ::  dlon, dlat
    1614 
    1615        INTEGER ::  i, j, imin, imax, jmin, jmax, l, nx, ny
    1616 
    1617 
    1618        SELECT CASE( TRIM(avg_grid % kind) )
     1613       CALL average_2d(avg_grid%cosmo_h, avg_grid%h(1,1,:),             &
     1614                       avg_grid%iii, avg_grid%jjj)
     1615
     1616    ENDIF
     1617
     1618!
     1619!-- Compute vertical weights and neighbours
     1620    CALL find_vertical_neighbours_and_weights_average(                      &
     1621       avg_grid, level_based_averaging                                      &
     1622    )
     1623
     1624 END SUBROUTINE init_averaging_grid
     1625
     1626
     1627 SUBROUTINE get_cosmo_averaging_region(avg_grid, cosmo_grid)
     1628    TYPE(grid_definition), INTENT(INOUT)         ::  avg_grid
     1629    TYPE(grid_definition), TARGET, INTENT(IN)    ::  cosmo_grid
     1630
     1631    REAL(wp), DIMENSION(:), POINTER              ::  cosmo_lon, cosmo_lat
     1632    REAL(wp)                                     ::  dlon, dlat
     1633
     1634    INTEGER ::  i, j, imin, imax, jmin, jmax, l, nx, ny
     1635
     1636
     1637    SELECT CASE( TRIM(avg_grid%kind) )
    16191638
    16201639       CASE('scalar', 'w')
    1621           cosmo_lon => cosmo_grid % lon
    1622           cosmo_lat => cosmo_grid % lat
     1640          cosmo_lon => cosmo_grid%lon
     1641          cosmo_lat => cosmo_grid%lat
    16231642
    16241643       CASE('u')
    1625           cosmo_lon => cosmo_grid % lonu
    1626           cosmo_lat => cosmo_grid % lat
     1644          cosmo_lon => cosmo_grid%lonu
     1645          cosmo_lat => cosmo_grid%lat
    16271646
    16281647       CASE('v')
    1629           cosmo_lon => cosmo_grid % lon
    1630           cosmo_lat => cosmo_grid % latv
     1648          cosmo_lon => cosmo_grid%lon
     1649          cosmo_lat => cosmo_grid%latv
    16311650
    16321651       CASE DEFAULT
    1633           message = "Averaging grid kind '" // TRIM(avg_grid % kind) // &
     1652          message = "Averaging grid kind '" // TRIM(avg_grid%kind) // &
    16341653                    "' is not supported. Use 'scalar', 'u', or 'v'."
    16351654          CALL inifor_abort('get_cosmo_averaging_region', message)
    16361655
    1637        END SELECT
    1638 
    1639 !
    1640 !--    FIXME: make dlon, dlat parameters of the grid_defintion type
    1641        dlon = cosmo_lon(1) - cosmo_lon(0)
    1642        dlat = cosmo_lat(1) - cosmo_lat(0)
    1643 
    1644        imin = FLOOR  ( (avg_grid % lon(1) - cosmo_lon(0)) / dlon )
    1645        imax = CEILING( (avg_grid % lon(2) - cosmo_lon(0)) / dlon )
    1646 
    1647        jmin = FLOOR  ( (avg_grid % lat(1) - cosmo_lat(0)) / dlat )
    1648        jmax = CEILING( (avg_grid % lat(2) - cosmo_lat(0)) / dlat )
    1649        
    1650        message = "Grid " // TRIM(avg_grid % name(1)) // " averages over " // &
    1651                  TRIM(str(imin)) // " <= i <= " // TRIM(str(imax)) //          &
    1652                  " and " //                                                    &
    1653                  TRIM(str(jmin)) // " <= j <= " // TRIM(str(jmax))
    1654        CALL report( 'get_cosmo_averaging_region', message )
    1655 
    1656        nx = imax - imin + 1
    1657        ny = jmax - jmin + 1
    1658        avg_grid % n_columns = nx * ny
    1659 
    1660        ALLOCATE( avg_grid % iii(avg_grid % n_columns),                         &
    1661                  avg_grid % jjj(avg_grid % n_columns) )
    1662 
    1663        l = 0
    1664        DO j = jmin, jmax
    1665        DO i = imin, imax
    1666           l = l + 1
    1667           avg_grid % iii(l) = i
    1668           avg_grid % jjj(l) = j
    1669        ENDDO
    1670        ENDDO
    1671 
    1672     END SUBROUTINE get_cosmo_averaging_region
     1656    END SELECT
     1657
     1658!
     1659!-- FIXME: make dlon, dlat parameters of the grid_defintion type
     1660    dlon = cosmo_lon(1) - cosmo_lon(0)
     1661    dlat = cosmo_lat(1) - cosmo_lat(0)
     1662
     1663    imin = FLOOR  ( (avg_grid%lon(1) - cosmo_lon(0)) / dlon )
     1664    imax = CEILING( (avg_grid%lon(2) - cosmo_lon(0)) / dlon )
     1665
     1666    jmin = FLOOR  ( (avg_grid%lat(1) - cosmo_lat(0)) / dlat )
     1667    jmax = CEILING( (avg_grid%lat(2) - cosmo_lat(0)) / dlat )
     1668   
     1669    message = "Grid " // TRIM(avg_grid%name(1)) // " averages over " // &
     1670              TRIM(str(imin)) // " <= i <= " // TRIM(str(imax)) //          &
     1671              " and " //                                                    &
     1672              TRIM(str(jmin)) // " <= j <= " // TRIM(str(jmax))
     1673    CALL report( 'get_cosmo_averaging_region', message )
     1674
     1675    nx = imax - imin + 1
     1676    ny = jmax - jmin + 1
     1677    avg_grid%n_columns = nx * ny
     1678
     1679    ALLOCATE( avg_grid%iii(avg_grid%n_columns),                         &
     1680              avg_grid%jjj(avg_grid%n_columns) )
     1681
     1682    l = 0
     1683    DO j = jmin, jmax
     1684    DO i = imin, imax
     1685       l = l + 1
     1686       avg_grid%iii(l) = i
     1687       avg_grid%jjj(l) = j
     1688    ENDDO
     1689    ENDDO
     1690
     1691 END SUBROUTINE get_cosmo_averaging_region
    16731692
    16741693
    … …  
    16831702!> 'modpoints'.
    16841703!------------------------------------------------------------------------------!
    1685     SUBROUTINE stretched_z(z, dz, dz_max, dz_stretch_factor, dz_stretch_level, &
    1686                            dz_stretch_level_start, dz_stretch_level_end,       &
    1687                            dz_stretch_factor_array)
    1688 
    1689        REAL(dp), DIMENSION(:), INTENT(INOUT) ::  z, dz, dz_stretch_factor_array
    1690        REAL(dp), DIMENSION(:), INTENT(INOUT) ::  dz_stretch_level_start, dz_stretch_level_end
    1691        REAL(dp), INTENT(IN) ::  dz_max, dz_stretch_factor, dz_stretch_level
    1692 
    1693        INTEGER ::  number_stretch_level_start        !< number of user-specified start levels for stretching
    1694        INTEGER ::  number_stretch_level_end          !< number of user-specified end levels for stretching
    1695 
    1696        REAL(dp), DIMENSION(:), ALLOCATABLE ::  min_dz_stretch_level_end
    1697        REAL(dp) ::  dz_level_end, dz_stretched
    1698 
    1699        INTEGER ::  dz_stretch_level_end_index(9)      !< vertical grid level index until which the vertical grid spacing is stretched
    1700        INTEGER ::  dz_stretch_level_start_index(9)    !< vertical grid level index above which the vertical grid spacing is stretched
    1701        INTEGER ::  dz_stretch_level_index = 0
    1702        INTEGER ::  k, n, number_dz
     1704 SUBROUTINE stretched_z(z, dz, dz_max, dz_stretch_factor, dz_stretch_level, &
     1705                        dz_stretch_level_start, dz_stretch_level_end,       &
     1706                        dz_stretch_factor_array)
     1707
     1708    REAL(wp), DIMENSION(:), INTENT(INOUT) ::  z, dz, dz_stretch_factor_array
     1709    REAL(wp), DIMENSION(:), INTENT(INOUT) ::  dz_stretch_level_start, dz_stretch_level_end
     1710    REAL(wp), INTENT(IN) ::  dz_max, dz_stretch_factor, dz_stretch_level
     1711
     1712    INTEGER ::  number_stretch_level_start        !< number of user-specified start levels for stretching
     1713    INTEGER ::  number_stretch_level_end          !< number of user-specified end levels for stretching
     1714
     1715    REAL(wp), DIMENSION(:), ALLOCATABLE ::  min_dz_stretch_level_end
     1716    REAL(wp) ::  dz_level_end, dz_stretched
     1717
     1718    INTEGER ::  dz_stretch_level_end_index(9)      !< vertical grid level index until which the vertical grid spacing is stretched
     1719    INTEGER ::  dz_stretch_level_start_index(9)    !< vertical grid level index above which the vertical grid spacing is stretched
     1720    INTEGER ::  dz_stretch_level_index = 0
     1721    INTEGER ::  k, n, number_dz
    17031722
    17041723!
    17051724!-- Compute height of u-levels from constant grid length and dz stretch factors
    1706        IF ( dz(1) == -1.0_dp )  THEN
    1707           message = 'missing dz'
    1708           CALL inifor_abort( 'stretched_z', message)
    1709        ELSEIF ( dz(1) <= 0.0_dp )  THEN
    1710           WRITE( message, * ) 'dz=', dz(1),' <= 0.0'
    1711           CALL inifor_abort( 'stretched_z', message)
    1712        ENDIF
     1725    IF ( dz(1) == -1.0_wp )  THEN
     1726       message = 'missing dz'
     1727       CALL inifor_abort( 'stretched_z', message)
     1728    ELSEIF ( dz(1) <= 0.0_wp )  THEN
     1729       WRITE( message, * ) 'dz=', dz(1),' <= 0.0'
     1730       CALL inifor_abort( 'stretched_z', message)
     1731    ENDIF
    17131732
    17141733!
    17151734!-- Initialize dz_stretch_level_start with the value of dz_stretch_level
    17161735!-- if it was set by the user
    1717        IF ( dz_stretch_level /= -9999999.9_dp ) THEN
    1718           dz_stretch_level_start(1) = dz_stretch_level
    1719        ENDIF
     1736    IF ( dz_stretch_level /= -9999999.9_wp ) THEN
     1737       dz_stretch_level_start(1) = dz_stretch_level
     1738    ENDIF
    17201739       
    17211740!
    … …  
    17271746!-- is used (Attention: The user is not allowed to specify a dz value equal
    17281747!-- to the default of dz_max = 999.0).
    1729        number_dz = COUNT( dz /= -1.0_dp .AND. dz /= dz_max )
    1730        number_stretch_level_start = COUNT( dz_stretch_level_start /=           &
    1731                                            -9999999.9_dp )
    1732        number_stretch_level_end = COUNT( dz_stretch_level_end /=               &
    1733                                          9999999.9_dp )
     1748    number_dz = COUNT( dz /= -1.0_wp .AND. dz /= dz_max )
     1749    number_stretch_level_start = COUNT( dz_stretch_level_start /=           &
     1750                                        -9999999.9_wp )
     1751    number_stretch_level_end = COUNT( dz_stretch_level_end /=               &
     1752                                      9999999.9_wp )
    17341753
    17351754!
    17361755!-- The number of specified end levels +1 has to be the same than the number
    17371756!-- of specified dz values
    1738        IF ( number_dz /= number_stretch_level_end + 1 ) THEN
    1739           WRITE( message, * ) 'The number of values for dz = ',                &
    1740                               number_dz, 'has to be the same than ',           &
    1741                               'the number of values for ',                     &
    1742                               'dz_stretch_level_end + 1 = ',                   &
    1743                               number_stretch_level_end+1
    1744           CALL inifor_abort( 'stretched_z', message)
    1745        ENDIF
     1757    IF ( number_dz /= number_stretch_level_end + 1 ) THEN
     1758       WRITE( message, * ) 'The number of values for dz = ',                &
     1759                           number_dz, 'has to be the same than ',           &
     1760                           'the number of values for ',                     &
     1761                           'dz_stretch_level_end + 1 = ',                   &
     1762                           number_stretch_level_end+1
     1763       CALL inifor_abort( 'stretched_z', message)
     1764    ENDIF
    17461765   
    17471766!
    1748 !--    The number of specified start levels has to be the same or one less than
    1749 !--    the number of specified dz values
    1750        IF ( number_dz /= number_stretch_level_start + 1 .AND.                  &
    1751             number_dz /= number_stretch_level_start ) THEN
    1752           WRITE( message, * ) 'The number of values for dz = ',         &
    1753                               number_dz, 'has to be the same or one ', &
    1754                               'more than& the number of values for ',  &
    1755                               'dz_stretch_level_start = ',             &
    1756                               number_stretch_level_start
    1757           CALL inifor_abort( 'stretched_z', message)
    1758        ENDIF
     1767!-- The number of specified start levels has to be the same or one less than
     1768!-- the number of specified dz values
     1769    IF ( number_dz /= number_stretch_level_start + 1 .AND.                  &
     1770         number_dz /= number_stretch_level_start ) THEN
     1771       WRITE( message, * ) 'The number of values for dz = ',         &
     1772                           number_dz, 'has to be the same or one ', &
     1773                           'more than& the number of values for ',  &
     1774                           'dz_stretch_level_start = ',             &
     1775                           number_stretch_level_start
     1776       CALL inifor_abort( 'stretched_z', message)
     1777    ENDIF
    17591778   
    1760 !--    The number of specified start levels has to be the same or one more than
    1761 !--    the number of specified end levels
    1762        IF ( number_stretch_level_start /= number_stretch_level_end + 1 .AND.   &
    1763             number_stretch_level_start /= number_stretch_level_end ) THEN
    1764           WRITE( message, * ) 'The number of values for ',              &
    1765                               'dz_stretch_level_start = ',              &
    1766                               dz_stretch_level_start, 'has to be the ',&
    1767                               'same or one more than& the number of ', &
    1768                               'values for dz_stretch_level_end = ',    &
    1769                               number_stretch_level_end
    1770           CALL inifor_abort( 'stretched_z', message)
    1771        ENDIF
     1779!-- The number of specified start levels has to be the same or one more than
     1780!-- the number of specified end levels
     1781    IF ( number_stretch_level_start /= number_stretch_level_end + 1 .AND.   &
     1782         number_stretch_level_start /= number_stretch_level_end ) THEN
     1783       WRITE( message, * ) 'The number of values for ',              &
     1784                           'dz_stretch_level_start = ',              &
     1785                           dz_stretch_level_start, 'has to be the ',&
     1786                           'same or one more than& the number of ', &
     1787                           'values for dz_stretch_level_end = ',    &
     1788                           number_stretch_level_end
     1789       CALL inifor_abort( 'stretched_z', message)
     1790    ENDIF
    17721791
    17731792!
    17741793!-- Initialize dz for the free atmosphere with the value of dz_max
    1775        IF ( dz(number_stretch_level_start+1) == -1.0_dp .AND.                     &
    1776             number_stretch_level_start /= 0 ) THEN
    1777           dz(number_stretch_level_start+1) = dz_max
    1778        ENDIF
     1794    IF ( dz(number_stretch_level_start+1) == -1.0_wp .AND.                     &
     1795         number_stretch_level_start /= 0 ) THEN
     1796       dz(number_stretch_level_start+1) = dz_max
     1797    ENDIF
    17791798       
    17801799!
    … …  
    17821801!-- atmosphere is desired (dz_stretch_level_end was not specified for the
    17831802!-- free atmosphere)
    1784        IF ( number_stretch_level_start == number_stretch_level_end + 1 ) THEN
    1785           dz_stretch_factor_array(number_stretch_level_start) =                   &
    1786           dz_stretch_factor
     1803    IF ( number_stretch_level_start == number_stretch_level_end + 1 )  THEN
     1804       dz_stretch_factor_array(number_stretch_level_start) =                   &
     1805       dz_stretch_factor
     1806    ENDIF
     1807
     1808!-- Allocation of arrays for stretching
     1809    ALLOCATE( min_dz_stretch_level_end(number_stretch_level_start) )
     1810
     1811!
     1812!-- The stretching region has to be large enough to allow for a smooth
     1813!-- transition between two different grid spacings
     1814    DO  n = 1, number_stretch_level_start
     1815       min_dz_stretch_level_end(n) = dz_stretch_level_start(n) +            &
     1816                                     4 * MAX( dz(n),dz(n+1) )
     1817    ENDDO
     1818
     1819    IF ( ANY( min_dz_stretch_level_end(1:number_stretch_level_start) >      &
     1820              dz_stretch_level_end(1:number_stretch_level_start) ) )  THEN
     1821    !IF ( ANY( min_dz_stretch_level_end >      &
     1822    !          dz_stretch_level_end ) ) THEN
     1823          message = 'Each dz_stretch_level_end has to be larger '  // &
     1824                    'than its corresponding value for ' //            &
     1825                    'dz_stretch_level_start + 4*MAX(dz(n),dz(n+1)) '//&
     1826                    'to allow for smooth grid stretching'
     1827          CALL inifor_abort('stretched_z', message)
     1828    ENDIF
     1829   
     1830!
     1831!-- Stretching must not be applied within the prandtl_layer
     1832!-- (first two grid points). For the default case dz_stretch_level_start
     1833!-- is negative. Therefore the absolut value is checked here.
     1834    IF ( ANY( ABS( dz_stretch_level_start ) < dz(1) * 1.5_wp ) )  THEN
     1835       WRITE( message, * ) 'Eeach dz_stretch_level_start has to be ',&
     1836                           'larger than ', dz(1) * 1.5
     1837          CALL inifor_abort( 'stretched_z', message)
     1838    ENDIF
     1839
     1840!
     1841!-- The stretching has to start and end on a grid level. Therefore
     1842!-- user-specified values have to ''interpolate'' to the next lowest level
     1843    IF ( number_stretch_level_start /= 0 )  THEN
     1844       dz_stretch_level_start(1) = INT( (dz_stretch_level_start(1) -        &
     1845                                         dz(1)/2.0) / dz(1) )               &
     1846                                   * dz(1) + dz(1)/2.0
     1847    ENDIF
     1848   
     1849    IF ( number_stretch_level_start > 1 )  THEN
     1850       DO  n = 2, number_stretch_level_start
     1851          dz_stretch_level_start(n) = INT( dz_stretch_level_start(n) /      &
     1852                                           dz(n) ) * dz(n)
     1853       ENDDO
     1854    ENDIF
     1855   
     1856    IF ( number_stretch_level_end /= 0 )  THEN
     1857       DO  n = 1, number_stretch_level_end
     1858          dz_stretch_level_end(n) = INT( dz_stretch_level_end(n) /          &
     1859                                         dz(n+1) ) * dz(n+1)
     1860       ENDDO
     1861    ENDIF
     1862 
     1863!
     1864!-- Determine stretching factor if necessary
     1865    IF ( number_stretch_level_end >= 1 )  THEN
     1866       CALL calculate_stretching_factor( number_stretch_level_end, dz,      &
     1867                                         dz_stretch_factor_array,           &   
     1868                                         dz_stretch_level_end,              &
     1869                                         dz_stretch_level_start )
     1870    ENDIF
     1871
     1872    z(1) = dz(1) * 0.5_wp
     1873!
     1874    dz_stretch_level_index = n
     1875    dz_stretched = dz(1)
     1876    DO  k = 2, n
     1877
     1878       IF ( dz_stretch_level <= z(k-1)  .AND.  dz_stretched < dz_max )  THEN
     1879
     1880          dz_stretched = dz_stretched * dz_stretch_factor
     1881          dz_stretched = MIN( dz_stretched, dz_max )
     1882
     1883          IF ( dz_stretch_level_index == n )  dz_stretch_level_index = k-1
     1884
    17871885       ENDIF
    17881886
    1789 !-- Allocation of arrays for stretching
    1790        ALLOCATE( min_dz_stretch_level_end(number_stretch_level_start) )
    1791 
    1792 !
    1793 !--    The stretching region has to be large enough to allow for a smooth
    1794 !--    transition between two different grid spacings
    1795        DO n = 1, number_stretch_level_start
    1796           min_dz_stretch_level_end(n) = dz_stretch_level_start(n) +            &
    1797                                         4 * MAX( dz(n),dz(n+1) )
    1798        ENDDO
    1799 
    1800        IF ( ANY( min_dz_stretch_level_end(1:number_stretch_level_start) >      &
    1801                  dz_stretch_level_end(1:number_stretch_level_start) ) ) THEN
    1802        !IF ( ANY( min_dz_stretch_level_end >      &
    1803        !          dz_stretch_level_end ) ) THEN
    1804              message = 'Each dz_stretch_level_end has to be larger '  // &
    1805                        'than its corresponding value for ' //            &
    1806                        'dz_stretch_level_start + 4*MAX(dz(n),dz(n+1)) '//&
    1807                        'to allow for smooth grid stretching'
    1808              CALL inifor_abort('stretched_z', message)
    1809        ENDIF
    1810        
    1811 !
    1812 !--    Stretching must not be applied within the prandtl_layer
    1813 !--    (first two grid points). For the default case dz_stretch_level_start
    1814 !--    is negative. Therefore the absolut value is checked here.
    1815        IF ( ANY( ABS( dz_stretch_level_start ) < dz(1) * 1.5_dp ) ) THEN
    1816           WRITE( message, * ) 'Eeach dz_stretch_level_start has to be ',&
    1817                               'larger than ', dz(1) * 1.5
    1818              CALL inifor_abort( 'stretched_z', message)
    1819        ENDIF
    1820 
    1821 !
    1822 !--    The stretching has to start and end on a grid level. Therefore
    1823 !--    user-specified values have to ''interpolate'' to the next lowest level
    1824        IF ( number_stretch_level_start /= 0 ) THEN
    1825           dz_stretch_level_start(1) = INT( (dz_stretch_level_start(1) -        &
    1826                                             dz(1)/2.0) / dz(1) )               &
    1827                                       * dz(1) + dz(1)/2.0
    1828        ENDIF
    1829        
    1830        IF ( number_stretch_level_start > 1 ) THEN
    1831           DO n = 2, number_stretch_level_start
    1832              dz_stretch_level_start(n) = INT( dz_stretch_level_start(n) /      &
    1833                                               dz(n) ) * dz(n)
    1834           ENDDO
    1835        ENDIF
    1836        
    1837        IF ( number_stretch_level_end /= 0 ) THEN
    1838           DO n = 1, number_stretch_level_end
    1839              dz_stretch_level_end(n) = INT( dz_stretch_level_end(n) /          &
    1840                                             dz(n+1) ) * dz(n+1)
    1841           ENDDO
    1842        ENDIF
    1843  
    1844 !
    1845 !--    Determine stretching factor if necessary
    1846        IF ( number_stretch_level_end >= 1 ) THEN
    1847           CALL calculate_stretching_factor( number_stretch_level_end, dz,      &
    1848                                             dz_stretch_factor_array,           &   
    1849                                             dz_stretch_level_end,              &
    1850                                             dz_stretch_level_start )
    1851        ENDIF
    1852 
    1853        z(1) = dz(1) * 0.5_dp
    1854 !
    1855        dz_stretch_level_index = n
    1856        dz_stretched = dz(1)
    1857        DO  k = 2, n
    1858 
    1859           IF ( dz_stretch_level <= z(k-1)  .AND.  dz_stretched < dz_max )  THEN
    1860 
    1861              dz_stretched = dz_stretched * dz_stretch_factor
    1862              dz_stretched = MIN( dz_stretched, dz_max )
    1863 
    1864              IF ( dz_stretch_level_index == n ) dz_stretch_level_index = k-1
    1865 
    1866           ENDIF
    1867 
    1868           z(k) = z(k-1) + dz_stretched
    1869 
    1870        ENDDO
    1871 !--    Determine u and v height levels considering the possibility of grid
    1872 !--    stretching in several heights.
    1873        n = 1
    1874        dz_stretch_level_start_index(:) = UBOUND(z, 1)
    1875        dz_stretch_level_end_index(:) = UBOUND(z, 1)
    1876        dz_stretched = dz(1)
    1877 
    1878 !--    The default value of dz_stretch_level_start is negative, thus the first
    1879 !--    condition is always true. Hence, the second condition is necessary.
    1880        DO  k = 2, UBOUND(z, 1)
    1881           IF ( dz_stretch_level_start(n) <= z(k-1) .AND.                      &
    1882                dz_stretch_level_start(n) /= -9999999.9_dp ) THEN
    1883              dz_stretched = dz_stretched * dz_stretch_factor_array(n)
    1884              
    1885              IF ( dz(n) > dz(n+1) ) THEN
    1886                 dz_stretched = MAX( dz_stretched, dz(n+1) ) !Restrict dz_stretched to the user-specified (higher) dz
    1887              ELSE
    1888                 dz_stretched = MIN( dz_stretched, dz(n+1) ) !Restrict dz_stretched to the user-specified (lower) dz
    1889              ENDIF
    1890              
    1891              IF ( dz_stretch_level_start_index(n) == UBOUND(z, 1) )            &
    1892              dz_stretch_level_start_index(n) = k-1
    1893              
     1887       z(k) = z(k-1) + dz_stretched
     1888
     1889    ENDDO
     1890!-- Determine u and v height levels considering the possibility of grid
     1891!-- stretching in several heights.
     1892    n = 1
     1893    dz_stretch_level_start_index(:) = UBOUND(z, 1)
     1894    dz_stretch_level_end_index(:) = UBOUND(z, 1)
     1895    dz_stretched = dz(1)
     1896
     1897!-- The default value of dz_stretch_level_start is negative, thus the first
     1898!-- condition is always true. Hence, the second condition is necessary.
     1899    DO  k = 2, UBOUND(z, 1)
     1900       IF ( dz_stretch_level_start(n) <= z(k-1) .AND.                      &
     1901            dz_stretch_level_start(n) /= -9999999.9_wp )  THEN
     1902          dz_stretched = dz_stretched * dz_stretch_factor_array(n)
     1903         
     1904          IF ( dz(n) > dz(n+1) )  THEN
     1905             dz_stretched = MAX( dz_stretched, dz(n+1) ) !Restrict dz_stretched to the user-specified (higher) dz
     1906          ELSE
     1907             dz_stretched = MIN( dz_stretched, dz(n+1) ) !Restrict dz_stretched to the user-specified (lower) dz
    18941908          ENDIF
    18951909         
    1896           z(k) = z(k-1) + dz_stretched
     1910          IF ( dz_stretch_level_start_index(n) == UBOUND(z, 1) )            &
     1911          dz_stretch_level_start_index(n) = k-1
    18971912         
    1898 !
    1899 !--       Make sure that the stretching ends exactly at dz_stretch_level_end
    1900           dz_level_end = ABS( z(k) - dz_stretch_level_end(n) )
    1901          
    1902           IF ( dz_level_end < dz(n+1)/3.0 ) THEN
    1903              z(k) = dz_stretch_level_end(n)
    1904              dz_stretched = dz(n+1)
    1905              dz_stretch_level_end_index(n) = k
    1906              n = n + 1             
    1907           ENDIF
    1908        ENDDO
    1909 
    1910        DEALLOCATE( min_dz_stretch_level_end )
    1911 
    1912     END SUBROUTINE stretched_z
     1913       ENDIF
     1914       
     1915       z(k) = z(k-1) + dz_stretched
     1916       
     1917!
     1918!--    Make sure that the stretching ends exactly at dz_stretch_level_end
     1919       dz_level_end = ABS( z(k) - dz_stretch_level_end(n) )
     1920       
     1921       IF ( dz_level_end < dz(n+1)/3.0 )  THEN
     1922          z(k) = dz_stretch_level_end(n)
     1923          dz_stretched = dz(n+1)
     1924          dz_stretch_level_end_index(n) = k
     1925          n = n + 1             
     1926       ENDIF
     1927    ENDDO
     1928
     1929    DEALLOCATE( min_dz_stretch_level_end )
     1930
     1931 END SUBROUTINE stretched_z
    19131932
    19141933
    … …  
    19281947                                         dz_stretch_level_start )
    19291948 
    1930     REAL(dp), DIMENSION(:), INTENT(IN)    ::  dz
    1931     REAL(dp), DIMENSION(:), INTENT(INOUT) ::  dz_stretch_factor_array
    1932     REAL(dp), DIMENSION(:), INTENT(IN)    ::  dz_stretch_level_end, dz_stretch_level_start
     1949    REAL(wp), DIMENSION(:), INTENT(IN)    ::  dz
     1950    REAL(wp), DIMENSION(:), INTENT(INOUT) ::  dz_stretch_factor_array
     1951    REAL(wp), DIMENSION(:), INTENT(IN)    ::  dz_stretch_level_end, dz_stretch_level_start
    19331952 
    19341953    INTEGER ::  iterations  !< number of iterations until stretch_factor_lower/upper_limit is reached 
    … …  
    19381957    INTEGER, INTENT(IN) ::  number_end !< number of user-specified end levels for stretching
    19391958       
    1940     REAL(dp) ::  delta_l               !< absolute difference between l and l_rounded
    1941     REAL(dp) ::  delta_stretch_factor  !< absolute difference between stretch_factor_1 and stretch_factor_2
    1942     REAL(dp) ::  delta_total_new       !< sum of delta_l and delta_stretch_factor for the next iteration (should be as small as possible)
    1943     REAL(dp) ::  delta_total_old       !< sum of delta_l and delta_stretch_factor for the last iteration
    1944     REAL(dp) ::  distance              !< distance between dz_stretch_level_start and dz_stretch_level_end (stretching region)
    1945     REAL(dp) ::  l                     !< value that fulfil Eq. (5) in the paper mentioned above together with stretch_factor_1 exactly
    1946     REAL(dp) ::  numerator             !< numerator of the quotient
    1947     REAL(dp) ::  stretch_factor_1      !< stretching factor that fulfil Eq. (5) togehter with l exactly
    1948     REAL(dp) ::  stretch_factor_2      !< stretching factor that fulfil Eq. (6) togehter with l_rounded exactly
     1959    REAL(wp) ::  delta_l               !< absolute difference between l and l_rounded
     1960    REAL(wp) ::  delta_stretch_factor  !< absolute difference between stretch_factor_1 and stretch_factor_2
     1961    REAL(wp) ::  delta_total_new       !< sum of delta_l and delta_stretch_factor for the next iteration (should be as small as possible)
     1962    REAL(wp) ::  delta_total_old       !< sum of delta_l and delta_stretch_factor for the last iteration
     1963    REAL(wp) ::  distance              !< distance between dz_stretch_level_start and dz_stretch_level_end (stretching region)
     1964    REAL(wp) ::  l                     !< value that fulfil Eq. (5) in the paper mentioned above together with stretch_factor_1 exactly
     1965    REAL(wp) ::  numerator             !< numerator of the quotient
     1966    REAL(wp) ::  stretch_factor_1      !< stretching factor that fulfil Eq. (5) togehter with l exactly
     1967    REAL(wp) ::  stretch_factor_2      !< stretching factor that fulfil Eq. (6) togehter with l_rounded exactly
    19491968   
    1950     REAL(dp) ::  dz_stretch_factor_array_2(9) = 1.08_dp  !< Array that contains all stretch_factor_2 that belongs to stretch_factor_1
     1969    REAL(wp) ::  dz_stretch_factor_array_2(9) = 1.08_wp  !< Array that contains all stretch_factor_2 that belongs to stretch_factor_1
    19511970   
    1952     REAL(dp), PARAMETER ::  stretch_factor_interval = 1.0E-06  !< interval for sampling possible stretching factors
    1953     REAL(dp), PARAMETER ::  stretch_factor_lower_limit = 0.88  !< lowest possible stretching factor
    1954     REAL(dp), PARAMETER ::  stretch_factor_upper_limit = 1.12  !< highest possible stretching factor
     1971    REAL(wp), PARAMETER ::  stretch_factor_interval = 1.0E-06  !< interval for sampling possible stretching factors
     1972    REAL(wp), PARAMETER ::  stretch_factor_lower_limit = 0.88  !< lowest possible stretching factor
     1973    REAL(wp), PARAMETER ::  stretch_factor_upper_limit = 1.12  !< highest possible stretching factor
    19551974 
    19561975 
    … …  
    19631982       delta_total_old = 1.0
    19641983       
    1965        IF ( dz(n) > dz(n+1) ) THEN
    1966           DO WHILE ( stretch_factor_1 >= stretch_factor_lower_limit )
     1984       IF ( dz(n) > dz(n+1) )  THEN
     1985          DO  WHILE ( stretch_factor_1 >= stretch_factor_lower_limit )
    19671986             
    19681987             stretch_factor_1 = 1.0 - iterations * stretch_factor_interval
    … …  
    19721991                         stretch_factor_1 - distance/dz(n)
    19731992             
    1974              IF ( numerator > 0.0 ) THEN
     1993             IF ( numerator > 0.0 )  THEN
    19751994                l = LOG( numerator ) / LOG( stretch_factor_1 ) - 1.0
    19761995                l_rounded = NINT( l )
    … …  
    19912010!--                stretch_factor_2 would guarantee that the stretched dz(n) is
    19922011!--                equal to dz(n+1) after l_rounded grid levels.
    1993              IF (delta_total_new < delta_total_old) THEN
     2012             IF (delta_total_new < delta_total_old)  THEN
    19942013                dz_stretch_factor_array(n) = stretch_factor_1
    19952014                dz_stretch_factor_array_2(n) = stretch_factor_2
    … …  
    20012020          ENDDO
    20022021             
    2003        ELSEIF ( dz(n) < dz(n+1) ) THEN
    2004           DO WHILE ( stretch_factor_1 <= stretch_factor_upper_limit )
     2022       ELSEIF ( dz(n) < dz(n+1) )  THEN
     2023          DO  WHILE ( stretch_factor_1 <= stretch_factor_upper_limit )
    20052024                     
    20062025             stretch_factor_1 = 1.0 + iterations * stretch_factor_interval
    … …  
    20272046!--          stretch_factor_2 would guarantee that the stretched dz(n) is
    20282047!--          equal to dz(n+1) after l_rounded grid levels.
    2029              IF (delta_total_new < delta_total_old) THEN
     2048             IF (delta_total_new < delta_total_old)  THEN
    20302049                dz_stretch_factor_array(n) = stretch_factor_1
    20312050                dz_stretch_factor_array_2(n) = stretch_factor_2
    … …  
    20452064!--    interval. If not, print a warning for the user.
    20462065       IF ( dz_stretch_factor_array_2(n) < stretch_factor_lower_limit .OR.     &
    2047             dz_stretch_factor_array_2(n) > stretch_factor_upper_limit ) THEN
    2048           WRITE( message, * ) 'stretch_factor_2 = ',                    &
     2066            dz_stretch_factor_array_2(n) > stretch_factor_upper_limit )  THEN
     2067          WRITE( message, * ) 'stretch_factor_2 = ',                           &
    20492068                                     dz_stretch_factor_array_2(n), ' which is',&
    20502069                                     ' responsible for exactly reaching& dz =',&
    … …  
    20672086!> coordinate vector 'z' and stores it in 'zw'.
    20682087!------------------------------------------------------------------------------!
    2069     SUBROUTINE midpoints(z, zw)
    2070 
    2071         REAL(dp), INTENT(IN)  ::  z(0:)
    2072         REAL(dp), INTENT(OUT) ::  zw(1:)
    2073 
    2074         INTEGER ::  k
    2075 
    2076         DO k = 1, UBOUND(zw, 1)
    2077            zw(k) = 0.5_dp * (z(k-1) + z(k))
    2078         ENDDO
    2079 
    2080     END SUBROUTINE midpoints
     2088 SUBROUTINE midpoints(z, zw)
     2089
     2090     REAL(wp), INTENT(IN)  ::  z(0:)
     2091     REAL(wp), INTENT(OUT) ::  zw(1:)
     2092
     2093     INTEGER ::  k
     2094
     2095     DO k = 1, UBOUND(zw, 1)
     2096        zw(k) = 0.5_wp * (z(k-1) + z(k))
     2097     ENDDO
     2098
     2099 END SUBROUTINE midpoints
    20812100
    20822101!------------------------------------------------------------------------------!
    … …  
    20852104!> Defines INFOR's IO groups.
    20862105!------------------------------------------------------------------------------!
    2087     SUBROUTINE setup_io_groups()
    2088 
    2089        INTEGER ::  ngroups
    2090 
    2091        ngroups = 16
    2092        ALLOCATE( io_group_list(ngroups) )
    2093 
    2094 !
    2095 !--    soil temp
    2096        io_group_list(1) = init_io_group(                                       &
    2097           in_files = soil_files,                                               &
    2098           out_vars = output_var_table(1:1),                                    &
    2099           in_var_list = input_var_table(1:1),                                  &
    2100           kind = 'soil-temperature'                                            &
    2101        )
    2102 
    2103 !
    2104 !--    soil water
    2105        io_group_list(2) = init_io_group(                                       &
    2106           in_files = soil_files,                                               &
    2107           out_vars = output_var_table(2:2),                                    &
    2108           in_var_list = input_var_table(2:2),                                  &
    2109           kind = 'soil-water'                                                  &
    2110        )
    2111 
    2112 !
    2113 !--    potential temperature, surface pressure, specific humidity including
    2114 !--    nudging and subsidence, and geostrophic winds ug, vg
    2115        io_group_list(3) = init_io_group(                                       &
    2116           in_files = flow_files,                                               &
    2117           out_vars = [output_var_table(56:64),                                 & ! internal averaged density and pressure profiles
    2118                       output_var_table(3:8), output_var_table(9:14),           &
    2119                       output_var_table(42:42), output_var_table(43:44),        &
    2120                       output_var_table(49:51), output_var_table(52:54)],       &
    2121           in_var_list = (/input_var_table(3), input_var_table(17),             & ! T, P, QV
    2122                           input_var_table(4) /),                               &
    2123           kind = 'thermodynamics',                                             &
    2124           n_output_quantities = 4                                              & ! P, Theta, Rho, qv
    2125        )
    2126 
    2127 !
    2128 !--    Moved to therodynamic io_group
    2129        !io_group_list(4) = init_io_group(                                       &
    2130        !   in_files = flow_files,                                               &
    2131        !   out_vars = [output_var_table(9:14), output_var_table(52:54)],        &
    2132        !   in_var_list = input_var_table(4:4),                                  &
    2133        !   kind = 'scalar'                                                      &
    2134        !)
    2135 
    2136 !
    2137 !--    u and v velocity
    2138        io_group_list(5) = init_io_group(                                       &
    2139           in_files = flow_files,                                               &
    2140           out_vars = [output_var_table(15:26), output_var_table(45:46)],       &
    2141           !out_vars = output_var_table(15:20),                                  &
    2142           in_var_list = input_var_table(5:6),                                  &
    2143           !in_var_list = input_var_table(5:5),                                  &
    2144           kind = 'velocities'                                                  &
    2145        )
     2106 SUBROUTINE setup_io_groups()
     2107
     2108    INTEGER ::  ngroups
     2109
     2110    ngroups = 16
     2111    ALLOCATE( io_group_list(ngroups) )
     2112
     2113!
     2114!-- soil temp
     2115    io_group_list(1) = init_io_group(                                       &
     2116       in_files = soil_files,                                               &
     2117       out_vars = output_var_table(1:1),                                    &
     2118       in_var_list = input_var_table(1:1),                                  &
     2119       kind = 'soil-temperature'                                            &
     2120    )
     2121
     2122!
     2123!-- soil water
     2124    io_group_list(2) = init_io_group(                                       &
     2125       in_files = soil_files,                                               &
     2126       out_vars = output_var_table(2:2),                                    &
     2127       in_var_list = input_var_table(2:2),                                  &
     2128       kind = 'soil-water'                                                  &
     2129    )
     2130
     2131!
     2132!-- potential temperature, surface pressure, specific humidity including
     2133!-- nudging and subsidence, and geostrophic winds ug, vg
     2134    io_group_list(3) = init_io_group(                                       &
     2135       in_files = flow_files,                                               &
     2136       out_vars = [output_var_table(56:64),                                 & ! internal averaged density and pressure profiles
     2137                   output_var_table(3:8), output_var_table(9:14),           &
     2138                   output_var_table(42:42), output_var_table(43:44),        &
     2139                   output_var_table(49:51), output_var_table(52:54)],       &
     2140       in_var_list = (/input_var_table(3), input_var_table(17),             & ! T, P, QV
     2141                       input_var_table(4) /),                               &
     2142       kind = 'thermodynamics',                                             &
     2143       n_output_quantities = 4                                              & ! P, Theta, Rho, qv
     2144    )
     2145
     2146!
     2147!-- Moved to therodynamic io_group
     2148    !io_group_list(4) = init_io_group(                                       &
     2149    !   in_files = flow_files,                                               &
     2150    !   out_vars = [output_var_table(9:14), output_var_table(52:54)],        &
     2151    !   in_var_list = input_var_table(4:4),                                  &
     2152    !   kind = 'scalar'                                                      &
     2153    !)
     2154
     2155!
     2156!-- u and v velocity
     2157    io_group_list(5) = init_io_group(                                       &
     2158       in_files = flow_files,                                               &
     2159       out_vars = [output_var_table(15:26), output_var_table(45:46)],       &
     2160       !out_vars = output_var_table(15:20),                                  &
     2161       in_var_list = input_var_table(5:6),                                  &
     2162       !in_var_list = input_var_table(5:5),                                  &
     2163       kind = 'velocities'                                                  &
     2164    )
    21462165   
    21472166!
    2148 !--    v velocity, deprecated!
    2149        !io_group_list(6) = init_io_group(                                       &
    2150        !   in_files = flow_files,                                               &
    2151        !   out_vars = output_var_table(21:26),                                  &
    2152        !   in_var_list = input_var_table(6:6),                                  &
    2153        !   kind = 'horizontal velocity'                                         &
    2154        !)
    2155        !io_group_list(6) % to_be_processed = .FALSE.
     2167!-- v velocity, deprecated!
     2168    !io_group_list(6) = init_io_group(                                       &
     2169    !   in_files = flow_files,                                               &
     2170    !   out_vars = output_var_table(21:26),                                  &
     2171    !   in_var_list = input_var_table(6:6),                                  &
     2172    !   kind = 'horizontal velocity'                                         &
     2173    !)
     2174    !io_group_list(6)%to_be_processed = .FALSE.
    21562175   
    21572176!
    2158 !--    w velocity and subsidence and w nudging
    2159        io_group_list(7) = init_io_group(                                       &
    2160           in_files = flow_files,                                               &
    2161           out_vars = [output_var_table(27:32), output_var_table(47:48)],       &
    2162           in_var_list = input_var_table(7:7),                                  &
    2163           kind = 'scalar'                                                      &
    2164        )
    2165 !
    2166 !--    rain
    2167        io_group_list(8) = init_io_group(                                       &
    2168           in_files = soil_moisture_files,                                      &
    2169           out_vars = output_var_table(33:33),                                  &
    2170           in_var_list = input_var_table(8:8),                                  &
    2171           kind = 'accumulated'                                                 &
    2172        )
    2173        io_group_list(8) % to_be_processed = .FALSE.
    2174 !
    2175 !--    snow
    2176        io_group_list(9) = init_io_group(                                       &
    2177           in_files = soil_moisture_files,                                      &
    2178           out_vars = output_var_table(34:34),                                  &
    2179           in_var_list = input_var_table(9:9),                                  &
    2180           kind = 'accumulated'                                                 &
    2181        )
    2182        io_group_list(9) % to_be_processed = .FALSE.
    2183 !
    2184 !--    graupel
    2185        io_group_list(10) = init_io_group(                                      &
    2186           in_files = soil_moisture_files,                                      &
    2187           out_vars = output_var_table(35:35),                                  &
    2188           in_var_list = input_var_table(10:10),                                &
    2189           kind = 'accumulated'                                                 &
    2190        )
    2191        io_group_list(10) % to_be_processed = .FALSE.
    2192 !
    2193 !--    evapotranspiration
    2194        io_group_list(11) = init_io_group(                                      &
    2195           in_files = soil_moisture_files,                                      &
    2196           out_vars = output_var_table(37:37),                                  &
    2197           in_var_list = input_var_table(11:11),                                &
    2198           kind = 'accumulated'                                                 &
    2199        )
    2200        io_group_list(11) % to_be_processed = .FALSE.
    2201 !
    2202 !--    2m air temperature
    2203        io_group_list(12) = init_io_group(                                      &
    2204           in_files = soil_moisture_files,                                      &
    2205           out_vars = output_var_table(36:36),                                  &
    2206           in_var_list = input_var_table(12:12),                                &
    2207           kind = 'surface'                                                     &
    2208        )
    2209        io_group_list(12) % to_be_processed = .FALSE.
    2210 !
    2211 !--    incoming diffusive sw flux
    2212        io_group_list(13) = init_io_group(                                      &
    2213           in_files = radiation_files,                                          &
    2214           out_vars = output_var_table(38:38),                                  &
    2215           in_var_list = input_var_table(13:13),                                &
    2216           kind = 'running average'                                             &
    2217        )
    2218        io_group_list(13) % to_be_processed = .FALSE.
    2219 !
    2220 !--    incoming direct sw flux
    2221        io_group_list(14) = init_io_group(                                      &
    2222           in_files = radiation_files,                                          &
    2223           out_vars = output_var_table(39:39),                                  &
    2224           in_var_list = input_var_table(14:14),                                &
    2225           kind = 'running average'                                             &
    2226        )
    2227        io_group_list(14) % to_be_processed = .FALSE.
    2228 !
    2229 !--    sw radiation balance
    2230        io_group_list(15) = init_io_group(                                      &
    2231           in_files = radiation_files,                                          &
    2232           out_vars = output_var_table(40:40),                                  &
    2233           in_var_list = input_var_table(15:15),                                &
    2234           kind = 'running average'                                             &
    2235        )
    2236        io_group_list(15) % to_be_processed = .FALSE.
    2237 !
    2238 !--    lw radiation balance
    2239        io_group_list(16) = init_io_group(                                      &
    2240           in_files = radiation_files,                                          &
    2241           out_vars = output_var_table(41:41),                                  &
    2242           in_var_list = input_var_table(16:16),                                &
    2243           kind = 'running average'                                             &
    2244        )
    2245        io_group_list(16) % to_be_processed = .FALSE.
    2246 
    2247     END SUBROUTINE setup_io_groups
     2177!-- w velocity and subsidence and w nudging
     2178    io_group_list(7) = init_io_group(                                       &
     2179       in_files = flow_files,                                               &
     2180       out_vars = [output_var_table(27:32), output_var_table(47:48)],       &
     2181       in_var_list = input_var_table(7:7),                                  &
     2182       kind = 'scalar'                                                      &
     2183    )
     2184!
     2185!-- rain
     2186    io_group_list(8) = init_io_group(                                       &
     2187       in_files = soil_moisture_files,                                      &
     2188       out_vars = output_var_table(33:33),                                  &
     2189       in_var_list = input_var_table(8:8),                                  &
     2190       kind = 'accumulated'                                                 &
     2191    )
     2192    io_group_list(8)%to_be_processed = .FALSE.
     2193!
     2194!-- snow
     2195    io_group_list(9) = init_io_group(                                       &
     2196       in_files = soil_moisture_files,                                      &
     2197       out_vars = output_var_table(34:34),                                  &
     2198       in_var_list = input_var_table(9:9),                                  &
     2199       kind = 'accumulated'                                                 &
     2200    )
     2201    io_group_list(9)%to_be_processed = .FALSE.
     2202!
     2203!-- graupel
     2204    io_group_list(10) = init_io_group(                                      &
     2205       in_files = soil_moisture_files,                                      &
     2206       out_vars = output_var_table(35:35),                                  &
     2207       in_var_list = input_var_table(10:10),                                &
     2208       kind = 'accumulated'                                                 &
     2209    )
     2210    io_group_list(10)%to_be_processed = .FALSE.
     2211!
     2212!-- evapotranspiration
     2213    io_group_list(11) = init_io_group(                                      &
     2214       in_files = soil_moisture_files,                                      &
     2215       out_vars = output_var_table(37:37),                                  &
     2216       in_var_list = input_var_table(11:11),                                &
     2217       kind = 'accumulated'                                                 &
     2218    )
     2219    io_group_list(11)%to_be_processed = .FALSE.
     2220!
     2221!-- 2m air temperature
     2222    io_group_list(12) = init_io_group(                                      &
     2223       in_files = soil_moisture_files,                                      &
     2224       out_vars = output_var_table(36:36),                                  &
     2225       in_var_list = input_var_table(12:12),                                &
     2226       kind = 'surface'                                                     &
     2227    )
     2228    io_group_list(12)%to_be_processed = .FALSE.
     2229!
     2230!-- incoming diffusive sw flux
     2231    io_group_list(13) = init_io_group(                                      &
     2232       in_files = radiation_files,                                          &
     2233       out_vars = output_var_table(38:38),                                  &
     2234       in_var_list = input_var_table(13:13),                                &
     2235       kind = 'running average'                                             &
     2236    )
     2237    io_group_list(13)%to_be_processed = .FALSE.
     2238!
     2239!-- incoming direct sw flux
     2240    io_group_list(14) = init_io_group(                                      &
     2241       in_files = radiation_files,                                          &
     2242       out_vars = output_var_table(39:39),                                  &
     2243       in_var_list = input_var_table(14:14),                                &
     2244       kind = 'running average'                                             &
     2245    )
     2246    io_group_list(14)%to_be_processed = .FALSE.
     2247!
     2248!-- sw radiation balance
     2249    io_group_list(15) = init_io_group(                                      &
     2250       in_files = radiation_files,                                          &
     2251       out_vars = output_var_table(40:40),                                  &
     2252       in_var_list = input_var_table(15:15),                                &
     2253       kind = 'running average'                                             &
     2254    )
     2255    io_group_list(15)%to_be_processed = .FALSE.
     2256!
     2257!-- lw radiation balance
     2258    io_group_list(16) = init_io_group(                                      &
     2259       in_files = radiation_files,                                          &
     2260       out_vars = output_var_table(41:41),                                  &
     2261       in_var_list = input_var_table(16:16),                                &
     2262       kind = 'running average'                                             &
     2263    )
     2264    io_group_list(16)%to_be_processed = .FALSE.
     2265
     2266 END SUBROUTINE setup_io_groups
    22482267
    22492268
    … …  
    22592278!> on the output quantity Theta.
    22602279!------------------------------------------------------------------------------!
    2261     FUNCTION init_io_group(in_files, out_vars, in_var_list, kind,              &
    2262                            n_output_quantities) RESULT(group)
    2263        CHARACTER(LEN=PATH), INTENT(IN) ::  in_files(:)
    2264        CHARACTER(LEN=*), INTENT(IN)    ::  kind
    2265        TYPE(nc_var), INTENT(IN)        ::  out_vars(:)
    2266        TYPE(nc_var), INTENT(IN)        ::  in_var_list(:)
    2267        INTEGER, OPTIONAL               ::  n_output_quantities
    2268 
    2269        TYPE(io_group)                  ::  group
    2270 
    2271        group % nt = SIZE(in_files)
    2272        group % nv = SIZE(out_vars)
    2273        group % n_inputs = SIZE(in_var_list)
    2274        group % kind = TRIM(kind)
    2275 !
    2276 !--    For the 'thermodynamics' IO group, one quantity more than input variables
    2277 !--    is needed to compute all output variables of the IO group. Concretely, in
    2278 !--    preprocess() the density is computed from T,P or PP,QV in adddition to
    2279 !--    the variables Theta, p, qv. In read_input_variables(),
    2280 !--    n_output_quantities is used to allocate the correct number of input
    2281 !--    buffers.
    2282        IF ( PRESENT(n_output_quantities) )  THEN
    2283           group % n_output_quantities = n_output_quantities
    2284        ELSE
    2285           group % n_output_quantities = group % n_inputs
    2286        ENDIF
    2287 
    2288        ALLOCATE(group % in_var_list(group % n_inputs))
    2289        ALLOCATE(group % in_files(group % nt))
    2290        ALLOCATE(group % out_vars(group % nv))
    2291 
    2292        group % in_var_list = in_var_list
    2293        group % in_files = in_files
    2294        group % out_vars = out_vars
    2295        group % to_be_processed = .TRUE.
    2296 
    2297     END FUNCTION init_io_group
     2280 FUNCTION init_io_group(in_files, out_vars, in_var_list, kind,              &
     2281                        n_output_quantities) RESULT(group)
     2282    CHARACTER(LEN=PATH), INTENT(IN) ::  in_files(:)
     2283    CHARACTER(LEN=*), INTENT(IN)    ::  kind
     2284    TYPE(nc_var), INTENT(IN)        ::  out_vars(:)
     2285    TYPE(nc_var), INTENT(IN)        ::  in_var_list(:)
     2286    INTEGER, OPTIONAL               ::  n_output_quantities
     2287
     2288    TYPE(io_group)                  ::  group
     2289
     2290    group%nt = SIZE(in_files)
     2291    group%nv = SIZE(out_vars)
     2292    group%n_inputs = SIZE(in_var_list)
     2293    group%kind = TRIM(kind)
     2294!
     2295!-- For the 'thermodynamics' IO group, one quantity more than input variables
     2296!-- is needed to compute all output variables of the IO group. Concretely, in
     2297!-- preprocess() the density is computed from T,P or PP,QV in adddition to
     2298!-- the variables Theta, p, qv. In read_input_variables(),
     2299!-- n_output_quantities is used to allocate the correct number of input
     2300!-- buffers.
     2301    IF ( PRESENT(n_output_quantities) )  THEN
     2302       group%n_output_quantities = n_output_quantities
     2303    ELSE
     2304       group%n_output_quantities = group%n_inputs
     2305    ENDIF
     2306
     2307    ALLOCATE(group%in_var_list(group%n_inputs))
     2308    ALLOCATE(group%in_files(group%nt))
     2309    ALLOCATE(group%out_vars(group%nv))
     2310
     2311    group%in_var_list = in_var_list
     2312    group%in_files = in_files
     2313    group%out_vars = out_vars
     2314    group%to_be_processed = .TRUE.
     2315
     2316 END FUNCTION init_io_group
    22982317
    22992318
    … …  
    23032322!> Deallocates all allocated variables.
    23042323!------------------------------------------------------------------------------!
    2305     SUBROUTINE fini_grids()
    2306 
    2307        CALL report('fini_grids', 'Deallocating grids', cfg % debug)
    2308        
    2309        DEALLOCATE(x, y, z, xu, yv, zw, z_column, zw_column)
    2310 
    2311        DEALLOCATE(palm_grid%x,  palm_grid%y,  palm_grid%z,                     &
    2312                   palm_grid%xu, palm_grid%yv, palm_grid%zw,                    &
    2313                   palm_grid%clon,  palm_grid%clat,                             &
    2314                   palm_grid%clonu, palm_grid%clatu)
    2315 
    2316        DEALLOCATE(palm_intermediate%x,  palm_intermediate%y,  palm_intermediate%z, &
    2317                   palm_intermediate%xu, palm_intermediate%yv, palm_intermediate%zw,&
    2318                   palm_intermediate%clon,  palm_intermediate%clat,             & 
    2319                   palm_intermediate%clonu, palm_intermediate%clatu)
    2320 
    2321        DEALLOCATE(cosmo_grid%lon,  cosmo_grid%lat,                             &
    2322                   cosmo_grid%lonu, cosmo_grid%latv,                            &
    2323                   cosmo_grid%hfl)
    2324 
    2325     END SUBROUTINE fini_grids
     2324 SUBROUTINE fini_grids()
     2325
     2326    CALL report('fini_grids', 'Deallocating grids', cfg%debug)
     2327   
     2328    DEALLOCATE(x, y, z, xu, yv, zw, z_column, zw_column)
     2329
     2330    DEALLOCATE(palm_grid%x,  palm_grid%y,  palm_grid%z,                        &
     2331               palm_grid%xu, palm_grid%yv, palm_grid%zw,                       &
     2332               palm_grid%clon,  palm_grid%clat,                                &
     2333               palm_grid%clonu, palm_grid%clatu)
     2334
     2335    DEALLOCATE(palm_intermediate%x,  palm_intermediate%y,  palm_intermediate%z, &
     2336               palm_intermediate%xu, palm_intermediate%yv, palm_intermediate%zw,&
     2337               palm_intermediate%clon,  palm_intermediate%clat,                & 
     2338               palm_intermediate%clonu, palm_intermediate%clatu)
     2339
     2340    DEALLOCATE(cosmo_grid%lon,  cosmo_grid%lat,                                &
     2341               cosmo_grid%lonu, cosmo_grid%latv,                               &
     2342               cosmo_grid%hfl)
     2343
     2344 END SUBROUTINE fini_grids
    23262345
    23272346
    … …  
    23312350!> Initializes the variable list.
    23322351!------------------------------------------------------------------------------!
    2333     SUBROUTINE setup_variable_tables(ic_mode)
    2334        CHARACTER(LEN=*), INTENT(IN) ::  ic_mode
    2335        INTEGER                      ::  n_invar = 0  !< number of variables in the input variable table
    2336        INTEGER                      ::  n_outvar = 0 !< number of variables in the output variable table
    2337        TYPE(nc_var), POINTER        ::  var
    2338 
    2339        IF (TRIM(cfg % start_date) == '')  THEN
    2340           message = 'Simulation start date has not been set.'
    2341           CALL inifor_abort('setup_variable_tables', message)
    2342        ENDIF
    2343 
    2344        nc_source_text = 'COSMO-DE analysis from ' // TRIM(cfg % start_date)
    2345 
    2346        n_invar = 17
    2347        n_outvar = 64
    2348        ALLOCATE( input_var_table(n_invar) )
    2349        ALLOCATE( output_var_table(n_outvar) )
     2352 SUBROUTINE setup_variable_tables(ic_mode)
     2353    CHARACTER(LEN=*), INTENT(IN) ::  ic_mode
     2354    INTEGER                      ::  n_invar = 0  !< number of variables in the input variable table
     2355    INTEGER                      ::  n_outvar = 0 !< number of variables in the output variable table
     2356    TYPE(nc_var), POINTER        ::  var
     2357
     2358    IF (TRIM(cfg%start_date) == '')  THEN
     2359       message = 'Simulation start date has not been set.'
     2360       CALL inifor_abort('setup_variable_tables', message)
     2361    ENDIF
     2362
     2363    nc_source_text = 'COSMO-DE