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

Last change on this file since 3655 was 3655, checked in by knoop, 3 years ago

Bugfix: made "unit" and "found" intend INOUT in module interface subroutines + automatic copyright update

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