Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |