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

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

chem_emissions: Replace global arrays also in mode_emis branch; diagnostic output: restructure initialization in order to work also when data output during spin-up is enabled; radiation: give informative message on raytracing distance only by core zero not by all cores

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