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

Last change on this file since 3261 was 3261, checked in by suehring, 3 years ago

Bugfix in raytrace_2d calls

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