Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |