Changeset 4756 for palm/trunk/UTIL/inifor/src/inifor_transform.f90

Timestamp:

Oct 26, 2020 10:05:58 AM (4 years ago)

Author:

eckhard

Message:

inifor: Fixed an error in surface pressure extrapolation

File:

• palm/trunk/UTIL/inifor/src/inifor_transform.f90 (modified) (8 diffs)

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

@@ -26 +26 @@
 ! -----------------
 ! $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 +387 @@
 !------------------------------------------------------------------------------!
  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 +444 @@
  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 +490 @@
                                            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 +629 @@
 ! 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 +705 @@
 
 
-!-----------------------------------------------------------------------------!
+!------------------------------------------------------------------------------!
 ! Description:
 ! -----------
@@ -689 +711 @@
 !> 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 +728 @@
 
 
-!-----------------------------------------------------------------------------!
+!------------------------------------------------------------------------------!
 ! Description:
 ! ------------