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

Last change on this file since 2932 was 2932, checked in by maronga, 4 years ago

renamed all Fortran NAMELISTS

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