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

Last change on this file since 3230 was 3230, checked in by schwenkel, 6 years ago

Bugfix in radiation_constant_surf

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