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

Last change on this file since 3705 was 3705, checked in by suehring, 4 years ago

last commit documented

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