Ignore:
Timestamp:
Oct 26, 2020 10:05:58 AM (3 years ago)
Author:
eckhard
Message:

inifor: Fixed an error in surface pressure extrapolation

Location:
palm/trunk/UTIL/inifor/src
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/UTIL/inifor/src/inifor.f90

    r4675 r4756  
    2626! -----------------
    2727! $Id$
     28! Improved variable naming and code readabilty
     29!
     30!
     31! 4675 2020-09-11 10:00:26Z eckhard
    2832! Support for soil profile initialization
    2933! Improved code formatting
     
    157161
    158162    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)     ::  output_arr !< array buffer for interpolated quantities
    159     REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_centre !< density profile of the centre averaging domain
     163    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_centre_cosmo !< density profile of the centre averaging domain
    160164    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  ug_cosmo   !< profile of the geostrophic wind in x direction on COSMO levels
    161165    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  vg_cosmo   !< profile of the geostrophic wind in y direction on COSMO levels
    162166    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  ug_palm    !< profile of the geostrophic wind in x direction interpolated onto PALM levels
    163167    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  vg_palm    !< profile of the geostrophic wind in y direction interpolated onto PALM levels
    164     REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_north  !< density profile of the northern averaging domain
    165     REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_south  !< density profile of the southern averaging domain
    166     REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_east   !< density profile of the eastern averaging domain
    167     REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_west   !< density profile of the western averaging domain
    168     REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_north    !< pressure profile of the northern averaging domain
    169     REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_south    !< pressure profile of the southern averaging domain
    170     REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_east     !< pressure profile of the eastern averaging domain
    171     REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_west     !< pressure profile of the western averaging domain
    172 
    173     REAL(wp), POINTER, DIMENSION(:) ::  internal_arr !< pointer to the currently processed internal array (density, pressure)
     168    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_north_cosmo  !< density profile of the northern averaging domain
     169    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_south_cosmo  !< density profile of the southern averaging domain
     170    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_east_cosmo   !< density profile of the eastern averaging domain
     171    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_west_cosmo   !< density profile of the western averaging domain
     172    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_north_cosmo    !< pressure profile of the northern averaging domain
     173    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_south_cosmo    !< pressure profile of the southern averaging domain
     174    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_east_cosmo     !< pressure profile of the eastern averaging domain
     175    REAL(wp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_west_cosmo     !< pressure profile of the western averaging domain
     176
     177    REAL(wp), POINTER, DIMENSION(:) ::  internal_cosmo_arr !< pointer to the currently processed internal array (density, pressure)
    174178    REAL(wp), POINTER, DIMENSION(:) ::  ug_vg_palm   !< pointer to the currently processed geostrophic wind component
    175179
     
    334338                            ELSE
    335339                               CALL get_surface_pressure(                      &
    336                                   output_arr(0,0,:), rho_centre,               &
     340                                  output_arr(0,0,:),                           &
     341                                  rho_centre_cosmo,                            &
    337342                                  output_var%averaging_grid )
    338343                            ENDIF
     
    362367   
    363368                         CASE( 'internal_density_centre' )
    364                             ALLOCATE( rho_centre(1:cosmo_grid%nz) )
    365                             internal_arr => rho_centre
     369                            ALLOCATE( rho_centre_cosmo(1:cosmo_grid%nz) )
     370                            internal_cosmo_arr => rho_centre_cosmo
    366371   
    367372                         CASE( 'internal_density_north' )
    368                             ALLOCATE( rho_north(1:cosmo_grid%nz) )
    369                             internal_arr => rho_north
     373                            ALLOCATE( rho_north_cosmo(1:cosmo_grid%nz) )
     374                            internal_cosmo_arr => rho_north_cosmo
    370375   
    371376                         CASE( 'internal_density_south' )
    372                             ALLOCATE( rho_south(1:cosmo_grid%nz) )
    373                             internal_arr => rho_south
     377                            ALLOCATE( rho_south_cosmo(1:cosmo_grid%nz) )
     378                            internal_cosmo_arr => rho_south_cosmo
    374379   
    375380                         CASE( 'internal_density_east' )
    376                             ALLOCATE( rho_east(1:cosmo_grid%nz) )
    377                             internal_arr => rho_east
     381                            ALLOCATE( rho_east_cosmo(1:cosmo_grid%nz) )
     382                            internal_cosmo_arr => rho_east_cosmo
    378383   
    379384                         CASE( 'internal_density_west' )
    380                             ALLOCATE( rho_west(1:cosmo_grid%nz) )
    381                             internal_arr => rho_west
     385                            ALLOCATE( rho_west_cosmo(1:cosmo_grid%nz) )
     386                            internal_cosmo_arr => rho_west_cosmo
    382387   
    383388                         CASE( 'internal_pressure_north' )
    384                             ALLOCATE( p_north(1:cosmo_grid%nz) )
    385                             internal_arr => p_north
     389                            ALLOCATE( p_north_cosmo(1:cosmo_grid%nz) )
     390                            internal_cosmo_arr => p_north_cosmo
    386391   
    387392                         CASE( 'internal_pressure_south' )
    388                             ALLOCATE( p_south(1:cosmo_grid%nz) )
    389                             internal_arr => p_south
     393                            ALLOCATE( p_south_cosmo(1:cosmo_grid%nz) )
     394                            internal_cosmo_arr => p_south_cosmo
    390395   
    391396                         CASE( 'internal_pressure_east' )
    392                             ALLOCATE( p_east(1:cosmo_grid%nz) )
    393                             internal_arr => p_east
     397                            ALLOCATE( p_east_cosmo(1:cosmo_grid%nz) )
     398                            internal_cosmo_arr => p_east_cosmo
    394399   
    395400                         CASE( 'internal_pressure_west' )
    396                             ALLOCATE( p_west(1:cosmo_grid%nz) )
    397                             internal_arr => p_west
     401                            ALLOCATE( p_west_cosmo(1:cosmo_grid%nz) )
     402                            internal_cosmo_arr => p_west_cosmo
    398403   
    399404                         CASE DEFAULT
     
    413418                            CALL average_pressure_perturbation(                &
    414419                               input_buffer(output_var%input_id)%array(:,:,:), &
    415                                internal_arr(:),                                &
     420                               internal_cosmo_arr(:),                          &
    416421                               cosmo_grid, output_var%averaging_grid           &
    417422                            )
     
    421426                            CALL average_profile(                              &
    422427                               input_buffer(output_var%input_id)%array(:,:,:), &
    423                                internal_arr(:),                                &
     428                               internal_cosmo_arr(:),                          &
    424429                               output_var%averaging_grid                       &
    425430                            )
     
    445450                            vg_palm = cfg%vg
    446451                         ELSE
    447                             CALL geostrophic_winds( p_north, p_south, p_east,  &
    448                                                     p_west, rho_centre, f3,    &
    449                                                     averaging_width_ew,        &
    450                                                     averaging_width_ns,        &
    451                                                     phi_n, lambda_n,           &
    452                                                     phi_centre, lam_centre,    &
    453                                                     ug_cosmo, vg_cosmo )
     452                            CALL geostrophic_winds(                            &
     453                               p_north_cosmo, p_south_cosmo, p_east_cosmo,     &
     454                               p_west_cosmo, rho_centre_cosmo, f3,             &
     455                               averaging_width_ew,                             &
     456                               averaging_width_ns,                             &
     457                               phi_n, lambda_n,                                &
     458                               phi_centre, lam_centre,                         &
     459                               ug_cosmo, vg_cosmo                              &
     460                            )
    454461
    455462                            CALL interpolate_1d( ug_cosmo, ug_palm,            &
     
    532539             ug_vg_have_been_computed = .FALSE.
    533540             IF ( group%kind == 'thermodynamics' )  THEN
    534                 DEALLOCATE( rho_centre )
     541                DEALLOCATE( rho_centre_cosmo )
    535542                DEALLOCATE( ug_palm )
    536543                DEALLOCATE( vg_palm )
     
    538545                DEALLOCATE( vg_cosmo )
    539546                IF ( .NOT. cfg%ug_defined_by_user )  THEN
    540                    DEALLOCATE( rho_north )
    541                    DEALLOCATE( rho_south )
    542                    DEALLOCATE( rho_east )
    543                    DEALLOCATE( rho_west )
    544                    DEALLOCATE( p_north )
    545                    DEALLOCATE( p_south )
    546                    DEALLOCATE( p_east )
    547                    DEALLOCATE( p_west )
     547                   DEALLOCATE( rho_north_cosmo )
     548                   DEALLOCATE( rho_south_cosmo )
     549                   DEALLOCATE( rho_east_cosmo )
     550                   DEALLOCATE( rho_west_cosmo )
     551                   DEALLOCATE( p_north_cosmo )
     552                   DEALLOCATE( p_south_cosmo )
     553                   DEALLOCATE( p_east_cosmo )
     554                   DEALLOCATE( p_west_cosmo )
    548555                ENDIF
    549556             ENDIF
  • palm/trunk/UTIL/inifor/src/inifor_defs.f90

    r4714 r4756  
    214214    ACHAR( 10 ) // ' Copyright 2017-2020 Deutscher Wetterdienst Offenbach' !< Copyright notice
    215215 CHARACTER(LEN=*), PARAMETER ::  LOG_FILE_NAME = 'inifor.log' !< Name of INIFOR's log file
    216  CHARACTER(LEN=*), PARAMETER ::  VERSION = '2.1.1'            !< INIFOR version number
     216 CHARACTER(LEN=*), PARAMETER ::  VERSION = '2.1.2'            !< INIFOR version number
    217217 
    218218 END MODULE inifor_defs
  • palm/trunk/UTIL/inifor/src/inifor_transform.f90

    r4714 r4756  
    2626! -----------------
    2727! $Id$
     28! Fixed an error in surface pressure extrapolation where the cosmo grid was
     29!    misinterpreted as palm grid
     30! Improved code readability and formatting
     31!
     32!
     33! 4714 2020-09-29 12:47:35Z eckhard
    2834! Fixed off-by-one indexing error for profile quantities
    2935!
     
    381387!------------------------------------------------------------------------------!
    382388 SUBROUTINE interp_average_profile(source_array, profile_array, avg_grid)
    383     TYPE(grid_definition), INTENT(IN)          ::  avg_grid
    384     REAL(wp), DIMENSION(:,:,:), INTENT(IN)     ::  source_array(0:,0:,:)
    385     REAL(wp), DIMENSION(:), INTENT(OUT)        ::  profile_array
     389    TYPE(grid_definition), INTENT(IN)      ::  avg_grid
     390    REAL(wp), DIMENSION(:,:,:), INTENT(IN) ::  source_array(0:,0:,:)
     391    REAL(wp), DIMENSION(:), INTENT(OUT)    ::  profile_array
    386392
    387393    INTEGER(iwp) ::  i_source, j_source, k_profile, k_source, l, m
     
    438444 SUBROUTINE average_profile( source_array, profile_array, avg_grid )
    439445
    440     TYPE(grid_definition), INTENT(IN)          ::  avg_grid
    441     REAL(wp), DIMENSION(:,:,:), INTENT(IN)     ::  source_array(0:,0:,:)
    442     REAL(wp), DIMENSION(:), INTENT(OUT)        ::  profile_array
     446    TYPE(grid_definition), INTENT(IN) ::  avg_grid
     447    REAL(wp), INTENT(IN)              ::  source_array(0:,0:,:)
     448    REAL(wp), INTENT(OUT)             ::  profile_array(:)
    443449
    444450    INTEGER(iwp) ::  i_source, j_source, l, nz, nlev
     
    484490                                           cosmo_grid, avg_grid )
    485491
    486     TYPE(grid_definition), INTENT(IN)          ::  cosmo_grid, avg_grid
    487     REAL(wp), DIMENSION(:,:,:), INTENT(IN)     ::  cosmo_pressure(0:,0:,:)
    488     REAL(wp), DIMENSION(:), INTENT(OUT)        ::  profile_array
    489 
    490     INTEGER(iwp) ::  i_source, j_source, l, nz, nlev
    491 
    492     REAL(wp)                            ::  ni_columns
    493     REAL(wp), DIMENSION(:), ALLOCATABLE ::  basic_state_pressure
     492    TYPE(grid_definition), INTENT(IN) ::  cosmo_grid, avg_grid
     493    REAL(wp), INTENT(IN)              ::  cosmo_pressure(0:,0:,:)
     494    REAL(wp), INTENT(OUT)             ::  profile_array(:)
     495
     496    INTEGER(iwp)          ::  i_source, j_source, l, nz, nlev
     497    REAL(wp)              ::  ni_columns
     498    REAL(wp), ALLOCATABLE ::  basic_state_pressure(:)
    494499
    495500    nlev = SIZE( cosmo_pressure, 3 )
     
    624629! Description:
    625630! ------------
    626 !> Takes the averaged pressure profile <p> and sets the lowest entry to the
    627 !> extrapolated pressure at the surface.
    628 !------------------------------------------------------------------------------!
    629  SUBROUTINE get_surface_pressure(p, rho, avg_grid)
    630     REAL(wp), DIMENSION(:), INTENT(IN)    ::  rho
    631     REAL(wp), DIMENSION(:), INTENT(INOUT) ::  p
     631!> Takes the averaged pressure profile p_palm and sets the lowest entry to the
     632!> extrapolated pressure at the surface, assuming a linear density profile.
     633!------------------------------------------------------------------------------!
     634 SUBROUTINE get_surface_pressure(p_palm, rho_cosmo, avg_grid)
     635    REAL(wp), DIMENSION(:), INTENT(IN)    ::  rho_cosmo
     636    REAL(wp), DIMENSION(:), INTENT(INOUT) ::  p_palm
    632637    TYPE(grid_definition), INTENT(IN)     ::  avg_grid
    633638
    634     REAL(wp)     ::  drhodz, dz, zk, rhok
    635     INTEGER(iwp) ::  k_min
    636 
    637     k_min = avg_grid%k_min
    638     zk    = avg_grid%z(k_min)
    639     rhok  = rho(k_min)
    640     dz    = avg_grid%z(k_min + 1) - avg_grid%z(k_min)
    641     drhodz = ( rho(k_min + 1) - rho(k_min) ) / dz
    642 
    643     p(1) = linear_density_pressure( p(k_min), zk, rhok, drhodz,                &
    644                                     z = 0.0_wp, g=G )
     639    REAL(wp)     ::  drhodz_surface_cosmo
     640    INTEGER(iwp) ::  k_min_palm
     641
     642    drhodz_surface_cosmo =                                                     &
     643       ( rho_cosmo(2) - rho_cosmo(1) ) /                                       &
     644       ( avg_grid%intermediate_h(1,1,2) - avg_grid%intermediate_h(1,1,1) )
     645
     646    k_min_palm = avg_grid%k_min
     647
     648    p_palm(1) = linear_density_pressure(                                       &
     649       p0 = p_palm(k_min_palm),                                                &
     650       z0 = avg_grid%z(k_min_palm),                                            &
     651       rho00 = rho_cosmo(1),                                                   &
     652       z00 = avg_grid%intermediate_h(1,1,1),                                   &
     653       drhodz = drhodz_surface_cosmo,                                          &
     654       g = G,                                                                  &
     655       z = 0.0_wp                                                              &
     656    )
    645657
    646658 END SUBROUTINE get_surface_pressure
    647659
    648660
    649  FUNCTION linear_density_pressure(pk, zk, rhok, drhodz, z, g)  RESULT(p)
    650 
    651     REAL(wp), INTENT(IN)  ::  pk, zk, rhok, drhodz, g, z
     661!------------------------------------------------------------------------------!
     662! Description:
     663! -----------
     664!> Computes the hydrostatic pressure p at height z given the pressure p0 at
     665!> height z0. The pressure is computed based on the solution of the hydrostatic
     666!> equation assuming a linear density profile with density rho00 at z00 and the
     667!> vertical density gradient drhodz.
     668!------------------------------------------------------------------------------!
     669 FUNCTION linear_density_pressure(p0, z0, rho00, z00, drhodz, g, z)  RESULT(p)
     670
     671    REAL(wp), INTENT(IN)  ::  p0, z0, rho00, z00, drhodz, g, z
    652672    REAL(wp) ::  p
    653673
    654     p = pk + ( zk - z ) * g * ( rhok - 0.5_wp * drhodz * (zk - z) )
     674    p = p0 - ( z - z0 ) * g * (                                                &
     675           rho00 + 0.5_wp * drhodz * ( z + z0 - 2.0_wp * z00 )                 &
     676        )
    655677
    656678 END FUNCTION linear_density_pressure
    657679
    658 !-----------------------------------------------------------------------------!
     680!------------------------------------------------------------------------------!
    659681! Description:
    660682! -----------
    661683!> This routine computes profiles of the two geostrophic wind components ug and
    662684!> vg.
    663 !-----------------------------------------------------------------------------!
     685!------------------------------------------------------------------------------!
    664686 SUBROUTINE geostrophic_winds(p_north, p_south, p_east, p_west, rho, f3,    &
    665687                              Lx, Ly, phi_n, lam_n, phi_g, lam_g, ug, vg)
     
    683705
    684706
    685 !-----------------------------------------------------------------------------!
     707!------------------------------------------------------------------------------!
    686708! Description:
    687709! -----------
     
    689711!> Cartesian coordinates (x,y) onto a sphere. It returns the latitude and
    690712!> lngitude of a geographical system centered at x0 and y0.
    691 !-----------------------------------------------------------------------------!
     713!------------------------------------------------------------------------------!
    692714 SUBROUTINE inv_plate_carree(x, y, x0, y0, r, lat, lon)
    693715    REAL(wp), INTENT(IN)  ::  x(:), y(:), x0, y0, r
     
    706728
    707729
    708 !-----------------------------------------------------------------------------!
     730!------------------------------------------------------------------------------!
    709731! Description:
    710732! ------------
Note: See TracChangeset for help on using the changeset viewer.