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

Last change on this file since 3846 was 3846, checked in by suehring, 3 years ago

Check for proper setting of dt_radiation

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