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

Last change on this file since 3769 was 3769, checked in by moh.hefny, 3 years ago

removed unused variables from urban_surface_mod radiation_model_mod and part of module_interface

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