Home

Context Navigation

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

Last change on this file since 1858 was 1857, checked in by maronga, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 89.2 KB

Rev	Line
[1826]	1	!> @file radiation_model_mod.f90
[1496]	2	!--------------------------------------------------------------------------------!
	3	! This file is part of PALM.
	4	!
	5	! PALM is free software: you can redistribute it and/or modify it under the terms
	6	! of the GNU General Public License as published by the Free Software Foundation,
	7	! either version 3 of the License, or (at your option) any later version.
	8	!
	9	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
	10	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
	11	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
	12	!
	13	! You should have received a copy of the GNU General Public License along with
	14	! PALM. If not, see <http://www.gnu.org/licenses/>.
	15	!
[1818]	16	! Copyright 1997-2016 Leibniz Universitaet Hannover
[1496]	17	!--------------------------------------------------------------------------------!
	18	!
	19	! Current revisions:
	20	! -----------------
[1854]	21	!
[1857]	22	!
[1496]	23	! Former revisions:
	24	! -----------------
	25	! $Id: radiation_model_mod.f90 1857 2016-04-13 12:56:38Z hoffmann $
	26	!
[1857]	27	! 1856 2016-04-13 12:56:17Z maronga
	28	! Bugfix: allocation of rad_lw_out for radiation_scheme = 'clear-sky'
	29	!
[1854]	30	! 1853 2016-04-11 09:00:35Z maronga
	31	! Added routine for radiation_scheme = constant.
	32	!
[1852]	33	! 1849 2016-04-08 11:33:18Z hoffmann
[1851]	34	! Adapted for modularization of microphysics
[1852]	35	!
[1827]	36	! 1826 2016-04-07 12:01:39Z maronga
	37	! Further modularization.
	38	!
[1789]	39	! 1788 2016-03-10 11:01:04Z maronga
	40	! Added new albedo class for pavements / roads.
	41	!
[1784]	42	! 1783 2016-03-06 18:36:17Z raasch
	43	! palm-netcdf-module removed in order to avoid a circular module dependency,
	44	! netcdf-variables moved to netcdf-module, new routine netcdf_handle_error_rad
	45	! added
	46	!
[1758]	47	! 1757 2016-02-22 15:49:32Z maronga
	48	! Added parameter unscheduled_radiation_calls. Bugfix: interpolation of sounding
	49	! profiles for pressure and temperature above the LES domain.
	50	!
[1710]	51	! 1709 2015-11-04 14:47:01Z maronga
	52	! Bugfix: set initial value for rrtm_lwuflx_dt to zero, small formatting
	53	! corrections
	54	!
[1702]	55	! 1701 2015-11-02 07:43:04Z maronga
	56	! Bugfixes: wrong index for output of timeseries, setting of nz_snd_end
	57	!
[1692]	58	! 1691 2015-10-26 16:17:44Z maronga
	59	! Added option for spin-up runs without radiation (skip_time_do_radiation). Bugfix
	60	! in calculation of pressure profiles. Bugfix in calculation of trace gas profiles.
	61	! Added output of radiative heating rates.
	62	!
[1683]	63	! 1682 2015-10-07 23:56:08Z knoop
	64	! Code annotations made doxygen readable
	65	!
[1607]	66	! 1606 2015-06-29 10:43:37Z maronga
	67	! Added preprocessor directive __netcdf to allow for compiling without netCDF.
	68	! Note, however, that RRTMG cannot be used without netCDF.
	69	!
[1591]	70	! 1590 2015-05-08 13:56:27Z maronga
	71	! Bugfix: definition of character strings requires same length for all elements
	72	!
[1588]	73	! 1587 2015-05-04 14:19:01Z maronga
	74	! Added albedo class for snow
	75	!
[1586]	76	! 1585 2015-04-30 07:05:52Z maronga
	77	! Added support for RRTMG
	78	!
[1572]	79	! 1571 2015-03-12 16:12:49Z maronga
	80	! Added missing KIND attribute. Removed upper-case variable names
	81	!
[1552]	82	! 1551 2015-03-03 14:18:16Z maronga
	83	! Added support for data output. Various variables have been renamed. Added
	84	! interface for different radiation schemes (currently: clear-sky, constant, and
	85	! RRTM (not yet implemented).
	86	!
[1497]	87	! 1496 2014-12-02 17:25:50Z maronga
	88	! Initial revision
	89	!
[1496]	90	!
	91	! Description:
	92	! ------------
[1682]	93	!> Radiation models and interfaces
[1826]	94	!> @todo move variable definitions used in radiation_init only to the subroutine
[1682]	95	!> as they are no longer required after initialization.
	96	!> @todo Output of full column vertical profiles used in RRTMG
	97	!> @todo Output of other rrtm arrays (such as volume mixing ratios)
	98	!> @todo Adapt for use with topography
	99	!>
	100	!> @note Many variables have a leading dummy dimension (0:0) in order to
	101	!> match the assume-size shape expected by the RRTMG model.
[1496]	102	!------------------------------------------------------------------------------!
[1682]	103	MODULE radiation_model_mod
	104
[1496]	105	USE arrays_3d, &
[1691]	106	ONLY: dzw, hyp, pt, q, ql, zw
[1496]	107
[1585]	108	USE cloud_parameters, &
[1849]	109	ONLY: cp, l_d_cp, rho_l
[1585]	110
	111	USE constants, &
	112	ONLY: pi
	113
[1496]	114	USE control_parameters, &
[1585]	115	ONLY: cloud_droplets, cloud_physics, g, initializing_actions, &
[1691]	116	large_scale_forcing, lsf_surf, phi, pt_surface, rho_surface, &
[1585]	117	surface_pressure, time_since_reference_point
[1496]	118
	119	USE indices, &
[1585]	120	ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb_s_inner, nzb, nzt
[1496]	121
	122	USE kinds
	123
[1849]	124	USE microphysics_mod, &
	125	ONLY: nc_const, sigma_gc
	126
[1606]	127	#if defined ( __netcdf )
[1783]	128	USE NETCDF
[1606]	129	#endif
[1585]	130
	131	#if defined ( __rrtmg )
[1788]	132
[1585]	133	USE parrrsw, &
	134	ONLY: naerec, nbndsw
[1551]	135
[1585]	136	USE parrrtm, &
	137	ONLY: nbndlw
	138
	139	USE rrtmg_lw_init, &
	140	ONLY: rrtmg_lw_ini
	141
	142	USE rrtmg_sw_init, &
	143	ONLY: rrtmg_sw_ini
	144
	145	USE rrtmg_lw_rad, &
	146	ONLY: rrtmg_lw
	147
	148	USE rrtmg_sw_rad, &
	149	ONLY: rrtmg_sw
	150	#endif
	151
	152
	153
[1496]	154	IMPLICIT NONE
	155
[1585]	156	CHARACTER(10) :: radiation_scheme = 'clear-sky' ! 'constant', 'clear-sky', or 'rrtmg'
[1551]	157
[1585]	158	!
	159	!-- Predefined Land surface classes (albedo_type) after Briegleb (1992)
[1788]	160	CHARACTER(37), DIMENSION(0:17), PARAMETER :: albedo_type_name = (/ &
[1590]	161	'user defined ', & ! 0
	162	'ocean ', & ! 1
	163	'mixed farming, tall grassland ', & ! 2
	164	'tall/medium grassland ', & ! 3
	165	'evergreen shrubland ', & ! 4
	166	'short grassland/meadow/shrubland ', & ! 5
	167	'evergreen needleleaf forest ', & ! 6
	168	'mixed deciduous evergreen forest ', & ! 7
	169	'deciduous forest ', & ! 8
	170	'tropical evergreen broadleaved forest', & ! 9
	171	'medium/tall grassland/woodland ', & ! 10
	172	'desert, sandy ', & ! 11
	173	'desert, rocky ', & ! 12
	174	'tundra ', & ! 13
	175	'land ice ', & ! 14
	176	'sea ice ', & ! 15
[1788]	177	'snow ', & ! 16
	178	'pavement/roads ' & ! 17
[1585]	179	/)
[1496]	180
[1682]	181	INTEGER(iwp) :: albedo_type = 5, & !< Albedo surface type (default: short grassland)
	182	day, & !< current day of the year
	183	day_init = 172, & !< day of the year at model start (21/06)
	184	dots_rad = 0 !< starting index for timeseries output
[1496]	185
[1757]	186	LOGICAL :: unscheduled_radiation_calls = .TRUE., & !< flag parameter indicating whether additional calls of the radiation code are allowed
	187	constant_albedo = .FALSE., & !< flag parameter indicating whether the albedo may change depending on zenith
	188	force_radiation_call = .FALSE., & !< flag parameter for unscheduled radiation calls
	189	lw_radiation = .TRUE., & !< flag parameter indicating whether longwave radiation shall be calculated
	190	radiation = .FALSE., & !< flag parameter indicating whether the radiation model is used
	191	sun_up = .TRUE., & !< flag parameter indicating whether the sun is up or down
	192	sw_radiation = .TRUE. !< flag parameter indicing whether shortwave radiation shall be calculated
[1585]	193
[1496]	194
[1691]	195	REAL(wp), PARAMETER :: d_seconds_hour = 0.000277777777778_wp, & !< inverse of seconds per hour (1/3600)
	196	d_hours_day = 0.0416666666667_wp, & !< inverse of hours per day (1/24)
	197	sigma_sb = 5.67037321E-8_wp, & !< Stefan-Boltzmann constant
	198	solar_constant = 1368.0_wp !< solar constant at top of atmosphere
[1585]	199
[1691]	200	REAL(wp) :: albedo = 9999999.9_wp, & !< NAMELIST alpha
	201	albedo_lw_dif = 9999999.9_wp, & !< NAMELIST aldif
	202	albedo_lw_dir = 9999999.9_wp, & !< NAMELIST aldir
	203	albedo_sw_dif = 9999999.9_wp, & !< NAMELIST asdif
	204	albedo_sw_dir = 9999999.9_wp, & !< NAMELIST asdir
	205	decl_1, & !< declination coef. 1
	206	decl_2, & !< declination coef. 2
	207	decl_3, & !< declination coef. 3
	208	dt_radiation = 0.0_wp, & !< radiation model timestep
	209	emissivity = 0.98_wp, & !< NAMELIST surface emissivity
	210	lambda = 0.0_wp, & !< longitude in degrees
	211	lon = 0.0_wp, & !< longitude in radians
	212	lat = 0.0_wp, & !< latitude in radians
	213	net_radiation = 0.0_wp, & !< net radiation at surface
	214	skip_time_do_radiation = 0.0_wp, & !< Radiation model is not called before this time
	215	sky_trans, & !< sky transmissivity
	216	time_radiation = 0.0_wp, & !< time since last call of radiation code
	217	time_utc, & !< current time in UTC
	218	time_utc_init = 43200.0_wp !< UTC time at model start (noon)
	219
[1682]	220	REAL(wp), DIMENSION(0:0) :: zenith !< solar zenith angle
[1585]	221
[1496]	222	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
[1682]	223	alpha, & !< surface broadband albedo (used for clear-sky scheme)
[1709]	224	rad_lw_out_change_0, & !< change in LW out due to change in surface temperature
[1682]	225	rad_net, & !< net radiation at the surface
	226	rad_net_av !< average of rad_net
[1496]	227
[1585]	228	!
	229	!-- Land surface albedos for solar zenith angle of 60Â° after Briegleb (1992)
	230	!-- (shortwave, longwave, broadband): sw, lw, bb,
[1788]	231	REAL(wp), DIMENSION(0:2,1:17), PARAMETER :: albedo_pars = RESHAPE( (/&
[1585]	232	0.06_wp, 0.06_wp, 0.06_wp, & ! 1
	233	0.09_wp, 0.28_wp, 0.19_wp, & ! 2
	234	0.11_wp, 0.33_wp, 0.23_wp, & ! 3
	235	0.11_wp, 0.33_wp, 0.23_wp, & ! 4
	236	0.14_wp, 0.34_wp, 0.25_wp, & ! 5
	237	0.06_wp, 0.22_wp, 0.14_wp, & ! 6
	238	0.06_wp, 0.27_wp, 0.17_wp, & ! 7
	239	0.06_wp, 0.31_wp, 0.19_wp, & ! 8
	240	0.06_wp, 0.22_wp, 0.14_wp, & ! 9
	241	0.06_wp, 0.28_wp, 0.18_wp, & ! 10
	242	0.35_wp, 0.51_wp, 0.43_wp, & ! 11
	243	0.24_wp, 0.40_wp, 0.32_wp, & ! 12
	244	0.10_wp, 0.27_wp, 0.19_wp, & ! 13
	245	0.90_wp, 0.65_wp, 0.77_wp, & ! 14
[1587]	246	0.90_wp, 0.65_wp, 0.77_wp, & ! 15
[1788]	247	0.95_wp, 0.70_wp, 0.82_wp, & ! 16
	248	0.08_wp, 0.08_wp, 0.08_wp & ! 17
	249	/), (/ 3, 17 /) )
[1496]	250
[1585]	251	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: &
[1691]	252	rad_lw_cs_hr, & !< longwave clear sky radiation heating rate (K/s)
	253	rad_lw_cs_hr_av, & !< average of rad_lw_cs_hr
	254	rad_lw_hr, & !< longwave radiation heating rate (K/s)
	255	rad_lw_hr_av, & !< average of rad_sw_hr
	256	rad_lw_in, & !< incoming longwave radiation (W/m2)
	257	rad_lw_in_av, & !< average of rad_lw_in
	258	rad_lw_out, & !< outgoing longwave radiation (W/m2)
	259	rad_lw_out_av, & !< average of rad_lw_out
	260	rad_sw_cs_hr, & !< shortwave clear sky radiation heating rate (K/s)
	261	rad_sw_cs_hr_av, & !< average of rad_sw_cs_hr
	262	rad_sw_hr, & !< shortwave radiation heating rate (K/s)
	263	rad_sw_hr_av, & !< average of rad_sw_hr
[1682]	264	rad_sw_in, & !< incoming shortwave radiation (W/m2)
	265	rad_sw_in_av, & !< average of rad_sw_in
	266	rad_sw_out, & !< outgoing shortwave radiation (W/m2)
[1691]	267	rad_sw_out_av !< average of rad_sw_out
[1585]	268
[1691]	269
[1585]	270	!
	271	!-- Variables and parameters used in RRTMG only
	272	#if defined ( __rrtmg )
[1682]	273	CHARACTER(LEN=12) :: rrtm_input_file = "RAD_SND_DATA" !< name of the NetCDF input file (sounding data)
[1585]	274
	275
	276	!
	277	!-- Flag parameters for RRTMGS (should not be changed)
[1682]	278	INTEGER(iwp), PARAMETER :: rrtm_inflglw = 2, & !< flag for lw cloud optical properties (0,1,2)
	279	rrtm_iceflglw = 0, & !< flag for lw ice particle specifications (0,1,2,3)
	280	rrtm_liqflglw = 1, & !< flag for lw liquid droplet specifications
	281	rrtm_inflgsw = 2, & !< flag for sw cloud optical properties (0,1,2)
	282	rrtm_iceflgsw = 0, & !< flag for sw ice particle specifications (0,1,2,3)
	283	rrtm_liqflgsw = 1 !< flag for sw liquid droplet specifications
[1585]	284
	285	!
	286	!-- The following variables should be only changed with care, as this will
	287	!-- require further setting of some variables, which is currently not
	288	!-- implemented (aerosols, ice phase).
[1682]	289	INTEGER(iwp) :: nzt_rad, & !< upper vertical limit for radiation calculations
	290	rrtm_icld = 0, & !< cloud flag (0: clear sky column, 1: cloudy column)
	291	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)
[1691]	292	rrtm_idrv = 1 !< longwave upward flux calculation option (0,1)
[1585]	293
[1788]	294	INTEGER(iwp) :: nc_stat !< local variable for storin the result of netCDF calls for error message handling
	295
[1682]	296	LOGICAL :: snd_exists = .FALSE. !< flag parameter to check whether a user-defined input files exists
[1585]	297
[1691]	298	REAL(wp), PARAMETER :: mol_mass_air_d_wv = 1.607793_wp !< molecular weight dry air / water vapor
[1585]	299
[1682]	300	REAL(wp), DIMENSION(:), ALLOCATABLE :: hyp_snd, & !< hypostatic pressure from sounding data (hPa)
	301	q_snd, & !< specific humidity from sounding data (kg/kg) - dummy at the moment
	302	rrtm_tsfc, & !< dummy array for storing surface temperature
	303	t_snd !< actual temperature from sounding data (hPa)
[1585]	304
[1691]	305	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: aldif, & !< longwave diffuse albedo solar angle of 60Â°
	306	aldir, & !< longwave direct albedo solar angle of 60Â°
	307	asdif, & !< shortwave diffuse albedo solar angle of 60Â°
	308	asdir, & !< shortwave direct albedo solar angle of 60Â°
	309	rrtm_ccl4vmr, & !< CCL4 volume mixing ratio (g/mol)
	310	rrtm_cfc11vmr, & !< CFC11 volume mixing ratio (g/mol)
	311	rrtm_cfc12vmr, & !< CFC12 volume mixing ratio (g/mol)
	312	rrtm_cfc22vmr, & !< CFC22 volume mixing ratio (g/mol)
	313	rrtm_ch4vmr, & !< CH4 volume mixing ratio
	314	rrtm_cicewp, & !< in-cloud ice water path (g/mÂ²)
	315	rrtm_cldfr, & !< cloud fraction (0,1)
	316	rrtm_cliqwp, & !< in-cloud liquid water path (g/mÂ²)
	317	rrtm_co2vmr, & !< CO2 volume mixing ratio (g/mol)
	318	rrtm_emis, & !< surface emissivity (0-1)
	319	rrtm_h2ovmr, & !< H2O volume mixing ratio
	320	rrtm_n2ovmr, & !< N2O volume mixing ratio
	321	rrtm_o2vmr, & !< O2 volume mixing ratio
	322	rrtm_o3vmr, & !< O3 volume mixing ratio
	323	rrtm_play, & !< pressure layers (hPa, zu-grid)
	324	rrtm_plev, & !< pressure layers (hPa, zw-grid)
	325	rrtm_reice, & !< cloud ice effective radius (microns)
	326	rrtm_reliq, & !< cloud water drop effective radius (microns)
	327	rrtm_tlay, & !< actual temperature (K, zu-grid)
	328	rrtm_tlev, & !< actual temperature (K, zw-grid)
	329	rrtm_lwdflx, & !< RRTM output of incoming longwave radiation flux (W/m2)
	330	rrtm_lwdflxc, & !< RRTM output of outgoing clear sky longwave radiation flux (W/m2)
	331	rrtm_lwuflx, & !< RRTM output of outgoing longwave radiation flux (W/m2)
	332	rrtm_lwuflxc, & !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
	333	rrtm_lwuflx_dt, & !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
	334	rrtm_lwuflxc_dt,& !< RRTM output of outgoing clear sky longwave radiation flux (W/m2)
	335	rrtm_lwhr, & !< RRTM output of longwave radiation heating rate (K/d)
	336	rrtm_lwhrc, & !< RRTM output of incoming longwave clear sky radiation heating rate (K/d)
	337	rrtm_swdflx, & !< RRTM output of incoming shortwave radiation flux (W/m2)
	338	rrtm_swdflxc, & !< RRTM output of outgoing clear sky shortwave radiation flux (W/m2)
	339	rrtm_swuflx, & !< RRTM output of outgoing shortwave radiation flux (W/m2)
	340	rrtm_swuflxc, & !< RRTM output of incoming clear sky shortwave radiation flux (W/m2)
	341	rrtm_swhr, & !< RRTM output of shortwave radiation heating rate (K/d)
	342	rrtm_swhrc !< RRTM output of incoming shortwave clear sky radiation heating rate (K/d)
[1585]	343
	344	!
	345	!-- Definition of arrays that are currently not used for calling RRTMG (due to setting of flag parameters)
[1682]	346	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: rad_lw_cs_in, & !< incoming clear sky longwave radiation (W/m2) (not used)
	347	rad_lw_cs_out, & !< outgoing clear sky longwave radiation (W/m2) (not used)
	348	rad_sw_cs_in, & !< incoming clear sky shortwave radiation (W/m2) (not used)
	349	rad_sw_cs_out, & !< outgoing clear sky shortwave radiation (W/m2) (not used)
	350	rrtm_aldif, & !< surface albedo for longwave diffuse radiation
	351	rrtm_aldir, & !< surface albedo for longwave direct radiation
	352	rrtm_asdif, & !< surface albedo for shortwave diffuse radiation
	353	rrtm_asdir, & !< surface albedo for shortwave direct radiation
	354	rrtm_lw_tauaer, & !< lw aerosol optical depth
	355	rrtm_lw_taucld, & !< lw in-cloud optical depth
	356	rrtm_sw_taucld, & !< sw in-cloud optical depth
	357	rrtm_sw_ssacld, & !< sw in-cloud single scattering albedo
	358	rrtm_sw_asmcld, & !< sw in-cloud asymmetry parameter
	359	rrtm_sw_fsfcld, & !< sw in-cloud forward scattering fraction
	360	rrtm_sw_tauaer, & !< sw aerosol optical depth
	361	rrtm_sw_ssaaer, & !< sw aerosol single scattering albedo
	362	rrtm_sw_asmaer, & !< sw aerosol asymmetry parameter
	363	rrtm_sw_ecaer !< sw aerosol optical detph at 0.55 microns (rrtm_iaer = 6 only)
[1691]	364
[1585]	365	#endif
	366
[1826]	367	INTERFACE radiation_check_data_output
	368	MODULE PROCEDURE radiation_check_data_output
	369	END INTERFACE radiation_check_data_output
[1496]	370
[1826]	371	INTERFACE radiation_check_data_output_pr
	372	MODULE PROCEDURE radiation_check_data_output_pr
	373	END INTERFACE radiation_check_data_output_pr
	374
	375	INTERFACE radiation_check_parameters
	376	MODULE PROCEDURE radiation_check_parameters
	377	END INTERFACE radiation_check_parameters
	378
[1551]	379	INTERFACE radiation_clearsky
	380	MODULE PROCEDURE radiation_clearsky
	381	END INTERFACE radiation_clearsky
[1853]	382
	383	INTERFACE radiation_constant
	384	MODULE PROCEDURE radiation_constant
	385	END INTERFACE radiation_constant
	386
[1826]	387	INTERFACE radiation_header
	388	MODULE PROCEDURE radiation_header
	389	END INTERFACE radiation_header
	390
	391	INTERFACE radiation_init
	392	MODULE PROCEDURE radiation_init
	393	END INTERFACE radiation_init
[1496]	394
[1826]	395	INTERFACE radiation_parin
	396	MODULE PROCEDURE radiation_parin
	397	END INTERFACE radiation_parin
	398
[1585]	399	INTERFACE radiation_rrtmg
	400	MODULE PROCEDURE radiation_rrtmg
	401	END INTERFACE radiation_rrtmg
[1551]	402
[1585]	403	INTERFACE radiation_tendency
	404	MODULE PROCEDURE radiation_tendency
	405	MODULE PROCEDURE radiation_tendency_ij
	406	END INTERFACE radiation_tendency
[1551]	407
[1496]	408	SAVE
	409
	410	PRIVATE
	411
[1826]	412	!
	413	!-- Public functions
	414	PUBLIC radiation_check_data_output, radiation_check_data_output_pr, &
[1853]	415	radiation_check_parameters, radiation_clearsky, radiation_constant, &
	416	radiation_header, radiation_init, radiation_parin, radiation_rrtmg, &
	417	radiation_tendency
[1826]	418
	419	!
	420	!-- Public variables and constants
	421	PUBLIC dots_rad, dt_radiation, force_radiation_call, &
	422	rad_net, rad_net_av, radiation, radiation_scheme, rad_lw_in, &
	423	rad_lw_in_av, rad_lw_out, rad_lw_out_av, rad_lw_out_change_0, &
	424	rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av, rad_sw_in, &
	425	rad_sw_in_av, rad_sw_out, rad_sw_out_av, rad_sw_cs_hr, &
	426	rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av, sigma_sb, &
	427	skip_time_do_radiation, time_radiation, unscheduled_radiation_calls
[1496]	428
[1691]	429
[1585]	430	#if defined ( __rrtmg )
[1709]	431	PUBLIC rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir, rrtm_idrv
[1585]	432	#endif
[1496]	433
	434	CONTAINS
	435
	436	!------------------------------------------------------------------------------!
	437	! Description:
	438	! ------------
[1826]	439	!> Check data output for radiation model
	440	!------------------------------------------------------------------------------!
	441	SUBROUTINE radiation_check_data_output( var, unit, i, ilen, k )
	442
	443
	444	USE control_parameters, &
	445	ONLY: data_output, message_string
	446
	447	IMPLICIT NONE
	448
	449	CHARACTER (LEN=*) :: unit !<
	450	CHARACTER (LEN=*) :: var !<
	451
	452	INTEGER(iwp) :: i
	453	INTEGER(iwp) :: ilen
	454	INTEGER(iwp) :: k
	455
	456	SELECT CASE ( TRIM( var ) )
	457
	458	CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_lw_cs_hr', 'rad_lw_hr', &
	459	'rad_sw_in', 'rad_sw_out', 'rad_sw_cs_hr', 'rad_sw_hr' )
	460	IF ( .NOT. radiation .OR. radiation_scheme /= 'rrtmg' ) THEN
	461	message_string = '"output of "' // TRIM( var ) // '" requi' // &
	462	'res radiation = .TRUE. and ' // &
	463	'radiation_scheme = "rrtmg"'
	464	CALL message( 'check_parameters', 'PA0406', 1, 2, 0, 6, 0 )
	465	ENDIF
	466	unit = 'W/m2'
	467
	468	CASE ( 'rad_net', 'rrtm_aldif', 'rrtm_aldir', 'rrtm_asdif', &
	469	'rrtm_asdir*' )
	470	IF ( k == 0 .OR. data_output(i)(ilen-2:ilen) /= '_xy' ) THEN
	471	message_string = 'illegal value for data_output: "' // &
	472	TRIM( var ) // '" & only 2d-horizontal ' // &
	473	'cross sections are allowed for this value'
	474	CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
	475	ENDIF
	476	IF ( .NOT. radiation .OR. radiation_scheme /= "rrtmg" ) THEN
	477	IF ( TRIM( var ) == 'rrtm_aldif*' .OR. &
	478	TRIM( var ) == 'rrtm_aldir*' .OR. &
	479	TRIM( var ) == 'rrtm_asdif*' .OR. &
	480	TRIM( var ) == 'rrtm_asdir*' ) &
	481	THEN
	482	message_string = 'output of "' // TRIM( var ) // '" require'&
	483	// 's radiation = .TRUE. and radiation_sch'&
	484	// 'eme = "rrtmg"'
	485	CALL message( 'check_parameters', 'PA0409', 1, 2, 0, 6, 0 )
	486	ENDIF
	487	ENDIF
	488
	489	IF ( TRIM( var ) == 'rad_net*' ) unit = 'W/m2'
	490	IF ( TRIM( var ) == 'rrtm_aldif*' ) unit = ''
	491	IF ( TRIM( var ) == 'rrtm_aldir*' ) unit = ''
	492	IF ( TRIM( var ) == 'rrtm_asdif*' ) unit = ''
	493	IF ( TRIM( var ) == 'rrtm_asdir*' ) unit = ''
	494
	495	CASE DEFAULT
	496	unit = 'illegal'
	497
	498	END SELECT
	499
	500
	501	END SUBROUTINE radiation_check_data_output
	502
	503	!------------------------------------------------------------------------------!
	504	! Description:
	505	! ------------
	506	!> Check data output of profiles for radiation model
	507	!------------------------------------------------------------------------------!
	508	SUBROUTINE radiation_check_data_output_pr( variable, var_count, unit, dopr_unit )
	509
	510	USE arrays_3d, &
	511	ONLY: zu
	512
	513	USE control_parameters, &
	514	ONLY: data_output_pr, message_string
	515
	516	USE indices
	517
	518	USE profil_parameter
	519
	520	USE statistics
	521
	522	IMPLICIT NONE
	523
	524	CHARACTER (LEN=*) :: unit !<
	525	CHARACTER (LEN=*) :: variable !<
	526	CHARACTER (LEN=*) :: dopr_unit !< local value of dopr_unit
	527
	528	INTEGER(iwp) :: user_pr_index !<
	529	INTEGER(iwp) :: var_count !<
	530
	531	SELECT CASE ( TRIM( variable ) )
	532
	533	CASE ( 'rad_net' )
	534	IF ( ( .NOT. radiation ) .OR. radiation_scheme == 'constant' )&
	535	THEN
	536	message_string = 'data_output_pr = ' // &
	537	TRIM( data_output_pr(var_count) ) // ' is' // &
	538	'not available for radiation = .FALSE. or ' //&
	539	'radiation_scheme = "constant"'
	540	CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
	541	ELSE
	542	dopr_index(var_count) = 101
	543	dopr_unit = 'W/m2'
	544	hom(:,2,101,:) = SPREAD( zw, 2, statistic_regions+1 )
	545	unit = dopr_unit
	546	ENDIF
	547
	548	CASE ( 'rad_lw_in' )
	549	IF ( ( .NOT. radiation) .OR. radiation_scheme == 'constant' ) &
	550	THEN
	551	message_string = 'data_output_pr = ' // &
	552	TRIM( data_output_pr(var_count) ) // ' is' // &
	553	'not available for radiation = .FALSE. or ' //&
	554	'radiation_scheme = "constant"'
	555	CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
	556	ELSE
	557	dopr_index(var_count) = 102
	558	dopr_unit = 'W/m2'
	559	hom(:,2,102,:) = SPREAD( zw, 2, statistic_regions+1 )
	560	unit = dopr_unit
	561	ENDIF
	562
	563	CASE ( 'rad_lw_out' )
	564	IF ( ( .NOT. radiation ) .OR. radiation_scheme == 'constant' ) &
	565	THEN
	566	message_string = 'data_output_pr = ' // &
	567	TRIM( data_output_pr(var_count) ) // ' is' // &
	568	'not available for radiation = .FALSE. or ' //&
	569	'radiation_scheme = "constant"'
	570	CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
	571	ELSE
	572	dopr_index(var_count) = 103
	573	dopr_unit = 'W/m2'
	574	hom(:,2,103,:) = SPREAD( zw, 2, statistic_regions+1 )
	575	unit = dopr_unit
	576	ENDIF
	577
	578	CASE ( 'rad_sw_in' )
	579	IF ( ( .NOT. radiation ) .OR. radiation_scheme == 'constant' ) &
	580	THEN
	581	message_string = 'data_output_pr = ' // &
	582	TRIM( data_output_pr(var_count) ) // ' is' // &
	583	'not available for radiation = .FALSE. or ' //&
	584	'radiation_scheme = "constant"'
	585	CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
	586	ELSE
	587	dopr_index(var_count) = 104
	588	dopr_unit = 'W/m2'
	589	hom(:,2,104,:) = SPREAD( zw, 2, statistic_regions+1 )
	590	unit = dopr_unit
	591	ENDIF
	592
	593	CASE ( 'rad_sw_out')
	594	IF ( ( .NOT. radiation ) .OR. radiation_scheme == 'constant' )&
	595	THEN
	596	message_string = 'data_output_pr = ' // &
	597	TRIM( data_output_pr(var_count) ) // ' is' // &
	598	'not available for radiation = .FALSE. or ' //&
	599	'radiation_scheme = "constant"'
	600	CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
	601	ELSE
	602	dopr_index(var_count) = 105
	603	dopr_unit = 'W/m2'
	604	hom(:,2,105,:) = SPREAD( zw, 2, statistic_regions+1 )
	605	unit = dopr_unit
	606	ENDIF
	607
	608	CASE ( 'rad_lw_cs_hr' )
	609	IF ( ( .NOT. radiation ) .OR. radiation_scheme /= 'rrtmg' ) &
	610	THEN
	611	message_string = 'data_output_pr = ' // &
	612	TRIM( data_output_pr(var_count) ) // ' is' // &
	613	'not available for radiation = .FALSE. or ' //&
	614	'radiation_scheme /= "rrtmg"'
	615	CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
	616	ELSE
	617	dopr_index(var_count) = 106
	618	dopr_unit = 'K/h'
	619	hom(:,2,106,:) = SPREAD( zu, 2, statistic_regions+1 )
	620	unit = dopr_unit
	621	ENDIF
	622
	623	CASE ( 'rad_lw_hr' )
	624	IF ( ( .NOT. radiation ) .OR. radiation_scheme /= 'rrtmg' ) &
	625	THEN
	626	message_string = 'data_output_pr = ' // &
	627	TRIM( data_output_pr(var_count) ) // ' is' // &
	628	'not available for radiation = .FALSE. or ' //&
	629	'radiation_scheme /= "rrtmg"'
	630	CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
	631	ELSE
	632	dopr_index(var_count) = 107
	633	dopr_unit = 'K/h'
	634	hom(:,2,107,:) = SPREAD( zu, 2, statistic_regions+1 )
	635	unit = dopr_unit
	636	ENDIF
	637
	638	CASE ( 'rad_sw_cs_hr' )
	639	IF ( ( .NOT. radiation ) .OR. radiation_scheme /= 'rrtmg' ) &
	640	THEN
	641	message_string = 'data_output_pr = ' // &
	642	TRIM( data_output_pr(var_count) ) // ' is' // &
	643	'not available for radiation = .FALSE. or ' //&
	644	'radiation_scheme /= "rrtmg"'
	645	CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
	646	ELSE
	647	dopr_index(var_count) = 108
	648	dopr_unit = 'K/h'
	649	hom(:,2,108,:) = SPREAD( zu, 2, statistic_regions+1 )
	650	unit = dopr_unit
	651	ENDIF
	652
	653	CASE ( 'rad_sw_hr' )
	654	IF ( ( .NOT. radiation ) .OR. radiation_scheme /= 'rrtmg' ) &
	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 /= "rrtmg"'
	660	CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
	661	ELSE
	662	dopr_index(var_count) = 109
	663	dopr_unit = 'K/h'
	664	hom(:,2,109,:) = SPREAD( zu, 2, statistic_regions+1 )
	665	unit = dopr_unit
	666	ENDIF
	667
	668
	669	CASE DEFAULT
	670	unit = 'illegal'
	671
	672	END SELECT
	673
	674
	675	END SUBROUTINE radiation_check_data_output_pr
	676
	677
	678	!------------------------------------------------------------------------------!
	679	! Description:
	680	! ------------
	681	!> Check parameters routine for radiation model
	682	!------------------------------------------------------------------------------!
	683	SUBROUTINE radiation_check_parameters
	684
	685	USE control_parameters, &
	686	ONLY: message_string, topography
	687
	688
	689	IMPLICIT NONE
	690
	691	IF ( radiation_scheme /= 'constant' .AND. &
	692	radiation_scheme /= 'clear-sky' .AND. &
	693	radiation_scheme /= 'rrtmg' ) THEN
	694	message_string = 'unknown radiation_scheme = '// &
	695	TRIM( radiation_scheme )
	696	CALL message( 'check_parameters', 'PA0405', 1, 2, 0, 6, 0 )
	697	ELSEIF ( radiation_scheme == 'rrtmg' ) THEN
	698	#if ! defined ( __rrtmg )
	699	message_string = 'radiation_scheme = "rrtmg" requires ' // &
	700	'compilation of PALM with pre-processor ' // &
	701	'directive -D__rrtmg'
	702	CALL message( 'check_parameters', 'PA0407', 1, 2, 0, 6, 0 )
	703	#endif
	704	#if defined ( __rrtmg ) && ! defined( __netcdf )
	705	message_string = 'radiation_scheme = "rrtmg" requires ' // &
	706	'the use of NetCDF (preprocessor directive ' // &
	707	'-D__netcdf'
	708	CALL message( 'check_parameters', 'PA0412', 1, 2, 0, 6, 0 )
	709	#endif
	710
	711	ENDIF
	712
	713	IF ( albedo_type == 0 .AND. albedo == 9999999.9_wp .AND. &
	714	radiation_scheme == 'clear-sky') THEN
	715	message_string = 'radiation_scheme = "clear-sky" in combination' // &
	716	'with albedo_type = 0 requires setting of albedo'// &
	717	' /= 9999999.9'
	718	CALL message( 'check_parameters', 'PA0410', 1, 2, 0, 6, 0 )
	719	ENDIF
	720
	721	IF ( albedo_type == 0 .AND. radiation_scheme == 'rrtmg' .AND. &
	722	( albedo_lw_dif == 9999999.9_wp .OR. albedo_lw_dir == 9999999.9_wp&
	723	.OR. albedo_sw_dif == 9999999.9_wp .OR. albedo_sw_dir == 9999999.9_wp&
	724	) ) THEN
	725	message_string = 'radiation_scheme = "rrtmg" in combination' // &
	726	'with albedo_type = 0 requires setting of ' // &
	727	'albedo_lw_dif /= 9999999.9' // &
	728	'albedo_lw_dir /= 9999999.9' // &
	729	'albedo_sw_dif /= 9999999.9 and' // &
	730	'albedo_sw_dir /= 9999999.9'
	731	CALL message( 'check_parameters', 'PA0411', 1, 2, 0, 6, 0 )
	732	ENDIF
	733
	734	IF ( topography /= 'flat' ) THEN
	735	message_string = 'radiation scheme cannot be used ' // &
	736	'in combination with topography /= "flat"'
	737	CALL message( 'check_parameters', 'PA0414', 1, 2, 0, 6, 0 )
	738	ENDIF
	739
	740	END SUBROUTINE radiation_check_parameters
	741
	742
	743	!------------------------------------------------------------------------------!
	744	! Description:
	745	! ------------
[1682]	746	!> Initialization of the radiation model
[1496]	747	!------------------------------------------------------------------------------!
[1826]	748	SUBROUTINE radiation_init
[1496]	749
	750	IMPLICIT NONE
	751
[1585]	752	!
	753	!-- Allocate array for storing the surface net radiation
	754	IF ( .NOT. ALLOCATED ( rad_net ) ) THEN
	755	ALLOCATE ( rad_net(nysg:nyng,nxlg:nxrg) )
	756	rad_net = 0.0_wp
	757	ENDIF
[1496]	758
	759	!
[1709]	760	!-- Allocate array for storing the surface net radiation
	761	IF ( .NOT. ALLOCATED ( rad_lw_out_change_0 ) ) THEN
	762	ALLOCATE ( rad_lw_out_change_0(nysg:nyng,nxlg:nxrg) )
	763	rad_lw_out_change_0 = 0.0_wp
	764	ENDIF
	765
	766	!
[1551]	767	!-- Fix net radiation in case of radiation_scheme = 'constant'
[1585]	768	IF ( radiation_scheme == 'constant' ) THEN
[1551]	769	rad_net = net_radiation
[1853]	770	! radiation = .FALSE.
[1551]	771	!
[1585]	772	!-- Calculate orbital constants
	773	ELSE
[1551]	774	decl_1 = SIN(23.45_wp * pi / 180.0_wp)
	775	decl_2 = 2.0_wp * pi / 365.0_wp
	776	decl_3 = decl_2 * 81.0_wp
[1585]	777	lat = phi * pi / 180.0_wp
	778	lon = lambda * pi / 180.0_wp
	779	ENDIF
	780
	781
[1853]	782	IF ( radiation_scheme == 'constant' ) THEN
	783
	784	IF ( .NOT. ALLOCATED ( rad_lw_out ) ) THEN
	785	ALLOCATE ( rad_lw_out(0:0,nysg:nyng,nxlg:nxrg) )
	786	ENDIF
	787
	788	ENDIF
	789
[1585]	790	IF ( radiation_scheme == 'clear-sky' ) THEN
	791
	792	ALLOCATE ( alpha(nysg:nyng,nxlg:nxrg) )
	793
	794	IF ( .NOT. ALLOCATED ( rad_sw_in ) ) THEN
	795	ALLOCATE ( rad_sw_in(0:0,nysg:nyng,nxlg:nxrg) )
	796	ENDIF
	797	IF ( .NOT. ALLOCATED ( rad_sw_out ) ) THEN
	798	ALLOCATE ( rad_sw_out(0:0,nysg:nyng,nxlg:nxrg) )
	799	ENDIF
	800
	801	IF ( .NOT. ALLOCATED ( rad_sw_in_av ) ) THEN
	802	ALLOCATE ( rad_sw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
	803	ENDIF
	804	IF ( .NOT. ALLOCATED ( rad_sw_out_av ) ) THEN
	805	ALLOCATE ( rad_sw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
	806	ENDIF
	807
	808	IF ( .NOT. ALLOCATED ( rad_lw_in ) ) THEN
	809	ALLOCATE ( rad_lw_in(0:0,nysg:nyng,nxlg:nxrg) )
	810	ENDIF
[1856]	811	IF ( .NOT. ALLOCATED ( rad_lw_out ) ) THEN
	812	ALLOCATE ( rad_lw_out(0:0,nysg:nyng,nxlg:nxrg) )
	813	ENDIF
[1585]	814
	815	IF ( .NOT. ALLOCATED ( rad_lw_in_av ) ) THEN
	816	ALLOCATE ( rad_lw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
	817	ENDIF
	818	IF ( .NOT. ALLOCATED ( rad_lw_out_av ) ) THEN
	819	ALLOCATE ( rad_lw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
	820	ENDIF
	821
	822	rad_sw_in = 0.0_wp
	823	rad_sw_out = 0.0_wp
	824	rad_lw_in = 0.0_wp
	825	rad_lw_out = 0.0_wp
	826
[1496]	827	!
[1585]	828	!-- Overwrite albedo if manually set in parameter file
	829	IF ( albedo_type /= 0 .AND. albedo == 9999999.9_wp ) THEN
	830	albedo = albedo_pars(2,albedo_type)
	831	ENDIF
	832
	833	alpha = albedo
	834
	835	!
	836	!-- Initialization actions for RRTMG
	837	ELSEIF ( radiation_scheme == 'rrtmg' ) THEN
	838	#if defined ( __rrtmg )
	839	!
	840	!-- Allocate albedos
	841	ALLOCATE ( rrtm_aldif(0:0,nysg:nyng,nxlg:nxrg) )
	842	ALLOCATE ( rrtm_aldir(0:0,nysg:nyng,nxlg:nxrg) )
	843	ALLOCATE ( rrtm_asdif(0:0,nysg:nyng,nxlg:nxrg) )
	844	ALLOCATE ( rrtm_asdir(0:0,nysg:nyng,nxlg:nxrg) )
	845	ALLOCATE ( aldif(nysg:nyng,nxlg:nxrg) )
	846	ALLOCATE ( aldir(nysg:nyng,nxlg:nxrg) )
	847	ALLOCATE ( asdif(nysg:nyng,nxlg:nxrg) )
	848	ALLOCATE ( asdir(nysg:nyng,nxlg:nxrg) )
	849
	850	IF ( albedo_type /= 0 ) THEN
	851	IF ( albedo_lw_dif == 9999999.9_wp ) THEN
	852	albedo_lw_dif = albedo_pars(0,albedo_type)
	853	albedo_lw_dir = albedo_lw_dif
	854	ENDIF
	855	IF ( albedo_sw_dif == 9999999.9_wp ) THEN
	856	albedo_sw_dif = albedo_pars(1,albedo_type)
	857	albedo_sw_dir = albedo_sw_dif
	858	ENDIF
	859	ENDIF
	860
	861	aldif(:,:) = albedo_lw_dif
	862	aldir(:,:) = albedo_lw_dir
	863	asdif(:,:) = albedo_sw_dif
	864	asdir(:,:) = albedo_sw_dir
	865	!
	866	!-- Calculate initial values of current (cosine of) the zenith angle and
	867	!-- whether the sun is up
	868	CALL calc_zenith
	869	!
	870	!-- Calculate initial surface albedo
	871	IF ( .NOT. constant_albedo ) THEN
	872	CALL calc_albedo
	873	ELSE
	874	rrtm_aldif(0,:,:) = aldif(:,:)
	875	rrtm_aldir(0,:,:) = aldir(:,:)
	876	rrtm_asdif(0,:,:) = asdif(:,:)
	877	rrtm_asdir(0,:,:) = asdir(:,:)
	878	ENDIF
	879
	880	!
	881	!-- Allocate surface emissivity
	882	ALLOCATE ( rrtm_emis(0:0,1:nbndlw+1) )
	883	rrtm_emis = emissivity
	884
	885	!
	886	!-- Allocate 3d arrays of radiative fluxes and heating rates
	887	IF ( .NOT. ALLOCATED ( rad_sw_in ) ) THEN
	888	ALLOCATE ( rad_sw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	889	rad_sw_in = 0.0_wp
	890	ENDIF
	891
	892	IF ( .NOT. ALLOCATED ( rad_sw_in_av ) ) THEN
	893	ALLOCATE ( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	894	ENDIF
	895
	896	IF ( .NOT. ALLOCATED ( rad_sw_out ) ) THEN
	897	ALLOCATE ( rad_sw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1691]	898	rad_sw_out = 0.0_wp
[1585]	899	ENDIF
	900
	901	IF ( .NOT. ALLOCATED ( rad_sw_out_av ) ) THEN
	902	ALLOCATE ( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	903	ENDIF
	904
[1691]	905	IF ( .NOT. ALLOCATED ( rad_sw_hr ) ) THEN
	906	ALLOCATE ( rad_sw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	907	rad_sw_hr = 0.0_wp
	908	ENDIF
[1585]	909
[1691]	910	IF ( .NOT. ALLOCATED ( rad_sw_hr_av ) ) THEN
	911	ALLOCATE ( rad_sw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	912	rad_sw_hr_av = 0.0_wp
	913	ENDIF
	914
	915	IF ( .NOT. ALLOCATED ( rad_sw_cs_hr ) ) THEN
	916	ALLOCATE ( rad_sw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	917	rad_sw_cs_hr = 0.0_wp
	918	ENDIF
	919
	920	IF ( .NOT. ALLOCATED ( rad_sw_cs_hr_av ) ) THEN
	921	ALLOCATE ( rad_sw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	922	rad_sw_cs_hr_av = 0.0_wp
	923	ENDIF
	924
[1585]	925	IF ( .NOT. ALLOCATED ( rad_lw_in ) ) THEN
	926	ALLOCATE ( rad_lw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	927	rad_lw_in = 0.0_wp
	928	ENDIF
	929
	930	IF ( .NOT. ALLOCATED ( rad_lw_in_av ) ) THEN
	931	ALLOCATE ( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	932	ENDIF
	933
	934	IF ( .NOT. ALLOCATED ( rad_lw_out ) ) THEN
	935	ALLOCATE ( rad_lw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	936	rad_lw_out = 0.0_wp
	937	ENDIF
	938
	939	IF ( .NOT. ALLOCATED ( rad_lw_out_av ) ) THEN
	940	ALLOCATE ( rad_lw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	941	ENDIF
	942
[1691]	943	IF ( .NOT. ALLOCATED ( rad_lw_hr ) ) THEN
	944	ALLOCATE ( rad_lw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	945	rad_lw_hr = 0.0_wp
	946	ENDIF
	947
	948	IF ( .NOT. ALLOCATED ( rad_lw_hr_av ) ) THEN
	949	ALLOCATE ( rad_lw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	950	rad_lw_hr_av = 0.0_wp
	951	ENDIF
	952
	953	IF ( .NOT. ALLOCATED ( rad_lw_cs_hr ) ) THEN
	954	ALLOCATE ( rad_lw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	955	rad_lw_cs_hr = 0.0_wp
	956	ENDIF
	957
	958	IF ( .NOT. ALLOCATED ( rad_lw_cs_hr_av ) ) THEN
	959	ALLOCATE ( rad_lw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	960	rad_lw_cs_hr_av = 0.0_wp
	961	ENDIF
	962
	963	ALLOCATE ( rad_sw_cs_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	964	ALLOCATE ( rad_sw_cs_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1585]	965	rad_sw_cs_in = 0.0_wp
	966	rad_sw_cs_out = 0.0_wp
	967
[1691]	968	ALLOCATE ( rad_lw_cs_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	969	ALLOCATE ( rad_lw_cs_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1585]	970	rad_lw_cs_in = 0.0_wp
	971	rad_lw_cs_out = 0.0_wp
	972
	973	!
	974	!-- Allocate dummy array for storing surface temperature
	975	ALLOCATE ( rrtm_tsfc(1) )
	976
	977	!
	978	!-- Initialize RRTMG
	979	IF ( lw_radiation ) CALL rrtmg_lw_ini ( cp )
	980	IF ( sw_radiation ) CALL rrtmg_sw_ini ( cp )
	981
	982	!
	983	!-- Set input files for RRTMG
	984	INQUIRE(FILE="RAD_SND_DATA", EXIST=snd_exists)
	985	IF ( .NOT. snd_exists ) THEN
	986	rrtm_input_file = "rrtmg_lw.nc"
	987	ENDIF
	988
	989	!
	990	!-- Read vertical layers for RRTMG from sounding data
	991	!-- The routine provides nzt_rad, hyp_snd(1:nzt_rad),
	992	!-- t_snd(nzt+2:nzt_rad), rrtm_play(1:nzt_rad), rrtm_plev(1_nzt_rad+1),
	993	!-- rrtm_tlay(nzt+2:nzt_rad), rrtm_tlev(nzt+2:nzt_rad+1)
	994	CALL read_sounding_data
	995
	996	!
	997	!-- Read trace gas profiles from file. This routine provides
	998	!-- the rrtm_ arrays (1:nzt_rad+1)
	999	CALL read_trace_gas_data
	1000	#endif
[1551]	1001	ENDIF
[1585]	1002
[1551]	1003	!
[1585]	1004	!-- Perform user actions if required
	1005	CALL user_init_radiation
	1006
	1007	!
	1008	!-- Calculate radiative fluxes at model start
	1009	IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
[1853]	1010
	1011	SELECT CASE ( radiation_scheme )
	1012	CASE ( 'rrtmg' )
	1013	CALL radiation_rrtmg
	1014	CASE ( 'clear-sky' )
	1015	CALL radiation_clearsky
	1016	CASE ( 'constant' )
	1017	CALL radiation_constant
	1018	CASE DEFAULT
	1019	END SELECT
	1020
[1585]	1021	ENDIF
	1022
[1496]	1023	RETURN
	1024
[1826]	1025	END SUBROUTINE radiation_init
[1496]	1026
	1027
	1028	!------------------------------------------------------------------------------!
	1029	! Description:
	1030	! ------------
[1682]	1031	!> A simple clear sky radiation model
[1496]	1032	!------------------------------------------------------------------------------!
[1551]	1033	SUBROUTINE radiation_clearsky
[1496]	1034
[1585]	1035
[1496]	1036	IMPLICIT NONE
	1037
[1691]	1038	INTEGER(iwp) :: i, j, k !< loop indices
	1039	REAL(wp) :: exn, & !< Exner functions at surface
[1709]	1040	exn1, & !< Exner functions at first grid level
	1041	pt1 !< potential temperature at first grid level
[1585]	1042
[1496]	1043	!
[1585]	1044	!-- Calculate current zenith angle
	1045	CALL calc_zenith
	1046
	1047	!
	1048	!-- Calculate sky transmissivity
	1049	sky_trans = 0.6_wp + 0.2_wp * zenith(0)
	1050
	1051	!
	1052	!-- Calculate value of the Exner function
	1053	exn = (surface_pressure / 1000.0_wp )**0.286_wp
	1054	!
	1055	!-- Calculate radiation fluxes and net radiation (rad_net) for each grid
	1056	!-- point
[1709]	1057	DO i = nxlg, nxrg
	1058	DO j = nysg, nyng
[1585]	1059	k = nzb_s_inner(j,i)
[1691]	1060
[1709]	1061	exn1 = (hyp(k+1) / 100000.0_wp )**0.286_wp
[1691]	1062
[1585]	1063	rad_sw_in(0,j,i) = solar_constant * sky_trans * zenith(0)
	1064	rad_sw_out(0,j,i) = alpha(j,i) * rad_sw_in(0,j,i)
[1691]	1065	rad_lw_out(0,j,i) = emissivity * sigma_sb * (pt(k,j,i) * exn)**4
[1585]	1066
[1691]	1067	IF ( cloud_physics ) THEN
[1709]	1068	pt1 = pt(k+1,j,i) + l_d_cp / exn1 * ql(k+1,j,i)
	1069	rad_lw_in(0,j,i) = 0.8_wp * sigma_sb * (pt1 * exn1)**4
[1691]	1070	ELSE
[1709]	1071	rad_lw_in(0,j,i) = 0.8_wp * sigma_sb * (pt(k+1,j,i) * exn1)**4
[1691]	1072	ENDIF
	1073
	1074	rad_net(j,i) = rad_sw_in(0,j,i) - rad_sw_out(0,j,i) &
	1075	+ rad_lw_in(0,j,i) - rad_lw_out(0,j,i)
	1076
[1585]	1077	ENDDO
	1078	ENDDO
	1079
	1080	END SUBROUTINE radiation_clearsky
	1081
	1082
	1083	!------------------------------------------------------------------------------!
	1084	! Description:
	1085	! ------------
[1853]	1086	!> This scheme keeps the prescribed net radiation constant during the run
	1087	!------------------------------------------------------------------------------!
	1088	SUBROUTINE radiation_constant
	1089
	1090
	1091	IMPLICIT NONE
	1092
	1093	INTEGER(iwp) :: i, j, k !< loop indices
	1094	REAL(wp) :: exn, & !< Exner functions at surface
	1095	pt1 !< potential temperature at first grid level
	1096
	1097	!
	1098	!-- Calculate value of the Exner function
	1099	exn = (surface_pressure / 1000.0_wp )**0.286_wp
	1100	!
	1101	!-- Prescribe net radiation and estimate a longwave outgoing radiative
	1102	!-- flux (needed in land surface model)
	1103	DO i = nxlg, nxrg
	1104	DO j = nysg, nyng
	1105	k = nzb_s_inner(j,i)
	1106
	1107	rad_net(j,i) = net_radiation
	1108	rad_lw_out(0,j,i) = emissivity * sigma_sb * (pt(k,j,i) * exn)**4
	1109
	1110	ENDDO
	1111	ENDDO
	1112
	1113	END SUBROUTINE radiation_constant
	1114
	1115	!------------------------------------------------------------------------------!
	1116	! Description:
	1117	! ------------
[1826]	1118	!> Header output for radiation model
	1119	!------------------------------------------------------------------------------!
	1120	SUBROUTINE radiation_header ( io )
	1121
	1122
	1123	IMPLICIT NONE
	1124
	1125	INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file
	1126
	1127
	1128
	1129	!
	1130	!-- Write radiation model header
	1131	WRITE( io, 3 )
	1132
	1133	IF ( radiation_scheme == "constant" ) THEN
	1134	WRITE( io, 4 ) net_radiation
	1135	ELSEIF ( radiation_scheme == "clear-sky" ) THEN
	1136	WRITE( io, 5 )
	1137	ELSEIF ( radiation_scheme == "rrtmg" ) THEN
	1138	WRITE( io, 6 )
	1139	IF ( .NOT. lw_radiation ) WRITE( io, 10 )
	1140	IF ( .NOT. sw_radiation ) WRITE( io, 11 )
	1141	ENDIF
	1142
	1143	IF ( albedo_type == 0 ) THEN
	1144	WRITE( io, 7 ) albedo
	1145	ELSE
	1146	WRITE( io, 8 ) TRIM( albedo_type_name(albedo_type) )
	1147	ENDIF
	1148	IF ( constant_albedo ) THEN
	1149	WRITE( io, 9 )
	1150	ENDIF
	1151
	1152	IF ( radiation .AND. radiation_scheme /= 'constant' ) THEN
	1153	WRITE ( io, 1 ) lambda
	1154	WRITE ( io, 2 ) day_init, time_utc_init
	1155	ENDIF
	1156
	1157	WRITE( io, 12 ) dt_radiation
	1158
	1159
	1160	1 FORMAT (' Geograph. longitude : lambda = ',F4.1,' degr')
	1161	2 FORMAT (' Day of the year at model start : day_init = ',I3 &
	1162	/' UTC time at model start : time_utc_init = ',F7.1' s')
	1163	3 FORMAT (//' Radiation model information:'/ &
	1164	' ----------------------------'/)
	1165	4 FORMAT (' --> Using constant net radiation: net_radiation = ', F6.2, &
	1166	// 'W/m**2')
	1167	5 FORMAT (' --> Simple radiation scheme for clear sky is used (no clouds,', &
	1168	' default)')
	1169	6 FORMAT (' --> RRTMG scheme is used')
	1170	7 FORMAT (/' User-specific surface albedo: albedo =', F6.3)
	1171	8 FORMAT (/' Albedo is set for land surface type: ', A)
	1172	9 FORMAT (/' --> Albedo is fixed during the run')
	1173	10 FORMAT (/' --> Longwave radiation is disabled')
	1174	11 FORMAT (/' --> Shortwave radiation is disabled.')
	1175	12 FORMAT (' Timestep: dt_radiation = ', F6.2, ' s')
	1176
	1177
	1178	END SUBROUTINE radiation_header
	1179
	1180
	1181	!------------------------------------------------------------------------------!
	1182	! Description:
	1183	! ------------
	1184	!> Parin for &radiation_par for radiation model
	1185	!------------------------------------------------------------------------------!
	1186	SUBROUTINE radiation_parin
	1187
	1188
	1189	IMPLICIT NONE
	1190
	1191	CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file
	1192
	1193	NAMELIST /radiation_par/ albedo, albedo_type, albedo_lw_dir, &
	1194	albedo_lw_dif, albedo_sw_dir, albedo_sw_dif, &
	1195	constant_albedo, day_init, dt_radiation, &
	1196	lambda, lw_radiation, net_radiation, &
	1197	radiation_scheme, skip_time_do_radiation, &
	1198	sw_radiation, time_utc_init, &
	1199	unscheduled_radiation_calls
	1200
	1201	line = ' '
	1202
	1203	!
	1204	!-- Try to find radiation model package
	1205	REWIND ( 11 )
	1206	line = ' '
	1207	DO WHILE ( INDEX( line, '&radiation_par' ) == 0 )
	1208	READ ( 11, '(A)', END=10 ) line
	1209	ENDDO
	1210	BACKSPACE ( 11 )
	1211
	1212	!
	1213	!-- Read user-defined namelist
	1214	READ ( 11, radiation_par )
	1215
	1216	!
	1217	!-- Set flag that indicates that the radiation model is switched on
	1218	radiation = .TRUE.
	1219
	1220	10 CONTINUE
	1221
	1222
	1223	END SUBROUTINE radiation_parin
	1224
	1225
	1226	!------------------------------------------------------------------------------!
	1227	! Description:
	1228	! ------------
[1682]	1229	!> Implementation of the RRTMG radiation_scheme
[1585]	1230	!------------------------------------------------------------------------------!
	1231	SUBROUTINE radiation_rrtmg
	1232
	1233	USE indices, &
	1234	ONLY: nbgp
	1235
	1236	USE particle_attributes, &
	1237	ONLY: grid_particles, number_of_particles, particles, &
	1238	particle_advection_start, prt_count
	1239
	1240	IMPLICIT NONE
	1241
	1242	#if defined ( __rrtmg )
	1243
[1691]	1244	INTEGER(iwp) :: i, j, k, n !< loop indices
[1585]	1245
[1691]	1246	REAL(wp) :: s_r2, & !< weighted sum over all droplets with r^2
	1247	s_r3 !< weighted sum over all droplets with r^3
[1585]	1248
	1249	!
	1250	!-- Calculate current (cosine of) zenith angle and whether the sun is up
	1251	CALL calc_zenith
	1252	!
	1253	!-- Calculate surface albedo
	1254	IF ( .NOT. constant_albedo ) THEN
	1255	CALL calc_albedo
	1256	ENDIF
	1257
	1258	!
	1259	!-- Prepare input data for RRTMG
	1260
	1261	!
	1262	!-- In case of large scale forcing with surface data, calculate new pressure
	1263	!-- profile. nzt_rad might be modified by these calls and all required arrays
	1264	!-- will then be re-allocated
[1691]	1265	IF ( large_scale_forcing .AND. lsf_surf ) THEN
[1585]	1266	CALL read_sounding_data
	1267	CALL read_trace_gas_data
	1268	ENDIF
	1269	!
	1270	!-- Loop over all grid points
	1271	DO i = nxl, nxr
	1272	DO j = nys, nyn
	1273
	1274	!
	1275	!-- Prepare profiles of temperature and H2O volume mixing ratio
[1691]	1276	rrtm_tlev(0,nzb+1) = pt(nzb,j,i) * ( surface_pressure &
	1277	/ 1000.0_wp )**0.286_wp
[1585]	1278
	1279	DO k = nzb+1, nzt+1
	1280	rrtm_tlay(0,k) = pt(k,j,i) * ( (hyp(k) ) / 100000.0_wp &
[1691]	1281	)*0.286_wp + l_d_cp ql(k,j,i)
	1282	rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * (q(k,j,i) - ql(k,j,i))
[1585]	1283
	1284	ENDDO
	1285
	1286	!
	1287	!-- Avoid temperature/humidity jumps at the top of the LES domain by
	1288	!-- linear interpolation from nzt+2 to nzt+7
	1289	DO k = nzt+2, nzt+7
	1290	rrtm_tlay(0,k) = rrtm_tlay(0,nzt+1) &
	1291	+ ( rrtm_tlay(0,nzt+8) - rrtm_tlay(0,nzt+1) ) &
	1292	/ ( rrtm_play(0,nzt+8) - rrtm_play(0,nzt+1) ) &
	1293	* ( rrtm_play(0,k) - rrtm_play(0,nzt+1) )
	1294
	1295	rrtm_h2ovmr(0,k) = rrtm_h2ovmr(0,nzt+1) &
	1296	+ ( rrtm_h2ovmr(0,nzt+8) - rrtm_h2ovmr(0,nzt+1) )&
	1297	/ ( rrtm_play(0,nzt+8) - rrtm_play(0,nzt+1) )&
	1298	* ( rrtm_play(0,k) - rrtm_play(0,nzt+1) )
	1299
	1300	ENDDO
	1301
	1302	!-- Linear interpolate to zw grid
	1303	DO k = nzb+2, nzt+8
	1304	rrtm_tlev(0,k) = rrtm_tlay(0,k-1) + (rrtm_tlay(0,k) - &
	1305	rrtm_tlay(0,k-1)) &
	1306	/ ( rrtm_play(0,k) - rrtm_play(0,k-1) ) &
	1307	* ( rrtm_plev(0,k) - rrtm_play(0,k-1) )
	1308	ENDDO
	1309
	1310
	1311	!
	1312	!-- Calculate liquid water path and cloud fraction for each column.
	1313	!-- Note that LWP is required in g/mÂ² instead of kg/kg m.
	1314	rrtm_cldfr = 0.0_wp
	1315	rrtm_reliq = 0.0_wp
	1316	rrtm_cliqwp = 0.0_wp
[1691]	1317	rrtm_icld = 0
[1585]	1318
	1319	DO k = nzb+1, nzt+1
[1691]	1320	rrtm_cliqwp(0,k) = ql(k,j,i) * 1000.0_wp * &
	1321	(rrtm_plev(0,k) - rrtm_plev(0,k+1)) &
	1322	* 100.0_wp / g
[1585]	1323
[1691]	1324	IF ( rrtm_cliqwp(0,k) > 0.0_wp ) THEN
[1585]	1325	rrtm_cldfr(0,k) = 1.0_wp
[1691]	1326	IF ( rrtm_icld == 0 ) rrtm_icld = 1
[1585]	1327
	1328	!
	1329	!-- Calculate cloud droplet effective radius
	1330	IF ( cloud_physics ) THEN
[1691]	1331	rrtm_reliq(0,k) = 1.0E6_wp * ( 3.0_wp * ql(k,j,i) &
	1332	* rho_surface &
	1333	/ ( 4.0_wp * pi * nc_const * rho_l ) &
	1334	)**0.33333333333333_wp &
	1335	* EXP( LOG( sigma_gc )**2 )
[1585]	1336
	1337	ELSEIF ( cloud_droplets ) THEN
	1338	number_of_particles = prt_count(k,j,i)
	1339
	1340	IF (number_of_particles <= 0) CYCLE
	1341	particles => grid_particles(k,j,i)%particles(1:number_of_particles)
	1342	s_r2 = 0.0_wp
	1343	s_r3 = 0.0_wp
	1344
	1345	DO n = 1, number_of_particles
	1346	IF ( particles(n)%particle_mask ) THEN
	1347	s_r2 = s_r2 + particles(n)%radius*2 &
	1348	particles(n)%weight_factor
	1349	s_r3 = s_r3 + particles(n)%radius*3 &
	1350	particles(n)%weight_factor
	1351	ENDIF
	1352	ENDDO
	1353
	1354	IF ( s_r2 > 0.0_wp ) rrtm_reliq(0,k) = s_r3 / s_r2
	1355
	1356	ENDIF
	1357
	1358	!
	1359	!-- Limit effective radius
[1691]	1360	IF ( rrtm_reliq(0,k) > 0.0_wp ) THEN
[1585]	1361	rrtm_reliq(0,k) = MAX(rrtm_reliq(0,k),2.5_wp)
	1362	rrtm_reliq(0,k) = MIN(rrtm_reliq(0,k),60.0_wp)
	1363	ENDIF
	1364	ENDIF
	1365	ENDDO
	1366
	1367	!
	1368	!-- Set surface temperature
	1369	rrtm_tsfc = pt(nzb,j,i) * (surface_pressure / 1000.0_wp )**0.286_wp
	1370
	1371	IF ( lw_radiation ) THEN
	1372	CALL rrtmg_lw( 1, nzt_rad , rrtm_icld , rrtm_idrv ,&
	1373	rrtm_play , rrtm_plev , rrtm_tlay , rrtm_tlev ,&
	1374	rrtm_tsfc , rrtm_h2ovmr , rrtm_o3vmr , rrtm_co2vmr ,&
	1375	rrtm_ch4vmr , rrtm_n2ovmr , rrtm_o2vmr , rrtm_cfc11vmr ,&
	1376	rrtm_cfc12vmr , rrtm_cfc22vmr, rrtm_ccl4vmr , rrtm_emis ,&
	1377	rrtm_inflglw , rrtm_iceflglw, rrtm_liqflglw, rrtm_cldfr ,&
	1378	rrtm_lw_taucld , rrtm_cicewp , rrtm_cliqwp , rrtm_reice ,&
	1379	rrtm_reliq , rrtm_lw_tauaer, &
	1380	rrtm_lwuflx , rrtm_lwdflx , rrtm_lwhr , &
[1691]	1381	rrtm_lwuflxc , rrtm_lwdflxc , rrtm_lwhrc , &
	1382	rrtm_lwuflx_dt , rrtm_lwuflxc_dt )
[1585]	1383
[1691]	1384	!
	1385	!-- Save fluxes
[1585]	1386	DO k = nzb, nzt+1
	1387	rad_lw_in(k,j,i) = rrtm_lwdflx(0,k)
	1388	rad_lw_out(k,j,i) = rrtm_lwuflx(0,k)
	1389	ENDDO
	1390
[1691]	1391	!
	1392	!-- Save heating rates (convert from K/d to K/h)
	1393	DO k = nzb+1, nzt+1
	1394	rad_lw_hr(k,j,i) = rrtm_lwhr(0,k) * d_hours_day
	1395	rad_lw_cs_hr(k,j,i) = rrtm_lwhrc(0,k) * d_hours_day
	1396	ENDDO
[1585]	1397
[1709]	1398	!
	1399	!-- Save change in LW heating rate
	1400	rad_lw_out_change_0(j,i) = rrtm_lwuflx_dt(0,nzb)
	1401
[1585]	1402	ENDIF
	1403
	1404	IF ( sw_radiation .AND. sun_up ) THEN
	1405	CALL rrtmg_sw( 1, nzt_rad , rrtm_icld , rrtm_iaer ,&
	1406	rrtm_play , rrtm_plev , rrtm_tlay , rrtm_tlev ,&
	1407	rrtm_tsfc , rrtm_h2ovmr , rrtm_o3vmr , rrtm_co2vmr ,&
	1408	rrtm_ch4vmr , rrtm_n2ovmr , rrtm_o2vmr , rrtm_asdir(:,j,i),&
	1409	rrtm_asdif(:,j,i), rrtm_aldir(:,j,i), rrtm_aldif(:,j,i), zenith,&
	1410	0.0_wp , day , solar_constant, rrtm_inflgsw,&
	1411	rrtm_iceflgsw , rrtm_liqflgsw, rrtm_cldfr , rrtm_sw_taucld ,&
	1412	rrtm_sw_ssacld , rrtm_sw_asmcld, rrtm_sw_fsfcld, rrtm_cicewp ,&
	1413	rrtm_cliqwp , rrtm_reice , rrtm_reliq , rrtm_sw_tauaer ,&
	1414	rrtm_sw_ssaaer , rrtm_sw_asmaer , rrtm_sw_ecaer , &
	1415	rrtm_swuflx , rrtm_swdflx , rrtm_swhr , &
	1416	rrtm_swuflxc , rrtm_swdflxc , rrtm_swhrc )
	1417
[1691]	1418	!
	1419	!-- Save fluxes
[1585]	1420	DO k = nzb, nzt+1
	1421	rad_sw_in(k,j,i) = rrtm_swdflx(0,k)
	1422	rad_sw_out(k,j,i) = rrtm_swuflx(0,k)
	1423	ENDDO
[1691]	1424
	1425	!
	1426	!-- Save heating rates (convert from K/d to K/s)
	1427	DO k = nzb+1, nzt+1
	1428	rad_sw_hr(k,j,i) = rrtm_swhr(0,k) * d_hours_day
	1429	rad_sw_cs_hr(k,j,i) = rrtm_swhrc(0,k) * d_hours_day
	1430	ENDDO
	1431
[1585]	1432	ENDIF
	1433
	1434	!
	1435	!-- Calculate surface net radiation
	1436	rad_net(j,i) = rad_sw_in(nzb,j,i) - rad_sw_out(nzb,j,i) &
	1437	+ rad_lw_in(nzb,j,i) - rad_lw_out(nzb,j,i)
	1438
	1439	ENDDO
	1440	ENDDO
	1441
	1442	CALL exchange_horiz( rad_lw_in, nbgp )
	1443	CALL exchange_horiz( rad_lw_out, nbgp )
[1691]	1444	CALL exchange_horiz( rad_lw_hr, nbgp )
	1445	CALL exchange_horiz( rad_lw_cs_hr, nbgp )
	1446
[1585]	1447	CALL exchange_horiz( rad_sw_in, nbgp )
	1448	CALL exchange_horiz( rad_sw_out, nbgp )
[1691]	1449	CALL exchange_horiz( rad_sw_hr, nbgp )
	1450	CALL exchange_horiz( rad_sw_cs_hr, nbgp )
	1451
[1585]	1452	CALL exchange_horiz_2d( rad_net, nbgp )
[1709]	1453	CALL exchange_horiz_2d( rad_lw_out_change_0, nbgp )
[1585]	1454	#endif
	1455
	1456	END SUBROUTINE radiation_rrtmg
	1457
	1458
	1459	!------------------------------------------------------------------------------!
	1460	! Description:
	1461	! ------------
[1682]	1462	!> Calculate the cosine of the zenith angle (variable is called zenith)
[1585]	1463	!------------------------------------------------------------------------------!
	1464	SUBROUTINE calc_zenith
	1465
	1466	IMPLICIT NONE
	1467
[1682]	1468	REAL(wp) :: declination, & !< solar declination angle
	1469	hour_angle !< solar hour angle
[1585]	1470	!
[1496]	1471	!-- Calculate current day and time based on the initial values and simulation
	1472	!-- time
[1585]	1473	day = day_init + INT(FLOOR( (time_utc_init + time_since_reference_point) &
	1474	/ 86400.0_wp ), KIND=iwp)
[1496]	1475	time_utc = MOD((time_utc_init + time_since_reference_point), 86400.0_wp)
	1476
	1477
	1478	!
	1479	!-- Calculate solar declination and hour angle
[1585]	1480	declination = ASIN( decl_1 * SIN(decl_2 * REAL(day, KIND=wp) - decl_3) )
[1496]	1481	hour_angle = 2.0_wp * pi * (time_utc / 86400.0_wp) + lon - pi
	1482
	1483	!
	1484	!-- Calculate zenith angle
[1585]	1485	zenith(0) = SIN(lat) * SIN(declination) + COS(lat) * COS(declination) &
[1496]	1486	* COS(hour_angle)
[1585]	1487	zenith(0) = MAX(0.0_wp,zenith(0))
[1496]	1488
	1489	!
[1585]	1490	!-- Check if the sun is up (otheriwse shortwave calculations can be skipped)
[1691]	1491	IF ( zenith(0) > 0.0_wp ) THEN
[1585]	1492	sun_up = .TRUE.
	1493	ELSE
	1494	sun_up = .FALSE.
	1495	END IF
[1496]	1496
[1585]	1497	END SUBROUTINE calc_zenith
	1498
[1606]	1499	#if defined ( __rrtmg ) && defined ( __netcdf )
[1585]	1500	!------------------------------------------------------------------------------!
	1501	! Description:
	1502	! ------------
[1682]	1503	!> Calculates surface albedo components based on Briegleb (1992) and
	1504	!> Briegleb et al. (1986)
[1585]	1505	!------------------------------------------------------------------------------!
	1506	SUBROUTINE calc_albedo
	1507
	1508	IMPLICIT NONE
	1509
	1510	IF ( sun_up ) THEN
[1496]	1511	!
[1585]	1512	!-- Ocean
	1513	IF ( albedo_type == 1 ) THEN
	1514	rrtm_aldir(0,:,:) = 0.026_wp / ( zenith(0)**1.7_wp + 0.065_wp ) &
	1515	+ 0.15_wp * ( zenith(0) - 0.1_wp ) &
	1516	* ( zenith(0) - 0.5_wp ) &
	1517	* ( zenith(0) - 1.0_wp )
	1518	rrtm_asdir(0,:,:) = rrtm_aldir(0,:,:)
	1519	!
	1520	!-- Snow
	1521	ELSEIF ( albedo_type == 16 ) THEN
	1522	IF ( zenith(0) < 0.5_wp ) THEN
	1523	rrtm_aldir(0,:,:) = 0.5_wp * (1.0_wp - aldif) &
	1524	* ( 3.0_wp / (1.0_wp + 4.0_wp &
	1525	* zenith(0))) - 1.0_wp
	1526	rrtm_asdir(0,:,:) = 0.5_wp * (1.0_wp - asdif) &
	1527	* ( 3.0_wp / (1.0_wp + 4.0_wp &
	1528	* zenith(0))) - 1.0_wp
[1496]	1529
[1585]	1530	rrtm_aldir(0,:,:) = MIN(0.98_wp, rrtm_aldir(0,:,:))
	1531	rrtm_asdir(0,:,:) = MIN(0.98_wp, rrtm_asdir(0,:,:))
	1532	ELSE
	1533	rrtm_aldir(0,:,:) = aldif
	1534	rrtm_asdir(0,:,:) = asdif
	1535	ENDIF
[1496]	1536	!
[1585]	1537	!-- Sea ice
	1538	ELSEIF ( albedo_type == 15 ) THEN
	1539	rrtm_aldir(0,:,:) = aldif
	1540	rrtm_asdir(0,:,:) = asdif
[1788]	1541
[1585]	1542	!
[1788]	1543	!-- Asphalt
	1544	ELSEIF ( albedo_type == 17 ) THEN
	1545	rrtm_aldir(0,:,:) = aldif
	1546	rrtm_asdir(0,:,:) = asdif
	1547	!
[1585]	1548	!-- Land surfaces
	1549	ELSE
	1550	SELECT CASE ( albedo_type )
[1496]	1551
[1585]	1552	!
	1553	!-- Surface types with strong zenith dependence
	1554	CASE ( 1, 2, 3, 4, 11, 12, 13 )
	1555	rrtm_aldir(0,:,:) = aldif * 1.4_wp / &
	1556	(1.0_wp + 0.8_wp * zenith(0))
	1557	rrtm_asdir(0,:,:) = asdif * 1.4_wp / &
	1558	(1.0_wp + 0.8_wp * zenith(0))
	1559	!
	1560	!-- Surface types with weak zenith dependence
	1561	CASE ( 5, 6, 7, 8, 9, 10, 14 )
	1562	rrtm_aldir(0,:,:) = aldif * 1.1_wp / &
	1563	(1.0_wp + 0.2_wp * zenith(0))
	1564	rrtm_asdir(0,:,:) = asdif * 1.1_wp / &
	1565	(1.0_wp + 0.2_wp * zenith(0))
[1496]	1566
[1585]	1567	CASE DEFAULT
	1568
	1569	END SELECT
	1570	ENDIF
	1571	!
	1572	!-- Diffusive albedo is taken from Table 2
	1573	rrtm_aldif(0,:,:) = aldif
	1574	rrtm_asdif(0,:,:) = asdif
	1575
	1576	ELSE
	1577
	1578	rrtm_aldir(0,:,:) = 0.0_wp
	1579	rrtm_asdir(0,:,:) = 0.0_wp
	1580	rrtm_aldif(0,:,:) = 0.0_wp
	1581	rrtm_asdif(0,:,:) = 0.0_wp
	1582	ENDIF
	1583	END SUBROUTINE calc_albedo
	1584
	1585	!------------------------------------------------------------------------------!
	1586	! Description:
	1587	! ------------
[1682]	1588	!> Read sounding data (pressure and temperature) from RADIATION_DATA.
[1585]	1589	!------------------------------------------------------------------------------!
	1590	SUBROUTINE read_sounding_data
	1591
	1592	IMPLICIT NONE
	1593
[1691]	1594	INTEGER(iwp) :: id, & !< NetCDF id of input file
	1595	id_dim_zrad, & !< pressure level id in the NetCDF file
	1596	id_var, & !< NetCDF variable id
	1597	k, & !< loop index
	1598	nz_snd, & !< number of vertical levels in the sounding data
	1599	nz_snd_start, & !< start vertical index for sounding data to be used
	1600	nz_snd_end !< end vertical index for souding data to be used
[1585]	1601
[1691]	1602	REAL(wp) :: t_surface !< actual surface temperature
[1585]	1603
[1691]	1604	REAL(wp), DIMENSION(:), ALLOCATABLE :: hyp_snd_tmp, & !< temporary hydrostatic pressure profile (sounding)
	1605	t_snd_tmp !< temporary temperature profile (sounding)
[1585]	1606
	1607	!
	1608	!-- In case of updates, deallocate arrays first (sufficient to check one
	1609	!-- array as the others are automatically allocated). This is required
	1610	!-- because nzt_rad might change during the update
	1611	IF ( ALLOCATED ( hyp_snd ) ) THEN
	1612	DEALLOCATE( hyp_snd )
	1613	DEALLOCATE( t_snd )
	1614	DEALLOCATE( q_snd )
	1615	DEALLOCATE ( rrtm_play )
	1616	DEALLOCATE ( rrtm_plev )
	1617	DEALLOCATE ( rrtm_tlay )
	1618	DEALLOCATE ( rrtm_tlev )
[1691]	1619
[1585]	1620	DEALLOCATE ( rrtm_h2ovmr )
	1621	DEALLOCATE ( rrtm_cicewp )
	1622	DEALLOCATE ( rrtm_cldfr )
	1623	DEALLOCATE ( rrtm_cliqwp )
	1624	DEALLOCATE ( rrtm_reice )
	1625	DEALLOCATE ( rrtm_reliq )
	1626	DEALLOCATE ( rrtm_lw_taucld )
	1627	DEALLOCATE ( rrtm_lw_tauaer )
[1691]	1628
[1585]	1629	DEALLOCATE ( rrtm_lwdflx )
[1691]	1630	DEALLOCATE ( rrtm_lwdflxc )
[1585]	1631	DEALLOCATE ( rrtm_lwuflx )
[1691]	1632	DEALLOCATE ( rrtm_lwuflxc )
	1633	DEALLOCATE ( rrtm_lwuflx_dt )
	1634	DEALLOCATE ( rrtm_lwuflxc_dt )
[1585]	1635	DEALLOCATE ( rrtm_lwhr )
	1636	DEALLOCATE ( rrtm_lwhrc )
[1691]	1637
[1585]	1638	DEALLOCATE ( rrtm_sw_taucld )
	1639	DEALLOCATE ( rrtm_sw_ssacld )
	1640	DEALLOCATE ( rrtm_sw_asmcld )
	1641	DEALLOCATE ( rrtm_sw_fsfcld )
	1642	DEALLOCATE ( rrtm_sw_tauaer )
	1643	DEALLOCATE ( rrtm_sw_ssaaer )
	1644	DEALLOCATE ( rrtm_sw_asmaer )
[1691]	1645	DEALLOCATE ( rrtm_sw_ecaer )
	1646
[1585]	1647	DEALLOCATE ( rrtm_swdflx )
[1691]	1648	DEALLOCATE ( rrtm_swdflxc )
[1585]	1649	DEALLOCATE ( rrtm_swuflx )
[1691]	1650	DEALLOCATE ( rrtm_swuflxc )
[1585]	1651	DEALLOCATE ( rrtm_swhr )
	1652	DEALLOCATE ( rrtm_swhrc )
[1691]	1653
[1585]	1654	ENDIF
	1655
	1656	!
	1657	!-- Open file for reading
	1658	nc_stat = NF90_OPEN( rrtm_input_file, NF90_NOWRITE, id )
[1783]	1659	CALL netcdf_handle_error_rad( 'read_sounding_data', 549 )
[1585]	1660
	1661	!
	1662	!-- Inquire dimension of z axis and save in nz_snd
	1663	nc_stat = NF90_INQ_DIMID( id, "Pressure", id_dim_zrad )
	1664	nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim_zrad, len = nz_snd )
[1783]	1665	CALL netcdf_handle_error_rad( 'read_sounding_data', 551 )
[1585]	1666
	1667	!
	1668	! !-- Allocate temporary array for storing pressure data
[1701]	1669	ALLOCATE( hyp_snd_tmp(1:nz_snd) )
[1585]	1670	hyp_snd_tmp = 0.0_wp
	1671
	1672
	1673	!-- Read pressure from file
	1674	nc_stat = NF90_INQ_VARID( id, "Pressure", id_var )
[1691]	1675	nc_stat = NF90_GET_VAR( id, id_var, hyp_snd_tmp(:), start = (/1/), &
[1585]	1676	count = (/nz_snd/) )
[1783]	1677	CALL netcdf_handle_error_rad( 'read_sounding_data', 552 )
[1585]	1678
	1679	!
	1680	!-- Allocate temporary array for storing temperature data
[1701]	1681	ALLOCATE( t_snd_tmp(1:nz_snd) )
[1585]	1682	t_snd_tmp = 0.0_wp
	1683
	1684	!
	1685	!-- Read temperature from file
	1686	nc_stat = NF90_INQ_VARID( id, "ReferenceTemperature", id_var )
[1691]	1687	nc_stat = NF90_GET_VAR( id, id_var, t_snd_tmp(:), start = (/1/), &
[1585]	1688	count = (/nz_snd/) )
[1783]	1689	CALL netcdf_handle_error_rad( 'read_sounding_data', 553 )
[1585]	1690
	1691	!
	1692	!-- Calculate start of sounding data
	1693	nz_snd_start = nz_snd + 1
[1701]	1694	nz_snd_end = nz_snd + 1
[1585]	1695
	1696	!
	1697	!-- Start filling vertical dimension at 10hPa above the model domain (hyp is
	1698	!-- in Pa, hyp_snd in hPa).
	1699	DO k = 1, nz_snd
[1691]	1700	IF ( hyp_snd_tmp(k) < ( hyp(nzt+1) - 1000.0_wp) * 0.01_wp ) THEN
[1585]	1701	nz_snd_start = k
	1702	EXIT
	1703	END IF
	1704	END DO
	1705
[1691]	1706	IF ( nz_snd_start <= nz_snd ) THEN
[1701]	1707	nz_snd_end = nz_snd
[1585]	1708	END IF
	1709
	1710
	1711	!
	1712	!-- Calculate of total grid points for RRTMG calculations
[1701]	1713	nzt_rad = nzt + nz_snd_end - nz_snd_start + 1
[1585]	1714
	1715	!
	1716	!-- Save data above LES domain in hyp_snd, t_snd and q_snd
	1717	!-- Note: q_snd_tmp is not calculated at the moment (dry residual atmosphere)
	1718	ALLOCATE( hyp_snd(nzb+1:nzt_rad) )
	1719	ALLOCATE( t_snd(nzb+1:nzt_rad) )
	1720	ALLOCATE( q_snd(nzb+1:nzt_rad) )
	1721	hyp_snd = 0.0_wp
	1722	t_snd = 0.0_wp
	1723	q_snd = 0.0_wp
	1724
[1757]	1725	hyp_snd(nzt+2:nzt_rad) = hyp_snd_tmp(nz_snd_start+1:nz_snd_end)
	1726	t_snd(nzt+2:nzt_rad) = t_snd_tmp(nz_snd_start+1:nz_snd_end)
[1585]	1727
	1728	nc_stat = NF90_CLOSE( id )
	1729
	1730	!
	1731	!-- Calculate pressure levels on zu and zw grid. Sounding data is added at
	1732	!-- top of the LES domain. This routine does not consider horizontal or
	1733	!-- vertical variability of pressure and temperature
	1734	ALLOCATE ( rrtm_play(0:0,nzb+1:nzt_rad+1) )
	1735	ALLOCATE ( rrtm_plev(0:0,nzb+1:nzt_rad+2) )
	1736
[1691]	1737	t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
[1585]	1738	DO k = nzb+1, nzt+1
	1739	rrtm_play(0,k) = hyp(k) * 0.01_wp
	1740	rrtm_plev(0,k) = surface_pressure * ( (t_surface - g/cp * zw(k-1)) / &
	1741	t_surface )**(1.0_wp/0.286_wp)
	1742	ENDDO
	1743
	1744	DO k = nzt+2, nzt_rad
	1745	rrtm_play(0,k) = hyp_snd(k)
	1746	rrtm_plev(0,k) = 0.5_wp * ( rrtm_play(0,k) + rrtm_play(0,k-1) )
	1747	ENDDO
	1748	rrtm_plev(0,nzt_rad+1) = MAX( 0.5 * hyp_snd(nzt_rad), &
	1749	1.5 * hyp_snd(nzt_rad) &
	1750	- 0.5 * hyp_snd(nzt_rad-1) )
	1751	rrtm_plev(0,nzt_rad+2) = MIN( 1.0E-4_wp, &
	1752	0.25_wp * rrtm_plev(0,nzt_rad+1) )
	1753
	1754	rrtm_play(0,nzt_rad+1) = 0.5 * rrtm_plev(0,nzt_rad+1)
	1755
	1756	!
	1757	!-- Calculate temperature/humidity levels at top of the LES domain.
	1758	!-- Currently, the temperature is taken from sounding data (might lead to a
	1759	!-- temperature jump at interface. To do: Humidity is currently not
	1760	!-- calculated above the LES domain.
	1761	ALLOCATE ( rrtm_tlay(0:0,nzb+1:nzt_rad+1) )
	1762	ALLOCATE ( rrtm_tlev(0:0,nzb+1:nzt_rad+2) )
	1763	ALLOCATE ( rrtm_h2ovmr(0:0,nzb+1:nzt_rad+1) )
	1764
	1765	DO k = nzt+8, nzt_rad
	1766	rrtm_tlay(0,k) = t_snd(k)
	1767	rrtm_h2ovmr(0,k) = q_snd(k)
	1768	ENDDO
[1691]	1769	rrtm_tlay(0,nzt_rad+1) = 2.0_wp * rrtm_tlay(0,nzt_rad) &
	1770	- rrtm_tlay(0,nzt_rad-1)
[1585]	1771	DO k = nzt+9, nzt_rad+1
	1772	rrtm_tlev(0,k) = rrtm_tlay(0,k-1) + (rrtm_tlay(0,k) &
	1773	- rrtm_tlay(0,k-1)) &
	1774	/ ( rrtm_play(0,k) - rrtm_play(0,k-1) ) &
	1775	* ( rrtm_plev(0,k) - rrtm_play(0,k-1) )
	1776	ENDDO
	1777	rrtm_h2ovmr(0,nzt_rad+1) = rrtm_h2ovmr(0,nzt_rad)
	1778
	1779	rrtm_tlev(0,nzt_rad+2) = 2.0_wp * rrtm_tlay(0,nzt_rad+1) &
	1780	- rrtm_tlev(0,nzt_rad)
	1781	!
	1782	!-- Allocate remaining RRTMG arrays
	1783	ALLOCATE ( rrtm_cicewp(0:0,nzb+1:nzt_rad+1) )
	1784	ALLOCATE ( rrtm_cldfr(0:0,nzb+1:nzt_rad+1) )
	1785	ALLOCATE ( rrtm_cliqwp(0:0,nzb+1:nzt_rad+1) )
	1786	ALLOCATE ( rrtm_reice(0:0,nzb+1:nzt_rad+1) )
	1787	ALLOCATE ( rrtm_reliq(0:0,nzb+1:nzt_rad+1) )
	1788	ALLOCATE ( rrtm_lw_taucld(1:nbndlw+1,0:0,nzb+1:nzt_rad+1) )
	1789	ALLOCATE ( rrtm_lw_tauaer(0:0,nzb+1:nzt_rad+1,1:nbndlw+1) )
	1790	ALLOCATE ( rrtm_sw_taucld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
	1791	ALLOCATE ( rrtm_sw_ssacld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
	1792	ALLOCATE ( rrtm_sw_asmcld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
	1793	ALLOCATE ( rrtm_sw_fsfcld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
	1794	ALLOCATE ( rrtm_sw_tauaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) )
	1795	ALLOCATE ( rrtm_sw_ssaaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) )
	1796	ALLOCATE ( rrtm_sw_asmaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) )
	1797	ALLOCATE ( rrtm_sw_ecaer(0:0,nzb+1:nzt_rad+1,1:naerec+1) )
	1798
	1799	!
	1800	!-- The ice phase is currently not considered in PALM
	1801	rrtm_cicewp = 0.0_wp
	1802	rrtm_reice = 0.0_wp
	1803
	1804	!
	1805	!-- Set other parameters (move to NAMELIST parameters in the future)
	1806	rrtm_lw_tauaer = 0.0_wp
	1807	rrtm_lw_taucld = 0.0_wp
	1808	rrtm_sw_taucld = 0.0_wp
	1809	rrtm_sw_ssacld = 0.0_wp
	1810	rrtm_sw_asmcld = 0.0_wp
	1811	rrtm_sw_fsfcld = 0.0_wp
	1812	rrtm_sw_tauaer = 0.0_wp
	1813	rrtm_sw_ssaaer = 0.0_wp
	1814	rrtm_sw_asmaer = 0.0_wp
	1815	rrtm_sw_ecaer = 0.0_wp
	1816
	1817
	1818	ALLOCATE ( rrtm_swdflx(0:0,nzb:nzt_rad+1) )
	1819	ALLOCATE ( rrtm_swuflx(0:0,nzb:nzt_rad+1) )
	1820	ALLOCATE ( rrtm_swhr(0:0,nzb+1:nzt_rad+1) )
	1821	ALLOCATE ( rrtm_swuflxc(0:0,nzb:nzt_rad+1) )
	1822	ALLOCATE ( rrtm_swdflxc(0:0,nzb:nzt_rad+1) )
	1823	ALLOCATE ( rrtm_swhrc(0:0,nzb+1:nzt_rad+1) )
	1824
	1825	rrtm_swdflx = 0.0_wp
	1826	rrtm_swuflx = 0.0_wp
	1827	rrtm_swhr = 0.0_wp
	1828	rrtm_swuflxc = 0.0_wp
	1829	rrtm_swdflxc = 0.0_wp
	1830	rrtm_swhrc = 0.0_wp
	1831
	1832	ALLOCATE ( rrtm_lwdflx(0:0,nzb:nzt_rad+1) )
	1833	ALLOCATE ( rrtm_lwuflx(0:0,nzb:nzt_rad+1) )
	1834	ALLOCATE ( rrtm_lwhr(0:0,nzb+1:nzt_rad+1) )
	1835	ALLOCATE ( rrtm_lwuflxc(0:0,nzb:nzt_rad+1) )
	1836	ALLOCATE ( rrtm_lwdflxc(0:0,nzb:nzt_rad+1) )
	1837	ALLOCATE ( rrtm_lwhrc(0:0,nzb+1:nzt_rad+1) )
	1838
	1839	rrtm_lwdflx = 0.0_wp
	1840	rrtm_lwuflx = 0.0_wp
	1841	rrtm_lwhr = 0.0_wp
	1842	rrtm_lwuflxc = 0.0_wp
	1843	rrtm_lwdflxc = 0.0_wp
	1844	rrtm_lwhrc = 0.0_wp
	1845
[1691]	1846	ALLOCATE ( rrtm_lwuflx_dt(0:0,nzb:nzt_rad+1) )
	1847	ALLOCATE ( rrtm_lwuflxc_dt(0:0,nzb:nzt_rad+1) )
[1585]	1848
[1709]	1849	rrtm_lwuflx_dt = 0.0_wp
[1691]	1850	rrtm_lwuflxc_dt = 0.0_wp
	1851
[1585]	1852	END SUBROUTINE read_sounding_data
	1853
	1854
	1855	!------------------------------------------------------------------------------!
	1856	! Description:
	1857	! ------------
[1682]	1858	!> Read trace gas data from file
[1585]	1859	!------------------------------------------------------------------------------!
	1860	SUBROUTINE read_trace_gas_data
	1861
	1862	USE rrsw_ncpar
	1863
	1864	IMPLICIT NONE
	1865
[1691]	1866	INTEGER(iwp), PARAMETER :: num_trace_gases = 9 !< number of trace gases (absorbers)
[1585]	1867
[1691]	1868	CHARACTER(LEN=5), DIMENSION(num_trace_gases), PARAMETER :: & !< trace gas names
[1585]	1869	trace_names = (/'O3 ', 'CO2 ', 'CH4 ', 'N2O ', 'O2 ', &
	1870	'CFC11', 'CFC12', 'CFC22', 'CCL4 '/)
	1871
[1691]	1872	INTEGER(iwp) :: id, & !< NetCDF id
	1873	k, & !< loop index
	1874	m, & !< loop index
	1875	n, & !< loop index
	1876	nabs, & !< number of absorbers
	1877	np, & !< number of pressure levels
	1878	id_abs, & !< NetCDF id of the respective absorber
	1879	id_dim, & !< NetCDF id of asborber's dimension
	1880	id_var !< NetCDf id ot the absorber
[1585]	1881
	1882	REAL(wp) :: p_mls_l, p_mls_u, p_wgt_l, p_wgt_u, p_mls_m
	1883
	1884
[1682]	1885	REAL(wp), DIMENSION(:), ALLOCATABLE :: p_mls, & !< pressure levels for the absorbers
	1886	rrtm_play_tmp, & !< temporary array for pressure zu-levels
	1887	rrtm_plev_tmp, & !< temporary array for pressure zw-levels
	1888	trace_path_tmp !< temporary array for storing trace gas path data
[1585]	1889
[1682]	1890	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: trace_mls, & !< array for storing the absorber amounts
	1891	trace_mls_path, & !< array for storing trace gas path data
	1892	trace_mls_tmp !< temporary array for storing trace gas data
[1585]	1893
	1894
	1895	!
	1896	!-- In case of updates, deallocate arrays first (sufficient to check one
	1897	!-- array as the others are automatically allocated)
	1898	IF ( ALLOCATED ( rrtm_o3vmr ) ) THEN
	1899	DEALLOCATE ( rrtm_o3vmr )
	1900	DEALLOCATE ( rrtm_co2vmr )
	1901	DEALLOCATE ( rrtm_ch4vmr )
	1902	DEALLOCATE ( rrtm_n2ovmr )
	1903	DEALLOCATE ( rrtm_o2vmr )
	1904	DEALLOCATE ( rrtm_cfc11vmr )
	1905	DEALLOCATE ( rrtm_cfc12vmr )
	1906	DEALLOCATE ( rrtm_cfc22vmr )
	1907	DEALLOCATE ( rrtm_ccl4vmr )
	1908	ENDIF
	1909
	1910	!
	1911	!-- Allocate trace gas profiles
	1912	ALLOCATE ( rrtm_o3vmr(0:0,1:nzt_rad+1) )
	1913	ALLOCATE ( rrtm_co2vmr(0:0,1:nzt_rad+1) )
	1914	ALLOCATE ( rrtm_ch4vmr(0:0,1:nzt_rad+1) )
	1915	ALLOCATE ( rrtm_n2ovmr(0:0,1:nzt_rad+1) )
	1916	ALLOCATE ( rrtm_o2vmr(0:0,1:nzt_rad+1) )
	1917	ALLOCATE ( rrtm_cfc11vmr(0:0,1:nzt_rad+1) )
	1918	ALLOCATE ( rrtm_cfc12vmr(0:0,1:nzt_rad+1) )
	1919	ALLOCATE ( rrtm_cfc22vmr(0:0,1:nzt_rad+1) )
	1920	ALLOCATE ( rrtm_ccl4vmr(0:0,1:nzt_rad+1) )
	1921
	1922	!
	1923	!-- Open file for reading
	1924	nc_stat = NF90_OPEN( rrtm_input_file, NF90_NOWRITE, id )
[1783]	1925	CALL netcdf_handle_error_rad( 'read_trace_gas_data', 549 )
[1585]	1926	!
	1927	!-- Inquire dimension ids and dimensions
	1928	nc_stat = NF90_INQ_DIMID( id, "Pressure", id_dim )
[1783]	1929	CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
[1585]	1930	nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, len = np)
[1783]	1931	CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
[1585]	1932
	1933	nc_stat = NF90_INQ_DIMID( id, "Absorber", id_dim )
[1783]	1934	CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
[1585]	1935	nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, len = nabs )
[1783]	1936	CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
[1585]	1937
	1938
	1939	!
	1940	!-- Allocate pressure, and trace gas arrays
	1941	ALLOCATE( p_mls(1:np) )
	1942	ALLOCATE( trace_mls(1:num_trace_gases,1:np) )
	1943	ALLOCATE( trace_mls_tmp(1:nabs,1:np) )
	1944
	1945
	1946	nc_stat = NF90_INQ_VARID( id, "Pressure", id_var )
[1783]	1947	CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
[1585]	1948	nc_stat = NF90_GET_VAR( id, id_var, p_mls )
[1783]	1949	CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
[1585]	1950
	1951	nc_stat = NF90_INQ_VARID( id, "AbsorberAmountMLS", id_var )
[1783]	1952	CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
[1585]	1953	nc_stat = NF90_GET_VAR( id, id_var, trace_mls_tmp )
[1783]	1954	CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
[1585]	1955
	1956
	1957	!
	1958	!-- Write absorber amounts (mls) to trace_mls
	1959	DO n = 1, num_trace_gases
	1960	CALL getAbsorberIndex( TRIM( trace_names(n) ), id_abs )
	1961
	1962	trace_mls(n,1:np) = trace_mls_tmp(id_abs,1:np)
	1963
	1964	!
	1965	!-- Replace missing values by zero
	1966	WHERE ( trace_mls(n,:) > 2.0_wp )
	1967	trace_mls(n,:) = 0.0_wp
	1968	END WHERE
	1969	END DO
	1970
	1971	DEALLOCATE ( trace_mls_tmp )
	1972
	1973	nc_stat = NF90_CLOSE( id )
[1783]	1974	CALL netcdf_handle_error_rad( 'read_trace_gas_data', 551 )
[1585]	1975
	1976	!
	1977	!-- Add extra pressure level for calculations of the trace gas paths
	1978	ALLOCATE ( rrtm_play_tmp(1:nzt_rad+1) )
	1979	ALLOCATE ( rrtm_plev_tmp(1:nzt_rad+2) )
	1980
	1981	rrtm_play_tmp(1:nzt_rad) = rrtm_play(0,1:nzt_rad)
	1982	rrtm_plev_tmp(1:nzt_rad+1) = rrtm_plev(0,1:nzt_rad+1)
	1983	rrtm_play_tmp(nzt_rad+1) = rrtm_plev(0,nzt_rad+1) * 0.5_wp
	1984	rrtm_plev_tmp(nzt_rad+2) = MIN( 1.0E-4_wp, 0.25_wp &
	1985	* rrtm_plev(0,nzt_rad+1) )
	1986
	1987	!
	1988	!-- Calculate trace gas path (zero at surface) with interpolation to the
	1989	!-- sounding levels
	1990	ALLOCATE ( trace_mls_path(1:nzt_rad+2,1:num_trace_gases) )
	1991
	1992	trace_mls_path(nzb+1,:) = 0.0_wp
	1993
	1994	DO k = nzb+2, nzt_rad+2
	1995	DO m = 1, num_trace_gases
	1996	trace_mls_path(k,m) = trace_mls_path(k-1,m)
	1997
	1998	!
	1999	!-- When the pressure level is higher than the trace gas pressure
	2000	!-- level, assume that
[1691]	2001	IF ( rrtm_plev_tmp(k-1) > p_mls(1) ) THEN
[1585]	2002
	2003	trace_mls_path(k,m) = trace_mls_path(k,m) + trace_mls(m,1) &
	2004	* ( rrtm_plev_tmp(k-1) &
	2005	- MAX( p_mls(1), rrtm_plev_tmp(k) ) &
	2006	) / g
	2007	ENDIF
	2008
	2009	!
	2010	!-- Integrate for each sounding level from the contributing p_mls
	2011	!-- levels
	2012	DO n = 2, np
	2013	!
	2014	!-- Limit p_mls so that it is within the model level
	2015	p_mls_u = MIN( rrtm_plev_tmp(k-1), &
	2016	MAX( rrtm_plev_tmp(k), p_mls(n) ) )
	2017	p_mls_l = MIN( rrtm_plev_tmp(k-1), &
	2018	MAX( rrtm_plev_tmp(k), p_mls(n-1) ) )
	2019
[1691]	2020	IF ( p_mls_l > p_mls_u ) THEN
[1585]	2021
	2022	!
	2023	!-- Calculate weights for interpolation
	2024	p_mls_m = 0.5_wp * (p_mls_l + p_mls_u)
	2025	p_wgt_u = (p_mls(n-1) - p_mls_m) / (p_mls(n-1) - p_mls(n))
	2026	p_wgt_l = (p_mls_m - p_mls(n)) / (p_mls(n-1) - p_mls(n))
	2027
	2028	!
	2029	!-- Add level to trace gas path
	2030	trace_mls_path(k,m) = trace_mls_path(k,m) &
	2031	+ ( p_wgt_u * trace_mls(m,n) &
	2032	+ p_wgt_l * trace_mls(m,n-1) ) &
[1691]	2033	* (p_mls_l - p_mls_u) / g
[1585]	2034	ENDIF
	2035	ENDDO
	2036
[1691]	2037	IF ( rrtm_plev_tmp(k) < p_mls(np) ) THEN
[1585]	2038	trace_mls_path(k,m) = trace_mls_path(k,m) + trace_mls(m,np) &
	2039	* ( MIN( rrtm_plev_tmp(k-1), p_mls(np) ) &
	2040	- rrtm_plev_tmp(k) &
	2041	) / g
	2042	ENDIF
[1496]	2043	ENDDO
	2044	ENDDO
	2045
	2046
[1585]	2047	!
	2048	!-- Prepare trace gas path profiles
	2049	ALLOCATE ( trace_path_tmp(1:nzt_rad+1) )
[1496]	2050
[1585]	2051	DO m = 1, num_trace_gases
	2052
	2053	trace_path_tmp(1:nzt_rad+1) = ( trace_mls_path(2:nzt_rad+2,m) &
	2054	- trace_mls_path(1:nzt_rad+1,m) ) * g &
	2055	/ ( rrtm_plev_tmp(1:nzt_rad+1) &
	2056	- rrtm_plev_tmp(2:nzt_rad+2) )
	2057
	2058	!
	2059	!-- Save trace gas paths to the respective arrays
	2060	SELECT CASE ( TRIM( trace_names(m) ) )
	2061
	2062	CASE ( 'O3' )
	2063
	2064	rrtm_o3vmr(0,:) = trace_path_tmp(:)
	2065
	2066	CASE ( 'CO2' )
	2067
	2068	rrtm_co2vmr(0,:) = trace_path_tmp(:)
	2069
	2070	CASE ( 'CH4' )
	2071
	2072	rrtm_ch4vmr(0,:) = trace_path_tmp(:)
	2073
	2074	CASE ( 'N2O' )
	2075
	2076	rrtm_n2ovmr(0,:) = trace_path_tmp(:)
	2077
	2078	CASE ( 'O2' )
	2079
	2080	rrtm_o2vmr(0,:) = trace_path_tmp(:)
	2081
	2082	CASE ( 'CFC11' )
	2083
	2084	rrtm_cfc11vmr(0,:) = trace_path_tmp(:)
	2085
	2086	CASE ( 'CFC12' )
	2087
	2088	rrtm_cfc12vmr(0,:) = trace_path_tmp(:)
	2089
	2090	CASE ( 'CFC22' )
	2091
	2092	rrtm_cfc22vmr(0,:) = trace_path_tmp(:)
	2093
	2094	CASE ( 'CCL4' )
	2095
	2096	rrtm_ccl4vmr(0,:) = trace_path_tmp(:)
	2097
	2098	CASE DEFAULT
	2099
	2100	END SELECT
	2101
	2102	ENDDO
	2103
	2104	DEALLOCATE ( trace_path_tmp )
	2105	DEALLOCATE ( trace_mls_path )
	2106	DEALLOCATE ( rrtm_play_tmp )
	2107	DEALLOCATE ( rrtm_plev_tmp )
	2108	DEALLOCATE ( trace_mls )
	2109	DEALLOCATE ( p_mls )
	2110
	2111	END SUBROUTINE read_trace_gas_data
	2112
[1826]	2113
[1783]	2114	SUBROUTINE netcdf_handle_error_rad( routine_name, errno )
	2115
	2116	USE control_parameters, &
	2117	ONLY: message_string
	2118
	2119	USE NETCDF
	2120
	2121	USE pegrid
	2122
	2123	IMPLICIT NONE
	2124
	2125	CHARACTER(LEN=6) :: message_identifier
	2126	CHARACTER(LEN=*) :: routine_name
	2127
	2128	INTEGER(iwp) :: errno
	2129
	2130	IF ( nc_stat /= NF90_NOERR ) THEN
	2131
	2132	WRITE( message_identifier, '(''NC'',I4.4)' ) errno
	2133	message_string = TRIM( NF90_STRERROR( nc_stat ) )
	2134
	2135	CALL message( routine_name, message_identifier, 2, 2, 0, 6, 1 )
	2136
	2137	ENDIF
	2138
	2139	END SUBROUTINE netcdf_handle_error_rad
[1585]	2140	#endif
	2141
	2142
[1551]	2143	!------------------------------------------------------------------------------!
	2144	! Description:
	2145	! ------------
[1682]	2146	!> Calculate temperature tendency due to radiative cooling/heating.
	2147	!> Cache-optimized version.
[1551]	2148	!------------------------------------------------------------------------------!
[1585]	2149	SUBROUTINE radiation_tendency_ij ( i, j, tend )
[1496]	2150
[1585]	2151	USE cloud_parameters, &
[1691]	2152	ONLY: pt_d_t
[1551]	2153
[1585]	2154	IMPLICIT NONE
	2155
[1691]	2156	INTEGER(iwp) :: i, j, k !< loop indices
[1585]	2157
[1691]	2158	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend !< pt tendency term
[1585]	2159
	2160	#if defined ( __rrtmg )
	2161	!
[1691]	2162	!-- Calculate tendency based on heating rate
[1585]	2163	DO k = nzb+1, nzt+1
[1691]	2164	tend(k,j,i) = tend(k,j,i) + (rad_lw_hr(k,j,i) + rad_sw_hr(k,j,i)) &
	2165	* pt_d_t(k) * d_seconds_hour
[1585]	2166	ENDDO
	2167
	2168	#endif
	2169
	2170	END SUBROUTINE radiation_tendency_ij
	2171
	2172
[1551]	2173	!------------------------------------------------------------------------------!
	2174	! Description:
	2175	! ------------
[1682]	2176	!> Calculate temperature tendency due to radiative cooling/heating.
	2177	!> Vector-optimized version
[1551]	2178	!------------------------------------------------------------------------------!
[1585]	2179	SUBROUTINE radiation_tendency ( tend )
[1551]	2180
[1585]	2181	USE cloud_parameters, &
[1691]	2182	ONLY: pt_d_t
[1551]	2183
[1585]	2184	USE indices, &
	2185	ONLY: nxl, nxr, nyn, nys
	2186
	2187	IMPLICIT NONE
	2188
[1691]	2189	INTEGER(iwp) :: i, j, k !< loop indices
[1585]	2190
[1691]	2191	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend !< pt tendency term
[1585]	2192
	2193	#if defined ( __rrtmg )
[1691]	2194	!
	2195	!-- Calculate tendency based on heating rate
[1585]	2196	DO i = nxl, nxr
	2197	DO j = nys, nyn
	2198	DO k = nzb+1, nzt+1
[1691]	2199	tend(k,j,i) = tend(k,j,i) + ( rad_lw_hr(k,j,i) &
	2200	+ rad_sw_hr(k,j,i) ) * pt_d_t(k) &
	2201	* d_seconds_hour
[1585]	2202	ENDDO
	2203	ENDDO
	2204	ENDDO
	2205	#endif
	2206
	2207	END SUBROUTINE radiation_tendency
	2208
[1496]	2209	END MODULE radiation_model_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |