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

Last change on this file since 3039 was 2718, checked in by maronga, 7 years ago

deleting of deprecated files; headers updated where needed

  • Property svn:keywords set to Id
File size: 14.1 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 2718 2018-01-02 08:49:38Z schwenkel $
27! Corrected "Former revisions" section
28!
29! 2696 2017-12-14 17:12:51Z kanani
30! Change in file header (GPL part)
31!
32! 2233 2017-05-30 18:08:54Z suehring
33!
34! 2232 2017-05-30 17:47:52Z suehring
35! Adjustments to new topography concept
36!
37! 2000 2016-08-20 18:09:15Z knoop
38! Forced header and separation lines into 80 columns
39!
40! 1873 2016-04-18 14:50:06Z maronga
41! Module renamed (removed _mod)
42!
43!
44! 1850 2016-04-08 13:29:27Z maronga
45! Module renamed
46!
47!
48! 1682 2015-10-07 23:56:08Z knoop
49! Code annotations made doxygen readable
50!
51! 1353 2014-04-08 15:21:23Z heinze
52! REAL constants provided with KIND-attribute
53!
54! 1322 2014-03-20 16:38:49Z raasch
55! exponent 4.0 changed to integer
56!
57! 1320 2014-03-20 08:40:49Z raasch
58! ONLY-attribute added to USE-statements,
59! kind-parameters added to all INTEGER and REAL declaration statements,
60! kinds are defined in new module kinds,
61! revision history before 2012 removed,
62! comment fields (!:) to be used for variable explanations added to
63! all variable declaration statements
64!
65! 1036 2012-10-22 13:43:42Z raasch
66! code put under GPL (PALM 3.9)
67!
68! Revision 1.1  2000/04/13 14:42:45  schroeter
69! Initial revision
70!
71!
72! Description:
73! -------------
74!> Calculation of the vertical divergences of the long-wave radiation-fluxes
75!> based on the parameterization of the cloud effective emissivity
76!------------------------------------------------------------------------------!
77 MODULE calc_radiation_mod
78 
79    USE kinds
80   
81    PRIVATE
82    PUBLIC calc_radiation
83   
84    LOGICAL, SAVE ::  first_call = .TRUE. !<
85    REAL(wp), SAVE ::  sigma = 5.67E-08_wp   !<
86
87    REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE ::  lwp_ground         !<
88    REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE ::  lwp_top            !<
89    REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE ::  blackbody_emission !<
90
91    INTERFACE calc_radiation
92       MODULE PROCEDURE calc_radiation
93       MODULE PROCEDURE calc_radiation_ij
94    END INTERFACE calc_radiation
95 
96 CONTAINS
97
98
99!------------------------------------------------------------------------------!
100! Description:
101! ------------
102!> Call for all grid points
103!------------------------------------------------------------------------------!
104    SUBROUTINE calc_radiation
105
106       USE arrays_3d,                                                          &
107           ONLY:  dzw, pt, ql, tend
108
109       USE cloud_parameters,                                                   &
110           ONLY:  cp, l_d_cp, pt_d_t, t_d_pt
111
112       USE control_parameters,                                                 &
113           ONLY:  rho_surface
114
115       USE indices,                                                            &
116           ONLY:  nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0
117
118       USE kinds
119
120       USE pegrid
121
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) * t_d_pt(nzb) + l_d_cp * ql(nzb,j,i)
162             blackbody_emission(nzb) = sigma * 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) * t_d_pt(k) + l_d_cp * ql(k,j,i)
174                blackbody_emission(k) = sigma * 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) * t_d_pt(nzt+1) + l_d_cp *      &
185                                 ql(nzt+1,j,i)
186             blackbody_emission(nzt+1) = sigma * 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                                ( pt_d_t(k) / ( rho_surface * cp ) *           &
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       USE arrays_3d,                                                          &
256           ONLY:  dzw, pt, ql, tend
257
258       USE cloud_parameters,                                                   &
259           ONLY:  cp, l_d_cp, pt_d_t, t_d_pt
260
261       USE control_parameters,                                                 &
262           ONLY:  rho_surface
263
264       USE indices,                                                            &
265           ONLY:  nzb, nzt, wall_flags_0
266
267       USE kinds
268
269       USE pegrid
270
271   
272       IMPLICIT NONE
273
274       INTEGER(iwp) ::  i      !<
275       INTEGER(iwp) ::  j      !<
276       INTEGER(iwp) ::  k      !<
277       INTEGER(iwp) ::  k_help !<
278
279       REAL(wp) :: df_p                      !<
280       REAL(wp) :: df_m                      !<
281       REAL(wp) :: effective_emission_up_m   !<
282       REAL(wp) :: effective_emission_up_p   !<
283       REAL(wp) :: effective_emission_down_m !<
284       REAL(wp) :: effective_emission_down_p !<
285       REAL(wp) :: f_up_m                    !<
286       REAL(wp) :: f_up_p                    !<
287       REAL(wp) :: f_down_m                  !<
288       REAL(wp) :: f_down_p                  !<
289       REAL(wp) :: impinging_flux_at_top     !<
290       REAL(wp) :: temperature               !<
291
292       
293!
294!--    On first call, allocate temporary arrays
295       IF ( first_call )  THEN
296          ALLOCATE( blackbody_emission(nzb:nzt+1), lwp_ground(nzb:nzt+1),      &
297                    lwp_top(nzb:nzt+1) )
298          first_call = .FALSE.
299       ENDIF
300
301!
302!--    Compute the liquid water path (LWP) and blackbody_emission
303!--    at all vertical levels
304       lwp_ground(nzb) = 0.0_wp
305       lwp_top(nzt+1)  = rho_surface * ql(nzt+1,j,i) * dzw(nzt+1)
306
307       temperature     = pt(nzb,j,i) * t_d_pt(nzb) + l_d_cp * ql(nzb,j,i)
308       blackbody_emission(nzb) = sigma * temperature**4
309
310       DO  k = nzb+1, nzt
311          k_help = ( nzt+nzb+1 ) - k
312          lwp_ground(k)   = lwp_ground(k-1) + rho_surface * ql(k,j,i) * dzw(k)
313
314          lwp_top(k_help) = lwp_top(k_help+1) +                                &
315                            rho_surface * ql(k_help,j,i) * dzw(k_help)
316
317          temperature     = pt(k,j,i) * t_d_pt(k) + l_d_cp * ql(k,j,i)
318          blackbody_emission(k) = sigma * temperature**4                       &
319                                      * MERGE( 1.0_wp, 0.0_wp,                 &
320                                               BTEST( wall_flags_0(k,j,i), 0 ) )
321
322       ENDDO
323       lwp_ground(nzt+1) = lwp_ground(nzt) +                                   &
324                           rho_surface * ql(nzt+1,j,i) * dzw(nzt+1)
325       lwp_top(nzb)      = lwp_top(nzb+1)
326
327       temperature       = pt(nzt+1,j,i) * t_d_pt(nzt+1) + l_d_cp *            &
328                           ql(nzt+1,j,i)
329       blackbody_emission(nzt+1) = sigma * temperature**4
330
331!
332!--    See Chlond '92, this is just a first guess
333       impinging_flux_at_top = blackbody_emission(nzb) - 100.0_wp
334
335       DO  k = nzb+1, nzt
336!
337!--       Store some computational time,
338!--       this may cause load imbalances if ql is not distributed uniformly
339          IF ( ql(k,j,i) /= 0.0_wp ) THEN
340!
341!--          Compute effective emissivities
342             effective_emission_up_p   = 1.0_wp -                              &
343                                         EXP( -130.0_wp * lwp_ground(k+1) )
344             effective_emission_up_m   = 1.0_wp -                              &
345                                         EXP( -130.0_wp * lwp_ground(k-1) )
346             effective_emission_down_p = 1.0_wp -                              &
347                                         EXP( -158.0_wp * lwp_top(k+1) )
348             effective_emission_down_m = 1.0_wp -                              &
349                                         EXP( -158.0_wp * lwp_top(k-1) ) 
350             
351!
352!--          Compute vertical long wave radiation fluxes
353             f_up_p = blackbody_emission(nzb) + effective_emission_up_p *      &
354                     ( blackbody_emission(k) - blackbody_emission(nzb) )
355
356             f_up_m = blackbody_emission(nzb) + effective_emission_up_m *      &
357                     ( blackbody_emission(k-1) - blackbody_emission(nzb) )
358
359             f_down_p = impinging_flux_at_top + effective_emission_down_p *    &
360                       ( blackbody_emission(k) - impinging_flux_at_top )
361
362             f_down_m = impinging_flux_at_top + effective_emission_down_m *    &
363                       ( blackbody_emission(k-1) - impinging_flux_at_top )
364
365!
366!-           Divergence of vertical long wave radiation fluxes
367             df_p = f_up_p - f_down_p
368             df_m = f_up_m - f_down_m
369
370!
371!--          Compute tendency term         
372             tend(k,j,i) = tend(k,j,i) - ( pt_d_t(k) / ( rho_surface * cp ) *  &
373                                         ( df_p - df_m ) / dzw(k) )            &
374                                      * MERGE( 1.0_wp, 0.0_wp,                 &
375                                               BTEST( wall_flags_0(k,j,i), 0 ) )
376
377          ENDIF
378
379       ENDDO
380
381    END SUBROUTINE calc_radiation_ij
382
383 END MODULE calc_radiation_mod
Note: See TracBrowser for help on using the repository browser.