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

Last change on this file since 70 was 57, checked in by raasch, 18 years ago

preliminary update of further changes, advec_particles is not running!

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