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

Last change on this file since 3847 was 3847, checked in by suehring, 3 years ago

Adjust message-call for checks that are especially carried out locally.

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