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

Last change on this file since 2232 was 2232, checked in by suehring, 4 years ago

Adjustments according new topography and surface-modelling concept implemented

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