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

Last change on this file since 3843 was 3843, checked in by pavelkrc, 3 years ago

Code review for radiation_model_mod.f90

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