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

Last change on this file since 2894 was 2894, checked in by Giersch, 4 years ago

Reading/Writing? data in case of restart runs revised

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