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

Last change on this file since 3704 was 3704, checked in by suehring, 5 years ago

Revision of virtual-measurement module and data output enabled. Further, post-processing tool added to merge distributed virtually sampled data and to output it into NetCDF files.

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