Home

Context Navigation

source: palm/trunk/SOURCE/radiation_model_mod.f90 @ 1854

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |