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

Last change on this file since 3861 was 3861, checked in by maronga, 2 years ago

bugfix in radiation model (rad_lw_out_change_0)

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