Home

Context Navigation

source: palm/trunk/UTIL/inifor/src/grid.f90 @ 3201

Last change on this file since 3201 was 3183, checked in by suehring, 6 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 178.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |