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

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

Correct level 2 initialization of spectral albedos in rrtmg branch, long- and shortwave albedos were mixed-up; Change order of albedo_pars so that it is now consistent with the defined order of albedo_pars in PIDS

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