Ignore:
Timestamp:
Mar 5, 2019 11:13:35 AM (4 years ago)
Author:
eckhard
Message:

inifor: bugfix: Fixes issue #815 with geostrophic wind profiles

File:
1 edited

Legend:

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

    r3716 r3779  
    1515! PALM. If not, see <http://www.gnu.org/licenses/>.
    1616!
    17 ! Copyright 2017-2018 Leibniz Universitaet Hannover
    18 ! Copyright 2017-2018 Deutscher Wetterdienst Offenbach
     17! Copyright 2017-2019 Leibniz Universitaet Hannover
     18! Copyright 2017-2019 Deutscher Wetterdienst Offenbach
    1919!------------------------------------------------------------------------------!
    2020!
     
    2626! -----------------
    2727! $Id$
     28! Remove basic state pressure before computing geostrophic wind
     29!  - Introduced new level-based profile averaging routine that does not rely on
     30!    external weights average_profile()
     31!  - Renamed original weights-based routine average_profile() ->
     32!    interp_average_profile()
     33! Average geostrophic wind components on coarse COSMO levels instead of fine PALM levels
     34!  - Introduced new profile interpolation routine for interpolating single
     35!    profiles from COSMO to PALM levels
     36!  - Renamed original array variant interpolate_1d() -> interpolate_1d_arr()
     37!
     38!
     39!
     40! 3716 2019-02-05 17:02:38Z eckhard
    2841! Include out-of-bounds error message in log
    2942!
     
    93106    USE inifor_control
    94107    USE inifor_defs,                                                           &
    95         ONLY: G, TO_DEGREES, TO_RADIANS, PI, dp
     108        ONLY: BETA, dp, G, P_SL, PI, RD, T_SL, TO_DEGREES, TO_RADIANS
    96109    USE inifor_types
    97110    USE inifor_util,                                                           &
    98         ONLY: real_to_str, str
     111        ONLY: get_basic_state, real_to_str, str
    99112
    100113    IMPLICIT NONE
    101114
    102115 CONTAINS
     116
     117
     118    SUBROUTINE interpolate_1d(in_arr, out_arr, outgrid)
     119       TYPE(grid_definition), INTENT(IN) ::  outgrid
     120       REAL(dp), INTENT(IN)              ::  in_arr(:)
     121       REAL(dp), INTENT(OUT)             ::  out_arr(:)
     122
     123       INTEGER :: i, j, k, l, nz
     124
     125       nz = UBOUND(out_arr, 1)
     126
     127       DO k = nz, LBOUND(out_arr, 1), -1
     128
     129!
     130!--       TODO: Remove IF clause and extrapolate based on a critical vertical
     131!--       TODO: index marking the lower bound of COSMO-DE data coverage.
     132!--       Check for negative interpolation weights indicating grid points
     133!--       below COSMO-DE domain and extrapolate from the top in such cells.
     134          IF (outgrid % w(1,k,1) < -1.0_dp .AND. k < nz)  THEN
     135             out_arr(k) = out_arr(k+1)
     136          ELSE
     137             out_arr(k) = 0.0_dp
     138             DO l = 1, 2
     139                out_arr(k) = out_arr(k) +                                      &
     140                    outgrid % w(1,k,l) * in_arr(outgrid % kkk(1,k,l) )
     141             END DO
     142          END IF
     143       END DO
     144
     145    END SUBROUTINE interpolate_1d
     146
    103147
    104148!------------------------------------------------------------------------------!
     
    125169!> outvar : Array of interpolated data
    126170!------------------------------------------------------------------------------!
    127     SUBROUTINE interpolate_1d(in_arr, out_arr, outgrid)
     171    SUBROUTINE interpolate_1d_arr(in_arr, out_arr, outgrid)
    128172       TYPE(grid_definition), INTENT(IN) ::  outgrid
    129173       REAL(dp), INTENT(IN)              ::  in_arr(0:,0:,0:)
     
    156200       END DO
    157201       END DO
    158     END SUBROUTINE interpolate_1d
     202    END SUBROUTINE interpolate_1d_arr
    159203
    160204
     
    174218!>     of PALM-4U point (i,j,k) on the input grid corresponding to the source
    175219!>     data invar. (The outgrid carries the relationship with the ingrid in the
    176 !      form of the interpoaltion weights.)
     220!      form of the interpolation weights.)
    177221!>
    178222!> outgrid % w_horiz: Array of weights for horizontal bi-linear interpolation
     
    298342!--    Interpolate from intermediate grid to palm_grid grid, includes
    299343!--    extrapolation for cells below COSMO domain.
    300        CALL interpolate_1d(intermediate_array, palm_array, palm_grid)
     344       CALL interpolate_1d_arr(intermediate_array, palm_array, palm_grid)
    301345
    302346       DEALLOCATE(intermediate_array)
     
    311355!> averaging grid 'avg_grid' and store the result in 'profile_array'.
    312356!------------------------------------------------------------------------------!
    313     SUBROUTINE average_profile(source_array, profile_array, avg_grid)
     357    SUBROUTINE interp_average_profile(source_array, profile_array, avg_grid)
    314358       TYPE(grid_definition), INTENT(IN)          ::  avg_grid
    315359       REAL(dp), DIMENSION(:,:,:), INTENT(IN)     ::  source_array
     
    358402       profile_array(1:avg_grid % k_min-1) = profile_array(avg_grid % k_min)
    359403
     404    END SUBROUTINE interp_average_profile
     405
     406
     407!------------------------------------------------------------------------------!
     408! Description:
     409! ------------
     410!> Average data horizontally from the source_array over the region given by the
     411!> averaging grid 'avg_grid' and store the result in 'profile_array'.
     412!------------------------------------------------------------------------------!
     413    SUBROUTINE average_profile( source_array, profile_array,                   &
     414                                source_grid, avg_grid )
     415
     416       TYPE(grid_definition), INTENT(IN)          ::  source_grid, avg_grid
     417       REAL(dp), DIMENSION(:,:,:), INTENT(IN)     ::  source_array
     418       REAL(dp), DIMENSION(:), INTENT(OUT)        ::  profile_array
     419
     420       INTEGER ::  i_source, j_source, k_profile, k_source, l, m, nz, nlev
     421
     422       REAL                            ::  ni_columns
     423
     424       nlev = SIZE( source_array, 3 )
     425       nz   = SIZE( profile_array, 1 )
     426
     427       IF ( nlev /= nz )  THEN
     428          message = "Lengths of input and output profiles do not match: " //   &
     429                    "cosmo_pressure(" // TRIM( str( nlev ) ) //                &
     430                    "), profile_array(" // TRIM( str( nz ) )  // ")."
     431          CALL inifor_abort('average_pressure_perturbation', message)
     432       ENDIF
     433     
     434       profile_array(:) = 0.0_dp
     435
     436       DO l = 1, avg_grid % n_columns
     437
     438          i_source = avg_grid % iii(l)
     439          j_source = avg_grid % jjj(l)
     440
     441          profile_array(:) = profile_array(:)                                  &
     442                           + source_array(i_source, j_source, :)
     443
     444       END DO
     445
     446       ni_columns = 1.0_dp / avg_grid % n_columns
     447       profile_array(:) = profile_array(:) * ni_columns
     448
    360449    END SUBROUTINE average_profile
     450
     451
     452!------------------------------------------------------------------------------!
     453! Description:
     454! ------------
     455!> This is a sister routine to average_profile() and differes from it in that
     456!> it removes the COSMO basic state pressure from the input array before
     457!> averaging.
     458!------------------------------------------------------------------------------!
     459    SUBROUTINE average_pressure_perturbation( cosmo_pressure, profile_array,   &
     460                                              cosmo_grid, avg_grid )
     461
     462       TYPE(grid_definition), INTENT(IN)          ::  cosmo_grid, avg_grid
     463       REAL(dp), DIMENSION(:,:,:), INTENT(IN)     ::  cosmo_pressure
     464       REAL(dp), DIMENSION(:), INTENT(OUT)        ::  profile_array
     465
     466       INTEGER ::  i_source, j_source, k_profile, k_source, l, m, nz, nlev
     467
     468       REAL(dp)                            ::  ni_columns
     469       REAL(dp), DIMENSION(:), ALLOCATABLE ::  basic_state_pressure
     470
     471       nlev = SIZE( cosmo_pressure, 3 )
     472       nz   = SIZE( profile_array, 1 )
     473
     474       IF ( nlev /= nz )  THEN
     475          message = "Lengths of input and output profiles do not match: " //   &
     476                    "cosmo_pressure(" // TRIM( str( nlev ) ) //                &
     477                    "), profile_array(" // TRIM( str( nz ) )  // ")."
     478          CALL inifor_abort('average_pressure_perturbation', message)
     479       ENDIF
     480
     481       ALLOCATE( basic_state_pressure(nz) )
     482       profile_array(:) = 0.0_dp
     483
     484       DO l = 1, avg_grid % n_columns
     485          i_source = avg_grid % iii(l)
     486          j_source = avg_grid % jjj(l)
     487
     488!
     489!--       Compute pressure perturbation by removing COSMO basic state pressure
     490          CALL get_basic_state( cosmo_grid % hfl(i_source,j_source,:), BETA,   &
     491                                P_SL, T_SL, RD, G, basic_state_pressure )
     492
     493          profile_array(:) = profile_array(:)                                  &
     494                           + cosmo_pressure(i_source, j_source, :)             &
     495                           - basic_state_pressure(:)
     496
     497!
     498!--    Loop over horizontal neighbours l
     499       END DO
     500
     501       DEALLOCATE( basic_state_pressure )
     502
     503       ni_columns = 1.0_dp / avg_grid % n_columns
     504       profile_array(:) = profile_array(:) * ni_columns
     505
     506    END SUBROUTINE average_pressure_perturbation
     507
     508
    361509
    362510
Note: See TracChangeset for help on using the changeset viewer.