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

Last change on this file since 3826 was 3826, checked in by raasch, 3 years ago

unused variable removed

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