Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |