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

Last change on this file since 4069 was 4069, checked in by Giersch, 2 years ago

Bugfix for masked output, compiler warning removed, test case for wind turbine model revised

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