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

Last change on this file since 3495 was 3495, checked in by kanani, 4 years ago

bugfix in calculating the apparent solar positions, resort ONLY list (radiation_model_mod)

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