Home

Context Navigation

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

Last change on this file since 2317 was 2317, checked in by suehring, 7 years ago
get topograpyh top index via function call
Property svn:keywords set to `Id`
File size: 145.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |