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

Last change on this file since 3633 was 3633, checked in by schwenkel, 4 years ago

your commit message

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