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

Last change on this file since 550 was 484, checked in by raasch, 14 years ago

typo in file headers removed

  • Property svn:keywords set to Id
File size: 9.3 KB
Line 
1 MODULE calc_radiation_mod
2
3!------------------------------------------------------------------------------!
4! Current revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: calc_radiation.f90 484 2010-02-05 07:36:54Z maronga $
11! RCS Log replace by Id keyword, revision history cleaned up
12!
13! Revision 1.6  2004/01/30 10:17:03  raasch
14! Scalar lower k index nzb replaced by 2d-array nzb_2d
15!
16! Revision 1.1  2000/04/13 14:42:45  schroeter
17! Initial revision
18!
19!
20! Description:
21! -------------
22! Calculation of the vertical divergences of the long-wave radiation-fluxes
23! based on the parameterization of the cloud effective emissivity
24!------------------------------------------------------------------------------!
25
26    PRIVATE
27    PUBLIC calc_radiation
28   
29    LOGICAL, SAVE ::  first_call = .TRUE.
30    REAL, SAVE    ::  sigma = 5.67E-08
31
32    REAL, DIMENSION(:), ALLOCATABLE, SAVE ::  lwp_ground, lwp_top, &
33                                              blackbody_emission
34
35    INTERFACE calc_radiation
36       MODULE PROCEDURE calc_radiation
37       MODULE PROCEDURE calc_radiation_ij
38    END INTERFACE calc_radiation
39 
40 CONTAINS
41
42
43!------------------------------------------------------------------------------!
44! Call for all grid points
45!------------------------------------------------------------------------------!
46    SUBROUTINE calc_radiation
47
48       USE arrays_3d
49       USE cloud_parameters
50       USE control_parameters
51       USE indices
52       USE pegrid
53
54       IMPLICIT NONE
55
56       INTEGER ::  i, j, k, k_help
57 
58       REAL :: df_p, df_m , effective_emission_up_m, effective_emission_up_p, &
59               effective_emission_down_m, effective_emission_down_p,          &
60               f_up_m, f_up_p, f_down_m, f_down_p, impinging_flux_at_top,     &
61               temperature
62
63
64!
65!--    On first call, allocate temporary arrays
66       IF ( first_call )  THEN
67          ALLOCATE( blackbody_emission(nzb:nzt+1), lwp_ground(nzb:nzt+1), &
68                    lwp_top(nzb:nzt+1) )
69          first_call = .FALSE.
70       ENDIF
71
72
73       DO  i = nxl, nxr
74          DO  j = nys, nyn
75!
76!--          Compute the liquid water path (LWP) and blackbody_emission
77!--          at all vertical levels
78             lwp_ground(nzb) = 0.0
79             lwp_top(nzt+1)  = rho_surface * ql(nzt+1,j,i) * dzw(nzt+1)
80
81             temperature     = pt(nzb,j,i) * t_d_pt(nzb) + l_d_cp * ql(nzb,j,i)
82             blackbody_emission(nzb) = sigma * temperature**4.0
83
84             DO  k = nzb_2d(j,i)+1, nzt
85
86                k_help = ( nzt+nzb+1 ) - k
87                lwp_ground(k)   = lwp_ground(k-1) + rho_surface * ql(k,j,i) * &
88                                  dzw(k)
89
90                lwp_top(k_help) = lwp_top(k_help+1) + &
91                                  rho_surface * ql(k_help,j,i) * dzw(k_help)
92
93                temperature     = pt(k,j,i) * t_d_pt(k) + l_d_cp * ql(k,j,i)
94                blackbody_emission(k) = sigma * temperature**4.0
95
96             ENDDO
97
98             lwp_ground(nzt+1) = lwp_ground(nzt) + &
99                                 rho_surface * ql(nzt+1,j,i) * dzw(nzt+1)
100             lwp_top(nzb)      = lwp_top(nzb+1)
101
102             temperature       = pt(nzt+1,j,i) * t_d_pt(nzt+1) + l_d_cp * &
103                                 ql(nzt+1,j,i)
104             blackbody_emission(nzt+1) = sigma * temperature**4.0
105
106!
107!--          See Chlond '92, this is just a first guess
108             impinging_flux_at_top = blackbody_emission(nzb) - 100.0
109
110             DO  k = nzb_2d(j,i)+1, nzt
111!
112!--             Save some computational time, but this may cause load
113!--             imbalances if ql is not distributed uniformly
114                IF ( ql(k,j,i) /= 0.0 ) THEN
115!
116!--                Compute effective emissivities
117                   effective_emission_up_p   = 1.0 - &
118                                               EXP( -130.0 * lwp_ground(k+1) )
119                   effective_emission_up_m   = 1.0 - &
120                                               EXP( -130.0 * lwp_ground(k-1) )
121                   effective_emission_down_p = 1.0 - &
122                                               EXP( -158.0 * lwp_top(k+1) )
123                   effective_emission_down_m = 1.0 - &
124                                               EXP( -158.0 * lwp_top(k-1) ) 
125
126!
127!--                Compute vertical long wave radiation fluxes
128                   f_up_p = blackbody_emission(nzb) + &
129                            effective_emission_up_p * &
130                           ( blackbody_emission(k) - blackbody_emission(nzb) )
131
132                   f_up_m = blackbody_emission(nzb) + &
133                            effective_emission_up_m * &
134                           ( blackbody_emission(k-1) - blackbody_emission(nzb) )
135
136                   f_down_p = impinging_flux_at_top + &
137                              effective_emission_down_p * &
138                             ( blackbody_emission(k) - impinging_flux_at_top )
139
140                   f_down_m = impinging_flux_at_top + &
141                              effective_emission_down_m * &
142                             ( blackbody_emission(k-1) - impinging_flux_at_top )
143
144!
145!--                Divergence of vertical long wave radiation fluxes
146                   df_p = f_up_p - f_down_p
147                   df_m = f_up_m - f_down_m
148
149!
150!--                Compute tendency term         
151                   tend(k,j,i) = tend(k,j,i) - &
152                                ( pt_d_t(k) / ( rho_surface * cp ) * &
153                                  ( df_p - df_m ) / dzw(k) )
154
155                ENDIF
156
157             ENDDO
158          ENDDO
159       ENDDO
160
161    END SUBROUTINE calc_radiation
162
163
164!------------------------------------------------------------------------------!
165! Call for grid point i,j
166!------------------------------------------------------------------------------!
167    SUBROUTINE calc_radiation_ij( i, j )
168
169       USE arrays_3d
170       USE cloud_parameters
171       USE control_parameters
172       USE indices
173       USE pegrid
174   
175       IMPLICIT NONE
176
177       INTEGER :: i, j, k, k_help
178
179       REAL :: df_p, df_m , effective_emission_up_m, effective_emission_up_p, &
180               effective_emission_down_m, effective_emission_down_p,          &
181               f_up_m, f_up_p, f_down_m, f_down_p, impinging_flux_at_top,     &
182               temperature
183
184!
185!--    On first call, allocate temporary arrays
186       IF ( first_call )  THEN
187          ALLOCATE( blackbody_emission(nzb:nzt+1), lwp_ground(nzb:nzt+1), &
188                    lwp_top(nzb:nzt+1) )
189          first_call = .FALSE.
190       ENDIF
191
192!
193!--    Compute the liquid water path (LWP) and blackbody_emission
194!--    at all vertical levels
195       lwp_ground(nzb) = 0.0
196       lwp_top(nzt+1)  = rho_surface * ql(nzt+1,j,i) * dzw(nzt+1)
197
198       temperature     = pt(nzb,j,i) * t_d_pt(nzb) + l_d_cp * ql(nzb,j,i)
199       blackbody_emission(nzb) = sigma * temperature**4.0
200
201       DO  k = nzb_2d(j,i)+1, nzt
202          k_help = ( nzt+nzb+1 ) - k
203          lwp_ground(k)   = lwp_ground(k-1) + rho_surface * ql(k,j,i) * dzw(k)
204
205          lwp_top(k_help) = lwp_top(k_help+1) + &
206                            rho_surface * ql(k_help,j,i) * dzw(k_help)
207
208          temperature     = pt(k,j,i) * t_d_pt(k) + l_d_cp * ql(k,j,i)
209          blackbody_emission(k) = sigma * temperature**4.0
210
211       ENDDO
212       lwp_ground(nzt+1) = lwp_ground(nzt) + &
213                           rho_surface * ql(nzt+1,j,i) * dzw(nzt+1)
214       lwp_top(nzb)      = lwp_top(nzb+1)
215
216       temperature       = pt(nzt+1,j,i) * t_d_pt(nzt+1) + l_d_cp * &
217                           ql(nzt+1,j,i)
218       blackbody_emission(nzt+1) = sigma * temperature**4.0
219
220!
221!--    See Chlond '92, this is just a first guess
222       impinging_flux_at_top = blackbody_emission(nzb) - 100.0
223
224       DO  k = nzb_2d(j,i)+1, nzt
225!
226!--       Store some computational time,
227!--       this may cause load imbalances if ql is not distributed uniformly
228          IF ( ql(k,j,i) /= 0.0 ) THEN
229!
230!--          Compute effective emissivities
231             effective_emission_up_p   = 1.0 - &
232                                         EXP( -130.0 * lwp_ground(k+1) )
233             effective_emission_up_m   = 1.0 - &
234                                         EXP( -130.0 * lwp_ground(k-1) )
235             effective_emission_down_p = 1.0 - &
236                                         EXP( -158.0 * lwp_top(k+1) )
237             effective_emission_down_m = 1.0 - &
238                                         EXP( -158.0 * lwp_top(k-1) ) 
239             
240!
241!--          Compute vertical long wave radiation fluxes
242             f_up_p = blackbody_emission(nzb) + effective_emission_up_p * &
243                     ( blackbody_emission(k) - blackbody_emission(nzb) )
244
245             f_up_m = blackbody_emission(nzb) + effective_emission_up_m * &
246                     ( blackbody_emission(k-1) - blackbody_emission(nzb) )
247
248             f_down_p = impinging_flux_at_top + effective_emission_down_p * &
249                       ( blackbody_emission(k) - impinging_flux_at_top )
250
251             f_down_m = impinging_flux_at_top + effective_emission_down_m * &
252                       ( blackbody_emission(k-1) - impinging_flux_at_top )
253
254!
255!-           Divergence of vertical long wave radiation fluxes
256             df_p = f_up_p - f_down_p
257             df_m = f_up_m - f_down_m
258
259!
260!--          Compute tendency term         
261             tend(k,j,i) = tend(k,j,i) - ( pt_d_t(k) / ( rho_surface * cp ) * &
262                                         ( df_p - df_m ) / dzw(k) )
263
264          ENDIF
265
266       ENDDO
267
268    END SUBROUTINE calc_radiation_ij
269
270 END MODULE calc_radiation_mod
Note: See TracBrowser for help on using the repository browser.