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

Last change on this file since 1827 was 1827, checked in by maronga, 5 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 87.5 KB
Line 
1!> @file radiation_model_mod.f90
2!--------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms
6! of the GNU General Public License as published by the Free Software Foundation,
7! either version 3 of the License, or (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
10! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
11! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with
14! PALM. If not, see <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2016 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------!
18!
19! Current revisions:
20! -----------------
21!
22!
23! Former revisions:
24! -----------------
25! $Id: radiation_model_mod.f90 1827 2016-04-07 12:12:23Z maronga $
26!
27! 1826 2016-04-07 12:01:39Z maronga
28! Further modularization.
29!
30! 1788 2016-03-10 11:01:04Z maronga
31! Added new albedo class for pavements / roads.
32!
33! 1783 2016-03-06 18:36:17Z raasch
34! palm-netcdf-module removed in order to avoid a circular module dependency,
35! netcdf-variables moved to netcdf-module, new routine netcdf_handle_error_rad
36! added
37!
38! 1757 2016-02-22 15:49:32Z maronga
39! Added parameter unscheduled_radiation_calls. Bugfix: interpolation of sounding
40! profiles for pressure and temperature above the LES domain.
41!
42! 1709 2015-11-04 14:47:01Z maronga
43! Bugfix: set initial value for rrtm_lwuflx_dt to zero, small formatting
44! corrections
45!
46! 1701 2015-11-02 07:43:04Z maronga
47! Bugfixes: wrong index for output of timeseries, setting of nz_snd_end
48!
49! 1691 2015-10-26 16:17:44Z maronga
50! Added option for spin-up runs without radiation (skip_time_do_radiation). Bugfix
51! in calculation of pressure profiles. Bugfix in calculation of trace gas profiles.
52! Added output of radiative heating rates.
53!
54! 1682 2015-10-07 23:56:08Z knoop
55! Code annotations made doxygen readable
56!
57! 1606 2015-06-29 10:43:37Z maronga
58! Added preprocessor directive __netcdf to allow for compiling without netCDF.
59! Note, however, that RRTMG cannot be used without netCDF.
60!
61! 1590 2015-05-08 13:56:27Z maronga
62! Bugfix: definition of character strings requires same length for all elements
63!
64! 1587 2015-05-04 14:19:01Z maronga
65! Added albedo class for snow
66!
67! 1585 2015-04-30 07:05:52Z maronga
68! Added support for RRTMG
69!
70! 1571 2015-03-12 16:12:49Z maronga
71! Added missing KIND attribute. Removed upper-case variable names
72!
73! 1551 2015-03-03 14:18:16Z maronga
74! Added support for data output. Various variables have been renamed. Added
75! interface for different radiation schemes (currently: clear-sky, constant, and
76! RRTM (not yet implemented).
77!
78! 1496 2014-12-02 17:25:50Z maronga
79! Initial revision
80!
81!
82! Description:
83! ------------
84!> Radiation models and interfaces
85!> @todo move variable definitions used in radiation_init only to the subroutine
86!>       as they are no longer required after initialization.
87!> @todo Output of full column vertical profiles used in RRTMG
88!> @todo Output of other rrtm arrays (such as volume mixing ratios)
89!> @todo Adapt for use with topography
90!>
91!> @note Many variables have a leading dummy dimension (0:0) in order to
92!>       match the assume-size shape expected by the RRTMG model.
93!------------------------------------------------------------------------------!
94 MODULE radiation_model_mod
95 
96    USE arrays_3d,                                                             &
97        ONLY:  dzw, hyp, pt, q, ql, zw
98
99    USE cloud_parameters,                                                      &
100        ONLY:  cp, l_d_cp, nc_const, rho_l, sigma_gc 
101
102    USE constants,                                                             &
103        ONLY:  pi
104
105    USE control_parameters,                                                    &
106        ONLY:  cloud_droplets, cloud_physics, g, initializing_actions,         &
107               large_scale_forcing, lsf_surf, phi, pt_surface, rho_surface,    &
108               surface_pressure, time_since_reference_point
109
110    USE indices,                                                               &
111        ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb_s_inner, nzb, nzt
112
113    USE kinds
114
115#if defined ( __netcdf )
116    USE NETCDF
117#endif
118
119#if defined ( __rrtmg )
120
121    USE parrrsw,                                                               &
122        ONLY:  naerec, nbndsw
123
124    USE parrrtm,                                                               &
125        ONLY:  nbndlw
126
127    USE rrtmg_lw_init,                                                         &
128        ONLY:  rrtmg_lw_ini
129
130    USE rrtmg_sw_init,                                                         &
131        ONLY:  rrtmg_sw_ini
132
133    USE rrtmg_lw_rad,                                                          &
134        ONLY:  rrtmg_lw
135
136    USE rrtmg_sw_rad,                                                          &
137        ONLY:  rrtmg_sw
138#endif
139
140
141
142    IMPLICIT NONE
143
144    CHARACTER(10) :: radiation_scheme = 'clear-sky' ! 'constant', 'clear-sky', or 'rrtmg'
145
146!
147!-- Predefined Land surface classes (albedo_type) after Briegleb (1992)
148    CHARACTER(37), DIMENSION(0:17), PARAMETER :: albedo_type_name = (/      &
149                                   'user defined                         ', & !  0
150                                   'ocean                                ', & !  1
151                                   'mixed farming, tall grassland        ', & !  2
152                                   'tall/medium grassland                ', & !  3
153                                   'evergreen shrubland                  ', & !  4
154                                   'short grassland/meadow/shrubland     ', & !  5
155                                   'evergreen needleleaf forest          ', & !  6
156                                   'mixed deciduous evergreen forest     ', & !  7
157                                   'deciduous forest                     ', & !  8
158                                   'tropical evergreen broadleaved forest', & !  9
159                                   'medium/tall grassland/woodland       ', & ! 10
160                                   'desert, sandy                        ', & ! 11
161                                   'desert, rocky                        ', & ! 12
162                                   'tundra                               ', & ! 13
163                                   'land ice                             ', & ! 14
164                                   'sea ice                              ', & ! 15
165                                   'snow                                 ', & ! 16
166                                   'pavement/roads                       '  & ! 17
167                                                         /)
168
169    INTEGER(iwp) :: albedo_type  = 5,    & !< Albedo surface type (default: short grassland)
170                    day,                 & !< current day of the year
171                    day_init     = 172,  & !< day of the year at model start (21/06)
172                    dots_rad     = 0       !< starting index for timeseries output
173
174    LOGICAL ::  unscheduled_radiation_calls = .TRUE., & !< flag parameter indicating whether additional calls of the radiation code are allowed
175                constant_albedo = .FALSE.,            & !< flag parameter indicating whether the albedo may change depending on zenith
176                force_radiation_call = .FALSE.,       & !< flag parameter for unscheduled radiation calls
177                lw_radiation = .TRUE.,                & !< flag parameter indicating whether longwave radiation shall be calculated
178                radiation = .FALSE.,                  & !< flag parameter indicating whether the radiation model is used
179                sun_up    = .TRUE.,                   & !< flag parameter indicating whether the sun is up or down
180                sw_radiation = .TRUE.                   !< flag parameter indicing whether shortwave radiation shall be calculated
181
182
183    REAL(wp), PARAMETER :: d_seconds_hour  = 0.000277777777778_wp,  & !< inverse of seconds per hour (1/3600)
184                           d_hours_day    = 0.0416666666667_wp,     & !< inverse of hours per day (1/24)
185                           sigma_sb       = 5.67037321E-8_wp,       & !< Stefan-Boltzmann constant
186                           solar_constant = 1368.0_wp                 !< solar constant at top of atmosphere
187
188    REAL(wp) :: albedo = 9999999.9_wp,           & !< NAMELIST alpha
189                albedo_lw_dif = 9999999.9_wp,    & !< NAMELIST aldif
190                albedo_lw_dir = 9999999.9_wp,    & !< NAMELIST aldir
191                albedo_sw_dif = 9999999.9_wp,    & !< NAMELIST asdif
192                albedo_sw_dir = 9999999.9_wp,    & !< NAMELIST asdir
193                decl_1,                          & !< declination coef. 1
194                decl_2,                          & !< declination coef. 2
195                decl_3,                          & !< declination coef. 3
196                dt_radiation = 0.0_wp,           & !< radiation model timestep
197                emissivity = 0.98_wp,            & !< NAMELIST surface emissivity
198                lambda = 0.0_wp,                 & !< longitude in degrees
199                lon = 0.0_wp,                    & !< longitude in radians
200                lat = 0.0_wp,                    & !< latitude in radians
201                net_radiation = 0.0_wp,          & !< net radiation at surface
202                skip_time_do_radiation = 0.0_wp, & !< Radiation model is not called before this time
203                sky_trans,                       & !< sky transmissivity
204                time_radiation = 0.0_wp,         & !< time since last call of radiation code
205                time_utc,                        & !< current time in UTC
206                time_utc_init = 43200.0_wp         !< UTC time at model start (noon)
207
208    REAL(wp), DIMENSION(0:0) ::  zenith        !< solar zenith angle
209
210    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
211                alpha,                       & !< surface broadband albedo (used for clear-sky scheme)
212                rad_lw_out_change_0,         & !< change in LW out due to change in surface temperature
213                rad_net,                     & !< net radiation at the surface
214                rad_net_av                     !< average of rad_net
215
216!
217!-- Land surface albedos for solar zenith angle of 60° after Briegleb (1992)     
218!-- (shortwave, longwave, broadband):   sw,      lw,      bb,
219    REAL(wp), DIMENSION(0:2,1:17), PARAMETER :: albedo_pars = RESHAPE( (/& 
220                                   0.06_wp, 0.06_wp, 0.06_wp,            & !  1
221                                   0.09_wp, 0.28_wp, 0.19_wp,            & !  2
222                                   0.11_wp, 0.33_wp, 0.23_wp,            & !  3
223                                   0.11_wp, 0.33_wp, 0.23_wp,            & !  4
224                                   0.14_wp, 0.34_wp, 0.25_wp,            & !  5
225                                   0.06_wp, 0.22_wp, 0.14_wp,            & !  6
226                                   0.06_wp, 0.27_wp, 0.17_wp,            & !  7
227                                   0.06_wp, 0.31_wp, 0.19_wp,            & !  8
228                                   0.06_wp, 0.22_wp, 0.14_wp,            & !  9
229                                   0.06_wp, 0.28_wp, 0.18_wp,            & ! 10
230                                   0.35_wp, 0.51_wp, 0.43_wp,            & ! 11
231                                   0.24_wp, 0.40_wp, 0.32_wp,            & ! 12
232                                   0.10_wp, 0.27_wp, 0.19_wp,            & ! 13
233                                   0.90_wp, 0.65_wp, 0.77_wp,            & ! 14
234                                   0.90_wp, 0.65_wp, 0.77_wp,            & ! 15
235                                   0.95_wp, 0.70_wp, 0.82_wp,            & ! 16
236                                   0.08_wp, 0.08_wp, 0.08_wp             & ! 17
237                                 /), (/ 3, 17 /) )
238
239    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: &
240                        rad_lw_cs_hr,                  & !< longwave clear sky radiation heating rate (K/s)
241                        rad_lw_cs_hr_av,               & !< average of rad_lw_cs_hr
242                        rad_lw_hr,                     & !< longwave radiation heating rate (K/s)
243                        rad_lw_hr_av,                  & !< average of rad_sw_hr
244                        rad_lw_in,                     & !< incoming longwave radiation (W/m2)
245                        rad_lw_in_av,                  & !< average of rad_lw_in
246                        rad_lw_out,                    & !< outgoing longwave radiation (W/m2)
247                        rad_lw_out_av,                 & !< average of rad_lw_out
248                        rad_sw_cs_hr,                  & !< shortwave clear sky radiation heating rate (K/s)
249                        rad_sw_cs_hr_av,               & !< average of rad_sw_cs_hr
250                        rad_sw_hr,                     & !< shortwave radiation heating rate (K/s)
251                        rad_sw_hr_av,                  & !< average of rad_sw_hr
252                        rad_sw_in,                     & !< incoming shortwave radiation (W/m2)
253                        rad_sw_in_av,                  & !< average of rad_sw_in
254                        rad_sw_out,                    & !< outgoing shortwave radiation (W/m2)
255                        rad_sw_out_av                    !< average of rad_sw_out
256
257
258!
259!-- Variables and parameters used in RRTMG only
260#if defined ( __rrtmg )
261    CHARACTER(LEN=12) :: rrtm_input_file = "RAD_SND_DATA" !< name of the NetCDF input file (sounding data)
262
263
264!
265!-- Flag parameters for RRTMGS (should not be changed)
266    INTEGER(iwp), PARAMETER :: rrtm_inflglw  = 2, & !< flag for lw cloud optical properties (0,1,2)
267                               rrtm_iceflglw = 0, & !< flag for lw ice particle specifications (0,1,2,3)
268                               rrtm_liqflglw = 1, & !< flag for lw liquid droplet specifications
269                               rrtm_inflgsw  = 2, & !< flag for sw cloud optical properties (0,1,2)
270                               rrtm_iceflgsw = 0, & !< flag for sw ice particle specifications (0,1,2,3)
271                               rrtm_liqflgsw = 1    !< flag for sw liquid droplet specifications
272
273!
274!-- The following variables should be only changed with care, as this will
275!-- require further setting of some variables, which is currently not
276!-- implemented (aerosols, ice phase).
277    INTEGER(iwp) :: nzt_rad,           & !< upper vertical limit for radiation calculations
278                    rrtm_icld = 0,     & !< cloud flag (0: clear sky column, 1: cloudy column)
279                    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)
280                    rrtm_idrv = 1        !< longwave upward flux calculation option (0,1)
281
282    INTEGER(iwp) :: nc_stat !< local variable for storin the result of netCDF calls for error message handling
283
284    LOGICAL :: snd_exists = .FALSE.      !< flag parameter to check whether a user-defined input files exists
285
286    REAL(wp), PARAMETER :: mol_mass_air_d_wv = 1.607793_wp !< molecular weight dry air / water vapor
287
288    REAL(wp), DIMENSION(:), ALLOCATABLE :: hyp_snd,     & !< hypostatic pressure from sounding data (hPa)
289                                           q_snd,       & !< specific humidity from sounding data (kg/kg) - dummy at the moment
290                                           rrtm_tsfc,   & !< dummy array for storing surface temperature
291                                           t_snd          !< actual temperature from sounding data (hPa)
292
293    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: aldif,          & !< longwave diffuse albedo solar angle of 60°
294                                             aldir,          & !< longwave direct albedo solar angle of 60°
295                                             asdif,          & !< shortwave diffuse albedo solar angle of 60°
296                                             asdir,          & !< shortwave direct albedo solar angle of 60°
297                                             rrtm_ccl4vmr,   & !< CCL4 volume mixing ratio (g/mol)
298                                             rrtm_cfc11vmr,  & !< CFC11 volume mixing ratio (g/mol)
299                                             rrtm_cfc12vmr,  & !< CFC12 volume mixing ratio (g/mol)
300                                             rrtm_cfc22vmr,  & !< CFC22 volume mixing ratio (g/mol)
301                                             rrtm_ch4vmr,    & !< CH4 volume mixing ratio
302                                             rrtm_cicewp,    & !< in-cloud ice water path (g/m²)
303                                             rrtm_cldfr,     & !< cloud fraction (0,1)
304                                             rrtm_cliqwp,    & !< in-cloud liquid water path (g/m²)
305                                             rrtm_co2vmr,    & !< CO2 volume mixing ratio (g/mol)
306                                             rrtm_emis,      & !< surface emissivity (0-1)   
307                                             rrtm_h2ovmr,    & !< H2O volume mixing ratio
308                                             rrtm_n2ovmr,    & !< N2O volume mixing ratio
309                                             rrtm_o2vmr,     & !< O2 volume mixing ratio
310                                             rrtm_o3vmr,     & !< O3 volume mixing ratio
311                                             rrtm_play,      & !< pressure layers (hPa, zu-grid)
312                                             rrtm_plev,      & !< pressure layers (hPa, zw-grid)
313                                             rrtm_reice,     & !< cloud ice effective radius (microns)
314                                             rrtm_reliq,     & !< cloud water drop effective radius (microns)
315                                             rrtm_tlay,      & !< actual temperature (K, zu-grid)
316                                             rrtm_tlev,      & !< actual temperature (K, zw-grid)
317                                             rrtm_lwdflx,    & !< RRTM output of incoming longwave radiation flux (W/m2)
318                                             rrtm_lwdflxc,   & !< RRTM output of outgoing clear sky longwave radiation flux (W/m2)
319                                             rrtm_lwuflx,    & !< RRTM output of outgoing longwave radiation flux (W/m2)
320                                             rrtm_lwuflxc,   & !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
321                                             rrtm_lwuflx_dt, & !< RRTM output of incoming clear sky longwave radiation flux (W/m2)
322                                             rrtm_lwuflxc_dt,& !< RRTM output of outgoing clear sky longwave radiation flux (W/m2)
323                                             rrtm_lwhr,      & !< RRTM output of longwave radiation heating rate (K/d)
324                                             rrtm_lwhrc,     & !< RRTM output of incoming longwave clear sky radiation heating rate (K/d)
325                                             rrtm_swdflx,    & !< RRTM output of incoming shortwave radiation flux (W/m2)
326                                             rrtm_swdflxc,   & !< RRTM output of outgoing clear sky shortwave radiation flux (W/m2)
327                                             rrtm_swuflx,    & !< RRTM output of outgoing shortwave radiation flux (W/m2)
328                                             rrtm_swuflxc,   & !< RRTM output of incoming clear sky shortwave radiation flux (W/m2)
329                                             rrtm_swhr,      & !< RRTM output of shortwave radiation heating rate (K/d)
330                                             rrtm_swhrc        !< RRTM output of incoming shortwave clear sky radiation heating rate (K/d)
331
332!
333!-- Definition of arrays that are currently not used for calling RRTMG (due to setting of flag parameters)
334    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  rad_lw_cs_in,   & !< incoming clear sky longwave radiation (W/m2) (not used)
335                                                rad_lw_cs_out,  & !< outgoing clear sky longwave radiation (W/m2) (not used)
336                                                rad_sw_cs_in,   & !< incoming clear sky shortwave radiation (W/m2) (not used)
337                                                rad_sw_cs_out,  & !< outgoing clear sky shortwave radiation (W/m2) (not used)
338                                                rrtm_aldif,     & !< surface albedo for longwave diffuse radiation
339                                                rrtm_aldir,     & !< surface albedo for longwave direct radiation
340                                                rrtm_asdif,     & !< surface albedo for shortwave diffuse radiation
341                                                rrtm_asdir,     & !< surface albedo for shortwave direct radiation
342                                                rrtm_lw_tauaer, & !< lw aerosol optical depth
343                                                rrtm_lw_taucld, & !< lw in-cloud optical depth
344                                                rrtm_sw_taucld, & !< sw in-cloud optical depth
345                                                rrtm_sw_ssacld, & !< sw in-cloud single scattering albedo
346                                                rrtm_sw_asmcld, & !< sw in-cloud asymmetry parameter
347                                                rrtm_sw_fsfcld, & !< sw in-cloud forward scattering fraction
348                                                rrtm_sw_tauaer, & !< sw aerosol optical depth
349                                                rrtm_sw_ssaaer, & !< sw aerosol single scattering albedo
350                                                rrtm_sw_asmaer, & !< sw aerosol asymmetry parameter
351                                                rrtm_sw_ecaer     !< sw aerosol optical detph at 0.55 microns (rrtm_iaer = 6 only)
352
353#endif
354
355    INTERFACE radiation_check_data_output
356       MODULE PROCEDURE radiation_check_data_output
357    END INTERFACE radiation_check_data_output
358
359    INTERFACE radiation_check_data_output_pr
360       MODULE PROCEDURE radiation_check_data_output_pr
361    END INTERFACE radiation_check_data_output_pr
362 
363    INTERFACE radiation_check_parameters
364       MODULE PROCEDURE radiation_check_parameters
365    END INTERFACE radiation_check_parameters
366 
367    INTERFACE radiation_clearsky
368       MODULE PROCEDURE radiation_clearsky
369    END INTERFACE radiation_clearsky
370 
371    INTERFACE radiation_header
372       MODULE PROCEDURE radiation_header
373    END INTERFACE radiation_header
374 
375    INTERFACE radiation_init
376       MODULE PROCEDURE radiation_init
377    END INTERFACE radiation_init
378
379    INTERFACE radiation_parin
380       MODULE PROCEDURE radiation_parin
381    END INTERFACE radiation_parin
382   
383    INTERFACE radiation_rrtmg
384       MODULE PROCEDURE radiation_rrtmg
385    END INTERFACE radiation_rrtmg
386
387    INTERFACE radiation_tendency
388       MODULE PROCEDURE radiation_tendency
389       MODULE PROCEDURE radiation_tendency_ij
390    END INTERFACE radiation_tendency
391
392    SAVE
393
394    PRIVATE
395
396!
397!-- Public functions
398    PUBLIC radiation_check_data_output, radiation_check_data_output_pr,        &
399           radiation_check_parameters, radiation_clearsky, radiation_header,   &
400           radiation_init, radiation_parin, radiation_rrtmg, radiation_tendency
401   
402!
403!-- Public variables and constants
404    PUBLIC dots_rad, dt_radiation, force_radiation_call,                       &
405           rad_net, rad_net_av, radiation, radiation_scheme, rad_lw_in,        &
406           rad_lw_in_av, rad_lw_out, rad_lw_out_av, rad_lw_out_change_0,       &
407           rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av, rad_sw_in,  &
408           rad_sw_in_av, rad_sw_out, rad_sw_out_av, rad_sw_cs_hr,              &
409           rad_sw_cs_hr_av, rad_sw_hr, rad_sw_hr_av, sigma_sb,                 &
410           skip_time_do_radiation, time_radiation, unscheduled_radiation_calls
411
412
413#if defined ( __rrtmg )
414    PUBLIC rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir, rrtm_idrv
415#endif
416
417 CONTAINS
418
419!------------------------------------------------------------------------------!
420! Description:
421! ------------
422!> Check data output for radiation model
423!------------------------------------------------------------------------------!
424    SUBROUTINE radiation_check_data_output( var, unit, i, ilen, k )
425 
426 
427       USE control_parameters,                                                 &
428           ONLY: data_output, message_string
429
430       IMPLICIT NONE
431
432       CHARACTER (LEN=*) ::  unit     !<
433       CHARACTER (LEN=*) ::  var      !<
434
435       INTEGER(iwp) :: i
436       INTEGER(iwp) :: ilen
437       INTEGER(iwp) :: k
438
439       SELECT CASE ( TRIM( var ) )
440
441         CASE ( 'rad_lw_in', 'rad_lw_out', 'rad_lw_cs_hr', 'rad_lw_hr',       &
442                 'rad_sw_in', 'rad_sw_out', 'rad_sw_cs_hr', 'rad_sw_hr' )
443             IF (  .NOT.  radiation  .OR.  radiation_scheme /= 'rrtmg' )  THEN
444                message_string = '"output of "' // TRIM( var ) // '" requi' // &
445                                 'res radiation = .TRUE. and ' //              &
446                                 'radiation_scheme = "rrtmg"'
447                CALL message( 'check_parameters', 'PA0406', 1, 2, 0, 6, 0 )
448             ENDIF
449             unit = 'W/m2'     
450
451          CASE ( 'rad_net*', 'rrtm_aldif*', 'rrtm_aldir*', 'rrtm_asdif*',      &
452                 'rrtm_asdir*' )
453             IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
454                message_string = 'illegal value for data_output: "' //         &
455                                 TRIM( var ) // '" & only 2d-horizontal ' //   &
456                                 'cross sections are allowed for this value'
457                CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
458             ENDIF
459             IF (  .NOT.  radiation  .OR.  radiation_scheme /= "rrtmg" )  THEN
460                IF ( TRIM( var ) == 'rrtm_aldif*'  .OR.                        &
461                     TRIM( var ) == 'rrtm_aldir*'  .OR.                        &
462                     TRIM( var ) == 'rrtm_asdif*'  .OR.                        &
463                     TRIM( var ) == 'rrtm_asdir*'      )                       &
464                THEN
465                   message_string = 'output of "' // TRIM( var ) // '" require'&
466                                    // 's radiation = .TRUE. and radiation_sch'&
467                                    // 'eme = "rrtmg"'
468                   CALL message( 'check_parameters', 'PA0409', 1, 2, 0, 6, 0 )
469                ENDIF
470             ENDIF
471
472             IF ( TRIM( var ) == 'rad_net*'      ) unit = 'W/m2'   
473             IF ( TRIM( var ) == 'rrtm_aldif*'   ) unit = ''   
474             IF ( TRIM( var ) == 'rrtm_aldir*'   ) unit = '' 
475             IF ( TRIM( var ) == 'rrtm_asdif*'   ) unit = '' 
476             IF ( TRIM( var ) == 'rrtm_asdir*'   ) unit = '' 
477
478          CASE DEFAULT
479             unit = 'illegal'
480
481       END SELECT
482
483
484    END SUBROUTINE radiation_check_data_output
485
486!------------------------------------------------------------------------------!
487! Description:
488! ------------
489!> Check data output of profiles for radiation model
490!------------------------------------------------------------------------------! 
491    SUBROUTINE radiation_check_data_output_pr( variable, var_count, unit, dopr_unit )
492 
493       USE arrays_3d,                                                          &
494           ONLY: zu
495
496       USE control_parameters,                                                 &
497           ONLY: data_output_pr, message_string
498
499       USE indices
500
501       USE profil_parameter
502
503       USE statistics
504
505       IMPLICIT NONE
506   
507       CHARACTER (LEN=*) ::  unit      !<
508       CHARACTER (LEN=*) ::  variable  !<
509       CHARACTER (LEN=*) ::  dopr_unit !< local value of dopr_unit
510 
511       INTEGER(iwp) ::  user_pr_index !<
512       INTEGER(iwp) ::  var_count     !<
513
514       SELECT CASE ( TRIM( variable ) )
515       
516         CASE ( 'rad_net' )
517             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme == 'constant' )&
518             THEN
519                message_string = 'data_output_pr = ' //                        &
520                                 TRIM( data_output_pr(var_count) ) // ' is' // &
521                                 'not available for radiation = .FALSE. or ' //&
522                                 'radiation_scheme = "constant"'
523                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
524             ELSE
525                dopr_index(var_count) = 101
526                dopr_unit  = 'W/m2'
527                hom(:,2,101,:)  = SPREAD( zw, 2, statistic_regions+1 )
528                unit = dopr_unit
529             ENDIF
530
531          CASE ( 'rad_lw_in' )
532             IF ( (  .NOT.  radiation)  .OR.  radiation_scheme == 'constant' ) &
533             THEN
534                message_string = 'data_output_pr = ' //                        &
535                                 TRIM( data_output_pr(var_count) ) // ' is' // &
536                                 'not available for radiation = .FALSE. or ' //&
537                                 'radiation_scheme = "constant"'
538                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
539             ELSE
540                dopr_index(var_count) = 102
541                dopr_unit  = 'W/m2'
542                hom(:,2,102,:)  = SPREAD( zw, 2, statistic_regions+1 )
543                unit = dopr_unit
544             ENDIF
545
546          CASE ( 'rad_lw_out' )
547             IF ( (  .NOT. radiation )  .OR.  radiation_scheme == 'constant' ) &
548             THEN
549                message_string = 'data_output_pr = ' //                        &
550                                 TRIM( data_output_pr(var_count) ) // ' is' // &
551                                 'not available for radiation = .FALSE. or ' //&
552                                 'radiation_scheme = "constant"'
553                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
554             ELSE
555                dopr_index(var_count) = 103
556                dopr_unit  = 'W/m2'
557                hom(:,2,103,:)  = SPREAD( zw, 2, statistic_regions+1 )
558                unit = dopr_unit   
559             ENDIF
560
561          CASE ( 'rad_sw_in' )
562             IF ( (  .NOT. radiation )  .OR.  radiation_scheme == 'constant' ) &
563             THEN
564                message_string = 'data_output_pr = ' //                        &
565                                 TRIM( data_output_pr(var_count) ) // ' is' // &
566                                 'not available for radiation = .FALSE. or ' //&
567                                 'radiation_scheme = "constant"'
568                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
569             ELSE
570                dopr_index(var_count) = 104
571                dopr_unit  = 'W/m2'
572                hom(:,2,104,:)  = SPREAD( zw, 2, statistic_regions+1 )
573                unit = dopr_unit
574             ENDIF
575
576          CASE ( 'rad_sw_out')
577             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme == 'constant' )&
578             THEN
579                message_string = 'data_output_pr = ' //                        &
580                                 TRIM( data_output_pr(var_count) ) // ' is' // &
581                                 'not available for radiation = .FALSE. or ' //&
582                                 'radiation_scheme = "constant"'
583                CALL message( 'check_parameters', 'PA0408', 1, 2, 0, 6, 0 )
584             ELSE
585                dopr_index(var_count) = 105
586                dopr_unit  = 'W/m2'
587                hom(:,2,105,:)  = SPREAD( zw, 2, statistic_regions+1 )
588                unit = dopr_unit
589             ENDIF
590
591          CASE ( 'rad_lw_cs_hr' )
592             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
593             THEN
594                message_string = 'data_output_pr = ' //                        &
595                                 TRIM( data_output_pr(var_count) ) // ' is' // &
596                                 'not available for radiation = .FALSE. or ' //&
597                                 'radiation_scheme /= "rrtmg"'
598                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
599             ELSE
600                dopr_index(var_count) = 106
601                dopr_unit  = 'K/h'
602                hom(:,2,106,:)  = SPREAD( zu, 2, statistic_regions+1 )
603                unit = dopr_unit
604             ENDIF
605
606          CASE ( 'rad_lw_hr' )
607             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
608             THEN
609                message_string = 'data_output_pr = ' //                        &
610                                 TRIM( data_output_pr(var_count) ) // ' is' // &
611                                 'not available for radiation = .FALSE. or ' //&
612                                 'radiation_scheme /= "rrtmg"'
613                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
614             ELSE
615                dopr_index(var_count) = 107
616                dopr_unit  = 'K/h'
617                hom(:,2,107,:)  = SPREAD( zu, 2, statistic_regions+1 )
618                unit = dopr_unit
619             ENDIF
620
621          CASE ( 'rad_sw_cs_hr' )
622             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
623             THEN
624                message_string = 'data_output_pr = ' //                        &
625                                 TRIM( data_output_pr(var_count) ) // ' is' // &
626                                 'not available for radiation = .FALSE. or ' //&
627                                 'radiation_scheme /= "rrtmg"'
628                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
629             ELSE
630                dopr_index(var_count) = 108
631                dopr_unit  = 'K/h'
632                hom(:,2,108,:)  = SPREAD( zu, 2, statistic_regions+1 )
633                unit = dopr_unit
634             ENDIF
635
636          CASE ( 'rad_sw_hr' )
637             IF ( (  .NOT.  radiation )  .OR.  radiation_scheme /= 'rrtmg' )   &
638             THEN
639                message_string = 'data_output_pr = ' //                        &
640                                 TRIM( data_output_pr(var_count) ) // ' is' // &
641                                 'not available for radiation = .FALSE. or ' //&
642                                 'radiation_scheme /= "rrtmg"'
643                CALL message( 'check_parameters', 'PA0413', 1, 2, 0, 6, 0 )
644             ELSE
645                dopr_index(var_count) = 109
646                dopr_unit  = 'K/h'
647                hom(:,2,109,:)  = SPREAD( zu, 2, statistic_regions+1 )
648                unit = dopr_unit
649             ENDIF
650
651
652          CASE DEFAULT
653             unit = 'illegal'
654
655       END SELECT
656
657
658    END SUBROUTINE radiation_check_data_output_pr
659 
660 
661!------------------------------------------------------------------------------!
662! Description:
663! ------------
664!> Check parameters routine for radiation model
665!------------------------------------------------------------------------------!
666    SUBROUTINE radiation_check_parameters
667
668       USE control_parameters,                                                 &
669           ONLY: message_string, topography
670                 
671   
672       IMPLICIT NONE
673
674       IF ( radiation_scheme /= 'constant'   .AND.                             &
675            radiation_scheme /= 'clear-sky'  .AND.                             &
676            radiation_scheme /= 'rrtmg' )  THEN
677          message_string = 'unknown radiation_scheme = '//                     &
678                           TRIM( radiation_scheme )
679          CALL message( 'check_parameters', 'PA0405', 1, 2, 0, 6, 0 )
680       ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
681#if ! defined ( __rrtmg )
682          message_string = 'radiation_scheme = "rrtmg" requires ' //           & 
683                           'compilation of PALM with pre-processor ' //        &
684                           'directive -D__rrtmg'
685          CALL message( 'check_parameters', 'PA0407', 1, 2, 0, 6, 0 )
686#endif
687#if defined ( __rrtmg ) && ! defined( __netcdf )
688          message_string = 'radiation_scheme = "rrtmg" requires ' //           & 
689                           'the use of NetCDF (preprocessor directive ' //     &
690                           '-D__netcdf'
691          CALL message( 'check_parameters', 'PA0412', 1, 2, 0, 6, 0 )
692#endif
693
694       ENDIF
695
696       IF ( albedo_type == 0  .AND.  albedo == 9999999.9_wp  .AND.             &
697            radiation_scheme == 'clear-sky')  THEN
698          message_string = 'radiation_scheme = "clear-sky" in combination' //  & 
699                           'with albedo_type = 0 requires setting of albedo'// &
700                           ' /= 9999999.9'
701          CALL message( 'check_parameters', 'PA0410', 1, 2, 0, 6, 0 )
702       ENDIF
703
704       IF ( albedo_type == 0  .AND.  radiation_scheme == 'rrtmg'  .AND.        &
705          (    albedo_lw_dif == 9999999.9_wp .OR. albedo_lw_dir == 9999999.9_wp&
706          .OR. albedo_sw_dif == 9999999.9_wp .OR. albedo_sw_dir == 9999999.9_wp& 
707          ) ) THEN
708          message_string = 'radiation_scheme = "rrtmg" in combination' //      & 
709                           'with albedo_type = 0 requires setting of ' //      &
710                           'albedo_lw_dif /= 9999999.9' //                     &
711                           'albedo_lw_dir /= 9999999.9' //                     &
712                           'albedo_sw_dif /= 9999999.9 and' //                 &
713                           'albedo_sw_dir /= 9999999.9'
714          CALL message( 'check_parameters', 'PA0411', 1, 2, 0, 6, 0 )
715       ENDIF
716
717       IF ( topography /= 'flat' )  THEN
718          message_string = 'radiation scheme cannot be used ' //               & 
719                           'in combination with  topography /= "flat"'
720          CALL message( 'check_parameters', 'PA0414', 1, 2, 0, 6, 0 )
721       ENDIF
722 
723    END SUBROUTINE radiation_check_parameters
724 
725 
726!------------------------------------------------------------------------------!
727! Description:
728! ------------
729!> Initialization of the radiation model
730!------------------------------------------------------------------------------!
731    SUBROUTINE radiation_init
732   
733       IMPLICIT NONE
734
735!
736!--    Allocate array for storing the surface net radiation
737       IF ( .NOT. ALLOCATED ( rad_net ) )  THEN
738          ALLOCATE ( rad_net(nysg:nyng,nxlg:nxrg) )
739          rad_net = 0.0_wp
740       ENDIF
741
742!
743!--    Allocate array for storing the surface net radiation
744       IF ( .NOT. ALLOCATED ( rad_lw_out_change_0 ) )  THEN
745          ALLOCATE ( rad_lw_out_change_0(nysg:nyng,nxlg:nxrg) )
746          rad_lw_out_change_0 = 0.0_wp
747       ENDIF
748
749!
750!--    Fix net radiation in case of radiation_scheme = 'constant'
751       IF ( radiation_scheme == 'constant' )  THEN
752          rad_net = net_radiation
753          radiation = .FALSE.
754!
755!--    Calculate orbital constants
756       ELSE
757          decl_1 = SIN(23.45_wp * pi / 180.0_wp)
758          decl_2 = 2.0_wp * pi / 365.0_wp
759          decl_3 = decl_2 * 81.0_wp
760          lat    = phi * pi / 180.0_wp
761          lon    = lambda * pi / 180.0_wp
762       ENDIF
763
764
765       IF ( radiation_scheme == 'clear-sky' )  THEN
766
767          ALLOCATE ( alpha(nysg:nyng,nxlg:nxrg) )
768
769          IF ( .NOT. ALLOCATED ( rad_sw_in ) )  THEN
770             ALLOCATE ( rad_sw_in(0:0,nysg:nyng,nxlg:nxrg) )
771          ENDIF
772          IF ( .NOT. ALLOCATED ( rad_sw_out ) )  THEN
773             ALLOCATE ( rad_sw_out(0:0,nysg:nyng,nxlg:nxrg) )
774          ENDIF
775
776          IF ( .NOT. ALLOCATED ( rad_sw_in_av ) )  THEN
777             ALLOCATE ( rad_sw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
778          ENDIF
779          IF ( .NOT. ALLOCATED ( rad_sw_out_av ) )  THEN
780             ALLOCATE ( rad_sw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
781          ENDIF
782
783          IF ( .NOT. ALLOCATED ( rad_lw_in ) )  THEN
784             ALLOCATE ( rad_lw_in(0:0,nysg:nyng,nxlg:nxrg) )
785          ENDIF
786          IF ( .NOT. ALLOCATED ( rad_lw_out ) )  THEN
787             ALLOCATE ( rad_lw_out(0:0,nysg:nyng,nxlg:nxrg) )
788          ENDIF
789
790          IF ( .NOT. ALLOCATED ( rad_lw_in_av ) )  THEN
791             ALLOCATE ( rad_lw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
792          ENDIF
793          IF ( .NOT. ALLOCATED ( rad_lw_out_av ) )  THEN
794             ALLOCATE ( rad_lw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
795          ENDIF
796
797          rad_sw_in  = 0.0_wp
798          rad_sw_out = 0.0_wp
799          rad_lw_in  = 0.0_wp
800          rad_lw_out = 0.0_wp
801
802!
803!--       Overwrite albedo if manually set in parameter file
804          IF ( albedo_type /= 0 .AND. albedo == 9999999.9_wp )  THEN
805             albedo = albedo_pars(2,albedo_type)
806          ENDIF
807   
808          alpha = albedo
809 
810!
811!--    Initialization actions for RRTMG
812       ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
813#if defined ( __rrtmg )
814!
815!--       Allocate albedos
816          ALLOCATE ( rrtm_aldif(0:0,nysg:nyng,nxlg:nxrg) )
817          ALLOCATE ( rrtm_aldir(0:0,nysg:nyng,nxlg:nxrg) )
818          ALLOCATE ( rrtm_asdif(0:0,nysg:nyng,nxlg:nxrg) )
819          ALLOCATE ( rrtm_asdir(0:0,nysg:nyng,nxlg:nxrg) )
820          ALLOCATE ( aldif(nysg:nyng,nxlg:nxrg) )
821          ALLOCATE ( aldir(nysg:nyng,nxlg:nxrg) )
822          ALLOCATE ( asdif(nysg:nyng,nxlg:nxrg) )
823          ALLOCATE ( asdir(nysg:nyng,nxlg:nxrg) )
824
825          IF ( albedo_type /= 0 )  THEN
826             IF ( albedo_lw_dif == 9999999.9_wp )  THEN
827                albedo_lw_dif = albedo_pars(0,albedo_type)
828                albedo_lw_dir = albedo_lw_dif
829             ENDIF
830             IF ( albedo_sw_dif == 9999999.9_wp )  THEN
831                albedo_sw_dif = albedo_pars(1,albedo_type)
832                albedo_sw_dir = albedo_sw_dif
833             ENDIF
834          ENDIF
835
836          aldif(:,:) = albedo_lw_dif
837          aldir(:,:) = albedo_lw_dir
838          asdif(:,:) = albedo_sw_dif
839          asdir(:,:) = albedo_sw_dir
840!
841!--       Calculate initial values of current (cosine of) the zenith angle and
842!--       whether the sun is up
843          CALL calc_zenith     
844!
845!--       Calculate initial surface albedo
846          IF ( .NOT. constant_albedo )  THEN
847             CALL calc_albedo
848          ELSE
849             rrtm_aldif(0,:,:) = aldif(:,:)
850             rrtm_aldir(0,:,:) = aldir(:,:)
851             rrtm_asdif(0,:,:) = asdif(:,:) 
852             rrtm_asdir(0,:,:) = asdir(:,:)   
853          ENDIF
854
855!
856!--       Allocate surface emissivity
857          ALLOCATE ( rrtm_emis(0:0,1:nbndlw+1) )
858          rrtm_emis = emissivity
859
860!
861!--       Allocate 3d arrays of radiative fluxes and heating rates
862          IF ( .NOT. ALLOCATED ( rad_sw_in ) )  THEN
863             ALLOCATE ( rad_sw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
864             rad_sw_in = 0.0_wp
865          ENDIF
866
867          IF ( .NOT. ALLOCATED ( rad_sw_in_av ) )  THEN
868             ALLOCATE ( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
869          ENDIF
870
871          IF ( .NOT. ALLOCATED ( rad_sw_out ) )  THEN
872             ALLOCATE ( rad_sw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
873             rad_sw_out = 0.0_wp
874          ENDIF
875
876          IF ( .NOT. ALLOCATED ( rad_sw_out_av ) )  THEN
877             ALLOCATE ( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
878          ENDIF
879
880          IF ( .NOT. ALLOCATED ( rad_sw_hr ) )  THEN
881             ALLOCATE ( rad_sw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
882             rad_sw_hr = 0.0_wp
883          ENDIF
884
885          IF ( .NOT. ALLOCATED ( rad_sw_hr_av ) )  THEN
886             ALLOCATE ( rad_sw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
887             rad_sw_hr_av = 0.0_wp
888          ENDIF
889
890          IF ( .NOT. ALLOCATED ( rad_sw_cs_hr ) )  THEN
891             ALLOCATE ( rad_sw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
892             rad_sw_cs_hr = 0.0_wp
893          ENDIF
894
895          IF ( .NOT. ALLOCATED ( rad_sw_cs_hr_av ) )  THEN
896             ALLOCATE ( rad_sw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
897             rad_sw_cs_hr_av = 0.0_wp
898          ENDIF
899
900          IF ( .NOT. ALLOCATED ( rad_lw_in ) )  THEN
901             ALLOCATE ( rad_lw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
902             rad_lw_in     = 0.0_wp
903          ENDIF
904
905          IF ( .NOT. ALLOCATED ( rad_lw_in_av ) )  THEN
906             ALLOCATE ( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
907          ENDIF
908
909          IF ( .NOT. ALLOCATED ( rad_lw_out ) )  THEN
910             ALLOCATE ( rad_lw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
911            rad_lw_out    = 0.0_wp
912          ENDIF
913
914          IF ( .NOT. ALLOCATED ( rad_lw_out_av ) )  THEN
915             ALLOCATE ( rad_lw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
916          ENDIF
917
918          IF ( .NOT. ALLOCATED ( rad_lw_hr ) )  THEN
919             ALLOCATE ( rad_lw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
920             rad_lw_hr = 0.0_wp
921          ENDIF
922
923          IF ( .NOT. ALLOCATED ( rad_lw_hr_av ) )  THEN
924             ALLOCATE ( rad_lw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
925             rad_lw_hr_av = 0.0_wp
926          ENDIF
927
928          IF ( .NOT. ALLOCATED ( rad_lw_cs_hr ) )  THEN
929             ALLOCATE ( rad_lw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
930             rad_lw_cs_hr = 0.0_wp
931          ENDIF
932
933          IF ( .NOT. ALLOCATED ( rad_lw_cs_hr_av ) )  THEN
934             ALLOCATE ( rad_lw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
935             rad_lw_cs_hr_av = 0.0_wp
936          ENDIF
937
938          ALLOCATE ( rad_sw_cs_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
939          ALLOCATE ( rad_sw_cs_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
940          rad_sw_cs_in  = 0.0_wp
941          rad_sw_cs_out = 0.0_wp
942
943          ALLOCATE ( rad_lw_cs_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
944          ALLOCATE ( rad_lw_cs_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
945          rad_lw_cs_in  = 0.0_wp
946          rad_lw_cs_out = 0.0_wp
947
948!
949!--       Allocate dummy array for storing surface temperature
950          ALLOCATE ( rrtm_tsfc(1) )
951
952!
953!--       Initialize RRTMG
954          IF ( lw_radiation )  CALL rrtmg_lw_ini ( cp )
955          IF ( sw_radiation )  CALL rrtmg_sw_ini ( cp )
956
957!
958!--       Set input files for RRTMG
959          INQUIRE(FILE="RAD_SND_DATA", EXIST=snd_exists) 
960          IF ( .NOT. snd_exists )  THEN
961             rrtm_input_file = "rrtmg_lw.nc"
962          ENDIF
963
964!
965!--       Read vertical layers for RRTMG from sounding data
966!--       The routine provides nzt_rad, hyp_snd(1:nzt_rad),
967!--       t_snd(nzt+2:nzt_rad), rrtm_play(1:nzt_rad), rrtm_plev(1_nzt_rad+1),
968!--       rrtm_tlay(nzt+2:nzt_rad), rrtm_tlev(nzt+2:nzt_rad+1)
969          CALL read_sounding_data
970
971!
972!--       Read trace gas profiles from file. This routine provides
973!--       the rrtm_ arrays (1:nzt_rad+1)
974          CALL read_trace_gas_data
975#endif
976       ENDIF
977
978!
979!--    Perform user actions if required
980       CALL user_init_radiation
981
982!
983!--    Calculate radiative fluxes at model start
984       IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
985          IF ( radiation_scheme == 'clear-sky' )  THEN
986             CALL radiation_clearsky
987          ELSEIF ( radiation_scheme == 'rrtmg' )  THEN
988             CALL radiation_rrtmg
989          ENDIF
990       ENDIF
991
992       RETURN
993
994    END SUBROUTINE radiation_init
995
996
997!------------------------------------------------------------------------------!
998! Description:
999! ------------
1000!> A simple clear sky radiation model
1001!------------------------------------------------------------------------------!
1002    SUBROUTINE radiation_clearsky
1003
1004       USE indices,                                                            &
1005           ONLY:  nbgp
1006
1007       IMPLICIT NONE
1008
1009       INTEGER(iwp) :: i, j, k   !< loop indices
1010       REAL(wp)     :: exn,   &  !< Exner functions at surface
1011                       exn1,  &  !< Exner functions at first grid level
1012                       pt1       !< potential temperature at first grid level
1013
1014!
1015!--    Calculate current zenith angle
1016       CALL calc_zenith
1017
1018!
1019!--    Calculate sky transmissivity
1020       sky_trans = 0.6_wp + 0.2_wp * zenith(0)
1021
1022!
1023!--    Calculate value of the Exner function
1024       exn = (surface_pressure / 1000.0_wp )**0.286_wp
1025!
1026!--    Calculate radiation fluxes and net radiation (rad_net) for each grid
1027!--    point
1028       DO i = nxlg, nxrg
1029          DO j = nysg, nyng
1030             k = nzb_s_inner(j,i)
1031
1032             exn1 = (hyp(k+1) / 100000.0_wp )**0.286_wp
1033
1034             rad_sw_in(0,j,i)  = solar_constant * sky_trans * zenith(0)
1035             rad_sw_out(0,j,i) = alpha(j,i) * rad_sw_in(0,j,i)
1036             rad_lw_out(0,j,i) = emissivity * sigma_sb * (pt(k,j,i) * exn)**4
1037
1038             IF ( cloud_physics )  THEN
1039                pt1 = pt(k+1,j,i) + l_d_cp / exn1 * ql(k+1,j,i)
1040                rad_lw_in(0,j,i)  = 0.8_wp * sigma_sb * (pt1 * exn1)**4
1041             ELSE
1042                rad_lw_in(0,j,i)  = 0.8_wp * sigma_sb * (pt(k+1,j,i) * exn1)**4
1043             ENDIF
1044
1045             rad_net(j,i) = rad_sw_in(0,j,i) - rad_sw_out(0,j,i)               &
1046                            + rad_lw_in(0,j,i) - rad_lw_out(0,j,i)
1047
1048          ENDDO
1049       ENDDO
1050
1051    END SUBROUTINE radiation_clearsky
1052
1053
1054!------------------------------------------------------------------------------!
1055! Description:
1056! ------------
1057!> Header output for radiation model
1058!------------------------------------------------------------------------------!
1059    SUBROUTINE radiation_header ( io )
1060
1061
1062       IMPLICIT NONE
1063 
1064       INTEGER(iwp), INTENT(IN) ::  io            !< Unit of the output file
1065   
1066
1067       
1068!
1069!--    Write radiation model header
1070       WRITE( io, 3 )
1071
1072       IF ( radiation_scheme == "constant" )  THEN
1073          WRITE( io, 4 ) net_radiation
1074       ELSEIF ( radiation_scheme == "clear-sky" )  THEN
1075          WRITE( io, 5 )
1076       ELSEIF ( radiation_scheme == "rrtmg" )  THEN
1077          WRITE( io, 6 )
1078          IF ( .NOT. lw_radiation )  WRITE( io, 10 )
1079          IF ( .NOT. sw_radiation )  WRITE( io, 11 )
1080       ENDIF
1081
1082       IF ( albedo_type == 0 )  THEN
1083          WRITE( io, 7 ) albedo
1084       ELSE
1085          WRITE( io, 8 ) TRIM( albedo_type_name(albedo_type) )
1086       ENDIF
1087       IF ( constant_albedo )  THEN
1088          WRITE( io, 9 )
1089       ENDIF
1090       
1091       IF ( radiation .AND. radiation_scheme /= 'constant' )  THEN
1092          WRITE ( io, 1 )  lambda
1093          WRITE ( io, 2 )  day_init, time_utc_init
1094       ENDIF
1095
1096       WRITE( io, 12 ) dt_radiation
1097 
1098
1099 1 FORMAT ('    Geograph. longitude            :   lambda = ',F4.1,' degr')
1100 2 FORMAT ('    Day of the year at model start :   day_init = ',I3   &
1101            /'    UTC time at model start        :   time_utc_init = ',F7.1' s')
1102 3 FORMAT (//' Radiation model information:'/                                  &
1103              ' ----------------------------'/)
1104 4 FORMAT ('    --> Using constant net radiation: net_radiation = ', F6.2,        &
1105           // 'W/m**2')
1106 5 FORMAT ('    --> Simple radiation scheme for clear sky is used (no clouds,',   &
1107                   ' default)')
1108 6 FORMAT ('    --> RRTMG scheme is used')
1109 7 FORMAT (/'    User-specific surface albedo: albedo =', F6.3)
1110 8 FORMAT (/'    Albedo is set for land surface type: ', A)
1111 9 FORMAT (/'    --> Albedo is fixed during the run')
111210 FORMAT (/'    --> Longwave radiation is disabled')
111311 FORMAT (/'    --> Shortwave radiation is disabled.')
111412 FORMAT  ('    Timestep: dt_radiation = ', F6.2, '  s')
1115
1116
1117    END SUBROUTINE radiation_header
1118   
1119
1120!------------------------------------------------------------------------------!
1121! Description:
1122! ------------
1123!> Parin for &radiation_par for radiation model
1124!------------------------------------------------------------------------------!
1125    SUBROUTINE radiation_parin
1126
1127
1128       IMPLICIT NONE
1129
1130       CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file
1131       
1132       NAMELIST /radiation_par/   albedo, albedo_type, albedo_lw_dir,          &
1133                                  albedo_lw_dif, albedo_sw_dir, albedo_sw_dif, &
1134                                  constant_albedo, day_init, dt_radiation,     &
1135                                  lambda, lw_radiation, net_radiation,         &
1136                                  radiation_scheme, skip_time_do_radiation,    &
1137                                  sw_radiation, time_utc_init,                 &
1138                                  unscheduled_radiation_calls
1139       
1140       line = ' '
1141       
1142!
1143!--    Try to find radiation model package
1144       REWIND ( 11 )
1145       line = ' '
1146       DO   WHILE ( INDEX( line, '&radiation_par' ) == 0 )
1147          READ ( 11, '(A)', END=10 )  line
1148       ENDDO
1149       BACKSPACE ( 11 )
1150
1151!
1152!--    Read user-defined namelist
1153       READ ( 11, radiation_par )
1154
1155!
1156!--    Set flag that indicates that the radiation model is switched on
1157       radiation = .TRUE.
1158
1159 10    CONTINUE
1160       
1161
1162    END SUBROUTINE radiation_parin
1163
1164
1165!------------------------------------------------------------------------------!
1166! Description:
1167! ------------
1168!> Implementation of the RRTMG radiation_scheme
1169!------------------------------------------------------------------------------!
1170    SUBROUTINE radiation_rrtmg
1171
1172       USE indices,                                                            &
1173           ONLY:  nbgp
1174
1175       USE particle_attributes,                                                &
1176           ONLY:  grid_particles, number_of_particles, particles,              &
1177                  particle_advection_start, prt_count
1178
1179       IMPLICIT NONE
1180
1181#if defined ( __rrtmg )
1182
1183       INTEGER(iwp) :: i, j, k, n !< loop indices
1184
1185       REAL(wp)     ::  s_r2, &   !< weighted sum over all droplets with r^2
1186                        s_r3      !< weighted sum over all droplets with r^3
1187
1188!
1189!--    Calculate current (cosine of) zenith angle and whether the sun is up
1190       CALL calc_zenith     
1191!
1192!--    Calculate surface albedo
1193       IF ( .NOT. constant_albedo )  THEN
1194          CALL calc_albedo
1195       ENDIF
1196
1197!
1198!--    Prepare input data for RRTMG
1199
1200!
1201!--    In case of large scale forcing with surface data, calculate new pressure
1202!--    profile. nzt_rad might be modified by these calls and all required arrays
1203!--    will then be re-allocated
1204       IF ( large_scale_forcing  .AND.  lsf_surf )  THEN
1205          CALL read_sounding_data
1206          CALL read_trace_gas_data
1207       ENDIF
1208!
1209!--    Loop over all grid points
1210       DO i = nxl, nxr
1211          DO j = nys, nyn
1212
1213!
1214!--          Prepare profiles of temperature and H2O volume mixing ratio
1215             rrtm_tlev(0,nzb+1) = pt(nzb,j,i) * ( surface_pressure             &
1216                                                  / 1000.0_wp )**0.286_wp
1217
1218             DO k = nzb+1, nzt+1
1219                rrtm_tlay(0,k) = pt(k,j,i) * ( (hyp(k) ) / 100000.0_wp         &
1220                                 )**0.286_wp + l_d_cp * ql(k,j,i)
1221                rrtm_h2ovmr(0,k) = mol_mass_air_d_wv * (q(k,j,i) - ql(k,j,i))
1222
1223             ENDDO
1224
1225!
1226!--          Avoid temperature/humidity jumps at the top of the LES domain by
1227!--          linear interpolation from nzt+2 to nzt+7
1228             DO k = nzt+2, nzt+7
1229                rrtm_tlay(0,k) = rrtm_tlay(0,nzt+1)                            &
1230                              + ( rrtm_tlay(0,nzt+8) - rrtm_tlay(0,nzt+1) )    &
1231                              / ( rrtm_play(0,nzt+8) - rrtm_play(0,nzt+1) )    &
1232                              * ( rrtm_play(0,k) - rrtm_play(0,nzt+1) )
1233
1234                rrtm_h2ovmr(0,k) = rrtm_h2ovmr(0,nzt+1)                        &
1235                              + ( rrtm_h2ovmr(0,nzt+8) - rrtm_h2ovmr(0,nzt+1) )&
1236                              / ( rrtm_play(0,nzt+8)   - rrtm_play(0,nzt+1)   )&
1237                              * ( rrtm_play(0,k) - rrtm_play(0,nzt+1) )
1238
1239             ENDDO
1240
1241!--          Linear interpolate to zw grid
1242             DO k = nzb+2, nzt+8
1243                rrtm_tlev(0,k)   = rrtm_tlay(0,k-1) + (rrtm_tlay(0,k) -        &
1244                                   rrtm_tlay(0,k-1))                           &
1245                                   / ( rrtm_play(0,k) - rrtm_play(0,k-1) )     &
1246                                   * ( rrtm_plev(0,k) - rrtm_play(0,k-1) )
1247             ENDDO
1248
1249
1250!
1251!--          Calculate liquid water path and cloud fraction for each column.
1252!--          Note that LWP is required in g/m² instead of kg/kg m.
1253             rrtm_cldfr  = 0.0_wp
1254             rrtm_reliq  = 0.0_wp
1255             rrtm_cliqwp = 0.0_wp
1256             rrtm_icld   = 0
1257
1258             DO k = nzb+1, nzt+1
1259                rrtm_cliqwp(0,k) =  ql(k,j,i) * 1000.0_wp *                    &
1260                                    (rrtm_plev(0,k) - rrtm_plev(0,k+1))        &
1261                                    * 100.0_wp / g
1262
1263                IF ( rrtm_cliqwp(0,k) > 0.0_wp )  THEN
1264                   rrtm_cldfr(0,k) = 1.0_wp
1265                   IF ( rrtm_icld == 0 )  rrtm_icld = 1
1266
1267!
1268!--                Calculate cloud droplet effective radius
1269                   IF ( cloud_physics )  THEN
1270                      rrtm_reliq(0,k) = 1.0E6_wp * ( 3.0_wp * ql(k,j,i)        &
1271                                        * rho_surface                          &
1272                                        / ( 4.0_wp * pi * nc_const * rho_l )   &
1273                                        )**0.33333333333333_wp                 &
1274                                        * EXP( LOG( sigma_gc )**2 )
1275
1276                   ELSEIF ( cloud_droplets )  THEN
1277                      number_of_particles = prt_count(k,j,i)
1278
1279                      IF (number_of_particles <= 0)  CYCLE
1280                      particles => grid_particles(k,j,i)%particles(1:number_of_particles)
1281                      s_r2 = 0.0_wp
1282                      s_r3 = 0.0_wp
1283
1284                      DO  n = 1, number_of_particles
1285                         IF ( particles(n)%particle_mask )  THEN
1286                            s_r2 = s_r2 + particles(n)%radius**2 * &
1287                                   particles(n)%weight_factor
1288                            s_r3 = s_r3 + particles(n)%radius**3 * &
1289                                   particles(n)%weight_factor
1290                         ENDIF
1291                      ENDDO
1292
1293                      IF ( s_r2 > 0.0_wp )  rrtm_reliq(0,k) = s_r3 / s_r2
1294
1295                   ENDIF
1296
1297!
1298!--                Limit effective radius
1299                   IF ( rrtm_reliq(0,k) > 0.0_wp )  THEN
1300                      rrtm_reliq(0,k) = MAX(rrtm_reliq(0,k),2.5_wp)
1301                      rrtm_reliq(0,k) = MIN(rrtm_reliq(0,k),60.0_wp)
1302                  ENDIF
1303                ENDIF
1304             ENDDO
1305
1306!
1307!--          Set surface temperature
1308             rrtm_tsfc = pt(nzb,j,i) * (surface_pressure / 1000.0_wp )**0.286_wp
1309
1310             IF ( lw_radiation )  THEN
1311               CALL rrtmg_lw( 1, nzt_rad      , rrtm_icld    , rrtm_idrv      ,&
1312               rrtm_play       , rrtm_plev    , rrtm_tlay    , rrtm_tlev      ,&
1313               rrtm_tsfc       , rrtm_h2ovmr  , rrtm_o3vmr   , rrtm_co2vmr    ,&
1314               rrtm_ch4vmr     , rrtm_n2ovmr  , rrtm_o2vmr   , rrtm_cfc11vmr  ,&
1315               rrtm_cfc12vmr   , rrtm_cfc22vmr, rrtm_ccl4vmr , rrtm_emis      ,&
1316               rrtm_inflglw    , rrtm_iceflglw, rrtm_liqflglw, rrtm_cldfr     ,&
1317               rrtm_lw_taucld  , rrtm_cicewp  , rrtm_cliqwp  , rrtm_reice     ,& 
1318               rrtm_reliq      , rrtm_lw_tauaer,                               &
1319               rrtm_lwuflx     , rrtm_lwdflx  , rrtm_lwhr  ,                   &
1320               rrtm_lwuflxc    , rrtm_lwdflxc , rrtm_lwhrc ,                   &
1321               rrtm_lwuflx_dt  ,  rrtm_lwuflxc_dt )
1322
1323!
1324!--             Save fluxes
1325                DO k = nzb, nzt+1
1326                   rad_lw_in(k,j,i)  = rrtm_lwdflx(0,k)
1327                   rad_lw_out(k,j,i) = rrtm_lwuflx(0,k)
1328                ENDDO
1329
1330!
1331!--             Save heating rates (convert from K/d to K/h)
1332                DO k = nzb+1, nzt+1
1333                   rad_lw_hr(k,j,i)     = rrtm_lwhr(0,k)  * d_hours_day
1334                   rad_lw_cs_hr(k,j,i)  = rrtm_lwhrc(0,k) * d_hours_day
1335                ENDDO
1336
1337!
1338!--             Save change in LW heating rate
1339                rad_lw_out_change_0(j,i) = rrtm_lwuflx_dt(0,nzb)
1340
1341             ENDIF
1342
1343             IF ( sw_radiation .AND. sun_up )  THEN
1344                CALL rrtmg_sw( 1, nzt_rad      , rrtm_icld  , rrtm_iaer       ,&
1345               rrtm_play       , rrtm_plev    , rrtm_tlay  , rrtm_tlev        ,&
1346               rrtm_tsfc       , rrtm_h2ovmr  , rrtm_o3vmr , rrtm_co2vmr      ,&
1347               rrtm_ch4vmr     , rrtm_n2ovmr  , rrtm_o2vmr , rrtm_asdir(:,j,i),&
1348               rrtm_asdif(:,j,i), rrtm_aldir(:,j,i), rrtm_aldif(:,j,i), zenith,&
1349               0.0_wp          , day          , solar_constant,   rrtm_inflgsw,&
1350               rrtm_iceflgsw   , rrtm_liqflgsw, rrtm_cldfr , rrtm_sw_taucld   ,&
1351               rrtm_sw_ssacld  , rrtm_sw_asmcld, rrtm_sw_fsfcld, rrtm_cicewp  ,&
1352               rrtm_cliqwp     , rrtm_reice   , rrtm_reliq , rrtm_sw_tauaer   ,&
1353               rrtm_sw_ssaaer     , rrtm_sw_asmaer  , rrtm_sw_ecaer ,          &
1354               rrtm_swuflx     , rrtm_swdflx  , rrtm_swhr  ,                   &
1355               rrtm_swuflxc    , rrtm_swdflxc , rrtm_swhrc )
1356 
1357!
1358!--             Save fluxes
1359                DO k = nzb, nzt+1
1360                   rad_sw_in(k,j,i)  = rrtm_swdflx(0,k)
1361                   rad_sw_out(k,j,i) = rrtm_swuflx(0,k)
1362                ENDDO
1363
1364!
1365!--             Save heating rates (convert from K/d to K/s)
1366                DO k = nzb+1, nzt+1
1367                   rad_sw_hr(k,j,i)     = rrtm_swhr(0,k)  * d_hours_day
1368                   rad_sw_cs_hr(k,j,i)  = rrtm_swhrc(0,k) * d_hours_day
1369                ENDDO
1370
1371             ENDIF
1372
1373!
1374!--          Calculate surface net radiation
1375             rad_net(j,i) = rad_sw_in(nzb,j,i) - rad_sw_out(nzb,j,i)           &
1376                            + rad_lw_in(nzb,j,i) - rad_lw_out(nzb,j,i)
1377
1378          ENDDO
1379       ENDDO
1380
1381       CALL exchange_horiz( rad_lw_in,  nbgp )
1382       CALL exchange_horiz( rad_lw_out, nbgp )
1383       CALL exchange_horiz( rad_lw_hr,    nbgp )
1384       CALL exchange_horiz( rad_lw_cs_hr, nbgp )
1385
1386       CALL exchange_horiz( rad_sw_in,  nbgp )
1387       CALL exchange_horiz( rad_sw_out, nbgp ) 
1388       CALL exchange_horiz( rad_sw_hr,    nbgp )
1389       CALL exchange_horiz( rad_sw_cs_hr, nbgp )
1390
1391       CALL exchange_horiz_2d( rad_net, nbgp )
1392       CALL exchange_horiz_2d( rad_lw_out_change_0, nbgp )
1393#endif
1394
1395    END SUBROUTINE radiation_rrtmg
1396
1397
1398!------------------------------------------------------------------------------!
1399! Description:
1400! ------------
1401!> Calculate the cosine of the zenith angle (variable is called zenith)
1402!------------------------------------------------------------------------------!
1403    SUBROUTINE calc_zenith
1404
1405       IMPLICIT NONE
1406
1407       REAL(wp) ::  declination,  & !< solar declination angle
1408                    hour_angle      !< solar hour angle
1409!
1410!--    Calculate current day and time based on the initial values and simulation
1411!--    time
1412       day = day_init + INT(FLOOR( (time_utc_init + time_since_reference_point)    &
1413                               / 86400.0_wp ), KIND=iwp)
1414       time_utc = MOD((time_utc_init + time_since_reference_point), 86400.0_wp)
1415
1416
1417!
1418!--    Calculate solar declination and hour angle   
1419       declination = ASIN( decl_1 * SIN(decl_2 * REAL(day, KIND=wp) - decl_3) )
1420       hour_angle  = 2.0_wp * pi * (time_utc / 86400.0_wp) + lon - pi
1421
1422!
1423!--    Calculate zenith angle
1424       zenith(0) = SIN(lat) * SIN(declination) + COS(lat) * COS(declination)      &
1425                                            * COS(hour_angle)
1426       zenith(0) = MAX(0.0_wp,zenith(0))
1427
1428!
1429!--    Check if the sun is up (otheriwse shortwave calculations can be skipped)
1430       IF ( zenith(0) > 0.0_wp )  THEN
1431          sun_up = .TRUE.
1432       ELSE
1433          sun_up = .FALSE.
1434       END IF
1435
1436    END SUBROUTINE calc_zenith
1437
1438#if defined ( __rrtmg ) && defined ( __netcdf )
1439!------------------------------------------------------------------------------!
1440! Description:
1441! ------------
1442!> Calculates surface albedo components based on Briegleb (1992) and
1443!> Briegleb et al. (1986)
1444!------------------------------------------------------------------------------!
1445    SUBROUTINE calc_albedo
1446
1447        IMPLICIT NONE
1448
1449        IF ( sun_up )  THEN
1450!
1451!--        Ocean
1452           IF ( albedo_type == 1 )  THEN
1453              rrtm_aldir(0,:,:) = 0.026_wp / ( zenith(0)**1.7_wp + 0.065_wp )  &
1454                                  + 0.15_wp * ( zenith(0) - 0.1_wp )           &
1455                                            * ( zenith(0) - 0.5_wp )           &
1456                                            * ( zenith(0) - 1.0_wp )
1457              rrtm_asdir(0,:,:) = rrtm_aldir(0,:,:)
1458!
1459!--        Snow
1460           ELSEIF ( albedo_type == 16 )  THEN
1461              IF ( zenith(0) < 0.5_wp )  THEN
1462                 rrtm_aldir(0,:,:) = 0.5_wp * (1.0_wp - aldif)                 &
1463                                     * ( 3.0_wp / (1.0_wp + 4.0_wp             &
1464                                     * zenith(0))) - 1.0_wp
1465                 rrtm_asdir(0,:,:) = 0.5_wp * (1.0_wp - asdif)                 &
1466                                     * ( 3.0_wp / (1.0_wp + 4.0_wp             &
1467                                     * zenith(0))) - 1.0_wp
1468
1469                 rrtm_aldir(0,:,:) = MIN(0.98_wp, rrtm_aldir(0,:,:))
1470                 rrtm_asdir(0,:,:) = MIN(0.98_wp, rrtm_asdir(0,:,:))
1471              ELSE
1472                 rrtm_aldir(0,:,:) = aldif
1473                 rrtm_asdir(0,:,:) = asdif
1474              ENDIF
1475!
1476!--        Sea ice
1477           ELSEIF ( albedo_type == 15 )  THEN
1478                 rrtm_aldir(0,:,:) = aldif
1479                 rrtm_asdir(0,:,:) = asdif
1480
1481!
1482!--        Asphalt
1483           ELSEIF ( albedo_type == 17 )  THEN
1484                 rrtm_aldir(0,:,:) = aldif
1485                 rrtm_asdir(0,:,:) = asdif
1486!
1487!--        Land surfaces
1488           ELSE
1489              SELECT CASE ( albedo_type )
1490
1491!
1492!--              Surface types with strong zenith dependence
1493                 CASE ( 1, 2, 3, 4, 11, 12, 13 )
1494                    rrtm_aldir(0,:,:) = aldif * 1.4_wp /                       &
1495                                        (1.0_wp + 0.8_wp * zenith(0))
1496                    rrtm_asdir(0,:,:) = asdif * 1.4_wp /                       &
1497                                        (1.0_wp + 0.8_wp * zenith(0))
1498!
1499!--              Surface types with weak zenith dependence
1500                 CASE ( 5, 6, 7, 8, 9, 10, 14 )
1501                    rrtm_aldir(0,:,:) = aldif * 1.1_wp /                       &
1502                                        (1.0_wp + 0.2_wp * zenith(0))
1503                    rrtm_asdir(0,:,:) = asdif * 1.1_wp /                       &
1504                                        (1.0_wp + 0.2_wp * zenith(0))
1505
1506                 CASE DEFAULT
1507
1508              END SELECT
1509           ENDIF
1510!
1511!--        Diffusive albedo is taken from Table 2
1512           rrtm_aldif(0,:,:) = aldif
1513           rrtm_asdif(0,:,:) = asdif
1514
1515        ELSE
1516
1517           rrtm_aldir(0,:,:) = 0.0_wp
1518           rrtm_asdir(0,:,:) = 0.0_wp
1519           rrtm_aldif(0,:,:) = 0.0_wp
1520           rrtm_asdif(0,:,:) = 0.0_wp
1521        ENDIF
1522    END SUBROUTINE calc_albedo
1523
1524!------------------------------------------------------------------------------!
1525! Description:
1526! ------------
1527!> Read sounding data (pressure and temperature) from RADIATION_DATA.
1528!------------------------------------------------------------------------------!
1529    SUBROUTINE read_sounding_data
1530
1531       IMPLICIT NONE
1532
1533       INTEGER(iwp) :: id,           & !< NetCDF id of input file
1534                       id_dim_zrad,  & !< pressure level id in the NetCDF file
1535                       id_var,       & !< NetCDF variable id
1536                       k,            & !< loop index
1537                       nz_snd,       & !< number of vertical levels in the sounding data
1538                       nz_snd_start, & !< start vertical index for sounding data to be used
1539                       nz_snd_end      !< end vertical index for souding data to be used
1540
1541       REAL(wp) :: t_surface           !< actual surface temperature
1542
1543       REAL(wp), DIMENSION(:), ALLOCATABLE ::  hyp_snd_tmp, & !< temporary hydrostatic pressure profile (sounding)
1544                                               t_snd_tmp      !< temporary temperature profile (sounding)
1545
1546!
1547!--    In case of updates, deallocate arrays first (sufficient to check one
1548!--    array as the others are automatically allocated). This is required
1549!--    because nzt_rad might change during the update
1550       IF ( ALLOCATED ( hyp_snd ) )  THEN
1551          DEALLOCATE( hyp_snd )
1552          DEALLOCATE( t_snd )
1553          DEALLOCATE( q_snd  )
1554          DEALLOCATE ( rrtm_play )
1555          DEALLOCATE ( rrtm_plev )
1556          DEALLOCATE ( rrtm_tlay )
1557          DEALLOCATE ( rrtm_tlev )
1558
1559          DEALLOCATE ( rrtm_h2ovmr )
1560          DEALLOCATE ( rrtm_cicewp )
1561          DEALLOCATE ( rrtm_cldfr )
1562          DEALLOCATE ( rrtm_cliqwp )
1563          DEALLOCATE ( rrtm_reice )
1564          DEALLOCATE ( rrtm_reliq )
1565          DEALLOCATE ( rrtm_lw_taucld )
1566          DEALLOCATE ( rrtm_lw_tauaer )
1567
1568          DEALLOCATE ( rrtm_lwdflx  )
1569          DEALLOCATE ( rrtm_lwdflxc )
1570          DEALLOCATE ( rrtm_lwuflx  )
1571          DEALLOCATE ( rrtm_lwuflxc )
1572          DEALLOCATE ( rrtm_lwuflx_dt )
1573          DEALLOCATE ( rrtm_lwuflxc_dt )
1574          DEALLOCATE ( rrtm_lwhr  )
1575          DEALLOCATE ( rrtm_lwhrc )
1576
1577          DEALLOCATE ( rrtm_sw_taucld )
1578          DEALLOCATE ( rrtm_sw_ssacld )
1579          DEALLOCATE ( rrtm_sw_asmcld )
1580          DEALLOCATE ( rrtm_sw_fsfcld )
1581          DEALLOCATE ( rrtm_sw_tauaer )
1582          DEALLOCATE ( rrtm_sw_ssaaer )
1583          DEALLOCATE ( rrtm_sw_asmaer ) 
1584          DEALLOCATE ( rrtm_sw_ecaer )   
1585 
1586          DEALLOCATE ( rrtm_swdflx  )
1587          DEALLOCATE ( rrtm_swdflxc )
1588          DEALLOCATE ( rrtm_swuflx  )
1589          DEALLOCATE ( rrtm_swuflxc )
1590          DEALLOCATE ( rrtm_swhr  )
1591          DEALLOCATE ( rrtm_swhrc )
1592
1593       ENDIF
1594
1595!
1596!--    Open file for reading
1597       nc_stat = NF90_OPEN( rrtm_input_file, NF90_NOWRITE, id )
1598       CALL netcdf_handle_error_rad( 'read_sounding_data', 549 )
1599
1600!
1601!--    Inquire dimension of z axis and save in nz_snd
1602       nc_stat = NF90_INQ_DIMID( id, "Pressure", id_dim_zrad )
1603       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim_zrad, len = nz_snd )
1604       CALL netcdf_handle_error_rad( 'read_sounding_data', 551 )
1605
1606!
1607! !--    Allocate temporary array for storing pressure data
1608       ALLOCATE( hyp_snd_tmp(1:nz_snd) )
1609       hyp_snd_tmp = 0.0_wp
1610
1611
1612!--    Read pressure from file
1613       nc_stat = NF90_INQ_VARID( id, "Pressure", id_var )
1614       nc_stat = NF90_GET_VAR( id, id_var, hyp_snd_tmp(:), start = (/1/),      &
1615                               count = (/nz_snd/) )
1616       CALL netcdf_handle_error_rad( 'read_sounding_data', 552 )
1617
1618!
1619!--    Allocate temporary array for storing temperature data
1620       ALLOCATE( t_snd_tmp(1:nz_snd) )
1621       t_snd_tmp = 0.0_wp
1622
1623!
1624!--    Read temperature from file
1625       nc_stat = NF90_INQ_VARID( id, "ReferenceTemperature", id_var )
1626       nc_stat = NF90_GET_VAR( id, id_var, t_snd_tmp(:), start = (/1/),        &
1627                               count = (/nz_snd/) )
1628       CALL netcdf_handle_error_rad( 'read_sounding_data', 553 )
1629
1630!
1631!--    Calculate start of sounding data
1632       nz_snd_start = nz_snd + 1
1633       nz_snd_end   = nz_snd + 1
1634
1635!
1636!--    Start filling vertical dimension at 10hPa above the model domain (hyp is
1637!--    in Pa, hyp_snd in hPa).
1638       DO  k = 1, nz_snd
1639          IF ( hyp_snd_tmp(k) < ( hyp(nzt+1) - 1000.0_wp) * 0.01_wp )  THEN
1640             nz_snd_start = k
1641             EXIT
1642          END IF
1643       END DO
1644
1645       IF ( nz_snd_start <= nz_snd )  THEN
1646          nz_snd_end = nz_snd
1647       END IF
1648
1649
1650!
1651!--    Calculate of total grid points for RRTMG calculations
1652       nzt_rad = nzt + nz_snd_end - nz_snd_start + 1
1653
1654!
1655!--    Save data above LES domain in hyp_snd, t_snd and q_snd
1656!--    Note: q_snd_tmp is not calculated at the moment (dry residual atmosphere)
1657       ALLOCATE( hyp_snd(nzb+1:nzt_rad) )
1658       ALLOCATE( t_snd(nzb+1:nzt_rad)   )
1659       ALLOCATE( q_snd(nzb+1:nzt_rad)   )
1660       hyp_snd = 0.0_wp
1661       t_snd = 0.0_wp
1662       q_snd = 0.0_wp
1663
1664       hyp_snd(nzt+2:nzt_rad) = hyp_snd_tmp(nz_snd_start+1:nz_snd_end)
1665       t_snd(nzt+2:nzt_rad)   = t_snd_tmp(nz_snd_start+1:nz_snd_end)
1666
1667       nc_stat = NF90_CLOSE( id )
1668
1669!
1670!--    Calculate pressure levels on zu and zw grid. Sounding data is added at
1671!--    top of the LES domain. This routine does not consider horizontal or
1672!--    vertical variability of pressure and temperature
1673       ALLOCATE ( rrtm_play(0:0,nzb+1:nzt_rad+1)   )
1674       ALLOCATE ( rrtm_plev(0:0,nzb+1:nzt_rad+2)   )
1675
1676       t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp
1677       DO k = nzb+1, nzt+1
1678          rrtm_play(0,k) = hyp(k) * 0.01_wp
1679          rrtm_plev(0,k) = surface_pressure * ( (t_surface - g/cp * zw(k-1)) / &
1680                         t_surface )**(1.0_wp/0.286_wp)
1681       ENDDO
1682
1683       DO k = nzt+2, nzt_rad
1684          rrtm_play(0,k) = hyp_snd(k)
1685          rrtm_plev(0,k) = 0.5_wp * ( rrtm_play(0,k) + rrtm_play(0,k-1) )
1686       ENDDO
1687       rrtm_plev(0,nzt_rad+1) = MAX( 0.5 * hyp_snd(nzt_rad),                   &
1688                                   1.5 * hyp_snd(nzt_rad)                      &
1689                                 - 0.5 * hyp_snd(nzt_rad-1) )
1690       rrtm_plev(0,nzt_rad+2)  = MIN( 1.0E-4_wp,                               &
1691                                      0.25_wp * rrtm_plev(0,nzt_rad+1) )
1692
1693       rrtm_play(0,nzt_rad+1) = 0.5 * rrtm_plev(0,nzt_rad+1)
1694
1695!
1696!--    Calculate temperature/humidity levels at top of the LES domain.
1697!--    Currently, the temperature is taken from sounding data (might lead to a
1698!--    temperature jump at interface. To do: Humidity is currently not
1699!--    calculated above the LES domain.
1700       ALLOCATE ( rrtm_tlay(0:0,nzb+1:nzt_rad+1)   )
1701       ALLOCATE ( rrtm_tlev(0:0,nzb+1:nzt_rad+2)   )
1702       ALLOCATE ( rrtm_h2ovmr(0:0,nzb+1:nzt_rad+1) )
1703
1704       DO k = nzt+8, nzt_rad
1705          rrtm_tlay(0,k)   = t_snd(k)
1706          rrtm_h2ovmr(0,k) = q_snd(k)
1707       ENDDO
1708       rrtm_tlay(0,nzt_rad+1) = 2.0_wp * rrtm_tlay(0,nzt_rad)                 &
1709                                - rrtm_tlay(0,nzt_rad-1)
1710       DO k = nzt+9, nzt_rad+1
1711          rrtm_tlev(0,k)   = rrtm_tlay(0,k-1) + (rrtm_tlay(0,k)                &
1712                             - rrtm_tlay(0,k-1))                               &
1713                             / ( rrtm_play(0,k) - rrtm_play(0,k-1) )           &
1714                             * ( rrtm_plev(0,k) - rrtm_play(0,k-1) )
1715       ENDDO
1716       rrtm_h2ovmr(0,nzt_rad+1) = rrtm_h2ovmr(0,nzt_rad)
1717
1718       rrtm_tlev(0,nzt_rad+2)   = 2.0_wp * rrtm_tlay(0,nzt_rad+1)              &
1719                                  - rrtm_tlev(0,nzt_rad)
1720!
1721!--    Allocate remaining RRTMG arrays
1722       ALLOCATE ( rrtm_cicewp(0:0,nzb+1:nzt_rad+1) )
1723       ALLOCATE ( rrtm_cldfr(0:0,nzb+1:nzt_rad+1) )
1724       ALLOCATE ( rrtm_cliqwp(0:0,nzb+1:nzt_rad+1) )
1725       ALLOCATE ( rrtm_reice(0:0,nzb+1:nzt_rad+1) )
1726       ALLOCATE ( rrtm_reliq(0:0,nzb+1:nzt_rad+1) )
1727       ALLOCATE ( rrtm_lw_taucld(1:nbndlw+1,0:0,nzb+1:nzt_rad+1) )
1728       ALLOCATE ( rrtm_lw_tauaer(0:0,nzb+1:nzt_rad+1,1:nbndlw+1) )
1729       ALLOCATE ( rrtm_sw_taucld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
1730       ALLOCATE ( rrtm_sw_ssacld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
1731       ALLOCATE ( rrtm_sw_asmcld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
1732       ALLOCATE ( rrtm_sw_fsfcld(1:nbndsw+1,0:0,nzb+1:nzt_rad+1) )
1733       ALLOCATE ( rrtm_sw_tauaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) )
1734       ALLOCATE ( rrtm_sw_ssaaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) )
1735       ALLOCATE ( rrtm_sw_asmaer(0:0,nzb+1:nzt_rad+1,1:nbndsw+1) ) 
1736       ALLOCATE ( rrtm_sw_ecaer(0:0,nzb+1:nzt_rad+1,1:naerec+1) )   
1737
1738!
1739!--    The ice phase is currently not considered in PALM
1740       rrtm_cicewp = 0.0_wp
1741       rrtm_reice  = 0.0_wp
1742
1743!
1744!--    Set other parameters (move to NAMELIST parameters in the future)
1745       rrtm_lw_tauaer = 0.0_wp
1746       rrtm_lw_taucld = 0.0_wp
1747       rrtm_sw_taucld = 0.0_wp
1748       rrtm_sw_ssacld = 0.0_wp
1749       rrtm_sw_asmcld = 0.0_wp
1750       rrtm_sw_fsfcld = 0.0_wp
1751       rrtm_sw_tauaer = 0.0_wp
1752       rrtm_sw_ssaaer = 0.0_wp
1753       rrtm_sw_asmaer = 0.0_wp
1754       rrtm_sw_ecaer  = 0.0_wp
1755
1756
1757       ALLOCATE ( rrtm_swdflx(0:0,nzb:nzt_rad+1)  )
1758       ALLOCATE ( rrtm_swuflx(0:0,nzb:nzt_rad+1)  )
1759       ALLOCATE ( rrtm_swhr(0:0,nzb+1:nzt_rad+1)  )
1760       ALLOCATE ( rrtm_swuflxc(0:0,nzb:nzt_rad+1) )
1761       ALLOCATE ( rrtm_swdflxc(0:0,nzb:nzt_rad+1) )
1762       ALLOCATE ( rrtm_swhrc(0:0,nzb+1:nzt_rad+1) )
1763
1764       rrtm_swdflx  = 0.0_wp
1765       rrtm_swuflx  = 0.0_wp
1766       rrtm_swhr    = 0.0_wp 
1767       rrtm_swuflxc = 0.0_wp
1768       rrtm_swdflxc = 0.0_wp
1769       rrtm_swhrc   = 0.0_wp
1770
1771       ALLOCATE ( rrtm_lwdflx(0:0,nzb:nzt_rad+1)  )
1772       ALLOCATE ( rrtm_lwuflx(0:0,nzb:nzt_rad+1)  )
1773       ALLOCATE ( rrtm_lwhr(0:0,nzb+1:nzt_rad+1)  )
1774       ALLOCATE ( rrtm_lwuflxc(0:0,nzb:nzt_rad+1) )
1775       ALLOCATE ( rrtm_lwdflxc(0:0,nzb:nzt_rad+1) )
1776       ALLOCATE ( rrtm_lwhrc(0:0,nzb+1:nzt_rad+1) )
1777
1778       rrtm_lwdflx  = 0.0_wp
1779       rrtm_lwuflx  = 0.0_wp
1780       rrtm_lwhr    = 0.0_wp 
1781       rrtm_lwuflxc = 0.0_wp
1782       rrtm_lwdflxc = 0.0_wp
1783       rrtm_lwhrc   = 0.0_wp
1784
1785       ALLOCATE ( rrtm_lwuflx_dt(0:0,nzb:nzt_rad+1) )
1786       ALLOCATE ( rrtm_lwuflxc_dt(0:0,nzb:nzt_rad+1) )
1787
1788       rrtm_lwuflx_dt = 0.0_wp
1789       rrtm_lwuflxc_dt = 0.0_wp
1790
1791    END SUBROUTINE read_sounding_data
1792
1793
1794!------------------------------------------------------------------------------!
1795! Description:
1796! ------------
1797!> Read trace gas data from file
1798!------------------------------------------------------------------------------!
1799    SUBROUTINE read_trace_gas_data
1800
1801       USE rrsw_ncpar
1802
1803       IMPLICIT NONE
1804
1805       INTEGER(iwp), PARAMETER :: num_trace_gases = 9 !< number of trace gases (absorbers)
1806
1807       CHARACTER(LEN=5), DIMENSION(num_trace_gases), PARAMETER ::              & !< trace gas names
1808           trace_names = (/'O3   ', 'CO2  ', 'CH4  ', 'N2O  ', 'O2   ',        &
1809                           'CFC11', 'CFC12', 'CFC22', 'CCL4 '/)
1810
1811       INTEGER(iwp) :: id,     & !< NetCDF id
1812                       k,      & !< loop index
1813                       m,      & !< loop index
1814                       n,      & !< loop index
1815                       nabs,   & !< number of absorbers
1816                       np,     & !< number of pressure levels
1817                       id_abs, & !< NetCDF id of the respective absorber
1818                       id_dim, & !< NetCDF id of asborber's dimension
1819                       id_var    !< NetCDf id ot the absorber
1820
1821       REAL(wp) :: p_mls_l, p_mls_u, p_wgt_l, p_wgt_u, p_mls_m
1822
1823
1824       REAL(wp), DIMENSION(:), ALLOCATABLE   ::  p_mls,         & !< pressure levels for the absorbers
1825                                                 rrtm_play_tmp, & !< temporary array for pressure zu-levels
1826                                                 rrtm_plev_tmp, & !< temporary array for pressure zw-levels
1827                                                 trace_path_tmp   !< temporary array for storing trace gas path data
1828
1829       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  trace_mls,      & !< array for storing the absorber amounts
1830                                                 trace_mls_path, & !< array for storing trace gas path data
1831                                                 trace_mls_tmp     !< temporary array for storing trace gas data
1832
1833
1834!
1835!--    In case of updates, deallocate arrays first (sufficient to check one
1836!--    array as the others are automatically allocated)
1837       IF ( ALLOCATED ( rrtm_o3vmr ) )  THEN
1838          DEALLOCATE ( rrtm_o3vmr  )
1839          DEALLOCATE ( rrtm_co2vmr )
1840          DEALLOCATE ( rrtm_ch4vmr )
1841          DEALLOCATE ( rrtm_n2ovmr )
1842          DEALLOCATE ( rrtm_o2vmr  )
1843          DEALLOCATE ( rrtm_cfc11vmr )
1844          DEALLOCATE ( rrtm_cfc12vmr )
1845          DEALLOCATE ( rrtm_cfc22vmr )
1846          DEALLOCATE ( rrtm_ccl4vmr  )
1847       ENDIF
1848
1849!
1850!--    Allocate trace gas profiles
1851       ALLOCATE ( rrtm_o3vmr(0:0,1:nzt_rad+1)  )
1852       ALLOCATE ( rrtm_co2vmr(0:0,1:nzt_rad+1) )
1853       ALLOCATE ( rrtm_ch4vmr(0:0,1:nzt_rad+1) )
1854       ALLOCATE ( rrtm_n2ovmr(0:0,1:nzt_rad+1) )
1855       ALLOCATE ( rrtm_o2vmr(0:0,1:nzt_rad+1)  )
1856       ALLOCATE ( rrtm_cfc11vmr(0:0,1:nzt_rad+1) )
1857       ALLOCATE ( rrtm_cfc12vmr(0:0,1:nzt_rad+1) )
1858       ALLOCATE ( rrtm_cfc22vmr(0:0,1:nzt_rad+1) )
1859       ALLOCATE ( rrtm_ccl4vmr(0:0,1:nzt_rad+1)  )
1860
1861!
1862!--    Open file for reading
1863       nc_stat = NF90_OPEN( rrtm_input_file, NF90_NOWRITE, id )
1864       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 549 )
1865!
1866!--    Inquire dimension ids and dimensions
1867       nc_stat = NF90_INQ_DIMID( id, "Pressure", id_dim )
1868       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
1869       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, len = np) 
1870       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
1871
1872       nc_stat = NF90_INQ_DIMID( id, "Absorber", id_dim )
1873       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
1874       nc_stat = NF90_INQUIRE_DIMENSION( id, id_dim, len = nabs ) 
1875       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
1876   
1877
1878!
1879!--    Allocate pressure, and trace gas arrays     
1880       ALLOCATE( p_mls(1:np) )
1881       ALLOCATE( trace_mls(1:num_trace_gases,1:np) ) 
1882       ALLOCATE( trace_mls_tmp(1:nabs,1:np) ) 
1883
1884
1885       nc_stat = NF90_INQ_VARID( id, "Pressure", id_var )
1886       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
1887       nc_stat = NF90_GET_VAR( id, id_var, p_mls )
1888       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
1889
1890       nc_stat = NF90_INQ_VARID( id, "AbsorberAmountMLS", id_var )
1891       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
1892       nc_stat = NF90_GET_VAR( id, id_var, trace_mls_tmp )
1893       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 550 )
1894
1895
1896!
1897!--    Write absorber amounts (mls) to trace_mls
1898       DO n = 1, num_trace_gases
1899          CALL getAbsorberIndex( TRIM( trace_names(n) ), id_abs )
1900
1901          trace_mls(n,1:np) = trace_mls_tmp(id_abs,1:np)
1902
1903!
1904!--       Replace missing values by zero
1905          WHERE ( trace_mls(n,:) > 2.0_wp ) 
1906             trace_mls(n,:) = 0.0_wp
1907          END WHERE
1908       END DO
1909
1910       DEALLOCATE ( trace_mls_tmp )
1911
1912       nc_stat = NF90_CLOSE( id )
1913       CALL netcdf_handle_error_rad( 'read_trace_gas_data', 551 )
1914
1915!
1916!--    Add extra pressure level for calculations of the trace gas paths
1917       ALLOCATE ( rrtm_play_tmp(1:nzt_rad+1) )
1918       ALLOCATE ( rrtm_plev_tmp(1:nzt_rad+2) )
1919
1920       rrtm_play_tmp(1:nzt_rad)   = rrtm_play(0,1:nzt_rad) 
1921       rrtm_plev_tmp(1:nzt_rad+1) = rrtm_plev(0,1:nzt_rad+1)
1922       rrtm_play_tmp(nzt_rad+1)   = rrtm_plev(0,nzt_rad+1) * 0.5_wp
1923       rrtm_plev_tmp(nzt_rad+2)   = MIN( 1.0E-4_wp, 0.25_wp                    &
1924                                         * rrtm_plev(0,nzt_rad+1) )
1925 
1926!
1927!--    Calculate trace gas path (zero at surface) with interpolation to the
1928!--    sounding levels
1929       ALLOCATE ( trace_mls_path(1:nzt_rad+2,1:num_trace_gases) )
1930
1931       trace_mls_path(nzb+1,:) = 0.0_wp
1932       
1933       DO k = nzb+2, nzt_rad+2
1934          DO m = 1, num_trace_gases
1935             trace_mls_path(k,m) = trace_mls_path(k-1,m)
1936
1937!
1938!--          When the pressure level is higher than the trace gas pressure
1939!--          level, assume that
1940             IF ( rrtm_plev_tmp(k-1) > p_mls(1) )  THEN             
1941               
1942                trace_mls_path(k,m) = trace_mls_path(k,m) + trace_mls(m,1)     &
1943                                      * ( rrtm_plev_tmp(k-1)                   &
1944                                          - MAX( p_mls(1), rrtm_plev_tmp(k) )  &
1945                                        ) / g
1946             ENDIF
1947
1948!
1949!--          Integrate for each sounding level from the contributing p_mls
1950!--          levels
1951             DO n = 2, np
1952!
1953!--             Limit p_mls so that it is within the model level
1954                p_mls_u = MIN( rrtm_plev_tmp(k-1),                             &
1955                          MAX( rrtm_plev_tmp(k), p_mls(n) ) )
1956                p_mls_l = MIN( rrtm_plev_tmp(k-1),                             &
1957                          MAX( rrtm_plev_tmp(k), p_mls(n-1) ) )
1958
1959                IF ( p_mls_l > p_mls_u )  THEN
1960
1961!
1962!--                Calculate weights for interpolation
1963                   p_mls_m = 0.5_wp * (p_mls_l + p_mls_u)
1964                   p_wgt_u = (p_mls(n-1) - p_mls_m) / (p_mls(n-1) - p_mls(n))
1965                   p_wgt_l = (p_mls_m - p_mls(n))   / (p_mls(n-1) - p_mls(n))
1966
1967!
1968!--                Add level to trace gas path
1969                   trace_mls_path(k,m) = trace_mls_path(k,m)                   &
1970                                         +  ( p_wgt_u * trace_mls(m,n)         &
1971                                            + p_wgt_l * trace_mls(m,n-1) )     &
1972                                         * (p_mls_l - p_mls_u) / g
1973                ENDIF
1974             ENDDO
1975
1976             IF ( rrtm_plev_tmp(k) < p_mls(np) )  THEN
1977                trace_mls_path(k,m) = trace_mls_path(k,m) + trace_mls(m,np)    &
1978                                      * ( MIN( rrtm_plev_tmp(k-1), p_mls(np) ) &
1979                                          - rrtm_plev_tmp(k)                   &
1980                                        ) / g
1981             ENDIF 
1982          ENDDO
1983       ENDDO
1984
1985
1986!
1987!--    Prepare trace gas path profiles
1988       ALLOCATE ( trace_path_tmp(1:nzt_rad+1) )
1989
1990       DO m = 1, num_trace_gases
1991
1992          trace_path_tmp(1:nzt_rad+1) = ( trace_mls_path(2:nzt_rad+2,m)        &
1993                                       - trace_mls_path(1:nzt_rad+1,m) ) * g   &
1994                                       / ( rrtm_plev_tmp(1:nzt_rad+1)          &
1995                                       - rrtm_plev_tmp(2:nzt_rad+2) )
1996
1997!
1998!--       Save trace gas paths to the respective arrays
1999          SELECT CASE ( TRIM( trace_names(m) ) )
2000
2001             CASE ( 'O3' )
2002
2003                rrtm_o3vmr(0,:) = trace_path_tmp(:)
2004
2005             CASE ( 'CO2' )
2006
2007                rrtm_co2vmr(0,:) = trace_path_tmp(:)
2008
2009             CASE ( 'CH4' )
2010
2011                rrtm_ch4vmr(0,:) = trace_path_tmp(:)
2012
2013             CASE ( 'N2O' )
2014
2015                rrtm_n2ovmr(0,:) = trace_path_tmp(:)
2016
2017             CASE ( 'O2' )
2018
2019                rrtm_o2vmr(0,:) = trace_path_tmp(:)
2020
2021             CASE ( 'CFC11' )
2022
2023                rrtm_cfc11vmr(0,:) = trace_path_tmp(:)
2024
2025             CASE ( 'CFC12' )
2026
2027                rrtm_cfc12vmr(0,:) = trace_path_tmp(:)
2028
2029             CASE ( 'CFC22' )
2030
2031                rrtm_cfc22vmr(0,:) = trace_path_tmp(:)
2032
2033             CASE ( 'CCL4' )
2034
2035                rrtm_ccl4vmr(0,:) = trace_path_tmp(:)
2036
2037             CASE DEFAULT
2038
2039          END SELECT
2040
2041       ENDDO
2042
2043       DEALLOCATE ( trace_path_tmp )
2044       DEALLOCATE ( trace_mls_path )
2045       DEALLOCATE ( rrtm_play_tmp )
2046       DEALLOCATE ( rrtm_plev_tmp )
2047       DEALLOCATE ( trace_mls )
2048       DEALLOCATE ( p_mls )
2049
2050    END SUBROUTINE read_trace_gas_data
2051
2052
2053    SUBROUTINE netcdf_handle_error_rad( routine_name, errno )
2054
2055       USE control_parameters,                                                 &
2056           ONLY:  message_string
2057
2058       USE NETCDF
2059
2060       USE pegrid
2061
2062       IMPLICIT NONE
2063
2064       CHARACTER(LEN=6) ::  message_identifier
2065       CHARACTER(LEN=*) ::  routine_name
2066
2067       INTEGER(iwp) ::  errno
2068
2069       IF ( nc_stat /= NF90_NOERR )  THEN
2070
2071          WRITE( message_identifier, '(''NC'',I4.4)' )  errno
2072          message_string = TRIM( NF90_STRERROR( nc_stat ) )
2073
2074          CALL message( routine_name, message_identifier, 2, 2, 0, 6, 1 )
2075
2076       ENDIF
2077
2078    END SUBROUTINE netcdf_handle_error_rad
2079#endif
2080
2081
2082!------------------------------------------------------------------------------!
2083! Description:
2084! ------------
2085!> Calculate temperature tendency due to radiative cooling/heating.
2086!> Cache-optimized version.
2087!------------------------------------------------------------------------------!
2088    SUBROUTINE radiation_tendency_ij ( i, j, tend )
2089
2090       USE cloud_parameters,                                                   &
2091           ONLY:  pt_d_t
2092
2093       IMPLICIT NONE
2094
2095       INTEGER(iwp) :: i, j, k !< loop indices
2096
2097       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend !< pt tendency term
2098
2099#if defined ( __rrtmg )
2100!
2101!--    Calculate tendency based on heating rate
2102       DO k = nzb+1, nzt+1
2103          tend(k,j,i) = tend(k,j,i) + (rad_lw_hr(k,j,i) + rad_sw_hr(k,j,i))    &
2104                                      * pt_d_t(k) * d_seconds_hour
2105       ENDDO
2106
2107#endif
2108
2109    END SUBROUTINE radiation_tendency_ij
2110
2111
2112!------------------------------------------------------------------------------!
2113! Description:
2114! ------------
2115!> Calculate temperature tendency due to radiative cooling/heating.
2116!> Vector-optimized version
2117!------------------------------------------------------------------------------!
2118    SUBROUTINE radiation_tendency ( tend )
2119
2120       USE cloud_parameters,                                                   &
2121           ONLY:  pt_d_t
2122
2123       USE indices,                                                            &
2124           ONLY:  nxl, nxr, nyn, nys
2125
2126       IMPLICIT NONE
2127
2128       INTEGER(iwp) :: i, j, k !< loop indices
2129
2130       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: tend !< pt tendency term
2131
2132#if defined ( __rrtmg )
2133!
2134!--    Calculate tendency based on heating rate
2135       DO  i = nxl, nxr
2136          DO  j = nys, nyn
2137             DO k = nzb+1, nzt+1
2138                tend(k,j,i) = tend(k,j,i) + ( rad_lw_hr(k,j,i)                 &
2139                                            +  rad_sw_hr(k,j,i) ) * pt_d_t(k)  &
2140                                            * d_seconds_hour
2141             ENDDO
2142         ENDDO
2143       ENDDO
2144#endif
2145
2146    END SUBROUTINE radiation_tendency
2147
2148 END MODULE radiation_model_mod
Note: See TracBrowser for help on using the repository browser.