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

Last change on this file since 2930 was 2930, checked in by suehring, 4 years ago

Remove default surfaces from radiation model and add check for this; revise checks for surface_fraction

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