Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |