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

Last change on this file since 3814 was 3814, checked in by pavelkrc, 3 years ago

Rename conflicting variable names

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