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

Last change on this file since 2200 was 2200, checked in by suehring, 6 years ago

Bugfixes in radiation_model; monotonic limiter in 5th order advection scheme removed

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