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

Last change on this file since 3337 was 3337, checked in by kanani, 3 years ago

reintegrate branch resler to trunk

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