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

Last change on this file since 1321 was 1321, checked in by raasch, 7 years ago

last commit documented

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