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

Last change on this file since 3435 was 3435, checked in by gronemeier, 4 years ago

new: terrain-following masked output; bugfixes: increase vertical dimension of gamma_w_green_sat by 1, add checks for masked output for chemistry_model and radiation_model, reordered calls to xxx_define_netcdf_grid in masked output part

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