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

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

indoor_model_mod: Revision of indoor model; urban_surface_mod: parameters for waste heat from cooling and heating are introduced to building data base; initialization of building data base moved to an earlier point of time before indoor model will be initialized; radiation_model_mod: minor improvement in some comments; synthetic_turbulence_generator: unused variable removed

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