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

Last change on this file since 3987 was 3987, checked in by kanani, 2 years ago

clean up location, debug and error messages

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