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

Last change on this file since 3421 was 3378, checked in by kanani, 3 years ago

merge fixes of radiation branch (r3362) to trunk

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