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

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

Offline nesting: data input modularized; time variable is defined relative to time_utc_init, so that input data is correctly mapped to actual model time; checks rephrased and new checks for the time dimension added; Netcdf input: routine to retrieve dimension length renamed

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