Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |