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

Last change on this file since 3065 was 3065, checked in by Giersch, 4 years ago

New vertical stretching procedure has been introduced

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