Home

Context Navigation

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

Last change on this file since 2224 was 2201, checked in by suehring, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 142.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |