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

Last change on this file since 3754 was 3754, checked in by kanani, 3 years ago

Bugfixes for calculation and I/O of view factors

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