Home

Context Navigation

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

Last change on this file since 2281 was 2270, checked in by maronga, 7 years ago
major revisions in land surface model
Property svn:keywords set to `Id`
File size: 144.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |