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

Last change on this file since 4178 was 4178, checked in by suehring, 2 years ago

Enable external time-dependent radiative forcing with downwelling short- and longwave radiation. Optionally, also downwelling diffuse radiation can be provided. Radiation data will be provided via dynamic input file.

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