source: palm/trunk/SOURCE/radiation_model.f90 @ 1710

Last change on this file since 1710 was 1710, checked in by maronga, 6 years ago

last commit documented

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