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

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

Forced header and separation lines into 80 columns

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