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

Last change on this file since 4134 was 4134, checked in by suehring, 23 months ago

Remove faulty debug write

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