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

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

Do not overwrite values of albedo in radiation_init in case albedo has been already initialized in the urban-surface model via ASCII input

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