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

Last change on this file since 3004 was 3004, checked in by Giersch, 3 years ago

precipitation_rate removed, further allocation checks for data output of averaged quantities implemented, double CALL of flow_statistics at the beginning of time_integration removed, further minor bugfixes, comments added

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