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

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

last commit documented

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