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

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

ONLY-attribute added to USE-statements,
kind-parameters added to all INTEGER and REAL declaration statements,
kinds are defined in new module kinds,
old module precision_kind is removed,
revision history before 2012 removed,
comment fields (!:) to be used for variable explanations added to all variable declaration statements

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