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

inifor: Added computation of geostrophic winds from COSMO input

Location:
palm/trunk/UTIL/inifor/src
Files:
1 deleted
8 edited

Legend:

Unmodified
Added
Removed
  • 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.