Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |