source: palm/trunk/SOURCE/diffusion_e.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: 11.7 KB
Line 
1 MODULE diffusion_e_mod
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Log: diffusion_e.f90,v $
11! Revision 1.18  2006/08/04 14:29:43  raasch
12! dissipation is stored in extra array diss if needed later on for calculating
13! the sgs particle velocities
14!
15! Revision 1.17  2006/02/23 10:31:46  raasch
16! nzb_2d replaced by nzb_s_inner
17!
18! Revision 1.16  2004/01/30 10:18:18  raasch
19! Scalar lower k index nzb replaced by 2d-array nzb_2d
20!
21! Revision 1.15  2003/03/14 13:39:33  raasch
22! Loop optimization for diffusion_e, l and ll are now automatic arrays
23!
24! Revision 1.14  2003/03/12 16:25:03  raasch
25! Full code replaced in the call for all gridpoints instead of calling the
26! _ij version (required by NEC, because otherwise no vectorization)
27!
28! Revision 1.13  2002/12/19 14:25:28  raasch
29! Correction of mixing length term (l(k)/ll(k)). The condition kh=3*km in
30! the unstable case is now also exactly met in the wall adjustment region.
31! Factor 0.7 in wall adjustment part replaced by variable
32! wall_adjustment_factor, which is set to ... in modules.f90.
33!
34! Revision 1.12  2002/06/11 12:51:59  raasch
35! Former subroutine changed to a module which allows to be called for all grid
36! points of a single vertical column with index i,j or for all grid points by
37! using function overloading.
38! 1D-array l is allocated only once in the first call.
39!
40! Revision 1.11  2001/08/21 08:24:34  raasch
41! Wall adjustment of mixing length to 0.7 z can be switched off
42!
43! Revision 1.10  2001/03/30 07:06:30  raasch
44! Near surface mixing length is limited to 0.7*zu,
45! e**1.5 replaced by e*SQRT(e) (more than 10% total increase in performance
46! of this routine),
47! Translation of remaining German identifiers (variables, subroutines, etc.)
48!
49! Revision 1.9  2001/01/22 06:05:28  raasch
50! Module test_variables removed
51!
52! Revision 1.8  2001/01/02 17:27:00  raasch
53! -dpt_dz_d, dpt_dz_u
54!
55! Revision 1.7  2000/07/03 12:56:34  raasch
56! array l changed from dummy argument to local allocatable array,
57! dummy arguments, whose corresponding actual arguments are pointers,
58! are now also defined as pointers
59! all comments translated into English
60!
61! Revision 1.6  2000/04/18 08:10:12  schroeter
62! Revision 1.4 wieder rueckgaengig gemacht, das Stabilitaets-
63! kriterium basiert nun wieder auf zentralen Differenzen
64!
65! Revision 1.5  2000/04/13 14:33:08  schroeter
66! je nach Initialisierungsmodus (trocken/feucht) fliesst in die
67! Berechnung des Mischungsweges pt oder vpt ein, wird durch
68! entsprechende Variablenuebergabe geregelt
69!
70! Revision 1.4  99/02/17  09:15:52  09:15:52  raasch (Siegfried Raasch)
71! Dissipation jetzt gemaess dem originalen Deardorff-Ansatz
72! Kriterium fuer reduzierten Mischungsweg im stabil geschichteten Fall enger
73! gefasst (Schichtung muss sowohl oberhalb als auch unterhalb des betrachteten
74! Gitterpunkts stabil sein)
75!
76! Revision 1.3  1998/07/06 12:10:56  raasch
77! + USE test_variables
78!
79! Revision 1.2  1998/03/11 11:48:59  raasch
80! Anpassung des Mischungsweges an den Prandtlschen Mischungsweg moeglich
81!
82! Revision 1.1  1997/09/19 07:40:24  raasch
83! Initial revision
84!
85!
86! Description:
87! ------------
88! Diffusion- and dissipation terms for the TKE
89!------------------------------------------------------------------------------!
90
91    PRIVATE
92    PUBLIC diffusion_e
93   
94
95    INTERFACE diffusion_e
96       MODULE PROCEDURE diffusion_e
97       MODULE PROCEDURE diffusion_e_ij
98    END INTERFACE diffusion_e
99 
100 CONTAINS
101
102
103!------------------------------------------------------------------------------!
104! Call for all grid points
105!------------------------------------------------------------------------------!
106    SUBROUTINE diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, theta, &
107                            rif, tend, zu )
108
109       USE control_parameters
110       USE grid_variables
111       USE indices
112       USE particle_attributes
113
114       IMPLICIT NONE
115
116       INTEGER ::  i, j, k
117       REAL            ::  dpt_dz, l_stable, phi_m
118       REAL            ::  ddzu(1:nzt+1), dd2zu(1:nzt), ddzw(1:nzt), &
119                           l_grid(1:nzt), zu(0:nzt+1)
120       REAL, DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) :: diss, tend
121       REAL, DIMENSION(:,:), POINTER   ::  rif
122       REAL, DIMENSION(:,:,:), POINTER ::  e, km, theta
123       REAL, DIMENSION(nzb+1:nzt-1,nys:nyn) ::  dissipation, l, ll
124 
125
126       DO  i = nxl, nxr
127          DO  j = nys, nyn
128!
129!--          First, calculate phi-function for eventually adjusting the &
130!--          mixing length to the prandtl mixing length
131             IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
132                IF ( rif(j,i) >= 0.0 )  THEN
133                   phi_m = 1.0 + 5.0 * rif(j,i)
134                ELSE
135                   phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
136                ENDIF
137             ENDIF
138
139             DO  k = nzb_s_inner(j,i)+1, nzt-1
140!
141!--             Calculate the mixing length (for dissipation)
142                dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k)
143                IF ( dpt_dz > 0.0 ) THEN
144                   l_stable = 0.76 * SQRT( e(k,j,i) ) / &
145                                     SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5
146                ELSE
147                   l_stable = l_grid(k)
148                ENDIF
149!
150!--             Adjustment of the mixing length
151                IF ( wall_adjustment )  THEN
152                   l(k,j)  = MIN( wall_adjustment_factor * zu(k), l_grid(k), &
153                                  l_stable )
154                   ll(k,j) = MIN( wall_adjustment_factor * zu(k), l_grid(k) )
155                ELSE
156                   l(k,j)  = MIN( l_grid(k), l_stable )
157                   ll(k,j) = l_grid(k)
158                ENDIF
159                IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
160                   l(k,j)  = MIN( l(k,j),  kappa * zu(k) / phi_m )
161                   ll(k,j) = MIN( ll(k,j), kappa * zu(k) / phi_m )
162                ENDIF
163
164             ENDDO
165          ENDDO
166!
167!--       Calculate the tendency terms
168          DO  j = nys, nyn
169             DO  k = nzb_s_inner(j,i)+1, nzt-1
170
171                 dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
172                                    e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
173
174                 tend(k,j,i) = tend(k,j,i)                                     &
175                                        + (                                    &
176                          ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
177                        - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
178                                          ) * ddx2                             &
179                                        + (                                    &
180                          ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
181                        - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
182                                          ) * ddy2                             &
183                                        + (                                    &
184               ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
185             - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
186                                          ) * ddzw(k)                          &
187                             - dissipation(k,j)
188
189             ENDDO
190          ENDDO
191
192!
193!--       Store dissipation if needed for calculating the sgs particle
194!--       velocities
195          IF ( use_sgs_for_particles )  THEN
196             DO  j = nys, nyn
197                DO  k = nzb_s_inner(j,i)+1, nzt-1
198                   diss(k,j,i) = dissipation(k,j)
199                ENDDO
200             ENDDO
201          ENDIF
202
203       ENDDO
204
205!
206!--    Boundary condition for dissipation
207       IF ( use_sgs_for_particles )  THEN
208          DO  i = nxl, nxr
209             DO  j = nys, nyn
210                diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
211             ENDDO
212          ENDDO
213       ENDIF
214
215    END SUBROUTINE diffusion_e
216
217
218!------------------------------------------------------------------------------!
219! Call for grid point i,j
220!------------------------------------------------------------------------------!
221    SUBROUTINE diffusion_e_ij( i, j, ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
222                               theta, rif, tend, zu )
223
224       USE control_parameters
225       USE grid_variables
226       USE indices
227       USE particle_attributes
228
229       IMPLICIT NONE
230
231       INTEGER         ::  i, j, k
232       REAL            ::  dpt_dz, l_stable, phi_m
233       REAL            ::  ddzu(1:nzt+1), dd2zu(1:nzt), ddzw(1:nzt), &
234                           l_grid(1:nzt), zu(0:nzt+1)
235       REAL, DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) ::  diss, tend
236       REAL, DIMENSION(:,:), POINTER   ::  rif
237       REAL, DIMENSION(:,:,:), POINTER ::  e, km, theta
238       REAL, DIMENSION(nzb+1:nzt-1)    ::  dissipation, l, ll
239
240
241!
242!--    First, calculate phi-function for eventually adjusting the mixing length
243!--    to the prandtl mixing length
244       IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
245          IF ( rif(j,i) >= 0.0 )  THEN
246             phi_m = 1.0 + 5.0 * rif(j,i)
247          ELSE
248             phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
249          ENDIF
250       ENDIF
251
252!
253!--    Calculate the mixing length (for dissipation)
254       DO  k = nzb_s_inner(j,i)+1, nzt-1
255          dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k)
256          IF ( dpt_dz > 0.0 ) THEN
257             l_stable = 0.76 * SQRT( e(k,j,i) ) / &
258                               SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5
259          ELSE
260             l_stable = l_grid(k)
261          ENDIF
262!
263!--       Adjustment of the mixing length
264          IF ( wall_adjustment )  THEN
265             l(k)  = MIN( wall_adjustment_factor * zu(k), l_grid(k), l_stable )
266             ll(k) = MIN( wall_adjustment_factor * zu(k), l_grid(k) )
267          ELSE
268             l(k)  = MIN( l_grid(k), l_stable )
269             ll(k) = l_grid(k)
270          ENDIF
271          IF ( adjust_mixing_length  .AND.  prandtl_layer )  THEN
272             l(k)  = MIN( l(k),  kappa * zu(k) / phi_m )
273             ll(k) = MIN( ll(k), kappa * zu(k) / phi_m )
274          ENDIF
275
276!
277!--       Calculate the tendency term
278          dissipation(k) = ( 0.19 + 0.74 * l(k) / ll(k) ) * e(k,j,i) * &
279                           SQRT( e(k,j,i) ) / l(k)
280
281          tend(k,j,i) = tend(k,j,i)                                           &
282                                       + (                                    &
283                         ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )  &
284                       - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )  &
285                                         ) * ddx2                             &
286                                       + (                                    &
287                         ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )  &
288                       - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )  &
289                                         ) * ddy2                             &
290                                       + (                                    &
291              ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
292            - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)   &
293                                         ) * ddzw(k)                          &
294                                       - dissipation(k)
295
296       ENDDO
297
298!
299!--    Store dissipation if needed for calculating the sgs particle velocities
300       IF ( use_sgs_for_particles )  THEN
301          DO  k = nzb_s_inner(j,i)+1, nzt-1
302             diss(k,j,i) = dissipation(k)
303          ENDDO
304!
305!--       Boundary condition for dissipation
306          diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
307       ENDIF
308
309    END SUBROUTINE diffusion_e_ij
310
311 END MODULE diffusion_e_mod
Note: See TracBrowser for help on using the repository browser.