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

Last change on this file since 4245 was 4245, checked in by pavelkrc, 2 years ago

Merge branch resler into trunk

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