Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |