Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |