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

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

Further revision of 2D surface output for radiation and chemistry quantities; bugfix for commit 3170

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