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

Last change on this file since 3667 was 3667, checked in by schwenkel, 3 years ago

Modified rrtmg input files check

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