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

Last change on this file since 3449 was 3449, checked in by suehring, 3 years ago

Branch resler -r 3439 re-integrated into current trunk: RTM 3.0, transpiration of plant canopy, output fixes in USM

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