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

Last change on this file since 3026 was 3026, checked in by schwenkel, 3 years ago

Changed the name specific humidity to mixing ratio

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