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

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

changes in the course of urban surface model implementation

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