Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |