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

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

get topograpyh top index via function call

  • Property svn:keywords set to Id
File size: 145.1 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! Get topography top index via Function call
23!
24! Former revisions:
25! -----------------
26! $Id: radiation_model_mod.f90 2317 2017-07-20 17:27:19Z suehring $
27! Improved syntax layout
28!
29! 2298 2017-06-29 09:28:18Z raasch
30! type of write_binary changed from CHARACTER to LOGICAL
31!
32! 2296 2017-06-28 07:53:56Z maronga
33! Added output of rad_sw_out for radiation_scheme = 'constant'
34!
35! 2270 2017-06-09 12:18:47Z maronga
36! Numbering changed (2 timeseries removed)
37!
38! 2249 2017-06-06 13:58:01Z sward
39! Allow for RRTMG runs without humidity/cloud physics
40!
41! 2248 2017-06-06 13:52:54Z sward
42! Error no changed
43!
44! 2233 2017-05-30 18:08:54Z suehring
45!
46! 2232 2017-05-30 17:47:52Z suehring
47! Adjustments to new topography concept
48! Bugfix in read restart
49!
50! 2200 2017-04-11 11:37:51Z suehring
51! Bugfix in call of exchange_horiz_2d and read restart data
52!
53! 2163 2017-03-01 13:23:15Z schwenkel
54! Bugfix in radiation_check_data_output
55!
56! 2157 2017-02-22 15:10:35Z suehring
57! Bugfix in read_restart data
58!
59! 2011 2016-09-19 17:29:57Z kanani
60! Removed CALL of auxiliary SUBROUTINE get_usm_info,
61! flag urban_surface is now defined in module control_parameters.
62!
63! 2007 2016-08-24 15:47:17Z kanani
64! Added calculation of solar directional vector for new urban surface
65! model,
66! accounted for urban_surface model in radiation_check_parameters,
67! correction of comments for zenith angle.
68!
69! 2000 2016-08-20 18:09:15Z knoop
70! Forced header and separation lines into 80 columns
71!
72! 1976 2016-07-27 13:28:04Z maronga
73! Output of 2D/3D/masked data is now directly done within this module. The
74! radiation schemes have been simplified for better usability so that
75! rad_lw_in, rad_lw_out, rad_sw_in, and rad_sw_out are available independent of
76! the radiation code used.
77!
78! 1856 2016-04-13 12:56:17Z maronga
79! Bugfix: allocation of rad_lw_out for radiation_scheme = 'clear-sky'
80!
81! 1853 2016-04-11 09:00:35Z maronga
82! Added routine for radiation_scheme = constant.
83
84! 1849 2016-04-08 11:33:18Z hoffmann
85! Adapted for modularization of microphysics
86!
87! 1826 2016-04-07 12:01:39Z maronga
88! Further modularization.
89!
90! 1788 2016-03-10 11:01:04Z maronga
91! Added new albedo class for pavements / roads.
92!
93! 1783 2016-03-06 18:36:17Z raasch
94! palm-netcdf-module removed in order to avoid a circular module dependency,
95! netcdf-variables moved to netcdf-module, new routine netcdf_handle_error_rad
96! added
97!
98! 1757 2016-02-22 15:49:32Z maronga
99! Added parameter unscheduled_radiation_calls. Bugfix: interpolation of sounding
100! profiles for pressure and temperature above the LES domain.
101!
102! 1709 2015-11-04 14:47:01Z maronga
103! Bugfix: set initial value for rrtm_lwuflx_dt to zero, small formatting
104! corrections
105!
106! 1701 2015-11-02 07:43:04Z maronga
107! Bugfixes: wrong index for output of timeseries, setting of nz_snd_end
108!
109! 1691 2015-10-26 16:17:44Z maronga
110! Added option for spin-up runs without radiation (skip_time_do_radiation). Bugfix
111! in calculation of pressure profiles. Bugfix in calculation of trace gas profiles.
112! Added output of radiative heating rates.
113!
114! 1682 2015-10-07 23:56:08Z knoop
115! Code annotations made doxygen readable
116!
117! 1606 2015-06-29 10:43:37Z maronga
118! Added preprocessor directive __netcdf to allow for compiling without netCDF.
119! Note, however, that RRTMG cannot be used without netCDF.
120!
121! 1590 2015-05-08 13:56:27Z maronga
122! Bugfix: definition of character strings requires same length for all elements
123!
124! 1587 2015-05-04 14:19:01Z maronga
125! Added albedo class for snow
126!
127! 1585 2015-04-30 07:05:52Z maronga
128! Added support for RRTMG
129!
130! 1571 2015-03-12 16:12:49Z maronga
131! Added missing KIND attribute. Removed upper-case variable names
132!
133! 1551 2015-03-03 14:18:16Z maronga
134! Added support for data output. Various variables have been renamed. Added
135! interface for different radiation schemes (currently: clear-sky, constant, and
136! RRTM (not yet implemented).
137!
138! 1496 2014-12-02 17:25:50Z maronga
139! Initial revision
140!
141!
142! Description:
143! ------------
144!> Radiation models and interfaces
145!> @todo move variable definitions used in radiation_init only to the subroutine
146!>       as they are no longer required after initialization.
147!> @todo Output of full column vertical profiles used in RRTMG
148!> @todo Output of other rrtm arrays (such as volume mixing ratios)
149!> @todo Adapt for use with topography
150!>
151!> @note Many variables have a leading dummy dimension (0:0) in order to
152!>       match the assume-size shape expected by the RRTMG model.
153!------------------------------------------------------------------------------!
154 MODULE radiation_model_mod
155 
156    USE arrays_3d,                                                             &
157        ONLY:  dzw, hyp, pt, q, ql, zu, zw
158
159    USE cloud_parameters,                                                      &
160        ONLY:  cp, l_d_cp, rho_l
161
162    USE constants,                                                             &
163        ONLY:  pi
164
165    USE control_parameters,                                                    &
166        ONLY:  cloud_droplets, cloud_physics, g, initializing_actions,         &
167               large_scale_forcing, lsf_surf, phi, pt_surface, rho_surface,    &
168               surface_pressure, time_since_reference_point
169
170    USE indices,                                                               &
171        ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt
172
173    USE kinds
174
175    USE microphysics_mod,                                                      &
176        ONLY:  nc_const, sigma_gc
177
178#if defined ( __netcdf )
179    USE NETCDF
180#endif
181
182#if defined ( __rrtmg )
183    USE parrrsw,                                                               &
184        ONLY:  naerec, nbndsw
185
186    USE parrrtm,                                                               &
187        ONLY:  nbndlw
188
189    USE rrtmg_lw_init,                                                         &
190        ONLY:  rrtmg_lw_ini
191
192    USE rrtmg_sw_init,                                                         &
193        ONLY:  rrtmg_sw_ini
194
195    USE rrtmg_lw_rad,                                                          &
196        ONLY:  rrtmg_lw
197
198    USE rrtmg_sw_rad,                                                          &
199        ONLY:  rrtmg_sw
200#endif
201    USE surface_mod,                                                           &
202        ONLY:  get_topography_top_index
203
204    IMPLICIT NONE
205
206    CHARACTER(10) :: radiation_scheme = 'clear-sky' ! 'constant', 'clear-sky', or 'rrtmg'
207
208!
209!-- Predefined Land surface classes (albedo_type) after Briegleb (1992)
210    CHARACTER(37), DIMENSION(0:17), PARAMETER :: albedo_type_name = (/      &
211                                   'user defined                         ', & !  0
212                                   'ocean                                ', & !  1
213                                   'mixed farming, tall grassland        ', & !  2
214                                   'tall/medium grassland                ', & !  3
215                                   'evergreen shrubland                  ', & !  4
216                                   'short grassland/meadow/shrubland     ', & !  5
217                                   'evergreen needleleaf forest          ', & !  6
218                                   'mixed deciduous evergreen forest     ', & !  7
219                                   'deciduous forest                     ', & !  8
220                                   'tropical evergreen broadleaved forest', & !  9
221                                   'medium/tall grassland/woodland       ', & ! 10
222                                   'desert, sandy                        ', & ! 11
223                                   'desert, rocky                        ', & ! 12
224                                   'tundra                               ', & ! 13
225                                   'land ice                             ', & ! 14
226                                   'sea ice                              ', & ! 15
227                                   'snow                                 ', & ! 16
228                                   'pavement/roads                       '  & ! 17
229                                                         /)
230
231    INTEGER(iwp) :: albedo_type  = 5,    & !< Albedo surface type (default: short grassland)
232                    day,                 & !< current day of the year
233                    day_init     = 172,  & !< day of the year at model start (21/06)
234                    dots_rad     = 0       !< starting index for timeseries output
235
236    LOGICAL ::  unscheduled_radiation_calls = .TRUE., & !< flag parameter indicating whether additional calls of the radiation code are allowed
237                constant_albedo = .FALSE.,            & !< flag parameter indicating whether the albedo may change depending on zenith
238                force_radiation_call = .FALSE.,       & !< flag parameter for unscheduled radiation calls
239                lw_radiation = .TRUE.,                & !< flag parameter indicating whether longwave radiation shall be calculated
240                radiation = .FALSE.,                  & !< flag parameter indicating whether the radiation model is used
241                sun_up    = .TRUE.,                   & !< flag parameter indicating whether the sun is up or down
242                sw_radiation = .TRUE.,                 & !< flag parameter indicing whether shortwave radiation shall be calculated
243                sun_direction = .FALSE.                 !< flag parameter indicing whether solar direction shall be calculated
244
245
246    REAL(wp), PARAMETER :: d_seconds_hour  = 0.000277777777778_wp,  & !< inverse of seconds per hour (1/3600)
247                           d_hours_day    = 0.0416666666667_wp,     & !< inverse of hours per day (1/24)
248                           sigma_sb       = 5.67037321E-8_wp,       & !< Stefan-Boltzmann constant
249                           solar_constant = 1368.0_wp                 !< solar constant at top of atmosphere
250
251    REAL(wp) :: albedo = 9999999.9_wp,           & !< NAMELIST alpha
252                albedo_lw_dif = 9999999.9_wp,    & !< NAMELIST aldif
253                albedo_lw_dir = 9999999.9_wp,    & !< NAMELIST aldir
254                albedo_sw_dif = 9999999.9_wp,    & !< NAMELIST asdif
255                albedo_sw_dir = 9999999.9_wp,    & !< NAMELIST asdir
256                decl_1,                          & !< declination coef. 1
257                decl_2,                          & !< declination coef. 2
258                decl_3,                          & !< declination coef. 3
259                dt_radiation = 0.0_wp,           & !< radiation model timestep
260                emissivity = 0.98_wp,            & !< NAMELIST surface emissivity
261                lambda = 0.0_wp,                 & !< longitude in degrees
262                lon = 0.0_wp,                    & !< longitude in radians
263                lat = 0.0_wp,                    & !< latitude in radians
264                net_radiation = 0.0_wp,          & !< net radiation at surface
265                skip_time_do_radiation = 0.0_wp, & !< Radiation model is not called before this time
266                sky_trans,                       & !< sky transmissivity
267                time_radiation = 0.0_wp,         & !< time since last call of radiation code
268                time_utc,                        & !< current time in UTC
269                time_utc_init = 43200.0_wp         !< UTC time at model start (noon)
270
271    REAL(wp), DIMENSION(0:0) ::  zenith,         & !< cosine of solar zenith angle
272                                 sun_dir_lat,    & !< solar directional vector in latitudes
273                                 sun_dir_lon       !< solar directional vector in longitudes
274
275    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
276                alpha,                       & !< surface broadband albedo (used for clear-sky scheme)
277                rad_lw_out_change_0,         & !< change in LW out due to change in surface temperature
278                rad_net,                     & !< net radiation at the surface
279                rad_net_av                     !< average of rad_net
280
281!
282!-- Land surface albedos for solar zenith angle of 60° after Briegleb (1992)     
283!-- (shortwave, longwave, broadband):   sw,      lw,      bb,
284    REAL(wp), DIMENSION(0:2,1:17), PARAMETER :: albedo_pars = RESHAPE( (/& 
285                                   0.06_wp, 0.06_wp, 0.06_wp,            & !  1
286                                   0.09_wp, 0.28_wp, 0.19_wp,            & !  2
287                                   0.11_wp, 0.33_wp, 0.23_wp,            & !  3
288                                   0.11_wp, 0.33_wp, 0.23_wp,            & !  4
289                                   0.14_wp, 0.34_wp, 0.25_wp,            & !  5
290                                   0.06_wp, 0.22_wp, 0.14_wp,            & !  6
291                                   0.06_wp, 0.27_wp, 0.17_wp,            & !  7
292                                   0.06_wp, 0.31_wp, 0.19_wp,            & !  8
293                                   0.06_wp, 0.22_wp, 0.14_wp,            & !  9
294                                   0.06_wp, 0.28_wp, 0.18_wp,            & ! 10
295                                   0.35_wp, 0.51_wp, 0.43_wp,            & ! 11
296                                   0.24_wp, 0.40_wp, 0.32_wp,            & ! 12
297                                   0.10_wp, 0.27_wp, 0.19_wp,            & ! 13
298                                   0.90_wp, 0.65_wp, 0.77_wp,            & ! 14
299                                   0.90_wp, 0.65_wp, 0.77_wp,            & ! 15
300                                   0.95_wp, 0.70_wp, 0.82_wp,            & ! 16
301                                   0.08_wp, 0.08_wp, 0.08_wp             & ! 17
302                                 /), (/ 3, 17 /) )
303
304    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: &
305                        rad_lw_cs_hr,                  & !< longwave clear sky radiation heating rate (K/s)
306                        rad_lw_cs_hr_av,               & !< average of rad_lw_cs_hr
307                        rad_lw_hr,                     & !< longwave radiation heating rate (K/s)
308                        rad_lw_hr_av,                  & !< average of rad_sw_hr
309                        rad_lw_in,                     & !< incoming longwave radiation (W/m2)
310                        rad_lw_in_av,                  & !< average of rad_lw_in
311                        rad_lw_out,                    & !< outgoing longwave radiation (W/m2)
312                        rad_lw_out_av,                 & !< average of rad_lw_out
313                        rad_sw_cs_hr,                  & !< shortwave clear sky radiation heating rate (K/s)
314                        rad_sw_cs_hr_av,               & !< average of rad_sw_cs_hr
315                        rad_sw_hr,                     & !< shortwave radiation heating rate (K/s)
316                        rad_sw_hr_av,                  & !< average of rad_sw_hr
317                        rad_sw_in,                     & !< incoming shortwave radiation (W/m2)
318                        rad_sw_in_av,                  & !< average of rad_sw_in
319                        rad_sw_out,                    & !< outgoing shortwave radiation (W/m2)
320                        rad_sw_out_av                    !< average of rad_sw_out
321
322
323!
324!-- Variables and parameters used in RRTMG only
325#if defined ( __rrtmg )
326    CHARACTER(LEN=12) :: rrtm_input_file = "RAD_SND_DATA" !< name of the NetCDF input file (sounding data)
327
328
329!
330!-- Flag parameters for RRTMGS (should not be changed)
331    INTEGER(iwp), PARAMETER :: rrtm_idrv     = 1, & !< flag for longwave upward flux calculation option (0,1)
332                               rrtm_inflglw  = 2, & !< flag for lw cloud optical properties (0,1,2)
333                               rrtm_iceflglw = 0, & !< flag for lw ice particle specifications (0,1,2,3)
334                               rrtm_liqflglw = 1, & !< flag for lw liquid droplet specifications
335                               rrtm_inflgsw  = 2, & !< flag for sw cloud optical properties (0,1,2)
336                               rrtm_iceflgsw = 0, & !< flag for sw ice particle specifications (0,1,2,3)
337                               rrtm_liqflgsw = 1    !< flag for sw liquid droplet specifications
338
339!
340!-- The following variables should be only changed with care, as this will
341!-- require further setting of some variables, which is currently not
342!-- implemented (aerosols, ice phase).
343    INTEGER(iwp) :: nzt_rad,           & !< upper vertical limit for radiation calculations
344                    rrtm_icld = 0,     & !< cloud flag (0: clear sky column, 1: cloudy column)
345                    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)
346
347    INTEGER(iwp) :: nc_stat !< local variable for storin the result of netCDF calls for error message handling
348
349    LOGICAL :: snd_exists = .FALSE.      !< flag parameter to check whether a user-defined input files exists
350
351    REAL(wp), PARAMETER :: mol_mass_air_d_wv = 1.607793_wp !< molecular weight dry air / water vapor
352
353    REAL(wp), DIMENSION(:), ALLOCATABLE :: hyp_snd,     & !< hypostatic pressure from sounding data (hPa)
354                                           q_snd,       & !< specific humidity from sounding data (kg/kg) - dummy at the moment
355                                           rrtm_tsfc,   & !< dummy array for storing surface temperature
356                                           t_snd          !< actual temperature from sounding data (hPa)
357
358    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: aldif,          & !< longwave diffuse albedo solar angle of 60°
359                                             aldir,          & !< longwave direct albedo solar angle of 60°
360                                             asdif,          & !< shortwave diffuse albedo solar angle of 60°
361                                             asdir,          & !< shortwave direct albedo solar angle of 60°
362                                             rrtm_ccl4vmr,   & !< CCL4 volume mixing ratio (g/mol)
363                                             rrtm_cfc11vmr,  & !< CFC11 volume mixing ratio (g/mol)
364                                             rrtm_cfc12vmr,  & !< CFC12 volume mixing ratio (g/mol)
365                                             rrtm_cfc22vmr,  & !< CFC22 volume mixing ratio (g/mol)
366                                             rrtm_ch4vmr,    & !< CH4 volume mixing ratio
367                                             rrtm_cicewp,    & !< in-cloud ice water path (g/m²)
368                                             rrtm_cldfr,     & !< cloud fraction (0,1)
369                                             rrtm_cliqwp,    & !< in-cloud liquid water path (g/m²)
370                                             rrtm_co2vmr,    & !< CO2 volume mixing ratio (g/mol)
371                                             rrtm_emis,      & !< surface emissivity (0-1)   
372                                             rrtm_h2ovmr,    & !< H2O volume mixing ratio
373                                             rrtm_n2ovmr,    & !< N2O volume mixing ratio
374                                             rrtm_o2vmr,     & !< O2 volume mixing ratio
375                                             rrtm_o3vmr,     & !< O3 volume mixing ratio
376                                             rrtm_play,      & !< pressure layers (hPa, zu-grid)
377                                             rrtm_plev,      & !< pressure layers (hPa, zw-grid)
378                                             rrtm_reice,     & !< cloud ice effective radius (microns)
379                                             rrtm_reliq,     & !< cloud water drop effective radius (microns)
380                                             rrtm_tlay,      & !< actual temperature (K, zu-grid)
381                                             rrtm_tlev,      & !< actual temperature (K, zw-grid)
382                                             rrtm_lwdflx,    & !< RRTM output of incoming longwave radiation flux (W/m2)
383                                             rrtm_lwdflxc,   & !< RRTM output of outgoing clear sky longwave radiation flux (W/m2)
384                                             rrtm_lwuflx,    & !< RRTM output of outgoing longwave radiation flux (W/m2)
385                                             rrtm_lwuflxc,   & !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
386                                             rrtm_lwuflx_dt, & !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
387                                             rrtm_lwuflxc_dt,& !< RRTM output of outgoing clear sky longwave radiation flux (W/m2)
388                                             rrtm_lwhr,      & !< RRTM output of longwave radiation heating rate (K/d)
389                                             rrtm_lwhrc,     & !< RRTM output of incoming longwave clear sky radiation heating rate (K/d)
390                                             rrtm_swdflx,    & !< RRTM output of incoming shortwave radiation flux (W/m2)
391                                             rrtm_swdflxc,   & !< RRTM output of outgoing clear sky shortwave radiation flux (W/m2)
392                                             rrtm_swuflx,    & !< RRTM output of outgoing shortwave radiation flux (W/m2)
393                                             rrtm_swuflxc,   & !< RRTM output of incoming clear sky shortwave radiation flux (W/m2)
394                                             rrtm_swhr,      & !< RRTM output of shortwave radiation heating rate (K/d)
395                                             rrtm_swhrc        !< RRTM output of incoming shortwave clear sky radiation heating rate (K/d)
396
397!
398!-- Definition of arrays that are currently not used for calling RRTMG (due to setting of flag parameters)
399    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  rad_lw_cs_in,   & !< incoming clear sky longwave radiation (W/m2) (not used)
400                                                rad_lw_cs_out,  & !< outgoing clear sky longwave radiation (W/m2) (not used)
401                                                rad_sw_cs_in,   & !< incoming clear sky shortwave radiation (W/m2) (not used)
402                                                rad_sw_cs_out,  & !< outgoing clear sky shortwave radiation (W/m2) (not used)
403                                                rrtm_aldif,     & !< surface albedo for longwave diffuse radiation
404                                                rrtm_aldir,     & !< surface albedo for longwave direct radiation
405                                                rrtm_asdif,     & !< surface albedo for shortwave diffuse radiation
406                                                rrtm_asdir,     & !< surface albedo for shortwave direct radiation
407                                                rrtm_lw_tauaer, & !< lw aerosol optical depth
408                                                rrtm_lw_taucld, & !< lw in-cloud optical depth
409                                                rrtm_sw_taucld, & !< sw in-cloud optical depth
410                                                rrtm_sw_ssacld, & !< sw in-cloud single scattering albedo
411                                                rrtm_sw_asmcld, & !< sw in-cloud asymmetry parameter
412                                                rrtm_sw_fsfcld, & !< sw in-cloud forward scattering fraction
413                                                rrtm_sw_tauaer, & !< sw aerosol optical depth
414                                                rrtm_sw_ssaaer, & !< sw aerosol single scattering albedo
415                                                rrtm_sw_asmaer, & !< sw aerosol asymmetry parameter
416                                                rrtm_sw_ecaer     !< sw aerosol optical detph at 0.55 microns (rrtm_iaer = 6 only)
417
418#endif
419
420    INTERFACE radiation_check_data_output
421       MODULE PROCEDURE radiation_check_data_output
422    END INTERFACE radiation_check_data_output
423
424    INTERFACE radiation_check_data_output_pr
425       MODULE PROCEDURE radiation_check_data_output_pr
426    END INTERFACE radiation_check_data_output_pr
427 
428    INTERFACE radiation_check_parameters
429       MODULE PROCEDURE radiation_check_parameters
430    END INTERFACE radiation_check_parameters
431 
432    INTERFACE radiation_clearsky
433       MODULE PROCEDURE radiation_clearsky
434    END INTERFACE radiation_clearsky
435 
436    INTERFACE radiation_constant
437       MODULE PROCEDURE radiation_constant
438    END INTERFACE radiation_constant
439 
440    INTERFACE radiation_control
441       MODULE PROCEDURE radiation_control
442    END INTERFACE radiation_control
443
444    INTERFACE radiation_3d_data_averaging
445       MODULE PROCEDURE radiation_3d_data_averaging
446    END INTERFACE radiation_3d_data_averaging
447
448    INTERFACE radiation_data_output_2d
449       MODULE PROCEDURE radiation_data_output_2d
450    END INTERFACE radiation_data_output_2d
451
452    INTERFACE radiation_data_output_3d
453       MODULE PROCEDURE radiation_data_output_3d
454    END INTERFACE radiation_data_output_3d
455
456    INTERFACE radiation_data_output_mask
457       MODULE PROCEDURE radiation_data_output_mask
458    END INTERFACE radiation_data_output_mask
459
460    INTERFACE radiation_define_netcdf_grid
461       MODULE PROCEDURE radiation_define_netcdf_grid
462    END INTERFACE radiation_define_netcdf_grid
463
464    INTERFACE radiation_header
465       MODULE PROCEDURE radiation_header
466    END INTERFACE radiation_header
467 
468    INTERFACE radiation_init
469       MODULE PROCEDURE radiation_init
470    END INTERFACE radiation_init
471
472    INTERFACE radiation_parin
473       MODULE PROCEDURE radiation_parin
474    END INTERFACE radiation_parin
475   
476    INTERFACE radiation_rrtmg
477       MODULE PROCEDURE radiation_rrtmg
478    END INTERFACE radiation_rrtmg
479
480    INTERFACE radiation_tendency
481       MODULE PROCEDURE radiation_tendency
482       MODULE PROCEDURE radiation_tendency_ij
483    END INTERFACE radiation_tendency
484
485    INTERFACE radiation_read_restart_data
486       MODULE PROCEDURE radiation_read_restart_data
487    END INTERFACE radiation_read_restart_data
488
489    INTERFACE radiation_last_actions
490       MODULE PROCEDURE radiation_last_actions
491    END INTERFACE radiation_last_actions
492
493    SAVE
494
495    PRIVATE
496
497!
498!-- Public functions / NEEDS SORTING
499    PUBLIC radiation_check_data_output, radiation_check_data_output_pr,        &
500           radiation_check_parameters, radiation_control,                      &
501           radiation_header, radiation_init, radiation_parin,                  &
502           radiation_3d_data_averaging, radiation_tendency,                    &
503           radiation_data_output_2d, radiation_data_output_3d,                 &
504           radiation_define_netcdf_grid, radiation_last_actions,               &
505           radiation_read_restart_data, radiation_data_output_mask
506   
507!
508!-- Public variables and constants / NEEDS SORTING
509    PUBLIC decl_1, decl_2, decl_3, dots_rad, dt_radiation, force_radiation_call,&
510           lat, lon, rad_net, rad_net_av, radiation, radiation_scheme, rad_lw_in,        &
511           rad_lw_in_av, rad_lw_out, rad_lw_out_av, rad_lw_out_change_0,       &
512           rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av, rad_sw_in,  &
513           rad_sw_in_av, rad_sw_out, rad_sw_out_av, rad_sw_cs_hr,              &
514           rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av, sigma_sb,                 &
515           skip_time_do_radiation, time_radiation, unscheduled_radiation_calls,&
516           zenith, calc_zenith, sun_direction, sun_dir_lat, sun_dir_lon,       &
517           day_init, time_utc_init
518
519
520#if defined ( __rrtmg )
521    PUBLIC rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir
522#endif
523
524 CONTAINS
525
526
527!------------------------------------------------------------------------------!
528! Description:
529! ------------
530!> This subroutine controls the calls of the radiation schemes
531!------------------------------------------------------------------------------!
532    SUBROUTINE radiation_control
533 
534 
535       IMPLICIT NONE
536
537
538       SELECT CASE ( TRIM( radiation_scheme ) )
539
540          CASE ( 'constant' )
541             CALL radiation_constant
542         
543          CASE ( 'clear-sky' )
544             CALL radiation_clearsky
545       
546          CASE ( 'rrtmg' )
547             CALL radiation_rrtmg
548
549          CASE DEFAULT
550
551       END SELECT
552
553
554    END SUBROUTINE radiation_control
555
556!------------------------------------------------------------------------------!
557! Description:
558! ------------
559!> Check data output for radiation model
560!------------------------------------------------------------------------------!
561    SUBROUTINE radiation_check_data_output( var, unit, i, ilen, k )
562 
563 
564       USE control_parameters,                                                 &
565           ONLY: data_output, message_string
566
567       IMPLICIT NONE
568
569       CHARACTER (LEN=*) ::  unit     !<
570       CHARACTER (LEN=*) ::  var      !<
571
572       INTEGER(iwp) :: i
573       INTEGER(iwp) :: ilen
574       INTEGER(iwp) :: k
575
576       SELECT CASE ( TRIM( var ) )
577
578         CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_sw_cs_hr', 'rad_sw_hr' )
579             IF (  .NOT.  radiation  .OR.  radiation_scheme /= 'rrtmg' )  THEN
580                message_string = '"output of "' // TRIM( var ) // '" requi' // &
581                                 'res radiation = .TRUE. and ' //              &
582                                 'radiation_scheme = "rrtmg"'
583                CALL message( 'check_parameters', 'PA0406', 1, 2, 0, 6, 0 )
584             ENDIF
585             unit = 'K/h'     
586
587         CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_sw_in', 'rad_sw_out' )
588             IF (  .NOT.  radiation  .OR.  radiation_scheme /= 'rrtmg' )  THEN
589                message_string = '"output of "' // TRIM( var ) // '" requi' // &
590                                 'res radiation = .TRUE. and ' //              &
591                                 'radiation_scheme = "rrtmg"'
592                CALL message( 'check_parameters', 'PA0406', 1, 2, 0, 6, 0 )
593             ENDIF
594             unit = 'W/m2'   
595
596          CASE ( 'rad_net*', 'rrtm_aldif*', 'rrtm_aldir*', 'rrtm_asdif*',      &
597                 'rrtm_asdir*' )
598             IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
599                message_string = 'illegal value for data_output: "' //         &
600                                 TRIM( var ) // '" & only 2d-horizontal ' //   &
601                                 'cross sections are allowed for this value'
602                CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
603             ENDIF
604             IF (  .NOT.  radiation  .OR.  radiation_scheme /= "rrtmg" )  THEN
605                IF ( TRIM( var ) == 'rrtm_aldif*'  .OR.                        &
606                     TRIM( var ) == 'rrtm_aldir*'  .OR.                        &
607                     TRIM( var ) == 'rrtm_asdif*'  .OR.                        &
608                     TRIM( var ) == 'rrtm_asdir*'      )                       &
609                THEN
610                   message_string = 'output of "' // TRIM( var ) // '" require'&
611                                    // 's radiation = .TRUE. and radiation_sch'&
612                                    // 'eme = "rrtmg"'
613                   CALL message( 'check_parameters', 'PA0409', 1, 2, 0, 6, 0 )
614                ENDIF
615             ENDIF
616
617             IF ( TRIM( var ) == 'rad_net*'      ) unit = 'W/m2'   
618             IF ( TRIM( var ) == 'rrtm_aldif*'   ) unit = ''   
619             IF ( TRIM( var ) == 'rrtm_aldir*'   ) unit = '' 
620             IF ( TRIM( var ) == 'rrtm_asdif*'   ) unit = '' 
621             IF ( TRIM( var ) == 'rrtm_asdir*'   ) unit = '' 
622
623          CASE DEFAULT
624             unit = 'illegal'
625
626       END SELECT
627
628
629    END SUBROUTINE radiation_check_data_output
630
631!------------------------------------------------------------------------------!
632! Description:
633! ------------
634!> Check data output of profiles for radiation model
635!------------------------------------------------------------------------------! 
636    SUBROUTINE radiation_check_data_output_pr( variable, var_count, unit,      &
637               dopr_unit )
638 
639       USE arrays_3d,                                                          &
640           ONLY: zu
641
642       USE control_parameters,                                                 &
643           ONLY: data_output_pr, message_string
644
645       USE indices
646
647       USE profil_parameter
648
649       USE statistics
650
651       IMPLICIT NONE
652   
653       CHARACTER (LEN=*) ::  unit      !<
654       CHARACTER (LEN=*) ::  variable  !<
655       CHARACTER (LEN=*) ::  dopr_unit !< local value of dopr_unit
656 
657       INTEGER(iwp) ::  user_pr_index !<
658       INTEGER(iwp) ::  var_count     !<
659
660       SELECT CASE ( TRIM( variable ) )
661       
662         CASE ( 'rad_net' )
663             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme == 'constant' )&
664             THEN
665                message_string = 'data_output_pr = ' //                        &
666                                 TRIM( data_output_pr(var_count) ) // ' is' // &
667                                 'not available for radiation = .FALSE. or ' //&
668                                 'radiation_scheme = "constant"'
669                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
670             ELSE
671                dopr_index(var_count) = 99
672                dopr_unit  = 'W/m2'
673                hom(:,2,99,:)  = SPREAD( zw, 2, statistic_regions+1 )
674                unit = dopr_unit
675             ENDIF
676
677          CASE ( 'rad_lw_in' )
678             IF ( (  .NOT.  radiation)  .OR.  radiation_scheme == 'constant' ) &
679             THEN
680                message_string = 'data_output_pr = ' //                        &
681                                 TRIM( data_output_pr(var_count) ) // ' is' // &
682                                 'not available for radiation = .FALSE. or ' //&
683                                 'radiation_scheme = "constant"'
684                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
685             ELSE
686                dopr_index(var_count) = 100
687                dopr_unit  = 'W/m2'
688                hom(:,2,100,:)  = SPREAD( zw, 2, statistic_regions+1 )
689                unit = dopr_unit
690             ENDIF
691
692          CASE ( 'rad_lw_out' )
693             IF ( (  .NOT. radiation )  .OR.  radiation_scheme == 'constant' ) &
694             THEN
695                message_string = 'data_output_pr = ' //                        &
696                                 TRIM( data_output_pr(var_count) ) // ' is' // &
697                                 'not available for radiation = .FALSE. or ' //&
698                                 'radiation_scheme = "constant"'
699                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
700             ELSE
701                dopr_index(var_count) = 101
702                dopr_unit  = 'W/m2'
703                hom(:,2,101,:)  = SPREAD( zw, 2, statistic_regions+1 )
704                unit = dopr_unit   
705             ENDIF
706
707          CASE ( 'rad_sw_in' )
708             IF ( (  .NOT. radiation )  .OR.  radiation_scheme == 'constant' ) &
709             THEN
710                message_string = 'data_output_pr = ' //                        &
711                                 TRIM( data_output_pr(var_count) ) // ' is' // &
712                                 'not available for radiation = .FALSE. or ' //&
713                                 'radiation_scheme = "constant"'
714                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
715             ELSE
716                dopr_index(var_count) = 102
717                dopr_unit  = 'W/m2'
718                hom(:,2,102,:)  = SPREAD( zw, 2, statistic_regions+1 )
719                unit = dopr_unit
720             ENDIF
721
722          CASE ( 'rad_sw_out')
723             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme == 'constant' )&
724             THEN
725                message_string = 'data_output_pr = ' //                        &
726                                 TRIM( data_output_pr(var_count) ) // ' is' // &
727                                 'not available for radiation = .FALSE. or ' //&
728                                 'radiation_scheme = "constant"'
729                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
730             ELSE
731                dopr_index(var_count) = 103
732                dopr_unit  = 'W/m2'
733                hom(:,2,103,:)  = SPREAD( zw, 2, statistic_regions+1 )
734                unit = dopr_unit
735             ENDIF
736
737          CASE ( 'rad_lw_cs_hr' )
738             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
739             THEN
740                message_string = 'data_output_pr = ' //                        &
741                                 TRIM( data_output_pr(var_count) ) // ' is' // &
742                                 'not available for radiation = .FALSE. or ' //&
743                                 'radiation_scheme /= "rrtmg"'
744                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
745             ELSE
746                dopr_index(var_count) = 104
747                dopr_unit  = 'K/h'
748                hom(:,2,104,:)  = SPREAD( zu, 2, statistic_regions+1 )
749                unit = dopr_unit
750             ENDIF
751
752          CASE ( 'rad_lw_hr' )
753             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
754             THEN
755                message_string = 'data_output_pr = ' //                        &
756                                 TRIM( data_output_pr(var_count) ) // ' is' // &
757                                 'not available for radiation = .FALSE. or ' //&
758                                 'radiation_scheme /= "rrtmg"'
759                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
760             ELSE
761                dopr_index(var_count) = 105
762                dopr_unit  = 'K/h'
763                hom(:,2,105,:)  = SPREAD( zu, 2, statistic_regions+1 )
764                unit = dopr_unit
765             ENDIF
766
767          CASE ( 'rad_sw_cs_hr' )
768             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
769             THEN
770                message_string = 'data_output_pr = ' //                        &
771                                 TRIM( data_output_pr(var_count) ) // ' is' // &
772                                 'not available for radiation = .FALSE. or ' //&
773                                 'radiation_scheme /= "rrtmg"'
774                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
775             ELSE
776                dopr_index(var_count) = 106
777                dopr_unit  = 'K/h'
778                hom(:,2,106,:)  = SPREAD( zu, 2, statistic_regions+1 )
779                unit = dopr_unit
780             ENDIF
781
782          CASE ( 'rad_sw_hr' )
783             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
784             THEN
785                message_string = 'data_output_pr = ' //                        &
786                                 TRIM( data_output_pr(var_count) ) // ' is' // &
787                                 'not available for radiation = .FALSE. or ' //&
788                                 'radiation_scheme /= "rrtmg"'
789                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
790             ELSE
791                dopr_index(var_count) = 107
792                dopr_unit  = 'K/h'
793                hom(:,2,107,:)  = SPREAD( zu, 2, statistic_regions+1 )
794                unit = dopr_unit
795             ENDIF
796
797
798          CASE DEFAULT
799             unit = 'illegal'
800
801       END SELECT
802
803
804    END SUBROUTINE radiation_check_data_output_pr
805 
806 
807!------------------------------------------------------------------------------!
808! Description:
809! ------------
810!> Check parameters routine for radiation model
811!------------------------------------------------------------------------------!
812    SUBROUTINE radiation_check_parameters
813
814       USE control_parameters,                                                 &
815           ONLY: message_string, topography, urban_surface
816                 
817   
818       IMPLICIT NONE
819       
820
821       IF ( radiation_scheme /= 'constant'   .AND.                             &
822            radiation_scheme /= 'clear-sky'  .AND.                             &
823            radiation_scheme /= 'rrtmg' )  THEN
824          message_string = 'unknown radiation_scheme = '//                     &
825                           TRIM( radiation_scheme )
826          CALL message( 'check_parameters', 'PA0405', 1, 2, 0, 6, 0 )
827       ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
828#if ! defined ( __rrtmg )
829          message_string = 'radiation_scheme = "rrtmg" requires ' //           & 
830                           'compilation of PALM with pre-processor ' //        &
831                           'directive -D__rrtmg'
832          CALL message( 'check_parameters', 'PA0407', 1, 2, 0, 6, 0 )
833#endif
834#if defined ( __rrtmg ) && ! defined( __netcdf )
835          message_string = 'radiation_scheme = "rrtmg" requires ' //           & 
836                           'the use of NetCDF (preprocessor directive ' //     &
837                           '-D__netcdf'
838          CALL message( 'check_parameters', 'PA0412', 1, 2, 0, 6, 0 )
839#endif
840
841       ENDIF
842
843       IF ( albedo_type == 0  .AND.  albedo == 9999999.9_wp  .AND.             &
844            radiation_scheme == 'clear-sky')  THEN
845          message_string = 'radiation_scheme = "clear-sky" in combination' //  & 
846                           'with albedo_type = 0 requires setting of albedo'// &
847                           ' /= 9999999.9'
848          CALL message( 'check_parameters', 'PA0410', 1, 2, 0, 6, 0 )
849       ENDIF
850
851       IF ( albedo_type == 0  .AND.  radiation_scheme == 'rrtmg'  .AND.        &
852          (    albedo_lw_dif == 9999999.9_wp .OR. albedo_lw_dir == 9999999.9_wp&
853          .OR. albedo_sw_dif == 9999999.9_wp .OR. albedo_sw_dir == 9999999.9_wp& 
854          ) ) THEN
855          message_string = 'radiation_scheme = "rrtmg" in combination' //      & 
856                           'with albedo_type = 0 requires setting of ' //      &
857                           'albedo_lw_dif /= 9999999.9' //                     &
858                           'albedo_lw_dir /= 9999999.9' //                     &
859                           'albedo_sw_dif /= 9999999.9 and' //                 &
860                           'albedo_sw_dir /= 9999999.9'
861          CALL message( 'check_parameters', 'PA0411', 1, 2, 0, 6, 0 )
862       ENDIF
863
864!
865!--    The following paramter check is temporarily extended by the urban_surface
866!--    flag, until a better solution comes up to omit this check in case of
867!--    urban surface model is used.
868       IF ( topography /= 'flat'  .AND.  .NOT.  urban_surface )  THEN
869          message_string = 'radiation scheme cannot be used ' //               & 
870                           'in combination with  topography /= "flat"'
871          CALL message( 'check_parameters', 'PA0414', 1, 2, 0, 6, 0 )
872       ENDIF
873 
874    END SUBROUTINE radiation_check_parameters
875 
876 
877!------------------------------------------------------------------------------!
878! Description:
879! ------------
880!> Initialization of the radiation model
881!------------------------------------------------------------------------------!
882    SUBROUTINE radiation_init
883   
884       IMPLICIT NONE
885
886!
887!--    Allocate array for storing the surface net radiation
888       IF ( .NOT. ALLOCATED ( rad_net ) )  THEN
889          ALLOCATE ( rad_net(nysg:nyng,nxlg:nxrg) )
890          rad_net = 0.0_wp
891       ENDIF
892
893!
894!--    Allocate array for storing the surface net radiation
895       IF ( .NOT. ALLOCATED ( rad_lw_out_change_0 ) )  THEN
896          ALLOCATE ( rad_lw_out_change_0(nysg:nyng,nxlg:nxrg) )
897          rad_lw_out_change_0 = 0.0_wp
898       ENDIF
899
900!
901!--    Fix net radiation in case of radiation_scheme = 'constant'
902       IF ( radiation_scheme == 'constant' )  THEN
903          rad_net = net_radiation
904!          radiation = .FALSE.
905!
906!--    Calculate orbital constants
907       ELSE
908          decl_1 = SIN(23.45_wp * pi / 180.0_wp)
909          decl_2 = 2.0_wp * pi / 365.0_wp
910          decl_3 = decl_2 * 81.0_wp
911          lat    = phi * pi / 180.0_wp
912          lon    = lambda * pi / 180.0_wp
913       ENDIF
914
915       IF ( radiation_scheme == 'clear-sky'  .OR.                              &
916            radiation_scheme == 'constant')  THEN
917
918          ALLOCATE ( alpha(nysg:nyng,nxlg:nxrg) )
919
920          IF ( .NOT. ALLOCATED ( rad_sw_in ) )  THEN
921             ALLOCATE ( rad_sw_in(0:0,nysg:nyng,nxlg:nxrg) )
922          ENDIF
923          IF ( .NOT. ALLOCATED ( rad_sw_out ) )  THEN
924             ALLOCATE ( rad_sw_out(0:0,nysg:nyng,nxlg:nxrg) )
925          ENDIF
926
927          IF ( .NOT. ALLOCATED ( rad_sw_in_av ) )  THEN
928             ALLOCATE ( rad_sw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
929          ENDIF
930          IF ( .NOT. ALLOCATED ( rad_sw_out_av ) )  THEN
931             ALLOCATE ( rad_sw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
932          ENDIF
933
934          IF ( .NOT. ALLOCATED ( rad_lw_in ) )  THEN
935             ALLOCATE ( rad_lw_in(0:0,nysg:nyng,nxlg:nxrg) )
936          ENDIF
937          IF ( .NOT. ALLOCATED ( rad_lw_out ) )  THEN
938             ALLOCATE ( rad_lw_out(0:0,nysg:nyng,nxlg:nxrg) )
939          ENDIF
940
941          IF ( .NOT. ALLOCATED ( rad_lw_in_av ) )  THEN
942             ALLOCATE ( rad_lw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
943          ENDIF
944          IF ( .NOT. ALLOCATED ( rad_lw_out_av ) )  THEN
945             ALLOCATE ( rad_lw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
946          ENDIF
947
948          rad_sw_in  = 0.0_wp
949          rad_sw_out = 0.0_wp
950          rad_lw_in  = 0.0_wp
951          rad_lw_out = 0.0_wp
952
953!
954!--       Overwrite albedo if manually set in parameter file
955          IF ( albedo_type /= 0 .AND. albedo == 9999999.9_wp )  THEN
956             albedo = albedo_pars(2,albedo_type)
957          ENDIF
958   
959          alpha = albedo
960 
961!
962!--    Initialization actions for RRTMG
963       ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
964#if defined ( __rrtmg )
965!
966!--       Allocate albedos
967          ALLOCATE ( rrtm_aldif(0:0,nysg:nyng,nxlg:nxrg) )
968          ALLOCATE ( rrtm_aldir(0:0,nysg:nyng,nxlg:nxrg) )
969          ALLOCATE ( rrtm_asdif(0:0,nysg:nyng,nxlg:nxrg) )
970          ALLOCATE ( rrtm_asdir(0:0,nysg:nyng,nxlg:nxrg) )
971          ALLOCATE ( aldif(nysg:nyng,nxlg:nxrg) )
972          ALLOCATE ( aldir(nysg:nyng,nxlg:nxrg) )
973          ALLOCATE ( asdif(nysg:nyng,nxlg:nxrg) )
974          ALLOCATE ( asdir(nysg:nyng,nxlg:nxrg) )
975
976          IF ( albedo_type /= 0 )  THEN
977             IF ( albedo_lw_dif == 9999999.9_wp )  THEN
978                albedo_lw_dif = albedo_pars(0,albedo_type)
979                albedo_lw_dir = albedo_lw_dif
980             ENDIF
981             IF ( albedo_sw_dif == 9999999.9_wp )  THEN
982                albedo_sw_dif = albedo_pars(1,albedo_type)
983                albedo_sw_dir = albedo_sw_dif
984             ENDIF
985          ENDIF
986
987          aldif(:,:) = albedo_lw_dif
988          aldir(:,:) = albedo_lw_dir
989          asdif(:,:) = albedo_sw_dif
990          asdir(:,:) = albedo_sw_dir
991!
992!--       Calculate initial values of current (cosine of) the zenith angle and
993!--       whether the sun is up
994          CALL calc_zenith     
995!
996!--       Calculate initial surface albedo
997          IF ( .NOT. constant_albedo )  THEN
998             CALL calc_albedo
999          ELSE
1000             rrtm_aldif(0,:,:) = aldif(:,:)
1001             rrtm_aldir(0,:,:) = aldir(:,:)
1002             rrtm_asdif(0,:,:) = asdif(:,:) 
1003             rrtm_asdir(0,:,:) = asdir(:,:)   
1004          ENDIF
1005
1006!
1007!--       Allocate surface emissivity
1008          ALLOCATE ( rrtm_emis(0:0,1:nbndlw+1) )
1009          rrtm_emis = emissivity
1010
1011!
1012!--       Allocate 3d arrays of radiative fluxes and heating rates
1013          IF ( .NOT. ALLOCATED ( rad_sw_in ) )  THEN
1014             ALLOCATE ( rad_sw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1015             rad_sw_in = 0.0_wp
1016          ENDIF
1017
1018          IF ( .NOT. ALLOCATED ( rad_sw_in_av ) )  THEN
1019             ALLOCATE ( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1020          ENDIF
1021
1022          IF ( .NOT. ALLOCATED ( rad_sw_out ) )  THEN
1023             ALLOCATE ( rad_sw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1024             rad_sw_out = 0.0_wp
1025          ENDIF
1026
1027          IF ( .NOT. ALLOCATED ( rad_sw_out_av ) )  THEN
1028             ALLOCATE ( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1029          ENDIF
1030
1031          IF ( .NOT. ALLOCATED ( rad_sw_hr ) )  THEN
1032             ALLOCATE ( rad_sw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1033             rad_sw_hr = 0.0_wp
1034          ENDIF
1035
1036          IF ( .NOT. ALLOCATED ( rad_sw_hr_av ) )  THEN
1037             ALLOCATE ( rad_sw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1038             rad_sw_hr_av = 0.0_wp
1039          ENDIF
1040
1041          IF ( .NOT. ALLOCATED ( rad_sw_cs_hr ) )  THEN
1042             ALLOCATE ( rad_sw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1043             rad_sw_cs_hr = 0.0_wp
1044          ENDIF
1045
1046          IF ( .NOT. ALLOCATED ( rad_sw_cs_hr_av ) )  THEN
1047             ALLOCATE ( rad_sw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1048             rad_sw_cs_hr_av = 0.0_wp
1049          ENDIF
1050
1051          IF ( .NOT. ALLOCATED ( rad_lw_in ) )  THEN
1052             ALLOCATE ( rad_lw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1053             rad_lw_in     = 0.0_wp
1054          ENDIF
1055
1056          IF ( .NOT. ALLOCATED ( rad_lw_in_av ) )  THEN
1057             ALLOCATE ( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1058          ENDIF
1059
1060          IF ( .NOT. ALLOCATED ( rad_lw_out ) )  THEN
1061             ALLOCATE ( rad_lw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1062            rad_lw_out    = 0.0_wp
1063          ENDIF
1064
1065          IF ( .NOT. ALLOCATED ( rad_lw_out_av ) )  THEN
1066             ALLOCATE ( rad_lw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1067          ENDIF
1068
1069          IF ( .NOT. ALLOCATED ( rad_lw_hr ) )  THEN
1070             ALLOCATE ( rad_lw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1071             rad_lw_hr = 0.0_wp
1072          ENDIF
1073
1074          IF ( .NOT. ALLOCATED ( rad_lw_hr_av ) )  THEN
1075             ALLOCATE ( rad_lw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1076             rad_lw_hr_av = 0.0_wp
1077          ENDIF
1078
1079          IF ( .NOT. ALLOCATED ( rad_lw_cs_hr ) )  THEN
1080             ALLOCATE ( rad_lw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1081             rad_lw_cs_hr = 0.0_wp
1082          ENDIF
1083
1084          IF ( .NOT. ALLOCATED ( rad_lw_cs_hr_av ) )  THEN
1085             ALLOCATE ( rad_lw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1086             rad_lw_cs_hr_av = 0.0_wp
1087          ENDIF
1088
1089          ALLOCATE ( rad_sw_cs_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1090          ALLOCATE ( rad_sw_cs_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1091          rad_sw_cs_in  = 0.0_wp
1092          rad_sw_cs_out = 0.0_wp
1093
1094          ALLOCATE ( rad_lw_cs_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1095          ALLOCATE ( rad_lw_cs_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1096          rad_lw_cs_in  = 0.0_wp
1097          rad_lw_cs_out = 0.0_wp
1098
1099!
1100!--       Allocate dummy array for storing surface temperature
1101          ALLOCATE ( rrtm_tsfc(1) )
1102
1103!
1104!--       Initialize RRTMG
1105          IF ( lw_radiation )  CALL rrtmg_lw_ini ( cp )
1106          IF ( sw_radiation )  CALL rrtmg_sw_ini ( cp )
1107
1108!
1109!--       Set input files for RRTMG
1110          INQUIRE(FILE="RAD_SND_DATA", EXIST=snd_exists) 
1111          IF ( .NOT. snd_exists )  THEN
1112             rrtm_input_file = "rrtmg_lw.nc"
1113          ENDIF
1114
1115!
1116!--       Read vertical layers for RRTMG from sounding data
1117!--       The routine provides nzt_rad, hyp_snd(1:nzt_rad),
1118!--       t_snd(nzt+2:nzt_rad), rrtm_play(1:nzt_rad), rrtm_plev(1_nzt_rad+1),
1119!--       rrtm_tlay(nzt+2:nzt_rad), rrtm_tlev(nzt+2:nzt_rad+1)
1120          CALL read_sounding_data
1121
1122!
1123!--       Read trace gas profiles from file. This routine provides
1124!--       the rrtm_ arrays (1:nzt_rad+1)
1125          CALL read_trace_gas_data
1126#endif
1127       ENDIF
1128
1129!
1130!--    Perform user actions if required
1131       CALL user_init_radiation
1132
1133!
1134!--    Calculate radiative fluxes at model start
1135       IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
1136
1137          SELECT CASE ( radiation_scheme )
1138             CASE ( 'rrtmg' )
1139                CALL radiation_rrtmg
1140             CASE ( 'clear-sky' )
1141                CALL radiation_clearsky
1142             CASE ( 'constant' )
1143                CALL radiation_constant
1144             CASE DEFAULT
1145          END SELECT
1146
1147       ENDIF
1148
1149       RETURN
1150
1151    END SUBROUTINE radiation_init
1152
1153
1154!------------------------------------------------------------------------------!
1155! Description:
1156! ------------
1157!> A simple clear sky radiation model
1158!------------------------------------------------------------------------------!
1159    SUBROUTINE radiation_clearsky
1160
1161
1162       IMPLICIT NONE
1163
1164       INTEGER(iwp) :: i, j, k   !< loop indices
1165       REAL(wp)     :: exn,   &  !< Exner functions at surface
1166                       exn1,  &  !< Exner functions at first grid level
1167                       pt1       !< potential temperature at first grid level
1168
1169!
1170!--    Calculate current zenith angle
1171       CALL calc_zenith
1172
1173!
1174!--    Calculate sky transmissivity
1175       sky_trans = 0.6_wp + 0.2_wp * zenith(0)
1176
1177!
1178!--    Calculate value of the Exner function
1179       exn = (surface_pressure / 1000.0_wp )**0.286_wp
1180!
1181!--    Calculate radiation fluxes and net radiation (rad_net) for each grid
1182!--    point
1183       DO i = nxlg, nxrg
1184          DO j = nysg, nyng
1185!
1186!--          Obtain vertical index of topography top
1187             k = get_topography_top_index( j, i, 's' )
1188
1189             exn1 = (hyp(k+1) / 100000.0_wp )**0.286_wp
1190
1191             rad_sw_in(0,j,i)  = solar_constant * sky_trans * zenith(0)
1192             rad_sw_out(0,j,i) = alpha(j,i) * rad_sw_in(0,j,i)
1193             rad_lw_out(0,j,i) = emissivity * sigma_sb * (pt(k,j,i) * exn)**4
1194
1195             IF ( cloud_physics )  THEN
1196                pt1 = pt(k+1,j,i) + l_d_cp / exn1 * ql(k+1,j,i)
1197                rad_lw_in(0,j,i)  = 0.8_wp * sigma_sb * (pt1 * exn1)**4
1198             ELSE
1199                rad_lw_in(0,j,i)  = 0.8_wp * sigma_sb * (pt(k+1,j,i) * exn1)**4
1200             ENDIF
1201
1202             rad_net(j,i) = rad_sw_in(0,j,i) - rad_sw_out(0,j,i)               &
1203                            + rad_lw_in(0,j,i) - rad_lw_out(0,j,i)
1204
1205
1206             rad_lw_out_change_0(j,i) = 3.0_wp * sigma_sb * emissivity         &
1207                                        * (pt(k,j,i) * exn) ** 3
1208
1209          ENDDO
1210       ENDDO
1211
1212    END SUBROUTINE radiation_clearsky
1213
1214
1215!------------------------------------------------------------------------------!
1216! Description:
1217! ------------
1218!> This scheme keeps the prescribed net radiation constant during the run
1219!------------------------------------------------------------------------------!
1220    SUBROUTINE radiation_constant
1221
1222
1223       IMPLICIT NONE
1224
1225       INTEGER(iwp) :: i, j, k   !< loop indices
1226       REAL(wp)     :: exn,   &  !< Exner functions at surface
1227                       exn1,  &  !< Exner functions at first grid level
1228                       pt1       !< potential temperature at first grid level
1229
1230!
1231!--    Calculate value of the Exner function
1232       exn = (surface_pressure / 1000.0_wp )**0.286_wp
1233!
1234!--    Prescribe net radiation and estimate the remaining radiative fluxes
1235       DO i = nxlg, nxrg
1236          DO j = nysg, nyng
1237!
1238!--          Obtain vertical index of topography top. So far it is identical to
1239!--          nzb.
1240             k = get_topography_top_index( j, i, 's' )
1241
1242             rad_net(j,i)      = net_radiation
1243
1244             exn1 = (hyp(k+1) / 100000.0_wp )**0.286_wp
1245
1246             IF ( cloud_physics )  THEN
1247                pt1 = pt(k+1,j,i) + l_d_cp / exn1 * ql(k+1,j,i)
1248                rad_lw_in(0,j,i)  = 0.8_wp * sigma_sb * (pt1 * exn1)**4
1249             ELSE
1250                rad_lw_in(0,j,i)  = 0.8_wp * sigma_sb * (pt(k+1,j,i) * exn1)**4
1251             ENDIF
1252
1253             rad_lw_out(0,j,i) = emissivity * sigma_sb * (pt(k,j,i) * exn)**4
1254
1255             rad_sw_in(0,j,i) = ( rad_net(j,i) - rad_lw_in(0,j,i)              &
1256                                  + rad_lw_out(0,j,i) )                        &
1257                                  / ( 1.0_wp - alpha(j,i) )
1258
1259             rad_sw_out(0,j,i) =  alpha(j,i) * rad_sw_in(0,j,i)
1260
1261          ENDDO
1262       ENDDO
1263
1264    END SUBROUTINE radiation_constant
1265
1266!------------------------------------------------------------------------------!
1267! Description:
1268! ------------
1269!> Header output for radiation model
1270!------------------------------------------------------------------------------!
1271    SUBROUTINE radiation_header ( io )
1272
1273
1274       IMPLICIT NONE
1275 
1276       INTEGER(iwp), INTENT(IN) ::  io            !< Unit of the output file
1277   
1278
1279       
1280!
1281!--    Write radiation model header
1282       WRITE( io, 3 )
1283
1284       IF ( radiation_scheme == "constant" )  THEN
1285          WRITE( io, 4 ) net_radiation
1286       ELSEIF ( radiation_scheme == "clear-sky" )  THEN
1287          WRITE( io, 5 )
1288       ELSEIF ( radiation_scheme == "rrtmg" )  THEN
1289          WRITE( io, 6 )
1290          IF ( .NOT. lw_radiation )  WRITE( io, 10 )
1291          IF ( .NOT. sw_radiation )  WRITE( io, 11 )
1292       ENDIF
1293
1294       IF ( albedo_type == 0 )  THEN
1295          WRITE( io, 7 ) albedo
1296       ELSE
1297          WRITE( io, 8 ) TRIM( albedo_type_name(albedo_type) )
1298       ENDIF
1299       IF ( constant_albedo )  THEN
1300          WRITE( io, 9 )
1301       ENDIF
1302       
1303       IF ( radiation .AND. radiation_scheme /= 'constant' )  THEN
1304          WRITE ( io, 1 )  lambda
1305          WRITE ( io, 2 )  day_init, time_utc_init
1306       ENDIF
1307
1308       WRITE( io, 12 ) dt_radiation
1309 
1310
1311 1 FORMAT ('    Geograph. longitude            :   lambda = ',F4.1,' degr')
1312 2 FORMAT ('    Day of the year at model start :   day_init = ',I3             &
1313            /'    UTC time at model start        :   time_utc_init = ',F7.1' s')
1314 3 FORMAT (//' Radiation model information:'/                                  &
1315              ' ----------------------------'/)
1316 4 FORMAT ('    --> Using constant net radiation: net_radiation = ', F6.2,     &
1317           // 'W/m**2')
1318 5 FORMAT ('    --> Simple radiation scheme for clear sky is used (no clouds,',&
1319                   ' default)')
1320 6 FORMAT ('    --> RRTMG scheme is used')
1321 7 FORMAT (/'    User-specific surface albedo: albedo =', F6.3)
1322 8 FORMAT (/'    Albedo is set for land surface type: ', A)
1323 9 FORMAT (/'    --> Albedo is fixed during the run')
132410 FORMAT (/'    --> Longwave radiation is disabled')
132511 FORMAT (/'    --> Shortwave radiation is disabled.')
132612 FORMAT  ('    Timestep: dt_radiation = ', F6.2, '  s')
1327
1328
1329    END SUBROUTINE radiation_header
1330   
1331
1332!------------------------------------------------------------------------------!
1333! Description:
1334! ------------
1335!> Parin for &radiation_par for radiation model
1336!------------------------------------------------------------------------------!
1337    SUBROUTINE radiation_parin
1338
1339
1340       IMPLICIT NONE
1341
1342       CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file
1343       
1344       NAMELIST /radiation_par/   albedo, albedo_type, albedo_lw_dir,          &
1345                                  albedo_lw_dif, albedo_sw_dir, albedo_sw_dif, &
1346                                  constant_albedo, day_init, dt_radiation,     &
1347                                  lambda, lw_radiation, net_radiation,         &
1348                                  radiation_scheme, skip_time_do_radiation,    &
1349                                  sw_radiation, time_utc_init,                 &
1350                                  unscheduled_radiation_calls
1351       
1352       line = ' '
1353       
1354!
1355!--    Try to find radiation model package
1356       REWIND ( 11 )
1357       line = ' '
1358       DO   WHILE ( INDEX( line, '&radiation_par' ) == 0 )
1359          READ ( 11, '(A)', END=10 )  line
1360       ENDDO
1361       BACKSPACE ( 11 )
1362
1363!
1364!--    Read user-defined namelist
1365       READ ( 11, radiation_par )
1366
1367!
1368!--    Set flag that indicates that the radiation model is switched on
1369       radiation = .TRUE.
1370
1371 10    CONTINUE
1372       
1373
1374    END SUBROUTINE radiation_parin
1375
1376
1377!------------------------------------------------------------------------------!
1378! Description:
1379! ------------
1380!> Implementation of the RRTMG radiation_scheme
1381!------------------------------------------------------------------------------!
1382    SUBROUTINE radiation_rrtmg
1383
1384       USE indices,                                                            &
1385           ONLY:  nbgp
1386
1387       USE particle_attributes,                                                &
1388           ONLY:  grid_particles, number_of_particles, particles,              &
1389                  particle_advection_start, prt_count
1390
1391       IMPLICIT NONE
1392
1393#if defined ( __rrtmg )
1394
1395       INTEGER(iwp) :: i, j, k, n !< loop indices
1396
1397       REAL(wp)     ::  s_r2, &   !< weighted sum over all droplets with r^2
1398                        s_r3      !< weighted sum over all droplets with r^3
1399
1400!
1401!--    Calculate current (cosine of) zenith angle and whether the sun is up
1402       CALL calc_zenith     
1403!
1404!--    Calculate surface albedo
1405       IF ( .NOT. constant_albedo )  THEN
1406          CALL calc_albedo
1407       ENDIF
1408
1409!
1410!--    Prepare input data for RRTMG
1411
1412!
1413!--    In case of large scale forcing with surface data, calculate new pressure
1414!--    profile. nzt_rad might be modified by these calls and all required arrays
1415!--    will then be re-allocated
1416       IF ( large_scale_forcing  .AND.  lsf_surf )  THEN
1417          CALL read_sounding_data
1418          CALL read_trace_gas_data
1419       ENDIF
1420!
1421!--    Loop over all grid points
1422       DO i = nxl, nxr
1423          DO j = nys, nyn
1424
1425!
1426!--          Prepare profiles of temperature and H2O volume mixing ratio
1427             rrtm_tlev(0,nzb+1) = pt(nzb,j,i) * ( surface_pressure             &
1428                                                  / 1000.0_wp )**0.286_wp
1429
1430
1431             IF ( cloud_physics )  THEN
1432                DO k = nzb+1, nzt+1
1433                   rrtm_tlay(0,k) = pt(k,j,i) * ( (hyp(k) ) / 100000.0_wp      &
1434                                    )**0.286_wp + l_d_cp * ql(k,j,i)
1435                   rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * (q(k,j,i) - ql(k,j,i))
1436                ENDDO
1437             ELSE
1438                DO k = nzb+1, nzt+1
1439                   rrtm_tlay(0,k) = pt(k,j,i) * ( (hyp(k) ) / 100000.0_wp      &
1440                                    )**0.286_wp
1441                   rrtm_h2ovmr(0,k) = 0.0_wp
1442                ENDDO
1443             ENDIF
1444
1445!
1446!--          Avoid temperature/humidity jumps at the top of the LES domain by
1447!--          linear interpolation from nzt+2 to nzt+7
1448             DO k = nzt+2, nzt+7
1449                rrtm_tlay(0,k) = rrtm_tlay(0,nzt+1)                            &
1450                              + ( rrtm_tlay(0,nzt+8) - rrtm_tlay(0,nzt+1) )    &
1451                              / ( rrtm_play(0,nzt+8) - rrtm_play(0,nzt+1) )    &
1452                              * ( rrtm_play(0,k) - rrtm_play(0,nzt+1) )
1453
1454                rrtm_h2ovmr(0,k) = rrtm_h2ovmr(0,nzt+1)                        &
1455                              + ( rrtm_h2ovmr(0,nzt+8) - rrtm_h2ovmr(0,nzt+1) )&
1456                              / ( rrtm_play(0,nzt+8)   - rrtm_play(0,nzt+1)   )&
1457                              * ( rrtm_play(0,k) - rrtm_play(0,nzt+1) )
1458
1459             ENDDO
1460
1461!--          Linear interpolate to zw grid
1462             DO k = nzb+2, nzt+8
1463                rrtm_tlev(0,k)   = rrtm_tlay(0,k-1) + (rrtm_tlay(0,k) -        &
1464                                   rrtm_tlay(0,k-1))                           &
1465                                   / ( rrtm_play(0,k) - rrtm_play(0,k-1) )     &
1466                                   * ( rrtm_plev(0,k) - rrtm_play(0,k-1) )
1467             ENDDO
1468
1469
1470!
1471!--          Calculate liquid water path and cloud fraction for each column.
1472!--          Note that LWP is required in g/m² instead of kg/kg m.
1473             rrtm_cldfr  = 0.0_wp
1474             rrtm_reliq  = 0.0_wp
1475             rrtm_cliqwp = 0.0_wp
1476             rrtm_icld   = 0
1477
1478             IF ( cloud_physics )  THEN
1479                DO k = nzb+1, nzt+1
1480                   rrtm_cliqwp(0,k) =  ql(k,j,i) * 1000.0_wp *                 &
1481                                       (rrtm_plev(0,k) - rrtm_plev(0,k+1))     &
1482                                       * 100.0_wp / g
1483
1484                   IF ( rrtm_cliqwp(0,k) > 0.0_wp )  THEN
1485                      rrtm_cldfr(0,k) = 1.0_wp
1486                      IF ( rrtm_icld == 0 )  rrtm_icld = 1
1487
1488!
1489!--                   Calculate cloud droplet effective radius
1490                      IF ( cloud_physics )  THEN
1491                         rrtm_reliq(0,k) = 1.0E6_wp * ( 3.0_wp * ql(k,j,i)     &
1492                                           * rho_surface                       &
1493                                           / ( 4.0_wp * pi * nc_const * rho_l )&
1494                                           )**0.33333333333333_wp              &
1495                                           * EXP( LOG( sigma_gc )**2 )
1496
1497                      ELSEIF ( cloud_droplets )  THEN
1498                         number_of_particles = prt_count(k,j,i)
1499
1500                         IF (number_of_particles <= 0)  CYCLE
1501                         particles => grid_particles(k,j,i)%particles(1:number_of_particles)
1502                         s_r2 = 0.0_wp
1503                         s_r3 = 0.0_wp
1504
1505                         DO  n = 1, number_of_particles
1506                            IF ( particles(n)%particle_mask )  THEN
1507                               s_r2 = s_r2 + particles(n)%radius**2 * &
1508                                      particles(n)%weight_factor
1509                               s_r3 = s_r3 + particles(n)%radius**3 * &
1510                                      particles(n)%weight_factor
1511                            ENDIF
1512                         ENDDO
1513
1514                         IF ( s_r2 > 0.0_wp )  rrtm_reliq(0,k) = s_r3 / s_r2
1515
1516                      ENDIF
1517
1518!
1519!--                   Limit effective radius
1520                      IF ( rrtm_reliq(0,k) > 0.0_wp )  THEN
1521                         rrtm_reliq(0,k) = MAX(rrtm_reliq(0,k),2.5_wp)
1522                         rrtm_reliq(0,k) = MIN(rrtm_reliq(0,k),60.0_wp)
1523                     ENDIF
1524                   ENDIF
1525                ENDDO
1526             ENDIF
1527
1528!
1529!--          Set surface temperature
1530             rrtm_tsfc = pt(nzb,j,i) * (surface_pressure / 1000.0_wp )**0.286_wp
1531
1532             IF ( lw_radiation )  THEN
1533               CALL rrtmg_lw( 1, nzt_rad      , rrtm_icld    , rrtm_idrv      ,&
1534               rrtm_play       , rrtm_plev    , rrtm_tlay    , rrtm_tlev      ,&
1535               rrtm_tsfc       , rrtm_h2ovmr  , rrtm_o3vmr   , rrtm_co2vmr    ,&
1536               rrtm_ch4vmr     , rrtm_n2ovmr  , rrtm_o2vmr   , rrtm_cfc11vmr  ,&
1537               rrtm_cfc12vmr   , rrtm_cfc22vmr, rrtm_ccl4vmr , rrtm_emis      ,&
1538               rrtm_inflglw    , rrtm_iceflglw, rrtm_liqflglw, rrtm_cldfr     ,&
1539               rrtm_lw_taucld  , rrtm_cicewp  , rrtm_cliqwp  , rrtm_reice     ,& 
1540               rrtm_reliq      , rrtm_lw_tauaer,                               &
1541               rrtm_lwuflx     , rrtm_lwdflx  , rrtm_lwhr  ,                   &
1542               rrtm_lwuflxc    , rrtm_lwdflxc , rrtm_lwhrc ,                   &
1543               rrtm_lwuflx_dt  ,  rrtm_lwuflxc_dt )
1544
1545!
1546!--             Save fluxes
1547                DO k = nzb, nzt+1
1548                   rad_lw_in(k,j,i)  = rrtm_lwdflx(0,k)
1549                   rad_lw_out(k,j,i) = rrtm_lwuflx(0,k)
1550                ENDDO
1551
1552!
1553!--             Save heating rates (convert from K/d to K/h)
1554                DO k = nzb+1, nzt+1
1555                   rad_lw_hr(k,j,i)     = rrtm_lwhr(0,k)  * d_hours_day
1556                   rad_lw_cs_hr(k,j,i)  = rrtm_lwhrc(0,k) * d_hours_day
1557                ENDDO
1558
1559!
1560!--             Save change in LW heating rate
1561                rad_lw_out_change_0(j,i) = rrtm_lwuflx_dt(0,nzb)
1562
1563             ENDIF
1564
1565             IF ( sw_radiation .AND. sun_up )  THEN
1566                CALL rrtmg_sw( 1, nzt_rad      , rrtm_icld  , rrtm_iaer       ,&
1567               rrtm_play       , rrtm_plev    , rrtm_tlay  , rrtm_tlev        ,&
1568               rrtm_tsfc       , rrtm_h2ovmr  , rrtm_o3vmr , rrtm_co2vmr      ,&
1569               rrtm_ch4vmr     , rrtm_n2ovmr  , rrtm_o2vmr , rrtm_asdir(:,j,i),&
1570               rrtm_asdif(:,j,i), rrtm_aldir(:,j,i), rrtm_aldif(:,j,i), zenith,&
1571               0.0_wp          , day          , solar_constant,   rrtm_inflgsw,&
1572               rrtm_iceflgsw   , rrtm_liqflgsw, rrtm_cldfr , rrtm_sw_taucld   ,&
1573               rrtm_sw_ssacld  , rrtm_sw_asmcld, rrtm_sw_fsfcld, rrtm_cicewp  ,&
1574               rrtm_cliqwp     , rrtm_reice   , rrtm_reliq , rrtm_sw_tauaer   ,&
1575               rrtm_sw_ssaaer     , rrtm_sw_asmaer  , rrtm_sw_ecaer ,          &
1576               rrtm_swuflx     , rrtm_swdflx  , rrtm_swhr  ,                   &
1577               rrtm_swuflxc    , rrtm_swdflxc , rrtm_swhrc )
1578 
1579!
1580!--             Save fluxes
1581                DO k = nzb, nzt+1
1582                   rad_sw_in(k,j,i)  = rrtm_swdflx(0,k)
1583                   rad_sw_out(k,j,i) = rrtm_swuflx(0,k)
1584                ENDDO
1585
1586!
1587!--             Save heating rates (convert from K/d to K/s)
1588                DO k = nzb+1, nzt+1
1589                   rad_sw_hr(k,j,i)     = rrtm_swhr(0,k)  * d_hours_day
1590                   rad_sw_cs_hr(k,j,i)  = rrtm_swhrc(0,k) * d_hours_day
1591                ENDDO
1592
1593             ENDIF
1594
1595!
1596!--          Calculate surface net radiation
1597             rad_net(j,i) = rad_sw_in(nzb,j,i) - rad_sw_out(nzb,j,i)           &
1598                            + rad_lw_in(nzb,j,i) - rad_lw_out(nzb,j,i)
1599
1600          ENDDO
1601       ENDDO
1602
1603       CALL exchange_horiz( rad_lw_in,  nbgp )
1604       CALL exchange_horiz( rad_lw_out, nbgp )
1605       CALL exchange_horiz( rad_lw_hr,    nbgp )
1606       CALL exchange_horiz( rad_lw_cs_hr, nbgp )
1607
1608       CALL exchange_horiz( rad_sw_in,  nbgp )
1609       CALL exchange_horiz( rad_sw_out, nbgp ) 
1610       CALL exchange_horiz( rad_sw_hr,    nbgp )
1611       CALL exchange_horiz( rad_sw_cs_hr, nbgp )
1612
1613       CALL exchange_horiz_2d( rad_net )
1614       CALL exchange_horiz_2d( rad_lw_out_change_0 )
1615#endif
1616
1617    END SUBROUTINE radiation_rrtmg
1618
1619
1620!------------------------------------------------------------------------------!
1621! Description:
1622! ------------
1623!> Calculate the cosine of the zenith angle (variable is called zenith)
1624!------------------------------------------------------------------------------!
1625    SUBROUTINE calc_zenith
1626
1627       IMPLICIT NONE
1628
1629       REAL(wp) ::  declination,  & !< solar declination angle
1630                    hour_angle      !< solar hour angle
1631!
1632!--    Calculate current day and time based on the initial values and simulation
1633!--    time
1634       day = day_init + INT(FLOOR( (time_utc_init + time_since_reference_point)&
1635                               / 86400.0_wp ), KIND=iwp)
1636       time_utc = MOD((time_utc_init + time_since_reference_point), 86400.0_wp)
1637
1638
1639!
1640!--    Calculate solar declination and hour angle   
1641       declination = ASIN( decl_1 * SIN(decl_2 * REAL(day, KIND=wp) - decl_3) )
1642       hour_angle  = 2.0_wp * pi * (time_utc / 86400.0_wp) + lon - pi
1643
1644!
1645!--    Calculate cosine of solar zenith angle
1646       zenith(0) = SIN(lat) * SIN(declination) + COS(lat) * COS(declination)   &
1647                                            * COS(hour_angle)
1648       zenith(0) = MAX(0.0_wp,zenith(0))
1649
1650!
1651!--    Calculate solar directional vector
1652       IF ( sun_direction )  THEN
1653
1654!
1655!--       Direction in longitudes equals to sin(solar_azimuth) * sin(zenith)
1656          sun_dir_lon(0) = -SIN(hour_angle) * COS(declination)
1657
1658!
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)              &
2500                                             + rad_lw_out(k,j,i)
2501                   ENDDO
2502                ENDDO
2503             ENDDO
2504
2505          CASE ( 'rad_lw_cs_hr' )
2506             DO  i = nxlg, nxrg
2507                DO  j = nysg, nyng
2508                   DO  k = nzb, nzt+1
2509                      rad_lw_cs_hr_av(k,j,i) = rad_lw_cs_hr_av(k,j,i)          &
2510                                               + rad_lw_cs_hr(k,j,i)
2511                   ENDDO
2512                ENDDO
2513             ENDDO
2514
2515          CASE ( 'rad_lw_hr' )
2516             DO  i = nxlg, nxrg
2517                DO  j = nysg, nyng
2518                   DO  k = nzb, nzt+1
2519                      rad_lw_hr_av(k,j,i) = rad_lw_hr_av(k,j,i)                &
2520                                            + rad_lw_hr(k,j,i)
2521                   ENDDO
2522                ENDDO
2523             ENDDO
2524
2525          CASE ( 'rad_sw_in' )
2526             DO  i = nxlg, nxrg
2527                DO  j = nysg, nyng
2528                   DO  k = nzb, nzt+1
2529                      rad_sw_in_av(k,j,i) = rad_sw_in_av(k,j,i)                &
2530                                            + rad_sw_in(k,j,i)
2531                   ENDDO
2532                ENDDO
2533             ENDDO
2534
2535          CASE ( 'rad_sw_out' )
2536             DO  i = nxlg, nxrg
2537                DO  j = nysg, nyng
2538                   DO  k = nzb, nzt+1
2539                      rad_sw_out_av(k,j,i) = rad_sw_out_av(k,j,i)              &
2540                                             + rad_sw_out(k,j,i)
2541                   ENDDO
2542                ENDDO
2543             ENDDO
2544
2545          CASE ( 'rad_sw_cs_hr' )
2546             DO  i = nxlg, nxrg
2547                DO  j = nysg, nyng
2548                   DO  k = nzb, nzt+1
2549                      rad_sw_cs_hr_av(k,j,i) = rad_sw_cs_hr_av(k,j,i)          &
2550                                               + rad_sw_cs_hr(k,j,i)
2551                   ENDDO
2552                ENDDO
2553             ENDDO
2554
2555          CASE ( 'rad_sw_hr' )
2556             DO  i = nxlg, nxrg
2557                DO  j = nysg, nyng
2558                   DO  k = nzb, nzt+1
2559                      rad_sw_hr_av(k,j,i) = rad_sw_hr_av(k,j,i)                &
2560                                            + rad_sw_hr(k,j,i)
2561                   ENDDO
2562                ENDDO
2563             ENDDO
2564
2565          CASE DEFAULT
2566             CONTINUE
2567
2568       END SELECT
2569
2570    ELSEIF ( mode == 'average' )  THEN
2571
2572       SELECT CASE ( TRIM( variable ) )
2573
2574         CASE ( 'rad_net*' )
2575             DO  i = nxlg, nxrg
2576                DO  j = nysg, nyng
2577                   rad_net_av(j,i) = rad_net_av(j,i) / REAL( average_count_3d, &
2578                                     KIND=wp )
2579                ENDDO
2580             ENDDO
2581
2582          CASE ( 'rad_lw_in' )
2583             DO  i = nxlg, nxrg
2584                DO  j = nysg, nyng
2585                   DO  k = nzb, nzt+1
2586                      rad_lw_in_av(k,j,i) = rad_lw_in_av(k,j,i)                &
2587                                            / REAL( average_count_3d, KIND=wp )
2588                   ENDDO
2589                ENDDO
2590             ENDDO
2591
2592          CASE ( 'rad_lw_out' )
2593             DO  i = nxlg, nxrg
2594                DO  j = nysg, nyng
2595                   DO  k = nzb, nzt+1
2596                      rad_lw_out_av(k,j,i) = rad_lw_out_av(k,j,i)              &
2597                                             / REAL( average_count_3d, KIND=wp )
2598                   ENDDO
2599                ENDDO
2600             ENDDO
2601
2602          CASE ( 'rad_lw_cs_hr' )
2603             DO  i = nxlg, nxrg
2604                DO  j = nysg, nyng
2605                   DO  k = nzb, nzt+1
2606                      rad_lw_cs_hr_av(k,j,i) = rad_lw_cs_hr_av(k,j,i)          &
2607                                             / REAL( average_count_3d, KIND=wp )
2608                   ENDDO
2609                ENDDO
2610             ENDDO
2611
2612          CASE ( 'rad_lw_hr' )
2613             DO  i = nxlg, nxrg
2614                DO  j = nysg, nyng
2615                   DO  k = nzb, nzt+1
2616                      rad_lw_hr_av(k,j,i) = rad_lw_hr_av(k,j,i)                &
2617                                            / REAL( average_count_3d, KIND=wp )
2618                   ENDDO
2619                ENDDO
2620             ENDDO
2621
2622          CASE ( 'rad_sw_in' )
2623             DO  i = nxlg, nxrg
2624                DO  j = nysg, nyng
2625                   DO  k = nzb, nzt+1
2626                      rad_sw_in_av(k,j,i) = rad_sw_in_av(k,j,i)                &
2627                                            / REAL( average_count_3d, KIND=wp )
2628                   ENDDO
2629                ENDDO
2630             ENDDO
2631
2632          CASE ( 'rad_sw_out' )
2633             DO  i = nxlg, nxrg
2634                DO  j = nysg, nyng
2635                   DO  k = nzb, nzt+1
2636                      rad_sw_out_av(k,j,i) = rad_sw_out_av(k,j,i)              &
2637                                             / REAL( average_count_3d, KIND=wp )
2638                   ENDDO
2639                ENDDO
2640             ENDDO
2641
2642          CASE ( 'rad_sw_cs_hr' )
2643             DO  i = nxlg, nxrg
2644                DO  j = nysg, nyng
2645                   DO  k = nzb, nzt+1
2646                      rad_sw_cs_hr_av(k,j,i) = rad_sw_cs_hr_av(k,j,i)          &
2647                                             / REAL( average_count_3d, KIND=wp )
2648                   ENDDO
2649                ENDDO
2650             ENDDO
2651
2652          CASE ( 'rad_sw_hr' )
2653             DO  i = nxlg, nxrg
2654                DO  j = nysg, nyng
2655                   DO  k = nzb, nzt+1
2656                      rad_sw_hr_av(k,j,i) = rad_sw_hr_av(k,j,i)                &
2657                                            / REAL( average_count_3d, KIND=wp )
2658                   ENDDO
2659                ENDDO
2660             ENDDO
2661
2662       END SELECT
2663
2664    ENDIF
2665
2666END SUBROUTINE radiation_3d_data_averaging
2667
2668
2669!------------------------------------------------------------------------------!
2670!
2671! Description:
2672! ------------
2673!> Subroutine defining appropriate grid for netcdf variables.
2674!> It is called out from subroutine netcdf.
2675!------------------------------------------------------------------------------!
2676SUBROUTINE radiation_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
2677   
2678    IMPLICIT NONE
2679
2680    CHARACTER (LEN=*), INTENT(IN)  ::  var         !<
2681    LOGICAL, INTENT(OUT)           ::  found       !<
2682    CHARACTER (LEN=*), INTENT(OUT) ::  grid_x      !<
2683    CHARACTER (LEN=*), INTENT(OUT) ::  grid_y      !<
2684    CHARACTER (LEN=*), INTENT(OUT) ::  grid_z      !<
2685
2686    found  = .TRUE.
2687
2688
2689!
2690!-- Check for the grid
2691    SELECT CASE ( TRIM( var ) )
2692
2693       CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_sw_cs_hr', 'rad_sw_hr',        &
2694              'rad_lw_cs_hr_xy', 'rad_lw_hr_xy', 'rad_sw_cs_hr_xy',            &
2695              'rad_sw_hr_xy', 'rad_lw_cs_hr_xz', 'rad_lw_hr_xz',               &
2696              'rad_sw_cs_hr_xz', 'rad_sw_hr_xz', 'rad_lw_cs_hr_yz',            &
2697              'rad_lw_hr_yz', 'rad_sw_cs_hr_yz', 'rad_sw_hr_yz' )
2698          grid_x = 'x'
2699          grid_y = 'y'
2700          grid_z = 'zu'
2701
2702       CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_sw_in', 'rad_sw_out',            &
2703              'rad_lw_in_xy', 'rad_lw_out_xy', 'rad_sw_in_xy','rad_sw_out_xy', &
2704              'rad_lw_in_xz', 'rad_lw_out_xz', 'rad_sw_in_xz','rad_sw_out_xz', &
2705              'rad_lw_in_yz', 'rad_lw_out_yz', 'rad_sw_in_yz','rad_sw_out_yz' )
2706          grid_x = 'x'
2707          grid_y = 'y'
2708          grid_z = 'zw'
2709
2710
2711       CASE DEFAULT
2712          found  = .FALSE.
2713          grid_x = 'none'
2714          grid_y = 'none'
2715          grid_z = 'none'
2716
2717        END SELECT
2718
2719    END SUBROUTINE radiation_define_netcdf_grid
2720
2721!------------------------------------------------------------------------------!
2722!
2723! Description:
2724! ------------
2725!> Subroutine defining 3D output variables
2726!------------------------------------------------------------------------------!
2727 SUBROUTINE radiation_data_output_2d( av, variable, found, grid, mode,         &
2728                                      local_pf, two_d )
2729 
2730    USE indices
2731
2732    USE kinds
2733
2734
2735    IMPLICIT NONE
2736
2737    CHARACTER (LEN=*) ::  grid     !<
2738    CHARACTER (LEN=*) ::  mode     !<
2739    CHARACTER (LEN=*) ::  variable !<
2740
2741    INTEGER(iwp) ::  av !<
2742    INTEGER(iwp) ::  i  !<
2743    INTEGER(iwp) ::  j  !<
2744    INTEGER(iwp) ::  k  !<
2745
2746    LOGICAL      ::  found !<
2747    LOGICAL      ::  two_d !< flag parameter that indicates 2D variables (horizontal cross sections)
2748
2749    REAL(wp), DIMENSION(nxlg:nxrg,nysg:nyng,nzb:nzt+1) ::  local_pf !<
2750
2751    found = .TRUE.
2752
2753    SELECT CASE ( TRIM( variable ) )
2754
2755       CASE ( 'rad_net*_xy' )        ! 2d-array
2756          IF ( av == 0 ) THEN
2757             DO  i = nxlg, nxrg
2758                DO  j = nysg, nyng
2759                   local_pf(i,j,nzb+1) = rad_net(j,i)
2760                ENDDO
2761             ENDDO
2762          ELSE
2763             DO  i = nxlg, nxrg
2764                DO  j = nysg, nyng
2765                   local_pf(i,j,nzb+1) = rad_net_av(j,i)
2766                ENDDO
2767             ENDDO
2768          ENDIF
2769          two_d = .TRUE.
2770          grid = 'zu1'
2771
2772 
2773       CASE ( 'rad_lw_in_xy', 'rad_lw_in_xz', 'rad_lw_in_yz' )
2774          IF ( av == 0 ) THEN
2775             DO  i = nxlg, nxrg
2776                DO  j = nysg, nyng
2777                   DO  k = nzb, nzt+1
2778                      local_pf(i,j,k) = rad_lw_in(k,j,i)
2779                   ENDDO
2780                ENDDO
2781             ENDDO
2782          ELSE
2783             DO  i = nxlg, nxrg
2784                DO  j = nysg, nyng
2785                   DO  k = nzb, nzt+1
2786                      local_pf(i,j,k) = rad_lw_in_av(k,j,i)
2787                   ENDDO
2788                ENDDO
2789             ENDDO
2790          ENDIF
2791          IF ( mode == 'xy' )  grid = 'zu'
2792
2793       CASE ( 'rad_lw_out_xy', 'rad_lw_out_xz', 'rad_lw_out_yz' )
2794          IF ( av == 0 ) THEN
2795             DO  i = nxlg, nxrg
2796                DO  j = nysg, nyng
2797                   DO  k = nzb, nzt+1
2798                      local_pf(i,j,k) = rad_lw_out(k,j,i)
2799                   ENDDO
2800                ENDDO
2801             ENDDO
2802          ELSE
2803             DO  i = nxlg, nxrg
2804                DO  j = nysg, nyng
2805                   DO  k = nzb, nzt+1
2806                      local_pf(i,j,k) = rad_lw_out_av(k,j,i)
2807                   ENDDO
2808                ENDDO
2809             ENDDO
2810          ENDIF   
2811          IF ( mode == 'xy' )  grid = 'zu'
2812
2813       CASE ( 'rad_lw_cs_hr_xy', 'rad_lw_cs_hr_xz', 'rad_lw_cs_hr_yz' )
2814          IF ( av == 0 ) THEN
2815             DO  i = nxlg, nxrg
2816                DO  j = nysg, nyng
2817                   DO  k = nzb, nzt+1
2818                      local_pf(i,j,k) = rad_lw_cs_hr(k,j,i)
2819                   ENDDO
2820                ENDDO
2821             ENDDO
2822          ELSE
2823             DO  i = nxlg, nxrg
2824                DO  j = nysg, nyng
2825                   DO  k = nzb, nzt+1
2826                      local_pf(i,j,k) = rad_lw_cs_hr_av(k,j,i)
2827                   ENDDO
2828                ENDDO
2829             ENDDO
2830          ENDIF
2831          IF ( mode == 'xy' )  grid = 'zw'
2832
2833       CASE ( 'rad_lw_hr_xy', 'rad_lw_hr_xz', 'rad_lw_hr_yz' )
2834          IF ( av == 0 ) THEN
2835             DO  i = nxlg, nxrg
2836                DO  j = nysg, nyng
2837                   DO  k = nzb, nzt+1
2838                      local_pf(i,j,k) = rad_lw_hr(k,j,i)
2839                   ENDDO
2840                ENDDO
2841             ENDDO
2842          ELSE
2843             DO  i = nxlg, nxrg
2844                DO  j = nysg, nyng
2845                   DO  k = nzb, nzt+1
2846                      local_pf(i,j,k) = rad_lw_hr_av(k,j,i)
2847                   ENDDO
2848                ENDDO
2849             ENDDO
2850          ENDIF
2851          IF ( mode == 'xy' )  grid = 'zw'
2852
2853       CASE ( 'rad_sw_in_xy', 'rad_sw_in_xz', 'rad_sw_in_yz' )
2854          IF ( av == 0 ) THEN
2855             DO  i = nxlg, nxrg
2856                DO  j = nysg, nyng
2857                   DO  k = nzb, nzt+1
2858                      local_pf(i,j,k) = rad_sw_in(k,j,i)
2859                   ENDDO
2860                ENDDO
2861             ENDDO
2862          ELSE
2863             DO  i = nxlg, nxrg
2864                DO  j = nysg, nyng
2865                   DO  k = nzb, nzt+1
2866                      local_pf(i,j,k) = rad_sw_in_av(k,j,i)
2867                   ENDDO
2868                ENDDO
2869             ENDDO
2870          ENDIF
2871          IF ( mode == 'xy' )  grid = 'zu'
2872
2873       CASE ( 'rad_sw_out_xy', 'rad_sw_out_xz', 'rad_sw_out_yz' )
2874          IF ( av == 0 ) THEN
2875             DO  i = nxlg, nxrg
2876                DO  j = nysg, nyng
2877                   DO  k = nzb, nzt+1
2878                      local_pf(i,j,k) = rad_sw_out(k,j,i)
2879                   ENDDO
2880                ENDDO
2881             ENDDO
2882          ELSE
2883             DO  i = nxlg, nxrg
2884                DO  j = nysg, nyng
2885                   DO  k = nzb, nzt+1
2886                      local_pf(i,j,k) = rad_sw_out_av(k,j,i)
2887                   ENDDO
2888                ENDDO
2889             ENDDO
2890          ENDIF
2891          IF ( mode == 'xy' )  grid = 'zu'
2892
2893       CASE ( 'rad_sw_cs_hr_xy', 'rad_sw_cs_hr_xz', 'rad_sw_cs_hr_yz' )
2894          IF ( av == 0 ) THEN
2895             DO  i = nxlg, nxrg
2896                DO  j = nysg, nyng
2897                   DO  k = nzb, nzt+1
2898                      local_pf(i,j,k) = rad_sw_cs_hr(k,j,i)
2899                   ENDDO
2900                ENDDO
2901             ENDDO
2902          ELSE
2903             DO  i = nxlg, nxrg
2904                DO  j = nysg, nyng
2905                   DO  k = nzb, nzt+1
2906                      local_pf(i,j,k) = rad_sw_cs_hr_av(k,j,i)
2907                   ENDDO
2908                ENDDO
2909             ENDDO
2910          ENDIF
2911          IF ( mode == 'xy' )  grid = 'zw'
2912
2913       CASE ( 'rad_sw_hr_xy', 'rad_sw_hr_xz', 'rad_sw_hr_yz' )
2914          IF ( av == 0 ) THEN
2915             DO  i = nxlg, nxrg
2916                DO  j = nysg, nyng
2917                   DO  k = nzb, nzt+1
2918                      local_pf(i,j,k) = rad_sw_hr(k,j,i)
2919                   ENDDO
2920                ENDDO
2921             ENDDO
2922          ELSE
2923             DO  i = nxlg, nxrg
2924                DO  j = nysg, nyng
2925                   DO  k = nzb, nzt+1
2926                      local_pf(i,j,k) = rad_sw_hr_av(k,j,i)
2927                   ENDDO
2928                ENDDO
2929             ENDDO
2930          ENDIF
2931          IF ( mode == 'xy' )  grid = 'zw'
2932
2933       CASE DEFAULT
2934          found = .FALSE.
2935          grid  = 'none'
2936
2937    END SELECT
2938 
2939 END SUBROUTINE radiation_data_output_2d
2940
2941
2942!------------------------------------------------------------------------------!
2943!
2944! Description:
2945! ------------
2946!> Subroutine defining 3D output variables
2947!------------------------------------------------------------------------------!
2948 SUBROUTINE radiation_data_output_3d( av, variable, found, local_pf )
2949 
2950
2951    USE indices
2952
2953    USE kinds
2954
2955
2956    IMPLICIT NONE
2957
2958    CHARACTER (LEN=*) ::  variable !<
2959
2960    INTEGER(iwp) ::  av    !<
2961    INTEGER(iwp) ::  i     !<
2962    INTEGER(iwp) ::  j     !<
2963    INTEGER(iwp) ::  k     !<
2964
2965    LOGICAL      ::  found !<
2966
2967    REAL(sp), DIMENSION(nxlg:nxrg,nysg:nyng,nzb:nzt+1) ::  local_pf !<
2968
2969
2970    found = .TRUE.
2971
2972
2973    SELECT CASE ( TRIM( variable ) )
2974
2975      CASE ( 'rad_sw_in' )
2976         IF ( av == 0 )  THEN
2977            DO  i = nxlg, nxrg
2978               DO  j = nysg, nyng
2979                  DO  k = nzb, nzt+1
2980                     local_pf(i,j,k) = rad_sw_in(k,j,i)
2981                  ENDDO
2982               ENDDO
2983            ENDDO
2984         ELSE
2985            DO  i = nxlg, nxrg
2986               DO  j = nysg, nyng
2987                  DO  k = nzb, nzt+1
2988                     local_pf(i,j,k) = rad_sw_in_av(k,j,i)
2989                  ENDDO
2990               ENDDO
2991            ENDDO
2992         ENDIF
2993
2994      CASE ( 'rad_sw_out' )
2995         IF ( av == 0 )  THEN
2996            DO  i = nxlg, nxrg
2997               DO  j = nysg, nyng
2998                  DO  k = nzb, nzt+1
2999                     local_pf(i,j,k) = rad_sw_out(k,j,i)
3000                  ENDDO
3001               ENDDO
3002            ENDDO
3003         ELSE
3004            DO  i = nxlg, nxrg
3005               DO  j = nysg, nyng
3006                  DO  k = nzb, nzt+1
3007                     local_pf(i,j,k) = rad_sw_out_av(k,j,i)
3008                  ENDDO
3009               ENDDO
3010            ENDDO
3011         ENDIF
3012
3013      CASE ( 'rad_sw_cs_hr' )
3014         IF ( av == 0 )  THEN
3015            DO  i = nxlg, nxrg
3016               DO  j = nysg, nyng
3017                  DO  k = nzb, nzt+1
3018                     local_pf(i,j,k) = rad_sw_cs_hr(k,j,i)
3019                  ENDDO
3020               ENDDO
3021            ENDDO
3022         ELSE
3023            DO  i = nxlg, nxrg
3024               DO  j = nysg, nyng
3025                  DO  k = nzb, nzt+1
3026                     local_pf(i,j,k) = rad_sw_cs_hr_av(k,j,i)
3027                  ENDDO
3028               ENDDO
3029            ENDDO
3030         ENDIF
3031
3032      CASE ( 'rad_sw_hr' )
3033         IF ( av == 0 )  THEN
3034            DO  i = nxlg, nxrg
3035               DO  j = nysg, nyng
3036                  DO  k = nzb, nzt+1
3037                     local_pf(i,j,k) = rad_sw_hr(k,j,i)
3038                  ENDDO
3039               ENDDO
3040            ENDDO
3041         ELSE
3042            DO  i = nxlg, nxrg
3043               DO  j = nysg, nyng
3044                  DO  k = nzb, nzt+1
3045                     local_pf(i,j,k) = rad_sw_hr_av(k,j,i)
3046                  ENDDO
3047               ENDDO
3048            ENDDO
3049         ENDIF
3050
3051      CASE ( 'rad_lw_in' )
3052         IF ( av == 0 )  THEN
3053            DO  i = nxlg, nxrg
3054               DO  j = nysg, nyng
3055                  DO  k = nzb, nzt+1
3056                     local_pf(i,j,k) = rad_lw_in(k,j,i)
3057                  ENDDO
3058               ENDDO
3059            ENDDO
3060         ELSE
3061            DO  i = nxlg, nxrg
3062               DO  j = nysg, nyng
3063                  DO  k = nzb, nzt+1
3064                     local_pf(i,j,k) = rad_lw_in_av(k,j,i)
3065                  ENDDO
3066               ENDDO
3067            ENDDO
3068         ENDIF
3069
3070      CASE ( 'rad_lw_out' )
3071         IF ( av == 0 )  THEN
3072            DO  i = nxlg, nxrg
3073               DO  j = nysg, nyng
3074                  DO  k = nzb, nzt+1
3075                     local_pf(i,j,k) = rad_lw_out(k,j,i)
3076                  ENDDO
3077               ENDDO
3078            ENDDO
3079         ELSE
3080            DO  i = nxlg, nxrg
3081               DO  j = nysg, nyng
3082                  DO  k = nzb, nzt+1
3083                     local_pf(i,j,k) = rad_lw_out_av(k,j,i)
3084                  ENDDO
3085               ENDDO
3086            ENDDO
3087         ENDIF
3088
3089      CASE ( 'rad_lw_cs_hr' )
3090         IF ( av == 0 )  THEN
3091            DO  i = nxlg, nxrg
3092               DO  j = nysg, nyng
3093                  DO  k = nzb, nzt+1
3094                     local_pf(i,j,k) = rad_lw_cs_hr(k,j,i)
3095                  ENDDO
3096               ENDDO
3097            ENDDO
3098         ELSE
3099            DO  i = nxlg, nxrg
3100               DO  j = nysg, nyng
3101                  DO  k = nzb, nzt+1
3102                     local_pf(i,j,k) = rad_lw_cs_hr_av(k,j,i)
3103                  ENDDO
3104               ENDDO
3105            ENDDO
3106         ENDIF
3107
3108      CASE ( 'rad_lw_hr' )
3109         IF ( av == 0 )  THEN
3110            DO  i = nxlg, nxrg
3111               DO  j = nysg, nyng
3112                  DO  k = nzb, nzt+1
3113                     local_pf(i,j,k) = rad_lw_hr(k,j,i)
3114                  ENDDO
3115               ENDDO
3116            ENDDO
3117         ELSE
3118            DO  i = nxlg, nxrg
3119               DO  j = nysg, nyng
3120                  DO  k = nzb, nzt+1
3121                     local_pf(i,j,k) = rad_lw_hr_av(k,j,i)
3122                  ENDDO
3123               ENDDO
3124            ENDDO
3125         ENDIF
3126
3127       CASE DEFAULT
3128          found = .FALSE.
3129
3130    END SELECT
3131
3132
3133 END SUBROUTINE radiation_data_output_3d
3134
3135!------------------------------------------------------------------------------!
3136!
3137! Description:
3138! ------------
3139!> Subroutine defining masked data output
3140!------------------------------------------------------------------------------!
3141 SUBROUTINE radiation_data_output_mask( av, variable, found, local_pf )
3142 
3143    USE control_parameters
3144       
3145    USE indices
3146   
3147    USE kinds
3148   
3149
3150    IMPLICIT NONE
3151
3152    CHARACTER (LEN=*) ::  variable   !<
3153
3154    INTEGER(iwp) ::  av   !<
3155    INTEGER(iwp) ::  i    !<
3156    INTEGER(iwp) ::  j    !<
3157    INTEGER(iwp) ::  k    !<
3158
3159    LOGICAL ::  found     !<
3160
3161    REAL(wp),                                                                  &
3162       DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) ::  &
3163          local_pf   !<
3164
3165
3166    found = .TRUE.
3167
3168    SELECT CASE ( TRIM( variable ) )
3169
3170
3171       CASE ( 'rad_lw_in' )
3172          IF ( av == 0 )  THEN
3173             DO  i = 1, mask_size_l(mid,1)
3174                DO  j = 1, mask_size_l(mid,2)
3175                   DO  k = 1, mask_size_l(mid,3)
3176                       local_pf(i,j,k) = rad_lw_in(mask_k(mid,k),              &
3177                                            mask_j(mid,j),mask_i(mid,i))
3178                    ENDDO
3179                 ENDDO
3180              ENDDO
3181          ELSE
3182             DO  i = 1, mask_size_l(mid,1)
3183                DO  j = 1, mask_size_l(mid,2)
3184                   DO  k = 1, mask_size_l(mid,3)
3185                       local_pf(i,j,k) = rad_lw_in_av(mask_k(mid,k),           &
3186                                               mask_j(mid,j),mask_i(mid,i))
3187                   ENDDO
3188                ENDDO
3189             ENDDO
3190          ENDIF
3191
3192       CASE ( 'rad_lw_out' )
3193          IF ( av == 0 )  THEN
3194             DO  i = 1, mask_size_l(mid,1)
3195                DO  j = 1, mask_size_l(mid,2)
3196                   DO  k = 1, mask_size_l(mid,3)
3197                       local_pf(i,j,k) = rad_lw_out(mask_k(mid,k),             &
3198                                            mask_j(mid,j),mask_i(mid,i))
3199                    ENDDO
3200                 ENDDO
3201              ENDDO
3202          ELSE
3203             DO  i = 1, mask_size_l(mid,1)
3204                DO  j = 1, mask_size_l(mid,2)
3205                   DO  k = 1, mask_size_l(mid,3)
3206                       local_pf(i,j,k) = rad_lw_out_av(mask_k(mid,k),          &
3207                                               mask_j(mid,j),mask_i(mid,i))
3208                   ENDDO
3209                ENDDO
3210             ENDDO
3211          ENDIF
3212
3213       CASE ( 'rad_lw_cs_hr' )
3214          IF ( av == 0 )  THEN
3215             DO  i = 1, mask_size_l(mid,1)
3216                DO  j = 1, mask_size_l(mid,2)
3217                   DO  k = 1, mask_size_l(mid,3)
3218                       local_pf(i,j,k) = rad_lw_cs_hr(mask_k(mid,k),           &
3219                                            mask_j(mid,j),mask_i(mid,i))
3220                    ENDDO
3221                 ENDDO
3222              ENDDO
3223          ELSE
3224             DO  i = 1, mask_size_l(mid,1)
3225                DO  j = 1, mask_size_l(mid,2)
3226                   DO  k = 1, mask_size_l(mid,3)
3227                       local_pf(i,j,k) = rad_lw_cs_hr_av(mask_k(mid,k),        &
3228                                               mask_j(mid,j),mask_i(mid,i))
3229                   ENDDO
3230                ENDDO
3231             ENDDO
3232          ENDIF
3233
3234       CASE ( 'rad_lw_hr' )
3235          IF ( av == 0 )  THEN
3236             DO  i = 1, mask_size_l(mid,1)
3237                DO  j = 1, mask_size_l(mid,2)
3238                   DO  k = 1, mask_size_l(mid,3)
3239                       local_pf(i,j,k) = rad_lw_hr(mask_k(mid,k),              &
3240                                            mask_j(mid,j),mask_i(mid,i))
3241                    ENDDO
3242                 ENDDO
3243              ENDDO
3244          ELSE
3245             DO  i = 1, mask_size_l(mid,1)
3246                DO  j = 1, mask_size_l(mid,2)
3247                   DO  k = 1, mask_size_l(mid,3)
3248                       local_pf(i,j,k) = rad_lw_hr_av(mask_k(mid,k),           &
3249                                               mask_j(mid,j),mask_i(mid,i))
3250                   ENDDO
3251                ENDDO
3252             ENDDO
3253          ENDIF
3254
3255       CASE ( 'rad_sw_in' )
3256          IF ( av == 0 )  THEN
3257             DO  i = 1, mask_size_l(mid,1)
3258                DO  j = 1, mask_size_l(mid,2)
3259                   DO  k = 1, mask_size_l(mid,3)
3260                       local_pf(i,j,k) = rad_sw_in(mask_k(mid,k),              &
3261                                            mask_j(mid,j),mask_i(mid,i))
3262                    ENDDO
3263                 ENDDO
3264              ENDDO
3265          ELSE
3266             DO  i = 1, mask_size_l(mid,1)
3267                DO  j = 1, mask_size_l(mid,2)
3268                   DO  k = 1, mask_size_l(mid,3)
3269                       local_pf(i,j,k) = rad_sw_in_av(mask_k(mid,k),           &
3270                                               mask_j(mid,j),mask_i(mid,i))
3271                   ENDDO
3272                ENDDO
3273             ENDDO
3274          ENDIF
3275
3276       CASE ( 'rad_sw_out' )
3277          IF ( av == 0 )  THEN
3278             DO  i = 1, mask_size_l(mid,1)
3279                DO  j = 1, mask_size_l(mid,2)
3280                   DO  k = 1, mask_size_l(mid,3)
3281                       local_pf(i,j,k) = rad_sw_out(mask_k(mid,k),             &
3282                                            mask_j(mid,j),mask_i(mid,i))
3283                    ENDDO
3284                 ENDDO
3285              ENDDO
3286          ELSE
3287             DO  i = 1, mask_size_l(mid,1)
3288                DO  j = 1, mask_size_l(mid,2)
3289                   DO  k = 1, mask_size_l(mid,3)
3290                       local_pf(i,j,k) = rad_sw_out_av(mask_k(mid,k),          &
3291                                               mask_j(mid,j),mask_i(mid,i))
3292                   ENDDO
3293                ENDDO
3294             ENDDO
3295          ENDIF
3296
3297       CASE ( 'rad_sw_cs_hr' )
3298          IF ( av == 0 )  THEN
3299             DO  i = 1, mask_size_l(mid,1)
3300                DO  j = 1, mask_size_l(mid,2)
3301                   DO  k = 1, mask_size_l(mid,3)
3302                       local_pf(i,j,k) = rad_sw_cs_hr(mask_k(mid,k),           &
3303                                            mask_j(mid,j),mask_i(mid,i))
3304                    ENDDO
3305                 ENDDO
3306              ENDDO
3307          ELSE
3308             DO  i = 1, mask_size_l(mid,1)
3309                DO  j = 1, mask_size_l(mid,2)
3310                   DO  k = 1, mask_size_l(mid,3)
3311                       local_pf(i,j,k) = rad_sw_cs_hr_av(mask_k(mid,k),        &
3312                                               mask_j(mid,j),mask_i(mid,i))
3313                   ENDDO
3314                ENDDO
3315             ENDDO
3316          ENDIF
3317
3318       CASE ( 'rad_sw_hr' )
3319          IF ( av == 0 )  THEN
3320             DO  i = 1, mask_size_l(mid,1)
3321                DO  j = 1, mask_size_l(mid,2)
3322                   DO  k = 1, mask_size_l(mid,3)
3323                       local_pf(i,j,k) = rad_sw_hr(mask_k(mid,k),              &
3324                                            mask_j(mid,j),mask_i(mid,i))
3325                    ENDDO
3326                 ENDDO
3327              ENDDO
3328          ELSE
3329             DO  i = 1, mask_size_l(mid,1)
3330                DO  j = 1, mask_size_l(mid,2)
3331                   DO  k = 1, mask_size_l(mid,3)
3332                       local_pf(i,j,k) = rad_sw_hr_av(mask_k(mid,k),           &
3333                                               mask_j(mid,j),mask_i(mid,i))
3334                   ENDDO
3335                ENDDO
3336             ENDDO
3337          ENDIF
3338
3339       CASE DEFAULT
3340          found = .FALSE.
3341
3342    END SELECT
3343
3344
3345 END SUBROUTINE radiation_data_output_mask
3346
3347
3348!------------------------------------------------------------------------------!
3349!
3350! Description:
3351! ------------
3352!> Subroutine defines masked output variables
3353!------------------------------------------------------------------------------!
3354 SUBROUTINE radiation_last_actions
3355 
3356
3357    USE control_parameters
3358       
3359    USE kinds
3360
3361    IMPLICIT NONE
3362
3363    IF ( write_binary )  THEN
3364       IF ( ALLOCATED( rad_net ) )  THEN
3365          WRITE ( 14 )  'rad_net             ';  WRITE ( 14 )  rad_net 
3366       ENDIF
3367       IF ( ALLOCATED( rad_net_av ) )  THEN
3368          WRITE ( 14 )  'rad_net_av          ';  WRITE ( 14 )  rad_net_av 
3369       ENDIF 
3370       IF ( ALLOCATED( rad_lw_in ) )  THEN
3371          WRITE ( 14 )  'rad_lw_in           ';  WRITE ( 14 )  rad_lw_in 
3372       ENDIF
3373       IF ( ALLOCATED( rad_lw_in_av ) )  THEN
3374          WRITE ( 14 )  'rad_lw_in_av        ';  WRITE ( 14 )  rad_lw_in_av 
3375       ENDIF
3376       IF ( ALLOCATED( rad_lw_out ) )  THEN
3377          WRITE ( 14 )  'rad_lw_out          ';  WRITE ( 14 )  rad_lw_out
3378       ENDIF
3379       IF ( ALLOCATED( rad_lw_out_av ) )  THEN
3380          WRITE ( 14 )  'rad_lw_out_av       ';  WRITE ( 14 )  rad_lw_out_av 
3381       ENDIF
3382       IF ( ALLOCATED( rad_lw_out_change_0 ) )  THEN
3383          WRITE ( 14 )  'rad_lw_out_change_0 '
3384          WRITE ( 14 )  rad_lw_out_change_0
3385       ENDIF
3386       IF ( ALLOCATED( rad_lw_cs_hr ) )  THEN
3387          WRITE ( 14 )  'rad_lw_cs_hr        ';  WRITE ( 14 )  rad_lw_cs_hr
3388       ENDIF
3389       IF ( ALLOCATED( rad_lw_cs_hr_av ) )  THEN
3390          WRITE ( 14 )  'rad_lw_cs_hr_av     ';  WRITE ( 14 )  rad_lw_cs_hr_av
3391       ENDIF
3392       IF ( ALLOCATED( rad_lw_hr ) )  THEN
3393          WRITE ( 14 )  'rad_lw_hr           ';  WRITE ( 14 )  rad_lw_hr
3394       ENDIF
3395       IF ( ALLOCATED( rad_lw_hr_av ) )  THEN
3396          WRITE ( 14 )  'rad_lw_hr_av        ';  WRITE ( 14 )  rad_lw_hr_av
3397       ENDIF
3398       IF ( ALLOCATED( rad_sw_in ) )  THEN
3399          WRITE ( 14 )  'rad_sw_in           ';  WRITE ( 14 )  rad_sw_in 
3400       ENDIF
3401       IF ( ALLOCATED( rad_sw_in_av ) )  THEN
3402          WRITE ( 14 )  'rad_sw_in_av        ';  WRITE ( 14 )  rad_sw_in_av 
3403       ENDIF
3404       IF ( ALLOCATED( rad_sw_out ) )  THEN
3405          WRITE ( 14 )  'rad_sw_out          ';  WRITE ( 14 )  rad_sw_out 
3406       ENDIF
3407       IF ( ALLOCATED( rad_sw_out_av ) )  THEN
3408          WRITE ( 14 )  'rad_sw_out_av       ';  WRITE ( 14 )  rad_sw_out_av 
3409       ENDIF
3410       IF ( ALLOCATED( rad_sw_cs_hr ) )  THEN
3411          WRITE ( 14 )  'rad_sw_cs_hr        ';  WRITE ( 14 )  rad_sw_cs_hr
3412       ENDIF
3413       IF ( ALLOCATED( rad_sw_cs_hr_av ) )  THEN
3414          WRITE ( 14 )  'rad_sw_cs_hr_av     ';  WRITE ( 14 )  rad_sw_cs_hr_av
3415       ENDIF
3416       IF ( ALLOCATED( rad_sw_hr ) )  THEN
3417          WRITE ( 14 )  'rad_sw_hr           ';  WRITE ( 14 )  rad_sw_hr
3418       ENDIF
3419       IF ( ALLOCATED( rad_sw_hr_av ) )  THEN
3420          WRITE ( 14 )  'rad_sw_hr_av        ';  WRITE ( 14 )  rad_sw_hr_av
3421       ENDIF
3422
3423       WRITE ( 14 )  '*** end rad ***     '
3424
3425    ENDIF
3426
3427 END SUBROUTINE radiation_last_actions
3428
3429
3430SUBROUTINE radiation_read_restart_data( i, nxlfa, nxl_on_file, nxrfa,          &
3431                                        nxr_on_file, nynfa, nyn_on_file, nysfa,&
3432                                        nys_on_file, offset_xa, offset_ya,     &
3433                                        overlap_count, tmp_2d, tmp_3d )
3434 
3435
3436    USE control_parameters
3437       
3438    USE indices
3439   
3440    USE kinds
3441   
3442    USE pegrid
3443
3444    IMPLICIT NONE
3445
3446    CHARACTER (LEN=20) :: field_char   !<
3447
3448    INTEGER(iwp) ::  i               !<
3449    INTEGER(iwp) ::  k               !<
3450    INTEGER(iwp) ::  nxlc            !<
3451    INTEGER(iwp) ::  nxlf            !<
3452    INTEGER(iwp) ::  nxl_on_file     !<
3453    INTEGER(iwp) ::  nxrc            !<
3454    INTEGER(iwp) ::  nxrf            !<
3455    INTEGER(iwp) ::  nxr_on_file     !<
3456    INTEGER(iwp) ::  nync            !<
3457    INTEGER(iwp) ::  nynf            !<
3458    INTEGER(iwp) ::  nyn_on_file     !<
3459    INTEGER(iwp) ::  nysc            !<
3460    INTEGER(iwp) ::  nysf            !<
3461    INTEGER(iwp) ::  nys_on_file     !<
3462    INTEGER(iwp) ::  overlap_count   !<
3463
3464    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  nxlfa       !<
3465    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  nxrfa       !<
3466    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  nynfa       !<
3467    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  nysfa       !<
3468    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  offset_xa   !<
3469    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  offset_ya   !<
3470
3471    REAL(wp),                                                                  &
3472       DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) ::&
3473          tmp_2d   !<
3474
3475    REAL(wp),                                                                  &
3476       DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) ::&
3477          tmp_3d   !<
3478
3479    REAL(wp),                                                                  &
3480       DIMENSION(0:0,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) ::&
3481          tmp_3d2   !<
3482
3483
3484
3485   IF ( initializing_actions == 'read_restart_data' )  THEN
3486      READ ( 13 )  field_char
3487
3488      DO  WHILE ( TRIM( field_char ) /= '*** end rad ***' )
3489
3490         DO  k = 1, overlap_count
3491
3492            nxlf = nxlfa(i,k)
3493            nxlc = nxlfa(i,k) + offset_xa(i,k)
3494            nxrf = nxrfa(i,k)
3495            nxrc = nxrfa(i,k) + offset_xa(i,k)
3496            nysf = nysfa(i,k)
3497            nysc = nysfa(i,k) + offset_ya(i,k)
3498            nynf = nynfa(i,k)
3499            nync = nynfa(i,k) + offset_ya(i,k)
3500
3501
3502            SELECT CASE ( TRIM( field_char ) )
3503
3504                CASE ( 'rad_net' )
3505                   IF ( .NOT. ALLOCATED( rad_net ) )  THEN
3506                      ALLOCATE( rad_net(nysg:nyng,nxlg:nxrg) )
3507                   ENDIF 
3508                   IF ( k == 1 )  READ ( 13 )  tmp_2d
3509                   rad_net(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =         &
3510                                 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3511
3512                CASE ( 'rad_net_av' )
3513                   IF ( .NOT. ALLOCATED( rad_net_av ) )  THEN
3514                      ALLOCATE( rad_net_av(nysg:nyng,nxlg:nxrg) )
3515                   ENDIF 
3516                   IF ( k == 1 )  READ ( 13 )  tmp_2d
3517                   rad_net_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =      &
3518                                 tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3519                CASE ( 'rad_lw_in' )
3520                   IF ( .NOT. ALLOCATED( rad_lw_in ) )  THEN
3521                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3522                           radiation_scheme == 'constant')  THEN
3523                         ALLOCATE( rad_lw_in(0:0,nysg:nyng,nxlg:nxrg) )
3524                      ELSE
3525                         ALLOCATE( rad_lw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3526                      ENDIF
3527                   ENDIF 
3528                   IF ( k == 1 )  THEN
3529                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3530                           radiation_scheme == 'constant')  THEN
3531                         READ ( 13 )  tmp_3d2
3532                         rad_lw_in(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3533                            tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3534                      ELSE
3535                         READ ( 13 )  tmp_3d
3536                         rad_lw_in(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3537                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3538                      ENDIF
3539                   ENDIF
3540
3541                CASE ( 'rad_lw_in_av' )
3542                   IF ( .NOT. ALLOCATED( rad_lw_in_av ) )  THEN
3543                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3544                           radiation_scheme == 'constant')  THEN
3545                         ALLOCATE( rad_lw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
3546                      ELSE
3547                         ALLOCATE( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3548                      ENDIF
3549                   ENDIF 
3550                   IF ( k == 1 )  THEN
3551                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3552                           radiation_scheme == 'constant')  THEN
3553                         READ ( 13 )  tmp_3d2
3554                         rad_lw_in_av(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3555                             tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3556                      ELSE
3557                         READ ( 13 )  tmp_3d
3558                         rad_lw_in_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3559                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3560                      ENDIF
3561                   ENDIF
3562
3563                CASE ( 'rad_lw_out' )
3564                   IF ( .NOT. ALLOCATED( rad_lw_out ) )  THEN
3565                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3566                           radiation_scheme == 'constant')  THEN
3567                         ALLOCATE( rad_lw_out(0:0,nysg:nyng,nxlg:nxrg) )
3568                      ELSE
3569                         ALLOCATE( rad_lw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3570                      ENDIF
3571                   ENDIF 
3572                   IF ( k == 1 )  THEN
3573                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3574                           radiation_scheme == 'constant')  THEN
3575                         READ ( 13 )  tmp_3d2
3576                         rad_lw_out(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3577                             tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3578                      ELSE
3579                         READ ( 13 )  tmp_3d
3580                         rad_lw_out(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3581                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3582                      ENDIF
3583                   ENDIF
3584
3585                CASE ( 'rad_lw_out_av' )
3586                   IF ( .NOT. ALLOCATED( rad_lw_out_av ) )  THEN
3587                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3588                           radiation_scheme == 'constant')  THEN
3589                         ALLOCATE( rad_lw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
3590                      ELSE
3591                         ALLOCATE( rad_lw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3592                      ENDIF
3593                   ENDIF 
3594                   IF ( k == 1 )  THEN
3595                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3596                           radiation_scheme == 'constant')  THEN
3597                         READ ( 13 )  tmp_3d2
3598                         rad_lw_out_av(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3599                             tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3600                      ELSE
3601                         READ ( 13 )  tmp_3d
3602                         rad_lw_out_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3603                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3604                      ENDIF
3605                   ENDIF
3606
3607                CASE ( 'rad_lw_out_change_0' )
3608                   IF ( .NOT. ALLOCATED( rad_lw_out_change_0 ) )  THEN
3609                      ALLOCATE( rad_lw_out_change_0(nysg:nyng,nxlg:nxrg) )
3610                   ENDIF
3611                   IF ( k == 1 )  READ ( 13 )  tmp_2d
3612                   rad_lw_out_change_0(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)&
3613                              = tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3614
3615                CASE ( 'rad_lw_cs_hr' )
3616                   IF ( .NOT. ALLOCATED( rad_lw_cs_hr ) )  THEN
3617                      ALLOCATE( rad_lw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3618                   ENDIF
3619                   IF ( k == 1 )  READ ( 13 )  tmp_3d
3620                   rad_lw_cs_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_cs_hr_av' )
3624                   IF ( .NOT. ALLOCATED( rad_lw_cs_hr_av ) )  THEN
3625                      ALLOCATE( rad_lw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3626                   ENDIF
3627                   IF ( k == 1 )  READ ( 13 )  tmp_3d
3628                   rad_lw_cs_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_lw_hr' )
3632                   IF ( .NOT. ALLOCATED( rad_lw_hr ) )  THEN
3633                      ALLOCATE( rad_lw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3634                   ENDIF
3635                   IF ( k == 1 )  READ ( 13 )  tmp_3d
3636                   rad_lw_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =      &
3637                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3638
3639                CASE ( 'rad_lw_hr_av' )
3640                   IF ( .NOT. ALLOCATED( rad_lw_hr_av ) )  THEN
3641                      ALLOCATE( rad_lw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3642                   ENDIF
3643                   IF ( k == 1 )  READ ( 13 )  tmp_3d
3644                   rad_lw_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =   &
3645                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3646
3647                CASE ( 'rad_sw_in' )
3648                   IF ( .NOT. ALLOCATED( rad_sw_in ) )  THEN
3649                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3650                           radiation_scheme == 'constant')  THEN
3651                         ALLOCATE( rad_sw_in(0:0,nysg:nyng,nxlg:nxrg) )
3652                      ELSE
3653                         ALLOCATE( rad_sw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3654                      ENDIF
3655                   ENDIF 
3656                   IF ( k == 1 )  THEN
3657                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3658                           radiation_scheme == 'constant')  THEN
3659                         READ ( 13 )  tmp_3d2
3660                         rad_sw_in(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3661                             tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3662                      ELSE
3663                         READ ( 13 )  tmp_3d
3664                         rad_sw_in(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3665                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3666                      ENDIF
3667                   ENDIF
3668
3669                CASE ( 'rad_sw_in_av' )
3670                   IF ( .NOT. ALLOCATED( rad_sw_in_av ) )  THEN
3671                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3672                           radiation_scheme == 'constant')  THEN
3673                         ALLOCATE( rad_sw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
3674                      ELSE
3675                         ALLOCATE( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3676                      ENDIF
3677                   ENDIF 
3678                   IF ( k == 1 )  THEN
3679                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3680                           radiation_scheme == 'constant')  THEN
3681                         READ ( 13 )  tmp_3d2
3682                         rad_sw_in_av(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3683                             tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3684                      ELSE
3685                         READ ( 13 )  tmp_3d
3686                         rad_sw_in_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3687                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3688                      ENDIF
3689                   ENDIF
3690
3691                CASE ( 'rad_sw_out' )
3692                   IF ( .NOT. ALLOCATED( rad_sw_out ) )  THEN
3693                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3694                           radiation_scheme == 'constant')  THEN
3695                         ALLOCATE( rad_sw_out(0:0,nysg:nyng,nxlg:nxrg) )
3696                      ELSE
3697                         ALLOCATE( rad_sw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3698                      ENDIF
3699                   ENDIF 
3700                   IF ( k == 1 )  THEN
3701                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3702                           radiation_scheme == 'constant')  THEN
3703                         READ ( 13 )  tmp_3d2
3704                         rad_sw_out(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3705                             tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3706                      ELSE
3707                         READ ( 13 )  tmp_3d
3708                         rad_sw_out(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3709                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3710                      ENDIF
3711                   ENDIF
3712
3713                CASE ( 'rad_sw_out_av' )
3714                   IF ( .NOT. ALLOCATED( rad_sw_out_av ) )  THEN
3715                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3716                           radiation_scheme == 'constant')  THEN
3717                         ALLOCATE( rad_sw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
3718                      ELSE
3719                         ALLOCATE( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3720                      ENDIF
3721                   ENDIF 
3722                   IF ( k == 1 )  THEN
3723                      IF ( radiation_scheme == 'clear-sky'  .OR.               &
3724                           radiation_scheme == 'constant')  THEN
3725                         READ ( 13 )  tmp_3d2
3726                         rad_sw_out_av(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3727                             tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3728                      ELSE
3729                         READ ( 13 )  tmp_3d
3730                         rad_sw_out_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
3731                             tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3732                      ENDIF
3733                   ENDIF
3734
3735                CASE ( 'rad_sw_cs_hr' )
3736                   IF ( .NOT. ALLOCATED( rad_sw_cs_hr ) )  THEN
3737                      ALLOCATE( rad_sw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3738                   ENDIF
3739                   IF ( k == 1 )  READ ( 13 )  tmp_3d
3740                   rad_sw_cs_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_cs_hr_av' )
3744                   IF ( .NOT. ALLOCATED( rad_sw_cs_hr_av ) )  THEN
3745                      ALLOCATE( rad_sw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3746                   ENDIF
3747                   IF ( k == 1 )  READ ( 13 )  tmp_3d
3748                   rad_sw_cs_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
3749                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3750
3751                CASE ( 'rad_sw_hr' )
3752                   IF ( .NOT. ALLOCATED( rad_sw_hr ) )  THEN
3753                      ALLOCATE( rad_sw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3754                   ENDIF
3755                   IF ( k == 1 )  READ ( 13 )  tmp_3d
3756                   rad_sw_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =      &
3757                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3758
3759                CASE ( 'rad_sw_hr_av' )
3760                   IF ( .NOT. ALLOCATED( rad_sw_hr_av ) )  THEN
3761                      ALLOCATE( rad_sw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
3762                   ENDIF
3763                   IF ( k == 1 )  READ ( 13 )  tmp_3d
3764                   rad_lw_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =   &
3765                           tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
3766
3767               CASE DEFAULT
3768                  WRITE( message_string, * ) 'unknown variable named "',       &
3769                                        TRIM( field_char ), '" found in',      &
3770                                        '&data from prior run on PE ', myid
3771                  CALL message( 'radiation_read_restart_data', 'PA0302', 1, 2, &
3772                                0, 6, 0 )
3773
3774            END SELECT
3775
3776         ENDDO
3777
3778         READ ( 13 )  field_char
3779
3780      ENDDO
3781   ENDIF
3782
3783 END SUBROUTINE radiation_read_restart_data
3784
3785
3786 END MODULE radiation_model_mod
Note: See TracBrowser for help on using the repository browser.