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

Last change on this file since 3274 was 3274, checked in by knoop, 3 years ago

Modularization of all bulk cloud physics code components

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