Home

Context Navigation

radiation_model_mod.f90 @ 3825

Last change on this file since 3825 was 3825, checked in by raasch, 6 years ago

unused variable removed

Property svn:keywords set to Id

Property svn:mergeinfo set to (toggle deleted branches)

/palm/branches/chemistry/SOURCE/radiation_model_mod.f90	2047-3190,3218-3297
/palm/branches/forwind/SOURCE/radiation_model_mod.f90	1564-1913
/palm/branches/mosaik_M2/radiation_model_mod.f90	2360-3471
/palm/branches/palm4u/SOURCE/radiation_model_mod.f90	2540-2692
/palm/branches/radiation/SOURCE/radiation_model_mod.f90	2081-3493
/palm/branches/rans/SOURCE/radiation_model_mod.f90	2078-3128
/palm/branches/resler/SOURCE/radiation_model_mod.f90	2023-3605
/palm/branches/salsa/SOURCE/radiation_model_mod.f90	2503-3460
/palm/branches/fricke/SOURCE/radiation_model_mod.f90	942-977
/palm/branches/hoffmann/SOURCE/radiation_model_mod.f90	989-1052
/palm/branches/letzel/masked_output/SOURCE/radiation_model_mod.f90	296-409
/palm/branches/suehring/radiation_model_mod.f90	423-666

File size: 496.8 KB

Line
1	!> @file radiation_model_mod.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of the PALM model system.
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 2015-2018 Institute of Computer Science of the
18	! Czech Academy of Sciences, Prague
19	! Copyright 2015-2018 Czech Technical University in Prague
20	! Copyright 1997-2019 Leibniz Universitaet Hannover
21	!------------------------------------------------------------------------------!
22	!
23	! Current revisions:
24	! ------------------
25	!
26	!
27	! Former revisions:
28	! -----------------
29	! $Id: radiation_model_mod.f90 3825 2019-03-27 16:01:25Z raasch $
30	<<<<<<< .mine
31	! unused variable removed
32	!
33	! 3814 2019-03-26 08:40:31Z pavelkrc
34	\|\|\|\|\|\|\| .r3822
35	=======
36	! Change zenith(0:0) and others to scalar.
37	! Code review.
38	!
39	! 3814 2019-03-26 08:40:31Z pavelkrc
40	>>>>>>> .r3824
41	! Rename exported nzu, nzp and related variables due to name conflict
42	!
43	! 3771 2019-02-28 12:19:33Z raasch
44	! rrtmg preprocessor for directives moved/added, save attribute added to temporary
45	! pointers to avoid compiler warnings about outlived pointer targets,
46	! statement added to avoid compiler warning about unused variable
47	!
48	! 3769 2019-02-28 10:16:49Z moh.hefny
49	! removed unused variables and subroutine radiation_radflux_gridbox
50	!
51	! 3767 2019-02-27 08:18:02Z raasch
52	! unused variable for file index removed from rrd-subroutines parameter list
53	!
54	! 3760 2019-02-21 18:47:35Z moh.hefny
55	! Bugfix: initialized simulated_time before calculating solar position
56	! to enable restart option with reading in SVF from file(s).
57	!
58	! 3754 2019-02-19 17:02:26Z kanani
59	! (resler, pavelkrc)
60	! Bugfixes: add further required MRT factors to read/write_svf,
61	! fix for aggregating view factors to eliminate local noise in reflected
62	! irradiance at mutually close surfaces (corners, presence of trees) in the
63	! angular discretization scheme.
64	!
65	! 3752 2019-02-19 09:37:22Z resler
66	! added read/write number of MRT factors to the respective routines
67	!
68	! 3705 2019-01-29 19:56:39Z suehring
69	! Make variables that are sampled in virtual measurement module public
70	!
71	! 3704 2019-01-29 19:51:41Z suehring
72	! Some interface calls moved to module_interface + cleanup
73	!
74	! 3667 2019-01-10 14:26:24Z schwenkel
75	! Modified check for rrtmg input files
76	!
77	! 3655 2019-01-07 16:51:22Z knoop
78	! nopointer option removed
79	!
80	! 3633 2018-12-17 16:17:57Z schwenkel
81	! Include check for rrtmg files
82	!
83	! 3630 2018-12-17 11:04:17Z knoop
84	! - fix initialization of date and time after calling zenith
85	! - fix a bug in radiation_solar_pos
86	!
87	! 3616 2018-12-10 09:44:36Z Salim
88	! fix manipulation of time variables in radiation_presimulate_solar_pos
89	!
90	! 3608 2018-12-07 12:59:57Z suehring $
91	! Bugfix radiation output
92	!
93	! 3607 2018-12-07 11:56:58Z suehring
94	! Output of radiation-related quantities migrated to radiation_model_mod.
95	!
96	! 3589 2018-11-30 15:09:51Z suehring
97	! Remove erroneous UTF encoding
98	!
99	! 3572 2018-11-28 11:40:28Z suehring
100	! Add filling the short- and longwave radiation flux arrays (e.g. diffuse,
101	! direct, reflected, resedual) for all surfaces. This is required to surface
102	! outputs in suface_output_mod. (M. Salim)
103	!
104	! 3571 2018-11-28 09:24:03Z moh.hefny
105	! Add an epsilon value to compare values in if statement to fix possible
106	! precsion related errors in raytrace routines.
107	!
108	! 3524 2018-11-14 13:36:44Z raasch
109	! missing cpp-directives added
110	!
111	! 3495 2018-11-06 15:22:17Z kanani
112	! Resort control_parameters ONLY list,
113	! From branch radiation@3491 moh.hefny:
114	! bugfix in calculating the apparent solar positions by updating
115	! the simulated time so that the actual time is correct.
116	!
117	! 3464 2018-10-30 18:08:55Z kanani
118	! From branch resler@3462, pavelkrc:
119	! add MRT shaping function for human
120	!
121	! 3449 2018-10-29 19:36:56Z suehring
122	! New RTM version 3.0: (Pavel Krc, Jaroslav Resler, ICS, Prague)
123	! - Interaction of plant canopy with LW radiation
124	! - Transpiration from resolved plant canopy dependent on radiation
125	! called from RTM
126	!
127	!
128	! 3435 2018-10-26 18:25:44Z gronemeier
129	! - workaround: return unit=illegal in check_data_output for certain variables
130	! when check called from init_masks
131	! - Use pointer in masked output to reduce code redundancies
132	! - Add terrain-following masked output
133	!
134	! 3424 2018-10-25 07:29:10Z gronemeier
135	! bugfix: add rad_lw_in, rad_lw_out, rad_sw_out to radiation_check_data_output
136	!
137	! 3378 2018-10-19 12:34:59Z kanani
138	! merge from radiation branch (r3362) into trunk
139	! (moh.hefny):
140	! - removed read/write_svf_on_init and read_dist_max_svf (not used anymore)
141	! - bugfix nzut > nzpt in calculating maxboxes
142	!
143	! 3372 2018-10-18 14:03:19Z raasch
144	! bugfix: kind type of 2nd argument of mpi_win_allocate changed, misplaced
145	! __parallel directive
146	!
147	! 3351 2018-10-15 18:40:42Z suehring
148	! Do not overwrite values of spectral and broadband albedo during initialization
149	! if they are already initialized in the urban-surface model via ASCII input.
150	!
151	! 3337 2018-10-12 15:17:09Z kanani
152	! - New RTM version 2.9: (Pavel Krc, Jaroslav Resler, ICS, Prague)
153	! added calculation of the MRT inside the RTM module
154	! MRT fluxes are consequently used in the new biometeorology module
155	! for calculation of biological indices (MRT, PET)
156	! Fixes of v. 2.5 and SVN trunk:
157	! - proper initialization of rad_net_l
158	! - make arrays nsurfs and surfstart TARGET to prevent some MPI problems
159	! - initialization of arrays used in MPI one-sided operation as 1-D arrays
160	! to prevent problems with some MPI/compiler combinations
161	! - fix indexing of target displacement in subroutine request_itarget to
162	! consider nzub
163	! - fix LAD dimmension range in PCB calculation
164	! - check ierr in all MPI calls
165	! - use proper per-gridbox sky and diffuse irradiance
166	! - fix shading for reflected irradiance
167	! - clear away the residuals of "atmospheric surfaces" implementation
168	! - fix rounding bug in raytrace_2d introduced in SVN trunk
169	! - New RTM version 2.5: (Pavel Krc, Jaroslav Resler, ICS, Prague)
170	! can use angular discretization for all SVF
171	! (i.e. reflected and emitted radiation in addition to direct and diffuse),
172	! allowing for much better scaling wih high resoltion and/or complex terrain
173	! - Unite array grow factors
174	! - Fix slightly shifted terrain height in raytrace_2d
175	! - Use more efficient MPI_Win_allocate for reverse gridsurf index
176	! - Fix random MPI RMA bugs on Intel compilers
177	! - Fix approx. double plant canopy sink values for reflected radiation
178	! - Fix mostly missing plant canopy sinks for direct radiation
179	! - Fix discretization errors for plant canopy sink in diffuse radiation
180	! - Fix rounding errors in raytrace_2d
181	!
182	! 3274 2018-09-24 15:42:55Z knoop
183	! Modularization of all bulk cloud physics code components
184	!
185	! 3272 2018-09-24 10:16:32Z suehring
186	! - split direct and diffusion shortwave radiation using RRTMG rather than using
187	! calc_diffusion_radiation, in case of RRTMG
188	! - removed the namelist variable split_diffusion_radiation. Now splitting depends
189	! on the choise of radiation radiation scheme
190	! - removed calculating the rdiation flux for surfaces at the radiation scheme
191	! in case of using RTM since it will be calculated anyway in the radiation
192	! interaction routine.
193	! - set SW radiation flux for surfaces to zero at night in case of no RTM is used
194	! - precalculate the unit vector yxdir of ray direction to avoid the temporarly
195	! array allocation during the subroutine call
196	! - fixed a bug in calculating the max number of boxes ray can cross in the domain
197	!
198	! 3264 2018-09-20 13:54:11Z moh.hefny
199	! Bugfix in raytrace_2d calls
200	!
201	! 3248 2018-09-14 09:42:06Z sward
202	! Minor formating changes
203	!
204	! 3246 2018-09-13 15:14:50Z sward
205	! Added error handling for input namelist via parin_fail_message
206	!
207	! 3241 2018-09-12 15:02:00Z raasch
208	! unused variables removed or commented
209	!
210	! 3233 2018-09-07 13:21:24Z schwenkel
211	! Adapted for the use of cloud_droplets
212	!
213	! 3230 2018-09-05 09:29:05Z schwenkel
214	! Bugfix in radiation_constant_surf: changed (10.0 - emissivity_urb) to
215	! (1.0 - emissivity_urb)
216	!
217	! 3226 2018-08-31 12:27:09Z suehring
218	! Bugfixes in calculation of sky-view factors and canopy-sink factors.
219	!
220	! 3186 2018-07-30 17:07:14Z suehring
221	! Remove print statement
222	!
223	! 3180 2018-07-27 11:00:56Z suehring
224	! Revise concept for calculation of effective radiative temperature and mapping
225	! of radiative heating
226	!
227	! 3175 2018-07-26 14:07:38Z suehring
228	! Bugfix for commit 3172
229	!
230	! 3173 2018-07-26 12:55:23Z suehring
231	! Revise output of surface radiation quantities in case of overhanging
232	! structures
233	!
234	! 3172 2018-07-26 12:06:06Z suehring
235	! Bugfixes:
236	! - temporal work-around for calculation of effective radiative surface
237	! temperature
238	! - prevent positive solar radiation during nighttime
239	!
240	! 3170 2018-07-25 15:19:37Z suehring
241	! Bugfix, map signle-column radiation forcing profiles on top of any topography
242	!
243	! 3156 2018-07-19 16:30:54Z knoop
244	! Bugfix: replaced usage of the pt array with the surf%pt_surface array
245	!
246	! 3137 2018-07-17 06:44:21Z maronga
247	! String length for trace_names fixed
248	!
249	! 3127 2018-07-15 08:01:25Z maronga
250	! A few pavement parameters updated.
251	!
252	! 3123 2018-07-12 16:21:53Z suehring
253	! Correct working precision for INTEGER number
254	!
255	! 3122 2018-07-11 21:46:41Z maronga
256	! Bugfix: maximum distance for raytracing was set to -999 m by default,
257	! effectively switching off all surface reflections when max_raytracing_dist
258	! was not explicitly set in namelist
259	!
260	! 3117 2018-07-11 09:59:11Z maronga
261	! Bugfix: water vapor was not transfered to RRTMG when bulk_cloud_model = .F.
262	! Bugfix: changed the calculation of RRTMG boundary conditions (Mohamed Salim)
263	! Bugfix: dry residual atmosphere is replaced by standard RRTMG atmosphere
264	!
265	! 3116 2018-07-10 14:31:58Z suehring
266	! Output of long/shortwave radiation at surface
267	!
268	! 3107 2018-07-06 15:55:51Z suehring
269	! Bugfix, missing index for dz
270	!
271	! 3066 2018-06-12 08:55:55Z Giersch
272	! Error message revised
273	!
274	! 3065 2018-06-12 07:03:02Z Giersch
275	! dz was replaced by dz(1), error message concerning vertical stretching was
276	! added
277	!
278	! 3049 2018-05-29 13:52:36Z Giersch
279	! Error messages revised
280	!
281	! 3045 2018-05-28 07:55:41Z Giersch
282	! Error message revised
283	!
284	! 3026 2018-05-22 10:30:53Z schwenkel
285	! Changed the name specific humidity to mixing ratio, since we are computing
286	! mixing ratios.
287	!
288	! 3016 2018-05-09 10:53:37Z Giersch
289	! Revised structure of reading svf data according to PALM coding standard:
290	! svf_code_field/len and fsvf removed, error messages PA0493 and PA0494 added,
291	! allocation status of output arrays checked.
292	!
293	! 3014 2018-05-09 08:42:38Z maronga
294	! Introduced plant canopy height similar to urban canopy height to limit
295	! the memory requirement to allocate lad.
296	! Deactivated automatic setting of minimum raytracing distance.
297	!
298	! 3004 2018-04-27 12:33:25Z Giersch
299	! Further allocation checks implemented (averaged data will be assigned to fill
300	! values if no allocation happened so far)
301	!
302	! 2995 2018-04-19 12:13:16Z Giersch
303	! IF-statement in radiation_init removed so that the calculation of radiative
304	! fluxes at model start is done in any case, bugfix in
305	! radiation_presimulate_solar_pos (end_time is the sum of end_time and the
306	! spinup_time specified in the p3d_file ), list of variables/fields that have
307	! to be written out or read in case of restarts has been extended
308	!
309	! 2977 2018-04-17 10:27:57Z kanani
310	! Implement changes from branch radiation (r2948-2971) with minor modifications,
311	! plus some formatting.
312	! (moh.hefny):
313	! - replaced plant_canopy by npcbl to check tree existence to avoid weird
314	! allocation of related arrays (after domain decomposition some domains
315	! contains no trees although plant_canopy (global parameter) is still TRUE).
316	! - added a namelist parameter to force RTM settings
317	! - enabled the option to switch radiation reflections off
318	! - renamed surf_reflections to surface_reflections
319	! - removed average_radiation flag from the namelist (now it is implicitly set
320	! in init_3d_model according to RTM)
321	! - edited read and write sky view factors and CSF routines to account for
322	! the sub-domains which may not contain any of them
323	!
324	! 2967 2018-04-13 11:22:08Z raasch
325	! bugfix: missing parallel cpp-directives added
326	!
327	! 2964 2018-04-12 16:04:03Z Giersch
328	! Error message PA0491 has been introduced which could be previously found in
329	! check_open. The variable numprocs_previous_run is only known in case of
330	! initializing_actions == read_restart_data
331	!
332	! 2963 2018-04-12 14:47:44Z suehring
333	! - Introduce index for vegetation/wall, pavement/green-wall and water/window
334	! surfaces, for clearer access of surface fraction, albedo, emissivity, etc. .
335	! - Minor bugfix in initialization of albedo for window surfaces
336	!
337	! 2944 2018-04-03 16:20:18Z suehring
338	! Fixed bad commit
339	!
340	! 2943 2018-04-03 16:17:10Z suehring
341	! No read of nsurfl from SVF file since it is calculated in
342	! radiation_interaction_init,
343	! allocation of arrays in radiation_read_svf only if not yet allocated,
344	! update of 2920 revision comment.
345	!
346	! 2932 2018-03-26 09:39:22Z maronga
347	! renamed radiation_par to radiation_parameters
348	!
349	! 2930 2018-03-23 16:30:46Z suehring
350	! Remove default surfaces from radiation model, does not make much sense to
351	! apply radiation model without energy-balance solvers; Further, add check for
352	! this.
353	!
354	! 2920 2018-03-22 11:22:01Z kanani
355	! - Bugfix: Initialize pcbl array (=-1)
356	! RTM version 2.0 (Jaroslav Resler, Pavel Krc, Mohamed Salim):
357	! - new major version of radiation interactions
358	! - substantially enhanced performance and scalability
359	! - processing of direct and diffuse solar radiation separated from reflected
360	! radiation, removed virtual surfaces
361	! - new type of sky discretization by azimuth and elevation angles
362	! - diffuse radiation processed cumulatively using sky view factor
363	! - used precalculated apparent solar positions for direct irradiance
364	! - added new 2D raytracing process for processing whole vertical column at once
365	! to increase memory efficiency and decrease number of MPI RMA operations
366	! - enabled limiting the number of view factors between surfaces by the distance
367	! and value
368	! - fixing issues induced by transferring radiation interactions from
369	! urban_surface_mod to radiation_mod
370	! - bugfixes and other minor enhancements
371	!
372	! 2906 2018-03-19 08:56:40Z Giersch
373	! NAMELIST paramter read/write_svf_on_init have been removed, functions
374	! check_open and close_file are used now for opening/closing files related to
375	! svf data, adjusted unit number and error numbers
376	!
377	! 2894 2018-03-15 09:17:58Z Giersch
378	! Calculations of the index range of the subdomain on file which overlaps with
379	! the current subdomain are already done in read_restart_data_mod
380	! radiation_read_restart_data was renamed to radiation_rrd_local and
381	! radiation_last_actions was renamed to radiation_wrd_local, variable named
382	! found has been introduced for checking if restart data was found, reading
383	! of restart strings has been moved completely to read_restart_data_mod,
384	! radiation_rrd_local is already inside the overlap loop programmed in
385	! read_restart_data_mod, the marker * end rad * is not necessary anymore,
386	! strings and their respective lengths are written out and read now in case of
387	! restart runs to get rid of prescribed character lengths (Giersch)
388	!
389	! 2809 2018-02-15 09:55:58Z suehring
390	! Bugfix for gfortran: Replace the function C_SIZEOF with STORAGE_SIZE
391	!
392	! 2753 2018-01-16 14:16:49Z suehring
393	! Tile approach for spectral albedo implemented.
394	!
395	! 2746 2018-01-15 12:06:04Z suehring
396	! Move flag plant canopy to modules
397	!
398	! 2724 2018-01-05 12:12:38Z maronga
399	! Set default of average_radiation to .FALSE.
400	!
401	! 2723 2018-01-05 09:27:03Z maronga
402	! Bugfix in calculation of rad_lw_out (clear-sky). First grid level was used
403	! instead of the surface value
404	!
405	! 2718 2018-01-02 08:49:38Z maronga
406	! Corrected "Former revisions" section
407	!
408	! 2707 2017-12-18 18:34:46Z suehring
409	! Changes from last commit documented
410	!
411	! 2706 2017-12-18 18:33:49Z suehring
412	! Bugfix, in average radiation case calculate exner function before using it.
413	!
414	! 2701 2017-12-15 15:40:50Z suehring
415	! Changes from last commit documented
416	!
417	! 2698 2017-12-14 18:46:24Z suehring
418	! Bugfix in get_topography_top_index
419	!
420	! 2696 2017-12-14 17:12:51Z kanani
421	! - Change in file header (GPL part)
422	! - Improved reading/writing of SVF from/to file (BM)
423	! - Bugfixes concerning RRTMG as well as average_radiation options (M. Salim)
424	! - Revised initialization of surface albedo and some minor bugfixes (MS)
425	! - Update net radiation after running radiation interaction routine (MS)
426	! - Revisions from M Salim included
427	! - Adjustment to topography and surface structure (MS)
428	! - Initialization of albedo and surface emissivity via input file (MS)
429	! - albedo_pars extended (MS)
430	!
431	! 2604 2017-11-06 13:29:00Z schwenkel
432	! bugfix for calculation of effective radius using morrison microphysics
433	!
434	! 2601 2017-11-02 16:22:46Z scharf
435	! added emissivity to namelist
436	!
437	! 2575 2017-10-24 09:57:58Z maronga
438	! Bugfix: calculation of shortwave and longwave albedos for RRTMG swapped
439	!
440	! 2547 2017-10-16 12:41:56Z schwenkel
441	! extended by cloud_droplets option, minor bugfix and correct calculation of
442	! cloud droplet number concentration
443	!
444	! 2544 2017-10-13 18:09:32Z maronga
445	! Moved date and time quantitis to separate module date_and_time_mod
446	!
447	! 2512 2017-10-04 08:26:59Z raasch
448	! upper bounds of cross section and 3d output changed from nx+1,ny+1 to nx,ny
449	! no output of ghost layer data
450	!
451	! 2504 2017-09-27 10:36:13Z maronga
452	! Updates pavement types and albedo parameters
453	!
454	! 2328 2017-08-03 12:34:22Z maronga
455	! Emissivity can now be set individually for each pixel.
456	! Albedo type can be inferred from land surface model.
457	! Added default albedo type for bare soil
458	!
459	! 2318 2017-07-20 17:27:44Z suehring
460	! Get topography top index via Function call
461	!
462	! 2317 2017-07-20 17:27:19Z suehring
463	! Improved syntax layout
464	!
465	! 2298 2017-06-29 09:28:18Z raasch
466	! type of write_binary changed from CHARACTER to LOGICAL
467	!
468	! 2296 2017-06-28 07:53:56Z maronga
469	! Added output of rad_sw_out for radiation_scheme = 'constant'
470	!
471	! 2270 2017-06-09 12:18:47Z maronga
472	! Numbering changed (2 timeseries removed)
473	!
474	! 2249 2017-06-06 13:58:01Z sward
475	! Allow for RRTMG runs without humidity/cloud physics
476	!
477	! 2248 2017-06-06 13:52:54Z sward
478	! Error no changed
479	!
480	! 2233 2017-05-30 18:08:54Z suehring
481	!
482	! 2232 2017-05-30 17:47:52Z suehring
483	! Adjustments to new topography concept
484	! Bugfix in read restart
485	!
486	! 2200 2017-04-11 11:37:51Z suehring
487	! Bugfix in call of exchange_horiz_2d and read restart data
488	!
489	! 2163 2017-03-01 13:23:15Z schwenkel
490	! Bugfix in radiation_check_data_output
491	!
492	! 2157 2017-02-22 15:10:35Z suehring
493	! Bugfix in read_restart data
494	!
495	! 2011 2016-09-19 17:29:57Z kanani
496	! Removed CALL of auxiliary SUBROUTINE get_usm_info,
497	! flag urban_surface is now defined in module control_parameters.
498	!
499	! 2007 2016-08-24 15:47:17Z kanani
500	! Added calculation of solar directional vector for new urban surface
501	! model,
502	! accounted for urban_surface model in radiation_check_parameters,
503	! correction of comments for zenith angle.
504	!
505	! 2000 2016-08-20 18:09:15Z knoop
506	! Forced header and separation lines into 80 columns
507	!
508	! 1976 2016-07-27 13:28:04Z maronga
509	! Output of 2D/3D/masked data is now directly done within this module. The
510	! radiation schemes have been simplified for better usability so that
511	! rad_lw_in, rad_lw_out, rad_sw_in, and rad_sw_out are available independent of
512	! the radiation code used.
513	!
514	! 1856 2016-04-13 12:56:17Z maronga
515	! Bugfix: allocation of rad_lw_out for radiation_scheme = 'clear-sky'
516	!
517	! 1853 2016-04-11 09:00:35Z maronga
518	! Added routine for radiation_scheme = constant.
519	!
520	! 1849 2016-04-08 11:33:18Z hoffmann
521	! Adapted for modularization of microphysics
522	!
523	! 1826 2016-04-07 12:01:39Z maronga
524	! Further modularization.
525	!
526	! 1788 2016-03-10 11:01:04Z maronga
527	! Added new albedo class for pavements / roads.
528	!
529	! 1783 2016-03-06 18:36:17Z raasch
530	! palm-netcdf-module removed in order to avoid a circular module dependency,
531	! netcdf-variables moved to netcdf-module, new routine netcdf_handle_error_rad
532	! added
533	!
534	! 1757 2016-02-22 15:49:32Z maronga
535	! Added parameter unscheduled_radiation_calls. Bugfix: interpolation of sounding
536	! profiles for pressure and temperature above the LES domain.
537	!
538	! 1709 2015-11-04 14:47:01Z maronga
539	! Bugfix: set initial value for rrtm_lwuflx_dt to zero, small formatting
540	! corrections
541	!
542	! 1701 2015-11-02 07:43:04Z maronga
543	! Bugfixes: wrong index for output of timeseries, setting of nz_snd_end
544	!
545	! 1691 2015-10-26 16:17:44Z maronga
546	! Added option for spin-up runs without radiation (skip_time_do_radiation). Bugfix
547	! in calculation of pressure profiles. Bugfix in calculation of trace gas profiles.
548	! Added output of radiative heating rates.
549	!
550	! 1682 2015-10-07 23:56:08Z knoop
551	! Code annotations made doxygen readable
552	!
553	! 1606 2015-06-29 10:43:37Z maronga
554	! Added preprocessor directive __netcdf to allow for compiling without netCDF.
555	! Note, however, that RRTMG cannot be used without netCDF.
556	!
557	! 1590 2015-05-08 13:56:27Z maronga
558	! Bugfix: definition of character strings requires same length for all elements
559	!
560	! 1587 2015-05-04 14:19:01Z maronga
561	! Added albedo class for snow
562	!
563	! 1585 2015-04-30 07:05:52Z maronga
564	! Added support for RRTMG
565	!
566	! 1571 2015-03-12 16:12:49Z maronga
567	! Added missing KIND attribute. Removed upper-case variable names
568	!
569	! 1551 2015-03-03 14:18:16Z maronga
570	! Added support for data output. Various variables have been renamed. Added
571	! interface for different radiation schemes (currently: clear-sky, constant, and
572	! RRTM (not yet implemented).
573	!
574	! 1496 2014-12-02 17:25:50Z maronga
575	! Initial revision
576	!
577	!
578	! Description:
579	! ------------
580	!> Radiation models and interfaces
581	!> @todo Replace dz(1) appropriatly to account for grid stretching
582	!> @todo move variable definitions used in radiation_init only to the subroutine
583	!> as they are no longer required after initialization.
584	!> @todo Output of full column vertical profiles used in RRTMG
585	!> @todo Output of other rrtm arrays (such as volume mixing ratios)
586	!> @todo Check for mis-used NINT() calls in raytrace_2d
587	!> RESULT: Original was correct (carefully verified formula), the change
588	!> to INT broke raytracing -- P. Krc
589	!> @todo Optimize radiation_tendency routines
590	!>
591	!> @note Many variables have a leading dummy dimension (0:0) in order to
592	!> match the assume-size shape expected by the RRTMG model.
593	!------------------------------------------------------------------------------!
594	MODULE radiation_model_mod
595
596	USE arrays_3d, &
597	ONLY: dzw, hyp, nc, pt, p, q, ql, u, v, w, zu, zw, exner, d_exner
598
599	USE basic_constants_and_equations_mod, &
600	ONLY: c_p, g, lv_d_cp, l_v, pi, r_d, rho_l, solar_constant, sigma_sb, &
601	barometric_formula
602
603	USE calc_mean_profile_mod, &
604	ONLY: calc_mean_profile
605
606	USE control_parameters, &
607	ONLY: cloud_droplets, coupling_char, dz, dt_spinup, end_time, &
608	humidity, &
609	initializing_actions, io_blocks, io_group, &
610	land_surface, large_scale_forcing, &
611	latitude, longitude, lsf_surf, &
612	message_string, plant_canopy, pt_surface, &
613	rho_surface, simulated_time, spinup_time, surface_pressure, &
614	read_svf, write_svf, &
615	time_since_reference_point, urban_surface, varnamelength
616
617	USE cpulog, &
618	ONLY: cpu_log, log_point, log_point_s
619
620	USE grid_variables, &
621	ONLY: ddx, ddy, dx, dy
622
623	USE date_and_time_mod, &
624	ONLY: calc_date_and_time, d_hours_day, d_seconds_hour, d_seconds_year,&
625	day_of_year, d_seconds_year, day_of_month, day_of_year_init, &
626	init_date_and_time, month_of_year, time_utc_init, time_utc
627
628	USE indices, &
629	ONLY: nnx, nny, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, &
630	nzb, nzt
631
632	USE, INTRINSIC :: iso_c_binding
633
634	USE kinds
635
636	USE bulk_cloud_model_mod, &
637	ONLY: bulk_cloud_model, microphysics_morrison, na_init, nc_const, sigma_gc
638
639	#if defined ( __netcdf )
640	USE NETCDF
641	#endif
642
643	USE netcdf_data_input_mod, &
644	ONLY: albedo_type_f, albedo_pars_f, building_type_f, pavement_type_f, &
645	vegetation_type_f, water_type_f
646
647	USE plant_canopy_model_mod, &
648	ONLY: lad_s, pc_heating_rate, pc_transpiration_rate, pc_latent_rate, &
649	plant_canopy_transpiration, pcm_calc_transpiration_rate
650
651	USE pegrid
652
653	#if defined ( __rrtmg )
654	USE parrrsw, &
655	ONLY: naerec, nbndsw
656
657	USE parrrtm, &
658	ONLY: nbndlw
659
660	USE rrtmg_lw_init, &
661	ONLY: rrtmg_lw_ini
662
663	USE rrtmg_sw_init, &
664	ONLY: rrtmg_sw_ini
665
666	USE rrtmg_lw_rad, &
667	ONLY: rrtmg_lw
668
669	USE rrtmg_sw_rad, &
670	ONLY: rrtmg_sw
671	#endif
672	USE statistics, &
673	ONLY: hom
674
675	USE surface_mod, &
676	ONLY: get_topography_top_index, get_topography_top_index_ji, &
677	ind_pav_green, ind_veg_wall, ind_wat_win, &
678	surf_lsm_h, surf_lsm_v, surf_type, surf_usm_h, surf_usm_v, &
679	vertical_surfaces_exist
680
681	IMPLICIT NONE
682
683	CHARACTER(10) :: radiation_scheme = 'clear-sky' ! 'constant', 'clear-sky', or 'rrtmg'
684
685	!
686	!-- Predefined Land surface classes (albedo_type) after Briegleb (1992)
687	CHARACTER(37), DIMENSION(0:33), PARAMETER :: albedo_type_name = (/ &
688	'user defined ', & ! 0
689	'ocean ', & ! 1
690	'mixed farming, tall grassland ', & ! 2
691	'tall/medium grassland ', & ! 3
692	'evergreen shrubland ', & ! 4
693	'short grassland/meadow/shrubland ', & ! 5
694	'evergreen needleleaf forest ', & ! 6
695	'mixed deciduous evergreen forest ', & ! 7
696	'deciduous forest ', & ! 8
697	'tropical evergreen broadleaved forest', & ! 9
698	'medium/tall grassland/woodland ', & ! 10
699	'desert, sandy ', & ! 11
700	'desert, rocky ', & ! 12
701	'tundra ', & ! 13
702	'land ice ', & ! 14
703	'sea ice ', & ! 15
704	'snow ', & ! 16
705	'bare soil ', & ! 17
706	'asphalt/concrete mix ', & ! 18
707	'asphalt (asphalt concrete) ', & ! 19
708	'concrete (Portland concrete) ', & ! 20
709	'sett ', & ! 21
710	'paving stones ', & ! 22
711	'cobblestone ', & ! 23
712	'metal ', & ! 24
713	'wood ', & ! 25
714	'gravel ', & ! 26
715	'fine gravel ', & ! 27
716	'pebblestone ', & ! 28
717	'woodchips ', & ! 29
718	'tartan (sports) ', & ! 30
719	'artifical turf (sports) ', & ! 31
720	'clay (sports) ', & ! 32
721	'building (dummy) ' & ! 33
722	/)
723
724	INTEGER(iwp) :: albedo_type = 9999999_iwp, & !< Albedo surface type
725	dots_rad = 0_iwp !< starting index for timeseries output
726
727	LOGICAL :: unscheduled_radiation_calls = .TRUE., & !< flag parameter indicating whether additional calls of the radiation code are allowed
728	constant_albedo = .FALSE., & !< flag parameter indicating whether the albedo may change depending on zenith
729	force_radiation_call = .FALSE., & !< flag parameter for unscheduled radiation calls
730	lw_radiation = .TRUE., & !< flag parameter indicating whether longwave radiation shall be calculated
731	radiation = .FALSE., & !< flag parameter indicating whether the radiation model is used
732	sun_up = .TRUE., & !< flag parameter indicating whether the sun is up or down
733	sw_radiation = .TRUE., & !< flag parameter indicating whether shortwave radiation shall be calculated
734	sun_direction = .FALSE., & !< flag parameter indicating whether solar direction shall be calculated
735	average_radiation = .FALSE., & !< flag to set the calculation of radiation averaging for the domain
736	radiation_interactions = .FALSE., & !< flag to activiate RTM (TRUE only if vertical urban/land surface and trees exist)
737	surface_reflections = .TRUE., & !< flag to switch the calculation of radiation interaction between surfaces.
738	!< When it switched off, only the effect of buildings and trees shadow
739	!< will be considered. However fewer SVFs are expected.
740	radiation_interactions_on = .TRUE. !< namelist flag to force RTM activiation regardless to vertical urban/land surface and trees
741
742	REAL(wp) :: albedo = 9999999.9_wp, & !< NAMELIST alpha
743	albedo_lw_dif = 9999999.9_wp, & !< NAMELIST aldif
744	albedo_lw_dir = 9999999.9_wp, & !< NAMELIST aldir
745	albedo_sw_dif = 9999999.9_wp, & !< NAMELIST asdif
746	albedo_sw_dir = 9999999.9_wp, & !< NAMELIST asdir
747	decl_1, & !< declination coef. 1
748	decl_2, & !< declination coef. 2
749	decl_3, & !< declination coef. 3
750	dt_radiation = 0.0_wp, & !< radiation model timestep
751	emissivity = 9999999.9_wp, & !< NAMELIST surface emissivity
752	lon = 0.0_wp, & !< longitude in radians
753	lat = 0.0_wp, & !< latitude in radians
754	net_radiation = 0.0_wp, & !< net radiation at surface
755	skip_time_do_radiation = 0.0_wp, & !< Radiation model is not called before this time
756	sky_trans, & !< sky transmissivity
757	time_radiation = 0.0_wp !< time since last call of radiation code
758
759
760	REAL(wp) :: cos_zenith !< cosine of solar zenith angle, also z-coordinate of solar unit vector
761	REAL(wp) :: sun_dir_lat !< y-coordinate of solar unit vector
762	REAL(wp) :: sun_dir_lon !< x-coordinate of solar unit vector
763
764	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rad_net_av !< average of net radiation (rad_net) at surface
765	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rad_lw_in_xy_av !< average of incoming longwave radiation at surface
766	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rad_lw_out_xy_av !< average of outgoing longwave radiation at surface
767	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rad_sw_in_xy_av !< average of incoming shortwave radiation at surface
768	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rad_sw_out_xy_av !< average of outgoing shortwave radiation at surface
769	!
770	!-- Land surface albedos for solar zenith angle of 60degree after Briegleb (1992)
771	!-- (shortwave, longwave, broadband): sw, lw, bb,
772	REAL(wp), DIMENSION(0:2,1:33), PARAMETER :: albedo_pars = RESHAPE( (/&
773	0.06_wp, 0.06_wp, 0.06_wp, & ! 1
774	0.09_wp, 0.28_wp, 0.19_wp, & ! 2
775	0.11_wp, 0.33_wp, 0.23_wp, & ! 3
776	0.11_wp, 0.33_wp, 0.23_wp, & ! 4
777	0.14_wp, 0.34_wp, 0.25_wp, & ! 5
778	0.06_wp, 0.22_wp, 0.14_wp, & ! 6
779	0.06_wp, 0.27_wp, 0.17_wp, & ! 7
780	0.06_wp, 0.31_wp, 0.19_wp, & ! 8
781	0.06_wp, 0.22_wp, 0.14_wp, & ! 9
782	0.06_wp, 0.28_wp, 0.18_wp, & ! 10
783	0.35_wp, 0.51_wp, 0.43_wp, & ! 11
784	0.24_wp, 0.40_wp, 0.32_wp, & ! 12
785	0.10_wp, 0.27_wp, 0.19_wp, & ! 13
786	0.90_wp, 0.65_wp, 0.77_wp, & ! 14
787	0.90_wp, 0.65_wp, 0.77_wp, & ! 15
788	0.95_wp, 0.70_wp, 0.82_wp, & ! 16
789	0.08_wp, 0.08_wp, 0.08_wp, & ! 17
790	0.17_wp, 0.17_wp, 0.17_wp, & ! 18
791	0.17_wp, 0.17_wp, 0.17_wp, & ! 19
792	0.30_wp, 0.30_wp, 0.30_wp, & ! 20
793	0.17_wp, 0.17_wp, 0.17_wp, & ! 21
794	0.17_wp, 0.17_wp, 0.17_wp, & ! 22
795	0.17_wp, 0.17_wp, 0.17_wp, & ! 23
796	0.17_wp, 0.17_wp, 0.17_wp, & ! 24
797	0.17_wp, 0.17_wp, 0.17_wp, & ! 25
798	0.17_wp, 0.17_wp, 0.17_wp, & ! 26
799	0.17_wp, 0.17_wp, 0.17_wp, & ! 27
800	0.17_wp, 0.17_wp, 0.17_wp, & ! 28
801	0.17_wp, 0.17_wp, 0.17_wp, & ! 29
802	0.17_wp, 0.17_wp, 0.17_wp, & ! 30
803	0.17_wp, 0.17_wp, 0.17_wp, & ! 31
804	0.17_wp, 0.17_wp, 0.17_wp, & ! 32
805	0.17_wp, 0.17_wp, 0.17_wp & ! 33
806	/), (/ 3, 33 /) )
807
808	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: &
809	rad_lw_cs_hr, & !< longwave clear sky radiation heating rate (K/s)
810	rad_lw_cs_hr_av, & !< average of rad_lw_cs_hr
811	rad_lw_hr, & !< longwave radiation heating rate (K/s)
812	rad_lw_hr_av, & !< average of rad_sw_hr
813	rad_lw_in, & !< incoming longwave radiation (W/m2)
814	rad_lw_in_av, & !< average of rad_lw_in
815	rad_lw_out, & !< outgoing longwave radiation (W/m2)
816	rad_lw_out_av, & !< average of rad_lw_out
817	rad_sw_cs_hr, & !< shortwave clear sky radiation heating rate (K/s)
818	rad_sw_cs_hr_av, & !< average of rad_sw_cs_hr
819	rad_sw_hr, & !< shortwave radiation heating rate (K/s)
820	rad_sw_hr_av, & !< average of rad_sw_hr
821	rad_sw_in, & !< incoming shortwave radiation (W/m2)
822	rad_sw_in_av, & !< average of rad_sw_in
823	rad_sw_out, & !< outgoing shortwave radiation (W/m2)
824	rad_sw_out_av !< average of rad_sw_out
825
826
827	!
828	!-- Variables and parameters used in RRTMG only
829	#if defined ( __rrtmg )
830	CHARACTER(LEN=12) :: rrtm_input_file = "RAD_SND_DATA" !< name of the NetCDF input file (sounding data)
831
832
833	!
834	!-- Flag parameters for RRTMGS (should not be changed)
835	INTEGER(iwp), PARAMETER :: rrtm_idrv = 1, & !< flag for longwave upward flux calculation option (0,1)
836	rrtm_inflglw = 2, & !< flag for lw cloud optical properties (0,1,2)
837	rrtm_iceflglw = 0, & !< flag for lw ice particle specifications (0,1,2,3)
838	rrtm_liqflglw = 1, & !< flag for lw liquid droplet specifications
839	rrtm_inflgsw = 2, & !< flag for sw cloud optical properties (0,1,2)
840	rrtm_iceflgsw = 0, & !< flag for sw ice particle specifications (0,1,2,3)
841	rrtm_liqflgsw = 1 !< flag for sw liquid droplet specifications
842
843	!
844	!-- The following variables should be only changed with care, as this will
845	!-- require further setting of some variables, which is currently not
846	!-- implemented (aerosols, ice phase).
847	INTEGER(iwp) :: nzt_rad, & !< upper vertical limit for radiation calculations
848	rrtm_icld = 0, & !< cloud flag (0: clear sky column, 1: cloudy column)
849	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)
850
851	INTEGER(iwp) :: nc_stat !< local variable for storin the result of netCDF calls for error message handling
852
853	LOGICAL :: snd_exists = .FALSE. !< flag parameter to check whether a user-defined input files exists
854	LOGICAL :: sw_exists = .FALSE. !< flag parameter to check whether that required rrtmg sw file exists
855	LOGICAL :: lw_exists = .FALSE. !< flag parameter to check whether that required rrtmg lw file exists
856
857
858	REAL(wp), PARAMETER :: mol_mass_air_d_wv = 1.607793_wp !< molecular weight dry air / water vapor
859
860	REAL(wp), DIMENSION(:), ALLOCATABLE :: hyp_snd, & !< hypostatic pressure from sounding data (hPa)
861	rrtm_tsfc, & !< dummy array for storing surface temperature
862	t_snd !< actual temperature from sounding data (hPa)
863
864	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rrtm_ccl4vmr, & !< CCL4 volume mixing ratio (g/mol)
865	rrtm_cfc11vmr, & !< CFC11 volume mixing ratio (g/mol)
866	rrtm_cfc12vmr, & !< CFC12 volume mixing ratio (g/mol)
867	rrtm_cfc22vmr, & !< CFC22 volume mixing ratio (g/mol)
868	rrtm_ch4vmr, & !< CH4 volume mixing ratio
869	rrtm_cicewp, & !< in-cloud ice water path (g/m2)
870	rrtm_cldfr, & !< cloud fraction (0,1)
871	rrtm_cliqwp, & !< in-cloud liquid water path (g/m2)
872	rrtm_co2vmr, & !< CO2 volume mixing ratio (g/mol)
873	rrtm_emis, & !< surface emissivity (0-1)
874	rrtm_h2ovmr, & !< H2O volume mixing ratio
875	rrtm_n2ovmr, & !< N2O volume mixing ratio
876	rrtm_o2vmr, & !< O2 volume mixing ratio
877	rrtm_o3vmr, & !< O3 volume mixing ratio
878	rrtm_play, & !< pressure layers (hPa, zu-grid)
879	rrtm_plev, & !< pressure layers (hPa, zw-grid)
880	rrtm_reice, & !< cloud ice effective radius (microns)
881	rrtm_reliq, & !< cloud water drop effective radius (microns)
882	rrtm_tlay, & !< actual temperature (K, zu-grid)
883	rrtm_tlev, & !< actual temperature (K, zw-grid)
884	rrtm_lwdflx, & !< RRTM output of incoming longwave radiation flux (W/m2)
885	rrtm_lwdflxc, & !< RRTM output of outgoing clear sky longwave radiation flux (W/m2)
886	rrtm_lwuflx, & !< RRTM output of outgoing longwave radiation flux (W/m2)
887	rrtm_lwuflxc, & !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
888	rrtm_lwuflx_dt, & !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
889	rrtm_lwuflxc_dt,& !< RRTM output of outgoing clear sky longwave radiation flux (W/m2)
890	rrtm_lwhr, & !< RRTM output of longwave radiation heating rate (K/d)
891	rrtm_lwhrc, & !< RRTM output of incoming longwave clear sky radiation heating rate (K/d)
892	rrtm_swdflx, & !< RRTM output of incoming shortwave radiation flux (W/m2)
893	rrtm_swdflxc, & !< RRTM output of outgoing clear sky shortwave radiation flux (W/m2)
894	rrtm_swuflx, & !< RRTM output of outgoing shortwave radiation flux (W/m2)
895	rrtm_swuflxc, & !< RRTM output of incoming clear sky shortwave radiation flux (W/m2)
896	rrtm_swhr, & !< RRTM output of shortwave radiation heating rate (K/d)
897	rrtm_swhrc, & !< RRTM output of incoming shortwave clear sky radiation heating rate (K/d)
898	rrtm_dirdflux, & !< RRTM output of incoming direct shortwave (W/m2)
899	rrtm_difdflux !< RRTM output of incoming diffuse shortwave (W/m2)
900
901	REAL(wp), DIMENSION(1) :: rrtm_aldif, & !< surface albedo for longwave diffuse radiation
902	rrtm_aldir, & !< surface albedo for longwave direct radiation
903	rrtm_asdif, & !< surface albedo for shortwave diffuse radiation
904	rrtm_asdir !< surface albedo for shortwave direct radiation
905
906	!
907	!-- Definition of arrays that are currently not used for calling RRTMG (due to setting of flag parameters)
908	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: rad_lw_cs_in, & !< incoming clear sky longwave radiation (W/m2) (not used)
909	rad_lw_cs_out, & !< outgoing clear sky longwave radiation (W/m2) (not used)
910	rad_sw_cs_in, & !< incoming clear sky shortwave radiation (W/m2) (not used)
911	rad_sw_cs_out, & !< outgoing clear sky shortwave radiation (W/m2) (not used)
912	rrtm_lw_tauaer, & !< lw aerosol optical depth
913	rrtm_lw_taucld, & !< lw in-cloud optical depth
914	rrtm_sw_taucld, & !< sw in-cloud optical depth
915	rrtm_sw_ssacld, & !< sw in-cloud single scattering albedo
916	rrtm_sw_asmcld, & !< sw in-cloud asymmetry parameter
917	rrtm_sw_fsfcld, & !< sw in-cloud forward scattering fraction
918	rrtm_sw_tauaer, & !< sw aerosol optical depth
919	rrtm_sw_ssaaer, & !< sw aerosol single scattering albedo
920	rrtm_sw_asmaer, & !< sw aerosol asymmetry parameter
921	rrtm_sw_ecaer !< sw aerosol optical detph at 0.55 microns (rrtm_iaer = 6 only)
922
923	#endif
924	!
925	!-- Parameters of urban and land surface models
926	INTEGER(iwp) :: nz_urban !< number of layers of urban surface (will be calculated)
927	INTEGER(iwp) :: nz_plant !< number of layers of plant canopy (will be calculated)
928	INTEGER(iwp) :: nz_urban_b !< bottom layer of urban surface (will be calculated)
929	INTEGER(iwp) :: nz_urban_t !< top layer of urban surface (will be calculated)
930	INTEGER(iwp) :: nz_plant_t !< top layer of plant canopy (will be calculated)
931	!-- parameters of urban and land surface models
932	INTEGER(iwp), PARAMETER :: nzut_free = 3 !< number of free layers above top of of topography
933	INTEGER(iwp), PARAMETER :: ndsvf = 2 !< number of dimensions of real values in SVF
934	INTEGER(iwp), PARAMETER :: idsvf = 2 !< number of dimensions of integer values in SVF
935	INTEGER(iwp), PARAMETER :: ndcsf = 1 !< number of dimensions of real values in CSF
936	INTEGER(iwp), PARAMETER :: idcsf = 2 !< number of dimensions of integer values in CSF
937	INTEGER(iwp), PARAMETER :: kdcsf = 4 !< number of dimensions of integer values in CSF calculation array
938	INTEGER(iwp), PARAMETER :: id = 1 !< position of d-index in surfl and surf
939	INTEGER(iwp), PARAMETER :: iz = 2 !< position of k-index in surfl and surf
940	INTEGER(iwp), PARAMETER :: iy = 3 !< position of j-index in surfl and surf
941	INTEGER(iwp), PARAMETER :: ix = 4 !< position of i-index in surfl and surf
942
943	INTEGER(iwp), PARAMETER :: nsurf_type = 10 !< number of surf types incl. phys.(land+urban) & (atm.,sky,boundary) surfaces - 1
944
945	INTEGER(iwp), PARAMETER :: iup_u = 0 !< 0 - index of urban upward surface (ground or roof)
946	INTEGER(iwp), PARAMETER :: idown_u = 1 !< 1 - index of urban downward surface (overhanging)
947	INTEGER(iwp), PARAMETER :: inorth_u = 2 !< 2 - index of urban northward facing wall
948	INTEGER(iwp), PARAMETER :: isouth_u = 3 !< 3 - index of urban southward facing wall
949	INTEGER(iwp), PARAMETER :: ieast_u = 4 !< 4 - index of urban eastward facing wall
950	INTEGER(iwp), PARAMETER :: iwest_u = 5 !< 5 - index of urban westward facing wall
951
952	INTEGER(iwp), PARAMETER :: iup_l = 6 !< 6 - index of land upward surface (ground or roof)
953	INTEGER(iwp), PARAMETER :: inorth_l = 7 !< 7 - index of land northward facing wall
954	INTEGER(iwp), PARAMETER :: isouth_l = 8 !< 8 - index of land southward facing wall
955	INTEGER(iwp), PARAMETER :: ieast_l = 9 !< 9 - index of land eastward facing wall
956	INTEGER(iwp), PARAMETER :: iwest_l = 10 !< 10- index of land westward facing wall
957
958	INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER :: idir = (/0, 0,0, 0,1,-1,0,0, 0,1,-1/) !< surface normal direction x indices
959	INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER :: jdir = (/0, 0,1,-1,0, 0,0,1,-1,0, 0/) !< surface normal direction y indices
960	INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER :: kdir = (/1,-1,0, 0,0, 0,1,0, 0,0, 0/) !< surface normal direction z indices
961	REAL(wp), DIMENSION(0:nsurf_type) :: facearea !< area of single face in respective
962	!< direction (will be calc'd)
963
964
965	!-- indices and sizes of urban and land surface models
966	INTEGER(iwp) :: startland !< start index of block of land and roof surfaces
967	INTEGER(iwp) :: endland !< end index of block of land and roof surfaces
968	INTEGER(iwp) :: nlands !< number of land and roof surfaces in local processor
969	INTEGER(iwp) :: startwall !< start index of block of wall surfaces
970	INTEGER(iwp) :: endwall !< end index of block of wall surfaces
971	INTEGER(iwp) :: nwalls !< number of wall surfaces in local processor
972
973	!-- indices needed for RTM netcdf output subroutines
974	INTEGER(iwp), PARAMETER :: nd = 5
975	CHARACTER(LEN=6), DIMENSION(0:nd-1), PARAMETER :: dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /)
976	INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER :: dirint_u = (/ iup_u, isouth_u, inorth_u, iwest_u, ieast_u /)
977	INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER :: dirint_l = (/ iup_l, isouth_l, inorth_l, iwest_l, ieast_l /)
978	INTEGER(iwp), DIMENSION(0:nd-1) :: dirstart
979	INTEGER(iwp), DIMENSION(0:nd-1) :: dirend
980
981	!-- indices and sizes of urban and land surface models
982	INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET :: surfl_l !< coordinates of i-th local surface in local grid - surfl[:,k] = [d, z, y, x]
983	INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET :: surf_l !< coordinates of i-th surface in grid - surf[:,k] = [d, z, y, x]
984	INTEGER(iwp), DIMENSION(:,:), POINTER :: surfl !< coordinates of i-th local surface in local grid - surfl[:,k] = [d, z, y, x]
985	INTEGER(iwp), DIMENSION(:,:), POINTER :: surf !< coordinates of i-th surface in grid - surf[:,k] = [d, z, y, x]
986	INTEGER(iwp) :: nsurfl !< number of all surfaces in local processor
987	INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET :: nsurfs !< array of number of all surfaces in individual processors
988	INTEGER(iwp) :: nsurf !< global number of surfaces in index array of surfaces (nsurf = proc nsurfs)
989	INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET :: surfstart !< starts of blocks of surfaces for individual processors in array surf (indexed from 1)
990	!< respective block for particular processor is surfstart[iproc+1]+1 : surfstart[iproc+1]+nsurfs[iproc+1]
991
992	!-- block variables needed for calculation of the plant canopy model inside the urban surface model
993	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: pct !< top layer of the plant canopy
994	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: pch !< heights of the plant canopy
995	INTEGER(iwp) :: npcbl = 0 !< number of the plant canopy gridboxes in local processor
996	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: pcbl !< k,j,i coordinates of l-th local plant canopy box pcbl[:,l] = [k, j, i]
997	REAL(wp), DIMENSION(:), ALLOCATABLE :: pcbinsw !< array of absorbed sw radiation for local plant canopy box
998	REAL(wp), DIMENSION(:), ALLOCATABLE :: pcbinswdir !< array of absorbed direct sw radiation for local plant canopy box
999	REAL(wp), DIMENSION(:), ALLOCATABLE :: pcbinswdif !< array of absorbed diffusion sw radiation for local plant canopy box
1000	REAL(wp), DIMENSION(:), ALLOCATABLE :: pcbinlw !< array of absorbed lw radiation for local plant canopy box
1001
1002	!-- configuration parameters (they can be setup in PALM config)
1003	LOGICAL :: raytrace_mpi_rma = .TRUE. !< use MPI RMA to access LAD and gridsurf from remote processes during raytracing
1004	LOGICAL :: rad_angular_discretization = .TRUE.!< whether to use fixed resolution discretization of view factors for
1005	!< reflected radiation (as opposed to all mutually visible pairs)
1006	LOGICAL :: plant_lw_interact = .TRUE. !< whether plant canopy interacts with LW radiation (in addition to SW)
1007	INTEGER(iwp) :: mrt_nlevels = 0 !< number of vertical boxes above surface for which to calculate MRT
1008	LOGICAL :: mrt_skip_roof = .TRUE. !< do not calculate MRT above roof surfaces
1009	LOGICAL :: mrt_include_sw = .TRUE. !< should MRT calculation include SW radiation as well?
1010	LOGICAL :: mrt_geom_human = .TRUE. !< MRT direction weights simulate human body instead of a sphere
1011	INTEGER(iwp) :: nrefsteps = 3 !< number of reflection steps to perform
1012	REAL(wp), PARAMETER :: ext_coef = 0.6_wp !< extinction coefficient (a.k.a. alpha)
1013	INTEGER(iwp), PARAMETER :: rad_version_len = 10 !< length of identification string of rad version
1014	CHARACTER(rad_version_len), PARAMETER :: rad_version = 'RAD v. 3.0' !< identification of version of binary svf and restart files
1015	INTEGER(iwp) :: raytrace_discrete_elevs = 40 !< number of discretization steps for elevation (nadir to zenith)
1016	INTEGER(iwp) :: raytrace_discrete_azims = 80 !< number of discretization steps for azimuth (out of 360 degrees)
1017	REAL(wp) :: max_raytracing_dist = -999.0_wp !< maximum distance for raytracing (in metres)
1018	REAL(wp) :: min_irrf_value = 1e-6_wp !< minimum potential irradiance factor value for raytracing
1019	REAL(wp), DIMENSION(1:30) :: svfnorm_report_thresh = 1e21_wp !< thresholds of SVF normalization values to report
1020	INTEGER(iwp) :: svfnorm_report_num !< number of SVF normalization thresholds to report
1021
1022	!-- radiation related arrays to be used in radiation_interaction routine
1023	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rad_sw_in_dir !< direct sw radiation
1024	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rad_sw_in_diff !< diffusion sw radiation
1025	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rad_lw_in_diff !< diffusion lw radiation
1026
1027	!-- parameters required for RRTMG lower boundary condition
1028	REAL(wp) :: albedo_urb !< albedo value retuned to RRTMG boundary cond.
1029	REAL(wp) :: emissivity_urb !< emissivity value retuned to RRTMG boundary cond.
1030	REAL(wp) :: t_rad_urb !< temperature value retuned to RRTMG boundary cond.
1031
1032	!-- type for calculation of svf
1033	TYPE t_svf
1034	INTEGER(iwp) :: isurflt !<
1035	INTEGER(iwp) :: isurfs !<
1036	REAL(wp) :: rsvf !<
1037	REAL(wp) :: rtransp !<
1038	END TYPE
1039
1040	!-- type for calculation of csf
1041	TYPE t_csf
1042	INTEGER(iwp) :: ip !<
1043	INTEGER(iwp) :: itx !<
1044	INTEGER(iwp) :: ity !<
1045	INTEGER(iwp) :: itz !<
1046	INTEGER(iwp) :: isurfs !< Idx of source face / -1 for sky
1047	REAL(wp) :: rcvf !< Canopy view factor for faces /
1048	!< canopy sink factor for sky (-1)
1049	END TYPE
1050
1051	!-- arrays storing the values of USM
1052	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: svfsurf !< svfsurf[:,isvf] = index of target and source surface for svf[isvf]
1053	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: svf !< array of shape view factors+direct irradiation factors for local surfaces
1054	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfins !< array of sw radiation falling to local surface after i-th reflection
1055	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinl !< array of lw radiation for local surface after i-th reflection
1056
1057	REAL(wp), DIMENSION(:), ALLOCATABLE :: skyvf !< array of sky view factor for each local surface
1058	REAL(wp), DIMENSION(:), ALLOCATABLE :: skyvft !< array of sky view factor including transparency for each local surface
1059	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: dsitrans !< dsidir[isvfl,i] = path transmittance of i-th
1060	!< direction of direct solar irradiance per target surface
1061	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: dsitransc !< dtto per plant canopy box
1062	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: dsidir !< dsidir[:,i] = unit vector of i-th
1063	!< direction of direct solar irradiance
1064	INTEGER(iwp) :: ndsidir !< number of apparent solar directions used
1065	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: dsidir_rev !< dsidir_rev[ielev,iazim] = i for dsidir or -1 if not present
1066
1067	INTEGER(iwp) :: nmrtbl !< No. of local grid boxes for which MRT is calculated
1068	INTEGER(iwp) :: nmrtf !< number of MRT factors for local processor
1069	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: mrtbl !< coordinates of i-th local MRT box - surfl[:,i] = [z, y, x]
1070	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: mrtfsurf !< mrtfsurf[:,imrtf] = index of target MRT box and source surface for mrtf[imrtf]
1071	REAL(wp), DIMENSION(:), ALLOCATABLE :: mrtf !< array of MRT factors for each local MRT box
1072	REAL(wp), DIMENSION(:), ALLOCATABLE :: mrtft !< array of MRT factors including transparency for each local MRT box
1073	REAL(wp), DIMENSION(:), ALLOCATABLE :: mrtsky !< array of sky view factor for each local MRT box
1074	REAL(wp), DIMENSION(:), ALLOCATABLE :: mrtskyt !< array of sky view factor including transparency for each local MRT box
1075	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: mrtdsit !< array of direct solar transparencies for each local MRT box
1076	REAL(wp), DIMENSION(:), ALLOCATABLE :: mrtinsw !< mean SW radiant flux for each MRT box
1077	REAL(wp), DIMENSION(:), ALLOCATABLE :: mrtinlw !< mean LW radiant flux for each MRT box
1078	REAL(wp), DIMENSION(:), ALLOCATABLE :: mrt !< mean radiant temperature for each MRT box
1079	REAL(wp), DIMENSION(:), ALLOCATABLE :: mrtinsw_av !< time average mean SW radiant flux for each MRT box
1080	REAL(wp), DIMENSION(:), ALLOCATABLE :: mrtinlw_av !< time average mean LW radiant flux for each MRT box
1081	REAL(wp), DIMENSION(:), ALLOCATABLE :: mrt_av !< time average mean radiant temperature for each MRT box
1082
1083	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinsw !< array of sw radiation falling to local surface including radiation from reflections
1084	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlw !< array of lw radiation falling to local surface including radiation from reflections
1085	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinswdir !< array of direct sw radiation falling to local surface
1086	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinswdif !< array of diffuse sw radiation from sky and model boundary falling to local surface
1087	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlwdif !< array of diffuse lw radiation from sky and model boundary falling to local surface
1088
1089	!< Outward radiation is only valid for nonvirtual surfaces
1090	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutsl !< array of reflected sw radiation for local surface in i-th reflection
1091	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutll !< array of reflected + emitted lw radiation for local surface in i-th reflection
1092	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfouts !< array of reflected sw radiation for all surfaces in i-th reflection
1093	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutl !< array of reflected + emitted lw radiation for all surfaces in i-th reflection
1094	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlg !< global array of incoming lw radiation from plant canopy
1095	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutsw !< array of total sw radiation outgoing from nonvirtual surfaces surfaces after all reflection
1096	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutlw !< array of total lw radiation outgoing from nonvirtual surfaces surfaces after all reflection
1097	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfemitlwl !< array of emitted lw radiation for local surface used to calculate effective surface temperature for radiation model
1098
1099	!-- block variables needed for calculation of the plant canopy model inside the urban surface model
1100	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: csfsurf !< csfsurf[:,icsf] = index of target surface and csf grid index for csf[icsf]
1101	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: csf !< array of plant canopy sink fators + direct irradiation factors (transparency)
1102	REAL(wp), DIMENSION(:,:,:), POINTER :: sub_lad !< subset of lad_s within urban surface, transformed to plain Z coordinate
1103	REAL(wp), DIMENSION(:), POINTER :: sub_lad_g !< sub_lad globalized (used to avoid MPI RMA calls in raytracing)
1104	REAL(wp) :: prototype_lad !< prototype leaf area density for computing effective optical depth
1105	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nzterr, plantt !< temporary global arrays for raytracing
1106	INTEGER(iwp) :: plantt_max
1107
1108	!-- arrays and variables for calculation of svf and csf
1109	TYPE(t_svf), DIMENSION(:), POINTER :: asvf !< pointer to growing svc array
1110	TYPE(t_csf), DIMENSION(:), POINTER :: acsf !< pointer to growing csf array
1111	TYPE(t_svf), DIMENSION(:), POINTER :: amrtf !< pointer to growing mrtf array
1112	TYPE(t_svf), DIMENSION(:), ALLOCATABLE, TARGET :: asvf1, asvf2 !< realizations of svf array
1113	TYPE(t_csf), DIMENSION(:), ALLOCATABLE, TARGET :: acsf1, acsf2 !< realizations of csf array
1114	TYPE(t_svf), DIMENSION(:), ALLOCATABLE, TARGET :: amrtf1, amrtf2 !< realizations of mftf array
1115	INTEGER(iwp) :: nsvfla !< dimmension of array allocated for storage of svf in local processor
1116	INTEGER(iwp) :: ncsfla !< dimmension of array allocated for storage of csf in local processor
1117	INTEGER(iwp) :: nmrtfa !< dimmension of array allocated for storage of mrt
1118	INTEGER(iwp) :: msvf, mcsf, mmrtf!< mod for swapping the growing array
1119	INTEGER(iwp), PARAMETER :: gasize = 100000_iwp !< initial size of growing arrays
1120	REAL(wp), PARAMETER :: grow_factor = 1.4_wp !< growth factor of growing arrays
1121	INTEGER(iwp) :: nsvfl !< number of svf for local processor
1122	INTEGER(iwp) :: ncsfl !< no. of csf in local processor
1123	!< needed only during calc_svf but must be here because it is
1124	!< shared between subroutines calc_svf and raytrace
1125	INTEGER(iwp), DIMENSION(:,:,:,:), POINTER :: gridsurf !< reverse index of local surfl[d,k,j,i] (for case rad_angular_discretization)
1126	INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE :: gridpcbl !< reverse index of local pcbl[k,j,i]
1127	INTEGER(iwp), PARAMETER :: nsurf_type_u = 6 !< number of urban surf types (used in gridsurf)
1128
1129	!-- temporary arrays for calculation of csf in raytracing
1130	INTEGER(iwp) :: maxboxesg !< max number of boxes ray can cross in the domain
1131	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: boxes !< coordinates of gridboxes being crossed by ray
1132	REAL(wp), DIMENSION(:), ALLOCATABLE :: crlens !< array of crossing lengths of ray for particular grid boxes
1133	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: lad_ip !< array of numbers of process where lad is stored
1134	#if defined( __parallel )
1135	INTEGER(kind=MPI_ADDRESS_KIND), &
1136	DIMENSION(:), ALLOCATABLE :: lad_disp !< array of displaycements of lad in local array of proc lad_ip
1137	INTEGER(iwp) :: win_lad !< MPI RMA window for leaf area density
1138	INTEGER(iwp) :: win_gridsurf !< MPI RMA window for reverse grid surface index
1139	#endif
1140	REAL(wp), DIMENSION(:), ALLOCATABLE :: lad_s_ray !< array of received lad_s for appropriate gridboxes crossed by ray
1141	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: target_surfl
1142	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: rt2_track
1143	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rt2_track_lad
1144	REAL(wp), DIMENSION(:), ALLOCATABLE :: rt2_track_dist
1145	REAL(wp), DIMENSION(:), ALLOCATABLE :: rt2_dist
1146
1147	!-- arrays for time averages
1148	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfradnet_av !< average of net radiation to local surface including radiation from reflections
1149	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinsw_av !< average of sw radiation falling to local surface including radiation from reflections
1150	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlw_av !< average of lw radiation falling to local surface including radiation from reflections
1151	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinswdir_av !< average of direct sw radiation falling to local surface
1152	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinswdif_av !< average of diffuse sw radiation from sky and model boundary falling to local surface
1153	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlwdif_av !< average of diffuse lw radiation from sky and model boundary falling to local surface
1154	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinswref_av !< average of sw radiation falling to surface from reflections
1155	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinlwref_av !< average of lw radiation falling to surface from reflections
1156	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutsw_av !< average of total sw radiation outgoing from nonvirtual surfaces surfaces after all reflection
1157	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfoutlw_av !< average of total lw radiation outgoing from nonvirtual surfaces surfaces after all reflection
1158	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfins_av !< average of array of residua of sw radiation absorbed in surface after last reflection
1159	REAL(wp), DIMENSION(:), ALLOCATABLE :: surfinl_av !< average of array of residua of lw radiation absorbed in surface after last reflection
1160	REAL(wp), DIMENSION(:), ALLOCATABLE :: pcbinlw_av !< Average of pcbinlw
1161	REAL(wp), DIMENSION(:), ALLOCATABLE :: pcbinsw_av !< Average of pcbinsw
1162	REAL(wp), DIMENSION(:), ALLOCATABLE :: pcbinswdir_av !< Average of pcbinswdir
1163	REAL(wp), DIMENSION(:), ALLOCATABLE :: pcbinswdif_av !< Average of pcbinswdif
1164	REAL(wp), DIMENSION(:), ALLOCATABLE :: pcbinswref_av !< Average of pcbinswref
1165
1166
1167	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1168	!-- Energy balance variables
1169	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1170	!-- parameters of the land, roof and wall surfaces
1171	REAL(wp), DIMENSION(:), ALLOCATABLE :: albedo_surf !< albedo of the surface
1172	REAL(wp), DIMENSION(:), ALLOCATABLE :: emiss_surf !< emissivity of the wall surface
1173
1174
1175	INTERFACE radiation_check_data_output
1176	MODULE PROCEDURE radiation_check_data_output
1177	END INTERFACE radiation_check_data_output
1178
1179	INTERFACE radiation_check_data_output_ts
1180	MODULE PROCEDURE radiation_check_data_output_ts
1181	END INTERFACE radiation_check_data_output_ts
1182
1183	INTERFACE radiation_check_data_output_pr
1184	MODULE PROCEDURE radiation_check_data_output_pr
1185	END INTERFACE radiation_check_data_output_pr
1186
1187	INTERFACE radiation_check_parameters
1188	MODULE PROCEDURE radiation_check_parameters
1189	END INTERFACE radiation_check_parameters
1190
1191	INTERFACE radiation_clearsky
1192	MODULE PROCEDURE radiation_clearsky
1193	END INTERFACE radiation_clearsky
1194
1195	INTERFACE radiation_constant
1196	MODULE PROCEDURE radiation_constant
1197	END INTERFACE radiation_constant
1198
1199	INTERFACE radiation_control
1200	MODULE PROCEDURE radiation_control
1201	END INTERFACE radiation_control
1202
1203	INTERFACE radiation_3d_data_averaging
1204	MODULE PROCEDURE radiation_3d_data_averaging
1205	END INTERFACE radiation_3d_data_averaging
1206
1207	INTERFACE radiation_data_output_2d
1208	MODULE PROCEDURE radiation_data_output_2d
1209	END INTERFACE radiation_data_output_2d
1210
1211	INTERFACE radiation_data_output_3d
1212	MODULE PROCEDURE radiation_data_output_3d
1213	END INTERFACE radiation_data_output_3d
1214
1215	INTERFACE radiation_data_output_mask
1216	MODULE PROCEDURE radiation_data_output_mask
1217	END INTERFACE radiation_data_output_mask
1218
1219	INTERFACE radiation_define_netcdf_grid
1220	MODULE PROCEDURE radiation_define_netcdf_grid
1221	END INTERFACE radiation_define_netcdf_grid
1222
1223	INTERFACE radiation_header
1224	MODULE PROCEDURE radiation_header
1225	END INTERFACE radiation_header
1226
1227	INTERFACE radiation_init
1228	MODULE PROCEDURE radiation_init
1229	END INTERFACE radiation_init
1230
1231	INTERFACE radiation_parin
1232	MODULE PROCEDURE radiation_parin
1233	END INTERFACE radiation_parin
1234
1235	INTERFACE radiation_rrtmg
1236	MODULE PROCEDURE radiation_rrtmg
1237	END INTERFACE radiation_rrtmg
1238
1239	#if defined( __rrtmg )
1240	INTERFACE radiation_tendency
1241	MODULE PROCEDURE radiation_tendency
1242	MODULE PROCEDURE radiation_tendency_ij
1243	END INTERFACE radiation_tendency
1244	#endif
1245
1246	INTERFACE radiation_rrd_local
1247	MODULE PROCEDURE radiation_rrd_local
1248	END INTERFACE radiation_rrd_local
1249
1250	INTERFACE radiation_wrd_local
1251	MODULE PROCEDURE radiation_wrd_local
1252	END INTERFACE radiation_wrd_local
1253
1254	INTERFACE radiation_interaction
1255	MODULE PROCEDURE radiation_interaction
1256	END INTERFACE radiation_interaction
1257
1258	INTERFACE radiation_interaction_init
1259	MODULE PROCEDURE radiation_interaction_init
1260	END INTERFACE radiation_interaction_init
1261
1262	INTERFACE radiation_presimulate_solar_pos
1263	MODULE PROCEDURE radiation_presimulate_solar_pos
1264	END INTERFACE radiation_presimulate_solar_pos
1265
1266	INTERFACE radiation_calc_svf
1267	MODULE PROCEDURE radiation_calc_svf
1268	END INTERFACE radiation_calc_svf
1269
1270	INTERFACE radiation_write_svf
1271	MODULE PROCEDURE radiation_write_svf
1272	END INTERFACE radiation_write_svf
1273
1274	INTERFACE radiation_read_svf
1275	MODULE PROCEDURE radiation_read_svf
1276	END INTERFACE radiation_read_svf
1277
1278
1279	SAVE
1280
1281	PRIVATE
1282
1283	!
1284	!-- Public functions / NEEDS SORTING
1285	PUBLIC radiation_check_data_output, radiation_check_data_output_pr, &
1286	radiation_check_data_output_ts, &
1287	radiation_check_parameters, radiation_control, &
1288	radiation_header, radiation_init, radiation_parin, &
1289	radiation_3d_data_averaging, &
1290	radiation_data_output_2d, radiation_data_output_3d, &
1291	radiation_define_netcdf_grid, radiation_wrd_local, &
1292	radiation_rrd_local, radiation_data_output_mask, &
1293	radiation_calc_svf, radiation_write_svf, &
1294	radiation_interaction, radiation_interaction_init, &
1295	radiation_read_svf, radiation_presimulate_solar_pos
1296
1297
1298	!
1299	!-- Public variables and constants / NEEDS SORTING
1300	PUBLIC albedo, albedo_type, decl_1, decl_2, decl_3, dots_rad, dt_radiation,&
1301	emissivity, force_radiation_call, lat, lon, mrt_geom_human, &
1302	mrt_include_sw, mrt_nlevels, mrtbl, mrtinsw, mrtinlw, nmrtbl, &
1303	rad_net_av, radiation, radiation_scheme, rad_lw_in, &
1304	rad_lw_in_av, rad_lw_out, rad_lw_out_av, &
1305	rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av, rad_sw_in, &
1306	rad_sw_in_av, rad_sw_out, rad_sw_out_av, rad_sw_cs_hr, &
1307	rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av, solar_constant, &
1308	skip_time_do_radiation, time_radiation, unscheduled_radiation_calls,&
1309	cos_zenith, calc_zenith, sun_direction, sun_dir_lat, sun_dir_lon, &
1310	idir, jdir, kdir, id, iz, iy, ix, &
1311	iup_u, inorth_u, isouth_u, ieast_u, iwest_u, &
1312	iup_l, inorth_l, isouth_l, ieast_l, iwest_l, &
1313	nsurf_type, nz_urban_b, nz_urban_t, nz_urban, pch, nsurf, &
1314	idsvf, ndsvf, idcsf, ndcsf, kdcsf, pct, &
1315	radiation_interactions, startwall, startland, endland, endwall, &
1316	skyvf, skyvft, radiation_interactions_on, average_radiation, &
1317	rad_sw_in_diff, rad_sw_in_dir
1318
1319
1320	#if defined ( __rrtmg )
1321	PUBLIC radiation_tendency, rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir
1322	#endif
1323
1324	CONTAINS
1325
1326
1327	!------------------------------------------------------------------------------!
1328	! Description:
1329	! ------------
1330	!> This subroutine controls the calls of the radiation schemes
1331	!------------------------------------------------------------------------------!
1332	SUBROUTINE radiation_control
1333
1334
1335	IMPLICIT NONE
1336
1337
1338	SELECT CASE ( TRIM( radiation_scheme ) )
1339
1340	CASE ( 'constant' )
1341	CALL radiation_constant
1342
1343	CASE ( 'clear-sky' )
1344	CALL radiation_clearsky
1345
1346	CASE ( 'rrtmg' )
1347	CALL radiation_rrtmg
1348
1349	CASE DEFAULT
1350
1351	END SELECT
1352
1353
1354	END SUBROUTINE radiation_control
1355
1356	!------------------------------------------------------------------------------!
1357	! Description:
1358	! ------------
1359	!> Check data output for radiation model
1360	!------------------------------------------------------------------------------!
1361	SUBROUTINE radiation_check_data_output( variable, unit, i, ilen, k )
1362
1363
1364	USE control_parameters, &
1365	ONLY: data_output, message_string
1366
1367	IMPLICIT NONE
1368
1369	CHARACTER (LEN=*) :: unit !<
1370	CHARACTER (LEN=*) :: variable !<
1371
1372	INTEGER(iwp) :: i, k
1373	INTEGER(iwp) :: ilen
1374	CHARACTER(LEN=varnamelength) :: var !< TRIM(variable)
1375
1376	var = TRIM(variable)
1377
1378	!-- first process diractional variables
1379	IF ( var(1:12) == 'rtm_rad_net_' .OR. var(1:13) == 'rtm_rad_insw_' .OR. &
1380	var(1:13) == 'rtm_rad_inlw_' .OR. var(1:16) == 'rtm_rad_inswdir_' .OR. &
1381	var(1:16) == 'rtm_rad_inswdif_' .OR. var(1:16) == 'rtm_rad_inswref_' .OR. &
1382	var(1:16) == 'rtm_rad_inlwdif_' .OR. var(1:16) == 'rtm_rad_inlwref_' .OR. &
1383	var(1:14) == 'rtm_rad_outsw_' .OR. var(1:14) == 'rtm_rad_outlw_' .OR. &
1384	var(1:14) == 'rtm_rad_ressw_' .OR. var(1:14) == 'rtm_rad_reslw_' ) THEN
1385	IF ( .NOT. radiation ) THEN
1386	message_string = 'output of "' // TRIM( var ) // '" require'&
1387	// 's radiation = .TRUE.'
1388	CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
1389	ENDIF
1390	unit = 'W/m2'
1391	ELSE IF ( var(1:7) == 'rtm_svf' .OR. var(1:7) == 'rtm_dif' .OR. &
1392	var(1:9) == 'rtm_skyvf' .OR. var(1:9) == 'rtm_skyvft' ) THEN
1393	IF ( .NOT. radiation ) THEN
1394	message_string = 'output of "' // TRIM( var ) // '" require'&
1395	// 's radiation = .TRUE.'
1396	CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
1397	ENDIF
1398	unit = '1'
1399	ELSE
1400	!-- non-directional variables
1401	SELECT CASE ( TRIM( var ) )
1402	CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_lw_in', 'rad_lw_out', &
1403	'rad_sw_cs_hr', 'rad_sw_hr', 'rad_sw_in', 'rad_sw_out' )
1404	IF ( .NOT. radiation .OR. radiation_scheme /= 'rrtmg' ) THEN
1405	message_string = '"output of "' // TRIM( var ) // '" requi' // &
1406	'res radiation = .TRUE. and ' // &
1407	'radiation_scheme = "rrtmg"'
1408	CALL message( 'check_parameters', 'PA0406', 1, 2, 0, 6, 0 )
1409	ENDIF
1410	unit = 'K/h'
1411
1412	CASE ( 'rad_net', 'rrtm_aldif', 'rrtm_aldir', 'rrtm_asdif', &
1413	'rrtm_asdir', 'rad_lw_in', 'rad_lw_out', 'rad_sw_in', &
1414	'rad_sw_out*')
1415	IF ( i == 0 .AND. ilen == 0 .AND. k == 0) THEN
1416	! Workaround for masked output (calls with i=ilen=k=0)
1417	unit = 'illegal'
1418	RETURN
1419	ENDIF
1420	IF ( k == 0 .OR. data_output(i)(ilen-2:ilen) /= '_xy' ) THEN
1421	message_string = 'illegal value for data_output: "' // &
1422	TRIM( var ) // '" & only 2d-horizontal ' // &
1423	'cross sections are allowed for this value'
1424	CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
1425	ENDIF
1426	IF ( .NOT. radiation .OR. radiation_scheme /= "rrtmg" ) THEN
1427	IF ( TRIM( var ) == 'rrtm_aldif*' .OR. &
1428	TRIM( var ) == 'rrtm_aldir*' .OR. &
1429	TRIM( var ) == 'rrtm_asdif*' .OR. &
1430	TRIM( var ) == 'rrtm_asdir*' ) &
1431	THEN
1432	message_string = 'output of "' // TRIM( var ) // '" require'&
1433	// 's radiation = .TRUE. and radiation_sch'&
1434	// 'eme = "rrtmg"'
1435	CALL message( 'check_parameters', 'PA0409', 1, 2, 0, 6, 0 )
1436	ENDIF
1437	ENDIF
1438
1439	IF ( TRIM( var ) == 'rad_net*' ) unit = 'W/m2'
1440	IF ( TRIM( var ) == 'rad_lw_in*' ) unit = 'W/m2'
1441	IF ( TRIM( var ) == 'rad_lw_out*' ) unit = 'W/m2'
1442	IF ( TRIM( var ) == 'rad_sw_in*' ) unit = 'W/m2'
1443	IF ( TRIM( var ) == 'rad_sw_out*' ) unit = 'W/m2'
1444	IF ( TRIM( var ) == 'rad_sw_in' ) unit = 'W/m2'
1445	IF ( TRIM( var ) == 'rrtm_aldif*' ) unit = ''
1446	IF ( TRIM( var ) == 'rrtm_aldir*' ) unit = ''
1447	IF ( TRIM( var ) == 'rrtm_asdif*' ) unit = ''
1448	IF ( TRIM( var ) == 'rrtm_asdir*' ) unit = ''
1449
1450	CASE ( 'rtm_rad_pc_inlw', 'rtm_rad_pc_insw', 'rtm_rad_pc_inswdir', &
1451	'rtm_rad_pc_inswdif', 'rtm_rad_pc_inswref')
1452	IF ( .NOT. radiation ) THEN
1453	message_string = 'output of "' // TRIM( var ) // '" require'&
1454	// 's radiation = .TRUE.'
1455	CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
1456	ENDIF
1457	unit = 'W'
1458
1459	CASE ( 'rtm_mrt', 'rtm_mrt_sw', 'rtm_mrt_lw' )
1460	IF ( i == 0 .AND. ilen == 0 .AND. k == 0) THEN
1461	! Workaround for masked output (calls with i=ilen=k=0)
1462	unit = 'illegal'
1463	RETURN
1464	ENDIF
1465
1466	IF ( .NOT. radiation ) THEN
1467	message_string = 'output of "' // TRIM( var ) // '" require'&
1468	// 's radiation = .TRUE.'
1469	CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
1470	ENDIF
1471	IF ( mrt_nlevels == 0 ) THEN
1472	message_string = 'output of "' // TRIM( var ) // '" require'&
1473	// 's mrt_nlevels > 0'
1474	CALL message( 'check_parameters', 'PA0510', 1, 2, 0, 6, 0 )
1475	ENDIF
1476	IF ( TRIM( var ) == 'rtm_mrt_sw' .AND. .NOT. mrt_include_sw ) THEN
1477	message_string = 'output of "' // TRIM( var ) // '" require'&
1478	// 's rtm_mrt_sw = .TRUE.'
1479	CALL message( 'check_parameters', 'PA0511', 1, 2, 0, 6, 0 )
1480	ENDIF
1481	IF ( TRIM( var ) == 'rtm_mrt' ) THEN
1482	unit = 'K'
1483	ELSE
1484	unit = 'W m-2'
1485	ENDIF
1486
1487	CASE DEFAULT
1488	unit = 'illegal'
1489
1490	END SELECT
1491	ENDIF
1492
1493	END SUBROUTINE radiation_check_data_output
1494
1495
1496	!------------------------------------------------------------------------------!
1497	! Description:
1498	! ------------
1499	!> Set module-specific timeseries units and labels
1500	!------------------------------------------------------------------------------!
1501	SUBROUTINE radiation_check_data_output_ts( dots_max, dots_num )
1502
1503
1504	INTEGER(iwp), INTENT(IN) :: dots_max
1505	INTEGER(iwp), INTENT(INOUT) :: dots_num
1506
1507	!
1508	!-- Next line is just to avoid compiler warning about unused variable.
1509	IF ( dots_max == 0 ) CONTINUE
1510
1511	!
1512	!-- Temporary solution to add LSM and radiation time series to the default
1513	!-- output
1514	IF ( land_surface .OR. radiation ) THEN
1515	IF ( TRIM( radiation_scheme ) == 'rrtmg' ) THEN
1516	dots_num = dots_num + 15
1517	ELSE
1518	dots_num = dots_num + 11
1519	ENDIF
1520	ENDIF
1521
1522
1523	END SUBROUTINE radiation_check_data_output_ts
1524
1525	!------------------------------------------------------------------------------!
1526	! Description:
1527	! ------------
1528	!> Check data output of profiles for radiation model
1529	!------------------------------------------------------------------------------!
1530	SUBROUTINE radiation_check_data_output_pr( variable, var_count, unit, &
1531	dopr_unit )
1532
1533	USE arrays_3d, &
1534	ONLY: zu
1535
1536	USE control_parameters, &
1537	ONLY: data_output_pr, message_string
1538
1539	USE indices
1540
1541	USE profil_parameter
1542
1543	USE statistics
1544
1545	IMPLICIT NONE
1546
1547	CHARACTER (LEN=*) :: unit !<
1548	CHARACTER (LEN=*) :: variable !<
1549	CHARACTER (LEN=*) :: dopr_unit !< local value of dopr_unit
1550
1551	INTEGER(iwp) :: var_count !<
1552
1553	SELECT CASE ( TRIM( variable ) )
1554
1555	CASE ( 'rad_net' )
1556	IF ( ( .NOT. radiation ) .OR. radiation_scheme == 'constant' )&
1557	THEN
1558	message_string = 'data_output_pr = ' // &
1559	TRIM( data_output_pr(var_count) ) // ' is' // &
1560	'not available for radiation = .FALSE. or ' //&
1561	'radiation_scheme = "constant"'
1562	CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
1563	ELSE
1564	dopr_index(var_count) = 99
1565	dopr_unit = 'W/m2'
1566	hom(:,2,99,:) = SPREAD( zw, 2, statistic_regions+1 )
1567	unit = dopr_unit
1568	ENDIF
1569
1570	CASE ( 'rad_lw_in' )
1571	IF ( ( .NOT. radiation) .OR. radiation_scheme == 'constant' ) &
1572	THEN
1573	message_string = 'data_output_pr = ' // &
1574	TRIM( data_output_pr(var_count) ) // ' is' // &
1575	'not available for radiation = .FALSE. or ' //&
1576	'radiation_scheme = "constant"'
1577	CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
1578	ELSE
1579	dopr_index(var_count) = 100
1580	dopr_unit = 'W/m2'
1581	hom(:,2,100,:) = SPREAD( zw, 2, statistic_regions+1 )
1582	unit = dopr_unit
1583	ENDIF
1584
1585	CASE ( 'rad_lw_out' )
1586	IF ( ( .NOT. radiation ) .OR. radiation_scheme == 'constant' ) &
1587	THEN
1588	message_string = 'data_output_pr = ' // &
1589	TRIM( data_output_pr(var_count) ) // ' is' // &
1590	'not available for radiation = .FALSE. or ' //&
1591	'radiation_scheme = "constant"'
1592	CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
1593	ELSE
1594	dopr_index(var_count) = 101
1595	dopr_unit = 'W/m2'
1596	hom(:,2,101,:) = SPREAD( zw, 2, statistic_regions+1 )
1597	unit = dopr_unit
1598	ENDIF
1599
1600	CASE ( 'rad_sw_in' )
1601	IF ( ( .NOT. radiation ) .OR. radiation_scheme == 'constant' ) &
1602	THEN
1603	message_string = 'data_output_pr = ' // &
1604	TRIM( data_output_pr(var_count) ) // ' is' // &
1605	'not available for radiation = .FALSE. or ' //&
1606	'radiation_scheme = "constant"'
1607	CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
1608	ELSE
1609	dopr_index(var_count) = 102
1610	dopr_unit = 'W/m2'
1611	hom(:,2,102,:) = SPREAD( zw, 2, statistic_regions+1 )
1612	unit = dopr_unit
1613	ENDIF
1614
1615	CASE ( 'rad_sw_out')
1616	IF ( ( .NOT. radiation ) .OR. radiation_scheme == 'constant' )&
1617	THEN
1618	message_string = 'data_output_pr = ' // &
1619	TRIM( data_output_pr(var_count) ) // ' is' // &
1620	'not available for radiation = .FALSE. or ' //&
1621	'radiation_scheme = "constant"'
1622	CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
1623	ELSE
1624	dopr_index(var_count) = 103
1625	dopr_unit = 'W/m2'
1626	hom(:,2,103,:) = SPREAD( zw, 2, statistic_regions+1 )
1627	unit = dopr_unit
1628	ENDIF
1629
1630	CASE ( 'rad_lw_cs_hr' )
1631	IF ( ( .NOT. radiation ) .OR. radiation_scheme /= 'rrtmg' ) &
1632	THEN
1633	message_string = 'data_output_pr = ' // &
1634	TRIM( data_output_pr(var_count) ) // ' is' // &
1635	'not available for radiation = .FALSE. or ' //&
1636	'radiation_scheme /= "rrtmg"'
1637	CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
1638	ELSE
1639	dopr_index(var_count) = 104
1640	dopr_unit = 'K/h'
1641	hom(:,2,104,:) = SPREAD( zu, 2, statistic_regions+1 )
1642	unit = dopr_unit
1643	ENDIF
1644
1645	CASE ( 'rad_lw_hr' )
1646	IF ( ( .NOT. radiation ) .OR. radiation_scheme /= 'rrtmg' ) &
1647	THEN
1648	message_string = 'data_output_pr = ' // &
1649	TRIM( data_output_pr(var_count) ) // ' is' // &
1650	'not available for radiation = .FALSE. or ' //&
1651	'radiation_scheme /= "rrtmg"'
1652	CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
1653	ELSE
1654	dopr_index(var_count) = 105
1655	dopr_unit = 'K/h'
1656	hom(:,2,105,:) = SPREAD( zu, 2, statistic_regions+1 )
1657	unit = dopr_unit
1658	ENDIF
1659
1660	CASE ( 'rad_sw_cs_hr' )
1661	IF ( ( .NOT. radiation ) .OR. radiation_scheme /= 'rrtmg' ) &
1662	THEN
1663	message_string = 'data_output_pr = ' // &
1664	TRIM( data_output_pr(var_count) ) // ' is' // &
1665	'not available for radiation = .FALSE. or ' //&
1666	'radiation_scheme /= "rrtmg"'
1667	CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
1668	ELSE
1669	dopr_index(var_count) = 106
1670	dopr_unit = 'K/h'
1671	hom(:,2,106,:) = SPREAD( zu, 2, statistic_regions+1 )
1672	unit = dopr_unit
1673	ENDIF
1674
1675	CASE ( 'rad_sw_hr' )
1676	IF ( ( .NOT. radiation ) .OR. radiation_scheme /= 'rrtmg' ) &
1677	THEN
1678	message_string = 'data_output_pr = ' // &
1679	TRIM( data_output_pr(var_count) ) // ' is' // &
1680	'not available for radiation = .FALSE. or ' //&
1681	'radiation_scheme /= "rrtmg"'
1682	CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
1683	ELSE
1684	dopr_index(var_count) = 107
1685	dopr_unit = 'K/h'
1686	hom(:,2,107,:) = SPREAD( zu, 2, statistic_regions+1 )
1687	unit = dopr_unit
1688	ENDIF
1689
1690
1691	CASE DEFAULT
1692	unit = 'illegal'
1693
1694	END SELECT
1695
1696
1697	END SUBROUTINE radiation_check_data_output_pr
1698
1699
1700	!------------------------------------------------------------------------------!
1701	! Description:
1702	! ------------
1703	!> Check parameters routine for radiation model
1704	!------------------------------------------------------------------------------!
1705	SUBROUTINE radiation_check_parameters
1706
1707	USE control_parameters, &
1708	ONLY: land_surface, message_string, urban_surface
1709
1710	USE netcdf_data_input_mod, &
1711	ONLY: input_pids_static
1712
1713	IMPLICIT NONE
1714
1715	!
1716	!-- In case no urban-surface or land-surface model is applied, usage of
1717	!-- a radiation model make no sense.
1718	IF ( .NOT. land_surface .AND. .NOT. urban_surface ) THEN
1719	message_string = 'Usage of radiation module is only allowed if ' // &
1720	'land-surface and/or urban-surface model is applied.'
1721	CALL message( 'check_parameters', 'PA0486', 1, 2, 0, 6, 0 )
1722	ENDIF
1723
1724	IF ( radiation_scheme /= 'constant' .AND. &
1725	radiation_scheme /= 'clear-sky' .AND. &
1726	radiation_scheme /= 'rrtmg' ) THEN
1727	message_string = 'unknown radiation_scheme = '// &
1728	TRIM( radiation_scheme )
1729	CALL message( 'check_parameters', 'PA0405', 1, 2, 0, 6, 0 )
1730	ELSEIF ( radiation_scheme == 'rrtmg' ) THEN
1731	#if ! defined ( __rrtmg )
1732	message_string = 'radiation_scheme = "rrtmg" requires ' // &
1733	'compilation of PALM with pre-processor ' // &
1734	'directive -D__rrtmg'
1735	CALL message( 'check_parameters', 'PA0407', 1, 2, 0, 6, 0 )
1736	#endif
1737	#if defined ( __rrtmg ) && ! defined( __netcdf )
1738	message_string = 'radiation_scheme = "rrtmg" requires ' // &
1739	'the use of NetCDF (preprocessor directive ' // &
1740	'-D__netcdf'
1741	CALL message( 'check_parameters', 'PA0412', 1, 2, 0, 6, 0 )
1742	#endif
1743
1744	ENDIF
1745	!
1746	!-- Checks performed only if data is given via namelist only.
1747	IF ( .NOT. input_pids_static ) THEN
1748	IF ( albedo_type == 0 .AND. albedo == 9999999.9_wp .AND. &
1749	radiation_scheme == 'clear-sky') THEN
1750	message_string = 'radiation_scheme = "clear-sky" in combination'//&
1751	'with albedo_type = 0 requires setting of'// &
1752	'albedo /= 9999999.9'
1753	CALL message( 'check_parameters', 'PA0410', 1, 2, 0, 6, 0 )
1754	ENDIF
1755
1756	IF ( albedo_type == 0 .AND. radiation_scheme == 'rrtmg' .AND. &
1757	( albedo_lw_dif == 9999999.9_wp .OR. albedo_lw_dir == 9999999.9_wp&
1758	.OR. albedo_sw_dif == 9999999.9_wp .OR. albedo_sw_dir == 9999999.9_wp&
1759	) ) THEN
1760	message_string = 'radiation_scheme = "rrtmg" in combination' // &
1761	'with albedo_type = 0 requires setting of ' // &
1762	'albedo_lw_dif /= 9999999.9' // &
1763	'albedo_lw_dir /= 9999999.9' // &
1764	'albedo_sw_dif /= 9999999.9 and' // &
1765	'albedo_sw_dir /= 9999999.9'
1766	CALL message( 'check_parameters', 'PA0411', 1, 2, 0, 6, 0 )
1767	ENDIF
1768	ENDIF
1769	!
1770	!-- Parallel rad_angular_discretization without raytrace_mpi_rma is not implemented
1771	#if defined( __parallel )
1772	IF ( rad_angular_discretization .AND. .NOT. raytrace_mpi_rma ) THEN
1773	message_string = 'rad_angular_discretization can only be used ' // &
1774	'together with raytrace_mpi_rma or when ' // &
1775	'no parallelization is applied.'
1776	CALL message( 'check_parameters', 'PA0486', 1, 2, 0, 6, 0 )
1777	ENDIF
1778	#endif
1779
1780	IF ( cloud_droplets .AND. radiation_scheme == 'rrtmg' .AND. &
1781	average_radiation ) THEN
1782	message_string = 'average_radiation = .T. with radiation_scheme'// &
1783	'= "rrtmg" in combination cloud_droplets = .T.'// &
1784	'is not implementd'
1785	CALL message( 'check_parameters', 'PA0560', 1, 2, 0, 6, 0 )
1786	ENDIF
1787
1788	!
1789	!-- Incialize svf normalization reporting histogram
1790	svfnorm_report_num = 1
1791	DO WHILE ( svfnorm_report_thresh(svfnorm_report_num) < 1e20_wp &
1792	.AND. svfnorm_report_num <= 30 )
1793	svfnorm_report_num = svfnorm_report_num + 1
1794	ENDDO
1795	svfnorm_report_num = svfnorm_report_num - 1
1796
1797
1798
1799	END SUBROUTINE radiation_check_parameters
1800
1801
1802	!------------------------------------------------------------------------------!
1803	! Description:
1804	! ------------
1805	!> Initialization of the radiation model
1806	!------------------------------------------------------------------------------!
1807	SUBROUTINE radiation_init
1808
1809	IMPLICIT NONE
1810
1811	INTEGER(iwp) :: i !< running index x-direction
1812	INTEGER(iwp) :: ioff !< offset in x between surface element reference grid point in atmosphere and actual surface
1813	INTEGER(iwp) :: j !< running index y-direction
1814	INTEGER(iwp) :: joff !< offset in y between surface element reference grid point in atmosphere and actual surface
1815	INTEGER(iwp) :: l !< running index for orientation of vertical surfaces
1816	INTEGER(iwp) :: m !< running index for surface elements
1817	#if defined( __rrtmg )
1818	INTEGER(iwp) :: ind_type !< running index for subgrid-surface tiles
1819	#endif
1820
1821	!
1822	!-- Activate radiation_interactions according to the existence of vertical surfaces and/or trees.
1823	!-- The namelist parameter radiation_interactions_on can override this behavior.
1824	!-- (This check cannot be performed in check_parameters, because vertical_surfaces_exist is first set in
1825	!-- init_surface_arrays.)
1826	IF ( radiation_interactions_on ) THEN
1827	IF ( vertical_surfaces_exist .OR. plant_canopy ) THEN
1828	radiation_interactions = .TRUE.
1829	average_radiation = .TRUE.
1830	ELSE
1831	radiation_interactions_on = .FALSE. !< reset namelist parameter: no interactions
1832	!< calculations necessary in case of flat surface
1833	ENDIF
1834	ELSEIF ( vertical_surfaces_exist .OR. plant_canopy ) THEN
1835	message_string = 'radiation_interactions_on is set to .FALSE. although ' // &
1836	'vertical surfaces and/or trees exist. The model will run ' // &
1837	'without RTM (no shadows, no radiation reflections)'
1838	CALL message( 'init_3d_model', 'PA0348', 0, 1, 0, 6, 0 )
1839	ENDIF
1840	!
1841	!-- If required, initialize radiation interactions between surfaces
1842	!-- via sky-view factors. This must be done before radiation is initialized.
1843	IF ( radiation_interactions ) CALL radiation_interaction_init
1844
1845	!
1846	!-- Initialize radiation model
1847	CALL location_message( 'initializing radiation model', .FALSE. )
1848
1849	!
1850	!-- Allocate array for storing the surface net radiation
1851	IF ( .NOT. ALLOCATED ( surf_lsm_h%rad_net ) .AND. &
1852	surf_lsm_h%ns > 0 ) THEN
1853	ALLOCATE( surf_lsm_h%rad_net(1:surf_lsm_h%ns) )
1854	surf_lsm_h%rad_net = 0.0_wp
1855	ENDIF
1856	IF ( .NOT. ALLOCATED ( surf_usm_h%rad_net ) .AND. &
1857	surf_usm_h%ns > 0 ) THEN
1858	ALLOCATE( surf_usm_h%rad_net(1:surf_usm_h%ns) )
1859	surf_usm_h%rad_net = 0.0_wp
1860	ENDIF
1861	DO l = 0, 3
1862	IF ( .NOT. ALLOCATED ( surf_lsm_v(l)%rad_net ) .AND. &
1863	surf_lsm_v(l)%ns > 0 ) THEN
1864	ALLOCATE( surf_lsm_v(l)%rad_net(1:surf_lsm_v(l)%ns) )
1865	surf_lsm_v(l)%rad_net = 0.0_wp
1866	ENDIF
1867	IF ( .NOT. ALLOCATED ( surf_usm_v(l)%rad_net ) .AND. &
1868	surf_usm_v(l)%ns > 0 ) THEN
1869	ALLOCATE( surf_usm_v(l)%rad_net(1:surf_usm_v(l)%ns) )
1870	surf_usm_v(l)%rad_net = 0.0_wp
1871	ENDIF
1872	ENDDO
1873
1874
1875	!
1876	!-- Allocate array for storing the surface longwave (out) radiation change
1877	IF ( .NOT. ALLOCATED ( surf_lsm_h%rad_lw_out_change_0 ) .AND. &
1878	surf_lsm_h%ns > 0 ) THEN
1879	ALLOCATE( surf_lsm_h%rad_lw_out_change_0(1:surf_lsm_h%ns) )
1880	surf_lsm_h%rad_lw_out_change_0 = 0.0_wp
1881	ENDIF
1882	IF ( .NOT. ALLOCATED ( surf_usm_h%rad_lw_out_change_0 ) .AND. &
1883	surf_usm_h%ns > 0 ) THEN
1884	ALLOCATE( surf_usm_h%rad_lw_out_change_0(1:surf_usm_h%ns) )
1885	surf_usm_h%rad_lw_out_change_0 = 0.0_wp
1886	ENDIF
1887	DO l = 0, 3
1888	IF ( .NOT. ALLOCATED ( surf_lsm_v(l)%rad_lw_out_change_0 ) .AND. &
1889	surf_lsm_v(l)%ns > 0 ) THEN
1890	ALLOCATE( surf_lsm_v(l)%rad_lw_out_change_0(1:surf_lsm_v(l)%ns) )
1891	surf_lsm_v(l)%rad_lw_out_change_0 = 0.0_wp
1892	ENDIF
1893	IF ( .NOT. ALLOCATED ( surf_usm_v(l)%rad_lw_out_change_0 ) .AND. &
1894	surf_usm_v(l)%ns > 0 ) THEN
1895	ALLOCATE( surf_usm_v(l)%rad_lw_out_change_0(1:surf_usm_v(l)%ns) )
1896	surf_usm_v(l)%rad_lw_out_change_0 = 0.0_wp
1897	ENDIF
1898	ENDDO
1899
1900	!
1901	!-- Allocate surface arrays for incoming/outgoing short/longwave radiation
1902	IF ( .NOT. ALLOCATED ( surf_lsm_h%rad_sw_in ) .AND. &
1903	surf_lsm_h%ns > 0 ) THEN
1904	ALLOCATE( surf_lsm_h%rad_sw_in(1:surf_lsm_h%ns) )
1905	ALLOCATE( surf_lsm_h%rad_sw_out(1:surf_lsm_h%ns) )
1906	ALLOCATE( surf_lsm_h%rad_sw_dir(1:surf_lsm_h%ns) )
1907	ALLOCATE( surf_lsm_h%rad_sw_dif(1:surf_lsm_h%ns) )
1908	ALLOCATE( surf_lsm_h%rad_sw_ref(1:surf_lsm_h%ns) )
1909	ALLOCATE( surf_lsm_h%rad_sw_res(1:surf_lsm_h%ns) )
1910	ALLOCATE( surf_lsm_h%rad_lw_in(1:surf_lsm_h%ns) )
1911	ALLOCATE( surf_lsm_h%rad_lw_out(1:surf_lsm_h%ns) )
1912	ALLOCATE( surf_lsm_h%rad_lw_dif(1:surf_lsm_h%ns) )
1913	ALLOCATE( surf_lsm_h%rad_lw_ref(1:surf_lsm_h%ns) )
1914	ALLOCATE( surf_lsm_h%rad_lw_res(1:surf_lsm_h%ns) )
1915	surf_lsm_h%rad_sw_in = 0.0_wp
1916	surf_lsm_h%rad_sw_out = 0.0_wp
1917	surf_lsm_h%rad_sw_dir = 0.0_wp
1918	surf_lsm_h%rad_sw_dif = 0.0_wp
1919	surf_lsm_h%rad_sw_ref = 0.0_wp
1920	surf_lsm_h%rad_sw_res = 0.0_wp
1921	surf_lsm_h%rad_lw_in = 0.0_wp
1922	surf_lsm_h%rad_lw_out = 0.0_wp
1923	surf_lsm_h%rad_lw_dif = 0.0_wp
1924	surf_lsm_h%rad_lw_ref = 0.0_wp
1925	surf_lsm_h%rad_lw_res = 0.0_wp
1926	ENDIF
1927	IF ( .NOT. ALLOCATED ( surf_usm_h%rad_sw_in ) .AND. &
1928	surf_usm_h%ns > 0 ) THEN
1929	ALLOCATE( surf_usm_h%rad_sw_in(1:surf_usm_h%ns) )
1930	ALLOCATE( surf_usm_h%rad_sw_out(1:surf_usm_h%ns) )
1931	ALLOCATE( surf_usm_h%rad_sw_dir(1:surf_usm_h%ns) )
1932	ALLOCATE( surf_usm_h%rad_sw_dif(1:surf_usm_h%ns) )
1933	ALLOCATE( surf_usm_h%rad_sw_ref(1:surf_usm_h%ns) )
1934	ALLOCATE( surf_usm_h%rad_sw_res(1:surf_usm_h%ns) )
1935	ALLOCATE( surf_usm_h%rad_lw_in(1:surf_usm_h%ns) )
1936	ALLOCATE( surf_usm_h%rad_lw_out(1:surf_usm_h%ns) )
1937	ALLOCATE( surf_usm_h%rad_lw_dif(1:surf_usm_h%ns) )
1938	ALLOCATE( surf_usm_h%rad_lw_ref(1:surf_usm_h%ns) )
1939	ALLOCATE( surf_usm_h%rad_lw_res(1:surf_usm_h%ns) )
1940	surf_usm_h%rad_sw_in = 0.0_wp
1941	surf_usm_h%rad_sw_out = 0.0_wp
1942	surf_usm_h%rad_sw_dir = 0.0_wp
1943	surf_usm_h%rad_sw_dif = 0.0_wp
1944	surf_usm_h%rad_sw_ref = 0.0_wp
1945	surf_usm_h%rad_sw_res = 0.0_wp
1946	surf_usm_h%rad_lw_in = 0.0_wp
1947	surf_usm_h%rad_lw_out = 0.0_wp
1948	surf_usm_h%rad_lw_dif = 0.0_wp
1949	surf_usm_h%rad_lw_ref = 0.0_wp
1950	surf_usm_h%rad_lw_res = 0.0_wp
1951	ENDIF
1952	DO l = 0, 3
1953	IF ( .NOT. ALLOCATED ( surf_lsm_v(l)%rad_sw_in ) .AND. &
1954	surf_lsm_v(l)%ns > 0 ) THEN
1955	ALLOCATE( surf_lsm_v(l)%rad_sw_in(1:surf_lsm_v(l)%ns) )
1956	ALLOCATE( surf_lsm_v(l)%rad_sw_out(1:surf_lsm_v(l)%ns) )
1957	ALLOCATE( surf_lsm_v(l)%rad_sw_dir(1:surf_lsm_v(l)%ns) )
1958	ALLOCATE( surf_lsm_v(l)%rad_sw_dif(1:surf_lsm_v(l)%ns) )
1959	ALLOCATE( surf_lsm_v(l)%rad_sw_ref(1:surf_lsm_v(l)%ns) )
1960	ALLOCATE( surf_lsm_v(l)%rad_sw_res(1:surf_lsm_v(l)%ns) )
1961
1962	ALLOCATE( surf_lsm_v(l)%rad_lw_in(1:surf_lsm_v(l)%ns) )
1963	ALLOCATE( surf_lsm_v(l)%rad_lw_out(1:surf_lsm_v(l)%ns) )
1964	ALLOCATE( surf_lsm_v(l)%rad_lw_dif(1:surf_lsm_v(l)%ns) )
1965	ALLOCATE( surf_lsm_v(l)%rad_lw_ref(1:surf_lsm_v(l)%ns) )
1966	ALLOCATE( surf_lsm_v(l)%rad_lw_res(1:surf_lsm_v(l)%ns) )
1967
1968	surf_lsm_v(l)%rad_sw_in = 0.0_wp
1969	surf_lsm_v(l)%rad_sw_out = 0.0_wp
1970	surf_lsm_v(l)%rad_sw_dir = 0.0_wp
1971	surf_lsm_v(l)%rad_sw_dif = 0.0_wp
1972	surf_lsm_v(l)%rad_sw_ref = 0.0_wp
1973	surf_lsm_v(l)%rad_sw_res = 0.0_wp
1974
1975	surf_lsm_v(l)%rad_lw_in = 0.0_wp
1976	surf_lsm_v(l)%rad_lw_out = 0.0_wp
1977	surf_lsm_v(l)%rad_lw_dif = 0.0_wp
1978	surf_lsm_v(l)%rad_lw_ref = 0.0_wp
1979	surf_lsm_v(l)%rad_lw_res = 0.0_wp
1980	ENDIF
1981	IF ( .NOT. ALLOCATED ( surf_usm_v(l)%rad_sw_in ) .AND. &
1982	surf_usm_v(l)%ns > 0 ) THEN
1983	ALLOCATE( surf_usm_v(l)%rad_sw_in(1:surf_usm_v(l)%ns) )
1984	ALLOCATE( surf_usm_v(l)%rad_sw_out(1:surf_usm_v(l)%ns) )
1985	ALLOCATE( surf_usm_v(l)%rad_sw_dir(1:surf_usm_v(l)%ns) )
1986	ALLOCATE( surf_usm_v(l)%rad_sw_dif(1:surf_usm_v(l)%ns) )
1987	ALLOCATE( surf_usm_v(l)%rad_sw_ref(1:surf_usm_v(l)%ns) )
1988	ALLOCATE( surf_usm_v(l)%rad_sw_res(1:surf_usm_v(l)%ns) )
1989	ALLOCATE( surf_usm_v(l)%rad_lw_in(1:surf_usm_v(l)%ns) )
1990	ALLOCATE( surf_usm_v(l)%rad_lw_out(1:surf_usm_v(l)%ns) )
1991	ALLOCATE( surf_usm_v(l)%rad_lw_dif(1:surf_usm_v(l)%ns) )
1992	ALLOCATE( surf_usm_v(l)%rad_lw_ref(1:surf_usm_v(l)%ns) )
1993	ALLOCATE( surf_usm_v(l)%rad_lw_res(1:surf_usm_v(l)%ns) )
1994	surf_usm_v(l)%rad_sw_in = 0.0_wp
1995	surf_usm_v(l)%rad_sw_out = 0.0_wp
1996	surf_usm_v(l)%rad_sw_dir = 0.0_wp
1997	surf_usm_v(l)%rad_sw_dif = 0.0_wp
1998	surf_usm_v(l)%rad_sw_ref = 0.0_wp
1999	surf_usm_v(l)%rad_sw_res = 0.0_wp
2000	surf_usm_v(l)%rad_lw_in = 0.0_wp
2001	surf_usm_v(l)%rad_lw_out = 0.0_wp
2002	surf_usm_v(l)%rad_lw_dif = 0.0_wp
2003	surf_usm_v(l)%rad_lw_ref = 0.0_wp
2004	surf_usm_v(l)%rad_lw_res = 0.0_wp
2005	ENDIF
2006	ENDDO
2007	!
2008	!-- Fix net radiation in case of radiation_scheme = 'constant'
2009	IF ( radiation_scheme == 'constant' ) THEN
2010	IF ( ALLOCATED( surf_lsm_h%rad_net ) ) &
2011	surf_lsm_h%rad_net = net_radiation
2012	IF ( ALLOCATED( surf_usm_h%rad_net ) ) &
2013	surf_usm_h%rad_net = net_radiation
2014	!
2015	!-- Todo: weight with inclination angle
2016	DO l = 0, 3
2017	IF ( ALLOCATED( surf_lsm_v(l)%rad_net ) ) &
2018	surf_lsm_v(l)%rad_net = net_radiation
2019	IF ( ALLOCATED( surf_usm_v(l)%rad_net ) ) &
2020	surf_usm_v(l)%rad_net = net_radiation
2021	ENDDO
2022	! radiation = .FALSE.
2023	!
2024	!-- Calculate orbital constants
2025	ELSE
2026	decl_1 = SIN(23.45_wp * pi / 180.0_wp)
2027	decl_2 = 2.0_wp * pi / 365.0_wp
2028	decl_3 = decl_2 * 81.0_wp
2029	lat = latitude * pi / 180.0_wp
2030	lon = longitude * pi / 180.0_wp
2031	ENDIF
2032
2033	IF ( radiation_scheme == 'clear-sky' .OR. &
2034	radiation_scheme == 'constant') THEN
2035
2036
2037	!
2038	!-- Allocate arrays for incoming/outgoing short/longwave radiation
2039	IF ( .NOT. ALLOCATED ( rad_sw_in ) ) THEN
2040	ALLOCATE ( rad_sw_in(0:0,nysg:nyng,nxlg:nxrg) )
2041	ENDIF
2042	IF ( .NOT. ALLOCATED ( rad_sw_out ) ) THEN
2043	ALLOCATE ( rad_sw_out(0:0,nysg:nyng,nxlg:nxrg) )
2044	ENDIF
2045
2046	IF ( .NOT. ALLOCATED ( rad_lw_in ) ) THEN
2047	ALLOCATE ( rad_lw_in(0:0,nysg:nyng,nxlg:nxrg) )
2048	ENDIF
2049	IF ( .NOT. ALLOCATED ( rad_lw_out ) ) THEN
2050	ALLOCATE ( rad_lw_out(0:0,nysg:nyng,nxlg:nxrg) )
2051	ENDIF
2052
2053	!
2054	!-- Allocate average arrays for incoming/outgoing short/longwave radiation
2055	IF ( .NOT. ALLOCATED ( rad_sw_in_av ) ) THEN
2056	ALLOCATE ( rad_sw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
2057	ENDIF
2058	IF ( .NOT. ALLOCATED ( rad_sw_out_av ) ) THEN
2059	ALLOCATE ( rad_sw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
2060	ENDIF
2061
2062	IF ( .NOT. ALLOCATED ( rad_lw_in_av ) ) THEN
2063	ALLOCATE ( rad_lw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
2064	ENDIF
2065	IF ( .NOT. ALLOCATED ( rad_lw_out_av ) ) THEN
2066	ALLOCATE ( rad_lw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
2067	ENDIF
2068	!
2069	!-- Allocate arrays for broadband albedo, and level 1 initialization
2070	!-- via namelist paramter, unless not already allocated.
2071	IF ( .NOT. ALLOCATED(surf_lsm_h%albedo) ) THEN
2072	ALLOCATE( surf_lsm_h%albedo(0:2,1:surf_lsm_h%ns) )
2073	surf_lsm_h%albedo = albedo
2074	ENDIF
2075	IF ( .NOT. ALLOCATED(surf_usm_h%albedo) ) THEN
2076	ALLOCATE( surf_usm_h%albedo(0:2,1:surf_usm_h%ns) )
2077	surf_usm_h%albedo = albedo
2078	ENDIF
2079
2080	DO l = 0, 3
2081	IF ( .NOT. ALLOCATED( surf_lsm_v(l)%albedo ) ) THEN
2082	ALLOCATE( surf_lsm_v(l)%albedo(0:2,1:surf_lsm_v(l)%ns) )
2083	surf_lsm_v(l)%albedo = albedo
2084	ENDIF
2085	IF ( .NOT. ALLOCATED( surf_usm_v(l)%albedo ) ) THEN
2086	ALLOCATE( surf_usm_v(l)%albedo(0:2,1:surf_usm_v(l)%ns) )
2087	surf_usm_v(l)%albedo = albedo
2088	ENDIF
2089	ENDDO
2090	!
2091	!-- Level 2 initialization of broadband albedo via given albedo_type.
2092	!-- Only if albedo_type is non-zero. In case of urban surface and
2093	!-- input data is read from ASCII file, albedo_type will be zero, so that
2094	!-- albedo won't be overwritten.
2095	DO m = 1, surf_lsm_h%ns
2096	IF ( surf_lsm_h%albedo_type(ind_veg_wall,m) /= 0 ) &
2097	surf_lsm_h%albedo(ind_veg_wall,m) = &
2098	albedo_pars(2,surf_lsm_h%albedo_type(ind_veg_wall,m))
2099	IF ( surf_lsm_h%albedo_type(ind_pav_green,m) /= 0 ) &
2100	surf_lsm_h%albedo(ind_pav_green,m) = &
2101	albedo_pars(2,surf_lsm_h%albedo_type(ind_pav_green,m))
2102	IF ( surf_lsm_h%albedo_type(ind_wat_win,m) /= 0 ) &
2103	surf_lsm_h%albedo(ind_wat_win,m) = &
2104	albedo_pars(2,surf_lsm_h%albedo_type(ind_wat_win,m))
2105	ENDDO
2106	DO m = 1, surf_usm_h%ns
2107	IF ( surf_usm_h%albedo_type(ind_veg_wall,m) /= 0 ) &
2108	surf_usm_h%albedo(ind_veg_wall,m) = &
2109	albedo_pars(2,surf_usm_h%albedo_type(ind_veg_wall,m))
2110	IF ( surf_usm_h%albedo_type(ind_pav_green,m) /= 0 ) &
2111	surf_usm_h%albedo(ind_pav_green,m) = &
2112	albedo_pars(2,surf_usm_h%albedo_type(ind_pav_green,m))
2113	IF ( surf_usm_h%albedo_type(ind_wat_win,m) /= 0 ) &
2114	surf_usm_h%albedo(ind_wat_win,m) = &
2115	albedo_pars(2,surf_usm_h%albedo_type(ind_wat_win,m))
2116	ENDDO
2117
2118	DO l = 0, 3
2119	DO m = 1, surf_lsm_v(l)%ns
2120	IF ( surf_lsm_v(l)%albedo_type(ind_veg_wall,m) /= 0 ) &
2121	surf_lsm_v(l)%albedo(ind_veg_wall,m) = &
2122	albedo_pars(2,surf_lsm_v(l)%albedo_type(ind_veg_wall,m))
2123	IF ( surf_lsm_v(l)%albedo_type(ind_pav_green,m) /= 0 ) &
2124	surf_lsm_v(l)%albedo(ind_pav_green,m) = &
2125	albedo_pars(2,surf_lsm_v(l)%albedo_type(ind_pav_green,m))
2126	IF ( surf_lsm_v(l)%albedo_type(ind_wat_win,m) /= 0 ) &
2127	surf_lsm_v(l)%albedo(ind_wat_win,m) = &
2128	albedo_pars(2,surf_lsm_v(l)%albedo_type(ind_wat_win,m))
2129	ENDDO
2130	DO m = 1, surf_usm_v(l)%ns
2131	IF ( surf_usm_v(l)%albedo_type(ind_veg_wall,m) /= 0 ) &
2132	surf_usm_v(l)%albedo(ind_veg_wall,m) = &
2133	albedo_pars(2,surf_usm_v(l)%albedo_type(ind_veg_wall,m))
2134	IF ( surf_usm_v(l)%albedo_type(ind_pav_green,m) /= 0 ) &
2135	surf_usm_v(l)%albedo(ind_pav_green,m) = &
2136	albedo_pars(2,surf_usm_v(l)%albedo_type(ind_pav_green,m))
2137	IF ( surf_usm_v(l)%albedo_type(ind_wat_win,m) /= 0 ) &
2138	surf_usm_v(l)%albedo(ind_wat_win,m) = &
2139	albedo_pars(2,surf_usm_v(l)%albedo_type(ind_wat_win,m))
2140	ENDDO
2141	ENDDO
2142
2143	!
2144	!-- Level 3 initialization at grid points where albedo type is zero.
2145	!-- This case, albedo is taken from file. In case of constant radiation
2146	!-- or clear sky, only broadband albedo is given.
2147	IF ( albedo_pars_f%from_file ) THEN
2148	!
2149	!-- Horizontal surfaces
2150	DO m = 1, surf_lsm_h%ns
2151	i = surf_lsm_h%i(m)
2152	j = surf_lsm_h%j(m)
2153	IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill ) THEN
2154	IF ( surf_lsm_h%albedo_type(ind_veg_wall,m) == 0 ) &
2155	surf_lsm_h%albedo(ind_veg_wall,m) = albedo_pars_f%pars_xy(0,j,i)
2156	IF ( surf_lsm_h%albedo_type(ind_pav_green,m) == 0 ) &
2157	surf_lsm_h%albedo(ind_pav_green,m) = albedo_pars_f%pars_xy(0,j,i)
2158	IF ( surf_lsm_h%albedo_type(ind_wat_win,m) == 0 ) &
2159	surf_lsm_h%albedo(ind_wat_win,m) = albedo_pars_f%pars_xy(0,j,i)
2160	ENDIF
2161	ENDDO
2162	DO m = 1, surf_usm_h%ns
2163	i = surf_usm_h%i(m)
2164	j = surf_usm_h%j(m)
2165	IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill ) THEN
2166	IF ( surf_usm_h%albedo_type(ind_veg_wall,m) == 0 ) &
2167	surf_usm_h%albedo(ind_veg_wall,m) = albedo_pars_f%pars_xy(0,j,i)
2168	IF ( surf_usm_h%albedo_type(ind_pav_green,m) == 0 ) &
2169	surf_usm_h%albedo(ind_pav_green,m) = albedo_pars_f%pars_xy(0,j,i)
2170	IF ( surf_usm_h%albedo_type(ind_wat_win,m) == 0 ) &
2171	surf_usm_h%albedo(ind_wat_win,m) = albedo_pars_f%pars_xy(0,j,i)
2172	ENDIF
2173	ENDDO
2174	!
2175	!-- Vertical surfaces
2176	DO l = 0, 3
2177
2178	ioff = surf_lsm_v(l)%ioff
2179	joff = surf_lsm_v(l)%joff
2180	DO m = 1, surf_lsm_v(l)%ns
2181	i = surf_lsm_v(l)%i(m) + ioff
2182	j = surf_lsm_v(l)%j(m) + joff
2183	IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill ) THEN
2184	IF ( surf_lsm_v(l)%albedo_type(ind_veg_wall,m) == 0 ) &
2185	surf_lsm_v(l)%albedo(ind_veg_wall,m) = albedo_pars_f%pars_xy(0,j,i)
2186	IF ( surf_lsm_v(l)%albedo_type(ind_pav_green,m) == 0 ) &
2187	surf_lsm_v(l)%albedo(ind_pav_green,m) = albedo_pars_f%pars_xy(0,j,i)
2188	IF ( surf_lsm_v(l)%albedo_type(ind_wat_win,m) == 0 ) &
2189	surf_lsm_v(l)%albedo(ind_wat_win,m) = albedo_pars_f%pars_xy(0,j,i)
2190	ENDIF
2191	ENDDO
2192
2193	ioff = surf_usm_v(l)%ioff
2194	joff = surf_usm_v(l)%joff
2195	DO m = 1, surf_usm_h%ns
2196	i = surf_usm_h%i(m) + joff
2197	j = surf_usm_h%j(m) + joff
2198	IF ( albedo_pars_f%pars_xy(0,j,i) /= albedo_pars_f%fill ) THEN
2199	IF ( surf_usm_v(l)%albedo_type(ind_veg_wall,m) == 0 ) &
2200	surf_usm_v(l)%albedo(ind_veg_wall,m) = albedo_pars_f%pars_xy(0,j,i)
2201	IF ( surf_usm_v(l)%albedo_type(ind_pav_green,m) == 0 ) &
2202	surf_usm_v(l)%albedo(ind_pav_green,m) = albedo_pars_f%pars_xy(0,j,i)
2203	IF ( surf_usm_v(l)%albedo_type(ind_wat_win,m) == 0 ) &
2204	surf_lsm_v(l)%albedo(ind_wat_win,m) = albedo_pars_f%pars_xy(0,j,i)
2205	ENDIF
2206	ENDDO
2207	ENDDO
2208
2209	ENDIF
2210	!
2211	!-- Initialization actions for RRTMG
2212	ELSEIF ( radiation_scheme == 'rrtmg' ) THEN
2213	#if defined ( __rrtmg )
2214	!
2215	!-- Allocate albedos for short/longwave radiation, horizontal surfaces
2216	!-- for wall/green/window (USM) or vegetation/pavement/water surfaces
2217	!-- (LSM).
2218	ALLOCATE ( surf_lsm_h%aldif(0:2,1:surf_lsm_h%ns) )
2219	ALLOCATE ( surf_lsm_h%aldir(0:2,1:surf_lsm_h%ns) )
2220	ALLOCATE ( surf_lsm_h%asdif(0:2,1:surf_lsm_h%ns) )
2221	ALLOCATE ( surf_lsm_h%asdir(0:2,1:surf_lsm_h%ns) )
2222	ALLOCATE ( surf_lsm_h%rrtm_aldif(0:2,1:surf_lsm_h%ns) )
2223	ALLOCATE ( surf_lsm_h%rrtm_aldir(0:2,1:surf_lsm_h%ns) )
2224	ALLOCATE ( surf_lsm_h%rrtm_asdif(0:2,1:surf_lsm_h%ns) )
2225	ALLOCATE ( surf_lsm_h%rrtm_asdir(0:2,1:surf_lsm_h%ns) )
2226
2227	ALLOCATE ( surf_usm_h%aldif(0:2,1:surf_usm_h%ns) )
2228	ALLOCATE ( surf_usm_h%aldir(0:2,1:surf_usm_h%ns) )
2229	ALLOCATE ( surf_usm_h%asdif(0:2,1:surf_usm_h%ns) )
2230	ALLOCATE ( surf_usm_h%asdir(0:2,1:surf_usm_h%ns) )
2231	ALLOCATE ( surf_usm_h%rrtm_aldif(0:2,1:surf_usm_h%ns) )
2232	ALLOCATE ( surf_usm_h%rrtm_aldir(0:2,1:surf_usm_h%ns) )
2233	ALLOCATE ( surf_usm_h%rrtm_asdif(0:2,1:surf_usm_h%ns) )
2234	ALLOCATE ( surf_usm_h%rrtm_asdir(0:2,1:surf_usm_h%ns) )
2235
2236	!
2237	!-- Allocate broadband albedo (temporary for the current radiation
2238	!-- implementations)
2239	IF ( .NOT. ALLOCATED(surf_lsm_h%albedo) ) &
2240	ALLOCATE( surf_lsm_h%albedo(0:2,1:surf_lsm_h%ns) )
2241	IF ( .NOT. ALLOCATED(surf_usm_h%albedo) ) &
2242	ALLOCATE( surf_usm_h%albedo(0:2,1:surf_usm_h%ns) )
2243
2244	!
2245	!-- Allocate albedos for short/longwave radiation, vertical surfaces
2246	DO l = 0, 3
2247
2248	ALLOCATE ( surf_lsm_v(l)%aldif(0:2,1:surf_lsm_v(l)%ns) )
2249	ALLOCATE ( surf_lsm_v(l)%aldir(0:2,1:surf_lsm_v(l)%ns) )
2250	ALLOCATE ( surf_lsm_v(l)%asdif(0:2,1:surf_lsm_v(l)%ns) )
2251	ALLOCATE ( surf_lsm_v(l)%asdir(0:2,1:surf_lsm_v(l)%ns) )
2252
2253	ALLOCATE ( surf_lsm_v(l)%rrtm_aldif(0:2,1:surf_lsm_v(l)%ns) )
2254	ALLOCATE ( surf_lsm_v(l)%rrtm_aldir(0:2,1:surf_lsm_v(l)%ns) )
2255	ALLOCATE ( surf_lsm_v(l)%rrtm_asdif(0:2,1:surf_lsm_v(l)%ns) )
2256	ALLOCATE ( surf_lsm_v(l)%rrtm_asdir(0:2,1:surf_lsm_v(l)%ns) )
2257
2258	ALLOCATE ( surf_usm_v(l)%aldif(0:2,1:surf_usm_v(l)%ns) )
2259	ALLOCATE ( surf_usm_v(l)%aldir(0:2,1:surf_usm_v(l)%ns) )
2260	ALLOCATE ( surf_usm_v(l)%asdif(0:2,1:surf_usm_v(l)%ns) )
2261	ALLOCATE ( surf_usm_v(l)%asdir(0:2,1:surf_usm_v(l)%ns) )
2262
2263	ALLOCATE ( surf_usm_v(l)%rrtm_aldif(0:2,1:surf_usm_v(l)%ns) )
2264	ALLOCATE ( surf_usm_v(l)%rrtm_aldir(0:2,1:surf_usm_v(l)%ns) )
2265	ALLOCATE ( surf_usm_v(l)%rrtm_asdif(0:2,1:surf_usm_v(l)%ns) )
2266	ALLOCATE ( surf_usm_v(l)%rrtm_asdir(0:2,1:surf_usm_v(l)%ns) )
2267	!
2268	!-- Allocate broadband albedo (temporary for the current radiation
2269	!-- implementations)
2270	IF ( .NOT. ALLOCATED( surf_lsm_v(l)%albedo ) ) &
2271	ALLOCATE( surf_lsm_v(l)%albedo(0:2,1:surf_lsm_v(l)%ns) )
2272	IF ( .NOT. ALLOCATED( surf_usm_v(l)%albedo ) ) &
2273	ALLOCATE( surf_usm_v(l)%albedo(0:2,1:surf_usm_v(l)%ns) )
2274
2275	ENDDO
2276	!
2277	!-- Level 1 initialization of spectral albedos via namelist
2278	!-- paramters. Please note, this case all surface tiles are initialized
2279	!-- the same.
2280	IF ( surf_lsm_h%ns > 0 ) THEN
2281	surf_lsm_h%aldif = albedo_lw_dif
2282	surf_lsm_h%aldir = albedo_lw_dir
2283	surf_lsm_h%asdif = albedo_sw_dif
2284	surf_lsm_h%asdir = albedo_sw_dir
2285	surf_lsm_h%albedo = albedo_sw_dif
2286	ENDIF
2287	IF ( surf_usm_h%ns > 0 ) THEN
2288	IF ( surf_usm_h%albedo_from_ascii ) THEN
2289	surf_usm_h%aldif = surf_usm_h%albedo
2290	surf_usm_h%aldir = surf_usm_h%albedo
2291	surf_usm_h%asdif = surf_usm_h%albedo
2292	surf_usm_h%asdir = surf_usm_h%albedo
2293	ELSE
2294	surf_usm_h%aldif = albedo_lw_dif
2295	surf_usm_h%aldir = albedo_lw_dir
2296	surf_usm_h%asdif = albedo_sw_dif
2297	surf_usm_h%asdir = albedo_sw_dir
2298	surf_usm_h%albedo = albedo_sw_dif
2299	ENDIF
2300	ENDIF
2301
2302	DO l = 0, 3
2303
2304	IF ( surf_lsm_v(l)%ns > 0 ) THEN
2305	surf_lsm_v(l)%aldif = albedo_lw_dif
2306	surf_lsm_v(l)%aldir = albedo_lw_dir
2307	surf_lsm_v(l)%asdif = albedo_sw_dif
2308	surf_lsm_v(l)%asdir = albedo_sw_dir
2309	surf_lsm_v(l)%albedo = albedo_sw_dif
2310	ENDIF
2311
2312	IF ( surf_usm_v(l)%ns > 0 ) THEN
2313	IF ( surf_usm_v(l)%albedo_from_ascii ) THEN
2314	surf_usm_v(l)%aldif = surf_usm_v(l)%albedo
2315	surf_usm_v(l)%aldir = surf_usm_v(l)%albedo
2316	surf_usm_v(l)%asdif = surf_usm_v(l)%albedo
2317	surf_usm_v(l)%asdir = surf_usm_v(l)%albedo
2318	ELSE
2319	surf_usm_v(l)%aldif = albedo_lw_dif
2320	surf_usm_v(l)%aldir = albedo_lw_dir
2321	surf_usm_v(l)%asdif = albedo_sw_dif
2322	surf_usm_v(l)%asdir = albedo_sw_dir
2323	ENDIF
2324	ENDIF
2325	ENDDO
2326
2327	!
2328	!-- Level 2 initialization of spectral albedos via albedo_type.
2329	!-- Please note, for natural- and urban-type surfaces, a tile approach
2330	!-- is applied so that the resulting albedo is calculated via the weighted
2331	!-- average of respective surface fractions.
2332	DO m = 1, surf_lsm_h%ns
2333	!
2334	!-- Spectral albedos for vegetation/pavement/water surfaces
2335	DO ind_type = 0, 2
2336	IF ( surf_lsm_h%albedo_type(ind_type,m) /= 0 ) THEN
2337	surf_lsm_h%aldif(ind_type,m) = &
2338	albedo_pars(0,surf_lsm_h%albedo_type(ind_type,m))
2339	surf_lsm_h%asdif(ind_type,m) = &
2340	albedo_pars(1,surf_lsm_h%albedo_type(ind_type,m))
2341	surf_lsm_h%aldir(ind_type,m) = &
2342	albedo_pars(0,surf_lsm_h%albedo_type(ind_type,m))
2343	surf_lsm_h%asdir(ind_type,m) = &
2344	albedo_pars(1,surf_lsm_h%albedo_type(ind_type,m))
2345	surf_lsm_h%albedo(ind_type,m) = &
2346	albedo_pars(2,surf_lsm_h%albedo_type(ind_type,m))
2347	ENDIF
2348	ENDDO
2349
2350	ENDDO
2351	!
2352	!-- For urban surface only if albedo has not been already initialized
2353	!-- in the urban-surface model via the ASCII file.
2354	IF ( .NOT. surf_usm_h%albedo_from_ascii ) THEN
2355	DO m = 1, surf_usm_h%ns
2356	!
2357	!-- Spectral albedos for wall/green/window surfaces
2358	DO ind_type = 0, 2
2359	IF ( surf_usm_h%albedo_type(ind_type,m) /= 0 ) THEN
2360	surf_usm_h%aldif(ind_type,m) = &
2361	albedo_pars(0,surf_usm_h%albedo_type(ind_type,m))
2362	surf_usm_h%asdif(ind_type,m) = &
2363	albedo_pars(1,surf_usm_h%albedo_type(ind_type,m))
2364	surf_usm_h%aldir(ind_type,m) = &
2365	albedo_pars(0,surf_usm_h%albedo_type(ind_type,m))
2366	surf_usm_h%asdir(ind_type,m) = &
2367	albedo_pars(1,surf_usm_h%albedo_type(ind_type,m))
2368	surf_usm_h%albedo(ind_type,m) = &
2369	albedo_pars(2,surf_usm_h%albedo_type(ind_type,m))
2370	ENDIF
2371	ENDDO
2372
2373	ENDDO
2374	ENDIF
2375
2376	DO l = 0, 3
2377
2378	DO m = 1, surf_lsm_v(l)%ns
2379	!
2380	!-- Spectral albedos for vegetation/pavement/water surfaces
2381	DO ind_type = 0, 2
2382	IF ( surf_lsm_v(l)%albedo_type(ind_type,m) /= 0 ) THEN
2383	surf_lsm_v(l)%aldif(ind_type,m) = &
2384	albedo_pars(0,surf_lsm_v(l)%albedo_type(ind_type,m))
2385	surf_lsm_v(l)%asdif(ind_type,m) = &
2386	albedo_pars(1,surf_lsm_v(l)%albedo_type(ind_type,m))
2387	surf_lsm_v(l)%aldir(ind_type,m) = &
2388	albedo_pars(0,surf_lsm_v(l)%albedo_type(ind_type,m))
2389	surf_lsm_v(l)%asdir(ind_type,m) = &
2390	albedo_pars(1,surf_lsm_v(l)%albedo_type(ind_type,m))
2391	surf_lsm_v(l)%albedo(ind_type,m) = &
2392	albedo_pars(2,surf_lsm_v(l)%albedo_type(ind_type,m))
2393	ENDIF
2394	ENDDO
2395	ENDDO
2396	!
2397	!-- For urban surface only if albedo has not been already initialized
2398	!-- in the urban-surface model via the ASCII file.
2399	IF ( .NOT. surf_usm_v(l)%albedo_from_ascii ) THEN
2400	DO m = 1, surf_usm_v(l)%ns
2401	!
2402	!-- Spectral albedos for wall/green/window surfaces
2403	DO ind_type = 0, 2
2404	IF ( surf_usm_v(l)%albedo_type(ind_type,m) /= 0 ) THEN
2405	surf_usm_v(l)%aldif(ind_type,m) = &
2406	albedo_pars(0,surf_usm_v(l)%albedo_type(ind_type,m))
2407	surf_usm_v(l)%asdif(ind_type,m) = &
2408	albedo_pars(1,surf_usm_v(l)%albedo_type(ind_type,m))
2409	surf_usm_v(l)%aldir(ind_type,m) = &
2410	albedo_pars(0,surf_usm_v(l)%albedo_type(ind_type,m))
2411	surf_usm_v(l)%asdir(ind_type,m) = &
2412	albedo_pars(1,surf_usm_v(l)%albedo_type(ind_type,m))
2413	surf_usm_v(l)%albedo(ind_type,m) = &
2414	albedo_pars(2,surf_usm_v(l)%albedo_type(ind_type,m))
2415	ENDIF
2416	ENDDO
2417
2418	ENDDO
2419	ENDIF
2420	ENDDO
2421	!
2422	!-- Level 3 initialization at grid points where albedo type is zero.
2423	!-- This case, spectral albedos are taken from file if available
2424	IF ( albedo_pars_f%from_file ) THEN
2425	!
2426	!-- Horizontal
2427	DO m = 1, surf_lsm_h%ns
2428	i = surf_lsm_h%i(m)
2429	j = surf_lsm_h%j(m)
2430	!
2431	!-- Spectral albedos for vegetation/pavement/water surfaces
2432	DO ind_type = 0, 2
2433	IF ( surf_lsm_h%albedo_type(ind_type,m) == 0 ) THEN
2434	IF ( albedo_pars_f%pars_xy(1,j,i) /= albedo_pars_f%fill )&
2435	surf_lsm_h%albedo(ind_type,m) = &
2436	albedo_pars_f%pars_xy(1,j,i)
2437	IF ( albedo_pars_f%pars_xy(1,j,i) /= albedo_pars_f%fill )&
2438	surf_lsm_h%aldir(ind_type,m) = &
2439	albedo_pars_f%pars_xy(1,j,i)
2440	IF ( albedo_pars_f%pars_xy(2,j,i) /= albedo_pars_f%fill )&
2441	surf_lsm_h%aldif(ind_type,m) = &
2442	albedo_pars_f%pars_xy(2,j,i)
2443	IF ( albedo_pars_f%pars_xy(3,j,i) /= albedo_pars_f%fill )&
2444	surf_lsm_h%asdir(ind_type,m) = &
2445	albedo_pars_f%pars_xy(3,j,i)
2446	IF ( albedo_pars_f%pars_xy(4,j,i) /= albedo_pars_f%fill )&
2447	surf_lsm_h%asdif(ind_type,m) = &
2448	albedo_pars_f%pars_xy(4,j,i)
2449	ENDIF
2450	ENDDO
2451	ENDDO
2452	!
2453	!-- For urban surface only if albedo has not been already initialized
2454	!-- in the urban-surface model via the ASCII file.
2455	IF ( .NOT. surf_usm_h%albedo_from_ascii ) THEN
2456	DO m = 1, surf_usm_h%ns
2457	i = surf_usm_h%i(m)
2458	j = surf_usm_h%j(m)
2459	!
2460	!-- Spectral albedos for wall/green/window surfaces
2461	DO ind_type = 0, 2
2462	IF ( surf_usm_h%albedo_type(ind_type,m) == 0 ) THEN
2463	IF ( albedo_pars_f%pars_xy(1,j,i) /= albedo_pars_f%fill )&
2464	surf_usm_h%albedo(ind_type,m) = &
2465	albedo_pars_f%pars_xy(1,j,i)
2466	IF ( albedo_pars_f%pars_xy(1,j,i) /= albedo_pars_f%fill )&
2467	surf_usm_h%aldir(ind_type,m) = &
2468	albedo_pars_f%pars_xy(1,j,i)
2469	IF ( albedo_pars_f%pars_xy(2,j,i) /= albedo_pars_f%fill )&
2470	surf_usm_h%aldif(ind_type,m) = &
2471	albedo_pars_f%pars_xy(2,j,i)
2472	IF ( albedo_pars_f%pars_xy(3,j,i) /= albedo_pars_f%fill )&
2473	surf_usm_h%asdir(ind_type,m) = &
2474	albedo_pars_f%pars_xy(3,j,i)
2475	IF ( albedo_pars_f%pars_xy(4,j,i) /= albedo_pars_f%fill )&
2476	surf_usm_h%asdif(ind_type,m) = &
2477	albedo_pars_f%pars_xy(4,j,i)
2478	ENDIF
2479	ENDDO
2480
2481	ENDDO
2482	ENDIF
2483	!
2484	!-- Vertical
2485	DO l = 0, 3
2486	ioff = surf_lsm_v(l)%ioff
2487	joff = surf_lsm_v(l)%joff
2488
2489	DO m = 1, surf_lsm_v(l)%ns
2490	i = surf_lsm_v(l)%i(m)
2491	j = surf_lsm_v(l)%j(m)
2492	!
2493	!-- Spectral albedos for vegetation/pavement/water surfaces
2494	DO ind_type = 0, 2
2495	IF ( surf_lsm_v(l)%albedo_type(ind_type,m) == 0 ) THEN
2496	IF ( albedo_pars_f%pars_xy(1,j+joff,i+ioff) /= &
2497	albedo_pars_f%fill ) &
2498	surf_lsm_v(l)%albedo(ind_type,m) = &
2499	albedo_pars_f%pars_xy(1,j+joff,i+ioff)
2500	IF ( albedo_pars_f%pars_xy(1,j+joff,i+ioff) /= &
2501	albedo_pars_f%fill ) &
2502	surf_lsm_v(l)%aldir(ind_type,m) = &
2503	albedo_pars_f%pars_xy(1,j+joff,i+ioff)
2504	IF ( albedo_pars_f%pars_xy(2,j+joff,i+ioff) /= &
2505	albedo_pars_f%fill ) &
2506	surf_lsm_v(l)%aldif(ind_type,m) = &
2507	albedo_pars_f%pars_xy(2,j+joff,i+ioff)
2508	IF ( albedo_pars_f%pars_xy(3,j+joff,i+ioff) /= &
2509	albedo_pars_f%fill ) &
2510	surf_lsm_v(l)%asdir(ind_type,m) = &
2511	albedo_pars_f%pars_xy(3,j+joff,i+ioff)
2512	IF ( albedo_pars_f%pars_xy(4,j+joff,i+ioff) /= &
2513	albedo_pars_f%fill ) &
2514	surf_lsm_v(l)%asdif(ind_type,m) = &
2515	albedo_pars_f%pars_xy(4,j+joff,i+ioff)
2516	ENDIF
2517	ENDDO
2518	ENDDO
2519	!
2520	!-- For urban surface only if albedo has not been already initialized
2521	!-- in the urban-surface model via the ASCII file.
2522	IF ( .NOT. surf_usm_v(l)%albedo_from_ascii ) THEN
2523	ioff = surf_usm_v(l)%ioff
2524	joff = surf_usm_v(l)%joff
2525
2526	DO m = 1, surf_usm_v(l)%ns
2527	i = surf_usm_v(l)%i(m)
2528	j = surf_usm_v(l)%j(m)
2529	!
2530	!-- Spectral albedos for wall/green/window surfaces
2531	DO ind_type = 0, 2
2532	IF ( surf_usm_v(l)%albedo_type(ind_type,m) == 0 ) THEN
2533	IF ( albedo_pars_f%pars_xy(1,j+joff,i+ioff) /= &
2534	albedo_pars_f%fill ) &
2535	surf_usm_v(l)%albedo(ind_type,m) = &
2536	albedo_pars_f%pars_xy(1,j+joff,i+ioff)
2537	IF ( albedo_pars_f%pars_xy(1,j+joff,i+ioff) /= &
2538	albedo_pars_f%fill ) &
2539	surf_usm_v(l)%aldir(ind_type,m) = &
2540	albedo_pars_f%pars_xy(1,j+joff,i+ioff)
2541	IF ( albedo_pars_f%pars_xy(2,j+joff,i+ioff) /= &
2542	albedo_pars_f%fill ) &
2543	surf_usm_v(l)%aldif(ind_type,m) = &
2544	albedo_pars_f%pars_xy(2,j+joff,i+ioff)
2545	IF ( albedo_pars_f%pars_xy(3,j+joff,i+ioff) /= &
2546	albedo_pars_f%fill ) &
2547	surf_usm_v(l)%asdir(ind_type,m) = &
2548	albedo_pars_f%pars_xy(3,j+joff,i+ioff)
2549	IF ( albedo_pars_f%pars_xy(4,j+joff,i+ioff) /= &
2550	albedo_pars_f%fill ) &
2551	surf_usm_v(l)%asdif(ind_type,m) = &
2552	albedo_pars_f%pars_xy(4,j+joff,i+ioff)
2553	ENDIF
2554	ENDDO
2555
2556	ENDDO
2557	ENDIF
2558	ENDDO
2559
2560	ENDIF
2561
2562	!
2563	!-- Calculate initial values of current (cosine of) the zenith angle and
2564	!-- whether the sun is up
2565	CALL calc_zenith
2566	! readjust date and time to its initial value
2567	CALL init_date_and_time
2568	!
2569	!-- Calculate initial surface albedo for different surfaces
2570	IF ( .NOT. constant_albedo ) THEN
2571	#if defined( __netcdf )
2572	!
2573	!-- Horizontally aligned natural and urban surfaces
2574	CALL calc_albedo( surf_lsm_h )
2575	CALL calc_albedo( surf_usm_h )
2576	!
2577	!-- Vertically aligned natural and urban surfaces
2578	DO l = 0, 3
2579	CALL calc_albedo( surf_lsm_v(l) )
2580	CALL calc_albedo( surf_usm_v(l) )
2581	ENDDO
2582	#endif
2583	ELSE
2584	!
2585	!-- Initialize sun-inclination independent spectral albedos
2586	!-- Horizontal surfaces
2587	IF ( surf_lsm_h%ns > 0 ) THEN
2588	surf_lsm_h%rrtm_aldir = surf_lsm_h%aldir
2589	surf_lsm_h%rrtm_asdir = surf_lsm_h%asdir
2590	surf_lsm_h%rrtm_aldif = surf_lsm_h%aldif
2591	surf_lsm_h%rrtm_asdif = surf_lsm_h%asdif
2592	ENDIF
2593	IF ( surf_usm_h%ns > 0 ) THEN
2594	surf_usm_h%rrtm_aldir = surf_usm_h%aldir
2595	surf_usm_h%rrtm_asdir = surf_usm_h%asdir
2596	surf_usm_h%rrtm_aldif = surf_usm_h%aldif
2597	surf_usm_h%rrtm_asdif = surf_usm_h%asdif
2598	ENDIF
2599	!
2600	!-- Vertical surfaces
2601	DO l = 0, 3
2602	IF ( surf_lsm_v(l)%ns > 0 ) THEN
2603	surf_lsm_v(l)%rrtm_aldir = surf_lsm_v(l)%aldir
2604	surf_lsm_v(l)%rrtm_asdir = surf_lsm_v(l)%asdir
2605	surf_lsm_v(l)%rrtm_aldif = surf_lsm_v(l)%aldif
2606	surf_lsm_v(l)%rrtm_asdif = surf_lsm_v(l)%asdif
2607	ENDIF
2608	IF ( surf_usm_v(l)%ns > 0 ) THEN
2609	surf_usm_v(l)%rrtm_aldir = surf_usm_v(l)%aldir
2610	surf_usm_v(l)%rrtm_asdir = surf_usm_v(l)%asdir
2611	surf_usm_v(l)%rrtm_aldif = surf_usm_v(l)%aldif
2612	surf_usm_v(l)%rrtm_asdif = surf_usm_v(l)%asdif
2613	ENDIF
2614	ENDDO
2615
2616	ENDIF
2617
2618	!
2619	!-- Allocate 3d arrays of radiative fluxes and heating rates
2620	IF ( .NOT. ALLOCATED ( rad_sw_in ) ) THEN
2621	ALLOCATE ( rad_sw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2622	rad_sw_in = 0.0_wp
2623	ENDIF
2624
2625	IF ( .NOT. ALLOCATED ( rad_sw_in_av ) ) THEN
2626	ALLOCATE ( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2627	ENDIF
2628
2629	IF ( .NOT. ALLOCATED ( rad_sw_out ) ) THEN
2630	ALLOCATE ( rad_sw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2631	rad_sw_out = 0.0_wp
2632	ENDIF
2633
2634	IF ( .NOT. ALLOCATED ( rad_sw_out_av ) ) THEN
2635	ALLOCATE ( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2636	ENDIF
2637
2638	IF ( .NOT. ALLOCATED ( rad_sw_hr ) ) THEN
2639	ALLOCATE ( rad_sw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2640	rad_sw_hr = 0.0_wp
2641	ENDIF
2642
2643	IF ( .NOT. ALLOCATED ( rad_sw_hr_av ) ) THEN
2644	ALLOCATE ( rad_sw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2645	rad_sw_hr_av = 0.0_wp
2646	ENDIF
2647
2648	IF ( .NOT. ALLOCATED ( rad_sw_cs_hr ) ) THEN
2649	ALLOCATE ( rad_sw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2650	rad_sw_cs_hr = 0.0_wp
2651	ENDIF
2652
2653	IF ( .NOT. ALLOCATED ( rad_sw_cs_hr_av ) ) THEN
2654	ALLOCATE ( rad_sw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2655	rad_sw_cs_hr_av = 0.0_wp
2656	ENDIF
2657
2658	IF ( .NOT. ALLOCATED ( rad_lw_in ) ) THEN
2659	ALLOCATE ( rad_lw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2660	rad_lw_in = 0.0_wp
2661	ENDIF
2662
2663	IF ( .NOT. ALLOCATED ( rad_lw_in_av ) ) THEN
2664	ALLOCATE ( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2665	ENDIF
2666
2667	IF ( .NOT. ALLOCATED ( rad_lw_out ) ) THEN
2668	ALLOCATE ( rad_lw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2669	rad_lw_out = 0.0_wp
2670	ENDIF
2671
2672	IF ( .NOT. ALLOCATED ( rad_lw_out_av ) ) THEN
2673	ALLOCATE ( rad_lw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2674	ENDIF
2675
2676	IF ( .NOT. ALLOCATED ( rad_lw_hr ) ) THEN
2677	ALLOCATE ( rad_lw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2678	rad_lw_hr = 0.0_wp
2679	ENDIF
2680
2681	IF ( .NOT. ALLOCATED ( rad_lw_hr_av ) ) THEN
2682	ALLOCATE ( rad_lw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2683	rad_lw_hr_av = 0.0_wp
2684	ENDIF
2685
2686	IF ( .NOT. ALLOCATED ( rad_lw_cs_hr ) ) THEN
2687	ALLOCATE ( rad_lw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2688	rad_lw_cs_hr = 0.0_wp
2689	ENDIF
2690
2691	IF ( .NOT. ALLOCATED ( rad_lw_cs_hr_av ) ) THEN
2692	ALLOCATE ( rad_lw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2693	rad_lw_cs_hr_av = 0.0_wp
2694	ENDIF
2695
2696	ALLOCATE ( rad_sw_cs_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2697	ALLOCATE ( rad_sw_cs_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2698	rad_sw_cs_in = 0.0_wp
2699	rad_sw_cs_out = 0.0_wp
2700
2701	ALLOCATE ( rad_lw_cs_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2702	ALLOCATE ( rad_lw_cs_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2703	rad_lw_cs_in = 0.0_wp
2704	rad_lw_cs_out = 0.0_wp
2705
2706	!
2707	!-- Allocate 1-element array for surface temperature
2708	!-- (RRTMG anticipates an array as passed argument).
2709	ALLOCATE ( rrtm_tsfc(1) )
2710	!
2711	!-- Allocate surface emissivity.
2712	!-- Values will be given directly before calling rrtm_lw.
2713	ALLOCATE ( rrtm_emis(0:0,1:nbndlw+1) )
2714
2715	!
2716	!-- Initialize RRTMG, before check if files are existent
2717	INQUIRE( FILE='rrtmg_lw.nc', EXIST=lw_exists )
2718	IF ( .NOT. lw_exists ) THEN
2719	message_string = 'Input file rrtmg_lw.nc' // &
2720	'&for rrtmg missing. ' // &
2721	'&Please provide <jobname>_lsw file in the INPUT directory.'
2722	CALL message( 'radiation_init', 'PA0583', 1, 2, 0, 6, 0 )
2723	ENDIF
2724	INQUIRE( FILE='rrtmg_sw.nc', EXIST=sw_exists )
2725	IF ( .NOT. sw_exists ) THEN
2726	message_string = 'Input file rrtmg_sw.nc' // &
2727	'&for rrtmg missing. ' // &
2728	'&Please provide <jobname>_rsw file in the INPUT directory.'
2729	CALL message( 'radiation_init', 'PA0584', 1, 2, 0, 6, 0 )
2730	ENDIF
2731
2732	IF ( lw_radiation ) CALL rrtmg_lw_ini ( c_p )
2733	IF ( sw_radiation ) CALL rrtmg_sw_ini ( c_p )
2734
2735	!
2736	!-- Set input files for RRTMG
2737	INQUIRE(FILE="RAD_SND_DATA", EXIST=snd_exists)
2738	IF ( .NOT. snd_exists ) THEN
2739	rrtm_input_file = "rrtmg_lw.nc"
2740	ENDIF
2741
2742	!
2743	!-- Read vertical layers for RRTMG from sounding data
2744	!-- The routine provides nzt_rad, hyp_snd(1:nzt_rad),
2745	!-- t_snd(nzt+2:nzt_rad), rrtm_play(1:nzt_rad), rrtm_plev(1_nzt_rad+1),
2746	!-- rrtm_tlay(nzt+2:nzt_rad), rrtm_tlev(nzt+2:nzt_rad+1)
2747	CALL read_sounding_data
2748
2749	!
2750	!-- Read trace gas profiles from file. This routine provides
2751	!-- the rrtm_ arrays (1:nzt_rad+1)
2752	CALL read_trace_gas_data
2753	#endif
2754	ENDIF
2755
2756	!
2757	!-- Perform user actions if required
2758	CALL user_init_radiation
2759
2760	!
2761	!-- Calculate radiative fluxes at model start
2762	SELECT CASE ( TRIM( radiation_scheme ) )
2763
2764	CASE ( 'rrtmg' )
2765	CALL radiation_rrtmg
2766
2767	CASE ( 'clear-sky' )
2768	CALL radiation_clearsky
2769
2770	CASE ( 'constant' )
2771	CALL radiation_constant
2772
2773	CASE DEFAULT
2774
2775	END SELECT
2776
2777	! readjust date and time to its initial value
2778	CALL init_date_and_time
2779
2780	CALL location_message( 'finished', .TRUE. )
2781
2782	!
2783	!-- Find all discretized apparent solar positions for radiation interaction.
2784	IF ( radiation_interactions ) CALL radiation_presimulate_solar_pos
2785
2786	!
2787	!-- If required, read or calculate and write out the SVF
2788	IF ( radiation_interactions .AND. read_svf) THEN
2789	!
2790	!-- Read sky-view factors and further required data from file
2791	CALL location_message( ' Start reading SVF from file', .FALSE. )
2792	CALL radiation_read_svf()
2793	CALL location_message( ' Reading SVF from file has finished', .TRUE. )
2794
2795	ELSEIF ( radiation_interactions .AND. .NOT. read_svf) THEN
2796	!
2797	!-- calculate SFV and CSF
2798	CALL location_message( ' Start calculation of SVF', .FALSE. )
2799	CALL radiation_calc_svf()
2800	CALL location_message( ' Calculation of SVF has finished', .TRUE. )
2801	ENDIF
2802
2803	IF ( radiation_interactions .AND. write_svf) THEN
2804	!
2805	!-- Write svf, csf svfsurf and csfsurf data to file
2806	CALL location_message( ' Start writing SVF in file', .FALSE. )
2807	CALL radiation_write_svf()
2808	CALL location_message( ' Writing SVF in file has finished', .TRUE. )
2809	ENDIF
2810
2811	!
2812	!-- Adjust radiative fluxes. In case of urban and land surfaces, also
2813	!-- call an initial interaction.
2814	IF ( radiation_interactions ) THEN
2815	CALL radiation_interaction
2816	ENDIF
2817
2818	RETURN
2819
2820	END SUBROUTINE radiation_init
2821
2822
2823	!------------------------------------------------------------------------------!
2824	! Description:
2825	! ------------
2826	!> A simple clear sky radiation model
2827	!------------------------------------------------------------------------------!
2828	SUBROUTINE radiation_clearsky
2829
2830
2831	IMPLICIT NONE
2832
2833	INTEGER(iwp) :: l !< running index for surface orientation
2834	REAL(wp) :: pt1 !< potential temperature at first grid level or mean value at urban layer top
2835	REAL(wp) :: pt1_l !< potential temperature at first grid level or mean value at urban layer top at local subdomain
2836	REAL(wp) :: ql1 !< liquid water mixing ratio at first grid level or mean value at urban layer top
2837	REAL(wp) :: ql1_l !< liquid water mixing ratio at first grid level or mean value at urban layer top at local subdomain
2838
2839	TYPE(surf_type), POINTER :: surf !< pointer on respective surface type, used to generalize routine
2840
2841	!
2842	!-- Calculate current zenith angle
2843	CALL calc_zenith
2844
2845	!
2846	!-- Calculate sky transmissivity
2847	sky_trans = 0.6_wp + 0.2_wp * cos_zenith
2848
2849	!
2850	!-- Calculate value of the Exner function at model surface
2851	!
2852	!-- In case averaged radiation is used, calculate mean temperature and
2853	!-- liquid water mixing ratio at the urban-layer top.
2854	IF ( average_radiation ) THEN
2855	pt1 = 0.0_wp
2856	IF ( bulk_cloud_model .OR. cloud_droplets ) ql1 = 0.0_wp
2857
2858	pt1_l = SUM( pt(nz_urban_t,nys:nyn,nxl:nxr) )
2859	IF ( bulk_cloud_model .OR. cloud_droplets ) ql1_l = SUM( ql(nz_urban_t,nys:nyn,nxl:nxr) )
2860
2861	#if defined( __parallel )
2862	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
2863	CALL MPI_ALLREDUCE( pt1_l, pt1, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
2864	IF ( ierr /= 0 ) THEN
2865	WRITE(9,*) 'Error MPI_AllReduce1:', ierr, pt1_l, pt1
2866	FLUSH(9)
2867	ENDIF
2868
2869	IF ( bulk_cloud_model .OR. cloud_droplets ) THEN
2870	CALL MPI_ALLREDUCE( ql1_l, ql1, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
2871	IF ( ierr /= 0 ) THEN
2872	WRITE(9,*) 'Error MPI_AllReduce2:', ierr, ql1_l, ql1
2873	FLUSH(9)
2874	ENDIF
2875	ENDIF
2876	#else
2877	pt1 = pt1_l
2878	IF ( bulk_cloud_model .OR. cloud_droplets ) ql1 = ql1_l
2879	#endif
2880
2881	IF ( bulk_cloud_model .OR. cloud_droplets ) pt1 = pt1 + lv_d_cp / exner(nz_urban_t) * ql1
2882	!
2883	!-- Finally, divide by number of grid points
2884	pt1 = pt1 / REAL( ( nx + 1 ) * ( ny + 1 ), KIND=wp )
2885	ENDIF
2886	!
2887	!-- Call clear-sky calculation for each surface orientation.
2888	!-- First, horizontal surfaces
2889	surf => surf_lsm_h
2890	CALL radiation_clearsky_surf
2891	surf => surf_usm_h
2892	CALL radiation_clearsky_surf
2893	!
2894	!-- Vertical surfaces
2895	DO l = 0, 3
2896	surf => surf_lsm_v(l)
2897	CALL radiation_clearsky_surf
2898	surf => surf_usm_v(l)
2899	CALL radiation_clearsky_surf
2900	ENDDO
2901
2902	CONTAINS
2903
2904	SUBROUTINE radiation_clearsky_surf
2905
2906	IMPLICIT NONE
2907
2908	INTEGER(iwp) :: i !< index x-direction
2909	INTEGER(iwp) :: j !< index y-direction
2910	INTEGER(iwp) :: k !< index z-direction
2911	INTEGER(iwp) :: m !< running index for surface elements
2912
2913	IF ( surf%ns < 1 ) RETURN
2914
2915	!
2916	!-- Calculate radiation fluxes and net radiation (rad_net) assuming
2917	!-- homogeneous urban radiation conditions.
2918	IF ( average_radiation ) THEN
2919
2920	k = nz_urban_t
2921
2922	surf%rad_sw_in = solar_constant * sky_trans * cos_zenith
2923	surf%rad_sw_out = albedo_urb * surf%rad_sw_in
2924
2925	surf%rad_lw_in = 0.8_wp * sigma_sb * (pt1 * exner(k+1))**4
2926
2927	surf%rad_lw_out = emissivity_urb * sigma_sb * (t_rad_urb)**4 &
2928	+ (1.0_wp - emissivity_urb) * surf%rad_lw_in
2929
2930	surf%rad_net = surf%rad_sw_in - surf%rad_sw_out &
2931	+ surf%rad_lw_in - surf%rad_lw_out
2932
2933	surf%rad_lw_out_change_0 = 3.0_wp * emissivity_urb * sigma_sb &
2934	* (t_rad_urb)**3
2935
2936	!
2937	!-- Calculate radiation fluxes and net radiation (rad_net) for each surface
2938	!-- element.
2939	ELSE
2940
2941	DO m = 1, surf%ns
2942	i = surf%i(m)
2943	j = surf%j(m)
2944	k = surf%k(m)
2945
2946	surf%rad_sw_in(m) = solar_constant * sky_trans * cos_zenith
2947
2948	!
2949	!-- Weighted average according to surface fraction.
2950	!-- ATTENTION: when radiation interactions are switched on the
2951	!-- calculated fluxes below are not actually used as they are
2952	!-- overwritten in radiation_interaction.
2953	surf%rad_sw_out(m) = ( surf%frac(ind_veg_wall,m) * &
2954	surf%albedo(ind_veg_wall,m) &
2955	+ surf%frac(ind_pav_green,m) * &
2956	surf%albedo(ind_pav_green,m) &
2957	+ surf%frac(ind_wat_win,m) * &
2958	surf%albedo(ind_wat_win,m) ) &
2959	* surf%rad_sw_in(m)
2960
2961	surf%rad_lw_out(m) = ( surf%frac(ind_veg_wall,m) * &
2962	surf%emissivity(ind_veg_wall,m) &
2963	+ surf%frac(ind_pav_green,m) * &
2964	surf%emissivity(ind_pav_green,m) &
2965	+ surf%frac(ind_wat_win,m) * &
2966	surf%emissivity(ind_wat_win,m) &
2967	) &
2968	* sigma_sb &
2969	* ( surf%pt_surface(m) * exner(nzb) )**4
2970
2971	surf%rad_lw_out_change_0(m) = &
2972	( surf%frac(ind_veg_wall,m) * &
2973	surf%emissivity(ind_veg_wall,m) &
2974	+ surf%frac(ind_pav_green,m) * &
2975	surf%emissivity(ind_pav_green,m) &
2976	+ surf%frac(ind_wat_win,m) * &
2977	surf%emissivity(ind_wat_win,m) &
2978	) * 3.0_wp * sigma_sb &
2979	* ( surf%pt_surface(m) * exner(nzb) )** 3
2980
2981
2982	IF ( bulk_cloud_model .OR. cloud_droplets ) THEN
2983	pt1 = pt(k,j,i) + lv_d_cp / exner(k) * ql(k,j,i)
2984	surf%rad_lw_in(m) = 0.8_wp * sigma_sb * (pt1 * exner(k))**4
2985	ELSE
2986	surf%rad_lw_in(m) = 0.8_wp * sigma_sb * (pt(k,j,i) * exner(k))**4
2987	ENDIF
2988
2989	surf%rad_net(m) = surf%rad_sw_in(m) - surf%rad_sw_out(m) &
2990	+ surf%rad_lw_in(m) - surf%rad_lw_out(m)
2991
2992	ENDDO
2993
2994	ENDIF
2995
2996	!
2997	!-- Fill out values in radiation arrays
2998	DO m = 1, surf%ns
2999	i = surf%i(m)
3000	j = surf%j(m)
3001	rad_sw_in(0,j,i) = surf%rad_sw_in(m)
3002	rad_sw_out(0,j,i) = surf%rad_sw_out(m)
3003	rad_lw_in(0,j,i) = surf%rad_lw_in(m)
3004	rad_lw_out(0,j,i) = surf%rad_lw_out(m)
3005	ENDDO
3006
3007	END SUBROUTINE radiation_clearsky_surf
3008
3009	END SUBROUTINE radiation_clearsky
3010
3011
3012	!------------------------------------------------------------------------------!
3013	! Description:
3014	! ------------
3015	!> This scheme keeps the prescribed net radiation constant during the run
3016	!------------------------------------------------------------------------------!
3017	SUBROUTINE radiation_constant
3018
3019
3020	IMPLICIT NONE
3021
3022	INTEGER(iwp) :: l !< running index for surface orientation
3023
3024	REAL(wp) :: pt1 !< potential temperature at first grid level or mean value at urban layer top
3025	REAL(wp) :: pt1_l !< potential temperature at first grid level or mean value at urban layer top at local subdomain
3026	REAL(wp) :: ql1 !< liquid water mixing ratio at first grid level or mean value at urban layer top
3027	REAL(wp) :: ql1_l !< liquid water mixing ratio at first grid level or mean value at urban layer top at local subdomain
3028
3029	TYPE(surf_type), POINTER :: surf !< pointer on respective surface type, used to generalize routine
3030
3031	!
3032	!-- In case averaged radiation is used, calculate mean temperature and
3033	!-- liquid water mixing ratio at the urban-layer top.
3034	IF ( average_radiation ) THEN
3035	pt1 = 0.0_wp
3036	IF ( bulk_cloud_model .OR. cloud_droplets ) ql1 = 0.0_wp
3037
3038	pt1_l = SUM( pt(nz_urban_t,nys:nyn,nxl:nxr) )
3039	IF ( bulk_cloud_model .OR. cloud_droplets ) ql1_l = SUM( ql(nz_urban_t,nys:nyn,nxl:nxr) )
3040
3041	#if defined( __parallel )
3042	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
3043	CALL MPI_ALLREDUCE( pt1_l, pt1, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
3044	IF ( ierr /= 0 ) THEN
3045	WRITE(9,*) 'Error MPI_AllReduce3:', ierr, pt1_l, pt1
3046	FLUSH(9)
3047	ENDIF
3048	IF ( bulk_cloud_model .OR. cloud_droplets ) THEN
3049	CALL MPI_ALLREDUCE( ql1_l, ql1, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
3050	IF ( ierr /= 0 ) THEN
3051	WRITE(9,*) 'Error MPI_AllReduce4:', ierr, ql1_l, ql1
3052	FLUSH(9)
3053	ENDIF
3054	ENDIF
3055	#else
3056	pt1 = pt1_l
3057	IF ( bulk_cloud_model .OR. cloud_droplets ) ql1 = ql1_l
3058	#endif
3059	IF ( bulk_cloud_model .OR. cloud_droplets ) pt1 = pt1 + lv_d_cp / exner(nz_urban_t+1) * ql1
3060	!
3061	!-- Finally, divide by number of grid points
3062	pt1 = pt1 / REAL( ( nx + 1 ) * ( ny + 1 ), KIND=wp )
3063	ENDIF
3064
3065	!
3066	!-- First, horizontal surfaces
3067	surf => surf_lsm_h
3068	CALL radiation_constant_surf
3069	surf => surf_usm_h
3070	CALL radiation_constant_surf
3071	!
3072	!-- Vertical surfaces
3073	DO l = 0, 3
3074	surf => surf_lsm_v(l)
3075	CALL radiation_constant_surf
3076	surf => surf_usm_v(l)
3077	CALL radiation_constant_surf
3078	ENDDO
3079
3080	CONTAINS
3081
3082	SUBROUTINE radiation_constant_surf
3083
3084	IMPLICIT NONE
3085
3086	INTEGER(iwp) :: i !< index x-direction
3087	INTEGER(iwp) :: ioff !< offset between surface element and adjacent grid point along x
3088	INTEGER(iwp) :: j !< index y-direction
3089	INTEGER(iwp) :: joff !< offset between surface element and adjacent grid point along y
3090	INTEGER(iwp) :: k !< index z-direction
3091	INTEGER(iwp) :: koff !< offset between surface element and adjacent grid point along z
3092	INTEGER(iwp) :: m !< running index for surface elements
3093
3094	IF ( surf%ns < 1 ) RETURN
3095
3096	!-- Calculate homogenoeus urban radiation fluxes
3097	IF ( average_radiation ) THEN
3098
3099	surf%rad_net = net_radiation
3100
3101	surf%rad_lw_in = 0.8_wp * sigma_sb * (pt1 * exner(nz_urban_t+1))**4
3102
3103	surf%rad_lw_out = emissivity_urb * sigma_sb * (t_rad_urb)**4 &
3104	+ ( 1.0_wp - emissivity_urb ) & ! shouldn't be this a bulk value -- emissivity_urb?
3105	* surf%rad_lw_in
3106
3107	surf%rad_lw_out_change_0 = 3.0_wp * emissivity_urb * sigma_sb &
3108	* t_rad_urb**3
3109
3110	surf%rad_sw_in = ( surf%rad_net - surf%rad_lw_in &
3111	+ surf%rad_lw_out ) &
3112	/ ( 1.0_wp - albedo_urb )
3113
3114	surf%rad_sw_out = albedo_urb * surf%rad_sw_in
3115
3116	!
3117	!-- Calculate radiation fluxes for each surface element
3118	ELSE
3119	!
3120	!-- Determine index offset between surface element and adjacent
3121	!-- atmospheric grid point
3122	ioff = surf%ioff
3123	joff = surf%joff
3124	koff = surf%koff
3125
3126	!
3127	!-- Prescribe net radiation and estimate the remaining radiative fluxes
3128	DO m = 1, surf%ns
3129	i = surf%i(m)
3130	j = surf%j(m)
3131	k = surf%k(m)
3132
3133	surf%rad_net(m) = net_radiation
3134
3135	IF ( bulk_cloud_model .OR. cloud_droplets ) THEN
3136	pt1 = pt(k,j,i) + lv_d_cp / exner(k) * ql(k,j,i)
3137	surf%rad_lw_in(m) = 0.8_wp * sigma_sb * (pt1 * exner(k))**4
3138	ELSE
3139	surf%rad_lw_in(m) = 0.8_wp * sigma_sb * &
3140	( pt(k,j,i) * exner(k) )**4
3141	ENDIF
3142
3143	!
3144	!-- Weighted average according to surface fraction.
3145	surf%rad_lw_out(m) = ( surf%frac(ind_veg_wall,m) * &
3146	surf%emissivity(ind_veg_wall,m) &
3147	+ surf%frac(ind_pav_green,m) * &
3148	surf%emissivity(ind_pav_green,m) &
3149	+ surf%frac(ind_wat_win,m) * &
3150	surf%emissivity(ind_wat_win,m) &
3151	) &
3152	* sigma_sb &
3153	* ( surf%pt_surface(m) * exner(nzb) )**4
3154
3155	surf%rad_sw_in(m) = ( surf%rad_net(m) - surf%rad_lw_in(m) &
3156	+ surf%rad_lw_out(m) ) &
3157	/ ( 1.0_wp - &
3158	( surf%frac(ind_veg_wall,m) * &
3159	surf%albedo(ind_veg_wall,m) &
3160	+ surf%frac(ind_pav_green,m) * &
3161	surf%albedo(ind_pav_green,m) &
3162	+ surf%frac(ind_wat_win,m) * &
3163	surf%albedo(ind_wat_win,m) ) &
3164	)
3165
3166	surf%rad_sw_out(m) = ( surf%frac(ind_veg_wall,m) * &
3167	surf%albedo(ind_veg_wall,m) &
3168	+ surf%frac(ind_pav_green,m) * &
3169	surf%albedo(ind_pav_green,m) &
3170	+ surf%frac(ind_wat_win,m) * &
3171	surf%albedo(ind_wat_win,m) ) &
3172	* surf%rad_sw_in(m)
3173
3174	ENDDO
3175
3176	ENDIF
3177
3178	!
3179	!-- Fill out values in radiation arrays
3180	DO m = 1, surf%ns
3181	i = surf%i(m)
3182	j = surf%j(m)
3183	rad_sw_in(0,j,i) = surf%rad_sw_in(m)
3184	rad_sw_out(0,j,i) = surf%rad_sw_out(m)
3185	rad_lw_in(0,j,i) = surf%rad_lw_in(m)
3186	rad_lw_out(0,j,i) = surf%rad_lw_out(m)
3187	ENDDO
3188
3189	END SUBROUTINE radiation_constant_surf
3190
3191
3192	END SUBROUTINE radiation_constant
3193
3194	!------------------------------------------------------------------------------!
3195	! Description:
3196	! ------------
3197	!> Header output for radiation model
3198	!------------------------------------------------------------------------------!
3199	SUBROUTINE radiation_header ( io )
3200
3201
3202	IMPLICIT NONE
3203
3204	INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file
3205
3206
3207
3208	!
3209	!-- Write radiation model header
3210	WRITE( io, 3 )
3211
3212	IF ( radiation_scheme == "constant" ) THEN
3213	WRITE( io, 4 ) net_radiation
3214	ELSEIF ( radiation_scheme == "clear-sky" ) THEN
3215	WRITE( io, 5 )
3216	ELSEIF ( radiation_scheme == "rrtmg" ) THEN
3217	WRITE( io, 6 )
3218	IF ( .NOT. lw_radiation ) WRITE( io, 10 )
3219	IF ( .NOT. sw_radiation ) WRITE( io, 11 )
3220	ENDIF
3221
3222	IF ( albedo_type_f%from_file .OR. vegetation_type_f%from_file .OR. &
3223	pavement_type_f%from_file .OR. water_type_f%from_file .OR. &
3224	building_type_f%from_file ) THEN
3225	WRITE( io, 13 )
3226	ELSE
3227	IF ( albedo_type == 0 ) THEN
3228	WRITE( io, 7 ) albedo
3229	ELSE
3230	WRITE( io, 8 ) TRIM( albedo_type_name(albedo_type) )
3231	ENDIF
3232	ENDIF
3233	IF ( constant_albedo ) THEN
3234	WRITE( io, 9 )
3235	ENDIF
3236
3237	WRITE( io, 12 ) dt_radiation
3238
3239
3240	3 FORMAT (//' Radiation model information:'/ &
3241	' ----------------------------'/)
3242	4 FORMAT (' --> Using constant net radiation: net_radiation = ', F6.2, &
3243	// 'W/m**2')
3244	5 FORMAT (' --> Simple radiation scheme for clear sky is used (no clouds,',&
3245	' default)')
3246	6 FORMAT (' --> RRTMG scheme is used')
3247	7 FORMAT (/' User-specific surface albedo: albedo =', F6.3)
3248	8 FORMAT (/' Albedo is set for land surface type: ', A)
3249	9 FORMAT (/' --> Albedo is fixed during the run')
3250	10 FORMAT (/' --> Longwave radiation is disabled')
3251	11 FORMAT (/' --> Shortwave radiation is disabled.')
3252	12 FORMAT (' Timestep: dt_radiation = ', F6.2, ' s')
3253	13 FORMAT (/' Albedo is set individually for each xy-location, according ', &
3254	'to given surface type.')
3255
3256
3257	END SUBROUTINE radiation_header
3258
3259
3260	!------------------------------------------------------------------------------!
3261	! Description:
3262	! ------------
3263	!> Parin for &radiation_parameters for radiation model
3264	!------------------------------------------------------------------------------!
3265	SUBROUTINE radiation_parin
3266
3267
3268	IMPLICIT NONE
3269
3270	CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file
3271
3272	NAMELIST /radiation_par/ albedo, albedo_lw_dif, albedo_lw_dir, &
3273	albedo_sw_dif, albedo_sw_dir, albedo_type, &
3274	constant_albedo, dt_radiation, emissivity, &
3275	lw_radiation, max_raytracing_dist, &
3276	min_irrf_value, mrt_geom_human, &
3277	mrt_include_sw, mrt_nlevels, &
3278	mrt_skip_roof, net_radiation, nrefsteps, &
3279	plant_lw_interact, rad_angular_discretization,&
3280	radiation_interactions_on, radiation_scheme, &
3281	raytrace_discrete_azims, &
3282	raytrace_discrete_elevs, raytrace_mpi_rma, &
3283	skip_time_do_radiation, surface_reflections, &
3284	svfnorm_report_thresh, sw_radiation, &
3285	unscheduled_radiation_calls
3286
3287
3288	NAMELIST /radiation_parameters/ albedo, albedo_lw_dif, albedo_lw_dir, &
3289	albedo_sw_dif, albedo_sw_dir, albedo_type, &
3290	constant_albedo, dt_radiation, emissivity, &
3291	lw_radiation, max_raytracing_dist, &
3292	min_irrf_value, mrt_geom_human, &
3293	mrt_include_sw, mrt_nlevels, &
3294	mrt_skip_roof, net_radiation, nrefsteps, &
3295	plant_lw_interact, rad_angular_discretization,&
3296	radiation_interactions_on, radiation_scheme, &
3297	raytrace_discrete_azims, &
3298	raytrace_discrete_elevs, raytrace_mpi_rma, &
3299	skip_time_do_radiation, surface_reflections, &
3300	svfnorm_report_thresh, sw_radiation, &
3301	unscheduled_radiation_calls
3302
3303	line = ' '
3304
3305	!
3306	!-- Try to find radiation model namelist
3307	REWIND ( 11 )
3308	line = ' '
3309	DO WHILE ( INDEX( line, '&radiation_parameters' ) == 0 )
3310	READ ( 11, '(A)', END=12 ) line
3311	ENDDO
3312	BACKSPACE ( 11 )
3313
3314	!
3315	!-- Read user-defined namelist
3316	READ ( 11, radiation_parameters, ERR = 10 )
3317
3318	!
3319	!-- Set flag that indicates that the radiation model is switched on
3320	radiation = .TRUE.
3321
3322	GOTO 14
3323
3324	10 BACKSPACE( 11 )
3325	READ( 11 , '(A)') line
3326	CALL parin_fail_message( 'radiation_parameters', line )
3327	!
3328	!-- Try to find old namelist
3329	12 REWIND ( 11 )
3330	line = ' '
3331	DO WHILE ( INDEX( line, '&radiation_par' ) == 0 )
3332	READ ( 11, '(A)', END=14 ) line
3333	ENDDO
3334	BACKSPACE ( 11 )
3335
3336	!
3337	!-- Read user-defined namelist
3338	READ ( 11, radiation_par, ERR = 13, END = 14 )
3339
3340	message_string = 'namelist radiation_par is deprecated and will be ' // &
3341	'removed in near future. Please use namelist ' // &
3342	'radiation_parameters instead'
3343	CALL message( 'radiation_parin', 'PA0487', 0, 1, 0, 6, 0 )
3344
3345	!
3346	!-- Set flag that indicates that the radiation model is switched on
3347	radiation = .TRUE.
3348
3349	IF ( .NOT. radiation_interactions_on .AND. surface_reflections ) THEN
3350	message_string = 'surface_reflections is allowed only when ' // &
3351	'radiation_interactions_on is set to TRUE'
3352	CALL message( 'radiation_parin', 'PA0293',1, 2, 0, 6, 0 )
3353	ENDIF
3354
3355	GOTO 14
3356
3357	13 BACKSPACE( 11 )
3358	READ( 11 , '(A)') line
3359	CALL parin_fail_message( 'radiation_par', line )
3360
3361	14 CONTINUE
3362
3363	END SUBROUTINE radiation_parin
3364
3365
3366	!------------------------------------------------------------------------------!
3367	! Description:
3368	! ------------
3369	!> Implementation of the RRTMG radiation_scheme
3370	!------------------------------------------------------------------------------!
3371	SUBROUTINE radiation_rrtmg
3372
3373	#if defined ( __rrtmg )
3374	USE indices, &
3375	ONLY: nbgp
3376
3377	USE particle_attributes, &
3378	ONLY: grid_particles, number_of_particles, particles, prt_count
3379
3380	IMPLICIT NONE
3381
3382
3383	INTEGER(iwp) :: i, j, k, l, m, n !< loop indices
3384	INTEGER(iwp) :: k_topo !< topography top index
3385
3386	REAL(wp) :: nc_rad, & !< number concentration of cloud droplets
3387	s_r2, & !< weighted sum over all droplets with r^2
3388	s_r3 !< weighted sum over all droplets with r^3
3389
3390	REAL(wp), DIMENSION(0:nzt+1) :: pt_av, q_av, ql_av
3391	REAL(wp), DIMENSION(0:0) :: zenith !< to provide indexed array
3392	!
3393	!-- Just dummy arguments
3394	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: rrtm_lw_taucld_dum, &
3395	rrtm_lw_tauaer_dum, &
3396	rrtm_sw_taucld_dum, &
3397	rrtm_sw_ssacld_dum, &
3398	rrtm_sw_asmcld_dum, &
3399	rrtm_sw_fsfcld_dum, &
3400	rrtm_sw_tauaer_dum, &
3401	rrtm_sw_ssaaer_dum, &
3402	rrtm_sw_asmaer_dum, &
3403	rrtm_sw_ecaer_dum
3404
3405	!
3406	!-- Calculate current (cosine of) zenith angle and whether the sun is up
3407	CALL calc_zenith
3408	zenith(0) = cos_zenith
3409	!
3410	!-- Calculate surface albedo. In case average radiation is applied,
3411	!-- this is not required.
3412	#if defined( __netcdf )
3413	IF ( .NOT. constant_albedo ) THEN
3414	!
3415	!-- Horizontally aligned default, natural and urban surfaces
3416	CALL calc_albedo( surf_lsm_h )
3417	CALL calc_albedo( surf_usm_h )
3418	!
3419	!-- Vertically aligned default, natural and urban surfaces
3420	DO l = 0, 3
3421	CALL calc_albedo( surf_lsm_v(l) )
3422	CALL calc_albedo( surf_usm_v(l) )
3423	ENDDO
3424	ENDIF
3425	#endif
3426
3427	!
3428	!-- Prepare input data for RRTMG
3429
3430	!
3431	!-- In case of large scale forcing with surface data, calculate new pressure
3432	!-- profile. nzt_rad might be modified by these calls and all required arrays
3433	!-- will then be re-allocated
3434	IF ( large_scale_forcing .AND. lsf_surf ) THEN
3435	CALL read_sounding_data
3436	CALL read_trace_gas_data
3437	ENDIF
3438
3439
3440	IF ( average_radiation ) THEN
3441
3442	rrtm_asdir(1) = albedo_urb
3443	rrtm_asdif(1) = albedo_urb
3444	rrtm_aldir(1) = albedo_urb
3445	rrtm_aldif(1) = albedo_urb
3446
3447	rrtm_emis = emissivity_urb
3448	!
3449	!-- Calculate mean pt profile. Actually, only one height level is required.
3450	CALL calc_mean_profile( pt, 4 )
3451	pt_av = hom(:, 1, 4, 0)
3452
3453	IF ( humidity ) THEN
3454	CALL calc_mean_profile( q, 41 )
3455	q_av = hom(:, 1, 41, 0)
3456	ENDIF
3457	!
3458	!-- Prepare profiles of temperature and H2O volume mixing ratio
3459	rrtm_tlev(0,nzb+1) = t_rad_urb
3460
3461	IF ( bulk_cloud_model ) THEN
3462
3463	CALL calc_mean_profile( ql, 54 )
3464	! average ql is now in hom(:, 1, 54, 0)
3465	ql_av = hom(:, 1, 54, 0)
3466
3467	DO k = nzb+1, nzt+1
3468	rrtm_tlay(0,k) = pt_av(k) * ( (hyp(k) ) / 100000._wp &
3469	)*.286_wp + lv_d_cp ql_av(k)
3470	rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * (q_av(k) - ql_av(k))
3471	ENDDO
3472	ELSE
3473	DO k = nzb+1, nzt+1
3474	rrtm_tlay(0,k) = pt_av(k) * ( (hyp(k) ) / 100000._wp &
3475	)**.286_wp
3476	ENDDO
3477
3478	IF ( humidity ) THEN
3479	DO k = nzb+1, nzt+1
3480	rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * q_av(k)
3481	ENDDO
3482	ELSE
3483	rrtm_h2ovmr(0,nzb+1:nzt+1) = 0.0_wp
3484	ENDIF
3485	ENDIF
3486
3487	!
3488	!-- Avoid temperature/humidity jumps at the top of the LES domain by
3489	!-- linear interpolation from nzt+2 to nzt+7
3490	DO k = nzt+2, nzt+7
3491	rrtm_tlay(0,k) = rrtm_tlay(0,nzt+1) &
3492	+ ( rrtm_tlay(0,nzt+8) - rrtm_tlay(0,nzt+1) ) &
3493	/ ( rrtm_play(0,nzt+8) - rrtm_play(0,nzt+1) ) &
3494	* ( rrtm_play(0,k) - rrtm_play(0,nzt+1) )
3495
3496	rrtm_h2ovmr(0,k) = rrtm_h2ovmr(0,nzt+1) &
3497	+ ( rrtm_h2ovmr(0,nzt+8) - rrtm_h2ovmr(0,nzt+1) )&
3498	/ ( rrtm_play(0,nzt+8) - rrtm_play(0,nzt+1) )&
3499	* ( rrtm_play(0,k) - rrtm_play(0,nzt+1) )
3500
3501	ENDDO
3502
3503	!-- Linear interpolate to zw grid
3504	DO k = nzb+2, nzt+8
3505	rrtm_tlev(0,k) = rrtm_tlay(0,k-1) + (rrtm_tlay(0,k) - &
3506	rrtm_tlay(0,k-1)) &
3507	/ ( rrtm_play(0,k) - rrtm_play(0,k-1) ) &
3508	* ( rrtm_plev(0,k) - rrtm_play(0,k-1) )
3509	ENDDO
3510
3511
3512	!
3513	!-- Calculate liquid water path and cloud fraction for each column.
3514	!-- Note that LWP is required in g/m2 instead of kg/kg m.
3515	rrtm_cldfr = 0.0_wp
3516	rrtm_reliq = 0.0_wp
3517	rrtm_cliqwp = 0.0_wp
3518	rrtm_icld = 0
3519
3520	IF ( bulk_cloud_model ) THEN
3521	DO k = nzb+1, nzt+1
3522	rrtm_cliqwp(0,k) = ql_av(k) * 1000._wp * &
3523	(rrtm_plev(0,k) - rrtm_plev(0,k+1)) &
3524	* 100._wp / g
3525
3526	IF ( rrtm_cliqwp(0,k) > 0._wp ) THEN
3527	rrtm_cldfr(0,k) = 1._wp
3528	IF ( rrtm_icld == 0 ) rrtm_icld = 1
3529
3530	!
3531	!-- Calculate cloud droplet effective radius
3532	rrtm_reliq(0,k) = 1.0E6_wp * ( 3.0_wp * ql_av(k) &
3533	* rho_surface &
3534	/ ( 4.0_wp * pi * nc_const * rho_l ) &
3535	)**0.33333333333333_wp &
3536	* EXP( LOG( sigma_gc )**2 )
3537	!
3538	!-- Limit effective radius
3539	IF ( rrtm_reliq(0,k) > 0.0_wp ) THEN
3540	rrtm_reliq(0,k) = MAX(rrtm_reliq(0,k),2.5_wp)
3541	rrtm_reliq(0,k) = MIN(rrtm_reliq(0,k),60.0_wp)
3542	ENDIF
3543	ENDIF
3544	ENDDO
3545	ENDIF
3546
3547	!
3548	!-- Set surface temperature
3549	rrtm_tsfc = t_rad_urb
3550
3551	IF ( lw_radiation ) THEN
3552
3553	CALL rrtmg_lw( 1, nzt_rad , rrtm_icld , rrtm_idrv ,&
3554	rrtm_play , rrtm_plev , rrtm_tlay , rrtm_tlev ,&
3555	rrtm_tsfc , rrtm_h2ovmr , rrtm_o3vmr , rrtm_co2vmr ,&
3556	rrtm_ch4vmr , rrtm_n2ovmr , rrtm_o2vmr , rrtm_cfc11vmr ,&
3557	rrtm_cfc12vmr , rrtm_cfc22vmr, rrtm_ccl4vmr , rrtm_emis ,&
3558	rrtm_inflglw , rrtm_iceflglw, rrtm_liqflglw, rrtm_cldfr ,&
3559	rrtm_lw_taucld , rrtm_cicewp , rrtm_cliqwp , rrtm_reice ,&
3560	rrtm_reliq , rrtm_lw_tauaer, &
3561	rrtm_lwuflx , rrtm_lwdflx , rrtm_lwhr , &
3562	rrtm_lwuflxc , rrtm_lwdflxc , rrtm_lwhrc , &
3563	rrtm_lwuflx_dt , rrtm_lwuflxc_dt )
3564
3565	!
3566	!-- Save fluxes
3567	DO k = nzb, nzt+1
3568	rad_lw_in(k,:,:) = rrtm_lwdflx(0,k)
3569	rad_lw_out(k,:,:) = rrtm_lwuflx(0,k)
3570	ENDDO
3571	rad_lw_in_diff(:,:) = rad_lw_in(0,:,:)
3572	!
3573	!-- Save heating rates (convert from K/d to K/h).
3574	!-- Further, even though an aggregated radiation is computed, map
3575	!-- signle-column profiles on top of any topography, in order to
3576	!-- obtain correct near surface radiation heating/cooling rates.
3577	DO i = nxl, nxr
3578	DO j = nys, nyn
3579	k_topo = get_topography_top_index_ji( j, i, 's' )
3580	DO k = k_topo+1, nzt+1
3581	rad_lw_hr(k,j,i) = rrtm_lwhr(0,k-k_topo) * d_hours_day
3582	rad_lw_cs_hr(k,j,i) = rrtm_lwhrc(0,k-k_topo) * d_hours_day
3583	ENDDO
3584	ENDDO
3585	ENDDO
3586
3587	ENDIF
3588
3589	IF ( sw_radiation .AND. sun_up ) THEN
3590	CALL rrtmg_sw( 1, nzt_rad , rrtm_icld , rrtm_iaer ,&
3591	rrtm_play , rrtm_plev , rrtm_tlay , rrtm_tlev ,&
3592	rrtm_tsfc , rrtm_h2ovmr , rrtm_o3vmr , rrtm_co2vmr ,&
3593	rrtm_ch4vmr , rrtm_n2ovmr , rrtm_o2vmr , rrtm_asdir ,&
3594	rrtm_asdif , rrtm_aldir , rrtm_aldif , zenith ,&
3595	0.0_wp , day_of_year , solar_constant, rrtm_inflgsw ,&
3596	rrtm_iceflgsw , rrtm_liqflgsw , rrtm_cldfr , rrtm_sw_taucld ,&
3597	rrtm_sw_ssacld , rrtm_sw_asmcld, rrtm_sw_fsfcld, rrtm_cicewp ,&
3598	rrtm_cliqwp , rrtm_reice , rrtm_reliq , rrtm_sw_tauaer ,&
3599	rrtm_sw_ssaaer , rrtm_sw_asmaer, rrtm_sw_ecaer , rrtm_swuflx ,&
3600	rrtm_swdflx , rrtm_swhr , rrtm_swuflxc , rrtm_swdflxc ,&
3601	rrtm_swhrc , rrtm_dirdflux , rrtm_difdflux )
3602
3603	!
3604	!-- Save fluxes:
3605	!-- - whole domain
3606	DO k = nzb, nzt+1
3607	rad_sw_in(k,:,:) = rrtm_swdflx(0,k)
3608	rad_sw_out(k,:,:) = rrtm_swuflx(0,k)
3609	ENDDO
3610	!-- - direct and diffuse SW at urban-surface-layer (required by RTM)
3611	rad_sw_in_dir(:,:) = rrtm_dirdflux(0,nzb)
3612	rad_sw_in_diff(:,:) = rrtm_difdflux(0,nzb)
3613
3614	!
3615	!-- Save heating rates (convert from K/d to K/s)
3616	DO k = nzb+1, nzt+1
3617	rad_sw_hr(k,:,:) = rrtm_swhr(0,k) * d_hours_day
3618	rad_sw_cs_hr(k,:,:) = rrtm_swhrc(0,k) * d_hours_day
3619	ENDDO
3620	!
3621	!-- Solar radiation is zero during night
3622	ELSE
3623	rad_sw_in = 0.0_wp
3624	rad_sw_out = 0.0_wp
3625	rad_sw_in_dir(:,:) = 0.0_wp
3626	rad_sw_in_diff(:,:) = 0.0_wp
3627	ENDIF
3628	!
3629	!-- RRTMG is called for each (j,i) grid point separately, starting at the
3630	!-- highest topography level. Here no RTM is used since average_radiation is false
3631	ELSE
3632	!
3633	!-- Loop over all grid points
3634	DO i = nxl, nxr
3635	DO j = nys, nyn
3636
3637	!
3638	!-- Prepare profiles of temperature and H2O volume mixing ratio
3639	DO m = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
3640	rrtm_tlev(0,nzb+1) = surf_lsm_h%pt_surface(m) * exner(nzb)
3641	ENDDO
3642	DO m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
3643	rrtm_tlev(0,nzb+1) = surf_usm_h%pt_surface(m) * exner(nzb)
3644	ENDDO
3645
3646
3647	IF ( bulk_cloud_model ) THEN
3648	DO k = nzb+1, nzt+1
3649	rrtm_tlay(0,k) = pt(k,j,i) * exner(k) &
3650	+ lv_d_cp * ql(k,j,i)
3651	rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * (q(k,j,i) - ql(k,j,i))
3652	ENDDO
3653	ELSEIF ( cloud_droplets ) THEN
3654	DO k = nzb+1, nzt+1
3655	rrtm_tlay(0,k) = pt(k,j,i) * exner(k) &
3656	+ lv_d_cp * ql(k,j,i)
3657	rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * q(k,j,i)
3658	ENDDO
3659	ELSE
3660	DO k = nzb+1, nzt+1
3661	rrtm_tlay(0,k) = pt(k,j,i) * exner(k)
3662	ENDDO
3663
3664	IF ( humidity ) THEN
3665	DO k = nzb+1, nzt+1
3666	rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * q(k,j,i)
3667	ENDDO
3668	ELSE
3669	rrtm_h2ovmr(0,nzb+1:nzt+1) = 0.0_wp
3670	ENDIF
3671	ENDIF
3672
3673	!
3674	!-- Avoid temperature/humidity jumps at the top of the LES domain by
3675	!-- linear interpolation from nzt+2 to nzt+7
3676	DO k = nzt+2, nzt+7
3677	rrtm_tlay(0,k) = rrtm_tlay(0,nzt+1) &
3678	+ ( rrtm_tlay(0,nzt+8) - rrtm_tlay(0,nzt+1) ) &
3679	/ ( rrtm_play(0,nzt+8) - rrtm_play(0,nzt+1) ) &
3680	* ( rrtm_play(0,k) - rrtm_play(0,nzt+1) )
3681
3682	rrtm_h2ovmr(0,k) = rrtm_h2ovmr(0,nzt+1) &
3683	+ ( rrtm_h2ovmr(0,nzt+8) - rrtm_h2ovmr(0,nzt+1) )&
3684	/ ( rrtm_play(0,nzt+8) - rrtm_play(0,nzt+1) )&
3685	* ( rrtm_play(0,k) - rrtm_play(0,nzt+1) )
3686
3687	ENDDO
3688
3689	!-- Linear interpolate to zw grid
3690	DO k = nzb+2, nzt+8
3691	rrtm_tlev(0,k) = rrtm_tlay(0,k-1) + (rrtm_tlay(0,k) - &
3692	rrtm_tlay(0,k-1)) &
3693	/ ( rrtm_play(0,k) - rrtm_play(0,k-1) ) &
3694	* ( rrtm_plev(0,k) - rrtm_play(0,k-1) )
3695	ENDDO
3696
3697
3698	!
3699	!-- Calculate liquid water path and cloud fraction for each column.
3700	!-- Note that LWP is required in g/m2 instead of kg/kg m.
3701	rrtm_cldfr = 0.0_wp
3702	rrtm_reliq = 0.0_wp
3703	rrtm_cliqwp = 0.0_wp
3704	rrtm_icld = 0
3705
3706	IF ( bulk_cloud_model .OR. cloud_droplets ) THEN
3707	DO k = nzb+1, nzt+1
3708	rrtm_cliqwp(0,k) = ql(k,j,i) * 1000.0_wp * &
3709	(rrtm_plev(0,k) - rrtm_plev(0,k+1)) &
3710	* 100.0_wp / g
3711
3712	IF ( rrtm_cliqwp(0,k) > 0.0_wp ) THEN
3713	rrtm_cldfr(0,k) = 1.0_wp
3714	IF ( rrtm_icld == 0 ) rrtm_icld = 1
3715
3716	!
3717	!-- Calculate cloud droplet effective radius
3718	IF ( bulk_cloud_model ) THEN
3719	!
3720	!-- Calculete effective droplet radius. In case of using
3721	!-- cloud_scheme = 'morrison' and a non reasonable number
3722	!-- of cloud droplets the inital aerosol number
3723	!-- concentration is considered.
3724	IF ( microphysics_morrison ) THEN
3725	IF ( nc(k,j,i) > 1.0E-20_wp ) THEN
3726	nc_rad = nc(k,j,i)
3727	ELSE
3728	nc_rad = na_init
3729	ENDIF
3730	ELSE
3731	nc_rad = nc_const
3732	ENDIF
3733
3734	rrtm_reliq(0,k) = 1.0E6_wp * ( 3.0_wp * ql(k,j,i) &
3735	* rho_surface &
3736	/ ( 4.0_wp * pi * nc_rad * rho_l ) &
3737	)**0.33333333333333_wp &
3738	* EXP( LOG( sigma_gc )**2 )
3739
3740	ELSEIF ( cloud_droplets ) THEN
3741	number_of_particles = prt_count(k,j,i)
3742
3743	IF (number_of_particles <= 0) CYCLE
3744	particles => grid_particles(k,j,i)%particles(1:number_of_particles)
3745	s_r2 = 0.0_wp
3746	s_r3 = 0.0_wp
3747
3748	DO n = 1, number_of_particles
3749	IF ( particles(n)%particle_mask ) THEN
3750	s_r2 = s_r2 + particles(n)%radius*2 &
3751	particles(n)%weight_factor
3752	s_r3 = s_r3 + particles(n)%radius*3 &
3753	particles(n)%weight_factor
3754	ENDIF
3755	ENDDO
3756
3757	IF ( s_r2 > 0.0_wp ) rrtm_reliq(0,k) = s_r3 / s_r2
3758
3759	ENDIF
3760
3761	!
3762	!-- Limit effective radius
3763	IF ( rrtm_reliq(0,k) > 0.0_wp ) THEN
3764	rrtm_reliq(0,k) = MAX(rrtm_reliq(0,k),2.5_wp)
3765	rrtm_reliq(0,k) = MIN(rrtm_reliq(0,k),60.0_wp)
3766	ENDIF
3767	ENDIF
3768	ENDDO
3769	ENDIF
3770
3771	!
3772	!-- Write surface emissivity and surface temperature at current
3773	!-- surface element on RRTMG-shaped array.
3774	!-- Please note, as RRTMG is a single column model, surface attributes
3775	!-- are only obtained from horizontally aligned surfaces (for
3776	!-- simplicity). Taking surface attributes from horizontal and
3777	!-- vertical walls would lead to multiple solutions.
3778	!-- Moreover, for natural- and urban-type surfaces, several surface
3779	!-- classes can exist at a surface element next to each other.
3780	!-- To obtain bulk parameters, apply a weighted average for these
3781	!-- surfaces.
3782	DO m = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
3783	rrtm_emis = surf_lsm_h%frac(ind_veg_wall,m) * &
3784	surf_lsm_h%emissivity(ind_veg_wall,m) + &
3785	surf_lsm_h%frac(ind_pav_green,m) * &
3786	surf_lsm_h%emissivity(ind_pav_green,m) + &
3787	surf_lsm_h%frac(ind_wat_win,m) * &
3788	surf_lsm_h%emissivity(ind_wat_win,m)
3789	rrtm_tsfc = surf_lsm_h%pt_surface(m) * exner(nzb)
3790	ENDDO
3791	DO m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
3792	rrtm_emis = surf_usm_h%frac(ind_veg_wall,m) * &
3793	surf_usm_h%emissivity(ind_veg_wall,m) + &
3794	surf_usm_h%frac(ind_pav_green,m) * &
3795	surf_usm_h%emissivity(ind_pav_green,m) + &
3796	surf_usm_h%frac(ind_wat_win,m) * &
3797	surf_usm_h%emissivity(ind_wat_win,m)
3798	rrtm_tsfc = surf_usm_h%pt_surface(m) * exner(nzb)
3799	ENDDO
3800	!
3801	!-- Obtain topography top index (lower bound of RRTMG)
3802	k_topo = get_topography_top_index_ji( j, i, 's' )
3803
3804	IF ( lw_radiation ) THEN
3805	!
3806	!-- Due to technical reasons, copy optical depth to dummy arguments
3807	!-- which are allocated on the exact size as the rrtmg_lw is called.
3808	!-- As one dimesion is allocated with zero size, compiler complains
3809	!-- that rank of the array does not match that of the
3810	!-- assumed-shaped arguments in the RRTMG library. In order to
3811	!-- avoid this, write to dummy arguments and give pass the entire
3812	!-- dummy array. Seems to be the only existing work-around.
3813	ALLOCATE( rrtm_lw_taucld_dum(1:nbndlw+1,0:0,k_topo+1:nzt_rad+1) )
3814	ALLOCATE( rrtm_lw_tauaer_dum(0:0,k_topo+1:nzt_rad+1,1:nbndlw+1) )
3815
3816	rrtm_lw_taucld_dum = &
3817	rrtm_lw_taucld(1:nbndlw+1,0:0,k_topo+1:nzt_rad+1)
3818	rrtm_lw_tauaer_dum = &
3819	rrtm_lw_tauaer(0:0,k_topo+1:nzt_rad+1,1:nbndlw+1)
3820
3821	CALL rrtmg_lw( 1, &
3822	nzt_rad-k_topo, &
3823	rrtm_icld, &
3824	rrtm_idrv, &
3825	rrtm_play(:,k_topo+1:nzt_rad+1), &
3826	rrtm_plev(:,k_topo+1:nzt_rad+2), &
3827	rrtm_tlay(:,k_topo+1:nzt_rad+1), &
3828	rrtm_tlev(:,k_topo+1:nzt_rad+2), &
3829	rrtm_tsfc, &
3830	rrtm_h2ovmr(:,k_topo+1:nzt_rad+1), &
3831	rrtm_o3vmr(:,k_topo+1:nzt_rad+1), &
3832	rrtm_co2vmr(:,k_topo+1:nzt_rad+1), &
3833	rrtm_ch4vmr(:,k_topo+1:nzt_rad+1), &
3834	rrtm_n2ovmr(:,k_topo+1:nzt_rad+1), &
3835	rrtm_o2vmr(:,k_topo+1:nzt_rad+1), &
3836	rrtm_cfc11vmr(:,k_topo+1:nzt_rad+1), &
3837	rrtm_cfc12vmr(:,k_topo+1:nzt_rad+1), &
3838	rrtm_cfc22vmr(:,k_topo+1:nzt_rad+1), &
3839	rrtm_ccl4vmr(:,k_topo+1:nzt_rad+1), &
3840	rrtm_emis, &
3841	rrtm_inflglw, &
3842	rrtm_iceflglw, &
3843	rrtm_liqflglw, &
3844	rrtm_cldfr(:,k_topo+1:nzt_rad+1), &
3845	rrtm_lw_taucld_dum, &
3846	rrtm_cicewp(:,k_topo+1:nzt_rad+1), &
3847	rrtm_cliqwp(:,k_topo+1:nzt_rad+1), &
3848	rrtm_reice(:,k_topo+1:nzt_rad+1), &
3849	rrtm_reliq(:,k_topo+1:nzt_rad+1), &
3850	rrtm_lw_tauaer_dum, &
3851	rrtm_lwuflx(:,k_topo:nzt_rad+1), &
3852	rrtm_lwdflx(:,k_topo:nzt_rad+1), &
3853	rrtm_lwhr(:,k_topo+1:nzt_rad+1), &
3854	rrtm_lwuflxc(:,k_topo:nzt_rad+1), &
3855	rrtm_lwdflxc(:,k_topo:nzt_rad+1), &
3856	rrtm_lwhrc(:,k_topo+1:nzt_rad+1), &
3857	rrtm_lwuflx_dt(:,k_topo:nzt_rad+1), &
3858	rrtm_lwuflxc_dt(:,k_topo:nzt_rad+1) )
3859
3860	DEALLOCATE ( rrtm_lw_taucld_dum )
3861	DEALLOCATE ( rrtm_lw_tauaer_dum )
3862	!
3863	!-- Save fluxes
3864	DO k = k_topo, nzt+1
3865	rad_lw_in(k,j,i) = rrtm_lwdflx(0,k)
3866	rad_lw_out(k,j,i) = rrtm_lwuflx(0,k)
3867	ENDDO
3868
3869	!
3870	!-- Save heating rates (convert from K/d to K/h)
3871	DO k = k_topo+1, nzt+1
3872	rad_lw_hr(k,j,i) = rrtm_lwhr(0,k-k_topo) * d_hours_day
3873	rad_lw_cs_hr(k,j,i) = rrtm_lwhrc(0,k-k_topo) * d_hours_day
3874	ENDDO
3875
3876	!
3877	!-- Save surface radiative fluxes and change in LW heating rate
3878	!-- onto respective surface elements
3879	!-- Horizontal surfaces
3880	DO m = surf_lsm_h%start_index(j,i), &
3881	surf_lsm_h%end_index(j,i)
3882	surf_lsm_h%rad_lw_in(m) = rrtm_lwdflx(0,k_topo)
3883	surf_lsm_h%rad_lw_out(m) = rrtm_lwuflx(0,k_topo)
3884	surf_lsm_h%rad_lw_out_change_0(m) = rrtm_lwuflx_dt(0,k_topo)
3885	ENDDO
3886	DO m = surf_usm_h%start_index(j,i), &
3887	surf_usm_h%end_index(j,i)
3888	surf_usm_h%rad_lw_in(m) = rrtm_lwdflx(0,k_topo)
3889	surf_usm_h%rad_lw_out(m) = rrtm_lwuflx(0,k_topo)
3890	surf_usm_h%rad_lw_out_change_0(m) = rrtm_lwuflx_dt(0,k_topo)
3891	ENDDO
3892	!
3893	!-- Vertical surfaces. Fluxes are obtain at vertical level of the
3894	!-- respective surface element
3895	DO l = 0, 3
3896	DO m = surf_lsm_v(l)%start_index(j,i), &
3897	surf_lsm_v(l)%end_index(j,i)
3898	k = surf_lsm_v(l)%k(m)
3899	surf_lsm_v(l)%rad_lw_in(m) = rrtm_lwdflx(0,k)
3900	surf_lsm_v(l)%rad_lw_out(m) = rrtm_lwuflx(0,k)
3901	surf_lsm_v(l)%rad_lw_out_change_0(m) = rrtm_lwuflx_dt(0,k)
3902	ENDDO
3903	DO m = surf_usm_v(l)%start_index(j,i), &
3904	surf_usm_v(l)%end_index(j,i)
3905	k = surf_usm_v(l)%k(m)
3906	surf_usm_v(l)%rad_lw_in(m) = rrtm_lwdflx(0,k)
3907	surf_usm_v(l)%rad_lw_out(m) = rrtm_lwuflx(0,k)
3908	surf_usm_v(l)%rad_lw_out_change_0(m) = rrtm_lwuflx_dt(0,k)
3909	ENDDO
3910	ENDDO
3911
3912	ENDIF
3913
3914	IF ( sw_radiation .AND. sun_up ) THEN
3915	!
3916	!-- Get albedo for direct/diffusive long/shortwave radiation at
3917	!-- current (y,x)-location from surface variables.
3918	!-- Only obtain it from horizontal surfaces, as RRTMG is a single
3919	!-- column model
3920	!-- (Please note, only one loop will entered, controlled by
3921	!-- start-end index.)
3922	DO m = surf_lsm_h%start_index(j,i), &
3923	surf_lsm_h%end_index(j,i)
3924	rrtm_asdir(1) = SUM( surf_lsm_h%frac(:,m) * &
3925	surf_lsm_h%rrtm_asdir(:,m) )
3926	rrtm_asdif(1) = SUM( surf_lsm_h%frac(:,m) * &
3927	surf_lsm_h%rrtm_asdif(:,m) )
3928	rrtm_aldir(1) = SUM( surf_lsm_h%frac(:,m) * &
3929	surf_lsm_h%rrtm_aldir(:,m) )
3930	rrtm_aldif(1) = SUM( surf_lsm_h%frac(:,m) * &
3931	surf_lsm_h%rrtm_aldif(:,m) )
3932	ENDDO
3933	DO m = surf_usm_h%start_index(j,i), &
3934	surf_usm_h%end_index(j,i)
3935	rrtm_asdir(1) = SUM( surf_usm_h%frac(:,m) * &
3936	surf_usm_h%rrtm_asdir(:,m) )
3937	rrtm_asdif(1) = SUM( surf_usm_h%frac(:,m) * &
3938	surf_usm_h%rrtm_asdif(:,m) )
3939	rrtm_aldir(1) = SUM( surf_usm_h%frac(:,m) * &
3940	surf_usm_h%rrtm_aldir(:,m) )
3941	rrtm_aldif(1) = SUM( surf_usm_h%frac(:,m) * &
3942	surf_usm_h%rrtm_aldif(:,m) )
3943	ENDDO
3944	!
3945	!-- Due to technical reasons, copy optical depths and other
3946	!-- to dummy arguments which are allocated on the exact size as the
3947	!-- rrtmg_sw is called.
3948	!-- As one dimesion is allocated with zero size, compiler complains
3949	!-- that rank of the array does not match that of the
3950	!-- assumed-shaped arguments in the RRTMG library. In order to
3951	!-- avoid this, write to dummy arguments and give pass the entire
3952	!-- dummy array. Seems to be the only existing work-around.
3953	ALLOCATE( rrtm_sw_taucld_dum(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1) )
3954	ALLOCATE( rrtm_sw_ssacld_dum(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1) )
3955	ALLOCATE( rrtm_sw_asmcld_dum(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1) )
3956	ALLOCATE( rrtm_sw_fsfcld_dum(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1) )
3957	ALLOCATE( rrtm_sw_tauaer_dum(0:0,k_topo+1:nzt_rad+1,1:nbndsw+1) )
3958	ALLOCATE( rrtm_sw_ssaaer_dum(0:0,k_topo+1:nzt_rad+1,1:nbndsw+1) )
3959	ALLOCATE( rrtm_sw_asmaer_dum(0:0,k_topo+1:nzt_rad+1,1:nbndsw+1) )
3960	ALLOCATE( rrtm_sw_ecaer_dum(0:0,k_topo+1:nzt_rad+1,1:naerec+1) )
3961
3962	rrtm_sw_taucld_dum = rrtm_sw_taucld(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1)
3963	rrtm_sw_ssacld_dum = rrtm_sw_ssacld(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1)
3964	rrtm_sw_asmcld_dum = rrtm_sw_asmcld(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1)
3965	rrtm_sw_fsfcld_dum = rrtm_sw_fsfcld(1:nbndsw+1,0:0,k_topo+1:nzt_rad+1)
3966	rrtm_sw_tauaer_dum = rrtm_sw_tauaer(0:0,k_topo+1:nzt_rad+1,1:nbndsw+1)
3967	rrtm_sw_ssaaer_dum = rrtm_sw_ssaaer(0:0,k_topo+1:nzt_rad+1,1:nbndsw+1)
3968	rrtm_sw_asmaer_dum = rrtm_sw_asmaer(0:0,k_topo+1:nzt_rad+1,1:nbndsw+1)
3969	rrtm_sw_ecaer_dum = rrtm_sw_ecaer(0:0,k_topo+1:nzt_rad+1,1:naerec+1)
3970
3971	CALL rrtmg_sw( 1, &
3972	nzt_rad-k_topo, &
3973	rrtm_icld, &
3974	rrtm_iaer, &
3975	rrtm_play(:,k_topo+1:nzt_rad+1), &
3976	rrtm_plev(:,k_topo+1:nzt_rad+2), &
3977	rrtm_tlay(:,k_topo+1:nzt_rad+1), &
3978	rrtm_tlev(:,k_topo+1:nzt_rad+2), &
3979	rrtm_tsfc, &
3980	rrtm_h2ovmr(:,k_topo+1:nzt_rad+1), &
3981	rrtm_o3vmr(:,k_topo+1:nzt_rad+1), &
3982	rrtm_co2vmr(:,k_topo+1:nzt_rad+1), &
3983	rrtm_ch4vmr(:,k_topo+1:nzt_rad+1), &
3984	rrtm_n2ovmr(:,k_topo+1:nzt_rad+1), &
3985	rrtm_o2vmr(:,k_topo+1:nzt_rad+1), &
3986	rrtm_asdir, &
3987	rrtm_asdif, &
3988	rrtm_aldir, &
3989	rrtm_aldif, &
3990	zenith, &
3991	0.0_wp, &
3992	day_of_year, &
3993	solar_constant, &
3994	rrtm_inflgsw, &
3995	rrtm_iceflgsw, &
3996	rrtm_liqflgsw, &
3997	rrtm_cldfr(:,k_topo+1:nzt_rad+1), &
3998	rrtm_sw_taucld_dum, &
3999	rrtm_sw_ssacld_dum, &
4000	rrtm_sw_asmcld_dum, &
4001	rrtm_sw_fsfcld_dum, &
4002	rrtm_cicewp(:,k_topo+1:nzt_rad+1), &
4003	rrtm_cliqwp(:,k_topo+1:nzt_rad+1), &
4004	rrtm_reice(:,k_topo+1:nzt_rad+1), &
4005	rrtm_reliq(:,k_topo+1:nzt_rad+1), &
4006	rrtm_sw_tauaer_dum, &
4007	rrtm_sw_ssaaer_dum, &
4008	rrtm_sw_asmaer_dum, &
4009	rrtm_sw_ecaer_dum, &
4010	rrtm_swuflx(:,k_topo:nzt_rad+1), &
4011	rrtm_swdflx(:,k_topo:nzt_rad+1), &
4012	rrtm_swhr(:,k_topo+1:nzt_rad+1), &
4013	rrtm_swuflxc(:,k_topo:nzt_rad+1), &
4014	rrtm_swdflxc(:,k_topo:nzt_rad+1), &
4015	rrtm_swhrc(:,k_topo+1:nzt_rad+1), &
4016	rrtm_dirdflux(:,k_topo:nzt_rad+1), &
4017	rrtm_difdflux(:,k_topo:nzt_rad+1) )
4018
4019	DEALLOCATE( rrtm_sw_taucld_dum )
4020	DEALLOCATE( rrtm_sw_ssacld_dum )
4021	DEALLOCATE( rrtm_sw_asmcld_dum )
4022	DEALLOCATE( rrtm_sw_fsfcld_dum )
4023	DEALLOCATE( rrtm_sw_tauaer_dum )
4024	DEALLOCATE( rrtm_sw_ssaaer_dum )
4025	DEALLOCATE( rrtm_sw_asmaer_dum )
4026	DEALLOCATE( rrtm_sw_ecaer_dum )
4027	!
4028	!-- Save fluxes
4029	DO k = nzb, nzt+1
4030	rad_sw_in(k,j,i) = rrtm_swdflx(0,k)
4031	rad_sw_out(k,j,i) = rrtm_swuflx(0,k)
4032	ENDDO
4033	!
4034	!-- Save heating rates (convert from K/d to K/s)
4035	DO k = nzb+1, nzt+1
4036	rad_sw_hr(k,j,i) = rrtm_swhr(0,k) * d_hours_day
4037	rad_sw_cs_hr(k,j,i) = rrtm_swhrc(0,k) * d_hours_day
4038	ENDDO
4039
4040	!
4041	!-- Save surface radiative fluxes onto respective surface elements
4042	!-- Horizontal surfaces
4043	DO m = surf_lsm_h%start_index(j,i), &
4044	surf_lsm_h%end_index(j,i)
4045	surf_lsm_h%rad_sw_in(m) = rrtm_swdflx(0,k_topo)
4046	surf_lsm_h%rad_sw_out(m) = rrtm_swuflx(0,k_topo)
4047	ENDDO
4048	DO m = surf_usm_h%start_index(j,i), &
4049	surf_usm_h%end_index(j,i)
4050	surf_usm_h%rad_sw_in(m) = rrtm_swdflx(0,k_topo)
4051	surf_usm_h%rad_sw_out(m) = rrtm_swuflx(0,k_topo)
4052	ENDDO
4053	!
4054	!-- Vertical surfaces. Fluxes are obtain at respective vertical
4055	!-- level of the surface element
4056	DO l = 0, 3
4057	DO m = surf_lsm_v(l)%start_index(j,i), &
4058	surf_lsm_v(l)%end_index(j,i)
4059	k = surf_lsm_v(l)%k(m)
4060	surf_lsm_v(l)%rad_sw_in(m) = rrtm_swdflx(0,k)
4061	surf_lsm_v(l)%rad_sw_out(m) = rrtm_swuflx(0,k)
4062	ENDDO
4063	DO m = surf_usm_v(l)%start_index(j,i), &
4064	surf_usm_v(l)%end_index(j,i)
4065	k = surf_usm_v(l)%k(m)
4066	surf_usm_v(l)%rad_sw_in(m