Home

Context Navigation

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

Last change on this file since 2296 was 2296, checked in by maronga, 8 years ago
added new spinup mechanism for surface/radiation models
Property svn:keywords set to `Id`
File size: 144.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |