Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |