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

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

Dollar sign added before Id; Revised structure of reading svf data according to PALM coding standard: svf_code_field/len and fsvf removed, error messages PA0493 and PA0494 added, allocation status of output arrays checked

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