Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |