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

Last change on this file since 3528 was 3528, checked in by suehring, 4 years ago

Bugfix in raytracing - add an epsilon value to overcome precision-related errors.

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