Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |