Home

Context Navigation

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

Last change on this file since 1989 was 1977, checked in by maronga, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 134.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |