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

Last change on this file since 2232 was 2232, checked in by suehring, 4 years ago

Adjustments according new topography and surface-modelling concept implemented

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