Home

Context Navigation

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

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |