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

Last change on this file since 2707 was 2696, checked in by kanani, 7 years ago

Merge of branch palm4u into trunk

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