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

Last change on this file since 4188 was 4188, checked in by suehring, 2 years ago

Minor adjustments in error messages and error numbers. Some typos are corrected.

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