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

Last change on this file since 2716 was 2716, checked in by kanani, 4 years ago

Correction of "Former revisions" section

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