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

Last change on this file since 3274 was 3274, checked in by knoop, 3 years ago

Modularization of all bulk cloud physics code components

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