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

Last change on this file since 4495 was 4495, checked in by raasch, 14 months ago

restart data handling with MPI-IO added, first part

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