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

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

Make level 3 initialization of urban-surfaces consistent to input data standard; revise flagging of ground-floor level surfaces at buidlings; bugfixes in level 3 initialization of albedo; bugfix in OpenMP directive; check for zero output timestep in surface output; assure that Obukhov lenght does not become zero

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