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

Last change on this file since 2967 was 2967, checked in by raasch, 3 years ago

bugfix: missing parallel cpp-directives added

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