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

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

Bugfixes in initialization and STG

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