Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |