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

Last change on this file since 2809 was 2809, checked in by schwenkel, 6 years ago

Bugfix for gfortran: Replace the function C_SIZEOF with STORAGE_SIZE

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