source: palm/trunk/SOURCE/calc_radiation_mod.f90 @ 1851

Last change on this file since 1851 was 1851, checked in by maronga, 5 years ago

last commit documented

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