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

Last change on this file since 2906 was 2906, checked in by Giersch, 4 years ago

new procedure for reading/writing svf data, initialization of local variable ids

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