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

Last change on this file since 1 was 1, checked in by raasch, 15 years ago

Initial repository layout and content

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