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

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

Location:
palm/trunk/UTIL/inifor
Files:
9 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/UTIL/inifor/README

    r3764 r3779  
    1 # INIFOR - Mesoscale Interface for Initializing and Forcing PALM-4U (v1.4.6)
     1# INIFOR - Mesoscale Interface for Initializing and Forcing PALM-4U (v1.4.7)
    22
    33INIFOR provides the meteorological fields required to initialize and drive the
  • palm/trunk/UTIL/inifor/src/inifor.f90

    r3680 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! Average geostrophic wind components on coarse COSMO levels instead of fine PALM levels
     29! Remove --debug netCDF output of internal pressure profiles
     30!
     31! 3680 2019-01-18 14:54:12Z knoop
    2832! Prefixed all INIFOR modules with inifor_
    2933!
     
    7478!> and forcing data for the urban climate model PALM-4U. The required
    7579!> meteorological fields are interpolated from output data of the mesoscale
    76 !> model COSMO-DE. This is the main program file.
     80!> model COSMO. This is the main program file.
    7781!------------------------------------------------------------------------------!
    7882 PROGRAM inifor
     
    9195    USE inifor_io
    9296    USE inifor_transform,                                                      &
    93         ONLY:  average_profile, interpolate_2d, interpolate_3d,                &
    94                geostrophic_winds, extrapolate_density, extrapolate_pressure,   &
    95                get_surface_pressure
     97        ONLY:  average_pressure_perturbation, average_profile, interpolate_1d, &
     98               interpolate_1d_arr, interpolate_2d, interpolate_3d,             &
     99               interp_average_profile, geostrophic_winds, extrapolate_density, &
     100               extrapolate_pressure, get_surface_pressure
    96101    USE inifor_types
    97102   
     
    104109    REAL(dp), ALLOCATABLE, DIMENSION(:,:,:)     ::  output_arr !< array buffer for interpolated quantities
    105110    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_centre !< density profile of the centre averaging domain
    106     REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  ug_arr     !< geostrophic wind in x direction
    107     REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  vg_arr     !< geostrophic wind in y direction
     111    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  ug_cosmo   !< profile of the geostrophic wind in x direction on COSMO levels
     112    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  vg_cosmo   !< profile of the geostrophic wind in y direction on COSMO levels
     113    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  ug_palm    !< profile of the geostrophic wind in x direction interpolated onto PALM levels
     114    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  vg_palm    !< profile of the geostrophic wind in y direction interpolated onto PALM levels
    108115    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_north  !< density profile of the northern averaging domain
    109116    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_south  !< density profile of the southern averaging domain
     
    116123
    117124    REAL(dp), POINTER, DIMENSION(:) ::  internal_arr !< pointer to the currently processed internal array (density, pressure)
    118     REAL(dp), POINTER, DIMENSION(:) ::  ug_vg_arr    !< pointer to the currently processed geostrophic wind component
     125    REAL(dp), POINTER, DIMENSION(:) ::  ug_vg_palm   !< pointer to the currently processed geostrophic wind component
    119126
    120127    TYPE(nc_var), POINTER        ::  output_var      !< pointer to the currently processed output variable
     
    255262 CALL run_control('time', 'alloc')
    256263                     
    257                       CALL average_profile(                                    &
     264                      CALL interp_average_profile(                             &
    258265                         input_buffer(output_var % input_id) % array(:,:,:),   &
    259266                         output_arr(0,0,:),                                    &
     
    282289
    283290                      CASE('internal_density_centre')
    284                          ALLOCATE( rho_centre( 1:output_var % grid % nz ) )
     291                         ALLOCATE( rho_centre( 1:cosmo_grid % nz) )
    285292                         internal_arr => rho_centre
    286293
    287294                      CASE('internal_density_north')
    288                          ALLOCATE( rho_north( 1:output_var % grid % nz ) )
     295                         ALLOCATE( rho_north( 1:cosmo_grid % nz) )
    289296                         internal_arr => rho_north
    290297
    291298                      CASE('internal_density_south')
    292                          ALLOCATE( rho_south( 1:output_var % grid % nz ) )
     299                         ALLOCATE( rho_south( 1:cosmo_grid % nz) )
    293300                         internal_arr => rho_south
    294301
    295302                      CASE('internal_density_east')
    296                          ALLOCATE( rho_east( 1:output_var % grid % nz) )
     303                         ALLOCATE( rho_east( 1:cosmo_grid % nz) )
    297304                         internal_arr => rho_east
    298305
    299306                      CASE('internal_density_west')
    300                          ALLOCATE( rho_west( 1:output_var % grid % nz ) )
     307                         ALLOCATE( rho_west( 1:cosmo_grid % nz) )
    301308                         internal_arr => rho_west
    302309
    303310                      CASE('internal_pressure_north')
    304                          ALLOCATE( p_north( 1:output_var % grid % nz ) )
     311                         ALLOCATE( p_north( 1:cosmo_grid % nz) )
    305312                         internal_arr => p_north
    306313
    307314                      CASE('internal_pressure_south')
    308                          ALLOCATE( p_south( 1:output_var % grid % nz ) )
     315                         ALLOCATE( p_south( 1:cosmo_grid % nz) )
    309316                         internal_arr => p_south
    310317
    311318                      CASE('internal_pressure_east')
    312                          ALLOCATE( p_east( 1:output_var % grid % nz) )
     319                         ALLOCATE( p_east( 1:cosmo_grid % nz) )
    313320                         internal_arr => p_east
    314321
    315322                      CASE('internal_pressure_west')
    316                          ALLOCATE( p_west( 1:output_var % grid % nz ) )
     323                         ALLOCATE( p_west( 1:cosmo_grid % nz) )
    317324                         internal_arr => p_west
    318325
     
    324331
    325332
    326                       CALL average_profile(                                 &
    327                          input_buffer(output_var % input_id) % array(:,:,:),&
    328                          internal_arr(:),                                   &
    329                          output_var % averaging_grid)
    330 
    331                       SELECT CASE (TRIM(output_var % name))
    332 
    333                       CASE('internal_density_centre',                          &
    334                            'internal_density_north',                           &
    335                            'internal_density_south',                           &
    336                            'internal_density_east',                            &
    337                            'internal_density_west')
    338                          CALL extrapolate_density(internal_arr,                &
    339                                                   output_var % averaging_grid)
    340 
    341                       CASE('internal_pressure_north')
    342                          CALL extrapolate_pressure(internal_arr, rho_north,    &
    343                                                    output_var % averaging_grid)
    344 
    345                       CASE('internal_pressure_south')
    346                          CALL extrapolate_pressure(internal_arr, rho_south,    &
    347                                                    output_var % averaging_grid)
    348 
    349                       CASE('internal_pressure_east')
    350                          CALL extrapolate_pressure(internal_arr, rho_east,     &
    351                                                    output_var % averaging_grid)
    352 
    353                       CASE('internal_pressure_west')
    354                          CALL extrapolate_pressure(internal_arr, rho_west,     &
    355                                                    output_var % averaging_grid)
     333                      SELECT CASE( TRIM( output_var % name ) )
     334
     335                      CASE( 'internal_pressure_north',                         &
     336                            'internal_pressure_south',                         &
     337                            'internal_pressure_east',                          &
     338                            'internal_pressure_west' )
     339
     340                         CALL average_pressure_perturbation(                   &
     341                            input_buffer(output_var % input_id) % array(:,:,:),&
     342                            internal_arr(:),                                   &
     343                            cosmo_grid, output_var % averaging_grid            &
     344                         )
    356345
    357346                      CASE DEFAULT
    358                          CALL inifor_abort('main loop', message)
     347
     348                         CALL average_profile(                                 &
     349                            input_buffer(output_var % input_id) % array(:,:,:),&
     350                            internal_arr(:),                                   &
     351                            cosmo_grid, output_var % averaging_grid            &
     352                         )
    359353
    360354                      END SELECT
    361355
    362                       IF (.TRUE.)  THEN
    363                          ALLOCATE( output_arr(1,1,1:output_var % grid % nz) )
    364                          output_arr(1,1,:) = internal_arr(:)
    365                       END IF
     356
     357!
     358!--                   Output of geostrophic pressure profiles (with --debug
     359!--                   option) is currently deactivated, since they are now
     360!--                   defined on averaged COSMO levels instead of PALM levels
     361!--                   (requires definiton of COSMO levels in netCDF output.)
     362                      !IF (.TRUE.)  THEN
     363                      !   ALLOCATE( output_arr(1,1,1:output_var % grid % nz) )
     364                      !   output_arr(1,1,:) = internal_arr(:)
     365                      !END IF
    366366 CALL run_control('time', 'comp')
    367367
     
    369369!--                This case gets called twice, the first time for ug, the
    370370!--                second time for vg. We compute ug and vg at the first call
    371 !--                and keep vg (and ug for that matter) around for the second
    372 !--                call.
     371!--                and keep both of them around for the second call.
    373372                   CASE ( 'geostrophic winds' )
    374373
    375374                      IF (.NOT. ug_vg_have_been_computed )  THEN
    376                          ALLOCATE( ug_arr(output_var % grid % nz) )
    377                          ALLOCATE( vg_arr(output_var % grid % nz) )
    378 
    379                          IF ( cfg % ug_is_set )  THEN
    380                             ug_arr = cfg % ug
    381                             vg_arr = cfg % vg
     375                         ALLOCATE( ug_palm(output_var % grid % nz) )
     376                         ALLOCATE( vg_palm(output_var % grid % nz) )
     377                         ALLOCATE( ug_cosmo(cosmo_grid % nz) )
     378                         ALLOCATE( vg_cosmo(cosmo_grid % nz) )
     379
     380                         IF ( cfg % ug_defined_by_user )  THEN
     381                            ug_palm = cfg % ug
     382                            vg_palm = cfg % vg
    382383                         ELSE
    383384                            CALL geostrophic_winds( p_north, p_south, p_east,  &
     
    387388                                                    phi_n, lambda_n,           &
    388389                                                    phi_centre, lam_centre,    &
    389                                                     ug_arr, vg_arr )
     390                                                    ug_cosmo, vg_cosmo )
     391
     392                            CALL interpolate_1d( ug_cosmo, ug_palm,             &
     393                                                 output_var % grid )
     394
     395                            CALL interpolate_1d( vg_cosmo, vg_palm,             &
     396                                                 output_var % grid )
    390397                         END IF
    391398
     
    395402
    396403!
    397 !--                   Prepare output of geostrophic winds
     404!--                   Select output array of current geostrophic wind component
    398405                      SELECT CASE(TRIM(output_var % name))
    399406                      CASE ('ls_forcing_ug')
    400                          ug_vg_arr => ug_arr
     407                         ug_vg_palm => ug_palm
    401408                      CASE ('ls_forcing_vg')
    402                          ug_vg_arr => vg_arr
     409                         ug_vg_palm => vg_palm
    403410                      END SELECT
    404411
    405412                      ALLOCATE( output_arr(1,1,output_var % grid % nz) )
    406                       output_arr(1,1,:) = ug_vg_arr(:)
     413                      output_arr(1,1,:) = ug_vg_palm(:)
    407414
    408415                   CASE ( 'average scalar' )
     
    447454!- Section 2.3: Write current time step of current variable
    448455!------------------------------------------------------------------------------
    449                    IF (.NOT. output_var % is_internal .OR. debugging_output)  THEN
     456!
     457!--                Output of geostrophic pressure profiles (with --debug
     458!--                option) is currently deactivated, since they are now
     459!--                defined on averaged COSMO levels instead of PALM levels
     460!--                (requires definiton of COSMO levels in netCDF output.)
     461                   !IF (.NOT. output_var % is_internal .OR. debugging_output)  THEN
     462
     463                   IF (.NOT. output_var % is_internal)  THEN
    450464                      message = "Writing variable '" // TRIM(output_var%name) // "'."
    451465                      CALL report('main loop', message)
     
    467481             IF ( group % kind == 'thermodynamics' )  THEN
    468482                DEALLOCATE( rho_centre )
    469                 DEALLOCATE( ug_arr )
    470                 DEALLOCATE( vg_arr )
    471                 IF ( .NOT. cfg % ug_is_set )  THEN
     483                DEALLOCATE( ug_palm )
     484                DEALLOCATE( vg_palm )
     485                DEALLOCATE( ug_cosmo )
     486                DEALLOCATE( vg_cosmo )
     487                IF ( .NOT. cfg % ug_defined_by_user )  THEN
    472488                   DEALLOCATE( rho_north )
    473489                   DEALLOCATE( rho_south )
  • palm/trunk/UTIL/inifor/src/inifor_control.f90

    r3680 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!
  • palm/trunk/UTIL/inifor/src/inifor_defs.f90

    r3764 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! Updated version number to 1.4.7, updated copyright note
     29!
     30!
     31! 3764 2019-02-26 13:42:09Z eckhard
    2832! Bumped version number
    2933!
     
    136140 INTEGER, PARAMETER          ::  FORCING_STEP = 1             !< Number of hours between forcing time steps [h]
    137141 REAL(dp), PARAMETER         ::  NUDGING_TAU = 21600.0_dp     !< Nudging relaxation time scale [s]
    138  CHARACTER(LEN=*), PARAMETER ::  VERSION = '1.4.6'            !< INIFOR version number
    139  CHARACTER(LEN=*), PARAMETER ::  COPYRIGHT = 'Copyright 2017-2018 Leibniz Universitaet Hannover' // &
    140      ACHAR( 10 ) // ' Copyright 2017-2018 Deutscher Wetterdienst Offenbach' !< Copyright notice
     142 CHARACTER(LEN=*), PARAMETER ::  VERSION = '1.4.7'            !< INIFOR version number
     143 CHARACTER(LEN=*), PARAMETER ::  COPYRIGHT = 'Copyright 2017-2019 Leibniz Universitaet Hannover' // &
     144     ACHAR( 10 ) // ' Copyright 2017-2019 Deutscher Wetterdienst Offenbach' !< Copyright notice
    141145
    142146 END MODULE inifor_defs
  • palm/trunk/UTIL/inifor/src/inifor_grid.f90

    r3765 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! Assigned names to averaging grids
     29! Improved variable naming and minor clean-up
     30!
     31!
     32! 3765 2019-02-26 13:45:46Z eckhard
    2833! Removed dependency on radiation input files
    2934!
     
    411416!--    Overwrite defaults with user configuration
    412417       CALL parse_command_line_arguments( cfg )
     418       CALL report('main_loop', 'Running INIFOR version ' // VERSION)
    413419
    414420       flow_prefix = TRIM(cfg % input_prefix)
     
    444450       CALL validate_config( cfg )
    445451
    446        CALL report('main_loop', 'Running INIFOR version ' // VERSION)
    447452       CALL report('setup_parameters', "initialization mode: " // TRIM(cfg % ic_mode))
    448453       CALL report('setup_parameters', "       forcing mode: " // TRIM(cfg % bc_mode))
    449454       CALL report('setup_parameters', "     averaging mode: " // TRIM(cfg % averaging_mode))
     455       CALL report('setup_parameters', "    averaging angle: " // real_to_str(cfg % averaging_angle))
    450456       CALL report('setup_parameters', "          data path: " // TRIM(cfg % input_path))
    451457       CALL report('setup_parameters', "           hhl file: " // TRIM(cfg % hhl_file))
     
    10511057               lonmin = lonmin_palm, lonmax = lonmax_palm,                     &
    10521058               latmin = latmin_palm, latmax = latmax_palm,                     &
    1053                kind='scalar')
     1059               kind='scalar', name='averaged initial scalar')
    10541060
    10551061       CALL init_averaging_grid(averaged_initial_w_profile, cosmo_grid,        &
     
    10571063               lonmin = lonmin_palm, lonmax = lonmax_palm,                     &
    10581064               latmin = latmin_palm, latmax = latmax_palm,                     &
    1059                kind='w')
     1065               kind='w', name='averaged initial w')
    10601066
    10611067       CALL init_averaging_grid(averaged_scalar_profile, cosmo_grid,           &
     
    10631069               lonmin = lam_west, lonmax = lam_east,                           &
    10641070               latmin = phi_south, latmax = phi_north,                         &
    1065                kind='scalar')
     1071               kind='scalar', name='centre geostrophic scalar')
    10661072
    10671073       CALL init_averaging_grid(averaged_w_profile, cosmo_grid,                &
     
    10691075               lonmin = lam_west, lonmax = lam_east,                           &
    10701076               latmin = phi_south, latmax = phi_north,                         &
    1071                kind='w')
     1077               kind='w', name='centre geostrophic w')
    10721078
    10731079       CALL init_averaging_grid(south_averaged_scalar_profile, cosmo_grid,     &
     
    10761082               latmin = phi_centre - averaging_angle,                          &
    10771083               latmax = phi_centre,                                            &
    1078                kind='scalar')
     1084               kind='scalar', name='south geostrophic scalar')
    10791085
    10801086       CALL init_averaging_grid(north_averaged_scalar_profile, cosmo_grid,     &
     
    10831089               latmin = phi_centre,                                            &
    10841090               latmax = phi_centre + averaging_angle,                          &
    1085                kind='scalar')
     1091               kind='scalar', name='north geostrophic scalar')
    10861092
    10871093       CALL init_averaging_grid(west_averaged_scalar_profile, cosmo_grid,      &
     
    10901096               lonmax = lam_centre,                                            &
    10911097               latmin = phi_south, latmax = phi_north,                         &
    1092                kind='scalar')
     1098               kind='scalar', name='west geostrophic scalar')
    10931099
    10941100       CALL init_averaging_grid(east_averaged_scalar_profile, cosmo_grid,      &
     
    10971103               lonmax = lam_centre + averaging_angle,                          &
    10981104               latmin = phi_south, latmax = phi_north,                         &
    1099                kind='scalar')
     1105               kind='scalar', name='east geostrophic scalar')
    11001106
    11011107!                                                                             
     
    13041310        grid % z0 = z0
    13051311
    1306         SELECT CASE( TRIM (kind) )
     1312        SELECT CASE( TRIM(kind) )
    13071313
    13081314        CASE('boundary')
     
    15291535!------------------------------------------------------------------------------!
    15301536    SUBROUTINE init_averaging_grid(avg_grid, cosmo_grid, x, y, z, z0,          &
    1531        lonmin, lonmax, latmin, latmax, kind)
     1537       lonmin, lonmax, latmin, latmax, kind, name)
    15321538
    15331539       TYPE(grid_definition), INTENT(INOUT) ::  avg_grid
     
    15411547
    15421548       CHARACTER(LEN=*), INTENT(IN)         ::  kind
    1543 
    1544        LOGICAL                              :: level_based_averaging
     1549       CHARACTER(LEN=*), INTENT(IN)         ::  name
     1550
     1551       LOGICAL                              ::  level_based_averaging
    15451552
    15461553       ALLOCATE( avg_grid % x(1) )
     
    15591566
    15601567       avg_grid % kind = TRIM(kind)
     1568       avg_grid % name(1) = TRIM(name)
    15611569
    15621570!
     
    15861594!
    15871595!--    For level-besed averaging, compute average heights
     1596       !level_based_averaging = ( TRIM(mode) == 'level' )
    15881597       level_based_averaging = ( TRIM(cfg % averaging_mode) == 'level' )
    15891598       IF (level_based_averaging)  THEN
     
    15921601          CALL average_2d(avg_grid % cosmo_h, avg_grid % h(1,1,:),             &
    15931602                          avg_grid % iii, avg_grid % jjj)
     1603
    15941604       END IF
    15951605
     
    16451655       jmax = CEILING( (avg_grid % lat(2) - cosmo_lat(0)) / dlat )
    16461656       
    1647        message = "Averaging '" // TRIM(avg_grid % kind) // "' over "//        &
     1657       message = "Grid " // TRIM(avg_grid % name(1)) // " averages over " // &
    16481658                 TRIM(str(imin)) // " <= i <= " // TRIM(str(imax)) //          &
    16491659                 " and " //                                                    &
     
    31673177       )
    31683178       output_var_table(57) % averaging_grid => north_averaged_scalar_profile
    3169        output_var_table(57) % to_be_processed = .NOT. cfg % ug_is_set
     3179       output_var_table(57) % to_be_processed = .NOT. cfg % ug_defined_by_user
    31703180
    31713181
     
    31823192       )
    31833193       output_var_table(58) % averaging_grid => south_averaged_scalar_profile
    3184        output_var_table(58) % to_be_processed = .NOT. cfg % ug_is_set
     3194       output_var_table(58) % to_be_processed = .NOT. cfg % ug_defined_by_user
    31853195
    31863196
     
    31973207       )
    31983208       output_var_table(59) % averaging_grid => east_averaged_scalar_profile
    3199        output_var_table(59) % to_be_processed = .NOT. cfg % ug_is_set
     3209       output_var_table(59) % to_be_processed = .NOT. cfg % ug_defined_by_user
    32003210
    32013211
     
    32123222       )
    32133223       output_var_table(60) % averaging_grid => west_averaged_scalar_profile
    3214        output_var_table(60) % to_be_processed = .NOT. cfg % ug_is_set
     3224       output_var_table(60) % to_be_processed = .NOT. cfg % ug_defined_by_user
    32153225
    32163226       output_var_table(61) = init_nc_var(                                     &
     
    32263236       )
    32273237       output_var_table(61) % averaging_grid => north_averaged_scalar_profile
    3228        output_var_table(61) % to_be_processed = .NOT. cfg % ug_is_set
     3238       output_var_table(61) % to_be_processed = .NOT. cfg % ug_defined_by_user
    32293239
    32303240
     
    32413251       )
    32423252       output_var_table(62) % averaging_grid => south_averaged_scalar_profile
    3243        output_var_table(62) % to_be_processed = .NOT. cfg % ug_is_set
     3253       output_var_table(62) % to_be_processed = .NOT. cfg % ug_defined_by_user
    32443254
    32453255
     
    32563266       )
    32573267       output_var_table(63) % averaging_grid => east_averaged_scalar_profile
    3258        output_var_table(63) % to_be_processed = .NOT. cfg % ug_is_set
     3268       output_var_table(63) % to_be_processed = .NOT. cfg % ug_defined_by_user
    32593269
    32603270
     
    32713281       )
    32723282       output_var_table(64) % averaging_grid => west_averaged_scalar_profile
    3273        output_var_table(64) % to_be_processed = .NOT. cfg % ug_is_set
     3283       output_var_table(64) % to_be_processed = .NOT. cfg % ug_defined_by_user
    32743284
    32753285!
     
    34163426          var % task            = "interpolate_2d"
    34173427
    3418        CASE( 'left scalar', 'right scalar') ! same as right
     3428       CASE( 'left scalar', 'right scalar')
    34193429          var % lod             = -1
    34203430          var % ndim            = 3
     
    34293439          var % task            = "interpolate_3d"
    34303440
    3431        CASE( 'north scalar', 'south scalar') ! same as south
     3441       CASE( 'north scalar', 'south scalar')
    34323442          var % lod             = -1
    34333443          var % ndim            = 3
  • palm/trunk/UTIL/inifor/src/inifor_io.f90

    r3764 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! Temporariliy disabled height-based geostrophic wind averaging
     29! Improved variable naming
     30!
     31!
     32! 3764 2019-02-26 13:42:09Z eckhard
    2833! Removed dependency on radiation input files
    2934!
     
    397402
    398403       cfg % p0_is_set = .FALSE.
    399        cfg % ug_is_set = .FALSE.
    400        cfg % vg_is_set = .FALSE.
     404       cfg % ug_defined_by_user = .FALSE.
     405       cfg % vg_defined_by_user = .FALSE.
    401406       cfg % flow_prefix_is_set = .FALSE.
    402407       cfg % input_prefix_is_set = .FALSE.
     
    442447
    443448             CASE( '-ug', '-u', '--geostrophic-u' )
    444                 cfg % ug_is_set = .TRUE.
     449                cfg % ug_defined_by_user = .TRUE.
    445450                CALL get_option_argument( i, arg )
    446451                READ(arg, *) cfg % ug
    447452
    448453             CASE( '-vg', '-v', '--geostrophic-v' )
    449                 cfg % vg_is_set = .TRUE.
     454                cfg % vg_defined_by_user = .TRUE.
    450455                CALL get_option_argument( i, arg )
    451456                READ(arg, *) cfg % vg
     
    712717
    713718      SELECT CASE( TRIM(cfg % averaging_mode) )
    714       CASE( 'level', 'height')
     719      CASE( 'level' )
     720      CASE( 'height' )
     721         message = "Averaging mode '" // TRIM(cfg % averaging_mode) //&
     722                   "' is currently not supported. " //&
     723                   "Please use level-based averaging by selecting 'level', " //&
     724                   "or by omitting the --averaging-mode option entirely."
     725         CALL inifor_abort( 'validate_config', message )
    715726      CASE DEFAULT
    716727         message = "Averaging mode '" // TRIM(cfg % averaging_mode) //&
    717728                   "' is not supported. " //&
    718                    "Please select either 'height' or 'level', " //&
    719                    "or omit the --averaging-mode option entirely, which corresponds "//&
    720                    "to the latter."
     729         !          "Please select either 'height' or 'level', " //&
     730         !          "or omit the --averaging-mode option entirely, which corresponds "//&
     731         !          "to the latter."
     732                   "Please use level-based averaging by selecting 'level', " //&
     733                   "or by omitting the --averaging-mode option entirely."
    721734         CALL inifor_abort( 'validate_config', message )
    722735      END SELECT
    723736
    724       IF ( cfg % ug_is_set .NEQV. cfg % vg_is_set )  THEN
     737      IF ( cfg % ug_defined_by_user .NEQV. cfg % vg_defined_by_user )  THEN
    725738         message = "You specified only one component of the geostrophic " // &
    726739                   "wind. Please specify either both or none."
     
    965978          var => output_variable_table(i)
    966979
    967           to_be_written = ( var % to_be_processed  .AND. .NOT. var % is_internal) .OR.                        &
    968                           ( var % is_internal  .AND.  debug )
     980          !to_be_written = ( var % to_be_processed  .AND. .NOT. var % is_internal) .OR. &
     981          !                ( var % is_internal  .AND.  debug )
     982          to_be_written = ( var % to_be_processed  .AND. .NOT. var % is_internal)
    969983
    970984          IF ( to_be_written )  THEN
  • 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
  • palm/trunk/UTIL/inifor/src/inifor_types.f90

    r3680 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! Improved variable naming
     29!
     30! 3680 2019-01-18 14:54:12Z knoop
    2831! Prefixed all INIFOR modules with inifor_
    2932!
     
    109112    LOGICAL              ::  debug                       !< indicates whether --debug option was given
    110113    LOGICAL              ::  p0_is_set                   !< indicates whether p0 was set manually
    111     LOGICAL              ::  ug_is_set                   !< indicates whether ug was set manually
    112     LOGICAL              ::  vg_is_set                   !< indicates whether vg was set manually
    113     LOGICAL              ::  flow_prefix_is_set          !<  indicates whether the flow prefix was set manually
    114     LOGICAL              ::  input_prefix_is_set         !<  indicates whether the input prefix was set manually
    115     LOGICAL              ::  radiation_prefix_is_set     !<  indicates whether the radiation prefix was set manually
    116     LOGICAL              ::  soil_prefix_is_set          !<  indicates whether the soil prefix was set manually
    117     LOGICAL              ::  soilmoisture_prefix_is_set  !<  indicates whether the soilmoisture prefix was set manually
     114    LOGICAL              ::  ug_defined_by_user          !< indicates whether ug was set manually
     115    LOGICAL              ::  vg_defined_by_user          !< indicates whether vg was set manually
     116    LOGICAL              ::  flow_prefix_is_set          !< indicates whether the flow prefix was set manually
     117    LOGICAL              ::  input_prefix_is_set         !< indicates whether the input prefix was set manually
     118    LOGICAL              ::  radiation_prefix_is_set     !< indicates whether the radiation prefix was set manually
     119    LOGICAL              ::  soil_prefix_is_set          !< indicates whether the soil prefix was set manually
     120    LOGICAL              ::  soilmoisture_prefix_is_set  !< indicates whether the soilmoisture prefix was set manually
    118121 END TYPE inifor_config
    119122
  • palm/trunk/UTIL/inifor/src/inifor_util.f90

    r3680 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!
Note: See TracChangeset for help on using the changeset viewer.