Home

Context Navigation

source: palm/trunk/UTIL/inifor/src/inifor_grid.f90 @ 3698

Last change on this file since 3698 was 3680, checked in by knoop, 5 years ago
Added cpp-option netcdf to inifor
Property svn:keywords set to `Id`
File size: 211.4 KB

Rev	Line
[3447]	1	!> @file src/inifor_grid.f90
[2696]	2	!------------------------------------------------------------------------------!
[2718]	3	! This file is part of the PALM model system.
[2696]	4	!
[2718]	5	! PALM is free software: you can redistribute it and/or modify it under the
	6	! terms of the GNU General Public License as published by the Free Software
	7	! Foundation, either version 3 of the License, or (at your option) any later
[2696]	8	! version.
	9	!
[2718]	10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
	11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
	12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
[2696]	13	!
	14	! You should have received a copy of the GNU General Public License along with
	15	! PALM. If not, see <http://www.gnu.org/licenses/>.
	16	!
[2718]	17	! Copyright 2017-2018 Leibniz Universitaet Hannover
	18	! Copyright 2017-2018 Deutscher Wetterdienst Offenbach
[2696]	19	!------------------------------------------------------------------------------!
	20	!
	21	! Current revisions:
	22	! -----------------
[3183]	23	!
	24	!
	25	! Former revisions:
	26	! -----------------
	27	! $Id: inifor_grid.f90 3680 2019-01-18 14:54:12Z suehring $
[3678]	28	! bugfix: Avoid empty averaging regions for small PALM domains
	29	! change: moved get_input_file_list() to io module
	30	!
	31	!
	32	! 3618 2018-12-10 13:25:22Z eckhard
[3618]	33	! Prefixed all INIFOR modules with inifor_, removed unused variables
	34	!
	35	!
	36	! 3615 2018-12-10 07:21:03Z raasch
[3615]	37	! bugfix: abort replaced by inifor_abort
	38	!
	39	! 3614 2018-12-10 07:05:46Z raasch
[3614]	40	! unused variables removed
	41	!
	42	! 3613 2018-12-07 18:20:37Z eckhard
[3613]	43	! Average initial profiles only over PALM domain region, not the
	44	! geostrophic-wind averaging region
	45	! Fix unintended modification of COSMO heights
	46	! Fix commenting out of setup_volumentric setting
	47	! Moved version output to setup_parameters()
	48	!
	49	!
	50	! 3557 2018-11-22 16:01:22Z eckhard
[3557]	51	! Updated documentation
	52	!
	53	!
	54	! 3537 2018-11-20 10:53:14Z eckhard
[3537]	55	! Read COSMO domain extents and soil depths from input files
	56	! Report averaging mode and debugging mode in log
	57	!
	58	!
	59	! 3456 2018-10-30 14:29:54Z eckhard
[3456]	60	! Remove surface forcing from netCDF output (radiation + precipitation)
	61	! NetCDf output of internal arrays only with --debug option
	62	!
	63	!
	64	! 3447 2018-10-29 15:52:54Z eckhard
[3447]	65	! Renamed source files for compatibilty with PALM build system
	66	!
	67	!
	68	! 3395 2018-10-22 17:32:49Z eckhard
[3395]	69	! Added computation of geostrophic winds form COSMO pressure fields
	70	! Introduced averaging grids and internal 'output' variables for computation of
	71	! geostrophic and large-scale forcings
	72	! Merged qv, uv and vg output variables and 'temperature' IO group into new
	73	! 'thermodynamics' IP group
	74	! Cleaned up terminal output, show some messages only with --debug option
	75	! Improved variable naming and readability
	76	!
	77	!
	78	! 3183 2018-07-27 14:25:55Z suehring
[3182]	79	! Introduced new PALM grid stretching
	80	! Updated variable names and metadata for PIDS v1.9 compatibility
	81	! Better compatibility with older Intel compilers:
	82	! - avoiding implicit array allocation with new get_netcdf_variable()
	83	! subroutine instead of function
	84	! Improved command line interface:
	85	! - Produce forcing variables when '-mode profile' is being used
	86	! - Renamend initial-condition mode variable 'mode' to 'ic_mode'
	87	! - Improved handling of the start date string
	88	! Removed unnecessary variables and routines
	89	!
[2696]	90	!
[3183]	91	! 3182 2018-07-27 13:36:03Z suehring
[2696]	92	! Initial revision
	93	!
	94	!
	95	!
	96	! Authors:
	97	! --------
[3557]	98	!> @author Eckhard Kadasch (Deutscher Wetterdienst, Offenbach)
[2696]	99	!
	100	!------------------------------------------------------------------------------!
	101	! Description:
	102	! ------------
	103	!> The grid module contains all variables and routines to specify and work with
	104	!> the numerical grids in INIFOR. By convention, all angles are stored in
	105	!> radians.
	106	!------------------------------------------------------------------------------!
	107
[3680]	108	#if defined ( __netcdf )
[3618]	109	MODULE inifor_grid
[2696]	110
[3618]	111	USE inifor_control
	112	USE inifor_defs, &
[2696]	113	ONLY: DATE, EARTH_RADIUS, TO_RADIANS, TO_DEGREES, PI, dp, hp, sp, &
[3182]	114	SNAME, LNAME, PATH, FORCING_STEP, WATER_ID, FILL_ITERATIONS, &
[3395]	115	BETA, P_SL, T_SL, BETA, RD, RV, G, P_REF, RD_PALM, CP_PALM, &
	116	RHO_L, OMEGA, HECTO
[3618]	117	USE inifor_io, &
[3537]	118	ONLY: get_netcdf_attribute, get_netcdf_dim_vector, &
	119	get_netcdf_variable, parse_command_line_arguments, &
[3678]	120	get_input_file_list, validate_config
[3618]	121	USE inifor_transform, &
[3395]	122	ONLY: average_2d, rotate_to_cosmo, find_horizontal_neighbours,&
[2696]	123	compute_horizontal_interp_weights, &
[3395]	124	find_vertical_neighbours_and_weights_interp, &
	125	find_vertical_neighbours_and_weights_average, interpolate_2d, &
	126	gamma_from_hemisphere, phic_to_phin, lamc_to_lamn, &
	127	project, centre_velocities, phi2phirot, rla2rlarot, uv2uvrot, &
	128	phirot2phi, rlarot2rla
[3618]	129	USE inifor_types
	130	USE inifor_util
	131	USE netcdf, &
	132	ONLY: NF90_MAX_NAME, NF90_MAX_VAR_DIMS
[2696]	133
	134	IMPLICIT NONE
	135
	136	SAVE
	137
[3395]	138	REAL(dp) :: averaging_angle = 0.0_dp !< latitudal and longitudal width of averaging regions [rad]
	139	REAL(dp) :: averaging_width_ns = 0.0_dp !< longitudal width of averaging regions [m]
	140	REAL(dp) :: averaging_width_ew = 0.0_dp !< latitudal width of averaging regions [m]
[3182]	141	REAL(dp) :: phi_equat = 0.0_dp !< latitude of rotated equator of COSMO-DE grid [rad]
	142	REAL(dp) :: phi_n = 0.0_dp !< latitude of rotated pole of COSMO-DE grid [rad]
	143	REAL(dp) :: lambda_n = 0.0_dp !< longitude of rotaded pole of COSMO-DE grid [rad]
	144	REAL(dp) :: phi_c = 0.0_dp !< rotated-grid latitude of the center of the PALM domain [rad]
	145	REAL(dp) :: lambda_c = 0.0_dp !< rotated-grid longitude of the centre of the PALM domain [rad]
	146	REAL(dp) :: phi_cn = 0.0_dp !< latitude of the rotated pole relative to the COSMO-DE grid [rad]
	147	REAL(dp) :: lambda_cn = 0.0_dp !< longitude of the rotated pole relative to the COSMO-DE grid [rad]
[3395]	148	REAL(dp) :: lam_centre = 0.0_dp !< longitude of the PLAM domain centre in the source (COSMO rotated-pole) system [rad]
	149	REAL(dp) :: phi_centre = 0.0_dp !< latitude of the PLAM domain centre in the source (COSMO rotated-pole) system [rad]
	150	REAL(dp) :: lam_east = 0.0_dp !< longitude of the east central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad]
	151	REAL(dp) :: lam_west = 0.0_dp !< longitude of the west central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad]
	152	REAL(dp) :: phi_north = 0.0_dp !< latitude of the north central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad]
	153	REAL(dp) :: phi_south = 0.0_dp !< latitude of the south central-averaging-domain boundary in the source (COSMO rotated-pole) system [rad]
[3182]	154	REAL(dp) :: gam = 0.0_dp !< angle for working around phirot2phi/rlarot2rla bug
	155	REAL(dp) :: dx = 0.0_dp !< PALM-4U grid spacing in x direction [m]
	156	REAL(dp) :: dy = 0.0_dp !< PALM-4U grid spacing in y direction [m]
	157	REAL(dp) :: dz(10) = -1.0_dp !< PALM-4U grid spacing in z direction [m]
	158	REAL(dp) :: dz_max = 1000.0_dp !< maximum vertical grid spacing [m]
	159	REAL(dp) :: dz_stretch_factor = 1.08_dp !< factor for vertical grid stretching [m]
	160	REAL(dp) :: dz_stretch_level = -9999999.9_dp!< height above which the vertical grid will be stretched [m]
	161	REAL(dp) :: dz_stretch_level_start(9) = -9999999.9_dp !< namelist parameter
	162	REAL(dp) :: dz_stretch_level_end(9) = 9999999.9_dp !< namelist parameter
	163	REAL(dp) :: dz_stretch_factor_array(9) = 1.08_dp !< namelist parameter
	164	REAL(dp) :: dxi = 0.0_dp !< inverse PALM-4U grid spacing in x direction [m^-1]
	165	REAL(dp) :: dyi = 0.0_dp !< inverse PALM-4U grid spacing in y direction [m^-1]
	166	REAL(dp) :: dzi = 0.0_dp !< inverse PALM-4U grid spacing in z direction [m^-1]
[3395]	167	REAL(dp) :: f3 = 0.0_dp !< Coriolis parameter
[3182]	168	REAL(dp) :: lx = 0.0_dp !< PALM-4U domain size in x direction [m]
	169	REAL(dp) :: ly = 0.0_dp !< PALM-4U domain size in y direction [m]
	170	REAL(dp) :: p0 = 0.0_dp !< PALM-4U surface pressure, at z0 [Pa]
	171	REAL(dp) :: x0 = 0.0_dp !< x coordinate of PALM-4U Earth tangent [m]
	172	REAL(dp) :: y0 = 0.0_dp !< y coordinate of PALM-4U Earth tangent [m]
	173	REAL(dp) :: z0 = 0.0_dp !< Elevation of the PALM-4U domain above sea level [m]
	174	REAL(dp) :: z_top = 0.0_dp !< height of the scalar top boundary [m]
	175	REAL(dp) :: zw_top = 0.0_dp !< height of the vertical velocity top boundary [m]
[3613]	176	REAL(dp) :: lonmin_cosmo = 0.0_dp !< Minimunm longitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
	177	REAL(dp) :: lonmax_cosmo = 0.0_dp !< Maximum longitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
	178	REAL(dp) :: latmin_cosmo = 0.0_dp !< Minimunm latitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
	179	REAL(dp) :: latmax_cosmo = 0.0_dp !< Maximum latitude of COSMO-DE's rotated-pole grid [COSMO rotated-pole rad]
	180	REAL(dp) :: lonmin_palm = 0.0_dp !< Minimunm longitude of PALM grid [COSMO rotated-pole rad]
	181	REAL(dp) :: lonmax_palm = 0.0_dp !< Maximum longitude of PALM grid [COSMO rotated-pole rad]
	182	REAL(dp) :: latmin_palm = 0.0_dp !< Minimunm latitude of PALM grid [COSMO rotated-pole rad]
	183	REAL(dp) :: latmax_palm = 0.0_dp !< Maximum latitude of PALM grid [COSMO rotated-pole rad]
[3182]	184	REAL(dp) :: latitude = 0.0_dp !< geographical latitude of the PALM-4U origin, from inipar namelist [deg]
	185	REAL(dp) :: longitude = 0.0_dp !< geographical longitude of the PALM-4U origin, from inipar namelist [deg]
	186	REAL(dp) :: origin_lat = 0.0_dp !< geographical latitude of the PALM-4U origin, from static driver netCDF file [deg]
	187	REAL(dp) :: origin_lon = 0.0_dp !< geographical longitude of the PALM-4U origin, from static driver netCDF file [deg]
	188	REAL(dp) :: rotation_angle = 0.0_dp !< clockwise angle the PALM-4U north is rotated away from geographical north [deg]
	189	REAL(dp) :: end_time = 0.0_dp !< PALM-4U simulation time [s]
[2696]	190
	191	REAL(dp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: hhl !< heights of half layers (cell faces) above sea level in COSMO-DE, read in from external file
	192	REAL(dp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: hfl !< heights of full layers (cell centres) above sea level in COSMO-DE, computed as arithmetic average of hhl
	193	REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET :: depths !< COSMO-DE's TERRA-ML soil layer depths
[3537]	194	REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET :: d_depth !< COSMO-DE's TERRA-ML soil layer thicknesses
	195	REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET :: d_depth_rho_inv !< inverted soil water mass
[2696]	196	REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET :: rlon !< longitudes of COSMO-DE's rotated-pole grid
	197	REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET :: rlat !< latitudes of COSMO-DE's rotated-pole grid
[3182]	198	REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET :: time !< output times
	199	REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET :: x !< base palm grid x coordinate vector pointed to by grid_definitions
	200	REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET :: xu !< base palm grid xu coordinate vector pointed to by grid_definitions
	201	REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET :: y !< base palm grid y coordinate vector pointed to by grid_definitions
	202	REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET :: yv !< base palm grid yv coordinate vector pointed to by grid_definitions
	203	REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET :: z_column !< base palm grid z coordinate vector including the top boundary coordinate (entire column)
	204	REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET :: zw_column !< base palm grid zw coordinate vector including the top boundary coordinate (entire column)
	205	REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET :: z !< base palm grid z coordinate vector pointed to by grid_definitions
	206	REAL(dp), DIMENSION(:), ALLOCATABLE, TARGET :: zw !< base palm grid zw coordinate vector pointed to by grid_definitions
[2696]	207
	208	INTEGER(hp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: soiltyp !< COSMO-DE soil type map
[3182]	209	INTEGER :: dz_stretch_level_end_index(9) !< vertical grid level index until which the vertical grid spacing is stretched
	210	INTEGER :: dz_stretch_level_start_index(9) !< vertical grid level index above which the vertical grid spacing is stretched
[2696]	211	INTEGER :: nt !< number of output time steps
	212	INTEGER :: nx !< number of PALM-4U grid points in x direction
	213	INTEGER :: ny !< number of PALM-4U grid points in y direction
	214	INTEGER :: nz !< number of PALM-4U grid points in z direction
	215	INTEGER :: nlon !< number of longitudal points in target grid (COSMO-DE)
	216	INTEGER :: nlat !< number of latitudal points in target grid (COSMO-DE)
	217	INTEGER :: nlev !< number of levels in target grid (COSMO-DE)
[3537]	218	INTEGER :: ndepths !< number of COSMO-DE soil layers
[2696]	219	INTEGER :: start_hour_flow !< start of flow forcing in number of hours relative to start_date
	220	INTEGER :: start_hour_soil !< start of soil forcing in number of hours relative to start_date, typically equals start_hour_flow
	221	INTEGER :: start_hour_radiation !< start of radiation forcing in number of hours relative to start_date, 0 to 2 hours before start_hour_flow to reconstruct hourly averages from one- to three hourly averages of the input data
	222	INTEGER :: start_hour_soilmoisture !< start of forcing for the soil moisture spin-up in number of hours relative to start_date, typically -672 (-4 weeks)
	223	INTEGER :: end_hour !< simulation time in hours
	224	INTEGER :: end_hour_soilmoisture !< end of soil moisture spin-up in hours relative to start_hour_flow
	225	INTEGER :: step_hour !< number of hours between forcing time steps
	226
[3557]	227	LOGICAL :: init_variables_required !< flag controlling whether init variables are to be processed
	228	LOGICAL :: boundary_variables_required !< flag controlling whether boundary grids are to be allocated and boundary variables are to be computed
	229	LOGICAL :: ls_forcing_variables_required !< flag controlling whether large-scale forcing variables are to be computed
	230	LOGICAL :: surface_forcing_required !< flag controlling whether surface forcing variables are to be computed
[2696]	231
	232	TYPE(nc_var), ALLOCATABLE, TARGET :: input_var_table(:) !< table of input variables
	233	TYPE(nc_var), ALLOCATABLE, TARGET :: output_var_table(:) !< table of input variables
[3557]	234	TYPE(nc_var) :: cosmo_var !< COSMO dummy variable, used for reading HHL, rlon, rlat
[2696]	235
	236	TYPE(grid_definition), TARGET :: palm_grid !< PALM-4U grid in the target system (COSMO-DE rotated-pole)
	237	TYPE(grid_definition), TARGET :: palm_intermediate !< PALM-4U grid with coarse vertical grid wiht levels interpolated from COSMO-DE grid
	238	TYPE(grid_definition), TARGET :: cosmo_grid !< target system (COSMO-DE rotated-pole)
[3557]	239	TYPE(grid_definition), TARGET :: scalars_east_grid !< grid for eastern scalar boundary condition
	240	TYPE(grid_definition), TARGET :: scalars_west_grid !< grid for western scalar boundary condition
	241	TYPE(grid_definition), TARGET :: scalars_north_grid !< grid for northern scalar boundary condition
	242	TYPE(grid_definition), TARGET :: scalars_south_grid !< grid for southern scalar boundary condition
	243	TYPE(grid_definition), TARGET :: scalars_top_grid !< grid for top scalar boundary condition
	244	TYPE(grid_definition), TARGET :: scalars_east_intermediate !< intermediate grid for eastern scalar boundary condition
	245	TYPE(grid_definition), TARGET :: scalars_west_intermediate !< intermediate grid for western scalar boundary condition
	246	TYPE(grid_definition), TARGET :: scalars_north_intermediate !< intermediate grid for northern scalar boundary condition
	247	TYPE(grid_definition), TARGET :: scalars_south_intermediate !< intermediate grid for southern scalar boundary condition
	248	TYPE(grid_definition), TARGET :: scalars_top_intermediate !< intermediate grid for top scalar boundary condition
	249	TYPE(grid_definition), TARGET :: u_initial_grid !< grid for u initial condition
	250	TYPE(grid_definition), TARGET :: u_east_grid !< grid for eastern u boundary condition
	251	TYPE(grid_definition), TARGET :: u_west_grid !< grid for western u boundary condition
	252	TYPE(grid_definition), TARGET :: u_north_grid !< grid for northern u boundary condition
	253	TYPE(grid_definition), TARGET :: u_south_grid !< grid for southern u boundary condition
	254	TYPE(grid_definition), TARGET :: u_top_grid !< grid for top u boundary condition
	255	TYPE(grid_definition), TARGET :: u_initial_intermediate !< intermediate grid for u initial condition
	256	TYPE(grid_definition), TARGET :: u_east_intermediate !< intermediate grid for eastern u boundary condition
	257	TYPE(grid_definition), TARGET :: u_west_intermediate !< intermediate grid for western u boundary condition
	258	TYPE(grid_definition), TARGET :: u_north_intermediate !< intermediate grid for northern u boundary condition
	259	TYPE(grid_definition), TARGET :: u_south_intermediate !< intermediate grid for southern u boundary condition
	260	TYPE(grid_definition), TARGET :: u_top_intermediate !< intermediate grid for top u boundary condition
	261	TYPE(grid_definition), TARGET :: v_initial_grid !< grid for v initial condition
	262	TYPE(grid_definition), TARGET :: v_east_grid !< grid for eastern v boundary condition
	263	TYPE(grid_definition), TARGET :: v_west_grid !< grid for western v boundary condition
	264	TYPE(grid_definition), TARGET :: v_north_grid !< grid for northern v boundary condition
	265	TYPE(grid_definition), TARGET :: v_south_grid !< grid for southern v boundary condition
	266	TYPE(grid_definition), TARGET :: v_top_grid !< grid for top v boundary condition
	267	TYPE(grid_definition), TARGET :: v_initial_intermediate !< intermediate grid for v initial condition
	268	TYPE(grid_definition), TARGET :: v_east_intermediate !< intermediate grid for eastern v boundary condition
	269	TYPE(grid_definition), TARGET :: v_west_intermediate !< intermediate grid for western v boundary condition
	270	TYPE(grid_definition), TARGET :: v_north_intermediate !< intermediate grid for northern v boundary condition
	271	TYPE(grid_definition), TARGET :: v_south_intermediate !< intermediate grid for southern v boundary condition
	272	TYPE(grid_definition), TARGET :: v_top_intermediate !< intermediate grid for top v boundary condition
	273	TYPE(grid_definition), TARGET :: w_initial_grid !< grid for w initial condition
	274	TYPE(grid_definition), TARGET :: w_east_grid !< grid for eastern w boundary condition
	275	TYPE(grid_definition), TARGET :: w_west_grid !< grid for western w boundary condition
	276	TYPE(grid_definition), TARGET :: w_north_grid !< grid for northern w boundary condition
	277	TYPE(grid_definition), TARGET :: w_south_grid !< grid for southern w boundary condition
	278	TYPE(grid_definition), TARGET :: w_top_grid !< grid for top w boundary condition
	279	TYPE(grid_definition), TARGET :: w_initial_intermediate !< intermediate grid for w initial condition
	280	TYPE(grid_definition), TARGET :: w_east_intermediate !< intermediate grid for eastern w boundary condition
	281	TYPE(grid_definition), TARGET :: w_west_intermediate !< intermediate grid for western w boundary condition
	282	TYPE(grid_definition), TARGET :: w_north_intermediate !< intermediate grid for northern w boundary condition
	283	TYPE(grid_definition), TARGET :: w_south_intermediate !< intermediate grid for southern w boundary condition
	284	TYPE(grid_definition), TARGET :: w_top_intermediate !< intermediate grid for top w boundary condition
[3613]	285	TYPE(grid_definition), TARGET :: north_averaged_scalar_profile !< grid of the northern geostrophic scalar averaging region
	286	TYPE(grid_definition), TARGET :: south_averaged_scalar_profile !< grid of the southern geostrophic scalar averaging region
	287	TYPE(grid_definition), TARGET :: west_averaged_scalar_profile !< grid of the western geostrophic scalar averaging region
	288	TYPE(grid_definition), TARGET :: east_averaged_scalar_profile !< grid of the eastern geostrophic scalar averaging region
	289	TYPE(grid_definition), TARGET :: averaged_scalar_profile !< grid of the central geostrophic scalar averaging region
	290	TYPE(grid_definition), TARGET :: averaged_w_profile !< grid of the central geostrophic w-velocity averaging region
	291	TYPE(grid_definition), TARGET :: averaged_initial_scalar_profile !< averaging grid for initial scalar profiles
	292	TYPE(grid_definition), TARGET :: averaged_initial_w_profile !< averaging grid for the initial w profile
[2696]	293
	294	TYPE(io_group), ALLOCATABLE, TARGET :: io_group_list(:) !< List of I/O groups, which group together output variables that share the same input variable
	295
[3182]	296	NAMELIST /inipar/ nx, ny, nz, dx, dy, dz, longitude, latitude, &
[3557]	297	dz_max, dz_stretch_factor, dz_stretch_level, &
	298	dz_stretch_level_start, dz_stretch_level_end
[2696]	299	NAMELIST /d3par/ end_time
	300
[3557]	301	CHARACTER(LEN=LNAME) :: nc_source_text = '' !< Text describing the source of the output data, e.g. 'COSMO-DE analysis from ...'
	302
	303	CHARACTER(LEN=PATH), ALLOCATABLE, DIMENSION(:) :: flow_files !< list of atmospheric input files (<prefix>YYYYMMDDHH-flow.nc)
	304	CHARACTER(LEN=PATH), ALLOCATABLE, DIMENSION(:) :: soil_moisture_files !< list of precipitation input files (<prefix>YYYYMMDDHH-soilmoisture.nc)
	305	CHARACTER(LEN=PATH), ALLOCATABLE, DIMENSION(:) :: soil_files !< list of soil input files (temperature, moisture, <prefix>YYYYMMDDHH-soil.nc)
	306	CHARACTER(LEN=PATH), ALLOCATABLE, DIMENSION(:) :: radiation_files !< list of radiation input files (<prefix>YYYYMMDDHH-rad.nc)
	307
	308	CHARACTER(LEN=SNAME) :: input_prefix !< prefix of input files, e.g. 'laf' for COSMO-DE analyses
	309	CHARACTER(LEN=SNAME) :: flow_prefix !< prefix of flow input files, e.g. 'laf' for COSMO-DE analyses
	310	CHARACTER(LEN=SNAME) :: soil_prefix !< prefix of soil input files, e.g. 'laf' for COSMO-DE analyses
	311	CHARACTER(LEN=SNAME) :: radiation_prefix !< prefix of radiation input files, e.g. 'laf' for COSMO-DE analyses
	312	CHARACTER(LEN=SNAME) :: soilmoisture_prefix !< prefix of input files for soil moisture spin-up, e.g. 'laf' for COSMO-DE analyses
	313	CHARACTER(LEN=SNAME) :: flow_suffix !< suffix of flow input files, e.g. 'flow'
	314	CHARACTER(LEN=SNAME) :: soil_suffix !< suffix of soil input files, e.g. 'soil'
	315	CHARACTER(LEN=SNAME) :: radiation_suffix !< suffix of radiation input files, e.g. 'radiation'
	316	CHARACTER(LEN=SNAME) :: soilmoisture_suffix !< suffix of input files for soil moisture spin-up, e.g. 'soilmoisture'
[2696]	317
[3557]	318	TYPE(nc_file) :: output_file !< metadata of the dynamic driver
[2696]	319
[3557]	320	TYPE(inifor_config) :: cfg !< container of the INIFOR command-line configuration
[3182]	321
[2696]	322	CONTAINS
	323
[3557]	324	!------------------------------------------------------------------------------!
	325	! Description:
	326	! ------------
	327	!> This routine initializes INIFOR. This includes parsing command-line options,
	328	!> setting the names of the input and output files, reading INIFOR's namelist
	329	!> file as well as reading in and setting grid parameters for defining
	330	!> interplation grids later in setup_grids().
	331	!------------------------------------------------------------------------------!
[2696]	332	SUBROUTINE setup_parameters()
[3557]	333	INTEGER :: k !< loop index
[2696]	334
	335	!
	336	!------------------------------------------------------------------------------
	337	! Section 1: Define default parameters
	338	!------------------------------------------------------------------------------
[3182]	339	cfg % start_date = '2013072100'
[2696]	340	end_hour = 2
	341	start_hour_soil = -2
	342	start_hour_soilmoisture = - (4 * 7 * 24) - 2
	343
[3557]	344	!
	345	!-- COSMO-DE and -D2 rotated pole position
[2696]	346	phi_n = 40.0_dp * TO_RADIANS
	347	phi_equat = 50.0_dp * TO_RADIANS
	348	lambda_n = -170.0_dp * TO_RADIANS
	349
[3557]	350	!
	351	!-- Defaultmain centre (_c) of the PALM-4U grid in the geographical system (_g)
[2696]	352	origin_lat = 52.325079_dp * TO_RADIANS ! south-west of Berlin, origin used for the Dec 2017 showcase simulation
	353	origin_lon = 13.082744_dp * TO_RADIANS
[3182]	354	cfg % z0 = 35.0_dp
[2696]	355
[3557]	356	!
	357	!-- Default atmospheric parameters
[3395]	358	cfg % ug = 0.0_dp
	359	cfg % vg = 0.0_dp
[3182]	360	cfg % p0 = P_SL
[2696]	361
[3557]	362	!
	363	!-- Parameters for file names
[2696]	364	start_hour_flow = 0
	365	start_hour_soil = 0
	366	start_hour_radiation = 0
	367	start_hour_soilmoisture = start_hour_flow - 2
	368	end_hour_soilmoisture = start_hour_flow
[3182]	369	step_hour = FORCING_STEP
	370
[3557]	371	input_prefix = 'laf'
[3182]	372	cfg % flow_prefix = input_prefix
[3395]	373	cfg % input_prefix = input_prefix
[3182]	374	cfg % soil_prefix = input_prefix
	375	cfg % radiation_prefix = input_prefix
	376	cfg % soilmoisture_prefix = input_prefix
	377
[2696]	378	flow_suffix = '-flow'
	379	soil_suffix = '-soil'
	380	radiation_suffix = '-rad'
	381	soilmoisture_suffix = '-soilmoisture'
[3395]	382
	383	cfg % debug = .FALSE.
	384	cfg % averaging_angle = 2.0_dp
[2696]	385	!
	386	!------------------------------------------------------------------------------
	387	! Section 2: Read command-line arguments, namelist, and grid parameters
	388	!------------------------------------------------------------------------------
	389
[3557]	390	!
	391	!-- Set default paths and modes
[3182]	392	cfg % input_path = './'
	393	cfg % hhl_file = ''
	394	cfg % soiltyp_file = ''
	395	cfg % namelist_file = './namelist'
	396	cfg % static_driver_file = ''
	397	cfg % output_file = './palm-4u-input.nc'
	398	cfg % ic_mode = 'volume'
	399	cfg % bc_mode = 'real'
[3395]	400	cfg % averaging_mode = 'level'
[2696]	401
[3557]	402	!
	403	!-- Overwrite defaults with user configuration
[3182]	404	CALL parse_command_line_arguments( cfg )
[2696]	405
[3395]	406	flow_prefix = TRIM(cfg % input_prefix)
	407	radiation_prefix = TRIM(cfg % input_prefix)
	408	soil_prefix = TRIM(cfg % input_prefix)
	409	soilmoisture_prefix = TRIM(cfg % input_prefix)
	410	IF (cfg % flow_prefix_is_set) flow_prefix = TRIM(cfg % flow_prefix)
	411	IF (cfg % radiation_prefix_is_set) radiation_prefix = TRIM(cfg % radiation_prefix)
	412	IF (cfg % soil_prefix_is_set) soil_prefix = TRIM(cfg % soil_prefix)
	413	IF (cfg % soilmoisture_prefix_is_set) soilmoisture_prefix = TRIM(cfg % soilmoisture_prefix)
	414
[3182]	415	output_file % name = cfg % output_file
	416	z0 = cfg % z0
	417	p0 = cfg % p0
[2696]	418
[3182]	419	init_variables_required = .TRUE.
	420	boundary_variables_required = TRIM( cfg % bc_mode ) == 'real'
	421	ls_forcing_variables_required = TRIM( cfg % bc_mode ) == 'ideal'
[3456]	422	surface_forcing_required = .FALSE.
[2696]	423
[3182]	424	IF ( ls_forcing_variables_required ) THEN
	425	message = "Averaging of large-scale forcing profiles " // &
	426	"has not been implemented, yet."
[3615]	427	CALL inifor_abort('setup_parameters', message)
[3182]	428	END IF
[2696]	429
[3557]	430	!
	431	!-- Set default file paths, if not specified by user.
[3182]	432	CALL normalize_path(cfg % input_path)
	433	IF (TRIM(cfg % hhl_file) == '') cfg % hhl_file = TRIM(cfg % input_path) // 'hhl.nc'
	434	IF (TRIM(cfg % soiltyp_file) == '') cfg % soiltyp_file = TRIM(cfg % input_path) // 'soil.nc'
	435
	436	CALL validate_config( cfg )
	437
[3613]	438	CALL report('main_loop', 'Running INIFOR version ' // VERSION)
[3182]	439	CALL report('setup_parameters', "initialization mode: " // TRIM(cfg % ic_mode))
	440	CALL report('setup_parameters', " forcing mode: " // TRIM(cfg % bc_mode))
[3537]	441	CALL report('setup_parameters', " averaging mode: " // TRIM(cfg % averaging_mode))
[3182]	442	CALL report('setup_parameters', " data path: " // TRIM(cfg % input_path))
	443	CALL report('setup_parameters', " hhl file: " // TRIM(cfg % hhl_file))
	444	CALL report('setup_parameters', " soiltyp file: " // TRIM(cfg % soiltyp_file))
	445	CALL report('setup_parameters', " namelist file: " // TRIM(cfg % namelist_file))
	446	CALL report('setup_parameters', " output data file: " // TRIM(output_file % name))
[3537]	447	IF (cfg % debug ) CALL report('setup_parameters', " debugging mode: enabled")
[3182]	448
[2696]	449	CALL run_control('time', 'init')
[3557]	450	!
	451	!-- Read in namelist parameters
[3182]	452	OPEN(10, FILE=cfg % namelist_file)
[2696]	453	READ(10, NML=inipar) ! nx, ny, nz, dx, dy, dz
	454	READ(10, NML=d3par) ! end_time
	455	CLOSE(10)
	456	CALL run_control('time', 'read')
	457
[3182]	458	end_hour = CEILING( end_time / 3600.0 * step_hour )
[2696]	459
[3557]	460	!
	461	!-- Generate input file lists
[3395]	462	CALL get_input_file_list( &
	463	cfg % start_date, start_hour_flow, end_hour, step_hour, &
	464	cfg % input_path, flow_prefix, flow_suffix, flow_files)
	465	CALL get_input_file_list( &
	466	cfg % start_date, start_hour_soil, end_hour, step_hour, &
	467	cfg % input_path, soil_prefix, soil_suffix, soil_files)
	468	CALL get_input_file_list( &
	469	cfg % start_date, start_hour_radiation, end_hour, step_hour, &
	470	cfg % input_path, radiation_prefix, radiation_suffix, radiation_files)
	471	CALL get_input_file_list( &
	472	cfg % start_date, start_hour_soilmoisture, end_hour_soilmoisture, step_hour, &
	473	cfg % input_path, soilmoisture_prefix, soilmoisture_suffix, soil_moisture_files)
[2696]	474
	475	!
	476	!------------------------------------------------------------------------------
	477	! Section 3: Check for consistency
	478	!------------------------------------------------------------------------------
[3182]	479
[2696]	480	!
	481	!------------------------------------------------------------------------------
	482	! Section 4: Compute additional parameters
	483	!------------------------------------------------------------------------------
	484	!------------------------------------------------------------------------------
	485	! Section 4.1: COSMO-DE parameters
	486	!------------------------------------------------------------------------------
	487
	488
	489	CALL run_control('time', 'init')
[3557]	490	!
	491	!-- Read COSMO soil type map
[2696]	492	cosmo_var % name = 'SOILTYP'
[3182]	493	CALL get_netcdf_variable(cfg % soiltyp_file, cosmo_var, soiltyp)
[2696]	494
	495	message = 'Reading PALM-4U origin from'
[3182]	496	IF (TRIM(cfg % static_driver_file) .NE. '') THEN
[2696]	497
[3182]	498	origin_lon = get_netcdf_attribute(cfg % static_driver_file, 'origin_lon')
	499	origin_lat = get_netcdf_attribute(cfg % static_driver_file, 'origin_lat')
[2696]	500
	501	message = TRIM(message) // " static driver file '" &
[3182]	502	// TRIM(cfg % static_driver_file) // "'"
[2696]	503
	504
	505	ELSE
	506
	507	origin_lon = longitude
	508	origin_lat = latitude
	509
	510	message = TRIM(message) // " namlist file '" &
[3182]	511	// TRIM(cfg % namelist_file) // "'"
[2696]	512
	513	END IF
	514	origin_lon = origin_lon * TO_RADIANS
	515	origin_lat = origin_lat * TO_RADIANS
	516
	517	CALL report('setup_parameters', message)
	518
[3557]	519	!
	520	!-- Read in COSMO heights of half layers (vertical cell faces)
[2696]	521	cosmo_var % name = 'HHL'
[3182]	522	CALL get_netcdf_variable(cfg % hhl_file, cosmo_var, hhl)
[3537]	523	CALL get_netcdf_dim_vector(cfg % hhl_file, 'rlon', rlon)
	524	CALL get_netcdf_dim_vector(cfg % hhl_file, 'rlat', rlat)
	525	CALL get_netcdf_dim_vector(soil_files(1), 'depth_2', depths)
[2696]	526	CALL run_control('time', 'read')
	527
	528	CALL reverse(hhl)
	529	nlon = SIZE(hhl, 1)
	530	nlat = SIZE(hhl, 2)
	531	nlev = SIZE(hhl, 3)
[3537]	532	ndepths = SIZE(depths)
[2696]	533
[3613]	534	lonmin_cosmo = MINVAL(rlon) * TO_RADIANS
	535	lonmax_cosmo = MAXVAL(rlon) * TO_RADIANS
	536	latmin_cosmo = MINVAL(rlat) * TO_RADIANS
	537	latmax_cosmo = MAXVAL(rlat) * TO_RADIANS
[2696]	538	CALL run_control('time', 'comp')
	539
[3557]	540	!
	541	!-- Appoximate COSMO-DE heights of full layers (cell centres)
[2696]	542	ALLOCATE( hfl(nlon, nlat, nlev-1) )
[3537]	543	ALLOCATE( d_depth(ndepths), d_depth_rho_inv(ndepths) )
	544
[2696]	545	CALL run_control('time', 'alloc')
[3537]	546	CALL get_soil_layer_thickness( depths, d_depth )
	547	d_depth_rho_inv = 1.0_dp / ( d_depth * RHO_L )
	548
[3557]	549	!
	550	!-- Appoximate COSMO-DE heights of full layers (cell centres)
[2696]	551	DO k = 1, nlev-1
	552	hfl(:,:,k) = 0.5_dp * ( hhl(:,:,k) + &
	553	hhl(:,:,k+1) )
	554	END DO
[3537]	555	CALL run_control('time', 'comp')
[2696]	556
[3537]	557
	558
[2696]	559	!------------------------------------------------------------------------------
	560	! Section 4.2: PALM-4U parameters
	561	!------------------------------------------------------------------------------
[3557]	562	!
	563	!-- PALM-4U domain extents
[2696]	564	lx = (nx+1) * dx
	565	ly = (ny+1) * dy
	566
[3557]	567	!
	568	!-- PALM-4U point of Earth tangency
[2696]	569	x0 = 0.0_dp
	570	y0 = 0.0_dp
	571
[3557]	572	!
	573	!-- time vector
[3182]	574	nt = CEILING(end_time / (step_hour * 3600.0_dp)) + 1
[2696]	575	ALLOCATE( time(nt) )
	576	CALL linspace(0.0_dp, 3600.0_dp * (nt-1), time)
	577	output_file % time => time
	578	CALL run_control('time', 'init')
	579
[3557]	580	!
	581	!-- Convert the PALM-4U origin coordinates to COSMO's rotated-pole grid
[2696]	582	phi_c = TO_RADIANS * &
	583	phi2phirot( origin_lat * TO_DEGREES, origin_lon * TO_DEGREES,&
	584	phi_n * TO_DEGREES, lambda_n * TO_DEGREES )
	585	lambda_c = TO_RADIANS * &
	586	rla2rlarot( origin_lat * TO_DEGREES, origin_lon * TO_DEGREES,&
	587	phi_n * TO_DEGREES, lambda_n * TO_DEGREES, &
	588	0.0_dp )
	589
[3557]	590	!
	591	!-- Set gamma according to whether PALM domain is in the northern or southern
	592	!-- hemisphere of the COSMO rotated-pole system. Gamma assumes either the
	593	!-- value 0 or PI and is needed to work around around a bug in the
	594	!-- rotated-pole coordinate transformations.
[2696]	595	gam = gamma_from_hemisphere(origin_lat, phi_equat)
	596
[3557]	597	!
	598	!-- Compute the north pole of the rotated-pole grid centred at the PALM-4U
	599	!-- domain centre. The resulting (phi_cn, lambda_cn) are coordinates in
	600	!-- COSMO-DE's rotated-pole grid.
[2696]	601	phi_cn = phic_to_phin(phi_c)
	602	lambda_cn = lamc_to_lamn(phi_c, lambda_c)
	603
	604	message = "PALM-4U origin:" // NEW_LINE('') // &
	605	" lon (lambda) = " // &
	606	TRIM(real_to_str_f(origin_lon * TO_DEGREES)) // " deg"// NEW_LINE(' ') //&
	607	" lat (phi ) = " // &
	608	TRIM(real_to_str_f(origin_lat * TO_DEGREES)) // " deg (geographical)" // NEW_LINE(' ') //&
	609	" lon (lambda) = " // &
	610	TRIM(real_to_str_f(lambda_c * TO_DEGREES)) // " deg" // NEW_LINE(' ') // &
	611	" lat (phi ) = " // &
	612	TRIM(real_to_str_f(phi_c * TO_DEGREES)) // " deg (COSMO-DE rotated-pole)"
	613	CALL report ('setup_parameters', message)
	614
	615	message = "North pole of the rotated COSMO-DE system:" // NEW_LINE(' ') // &
	616	" lon (lambda) = " // &
	617	TRIM(real_to_str_f(lambda_n * TO_DEGREES)) // " deg" // NEW_LINE(' ') //&
	618	" lat (phi ) = " // &
	619	TRIM(real_to_str_f(phi_n * TO_DEGREES)) // " deg (geographical)"
	620	CALL report ('setup_parameters', message)
	621
	622	message = "North pole of the rotated palm system:" // NEW_LINE(' ') // &
	623	" lon (lambda) = " // &
	624	TRIM(real_to_str_f(lambda_cn * TO_DEGREES)) // " deg" // NEW_LINE(' ') // &
	625	" lat (phi ) = " // &
	626	TRIM(real_to_str_f(phi_cn * TO_DEGREES)) // " deg (COSMO-DE rotated-pole)"
	627	CALL report ('setup_parameters', message)
	628
	629	CALL run_control('time', 'comp')
	630
[3395]	631	!------------------------------------------------------------------------------
	632	! Section 4.3: INIFOR averaging domains
	633	!------------------------------------------------------------------------------
	634
[3557]	635	!
	636	!-- Compute coordiantes of the PALM centre in the source (COSMO) system
[3395]	637	phi_centre = phirot2phi( &
	638	phirot = project(0.5_dply, y0, EARTH_RADIUS) TO_DEGREES, &
	639	rlarot = project(0.5_dplx, x0, EARTH_RADIUS) TO_DEGREES, &
[3618]	640	polphi = phi_cn * TO_DEGREES, &
[3395]	641	polgam = gam * TO_DEGREES &
	642	) * TO_RADIANS
	643
	644	lam_centre = rlarot2rla( &
	645	phirot = project(0.5_dply, y0, EARTH_RADIUS) TO_DEGREES, &
	646	rlarot = project(0.5_dplx, x0, EARTH_RADIUS) TO_DEGREES, &
	647	polphi = phi_cn * TO_DEGREES, pollam = lambda_cn * TO_DEGREES, &
	648	polgam = gam * TO_DEGREES &
	649	) * TO_RADIANS
	650
	651	message = "PALM-4U centre:" // NEW_LINE('') // &
	652	" lon (lambda) = " // &
	653	TRIM(real_to_str_f(lam_centre * TO_DEGREES)) // " deg" // NEW_LINE(' ') // &
	654	" lat (phi ) = " // &
	655	TRIM(real_to_str_f(phi_centre * TO_DEGREES)) // " deg (COSMO-DE rotated-pole)"
	656	CALL report( 'setup_parameters', message )
	657
[3557]	658	!
	659	!-- Compute boundaries of the central averaging box
[3395]	660	averaging_angle = cfg % averaging_angle * TO_RADIANS
	661	lam_east = lam_centre + 0.5_dp * averaging_angle
	662	lam_west = lam_centre - 0.5_dp * averaging_angle
	663	phi_north = phi_centre + 0.5_dp * averaging_angle
	664	phi_south = phi_centre - 0.5_dp * averaging_angle
	665	averaging_width_ew = averaging_angle * COS(phi_centre) * EARTH_RADIUS
	666	averaging_width_ns = averaging_angle * EARTH_RADIUS
	667
[3557]	668	!
	669	!-- Coriolis parameter
[3395]	670	f3 = 2.0_dp * OMEGA * SIN( &
	671	TO_RADIANSphirot2phi( phi_centre TO_DEGREES, lam_centre * TO_DEGREES,&
[3618]	672	phi_n * TO_DEGREES, &
[3395]	673	gam * TO_DEGREES ) &
	674	)
	675
[2696]	676	END SUBROUTINE setup_parameters
	677
	678
	679	!------------------------------------------------------------------------------!
	680	! Description:
	681	! ------------
[3557]	682	!> Defines the COSMO, PALM-4U, PALM-4U boundary grids, in partucular their
	683	!> coordinates and interpolation weights
[2696]	684	!------------------------------------------------------------------------------!
[3557]	685	SUBROUTINE setup_grids()
[2696]	686	CHARACTER :: interp_mode
[3182]	687
[2696]	688	!------------------------------------------------------------------------------
[3182]	689	! Section 0: Define base PALM-4U grid coordinate vectors
[2696]	690	!------------------------------------------------------------------------------
[3557]	691	!
	692	!-- palm x y z, we allocate the column to nz+1 in order to include the top
	693	!-- scalar boundary. The interpolation grids will be associated with
	694	!-- a shorter column that omits the top element.
[3182]	695	ALLOCATE( x(0:nx), y(0:ny), z(1:nz), z_column(1:nz+1) )
	696	CALL linspace(0.5_dp * dx, lx - 0.5_dp * dx, x)
	697	CALL linspace(0.5_dp * dy, ly - 0.5_dp * dy, y)
	698	CALL stretched_z(z_column, dz, dz_max=dz_max, &
	699	dz_stretch_factor=dz_stretch_factor, &
	700	dz_stretch_level=dz_stretch_level, &
	701	dz_stretch_level_start=dz_stretch_level_start, &
	702	dz_stretch_level_end=dz_stretch_level_end, &
	703	dz_stretch_factor_array=dz_stretch_factor_array)
	704	z(1:nz) = z_column(1:nz)
	705	z_top = z_column(nz+1)
	706
[3557]	707	!
	708	!-- palm xu yv zw, compared to the scalar grid, velocity coordinates
	709	!-- contain one element less.
[3182]	710	ALLOCATE( xu(1:nx), yv(1:ny), zw(1:nz-1), zw_column(1:nz))
	711	CALL linspace(dx, lx - dx, xu)
	712	CALL linspace(dy, ly - dy, yv)
	713	CALL midpoints(z_column, zw_column)
	714	zw(1:nz-1) = zw_column(1:nz-1)
	715	zw_top = zw_column(nz)
	716
	717
	718	!------------------------------------------------------------------------------
	719	! Section 1: Define initialization and boundary grids
	720	!------------------------------------------------------------------------------
[2696]	721	CALL init_grid_definition('palm', grid=palm_grid, &
	722	xmin=0.0_dp, xmax=lx, &
	723	ymin=0.0_dp, ymax=ly, &
[3613]	724	x0=x0, y0=y0, z0=z0, &
[3182]	725	nx=nx, ny=ny, nz=nz, z=z, zw=zw, ic_mode=cfg % ic_mode)
[2696]	726
[3557]	727	!
	728	!-- Subtracting 1 because arrays will be allocated with nlon + 1 elements.
[2696]	729	CALL init_grid_definition('cosmo-de', grid=cosmo_grid, &
[3613]	730	xmin=lonmin_cosmo, xmax=lonmax_cosmo, &
	731	ymin=latmin_cosmo, ymax=latmax_cosmo, &
	732	x0=x0, y0=y0, z0=0.0_dp, &
[2696]	733	nx=nlon-1, ny=nlat-1, nz=nlev-1)
	734
[3557]	735	!
	736	!-- Define intermediate grid. This is the same as palm_grid except with a
	737	!-- much coarser vertical grid. The vertical levels are interpolated in each
	738	!-- PALM column from COSMO's secondary levels. The main levels are then
	739	!-- computed as the averages of the bounding secondary levels.
[2696]	740	CALL init_grid_definition('palm intermediate', grid=palm_intermediate, &
	741	xmin=0.0_dp, xmax=lx, &
	742	ymin=0.0_dp, ymax=ly, &
[3613]	743	x0=x0, y0=y0, z0=z0, &
[2696]	744	nx=nx, ny=ny, nz=nlev-2)
	745
	746	CALL init_grid_definition('boundary', grid=u_initial_grid, &
	747	xmin = dx, xmax = lx - dx, &
	748	ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, &
	749	x0=x0, y0=y0, z0 = z0, &
	750	nx = nx-1, ny = ny, nz = nz, &
[3182]	751	z=z, ic_mode=cfg % ic_mode)
[2696]	752
	753	CALL init_grid_definition('boundary', grid=v_initial_grid, &
	754	xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, &
	755	ymin = dy, ymax = ly - dy, &
	756	x0=x0, y0=y0, z0 = z0, &
	757	nx = nx, ny = ny-1, nz = nz, &
[3182]	758	z=z, ic_mode=cfg % ic_mode)
[2696]	759
	760	CALL init_grid_definition('boundary', grid=w_initial_grid, &
	761	xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, &
	762	ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, &
	763	x0=x0, y0=y0, z0 = z0, &
	764	nx = nx, ny = ny, nz = nz-1, &
[3182]	765	z=zw, ic_mode=cfg % ic_mode)
[2696]	766
	767	CALL init_grid_definition('boundary intermediate', grid=u_initial_intermediate, &
	768	xmin = dx, xmax = lx - dx, &
	769	ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, &
	770	x0=x0, y0=y0, z0 = z0, &
[3182]	771	nx = nx-1, ny = ny, nz = nlev - 2)
[2696]	772
	773	CALL init_grid_definition('boundary intermediate', grid=v_initial_intermediate, &
	774	xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, &
	775	ymin = dy, ymax = ly - dy, &
	776	x0=x0, y0=y0, z0 = z0, &
[3182]	777	nx = nx, ny = ny-1, nz = nlev - 2)
[2696]	778
	779	CALL init_grid_definition('boundary intermediate', grid=w_initial_intermediate, &
	780	xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, &
	781	ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, &
	782	x0=x0, y0=y0, z0 = z0, &
[3395]	783	nx = nx, ny = ny, nz = nlev - 1)
[2696]	784
	785	IF (boundary_variables_required) THEN
	786	!
	787	!------------------------------------------------------------------------------
	788	! Section 2: Define PALM-4U boundary grids
	789	!------------------------------------------------------------------------------
	790	CALL init_grid_definition('boundary', grid=scalars_east_grid, &
	791	xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx, &
	792	ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, &
	793	x0=x0, y0=y0, z0 = z0, &
[3182]	794	nx = 0, ny = ny, nz = nz, z=z)
[2696]	795
	796	CALL init_grid_definition('boundary', grid=scalars_west_grid, &
	797	xmin = -0.5_dp * dx, xmax = -0.5_dp * dx, &
	798	ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, &
	799	x0=x0, y0=y0, z0 = z0, &
[3182]	800	nx = 0, ny = ny, nz = nz, z=z)
[2696]	801
	802	CALL init_grid_definition('boundary', grid=scalars_north_grid, &
	803	xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, &
	804	ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy, &
	805	x0=x0, y0=y0, z0 = z0, &
[3182]	806	nx = nx, ny = 0, nz = nz, z=z)
[2696]	807
	808	CALL init_grid_definition('boundary', grid=scalars_south_grid, &
	809	xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, &
	810	ymin = -0.5_dp * dy, ymax = -0.5_dp * dy, &
	811	x0=x0, y0=y0, z0 = z0, &
[3182]	812	nx = nx, ny = 0, nz = nz, z=z)
[2696]	813
	814	CALL init_grid_definition('boundary', grid=scalars_top_grid, &
	815	xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, &
	816	ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, &
	817	x0=x0, y0=y0, z0 = z0, &
[3182]	818	nx = nx, ny = ny, nz = 1, z=(/z_top/))
[2696]	819
	820	CALL init_grid_definition('boundary', grid=u_east_grid, &
	821	xmin = lx, xmax = lx, &
	822	ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, &
	823	x0=x0, y0=y0, z0 = z0, &
[3182]	824	nx = 0, ny = ny, nz = nz, z=z)
[2696]	825
	826	CALL init_grid_definition('boundary', grid=u_west_grid, &
	827	xmin = 0.0_dp, xmax = 0.0_dp, &
	828	ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, &
	829	x0=x0, y0=y0, z0 = z0, &
[3182]	830	nx = 0, ny = ny, nz = nz, z=z)
[2696]	831
	832	CALL init_grid_definition('boundary', grid=u_north_grid, &
	833	xmin = dx, xmax = lx - dx, &
	834	ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy, &
	835	x0=x0, y0=y0, z0 = z0, &
[3182]	836	nx = nx-1, ny = 0, nz = nz, z=z)
	837
[2696]	838	CALL init_grid_definition('boundary', grid=u_south_grid, &
	839	xmin = dx, xmax = lx - dx, &
	840	ymin = -0.5_dp * dy, ymax = -0.5_dp * dy, &
	841	x0=x0, y0=y0, z0 = z0, &
[3182]	842	nx = nx-1, ny = 0, nz = nz, z=z)
[2696]	843
	844	CALL init_grid_definition('boundary', grid=u_top_grid, &
	845	xmin = dx, xmax = lx - dx, &
	846	ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, &
	847	x0=x0, y0=y0, z0 = z0, &
[3182]	848	nx = nx-1, ny = ny, nz = 1, z=(/z_top/))
[2696]	849
	850	CALL init_grid_definition('boundary', grid=v_east_grid, &
	851	xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx, &
	852	ymin = dy, ymax = ly - dy, &
	853	x0=x0, y0=y0, z0 = z0, &
[3182]	854	nx = 0, ny = ny-1, nz = nz, z=z)
[2696]	855
	856	CALL init_grid_definition('boundary', grid=v_west_grid, &
	857	xmin = -0.5_dp * dx, xmax = -0.5_dp * dx, &
	858	ymin = dy, ymax = ly - dy, &
	859	x0=x0, y0=y0, z0 = z0, &
[3182]	860	nx = 0, ny = ny-1, nz = nz, z=z)
[2696]	861
	862	CALL init_grid_definition('boundary', grid=v_north_grid, &
	863	xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, &
	864	ymin = ly, ymax = ly, &
	865	x0=x0, y0=y0, z0 = z0, &
[3182]	866	nx = nx, ny = 0, nz = nz, z=z)
[2696]	867
	868	CALL init_grid_definition('boundary', grid=v_south_grid, &
	869	xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, &
	870	ymin = 0.0_dp, ymax = 0.0_dp, &
	871	x0=x0, y0=y0, z0 = z0, &
[3182]	872	nx = nx, ny = 0, nz = nz, z=z)
[2696]	873
	874	CALL init_grid_definition('boundary', grid=v_top_grid, &
	875	xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, &
	876	ymin = dy, ymax = ly - dy, &
	877	x0=x0, y0=y0, z0 = z0, &
[3182]	878	nx = nx, ny = ny-1, nz = 1, z=(/z_top/))
[2696]	879
	880	CALL init_grid_definition('boundary', grid=w_east_grid, &
	881	xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx, &
	882	ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, &
	883	x0=x0, y0=y0, z0 = z0, &
[3182]	884	nx = 0, ny = ny, nz = nz - 1, z=zw)
[2696]	885
	886	CALL init_grid_definition('boundary', grid=w_west_grid, &
	887	xmin = -0.5_dp * dx, xmax = -0.5_dp * dx, &
	888	ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, &
	889	x0=x0, y0=y0, z0 = z0, &
[3182]	890	nx = 0, ny = ny, nz = nz - 1, z=zw)
[2696]	891
	892	CALL init_grid_definition('boundary', grid=w_north_grid, &
	893	xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, &
	894	ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy, &
	895	x0=x0, y0=y0, z0 = z0, &
[3182]	896	nx = nx, ny = 0, nz = nz - 1, z=zw)
[2696]	897
	898	CALL init_grid_definition('boundary', grid=w_south_grid, &
	899	xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, &
	900	ymin = -0.5_dp * dy, ymax = -0.5_dp * dy, &
	901	x0=x0, y0=y0, z0 = z0, &
[3182]	902	nx = nx, ny = 0, nz = nz - 1, z=zw)
[2696]	903
	904	CALL init_grid_definition('boundary', grid=w_top_grid, &
	905	xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, &
	906	ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, &
	907	x0=x0, y0=y0, z0 = z0, &
[3182]	908	nx = nx, ny = ny, nz = 1, z=(/zw_top/))
[2696]	909
	910	CALL init_grid_definition('boundary intermediate', grid=scalars_east_intermediate, &
	911	xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx, &
	912	ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, &
	913	x0=x0, y0=y0, z0 = z0, &
[3182]	914	nx = 0, ny = ny, nz = nlev - 2)
[2696]	915
	916	CALL init_grid_definition('boundary intermediate', grid=scalars_west_intermediate, &
	917	xmin = -0.5_dp * dx, xmax = -0.5_dp * dx, &
	918	ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, &
	919	x0=x0, y0=y0, z0 = z0, &
[3182]	920	nx = 0, ny = ny, nz = nlev - 2)
[2696]	921
	922	CALL init_grid_definition('boundary intermediate', grid=scalars_north_intermediate, &
	923	xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, &
	924	ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy, &
	925	x0=x0, y0=y0, z0 = z0, &
[3182]	926	nx = nx, ny = 0, nz = nlev - 2)
[2696]	927
	928	CALL init_grid_definition('boundary intermediate', grid=scalars_south_intermediate, &
	929	xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, &
	930	ymin = -0.5_dp * dy, ymax = -0.5_dp * dy, &
	931	x0=x0, y0=y0, z0 = z0, &
[3182]	932	nx = nx, ny = 0, nz = nlev - 2)
[2696]	933
	934	CALL init_grid_definition('boundary intermediate', grid=scalars_top_intermediate, &
	935	xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, &
	936	ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, &
	937	x0=x0, y0=y0, z0 = z0, &
[3182]	938	nx = nx, ny = ny, nz = nlev - 2)
[2696]	939
	940	CALL init_grid_definition('boundary intermediate', grid=u_east_intermediate, &
	941	xmin = lx, xmax = lx, &
	942	ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, &
	943	x0=x0, y0=y0, z0 = z0, &
[3182]	944	nx = 0, ny = ny, nz = nlev - 2)
[2696]	945
	946	CALL init_grid_definition('boundary intermediate', grid=u_west_intermediate, &
	947	xmin = 0.0_dp, xmax = 0.0_dp, &
	948	ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, &
	949	x0=x0, y0=y0, z0 = z0, &
[3182]	950	nx = 0, ny = ny, nz = nlev - 2)
[2696]	951
	952	CALL init_grid_definition('boundary intermediate', grid=u_north_intermediate, &
	953	xmin = dx, xmax = lx - dx, &
	954	ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy, &
	955	x0=x0, y0=y0, z0 = z0, &
[3182]	956	nx = nx-1, ny = 0, nz = nlev - 2)
[2696]	957
	958	CALL init_grid_definition('boundary intermediate', grid=u_south_intermediate, &
	959	xmin = dx, xmax = lx - dx, &
	960	ymin = -0.5_dp * dy, ymax = -0.5_dp * dy, &
	961	x0=x0, y0=y0, z0 = z0, &
[3182]	962	nx = nx-1, ny = 0, nz = nlev - 2)
[2696]	963
	964	CALL init_grid_definition('boundary intermediate', grid=u_top_intermediate, &
	965	xmin = dx, xmax = lx - dx, &
	966	ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, &
	967	x0=x0, y0=y0, z0 = z0, &
[3182]	968	nx = nx-1, ny = ny, nz = nlev - 2)
[2696]	969
	970	CALL init_grid_definition('boundary intermediate', grid=v_east_intermediate, &
	971	xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx, &
	972	ymin = dy, ymax = ly - dy, &
	973	x0=x0, y0=y0, z0 = z0, &
[3182]	974	nx = 0, ny = ny-1, nz = nlev - 2)
[2696]	975
	976	CALL init_grid_definition('boundary intermediate', grid=v_west_intermediate, &
	977	xmin = -0.5_dp * dx, xmax = -0.5_dp * dx, &
	978	ymin = dy, ymax = ly - dy, &
	979	x0=x0, y0=y0, z0 = z0, &
[3182]	980	nx = 0, ny = ny-1, nz = nlev - 2)
[2696]	981
	982	CALL init_grid_definition('boundary intermediate', grid=v_north_intermediate, &
	983	xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, &
	984	ymin = ly, ymax = ly, &
	985	x0=x0, y0=y0, z0 = z0, &
[3182]	986	nx = nx, ny = 0, nz = nlev - 2)
[2696]	987
	988	CALL init_grid_definition('boundary intermediate', grid=v_south_intermediate, &
	989	xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, &
	990	ymin = 0.0_dp, ymax = 0.0_dp, &
	991	x0=x0, y0=y0, z0 = z0, &
[3182]	992	nx = nx, ny = 0, nz = nlev - 2)
[2696]	993
	994	CALL init_grid_definition('boundary intermediate', grid=v_top_intermediate, &
	995	xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, &
	996	ymin = dy, ymax = ly - dy, &
	997	x0=x0, y0=y0, z0 = z0, &
[3182]	998	nx = nx, ny = ny-1, nz = nlev - 2)
[2696]	999
	1000	CALL init_grid_definition('boundary intermediate', grid=w_east_intermediate, &
	1001	xmin = lx + 0.5_dp * dx, xmax = lx + 0.5_dp * dx, &
	1002	ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, &
	1003	x0=x0, y0=y0, z0 = z0, &
[3395]	1004	nx = 0, ny = ny, nz = nlev - 1)
[2696]	1005
	1006	CALL init_grid_definition('boundary intermediate', grid=w_west_intermediate, &
	1007	xmin = -0.5_dp * dx, xmax = -0.5_dp * dx, &
	1008	ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, &
	1009	x0=x0, y0=y0, z0 = z0, &
[3395]	1010	nx = 0, ny = ny, nz = nlev - 1)
[2696]	1011
	1012	CALL init_grid_definition('boundary intermediate', grid=w_north_intermediate, &
	1013	xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, &
	1014	ymin = ly + 0.5_dp * dy, ymax = ly + 0.5_dp * dy, &
	1015	x0=x0, y0=y0, z0 = z0, &
[3395]	1016	nx = nx, ny = 0, nz = nlev - 1)
[2696]	1017
	1018	CALL init_grid_definition('boundary intermediate', grid=w_south_intermediate, &
	1019	xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, &
	1020	ymin = -0.5_dp * dy, ymax = -0.5_dp * dy, &
	1021	x0=x0, y0=y0, z0 = z0, &
[3395]	1022	nx = nx, ny = 0, nz = nlev - 1)
[2696]	1023
	1024	CALL init_grid_definition('boundary intermediate', grid=w_top_intermediate, &
	1025	xmin = 0.5_dp * dx, xmax = lx - 0.5_dp * dx, &
	1026	ymin = 0.5_dp * dy, ymax = ly - 0.5_dp * dy, &
	1027	x0=x0, y0=y0, z0 = z0, &
[3395]	1028	nx = nx, ny = ny, nz = nlev - 1)
[2696]	1029	END IF
	1030
	1031	!
	1032	!------------------------------------------------------------------------------
	1033	! Section 3: Define profile grids
	1034	!------------------------------------------------------------------------------
	1035
[3613]	1036	lonmin_palm = MINVAL(palm_intermediate % clon)
	1037	lonmax_palm = MAXVAL(palm_intermediate % clon)
	1038	latmin_palm = MINVAL(palm_intermediate % clat)
	1039	latmax_palm = MAXVAL(palm_intermediate % clat)
	1040
	1041	CALL init_averaging_grid(averaged_initial_scalar_profile, cosmo_grid, &
	1042	x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0, &
	1043	lonmin = lonmin_palm, lonmax = lonmax_palm, &
	1044	latmin = latmin_palm, latmax = latmax_palm, &
	1045	kind='scalar')
	1046
	1047	CALL init_averaging_grid(averaged_initial_w_profile, cosmo_grid, &
	1048	x = 0.5_dp * lx, y = 0.5_dp * ly, z = zw, z0 = z0, &
	1049	lonmin = lonmin_palm, lonmax = lonmax_palm, &
	1050	latmin = latmin_palm, latmax = latmax_palm, &
	1051	kind='w')
	1052
[3395]	1053	CALL init_averaging_grid(averaged_scalar_profile, cosmo_grid, &
	1054	x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0, &
	1055	lonmin = lam_west, lonmax = lam_east, &
	1056	latmin = phi_south, latmax = phi_north, &
	1057	kind='scalar')
	1058
	1059	CALL init_averaging_grid(averaged_w_profile, cosmo_grid, &
	1060	x = 0.5_dp * lx, y = 0.5_dp * ly, z = zw, z0 = z0, &
	1061	lonmin = lam_west, lonmax = lam_east, &
	1062	latmin = phi_south, latmax = phi_north, &
	1063	kind='w')
	1064
	1065	CALL init_averaging_grid(south_averaged_scalar_profile, cosmo_grid, &
	1066	x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0, &
	1067	lonmin = lam_west, lonmax = lam_east, &
	1068	latmin = phi_centre - averaging_angle, &
	1069	latmax = phi_centre, &
	1070	kind='scalar')
	1071
	1072	CALL init_averaging_grid(north_averaged_scalar_profile, cosmo_grid, &
	1073	x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0, &
	1074	lonmin = lam_west, lonmax = lam_east, &
	1075	latmin = phi_centre, &
	1076	latmax = phi_centre + averaging_angle, &
	1077	kind='scalar')
	1078
	1079	CALL init_averaging_grid(west_averaged_scalar_profile, cosmo_grid, &
	1080	x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0, &
	1081	lonmin = lam_centre - averaging_angle, &
	1082	lonmax = lam_centre, &
	1083	latmin = phi_south, latmax = phi_north, &
	1084	kind='scalar')
	1085
	1086	CALL init_averaging_grid(east_averaged_scalar_profile, cosmo_grid, &
	1087	x = 0.5_dp * lx, y = 0.5_dp * ly, z = z, z0 = z0, &
	1088	lonmin = lam_centre, &
	1089	lonmax = lam_centre + averaging_angle, &
	1090	latmin = phi_south, latmax = phi_north, &
	1091	kind='scalar')
	1092
[2696]	1093	!
	1094	!------------------------------------------------------------------------------
	1095	! Section 4: Precompute neighbours and weights for interpolation
	1096	!------------------------------------------------------------------------------
	1097	interp_mode = 's'
[3182]	1098	CALL setup_interpolation(cosmo_grid, palm_grid, palm_intermediate, interp_mode, ic_mode=cfg % ic_mode)
[2696]	1099	IF (boundary_variables_required) THEN
	1100	CALL setup_interpolation(cosmo_grid, scalars_east_grid, scalars_east_intermediate, interp_mode)
	1101	CALL setup_interpolation(cosmo_grid, scalars_west_grid, scalars_west_intermediate, interp_mode)
	1102	CALL setup_interpolation(cosmo_grid, scalars_north_grid, scalars_north_intermediate, interp_mode)
	1103	CALL setup_interpolation(cosmo_grid, scalars_south_grid, scalars_south_intermediate, interp_mode)
	1104	CALL setup_interpolation(cosmo_grid, scalars_top_grid, scalars_top_intermediate, interp_mode)
	1105	END IF
	1106
	1107	interp_mode = 'u'
[3182]	1108	CALL setup_interpolation(cosmo_grid, u_initial_grid, u_initial_intermediate, interp_mode, ic_mode=cfg % ic_mode)
[2696]	1109	IF (boundary_variables_required) THEN
	1110	CALL setup_interpolation(cosmo_grid, u_east_grid, u_east_intermediate, interp_mode)
	1111	CALL setup_interpolation(cosmo_grid, u_west_grid, u_west_intermediate, interp_mode)
	1112	CALL setup_interpolation(cosmo_grid, u_north_grid, u_north_intermediate, interp_mode)
	1113	CALL setup_interpolation(cosmo_grid, u_south_grid, u_south_intermediate, interp_mode)
	1114	CALL setup_interpolation(cosmo_grid, u_top_grid, u_top_intermediate, interp_mode)
	1115	END IF
	1116
	1117	interp_mode = 'v'
[3182]	1118	CALL setup_interpolation(cosmo_grid, v_initial_grid, v_initial_intermediate, interp_mode, ic_mode=cfg % ic_mode)
[2696]	1119	IF (boundary_variables_required) THEN
	1120	CALL setup_interpolation(cosmo_grid, v_east_grid, v_east_intermediate, interp_mode)
	1121	CALL setup_interpolation(cosmo_grid, v_west_grid, v_west_intermediate, interp_mode)
	1122	CALL setup_interpolation(cosmo_grid, v_north_grid, v_north_intermediate, interp_mode)
	1123	CALL setup_interpolation(cosmo_grid, v_south_grid, v_south_intermediate, interp_mode)
	1124	CALL setup_interpolation(cosmo_grid, v_top_grid, v_top_intermediate, interp_mode)
	1125	END IF
	1126
	1127	interp_mode = 'w'
[3182]	1128	CALL setup_interpolation(cosmo_grid, w_initial_grid, w_initial_intermediate, interp_mode, ic_mode=cfg % ic_mode)
[2696]	1129	IF (boundary_variables_required) THEN
	1130	CALL setup_interpolation(cosmo_grid, w_east_grid, w_east_intermediate, interp_mode)
	1131	CALL setup_interpolation(cosmo_grid, w_west_grid, w_west_intermediate, interp_mode)
	1132	CALL setup_interpolation(cosmo_grid, w_north_grid, w_north_intermediate, interp_mode)
	1133	CALL setup_interpolation(cosmo_grid, w_south_grid, w_south_intermediate, interp_mode)
	1134	CALL setup_interpolation(cosmo_grid, w_top_grid, w_top_intermediate, interp_mode)
	1135	END IF
	1136
[3182]	1137	IF (TRIM(cfg % ic_mode) == 'profile') THEN
[3395]	1138	!TODO: remove this conditional if not needed.
[2696]	1139	END IF
	1140
	1141
	1142	END SUBROUTINE setup_grids
	1143
	1144
[3557]	1145	!------------------------------------------------------------------------------!
	1146	! Description:
	1147	! ------------
	1148	!> Driver for computing neighbour indices and weights for horizontal and
	1149	!> vertical interpolation.
	1150	!------------------------------------------------------------------------------!
[3395]	1151	SUBROUTINE setup_interpolation(cosmo_grid, grid, intermediate_grid, kind, ic_mode)
[2696]	1152
	1153	TYPE(grid_definition), INTENT(IN), TARGET :: cosmo_grid
[3395]	1154	TYPE(grid_definition), INTENT(INOUT), TARGET :: grid, intermediate_grid
[2696]	1155	CHARACTER, INTENT(IN) :: kind
[3182]	1156	CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: ic_mode
[2696]	1157
[3395]	1158	REAL(dp), DIMENSION(:), POINTER :: cosmo_lat, cosmo_lon
	1159	REAL(dp), DIMENSION(:,:,:), POINTER :: cosmo_h
[2696]	1160
[3395]	1161	LOGICAL :: setup_volumetric
[2696]	1162
	1163	!------------------------------------------------------------------------------
	1164	! Section 1: Horizontal interpolation
	1165	!------------------------------------------------------------------------------
[3557]	1166	!
	1167	!-- Select horizontal coordinates according to kind of points (s/w, u, v)
[2696]	1168	SELECT CASE(kind)
	1169
[3557]	1170	!
	1171	!-- scalars
	1172	CASE('s')
[2696]	1173
[3395]	1174	cosmo_lat => cosmo_grid % lat
	1175	cosmo_lon => cosmo_grid % lon
	1176	cosmo_h => cosmo_grid % hfl
[3557]	1177	!
	1178	!-- vertical velocity
	1179	CASE('w')
[2696]	1180
[3395]	1181	cosmo_lat => cosmo_grid % lat
	1182	cosmo_lon => cosmo_grid % lon
	1183	cosmo_h => cosmo_grid % hhl
[3557]	1184	!
	1185	!-- x velocity
	1186	CASE('u')
[2696]	1187
[3395]	1188	cosmo_lat => cosmo_grid % lat
	1189	cosmo_lon => cosmo_grid % lonu
	1190	cosmo_h => cosmo_grid % hfl
[2696]	1191
[3557]	1192	!
	1193	!-- y velocity
	1194	CASE('v')
[2696]	1195
[3395]	1196	cosmo_lat => cosmo_grid % latv
	1197	cosmo_lon => cosmo_grid % lon
	1198	cosmo_h => cosmo_grid % hfl
[2696]	1199
	1200	CASE DEFAULT
	1201
[3182]	1202	message = "Interpolation quantity '" // kind // "' is not supported."
[3615]	1203	CALL inifor_abort('setup_interpolation', message)
[2696]	1204
	1205	END SELECT
	1206
[3395]	1207	CALL find_horizontal_neighbours(cosmo_lat, cosmo_lon, &
	1208	intermediate_grid % clat, intermediate_grid % clon, &
	1209	intermediate_grid % ii, intermediate_grid % jj)
[2696]	1210
[3395]	1211	CALL compute_horizontal_interp_weights(cosmo_lat, cosmo_lon, &
	1212	intermediate_grid % clat, intermediate_grid % clon, &
	1213	intermediate_grid % ii, intermediate_grid % jj, &
	1214	intermediate_grid % w_horiz)
[2696]	1215
	1216	!------------------------------------------------------------------------------
	1217	! Section 2: Vertical interpolation
	1218	!------------------------------------------------------------------------------
	1219
[3557]	1220	!
	1221	!-- If profile initialization is chosen, we--somewhat counterintuitively--
	1222	!-- don't need to compute vertical interpolation weights. At least, we
	1223	!-- don't need them on the intermediate grid, which fills the entire PALM
	1224	!-- domain volume. Instead we need vertical weights for the intermediate
	1225	!-- profile grids, which get computed in setup_averaging().
[3613]	1226	setup_volumetric = .TRUE.
[3395]	1227	IF (PRESENT(ic_mode)) THEN
	1228	IF (TRIM(ic_mode) == 'profile') setup_volumetric = .FALSE.
	1229	END IF
[2696]	1230
[3395]	1231	IF (setup_volumetric) THEN
	1232	ALLOCATE( intermediate_grid % h(0:intermediate_grid % nx, &
	1233	0:intermediate_grid % ny, &
	1234	0:intermediate_grid % nz) )
	1235	intermediate_grid % h(:,:,:) = - EARTH_RADIUS
	1236
[3557]	1237	!
	1238	!-- For w points, use hhl, for scalars use hfl
	1239	!-- compute the full heights for the intermediate grids
[3395]	1240	CALL interpolate_2d(cosmo_h, intermediate_grid % h, intermediate_grid)
	1241	CALL find_vertical_neighbours_and_weights_interp(grid, intermediate_grid)
[2696]	1242	END IF
	1243
	1244	END SUBROUTINE setup_interpolation
	1245
	1246	!------------------------------------------------------------------------------!
	1247	! Description:
	1248	! ------------
	1249	!> Initializes grid_definition-type variables.
	1250	!>
	1251	!> Input parameters:
	1252	!> -----------------
[3182]	1253	!> kind : Grid kind, distinguishes between PALM-4U and COSMO-DE grids
	1254	!> as well as grids covering the boundary surfaces. Valid kinds are:
[2696]	1255	!> - 'palm'
	1256	!> - 'cosmo-de'
	1257	!> - 'eastwest-scalar'
	1258	!>
	1259	!> <xyx>min, <xyz>max : Domain minima and maxima in x, y, and z direction. Note
	1260	!> that these values do not necessarily translate to the outmost coordinates
	1261	!> of the generated grid but rather refer to the extent of the underlying
	1262	!> PALM-4U computational domain (i.e. the outer cell faces). The coordinates
	1263	!> of the generated grid will be inferred from this information taking into
[3182]	1264	!> account the initialization mode ic_mode. For example, the coordinates of a
[2696]	1265	!> boundary grid initialized using mode 'eastwest-scalar' will be located in
	1266	!> planes one half grid point outwards of xmin and xmax.
	1267	!>
	1268	!> z0 : Elevation of the PALM-4U domain above sea level [m]
	1269	!>
	1270	!> n<xyz> : Number of grod points in x, y, and z direction
	1271	!>
	1272	!> Output parameters:
	1273	!> ------------------
	1274	!> grid : Grid variable to be initialized.
	1275	!------------------------------------------------------------------------------!
[3182]	1276	SUBROUTINE init_grid_definition(kind, xmin, xmax, ymin, ymax, &
	1277	x0, y0, z0, nx, ny, nz, z, zw, grid, ic_mode)
[2696]	1278	CHARACTER(LEN=*), INTENT(IN) :: kind
[3182]	1279	CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: ic_mode
[2696]	1280	INTEGER, INTENT(IN) :: nx, ny, nz
[3182]	1281	REAL(dp), INTENT(IN) :: xmin, xmax, ymin, ymax
[2696]	1282	REAL(dp), INTENT(IN) :: x0, y0, z0
[3182]	1283	REAL(dp), INTENT(IN), TARGET, OPTIONAL :: z(:)
	1284	REAL(dp), INTENT(IN), TARGET, OPTIONAL :: zw(:)
[2696]	1285	TYPE(grid_definition), INTENT(INOUT) :: grid
	1286
	1287	grid % nx = nx
	1288	grid % ny = ny
	1289	grid % nz = nz
	1290
	1291	grid % lx = xmax - xmin
	1292	grid % ly = ymax - ymin
	1293
	1294	grid % x0 = x0
	1295	grid % y0 = y0
	1296	grid % z0 = z0
	1297
	1298	SELECT CASE( TRIM (kind) )
	1299
	1300	CASE('boundary')
[3182]	1301
	1302	IF (.NOT.PRESENT(z)) THEN
	1303	message = "z has not been passed but is required for 'boundary' grids"
[3615]	1304	CALL inifor_abort('init_grid_definition', message)
[2696]	1305	END IF
	1306
	1307	ALLOCATE( grid % x(0:nx) )
	1308	CALL linspace(xmin, xmax, grid % x)
	1309
	1310	ALLOCATE( grid % y(0:ny) )
	1311	CALL linspace(ymin, ymax, grid % y)
	1312
[3182]	1313	grid % z => z
[2696]	1314
[3557]	1315	!
	1316	!-- Allocate neighbour indices and weights
[3182]	1317	IF (TRIM(ic_mode) .NE. 'profile') THEN
	1318	ALLOCATE( grid % kk(0:nx, 0:ny, 1:nz, 2) )
[2696]	1319	grid % kk(:,:,:,:) = -1
	1320
[3182]	1321	ALLOCATE( grid % w_verti(0:nx, 0:ny, 1:nz, 2) )
[2696]	1322	grid % w_verti(:,:,:,:) = 0.0_dp
	1323	END IF
	1324
	1325	CASE('boundary intermediate')
	1326
	1327	ALLOCATE( grid % x(0:nx) )
	1328	CALL linspace(xmin, xmax, grid % x)
	1329
	1330	ALLOCATE( grid % y(0:ny) )
	1331	CALL linspace(ymin, ymax, grid % y)
	1332
[3182]	1333	ALLOCATE( grid % clon(0:nx, 0:ny), grid % clat(0:nx, 0:ny) )
[2696]	1334
	1335	CALL rotate_to_cosmo( &
[3182]	1336	phir = project( grid % y, y0, EARTH_RADIUS ) , & ! = plate-carree latitude
	1337	lamr = project( grid % x, x0, EARTH_RADIUS ) , & ! = plate-carree longitude
[2696]	1338	phip = phi_cn, lamp = lambda_cn, &
	1339	phi = grid % clat, &
	1340	lam = grid % clon, &
	1341	gam = gam &
	1342	)
	1343
[3557]	1344	!
	1345	!-- Allocate neighbour indices and weights
[2696]	1346	ALLOCATE( grid % ii(0:nx, 0:ny, 4), &
	1347	grid % jj(0:nx, 0:ny, 4) )
	1348	grid % ii(:,:,:) = -1
	1349	grid % jj(:,:,:) = -1
	1350
	1351	ALLOCATE( grid % w_horiz(0:nx, 0:ny, 4) )
	1352	grid % w_horiz(:,:,:) = 0.0_dp
	1353
[3557]	1354	!
	1355	!-- This mode initializes a Cartesian PALM-4U grid and adds the
	1356	!-- corresponding latitudes and longitudes of the rotated pole grid.
[2696]	1357	CASE('palm')
[3182]	1358
	1359	IF (.NOT.PRESENT(z)) THEN
	1360	message = "z has not been passed but is required for 'palm' grids"
[3615]	1361	CALL inifor_abort('init_grid_definition', message)
[3182]	1362	END IF
	1363
	1364	IF (.NOT.PRESENT(zw)) THEN
	1365	message = "zw has not been passed but is required for 'palm' grids"
[3615]	1366	CALL inifor_abort('init_grid_definition', message)
[3182]	1367	END IF
	1368
[2696]	1369	grid % name(1) = 'x and lon'
	1370	grid % name(2) = 'y and lat'
	1371	grid % name(3) = 'z'
	1372
[3557]	1373	!
	1374	!-- TODO: Remove use of global dx, dy, dz variables. Consider
	1375	!-- TODO: associating global x,y, and z arrays.
[3182]	1376	ALLOCATE( grid % x(0:nx), grid % y(0:ny) )
	1377	ALLOCATE( grid % xu(1:nx), grid % yv(1:ny) )
	1378	CALL linspace(xmin + 0.5_dp* dx, xmax - 0.5_dp* dx, grid % x)
	1379	CALL linspace(ymin + 0.5_dp* dy, ymax - 0.5_dp* dy, grid % y)
	1380	grid % z => z
	1381	CALL linspace(xmin + dx, xmax - dx, grid % xu)
	1382	CALL linspace(ymin + dy, ymax - dy, grid % yv)
	1383	grid % zw => zw
[2696]	1384
	1385	grid % depths => depths
	1386
[3557]	1387	!
	1388	!-- Allocate neighbour indices and weights
[3182]	1389	IF (TRIM(ic_mode) .NE. 'profile') THEN
	1390	ALLOCATE( grid % kk(0:nx, 0:ny, 1:nz, 2) )
[2696]	1391	grid % kk(:,:,:,:) = -1
	1392
[3182]	1393	ALLOCATE( grid % w_verti(0:nx, 0:ny, 1:nz, 2) )
[2696]	1394	grid % w_verti(:,:,:,:) = 0.0_dp
	1395	END IF
	1396
	1397	CASE('palm intermediate')
[3182]	1398
[2696]	1399	grid % name(1) = 'x and lon'
	1400	grid % name(2) = 'y and lat'
[3182]	1401	grid % name(3) = 'interpolated hhl or hfl'
[2696]	1402
[3557]	1403	!
	1404	!-- TODO: Remove use of global dx, dy, dz variables. Consider
	1405	!-- TODO: associating global x,y, and z arrays.
[3182]	1406	ALLOCATE( grid % x(0:nx), grid % y(0:ny) )
	1407	ALLOCATE( grid % xu(1:nx), grid % yv(1:ny) )
	1408	CALL linspace(xmin + 0.5_dpdx, xmax - 0.5_dpdx, grid % x)
	1409	CALL linspace(ymin + 0.5_dpdy, ymax - 0.5_dpdy, grid % y)
	1410	CALL linspace(xmin + dx, xmax - dx, grid % xu)
	1411	CALL linspace(ymin + dy, ymax - dy, grid % yv)
[2696]	1412
	1413	grid % depths => depths
	1414
[3557]	1415	!
	1416	!-- Allocate rotated-pole coordinates, clon is for (c)osmo-de (lon)gitude
[2696]	1417	ALLOCATE( grid % clon(0:nx, 0:ny), grid % clat(0:nx, 0:ny) )
	1418	ALLOCATE( grid % clonu(1:nx, 0:ny), grid % clatu(1:nx, 0:ny) )
	1419	ALLOCATE( grid % clonv(0:nx, 1:ny), grid % clatv(0:nx, 1:ny) )
	1420
[3557]	1421	!
	1422	!-- Compute rotated-pole coordinates of...
	1423	!-- ... PALM-4U centres
[2696]	1424	CALL rotate_to_cosmo( &
	1425	phir = project( grid % y, y0, EARTH_RADIUS ) , & ! = plate-carree latitude
	1426	lamr = project( grid % x, x0, EARTH_RADIUS ) , & ! = plate-carree longitude
	1427	phip = phi_cn, lamp = lambda_cn, &
	1428	phi = grid % clat, &
	1429	lam = grid % clon, &
	1430	gam = gam &
	1431	)
	1432
[3557]	1433	!
	1434	!-- ... PALM-4U u winds
[2696]	1435	CALL rotate_to_cosmo( &
	1436	phir = project( grid % y, y0, EARTH_RADIUS ), & ! = plate-carree latitude
	1437	lamr = project( grid % xu, x0, EARTH_RADIUS ), & ! = plate-carree longitude
	1438	phip = phi_cn, lamp = lambda_cn, &
	1439	phi = grid % clatu, &
	1440	lam = grid % clonu, &
	1441	gam = gam &
	1442	)
	1443
[3557]	1444	!
	1445	!-- ... PALM-4U v winds
[2696]	1446	CALL rotate_to_cosmo( &
	1447	phir = project( grid % yv, y0, EARTH_RADIUS ), & ! = plate-carree latitude
	1448	lamr = project( grid % x, x0, EARTH_RADIUS ), & ! = plate-carree longitude
	1449	phip = phi_cn, lamp = lambda_cn, &
	1450	phi = grid % clatv, &
	1451	lam = grid % clonv, &
	1452	gam = gam &
	1453	)
	1454
[3557]	1455	!
	1456	!-- Allocate neighbour indices and weights
[2696]	1457	ALLOCATE( grid % ii(0:nx, 0:ny, 4), &
	1458	grid % jj(0:nx, 0:ny, 4) )
	1459	grid % ii(:,:,:) = -1
	1460	grid % jj(:,:,:) = -1
	1461
	1462	ALLOCATE( grid % w_horiz(0:nx, 0:ny, 4) )
	1463	grid % w_horiz(:,:,:) = 0.0_dp
	1464
	1465	CASE('cosmo-de')
	1466	grid % name(1) = 'rlon' ! of COMSO-DE cell centres (scalars)
	1467	grid % name(2) = 'rlat' ! of COMSO-DE cell centres (scalars)
	1468	grid % name(3) = 'height'
	1469
[3182]	1470	ALLOCATE( grid % lon(0:nx), grid % lat(0:ny) )
	1471	ALLOCATE( grid % lonu(0:nx), grid % latv(0:ny) )
[2696]	1472
	1473	CALL linspace(xmin, xmax, grid % lon)
	1474	CALL linspace(ymin, ymax, grid % lat)
[3182]	1475	grid % lonu(:) = grid % lon + 0.5_dp * (grid % lx / grid % nx)
	1476	grid % latv(:) = grid % lat + 0.5_dp * (grid % ly / grid % ny)
[2696]	1477
[3557]	1478	!
	1479	!-- Point to heights of half levels (hhl) and compute heights of full
	1480	!-- levels (hfl) as arithmetic averages
[2696]	1481	grid % hhl => hhl
	1482	grid % hfl => hfl
	1483	grid % depths => depths
	1484
	1485	CASE DEFAULT
	1486	message = "Grid kind '" // TRIM(kind) // "' is not recognized."
[3615]	1487	CALL inifor_abort('init_grid_definition', message)
[2696]	1488
	1489	END SELECT
	1490
	1491	END SUBROUTINE init_grid_definition
	1492
	1493
[3182]	1494	!------------------------------------------------------------------------------!
	1495	! Description:
	1496	! ------------
[3613]	1497	!> Initializes averaging grids
[3395]	1498	!>
	1499	!> Input parameters:
	1500	!> -----------------
	1501	!>
	1502	!> cosmo_grid : grid_definition-type varieble initialized as 'cosmo-de' grid
	1503	!> providing longitudes and latitudes of the parent grid
	1504	!>
	1505	!> x, y : location of the profile to be averaged in the PALM system [m]
	1506	!>
	1507	!> z : vertical grid coordinates of the profile in the PALM system [m]
	1508	!>
	1509	!> <lat/lon><min/max>: boundaries of the averaging region in the parent system,
	1510	!> i.e. the COSMO-DE rotated-pole system. [rad]
	1511	!>
	1512	!> kind : kind of quantity to be averaged using this averaging grid.
	1513	!> Destinguishes COSMO-DE scalar and w-velocity levels. Note that finding the
	1514	!> parent/COSMO columns for the region in get_cosmo_averaging_region() is
	1515	!> independent of 'kind' b/ecause INIFOR uses column-centred u and v velocity
	1516	!> components, which are computed in the preprocessing step.
	1517	!>
	1518	!> Output parameters:
	1519	!> ------------------
[3613]	1520	!> avg_grid : averagin grid to be initialized
[3395]	1521	!------------------------------------------------------------------------------!
	1522	SUBROUTINE init_averaging_grid(avg_grid, cosmo_grid, x, y, z, z0, &
	1523	lonmin, lonmax, latmin, latmax, kind)
	1524
	1525	TYPE(grid_definition), INTENT(INOUT) :: avg_grid
	1526	TYPE(grid_definition), INTENT(IN) :: cosmo_grid
	1527	REAL(dp), INTENT(IN) :: x, y, z0
	1528	REAL(dp), INTENT(IN), TARGET :: z(:)
[3613]	1529	REAL(dp), INTENT(IN) :: lonmin !< lower longitude bound of the averaging grid region [COSMO rotated-pole rad]
	1530	REAL(dp), INTENT(IN) :: lonmax !< upper longitude bound of the averaging grid region [COSMO rotated-pole rad]
	1531	REAL(dp), INTENT(IN) :: latmin !< lower latitude bound of the averaging grid region [COSMO rotated-pole rad]
	1532	REAL(dp), INTENT(IN) :: latmax !< lower latitude bound of the averaging grid region [COSMO rotated-pole rad]
[3395]	1533
	1534	CHARACTER(LEN=*), INTENT(IN) :: kind
	1535
	1536	LOGICAL :: level_based_averaging
	1537
	1538	ALLOCATE( avg_grid % x(1) )
	1539	ALLOCATE( avg_grid % y(1) )
	1540	avg_grid % x(1) = x
	1541	avg_grid % y(1) = y
	1542	avg_grid % z => z
	1543	avg_grid % z0 = z0
	1544
	1545	avg_grid % nz = SIZE(z, 1)
	1546
	1547	ALLOCATE( avg_grid % lon(2) )
	1548	ALLOCATE( avg_grid % lat(2) )
	1549	avg_grid % lon(1:2) = (/lonmin, lonmax/)
	1550	avg_grid % lat(1:2) = (/latmin, latmax/)
	1551
	1552	avg_grid % kind = TRIM(kind)
	1553
[3557]	1554	!
	1555	!-- Find and store COSMO columns that fall into the coordinate range
	1556	!-- given by avg_grid % clon, % clat
[3395]	1557	CALL get_cosmo_averaging_region(avg_grid, cosmo_grid)
	1558
	1559	ALLOCATE (avg_grid % kkk(avg_grid % n_columns, avg_grid % nz, 2) )
	1560	ALLOCATE (avg_grid % w(avg_grid % n_columns, avg_grid % nz, 2) )
[3557]	1561	!
	1562	!-- Compute average COSMO levels in the averaging region
[3395]	1563	SELECT CASE(avg_grid % kind)
	1564
	1565	CASE('scalar', 'u', 'v')
	1566	avg_grid % cosmo_h => cosmo_grid % hfl
	1567
	1568	CASE('w')
	1569	avg_grid % cosmo_h => cosmo_grid % hhl
	1570
	1571	CASE DEFAULT
	1572	message = "Averaging grid kind '" // TRIM(avg_grid % kind) // &
	1573	"' is not supported. Use 'scalar', 'u', or 'v'."
[3615]	1574	CALL inifor_abort('get_cosmo_averaging_region', message)
[3395]	1575
	1576	END SELECT
	1577
[3557]	1578	!
	1579	!-- For level-besed averaging, compute average heights
[3613]	1580	level_based_averaging = ( TRIM(cfg % averaging_mode) == 'level' )
[3395]	1581	IF (level_based_averaging) THEN
	1582	ALLOCATE(avg_grid % h(1,1,SIZE(avg_grid % cosmo_h, 3)) )
	1583
[3613]	1584	CALL average_2d(avg_grid % cosmo_h, avg_grid % h(1,1,:), &
[3395]	1585	avg_grid % iii, avg_grid % jjj)
	1586	END IF
	1587
[3557]	1588	!
	1589	!-- Compute vertical weights and neighbours
[3613]	1590	CALL find_vertical_neighbours_and_weights_average( &
	1591	avg_grid, level_based_averaging &
	1592	)
[3395]	1593
	1594	END SUBROUTINE init_averaging_grid
	1595
	1596
	1597	SUBROUTINE get_cosmo_averaging_region(avg_grid, cosmo_grid)
	1598	TYPE(grid_definition), INTENT(INOUT) :: avg_grid
	1599	TYPE(grid_definition), TARGET, INTENT(IN) :: cosmo_grid
	1600
	1601	REAL(dp), DIMENSION(:), POINTER :: cosmo_lon, cosmo_lat
	1602	REAL(dp) :: dlon, dlat
	1603
	1604	INTEGER :: i, j, imin, imax, jmin, jmax, l, nx, ny
	1605
	1606
	1607	SELECT CASE( TRIM(avg_grid % kind) )
	1608
	1609	CASE('scalar', 'w')
	1610	cosmo_lon => cosmo_grid % lon
	1611	cosmo_lat => cosmo_grid % lat
	1612
	1613	CASE('u')
	1614	cosmo_lon => cosmo_grid % lonu
	1615	cosmo_lat => cosmo_grid % lat
	1616
	1617	CASE('v')
	1618	cosmo_lon => cosmo_grid % lon
	1619	cosmo_lat => cosmo_grid % latv
	1620
	1621	CASE DEFAULT
	1622	message = "Averaging grid kind '" // TRIM(avg_grid % kind) // &
	1623	"' is not supported. Use 'scalar', 'u', or 'v'."
[3615]	1624	CALL inifor_abort('get_cosmo_averaging_region', message)
[3395]	1625
	1626	END SELECT
	1627
[3557]	1628	!
	1629	!-- FIXME: make dlon, dlat parameters of the grid_defintion type
[3395]	1630	dlon = cosmo_lon(1) - cosmo_lon(0)
	1631	dlat = cosmo_lat(1) - cosmo_lat(0)
	1632
[3678]	1633	imin = FLOOR ( (avg_grid % lon(1) - cosmo_lon(0)) / dlon )
	1634	imax = CEILING( (avg_grid % lon(2) - cosmo_lon(0)) / dlon )
[3395]	1635
[3678]	1636	jmin = FLOOR ( (avg_grid % lat(1) - cosmo_lat(0)) / dlat )
	1637	jmax = CEILING( (avg_grid % lat(2) - cosmo_lat(0)) / dlat )
[3395]	1638
[3537]	1639	message = "Averaging '" // TRIM(avg_grid % kind) // "' over "// &
	1640	TRIM(str(imin)) // " <= i <= " // TRIM(str(imax)) // &
[3395]	1641	" and " // &
[3537]	1642	TRIM(str(jmin)) // " <= j <= " // TRIM(str(jmax))
[3395]	1643	CALL report( 'get_cosmo_averaging_region', message )
	1644
	1645	nx = imax - imin + 1
	1646	ny = jmax - jmin + 1
	1647	avg_grid % n_columns = nx * ny
	1648
[3537]	1649	ALLOCATE( avg_grid % iii(avg_grid % n_columns), &
[3395]	1650	avg_grid % jjj(avg_grid % n_columns) )
	1651
	1652	l = 0
	1653	DO j = jmin, jmax
	1654	DO i = imin, imax
	1655	l = l + 1
	1656	avg_grid % iii(l) = i
	1657	avg_grid % jjj(l) = j
	1658	END DO
	1659	END DO
	1660
	1661	END SUBROUTINE get_cosmo_averaging_region
	1662
	1663
	1664	!------------------------------------------------------------------------------!
	1665	! Description:
	1666	! ------------
[3182]	1667	!> PALM's stretched vertical grid generator. Forked from PALM revision 3139, see
	1668	!> https://palm.muk.uni-hannover.de/trac/browser/palm/trunk/SOURCE/init_grid.f90?rev=3139
	1669	!>
	1670	!> This routine computes the levels of scalar points. The levels of the velocity
	1671	!> points are then obtained as the midpoints inbetween using the INIFOR routine
	1672	!> 'modpoints'.
	1673	!------------------------------------------------------------------------------!
	1674	SUBROUTINE stretched_z(z, dz, dz_max, dz_stretch_factor, dz_stretch_level, &
	1675	dz_stretch_level_start, dz_stretch_level_end, &
	1676	dz_stretch_factor_array)
	1677
	1678	REAL(dp), DIMENSION(:), INTENT(INOUT) :: z, dz, dz_stretch_factor_array
	1679	REAL(dp), DIMENSION(:), INTENT(INOUT) :: dz_stretch_level_start, dz_stretch_level_end
	1680	REAL(dp), INTENT(IN) :: dz_max, dz_stretch_factor, dz_stretch_level
	1681
	1682	INTEGER :: number_stretch_level_start !< number of user-specified start levels for stretching
	1683	INTEGER :: number_stretch_level_end !< number of user-specified end levels for stretching
	1684
	1685	REAL(dp), DIMENSION(:), ALLOCATABLE :: min_dz_stretch_level_end
	1686	REAL(dp) :: dz_level_end, dz_stretched
	1687
[3557]	1688	INTEGER :: dz_stretch_level_end_index(9) !< vertical grid level index until which the vertical grid spacing is stretched
	1689	INTEGER :: dz_stretch_level_start_index(9) !< vertical grid level index above which the vertical grid spacing is stretched
[3182]	1690	INTEGER :: dz_stretch_level_index = 0
	1691	INTEGER :: k, n, number_dz
[3557]	1692
[3182]	1693	!
	1694	!-- Compute height of u-levels from constant grid length and dz stretch factors
	1695	IF ( dz(1) == -1.0_dp ) THEN
	1696	message = 'missing dz'
[3615]	1697	CALL inifor_abort( 'stretched_z', message)
[3182]	1698	ELSEIF ( dz(1) <= 0.0_dp ) THEN
	1699	WRITE( message, * ) 'dz=', dz(1),' <= 0.0'
[3615]	1700	CALL inifor_abort( 'stretched_z', message)
[3182]	1701	ENDIF
	1702
	1703	!
	1704	!-- Initialize dz_stretch_level_start with the value of dz_stretch_level
	1705	!-- if it was set by the user
	1706	IF ( dz_stretch_level /= -9999999.9_dp ) THEN
	1707	dz_stretch_level_start(1) = dz_stretch_level
	1708	ENDIF
	1709
	1710	!
	1711	!-- Determine number of dz values and stretching levels specified by the
	1712	!-- user to allow right controlling of the stretching mechanism and to
	1713	!-- perform error checks. The additional requirement that dz /= dz_max
	1714	!-- for counting number of user-specified dz values is necessary. Otherwise
	1715	!-- restarts would abort if the old stretching mechanism with dz_stretch_level
	1716	!-- is used (Attention: The user is not allowed to specify a dz value equal
	1717	!-- to the default of dz_max = 999.0).
	1718	number_dz = COUNT( dz /= -1.0_dp .AND. dz /= dz_max )
	1719	number_stretch_level_start = COUNT( dz_stretch_level_start /= &
	1720	-9999999.9_dp )
	1721	number_stretch_level_end = COUNT( dz_stretch_level_end /= &
	1722	9999999.9_dp )
	1723
	1724	!
	1725	!-- The number of specified end levels +1 has to be the same than the number
	1726	!-- of specified dz values
	1727	IF ( number_dz /= number_stretch_level_end + 1 ) THEN
	1728	WRITE( message, * ) 'The number of values for dz = ', &
	1729	number_dz, 'has to be the same than ', &
	1730	'the number of values for ', &
	1731	'dz_stretch_level_end + 1 = ', &
	1732	number_stretch_level_end+1
[3615]	1733	CALL inifor_abort( 'stretched_z', message)
[3182]	1734	ENDIF
	1735
	1736	!
	1737	!-- The number of specified start levels has to be the same or one less than
	1738	!-- the number of specified dz values
	1739	IF ( number_dz /= number_stretch_level_start + 1 .AND. &
	1740	number_dz /= number_stretch_level_start ) THEN
	1741	WRITE( message, * ) 'The number of values for dz = ', &
	1742	number_dz, 'has to be the same or one ', &
	1743	'more than& the number of values for ', &
	1744	'dz_stretch_level_start = ', &
	1745	number_stretch_level_start
[3615]	1746	CALL inifor_abort( 'stretched_z', message)
[3182]	1747	ENDIF
	1748
	1749	!-- The number of specified start levels has to be the same or one more than
	1750	!-- the number of specified end levels
	1751	IF ( number_stretch_level_start /= number_stretch_level_end + 1 .AND. &
	1752	number_stretch_level_start /= number_stretch_level_end ) THEN
	1753	WRITE( message, * ) 'The number of values for ', &
	1754	'dz_stretch_level_start = ', &
	1755	dz_stretch_level_start, 'has to be the ',&
	1756	'same or one more than& the number of ', &
	1757	'values for dz_stretch_level_end = ', &
	1758	number_stretch_level_end
[3615]	1759	CALL inifor_abort( 'stretched_z', message)
[3182]	1760	ENDIF
	1761
	1762	!
	1763	!-- Initialize dz for the free atmosphere with the value of dz_max
	1764	IF ( dz(number_stretch_level_start+1) == -1.0_dp .AND. &
	1765	number_stretch_level_start /= 0 ) THEN
	1766	dz(number_stretch_level_start+1) = dz_max
	1767	ENDIF
	1768
	1769	!
	1770	!-- Initialize the stretching factor if (infinitely) stretching in the free
	1771	!-- atmosphere is desired (dz_stretch_level_end was not specified for the
	1772	!-- free atmosphere)
	1773	IF ( number_stretch_level_start == number_stretch_level_end + 1 ) THEN
	1774	dz_stretch_factor_array(number_stretch_level_start) = &
	1775	dz_stretch_factor
	1776	ENDIF
	1777
	1778	!-- Allocation of arrays for stretching
	1779	ALLOCATE( min_dz_stretch_level_end(number_stretch_level_start) )
	1780
	1781	!
	1782	!-- The stretching region has to be large enough to allow for a smooth
	1783	!-- transition between two different grid spacings
	1784	DO n = 1, number_stretch_level_start
	1785	min_dz_stretch_level_end(n) = dz_stretch_level_start(n) + &
	1786	4 * MAX( dz(n),dz(n+1) )
	1787	ENDDO
	1788
	1789	IF ( ANY( min_dz_stretch_level_end(1:number_stretch_level_start) > &
	1790	dz_stretch_level_end(1:number_stretch_level_start) ) ) THEN
	1791	!IF ( ANY( min_dz_stretch_level_end > &
	1792	! dz_stretch_level_end ) ) THEN
	1793	message = 'Each dz_stretch_level_end has to be larger ' // &
	1794	'than its corresponding value for ' // &
	1795	'dz_stretch_level_start + 4*MAX(dz(n),dz(n+1)) '//&
	1796	'to allow for smooth grid stretching'
[3615]	1797	CALL inifor_abort('stretched_z', message)
[3182]	1798	ENDIF
	1799
	1800	!
	1801	!-- Stretching must not be applied within the prandtl_layer
	1802	!-- (first two grid points). For the default case dz_stretch_level_start
	1803	!-- is negative. Therefore the absolut value is checked here.
	1804	IF ( ANY( ABS( dz_stretch_level_start ) < dz(1) * 1.5_dp ) ) THEN
	1805	WRITE( message, * ) 'Eeach dz_stretch_level_start has to be ',&
	1806	'larger than ', dz(1) * 1.5
[3615]	1807	CALL inifor_abort( 'stretched_z', message)
[3182]	1808	ENDIF
	1809
	1810	!
	1811	!-- The stretching has to start and end on a grid level. Therefore
	1812	!-- user-specified values have to ''interpolate'' to the next lowest level
	1813	IF ( number_stretch_level_start /= 0 ) THEN
	1814	dz_stretch_level_start(1) = INT( (dz_stretch_level_start(1) - &
	1815	dz(1)/2.0) / dz(1) ) &
	1816	* dz(1) + dz(1)/2.0
	1817	ENDIF
	1818
	1819	IF ( number_stretch_level_start > 1 ) THEN
	1820	DO n = 2, number_stretch_level_start
	1821	dz_stretch_level_start(n) = INT( dz_stretch_level_start(n) / &
	1822	dz(n) ) * dz(n)
	1823	ENDDO
	1824	ENDIF
	1825
	1826	IF ( number_stretch_level_end /= 0 ) THEN
	1827	DO n = 1, number_stretch_level_end
	1828	dz_stretch_level_end(n) = INT( dz_stretch_level_end(n) / &
	1829	dz(n+1) ) * dz(n+1)
	1830	ENDDO
	1831	ENDIF
	1832
	1833	!
	1834	!-- Determine stretching factor if necessary
	1835	IF ( number_stretch_level_end >= 1 ) THEN
	1836	CALL calculate_stretching_factor( number_stretch_level_end, dz, &
	1837	dz_stretch_factor_array, &
	1838	dz_stretch_level_end, &
	1839	dz_stretch_level_start )
	1840	ENDIF
	1841
	1842	z(1) = dz(1) * 0.5_dp
	1843	!
	1844	dz_stretch_level_index = n
	1845	dz_stretched = dz(1)
	1846	DO k = 2, n
	1847
	1848	IF ( dz_stretch_level <= z(k-1) .AND. dz_stretched < dz_max ) THEN
	1849
	1850	dz_stretched = dz_stretched * dz_stretch_factor
	1851	dz_stretched = MIN( dz_stretched, dz_max )
	1852
	1853	IF ( dz_stretch_level_index == n ) dz_stretch_level_index = k-1
	1854
	1855	ENDIF
	1856
	1857	z(k) = z(k-1) + dz_stretched
	1858
	1859	ENDDO
	1860	!-- Determine u and v height levels considering the possibility of grid
	1861	!-- stretching in several heights.
	1862	n = 1
	1863	dz_stretch_level_start_index(:) = UBOUND(z, 1)
	1864	dz_stretch_level_end_index(:) = UBOUND(z, 1)
	1865	dz_stretched = dz(1)
	1866
	1867	!-- The default value of dz_stretch_level_start is negative, thus the first
	1868	!-- condition is always true. Hence, the second condition is necessary.
	1869	DO k = 2, UBOUND(z, 1)
	1870	IF ( dz_stretch_level_start(n) <= z(k-1) .AND. &
	1871	dz_stretch_level_start(n) /= -9999999.9_dp ) THEN
	1872	dz_stretched = dz_stretched * dz_stretch_factor_array(n)
	1873
	1874	IF ( dz(n) > dz(n+1) ) THEN
	1875	dz_stretched = MAX( dz_stretched, dz(n+1) ) !Restrict dz_stretched to the user-specified (higher) dz
	1876	ELSE
	1877	dz_stretched = MIN( dz_stretched, dz(n+1) ) !Restrict dz_stretched to the user-specified (lower) dz
	1878	ENDIF
	1879
	1880	IF ( dz_stretch_level_start_index(n) == UBOUND(z, 1) ) &
	1881	dz_stretch_level_start_index(n) = k-1
	1882
	1883	ENDIF
	1884
	1885	z(k) = z(k-1) + dz_stretched
	1886
	1887	!
	1888	!-- Make sure that the stretching ends exactly at dz_stretch_level_end
	1889	dz_level_end = ABS( z(k) - dz_stretch_level_end(n) )
	1890
	1891	IF ( dz_level_end < dz(n+1)/3.0 ) THEN
	1892	z(k) = dz_stretch_level_end(n)
	1893	dz_stretched = dz(n+1)
	1894	dz_stretch_level_end_index(n) = k
	1895	n = n + 1
	1896	ENDIF
	1897	ENDDO
	1898
	1899	DEALLOCATE( min_dz_stretch_level_end )
	1900
	1901	END SUBROUTINE stretched_z
	1902
	1903
[3557]	1904	!------------------------------------------------------------------------------!
[3182]	1905	! Description: [PALM subroutine]
	1906	! -----------------------------------------------------------------------------!
	1907	!> Calculation of the stretching factor through an iterative method. Ideas were
	1908	!> taken from the paper "Regional stretched grid generation and its application
	1909	!> to the NCAR RegCM (1999)". Normally, no analytic solution exists because the
	1910	!> system of equations has two variables (r,l) but four requirements
	1911	!> (l=integer, r=[0,88;1,2], Eq(6), Eq(5) starting from index j=1) which
	1912	!> results into an overdetermined system.
	1913	!------------------------------------------------------------------------------!
[3618]	1914	SUBROUTINE calculate_stretching_factor( number_end, dz, &
[3182]	1915	dz_stretch_factor_array, &
	1916	dz_stretch_level_end, &
	1917	dz_stretch_level_start )
	1918
	1919	REAL(dp), DIMENSION(:), INTENT(IN) :: dz
	1920	REAL(dp), DIMENSION(:), INTENT(INOUT) :: dz_stretch_factor_array
	1921	REAL(dp), DIMENSION(:), INTENT(IN) :: dz_stretch_level_end, dz_stretch_level_start
	1922
	1923	INTEGER :: iterations !< number of iterations until stretch_factor_lower/upper_limit is reached
	1924	INTEGER :: l_rounded !< after l_rounded grid levels dz(n) is strechted to dz(n+1) with stretch_factor_2
	1925	INTEGER :: n !< loop variable for stretching
	1926
	1927	INTEGER, INTENT(IN) :: number_end !< number of user-specified end levels for stretching
	1928
	1929	REAL(dp) :: delta_l !< absolute difference between l and l_rounded
	1930	REAL(dp) :: delta_stretch_factor !< absolute difference between stretch_factor_1 and stretch_factor_2
	1931	REAL(dp) :: delta_total_new !< sum of delta_l and delta_stretch_factor for the next iteration (should be as small as possible)
	1932	REAL(dp) :: delta_total_old !< sum of delta_l and delta_stretch_factor for the last iteration
	1933	REAL(dp) :: distance !< distance between dz_stretch_level_start and dz_stretch_level_end (stretching region)
	1934	REAL(dp) :: l !< value that fulfil Eq. (5) in the paper mentioned above together with stretch_factor_1 exactly
	1935	REAL(dp) :: numerator !< numerator of the quotient
	1936	REAL(dp) :: stretch_factor_1 !< stretching factor that fulfil Eq. (5) togehter with l exactly
	1937	REAL(dp) :: stretch_factor_2 !< stretching factor that fulfil Eq. (6) togehter with l_rounded exactly
	1938
	1939	REAL(dp) :: dz_stretch_factor_array_2(9) = 1.08_dp !< Array that contains all stretch_factor_2 that belongs to stretch_factor_1
	1940
	1941	REAL(dp), PARAMETER :: stretch_factor_interval = 1.0E-06 !< interval for sampling possible stretching factors
	1942	REAL(dp), PARAMETER :: stretch_factor_lower_limit = 0.88 !< lowest possible stretching factor
	1943	REAL(dp), PARAMETER :: stretch_factor_upper_limit = 1.12 !< highest possible stretching factor
	1944
	1945
	1946	l = 0
	1947	DO n = 1, number_end
	1948
	1949	iterations = 1
	1950	stretch_factor_1 = 1.0
	1951	stretch_factor_2 = 1.0
	1952	delta_total_old = 1.0
	1953
	1954	IF ( dz(n) > dz(n+1) ) THEN
	1955	DO WHILE ( stretch_factor_1 >= stretch_factor_lower_limit )
	1956
	1957	stretch_factor_1 = 1.0 - iterations * stretch_factor_interval
	1958	distance = ABS( dz_stretch_level_end(n) - &
	1959	dz_stretch_level_start(n) )
	1960	numerator = distance*stretch_factor_1/dz(n) + &
	1961	stretch_factor_1 - distance/dz(n)
	1962
	1963	IF ( numerator > 0.0 ) THEN
	1964	l = LOG( numerator ) / LOG( stretch_factor_1 ) - 1.0
	1965	l_rounded = NINT( l )
	1966	delta_l = ABS( l_rounded - l ) / l
	1967	ENDIF
	1968
	1969	stretch_factor_2 = EXP( LOG( dz(n+1)/dz(n) ) / (l_rounded) )
	1970
	1971	delta_stretch_factor = ABS( stretch_factor_1 - &
	1972	stretch_factor_2 ) / &
	1973	stretch_factor_2
	1974
	1975	delta_total_new = delta_l + delta_stretch_factor
	1976
	1977	!
	1978	!-- stretch_factor_1 is taken to guarantee that the stretching
	1979	!-- procedure ends as close as possible to dz_stretch_level_end.
	1980	!-- stretch_factor_2 would guarantee that the stretched dz(n) is
	1981	!-- equal to dz(n+1) after l_rounded grid levels.
	1982	IF (delta_total_new < delta_total_old) THEN
	1983	dz_stretch_factor_array(n) = stretch_factor_1
	1984	dz_stretch_factor_array_2(n) = stretch_factor_2
	1985	delta_total_old = delta_total_new
	1986	ENDIF
	1987
	1988	iterations = iterations + 1
	1989
	1990	ENDDO
	1991
	1992	ELSEIF ( dz(n) < dz(n+1) ) THEN
	1993	DO WHILE ( stretch_factor_1 <= stretch_factor_upper_limit )
	1994
	1995	stretch_factor_1 = 1.0 + iterations * stretch_factor_interval
	1996	distance = ABS( dz_stretch_level_end(n) - &
	1997	dz_stretch_level_start(n) )
	1998	numerator = distance*stretch_factor_1/dz(n) + &
	1999	stretch_factor_1 - distance/dz(n)
	2000
	2001	l = LOG( numerator ) / LOG( stretch_factor_1 ) - 1.0
	2002	l_rounded = NINT( l )
	2003	delta_l = ABS( l_rounded - l ) / l
	2004
	2005	stretch_factor_2 = EXP( LOG( dz(n+1)/dz(n) ) / (l_rounded) )
	2006
	2007	delta_stretch_factor = ABS( stretch_factor_1 - &
	2008	stretch_factor_2 ) / &
	2009	stretch_factor_2
	2010
	2011	delta_total_new = delta_l + delta_stretch_factor
	2012
	2013	!
[3557]	2014	!-- stretch_factor_1 is taken to guarantee that the stretching
	2015	!-- procedure ends as close as possible to dz_stretch_level_end.
	2016	!-- stretch_factor_2 would guarantee that the stretched dz(n) is
	2017	!-- equal to dz(n+1) after l_rounded grid levels.
[3182]	2018	IF (delta_total_new < delta_total_old) THEN
	2019	dz_stretch_factor_array(n) = stretch_factor_1
	2020	dz_stretch_factor_array_2(n) = stretch_factor_2
	2021	delta_total_old = delta_total_new
	2022	ENDIF
	2023
	2024	iterations = iterations + 1
	2025	ENDDO
	2026
	2027	ELSE
	2028	message = 'Two adjacent values of dz must be different'
[3615]	2029	CALL inifor_abort( 'calculate_stretching_factor', message)
[3182]	2030	ENDIF
	2031
	2032	!
	2033	!-- Check if also the second stretching factor fits into the allowed
	2034	!-- interval. If not, print a warning for the user.
	2035	IF ( dz_stretch_factor_array_2(n) < stretch_factor_lower_limit .OR. &
	2036	dz_stretch_factor_array_2(n) > stretch_factor_upper_limit ) THEN
	2037	WRITE( message, * ) 'stretch_factor_2 = ', &
	2038	dz_stretch_factor_array_2(n), ' which is',&
	2039	' responsible for exactly reaching& dz =',&
	2040	dz(n+1), 'after a specific amount of', &
	2041	' grid levels& exceeds the upper', &
	2042	' limit =', stretch_factor_upper_limit, &
	2043	' &or lower limit = ', &
	2044	stretch_factor_lower_limit
[3615]	2045	CALL inifor_abort( 'calculate_stretching_factor', message )
[3182]	2046
	2047	ENDIF
	2048	ENDDO
	2049
	2050	END SUBROUTINE calculate_stretching_factor
	2051
[3557]	2052	!------------------------------------------------------------------------------!
	2053	! Description:
	2054	! ------------
	2055	!> Computes the midpoint between two neighbouring coordinates in the given
	2056	!> coordinate vector 'z' and stores it in 'zw'.
	2057	!------------------------------------------------------------------------------!
[3182]	2058	SUBROUTINE midpoints(z, zw)
	2059
	2060	REAL(dp), INTENT(IN) :: z(0:)
	2061	REAL(dp), INTENT(OUT) :: zw(1:)
	2062
	2063	INTEGER :: k
	2064
	2065	DO k = 1, UBOUND(zw, 1)
	2066	zw(k) = 0.5_dp * (z(k-1) + z(k))
	2067	END DO
	2068
	2069	END SUBROUTINE midpoints
	2070
[3557]	2071	!------------------------------------------------------------------------------!
	2072	! Description:
	2073	! ------------
	2074	!> Defines INFOR's IO groups.
	2075	!------------------------------------------------------------------------------!
[2696]	2076	SUBROUTINE setup_io_groups()
	2077
	2078	INTEGER :: ngroups
	2079
	2080	ngroups = 16
	2081	ALLOCATE( io_group_list(ngroups) )
[3557]	2082
	2083	!
	2084	!-- soil temp
[2696]	2085	io_group_list(1) = init_io_group( &
	2086	in_files = soil_files, &
	2087	out_vars = output_var_table(1:1), &
	2088	in_var_list = input_var_table(1:1), &
	2089	kind = 'soil-temperature' &
	2090	)
	2091
[3557]	2092	!
	2093	!-- soil water
[2696]	2094	io_group_list(2) = init_io_group( &
	2095	in_files = soil_files, &
	2096	out_vars = output_var_table(2:2), &
	2097	in_var_list = input_var_table(2:2), &
	2098	kind = 'soil-water' &
	2099	)
	2100
[3557]	2101	!
	2102	!-- potential temperature, surface pressure, specific humidity including
	2103	!-- nudging and subsidence, and geostrophic winds ug, vg
[2696]	2104	io_group_list(3) = init_io_group( &
	2105	in_files = flow_files, &
[3395]	2106	out_vars = [output_var_table(56:64), & ! internal averaged density and pressure profiles
	2107	output_var_table(3:8), output_var_table(9:14), &
	2108	output_var_table(42:42), output_var_table(43:44), &
	2109	output_var_table(49:51), output_var_table(52:54)], &
	2110	in_var_list = (/input_var_table(3), input_var_table(17), & ! T, P, QV
	2111	input_var_table(4) /), &
	2112	kind = 'thermodynamics', &
	2113	n_output_quantities = 4 & ! P, Theta, Rho, qv
[2696]	2114	)
	2115
[3557]	2116	!
	2117	!-- Moved to therodynamic io_group
[3395]	2118	!io_group_list(4) = init_io_group( &
	2119	! in_files = flow_files, &
	2120	! out_vars = [output_var_table(9:14), output_var_table(52:54)], &
	2121	! in_var_list = input_var_table(4:4), &
	2122	! kind = 'scalar' &
	2123	!)
[2696]	2124
[3557]	2125	!
	2126	!-- u and v velocity
[2696]	2127	io_group_list(5) = init_io_group( &
	2128	in_files = flow_files, &
[3395]	2129	out_vars = [output_var_table(15:26), output_var_table(45:46)], &
[2696]	2130	!out_vars = output_var_table(15:20), &
	2131	in_var_list = input_var_table(5:6), &
	2132	!in_var_list = input_var_table(5:5), &
	2133	kind = 'velocities' &
	2134	)
	2135
[3557]	2136	!
	2137	!-- v velocity, deprecated!
[2696]	2138	!io_group_list(6) = init_io_group( &
	2139	! in_files = flow_files, &
	2140	! out_vars = output_var_table(21:26), &
	2141	! in_var_list = input_var_table(6:6), &
	2142	! kind = 'horizontal velocity' &
	2143	!)
	2144	!io_group_list(6) % to_be_processed = .FALSE.
	2145
[3557]	2146	!
	2147	!-- w velocity and subsidence and w nudging
[2696]	2148	io_group_list(7) = init_io_group( &
	2149	in_files = flow_files, &
[3182]	2150	out_vars = [output_var_table(27:32), output_var_table(47:48)], &
[2696]	2151	in_var_list = input_var_table(7:7), &
	2152	kind = 'scalar' &
	2153	)
[3557]	2154	!
	2155	!-- rain
[2696]	2156	io_group_list(8) = init_io_group( &
	2157	in_files = soil_moisture_files, &
	2158	out_vars = output_var_table(33:33), &
	2159	in_var_list = input_var_table(8:8), &
	2160	kind = 'accumulated' &
	2161	)
[3182]	2162	io_group_list(8) % to_be_processed = .FALSE.
[3557]	2163	!
	2164	!-- snow
[2696]	2165	io_group_list(9) = init_io_group( &
	2166	in_files = soil_moisture_files, &
	2167	out_vars = output_var_table(34:34), &
	2168	in_var_list = input_var_table(9:9), &
	2169	kind = 'accumulated' &
	2170	)
[3182]	2171	io_group_list(9) % to_be_processed = .FALSE.
[3557]	2172	!
	2173	!-- graupel
[2696]	2174	io_group_list(10) = init_io_group( &
	2175	in_files = soil_moisture_files, &
	2176	out_vars = output_var_table(35:35), &
	2177	in_var_list = input_var_table(10:10), &
	2178	kind = 'accumulated' &
	2179	)
[3182]	2180	io_group_list(10) % to_be_processed = .FALSE.
[3557]	2181	!
	2182	!-- evapotranspiration
[2696]	2183	io_group_list(11) = init_io_group( &
	2184	in_files = soil_moisture_files, &
	2185	out_vars = output_var_table(37:37), &
	2186	in_var_list = input_var_table(11:11), &
	2187	kind = 'accumulated' &
	2188	)
[3182]	2189	io_group_list(11) % to_be_processed = .FALSE.
[3557]	2190	!
	2191	!-- 2m air temperature
[2696]	2192	io_group_list(12) = init_io_group( &
	2193	in_files = soil_moisture_files, &
	2194	out_vars = output_var_table(36:36), &
	2195	in_var_list = input_var_table(12:12), &
	2196	kind = 'surface' &
	2197	)
[3557]	2198	!
	2199	!-- incoming diffusive sw flux
[2696]	2200	io_group_list(13) = init_io_group( &
	2201	in_files = radiation_files, &
	2202	out_vars = output_var_table(38:38), &
	2203	in_var_list = input_var_table(13:13), &
	2204	kind = 'running average' &
	2205	)
	2206	io_group_list(13) % to_be_processed = .FALSE.
[3557]	2207	!
	2208	!-- incoming direct sw flux
[2696]	2209	io_group_list(14) = init_io_group( &
	2210	in_files = radiation_files, &
	2211	out_vars = output_var_table(39:39), &
	2212	in_var_list = input_var_table(14:14), &
	2213	kind = 'running average' &
	2214	)
	2215	io_group_list(14) % to_be_processed = .FALSE.
[3557]	2216	!
	2217	!-- sw radiation balance
[2696]	2218	io_group_list(15) = init_io_group( &
	2219	in_files = radiation_files, &
	2220	out_vars = output_var_table(40:40), &
	2221	in_var_list = input_var_table(15:15), &
	2222	kind = 'running average' &
	2223	)
[3182]	2224	io_group_list(15) % to_be_processed = .FALSE.
[3557]	2225	!
	2226	!-- lw radiation balance
[2696]	2227	io_group_list(16) = init_io_group( &
	2228	in_files = radiation_files, &
	2229	out_vars = output_var_table(41:41), &
	2230	in_var_list = input_var_table(16:16), &
	2231	kind = 'running average' &
	2232	)
[3182]	2233	io_group_list(16) % to_be_processed = .FALSE.
[2696]	2234
	2235	END SUBROUTINE setup_io_groups
	2236
	2237
[3557]	2238	!------------------------------------------------------------------------------!
	2239	! Description:
	2240	! ------------
	2241	!> Initializes an IO group with the list of input files, the lists of
	2242	!> input and output variables, and sets the number of output quantities.
	2243	!>
	2244	!> In this context, output quantities refers to the physical quantities required
	2245	!> to interpolate netCDF output variables, not the output variables itsleft. For
	2246	!> instance, the output variables ls_forcing_left/_right/_north/_south all rely
	2247	!> on the output quantity Theta.
	2248	!------------------------------------------------------------------------------!
[3395]	2249	FUNCTION init_io_group(in_files, out_vars, in_var_list, kind, &
	2250	n_output_quantities) RESULT(group)
[2696]	2251	CHARACTER(LEN=PATH), INTENT(IN) :: in_files(:)
	2252	CHARACTER(LEN=*), INTENT(IN) :: kind
	2253	TYPE(nc_var), INTENT(IN) :: out_vars(:)
	2254	TYPE(nc_var), INTENT(IN) :: in_var_list(:)
[3395]	2255	INTEGER, OPTIONAL :: n_output_quantities
[2696]	2256
	2257	TYPE(io_group) :: group
	2258
	2259	group % nt = SIZE(in_files)
	2260	group % nv = SIZE(out_vars)
[3395]	2261	group % n_inputs = SIZE(in_var_list)
[2696]	2262	group % kind = TRIM(kind)
[3557]	2263	!
	2264	!-- For the 'thermodynamics' IO group, one quantity more than input variables
	2265	!-- is needed to compute all output variables of the IO group. Concretely, in
	2266	!-- preprocess() the density is computed from T,P or PP,QV in adddition to
	2267	!-- the variables Theta, p, qv. In read_input_variables(),
	2268	!-- n_output_quantities is used to allocate the correct number of input
	2269	!-- buffers.
[3395]	2270	IF ( PRESENT(n_output_quantities) ) THEN
	2271	group % n_output_quantities = n_output_quantities
	2272	ELSE
	2273	group % n_output_quantities = group % n_inputs
	2274	END IF
	2275
	2276	ALLOCATE(group % in_var_list(group % n_inputs))
[2696]	2277	ALLOCATE(group % in_files(group % nt))
	2278	ALLOCATE(group % out_vars(group % nv))
	2279
	2280	group % in_var_list = in_var_list
	2281	group % in_files = in_files
	2282	group % out_vars = out_vars
	2283	group % to_be_processed = .TRUE.
	2284
	2285	END FUNCTION init_io_group
	2286
	2287
	2288	!------------------------------------------------------------------------------!
	2289	! Description:
	2290	! ------------
	2291	!> Deallocates all allocated variables.
	2292	!------------------------------------------------------------------------------!
	2293	SUBROUTINE fini_grids()
	2294
[3395]	2295	CALL report('fini_grids', 'Deallocating grids', cfg % debug)
[2696]	2296
[3182]	2297	DEALLOCATE(x, y, z, xu, yv, zw, z_column, zw_column)
	2298
[2696]	2299	DEALLOCATE(palm_grid%x, palm_grid%y, palm_grid%z, &
	2300	palm_grid%xu, palm_grid%yv, palm_grid%zw, &
	2301	palm_grid%clon, palm_grid%clat, &
	2302	palm_grid%clonu, palm_grid%clatu)
	2303
	2304	DEALLOCATE(palm_intermediate%x, palm_intermediate%y, palm_intermediate%z, &
	2305	palm_intermediate%xu, palm_intermediate%yv, palm_intermediate%zw,&
	2306	palm_intermediate%clon, palm_intermediate%clat, &
	2307	palm_intermediate%clonu, palm_intermediate%clatu)
	2308
[3182]	2309	DEALLOCATE(cosmo_grid%lon, cosmo_grid%lat, &
	2310	cosmo_grid%lonu, cosmo_grid%latv, &
[2696]	2311	cosmo_grid%hfl)
	2312
	2313	END SUBROUTINE fini_grids
	2314
	2315
	2316	!------------------------------------------------------------------------------!
	2317	! Description:
	2318	! ------------
[3557]	2319	!> Initializes the variable list.
[2696]	2320	!------------------------------------------------------------------------------!
[3182]	2321	SUBROUTINE setup_variable_tables(ic_mode)
	2322	CHARACTER(LEN=*), INTENT(IN) :: ic_mode
[3395]	2323	INTEGER :: n_invar = 0 !< number of variables in the input variable table
	2324	INTEGER :: n_outvar = 0 !< number of variables in the output variable table
[2696]	2325	TYPE(nc_var), POINTER :: var
	2326
[3182]	2327	IF (TRIM(cfg % start_date) == '') THEN
[2696]	2328	message = 'Simulation start date has not been set.'
[3615]	2329	CALL inifor_abort('setup_variable_tables', message)
[2696]	2330	END IF
	2331
[3182]	2332	nc_source_text = 'COSMO-DE analysis from ' // TRIM(cfg % start_date)
[2696]	2333
	2334	n_invar = 17
[3395]	2335	n_outvar = 64
[2696]	2336	ALLOCATE( input_var_table(n_invar) )
	2337	ALLOCATE( output_var_table(n_outvar) )
	2338
	2339	!
	2340	!------------------------------------------------------------------------------
	2341	!- Section 1: NetCDF input variables
	2342	!------------------------------------------------------------------------------
	2343	var => input_var_table(1)
	2344	var % name = 'T_SO'
	2345	var % to_be_processed = .TRUE.
	2346	var % is_upside_down = .FALSE.
	2347
	2348	var => input_var_table(2)
	2349	var % name = 'W_SO'
	2350	var % to_be_processed = .TRUE.
	2351	var % is_upside_down = .FALSE.
	2352
	2353	var => input_var_table(3)
	2354	var % name = 'T'
	2355	var % to_be_processed = .TRUE.
	2356	var % is_upside_down = .TRUE.
	2357
	2358	var => input_var_table(4)
	2359	var % name = 'QV'
	2360	var % to_be_processed = .TRUE.
	2361	var % is_upside_down = .TRUE.
	2362
	2363	var => input_var_table(5)
	2364	var % name = 'U'
	2365	var % to_be_processed = .TRUE.
	2366	var % is_upside_down = .TRUE.
	2367
	2368	var => input_var_table(6)
	2369	var % name = 'V'
	2370	var % to_be_processed = .TRUE.
	2371	var % is_upside_down = .TRUE.
	2372
	2373	var => input_var_table(7)
	2374	var % name = 'W'
	2375	var % to_be_processed = .TRUE.
	2376	var % is_upside_down = .TRUE.
	2377
	2378	var => input_var_table(8)
	2379	var % name = 'RAIN_GSP'
	2380	var % to_be_processed = .TRUE.
	2381	var % is_upside_down = .FALSE.
	2382
	2383	var => input_var_table(9)
	2384	var % name = 'SNOW_GSP'
	2385	var % to_be_processed = .TRUE.
	2386	var % is_upside_down = .FALSE.
	2387
	2388	var => input_var_table(10)
	2389	var % name = 'GRAU_GSP'
	2390	var % to_be_processed = .TRUE.
	2391	var % is_upside_down = .FALSE.
	2392
	2393	var => input_var_table(11)
	2394	var % name = 'AEVAP_S'
	2395	var % to_be_processed = .TRUE.
	2396	var % is_upside_down = .FALSE.
	2397
	2398	var => input_var_table(12)
	2399	var % name = 'T_2M'
	2400	var % to_be_processed = .TRUE.
	2401	var % is_upside_down = .FALSE.
	2402
	2403	var => input_var_table(13)
	2404	var % name = 'ASWDIFD_S'
	2405	var % to_be_processed = .TRUE.
	2406	var % is_upside_down = .FALSE.
	2407
	2408	var => input_var_table(14)
	2409	var % name = 'ASWDIR_S'
	2410	var % to_be_processed = .TRUE.
	2411	var % is_upside_down = .FALSE.
	2412
	2413	var => input_var_table(15)
	2414	var % name = 'ASOB_S'
	2415	var % to_be_processed = .TRUE.
	2416	var % is_upside_down = .FALSE.
	2417
	2418	var => input_var_table(16)
	2419	var % name = 'ATHB_S'
	2420	var % to_be_processed = .TRUE.
	2421	var % is_upside_down = .FALSE.
	2422
	2423	var => input_var_table(17)
	2424	var % name = 'P'
	2425	var % to_be_processed = .TRUE.
	2426	var % is_upside_down = .TRUE.
	2427
	2428	!
	2429	!------------------------------------------------------------------------------
	2430	!- Section 2: NetCDF output variables
	2431	!------------------------------------------------------------------------------
[3182]	2432	!
	2433	!------------------------------------------------------------------------------
	2434	! Section 2.1: Realistic forcings, i.e. 3D initial and boundary conditions
	2435	!------------------------------------------------------------------------------
[2696]	2436	output_var_table(1) = init_nc_var( &
	2437	name = 'init_soil_t', &
	2438	std_name = "", &
	2439	long_name = "initial soil temperature", &
	2440	units = "K", &
	2441	kind = "init soil", &
	2442	input_id = 1, &
	2443	output_file = output_file, &
	2444	grid = palm_grid, &
	2445	intermediate_grid = palm_intermediate &
	2446	)
	2447
	2448	output_var_table(2) = init_nc_var( &
	2449	name = 'init_soil_m', &
	2450	std_name = "", &
	2451	long_name = "initial soil moisture", &
	2452	units = "m^3/m^3", &
	2453	kind = "init soil", &
	2454	input_id = 1, &
	2455	output_file = output_file, &
	2456	grid = palm_grid, &
	2457	intermediate_grid = palm_intermediate &
	2458	)
	2459
	2460	output_var_table(3) = init_nc_var( &
[3182]	2461	name = 'init_atmosphere_pt', &
[2696]	2462	std_name = "", &
	2463	long_name = "initial potential temperature", &
	2464	units = "K", &
	2465	kind = "init scalar", &
	2466	input_id = 1, & ! first in (T, p) IO group
	2467	output_file = output_file, &
	2468	grid = palm_grid, &
	2469	intermediate_grid = palm_intermediate, &
[3182]	2470	is_profile = (TRIM(ic_mode) == 'profile') &
[2696]	2471	)
[3182]	2472	IF (TRIM(ic_mode) == 'profile') THEN
[3613]	2473	output_var_table(3) % averaging_grid => averaged_initial_scalar_profile
[2696]	2474	END IF
	2475
	2476	output_var_table(4) = init_nc_var( &
	2477	name = 'ls_forcing_left_pt', &
	2478	std_name = "", &
	2479	long_name = "large-scale forcing for left model boundary for the potential temperature", &
	2480	units = "K", &
	2481	kind = "left scalar", &
	2482	input_id = 1, &
	2483	grid = scalars_west_grid, &
	2484	intermediate_grid = scalars_west_intermediate, &
	2485	output_file = output_file &
	2486	)
	2487
	2488	output_var_table(5) = init_nc_var( &
	2489	name = 'ls_forcing_right_pt', &
	2490	std_name = "", &
	2491	long_name = "large-scale forcing for right model boundary for the potential temperature", &
	2492	units = "K", &
	2493	kind = "right scalar", &
	2494	input_id = 1, &
	2495	grid = scalars_east_grid, &
	2496	intermediate_grid = scalars_east_intermediate, &
	2497	output_file = output_file &
	2498	)
	2499
	2500	output_var_table(6) = init_nc_var( &
	2501	name = 'ls_forcing_north_pt', &
	2502	std_name = "", &
	2503	long_name = "large-scale forcing for north model boundary for the potential temperature", &
	2504	units = "K", &
	2505	kind = "north scalar", &
	2506	input_id = 1, &
	2507	grid = scalars_north_grid, &
	2508	intermediate_grid = scalars_north_intermediate, &
	2509	output_file = output_file &
	2510	)
	2511
	2512	output_var_table(7) = init_nc_var( &
	2513	name = 'ls_forcing_south_pt', &
	2514	std_name = "", &
	2515	long_name = "large-scale forcing for south model boundary for the potential temperature", &
	2516	units = "K", &
	2517	kind = "south scalar", &
	2518	input_id = 1, &
	2519	grid = scalars_south_grid, &
	2520	intermediate_grid = scalars_south_intermediate, &
	2521	output_file = output_file &
	2522	)
	2523
	2524	output_var_table(8) = init_nc_var( &
	2525	name = 'ls_forcing_top_pt', &
	2526	std_name = "", &
	2527	long_name = "large-scale forcing for top model boundary for the potential temperature", &
	2528	units = "K", &
	2529	kind = "top scalar", &
	2530	input_id = 1, &
	2531	grid = scalars_top_grid, &
	2532	intermediate_grid = scalars_top_intermediate, &
	2533	output_file = output_file &
	2534	)
	2535
	2536	output_var_table(9) = init_nc_var( &
[3182]	2537	name = 'init_atmosphere_qv', &
[2696]	2538	std_name = "", &
	2539	long_name = "initial specific humidity", &
	2540	units = "kg/kg", &
	2541	kind = "init scalar", &
[3395]	2542	input_id = 3, &
[2696]	2543	output_file = output_file, &
	2544	grid = palm_grid, &
	2545	intermediate_grid = palm_intermediate, &
[3182]	2546	is_profile = (TRIM(ic_mode) == 'profile') &
[2696]	2547	)
[3182]	2548	IF (TRIM(ic_mode) == 'profile') THEN
[3613]	2549	output_var_table(9) % averaging_grid => averaged_initial_scalar_profile
[2696]	2550	END IF
	2551
	2552	output_var_table(10) = init_nc_var( &
	2553	name = 'ls_forcing_left_qv', &
	2554	std_name = "", &
	2555	long_name = "large-scale forcing for left model boundary for the specific humidity", &
	2556	units = "kg/kg", &
	2557	kind = "left scalar", &
[3395]	2558	input_id = 3, &
[2696]	2559	output_file = output_file, &
	2560	grid = scalars_west_grid, &
	2561	intermediate_grid = scalars_west_intermediate &
	2562	)
	2563
	2564	output_var_table(11) = init_nc_var( &
	2565	name = 'ls_forcing_right_qv', &
	2566	std_name = "", &
	2567	long_name = "large-scale forcing for right model boundary for the specific humidity", &
	2568	units = "kg/kg", &
	2569	kind = "right scalar", &
[3395]	2570	input_id = 3, &
[2696]	2571	output_file = output_file, &
	2572	grid = scalars_east_grid, &
	2573	intermediate_grid = scalars_east_intermediate &
	2574	)
	2575
	2576	output_var_table(12) = init_nc_var( &
	2577	name = 'ls_forcing_north_qv', &
	2578	std_name = "", &
	2579	long_name = "large-scale forcing for north model boundary for the specific humidity", &
	2580	units = "kg/kg", &
	2581	kind = "north scalar", &
[3395]	2582	input_id = 3, &
[2696]	2583	output_file = output_file, &
	2584	grid = scalars_north_grid, &
	2585	intermediate_grid = scalars_north_intermediate &
	2586	)
	2587
	2588	output_var_table(13) = init_nc_var( &
	2589	name = 'ls_forcing_south_qv', &
	2590	std_name = "", &
	2591	long_name = "large-scale forcing for south model boundary for the specific humidity", &
	2592	units = "kg/kg", &
	2593	kind = "south scalar", &
[3395]	2594	input_id = 3, &
[2696]	2595	output_file = output_file, &
	2596	grid = scalars_south_grid, &
	2597	intermediate_grid = scalars_south_intermediate &
	2598	)
	2599
	2600	output_var_table(14) = init_nc_var( &
	2601	name = 'ls_forcing_top_qv', &
	2602	std_name = "", &
	2603	long_name = "large-scale forcing for top model boundary for the specific humidity", &
	2604	units = "kg/kg", &
	2605	kind = "top scalar", &
[3395]	2606	input_id = 3, &
[2696]	2607	output_file = output_file, &
	2608	grid = scalars_top_grid, &
	2609	intermediate_grid = scalars_top_intermediate &
	2610	)
	2611
	2612	output_var_table(15) = init_nc_var( &
[3182]	2613	name = 'init_atmosphere_u', &
[2696]	2614	std_name = "", &
	2615	long_name = "initial wind component in x direction", &
	2616	units = "m/s", &
	2617	kind = "init u", &
	2618	input_id = 1, & ! first in (U, V) I/O group
	2619	output_file = output_file, &
	2620	grid = u_initial_grid, &
	2621	intermediate_grid = u_initial_intermediate, &
[3182]	2622	is_profile = (TRIM(ic_mode) == 'profile') &
[2696]	2623	)
[3182]	2624	IF (TRIM(ic_mode) == 'profile') THEN
[3613]	2625	output_var_table(15) % averaging_grid => averaged_initial_scalar_profile
[2696]	2626	END IF
	2627
	2628	output_var_table(16) = init_nc_var( &
	2629	name = 'ls_forcing_left_u', &
	2630	std_name = "", &
	2631	long_name = "large-scale forcing for left model boundary for the wind component in x direction", &
	2632	units = "m/s", &
	2633	kind = "left u", &
	2634	input_id = 1, & ! first in (U, V) I/O group
	2635	output_file = output_file, &
	2636	grid = u_west_grid, &
	2637	intermediate_grid = u_west_intermediate &
	2638	)
	2639
	2640	output_var_table(17) = init_nc_var( &
	2641	name = 'ls_forcing_right_u', &
	2642	std_name = "", &
	2643	long_name = "large-scale forcing for right model boundary for the wind component in x direction", &
	2644	units = "m/s", &
	2645	kind = "right u", &
	2646	input_id = 1, & ! first in (U, V) I/O group
	2647	output_file = output_file, &
	2648	grid = u_east_grid, &
	2649	intermediate_grid = u_east_intermediate &
	2650	)
	2651
	2652	output_var_table(18) = init_nc_var( &
	2653	name = 'ls_forcing_north_u', &
	2654	std_name = "", &
	2655	long_name = "large-scale forcing for north model boundary for the wind component in x direction", &
	2656	units = "m/s", &
	2657	kind = "north u", &
	2658	input_id = 1, & ! first in (U, V) I/O group
	2659	output_file = output_file, &
	2660	grid = u_north_grid, &
	2661	intermediate_grid = u_north_intermediate &
	2662	)
	2663
	2664	output_var_table(19) = init_nc_var( &
	2665	name = 'ls_forcing_south_u', &
	2666	std_name = "", &
	2667	long_name = "large-scale forcing for south model boundary for the wind component in x direction", &
	2668	units = "m/s", &
	2669	kind = "south u", &
	2670	input_id = 1, & ! first in (U, V) I/O group
	2671	output_file = output_file, &
	2672	grid = u_south_grid, &
	2673	intermediate_grid = u_south_intermediate &
	2674	)
	2675
	2676	output_var_table(20) = init_nc_var( &
	2677	name = 'ls_forcing_top_u', &
	2678	std_name = "", &
	2679	long_name = "large-scale forcing for top model boundary for the wind component in x direction", &
	2680	units = "m/s", &
	2681	kind = "top u", &
	2682	input_id = 1, & ! first in (U, V) I/O group
	2683	output_file = output_file, &
	2684	grid = u_top_grid, &
	2685	intermediate_grid = u_top_intermediate &
	2686	)
	2687
	2688	output_var_table(21) = init_nc_var( &
[3182]	2689	name = 'init_atmosphere_v', &
[2696]	2690	std_name = "", &
	2691	long_name = "initial wind component in y direction", &
	2692	units = "m/s", &
	2693	kind = "init v", &
	2694	input_id = 2, & ! second in (U, V) I/O group
	2695	output_file = output_file, &
	2696	grid = v_initial_grid, &
	2697	intermediate_grid = v_initial_intermediate, &
[3182]	2698	is_profile = (TRIM(ic_mode) == 'profile') &
[2696]	2699	)
[3182]	2700	IF (TRIM(ic_mode) == 'profile') THEN
[3613]	2701	output_var_table(21) % averaging_grid => averaged_initial_scalar_profile
[2696]	2702	END IF
	2703
	2704	output_var_table(22) = init_nc_var( &
	2705	name = 'ls_forcing_left_v', &
	2706	std_name = "", &
	2707	long_name = "large-scale forcing for left model boundary for the wind component in y direction", &
	2708	units = "m/s", &
	2709	kind = "right v", &
	2710	input_id = 2, & ! second in (U, V) I/O group
	2711	output_file = output_file, &
	2712	grid = v_west_grid, &
	2713	intermediate_grid = v_west_intermediate &
	2714	)
	2715
	2716	output_var_table(23) = init_nc_var( &
	2717	name = 'ls_forcing_right_v', &
	2718	std_name = "", &
	2719	long_name = "large-scale forcing for right model boundary for the wind component in y direction", &
	2720	units = "m/s", &
	2721	kind = "right v", &
	2722	input_id = 2, & ! second in (U, V) I/O group
	2723	output_file = output_file, &
	2724	grid = v_east_grid, &
	2725	intermediate_grid = v_east_intermediate &
	2726	)
	2727
	2728	output_var_table(24) = init_nc_var( &
	2729	name = 'ls_forcing_north_v', &
	2730	std_name = "", &
	2731	long_name = "large-scale forcing for north model boundary for the wind component in y direction", &
	2732	units = "m/s", &
	2733	kind = "north v", &
	2734	input_id = 2, & ! second in (U, V) I/O group
	2735	output_file = output_file, &
	2736	grid = v_north_grid, &
	2737	intermediate_grid = v_north_intermediate &
	2738	)
	2739
	2740	output_var_table(25) = init_nc_var( &
	2741	name = 'ls_forcing_south_v', &
	2742	std_name = "", &
	2743	long_name = "large-scale forcing for south model boundary for the wind component in y direction", &
	2744	units = "m/s", &
	2745	kind = "south v", &
	2746	input_id = 2, & ! second in (U, V) I/O group
	2747	output_file = output_file, &
	2748	grid = v_south_grid, &
	2749	intermediate_grid = v_south_intermediate &
	2750	)
	2751
	2752	output_var_table(26) = init_nc_var( &
	2753	name = 'ls_forcing_top_v', &
	2754	std_name = "", &
	2755	long_name = "large-scale forcing for top model boundary for the wind component in y direction", &
	2756	units = "m/s", &
	2757	kind = "top v", &
	2758	input_id = 2, & ! second in (U, V) I/O group
	2759	output_file = output_file, &
	2760	grid = v_top_grid, &
	2761	intermediate_grid = v_top_intermediate &
	2762	)
	2763
	2764	output_var_table(27) = init_nc_var( &
[3182]	2765	name = 'init_atmosphere_w', &
[2696]	2766	std_name = "", &
	2767	long_name = "initial wind component in z direction", &
	2768	units = "m/s", &
	2769	kind = "init w", &
	2770	input_id = 1, &
	2771	output_file = output_file, &
	2772	grid = w_initial_grid, &
	2773	intermediate_grid = w_initial_intermediate, &
[3182]	2774	is_profile = (TRIM(ic_mode) == 'profile') &
[2696]	2775	)
[3182]	2776	IF (TRIM(ic_mode) == 'profile') THEN
[3613]	2777	output_var_table(27) % averaging_grid => averaged_initial_w_profile
[2696]	2778	END IF
	2779
	2780	output_var_table(28) = init_nc_var( &
	2781	name = 'ls_forcing_left_w', &
	2782	std_name = "", &
	2783	long_name = "large-scale forcing for left model boundary for the wind component in z direction", &
	2784	units = "m/s", &
	2785	kind = "left w", &
	2786	input_id = 1, &
	2787	output_file = output_file, &
	2788	grid = w_west_grid, &
	2789	intermediate_grid = w_west_intermediate &
	2790	)
	2791
	2792	output_var_table(29) = init_nc_var( &
	2793	name = 'ls_forcing_right_w', &
	2794	std_name = "", &
	2795	long_name = "large-scale forcing for right model boundary for the wind component in z direction", &
	2796	units = "m/s", &
	2797	kind = "right w", &
	2798	input_id = 1, &
	2799	output_file = output_file, &
	2800	grid = w_east_grid, &
	2801	intermediate_grid = w_east_intermediate &
	2802	)
	2803
	2804	output_var_table(30) = init_nc_var( &
	2805	name = 'ls_forcing_north_w', &
	2806	std_name = "", &
	2807	long_name = "large-scale forcing for north model boundary for the wind component in z direction", &
	2808	units = "m/s", &
	2809	kind = "north w", &
	2810	input_id = 1, &
	2811	output_file = output_file, &
	2812	grid = w_north_grid, &
	2813	intermediate_grid = w_north_intermediate &
	2814	)
	2815
	2816	output_var_table(31) = init_nc_var( &
	2817	name = 'ls_forcing_south_w', &
	2818	std_name = "", &
	2819	long_name = "large-scale forcing for south model boundary for the wind component in z direction", &
	2820	units = "m/s", &
	2821	kind = "south w", &
	2822	input_id = 1, &
	2823	output_file = output_file, &
	2824	grid = w_south_grid, &
	2825	intermediate_grid = w_south_intermediate &
	2826	)
	2827
	2828	output_var_table(32) = init_nc_var( &
	2829	name = 'ls_forcing_top_w', &
	2830	std_name = "", &
	2831	long_name = "large-scale forcing for top model boundary for the wind component in z direction", &
	2832	units = "m/s", &
	2833	kind = "top w", &
	2834	input_id = 1, &
	2835	output_file = output_file, &
	2836	grid = w_top_grid, &
	2837	intermediate_grid = w_top_intermediate &
	2838	)
	2839
	2840	output_var_table(33) = init_nc_var( &
	2841	name = 'ls_forcing_soil_rain', &
	2842	std_name = "", &
	2843	long_name = "large-scale forcing rain", &
	2844	units = "kg/m2", &
	2845	kind = "surface forcing", &
	2846	input_id = 1, &
	2847	output_file = output_file, &
	2848	grid = palm_grid, &
	2849	intermediate_grid = palm_intermediate &
	2850	)
	2851
	2852	output_var_table(34) = init_nc_var( &
	2853	name = 'ls_forcing_soil_snow', &
	2854	std_name = "", &
	2855	long_name = "large-scale forcing snow", &
	2856	units = "kg/m2", &
	2857	kind = "surface forcing", &
	2858	input_id = 1, &
	2859	output_file = output_file, &
	2860	grid = palm_grid, &
	2861	intermediate_grid = palm_intermediate &
	2862	)
	2863
	2864	output_var_table(35) = init_nc_var( &
	2865	name = 'ls_forcing_soil_graupel', &
	2866	std_name = "", &
	2867	long_name = "large-scale forcing graupel", &
	2868	units = "kg/m2", &
	2869	kind = "surface forcing", &
	2870	input_id = 1, &
	2871	output_file = output_file, &
	2872	grid = palm_grid, &
	2873	intermediate_grid = palm_intermediate &
	2874	)
	2875
	2876	output_var_table(36) = init_nc_var( &
	2877	name = 'ls_forcing_soil_t_2m', &
	2878	std_name = "", &
	2879	long_name = "large-scale forcing 2m air temperature", &
	2880	units = "kg/m2", &
	2881	kind = "surface forcing", &
	2882	input_id = 1, &
	2883	output_file = output_file, &
	2884	grid = palm_grid, &
	2885	intermediate_grid = palm_intermediate &
	2886	)
	2887
	2888	output_var_table(37) = init_nc_var( &
	2889	name = 'ls_forcing_soil_evap', &
	2890	std_name = "", &
	2891	long_name = "large-scale forcing evapo-transpiration", &
	2892	units = "kg/m2", &
	2893	kind = "surface forcing", &
	2894	input_id = 1, &
	2895	output_file = output_file, &
	2896	grid = palm_grid, &
	2897	intermediate_grid = palm_intermediate &
	2898	)
[3557]	2899
[2696]	2900	output_var_table(38) = init_nc_var( &
	2901	name = 'rad_swd_dif_0', &
	2902	std_name = "", &
	2903	long_name = "incoming diffuse shortwave radiative flux at the surface", &
	2904	units = "W/m2", &
	2905	kind = "surface forcing", &
	2906	input_id = 1, &
	2907	output_file = output_file, &
	2908	grid = palm_grid, &
	2909	intermediate_grid = palm_intermediate &
	2910	)
	2911
	2912	output_var_table(39) = init_nc_var( &
	2913	name = 'rad_swd_dir_0', &
	2914	std_name = "", &
	2915	long_name = "incoming direct shortwave radiative flux at the surface", &
	2916	units = "W/m2", &
	2917	kind = "surface forcing", &
	2918	input_id = 1, &
	2919	output_file = output_file, &
	2920	grid = palm_grid, &
	2921	intermediate_grid = palm_intermediate &
	2922	)
	2923
	2924	output_var_table(40) = init_nc_var( &
	2925	name = 'rad_sw_bal_0', &
	2926	std_name = "", &
	2927	long_name = "shortwave radiation balance at the surface", &
	2928	units = "W/m2", &
	2929	kind = "surface forcing", &
	2930	input_id = 1, &
	2931	output_file = output_file, &
	2932	grid = palm_grid, &
	2933	intermediate_grid = palm_intermediate &
	2934	)
	2935
	2936	output_var_table(41) = init_nc_var( &
	2937	name = 'rad_lw_bal_0', &
	2938	std_name = "", &
	2939	long_name = "longwave radiation balance at the surface", &
	2940	units = "W/m2", &
	2941	kind = "surface forcing", &
	2942	input_id = 1, &
	2943	output_file = output_file, &
	2944	grid = palm_grid, &
	2945	intermediate_grid = palm_intermediate &
	2946	)
[3182]	2947	!
	2948	!------------------------------------------------------------------------------
	2949	! Section 2.2: Idealized large-scale forcings
	2950	!------------------------------------------------------------------------------
[2696]	2951	output_var_table(42) = init_nc_var( &
	2952	name = 'surface_forcing_surface_pressure', &
	2953	std_name = "", &
	2954	long_name = "surface pressure", &
	2955	units = "Pa", &
	2956	kind = "time series", &
	2957	input_id = 2, & ! second in (T, p) I/O group
	2958	output_file = output_file, &
	2959	grid = palm_grid, &
	2960	intermediate_grid = palm_intermediate &
	2961	)
[3395]	2962	output_var_table(42) % averaging_grid => averaged_scalar_profile
[2696]	2963
	2964	output_var_table(43) = init_nc_var( &
	2965	name = 'ls_forcing_ug', &
	2966	std_name = "", &
	2967	long_name = "geostrophic wind (u component)", &
	2968	units = "m/s", &
[3395]	2969	kind = "geostrophic", &
[2696]	2970	input_id = 1, &
	2971	output_file = output_file, &
[3395]	2972	grid = averaged_scalar_profile, &
	2973	intermediate_grid = averaged_scalar_profile &
[2696]	2974	)
	2975
	2976	output_var_table(44) = init_nc_var( &
	2977	name = 'ls_forcing_vg', &
	2978	std_name = "", &
	2979	long_name = "geostrophic wind (v component)", &
	2980	units = "m/s", &
[3395]	2981	kind = "geostrophic", &
[2696]	2982	input_id = 1, &
	2983	output_file = output_file, &
[3395]	2984	grid = averaged_scalar_profile, &
	2985	intermediate_grid = averaged_scalar_profile &
[2696]	2986	)
	2987
[3182]	2988	output_var_table(45) = init_nc_var( &
	2989	name = 'nudging_u', &
	2990	std_name = "", &
	2991	long_name = "wind component in x direction", &
	2992	units = "m/s", &
[3395]	2993	kind = "geostrophic", &
[3182]	2994	input_id = 1, &
	2995	output_file = output_file, &
[3395]	2996	grid = averaged_scalar_profile, &
	2997	intermediate_grid = averaged_scalar_profile &
[3182]	2998	)
	2999	output_var_table(45) % to_be_processed = ls_forcing_variables_required
	3000
	3001	output_var_table(46) = init_nc_var( &
	3002	name = 'nudging_v', &
	3003	std_name = "", &
	3004	long_name = "wind component in y direction", &
	3005	units = "m/s", &
	3006	kind = "large-scale scalar forcing", &
	3007	input_id = 1, &
	3008	output_file = output_file, &
[3395]	3009	grid = averaged_scalar_profile, &
	3010	intermediate_grid = averaged_scalar_profile &
[3182]	3011	)
	3012	output_var_table(46) % to_be_processed = ls_forcing_variables_required
	3013
	3014	output_var_table(47) = init_nc_var( &
	3015	name = 'ls_forcing_sub_w', &
	3016	std_name = "", &
	3017	long_name = "subsidence velocity of w", &
	3018	units = "m/s", &
	3019	kind = "large-scale w forcing", &
	3020	input_id = 1, &
	3021	output_file = output_file, &
[3395]	3022	grid = averaged_scalar_profile, &
	3023	intermediate_grid = averaged_scalar_profile &
[3182]	3024	)
	3025	output_var_table(47) % to_be_processed = ls_forcing_variables_required
	3026
	3027	output_var_table(48) = init_nc_var( &
	3028	name = 'nudging_w', &
	3029	std_name = "", &
	3030	long_name = "wind component in w direction", &
	3031	units = "m/s", &
	3032	kind = "large-scale w forcing", &
	3033	input_id = 1, &
	3034	output_file = output_file, &
[3395]	3035	grid = averaged_w_profile, &
	3036	intermediate_grid = averaged_w_profile &
[3182]	3037	)
	3038	output_var_table(48) % to_be_processed = ls_forcing_variables_required
	3039
	3040
	3041	output_var_table(49) = init_nc_var( &
	3042	name = 'ls_forcing_adv_pt', &
	3043	std_name = "", &
	3044	long_name = "advection of potential temperature", &
	3045	units = "K/s", &
	3046	kind = "large-scale scalar forcing", &
	3047	input_id = 1, &
	3048	output_file = output_file, &
[3395]	3049	grid = averaged_scalar_profile, &
	3050	intermediate_grid = averaged_scalar_profile &
[3182]	3051	)
	3052	output_var_table(49) % to_be_processed = ls_forcing_variables_required
	3053
	3054	output_var_table(50) = init_nc_var( &
	3055	name = 'ls_forcing_sub_pt', &
	3056	std_name = "", &
	3057	long_name = "subsidence velocity of potential temperature", &
	3058	units = "K/s", &
	3059	kind = "large-scale scalar forcing", &
	3060	input_id = 1, &
	3061	output_file = output_file, &
[3395]	3062	grid = averaged_scalar_profile, &
	3063	intermediate_grid = averaged_scalar_profile &
[3182]	3064	)
	3065	output_var_table(50) % to_be_processed = ls_forcing_variables_required
	3066
	3067	output_var_table(51) = init_nc_var( &
	3068	name = 'nudging_pt', &
	3069	std_name = "", &
	3070	long_name = "potential temperature", &
	3071	units = "K", &
	3072	kind = "large-scale scalar forcing", &
	3073	input_id = 1, &
	3074	output_file = output_file, &
[3395]	3075	grid = averaged_scalar_profile, &
	3076	intermediate_grid = averaged_scalar_profile &
[3182]	3077	)
	3078	output_var_table(51) % to_be_processed = ls_forcing_variables_required
	3079
	3080	output_var_table(52) = init_nc_var( &
	3081	name = 'ls_forcing_adv_qv', &
	3082	std_name = "", &
	3083	long_name = "advection of specific humidity", &
	3084	units = "kg/kg/s", &
	3085	kind = "large-scale scalar forcing", &
[3395]	3086	input_id = 3, &
[3182]	3087	output_file = output_file, &
[3395]	3088	grid = averaged_scalar_profile, &
	3089	intermediate_grid = averaged_scalar_profile &
[3182]	3090	)
	3091	output_var_table(52) % to_be_processed = ls_forcing_variables_required
	3092
	3093
	3094	output_var_table(53) = init_nc_var( &
	3095	name = 'ls_forcing_sub_qv', &
	3096	std_name = "", &
	3097	long_name = "subsidence velocity of specific humidity", &
	3098	units = "kg/kg/s", &
	3099	kind = "large-scale scalar forcing", &
[3395]	3100	input_id = 3, &
[3182]	3101	output_file = output_file, &
[3395]	3102	grid = averaged_scalar_profile, &
	3103	intermediate_grid = averaged_scalar_profile &
[3182]	3104	)
	3105	output_var_table(53) % to_be_processed = ls_forcing_variables_required
	3106
	3107	output_var_table(54) = init_nc_var( &
	3108	name = 'nudging_qv', &
	3109	std_name = "", &
	3110	long_name = "specific humidity", &
	3111	units = "kg/kg", &
	3112	kind = "large-scale scalar forcing", &
[3395]	3113	input_id = 3, &
[3182]	3114	output_file = output_file, &
[3395]	3115	grid = averaged_scalar_profile, &
	3116	intermediate_grid = averaged_scalar_profile &
[3182]	3117	)
	3118	output_var_table(54) % to_be_processed = ls_forcing_variables_required
	3119
	3120	output_var_table(55) = init_nc_var( &
	3121	name = 'nudging_tau', &
	3122	std_name = "", &
	3123	long_name = "nudging relaxation time scale", &
	3124	units = "s", &
	3125	kind = "constant scalar profile", &
	3126	input_id = 1, &
	3127	output_file = output_file, &
[3395]	3128	grid = averaged_scalar_profile, &
	3129	intermediate_grid = averaged_scalar_profile &
[3182]	3130	)
	3131	output_var_table(55) % to_be_processed = ls_forcing_variables_required
	3132
	3133
[3395]	3134	output_var_table(56) = init_nc_var( &
	3135	name = 'internal_density_centre', &
	3136	std_name = "", &
	3137	long_name = "", &
	3138	units = "", &
	3139	kind = "internal profile", &
	3140	input_id = 4, &
	3141	output_file = output_file, &
	3142	grid = averaged_scalar_profile, &
	3143	intermediate_grid = averaged_scalar_profile &
	3144	)
	3145	output_var_table(56) % averaging_grid => averaged_scalar_profile
	3146
	3147
	3148	output_var_table(57) = init_nc_var( &
	3149	name = 'internal_density_north', &
	3150	std_name = "", &
	3151	long_name = "", &
	3152	units = "", &
	3153	kind = "internal profile", &
	3154	input_id = 4, &
	3155	output_file = output_file, &
	3156	grid = north_averaged_scalar_profile, &
	3157	intermediate_grid = north_averaged_scalar_profile &
	3158	)
	3159	output_var_table(57) % averaging_grid => north_averaged_scalar_profile
	3160	output_var_table(57) % to_be_processed = .NOT. cfg % ug_is_set
	3161
	3162
	3163	output_var_table(58) = init_nc_var( &
	3164	name = 'internal_density_south', &
	3165	std_name = "", &
	3166	long_name = "", &
	3167	units = "", &
	3168	kind = "internal profile", &
	3169	input_id = 4, &
	3170	output_file = output_file, &
	3171	grid = south_averaged_scalar_profile, &
	3172	intermediate_grid = south_averaged_scalar_profile &
	3173	)
	3174	output_var_table(58) % averaging_grid => south_averaged_scalar_profile
	3175	output_var_table(58) % to_be_processed = .NOT. cfg % ug_is_set
	3176
	3177
	3178	output_var_table(59) = init_nc_var( &
	3179	name = 'internal_density_east', &
	3180	std_name = "", &
	3181	long_name = "", &
	3182	units = "", &
	3183	kind = "internal profile", &
	3184	input_id = 4, &
	3185	output_file = output_file, &
	3186	grid = east_averaged_scalar_profile, &
	3187	intermediate_grid = east_averaged_scalar_profile &
	3188	)
	3189	output_var_table(59) % averaging_grid => east_averaged_scalar_profile
	3190	output_var_table(59) % to_be_processed = .NOT. cfg % ug_is_set
	3191
	3192
	3193	output_var_table(60) = init_nc_var( &
	3194	name = 'internal_density_west', &
	3195	std_name = "", &
	3196	long_name = "", &
	3197	units = "", &
	3198	kind = "internal profile", &
	3199	input_id = 4, &
	3200	output_file = output_file, &
	3201	grid = west_averaged_scalar_profile, &
	3202	intermediate_grid = west_averaged_scalar_profile &
	3203	)
	3204	output_var_table(60) % averaging_grid => west_averaged_scalar_profile
	3205	output_var_table(60) % to_be_processed = .NOT. cfg % ug_is_set
	3206
	3207	output_var_table(61) = init_nc_var( &
	3208	name = 'internal_pressure_north', &
	3209	std_name = "", &
	3210	long_name = "", &
	3211	units = "", &
	3212	kind = "internal profile", &
	3213	input_id = 2, &
	3214	output_file = output_file, &
	3215	grid = north_averaged_scalar_profile, &
	3216	intermediate_grid = north_averaged_scalar_profile &
	3217	)
	3218	output_var_table(61) % averaging_grid => north_averaged_scalar_profile
	3219	output_var_table(61) % to_be_processed = .NOT. cfg % ug_is_set
	3220
	3221
	3222	output_var_table(62) = init_nc_var( &
	3223	name = 'internal_pressure_south', &
	3224	std_name = "", &
	3225	long_name = "", &
	3226	units = "", &
	3227	kind = "internal profile", &
	3228	input_id = 2, &
	3229	output_file = output_file, &
	3230	grid = south_averaged_scalar_profile, &
	3231	intermediate_grid = south_averaged_scalar_profile &
	3232	)
	3233	output_var_table(62) % averaging_grid => south_averaged_scalar_profile
	3234	output_var_table(62) % to_be_processed = .NOT. cfg % ug_is_set
	3235
	3236
	3237	output_var_table(63) = init_nc_var( &
	3238	name = 'internal_pressure_east', &
	3239	std_name = "", &
	3240	long_name = "", &
	3241	units = "", &
	3242	kind = "internal profile", &
	3243	input_id = 2, &
	3244	output_file = output_file, &
	3245	grid = east_averaged_scalar_profile, &
	3246	intermediate_grid = east_averaged_scalar_profile &
	3247	)
	3248	output_var_table(63) % averaging_grid => east_averaged_scalar_profile
	3249	output_var_table(63) % to_be_processed = .NOT. cfg % ug_is_set
	3250
	3251
	3252	output_var_table(64) = init_nc_var( &
	3253	name = 'internal_pressure_west', &
	3254	std_name = "", &
	3255	long_name = "", &
	3256	units = "", &
	3257	kind = "internal profile", &
	3258	input_id = 2, &
	3259	output_file = output_file, &
	3260	grid = west_averaged_scalar_profile, &
	3261	intermediate_grid = west_averaged_scalar_profile &
	3262	)
	3263	output_var_table(64) % averaging_grid => west_averaged_scalar_profile
	3264	output_var_table(64) % to_be_processed = .NOT. cfg % ug_is_set
	3265
[3557]	3266	!
	3267	!-- Attributes shared among all variables
[2696]	3268	output_var_table(:) % source = nc_source_text
	3269
[3182]	3270
[2696]	3271	END SUBROUTINE setup_variable_tables
	3272
	3273
	3274	!------------------------------------------------------------------------------!
	3275	! Description:
	3276	! ------------
[3557]	3277	!> Initializes an nc_var varible with the given parameters. The 'kind'
[2696]	3278	!> parameter is used to infer the correct netCDF IDs and the level of detail,
	3279	!> 'lod', as defined by the PALM-4U input data standard.
	3280	!------------------------------------------------------------------------------!
	3281	FUNCTION init_nc_var(name, std_name, long_name, units, kind, input_id, &
[3395]	3282	grid, intermediate_grid, output_file, is_profile) &
	3283	RESULT(var)
[2696]	3284
	3285	CHARACTER(LEN=*), INTENT(IN) :: name, std_name, long_name, units, kind
	3286	INTEGER, INTENT(IN) :: input_id
	3287	TYPE(grid_definition), INTENT(IN), TARGET :: grid, intermediate_grid
	3288	TYPE(nc_file), INTENT(IN) :: output_file
	3289	LOGICAL, INTENT(IN), OPTIONAL :: is_profile
	3290
	3291	CHARACTER(LEN=LNAME) :: out_var_kind
	3292	TYPE(nc_var) :: var
	3293
	3294	out_var_kind = TRIM(kind)
	3295
	3296	IF (PRESENT(is_profile)) THEN
	3297	IF (is_profile) out_var_kind = TRIM(kind) // ' profile'
	3298	END IF
	3299
	3300	var % name = name
	3301	var % standard_name = std_name
	3302	var % long_name = long_name
	3303	var % units = units
	3304	var % kind = TRIM(out_var_kind)
	3305	var % input_id = input_id
	3306	var % nt = SIZE (output_file % time)
	3307	var % grid => grid
	3308	var % intermediate_grid => intermediate_grid
	3309
	3310	SELECT CASE( TRIM(out_var_kind) )
	3311
[3557]	3312	!
	3313	!-- TODO: Using global module variables 'init_variables_required' and
	3314	!-- TODO: 'boundary_variables_required'. Encapsulate in settings type
	3315	!-- TODO: and pass into init_nc_var.
[2696]	3316	CASE( 'init soil' )
	3317	var % nt = 1
	3318	var % lod = 2
	3319	var % ndim = 3
	3320	var % dimids(1:3) = output_file % dimids_soil
	3321	var % dimvarids(1:3) = output_file % dimvarids_soil
	3322	var % to_be_processed = init_variables_required
[3395]	3323	var % is_internal = .FALSE.
[2696]	3324	var % task = "interpolate_2d"
	3325
	3326	CASE( 'init scalar' )
	3327	var % nt = 1
	3328	var % lod = 2
	3329	var % ndim = 3
	3330	var % dimids(1:3) = output_file % dimids_scl
	3331	var % dimvarids(1:3) = output_file % dimvarids_scl
	3332	var % to_be_processed = init_variables_required
[3395]	3333	var % is_internal = .FALSE.
[2696]	3334	var % task = "interpolate_3d"
	3335
	3336	CASE( 'init u' )
	3337	var % nt = 1
	3338	var % lod = 2
	3339	var % ndim = 3
	3340	var % dimids(1) = output_file % dimids_vel(1)
	3341	var % dimids(2) = output_file % dimids_scl(2)
	3342	var % dimids(3) = output_file % dimids_scl(3)
	3343	var % dimvarids(1) = output_file % dimvarids_vel(1)
	3344	var % dimvarids(2) = output_file % dimvarids_scl(2)
	3345	var % dimvarids(3) = output_file % dimvarids_scl(3)
	3346	var % to_be_processed = init_variables_required
[3395]	3347	var % is_internal = .FALSE.
[2696]	3348	var % task = "interpolate_3d"
	3349
	3350	CASE( 'init v' )
	3351	var % nt = 1
	3352	var % lod = 2
	3353	var % ndim = 3
	3354	var % dimids(1) = output_file % dimids_scl(1)
	3355	var % dimids(2) = output_file % dimids_vel(2)
	3356	var % dimids(3) = output_file % dimids_scl(3)
	3357	var % dimvarids(1) = output_file % dimvarids_scl(1)
	3358	var % dimvarids(2) = output_file % dimvarids_vel(2)
	3359	var % dimvarids(3) = output_file % dimvarids_scl(3)
	3360	var % to_be_processed = init_variables_required
[3395]	3361	var % is_internal = .FALSE.
[2696]	3362	var % task = "interpolate_3d"
	3363
	3364	CASE( 'init w' )
	3365	var % nt = 1
	3366	var % lod = 2
	3367	var % ndim = 3
	3368	var % dimids(1) = output_file % dimids_scl(1)
	3369	var % dimids(2) = output_file % dimids_scl(2)
	3370	var % dimids(3) = output_file % dimids_vel(3)
	3371	var % dimvarids(1) = output_file % dimvarids_scl(1)
	3372	var % dimvarids(2) = output_file % dimvarids_scl(2)
	3373	var % dimvarids(3) = output_file % dimvarids_vel(3)
	3374	var % to_be_processed = init_variables_required
[3395]	3375	var % is_internal = .FALSE.
[2696]	3376	var % task = "interpolate_3d"
	3377
	3378	CASE( 'init scalar profile', 'init u profile', 'init v profile')
	3379	var % nt = 1
	3380	var % lod = 1
	3381	var % ndim = 1
	3382	var % dimids(1) = output_file % dimids_scl(3) !z
	3383	var % dimvarids(1) = output_file % dimvarids_scl(3) !z
	3384	var % to_be_processed = init_variables_required
[3395]	3385	var % is_internal = .FALSE.
[2696]	3386	var % task = "average profile"
	3387
	3388	CASE( 'init w profile')
	3389	var % nt = 1
	3390	var % lod = 1
	3391	var % ndim = 1
	3392	var % dimids(1) = output_file % dimids_vel(3) !z
	3393	var % dimvarids(1) = output_file % dimvarids_vel(3) !z
	3394	var % to_be_processed = init_variables_required
[3395]	3395	var % is_internal = .FALSE.
[2696]	3396	var % task = "average profile"
	3397
[3182]	3398	CASE( 'surface forcing' )
	3399	var % lod = -1
[2696]	3400	var % ndim = 3
	3401	var % dimids(3) = output_file % dimid_time
	3402	var % dimids(1:2) = output_file % dimids_soil(1:2)
	3403	var % dimvarids(3) = output_file % dimvarid_time
	3404	var % dimvarids(1:2) = output_file % dimvarids_soil(1:2)
[3456]	3405	var % to_be_processed = surface_forcing_required
[3395]	3406	var % is_internal = .FALSE.
[2696]	3407	var % task = "interpolate_2d"
	3408
[3182]	3409	CASE( 'left scalar', 'right scalar') ! same as right
	3410	var % lod = -1
[2696]	3411	var % ndim = 3
	3412	var % dimids(3) = output_file % dimid_time
	3413	var % dimids(1) = output_file % dimids_scl(2)
	3414	var % dimids(2) = output_file % dimids_scl(3)
	3415	var % dimvarids(3) = output_file % dimvarid_time
	3416	var % dimvarids(1) = output_file % dimvarids_scl(2)
	3417	var % dimvarids(2) = output_file % dimvarids_scl(3)
	3418	var % to_be_processed = boundary_variables_required
[3395]	3419	var % is_internal = .FALSE.
[2696]	3420	var % task = "interpolate_3d"
	3421
[3182]	3422	CASE( 'north scalar', 'south scalar') ! same as south
	3423	var % lod = -1
[2696]	3424	var % ndim = 3
	3425	var % dimids(3) = output_file % dimid_time
	3426	var % dimids(1) = output_file % dimids_scl(1)
	3427	var % dimids(2) = output_file % dimids_scl(3)
	3428	var % dimvarids(3) = output_file % dimvarid_time
	3429	var % dimvarids(1) = output_file % dimvarids_scl(1)
	3430	var % dimvarids(2) = output_file % dimvarids_scl(3)
	3431	var % to_be_processed = boundary_variables_required
[3395]	3432	var % is_internal = .FALSE.
[2696]	3433	var % task = "interpolate_3d"
	3434
[3182]	3435	CASE( 'top scalar', 'top w' )
	3436	var % lod = -1
[2696]	3437	var % ndim = 3
	3438	var % dimids(3) = output_file % dimid_time
	3439	var % dimids(1) = output_file % dimids_scl(1)
	3440	var % dimids(2) = output_file % dimids_scl(2)
	3441	var % dimvarids(3) = output_file % dimvarid_time
	3442	var % dimvarids(1) = output_file % dimvarids_scl(1)
	3443	var % dimvarids(2) = output_file % dimvarids_scl(2)
	3444	var % to_be_processed = boundary_variables_required
[3395]	3445	var % is_internal = .FALSE.
[2696]	3446	var % task = "interpolate_3d"
	3447
[3182]	3448	CASE( 'left u', 'right u' )
	3449	var % lod = -1
[2696]	3450	var % ndim = 3
	3451	var % dimids(3) = output_file % dimid_time
	3452	var % dimids(1) = output_file % dimids_scl(2)
	3453	var % dimids(2) = output_file % dimids_scl(3)
	3454	var % dimvarids(3) = output_file % dimvarid_time
	3455	var % dimvarids(1) = output_file % dimvarids_scl(2)
	3456	var % dimvarids(2) = output_file % dimvarids_scl(3)
	3457	var % to_be_processed = boundary_variables_required
[3395]	3458	var % is_internal = .FALSE.
[2696]	3459	var % task = "interpolate_3d"
	3460
[3182]	3461	CASE( 'north u', 'south u' )
	3462	var % lod = -1
[2696]	3463	var % ndim = 3
	3464	var % dimids(3) = output_file % dimid_time !t
	3465	var % dimids(1) = output_file % dimids_vel(1) !x
	3466	var % dimids(2) = output_file % dimids_scl(3) !z
	3467	var % dimvarids(3) = output_file % dimvarid_time
	3468	var % dimvarids(1) = output_file % dimvarids_vel(1)
	3469	var % dimvarids(2) = output_file % dimvarids_scl(3)
	3470	var % to_be_processed = boundary_variables_required
[3395]	3471	var % is_internal = .FALSE.
[2696]	3472	var % task = "interpolate_3d"
	3473
[3182]	3474	CASE( 'top u' )
	3475	var % lod = -1
[2696]	3476	var % ndim = 3
	3477	var % dimids(3) = output_file % dimid_time !t
	3478	var % dimids(1) = output_file % dimids_vel(1) !x
	3479	var % dimids(2) = output_file % dimids_scl(2) !z
	3480	var % dimvarids(3) = output_file % dimvarid_time
	3481	var % dimvarids(1) = output_file % dimvarids_vel(1)
	3482	var % dimvarids(2) = output_file % dimvarids_scl(2)
	3483	var % to_be_processed = boundary_variables_required
[3395]	3484	var % is_internal = .FALSE.
[2696]	3485	var % task = "interpolate_3d"
	3486
[3182]	3487	CASE( 'left v', 'right v' )
	3488	var % lod = -1
[2696]	3489	var % ndim = 3
	3490	var % dimids(3) = output_file % dimid_time
	3491	var % dimids(1) = output_file % dimids_vel(2)
	3492	var % dimids(2) = output_file % dimids_scl(3)
	3493	var % dimvarids(3) = output_file % dimvarid_time
	3494	var % dimvarids(1) = output_file % dimvarids_vel(2)
	3495	var % dimvarids(2) = output_file % dimvarids_scl(3)
	3496	var % to_be_processed = boundary_variables_required
[3395]	3497	var % is_internal = .FALSE.
[2696]	3498	var % task = "interpolate_3d"
	3499
[3182]	3500	CASE( 'north v', 'south v' )
	3501	var % lod = -1
[2696]	3502	var % ndim = 3
	3503	var % dimids(3) = output_file % dimid_time !t
	3504	var % dimids(1) = output_file % dimids_scl(1) !x
	3505	var % dimids(2) = output_file % dimids_scl(3) !z
	3506	var % dimvarids(3) = output_file % dimvarid_time
	3507	var % dimvarids(1) = output_file % dimvarids_scl(1)
	3508	var % dimvarids(2) = output_file % dimvarids_scl(3)
	3509	var % to_be_processed = boundary_variables_required
[3395]	3510	var % is_internal = .FALSE.
[2696]	3511	var % task = "interpolate_3d"
	3512
[3182]	3513	CASE( 'top v' )
	3514	var % lod = -1
[2696]	3515	var % ndim = 3
	3516	var % dimids(3) = output_file % dimid_time !t
	3517	var % dimids(1) = output_file % dimids_scl(1) !x
	3518	var % dimids(2) = output_file % dimids_vel(2) !z
	3519	var % dimvarids(3) = output_file % dimvarid_time
	3520	var % dimvarids(1) = output_file % dimvarids_scl(1)
	3521	var % dimvarids(2) = output_file % dimvarids_vel(2)
	3522	var % to_be_processed = boundary_variables_required
[3395]	3523	var % is_internal = .FALSE.
[2696]	3524	var % task = "interpolate_3d"
	3525
[3182]	3526	CASE( 'left w', 'right w')
	3527	var % lod = -1
[2696]	3528	var % ndim = 3
	3529	var % dimids(3) = output_file % dimid_time
	3530	var % dimids(1) = output_file % dimids_scl(2)
	3531	var % dimids(2) = output_file % dimids_vel(3)
	3532	var % dimvarids(3) = output_file % dimvarid_time
	3533	var % dimvarids(1) = output_file % dimvarids_scl(2)
	3534	var % dimvarids(2) = output_file % dimvarids_vel(3)
	3535	var % to_be_processed = boundary_variables_required
[3395]	3536	var % is_internal = .FALSE.
[2696]	3537	var % task = "interpolate_3d"
	3538
[3182]	3539	CASE( 'north w', 'south w' )
	3540	var % lod = -1
[2696]	3541	var % ndim = 3
	3542	var % dimids(3) = output_file % dimid_time !t
	3543	var % dimids(1) = output_file % dimids_scl(1) !x
	3544	var % dimids(2) = output_file % dimids_vel(3) !z
	3545	var % dimvarids(3) = output_file % dimvarid_time
	3546	var % dimvarids(1) = output_file % dimvarids_scl(1)
	3547	var % dimvarids(2) = output_file % dimvarids_vel(3)
	3548	var % to_be_processed = boundary_variables_required
[3395]	3549	var % is_internal = .FALSE.
[2696]	3550	var % task = "interpolate_3d"
	3551
[3182]	3552	CASE( 'time series' )
[2696]	3553	var % lod = 0
	3554	var % ndim = 1
	3555	var % dimids(1) = output_file % dimid_time !t
	3556	var % dimvarids(1) = output_file % dimvarid_time
	3557	var % to_be_processed = .TRUE.
[3395]	3558	var % is_internal = .FALSE.
	3559	var % task = "average profile"
[2696]	3560
[3182]	3561	CASE( 'constant scalar profile' )
	3562	var % lod = -1
[2696]	3563	var % ndim = 2
	3564	var % dimids(2) = output_file % dimid_time !t
	3565	var % dimids(1) = output_file % dimids_scl(3) !z
	3566	var % dimvarids(2) = output_file % dimvarid_time
	3567	var % dimvarids(1) = output_file % dimvarids_scl(3)
	3568	var % to_be_processed = .TRUE.
[3395]	3569	var % is_internal = .FALSE.
[3182]	3570	var % task = "set profile"
[2696]	3571
[3182]	3572	CASE( 'large-scale scalar forcing' )
	3573	var % lod = -1
	3574	var % ndim = 2
	3575	var % dimids(2) = output_file % dimid_time !t
	3576	var % dimids(1) = output_file % dimids_scl(3) !z
	3577	var % dimvarids(2) = output_file % dimvarid_time
	3578	var % dimvarids(1) = output_file % dimvarids_scl(3)
	3579	var % to_be_processed = ls_forcing_variables_required
[3395]	3580	var % is_internal = .FALSE.
[3182]	3581	var % task = "average large-scale profile"
	3582
[3395]	3583	CASE( 'geostrophic' )
	3584	var % lod = -1
	3585	var % ndim = 2
	3586	var % dimids(2) = output_file % dimid_time !t
	3587	var % dimids(1) = output_file % dimids_scl(3) !z
	3588	var % dimvarids(2) = output_file % dimvarid_time
	3589	var % dimvarids(1) = output_file % dimvarids_scl(3)
	3590	var % to_be_processed = .TRUE.
	3591	var % is_internal = .FALSE.
	3592	var % task = "geostrophic winds"
	3593
[3182]	3594	CASE( 'large-scale w forcing' )
	3595	var % lod = -1
	3596	var % ndim = 2
	3597	var % dimids(2) = output_file % dimid_time !t
	3598	var % dimids(1) = output_file % dimids_vel(3) !z
	3599	var % dimvarids(2) = output_file % dimvarid_time
	3600	var % dimvarids(1) = output_file % dimvarids_vel(3)
	3601	var % to_be_processed = ls_forcing_variables_required
[3395]	3602	var % is_internal = .FALSE.
[3182]	3603	var % task = "average large-scale profile"
	3604
[3395]	3605	CASE( 'internal profile' )
	3606	var % lod = -1
	3607	var % ndim = 2
	3608	var % dimids(2) = output_file % dimid_time !t
	3609	var % dimids(1) = output_file % dimids_scl(3) !z
	3610	var % dimvarids(2) = output_file % dimvarid_time
	3611	var % dimvarids(1) = output_file % dimvarids_scl(3)
	3612	var % to_be_processed = .TRUE.
	3613	var % is_internal = .TRUE.
	3614	var % task = "internal profile"
	3615
[2696]	3616	CASE DEFAULT
	3617	message = "Variable kind '" // TRIM(kind) // "' not recognized."
[3615]	3618	CALL inifor_abort ('init_nc_var', message)
[2696]	3619
	3620	END SELECT
	3621
	3622	END FUNCTION init_nc_var
	3623
	3624
	3625	SUBROUTINE fini_variables()
	3626
[3395]	3627	CALL report('fini_variables', 'Deallocating variable table', cfg % debug)
[2696]	3628	DEALLOCATE( input_var_table )
	3629
	3630	END SUBROUTINE fini_variables
	3631
	3632
	3633	SUBROUTINE fini_io_groups()
	3634
[3395]	3635	CALL report('fini_io_groups', 'Deallocating IO groups', cfg % debug)
[2696]	3636	DEALLOCATE( io_group_list )
	3637
	3638	END SUBROUTINE fini_io_groups
	3639
	3640
	3641	SUBROUTINE fini_file_lists()
	3642
[3395]	3643	CALL report('fini_file_lists', 'Deallocating file lists', cfg % debug)
[2696]	3644	DEALLOCATE( flow_files, soil_files, radiation_files, soil_moisture_files )
	3645
	3646	END SUBROUTINE fini_file_lists
	3647
	3648
	3649	!------------------------------------------------------------------------------!
	3650	! Description:
	3651	! ------------
	3652	!> Carries out any physical conversion of the quantities in the given input
	3653	!> buffer needed to obtain the quantity required by PALM-4U. For instance,
	3654	!> velocities are rotated to the PALM-4U coordinate system and the potential
	3655	!> temperature is computed from the absolute temperature and pressure.
	3656	!>
	3657	!> Note, that the preprocessing does not include any grid change. The result
	3658	!> array will match a COSMO-DE scalar array.
	3659	!------------------------------------------------------------------------------!
	3660	SUBROUTINE preprocess(group, input_buffer, cosmo_grid, iter)
	3661
	3662	TYPE(io_group), INTENT(INOUT), TARGET :: group
	3663	TYPE(container), INTENT(INOUT), ALLOCATABLE :: input_buffer(:)
	3664	TYPE(grid_definition), INTENT(IN) :: cosmo_grid
	3665	INTEGER, INTENT(IN) :: iter
	3666
	3667	REAL(dp), ALLOCATABLE :: basic_state_pressure(:)
[3182]	3668	TYPE(container), ALLOCATABLE :: preprocess_buffer(:)
[2696]	3669	INTEGER :: hour, dt
	3670	INTEGER :: i, j, k
	3671	INTEGER :: nx, ny, nz
	3672
	3673	input_buffer(:) % is_preprocessed = .FALSE.
	3674
	3675	SELECT CASE( group % kind )
	3676
	3677	CASE( 'velocities' )
[3557]	3678	!
	3679	!-- Allocate a compute buffer with the same number of arrays as the input
[3182]	3680	ALLOCATE( preprocess_buffer( SIZE(input_buffer) ) )
[2696]	3681
[3557]	3682	!
	3683	!-- Allocate u and v arrays with scalar dimensions
[2696]	3684	nx = SIZE(input_buffer(1) % array, 1)
	3685	ny = SIZE(input_buffer(1) % array, 2)
	3686	nz = SIZE(input_buffer(1) % array, 3)
[3182]	3687	ALLOCATE( preprocess_buffer(1) % array(nx, ny, nz) ) ! u buffer
	3688	ALLOCATE( preprocess_buffer(2) % array(nx, ny, nz) ) ! v buffer
[2696]	3689
	3690	CALL run_control('time', 'alloc')
	3691
[3557]	3692	!
	3693	!-- interpolate U and V to centres
[2696]	3694	CALL centre_velocities( u_face = input_buffer(1) % array, &
	3695	v_face = input_buffer(2) % array, &
[3182]	3696	u_centre = preprocess_buffer(1) % array, &
	3697	v_centre = preprocess_buffer(2) % array )
[2696]	3698
[3182]	3699	cfg % rotation_method = 'rotated-pole'
	3700	SELECT CASE(cfg % rotation_method)
	3701
	3702	CASE('rotated-pole')
[3557]	3703	!
	3704	!-- rotate U and V to PALM-4U orientation and overwrite U and V with
	3705	!-- rotated velocities
[3182]	3706	DO k = 1, nz
[3395]	3707	DO j = 1, ny
	3708	DO i = 1, nx
[3182]	3709	CALL uv2uvrot( urot = preprocess_buffer(1) % array(i,j,k), &
	3710	vrot = preprocess_buffer(2) % array(i,j,k), &
	3711	rlat = cosmo_grid % lat(j-1), &
	3712	rlon = cosmo_grid % lon(i-1), &
	3713	pollat = phi_cn, &
	3714	pollon = lambda_cn, &
	3715	u = input_buffer(1) % array(i,j,k), &
	3716	v = input_buffer(2) % array(i,j,k) )
	3717	END DO
	3718	END DO
	3719	END DO
	3720
	3721	CASE DEFAULT
	3722	message = "Rotation method '" // TRIM(cfg % rotation_method) // &
	3723	"' not recognized."
[3615]	3724	CALL inifor_abort('preprocess', message)
[3182]	3725
	3726	END SELECT
	3727
[2696]	3728	input_buffer(1) % array(1,:,:) = 0.0_dp
	3729	input_buffer(2) % array(1,:,:) = 0.0_dp
	3730	input_buffer(1) % array(:,1,:) = 0.0_dp
	3731	input_buffer(2) % array(:,1,:) = 0.0_dp
	3732
	3733	input_buffer(1:2) % is_preprocessed = .TRUE.
	3734	CALL run_control('time', 'comp')
	3735
[3182]	3736	DEALLOCATE( preprocess_buffer )
[2696]	3737	CALL run_control('time', 'alloc')
	3738
	3739	message = "Input buffers for group '" // TRIM(group % kind) // "'"//&
	3740	" preprocessed sucessfully."
	3741	CALL report('preprocess', message)
	3742
[3395]	3743	CASE( 'thermodynamics' ) ! T, P, QV
[2696]	3744	nx = SIZE(input_buffer(1) % array, 1)
	3745	ny = SIZE(input_buffer(1) % array, 2)
	3746	nz = SIZE(input_buffer(1) % array, 3)
	3747
[3557]	3748	!
	3749	!-- Compute absolute pressure if presure perturbation has been read in.
[2696]	3750	IF ( TRIM(group % in_var_list(2) % name) == 'PP' ) THEN
	3751	message = "Absolute pressure, P, not available, " // &
	3752	"computing from pressure preturbation PP."
	3753	CALL report('preprocess', message)
	3754
	3755	ALLOCATE( basic_state_pressure(1:nz) )
	3756	CALL run_control('time', 'alloc')
	3757
	3758	DO j = 1, ny
	3759	DO i = 1, nx
	3760
	3761	CALL get_basic_state(cosmo_grid % hfl(i,j,:), BETA, P_SL, T_SL,&
	3762	RD, G, basic_state_pressure)
	3763
[3557]	3764	!
	3765	!-- Overwrite pressure perturbation with absolute pressure. HECTO
	3766	!-- converts pressure perturbation from hPa to Pa.
[2696]	3767	input_buffer (2) % array(i,j,:) = &
[3395]	3768	HECTO * input_buffer (2) % array(i,j,:) + &
	3769	basic_state_pressure(:)
[2696]	3770
	3771	END DO
	3772	END DO
	3773	CALL run_control('time', 'comp')
	3774
	3775	DEALLOCATE( basic_state_pressure )
	3776	CALL run_control('time', 'alloc')
	3777
	3778	group % in_var_list(2) % name = 'P'
	3779
	3780	END IF
[3557]	3781	!
	3782	!-- mark pressure as preprocessed
[3395]	3783	input_buffer(2) % is_preprocessed = .TRUE.
[2696]	3784
[3557]	3785	!
	3786	!-- Copy temperature to the last input buffer array
[3395]	3787	ALLOCATE( &
	3788	input_buffer( group % n_output_quantities ) % array (nx, ny, nz) &
	3789	)
	3790	input_buffer(group % n_output_quantities) % array(:,:,:) = &
	3791	input_buffer(1) % array(:,:,:)
	3792
[3557]	3793	!
	3794	!-- Convert absolute in place to potential temperature
[3395]	3795	CALL potential_temperature( &
	3796	t = input_buffer(1) % array(:,:,:), &
	3797	p = input_buffer(2) % array(:,:,:), &
	3798	p_ref = P_REF, &
	3799	r = RD_PALM, &
	3800	cp = CP_PALM &
	3801	)
[2696]	3802
[3557]	3803	!
	3804	!-- mark potential temperature as preprocessed
[3395]	3805	input_buffer(1) % is_preprocessed = .TRUE.
	3806
[3557]	3807	!
	3808	!-- Convert temperature copy to density
[3395]	3809	CALL moist_density( &
	3810	t_rho = input_buffer(group % n_output_quantities) % array(:,:,:), &
	3811	p = input_buffer(2) % array(:,:,:), &
	3812	qv = input_buffer(3) % array(:,:,:), &
	3813	rd = RD, &
	3814	rv = RV &
	3815	)
	3816
[3557]	3817	!
	3818	!-- mark qv as preprocessed
[3395]	3819	input_buffer(3) % is_preprocessed = .TRUE.
	3820
[3557]	3821	!
	3822	!-- mark density as preprocessed
[3395]	3823	input_buffer(group % n_output_quantities) % is_preprocessed = .TRUE.
	3824
	3825
[2696]	3826	message = "Input buffers for group '" // TRIM(group % kind) // "'"//&
	3827	" preprocessed sucessfully."
	3828	CALL report('preprocess', message)
	3829	CALL run_control('time', 'comp')
	3830
	3831	CASE( 'scalar' ) ! S or W
	3832	input_buffer(:) % is_preprocessed = .TRUE.
	3833	CALL run_control('time', 'comp')
	3834
	3835	CASE( 'soil-temperature' ) !
	3836
	3837	CALL fill_water_cells(soiltyp, input_buffer(1) % array, &
	3838	SIZE(input_buffer(1) % array, 3), &
	3839	FILL_ITERATIONS)
	3840	input_buffer(:) % is_preprocessed = .TRUE.
	3841	CALL run_control('time', 'comp')
	3842
	3843	CASE( 'soil-water' ) !
	3844
	3845	CALL fill_water_cells(soiltyp, input_buffer(1) % array, &
	3846	SIZE(input_buffer(1) % array, 3), &
	3847	FILL_ITERATIONS)
	3848
	3849	nx = SIZE(input_buffer(1) % array, 1)
	3850	ny = SIZE(input_buffer(1) % array, 2)
	3851	nz = SIZE(input_buffer(1) % array, 3)
	3852
	3853	DO k = 1, nz
	3854	DO j = 1, ny
	3855	DO i = 1, nx
	3856	input_buffer(1) % array(i,j,k) = &
	3857	input_buffer(1) % array(i,j,k) * d_depth_rho_inv(k)
	3858	END DO
	3859	END DO
	3860	END DO
	3861
	3862	message = "Converted soil water from [kg/m^2] to [m^3/m^3]"
	3863	CALL report('preprocess', message)
	3864
	3865	input_buffer(:) % is_preprocessed = .TRUE.
	3866	CALL run_control('time', 'comp')
	3867
	3868	CASE( 'surface' ) !
	3869	input_buffer(:) % is_preprocessed = .TRUE.
	3870	CALL run_control('time', 'comp')
	3871
	3872	CASE( 'accumulated' ) !
	3873	message = "De-accumulating '" // TRIM(group % in_var_list(1) % name) //&
	3874	"' in iteration " // TRIM(str(iter))
	3875	CALL report('preprocess', message)
	3876
	3877	hour = iter - 1
	3878	dt = MODULO(hour, 3) + 1 ! averaging period
	3879	SELECT CASE(dt)
	3880
[3557]	3881	!
	3882	!-- input has been accumulated over one hour. Leave as is
	3883	!-- input_buffer(1) % array(:,:,:) carrries one-hour integral
[2696]	3884	CASE(1)
	3885
[3557]	3886	!
	3887	!-- input has been accumulated over two hours. Subtract previous step
	3888	!-- input_buffer(1) % array(:,:,:) carrries one-hour integral
	3889	!-- input_buffer(2) % array(:,:,:) carrries two-hour integral
[2696]	3890	CASE(2)
	3891	CALL deaverage( &
	3892	avg_1 = input_buffer(1) % array(:,:,:), t1 = 1.0_dp, &
	3893	avg_2 = input_buffer(2) % array(:,:,:), t2 = 1.0_dp, &
	3894	avg_3 = input_buffer(1) % array(:,:,:), t3 = 1.0_dp )
[3557]	3895	!
	3896	!-- input_buffer(1) % array(:,:,:) carrries one-hour integral of second hour
[2696]	3897
[3557]	3898	!
	3899	!-- input has been accumulated over three hours. Subtract previous step
	3900	!-- input_buffer(1) % array(:,:,:) carrries three-hour integral
	3901	!-- input_buffer(2) % array(:,:,:) still carrries two-hour integral
[2696]	3902	CASE(3)
	3903	CALL deaverage( &
	3904	avg_1 = input_buffer(2) % array(:,:,:), t1 = 1.0_dp, &
	3905	avg_2 = input_buffer(1) % array(:,:,:), t2 = 1.0_dp, &
	3906	avg_3 = input_buffer(1) % array(:,:,:), t3 = 1.0_dp )
[3557]	3907	!
	3908	!-- input_buffer(1) % array(:,:,:) carrries one-hour integral of third hourA
[2696]	3909
	3910	CASE DEFAULT
	3911	message = "Invalid averaging period '" // TRIM(str(dt)) // " hours"
[3615]	3912	CALL inifor_abort('preprocess', message)
[2696]	3913
	3914	END SELECT
	3915	input_buffer(:) % is_preprocessed = .TRUE.
	3916	CALL run_control('time', 'comp')
	3917
	3918	CASE( 'running average' ) !
	3919	message = "De-averaging '" // TRIM(group % in_var_list(1) % name) // &
	3920	"' in iteration " // TRIM(str(iter))
	3921	CALL report('preprocess', message)
	3922
	3923	hour = iter - 1
[3557]	3924	!
	3925	!-- averaging period
	3926	dt = MODULO(hour, 3) + 1
[2696]	3927	SELECT CASE(dt)
[3557]	3928	!
	3929	!-- input has been accumulated over one hour. Leave as is
	3930	!-- input_buffer(1) % array(:,:,:) carrries one-hour integral
[2696]	3931	CASE(1)
	3932
[3557]	3933	!
	3934	!-- input has been accumulated over two hours. Subtract previous step
	3935	!-- input_buffer(1) % array(:,:,:) carrries one-hour integral
	3936	!-- input_buffer(2) % array(:,:,:) carrries two-hour integral
[2696]	3937	CASE(2)
	3938	CALL deaverage( input_buffer(1) % array(:,:,:), 1.0_dp, &
	3939	input_buffer(2) % array(:,:,:), 2.0_dp, &
	3940	input_buffer(1) % array(:,:,:), 1.0_dp)
[3557]	3941	!
	3942	!-- input_buffer(1) % array(:,:,:) carrries one-hour integral of second hour
[2696]	3943
[3557]	3944	!
	3945	!-- input has been accumulated over three hours. Subtract previous step
	3946	!-- input_buffer(1) % array(:,:,:) carrries three-hour integral
	3947	!-- input_buffer(2) % array(:,:,:) still carrries two-hour integral
[2696]	3948	CASE(3)
	3949	CALL deaverage( input_buffer(2) % array(:,:,:), 2.0_dp, &
	3950	input_buffer(1) % array(:,:,:), 3.0_dp, &
	3951	input_buffer(1) % array(:,:,:), 1.0_dp)
[3557]	3952	!
	3953	!-- input_buffer(1) % array(:,:,:) carrries one-hour integral of third hourA
[2696]	3954
	3955	CASE DEFAULT
	3956	message = "Invalid averaging period '" // TRIM(str(dt)) // " hours"
[3615]	3957	CALL inifor_abort('preprocess', message)
[2696]	3958
	3959	END SELECT
	3960	input_buffer(:) % is_preprocessed = .TRUE.
	3961
	3962	CASE DEFAULT
[3395]	3963	message = "IO group kind '" // TRIM(group % kind) // "' is not supported."
[3615]	3964	CALL inifor_abort('prerpocess', message)
[2696]	3965
	3966	END SELECT
	3967	CALL run_control('time', 'comp')
	3968
	3969	END SUBROUTINE preprocess
	3970
	3971
[3557]	3972	!------------------------------------------------------------------------------!
[2696]	3973	! Description:
	3974	! ------------
	3975	!> Computes average soil values in COSMO-DE water cells from neighbouring
	3976	!> non-water cells. This is done as a preprocessing step for the COSMO-DE
	3977	!> soil input arrays, which contain unphysical values for water cells.
	3978	!>
	3979	!> This routine computes the average of up to all nine neighbouring cells
	3980	!> or keeps the original value, if not at least one non-water neightbour
	3981	!> is available.
	3982	!>
	3983	!> By repeatedly applying this step, soil data can be extrapolated into
	3984	!> 'water' regions occupying multiple cells, one cell per iteration.
	3985	!>
	3986	!> Input parameters:
	3987	!> -----------------
	3988	!> soiltyp : 2d map of COSMO-DE soil types
	3989	!> nz : number of layers in the COSMO-DE soil
	3990	!> niter : number iterations
	3991	!>
	3992	!> Output parameters:
	3993	!> ------------------
	3994	!> array : the soil array (i.e. water content or temperature)
	3995	!------------------------------------------------------------------------------!
	3996	SUBROUTINE fill_water_cells(soiltyp, array, nz, niter)
	3997	INTEGER(hp), DIMENSION(:,:,:), INTENT(IN) :: soiltyp
	3998	REAL(dp), DIMENSION(:,:,:), INTENT(INOUT) :: array
	3999	INTEGER, INTENT(IN) :: nz, niter
	4000
	4001	REAL(dp), DIMENSION(nz) :: column
	4002	INTEGER(hp), DIMENSION(:,:), ALLOCATABLE :: old_soiltyp, new_soiltyp
	4003	INTEGER :: l, i, j, nx, ny, n_cells, ii, jj, iter
	4004	INTEGER, DIMENSION(8) :: di, dj
	4005
	4006	nx = SIZE(array, 1)
	4007	ny = SIZE(array, 2)
	4008	di = (/ -1, -1, -1, 0, 0, 1, 1, 1 /)
	4009	dj = (/ -1, 0, 1, -1, 1, -1, 0, 1 /)
	4010
	4011	ALLOCATE(old_soiltyp(SIZE(soiltyp,1), &
	4012	SIZE(soiltyp,2) ))
	4013
	4014	ALLOCATE(new_soiltyp(SIZE(soiltyp,1), &
	4015	SIZE(soiltyp,2) ))
	4016
	4017	old_soiltyp(:,:) = soiltyp(:,:,1)
	4018	new_soiltyp(:,:) = soiltyp(:,:,1)
	4019
	4020	DO iter = 1, niter
	4021
	4022	DO j = 1, ny
	4023	DO i = 1, nx
	4024
	4025	IF (old_soiltyp(i,j) == WATER_ID) THEN
	4026
	4027	n_cells = 0
	4028	column(:) = 0.0_dp
	4029	DO l = 1, SIZE(di)
	4030
	4031	ii = MIN(nx, MAX(1, i + di(l)))
	4032	jj = MIN(ny, MAX(1, j + dj(l)))
	4033
	4034	IF (old_soiltyp(ii,jj) .NE. WATER_ID) THEN
	4035	n_cells = n_cells + 1
	4036	column(:) = column(:) + array(ii,jj,:)
	4037	END IF
	4038
	4039	END DO
	4040
[3557]	4041	!
	4042	!-- Overwrite if at least one non-water neighbour cell is available
[2696]	4043	IF (n_cells > 0) THEN
	4044	array(i,j,:) = column(:) / n_cells
	4045	new_soiltyp(i,j) = 0
	4046	END IF
	4047
	4048	END IF
	4049
	4050	END DO
	4051	END DO
	4052
	4053	old_soiltyp(:,:) = new_soiltyp(:,:)
	4054
	4055	END DO
	4056
	4057	DEALLOCATE(old_soiltyp, new_soiltyp)
	4058
	4059	END SUBROUTINE fill_water_cells
	4060
[3537]	4061
	4062	!------------------------------------------------------------------------------!
	4063	! Description:
	4064	! ------------
	4065	!> Fills the thickness array of the COSMO soil layers. Since COSMO's (i.e.
	4066	!> TERRA_ML's [1]) soil layer boundaries follow the rule
	4067	!>
	4068	!> depth(0) = 0.0, and
	4069	!> depth(k) = 0.01 * 3**(k-1), k in [1,2,3,...,7]
	4070	!>
	4071	!> and full levels are defined as the midpoints between two layer boundaries,
	4072	!> all except the first layer thicknesses equal the depth of the midpoint.
	4073	!>
	4074	!> [1] A Description of the Nonhydrostatic Regional COSMO Model Part II :
	4075	!> Physical Parameterization*, Sect. 11 TERRA_ML.
	4076	!> http://www.cosmo-model.org/content/model/documentation/core/cosmoPhysParamtr.pdf)
	4077	!>
	4078	!> Input parameters:
	4079	!> -----------------
	4080	!>
	4081	!> depths: array of full soil layer depths (cell midpoints)
	4082	!>
	4083	!>
	4084	!> Output parameters:
	4085	!> ------------------
	4086	!>
	4087	!> d_depth: array of soil layer thicknesses
	4088	!>
	4089	!------------------------------------------------------------------------------!
	4090	SUBROUTINE get_soil_layer_thickness(depths, d_depth)
	4091
	4092	REAL(dp), INTENT(IN) :: depths(:)
	4093	REAL(dp), INTENT(OUT) :: d_depth(:)
	4094
	4095	d_depth(:) = depths(:)
	4096	d_depth(1) = 2.0_dp * depths(1)
	4097
	4098	END SUBROUTINE get_soil_layer_thickness
	4099
[3618]	4100	END MODULE inifor_grid
[3680]	4101	#endif

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |