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

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

File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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 analysis from ' // TRIM(cfg%start_date)
     2364
     2365    n_invar = 17
     2366    n_outvar = 64
     2367    ALLOCATE( input_var_table(n_invar) )
     2368    ALLOCATE( output_var_table(n_outvar) )
    23502369
    23512370!
     
    23532372!- Section 1: NetCDF input variables
    23542373!------------------------------------------------------------------------------
    2355        var => input_var_table(1)
    2356        var % name = 'T_SO'
    2357        var % to_be_processed = .TRUE.
    2358        var % is_upside_down = .FALSE.
    2359 
    2360        var => input_var_table(2)
    2361        var % name = 'W_SO'
    2362        var % to_be_processed = .TRUE.
    2363        var % is_upside_down = .FALSE.
    2364 
    2365        var => input_var_table(3)
    2366        var % name = 'T'
    2367        var % to_be_processed = .TRUE.
    2368        var % is_upside_down = .TRUE.
    2369 
    2370        var => input_var_table(4)
    2371        var % name = 'QV'
    2372        var % to_be_processed = .TRUE.
    2373        var % is_upside_down = .TRUE.
    2374 
    2375        var => input_var_table(5)
    2376        var % name = 'U'
    2377        var % to_be_processed = .TRUE.
    2378        var % is_upside_down = .TRUE.
    2379 
    2380        var => input_var_table(6)
    2381        var % name = 'V'
    2382        var % to_be_processed = .TRUE.
    2383        var % is_upside_down = .TRUE.
    2384 
    2385        var => input_var_table(7)
    2386        var % name = 'W'
    2387        var % to_be_processed = .TRUE.
    2388        var % is_upside_down = .TRUE.
    2389 
    2390        var => input_var_table(8)
    2391        var % name = 'RAIN_GSP'
    2392        var % to_be_processed = .TRUE.
    2393        var % is_upside_down = .FALSE.
    2394 
    2395        var => input_var_table(9)
    2396        var % name = 'SNOW_GSP'
    2397        var % to_be_processed = .TRUE.
    2398        var % is_upside_down = .FALSE.
    2399 
    2400        var => input_var_table(10)
    2401        var % name = 'GRAU_GSP'
    2402        var % to_be_processed = .TRUE.
    2403        var % is_upside_down = .FALSE.
    2404 
    2405        var => input_var_table(11)
    2406        var % name = 'AEVAP_S'
    2407        var % to_be_processed = .TRUE.
    2408        var % is_upside_down = .FALSE.
    2409 
    2410        var => input_var_table(12)
    2411        var % name = 'T_2M'
    2412        var % to_be_processed = .TRUE.
    2413        var % is_upside_down = .FALSE.
    2414 
    2415        var => input_var_table(13)
    2416        var % name = 'ASWDIFD_S'
    2417        var % to_be_processed = .TRUE.
    2418        var % is_upside_down = .FALSE.
    2419 
    2420        var => input_var_table(14)
    2421        var % name = 'ASWDIR_S'
    2422        var % to_be_processed = .TRUE.
    2423        var % is_upside_down = .FALSE.
    2424 
    2425        var => input_var_table(15)
    2426        var % name = 'ASOB_S'
    2427        var % to_be_processed = .TRUE.
    2428        var % is_upside_down = .FALSE.
    2429 
    2430        var => input_var_table(16)
    2431        var % name = 'ATHB_S'
    2432        var % to_be_processed = .TRUE.
    2433        var % is_upside_down = .FALSE.
    2434 
    2435        var => input_var_table(17)
    2436        var % name = 'P'
    2437        var % to_be_processed = .TRUE.
    2438        var % is_upside_down = .TRUE.
     2374    var => input_var_table(1)
     2375    var%name = 'T_SO'
     2376    var%to_be_processed = .TRUE.
     2377    var%is_upside_down = .FALSE.
     2378
     2379    var => input_var_table(2)
     2380    var%name = 'W_SO'
     2381    var%to_be_processed = .TRUE.
     2382    var%is_upside_down = .FALSE.
     2383
     2384    var => input_var_table(3)
     2385    var%name = 'T'
     2386    var%to_be_processed = .TRUE.
     2387    var%is_upside_down = .TRUE.
     2388
     2389    var => input_var_table(4)
     2390    var%name = 'QV'
     2391    var%to_be_processed = .TRUE.
     2392    var%is_upside_down = .TRUE.
     2393
     2394    var => input_var_table(5)
     2395    var%name = 'U'
     2396    var%to_be_processed = .TRUE.
     2397    var%is_upside_down = .TRUE.
     2398
     2399    var => input_var_table(6)
     2400    var%name = 'V'
     2401    var%to_be_processed = .TRUE.
     2402    var%is_upside_down = .TRUE.
     2403
     2404    var => input_var_table(7)
     2405    var%name = 'W'
     2406    var%to_be_processed = .TRUE.
     2407    var%is_upside_down = .TRUE.
     2408
     2409    var => input_var_table(8)
     2410    var%name = 'RAIN_GSP'
     2411    var%to_be_processed = .TRUE.
     2412    var%is_upside_down = .FALSE.
     2413
     2414    var => input_var_table(9)
     2415    var%name = 'SNOW_GSP'
     2416    var%to_be_processed = .TRUE.
     2417    var%is_upside_down = .FALSE.
     2418
     2419    var => input_var_table(10)
     2420    var%name = 'GRAU_GSP'
     2421    var%to_be_processed = .TRUE.
     2422    var%is_upside_down = .FALSE.
     2423
     2424    var => input_var_table(11)
     2425    var%name = 'AEVAP_S'
     2426    var%to_be_processed = .TRUE.
     2427    var%is_upside_down = .FALSE.
     2428
     2429    var => input_var_table(12)
     2430    var%name = 'T_2M'
     2431    var%to_be_processed = .TRUE.
     2432    var%is_upside_down = .FALSE.
     2433
     2434    var => input_var_table(13)
     2435    var%name = 'ASWDIFD_S'
     2436    var%to_be_processed = .TRUE.
     2437    var%is_upside_down = .FALSE.
     2438
     2439    var => input_var_table(14)
     2440    var%name = 'ASWDIR_S'
     2441    var%to_be_processed = .TRUE.
     2442    var%is_upside_down = .FALSE.
     2443
     2444    var => input_var_table(15)
     2445    var%name = 'ASOB_S'
     2446    var%to_be_processed = .TRUE.
     2447    var%is_upside_down = .FALSE.
     2448
     2449    var => input_var_table(16)
     2450    var%name = 'ATHB_S'
     2451    var%to_be_processed = .TRUE.
     2452    var%is_upside_down = .FALSE.
     2453
     2454    var => input_var_table(17)
     2455    var%name = 'P'
     2456    var%to_be_processed = .TRUE.
     2457    var%is_upside_down = .TRUE.
    24392458
    24402459!
     
    24462465! Section 2.1: Realistic forcings, i.e. 3D initial and boundary conditions
    24472466!------------------------------------------------------------------------------
    2448        output_var_table(1) = init_nc_var(                                      &
    2449           name              = 'init_soil_t',                                   &
    2450           std_name          = "",                                              &
    2451           long_name         = "initial soil temperature",                      &
    2452           units             = "K",                                             &
    2453           kind              = "init soil",                                     &
    2454           input_id          = 1,                                               &
    2455           output_file       = output_file,                                     &
    2456           grid              = palm_grid,                                       &
    2457           intermediate_grid = palm_intermediate                                &
    2458        )
    2459 
    2460        output_var_table(2) = init_nc_var(                                      &
    2461           name              = 'init_soil_m',                                   &
    2462           std_name          = "",                                              &
    2463           long_name         = "initial soil moisture",                         &
    2464           units             = "m^3/m^3",                                       &
    2465           kind              = "init soil",                                     &
    2466           input_id          = 1,                                               &
    2467           output_file       = output_file,                                     &
    2468           grid              = palm_grid,                                       &
    2469           intermediate_grid = palm_intermediate                                &
    2470        )
    2471 
    2472        output_var_table(3) = init_nc_var(                                      &
    2473           name              = 'init_atmosphere_pt',                            &
    2474           std_name          = "",                                              &
    2475           long_name         = "initial potential temperature",                 &
    2476           units             = "K",                                             &
    2477           kind              = "init scalar",                                   &
    2478           input_id          = 1,                                               & ! first in (T, p) IO group
    2479           output_file       = output_file,                                     &
    2480           grid              = palm_grid,                                       &
    2481           intermediate_grid = palm_intermediate,                               &
    2482           is_profile = (TRIM(ic_mode) == 'profile')                            &
    2483        )
    2484        IF (TRIM(ic_mode) == 'profile')  THEN
    2485           output_var_table(3) % averaging_grid => averaged_initial_scalar_profile
    2486        ENDIF
    2487 
    2488        output_var_table(4) = init_nc_var(                                      &
    2489           name              = 'ls_forcing_left_pt',                            &
    2490           std_name          = "",                                              &
    2491           long_name         = "large-scale forcing for left model boundary for the potential temperature", &
    2492           units             = "K",                                             &
    2493           kind              = "left scalar",                                   &
    2494           input_id          = 1,                                               &
    2495           grid              = scalars_west_grid,                               &
    2496           intermediate_grid = scalars_west_intermediate,                       &
    2497           output_file = output_file                                            &
    2498        )
    2499 
    2500        output_var_table(5) = init_nc_var(                                      &
    2501           name              = 'ls_forcing_right_pt',                           &
    2502           std_name          = "",                                              &
    2503           long_name         = "large-scale forcing for right model boundary for the potential temperature", &
    2504           units             = "K",                                             &
    2505           kind              = "right scalar",                                  &
    2506           input_id          = 1,                                               &
    2507           grid              = scalars_east_grid,                               &
    2508           intermediate_grid = scalars_east_intermediate,                       &
    2509           output_file = output_file                                            &
    2510        )
    2511 
    2512        output_var_table(6) = init_nc_var(                                      &
    2513           name              = 'ls_forcing_north_pt',                           &
    2514           std_name          = "",                                              &
    2515           long_name         = "large-scale forcing for north model boundary for the potential temperature", &
    2516           units             = "K",                                             &
    2517           kind              = "north scalar",                                  &
    2518           input_id          = 1,                                               &
    2519           grid              = scalars_north_grid,                              &
    2520           intermediate_grid = scalars_north_intermediate,                      &
    2521           output_file = output_file                                            &
    2522        )
    2523 
    2524        output_var_table(7) = init_nc_var(                                      &
    2525           name              = 'ls_forcing_south_pt',                           &
    2526           std_name          = "",                                              &
    2527           long_name         = "large-scale forcing for south model boundary for the potential temperature", &
    2528           units             = "K",                                             &
    2529           kind              = "south scalar",                                  &
    2530           input_id          = 1,                                               &
    2531           grid              = scalars_south_grid,                              &
    2532           intermediate_grid = scalars_south_intermediate,                      &
    2533           output_file = output_file                                            &
    2534        )
    2535 
    2536        output_var_table(8) = init_nc_var(                                      &
    2537           name              = 'ls_forcing_top_pt',                             &
    2538           std_name          = "",                                              &
    2539           long_name         = "large-scale forcing for top model boundary for the potential temperature", &
    2540           units             = "K",                                             &
    2541           kind              = "top scalar",                                    &
    2542           input_id          = 1,                                               &
    2543           grid              = scalars_top_grid,                                &
    2544           intermediate_grid = scalars_top_intermediate,                        &
    2545           output_file = output_file                                            &
    2546        )
    2547 
    2548        output_var_table(9) = init_nc_var(                                      &
    2549           name              = 'init_atmosphere_qv',                            &
    2550           std_name          = "",                                              &
    2551           long_name         = "initial specific humidity",                     &
    2552           units             = "kg/kg",                                         &
    2553           kind              = "init scalar",                                   &
    2554           input_id          = 3,                                               &
    2555           output_file       = output_file,                                     &
    2556           grid              = palm_grid,                                       &
    2557           intermediate_grid = palm_intermediate,                               &
    2558           is_profile = (TRIM(ic_mode) == 'profile')                            &
    2559        )
    2560        IF (TRIM(ic_mode) == 'profile')  THEN
    2561           output_var_table(9) % averaging_grid => averaged_initial_scalar_profile
    2562        ENDIF
    2563 
    2564        output_var_table(10) = init_nc_var(                                     &
    2565           name              = 'ls_forcing_left_qv',                            &
    2566           std_name          = "",                                              &
    2567           long_name         = "large-scale forcing for left model boundary for the specific humidity", &
    2568           units             = "kg/kg",                                         &
    2569           kind              = "left scalar",                                   &
    2570           input_id          = 3,                                               &
    2571           output_file       = output_file,                                     &
    2572           grid              = scalars_west_grid,                               &
    2573           intermediate_grid = scalars_west_intermediate                        &
    2574        )
    2575 
    2576        output_var_table(11) = init_nc_var(                                     &
    2577           name              = 'ls_forcing_right_qv',                           &
    2578           std_name          = "",                                              &
    2579           long_name         = "large-scale forcing for right model boundary for the specific humidity", &
    2580           units             = "kg/kg",                                         &
    2581           kind              = "right scalar",                                  &
    2582           input_id          = 3,                                               &
    2583           output_file       = output_file,                                     &
    2584           grid              = scalars_east_grid,                               &
    2585           intermediate_grid = scalars_east_intermediate                        &
    2586        )
    2587 
    2588        output_var_table(12) = init_nc_var(                                     &
    2589           name              = 'ls_forcing_north_qv',                           &
    2590           std_name          = "",                                              &
    2591           long_name         = "large-scale forcing for north model boundary for the specific humidity", &
    2592           units             = "kg/kg",                                         &
    2593           kind              = "north scalar",                                  &
    2594           input_id          = 3,                                               &
    2595           output_file       = output_file,                                     &
    2596           grid              = scalars_north_grid,                              &
    2597           intermediate_grid = scalars_north_intermediate                       &
    2598        )
    2599 
    2600        output_var_table(13) = init_nc_var(                                     &
    2601           name              = 'ls_forcing_south_qv',                           &
    2602           std_name          = "",                                              &
    2603           long_name         = "large-scale forcing for south model boundary for the specific humidity", &
    2604           units             = "kg/kg",                                         &
    2605           kind              = "south scalar",                                  &
    2606           input_id          = 3,                                               &
    2607           output_file       = output_file,                                     &
    2608           grid              = scalars_south_grid,                              &
    2609           intermediate_grid = scalars_south_intermediate                       &
    2610        )
    2611 
    2612        output_var_table(14) = init_nc_var(                                     &
    2613           name              = 'ls_forcing_top_qv',                             &
    2614           std_name          = "",                                              &
    2615           long_name         = "large-scale forcing for top model boundary for the specific humidity", &
    2616           units             = "kg/kg",                                         &
    2617           kind              = "top scalar",                                    &
    2618           input_id          = 3,                                               &
    2619           output_file       = output_file,                                     &
    2620           grid              = scalars_top_grid,                                &
    2621           intermediate_grid = scalars_top_intermediate                         &
    2622        )
    2623 
    2624        output_var_table(15) = init_nc_var(                                     &
    2625           name              = 'init_atmosphere_u',                             &
    2626           std_name          = "",                                              &
    2627           long_name         = "initial wind component in x direction",         &
    2628           units             = "m/s",                                           &
    2629           kind              = "init u",                                        &
    2630           input_id          = 1,                                               & ! first in (U, V) I/O group
    2631           output_file       = output_file,                                     &
    2632           grid              = u_initial_grid,                                  &
    2633           intermediate_grid = u_initial_intermediate,                          &
    2634           is_profile = (TRIM(ic_mode) == 'profile')                            &
    2635        )
    2636        IF (TRIM(ic_mode) == 'profile')  THEN
    2637           output_var_table(15) % averaging_grid => averaged_initial_scalar_profile
    2638        ENDIF
    2639 
    2640        output_var_table(16) = init_nc_var(                                     &
    2641           name              = 'ls_forcing_left_u',                             &
    2642           std_name          = "",                                              &
    2643           long_name         = "large-scale forcing for left model boundary for the wind component in x direction", &
    2644           units             = "m/s",                                           &
    2645           kind              = "left u",                                        &
    2646           input_id          = 1,                                               & ! first in (U, V) I/O group
    2647           output_file       = output_file,                                     &
    2648           grid              = u_west_grid,                                     &
    2649           intermediate_grid = u_west_intermediate                              &
    2650        )
    2651 
    2652        output_var_table(17) = init_nc_var(                                     &
    2653           name              = 'ls_forcing_right_u',                            &
    2654           std_name          = "",                                              &
    2655           long_name         = "large-scale forcing for right model boundary for the wind component in x direction", &
    2656           units             = "m/s",                                           &
    2657           kind              = "right u",                                       &
    2658           input_id          = 1,                                               & ! first in (U, V) I/O group
    2659           output_file       = output_file,                                     &
    2660           grid              = u_east_grid,                                     &
    2661           intermediate_grid = u_east_intermediate                              &
    2662        )
    2663 
    2664        output_var_table(18) = init_nc_var(                                     &
    2665           name              = 'ls_forcing_north_u',                            &
    2666           std_name          = "",                                              &
    2667           long_name         = "large-scale forcing for north model boundary for the wind component in x direction", &
    2668           units             = "m/s",                                           &
    2669           kind              = "north u",                                       &
    2670           input_id          = 1,                                               & ! first in (U, V) I/O group
    2671           output_file       = output_file,                                     &
    2672           grid              = u_north_grid,                                    &
    2673           intermediate_grid = u_north_intermediate                             &
    2674        )
    2675 
    2676        output_var_table(19) = init_nc_var(                                     &
    2677           name              = 'ls_forcing_south_u',                            &
    2678           std_name          = "",                                              &
    2679           long_name         = "large-scale forcing for south model boundary for the wind component in x direction", &
    2680           units             = "m/s",                                           &
    2681           kind              = "south u",                                       &
    2682           input_id          = 1,                                               & ! first in (U, V) I/O group
    2683           output_file       = output_file,                                     &
    2684           grid              = u_south_grid,                                    &
    2685           intermediate_grid = u_south_intermediate                             &
    2686        )
    2687 
    2688        output_var_table(20) = init_nc_var(                                     &
    2689           name              = 'ls_forcing_top_u',                              &
    2690           std_name          = "",                                              &
    2691           long_name         = "large-scale forcing for top model boundary for the wind component in x direction", &
    2692           units             = "m/s",                                           &
    2693           kind              = "top u",                                         &
    2694           input_id          = 1,                                               & ! first in (U, V) I/O group
    2695           output_file       = output_file,                                     &
    2696           grid              = u_top_grid,                                      &
    2697           intermediate_grid = u_top_intermediate                               &
    2698        )
    2699 
    2700        output_var_table(21) = init_nc_var(                                     &
    2701           name              = 'init_atmosphere_v',                             &
    2702           std_name          = "",                                              &
    2703           long_name         = "initial wind component in y direction",         &
    2704           units             = "m/s",                                           &
    2705           kind              = "init v",                                        &
    2706           input_id          = 2,                                               & ! second in (U, V) I/O group
    2707           output_file       = output_file,                                     &
    2708           grid              = v_initial_grid,                                  &
    2709           intermediate_grid = v_initial_intermediate,                          &
    2710           is_profile = (TRIM(ic_mode) == 'profile')                            &
    2711        )
    2712        IF (TRIM(ic_mode) == 'profile')  THEN
    2713           output_var_table(21) % averaging_grid => averaged_initial_scalar_profile
    2714        ENDIF
    2715 
    2716        output_var_table(22) = init_nc_var(                                     &
    2717           name              = 'ls_forcing_left_v',                             &
    2718           std_name          = "",                                              &
    2719           long_name         = "large-scale forcing for left model boundary for the wind component in y direction", &
    2720           units             = "m/s",                                           &
    2721           kind              = "right v",                                       &
    2722           input_id          = 2,                                               & ! second in (U, V) I/O group
    2723           output_file       = output_file,                                     &
    2724           grid              = v_west_grid,                                     &
    2725           intermediate_grid = v_west_intermediate                              &
    2726        )
    2727 
    2728        output_var_table(23) = init_nc_var(                                     &
    2729           name              = 'ls_forcing_right_v',                            &
    2730           std_name          = "",                                              &
    2731           long_name         = "large-scale forcing for right model boundary for the wind component in y direction", &
    2732           units             = "m/s",                                           &
    2733           kind              = "right v",                                       &
    2734           input_id          = 2,                                               & ! second in (U, V) I/O group
    2735           output_file       = output_file,                                     &
    2736           grid              = v_east_grid,                                     &
    2737           intermediate_grid = v_east_intermediate                              &
    2738        )
    2739 
    2740        output_var_table(24) = init_nc_var(                                     &
    2741           name              = 'ls_forcing_north_v',                            &
    2742           std_name          = "",                                              &
    2743           long_name         = "large-scale forcing for north model boundary for the wind component in y direction", &
    2744           units             = "m/s",                                           &
    2745           kind              = "north v",                                       &
    2746           input_id          = 2,                                               & ! second in (U, V) I/O group
    2747           output_file       = output_file,                                     &
    2748           grid              = v_north_grid,                                    &
    2749           intermediate_grid = v_north_intermediate                             &
    2750        )
    2751 
    2752        output_var_table(25) = init_nc_var(                                     &
    2753           name              = 'ls_forcing_south_v',                            &
    2754           std_name          = "",                                              &
    2755           long_name         = "large-scale forcing for south model boundary for the wind component in y direction", &
    2756           units             = "m/s",                                           &
    2757           kind              = "south v",                                       &
    2758           input_id          = 2,                                               & ! second in (U, V) I/O group
    2759           output_file       = output_file,                                     &
    2760           grid              = v_south_grid,                                    &
    2761           intermediate_grid = v_south_intermediate                             &
    2762        )
    2763 
    2764        output_var_table(26) = init_nc_var(                                     &
    2765           name              = 'ls_forcing_top_v',                              &
    2766           std_name          = "",                                              &
    2767           long_name         = "large-scale forcing for top model boundary for the wind component in y direction", &
    2768           units             = "m/s",                                           &
    2769           kind              = "top v",                                         &
    2770           input_id          = 2,                                               & ! second in (U, V) I/O group
    2771           output_file       = output_file,                                     &
    2772           grid              = v_top_grid,                                      &
    2773           intermediate_grid = v_top_intermediate                               &
    2774        )
    2775 
    2776        output_var_table(27) = init_nc_var(                                     &
    2777           name              = 'init_atmosphere_w',                             &
    2778           std_name          = "",                                              &
    2779           long_name         = "initial wind component in z direction",         &
    2780           units             = "m/s",                                           &
    2781           kind              = "init w",                                        &
    2782           input_id          = 1,                                               &
    2783           output_file       = output_file,                                     &
    2784           grid              = w_initial_grid,                                  &
    2785           intermediate_grid = w_initial_intermediate,                          &
    2786           is_profile = (TRIM(ic_mode) == 'profile')                            &
    2787        )
    2788        IF (TRIM(ic_mode) == 'profile')  THEN
    2789           output_var_table(27) % averaging_grid => averaged_initial_w_profile
    2790        ENDIF
    2791 
    2792        output_var_table(28) = init_nc_var(                                     &
    2793           name              = 'ls_forcing_left_w',                             &
    2794           std_name          = "",                                              &
    2795           long_name         = "large-scale forcing for left model boundary for the wind component in z direction", &
    2796           units             = "m/s",                                           &
    2797           kind              = "left w",                                        &
    2798           input_id          = 1,                                               &
    2799           output_file       = output_file,                                     &
    2800           grid              = w_west_grid,                                     &
    2801           intermediate_grid = w_west_intermediate                              &
    2802        )
    2803 
    2804        output_var_table(29) = init_nc_var(                                     &
    2805           name              = 'ls_forcing_right_w',                            &
    2806           std_name          = "",                                              &
    2807           long_name         = "large-scale forcing for right model boundary for the wind component in z direction", &
    2808           units             = "m/s",                                           &
    2809           kind              = "right w",                                       &
    2810           input_id          = 1,                                               &
    2811           output_file       = output_file,                                     &
    2812           grid              = w_east_grid,                                     &
    2813           intermediate_grid = w_east_intermediate                              &
    2814        )
    2815 
    2816        output_var_table(30) = init_nc_var(                                     &
    2817           name              = 'ls_forcing_north_w',                            &
    2818           std_name          = "",                                              &
    2819           long_name         = "large-scale forcing for north model boundary for the wind component in z direction", &
    2820           units             = "m/s",                                           &
    2821           kind              = "north w",                                       &
    2822           input_id          = 1,                                               &
    2823           output_file       = output_file,                                     &
    2824           grid              = w_north_grid,                                    &
    2825           intermediate_grid = w_north_intermediate                             &
    2826        )
    2827 
    2828        output_var_table(31) = init_nc_var(                                     &
    2829           name              = 'ls_forcing_south_w',                            &
    2830           std_name          = "",                                              &
    2831           long_name         = "large-scale forcing for south model boundary for the wind component in z direction", &
    2832           units             = "m/s",                                           &
    2833           kind              = "south w",                                       &
    2834           input_id          = 1,                                               &
    2835           output_file       = output_file,                                     &
    2836           grid              = w_south_grid,                                    &
    2837           intermediate_grid = w_south_intermediate                             &
    2838        )
    2839 
    2840        output_var_table(32) = init_nc_var(                                     &
    2841           name              = 'ls_forcing_top_w',                              &
    2842           std_name          = "",                                              &
    2843           long_name         = "large-scale forcing for top model boundary for the wind component in z direction", &
    2844           units             = "m/s",                                           &
    2845           kind              = "top w",                                         &
    2846           input_id          = 1,                                               &
    2847           output_file       = output_file,                                     &
    2848           grid              = w_top_grid,                                      &
    2849           intermediate_grid = w_top_intermediate                               &
    2850        )
    2851 
    2852        output_var_table(33) = init_nc_var(                                     &
    2853           name              = 'ls_forcing_soil_rain',                          &
    2854           std_name          = "",                                              &
    2855           long_name         = "large-scale forcing rain",                      &
    2856           units             = "kg/m2",                                         &
    2857           kind              = "surface forcing",                               &
    2858           input_id          = 1,                                               &
    2859           output_file       = output_file,                                     &
    2860           grid              = palm_grid,                                       &
    2861           intermediate_grid = palm_intermediate                                &
    2862        )
    2863 
    2864        output_var_table(34) = init_nc_var(                                     &
    2865           name              = 'ls_forcing_soil_snow',                          &
    2866           std_name          = "",                                              &
    2867           long_name         = "large-scale forcing snow",                      &
    2868           units             = "kg/m2",                                         &
    2869           kind              = "surface forcing",                               &
    2870           input_id          = 1,                                               &
    2871           output_file       = output_file,                                     &
    2872           grid              = palm_grid,                                       &
    2873           intermediate_grid = palm_intermediate                                &
    2874        )
    2875 
    2876        output_var_table(35) = init_nc_var(                                     &
    2877           name              = 'ls_forcing_soil_graupel',                       &
    2878           std_name          = "",                                              &
    2879           long_name         = "large-scale forcing graupel",                   &
    2880           units             = "kg/m2",                                         &
    2881           kind              = "surface forcing",                               &
    2882           input_id          = 1,                                               &
    2883           output_file       = output_file,                                     &
    2884           grid              = palm_grid,                                       &
    2885           intermediate_grid = palm_intermediate                                &
    2886        )
    2887 
    2888        output_var_table(36) = init_nc_var(                                     &
    2889           name              = 'ls_forcing_soil_t_2m',                          &
    2890           std_name          = "",                                              &
    2891           long_name         = "large-scale forcing 2m air temperature",        &
    2892           units             = "kg/m2",                                         &
    2893           kind              = "surface forcing",                               &
    2894           input_id          = 1,                                               &
    2895           output_file       = output_file,                                     &
    2896           grid              = palm_grid,                                       &
    2897           intermediate_grid = palm_intermediate                                &
    2898        )
    2899 
    2900        output_var_table(37) = init_nc_var(                                     &
    2901           name              = 'ls_forcing_soil_evap',                          &
    2902           std_name          = "",                                              &
    2903           long_name         = "large-scale forcing evapo-transpiration",       &
    2904           units             = "kg/m2",                                         &
    2905           kind              = "surface forcing",                               &
    2906           input_id          = 1,                                               &
    2907           output_file       = output_file,                                     &
    2908           grid              = palm_grid,                                       &
    2909           intermediate_grid = palm_intermediate                                &
    2910        )
    2911 
    2912        output_var_table(38) = init_nc_var(                                     &
    2913           name              = 'rad_swd_dif_0',                                 &
    2914           std_name          = "",                                              &
    2915           long_name         = "incoming diffuse shortwave radiative flux at the surface", &
    2916           units             = "W/m2",                                          &
    2917           kind              = "surface forcing",                               &
    2918           input_id          = 1,                                               &
    2919           output_file       = output_file,                                     &
    2920           grid              = palm_grid,                                       &
    2921           intermediate_grid = palm_intermediate                                &
    2922        )
    2923 
    2924        output_var_table(39) = init_nc_var(                                     &
    2925           name              = 'rad_swd_dir_0',                                 &
    2926           std_name          = "",                                              &
    2927           long_name         = "incoming direct shortwave radiative flux at the surface", &
    2928           units             = "W/m2",                                          &
    2929           kind              = "surface forcing",                               &
    2930           input_id          = 1,                                               &
    2931           output_file       = output_file,                                     &
    2932           grid              = palm_grid,                                       &
    2933           intermediate_grid = palm_intermediate                                &
    2934        )
    2935 
    2936        output_var_table(40) = init_nc_var(                                     &
    2937           name              = 'rad_sw_bal_0',                                  &
    2938           std_name          = "",                                              &
    2939           long_name         = "shortwave radiation balance at the surface",    &
    2940           units             = "W/m2",                                          &
    2941           kind              = "surface forcing",                               &
    2942           input_id          = 1,                                               &
    2943           output_file       = output_file,                                     &
    2944           grid              = palm_grid,                                       &
    2945           intermediate_grid = palm_intermediate                                &
    2946        )
    2947 
    2948        output_var_table(41) = init_nc_var(                                     &
    2949           name              = 'rad_lw_bal_0',                                  &
    2950           std_name          = "",                                              &
    2951           long_name         = "longwave radiation balance at the surface",     &
    2952           units             = "W/m2",                                          &
    2953           kind              = "surface forcing",                               &
    2954           input_id          = 1,                                               &
    2955           output_file       = output_file,                                     &
    2956           grid              = palm_grid,                                       &
    2957           intermediate_grid = palm_intermediate                                &
    2958        )
     2467    output_var_table(1) = init_nc_var(                                      &
     2468       name              = 'init_soil_t',                                   &
     2469       std_name          = "",                                              &
     2470       long_name         = "initial soil temperature",                      &
     2471       units             = "K",                                             &
     2472       kind              = "init soil",                                     &
     2473       input_id          = 1,                                               &
     2474       output_file       = output_file,                                     &
     2475       grid              = palm_grid,                                       &
     2476       intermediate_grid = palm_intermediate                                &
     2477    )
     2478
     2479    output_var_table(2) = init_nc_var(                                      &
     2480       name              = 'init_soil_m',                                   &
     2481       std_name          = "",                                              &
     2482       long_name         = "initial soil moisture",                         &
     2483       units             = "m^3/m^3",                                       &
     2484       kind              = "init soil",                                     &
     2485       input_id          = 1,                                               &
     2486       output_file       = output_file,                                     &
     2487       grid              = palm_grid,                                       &
     2488       intermediate_grid = palm_intermediate                                &
     2489    )
     2490
     2491    output_var_table(3) = init_nc_var(                                      &
     2492       name              = 'init_atmosphere_pt',                            &
     2493       std_name          = "",                                              &
     2494       long_name         = "initial potential temperature",                 &
     2495       units             = "K",                                             &
     2496       kind              = "init scalar",                                   &
     2497       input_id          = 1,                                               & ! first in (T, p) IO group
     2498       output_file       = output_file,                                     &
     2499       grid              = palm_grid,                                       &
     2500       intermediate_grid = palm_intermediate,                               &
     2501       is_profile = (TRIM(ic_mode) == 'profile')                            &
     2502    )
     2503    IF (TRIM(ic_mode) == 'profile')  THEN
     2504       output_var_table(3)%averaging_grid => averaged_initial_scalar_profile
     2505    ENDIF
     2506
     2507    output_var_table(4) = init_nc_var(                                      &
     2508       name              = 'ls_forcing_left_pt',                            &
     2509       std_name          = "",                                              &
     2510       long_name         = "large-scale forcing for left model boundary for the potential temperature", &
     2511       units             = "K",                                             &
     2512       kind              = "left scalar",                                   &
     2513       input_id          = 1,                                               &
     2514       grid              = scalars_west_grid,                               &
     2515       intermediate_grid = scalars_west_intermediate,                       &
     2516       output_file = output_file                                            &
     2517    )
     2518
     2519    output_var_table(5) = init_nc_var(                                      &
     2520       name              = 'ls_forcing_right_pt',                           &
     2521       std_name          = "",                                              &
     2522       long_name         = "large-scale forcing for right model boundary for the potential temperature", &
     2523       units             = "K",                                             &
     2524       kind              = "right scalar",                                  &
     2525       input_id          = 1,                                               &
     2526       grid              = scalars_east_grid,                               &
     2527       intermediate_grid = scalars_east_intermediate,                       &
     2528       output_file = output_file                                            &
     2529    )
     2530
     2531    output_var_table(6) = init_nc_var(                                      &
     2532       name              = 'ls_forcing_north_pt',                           &
     2533       std_name          = "",                                              &
     2534       long_name         = "large-scale forcing for north model boundary for the potential temperature", &
     2535       units             = "K",                                             &
     2536       kind              = "north scalar",                                  &
     2537       input_id          = 1,                                               &
     2538       grid              = scalars_north_grid,                              &
     2539       intermediate_grid = scalars_north_intermediate,                      &
     2540       output_file = output_file                                            &
     2541    )
     2542
     2543    output_var_table(7) = init_nc_var(                                      &
     2544       name              = 'ls_forcing_south_pt',                           &
     2545       std_name          = "",                                              &
     2546       long_name         = "large-scale forcing for south model boundary for the potential temperature", &
     2547       units             = "K",                                             &
     2548       kind              = "south scalar",                                  &
     2549       input_id          = 1,                                               &
     2550       grid              = scalars_south_grid,                              &
     2551       intermediate_grid = scalars_south_intermediate,                      &
     2552       output_file = output_file                                            &
     2553    )
     2554
     2555    output_var_table(8) = init_nc_var(                                      &
     2556       name              = 'ls_forcing_top_pt',                             &
     2557       std_name          = "",                                              &
     2558       long_name         = "large-scale forcing for top model boundary for the potential temperature", &
     2559       units             = "K",                                             &
     2560       kind              = "top scalar",                                    &
     2561       input_id          = 1,                                               &
     2562       grid              = scalars_top_grid,                                &
     2563       intermediate_grid = scalars_top_intermediate,                        &
     2564       output_file = output_file                                            &
     2565    )
     2566
     2567    output_var_table(9) = init_nc_var(                                      &
     2568       name              = 'init_atmosphere_qv',                            &
     2569       std_name          = "",                                              &
     2570       long_name         = "initial specific humidity",                     &
     2571       units             = "kg/kg",                                         &
     2572       kind              = "init scalar",                                   &
     2573       input_id          = 3,                                               &
     2574       output_file       = output_file,                                     &
     2575       grid              = palm_grid,                                       &
     2576       intermediate_grid = palm_intermediate,                               &
     2577       is_profile = (TRIM(ic_mode) == 'profile')                            &
     2578    )
     2579    IF (TRIM(ic_mode) == 'profile')  THEN
     2580       output_var_table(9)%averaging_grid => averaged_initial_scalar_profile
     2581    ENDIF
     2582
     2583    output_var_table(10) = init_nc_var(                                     &
     2584       name              = 'ls_forcing_left_qv',                            &
     2585       std_name          = "",                                              &
     2586       long_name         = "large-scale forcing for left model boundary for the specific humidity", &
     2587       units             = "kg/kg",                                         &
     2588       kind              = "left scalar",                                   &
     2589       input_id          = 3,                                               &
     2590       output_file       = output_file,                                     &
     2591       grid              = scalars_west_grid,                               &
     2592       intermediate_grid = scalars_west_intermediate                        &
     2593    )
     2594
     2595    output_var_table(11) = init_nc_var(                                     &
     2596       name              = 'ls_forcing_right_qv',                           &
     2597       std_name          = "",                                              &
     2598       long_name         = "large-scale forcing for right model boundary for the specific humidity", &
     2599       units             = "kg/kg",                                         &
     2600       kind              = "right scalar",                                  &
     2601       input_id          = 3,                                               &
     2602       output_file       = output_file,                                     &
     2603       grid              = scalars_east_grid,                               &
     2604       intermediate_grid = scalars_east_intermediate                        &
     2605    )
     2606
     2607    output_var_table(12) = init_nc_var(                                     &
     2608       name              = 'ls_forcing_north_qv',                           &
     2609       std_name          = "",                                              &
     2610       long_name         = "large-scale forcing for north model boundary for the specific humidity", &
     2611       units             = "kg/kg",                                         &
     2612       kind              = "north scalar",                                  &
     2613       input_id          = 3,                                               &
     2614       output_file       = output_file,                                     &
     2615       grid              = scalars_north_grid,                              &
     2616       intermediate_grid = scalars_north_intermediate                       &
     2617    )
     2618
     2619    output_var_table(13) = init_nc_var(                                     &
     2620       name              = 'ls_forcing_south_qv',                           &
     2621       std_name          = "",                                              &
     2622       long_name         = "large-scale forcing for south model boundary for the specific humidity", &
     2623       units             = "kg/kg",                                         &
     2624       kind              = "south scalar",                                  &
     2625       input_id          = 3,                                               &
     2626       output_file       = output_file,                                     &
     2627       grid              = scalars_south_grid,                              &
     2628       intermediate_grid = scalars_south_intermediate                       &
     2629    )
     2630
     2631    output_var_table(14) = init_nc_var(                                     &
     2632       name              = 'ls_forcing_top_qv',                             &
     2633       std_name          = "",                                              &
     2634       long_name         = "large-scale forcing for top model boundary for the specific humidity", &
     2635       units             = "kg/kg",                                         &
     2636       kind              = "top scalar",                                    &
     2637       input_id          = 3,                                               &
     2638       output_file       = output_file,                                     &
     2639       grid              = scalars_top_grid,                                &
     2640       intermediate_grid = scalars_top_intermediate                         &
     2641    )
     2642
     2643    output_var_table(15) = init_nc_var(                                     &
     2644       name              = 'init_atmosphere_u',                             &
     2645       std_name          = "",                                              &
     2646       long_name         = "initial wind component in x direction",         &
     2647       units             = "m/s",                                           &
     2648       kind              = "init u",                                        &
     2649       input_id          = 1,                                               & ! first in (U, V) I/O group
     2650       output_file       = output_file,                                     &
     2651       grid              = u_initial_grid,                                  &
     2652       intermediate_grid = u_initial_intermediate,                          &
     2653       is_profile = (TRIM(ic_mode) == 'profile')                            &
     2654    )
     2655    IF (TRIM(ic_mode) == 'profile')  THEN
     2656       output_var_table(15)%averaging_grid => averaged_initial_scalar_profile
     2657    ENDIF
     2658
     2659    output_var_table(16) = init_nc_var(                                     &
     2660       name              = 'ls_forcing_left_u',                             &
     2661       std_name          = "",                                              &
     2662       long_name         = "large-scale forcing for left model boundary for the wind component in x direction", &
     2663       units             = "m/s",                                           &
     2664       kind              = "left u",                                        &
     2665       input_id          = 1,                                               & ! first in (U, V) I/O group
     2666       output_file       = output_file,                                     &
     2667       grid              = u_west_grid,                                     &
     2668       intermediate_grid = u_west_intermediate                              &
     2669    )
     2670
     2671    output_var_table(17) = init_nc_var(                                     &
     2672       name              = 'ls_forcing_right_u',                            &
     2673       std_name          = "",                                              &
     2674       long_name         = "large-scale forcing for right model boundary for the wind component in x direction", &
     2675       units             = "m/s",                                           &
     2676       kind              = "right u",                                       &
     2677       input_id          = 1,                                               & ! first in (U, V) I/O group
     2678       output_file       = output_file,                                     &
     2679       grid              = u_east_grid,                                     &
     2680       intermediate_grid = u_east_intermediate                              &
     2681    )
     2682
     2683    output_var_table(18) = init_nc_var(                                     &
     2684       name              = 'ls_forcing_north_u',                            &
     2685       std_name          = "",                                              &
     2686       long_name         = "large-scale forcing for north model boundary for the wind component in x direction", &
     2687       units             = "m/s",                                           &
     2688       kind              = "north u",                                       &
     2689       input_id          = 1,                                               & ! first in (U, V) I/O group
     2690       output_file       = output_file,                                     &
     2691       grid              = u_north_grid,                                    &
     2692       intermediate_grid = u_north_intermediate                             &
     2693    )
     2694
     2695    output_var_table(19) = init_nc_var(                                     &
     2696       name              = 'ls_forcing_south_u',                            &
     2697       std_name          = "",                                              &
     2698       long_name         = "large-scale forcing for south model boundary for the wind component in x direction", &
     2699       units             = "m/s",                                           &
     2700       kind              = "south u",                                       &
     2701       input_id          = 1,                                               & ! first in (U, V) I/O group
     2702       output_file       = output_file,                                     &
     2703       grid              = u_south_grid,                                    &
     2704       intermediate_grid = u_south_intermediate                             &
     2705    )
     2706
     2707    output_var_table(20) = init_nc_var(                                     &
     2708       name              = 'ls_forcing_top_u',                              &
     2709       std_name          = "",                                              &
     2710       long_name         = "large-scale forcing for top model boundary for the wind component in x direction", &
     2711       units             = "m/s",                                           &
     2712       kind              = "top u",                                         &
     2713       input_id          = 1,                                               & ! first in (U, V) I/O group
     2714       output_file       = output_file,                                     &
     2715       grid              = u_top_grid,                                      &
     2716       intermediate_grid = u_top_intermediate                               &
     2717    )
     2718
     2719    output_var_table(21) = init_nc_var(                                     &
     2720       name              = 'init_atmosphere_v',                             &
     2721       std_name          = "",                                              &
     2722       long_name         = "initial wind component in y direction",         &
     2723       units             = "m/s",                                           &
     2724       kind              = "init v",                                        &
     2725       input_id          = 2,                                               & ! second in (U, V) I/O group
     2726&nb