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

Last change on this file since 1976 was 1976, checked in by maronga, 5 years ago

further modularization of land surface model (2D/3D output and restart data). Bugfix for restart runs without land surface model

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