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

Last change on this file since 1325 was 1323, checked in by raasch, 10 years ago

last commit documented

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