Home

Context Navigation

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

Last change on this file since 2065 was 2012, checked in by kanani, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 135.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |