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

Last change on this file since 3825 was 3825, checked in by raasch, 3 years ago

unused variable removed

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