Index: palm/trunk/UTIL/inifor/src/grid.f90
===================================================================
--- palm/trunk/UTIL/inifor/src/grid.f90	(revision 3183)
+++ palm/trunk/UTIL/inifor/src/grid.f90	(revision 3395)
@@ -26,4 +26,14 @@
 ! -----------------
 ! $Id$
+! Added computation of geostrophic winds form COSMO pressure fields
+! Introduced averaging grids and internal 'output' variables for computation of
+!     geostrophic and large-scale forcings
+! Merged qv, uv and vg output variables and 'temperature' IO group into new
+!     'thermodynamics' IP group
+! Cleaned up terminal output, show some messages only with --debug option
+! Improved variable naming and readability
+! 
+! 
+! 3183 2018-07-27 14:25:55Z suehring
 ! Introduced new PALM grid stretching
 ! Updated variable names and metadata for PIDS v1.9 compatibility
@@ -61,5 +71,6 @@
         ONLY:  DATE, EARTH_RADIUS, TO_RADIANS, TO_DEGREES, PI, dp, hp, sp,     &
                SNAME, LNAME, PATH, FORCING_STEP, WATER_ID, FILL_ITERATIONS,    &
-               BETA, P_SL, T_SL, BETA, RD, G, P_REF, RD_PALM, CP_PALM, RHO_L
+               BETA, P_SL, T_SL, BETA, RD, RV, G, P_REF, RD_PALM, CP_PALM,     &
+               RHO_L, OMEGA, HECTO
     USE io,                                                                    &
         ONLY:  get_netcdf_variable, get_netcdf_attribute,                      &
@@ -68,9 +79,11 @@
         ONLY:  NF90_MAX_NAME, NF90_MAX_VAR_DIMS
     USE transform,                                                             &
-        ONLY:  rotate_to_cosmo, find_horizontal_neighbours,                    &
+        ONLY:  average_2d, rotate_to_cosmo, find_horizontal_neighbours,&
                compute_horizontal_interp_weights,                              &
-               find_vertical_neighbours_and_weights, interpolate_2d,           &
-               gamma_from_hemisphere, phic_to_phin, lamc_to_lamn, average_2d,  &
-               project, centre_velocities, phi2phirot, rla2rlarot, uv2uvrot
+               find_vertical_neighbours_and_weights_interp,                    &
+               find_vertical_neighbours_and_weights_average, interpolate_2d,   &
+               gamma_from_hemisphere, phic_to_phin, lamc_to_lamn,              &
+               project, centre_velocities, phi2phirot, rla2rlarot, uv2uvrot,   &
+               phirot2phi, rlarot2rla
     USE types
     USE util
@@ -80,4 +93,7 @@
     SAVE
     
+    REAL(dp) ::  averaging_angle   = 0.0_dp       !< latitudal and longitudal width of averaging regions [rad]
+    REAL(dp) ::  averaging_width_ns = 0.0_dp       !< longitudal width of averaging regions [m]
+    REAL(dp) ::  averaging_width_ew = 0.0_dp       !< latitudal width of averaging regions [m]
     REAL(dp) ::  phi_equat         = 0.0_dp       !< latitude of rotated equator of COSMO-DE grid [rad]
     REAL(dp) ::  phi_n             = 0.0_dp       !< latitude of rotated pole of COSMO-DE grid [rad]
@@ -87,4 +103,10 @@
     REAL(dp) ::  phi_cn            = 0.0_dp       !< latitude of the rotated pole relative to the COSMO-DE grid [rad]
     REAL(dp) ::  lambda_cn         = 0.0_dp       !< longitude of the rotated pole relative to the COSMO-DE grid [rad]
+    REAL(dp) ::  lam_centre        = 0.0_dp       !< longitude of the PLAM domain centre in the source (COSMO rotated-pole) system [rad]
+    REAL(dp) ::  phi_centre        = 0.0_dp       !< latitude of the PLAM domain centre in the source (COSMO rotated-pole) system [rad]
+    REAL(dp) ::  lam_east          = 0.0_dp       !< longitude of the east central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad]
+    REAL(dp) ::  lam_west          = 0.0_dp       !< longitude of the west central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad]
+    REAL(dp) ::  phi_north         = 0.0_dp       !< latitude of the north central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad]
+    REAL(dp) ::  phi_south         = 0.0_dp       !< latitude of the south central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad]
     REAL(dp) ::  gam               = 0.0_dp       !< angle for working around phirot2phi/rlarot2rla bug
     REAL(dp) ::  dx                = 0.0_dp       !< PALM-4U grid spacing in x direction [m]
@@ -100,8 +122,7 @@
     REAL(dp) ::  dyi               = 0.0_dp       !< inverse PALM-4U grid spacing in y direction [m^-1]
     REAL(dp) ::  dzi               = 0.0_dp       !< inverse PALM-4U grid spacing in z direction [m^-1]
+    REAL(dp) ::  f3                = 0.0_dp       !< Coriolis parameter
     REAL(dp) ::  lx                = 0.0_dp       !< PALM-4U domain size in x direction [m]
     REAL(dp) ::  ly                = 0.0_dp       !< PALM-4U domain size in y direction [m]
-    REAL(dp) ::  ug                = 0.0_dp       !< geostrophic wind in x direction [m/s]
-    REAL(dp) ::  vg                = 0.0_dp       !< geostrophic wind in y direction [m/s]
     REAL(dp) ::  p0                = 0.0_dp       !< PALM-4U surface pressure, at z0 [Pa]
     REAL(dp) ::  x0                = 0.0_dp       !< x coordinate of PALM-4U Earth tangent [m] 
@@ -140,7 +161,4 @@
     INTEGER ::  dz_stretch_level_end_index(9)               !< vertical grid level index until which the vertical grid spacing is stretched
     INTEGER ::  dz_stretch_level_start_index(9)             !< vertical grid level index above which the vertical grid spacing is stretched
-    INTEGER ::  i     !< indexing variable
-    INTEGER ::  average_imin, average_imax, average_jmin, average_jmax !< index ranges for profile averaging
-    INTEGER ::  k     !< indexing variable
     INTEGER ::  nt    !< number of output time steps
     INTEGER ::  nx    !< number of PALM-4U grid points in x direction
@@ -164,6 +182,4 @@
     LOGICAL ::  profile_grids_required
 
-    INTEGER ::  n_invar = 0 !< number of variables in the input variable table
-    INTEGER ::  n_outvar = 0 !< number of variables in the output variable table
     TYPE(nc_var), ALLOCATABLE, TARGET ::  input_var_table(:)  !< table of input variables
     TYPE(nc_var), ALLOCATABLE, TARGET ::  output_var_table(:) !< table of input variables
@@ -219,8 +235,11 @@
     TYPE(grid_definition), TARGET ::  w_south_intermediate            !< 
     TYPE(grid_definition), TARGET ::  w_top_intermediate              !< 
-    TYPE(grid_definition), TARGET ::  scalar_profile_grid             !< 
-    TYPE(grid_definition), TARGET ::  scalar_profile_intermediate     !< 
-    TYPE(grid_definition), TARGET ::  w_profile_grid                  !< 
-    TYPE(grid_definition), TARGET ::  w_profile_intermediate          !< 
+
+    TYPE(grid_definition), TARGET ::  north_averaged_scalar_profile
+    TYPE(grid_definition), TARGET ::  south_averaged_scalar_profile
+    TYPE(grid_definition), TARGET ::  west_averaged_scalar_profile
+    TYPE(grid_definition), TARGET ::  east_averaged_scalar_profile
+    TYPE(grid_definition), TARGET ::  averaged_scalar_profile
+    TYPE(grid_definition), TARGET ::  averaged_w_profile
 
     TYPE(io_group), ALLOCATABLE, TARGET ::  io_group_list(:)  !< List of I/O groups, which group together output variables that share the same input variable
@@ -237,4 +256,8 @@
     CHARACTER(LEN=PATH), ALLOCATABLE, DIMENSION(:) ::  radiation_files
     CHARACTER(LEN=SNAME) ::  input_prefix         !< Prefix of input files, e.g. 'laf' for COSMO-DE analyses
+    CHARACTER(LEN=SNAME) ::  flow_prefix          !< Prefix of flow input files, e.g. 'laf' for COSMO-DE analyses
+    CHARACTER(LEN=SNAME) ::  soil_prefix          !< Prefix of soil input files, e.g. 'laf' for COSMO-DE analyses
+    CHARACTER(LEN=SNAME) ::  radiation_prefix     !< Prefix of radiation input files, e.g. 'laf' for COSMO-DE analyses
+    CHARACTER(LEN=SNAME) ::  soilmoisture_prefix  !< Prefix of input files for soil moisture spin-up, e.g. 'laf' for COSMO-DE analyses
     CHARACTER(LEN=SNAME) ::  flow_suffix          !< Suffix of flow input files, e.g. 'flow'
     CHARACTER(LEN=SNAME) ::  soil_suffix          !< Suffix of soil input files, e.g. 'soil'
@@ -249,4 +272,5 @@
     
     SUBROUTINE setup_parameters()
+       INTEGER ::  k
 
 !
@@ -282,6 +306,6 @@
 
        ! Default atmospheric parameters
-       ug = 0.0_dp
-       vg = 0.0_dp
+       cfg % ug = 0.0_dp
+       cfg % vg = 0.0_dp
        cfg % p0 = P_SL
 
@@ -296,4 +320,5 @@
        input_prefix = 'laf'  ! 'laf' for COSMO-DE analyses
        cfg % flow_prefix = input_prefix
+       cfg % input_prefix = input_prefix
        cfg % soil_prefix = input_prefix
        cfg % radiation_prefix = input_prefix
@@ -304,4 +329,7 @@
        radiation_suffix = '-rad'
        soilmoisture_suffix = '-soilmoisture'
+
+       cfg % debug = .FALSE.
+       cfg % averaging_angle = 2.0_dp
 !
 !------------------------------------------------------------------------------
@@ -318,7 +346,17 @@
        cfg % ic_mode = 'volume'
        cfg % bc_mode = 'real'
-
-       ! Update default file names and domain centre
+       cfg % averaging_mode = 'level'
+
+       ! Overwrite defaults with user configuration
        CALL parse_command_line_arguments( cfg )
+
+       flow_prefix = TRIM(cfg % input_prefix)
+       radiation_prefix = TRIM(cfg % input_prefix)
+       soil_prefix = TRIM(cfg % input_prefix)
+       soilmoisture_prefix = TRIM(cfg % input_prefix)
+       IF (cfg % flow_prefix_is_set)  flow_prefix = TRIM(cfg % flow_prefix)
+       IF (cfg % radiation_prefix_is_set)  radiation_prefix = TRIM(cfg % radiation_prefix)
+       IF (cfg % soil_prefix_is_set)  soil_prefix = TRIM(cfg % soil_prefix)
+       IF (cfg % soilmoisture_prefix_is_set)  soilmoisture_prefix = TRIM(cfg % soilmoisture_prefix)
 
        output_file % name = cfg % output_file
@@ -336,4 +374,6 @@
        END IF
 
+
+       ! Set default file paths, if not specified by user.
        CALL normalize_path(cfg % input_path)
        IF (TRIM(cfg % hhl_file) == '')  cfg % hhl_file = TRIM(cfg % input_path) // 'hhl.nc'
@@ -361,12 +401,16 @@
 
        ! Generate input file lists
-       CALL get_input_file_list(cfg % start_date, start_hour_flow, end_hour, step_hour,  &
-          cfg % input_path, cfg % flow_prefix, flow_suffix, flow_files)
-       CALL get_input_file_list(cfg % start_date, start_hour_soil, end_hour, step_hour,  &
-          cfg % input_path, cfg % soil_prefix, soil_suffix, soil_files)
-       CALL get_input_file_list(cfg % start_date, start_hour_radiation, end_hour, step_hour, &
-          cfg % input_path, cfg % radiation_prefix, radiation_suffix, radiation_files)
-       CALL get_input_file_list(cfg % start_date, start_hour_soilmoisture, end_hour_soilmoisture, step_hour, &
-          cfg % input_path, cfg % soilmoisture_prefix, soilmoisture_suffix, soil_moisture_files)
+       CALL get_input_file_list(                                               &
+          cfg % start_date, start_hour_flow, end_hour, step_hour,              &
+          cfg % input_path, flow_prefix, flow_suffix, flow_files)
+       CALL get_input_file_list(                                               &
+          cfg % start_date, start_hour_soil, end_hour, step_hour,              &
+          cfg % input_path, soil_prefix, soil_suffix, soil_files)
+       CALL get_input_file_list(                                               &
+          cfg % start_date, start_hour_radiation, end_hour, step_hour,         &
+          cfg % input_path, radiation_prefix, radiation_suffix, radiation_files)
+       CALL get_input_file_list(                                               &
+          cfg % start_date, start_hour_soilmoisture, end_hour_soilmoisture, step_hour, &
+          cfg % input_path, soilmoisture_prefix, soilmoisture_suffix, soil_moisture_files)
 
 !
@@ -499,4 +543,47 @@
 
  CALL run_control('time', 'comp')
+
+!------------------------------------------------------------------------------
+! Section 4.3: INIFOR averaging domains
+!------------------------------------------------------------------------------
+
+    ! Compute coordiantes of the PALM centre in the source (COSMO) system
+    phi_centre = phirot2phi(                                                   &
+       phirot = project(0.5_dp*ly, y0, EARTH_RADIUS) * TO_DEGREES,             &
+       rlarot = project(0.5_dp*lx, x0, EARTH_RADIUS) * TO_DEGREES,             &
+       polphi = phi_cn * TO_DEGREES, pollam = lambda_cn * TO_DEGREES,          &
+       polgam = gam * TO_DEGREES                                               &
+    ) * TO_RADIANS
+
+    lam_centre = rlarot2rla(                                                   &
+       phirot = project(0.5_dp*ly, y0, EARTH_RADIUS) * TO_DEGREES,             &
+       rlarot = project(0.5_dp*lx, x0, EARTH_RADIUS) * TO_DEGREES,             &
+       polphi = phi_cn * TO_DEGREES, pollam = lambda_cn * TO_DEGREES,          &
+       polgam = gam * TO_DEGREES                                               &
+    ) * TO_RADIANS
+
+    message = "PALM-4U centre:" // NEW_LINE('') // &
+       "           lon (lambda) = " // &
+       TRIM(real_to_str_f(lam_centre * TO_DEGREES)) // " deg" // NEW_LINE(' ') // &
+       "           lat (phi   ) = " // &
+       TRIM(real_to_str_f(phi_centre * TO_DEGREES)) // " deg (COSMO-DE rotated-pole)"
+    CALL report( 'setup_parameters', message )
+
+
+    ! Compute boundaries of the central averaging box
+    averaging_angle = cfg % averaging_angle * TO_RADIANS
+    lam_east = lam_centre + 0.5_dp * averaging_angle
+    lam_west = lam_centre - 0.5_dp * averaging_angle
+    phi_north = phi_centre + 0.5_dp * averaging_angle
+    phi_south = phi_centre - 0.5_dp * averaging_angle
+    averaging_width_ew = averaging_angle * COS(phi_centre) * EARTH_RADIUS
+    averaging_width_ns = averaging_angle * EARTH_RADIUS
+
+    ! Coriolis parameter
+    f3 = 2.0_dp * OMEGA * SIN(                                                 &
+       TO_RADIANS*phirot2phi( phi_centre * TO_DEGREES, lam_centre * TO_DEGREES,&
+                              phi_n * TO_DEGREES, lambda_n * TO_DEGREES,       &
+                              gam * TO_DEGREES )                               &
+       )
 
     END SUBROUTINE setup_parameters
@@ -603,5 +690,5 @@
                ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy,                    &
                x0=x0, y0=y0, z0 = z0,                                          &
-               nx = nx, ny = ny, nz = nlev - 2)
+               nx = nx, ny = ny, nz = nlev - 1)
 
       IF (boundary_variables_required)  THEN
@@ -824,5 +911,5 @@
                   ymin =  0.5_dp * dy, ymax = ly - 0.5_dp * dy,                   &
                   x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = 0, ny = ny, nz = nlev - 2)
+                  nx = 0, ny = ny, nz = nlev - 1)
 
           CALL init_grid_definition('boundary intermediate', grid=w_west_intermediate,         &
@@ -830,5 +917,5 @@
                   ymin =  0.5_dp * dy, ymax = ly - 0.5_dp * dy,                   &
                   x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = 0, ny = ny, nz = nlev - 2)
+                  nx = 0, ny = ny, nz = nlev - 1)
 
           CALL init_grid_definition('boundary intermediate', grid=w_north_intermediate,        &
@@ -836,5 +923,5 @@
                   ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy,               &
                   x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = nx, ny = 0, nz = nlev - 2)
+                  nx = nx, ny = 0, nz = nlev - 1)
 
           CALL init_grid_definition('boundary intermediate', grid=w_south_intermediate,        &
@@ -842,5 +929,5 @@
                   ymin = -0.5_dp * dy, ymax = -0.5_dp * dy,                       &
                   x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = nx, ny = 0, nz = nlev - 2)
+                  nx = nx, ny = 0, nz = nlev - 1)
 
           CALL init_grid_definition('boundary intermediate', grid=w_top_intermediate,          &
@@ -848,5 +935,5 @@
                   ymin =  0.5_dp * dy, ymax = ly - 0.5_dp * dy,                   &
                   x0=x0, y0=y0, z0 = z0,                                          &
-                  nx = nx, ny = ny, nz = nlev - 2)
+                  nx = nx, ny = ny, nz = nlev - 1)
        END IF
 
@@ -856,32 +943,74 @@
 !------------------------------------------------------------------------------
 
-       profile_grids_required = ( TRIM(cfg % ic_mode) == 'profile' .OR.              &
+       CALL init_averaging_grid(averaged_scalar_profile, cosmo_grid,           &
+               x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0,               &
+               lonmin = lam_west, lonmax = lam_east,                           &
+               latmin = phi_south, latmax = phi_north,                         &
+               kind='scalar')
+
+       CALL init_averaging_grid(averaged_w_profile, cosmo_grid,                &
+               x = 0.5_dp * lx, y = 0.5_dp * ly, z = zw, z0 = z0,              &
+               lonmin = lam_west, lonmax = lam_east,                           &
+               latmin = phi_south, latmax = phi_north,                         &
+               kind='w')
+
+       CALL init_averaging_grid(south_averaged_scalar_profile, cosmo_grid,     &
+               x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0,               &
+               lonmin = lam_west, lonmax = lam_east,                           &
+               latmin = phi_centre - averaging_angle,                          &
+               latmax = phi_centre,                                            &
+               kind='scalar')
+
+       CALL init_averaging_grid(north_averaged_scalar_profile, cosmo_grid,     &
+               x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0,               &
+               lonmin = lam_west, lonmax = lam_east,                           &
+               latmin = phi_centre,                                            &
+               latmax = phi_centre + averaging_angle,                          &
+               kind='scalar')
+
+       CALL init_averaging_grid(west_averaged_scalar_profile, cosmo_grid,      &
+               x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0,               &
+               lonmin = lam_centre - averaging_angle,                          &
+               lonmax = lam_centre,                                            &
+               latmin = phi_south, latmax = phi_north,                         &
+               kind='scalar')
+
+       CALL init_averaging_grid(east_averaged_scalar_profile, cosmo_grid,      &
+               x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0,               &
+               lonmin = lam_centre,                                            &
+               lonmax = lam_centre + averaging_angle,                          &
+               latmin = phi_south, latmax = phi_north,                         &
+               kind='scalar')
+
+       profile_grids_required = ( TRIM(cfg % ic_mode) == 'profile' .OR.        &
                                   TRIM(cfg % bc_mode) == 'ideal' )
 
-       IF (profile_grids_required)  THEN
-          CALL init_grid_definition('boundary', grid=scalar_profile_grid,      &
-                  xmin = 0.5_dp * lx, xmax = 0.5_dp * lx,                      &
-                  ymin = 0.5_dp * ly, ymax = 0.5_dp * ly,                      &
-                  x0=x0, y0=y0, z0 = z0,                                       &
-                  nx = 0, ny = 0, nz = nz, z=z)
-
-          CALL init_grid_definition('boundary', grid=w_profile_grid,           &
-                  xmin = 0.5_dp * lx, xmax = 0.5_dp * lx,                      &
-                  ymin = 0.5_dp * ly, ymax = 0.5_dp * ly,                      &
-                  x0=x0, y0=y0, z0 = z0,                                       &
-                  nx = 0, ny = 0, nz = nz - 1, z=zw)
-
-          CALL init_grid_definition('boundary', grid=scalar_profile_intermediate,&
-                  xmin = 0.5_dp * lx, xmax = 0.5_dp * lx,                      &
-                  ymin = 0.5_dp * ly, ymax = 0.5_dp * ly,                      &
-                  x0=x0, y0=y0, z0 = z0,                                       &
-                  nx = 0, ny = 0, nz = nlev - 2, z=z)
-
-          CALL init_grid_definition('boundary', grid=w_profile_intermediate,   &
-                  xmin = 0.5_dp * lx, xmax = 0.5_dp * lx,                      &
-                  ymin = 0.5_dp * ly, ymax = 0.5_dp * ly,                      &
-                  x0=x0, y0=y0, z0 = z0,                                       &
-                  nx = 0, ny = 0, nz = nlev - 2, z=zw)
-       END IF
+       !IF (profile_grids_required)  THEN
+       !   ! Here, I use the 'boundary' kind, so the vertical coordinate uses the
+       !   ! PALM z coordinate.
+       !   CALL init_grid_definition('boundary', grid=scalar_profile_grid,      &
+       !           xmin = 0.5_dp * lx, xmax = 0.5_dp * lx,                      &
+       !           ymin = 0.5_dp * ly, ymax = 0.5_dp * ly,                      &
+       !           x0=x0, y0=y0, z0 = z0,                                       &
+       !           nx = 0, ny = 0, nz = nz, z=z)
+
+       !   CALL init_grid_definition('boundary', grid=w_profile_grid,           &
+       !           xmin = 0.5_dp * lx, xmax = 0.5_dp * lx,                      &
+       !           ymin = 0.5_dp * ly, ymax = 0.5_dp * ly,                      &
+       !           x0=x0, y0=y0, z0 = z0,                                       &
+       !           nx = 0, ny = 0, nz = nz - 1, z=zw)
+
+       !   CALL init_grid_definition('boundary', grid=scalar_profile_intermediate,&
+       !           xmin = 0.5_dp * lx, xmax = 0.5_dp * lx,                      &
+       !           ymin = 0.5_dp * ly, ymax = 0.5_dp * ly,                      &
+       !           x0=x0, y0=y0, z0 = z0,                                       &
+       !           nx = 0, ny = 0, nz = nlev - 2, z=z)
+
+       !   CALL init_grid_definition('boundary', grid=w_profile_intermediate,   &
+       !           xmin = 0.5_dp * lx, xmax = 0.5_dp * lx,                      &
+       !           ymin = 0.5_dp * ly, ymax = 0.5_dp * ly,                      &
+       !           x0=x0, y0=y0, z0 = z0,                                       &
+       !           nx = 0, ny = 0, nz = nlev - 2, z=zw)
+       !END IF
 
 !                                                                              
@@ -930,6 +1059,5 @@
 
        IF (TRIM(cfg % ic_mode) == 'profile')  THEN
-           CALL setup_averaging(cosmo_grid, palm_intermediate,                 &
-                average_imin, average_imax, average_jmin, average_jmax) 
+           !TODO: remove this conditional if not needed. 
        END IF
         
@@ -938,21 +1066,15 @@
 
 
-    SUBROUTINE setup_interpolation(cosmo_grid, palm_grid, palm_intermediate, kind, ic_mode)
+    SUBROUTINE setup_interpolation(cosmo_grid, grid, intermediate_grid, kind, ic_mode)
 
        TYPE(grid_definition), INTENT(IN), TARGET    ::  cosmo_grid
-       TYPE(grid_definition), INTENT(INOUT), TARGET ::  palm_grid, palm_intermediate
+       TYPE(grid_definition), INTENT(INOUT), TARGET ::  grid, intermediate_grid
        CHARACTER, INTENT(IN)                        ::  kind
        CHARACTER(LEN=*), INTENT(IN), OPTIONAL       ::  ic_mode
 
-       TYPE(grid_definition), POINTER      ::  grid
-       REAL(dp), DIMENSION(:), POINTER     ::  lat, lon
-       REAL(dp), DIMENSION(:,:,:), POINTER ::  h
-
-       LOGICAL :: setup_vertical
-
-       setup_vertical = .TRUE.
-       IF (PRESENT(ic_mode))  THEN
-          IF (TRIM(ic_mode) == 'profile')  setup_vertical = .FALSE.
-       END IF
+       REAL(dp), DIMENSION(:), POINTER     ::  cosmo_lat, cosmo_lon
+       REAL(dp), DIMENSION(:,:,:), POINTER ::  cosmo_h
+
+       LOGICAL :: setup_volumetric
 
 !------------------------------------------------------------------------------
@@ -964,25 +1086,25 @@
        CASE('s') ! scalars
 
-          lat => cosmo_grid % lat
-          lon => cosmo_grid % lon
-          h   => cosmo_grid % hfl
+          cosmo_lat => cosmo_grid % lat
+          cosmo_lon => cosmo_grid % lon
+          cosmo_h   => cosmo_grid % hfl
 
        CASE('w') ! vertical velocity
 
-          lat => cosmo_grid % lat
-          lon => cosmo_grid % lon
-          h   => cosmo_grid % hhl
+          cosmo_lat => cosmo_grid % lat
+          cosmo_lon => cosmo_grid % lon
+          cosmo_h   => cosmo_grid % hhl
 
        CASE('u') ! x velocity
 
-          lat => cosmo_grid % lat
-          lon => cosmo_grid % lonu
-          h   => cosmo_grid % hfl
+          cosmo_lat => cosmo_grid % lat
+          cosmo_lon => cosmo_grid % lonu
+          cosmo_h   => cosmo_grid % hfl
 
        CASE('v') ! y velocity
 
-          lat => cosmo_grid % latv
-          lon => cosmo_grid % lon
-          h   => cosmo_grid % hfl
+          cosmo_lat => cosmo_grid % latv
+          cosmo_lon => cosmo_grid % lon
+          cosmo_h   => cosmo_grid % hfl
 
        CASE DEFAULT
@@ -993,11 +1115,12 @@
        END SELECT
 
-       grid => palm_intermediate
-
-       CALL find_horizontal_neighbours(lat, lon,                               &
-          grid % clat, grid % clon, grid % ii, grid % jj)
-
-       CALL compute_horizontal_interp_weights(lat, lon,                        &
-          grid % clat, grid % clon, grid % ii, grid % jj, grid % w_horiz)
+       CALL find_horizontal_neighbours(cosmo_lat, cosmo_lon,                   &
+          intermediate_grid % clat, intermediate_grid % clon,                  &
+          intermediate_grid % ii, intermediate_grid % jj)
+
+       CALL compute_horizontal_interp_weights(cosmo_lat, cosmo_lon,            &
+          intermediate_grid % clat, intermediate_grid % clon,                  &
+          intermediate_grid % ii, intermediate_grid % jj,                      &
+          intermediate_grid % w_horiz)
 
 !------------------------------------------------------------------------------
@@ -1005,55 +1128,27 @@
 !------------------------------------------------------------------------------
 
-       IF (setup_vertical)  THEN
-          ALLOCATE( grid % h(0:grid % nx, 0:grid % ny, 0:grid % nz) ) 
-          grid % h(:,:,:) = - EARTH_RADIUS
+       ! If profile initialization is chosen, we--somewhat counterintuitively--
+       ! don't need to compute vertical interpolation weights. At least, we
+       ! don't need them on the intermediate grid, which fills the entire PALM
+       ! domain volume. Instead we need vertical weights for the intermediate
+       ! profile grids, which get computed in setup_averaging().
+       setup_volumetric = .TRUE.
+       IF (PRESENT(ic_mode))  THEN
+          IF (TRIM(ic_mode) == 'profile')  setup_volumetric = .FALSE.
+       END IF
+
+       IF (setup_volumetric)  THEN
+          ALLOCATE( intermediate_grid % h(0:intermediate_grid % nx,            &
+                                          0:intermediate_grid % ny,            &
+                                          0:intermediate_grid % nz) ) 
+          intermediate_grid % h(:,:,:) = - EARTH_RADIUS
 
           ! For w points, use hhl, for scalars use hfl
           ! compute the full heights for the intermediate grids
-          CALL interpolate_2d(cosmo_grid % hfl, grid % h, grid)
-          CALL find_vertical_neighbours_and_weights(palm_grid, grid)
+          CALL interpolate_2d(cosmo_h, intermediate_grid % h, intermediate_grid)
+          CALL find_vertical_neighbours_and_weights_interp(grid, intermediate_grid)
        END IF
        
     END SUBROUTINE setup_interpolation
-
-
-    SUBROUTINE setup_averaging(cosmo_grid, palm_intermediate, imin, imax, jmin, jmax)
-
-       TYPE(grid_definition), INTENT(IN) ::  cosmo_grid, palm_intermediate
-       INTEGER, INTENT(INOUT)            ::  imin, imax, jmin, jmax
-
-       TYPE(grid_definition), POINTER    ::  grid
-       REAL(dp)                          ::  lonmin_pos,lonmax_pos, latmin_pos, latmax_pos
-       REAL(dp)                          ::  cosmo_dxi, cosmo_dyi
-
-       cosmo_dxi = 1.0_dp / (cosmo_grid % lon(1) - cosmo_grid % lon(0))
-       cosmo_dyi = 1.0_dp / (cosmo_grid % lat(1) - cosmo_grid % lat(0))
-
-       ! find horizontal index ranges for profile averaging
-       lonmin_pos = (MINVAL(palm_intermediate % clon(:,:)) - cosmo_grid % lon(0)) * cosmo_dxi
-       lonmax_pos = (MAXVAL(palm_intermediate % clon(:,:)) - cosmo_grid % lon(0)) * cosmo_dxi
-       latmin_pos = (MINVAL(palm_intermediate % clat(:,:)) - cosmo_grid % lat(0)) * cosmo_dyi
-       latmax_pos = (MAXVAL(palm_intermediate % clat(:,:)) - cosmo_grid % lat(0)) * cosmo_dyi
-
-       imin = FLOOR(lonmin_pos)
-       imax = CEILING(lonmax_pos)
-       jmin = FLOOR(latmin_pos)
-       jmax = CEILING(latmax_pos)
-       
-       ! average heights for intermediate scalar and w profile grids
-       grid => scalar_profile_intermediate
-       ALLOCATE( grid % h(0:grid % nx, 0:grid % ny, 0:grid % nz) ) 
-       grid % h(:,:,:) = - EARTH_RADIUS
-       CALL average_2d(cosmo_grid % hfl, grid % h(0,0,:), imin, imax, jmin, jmax)
-       CALL find_vertical_neighbours_and_weights(scalar_profile_grid, grid)
-
-       grid => w_profile_intermediate
-       ALLOCATE( grid % h(0:grid % nx, 0:grid % ny, 0:grid % nz) ) 
-       grid % h(:,:,:) = - EARTH_RADIUS
-       CALL average_2d(cosmo_grid % hhl, grid % h(0,0,:), imin, imax, jmin, jmax)
-       CALL find_vertical_neighbours_and_weights(w_profile_grid, grid)
-
-    END SUBROUTINE setup_averaging
-
 
 !------------------------------------------------------------------------------!
@@ -1291,4 +1386,177 @@
 
     END SUBROUTINE init_grid_definition
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Initializes averagin_grid-type variables.
+!> 
+!> Input parameters:
+!> -----------------
+!>
+!> cosmo_grid : grid_definition-type varieble initialized as 'cosmo-de' grid
+!>    providing longitudes and latitudes of the parent grid
+!>
+!> x, y : location of the profile to be averaged in the PALM system [m]
+!>
+!> z : vertical grid coordinates of the profile in the PALM system [m]
+!>
+!> <lat/lon><min/max>: boundaries of the averaging region in the parent system,
+!>    i.e. the COSMO-DE rotated-pole system. [rad]
+!>
+!> kind : kind of quantity to be averaged using this averaging grid.
+!>    Destinguishes COSMO-DE scalar and w-velocity levels. Note that finding the
+!>    parent/COSMO columns for the region in get_cosmo_averaging_region() is
+!>    independent of 'kind' b/ecause INIFOR uses column-centred u and v velocity
+!>    components, which are computed in the preprocessing step.
+!>
+!> Output parameters:
+!> ------------------
+!> avg_grid : averaging_grid variable to be initialized.
+!------------------------------------------------------------------------------!
+    SUBROUTINE init_averaging_grid(avg_grid, cosmo_grid, x, y, z, z0,          &
+       lonmin, lonmax, latmin, latmax, kind)
+
+       TYPE(grid_definition), INTENT(INOUT) ::  avg_grid
+       TYPE(grid_definition), INTENT(IN)    ::  cosmo_grid
+       REAL(dp), INTENT(IN)                 ::  x, y, z0
+       REAL(dp), INTENT(IN), TARGET         ::  z(:)
+       REAL(dp), INTENT(IN)                 ::  lonmin, lonmax, latmin, latmax
+
+       CHARACTER(LEN=*), INTENT(IN)         ::  kind
+
+       LOGICAL                              :: level_based_averaging
+
+       INTEGER :: i, j
+
+       ALLOCATE( avg_grid % x(1) )
+       ALLOCATE( avg_grid % y(1) )
+       avg_grid % x(1) = x
+       avg_grid % y(1) = y
+       avg_grid % z => z
+       avg_grid % z0 = z0
+
+       avg_grid % nz = SIZE(z, 1)
+
+       ALLOCATE( avg_grid % lon(2) )
+       ALLOCATE( avg_grid % lat(2) )
+       avg_grid % lon(1:2) = (/lonmin, lonmax/)
+       avg_grid % lat(1:2) = (/latmin, latmax/)
+
+       avg_grid % kind = TRIM(kind)
+
+       ! Find and store COSMO columns that fall into the coordinate range
+       ! given by avg_grid % clon, % clat
+       CALL get_cosmo_averaging_region(avg_grid, cosmo_grid)
+
+       ALLOCATE (avg_grid % kkk(avg_grid % n_columns, avg_grid % nz, 2) )
+       ALLOCATE (avg_grid % w(avg_grid % n_columns, avg_grid % nz, 2) )
+       ! Compute average COSMO levels in the averaging region
+       SELECT CASE(avg_grid % kind)
+
+       CASE('scalar', 'u', 'v')
+          avg_grid % cosmo_h => cosmo_grid % hfl
+
+       CASE('w')
+          avg_grid % cosmo_h => cosmo_grid % hhl
+
+       CASE DEFAULT
+          message = "Averaging grid kind '" // TRIM(avg_grid % kind) // &
+                    "' is not supported. Use 'scalar', 'u', or 'v'."
+          CALL abort('get_cosmo_averaging_region', message)
+
+       END SELECT
+
+       ! For level-besed averaging, compute average heights
+       level_based_averaging = .TRUE.
+       IF (level_based_averaging)  THEN
+          ALLOCATE(avg_grid % h(1,1,SIZE(avg_grid % cosmo_h, 3)) )
+ 
+          CALL average_2d(avg_grid % cosmo_h, avg_grid % h(1,1,:),                            &
+                          avg_grid % iii, avg_grid % jjj)
+       END IF
+
+       ! Copy averaged grid to all COSMO columnts, leads to computing the same
+       ! vertical interpolation weights for all columns and to COSMO grid level
+       ! based averaging onto the averaged COSMO heights.
+       IF ( TRIM(cfg % averaging_mode) == 'level' )  THEN
+          FORALL(                                                              &
+             i=1:SIZE(avg_grid % cosmo_h, 1),                                  &
+             j=1:SIZE(avg_grid % cosmo_h, 2)                                   &
+          )   avg_grid % cosmo_h(i,j,:) = avg_grid % h(1,1,:)
+       END IF
+
+       ! Compute vertical weights and neighbours
+       CALL find_vertical_neighbours_and_weights_average(avg_grid)
+
+    END SUBROUTINE init_averaging_grid
+
+
+    SUBROUTINE get_cosmo_averaging_region(avg_grid, cosmo_grid)
+       TYPE(grid_definition), INTENT(INOUT)         ::  avg_grid
+       TYPE(grid_definition), TARGET, INTENT(IN)    ::  cosmo_grid
+
+       REAL(dp), DIMENSION(:), POINTER              ::  cosmo_lon, cosmo_lat
+       REAL(dp)                                     ::  dlon, dlat
+
+       INTEGER ::  i, j, imin, imax, jmin, jmax, l, nx, ny
+
+
+       SELECT CASE( TRIM(avg_grid % kind) )
+
+       CASE('scalar', 'w')
+          cosmo_lon => cosmo_grid % lon
+          cosmo_lat => cosmo_grid % lat
+
+       CASE('u')
+          cosmo_lon => cosmo_grid % lonu
+          cosmo_lat => cosmo_grid % lat
+
+       CASE('v')
+          cosmo_lon => cosmo_grid % lon
+          cosmo_lat => cosmo_grid % latv
+
+       CASE DEFAULT
+          message = "Averaging grid kind '" // TRIM(avg_grid % kind) // &
+                    "' is not supported. Use 'scalar', 'u', or 'v'."
+          CALL abort('get_cosmo_averaging_region', message)
+
+       END SELECT
+
+
+       ! FIXME: make dlon, dlat parameters of the grid_defintion type
+       dlon = cosmo_lon(1) - cosmo_lon(0)
+       dlat = cosmo_lat(1) - cosmo_lat(0)
+
+       imin = CEILING( (avg_grid % lon(1) - cosmo_lon(0)) / dlon )
+       imax = FLOOR(   (avg_grid % lon(2) - cosmo_lon(0)) / dlon )
+
+       jmin = CEILING( (avg_grid % lat(1) - cosmo_lat(0)) / dlat )
+       jmax = FLOOR(   (avg_grid % lat(2) - cosmo_lat(0)) / dlat )
+       
+       message = "Averaging over "//                                           &
+                 TRIM(str(imin)) // " < i < " // TRIM(str(imax)) //            &
+                 " and " //                                                    &
+                 TRIM(str(jmin)) // " < j < " // TRIM(str(jmax))
+       CALL report( 'get_cosmo_averaging_region', message )
+
+       nx = imax - imin + 1
+       ny = jmax - jmin + 1
+       avg_grid % n_columns = nx * ny
+
+       ALLOCATE( avg_grid % iii(avg_grid % n_columns),                                 &
+                 avg_grid % jjj(avg_grid % n_columns) )
+
+       l = 0
+       DO j = jmin, jmax
+       DO i = imin, imax
+          l = l + 1
+          avg_grid % iii(l) = i
+          avg_grid % jjj(l) = j
+       END DO
+       END DO
+
+    END SUBROUTINE get_cosmo_averaging_region
 
 
@@ -1718,25 +1986,30 @@
        )
 
-       !potential temperature, surface pressure, including nudging and subsidence
+       !potential temperature, surface pressure, specific humidity including
+       !nudging and subsidence, and geostrophic winds ug, vg
        io_group_list(3) = init_io_group(                                       &
           in_files = flow_files,                                               &
-          out_vars = [output_var_table(3:8), output_var_table(42:42),          &
-                      output_var_table(49:51)],                                &
-          in_var_list = (/input_var_table(3), input_var_table(17)/),           &
-          kind = 'temperature'                                                 &
-       )
-
-       !specific humidity including nudging and subsidence
-       io_group_list(4) = init_io_group(                                       &
-          in_files = flow_files,                                               &
-          out_vars = [output_var_table(9:14), output_var_table(52:54)],        &
-          in_var_list = input_var_table(4:4),                                  &
-          kind = 'scalar'                                                      &
-       )
-
-       !u and v velocity and geostrophic winds ug, vg
+          out_vars = [output_var_table(56:64),                                 & ! internal averaged density and pressure profiles
+                      output_var_table(3:8), output_var_table(9:14),           &
+                      output_var_table(42:42), output_var_table(43:44),        &
+                      output_var_table(49:51), output_var_table(52:54)],       &
+          in_var_list = (/input_var_table(3), input_var_table(17),             & ! T, P, QV
+                          input_var_table(4) /),                               &
+          kind = 'thermodynamics',                                             &
+          n_output_quantities = 4                                              & ! P, Theta, Rho, qv
+       )
+
+       !Moved to therodynamic io_group
+       !io_group_list(4) = init_io_group(                                       &
+       !   in_files = flow_files,                                               &
+       !   out_vars = [output_var_table(9:14), output_var_table(52:54)],        &
+       !   in_var_list = input_var_table(4:4),                                  &
+       !   kind = 'scalar'                                                      &
+       !)
+
+       !u and v velocity, ug anv vg moved to thermodynamic io group. 
        io_group_list(5) = init_io_group(                                       &
           in_files = flow_files,                                               &
-          out_vars = [output_var_table(15:26), output_var_table(43:46)],       &
+          out_vars = [output_var_table(15:26), output_var_table(45:46)],       &
           !out_vars = output_var_table(15:20),                                  &
           in_var_list = input_var_table(5:6),                                  &
@@ -1845,9 +2118,11 @@
 
 
-    FUNCTION init_io_group(in_files, out_vars, in_var_list, kind) RESULT(group)
+    FUNCTION init_io_group(in_files, out_vars, in_var_list, kind,              &
+                           n_output_quantities) RESULT(group)
        CHARACTER(LEN=PATH), INTENT(IN) ::  in_files(:)
        CHARACTER(LEN=*), INTENT(IN)    ::  kind
        TYPE(nc_var), INTENT(IN)        ::  out_vars(:)
        TYPE(nc_var), INTENT(IN)        ::  in_var_list(:)
+       INTEGER, OPTIONAL               ::  n_output_quantities
 
        TYPE(io_group)                  ::  group
@@ -1855,7 +2130,14 @@
        group % nt = SIZE(in_files)
        group % nv = SIZE(out_vars)
+       group % n_inputs = SIZE(in_var_list)
        group % kind = TRIM(kind)
 
-       ALLOCATE(group % in_var_list( SIZE(in_var_list) ))
+       IF ( PRESENT(n_output_quantities) )  THEN
+          group % n_output_quantities = n_output_quantities
+       ELSE
+          group % n_output_quantities = group % n_inputs
+       END IF
+
+       ALLOCATE(group % in_var_list(group % n_inputs))
        ALLOCATE(group % in_files(group % nt))
        ALLOCATE(group % out_vars(group % nv))
@@ -1876,5 +2158,5 @@
     SUBROUTINE fini_grids()
 
-       CALL report('fini_grids', 'Deallocating grids')
+       CALL report('fini_grids', 'Deallocating grids', cfg % debug)
        
        DEALLOCATE(x, y, z, xu, yv, zw, z_column, zw_column)
@@ -1904,4 +2186,6 @@
     SUBROUTINE setup_variable_tables(ic_mode)
        CHARACTER(LEN=*), INTENT(IN) ::  ic_mode
+       INTEGER                      ::  n_invar = 0  !< number of variables in the input variable table
+       INTEGER                      ::  n_outvar = 0 !< number of variables in the output variable table
        TYPE(nc_var), POINTER        ::  var
 
@@ -1914,5 +2198,5 @@
 
        n_invar = 17
-       n_outvar = 55
+       n_outvar = 64
        ALLOCATE( input_var_table(n_invar) )
        ALLOCATE( output_var_table(n_outvar) )
@@ -2052,6 +2336,5 @@
        )
        IF (TRIM(ic_mode) == 'profile')  THEN
-          output_var_table(3) % grid => scalar_profile_grid
-          output_var_table(3) % intermediate_grid => scalar_profile_intermediate
+          output_var_table(3) % averaging_grid => averaged_scalar_profile
        END IF
 
@@ -2122,5 +2405,5 @@
           units             = "kg/kg",                                         &
           kind              = "init scalar",                                   &
-          input_id          = 1,                                               &
+          input_id          = 3,                                               &
           output_file       = output_file,                                     &
           grid              = palm_grid,                                       &
@@ -2129,6 +2412,5 @@
        )
        IF (TRIM(ic_mode) == 'profile')  THEN
-          output_var_table(9) % grid => scalar_profile_grid
-          output_var_table(9) % intermediate_grid => scalar_profile_intermediate
+          output_var_table(9) % averaging_grid => averaged_scalar_profile
        END IF
 
@@ -2139,5 +2421,5 @@
           units             = "kg/kg",                                         &
           kind              = "left scalar",                                   &
-          input_id          = 1,                                               &
+          input_id          = 3,                                               &
           output_file       = output_file,                                     &
           grid              = scalars_west_grid,                               &
@@ -2151,5 +2433,5 @@
           units             = "kg/kg",                                         &
           kind              = "right scalar",                                  &
-          input_id          = 1,                                               &
+          input_id          = 3,                                               &
           output_file       = output_file,                                     &
           grid              = scalars_east_grid,                               &
@@ -2163,5 +2445,5 @@
           units             = "kg/kg",                                         &
           kind              = "north scalar",                                  &
-          input_id          = 1,                                               &
+          input_id          = 3,                                               &
           output_file       = output_file,                                     &
           grid              = scalars_north_grid,                              &
@@ -2175,5 +2457,5 @@
           units             = "kg/kg",                                         &
           kind              = "south scalar",                                  &
-          input_id          = 1,                                               &
+          input_id          = 3,                                               &
           output_file       = output_file,                                     &
           grid              = scalars_south_grid,                              &
@@ -2187,5 +2469,5 @@
           units             = "kg/kg",                                         &
           kind              = "top scalar",                                    &
-          input_id          = 1,                                               &
+          input_id          = 3,                                               &
           output_file       = output_file,                                     &
           grid              = scalars_top_grid,                                &
@@ -2206,6 +2488,5 @@
        )
        IF (TRIM(ic_mode) == 'profile')  THEN
-          output_var_table(15) % grid => scalar_profile_grid
-          output_var_table(15) % intermediate_grid => scalar_profile_intermediate
+          output_var_table(15) % averaging_grid => averaged_scalar_profile
        END IF 
 
@@ -2283,6 +2564,5 @@
        )
        IF (TRIM(ic_mode) == 'profile')  THEN
-          output_var_table(21) % grid => scalar_profile_grid
-          output_var_table(21) % intermediate_grid => scalar_profile_intermediate
+          output_var_table(21) % averaging_grid => averaged_scalar_profile
        END IF
 
@@ -2360,6 +2640,5 @@
        )
        IF (TRIM(ic_mode) == 'profile')  THEN
-          output_var_table(27) % grid => w_profile_grid
-          output_var_table(27) % intermediate_grid => w_profile_intermediate
+          output_var_table(27) % averaging_grid => averaged_w_profile
        END IF
 
@@ -2546,4 +2825,5 @@
           intermediate_grid = palm_intermediate                                &
        )
+       output_var_table(42) % averaging_grid => averaged_scalar_profile
 
        output_var_table(43) = init_nc_var(                                     &
@@ -2552,9 +2832,9 @@
           long_name         = "geostrophic wind (u component)",                &
           units             = "m/s",                                           &
-          kind              = "constant scalar profile",                       &
+          kind              = "geostrophic",                                   &
           input_id          = 1,                                               &
           output_file       = output_file,                                     &
-          grid              = scalar_profile_grid,                             &
-          intermediate_grid = scalar_profile_intermediate                      &
+          grid              = averaged_scalar_profile,                         &
+          intermediate_grid = averaged_scalar_profile                          &
        )
 
@@ -2564,9 +2844,9 @@
           long_name         = "geostrophic wind (v component)",                &
           units             = "m/s",                                           &
-          kind              = "constant scalar profile",                       &
+          kind              = "geostrophic",                                   &
           input_id          = 1,                                               &
           output_file       = output_file,                                     &
-          grid              = scalar_profile_grid,                             &
-          intermediate_grid = scalar_profile_intermediate                      &
+          grid              = averaged_scalar_profile,                         &
+          intermediate_grid = averaged_scalar_profile                          &
        )
 
@@ -2576,9 +2856,9 @@
           long_name         = "wind component in x direction",                 &
           units             = "m/s",                                           &
-          kind              = "large-scale scalar forcing",                    &
+          kind              = "geostrophic",                                   &
           input_id          = 1,                                               &
           output_file       = output_file,                                     &
-          grid              = scalar_profile_grid,                             &
-          intermediate_grid = scalar_profile_intermediate                      &
+          grid              = averaged_scalar_profile,                         &
+          intermediate_grid = averaged_scalar_profile                          &
        )
        output_var_table(45) % to_be_processed = ls_forcing_variables_required
@@ -2592,6 +2872,6 @@
           input_id          = 1,                                               &
           output_file       = output_file,                                     &
-          grid              = scalar_profile_grid,                             &
-          intermediate_grid = scalar_profile_intermediate                      &
+          grid              = averaged_scalar_profile,                         &
+          intermediate_grid = averaged_scalar_profile                          &
        )
        output_var_table(46) % to_be_processed = ls_forcing_variables_required
@@ -2605,6 +2885,6 @@
           input_id          = 1,                                               &
           output_file       = output_file,                                     &
-          grid              = w_profile_grid,                                  &
-          intermediate_grid = w_profile_intermediate                           &
+          grid              = averaged_scalar_profile,                         &
+          intermediate_grid = averaged_scalar_profile                          &
        )
        output_var_table(47) % to_be_processed = ls_forcing_variables_required
@@ -2618,6 +2898,6 @@
           input_id          = 1,                                               &
           output_file       = output_file,                                     &
-          grid              = w_profile_grid,                                  &
-          intermediate_grid = w_profile_intermediate                           &
+          grid              = averaged_w_profile,                              &
+          intermediate_grid = averaged_w_profile                               &
        )
        output_var_table(48) % to_be_processed = ls_forcing_variables_required
@@ -2632,6 +2912,6 @@
           input_id          = 1,                                               &
           output_file       = output_file,                                     &
-          grid              = scalar_profile_grid,                             &
-          intermediate_grid = scalar_profile_intermediate                      &
+          grid              = averaged_scalar_profile,                         &
+          intermediate_grid = averaged_scalar_profile                          &
        )
        output_var_table(49) % to_be_processed = ls_forcing_variables_required
@@ -2645,6 +2925,6 @@
           input_id          = 1,                                               &
           output_file       = output_file,                                     &
-          grid              = scalar_profile_grid,                             &
-          intermediate_grid = scalar_profile_intermediate                      &
+          grid              = averaged_scalar_profile,                         &
+          intermediate_grid = averaged_scalar_profile                          &
        )
        output_var_table(50) % to_be_processed = ls_forcing_variables_required
@@ -2658,6 +2938,6 @@
           input_id          = 1,                                               &
           output_file       = output_file,                                     &
-          grid              = scalar_profile_grid,                             &
-          intermediate_grid = scalar_profile_intermediate                      &
+          grid              = averaged_scalar_profile,                         &
+          intermediate_grid = averaged_scalar_profile                          &
        )
        output_var_table(51) % to_be_processed = ls_forcing_variables_required
@@ -2669,8 +2949,8 @@
           units             = "kg/kg/s",                                       &
           kind              = "large-scale scalar forcing",                    &
-          input_id          = 1,                                               &
-          output_file       = output_file,                                     &
-          grid              = scalar_profile_grid,                             &
-          intermediate_grid = scalar_profile_intermediate                      &
+          input_id          = 3,                                               &
+          output_file       = output_file,                                     &
+          grid              = averaged_scalar_profile,                         &
+          intermediate_grid = averaged_scalar_profile                          &
        )
        output_var_table(52) % to_be_processed = ls_forcing_variables_required
@@ -2683,8 +2963,8 @@
           units             = "kg/kg/s",                                       &
           kind              = "large-scale scalar forcing",                    &
-          input_id          = 1,                                               &
-          output_file       = output_file,                                     &
-          grid              = scalar_profile_grid,                             &
-          intermediate_grid = scalar_profile_intermediate                      &
+          input_id          = 3,                                               &
+          output_file       = output_file,                                     &
+          grid              = averaged_scalar_profile,                         &
+          intermediate_grid = averaged_scalar_profile                          &
        )
        output_var_table(53) % to_be_processed = ls_forcing_variables_required
@@ -2696,8 +2976,8 @@
           units             = "kg/kg",                                         &
           kind              = "large-scale scalar forcing",                    &
-          input_id          = 1,                                               &
-          output_file       = output_file,                                     &
-          grid              = scalar_profile_grid,                             &
-          intermediate_grid = scalar_profile_intermediate                      &
+          input_id          = 3,                                               &
+          output_file       = output_file,                                     &
+          grid              = averaged_scalar_profile,                         &
+          intermediate_grid = averaged_scalar_profile                          &
        )
        output_var_table(54) % to_be_processed = ls_forcing_variables_required
@@ -2711,9 +2991,141 @@
           input_id          = 1,                                               &
           output_file       = output_file,                                     &
-          grid              = scalar_profile_grid,                             &
-          intermediate_grid = scalar_profile_intermediate                      &
+          grid              = averaged_scalar_profile,                         &
+          intermediate_grid = averaged_scalar_profile                          &
        )
        output_var_table(55) % to_be_processed = ls_forcing_variables_required
 
+
+       output_var_table(56) = init_nc_var(                                     &
+          name              = 'internal_density_centre',                              &
+          std_name          = "",                                              &
+          long_name         = "",                                              &
+          units             = "",                                              &
+          kind              = "internal profile",                              &
+          input_id          = 4,                                               &
+          output_file       = output_file,                                     &
+          grid              = averaged_scalar_profile,                         &
+          intermediate_grid = averaged_scalar_profile                          &
+       )
+       output_var_table(56) % averaging_grid => averaged_scalar_profile
+
+
+       output_var_table(57) = init_nc_var(                                     &
+          name              = 'internal_density_north',                       &
+          std_name          = "",                                              &
+          long_name         = "",                                              &
+          units             = "",                                              &
+          kind              = "internal profile",                              &
+          input_id          = 4,                                               &
+          output_file       = output_file,                                     &
+          grid              = north_averaged_scalar_profile,                   &
+          intermediate_grid = north_averaged_scalar_profile                    &
+       )
+       output_var_table(57) % averaging_grid => north_averaged_scalar_profile
+       output_var_table(57) % to_be_processed = .NOT. cfg % ug_is_set
+
+
+       output_var_table(58) = init_nc_var(                                     &
+          name              = 'internal_density_south',                       &
+          std_name          = "",                                              &
+          long_name         = "",                                              &
+          units             = "",                                              &
+          kind              = "internal profile",                              &
+          input_id          = 4,                                               &
+          output_file       = output_file,                                     &
+          grid              = south_averaged_scalar_profile,                   &
+          intermediate_grid = south_averaged_scalar_profile                    &
+       )
+       output_var_table(58) % averaging_grid => south_averaged_scalar_profile
+       output_var_table(58) % to_be_processed = .NOT. cfg % ug_is_set
+
+
+       output_var_table(59) = init_nc_var(                                     &
+          name              = 'internal_density_east',                        &
+          std_name          = "",                                              &
+          long_name         = "",                                              &
+          units             = "",                                              &
+          kind              = "internal profile",                              &
+          input_id          = 4,                                               &
+          output_file       = output_file,                                     &
+          grid              = east_averaged_scalar_profile,                    &
+          intermediate_grid = east_averaged_scalar_profile                     &
+       )
+       output_var_table(59) % averaging_grid => east_averaged_scalar_profile
+       output_var_table(59) % to_be_processed = .NOT. cfg % ug_is_set
+
+
+       output_var_table(60) = init_nc_var(                                     &
+          name              = 'internal_density_west',                        &
+          std_name          = "",                                              &
+          long_name         = "",                                              &
+          units             = "",                                              &
+          kind              = "internal profile",                              &
+          input_id          = 4,                                               &
+          output_file       = output_file,                                     &
+          grid              = west_averaged_scalar_profile,                    &
+          intermediate_grid = west_averaged_scalar_profile                     &
+       )
+       output_var_table(60) % averaging_grid => west_averaged_scalar_profile
+       output_var_table(60) % to_be_processed = .NOT. cfg % ug_is_set
+
+       output_var_table(61) = init_nc_var(                                     &
+          name              = 'internal_pressure_north',                       &
+          std_name          = "",                                              &
+          long_name         = "",                                              &
+          units             = "",                                              &
+          kind              = "internal profile",                              &
+          input_id          = 2,                                               &
+          output_file       = output_file,                                     &
+          grid              = north_averaged_scalar_profile,                   &
+          intermediate_grid = north_averaged_scalar_profile                    &
+       )
+       output_var_table(61) % averaging_grid => north_averaged_scalar_profile
+       output_var_table(61) % to_be_processed = .NOT. cfg % ug_is_set
+
+
+       output_var_table(62) = init_nc_var(                                     &
+          name              = 'internal_pressure_south',                       &
+          std_name          = "",                                              &
+          long_name         = "",                                              &
+          units             = "",                                              &
+          kind              = "internal profile",                              &
+          input_id          = 2,                                               &
+          output_file       = output_file,                                     &
+          grid              = south_averaged_scalar_profile,                   &
+          intermediate_grid = south_averaged_scalar_profile                    &
+       )
+       output_var_table(62) % averaging_grid => south_averaged_scalar_profile
+       output_var_table(62) % to_be_processed = .NOT. cfg % ug_is_set
+
+
+       output_var_table(63) = init_nc_var(                                     &
+          name              = 'internal_pressure_east',                        &
+          std_name          = "",                                              &
+          long_name         = "",                                              &
+          units             = "",                                              &
+          kind              = "internal profile",                              &
+          input_id          = 2,                                               &
+          output_file       = output_file,                                     &
+          grid              = east_averaged_scalar_profile,                    &
+          intermediate_grid = east_averaged_scalar_profile                     &
+       )
+       output_var_table(63) % averaging_grid => east_averaged_scalar_profile
+       output_var_table(63) % to_be_processed = .NOT. cfg % ug_is_set
+
+
+       output_var_table(64) = init_nc_var(                                     &
+          name              = 'internal_pressure_west',                        &
+          std_name          = "",                                              &
+          long_name         = "",                                              &
+          units             = "",                                              &
+          kind              = "internal profile",                              &
+          input_id          = 2,                                               &
+          output_file       = output_file,                                     &
+          grid              = west_averaged_scalar_profile,                    &
+          intermediate_grid = west_averaged_scalar_profile                     &
+       )
+       output_var_table(64) % averaging_grid => west_averaged_scalar_profile
+       output_var_table(64) % to_be_processed = .NOT. cfg % ug_is_set
 
        ! Attributes shared among all variables
@@ -2732,6 +3144,6 @@
 !------------------------------------------------------------------------------!
     FUNCTION init_nc_var(name, std_name, long_name, units, kind, input_id,     &
-                         grid, intermediate_grid, output_file, is_profile      &
-             ) RESULT(var)
+                         grid, intermediate_grid, output_file, is_profile)     &
+       RESULT(var)
 
        CHARACTER(LEN=*), INTENT(IN)      ::  name, std_name, long_name, units, kind
@@ -2772,4 +3184,5 @@
           var % dimvarids(1:3)  = output_file % dimvarids_soil
           var % to_be_processed = init_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "interpolate_2d"
 
@@ -2781,4 +3194,5 @@
           var % dimvarids(1:3)  = output_file % dimvarids_scl
           var % to_be_processed = init_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "interpolate_3d"
 
@@ -2794,4 +3208,5 @@
           var % dimvarids(3)    = output_file % dimvarids_scl(3)
           var % to_be_processed = init_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "interpolate_3d"
 
@@ -2807,4 +3222,5 @@
           var % dimvarids(3)    = output_file % dimvarids_scl(3)
           var % to_be_processed = init_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "interpolate_3d"
 
@@ -2820,4 +3236,5 @@
           var % dimvarids(3)    = output_file % dimvarids_vel(3)
           var % to_be_processed = init_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "interpolate_3d"
 
@@ -2829,4 +3246,5 @@
           var % dimvarids(1)    = output_file % dimvarids_scl(3) !z
           var % to_be_processed = init_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "average profile"
 
@@ -2838,4 +3256,5 @@
           var % dimvarids(1)    = output_file % dimvarids_vel(3) !z
           var % to_be_processed = init_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "average profile"
 
@@ -2848,4 +3267,5 @@
           var % dimvarids(1:2)  = output_file % dimvarids_soil(1:2)
           var % to_be_processed = boundary_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "interpolate_2d"
 
@@ -2860,4 +3280,5 @@
           var % dimvarids(2)    = output_file % dimvarids_scl(3)
           var % to_be_processed = boundary_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "interpolate_3d"
 
@@ -2872,4 +3293,5 @@
           var % dimvarids(2)    = output_file % dimvarids_scl(3)
           var % to_be_processed = boundary_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "interpolate_3d"
 
@@ -2884,4 +3306,5 @@
           var % dimvarids(2)    = output_file % dimvarids_scl(2)
           var % to_be_processed = boundary_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "interpolate_3d"
 
@@ -2896,4 +3319,5 @@
           var % dimvarids(2)    = output_file % dimvarids_scl(3)
           var % to_be_processed = boundary_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "interpolate_3d"
 
@@ -2908,4 +3332,5 @@
           var % dimvarids(2)    = output_file % dimvarids_scl(3)
           var % to_be_processed = boundary_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "interpolate_3d"
 
@@ -2920,4 +3345,5 @@
           var % dimvarids(2)    = output_file % dimvarids_scl(2)
           var % to_be_processed = boundary_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "interpolate_3d"
 
@@ -2932,4 +3358,5 @@
           var % dimvarids(2)    = output_file % dimvarids_scl(3)
           var % to_be_processed = boundary_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "interpolate_3d"
 
@@ -2944,4 +3371,5 @@
           var % dimvarids(2)    = output_file % dimvarids_scl(3)
           var % to_be_processed = boundary_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "interpolate_3d"
 
@@ -2956,4 +3384,5 @@
           var % dimvarids(2)    = output_file % dimvarids_vel(2)
           var % to_be_processed = boundary_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "interpolate_3d"
 
@@ -2968,4 +3397,5 @@
           var % dimvarids(2)    = output_file % dimvarids_vel(3)
           var % to_be_processed = boundary_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "interpolate_3d"
 
@@ -2980,4 +3410,5 @@
           var % dimvarids(2)    = output_file % dimvarids_vel(3)
           var % to_be_processed = boundary_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "interpolate_3d"
 
@@ -2988,5 +3419,6 @@
           var % dimvarids(1)    = output_file % dimvarid_time
           var % to_be_processed = .TRUE.
-          var % task            = "average scalar"
+          var % is_internal     = .FALSE.
+          var % task            = "average profile"
 
        CASE( 'constant scalar profile' )
@@ -2998,4 +3430,5 @@
           var % dimvarids(1)    = output_file % dimvarids_scl(3)
           var % to_be_processed = .TRUE.
+          var % is_internal     = .FALSE.
           var % task            = "set profile"
 
@@ -3008,5 +3441,17 @@
           var % dimvarids(1)    = output_file % dimvarids_scl(3)
           var % to_be_processed = ls_forcing_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "average large-scale profile"
+
+       CASE( 'geostrophic' )
+          var % lod             = -1
+          var % ndim            = 2
+          var % dimids(2)       = output_file % dimid_time    !t
+          var % dimids(1)       = output_file % dimids_scl(3) !z
+          var % dimvarids(2)    = output_file % dimvarid_time
+          var % dimvarids(1)    = output_file % dimvarids_scl(3)
+          var % to_be_processed = .TRUE.
+          var % is_internal     = .FALSE.
+          var % task            = "geostrophic winds"
 
        CASE( 'large-scale w forcing' )
@@ -3018,5 +3463,17 @@
           var % dimvarids(1)    = output_file % dimvarids_vel(3)
           var % to_be_processed = ls_forcing_variables_required
+          var % is_internal     = .FALSE.
           var % task            = "average large-scale profile"
+
+       CASE( 'internal profile' )
+          var % lod             = -1
+          var % ndim            = 2
+          var % dimids(2)       = output_file % dimid_time    !t
+          var % dimids(1)       = output_file % dimids_scl(3) !z
+          var % dimvarids(2)    = output_file % dimvarid_time
+          var % dimvarids(1)    = output_file % dimvarids_scl(3)
+          var % to_be_processed = .TRUE.
+          var % is_internal     = .TRUE.
+          var % task            = "internal profile"
 
        CASE DEFAULT
@@ -3031,5 +3488,5 @@
     SUBROUTINE fini_variables()
 
-       CALL report('fini_variables', 'Deallocating variable table')
+       CALL report('fini_variables', 'Deallocating variable table', cfg % debug)
        DEALLOCATE( input_var_table )
 
@@ -3039,5 +3496,5 @@
     SUBROUTINE fini_io_groups()
 
-       CALL report('fini_io_groups', 'Deallocating IO groups')
+       CALL report('fini_io_groups', 'Deallocating IO groups', cfg % debug)
        DEALLOCATE( io_group_list )
 
@@ -3047,5 +3504,5 @@
     SUBROUTINE fini_file_lists()
        
-       CALL report('fini_file_lists', 'Deallocating file lists')
+       CALL report('fini_file_lists', 'Deallocating file lists', cfg % debug)
        DEALLOCATE( flow_files, soil_files, radiation_files, soil_moisture_files )
 
@@ -3134,6 +3591,6 @@
              ! rotated velocities
              DO k = 1, nz
-             DO j = 2, ny
-             DO i = 2, nx
+             DO j = 1, ny
+             DO i = 1, nx
                 CALL uv2uvrot( urot = preprocess_buffer(1) % array(i,j,k),     &
                                vrot = preprocess_buffer(2) % array(i,j,k),     &
@@ -3171,5 +3628,5 @@
           CALL report('preprocess', message)
        
-       CASE( 'temperature' ) ! P and T
+       CASE( 'thermodynamics' ) ! T, P, QV
           nx = SIZE(input_buffer(1) % array, 1)
           ny = SIZE(input_buffer(1) % array, 2)
@@ -3192,5 +3649,6 @@
 
                 input_buffer (2) % array(i,j,:) =                              &
-                   input_buffer (2) % array(i,j,:) + basic_state_pressure(:)
+                   HECTO * input_buffer (2) % array(i,j,:) +                   &
+                   basic_state_pressure(:)
 
              END DO
@@ -3204,10 +3662,42 @@
 
           END IF
+          ! pressure
+          input_buffer(2) % is_preprocessed = .TRUE.
+
+          ! Copy temperature to last input buffer array
+          ALLOCATE(                                                            &
+              input_buffer( group % n_output_quantities ) % array (nx, ny, nz) &
+          )
+          input_buffer(group % n_output_quantities) % array(:,:,:) =           &
+              input_buffer(1) % array(:,:,:)
 
           ! Convert absolute to potential temperature
-          input_buffer(1) % array(:,:,:) = input_buffer(1) % array(:,:,:) *    &
-             EXP( RD_PALM/CP_PALM * LOG( P_REF / input_buffer(2) % array(:,:,:) ))
-
-          input_buffer(1:2) % is_preprocessed = .TRUE.
+          CALL potential_temperature(                                          &
+             t = input_buffer(1) % array(:,:,:),                               &
+             p = input_buffer(2) % array(:,:,:),                               &
+             p_ref = P_REF,                                                    &
+             r = RD_PALM,                                                      &
+             cp = CP_PALM                                                      &
+          )
+
+          ! potential temperature
+          input_buffer(1) % is_preprocessed = .TRUE.
+
+          ! Convert temperature copy to density
+          CALL moist_density(                                                  &
+             t_rho = input_buffer(group % n_output_quantities) % array(:,:,:), &
+             p = input_buffer(2) % array(:,:,:),                               &
+             qv = input_buffer(3) % array(:,:,:),                              &
+             rd = RD,                                                          &
+             rv = RV                                                           &
+          )
+
+          ! qv
+          input_buffer(3) % is_preprocessed = .TRUE.
+
+          ! density
+          input_buffer(group % n_output_quantities) % is_preprocessed = .TRUE.
+
+
           message = "Input buffers for group '" // TRIM(group % kind) // "'"//&
              " preprocessed sucessfully."
@@ -3336,5 +3826,5 @@
 
        CASE DEFAULT
-          message = "Buffer kind '" // TRIM(group % kind) // "' is not supported."
+          message = "IO group kind '" // TRIM(group % kind) // "' is not supported."
           CALL abort('prerpocess', message)
 
