Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |