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

Last change on this file since 1850 was 1850, checked in by maronga, 6 years ago

added _mod string to several filenames to meet the naming convection for modules

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