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

Last change on this file since 3172 was 3172, checked in by suehring, 7 years ago

temporal work-around for calculation of effective radiative surface temperature; prevent positive solar radiation during nighttime

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