Home

Context Navigation

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

Last change on this file since 1322 was 1322, checked in by raasch, 10 years ago
REAL functions and a lot of REAL constants provided with KIND-attribute, some small bugfixes
Property svn:keywords set to `Id`
File size: 12.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |