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

Last change on this file since 2920 was 2920, checked in by kanani, 5 years ago

Optimize SVF calculation, clean-up, bugfixes

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