Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |