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

Last change on this file since 2296 was 2296, checked in by maronga, 7 years ago

added new spinup mechanism for surface/radiation models

  • Property svn:keywords set to Id
File size: 144.5 KB
Line 
1!> @file radiation_model_mod.f90
2!------------------------------------------------------------------------------!
3! This file is part of PALM.
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 1997-2017 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: radiation_model_mod.f90 2296 2017-06-28 07:53:56Z maronga $
27! Added output of rad_sw_out for radiation_scheme = 'constant'
28!
29! 2270 2017-06-09 12:18:47Z maronga
30! Numbering changed (2 timeseries removed)
31!
32! 2249 2017-06-06 13:58:01Z sward
33! Allow for RRTMG runs without humidity/cloud physics
34!
35! 2248 2017-06-06 13:52:54Z sward
36! Error no changed
37!
38! 2233 2017-05-30 18:08:54Z suehring
39!
40! 2232 2017-05-30 17:47:52Z suehring
41! Adjustments to new topography concept
42! Bugfix in read restart
43!
44! 2200 2017-04-11 11:37:51Z suehring
45! Bugfix in call of exchange_horiz_2d and read restart data
46!
47! 2163 2017-03-01 13:23:15Z schwenkel
48! Bugfix in radiation_check_data_output
49!
50! 2157 2017-02-22 15:10:35Z suehring
51! Bugfix in read_restart data
52!
53! 2011 2016-09-19 17:29:57Z kanani
54! Removed CALL of auxiliary SUBROUTINE get_usm_info,
55! flag urban_surface is now defined in module control_parameters.
56!
57! 2007 2016-08-24 15:47:17Z kanani
58! Added calculation of solar directional vector for new urban surface
59! model,
60! accounted for urban_surface model in radiation_check_parameters,
61! correction of comments for zenith angle.
62!
63! 2000 2016-08-20 18:09:15Z knoop
64! Forced header and separation lines into 80 columns
65!
66! 1976 2016-07-27 13:28:04Z maronga
67! Output of 2D/3D/masked data is now directly done within this module. The
68! radiation schemes have been simplified for better usability so that
69! rad_lw_in, rad_lw_out, rad_sw_in, and rad_sw_out are available independent of
70! the radiation code used.
71!
72! 1856 2016-04-13 12:56:17Z maronga
73! Bugfix: allocation of rad_lw_out for radiation_scheme = 'clear-sky'
74!
75! 1853 2016-04-11 09:00:35Z maronga
76! Added routine for radiation_scheme = constant.
77
78! 1849 2016-04-08 11:33:18Z hoffmann
79! Adapted for modularization of microphysics
80!
81! 1826 2016-04-07 12:01:39Z maronga
82! Further modularization.
83!
84! 1788 2016-03-10 11:01:04Z maronga
85! Added new albedo class for pavements / roads.
86!
87! 1783 2016-03-06 18:36:17Z raasch
88! palm-netcdf-module removed in order to avoid a circular module dependency,
89! netcdf-variables moved to netcdf-module, new routine netcdf_handle_error_rad
90! added
91!
92! 1757 2016-02-22 15:49:32Z maronga
93! Added parameter unscheduled_radiation_calls. Bugfix: interpolation of sounding
94! profiles for pressure and temperature above the LES domain.
95!
96! 1709 2015-11-04 14:47:01Z maronga
97! Bugfix: set initial value for rrtm_lwuflx_dt to zero, small formatting
98! corrections
99!
100! 1701 2015-11-02 07:43:04Z maronga
101! Bugfixes: wrong index for output of timeseries, setting of nz_snd_end
102!
103! 1691 2015-10-26 16:17:44Z maronga
104! Added option for spin-up runs without radiation (skip_time_do_radiation). Bugfix
105! in calculation of pressure profiles. Bugfix in calculation of trace gas profiles.
106! Added output of radiative heating rates.
107!
108! 1682 2015-10-07 23:56:08Z knoop
109! Code annotations made doxygen readable
110!
111! 1606 2015-06-29 10:43:37Z maronga
112! Added preprocessor directive __netcdf to allow for compiling without netCDF.
113! Note, however, that RRTMG cannot be used without netCDF.
114!
115! 1590 2015-05-08 13:56:27Z maronga
116! Bugfix: definition of character strings requires same length for all elements
117!
118! 1587 2015-05-04 14:19:01Z maronga
119! Added albedo class for snow
120!
121! 1585 2015-04-30 07:05:52Z maronga
122! Added support for RRTMG
123!
124! 1571 2015-03-12 16:12:49Z maronga
125! Added missing KIND attribute. Removed upper-case variable names
126!
127! 1551 2015-03-03 14:18:16Z maronga
128! Added support for data output. Various variables have been renamed. Added
129! interface for different radiation schemes (currently: clear-sky, constant, and
130! RRTM (not yet implemented).
131!
132! 1496 2014-12-02 17:25:50Z maronga
133! Initial revision
134!
135!
136! Description:
137! ------------
138!> Radiation models and interfaces
139!> @todo move variable definitions used in radiation_init only to the subroutine
140!>       as they are no longer required after initialization.
141!> @todo Output of full column vertical profiles used in RRTMG
142!> @todo Output of other rrtm arrays (such as volume mixing ratios)
143!> @todo Adapt for use with topography
144!>
145!> @note Many variables have a leading dummy dimension (0:0) in order to
146!>       match the assume-size shape expected by the RRTMG model.
147!------------------------------------------------------------------------------!
148 MODULE radiation_model_mod
149 
150    USE arrays_3d,                                                             &
151        ONLY:  dzw, hyp, pt, q, ql, zu, zw
152
153    USE cloud_parameters,                                                      &
154        ONLY:  cp, l_d_cp, rho_l
155
156    USE constants,                                                             &
157        ONLY:  pi
158
159    USE control_parameters,                                                    &
160        ONLY:  cloud_droplets, cloud_physics, g, initializing_actions,         &
161               large_scale_forcing, lsf_surf, phi, pt_surface, rho_surface,    &
162               surface_pressure, time_since_reference_point
163
164    USE indices,                                                               &
165        ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt,           &
166               wall_flags_0
167
168    USE kinds
169
170    USE microphysics_mod,                                                      &
171        ONLY:  nc_const, sigma_gc
172
173#if defined ( __netcdf )
174    USE NETCDF
175#endif
176
177#if defined ( __rrtmg )
178    USE parrrsw,                                                               &
179        ONLY:  naerec, nbndsw
180
181    USE parrrtm,                                                               &
182        ONLY:  nbndlw
183
184    USE rrtmg_lw_init,                                                         &
185        ONLY:  rrtmg_lw_ini
186
187    USE rrtmg_sw_init,                                                         &
188        ONLY:  rrtmg_sw_ini
189
190    USE rrtmg_lw_rad,                                                          &
191        ONLY:  rrtmg_lw
192
193    USE rrtmg_sw_rad,                                                          &
194        ONLY:  rrtmg_sw
195#endif
196
197
198
199    IMPLICIT NONE
200
201    CHARACTER(10) :: radiation_scheme = 'clear-sky' ! 'constant', 'clear-sky', or 'rrtmg'
202
203!
204!-- Predefined Land surface classes (albedo_type) after Briegleb (1992)
205    CHARACTER(37), DIMENSION(0:17), PARAMETER :: albedo_type_name = (/      &
206                                   'user defined                         ', & !  0
207                                   'ocean                                ', & !  1
208                                   'mixed farming, tall grassland        ', & !  2
209                                   'tall/medium grassland                ', & !  3
210                                   'evergreen shrubland                  ', & !  4
211                                   'short grassland/meadow/shrubland     ', & !  5
212                                   'evergreen needleleaf forest          ', & !  6
213                                   'mixed deciduous evergreen forest     ', & !  7
214                                   'deciduous forest                     ', & !  8
215                                   'tropical evergreen broadleaved forest', & !  9
216                                   'medium/tall grassland/woodland       ', & ! 10
217                                   'desert, sandy                        ', & ! 11
218                                   'desert, rocky                        ', & ! 12
219                                   'tundra                               ', & ! 13
220                                   'land ice                             ', & ! 14
221                                   'sea ice                              ', & ! 15
222                                   'snow                                 ', & ! 16
223                                   'pavement/roads                       '  & ! 17
224                                                         /)
225
226    INTEGER(iwp) :: albedo_type  = 5,    & !< Albedo surface type (default: short grassland)
227                    day,                 & !< current day of the year
228                    day_init     = 172,  & !< day of the year at model start (21/06)
229                    dots_rad     = 0       !< starting index for timeseries output
230
231    LOGICAL ::  unscheduled_radiation_calls = .TRUE., & !< flag parameter indicating whether additional calls of the radiation code are allowed
232                constant_albedo = .FALSE.,            & !< flag parameter indicating whether the albedo may change depending on zenith
233                force_radiation_call = .FALSE.,       & !< flag parameter for unscheduled radiation calls
234                lw_radiation = .TRUE.,                & !< flag parameter indicating whether longwave radiation shall be calculated
235                radiation = .FALSE.,                  & !< flag parameter indicating whether the radiation model is used
236                sun_up    = .TRUE.,                   & !< flag parameter indicating whether the sun is up or down
237                sw_radiation = .TRUE.,                 & !< flag parameter indicing whether shortwave radiation shall be calculated
238                sun_direction = .FALSE.                 !< flag parameter indicing whether solar direction shall be calculated
239
240
241    REAL(wp), PARAMETER :: d_seconds_hour  = 0.000277777777778_wp,  & !< inverse of seconds per hour (1/3600)
242                           d_hours_day    = 0.0416666666667_wp,     & !< inverse of hours per day (1/24)
243                           sigma_sb       = 5.67037321E-8_wp,       & !< Stefan-Boltzmann constant
244                           solar_constant = 1368.0_wp                 !< solar constant at top of atmosphere
245
246    REAL(wp) :: albedo = 9999999.9_wp,           & !< NAMELIST alpha
247                albedo_lw_dif = 9999999.9_wp,    & !< NAMELIST aldif
248                albedo_lw_dir = 9999999.9_wp,    & !< NAMELIST aldir
249                albedo_sw_dif = 9999999.9_wp,    & !< NAMELIST asdif
250                albedo_sw_dir = 9999999.9_wp,    & !< NAMELIST asdir
251                decl_1,                          & !< declination coef. 1
252                decl_2,                          & !< declination coef. 2
253                decl_3,                          & !< declination coef. 3
254                dt_radiation = 0.0_wp,           & !< radiation model timestep
255                emissivity = 0.98_wp,            & !< NAMELIST surface emissivity
256                lambda = 0.0_wp,                 & !< longitude in degrees
257                lon = 0.0_wp,                    & !< longitude in radians
258                lat = 0.0_wp,                    & !< latitude in radians
259                net_radiation = 0.0_wp,          & !< net radiation at surface
260                skip_time_do_radiation = 0.0_wp, & !< Radiation model is not called before this time
261                sky_trans,                       & !< sky transmissivity
262                time_radiation = 0.0_wp,         & !< time since last call of radiation code
263                time_utc,                        & !< current time in UTC
264                time_utc_init = 43200.0_wp         !< UTC time at model start (noon)
265
266    REAL(wp), DIMENSION(0:0) ::  zenith,         & !< cosine of solar zenith angle
267                                 sun_dir_lat,    & !< solar directional vector in latitudes
268                                 sun_dir_lon       !< solar directional vector in longitudes
269
270    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
271                alpha,                       & !< surface broadband albedo (used for clear-sky scheme)
272                rad_lw_out_change_0,         & !< change in LW out due to change in surface temperature
273                rad_net,                     & !< net radiation at the surface
274                rad_net_av                     !< average of rad_net
275
276!
277!-- Land surface albedos for solar zenith angle of 60° after Briegleb (1992)     
278!-- (shortwave, longwave, broadband):   sw,      lw,      bb,
279    REAL(wp), DIMENSION(0:2,1:17), PARAMETER :: albedo_pars = RESHAPE( (/& 
280                                   0.06_wp, 0.06_wp, 0.06_wp,            & !  1
281                                   0.09_wp, 0.28_wp, 0.19_wp,            & !  2
282                                   0.11_wp, 0.33_wp, 0.23_wp,            & !  3
283                                   0.11_wp, 0.33_wp, 0.23_wp,            & !  4
284                                   0.14_wp, 0.34_wp, 0.25_wp,            & !  5
285                                   0.06_wp, 0.22_wp, 0.14_wp,            & !  6
286                                   0.06_wp, 0.27_wp, 0.17_wp,            & !  7
287                                   0.06_wp, 0.31_wp, 0.19_wp,            & !  8
288                                   0.06_wp, 0.22_wp, 0.14_wp,            & !  9
289                                   0.06_wp, 0.28_wp, 0.18_wp,            & ! 10
290                                   0.35_wp, 0.51_wp, 0.43_wp,            & ! 11
291                                   0.24_wp, 0.40_wp, 0.32_wp,            & ! 12
292                                   0.10_wp, 0.27_wp, 0.19_wp,            & ! 13
293                                   0.90_wp, 0.65_wp, 0.77_wp,            & ! 14
294                                   0.90_wp, 0.65_wp, 0.77_wp,            & ! 15
295                                   0.95_wp, 0.70_wp, 0.82_wp,            & ! 16
296                                   0.08_wp, 0.08_wp, 0.08_wp             & ! 17
297                                 /), (/ 3, 17 /) )
298
299    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: &
300                        rad_lw_cs_hr,                  & !< longwave clear sky radiation heating rate (K/s)
301                        rad_lw_cs_hr_av,               & !< average of rad_lw_cs_hr
302                        rad_lw_hr,                     & !< longwave radiation heating rate (K/s)
303                        rad_lw_hr_av,                  & !< average of rad_sw_hr
304                        rad_lw_in,                     & !< incoming longwave radiation (W/m2)
305                        rad_lw_in_av,                  & !< average of rad_lw_in
306                        rad_lw_out,                    & !< outgoing longwave radiation (W/m2)
307                        rad_lw_out_av,                 & !< average of rad_lw_out
308                        rad_sw_cs_hr,                  & !< shortwave clear sky radiation heating rate (K/s)
309                        rad_sw_cs_hr_av,               & !< average of rad_sw_cs_hr
310                        rad_sw_hr,                     & !< shortwave radiation heating rate (K/s)
311                        rad_sw_hr_av,                  & !< average of rad_sw_hr
312                        rad_sw_in,                     & !< incoming shortwave radiation (W/m2)
313                        rad_sw_in_av,                  & !< average of rad_sw_in
314                        rad_sw_out,                    & !< outgoing shortwave radiation (W/m2)
315                        rad_sw_out_av                    !< average of rad_sw_out
316
317
318!
319!-- Variables and parameters used in RRTMG only
320#if defined ( __rrtmg )
321    CHARACTER(LEN=12) :: rrtm_input_file = "RAD_SND_DATA" !< name of the NetCDF input file (sounding data)
322
323
324!
325!-- Flag parameters for RRTMGS (should not be changed)
326    INTEGER(iwp), PARAMETER :: rrtm_idrv     = 1, & !< flag for longwave upward flux calculation option (0,1)
327                               rrtm_inflglw  = 2, & !< flag for lw cloud optical properties (0,1,2)
328                               rrtm_iceflglw = 0, & !< flag for lw ice particle specifications (0,1,2,3)
329                               rrtm_liqflglw = 1, & !< flag for lw liquid droplet specifications
330                               rrtm_inflgsw  = 2, & !< flag for sw cloud optical properties (0,1,2)
331                               rrtm_iceflgsw = 0, & !< flag for sw ice particle specifications (0,1,2,3)
332                               rrtm_liqflgsw = 1    !< flag for sw liquid droplet specifications
333
334!
335!-- The following variables should be only changed with care, as this will
336!-- require further setting of some variables, which is currently not
337!-- implemented (aerosols, ice phase).
338    INTEGER(iwp) :: nzt_rad,           & !< upper vertical limit for radiation calculations
339                    rrtm_icld = 0,     & !< cloud flag (0: clear sky column, 1: cloudy column)
340                    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)
341
342    INTEGER(iwp) :: nc_stat !< local variable for storin the result of netCDF calls for error message handling
343
344    LOGICAL :: snd_exists = .FALSE.      !< flag parameter to check whether a user-defined input files exists
345
346    REAL(wp), PARAMETER :: mol_mass_air_d_wv = 1.607793_wp !< molecular weight dry air / water vapor
347
348    REAL(wp), DIMENSION(:), ALLOCATABLE :: hyp_snd,     & !< hypostatic pressure from sounding data (hPa)
349                                           q_snd,       & !< specific humidity from sounding data (kg/kg) - dummy at the moment
350                                           rrtm_tsfc,   & !< dummy array for storing surface temperature
351                                           t_snd          !< actual temperature from sounding data (hPa)
352
353    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: aldif,          & !< longwave diffuse albedo solar angle of 60°
354                                             aldir,          & !< longwave direct albedo solar angle of 60°
355                                             asdif,          & !< shortwave diffuse albedo solar angle of 60°
356                                             asdir,          & !< shortwave direct albedo solar angle of 60°
357                                             rrtm_ccl4vmr,   & !< CCL4 volume mixing ratio (g/mol)
358                                             rrtm_cfc11vmr,  & !< CFC11 volume mixing ratio (g/mol)
359                                             rrtm_cfc12vmr,  & !< CFC12 volume mixing ratio (g/mol)
360                                             rrtm_cfc22vmr,  & !< CFC22 volume mixing ratio (g/mol)
361                                             rrtm_ch4vmr,    & !< CH4 volume mixing ratio
362                                             rrtm_cicewp,    & !< in-cloud ice water path (g/m²)
363                                             rrtm_cldfr,     & !< cloud fraction (0,1)
364                                             rrtm_cliqwp,    & !< in-cloud liquid water path (g/m²)
365                                             rrtm_co2vmr,    & !< CO2 volume mixing ratio (g/mol)
366                                             rrtm_emis,      & !< surface emissivity (0-1)   
367                                             rrtm_h2ovmr,    & !< H2O volume mixing ratio
368                                             rrtm_n2ovmr,    & !< N2O volume mixing ratio
369                                             rrtm_o2vmr,     & !< O2 volume mixing ratio
370                                             rrtm_o3vmr,     & !< O3 volume mixing ratio
371                                             rrtm_play,      & !< pressure layers (hPa, zu-grid)
372                                             rrtm_plev,      & !< pressure layers (hPa, zw-grid)
373                                             rrtm_reice,     & !< cloud ice effective radius (microns)
374                                             rrtm_reliq,     & !< cloud water drop effective radius (microns)
375                                             rrtm_tlay,      & !< actual temperature (K, zu-grid)
376                                             rrtm_tlev,      & !< actual temperature (K, zw-grid)
377                                             rrtm_lwdflx,    & !< RRTM output of incoming longwave radiation flux (W/m2)
378                                             rrtm_lwdflxc,   & !< RRTM output of outgoing clear sky longwave radiation flux (W/m2)
379                                             rrtm_lwuflx,    & !< RRTM output of outgoing longwave radiation flux (W/m2)
380                                             rrtm_lwuflxc,   & !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
381                                             rrtm_lwuflx_dt, & !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
382                                             rrtm_lwuflxc_dt,& !< RRTM output of outgoing clear sky longwave radiation flux (W/m2)
383                                             rrtm_lwhr,      & !< RRTM output of longwave radiation heating rate (K/d)
384                                             rrtm_lwhrc,     & !< RRTM output of incoming longwave clear sky radiation heating rate (K/d)
385                                             rrtm_swdflx,    & !< RRTM output of incoming shortwave radiation flux (W/m2)
386                                             rrtm_swdflxc,   & !< RRTM output of outgoing clear sky shortwave radiation flux (W/m2)
387                                             rrtm_swuflx,    & !< RRTM output of outgoing shortwave radiation flux (W/m2)
388                                             rrtm_swuflxc,   & !< RRTM output of incoming clear sky shortwave radiation flux (W/m2)
389                                             rrtm_swhr,      & !< RRTM output of shortwave radiation heating rate (K/d)
390                                             rrtm_swhrc        !< RRTM output of incoming shortwave clear sky radiation heating rate (K/d)
391
392!
393!-- Definition of arrays that are currently not used for calling RRTMG (due to setting of flag parameters)
394    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  rad_lw_cs_in,   & !< incoming clear sky longwave radiation (W/m2) (not used)
395                                                rad_lw_cs_out,  & !< outgoing clear sky longwave radiation (W/m2) (not used)
396                                                rad_sw_cs_in,   & !< incoming clear sky shortwave radiation (W/m2) (not used)
397                                                rad_sw_cs_out,  & !< outgoing clear sky shortwave radiation (W/m2) (not used)
398                                                rrtm_aldif,     & !< surface albedo for longwave diffuse radiation
399                                                rrtm_aldir,     & !< surface albedo for longwave direct radiation
400                                                rrtm_asdif,     & !< surface albedo for shortwave diffuse radiation
401                                                rrtm_asdir,     & !< surface albedo for shortwave direct radiation
402                                                rrtm_lw_tauaer, & !< lw aerosol optical depth
403                                                rrtm_lw_taucld, & !< lw in-cloud optical depth
404                                                rrtm_sw_taucld, & !< sw in-cloud optical depth
405                                                rrtm_sw_ssacld, & !< sw in-cloud single scattering albedo
406                                                rrtm_sw_asmcld, & !< sw in-cloud asymmetry parameter
407                                                rrtm_sw_fsfcld, & !< sw in-cloud forward scattering fraction
408                                                rrtm_sw_tauaer, & !< sw aerosol optical depth
409                                                rrtm_sw_ssaaer, & !< sw aerosol single scattering albedo
410                                                rrtm_sw_asmaer, & !< sw aerosol asymmetry parameter
411                                                rrtm_sw_ecaer     !< sw aerosol optical detph at 0.55 microns (rrtm_iaer = 6 only)
412
413#endif
414
415    INTERFACE radiation_check_data_output
416       MODULE PROCEDURE radiation_check_data_output
417    END INTERFACE radiation_check_data_output
418
419    INTERFACE radiation_check_data_output_pr
420       MODULE PROCEDURE radiation_check_data_output_pr
421    END INTERFACE radiation_check_data_output_pr
422 
423    INTERFACE radiation_check_parameters
424       MODULE PROCEDURE radiation_check_parameters
425    END INTERFACE radiation_check_parameters
426 
427    INTERFACE radiation_clearsky
428       MODULE PROCEDURE radiation_clearsky
429    END INTERFACE radiation_clearsky
430 
431    INTERFACE radiation_constant
432       MODULE PROCEDURE radiation_constant
433    END INTERFACE radiation_constant
434 
435    INTERFACE radiation_control
436       MODULE PROCEDURE radiation_control
437    END INTERFACE radiation_control
438
439    INTERFACE radiation_3d_data_averaging
440       MODULE PROCEDURE radiation_3d_data_averaging
441    END INTERFACE radiation_3d_data_averaging
442
443    INTERFACE radiation_data_output_2d
444       MODULE PROCEDURE radiation_data_output_2d
445    END INTERFACE radiation_data_output_2d
446
447    INTERFACE radiation_data_output_3d
448       MODULE PROCEDURE radiation_data_output_3d
449    END INTERFACE radiation_data_output_3d
450
451    INTERFACE radiation_data_output_mask
452       MODULE PROCEDURE radiation_data_output_mask
453    END INTERFACE radiation_data_output_mask
454
455    INTERFACE radiation_define_netcdf_grid
456       MODULE PROCEDURE radiation_define_netcdf_grid
457    END INTERFACE radiation_define_netcdf_grid
458
459    INTERFACE radiation_header
460       MODULE PROCEDURE radiation_header
461    END INTERFACE radiation_header 
462 
463    INTERFACE radiation_init
464       MODULE PROCEDURE radiation_init
465    END INTERFACE radiation_init
466
467    INTERFACE radiation_parin
468       MODULE PROCEDURE radiation_parin
469    END INTERFACE radiation_parin
470   
471    INTERFACE radiation_rrtmg
472       MODULE PROCEDURE radiation_rrtmg
473    END INTERFACE radiation_rrtmg
474
475    INTERFACE radiation_tendency
476       MODULE PROCEDURE radiation_tendency
477       MODULE PROCEDURE radiation_tendency_ij
478    END INTERFACE radiation_tendency
479
480    INTERFACE radiation_read_restart_data
481       MODULE PROCEDURE radiation_read_restart_data
482    END INTERFACE radiation_read_restart_data
483
484    INTERFACE radiation_last_actions
485       MODULE PROCEDURE radiation_last_actions
486    END INTERFACE radiation_last_actions
487
488    SAVE
489
490    PRIVATE
491
492!
493!-- Public functions / NEEDS SORTING
494    PUBLIC radiation_check_data_output, radiation_check_data_output_pr,        &
495           radiation_check_parameters, radiation_control,                      &
496           radiation_header, radiation_init, radiation_parin,                  &
497           radiation_3d_data_averaging, radiation_tendency,                    &
498           radiation_data_output_2d, radiation_data_output_3d,                 &
499           radiation_define_netcdf_grid, radiation_last_actions,               &
500           radiation_read_restart_data, radiation_data_output_mask
501   
502!
503!-- Public variables and constants / NEEDS SORTING
504    PUBLIC dots_rad, dt_radiation, force_radiation_call,                       &
505           rad_net, rad_net_av, radiation, radiation_scheme, rad_lw_in,        &
506           rad_lw_in_av, rad_lw_out, rad_lw_out_av, rad_lw_out_change_0,       &
507           rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av, rad_sw_in,  &
508           rad_sw_in_av, rad_sw_out, rad_sw_out_av, rad_sw_cs_hr,              &
509           rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av, sigma_sb,                 &
510           skip_time_do_radiation, time_radiation, unscheduled_radiation_calls,&
511           zenith, calc_zenith, sun_direction, sun_dir_lat, sun_dir_lon,       &
512           day_init, time_utc_init
513
514
515#if defined ( __rrtmg )
516    PUBLIC rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir
517#endif
518
519 CONTAINS
520
521
522!------------------------------------------------------------------------------!
523! Description:
524! ------------
525!> This subroutine controls the calls of the radiation schemes
526!------------------------------------------------------------------------------!
527    SUBROUTINE radiation_control
528 
529 
530       IMPLICIT NONE
531
532
533       SELECT CASE ( TRIM( radiation_scheme ) )
534
535          CASE ( 'constant' )
536             CALL radiation_constant
537         
538          CASE ( 'clear-sky' )
539             CALL radiation_clearsky
540       
541          CASE ( 'rrtmg' )
542             CALL radiation_rrtmg
543
544          CASE DEFAULT
545
546       END SELECT
547
548
549    END SUBROUTINE radiation_control
550
551!------------------------------------------------------------------------------!
552! Description:
553! ------------
554!> Check data output for radiation model
555!------------------------------------------------------------------------------!
556    SUBROUTINE radiation_check_data_output( var, unit, i, ilen, k )
557 
558 
559       USE control_parameters,                                                 &
560           ONLY: data_output, message_string
561
562       IMPLICIT NONE
563
564       CHARACTER (LEN=*) ::  unit     !<
565       CHARACTER (LEN=*) ::  var      !<
566
567       INTEGER(iwp) :: i
568       INTEGER(iwp) :: ilen
569       INTEGER(iwp) :: k
570
571       SELECT CASE ( TRIM( var ) )
572
573         CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_sw_cs_hr', 'rad_sw_hr' )
574             IF (  .NOT.  radiation  .OR.  radiation_scheme /= 'rrtmg' )  THEN
575                message_string = '"output of "' // TRIM( var ) // '" requi' // &
576                                 'res radiation = .TRUE. and ' //              &
577                                 'radiation_scheme = "rrtmg"'
578                CALL message( 'check_parameters', 'PA0406', 1, 2, 0, 6, 0 )
579             ENDIF
580             unit = 'K/h'     
581
582         CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_sw_in', 'rad_sw_out' )
583             IF (  .NOT.  radiation  .OR.  radiation_scheme /= 'rrtmg' )  THEN
584                message_string = '"output of "' // TRIM( var ) // '" requi' // &
585                                 'res radiation = .TRUE. and ' //              &
586                                 'radiation_scheme = "rrtmg"'
587                CALL message( 'check_parameters', 'PA0406', 1, 2, 0, 6, 0 )
588             ENDIF
589             unit = 'W/m2'   
590
591          CASE ( 'rad_net*', 'rrtm_aldif*', 'rrtm_aldir*', 'rrtm_asdif*',      &
592                 'rrtm_asdir*' )
593             IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
594                message_string = 'illegal value for data_output: "' //         &
595                                 TRIM( var ) // '" & only 2d-horizontal ' //   &
596                                 'cross sections are allowed for this value'
597                CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
598             ENDIF
599             IF (  .NOT.  radiation  .OR.  radiation_scheme /= "rrtmg" )  THEN
600                IF ( TRIM( var ) == 'rrtm_aldif*'  .OR.                        &
601                     TRIM( var ) == 'rrtm_aldir*'  .OR.                        &
602                     TRIM( var ) == 'rrtm_asdif*'  .OR.                        &
603                     TRIM( var ) == 'rrtm_asdir*'      )                       &
604                THEN
605                   message_string = 'output of "' // TRIM( var ) // '" require'&
606                                    // 's radiation = .TRUE. and radiation_sch'&
607                                    // 'eme = "rrtmg"'
608                   CALL message( 'check_parameters', 'PA0409', 1, 2, 0, 6, 0 )
609                ENDIF
610             ENDIF
611
612             IF ( TRIM( var ) == 'rad_net*'      ) unit = 'W/m2'   
613             IF ( TRIM( var ) == 'rrtm_aldif*'   ) unit = ''   
614             IF ( TRIM( var ) == 'rrtm_aldir*'   ) unit = '' 
615             IF ( TRIM( var ) == 'rrtm_asdif*'   ) unit = '' 
616             IF ( TRIM( var ) == 'rrtm_asdir*'   ) unit = '' 
617
618          CASE DEFAULT
619             unit = 'illegal'
620
621       END SELECT
622
623
624    END SUBROUTINE radiation_check_data_output
625
626!------------------------------------------------------------------------------!
627! Description:
628! ------------
629!> Check data output of profiles for radiation model
630!------------------------------------------------------------------------------! 
631    SUBROUTINE radiation_check_data_output_pr( variable, var_count, unit, dopr_unit )
632 
633       USE arrays_3d,                                                          &
634           ONLY: zu
635
636       USE control_parameters,                                                 &
637           ONLY: data_output_pr, message_string
638
639       USE indices
640
641       USE profil_parameter
642
643       USE statistics
644
645       IMPLICIT NONE
646   
647       CHARACTER (LEN=*) ::  unit      !<
648       CHARACTER (LEN=*) ::  variable  !<
649       CHARACTER (LEN=*) ::  dopr_unit !< local value of dopr_unit
650 
651       INTEGER(iwp) ::  user_pr_index !<
652       INTEGER(iwp) ::  var_count     !<
653
654       SELECT CASE ( TRIM( variable ) )
655       
656         CASE ( 'rad_net' )
657             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme == 'constant' )&
658             THEN
659                message_string = 'data_output_pr = ' //                        &
660                                 TRIM( data_output_pr(var_count) ) // ' is' // &
661                                 'not available for radiation = .FALSE. or ' //&
662                                 'radiation_scheme = "constant"'
663                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
664             ELSE
665                dopr_index(var_count) = 99
666                dopr_unit  = 'W/m2'
667                hom(:,2,99,:)  = SPREAD( zw, 2, statistic_regions+1 )
668                unit = dopr_unit
669             ENDIF
670
671          CASE ( 'rad_lw_in' )
672             IF ( (  .NOT.  radiation)  .OR.  radiation_scheme == 'constant' ) &
673             THEN
674                message_string = 'data_output_pr = ' //                        &
675                                 TRIM( data_output_pr(var_count) ) // ' is' // &
676                                 'not available for radiation = .FALSE. or ' //&
677                                 'radiation_scheme = "constant"'
678                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
679             ELSE
680                dopr_index(var_count) = 100
681                dopr_unit  = 'W/m2'
682                hom(:,2,100,:)  = SPREAD( zw, 2, statistic_regions+1 )
683                unit = dopr_unit 
684             ENDIF
685
686          CASE ( 'rad_lw_out' )
687             IF ( (  .NOT. radiation )  .OR.  radiation_scheme == 'constant' ) &
688             THEN
689                message_string = 'data_output_pr = ' //                        &
690                                 TRIM( data_output_pr(var_count) ) // ' is' // &
691                                 'not available for radiation = .FALSE. or ' //&
692                                 'radiation_scheme = "constant"'
693                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
694             ELSE
695                dopr_index(var_count) = 101
696                dopr_unit  = 'W/m2'
697                hom(:,2,101,:)  = SPREAD( zw, 2, statistic_regions+1 )
698                unit = dopr_unit   
699             ENDIF
700
701          CASE ( 'rad_sw_in' )
702             IF ( (  .NOT. radiation )  .OR.  radiation_scheme == 'constant' ) &
703             THEN
704                message_string = 'data_output_pr = ' //                        &
705                                 TRIM( data_output_pr(var_count) ) // ' is' // &
706                                 'not available for radiation = .FALSE. or ' //&
707                                 'radiation_scheme = "constant"'
708                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
709             ELSE
710                dopr_index(var_count) = 102
711                dopr_unit  = 'W/m2'
712                hom(:,2,102,:)  = SPREAD( zw, 2, statistic_regions+1 )
713                unit = dopr_unit
714             ENDIF
715
716          CASE ( 'rad_sw_out')
717             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme == 'constant' )&
718             THEN
719                message_string = 'data_output_pr = ' //                        &
720                                 TRIM( data_output_pr(var_count) ) // ' is' // &
721                                 'not available for radiation = .FALSE. or ' //&
722                                 'radiation_scheme = "constant"'
723                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
724             ELSE
725                dopr_index(var_count) = 103
726                dopr_unit  = 'W/m2'
727                hom(:,2,103,:)  = SPREAD( zw, 2, statistic_regions+1 )
728                unit = dopr_unit
729             ENDIF
730
731          CASE ( 'rad_lw_cs_hr' )
732             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
733             THEN
734                message_string = 'data_output_pr = ' //                        &
735                                 TRIM( data_output_pr(var_count) ) // ' is' // &
736                                 'not available for radiation = .FALSE. or ' //&
737                                 'radiation_scheme /= "rrtmg"'
738                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
739             ELSE
740                dopr_index(var_count) = 104
741                dopr_unit  = 'K/h'
742                hom(:,2,104,:)  = SPREAD( zu, 2, statistic_regions+1 )
743                unit = dopr_unit
744             ENDIF
745
746          CASE ( 'rad_lw_hr' )
747             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
748             THEN
749                message_string = 'data_output_pr = ' //                        &
750                                 TRIM( data_output_pr(var_count) ) // ' is' // &
751                                 'not available for radiation = .FALSE. or ' //&
752                                 'radiation_scheme /= "rrtmg"'
753                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
754             ELSE
755                dopr_index(var_count) = 105
756                dopr_unit  = 'K/h'
757                hom(:,2,105,:)  = SPREAD( zu, 2, statistic_regions+1 )
758                unit = dopr_unit
759             ENDIF
760
761          CASE ( 'rad_sw_cs_hr' )
762             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
763             THEN
764                message_string = 'data_output_pr = ' //                        &
765                                 TRIM( data_output_pr(var_count) ) // ' is' // &
766                                 'not available for radiation = .FALSE. or ' //&
767                                 'radiation_scheme /= "rrtmg"'
768                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
769             ELSE
770                dopr_index(var_count) = 106
771                dopr_unit  = 'K/h'
772                hom(:,2,106,:)  = SPREAD( zu, 2, statistic_regions+1 )
773                unit = dopr_unit
774             ENDIF
775
776          CASE ( 'rad_sw_hr' )
777             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
778             THEN
779                message_string = 'data_output_pr = ' //                        &
780                                 TRIM( data_output_pr(var_count) ) // ' is' // &
781                                 'not available for radiation = .FALSE. or ' //&
782                                 'radiation_scheme /= "rrtmg"'
783                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
784             ELSE
785                dopr_index(var_count) = 107
786                dopr_unit  = 'K/h'
787                hom(:,2,107,:)  = SPREAD( zu, 2, statistic_regions+1 )
788                unit = dopr_unit
789             ENDIF
790
791
792          CASE DEFAULT
793             unit = 'illegal'
794
795       END SELECT
796
797
798    END SUBROUTINE radiation_check_data_output_pr
799 
800 
801!------------------------------------------------------------------------------!
802! Description:
803! ------------
804!> Check parameters routine for radiation model
805!------------------------------------------------------------------------------!
806    SUBROUTINE radiation_check_parameters
807
808       USE control_parameters,                                                 &
809           ONLY: message_string, topography, urban_surface
810                 
811   
812       IMPLICIT NONE
813       
814
815       IF ( radiation_scheme /= 'constant'   .AND.                             &
816            radiation_scheme /= 'clear-sky'  .AND.                             &
817            radiation_scheme /= 'rrtmg' )  THEN
818          message_string = 'unknown radiation_scheme = '//                     &
819                           TRIM( radiation_scheme )
820          CALL message( 'check_parameters', 'PA0405', 1, 2, 0, 6, 0 )
821       ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
822#if ! defined ( __rrtmg )
823          message_string = 'radiation_scheme = "rrtmg" requires ' //           & 
824                           'compilation of PALM with pre-processor ' //        &
825                           'directive -D__rrtmg'
826          CALL message( 'check_parameters', 'PA0407', 1, 2, 0, 6, 0 )
827#endif
828#if defined ( __rrtmg ) && ! defined( __netcdf )
829          message_string = 'radiation_scheme = "rrtmg" requires ' //           & 
830                           'the use of NetCDF (preprocessor directive ' //     &
831                           '-D__netcdf'
832          CALL message( 'check_parameters', 'PA0412', 1, 2, 0, 6, 0 )
833#endif
834
835       ENDIF
836
837       IF ( albedo_type == 0  .AND.  albedo == 9999999.9_wp  .AND.             &
838            radiation_scheme == 'clear-sky')  THEN
839          message_string = 'radiation_scheme = "clear-sky" in combination' //  & 
840                           'with albedo_type = 0 requires setting of albedo'// &
841                           ' /= 9999999.9'
842          CALL message( 'check_parameters', 'PA0410', 1, 2, 0, 6, 0 )
843       ENDIF
844
845       IF ( albedo_type == 0  .AND.  radiation_scheme == 'rrtmg'  .AND.        &
846          (    albedo_lw_dif == 9999999.9_wp .OR. albedo_lw_dir == 9999999.9_wp&
847          .OR. albedo_sw_dif == 9999999.9_wp .OR. albedo_sw_dir == 9999999.9_wp& 
848          ) ) THEN
849          message_string = 'radiation_scheme = "rrtmg" in combination' //      & 
850                           'with albedo_type = 0 requires setting of ' //      &
851                           'albedo_lw_dif /= 9999999.9' //                     &
852                           'albedo_lw_dir /= 9999999.9' //                     &
853                           'albedo_sw_dif /= 9999999.9 and' //                 &
854                           'albedo_sw_dir /= 9999999.9'
855          CALL message( 'check_parameters', 'PA0411', 1, 2, 0, 6, 0 )
856       ENDIF
857
858!
859!--    The following paramter check is temporarily extended by the urban_surface
860!--    flag, until a better solution comes up to omit this check in case of
861!--    urban surface model is used.
862       IF ( topography /= 'flat'  .AND.  .NOT.  urban_surface )  THEN
863          message_string = 'radiation scheme cannot be used ' //               & 
864                           'in combination with  topography /= "flat"'
865          CALL message( 'check_parameters', 'PA0414', 1, 2, 0, 6, 0 )
866       ENDIF
867 
868    END SUBROUTINE radiation_check_parameters 
869 
870 
871!------------------------------------------------------------------------------!
872! Description:
873! ------------
874!> Initialization of the radiation model
875!------------------------------------------------------------------------------!
876    SUBROUTINE radiation_init
877   
878       IMPLICIT NONE
879
880!
881!--    Allocate array for storing the surface net radiation
882       IF ( .NOT. ALLOCATED ( rad_net ) )  THEN
883          ALLOCATE ( rad_net(nysg:nyng,nxlg:nxrg) )
884          rad_net = 0.0_wp
885       ENDIF
886
887!
888!--    Allocate array for storing the surface net radiation
889       IF ( .NOT. ALLOCATED ( rad_lw_out_change_0 ) )  THEN
890          ALLOCATE ( rad_lw_out_change_0(nysg:nyng,nxlg:nxrg) )
891          rad_lw_out_change_0 = 0.0_wp
892       ENDIF
893
894!
895!--    Fix net radiation in case of radiation_scheme = 'constant'
896       IF ( radiation_scheme == 'constant' )  THEN
897          rad_net = net_radiation
898!          radiation = .FALSE.
899!
900!--    Calculate orbital constants
901       ELSE
902          decl_1 = SIN(23.45_wp * pi / 180.0_wp)
903          decl_2 = 2.0_wp * pi / 365.0_wp
904          decl_3 = decl_2 * 81.0_wp
905          lat    = phi * pi / 180.0_wp
906          lon    = lambda * pi / 180.0_wp
907       ENDIF
908
909       IF ( radiation_scheme == 'clear-sky'  .OR.                              &
910            radiation_scheme == 'constant')  THEN
911
912          ALLOCATE ( alpha(nysg:nyng,nxlg:nxrg) )
913
914          IF ( .NOT. ALLOCATED ( rad_sw_in ) )  THEN
915             ALLOCATE ( rad_sw_in(0:0,nysg:nyng,nxlg:nxrg) )
916          ENDIF
917          IF ( .NOT. ALLOCATED ( rad_sw_out ) )  THEN
918             ALLOCATE ( rad_sw_out(0:0,nysg:nyng,nxlg:nxrg) )
919          ENDIF
920
921          IF ( .NOT. ALLOCATED ( rad_sw_in_av ) )  THEN
922             ALLOCATE ( rad_sw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
923          ENDIF
924          IF ( .NOT. ALLOCATED ( rad_sw_out_av ) )  THEN
925             ALLOCATE ( rad_sw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
926          ENDIF
927
928          IF ( .NOT. ALLOCATED ( rad_lw_in ) )  THEN
929             ALLOCATE ( rad_lw_in(0:0,nysg:nyng,nxlg:nxrg) )
930          ENDIF
931          IF ( .NOT. ALLOCATED ( rad_lw_out ) )  THEN
932             ALLOCATE ( rad_lw_out(0:0,nysg:nyng,nxlg:nxrg) )
933          ENDIF
934
935          IF ( .NOT. ALLOCATED ( rad_lw_in_av ) )  THEN
936             ALLOCATE ( rad_lw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
937          ENDIF
938          IF ( .NOT. ALLOCATED ( rad_lw_out_av ) )  THEN
939             ALLOCATE ( rad_lw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
940          ENDIF
941
942          rad_sw_in  = 0.0_wp
943          rad_sw_out = 0.0_wp
944          rad_lw_in  = 0.0_wp
945          rad_lw_out = 0.0_wp
946
947!
948!--       Overwrite albedo if manually set in parameter file
949          IF ( albedo_type /= 0 .AND. albedo == 9999999.9_wp )  THEN
950             albedo = albedo_pars(2,albedo_type)
951          ENDIF
952   
953          alpha = albedo
954 
955!
956!--    Initialization actions for RRTMG
957       ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
958#if defined ( __rrtmg )
959!
960!--       Allocate albedos
961          ALLOCATE ( rrtm_aldif(0:0,nysg:nyng,nxlg:nxrg) )
962          ALLOCATE ( rrtm_aldir(0:0,nysg:nyng,nxlg:nxrg) )
963          ALLOCATE ( rrtm_asdif(0:0,nysg:nyng,nxlg:nxrg) )
964          ALLOCATE ( rrtm_asdir(0:0,nysg:nyng,nxlg:nxrg) )
965          ALLOCATE ( aldif(nysg:nyng,nxlg:nxrg) )
966          ALLOCATE ( aldir(nysg:nyng,nxlg:nxrg) )
967          ALLOCATE ( asdif(nysg:nyng,nxlg:nxrg) )
968          ALLOCATE ( asdir(nysg:nyng,nxlg:nxrg) )
969
970          IF ( albedo_type /= 0 )  THEN
971             IF ( albedo_lw_dif == 9999999.9_wp )  THEN
972                albedo_lw_dif = albedo_pars(0,albedo_type)
973                albedo_lw_dir = albedo_lw_dif
974             ENDIF
975             IF ( albedo_sw_dif == 9999999.9_wp )  THEN
976                albedo_sw_dif = albedo_pars(1,albedo_type)
977                albedo_sw_dir = albedo_sw_dif
978             ENDIF
979          ENDIF
980
981          aldif(:,:) = albedo_lw_dif
982          aldir(:,:) = albedo_lw_dir
983          asdif(:,:) = albedo_sw_dif
984          asdir(:,:) = albedo_sw_dir
985!
986!--       Calculate initial values of current (cosine of) the zenith angle and
987!--       whether the sun is up
988          CALL calc_zenith     
989!
990!--       Calculate initial surface albedo
991          IF ( .NOT. constant_albedo )  THEN
992             CALL calc_albedo
993          ELSE
994             rrtm_aldif(0,:,:) = aldif(:,:)
995             rrtm_aldir(0,:,:) = aldir(:,:)
996             rrtm_asdif(0,:,:) = asdif(:,:) 
997             rrtm_asdir(0,:,:) = asdir(:,:)   
998          ENDIF
999
1000!
1001!--       Allocate surface emissivity
1002          ALLOCATE ( rrtm_emis(0:0,1:nbndlw+1) )
1003          rrtm_emis = emissivity
1004
1005!
1006!--       Allocate 3d arrays of radiative fluxes and heating rates
1007          IF ( .NOT. ALLOCATED ( rad_sw_in ) )  THEN
1008             ALLOCATE ( rad_sw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1009             rad_sw_in = 0.0_wp
1010          ENDIF
1011
1012          IF ( .NOT. ALLOCATED ( rad_sw_in_av ) )  THEN
1013             ALLOCATE ( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1014          ENDIF
1015
1016          IF ( .NOT. ALLOCATED ( rad_sw_out ) )  THEN
1017             ALLOCATE ( rad_sw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1018             rad_sw_out = 0.0_wp
1019          ENDIF
1020
1021          IF ( .NOT. ALLOCATED ( rad_sw_out_av ) )  THEN
1022             ALLOCATE ( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1023          ENDIF
1024
1025          IF ( .NOT. ALLOCATED ( rad_sw_hr ) )  THEN
1026             ALLOCATE ( rad_sw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1027             rad_sw_hr = 0.0_wp
1028          ENDIF
1029
1030          IF ( .NOT. ALLOCATED ( rad_sw_hr_av ) )  THEN
1031             ALLOCATE ( rad_sw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1032             rad_sw_hr_av = 0.0_wp
1033          ENDIF
1034
1035          IF ( .NOT. ALLOCATED ( rad_sw_cs_hr ) )  THEN
1036             ALLOCATE ( rad_sw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1037             rad_sw_cs_hr = 0.0_wp
1038          ENDIF
1039
1040          IF ( .NOT. ALLOCATED ( rad_sw_cs_hr_av ) )  THEN
1041             ALLOCATE ( rad_sw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1042             rad_sw_cs_hr_av = 0.0_wp
1043          ENDIF
1044
1045          IF ( .NOT. ALLOCATED ( rad_lw_in ) )  THEN
1046             ALLOCATE ( rad_lw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1047             rad_lw_in     = 0.0_wp
1048          ENDIF
1049
1050          IF ( .NOT. ALLOCATED ( rad_lw_in_av ) )  THEN
1051             ALLOCATE ( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1052          ENDIF
1053
1054          IF ( .NOT. ALLOCATED ( rad_lw_out ) )  THEN
1055             ALLOCATE ( rad_lw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1056            rad_lw_out    = 0.0_wp
1057          ENDIF
1058
1059          IF ( .NOT. ALLOCATED ( rad_lw_out_av ) )  THEN
1060             ALLOCATE ( rad_lw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1061          ENDIF
1062
1063          IF ( .NOT. ALLOCATED ( rad_lw_hr ) )  THEN
1064             ALLOCATE ( rad_lw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1065             rad_lw_hr = 0.0_wp
1066          ENDIF
1067
1068          IF ( .NOT. ALLOCATED ( rad_lw_hr_av ) )  THEN
1069             ALLOCATE ( rad_lw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1070             rad_lw_hr_av = 0.0_wp
1071          ENDIF
1072
1073          IF ( .NOT. ALLOCATED ( rad_lw_cs_hr ) )  THEN
1074             ALLOCATE ( rad_lw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1075             rad_lw_cs_hr = 0.0_wp
1076          ENDIF
1077
1078          IF ( .NOT. ALLOCATED ( rad_lw_cs_hr_av ) )  THEN
1079             ALLOCATE ( rad_lw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1080             rad_lw_cs_hr_av = 0.0_wp
1081          ENDIF
1082
1083          ALLOCATE ( rad_sw_cs_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1084          ALLOCATE ( rad_sw_cs_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1085          rad_sw_cs_in  = 0.0_wp
1086          rad_sw_cs_out = 0.0_wp
1087
1088          ALLOCATE ( rad_lw_cs_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1089          ALLOCATE ( rad_lw_cs_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1090          rad_lw_cs_in  = 0.0_wp
1091          rad_lw_cs_out = 0.0_wp
1092
1093!
1094!--       Allocate dummy array for storing surface temperature
1095          ALLOCATE ( rrtm_tsfc(1) )
1096
1097!
1098!--       Initialize RRTMG
1099          IF ( lw_radiation )  CALL rrtmg_lw_ini ( cp )
1100          IF ( sw_radiation )  CALL rrtmg_sw_ini ( cp )
1101
1102!
1103!--       Set input files for RRTMG
1104          INQUIRE(FILE="RAD_SND_DATA", EXIST=snd_exists) 
1105          IF ( .NOT. snd_exists )  THEN
1106             rrtm_input_file = "rrtmg_lw.nc"
1107          ENDIF
1108
1109!
1110!--       Read vertical layers for RRTMG from sounding data
1111!--       The routine provides nzt_rad, hyp_snd(1:nzt_rad),
1112!--       t_snd(nzt+2:nzt_rad), rrtm_play(1:nzt_rad), rrtm_plev(1_nzt_rad+1),
1113!--       rrtm_tlay(nzt+2:nzt_rad), rrtm_tlev(nzt+2:nzt_rad+1)
1114          CALL read_sounding_data
1115
1116!
1117!--       Read trace gas profiles from file. This routine provides
1118!--       the rrtm_ arrays (1:nzt_rad+1)
1119          CALL read_trace_gas_data
1120#endif
1121       ENDIF
1122
1123!
1124!--    Perform user actions if required
1125       CALL user_init_radiation
1126
1127!
1128!--    Calculate radiative fluxes at model start
1129       IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
1130
1131          SELECT CASE ( radiation_scheme )
1132             CASE ( 'rrtmg' )
1133                CALL radiation_rrtmg
1134             CASE ( 'clear-sky' )
1135                CALL radiation_clearsky
1136             CASE ( 'constant' )
1137                CALL radiation_constant
1138             CASE DEFAULT
1139          END SELECT
1140
1141       ENDIF
1142
1143       RETURN
1144
1145    END SUBROUTINE radiation_init
1146
1147
1148!------------------------------------------------------------------------------!
1149! Description:
1150! ------------
1151!> A simple clear sky radiation model
1152!------------------------------------------------------------------------------!
1153    SUBROUTINE radiation_clearsky
1154
1155
1156       IMPLICIT NONE
1157
1158       INTEGER(iwp) :: i, j, k   !< loop indices
1159       REAL(wp)     :: exn,   &  !< Exner functions at surface
1160                       exn1,  &  !< Exner functions at first grid level
1161                       pt1       !< potential temperature at first grid level
1162
1163!
1164!--    Calculate current zenith angle
1165       CALL calc_zenith
1166
1167!
1168!--    Calculate sky transmissivity
1169       sky_trans = 0.6_wp + 0.2_wp * zenith(0)
1170
1171!
1172!--    Calculate value of the Exner function
1173       exn = (surface_pressure / 1000.0_wp )**0.286_wp
1174!
1175!--    Calculate radiation fluxes and net radiation (rad_net) for each grid
1176!--    point
1177       DO i = nxlg, nxrg
1178          DO j = nysg, nyng
1179!
1180!--          Obtain vertical index of topography top
1181             k = MAXLOC(                                                       &
1182                        MERGE( 1, 0,                                           &
1183                               BTEST( wall_flags_0(:,j,i), 12 )                &
1184                             ), DIM = 1                                        &
1185                       ) - 1
1186
1187             exn1 = (hyp(k+1) / 100000.0_wp )**0.286_wp
1188
1189             rad_sw_in(0,j,i)  = solar_constant * sky_trans * zenith(0)
1190             rad_sw_out(0,j,i) = alpha(j,i) * rad_sw_in(0,j,i)
1191             rad_lw_out(0,j,i) = emissivity * sigma_sb * (pt(k,j,i) * exn)**4
1192
1193             IF ( cloud_physics )  THEN
1194                pt1 = pt(k+1,j,i) + l_d_cp / exn1 * ql(k+1,j,i)
1195                rad_lw_in(0,j,i)  = 0.8_wp * sigma_sb * (pt1 * exn1)**4
1196             ELSE
1197                rad_lw_in(0,j,i)  = 0.8_wp * sigma_sb * (pt(k+1,j,i) * exn1)**4
1198             ENDIF
1199
1200             rad_net(j,i) = rad_sw_in(0,j,i) - rad_sw_out(0,j,i)               &
1201                            + rad_lw_in(0,j,i) - rad_lw_out(0,j,i)
1202
1203
1204             rad_lw_out_change_0(j,i) = 3.0_wp * sigma_sb * emissivity         &
1205                                        * (pt(k,j,i) * exn) ** 3
1206
1207          ENDDO
1208       ENDDO
1209
1210    END SUBROUTINE radiation_clearsky
1211
1212
1213!------------------------------------------------------------------------------!
1214! Description:
1215! ------------
1216!> This scheme keeps the prescribed net radiation constant during the run
1217!------------------------------------------------------------------------------!
1218    SUBROUTINE radiation_constant
1219
1220
1221       IMPLICIT NONE
1222
1223       INTEGER(iwp) :: i, j, k   !< loop indices
1224       REAL(wp)     :: exn,   &  !< Exner functions at surface
1225                       exn1,  &  !< Exner functions at first grid level
1226                       pt1       !< potential temperature at first grid level
1227
1228!
1229!--    Calculate value of the Exner function
1230       exn = (surface_pressure / 1000.0_wp )**0.286_wp
1231!
1232!--    Prescribe net radiation and estimate the remaining radiative fluxes
1233       DO i = nxlg, nxrg
1234          DO j = nysg, nyng
1235!
1236!--          Obtain vertical index of topography top. So far it is identical to nzb.
1237             k = MAXLOC(                                                       &
1238                        MERGE( 1, 0,                                           &
1239                               BTEST( wall_flags_0(:,j,i), 12 )                &
1240                             ), DIM = 1                                        &
1241                       ) - 1
1242
1243             rad_net(j,i)      = net_radiation
1244
1245             exn1 = (hyp(k+1) / 100000.0_wp )**0.286_wp
1246
1247             IF ( cloud_physics )  THEN
1248                pt1 = pt(k+1,j,i) + l_d_cp / exn1 * ql(k+1,j,i)
1249                rad_lw_in(0,j,i)  = 0.8_wp * sigma_sb * (pt1 * exn1)**4
1250             ELSE
1251                rad_lw_in(0,j,i)  = 0.8_wp * sigma_sb * (pt(k+1,j,i) * exn1)**4
1252             ENDIF
1253
1254             rad_lw_out(0,j,i) = emissivity * sigma_sb * (pt(k,j,i) * exn)**4
1255
1256             rad_sw_in(0,j,i) = ( rad_net(j,i) - rad_lw_in(0,j,i)              &
1257                                  + rad_lw_out(0,j,i) )                        &
1258                                  / ( 1.0_wp - alpha(j,i) )
1259
1260             rad_sw_out(0,j,i) =  alpha(j,i) * rad_sw_in(0,j,i)
1261
1262          ENDDO
1263       ENDDO
1264
1265    END SUBROUTINE radiation_constant
1266
1267!------------------------------------------------------------------------------!
1268! Description:
1269! ------------
1270!> Header output for radiation model
1271!------------------------------------------------------------------------------!
1272    SUBROUTINE radiation_header ( io )
1273
1274
1275       IMPLICIT NONE
1276 
1277       INTEGER(iwp), INTENT(IN) ::  io            !< Unit of the output file
1278   
1279
1280       
1281!
1282!--    Write radiation model header
1283       WRITE( io, 3 )
1284
1285       IF ( radiation_scheme == "constant" )  THEN
1286          WRITE( io, 4 ) net_radiation
1287       ELSEIF ( radiation_scheme == "clear-sky" )  THEN
1288          WRITE( io, 5 )
1289       ELSEIF ( radiation_scheme == "rrtmg" )  THEN
1290          WRITE( io, 6 )
1291          IF ( .NOT. lw_radiation )  WRITE( io, 10 )
1292          IF ( .NOT. sw_radiation )  WRITE( io, 11 )
1293       ENDIF
1294
1295       IF ( albedo_type == 0 )  THEN
1296          WRITE( io, 7 ) albedo
1297       ELSE
1298          WRITE( io, 8 ) TRIM( albedo_type_name(albedo_type) )
1299       ENDIF
1300       IF ( constant_albedo )  THEN
1301          WRITE( io, 9 )
1302       ENDIF
1303       
1304       IF ( radiation .AND. radiation_scheme /= 'constant' )  THEN
1305          WRITE ( io, 1 )  lambda
1306          WRITE ( io, 2 )  day_init, time_utc_init
1307       ENDIF
1308
1309       WRITE( io, 12 ) dt_radiation
1310 
1311
1312 1 FORMAT ('    Geograph. longitude            :   lambda = ',F4.1,' degr')
1313 2 FORMAT ('    Day of the year at model start :   day_init = ',I3   &
1314            /'    UTC time at model start        :   time_utc_init = ',F7.1' s')
1315 3 FORMAT (//' Radiation model information:'/                                  &
1316              ' ----------------------------'/)
1317 4 FORMAT ('    --> Using constant net radiation: net_radiation = ', F6.2,        &
1318           // 'W/m**2')
1319 5 FORMAT ('    --> Simple radiation scheme for clear sky is used (no clouds,',   &
1320                   ' default)')
1321 6 FORMAT ('    --> RRTMG scheme is used')
1322 7 FORMAT (/'    User-specific surface albedo: albedo =', F6.3)
1323 8 FORMAT (/'    Albedo is set for land surface type: ', A)
1324 9 FORMAT (/'    --> Albedo is fixed during the run')
132510 FORMAT (/'    --> Longwave radiation is disabled')
132611 FORMAT (/'    --> Shortwave radiation is disabled.')
132712 FORMAT  ('    Timestep: dt_radiation = ', F6.2, '  s')
1328
1329
1330    END SUBROUTINE radiation_header
1331   
1332
1333!------------------------------------------------------------------------------!
1334! Description:
1335! ------------
1336!> Parin for &radiation_par for radiation model
1337!------------------------------------------------------------------------------!
1338    SUBROUTINE radiation_parin
1339
1340
1341       IMPLICIT NONE
1342
1343       CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file
1344       
1345       NAMELIST /radiation_par/   albedo, albedo_type, albedo_lw_dir,          &
1346                                  albedo_lw_dif, albedo_sw_dir, albedo_sw_dif, &
1347                                  constant_albedo, day_init, dt_radiation,     &
1348                                  lambda, lw_radiation, net_radiation,         &
1349                                  radiation_scheme, skip_time_do_radiation,    &
1350                                  sw_radiation, time_utc_init,                 &
1351                                  unscheduled_radiation_calls
1352       
1353       line = ' '
1354       
1355!
1356!--    Try to find radiation model package
1357       REWIND ( 11 )
1358       line = ' '
1359       DO   WHILE ( INDEX( line, '&radiation_par' ) == 0 )
1360          READ ( 11, '(A)', END=10 )  line
1361       ENDDO
1362       BACKSPACE ( 11 )
1363
1364!
1365!--    Read user-defined namelist
1366       READ ( 11, radiation_par )
1367
1368!
1369!--    Set flag that indicates that the radiation model is switched on
1370       radiation = .TRUE.
1371
1372 10    CONTINUE
1373       
1374
1375    END SUBROUTINE radiation_parin
1376
1377
1378!------------------------------------------------------------------------------!
1379! Description:
1380! ------------
1381!> Implementation of the RRTMG radiation_scheme
1382!------------------------------------------------------------------------------!
1383    SUBROUTINE radiation_rrtmg
1384
1385       USE indices,                                                            &
1386           ONLY:  nbgp
1387
1388       USE particle_attributes,                                                &
1389           ONLY:  grid_particles, number_of_particles, particles,              &
1390                  particle_advection_start, prt_count
1391
1392       IMPLICIT NONE
1393
1394#if defined ( __rrtmg )
1395
1396       INTEGER(iwp) :: i, j, k, n !< loop indices
1397
1398       REAL(wp)     ::  s_r2, &   !< weighted sum over all droplets with r^2
1399                        s_r3      !< weighted sum over all droplets with r^3
1400
1401!
1402!--    Calculate current (cosine of) zenith angle and whether the sun is up
1403       CALL calc_zenith     
1404!
1405!--    Calculate surface albedo
1406       IF ( .NOT. constant_albedo )  THEN
1407          CALL calc_albedo
1408       ENDIF
1409
1410!
1411!--    Prepare input data for RRTMG
1412
1413!
1414!--    In case of large scale forcing with surface data, calculate new pressure
1415!--    profile. nzt_rad might be modified by these calls and all required arrays
1416!--    will then be re-allocated
1417       IF ( large_scale_forcing  .AND.  lsf_surf )  THEN
1418          CALL read_sounding_data
1419          CALL read_trace_gas_data
1420       ENDIF
1421!
1422!--    Loop over all grid points
1423       DO i = nxl, nxr
1424          DO j = nys, nyn
1425
1426!
1427!--          Prepare profiles of temperature and H2O volume mixing ratio
1428             rrtm_tlev(0,nzb+1) = pt(nzb,j,i) * ( surface_pressure             &
1429                                                  / 1000.0_wp )**0.286_wp
1430
1431
1432             IF ( cloud_physics )  THEN
1433                DO k = nzb+1, nzt+1
1434                   rrtm_tlay(0,k) = pt(k,j,i) * ( (hyp(k) ) / 100000.0_wp      &
1435                                    )**0.286_wp + l_d_cp * ql(k,j,i)
1436                   rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * (q(k,j,i) - ql(k,j,i))
1437                ENDDO
1438             ELSE
1439                DO k = nzb+1, nzt+1
1440                   rrtm_tlay(0,k) = pt(k,j,i) * ( (hyp(k) ) / 100000.0_wp      &
1441                                    )**0.286_wp
1442                   rrtm_h2ovmr(0,k) = 0.0_wp
1443                ENDDO
1444             ENDIF
1445
1446!
1447!--          Avoid temperature/humidity jumps at the top of the LES domain by
1448!--          linear interpolation from nzt+2 to nzt+7
1449             DO k = nzt+2, nzt+7
1450                rrtm_tlay(0,k) = rrtm_tlay(0,nzt+1)                            &
1451                              + ( rrtm_tlay(0,nzt+8) - rrtm_tlay(0,nzt+1) )    &
1452                              / ( rrtm_play(0,nzt+8) - rrtm_play(0,nzt+1) )    &
1453                              * ( rrtm_play(0,k) - rrtm_play(0,nzt+1) )
1454
1455                rrtm_h2ovmr(0,k) = rrtm_h2ovmr(0,nzt+1)                        &
1456                              + ( rrtm_h2ovmr(0,nzt+8) - rrtm_h2ovmr(0,nzt+1) )&
1457                              / ( rrtm_play(0,nzt+8)   - rrtm_play(0,nzt+1)   )&
1458                              * ( rrtm_play(0,k) - rrtm_play(0,nzt+1) )
1459
1460             ENDDO
1461
1462!--          Linear interpolate to zw grid
1463             DO k = nzb+2, nzt+8
1464                rrtm_tlev(0,k)   = rrtm_tlay(0,k-1) + (rrtm_tlay(0,k) -        &
1465                                   rrtm_tlay(0,k-1))                           &
1466                                   / ( rrtm_play(0,k) - rrtm_play(0,k-1) )     &
1467                                   * ( rrtm_plev(0,k) - rrtm_play(0,k-1) )
1468             ENDDO
1469
1470
1471!
1472!--          Calculate liquid water path and cloud fraction for each column.
1473!--          Note that LWP is required in g/m² instead of kg/kg m.
1474             rrtm_cldfr  = 0.0_wp
1475             rrtm_reliq  = 0.0_wp
1476             rrtm_cliqwp = 0.0_wp
1477             rrtm_icld   = 0
1478
1479             IF ( cloud_physics )  THEN
1480                DO k = nzb+1, nzt+1
1481                   rrtm_cliqwp(0,k) =  ql(k,j,i) * 1000.0_wp *                 &
1482                                       (rrtm_plev(0,k) - rrtm_plev(0,k+1))     &
1483                                       * 100.0_wp / g 
1484
1485                   IF ( rrtm_cliqwp(0,k) > 0.0_wp )  THEN
1486                      rrtm_cldfr(0,k) = 1.0_wp
1487                      IF ( rrtm_icld == 0 )  rrtm_icld = 1
1488
1489!
1490!--                   Calculate cloud droplet effective radius
1491                      IF ( cloud_physics )  THEN
1492                         rrtm_reliq(0,k) = 1.0E6_wp * ( 3.0_wp * ql(k,j,i)        &
1493                                           * rho_surface                          &
1494                                           / ( 4.0_wp * pi * nc_const * rho_l )   &
1495                                           )**0.33333333333333_wp                 &
1496                                           * EXP( LOG( sigma_gc )**2 )
1497
1498                      ELSEIF ( cloud_droplets )  THEN
1499                         number_of_particles = prt_count(k,j,i)
1500
1501                         IF (number_of_particles <= 0)  CYCLE
1502                         particles => grid_particles(k,j,i)%particles(1:number_of_particles)
1503                         s_r2 = 0.0_wp
1504                         s_r3 = 0.0_wp
1505
1506                         DO  n = 1, number_of_particles
1507                            IF ( particles(n)%particle_mask )  THEN
1508                               s_r2 = s_r2 + particles(n)%radius**2 * &
1509                                      particles(n)%weight_factor
1510                               s_r3 = s_r3 + particles(n)%radius**3 * &
1511                                      particles(n)%weight_factor
1512                            ENDIF
1513                         ENDDO
1514
1515                         IF ( s_r2 > 0.0_wp )  rrtm_reliq(0,k) = s_r3 / s_r2
1516
1517                      ENDIF
1518
1519!
1520!--                   Limit effective radius
1521                      IF ( rrtm_reliq(0,k) > 0.0_wp )  THEN
1522                         rrtm_reliq(0,k) = MAX(rrtm_reliq(0,k),2.5_wp)
1523                         rrtm_reliq(0,k) = MIN(rrtm_reliq(0,k),60.0_wp)
1524                     ENDIF
1525                   ENDIF
1526                ENDDO
1527             ENDIF
1528
1529!
1530!--          Set surface temperature
1531             rrtm_tsfc = pt(nzb,j,i) * (surface_pressure / 1000.0_wp )**0.286_wp
1532
1533             IF ( lw_radiation )  THEN
1534               CALL rrtmg_lw( 1, nzt_rad      , rrtm_icld    , rrtm_idrv      ,&
1535               rrtm_play       , rrtm_plev    , rrtm_tlay    , rrtm_tlev      ,&
1536               rrtm_tsfc       , rrtm_h2ovmr  , rrtm_o3vmr   , rrtm_co2vmr    ,&
1537               rrtm_ch4vmr     , rrtm_n2ovmr  , rrtm_o2vmr   , rrtm_cfc11vmr  ,&
1538               rrtm_cfc12vmr   , rrtm_cfc22vmr, rrtm_ccl4vmr , rrtm_emis      ,&
1539               rrtm_inflglw    , rrtm_iceflglw, rrtm_liqflglw, rrtm_cldfr     ,&
1540               rrtm_lw_taucld  , rrtm_cicewp  , rrtm_cliqwp  , rrtm_reice     ,& 
1541               rrtm_reliq      , rrtm_lw_tauaer,                               &
1542               rrtm_lwuflx     , rrtm_lwdflx  , rrtm_lwhr  ,                   &
1543               rrtm_lwuflxc    , rrtm_lwdflxc , rrtm_lwhrc ,                   &
1544               rrtm_lwuflx_dt  ,  rrtm_lwuflxc_dt )
1545
1546!
1547!--             Save fluxes
1548                DO k = nzb, nzt+1
1549                   rad_lw_in(k,j,i)  = rrtm_lwdflx(0,k)
1550                   rad_lw_out(k,j,i) = rrtm_lwuflx(0,k)
1551                ENDDO
1552
1553!
1554!--             Save heating rates (convert from K/d to K/h)
1555                DO k = nzb+1, nzt+1
1556                   rad_lw_hr(k,j,i)     = rrtm_lwhr(0,k)  * d_hours_day
1557                   rad_lw_cs_hr(k,j,i)  = rrtm_lwhrc(0,k) * d_hours_day
1558                ENDDO
1559
1560!
1561!--             Save change in LW heating rate
1562                rad_lw_out_change_0(j,i) = rrtm_lwuflx_dt(0,nzb)
1563
1564             ENDIF
1565
1566             IF ( sw_radiation .AND. sun_up )  THEN
1567                CALL rrtmg_sw( 1, nzt_rad      , rrtm_icld  , rrtm_iaer       ,&
1568               rrtm_play       , rrtm_plev    , rrtm_tlay  , rrtm_tlev        ,&
1569               rrtm_tsfc       , rrtm_h2ovmr  , rrtm_o3vmr , rrtm_co2vmr      ,&
1570               rrtm_ch4vmr     , rrtm_n2ovmr  , rrtm_o2vmr , rrtm_asdir(:,j,i),&
1571               rrtm_asdif(:,j,i), rrtm_aldir(:,j,i), rrtm_aldif(:,j,i), zenith,&
1572               0.0_wp          , day          , solar_constant,   rrtm_inflgsw,&
1573               rrtm_iceflgsw   , rrtm_liqflgsw, rrtm_cldfr , rrtm_sw_taucld   ,&
1574               rrtm_sw_ssacld  , rrtm_sw_asmcld, rrtm_sw_fsfcld, rrtm_cicewp  ,&
1575               rrtm_cliqwp     , rrtm_reice   , rrtm_reliq , rrtm_sw_tauaer   ,&
1576               rrtm_sw_ssaaer     , rrtm_sw_asmaer  , rrtm_sw_ecaer ,          &
1577               rrtm_swuflx     , rrtm_swdflx  , rrtm_swhr  ,                   &
1578               rrtm_swuflxc    , rrtm_swdflxc , rrtm_swhrc )
1579 
1580!
1581!--             Save fluxes
1582                DO k = nzb, nzt+1
1583                   rad_sw_in(k,j,i)  = rrtm_swdflx(0,k)
1584                   rad_sw_out(k,j,i) = rrtm_swuflx(0,k)
1585                ENDDO
1586
1587!
1588!--             Save heating rates (convert from K/d to K/s)
1589                DO k = nzb+1, nzt+1
1590                   rad_sw_hr(k,j,i)     = rrtm_swhr(0,k)  * d_hours_day
1591                   rad_sw_cs_hr(k,j,i)  = rrtm_swhrc(0,k) * d_hours_day
1592                ENDDO
1593
1594             ENDIF
1595
1596!
1597!--          Calculate surface net radiation
1598             rad_net(j,i) = rad_sw_in(nzb,j,i) - rad_sw_out(nzb,j,i)           &
1599                            + rad_lw_in(nzb,j,i) - rad_lw_out(nzb,j,i)
1600
1601          ENDDO
1602       ENDDO
1603
1604       CALL exchange_horiz( rad_lw_in,  nbgp )
1605       CALL exchange_horiz( rad_lw_out, nbgp )
1606       CALL exchange_horiz( rad_lw_hr,    nbgp )
1607       CALL exchange_horiz( rad_lw_cs_hr, nbgp )
1608
1609       CALL exchange_horiz( rad_sw_in,  nbgp )
1610       CALL exchange_horiz( rad_sw_out, nbgp ) 
1611       CALL exchange_horiz( rad_sw_hr,    nbgp )
1612       CALL exchange_horiz( rad_sw_cs_hr, nbgp )
1613
1614       CALL exchange_horiz_2d( rad_net )
1615       CALL exchange_horiz_2d( rad_lw_out_change_0 )
1616#endif
1617
1618    END SUBROUTINE radiation_rrtmg
1619
1620
1621!------------------------------------------------------------------------------!
1622! Description:
1623! ------------
1624!> Calculate the cosine of the zenith angle (variable is called zenith)
1625!------------------------------------------------------------------------------!
1626    SUBROUTINE calc_zenith
1627
1628       IMPLICIT NONE
1629
1630       REAL(wp) ::  declination,  & !< solar declination angle
1631                    hour_angle      !< solar hour angle
1632!
1633!--    Calculate current day and time based on the initial values and simulation
1634!--    time
1635       day = day_init + INT(FLOOR( (time_utc_init + time_since_reference_point)    &
1636                               / 86400.0_wp ), KIND=iwp)
1637       time_utc = MOD((time_utc_init + time_since_reference_point), 86400.0_wp)
1638
1639
1640!
1641!--    Calculate solar declination and hour angle   
1642       declination = ASIN( decl_1 * SIN(decl_2 * REAL(day, KIND=wp) - decl_3) )
1643       hour_angle  = 2.0_wp * pi * (time_utc / 86400.0_wp) + lon - pi
1644
1645!
1646!--    Calculate cosine of solar zenith angle
1647       zenith(0) = SIN(lat) * SIN(declination) + COS(lat) * COS(declination)      &
1648                                            * COS(hour_angle)
1649       zenith(0) = MAX(0.0_wp,zenith(0))
1650
1651!
1652!--    Calculate solar directional vector
1653       IF ( sun_direction )  THEN
1654!--       Direction in longitudes equals to sin(solar_azimuth) * sin(zenith)
1655          sun_dir_lon(0) = -SIN(hour_angle) * COS(declination)
1656!--       Direction in latitues equals to cos(solar_azimuth) * sin(zenith)
1657          sun_dir_lat(0) = SIN(declination) * COS(lat) - COS(hour_angle) &
1658                              * COS(declination) * SIN(lat)
1659       ENDIF
1660
1661!
1662!--    Check if the sun is up (otheriwse shortwave calculations can be skipped)
1663       IF ( zenith(0) > 0.0_wp )  THEN
1664          sun_up = .TRUE.
1665       ELSE
1666          sun_up = .FALSE.
1667       END IF
1668
1669    END SUBROUTINE calc_zenith
1670
1671#if defined ( __rrtmg ) && defined ( __netcdf )
1672!------------------------------------------------------------------------------!
1673! Description:
1674! ------------
1675!> Calculates surface albedo components based on Briegleb (1992) and
1676!> Briegleb et al. (1986)
1677!------------------------------------------------------------------------------!
1678    SUBROUTINE calc_albedo
1679
1680        IMPLICIT NONE
1681
1682        IF ( sun_up )  THEN
1683!
1684!--        Ocean
1685           IF ( albedo_type == 1 )  THEN
1686              rrtm_aldir(0,:,:) = 0.026_wp / ( zenith(0)**1.7_wp + 0.065_wp )  &
1687                                  + 0.15_wp * ( zenith(0) - 0.1_wp )           &
1688                                            * ( zenith(0) - 0.5_wp )           &
1689                                            * ( zenith(0) - 1.0_wp )
1690              rrtm_asdir(0,:,:) = rrtm_aldir(0,:,:)
1691!
1692!--        Snow
1693           ELSEIF ( albedo_type == 16 )  THEN
1694              IF ( zenith(0) < 0.5_wp )  THEN
1695                 rrtm_aldir(0,:,:) = 0.5_wp * (1.0_wp - aldif)                 &
1696                                     * ( 3.0_wp / (1.0_wp + 4.0_wp             &
1697                                     * zenith(0))) - 1.0_wp
1698                 rrtm_asdir(0,:,:) = 0.5_wp * (1.0_wp - asdif)                 &
1699                                     * ( 3.0_wp / (1.0_wp + 4.0_wp             &
1700                                     * zenith(0))) - 1.0_wp
1701
1702                 rrtm_aldir(0,:,:) = MIN(0.98_wp, rrtm_aldir(0,:,:))
1703                 rrtm_asdir(0,:,:) = MIN(0.98_wp, rrtm_asdir(0,:,:))
1704              ELSE
1705                 rrtm_aldir(0,:,:) = aldif
1706                 rrtm_asdir(0,:,:) = asdif
1707              ENDIF
1708!
1709!--        Sea ice
1710           ELSEIF ( albedo_type == 15 )  THEN
1711                 rrtm_aldir(0,:,:) = aldif
1712                 rrtm_asdir(0,:,:) = asdif
1713
1714!
1715!--        Asphalt
1716           ELSEIF ( albedo_type == 17 )  THEN
1717                 rrtm_aldir(0,:,:) = aldif
1718                 rrtm_asdir(0,:,:) = asdif
1719!
1720!--        Land surfaces
1721           ELSE
1722              SELECT CASE ( albedo_type )
1723
1724!
1725!--              Surface types with strong zenith dependence
1726                 CASE ( 1, 2, 3, 4, 11, 12, 13 )
1727                    rrtm_aldir(0,:,:) = aldif * 1.4_wp /                       &
1728                                        (1.0_wp + 0.8_wp * zenith(0))
1729                    rrtm_asdir(0,:,:) = asdif * 1.4_wp /                       &
1730                                        (1.0_wp + 0.8_wp * zenith(0))
1731!
1732!--              Surface types with weak zenith dependence
1733                 CASE ( 5, 6, 7, 8, 9, 10, 14 )
1734                    rrtm_aldir(0,:,:) = aldif * 1.1_wp /                       &
1735                                        (1.0_wp + 0.2_wp * zenith(0))
1736                    rrtm_asdir(0,:,:) = asdif * 1.1_wp /                       &
1737                                        (1.0_wp + 0.2_wp * zenith(0))
1738
1739                 CASE DEFAULT
1740
1741              END SELECT
1742           ENDIF
1743!
1744!--        Diffusive albedo is taken from Table 2
1745           rrtm_aldif(0,:,:) = aldif
1746           rrtm_asdif(0,:,:) = asdif
1747
1748        ELSE
1749
1750           rrtm_aldir(0,:,:) = 0.0_wp
1751           rrtm_asdir(0,:,:) = 0.0_wp
1752           rrtm_aldif(0,:,:) = 0.0_wp
1753           rrtm_asdif(0,:,:) = 0.0_wp
1754        ENDIF
1755    END SUBROUTINE calc_albedo
1756
1757!------------------------------------------------------------------------------!
1758! Description:
1759! ------------
1760!> Read sounding data (pressure and temperature) from RADIATION_DATA.
1761!------------------------------------------------------------------------------!
1762    SUBROUTINE read_sounding_data
1763
1764       IMPLICIT NONE
1765
1766       INTEGER(iwp) :: id,           & !< NetCDF id of input file
1767                       id_dim_zrad,  & !< pressure level id in the NetCDF file
1768                       id_var,       & !< NetCDF variable id
1769                       k,            & !< loop index
1770                       nz_snd,       & !< number of vertical levels in the sounding data
1771                       nz_snd_start, & !< start vertical index for sounding data to be used
1772                       nz_snd_end      !< end vertical index for souding data to be used
1773
1774       REAL(wp) :: t_surface           !< actual surface temperature
1775
1776       REAL(wp), DIMENSION(:), ALLOCATABLE ::  hyp_snd_tmp, & !< temporary hydrostatic pressure profile (sounding)
1777                                               t_snd_tmp      !< temporary temperature profile (sounding)
1778
1779!
1780!--    In case of updates, deallocate arrays first (sufficient to check one
1781!--    array as the others are automatically allocated). This is required
1782!--    because nzt_rad might change during the update
1783       IF ( ALLOCATED ( hyp_snd ) )  THEN
1784          DEALLOCATE( hyp_snd )
1785          DEALLOCATE( t_snd )
1786          DEALLOCATE( q_snd  )
1787          DEALLOCATE ( rrtm_play )
1788          DEALLOCATE ( rrtm_plev )
1789          DEALLOCATE ( rrtm_tlay )
1790          DEALLOCATE ( rrtm_tlev )
1791
1792          DEALLOCATE ( rrtm_h2ovmr )
1793          DEALLOCATE ( rrtm_cicewp )
1794          DEALLOCATE ( rrtm_cldfr )
1795          DEALLOCATE ( rrtm_cliqwp )
1796          DEALLOCATE ( rrtm_reice )
1797          DEALLOCATE ( rrtm_reliq )
1798          DEALLOCATE ( rrtm_lw_taucld )
1799          DEALLOCATE ( rrtm_lw_tauaer )
1800
1801          DEALLOCATE ( rrtm_lwdflx  )
1802          DEALLOCATE ( rrtm_lwdflxc )
1803          DEALLOCATE ( rrtm_lwuflx  )
1804          DEALLOCATE ( rrtm_lwuflxc )
1805          DEALLOCATE ( rrtm_lwuflx_dt )
1806          DEALLOCATE ( rrtm_lwuflxc_dt )
1807          DEALLOCATE ( rrtm_lwhr  )
1808          DEALLOCATE ( rrtm_lwhrc )
1809
1810          DEALLOCATE ( rrtm_sw_taucld )
1811          DEALLOCATE ( rrtm_sw_ssacld )
1812          DEALLOCATE ( rrtm_sw_asmcld )
1813          DEALLOCATE ( rrtm_sw_fsfcld )
1814          DEALLOCATE ( rrtm_sw_tauaer )
1815          DEALLOCATE ( rrtm_sw_ssaaer )
1816          DEALLOCATE ( rrtm_sw_asmaer ) 
1817          DEALLOCATE ( rrtm_sw_ecaer )   
1818 
1819          DEALLOCATE ( rrtm_swdflx  )
1820          DEALLOCATE ( rrtm_swdflxc )
1821          DEALLOCATE ( rrtm_swuflx  )
1822          DEALLOCATE ( rrtm_swuflxc )
1823          DEALLOCATE ( rrtm_swhr  )
1824          DEALLOCATE ( rrtm_swhrc )
1825
1826       ENDIF
1827
1828!
1829!--    Open file for reading
1830       nc_stat = NF90_OPEN( rrtm_input_file, NF90_NOWRITE, id )
1831       CALL netcdf_handle_error_rad( 'read_sounding_data', 549 )
1832
1833!
1834!--    Inquire dimension of z axis and save in nz_snd
1835       nc_stat = NF90_INQ_DIMID( id, "Pressure", id_dim_zrad )
1836       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim_zrad, len = nz_snd )
1837       CALL netcdf_handle_error_rad( 'read_sounding_data', 551 )
1838
1839!
1840! !--    Allocate temporary array for storing pressure data
1841       ALLOCATE( hyp_snd_tmp(1:nz_snd) )
1842       hyp_snd_tmp = 0.0_wp
1843
1844
1845!--    Read pressure from file
1846       nc_stat = NF90_INQ_VARID( id, "Pressure", id_var )
1847       nc_stat = NF90_GET_VAR( id, id_var, hyp_snd_tmp(:), start = (/1/),      &
1848                               count = (/nz_snd/) )
1849       CALL netcdf_handle_error_rad( 'read_sounding_data', 552 )
1850
1851!
1852!--    Allocate temporary array for storing temperature data
1853       ALLOCATE( t_snd_tmp(1:nz_snd) )
1854       t_snd_tmp = 0.0_wp
1855
1856!
1857!--    Read temperature from file
1858       nc_stat = NF90_INQ_VARID( id, "ReferenceTemperature", id_var )
1859       nc_stat = NF90_GET_VAR( id, id_var, t_snd_tmp(:), start = (/1/),        &
1860                               count = (/nz_snd/) )
1861       CALL netcdf_handle_error_rad( 'read_sounding_data', 553 )
1862
1863!
1864!--    Calculate start of sounding data
1865       nz_snd_start = nz_snd + 1
1866       nz_snd_end   = nz_snd + 1
1867
1868!
1869!--    Start filling vertical dimension at 10hPa above the model domain (hyp is
1870!--    in Pa, hyp_snd in hPa).
1871       DO  k = 1, nz_snd
1872          IF ( hyp_snd_tmp(k) < ( hyp(nzt+1) - 1000.0_wp) * 0.01_wp )  THEN
1873             nz_snd_start = k
1874             EXIT
1875          END IF
1876       END DO
1877
1878       IF ( nz_snd_start <= nz_snd )  THEN
1879          nz_snd_end = nz_snd
1880       END IF
1881
1882
1883!
1884!--    Calculate of total grid points for RRTMG calculations
1885       nzt_rad = nzt + nz_snd_end - nz_snd_start + 1
1886
1887!
1888!--    Save data above LES domain in hyp_snd, t_snd and q_snd
1889!--    Note: q_snd_tmp is not calculated at the moment (dry residual atmosphere)
1890       ALLOCATE( hyp_snd(nzb+1:nzt_rad) )
1891       ALLOCATE( t_snd(nzb+1:nzt_rad)   )
1892       ALLOCATE( q_snd(nzb+1:nzt_rad)   )
1893       hyp_snd = 0.0_wp
1894       t_snd = 0.0_wp
1895       q_snd = 0.0_wp
1896
1897       hyp_snd(nzt+2:nzt_rad) = hyp_snd_tmp(nz_snd_start+1:nz_snd_end)
1898       t_snd(nzt+2:nzt_rad)   = t_snd_tmp(nz_snd_start+1:nz_snd_end)
1899
1900       nc_stat = NF90_CLOSE( id )
1901
1902!
1903!--    Calculate pressure levels on zu and zw grid. Sounding data is added at
1904!--    top of the LES domain. This routine does not consider horizontal or
1905!--    vertical variability of pressure and temperature
1906       ALLOCATE ( rrtm_play(0:0,nzb+1:nzt_rad+1)   )
1907       ALLOCATE ( rrtm_plev(0:0,nzb+1:nzt_rad+2)   )
1908
1909       t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
1910       DO k = nzb+1, nzt+1
1911          rrtm_play(0,k) = hyp(k) * 0.01_wp
1912          rrtm_plev(0,k) = surface_pressure * ( (t_surface - g/cp * zw(k-1)) / &
1913                         t_surface )**(1.0_wp/0.286_wp)
1914       ENDDO
1915
1916       DO k = nzt+2, nzt_rad
1917          rrtm_play(0,k) = hyp_snd(k)
1918          rrtm_plev(0,k) = 0.5_wp * ( rrtm_play(0,k) + rrtm_play(0,k-1) )
1919       ENDDO
1920       rrtm_plev(0,nzt_rad+1) = MAX( 0.5 * hyp_snd(nzt_rad),                   &
1921                                   1.5 * hyp_snd(nzt_rad)                      &
1922                                 - 0.5 * hyp_snd(nzt_rad-1) )
1923       rrtm_plev(0,nzt_rad+2)  = MIN( 1.0E-4_wp,                               &
1924                                      0.25_wp * rrtm_plev(0,nzt_rad+1) )
1925
1926       rrtm_play(0,nzt_rad+1) = 0.5 * rrtm_plev(0,nzt_rad+1)
1927
1928!
1929!--    Calculate temperature/humidity levels at top of the LES domain.
1930!--    Currently, the temperature is taken from sounding data (might lead to a
1931!--    temperature jump at interface. To do: Humidity is currently not
1932!--    calculated above the LES domain.
1933       ALLOCATE ( rrtm_tlay(0:0,nzb+1:nzt_rad+1)   )
1934       ALLOCATE ( rrtm_tlev(0:0,nzb+1:nzt_rad+2)   )
1935       ALLOCATE ( rrtm_h2ovmr(0:0,nzb+1:nzt_rad+1) )
1936
1937       DO k = nzt+8, nzt_rad
1938          rrtm_tlay(0,k)   = t_snd(k)
1939          rrtm_h2ovmr(0,k) = q_snd(k)
1940       ENDDO
1941       rrtm_tlay(0,nzt_rad+1) = 2.0_wp * rrtm_tlay(0,nzt_rad)                 &
1942                                - rrtm_tlay(0,nzt_rad-1)
1943       DO k = nzt+9, nzt_rad+1
1944          rrtm_tlev(0,k)   = rrtm_tlay(0,k-1) + (rrtm_tlay(0,k)                &
1945                             - rrtm_tlay(0,k-1))                               &
1946                             / ( rrtm_play(0,k) - rrtm_play(0,k-1) )           &
1947                             * ( rrtm_plev(0,k) - rrtm_play(0,k-1) )
1948       ENDDO
1949       rrtm_h2ovmr(0,nzt_rad+1) = rrtm_h2ovmr(0,nzt_rad)
1950
1951       rrtm_tlev(0,nzt_rad+2)   = 2.0_wp * rrtm_tlay(0,nzt_rad+1)              &
1952                                  - rrtm_tlev(0,nzt_rad)
1953!
1954!--    Allocate remaining RRTMG arrays
1955       ALLOCATE ( rrtm_cicewp(0:0,nzb+1:nzt_rad+1) )
1956       ALLOCATE ( rrtm_cldfr(0:0,nzb+1:nzt_rad+1) )
1957       ALLOCATE ( rrtm_cliqwp(0:0,nzb+1:nzt_rad+1) )
1958       ALLOCATE ( rrtm_reice(0:0,nzb+1:nzt_rad+1) )
1959       ALLOCATE ( rrtm_reliq(0:0,nzb+1:nzt_rad+1) )
1960       ALLOCATE ( rrtm_lw_taucld(1:nbndlw+1,0:0,nzb+1:nzt_rad+1) )
1961       ALLOCATE ( rrtm_lw_tauaer(0:0,nzb+1:nzt_rad+1,1:nbndlw+1) )
1962       ALLOCATE ( rrtm_sw_taucld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
1963       ALLOCATE ( rrtm_sw_ssacld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
1964       ALLOCATE ( rrtm_sw_asmcld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
1965       ALLOCATE ( rrtm_sw_fsfcld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
1966       ALLOCATE ( rrtm_sw_tauaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) )
1967       ALLOCATE ( rrtm_sw_ssaaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) )
1968       ALLOCATE ( rrtm_sw_asmaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) ) 
1969       ALLOCATE ( rrtm_sw_ecaer(0:0,nzb+1:nzt_rad+1,1:naerec+1) )   
1970
1971!
1972!--    The ice phase is currently not considered in PALM
1973       rrtm_cicewp = 0.0_wp
1974       rrtm_reice  = 0.0_wp
1975
1976!
1977!--    Set other parameters (move to NAMELIST parameters in the future)
1978       rrtm_lw_tauaer = 0.0_wp
1979       rrtm_lw_taucld = 0.0_wp
1980       rrtm_sw_taucld = 0.0_wp
1981       rrtm_sw_ssacld = 0.0_wp
1982       rrtm_sw_asmcld = 0.0_wp
1983       rrtm_sw_fsfcld = 0.0_wp
1984       rrtm_sw_tauaer = 0.0_wp
1985       rrtm_sw_ssaaer = 0.0_wp
1986       rrtm_sw_asmaer = 0.0_wp
1987       rrtm_sw_ecaer  = 0.0_wp
1988
1989
1990       ALLOCATE ( rrtm_swdflx(0:0,nzb:nzt_rad+1)  )
1991       ALLOCATE ( rrtm_swuflx(0:0,nzb:nzt_rad+1)  )
1992       ALLOCATE ( rrtm_swhr(0:0,nzb+1:nzt_rad+1)  )
1993       ALLOCATE ( rrtm_swuflxc(0:0,nzb:nzt_rad+1) )
1994       ALLOCATE ( rrtm_swdflxc(0:0,nzb:nzt_rad+1) )
1995       ALLOCATE ( rrtm_swhrc(0:0,nzb+1:nzt_rad+1) )
1996
1997       rrtm_swdflx  = 0.0_wp
1998       rrtm_swuflx  = 0.0_wp
1999       rrtm_swhr    = 0.0_wp 
2000       rrtm_swuflxc = 0.0_wp
2001       rrtm_swdflxc = 0.0_wp
2002       rrtm_swhrc   = 0.0_wp
2003
2004       ALLOCATE ( rrtm_lwdflx(0:0,nzb:nzt_rad+1)  )
2005       ALLOCATE ( rrtm_lwuflx(0:0,nzb:nzt_rad+1)  )
2006       ALLOCATE ( rrtm_lwhr(0:0,nzb+1:nzt_rad+1)  )
2007       ALLOCATE ( rrtm_lwuflxc(0:0,nzb:nzt_rad+1) )
2008       ALLOCATE ( rrtm_lwdflxc(0:0,nzb:nzt_rad+1) )
2009       ALLOCATE ( rrtm_lwhrc(0:0,nzb+1:nzt_rad+1) )
2010
2011       rrtm_lwdflx  = 0.0_wp
2012       rrtm_lwuflx  = 0.0_wp
2013       rrtm_lwhr    = 0.0_wp 
2014       rrtm_lwuflxc = 0.0_wp
2015       rrtm_lwdflxc = 0.0_wp
2016       rrtm_lwhrc   = 0.0_wp
2017
2018       ALLOCATE ( rrtm_lwuflx_dt(0:0,nzb:nzt_rad+1) )
2019       ALLOCATE ( rrtm_lwuflxc_dt(0:0,nzb:nzt_rad+1) )
2020
2021       rrtm_lwuflx_dt = 0.0_wp
2022       rrtm_lwuflxc_dt = 0.0_wp
2023
2024    END SUBROUTINE read_sounding_data
2025
2026
2027!------------------------------------------------------------------------------!
2028! Description:
2029! ------------
2030!> Read trace gas data from file
2031!------------------------------------------------------------------------------!
2032    SUBROUTINE read_trace_gas_data
2033
2034       USE rrsw_ncpar
2035
2036       IMPLICIT NONE
2037
2038       INTEGER(iwp), PARAMETER :: num_trace_gases = 9 !< number of trace gases (absorbers)
2039
2040       CHARACTER(LEN=5), DIMENSION(num_trace_gases), PARAMETER ::              & !< trace gas names
2041           trace_names = (/'O3   ', 'CO2  ', 'CH4  ', 'N2O  ', 'O2   ',        &
2042                           'CFC11', 'CFC12', 'CFC22', 'CCL4 '/)
2043
2044       INTEGER(iwp) :: id,     & !< NetCDF id
2045                       k,      & !< loop index
2046                       m,      & !< loop index
2047                       n,      & !< loop index
2048                       nabs,   & !< number of absorbers
2049                       np,     & !< number of pressure levels
2050                       id_abs, & !< NetCDF id of the respective absorber
2051                       id_dim, & !< NetCDF id of asborber's dimension
2052                       id_var    !< NetCDf id ot the absorber
2053
2054       REAL(wp) :: p_mls_l, p_mls_u, p_wgt_l, p_wgt_u, p_mls_m
2055
2056
2057       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  p_mls,         & !< pressure levels for the absorbers
2058                                                 rrtm_play_tmp, & !< temporary array for pressure zu-levels
2059                                                 rrtm_plev_tmp, & !< temporary array for pressure zw-levels
2060                                                 trace_path_tmp   !< temporary array for storing trace gas path data
2061
2062       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  trace_mls,      & !< array for storing the absorber amounts
2063                                                 trace_mls_path, & !< array for storing trace gas path data
2064                                                 trace_mls_tmp     !< temporary array for storing trace gas data
2065
2066
2067!
2068!--    In case of updates, deallocate arrays first (sufficient to check one
2069!--    array as the others are automatically allocated)
2070       IF ( ALLOCATED ( rrtm_o3vmr ) )  THEN
2071          DEALLOCATE ( rrtm_o3vmr  )
2072          DEALLOCATE ( rrtm_co2vmr )
2073          DEALLOCATE ( rrtm_ch4vmr )
2074          DEALLOCATE ( rrtm_n2ovmr )
2075          DEALLOCATE ( rrtm_o2vmr  )
2076          DEALLOCATE ( rrtm_cfc11vmr )
2077          DEALLOCATE ( rrtm_cfc12vmr )
2078          DEALLOCATE ( rrtm_cfc22vmr )
2079          DEALLOCATE ( rrtm_ccl4vmr  )
2080       ENDIF
2081
2082!
2083!--    Allocate trace gas profiles
2084       ALLOCATE ( rrtm_o3vmr(0:0,1:nzt_rad+1)  )
2085       ALLOCATE ( rrtm_co2vmr(0:0,1:nzt_rad+1) )
2086       ALLOCATE ( rrtm_ch4vmr(0:0,1:nzt_rad+1) )
2087       ALLOCATE ( rrtm_n2ovmr(0:0,1:nzt_rad+1) )
2088       ALLOCATE ( rrtm_o2vmr(0:0,1:nzt_rad+1)  )
2089       ALLOCATE ( rrtm_cfc11vmr(0:0,1:nzt_rad+1) )
2090       ALLOCATE ( rrtm_cfc12vmr(0:0,1:nzt_rad+1) )
2091       ALLOCATE ( rrtm_cfc22vmr(0:0,1:nzt_rad+1) )
2092       ALLOCATE ( rrtm_ccl4vmr(0:0,1:nzt_rad+1)  )
2093
2094!
2095!--    Open file for reading
2096       nc_stat = NF90_OPEN( rrtm_input_file, NF90_NOWRITE, id )
2097       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 549 )
2098!
2099!--    Inquire dimension ids and dimensions
2100       nc_stat = NF90_INQ_DIMID( id, "Pressure", id_dim )
2101       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
2102       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, len = np) 
2103       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
2104
2105       nc_stat = NF90_INQ_DIMID( id, "Absorber", id_dim )
2106       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
2107       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, len = nabs ) 
2108       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
2109   
2110
2111!
2112!--    Allocate pressure, and trace gas arrays     
2113       ALLOCATE( p_mls(1:np) )
2114       ALLOCATE( trace_mls(1:num_trace_gases,1:np) ) 
2115       ALLOCATE( trace_mls_tmp(1:nabs,1:np) ) 
2116
2117
2118       nc_stat = NF90_INQ_VARID( id, "Pressure", id_var )
2119       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
2120       nc_stat = NF90_GET_VAR( id, id_var, p_mls )
2121       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
2122
2123       nc_stat = NF90_INQ_VARID( id, "AbsorberAmountMLS", id_var )
2124       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
2125       nc_stat = NF90_GET_VAR( id, id_var, trace_mls_tmp )
2126       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
2127
2128
2129!
2130!--    Write absorber amounts (mls) to trace_mls
2131       DO n = 1, num_trace_gases
2132          CALL getAbsorberIndex( TRIM( trace_names(n) ), id_abs )
2133
2134          trace_mls(n,1:np) = trace_mls_tmp(id_abs,1:np)
2135
2136!
2137!--       Replace missing values by zero
2138          WHERE ( trace_mls(n,:) > 2.0_wp ) 
2139             trace_mls(n,:) = 0.0_wp
2140          END WHERE
2141       END DO
2142
2143       DEALLOCATE ( trace_mls_tmp )
2144
2145       nc_stat = NF90_CLOSE( id )
2146       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 551 )
2147
2148!
2149!--    Add extra pressure level for calculations of the trace gas paths
2150       ALLOCATE ( rrtm_play_tmp(1:nzt_rad+1) )
2151       ALLOCATE ( rrtm_plev_tmp(1:nzt_rad+2) )
2152
2153       rrtm_play_tmp(1:nzt_rad)   = rrtm_play(0,1:nzt_rad) 
2154       rrtm_plev_tmp(1:nzt_rad+1) = rrtm_plev(0,1:nzt_rad+1)
2155       rrtm_play_tmp(nzt_rad+1)   = rrtm_plev(0,nzt_rad+1) * 0.5_wp
2156       rrtm_plev_tmp(nzt_rad+2)   = MIN( 1.0E-4_wp, 0.25_wp                    &
2157                                         * rrtm_plev(0,nzt_rad+1) )
2158 
2159!
2160!--    Calculate trace gas path (zero at surface) with interpolation to the
2161!--    sounding levels
2162       ALLOCATE ( trace_mls_path(1:nzt_rad+2,1:num_trace_gases) )
2163
2164       trace_mls_path(nzb+1,:) = 0.0_wp
2165       
2166       DO k = nzb+2, nzt_rad+2
2167          DO m = 1, num_trace_gases
2168             trace_mls_path(k,m) = trace_mls_path(k-1,m)
2169
2170!
2171!--          When the pressure level is higher than the trace gas pressure
2172!--          level, assume that
2173             IF ( rrtm_plev_tmp(k-1) > p_mls(1) )  THEN             
2174               
2175                trace_mls_path(k,m) = trace_mls_path(k,m) + trace_mls(m,1)     &
2176                                      * ( rrtm_plev_tmp(k-1)                   &
2177                                          - MAX( p_mls(1), rrtm_plev_tmp(k) )  &
2178                                        ) / g
2179             ENDIF
2180
2181!
2182!--          Integrate for each sounding level from the contributing p_mls
2183!--          levels
2184             DO n = 2, np
2185!
2186!--             Limit p_mls so that it is within the model level
2187                p_mls_u = MIN( rrtm_plev_tmp(k-1),                             &
2188                          MAX( rrtm_plev_tmp(k), p_mls(n) ) )
2189                p_mls_l = MIN( rrtm_plev_tmp(k-1),                             &
2190                          MAX( rrtm_plev_tmp(k), p_mls(n-1) ) )
2191
2192                IF ( p_mls_l > p_mls_u )  THEN
2193
2194!
2195!--                Calculate weights for interpolation
2196                   p_mls_m = 0.5_wp * (p_mls_l + p_mls_u)
2197                   p_wgt_u = (p_mls(n-1) - p_mls_m) / (p_mls(n-1) - p_mls(n))
2198                   p_wgt_l = (p_mls_m - p_mls(n))   / (p_mls(n-1) - p_mls(n))
2199
2200!
2201!--                Add level to trace gas path
2202                   trace_mls_path(k,m) = trace_mls_path(k,m)                   &
2203                                         +  ( p_wgt_u * trace_mls(m,n)         &
2204                                            + p_wgt_l * trace_mls(m,n-1) )     &
2205                                         * (p_mls_l - p_mls_u) / g
2206                ENDIF
2207             ENDDO
2208
2209             IF ( rrtm_plev_tmp(k) < p_mls(np) )  THEN
2210                trace_mls_path(k,m) = trace_mls_path(k,m) + trace_mls(m,np)    &
2211                                      * ( MIN( rrtm_plev_tmp(k-1), p_mls(np) ) &
2212                                          - rrtm_plev_tmp(k)                   &
2213                                        ) / g 
2214             ENDIF 
2215          ENDDO
2216       ENDDO
2217
2218
2219!
2220!--    Prepare trace gas path profiles
2221       ALLOCATE ( trace_path_tmp(1:nzt_rad+1) )
2222
2223       DO m = 1, num_trace_gases
2224
2225          trace_path_tmp(1:nzt_rad+1) = ( trace_mls_path(2:nzt_rad+2,m)        &
2226                                       - trace_mls_path(1:nzt_rad+1,m) ) * g   &
2227                                       / ( rrtm_plev_tmp(1:nzt_rad+1)          &
2228                                       - rrtm_plev_tmp(2:nzt_rad+2) )
2229
2230!
2231!--       Save trace gas paths to the respective arrays
2232          SELECT CASE ( TRIM( trace_names(m) ) )
2233
2234             CASE ( 'O3' )
2235
2236                rrtm_o3vmr(0,:) = trace_path_tmp(:)
2237
2238             CASE ( 'CO2' )
2239
2240                rrtm_co2vmr(0,:) = trace_path_tmp(:)
2241
2242             CASE ( 'CH4' )
2243
2244                rrtm_ch4vmr(0,:) = trace_path_tmp(:)
2245
2246             CASE ( 'N2O' )
2247
2248                rrtm_n2ovmr(0,:) = trace_path_tmp(:)
2249
2250             CASE ( 'O2' )
2251
2252                rrtm_o2vmr(0,:) = trace_path_tmp(:)
2253
2254             CASE ( 'CFC11' )
2255
2256                rrtm_cfc11vmr(0,:) = trace_path_tmp(:)
2257
2258             CASE ( 'CFC12' )
2259
2260                rrtm_cfc12vmr(0,:) = trace_path_tmp(:)
2261
2262             CASE ( 'CFC22' )
2263
2264                rrtm_cfc22vmr(0,:) = trace_path_tmp(:)
2265
2266             CASE ( 'CCL4' )
2267
2268                rrtm_ccl4vmr(0,:) = trace_path_tmp(:)
2269
2270             CASE DEFAULT
2271
2272          END SELECT
2273
2274       ENDDO
2275
2276       DEALLOCATE ( trace_path_tmp )
2277       DEALLOCATE ( trace_mls_path )
2278       DEALLOCATE ( rrtm_play_tmp )
2279       DEALLOCATE ( rrtm_plev_tmp )
2280       DEALLOCATE ( trace_mls )
2281       DEALLOCATE ( p_mls )
2282
2283    END SUBROUTINE read_trace_gas_data
2284
2285
2286    SUBROUTINE netcdf_handle_error_rad( routine_name, errno )
2287
2288       USE control_parameters,                                                 &
2289           ONLY:  message_string
2290
2291       USE NETCDF
2292
2293       USE pegrid
2294
2295       IMPLICIT NONE
2296
2297       CHARACTER(LEN=6) ::  message_identifier
2298       CHARACTER(LEN=*) ::  routine_name
2299
2300       INTEGER(iwp) ::  errno
2301
2302       IF ( nc_stat /= NF90_NOERR )  THEN
2303
2304          WRITE( message_identifier, '(''NC'',I4.4)' )  errno
2305          message_string = TRIM( NF90_STRERROR( nc_stat ) )
2306
2307          CALL message( routine_name, message_identifier, 2, 2, 0, 6, 1 )
2308
2309       ENDIF
2310
2311    END SUBROUTINE netcdf_handle_error_rad
2312#endif
2313
2314
2315!------------------------------------------------------------------------------!
2316! Description:
2317! ------------
2318!> Calculate temperature tendency due to radiative cooling/heating.
2319!> Cache-optimized version.
2320!------------------------------------------------------------------------------!
2321 SUBROUTINE radiation_tendency_ij ( i, j, tend )
2322
2323    USE cloud_parameters,                                                      &
2324        ONLY:  pt_d_t
2325
2326    IMPLICIT NONE
2327
2328    INTEGER(iwp) :: i, j, k !< loop indices
2329
2330    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend !< pt tendency term
2331
2332    IF ( radiation_scheme == 'rrtmg' )  THEN
2333#if defined  ( __rrtmg )
2334!
2335!--    Calculate tendency based on heating rate
2336       DO k = nzb+1, nzt+1
2337          tend(k,j,i) = tend(k,j,i) + (rad_lw_hr(k,j,i) + rad_sw_hr(k,j,i))    &
2338                                         * pt_d_t(k) * d_seconds_hour
2339       ENDDO
2340#endif
2341    ENDIF
2342
2343    END SUBROUTINE radiation_tendency_ij
2344
2345
2346!------------------------------------------------------------------------------!
2347! Description:
2348! ------------
2349!> Calculate temperature tendency due to radiative cooling/heating.
2350!> Vector-optimized version
2351!------------------------------------------------------------------------------!
2352 SUBROUTINE radiation_tendency ( tend )
2353
2354    USE cloud_parameters,                                                      &
2355        ONLY:  pt_d_t
2356
2357    USE indices,                                                               &
2358        ONLY:  nxl, nxr, nyn, nys
2359
2360    IMPLICIT NONE
2361
2362    INTEGER(iwp) :: i, j, k !< loop indices
2363
2364    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend !< pt tendency term
2365
2366    IF ( radiation_scheme == 'rrtmg' )  THEN
2367#if defined  ( __rrtmg )
2368!
2369!--    Calculate tendency based on heating rate
2370       DO  i = nxl, nxr
2371          DO  j = nys, nyn
2372             DO k = nzb+1, nzt+1
2373                tend(k,j,i) = tend(k,j,i) + ( rad_lw_hr(k,j,i)                 &
2374                                          +  rad_sw_hr(k,j,i) ) * pt_d_t(k)    &
2375                                          * d_seconds_hour
2376             ENDDO
2377          ENDDO
2378       ENDDO
2379#endif
2380    ENDIF
2381
2382
2383 END SUBROUTINE radiation_tendency
2384
2385!------------------------------------------------------------------------------!
2386!
2387! Description:
2388! ------------
2389!> Subroutine for averaging 3D data
2390!------------------------------------------------------------------------------!
2391SUBROUTINE radiation_3d_data_averaging( mode, variable )
2392 
2393
2394    USE control_parameters
2395
2396    USE indices
2397
2398    USE kinds
2399
2400    IMPLICIT NONE
2401
2402    CHARACTER (LEN=*) ::  mode    !<
2403    CHARACTER (LEN=*) :: variable !<
2404
2405    INTEGER(iwp) ::  i !<
2406    INTEGER(iwp) ::  j !<
2407    INTEGER(iwp) ::  k !<
2408
2409    IF ( mode == 'allocate' )  THEN
2410
2411       SELECT CASE ( TRIM( variable ) )
2412
2413             CASE ( 'rad_net*' )
2414                IF ( .NOT. ALLOCATED( rad_net_av ) )  THEN
2415                   ALLOCATE( rad_net_av(nysg:nyng,nxlg:nxrg) )
2416                ENDIF
2417                rad_net_av = 0.0_wp
2418
2419             CASE ( 'rad_lw_in' )
2420                IF ( .NOT. ALLOCATED( rad_lw_in_av ) )  THEN
2421                   ALLOCATE( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2422                ENDIF
2423                rad_lw_in_av = 0.0_wp
2424
2425             CASE ( 'rad_lw_out' )
2426                IF ( .NOT. ALLOCATED( rad_lw_out_av ) )  THEN
2427                   ALLOCATE( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2428                ENDIF
2429                rad_lw_out_av = 0.0_wp
2430
2431             CASE ( 'rad_lw_cs_hr' )
2432                IF ( .NOT. ALLOCATED( rad_lw_cs_hr_av ) )  THEN
2433                   ALLOCATE( rad_lw_cs_hr_av(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
2434                ENDIF
2435                rad_lw_cs_hr_av = 0.0_wp
2436
2437             CASE ( 'rad_lw_hr' )
2438                IF ( .NOT. ALLOCATED( rad_lw_hr_av ) )  THEN
2439                   ALLOCATE( rad_lw_hr_av(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
2440                ENDIF
2441                rad_lw_hr_av = 0.0_wp
2442
2443             CASE ( 'rad_sw_in' )
2444                IF ( .NOT. ALLOCATED( rad_sw_in_av ) )  THEN
2445                   ALLOCATE( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2446                ENDIF
2447                rad_sw_in_av = 0.0_wp
2448
2449             CASE ( 'rad_sw_out' )
2450                IF ( .NOT. ALLOCATED( rad_sw_out_av ) )  THEN
2451                   ALLOCATE( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2452                ENDIF
2453                rad_sw_out_av = 0.0_wp
2454
2455             CASE ( 'rad_sw_cs_hr' )
2456                IF ( .NOT. ALLOCATED( rad_sw_cs_hr_av ) )  THEN
2457                   ALLOCATE( rad_sw_cs_hr_av(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
2458                ENDIF
2459                rad_sw_cs_hr_av = 0.0_wp
2460
2461             CASE ( 'rad_sw_hr' )
2462                IF ( .NOT. ALLOCATED( rad_sw_hr_av ) )  THEN
2463                   ALLOCATE( rad_sw_hr_av(nzb+1:nzt+1,nysg:nyng,nxlg:nxrg) )
2464                ENDIF
2465                rad_sw_hr_av = 0.0_wp
2466
2467          CASE DEFAULT
2468             CONTINUE
2469
2470       END SELECT
2471
2472    ELSEIF ( mode == 'sum' )  THEN
2473
2474       SELECT CASE ( TRIM( variable ) )
2475
2476          CASE ( 'rad_net*' )
2477             DO  i = nxlg, nxrg
2478                DO  j = nysg, nyng
2479                   rad_net_av(j,i) = rad_net_av(j,i) + rad_net(j,i)
2480                ENDDO
2481             ENDDO
2482
2483          CASE ( 'rad_lw_in' )
2484             DO  i = nxlg, nxrg
2485                DO  j = nysg, nyng
2486                   DO  k = nzb, nzt+1
2487                      rad_lw_in_av(k,j,i) = rad_lw_in_av(k,j,i) + rad_lw_in(k,j,i)
2488                   ENDDO
2489                ENDDO
2490             ENDDO
2491
2492          CASE ( 'rad_lw_out' )
2493             DO  i = nxlg, nxrg
2494                DO  j = nysg, nyng
2495                   DO  k = nzb, nzt+1
2496                      rad_lw_out_av(k,j,i) = rad_lw_out_av(k,j,i) + rad_lw_out(k,j,i)
2497                   ENDDO
2498                ENDDO
2499             ENDDO
2500
2501          CASE ( 'rad_lw_cs_hr' )
2502             DO  i = nxlg, nxrg
2503                DO  j = nysg, nyng
2504                   DO  k = nzb, nzt+1
2505                      rad_lw_cs_hr_av(k,j,i) = rad_lw_cs_hr_av(k,j,i) + rad_lw_cs_hr(k,j,i)
2506                   ENDDO
2507                ENDDO
2508             ENDDO
2509
2510          CASE ( 'rad_lw_hr' )
2511             DO  i = nxlg, nxrg
2512                DO  j = nysg, nyng
2513                   DO  k = nzb, nzt+1
2514                      rad_lw_hr_av(k,j,i) = rad_lw_hr_av(k,j,i) + rad_lw_hr(k,j,i)
2515                   ENDDO
2516                ENDDO
2517             ENDDO
2518
2519          CASE ( 'rad_sw_in' )
2520             DO  i = nxlg, nxrg
2521                DO  j = nysg, nyng
2522                   DO  k = nzb, nzt+1
2523                      rad_sw_in_av(k,j,i) = rad_sw_in_av(k,j,i) + rad_sw_in(k,j,i)
2524                   ENDDO
2525                ENDDO
2526             ENDDO
2527
2528          CASE ( 'rad_sw_out' )
2529             DO  i = nxlg, nxrg
2530                DO  j = nysg, nyng
2531                   DO  k = nzb, nzt+1
2532                      rad_sw_out_av(k,j,i) = rad_sw_out_av(k,j,i) + rad_sw_out(k,j,i)
2533                   ENDDO
2534                ENDDO
2535             ENDDO
2536
2537          CASE ( 'rad_sw_cs_hr' )
2538             DO  i = nxlg, nxrg
2539                DO  j = nysg, nyng
2540                   DO  k = nzb, nzt+1
2541                      rad_sw_cs_hr_av(k,j,i) = rad_sw_cs_hr_av(k,j,i) + rad_sw_cs_hr(k,j,i)
2542                   ENDDO
2543                ENDDO
2544             ENDDO
2545
2546          CASE ( 'rad_sw_hr' )
2547             DO  i = nxlg, nxrg
2548                DO  j = nysg, nyng
2549                   DO  k = nzb, nzt+1
2550                      rad_sw_hr_av(k,j,i) = rad_sw_hr_av(k,j,i) + rad_sw_hr(k,j,i)
2551                   ENDDO
2552                ENDDO
2553             ENDDO
2554
2555          CASE DEFAULT
2556             CONTINUE
2557
2558       END SELECT
2559
2560    ELSEIF ( mode == 'average' )  THEN
2561
2562       SELECT CASE ( TRIM( variable ) )
2563
2564         CASE ( 'rad_net*' )
2565             DO  i = nxlg, nxrg
2566                DO  j = nysg, nyng
2567                   rad_net_av(j,i) = rad_net_av(j,i) / REAL( average_count_3d, KIND=wp )
2568                ENDDO
2569             ENDDO
2570
2571          CASE ( 'rad_lw_in' )
2572             DO  i = nxlg, nxrg
2573                DO  j = nysg, nyng
2574                   DO  k = nzb, nzt+1
2575                      rad_lw_in_av(k,j,i) = rad_lw_in_av(k,j,i) / REAL( average_count_3d, KIND=wp )
2576                   ENDDO
2577                ENDDO
2578             ENDDO
2579
2580          CASE ( 'rad_lw_out' )
2581             DO  i = nxlg, nxrg
2582                DO  j = nysg, nyng
2583                   DO  k = nzb, nzt+1
2584                      rad_lw_out_av(k,j,i) = rad_lw_out_av(k,j,i) / REAL( average_count_3d, KIND=wp )
2585                   ENDDO
2586                ENDDO
2587             ENDDO
2588
2589          CASE ( 'rad_lw_cs_hr' )
2590             DO  i = nxlg, nxrg
2591                DO  j = nysg, nyng
2592                   DO  k = nzb, nzt+1
2593                      rad_lw_cs_hr_av(k,j,i) = rad_lw_cs_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
2594                   ENDDO
2595                ENDDO
2596             ENDDO
2597
2598          CASE ( 'rad_lw_hr' )
2599             DO  i = nxlg, nxrg
2600                DO  j = nysg, nyng
2601                   DO  k = nzb, nzt+1
2602                      rad_lw_hr_av(k,j,i) = rad_lw_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
2603                   ENDDO
2604                ENDDO
2605             ENDDO
2606
2607          CASE ( 'rad_sw_in' )
2608             DO  i = nxlg, nxrg
2609                DO  j = nysg, nyng
2610                   DO  k = nzb, nzt+1
2611                      rad_sw_in_av(k,j,i) = rad_sw_in_av(k,j,i) / REAL( average_count_3d, KIND=wp )
2612                   ENDDO
2613                ENDDO
2614             ENDDO
2615
2616          CASE ( 'rad_sw_out' )
2617             DO  i = nxlg, nxrg
2618                DO  j = nysg, nyng
2619                   DO  k = nzb, nzt+1
2620                      rad_sw_out_av(k,j,i) = rad_sw_out_av(k,j,i) / REAL( average_count_3d, KIND=wp )
2621                   ENDDO
2622                ENDDO
2623             ENDDO
2624
2625          CASE ( 'rad_sw_cs_hr' )
2626             DO  i = nxlg, nxrg
2627                DO  j = nysg, nyng
2628                   DO  k = nzb, nzt+1
2629                      rad_sw_cs_hr_av(k,j,i) = rad_sw_cs_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
2630                   ENDDO
2631                ENDDO
2632             ENDDO
2633
2634          CASE ( 'rad_sw_hr' )
2635             DO  i = nxlg, nxrg
2636                DO  j = nysg, nyng
2637                   DO  k = nzb, nzt+1
2638                      rad_sw_hr_av(k,j,i) = rad_sw_hr_av(k,j,i) / REAL( average_count_3d, KIND=wp )
2639                   ENDDO
2640                ENDDO
2641             ENDDO
2642
2643       END SELECT
2644
2645    ENDIF
2646
2647END SUBROUTINE radiation_3d_data_averaging
2648
2649
2650!------------------------------------------------------------------------------!
2651!
2652! Description:
2653! ------------
2654!> Subroutine defining appropriate grid for netcdf variables.
2655!> It is called out from subroutine netcdf.
2656!------------------------------------------------------------------------------!
2657SUBROUTINE radiation_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
2658   
2659    IMPLICIT NONE
2660
2661    CHARACTER (LEN=*), INTENT(IN)  ::  var         !<
2662    LOGICAL, INTENT(OUT)           ::  found       !<
2663    CHARACTER (LEN=*), INTENT(OUT) ::  grid_x      !<
2664    CHARACTER (LEN=*), INTENT(OUT) ::  grid_y      !<
2665    CHARACTER (LEN=*), INTENT(OUT) ::  grid_z      !<
2666
2667    found  = .TRUE.
2668
2669
2670!
2671!-- Check for the grid
2672    SELECT CASE ( TRIM( var ) )
2673
2674       CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_sw_cs_hr', 'rad_sw_hr',        &
2675              'rad_lw_cs_hr_xy', 'rad_lw_hr_xy', 'rad_sw_cs_hr_xy',            &
2676              'rad_sw_hr_xy', 'rad_lw_cs_hr_xz', 'rad_lw_hr_xz',               &
2677              'rad_sw_cs_hr_xz', 'rad_sw_hr_xz', 'rad_lw_cs_hr_yz',            &
2678              'rad_lw_hr_yz', 'rad_sw_cs_hr_yz', 'rad_sw_hr_yz' )
2679          grid_x = 'x'
2680          grid_y = 'y'
2681          grid_z = 'zu'
2682
2683       CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_sw_in', 'rad_sw_out',            &
2684              'rad_lw_in_xy', 'rad_lw_out_xy', 'rad_sw_in_xy','rad_sw_out_xy', &
2685              'rad_lw_in_xz', 'rad_lw_out_xz', 'rad_sw_in_xz','rad_sw_out_xz', &
2686              'rad_lw_in_yz', 'rad_lw_out_yz', 'rad_sw_in_yz','rad_sw_out_yz' )
2687          grid_x = 'x'
2688          grid_y = 'y'
2689          grid_z = 'zw'
2690
2691
2692       CASE DEFAULT
2693          found  = .FALSE.
2694          grid_x = 'none'
2695          grid_y = 'none'
2696          grid_z = 'none'
2697
2698        END SELECT
2699
2700    END SUBROUTINE radiation_define_netcdf_grid
2701
2702!------------------------------------------------------------------------------!
2703!
2704! Description:
2705! ------------
2706!> Subroutine defining 3D output variables
2707!------------------------------------------------------------------------------!
2708 SUBROUTINE radiation_data_output_2d( av, variable, found, grid, mode,         &
2709                                      local_pf, two_d )
2710 
2711    USE indices
2712
2713    USE kinds
2714
2715
2716    IMPLICIT NONE
2717
2718    CHARACTER (LEN=*) ::  grid     !<
2719    CHARACTER (LEN=*) ::  mode     !<
2720    CHARACTER (LEN=*) ::  variable !<
2721
2722    INTEGER(iwp) ::  av !<
2723    INTEGER(iwp) ::  i  !<
2724    INTEGER(iwp) ::  j  !<
2725    INTEGER(iwp) ::  k  !<
2726
2727    LOGICAL      ::  found !<
2728    LOGICAL      ::  two_d !< flag parameter that indicates 2D variables (horizontal cross sections)
2729
2730    REAL(wp), DIMENSION(nxlg:nxrg,nysg:nyng,nzb:nzt+1) ::  local_pf !<
2731
2732    found = .TRUE.
2733
2734    SELECT CASE ( TRIM( variable ) )
2735
2736       CASE ( 'rad_net*_xy' )        ! 2d-array
2737          IF ( av == 0 ) THEN
2738             DO  i = nxlg, nxrg
2739                DO  j = nysg, nyng
2740                   local_pf(i,j,nzb+1) = rad_net(j,i)
2741                ENDDO
2742             ENDDO
2743          ELSE
2744             DO  i = nxlg, nxrg
2745                DO  j = nysg, nyng 
2746                   local_pf(i,j,nzb+1) = rad_net_av(j,i)
2747                ENDDO
2748             ENDDO
2749          ENDIF
2750          two_d = .TRUE.
2751          grid = 'zu1'
2752
2753 
2754       CASE ( 'rad_lw_in_xy', 'rad_lw_in_xz', 'rad_lw_in_yz' )
2755          IF ( av == 0 ) THEN
2756             DO  i = nxlg, nxrg
2757                DO  j = nysg, nyng
2758                   DO  k = nzb, nzt+1
2759                      local_pf(i,j,k) = rad_lw_in(k,j,i)
2760                   ENDDO
2761                ENDDO
2762             ENDDO
2763          ELSE
2764             DO  i = nxlg, nxrg
2765                DO  j = nysg, nyng 
2766                   DO  k = nzb, nzt+1
2767                      local_pf(i,j,k) = rad_lw_in_av(k,j,i)
2768                   ENDDO
2769                ENDDO
2770             ENDDO
2771          ENDIF
2772          IF ( mode == 'xy' )  grid = 'zu'
2773
2774       CASE ( 'rad_lw_out_xy', 'rad_lw_out_xz', 'rad_lw_out_yz' )
2775          IF ( av == 0 ) THEN
2776             DO  i = nxlg, nxrg
2777                DO  j = nysg, nyng
2778                   DO  k = nzb, nzt+1
2779                      local_pf(i,j,k) = rad_lw_out(k,j,i)
2780                   ENDDO
2781                ENDDO
2782             ENDDO
2783          ELSE
2784             DO  i = nxlg, nxrg
2785                DO  j = nysg, nyng 
2786                   DO  k = nzb, nzt+1
2787                      local_pf(i,j,k) = rad_lw_out_av(k,j,i)
2788                   ENDDO
2789                ENDDO
2790             ENDDO
2791          ENDIF   
2792          IF ( mode == 'xy' )  grid = 'zu'
2793
2794       CASE ( 'rad_lw_cs_hr_xy', 'rad_lw_cs_hr_xz', 'rad_lw_cs_hr_yz' )
2795          IF ( av == 0 ) THEN
2796             DO  i = nxlg, nxrg
2797                DO  j = nysg, nyng
2798                   DO  k = nzb, nzt+1
2799                      local_pf(i,j,k) = rad_lw_cs_hr(k,j,i)
2800                   ENDDO
2801                ENDDO
2802             ENDDO
2803          ELSE
2804             DO  i = nxlg, nxrg
2805                DO  j = nysg, nyng 
2806                   DO  k = nzb, nzt+1
2807                      local_pf(i,j,k) = rad_lw_cs_hr_av(k,j,i)
2808                   ENDDO
2809                ENDDO
2810             ENDDO
2811          ENDIF
2812          IF ( mode == 'xy' )  grid = 'zw'
2813
2814       CASE ( 'rad_lw_hr_xy', 'rad_lw_hr_xz', 'rad_lw_hr_yz' )
2815          IF ( av == 0 ) THEN
2816             DO  i = nxlg, nxrg
2817                DO  j = nysg, nyng
2818                   DO  k = nzb, nzt+1
2819                      local_pf(i,j,k) = rad_lw_hr(k,j,i)
2820                   ENDDO
2821                ENDDO
2822             ENDDO
2823          ELSE
2824             DO  i = nxlg, nxrg
2825                DO  j = nysg, nyng 
2826                   DO  k = nzb, nzt+1
2827                      local_pf(i,j,k) = rad_lw_hr_av(k,j,i)
2828                   ENDDO
2829                ENDDO
2830             ENDDO
2831          ENDIF
2832          IF ( mode == 'xy' )  grid = 'zw'
2833
2834       CASE ( 'rad_sw_in_xy', 'rad_sw_in_xz', 'rad_sw_in_yz' )
2835          IF ( av == 0 ) THEN
2836             DO  i = nxlg, nxrg
2837                DO  j = nysg, nyng
2838                   DO  k = nzb, nzt+1
2839                      local_pf(i,j,k) = rad_sw_in(k,j,i)
2840                   ENDDO
2841                ENDDO
2842             ENDDO
2843          ELSE
2844             DO  i = nxlg, nxrg
2845                DO  j = nysg, nyng 
2846                   DO  k = nzb, nzt+1
2847                      local_pf(i,j,k) = rad_sw_in_av(k,j,i)
2848                   ENDDO
2849                ENDDO
2850             ENDDO
2851          ENDIF
2852          IF ( mode == 'xy' )  grid = 'zu'
2853
2854       CASE ( 'rad_sw_out_xy', 'rad_sw_out_xz', 'rad_sw_out_yz' )
2855          IF ( av == 0 ) THEN
2856             DO  i = nxlg, nxrg
2857                DO  j = nysg, nyng
2858                   DO  k = nzb, nzt+1
2859                      local_pf(i,j,k) = rad_sw_out(k,j,i)
2860                   ENDDO
2861                ENDDO
2862             ENDDO
2863          ELSE
2864             DO  i = nxlg, nxrg
2865                DO  j = nysg, nyng 
2866                   DO  k = nzb, nzt+1
2867                      local_pf(i,j,k) = rad_sw_out_av(k,j,i)
2868                   ENDDO
2869                ENDDO
2870             ENDDO
2871          ENDIF
2872          IF ( mode == 'xy' )  grid = 'zu'
2873
2874       CASE ( 'rad_sw_cs_hr_xy', 'rad_sw_cs_hr_xz', 'rad_sw_cs_hr_yz' )
2875          IF ( av == 0 ) THEN
2876             DO  i = nxlg, nxrg
2877                DO  j = nysg, nyng
2878                   DO  k = nzb, nzt+1
2879                      local_pf(i,j,k) = rad_sw_cs_hr(k,j,i)
2880                   ENDDO
2881                ENDDO
2882             ENDDO
2883          ELSE
2884             DO  i = nxlg, nxrg
2885                DO  j = nysg, nyng 
2886                   DO  k = nzb, nzt+1
2887                      local_pf(i,j,k) = rad_sw_cs_hr_av(k,j,i)
2888                   ENDDO
2889                ENDDO
2890             ENDDO
2891          ENDIF
2892          IF ( mode == 'xy' )  grid = 'zw'
2893
2894       CASE ( 'rad_sw_hr_xy', 'rad_sw_hr_xz', 'rad_sw_hr_yz' )
2895          IF ( av == 0 ) THEN
2896             DO  i = nxlg, nxrg
2897                DO  j = nysg, nyng
2898                   DO  k = nzb, nzt+1
2899                      local_pf(i,j,k) = rad_sw_hr(k,j,i)
2900                   ENDDO
2901                ENDDO
2902             ENDDO
2903          ELSE
2904             DO  i = nxlg, nxrg
2905                DO  j = nysg, nyng 
2906                   DO  k = nzb, nzt+1
2907                      local_pf(i,j,k) = rad_sw_hr_av(k,j,i)
2908                   ENDDO
2909                ENDDO
2910             ENDDO
2911          ENDIF
2912          IF ( mode == 'xy' )  grid = 'zw'
2913
2914       CASE DEFAULT
2915          found = .FALSE.
2916          grid  = 'none'
2917
2918    END SELECT
2919 
2920 END SUBROUTINE radiation_data_output_2d
2921
2922
2923!------------------------------------------------------------------------------!
2924!
2925! Description:
2926! ------------
2927!> Subroutine defining 3D output variables
2928!------------------------------------------------------------------------------!
2929 SUBROUTINE radiation_data_output_3d( av, variable, found, local_pf )
2930 
2931
2932    USE indices
2933
2934    USE kinds
2935
2936
2937    IMPLICIT NONE
2938
2939    CHARACTER (LEN=*) ::  variable !<
2940
2941    INTEGER(iwp) ::  av    !<
2942    INTEGER(iwp) ::  i     !<
2943    INTEGER(iwp) ::  j     !<
2944    INTEGER(iwp) ::  k     !<
2945
2946    LOGICAL      ::  found !<
2947
2948    REAL(sp), DIMENSION(nxlg:nxrg,nysg:nyng,nzb:nzt+1) ::  local_pf !<
2949
2950
2951    found = .TRUE.
2952
2953
2954    SELECT CASE ( TRIM( variable ) )
2955
2956      CASE ( 'rad_sw_in' )
2957         IF ( av == 0 )  THEN
2958            DO  i = nxlg, nxrg
2959               DO  j = nysg, nyng
2960                  DO  k = nzb, nzt+1
2961                     local_pf(i,j,k) = rad_sw_in(k,j,i)
2962                  ENDDO
2963               ENDDO
2964            ENDDO
2965         ELSE
2966            DO  i = nxlg, nxrg
2967               DO  j = nysg, nyng
2968                  DO  k = nzb, nzt+1
2969                     local_pf(i,j,k) = rad_sw_in_av(k,j,i)
2970                  ENDDO
2971               ENDDO
2972            ENDDO
2973         ENDIF
2974
2975      CASE ( 'rad_sw_out' )
2976         IF ( av == 0 )  THEN
2977            DO  i = nxlg, nxrg
2978               DO  j = nysg, nyng
2979                  DO  k = nzb, nzt+1
2980                     local_pf(i,j,k) = rad_sw_out(k,j,i)
2981                  ENDDO
2982               ENDDO
2983            ENDDO
2984         ELSE
2985            DO  i = nxlg, nxrg
2986               DO  j = nysg, nyng
2987                  DO  k = nzb, nzt+1
2988                     local_pf(i,j,k) = rad_sw_out_av(k,j,i)
2989                  ENDDO
2990               ENDDO
2991            ENDDO
2992         ENDIF
2993
2994      CASE ( 'rad_sw_cs_hr' )
2995         IF ( av == 0 )  THEN
2996            DO  i = nxlg, nxrg
2997               DO  j = nysg, nyng
2998                  DO  k = nzb, nzt+1
2999                     local_pf(i,j,k) = rad_sw_cs_hr(k,j,i)
3000                  ENDDO
3001               ENDDO
3002            ENDDO
3003         ELSE
3004            DO  i = nxlg, nxrg
3005               DO  j = nysg, nyng
3006                  DO  k = nzb, nzt+1
3007                     local_pf(i,j,k) = rad_sw_cs_hr_av(k,j,i)
3008                  ENDDO
3009               ENDDO
3010            ENDDO
3011         ENDIF
3012
3013      CASE ( 'rad_sw_hr' )
3014         IF ( av == 0 )  THEN
3015            DO  i = nxlg, nxrg
3016               DO  j = nysg, nyng
3017                  DO  k = nzb, nzt+1
3018                     local_pf(i,j,k) = rad_sw_hr(k,j,i)
3019                  ENDDO
3020               ENDDO
3021            ENDDO
3022         ELSE
3023            DO  i = nxlg, nxrg
3024               DO  j = nysg, nyng
3025                  DO  k = nzb, nzt+1
3026                     local_pf(i,j,k) = rad_sw_hr_av(k,j,i)
3027                  ENDDO
3028               ENDDO
3029            ENDDO
3030         ENDIF
3031
3032      CASE ( 'rad_lw_in' )
3033         IF ( av == 0 )  THEN
3034            DO  i = nxlg, nxrg
3035               DO  j = nysg, nyng
3036                  DO  k = nzb, nzt+1
3037                     local_pf(i,j,k) = rad_lw_in(k,j,i)
3038                  ENDDO
3039               ENDDO
3040            ENDDO
3041         ELSE
3042            DO  i = nxlg, nxrg
3043               DO  j = nysg, nyng
3044                  DO  k = nzb, nzt+1
3045                     local_pf(i,j,k) = rad_lw_in_av(k,j,i)
3046                  ENDDO
3047               ENDDO
3048            ENDDO
3049         ENDIF
3050
3051      CASE ( 'rad_lw_out' )
3052         IF ( av == 0 )  THEN
3053            DO  i = nxlg, nxrg
3054               DO  j = nysg, nyng
3055                  DO  k = nzb, nzt+1
3056                     local_pf(i,j,k) = rad_lw_out(k,j,i)
3057                  ENDDO
3058               ENDDO
3059            ENDDO
3060         ELSE
3061            DO  i = nxlg, nxrg
3062               DO  j = nysg, nyng
3063                  DO  k = nzb, nzt+1
3064                     local_pf(i,j,k) = rad_lw_out_av(k,j,i)
3065                  ENDDO
3066               ENDDO
3067            ENDDO
3068         ENDIF
3069
3070      CASE ( 'rad_lw_cs_hr' )
3071         IF ( av == 0 )  THEN
3072            DO  i = nxlg, nxrg
3073               DO  j = nysg, nyng
3074                  DO  k = nzb, nzt+1
3075                     local_pf(i,j,k) = rad_lw_cs_hr(k,j,i)
3076                  ENDDO
3077               ENDDO
3078            ENDDO
3079         ELSE
3080            DO  i = nxlg, nxrg
3081               DO  j = nysg, nyng
3082                  DO  k = nzb, nzt+1
3083                     local_pf(i,j,k) = rad_lw_cs_hr_av(k,j,i)
3084                  ENDDO
3085               ENDDO
3086            ENDDO
3087         ENDIF
3088
3089      CASE ( 'rad_lw_hr' )
3090         IF ( av == 0 )  THEN
3091            DO  i = nxlg, nxrg
3092               DO  j = nysg, nyng
3093                  DO  k = nzb, nzt+1
3094                     local_pf(i,j,k) = rad_lw_hr(k,j,i)
3095                  ENDDO
3096               ENDDO
3097            ENDDO
3098         ELSE
3099            DO  i = nxlg, nxrg
3100               DO  j = nysg, nyng
3101                  DO  k = nzb, nzt+1
3102                     local_pf(i,j,k) = rad_lw_hr_av(k,j,i)
3103                  ENDDO
3104               ENDDO
3105            ENDDO
3106         ENDIF
3107
3108       CASE DEFAULT
3109          found = .FALSE.
3110
3111    END SELECT
3112
3113
3114 END SUBROUTINE radiation_data_output_3d
3115
3116!------------------------------------------------------------------------------!
3117!
3118! Description:
3119! ------------
3120!> Subroutine defining masked data output
3121!------------------------------------------------------------------------------!
3122 SUBROUTINE radiation_data_output_mask( av, variable, found, local_pf )
3123 
3124    USE control_parameters
3125       
3126    USE indices
3127   
3128    USE kinds
3129   
3130
3131    IMPLICIT NONE
3132
3133    CHARACTER (LEN=*) ::  variable   !<
3134
3135    INTEGER(iwp) ::  av   !<
3136    INTEGER(iwp) ::  i    !<
3137    INTEGER(iwp) ::  j    !<
3138    INTEGER(iwp) ::  k    !<
3139
3140    LOGICAL ::  found     !<
3141
3142    REAL(wp),                                                                  &
3143       DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) ::  &
3144          local_pf   !<
3145
3146
3147    found = .TRUE.
3148
3149    SELECT CASE ( TRIM( variable ) )
3150
3151
3152       CASE ( 'rad_lw_in' )
3153          IF ( av == 0 )  THEN
3154             DO  i = 1, mask_size_l(mid,1)
3155                DO  j = 1, mask_size_l(mid,2)
3156                   DO  k = 1, mask_size_l(mid,3)
3157                       local_pf(i,j,k) = rad_lw_in(mask_k(mid,k),              &
3158                                            mask_j(mid,j),mask_i(mid,i))
3159                    ENDDO
3160                 ENDDO
3161              ENDDO
3162          ELSE
3163             DO  i = 1, mask_size_l(mid,1)
3164                DO  j = 1, mask_size_l(mid,2)
3165                   DO  k = 1, mask_size_l(mid,3)
3166                       local_pf(i,j,k) = rad_lw_in_av(mask_k(mid,k),           &
3167                                               mask_j(mid,j),mask_i(mid,i))
3168                   ENDDO
3169                ENDDO
3170             ENDDO
3171          ENDIF
3172
3173       CASE ( 'rad_lw_out' )
3174          IF ( av == 0 )  THEN
3175             DO  i = 1, mask_size_l(mid,1)
3176                DO  j = 1, mask_size_l(mid,2)
3177                   DO  k = 1, mask_size_l(mid,3)
3178                       local_pf(i,j,k) = rad_lw_out(mask_k(mid,k),             &
3179                                            mask_j(mid,j),mask_i(mid,i))
3180                    ENDDO
3181                 ENDDO
3182              ENDDO
3183          ELSE
3184             DO  i = 1, mask_size_l(mid,1)
3185                DO  j = 1, mask_size_l(mid,2)
3186                   DO  k = 1, mask_size_l(mid,3)
3187                       local_pf(i,j,k) = rad_lw_out_av(mask_k(mid,k),          &
3188                                               mask_j(mid,j),mask_i(mid,i))
3189                   ENDDO
3190                ENDDO
3191             ENDDO
3192          ENDIF
3193
3194       CASE ( 'rad_lw_cs_hr' )
3195          IF ( av == 0 )  THEN
3196             DO  i = 1, mask_size_l(mid,1)
3197                DO  j = 1, mask_size_l(mid,2)
3198                   DO  k = 1, mask_size_l(mid,3)
3199                       local_pf(i,j,k) = rad_lw_cs_hr(mask_k(mid,k),           &
3200                                            mask_j(mid,j),mask_i(mid,i))
3201                    ENDDO
3202                 ENDDO
3203              ENDDO
3204          ELSE
3205             DO  i = 1, mask_size_l(mid,1)
3206                DO  j = 1, mask_size_l(mid,2)
3207                   DO  k = 1, mask_size_l(mid,3)
3208                       local_pf(i,j,k) = rad_lw_cs_hr_av(mask_k(mid,k),        &
3209                                               mask_j(mid,j),mask_i(mid,i))
3210                   ENDDO
3211                ENDDO
3212             ENDDO
3213          ENDIF
3214
3215       CASE ( 'rad_lw_hr' )
3216          IF ( av == 0 )  THEN
3217             DO  i = 1, mask_size_l(mid,1)
3218                DO  j = 1, mask_size_l(mid,2)
3219                   DO  k = 1, mask_size_l(mid,3)
3220                       local_pf(i,j,k) = rad_lw_hr(mask_k(mid,k),              &
3221                                            mask_j(mid,j),mask_i(mid,i))
3222                    ENDDO
3223                 ENDDO
3224              ENDDO
3225          ELSE
3226             DO  i = 1, mask_size_l(mid,1)
3227                DO  j = 1, mask_size_l(mid,2)
3228                   DO  k = 1, mask_size_l(mid,3)
3229                       local_pf(i,j,k) = rad_lw_hr_av(mask_k(mid,k),           &
3230                                               mask_j(mid,j),mask_i(mid,i))
3231                   ENDDO
3232                ENDDO
3233             ENDDO
3234          ENDIF
3235
3236       CASE ( 'rad_sw_in' )
3237          IF ( av == 0 )  THEN
3238             DO  i = 1, mask_size_l(mid,1)
3239                DO  j = 1, mask_size_l(mid,2)
3240                   DO  k = 1, mask_size_l(mid,3)
3241                       local_pf(i,j,k) = rad_sw_in(mask_k(mid,k),              &
3242                                            mask_j(mid,j),mask_i(mid,i))
3243                    ENDDO
3244                 ENDDO
3245              ENDDO
3246          ELSE
3247             DO  i = 1, mask_size_l(mid,1)
3248                DO  j = 1, mask_size_l(mid,2)
3249                   DO  k = 1, mask_size_l(mid,3)
3250                       local_pf(i,j,k) = rad_sw_in_av(mask_k(mid,k),           &
3251                                               mask_j(mid,j),mask_i(mid,i))
3252                   ENDDO
3253                ENDDO
3254             ENDDO
3255          ENDIF
3256
3257       CASE ( 'rad_sw_out' )
3258          IF ( av == 0 )  THEN
3259             DO  i = 1, mask_size_l(mid,1)
3260                DO  j = 1, mask_size_l(mid,2)
3261                   DO  k = 1, mask_size_l(mid,3)
3262                       local_pf(i,j,k) = rad_sw_out(mask_k(mid,k),             &
3263                                            mask_j(mid,j),mask_i(mid,i))
3264                    ENDDO
3265                 ENDDO
3266              ENDDO
3267          ELSE
3268             DO  i = 1, mask_size_l(mid,1)
3269                DO  j = 1, mask_size_l(mid,2)
3270                   DO  k = 1, mask_size_l(mid,3)
3271                       local_pf(i,j,k) = rad_sw_out_av(mask_k(mid,k),          &
3272                                               mask_j(mid,j),mask_i(mid,i))
3273                   ENDDO
3274                ENDDO
3275             ENDDO
3276          ENDIF
3277
3278       CASE ( 'rad_sw_cs_hr' )
3279          IF ( av == 0 )  THEN
3280             DO  i = 1, mask_size_l(mid,1)
3281                DO  j = 1, mask_size_l(mid,2)
3282                   DO  k = 1, mask_size_l(mid,3)
3283                       local_pf(i,j,k) = rad_sw_cs_hr(mask_k(mid,k),           &
3284                                            mask_j(mid,j),mask_i(mid,i))
3285                    ENDDO
3286                 ENDDO
3287              ENDDO
3288          ELSE
3289             DO  i = 1, mask_size_l(mid,1)
3290                DO  j = 1, mask_size_l(mid,2)
3291                   DO  k = 1, mask_size_l(mid,3)
3292                       local_pf(i,j,k) = rad_sw_cs_hr_av(mask_k(mid,k),        &
3293                                               mask_j(mid,j),mask_i(mid,i))
3294                   ENDDO
3295                ENDDO
3296             ENDDO
3297          ENDIF
3298
3299       CASE ( 'rad_sw_hr' )
3300          IF ( av == 0 )  THEN
3301             DO  i = 1, mask_size_l(mid,1)
3302                DO  j = 1, mask_size_l(mid,2)
3303                   DO  k = 1, mask_size_l(mid,3)
3304                       local_pf(i,j,k) = rad_sw_hr(mask_k(mid,k),              &
3305                                            mask_j(mid,j),mask_i(mid,i))
3306                    ENDDO
3307                 ENDDO
3308              ENDDO
3309          ELSE
3310             DO  i = 1, mask_size_l(mid,1)
3311                DO  j = 1, mask_size_l(mid,2)
3312                   DO  k = 1, mask_size_l(mid,3)
3313                       local_pf(i,j,k) = rad_sw_hr_av(mask_k(mid,k),           &
3314                                               mask_j(mid,j),mask_i(mid,i))
3315                   ENDDO
3316                ENDDO
3317             ENDDO
3318          ENDIF
3319
3320       CASE DEFAULT
3321          found = .FALSE.
3322
3323    END SELECT
3324
3325
3326 END SUBROUTINE radiation_data_output_mask
3327
3328
3329!------------------------------------------------------------------------------!
3330!
3331! Description:
3332! ------------
3333!> Subroutine defines masked output variables
3334!------------------------------------------------------------------------------!
3335 SUBROUTINE radiation_last_actions
3336 
3337
3338    USE control_parameters
3339       
3340    USE kinds
3341
3342    IMPLICIT NONE
3343
3344    IF ( write_binary(1:4) == 'true' )  THEN
3345       IF ( ALLOCATED( rad_net ) )  THEN
3346          WRITE ( 14 )  'rad_net             ';  WRITE ( 14 )  rad_net 
3347       ENDIF
3348       IF ( ALLOCATED( rad_net_av ) )  THEN
3349          WRITE ( 14 )  'rad_net_av          ';  WRITE ( 14 )  rad_net_av 
3350       ENDIF 
3351       IF ( ALLOCATED( rad_lw_in ) )  THEN
3352          WRITE ( 14 )  'rad_lw_in           ';  WRITE ( 14 )  rad_lw_in 
3353       ENDIF
3354       IF ( ALLOCATED( rad_lw_in_av ) )  THEN
3355          WRITE ( 14 )  'rad_lw_in_av        ';  WRITE ( 14 )  rad_lw_in_av 
3356       ENDIF
3357       IF ( ALLOCATED( rad_lw_out ) )  THEN
3358          WRITE ( 14 )  'rad_lw_out          ';  WRITE ( 14 )  rad_lw_out 
3359       ENDIF
3360       IF ( ALLOCATED( rad_lw_out_av ) )  THEN
3361          WRITE ( 14 )  'rad_lw_out_av       ';  WRITE ( 14 )  rad_lw_out_av 
3362       ENDIF
3363       IF ( ALLOCATED( rad_lw_out_change_0 ) )  THEN
3364          WRITE ( 14 )  'rad_lw_out_change_0 '
3365          WRITE ( 14 )  rad_lw_out_change_0
3366       ENDIF
3367       IF ( ALLOCATED( rad_lw_cs_hr ) )  THEN
3368          WRITE ( 14 )  'rad_lw_cs_hr        ';  WRITE ( 14 )  rad_lw_cs_hr
3369       ENDIF
3370       IF ( ALLOCATED( rad_lw_cs_hr_av ) )  THEN
3371          WRITE ( 14 )  'rad_lw_cs_hr_av     ';  WRITE ( 14 )  rad_lw_cs_hr_av
3372       ENDIF
3373       IF ( ALLOCATED( rad_lw_hr ) )  THEN
3374          WRITE ( 14 )  'rad_lw_hr           ';  WRITE ( 14 )  rad_lw_hr
3375       ENDIF
3376       IF ( ALLOCATED( rad_lw_hr_av ) )  THEN
3377          WRITE ( 14 )  'rad_lw_hr_av        ';  WRITE ( 14 )  rad_lw_hr_av
3378       ENDIF
3379       IF ( ALLOCATED( rad_sw_in ) )  THEN
3380          WRITE ( 14 )  'rad_sw_in           ';  WRITE ( 14 )  rad_sw_in 
3381       ENDIF
3382       IF ( ALLOCATED( rad_sw_in_av ) )  THEN
3383          WRITE ( 14 )  'rad_sw_in_av        ';  WRITE ( 14 )  rad_sw_in_av 
3384       ENDIF
3385       IF ( ALLOCATED( rad_sw_out ) )  THEN
3386          WRITE ( 14 )  'rad_sw_out          ';  WRITE ( 14 )  rad_sw_out 
3387       ENDIF
3388       IF ( ALLOCATED( rad_sw_out_av ) )  THEN
3389          WRITE ( 14 )  'rad_sw_out_av       ';  WRITE ( 14 )  rad_sw_out_av 
3390       ENDIF
3391       IF ( ALLOCATED( rad_sw_cs_hr ) )  THEN
3392          WRITE ( 14 )  'rad_sw_cs_hr        ';  WRITE ( 14 )  rad_sw_cs_hr
3393       ENDIF
3394       IF ( ALLOCATED( rad_sw_cs_hr_av ) )  THEN
3395          WRITE ( 14 )  'rad_sw_cs_hr_av     ';  WRITE ( 14 )  rad_sw_cs_hr_av
3396       ENDIF
3397       IF ( ALLOCATED( rad_sw_hr ) )  THEN
3398          WRITE ( 14 )  'rad_sw_hr           ';  WRITE ( 14 )  rad_sw_hr
3399       ENDIF
3400       IF ( ALLOCATED( rad_sw_hr_av ) )  THEN
3401          WRITE ( 14 )  'rad_sw_hr_av        ';  WRITE ( 14 )  rad_sw_hr_av
3402       ENDIF
3403
3404       WRITE ( 14 )  '*** end rad ***     '
3405
3406    ENDIF
3407
3408 END SUBROUTINE radiation_last_actions
3409
3410
3411SUBROUTINE radiation_read_restart_data( i, nxlfa, nxl_on_file, nxrfa, nxr_on_file,   &
3412                                     nynfa, nyn_on_file, nysfa, nys_on_file,   &
3413                                     offset_xa, offset_ya, overlap_count,      &
3414                                     tmp_2d, tmp_3d )
3415 
3416
3417    USE control_parameters
3418       
3419    USE indices
3420   
3421    USE kinds
3422   
3423    USE pegrid
3424
3425    IMPLICIT NONE
3426
3427    CHARACTER (LEN=20) :: field_char   !<
3428
3429    INTEGER(iwp) ::  i               !<
3430    INTEGER(iwp) ::  k               !<
3431    INTEGER(iwp) ::  nxlc            !<
3432    INTEGER(iwp) ::  nxlf            !<
3433    INTEGER(iwp) ::  nxl_on_file     !<
3434    INTEGER(iwp) ::  nxrc            !<
3435    INTEGER(iwp) ::  nxrf            !<
3436    INTEGER(iwp) ::  nxr_on_file     !<
3437    INTEGER(iwp) ::  nync            !<
3438    INTEGER(iwp) ::  nynf            !<
3439    INTEGER(iwp) ::  nyn_on_file     !<
3440    INTEGER(iwp) ::  nysc            !<
3441    INTEGER(iwp) ::  nysf            !<
3442    INTEGER(iwp) ::  nys_on_file     !<
3443    INTEGER(iwp) ::  overlap_count   !<
3444
3445    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  nxlfa       !<
3446    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  nxrfa       !<
3447    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  nynfa       !<
3448    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  nysfa       !<
3449    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  offset_xa   !<
3450    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  offset_ya   !<
3451
3452    REAL(wp),                                                                  &
3453       DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) ::&
3454          tmp_2d   !<
3455
3456    REAL(wp),                                                                  &
3457       DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) ::&
3458          tmp_3d   !<
3459
3460    REAL(wp),                                                                  &
3461       DIMENSION(0:0,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) ::&
3462          tmp_3d2   !<
3463
3464
3465
3466   IF ( initializing_actions == 'read_restart_data' )  THEN
3467      READ ( 13 )  field_char
3468
3469      DO  WHILE ( TRIM( field_char ) /= '*** end rad ***' )
3470
3471         DO  k = 1, overlap_count
3472
3473            nxlf = nxlfa(i,k)
3474            nxlc = nxlfa(i,k) + offset_xa(i,k)
3475            nxrf = nxrfa(i,k)
3476            nxrc = nxrfa(i,k) + offset_xa(i,k)
3477            nysf = nysfa(i,k)
3478            nysc = nysfa(i,k) + offset_ya(i,k)
3479            nynf = nynfa(i,k)
3480            nync = nynfa(i,k) + offset_ya(i,k)
3481
3482
3483            SELECT CASE ( TRIM( field_char ) )
3484
3485                CASE ( 'rad_net' )
3486                   IF ( .NOT. ALLOCATED( rad_net ) )  THEN
3487                      ALLOCATE( rad_net(nysg:nyng,nxlg:nxrg) )
3488                   ENDIF 
3489                   IF ( k == 1 )  READ ( 13 )  tmp_2d
3490                   rad_net(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
3491                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3492
3493                CASE ( 'rad_net_av' )
3494                   IF ( .NOT. ALLOCATED( rad_net_av ) )  THEN
3495                      ALLOCATE( rad_net_av(nysg:nyng,nxlg:nxrg) )
3496                   ENDIF 
3497                   IF ( k == 1 )  READ ( 13 )  tmp_2d
3498                   rad_net_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
3499                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3500                CASE ( 'rad_lw_in' )
3501                   IF ( .NOT. ALLOCATED( rad_lw_in ) )  THEN
3502                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3503                           radiation_scheme == 'constant')  THEN
3504                         ALLOCATE( rad_lw_in(0:0,nysg:nyng,nxlg:nxrg) )
3505                      ELSE
3506                         ALLOCATE( rad_lw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3507                      ENDIF
3508                   ENDIF 
3509                   IF ( k == 1 )  THEN
3510                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3511                           radiation_scheme == 'constant')  THEN
3512                         READ ( 13 )  tmp_3d2
3513                         rad_lw_in(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3514                             tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3515                      ELSE
3516                         READ ( 13 )  tmp_3d
3517                         rad_lw_in(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3518                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3519                      ENDIF
3520                   ENDIF
3521
3522                CASE ( 'rad_lw_in_av' )
3523                   IF ( .NOT. ALLOCATED( rad_lw_in_av ) )  THEN
3524                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3525                           radiation_scheme == 'constant')  THEN
3526                         ALLOCATE( rad_lw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
3527                      ELSE
3528                         ALLOCATE( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3529                      ENDIF
3530                   ENDIF 
3531                   IF ( k == 1 )  THEN
3532                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3533                           radiation_scheme == 'constant')  THEN
3534                         READ ( 13 )  tmp_3d2
3535                         rad_lw_in_av(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3536                             tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3537                      ELSE
3538                         READ ( 13 )  tmp_3d
3539                         rad_lw_in_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3540                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3541                      ENDIF
3542                   ENDIF
3543
3544                CASE ( 'rad_lw_out' )
3545                   IF ( .NOT. ALLOCATED( rad_lw_out ) )  THEN
3546                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3547                           radiation_scheme == 'constant')  THEN
3548                         ALLOCATE( rad_lw_out(0:0,nysg:nyng,nxlg:nxrg) )
3549                      ELSE
3550                         ALLOCATE( rad_lw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3551                      ENDIF
3552                   ENDIF 
3553                   IF ( k == 1 )  THEN
3554                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3555                           radiation_scheme == 'constant')  THEN
3556                         READ ( 13 )  tmp_3d2
3557                         rad_lw_out(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3558                             tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3559                      ELSE
3560                         READ ( 13 )  tmp_3d
3561                         rad_lw_out(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3562                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3563                      ENDIF
3564                   ENDIF
3565
3566                CASE ( 'rad_lw_out_av' )
3567                   IF ( .NOT. ALLOCATED( rad_lw_out_av ) )  THEN
3568                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3569                           radiation_scheme == 'constant')  THEN
3570                         ALLOCATE( rad_lw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
3571                      ELSE
3572                         ALLOCATE( rad_lw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3573                      ENDIF
3574                   ENDIF 
3575                   IF ( k == 1 )  THEN
3576                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3577                           radiation_scheme == 'constant')  THEN
3578                         READ ( 13 )  tmp_3d2
3579                         rad_lw_out_av(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3580                             tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3581                      ELSE
3582                         READ ( 13 )  tmp_3d
3583                         rad_lw_out_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3584                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3585                      ENDIF
3586                   ENDIF
3587
3588                CASE ( 'rad_lw_out_change_0' )
3589                   IF ( .NOT. ALLOCATED( rad_lw_out_change_0 ) )  THEN
3590                      ALLOCATE( rad_lw_out_change_0(nysg:nyng,nxlg:nxrg) )
3591                   ENDIF
3592                   IF ( k == 1 )  READ ( 13 )  tmp_2d
3593                   rad_lw_out_change_0(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)&
3594                              = tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3595
3596                CASE ( 'rad_lw_cs_hr' )
3597                   IF ( .NOT. ALLOCATED( rad_lw_cs_hr ) )  THEN
3598                      ALLOCATE( rad_lw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3599                   ENDIF
3600                   IF ( k == 1 )  READ ( 13 )  tmp_3d
3601                   rad_lw_cs_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3602                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3603
3604                CASE ( 'rad_lw_cs_hr_av' )
3605                   IF ( .NOT. ALLOCATED( rad_lw_cs_hr_av ) )  THEN
3606                      ALLOCATE( rad_lw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3607                   ENDIF
3608                   IF ( k == 1 )  READ ( 13 )  tmp_3d
3609                   rad_lw_cs_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3610                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3611
3612                CASE ( 'rad_lw_hr' )
3613                   IF ( .NOT. ALLOCATED( rad_lw_hr ) )  THEN
3614                      ALLOCATE( rad_lw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3615                   ENDIF
3616                   IF ( k == 1 )  READ ( 13 )  tmp_3d
3617                   rad_lw_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3618                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3619
3620                CASE ( 'rad_lw_hr_av' )
3621                   IF ( .NOT. ALLOCATED( rad_lw_hr_av ) )  THEN
3622                      ALLOCATE( rad_lw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3623                   ENDIF
3624                   IF ( k == 1 )  READ ( 13 )  tmp_3d
3625                   rad_lw_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3626                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3627
3628                CASE ( 'rad_sw_in' )
3629                   IF ( .NOT. ALLOCATED( rad_sw_in ) )  THEN
3630                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3631                           radiation_scheme == 'constant')  THEN
3632                         ALLOCATE( rad_sw_in(0:0,nysg:nyng,nxlg:nxrg) )
3633                      ELSE
3634                         ALLOCATE( rad_sw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3635                      ENDIF
3636                   ENDIF 
3637                   IF ( k == 1 )  THEN
3638                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3639                           radiation_scheme == 'constant')  THEN
3640                         READ ( 13 )  tmp_3d2
3641                         rad_sw_in(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3642                             tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3643                      ELSE
3644                         READ ( 13 )  tmp_3d
3645                         rad_sw_in(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3646                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3647                      ENDIF
3648                   ENDIF
3649
3650                CASE ( 'rad_sw_in_av' )
3651                   IF ( .NOT. ALLOCATED( rad_sw_in_av ) )  THEN
3652                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3653                           radiation_scheme == 'constant')  THEN
3654                         ALLOCATE( rad_sw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
3655                      ELSE
3656                         ALLOCATE( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3657                      ENDIF
3658                   ENDIF 
3659                   IF ( k == 1 )  THEN
3660                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3661                           radiation_scheme == 'constant')  THEN
3662                         READ ( 13 )  tmp_3d2
3663                         rad_sw_in_av(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3664                             tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3665                      ELSE
3666                         READ ( 13 )  tmp_3d
3667                         rad_sw_in_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3668                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3669                      ENDIF
3670                   ENDIF
3671
3672                CASE ( 'rad_sw_out' )
3673                   IF ( .NOT. ALLOCATED( rad_sw_out ) )  THEN
3674                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3675                           radiation_scheme == 'constant')  THEN
3676                         ALLOCATE( rad_sw_out(0:0,nysg:nyng,nxlg:nxrg) )
3677                      ELSE
3678                         ALLOCATE( rad_sw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3679                      ENDIF
3680                   ENDIF 
3681                   IF ( k == 1 )  THEN
3682                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3683                           radiation_scheme == 'constant')  THEN
3684                         READ ( 13 )  tmp_3d2
3685                         rad_sw_out(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3686                             tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3687                      ELSE
3688                         READ ( 13 )  tmp_3d
3689                         rad_sw_out(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3690                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3691                      ENDIF
3692                   ENDIF
3693
3694                CASE ( 'rad_sw_out_av' )
3695                   IF ( .NOT. ALLOCATED( rad_sw_out_av ) )  THEN
3696                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3697                           radiation_scheme == 'constant')  THEN
3698                         ALLOCATE( rad_sw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
3699                      ELSE
3700                         ALLOCATE( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3701                      ENDIF
3702                   ENDIF 
3703                   IF ( k == 1 )  THEN
3704                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3705                           radiation_scheme == 'constant')  THEN
3706                         READ ( 13 )  tmp_3d2
3707                         rad_sw_out_av(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3708                             tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3709                      ELSE
3710                         READ ( 13 )  tmp_3d
3711                         rad_sw_out_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3712                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3713                      ENDIF
3714                   ENDIF
3715
3716                CASE ( 'rad_sw_cs_hr' )
3717                   IF ( .NOT. ALLOCATED( rad_sw_cs_hr ) )  THEN
3718                      ALLOCATE( rad_sw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3719                   ENDIF
3720                   IF ( k == 1 )  READ ( 13 )  tmp_3d
3721                   rad_sw_cs_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3722                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3723
3724                CASE ( 'rad_sw_cs_hr_av' )
3725                   IF ( .NOT. ALLOCATED( rad_sw_cs_hr_av ) )  THEN
3726                      ALLOCATE( rad_sw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3727                   ENDIF
3728                   IF ( k == 1 )  READ ( 13 )  tmp_3d
3729                   rad_sw_cs_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3730                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3731
3732                CASE ( 'rad_sw_hr' )
3733                   IF ( .NOT. ALLOCATED( rad_sw_hr ) )  THEN
3734                      ALLOCATE( rad_sw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3735                   ENDIF
3736                   IF ( k == 1 )  READ ( 13 )  tmp_3d
3737                   rad_sw_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3738                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3739
3740                CASE ( 'rad_sw_hr_av' )
3741                   IF ( .NOT. ALLOCATED( rad_sw_hr_av ) )  THEN
3742                      ALLOCATE( rad_sw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3743                   ENDIF
3744                   IF ( k == 1 )  READ ( 13 )  tmp_3d
3745                   rad_lw_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3746                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3747
3748               CASE DEFAULT
3749                  WRITE( message_string, * ) 'unknown variable named "',       &
3750                                        TRIM( field_char ), '" found in',      &
3751                                        '&data from prior run on PE ', myid
3752                  CALL message( 'radiation_read_restart_data', 'PA0302', 1, 2, 0, 6, &
3753                                 0 )
3754
3755            END SELECT
3756
3757         ENDDO
3758
3759         READ ( 13 )  field_char
3760
3761      ENDDO
3762   ENDIF
3763
3764 END SUBROUTINE radiation_read_restart_data
3765
3766
3767 END MODULE radiation_model_mod
Note: See TracBrowser for help on using the repository browser.