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

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

radiation: Take external radiation input from root domain dynamic file if no dynamic input file is provided for each nested domain; radiation: Combine MPI_ALLREDUCE calls to reduce latency; land_surface + plant_canopy: Give specific error numbers; land-surface: Adjust error messages for local checks; init_3d_model: Deallocate temporary string array for netcdf-data input since it may be re-used in other modules for different input files

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