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

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

Remove erroneous UTF encoding; last commit documented

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