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

Last change on this file since 3685 was 3685, checked in by knoop, 3 years ago

Some interface calls moved to module_interface + cleanup

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