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

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

Bugfix in the calculation of fixed number of output time levels for parallel netcdf output, error message related to reading sky view factors revised

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