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

Last change on this file since 4008 was 4008, checked in by moh.hefny, 2 years ago

Bugfix in check radiation related variables

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