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

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

Merge with branch radition - improved representation of diffuse radiation

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