Home

Context Navigation

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

Last change on this file since 2502 was 2328, checked in by maronga, 7 years ago
some improvements in land surface model
Property svn:keywords set to `Id`
File size: 146.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |