Home

Context Navigation

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

Last change on this file since 2209 was 2201, checked in by suehring, 7 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 142.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |