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

Last change on this file since 3246 was 3246, checked in by sward, 3 years ago

Added error handling for wrong input parameters

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