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

Last change on this file since 2723 was 2723, checked in by maronga, 7 years ago

bugfixes for spinup mechanism to work with lsm+usm+radiation

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