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

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

last commit documented / copyright update

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