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

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

major revisions in land surface model

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