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

Last change on this file since 3824 was 3824, checked in by pavelkrc, 3 years ago

Code review of radiation_model_mod.f90

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