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

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

Merge from branch resler: Changed behaviour of masked output over surface to follow terrain and ignore buildings

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