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

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

fix time variables during solar positions calc

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