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

Last change on this file since 3743 was 3743, checked in by moh.hefny, 3 years ago

added read/write number of MRT factors to the respective routines

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