Home

Context Navigation

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

Last change on this file since 1587 was 1587, checked in by maronga, 10 years ago
small adjustments for RRTMG
Property svn:keywords set to `Id`
File size: 63.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |