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

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

Bugfix in initialization of surfaces in cyclic-fill case + bugfix in radiation output

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