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

Last change on this file since 3178 was 3178, checked in by suehring, 6 years ago

Revise concept for calculation of radiative temperature

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