Index: /palm/trunk/UTIL/inifor/src/inifor.f90
===================================================================
--- /palm/trunk/UTIL/inifor/src/inifor.f90	(revision 4755)
+++ /palm/trunk/UTIL/inifor/src/inifor.f90	(revision 4756)
@@ -26,4 +26,8 @@
 ! -----------------
 ! $Id$
+! Improved variable naming and code readabilty
+!
+!
+! 4675 2020-09-11 10:00:26Z eckhard
 ! Support for soil profile initialization
 ! Improved code formatting
@@ -157,19 +161,19 @@
 
     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)     ::  output_arr !< array buffer for interpolated quantities
-    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_centre !< density profile of the centre averaging domain
+    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_centre_cosmo !< density profile of the centre averaging domain
     REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  ug_cosmo   !< profile of the geostrophic wind in x direction on COSMO levels
     REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  vg_cosmo   !< profile of the geostrophic wind in y direction on COSMO levels
     REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  ug_palm    !< profile of the geostrophic wind in x direction interpolated onto PALM levels
     REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  vg_palm    !< profile of the geostrophic wind in y direction interpolated onto PALM levels
-    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_north  !< density profile of the northern averaging domain
-    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_south  !< density profile of the southern averaging domain
-    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_east   !< density profile of the eastern averaging domain
-    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_west   !< density profile of the western averaging domain
-    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_north    !< pressure profile of the northern averaging domain
-    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_south    !< pressure profile of the southern averaging domain
-    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_east     !< pressure profile of the eastern averaging domain
-    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_west     !< pressure profile of the western averaging domain
-
-    REAL(wp), POINTER, DIMENSION(:) ::  internal_arr !< pointer to the currently processed internal array (density, pressure)
+    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_north_cosmo  !< density profile of the northern averaging domain
+    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_south_cosmo  !< density profile of the southern averaging domain
+    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_east_cosmo   !< density profile of the eastern averaging domain
+    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_west_cosmo   !< density profile of the western averaging domain
+    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_north_cosmo    !< pressure profile of the northern averaging domain
+    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_south_cosmo    !< pressure profile of the southern averaging domain
+    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_east_cosmo     !< pressure profile of the eastern averaging domain
+    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_west_cosmo     !< pressure profile of the western averaging domain
+
+    REAL(wp), POINTER, DIMENSION(:) ::  internal_cosmo_arr !< pointer to the currently processed internal array (density, pressure)
     REAL(wp), POINTER, DIMENSION(:) ::  ug_vg_palm   !< pointer to the currently processed geostrophic wind component
 
@@ -334,5 +338,6 @@
                             ELSE
                                CALL get_surface_pressure(                      &
-                                  output_arr(0,0,:), rho_centre,               &
+                                  output_arr(0,0,:),                           &
+                                  rho_centre_cosmo,                            &
                                   output_var%averaging_grid )
                             ENDIF
@@ -362,38 +367,38 @@
    
                          CASE( 'internal_density_centre' )
-                            ALLOCATE( rho_centre(1:cosmo_grid%nz) )
-                            internal_arr => rho_centre
+                            ALLOCATE( rho_centre_cosmo(1:cosmo_grid%nz) )
+                            internal_cosmo_arr => rho_centre_cosmo
    
                          CASE( 'internal_density_north' )
-                            ALLOCATE( rho_north(1:cosmo_grid%nz) )
-                            internal_arr => rho_north
+                            ALLOCATE( rho_north_cosmo(1:cosmo_grid%nz) )
+                            internal_cosmo_arr => rho_north_cosmo
    
                          CASE( 'internal_density_south' )
-                            ALLOCATE( rho_south(1:cosmo_grid%nz) )
-                            internal_arr => rho_south
+                            ALLOCATE( rho_south_cosmo(1:cosmo_grid%nz) )
+                            internal_cosmo_arr => rho_south_cosmo
    
                          CASE( 'internal_density_east' )
-                            ALLOCATE( rho_east(1:cosmo_grid%nz) )
-                            internal_arr => rho_east
+                            ALLOCATE( rho_east_cosmo(1:cosmo_grid%nz) )
+                            internal_cosmo_arr => rho_east_cosmo
    
                          CASE( 'internal_density_west' )
-                            ALLOCATE( rho_west(1:cosmo_grid%nz) )
-                            internal_arr => rho_west
+                            ALLOCATE( rho_west_cosmo(1:cosmo_grid%nz) )
+                            internal_cosmo_arr => rho_west_cosmo
    
                          CASE( 'internal_pressure_north' )
-                            ALLOCATE( p_north(1:cosmo_grid%nz) )
-                            internal_arr => p_north
+                            ALLOCATE( p_north_cosmo(1:cosmo_grid%nz) )
+                            internal_cosmo_arr => p_north_cosmo
    
                          CASE( 'internal_pressure_south' )
-                            ALLOCATE( p_south(1:cosmo_grid%nz) )
-                            internal_arr => p_south
+                            ALLOCATE( p_south_cosmo(1:cosmo_grid%nz) )
+                            internal_cosmo_arr => p_south_cosmo
    
                          CASE( 'internal_pressure_east' )
-                            ALLOCATE( p_east(1:cosmo_grid%nz) )
-                            internal_arr => p_east
+                            ALLOCATE( p_east_cosmo(1:cosmo_grid%nz) )
+                            internal_cosmo_arr => p_east_cosmo
    
                          CASE( 'internal_pressure_west' )
-                            ALLOCATE( p_west(1:cosmo_grid%nz) )
-                            internal_arr => p_west
+                            ALLOCATE( p_west_cosmo(1:cosmo_grid%nz) )
+                            internal_cosmo_arr => p_west_cosmo
    
                          CASE DEFAULT
@@ -413,5 +418,5 @@
                             CALL average_pressure_perturbation(                &
                                input_buffer(output_var%input_id)%array(:,:,:), &
-                               internal_arr(:),                                &
+                               internal_cosmo_arr(:),                          &
                                cosmo_grid, output_var%averaging_grid           &
                             )
@@ -421,5 +426,5 @@
                             CALL average_profile(                              &
                                input_buffer(output_var%input_id)%array(:,:,:), &
-                               internal_arr(:),                                &
+                               internal_cosmo_arr(:),                          &
                                output_var%averaging_grid                       &
                             )
@@ -445,11 +450,13 @@
                             vg_palm = cfg%vg
                          ELSE
-                            CALL geostrophic_winds( p_north, p_south, p_east,  &
-                                                    p_west, rho_centre, f3,    &
-                                                    averaging_width_ew,        &
-                                                    averaging_width_ns,        &
-                                                    phi_n, lambda_n,           &
-                                                    phi_centre, lam_centre,    &
-                                                    ug_cosmo, vg_cosmo )
+                            CALL geostrophic_winds(                            &
+                               p_north_cosmo, p_south_cosmo, p_east_cosmo,     &
+                               p_west_cosmo, rho_centre_cosmo, f3,             &
+                               averaging_width_ew,                             &
+                               averaging_width_ns,                             &
+                               phi_n, lambda_n,                                &
+                               phi_centre, lam_centre,                         &
+                               ug_cosmo, vg_cosmo                              &
+                            )
 
                             CALL interpolate_1d( ug_cosmo, ug_palm,            &
@@ -532,5 +539,5 @@
              ug_vg_have_been_computed = .FALSE.
              IF ( group%kind == 'thermodynamics' )  THEN
-                DEALLOCATE( rho_centre )
+                DEALLOCATE( rho_centre_cosmo )
                 DEALLOCATE( ug_palm )
                 DEALLOCATE( vg_palm )
@@ -538,12 +545,12 @@
                 DEALLOCATE( vg_cosmo )
                 IF ( .NOT. cfg%ug_defined_by_user )  THEN
-                   DEALLOCATE( rho_north )
-                   DEALLOCATE( rho_south )
-                   DEALLOCATE( rho_east )
-                   DEALLOCATE( rho_west )
-                   DEALLOCATE( p_north )
-                   DEALLOCATE( p_south )
-                   DEALLOCATE( p_east )
-                   DEALLOCATE( p_west )
+                   DEALLOCATE( rho_north_cosmo )
+                   DEALLOCATE( rho_south_cosmo )
+                   DEALLOCATE( rho_east_cosmo )
+                   DEALLOCATE( rho_west_cosmo )
+                   DEALLOCATE( p_north_cosmo )
+                   DEALLOCATE( p_south_cosmo )
+                   DEALLOCATE( p_east_cosmo )
+                   DEALLOCATE( p_west_cosmo )
                 ENDIF
              ENDIF
Index: /palm/trunk/UTIL/inifor/src/inifor_defs.f90
===================================================================
--- /palm/trunk/UTIL/inifor/src/inifor_defs.f90	(revision 4755)
+++ /palm/trunk/UTIL/inifor/src/inifor_defs.f90	(revision 4756)
@@ -214,5 +214,5 @@
     ACHAR( 10 ) // ' Copyright 2017-2020 Deutscher Wetterdienst Offenbach' !< Copyright notice
  CHARACTER(LEN=*), PARAMETER ::  LOG_FILE_NAME = 'inifor.log' !< Name of INIFOR's log file
- CHARACTER(LEN=*), PARAMETER ::  VERSION = '2.1.1'            !< INIFOR version number
+ CHARACTER(LEN=*), PARAMETER ::  VERSION = '2.1.2'            !< INIFOR version number
  
  END MODULE inifor_defs
Index: /palm/trunk/UTIL/inifor/src/inifor_transform.f90
===================================================================
--- /palm/trunk/UTIL/inifor/src/inifor_transform.f90	(revision 4755)
+++ /palm/trunk/UTIL/inifor/src/inifor_transform.f90	(revision 4756)
@@ -26,4 +26,10 @@
 ! -----------------
 ! $Id$
+! Fixed an error in surface pressure extrapolation where the cosmo grid was
+!    misinterpreted as palm grid
+! Improved code readability and formatting
+!
+!
+! 4714 2020-09-29 12:47:35Z eckhard
 ! Fixed off-by-one indexing error for profile quantities
 !
@@ -381,7 +387,7 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE interp_average_profile(source_array, profile_array, avg_grid)
-    TYPE(grid_definition), INTENT(IN)          ::  avg_grid
-    REAL(wp), DIMENSION(:,:,:), INTENT(IN)     ::  source_array(0:,0:,:)
-    REAL(wp), DIMENSION(:), INTENT(OUT)        ::  profile_array
+    TYPE(grid_definition), INTENT(IN)      ::  avg_grid
+    REAL(wp), DIMENSION(:,:,:), INTENT(IN) ::  source_array(0:,0:,:)
+    REAL(wp), DIMENSION(:), INTENT(OUT)    ::  profile_array
 
     INTEGER(iwp) ::  i_source, j_source, k_profile, k_source, l, m
@@ -438,7 +444,7 @@
  SUBROUTINE average_profile( source_array, profile_array, avg_grid )
 
-    TYPE(grid_definition), INTENT(IN)          ::  avg_grid
-    REAL(wp), DIMENSION(:,:,:), INTENT(IN)     ::  source_array(0:,0:,:)
-    REAL(wp), DIMENSION(:), INTENT(OUT)        ::  profile_array
+    TYPE(grid_definition), INTENT(IN) ::  avg_grid
+    REAL(wp), INTENT(IN)              ::  source_array(0:,0:,:)
+    REAL(wp), INTENT(OUT)             ::  profile_array(:)
 
     INTEGER(iwp) ::  i_source, j_source, l, nz, nlev
@@ -484,12 +490,11 @@
                                            cosmo_grid, avg_grid )
 
-    TYPE(grid_definition), INTENT(IN)          ::  cosmo_grid, avg_grid
-    REAL(wp), DIMENSION(:,:,:), INTENT(IN)     ::  cosmo_pressure(0:,0:,:)
-    REAL(wp), DIMENSION(:), INTENT(OUT)        ::  profile_array
-
-    INTEGER(iwp) ::  i_source, j_source, l, nz, nlev
-
-    REAL(wp)                            ::  ni_columns
-    REAL(wp), DIMENSION(:), ALLOCATABLE ::  basic_state_pressure
+    TYPE(grid_definition), INTENT(IN) ::  cosmo_grid, avg_grid
+    REAL(wp), INTENT(IN)              ::  cosmo_pressure(0:,0:,:)
+    REAL(wp), INTENT(OUT)             ::  profile_array(:)
+
+    INTEGER(iwp)          ::  i_source, j_source, l, nz, nlev
+    REAL(wp)              ::  ni_columns
+    REAL(wp), ALLOCATABLE ::  basic_state_pressure(:)
 
     nlev = SIZE( cosmo_pressure, 3 )
@@ -624,42 +629,59 @@
 ! Description:
 ! ------------
-!> Takes the averaged pressure profile <p> and sets the lowest entry to the
-!> extrapolated pressure at the surface.
-!------------------------------------------------------------------------------!
- SUBROUTINE get_surface_pressure(p, rho, avg_grid)
-    REAL(wp), DIMENSION(:), INTENT(IN)    ::  rho
-    REAL(wp), DIMENSION(:), INTENT(INOUT) ::  p
+!> Takes the averaged pressure profile p_palm and sets the lowest entry to the
+!> extrapolated pressure at the surface, assuming a linear density profile.
+!------------------------------------------------------------------------------!
+ SUBROUTINE get_surface_pressure(p_palm, rho_cosmo, avg_grid)
+    REAL(wp), DIMENSION(:), INTENT(IN)    ::  rho_cosmo
+    REAL(wp), DIMENSION(:), INTENT(INOUT) ::  p_palm
     TYPE(grid_definition), INTENT(IN)     ::  avg_grid
 
-    REAL(wp)     ::  drhodz, dz, zk, rhok
-    INTEGER(iwp) ::  k_min
-
-    k_min = avg_grid%k_min
-    zk    = avg_grid%z(k_min)
-    rhok  = rho(k_min)
-    dz    = avg_grid%z(k_min + 1) - avg_grid%z(k_min)
-    drhodz = ( rho(k_min + 1) - rho(k_min) ) / dz
-
-    p(1) = linear_density_pressure( p(k_min), zk, rhok, drhodz,                &
-                                    z = 0.0_wp, g=G )
+    REAL(wp)     ::  drhodz_surface_cosmo
+    INTEGER(iwp) ::  k_min_palm
+
+    drhodz_surface_cosmo =                                                     &
+       ( rho_cosmo(2) - rho_cosmo(1) ) /                                       &
+       ( avg_grid%intermediate_h(1,1,2) - avg_grid%intermediate_h(1,1,1) )
+
+    k_min_palm = avg_grid%k_min
+
+    p_palm(1) = linear_density_pressure(                                       &
+       p0 = p_palm(k_min_palm),                                                &
+       z0 = avg_grid%z(k_min_palm),                                            &
+       rho00 = rho_cosmo(1),                                                   &
+       z00 = avg_grid%intermediate_h(1,1,1),                                   &
+       drhodz = drhodz_surface_cosmo,                                          &
+       g = G,                                                                  &
+       z = 0.0_wp                                                              &
+    )
 
  END SUBROUTINE get_surface_pressure
 
 
- FUNCTION linear_density_pressure(pk, zk, rhok, drhodz, z, g)  RESULT(p)
-
-    REAL(wp), INTENT(IN)  ::  pk, zk, rhok, drhodz, g, z
+!------------------------------------------------------------------------------!
+! Description:
+! -----------
+!> Computes the hydrostatic pressure p at height z given the pressure p0 at
+!> height z0. The pressure is computed based on the solution of the hydrostatic
+!> equation assuming a linear density profile with density rho00 at z00 and the
+!> vertical density gradient drhodz. 
+!------------------------------------------------------------------------------!
+ FUNCTION linear_density_pressure(p0, z0, rho00, z00, drhodz, g, z)  RESULT(p)
+
+    REAL(wp), INTENT(IN)  ::  p0, z0, rho00, z00, drhodz, g, z
     REAL(wp) ::  p
 
-    p = pk + ( zk - z ) * g * ( rhok - 0.5_wp * drhodz * (zk - z) )
+    p = p0 - ( z - z0 ) * g * (                                                &
+           rho00 + 0.5_wp * drhodz * ( z + z0 - 2.0_wp * z00 )                 &
+        )
 
  END FUNCTION linear_density_pressure
 
-!-----------------------------------------------------------------------------!
+!------------------------------------------------------------------------------!
 ! Description:
 ! -----------
 !> This routine computes profiles of the two geostrophic wind components ug and
 !> vg.
-!-----------------------------------------------------------------------------!
+!------------------------------------------------------------------------------!
  SUBROUTINE geostrophic_winds(p_north, p_south, p_east, p_west, rho, f3,    &
                               Lx, Ly, phi_n, lam_n, phi_g, lam_g, ug, vg)
@@ -683,5 +705,5 @@
 
 
-!-----------------------------------------------------------------------------!
+!------------------------------------------------------------------------------!
 ! Description:
 ! -----------
@@ -689,5 +711,5 @@
 !> Cartesian coordinates (x,y) onto a sphere. It returns the latitude and
 !> lngitude of a geographical system centered at x0 and y0.
-!-----------------------------------------------------------------------------!
+!------------------------------------------------------------------------------!
  SUBROUTINE inv_plate_carree(x, y, x0, y0, r, lat, lon)
     REAL(wp), INTENT(IN)  ::  x(:), y(:), x0, y0, r
@@ -706,5 +728,5 @@
 
 
-!-----------------------------------------------------------------------------!
+!------------------------------------------------------------------------------!
 ! Description:
 ! ------------
