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

Last change on this file since 3049 was 3049, checked in by Giersch, 3 years ago

Revision history corrected

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