Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |