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

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

Merge with branch resler: biomet- output of bio_mrt added; plant_canopy - separate vertical dimension for 3D output (to save disk space); radiation - remove unused plant canopy variables; urban-surface model - do not add anthropogenic heat during wall spin-up

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