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

Last change on this file since 4179 was 4179, checked in by suehring, 2 years ago

remove debug prints

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