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

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

Bugfix, pass dummy string to MPI_INFO_SET

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