Home

Context Navigation

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

Last change on this file since 2542 was 2512, checked in by raasch, 7 years ago
upper bounds of cross section and 3d output changed from nx+1,ny+1 to nx,ny; no output if redundant ghost layer data to NetCDF files
Property svn:keywords set to `Id`
File size: 148.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |