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

Last change on this file since 3768 was 3767, checked in by raasch, 2 years ago

unused variables removed from rrd-subroutines parameter list

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