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

Last change on this file since 3066 was 3066, checked in by Giersch, 3 years ago

Error messages related to the new vertical grid stretching revised, Bugfix in IF statement before error message

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