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

Last change on this file since 2995 was 2995, checked in by Giersch, 3 years ago

Bugfixes related to spinup and radiation model

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