source: palm/trunk/SOURCE/production_e.f90 @ 75

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

preliminary update for changes concerning non-cyclic boundary conditions

  • Property svn:keywords set to Id
File size: 37.4 KB
Line 
1 MODULE production_e_mod
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6! Wall functions now include diabatic conditions, call of routine wall_fluxes_e,
7! reference temperature pt_reference can be used in buoyancy term,
8! moisture renamed humidity
9!
10! Former revisions:
11! -----------------
12! $Id: production_e.f90 75 2007-03-22 09:54:05Z raasch $
13!
14! 37 2007-03-01 08:33:54Z raasch
15! Calculation extended for gridpoint nzt, extended for given temperature /
16! humidity fluxes at the top, wall-part is now executed in case that a
17! Prandtl-layer is switched on (instead of surfaces fluxes switched on)
18!
19! RCS Log replace by Id keyword, revision history cleaned up
20!
21! Revision 1.21  2006/04/26 12:45:35  raasch
22! OpenMP parallelization of production_e_init
23!
24! Revision 1.1  1997/09/19 07:45:35  raasch
25! Initial revision
26!
27!
28! Description:
29! ------------
30! Production terms (shear + buoyancy) of the TKE
31! WARNING: the case with prandtl_layer = F and use_surface_fluxes = T is
32!          not considered well!
33!------------------------------------------------------------------------------!
34
35    USE wall_fluxes_mod
36
37    PRIVATE
38    PUBLIC production_e, production_e_init
39
40    LOGICAL, SAVE ::  first_call = .TRUE.
41
42    REAL, DIMENSION(:,:), ALLOCATABLE, SAVE ::  u_0, v_0
43
44    INTERFACE production_e
45       MODULE PROCEDURE production_e
46       MODULE PROCEDURE production_e_ij
47    END INTERFACE production_e
48   
49    INTERFACE production_e_init
50       MODULE PROCEDURE production_e_init
51    END INTERFACE production_e_init
52 
53 CONTAINS
54
55
56!------------------------------------------------------------------------------!
57! Call for all grid points
58!------------------------------------------------------------------------------!
59    SUBROUTINE production_e
60
61       USE arrays_3d
62       USE cloud_parameters
63       USE control_parameters
64       USE grid_variables
65       USE indices
66       USE statistics
67
68       IMPLICIT NONE
69
70       INTEGER ::  i, j, k
71
72       REAL    ::  def, dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz, &
73                   k1, k2, theta, temp
74
75!       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  usvs, vsus, wsus, wsvs
76       REAL, DIMENSION(nzb:nzt+1) ::   usvs, vsus, wsus, wsvs
77
78!
79!--    First calculate horizontal momentum flux u'v', w'v', v'u', w'u' at
80!--    vertical walls, if neccessary
81!--    So far, results are slightly different from the ij-Version.
82!--    Therefore, ij-Version is called further below within the ij-loops.
83!       IF ( topography /= 'flat' )  THEN
84!          CALL wall_fluxes_e( usvs, 1.0, 0.0, 0.0, 0.0, wall_e_y )
85!          CALL wall_fluxes_e( wsvs, 0.0, 0.0, 1.0, 0.0, wall_e_y )
86!          CALL wall_fluxes_e( vsus, 0.0, 1.0, 0.0, 0.0, wall_e_x )
87!          CALL wall_fluxes_e( wsus, 0.0, 0.0, 0.0, 1.0, wall_e_x )
88!       ENDIF
89
90!
91!--    Calculate TKE production by shear
92       DO  i = nxl, nxr
93
94          DO  j = nys, nyn
95             DO  k = nzb_diff_s_outer(j,i), nzt
96
97                dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
98                dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
99                                 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
100                dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
101                                 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
102
103                dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
104                                 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
105                dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
106                dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
107                                 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
108
109                dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
110                                 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
111                dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
112                                 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
113                dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
114
115                def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
116                      dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
117                      dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
118
119                IF ( def < 0.0 )  def = 0.0
120
121                tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
122               
123             ENDDO
124          ENDDO
125
126          IF ( prandtl_layer )  THEN
127
128!
129!--          Position beneath wall
130!--          (2) - Will allways be executed.
131!--          'bottom and wall: use u_0,v_0 and wall functions'
132             DO  j = nys, nyn
133
134                IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) ) &
135                THEN
136
137                   k = nzb_diff_s_inner(j,i) - 1
138                   dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
139                   dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
140                                  u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
141                   dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
142                   dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
143                                  v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
144                   dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
145
146                   IF ( wall_e_y(j,i) /= 0.0 )  THEN
147                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
148                                          usvs, 1.0, 0.0, 0.0, 0.0 )
149                      dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i)
150!                      dudy = wall_e_y(j,i) * usvs(k,j,i) / km(k,j,i)
151                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
152                                          wsvs, 0.0, 0.0, 1.0, 0.0 )
153                      dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i)
154!                      dwdy = wall_e_y(j,i) * wsvs(k,j,i) / km(k,j,i)
155                   ELSE
156                      dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
157                                      u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
158                      dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
159                                      w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
160                   ENDIF
161
162                   IF ( wall_e_x(j,i) /= 0.0 )  THEN
163                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
164                                          vsus, 0.0, 1.0, 0.0, 0.0 )
165                      dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i)
166!                      dvdx = wall_e_x(j,i) * vsus(k,j,i) / km(k,j,i)
167                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
168                                          wsus, 0.0, 0.0, 0.0, 1.0 )
169                      dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i)
170!                      dwdx = wall_e_x(j,i) * wsus(k,j,i) / km(k,j,i)
171                   ELSE
172                      dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
173                                      v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
174                      dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
175                                      w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
176                   ENDIF
177
178                   def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
179                         dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
180                         dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
181
182                   IF ( def < 0.0 )  def = 0.0
183
184                   tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
185
186
187!
188!--                (3) - will be executed only, if there is at least one level
189!--                between (2) and (4), i.e. the topography must have a
190!--                minimum height of 2 dz. Wall fluxes for this case have
191!--                already been calculated for (2).
192!--                'wall only: use wall functions'
193
194                   DO  k = nzb_diff_s_inner(j,i), nzb_diff_s_outer(j,i)-2
195
196                      dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
197                      dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
198                                     u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
199                      dvdy =       ( v(k,j+1,i) - v(k,j,i)     ) * ddy
200                      dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
201                                     v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
202                      dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
203
204                      IF ( wall_e_y(j,i) /= 0.0 )  THEN
205                         dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i)
206!                         dudy = wall_e_y(j,i) * usvs(k,j,i) / km(k,j,i)
207                         dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i)
208!                         dwdy = wall_e_y(j,i) * wsvs(k,j,i) / km(k,j,i)
209                      ELSE
210                         dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
211                                         u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
212                         dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
213                                         w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
214                      ENDIF
215
216                      IF ( wall_e_x(j,i) /= 0.0 )  THEN
217                         dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i)
218!                         dvdx = wall_e_x(j,i) * vsus(k,j,i) / km(k,j,i)
219                         dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i)
220!                         dwdx = wall_e_x(j,i) * wsus(k,j,i) / km(k,j,i)
221                      ELSE
222                         dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
223                                         v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
224                         dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
225                                         w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
226                      ENDIF
227
228                      def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
229                           dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 +  &
230                           dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
231
232                      IF ( def < 0.0 )  def = 0.0
233
234                      tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
235
236                   ENDDO
237
238                ENDIF
239
240             ENDDO
241
242!
243!--          (4) - will allways be executed.
244!--          'special case: free atmosphere' (as for case (0))
245             DO  j = nys, nyn
246
247                IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) ) &
248                THEN
249
250                   k = nzb_diff_s_outer(j,i)-1
251
252                   dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
253                   dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
254                                    u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
255                   dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
256                                    u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
257
258                   dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
259                                    v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
260                   dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
261                   dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
262                                    v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
263
264                   dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
265                                    w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
266                   dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
267                                    w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
268                   dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
269
270                   def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
271                         dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
272                         dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
273
274                   IF ( def < 0.0 )  def = 0.0
275
276                   tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
277
278                ENDIF
279
280             ENDDO
281
282!
283!--          Position without adjacent wall
284!--          (1) - will allways be executed.
285!--          'bottom only: use u_0,v_0'
286             DO  j = nys, nyn
287
288                IF ( ( wall_e_x(j,i) == 0.0 ) .AND. ( wall_e_y(j,i) == 0.0 ) ) &
289                THEN
290
291                   k = nzb_diff_s_inner(j,i)-1
292
293                   dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
294                   dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
295                                    u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
296                   dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
297                                    u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
298
299                   dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
300                                    v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
301                   dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
302                   dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
303                                    v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
304
305                   dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
306                                    w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
307                   dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
308                                    w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
309                   dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
310
311                   def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
312                         dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
313                         dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
314
315                   IF ( def < 0.0 )  def = 0.0
316
317                   tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
318               
319                ENDIF
320
321             ENDDO
322
323          ELSEIF ( use_surface_fluxes )  THEN
324
325             DO  j = nys, nyn
326
327                k = nzb_diff_s_outer(j,i)-1
328
329                dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
330                dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
331                                 u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
332                dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
333                                 u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
334
335                dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
336                                 v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
337                dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
338                dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
339                                 v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
340
341                dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
342                                 w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
343                dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
344                                 w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
345                dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
346
347                def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
348                      dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
349                      dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
350
351                IF ( def < 0.0 )  def = 0.0
352
353                tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
354
355             ENDDO
356
357          ENDIF
358
359!
360!--       Calculate TKE production by buoyancy
361          IF ( .NOT. humidity )  THEN
362
363             IF ( use_pt_reference )  THEN
364
365                DO  j = nys, nyn
366                   DO  k = nzb_diff_s_inner(j,i), nzt_diff
367                      tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g/pt_reference * &
368                                       ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
369                   ENDDO
370
371                   IF ( use_surface_fluxes )  THEN
372                      k = nzb_diff_s_inner(j,i)-1
373                      tend(k,j,i) = tend(k,j,i) + g / pt_reference * shf(j,i)
374                   ENDIF
375
376                   IF ( use_top_fluxes )  THEN
377                      k = nzt
378                      tend(k,j,i) = tend(k,j,i) + g / pt_reference * tswst(j,i)
379                   ENDIF
380                ENDDO
381
382             ELSE
383
384                DO  j = nys, nyn
385                   DO  k = nzb_diff_s_inner(j,i), nzt_diff
386                      tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt(k,j,i) * &
387                                       ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
388                   ENDDO
389
390                   IF ( use_surface_fluxes )  THEN
391                      k = nzb_diff_s_inner(j,i)-1
392                      tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
393                   ENDIF
394
395                   IF ( use_top_fluxes )  THEN
396                      k = nzt
397                      tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
398                   ENDIF
399                ENDDO
400
401             ENDIF
402
403          ELSE
404
405             DO  j = nys, nyn
406
407                DO  k = nzb_diff_s_inner(j,i), nzt_diff
408
409                   IF ( .NOT. cloud_physics )  THEN
410                      k1 = 1.0 + 0.61 * q(k,j,i)
411                      k2 = 0.61 * pt(k,j,i)
412                   ELSE
413                      IF ( ql(k,j,i) == 0.0 )  THEN
414                         k1 = 1.0 + 0.61 * q(k,j,i)
415                         k2 = 0.61 * pt(k,j,i)
416                      ELSE
417                         theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
418                         temp  = theta * t_d_pt(k)
419                         k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
420                                    ( q(k,j,i) - ql(k,j,i) ) *          &
421                              ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
422                              ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
423                              ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
424                         k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
425                      ENDIF
426                   ENDIF
427
428                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * &
429                                      ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +  &
430                                        k2 * ( q(k+1,j,i) - q(k-1,j,i) )    &
431                                      ) * dd2zu(k)
432                ENDDO
433
434             ENDDO
435
436             IF ( use_surface_fluxes )  THEN
437
438                DO  j = nys, nyn
439
440                   k = nzb_diff_s_inner(j,i)-1
441
442                   IF ( .NOT. cloud_physics )  THEN
443                      k1 = 1.0 + 0.61 * q(k,j,i)
444                      k2 = 0.61 * pt(k,j,i)
445                   ELSE
446                      IF ( ql(k,j,i) == 0.0 )  THEN
447                         k1 = 1.0 + 0.61 * q(k,j,i)
448                         k2 = 0.61 * pt(k,j,i)
449                      ELSE
450                         theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
451                         temp  = theta * t_d_pt(k)
452                         k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
453                                    ( q(k,j,i) - ql(k,j,i) ) *          &
454                              ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
455                              ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
456                              ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
457                         k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
458                      ENDIF
459                   ENDIF
460
461                   tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
462                                         ( k1* shf(j,i) + k2 * qsws(j,i) )
463                ENDDO
464
465             ENDIF
466
467             IF ( use_top_fluxes )  THEN
468
469                DO  j = nys, nyn
470
471                   k = nzt
472
473                   IF ( .NOT. cloud_physics )  THEN
474                      k1 = 1.0 + 0.61 * q(k,j,i)
475                      k2 = 0.61 * pt(k,j,i)
476                   ELSE
477                      IF ( ql(k,j,i) == 0.0 )  THEN
478                         k1 = 1.0 + 0.61 * q(k,j,i)
479                         k2 = 0.61 * pt(k,j,i)
480                      ELSE
481                         theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
482                         temp  = theta * t_d_pt(k)
483                         k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
484                                    ( q(k,j,i) - ql(k,j,i) ) *          &
485                              ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
486                              ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
487                              ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
488                         k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
489                      ENDIF
490                   ENDIF
491
492                   tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
493                                         ( k1* tswst(j,i) + k2 * qswst(j,i) )
494                ENDDO
495
496             ENDIF
497
498          ENDIF
499
500       ENDDO
501
502    END SUBROUTINE production_e
503
504
505!------------------------------------------------------------------------------!
506! Call for grid point i,j
507!------------------------------------------------------------------------------!
508    SUBROUTINE production_e_ij( i, j )
509
510       USE arrays_3d
511       USE cloud_parameters
512       USE control_parameters
513       USE grid_variables
514       USE indices
515       USE statistics
516
517       IMPLICIT NONE
518
519       INTEGER ::  i, j, k
520
521       REAL    ::  def, dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz, &
522                   k1, k2, theta, temp
523
524       REAL, DIMENSION(nzb:nzt+1) ::  usvs, vsus, wsus, wsvs
525
526!
527!--    Calculate TKE production by shear
528       DO  k = nzb_diff_s_outer(j,i), nzt
529
530          dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
531          dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
532                           u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
533          dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
534                           u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
535
536          dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
537                           v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
538          dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
539          dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
540                           v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
541
542          dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
543                           w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
544          dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
545                           w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
546          dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
547
548          def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 )                       &
549                + dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + dvdz**2 &
550                + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
551
552          IF ( def < 0.0 )  def = 0.0
553
554          tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
555               
556       ENDDO
557
558       IF ( prandtl_layer )  THEN
559
560          IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) )  THEN
561
562!
563!--          Position beneath wall
564!--          (2) - Will allways be executed.
565!--          'bottom and wall: use u_0,v_0 and wall functions'
566             k = nzb_diff_s_inner(j,i)-1
567
568             dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
569             dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
570                            u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
571             dvdy =       ( v(k,j+1,i) - v(k,j,i)     ) * ddy
572             dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
573                            v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
574             dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
575
576             IF ( wall_e_y(j,i) /= 0.0 )  THEN
577                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
578                                    usvs, 1.0, 0.0, 0.0, 0.0 )
579                dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i)
580                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
581                                    wsvs, 0.0, 0.0, 1.0, 0.0 )
582                dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i)
583             ELSE
584                dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
585                                u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
586                dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
587                                w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
588             ENDIF
589
590             IF ( wall_e_x(j,i) /= 0.0 )  THEN
591                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
592                                    vsus, 0.0, 1.0, 0.0, 0.0 )
593                dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i)
594                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
595                                    wsus, 0.0, 0.0, 0.0, 1.0 )
596                dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i)
597             ELSE
598                dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
599                                v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
600                dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
601                                w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
602             ENDIF
603
604             def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
605                   dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
606                   dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
607
608             IF ( def < 0.0 )  def = 0.0
609
610             tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
611
612!
613!--          (3) - will be executed only, if there is at least one level
614!--          between (2) and (4), i.e. the topography must have a
615!--          minimum height of 2 dz. Wall fluxes for this case have
616!--          already been calculated for (2).
617!--          'wall only: use wall functions'
618             DO  k = nzb_diff_s_inner(j,i), nzb_diff_s_outer(j,i)-2
619
620                dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
621                dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
622                               u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
623                dvdy =       ( v(k,j+1,i) - v(k,j,i)     ) * ddy
624                dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
625                               v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
626                dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
627
628                IF ( wall_e_y(j,i) /= 0.0 )  THEN
629                   dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i)
630                   dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i)
631                ELSE
632                   dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
633                                   u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
634                   dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
635                                   w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
636                ENDIF
637
638                IF ( wall_e_x(j,i) /= 0.0 )  THEN
639                   dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i)
640                   dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i)
641                ELSE
642                   dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
643                                   v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
644                   dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
645                                   w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
646                ENDIF
647
648                def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
649                      dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
650                      dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
651
652                IF ( def < 0.0 )  def = 0.0
653
654                tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
655
656             ENDDO
657
658!
659!--          (4) - will allways be executed.
660!--          'special case: free atmosphere' (as for case (0))
661             k = nzb_diff_s_outer(j,i)-1
662
663             dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
664             dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
665                              u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
666             dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
667                              u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
668
669             dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
670                              v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
671             dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
672             dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
673                              v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
674
675             dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
676                              w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
677             dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
678                              w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
679             dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
680
681             def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
682                   dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
683                   dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
684
685             IF ( def < 0.0 )  def = 0.0
686
687             tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
688
689          ELSE
690
691!
692!--          Position without adjacent wall
693!--          (1) - will allways be executed.
694!--          'bottom only: use u_0,v_0'
695             k = nzb_diff_s_inner(j,i)-1
696
697             dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
698             dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
699                              u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
700             dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
701                              u_0(j,i)   - u_0(j,i+1)   ) * dd2zu(k)
702
703             dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
704                              v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
705             dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
706             dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
707                              v_0(j,i)   - v_0(j+1,i)   ) * dd2zu(k)
708
709             dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
710                              w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
711             dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
712                              w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
713             dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
714
715             def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 )                       &
716                   + dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + dvdz**2 &
717                   + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
718
719             IF ( def < 0.0 )  def = 0.0
720
721             tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
722
723          ENDIF
724
725       ELSEIF ( use_surface_fluxes )  THEN
726
727          k = nzb_diff_s_outer(j,i)-1
728
729          dudx  =        ( u(k,j,i+1) - u(k,j,i)     ) * ddx
730          dudy  = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
731                           u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
732          dudz  = 0.5  * ( u(k+1,j,i) + u(k+1,j,i+1) - &
733                           u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
734
735          dvdx  = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
736                           v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
737          dvdy  =        ( v(k,j+1,i) - v(k,j,i)     ) * ddy
738          dvdz  = 0.5  * ( v(k+1,j,i) + v(k+1,j+1,i) - &
739                           v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
740
741          dwdx  = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
742                           w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
743          dwdy  = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
744                           w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
745          dwdz  =        ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
746
747          def = 2.0 * ( dudx**2 + dvdy**2 + dwdz**2 ) +           &
748                dudy**2 + dvdx**2 + dwdx**2 + dwdy**2 + dudz**2 + &
749                dvdz**2 + 2.0 * ( dvdx*dudy + dwdx*dudz + dwdy*dvdz )
750
751          IF ( def < 0.0 )  def = 0.0
752
753          tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
754
755       ENDIF
756
757!
758!--    Calculate TKE production by buoyancy
759       IF ( .NOT. humidity )  THEN
760
761          IF ( use_pt_reference )  THEN
762
763             DO  k = nzb_diff_s_inner(j,i), nzt_diff
764                tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt_reference * &
765                                       ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
766             ENDDO
767
768             IF ( use_surface_fluxes )  THEN
769                k = nzb_diff_s_inner(j,i)-1
770                tend(k,j,i) = tend(k,j,i) + g / pt_reference * shf(j,i)
771             ENDIF
772
773             IF ( use_top_fluxes )  THEN
774                k = nzt
775                tend(k,j,i) = tend(k,j,i) + g / pt_reference * tswst(j,i)
776             ENDIF
777
778          ELSE
779
780             DO  k = nzb_diff_s_inner(j,i), nzt_diff
781                tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / pt(k,j,i) * &
782                                       ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
783             ENDDO
784
785             IF ( use_surface_fluxes )  THEN
786                k = nzb_diff_s_inner(j,i)-1
787                tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
788             ENDIF
789
790             IF ( use_top_fluxes )  THEN
791                k = nzt
792                tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
793             ENDIF
794
795         ENDIF
796
797       ELSE
798
799          DO  k = nzb_diff_s_inner(j,i), nzt_diff
800
801             IF ( .NOT. cloud_physics )  THEN
802                k1 = 1.0 + 0.61 * q(k,j,i)
803                k2 = 0.61 * pt(k,j,i)
804             ELSE
805                IF ( ql(k,j,i) == 0.0 )  THEN
806                   k1 = 1.0 + 0.61 * q(k,j,i)
807                   k2 = 0.61 * pt(k,j,i)
808                ELSE
809                   theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
810                   temp  = theta * t_d_pt(k)
811                   k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
812                              ( q(k,j,i) - ql(k,j,i) ) *          &
813                        ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
814                        ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
815                        ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
816                   k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
817                ENDIF
818             ENDIF
819
820             tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) *      &
821                                      ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
822                                        k2 * ( q(k+1,j,i) - q(k-1,j,i) )   &
823                                      ) * dd2zu(k)
824          ENDDO
825
826          IF ( use_surface_fluxes )  THEN
827             k = nzb_diff_s_inner(j,i)-1
828
829             IF ( .NOT. cloud_physics ) THEN
830                k1 = 1.0 + 0.61 * q(k,j,i)
831                k2 = 0.61 * pt(k,j,i)
832             ELSE
833                IF ( ql(k,j,i) == 0.0 )  THEN
834                   k1 = 1.0 + 0.61 * q(k,j,i)
835                   k2 = 0.61 * pt(k,j,i)
836                ELSE
837                   theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
838                   temp  = theta * t_d_pt(k)
839                   k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
840                              ( q(k,j,i) - ql(k,j,i) ) *          &
841                        ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
842                        ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
843                        ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
844                   k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
845                ENDIF
846             ENDIF
847
848             tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
849                                         ( k1* shf(j,i) + k2 * qsws(j,i) )
850          ENDIF
851
852          IF ( use_top_fluxes )  THEN
853             k = nzt
854
855             IF ( .NOT. cloud_physics ) THEN
856                k1 = 1.0 + 0.61 * q(k,j,i)
857                k2 = 0.61 * pt(k,j,i)
858             ELSE
859                IF ( ql(k,j,i) == 0.0 )  THEN
860                   k1 = 1.0 + 0.61 * q(k,j,i)
861                   k2 = 0.61 * pt(k,j,i)
862                ELSE
863                   theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
864                   temp  = theta * t_d_pt(k)
865                   k1 = ( 1.0 - q(k,j,i) + 1.61 *                 &
866                              ( q(k,j,i) - ql(k,j,i) ) *          &
867                        ( 1.0 + 0.622 * l_d_r / temp ) ) /        &
868                        ( 1.0 + 0.622 * l_d_r * l_d_cp *          &
869                        ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
870                   k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
871                ENDIF
872             ENDIF
873
874             tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
875                                         ( k1* tswst(j,i) + k2 * qswst(j,i) )
876          ENDIF
877
878       ENDIF
879
880    END SUBROUTINE production_e_ij
881
882
883    SUBROUTINE production_e_init
884
885       USE arrays_3d
886       USE control_parameters
887       USE grid_variables
888       USE indices
889
890       IMPLICIT NONE
891
892       INTEGER ::  i, j, ku, kv
893
894       IF ( prandtl_layer )  THEN
895
896          IF ( first_call )  THEN
897             ALLOCATE( u_0(nys-1:nyn+1,nxl-1:nxr+1), &
898                       v_0(nys-1:nyn+1,nxl-1:nxr+1) )
899             first_call = .FALSE.
900          ENDIF
901
902!
903!--       Calculate a virtual velocity at the surface in a way that the
904!--       vertical velocity gradient at k = 1 (u(k+1)-u_0) matches the
905!--       Prandtl law (-w'u'/km). This gradient is used in the TKE shear
906!--       production term at k=1 (see production_e_ij).
907!--       The velocity gradient has to be limited in case of too small km
908!--       (otherwise the timestep may be significantly reduced by large
909!--       surface winds).
910!--       WARNING: the exact analytical solution would require the determination
911!--                of the eddy diffusivity by km = u* * kappa * zp / phi_m.
912          !$OMP PARALLEL DO PRIVATE( ku, kv )
913          DO  i = nxl, nxr
914             DO  j = nys, nyn
915
916                ku = nzb_u_inner(j,i)+1
917                kv = nzb_v_inner(j,i)+1
918
919                u_0(j,i) = u(ku+1,j,i) + usws(j,i) * ( zu(ku+1) - zu(ku-1) ) / &
920                                 ( 0.5 * ( km(ku,j,i) + km(ku,j,i-1) ) +       &
921                                   1.0E-20 )
922!                                  ( us(j,i) * kappa * zu(1) )
923                v_0(j,i) = v(kv+1,j,i) + vsws(j,i) * ( zu(kv+1) - zu(kv-1) ) / &
924                                 ( 0.5 * ( km(kv,j,i) + km(kv,j-1,i) ) +       &
925                                   1.0E-20 )
926!                                  ( us(j,i) * kappa * zu(1) )
927
928                IF ( ABS( u(ku+1,j,i) - u_0(j,i) )  > &
929                     ABS( u(ku+1,j,i) - u(ku-1,j,i) ) )  u_0(j,i) = u(ku-1,j,i)
930                IF ( ABS( v(kv+1,j,i) - v_0(j,i) )  > &
931                     ABS( v(kv+1,j,i) - v(kv-1,j,i) ) )  v_0(j,i) = v(kv-1,j,i)
932
933             ENDDO
934          ENDDO
935
936          CALL exchange_horiz_2d( u_0 )
937          CALL exchange_horiz_2d( v_0 )
938
939       ENDIF
940
941    END SUBROUTINE production_e_init
942
943 END MODULE production_e_mod
Note: See TracBrowser for help on using the repository browser.