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

Last change on this file since 2298 was 2298, checked in by raasch, 5 years ago

write_binary is of type LOGICAL now, MPI2-related code removed, obsolete variables removed, sendrecv_in_background related parts removed, missing variable descriptions added

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