Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |