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

Last change on this file since 1310 was 1310, checked in by raasch, 7 years ago

update of GPL copyright

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