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

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

Bugfixes in calculation of sky-view and canopy-sink factors

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