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

Last change on this file since 3116 was 3116, checked in by suehring, 6 years ago

Output of radiation terms at surface

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