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

Last change on this file since 3992 was 3992, checked in by suehring, 2 years ago

Bugfix in rrtmg radiation branch in a nested run when the lowest prognistic grid points in a child domain are all inside topography

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