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

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

rrtmg preprocessor for directives moved/added, save attribute added to temporary pointers to avoid compiler warnings about outlived pointer targets, statement added to avoid compiler warning about unused variable

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