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

Last change on this file since 2941 was 2941, checked in by kanani, 7 years ago

Fixes for spinup and sky-view-factor (SVF) treatment

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