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

Last change on this file since 2000 was 2000, checked in by knoop, 5 years ago

Forced header and separation lines into 80 columns

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