Home

Context Navigation

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

Last change on this file since 3447 was 3447, checked in by eckhard, 5 years ago
inifor: Renamed source files for compatibilty with PALM build system
Property svn:keywords set to `Id`
File size: 204.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |