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

Last change on this file since 2724 was 2724, checked in by maronga, 7 years ago

some changes in spinup mechanism and additional check in land surface model

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