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

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

Code adjusted according to coding standards, renamed namelists, error messages revised until PA0347, output CASE 108 disabled

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