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

Last change on this file since 1709 was 1709, checked in by maronga, 7 years ago

several bugfixes related to the new surface layer routine and land-surface-radiation interaction

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