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

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

Bugfix for commit 3172

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