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

Last change on this file since 1937 was 1874, checked in by maronga, 8 years ago

last commit documented

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