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

Last change on this file since 3117 was 3117, checked in by maronga, 6 years ago

bugfixes in RRTMG interface

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