source: palm/tags/release-4.0/SOURCE/calc_radiation.f90 @ 3999

Last change on this file since 3999 was 1354, checked in by heinze, 10 years ago

last commit documented

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