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

Last change on this file since 2701 was 2701, checked in by suehring, 4 years ago

changes from last commit documented

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