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

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

Mask topography while imposing inflow perturbations at the boundaries; do not impose perturbations at top boundary as well as ghost points; Remove print statement; Read zsoil dimension lenght only if soil variables are provided in dynamic driver

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