Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |