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

Last change on this file since 3859 was 3859, checked in by maronga, 4 years ago

comments in radiation model updated, minor bugfix in palm_csd

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