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

Last change on this file since 2746 was 2746, checked in by suehring, 7 years ago

Read information from statitic driver for resolved vegetation independently from land- or urban-surface model

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