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

Last change on this file since 3630 was 3630, checked in by knoop, 4 years ago

Bugfix: changed ERROR STOP to STOP 1 as PGI currently can not handle it.

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