Home

Context Navigation

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

Last change on this file since 2246 was 2242, checked in by maronga, 7 years ago
revised number of soil layers. added support for RRTMG runs with dry atmosphere/soil
Property svn:keywords set to `Id`
File size: 144.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
6	! terms of the GNU General Public License as published by the Free Software
7	! Foundation, either version 3 of the License, or (at your option) any later
8	! version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 1997-2017 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: radiation_model_mod.f90 2242 2017-06-01 14:12:31Z sward $
27	! Allow for RRTMG runs without humidity/cloud physics
28	!
29	! 2233 2017-05-30 18:08:54Z suehring
30	!
31	! 2232 2017-05-30 17:47:52Z suehring
32	! Adjustments to new topography concept
33	! Bugfix in read restart
34	!
35	! 2200 2017-04-11 11:37:51Z suehring
36	! Bugfix in call of exchange_horiz_2d and read restart data
37	!
38	! 2163 2017-03-01 13:23:15Z schwenkel
39	! Bugfix in radiation_check_data_output
40	!
41	! 2157 2017-02-22 15:10:35Z suehring
42	! Bugfix in read_restart data
43	!
44	! 2011 2016-09-19 17:29:57Z kanani
45	! Removed CALL of auxiliary SUBROUTINE get_usm_info,
46	! flag urban_surface is now defined in module control_parameters.
47	!
48	! 2007 2016-08-24 15:47:17Z kanani
49	! Added calculation of solar directional vector for new urban surface
50	! model,
51	! accounted for urban_surface model in radiation_check_parameters,
52	! correction of comments for zenith angle.
53	!
54	! 2000 2016-08-20 18:09:15Z knoop
55	! Forced header and separation lines into 80 columns
56	!
57	! 1976 2016-07-27 13:28:04Z maronga
58	! Output of 2D/3D/masked data is now directly done within this module. The
59	! radiation schemes have been simplified for better usability so that
60	! rad_lw_in, rad_lw_out, rad_sw_in, and rad_sw_out are available independent of
61	! the radiation code used.
62	!
63	! 1856 2016-04-13 12:56:17Z maronga
64	! Bugfix: allocation of rad_lw_out for radiation_scheme = 'clear-sky'
65	!
66	! 1853 2016-04-11 09:00:35Z maronga
67	! Added routine for radiation_scheme = constant.
68	!
69	! 1849 2016-04-08 11:33:18Z hoffmann
70	! Adapted for modularization of microphysics
71	!
72	! 1826 2016-04-07 12:01:39Z maronga
73	! Further modularization.
74	!
75	! 1788 2016-03-10 11:01:04Z maronga
76	! Added new albedo class for pavements / roads.
77	!
78	! 1783 2016-03-06 18:36:17Z raasch
79	! palm-netcdf-module removed in order to avoid a circular module dependency,
80	! netcdf-variables moved to netcdf-module, new routine netcdf_handle_error_rad
81	! added
82	!
83	! 1757 2016-02-22 15:49:32Z maronga
84	! Added parameter unscheduled_radiation_calls. Bugfix: interpolation of sounding
85	! profiles for pressure and temperature above the LES domain.
86	!
87	! 1709 2015-11-04 14:47:01Z maronga
88	! Bugfix: set initial value for rrtm_lwuflx_dt to zero, small formatting
89	! corrections
90	!
91	! 1701 2015-11-02 07:43:04Z maronga
92	! Bugfixes: wrong index for output of timeseries, setting of nz_snd_end
93	!
94	! 1691 2015-10-26 16:17:44Z maronga
95	! Added option for spin-up runs without radiation (skip_time_do_radiation). Bugfix
96	! in calculation of pressure profiles. Bugfix in calculation of trace gas profiles.
97	! Added output of radiative heating rates.
98	!
99	! 1682 2015-10-07 23:56:08Z knoop
100	! Code annotations made doxygen readable
101	!
102	! 1606 2015-06-29 10:43:37Z maronga
103	! Added preprocessor directive __netcdf to allow for compiling without netCDF.
104	! Note, however, that RRTMG cannot be used without netCDF.
105	!
106	! 1590 2015-05-08 13:56:27Z maronga
107	! Bugfix: definition of character strings requires same length for all elements
108	!
109	! 1587 2015-05-04 14:19:01Z maronga
110	! Added albedo class for snow
111	!
112	! 1585 2015-04-30 07:05:52Z maronga
113	! Added support for RRTMG
114	!
115	! 1571 2015-03-12 16:12:49Z maronga
116	! Added missing KIND attribute. Removed upper-case variable names
117	!
118	! 1551 2015-03-03 14:18:16Z maronga
119	! Added support for data output. Various variables have been renamed. Added
120	! interface for different radiation schemes (currently: clear-sky, constant, and
121	! RRTM (not yet implemented).
122	!
123	! 1496 2014-12-02 17:25:50Z maronga
124	! Initial revision
125	!
126	!
127	! Description:
128	! ------------
129	!> Radiation models and interfaces
130	!> @todo move variable definitions used in radiation_init only to the subroutine
131	!> as they are no longer required after initialization.
132	!> @todo Output of full column vertical profiles used in RRTMG
133	!> @todo Output of other rrtm arrays (such as volume mixing ratios)
134	!> @todo Adapt for use with topography
135	!>
136	!> @note Many variables have a leading dummy dimension (0:0) in order to
137	!> match the assume-size shape expected by the RRTMG model.
138	!------------------------------------------------------------------------------!
139	MODULE radiation_model_mod
140
141	USE arrays_3d, &
142	ONLY: dzw, hyp, pt, q, ql, zu, zw
143
144	USE cloud_parameters, &
145	ONLY: cp, l_d_cp, rho_l
146
147	USE constants, &
148	ONLY: pi
149
150	USE control_parameters, &
151	ONLY: cloud_droplets, cloud_physics, g, initializing_actions, &
152	large_scale_forcing, lsf_surf, phi, pt_surface, rho_surface, &
153	surface_pressure, time_since_reference_point
154
155	USE indices, &
156	ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt, &
157	wall_flags_0
158
159	USE kinds
160
161	USE microphysics_mod, &
162	ONLY: nc_const, sigma_gc
163
164	#if defined ( __netcdf )
165	USE NETCDF
166	#endif
167
168	#if defined ( __rrtmg )
169	USE parrrsw, &
170	ONLY: naerec, nbndsw
171
172	USE parrrtm, &
173	ONLY: nbndlw
174
175	USE rrtmg_lw_init, &
176	ONLY: rrtmg_lw_ini
177
178	USE rrtmg_sw_init, &
179	ONLY: rrtmg_sw_ini
180
181	USE rrtmg_lw_rad, &
182	ONLY: rrtmg_lw
183
184	USE rrtmg_sw_rad, &
185	ONLY: rrtmg_sw
186	#endif
187
188
189
190	IMPLICIT NONE
191
192	CHARACTER(10) :: radiation_scheme = 'clear-sky' ! 'constant', 'clear-sky', or 'rrtmg'
193
194	!
195	!-- Predefined Land surface classes (albedo_type) after Briegleb (1992)
196	CHARACTER(37), DIMENSION(0:17), PARAMETER :: albedo_type_name = (/ &
197	'user defined ', & ! 0
198	'ocean ', & ! 1
199	'mixed farming, tall grassland ', & ! 2
200	'tall/medium grassland ', & ! 3
201	'evergreen shrubland ', & ! 4
202	'short grassland/meadow/shrubland ', & ! 5
203	'evergreen needleleaf forest ', & ! 6
204	'mixed deciduous evergreen forest ', & ! 7
205	'deciduous forest ', & ! 8
206	'tropical evergreen broadleaved forest', & ! 9
207	'medium/tall grassland/woodland ', & ! 10
208	'desert, sandy ', & ! 11
209	'desert, rocky ', & ! 12
210	'tundra ', & ! 13
211	'land ice ', & ! 14
212	'sea ice ', & ! 15
213	'snow ', & ! 16
214	'pavement/roads ' & ! 17
215	/)
216
217	INTEGER(iwp) :: albedo_type = 5, & !< Albedo surface type (default: short grassland)
218	day, & !< current day of the year
219	day_init = 172, & !< day of the year at model start (21/06)
220	dots_rad = 0 !< starting index for timeseries output
221
222	LOGICAL :: unscheduled_radiation_calls = .TRUE., & !< flag parameter indicating whether additional calls of the radiation code are allowed
223	constant_albedo = .FALSE., & !< flag parameter indicating whether the albedo may change depending on zenith
224	force_radiation_call = .FALSE., & !< flag parameter for unscheduled radiation calls
225	lw_radiation = .TRUE., & !< flag parameter indicating whether longwave radiation shall be calculated
226	radiation = .FALSE., & !< flag parameter indicating whether the radiation model is used
227	sun_up = .TRUE., & !< flag parameter indicating whether the sun is up or down
228	sw_radiation = .TRUE., & !< flag parameter indicing whether shortwave radiation shall be calculated
229	sun_direction = .FALSE. !< flag parameter indicing whether solar direction shall be calculated
230
231
232	REAL(wp), PARAMETER :: d_seconds_hour = 0.000277777777778_wp, & !< inverse of seconds per hour (1/3600)
233	d_hours_day = 0.0416666666667_wp, & !< inverse of hours per day (1/24)
234	sigma_sb = 5.67037321E-8_wp, & !< Stefan-Boltzmann constant
235	solar_constant = 1368.0_wp !< solar constant at top of atmosphere
236
237	REAL(wp) :: albedo = 9999999.9_wp, & !< NAMELIST alpha
238	albedo_lw_dif = 9999999.9_wp, & !< NAMELIST aldif
239	albedo_lw_dir = 9999999.9_wp, & !< NAMELIST aldir
240	albedo_sw_dif = 9999999.9_wp, & !< NAMELIST asdif
241	albedo_sw_dir = 9999999.9_wp, & !< NAMELIST asdir
242	decl_1, & !< declination coef. 1
243	decl_2, & !< declination coef. 2
244	decl_3, & !< declination coef. 3
245	dt_radiation = 0.0_wp, & !< radiation model timestep
246	emissivity = 0.98_wp, & !< NAMELIST surface emissivity
247	lambda = 0.0_wp, & !< longitude in degrees
248	lon = 0.0_wp, & !< longitude in radians
249	lat = 0.0_wp, & !< latitude in radians
250	net_radiation = 0.0_wp, & !< net radiation at surface
251	skip_time_do_radiation = 0.0_wp, & !< Radiation model is not called before this time
252	sky_trans, & !< sky transmissivity
253	time_radiation = 0.0_wp, & !< time since last call of radiation code
254	time_utc, & !< current time in UTC
255	time_utc_init = 43200.0_wp !< UTC time at model start (noon)
256
257	REAL(wp), DIMENSION(0:0) :: zenith, & !< cosine of solar zenith angle
258	sun_dir_lat, & !< solar directional vector in latitudes
259	sun_dir_lon !< solar directional vector in longitudes
260
261	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
262	alpha, & !< surface broadband albedo (used for clear-sky scheme)
263	rad_lw_out_change_0, & !< change in LW out due to change in surface temperature
264	rad_net, & !< net radiation at the surface
265	rad_net_av !< average of rad_net
266
267	!
268	!-- Land surface albedos for solar zenith angle of 60Â° after Briegleb (1992)
269	!-- (shortwave, longwave, broadband): sw, lw, bb,
270	REAL(wp), DIMENSION(0:2,1:17), PARAMETER :: albedo_pars = RESHAPE( (/&
271	0.06_wp, 0.06_wp, 0.06_wp, & ! 1
272	0.09_wp, 0.28_wp, 0.19_wp, & ! 2
273	0.11_wp, 0.33_wp, 0.23_wp, & ! 3
274	0.11_wp, 0.33_wp, 0.23_wp, & ! 4
275	0.14_wp, 0.34_wp, 0.25_wp, & ! 5
276	0.06_wp, 0.22_wp, 0.14_wp, & ! 6
277	0.06_wp, 0.27_wp, 0.17_wp, & ! 7
278	0.06_wp, 0.31_wp, 0.19_wp, & ! 8
279	0.06_wp, 0.22_wp, 0.14_wp, & ! 9
280	0.06_wp, 0.28_wp, 0.18_wp, & ! 10
281	0.35_wp, 0.51_wp, 0.43_wp, & ! 11
282	0.24_wp, 0.40_wp, 0.32_wp, & ! 12
283	0.10_wp, 0.27_wp, 0.19_wp, & ! 13
284	0.90_wp, 0.65_wp, 0.77_wp, & ! 14
285	0.90_wp, 0.65_wp, 0.77_wp, & ! 15
286	0.95_wp, 0.70_wp, 0.82_wp, & ! 16
287	0.08_wp, 0.08_wp, 0.08_wp & ! 17
288	/), (/ 3, 17 /) )
289
290	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: &
291	rad_lw_cs_hr, & !< longwave clear sky radiation heating rate (K/s)
292	rad_lw_cs_hr_av, & !< average of rad_lw_cs_hr
293	rad_lw_hr, & !< longwave radiation heating rate (K/s)
294	rad_lw_hr_av, & !< average of rad_sw_hr
295	rad_lw_in, & !< incoming longwave radiation (W/m2)
296	rad_lw_in_av, & !< average of rad_lw_in
297	rad_lw_out, & !< outgoing longwave radiation (W/m2)
298	rad_lw_out_av, & !< average of rad_lw_out
299	rad_sw_cs_hr, & !< shortwave clear sky radiation heating rate (K/s)
300	rad_sw_cs_hr_av, & !< average of rad_sw_cs_hr
301	rad_sw_hr, & !< shortwave radiation heating rate (K/s)
302	rad_sw_hr_av, & !< average of rad_sw_hr
303	rad_sw_in, & !< incoming shortwave radiation (W/m2)
304	rad_sw_in_av, & !< average of rad_sw_in
305	rad_sw_out, & !< outgoing shortwave radiation (W/m2)
306	rad_sw_out_av !< average of rad_sw_out
307
308
309	!
310	!-- Variables and parameters used in RRTMG only
311	#if defined ( __rrtmg )
312	CHARACTER(LEN=12) :: rrtm_input_file = "RAD_SND_DATA" !< name of the NetCDF input file (sounding data)
313
314
315	!
316	!-- Flag parameters for RRTMGS (should not be changed)
317	INTEGER(iwp), PARAMETER :: rrtm_idrv = 1, & !< flag for longwave upward flux calculation option (0,1)
318	rrtm_inflglw = 2, & !< flag for lw cloud optical properties (0,1,2)
319	rrtm_iceflglw = 0, & !< flag for lw ice particle specifications (0,1,2,3)
320	rrtm_liqflglw = 1, & !< flag for lw liquid droplet specifications
321	rrtm_inflgsw = 2, & !< flag for sw cloud optical properties (0,1,2)
322	rrtm_iceflgsw = 0, & !< flag for sw ice particle specifications (0,1,2,3)
323	rrtm_liqflgsw = 1 !< flag for sw liquid droplet specifications
324
325	!
326	!-- The following variables should be only changed with care, as this will
327	!-- require further setting of some variables, which is currently not
328	!-- implemented (aerosols, ice phase).
329	INTEGER(iwp) :: nzt_rad, & !< upper vertical limit for radiation calculations
330	rrtm_icld = 0, & !< cloud flag (0: clear sky column, 1: cloudy column)
331	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)
332
333	INTEGER(iwp) :: nc_stat !< local variable for storin the result of netCDF calls for error message handling
334
335	LOGICAL :: snd_exists = .FALSE. !< flag parameter to check whether a user-defined input files exists
336
337	REAL(wp), PARAMETER :: mol_mass_air_d_wv = 1.607793_wp !< molecular weight dry air / water vapor
338
339	REAL(wp), DIMENSION(:), ALLOCATABLE :: hyp_snd, & !< hypostatic pressure from sounding data (hPa)
340	q_snd, & !< specific humidity from sounding data (kg/kg) - dummy at the moment
341	rrtm_tsfc, & !< dummy array for storing surface temperature
342	t_snd !< actual temperature from sounding data (hPa)
343
344	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: aldif, & !< longwave diffuse albedo solar angle of 60Â°
345	aldir, & !< longwave direct albedo solar angle of 60Â°
346	asdif, & !< shortwave diffuse albedo solar angle of 60Â°
347	asdir, & !< shortwave direct albedo solar angle of 60Â°
348	rrtm_ccl4vmr, & !< CCL4 volume mixing ratio (g/mol)
349	rrtm_cfc11vmr, & !< CFC11 volume mixing ratio (g/mol)
350	rrtm_cfc12vmr, & !< CFC12 volume mixing ratio (g/mol)
351	rrtm_cfc22vmr, & !< CFC22 volume mixing ratio (g/mol)
352	rrtm_ch4vmr, & !< CH4 volume mixing ratio
353	rrtm_cicewp, & !< in-cloud ice water path (g/mÂ²)
354	rrtm_cldfr, & !< cloud fraction (0,1)
355	rrtm_cliqwp, & !< in-cloud liquid water path (g/mÂ²)
356	rrtm_co2vmr, & !< CO2 volume mixing ratio (g/mol)
357	rrtm_emis, & !< surface emissivity (0-1)
358	rrtm_h2ovmr, & !< H2O volume mixing ratio
359	rrtm_n2ovmr, & !< N2O volume mixing ratio
360	rrtm_o2vmr, & !< O2 volume mixing ratio
361	rrtm_o3vmr, & !< O3 volume mixing ratio
362	rrtm_play, & !< pressure layers (hPa, zu-grid)
363	rrtm_plev, & !< pressure layers (hPa, zw-grid)
364	rrtm_reice, & !< cloud ice effective radius (microns)
365	rrtm_reliq, & !< cloud water drop effective radius (microns)
366	rrtm_tlay, & !< actual temperature (K, zu-grid)
367	rrtm_tlev, & !< actual temperature (K, zw-grid)
368	rrtm_lwdflx, & !< RRTM output of incoming longwave radiation flux (W/m2)
369	rrtm_lwdflxc, & !< RRTM output of outgoing clear sky longwave radiation flux (W/m2)
370	rrtm_lwuflx, & !< RRTM output of outgoing longwave radiation flux (W/m2)
371	rrtm_lwuflxc, & !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
372	rrtm_lwuflx_dt, & !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
373	rrtm_lwuflxc_dt,& !< RRTM output of outgoing clear sky longwave radiation flux (W/m2)
374	rrtm_lwhr, & !< RRTM output of longwave radiation heating rate (K/d)
375	rrtm_lwhrc, & !< RRTM output of incoming longwave clear sky radiation heating rate (K/d)
376	rrtm_swdflx, & !< RRTM output of incoming shortwave radiation flux (W/m2)
377	rrtm_swdflxc, & !< RRTM output of outgoing clear sky shortwave radiation flux (W/m2)
378	rrtm_swuflx, & !< RRTM output of outgoing shortwave radiation flux (W/m2)
379	rrtm_swuflxc, & !< RRTM output of incoming clear sky shortwave radiation flux (W/m2)
380	rrtm_swhr, & !< RRTM output of shortwave radiation heating rate (K/d)
381	rrtm_swhrc !< RRTM output of incoming shortwave clear sky radiation heating rate (K/d)
382
383	!
384	!-- Definition of arrays that are currently not used for calling RRTMG (due to setting of flag parameters)
385	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: rad_lw_cs_in, & !< incoming clear sky longwave radiation (W/m2) (not used)
386	rad_lw_cs_out, & !< outgoing clear sky longwave radiation (W/m2) (not used)
387	rad_sw_cs_in, & !< incoming clear sky shortwave radiation (W/m2) (not used)
388	rad_sw_cs_out, & !< outgoing clear sky shortwave radiation (W/m2) (not used)
389	rrtm_aldif, & !< surface albedo for longwave diffuse radiation
390	rrtm_aldir, & !< surface albedo for longwave direct radiation
391	rrtm_asdif, & !< surface albedo for shortwave diffuse radiation
392	rrtm_asdir, & !< surface albedo for shortwave direct radiation
393	rrtm_lw_tauaer, & !< lw aerosol optical depth
394	rrtm_lw_taucld, & !< lw in-cloud optical depth
395	rrtm_sw_taucld, & !< sw in-cloud optical depth
396	rrtm_sw_ssacld, & !< sw in-cloud single scattering albedo
397	rrtm_sw_asmcld, & !< sw in-cloud asymmetry parameter
398	rrtm_sw_fsfcld, & !< sw in-cloud forward scattering fraction
399	rrtm_sw_tauaer, & !< sw aerosol optical depth
400	rrtm_sw_ssaaer, & !< sw aerosol single scattering albedo
401	rrtm_sw_asmaer, & !< sw aerosol asymmetry parameter
402	rrtm_sw_ecaer !< sw aerosol optical detph at 0.55 microns (rrtm_iaer = 6 only)
403
404	#endif
405
406	INTERFACE radiation_check_data_output
407	MODULE PROCEDURE radiation_check_data_output
408	END INTERFACE radiation_check_data_output
409
410	INTERFACE radiation_check_data_output_pr
411	MODULE PROCEDURE radiation_check_data_output_pr
412	END INTERFACE radiation_check_data_output_pr
413
414	INTERFACE radiation_check_parameters
415	MODULE PROCEDURE radiation_check_parameters
416	END INTERFACE radiation_check_parameters
417
418	INTERFACE radiation_clearsky
419	MODULE PROCEDURE radiation_clearsky
420	END INTERFACE radiation_clearsky
421
422	INTERFACE radiation_constant
423	MODULE PROCEDURE radiation_constant
424	END INTERFACE radiation_constant
425
426	INTERFACE radiation_control
427	MODULE PROCEDURE radiation_control
428	END INTERFACE radiation_control
429
430	INTERFACE radiation_3d_data_averaging
431	MODULE PROCEDURE radiation_3d_data_averaging
432	END INTERFACE radiation_3d_data_averaging
433
434	INTERFACE radiation_data_output_2d
435	MODULE PROCEDURE radiation_data_output_2d
436	END INTERFACE radiation_data_output_2d
437
438	INTERFACE radiation_data_output_3d
439	MODULE PROCEDURE radiation_data_output_3d
440	END INTERFACE radiation_data_output_3d
441
442	INTERFACE radiation_data_output_mask
443	MODULE PROCEDURE radiation_data_output_mask
444	END INTERFACE radiation_data_output_mask
445
446	INTERFACE radiation_define_netcdf_grid
447	MODULE PROCEDURE radiation_define_netcdf_grid
448	END INTERFACE radiation_define_netcdf_grid
449
450	INTERFACE radiation_header
451	MODULE PROCEDURE radiation_header
452	END INTERFACE radiation_header
453
454	INTERFACE radiation_init
455	MODULE PROCEDURE radiation_init
456	END INTERFACE radiation_init
457
458	INTERFACE radiation_parin
459	MODULE PROCEDURE radiation_parin
460	END INTERFACE radiation_parin
461
462	INTERFACE radiation_rrtmg
463	MODULE PROCEDURE radiation_rrtmg
464	END INTERFACE radiation_rrtmg
465
466	INTERFACE radiation_tendency
467	MODULE PROCEDURE radiation_tendency
468	MODULE PROCEDURE radiation_tendency_ij
469	END INTERFACE radiation_tendency
470
471	INTERFACE radiation_read_restart_data
472	MODULE PROCEDURE radiation_read_restart_data
473	END INTERFACE radiation_read_restart_data
474
475	INTERFACE radiation_last_actions
476	MODULE PROCEDURE radiation_last_actions
477	END INTERFACE radiation_last_actions
478
479	SAVE
480
481	PRIVATE
482
483	!
484	!-- Public functions / NEEDS SORTING
485	PUBLIC radiation_check_data_output, radiation_check_data_output_pr, &
486	radiation_check_parameters, radiation_control, &
487	radiation_header, radiation_init, radiation_parin, &
488	radiation_3d_data_averaging, radiation_tendency, &
489	radiation_data_output_2d, radiation_data_output_3d, &
490	radiation_define_netcdf_grid, radiation_last_actions, &
491	radiation_read_restart_data, radiation_data_output_mask
492
493	!
494	!-- Public variables and constants / NEEDS SORTING
495	PUBLIC dots_rad, dt_radiation, force_radiation_call, &
496	rad_net, rad_net_av, radiation, radiation_scheme, rad_lw_in, &
497	rad_lw_in_av, rad_lw_out, rad_lw_out_av, rad_lw_out_change_0, &
498	rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av, rad_sw_in, &
499	rad_sw_in_av, rad_sw_out, rad_sw_out_av, rad_sw_cs_hr, &
500	rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av, sigma_sb, &
501	skip_time_do_radiation, time_radiation, unscheduled_radiation_calls,&
502	zenith, calc_zenith, sun_direction, sun_dir_lat, sun_dir_lon, &
503	day_init, time_utc_init
504
505
506	#if defined ( __rrtmg )
507	PUBLIC rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir
508	#endif
509
510	CONTAINS
511
512
513	!------------------------------------------------------------------------------!
514	! Description:
515	! ------------
516	!> This subroutine controls the calls of the radiation schemes
517	!------------------------------------------------------------------------------!
518	SUBROUTINE radiation_control
519
520
521	IMPLICIT NONE
522
523
524	SELECT CASE ( TRIM( radiation_scheme ) )
525
526	CASE ( 'constant' )
527	CALL radiation_constant
528
529	CASE ( 'clear-sky' )
530	CALL radiation_clearsky
531
532	CASE ( 'rrtmg' )
533	CALL radiation_rrtmg
534
535	CASE DEFAULT
536
537	END SELECT
538
539
540	END SUBROUTINE radiation_control
541
542	!------------------------------------------------------------------------------!
543	! Description:
544	! ------------
545	!> Check data output for radiation model
546	!------------------------------------------------------------------------------!
547	SUBROUTINE radiation_check_data_output( var, unit, i, ilen, k )
548
549
550	USE control_parameters, &
551	ONLY: data_output, message_string
552
553	IMPLICIT NONE
554
555	CHARACTER (LEN=*) :: unit !<
556	CHARACTER (LEN=*) :: var !<
557
558	INTEGER(iwp) :: i
559	INTEGER(iwp) :: ilen
560	INTEGER(iwp) :: k
561
562	SELECT CASE ( TRIM( var ) )
563
564	CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_sw_cs_hr', 'rad_sw_hr' )
565	IF ( .NOT. radiation .OR. radiation_scheme /= 'rrtmg' ) THEN
566	message_string = '"output of "' // TRIM( var ) // '" requi' // &
567	'res radiation = .TRUE. and ' // &
568	'radiation_scheme = "rrtmg"'
569	CALL message( 'check_parameters', 'PA0406', 1, 2, 0, 6, 0 )
570	ENDIF
571	unit = 'K/h'
572
573	CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_sw_in', 'rad_sw_out' )
574	IF ( .NOT. radiation .OR. radiation_scheme /= 'rrtmg' ) THEN
575	message_string = '"output of "' // TRIM( var ) // '" requi' // &
576	'res radiation = .TRUE. and ' // &
577	'radiation_scheme = "rrtmg"'
578	CALL message( 'check_parameters', 'PA0406', 1, 2, 0, 6, 0 )
579	ENDIF
580	unit = 'W/m2'
581
582	CASE ( 'rad_net', 'rrtm_aldif', 'rrtm_aldir', 'rrtm_asdif', &
583	'rrtm_asdir*' )
584	IF ( k == 0 .OR. data_output(i)(ilen-2:ilen) /= '_xy' ) THEN
585	message_string = 'illegal value for data_output: "' // &
586	TRIM( var ) // '" & only 2d-horizontal ' // &
587	'cross sections are allowed for this value'
588	CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
589	ENDIF
590	IF ( .NOT. radiation .OR. radiation_scheme /= "rrtmg" ) THEN
591	IF ( TRIM( var ) == 'rrtm_aldif*' .OR. &
592	TRIM( var ) == 'rrtm_aldir*' .OR. &
593	TRIM( var ) == 'rrtm_asdif*' .OR. &
594	TRIM( var ) == 'rrtm_asdir*' ) &
595	THEN
596	message_string = 'output of "' // TRIM( var ) // '" require'&
597	// 's radiation = .TRUE. and radiation_sch'&
598	// 'eme = "rrtmg"'
599	CALL message( 'check_parameters', 'PA0409', 1, 2, 0, 6, 0 )
600	ENDIF
601	ENDIF
602
603	IF ( TRIM( var ) == 'rad_net*' ) unit = 'W/m2'
604	IF ( TRIM( var ) == 'rrtm_aldif*' ) unit = ''
605	IF ( TRIM( var ) == 'rrtm_aldir*' ) unit = ''
606	IF ( TRIM( var ) == 'rrtm_asdif*' ) unit = ''
607	IF ( TRIM( var ) == 'rrtm_asdir*' ) unit = ''
608
609	CASE DEFAULT
610	unit = 'illegal'
611
612	END SELECT
613
614
615	END SUBROUTINE radiation_check_data_output
616
617	!------------------------------------------------------------------------------!
618	! Description:
619	! ------------
620	!> Check data output of profiles for radiation model
621	!------------------------------------------------------------------------------!
622	SUBROUTINE radiation_check_data_output_pr( variable, var_count, unit, dopr_unit )
623
624	USE arrays_3d, &
625	ONLY: zu
626
627	USE control_parameters, &
628	ONLY: data_output_pr, message_string
629
630	USE indices
631
632	USE profil_parameter
633
634	USE statistics
635
636	IMPLICIT NONE
637
638	CHARACTER (LEN=*) :: unit !<
639	CHARACTER (LEN=*) :: variable !<
640	CHARACTER (LEN=*) :: dopr_unit !< local value of dopr_unit
641
642	INTEGER(iwp) :: user_pr_index !<
643	INTEGER(iwp) :: var_count !<
644
645	SELECT CASE ( TRIM( variable ) )
646
647	CASE ( 'rad_net' )
648	IF ( ( .NOT. radiation ) .OR. radiation_scheme == 'constant' )&
649	THEN
650	message_string = 'data_output_pr = ' // &
651	TRIM( data_output_pr(var_count) ) // ' is' // &
652	'not available for radiation = .FALSE. or ' //&
653	'radiation_scheme = "constant"'
654	CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
655	ELSE
656	dopr_index(var_count) = 101
657	dopr_unit = 'W/m2'
658	hom(:,2,101,:) = SPREAD( zw, 2, statistic_regions+1 )
659	unit = dopr_unit
660	ENDIF
661
662	CASE ( 'rad_lw_in' )
663	IF ( ( .NOT. radiation) .OR. radiation_scheme == 'constant' ) &
664	THEN
665	message_string = 'data_output_pr = ' // &
666	TRIM( data_output_pr(var_count) ) // ' is' // &
667	'not available for radiation = .FALSE. or ' //&
668	'radiation_scheme = "constant"'
669	CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
670	ELSE
671	dopr_index(var_count) = 102
672	dopr_unit = 'W/m2'
673	hom(:,2,102,:) = SPREAD( zw, 2, statistic_regions+1 )
674	unit = dopr_unit
675	ENDIF
676
677	CASE ( 'rad_lw_out' )
678	IF ( ( .NOT. radiation ) .OR. radiation_scheme == 'constant' ) &
679	THEN
680	message_string = 'data_output_pr = ' // &
681	TRIM( data_output_pr(var_count) ) // ' is' // &
682	'not available for radiation = .FALSE. or ' //&
683	'radiation_scheme = "constant"'
684	CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
685	ELSE
686	dopr_index(var_count) = 103
687	dopr_unit = 'W/m2'
688	hom(:,2,103,:) = SPREAD( zw, 2, statistic_regions+1 )
689	unit = dopr_unit
690	ENDIF
691
692	CASE ( 'rad_sw_in' )
693	IF ( ( .NOT. radiation ) .OR. radiation_scheme == 'constant' ) &
694	THEN
695	message_string = 'data_output_pr = ' // &
696	TRIM( data_output_pr(var_count) ) // ' is' // &
697	'not available for radiation = .FALSE. or ' //&
698	'radiation_scheme = "constant"'
699	CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
700	ELSE
701	dopr_index(var_count) = 104
702	dopr_unit = 'W/m2'
703	hom(:,2,104,:) = SPREAD( zw, 2, statistic_regions+1 )
704	unit = dopr_unit
705	ENDIF
706
707	CASE ( 'rad_sw_out')
708	IF ( ( .NOT. radiation ) .OR. radiation_scheme == 'constant' )&
709	THEN
710	message_string = 'data_output_pr = ' // &
711	TRIM( data_output_pr(var_count) ) // ' is' // &
712	'not available for radiation = .FALSE. or ' //&
713	'radiation_scheme = "constant"'
714	CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
715	ELSE
716	dopr_index(var_count) = 105
717	dopr_unit = 'W/m2'
718	hom(:,2,105,:) = SPREAD( zw, 2, statistic_regions+1 )
719	unit = dopr_unit
720	ENDIF
721
722	CASE ( 'rad_lw_cs_hr' )
723	IF ( ( .NOT. radiation ) .OR. radiation_scheme /= 'rrtmg' ) &
724	THEN
725	message_string = 'data_output_pr = ' // &
726	TRIM( data_output_pr(var_count) ) // ' is' // &
727	'not available for radiation = .FALSE. or ' //&
728	'radiation_scheme /= "rrtmg"'
729	CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
730	ELSE
731	dopr_index(var_count) = 106
732	dopr_unit = 'K/h'
733	hom(:,2,106,:) = SPREAD( zu, 2, statistic_regions+1 )
734	unit = dopr_unit
735	ENDIF
736
737	CASE ( 'rad_lw_hr' )
738	IF ( ( .NOT. radiation ) .OR. radiation_scheme /= 'rrtmg' ) &
739	THEN
740	message_string = 'data_output_pr = ' // &
741	TRIM( data_output_pr(var_count) ) // ' is' // &
742	'not available for radiation = .FALSE. or ' //&
743	'radiation_scheme /= "rrtmg"'
744	CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
745	ELSE
746	dopr_index(var_count) = 107
747	dopr_unit = 'K/h'
748	hom(:,2,107,:) = SPREAD( zu, 2, statistic_regions+1 )
749	unit = dopr_unit
750	ENDIF
751
752	CASE ( 'rad_sw_cs_hr' )
753	IF ( ( .NOT. radiation ) .OR. radiation_scheme /= 'rrtmg' ) &
754	THEN
755	message_string = 'data_output_pr = ' // &
756	TRIM( data_output_pr(var_count) ) // ' is' // &
757	'not available for radiation = .FALSE. or ' //&
758	'radiation_scheme /= "rrtmg"'
759	CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
760	ELSE
761	dopr_index(var_count) = 108
762	dopr_unit = 'K/h'
763	hom(:,2,108,:) = SPREAD( zu, 2, statistic_regions+1 )
764	unit = dopr_unit
765	ENDIF
766
767	CASE ( 'rad_sw_hr' )
768	IF ( ( .NOT. radiation ) .OR. radiation_scheme /= 'rrtmg' ) &
769	THEN
770	message_string = 'data_output_pr = ' // &
771	TRIM( data_output_pr(var_count) ) // ' is' // &
772	'not available for radiation = .FALSE. or ' //&
773	'radiation_scheme /= "rrtmg"'
774	CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
775	ELSE
776	dopr_index(var_count) = 109
777	dopr_unit = 'K/h'
778	hom(:,2,109,:) = SPREAD( zu, 2, statistic_regions+1 )
779	unit = dopr_unit
780	ENDIF
781
782
783	CASE DEFAULT
784	unit = 'illegal'
785
786	END SELECT
787
788
789	END SUBROUTINE radiation_check_data_output_pr
790
791
792	!------------------------------------------------------------------------------!
793	! Description:
794	! ------------
795	!> Check parameters routine for radiation model
796	!------------------------------------------------------------------------------!
797	SUBROUTINE radiation_check_parameters
798
799	USE control_parameters, &
800	ONLY: message_string, topography, urban_surface
801
802
803	IMPLICIT NONE
804
805
806	IF ( radiation_scheme /= 'constant' .AND. &
807	radiation_scheme /= 'clear-sky' .AND. &
808	radiation_scheme /= 'rrtmg' ) THEN
809	message_string = 'unknown radiation_scheme = '// &
810	TRIM( radiation_scheme )
811	CALL message( 'check_parameters', 'PA0405', 1, 2, 0, 6, 0 )
812	ELSEIF ( radiation_scheme == 'rrtmg' ) THEN
813	#if ! defined ( __rrtmg )
814	message_string = 'radiation_scheme = "rrtmg" requires ' // &
815	'compilation of PALM with pre-processor ' // &
816	'directive -D__rrtmg'
817	CALL message( 'check_parameters', 'PA0407', 1, 2, 0, 6, 0 )
818	#endif
819	#if defined ( __rrtmg ) && ! defined( __netcdf )
820	message_string = 'radiation_scheme = "rrtmg" requires ' // &
821	'the use of NetCDF (preprocessor directive ' // &
822	'-D__netcdf'
823	CALL message( 'check_parameters', 'PA0412', 1, 2, 0, 6, 0 )
824	#endif
825
826	ENDIF
827
828	IF ( albedo_type == 0 .AND. albedo == 9999999.9_wp .AND. &
829	radiation_scheme == 'clear-sky') THEN
830	message_string = 'radiation_scheme = "clear-sky" in combination' // &
831	'with albedo_type = 0 requires setting of albedo'// &
832	' /= 9999999.9'
833	CALL message( 'check_parameters', 'PA0410', 1, 2, 0, 6, 0 )
834	ENDIF
835
836	IF ( albedo_type == 0 .AND. radiation_scheme == 'rrtmg' .AND. &
837	( albedo_lw_dif == 9999999.9_wp .OR. albedo_lw_dir == 9999999.9_wp&
838	.OR. albedo_sw_dif == 9999999.9_wp .OR. albedo_sw_dir == 9999999.9_wp&
839	) ) THEN
840	message_string = 'radiation_scheme = "rrtmg" in combination' // &
841	'with albedo_type = 0 requires setting of ' // &
842	'albedo_lw_dif /= 9999999.9' // &
843	'albedo_lw_dir /= 9999999.9' // &
844	'albedo_sw_dif /= 9999999.9 and' // &
845	'albedo_sw_dir /= 9999999.9'
846	CALL message( 'check_parameters', 'PA0411', 1, 2, 0, 6, 0 )
847	ENDIF
848
849	!
850	!-- The following paramter check is temporarily extended by the urban_surface
851	!-- flag, until a better solution comes up to omit this check in case of
852	!-- urban surface model is used.
853	IF ( topography /= 'flat' .AND. .NOT. urban_surface ) THEN
854	message_string = 'radiation scheme cannot be used ' // &
855	'in combination with topography /= "flat"'
856	CALL message( 'check_parameters', 'PA0414', 1, 2, 0, 6, 0 )
857	ENDIF
858
859	END SUBROUTINE radiation_check_parameters
860
861
862	!------------------------------------------------------------------------------!
863	! Description:
864	! ------------
865	!> Initialization of the radiation model
866	!------------------------------------------------------------------------------!
867	SUBROUTINE radiation_init
868
869	IMPLICIT NONE
870
871	!
872	!-- Allocate array for storing the surface net radiation
873	IF ( .NOT. ALLOCATED ( rad_net ) ) THEN
874	ALLOCATE ( rad_net(nysg:nyng,nxlg:nxrg) )
875	rad_net = 0.0_wp
876	ENDIF
877
878	!
879	!-- Allocate array for storing the surface net radiation
880	IF ( .NOT. ALLOCATED ( rad_lw_out_change_0 ) ) THEN
881	ALLOCATE ( rad_lw_out_change_0(nysg:nyng,nxlg:nxrg) )
882	rad_lw_out_change_0 = 0.0_wp
883	ENDIF
884
885	!
886	!-- Fix net radiation in case of radiation_scheme = 'constant'
887	IF ( radiation_scheme == 'constant' ) THEN
888	rad_net = net_radiation
889	! radiation = .FALSE.
890	!
891	!-- Calculate orbital constants
892	ELSE
893	decl_1 = SIN(23.45_wp * pi / 180.0_wp)
894	decl_2 = 2.0_wp * pi / 365.0_wp
895	decl_3 = decl_2 * 81.0_wp
896	lat = phi * pi / 180.0_wp
897	lon = lambda * pi / 180.0_wp
898	ENDIF
899
900	IF ( radiation_scheme == 'clear-sky' .OR. &
901	radiation_scheme == 'constant') THEN
902
903	ALLOCATE ( alpha(nysg:nyng,nxlg:nxrg) )
904
905	IF ( .NOT. ALLOCATED ( rad_sw_in ) ) THEN
906	ALLOCATE ( rad_sw_in(0:0,nysg:nyng,nxlg:nxrg) )
907	ENDIF
908	IF ( .NOT. ALLOCATED ( rad_sw_out ) ) THEN
909	ALLOCATE ( rad_sw_out(0:0,nysg:nyng,nxlg:nxrg) )
910	ENDIF
911
912	IF ( .NOT. ALLOCATED ( rad_sw_in_av ) ) THEN
913	ALLOCATE ( rad_sw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
914	ENDIF
915	IF ( .NOT. ALLOCATED ( rad_sw_out_av ) ) THEN
916	ALLOCATE ( rad_sw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
917	ENDIF
918
919	IF ( .NOT. ALLOCATED ( rad_lw_in ) ) THEN
920	ALLOCATE ( rad_lw_in(0:0,nysg:nyng,nxlg:nxrg) )
921	ENDIF
922	IF ( .NOT. ALLOCATED ( rad_lw_out ) ) THEN
923	ALLOCATE ( rad_lw_out(0:0,nysg:nyng,nxlg:nxrg) )
924	ENDIF
925
926	IF ( .NOT. ALLOCATED ( rad_lw_in_av ) ) THEN
927	ALLOCATE ( rad_lw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
928	ENDIF
929	IF ( .NOT. ALLOCATED ( rad_lw_out_av ) ) THEN
930	ALLOCATE ( rad_lw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
931	ENDIF
932
933	rad_sw_in = 0.0_wp
934	rad_sw_out = 0.0_wp
935	rad_lw_in = 0.0_wp
936	rad_lw_out = 0.0_wp
937
938	!
939	!-- Overwrite albedo if manually set in parameter file
940	IF ( albedo_type /= 0 .AND. albedo == 9999999.9_wp ) THEN
941	albedo = albedo_pars(2,albedo_type)
942	ENDIF
943
944	alpha = albedo
945
946	!
947	!-- Initialization actions for RRTMG
948	ELSEIF ( radiation_scheme == 'rrtmg' ) THEN
949	#if defined ( __rrtmg )
950	!
951	!-- Allocate albedos
952	ALLOCATE ( rrtm_aldif(0:0,nysg:nyng,nxlg:nxrg) )
953	ALLOCATE ( rrtm_aldir(0:0,nysg:nyng,nxlg:nxrg) )
954	ALLOCATE ( rrtm_asdif(0:0,nysg:nyng,nxlg:nxrg) )
955	ALLOCATE ( rrtm_asdir(0:0,nysg:nyng,nxlg:nxrg) )
956	ALLOCATE ( aldif(nysg:nyng,nxlg:nxrg) )
957	ALLOCATE ( aldir(nysg:nyng,nxlg:nxrg) )
958	ALLOCATE ( asdif(nysg:nyng,nxlg:nxrg) )
959	ALLOCATE ( asdir(nysg:nyng,nxlg:nxrg) )
960
961	IF ( albedo_type /= 0 ) THEN
962	IF ( albedo_lw_dif == 9999999.9_wp ) THEN
963	albedo_lw_dif = albedo_pars(0,albedo_type)
964	albedo_lw_dir = albedo_lw_dif
965	ENDIF
966	IF ( albedo_sw_dif == 9999999.9_wp ) THEN
967	albedo_sw_dif = albedo_pars(1,albedo_type)
968	albedo_sw_dir = albedo_sw_dif
969	ENDIF
970	ENDIF
971
972	aldif(:,:) = albedo_lw_dif
973	aldir(:,:) = albedo_lw_dir
974	asdif(:,:) = albedo_sw_dif
975	asdir(:,:) = albedo_sw_dir
976	!
977	!-- Calculate initial values of current (cosine of) the zenith angle and
978	!-- whether the sun is up
979	CALL calc_zenith
980	!
981	!-- Calculate initial surface albedo
982	IF ( .NOT. constant_albedo ) THEN
983	CALL calc_albedo
984	ELSE
985	rrtm_aldif(0,:,:) = aldif(:,:)
986	rrtm_aldir(0,:,:) = aldir(:,:)
987	rrtm_asdif(0,:,:) = asdif(:,:)
988	rrtm_asdir(0,:,:) = asdir(:,:)
989	ENDIF
990
991	!
992	!-- Allocate surface emissivity
993	ALLOCATE ( rrtm_emis(0:0,1:nbndlw+1) )
994	rrtm_emis = emissivity
995
996	!
997	!-- Allocate 3d arrays of radiative fluxes and heating rates
998	IF ( .NOT. ALLOCATED ( rad_sw_in ) ) THEN
999	ALLOCATE ( rad_sw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1000	rad_sw_in = 0.0_wp
1001	ENDIF
1002
1003	IF ( .NOT. ALLOCATED ( rad_sw_in_av ) ) THEN
1004	ALLOCATE ( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1005	ENDIF
1006
1007	IF ( .NOT. ALLOCATED ( rad_sw_out ) ) THEN
1008	ALLOCATE ( rad_sw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1009	rad_sw_out = 0.0_wp
1010	ENDIF
1011
1012	IF ( .NOT. ALLOCATED ( rad_sw_out_av ) ) THEN
1013	ALLOCATE ( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1014	ENDIF
1015
1016	IF ( .NOT. ALLOCATED ( rad_sw_hr ) ) THEN
1017	ALLOCATE ( rad_sw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1018	rad_sw_hr = 0.0_wp
1019	ENDIF
1020
1021	IF ( .NOT. ALLOCATED ( rad_sw_hr_av ) ) THEN
1022	ALLOCATE ( rad_sw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1023	rad_sw_hr_av = 0.0_wp
1024	ENDIF
1025
1026	IF ( .NOT. ALLOCATED ( rad_sw_cs_hr ) ) THEN
1027	ALLOCATE ( rad_sw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1028	rad_sw_cs_hr = 0.0_wp
1029	ENDIF
1030
1031	IF ( .NOT. ALLOCATED ( rad_sw_cs_hr_av ) ) THEN
1032	ALLOCATE ( rad_sw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1033	rad_sw_cs_hr_av = 0.0_wp
1034	ENDIF
1035
1036	IF ( .NOT. ALLOCATED ( rad_lw_in ) ) THEN
1037	ALLOCATE ( rad_lw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1038	rad_lw_in = 0.0_wp
1039	ENDIF
1040
1041	IF ( .NOT. ALLOCATED ( rad_lw_in_av ) ) THEN
1042	ALLOCATE ( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1043	ENDIF
1044
1045	IF ( .NOT. ALLOCATED ( rad_lw_out ) ) THEN
1046	ALLOCATE ( rad_lw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1047	rad_lw_out = 0.0_wp
1048	ENDIF
1049
1050	IF ( .NOT. ALLOCATED ( rad_lw_out_av ) ) THEN
1051	ALLOCATE ( rad_lw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1052	ENDIF
1053
1054	IF ( .NOT. ALLOCATED ( rad_lw_hr ) ) THEN
1055	ALLOCATE ( rad_lw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1056	rad_lw_hr = 0.0_wp
1057	ENDIF
1058
1059	IF ( .NOT. ALLOCATED ( rad_lw_hr_av ) ) THEN
1060	ALLOCATE ( rad_lw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1061	rad_lw_hr_av = 0.0_wp
1062	ENDIF
1063
1064	IF ( .NOT. ALLOCATED ( rad_lw_cs_hr ) ) THEN
1065	ALLOCATE ( rad_lw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1066	rad_lw_cs_hr = 0.0_wp
1067	ENDIF
1068
1069	IF ( .NOT. ALLOCATED ( rad_lw_cs_hr_av ) ) THEN
1070	ALLOCATE ( rad_lw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1071	rad_lw_cs_hr_av = 0.0_wp
1072	ENDIF
1073
1074	ALLOCATE ( rad_sw_cs_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1075	ALLOCATE ( rad_sw_cs_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1076	rad_sw_cs_in = 0.0_wp
1077	rad_sw_cs_out = 0.0_wp
1078
1079	ALLOCATE ( rad_lw_cs_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1080	ALLOCATE ( rad_lw_cs_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1081	rad_lw_cs_in = 0.0_wp
1082	rad_lw_cs_out = 0.0_wp
1083
1084	!
1085	!-- Allocate dummy array for storing surface temperature
1086	ALLOCATE ( rrtm_tsfc(1) )
1087
1088	!
1089	!-- Initialize RRTMG
1090	IF ( lw_radiation ) CALL rrtmg_lw_ini ( cp )
1091	IF ( sw_radiation ) CALL rrtmg_sw_ini ( cp )
1092
1093	!
1094	!-- Set input files for RRTMG
1095	INQUIRE(FILE="RAD_SND_DATA", EXIST=snd_exists)
1096	IF ( .NOT. snd_exists ) THEN
1097	rrtm_input_file = "rrtmg_lw.nc"
1098	ENDIF
1099
1100	!
1101	!-- Read vertical layers for RRTMG from sounding data
1102	!-- The routine provides nzt_rad, hyp_snd(1:nzt_rad),
1103	!-- t_snd(nzt+2:nzt_rad), rrtm_play(1:nzt_rad), rrtm_plev(1_nzt_rad+1),
1104	!-- rrtm_tlay(nzt+2:nzt_rad), rrtm_tlev(nzt+2:nzt_rad+1)
1105	CALL read_sounding_data
1106
1107	!
1108	!-- Read trace gas profiles from file. This routine provides
1109	!-- the rrtm_ arrays (1:nzt_rad+1)
1110	CALL read_trace_gas_data
1111	#endif
1112	ENDIF
1113
1114	!
1115	!-- Perform user actions if required
1116	CALL user_init_radiation
1117
1118	!
1119	!-- Calculate radiative fluxes at model start
1120	IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
1121
1122	SELECT CASE ( radiation_scheme )
1123	CASE ( 'rrtmg' )
1124	CALL radiation_rrtmg
1125	CASE ( 'clear-sky' )
1126	CALL radiation_clearsky
1127	CASE ( 'constant' )
1128	CALL radiation_constant
1129	CASE DEFAULT
1130	END SELECT
1131
1132	ENDIF
1133
1134	RETURN
1135
1136	END SUBROUTINE radiation_init
1137
1138
1139	!------------------------------------------------------------------------------!
1140	! Description:
1141	! ------------
1142	!> A simple clear sky radiation model
1143	!------------------------------------------------------------------------------!
1144	SUBROUTINE radiation_clearsky
1145
1146
1147	IMPLICIT NONE
1148
1149	INTEGER(iwp) :: i, j, k !< loop indices
1150	REAL(wp) :: exn, & !< Exner functions at surface
1151	exn1, & !< Exner functions at first grid level
1152	pt1 !< potential temperature at first grid level
1153
1154	!
1155	!-- Calculate current zenith angle
1156	CALL calc_zenith
1157
1158	!
1159	!-- Calculate sky transmissivity
1160	sky_trans = 0.6_wp + 0.2_wp * zenith(0)
1161
1162	!
1163	!-- Calculate value of the Exner function
1164	exn = (surface_pressure / 1000.0_wp )**0.286_wp
1165	!
1166	!-- Calculate radiation fluxes and net radiation (rad_net) for each grid
1167	!-- point
1168	DO i = nxlg, nxrg
1169	DO j = nysg, nyng
1170	!
1171	!-- Obtain vertical index of topography top
1172	k = MAXLOC( &
1173	MERGE( 1, 0, &
1174	BTEST( wall_flags_0(:,j,i), 12 ) &
1175	), DIM = 1 &
1176	) - 1
1177
1178	exn1 = (hyp(k+1) / 100000.0_wp )**0.286_wp
1179
1180	rad_sw_in(0,j,i) = solar_constant * sky_trans * zenith(0)
1181	rad_sw_out(0,j,i) = alpha(j,i) * rad_sw_in(0,j,i)
1182	rad_lw_out(0,j,i) = emissivity * sigma_sb * (pt(k,j,i) * exn)**4
1183
1184	IF ( cloud_physics ) THEN
1185	pt1 = pt(k+1,j,i) + l_d_cp / exn1 * ql(k+1,j,i)
1186	rad_lw_in(0,j,i) = 0.8_wp * sigma_sb * (pt1 * exn1)**4
1187	ELSE
1188	rad_lw_in(0,j,i) = 0.8_wp * sigma_sb * (pt(k+1,j,i) * exn1)**4
1189	ENDIF
1190
1191	rad_net(j,i) = rad_sw_in(0,j,i) - rad_sw_out(0,j,i) &
1192	+ rad_lw_in(0,j,i) - rad_lw_out(0,j,i)
1193
1194
1195	rad_lw_out_change_0(j,i) = 3.0_wp * sigma_sb * emissivity &
1196	* (pt(k,j,i) * exn) ** 3
1197
1198	ENDDO
1199	ENDDO
1200
1201	END SUBROUTINE radiation_clearsky
1202
1203
1204	!------------------------------------------------------------------------------!
1205	! Description:
1206	! ------------
1207	!> This scheme keeps the prescribed net radiation constant during the run
1208	!------------------------------------------------------------------------------!
1209	SUBROUTINE radiation_constant
1210
1211
1212	IMPLICIT NONE
1213
1214	INTEGER(iwp) :: i, j, k !< loop indices
1215	REAL(wp) :: exn, & !< Exner functions at surface
1216	exn1, & !< Exner functions at first grid level
1217	pt1 !< potential temperature at first grid level
1218
1219	!
1220	!-- Calculate value of the Exner function
1221	exn = (surface_pressure / 1000.0_wp )**0.286_wp
1222	!
1223	!-- Prescribe net radiation and estimate the remaining radiative fluxes
1224	DO i = nxlg, nxrg
1225	DO j = nysg, nyng
1226	!
1227	!-- Obtain vertical index of topography top. So far it is identical to nzb.
1228	k = MAXLOC( &
1229	MERGE( 1, 0, &
1230	BTEST( wall_flags_0(:,j,i), 12 ) &
1231	), DIM = 1 &
1232	) - 1
1233
1234	rad_net(j,i) = net_radiation
1235
1236	exn1 = (hyp(k+1) / 100000.0_wp )**0.286_wp
1237
1238	IF ( cloud_physics ) THEN
1239	pt1 = pt(k+1,j,i) + l_d_cp / exn1 * ql(k+1,j,i)
1240	rad_lw_in(0,j,i) = 0.8_wp * sigma_sb * (pt1 * exn1)**4
1241	ELSE
1242	rad_lw_in(0,j,i) = 0.8_wp * sigma_sb * (pt(k+1,j,i) * exn1)**4
1243	ENDIF
1244
1245	rad_lw_out(0,j,i) = emissivity * sigma_sb * (pt(k,j,i) * exn)**4
1246
1247	rad_sw_in(0,j,i) = ( rad_net(j,i) - rad_lw_in(0,j,i) &
1248	+ rad_lw_out(0,j,i) ) &
1249	/ ( 1.0_wp - alpha(j,i) )
1250
1251	ENDDO
1252	ENDDO
1253
1254	END SUBROUTINE radiation_constant
1255
1256	!------------------------------------------------------------------------------!
1257	! Description:
1258	! ------------
1259	!> Header output for radiation model
1260	!------------------------------------------------------------------------------!
1261	SUBROUTINE radiation_header ( io )
1262
1263
1264	IMPLICIT NONE
1265
1266	INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file
1267
1268
1269
1270	!
1271	!-- Write radiation model header
1272	WRITE( io, 3 )
1273
1274	IF ( radiation_scheme == "constant" ) THEN
1275	WRITE( io, 4 ) net_radiation
1276	ELSEIF ( radiation_scheme == "clear-sky" ) THEN
1277	WRITE( io, 5 )
1278	ELSEIF ( radiation_scheme == "rrtmg" ) THEN
1279	WRITE( io, 6 )
1280	IF ( .NOT. lw_radiation ) WRITE( io, 10 )
1281	IF ( .NOT. sw_radiation ) WRITE( io, 11 )
1282	ENDIF
1283
1284	IF ( albedo_type == 0 ) THEN
1285	WRITE( io, 7 ) albedo
1286	ELSE
1287	WRITE( io, 8 ) TRIM( albedo_type_name(albedo_type) )
1288	ENDIF
1289	IF ( constant_albedo ) THEN
1290	WRITE( io, 9 )
1291	ENDIF
1292
1293	IF ( radiation .AND. radiation_scheme /= 'constant' ) THEN
1294	WRITE ( io, 1 ) lambda
1295	WRITE ( io, 2 ) day_init, time_utc_init
1296	ENDIF
1297
1298	WRITE( io, 12 ) dt_radiation
1299
1300
1301	1 FORMAT (' Geograph. longitude : lambda = ',F4.1,' degr')
1302	2 FORMAT (' Day of the year at model start : day_init = ',I3 &
1303	/' UTC time at model start : time_utc_init = ',F7.1' s')
1304	3 FORMAT (//' Radiation model information:'/ &
1305	' ----------------------------'/)
1306	4 FORMAT (' --> Using constant net radiation: net_radiation = ', F6.2, &
1307	// 'W/m**2')
1308	5 FORMAT (' --> Simple radiation scheme for clear sky is used (no clouds,', &
1309	' default)')
1310	6 FORMAT (' --> RRTMG scheme is used')
1311	7 FORMAT (/' User-specific surface albedo: albedo =', F6.3)
1312	8 FORMAT (/' Albedo is set for land surface type: ', A)
1313	9 FORMAT (/' --> Albedo is fixed during the run')
1314	10 FORMAT (/' --> Longwave radiation is disabled')
1315	11 FORMAT (/' --> Shortwave radiation is disabled.')
1316	12 FORMAT (' Timestep: dt_radiation = ', F6.2, ' s')
1317
1318
1319	END SUBROUTINE radiation_header
1320
1321
1322	!------------------------------------------------------------------------------!
1323	! Description:
1324	! ------------
1325	!> Parin for &radiation_par for radiation model
1326	!------------------------------------------------------------------------------!
1327	SUBROUTINE radiation_parin
1328
1329
1330	IMPLICIT NONE
1331
1332	CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file
1333
1334	NAMELIST /radiation_par/ albedo, albedo_type, albedo_lw_dir, &
1335	albedo_lw_dif, albedo_sw_dir, albedo_sw_dif, &
1336	constant_albedo, day_init, dt_radiation, &
1337	lambda, lw_radiation, net_radiation, &
1338	radiation_scheme, skip_time_do_radiation, &
1339	sw_radiation, time_utc_init, &
1340	unscheduled_radiation_calls
1341
1342	line = ' '
1343
1344	!
1345	!-- Try to find radiation model package
1346	REWIND ( 11 )
1347	line = ' '
1348	DO WHILE ( INDEX( line, '&radiation_par' ) == 0 )
1349	READ ( 11, '(A)', END=10 ) line
1350	ENDDO
1351	BACKSPACE ( 11 )
1352
1353	!
1354	!-- Read user-defined namelist
1355	READ ( 11, radiation_par )
1356
1357	!
1358	!-- Set flag that indicates that the radiation model is switched on
1359	radiation = .TRUE.
1360
1361	10 CONTINUE
1362
1363
1364	END SUBROUTINE radiation_parin
1365
1366
1367	!------------------------------------------------------------------------------!
1368	! Description:
1369	! ------------
1370	!> Implementation of the RRTMG radiation_scheme
1371	!------------------------------------------------------------------------------!
1372	SUBROUTINE radiation_rrtmg
1373
1374	USE indices, &
1375	ONLY: nbgp
1376
1377	USE particle_attributes, &
1378	ONLY: grid_particles, number_of_particles, particles, &
1379	particle_advection_start, prt_count
1380
1381	IMPLICIT NONE
1382
1383	#if defined ( __rrtmg )
1384
1385	INTEGER(iwp) :: i, j, k, n !< loop indices
1386
1387	REAL(wp) :: s_r2, & !< weighted sum over all droplets with r^2
1388	s_r3 !< weighted sum over all droplets with r^3
1389
1390	!
1391	!-- Calculate current (cosine of) zenith angle and whether the sun is up
1392	CALL calc_zenith
1393	!
1394	!-- Calculate surface albedo
1395	IF ( .NOT. constant_albedo ) THEN
1396	CALL calc_albedo
1397	ENDIF
1398
1399	!
1400	!-- Prepare input data for RRTMG
1401
1402	!
1403	!-- In case of large scale forcing with surface data, calculate new pressure
1404	!-- profile. nzt_rad might be modified by these calls and all required arrays
1405	!-- will then be re-allocated
1406	IF ( large_scale_forcing .AND. lsf_surf ) THEN
1407	CALL read_sounding_data
1408	CALL read_trace_gas_data
1409	ENDIF
1410	!
1411	!-- Loop over all grid points
1412	DO i = nxl, nxr
1413	DO j = nys, nyn
1414
1415	!
1416	!-- Prepare profiles of temperature and H2O volume mixing ratio
1417	rrtm_tlev(0,nzb+1) = pt(nzb,j,i) * ( surface_pressure &
1418	/ 1000.0_wp )**0.286_wp
1419
1420
1421	IF ( cloud_physics ) THEN
1422	DO k = nzb+1, nzt+1
1423	rrtm_tlay(0,k) = pt(k,j,i) * ( (hyp(k) ) / 100000.0_wp &
1424	)*0.286_wp + l_d_cp ql(k,j,i)
1425	rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * (q(k,j,i) - ql(k,j,i))
1426	ENDDO
1427	ELSE
1428	DO k = nzb+1, nzt+1
1429	rrtm_tlay(0,k) = pt(k,j,i) * ( (hyp(k) ) / 100000.0_wp &
1430	)**0.286_wp
1431	rrtm_h2ovmr(0,k) = 0.0_wp
1432	ENDDO
1433	ENDIF
1434
1435	!
1436	!-- Avoid temperature/humidity jumps at the top of the LES domain by
1437	!-- linear interpolation from nzt+2 to nzt+7
1438	DO k = nzt+2, nzt+7
1439	rrtm_tlay(0,k) = rrtm_tlay(0,nzt+1) &
1440	+ ( rrtm_tlay(0,nzt+8) - rrtm_tlay(0,nzt+1) ) &
1441	/ ( rrtm_play(0,nzt+8) - rrtm_play(0,nzt+1) ) &
1442	* ( rrtm_play(0,k) - rrtm_play(0,nzt+1) )
1443
1444	rrtm_h2ovmr(0,k) = rrtm_h2ovmr(0,nzt+1) &
1445	+ ( rrtm_h2ovmr(0,nzt+8) - rrtm_h2ovmr(0,nzt+1) )&
1446	/ ( rrtm_play(0,nzt+8) - rrtm_play(0,nzt+1) )&
1447	* ( rrtm_play(0,k) - rrtm_play(0,nzt+1) )
1448
1449	ENDDO
1450
1451	!-- Linear interpolate to zw grid
1452	DO k = nzb+2, nzt+8
1453	rrtm_tlev(0,k) = rrtm_tlay(0,k-1) + (rrtm_tlay(0,k) - &
1454	rrtm_tlay(0,k-1)) &
1455	/ ( rrtm_play(0,k) - rrtm_play(0,k-1) ) &
1456	* ( rrtm_plev(0,k) - rrtm_play(0,k-1) )
1457	ENDDO
1458
1459
1460	!
1461	!-- Calculate liquid water path and cloud fraction for each column.
1462	!-- Note that LWP is required in g/mÂ² instead of kg/kg m.
1463	rrtm_cldfr = 0.0_wp
1464	rrtm_reliq = 0.0_wp
1465	rrtm_cliqwp = 0.0_wp
1466	rrtm_icld = 0
1467
1468	IF ( cloud_physics ) THEN
1469	DO k = nzb+1, nzt+1
1470	rrtm_cliqwp(0,k) = ql(k,j,i) * 1000.0_wp * &
1471	(rrtm_plev(0,k) - rrtm_plev(0,k+1)) &
1472	* 100.0_wp / g
1473
1474	IF ( rrtm_cliqwp(0,k) > 0.0_wp ) THEN
1475	rrtm_cldfr(0,k) = 1.0_wp
1476	IF ( rrtm_icld == 0 ) rrtm_icld = 1
1477
1478	!
1479	!-- Calculate cloud droplet effective radius
1480	IF ( cloud_physics ) THEN
1481	rrtm_reliq(0,k) = 1.0E6_wp * ( 3.0_wp * ql(k,j,i) &
1482	* rho_surface &
1483	/ ( 4.0_wp * pi * nc_const * rho_l ) &
1484	)**0.33333333333333_wp &
1485	* EXP( LOG( sigma_gc )**2 )
1486
1487	ELSEIF ( cloud_droplets ) THEN
1488	number_of_particles = prt_count(k,j,i)
1489
1490	IF (number_of_particles <= 0) CYCLE
1491	particles => grid_particles(k,j,i)%particles(1:number_of_particles)
1492	s_r2 = 0.0_wp
1493	s_r3 = 0.0_wp
1494
1495	DO n = 1, number_of_particles
1496	IF ( particles(n)%particle_mask ) THEN
1497	s_r2 = s_r2 + particles(n)%radius*2 &
1498	particles(n)%weight_factor
1499	s_r3 = s_r3 + particles(n)%radius*3 &
1500	particles(n)%weight_factor
1501	ENDIF
1502	ENDDO
1503
1504	IF ( s_r2 > 0.0_wp ) rrtm_reliq(0,k) = s_r3 / s_r2
1505
1506	ENDIF
1507
1508	!
1509	!-- Limit effective radius
1510	IF ( rrtm_reliq(0,k) > 0.0_wp ) THEN
1511	rrtm_reliq(0,k) = MAX(rrtm_reliq(0,k),2.5_wp)
1512	rrtm_reliq(0,k) = MIN(rrtm_reliq(0,k),60.0_wp)
1513	ENDIF
1514	ENDIF
1515	ENDDO
1516	ENDIF
1517
1518	!
1519	!-- Set surface temperature
1520	rrtm_tsfc = pt(nzb,j,i) * (surface_pressure / 1000.0_wp )**0.286_wp
1521
1522	IF ( lw_radiation ) THEN
1523	CALL rrtmg_lw( 1, nzt_rad , rrtm_icld , rrtm_idrv ,&
1524	rrtm_play , rrtm_plev , rrtm_tlay , rrtm_tlev ,&
1525	rrtm_tsfc , rrtm_h2ovmr , rrtm_o3vmr , rrtm_co2vmr ,&
1526	rrtm_ch4vmr , rrtm_n2ovmr , rrtm_o2vmr , rrtm_cfc11vmr ,&
1527	rrtm_cfc12vmr , rrtm_cfc22vmr, rrtm_ccl4vmr , rrtm_emis ,&
1528	rrtm_inflglw , rrtm_iceflglw, rrtm_liqflglw, rrtm_cldfr ,&
1529	rrtm_lw_taucld , rrtm_cicewp , rrtm_cliqwp , rrtm_reice ,&
1530	rrtm_reliq , rrtm_lw_tauaer, &
1531	rrtm_lwuflx , rrtm_lwdflx , rrtm_lwhr , &
1532	rrtm_lwuflxc , rrtm_lwdflxc , rrtm_lwhrc , &
1533	rrtm_lwuflx_dt , rrtm_lwuflxc_dt )
1534
1535	!
1536	!-- Save fluxes
1537	DO k = nzb, nzt+1
1538	rad_lw_in(k,j,i) = rrtm_lwdflx(0,k)
1539	rad_lw_out(k,j,i) = rrtm_lwuflx(0,k)
1540	ENDDO
1541
1542	!
1543	!-- Save heating rates (convert from K/d to K/h)
1544	DO k = nzb+1, nzt+1
1545	rad_lw_hr(k,j,i) = rrtm_lwhr(0,k) * d_hours_day
1546	rad_lw_cs_hr(k,j,i) = rrtm_lwhrc(0,k) * d_hours_day
1547	ENDDO
1548
1549	!
1550	!-- Save change in LW heating rate
1551	rad_lw_out_change_0(j,i) = rrtm_lwuflx_dt(0,nzb)
1552
1553	ENDIF
1554
1555	IF ( sw_radiation .AND. sun_up ) THEN
1556	CALL rrtmg_sw( 1, nzt_rad , rrtm_icld , rrtm_iaer ,&
1557	rrtm_play , rrtm_plev , rrtm_tlay , rrtm_tlev ,&
1558	rrtm_tsfc , rrtm_h2ovmr , rrtm_o3vmr , rrtm_co2vmr ,&
1559	rrtm_ch4vmr , rrtm_n2ovmr , rrtm_o2vmr , rrtm_asdir(:,j,i),&
1560	rrtm_asdif(:,j,i), rrtm_aldir(:,j,i), rrtm_aldif(:,j,i), zenith,&
1561	0.0_wp , day , solar_constant, rrtm_inflgsw,&
1562	rrtm_iceflgsw , rrtm_liqflgsw, rrtm_cldfr , rrtm_sw_taucld ,&
1563	rrtm_sw_ssacld , rrtm_sw_asmcld, rrtm_sw_fsfcld, rrtm_cicewp ,&
1564	rrtm_cliqwp , rrtm_reice , rrtm_reliq , rrtm_sw_tauaer ,&
1565	rrtm_sw_ssaaer , rrtm_sw_asmaer , rrtm_sw_ecaer , &
1566	rrtm_swuflx , rrtm_swdflx , rrtm_swhr , &
1567	rrtm_swuflxc , rrtm_swdflxc , rrtm_swhrc )
1568
1569	!
1570	!-- Save fluxes
1571	DO k = nzb, nzt+1
1572	rad_sw_in(k,j,i) = rrtm_swdflx(0,k)
1573	rad_sw_out(k,j,i) = rrtm_swuflx(0,k)
1574	ENDDO
1575
1576	!
1577	!-- Save heating rates (convert from K/d to K/s)
1578	DO k = nzb+1, nzt+1
1579	rad_sw_hr(k,j,i) = rrtm_swhr(0,k) * d_hours_day
1580	rad_sw_cs_hr(k,j,i) = rrtm_swhrc(0,k) * d_hours_day
1581	ENDDO
1582
1583	ENDIF
1584
1585	!
1586	!-- Calculate surface net radiation
1587	rad_net(j,i) = rad_sw_in(nzb,j,i) - rad_sw_out(nzb,j,i) &
1588	+ rad_lw_in(nzb,j,i) - rad_lw_out(nzb,j,i)
1589
1590	ENDDO
1591	ENDDO
1592
1593	CALL exchange_horiz( rad_lw_in, nbgp )
1594	CALL exchange_horiz( rad_lw_out, nbgp )
1595	CALL exchange_horiz( rad_lw_hr, nbgp )
1596	CALL exchange_horiz( rad_lw_cs_hr, nbgp )
1597
1598	CALL exchange_horiz( rad_sw_in, nbgp )
1599	CALL exchange_horiz( rad_sw_out, nbgp )
1600	CALL exchange_horiz( rad_sw_hr, nbgp )
1601	CALL exchange_horiz( rad_sw_cs_hr, nbgp )
1602
1603	CALL exchange_horiz_2d( rad_net )
1604	CALL exchange_horiz_2d( rad_lw_out_change_0 )
1605	#endif
1606
1607	END SUBROUTINE radiation_rrtmg
1608
1609
1610	!------------------------------------------------------------------------------!
1611	! Description:
1612	! ------------
1613	!> Calculate the cosine of the zenith angle (variable is called zenith)
1614	!------------------------------------------------------------------------------!
1615	SUBROUTINE calc_zenith
1616
1617	IMPLICIT NONE
1618
1619	REAL(wp) :: declination, & !< solar declination angle
1620	hour_angle !< solar hour angle
1621	!
1622	!-- Calculate current day and time based on the initial values and simulation
1623	!-- time
1624	day = day_init + INT(FLOOR( (time_utc_init + time_since_reference_point) &
1625	/ 86400.0_wp ), KIND=iwp)
1626	time_utc = MOD((time_utc_init + time_since_reference_point), 86400.0_wp)
1627
1628
1629	!
1630	!-- Calculate solar declination and hour angle
1631	declination = ASIN( decl_1 * SIN(decl_2 * REAL(day, KIND=wp) - decl_3) )
1632	hour_angle = 2.0_wp * pi * (time_utc / 86400.0_wp) + lon - pi
1633
1634	!
1635	!-- Calculate cosine of solar zenith angle
1636	zenith(0) = SIN(lat) * SIN(declination) + COS(lat) * COS(declination) &
1637	* COS(hour_angle)
1638	zenith(0) = MAX(0.0_wp,zenith(0))
1639
1640	!
1641	!-- Calculate solar directional vector
1642	IF ( sun_direction ) THEN
1643	!-- Direction in longitudes equals to sin(solar_azimuth) * sin(zenith)
1644	sun_dir_lon(0) = -SIN(hour_angle) * COS(declination)
1645	!-- Direction in latitues equals to cos(solar_azimuth) * sin(zenith)
1646	sun_dir_lat(0) = SIN(declination) * COS(lat) - COS(hour_angle) &
1647	* COS(declination) * SIN(lat)
1648	ENDIF
1649
1650	!
1651	!-- Check if the sun is up (otheriwse shortwave calculations can be skipped)
1652	IF ( zenith(0) > 0.0_wp ) THEN
1653	sun_up = .TRUE.
1654	ELSE
1655	sun_up = .FALSE.
1656	END IF
1657
1658	END SUBROUTINE calc_zenith
1659
1660	#if defined ( __rrtmg ) && defined ( __netcdf )
1661	!------------------------------------------------------------------------------!
1662	! Description:
1663	! ------------
1664	!> Calculates surface albedo components based on Briegleb (1992) and
1665	!> Briegleb et al. (1986)
1666	!------------------------------------------------------------------------------!
1667	SUBROUTINE calc_albedo
1668
1669	IMPLICIT NONE
1670
1671	IF ( sun_up ) THEN
1672	!
1673	!-- Ocean
1674	IF ( albedo_type == 1 ) THEN
1675	rrtm_aldir(0,:,:) = 0.026_wp / ( zenith(0)**1.7_wp + 0.065_wp ) &
1676	+ 0.15_wp * ( zenith(0) - 0.1_wp ) &
1677	* ( zenith(0) - 0.5_wp ) &
1678	* ( zenith(0) - 1.0_wp )
1679	rrtm_asdir(0,:,:) = rrtm_aldir(0,:,:)
1680	!
1681	!-- Snow
1682	ELSEIF ( albedo_type == 16 ) THEN
1683	IF ( zenith(0) < 0.5_wp ) THEN
1684	rrtm_aldir(0,:,:) = 0.5_wp * (1.0_wp - aldif) &
1685	* ( 3.0_wp / (1.0_wp + 4.0_wp &
1686	* zenith(0))) - 1.0_wp
1687	rrtm_asdir(0,:,:) = 0.5_wp * (1.0_wp - asdif) &
1688	* ( 3.0_wp / (1.0_wp + 4.0_wp &
1689	* zenith(0))) - 1.0_wp
1690
1691	rrtm_aldir(0,:,:) = MIN(0.98_wp, rrtm_aldir(0,:,:))
1692	rrtm_asdir(0,:,:) = MIN(0.98_wp, rrtm_asdir(0,:,:))
1693	ELSE
1694	rrtm_aldir(0,:,:) = aldif
1695	rrtm_asdir(0,:,:) = asdif
1696	ENDIF
1697	!
1698	!-- Sea ice
1699	ELSEIF ( albedo_type == 15 ) THEN
1700	rrtm_aldir(0,:,:) = aldif
1701	rrtm_asdir(0,:,:) = asdif
1702
1703	!
1704	!-- Asphalt
1705	ELSEIF ( albedo_type == 17 ) THEN
1706	rrtm_aldir(0,:,:) = aldif
1707	rrtm_asdir(0,:,:) = asdif
1708	!
1709	!-- Land surfaces
1710	ELSE
1711	SELECT CASE ( albedo_type )
1712
1713	!
1714	!-- Surface types with strong zenith dependence
1715	CASE ( 1, 2, 3, 4, 11, 12, 13 )
1716	rrtm_aldir(0,:,:) = aldif * 1.4_wp / &
1717	(1.0_wp + 0.8_wp * zenith(0))
1718	rrtm_asdir(0,:,:) = asdif * 1.4_wp / &
1719	(1.0_wp + 0.8_wp * zenith(0))
1720	!
1721	!-- Surface types with weak zenith dependence
1722	CASE ( 5, 6, 7, 8, 9, 10, 14 )
1723	rrtm_aldir(0,:,:) = aldif * 1.1_wp / &
1724	(1.0_wp + 0.2_wp * zenith(0))
1725	rrtm_asdir(0,:,:) = asdif * 1.1_wp / &
1726	(1.0_wp + 0.2_wp * zenith(0))
1727
1728	CASE DEFAULT
1729
1730	END SELECT
1731	ENDIF
1732	!
1733	!-- Diffusive albedo is taken from Table 2
1734	rrtm_aldif(0,:,:) = aldif
1735	rrtm_asdif(0,:,:) = asdif
1736
1737	ELSE
1738
1739	rrtm_aldir(0,:,:) = 0.0_wp
1740	rrtm_asdir(0,:,:) = 0.0_wp
1741	rrtm_aldif(0,:,:) = 0.0_wp
1742	rrtm_asdif(0,:,:) = 0.0_wp
1743	ENDIF
1744	END SUBROUTINE calc_albedo
1745
1746	!------------------------------------------------------------------------------!
1747	! Description:
1748	! ------------
1749	!> Read sounding data (pressure and temperature) from RADIATION_DATA.
1750	!------------------------------------------------------------------------------!
1751	SUBROUTINE read_sounding_data
1752
1753	IMPLICIT NONE
1754
1755	INTEGER(iwp) :: id, & !< NetCDF id of input file
1756	id_dim_zrad, & !< pressure level id in the NetCDF file
1757	id_var, & !< NetCDF variable id
1758	k, & !< loop index
1759	nz_snd, & !< number of vertical levels in the sounding data
1760	nz_snd_start, & !< start vertical index for sounding data to be used
1761	nz_snd_end !< end vertical index for souding data to be used
1762
1763	REAL(wp) :: t_surface !< actual surface temperature
1764
1765	REAL(wp), DIMENSION(:), ALLOCATABLE :: hyp_snd_tmp, & !< temporary hydrostatic pressure profile (sounding)
1766	t_snd_tmp !< temporary temperature profile (sounding)
1767
1768	!
1769	!-- In case of updates, deallocate arrays first (sufficient to check one
1770	!-- array as the others are automatically allocated). This is required
1771	!-- because nzt_rad might change during the update
1772	IF ( ALLOCATED ( hyp_snd ) ) THEN
1773	DEALLOCATE( hyp_snd )
1774	DEALLOCATE( t_snd )
1775	DEALLOCATE( q_snd )
1776	DEALLOCATE ( rrtm_play )
1777	DEALLOCATE ( rrtm_plev )
1778	DEALLOCATE ( rrtm_tlay )
1779	DEALLOCATE ( rrtm_tlev )
1780
1781	DEALLOCATE ( rrtm_h2ovmr )
1782	DEALLOCATE ( rrtm_cicewp )
1783	DEALLOCATE ( rrtm_cldfr )
1784	DEALLOCATE ( rrtm_cliqwp )
1785	DEALLOCATE ( rrtm_reice )
1786	DEALLOCATE ( rrtm_reliq )
1787	DEALLOCATE ( rrtm_lw_taucld )
1788	DEALLOCATE ( rrtm_lw_tauaer )
1789
1790	DEALLOCATE ( rrtm_lwdflx )
1791	DEALLOCATE ( rrtm_lwdflxc )
1792	DEALLOCATE ( rrtm_lwuflx )
1793	DEALLOCATE ( rrtm_lwuflxc )
1794	DEALLOCATE ( rrtm_lwuflx_dt )
1795	DEALLOCATE ( rrtm_lwuflxc_dt )
1796	DEALLOCATE ( rrtm_lwhr )
1797	DEALLOCATE ( rrtm_lwhrc )
1798
1799	DEALLOCATE ( rrtm_sw_taucld )
1800	DEALLOCATE ( rrtm_sw_ssacld )
1801	DEALLOCATE ( rrtm_sw_asmcld )
1802	DEALLOCATE ( rrtm_sw_fsfcld )
1803	DEALLOCATE ( rrtm_sw_tauaer )
1804	DEALLOCATE ( rrtm_sw_ssaaer )
1805	DEALLOCATE ( rrtm_sw_asmaer )
1806	DEALLOCATE ( rrtm_sw_ecaer )
1807
1808	DEALLOCATE ( rrtm_swdflx )
1809	DEALLOCATE ( rrtm_swdflxc )
1810	DEALLOCATE ( rrtm_swuflx )
1811	DEALLOCATE ( rrtm_swuflxc )
1812	DEALLOCATE ( rrtm_swhr )
1813	DEALLOCATE ( rrtm_swhrc )
1814
1815	ENDIF
1816
1817	!
1818	!-- Open file for reading
1819	nc_stat = NF90_OPEN( rrtm_input_file, NF90_NOWRITE, id )
1820	CALL netcdf_handle_error_rad( 'read_sounding_data', 549 )
1821
1822	!
1823	!-- Inquire dimension of z axis and save in nz_snd
1824	nc_stat = NF90_INQ_DIMID( id, "Pressure", id_dim_zrad )
1825	nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim_zrad, len = nz_snd )
1826	CALL netcdf_handle_error_rad( 'read_sounding_data', 551 )
1827
1828	!
1829	! !-- Allocate temporary array for storing pressure data
1830	ALLOCATE( hyp_snd_tmp(1:nz_snd) )
1831	hyp_snd_tmp = 0.0_wp
1832
1833
1834	!-- Read pressure from file
1835	nc_stat = NF90_INQ_VARID( id, "Pressure", id_var )
1836	nc_stat = NF90_GET_VAR( id, id_var, hyp_snd_tmp(:), start = (/1/), &
1837	count = (/nz_snd/) )
1838	CALL netcdf_handle_error_rad( 'read_sounding_data', 552 )
1839
1840	!
1841	!-- Allocate temporary array for storing temperature data
1842	ALLOCATE( t_snd_tmp(1:nz_snd) )
1843	t_snd_tmp = 0.0_wp
1844
1845	!
1846	!-- Read temperature from file
1847	nc_stat = NF90_INQ_VARID( id, "ReferenceTemperature", id_var )
1848	nc_stat = NF90_GET_VAR( id, id_var, t_snd_tmp(:), start = (/1/), &
1849	count = (/nz_snd/) )
1850	CALL netcdf_handle_error_rad( 'read_sounding_data', 553 )
1851
1852	!
1853	!-- Calculate start of sounding data
1854	nz_snd_start = nz_snd + 1
1855	nz_snd_end = nz_snd + 1
1856
1857	!
1858	!-- Start filling vertical dimension at 10hPa above the model domain (hyp is
1859	!-- in Pa, hyp_snd in hPa).
1860	DO k = 1, nz_snd
1861	IF ( hyp_snd_tmp(k) < ( hyp(nzt+1) - 1000.0_wp) * 0.01_wp ) THEN
1862	nz_snd_start = k
1863	EXIT
1864	END IF
1865	END DO
1866
1867	IF ( nz_snd_start <= nz_snd ) THEN
1868	nz_snd_end = nz_snd
1869	END IF
1870
1871
1872	!
1873	!-- Calculate of total grid points for RRTMG calculations
1874	nzt_rad = nzt + nz_snd_end - nz_snd_start + 1
1875
1876	!
1877	!-- Save data above LES domain in hyp_snd, t_snd and q_snd
1878	!-- Note: q_snd_tmp is not calculated at the moment (dry residual atmosphere)
1879	ALLOCATE( hyp_snd(nzb+1:nzt_rad) )
1880	ALLOCATE( t_snd(nzb+1:nzt_rad) )
1881	ALLOCATE( q_snd(nzb+1:nzt_rad) )
1882	hyp_snd = 0.0_wp
1883	t_snd = 0.0_wp
1884	q_snd = 0.0_wp
1885
1886	hyp_snd(nzt+2:nzt_rad) = hyp_snd_tmp(nz_snd_start+1:nz_snd_end)
1887	t_snd(nzt+2:nzt_rad) = t_snd_tmp(nz_snd_start+1:nz_snd_end)
1888
1889	nc_stat = NF90_CLOSE( id )
1890
1891	!
1892	!-- Calculate pressure levels on zu and zw grid. Sounding data is added at
1893	!-- top of the LES domain. This routine does not consider horizontal or
1894	!-- vertical variability of pressure and temperature
1895	ALLOCATE ( rrtm_play(0:0,nzb+1:nzt_rad+1) )
1896	ALLOCATE ( rrtm_plev(0:0,nzb+1:nzt_rad+2) )
1897
1898	t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
1899	DO k = nzb+1, nzt+1
1900	rrtm_play(0,k) = hyp(k) * 0.01_wp
1901	rrtm_plev(0,k) = surface_pressure * ( (t_surface - g/cp * zw(k-1)) / &
1902	t_surface )**(1.0_wp/0.286_wp)
1903	ENDDO
1904
1905	DO k = nzt+2, nzt_rad
1906	rrtm_play(0,k) = hyp_snd(k)
1907	rrtm_plev(0,k) = 0.5_wp * ( rrtm_play(0,k) + rrtm_play(0,k-1) )
1908	ENDDO
1909	rrtm_plev(0,nzt_rad+1) = MAX( 0.5 * hyp_snd(nzt_rad), &
1910	1.5 * hyp_snd(nzt_rad) &
1911	- 0.5 * hyp_snd(nzt_rad-1) )
1912	rrtm_plev(0,nzt_rad+2) = MIN( 1.0E-4_wp, &
1913	0.25_wp * rrtm_plev(0,nzt_rad+1) )
1914
1915	rrtm_play(0,nzt_rad+1) = 0.5 * rrtm_plev(0,nzt_rad+1)
1916
1917	!
1918	!-- Calculate temperature/humidity levels at top of the LES domain.
1919	!-- Currently, the temperature is taken from sounding data (might lead to a
1920	!-- temperature jump at interface. To do: Humidity is currently not
1921	!-- calculated above the LES domain.
1922	ALLOCATE ( rrtm_tlay(0:0,nzb+1:nzt_rad+1) )
1923	ALLOCATE ( rrtm_tlev(0:0,nzb+1:nzt_rad+2) )
1924	ALLOCATE ( rrtm_h2ovmr(0:0,nzb+1:nzt_rad+1) )
1925
1926	DO k = nzt+8, nzt_rad
1927	rrtm_tlay(0,k) = t_snd(k)
1928	rrtm_h2ovmr(0,k) = q_snd(k)
1929	ENDDO
1930	rrtm_tlay(0,nzt_rad+1) = 2.0_wp * rrtm_tlay(0,nzt_rad) &
1931	- rrtm_tlay(0,nzt_rad-1)
1932	DO k = nzt+9, nzt_rad+1
1933	rrtm_tlev(0,k) = rrtm_tlay(0,k-1) + (rrtm_tlay(0,k) &
1934	- rrtm_tlay(0,k-1)) &
1935	/ ( rrtm_play(0,k) - rrtm_play(0,k-1) ) &
1936	* ( rrtm_plev(0,k) - rrtm_play(0,k-1) )
1937	ENDDO
1938	rrtm_h2ovmr(0,nzt_rad+1) = rrtm_h2ovmr(0,nzt_rad)
1939
1940	rrtm_tlev(0,nzt_rad+2) = 2.0_wp * rrtm_tlay(0,nzt_rad+1) &
1941	- rrtm_tlev(0,nzt_rad)
1942	!
1943	!-- Allocate remaining RRTMG arrays
1944	ALLOCATE ( rrtm_cicewp(0:0,nzb+1:nzt_rad+1) )
1945	ALLOCATE ( rrtm_cldfr(0:0,nzb+1:nzt_rad+1) )
1946	ALLOCATE ( rrtm_cliqwp(0:0,nzb+1:nzt_rad+1) )
1947	ALLOCATE ( rrtm_reice(0:0,nzb+1:nzt_rad+1) )
1948	ALLOCATE ( rrtm_reliq(0:0,nzb+1:nzt_rad+1) )
1949	ALLOCATE ( rrtm_lw_taucld(1:nbndlw+1,0:0,nzb+1:nzt_rad+1) )
1950	ALLOCATE ( rrtm_lw_tauaer(0:0,nzb+1:nzt_rad+1,1:nbndlw+1) )
1951	ALLOCATE ( rrtm_sw_taucld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
1952	ALLOCATE ( rrtm_sw_ssacld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
1953	ALLOCATE ( rrtm_sw_asmcld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
1954	ALLOCATE ( rrtm_sw_fsfcld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
1955	ALLOCATE ( rrtm_sw_tauaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) )
1956	ALLOCATE ( rrtm_sw_ssaaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) )
1957	ALLOCATE ( rrtm_sw_asmaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) )
1958	ALLOCATE ( rrtm_sw_ecaer(0:0,nzb+1:nzt_rad+1,1:naerec+1) )
1959
1960	!
1961	!-- The ice phase is currently not considered in PALM
1962	rrtm_cicewp = 0.0_wp
1963	rrtm_reice = 0.0_wp
1964
1965	!
1966	!-- Set other parameters (move to NAMELIST parameters in the future)
1967	rrtm_lw_tauaer = 0.0_wp
1968	rrtm_lw_taucld = 0.0_wp
1969	rrtm_sw_taucld = 0.0_wp
1970	rrtm_sw_ssacld = 0.0_wp
1971	rrtm_sw_asmcld = 0.0_wp
1972	rrtm_sw_fsfcld = 0.0_wp
1973	rrtm_sw_tauaer = 0.0_wp
1974	rrtm_sw_ssaaer = 0.0_wp
1975	rrtm_sw_asmaer = 0.0_wp
1976	rrtm_sw_ecaer = 0.0_wp
1977
1978
1979	ALLOCATE ( rrtm_swdflx(0:0,nzb:nzt_rad+1) )
1980	ALLOCATE ( rrtm_swuflx(0:0,nzb:nzt_rad+1) )
1981	ALLOCATE ( rrtm_swhr(0:0,nzb+1:nzt_rad+1) )
1982	ALLOCATE ( rrtm_swuflxc(0:0,nzb:nzt_rad+1) )
1983	ALLOCATE ( rrtm_swdflxc(0:0,nzb:nzt_rad+1) )
1984	ALLOCATE ( rrtm_swhrc(0:0,nzb+1:nzt_rad+1) )
1985
1986	rrtm_swdflx = 0.0_wp
1987	rrtm_swuflx = 0.0_wp
1988	rrtm_swhr = 0.0_wp
1989	rrtm_swuflxc = 0.0_wp
1990	rrtm_swdflxc = 0.0_wp
1991	rrtm_swhrc = 0.0_wp
1992
1993	ALLOCATE ( rrtm_lwdflx(0:0,nzb:nzt_rad+1) )
1994	ALLOCATE ( rrtm_lwuflx(0:0,nzb:nzt_rad+1) )
1995	ALLOCATE ( rrtm_lwhr(0:0,nzb+1:nzt_rad+1) )
1996	ALLOCATE ( rrtm_lwuflxc(0:0,nzb:nzt_rad+1) )
1997	ALLOCATE ( rrtm_lwdflxc(0:0,nzb:nzt_rad+1) )
1998	ALLOCATE ( rrtm_lwhrc(0:0,nzb+1:nzt_rad+1) )
1999
2000	rrtm_lwdflx = 0.0_wp
2001	rrtm_lwuflx = 0.0_wp
2002	rrtm_lwhr = 0.0_wp
2003	rrtm_lwuflxc = 0.0_wp
2004	rrtm_lwdflxc = 0.0_wp
2005	rrtm_lwhrc = 0.0_wp
2006
2007	ALLOCATE ( rrtm_lwuflx_dt(0:0,nzb:nzt_rad+1) )
2008	ALLOCATE ( rrtm_lwuflxc_dt(0:0,nzb:nzt_rad+1) )
2009
2010	rrtm_lwuflx_dt = 0.0_wp
2011	rrtm_lwuflxc_dt = 0.0_wp
2012
2013	END SUBROUTINE read_sounding_data
2014
2015
2016	!------------------------------------------------------------------------------!
2017	! Description:
2018	! ------------
2019	!> Read trace gas data from file
2020	!------------------------------------------------------------------------------!
2021	SUBROUTINE read_trace_gas_data
2022
2023	USE rrsw_ncpar
2024
2025	IMPLICIT NONE
2026
2027	INTEGER(iwp), PARAMETER :: num_trace_gases = 9 !< number of trace gases (absorbers)
2028
2029	CHARACTER(LEN=5), DIMENSION(num_trace_gases), PARAMETER :: & !< trace gas names
2030	trace_names = (/'O3 ', 'CO2 ', 'CH4 ', 'N2O ', 'O2 ', &
2031	'CFC11', 'CFC12', 'CFC22', 'CCL4 '/)
2032
2033	INTEGER(iwp) :: id, & !< NetCDF id
2034	k, & !< loop index
2035	m, & !< loop index
2036	n, & !< loop index
2037	nabs, & !< number of absorbers
2038	np, & !< number of pressure levels
2039	id_abs, & !< NetCDF id of the respective absorber
2040	id_dim, & !< NetCDF id of asborber's dimension
2041	id_var !< NetCDf id ot the absorber
2042
2043	REAL(wp) :: p_mls_l, p_mls_u, p_wgt_l, p_wgt_u, p_mls_m
2044
2045
2046	REAL(wp), DIMENSION(:), ALLOCATABLE :: p_mls, & !< pressure levels for the absorbers
2047	rrtm_play_tmp, & !< temporary array for pressure zu-levels
2048	rrtm_plev_tmp, & !< temporary array for pressure zw-levels
2049	trace_path_tmp !< temporary array for storing trace gas path data
2050
2051	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: trace_mls, & !< array for storing the absorber amounts
2052	trace_mls_path, & !< array for storing trace gas path data
2053	trace_mls_tmp !< temporary array for storing trace gas data
2054
2055
2056	!
2057	!-- In case of updates, deallocate arrays first (sufficient to check one
2058	!-- array as the others are automatically allocated)
2059	IF ( ALLOCATED ( rrtm_o3vmr ) ) THEN
2060	DEALLOCATE ( rrtm_o3vmr )
2061	DEALLOCATE ( rrtm_co2vmr )
2062	DEALLOCATE ( rrtm_ch4vmr )
2063	DEALLOCATE ( rrtm_n2ovmr )
2064	DEALLOCATE ( rrtm_o2vmr )
2065	DEALLOCATE ( rrtm_cfc11vmr )
2066	DEALLOCATE ( rrtm_cfc12vmr )
2067	DEALLOCATE ( rrtm_cfc22vmr )
2068	DEALLOCATE ( rrtm_ccl4vmr )
2069	ENDIF
2070
2071	!
2072	!-- Allocate trace gas profiles
2073	ALLOCATE ( rrtm_o3vmr(0:0,1:nzt_rad+1) )
2074	ALLOCATE ( rrtm_co2vmr(0:0,1:nzt_rad+1) )
2075	ALLOCATE ( rrtm_ch4vmr(0:0,1:nzt_rad+1) )
2076	ALLOCATE ( rrtm_n2ovmr(0:0,1:nzt_rad+1) )
2077	ALLOCATE ( rrtm_o2vmr(0:0,1:nzt_rad+1) )
2078	ALLOCATE ( rrtm_cfc11vmr(0:0,1:nzt_rad+1) )
2079	ALLOCATE ( rrtm_cfc12vmr(0:0,1:nzt_rad+1) )
2080	ALLOCATE ( rrtm_cfc22vmr(0:0,1:nzt_rad+1) )
2081	ALLOCATE ( rrtm_ccl4vmr(0:0,1:nzt_rad+1) )
2082
2083	!
2084	!-- Open file for reading
2085	nc_stat = NF90_OPEN( rrtm_input_file, NF90_NOWRITE, id )
2086	CALL netcdf_handle_error_rad( 'read_trace_gas_data', 549 )
2087	!
2088	!-- Inquire dimension ids and dimensions
2089	nc_stat = NF90_INQ_DIMID( id, "Pressure", id_dim )
2090	CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
2091	nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, len = np)
2092	CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
2093
2094	nc_stat = NF90_INQ_DIMID( id, "Absorber", id_dim )
2095	CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
2096	nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, len = nabs )
2097	CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
2098
2099
2100	!
2101	!-- Allocate pressure, and trace gas arrays
2102	ALLOCATE( p_mls(1:np) )
2103	ALLOCATE( trace_mls(1:num_trace_gases,1:np) )
2104	ALLOCATE( trace_mls_tmp(1:nabs,1:np) )
2105
2106
2107	nc_stat = NF90_INQ_VARID( id, "Pressure", id_var )
2108	CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
2109	nc_stat = NF90_GET_VAR( id, id_var, p_mls )
2110	CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
2111
2112	nc_stat = NF90_INQ_VARID( id, "AbsorberAmountMLS", id_var )
2113	CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
2114	nc_stat = NF90_GET_VAR( id, id_var, trace_mls_tmp )
2115	CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
2116
2117
2118	!
2119	!-- Write absorber amounts (mls) to trace_mls
2120	DO n = 1, num_trace_gases
2121	CALL getAbsorberIndex( TRIM( trace_names(n) ), id_abs )
2122
2123	trace_mls(n,1:np) = trace_mls_tmp(id_abs,1:np)
2124
2125	!
2126	!-- Replace missing values by zero
2127	WHERE ( trace_mls(n,:) > 2.0_wp )
2128	trace_mls(n,:) = 0.0_wp
2129	END WHERE
2130	END DO
2131
2132	DEALLOCATE ( trace_mls_tmp )
2133
2134	nc_stat = NF90_CLOSE( id )
2135	CALL netcdf_handle_error_rad( 'read_trace_gas_data', 551 )
2136
2137	!
2138	!-- Add extra pressure level for calculations of the trace gas paths
2139	ALLOCATE ( rrtm_play_tmp(1:nzt_rad+1) )
2140	ALLOCATE ( rrtm_plev_tmp(1:nzt_rad+2) )
2141
2142	rrtm_play_tmp(1:nzt_rad) = rrtm_play(0,1:nzt_rad)
2143	rrtm_plev_tmp(1:nzt_rad+1) = rrtm_plev(0,1:nzt_rad+1)
2144	rrtm_play_tmp(nzt_rad+1) = rrtm_plev(0,nzt_rad+1) * 0.5_wp
2145	rrtm_plev_tmp(nzt_rad+2) = MIN( 1.0E-4_wp, 0.25_wp &
2146	* rrtm_plev(0,nzt_rad+1) )
2147
2148	!
2149	!-- Calculate trace gas path (zero at surface) with interpolation to the
2150	!-- sounding levels
2151	ALLOCATE ( trace_mls_path(1:nzt_rad+2,1:num_trace_gases) )
2152
2153	trace_mls_path(nzb+1,:) = 0.0_wp
2154
2155	DO k = nzb+2, nzt_rad+2
2156	DO m = 1, num_trace_gases
2157	trace_mls_path(k,m) = trace_mls_path(k-1,m)
2158
2159	!
2160	!-- When the pressure level is higher than the trace gas pressure
2161	!-- level, assume that
2162	IF ( rrtm_plev_tmp(k-1) > p_mls(1) ) THEN
2163
2164	trace_mls_path(k,m) = trace_mls_path(k,m) + trace_mls(m,1) &
2165	* ( rrtm_plev_tmp(k-1) &
2166	- MAX( p_mls(1), rrtm_plev_tmp(k) ) &
2167	) / g
2168	ENDIF
2169
2170	!
2171	!-- Integrate for each sounding level from the contributing p_mls
2172	!-- levels
2173	DO n = 2, np
2174	!
2175	!-- Limit p_mls so that it is within the model level
2176	p_mls_u = MIN( rrtm_plev_tmp(k-1), &
2177	MAX( rrtm_plev_tmp(k), p_mls(n) ) )
2178	p_mls_l = MIN( rrtm_plev_tmp(k-1), &
2179	MAX( rrtm_plev_tmp(k), p_mls(n-1) ) )
2180
2181	IF ( p_mls_l > p_mls_u ) THEN
2182
2183	!
2184	!-- Calculate weights for interpolation
2185	p_mls_m = 0.5_wp * (p_mls_l + p_mls_u)
2186	p_wgt_u = (p_mls(n-1) - p_mls_m) / (p_mls(n-1) - p_mls(n))
2187	p_wgt_l = (p_mls_m - p_mls(n)) / (p_mls(n-1) - p_mls(n))
2188
2189	!
2190	!-- Add level to trace gas path
2191	trace_mls_path(k,m) = trace_mls_path(k,m) &
2192	+ ( p_wgt_u * trace_mls(m,n) &
2193	+ p_wgt_l * trace_mls(m,n-1) ) &
2194	* (p_mls_l - p_mls_u) / g
2195	ENDIF
2196	ENDDO
2197
2198	IF ( rrtm_plev_tmp(k) < p_mls(np) ) THEN
2199	trace_mls_path(k,m) = trace_mls_path(k,m) + trace_mls(m,np) &
2200	* ( MIN( rrtm_plev_tmp(k-1), p_mls(np) ) &
2201	- rrtm_plev_tmp(k) &
2202	) / g
2203	ENDIF
2204	ENDDO
2205	ENDDO
2206
2207
2208	!
2209	!-- Prepare trace gas path profiles
2210	ALLOCATE ( trace_path_tmp(1:nzt_rad+1) )
2211
2212	DO m = 1, num_trace_gases
2213
2214	trace_path_tmp(1:nzt_rad+1) = ( trace_mls_path(2:nzt_rad+2,m) &
2215	- trace_mls_path(1:nzt_rad+1,m) ) * g &
2216	/ ( rrtm_plev_tmp(1:nzt_rad+1) &
2217	- rrtm_plev_tmp(2:nzt_rad+2) )
2218
2219	!
2220	!-- Save trace gas paths to the respective arrays
2221	SELECT CASE ( TRIM( trace_names(m) ) )
2222
2223	CASE ( 'O3' )
2224
2225	rrtm_o3vmr(0,:) = trace_path_tmp(:)
2226
2227	CASE ( 'CO2' )
2228
2229	rrtm_co2vmr(0,:) = trace_path_tmp(:)
2230
2231	CASE ( 'CH4' )
2232
2233	rrtm_ch4vmr(0,:) = trace_path_tmp(:)
2234
2235	CASE ( 'N2O' )
2236
2237	rrtm_n2ovmr(0,:) = trace_path_tmp(:)
2238
2239	CASE ( 'O2' )
2240
2241	rrtm_o2vmr(0,:) = trace_path_tmp(:)
2242
2243	CASE ( 'CFC11' )
2244
2245	rrtm_cfc11vmr(0,:) = trace_path_tmp(:)
2246
2247	CASE ( 'CFC12' )
2248
2249	rrtm_cfc12vmr(0,:) = trace_path_tmp(:)
2250
2251	CASE ( 'CFC22' )
2252
2253	rrtm_cfc22vmr(0,:) = trace_path_tmp(:)
2254
2255	CASE ( 'CCL4' )
2256
2257	rrtm_ccl4vmr(0,:) = trace_path_tmp(:)
2258
2259	CASE DEFAULT
2260
2261	END SELECT
2262
2263	ENDDO
2264
2265	DEALLOCATE ( trace_path_tmp )
2266	DEALLOCATE ( trace_mls_path )
2267	DEALLOCATE ( rrtm_play_tmp )
2268	DEALLOCATE ( rrtm_plev_tmp )
2269	DEALLOCATE ( trace_mls )
2270	DEALLOCATE ( p_mls )
2271
2272	END SUBROUTINE read_trace_gas_data
2273
2274
2275	SUBROUTINE netcdf_handle_error_rad( routine_name, errno )
2276
2277	USE control_parameters, &
2278	ONLY: message_string
2279
2280	USE NETCDF
2281
2282	USE pegrid
2283
2284	IMPLICIT NONE
2285
2286	CHARACTER(LEN=6) :: message_identifier
2287	CHARACTER(LEN=*) :: routine_name
2288
2289	INTEGER(iwp) :: errno
2290
2291	IF ( nc_stat /= NF90_NOERR ) THEN
2292
2293	WRITE( message_identifier, '(''NC'',I4.4)' ) errno
2294	message_string = TRIM( NF90_STRERROR( nc_stat ) )
2295
2296	CALL message( routine_name, message_identifier, 2, 2, 0, 6, 1 )
2297
2298	ENDIF
2299
2300	END SUBROUTINE netcdf_handle_error_rad
2301	#endif
2302
2303
2304	!------------------------------------------------------------------------------!
2305	! Description:
2306	! ------------
2307	!> Calculate temperature tendency due to radiative cooling/heating.
2308	!> Cache-optimized version.
2309	!------------------------------------------------------------------------------!
2310	SUBROUTINE radiation_tendency_ij ( i, j, tend )
2311
2312	USE cloud_parameters, &
2313	ONLY: pt_d_t
2314
2315	IMPLICIT NONE
2316
2317	INTEGER(iwp) :: i, j, k !< loop indices
2318
2319	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend !< pt tendency term
2320
2321	IF ( radiation_scheme == 'rrtmg' ) THEN
2322	#if defined ( __rrtmg )
2323	!
2324	!-- Calculate tendency based on heating rate
2325	DO k = nzb+1, nzt+1
2326	tend(k,j,i) = tend(k,j,i) + (rad_lw_hr(k,j,i) + rad_sw_hr(k,j,i)) &
2327	* pt_d_t(k) * d_seconds_hour
2328	ENDDO
2329	#endif
2330	ENDIF
2331
2332	END SUBROUTINE radiation_tendency_ij
2333
2334
2335	!------------------------------------------------------------------------------!
2336	! Description:
2337	! ------------
2338	!> Calculate temperature tendency due to radiative cooling/heating.
2339	!> Vector-optimized version
2340	!------------------------------------------------------------------------------!
2341	SUBROUTINE radiation_tendency ( tend )
2342
2343	USE cloud_parameters, &
2344	ONLY: pt_d_t
2345
2346	USE indices, &
2347	ONLY: nxl, nxr, nyn, nys
2348
2349	IMPLICIT NONE
2350
2351	INTEGER(iwp) :: i, j, k !< loop indices
2352
2353	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend !< pt tendency term
2354
2355	IF ( radiation_scheme == 'rrtmg' ) THEN
2356	#if defined ( __rrtmg )
2357	!
2358	!-- Calculate tendency based on heating rate
2359	DO i = nxl, nxr
2360	DO j = nys, nyn
2361	DO k = nzb+1, nzt+1
2362	tend(k,j,i) = tend(k,j,i) + ( rad_lw_hr(k,j,i) &
2363	+ rad_sw_hr(k,j,i) ) * pt_d_t(k) &
2364	* d_seconds_hour
2365	ENDDO
2366	ENDDO
2367	ENDDO
2368	#endif
2369	ENDIF
2370
2371
2372	END SUBROUTINE radiation_tendency
2373
2374	!------------------------------------------------------------------------------!
2375	!
2376	! Description:
2377	! ------------
2378	!> Subroutine for averaging 3D data
2379	!------------------------------------------------------------------------------!
2380	SUBROUTINE radiation_3d_data_averaging( mode, variable )
2381
2382
2383	USE control_parameters
2384
2385	USE indices
2386
2387	USE kinds
2388
2389	IMPLICIT NONE
2390
2391	CHARACTER (LEN=*) :: mode !<
2392	CHARACTER (LEN=*) :: variable !<
2393
2394	INTEGER(iwp) :: i !<
2395	INTEGER(iwp) :: j !<
2396	INTEGER(iwp) :: k !<
2397
2398	IF ( mode == 'allocate' ) THEN
2399
2400	SELECT CASE ( TRIM( variable ) )
2401
2402	CASE ( 'rad_net*' )
2403	IF ( .NOT. ALLOCATED( rad_net_av ) ) THEN
2404	ALLOCATE( rad_net_av(nysg:nyng,nxlg:nxrg) )
2405	ENDIF
2406	rad_net_av = 0.0_wp
2407
2408	CASE ( 'rad_lw_in' )
2409	IF ( .NOT. ALLOCATED( rad_lw_in_av ) ) THEN
2410	ALLOCATE( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2411	ENDIF
2412	rad_lw_in_av = 0.0_wp
2413
2414	CASE ( 'rad_lw_out' )
2415	IF ( .NOT. ALLOCATED( rad_lw_out_av ) ) THEN
2416	ALLOCATE( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2417	ENDIF
2418	rad_lw_out_av = 0.0_wp
2419
2420	CASE ( 'rad_lw_cs_hr' )
2421	IF ( .NOT. ALLOCATED( rad_lw_cs_hr_av ) ) THEN
2422	ALLOCATE( rad_lw_cs_hr_av(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
2423	ENDIF
2424	rad_lw_cs_hr_av = 0.0_wp
2425
2426	CASE ( 'rad_lw_hr' )
2427	IF ( .NOT. ALLOCATED( rad_lw_hr_av ) ) THEN
2428	ALLOCATE( rad_lw_hr_av(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
2429	ENDIF
2430	rad_lw_hr_av = 0.0_wp
2431
2432	CASE ( 'rad_sw_in' )
2433	IF ( .NOT. ALLOCATED( rad_sw_in_av ) ) THEN
2434	ALLOCATE( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2435	ENDIF
2436	rad_sw_in_av = 0.0_wp
2437
2438	CASE ( 'rad_sw_out' )
2439	IF ( .NOT. ALLOCATED( rad_sw_out_av ) ) THEN
2440	ALLOCATE( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2441	ENDIF
2442	rad_sw_out_av = 0.0_wp
2443
2444	CASE ( 'rad_sw_cs_hr' )
2445	IF ( .NOT. ALLOCATED( rad_sw_cs_hr_av ) ) THEN
2446	ALLOCATE( rad_sw_cs_hr_av(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
2447	ENDIF
2448	rad_sw_cs_hr_av = 0.0_wp
2449
2450	CASE ( 'rad_sw_hr' )
2451	IF ( .NOT. ALLOCATED( rad_sw_hr_av ) ) THEN
2452	ALLOCATE( rad_sw_hr_av(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
2453	ENDIF
2454	rad_sw_hr_av = 0.0_wp
2455
2456	CASE DEFAULT
2457	CONTINUE
2458
2459	END SELECT
2460
2461	ELSEIF ( mode == 'sum' ) THEN
2462
2463	SELECT CASE ( TRIM( variable ) )
2464
2465	CASE ( 'rad_net*' )
2466	DO i = nxlg, nxrg
2467	DO j = nysg, nyng
2468	rad_net_av(j,i) = rad_net_av(j,i) + rad_net(j,i)
2469	ENDDO
2470	ENDDO
2471
2472	CASE ( 'rad_lw_in' )
2473	DO i = nxlg, nxrg
2474	DO j = nysg, nyng
2475	DO k = nzb, nzt+1
2476	rad_lw_in_av(k,j,i) = rad_lw_in_av(k,j,i) + rad_lw_in(k,j,i)
2477	ENDDO
2478	ENDDO
2479	ENDDO
2480
2481	CASE ( 'rad_lw_out' )
2482	DO i = nxlg, nxrg
2483	DO j = nysg, nyng
2484	DO k = nzb, nzt+1
2485	rad_lw_out_av(k,j,i) = rad_lw_out_av(k,j,i) + rad_lw_out(k,j,i)
2486	ENDDO
2487	ENDDO
2488	ENDDO
2489
2490	CASE ( 'rad_lw_cs_hr' )
2491	DO i = nxlg, nxrg
2492	DO j = nysg, nyng
2493	DO k = nzb, nzt+1
2494	rad_lw_cs_hr_av(k,j,i) = rad_lw_cs_hr_av(k,j,i) + rad_lw_cs_hr(k,j,i)
2495	ENDDO
2496	ENDDO
2497	ENDDO
2498
2499	CASE ( 'rad_lw_hr' )
2500	DO i = nxlg, nxrg
2501	DO j = nysg, nyng
2502	DO k = nzb, nzt+1
2503	rad_lw_hr_av(k,j,i) = rad_lw_hr_av(k,j,i) + rad_lw_hr(k,j,i)
2504	ENDDO
2505	ENDDO
2506	ENDDO
2507
2508	CASE ( 'rad_sw_in' )
2509	DO i = nxlg, nxrg
2510	DO j = nysg, nyng
2511	DO k = nzb, nzt+1
2512	rad_sw_in_av(k,j,i) = rad_sw_in_av(k,j,i) + rad_sw_in(k,j,i)
2513	ENDDO
2514	ENDDO
2515	ENDDO
2516
2517	CASE ( 'rad_sw_out' )
2518	DO i = nxlg, nxrg
2519	DO j = nysg, nyng
2520	DO k = nzb, nzt+1
2521	rad_sw_out_av(k,j,i) = rad_sw_out_av(k,j,i) + rad_sw_out(k,j,i)
2522	ENDDO
2523	ENDDO
2524	ENDDO
2525
2526	CASE ( 'rad_sw_cs_hr' )
2527	DO i = nxlg, nxrg
2528	DO j = nysg, nyng
2529	DO k = nzb, nzt+1
2530	rad_sw_cs_hr_av(k,j,i) = rad_sw_cs_hr_av(k,j,i) + rad_sw_cs_hr(k,j,i)
2531	ENDDO
2532	ENDDO
2533	ENDDO
2534
2535	CASE ( 'rad_sw_hr' )
2536	DO i = nxlg, nxrg
2537	DO j = nysg, nyng
2538	DO k = nzb, nzt+1
2539	rad_sw_hr_av(k,j,i) = rad_sw_hr_av(k,j,i) + rad_sw_hr(k,j,i)
2540	ENDDO
2541	ENDDO
2542	ENDDO
2543
2544	CASE DEFAULT
2545	CONTINUE
2546
2547	END SELECT
2548
2549	ELSEIF ( mode == 'average' ) THEN
2550
2551	SELECT CASE ( TRIM( variable ) )
2552
2553	CASE ( 'rad_net*' )
2554	DO i = nxlg, nxrg
2555	DO j = nysg, nyng
2556	rad_net_av(j,i) = rad_net_av(j,i) / REAL( average_count_3d, KIND=wp )
2557	ENDDO
2558	ENDDO
2559
2560	CASE ( 'rad_lw_in' )
2561	DO i = nxlg, nxrg
2562	DO j = nysg, nyng
2563	DO k = nzb, nzt+1
2564	rad_lw_in_av(k,j,i) = rad_lw_in_av(k,j,i) / REAL( average_count_3d, KIND=wp )
2565	ENDDO
2566	ENDDO
2567	ENDDO
2568
2569	CASE ( 'rad_lw_out' )
2570	DO i = nxlg, nxrg
2571	DO j = nysg, nyng
2572	DO k = nzb, nzt+1
2573	rad_lw_out_av(k,j,i) = rad_lw_out_av(k,j,i) / REAL( average_count_3d, KIND=wp )
2574	ENDDO
2575	ENDDO
2576	ENDDO
2577
2578	CASE ( 'rad_lw_cs_hr' )
2579	DO i = nxlg, nxrg
2580	DO j = nysg, nyng
2581	DO k = nzb, nzt+1
2582	rad_lw_cs_hr_av(k,j,i) = rad_lw_cs_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
2583	ENDDO
2584	ENDDO
2585	ENDDO
2586
2587	CASE ( 'rad_lw_hr' )
2588	DO i = nxlg, nxrg
2589	DO j = nysg, nyng
2590	DO k = nzb, nzt+1
2591	rad_lw_hr_av(k,j,i) = rad_lw_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
2592	ENDDO
2593	ENDDO
2594	ENDDO
2595
2596	CASE ( 'rad_sw_in' )
2597	DO i = nxlg, nxrg
2598	DO j = nysg, nyng
2599	DO k = nzb, nzt+1
2600	rad_sw_in_av(k,j,i) = rad_sw_in_av(k,j,i) / REAL( average_count_3d, KIND=wp )
2601	ENDDO
2602	ENDDO
2603	ENDDO
2604
2605	CASE ( 'rad_sw_out' )
2606	DO i = nxlg, nxrg
2607	DO j = nysg, nyng
2608	DO k = nzb, nzt+1
2609	rad_sw_out_av(k,j,i) = rad_sw_out_av(k,j,i) / REAL( average_count_3d, KIND=wp )
2610	ENDDO
2611	ENDDO
2612	ENDDO
2613
2614	CASE ( 'rad_sw_cs_hr' )
2615	DO i = nxlg, nxrg
2616	DO j = nysg, nyng
2617	DO k = nzb, nzt+1
2618	rad_sw_cs_hr_av(k,j,i) = rad_sw_cs_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
2619	ENDDO
2620	ENDDO
2621	ENDDO
2622
2623	CASE ( 'rad_sw_hr' )
2624	DO i = nxlg, nxrg
2625	DO j = nysg, nyng
2626	DO k = nzb, nzt+1
2627	rad_sw_hr_av(k,j,i) = rad_sw_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
2628	ENDDO
2629	ENDDO
2630	ENDDO
2631
2632	END SELECT
2633
2634	ENDIF
2635
2636	END SUBROUTINE radiation_3d_data_averaging
2637
2638
2639	!------------------------------------------------------------------------------!
2640	!
2641	! Description:
2642	! ------------
2643	!> Subroutine defining appropriate grid for netcdf variables.
2644	!> It is called out from subroutine netcdf.
2645	!------------------------------------------------------------------------------!
2646	SUBROUTINE radiation_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
2647
2648	IMPLICIT NONE
2649
2650	CHARACTER (LEN=*), INTENT(IN) :: var !<
2651	LOGICAL, INTENT(OUT) :: found !<
2652	CHARACTER (LEN=*), INTENT(OUT) :: grid_x !<
2653	CHARACTER (LEN=*), INTENT(OUT) :: grid_y !<
2654	CHARACTER (LEN=*), INTENT(OUT) :: grid_z !<
2655
2656	found = .TRUE.
2657
2658
2659	!
2660	!-- Check for the grid
2661	SELECT CASE ( TRIM( var ) )
2662
2663	CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_sw_cs_hr', 'rad_sw_hr', &
2664	'rad_lw_cs_hr_xy', 'rad_lw_hr_xy', 'rad_sw_cs_hr_xy', &
2665	'rad_sw_hr_xy', 'rad_lw_cs_hr_xz', 'rad_lw_hr_xz', &
2666	'rad_sw_cs_hr_xz', 'rad_sw_hr_xz', 'rad_lw_cs_hr_yz', &
2667	'rad_lw_hr_yz', 'rad_sw_cs_hr_yz', 'rad_sw_hr_yz' )
2668	grid_x = 'x'
2669	grid_y = 'y'
2670	grid_z = 'zu'
2671
2672	CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_sw_in', 'rad_sw_out', &
2673	'rad_lw_in_xy', 'rad_lw_out_xy', 'rad_sw_in_xy','rad_sw_out_xy', &
2674	'rad_lw_in_xz', 'rad_lw_out_xz', 'rad_sw_in_xz','rad_sw_out_xz', &
2675	'rad_lw_in_yz', 'rad_lw_out_yz', 'rad_sw_in_yz','rad_sw_out_yz' )
2676	grid_x = 'x'
2677	grid_y = 'y'
2678	grid_z = 'zw'
2679
2680
2681	CASE DEFAULT
2682	found = .FALSE.
2683	grid_x = 'none'
2684	grid_y = 'none'
2685	grid_z = 'none'
2686
2687	END SELECT
2688
2689	END SUBROUTINE radiation_define_netcdf_grid
2690
2691	!------------------------------------------------------------------------------!
2692	!
2693	! Description:
2694	! ------------
2695	!> Subroutine defining 3D output variables
2696	!------------------------------------------------------------------------------!
2697	SUBROUTINE radiation_data_output_2d( av, variable, found, grid, mode, &
2698	local_pf, two_d )
2699
2700	USE indices
2701
2702	USE kinds
2703
2704
2705	IMPLICIT NONE
2706
2707	CHARACTER (LEN=*) :: grid !<
2708	CHARACTER (LEN=*) :: mode !<
2709	CHARACTER (LEN=*) :: variable !<
2710
2711	INTEGER(iwp) :: av !<
2712	INTEGER(iwp) :: i !<
2713	INTEGER(iwp) :: j !<
2714	INTEGER(iwp) :: k !<
2715
2716	LOGICAL :: found !<
2717	LOGICAL :: two_d !< flag parameter that indicates 2D variables (horizontal cross sections)
2718
2719	REAL(wp), DIMENSION(nxlg:nxrg,nysg:nyng,nzb:nzt+1) :: local_pf !<
2720
2721	found = .TRUE.
2722
2723	SELECT CASE ( TRIM( variable ) )
2724
2725	CASE ( 'rad_net*_xy' ) ! 2d-array
2726	IF ( av == 0 ) THEN
2727	DO i = nxlg, nxrg
2728	DO j = nysg, nyng
2729	local_pf(i,j,nzb+1) = rad_net(j,i)
2730	ENDDO
2731	ENDDO
2732	ELSE
2733	DO i = nxlg, nxrg
2734	DO j = nysg, nyng
2735	local_pf(i,j,nzb+1) = rad_net_av(j,i)
2736	ENDDO
2737	ENDDO
2738	ENDIF
2739	two_d = .TRUE.
2740	grid = 'zu1'
2741
2742
2743	CASE ( 'rad_lw_in_xy', 'rad_lw_in_xz', 'rad_lw_in_yz' )
2744	IF ( av == 0 ) THEN
2745	DO i = nxlg, nxrg
2746	DO j = nysg, nyng
2747	DO k = nzb, nzt+1
2748	local_pf(i,j,k) = rad_lw_in(k,j,i)
2749	ENDDO
2750	ENDDO
2751	ENDDO
2752	ELSE
2753	DO i = nxlg, nxrg
2754	DO j = nysg, nyng
2755	DO k = nzb, nzt+1
2756	local_pf(i,j,k) = rad_lw_in_av(k,j,i)
2757	ENDDO
2758	ENDDO
2759	ENDDO
2760	ENDIF
2761	IF ( mode == 'xy' ) grid = 'zu'
2762
2763	CASE ( 'rad_lw_out_xy', 'rad_lw_out_xz', 'rad_lw_out_yz' )
2764	IF ( av == 0 ) THEN
2765	DO i = nxlg, nxrg
2766	DO j = nysg, nyng
2767	DO k = nzb, nzt+1
2768	local_pf(i,j,k) = rad_lw_out(k,j,i)
2769	ENDDO
2770	ENDDO
2771	ENDDO
2772	ELSE
2773	DO i = nxlg, nxrg
2774	DO j = nysg, nyng
2775	DO k = nzb, nzt+1
2776	local_pf(i,j,k) = rad_lw_out_av(k,j,i)
2777	ENDDO
2778	ENDDO
2779	ENDDO
2780	ENDIF
2781	IF ( mode == 'xy' ) grid = 'zu'
2782
2783	CASE ( 'rad_lw_cs_hr_xy', 'rad_lw_cs_hr_xz', 'rad_lw_cs_hr_yz' )
2784	IF ( av == 0 ) THEN
2785	DO i = nxlg, nxrg
2786	DO j = nysg, nyng
2787	DO k = nzb, nzt+1
2788	local_pf(i,j,k) = rad_lw_cs_hr(k,j,i)
2789	ENDDO
2790	ENDDO
2791	ENDDO
2792	ELSE
2793	DO i = nxlg, nxrg
2794	DO j = nysg, nyng
2795	DO k = nzb, nzt+1
2796	local_pf(i,j,k) = rad_lw_cs_hr_av(k,j,i)
2797	ENDDO
2798	ENDDO
2799	ENDDO
2800	ENDIF
2801	IF ( mode == 'xy' ) grid = 'zw'
2802
2803	CASE ( 'rad_lw_hr_xy', 'rad_lw_hr_xz', 'rad_lw_hr_yz' )
2804	IF ( av == 0 ) THEN
2805	DO i = nxlg, nxrg
2806	DO j = nysg, nyng
2807	DO k = nzb, nzt+1
2808	local_pf(i,j,k) = rad_lw_hr(k,j,i)
2809	ENDDO
2810	ENDDO
2811	ENDDO
2812	ELSE
2813	DO i = nxlg, nxrg
2814	DO j = nysg, nyng
2815	DO k = nzb, nzt+1
2816	local_pf(i,j,k) = rad_lw_hr_av(k,j,i)
2817	ENDDO
2818	ENDDO
2819	ENDDO
2820	ENDIF
2821	IF ( mode == 'xy' ) grid = 'zw'
2822
2823	CASE ( 'rad_sw_in_xy', 'rad_sw_in_xz', 'rad_sw_in_yz' )
2824	IF ( av == 0 ) THEN
2825	DO i = nxlg, nxrg
2826	DO j = nysg, nyng
2827	DO k = nzb, nzt+1
2828	local_pf(i,j,k) = rad_sw_in(k,j,i)
2829	ENDDO
2830	ENDDO
2831	ENDDO
2832	ELSE
2833	DO i = nxlg, nxrg
2834	DO j = nysg, nyng
2835	DO k = nzb, nzt+1
2836	local_pf(i,j,k) = rad_sw_in_av(k,j,i)
2837	ENDDO
2838	ENDDO
2839	ENDDO
2840	ENDIF
2841	IF ( mode == 'xy' ) grid = 'zu'
2842
2843	CASE ( 'rad_sw_out_xy', 'rad_sw_out_xz', 'rad_sw_out_yz' )
2844	IF ( av == 0 ) THEN
2845	DO i = nxlg, nxrg
2846	DO j = nysg, nyng
2847	DO k = nzb, nzt+1
2848	local_pf(i,j,k) = rad_sw_out(k,j,i)
2849	ENDDO
2850	ENDDO
2851	ENDDO
2852	ELSE
2853	DO i = nxlg, nxrg
2854	DO j = nysg, nyng
2855	DO k = nzb, nzt+1
2856	local_pf(i,j,k) = rad_sw_out_av(k,j,i)
2857	ENDDO
2858	ENDDO
2859	ENDDO
2860	ENDIF
2861	IF ( mode == 'xy' ) grid = 'zu'
2862
2863	CASE ( 'rad_sw_cs_hr_xy', 'rad_sw_cs_hr_xz', 'rad_sw_cs_hr_yz' )
2864	IF ( av == 0 ) THEN
2865	DO i = nxlg, nxrg
2866	DO j = nysg, nyng
2867	DO k = nzb, nzt+1
2868	local_pf(i,j,k) = rad_sw_cs_hr(k,j,i)
2869	ENDDO
2870	ENDDO
2871	ENDDO
2872	ELSE
2873	DO i = nxlg, nxrg
2874	DO j = nysg, nyng
2875	DO k = nzb, nzt+1
2876	local_pf(i,j,k) = rad_sw_cs_hr_av(k,j,i)
2877	ENDDO
2878	ENDDO
2879	ENDDO
2880	ENDIF
2881	IF ( mode == 'xy' ) grid = 'zw'
2882
2883	CASE ( 'rad_sw_hr_xy', 'rad_sw_hr_xz', 'rad_sw_hr_yz' )
2884	IF ( av == 0 ) THEN
2885	DO i = nxlg, nxrg
2886	DO j = nysg, nyng
2887	DO k = nzb, nzt+1
2888	local_pf(i,j,k) = rad_sw_hr(k,j,i)
2889	ENDDO
2890	ENDDO
2891	ENDDO
2892	ELSE
2893	DO i = nxlg, nxrg
2894	DO j = nysg, nyng
2895	DO k = nzb, nzt+1
2896	local_pf(i,j,k) = rad_sw_hr_av(k,j,i)
2897	ENDDO
2898	ENDDO
2899	ENDDO
2900	ENDIF
2901	IF ( mode == 'xy' ) grid = 'zw'
2902
2903	CASE DEFAULT
2904	found = .FALSE.
2905	grid = 'none'
2906
2907	END SELECT
2908
2909	END SUBROUTINE radiation_data_output_2d
2910
2911
2912	!------------------------------------------------------------------------------!
2913	!
2914	! Description:
2915	! ------------
2916	!> Subroutine defining 3D output variables
2917	!------------------------------------------------------------------------------!
2918	SUBROUTINE radiation_data_output_3d( av, variable, found, local_pf )
2919
2920
2921	USE indices
2922
2923	USE kinds
2924
2925
2926	IMPLICIT NONE
2927
2928	CHARACTER (LEN=*) :: variable !<
2929
2930	INTEGER(iwp) :: av !<
2931	INTEGER(iwp) :: i !<
2932	INTEGER(iwp) :: j !<
2933	INTEGER(iwp) :: k !<
2934
2935	LOGICAL :: found !<
2936
2937	REAL(sp), DIMENSION(nxlg:nxrg,nysg:nyng,nzb:nzt+1) :: local_pf !<
2938
2939
2940	found = .TRUE.
2941
2942
2943	SELECT CASE ( TRIM( variable ) )
2944
2945	CASE ( 'rad_sw_in' )
2946	IF ( av == 0 ) THEN
2947	DO i = nxlg, nxrg
2948	DO j = nysg, nyng
2949	DO k = nzb, nzt+1
2950	local_pf(i,j,k) = rad_sw_in(k,j,i)
2951	ENDDO
2952	ENDDO
2953	ENDDO
2954	ELSE
2955	DO i = nxlg, nxrg
2956	DO j = nysg, nyng
2957	DO k = nzb, nzt+1
2958	local_pf(i,j,k) = rad_sw_in_av(k,j,i)
2959	ENDDO
2960	ENDDO
2961	ENDDO
2962	ENDIF
2963
2964	CASE ( 'rad_sw_out' )
2965	IF ( av == 0 ) THEN
2966	DO i = nxlg, nxrg
2967	DO j = nysg, nyng
2968	DO k = nzb, nzt+1
2969	local_pf(i,j,k) = rad_sw_out(k,j,i)
2970	ENDDO
2971	ENDDO
2972	ENDDO
2973	ELSE
2974	DO i = nxlg, nxrg
2975	DO j = nysg, nyng
2976	DO k = nzb, nzt+1
2977	local_pf(i,j,k) = rad_sw_out_av(k,j,i)
2978	ENDDO
2979	ENDDO
2980	ENDDO
2981	ENDIF
2982
2983	CASE ( 'rad_sw_cs_hr' )
2984	IF ( av == 0 ) THEN
2985	DO i = nxlg, nxrg
2986	DO j = nysg, nyng
2987	DO k = nzb, nzt+1
2988	local_pf(i,j,k) = rad_sw_cs_hr(k,j,i)
2989	ENDDO
2990	ENDDO
2991	ENDDO
2992	ELSE
2993	DO i = nxlg, nxrg
2994	DO j = nysg, nyng
2995	DO k = nzb, nzt+1
2996	local_pf(i,j,k) = rad_sw_cs_hr_av(k,j,i)
2997	ENDDO
2998	ENDDO
2999	ENDDO
3000	ENDIF
3001
3002	CASE ( 'rad_sw_hr' )
3003	IF ( av == 0 ) THEN
3004	DO i = nxlg, nxrg
3005	DO j = nysg, nyng
3006	DO k = nzb, nzt+1
3007	local_pf(i,j,k) = rad_sw_hr(k,j,i)
3008	ENDDO
3009	ENDDO
3010	ENDDO
3011	ELSE
3012	DO i = nxlg, nxrg
3013	DO j = nysg, nyng
3014	DO k = nzb, nzt+1
3015	local_pf(i,j,k) = rad_sw_hr_av(k,j,i)
3016	ENDDO
3017	ENDDO
3018	ENDDO
3019	ENDIF
3020
3021	CASE ( 'rad_lw_in' )
3022	IF ( av == 0 ) THEN
3023	DO i = nxlg, nxrg
3024	DO j = nysg, nyng
3025	DO k = nzb, nzt+1
3026	local_pf(i,j,k) = rad_lw_in(k,j,i)
3027	ENDDO
3028	ENDDO
3029	ENDDO
3030	ELSE
3031	DO i = nxlg, nxrg
3032	DO j = nysg, nyng
3033	DO k = nzb, nzt+1
3034	local_pf(i,j,k) = rad_lw_in_av(k,j,i)
3035	ENDDO
3036	ENDDO
3037	ENDDO
3038	ENDIF
3039
3040	CASE ( 'rad_lw_out' )
3041	IF ( av == 0 ) THEN
3042	DO i = nxlg, nxrg
3043	DO j = nysg, nyng
3044	DO k = nzb, nzt+1
3045	local_pf(i,j,k) = rad_lw_out(k,j,i)
3046	ENDDO
3047	ENDDO
3048	ENDDO
3049	ELSE
3050	DO i = nxlg, nxrg
3051	DO j = nysg, nyng
3052	DO k = nzb, nzt+1
3053	local_pf(i,j,k) = rad_lw_out_av(k,j,i)
3054	ENDDO
3055	ENDDO
3056	ENDDO
3057	ENDIF
3058
3059	CASE ( 'rad_lw_cs_hr' )
3060	IF ( av == 0 ) THEN
3061	DO i = nxlg, nxrg
3062	DO j = nysg, nyng
3063	DO k = nzb, nzt+1
3064	local_pf(i,j,k) = rad_lw_cs_hr(k,j,i)
3065	ENDDO
3066	ENDDO
3067	ENDDO
3068	ELSE
3069	DO i = nxlg, nxrg
3070	DO j = nysg, nyng
3071	DO k = nzb, nzt+1
3072	local_pf(i,j,k) = rad_lw_cs_hr_av(k,j,i)
3073	ENDDO
3074	ENDDO
3075	ENDDO
3076	ENDIF
3077
3078	CASE ( 'rad_lw_hr' )
3079	IF ( av == 0 ) THEN
3080	DO i = nxlg, nxrg
3081	DO j = nysg, nyng
3082	DO k = nzb, nzt+1
3083	local_pf(i,j,k) = rad_lw_hr(k,j,i)
3084	ENDDO
3085	ENDDO
3086	ENDDO
3087	ELSE
3088	DO i = nxlg, nxrg
3089	DO j = nysg, nyng
3090	DO k = nzb, nzt+1
3091	local_pf(i,j,k) = rad_lw_hr_av(k,j,i)
3092	ENDDO
3093	ENDDO
3094	ENDDO
3095	ENDIF
3096
3097	CASE DEFAULT
3098	found = .FALSE.
3099
3100	END SELECT
3101
3102
3103	END SUBROUTINE radiation_data_output_3d
3104
3105	!------------------------------------------------------------------------------!
3106	!
3107	! Description:
3108	! ------------
3109	!> Subroutine defining masked data output
3110	!------------------------------------------------------------------------------!
3111	SUBROUTINE radiation_data_output_mask( av, variable, found, local_pf )
3112
3113	USE control_parameters
3114
3115	USE indices
3116
3117	USE kinds
3118
3119
3120	IMPLICIT NONE
3121
3122	CHARACTER (LEN=*) :: variable !<
3123
3124	INTEGER(iwp) :: av !<
3125	INTEGER(iwp) :: i !<
3126	INTEGER(iwp) :: j !<
3127	INTEGER(iwp) :: k !<
3128
3129	LOGICAL :: found !<
3130
3131	REAL(wp), &
3132	DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) :: &
3133	local_pf !<
3134
3135
3136	found = .TRUE.
3137
3138	SELECT CASE ( TRIM( variable ) )
3139
3140
3141	CASE ( 'rad_lw_in' )
3142	IF ( av == 0 ) THEN
3143	DO i = 1, mask_size_l(mid,1)
3144	DO j = 1, mask_size_l(mid,2)
3145	DO k = 1, mask_size_l(mid,3)
3146	local_pf(i,j,k) = rad_lw_in(mask_k(mid,k), &
3147	mask_j(mid,j),mask_i(mid,i))
3148	ENDDO
3149	ENDDO
3150	ENDDO
3151	ELSE
3152	DO i = 1, mask_size_l(mid,1)
3153	DO j = 1, mask_size_l(mid,2)
3154	DO k = 1, mask_size_l(mid,3)
3155	local_pf(i,j,k) = rad_lw_in_av(mask_k(mid,k), &
3156	mask_j(mid,j),mask_i(mid,i))
3157	ENDDO
3158	ENDDO
3159	ENDDO
3160	ENDIF
3161
3162	CASE ( 'rad_lw_out' )
3163	IF ( av == 0 ) THEN
3164	DO i = 1, mask_size_l(mid,1)
3165	DO j = 1, mask_size_l(mid,2)
3166	DO k = 1, mask_size_l(mid,3)
3167	local_pf(i,j,k) = rad_lw_out(mask_k(mid,k), &
3168	mask_j(mid,j),mask_i(mid,i))
3169	ENDDO
3170	ENDDO
3171	ENDDO
3172	ELSE
3173	DO i = 1, mask_size_l(mid,1)
3174	DO j = 1, mask_size_l(mid,2)
3175	DO k = 1, mask_size_l(mid,3)
3176	local_pf(i,j,k) = rad_lw_out_av(mask_k(mid,k), &
3177	mask_j(mid,j),mask_i(mid,i))
3178	ENDDO
3179	ENDDO
3180	ENDDO
3181	ENDIF
3182
3183	CASE ( 'rad_lw_cs_hr' )
3184	IF ( av == 0 ) THEN
3185	DO i = 1, mask_size_l(mid,1)
3186	DO j = 1, mask_size_l(mid,2)
3187	DO k = 1, mask_size_l(mid,3)
3188	local_pf(i,j,k) = rad_lw_cs_hr(mask_k(mid,k), &
3189	mask_j(mid,j),mask_i(mid,i))
3190	ENDDO
3191	ENDDO
3192	ENDDO
3193	ELSE
3194	DO i = 1, mask_size_l(mid,1)
3195	DO j = 1, mask_size_l(mid,2)
3196	DO k = 1, mask_size_l(mid,3)
3197	local_pf(i,j,k) = rad_lw_cs_hr_av(mask_k(mid,k), &
3198	mask_j(mid,j),mask_i(mid,i))
3199	ENDDO
3200	ENDDO
3201	ENDDO
3202	ENDIF
3203
3204	CASE ( 'rad_lw_hr' )
3205	IF ( av == 0 ) THEN
3206	DO i = 1, mask_size_l(mid,1)
3207	DO j = 1, mask_size_l(mid,2)
3208	DO k = 1, mask_size_l(mid,3)
3209	local_pf(i,j,k) = rad_lw_hr(mask_k(mid,k), &
3210	mask_j(mid,j),mask_i(mid,i))
3211	ENDDO
3212	ENDDO
3213	ENDDO
3214	ELSE
3215	DO i = 1, mask_size_l(mid,1)
3216	DO j = 1, mask_size_l(mid,2)
3217	DO k = 1, mask_size_l(mid,3)
3218	local_pf(i,j,k) = rad_lw_hr_av(mask_k(mid,k), &
3219	mask_j(mid,j),mask_i(mid,i))
3220	ENDDO
3221	ENDDO
3222	ENDDO
3223	ENDIF
3224
3225	CASE ( 'rad_sw_in' )
3226	IF ( av == 0 ) THEN
3227	DO i = 1, mask_size_l(mid,1)
3228	DO j = 1, mask_size_l(mid,2)
3229	DO k = 1, mask_size_l(mid,3)
3230	local_pf(i,j,k) = rad_sw_in(mask_k(mid,k), &
3231	mask_j(mid,j),mask_i(mid,i))
3232	ENDDO
3233	ENDDO
3234	ENDDO
3235	ELSE
3236	DO i = 1, mask_size_l(mid,1)
3237	DO j = 1, mask_size_l(mid,2)
3238	DO k = 1, mask_size_l(mid,3)
3239	local_pf(i,j,k) = rad_sw_in_av(mask_k(mid,k), &
3240	mask_j(mid,j),mask_i(mid,i))
3241	ENDDO
3242	ENDDO
3243	ENDDO
3244	ENDIF
3245
3246	CASE ( 'rad_sw_out' )
3247	IF ( av == 0 ) THEN
3248	DO i = 1, mask_size_l(mid,1)
3249	DO j = 1, mask_size_l(mid,2)
3250	DO k = 1, mask_size_l(mid,3)
3251	local_pf(i,j,k) = rad_sw_out(mask_k(mid,k), &
3252	mask_j(mid,j),mask_i(mid,i))
3253	ENDDO
3254	ENDDO
3255	ENDDO
3256	ELSE
3257	DO i = 1, mask_size_l(mid,1)
3258	DO j = 1, mask_size_l(mid,2)
3259	DO k = 1, mask_size_l(mid,3)
3260	local_pf(i,j,k) = rad_sw_out_av(mask_k(mid,k), &
3261	mask_j(mid,j),mask_i(mid,i))
3262	ENDDO
3263	ENDDO
3264	ENDDO
3265	ENDIF
3266
3267	CASE ( 'rad_sw_cs_hr' )
3268	IF ( av == 0 ) THEN
3269	DO i = 1, mask_size_l(mid,1)
3270	DO j = 1, mask_size_l(mid,2)
3271	DO k = 1, mask_size_l(mid,3)
3272	local_pf(i,j,k) = rad_sw_cs_hr(mask_k(mid,k), &
3273	mask_j(mid,j),mask_i(mid,i))
3274	ENDDO
3275	ENDDO
3276	ENDDO
3277	ELSE
3278	DO i = 1, mask_size_l(mid,1)
3279	DO j = 1, mask_size_l(mid,2)
3280	DO k = 1, mask_size_l(mid,3)
3281	local_pf(i,j,k) = rad_sw_cs_hr_av(mask_k(mid,k), &
3282	mask_j(mid,j),mask_i(mid,i))
3283	ENDDO
3284	ENDDO
3285	ENDDO
3286	ENDIF
3287
3288	CASE ( 'rad_sw_hr' )
3289	IF ( av == 0 ) THEN
3290	DO i = 1, mask_size_l(mid,1)
3291	DO j = 1, mask_size_l(mid,2)
3292	DO k = 1, mask_size_l(mid,3)
3293	local_pf(i,j,k) = rad_sw_hr(mask_k(mid,k), &
3294	mask_j(mid,j),mask_i(mid,i))
3295	ENDDO
3296	ENDDO
3297	ENDDO
3298	ELSE
3299	DO i = 1, mask_size_l(mid,1)
3300	DO j = 1, mask_size_l(mid,2)
3301	DO k = 1, mask_size_l(mid,3)
3302	local_pf(i,j,k) = rad_sw_hr_av(mask_k(mid,k), &
3303	mask_j(mid,j),mask_i(mid,i))
3304	ENDDO
3305	ENDDO
3306	ENDDO
3307	ENDIF
3308
3309	CASE DEFAULT
3310	found = .FALSE.
3311
3312	END SELECT
3313
3314
3315	END SUBROUTINE radiation_data_output_mask
3316
3317
3318	!------------------------------------------------------------------------------!
3319	!
3320	! Description:
3321	! ------------
3322	!> Subroutine defines masked output variables
3323	!------------------------------------------------------------------------------!
3324	SUBROUTINE radiation_last_actions
3325
3326
3327	USE control_parameters
3328
3329	USE kinds
3330
3331	IMPLICIT NONE
3332
3333	IF ( write_binary(1:4) == 'true' ) THEN
3334	IF ( ALLOCATED( rad_net ) ) THEN
3335	WRITE ( 14 ) 'rad_net '; WRITE ( 14 ) rad_net
3336	ENDIF
3337	IF ( ALLOCATED( rad_net_av ) ) THEN
3338	WRITE ( 14 ) 'rad_net_av '; WRITE ( 14 ) rad_net_av
3339	ENDIF
3340	IF ( ALLOCATED( rad_lw_in ) ) THEN
3341	WRITE ( 14 ) 'rad_lw_in '; WRITE ( 14 ) rad_lw_in
3342	ENDIF
3343	IF ( ALLOCATED( rad_lw_in_av ) ) THEN
3344	WRITE ( 14 ) 'rad_lw_in_av '; WRITE ( 14 ) rad_lw_in_av
3345	ENDIF
3346	IF ( ALLOCATED( rad_lw_out ) ) THEN
3347	WRITE ( 14 ) 'rad_lw_out '; WRITE ( 14 ) rad_lw_out
3348	ENDIF
3349	IF ( ALLOCATED( rad_lw_out_av ) ) THEN
3350	WRITE ( 14 ) 'rad_lw_out_av '; WRITE ( 14 ) rad_lw_out_av
3351	ENDIF
3352	IF ( ALLOCATED( rad_lw_out_change_0 ) ) THEN
3353	WRITE ( 14 ) 'rad_lw_out_change_0 '
3354	WRITE ( 14 ) rad_lw_out_change_0
3355	ENDIF
3356	IF ( ALLOCATED( rad_lw_cs_hr ) ) THEN
3357	WRITE ( 14 ) 'rad_lw_cs_hr '; WRITE ( 14 ) rad_lw_cs_hr
3358	ENDIF
3359	IF ( ALLOCATED( rad_lw_cs_hr_av ) ) THEN
3360	WRITE ( 14 ) 'rad_lw_cs_hr_av '; WRITE ( 14 ) rad_lw_cs_hr_av
3361	ENDIF
3362	IF ( ALLOCATED( rad_lw_hr ) ) THEN
3363	WRITE ( 14 ) 'rad_lw_hr '; WRITE ( 14 ) rad_lw_hr
3364	ENDIF
3365	IF ( ALLOCATED( rad_lw_hr_av ) ) THEN
3366	WRITE ( 14 ) 'rad_lw_hr_av '; WRITE ( 14 ) rad_lw_hr_av
3367	ENDIF
3368	IF ( ALLOCATED( rad_sw_in ) ) THEN
3369	WRITE ( 14 ) 'rad_sw_in '; WRITE ( 14 ) rad_sw_in
3370	ENDIF
3371	IF ( ALLOCATED( rad_sw_in_av ) ) THEN
3372	WRITE ( 14 ) 'rad_sw_in_av '; WRITE ( 14 ) rad_sw_in_av
3373	ENDIF
3374	IF ( ALLOCATED( rad_sw_out ) ) THEN
3375	WRITE ( 14 ) 'rad_sw_out '; WRITE ( 14 ) rad_sw_out
3376	ENDIF
3377	IF ( ALLOCATED( rad_sw_out_av ) ) THEN
3378	WRITE ( 14 ) 'rad_sw_out_av '; WRITE ( 14 ) rad_sw_out_av
3379	ENDIF
3380	IF ( ALLOCATED( rad_sw_cs_hr ) ) THEN
3381	WRITE ( 14 ) 'rad_sw_cs_hr '; WRITE ( 14 ) rad_sw_cs_hr
3382	ENDIF
3383	IF ( ALLOCATED( rad_sw_cs_hr_av ) ) THEN
3384	WRITE ( 14 ) 'rad_sw_cs_hr_av '; WRITE ( 14 ) rad_sw_cs_hr_av
3385	ENDIF
3386	IF ( ALLOCATED( rad_sw_hr ) ) THEN
3387	WRITE ( 14 ) 'rad_sw_hr '; WRITE ( 14 ) rad_sw_hr
3388	ENDIF
3389	IF ( ALLOCATED( rad_sw_hr_av ) ) THEN
3390	WRITE ( 14 ) 'rad_sw_hr_av '; WRITE ( 14 ) rad_sw_hr_av
3391	ENDIF
3392
3393	WRITE ( 14 ) '* end rad * '
3394
3395	ENDIF
3396
3397	END SUBROUTINE radiation_last_actions
3398
3399
3400	SUBROUTINE radiation_read_restart_data( i, nxlfa, nxl_on_file, nxrfa, nxr_on_file, &
3401	nynfa, nyn_on_file, nysfa, nys_on_file, &
3402	offset_xa, offset_ya, overlap_count, &
3403	tmp_2d, tmp_3d )
3404
3405
3406	USE control_parameters
3407
3408	USE indices
3409
3410	USE kinds
3411
3412	USE pegrid
3413
3414	IMPLICIT NONE
3415
3416	CHARACTER (LEN=20) :: field_char !<
3417
3418	INTEGER(iwp) :: i !<
3419	INTEGER(iwp) :: k !<
3420	INTEGER(iwp) :: nxlc !<
3421	INTEGER(iwp) :: nxlf !<
3422	INTEGER(iwp) :: nxl_on_file !<
3423	INTEGER(iwp) :: nxrc !<
3424	INTEGER(iwp) :: nxrf !<
3425	INTEGER(iwp) :: nxr_on_file !<
3426	INTEGER(iwp) :: nync !<
3427	INTEGER(iwp) :: nynf !<
3428	INTEGER(iwp) :: nyn_on_file !<
3429	INTEGER(iwp) :: nysc !<
3430	INTEGER(iwp) :: nysf !<
3431	INTEGER(iwp) :: nys_on_file !<
3432	INTEGER(iwp) :: overlap_count !<
3433
3434	INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: nxlfa !<
3435	INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: nxrfa !<
3436	INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: nynfa !<
3437	INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: nysfa !<
3438	INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: offset_xa !<
3439	INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) :: offset_ya !<
3440
3441	REAL(wp), &
3442	DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) ::&
3443	tmp_2d !<
3444
3445	REAL(wp), &
3446	DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) ::&
3447	tmp_3d !<
3448
3449	REAL(wp), &
3450	DIMENSION(0:0,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) ::&
3451	tmp_3d2 !<
3452
3453
3454
3455	IF ( initializing_actions == 'read_restart_data' ) THEN
3456	READ ( 13 ) field_char
3457
3458	DO WHILE ( TRIM( field_char ) /= '* end rad *' )
3459
3460	DO k = 1, overlap_count
3461
3462	nxlf = nxlfa(i,k)
3463	nxlc = nxlfa(i,k) + offset_xa(i,k)
3464	nxrf = nxrfa(i,k)
3465	nxrc = nxrfa(i,k) + offset_xa(i,k)
3466	nysf = nysfa(i,k)
3467	nysc = nysfa(i,k) + offset_ya(i,k)
3468	nynf = nynfa(i,k)
3469	nync = nynfa(i,k) + offset_ya(i,k)
3470
3471
3472	SELECT CASE ( TRIM( field_char ) )
3473
3474	CASE ( 'rad_net' )
3475	IF ( .NOT. ALLOCATED( rad_net ) ) THEN
3476	ALLOCATE( rad_net(nysg:nyng,nxlg:nxrg) )
3477	ENDIF
3478	IF ( k == 1 ) READ ( 13 ) tmp_2d
3479	rad_net(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3480	tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3481
3482	CASE ( 'rad_net_av' )
3483	IF ( .NOT. ALLOCATED( rad_net_av ) ) THEN
3484	ALLOCATE( rad_net_av(nysg:nyng,nxlg:nxrg) )
3485	ENDIF
3486	IF ( k == 1 ) READ ( 13 ) tmp_2d
3487	rad_net_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3488	tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3489	CASE ( 'rad_lw_in' )
3490	IF ( .NOT. ALLOCATED( rad_lw_in ) ) THEN
3491	IF ( radiation_scheme == 'clear-sky' .OR. &
3492	radiation_scheme == 'constant') THEN
3493	ALLOCATE( rad_lw_in(0:0,nysg:nyng,nxlg:nxrg) )
3494	ELSE
3495	ALLOCATE( rad_lw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3496	ENDIF
3497	ENDIF
3498	IF ( k == 1 ) THEN
3499	IF ( radiation_scheme == 'clear-sky' .OR. &
3500	radiation_scheme == 'constant') THEN
3501	READ ( 13 ) tmp_3d2
3502	rad_lw_in(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3503	tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3504	ELSE
3505	READ ( 13 ) tmp_3d
3506	rad_lw_in(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3507	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3508	ENDIF
3509	ENDIF
3510
3511	CASE ( 'rad_lw_in_av' )
3512	IF ( .NOT. ALLOCATED( rad_lw_in_av ) ) THEN
3513	IF ( radiation_scheme == 'clear-sky' .OR. &
3514	radiation_scheme == 'constant') THEN
3515	ALLOCATE( rad_lw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
3516	ELSE
3517	ALLOCATE( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3518	ENDIF
3519	ENDIF
3520	IF ( k == 1 ) THEN
3521	IF ( radiation_scheme == 'clear-sky' .OR. &
3522	radiation_scheme == 'constant') THEN
3523	READ ( 13 ) tmp_3d2
3524	rad_lw_in_av(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3525	tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3526	ELSE
3527	READ ( 13 ) tmp_3d
3528	rad_lw_in_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3529	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3530	ENDIF
3531	ENDIF
3532
3533	CASE ( 'rad_lw_out' )
3534	IF ( .NOT. ALLOCATED( rad_lw_out ) ) THEN
3535	IF ( radiation_scheme == 'clear-sky' .OR. &
3536	radiation_scheme == 'constant') THEN
3537	ALLOCATE( rad_lw_out(0:0,nysg:nyng,nxlg:nxrg) )
3538	ELSE
3539	ALLOCATE( rad_lw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3540	ENDIF
3541	ENDIF
3542	IF ( k == 1 ) THEN
3543	IF ( radiation_scheme == 'clear-sky' .OR. &
3544	radiation_scheme == 'constant') THEN
3545	READ ( 13 ) tmp_3d2
3546	rad_lw_out(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3547	tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3548	ELSE
3549	READ ( 13 ) tmp_3d
3550	rad_lw_out(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3551	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3552	ENDIF
3553	ENDIF
3554
3555	CASE ( 'rad_lw_out_av' )
3556	IF ( .NOT. ALLOCATED( rad_lw_out_av ) ) THEN
3557	IF ( radiation_scheme == 'clear-sky' .OR. &
3558	radiation_scheme == 'constant') THEN
3559	ALLOCATE( rad_lw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
3560	ELSE
3561	ALLOCATE( rad_lw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3562	ENDIF
3563	ENDIF
3564	IF ( k == 1 ) THEN
3565	IF ( radiation_scheme == 'clear-sky' .OR. &
3566	radiation_scheme == 'constant') THEN
3567	READ ( 13 ) tmp_3d2
3568	rad_lw_out_av(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3569	tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3570	ELSE
3571	READ ( 13 ) tmp_3d
3572	rad_lw_out_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3573	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3574	ENDIF
3575	ENDIF
3576
3577	CASE ( 'rad_lw_out_change_0' )
3578	IF ( .NOT. ALLOCATED( rad_lw_out_change_0 ) ) THEN
3579	ALLOCATE( rad_lw_out_change_0(nysg:nyng,nxlg:nxrg) )
3580	ENDIF
3581	IF ( k == 1 ) READ ( 13 ) tmp_2d
3582	rad_lw_out_change_0(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)&
3583	= tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3584
3585	CASE ( 'rad_lw_cs_hr' )
3586	IF ( .NOT. ALLOCATED( rad_lw_cs_hr ) ) THEN
3587	ALLOCATE( rad_lw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3588	ENDIF
3589	IF ( k == 1 ) READ ( 13 ) tmp_3d
3590	rad_lw_cs_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3591	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3592
3593	CASE ( 'rad_lw_cs_hr_av' )
3594	IF ( .NOT. ALLOCATED( rad_lw_cs_hr_av ) ) THEN
3595	ALLOCATE( rad_lw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3596	ENDIF
3597	IF ( k == 1 ) READ ( 13 ) tmp_3d
3598	rad_lw_cs_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3599	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3600
3601	CASE ( 'rad_lw_hr' )
3602	IF ( .NOT. ALLOCATED( rad_lw_hr ) ) THEN
3603	ALLOCATE( rad_lw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3604	ENDIF
3605	IF ( k == 1 ) READ ( 13 ) tmp_3d
3606	rad_lw_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3607	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3608
3609	CASE ( 'rad_lw_hr_av' )
3610	IF ( .NOT. ALLOCATED( rad_lw_hr_av ) ) THEN
3611	ALLOCATE( rad_lw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3612	ENDIF
3613	IF ( k == 1 ) READ ( 13 ) tmp_3d
3614	rad_lw_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3615	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3616
3617	CASE ( 'rad_sw_in' )
3618	IF ( .NOT. ALLOCATED( rad_sw_in ) ) THEN
3619	IF ( radiation_scheme == 'clear-sky' .OR. &
3620	radiation_scheme == 'constant') THEN
3621	ALLOCATE( rad_sw_in(0:0,nysg:nyng,nxlg:nxrg) )
3622	ELSE
3623	ALLOCATE( rad_sw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3624	ENDIF
3625	ENDIF
3626	IF ( k == 1 ) THEN
3627	IF ( radiation_scheme == 'clear-sky' .OR. &
3628	radiation_scheme == 'constant') THEN
3629	READ ( 13 ) tmp_3d2
3630	rad_sw_in(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3631	tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3632	ELSE
3633	READ ( 13 ) tmp_3d
3634	rad_sw_in(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3635	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3636	ENDIF
3637	ENDIF
3638
3639	CASE ( 'rad_sw_in_av' )
3640	IF ( .NOT. ALLOCATED( rad_sw_in_av ) ) THEN
3641	IF ( radiation_scheme == 'clear-sky' .OR. &
3642	radiation_scheme == 'constant') THEN
3643	ALLOCATE( rad_sw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
3644	ELSE
3645	ALLOCATE( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3646	ENDIF
3647	ENDIF
3648	IF ( k == 1 ) THEN
3649	IF ( radiation_scheme == 'clear-sky' .OR. &
3650	radiation_scheme == 'constant') THEN
3651	READ ( 13 ) tmp_3d2
3652	rad_sw_in_av(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3653	tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3654	ELSE
3655	READ ( 13 ) tmp_3d
3656	rad_sw_in_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3657	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3658	ENDIF
3659	ENDIF
3660
3661	CASE ( 'rad_sw_out' )
3662	IF ( .NOT. ALLOCATED( rad_sw_out ) ) THEN
3663	IF ( radiation_scheme == 'clear-sky' .OR. &
3664	radiation_scheme == 'constant') THEN
3665	ALLOCATE( rad_sw_out(0:0,nysg:nyng,nxlg:nxrg) )
3666	ELSE
3667	ALLOCATE( rad_sw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3668	ENDIF
3669	ENDIF
3670	IF ( k == 1 ) THEN
3671	IF ( radiation_scheme == 'clear-sky' .OR. &
3672	radiation_scheme == 'constant') THEN
3673	READ ( 13 ) tmp_3d2
3674	rad_sw_out(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3675	tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3676	ELSE
3677	READ ( 13 ) tmp_3d
3678	rad_sw_out(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3679	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3680	ENDIF
3681	ENDIF
3682
3683	CASE ( 'rad_sw_out_av' )
3684	IF ( .NOT. ALLOCATED( rad_sw_out_av ) ) THEN
3685	IF ( radiation_scheme == 'clear-sky' .OR. &
3686	radiation_scheme == 'constant') THEN
3687	ALLOCATE( rad_sw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
3688	ELSE
3689	ALLOCATE( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3690	ENDIF
3691	ENDIF
3692	IF ( k == 1 ) THEN
3693	IF ( radiation_scheme == 'clear-sky' .OR. &
3694	radiation_scheme == 'constant') THEN
3695	READ ( 13 ) tmp_3d2
3696	rad_sw_out_av(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3697	tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3698	ELSE
3699	READ ( 13 ) tmp_3d
3700	rad_sw_out_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3701	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3702	ENDIF
3703	ENDIF
3704
3705	CASE ( 'rad_sw_cs_hr' )
3706	IF ( .NOT. ALLOCATED( rad_sw_cs_hr ) ) THEN
3707	ALLOCATE( rad_sw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3708	ENDIF
3709	IF ( k == 1 ) READ ( 13 ) tmp_3d
3710	rad_sw_cs_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3711	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3712
3713	CASE ( 'rad_sw_cs_hr_av' )
3714	IF ( .NOT. ALLOCATED( rad_sw_cs_hr_av ) ) THEN
3715	ALLOCATE( rad_sw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3716	ENDIF
3717	IF ( k == 1 ) READ ( 13 ) tmp_3d
3718	rad_sw_cs_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3719	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3720
3721	CASE ( 'rad_sw_hr' )
3722	IF ( .NOT. ALLOCATED( rad_sw_hr ) ) THEN
3723	ALLOCATE( rad_sw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3724	ENDIF
3725	IF ( k == 1 ) READ ( 13 ) tmp_3d
3726	rad_sw_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3727	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3728
3729	CASE ( 'rad_sw_hr_av' )
3730	IF ( .NOT. ALLOCATED( rad_sw_hr_av ) ) THEN
3731	ALLOCATE( rad_sw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3732	ENDIF
3733	IF ( k == 1 ) READ ( 13 ) tmp_3d
3734	rad_lw_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3735	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3736
3737	CASE DEFAULT
3738	WRITE( message_string, * ) 'unknown variable named "', &
3739	TRIM( field_char ), '" found in', &
3740	'&data from prior run on PE ', myid
3741	CALL message( 'radiation_read_restart_data', 'PA0441', 1, 2, 0, 6, &
3742	0 )
3743
3744	END SELECT
3745
3746	ENDDO
3747
3748	READ ( 13 ) field_char
3749
3750	ENDDO
3751	ENDIF
3752
3753	END SUBROUTINE radiation_read_restart_data
3754
3755
3756	END MODULE radiation_model_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |