Changeset 3395 for palm/trunk/UTIL


Ignore:
Timestamp:
Oct 22, 2018 5:32:49 PM (6 years ago)
Author:
eckhard
Message:

inifor: Added computation of geostrophic winds from COSMO input

Location:
palm/trunk/UTIL/inifor
Files:
2 deleted
12 edited

Legend:

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

    r3183 r3395  
    2525# -----------------
    2626# $Id$
     27# Updated build order
     28#
     29# 3183 2018-07-27 14:25:55Z suehring
    2730# Added __netcdf4 preprocessor flag
    2831#
     
    4750TEST_PATH = $(PROJECT_PATH)/tests
    4851
    49 MODULES = $(SRC_PATH)/defs.mod $(SRC_PATH)/util.mod $(SRC_PATH)/control.mod $(SRC_PATH)/types.mod $(SRC_PATH)/transform.mod $(SRC_PATH)/io.mod $(SRC_PATH)/grid.mod
     52MODULES = $(SRC_PATH)/defs.mod $(SRC_PATH)/types.mod $(SRC_PATH)/util.mod $(SRC_PATH)/control.mod $(SRC_PATH)/transform.mod $(SRC_PATH)/io.mod $(SRC_PATH)/grid.mod
    5053SOURCES = $(MODULES:%.mod=%.f90) $(SRC_PATH)/$(PROJECT).f90
    5154OBJECTS = $(SOURCES:%.f90=%.o)
  • palm/trunk/UTIL/inifor/Makefile.gnu

    r3183 r3395  
    2525# -----------------
    2626# $Id$
     27# Updated build order
     28#
     29# 3183 2018-07-27 14:25:55Z suehring
    2730# Added __netcdf4 preprocessor flag
    2831# Corrected compilation order
     
    4447PROJECT = inifor
    4548PROJECT_PATH = .
    46 BIN_PATH  = $(PROJECT_PATH)/bin
     49BIN_PATH  = $(PROJECT_PATH)/../../SCRIPTS
    4750SRC_PATH  = $(PROJECT_PATH)/src
    4851TEST_PATH = $(PROJECT_PATH)/tests
    4952
    50 MODULES = $(SRC_PATH)/defs.mod $(SRC_PATH)/util.mod $(SRC_PATH)/control.mod $(SRC_PATH)/types.mod $(SRC_PATH)/transform.mod $(SRC_PATH)/io.mod $(SRC_PATH)/grid.mod
     53MODULES = $(SRC_PATH)/defs.mod $(SRC_PATH)/types.mod $(SRC_PATH)/util.mod $(SRC_PATH)/control.mod $(SRC_PATH)/transform.mod $(SRC_PATH)/io.mod $(SRC_PATH)/grid.mod
    5154SOURCES = $(MODULES:%.mod=%.f90) $(SRC_PATH)/$(PROJECT).f90
    5255OBJECTS = $(SOURCES:%.f90=%.o)
  • palm/trunk/UTIL/inifor/Makefile.ifort

    r3183 r3395  
    2525# -----------------
    2626# $Id$
     27# Updated build order
     28#
     29# 3183 2018-07-27 14:25:55Z suehring
    2730# Added __netcdf4 preprocessor flag
    2831# Corrected compilation order
     
    4447PROJECT = inifor
    4548PROJECT_PATH = .
    46 BIN_PATH  = $(PROJECT_PATH)/bin
     49BIN_PATH  = $(PROJECT_PATH)/../../SCRIPTS
    4750SRC_PATH  = $(PROJECT_PATH)/src
    4851TEST_PATH = $(PROJECT_PATH)/tests
    4952
    50 MODULES = $(SRC_PATH)/defs.mod $(SRC_PATH)/util.mod $(SRC_PATH)/control.mod $(SRC_PATH)/types.mod $(SRC_PATH)/transform.mod $(SRC_PATH)/io.mod $(SRC_PATH)/grid.mod
     53MODULES = $(SRC_PATH)/defs.mod $(SRC_PATH)/types.mod $(SRC_PATH)/util.mod $(SRC_PATH)/control.mod $(SRC_PATH)/transform.mod $(SRC_PATH)/io.mod $(SRC_PATH)/grid.mod
    5154SOURCES = $(MODULES:%.mod=%.f90) $(SRC_PATH)/$(PROJECT).f90
    5255OBJECTS = $(SOURCES:%.f90=%.o)
  • palm/trunk/UTIL/inifor/README

    r3182 r3395  
    1 # INIFOR - Mesoscale Interface for Initializing and Forcing PALM-4U (v1.3.0)
     1# INIFOR - Mesoscale Interface for Initializing and Forcing PALM-4U (v1.4.0)
    22
    33INIFOR provides the meteorological fields required to initialize and drive the
     
    1414
    1515
     16## MINIMAL USAGE EXAMPLE
     17
     181. Customize `./namelist` (number of grid points and spacings, end_time)
     192. Run `current_version/trunk/SCRIPTS/inifor --path <scenario path> --date <YYYYMMDD>`
     20
     21
    1622## USAGE
    1723
    18 1. Customize `./namelist` (number or grid points and spacings, end_time)
    19 2. Run `current_version/trunk/SCRIPTS/inifor -path <scenario path> -date <YYYYMMDD>`
     24After compilation, the `inifor` binary resides in the `$PALM_BIN` path, i.e. in
     25`<PALM path>/current_version/trunk/SCRIPTS/`.
     26
     27In order to run, INIFOR requires three kinds of inputs:
     28
     291. hourly COSMO model output,
     302. a steering namelist file, and
     313. command-line options.
     32
     33In addition, a static driver file may be supplied (see `--static-driver` option
     34below) in order to pass the coordinates of the PALM origin to INIFOR. If no
     35static driver is passed to INIFOR, origin coordinates are read from the namelist
     36file.
     37
     38A typical `inifor` call looks like this:
     39
     40    inifor --path /data/evaluation/20170729 --date 2017073006 \
     41           --init-mode profile -n namelist -o dynamic_driver.nc \
     42           --elevation 110 --input-prefix lff0
     43
     44### INPUT DATA: COSMO MODEL OUTPUT
     45
     46INIFOR processes COSMO model output which it requires to be stored in a set of
     47netCDF files located in a user-specified path (see `--path` option). These are:
     48
     49- hhl.nc: This file provides the COSMO numerical grid. (*hhl* abbreviates
     50  *height of half layers*, i.e. the heights of the vertical cell boundaries.)
     51- soil.nc: This file provides the COSMO soil map and is used to destinguish land
     52  from water cells.
     53- `<prefix>YYYYMMDDHH-<suffix>.nc`: Each of these files contains COSMO model
     54  of one time step at the given time in UTC (Y..year, M..month, D..day, H..hour).
     55    - The `<prefix>` distinguishes different DWD products, for instance COSMO
     56      analyses (`laf`) or forecasts (`lff`).
     57    - The `<suffix>` distinguishes four kinds of COSMO model output data, namely
     58        - `flow` (atmospheric fields)
     59        - `rad` (radiation)
     60        - `soil` (soil moisture and temperature)
     61        - `soilmoisture` (precipitation and evaporation)
     62    - For example, laf2016010100-flow.nc contains the atmospheric fields of the
     63      COSMO analysis of Januray 1st, 2016 for 0:00 UTC.
    2064
    2165
    22 ## AVAILABLE NAMELIST PARAMETERS
     66### AVAILABLE NAMELIST PARAMETERS
    2367
    2468INIFOR mirrors a subset of the PALM-4U's Fortran namelists `inipar` and `d3par`
     
    2670
    2771
    28 ### inipar
     72#### inipar
    2973
    3074    nx, ny, nz - number of PALM-4U grid points in x, y, and z direction
     
    3478    dz_max - maximum vertical grid spacing [m]
    3579    dz_stretch_level_start(9) - array of height levels above which the grid is
    36     to be stretched vertically [m]
     80        to be stretched vertically [m]
    3781    dz_stretch_level_end(9) - array of height levels until which the grid is to
    38     be stretched vertically [m]
     82        be stretched vertically [m]
    3983    longitude, latitude - geographical coordinates of the PALM-4U origin [deg]
    4084
    4185
    42 ### d3par
     86#### d3par
    4387
    4488    end_time - PALM-4U simulation time. INIFOR will produce hourly forcing data
     
    4791
    4892
    49 ### EXAMPLE NAMELIST FILE
     93#### EXAMPLE NAMELIST FILE
    5094
    5195    &inipar nx = 4679, ny = 3939, nz = 360
     
    59103
    60104
    61 ## AVAILABLE COMMAND-LINE PARAMETERS
     105### AVAILABLE COMMAND-LINE PARAMETERS
     106    --averaging-mode <mode>:
     107        Selects how averaged quantities (large-scale forcing terms) are computed.
     108        INIFOR supports averaging along input model levels ('level') and along
     109        constant heights ('height'). Default: level
     110
     111    -a, --averaging-angle <angle>:
     112        Width of the averaging box in longitudal and latitudal direction in the
     113        source coordinate system (COSMO rotated-pole) [deg]. Default: 2.0
    62114
    63115    -d, --date, -date <date>:
     
    65117        hours (HH) are given, INIFOR assumes that the simulation starts at O UTC
    66118        on that day. Default: 20130721
    67        
     119
    68120    -i, --init-mode, -mode <mode>:
    69121        Set the PALM-4U initialization mode. INIFOR can provide initial conditions
     
    72124
    73125    -l, --hhl-file, -hhl <netCDF file>:
    74         Location of the netCDF file containing the vertical COSMO-DE grid levels,
     126        Location of the netCDF file containing the vertical COSMO grid levels,
    75127        specifically the heights of half levels (hhl, i.e. vertical cell faces).
    76128        Default: <scenario path>/hhl.nc
     
    83135    -o, --output <output file>:
    84136        Name of the INIFOR output file, i.e. the PALM-4U dynamic driver.
    85         Default: ./palm-4u-input.nc
     137        Default: ./dynamic_driver.nc
    86138
    87139    -p, --path, -path <scenario path>:
     
    89141
    90142    -r, --surface-pressure, -p0 <pressure>:
    91         Surface pressure at z=0 in the PALM-4U domain [Pa].
    92         Default: 1e5
     143        Manually set the pressure at z=0 in the PALM-4U domain [Pa]. If not
     144        given, the surface pressure is computed from the COSMO pressure field.
    93145
    94146    -s, --soil-file, -soil <netCDF file>:
    95         Location of the netCDF file containing the COSMO-DE soil type map.
     147        Location of the netCDF file containing the COSMO soil type map.
    96148        Default: <scenario path>/soil.nc
    97149
    98150    -t, --static-driver, -static <netCDF file>:
    99         Location of the netCDF file containing the static driver file for the case
     151        Location of the netCDF file containing the static driver for the case
    100152        to be simulated with PALM-4U. Optional parameter. Default: None
    101153
    102154    -u, --geostrophic-u, -ug <velocity>:
    103         Specifies the geostrophic wind in x direction [m/s]. Default: 0
     155        Manually specify the geostrophic wind in x direction [m/s]. If not
     156        given, the geostrophic wind is computed from the COSMO pressure field.
    104157
    105158    -v, --geostrophic-v, -vg <velocity>:
    106         Specifies the geostrophic wind in y direction [m/s]. Default: 0
    107 
    108     --version:
    109         Output version number and exit.
     159        Manually specify the geostrophic wind in y direction [m/s]. If not
     160        given, the geostrophic wind is computed from the COSMO pressure field.
    110161
    111162    -z, --elevation, -z0 <height>: Specifies the elevation of the PALM-4U domain
     
    113164
    114165
    115 ## ADDITIONAL COMMAND-LINE OPTIONS
     166#### ADDITIONAL COMMAND-LINE OPTIONS
     167
     168    --input-prefix <prefix>:
     169        Set the file prefixes for all input files. Individual prefixes can be
     170        overwritten with the options below. Default: laf
    116171
    117172    --flow-prefix <prefix>:
     
    126181    --soilmoisture-prefix <prefix>:
    127182        Set the file prefix of soil moisture input files. Default: laf
     183
     184    --debug:
     185        Enable debugging messages.
     186
     187    --version:
     188        Output version number and exit.
     189
  • palm/trunk/UTIL/inifor/src/control.f90

    r3183 r3395  
    2626! -----------------
    2727! $Id$
     28! Suppress debugging messages unless --debug option is given
     29!
     30!
     31! 3183 2018-07-27 14:25:55Z suehring
    2832! Added version and copyright output
    2933!
     
    5761 CONTAINS
    5862
    59     SUBROUTINE report(routine, message)
    60 
    61        CHARACTER(LEN=*), INTENT(IN) ::  routine
    62        CHARACTER(LEN=*), INTENT(IN) ::  message
    63        INTEGER                      ::  u
    64        LOGICAL, SAVE                ::  is_first_run = .TRUE.
    65 
    66        PRINT *, "inifor: " // TRIM(message) // "  [ " // TRIM(routine) // " ]"
     63    SUBROUTINE report(routine, message, debug)
     64
     65       CHARACTER(LEN=*), INTENT(IN)  ::  routine
     66       CHARACTER(LEN=*), INTENT(IN)  ::  message
     67       LOGICAL, OPTIONAL, INTENT(IN) ::  debug
     68       INTEGER                       ::  u
     69       LOGICAL, SAVE                 ::  is_first_run = .TRUE.
     70       LOGICAL                       ::  suppress_message
     71
    6772
    6873       IF (is_first_run)  THEN
     
    7378       END IF
    7479         
    75        WRITE(u, *)  TRIM(message) // "  [ " // TRIM(routine) // " ]"
     80
     81       suppress_message = .FALSE.
     82       IF (PRESENT(debug))  THEN
     83          IF (.NOT. debug)  suppress_message = .TRUE.
     84       END IF
     85
     86       IF (.NOT. suppress_message)  THEN
     87          PRINT *, "inifor: " // TRIM(message) // "  [ " // TRIM(routine) // " ]"
     88          WRITE(u, *)  TRIM(message) // "  [ " // TRIM(routine) // " ]"
     89       END IF
    7690
    7791       CLOSE(u)
  • palm/trunk/UTIL/inifor/src/defs.f90

    r3183 r3395  
    2626! -----------------
    2727! $Id$
     28! New parameters for computation of geostrophic winds
     29! Bumped INIFOR version number
     30!
     31!
     32! 3183 2018-07-27 14:25:55Z suehring
    2833! Updated defaults for soil extrapolation steps and nudging time-scale
    2934! Improved handling of the start date string
     
    7378 REAL(dp), PARAMETER ::  G    = 9.80665_dp                    !< acceleration of Earth's gravity, used in computation of COSMO-DE's basic state [m/s/s]
    7479 REAL(dp), PARAMETER ::  RHO_L = 1e3_dp                       !< density of liquid water, used to convert W_SO from [kg/m^2] to [m^3/m^3], in [kg/m^3]
     80 REAL(dp), PARAMETER ::  HECTO = 100_dp                       !< unit conversion factor from hPa to Pa
    7581
    7682 ! PALM-4U parameters
     83 REAL(dp), PARAMETER ::  OMEGA   = 7.29e-5_dp                 !< angular velocity of Earth's rotation [s^-1]
    7784 REAL(dp), PARAMETER ::  P_REF   = 1e5_dp                     !< Reference pressure for potential temperature [Pa]
    7885 REAL(dp), PARAMETER ::  RD_PALM = 287.0_dp                   !< specific gas constant of dry air, used in computation of PALM-4U's potential temperature [J/kg/K]
     
    8390 INTEGER, PARAMETER          ::  FORCING_STEP = 1             !< Number of hours between forcing time steps [h]
    8491 REAL(dp), PARAMETER         ::  NUDGING_TAU = 21600.0_dp     !< Nudging relaxation time scale [s]
    85  CHARACTER(LEN=*), PARAMETER ::  VERSION = '1.3.0'            !< INIFOR version number
     92 CHARACTER(LEN=*), PARAMETER ::  VERSION = '1.4.0'            !< INIFOR version number
    8693 CHARACTER(LEN=*), PARAMETER ::  COPYRIGHT = 'Copyright 2017-2018 Leibniz Universitaet Hannover' // &
    8794     NEW_LINE(' ') // ' Copyright 2017-2018 Deutscher Wetterdienst Offenbach' !< Copyright notice
  • palm/trunk/UTIL/inifor/src/grid.f90

    r3183 r3395  
    2626! -----------------
    2727! $Id$
     28! Added computation of geostrophic winds form COSMO pressure fields
     29! Introduced averaging grids and internal 'output' variables for computation of
     30!     geostrophic and large-scale forcings
     31! Merged qv, uv and vg output variables and 'temperature' IO group into new
     32!     'thermodynamics' IP group
     33! Cleaned up terminal output, show some messages only with --debug option
     34! Improved variable naming and readability
     35!
     36!
     37! 3183 2018-07-27 14:25:55Z suehring
    2838! Introduced new PALM grid stretching
    2939! Updated variable names and metadata for PIDS v1.9 compatibility
     
    6171        ONLY:  DATE, EARTH_RADIUS, TO_RADIANS, TO_DEGREES, PI, dp, hp, sp,     &
    6272               SNAME, LNAME, PATH, FORCING_STEP, WATER_ID, FILL_ITERATIONS,    &
    63                BETA, P_SL, T_SL, BETA, RD, G, P_REF, RD_PALM, CP_PALM, RHO_L
     73               BETA, P_SL, T_SL, BETA, RD, RV, G, P_REF, RD_PALM, CP_PALM,     &
     74               RHO_L, OMEGA, HECTO
    6475    USE io,                                                                    &
    6576        ONLY:  get_netcdf_variable, get_netcdf_attribute,                      &
     
    6879        ONLY:  NF90_MAX_NAME, NF90_MAX_VAR_DIMS
    6980    USE transform,                                                             &
    70         ONLY:  rotate_to_cosmo, find_horizontal_neighbours,                    &
     81        ONLY:  average_2d, rotate_to_cosmo, find_horizontal_neighbours,&
    7182               compute_horizontal_interp_weights,                              &
    72                find_vertical_neighbours_and_weights, interpolate_2d,           &
    73                gamma_from_hemisphere, phic_to_phin, lamc_to_lamn, average_2d,  &
    74                project, centre_velocities, phi2phirot, rla2rlarot, uv2uvrot
     83               find_vertical_neighbours_and_weights_interp,                    &
     84               find_vertical_neighbours_and_weights_average, interpolate_2d,   &
     85               gamma_from_hemisphere, phic_to_phin, lamc_to_lamn,              &
     86               project, centre_velocities, phi2phirot, rla2rlarot, uv2uvrot,   &
     87               phirot2phi, rlarot2rla
    7588    USE types
    7689    USE util
     
    8093    SAVE
    8194   
     95    REAL(dp) ::  averaging_angle   = 0.0_dp       !< latitudal and longitudal width of averaging regions [rad]
     96    REAL(dp) ::  averaging_width_ns = 0.0_dp       !< longitudal width of averaging regions [m]
     97    REAL(dp) ::  averaging_width_ew = 0.0_dp       !< latitudal width of averaging regions [m]
    8298    REAL(dp) ::  phi_equat         = 0.0_dp       !< latitude of rotated equator of COSMO-DE grid [rad]
    8399    REAL(dp) ::  phi_n             = 0.0_dp       !< latitude of rotated pole of COSMO-DE grid [rad]
     
    87103    REAL(dp) ::  phi_cn            = 0.0_dp       !< latitude of the rotated pole relative to the COSMO-DE grid [rad]
    88104    REAL(dp) ::  lambda_cn         = 0.0_dp       !< longitude of the rotated pole relative to the COSMO-DE grid [rad]
     105    REAL(dp) ::  lam_centre        = 0.0_dp       !< longitude of the PLAM domain centre in the source (COSMO rotated-pole) system [rad]
     106    REAL(dp) ::  phi_centre        = 0.0_dp       !< latitude of the PLAM domain centre in the source (COSMO rotated-pole) system [rad]
     107    REAL(dp) ::  lam_east          = 0.0_dp       !< longitude of the east central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad]
     108    REAL(dp) ::  lam_west          = 0.0_dp       !< longitude of the west central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad]
     109    REAL(dp) ::  phi_north         = 0.0_dp       !< latitude of the north central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad]
     110    REAL(dp) ::  phi_south         = 0.0_dp       !< latitude of the south central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad]
    89111    REAL(dp) ::  gam               = 0.0_dp       !< angle for working around phirot2phi/rlarot2rla bug
    90112    REAL(dp) ::  dx                = 0.0_dp       !< PALM-4U grid spacing in x direction [m]
     
    100122    REAL(dp) ::  dyi               = 0.0_dp       !< inverse PALM-4U grid spacing in y direction [m^-1]
    101123    REAL(dp) ::  dzi               = 0.0_dp       !< inverse PALM-4U grid spacing in z direction [m^-1]
     124    REAL(dp) ::  f3                = 0.0_dp       !< Coriolis parameter
    102125    REAL(dp) ::  lx                = 0.0_dp       !< PALM-4U domain size in x direction [m]
    103126    REAL(dp) ::  ly                = 0.0_dp       !< PALM-4U domain size in y direction [m]
    104     REAL(dp) ::  ug                = 0.0_dp       !< geostrophic wind in x direction [m/s]
    105     REAL(dp) ::  vg                = 0.0_dp       !< geostrophic wind in y direction [m/s]
    106127    REAL(dp) ::  p0                = 0.0_dp       !< PALM-4U surface pressure, at z0 [Pa]
    107128    REAL(dp) ::  x0                = 0.0_dp       !< x coordinate of PALM-4U Earth tangent [m]
     
    140161    INTEGER ::  dz_stretch_level_end_index(9)               !< vertical grid level index until which the vertical grid spacing is stretched
    141162    INTEGER ::  dz_stretch_level_start_index(9)             !< vertical grid level index above which the vertical grid spacing is stretched
    142     INTEGER ::  i     !< indexing variable
    143     INTEGER ::  average_imin, average_imax, average_jmin, average_jmax !< index ranges for profile averaging
    144     INTEGER ::  k     !< indexing variable
    145163    INTEGER ::  nt    !< number of output time steps
    146164    INTEGER ::  nx    !< number of PALM-4U grid points in x direction
     
    164182    LOGICAL ::  profile_grids_required
    165183
    166     INTEGER ::  n_invar = 0 !< number of variables in the input variable table
    167     INTEGER ::  n_outvar = 0 !< number of variables in the output variable table
    168184    TYPE(nc_var), ALLOCATABLE, TARGET ::  input_var_table(:)  !< table of input variables
    169185    TYPE(nc_var), ALLOCATABLE, TARGET ::  output_var_table(:) !< table of input variables
     
    219235    TYPE(grid_definition), TARGET ::  w_south_intermediate            !<
    220236    TYPE(grid_definition), TARGET ::  w_top_intermediate              !<
    221     TYPE(grid_definition), TARGET ::  scalar_profile_grid             !<
    222     TYPE(grid_definition), TARGET ::  scalar_profile_intermediate     !<
    223     TYPE(grid_definition), TARGET ::  w_profile_grid                  !<
    224     TYPE(grid_definition), TARGET ::  w_profile_intermediate          !<
     237
     238    TYPE(grid_definition), TARGET ::  north_averaged_scalar_profile
     239    TYPE(grid_definition), TARGET ::  south_averaged_scalar_profile
     240    TYPE(grid_definition), TARGET ::  west_averaged_scalar_profile
     241    TYPE(grid_definition), TARGET ::  east_averaged_scalar_profile
     242    TYPE(grid_definition), TARGET ::  averaged_scalar_profile
     243    TYPE(grid_definition), TARGET ::  averaged_w_profile
    225244
    226245    TYPE(io_group), ALLOCATABLE, TARGET ::  io_group_list(:)  !< List of I/O groups, which group together output variables that share the same input variable
     
    237256    CHARACTER(LEN=PATH), ALLOCATABLE, DIMENSION(:) ::  radiation_files
    238257    CHARACTER(LEN=SNAME) ::  input_prefix         !< Prefix of input files, e.g. 'laf' for COSMO-DE analyses
     258    CHARACTER(LEN=SNAME) ::  flow_prefix          !< Prefix of flow input files, e.g. 'laf' for COSMO-DE analyses
     259    CHARACTER(LEN=SNAME) ::  soil_prefix          !< Prefix of soil input files, e.g. 'laf' for COSMO-DE analyses
     260    CHARACTER(LEN=SNAME) ::  radiation_prefix     !< Prefix of radiation input files, e.g. 'laf' for COSMO-DE analyses
     261    CHARACTER(LEN=SNAME) ::  soilmoisture_prefix  !< Prefix of input files for soil moisture spin-up, e.g. 'laf' for COSMO-DE analyses
    239262    CHARACTER(LEN=SNAME) ::  flow_suffix          !< Suffix of flow input files, e.g. 'flow'
    240263    CHARACTER(LEN=SNAME) ::  soil_suffix          !< Suffix of soil input files, e.g. 'soil'
     
    249272   
    250273    SUBROUTINE setup_parameters()
     274       INTEGER ::  k
    251275
    252276!
     
    282306
    283307       ! Default atmospheric parameters
    284        ug = 0.0_dp
    285        vg = 0.0_dp
     308       cfg % ug = 0.0_dp
     309       cfg % vg = 0.0_dp
    286310       cfg % p0 = P_SL
    287311
     
    296320       input_prefix = 'laf'  ! 'laf' for COSMO-DE analyses
    297321       cfg % flow_prefix = input_prefix
     322       cfg % input_prefix = input_prefix
    298323       cfg % soil_prefix = input_prefix
    299324       cfg % radiation_prefix = input_prefix
     
    304329       radiation_suffix = '-rad'
    305330       soilmoisture_suffix = '-soilmoisture'
     331
     332       cfg % debug = .FALSE.
     333       cfg % averaging_angle = 2.0_dp
    306334!
    307335!------------------------------------------------------------------------------
     
    318346       cfg % ic_mode = 'volume'
    319347       cfg % bc_mode = 'real'
    320 
    321        ! Update default file names and domain centre
     348       cfg % averaging_mode = 'level'
     349
     350       ! Overwrite defaults with user configuration
    322351       CALL parse_command_line_arguments( cfg )
     352
     353       flow_prefix = TRIM(cfg % input_prefix)
     354       radiation_prefix = TRIM(cfg % input_prefix)
     355       soil_prefix = TRIM(cfg % input_prefix)
     356       soilmoisture_prefix = TRIM(cfg % input_prefix)
     357       IF (cfg % flow_prefix_is_set)  flow_prefix = TRIM(cfg % flow_prefix)
     358       IF (cfg % radiation_prefix_is_set)  radiation_prefix = TRIM(cfg % radiation_prefix)
     359       IF (cfg % soil_prefix_is_set)  soil_prefix = TRIM(cfg % soil_prefix)
     360       IF (cfg % soilmoisture_prefix_is_set)  soilmoisture_prefix = TRIM(cfg % soilmoisture_prefix)
    323361
    324362       output_file % name = cfg % output_file
     
    336374       END IF
    337375
     376
     377       ! Set default file paths, if not specified by user.
    338378       CALL normalize_path(cfg % input_path)
    339379       IF (TRIM(cfg % hhl_file) == '')  cfg % hhl_file = TRIM(cfg % input_path) // 'hhl.nc'
     
    361401
    362402       ! Generate input file lists
    363        CALL get_input_file_list(cfg % start_date, start_hour_flow, end_hour, step_hour,  &
    364           cfg % input_path, cfg % flow_prefix, flow_suffix, flow_files)
    365        CALL get_input_file_list(cfg % start_date, start_hour_soil, end_hour, step_hour,  &
    366           cfg % input_path, cfg % soil_prefix, soil_suffix, soil_files)
    367        CALL get_input_file_list(cfg % start_date, start_hour_radiation, end_hour, step_hour, &
    368           cfg % input_path, cfg % radiation_prefix, radiation_suffix, radiation_files)
    369        CALL get_input_file_list(cfg % start_date, start_hour_soilmoisture, end_hour_soilmoisture, step_hour, &
    370           cfg % input_path, cfg % soilmoisture_prefix, soilmoisture_suffix, soil_moisture_files)
     403       CALL get_input_file_list(                                               &
     404          cfg % start_date, start_hour_flow, end_hour, step_hour,              &
     405          cfg % input_path, flow_prefix, flow_suffix, flow_files)
     406       CALL get_input_file_list(                                               &
     407          cfg % start_date, start_hour_soil, end_hour, step_hour,              &
     408          cfg % input_path, soil_prefix, soil_suffix, soil_files)
     409       CALL get_input_file_list(                                               &
     410          cfg % start_date, start_hour_radiation, end_hour, step_hour,         &
     411          cfg % input_path, radiation_prefix, radiation_suffix, radiation_files)
     412       CALL get_input_file_list(                                               &
     413          cfg % start_date, start_hour_soilmoisture, end_hour_soilmoisture, step_hour, &
     414          cfg % input_path, soilmoisture_prefix, soilmoisture_suffix, soil_moisture_files)
    371415
    372416!
     
    499543
    500544 CALL run_control('time', 'comp')
     545
     546!------------------------------------------------------------------------------
     547! Section 4.3: INIFOR averaging domains
     548!------------------------------------------------------------------------------
     549
     550    ! Compute coordiantes of the PALM centre in the source (COSMO) system
     551    phi_centre = phirot2phi(                                                   &
     552       phirot = project(0.5_dp*ly, y0, EARTH_RADIUS) * TO_DEGREES,             &
     553       rlarot = project(0.5_dp*lx, x0, EARTH_RADIUS) * TO_DEGREES,             &
     554       polphi = phi_cn * TO_DEGREES, pollam = lambda_cn * TO_DEGREES,          &
     555       polgam = gam * TO_DEGREES                                               &
     556    ) * TO_RADIANS
     557
     558    lam_centre = rlarot2rla(                                                   &
     559       phirot = project(0.5_dp*ly, y0, EARTH_RADIUS) * TO_DEGREES,             &
     560       rlarot = project(0.5_dp*lx, x0, EARTH_RADIUS) * TO_DEGREES,             &
     561       polphi = phi_cn * TO_DEGREES, pollam = lambda_cn * TO_DEGREES,          &
     562       polgam = gam * TO_DEGREES                                               &
     563    ) * TO_RADIANS
     564
     565    message = "PALM-4U centre:" // NEW_LINE('') // &
     566       "           lon (lambda) = " // &
     567       TRIM(real_to_str_f(lam_centre * TO_DEGREES)) // " deg" // NEW_LINE(' ') // &
     568       "           lat (phi   ) = " // &
     569       TRIM(real_to_str_f(phi_centre * TO_DEGREES)) // " deg (COSMO-DE rotated-pole)"
     570    CALL report( 'setup_parameters', message )
     571
     572
     573    ! Compute boundaries of the central averaging box
     574    averaging_angle = cfg % averaging_angle * TO_RADIANS
     575    lam_east = lam_centre + 0.5_dp * averaging_angle
     576    lam_west = lam_centre - 0.5_dp * averaging_angle
     577    phi_north = phi_centre + 0.5_dp * averaging_angle
     578    phi_south = phi_centre - 0.5_dp * averaging_angle
     579    averaging_width_ew = averaging_angle * COS(phi_centre) * EARTH_RADIUS
     580    averaging_width_ns = averaging_angle * EARTH_RADIUS
     581
     582    ! Coriolis parameter
     583    f3 = 2.0_dp * OMEGA * SIN(                                                 &
     584       TO_RADIANS*phirot2phi( phi_centre * TO_DEGREES, lam_centre * TO_DEGREES,&
     585                              phi_n * TO_DEGREES, lambda_n * TO_DEGREES,       &
     586                              gam * TO_DEGREES )                               &
     587       )
    501588
    502589    END SUBROUTINE setup_parameters
     
    603690               ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy,                    &
    604691               x0=x0, y0=y0, z0 = z0,                                          &
    605                nx = nx, ny = ny, nz = nlev - 2)
     692               nx = nx, ny = ny, nz = nlev - 1)
    606693
    607694      IF (boundary_variables_required)  THEN
     
    824911                  ymin =  0.5_dp * dy, ymax = ly - 0.5_dp * dy,                   &
    825912                  x0=x0, y0=y0, z0 = z0,                                          &
    826                   nx = 0, ny = ny, nz = nlev - 2)
     913                  nx = 0, ny = ny, nz = nlev - 1)
    827914
    828915          CALL init_grid_definition('boundary intermediate', grid=w_west_intermediate,         &
     
    830917                  ymin =  0.5_dp * dy, ymax = ly - 0.5_dp * dy,                   &
    831918                  x0=x0, y0=y0, z0 = z0,                                          &
    832                   nx = 0, ny = ny, nz = nlev - 2)
     919                  nx = 0, ny = ny, nz = nlev - 1)
    833920
    834921          CALL init_grid_definition('boundary intermediate', grid=w_north_intermediate,        &
     
    836923                  ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy,               &
    837924                  x0=x0, y0=y0, z0 = z0,                                          &
    838                   nx = nx, ny = 0, nz = nlev - 2)
     925                  nx = nx, ny = 0, nz = nlev - 1)
    839926
    840927          CALL init_grid_definition('boundary intermediate', grid=w_south_intermediate,        &
     
    842929                  ymin = -0.5_dp * dy, ymax = -0.5_dp * dy,                       &
    843930                  x0=x0, y0=y0, z0 = z0,                                          &
    844                   nx = nx, ny = 0, nz = nlev - 2)
     931                  nx = nx, ny = 0, nz = nlev - 1)
    845932
    846933          CALL init_grid_definition('boundary intermediate', grid=w_top_intermediate,          &
     
    848935                  ymin =  0.5_dp * dy, ymax = ly - 0.5_dp * dy,                   &
    849936                  x0=x0, y0=y0, z0 = z0,                                          &
    850                   nx = nx, ny = ny, nz = nlev - 2)
     937                  nx = nx, ny = ny, nz = nlev - 1)
    851938       END IF
    852939
     
    856943!------------------------------------------------------------------------------
    857944
    858        profile_grids_required = ( TRIM(cfg % ic_mode) == 'profile' .OR.              &
     945       CALL init_averaging_grid(averaged_scalar_profile, cosmo_grid,           &
     946               x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0,               &
     947               lonmin = lam_west, lonmax = lam_east,                           &
     948               latmin = phi_south, latmax = phi_north,                         &
     949               kind='scalar')
     950
     951       CALL init_averaging_grid(averaged_w_profile, cosmo_grid,                &
     952               x = 0.5_dp * lx, y = 0.5_dp * ly, z = zw, z0 = z0,              &
     953               lonmin = lam_west, lonmax = lam_east,                           &
     954               latmin = phi_south, latmax = phi_north,                         &
     955               kind='w')
     956
     957       CALL init_averaging_grid(south_averaged_scalar_profile, cosmo_grid,     &
     958               x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0,               &
     959               lonmin = lam_west, lonmax = lam_east,                           &
     960               latmin = phi_centre - averaging_angle,                          &
     961               latmax = phi_centre,                                            &
     962               kind='scalar')
     963
     964       CALL init_averaging_grid(north_averaged_scalar_profile, cosmo_grid,     &
     965               x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0,               &
     966               lonmin = lam_west, lonmax = lam_east,                           &
     967               latmin = phi_centre,                                            &
     968               latmax = phi_centre + averaging_angle,                          &
     969               kind='scalar')
     970
     971       CALL init_averaging_grid(west_averaged_scalar_profile, cosmo_grid,      &
     972               x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0,               &
     973               lonmin = lam_centre - averaging_angle,                          &
     974               lonmax = lam_centre,                                            &
     975               latmin = phi_south, latmax = phi_north,                         &
     976               kind='scalar')
     977
     978       CALL init_averaging_grid(east_averaged_scalar_profile, cosmo_grid,      &
     979               x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0,               &
     980               lonmin = lam_centre,                                            &
     981               lonmax = lam_centre + averaging_angle,                          &
     982               latmin = phi_south, latmax = phi_north,                         &
     983               kind='scalar')
     984
     985       profile_grids_required = ( TRIM(cfg % ic_mode) == 'profile' .OR.        &
    859986                                  TRIM(cfg % bc_mode) == 'ideal' )
    860987
    861        IF (profile_grids_required)  THEN
    862           CALL init_grid_definition('boundary', grid=scalar_profile_grid,      &
    863                   xmin = 0.5_dp * lx, xmax = 0.5_dp * lx,                      &
    864                   ymin = 0.5_dp * ly, ymax = 0.5_dp * ly,                      &
    865                   x0=x0, y0=y0, z0 = z0,                                       &
    866                   nx = 0, ny = 0, nz = nz, z=z)
    867 
    868           CALL init_grid_definition('boundary', grid=w_profile_grid,           &
    869                   xmin = 0.5_dp * lx, xmax = 0.5_dp * lx,                      &
    870                   ymin = 0.5_dp * ly, ymax = 0.5_dp * ly,                      &
    871                   x0=x0, y0=y0, z0 = z0,                                       &
    872                   nx = 0, ny = 0, nz = nz - 1, z=zw)
    873 
    874           CALL init_grid_definition('boundary', grid=scalar_profile_intermediate,&
    875                   xmin = 0.5_dp * lx, xmax = 0.5_dp * lx,                      &
    876                   ymin = 0.5_dp * ly, ymax = 0.5_dp * ly,                      &
    877                   x0=x0, y0=y0, z0 = z0,                                       &
    878                   nx = 0, ny = 0, nz = nlev - 2, z=z)
    879 
    880           CALL init_grid_definition('boundary', grid=w_profile_intermediate,   &
    881                   xmin = 0.5_dp * lx, xmax = 0.5_dp * lx,                      &
    882                   ymin = 0.5_dp * ly, ymax = 0.5_dp * ly,                      &
    883                   x0=x0, y0=y0, z0 = z0,                                       &
    884                   nx = 0, ny = 0, nz = nlev - 2, z=zw)
    885        END IF
     988       !IF (profile_grids_required)  THEN
     989       !   ! Here, I use the 'boundary' kind, so the vertical coordinate uses the
     990       !   ! PALM z coordinate.
     991       !   CALL init_grid_definition('boundary', grid=scalar_profile_grid,      &
     992       !           xmin = 0.5_dp * lx, xmax = 0.5_dp * lx,                      &
     993       !           ymin = 0.5_dp * ly, ymax = 0.5_dp * ly,                      &
     994       !           x0=x0, y0=y0, z0 = z0,                                       &
     995       !           nx = 0, ny = 0, nz = nz, z=z)
     996
     997       !   CALL init_grid_definition('boundary', grid=w_profile_grid,           &
     998       !           xmin = 0.5_dp * lx, xmax = 0.5_dp * lx,                      &
     999       !           ymin = 0.5_dp * ly, ymax = 0.5_dp * ly,                      &
     1000       !           x0=x0, y0=y0, z0 = z0,                                       &
     1001       !           nx = 0, ny = 0, nz = nz - 1, z=zw)
     1002
     1003       !   CALL init_grid_definition('boundary', grid=scalar_profile_intermediate,&
     1004       !           xmin = 0.5_dp * lx, xmax = 0.5_dp * lx,                      &
     1005       !           ymin = 0.5_dp * ly, ymax = 0.5_dp * ly,                      &
     1006       !           x0=x0, y0=y0, z0 = z0,                                       &
     1007       !           nx = 0, ny = 0, nz = nlev - 2, z=z)
     1008
     1009       !   CALL init_grid_definition('boundary', grid=w_profile_intermediate,   &
     1010       !           xmin = 0.5_dp * lx, xmax = 0.5_dp * lx,                      &
     1011       !           ymin = 0.5_dp * ly, ymax = 0.5_dp * ly,                      &
     1012       !           x0=x0, y0=y0, z0 = z0,                                       &
     1013       !           nx = 0, ny = 0, nz = nlev - 2, z=zw)
     1014       !END IF
    8861015
    8871016!                                                                             
     
    9301059
    9311060       IF (TRIM(cfg % ic_mode) == 'profile')  THEN
    932            CALL setup_averaging(cosmo_grid, palm_intermediate,                 &
    933                 average_imin, average_imax, average_jmin, average_jmax)
     1061           !TODO: remove this conditional if not needed.
    9341062       END IF
    9351063       
     
    9381066
    9391067
    940     SUBROUTINE setup_interpolation(cosmo_grid, palm_grid, palm_intermediate, kind, ic_mode)
     1068    SUBROUTINE setup_interpolation(cosmo_grid, grid, intermediate_grid, kind, ic_mode)
    9411069
    9421070       TYPE(grid_definition), INTENT(IN), TARGET    ::  cosmo_grid
    943        TYPE(grid_definition), INTENT(INOUT), TARGET ::  palm_grid, palm_intermediate
     1071       TYPE(grid_definition), INTENT(INOUT), TARGET ::  grid, intermediate_grid
    9441072       CHARACTER, INTENT(IN)                        ::  kind
    9451073       CHARACTER(LEN=*), INTENT(IN), OPTIONAL       ::  ic_mode
    9461074
    947        TYPE(grid_definition), POINTER      ::  grid
    948        REAL(dp), DIMENSION(:), POINTER     ::  lat, lon
    949        REAL(dp), DIMENSION(:,:,:), POINTER ::  h
    950 
    951        LOGICAL :: setup_vertical
    952 
    953        setup_vertical = .TRUE.
    954        IF (PRESENT(ic_mode))  THEN
    955           IF (TRIM(ic_mode) == 'profile')  setup_vertical = .FALSE.
    956        END IF
     1075       REAL(dp), DIMENSION(:), POINTER     ::  cosmo_lat, cosmo_lon
     1076       REAL(dp), DIMENSION(:,:,:), POINTER ::  cosmo_h
     1077
     1078       LOGICAL :: setup_volumetric
    9571079
    9581080!------------------------------------------------------------------------------
     
    9641086       CASE('s') ! scalars
    9651087
    966           lat => cosmo_grid % lat
    967           lon => cosmo_grid % lon
    968           h   => cosmo_grid % hfl
     1088          cosmo_lat => cosmo_grid % lat
     1089          cosmo_lon => cosmo_grid % lon
     1090          cosmo_h   => cosmo_grid % hfl
    9691091
    9701092       CASE('w') ! vertical velocity
    9711093
    972           lat => cosmo_grid % lat
    973           lon => cosmo_grid % lon
    974           h   => cosmo_grid % hhl
     1094          cosmo_lat => cosmo_grid % lat
     1095          cosmo_lon => cosmo_grid % lon
     1096          cosmo_h   => cosmo_grid % hhl
    9751097
    9761098       CASE('u') ! x velocity
    9771099
    978           lat => cosmo_grid % lat
    979           lon => cosmo_grid % lonu
    980           h   => cosmo_grid % hfl
     1100          cosmo_lat => cosmo_grid % lat
     1101          cosmo_lon => cosmo_grid % lonu
     1102          cosmo_h   => cosmo_grid % hfl
    9811103
    9821104       CASE('v') ! y velocity
    9831105
    984           lat => cosmo_grid % latv
    985           lon => cosmo_grid % lon
    986           h   => cosmo_grid % hfl
     1106          cosmo_lat => cosmo_grid % latv
     1107          cosmo_lon => cosmo_grid % lon
     1108          cosmo_h   => cosmo_grid % hfl
    9871109
    9881110       CASE DEFAULT
     
    9931115       END SELECT
    9941116
    995        grid => palm_intermediate
    996 
    997        CALL find_horizontal_neighbours(lat, lon,                               &
    998           grid % clat, grid % clon, grid % ii, grid % jj)
    999 
    1000        CALL compute_horizontal_interp_weights(lat, lon,                        &
    1001           grid % clat, grid % clon, grid % ii, grid % jj, grid % w_horiz)
     1117       CALL find_horizontal_neighbours(cosmo_lat, cosmo_lon,                   &
     1118          intermediate_grid % clat, intermediate_grid % clon,                  &
     1119          intermediate_grid % ii, intermediate_grid % jj)
     1120
     1121       CALL compute_horizontal_interp_weights(cosmo_lat, cosmo_lon,            &
     1122          intermediate_grid % clat, intermediate_grid % clon,                  &
     1123          intermediate_grid % ii, intermediate_grid % jj,                      &
     1124          intermediate_grid % w_horiz)
    10021125
    10031126!------------------------------------------------------------------------------
     
    10051128!------------------------------------------------------------------------------
    10061129
    1007        IF (setup_vertical)  THEN
    1008           ALLOCATE( grid % h(0:grid % nx, 0:grid % ny, 0:grid % nz) )
    1009           grid % h(:,:,:) = - EARTH_RADIUS
     1130       ! If profile initialization is chosen, we--somewhat counterintuitively--
     1131       ! don't need to compute vertical interpolation weights. At least, we
     1132       ! don't need them on the intermediate grid, which fills the entire PALM
     1133       ! domain volume. Instead we need vertical weights for the intermediate
     1134       ! profile grids, which get computed in setup_averaging().
     1135       setup_volumetric = .TRUE.
     1136       IF (PRESENT(ic_mode))  THEN
     1137          IF (TRIM(ic_mode) == 'profile')  setup_volumetric = .FALSE.
     1138       END IF
     1139
     1140       IF (setup_volumetric)  THEN
     1141          ALLOCATE( intermediate_grid % h(0:intermediate_grid % nx,            &
     1142                                          0:intermediate_grid % ny,            &
     1143                                          0:intermediate_grid % nz) )
     1144          intermediate_grid % h(:,:,:) = - EARTH_RADIUS
    10101145
    10111146          ! For w points, use hhl, for scalars use hfl
    10121147          ! compute the full heights for the intermediate grids
    1013           CALL interpolate_2d(cosmo_grid % hfl, grid % h, grid)
    1014           CALL find_vertical_neighbours_and_weights(palm_grid, grid)
     1148          CALL interpolate_2d(cosmo_h, intermediate_grid % h, intermediate_grid)
     1149          CALL find_vertical_neighbours_and_weights_interp(grid, intermediate_grid)
    10151150       END IF
    10161151       
    10171152    END SUBROUTINE setup_interpolation
    1018 
    1019 
    1020     SUBROUTINE setup_averaging(cosmo_grid, palm_intermediate, imin, imax, jmin, jmax)
    1021 
    1022        TYPE(grid_definition), INTENT(IN) ::  cosmo_grid, palm_intermediate
    1023        INTEGER, INTENT(INOUT)            ::  imin, imax, jmin, jmax
    1024 
    1025        TYPE(grid_definition), POINTER    ::  grid
    1026        REAL(dp)                          ::  lonmin_pos,lonmax_pos, latmin_pos, latmax_pos
    1027        REAL(dp)                          ::  cosmo_dxi, cosmo_dyi
    1028 
    1029        cosmo_dxi = 1.0_dp / (cosmo_grid % lon(1) - cosmo_grid % lon(0))
    1030        cosmo_dyi = 1.0_dp / (cosmo_grid % lat(1) - cosmo_grid % lat(0))
    1031 
    1032        ! find horizontal index ranges for profile averaging
    1033        lonmin_pos = (MINVAL(palm_intermediate % clon(:,:)) - cosmo_grid % lon(0)) * cosmo_dxi
    1034        lonmax_pos = (MAXVAL(palm_intermediate % clon(:,:)) - cosmo_grid % lon(0)) * cosmo_dxi
    1035        latmin_pos = (MINVAL(palm_intermediate % clat(:,:)) - cosmo_grid % lat(0)) * cosmo_dyi
    1036        latmax_pos = (MAXVAL(palm_intermediate % clat(:,:)) - cosmo_grid % lat(0)) * cosmo_dyi
    1037 
    1038        imin = FLOOR(lonmin_pos)
    1039        imax = CEILING(lonmax_pos)
    1040        jmin = FLOOR(latmin_pos)
    1041        jmax = CEILING(latmax_pos)
    1042        
    1043        ! average heights for intermediate scalar and w profile grids
    1044        grid => scalar_profile_intermediate
    1045        ALLOCATE( grid % h(0:grid % nx, 0:grid % ny, 0:grid % nz) )
    1046        grid % h(:,:,:) = - EARTH_RADIUS
    1047        CALL average_2d(cosmo_grid % hfl, grid % h(0,0,:), imin, imax, jmin, jmax)
    1048        CALL find_vertical_neighbours_and_weights(scalar_profile_grid, grid)
    1049 
    1050        grid => w_profile_intermediate
    1051        ALLOCATE( grid % h(0:grid % nx, 0:grid % ny, 0:grid % nz) )
    1052        grid % h(:,:,:) = - EARTH_RADIUS
    1053        CALL average_2d(cosmo_grid % hhl, grid % h(0,0,:), imin, imax, jmin, jmax)
    1054        CALL find_vertical_neighbours_and_weights(w_profile_grid, grid)
    1055 
    1056     END SUBROUTINE setup_averaging
    1057 
    10581153
    10591154!------------------------------------------------------------------------------!
     
    12911386
    12921387    END SUBROUTINE init_grid_definition
     1388
     1389
     1390!------------------------------------------------------------------------------!
     1391! Description:
     1392! ------------
     1393!> Initializes averagin_grid-type variables.
     1394!>
     1395!> Input parameters:
     1396!> -----------------
     1397!>
     1398!> cosmo_grid : grid_definition-type varieble initialized as 'cosmo-de' grid
     1399!>    providing longitudes and latitudes of the parent grid
     1400!>
     1401!> x, y : location of the profile to be averaged in the PALM system [m]
     1402!>
     1403!> z : vertical grid coordinates of the profile in the PALM system [m]
     1404!>
     1405!> <lat/lon><min/max>: boundaries of the averaging region in the parent system,
     1406!>    i.e. the COSMO-DE rotated-pole system. [rad]
     1407!>
     1408!> kind : kind of quantity to be averaged using this averaging grid.
     1409!>    Destinguishes COSMO-DE scalar and w-velocity levels. Note that finding the
     1410!>    parent/COSMO columns for the region in get_cosmo_averaging_region() is
     1411!>    independent of 'kind' b/ecause INIFOR uses column-centred u and v velocity
     1412!>    components, which are computed in the preprocessing step.
     1413!>
     1414!> Output parameters:
     1415!> ------------------
     1416!> avg_grid : averaging_grid variable to be initialized.
     1417!------------------------------------------------------------------------------!
     1418    SUBROUTINE init_averaging_grid(avg_grid, cosmo_grid, x, y, z, z0,          &
     1419       lonmin, lonmax, latmin, latmax, kind)
     1420
     1421       TYPE(grid_definition), INTENT(INOUT) ::  avg_grid
     1422       TYPE(grid_definition), INTENT(IN)    ::  cosmo_grid
     1423       REAL(dp), INTENT(IN)                 ::  x, y, z0
     1424       REAL(dp), INTENT(IN), TARGET         ::  z(:)
     1425       REAL(dp), INTENT(IN)                 ::  lonmin, lonmax, latmin, latmax
     1426
     1427       CHARACTER(LEN=*), INTENT(IN)         ::  kind
     1428
     1429       LOGICAL                              :: level_based_averaging
     1430
     1431       INTEGER :: i, j
     1432
     1433       ALLOCATE( avg_grid % x(1) )
     1434       ALLOCATE( avg_grid % y(1) )
     1435       avg_grid % x(1) = x
     1436       avg_grid % y(1) = y
     1437       avg_grid % z => z
     1438       avg_grid % z0 = z0
     1439
     1440       avg_grid % nz = SIZE(z, 1)
     1441
     1442       ALLOCATE( avg_grid % lon(2) )
     1443       ALLOCATE( avg_grid % lat(2) )
     1444       avg_grid % lon(1:2) = (/lonmin, lonmax/)
     1445       avg_grid % lat(1:2) = (/latmin, latmax/)
     1446
     1447       avg_grid % kind = TRIM(kind)
     1448
     1449       ! Find and store COSMO columns that fall into the coordinate range
     1450       ! given by avg_grid % clon, % clat
     1451       CALL get_cosmo_averaging_region(avg_grid, cosmo_grid)
     1452
     1453       ALLOCATE (avg_grid % kkk(avg_grid % n_columns, avg_grid % nz, 2) )
     1454       ALLOCATE (avg_grid % w(avg_grid % n_columns, avg_grid % nz, 2) )
     1455       ! Compute average COSMO levels in the averaging region
     1456       SELECT CASE(avg_grid % kind)
     1457
     1458       CASE('scalar', 'u', 'v')
     1459          avg_grid % cosmo_h => cosmo_grid % hfl
     1460
     1461       CASE('w')
     1462          avg_grid % cosmo_h => cosmo_grid % hhl
     1463
     1464       CASE DEFAULT
     1465          message = "Averaging grid kind '" // TRIM(avg_grid % kind) // &
     1466                    "' is not supported. Use 'scalar', 'u', or 'v'."
     1467          CALL abort('get_cosmo_averaging_region', message)
     1468
     1469       END SELECT
     1470
     1471       ! For level-besed averaging, compute average heights
     1472       level_based_averaging = .TRUE.
     1473       IF (level_based_averaging)  THEN
     1474          ALLOCATE(avg_grid % h(1,1,SIZE(avg_grid % cosmo_h, 3)) )
     1475 
     1476          CALL average_2d(avg_grid % cosmo_h, avg_grid % h(1,1,:),                            &
     1477                          avg_grid % iii, avg_grid % jjj)
     1478       END IF
     1479
     1480       ! Copy averaged grid to all COSMO columnts, leads to computing the same
     1481       ! vertical interpolation weights for all columns and to COSMO grid level
     1482       ! based averaging onto the averaged COSMO heights.
     1483       IF ( TRIM(cfg % averaging_mode) == 'level' )  THEN
     1484          FORALL(                                                              &
     1485             i=1:SIZE(avg_grid % cosmo_h, 1),                                  &
     1486             j=1:SIZE(avg_grid % cosmo_h, 2)                                   &
     1487          )   avg_grid % cosmo_h(i,j,:) = avg_grid % h(1,1,:)
     1488       END IF
     1489
     1490       ! Compute vertical weights and neighbours
     1491       CALL find_vertical_neighbours_and_weights_average(avg_grid)
     1492
     1493    END SUBROUTINE init_averaging_grid
     1494
     1495
     1496    SUBROUTINE get_cosmo_averaging_region(avg_grid, cosmo_grid)
     1497       TYPE(grid_definition), INTENT(INOUT)         ::  avg_grid
     1498       TYPE(grid_definition), TARGET, INTENT(IN)    ::  cosmo_grid
     1499
     1500       REAL(dp), DIMENSION(:), POINTER              ::  cosmo_lon, cosmo_lat
     1501       REAL(dp)                                     ::  dlon, dlat
     1502
     1503       INTEGER ::  i, j, imin, imax, jmin, jmax, l, nx, ny
     1504
     1505
     1506       SELECT CASE( TRIM(avg_grid % kind) )
     1507
     1508       CASE('scalar', 'w')
     1509          cosmo_lon => cosmo_grid % lon
     1510          cosmo_lat => cosmo_grid % lat
     1511
     1512       CASE('u')
     1513          cosmo_lon => cosmo_grid % lonu
     1514          cosmo_lat => cosmo_grid % lat
     1515
     1516       CASE('v')
     1517          cosmo_lon => cosmo_grid % lon
     1518          cosmo_lat => cosmo_grid % latv
     1519
     1520       CASE DEFAULT
     1521          message = "Averaging grid kind '" // TRIM(avg_grid % kind) // &
     1522                    "' is not supported. Use 'scalar', 'u', or 'v'."
     1523          CALL abort('get_cosmo_averaging_region', message)
     1524
     1525       END SELECT
     1526
     1527
     1528       ! FIXME: make dlon, dlat parameters of the grid_defintion type
     1529       dlon = cosmo_lon(1) - cosmo_lon(0)
     1530       dlat = cosmo_lat(1) - cosmo_lat(0)
     1531
     1532       imin = CEILING( (avg_grid % lon(1) - cosmo_lon(0)) / dlon )
     1533       imax = FLOOR(   (avg_grid % lon(2) - cosmo_lon(0)) / dlon )
     1534
     1535       jmin = CEILING( (avg_grid % lat(1) - cosmo_lat(0)) / dlat )
     1536       jmax = FLOOR(   (avg_grid % lat(2) - cosmo_lat(0)) / dlat )
     1537       
     1538       message = "Averaging over "//                                           &
     1539                 TRIM(str(imin)) // " < i < " // TRIM(str(imax)) //            &
     1540                 " and " //                                                    &
     1541                 TRIM(str(jmin)) // " < j < " // TRIM(str(jmax))
     1542       CALL report( 'get_cosmo_averaging_region', message )
     1543
     1544       nx = imax - imin + 1
     1545       ny = jmax - jmin + 1
     1546       avg_grid % n_columns = nx * ny
     1547
     1548       ALLOCATE( avg_grid % iii(avg_grid % n_columns),                                 &
     1549                 avg_grid % jjj(avg_grid % n_columns) )
     1550
     1551       l = 0
     1552       DO j = jmin, jmax
     1553       DO i = imin, imax
     1554          l = l + 1
     1555          avg_grid % iii(l) = i
     1556          avg_grid % jjj(l) = j
     1557       END DO
     1558       END DO
     1559
     1560    END SUBROUTINE get_cosmo_averaging_region
    12931561
    12941562
     
    17181986       )
    17191987
    1720        !potential temperature, surface pressure, including nudging and subsidence
     1988       !potential temperature, surface pressure, specific humidity including
     1989       !nudging and subsidence, and geostrophic winds ug, vg
    17211990       io_group_list(3) = init_io_group(                                       &
    17221991          in_files = flow_files,                                               &
    1723           out_vars = [output_var_table(3:8), output_var_table(42:42),          &
    1724                       output_var_table(49:51)],                                &
    1725           in_var_list = (/input_var_table(3), input_var_table(17)/),           &
    1726           kind = 'temperature'                                                 &
    1727        )
    1728 
    1729        !specific humidity including nudging and subsidence
    1730        io_group_list(4) = init_io_group(                                       &
    1731           in_files = flow_files,                                               &
    1732           out_vars = [output_var_table(9:14), output_var_table(52:54)],        &
    1733           in_var_list = input_var_table(4:4),                                  &
    1734           kind = 'scalar'                                                      &
    1735        )
    1736 
    1737        !u and v velocity and geostrophic winds ug, vg
     1992          out_vars = [output_var_table(56:64),                                 & ! internal averaged density and pressure profiles
     1993                      output_var_table(3:8), output_var_table(9:14),           &
     1994                      output_var_table(42:42), output_var_table(43:44),        &
     1995                      output_var_table(49:51), output_var_table(52:54)],       &
     1996          in_var_list = (/input_var_table(3), input_var_table(17),             & ! T, P, QV
     1997                          input_var_table(4) /),                               &
     1998          kind = 'thermodynamics',                                             &
     1999          n_output_quantities = 4                                              & ! P, Theta, Rho, qv
     2000       )
     2001
     2002       !Moved to therodynamic io_group
     2003       !io_group_list(4) = init_io_group(                                       &
     2004       !   in_files = flow_files,                                               &
     2005       !   out_vars = [output_var_table(9:14), output_var_table(52:54)],        &
     2006       !   in_var_list = input_var_table(4:4),                                  &
     2007       !   kind = 'scalar'                                                      &
     2008       !)
     2009
     2010       !u and v velocity, ug anv vg moved to thermodynamic io group.
    17382011       io_group_list(5) = init_io_group(                                       &
    17392012          in_files = flow_files,                                               &
    1740           out_vars = [output_var_table(15:26), output_var_table(43:46)],       &
     2013          out_vars = [output_var_table(15:26), output_var_table(45:46)],       &
    17412014          !out_vars = output_var_table(15:20),                                  &
    17422015          in_var_list = input_var_table(5:6),                                  &
     
    18452118
    18462119
    1847     FUNCTION init_io_group(in_files, out_vars, in_var_list, kind) RESULT(group)
     2120    FUNCTION init_io_group(in_files, out_vars, in_var_list, kind,              &
     2121                           n_output_quantities) RESULT(group)
    18482122       CHARACTER(LEN=PATH), INTENT(IN) ::  in_files(:)
    18492123       CHARACTER(LEN=*), INTENT(IN)    ::  kind
    18502124       TYPE(nc_var), INTENT(IN)        ::  out_vars(:)
    18512125       TYPE(nc_var), INTENT(IN)        ::  in_var_list(:)
     2126       INTEGER, OPTIONAL               ::  n_output_quantities
    18522127
    18532128       TYPE(io_group)                  ::  group
     
    18552130       group % nt = SIZE(in_files)
    18562131       group % nv = SIZE(out_vars)
     2132       group % n_inputs = SIZE(in_var_list)
    18572133       group % kind = TRIM(kind)
    18582134
    1859        ALLOCATE(group % in_var_list( SIZE(in_var_list) ))
     2135       IF ( PRESENT(n_output_quantities) )  THEN
     2136          group % n_output_quantities = n_output_quantities
     2137       ELSE
     2138          group % n_output_quantities = group % n_inputs
     2139       END IF
     2140
     2141       ALLOCATE(group % in_var_list(group % n_inputs))
    18602142       ALLOCATE(group % in_files(group % nt))
    18612143       ALLOCATE(group % out_vars(group % nv))
     
    18762158    SUBROUTINE fini_grids()
    18772159
    1878        CALL report('fini_grids', 'Deallocating grids')
     2160       CALL report('fini_grids', 'Deallocating grids', cfg % debug)
    18792161       
    18802162       DEALLOCATE(x, y, z, xu, yv, zw, z_column, zw_column)
     
    19042186    SUBROUTINE setup_variable_tables(ic_mode)
    19052187       CHARACTER(LEN=*), INTENT(IN) ::  ic_mode
     2188       INTEGER                      ::  n_invar = 0  !< number of variables in the input variable table
     2189       INTEGER                      ::  n_outvar = 0 !< number of variables in the output variable table
    19062190       TYPE(nc_var), POINTER        ::  var
    19072191
     
    19142198
    19152199       n_invar = 17
    1916        n_outvar = 55
     2200       n_outvar = 64
    19172201       ALLOCATE( input_var_table(n_invar) )
    19182202       ALLOCATE( output_var_table(n_outvar) )
     
    20522336       )
    20532337       IF (TRIM(ic_mode) == 'profile')  THEN
    2054           output_var_table(3) % grid => scalar_profile_grid
    2055           output_var_table(3) % intermediate_grid => scalar_profile_intermediate
     2338          output_var_table(3) % averaging_grid => averaged_scalar_profile
    20562339       END IF
    20572340
     
    21222405          units             = "kg/kg",                                         &
    21232406          kind              = "init scalar",                                   &
    2124           input_id          = 1,                                               &
     2407          input_id          = 3,                                               &
    21252408          output_file       = output_file,                                     &
    21262409          grid              = palm_grid,                                       &
     
    21292412       )
    21302413       IF (TRIM(ic_mode) == 'profile')  THEN
    2131           output_var_table(9) % grid => scalar_profile_grid
    2132           output_var_table(9) % intermediate_grid => scalar_profile_intermediate
     2414          output_var_table(9) % averaging_grid => averaged_scalar_profile
    21332415       END IF
    21342416
     
    21392421          units             = "kg/kg",                                         &
    21402422          kind              = "left scalar",                                   &
    2141           input_id          = 1,                                               &
     2423          input_id          = 3,                                               &
    21422424          output_file       = output_file,                                     &
    21432425          grid              = scalars_west_grid,                               &
     
    21512433          units             = "kg/kg",                                         &
    21522434          kind              = "right scalar",                                  &
    2153           input_id          = 1,                                               &
     2435          input_id          = 3,                                               &
    21542436          output_file       = output_file,                                     &
    21552437          grid              = scalars_east_grid,                               &
     
    21632445          units             = "kg/kg",                                         &
    21642446          kind              = "north scalar",                                  &
    2165           input_id          = 1,                                               &
     2447          input_id          = 3,                                               &
    21662448          output_file       = output_file,                                     &
    21672449          grid              = scalars_north_grid,                              &
     
    21752457          units             = "kg/kg",                                         &
    21762458          kind              = "south scalar",                                  &
    2177           input_id          = 1,                                               &
     2459          input_id          = 3,                                               &
    21782460          output_file       = output_file,                                     &
    21792461          grid              = scalars_south_grid,                              &
     
    21872469          units             = "kg/kg",                                         &
    21882470          kind              = "top scalar",                                    &
    2189           input_id          = 1,                                               &
     2471          input_id          = 3,                                               &
    21902472          output_file       = output_file,                                     &
    21912473          grid              = scalars_top_grid,                                &
     
    22062488       )
    22072489       IF (TRIM(ic_mode) == 'profile')  THEN
    2208           output_var_table(15) % grid => scalar_profile_grid
    2209           output_var_table(15) % intermediate_grid => scalar_profile_intermediate
     2490          output_var_table(15) % averaging_grid => averaged_scalar_profile
    22102491       END IF
    22112492
     
    22832564       )
    22842565       IF (TRIM(ic_mode) == 'profile')  THEN
    2285           output_var_table(21) % grid => scalar_profile_grid
    2286           output_var_table(21) % intermediate_grid => scalar_profile_intermediate
     2566          output_var_table(21) % averaging_grid => averaged_scalar_profile
    22872567       END IF
    22882568
     
    23602640       )
    23612641       IF (TRIM(ic_mode) == 'profile')  THEN
    2362           output_var_table(27) % grid => w_profile_grid
    2363           output_var_table(27) % intermediate_grid => w_profile_intermediate
     2642          output_var_table(27) % averaging_grid => averaged_w_profile
    23642643       END IF
    23652644
     
    25462825          intermediate_grid = palm_intermediate                                &
    25472826       )
     2827       output_var_table(42) % averaging_grid => averaged_scalar_profile
    25482828
    25492829       output_var_table(43) = init_nc_var(                                     &
     
    25522832          long_name         = "geostrophic wind (u component)",                &
    25532833          units             = "m/s",                                           &
    2554           kind              = "constant scalar profile",                       &
     2834          kind              = "geostrophic",                                   &
    25552835          input_id          = 1,                                               &
    25562836          output_file       = output_file,                                     &
    2557           grid              = scalar_profile_grid,                             &
    2558           intermediate_grid = scalar_profile_intermediate                      &
     2837          grid              = averaged_scalar_profile,                         &
     2838          intermediate_grid = averaged_scalar_profile                          &
    25592839       )
    25602840
     
    25642844          long_name         = "geostrophic wind (v component)",                &
    25652845          units             = "m/s",                                           &
    2566           kind              = "constant scalar profile",                       &
     2846          kind              = "geostrophic",                                   &
    25672847          input_id          = 1,                                               &
    25682848          output_file       = output_file,                                     &
    2569           grid              = scalar_profile_grid,                             &
    2570           intermediate_grid = scalar_profile_intermediate                      &
     2849          grid              = averaged_scalar_profile,                         &
     2850          intermediate_grid = averaged_scalar_profile                          &
    25712851       )
    25722852
     
    25762856          long_name         = "wind component in x direction",                 &
    25772857          units             = "m/s",                                           &
    2578           kind              = "large-scale scalar forcing",                    &
     2858          kind              = "geostrophic",                                   &
    25792859          input_id          = 1,                                               &
    25802860          output_file       = output_file,                                     &
    2581           grid              = scalar_profile_grid,                             &
    2582           intermediate_grid = scalar_profile_intermediate                      &
     2861          grid              = averaged_scalar_profile,                         &
     2862          intermediate_grid = averaged_scalar_profile                          &
    25832863       )
    25842864       output_var_table(45) % to_be_processed = ls_forcing_variables_required
     
    25922872          input_id          = 1,                                               &
    25932873          output_file       = output_file,                                     &
    2594           grid              = scalar_profile_grid,                             &
    2595           intermediate_grid = scalar_profile_intermediate                      &
     2874          grid              = averaged_scalar_profile,                         &
     2875          intermediate_grid = averaged_scalar_profile                          &
    25962876       )
    25972877       output_var_table(46) % to_be_processed = ls_forcing_variables_required
     
    26052885          input_id          = 1,                                               &
    26062886          output_file       = output_file,                                     &
    2607           grid              = w_profile_grid,                                  &
    2608           intermediate_grid = w_profile_intermediate                           &
     2887          grid              = averaged_scalar_profile,                         &
     2888          intermediate_grid = averaged_scalar_profile                          &
    26092889       )
    26102890       output_var_table(47) % to_be_processed = ls_forcing_variables_required
     
    26182898          input_id          = 1,                                               &
    26192899          output_file       = output_file,                                     &
    2620           grid              = w_profile_grid,                                  &
    2621           intermediate_grid = w_profile_intermediate                           &
     2900          grid              = averaged_w_profile,                              &
     2901          intermediate_grid = averaged_w_profile                               &
    26222902       )
    26232903       output_var_table(48) % to_be_processed = ls_forcing_variables_required
     
    26322912          input_id          = 1,                                               &
    26332913          output_file       = output_file,                                     &
    2634           grid              = scalar_profile_grid,                             &
    2635           intermediate_grid = scalar_profile_intermediate                      &
     2914          grid              = averaged_scalar_profile,                         &
     2915          intermediate_grid = averaged_scalar_profile                          &
    26362916       )
    26372917       output_var_table(49) % to_be_processed = ls_forcing_variables_required
     
    26452925          input_id          = 1,                                               &
    26462926          output_file       = output_file,                                     &
    2647           grid              = scalar_profile_grid,                             &
    2648           intermediate_grid = scalar_profile_intermediate                      &
     2927          grid              = averaged_scalar_profile,                         &
     2928          intermediate_grid = averaged_scalar_profile                          &
    26492929       )
    26502930       output_var_table(50) % to_be_processed = ls_forcing_variables_required
     
    26582938          input_id          = 1,                                               &
    26592939          output_file       = output_file,                                     &
    2660           grid              = scalar_profile_grid,                             &
    2661           intermediate_grid = scalar_profile_intermediate                      &
     2940          grid              = averaged_scalar_profile,                         &
     2941          intermediate_grid = averaged_scalar_profile                          &
    26622942       )
    26632943       output_var_table(51) % to_be_processed = ls_forcing_variables_required
     
    26692949          units             = "kg/kg/s",                                       &
    26702950          kind              = "large-scale scalar forcing",                    &
    2671           input_id          = 1,                                               &
    2672           output_file       = output_file,                                     &
    2673           grid              = scalar_profile_grid,                             &
    2674           intermediate_grid = scalar_profile_intermediate                      &
     2951          input_id          = 3,                                               &
     2952          output_file       = output_file,                                     &
     2953          grid              = averaged_scalar_profile,                         &
     2954          intermediate_grid = averaged_scalar_profile                          &
    26752955       )
    26762956       output_var_table(52) % to_be_processed = ls_forcing_variables_required
     
    26832963          units             = "kg/kg/s",                                       &
    26842964          kind              = "large-scale scalar forcing",                    &
    2685           input_id          = 1,                                               &
    2686           output_file       = output_file,                                     &
    2687           grid              = scalar_profile_grid,                             &
    2688           intermediate_grid = scalar_profile_intermediate                      &
     2965          input_id          = 3,                                               &
     2966          output_file       = output_file,                                     &
     2967          grid              = averaged_scalar_profile,                         &
     2968          intermediate_grid = averaged_scalar_profile                          &
    26892969       )
    26902970       output_var_table(53) % to_be_processed = ls_forcing_variables_required
     
    26962976          units             = "kg/kg",                                         &
    26972977          kind              = "large-scale scalar forcing",                    &
    2698           input_id          = 1,                                               &
    2699           output_file       = output_file,                                     &
    2700           grid              = scalar_profile_grid,                             &
    2701           intermediate_grid = scalar_profile_intermediate                      &
     2978          input_id          = 3,                                               &
     2979          output_file       = output_file,                                     &
     2980          grid              = averaged_scalar_profile,                         &
     2981          intermediate_grid = averaged_scalar_profile                          &
    27022982       )
    27032983       output_var_table(54) % to_be_processed = ls_forcing_variables_required
     
    27112991          input_id          = 1,                                               &
    27122992          output_file       = output_file,                                     &
    2713           grid              = scalar_profile_grid,                             &
    2714           intermediate_grid = scalar_profile_intermediate                      &
     2993          grid              = averaged_scalar_profile,                         &
     2994          intermediate_grid = averaged_scalar_profile                          &
    27152995       )
    27162996       output_var_table(55) % to_be_processed = ls_forcing_variables_required
    27172997
     2998
     2999       output_var_table(56) = init_nc_var(                                     &
     3000          name              = 'internal_density_centre',                              &
     3001          std_name          = "",                                              &
     3002          long_name         = "",                                              &
     3003          units             = "",                                              &
     3004          kind              = "internal profile",                              &
     3005          input_id          = 4,                                               &
     3006          output_file       = output_file,                                     &
     3007          grid              = averaged_scalar_profile,                         &
     3008          intermediate_grid = averaged_scalar_profile                          &
     3009       )
     3010       output_var_table(56) % averaging_grid => averaged_scalar_profile
     3011
     3012
     3013       output_var_table(57) = init_nc_var(                                     &
     3014          name              = 'internal_density_north',                       &
     3015          std_name          = "",                                              &
     3016          long_name         = "",                                              &
     3017          units             = "",                                              &
     3018          kind              = "internal profile",                              &
     3019          input_id          = 4,                                               &
     3020          output_file       = output_file,                                     &
     3021          grid              = north_averaged_scalar_profile,                   &
     3022          intermediate_grid = north_averaged_scalar_profile                    &
     3023       )
     3024       output_var_table(57) % averaging_grid => north_averaged_scalar_profile
     3025       output_var_table(57) % to_be_processed = .NOT. cfg % ug_is_set
     3026
     3027
     3028       output_var_table(58) = init_nc_var(                                     &
     3029          name              = 'internal_density_south',                       &
     3030          std_name          = "",                                              &
     3031          long_name         = "",                                              &
     3032          units             = "",                                              &
     3033          kind              = "internal profile",                              &
     3034          input_id          = 4,                                               &
     3035          output_file       = output_file,                                     &
     3036          grid              = south_averaged_scalar_profile,                   &
     3037          intermediate_grid = south_averaged_scalar_profile                    &
     3038       )
     3039       output_var_table(58) % averaging_grid => south_averaged_scalar_profile
     3040       output_var_table(58) % to_be_processed = .NOT. cfg % ug_is_set
     3041
     3042
     3043       output_var_table(59) = init_nc_var(                                     &
     3044          name              = 'internal_density_east',                        &
     3045          std_name          = "",                                              &
     3046          long_name         = "",                                              &
     3047          units             = "",                                              &
     3048          kind              = "internal profile",                              &
     3049          input_id          = 4,                                               &
     3050          output_file       = output_file,                                     &
     3051          grid              = east_averaged_scalar_profile,                    &
     3052          intermediate_grid = east_averaged_scalar_profile                     &
     3053       )
     3054       output_var_table(59) % averaging_grid => east_averaged_scalar_profile
     3055       output_var_table(59) % to_be_processed = .NOT. cfg % ug_is_set
     3056
     3057
     3058       output_var_table(60) = init_nc_var(                                     &
     3059          name              = 'internal_density_west',                        &
     3060          std_name          = "",                                              &
     3061          long_name         = "",                                              &
     3062          units             = "",                                              &
     3063          kind              = "internal profile",                              &
     3064          input_id          = 4,                                               &
     3065          output_file       = output_file,                                     &
     3066          grid              = west_averaged_scalar_profile,                    &
     3067          intermediate_grid = west_averaged_scalar_profile                     &
     3068       )
     3069       output_var_table(60) % averaging_grid => west_averaged_scalar_profile
     3070       output_var_table(60) % to_be_processed = .NOT. cfg % ug_is_set
     3071
     3072       output_var_table(61) = init_nc_var(                                     &
     3073          name              = 'internal_pressure_north',                       &
     3074          std_name          = "",                                              &
     3075          long_name         = "",                                              &
     3076          units             = "",                                              &
     3077          kind              = "internal profile",                              &
     3078          input_id          = 2,                                               &
     3079          output_file       = output_file,                                     &
     3080          grid              = north_averaged_scalar_profile,                   &
     3081          intermediate_grid = north_averaged_scalar_profile                    &
     3082       )
     3083       output_var_table(61) % averaging_grid => north_averaged_scalar_profile
     3084       output_var_table(61) % to_be_processed = .NOT. cfg % ug_is_set
     3085
     3086
     3087       output_var_table(62) = init_nc_var(                                     &
     3088          name              = 'internal_pressure_south',                       &
     3089          std_name          = "",                                              &
     3090          long_name         = "",                                              &
     3091          units             = "",                                              &
     3092          kind              = "internal profile",                              &
     3093          input_id          = 2,                                               &
     3094          output_file       = output_file,                                     &
     3095          grid              = south_averaged_scalar_profile,                   &
     3096          intermediate_grid = south_averaged_scalar_profile                    &
     3097       )
     3098       output_var_table(62) % averaging_grid => south_averaged_scalar_profile
     3099       output_var_table(62) % to_be_processed = .NOT. cfg % ug_is_set
     3100
     3101
     3102       output_var_table(63) = init_nc_var(                                     &
     3103          name              = 'internal_pressure_east',                        &
     3104          std_name          = "",                                              &
     3105          long_name         = "",                                              &
     3106          units             = "",                                              &
     3107          kind              = "internal profile",                              &
     3108          input_id          = 2,                                               &
     3109          output_file       = output_file,                                     &
     3110          grid              = east_averaged_scalar_profile,                    &
     3111          intermediate_grid = east_averaged_scalar_profile                     &
     3112       )
     3113       output_var_table(63) % averaging_grid => east_averaged_scalar_profile
     3114       output_var_table(63) % to_be_processed = .NOT. cfg % ug_is_set
     3115
     3116
     3117       output_var_table(64) = init_nc_var(                                     &
     3118          name              = 'internal_pressure_west',                        &
     3119          std_name          = "",                                              &
     3120          long_name         = "",                                              &
     3121          units             = "",                                              &
     3122          kind              = "internal profile",                              &
     3123          input_id          = 2,                                               &
     3124          output_file       = output_file,                                     &
     3125          grid              = west_averaged_scalar_profile,                    &
     3126          intermediate_grid = west_averaged_scalar_profile                     &
     3127       )
     3128       output_var_table(64) % averaging_grid => west_averaged_scalar_profile
     3129       output_var_table(64) % to_be_processed = .NOT. cfg % ug_is_set
    27183130
    27193131       ! Attributes shared among all variables
     
    27323144!------------------------------------------------------------------------------!
    27333145    FUNCTION init_nc_var(name, std_name, long_name, units, kind, input_id,     &
    2734                          grid, intermediate_grid, output_file, is_profile      &
    2735              ) RESULT(var)
     3146                         grid, intermediate_grid, output_file, is_profile)     &
     3147       RESULT(var)
    27363148
    27373149       CHARACTER(LEN=*), INTENT(IN)      ::  name, std_name, long_name, units, kind
     
    27723184          var % dimvarids(1:3)  = output_file % dimvarids_soil
    27733185          var % to_be_processed = init_variables_required
     3186          var % is_internal     = .FALSE.
    27743187          var % task            = "interpolate_2d"
    27753188
     
    27813194          var % dimvarids(1:3)  = output_file % dimvarids_scl
    27823195          var % to_be_processed = init_variables_required
     3196          var % is_internal     = .FALSE.
    27833197          var % task            = "interpolate_3d"
    27843198
     
    27943208          var % dimvarids(3)    = output_file % dimvarids_scl(3)
    27953209          var % to_be_processed = init_variables_required
     3210          var % is_internal     = .FALSE.
    27963211          var % task            = "interpolate_3d"
    27973212
     
    28073222          var % dimvarids(3)    = output_file % dimvarids_scl(3)
    28083223          var % to_be_processed = init_variables_required
     3224          var % is_internal     = .FALSE.
    28093225          var % task            = "interpolate_3d"
    28103226
     
    28203236          var % dimvarids(3)    = output_file % dimvarids_vel(3)
    28213237          var % to_be_processed = init_variables_required
     3238          var % is_internal     = .FALSE.
    28223239          var % task            = "interpolate_3d"
    28233240
     
    28293246          var % dimvarids(1)    = output_file % dimvarids_scl(3) !z
    28303247          var % to_be_processed = init_variables_required
     3248          var % is_internal     = .FALSE.
    28313249          var % task            = "average profile"
    28323250
     
    28383256          var % dimvarids(1)    = output_file % dimvarids_vel(3) !z
    28393257          var % to_be_processed = init_variables_required
     3258          var % is_internal     = .FALSE.
    28403259          var % task            = "average profile"
    28413260
     
    28483267          var % dimvarids(1:2)  = output_file % dimvarids_soil(1:2)
    28493268          var % to_be_processed = boundary_variables_required
     3269          var % is_internal     = .FALSE.
    28503270          var % task            = "interpolate_2d"
    28513271
     
    28603280          var % dimvarids(2)    = output_file % dimvarids_scl(3)
    28613281          var % to_be_processed = boundary_variables_required
     3282          var % is_internal     = .FALSE.
    28623283          var % task            = "interpolate_3d"
    28633284
     
    28723293          var % dimvarids(2)    = output_file % dimvarids_scl(3)
    28733294          var % to_be_processed = boundary_variables_required
     3295          var % is_internal     = .FALSE.
    28743296          var % task            = "interpolate_3d"
    28753297
     
    28843306          var % dimvarids(2)    = output_file % dimvarids_scl(2)
    28853307          var % to_be_processed = boundary_variables_required
     3308          var % is_internal     = .FALSE.
    28863309          var % task            = "interpolate_3d"
    28873310
     
    28963319          var % dimvarids(2)    = output_file % dimvarids_scl(3)
    28973320          var % to_be_processed = boundary_variables_required
     3321          var % is_internal     = .FALSE.
    28983322          var % task            = "interpolate_3d"
    28993323
     
    29083332          var % dimvarids(2)    = output_file % dimvarids_scl(3)
    29093333          var % to_be_processed = boundary_variables_required
     3334          var % is_internal     = .FALSE.
    29103335          var % task            = "interpolate_3d"
    29113336
     
    29203345          var % dimvarids(2)    = output_file % dimvarids_scl(2)
    29213346          var % to_be_processed = boundary_variables_required
     3347          var % is_internal     = .FALSE.
    29223348          var % task            = "interpolate_3d"
    29233349
     
    29323358          var % dimvarids(2)    = output_file % dimvarids_scl(3)
    29333359          var % to_be_processed = boundary_variables_required
     3360          var % is_internal     = .FALSE.
    29343361          var % task            = "interpolate_3d"
    29353362
     
    29443371          var % dimvarids(2)    = output_file % dimvarids_scl(3)
    29453372          var % to_be_processed = boundary_variables_required
     3373          var % is_internal     = .FALSE.
    29463374          var % task            = "interpolate_3d"
    29473375
     
    29563384          var % dimvarids(2)    = output_file % dimvarids_vel(2)
    29573385          var % to_be_processed = boundary_variables_required
     3386          var % is_internal     = .FALSE.
    29583387          var % task            = "interpolate_3d"
    29593388
     
    29683397          var % dimvarids(2)    = output_file % dimvarids_vel(3)
    29693398          var % to_be_processed = boundary_variables_required
     3399          var % is_internal     = .FALSE.
    29703400          var % task            = "interpolate_3d"
    29713401
     
    29803410          var % dimvarids(2)    = output_file % dimvarids_vel(3)
    29813411          var % to_be_processed = boundary_variables_required
     3412          var % is_internal     = .FALSE.
    29823413          var % task            = "interpolate_3d"
    29833414
     
    29883419          var % dimvarids(1)    = output_file % dimvarid_time
    29893420          var % to_be_processed = .TRUE.
    2990           var % task            = "average scalar"
     3421          var % is_internal     = .FALSE.
     3422          var % task            = "average profile"
    29913423
    29923424       CASE( 'constant scalar profile' )
     
    29983430          var % dimvarids(1)    = output_file % dimvarids_scl(3)
    29993431          var % to_be_processed = .TRUE.
     3432          var % is_internal     = .FALSE.
    30003433          var % task            = "set profile"
    30013434
     
    30083441          var % dimvarids(1)    = output_file % dimvarids_scl(3)
    30093442          var % to_be_processed = ls_forcing_variables_required
     3443          var % is_internal     = .FALSE.
    30103444          var % task            = "average large-scale profile"
     3445
     3446       CASE( 'geostrophic' )
     3447          var % lod             = -1
     3448          var % ndim            = 2
     3449          var % dimids(2)       = output_file % dimid_time    !t
     3450          var % dimids(1)       = output_file % dimids_scl(3) !z
     3451          var % dimvarids(2)    = output_file % dimvarid_time
     3452          var % dimvarids(1)    = output_file % dimvarids_scl(3)
     3453          var % to_be_processed = .TRUE.
     3454          var % is_internal     = .FALSE.
     3455          var % task            = "geostrophic winds"
    30113456
    30123457       CASE( 'large-scale w forcing' )
     
    30183463          var % dimvarids(1)    = output_file % dimvarids_vel(3)
    30193464          var % to_be_processed = ls_forcing_variables_required
     3465          var % is_internal     = .FALSE.
    30203466          var % task            = "average large-scale profile"
     3467
     3468       CASE( 'internal profile' )
     3469          var % lod             = -1
     3470          var % ndim            = 2
     3471          var % dimids(2)       = output_file % dimid_time    !t
     3472          var % dimids(1)       = output_file % dimids_scl(3) !z
     3473          var % dimvarids(2)    = output_file % dimvarid_time
     3474          var % dimvarids(1)    = output_file % dimvarids_scl(3)
     3475          var % to_be_processed = .TRUE.
     3476          var % is_internal     = .TRUE.
     3477          var % task            = "internal profile"
    30213478
    30223479       CASE DEFAULT
     
    30313488    SUBROUTINE fini_variables()
    30323489
    3033        CALL report('fini_variables', 'Deallocating variable table')
     3490       CALL report('fini_variables', 'Deallocating variable table', cfg % debug)
    30343491       DEALLOCATE( input_var_table )
    30353492
     
    30393496    SUBROUTINE fini_io_groups()
    30403497
    3041        CALL report('fini_io_groups', 'Deallocating IO groups')
     3498       CALL report('fini_io_groups', 'Deallocating IO groups', cfg % debug)
    30423499       DEALLOCATE( io_group_list )
    30433500
     
    30473504    SUBROUTINE fini_file_lists()
    30483505       
    3049        CALL report('fini_file_lists', 'Deallocating file lists')
     3506       CALL report('fini_file_lists', 'Deallocating file lists', cfg % debug)
    30503507       DEALLOCATE( flow_files, soil_files, radiation_files, soil_moisture_files )
    30513508
     
    31343591             ! rotated velocities
    31353592             DO k = 1, nz
    3136              DO j = 2, ny
    3137              DO i = 2, nx
     3593             DO j = 1, ny
     3594             DO i = 1, nx
    31383595                CALL uv2uvrot( urot = preprocess_buffer(1) % array(i,j,k),     &
    31393596                               vrot = preprocess_buffer(2) % array(i,j,k),     &
     
    31713628          CALL report('preprocess', message)
    31723629       
    3173        CASE( 'temperature' ) ! P and T
     3630       CASE( 'thermodynamics' ) ! T, P, QV
    31743631          nx = SIZE(input_buffer(1) % array, 1)
    31753632          ny = SIZE(input_buffer(1) % array, 2)
     
    31923649
    31933650                input_buffer (2) % array(i,j,:) =                              &
    3194                    input_buffer (2) % array(i,j,:) + basic_state_pressure(:)
     3651                   HECTO * input_buffer (2) % array(i,j,:) +                   &
     3652                   basic_state_pressure(:)
    31953653
    31963654             END DO
     
    32043662
    32053663          END IF
     3664          ! pressure
     3665          input_buffer(2) % is_preprocessed = .TRUE.
     3666
     3667          ! Copy temperature to last input buffer array
     3668          ALLOCATE(                                                            &
     3669              input_buffer( group % n_output_quantities ) % array (nx, ny, nz) &
     3670          )
     3671          input_buffer(group % n_output_quantities) % array(:,:,:) =           &
     3672              input_buffer(1) % array(:,:,:)
    32063673
    32073674          ! Convert absolute to potential temperature
    3208           input_buffer(1) % array(:,:,:) = input_buffer(1) % array(:,:,:) *    &
    3209              EXP( RD_PALM/CP_PALM * LOG( P_REF / input_buffer(2) % array(:,:,:) ))
    3210 
    3211           input_buffer(1:2) % is_preprocessed = .TRUE.
     3675          CALL potential_temperature(                                          &
     3676             t = input_buffer(1) % array(:,:,:),                               &
     3677             p = input_buffer(2) % array(:,:,:),                               &
     3678             p_ref = P_REF,                                                    &
     3679             r = RD_PALM,                                                      &
     3680             cp = CP_PALM                                                      &
     3681          )
     3682
     3683          ! potential temperature
     3684          input_buffer(1) % is_preprocessed = .TRUE.
     3685
     3686          ! Convert temperature copy to density
     3687          CALL moist_density(                                                  &
     3688             t_rho = input_buffer(group % n_output_quantities) % array(:,:,:), &
     3689             p = input_buffer(2) % array(:,:,:),                               &
     3690             qv = input_buffer(3) % array(:,:,:),                              &
     3691             rd = RD,                                                          &
     3692             rv = RV                                                           &
     3693          )
     3694
     3695          ! qv
     3696          input_buffer(3) % is_preprocessed = .TRUE.
     3697
     3698          ! density
     3699          input_buffer(group % n_output_quantities) % is_preprocessed = .TRUE.
     3700
     3701
    32123702          message = "Input buffers for group '" // TRIM(group % kind) // "'"//&
    32133703             " preprocessed sucessfully."
     
    33363826
    33373827       CASE DEFAULT
    3338           message = "Buffer kind '" // TRIM(group % kind) // "' is not supported."
     3828          message = "IO group kind '" // TRIM(group % kind) // "' is not supported."
    33393829          CALL abort('prerpocess', message)
    33403830
  • palm/trunk/UTIL/inifor/src/inifor.f90

    r3183 r3395  
    2626! -----------------
    2727! $Id$
     28! Added main loop support for computation of geostrophic winds and surface
     29!     pressure
     30! Cleaned up terminal output, show some meVssages only with --debug option
     31!
     32! 3183 2018-07-27 14:25:55Z suehring
    2833! Introduced new PALM grid stretching
    2934! Renamend initial-condition mode variable 'mode' to 'ic_mode'
     
    5661               fini_file_lists, preprocess, origin_lon, origin_lat,            &
    5762               output_file, io_group_list, output_var_table,                   &
    58                cosmo_grid, palm_grid, nx, ny, nz, ug, vg, p0, cfg,             &
    59                average_imin, average_imax, average_jmin, average_jmax
    60 
     63               cosmo_grid, palm_grid, nx, ny, nz, p0, cfg, f3,                 &
     64               averaging_width_ns, averaging_width_ew, phi_n, lambda_n,        &
     65               lam_centre, phi_centre
    6166    USE io
    6267    USE transform,                                                             &
    63         ONLY:  average_profile, average_2d, interpolate_2d, interpolate_3d
     68        ONLY:  average_profile, interpolate_2d, interpolate_3d,                &
     69               geostrophic_winds, extrapolate_density, extrapolate_pressure,   &
     70               get_surface_pressure
    6471    USE types
    6572   
     
    7077    INTEGER                                 ::  iter
    7178    REAL(dp), ALLOCATABLE, DIMENSION(:,:,:) ::  output_arr
     79    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_centre
     80    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  ug_arr
     81    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  vg_arr
     82    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_north
     83    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_south
     84    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_east
     85    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  rho_west
     86    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_north
     87    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_south
     88    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_east
     89    REAL(dp), ALLOCATABLE, DIMENSION(:), TARGET ::  p_west
     90    REAL(dp), POINTER, DIMENSION(:)         ::  internal_arr
     91    REAL(dp), POINTER, DIMENSION(:)         ::  ug_vg_arr
    7292    TYPE(nc_var), POINTER                   ::  output_var
    7393    TYPE(io_group), POINTER                 ::  group
    7494    TYPE(container), ALLOCATABLE            ::  input_buffer(:)
     95    LOGICAL, SAVE                           ::  ug_vg_have_been_computed = .FALSE.
     96    LOGICAL, SAVE                           ::  debugging_output = .TRUE.
    7597   
    7698!> \mainpage About INIFOR
     
    199221 CALL run_control('time', 'alloc')
    200222                     
    201 
    202223                      CALL average_profile(                                    &
    203224                         input_buffer(output_var % input_id) % array(:,:,:),   &
    204                          output_arr(:,:,:), average_imin, average_imax,        &
    205                          average_jmin, average_jmax,                           &
    206                          output_var % intermediate_grid,                       &
    207                          output_var % grid)
    208  CALL run_control('time', 'comp')
     225                         output_arr(0,0,:),                                    &
     226                         output_var % averaging_grid)
     227
     228                      IF ( TRIM(output_var % name) ==                          &
     229                           'surface_forcing_surface_pressure' )  THEN
     230
     231                         IF ( cfg % p0_is_set )  THEN
     232                            output_arr(0,0,1) = p0
     233                         ELSE
     234                            CALL get_surface_pressure(                         &
     235                               output_arr(0,0,:), rho_centre,                  &
     236                               output_var % averaging_grid)
     237                         END IF
     238
     239                      END IF
     240 CALL run_control('time', 'comp')
     241
     242                   CASE ( 'internal profile' )
     243
     244                      message = "Averaging of internal profile for variable '" //&
     245                         TRIM(output_var % name) // "' is not supported."
     246
     247                      SELECT CASE (TRIM(output_var % name))
     248
     249                      CASE('internal_density_centre')
     250                         ALLOCATE( rho_centre( 1:output_var % grid % nz ) )
     251                         internal_arr => rho_centre
     252
     253                      CASE('internal_density_north')
     254                         ALLOCATE( rho_north( 1:output_var % grid % nz ) )
     255                         internal_arr => rho_north
     256
     257                      CASE('internal_density_south')
     258                         ALLOCATE( rho_south( 1:output_var % grid % nz ) )
     259                         internal_arr => rho_south
     260
     261                      CASE('internal_density_east')
     262                         ALLOCATE( rho_east( 1:output_var % grid % nz) )
     263                         internal_arr => rho_east
     264
     265                      CASE('internal_density_west')
     266                         ALLOCATE( rho_west( 1:output_var % grid % nz ) )
     267                         internal_arr => rho_west
     268
     269                      CASE('internal_pressure_north')
     270                         ALLOCATE( p_north( 1:output_var % grid % nz ) )
     271                         internal_arr => p_north
     272
     273                      CASE('internal_pressure_south')
     274                         ALLOCATE( p_south( 1:output_var % grid % nz ) )
     275                         internal_arr => p_south
     276
     277                      CASE('internal_pressure_east')
     278                         ALLOCATE( p_east( 1:output_var % grid % nz) )
     279                         internal_arr => p_east
     280
     281                      CASE('internal_pressure_west')
     282                         ALLOCATE( p_west( 1:output_var % grid % nz ) )
     283                         internal_arr => p_west
     284
     285                      CASE DEFAULT
     286                         CALL abort('main loop', message)
     287
     288                      END SELECT
     289 CALL run_control('time', 'alloc')
     290
     291
     292                      CALL average_profile(                                 &
     293                         input_buffer(output_var % input_id) % array(:,:,:),&
     294                         internal_arr(:),                                   &
     295                         output_var % averaging_grid)
     296
     297                      SELECT CASE (TRIM(output_var % name))
     298
     299                      CASE('internal_density_centre',                          &
     300                           'internal_density_north',                           &
     301                           'internal_density_south',                           &
     302                           'internal_density_east',                            &
     303                           'internal_density_west')
     304                         CALL extrapolate_density(internal_arr,                &
     305                                                  output_var % averaging_grid)
     306
     307                      CASE('internal_pressure_north')
     308                         CALL extrapolate_pressure(internal_arr, rho_north,    &
     309                                                   output_var % averaging_grid)
     310
     311                      CASE('internal_pressure_south')
     312                         CALL extrapolate_pressure(internal_arr, rho_south,    &
     313                                                   output_var % averaging_grid)
     314
     315                      CASE('internal_pressure_east')
     316                         CALL extrapolate_pressure(internal_arr, rho_east,     &
     317                                                   output_var % averaging_grid)
     318
     319                      CASE('internal_pressure_west')
     320                         CALL extrapolate_pressure(internal_arr, rho_west,     &
     321                                                   output_var % averaging_grid)
     322
     323                      CASE DEFAULT
     324                         CALL abort('main loop', message)
     325
     326                      END SELECT
     327
     328                      IF (.TRUE.)  THEN
     329                         ALLOCATE( output_arr(1,1,1:output_var % grid % nz) )
     330                         output_arr(1,1,:) = internal_arr(:)
     331                      END IF
     332 CALL run_control('time', 'comp')
     333
     334
     335                   ! This case gets called twice, the first time for ug, the
     336                   ! second time for vg. We compute ug and vg at the first call
     337                   ! and keep vg (and ug for that matter) around for the second
     338                   ! call.
     339                   CASE ( 'geostrophic winds' )
     340
     341                      IF (.NOT. ug_vg_have_been_computed )  THEN
     342                         ALLOCATE( ug_arr(output_var % grid % nz) )
     343                         ALLOCATE( vg_arr(output_var % grid % nz) )
     344
     345                         IF ( cfg % ug_is_set )  THEN
     346                            ug_arr = cfg % ug
     347                            vg_arr = cfg % vg
     348                         ELSE
     349                            CALL geostrophic_winds( p_north, p_south, p_east,     &
     350                                                    p_west, rho_centre, f3,      &
     351                                                    averaging_width_ew,           &
     352                                                    averaging_width_ns,           &
     353                                                    phi_n, lambda_n, phi_centre, lam_centre, &
     354                                                    ug_arr, vg_arr )
     355                         END IF
     356
     357                         ug_vg_have_been_computed = .TRUE.
     358
     359                      END IF
     360
     361                      ! Prepare output of geostrophic winds
     362                      SELECT CASE(TRIM(output_var % name))
     363                      CASE ('ls_forcing_ug')
     364                         ug_vg_arr => ug_arr
     365                      CASE ('ls_forcing_vg')
     366                         ug_vg_arr => vg_arr
     367                      END SELECT
     368
     369                      ALLOCATE( output_arr(1,1,output_var % grid % nz) )
     370                      output_arr(1,1,:) = ug_vg_arr(:)
    209371
    210372                   CASE ( 'average scalar' )
     
    222384                      SELECT CASE (TRIM(output_var % name))
    223385
    224                       CASE('ls_forcing_ug')
    225                           output_arr(1, 1, :) = ug
    226 
    227                       CASE('ls_forcing_vg')
    228                           output_arr(1, 1, :) = vg
     386                      !CASE('ls_forcing_ug')
     387                      !    output_arr(1, 1, :) = ug
     388
     389                      !CASE('ls_forcing_vg')
     390                      !    output_arr(1, 1, :) = vg
    229391
    230392                      CASE('nudging_tau')
     
    256418!- Section 2.3: Write current time step of current variable
    257419!------------------------------------------------------------------------------
    258                    message = "Writing variable '" // TRIM(output_var%name) // "'."
    259                    CALL report('main loop', message)
    260                    CALL update_output(output_var, output_arr, iter, output_file)
     420                   IF (.NOT. output_var % is_internal .OR. debugging_output)  THEN
     421                      message = "Writing variable '" // TRIM(output_var%name) // "'."
     422                      CALL report('main loop', message)
     423                      CALL update_output(output_var, output_arr, iter, output_file)
    261424 CALL run_control('time', 'write')
    262 
    263                    DEALLOCATE(output_arr)
     425                   END IF
     426
     427                   IF (ALLOCATED(output_arr))  DEALLOCATE(output_arr)
    264428 CALL run_control('time', 'alloc')
    265429
     
    267431
    268432             END DO ! ouput variables
     433
     434             IF ( group % kind == 'thermodynamics' )  THEN
     435                DEALLOCATE( rho_centre )
     436                IF ( .NOT. cfg % ug_is_set )  THEN
     437                   DEALLOCATE( rho_north )
     438                   DEALLOCATE( rho_south )
     439                   DEALLOCATE( rho_east )
     440                   DEALLOCATE( rho_west )
     441                   DEALLOCATE( p_north )
     442                   DEALLOCATE( p_south )
     443                   DEALLOCATE( p_east )
     444                   DEALLOCATE( p_west )
     445                END IF
     446             END IF
    269447
    270448             IF ( group % kind == 'running average' .OR. &
     
    273451                ! accumulated COSMO-DE quantities (precipitation).
    274452             ELSE
    275                 CALL report('main loop', 'Deallocating input buffer')
     453                CALL report('main loop', 'Deallocating input buffer', cfg % debug)
    276454                DEALLOCATE(input_buffer)
    277455             END IF
     
    281459
    282460          IF (ALLOCATED(input_buffer))  THEN
    283              CALL report('main loop', 'Deallocating input buffer')
     461             CALL report('main loop', 'Deallocating input buffer', cfg % debug)
    284462             DEALLOCATE(input_buffer)
    285463          END IF
     
    294472          END IF
    295473
    296           CALL report('main loop', message)
     474          CALL report('main loop', message, cfg % debug)
    297475
    298476       END IF ! IO group % to_be_processed
  • palm/trunk/UTIL/inifor/src/io.f90

    r3262 r3395  
    2626! -----------------
    2727! $Id$
    28 ! Removed unnecessary check for output file
     28! Added command-line options for configuring the computation of geostrophic
     29!     winds (--averagin-mode, --averaging-angle)
     30! Added command-line option --input-prefix for setting input file prefixes all
     31!     at once
     32! Added --debug option for more verbose terminal output
     33! Added option-specific *_is_set LOGICALs to indicate invocation from the
     34!     command-line
     35! Improved error messages in case of empty file-name strings
     36! Improved routine naming
    2937!
    3038! 3183 2018-07-27 14:25:55Z suehring
     
    148156       INTEGER                            ::  arg_count, i
    149157
     158       cfg % p0_is_set = .FALSE.
     159       cfg % ug_is_set = .FALSE.
     160       cfg % vg_is_set = .FALSE.
     161       cfg % flow_prefix_is_set = .FALSE.
     162       cfg % input_prefix_is_set = .FALSE.
     163       cfg % radiation_prefix_is_set = .FALSE.
     164       cfg % soil_prefix_is_set = .FALSE.
     165       cfg % soilmoisture_prefix_is_set = .FALSE.
     166
    150167       arg_count = COMMAND_ARGUMENT_COUNT()
    151168       IF (arg_count .GT. 0)  THEN
     
    170187             SELECT CASE( TRIM(option) )
    171188
     189             CASE( '--averaging-mode' )
     190                CALL get_option_argument( i, arg )
     191                cfg % averaging_mode = TRIM(arg)
     192
    172193             CASE( '-date', '-d', '--date' )
    173194                CALL get_option_argument( i, arg )
    174195                cfg % start_date = TRIM(arg)
     196
     197             CASE( '--debug' )
     198                cfg % debug = .TRUE.
    175199
    176200             ! Elevation of the PALM-4U domain above sea level
     
    181205             ! surface pressure, at z0
    182206             CASE( '-p0', '-r', '--surface-pressure' )
     207                cfg % p0_is_set = .TRUE.
    183208                CALL get_option_argument( i, arg )
    184209                READ(arg, *) cfg % p0
     
    186211             ! geostrophic wind in x direction
    187212             CASE( '-ug', '-u', '--geostrophic-u' )
     213                cfg % ug_is_set = .TRUE.
    188214                CALL get_option_argument( i, arg )
    189215                READ(arg, *) cfg % ug
     
    191217             ! geostrophic wind in y direction
    192218             CASE( '-vg', '-v', '--geostrophic-v' )
     219                cfg % vg_is_set = .TRUE.
    193220                CALL get_option_argument( i, arg )
    194221                READ(arg, *) cfg % vg
     
    208235             CASE( '-hhl', '-l', '--hhl-file' )
    209236                CALL get_option_argument( i, arg )
    210                  cfg % hhl_file = TRIM(arg)
     237                cfg % hhl_file = TRIM(arg)
     238
     239             CASE( '--input-prefix')
     240                CALL get_option_argument( i, arg )
     241                cfg % input_prefix = TRIM(arg)
     242                cfg % input_prefix_is_set = .TRUE.
     243   
     244             CASE( '-a', '--averaging-angle' )
     245                CALL get_option_argument( i, arg )
     246                READ(arg, *) cfg % averaging_angle
    211247
    212248             CASE( '-static', '-t', '--static-driver' )
    213249                CALL get_option_argument( i, arg )
    214                  cfg % static_driver_file = TRIM(arg)
     250                cfg % static_driver_file = TRIM(arg)
    215251
    216252             CASE( '-soil', '-s', '--soil-file')
    217253                CALL get_option_argument( i, arg )
    218                  cfg % soiltyp_file = TRIM(arg)
     254                cfg % soiltyp_file = TRIM(arg)
    219255
    220256             CASE( '--flow-prefix')
    221257                CALL get_option_argument( i, arg )
    222                  cfg % flow_prefix = TRIM(arg)
    223 
     258                cfg % flow_prefix = TRIM(arg)
     259                cfg % flow_prefix_is_set = .TRUE.
     260   
    224261             CASE( '--radiation-prefix')
    225262                CALL get_option_argument( i, arg )
    226                  cfg % radiation_prefix = TRIM(arg)
    227 
     263                cfg % radiation_prefix = TRIM(arg)
     264                cfg % radiation_prefix_is_set = .TRUE.
     265   
    228266             CASE( '--soil-prefix')
    229267                CALL get_option_argument( i, arg )
    230                  cfg % soil_prefix = TRIM(arg)
    231 
     268                cfg % soil_prefix = TRIM(arg)
     269                cfg % soil_prefix_is_set = .TRUE.
     270   
    232271             CASE( '--soilmoisture-prefix')
    233272                CALL get_option_argument( i, arg )
    234                  cfg % soilmoisture_prefix = TRIM(arg)
     273                cfg % soilmoisture_prefix = TRIM(arg)
     274                cfg % soilmoisture_prefix_is_set = .TRUE.
    235275
    236276             CASE( '-o', '--output' )
     
    297337
    298338      all_files_present = .TRUE.
    299       all_files_present = all_files_present .AND. file_present(cfg % hhl_file)
    300       all_files_present = all_files_present .AND. file_present(cfg % namelist_file)
    301       all_files_present = all_files_present .AND. file_present(cfg % soiltyp_file)
     339      all_files_present = all_files_present .AND. file_present(cfg % hhl_file, 'HHL')
     340      all_files_present = all_files_present .AND. file_present(cfg % namelist_file, 'NAMELIST')
     341      all_files_present = all_files_present .AND. file_present(cfg % soiltyp_file, 'SOILTYP')
    302342
    303343      ! Only check optional static driver file name, if it has been given.
    304344      IF (TRIM(cfg % static_driver_file) .NE. '')  THEN
    305          all_files_present = all_files_present .AND. file_present(cfg % static_driver_file)
     345         all_files_present = all_files_present .AND. file_present(cfg % static_driver_file, 'static driver')
    306346      END IF
    307347
     
    335375      END SELECT
    336376
     377      SELECT CASE( TRIM(cfg % averaging_mode) )
     378      CASE( 'level', 'height')
     379      CASE DEFAULT
     380         message = "Averaging mode '" // TRIM(cfg % averaging_mode) //&
     381                   "' is not supported. " //&
     382                   "Please select either 'height' or 'level', " //&
     383                   "or omit the --averaging-mode option entirely, which corresponds "//&
     384                   "to the latter."
     385         CALL abort( 'validate_config', message )
     386      END SELECT
     387
     388      IF ( cfg % ug_is_set .NEQV. cfg % vg_is_set )  THEN
     389         message = "You specified only one component of the geostrophic " // &
     390                   "wind. Please specify either both or none."
     391         CALL abort( 'validate_config', message )
     392      END IF
    337393
    338394   END SUBROUTINE validate_config
    339395
    340396
    341    LOGICAL FUNCTION file_present(filename)
     397   LOGICAL FUNCTION file_present(filename, kind)
    342398      CHARACTER(LEN=PATH), INTENT(IN) ::  filename
    343 
    344       INQUIRE(FILE=filename, EXIST=file_present)
    345 
    346       IF (.NOT. file_present)  THEN
    347          message = "The given file '" // "' does not exist."
     399      CHARACTER(LEN=*), INTENT(IN)    ::  kind
     400
     401      IF (LEN(TRIM(filename))==0)  THEN
     402
     403         file_present = .FALSE.
     404         message = "No name was given for the " // TRIM(kind) // " file."
    348405         CALL report('file_present', message)
     406
     407      ELSE
     408
     409         INQUIRE(FILE=filename, EXIST=file_present)
     410
     411         IF (.NOT. file_present)  THEN
     412            message = "The given " // TRIM(kind) // " file '" //               &
     413                      TRIM(filename) // "' does not exist."
     414            CALL report('file_present', message)
     415         END IF
     416
    349417      END IF
    350418
     
    388456#endif
    389457
    390 !
    391 !------------------------------------------------------------------------------
    392 !- Section 1: Write global NetCDF attributes
     458!------------------------------------------------------------------------------
     459!- Section 1: Define NetCDF dimensions and coordinates
     460!------------------------------------------------------------------------------
     461       nt = SIZE(output_file % time)
     462       nx = palm_grid % nx
     463       ny = palm_grid % ny
     464       nz = palm_grid % nz
     465       z0 = palm_grid % z0
     466
     467
     468!
     469!------------------------------------------------------------------------------
     470!- Section 2: Write global NetCDF attributes
    393471!------------------------------------------------------------------------------
    394472       CALL date_and_time(DATE=date_string, TIME=time_string, ZONE=zone_string)
     
    406484       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lat',     TRIM(real_to_str(origin_lat*TO_DEGREES, '(F18.13)'))))
    407485       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_lon',     TRIM(real_to_str(origin_lon*TO_DEGREES, '(F18.13)'))))
     486       ! FIXME: This is the elevation relative to COSMO-DE/D2 sea level and does
     487       ! FIXME: not necessarily comply with DHHN2016 (c.f. PALM Input Data
     488       ! FIXME: Standard v1.9., origin_z)
     489       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'origin_z',       TRIM(real_to_str(z0, '(F18.13)'))))
    408490       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'inifor_version', TRIM(VERSION)))
    409491       CALL check(nf90_put_att(ncid, NF90_GLOBAL, 'palm_version',   '--'))
    410492
    411493!
    412 !------------------------------------------------------------------------------
    413 !- Section 2: Define NetCDF dimensions and coordinates
    414 !------------------------------------------------------------------------------
    415        nt = SIZE(output_file % time)
    416        nx = palm_grid % nx
    417        ny = palm_grid % ny
    418        nz = palm_grid % nz
    419        z0 = palm_grid % z0
    420 
    421494!
    422495!------------------------------------------------------------------------------
     
    618691       ELSE
    619692
    620           nbuffers = SIZE( group % in_var_list )
     693          ! Allocate one input buffer per input_variable. If more quantities
     694          ! have to be computed than input variables exist in this group,
     695          ! allocate more buffers. For instance, in the thermodynamics group,
     696          ! there are three input variabels (PP, T, Qv) and four quantities
     697          ! necessart (P, Theta, Rho, qv) for the corresponding output fields
     698          ! (p0, Theta, qv, ug, and vg)
     699          nbuffers = MAX( group % n_inputs, group % n_output_quantities )
    621700          ALLOCATE( buffer(nbuffers) )
    622701 CALL run_control('time', 'alloc')
    623702         
    624           DO ivar = 1, nbuffers
     703          ! Read in all input variables, leave extra buffers-if any-untouched.
     704          DO ivar = 1, group % n_inputs
    625705
    626706             input_var => group % in_var_list(ivar)
     
    628708             ! Check wheather P or PP is present in input file
    629709             IF (input_var % name == 'P')  THEN
    630                 input_var % name = TRIM( get_pressure_var(input_file) )
     710                input_var % name = TRIM( get_pressure_varname(input_file) )
    631711 CALL run_control('time', 'read')
    632712             END IF
     
    668748!> perturbation, 'PP', and returns the appropriate string.
    669749!------------------------------------------------------------------------------!
    670     CHARACTER(LEN=2) FUNCTION get_pressure_var(input_file) RESULT(var)
     750    CHARACTER(LEN=2) FUNCTION get_pressure_varname(input_file) RESULT(var)
    671751       CHARACTER(LEN=*) ::  input_file
    672752       INTEGER          ::  ncid, varid
     
    692772       CALL check(nf90_close(ncid))
    693773
    694     END FUNCTION get_pressure_var
     774    END FUNCTION get_pressure_varname
    695775
    696776
     
    805885                                   start=start(1:ndim+1) ) )
    806886
    807        CASE ( 'constant scalar profile' )
     887       CASE ( 'constant scalar profile', 'geostrophic', 'internal profile' )
    808888
    809889          CALL check(nf90_put_var( ncid, var%varid, array(1,1,:),              &
  • palm/trunk/UTIL/inifor/src/transform.f90

    r3183 r3395  
    2626! -----------------
    2727! $Id$
     28! Switched addressing of averaging regions from index bounds to list of columns
     29! Added routines for the computation of geostrophic winds including:
     30!    - the downward extrapolation of density (linear) and
     31!    - pressure (hydrostatic equation) and
     32!    - rotation of geostrophic wind components to PALM frame of reference
     33!
     34! 3183 2018-07-27 14:25:55Z suehring
    2835! Introduced new PALM grid stretching
    2936! Removed unnecessary subroutine parameters
     
    5158    USE control
    5259    USE defs,                                                                  &
    53         ONLY: TO_DEGREES, TO_RADIANS, PI, dp
     60        ONLY: G, TO_DEGREES, TO_RADIANS, PI, dp
    5461    USE types
    5562    USE util,                                                                  &       
     
    176183
    177184
    178     SUBROUTINE average_2d(in_arr, out_arr, imin, imax, jmin, jmax)
    179        REAL(dp), INTENT(IN)  ::  in_arr(0:,0:,0:)
    180        REAL(dp), INTENT(OUT) ::  out_arr(0:)
    181        INTEGER, INTENT(IN)   ::  imin, imax, jmin, jmax
    182 
    183        INTEGER  ::  i, j, k
     185    SUBROUTINE average_2d(in_arr, out_arr, ii, jj)
     186       REAL(dp), INTENT(IN)              ::  in_arr(0:,0:,0:)
     187       REAL(dp), INTENT(OUT)             ::  out_arr(0:)
     188       INTEGER, INTENT(IN), DIMENSION(:) ::  ii, jj
     189
     190       INTEGER  ::  i, j, k, l
    184191       REAL(dp) ::  ni
    185        
    186        IF (imin < 0)  CALL abort('average_2d', "imin < 0.")
    187        IF (jmin < 0)  CALL abort('average_2d', "jmin < 0.")
    188        IF (imax > UBOUND(in_arr, 1))  CALL abort('average_2d', "imax out of i bound.")
    189        IF (jmax > UBOUND(in_arr, 2))  CALL abort('average_2d', "jmax out of j bound.")
     192
     193       IF (SIZE(ii) .NE. SIZE(jj))  THEN
     194          message = "Length of 'ii' and 'jj' index lists do not match." //     &
     195             NEW_LINE(' ') // "ii has " // str(SIZE(ii)) // " elements, " //        &
     196             NEW_LINE(' ') // "jj has " // str(SIZE(jj)) // "."
     197          CALL abort('average_2d', message)
     198       END IF
    190199
    191200       DO k = 0, UBOUND(out_arr, 1)
    192201
    193202          out_arr(k) = 0.0_dp
    194           DO j = jmin, jmax
    195           DO i = imin, imax
    196              out_arr(k) = out_arr(k) + in_arr(i, j, k)
     203          DO l = 1, UBOUND(ii, 1)
     204             i = ii(l)
     205             j = jj(l)
     206             out_arr(k) = out_arr(k) +&
     207                         in_arr(i, j, k)
    197208          END DO
    198           END DO
    199 
    200        END DO
    201    
    202        ! devide by number of grid points
    203        ni = 1.0_dp / ( (imax - imin + 1) * (jmax - jmin + 1) )
     209
     210       END DO
     211
     212       ni = 1.0_dp / SIZE(ii)
    204213       out_arr(:) = out_arr(:) * ni
    205214
     
    212221       REAL(dp), DIMENSION(:,:,:), INTENT(OUT) ::  palm_array
    213222       REAL(dp), DIMENSION(:,:,:), ALLOCATABLE ::  intermediate_array
    214        INTEGER ::  nx, ny, nz
     223       INTEGER ::  nx, ny, nlev
    215224
    216225       nx = palm_intermediate % nx
    217226       ny = palm_intermediate % ny
    218        nz = palm_intermediate % nz ! nlev
     227       nlev = palm_intermediate % nz ! nlev
    219228
    220229       ! Interpolate from COSMO-DE to intermediate grid. Allocating with one
    221230       ! less point in the vertical, since scalars like T have 50 instead of 51
    222231       ! points in COSMO-DE.
    223        ALLOCATE(intermediate_array(0:nx, 0:ny, 0:nz-1)) !
     232       ALLOCATE(intermediate_array(0:nx, 0:ny, 0:nlev-1)) !
    224233
    225234       CALL interpolate_2d(source_array, intermediate_array, palm_intermediate)
     
    234243
    235244
    236     SUBROUTINE average_profile(source_array, profile_array, imin, imax, jmin, jmax,&
    237                                palm_intermediate, palm_grid)
    238        TYPE(grid_definition), INTENT(IN)          ::  palm_intermediate, palm_grid
     245    SUBROUTINE average_profile(source_array, profile_array, avg_grid)
     246       TYPE(grid_definition), INTENT(IN)          ::  avg_grid
    239247       REAL(dp), DIMENSION(:,:,:), INTENT(IN)     ::  source_array
    240        INTEGER, INTENT(IN)                        ::  imin, imax, jmin, jmax
    241        REAL(dp), DIMENSION(0:,0:,0:), INTENT(OUT) ::  profile_array
    242        REAL(dp), DIMENSION(:,:,:), ALLOCATABLE    ::  intermediate_array
    243        INTEGER                                    ::  nx, ny, nz
    244 
    245        nx = palm_intermediate % nx
    246        ny = palm_intermediate % ny
    247        nz = palm_intermediate % nz
    248        ALLOCATE(intermediate_array(0:nx, 0:ny, 0:nz-1))
    249        intermediate_array(:,:,:) = 0.0_dp
    250 
    251        ! average input array to intermediate profile
    252        CALL average_2d(source_array, intermediate_array(0,0,:), imin, imax, jmin, jmax)
    253 
    254        ! vertically interpolate to ouput array
    255        CALL interpolate_1d(intermediate_array, profile_array, palm_grid)
    256 
    257        DEALLOCATE(intermediate_array)
     248       REAL(dp), DIMENSION(:), INTENT(OUT)        ::  profile_array
     249
     250       INTEGER ::  i_source, j_source, k_profile, k_source, l, m
     251
     252       REAL ::  ni_columns
     253
     254       profile_array(:) = 0.0_dp
     255
     256       DO l = 1, avg_grid % n_columns
     257          i_source = avg_grid % iii(l)
     258          j_source = avg_grid % jjj(l)
     259
     260          DO k_profile = avg_grid % k_min, UBOUND(profile_array, 1) ! PALM levels
     261
     262             DO m = 1, 2 ! vertical interpolation neighbours
     263
     264                k_source = avg_grid % kkk(l, k_profile, m)
     265
     266                profile_array(k_profile) = profile_array(k_profile)    &
     267                   + avg_grid % w(l, k_profile, m)                             &
     268                   * source_array(i_source, j_source, k_source)
     269
     270             END DO ! m, vertical interpolation neighbours
     271
     272          END DO    ! k_profile, PALM levels
     273
     274       END DO       ! l, horizontal neighbours
     275
     276       ni_columns = 1.0_dp / avg_grid % n_columns
     277       profile_array(:) = profile_array(:) * ni_columns
     278
     279       ! Extrapolate constant to the bottom
     280       profile_array(1:avg_grid % k_min-1) = profile_array(avg_grid % k_min)
    258281
    259282    END SUBROUTINE average_profile
    260283
     284
     285    SUBROUTINE extrapolate_density(rho, avg_grid)
     286       REAL(dp), DIMENSION(:), INTENT(INOUT) ::  rho
     287       TYPE(grid_definition), INTENT(IN)     ::  avg_grid
     288
     289       REAL(dp) ::  drhodz, dz, zk, rhok
     290       INTEGER  ::  k_min
     291
     292       k_min  = avg_grid % k_min
     293       zk     = avg_grid % z(k_min)
     294       rhok   = rho(k_min)
     295       dz     = avg_grid % z(k_min + 1) - avg_grid % z(k_min)
     296       drhodz = (rho(k_min + 1) - rho(k_min)) / dz
     297
     298       rho(1:k_min-1) = rhok + drhodz * (avg_grid % z(1:k_min-1) - zk)
     299
     300    END SUBROUTINE extrapolate_density
     301
     302
     303    SUBROUTINE extrapolate_pressure(p, rho, avg_grid)
     304       REAL(dp), DIMENSION(:), INTENT(IN)    ::  rho
     305       REAL(dp), DIMENSION(:), INTENT(INOUT) ::  p
     306       TYPE(grid_definition), INTENT(IN)     ::  avg_grid
     307
     308       REAL(dp) ::  drhodz, dz, zk, rhok
     309       INTEGER  ::  k, k_min
     310
     311       k_min = avg_grid % k_min
     312       zk    = avg_grid % z(k_min)
     313       rhok  = rho(k_min)
     314       dz    = avg_grid % z(k_min + 1) - avg_grid % z(k_min)
     315       drhodz = 0.5_dp * (rho(k_min + 1) - rho(k_min)) / dz
     316
     317       DO k = 1, k_min-1
     318          p(k) = constant_density_pressure(p(k_min), zk, rhok, drhodz,         &
     319                                           avg_grid % z(k), G)
     320       END DO
     321
     322    END SUBROUTINE extrapolate_pressure
     323
     324
     325!------------------------------------------------------------------------------!
     326! Description:
     327! ------------
     328!> Takes the averaged pressure profile <p> and sets the lowest entry to the
     329!> extrapolated pressure at the surface.
     330!------------------------------------------------------------------------------!
     331    SUBROUTINE get_surface_pressure(p, rho, avg_grid)
     332       REAL(dp), DIMENSION(:), INTENT(IN)    ::  rho
     333       REAL(dp), DIMENSION(:), INTENT(INOUT) ::  p
     334       TYPE(grid_definition), INTENT(IN)     ::  avg_grid
     335
     336       REAL(dp) ::  drhodz, dz, zk, rhok
     337       INTEGER  ::  k, k_min
     338
     339       k_min = avg_grid % k_min
     340       zk    = avg_grid % z(k_min)
     341       rhok  = rho(k_min)
     342       dz    = avg_grid % z(k_min + 1) - avg_grid % z(k_min)
     343       drhodz = 0.5_dp * (rho(k_min + 1) - rho(k_min)) / dz
     344
     345       p(1) = constant_density_pressure(p(k_min), zk, rhok, drhodz,            &
     346                                        0.0, G)
     347
     348    END SUBROUTINE get_surface_pressure
     349
     350
     351    FUNCTION constant_density_pressure(pk, zk, rhok, drhodz, z, g)  RESULT(p)
     352
     353       REAL(dp), INTENT(IN)  ::  pk, zk, rhok, drhodz, g, z
     354       REAL(dp) ::  p
     355
     356       p = pk + ( zk - z ) * g * ( rhok + 0.5*drhodz * (zk - z) )
     357
     358    END FUNCTION constant_density_pressure
     359
     360!-----------------------------------------------------------------------------!
     361! Description:
     362! -----------
     363!> This routine computes profiles of the two geostrophic wind components ug and
     364!> vg.
     365!-----------------------------------------------------------------------------!
     366    SUBROUTINE geostrophic_winds(p_north, p_south, p_east, p_west, rho, f3,    &
     367                                 Lx, Ly, phi_n, lam_n, phi_g, lam_g, ug, vg)
     368
     369    REAL(dp), DIMENSION(:), INTENT(IN)  ::  p_north, p_south, p_east, p_west,  &
     370                                            rho
     371    REAL(dp), INTENT(IN)                ::  f3, Lx, Ly, phi_n, lam_n, phi_g, lam_g
     372    REAL(dp), DIMENSION(:), INTENT(OUT) ::  ug, vg
     373    REAL(dp)                            ::  facx, facy
     374
     375    facx = 1.0_dp / (Lx * f3)
     376    facy = 1.0_dp / (Ly * f3)
     377    ug(:) = - facy / rho(:) * (p_north(:) - p_south(:))
     378    vg(:) =   facx / rho(:) * (p_east(:) - p_west(:))
     379
     380    CALL rotate_vector_field(                                                  &
     381       ug, vg, angle = meridian_convergence_rotated(phi_n, lam_n, phi_g, lam_g)&
     382    )
     383
     384    END SUBROUTINE geostrophic_winds
    261385
    262386
     
    419543       
    420544
     545    SUBROUTINE rotate_vector_field(x, y, angle)
     546       REAL(dp), DIMENSION(:), INTENT(INOUT) :: x, y  !< x and y coodrinate in arbitrary units
     547       REAL(dp), INTENT(IN)                  :: angle !< rotation angle [deg]
     548
     549       INTEGER  :: i
     550       REAL(dp) :: sine, cosine, v_rot(2), rotation(2,2)
     551
     552       sine = SIN(angle * TO_RADIANS)
     553       cosine = COS(angle * TO_RADIANS)
     554       ! RESAHPE() fills columns first, so the rotation matrix becomes
     555       !
     556       ! rotation = [ cosine   -sine  ]
     557       !            [  sine    cosine ]
     558       rotation = RESHAPE( (/cosine, sine, -sine, cosine/), (/2, 2/) )
     559
     560       DO i = LBOUND(x, 1), UBOUND(x, 1)
     561
     562          v_rot(:) = MATMUL(rotation, (/x(i), y(i)/))
     563
     564          x(i) = v_rot(1)
     565          y(i) = v_rot(2)
     566
     567       END DO
     568
     569    END SUBROUTINE rotate_vector_field
     570
     571
     572
     573!------------------------------------------------------------------------------!
     574! Description:
     575! ------------
     576!> This routine computes the local meridian convergence between a rotated-pole
     577!> and a geographical system using the Eq. (6) given in the DWD manual
     578!>
     579!>    Baldauf et al. (2018), Beschreibung des operationelle Kürzestfrist-
     580!>    vorhersagemodells COSMO-D2 und COSMO-D2-EPS und seiner Ausgabe in die
     581!>    Datenbanken des DWD.
     582!>    https://www.dwd.de/SharedDocs/downloads/DE/modelldokumentationen/nwv/cosmo_d2/cosmo_d2_dbbeschr_aktuell.pdf?__blob=publicationFile&v=2
     583!>
     584    FUNCTION meridian_convergence_rotated(phi_n, lam_n, phi_g, lam_g)          &
     585       RESULT(delta)
     586
     587       REAL(dp), INTENT(IN) ::  phi_n, lam_n, phi_g, lam_g
     588       REAL(dp)             ::  delta
     589
     590       delta = atan2( COS(phi_n) * SIN(lam_n - lam_g),                         &
     591                      COS(phi_g) * SIN(phi_n) -                                &
     592                      SIN(phi_g) * COS(phi_n) * COS(lam_n - lam_g) )
     593
     594    END FUNCTION meridian_convergence_rotated
    421595
    422596!------------------------------------------------------------------------------!
     
    512686
    513687   
    514     SUBROUTINE find_vertical_neighbours_and_weights(palm_grid, palm_intermediate)
     688    SUBROUTINE find_vertical_neighbours_and_weights_interp( palm_grid,         &
     689                                                            palm_intermediate )
    515690       TYPE(grid_definition), INTENT(INOUT) ::  palm_grid
    516691       TYPE(grid_definition), INTENT(IN)    ::  palm_intermediate
     
    528703
    529704       ! in each column of the fine grid, find vertical neighbours of every cell
     705       DO j = 0, ny
    530706       DO i = 0, nx
    531        DO j = 0, ny
    532707
    533708          k_intermediate = 0
     
    536711          column_top  = palm_intermediate % h(i,j,nlev)
    537712
    538           ! scan through palm_grid column until and set neighbour indices in
     713          ! scan through palm_grid column and set neighbour indices in
    539714          ! case current_height is either below column_base, in the current
    540715          ! cell, or above column_top. Keep increasing current cell index until
     
    580755                   h_top    = palm_intermediate % h(i,j,k_intermediate+1)
    581756                   h_bottom = palm_intermediate % h(i,j,k_intermediate)
     757                   point_is_in_current_cell = (                                &
     758                      current_height >= h_bottom .AND.                         &
     759                      current_height <  h_top                                  &
     760                   )
     761                END DO
     762
     763                IF (k_intermediate > nlev-1)  THEN
     764                   message = "Index " // TRIM(str(k_intermediate)) //          &
     765                             " is above intermediate grid range."
     766                   CALL abort('find_vertical_neighbours', message)
     767                END IF
     768   
     769                palm_grid % kk(i,j,k,1) = k_intermediate
     770                palm_grid % kk(i,j,k,2) = k_intermediate + 1
     771
     772                ! copmute vertical weights
     773                weight = (h_top - current_height) / (h_top - h_bottom)
     774                palm_grid % w_verti(i,j,k,1) = weight
     775                palm_grid % w_verti(i,j,k,2) = 1.0_dp - weight
     776             END IF
     777
     778          END DO
     779
     780       END DO
     781       END DO
     782
     783    END SUBROUTINE find_vertical_neighbours_and_weights_interp
     784
     785
     786    SUBROUTINE find_vertical_neighbours_and_weights_average( avg_grid )
     787       TYPE(grid_definition), INTENT(INOUT) ::  avg_grid
     788
     789       INTEGER  ::  i, j, k_palm, k_intermediate, l, nlev
     790       LOGICAL  ::  point_is_below_grid, point_is_above_grid,                  &
     791                    point_is_in_current_cell
     792       REAL(dp) ::  current_height, column_base, column_top, h_top, h_bottom,  &
     793                    weight
     794
     795
     796       avg_grid % k_min = LBOUND(avg_grid % z, 1)
     797
     798       nlev = SIZE(avg_grid % cosmo_h, 3)
     799
     800       ! in each column of the fine grid, find vertical neighbours of every cell
     801       DO l = 1, avg_grid % n_columns
     802
     803          i = avg_grid % iii(l)
     804          j = avg_grid % jjj(l)
     805
     806          column_base = avg_grid % cosmo_h(i,j,1)
     807          column_top  = avg_grid % cosmo_h(i,j,nlev)
     808
     809          ! scan through avg_grid column until and set neighbour indices in
     810          ! case current_height is either below column_base, in the current
     811          ! cell, or above column_top. Keep increasing current cell index until
     812          ! the current cell overlaps with the current_height.
     813          k_intermediate = 1 !avg_grid % cosmo_h is indezed 1-based.
     814          DO k_palm = 1, avg_grid % nz
     815
     816             ! Memorize the top and bottom boundaries of the coarse cell and the
     817             ! current height within it
     818             current_height = avg_grid % z(k_palm) + avg_grid % z0
     819             h_top    = avg_grid % cosmo_h(i,j,k_intermediate+1)
     820             h_bottom = avg_grid % cosmo_h(i,j,k_intermediate)
     821
     822             point_is_above_grid = (current_height > column_top) !22000m, very unlikely
     823             point_is_below_grid = (current_height < column_base)
     824
     825             point_is_in_current_cell = (                                      &
     826                current_height >= h_bottom .AND.                               &
     827                current_height <  h_top                                        &
     828             )
     829
     830             ! set default weights
     831             avg_grid % w(l,k_palm,1:2) = 0.0_dp
     832
     833             IF (point_is_above_grid)  THEN
     834
     835                avg_grid % kkk(l,k_palm,1:2) = nlev
     836                avg_grid % w(l,k_palm,1:2) = - 2.0_dp
     837
     838                message = "PALM-4U grid extends above COSMO-DE model top."
     839                CALL abort('find_vertical_neighbours_and_weights_average', message)
     840
     841             ELSE IF (point_is_below_grid)  THEN
     842
     843                avg_grid % kkk(l,k_palm,1:2) = 0
     844                avg_grid % w(l,k_palm,1:2) = - 2.0_dp
     845                avg_grid % k_min = MAX(k_palm + 1, avg_grid % k_min)
     846             ELSE
     847                ! cycle through intermediate levels until current
     848                ! intermediate-grid cell overlaps with current_height
     849                DO WHILE (.NOT. point_is_in_current_cell .AND. k_intermediate <= nlev-1)
     850                   k_intermediate = k_intermediate + 1
     851
     852                   h_top    = avg_grid % cosmo_h(i,j,k_intermediate+1)
     853                   h_bottom = avg_grid % cosmo_h(i,j,k_intermediate)
    582854                   point_is_in_current_cell = (                                &
    583855                      current_height >= h_bottom .AND.                         &
     
    594866                END IF
    595867   
    596                 palm_grid % kk(i,j,k,1) = k_intermediate
    597                 palm_grid % kk(i,j,k,2) = k_intermediate + 1
     868                avg_grid % kkk(l,k_palm,1) = k_intermediate
     869                avg_grid % kkk(l,k_palm,2) = k_intermediate + 1
    598870
    599871                ! copmute vertical weights
    600872                weight = (h_top - current_height) / (h_top - h_bottom)
    601                 palm_grid % w_verti(i,j,k,1) = weight
    602                 palm_grid % w_verti(i,j,k,2) = 1.0_dp - weight
     873                avg_grid % w(l,k_palm,1) = weight
     874                avg_grid % w(l,k_palm,2) = 1.0_dp - weight
    603875             END IF
    604876
    605           END DO
    606 
    607        END DO
    608        END DO
    609 
    610     END SUBROUTINE find_vertical_neighbours_and_weights
     877          END DO ! k, PALM levels
     878       END DO ! l, averaging columns
     879 
     880    END SUBROUTINE find_vertical_neighbours_and_weights_average
    611881
    612882!------------------------------------------------------------------------------!
  • palm/trunk/UTIL/inifor/src/types.f90

    r3183 r3395  
    2626! -----------------
    2727! $Id$
     28! Added *_is_set LOGICALs to inifor_config type to indicate option invocation
     29!     from the command-line
     30! Added 1D index vertical weights lists to support addressing averaging regions
     31!     by list of columns instead of index bounds
     32!
     33!
     34! 3183 2018-07-27 14:25:55Z suehring
    2835! Introduced new PALM grid stretching:
    2936! - Converted vertical grid_definition coordinte variables to pointers
     
    6673
    6774    CHARACTER(LEN=SNAME) ::  flow_prefix          !< Prefix of flow input files, e.g. 'laf' for COSMO-DE analyses
     75    CHARACTER(LEN=SNAME) ::  input_prefix         !< Prefix of all input files, e.g. 'laf' for COSMO-DE analyses
     76    CHARACTER(LEN=SNAME) ::  radiation_prefix     !< Prefix of radiation input files, e.g 'laf' for COSMO-DE analyses
    6877    CHARACTER(LEN=SNAME) ::  soil_prefix          !< Prefix of soil input files, e.g. 'laf' for COSMO-DE analyses
    69     CHARACTER(LEN=SNAME) ::  radiation_prefix     !< Prefix of radiation input files, e.g 'laf' for COSMO-DE analyses
    7078    CHARACTER(LEN=SNAME) ::  soilmoisture_prefix  !< Prefix of input files for soil moisture spin-up, e.g 'laf' for COSMO-DE analyses
    7179
     80    CHARACTER(LEN=SNAME) ::  averaging_mode
    7281    CHARACTER(LEN=SNAME) ::  bc_mode
    7382    CHARACTER(LEN=SNAME) ::  ic_mode
     
    7786    REAL(dp)             ::  ug
    7887    REAL(dp)             ::  vg
    79     REAL(dp)             ::  z0 !< Elevation of the PALM-4U domain above sea level [m]
     88    REAL(dp)             ::  z0                   !< Elevation of the PALM-4U domain above sea level [m]
     89    REAL(dp)             ::  averaging_angle      !< latitudal and longitudal width of averaging regions [deg]
     90   
     91
     92    LOGICAL              ::  debug
     93    LOGICAL              ::  p0_is_set
     94    LOGICAL              ::  ug_is_set
     95    LOGICAL              ::  vg_is_set
     96    LOGICAL              ::  flow_prefix_is_set          !<
     97    LOGICAL              ::  input_prefix_is_set         !<
     98    LOGICAL              ::  radiation_prefix_is_set     !<
     99    LOGICAL              ::  soil_prefix_is_set          !<
     100    LOGICAL              ::  soilmoisture_prefix_is_set  !<
    80101 END TYPE inifor_config
    81102
    82103 TYPE grid_definition
    83104    CHARACTER(LEN=SNAME)  ::  name(3)       !< names of the grid dimensions, e.g. (/'x', 'y', 'z'/) or (/'latitude', 'longitude', 'height'/)
     105    CHARACTER(LEN=SNAME)  ::  kind          !< names of the grid dimensions, e.g. (/'x', 'y', 'z'/) or (/'latitude', 'longitude', 'height'/)
     106    INTEGER               ::  k_min         !< Index of lowest PALM grid level that is not cut by local COSMO orography; vertically separates interpolation and extrapolation region.
    84107    INTEGER               ::  nx            !< number of gridpoints in the first dimension
    85108    INTEGER               ::  ny            !< number of gridpoints in the second dimension
    86     INTEGER               ::  nz            !< number of gridpoints in the third dimension
     109    INTEGER               ::  nz            !< number of gridpoints in the third dimension, used for PALM points
     110    INTEGER               ::  nlev          !< number of COSMO grid levels
     111    INTEGER               ::  n_columns     !< number of averaging columns of the source grid
    87112    INTEGER, ALLOCATABLE  ::  ii(:,:,:)     !< Given a point (i,j,k) in the PALM-4U grid, ii(i,j,l) gives the x index of the l'th horizontl neighbour on the COSMO-DE grid.
    88113    INTEGER, ALLOCATABLE  ::  jj(:,:,:)     !< Given a point (i,j,k) in the PALM-4U grid, jj(i,j,l) gives the y index of the l'th horizontl neighbour on the COSMO-DE grid.
    89114    INTEGER, ALLOCATABLE  ::  kk(:,:,:,:)   !< Given a point (i,j,k) in the PALM-4U grid, kk(i,j,k,l) gives the z index of the l'th vertical neighbour in the intermediate grid.
     115    INTEGER, ALLOCATABLE  ::  iii(:)        !< profile averaging neighbour indices
     116    INTEGER, ALLOCATABLE  ::  jjj(:)        !< profile averaging neighbour indices
     117    INTEGER, ALLOCATABLE  ::  kkk(:,:,:)    !< indices of vertical interpolation neightbours, kkk(<source column>, <PALM k level>, <neighbour index>)
    90118    REAL(dp)              ::  lx            !< domain length in the first dimension [m]
    91119    REAL(dp)              ::  ly            !< domain length in the second dimension [m]
     
    97125    REAL(dp), POINTER     ::  z(:)          !< coordinates of cell centers in z direction [m]
    98126    REAL(dp), ALLOCATABLE ::  h(:,:,:)      !< heights grid point for intermediate grids [m]
     127    REAL(dp), POINTER     ::  cosmo_h(:,:,:)!< pointer to appropriate COSMO level heights (scalar/w) [m]
    99128    REAL(dp), POINTER     ::  hhl(:,:,:)    !< heights of half layers (cell faces) above sea level in COSMO-DE, read in from
    100129    REAL(dp), POINTER     ::  hfl(:,:,:)    !< heights of full layers (cell centres) above sea level in COSMO-DE, computed as arithmetic average of hhl
     
    115144    REAL(dp), ALLOCATABLE ::  w_horiz(:,:,:)   !< weights for bilinear horizontal interpolation
    116145    REAL(dp), ALLOCATABLE ::  w_verti(:,:,:,:) !< weights for linear vertical interpolation
     146    REAL(dp), ALLOCATABLE ::  w(:,:,:)      !< vertical interpolation weights, w(<source_column>, <PALM k level>, <neighbour index>) [-]
    117147 END TYPE grid_definition
    118148
     
    150180    CHARACTER(LEN=SNAME)                  ::  task                      !< Processing task that generates this variable, e.g. 'interpolate_2d' or 'average profile'
    151181    LOGICAL                               ::  to_be_processed = .FALSE. !< INIFOR flag indicating whether variable shall be processed
     182    LOGICAL                               ::  is_internal = .FALSE.     !< INIFOR flag indicating whether variable shall be written to netCDF file (.FALSE.) or kept for later (.TRUE.)
    152183    LOGICAL                               ::  is_read = .FALSE.         !< INIFOR flag indicating whether variable has been read
    153     LOGICAL                               ::  is_upside_down  = .FALSE. !< INIFOR flag indicating whether variable shall be processed
     184    LOGICAL                               ::  is_upside_down  = .FALSE. !< INIFOR flag indicating whether vertical dimension is reversed (typically the case with COSMO-DE atmospheric fields)
    154185    TYPE(grid_definition), POINTER        ::  grid                      !< Pointer to the corresponding output grid
    155186    TYPE(grid_definition), POINTER        ::  intermediate_grid         !< Pointer to the corresponding intermediate grid
     187    TYPE(grid_definition), POINTER        ::  averaging_grid         !< Pointer to the corresponding intermediate grid
    156188 END TYPE nc_var
    157189
     
    159191 TYPE io_group                                          !< Input/Output group, groups together output variabels that share their input variables. For instance, all boundary surfaces and initialization fields of the potential temperature are base on T and p.
    160192    INTEGER                          ::  nt             !< maximum number of output time steps across all output variables
    161     INTEGER                          ::  nv             !< number of output variables
     193    INTEGER                          ::  nv             !< number of netCDF output variables
     194    INTEGER                          ::  n_inputs       !< number of input variables
     195    INTEGER                          ::  n_output_quantities !< number of physical quantities required for computing netCDF output variables
    162196    CHARACTER(LEN=SNAME)             ::  kind           !< kind of I/O group
    163197    CHARACTER(LEN=PATH), ALLOCATABLE ::  in_files(:)    !< list of nt input files
  • palm/trunk/UTIL/inifor/src/util.f90

    r3183 r3395  
    2626! -----------------
    2727! $Id$
     28! New routines for computing potential temperature and moist air density
     29! Increased number of digits in real-to-str conversion
     30!
     31! 3183 2018-07-27 14:25:55Z suehring
    2832! Improved real-to-string conversion
    2933!
     
    4852    USE defs,                                                                  &
    4953        ONLY :  dp, PI, DATE, SNAME
     54    USE types,                                                                 &
     55        ONLY :  grid_definition
    5056
    5157    IMPLICIT NONE
     
    281287
    282288
     289!------------------------------------------------------------------------------!
     290! Description:
     291! ------------
     292!> Converts the absolute temperature to the potential temperature in place using
     293!> the identity a^b = e^(b ln(a)).
     294!>
     295!>     theta = T * (p_ref/p)^(R/c_p) = T * e^( R/c_p * ln(p_ref/p) )
     296!------------------------------------------------------------------------------!
     297    SUBROUTINE potential_temperature(t, p, p_ref, r, cp)
     298       REAL(dp), DIMENSION(:,:,:), INTENT(INOUT) ::  t
     299       REAL(dp), DIMENSION(:,:,:), INTENT(IN)    ::  p
     300       REAL(dp), INTENT(IN)                      ::  p_ref, r, cp
     301       REAL(dp)                                  ::  rcp
     302
     303       rcp = r/cp
     304       t(:,:,:) =  t(:,:,:) * EXP( rcp * LOG(p_ref / p(:,:,:)) )
     305
     306    END SUBROUTINE potential_temperature
     307
     308
     309!------------------------------------------------------------------------------!
     310! Description:
     311! ------------
     312!> Compute the density in place of the given temperature (t_rho).
     313!------------------------------------------------------------------------------!
     314   SUBROUTINE moist_density(t_rho, p, qv, rd, rv)
     315       REAL(dp), DIMENSION(:,:,:), INTENT(INOUT) ::  t_rho
     316       REAL(dp), DIMENSION(:,:,:), INTENT(IN)    ::  p, qv
     317       REAL(dp), INTENT(IN)                      ::  rd, rv
     318
     319       t_rho(:,:,:) = p(:,:,:) / (                                             &
     320          (rv * qv(:,:,:) + rd * (1.0_dp - qv(:,:,:))) * t_rho(:,:,:)          &
     321       )
     322
     323    END SUBROUTINE moist_density
     324
     325
    283326    ! Convert a real number to a string in scientific notation
    284327    ! showing four significant digits.
     
    298341
    299342
    300     CHARACTER(LEN=12) FUNCTION real_to_str_f(val)
     343    CHARACTER(LEN=16) FUNCTION real_to_str_f(val)
    301344
    302345        REAL(dp), INTENT(IN) ::  val
    303346
    304         WRITE(real_to_str_f, '(F12.4)') val
     347        WRITE(real_to_str_f, '(F16.8)') val
    305348        real_to_str_f = ADJUSTL(real_to_str_f)
    306349
Note: See TracChangeset for help on using the changeset viewer.