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

Last change on this file since 2008 was 2008, checked in by kanani, 5 years ago

last commit documented

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