source: palm/trunk/SOURCE/calc_radiation.f90 @ 1873

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

revised renaming of modules

  • Property svn:keywords set to Id
File size: 13.1 KB
Line 
1!> @file calc_radiation.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! Module renamed (removed _mod)
22!
23!
24! Former revisions:
25! -----------------
26! $Id: calc_radiation.f90 1873 2016-04-18 14:50:06Z maronga $
27!
28! 1850 2016-04-08 13:29:27Z maronga
29! Module renamed
30!
31!
32! 1682 2015-10-07 23:56:08Z knoop
33! Code annotations made doxygen readable
34!
35! 1353 2014-04-08 15:21:23Z heinze
36! REAL constants provided with KIND-attribute
37!
38! 1322 2014-03-20 16:38:49Z raasch
39! exponent 4.0 changed to integer
40!
41! 1320 2014-03-20 08:40:49Z raasch
42! ONLY-attribute added to USE-statements,
43! kind-parameters added to all INTEGER and REAL declaration statements,
44! kinds are defined in new module kinds,
45! revision history before 2012 removed,
46! comment fields (!:) to be used for variable explanations added to
47! all variable declaration statements
48!
49! 1036 2012-10-22 13:43:42Z raasch
50! code put under GPL (PALM 3.9)
51!
52! Revision 1.1  2000/04/13 14:42:45  schroeter
53! Initial revision
54!
55!
56! Description:
57! -------------
58!> Calculation of the vertical divergences of the long-wave radiation-fluxes
59!> based on the parameterization of the cloud effective emissivity
60!------------------------------------------------------------------------------!
61 MODULE calc_radiation_mod
62 
63    USE kinds
64   
65    PRIVATE
66    PUBLIC calc_radiation
67   
68    LOGICAL, SAVE ::  first_call = .TRUE. !<
69    REAL(wp), SAVE ::  sigma = 5.67E-08_wp   !<
70
71    REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE ::  lwp_ground         !<
72    REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE ::  lwp_top            !<
73    REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE ::  blackbody_emission !<
74
75    INTERFACE calc_radiation
76       MODULE PROCEDURE calc_radiation
77       MODULE PROCEDURE calc_radiation_ij
78    END INTERFACE calc_radiation
79 
80 CONTAINS
81
82
83!------------------------------------------------------------------------------!
84! Description:
85! ------------
86!> Call for all grid points
87!------------------------------------------------------------------------------!
88    SUBROUTINE calc_radiation
89
90       USE arrays_3d,                                                          &
91           ONLY:  dzw, pt, ql, tend
92
93       USE cloud_parameters,                                                   &
94           ONLY:  cp, l_d_cp, pt_d_t, t_d_pt
95
96       USE control_parameters,                                                 &
97           ONLY:  rho_surface
98
99       USE indices,                                                            &
100           ONLY:  nxl, nxr, nyn, nys, nzb, nzb_s_inner, nzt
101
102       USE kinds
103
104       USE pegrid
105
106
107       IMPLICIT NONE
108
109       INTEGER(iwp) ::  i      !<
110       INTEGER(iwp) ::  j      !<
111       INTEGER(iwp) ::  k      !<
112       INTEGER(iwp) ::  k_help !<
113 
114       REAL(wp) :: df_p                      !<
115       REAL(wp) :: df_m                      !<
116       REAL(wp) :: effective_emission_up_m   !<
117       REAL(wp) :: effective_emission_up_p   !<
118       REAL(wp) :: effective_emission_down_m !<
119       REAL(wp) :: effective_emission_down_p !<
120       REAL(wp) :: f_up_m                    !<
121       REAL(wp) :: f_up_p                    !<
122       REAL(wp) :: f_down_m                  !<
123       REAL(wp) :: f_down_p                  !<
124       REAL(wp) :: impinging_flux_at_top     !<
125       REAL(wp) :: temperature               !<
126
127
128!
129!--    On first call, allocate temporary arrays
130       IF ( first_call )  THEN
131          ALLOCATE( blackbody_emission(nzb:nzt+1), lwp_ground(nzb:nzt+1),      &
132                    lwp_top(nzb:nzt+1) )
133          first_call = .FALSE.
134       ENDIF
135
136
137       DO  i = nxl, nxr
138          DO  j = nys, nyn
139!
140!--          Compute the liquid water path (LWP) and blackbody_emission
141!--          at all vertical levels
142             lwp_ground(nzb) = 0.0_wp
143             lwp_top(nzt+1)  = rho_surface * ql(nzt+1,j,i) * dzw(nzt+1)
144
145             temperature     = pt(nzb,j,i) * t_d_pt(nzb) + l_d_cp * ql(nzb,j,i)
146             blackbody_emission(nzb) = sigma * temperature**4
147
148             DO  k = nzb_s_inner(j,i)+1, nzt
149
150                k_help = ( nzt+nzb+1 ) - k
151                lwp_ground(k)   = lwp_ground(k-1) + rho_surface * ql(k,j,i) *  &
152                                  dzw(k)
153
154                lwp_top(k_help) = lwp_top(k_help+1) +                          &
155                                  rho_surface * ql(k_help,j,i) * dzw(k_help)
156
157                temperature     = pt(k,j,i) * t_d_pt(k) + l_d_cp * ql(k,j,i)
158                blackbody_emission(k) = sigma * temperature**4
159
160             ENDDO
161
162             lwp_ground(nzt+1) = lwp_ground(nzt) +                             &
163                                 rho_surface * ql(nzt+1,j,i) * dzw(nzt+1)
164             lwp_top(nzb)      = lwp_top(nzb+1)
165
166             temperature       = pt(nzt+1,j,i) * t_d_pt(nzt+1) + l_d_cp *      &
167                                 ql(nzt+1,j,i)
168             blackbody_emission(nzt+1) = sigma * temperature**4
169
170!
171!--          See Chlond '92, this is just a first guess
172             impinging_flux_at_top = blackbody_emission(nzb) - 100.0_wp
173
174             DO  k = nzb_s_inner(j,i)+1, nzt
175!
176!--             Save some computational time, but this may cause load
177!--             imbalances if ql is not distributed uniformly
178                IF ( ql(k,j,i) /= 0.0_wp ) THEN
179!
180!--                Compute effective emissivities
181                   effective_emission_up_p   = 1.0_wp -                        &
182                                               EXP( -130.0_wp * lwp_ground(k+1) )
183                   effective_emission_up_m   = 1.0_wp -                        &
184                                               EXP( -130.0_wp * lwp_ground(k-1) )
185                   effective_emission_down_p = 1.0_wp -                        &
186                                               EXP( -158.0_wp * lwp_top(k+1) )
187                   effective_emission_down_m = 1.0_wp -                        &
188                                               EXP( -158.0_wp * lwp_top(k-1) ) 
189
190!
191!--                Compute vertical long wave radiation fluxes
192                   f_up_p = blackbody_emission(nzb) +                          &
193                            effective_emission_up_p *                          &
194                           ( blackbody_emission(k) - blackbody_emission(nzb) )
195
196                   f_up_m = blackbody_emission(nzb) +                          &
197                            effective_emission_up_m *                          &
198                           ( blackbody_emission(k-1) - blackbody_emission(nzb) )
199
200                   f_down_p = impinging_flux_at_top +                          &
201                              effective_emission_down_p *                      &
202                             ( blackbody_emission(k) - impinging_flux_at_top )
203
204                   f_down_m = impinging_flux_at_top +                          &
205                              effective_emission_down_m *                      &
206                             ( blackbody_emission(k-1) - impinging_flux_at_top )
207
208!
209!--                Divergence of vertical long wave radiation fluxes
210                   df_p = f_up_p - f_down_p
211                   df_m = f_up_m - f_down_m
212
213!
214!--                Compute tendency term         
215                   tend(k,j,i) = tend(k,j,i) -                                 &
216                                ( pt_d_t(k) / ( rho_surface * cp ) *           &
217                                  ( df_p - df_m ) / dzw(k) )
218
219                ENDIF
220
221             ENDDO
222          ENDDO
223       ENDDO
224
225    END SUBROUTINE calc_radiation
226
227
228!------------------------------------------------------------------------------!
229! Description:
230! ------------
231!> Call for grid point i,j
232!------------------------------------------------------------------------------!
233    SUBROUTINE calc_radiation_ij( i, j )
234
235       USE arrays_3d,                                                          &
236           ONLY:  dzw, pt, ql, tend
237
238       USE cloud_parameters,                                                   &
239           ONLY:  cp, l_d_cp, pt_d_t, t_d_pt
240
241       USE control_parameters,                                                 &
242           ONLY:  rho_surface
243
244       USE indices,                                                            &
245           ONLY:  nzb, nzb_s_inner, nzt
246
247       USE kinds
248
249       USE pegrid
250
251   
252       IMPLICIT NONE
253
254       INTEGER(iwp) ::  i      !<
255       INTEGER(iwp) ::  j      !<
256       INTEGER(iwp) ::  k      !<
257       INTEGER(iwp) ::  k_help !<
258
259       REAL(wp) :: df_p                      !<
260       REAL(wp) :: df_m                      !<
261       REAL(wp) :: effective_emission_up_m   !<
262       REAL(wp) :: effective_emission_up_p   !<
263       REAL(wp) :: effective_emission_down_m !<
264       REAL(wp) :: effective_emission_down_p !<
265       REAL(wp) :: f_up_m                    !<
266       REAL(wp) :: f_up_p                    !<
267       REAL(wp) :: f_down_m                  !<
268       REAL(wp) :: f_down_p                  !<
269       REAL(wp) :: impinging_flux_at_top     !<
270       REAL(wp) :: temperature               !<
271
272       
273!
274!--    On first call, allocate temporary arrays
275       IF ( first_call )  THEN
276          ALLOCATE( blackbody_emission(nzb:nzt+1), lwp_ground(nzb:nzt+1),      &
277                    lwp_top(nzb:nzt+1) )
278          first_call = .FALSE.
279       ENDIF
280
281!
282!--    Compute the liquid water path (LWP) and blackbody_emission
283!--    at all vertical levels
284       lwp_ground(nzb) = 0.0_wp
285       lwp_top(nzt+1)  = rho_surface * ql(nzt+1,j,i) * dzw(nzt+1)
286
287       temperature     = pt(nzb,j,i) * t_d_pt(nzb) + l_d_cp * ql(nzb,j,i)
288       blackbody_emission(nzb) = sigma * temperature**4
289
290       DO  k = nzb_s_inner(j,i)+1, nzt
291          k_help = ( nzt+nzb+1 ) - k
292          lwp_ground(k)   = lwp_ground(k-1) + rho_surface * ql(k,j,i) * dzw(k)
293
294          lwp_top(k_help) = lwp_top(k_help+1) +                                &
295                            rho_surface * ql(k_help,j,i) * dzw(k_help)
296
297          temperature     = pt(k,j,i) * t_d_pt(k) + l_d_cp * ql(k,j,i)
298          blackbody_emission(k) = sigma * temperature**4
299
300       ENDDO
301       lwp_ground(nzt+1) = lwp_ground(nzt) +                                   &
302                           rho_surface * ql(nzt+1,j,i) * dzw(nzt+1)
303       lwp_top(nzb)      = lwp_top(nzb+1)
304
305       temperature       = pt(nzt+1,j,i) * t_d_pt(nzt+1) + l_d_cp *            &
306                           ql(nzt+1,j,i)
307       blackbody_emission(nzt+1) = sigma * temperature**4
308
309!
310!--    See Chlond '92, this is just a first guess
311       impinging_flux_at_top = blackbody_emission(nzb) - 100.0_wp
312
313       DO  k = nzb_s_inner(j,i)+1, nzt
314!
315!--       Store some computational time,
316!--       this may cause load imbalances if ql is not distributed uniformly
317          IF ( ql(k,j,i) /= 0.0_wp ) THEN
318!
319!--          Compute effective emissivities
320             effective_emission_up_p   = 1.0_wp -                              &
321                                         EXP( -130.0_wp * lwp_ground(k+1) )
322             effective_emission_up_m   = 1.0_wp -                              &
323                                         EXP( -130.0_wp * lwp_ground(k-1) )
324             effective_emission_down_p = 1.0_wp -                              &
325                                         EXP( -158.0_wp * lwp_top(k+1) )
326             effective_emission_down_m = 1.0_wp -                              &
327                                         EXP( -158.0_wp * lwp_top(k-1) ) 
328             
329!
330!--          Compute vertical long wave radiation fluxes
331             f_up_p = blackbody_emission(nzb) + effective_emission_up_p *      &
332                     ( blackbody_emission(k) - blackbody_emission(nzb) )
333
334             f_up_m = blackbody_emission(nzb) + effective_emission_up_m *      &
335                     ( blackbody_emission(k-1) - blackbody_emission(nzb) )
336
337             f_down_p = impinging_flux_at_top + effective_emission_down_p *    &
338                       ( blackbody_emission(k) - impinging_flux_at_top )
339
340             f_down_m = impinging_flux_at_top + effective_emission_down_m *    &
341                       ( blackbody_emission(k-1) - impinging_flux_at_top )
342
343!
344!-           Divergence of vertical long wave radiation fluxes
345             df_p = f_up_p - f_down_p
346             df_m = f_up_m - f_down_m
347
348!
349!--          Compute tendency term         
350             tend(k,j,i) = tend(k,j,i) - ( pt_d_t(k) / ( rho_surface * cp ) *  &
351                                         ( df_p - df_m ) / dzw(k) )
352
353          ENDIF
354
355       ENDDO
356
357    END SUBROUTINE calc_radiation_ij
358
359 END MODULE calc_radiation_mod
Note: See TracBrowser for help on using the repository browser.