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

Last change on this file since 3464 was 3464, checked in by kanani, 3 years ago

from branch resler@3462: add MRT shaping function (radiation_model_mod), use basic constants (biometeorology_mod), adjust precision to wp (biometeorology_pt_mod)

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