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

Last change on this file since 4168 was 4168, checked in by suehring, 3 years ago

Replace get_topography_top_index functions by pre-calculated arrays in order to save computational resources

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