Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |