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

Last change on this file since 4208 was 4208, checked in by suehring, 21 months ago

Bugfix in accessing albedo_pars in the clear-sky branch (merge from branch)

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