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

Last change on this file since 1036 was 1036, checked in by raasch, 9 years ago

code has been put under the GNU General Public License (v3)

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