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

Last change on this file since 3524 was 3524, checked in by raasch, 4 years ago

unused variables removed, missing working precision added, missing preprocessor directives added, bugfix concerning allocation of t_surf_wall_v in nopointer case, declaration statements rearranged to avoid compile time errors, mpi_abort arguments replaced to avoid compile errors

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