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

Last change on this file since 1682 was 1682, checked in by knoop, 8 years ago

Code annotations made doxygen readable

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