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

Last change on this file since 3248 was 3248, checked in by sward, 3 years ago

Minor format changes

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