Home

Context Navigation

source: palm/trunk/SOURCE/radiation_model_mod.f90 @ 2008

Last change on this file since 2008 was 2008, checked in by kanani, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 136.1 KB

Rev	Line
[1826]	1	!> @file radiation_model_mod.f90
[2000]	2	!------------------------------------------------------------------------------!
[1496]	3	! This file is part of PALM.
	4	!
[2000]	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
	8	! version.
[1496]	9	!
	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.
	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	!
[1818]	17	! Copyright 1997-2016 Leibniz Universitaet Hannover
[2000]	18	!------------------------------------------------------------------------------!
[1496]	19	!
	20	! Current revisions:
	21	! -----------------
[2008]	22	!
	23	!
	24	! Former revisions:
	25	! -----------------
	26	! $Id: radiation_model_mod.f90 2008 2016-08-24 15:53:00Z kanani $
	27	!
	28	! 2007 2016-08-24 15:47:17Z kanani
[2007]	29	! Added calculation of solar directional vector for new urban surface
	30	! model,
	31	! accounted for urban_surface model in radiation_check_parameters,
	32	! correction of comments for zenith angle.
[1977]	33	!
[2001]	34	! 2000 2016-08-20 18:09:15Z knoop
	35	! Forced header and separation lines into 80 columns
	36	!
[1977]	37	! 1976 2016-07-27 13:28:04Z maronga
[1976]	38	! Output of 2D/3D/masked data is now directly done within this module. The
	39	! radiation schemes have been simplified for better usability so that
	40	! rad_lw_in, rad_lw_out, rad_sw_in, and rad_sw_out are available independent of
	41	! the radiation code used.
[1854]	42	!
[1857]	43	! 1856 2016-04-13 12:56:17Z maronga
	44	! Bugfix: allocation of rad_lw_out for radiation_scheme = 'clear-sky'
	45	!
[1854]	46	! 1853 2016-04-11 09:00:35Z maronga
	47	! Added routine for radiation_scheme = constant.
	48	!
[1852]	49	! 1849 2016-04-08 11:33:18Z hoffmann
[1851]	50	! Adapted for modularization of microphysics
[1852]	51	!
[1827]	52	! 1826 2016-04-07 12:01:39Z maronga
	53	! Further modularization.
	54	!
[1789]	55	! 1788 2016-03-10 11:01:04Z maronga
	56	! Added new albedo class for pavements / roads.
	57	!
[1784]	58	! 1783 2016-03-06 18:36:17Z raasch
	59	! palm-netcdf-module removed in order to avoid a circular module dependency,
	60	! netcdf-variables moved to netcdf-module, new routine netcdf_handle_error_rad
	61	! added
	62	!
[1758]	63	! 1757 2016-02-22 15:49:32Z maronga
	64	! Added parameter unscheduled_radiation_calls. Bugfix: interpolation of sounding
	65	! profiles for pressure and temperature above the LES domain.
	66	!
[1710]	67	! 1709 2015-11-04 14:47:01Z maronga
	68	! Bugfix: set initial value for rrtm_lwuflx_dt to zero, small formatting
	69	! corrections
	70	!
[1702]	71	! 1701 2015-11-02 07:43:04Z maronga
	72	! Bugfixes: wrong index for output of timeseries, setting of nz_snd_end
	73	!
[1692]	74	! 1691 2015-10-26 16:17:44Z maronga
	75	! Added option for spin-up runs without radiation (skip_time_do_radiation). Bugfix
	76	! in calculation of pressure profiles. Bugfix in calculation of trace gas profiles.
	77	! Added output of radiative heating rates.
	78	!
[1683]	79	! 1682 2015-10-07 23:56:08Z knoop
	80	! Code annotations made doxygen readable
	81	!
[1607]	82	! 1606 2015-06-29 10:43:37Z maronga
	83	! Added preprocessor directive __netcdf to allow for compiling without netCDF.
	84	! Note, however, that RRTMG cannot be used without netCDF.
	85	!
[1591]	86	! 1590 2015-05-08 13:56:27Z maronga
	87	! Bugfix: definition of character strings requires same length for all elements
	88	!
[1588]	89	! 1587 2015-05-04 14:19:01Z maronga
	90	! Added albedo class for snow
	91	!
[1586]	92	! 1585 2015-04-30 07:05:52Z maronga
	93	! Added support for RRTMG
	94	!
[1572]	95	! 1571 2015-03-12 16:12:49Z maronga
	96	! Added missing KIND attribute. Removed upper-case variable names
	97	!
[1552]	98	! 1551 2015-03-03 14:18:16Z maronga
	99	! Added support for data output. Various variables have been renamed. Added
	100	! interface for different radiation schemes (currently: clear-sky, constant, and
	101	! RRTM (not yet implemented).
	102	!
[1497]	103	! 1496 2014-12-02 17:25:50Z maronga
	104	! Initial revision
	105	!
[1496]	106	!
	107	! Description:
	108	! ------------
[1682]	109	!> Radiation models and interfaces
[1826]	110	!> @todo move variable definitions used in radiation_init only to the subroutine
[1682]	111	!> as they are no longer required after initialization.
	112	!> @todo Output of full column vertical profiles used in RRTMG
	113	!> @todo Output of other rrtm arrays (such as volume mixing ratios)
	114	!> @todo Adapt for use with topography
	115	!>
	116	!> @note Many variables have a leading dummy dimension (0:0) in order to
	117	!> match the assume-size shape expected by the RRTMG model.
[1496]	118	!------------------------------------------------------------------------------!
[1682]	119	MODULE radiation_model_mod
	120
[1496]	121	USE arrays_3d, &
[2007]	122	ONLY: dzw, hyp, pt, q, ql, zu, zw
[1496]	123
[1585]	124	USE cloud_parameters, &
[1849]	125	ONLY: cp, l_d_cp, rho_l
[1585]	126
	127	USE constants, &
	128	ONLY: pi
	129
[1496]	130	USE control_parameters, &
[1585]	131	ONLY: cloud_droplets, cloud_physics, g, initializing_actions, &
[1691]	132	large_scale_forcing, lsf_surf, phi, pt_surface, rho_surface, &
[1585]	133	surface_pressure, time_since_reference_point
[1496]	134
	135	USE indices, &
[1585]	136	ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb_s_inner, nzb, nzt
[1496]	137
	138	USE kinds
	139
[1849]	140	USE microphysics_mod, &
	141	ONLY: nc_const, sigma_gc
	142
[1606]	143	#if defined ( __netcdf )
[1783]	144	USE NETCDF
[1606]	145	#endif
[1585]	146
	147	#if defined ( __rrtmg )
	148	USE parrrsw, &
	149	ONLY: naerec, nbndsw
[1551]	150
[1585]	151	USE parrrtm, &
	152	ONLY: nbndlw
	153
	154	USE rrtmg_lw_init, &
	155	ONLY: rrtmg_lw_ini
	156
	157	USE rrtmg_sw_init, &
	158	ONLY: rrtmg_sw_ini
	159
	160	USE rrtmg_lw_rad, &
	161	ONLY: rrtmg_lw
	162
	163	USE rrtmg_sw_rad, &
	164	ONLY: rrtmg_sw
	165	#endif
	166
	167
	168
[1496]	169	IMPLICIT NONE
	170
[1585]	171	CHARACTER(10) :: radiation_scheme = 'clear-sky' ! 'constant', 'clear-sky', or 'rrtmg'
[1551]	172
[1585]	173	!
	174	!-- Predefined Land surface classes (albedo_type) after Briegleb (1992)
[1788]	175	CHARACTER(37), DIMENSION(0:17), PARAMETER :: albedo_type_name = (/ &
[1590]	176	'user defined ', & ! 0
	177	'ocean ', & ! 1
	178	'mixed farming, tall grassland ', & ! 2
	179	'tall/medium grassland ', & ! 3
	180	'evergreen shrubland ', & ! 4
	181	'short grassland/meadow/shrubland ', & ! 5
	182	'evergreen needleleaf forest ', & ! 6
	183	'mixed deciduous evergreen forest ', & ! 7
	184	'deciduous forest ', & ! 8
	185	'tropical evergreen broadleaved forest', & ! 9
	186	'medium/tall grassland/woodland ', & ! 10
	187	'desert, sandy ', & ! 11
	188	'desert, rocky ', & ! 12
	189	'tundra ', & ! 13
	190	'land ice ', & ! 14
	191	'sea ice ', & ! 15
[1788]	192	'snow ', & ! 16
	193	'pavement/roads ' & ! 17
[1585]	194	/)
[1496]	195
[1682]	196	INTEGER(iwp) :: albedo_type = 5, & !< Albedo surface type (default: short grassland)
	197	day, & !< current day of the year
	198	day_init = 172, & !< day of the year at model start (21/06)
	199	dots_rad = 0 !< starting index for timeseries output
[1496]	200
[1757]	201	LOGICAL :: unscheduled_radiation_calls = .TRUE., & !< flag parameter indicating whether additional calls of the radiation code are allowed
	202	constant_albedo = .FALSE., & !< flag parameter indicating whether the albedo may change depending on zenith
	203	force_radiation_call = .FALSE., & !< flag parameter for unscheduled radiation calls
	204	lw_radiation = .TRUE., & !< flag parameter indicating whether longwave radiation shall be calculated
	205	radiation = .FALSE., & !< flag parameter indicating whether the radiation model is used
	206	sun_up = .TRUE., & !< flag parameter indicating whether the sun is up or down
[2007]	207	sw_radiation = .TRUE., & !< flag parameter indicing whether shortwave radiation shall be calculated
	208	sun_direction = .FALSE. !< flag parameter indicing whether solar direction shall be calculated
[1585]	209
[1496]	210
[1691]	211	REAL(wp), PARAMETER :: d_seconds_hour = 0.000277777777778_wp, & !< inverse of seconds per hour (1/3600)
	212	d_hours_day = 0.0416666666667_wp, & !< inverse of hours per day (1/24)
	213	sigma_sb = 5.67037321E-8_wp, & !< Stefan-Boltzmann constant
	214	solar_constant = 1368.0_wp !< solar constant at top of atmosphere
[1585]	215
[1691]	216	REAL(wp) :: albedo = 9999999.9_wp, & !< NAMELIST alpha
	217	albedo_lw_dif = 9999999.9_wp, & !< NAMELIST aldif
	218	albedo_lw_dir = 9999999.9_wp, & !< NAMELIST aldir
	219	albedo_sw_dif = 9999999.9_wp, & !< NAMELIST asdif
	220	albedo_sw_dir = 9999999.9_wp, & !< NAMELIST asdir
	221	decl_1, & !< declination coef. 1
	222	decl_2, & !< declination coef. 2
	223	decl_3, & !< declination coef. 3
	224	dt_radiation = 0.0_wp, & !< radiation model timestep
	225	emissivity = 0.98_wp, & !< NAMELIST surface emissivity
	226	lambda = 0.0_wp, & !< longitude in degrees
	227	lon = 0.0_wp, & !< longitude in radians
	228	lat = 0.0_wp, & !< latitude in radians
	229	net_radiation = 0.0_wp, & !< net radiation at surface
	230	skip_time_do_radiation = 0.0_wp, & !< Radiation model is not called before this time
	231	sky_trans, & !< sky transmissivity
	232	time_radiation = 0.0_wp, & !< time since last call of radiation code
	233	time_utc, & !< current time in UTC
	234	time_utc_init = 43200.0_wp !< UTC time at model start (noon)
	235
[2007]	236	REAL(wp), DIMENSION(0:0) :: zenith, & !< cosine of solar zenith angle
	237	sun_dir_lat, & !< solar directional vector in latitudes
	238	sun_dir_lon !< solar directional vector in longitudes
[1585]	239
[1496]	240	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
[1682]	241	alpha, & !< surface broadband albedo (used for clear-sky scheme)
[1709]	242	rad_lw_out_change_0, & !< change in LW out due to change in surface temperature
[1682]	243	rad_net, & !< net radiation at the surface
	244	rad_net_av !< average of rad_net
[1496]	245
[1585]	246	!
	247	!-- Land surface albedos for solar zenith angle of 60Â° after Briegleb (1992)
	248	!-- (shortwave, longwave, broadband): sw, lw, bb,
[1788]	249	REAL(wp), DIMENSION(0:2,1:17), PARAMETER :: albedo_pars = RESHAPE( (/&
[1585]	250	0.06_wp, 0.06_wp, 0.06_wp, & ! 1
	251	0.09_wp, 0.28_wp, 0.19_wp, & ! 2
	252	0.11_wp, 0.33_wp, 0.23_wp, & ! 3
	253	0.11_wp, 0.33_wp, 0.23_wp, & ! 4
	254	0.14_wp, 0.34_wp, 0.25_wp, & ! 5
	255	0.06_wp, 0.22_wp, 0.14_wp, & ! 6
	256	0.06_wp, 0.27_wp, 0.17_wp, & ! 7
	257	0.06_wp, 0.31_wp, 0.19_wp, & ! 8
	258	0.06_wp, 0.22_wp, 0.14_wp, & ! 9
	259	0.06_wp, 0.28_wp, 0.18_wp, & ! 10
	260	0.35_wp, 0.51_wp, 0.43_wp, & ! 11
	261	0.24_wp, 0.40_wp, 0.32_wp, & ! 12
	262	0.10_wp, 0.27_wp, 0.19_wp, & ! 13
	263	0.90_wp, 0.65_wp, 0.77_wp, & ! 14
[1587]	264	0.90_wp, 0.65_wp, 0.77_wp, & ! 15
[1788]	265	0.95_wp, 0.70_wp, 0.82_wp, & ! 16
	266	0.08_wp, 0.08_wp, 0.08_wp & ! 17
	267	/), (/ 3, 17 /) )
[1496]	268
[1585]	269	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: &
[1691]	270	rad_lw_cs_hr, & !< longwave clear sky radiation heating rate (K/s)
	271	rad_lw_cs_hr_av, & !< average of rad_lw_cs_hr
	272	rad_lw_hr, & !< longwave radiation heating rate (K/s)
	273	rad_lw_hr_av, & !< average of rad_sw_hr
	274	rad_lw_in, & !< incoming longwave radiation (W/m2)
	275	rad_lw_in_av, & !< average of rad_lw_in
	276	rad_lw_out, & !< outgoing longwave radiation (W/m2)
	277	rad_lw_out_av, & !< average of rad_lw_out
	278	rad_sw_cs_hr, & !< shortwave clear sky radiation heating rate (K/s)
	279	rad_sw_cs_hr_av, & !< average of rad_sw_cs_hr
	280	rad_sw_hr, & !< shortwave radiation heating rate (K/s)
	281	rad_sw_hr_av, & !< average of rad_sw_hr
[1682]	282	rad_sw_in, & !< incoming shortwave radiation (W/m2)
	283	rad_sw_in_av, & !< average of rad_sw_in
	284	rad_sw_out, & !< outgoing shortwave radiation (W/m2)
[1691]	285	rad_sw_out_av !< average of rad_sw_out
[1585]	286
[1691]	287
[1585]	288	!
	289	!-- Variables and parameters used in RRTMG only
	290	#if defined ( __rrtmg )
[1682]	291	CHARACTER(LEN=12) :: rrtm_input_file = "RAD_SND_DATA" !< name of the NetCDF input file (sounding data)
[1585]	292
	293
	294	!
	295	!-- Flag parameters for RRTMGS (should not be changed)
[1976]	296	INTEGER(iwp), PARAMETER :: rrtm_idrv = 1, & !< flag for longwave upward flux calculation option (0,1)
	297	rrtm_inflglw = 2, & !< flag for lw cloud optical properties (0,1,2)
[1682]	298	rrtm_iceflglw = 0, & !< flag for lw ice particle specifications (0,1,2,3)
	299	rrtm_liqflglw = 1, & !< flag for lw liquid droplet specifications
	300	rrtm_inflgsw = 2, & !< flag for sw cloud optical properties (0,1,2)
	301	rrtm_iceflgsw = 0, & !< flag for sw ice particle specifications (0,1,2,3)
	302	rrtm_liqflgsw = 1 !< flag for sw liquid droplet specifications
[1585]	303
	304	!
	305	!-- The following variables should be only changed with care, as this will
	306	!-- require further setting of some variables, which is currently not
	307	!-- implemented (aerosols, ice phase).
[1682]	308	INTEGER(iwp) :: nzt_rad, & !< upper vertical limit for radiation calculations
	309	rrtm_icld = 0, & !< cloud flag (0: clear sky column, 1: cloudy column)
[1976]	310	rrtm_iaer = 0 !< aerosol option flag (0: no aerosol layers, for lw only: 6 (requires setting of rrtm_sw_ecaer), 10: one or more aerosol layers (not implemented)
[1585]	311
[1788]	312	INTEGER(iwp) :: nc_stat !< local variable for storin the result of netCDF calls for error message handling
	313
[1682]	314	LOGICAL :: snd_exists = .FALSE. !< flag parameter to check whether a user-defined input files exists
[1585]	315
[1691]	316	REAL(wp), PARAMETER :: mol_mass_air_d_wv = 1.607793_wp !< molecular weight dry air / water vapor
[1585]	317
[1682]	318	REAL(wp), DIMENSION(:), ALLOCATABLE :: hyp_snd, & !< hypostatic pressure from sounding data (hPa)
	319	q_snd, & !< specific humidity from sounding data (kg/kg) - dummy at the moment
	320	rrtm_tsfc, & !< dummy array for storing surface temperature
	321	t_snd !< actual temperature from sounding data (hPa)
[1585]	322
[1691]	323	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: aldif, & !< longwave diffuse albedo solar angle of 60Â°
	324	aldir, & !< longwave direct albedo solar angle of 60Â°
	325	asdif, & !< shortwave diffuse albedo solar angle of 60Â°
	326	asdir, & !< shortwave direct albedo solar angle of 60Â°
	327	rrtm_ccl4vmr, & !< CCL4 volume mixing ratio (g/mol)
	328	rrtm_cfc11vmr, & !< CFC11 volume mixing ratio (g/mol)
	329	rrtm_cfc12vmr, & !< CFC12 volume mixing ratio (g/mol)
	330	rrtm_cfc22vmr, & !< CFC22 volume mixing ratio (g/mol)
	331	rrtm_ch4vmr, & !< CH4 volume mixing ratio
	332	rrtm_cicewp, & !< in-cloud ice water path (g/mÂ²)
	333	rrtm_cldfr, & !< cloud fraction (0,1)
	334	rrtm_cliqwp, & !< in-cloud liquid water path (g/mÂ²)
	335	rrtm_co2vmr, & !< CO2 volume mixing ratio (g/mol)
	336	rrtm_emis, & !< surface emissivity (0-1)
	337	rrtm_h2ovmr, & !< H2O volume mixing ratio
	338	rrtm_n2ovmr, & !< N2O volume mixing ratio
	339	rrtm_o2vmr, & !< O2 volume mixing ratio
	340	rrtm_o3vmr, & !< O3 volume mixing ratio
	341	rrtm_play, & !< pressure layers (hPa, zu-grid)
	342	rrtm_plev, & !< pressure layers (hPa, zw-grid)
	343	rrtm_reice, & !< cloud ice effective radius (microns)
	344	rrtm_reliq, & !< cloud water drop effective radius (microns)
	345	rrtm_tlay, & !< actual temperature (K, zu-grid)
	346	rrtm_tlev, & !< actual temperature (K, zw-grid)
	347	rrtm_lwdflx, & !< RRTM output of incoming longwave radiation flux (W/m2)
	348	rrtm_lwdflxc, & !< RRTM output of outgoing clear sky longwave radiation flux (W/m2)
	349	rrtm_lwuflx, & !< RRTM output of outgoing longwave radiation flux (W/m2)
	350	rrtm_lwuflxc, & !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
	351	rrtm_lwuflx_dt, & !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
	352	rrtm_lwuflxc_dt,& !< RRTM output of outgoing clear sky longwave radiation flux (W/m2)
	353	rrtm_lwhr, & !< RRTM output of longwave radiation heating rate (K/d)
	354	rrtm_lwhrc, & !< RRTM output of incoming longwave clear sky radiation heating rate (K/d)
	355	rrtm_swdflx, & !< RRTM output of incoming shortwave radiation flux (W/m2)
	356	rrtm_swdflxc, & !< RRTM output of outgoing clear sky shortwave radiation flux (W/m2)
	357	rrtm_swuflx, & !< RRTM output of outgoing shortwave radiation flux (W/m2)
	358	rrtm_swuflxc, & !< RRTM output of incoming clear sky shortwave radiation flux (W/m2)
	359	rrtm_swhr, & !< RRTM output of shortwave radiation heating rate (K/d)
	360	rrtm_swhrc !< RRTM output of incoming shortwave clear sky radiation heating rate (K/d)
[1585]	361
	362	!
	363	!-- Definition of arrays that are currently not used for calling RRTMG (due to setting of flag parameters)
[1682]	364	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: rad_lw_cs_in, & !< incoming clear sky longwave radiation (W/m2) (not used)
	365	rad_lw_cs_out, & !< outgoing clear sky longwave radiation (W/m2) (not used)
	366	rad_sw_cs_in, & !< incoming clear sky shortwave radiation (W/m2) (not used)
	367	rad_sw_cs_out, & !< outgoing clear sky shortwave radiation (W/m2) (not used)
	368	rrtm_aldif, & !< surface albedo for longwave diffuse radiation
	369	rrtm_aldir, & !< surface albedo for longwave direct radiation
	370	rrtm_asdif, & !< surface albedo for shortwave diffuse radiation
	371	rrtm_asdir, & !< surface albedo for shortwave direct radiation
	372	rrtm_lw_tauaer, & !< lw aerosol optical depth
	373	rrtm_lw_taucld, & !< lw in-cloud optical depth
	374	rrtm_sw_taucld, & !< sw in-cloud optical depth
	375	rrtm_sw_ssacld, & !< sw in-cloud single scattering albedo
	376	rrtm_sw_asmcld, & !< sw in-cloud asymmetry parameter
	377	rrtm_sw_fsfcld, & !< sw in-cloud forward scattering fraction
	378	rrtm_sw_tauaer, & !< sw aerosol optical depth
	379	rrtm_sw_ssaaer, & !< sw aerosol single scattering albedo
	380	rrtm_sw_asmaer, & !< sw aerosol asymmetry parameter
	381	rrtm_sw_ecaer !< sw aerosol optical detph at 0.55 microns (rrtm_iaer = 6 only)
[1691]	382
[1585]	383	#endif
	384
[1826]	385	INTERFACE radiation_check_data_output
	386	MODULE PROCEDURE radiation_check_data_output
	387	END INTERFACE radiation_check_data_output
[1496]	388
[1826]	389	INTERFACE radiation_check_data_output_pr
	390	MODULE PROCEDURE radiation_check_data_output_pr
	391	END INTERFACE radiation_check_data_output_pr
	392
	393	INTERFACE radiation_check_parameters
	394	MODULE PROCEDURE radiation_check_parameters
	395	END INTERFACE radiation_check_parameters
	396
[1551]	397	INTERFACE radiation_clearsky
	398	MODULE PROCEDURE radiation_clearsky
	399	END INTERFACE radiation_clearsky
[1853]	400
	401	INTERFACE radiation_constant
	402	MODULE PROCEDURE radiation_constant
	403	END INTERFACE radiation_constant
	404
[1976]	405	INTERFACE radiation_control
	406	MODULE PROCEDURE radiation_control
	407	END INTERFACE radiation_control
	408
	409	INTERFACE radiation_3d_data_averaging
	410	MODULE PROCEDURE radiation_3d_data_averaging
	411	END INTERFACE radiation_3d_data_averaging
	412
	413	INTERFACE radiation_data_output_2d
	414	MODULE PROCEDURE radiation_data_output_2d
	415	END INTERFACE radiation_data_output_2d
	416
	417	INTERFACE radiation_data_output_3d
	418	MODULE PROCEDURE radiation_data_output_3d
	419	END INTERFACE radiation_data_output_3d
	420
	421	INTERFACE radiation_data_output_mask
	422	MODULE PROCEDURE radiation_data_output_mask
	423	END INTERFACE radiation_data_output_mask
	424
	425	INTERFACE radiation_define_netcdf_grid
	426	MODULE PROCEDURE radiation_define_netcdf_grid
	427	END INTERFACE radiation_define_netcdf_grid
	428
[1826]	429	INTERFACE radiation_header
	430	MODULE PROCEDURE radiation_header
	431	END INTERFACE radiation_header
	432
	433	INTERFACE radiation_init
	434	MODULE PROCEDURE radiation_init
	435	END INTERFACE radiation_init
[1496]	436
[1826]	437	INTERFACE radiation_parin
	438	MODULE PROCEDURE radiation_parin
	439	END INTERFACE radiation_parin
	440
[1585]	441	INTERFACE radiation_rrtmg
	442	MODULE PROCEDURE radiation_rrtmg
	443	END INTERFACE radiation_rrtmg
[1551]	444
[1585]	445	INTERFACE radiation_tendency
	446	MODULE PROCEDURE radiation_tendency
	447	MODULE PROCEDURE radiation_tendency_ij
	448	END INTERFACE radiation_tendency
[1551]	449
[1976]	450	INTERFACE radiation_read_restart_data
	451	MODULE PROCEDURE radiation_read_restart_data
	452	END INTERFACE radiation_read_restart_data
	453
	454	INTERFACE radiation_last_actions
	455	MODULE PROCEDURE radiation_last_actions
	456	END INTERFACE radiation_last_actions
	457
[1496]	458	SAVE
	459
	460	PRIVATE
	461
[1826]	462	!
[1976]	463	!-- Public functions / NEEDS SORTING
[1826]	464	PUBLIC radiation_check_data_output, radiation_check_data_output_pr, &
[1976]	465	radiation_check_parameters, radiation_control, &
	466	radiation_header, radiation_init, radiation_parin, &
	467	radiation_3d_data_averaging, radiation_tendency, &
	468	radiation_data_output_2d, radiation_data_output_3d, &
	469	radiation_define_netcdf_grid, radiation_last_actions, &
	470	radiation_read_restart_data, radiation_data_output_mask
[1826]	471
	472	!
[1976]	473	!-- Public variables and constants / NEEDS SORTING
[1826]	474	PUBLIC dots_rad, dt_radiation, force_radiation_call, &
	475	rad_net, rad_net_av, radiation, radiation_scheme, rad_lw_in, &
	476	rad_lw_in_av, rad_lw_out, rad_lw_out_av, rad_lw_out_change_0, &
	477	rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av, rad_sw_in, &
	478	rad_sw_in_av, rad_sw_out, rad_sw_out_av, rad_sw_cs_hr, &
[2007]	479	rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av, sigma_sb, &
	480	skip_time_do_radiation, time_radiation, unscheduled_radiation_calls,&
	481	zenith, calc_zenith, sun_direction, sun_dir_lat, sun_dir_lon, &
	482	day_init, time_utc_init
[1496]	483
[1691]	484
[1585]	485	#if defined ( __rrtmg )
[1976]	486	PUBLIC rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir
[1585]	487	#endif
[1496]	488
	489	CONTAINS
	490
[1976]	491
[1496]	492	!------------------------------------------------------------------------------!
	493	! Description:
	494	! ------------
[1976]	495	!> This subroutine controls the calls of the radiation schemes
	496	!------------------------------------------------------------------------------!
	497	SUBROUTINE radiation_control
	498
	499
	500	IMPLICIT NONE
	501
	502
	503	SELECT CASE ( TRIM( radiation_scheme ) )
	504
	505	CASE ( 'constant' )
	506	CALL radiation_constant
	507
	508	CASE ( 'clear-sky' )
	509	CALL radiation_clearsky
	510
	511	CASE ( 'rrtmg' )
	512	CALL radiation_rrtmg
	513
	514	CASE DEFAULT
	515
	516	END SELECT
	517
	518
	519	END SUBROUTINE radiation_control
	520
	521	!------------------------------------------------------------------------------!
	522	! Description:
	523	! ------------
[1826]	524	!> Check data output for radiation model
	525	!------------------------------------------------------------------------------!
	526	SUBROUTINE radiation_check_data_output( var, unit, i, ilen, k )
	527
	528
	529	USE control_parameters, &
	530	ONLY: data_output, message_string
	531
	532	IMPLICIT NONE
	533
	534	CHARACTER (LEN=*) :: unit !<
	535	CHARACTER (LEN=*) :: var !<
	536
	537	INTEGER(iwp) :: i
	538	INTEGER(iwp) :: ilen
	539	INTEGER(iwp) :: k
	540
	541	SELECT CASE ( TRIM( var ) )
	542
[1976]	543	CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_sw_cs_hr', 'rad_sw_hr' )
[1826]	544	IF ( .NOT. radiation .OR. radiation_scheme /= 'rrtmg' ) THEN
	545	message_string = '"output of "' // TRIM( var ) // '" requi' // &
	546	'res radiation = .TRUE. and ' // &
	547	'radiation_scheme = "rrtmg"'
	548	CALL message( 'check_parameters', 'PA0406', 1, 2, 0, 6, 0 )
	549	ENDIF
	550	unit = 'W/m2'
	551
	552	CASE ( 'rad_net', 'rrtm_aldif', 'rrtm_aldir', 'rrtm_asdif', &
	553	'rrtm_asdir*' )
	554	IF ( k == 0 .OR. data_output(i)(ilen-2:ilen) /= '_xy' ) THEN
	555	message_string = 'illegal value for data_output: "' // &
	556	TRIM( var ) // '" & only 2d-horizontal ' // &
	557	'cross sections are allowed for this value'
	558	CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
	559	ENDIF
	560	IF ( .NOT. radiation .OR. radiation_scheme /= "rrtmg" ) THEN
	561	IF ( TRIM( var ) == 'rrtm_aldif*' .OR. &
	562	TRIM( var ) == 'rrtm_aldir*' .OR. &
	563	TRIM( var ) == 'rrtm_asdif*' .OR. &
	564	TRIM( var ) == 'rrtm_asdir*' ) &
	565	THEN
	566	message_string = 'output of "' // TRIM( var ) // '" require'&
	567	// 's radiation = .TRUE. and radiation_sch'&
	568	// 'eme = "rrtmg"'
	569	CALL message( 'check_parameters', 'PA0409', 1, 2, 0, 6, 0 )
	570	ENDIF
	571	ENDIF
	572
	573	IF ( TRIM( var ) == 'rad_net*' ) unit = 'W/m2'
	574	IF ( TRIM( var ) == 'rrtm_aldif*' ) unit = ''
	575	IF ( TRIM( var ) == 'rrtm_aldir*' ) unit = ''
	576	IF ( TRIM( var ) == 'rrtm_asdif*' ) unit = ''
	577	IF ( TRIM( var ) == 'rrtm_asdir*' ) unit = ''
	578
	579	CASE DEFAULT
	580	unit = 'illegal'
	581
	582	END SELECT
	583
	584
	585	END SUBROUTINE radiation_check_data_output
	586
	587	!------------------------------------------------------------------------------!
	588	! Description:
	589	! ------------
	590	!> Check data output of profiles for radiation model
	591	!------------------------------------------------------------------------------!
	592	SUBROUTINE radiation_check_data_output_pr( variable, var_count, unit, dopr_unit )
	593
	594	USE arrays_3d, &
	595	ONLY: zu
	596
	597	USE control_parameters, &
	598	ONLY: data_output_pr, message_string
	599
	600	USE indices
	601
	602	USE profil_parameter
	603
	604	USE statistics
	605
	606	IMPLICIT NONE
	607
	608	CHARACTER (LEN=*) :: unit !<
	609	CHARACTER (LEN=*) :: variable !<
	610	CHARACTER (LEN=*) :: dopr_unit !< local value of dopr_unit
	611
	612	INTEGER(iwp) :: user_pr_index !<
	613	INTEGER(iwp) :: var_count !<
	614
	615	SELECT CASE ( TRIM( variable ) )
	616
	617	CASE ( 'rad_net' )
	618	IF ( ( .NOT. radiation ) .OR. radiation_scheme == 'constant' )&
	619	THEN
	620	message_string = 'data_output_pr = ' // &
	621	TRIM( data_output_pr(var_count) ) // ' is' // &
	622	'not available for radiation = .FALSE. or ' //&
	623	'radiation_scheme = "constant"'
	624	CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
	625	ELSE
	626	dopr_index(var_count) = 101
	627	dopr_unit = 'W/m2'
	628	hom(:,2,101,:) = SPREAD( zw, 2, statistic_regions+1 )
	629	unit = dopr_unit
	630	ENDIF
	631
	632	CASE ( 'rad_lw_in' )
	633	IF ( ( .NOT. radiation) .OR. radiation_scheme == 'constant' ) &
	634	THEN
	635	message_string = 'data_output_pr = ' // &
	636	TRIM( data_output_pr(var_count) ) // ' is' // &
	637	'not available for radiation = .FALSE. or ' //&
	638	'radiation_scheme = "constant"'
	639	CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
	640	ELSE
	641	dopr_index(var_count) = 102
	642	dopr_unit = 'W/m2'
	643	hom(:,2,102,:) = SPREAD( zw, 2, statistic_regions+1 )
	644	unit = dopr_unit
	645	ENDIF
	646
	647	CASE ( 'rad_lw_out' )
	648	IF ( ( .NOT. radiation ) .OR. radiation_scheme == 'constant' ) &
	649	THEN
	650	message_string = 'data_output_pr = ' // &
	651	TRIM( data_output_pr(var_count) ) // ' is' // &
	652	'not available for radiation = .FALSE. or ' //&
	653	'radiation_scheme = "constant"'
	654	CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
	655	ELSE
	656	dopr_index(var_count) = 103
	657	dopr_unit = 'W/m2'
	658	hom(:,2,103,:) = SPREAD( zw, 2, statistic_regions+1 )
	659	unit = dopr_unit
	660	ENDIF
	661
	662	CASE ( 'rad_sw_in' )
	663	IF ( ( .NOT. radiation ) .OR. radiation_scheme == 'constant' ) &
	664	THEN
	665	message_string = 'data_output_pr = ' // &
	666	TRIM( data_output_pr(var_count) ) // ' is' // &
	667	'not available for radiation = .FALSE. or ' //&
	668	'radiation_scheme = "constant"'
	669	CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
	670	ELSE
	671	dopr_index(var_count) = 104
	672	dopr_unit = 'W/m2'
	673	hom(:,2,104,:) = SPREAD( zw, 2, statistic_regions+1 )
	674	unit = dopr_unit
	675	ENDIF
	676
	677	CASE ( 'rad_sw_out')
	678	IF ( ( .NOT. radiation ) .OR. radiation_scheme == 'constant' )&
	679	THEN
	680	message_string = 'data_output_pr = ' // &
	681	TRIM( data_output_pr(var_count) ) // ' is' // &
	682	'not available for radiation = .FALSE. or ' //&
	683	'radiation_scheme = "constant"'
	684	CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
	685	ELSE
	686	dopr_index(var_count) = 105
	687	dopr_unit = 'W/m2'
	688	hom(:,2,105,:) = SPREAD( zw, 2, statistic_regions+1 )
	689	unit = dopr_unit
	690	ENDIF
	691
	692	CASE ( 'rad_lw_cs_hr' )
	693	IF ( ( .NOT. radiation ) .OR. radiation_scheme /= 'rrtmg' ) &
	694	THEN
	695	message_string = 'data_output_pr = ' // &
	696	TRIM( data_output_pr(var_count) ) // ' is' // &
	697	'not available for radiation = .FALSE. or ' //&
	698	'radiation_scheme /= "rrtmg"'
	699	CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
	700	ELSE
	701	dopr_index(var_count) = 106
	702	dopr_unit = 'K/h'
	703	hom(:,2,106,:) = SPREAD( zu, 2, statistic_regions+1 )
	704	unit = dopr_unit
	705	ENDIF
	706
	707	CASE ( 'rad_lw_hr' )
	708	IF ( ( .NOT. radiation ) .OR. radiation_scheme /= 'rrtmg' ) &
	709	THEN
	710	message_string = 'data_output_pr = ' // &
	711	TRIM( data_output_pr(var_count) ) // ' is' // &
	712	'not available for radiation = .FALSE. or ' //&
	713	'radiation_scheme /= "rrtmg"'
	714	CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
	715	ELSE
	716	dopr_index(var_count) = 107
	717	dopr_unit = 'K/h'
	718	hom(:,2,107,:) = SPREAD( zu, 2, statistic_regions+1 )
	719	unit = dopr_unit
	720	ENDIF
	721
	722	CASE ( 'rad_sw_cs_hr' )
	723	IF ( ( .NOT. radiation ) .OR. radiation_scheme /= 'rrtmg' ) &
	724	THEN
	725	message_string = 'data_output_pr = ' // &
	726	TRIM( data_output_pr(var_count) ) // ' is' // &
	727	'not available for radiation = .FALSE. or ' //&
	728	'radiation_scheme /= "rrtmg"'
	729	CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
	730	ELSE
	731	dopr_index(var_count) = 108
	732	dopr_unit = 'K/h'
	733	hom(:,2,108,:) = SPREAD( zu, 2, statistic_regions+1 )
	734	unit = dopr_unit
	735	ENDIF
	736
	737	CASE ( 'rad_sw_hr' )
	738	IF ( ( .NOT. radiation ) .OR. radiation_scheme /= 'rrtmg' ) &
	739	THEN
	740	message_string = 'data_output_pr = ' // &
	741	TRIM( data_output_pr(var_count) ) // ' is' // &
	742	'not available for radiation = .FALSE. or ' //&
	743	'radiation_scheme /= "rrtmg"'
	744	CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
	745	ELSE
	746	dopr_index(var_count) = 109
	747	dopr_unit = 'K/h'
	748	hom(:,2,109,:) = SPREAD( zu, 2, statistic_regions+1 )
	749	unit = dopr_unit
	750	ENDIF
	751
	752
	753	CASE DEFAULT
	754	unit = 'illegal'
	755
	756	END SELECT
	757
	758
	759	END SUBROUTINE radiation_check_data_output_pr
	760
	761
	762	!------------------------------------------------------------------------------!
	763	! Description:
	764	! ------------
	765	!> Check parameters routine for radiation model
	766	!------------------------------------------------------------------------------!
	767	SUBROUTINE radiation_check_parameters
	768
	769	USE control_parameters, &
	770	ONLY: message_string, topography
	771
	772
	773	IMPLICIT NONE
[2007]	774
	775	LOGICAL :: urban_surface_af = .FALSE. !< auxiliary flag used for parameter check
[1826]	776
[2007]	777
[1826]	778	IF ( radiation_scheme /= 'constant' .AND. &
	779	radiation_scheme /= 'clear-sky' .AND. &
	780	radiation_scheme /= 'rrtmg' ) THEN
	781	message_string = 'unknown radiation_scheme = '// &
	782	TRIM( radiation_scheme )
	783	CALL message( 'check_parameters', 'PA0405', 1, 2, 0, 6, 0 )
	784	ELSEIF ( radiation_scheme == 'rrtmg' ) THEN
	785	#if ! defined ( __rrtmg )
	786	message_string = 'radiation_scheme = "rrtmg" requires ' // &
	787	'compilation of PALM with pre-processor ' // &
	788	'directive -D__rrtmg'
	789	CALL message( 'check_parameters', 'PA0407', 1, 2, 0, 6, 0 )
	790	#endif
	791	#if defined ( __rrtmg ) && ! defined( __netcdf )
	792	message_string = 'radiation_scheme = "rrtmg" requires ' // &
	793	'the use of NetCDF (preprocessor directive ' // &
	794	'-D__netcdf'
	795	CALL message( 'check_parameters', 'PA0412', 1, 2, 0, 6, 0 )
	796	#endif
	797
	798	ENDIF
	799
	800	IF ( albedo_type == 0 .AND. albedo == 9999999.9_wp .AND. &
	801	radiation_scheme == 'clear-sky') THEN
	802	message_string = 'radiation_scheme = "clear-sky" in combination' // &
	803	'with albedo_type = 0 requires setting of albedo'// &
	804	' /= 9999999.9'
	805	CALL message( 'check_parameters', 'PA0410', 1, 2, 0, 6, 0 )
	806	ENDIF
	807
	808	IF ( albedo_type == 0 .AND. radiation_scheme == 'rrtmg' .AND. &
	809	( albedo_lw_dif == 9999999.9_wp .OR. albedo_lw_dir == 9999999.9_wp&
	810	.OR. albedo_sw_dif == 9999999.9_wp .OR. albedo_sw_dir == 9999999.9_wp&
	811	) ) THEN
	812	message_string = 'radiation_scheme = "rrtmg" in combination' // &
	813	'with albedo_type = 0 requires setting of ' // &
	814	'albedo_lw_dif /= 9999999.9' // &
	815	'albedo_lw_dir /= 9999999.9' // &
	816	'albedo_sw_dif /= 9999999.9 and' // &
	817	'albedo_sw_dir /= 9999999.9'
	818	CALL message( 'check_parameters', 'PA0411', 1, 2, 0, 6, 0 )
	819	ENDIF
	820
[2007]	821	!
	822	!-- The following paramter check is temporarily extended by the urban_surface
	823	!-- flag, until a better solution comes up to omit this check in case of
	824	!-- urban surface model is used.
	825	!-- Routine get_usm_info provides the value for the urban_surface flag,
	826	!-- because the value cannot be retrieved via USE due to circular dependencies
	827	!-- between modules radiation_model_mod and urban_surface_mod.
	828	CALL get_usm_info( urban_surface_af )
	829	IF ( topography /= 'flat' .AND. .NOT. urban_surface_af ) THEN
[1826]	830	message_string = 'radiation scheme cannot be used ' // &
	831	'in combination with topography /= "flat"'
	832	CALL message( 'check_parameters', 'PA0414', 1, 2, 0, 6, 0 )
	833	ENDIF
	834
	835	END SUBROUTINE radiation_check_parameters
	836
	837
	838	!------------------------------------------------------------------------------!
	839	! Description:
	840	! ------------
[1682]	841	!> Initialization of the radiation model
[1496]	842	!------------------------------------------------------------------------------!
[1826]	843	SUBROUTINE radiation_init
[1496]	844
	845	IMPLICIT NONE
	846
[1585]	847	!
	848	!-- Allocate array for storing the surface net radiation
	849	IF ( .NOT. ALLOCATED ( rad_net ) ) THEN
	850	ALLOCATE ( rad_net(nysg:nyng,nxlg:nxrg) )
	851	rad_net = 0.0_wp
	852	ENDIF
[1496]	853
	854	!
[1709]	855	!-- Allocate array for storing the surface net radiation
	856	IF ( .NOT. ALLOCATED ( rad_lw_out_change_0 ) ) THEN
	857	ALLOCATE ( rad_lw_out_change_0(nysg:nyng,nxlg:nxrg) )
	858	rad_lw_out_change_0 = 0.0_wp
	859	ENDIF
	860
	861	!
[1551]	862	!-- Fix net radiation in case of radiation_scheme = 'constant'
[1585]	863	IF ( radiation_scheme == 'constant' ) THEN
[1551]	864	rad_net = net_radiation
[1853]	865	! radiation = .FALSE.
[1551]	866	!
[1585]	867	!-- Calculate orbital constants
	868	ELSE
[1551]	869	decl_1 = SIN(23.45_wp * pi / 180.0_wp)
	870	decl_2 = 2.0_wp * pi / 365.0_wp
	871	decl_3 = decl_2 * 81.0_wp
[1585]	872	lat = phi * pi / 180.0_wp
	873	lon = lambda * pi / 180.0_wp
	874	ENDIF
	875
[1976]	876	IF ( radiation_scheme == 'clear-sky' .OR. &
	877	radiation_scheme == 'constant') THEN
[1585]	878
	879	ALLOCATE ( alpha(nysg:nyng,nxlg:nxrg) )
	880
	881	IF ( .NOT. ALLOCATED ( rad_sw_in ) ) THEN
	882	ALLOCATE ( rad_sw_in(0:0,nysg:nyng,nxlg:nxrg) )
	883	ENDIF
	884	IF ( .NOT. ALLOCATED ( rad_sw_out ) ) THEN
	885	ALLOCATE ( rad_sw_out(0:0,nysg:nyng,nxlg:nxrg) )
	886	ENDIF
	887
	888	IF ( .NOT. ALLOCATED ( rad_sw_in_av ) ) THEN
	889	ALLOCATE ( rad_sw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
	890	ENDIF
	891	IF ( .NOT. ALLOCATED ( rad_sw_out_av ) ) THEN
	892	ALLOCATE ( rad_sw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
	893	ENDIF
	894
	895	IF ( .NOT. ALLOCATED ( rad_lw_in ) ) THEN
	896	ALLOCATE ( rad_lw_in(0:0,nysg:nyng,nxlg:nxrg) )
	897	ENDIF
[1856]	898	IF ( .NOT. ALLOCATED ( rad_lw_out ) ) THEN
	899	ALLOCATE ( rad_lw_out(0:0,nysg:nyng,nxlg:nxrg) )
	900	ENDIF
[1585]	901
	902	IF ( .NOT. ALLOCATED ( rad_lw_in_av ) ) THEN
	903	ALLOCATE ( rad_lw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
	904	ENDIF
	905	IF ( .NOT. ALLOCATED ( rad_lw_out_av ) ) THEN
	906	ALLOCATE ( rad_lw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
	907	ENDIF
	908
	909	rad_sw_in = 0.0_wp
	910	rad_sw_out = 0.0_wp
	911	rad_lw_in = 0.0_wp
	912	rad_lw_out = 0.0_wp
	913
[1496]	914	!
[1585]	915	!-- Overwrite albedo if manually set in parameter file
	916	IF ( albedo_type /= 0 .AND. albedo == 9999999.9_wp ) THEN
	917	albedo = albedo_pars(2,albedo_type)
	918	ENDIF
	919
	920	alpha = albedo
	921
	922	!
	923	!-- Initialization actions for RRTMG
	924	ELSEIF ( radiation_scheme == 'rrtmg' ) THEN
	925	#if defined ( __rrtmg )
	926	!
	927	!-- Allocate albedos
	928	ALLOCATE ( rrtm_aldif(0:0,nysg:nyng,nxlg:nxrg) )
	929	ALLOCATE ( rrtm_aldir(0:0,nysg:nyng,nxlg:nxrg) )
	930	ALLOCATE ( rrtm_asdif(0:0,nysg:nyng,nxlg:nxrg) )
	931	ALLOCATE ( rrtm_asdir(0:0,nysg:nyng,nxlg:nxrg) )
	932	ALLOCATE ( aldif(nysg:nyng,nxlg:nxrg) )
	933	ALLOCATE ( aldir(nysg:nyng,nxlg:nxrg) )
	934	ALLOCATE ( asdif(nysg:nyng,nxlg:nxrg) )
	935	ALLOCATE ( asdir(nysg:nyng,nxlg:nxrg) )
	936
	937	IF ( albedo_type /= 0 ) THEN
	938	IF ( albedo_lw_dif == 9999999.9_wp ) THEN
	939	albedo_lw_dif = albedo_pars(0,albedo_type)
	940	albedo_lw_dir = albedo_lw_dif
	941	ENDIF
	942	IF ( albedo_sw_dif == 9999999.9_wp ) THEN
	943	albedo_sw_dif = albedo_pars(1,albedo_type)
	944	albedo_sw_dir = albedo_sw_dif
	945	ENDIF
	946	ENDIF
	947
	948	aldif(:,:) = albedo_lw_dif
	949	aldir(:,:) = albedo_lw_dir
	950	asdif(:,:) = albedo_sw_dif
	951	asdir(:,:) = albedo_sw_dir
	952	!
	953	!-- Calculate initial values of current (cosine of) the zenith angle and
	954	!-- whether the sun is up
	955	CALL calc_zenith
	956	!
	957	!-- Calculate initial surface albedo
	958	IF ( .NOT. constant_albedo ) THEN
	959	CALL calc_albedo
	960	ELSE
	961	rrtm_aldif(0,:,:) = aldif(:,:)
	962	rrtm_aldir(0,:,:) = aldir(:,:)
	963	rrtm_asdif(0,:,:) = asdif(:,:)
	964	rrtm_asdir(0,:,:) = asdir(:,:)
	965	ENDIF
	966
	967	!
	968	!-- Allocate surface emissivity
	969	ALLOCATE ( rrtm_emis(0:0,1:nbndlw+1) )
	970	rrtm_emis = emissivity
	971
	972	!
	973	!-- Allocate 3d arrays of radiative fluxes and heating rates
	974	IF ( .NOT. ALLOCATED ( rad_sw_in ) ) THEN
	975	ALLOCATE ( rad_sw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	976	rad_sw_in = 0.0_wp
	977	ENDIF
	978
	979	IF ( .NOT. ALLOCATED ( rad_sw_in_av ) ) THEN
	980	ALLOCATE ( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	981	ENDIF
	982
	983	IF ( .NOT. ALLOCATED ( rad_sw_out ) ) THEN
	984	ALLOCATE ( rad_sw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1691]	985	rad_sw_out = 0.0_wp
[1585]	986	ENDIF
	987
	988	IF ( .NOT. ALLOCATED ( rad_sw_out_av ) ) THEN
	989	ALLOCATE ( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	990	ENDIF
	991
[1691]	992	IF ( .NOT. ALLOCATED ( rad_sw_hr ) ) THEN
	993	ALLOCATE ( rad_sw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	994	rad_sw_hr = 0.0_wp
	995	ENDIF
[1585]	996
[1691]	997	IF ( .NOT. ALLOCATED ( rad_sw_hr_av ) ) THEN
	998	ALLOCATE ( rad_sw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	999	rad_sw_hr_av = 0.0_wp
	1000	ENDIF
	1001
	1002	IF ( .NOT. ALLOCATED ( rad_sw_cs_hr ) ) THEN
	1003	ALLOCATE ( rad_sw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	1004	rad_sw_cs_hr = 0.0_wp
	1005	ENDIF
	1006
	1007	IF ( .NOT. ALLOCATED ( rad_sw_cs_hr_av ) ) THEN
	1008	ALLOCATE ( rad_sw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	1009	rad_sw_cs_hr_av = 0.0_wp
	1010	ENDIF
	1011
[1585]	1012	IF ( .NOT. ALLOCATED ( rad_lw_in ) ) THEN
	1013	ALLOCATE ( rad_lw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	1014	rad_lw_in = 0.0_wp
	1015	ENDIF
	1016
	1017	IF ( .NOT. ALLOCATED ( rad_lw_in_av ) ) THEN
	1018	ALLOCATE ( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	1019	ENDIF
	1020
	1021	IF ( .NOT. ALLOCATED ( rad_lw_out ) ) THEN
	1022	ALLOCATE ( rad_lw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	1023	rad_lw_out = 0.0_wp
	1024	ENDIF
	1025
	1026	IF ( .NOT. ALLOCATED ( rad_lw_out_av ) ) THEN
	1027	ALLOCATE ( rad_lw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	1028	ENDIF
	1029
[1691]	1030	IF ( .NOT. ALLOCATED ( rad_lw_hr ) ) THEN
	1031	ALLOCATE ( rad_lw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	1032	rad_lw_hr = 0.0_wp
	1033	ENDIF
	1034
	1035	IF ( .NOT. ALLOCATED ( rad_lw_hr_av ) ) THEN
	1036	ALLOCATE ( rad_lw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	1037	rad_lw_hr_av = 0.0_wp
	1038	ENDIF
	1039
	1040	IF ( .NOT. ALLOCATED ( rad_lw_cs_hr ) ) THEN
	1041	ALLOCATE ( rad_lw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	1042	rad_lw_cs_hr = 0.0_wp
	1043	ENDIF
	1044
	1045	IF ( .NOT. ALLOCATED ( rad_lw_cs_hr_av ) ) THEN
	1046	ALLOCATE ( rad_lw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	1047	rad_lw_cs_hr_av = 0.0_wp
	1048	ENDIF
	1049
	1050	ALLOCATE ( rad_sw_cs_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	1051	ALLOCATE ( rad_sw_cs_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1585]	1052	rad_sw_cs_in = 0.0_wp
	1053	rad_sw_cs_out = 0.0_wp
	1054
[1691]	1055	ALLOCATE ( rad_lw_cs_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	1056	ALLOCATE ( rad_lw_cs_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1585]	1057	rad_lw_cs_in = 0.0_wp
	1058	rad_lw_cs_out = 0.0_wp
	1059
	1060	!
	1061	!-- Allocate dummy array for storing surface temperature
	1062	ALLOCATE ( rrtm_tsfc(1) )
	1063
	1064	!
	1065	!-- Initialize RRTMG
	1066	IF ( lw_radiation ) CALL rrtmg_lw_ini ( cp )
	1067	IF ( sw_radiation ) CALL rrtmg_sw_ini ( cp )
	1068
	1069	!
	1070	!-- Set input files for RRTMG
	1071	INQUIRE(FILE="RAD_SND_DATA", EXIST=snd_exists)
	1072	IF ( .NOT. snd_exists ) THEN
	1073	rrtm_input_file = "rrtmg_lw.nc"
	1074	ENDIF
	1075
	1076	!
	1077	!-- Read vertical layers for RRTMG from sounding data
	1078	!-- The routine provides nzt_rad, hyp_snd(1:nzt_rad),
	1079	!-- t_snd(nzt+2:nzt_rad), rrtm_play(1:nzt_rad), rrtm_plev(1_nzt_rad+1),
	1080	!-- rrtm_tlay(nzt+2:nzt_rad), rrtm_tlev(nzt+2:nzt_rad+1)
	1081	CALL read_sounding_data
	1082
	1083	!
	1084	!-- Read trace gas profiles from file. This routine provides
	1085	!-- the rrtm_ arrays (1:nzt_rad+1)
	1086	CALL read_trace_gas_data
	1087	#endif
[1551]	1088	ENDIF
[1585]	1089
[1551]	1090	!
[1585]	1091	!-- Perform user actions if required
	1092	CALL user_init_radiation
	1093
	1094	!
	1095	!-- Calculate radiative fluxes at model start
	1096	IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
[1853]	1097
	1098	SELECT CASE ( radiation_scheme )
	1099	CASE ( 'rrtmg' )
	1100	CALL radiation_rrtmg
	1101	CASE ( 'clear-sky' )
	1102	CALL radiation_clearsky
	1103	CASE ( 'constant' )
	1104	CALL radiation_constant
	1105	CASE DEFAULT
	1106	END SELECT
	1107
[1585]	1108	ENDIF
	1109
[1496]	1110	RETURN
	1111
[1826]	1112	END SUBROUTINE radiation_init
[1496]	1113
	1114
	1115	!------------------------------------------------------------------------------!
	1116	! Description:
	1117	! ------------
[1682]	1118	!> A simple clear sky radiation model
[1496]	1119	!------------------------------------------------------------------------------!
[1551]	1120	SUBROUTINE radiation_clearsky
[1496]	1121
[1585]	1122
[1496]	1123	IMPLICIT NONE
	1124
[1691]	1125	INTEGER(iwp) :: i, j, k !< loop indices
	1126	REAL(wp) :: exn, & !< Exner functions at surface
[1709]	1127	exn1, & !< Exner functions at first grid level
	1128	pt1 !< potential temperature at first grid level
[1585]	1129
[1496]	1130	!
[1585]	1131	!-- Calculate current zenith angle
	1132	CALL calc_zenith
	1133
	1134	!
	1135	!-- Calculate sky transmissivity
	1136	sky_trans = 0.6_wp + 0.2_wp * zenith(0)
	1137
	1138	!
	1139	!-- Calculate value of the Exner function
	1140	exn = (surface_pressure / 1000.0_wp )**0.286_wp
	1141	!
	1142	!-- Calculate radiation fluxes and net radiation (rad_net) for each grid
	1143	!-- point
[1709]	1144	DO i = nxlg, nxrg
	1145	DO j = nysg, nyng
[1585]	1146	k = nzb_s_inner(j,i)
[1691]	1147
[1709]	1148	exn1 = (hyp(k+1) / 100000.0_wp )**0.286_wp
[1691]	1149
[1585]	1150	rad_sw_in(0,j,i) = solar_constant * sky_trans * zenith(0)
	1151	rad_sw_out(0,j,i) = alpha(j,i) * rad_sw_in(0,j,i)
[1691]	1152	rad_lw_out(0,j,i) = emissivity * sigma_sb * (pt(k,j,i) * exn)**4
[1585]	1153
[1691]	1154	IF ( cloud_physics ) THEN
[1709]	1155	pt1 = pt(k+1,j,i) + l_d_cp / exn1 * ql(k+1,j,i)
	1156	rad_lw_in(0,j,i) = 0.8_wp * sigma_sb * (pt1 * exn1)**4
[1691]	1157	ELSE
[1709]	1158	rad_lw_in(0,j,i) = 0.8_wp * sigma_sb * (pt(k+1,j,i) * exn1)**4
[1691]	1159	ENDIF
	1160
	1161	rad_net(j,i) = rad_sw_in(0,j,i) - rad_sw_out(0,j,i) &
	1162	+ rad_lw_in(0,j,i) - rad_lw_out(0,j,i)
	1163
[1976]	1164
	1165	rad_lw_out_change_0(j,i) = 3.0_wp * sigma_sb * emissivity &
	1166	* (pt(k,j,i) * exn) ** 3
	1167
[1585]	1168	ENDDO
	1169	ENDDO
	1170
	1171	END SUBROUTINE radiation_clearsky
	1172
	1173
	1174	!------------------------------------------------------------------------------!
	1175	! Description:
	1176	! ------------
[1853]	1177	!> This scheme keeps the prescribed net radiation constant during the run
	1178	!------------------------------------------------------------------------------!
	1179	SUBROUTINE radiation_constant
	1180
	1181
	1182	IMPLICIT NONE
	1183
	1184	INTEGER(iwp) :: i, j, k !< loop indices
	1185	REAL(wp) :: exn, & !< Exner functions at surface
[1976]	1186	exn1, & !< Exner functions at first grid level
[1853]	1187	pt1 !< potential temperature at first grid level
	1188
	1189	!
	1190	!-- Calculate value of the Exner function
	1191	exn = (surface_pressure / 1000.0_wp )**0.286_wp
	1192	!
[1976]	1193	!-- Prescribe net radiation and estimate the remaining radiative fluxes
[1853]	1194	DO i = nxlg, nxrg
	1195	DO j = nysg, nyng
	1196	k = nzb_s_inner(j,i)
	1197
	1198	rad_net(j,i) = net_radiation
[1976]	1199
	1200	exn1 = (hyp(k+1) / 100000.0_wp )**0.286_wp
	1201
	1202	IF ( cloud_physics ) THEN
	1203	pt1 = pt(k+1,j,i) + l_d_cp / exn1 * ql(k+1,j,i)
	1204	rad_lw_in(0,j,i) = 0.8_wp * sigma_sb * (pt1 * exn1)**4
	1205	ELSE
	1206	rad_lw_in(0,j,i) = 0.8_wp * sigma_sb * (pt(k+1,j,i) * exn1)**4
	1207	ENDIF
	1208
[1853]	1209	rad_lw_out(0,j,i) = emissivity * sigma_sb * (pt(k,j,i) * exn)**4
	1210
[1976]	1211	rad_sw_in(0,j,i) = ( rad_net(j,i) - rad_lw_in(0,j,i) &
	1212	+ rad_lw_out(0,j,i) ) &
	1213	/ ( 1.0_wp - alpha(j,i) )
	1214
[1853]	1215	ENDDO
	1216	ENDDO
	1217
	1218	END SUBROUTINE radiation_constant
	1219
	1220	!------------------------------------------------------------------------------!
	1221	! Description:
	1222	! ------------
[1826]	1223	!> Header output for radiation model
	1224	!------------------------------------------------------------------------------!
	1225	SUBROUTINE radiation_header ( io )
	1226
	1227
	1228	IMPLICIT NONE
	1229
	1230	INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file
	1231
	1232
	1233
	1234	!
	1235	!-- Write radiation model header
	1236	WRITE( io, 3 )
	1237
	1238	IF ( radiation_scheme == "constant" ) THEN
	1239	WRITE( io, 4 ) net_radiation
	1240	ELSEIF ( radiation_scheme == "clear-sky" ) THEN
	1241	WRITE( io, 5 )
	1242	ELSEIF ( radiation_scheme == "rrtmg" ) THEN
	1243	WRITE( io, 6 )
	1244	IF ( .NOT. lw_radiation ) WRITE( io, 10 )
	1245	IF ( .NOT. sw_radiation ) WRITE( io, 11 )
	1246	ENDIF
	1247
	1248	IF ( albedo_type == 0 ) THEN
	1249	WRITE( io, 7 ) albedo
	1250	ELSE
	1251	WRITE( io, 8 ) TRIM( albedo_type_name(albedo_type) )
	1252	ENDIF
	1253	IF ( constant_albedo ) THEN
	1254	WRITE( io, 9 )
	1255	ENDIF
	1256
	1257	IF ( radiation .AND. radiation_scheme /= 'constant' ) THEN
	1258	WRITE ( io, 1 ) lambda
	1259	WRITE ( io, 2 ) day_init, time_utc_init
	1260	ENDIF
	1261
	1262	WRITE( io, 12 ) dt_radiation
	1263
	1264
	1265	1 FORMAT (' Geograph. longitude : lambda = ',F4.1,' degr')
	1266	2 FORMAT (' Day of the year at model start : day_init = ',I3 &
	1267	/' UTC time at model start : time_utc_init = ',F7.1' s')
	1268	3 FORMAT (//' Radiation model information:'/ &
	1269	' ----------------------------'/)
	1270	4 FORMAT (' --> Using constant net radiation: net_radiation = ', F6.2, &
	1271	// 'W/m**2')
	1272	5 FORMAT (' --> Simple radiation scheme for clear sky is used (no clouds,', &
	1273	' default)')
	1274	6 FORMAT (' --> RRTMG scheme is used')
	1275	7 FORMAT (/' User-specific surface albedo: albedo =', F6.3)
	1276	8 FORMAT (/' Albedo is set for land surface type: ', A)
	1277	9 FORMAT (/' --> Albedo is fixed during the run')
	1278	10 FORMAT (/' --> Longwave radiation is disabled')
	1279	11 FORMAT (/' --> Shortwave radiation is disabled.')
	1280	12 FORMAT (' Timestep: dt_radiation = ', F6.2, ' s')
	1281
	1282
	1283	END SUBROUTINE radiation_header
	1284
	1285
	1286	!------------------------------------------------------------------------------!
	1287	! Description:
	1288	! ------------
	1289	!> Parin for &radiation_par for radiation model
	1290	!------------------------------------------------------------------------------!
	1291	SUBROUTINE radiation_parin
	1292
	1293
	1294	IMPLICIT NONE
	1295
	1296	CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file
	1297
	1298	NAMELIST /radiation_par/ albedo, albedo_type, albedo_lw_dir, &
	1299	albedo_lw_dif, albedo_sw_dir, albedo_sw_dif, &
	1300	constant_albedo, day_init, dt_radiation, &
	1301	lambda, lw_radiation, net_radiation, &
	1302	radiation_scheme, skip_time_do_radiation, &
	1303	sw_radiation, time_utc_init, &
	1304	unscheduled_radiation_calls
	1305
	1306	line = ' '
	1307
	1308	!
	1309	!-- Try to find radiation model package
	1310	REWIND ( 11 )
	1311	line = ' '
	1312	DO WHILE ( INDEX( line, '&radiation_par' ) == 0 )
	1313	READ ( 11, '(A)', END=10 ) line
	1314	ENDDO
	1315	BACKSPACE ( 11 )
	1316
	1317	!
	1318	!-- Read user-defined namelist
	1319	READ ( 11, radiation_par )
	1320
	1321	!
	1322	!-- Set flag that indicates that the radiation model is switched on
	1323	radiation = .TRUE.
	1324
	1325	10 CONTINUE
	1326
	1327
	1328	END SUBROUTINE radiation_parin
	1329
	1330
	1331	!------------------------------------------------------------------------------!
	1332	! Description:
	1333	! ------------
[1682]	1334	!> Implementation of the RRTMG radiation_scheme
[1585]	1335	!------------------------------------------------------------------------------!
	1336	SUBROUTINE radiation_rrtmg
	1337
	1338	USE indices, &
	1339	ONLY: nbgp
	1340
	1341	USE particle_attributes, &
	1342	ONLY: grid_particles, number_of_particles, particles, &
	1343	particle_advection_start, prt_count
	1344
	1345	IMPLICIT NONE
	1346
	1347	#if defined ( __rrtmg )
	1348
[1691]	1349	INTEGER(iwp) :: i, j, k, n !< loop indices
[1585]	1350
[1691]	1351	REAL(wp) :: s_r2, & !< weighted sum over all droplets with r^2
	1352	s_r3 !< weighted sum over all droplets with r^3
[1585]	1353
	1354	!
	1355	!-- Calculate current (cosine of) zenith angle and whether the sun is up
	1356	CALL calc_zenith
	1357	!
	1358	!-- Calculate surface albedo
	1359	IF ( .NOT. constant_albedo ) THEN
	1360	CALL calc_albedo
	1361	ENDIF
	1362
	1363	!
	1364	!-- Prepare input data for RRTMG
	1365
	1366	!
	1367	!-- In case of large scale forcing with surface data, calculate new pressure
	1368	!-- profile. nzt_rad might be modified by these calls and all required arrays
	1369	!-- will then be re-allocated
[1691]	1370	IF ( large_scale_forcing .AND. lsf_surf ) THEN
[1585]	1371	CALL read_sounding_data
	1372	CALL read_trace_gas_data
	1373	ENDIF
	1374	!
	1375	!-- Loop over all grid points
	1376	DO i = nxl, nxr
	1377	DO j = nys, nyn
	1378
	1379	!
	1380	!-- Prepare profiles of temperature and H2O volume mixing ratio
[1691]	1381	rrtm_tlev(0,nzb+1) = pt(nzb,j,i) * ( surface_pressure &
	1382	/ 1000.0_wp )**0.286_wp
[1585]	1383
	1384	DO k = nzb+1, nzt+1
	1385	rrtm_tlay(0,k) = pt(k,j,i) * ( (hyp(k) ) / 100000.0_wp &
[1691]	1386	)*0.286_wp + l_d_cp ql(k,j,i)
	1387	rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * (q(k,j,i) - ql(k,j,i))
[1585]	1388
	1389	ENDDO
	1390
	1391	!
	1392	!-- Avoid temperature/humidity jumps at the top of the LES domain by
	1393	!-- linear interpolation from nzt+2 to nzt+7
	1394	DO k = nzt+2, nzt+7
	1395	rrtm_tlay(0,k) = rrtm_tlay(0,nzt+1) &
	1396	+ ( rrtm_tlay(0,nzt+8) - rrtm_tlay(0,nzt+1) ) &
	1397	/ ( rrtm_play(0,nzt+8) - rrtm_play(0,nzt+1) ) &
	1398	* ( rrtm_play(0,k) - rrtm_play(0,nzt+1) )
	1399
	1400	rrtm_h2ovmr(0,k) = rrtm_h2ovmr(0,nzt+1) &
	1401	+ ( rrtm_h2ovmr(0,nzt+8) - rrtm_h2ovmr(0,nzt+1) )&
	1402	/ ( rrtm_play(0,nzt+8) - rrtm_play(0,nzt+1) )&
	1403	* ( rrtm_play(0,k) - rrtm_play(0,nzt+1) )
	1404
	1405	ENDDO
	1406
	1407	!-- Linear interpolate to zw grid
	1408	DO k = nzb+2, nzt+8
	1409	rrtm_tlev(0,k) = rrtm_tlay(0,k-1) + (rrtm_tlay(0,k) - &
	1410	rrtm_tlay(0,k-1)) &
	1411	/ ( rrtm_play(0,k) - rrtm_play(0,k-1) ) &
	1412	* ( rrtm_plev(0,k) - rrtm_play(0,k-1) )
	1413	ENDDO
	1414
	1415
	1416	!
	1417	!-- Calculate liquid water path and cloud fraction for each column.
	1418	!-- Note that LWP is required in g/mÂ² instead of kg/kg m.
	1419	rrtm_cldfr = 0.0_wp
	1420	rrtm_reliq = 0.0_wp
	1421	rrtm_cliqwp = 0.0_wp
[1691]	1422	rrtm_icld = 0
[1585]	1423
	1424	DO k = nzb+1, nzt+1
[1691]	1425	rrtm_cliqwp(0,k) = ql(k,j,i) * 1000.0_wp * &
	1426	(rrtm_plev(0,k) - rrtm_plev(0,k+1)) &
	1427	* 100.0_wp / g
[1585]	1428
[1691]	1429	IF ( rrtm_cliqwp(0,k) > 0.0_wp ) THEN
[1585]	1430	rrtm_cldfr(0,k) = 1.0_wp
[1691]	1431	IF ( rrtm_icld == 0 ) rrtm_icld = 1
[1585]	1432
	1433	!
	1434	!-- Calculate cloud droplet effective radius
	1435	IF ( cloud_physics ) THEN
[1691]	1436	rrtm_reliq(0,k) = 1.0E6_wp * ( 3.0_wp * ql(k,j,i) &
	1437	* rho_surface &
	1438	/ ( 4.0_wp * pi * nc_const * rho_l ) &
	1439	)**0.33333333333333_wp &
	1440	* EXP( LOG( sigma_gc )**2 )
[1585]	1441
	1442	ELSEIF ( cloud_droplets ) THEN
	1443	number_of_particles = prt_count(k,j,i)
	1444
	1445	IF (number_of_particles <= 0) CYCLE
	1446	particles => grid_particles(k,j,i)%particles(1:number_of_particles)
	1447	s_r2 = 0.0_wp
	1448	s_r3 = 0.0_wp
	1449
	1450	DO n = 1, number_of_particles
	1451	IF ( particles(n)%particle_mask ) THEN
	1452	s_r2 = s_r2 + particles(n)%radius*2 &
	1453	particles(n)%weight_factor
	1454	s_r3 = s_r3 + particles(n)%radius*3 &
	1455	particles(n)%weight_factor
	1456	ENDIF
	1457	ENDDO
	1458
	1459	IF ( s_r2 > 0.0_wp ) rrtm_reliq(0,k) = s_r3 / s_r2
	1460
	1461	ENDIF
	1462
	1463	!
	1464	!-- Limit effective radius
[1691]	1465	IF ( rrtm_reliq(0,k) > 0.0_wp ) THEN
[1585]	1466	rrtm_reliq(0,k) = MAX(rrtm_reliq(0,k),2.5_wp)
	1467	rrtm_reliq(0,k) = MIN(rrtm_reliq(0,k),60.0_wp)
	1468	ENDIF
	1469	ENDIF
	1470	ENDDO
	1471
	1472	!
	1473	!-- Set surface temperature
	1474	rrtm_tsfc = pt(nzb,j,i) * (surface_pressure / 1000.0_wp )**0.286_wp
	1475
	1476	IF ( lw_radiation ) THEN
	1477	CALL rrtmg_lw( 1, nzt_rad , rrtm_icld , rrtm_idrv ,&
	1478	rrtm_play , rrtm_plev , rrtm_tlay , rrtm_tlev ,&
	1479	rrtm_tsfc , rrtm_h2ovmr , rrtm_o3vmr , rrtm_co2vmr ,&
	1480	rrtm_ch4vmr , rrtm_n2ovmr , rrtm_o2vmr , rrtm_cfc11vmr ,&
	1481	rrtm_cfc12vmr , rrtm_cfc22vmr, rrtm_ccl4vmr , rrtm_emis ,&
	1482	rrtm_inflglw , rrtm_iceflglw, rrtm_liqflglw, rrtm_cldfr ,&
	1483	rrtm_lw_taucld , rrtm_cicewp , rrtm_cliqwp , rrtm_reice ,&
	1484	rrtm_reliq , rrtm_lw_tauaer, &
	1485	rrtm_lwuflx , rrtm_lwdflx , rrtm_lwhr , &
[1691]	1486	rrtm_lwuflxc , rrtm_lwdflxc , rrtm_lwhrc , &
	1487	rrtm_lwuflx_dt , rrtm_lwuflxc_dt )
[1585]	1488
[1691]	1489	!
	1490	!-- Save fluxes
[1585]	1491	DO k = nzb, nzt+1
	1492	rad_lw_in(k,j,i) = rrtm_lwdflx(0,k)
	1493	rad_lw_out(k,j,i) = rrtm_lwuflx(0,k)
	1494	ENDDO
	1495
[1691]	1496	!
	1497	!-- Save heating rates (convert from K/d to K/h)
	1498	DO k = nzb+1, nzt+1
	1499	rad_lw_hr(k,j,i) = rrtm_lwhr(0,k) * d_hours_day
	1500	rad_lw_cs_hr(k,j,i) = rrtm_lwhrc(0,k) * d_hours_day
	1501	ENDDO
[1585]	1502
[1709]	1503	!
	1504	!-- Save change in LW heating rate
	1505	rad_lw_out_change_0(j,i) = rrtm_lwuflx_dt(0,nzb)
	1506
[1585]	1507	ENDIF
	1508
	1509	IF ( sw_radiation .AND. sun_up ) THEN
	1510	CALL rrtmg_sw( 1, nzt_rad , rrtm_icld , rrtm_iaer ,&
	1511	rrtm_play , rrtm_plev , rrtm_tlay , rrtm_tlev ,&
	1512	rrtm_tsfc , rrtm_h2ovmr , rrtm_o3vmr , rrtm_co2vmr ,&
	1513	rrtm_ch4vmr , rrtm_n2ovmr , rrtm_o2vmr , rrtm_asdir(:,j,i),&
	1514	rrtm_asdif(:,j,i), rrtm_aldir(:,j,i), rrtm_aldif(:,j,i), zenith,&
	1515	0.0_wp , day , solar_constant, rrtm_inflgsw,&
	1516	rrtm_iceflgsw , rrtm_liqflgsw, rrtm_cldfr , rrtm_sw_taucld ,&
	1517	rrtm_sw_ssacld , rrtm_sw_asmcld, rrtm_sw_fsfcld, rrtm_cicewp ,&
	1518	rrtm_cliqwp , rrtm_reice , rrtm_reliq , rrtm_sw_tauaer ,&
	1519	rrtm_sw_ssaaer , rrtm_sw_asmaer , rrtm_sw_ecaer , &
	1520	rrtm_swuflx , rrtm_swdflx , rrtm_swhr , &
	1521	rrtm_swuflxc , rrtm_swdflxc , rrtm_swhrc )
	1522
[1691]	1523	!
	1524	!-- Save fluxes
[1585]	1525	DO k = nzb, nzt+1
	1526	rad_sw_in(k,j,i) = rrtm_swdflx(0,k)
	1527	rad_sw_out(k,j,i) = rrtm_swuflx(0,k)
	1528	ENDDO
[1691]	1529
	1530	!
	1531	!-- Save heating rates (convert from K/d to K/s)
	1532	DO k = nzb+1, nzt+1
	1533	rad_sw_hr(k,j,i) = rrtm_swhr(0,k) * d_hours_day
	1534	rad_sw_cs_hr(k,j,i) = rrtm_swhrc(0,k) * d_hours_day
	1535	ENDDO
	1536
[1585]	1537	ENDIF
	1538
	1539	!
	1540	!-- Calculate surface net radiation
	1541	rad_net(j,i) = rad_sw_in(nzb,j,i) - rad_sw_out(nzb,j,i) &
	1542	+ rad_lw_in(nzb,j,i) - rad_lw_out(nzb,j,i)
	1543
	1544	ENDDO
	1545	ENDDO
	1546
	1547	CALL exchange_horiz( rad_lw_in, nbgp )
	1548	CALL exchange_horiz( rad_lw_out, nbgp )
[1691]	1549	CALL exchange_horiz( rad_lw_hr, nbgp )
	1550	CALL exchange_horiz( rad_lw_cs_hr, nbgp )
	1551
[1585]	1552	CALL exchange_horiz( rad_sw_in, nbgp )
	1553	CALL exchange_horiz( rad_sw_out, nbgp )
[1691]	1554	CALL exchange_horiz( rad_sw_hr, nbgp )
	1555	CALL exchange_horiz( rad_sw_cs_hr, nbgp )
	1556
[1585]	1557	CALL exchange_horiz_2d( rad_net, nbgp )
[1709]	1558	CALL exchange_horiz_2d( rad_lw_out_change_0, nbgp )
[1585]	1559	#endif
	1560
	1561	END SUBROUTINE radiation_rrtmg
	1562
	1563
	1564	!------------------------------------------------------------------------------!
	1565	! Description:
	1566	! ------------
[1682]	1567	!> Calculate the cosine of the zenith angle (variable is called zenith)
[1585]	1568	!------------------------------------------------------------------------------!
	1569	SUBROUTINE calc_zenith
	1570
	1571	IMPLICIT NONE
	1572
[1682]	1573	REAL(wp) :: declination, & !< solar declination angle
	1574	hour_angle !< solar hour angle
[1585]	1575	!
[1496]	1576	!-- Calculate current day and time based on the initial values and simulation
	1577	!-- time
[1585]	1578	day = day_init + INT(FLOOR( (time_utc_init + time_since_reference_point) &
	1579	/ 86400.0_wp ), KIND=iwp)
[1496]	1580	time_utc = MOD((time_utc_init + time_since_reference_point), 86400.0_wp)
	1581
	1582
	1583	!
	1584	!-- Calculate solar declination and hour angle
[1585]	1585	declination = ASIN( decl_1 * SIN(decl_2 * REAL(day, KIND=wp) - decl_3) )
[1496]	1586	hour_angle = 2.0_wp * pi * (time_utc / 86400.0_wp) + lon - pi
	1587
	1588	!
[2007]	1589	!-- Calculate cosine of solar zenith angle
[1585]	1590	zenith(0) = SIN(lat) * SIN(declination) + COS(lat) * COS(declination) &
[1496]	1591	* COS(hour_angle)
[1585]	1592	zenith(0) = MAX(0.0_wp,zenith(0))
[1496]	1593
	1594	!
[2007]	1595	!-- Calculate solar directional vector
	1596	IF ( sun_direction ) THEN
	1597	!-- Direction in longitudes equals to sin(solar_azimuth) * sin(zenith)
	1598	sun_dir_lon(0) = -SIN(hour_angle) * COS(declination)
	1599	!-- Direction in latitues equals to cos(solar_azimuth) * sin(zenith)
	1600	sun_dir_lat(0) = SIN(declination) * COS(lat) - COS(hour_angle) &
	1601	* COS(declination) * SIN(lat)
	1602	ENDIF
	1603
	1604	!
[1585]	1605	!-- Check if the sun is up (otheriwse shortwave calculations can be skipped)
[1691]	1606	IF ( zenith(0) > 0.0_wp ) THEN
[1585]	1607	sun_up = .TRUE.
	1608	ELSE
	1609	sun_up = .FALSE.
	1610	END IF
[1496]	1611
[1585]	1612	END SUBROUTINE calc_zenith
	1613
[1606]	1614	#if defined ( __rrtmg ) && defined ( __netcdf )
[1585]	1615	!------------------------------------------------------------------------------!
	1616	! Description:
	1617	! ------------
[1682]	1618	!> Calculates surface albedo components based on Briegleb (1992) and
	1619	!> Briegleb et al. (1986)
[1585]	1620	!------------------------------------------------------------------------------!
	1621	SUBROUTINE calc_albedo
	1622
	1623	IMPLICIT NONE
	1624
	1625	IF ( sun_up ) THEN
[1496]	1626	!
[1585]	1627	!-- Ocean
	1628	IF ( albedo_type == 1 ) THEN
	1629	rrtm_aldir(0,:,:) = 0.026_wp / ( zenith(0)**1.7_wp + 0.065_wp ) &
	1630	+ 0.15_wp * ( zenith(0) - 0.1_wp ) &
	1631	* ( zenith(0) - 0.5_wp ) &
	1632	* ( zenith(0) - 1.0_wp )
	1633	rrtm_asdir(0,:,:) = rrtm_aldir(0,:,:)
	1634	!
	1635	!-- Snow
	1636	ELSEIF ( albedo_type == 16 ) THEN
	1637	IF ( zenith(0) < 0.5_wp ) THEN
	1638	rrtm_aldir(0,:,:) = 0.5_wp * (1.0_wp - aldif) &
	1639	* ( 3.0_wp / (1.0_wp + 4.0_wp &
	1640	* zenith(0))) - 1.0_wp
	1641	rrtm_asdir(0,:,:) = 0.5_wp * (1.0_wp - asdif) &
	1642	* ( 3.0_wp / (1.0_wp + 4.0_wp &
	1643	* zenith(0))) - 1.0_wp
[1496]	1644
[1585]	1645	rrtm_aldir(0,:,:) = MIN(0.98_wp, rrtm_aldir(0,:,:))
	1646	rrtm_asdir(0,:,:) = MIN(0.98_wp, rrtm_asdir(0,:,:))
	1647	ELSE
	1648	rrtm_aldir(0,:,:) = aldif
	1649	rrtm_asdir(0,:,:) = asdif
	1650	ENDIF
[1496]	1651	!
[1585]	1652	!-- Sea ice
	1653	ELSEIF ( albedo_type == 15 ) THEN
	1654	rrtm_aldir(0,:,:) = aldif
	1655	rrtm_asdir(0,:,:) = asdif
[1788]	1656
[1585]	1657	!
[1788]	1658	!-- Asphalt
	1659	ELSEIF ( albedo_type == 17 ) THEN
	1660	rrtm_aldir(0,:,:) = aldif
	1661	rrtm_asdir(0,:,:) = asdif
	1662	!
[1585]	1663	!-- Land surfaces
	1664	ELSE
	1665	SELECT CASE ( albedo_type )
[1496]	1666
[1585]	1667	!
	1668	!-- Surface types with strong zenith dependence
	1669	CASE ( 1, 2, 3, 4, 11, 12, 13 )
	1670	rrtm_aldir(0,:,:) = aldif * 1.4_wp / &
	1671	(1.0_wp + 0.8_wp * zenith(0))
	1672	rrtm_asdir(0,:,:) = asdif * 1.4_wp / &
	1673	(1.0_wp + 0.8_wp * zenith(0))
	1674	!
	1675	!-- Surface types with weak zenith dependence
	1676	CASE ( 5, 6, 7, 8, 9, 10, 14 )
	1677	rrtm_aldir(0,:,:) = aldif * 1.1_wp / &
	1678	(1.0_wp + 0.2_wp * zenith(0))
	1679	rrtm_asdir(0,:,:) = asdif * 1.1_wp / &
	1680	(1.0_wp + 0.2_wp * zenith(0))
[1496]	1681
[1585]	1682	CASE DEFAULT
	1683
	1684	END SELECT
	1685	ENDIF
	1686	!
	1687	!-- Diffusive albedo is taken from Table 2
	1688	rrtm_aldif(0,:,:) = aldif
	1689	rrtm_asdif(0,:,:) = asdif
	1690
	1691	ELSE
	1692
	1693	rrtm_aldir(0,:,:) = 0.0_wp
	1694	rrtm_asdir(0,:,:) = 0.0_wp
	1695	rrtm_aldif(0,:,:) = 0.0_wp
	1696	rrtm_asdif(0,:,:) = 0.0_wp
	1697	ENDIF
	1698	END SUBROUTINE calc_albedo
	1699
	1700	!------------------------------------------------------------------------------!
	1701	! Description:
	1702	! ------------
[1682]	1703	!> Read sounding data (pressure and temperature) from RADIATION_DATA.
[1585]	1704	!------------------------------------------------------------------------------!
	1705	SUBROUTINE read_sounding_data
	1706
	1707	IMPLICIT NONE
	1708
[1691]	1709	INTEGER(iwp) :: id, & !< NetCDF id of input file
	1710	id_dim_zrad, & !< pressure level id in the NetCDF file
	1711	id_var, & !< NetCDF variable id
	1712	k, & !< loop index
	1713	nz_snd, & !< number of vertical levels in the sounding data
	1714	nz_snd_start, & !< start vertical index for sounding data to be used
	1715	nz_snd_end !< end vertical index for souding data to be used
[1585]	1716
[1691]	1717	REAL(wp) :: t_surface !< actual surface temperature
[1585]	1718
[1691]	1719	REAL(wp), DIMENSION(:), ALLOCATABLE :: hyp_snd_tmp, & !< temporary hydrostatic pressure profile (sounding)
	1720	t_snd_tmp !< temporary temperature profile (sounding)
[1585]	1721
	1722	!
	1723	!-- In case of updates, deallocate arrays first (sufficient to check one
	1724	!-- array as the others are automatically allocated). This is required
	1725	!-- because nzt_rad might change during the update
	1726	IF ( ALLOCATED ( hyp_snd ) ) THEN
	1727	DEALLOCATE( hyp_snd )
	1728	DEALLOCATE( t_snd )
	1729	DEALLOCATE( q_snd )
	1730	DEALLOCATE ( rrtm_play )
	1731	DEALLOCATE ( rrtm_plev )
	1732	DEALLOCATE ( rrtm_tlay )
	1733	DEALLOCATE ( rrtm_tlev )
[1691]	1734
[1585]	1735	DEALLOCATE ( rrtm_h2ovmr )
	1736	DEALLOCATE ( rrtm_cicewp )
	1737	DEALLOCATE ( rrtm_cldfr )
	1738	DEALLOCATE ( rrtm_cliqwp )
	1739	DEALLOCATE ( rrtm_reice )
	1740	DEALLOCATE ( rrtm_reliq )
	1741	DEALLOCATE ( rrtm_lw_taucld )
	1742	DEALLOCATE ( rrtm_lw_tauaer )
[1691]	1743
[1585]	1744	DEALLOCATE ( rrtm_lwdflx )
[1691]	1745	DEALLOCATE ( rrtm_lwdflxc )
[1585]	1746	DEALLOCATE ( rrtm_lwuflx )
[1691]	1747	DEALLOCATE ( rrtm_lwuflxc )
	1748	DEALLOCATE ( rrtm_lwuflx_dt )
	1749	DEALLOCATE ( rrtm_lwuflxc_dt )
[1585]	1750	DEALLOCATE ( rrtm_lwhr )
	1751	DEALLOCATE ( rrtm_lwhrc )
[1691]	1752
[1585]	1753	DEALLOCATE ( rrtm_sw_taucld )
	1754	DEALLOCATE ( rrtm_sw_ssacld )
	1755	DEALLOCATE ( rrtm_sw_asmcld )
	1756	DEALLOCATE ( rrtm_sw_fsfcld )
	1757	DEALLOCATE ( rrtm_sw_tauaer )
	1758	DEALLOCATE ( rrtm_sw_ssaaer )
	1759	DEALLOCATE ( rrtm_sw_asmaer )
[1691]	1760	DEALLOCATE ( rrtm_sw_ecaer )
	1761
[1585]	1762	DEALLOCATE ( rrtm_swdflx )
[1691]	1763	DEALLOCATE ( rrtm_swdflxc )
[1585]	1764	DEALLOCATE ( rrtm_swuflx )
[1691]	1765	DEALLOCATE ( rrtm_swuflxc )
[1585]	1766	DEALLOCATE ( rrtm_swhr )
	1767	DEALLOCATE ( rrtm_swhrc )
[1691]	1768
[1585]	1769	ENDIF
	1770
	1771	!
	1772	!-- Open file for reading
	1773	nc_stat = NF90_OPEN( rrtm_input_file, NF90_NOWRITE, id )
[1783]	1774	CALL netcdf_handle_error_rad( 'read_sounding_data', 549 )
[1585]	1775
	1776	!
	1777	!-- Inquire dimension of z axis and save in nz_snd
	1778	nc_stat = NF90_INQ_DIMID( id, "Pressure", id_dim_zrad )
	1779	nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim_zrad, len = nz_snd )
[1783]	1780	CALL netcdf_handle_error_rad( 'read_sounding_data', 551 )
[1585]	1781
	1782	!
	1783	! !-- Allocate temporary array for storing pressure data
[1701]	1784	ALLOCATE( hyp_snd_tmp(1:nz_snd) )
[1585]	1785	hyp_snd_tmp = 0.0_wp
	1786
	1787
	1788	!-- Read pressure from file
	1789	nc_stat = NF90_INQ_VARID( id, "Pressure", id_var )
[1691]	1790	nc_stat = NF90_GET_VAR( id, id_var, hyp_snd_tmp(:), start = (/1/), &
[1585]	1791	count = (/nz_snd/) )
[1783]	1792	CALL netcdf_handle_error_rad( 'read_sounding_data', 552 )
[1585]	1793
	1794	!
	1795	!-- Allocate temporary array for storing temperature data
[1701]	1796	ALLOCATE( t_snd_tmp(1:nz_snd) )
[1585]	1797	t_snd_tmp = 0.0_wp
	1798
	1799	!
	1800	!-- Read temperature from file
	1801	nc_stat = NF90_INQ_VARID( id, "ReferenceTemperature", id_var )
[1691]	1802	nc_stat = NF90_GET_VAR( id, id_var, t_snd_tmp(:), start = (/1/), &
[1585]	1803	count = (/nz_snd/) )
[1783]	1804	CALL netcdf_handle_error_rad( 'read_sounding_data', 553 )
[1585]	1805
	1806	!
	1807	!-- Calculate start of sounding data
	1808	nz_snd_start = nz_snd + 1
[1701]	1809	nz_snd_end = nz_snd + 1
[1585]	1810
	1811	!
	1812	!-- Start filling vertical dimension at 10hPa above the model domain (hyp is
	1813	!-- in Pa, hyp_snd in hPa).
	1814	DO k = 1, nz_snd
[1691]	1815	IF ( hyp_snd_tmp(k) < ( hyp(nzt+1) - 1000.0_wp) * 0.01_wp ) THEN
[1585]	1816	nz_snd_start = k
	1817	EXIT
	1818	END IF
	1819	END DO
	1820
[1691]	1821	IF ( nz_snd_start <= nz_snd ) THEN
[1701]	1822	nz_snd_end = nz_snd
[1585]	1823	END IF
	1824
	1825
	1826	!
	1827	!-- Calculate of total grid points for RRTMG calculations
[1701]	1828	nzt_rad = nzt + nz_snd_end - nz_snd_start + 1
[1585]	1829
	1830	!
	1831	!-- Save data above LES domain in hyp_snd, t_snd and q_snd
	1832	!-- Note: q_snd_tmp is not calculated at the moment (dry residual atmosphere)
	1833	ALLOCATE( hyp_snd(nzb+1:nzt_rad) )
	1834	ALLOCATE( t_snd(nzb+1:nzt_rad) )
	1835	ALLOCATE( q_snd(nzb+1:nzt_rad) )
	1836	hyp_snd = 0.0_wp
	1837	t_snd = 0.0_wp
	1838	q_snd = 0.0_wp
	1839
[1757]	1840	hyp_snd(nzt+2:nzt_rad) = hyp_snd_tmp(nz_snd_start+1:nz_snd_end)
	1841	t_snd(nzt+2:nzt_rad) = t_snd_tmp(nz_snd_start+1:nz_snd_end)
[1585]	1842
	1843	nc_stat = NF90_CLOSE( id )
	1844
	1845	!
	1846	!-- Calculate pressure levels on zu and zw grid. Sounding data is added at
	1847	!-- top of the LES domain. This routine does not consider horizontal or
	1848	!-- vertical variability of pressure and temperature
	1849	ALLOCATE ( rrtm_play(0:0,nzb+1:nzt_rad+1) )
	1850	ALLOCATE ( rrtm_plev(0:0,nzb+1:nzt_rad+2) )
	1851
[1691]	1852	t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
[1585]	1853	DO k = nzb+1, nzt+1
	1854	rrtm_play(0,k) = hyp(k) * 0.01_wp
	1855	rrtm_plev(0,k) = surface_pressure * ( (t_surface - g/cp * zw(k-1)) / &
	1856	t_surface )**(1.0_wp/0.286_wp)
	1857	ENDDO
	1858
	1859	DO k = nzt+2, nzt_rad
	1860	rrtm_play(0,k) = hyp_snd(k)
	1861	rrtm_plev(0,k) = 0.5_wp * ( rrtm_play(0,k) + rrtm_play(0,k-1) )
	1862	ENDDO
	1863	rrtm_plev(0,nzt_rad+1) = MAX( 0.5 * hyp_snd(nzt_rad), &
	1864	1.5 * hyp_snd(nzt_rad) &
	1865	- 0.5 * hyp_snd(nzt_rad-1) )
	1866	rrtm_plev(0,nzt_rad+2) = MIN( 1.0E-4_wp, &
	1867	0.25_wp * rrtm_plev(0,nzt_rad+1) )
	1868
	1869	rrtm_play(0,nzt_rad+1) = 0.5 * rrtm_plev(0,nzt_rad+1)
	1870
	1871	!
	1872	!-- Calculate temperature/humidity levels at top of the LES domain.
	1873	!-- Currently, the temperature is taken from sounding data (might lead to a
	1874	!-- temperature jump at interface. To do: Humidity is currently not
	1875	!-- calculated above the LES domain.
	1876	ALLOCATE ( rrtm_tlay(0:0,nzb+1:nzt_rad+1) )
	1877	ALLOCATE ( rrtm_tlev(0:0,nzb+1:nzt_rad+2) )
	1878	ALLOCATE ( rrtm_h2ovmr(0:0,nzb+1:nzt_rad+1) )
	1879
	1880	DO k = nzt+8, nzt_rad
	1881	rrtm_tlay(0,k) = t_snd(k)
	1882	rrtm_h2ovmr(0,k) = q_snd(k)
	1883	ENDDO
[1691]	1884	rrtm_tlay(0,nzt_rad+1) = 2.0_wp * rrtm_tlay(0,nzt_rad) &
	1885	- rrtm_tlay(0,nzt_rad-1)
[1585]	1886	DO k = nzt+9, nzt_rad+1
	1887	rrtm_tlev(0,k) = rrtm_tlay(0,k-1) + (rrtm_tlay(0,k) &
	1888	- rrtm_tlay(0,k-1)) &
	1889	/ ( rrtm_play(0,k) - rrtm_play(0,k-1) ) &
	1890	* ( rrtm_plev(0,k) - rrtm_play(0,k-1) )
	1891	ENDDO
	1892	rrtm_h2ovmr(0,nzt_rad+1) = rrtm_h2ovmr(0,nzt_rad)
	1893
	1894	rrtm_tlev(0,nzt_rad+2) = 2.0_wp * rrtm_tlay(0,nzt_rad+1) &
	1895	- rrtm_tlev(0,nzt_rad)
	1896	!
	1897	!-- Allocate remaining RRTMG arrays
	1898	ALLOCATE ( rrtm_cicewp(0:0,nzb+1:nzt_rad+1) )
	1899	ALLOCATE ( rrtm_cldfr(0:0,nzb+1:nzt_rad+1) )
	1900	ALLOCATE ( rrtm_cliqwp(0:0,nzb+1:nzt_rad+1) )
	1901	ALLOCATE ( rrtm_reice(0:0,nzb+1:nzt_rad+1) )
	1902	ALLOCATE ( rrtm_reliq(0:0,nzb+1:nzt_rad+1) )
	1903	ALLOCATE ( rrtm_lw_taucld(1:nbndlw+1,0:0,nzb+1:nzt_rad+1) )
	1904	ALLOCATE ( rrtm_lw_tauaer(0:0,nzb+1:nzt_rad+1,1:nbndlw+1) )
	1905	ALLOCATE ( rrtm_sw_taucld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
	1906	ALLOCATE ( rrtm_sw_ssacld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
	1907	ALLOCATE ( rrtm_sw_asmcld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
	1908	ALLOCATE ( rrtm_sw_fsfcld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
	1909	ALLOCATE ( rrtm_sw_tauaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) )
	1910	ALLOCATE ( rrtm_sw_ssaaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) )
	1911	ALLOCATE ( rrtm_sw_asmaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) )
	1912	ALLOCATE ( rrtm_sw_ecaer(0:0,nzb+1:nzt_rad+1,1:naerec+1) )
	1913
	1914	!
	1915	!-- The ice phase is currently not considered in PALM
	1916	rrtm_cicewp = 0.0_wp
	1917	rrtm_reice = 0.0_wp
	1918
	1919	!
	1920	!-- Set other parameters (move to NAMELIST parameters in the future)
	1921	rrtm_lw_tauaer = 0.0_wp
	1922	rrtm_lw_taucld = 0.0_wp
	1923	rrtm_sw_taucld = 0.0_wp
	1924	rrtm_sw_ssacld = 0.0_wp
	1925	rrtm_sw_asmcld = 0.0_wp
	1926	rrtm_sw_fsfcld = 0.0_wp
	1927	rrtm_sw_tauaer = 0.0_wp
	1928	rrtm_sw_ssaaer = 0.0_wp
	1929	rrtm_sw_asmaer = 0.0_wp
	1930	rrtm_sw_ecaer = 0.0_wp
	1931
	1932
	1933	ALLOCATE ( rrtm_swdflx(0:0,nzb:nzt_rad+1) )
	1934	ALLOCATE ( rrtm_swuflx(0:0,nzb:nzt_rad+1) )
	1935	ALLOCATE ( rrtm_swhr(0:0,nzb+1:nzt_rad+1) )
	1936	ALLOCATE ( rrtm_swuflxc(0:0,nzb:nzt_rad+1) )
	1937	ALLOCATE ( rrtm_swdflxc(0:0,nzb:nzt_rad+1) )
	1938	ALLOCATE ( rrtm_swhrc(0:0,nzb+1:nzt_rad+1) )
	1939
	1940	rrtm_swdflx = 0.0_wp
	1941	rrtm_swuflx = 0.0_wp
	1942	rrtm_swhr = 0.0_wp
	1943	rrtm_swuflxc = 0.0_wp
	1944	rrtm_swdflxc = 0.0_wp
	1945	rrtm_swhrc = 0.0_wp
	1946
	1947	ALLOCATE ( rrtm_lwdflx(0:0,nzb:nzt_rad+1) )
	1948	ALLOCATE ( rrtm_lwuflx(0:0,nzb:nzt_rad+1) )
	1949	ALLOCATE ( rrtm_lwhr(0:0,nzb+1:nzt_rad+1) )
	1950	ALLOCATE ( rrtm_lwuflxc(0:0,nzb:nzt_rad+1) )
	1951	ALLOCATE ( rrtm_lwdflxc(0:0,nzb:nzt_rad+1) )
	1952	ALLOCATE ( rrtm_lwhrc(0:0,nzb+1:nzt_rad+1) )
	1953
	1954	rrtm_lwdflx = 0.0_wp
	1955	rrtm_lwuflx = 0.0_wp
	1956	rrtm_lwhr = 0.0_wp
	1957	rrtm_lwuflxc = 0.0_wp
	1958	rrtm_lwdflxc = 0.0_wp
	1959	rrtm_lwhrc = 0.0_wp
	1960
[1691]	1961	ALLOCATE ( rrtm_lwuflx_dt(0:0,nzb:nzt_rad+1) )
	1962	ALLOCATE ( rrtm_lwuflxc_dt(0:0,nzb:nzt_rad+1) )
[1585]	1963
[1709]	1964	rrtm_lwuflx_dt = 0.0_wp
[1691]	1965	rrtm_lwuflxc_dt = 0.0_wp
	1966
[1585]	1967	END SUBROUTINE read_sounding_data
	1968
	1969
	1970	!------------------------------------------------------------------------------!
	1971	! Description:
	1972	! ------------
[1682]	1973	!> Read trace gas data from file
[1585]	1974	!------------------------------------------------------------------------------!
	1975	SUBROUTINE read_trace_gas_data
	1976
	1977	USE rrsw_ncpar
	1978
	1979	IMPLICIT NONE
	1980
[1691]	1981	INTEGER(iwp), PARAMETER :: num_trace_gases = 9 !< number of trace gases (absorbers)
[1585]	1982
[1691]	1983	CHARACTER(LEN=5), DIMENSION(num_trace_gases), PARAMETER :: & !< trace gas names
[1585]	1984	trace_names = (/'O3 ', 'CO2 ', 'CH4 ', 'N2O ', 'O2 ', &
	1985	'CFC11', 'CFC12', 'CFC22', 'CCL4 '/)
	1986
[1691]	1987	INTEGER(iwp) :: id, & !< NetCDF id
	1988	k, & !< loop index
	1989	m, & !< loop index
	1990	n, & !< loop index
	1991	nabs, & !< number of absorbers
	1992	np, & !< number of pressure levels
	1993	id_abs, & !< NetCDF id of the respective absorber
	1994	id_dim, & !< NetCDF id of asborber's dimension
	1995	id_var !< NetCDf id ot the absorber
[1585]	1996
	1997	REAL(wp) :: p_mls_l, p_mls_u, p_wgt_l, p_wgt_u, p_mls_m
	1998
	1999
[1682]	2000	REAL(wp), DIMENSION(:), ALLOCATABLE :: p_mls, & !< pressure levels for the absorbers
	2001	rrtm_play_tmp, & !< temporary array for pressure zu-levels
	2002	rrtm_plev_tmp, & !< temporary array for pressure zw-levels
	2003	trace_path_tmp !< temporary array for storing trace gas path data
[1585]	2004
[1682]	2005	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: trace_mls, & !< array for storing the absorber amounts
	2006	trace_mls_path, & !< array for storing trace gas path data
	2007	trace_mls_tmp !< temporary array for storing trace gas data
[1585]	2008
	2009
	2010	!
	2011	!-- In case of updates, deallocate arrays first (sufficient to check one
	2012	!-- array as the others are automatically allocated)
	2013	IF ( ALLOCATED ( rrtm_o3vmr ) ) THEN
	2014	DEALLOCATE ( rrtm_o3vmr )
	2015	DEALLOCATE ( rrtm_co2vmr )
	2016	DEALLOCATE ( rrtm_ch4vmr )
	2017	DEALLOCATE ( rrtm_n2ovmr )
	2018	DEALLOCATE ( rrtm_o2vmr )
	2019	DEALLOCATE ( rrtm_cfc11vmr )
	2020	DEALLOCATE ( rrtm_cfc12vmr )
	2021	DEALLOCATE ( rrtm_cfc22vmr )
	2022	DEALLOCATE ( rrtm_ccl4vmr )
	2023	ENDIF
	2024
	2025	!
	2026	!-- Allocate trace gas profiles
	2027	ALLOCATE ( rrtm_o3vmr(0:0,1:nzt_rad+1) )
	2028	ALLOCATE ( rrtm_co2vmr(0:0,1:nzt_rad+1) )
	2029	ALLOCATE ( rrtm_ch4vmr(0:0,1:nzt_rad+1) )
	2030	ALLOCATE ( rrtm_n2ovmr(0:0,1:nzt_rad+1) )
	2031	ALLOCATE ( rrtm_o2vmr(0:0,1:nzt_rad+1) )
	2032	ALLOCATE ( rrtm_cfc11vmr(0:0,1:nzt_rad+1) )
	2033	ALLOCATE ( rrtm_cfc12vmr(0:0,1:nzt_rad+1) )
	2034	ALLOCATE ( rrtm_cfc22vmr(0:0,1:nzt_rad+1) )
	2035	ALLOCATE ( rrtm_ccl4vmr(0:0,1:nzt_rad+1) )
	2036
	2037	!
	2038	!-- Open file for reading
	2039	nc_stat = NF90_OPEN( rrtm_input_file, NF90_NOWRITE, id )
[1783]	2040	CALL netcdf_handle_error_rad( 'read_trace_gas_data', 549 )
[1585]	2041	!
	2042	!-- Inquire dimension ids and dimensions
	2043	nc_stat = NF90_INQ_DIMID( id, "Pressure", id_dim )
[1783]	2044	CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
[1585]	2045	nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, len = np)
[1783]	2046	CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
[1585]	2047
	2048	nc_stat = NF90_INQ_DIMID( id, "Absorber", id_dim )
[1783]	2049	CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
[1585]	2050	nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, len = nabs )
[1783]	2051	CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
[1585]	2052
	2053
	2054	!
	2055	!-- Allocate pressure, and trace gas arrays
	2056	ALLOCATE( p_mls(1:np) )
	2057	ALLOCATE( trace_mls(1:num_trace_gases,1:np) )
	2058	ALLOCATE( trace_mls_tmp(1:nabs,1:np) )
	2059
	2060
	2061	nc_stat = NF90_INQ_VARID( id, "Pressure", id_var )
[1783]	2062	CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
[1585]	2063	nc_stat = NF90_GET_VAR( id, id_var, p_mls )
[1783]	2064	CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
[1585]	2065
	2066	nc_stat = NF90_INQ_VARID( id, "AbsorberAmountMLS", id_var )
[1783]	2067	CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
[1585]	2068	nc_stat = NF90_GET_VAR( id, id_var, trace_mls_tmp )
[1783]	2069	CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
[1585]	2070
	2071
	2072	!
	2073	!-- Write absorber amounts (mls) to trace_mls
	2074	DO n = 1, num_trace_gases
	2075	CALL getAbsorberIndex( TRIM( trace_names(n) ), id_abs )
	2076
	2077	trace_mls(n,1:np) = trace_mls_tmp(id_abs,1:np)
	2078
	2079	!
	2080	!-- Replace missing values by zero
	2081	WHERE ( trace_mls(n,:) > 2.0_wp )
	2082	trace_mls(n,:) = 0.0_wp
	2083	END WHERE
	2084	END DO
	2085
	2086	DEALLOCATE ( trace_mls_tmp )
	2087
	2088	nc_stat = NF90_CLOSE( id )
[1783]	2089	CALL netcdf_handle_error_rad( 'read_trace_gas_data', 551 )
[1585]	2090
	2091	!
	2092	!-- Add extra pressure level for calculations of the trace gas paths
	2093	ALLOCATE ( rrtm_play_tmp(1:nzt_rad+1) )
	2094	ALLOCATE ( rrtm_plev_tmp(1:nzt_rad+2) )
	2095
	2096	rrtm_play_tmp(1:nzt_rad) = rrtm_play(0,1:nzt_rad)
	2097	rrtm_plev_tmp(1:nzt_rad+1) = rrtm_plev(0,1:nzt_rad+1)
	2098	rrtm_play_tmp(nzt_rad+1) = rrtm_plev(0,nzt_rad+1) * 0.5_wp
	2099	rrtm_plev_tmp(nzt_rad+2) = MIN( 1.0E-4_wp, 0.25_wp &
	2100	* rrtm_plev(0,nzt_rad+1) )
	2101
	2102	!
	2103	!-- Calculate trace gas path (zero at surface) with interpolation to the
	2104	!-- sounding levels
	2105	ALLOCATE ( trace_mls_path(1:nzt_rad+2,1:num_trace_gases) )
	2106
	2107	trace_mls_path(nzb+1,:) = 0.0_wp
	2108
	2109	DO k = nzb+2, nzt_rad+2
	2110	DO m = 1, num_trace_gases
	2111	trace_mls_path(k,m) = trace_mls_path(k-1,m)
	2112
	2113	!
	2114	!-- When the pressure level is higher than the trace gas pressure
	2115	!-- level, assume that
[1691]	2116	IF ( rrtm_plev_tmp(k-1) > p_mls(1) ) THEN
[1585]	2117
	2118	trace_mls_path(k,m) = trace_mls_path(k,m) + trace_mls(m,1) &
	2119	* ( rrtm_plev_tmp(k-1) &
	2120	- MAX( p_mls(1), rrtm_plev_tmp(k) ) &
	2121	) / g
	2122	ENDIF
	2123
	2124	!
	2125	!-- Integrate for each sounding level from the contributing p_mls
	2126	!-- levels
	2127	DO n = 2, np
	2128	!
	2129	!-- Limit p_mls so that it is within the model level
	2130	p_mls_u = MIN( rrtm_plev_tmp(k-1), &
	2131	MAX( rrtm_plev_tmp(k), p_mls(n) ) )
	2132	p_mls_l = MIN( rrtm_plev_tmp(k-1), &
	2133	MAX( rrtm_plev_tmp(k), p_mls(n-1) ) )
	2134
[1691]	2135	IF ( p_mls_l > p_mls_u ) THEN
[1585]	2136
	2137	!
	2138	!-- Calculate weights for interpolation
	2139	p_mls_m = 0.5_wp * (p_mls_l + p_mls_u)
	2140	p_wgt_u = (p_mls(n-1) - p_mls_m) / (p_mls(n-1) - p_mls(n))
	2141	p_wgt_l = (p_mls_m - p_mls(n)) / (p_mls(n-1) - p_mls(n))
	2142
	2143	!
	2144	!-- Add level to trace gas path
	2145	trace_mls_path(k,m) = trace_mls_path(k,m) &
	2146	+ ( p_wgt_u * trace_mls(m,n) &
	2147	+ p_wgt_l * trace_mls(m,n-1) ) &
[1691]	2148	* (p_mls_l - p_mls_u) / g
[1585]	2149	ENDIF
	2150	ENDDO
	2151
[1691]	2152	IF ( rrtm_plev_tmp(k) < p_mls(np) ) THEN
[1585]	2153	trace_mls_path(k,m) = trace_mls_path(k,m) + trace_mls(m,np) &
	2154	* ( MIN( rrtm_plev_tmp(k-1), p_mls(np) ) &
	2155	- rrtm_plev_tmp(k) &
	2156	) / g
	2157	ENDIF
[1496]	2158	ENDDO
	2159	ENDDO
	2160
	2161
[1585]	2162	!
	2163	!-- Prepare trace gas path profiles
	2164	ALLOCATE ( trace_path_tmp(1:nzt_rad+1) )
[1496]	2165
[1585]	2166	DO m = 1, num_trace_gases
	2167
	2168	trace_path_tmp(1:nzt_rad+1) = ( trace_mls_path(2:nzt_rad+2,m) &
	2169	- trace_mls_path(1:nzt_rad+1,m) ) * g &
	2170	/ ( rrtm_plev_tmp(1:nzt_rad+1) &
	2171	- rrtm_plev_tmp(2:nzt_rad+2) )
	2172
	2173	!
	2174	!-- Save trace gas paths to the respective arrays
	2175	SELECT CASE ( TRIM( trace_names(m) ) )
	2176
	2177	CASE ( 'O3' )
	2178
	2179	rrtm_o3vmr(0,:) = trace_path_tmp(:)
	2180
	2181	CASE ( 'CO2' )
	2182
	2183	rrtm_co2vmr(0,:) = trace_path_tmp(:)
	2184
	2185	CASE ( 'CH4' )
	2186
	2187	rrtm_ch4vmr(0,:) = trace_path_tmp(:)
	2188
	2189	CASE ( 'N2O' )
	2190
	2191	rrtm_n2ovmr(0,:) = trace_path_tmp(:)
	2192
	2193	CASE ( 'O2' )
	2194
	2195	rrtm_o2vmr(0,:) = trace_path_tmp(:)
	2196
	2197	CASE ( 'CFC11' )
	2198
	2199	rrtm_cfc11vmr(0,:) = trace_path_tmp(:)
	2200
	2201	CASE ( 'CFC12' )
	2202
	2203	rrtm_cfc12vmr(0,:) = trace_path_tmp(:)
	2204
	2205	CASE ( 'CFC22' )
	2206
	2207	rrtm_cfc22vmr(0,:) = trace_path_tmp(:)
	2208
	2209	CASE ( 'CCL4' )
	2210
	2211	rrtm_ccl4vmr(0,:) = trace_path_tmp(:)
	2212
	2213	CASE DEFAULT
	2214
	2215	END SELECT
	2216
	2217	ENDDO
	2218
	2219	DEALLOCATE ( trace_path_tmp )
	2220	DEALLOCATE ( trace_mls_path )
	2221	DEALLOCATE ( rrtm_play_tmp )
	2222	DEALLOCATE ( rrtm_plev_tmp )
	2223	DEALLOCATE ( trace_mls )
	2224	DEALLOCATE ( p_mls )
	2225
	2226	END SUBROUTINE read_trace_gas_data
	2227
[1826]	2228
[1783]	2229	SUBROUTINE netcdf_handle_error_rad( routine_name, errno )
	2230
	2231	USE control_parameters, &
	2232	ONLY: message_string
	2233
	2234	USE NETCDF
	2235
	2236	USE pegrid
	2237
	2238	IMPLICIT NONE
	2239
	2240	CHARACTER(LEN=6) :: message_identifier
	2241	CHARACTER(LEN=*) :: routine_name
	2242
	2243	INTEGER(iwp) :: errno
	2244
	2245	IF ( nc_stat /= NF90_NOERR ) THEN
	2246
	2247	WRITE( message_identifier, '(''NC'',I4.4)' ) errno
	2248	message_string = TRIM( NF90_STRERROR( nc_stat ) )
	2249
	2250	CALL message( routine_name, message_identifier, 2, 2, 0, 6, 1 )
	2251
	2252	ENDIF
	2253
	2254	END SUBROUTINE netcdf_handle_error_rad
[1585]	2255	#endif
	2256
	2257
[1551]	2258	!------------------------------------------------------------------------------!
	2259	! Description:
	2260	! ------------
[1682]	2261	!> Calculate temperature tendency due to radiative cooling/heating.
	2262	!> Cache-optimized version.
[1551]	2263	!------------------------------------------------------------------------------!
[1976]	2264	SUBROUTINE radiation_tendency_ij ( i, j, tend )
[1496]	2265
[1976]	2266	USE cloud_parameters, &
	2267	ONLY: pt_d_t
[1551]	2268
[1976]	2269	IMPLICIT NONE
[1585]	2270
[1976]	2271	INTEGER(iwp) :: i, j, k !< loop indices
[1585]	2272
[1976]	2273	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend !< pt tendency term
[1585]	2274
[1976]	2275	IF ( radiation_scheme == 'rrtmg' ) THEN
	2276	#if defined ( __rrtmg )
[1585]	2277	!
[1691]	2278	!-- Calculate tendency based on heating rate
[1585]	2279	DO k = nzb+1, nzt+1
[1691]	2280	tend(k,j,i) = tend(k,j,i) + (rad_lw_hr(k,j,i) + rad_sw_hr(k,j,i)) &
[1976]	2281	* pt_d_t(k) * d_seconds_hour
[1585]	2282	ENDDO
	2283	#endif
[1976]	2284	ENDIF
[1585]	2285
	2286	END SUBROUTINE radiation_tendency_ij
	2287
	2288
[1551]	2289	!------------------------------------------------------------------------------!
	2290	! Description:
	2291	! ------------
[1682]	2292	!> Calculate temperature tendency due to radiative cooling/heating.
	2293	!> Vector-optimized version
[1551]	2294	!------------------------------------------------------------------------------!
[1976]	2295	SUBROUTINE radiation_tendency ( tend )
[1551]	2296
[1976]	2297	USE cloud_parameters, &
	2298	ONLY: pt_d_t
[1551]	2299
[1976]	2300	USE indices, &
	2301	ONLY: nxl, nxr, nyn, nys
[1585]	2302
[1976]	2303	IMPLICIT NONE
[1585]	2304
[1976]	2305	INTEGER(iwp) :: i, j, k !< loop indices
[1585]	2306
[1976]	2307	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend !< pt tendency term
[1585]	2308
[1976]	2309	IF ( radiation_scheme == 'rrtmg' ) THEN
	2310	#if defined ( __rrtmg )
[1691]	2311	!
	2312	!-- Calculate tendency based on heating rate
[1585]	2313	DO i = nxl, nxr
	2314	DO j = nys, nyn
	2315	DO k = nzb+1, nzt+1
[1691]	2316	tend(k,j,i) = tend(k,j,i) + ( rad_lw_hr(k,j,i) &
[1976]	2317	+ rad_sw_hr(k,j,i) ) * pt_d_t(k) &
	2318	* d_seconds_hour
[1585]	2319	ENDDO
[1976]	2320	ENDDO
[1585]	2321	ENDDO
	2322	#endif
[1976]	2323	ENDIF
[1585]	2324
	2325
[1976]	2326	END SUBROUTINE radiation_tendency
	2327
	2328	!------------------------------------------------------------------------------!
	2329	!
	2330	! Description:
	2331	! ------------
	2332	!> Subroutine for averaging 3D data
	2333	!------------------------------------------------------------------------------!
	2334	SUBROUTINE radiation_3d_data_averaging( mode, variable )
	2335
	2336
	2337	USE control_parameters
	2338
	2339	USE indices
	2340
	2341	USE kinds
	2342
	2343	IMPLICIT NONE
	2344
	2345	CHARACTER (LEN=*) :: mode !<
	2346	CHARACTER (LEN=*) :: variable !<
	2347
	2348	INTEGER(iwp) :: i !<
	2349	INTEGER(iwp) :: j !<
	2350	INTEGER(iwp) :: k !<
	2351
	2352	IF ( mode == 'allocate' ) THEN
	2353
	2354	SELECT CASE ( TRIM( variable ) )
	2355
	2356	CASE ( 'rad_net*' )
	2357	IF ( .NOT. ALLOCATED( rad_net_av ) ) THEN
	2358	ALLOCATE( rad_net_av(nysg:nyng,nxlg:nxrg) )
	2359	ENDIF
	2360	rad_net_av = 0.0_wp
	2361
	2362	CASE ( 'rad_lw_in' )
	2363	IF ( .NOT. ALLOCATED( rad_lw_in_av ) ) THEN
	2364	ALLOCATE( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	2365	ENDIF
	2366	rad_lw_in_av = 0.0_wp
	2367
	2368	CASE ( 'rad_lw_out' )
	2369	IF ( .NOT. ALLOCATED( rad_lw_out_av ) ) THEN
	2370	ALLOCATE( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	2371	ENDIF
	2372	rad_lw_out_av = 0.0_wp
	2373
	2374	CASE ( 'rad_lw_cs_hr' )
	2375	IF ( .NOT. ALLOCATED( rad_lw_cs_hr_av ) ) THEN
	2376	ALLOCATE( rad_lw_cs_hr_av(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
	2377	ENDIF
	2378	rad_lw_cs_hr_av = 0.0_wp
	2379
	2380	CASE ( 'rad_lw_hr' )
	2381	IF ( .NOT. ALLOCATED( rad_lw_hr_av ) ) THEN
	2382	ALLOCATE( rad_lw_hr_av(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
	2383	ENDIF
	2384	rad_lw_hr_av = 0.0_wp
	2385
	2386	CASE ( 'rad_sw_in' )
	2387	IF ( .NOT. ALLOCATED( rad_sw_in_av ) ) THEN
	2388	ALLOCATE( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	2389	ENDIF
	2390	rad_sw_in_av = 0.0_wp
	2391
	2392	CASE ( 'rad_sw_out' )
	2393	IF ( .NOT. ALLOCATED( rad_sw_out_av ) ) THEN
	2394	ALLOCATE( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	2395	ENDIF
	2396	rad_sw_out_av = 0.0_wp
	2397
	2398	CASE ( 'rad_sw_cs_hr' )
	2399	IF ( .NOT. ALLOCATED( rad_sw_cs_hr_av ) ) THEN
	2400	ALLOCATE( rad_sw_cs_hr_av(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
	2401	ENDIF
	2402	rad_sw_cs_hr_av = 0.0_wp
	2403
	2404	CASE ( 'rad_sw_hr' )
	2405	IF ( .NOT. ALLOCATED( rad_sw_hr_av ) ) THEN
	2406	ALLOCATE( rad_sw_hr_av(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
	2407	ENDIF
	2408	rad_sw_hr_av = 0.0_wp
	2409
	2410	CASE DEFAULT
	2411	CONTINUE
	2412
	2413	END SELECT
	2414
	2415	ELSEIF ( mode == 'sum' ) THEN
	2416
	2417	SELECT CASE ( TRIM( variable ) )
	2418
	2419	CASE ( 'rad_net*' )
	2420	DO i = nxlg, nxrg
	2421	DO j = nysg, nyng
	2422	rad_net_av(j,i) = rad_net_av(j,i) + rad_net(j,i)
	2423	ENDDO
	2424	ENDDO
	2425
	2426	CASE ( 'rad_lw_in' )
	2427	DO i = nxlg, nxrg
	2428	DO j = nysg, nyng
	2429	DO k = nzb, nzt+1
	2430	rad_lw_in_av(k,j,i) = rad_lw_in_av(k,j,i) + rad_lw_in(k,j,i)
	2431	ENDDO
	2432	ENDDO
	2433	ENDDO
	2434
	2435	CASE ( 'rad_lw_out' )
	2436	DO i = nxlg, nxrg
	2437	DO j = nysg, nyng
	2438	DO k = nzb, nzt+1
	2439	rad_lw_out_av(k,j,i) = rad_lw_out_av(k,j,i) + rad_lw_out(k,j,i)
	2440	ENDDO
	2441	ENDDO
	2442	ENDDO
	2443
	2444	CASE ( 'rad_lw_cs_hr' )
	2445	DO i = nxlg, nxrg
	2446	DO j = nysg, nyng
	2447	DO k = nzb, nzt+1
	2448	rad_lw_cs_hr_av(k,j,i) = rad_lw_cs_hr_av(k,j,i) + rad_lw_cs_hr(k,j,i)
	2449	ENDDO
	2450	ENDDO
	2451	ENDDO
	2452
	2453	CASE ( 'rad_lw_hr' )
	2454	DO i = nxlg, nxrg
	2455	DO j = nysg, nyng
	2456	DO k = nzb, nzt+1
	2457	rad_lw_hr_av(k,j,i) = rad_lw_hr_av(k,j,i) + rad_lw_hr(k,j,i)
	2458	ENDDO
	2459	ENDDO
	2460	ENDDO
	2461
	2462	CASE ( 'rad_sw_in' )
	2463	DO i = nxlg, nxrg
	2464	DO j = nysg, nyng
	2465	DO k = nzb, nzt+1
	2466	rad_sw_in_av(k,j,i) = rad_sw_in_av(k,j,i) + rad_sw_in(k,j,i)
	2467	ENDDO
	2468	ENDDO
	2469	ENDDO
	2470
	2471	CASE ( 'rad_sw_out' )
	2472	DO i = nxlg, nxrg
	2473	DO j = nysg, nyng
	2474	DO k = nzb, nzt+1
	2475	rad_sw_out_av(k,j,i) = rad_sw_out_av(k,j,i) + rad_sw_out(k,j,i)
	2476	ENDDO
	2477	ENDDO
	2478	ENDDO
	2479
	2480	CASE ( 'rad_sw_cs_hr' )
	2481	DO i = nxlg, nxrg
	2482	DO j = nysg, nyng
	2483	DO k = nzb, nzt+1
	2484	rad_sw_cs_hr_av(k,j,i) = rad_sw_cs_hr_av(k,j,i) + rad_sw_cs_hr(k,j,i)
	2485	ENDDO
	2486	ENDDO
	2487	ENDDO
	2488
	2489	CASE ( 'rad_sw_hr' )
	2490	DO i = nxlg, nxrg
	2491	DO j = nysg, nyng
	2492	DO k = nzb, nzt+1
	2493	rad_sw_hr_av(k,j,i) = rad_sw_hr_av(k,j,i) + rad_sw_hr(k,j,i)
	2494	ENDDO
	2495	ENDDO
	2496	ENDDO
	2497
	2498	CASE DEFAULT
	2499	CONTINUE
	2500
	2501	END SELECT
	2502
	2503	ELSEIF ( mode == 'average' ) THEN
	2504
	2505	SELECT CASE ( TRIM( variable ) )
	2506
	2507	CASE ( 'rad_net*' )
	2508	DO i = nxlg, nxrg
	2509	DO j = nysg, nyng
	2510	rad_net_av(j,i) = rad_net_av(j,i) / REAL( average_count_3d, KIND=wp )
	2511	ENDDO
	2512	ENDDO
	2513
	2514	CASE ( 'rad_lw_in' )
	2515	DO i = nxlg, nxrg
	2516	DO j = nysg, nyng
	2517	DO k = nzb, nzt+1
	2518	rad_lw_in_av(k,j,i) = rad_lw_in_av(k,j,i) / REAL( average_count_3d, KIND=wp )
	2519	ENDDO
	2520	ENDDO
	2521	ENDDO
	2522
	2523	CASE ( 'rad_lw_out' )
	2524	DO i = nxlg, nxrg
	2525	DO j = nysg, nyng
	2526	DO k = nzb, nzt+1
	2527	rad_lw_out_av(k,j,i) = rad_lw_out_av(k,j,i) / REAL( average_count_3d, KIND=wp )
	2528	ENDDO
	2529	ENDDO
	2530	ENDDO
	2531
	2532	CASE ( 'rad_lw_cs_hr' )
	2533	DO i = nxlg, nxrg
	2534	DO j = nysg, nyng
	2535	DO k = nzb, nzt+1
	2536	rad_lw_cs_hr_av(k,j,i) = rad_lw_cs_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
	2537	ENDDO
	2538	ENDDO
	2539	ENDDO
	2540
	2541	CASE ( 'rad_lw_hr' )
	2542	DO i = nxlg, nxrg
	2543	DO j = nysg, nyng
	2544	DO k = nzb, nzt+1
	2545	rad_lw_hr_av(k,j,i) = rad_lw_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
	2546	ENDDO
	2547	ENDDO
	2548	ENDDO
	2549
	2550	CASE ( 'rad_sw_in' )
	2551	DO i = nxlg, nxrg
	2552	DO j = nysg, nyng
	2553	DO k = nzb, nzt+1
	2554	rad_sw_in_av(k,j,i) = rad_sw_in_av(k,j,i) / REAL( average_count_3d, KIND=wp )
	2555	ENDDO
	2556	ENDDO
	2557	ENDDO
	2558
	2559	CASE ( 'rad_sw_out' )
	2560	DO i = nxlg, nxrg
	2561	DO j = nysg, nyng
	2562	DO k = nzb, nzt+1
	2563	rad_sw_out_av(k,j,i) = rad_sw_out_av(k,j,i) / REAL( average_count_3d, KIND=wp )
	2564	ENDDO
	2565	ENDDO
	2566	ENDDO
	2567
	2568	CASE ( 'rad_sw_cs_hr' )
	2569	DO i = nxlg, nxrg
	2570	DO j = nysg, nyng
	2571	DO k = nzb, nzt+1
	2572	rad_sw_cs_hr_av(k,j,i) = rad_sw_cs_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
	2573	ENDDO
	2574	ENDDO
	2575	ENDDO
	2576
	2577	CASE ( 'rad_sw_hr' )
	2578	DO i = nxlg, nxrg
	2579	DO j = nysg, nyng
	2580	DO k = nzb, nzt+1
	2581	rad_sw_hr_av(k,j,i) = rad_sw_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
	2582	ENDDO
	2583	ENDDO
	2584	ENDDO
	2585
	2586	END SELECT
	2587
	2588	ENDIF
	2589
	2590	END SUBROUTINE radiation_3d_data_averaging
	2591
	2592
	2593	!------------------------------------------------------------------------------!
	2594	!
	2595	! Description:
	2596	! ------------
	2597	!> Subroutine defining appropriate grid for netcdf variables.
	2598	!> It is called out from subroutine netcdf.
	2599	!------------------------------------------------------------------------------!
	2600	SUBROUTINE radiation_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
	2601
	2602	IMPLICIT NONE
	2603
	2604	CHARACTER (LEN=*), INTENT(IN) :: var !<
	2605	LOGICAL, INTENT(OUT) :: found !<
	2606	CHARACTER (LEN=*), INTENT(OUT) :: grid_x !<
	2607	CHARACTER (LEN=*), INTENT(OUT) :: grid_y !<
	2608	CHARACTER (LEN=*), INTENT(OUT) :: grid_z !<
	2609
	2610	found = .TRUE.
	2611
	2612
	2613	!
	2614	!-- Check for the grid
	2615	SELECT CASE ( TRIM( var ) )
	2616
	2617	CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_sw_cs_hr', 'rad_sw_hr', &
	2618	'rad_lw_cs_hr_xy', 'rad_lw_hr_xy', 'rad_sw_cs_hr_xy', &
	2619	'rad_sw_hr_xy', 'rad_lw_cs_hr_xz', 'rad_lw_hr_xz', &
	2620	'rad_sw_cs_hr_xz', 'rad_sw_hr_xz', 'rad_lw_cs_hr_yz', &
	2621	'rad_lw_hr_yz', 'rad_sw_cs_hr_yz', 'rad_sw_hr_yz' )
	2622	grid_x = 'x'
	2623	grid_y = 'y'
	2624	grid_z = 'zu'
	2625
	2626	CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_sw_in', 'rad_sw_out', &
	2627	'rad_lw_in_xy', 'rad_lw_out_xy', 'rad_sw_in_xy','rad_sw_out_xy', &
	2628	'rad_lw_in_xz', 'rad_lw_out_xz', 'rad_sw_in_xz','rad_sw_out_xz', &
	2629	'rad_lw_in_yz', 'rad_lw_out_yz', 'rad_sw_in_yz','rad_sw_out_yz' )
	2630	grid_x = 'x'
	2631	grid_y = 'y'
	2632	grid_z = 'zw'
	2633
	2634
	2635	CASE DEFAULT
	2636	found = .FALSE.
	2637	grid_x = 'none'
	2638	grid_y = 'none'
	2639	grid_z = 'none'
	2640
	2641	END SELECT
	2642
	2643	END SUBROUTINE radiation_define_netcdf_grid
	2644
	2645	!------------------------------------------------------------------------------!
	2646	!
	2647	! Description:
	2648	! ------------
	2649	!> Subroutine defining 3D output variables
	2650	!------------------------------------------------------------------------------!
	2651	SUBROUTINE radiation_data_output_2d( av, variable, found, grid, mode, &
	2652	local_pf, two_d )
	2653
	2654	USE indices
	2655
	2656	USE kinds
	2657
	2658
	2659	IMPLICIT NONE
	2660
	2661	CHARACTER (LEN=*) :: grid !<
	2662	CHARACTER (LEN=*) :: mode !<
	2663	CHARACTER (LEN=*) :: variable !<
	2664
	2665	INTEGER(iwp) :: av !<
	2666	INTEGER(iwp) :: i !<
	2667	INTEGER(iwp) :: j !<
	2668	INTEGER(iwp) :: k !<
	2669
	2670	LOGICAL :: found !<
	2671	LOGICAL :: two_d !< flag parameter that indicates 2D variables (horizontal cross sections)
	2672
	2673	REAL(wp), DIMENSION(nxlg:nxrg,nysg:nyng,nzb:nzt+1) :: local_pf !<
	2674
	2675	found = .TRUE.
	2676
	2677	SELECT CASE ( TRIM( variable ) )
	2678
	2679	CASE ( 'rad_net*_xy' ) ! 2d-array
	2680	IF ( av == 0 ) THEN
	2681	DO i = nxlg, nxrg
	2682	DO j = nysg, nyng
	2683	local_pf(i,j,nzb+1) = rad_net(j,i)
	2684	ENDDO
	2685	ENDDO
	2686	ELSE
	2687	DO i = nxlg, nxrg
	2688	DO j = nysg, nyng
	2689	local_pf(i,j,nzb+1) = rad_net_av(j,i)
	2690	ENDDO
	2691	ENDDO
	2692	ENDIF
	2693	two_d = .TRUE.
	2694	grid = 'zu1'
	2695
	2696
	2697	CASE ( 'rad_lw_in_xy', 'rad_lw_in_xz', 'rad_lw_in_yz' )
	2698	IF ( av == 0 ) THEN
	2699	DO i = nxlg, nxrg
	2700	DO j = nysg, nyng
	2701	DO k = nzb, nzt+1
	2702	local_pf(i,j,k) = rad_lw_in(k,j,i)
	2703	ENDDO
	2704	ENDDO
	2705	ENDDO
	2706	ELSE
	2707	DO i = nxlg, nxrg
	2708	DO j = nysg, nyng
	2709	DO k = nzb, nzt+1
	2710	local_pf(i,j,k) = rad_lw_in_av(k,j,i)
	2711	ENDDO
	2712	ENDDO
	2713	ENDDO
	2714	ENDIF
	2715	IF ( mode == 'xy' ) grid = 'zu'
	2716
	2717	CASE ( 'rad_lw_out_xy', 'rad_lw_out_xz', 'rad_lw_out_yz' )
	2718	IF ( av == 0 ) THEN
	2719	DO i = nxlg, nxrg
	2720	DO j = nysg, nyng
	2721	DO k = nzb, nzt+1
	2722	local_pf(i,j,k) = rad_lw_out(k,j,i)
	2723	ENDDO
	2724	ENDDO
	2725	ENDDO
	2726	ELSE
	2727	DO i = nxlg, nxrg
	2728	DO j = nysg, nyng
	2729	DO k = nzb, nzt+1
	2730	local_pf(i,j,k) = rad_lw_out_av(k,j,i)
	2731	ENDDO
	2732	ENDDO
	2733	ENDDO
	2734	ENDIF
	2735	IF ( mode == 'xy' ) grid = 'zu'
	2736
	2737	CASE ( 'rad_lw_cs_hr_xy', 'rad_lw_cs_hr_xz', 'rad_lw_cs_hr_yz' )
	2738	IF ( av == 0 ) THEN
	2739	DO i = nxlg, nxrg
	2740	DO j = nysg, nyng
	2741	DO k = nzb, nzt+1
	2742	local_pf(i,j,k) = rad_lw_cs_hr(k,j,i)
	2743	ENDDO
	2744	ENDDO
	2745	ENDDO
	2746	ELSE
	2747	DO i = nxlg, nxrg
	2748	DO j = nysg, nyng
	2749	DO k = nzb, nzt+1
	2750	local_pf(i,j,k) = rad_lw_cs_hr_av(k,j,i)
	2751	ENDDO
	2752	ENDDO
	2753	ENDDO
	2754	ENDIF
	2755	IF ( mode == 'xy' ) grid = 'zw'
	2756
	2757	CASE ( 'rad_lw_hr_xy', 'rad_lw_hr_xz', 'rad_lw_hr_yz' )
	2758	IF ( av == 0 ) THEN
	2759	DO i = nxlg, nxrg
	2760	DO j = nysg, nyng
	2761	DO k = nzb, nzt+1
	2762	local_pf(i,j,k) = rad_lw_hr(k,j,i)
	2763	ENDDO
	2764	ENDDO
	2765	ENDDO
	2766	ELSE
	2767	DO i = nxlg, nxrg
	2768	DO j = nysg, nyng
	2769	DO k = nzb, nzt+1
	2770	local_pf(i,j,k) = rad_lw_hr_av(k,j,i)
	2771	ENDDO
	2772	ENDDO
	2773	ENDDO
	2774	ENDIF
	2775	IF ( mode == 'xy' ) grid = 'zw'
	2776
	2777	CASE ( 'rad_sw_in_xy', 'rad_sw_in_xz', 'rad_sw_in_yz' )
	2778	IF ( av == 0 ) THEN
	2779	DO i = nxlg, nxrg
	2780	DO j = nysg, nyng
	2781	DO k = nzb, nzt+1
	2782	local_pf(i,j,k) = rad_sw_in(k,j,i)
	2783	ENDDO
	2784	ENDDO
	2785	ENDDO
	2786	ELSE
	2787	DO i = nxlg, nxrg
	2788	DO j = nysg, nyng
	2789	DO k = nzb, nzt+1
	2790	local_pf(i,j,k) = rad_sw_in_av(k,j,i)
	2791	ENDDO
	2792	ENDDO
	2793	ENDDO
	2794	ENDIF
	2795	IF ( mode == 'xy' ) grid = 'zu'
	2796
	2797	CASE ( 'rad_sw_out_xy', 'rad_sw_out_xz', 'rad_sw_out_yz' )
	2798	IF ( av == 0 ) THEN
	2799	DO i = nxlg, nxrg
	2800	DO j = nysg, nyng
	2801	DO k = nzb, nzt+1
	2802	local_pf(i,j,k) = rad_sw_out(k,j,i)
	2803	ENDDO
	2804	ENDDO
	2805	ENDDO
	2806	ELSE
	2807	DO i = nxlg, nxrg
	2808	DO j = nysg, nyng
	2809	DO k = nzb, nzt+1
	2810	local_pf(i,j,k) = rad_sw_out_av(k,j,i)
	2811	ENDDO
	2812	ENDDO
	2813	ENDDO
	2814	ENDIF
	2815	IF ( mode == 'xy' ) grid = 'zu'
	2816
	2817	CASE ( 'rad_sw_cs_hr_xy', 'rad_sw_cs_hr_xz', 'rad_sw_cs_hr_yz' )
	2818	IF ( av == 0 ) THEN
	2819	DO i = nxlg, nxrg
	2820	DO j = nysg, nyng
	2821	DO k = nzb, nzt+1
	2822	local_pf(i,j,k) = rad_sw_cs_hr(k,j,i)
	2823	ENDDO
	2824	ENDDO
	2825	ENDDO
	2826	ELSE
	2827	DO i = nxlg, nxrg
	2828	DO j = nysg, nyng
	2829	DO k = nzb, nzt+1
	2830	local_pf(i,j,k) = rad_sw_cs_hr_av(k,j,i)
	2831	ENDDO
	2832	ENDDO
	2833	ENDDO
	2834	ENDIF
	2835	IF ( mode == 'xy' ) grid = 'zw'
	2836
	2837	CASE ( 'rad_sw_hr_xy', 'rad_sw_hr_xz', 'rad_sw_hr_yz' )
	2838	IF ( av == 0 ) THEN
	2839	DO i = nxlg, nxrg
	2840	DO j = nysg, nyng
	2841	DO k = nzb, nzt+1
	2842	local_pf(i,j,k) = rad_sw_hr(k,j,i)
	2843	ENDDO
	2844	ENDDO
	2845	ENDDO
	2846	ELSE
	2847	DO i = nxlg, nxrg
	2848	DO j = nysg, nyng
	2849	DO k = nzb, nzt+1
	2850	local_pf(i,j,k) = rad_sw_hr_av(k,j,i)
	2851	ENDDO
	2852	ENDDO
	2853	ENDDO
	2854	ENDIF
	2855	IF ( mode == 'xy' ) grid = 'zw'
	2856
	2857	CASE DEFAULT
	2858	found = .FALSE.
	2859	grid = 'none'
	2860
	2861	END SELECT
	2862
	2863	END SUBROUTINE radiation_data_output_2d
	2864
	2865
	2866	!------------------------------------------------------------------------------!
	2867	!
	2868	! Description:
	2869	! ------------
	2870	!> Subroutine defining 3D output variables
	2871	!------------------------------------------------------------------------------!
	2872	SUBROUTINE radiation_data_output_3d( av, variable, found, local_pf )
	2873
	2874
	2875	USE indices
	2876
	2877	USE kinds
	2878
	2879
	2880	IMPLICIT NONE
	2881
	2882	CHARACTER (LEN=*) :: variable !<
	2883
	2884	INTEGER(iwp) :: av !<
	2885	INTEGER(iwp) :: i !<
	2886	INTEGER(iwp) :: j !<
	2887	INTEGER(iwp) :: k !<
	2888
	2889	LOGICAL :: found !<
	2890
	2891	REAL(sp), DIMENSION(nxlg:nxrg,nysg:nyng,nzb:nzt+1) :: local_pf !<
	2892
	2893
	2894	found = .TRUE.
	2895
	2896
	2897	SELECT CASE ( TRIM( variable ) )
	2898
	2899	CASE ( 'rad_sw_in' )
	2900	IF ( av == 0 ) THEN
	2901	DO i = nxlg, nxrg
	2902	DO j = nysg, nyng
	2903	DO k = nzb, nzt+1
	2904	local_pf(i,j,k) = rad_sw_in(k,j,i)
	2905	ENDDO
	2906	ENDDO
	2907	ENDDO
	2908	ELSE
	2909	DO i = nxlg, nxrg
	2910	DO j = nysg, nyng
	2911	DO k = nzb, nzt+1
	2912	local_pf(i,j,k) = rad_sw_in_av(k,j,i)
	2913	ENDDO
	2914	ENDDO
	2915	ENDDO
	2916	ENDIF
	2917
	2918	CASE ( 'rad_sw_out' )
	2919	IF ( av == 0 ) THEN
	2920	DO i = nxlg, nxrg
	2921	DO j = nysg, nyng
	2922	DO k = nzb, nzt+1
	2923	local_pf(i,j,k) = rad_sw_out(k,j,i)
	2924	ENDDO
	2925	ENDDO
	2926	ENDDO
	2927	ELSE
	2928	DO i = nxlg, nxrg
	2929	DO j = nysg, nyng
	2930	DO k = nzb, nzt+1
	2931	local_pf(i,j,k) = rad_sw_out_av(k,j,i)
	2932	ENDDO
	2933	ENDDO
	2934	ENDDO
	2935	ENDIF
	2936
	2937	CASE ( 'rad_sw_cs_hr' )
	2938	IF ( av == 0 ) THEN
	2939	DO i = nxlg, nxrg
	2940	DO j = nysg, nyng
	2941	DO k = nzb, nzt+1
	2942	local_pf(i,j,k) = rad_sw_cs_hr(k,j,i)
	2943	ENDDO
	2944	ENDDO
	2945	ENDDO
	2946	ELSE
	2947	DO i = nxlg, nxrg
	2948	DO j = nysg, nyng
	2949	DO k = nzb, nzt+1
	2950	local_pf(i,j,k) = rad_sw_cs_hr_av(k,j,i)
	2951	ENDDO
	2952	ENDDO
	2953	ENDDO
	2954	ENDIF
	2955
	2956	CASE ( 'rad_sw_hr' )
	2957	IF ( av == 0 ) THEN
	2958	DO i = nxlg, nxrg
	2959	DO j = nysg, nyng
	2960	DO k = nzb, nzt+1
	2961	local_pf(i,j,k) = rad_sw_hr(k,j,i)
	2962	ENDDO
	2963	ENDDO
	2964	ENDDO
	2965	ELSE
	2966	DO i = nxlg, nxrg
	2967	DO j = nysg, nyng
	2968	DO k = nzb, nzt+1
	2969	local_pf(i,j,k) = rad_sw_hr_av(k,j,i)
	2970	ENDDO
	2971	ENDDO
	2972	ENDDO
	2973	ENDIF
	2974
	2975	CASE ( 'rad_lw_in' )
	2976	IF ( av == 0 ) THEN
	2977	DO i = nxlg, nxrg
	2978	DO j = nysg, nyng
	2979	DO k = nzb, nzt+1
	2980	local_pf(i,j,k) = rad_lw_in(k,j,i)
	2981	ENDDO
	2982	ENDDO
	2983	ENDDO
	2984	ELSE
	2985	DO i = nxlg, nxrg
	2986	DO j = nysg, nyng
	2987	DO k = nzb, nzt+1
	2988	local_pf(i,j,k) = rad_lw_in_av(k,j,i)
	2989	ENDDO
	2990	ENDDO
	2991	ENDDO
	2992	ENDIF
	2993
	2994	CASE ( 'rad_lw_out' )
	2995	IF ( av == 0 ) THEN
	2996	DO i = nxlg, nxrg
	2997	DO j = nysg, nyng
	2998	DO k = nzb, nzt+1
	2999	local_pf(i,j,k) = rad_lw_out(k,j,i)
	3000	ENDDO
	3001	ENDDO
	3002	ENDDO
	3003	ELSE
	3004	DO i = nxlg, nxrg
	3005	DO j = nysg, nyng
	3006	DO k = nzb, nzt+1
	3007	local_pf(i,j,k) = rad_lw_out_av(k,j,i)
	3008	ENDDO
	3009	ENDDO
	3010	ENDDO
	3011	ENDIF
	3012
	3013	CASE ( 'rad_lw_cs_hr' )
	3014	IF ( av == 0 ) THEN
	3015	DO i = nxlg, nxrg
	3016	DO j = nysg, nyng
	3017	DO k = nzb, nzt+1
	3018	local_pf(i,j,k) = rad_lw_cs_hr(k,j,i)
	3019	ENDDO
	3020	ENDDO
	3021	ENDDO
	3022	ELSE
	3023	DO i = nxlg, nxrg
	3024	DO j = nysg, nyng
	3025	DO k = nzb, nzt+1
	3026	local_pf(i,j,k) = rad_lw_cs_hr_av(k,j,i)
	3027	ENDDO
	3028	ENDDO
	3029	ENDDO
	3030	ENDIF
	3031
	3032	CASE ( 'rad_lw_hr' )
	3033	IF ( av == 0 ) THEN
	3034	DO i = nxlg, nxrg
	3035	DO j = nysg, nyng
	3036	DO k = nzb, nzt+1
	3037	local_pf(i,j,k) = rad_lw_hr(k,j,i)
	3038	ENDDO
	3039	ENDDO
	3040	ENDDO
	3041	ELSE
	3042	DO i = nxlg, nxrg
	3043	DO j = nysg, nyng
	3044	DO k = nzb, nzt+1
	3045	local_pf(i,j,k) = rad_lw_hr_av(k,j,i)
	3046	ENDDO
	3047	ENDDO
	3048	ENDDO
	3049	ENDIF
	3050
	3051	CASE DEFAULT
	3052	found = .FALSE.
	3053
	3054	END SELECT
	3055
	3056
	3057	END SUBROUTINE radiation_data_output_3d
	3058
	3059	!------------------------------------------------------------------------------!
	3060	!
	3061	! Description:
	3062	! ------------
	3063	!> Subroutine defining masked data output
	3064	!------------------------------------------------------------------------------!
	3065	SUBROUTINE radiation_data_output_mask( av, variable, found, local_pf )
	3066
	3067	USE control_parameters
	3068
	3069	USE indices
	3070
	3071	USE kinds
	3072
	3073
	3074	IMPLICIT NONE
	3075
	3076	CHARACTER (LEN=*) :: variable !<
	3077
	3078	INTEGER(iwp) :: av !<
	3079	INTEGER(iwp) :: i !<
	3080	INTEGER(iwp) :: j !<
	3081	INTEGER(iwp) :: k !<
	3082
	3083	LOGICAL :: found !<
	3084
	3085	REAL(wp), &
	3086	DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) :: &
	3087	local_pf !<
	3088
	3089
	3090	found = .TRUE.
	3091
	3092	SELECT CASE ( TRIM( variable ) )
	3093
	3094
	3095	CASE ( 'rad_lw_in' )
	3096	IF ( av == 0 ) THEN
	3097	DO i = 1, mask_size_l(mid,1)
	3098	DO j = 1, mask_size_l(mid,2)
	3099	DO k = 1, mask_size_l(mid,3)
	3100	local_pf(i,j,k) = rad_lw_in(mask_k(mid,k), &
	3101	mask_j(mid,j),mask_i(mid,i))
	3102	ENDDO
	3103	ENDDO
	3104	ENDDO
	3105	ELSE
	3106	DO i = 1, mask_size_l(mid,1)
	3107	DO j = 1, mask_size_l(mid,2)
	3108	DO k = 1, mask_size_l(mid,3)
	3109	local_pf(i,j,k) = rad_lw_in_av(mask_k(mid,k), &
	3110	mask_j(mid,j),mask_i(mid,i))
	3111	ENDDO
	3112	ENDDO
	3113	ENDDO
	3114	ENDIF
	3115
	3116	CASE ( 'rad_lw_out' )
	3117	IF ( av == 0 ) THEN
	3118	DO i = 1, mask_size_l(mid,1)
	3119	DO j = 1, mask_size_l(mid,2)
	3120	DO k = 1, mask_size_l(mid,3)
	3121	local_pf(i,j,k) = rad_lw_out(mask_k(mid,k), &
	3122	mask_j(mid,j),mask_i(mid,i))
	3123	ENDDO
	3124	ENDDO
	3125	ENDDO
	3126	ELSE
	3127	DO i = 1, mask_size_l(mid,1)
	3128	DO j = 1, mask_size_l(mid,2)
	3129	DO k = 1, mask_size_l(mid,3)
	3130	local_pf(i,j,k) = rad_lw_out_av(mask_k(mid,k), &
	3131	mask_j(mid,j),mask_i(mid,i))
	3132	ENDDO
	3133	ENDDO
	3134	ENDDO
	3135	ENDIF
	3136
	3137	CASE ( 'rad_lw_cs_hr' )
	3138	IF ( av == 0 ) THEN
	3139	DO i = 1, mask_size_l(mid,1)
	3140	DO j = 1, mask_size_l(mid,2)
	3141	DO k = 1, mask_size_l(mid,3)
	3142	local_pf(i,j,k) = rad_lw_cs_hr(mask_k(mid,k), &
	3143	mask_j(mid,j),mask_i(mid,i))
	3144	ENDDO
	3145	ENDDO
	3146	ENDDO
	3147	ELSE
	3148	DO i = 1, mask_size_l(mid,1)
	3149	DO j = 1, mask_size_l(mid,2)
	3150	DO k = 1, mask_size_l(mid,3)
	3151	local_pf(i,j,k) = rad_lw_cs_hr_av(mask_k(mid,k), &
	3152	mask_j(mid,j),mask_i(mid,i))
	3153	ENDDO
	3154	ENDDO
	3155	ENDDO
	3156	ENDIF
	3157
	3158	CASE ( 'rad_lw_hr' )
	3159	IF ( av == 0 ) THEN
	3160	DO i = 1, mask_size_l(mid,1)
	3161	DO j = 1, mask_size_l(mid,2)
	3162	DO k = 1, mask_size_l(mid,3)
	3163	local_pf(i,j,k) = rad_lw_hr(mask_k(mid,k), &
	3164	mask_j(mid,j),mask_i(mid,i))
	3165	ENDDO
	3166	ENDDO
	3167	ENDDO
	3168	ELSE
	3169	DO i = 1, mask_size_l(mid,1)
	3170	DO j = 1, mask_size_l(mid,2)
	3171	DO k = 1, mask_size_l(mid,3)
	3172	local_pf(i,j,k) = rad_lw_hr_av(mask_k(mid,k), &
	3173	mask_j(mid,j),mask_i(mid,i))
	3174	ENDDO
	3175	ENDDO
	3176	ENDDO
	3177	ENDIF
	3178
	3179	CASE ( 'rad_sw_in' )
	3180	IF ( av == 0 ) THEN
	3181	DO i = 1, mask_size_l(mid,1)
	3182	DO j = 1, mask_size_l(mid,2)
	3183	DO k = 1, mask_size_l(mid,3)
	3184	local_pf(i,j,k) = rad_sw_in(mask_k(mid,k), &
	3185	mask_j(mid,j),mask_i(mid,i))
	3186	ENDDO
	3187	ENDDO
	3188	ENDDO
	3189	ELSE
	3190	DO i = 1, mask_size_l(mid,1)
	3191	DO j = 1, mask_size_l(mid,2)
	3192	DO k = 1, mask_size_l(mid,3)
	3193	local_pf(i,j,k) = rad_sw_in_av(mask_k(mid,k), &
	3194	mask_j(mid,j),mask_i(mid,i))
	3195	ENDDO
	3196	ENDDO
	3197	ENDDO
	3198	ENDIF
	3199
	3200	CASE ( 'rad_sw_out' )
	3201	IF ( av == 0 ) THEN
	3202	DO i = 1, mask_size_l(mid,1)
	3203	DO j = 1, mask_size_l(mid,2)
	3204	DO k = 1, mask_size_l(mid,3)
	3205	local_pf(i,j,k) = rad_sw_out(mask_k(mid,k), &
	3206	mask_j(mid,j),mask_i(mid,i))
	3207	ENDDO
	3208	ENDDO
	3209	ENDDO
	3210	ELSE
	3211	DO i = 1, mask_size_l(mid,1)
	3212	DO j = 1, mask_size_l(mid,2)
	3213	DO k = 1, mask_size_l(mid,3)
	3214	local_pf(i,j,k) = rad_sw_out_av(mask_k(mid,k), &
	3215	mask_j(mid,j),mask_i(mid,i))
	3216	ENDDO
	3217	ENDDO
	3218	ENDDO
	3219	ENDIF
	3220
	3221	CASE ( 'rad_sw_cs_hr' )
	3222	IF ( av == 0 ) THEN
	3223	DO i = 1, mask_size_l(mid,1)
	3224	DO j = 1, mask_size_l(mid,2)
	3225	DO k = 1, mask_size_l(mid,3)
	3226	local_pf(i,j,k) = rad_sw_cs_hr(mask_k(mid,k), &
	3227	mask_j(mid,j),mask_i(mid,i))
	3228	ENDDO
	3229	ENDDO
	3230	ENDDO
	3231	ELSE
	3232	DO i = 1, mask_size_l(mid,1)
	3233	DO j = 1, mask_size_l(mid,2)
	3234	DO k = 1, mask_size_l(mid,3)
	3235	local_pf(i,j,k) = rad_sw_cs_hr_av(mask_k(mid,k), &
	3236	mask_j(mid,j),mask_i(mid,i))
	3237	ENDDO
	3238	ENDDO
	3239	ENDDO
	3240	ENDIF
	3241
	3242	CASE ( 'rad_sw_hr' )
	3243	IF ( av == 0 ) THEN
	3244	DO i = 1, mask_size_l(mid,1)
	3245	DO j = 1, mask_size_l(mid,2)
	3246	DO k = 1, mask_size_l(mid,3)
	3247	local_pf(i,j,k) = rad_sw_hr(mask_k(mid,k), &
	3248	mask_j(mid,j),mask_i(mid,i))
	3249	ENDDO
	3250	ENDDO
	3251	ENDDO
	3252	ELSE
	3253	DO i = 1, mask_size_l(mid,1)
	3254	DO j = 1, mask_size_l(mid,2)
	3255	DO k = 1, mask_size_l(mid,3)
	3256	local_pf(i,j,k) = rad_sw_hr_av(mask_k(mid,k), &
	3257	mask_j(mid,j),mask_i(mid,i))
	3258	ENDDO
	3259	ENDDO
	3260	ENDDO
	3261	ENDIF
	3262
	3263	CASE DEFAULT
	3264	found = .FALSE.
	3265
	3266	END SELECT
	3267
	3268
	3269	END SUBROUTINE radiation_data_output_mask
	3270
	3271
	3272	!------------------------------------------------------------------------------!
	3273	!
	3274	! Description:
	3275	! ------------
	3276	!> Subroutine defines masked output variables
	3277	!------------------------------------------------------------------------------!
	3278	SUBROUTINE radiation_last_actions
	3279
	3280
	3281	USE control_parameters
	3282
	3283	USE kinds
	3284
	3285	IMPLICIT NONE
	3286
	3287	IF ( write_binary(1:4) == 'true' ) THEN
	3288	IF ( ALLOCATED( rad_net ) ) THEN
	3289	WRITE ( 14 ) 'rad_net '; WRITE ( 14 ) rad_net
	3290	ENDIF
	3291	IF ( ALLOCATED( rad_net_av ) ) THEN
	3292	WRITE ( 14 ) 'rad_net_av '; WRITE ( 14 ) rad_net_av
	3293	ENDIF
	3294	IF ( ALLOCATED( rad_lw_in ) ) THEN
	3295	WRITE ( 14 ) 'rad_lw_in '; WRITE ( 14 ) rad_lw_in
	3296	ENDIF
	3297	IF ( ALLOCATED( rad_lw_in_av ) ) THEN
	3298	WRITE ( 14 ) 'rad_lw_in_av '; WRITE ( 14 ) rad_lw_in_av
	3299	ENDIF
	3300	IF ( ALLOCATED( rad_lw_out ) ) THEN
	3301	WRITE ( 14 ) 'rad_lw_out '; WRITE ( 14 ) rad_lw_out
	3302	ENDIF
	3303	IF ( ALLOCATED( rad_lw_out_av ) ) THEN
	3304	WRITE ( 14 ) 'rad_lw_out_av '; WRITE ( 14 ) rad_lw_out_av
	3305	ENDIF
	3306	IF ( ALLOCATED( rad_lw_out_change_0 ) ) THEN
	3307	WRITE ( 14 ) 'rad_lw_out_change_0 '
	3308	WRITE ( 14 ) rad_lw_out_change_0
	3309	ENDIF
	3310	IF ( ALLOCATED( rad_lw_cs_hr ) ) THEN
	3311	WRITE ( 14 ) 'rad_lw_cs_hr '; WRITE ( 14 ) rad_lw_cs_hr
	3312	ENDIF
	3313	IF ( ALLOCATED( rad_lw_cs_hr_av ) ) THEN
	3314	WRITE ( 14 ) 'rad_lw_cs_hr_av '; WRITE ( 14 ) rad_lw_cs_hr_av
	3315	ENDIF
	3316	IF ( ALLOCATED( rad_lw_hr ) ) THEN
	3317	WRITE ( 14 ) 'rad_lw_hr '; WRITE ( 14 ) rad_lw_hr
	3318	ENDIF
	3319	IF ( ALLOCATED( rad_lw_hr_av ) ) THEN
	3320	WRITE ( 14 ) 'rad_lw_hr_av '; WRITE ( 14 ) rad_lw_hr_av
	3321	ENDIF
	3322	IF ( ALLOCATED( rad_sw_in ) ) THEN
	3323	WRITE ( 14 ) 'rad_sw_in '; WRITE ( 14 ) rad_sw_in
	3324	ENDIF
	3325	IF ( ALLOCATED( rad_sw_in_av ) ) THEN
	3326	WRITE ( 14 ) 'rad_sw_in_av '; WRITE ( 14 ) rad_sw_in_av
	3327	ENDIF
	3328	IF ( ALLOCATED( rad_sw_out ) ) THEN
	3329	WRITE ( 14 ) 'rad_sw_out '; WRITE ( 14 ) rad_sw_out
	3330	ENDIF
	3331	IF ( ALLOCATED( rad_sw_out_av ) ) THEN
	3332	WRITE ( 14 ) 'rad_sw_out_av '; WRITE ( 14 ) rad_sw_out_av
	3333	ENDIF
	3334	IF ( ALLOCATED( rad_sw_cs_hr ) ) THEN
	3335	WRITE ( 14 ) 'rad_sw_cs_hr '; WRITE ( 14 ) rad_sw_cs_hr
	3336	ENDIF
	3337	IF ( ALLOCATED( rad_sw_cs_hr_av ) ) THEN
	3338	WRITE ( 14 ) 'rad_sw_cs_hr_av '; WRITE ( 14 ) rad_sw_cs_hr_av
	3339	ENDIF
	3340	IF ( ALLOCATED( rad_sw_hr ) ) THEN
	3341	WRITE ( 14 ) 'rad_sw_hr '; WRITE ( 14 ) rad_sw_hr
	3342	ENDIF
	3343	IF ( ALLOCATED( rad_sw_hr_av ) ) THEN
	3344	WRITE ( 14 ) 'rad_sw_hr_av '; WRITE ( 14 ) rad_sw_hr_av
	3345	ENDIF
	3346
	3347	WRITE ( 14 ) '* end rad * '
	3348
	3349	ENDIF
	3350
	3351	END SUBROUTINE radiation_last_actions
	3352
	3353
	3354	SUBROUTINE radiation_read_restart_data( i, nxlfa, nxl_on_file, nxrfa, nxr_on_file, &
	3355	nynfa, nyn_on_file, nysfa, nys_on_file, &
	3356	offset_xa, offset_ya, overlap_count, &
	3357	tmp_2d, tmp_3d )
	3358
	3359
	3360	USE control_parameters
	3361
	3362	USE indices
	3363
	3364	USE kinds
	3365
	3366	USE pegrid
	3367
	3368	IMPLICIT NONE
	3369
	3370	CHARACTER (LEN=20) :: field_char !<
	3371
	3372	INTEGER(iwp) :: i !<
	3373	INTEGER(iwp) :: k !<
	3374	INTEGER(iwp) :: nxlc !<
	3375	INTEGER(iwp) :: nxlf !<
	3376	INTEGER(iwp) :: nxl_on_file !<
	3377	INTEGER(iwp) :: nxrc !<
	3378	INTEGER(iwp) :: nxrf !<
	3379	INTEGER(iwp) :: nxr_on_file !<
	3380	INTEGER(iwp) :: nync !<
	3381	INTEGER(iwp) :: nynf !<
	3382	INTEGER(iwp) :: nyn_on_file !<
	3383	INTEGER(iwp) :: nysc !<
	3384	INTEGER(iwp) :: nysf !<
	3385	INTEGER(iwp) :: nys_on_file !<
	3386	INTEGER(iwp) :: overlap_count !<
	3387
	3388	INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: nxlfa !<
	3389	INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: nxrfa !<
	3390	INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: nynfa !<
	3391	INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: nysfa !<
	3392	INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: offset_xa !<
	3393	INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: offset_ya !<
	3394
	3395	REAL(wp), &
	3396	DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) ::&
	3397	tmp_2d !<
	3398
	3399	REAL(wp), &
	3400	DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) ::&
	3401	tmp_3d !<
	3402
	3403
	3404
	3405	IF ( initializing_actions == 'read_restart_data' ) THEN
	3406	READ ( 13 ) field_char
	3407
	3408	DO WHILE ( TRIM( field_char ) /= '* end rad *' )
	3409
	3410	DO k = 1, overlap_count
	3411
	3412	nxlf = nxlfa(i,k)
	3413	nxlc = nxlfa(i,k) + offset_xa(i,k)
	3414	nxrf = nxrfa(i,k)
	3415	nxrc = nxrfa(i,k) + offset_xa(i,k)
	3416	nysf = nysfa(i,k)
	3417	nysc = nysfa(i,k) + offset_ya(i,k)
	3418	nynf = nynfa(i,k)
	3419	nync = nynfa(i,k) + offset_ya(i,k)
	3420
	3421
	3422	SELECT CASE ( TRIM( field_char ) )
	3423
	3424	CASE ( 'rad_net' )
	3425	IF ( .NOT. ALLOCATED( rad_net ) ) THEN
	3426	ALLOCATE( rad_net(nysg:nyng,nxlg:nxrg) )
	3427	ENDIF
	3428	IF ( k == 1 ) READ ( 13 ) tmp_2d
	3429	rad_net(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
	3430	tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
	3431
	3432	CASE ( 'rad_net_av' )
	3433	IF ( .NOT. ALLOCATED( rad_net_av ) ) THEN
	3434	ALLOCATE( rad_net_av(nysg:nyng,nxlg:nxrg) )
	3435	ENDIF
	3436	IF ( k == 1 ) READ ( 13 ) tmp_2d
	3437	rad_net_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
	3438	tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
	3439	CASE ( 'rad_lw_in' )
	3440	IF ( .NOT. ALLOCATED( rad_lw_in ) ) THEN
	3441	ALLOCATE( rad_lw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	3442	ENDIF
	3443	IF ( k == 1 ) READ ( 13 ) tmp_3d
	3444	rad_lw_in(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
	3445	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
	3446
	3447	CASE ( 'rad_lw_in_av' )
	3448	IF ( .NOT. ALLOCATED( rad_lw_in_av ) ) THEN
	3449	ALLOCATE( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	3450	ENDIF
	3451	IF ( k == 1 ) READ ( 13 ) tmp_3d
	3452	rad_lw_in_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
	3453	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
	3454
	3455	CASE ( 'rad_lw_out' )
	3456	IF ( .NOT. ALLOCATED( rad_lw_out ) ) THEN
	3457	ALLOCATE( rad_lw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	3458	ENDIF
	3459	IF ( k == 1 ) READ ( 13 ) tmp_3d
	3460	rad_lw_out(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
	3461	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
	3462
	3463	CASE ( 'rad_lw_out_av' )
	3464	IF ( .NOT. ALLOCATED( rad_lw_out_av ) ) THEN
	3465	ALLOCATE( rad_lw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	3466	ENDIF
	3467	IF ( k == 1 ) READ ( 13 ) tmp_3d
	3468	rad_lw_out_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
	3469	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
	3470
	3471	CASE ( 'rad_lw_out_change_0' )
	3472	IF ( .NOT. ALLOCATED( rad_lw_out_change_0 ) ) THEN
	3473	ALLOCATE( rad_lw_out_change_0(nysg:nyng,nxlg:nxrg) )
	3474	ENDIF
	3475	IF ( k == 1 ) READ ( 13 ) tmp_2d
	3476	rad_lw_out_change_0(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)&
	3477	= tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
	3478
	3479	CASE ( 'rad_lw_cs_hr' )
	3480	IF ( .NOT. ALLOCATED( rad_lw_cs_hr ) ) THEN
	3481	ALLOCATE( rad_lw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	3482	ENDIF
	3483	IF ( k == 1 ) READ ( 13 ) tmp_3d
	3484	rad_lw_cs_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
	3485	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
	3486
	3487	CASE ( 'rad_lw_cs_hr_av' )
	3488	IF ( .NOT. ALLOCATED( rad_lw_cs_hr_av ) ) THEN
	3489	ALLOCATE( rad_lw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	3490	ENDIF
	3491	IF ( k == 1 ) READ ( 13 ) tmp_3d
	3492	rad_lw_cs_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
	3493	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
	3494
	3495	CASE ( 'rad_lw_hr' )
	3496	IF ( .NOT. ALLOCATED( rad_lw_hr ) ) THEN
	3497	ALLOCATE( rad_lw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	3498	ENDIF
	3499	IF ( k == 1 ) READ ( 13 ) tmp_3d
	3500	rad_lw_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
	3501	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
	3502
	3503	CASE ( 'rad_lw_hr_av' )
	3504	IF ( .NOT. ALLOCATED( rad_lw_hr_av ) ) THEN
	3505	ALLOCATE( rad_lw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	3506	ENDIF
	3507	IF ( k == 1 ) READ ( 13 ) tmp_3d
	3508	rad_lw_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
	3509	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
	3510
	3511	CASE ( 'rad_sw_in' )
	3512	IF ( .NOT. ALLOCATED( rad_sw_in ) ) THEN
	3513	ALLOCATE( rad_sw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	3514	ENDIF
	3515	IF ( k == 1 ) READ ( 13 ) tmp_3d
	3516	rad_sw_in(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
	3517	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
	3518
	3519	CASE ( 'rad_sw_in_av' )
	3520	IF ( .NOT. ALLOCATED( rad_sw_in_av ) ) THEN
	3521	ALLOCATE( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	3522	ENDIF
	3523	IF ( k == 1 ) READ ( 13 ) tmp_3d
	3524	rad_sw_in_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
	3525	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
	3526
	3527	CASE ( 'rad_sw_out' )
	3528	IF ( .NOT. ALLOCATED( rad_sw_out ) ) THEN
	3529	ALLOCATE( rad_sw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	3530	ENDIF
	3531	IF ( k == 1 ) READ ( 13 ) tmp_3d
	3532	rad_sw_out(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
	3533	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
	3534
	3535	CASE ( 'rad_sw_out_av' )
	3536	IF ( .NOT. ALLOCATED( rad_sw_out_av ) ) THEN
	3537	ALLOCATE( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	3538	ENDIF
	3539	IF ( k == 1 ) READ ( 13 ) tmp_3d
	3540	rad_sw_out_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
	3541	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
	3542
	3543	CASE ( 'rad_sw_cs_hr' )
	3544	IF ( .NOT. ALLOCATED( rad_sw_cs_hr ) ) THEN
	3545	ALLOCATE( rad_sw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	3546	ENDIF
	3547	IF ( k == 1 ) READ ( 13 ) tmp_3d
	3548	rad_sw_cs_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
	3549	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
	3550
	3551	CASE ( 'rad_sw_cs_hr_av' )
	3552	IF ( .NOT. ALLOCATED( rad_sw_cs_hr_av ) ) THEN
	3553	ALLOCATE( rad_sw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	3554	ENDIF
	3555	IF ( k == 1 ) READ ( 13 ) tmp_3d
	3556	rad_sw_cs_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
	3557	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
	3558
	3559	CASE ( 'rad_sw_hr' )
	3560	IF ( .NOT. ALLOCATED( rad_sw_hr ) ) THEN
	3561	ALLOCATE( rad_sw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	3562	ENDIF
	3563	IF ( k == 1 ) READ ( 13 ) tmp_3d
	3564	rad_sw_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
	3565	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
	3566
	3567	CASE ( 'rad_sw_hr_av' )
	3568	IF ( .NOT. ALLOCATED( rad_sw_hr_av ) ) THEN
	3569	ALLOCATE( rad_sw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	3570	ENDIF
	3571	IF ( k == 1 ) READ ( 13 ) tmp_3d
	3572	rad_lw_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
	3573	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
	3574
	3575	CASE DEFAULT
	3576	WRITE( message_string, * ) 'unknown variable named "', &
	3577	TRIM( field_char ), '" found in', &
	3578	'&data from prior run on PE ', myid
	3579	CALL message( 'radiation_read_restart_data', 'PA0441', 1, 2, 0, 6, &
	3580	0 )
	3581
	3582	END SELECT
	3583
	3584	ENDDO
	3585
	3586	READ ( 13 ) field_char
	3587
	3588	ENDDO
	3589	ENDIF
	3590
	3591	END SUBROUTINE radiation_read_restart_data
	3592
	3593
[1496]	3594	END MODULE radiation_model_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |