Home

Context Navigation

source: palm/trunk/SOURCE/radiation_model.f90 @ 1684

Last change on this file since 1684 was 1683, checked in by knoop, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 64.3 KB

Rev	Line
[1682]	1	!> @file radiation_model.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	!
[1585]	16	! Copyright 1997-2015 Leibniz Universitaet Hannover
[1496]	17	!--------------------------------------------------------------------------------!
	18	!
	19	! Current revisions:
	20	! -----------------
	21	!
[1683]	22	!
[1496]	23	! Former revisions:
	24	! -----------------
	25	! $Id: radiation_model.f90 1683 2015-10-07 23:57:51Z knoop $
	26	!
[1683]	27	! 1682 2015-10-07 23:56:08Z knoop
	28	! Code annotations made doxygen readable
	29	!
[1607]	30	! 1606 2015-06-29 10:43:37Z maronga
	31	! Added preprocessor directive __netcdf to allow for compiling without netCDF.
	32	! Note, however, that RRTMG cannot be used without netCDF.
	33	!
[1591]	34	! 1590 2015-05-08 13:56:27Z maronga
	35	! Bugfix: definition of character strings requires same length for all elements
	36	!
[1588]	37	! 1587 2015-05-04 14:19:01Z maronga
	38	! Added albedo class for snow
	39	!
[1586]	40	! 1585 2015-04-30 07:05:52Z maronga
	41	! Added support for RRTMG
	42	!
[1572]	43	! 1571 2015-03-12 16:12:49Z maronga
	44	! Added missing KIND attribute. Removed upper-case variable names
	45	!
[1552]	46	! 1551 2015-03-03 14:18:16Z maronga
	47	! Added support for data output. Various variables have been renamed. Added
	48	! interface for different radiation schemes (currently: clear-sky, constant, and
	49	! RRTM (not yet implemented).
	50	!
[1497]	51	! 1496 2014-12-02 17:25:50Z maronga
	52	! Initial revision
	53	!
[1496]	54	!
	55	! Description:
	56	! ------------
[1682]	57	!> Radiation models and interfaces
	58	!> @todo move variable definitions used in init_radiation only to the subroutine
	59	!> as they are no longer required after initialization.
	60	!> @todo Output of full column vertical profiles used in RRTMG
	61	!> @todo Output of other rrtm arrays (such as volume mixing ratios)
	62	!> @todo Adapt for use with topography
	63	!> @todo Remove 3D dummy arrays (such as clear-sky output)
	64	!>
	65	!> @note Many variables have a leading dummy dimension (0:0) in order to
	66	!> match the assume-size shape expected by the RRTMG model.
[1496]	67	!------------------------------------------------------------------------------!
[1682]	68	MODULE radiation_model_mod
	69
[1496]	70
	71	USE arrays_3d, &
[1585]	72	ONLY: hyp, pt, q, ql, zw
[1496]	73
[1585]	74	USE cloud_parameters, &
	75	ONLY: cp, l_v, nc_const, rho_l, sigma_gc
	76
	77	USE constants, &
	78	ONLY: pi
	79
[1496]	80	USE control_parameters, &
[1585]	81	ONLY: cloud_droplets, cloud_physics, g, initializing_actions, &
	82	large_scale_forcing, lsf_surf, phi, pt_surface, &
	83	surface_pressure, time_since_reference_point
[1496]	84
	85	USE indices, &
[1585]	86	ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb_s_inner, nzb, nzt
[1496]	87
	88	USE kinds
	89
[1606]	90	#if defined ( __netcdf )
[1585]	91	USE netcdf
[1606]	92	#endif
[1585]	93
[1551]	94	USE netcdf_control, &
	95	ONLY: dots_label, dots_num, dots_unit
[1496]	96
[1585]	97	#if defined ( __rrtmg )
	98	USE parrrsw, &
	99	ONLY: naerec, nbndsw
[1551]	100
[1585]	101	USE parrrtm, &
	102	ONLY: nbndlw
	103
	104	USE rrtmg_lw_init, &
	105	ONLY: rrtmg_lw_ini
	106
	107	USE rrtmg_sw_init, &
	108	ONLY: rrtmg_sw_ini
	109
	110	USE rrtmg_lw_rad, &
	111	ONLY: rrtmg_lw
	112
	113	USE rrtmg_sw_rad, &
	114	ONLY: rrtmg_sw
	115	#endif
	116
	117
	118
[1496]	119	IMPLICIT NONE
	120
[1585]	121	CHARACTER(10) :: radiation_scheme = 'clear-sky' ! 'constant', 'clear-sky', or 'rrtmg'
[1551]	122
[1585]	123	!
	124	!-- Predefined Land surface classes (albedo_type) after Briegleb (1992)
[1590]	125	CHARACTER(37), DIMENSION(0:16), PARAMETER :: albedo_type_name = (/ &
	126	'user defined ', & ! 0
	127	'ocean ', & ! 1
	128	'mixed farming, tall grassland ', & ! 2
	129	'tall/medium grassland ', & ! 3
	130	'evergreen shrubland ', & ! 4
	131	'short grassland/meadow/shrubland ', & ! 5
	132	'evergreen needleleaf forest ', & ! 6
	133	'mixed deciduous evergreen forest ', & ! 7
	134	'deciduous forest ', & ! 8
	135	'tropical evergreen broadleaved forest', & ! 9
	136	'medium/tall grassland/woodland ', & ! 10
	137	'desert, sandy ', & ! 11
	138	'desert, rocky ', & ! 12
	139	'tundra ', & ! 13
	140	'land ice ', & ! 14
	141	'sea ice ', & ! 15
	142	'snow ' & ! 16
[1585]	143	/)
[1496]	144
[1682]	145	INTEGER(iwp) :: albedo_type = 5, & !< Albedo surface type (default: short grassland)
	146	day, & !< current day of the year
	147	day_init = 172, & !< day of the year at model start (21/06)
	148	dots_rad = 0 !< starting index for timeseries output
[1496]	149
	150
	151
[1585]	152
	153
	154
[1682]	155	LOGICAL :: constant_albedo = .FALSE., & !< flag parameter indicating whether the albedo may change depending on zenith
	156	lw_radiation = .TRUE., & !< flag parameter indicing whether longwave radiation shall be calculated
	157	radiation = .FALSE., & !< flag parameter indicating whether the radiation model is used
	158	sun_up = .TRUE., & !< flag parameter indicating whether the sun is up or down
	159	sw_radiation = .TRUE. !< flag parameter indicing whether shortwave radiation shall be calculated
[1585]	160
[1682]	161	REAL(wp), PARAMETER :: sigma_sb = 5.67E-8_wp, & !< Stefan-Boltzmann constant
	162	solar_constant = 1368.0_wp !< solar constant at top of atmosphere
[1496]	163
[1682]	164	REAL(wp) :: albedo = 9999999.9_wp, & !< NAMELIST alpha
	165	albedo_lw_dif = 9999999.9_wp, & !< NAMELIST aldif
	166	albedo_lw_dir = 9999999.9_wp, & !< NAMELIST aldir
	167	albedo_sw_dif = 9999999.9_wp, & !< NAMELIST asdif
	168	albedo_sw_dir = 9999999.9_wp, & !< NAMELIST asdir
	169	dt_radiation = 0.0_wp, & !< radiation model timestep
	170	emissivity = 0.95_wp, & !< NAMELIST surface emissivity
	171	exn, & !< Exner function
	172	lon = 0.0_wp, & !< longitude in radians
	173	lat = 0.0_wp, & !< latitude in radians
	174	decl_1, & !< declination coef. 1
	175	decl_2, & !< declination coef. 2
	176	decl_3, & !< declination coef. 3
	177	time_utc, & !< current time in UTC
	178	time_utc_init = 43200.0_wp, & !< UTC time at model start (noon)
	179	lambda = 0.0_wp, & !< longitude in degrees
	180	net_radiation = 0.0_wp, & !< net radiation at surface
	181	time_radiation = 0.0_wp, & !< time since last call of radiation code
	182	sky_trans !< sky transmissivity
[1496]	183
[1585]	184
[1682]	185	REAL(wp), DIMENSION(0:0) :: zenith !< solar zenith angle
[1585]	186
[1496]	187	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
[1682]	188	alpha, & !< surface broadband albedo (used for clear-sky scheme)
	189	rad_net, & !< net radiation at the surface
	190	rad_net_av !< average of rad_net
[1496]	191
[1585]	192	!
	193	!-- Land surface albedos for solar zenith angle of 60Â° after Briegleb (1992)
	194	!-- (shortwave, longwave, broadband): sw, lw, bb,
[1587]	195	REAL(wp), DIMENSION(0:2,1:16), PARAMETER :: albedo_pars = RESHAPE( (/&
[1585]	196	0.06_wp, 0.06_wp, 0.06_wp, & ! 1
	197	0.09_wp, 0.28_wp, 0.19_wp, & ! 2
	198	0.11_wp, 0.33_wp, 0.23_wp, & ! 3
	199	0.11_wp, 0.33_wp, 0.23_wp, & ! 4
	200	0.14_wp, 0.34_wp, 0.25_wp, & ! 5
	201	0.06_wp, 0.22_wp, 0.14_wp, & ! 6
	202	0.06_wp, 0.27_wp, 0.17_wp, & ! 7
	203	0.06_wp, 0.31_wp, 0.19_wp, & ! 8
	204	0.06_wp, 0.22_wp, 0.14_wp, & ! 9
	205	0.06_wp, 0.28_wp, 0.18_wp, & ! 10
	206	0.35_wp, 0.51_wp, 0.43_wp, & ! 11
	207	0.24_wp, 0.40_wp, 0.32_wp, & ! 12
	208	0.10_wp, 0.27_wp, 0.19_wp, & ! 13
	209	0.90_wp, 0.65_wp, 0.77_wp, & ! 14
[1587]	210	0.90_wp, 0.65_wp, 0.77_wp, & ! 15
	211	0.95_wp, 0.70_wp, 0.82_wp & ! 16
	212	/), (/ 3, 16 /) )
[1496]	213
[1585]	214	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: &
[1682]	215	rad_sw_in, & !< incoming shortwave radiation (W/m2)
	216	rad_sw_in_av, & !< average of rad_sw_in
	217	rad_sw_out, & !< outgoing shortwave radiation (W/m2)
	218	rad_sw_out_av, & !< average of rad_sw_out
	219	rad_lw_in, & !< incoming longwave radiation (W/m2)
	220	rad_lw_in_av, & !< average of rad_lw_in
	221	rad_lw_out, & !< outgoing longwave radiation (W/m2)
	222	rad_lw_out_av !< average of rad_lw_out
[1585]	223
	224	!
	225	!-- Variables and parameters used in RRTMG only
	226	#if defined ( __rrtmg )
[1682]	227	CHARACTER(LEN=12) :: rrtm_input_file = "RAD_SND_DATA" !< name of the NetCDF input file (sounding data)
[1585]	228
	229
	230	!
	231	!-- Flag parameters for RRTMGS (should not be changed)
[1682]	232	INTEGER(iwp), PARAMETER :: rrtm_inflglw = 2, & !< flag for lw cloud optical properties (0,1,2)
	233	rrtm_iceflglw = 0, & !< flag for lw ice particle specifications (0,1,2,3)
	234	rrtm_liqflglw = 1, & !< flag for lw liquid droplet specifications
	235	rrtm_inflgsw = 2, & !< flag for sw cloud optical properties (0,1,2)
	236	rrtm_iceflgsw = 0, & !< flag for sw ice particle specifications (0,1,2,3)
	237	rrtm_liqflgsw = 1 !< flag for sw liquid droplet specifications
[1585]	238
	239	!
	240	!-- The following variables should be only changed with care, as this will
	241	!-- require further setting of some variables, which is currently not
	242	!-- implemented (aerosols, ice phase).
[1682]	243	INTEGER(iwp) :: nzt_rad, & !< upper vertical limit for radiation calculations
	244	rrtm_icld = 0, & !< cloud flag (0: clear sky column, 1: cloudy column)
	245	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)
	246	rrtm_idrv = 0 !< longwave upward flux calculation option (0,1)
[1585]	247
[1682]	248	LOGICAL :: snd_exists = .FALSE. !< flag parameter to check whether a user-defined input files exists
[1585]	249
[1682]	250	REAL(wp), PARAMETER :: mol_weight_air_d_wv = 1.607793_wp !< molecular weight dry air / water vapor
[1585]	251
[1682]	252	REAL(wp), DIMENSION(:), ALLOCATABLE :: hyp_snd, & !< hypostatic pressure from sounding data (hPa)
	253	q_snd, & !< specific humidity from sounding data (kg/kg) - dummy at the moment
	254	rrtm_tsfc, & !< dummy array for storing surface temperature
	255	t_snd !< actual temperature from sounding data (hPa)
[1585]	256
[1682]	257	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: aldif, & !< longwave diffuse albedo solar angle of 60Â°
	258	aldir, & !< longwave direct albedo solar angle of 60Â°
	259	asdif, & !< shortwave diffuse albedo solar angle of 60Â°
	260	asdir, & !< shortwave direct albedo solar angle of 60Â°
	261	rrtm_ccl4vmr, & !< CCL4 volume mixing ratio (g/mol)
	262	rrtm_cfc11vmr, & !< CFC11 volume mixing ratio (g/mol)
	263	rrtm_cfc12vmr, & !< CFC12 volume mixing ratio (g/mol)
	264	rrtm_cfc22vmr, & !< CFC22 volume mixing ratio (g/mol)
	265	rrtm_ch4vmr, & !< CH4 volume mixing ratio
	266	rrtm_cicewp, & !< in-cloud ice water path (g/mÂ²)
	267	rrtm_cldfr, & !< cloud fraction (0,1)
	268	rrtm_cliqwp, & !< in-cloud liquid water path (g/mÂ²)
	269	rrtm_co2vmr, & !< CO2 volume mixing ratio (g/mol)
	270	rrtm_emis, & !< surface emissivity (0-1)
	271	rrtm_h2ovmr, & !< H2O volume mixing ratio
	272	rrtm_n2ovmr, & !< N2O volume mixing ratio
	273	rrtm_o2vmr, & !< O2 volume mixing ratio
	274	rrtm_o3vmr, & !< O3 volume mixing ratio
	275	rrtm_play, & !< pressure layers (hPa, zu-grid)
	276	rrtm_plev, & !< pressure layers (hPa, zw-grid)
	277	rrtm_reice, & !< cloud ice effective radius (microns)
	278	rrtm_reliq, & !< cloud water drop effective radius (microns)
	279	rrtm_tlay, & !< actual temperature (K, zu-grid)
	280	rrtm_tlev, & !< actual temperature (K, zw-grid)
	281	rrtm_lwdflx, & !< RRTM output of incoming longwave radiation flux (W/m2)
	282	rrtm_lwuflx, & !< RRTM output of outgoing longwave radiation flux (W/m2)
	283	rrtm_lwhr, & !< RRTM output of longwave radiation heating rate (K/d)
	284	rrtm_lwuflxc, & !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
	285	rrtm_lwdflxc, & !< RRTM output of outgoing clear sky longwave radiation flux (W/m2)
	286	rrtm_lwhrc, & !< RRTM output of incoming longwave clear sky radiation heating rate (K/d)
	287	rrtm_swdflx, & !< RRTM output of incoming shortwave radiation flux (W/m2)
	288	rrtm_swuflx, & !< RRTM output of outgoing shortwave radiation flux (W/m2)
	289	rrtm_swhr, & !< RRTM output of shortwave radiation heating rate (K/d)
	290	rrtm_swuflxc, & !< RRTM output of incoming clear sky shortwave radiation flux (W/m2)
	291	rrtm_swdflxc, & !< RRTM output of outgoing clear sky shortwave radiation flux (W/m2)
	292	rrtm_swhrc !< RRTM output of incoming shortwave clear sky radiation heating rate (K/d)
[1585]	293
	294	!
	295	!-- Definition of arrays that are currently not used for calling RRTMG (due to setting of flag parameters)
[1682]	296	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: rad_lw_cs_in, & !< incoming clear sky longwave radiation (W/m2) (not used)
	297	rad_lw_cs_out, & !< outgoing clear sky longwave radiation (W/m2) (not used)
	298	rad_lw_cs_hr, & !< longwave clear sky radiation heating rate (K/d) (not used)
	299	rad_lw_hr, & !< longwave radiation heating rate (K/d)
	300	rad_sw_cs_in, & !< incoming clear sky shortwave radiation (W/m2) (not used)
	301	rad_sw_cs_out, & !< outgoing clear sky shortwave radiation (W/m2) (not used)
	302	rad_sw_cs_hr, & !< shortwave clear sky radiation heating rate (K/d) (not used)
	303	rad_sw_hr, & !< shortwave radiation heating rate (K/d)
	304	rrtm_aldif, & !< surface albedo for longwave diffuse radiation
	305	rrtm_aldir, & !< surface albedo for longwave direct radiation
	306	rrtm_asdif, & !< surface albedo for shortwave diffuse radiation
	307	rrtm_asdir, & !< surface albedo for shortwave direct radiation
	308	rrtm_lw_tauaer, & !< lw aerosol optical depth
	309	rrtm_lw_taucld, & !< lw in-cloud optical depth
	310	rrtm_sw_taucld, & !< sw in-cloud optical depth
	311	rrtm_sw_ssacld, & !< sw in-cloud single scattering albedo
	312	rrtm_sw_asmcld, & !< sw in-cloud asymmetry parameter
	313	rrtm_sw_fsfcld, & !< sw in-cloud forward scattering fraction
	314	rrtm_sw_tauaer, & !< sw aerosol optical depth
	315	rrtm_sw_ssaaer, & !< sw aerosol single scattering albedo
	316	rrtm_sw_asmaer, & !< sw aerosol asymmetry parameter
	317	rrtm_sw_ecaer !< sw aerosol optical detph at 0.55 microns (rrtm_iaer = 6 only)
[1585]	318	#endif
	319
[1496]	320	INTERFACE init_radiation
	321	MODULE PROCEDURE init_radiation
	322	END INTERFACE init_radiation
	323
[1551]	324	INTERFACE radiation_clearsky
	325	MODULE PROCEDURE radiation_clearsky
	326	END INTERFACE radiation_clearsky
[1496]	327
[1585]	328	INTERFACE radiation_rrtmg
	329	MODULE PROCEDURE radiation_rrtmg
	330	END INTERFACE radiation_rrtmg
[1551]	331
[1585]	332	INTERFACE radiation_tendency
	333	MODULE PROCEDURE radiation_tendency
	334	MODULE PROCEDURE radiation_tendency_ij
	335	END INTERFACE radiation_tendency
[1551]	336
[1496]	337	SAVE
	338
	339	PRIVATE
	340
[1585]	341	PUBLIC albedo, albedo_type, albedo_type_name, albedo_lw_dif, albedo_lw_dir,&
	342	albedo_sw_dif, albedo_sw_dir, constant_albedo, day_init, dots_rad, &
	343	dt_radiation, init_radiation, lambda, lw_radiation, net_radiation, &
	344	rad_net, rad_net_av, radiation, radiation_clearsky, &
	345	radiation_rrtmg, radiation_scheme, radiation_tendency, rad_lw_in, &
	346	rad_lw_in_av, rad_lw_out, rad_lw_out_av, rad_sw_in, rad_sw_in_av, &
	347	rad_sw_out, rad_sw_out_av, sigma_sb, sw_radiation, time_radiation, &
	348	time_utc_init
[1496]	349
[1585]	350	#if defined ( __rrtmg )
	351	PUBLIC rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir
	352	#endif
[1496]	353
	354	CONTAINS
	355
	356	!------------------------------------------------------------------------------!
	357	! Description:
	358	! ------------
[1682]	359	!> Initialization of the radiation model
[1496]	360	!------------------------------------------------------------------------------!
	361	SUBROUTINE init_radiation
	362
	363	IMPLICIT NONE
	364
[1585]	365	!
	366	!-- Allocate array for storing the surface net radiation
	367	IF ( .NOT. ALLOCATED ( rad_net ) ) THEN
	368	ALLOCATE ( rad_net(nysg:nyng,nxlg:nxrg) )
	369	rad_net = 0.0_wp
	370	ENDIF
[1496]	371
	372	!
[1551]	373	!-- Fix net radiation in case of radiation_scheme = 'constant'
[1585]	374	IF ( radiation_scheme == 'constant' ) THEN
[1551]	375	rad_net = net_radiation
[1585]	376	radiation = .FALSE.
[1551]	377	!
[1585]	378	!-- Calculate orbital constants
	379	ELSE
[1551]	380	decl_1 = SIN(23.45_wp * pi / 180.0_wp)
	381	decl_2 = 2.0_wp * pi / 365.0_wp
	382	decl_3 = decl_2 * 81.0_wp
[1585]	383	lat = phi * pi / 180.0_wp
	384	lon = lambda * pi / 180.0_wp
	385	ENDIF
	386
	387
	388	IF ( radiation_scheme == 'clear-sky' ) THEN
	389
	390	ALLOCATE ( alpha(nysg:nyng,nxlg:nxrg) )
	391
	392	IF ( .NOT. ALLOCATED ( rad_sw_in ) ) THEN
	393	ALLOCATE ( rad_sw_in(0:0,nysg:nyng,nxlg:nxrg) )
	394	ENDIF
	395	IF ( .NOT. ALLOCATED ( rad_sw_out ) ) THEN
	396	ALLOCATE ( rad_sw_out(0:0,nysg:nyng,nxlg:nxrg) )
	397	ENDIF
	398
	399	IF ( .NOT. ALLOCATED ( rad_sw_in_av ) ) THEN
	400	ALLOCATE ( rad_sw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
	401	ENDIF
	402	IF ( .NOT. ALLOCATED ( rad_sw_out_av ) ) THEN
	403	ALLOCATE ( rad_sw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
	404	ENDIF
	405
	406	IF ( .NOT. ALLOCATED ( rad_lw_in ) ) THEN
	407	ALLOCATE ( rad_lw_in(0:0,nysg:nyng,nxlg:nxrg) )
	408	ENDIF
	409	IF ( .NOT. ALLOCATED ( rad_lw_out ) ) THEN
	410	ALLOCATE ( rad_lw_out(0:0,nysg:nyng,nxlg:nxrg) )
	411	ENDIF
	412
	413	IF ( .NOT. ALLOCATED ( rad_lw_in_av ) ) THEN
	414	ALLOCATE ( rad_lw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
	415	ENDIF
	416	IF ( .NOT. ALLOCATED ( rad_lw_out_av ) ) THEN
	417	ALLOCATE ( rad_lw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
	418	ENDIF
	419
	420	rad_sw_in = 0.0_wp
	421	rad_sw_out = 0.0_wp
	422	rad_lw_in = 0.0_wp
	423	rad_lw_out = 0.0_wp
	424
[1496]	425	!
[1585]	426	!-- Overwrite albedo if manually set in parameter file
	427	IF ( albedo_type /= 0 .AND. albedo == 9999999.9_wp ) THEN
	428	albedo = albedo_pars(2,albedo_type)
	429	ENDIF
	430
	431	alpha = albedo
	432
	433	!
	434	!-- Initialization actions for RRTMG
	435	ELSEIF ( radiation_scheme == 'rrtmg' ) THEN
	436	#if defined ( __rrtmg )
	437	!
	438	!-- Allocate albedos
	439	ALLOCATE ( rrtm_aldif(0:0,nysg:nyng,nxlg:nxrg) )
	440	ALLOCATE ( rrtm_aldir(0:0,nysg:nyng,nxlg:nxrg) )
	441	ALLOCATE ( rrtm_asdif(0:0,nysg:nyng,nxlg:nxrg) )
	442	ALLOCATE ( rrtm_asdir(0:0,nysg:nyng,nxlg:nxrg) )
	443	ALLOCATE ( aldif(nysg:nyng,nxlg:nxrg) )
	444	ALLOCATE ( aldir(nysg:nyng,nxlg:nxrg) )
	445	ALLOCATE ( asdif(nysg:nyng,nxlg:nxrg) )
	446	ALLOCATE ( asdir(nysg:nyng,nxlg:nxrg) )
	447
	448	IF ( albedo_type /= 0 ) THEN
	449	IF ( albedo_lw_dif == 9999999.9_wp ) THEN
	450	albedo_lw_dif = albedo_pars(0,albedo_type)
	451	albedo_lw_dir = albedo_lw_dif
	452	ENDIF
	453	IF ( albedo_sw_dif == 9999999.9_wp ) THEN
	454	albedo_sw_dif = albedo_pars(1,albedo_type)
	455	albedo_sw_dir = albedo_sw_dif
	456	ENDIF
	457	ENDIF
	458
	459	aldif(:,:) = albedo_lw_dif
	460	aldir(:,:) = albedo_lw_dir
	461	asdif(:,:) = albedo_sw_dif
	462	asdir(:,:) = albedo_sw_dir
	463	!
	464	!-- Calculate initial values of current (cosine of) the zenith angle and
	465	!-- whether the sun is up
	466	CALL calc_zenith
	467	!
	468	!-- Calculate initial surface albedo
	469	IF ( .NOT. constant_albedo ) THEN
	470	CALL calc_albedo
	471	ELSE
	472	rrtm_aldif(0,:,:) = aldif(:,:)
	473	rrtm_aldir(0,:,:) = aldir(:,:)
	474	rrtm_asdif(0,:,:) = asdif(:,:)
	475	rrtm_asdir(0,:,:) = asdir(:,:)
	476	ENDIF
	477
	478	!
	479	!-- Allocate surface emissivity
	480	ALLOCATE ( rrtm_emis(0:0,1:nbndlw+1) )
	481	rrtm_emis = emissivity
	482
	483	!
	484	!-- Allocate 3d arrays of radiative fluxes and heating rates
	485	ALLOCATE ( rad_sw_hr(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
	486	ALLOCATE ( rad_sw_cs_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	487	ALLOCATE ( rad_sw_cs_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	488	ALLOCATE ( rad_sw_cs_hr(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
	489
	490	IF ( .NOT. ALLOCATED ( rad_sw_in ) ) THEN
	491	ALLOCATE ( rad_sw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	492	rad_sw_in = 0.0_wp
	493
	494	ENDIF
	495
	496	IF ( .NOT. ALLOCATED ( rad_sw_in_av ) ) THEN
	497	ALLOCATE ( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	498	rad_sw_out = 0.0_wp
	499	ENDIF
	500
	501	IF ( .NOT. ALLOCATED ( rad_sw_out ) ) THEN
	502	ALLOCATE ( rad_sw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	503	ENDIF
	504
	505	IF ( .NOT. ALLOCATED ( rad_sw_out_av ) ) THEN
	506	ALLOCATE ( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	507	ENDIF
	508
	509	ALLOCATE ( rad_lw_hr(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
	510	ALLOCATE ( rad_lw_cs_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	511	ALLOCATE ( rad_lw_cs_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	512	ALLOCATE ( rad_lw_cs_hr(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
	513
	514	IF ( .NOT. ALLOCATED ( rad_lw_in ) ) THEN
	515	ALLOCATE ( rad_lw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	516	rad_lw_in = 0.0_wp
	517	ENDIF
	518
	519	IF ( .NOT. ALLOCATED ( rad_lw_in_av ) ) THEN
	520	ALLOCATE ( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	521	ENDIF
	522
	523	IF ( .NOT. ALLOCATED ( rad_lw_out ) ) THEN
	524	ALLOCATE ( rad_lw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	525	rad_lw_out = 0.0_wp
	526	ENDIF
	527
	528	IF ( .NOT. ALLOCATED ( rad_lw_out_av ) ) THEN
	529	ALLOCATE ( rad_lw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	530	ENDIF
	531
	532	rad_sw_hr = 0.0_wp
	533	rad_sw_cs_in = 0.0_wp
	534	rad_sw_cs_out = 0.0_wp
	535	rad_sw_cs_hr = 0.0_wp
	536
	537	rad_lw_hr = 0.0_wp
	538	rad_lw_cs_in = 0.0_wp
	539	rad_lw_cs_out = 0.0_wp
	540	rad_lw_cs_hr = 0.0_wp
	541
	542	!
	543	!-- Allocate dummy array for storing surface temperature
	544	ALLOCATE ( rrtm_tsfc(1) )
	545
	546	!
	547	!-- Initialize RRTMG
	548	IF ( lw_radiation ) CALL rrtmg_lw_ini ( cp )
	549	IF ( sw_radiation ) CALL rrtmg_sw_ini ( cp )
	550
	551	!
	552	!-- Set input files for RRTMG
	553	INQUIRE(FILE="RAD_SND_DATA", EXIST=snd_exists)
	554	IF ( .NOT. snd_exists ) THEN
	555	rrtm_input_file = "rrtmg_lw.nc"
	556	ENDIF
	557
	558	!
	559	!-- Read vertical layers for RRTMG from sounding data
	560	!-- The routine provides nzt_rad, hyp_snd(1:nzt_rad),
	561	!-- t_snd(nzt+2:nzt_rad), rrtm_play(1:nzt_rad), rrtm_plev(1_nzt_rad+1),
	562	!-- rrtm_tlay(nzt+2:nzt_rad), rrtm_tlev(nzt+2:nzt_rad+1)
	563	CALL read_sounding_data
	564
	565	!
	566	!-- Read trace gas profiles from file. This routine provides
	567	!-- the rrtm_ arrays (1:nzt_rad+1)
	568	CALL read_trace_gas_data
	569	#endif
[1551]	570	ENDIF
[1585]	571
[1551]	572	!
[1585]	573	!-- Perform user actions if required
	574	CALL user_init_radiation
	575
	576
	577	!
[1551]	578	!-- Add timeseries for radiation model
[1585]	579	dots_rad = dots_num + 1
[1551]	580	dots_label(dots_num+1) = "rad_net"
[1585]	581	dots_label(dots_num+2) = "rad_lw_in"
	582	dots_label(dots_num+3) = "rad_lw_out"
	583	dots_label(dots_num+4) = "rad_sw_in"
	584	dots_label(dots_num+5) = "rad_sw_out"
	585	dots_unit(dots_num+1:dots_num+5) = "W/m2"
	586	dots_num = dots_num + 5
[1496]	587
[1585]	588	!
	589	!-- Output of albedos is only required for RRTMG
	590	IF ( radiation_scheme == 'rrtmg' ) THEN
	591	dots_label(dots_num+1) = "rrtm_aldif"
	592	dots_label(dots_num+2) = "rrtm_aldir"
	593	dots_label(dots_num+3) = "rrtm_asdif"
	594	dots_label(dots_num+4) = "rrtm_asdir"
	595	dots_unit(dots_num+1:dots_num+4) = ""
	596	dots_num = dots_num + 4
	597	ENDIF
[1551]	598
[1585]	599	!
	600	!-- Calculate radiative fluxes at model start
	601	IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
	602	IF ( radiation_scheme == 'clear-sky' ) THEN
	603	CALL radiation_clearsky
	604	ELSEIF ( radiation_scheme == 'rrtmg' ) THEN
	605	CALL radiation_rrtmg
	606	ENDIF
	607	ENDIF
	608
[1496]	609	RETURN
	610
	611	END SUBROUTINE init_radiation
	612
	613
	614	!------------------------------------------------------------------------------!
	615	! Description:
	616	! ------------
[1682]	617	!> A simple clear sky radiation model
[1496]	618	!------------------------------------------------------------------------------!
[1551]	619	SUBROUTINE radiation_clearsky
[1496]	620
[1585]	621	USE indices, &
	622	ONLY: nbgp
	623
[1496]	624	IMPLICIT NONE
	625
[1585]	626	INTEGER(iwp) :: i, j, k
	627
[1496]	628	!
[1585]	629	!-- Calculate current zenith angle
	630	CALL calc_zenith
	631
	632	!
	633	!-- Calculate sky transmissivity
	634	sky_trans = 0.6_wp + 0.2_wp * zenith(0)
	635
	636	!
	637	!-- Calculate value of the Exner function
	638	exn = (surface_pressure / 1000.0_wp )**0.286_wp
	639	!
	640	!-- Calculate radiation fluxes and net radiation (rad_net) for each grid
	641	!-- point
	642	DO i = nxl, nxr
	643	DO j = nys, nyn
	644	k = nzb_s_inner(j,i)
	645	rad_sw_in(0,j,i) = solar_constant * sky_trans * zenith(0)
	646	rad_sw_out(0,j,i) = alpha(j,i) * rad_sw_in(0,j,i)
	647	rad_lw_out(0,j,i) = sigma_sb * (pt(k,j,i) * exn)**4
	648	rad_lw_in(0,j,i) = 0.8_wp * sigma_sb * (pt(k+1,j,i) * exn)**4
	649	rad_net(j,i) = rad_sw_in(0,j,i) - rad_sw_out(0,j,i) &
	650	+ rad_lw_in(0,j,i) - rad_lw_out(0,j,i)
	651
	652	ENDDO
	653	ENDDO
	654
	655	CALL exchange_horiz_2d( rad_lw_in, nbgp )
	656	CALL exchange_horiz_2d( rad_lw_out, nbgp )
	657	CALL exchange_horiz_2d( rad_sw_in, nbgp )
	658	CALL exchange_horiz_2d( rad_sw_out, nbgp )
	659	CALL exchange_horiz_2d( rad_net, nbgp )
	660
	661	RETURN
	662
	663	END SUBROUTINE radiation_clearsky
	664
	665
	666	!------------------------------------------------------------------------------!
	667	! Description:
	668	! ------------
[1682]	669	!> Implementation of the RRTMG radiation_scheme
[1585]	670	!------------------------------------------------------------------------------!
	671	SUBROUTINE radiation_rrtmg
	672
	673	USE indices, &
	674	ONLY: nbgp
	675
	676	USE particle_attributes, &
	677	ONLY: grid_particles, number_of_particles, particles, &
	678	particle_advection_start, prt_count
	679
	680	IMPLICIT NONE
	681
	682	#if defined ( __rrtmg )
	683
[1682]	684	INTEGER(iwp) :: i, j, k, n !<
[1585]	685
[1682]	686	REAL(wp) :: s_r2 !< weighted sum over all droplets with r^2
	687	REAL(wp) :: s_r3 !< weighted sum over all droplets with r^3
[1585]	688
	689	!
	690	!-- Calculate current (cosine of) zenith angle and whether the sun is up
	691	CALL calc_zenith
	692	!
	693	!-- Calculate surface albedo
	694	IF ( .NOT. constant_albedo ) THEN
	695	CALL calc_albedo
	696	ENDIF
	697
	698	!
	699	!-- Prepare input data for RRTMG
	700
	701	!
	702	!-- In case of large scale forcing with surface data, calculate new pressure
	703	!-- profile. nzt_rad might be modified by these calls and all required arrays
	704	!-- will then be re-allocated
	705	IF ( large_scale_forcing .AND. lsf_surf ) THEN
	706	CALL read_sounding_data
	707	CALL read_trace_gas_data
	708	ENDIF
	709	!
	710	!-- Loop over all grid points
	711	DO i = nxl, nxr
	712	DO j = nys, nyn
	713
	714	!
	715	!-- Prepare profiles of temperature and H2O volume mixing ratio
	716	rrtm_tlev(0,nzb+1) = pt(nzb,j,i)
	717
	718	DO k = nzb+1, nzt+1
	719	rrtm_tlay(0,k) = pt(k,j,i) * ( (hyp(k) ) / 100000.0_wp &
	720	)*0.286_wp + ( l_v / cp ) ql(k,j,i)
	721	rrtm_h2ovmr(0,k) = mol_weight_air_d_wv * (q(k,j,i) - ql(k,j,i))
	722
	723	ENDDO
	724
	725	!
	726	!-- Avoid temperature/humidity jumps at the top of the LES domain by
	727	!-- linear interpolation from nzt+2 to nzt+7
	728	DO k = nzt+2, nzt+7
	729	rrtm_tlay(0,k) = rrtm_tlay(0,nzt+1) &
	730	+ ( rrtm_tlay(0,nzt+8) - rrtm_tlay(0,nzt+1) ) &
	731	/ ( rrtm_play(0,nzt+8) - rrtm_play(0,nzt+1) ) &
	732	* ( rrtm_play(0,k) - rrtm_play(0,nzt+1) )
	733
	734	rrtm_h2ovmr(0,k) = rrtm_h2ovmr(0,nzt+1) &
	735	+ ( rrtm_h2ovmr(0,nzt+8) - rrtm_h2ovmr(0,nzt+1) )&
	736	/ ( rrtm_play(0,nzt+8) - rrtm_play(0,nzt+1) )&
	737	* ( rrtm_play(0,k) - rrtm_play(0,nzt+1) )
	738
	739	ENDDO
	740
	741	!-- Linear interpolate to zw grid
	742	DO k = nzb+2, nzt+8
	743	rrtm_tlev(0,k) = rrtm_tlay(0,k-1) + (rrtm_tlay(0,k) - &
	744	rrtm_tlay(0,k-1)) &
	745	/ ( rrtm_play(0,k) - rrtm_play(0,k-1) ) &
	746	* ( rrtm_plev(0,k) - rrtm_play(0,k-1) )
	747	ENDDO
	748
	749
	750	!
	751	!-- Calculate liquid water path and cloud fraction for each column.
	752	!-- Note that LWP is required in g/mÂ² instead of kg/kg m.
	753	rrtm_cldfr = 0.0_wp
	754	rrtm_reliq = 0.0_wp
	755	rrtm_cliqwp = 0.0_wp
	756
	757	DO k = nzb+1, nzt+1
	758	rrtm_cliqwp(0,k) = ql(k,j,i) * 1000.0_wp * &
	759	(rrtm_plev(0,k) - rrtm_plev(0,k+1)) * 100.0_wp / g
	760
	761	IF ( rrtm_cliqwp(0,k) .GT. 0 ) THEN
	762	rrtm_cldfr(0,k) = 1.0_wp
	763	rrtm_icld = 1
	764
	765	!
	766	!-- Calculate cloud droplet effective radius
	767	IF ( cloud_physics ) THEN
	768	rrtm_reliq(0,k) = 1.0E6_wp * ( 3.0_wp * ql(k,j,i) &
	769	/ ( 4.0_wp * pi * nc_const * rho_l ) &
	770	)**0.33333333333333_wp &
	771	* EXP( LOG( sigma_gc )**2 )
	772
	773	ELSEIF ( cloud_droplets ) THEN
	774	number_of_particles = prt_count(k,j,i)
	775
	776	IF (number_of_particles <= 0) CYCLE
	777	particles => grid_particles(k,j,i)%particles(1:number_of_particles)
	778	s_r2 = 0.0_wp
	779	s_r3 = 0.0_wp
	780
	781	DO n = 1, number_of_particles
	782	IF ( particles(n)%particle_mask ) THEN
	783	s_r2 = s_r2 + particles(n)%radius*2 &
	784	particles(n)%weight_factor
	785	s_r3 = s_r3 + particles(n)%radius*3 &
	786	particles(n)%weight_factor
	787	ENDIF
	788	ENDDO
	789
	790	IF ( s_r2 > 0.0_wp ) rrtm_reliq(0,k) = s_r3 / s_r2
	791
	792	ENDIF
	793
	794	!
	795	!-- Limit effective radius
	796	IF ( rrtm_reliq(0,k) .GT. 0.0_wp ) THEN
	797	rrtm_reliq(0,k) = MAX(rrtm_reliq(0,k),2.5_wp)
	798	rrtm_reliq(0,k) = MIN(rrtm_reliq(0,k),60.0_wp)
	799	ENDIF
	800	ELSE
	801	rrtm_icld = 0
	802	ENDIF
	803	ENDDO
	804
	805	!
	806	!-- Set surface temperature
	807	rrtm_tsfc = pt(nzb,j,i) * (surface_pressure / 1000.0_wp )**0.286_wp
	808
	809	IF ( lw_radiation ) THEN
	810	CALL rrtmg_lw( 1, nzt_rad , rrtm_icld , rrtm_idrv ,&
	811	rrtm_play , rrtm_plev , rrtm_tlay , rrtm_tlev ,&
	812	rrtm_tsfc , rrtm_h2ovmr , rrtm_o3vmr , rrtm_co2vmr ,&
	813	rrtm_ch4vmr , rrtm_n2ovmr , rrtm_o2vmr , rrtm_cfc11vmr ,&
	814	rrtm_cfc12vmr , rrtm_cfc22vmr, rrtm_ccl4vmr , rrtm_emis ,&
	815	rrtm_inflglw , rrtm_iceflglw, rrtm_liqflglw, rrtm_cldfr ,&
	816	rrtm_lw_taucld , rrtm_cicewp , rrtm_cliqwp , rrtm_reice ,&
	817	rrtm_reliq , rrtm_lw_tauaer, &
	818	rrtm_lwuflx , rrtm_lwdflx , rrtm_lwhr , &
	819	rrtm_lwuflxc , rrtm_lwdflxc , rrtm_lwhrc )
	820
	821	DO k = nzb, nzt+1
	822	rad_lw_in(k,j,i) = rrtm_lwdflx(0,k)
	823	rad_lw_out(k,j,i) = rrtm_lwuflx(0,k)
	824	ENDDO
	825
	826
	827	ENDIF
	828
	829	IF ( sw_radiation .AND. sun_up ) THEN
	830	CALL rrtmg_sw( 1, nzt_rad , rrtm_icld , rrtm_iaer ,&
	831	rrtm_play , rrtm_plev , rrtm_tlay , rrtm_tlev ,&
	832	rrtm_tsfc , rrtm_h2ovmr , rrtm_o3vmr , rrtm_co2vmr ,&
	833	rrtm_ch4vmr , rrtm_n2ovmr , rrtm_o2vmr , rrtm_asdir(:,j,i),&
	834	rrtm_asdif(:,j,i), rrtm_aldir(:,j,i), rrtm_aldif(:,j,i), zenith,&
	835	0.0_wp , day , solar_constant, rrtm_inflgsw,&
	836	rrtm_iceflgsw , rrtm_liqflgsw, rrtm_cldfr , rrtm_sw_taucld ,&
	837	rrtm_sw_ssacld , rrtm_sw_asmcld, rrtm_sw_fsfcld, rrtm_cicewp ,&
	838	rrtm_cliqwp , rrtm_reice , rrtm_reliq , rrtm_sw_tauaer ,&
	839	rrtm_sw_ssaaer , rrtm_sw_asmaer , rrtm_sw_ecaer , &
	840	rrtm_swuflx , rrtm_swdflx , rrtm_swhr , &
	841	rrtm_swuflxc , rrtm_swdflxc , rrtm_swhrc )
	842
	843	DO k = nzb, nzt+1
	844	rad_sw_in(k,j,i) = rrtm_swdflx(0,k)
	845	rad_sw_out(k,j,i) = rrtm_swuflx(0,k)
	846	ENDDO
	847	ENDIF
	848
	849	!
	850	!-- Calculate surface net radiation
	851	rad_net(j,i) = rad_sw_in(nzb,j,i) - rad_sw_out(nzb,j,i) &
	852	+ rad_lw_in(nzb,j,i) - rad_lw_out(nzb,j,i)
	853
	854	ENDDO
	855	ENDDO
	856
	857	CALL exchange_horiz( rad_lw_in, nbgp )
	858	CALL exchange_horiz( rad_lw_out, nbgp )
	859	CALL exchange_horiz( rad_sw_in, nbgp )
	860	CALL exchange_horiz( rad_sw_out, nbgp )
	861	CALL exchange_horiz_2d( rad_net, nbgp )
	862	#endif
	863
	864	END SUBROUTINE radiation_rrtmg
	865
	866
	867	!------------------------------------------------------------------------------!
	868	! Description:
	869	! ------------
[1682]	870	!> Calculate the cosine of the zenith angle (variable is called zenith)
[1585]	871	!------------------------------------------------------------------------------!
	872	SUBROUTINE calc_zenith
	873
	874	IMPLICIT NONE
	875
[1682]	876	REAL(wp) :: declination, & !< solar declination angle
	877	hour_angle !< solar hour angle
[1585]	878	!
[1496]	879	!-- Calculate current day and time based on the initial values and simulation
	880	!-- time
[1585]	881	day = day_init + INT(FLOOR( (time_utc_init + time_since_reference_point) &
	882	/ 86400.0_wp ), KIND=iwp)
[1496]	883	time_utc = MOD((time_utc_init + time_since_reference_point), 86400.0_wp)
	884
	885
	886	!
	887	!-- Calculate solar declination and hour angle
[1585]	888	declination = ASIN( decl_1 * SIN(decl_2 * REAL(day, KIND=wp) - decl_3) )
[1496]	889	hour_angle = 2.0_wp * pi * (time_utc / 86400.0_wp) + lon - pi
	890
	891	!
	892	!-- Calculate zenith angle
[1585]	893	zenith(0) = SIN(lat) * SIN(declination) + COS(lat) * COS(declination) &
[1496]	894	* COS(hour_angle)
[1585]	895	zenith(0) = MAX(0.0_wp,zenith(0))
[1496]	896
	897	!
[1585]	898	!-- Check if the sun is up (otheriwse shortwave calculations can be skipped)
	899	IF ( zenith(0) .GT. 0.0_wp ) THEN
	900	sun_up = .TRUE.
	901	ELSE
	902	sun_up = .FALSE.
	903	END IF
[1496]	904
[1585]	905	END SUBROUTINE calc_zenith
	906
[1606]	907	#if defined ( __rrtmg ) && defined ( __netcdf )
[1585]	908	!------------------------------------------------------------------------------!
	909	! Description:
	910	! ------------
[1682]	911	!> Calculates surface albedo components based on Briegleb (1992) and
	912	!> Briegleb et al. (1986)
[1585]	913	!------------------------------------------------------------------------------!
	914	SUBROUTINE calc_albedo
	915
	916	IMPLICIT NONE
	917
	918	IF ( sun_up ) THEN
[1496]	919	!
[1585]	920	!-- Ocean
	921	IF ( albedo_type == 1 ) THEN
	922	rrtm_aldir(0,:,:) = 0.026_wp / ( zenith(0)**1.7_wp + 0.065_wp ) &
	923	+ 0.15_wp * ( zenith(0) - 0.1_wp ) &
	924	* ( zenith(0) - 0.5_wp ) &
	925	* ( zenith(0) - 1.0_wp )
	926	rrtm_asdir(0,:,:) = rrtm_aldir(0,:,:)
	927	!
	928	!-- Snow
	929	ELSEIF ( albedo_type == 16 ) THEN
	930	IF ( zenith(0) < 0.5_wp ) THEN
	931	rrtm_aldir(0,:,:) = 0.5_wp * (1.0_wp - aldif) &
	932	* ( 3.0_wp / (1.0_wp + 4.0_wp &
	933	* zenith(0))) - 1.0_wp
	934	rrtm_asdir(0,:,:) = 0.5_wp * (1.0_wp - asdif) &
	935	* ( 3.0_wp / (1.0_wp + 4.0_wp &
	936	* zenith(0))) - 1.0_wp
[1496]	937
[1585]	938	rrtm_aldir(0,:,:) = MIN(0.98_wp, rrtm_aldir(0,:,:))
	939	rrtm_asdir(0,:,:) = MIN(0.98_wp, rrtm_asdir(0,:,:))
	940	ELSE
	941	rrtm_aldir(0,:,:) = aldif
	942	rrtm_asdir(0,:,:) = asdif
	943	ENDIF
[1496]	944	!
[1585]	945	!-- Sea ice
	946	ELSEIF ( albedo_type == 15 ) THEN
	947	rrtm_aldir(0,:,:) = aldif
	948	rrtm_asdir(0,:,:) = asdif
	949	!
	950	!-- Land surfaces
	951	ELSE
	952	SELECT CASE ( albedo_type )
[1496]	953
[1585]	954	!
	955	!-- Surface types with strong zenith dependence
	956	CASE ( 1, 2, 3, 4, 11, 12, 13 )
	957	rrtm_aldir(0,:,:) = aldif * 1.4_wp / &
	958	(1.0_wp + 0.8_wp * zenith(0))
	959	rrtm_asdir(0,:,:) = asdif * 1.4_wp / &
	960	(1.0_wp + 0.8_wp * zenith(0))
	961	!
	962	!-- Surface types with weak zenith dependence
	963	CASE ( 5, 6, 7, 8, 9, 10, 14 )
	964	rrtm_aldir(0,:,:) = aldif * 1.1_wp / &
	965	(1.0_wp + 0.2_wp * zenith(0))
	966	rrtm_asdir(0,:,:) = asdif * 1.1_wp / &
	967	(1.0_wp + 0.2_wp * zenith(0))
[1496]	968
[1585]	969	CASE DEFAULT
	970
	971	END SELECT
	972	ENDIF
	973	!
	974	!-- Diffusive albedo is taken from Table 2
	975	rrtm_aldif(0,:,:) = aldif
	976	rrtm_asdif(0,:,:) = asdif
	977
	978	ELSE
	979
	980	rrtm_aldir(0,:,:) = 0.0_wp
	981	rrtm_asdir(0,:,:) = 0.0_wp
	982	rrtm_aldif(0,:,:) = 0.0_wp
	983	rrtm_asdif(0,:,:) = 0.0_wp
	984	ENDIF
	985	END SUBROUTINE calc_albedo
	986
	987	!------------------------------------------------------------------------------!
	988	! Description:
	989	! ------------
[1682]	990	!> Read sounding data (pressure and temperature) from RADIATION_DATA.
[1585]	991	!------------------------------------------------------------------------------!
	992	SUBROUTINE read_sounding_data
	993
	994	USE netcdf_control
	995
	996	IMPLICIT NONE
	997
	998	INTEGER(iwp) :: id, id_dim_zrad, id_var, k, nz_snd, nz_snd_start, nz_snd_end
	999
	1000	REAL(wp) :: t_surface
	1001
	1002	REAL(wp), DIMENSION(:), ALLOCATABLE :: hyp_snd_tmp, t_snd_tmp
	1003
	1004	!
	1005	!-- In case of updates, deallocate arrays first (sufficient to check one
	1006	!-- array as the others are automatically allocated). This is required
	1007	!-- because nzt_rad might change during the update
	1008	IF ( ALLOCATED ( hyp_snd ) ) THEN
	1009	DEALLOCATE( hyp_snd )
	1010	DEALLOCATE( t_snd )
	1011	DEALLOCATE( q_snd )
	1012	DEALLOCATE ( rrtm_play )
	1013	DEALLOCATE ( rrtm_plev )
	1014	DEALLOCATE ( rrtm_tlay )
	1015	DEALLOCATE ( rrtm_tlev )
	1016	DEALLOCATE ( rrtm_h2ovmr )
	1017	DEALLOCATE ( rrtm_cicewp )
	1018	DEALLOCATE ( rrtm_cldfr )
	1019	DEALLOCATE ( rrtm_cliqwp )
	1020	DEALLOCATE ( rrtm_reice )
	1021	DEALLOCATE ( rrtm_reliq )
	1022	DEALLOCATE ( rrtm_lw_taucld )
	1023	DEALLOCATE ( rrtm_lw_tauaer )
	1024	DEALLOCATE ( rrtm_lwdflx )
	1025	DEALLOCATE ( rrtm_lwuflx )
	1026	DEALLOCATE ( rrtm_lwhr )
	1027	DEALLOCATE ( rrtm_lwuflxc )
	1028	DEALLOCATE ( rrtm_lwdflxc )
	1029	DEALLOCATE ( rrtm_lwhrc )
	1030	DEALLOCATE ( rrtm_sw_taucld )
	1031	DEALLOCATE ( rrtm_sw_ssacld )
	1032	DEALLOCATE ( rrtm_sw_asmcld )
	1033	DEALLOCATE ( rrtm_sw_fsfcld )
	1034	DEALLOCATE ( rrtm_sw_tauaer )
	1035	DEALLOCATE ( rrtm_sw_ssaaer )
	1036	DEALLOCATE ( rrtm_sw_asmaer )
	1037	DEALLOCATE ( rrtm_sw_ecaer )
	1038	DEALLOCATE ( rrtm_swdflx )
	1039	DEALLOCATE ( rrtm_swuflx )
	1040	DEALLOCATE ( rrtm_swhr )
	1041	DEALLOCATE ( rrtm_swuflxc )
	1042	DEALLOCATE ( rrtm_swdflxc )
	1043	DEALLOCATE ( rrtm_swhrc )
	1044	ENDIF
	1045
	1046	!
	1047	!-- Open file for reading
	1048	nc_stat = NF90_OPEN( rrtm_input_file, NF90_NOWRITE, id )
	1049	CALL handle_netcdf_error( 'netcdf', 549 )
	1050
	1051	!
	1052	!-- Inquire dimension of z axis and save in nz_snd
	1053	nc_stat = NF90_INQ_DIMID( id, "Pressure", id_dim_zrad )
	1054	nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim_zrad, len = nz_snd )
	1055	CALL handle_netcdf_error( 'netcdf', 551 )
	1056
	1057	!
	1058	! !-- Allocate temporary array for storing pressure data
	1059	ALLOCATE( hyp_snd_tmp(nzb+1:nz_snd) )
	1060	hyp_snd_tmp = 0.0_wp
	1061
	1062
	1063	!-- Read pressure from file
	1064	nc_stat = NF90_INQ_VARID( id, "Pressure", id_var )
	1065	nc_stat = NF90_GET_VAR( id, id_var, hyp_snd_tmp(:), start = (/1/), &
	1066	count = (/nz_snd/) )
	1067	CALL handle_netcdf_error( 'netcdf', 552 )
	1068
	1069	!
	1070	!-- Allocate temporary array for storing temperature data
	1071	ALLOCATE( t_snd_tmp(nzb+1:nz_snd) )
	1072	t_snd_tmp = 0.0_wp
	1073
	1074	!
	1075	!-- Read temperature from file
	1076	nc_stat = NF90_INQ_VARID( id, "ReferenceTemperature", id_var )
	1077	nc_stat = NF90_GET_VAR( id, id_var, t_snd_tmp(:), start = (/1/), &
	1078	count = (/nz_snd/) )
	1079	CALL handle_netcdf_error( 'netcdf', 553 )
	1080
	1081	!
	1082	!-- Calculate start of sounding data
	1083	nz_snd_start = nz_snd + 1
	1084	nz_snd_end = nz_snd_end
	1085
	1086	!
	1087	!-- Start filling vertical dimension at 10hPa above the model domain (hyp is
	1088	!-- in Pa, hyp_snd in hPa).
	1089	DO k = 1, nz_snd
	1090	IF ( hyp_snd_tmp(k) .LT. (hyp(nzt+1) - 1000.0_wp) * 0.01_wp ) THEN
	1091	nz_snd_start = k
	1092	EXIT
	1093	END IF
	1094	END DO
	1095
	1096	IF ( nz_snd_start .LE. nz_snd ) THEN
	1097	nz_snd_end = nz_snd - 1
	1098	END IF
	1099
	1100
	1101	!
	1102	!-- Calculate of total grid points for RRTMG calculations
	1103	nzt_rad = nzt + nz_snd_end - nz_snd_start + 2
	1104
	1105	!
	1106	!-- Save data above LES domain in hyp_snd, t_snd and q_snd
	1107	!-- Note: q_snd_tmp is not calculated at the moment (dry residual atmosphere)
	1108	ALLOCATE( hyp_snd(nzb+1:nzt_rad) )
	1109	ALLOCATE( t_snd(nzb+1:nzt_rad) )
	1110	ALLOCATE( q_snd(nzb+1:nzt_rad) )
	1111	hyp_snd = 0.0_wp
	1112	t_snd = 0.0_wp
	1113	q_snd = 0.0_wp
	1114
	1115	hyp_snd(nzt+2:nzt_rad) = hyp_snd_tmp(nz_snd_start:nz_snd_end)
	1116	t_snd(nzt+2:nzt_rad) = t_snd_tmp(nz_snd_start:nz_snd_end)
	1117
	1118	DEALLOCATE ( hyp_snd_tmp )
	1119	DEALLOCATE ( t_snd_tmp )
	1120
	1121	nc_stat = NF90_CLOSE( id )
	1122
	1123	!
	1124	!-- Calculate pressure levels on zu and zw grid. Sounding data is added at
	1125	!-- top of the LES domain. This routine does not consider horizontal or
	1126	!-- vertical variability of pressure and temperature
	1127	ALLOCATE ( rrtm_play(0:0,nzb+1:nzt_rad+1) )
	1128	ALLOCATE ( rrtm_plev(0:0,nzb+1:nzt_rad+2) )
	1129
	1130	t_surface = pt_surface * (surface_pressure / 1000.0_wp )**0.286_wp
	1131	DO k = nzb+1, nzt+1
	1132	rrtm_play(0,k) = hyp(k) * 0.01_wp
	1133	rrtm_plev(0,k) = surface_pressure * ( (t_surface - g/cp * zw(k-1)) / &
	1134	t_surface )**(1.0_wp/0.286_wp)
	1135	ENDDO
	1136
	1137	DO k = nzt+2, nzt_rad
	1138	rrtm_play(0,k) = hyp_snd(k)
	1139	rrtm_plev(0,k) = 0.5_wp * ( rrtm_play(0,k) + rrtm_play(0,k-1) )
	1140	ENDDO
	1141	rrtm_plev(0,nzt_rad+1) = MAX( 0.5 * hyp_snd(nzt_rad), &
	1142	1.5 * hyp_snd(nzt_rad) &
	1143	- 0.5 * hyp_snd(nzt_rad-1) )
	1144	rrtm_plev(0,nzt_rad+2) = MIN( 1.0E-4_wp, &
	1145	0.25_wp * rrtm_plev(0,nzt_rad+1) )
	1146
	1147	rrtm_play(0,nzt_rad+1) = 0.5 * rrtm_plev(0,nzt_rad+1)
	1148
	1149	!
	1150	!-- Calculate temperature/humidity levels at top of the LES domain.
	1151	!-- Currently, the temperature is taken from sounding data (might lead to a
	1152	!-- temperature jump at interface. To do: Humidity is currently not
	1153	!-- calculated above the LES domain.
	1154	ALLOCATE ( rrtm_tlay(0:0,nzb+1:nzt_rad+1) )
	1155	ALLOCATE ( rrtm_tlev(0:0,nzb+1:nzt_rad+2) )
	1156	ALLOCATE ( rrtm_h2ovmr(0:0,nzb+1:nzt_rad+1) )
	1157
	1158	DO k = nzt+8, nzt_rad
	1159	rrtm_tlay(0,k) = t_snd(k)
	1160	rrtm_h2ovmr(0,k) = q_snd(k)
	1161	ENDDO
	1162	rrtm_tlay(0,nzt_rad+1) = 2.0_wp * rrtm_tlay(0,nzt_rad) &
	1163	- rrtm_tlay(0,nzt_rad-1)
	1164	DO k = nzt+9, nzt_rad+1
	1165	rrtm_tlev(0,k) = rrtm_tlay(0,k-1) + (rrtm_tlay(0,k) &
	1166	- rrtm_tlay(0,k-1)) &
	1167	/ ( rrtm_play(0,k) - rrtm_play(0,k-1) ) &
	1168	* ( rrtm_plev(0,k) - rrtm_play(0,k-1) )
	1169	ENDDO
	1170	rrtm_h2ovmr(0,nzt_rad+1) = rrtm_h2ovmr(0,nzt_rad)
	1171
	1172	rrtm_tlev(0,nzt_rad+2) = 2.0_wp * rrtm_tlay(0,nzt_rad+1) &
	1173	- rrtm_tlev(0,nzt_rad)
	1174	!
	1175	!-- Allocate remaining RRTMG arrays
	1176	ALLOCATE ( rrtm_cicewp(0:0,nzb+1:nzt_rad+1) )
	1177	ALLOCATE ( rrtm_cldfr(0:0,nzb+1:nzt_rad+1) )
	1178	ALLOCATE ( rrtm_cliqwp(0:0,nzb+1:nzt_rad+1) )
	1179	ALLOCATE ( rrtm_reice(0:0,nzb+1:nzt_rad+1) )
	1180	ALLOCATE ( rrtm_reliq(0:0,nzb+1:nzt_rad+1) )
	1181	ALLOCATE ( rrtm_lw_taucld(1:nbndlw+1,0:0,nzb+1:nzt_rad+1) )
	1182	ALLOCATE ( rrtm_lw_tauaer(0:0,nzb+1:nzt_rad+1,1:nbndlw+1) )
	1183	ALLOCATE ( rrtm_sw_taucld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
	1184	ALLOCATE ( rrtm_sw_ssacld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
	1185	ALLOCATE ( rrtm_sw_asmcld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
	1186	ALLOCATE ( rrtm_sw_fsfcld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
	1187	ALLOCATE ( rrtm_sw_tauaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) )
	1188	ALLOCATE ( rrtm_sw_ssaaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) )
	1189	ALLOCATE ( rrtm_sw_asmaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) )
	1190	ALLOCATE ( rrtm_sw_ecaer(0:0,nzb+1:nzt_rad+1,1:naerec+1) )
	1191
	1192	!
	1193	!-- The ice phase is currently not considered in PALM
	1194	rrtm_cicewp = 0.0_wp
	1195	rrtm_reice = 0.0_wp
	1196
	1197	!
	1198	!-- Set other parameters (move to NAMELIST parameters in the future)
	1199	rrtm_lw_tauaer = 0.0_wp
	1200	rrtm_lw_taucld = 0.0_wp
	1201	rrtm_sw_taucld = 0.0_wp
	1202	rrtm_sw_ssacld = 0.0_wp
	1203	rrtm_sw_asmcld = 0.0_wp
	1204	rrtm_sw_fsfcld = 0.0_wp
	1205	rrtm_sw_tauaer = 0.0_wp
	1206	rrtm_sw_ssaaer = 0.0_wp
	1207	rrtm_sw_asmaer = 0.0_wp
	1208	rrtm_sw_ecaer = 0.0_wp
	1209
	1210
	1211	ALLOCATE ( rrtm_swdflx(0:0,nzb:nzt_rad+1) )
	1212	ALLOCATE ( rrtm_swuflx(0:0,nzb:nzt_rad+1) )
	1213	ALLOCATE ( rrtm_swhr(0:0,nzb+1:nzt_rad+1) )
	1214	ALLOCATE ( rrtm_swuflxc(0:0,nzb:nzt_rad+1) )
	1215	ALLOCATE ( rrtm_swdflxc(0:0,nzb:nzt_rad+1) )
	1216	ALLOCATE ( rrtm_swhrc(0:0,nzb+1:nzt_rad+1) )
	1217
	1218	rrtm_swdflx = 0.0_wp
	1219	rrtm_swuflx = 0.0_wp
	1220	rrtm_swhr = 0.0_wp
	1221	rrtm_swuflxc = 0.0_wp
	1222	rrtm_swdflxc = 0.0_wp
	1223	rrtm_swhrc = 0.0_wp
	1224
	1225	ALLOCATE ( rrtm_lwdflx(0:0,nzb:nzt_rad+1) )
	1226	ALLOCATE ( rrtm_lwuflx(0:0,nzb:nzt_rad+1) )
	1227	ALLOCATE ( rrtm_lwhr(0:0,nzb+1:nzt_rad+1) )
	1228	ALLOCATE ( rrtm_lwuflxc(0:0,nzb:nzt_rad+1) )
	1229	ALLOCATE ( rrtm_lwdflxc(0:0,nzb:nzt_rad+1) )
	1230	ALLOCATE ( rrtm_lwhrc(0:0,nzb+1:nzt_rad+1) )
	1231
	1232	rrtm_lwdflx = 0.0_wp
	1233	rrtm_lwuflx = 0.0_wp
	1234	rrtm_lwhr = 0.0_wp
	1235	rrtm_lwuflxc = 0.0_wp
	1236	rrtm_lwdflxc = 0.0_wp
	1237	rrtm_lwhrc = 0.0_wp
	1238
	1239
	1240	END SUBROUTINE read_sounding_data
	1241
	1242
	1243	!------------------------------------------------------------------------------!
	1244	! Description:
	1245	! ------------
[1682]	1246	!> Read trace gas data from file
[1585]	1247	!------------------------------------------------------------------------------!
	1248	SUBROUTINE read_trace_gas_data
	1249
	1250	USE netcdf_control
	1251	USE rrsw_ncpar
	1252
	1253	IMPLICIT NONE
	1254
[1682]	1255	INTEGER(iwp), PARAMETER :: num_trace_gases = 9 !< number of trace gases
[1585]	1256
	1257	CHARACTER(LEN=5), DIMENSION(num_trace_gases), PARAMETER :: &
	1258	trace_names = (/'O3 ', 'CO2 ', 'CH4 ', 'N2O ', 'O2 ', &
	1259	'CFC11', 'CFC12', 'CFC22', 'CCL4 '/)
	1260
	1261	INTEGER(iwp) :: id, k, m, n, nabs, np, id_abs, id_dim, id_var
	1262
	1263	REAL(wp) :: p_mls_l, p_mls_u, p_wgt_l, p_wgt_u, p_mls_m
	1264
	1265
[1682]	1266	REAL(wp), DIMENSION(:), ALLOCATABLE :: p_mls, & !< pressure levels for the absorbers
	1267	rrtm_play_tmp, & !< temporary array for pressure zu-levels
	1268	rrtm_plev_tmp, & !< temporary array for pressure zw-levels
	1269	trace_path_tmp !< temporary array for storing trace gas path data
[1585]	1270
[1682]	1271	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: trace_mls, & !< array for storing the absorber amounts
	1272	trace_mls_path, & !< array for storing trace gas path data
	1273	trace_mls_tmp !< temporary array for storing trace gas data
[1585]	1274
	1275
	1276	!
	1277	!-- In case of updates, deallocate arrays first (sufficient to check one
	1278	!-- array as the others are automatically allocated)
	1279	IF ( ALLOCATED ( rrtm_o3vmr ) ) THEN
	1280	DEALLOCATE ( rrtm_o3vmr )
	1281	DEALLOCATE ( rrtm_co2vmr )
	1282	DEALLOCATE ( rrtm_ch4vmr )
	1283	DEALLOCATE ( rrtm_n2ovmr )
	1284	DEALLOCATE ( rrtm_o2vmr )
	1285	DEALLOCATE ( rrtm_cfc11vmr )
	1286	DEALLOCATE ( rrtm_cfc12vmr )
	1287	DEALLOCATE ( rrtm_cfc22vmr )
	1288	DEALLOCATE ( rrtm_ccl4vmr )
	1289	ENDIF
	1290
	1291	!
	1292	!-- Allocate trace gas profiles
	1293	ALLOCATE ( rrtm_o3vmr(0:0,1:nzt_rad+1) )
	1294	ALLOCATE ( rrtm_co2vmr(0:0,1:nzt_rad+1) )
	1295	ALLOCATE ( rrtm_ch4vmr(0:0,1:nzt_rad+1) )
	1296	ALLOCATE ( rrtm_n2ovmr(0:0,1:nzt_rad+1) )
	1297	ALLOCATE ( rrtm_o2vmr(0:0,1:nzt_rad+1) )
	1298	ALLOCATE ( rrtm_cfc11vmr(0:0,1:nzt_rad+1) )
	1299	ALLOCATE ( rrtm_cfc12vmr(0:0,1:nzt_rad+1) )
	1300	ALLOCATE ( rrtm_cfc22vmr(0:0,1:nzt_rad+1) )
	1301	ALLOCATE ( rrtm_ccl4vmr(0:0,1:nzt_rad+1) )
	1302
	1303	!
	1304	!-- Open file for reading
	1305	nc_stat = NF90_OPEN( rrtm_input_file, NF90_NOWRITE, id )
	1306	CALL handle_netcdf_error( 'netcdf', 549 )
	1307	!
	1308	!-- Inquire dimension ids and dimensions
	1309	nc_stat = NF90_INQ_DIMID( id, "Pressure", id_dim )
	1310	CALL handle_netcdf_error( 'netcdf', 550 )
	1311	nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, len = np)
	1312	CALL handle_netcdf_error( 'netcdf', 550 )
	1313
	1314	nc_stat = NF90_INQ_DIMID( id, "Absorber", id_dim )
	1315	CALL handle_netcdf_error( 'netcdf', 550 )
	1316	nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, len = nabs )
	1317	CALL handle_netcdf_error( 'netcdf', 550 )
	1318
	1319
	1320	!
	1321	!-- Allocate pressure, and trace gas arrays
	1322	ALLOCATE( p_mls(1:np) )
	1323	ALLOCATE( trace_mls(1:num_trace_gases,1:np) )
	1324	ALLOCATE( trace_mls_tmp(1:nabs,1:np) )
	1325
	1326
	1327	nc_stat = NF90_INQ_VARID( id, "Pressure", id_var )
	1328	CALL handle_netcdf_error( 'netcdf', 550 )
	1329	nc_stat = NF90_GET_VAR( id, id_var, p_mls )
	1330	CALL handle_netcdf_error( 'netcdf', 550 )
	1331
	1332	nc_stat = NF90_INQ_VARID( id, "AbsorberAmountMLS", id_var )
	1333	CALL handle_netcdf_error( 'netcdf', 550 )
	1334	nc_stat = NF90_GET_VAR( id, id_var, trace_mls_tmp )
	1335	CALL handle_netcdf_error( 'netcdf', 550 )
	1336
	1337
	1338	!
	1339	!-- Write absorber amounts (mls) to trace_mls
	1340	DO n = 1, num_trace_gases
	1341	CALL getAbsorberIndex( TRIM( trace_names(n) ), id_abs )
	1342
	1343	trace_mls(n,1:np) = trace_mls_tmp(id_abs,1:np)
	1344
	1345	!
	1346	!-- Replace missing values by zero
	1347	WHERE ( trace_mls(n,:) > 2.0_wp )
	1348	trace_mls(n,:) = 0.0_wp
	1349	END WHERE
	1350	END DO
	1351
	1352	DEALLOCATE ( trace_mls_tmp )
	1353
	1354	nc_stat = NF90_CLOSE( id )
	1355	CALL handle_netcdf_error( 'netcdf', 551 )
	1356
	1357	!
	1358	!-- Add extra pressure level for calculations of the trace gas paths
	1359	ALLOCATE ( rrtm_play_tmp(1:nzt_rad+1) )
	1360	ALLOCATE ( rrtm_plev_tmp(1:nzt_rad+2) )
	1361
	1362	rrtm_play_tmp(1:nzt_rad) = rrtm_play(0,1:nzt_rad)
	1363	rrtm_plev_tmp(1:nzt_rad+1) = rrtm_plev(0,1:nzt_rad+1)
	1364	rrtm_play_tmp(nzt_rad+1) = rrtm_plev(0,nzt_rad+1) * 0.5_wp
	1365	rrtm_plev_tmp(nzt_rad+2) = MIN( 1.0E-4_wp, 0.25_wp &
	1366	* rrtm_plev(0,nzt_rad+1) )
	1367
	1368	!
	1369	!-- Calculate trace gas path (zero at surface) with interpolation to the
	1370	!-- sounding levels
	1371	ALLOCATE ( trace_mls_path(1:nzt_rad+2,1:num_trace_gases) )
	1372
	1373	trace_mls_path(nzb+1,:) = 0.0_wp
	1374
	1375	DO k = nzb+2, nzt_rad+2
	1376	DO m = 1, num_trace_gases
	1377	trace_mls_path(k,m) = trace_mls_path(k-1,m)
	1378
	1379	!
	1380	!-- When the pressure level is higher than the trace gas pressure
	1381	!-- level, assume that
	1382	IF ( rrtm_plev_tmp(k-1) .GT. p_mls(1) ) THEN
	1383
	1384	trace_mls_path(k,m) = trace_mls_path(k,m) + trace_mls(m,1) &
	1385	* ( rrtm_plev_tmp(k-1) &
	1386	- MAX( p_mls(1), rrtm_plev_tmp(k) ) &
	1387	) / g
	1388	ENDIF
	1389
	1390	!
	1391	!-- Integrate for each sounding level from the contributing p_mls
	1392	!-- levels
	1393	DO n = 2, np
	1394	!
	1395	!-- Limit p_mls so that it is within the model level
	1396	p_mls_u = MIN( rrtm_plev_tmp(k-1), &
	1397	MAX( rrtm_plev_tmp(k), p_mls(n) ) )
	1398	p_mls_l = MIN( rrtm_plev_tmp(k-1), &
	1399	MAX( rrtm_plev_tmp(k), p_mls(n-1) ) )
	1400
	1401	IF ( p_mls_l .GT. p_mls_u ) THEN
	1402
	1403	!
	1404	!-- Calculate weights for interpolation
	1405	p_mls_m = 0.5_wp * (p_mls_l + p_mls_u)
	1406	p_wgt_u = (p_mls(n-1) - p_mls_m) / (p_mls(n-1) - p_mls(n))
	1407	p_wgt_l = (p_mls_m - p_mls(n)) / (p_mls(n-1) - p_mls(n))
	1408
	1409	!
	1410	!-- Add level to trace gas path
	1411	trace_mls_path(k,m) = trace_mls_path(k,m) &
	1412	+ ( p_wgt_u * trace_mls(m,n) &
	1413	+ p_wgt_l * trace_mls(m,n-1) ) &
	1414	* (p_mls_l + p_mls_u) / g
	1415	ENDIF
	1416	ENDDO
	1417
	1418	IF ( rrtm_plev_tmp(k) .LT. p_mls(np) ) THEN
	1419	trace_mls_path(k,m) = trace_mls_path(k,m) + trace_mls(m,np) &
	1420	* ( MIN( rrtm_plev_tmp(k-1), p_mls(np) ) &
	1421	- rrtm_plev_tmp(k) &
	1422	) / g
	1423	ENDIF
[1496]	1424	ENDDO
	1425	ENDDO
	1426
	1427
[1585]	1428	!
	1429	!-- Prepare trace gas path profiles
	1430	ALLOCATE ( trace_path_tmp(1:nzt_rad+1) )
[1496]	1431
[1585]	1432	DO m = 1, num_trace_gases
	1433
	1434	trace_path_tmp(1:nzt_rad+1) = ( trace_mls_path(2:nzt_rad+2,m) &
	1435	- trace_mls_path(1:nzt_rad+1,m) ) * g &
	1436	/ ( rrtm_plev_tmp(1:nzt_rad+1) &
	1437	- rrtm_plev_tmp(2:nzt_rad+2) )
	1438
	1439	!
	1440	!-- Save trace gas paths to the respective arrays
	1441	SELECT CASE ( TRIM( trace_names(m) ) )
	1442
	1443	CASE ( 'O3' )
	1444
	1445	rrtm_o3vmr(0,:) = trace_path_tmp(:)
	1446
	1447	CASE ( 'CO2' )
	1448
	1449	rrtm_co2vmr(0,:) = trace_path_tmp(:)
	1450
	1451	CASE ( 'CH4' )
	1452
	1453	rrtm_ch4vmr(0,:) = trace_path_tmp(:)
	1454
	1455	CASE ( 'N2O' )
	1456
	1457	rrtm_n2ovmr(0,:) = trace_path_tmp(:)
	1458
	1459	CASE ( 'O2' )
	1460
	1461	rrtm_o2vmr(0,:) = trace_path_tmp(:)
	1462
	1463	CASE ( 'CFC11' )
	1464
	1465	rrtm_cfc11vmr(0,:) = trace_path_tmp(:)
	1466
	1467	CASE ( 'CFC12' )
	1468
	1469	rrtm_cfc12vmr(0,:) = trace_path_tmp(:)
	1470
	1471	CASE ( 'CFC22' )
	1472
	1473	rrtm_cfc22vmr(0,:) = trace_path_tmp(:)
	1474
	1475	CASE ( 'CCL4' )
	1476
	1477	rrtm_ccl4vmr(0,:) = trace_path_tmp(:)
	1478
	1479	CASE DEFAULT
	1480
	1481	END SELECT
	1482
	1483	ENDDO
	1484
	1485	DEALLOCATE ( trace_path_tmp )
	1486	DEALLOCATE ( trace_mls_path )
	1487	DEALLOCATE ( rrtm_play_tmp )
	1488	DEALLOCATE ( rrtm_plev_tmp )
	1489	DEALLOCATE ( trace_mls )
	1490	DEALLOCATE ( p_mls )
	1491
	1492	END SUBROUTINE read_trace_gas_data
	1493
	1494	#endif
	1495
	1496
[1551]	1497	!------------------------------------------------------------------------------!
	1498	! Description:
	1499	! ------------
[1682]	1500	!> Calculate temperature tendency due to radiative cooling/heating.
	1501	!> Cache-optimized version.
[1551]	1502	!------------------------------------------------------------------------------!
[1585]	1503	SUBROUTINE radiation_tendency_ij ( i, j, tend )
[1496]	1504
[1585]	1505	USE arrays_3d, &
	1506	ONLY: dzw
[1496]	1507
[1585]	1508	USE cloud_parameters, &
	1509	ONLY: pt_d_t, cp
[1551]	1510
[1585]	1511	USE control_parameters, &
	1512	ONLY: rho_surface
	1513
	1514	IMPLICIT NONE
	1515
	1516	INTEGER(iwp) :: i, j, k
	1517
	1518	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend
	1519
	1520	#if defined ( __rrtmg )
	1521
	1522	REAL(wp) :: rad_sw_net_l, rad_sw_net_u, rad_lw_net_l, rad_lw_net_u
	1523
	1524	rad_sw_net_l = 0.0_wp
	1525	rad_sw_net_u = 0.0_wp
	1526	rad_lw_net_l = 0.0_wp
	1527	rad_lw_net_u = 0.0_wp
	1528
	1529	!
	1530	!-- Calculate radiative flux divergence
	1531	DO k = nzb+1, nzt+1
	1532
	1533	rad_sw_net_l = rad_sw_in(k-1,j,i) - rad_sw_out(k-1,j,i)
	1534	rad_sw_net_u = rad_sw_in(k,j,i) - rad_sw_out(k,j,i)
	1535	rad_lw_net_l = rad_lw_in(k-1,j,i) - rad_lw_out(k-1,j,i)
	1536	rad_lw_net_u = rad_lw_in(k,j,i) - rad_lw_out(k,j,i)
	1537
	1538	tend(k,j,i) = tend(k,j,i) + ( rad_sw_net_u - rad_sw_net_l &
	1539	+ rad_lw_net_u - rad_lw_net_l ) / &
	1540	( dzw(k) * cp * rho_surface ) * pt_d_t(k)
	1541	ENDDO
	1542
	1543	#endif
	1544
	1545	END SUBROUTINE radiation_tendency_ij
	1546
	1547
[1551]	1548	!------------------------------------------------------------------------------!
	1549	! Description:
	1550	! ------------
[1682]	1551	!> Calculate temperature tendency due to radiative cooling/heating.
	1552	!> Vector-optimized version
[1551]	1553	!------------------------------------------------------------------------------!
[1585]	1554	SUBROUTINE radiation_tendency ( tend )
[1551]	1555
[1585]	1556	USE arrays_3d, &
	1557	ONLY: dzw
[1551]	1558
[1585]	1559	USE cloud_parameters, &
	1560	ONLY: pt_d_t, cp
[1551]	1561
[1585]	1562	USE indices, &
	1563	ONLY: nxl, nxr, nyn, nys
	1564
	1565	USE control_parameters, &
	1566	ONLY: rho_surface
	1567
	1568	IMPLICIT NONE
	1569
	1570	INTEGER(iwp) :: i, j, k
	1571
	1572	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend
	1573
	1574	#if defined ( __rrtmg )
	1575
	1576	REAL(wp) :: rad_sw_net_l, rad_sw_net_u, rad_lw_net_l, rad_lw_net_u
	1577
	1578	rad_sw_net_l = 0.0_wp
	1579	rad_sw_net_u = 0.0_wp
	1580	rad_lw_net_l = 0.0_wp
	1581	rad_lw_net_u = 0.0_wp
	1582
	1583	DO i = nxl, nxr
	1584	DO j = nys, nyn
	1585
	1586	!
	1587	!-- Calculate radiative flux divergence
	1588	DO k = nzb+1, nzt+1
	1589
	1590	rad_sw_net_l = rad_sw_in(k-1,j,i) - rad_sw_out(k-1,j,i)
	1591	rad_sw_net_u = rad_sw_in(k ,j,i) - rad_sw_out(k ,j,i)
	1592	rad_lw_net_l = rad_lw_in(k-1,j,i) - rad_lw_out(k-1,j,i)
	1593	rad_lw_net_u = rad_lw_in(k ,j,i) - rad_lw_out(k ,j,i)
	1594
	1595	tend(k,j,i) = tend(k,j,i) + ( rad_sw_net_u - rad_sw_net_l &
	1596	+ rad_lw_net_u - rad_lw_net_l ) / &
	1597	( dzw(k) * cp * rho_surface ) * pt_d_t(k)
	1598	ENDDO
	1599	ENDDO
	1600	ENDDO
	1601	#endif
	1602
	1603	END SUBROUTINE radiation_tendency
	1604
[1496]	1605	END MODULE radiation_model_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |