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

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

diagnostic output: Modularize diagnostic output, rename subroutines; formatting adjustments; allocate arrays only when required; add output of uu, vv, ww to enable variance calculation via temporal EC method; radiation: bugfix in masked data output; flow_statistics: Correct conversion to kinematic vertical scalar fluxes in case of pw-scheme and statistic regions

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