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

Last change on this file since 3607 was 3607, checked in by suehring, 4 years ago

Output of radiation-related quantities migrated from urban_surface_model_mod to radiation_model_mod

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