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

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

last commit documented

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