source: palm/trunk/SOURCE/flow_statistics.f90 @ 3

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

RCS Log replace by Id keyword, revision history cleaned up

File size: 33.5 KB
Line 
1 SUBROUTINE flow_statistics
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id$
11! RCS Log replace by Id keyword, revision history cleaned up
12!
13! Revision 1.41  2006/08/04 14:37:50  raasch
14! Error removed in non-parallel part (sums_l)
15!
16! Revision 1.1  1997/08/11 06:15:17  raasch
17! Initial revision
18!
19!
20! Description:
21! ------------
22! Compute average profiles and further average flow quantities for the different
23! user-defined (sub-)regions. The region indexed 0 is the total model domain.
24!
25! NOTE: For simplicity, nzb_s_outer and nzb_diff_s_outer are being used as a
26! ----  lower vertical index for k-loops for all variables so that regardless
27! of the variable and its respective staggered grid always the same number of
28! grid points is used for 2D averages. The disadvantage: depending on the
29! variable, up to one grid layer adjacent to the (vertical walls of the)
30! topography is missed out by this simplification.
31!------------------------------------------------------------------------------!
32
33    USE arrays_3d
34    USE cloud_parameters
35    USE cpulog
36    USE grid_variables
37    USE indices
38    USE interfaces
39    USE pegrid
40    USE statistics
41    USE control_parameters
42
43    IMPLICIT NONE
44
45    INTEGER ::  i, j, k, omp_get_thread_num, sr, tn
46    LOGICAL ::  first
47    REAL    ::  height, pts, sums_l_eper, sums_l_etot, ust, ust2, u2, vst, &
48                vst2, v2, w2, z_i(2)
49    REAL    ::  sums_ll(nzb:nzt+1,2)
50
51
52    CALL cpu_log( log_point(10), 'flow_statistics', 'start' )
53
54!
55!-- To be on the safe side, check whether flow_statistics has already been
56!-- called once after the current time step
57    IF ( flow_statistics_called )  THEN
58       IF ( myid == 0 )  PRINT*, '+++ WARNING: flow_statistics is called two', &
59                                 ' times within one timestep'
60       CALL local_stop
61    ENDIF
62
63!
64!-- Compute statistics for each (sub-)region
65    DO  sr = 0, statistic_regions
66
67!
68!--    Initialize (local) summation array
69       sums_l = 0.0
70
71!
72!--    Store sums that have been computed in other subroutines in summation
73!--    array
74       sums_l(:,11,:) = sums_l_l(:,sr,:)      ! mixing length from diffusivities
75!--    WARNING: next line still has to be adjusted for OpenMP
76       sums_l(:,21,0) = sums_wsts_bc_l(:,sr)  ! heat flux from advec_s_bc
77       sums_l(nzb+9,var_sum,0)  = sums_divold_l(sr)  ! old divergence from pres
78       sums_l(nzb+10,var_sum,0) = sums_divnew_l(sr)  ! new divergence from pres
79!--    WARNING: next four lines still may have to be adjusted for OpenMP
80       sums_l(nzb:nzb+2,var_sum-1,0)    = sums_up_fraction_l(1,1:3,sr)! upstream
81       sums_l(nzb+3:nzb+5,var_sum-1,0)  = sums_up_fraction_l(2,1:3,sr)! parts
82       sums_l(nzb+6:nzb+8,var_sum-1,0)  = sums_up_fraction_l(3,1:3,sr)! from
83       sums_l(nzb+9:nzb+11,var_sum-1,0) = sums_up_fraction_l(4,1:3,sr)! spline
84
85!
86!--    Horizontally averaged profiles of horizontal velocities and temperature.
87!--    They must have been computed before, because they are already required
88!--    for other horizontal averages.
89       tn = 0
90       !$OMP PARALLEL PRIVATE( i, j, k, tn )
91#if defined( __lcmuk )
92       tn = omp_get_thread_num()
93#else
94!$     tn = omp_get_thread_num()
95#endif
96
97       !$OMP DO
98       DO  i = nxl, nxr
99          DO  j =  nys, nyn
100             DO  k = nzb_s_outer(j,i), nzt+1
101                sums_l(k,1,tn)  = sums_l(k,1,tn)  + u(k,j,i)  * rmask(j,i,sr)
102                sums_l(k,2,tn)  = sums_l(k,2,tn)  + v(k,j,i)  * rmask(j,i,sr)
103                sums_l(k,4,tn)  = sums_l(k,4,tn)  + pt(k,j,i) * rmask(j,i,sr)
104             ENDDO
105          ENDDO
106       ENDDO
107
108!
109!--    Horizontally averaged profiles of virtual potential temperature,
110!--    total water content, specific humidity and liquid water potential
111!--    temperature
112       IF ( moisture )  THEN
113          !$OMP DO
114          DO  i = nxl, nxr
115             DO  j =  nys, nyn
116                DO  k = nzb_s_outer(j,i), nzt+1
117                   sums_l(k,44,tn)  = sums_l(k,44,tn) + &
118                                      vpt(k,j,i) * rmask(j,i,sr)
119                   sums_l(k,41,tn)  = sums_l(k,41,tn) + &
120                                      q(k,j,i) * rmask(j,i,sr)
121                ENDDO
122             ENDDO
123          ENDDO
124          IF ( cloud_physics )  THEN
125             !$OMP DO
126             DO  i = nxl, nxr
127                DO  j =  nys, nyn
128                   DO  k = nzb_s_outer(j,i), nzt+1
129                      sums_l(k,42,tn) = sums_l(k,42,tn) + &
130                                      ( q(k,j,i) - ql(k,j,i) ) * rmask(j,i,sr)
131                      sums_l(k,43,tn) = sums_l(k,43,tn) + ( &
132                                      pt(k,j,i) + l_d_cp*pt_d_t(k) * ql(k,j,i) &
133                                                          ) * rmask(j,i,sr)
134                   ENDDO
135                ENDDO
136             ENDDO
137          ENDIF
138       ENDIF
139
140!
141!--    Horizontally averaged profiles of passive scalar
142       IF ( passive_scalar )  THEN
143          !$OMP DO
144          DO  i = nxl, nxr
145             DO  j =  nys, nyn
146                DO  k = nzb_s_outer(j,i), nzt+1
147                   sums_l(k,41,tn)  = sums_l(k,41,tn) + q(k,j,i) * rmask(j,i,sr)
148                ENDDO
149             ENDDO
150          ENDDO
151       ENDIF
152       !$OMP END PARALLEL
153
154!
155!--    Summation of thread sums
156       IF ( threads_per_task > 1 )  THEN
157          DO  i = 1, threads_per_task-1
158             sums_l(:,1,0) = sums_l(:,1,0) + sums_l(:,1,i)
159             sums_l(:,2,0) = sums_l(:,2,0) + sums_l(:,2,i)
160             sums_l(:,4,0) = sums_l(:,4,0) + sums_l(:,4,i)
161             IF ( moisture )  THEN
162                sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
163                sums_l(:,44,0) = sums_l(:,44,0) + sums_l(:,44,i)
164                IF ( cloud_physics )  THEN
165                   sums_l(:,42,0) = sums_l(:,42,0) + sums_l(:,42,i)
166                   sums_l(:,43,0) = sums_l(:,43,0) + sums_l(:,43,i)
167                ENDIF
168             ENDIF
169             IF ( passive_scalar )  THEN
170                sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
171             ENDIF
172          ENDDO
173       ENDIF
174
175#if defined( __parallel )
176!
177!--    Compute total sum from local sums
178       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, MPI_REAL, &
179                           MPI_SUM, comm2d, ierr )
180       CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, MPI_REAL, &
181                           MPI_SUM, comm2d, ierr )
182       CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, MPI_REAL, &
183                           MPI_SUM, comm2d, ierr )
184       IF ( moisture ) THEN
185          CALL MPI_ALLREDUCE( sums_l(nzb,44,0), sums(nzb,44), nzt+2-nzb, &
186                              MPI_REAL, MPI_SUM, comm2d, ierr )
187          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb, &
188                              MPI_REAL, MPI_SUM, comm2d, ierr )
189          IF ( cloud_physics ) THEN
190             CALL MPI_ALLREDUCE( sums_l(nzb,42,0), sums(nzb,42), nzt+2-nzb, &
191                                 MPI_REAL, MPI_SUM, comm2d, ierr )
192             CALL MPI_ALLREDUCE( sums_l(nzb,43,0), sums(nzb,43), nzt+2-nzb, &
193                                 MPI_REAL, MPI_SUM, comm2d, ierr )
194          ENDIF
195       ENDIF
196
197       IF ( passive_scalar )  THEN
198          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb, &
199                              MPI_REAL, MPI_SUM, comm2d, ierr )
200       ENDIF
201#else
202       sums(:,1) = sums_l(:,1,0)
203       sums(:,2) = sums_l(:,2,0)
204       sums(:,4) = sums_l(:,4,0)
205       IF ( moisture ) THEN
206          sums(:,44) = sums_l(:,44,0)
207          sums(:,41) = sums_l(:,41,0)
208          IF ( cloud_physics ) THEN
209             sums(:,42) = sums_l(:,42,0)
210             sums(:,43) = sums_l(:,43,0)
211          ENDIF
212       ENDIF
213       IF ( passive_scalar )  sums(:,41) = sums_l(:,41,0)
214#endif
215
216!
217!--    Final values are obtained by division by the total number of grid points
218!--    used for summation. After that store profiles.
219       sums(:,1) = sums(:,1) / ngp_2dh_outer(:,sr)
220       sums(:,2) = sums(:,2) / ngp_2dh_outer(:,sr)
221       sums(:,4) = sums(:,4) / ngp_2dh_outer(:,sr)
222       hom(:,1,1,sr) = sums(:,1)             ! u
223       hom(:,1,2,sr) = sums(:,2)             ! v
224       hom(:,1,4,sr) = sums(:,4)             ! pt
225
226!
227!--    Humidity and cloud parameters
228       IF ( moisture ) THEN
229          sums(:,44) = sums(:,44) / ngp_2dh_outer(:,sr)
230          sums(:,41) = sums(:,41) / ngp_2dh_outer(:,sr)
231          hom(:,1,44,sr) = sums(:,44)             ! vpt
232          hom(:,1,41,sr) = sums(:,41)             ! qv (q)
233          IF ( cloud_physics ) THEN
234             sums(:,42) = sums(:,42) / ngp_2dh_outer(:,sr)
235             sums(:,43) = sums(:,43) / ngp_2dh_outer(:,sr)
236             hom(:,1,42,sr) = sums(:,42)             ! qv
237             hom(:,1,43,sr) = sums(:,43)             ! pt
238          ENDIF
239       ENDIF
240
241!
242!--    Passive scalar
243       IF ( passive_scalar )  hom(:,1,41,sr) = sums(:,41) / ngp_2dh_outer(:,sr)
244
245!
246!--    Horizontally averaged profiles of the remaining prognostic variables,
247!--    variances, the total and the perturbation energy (single values in last
248!--    column of sums_l) and some diagnostic quantities.
249!--    NOTE: for simplicity, nzb_s_outer is used below, although strictly
250!--    ----  speaking the following k-loop would have to be split up and
251!--          rearranged according to the staggered grid.
252       tn = 0
253#if defined( __lcmuk )
254       !$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper, sums_l_etot, &
255       !$OMP                    tn, ust, ust2, u2, vst, vst2, v2, w2 )
256       tn = omp_get_thread_num()
257#else
258       !$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper, sums_l_etot, tn, ust, ust2, u2, vst, vst2, v2, w2 )
259!$     tn = omp_get_thread_num()
260#endif
261       !$OMP DO
262       DO  i = nxl, nxr
263          DO  j =  nys, nyn
264             sums_l_etot = 0.0
265             sums_l_eper = 0.0
266             DO  k = nzb_s_outer(j,i), nzt+1
267                u2   = u(k,j,i)**2
268                v2   = v(k,j,i)**2
269                w2   = w(k,j,i)**2
270                ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2
271                vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2
272!
273!--             Prognostic and diagnostic variables
274                sums_l(k,3,tn)  = sums_l(k,3,tn)  + w(k,j,i)  * rmask(j,i,sr)
275                sums_l(k,8,tn)  = sums_l(k,8,tn)  + e(k,j,i)  * rmask(j,i,sr)
276                sums_l(k,9,tn)  = sums_l(k,9,tn)  + km(k,j,i) * rmask(j,i,sr)
277                sums_l(k,10,tn) = sums_l(k,10,tn) + kh(k,j,i) * rmask(j,i,sr)
278                sums_l(k,40,tn) = sums_l(k,40,tn) + p(k,j,i)
279
280!
281!--             Variances
282                sums_l(k,30,tn) = sums_l(k,30,tn) + ust2 * rmask(j,i,sr)
283                sums_l(k,31,tn) = sums_l(k,31,tn) + vst2 * rmask(j,i,sr)
284                sums_l(k,32,tn) = sums_l(k,32,tn) + w2   * rmask(j,i,sr)
285                sums_l(k,33,tn) = sums_l(k,33,tn) + &
286                                  ( pt(k,j,i)-hom(k,1,4,sr) )**2 * rmask(j,i,sr)
287!
288!--             Higher moments
289!--             (Computation of the skewness of w further below)
290                sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i) * w2 * &
291                                                    rmask(j,i,sr)
292!
293!--             Perturbation energy
294                sums_l(k,34,tn) = sums_l(k,34,tn) + 0.5 * ( ust2 + vst2 + w2 ) &
295                                                    * rmask(j,i,sr)
296                sums_l_etot  = sums_l_etot + &
297                                        0.5 * ( u2 + v2 + w2 ) * rmask(j,i,sr)
298                sums_l_eper  = sums_l_eper + &
299                                        0.5 * ( ust2+vst2+w2 ) * rmask(j,i,sr)
300             ENDDO
301!
302!--          Total and perturbation energy for the total domain (being
303!--          collected in the last column of sums_l). Summation of these
304!--          quantities is seperated from the previous loop in order to
305!--          allow vectorization of that loop.
306             sums_l(nzb+4,var_sum,tn) = sums_l(nzb+4,var_sum,tn) + sums_l_etot
307             sums_l(nzb+5,var_sum,tn) = sums_l(nzb+5,var_sum,tn) + sums_l_eper
308!
309!--          2D-arrays (being collected in the last column of sums_l)
310             sums_l(nzb,var_sum,tn)   = sums_l(nzb,var_sum,tn) +   &
311                                        us(j,i)   * rmask(j,i,sr)
312             sums_l(nzb+1,var_sum,tn) = sums_l(nzb+1,var_sum,tn) + &
313                                        usws(j,i) * rmask(j,i,sr)
314             sums_l(nzb+2,var_sum,tn) = sums_l(nzb+2,var_sum,tn) + &
315                                        vsws(j,i) * rmask(j,i,sr)
316             sums_l(nzb+3,var_sum,tn) = sums_l(nzb+3,var_sum,tn) + &
317                                        ts(j,i)   * rmask(j,i,sr)
318          ENDDO
319       ENDDO
320
321!
322!--    Horizontally averaged profiles of the vertical fluxes
323       !$OMP DO
324       DO  i = nxl, nxr
325          DO  j = nys, nyn
326!
327!--          Subgridscale fluxes (without Prandtl layer from k=nzb,
328!--          oterwise from k=nzb+1)
329!--          NOTE: for simplicity, nzb_diff_s_outer is used below, although
330!--          ----  strictly speaking the following k-loop would have to be
331!--                split up according to the staggered grid.
332             DO  k = nzb_diff_s_outer(j,i)-1, nzt
333!
334!--             Momentum flux w"u"
335                sums_l(k,12,tn) = sums_l(k,12,tn) - 0.25 * (                   &
336                               km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) &
337                                                           ) * (               &
338                                   ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
339                                 + ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
340                                                               ) * rmask(j,i,sr)
341!
342!--             Momentum flux w"v"
343                sums_l(k,14,tn) = sums_l(k,14,tn) - 0.25 * (                   &
344                               km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) &
345                                                           ) * (               &
346                                   ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
347                                 + ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
348                                                               ) * rmask(j,i,sr)
349!
350!--             Heat flux w"pt"
351                sums_l(k,16,tn) = sums_l(k,16,tn)                              &
352                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
353                                               * ( pt(k+1,j,i) - pt(k,j,i) )   &
354                                               * ddzu(k+1) * rmask(j,i,sr)
355
356
357!
358!--             Buoyancy flux, water flux (humidity flux) w"q"
359                IF ( moisture ) THEN
360                   sums_l(k,45,tn) = sums_l(k,45,tn)                           &
361                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
362                                               * ( vpt(k+1,j,i) - vpt(k,j,i) ) &
363                                               * ddzu(k+1) * rmask(j,i,sr)
364                   sums_l(k,48,tn) = sums_l(k,48,tn)                           &
365                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
366                                               * ( q(k+1,j,i) - q(k,j,i) )     &
367                                               * ddzu(k+1) * rmask(j,i,sr)
368                   IF ( cloud_physics ) THEN
369                      sums_l(k,51,tn) = sums_l(k,51,tn)                        &
370                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
371                                               * ( ( q(k+1,j,i) - ql(k+1,j,i) )&
372                                                - ( q(k,j,i) - ql(k,j,i) ) )   &
373                                               * ddzu(k+1) * rmask(j,i,sr) 
374                   ENDIF
375                ENDIF
376
377!
378!--             Passive scalar flux
379                IF ( passive_scalar )  THEN
380                   sums_l(k,48,tn) = sums_l(k,48,tn)                           &
381                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
382                                               * ( q(k+1,j,i) - q(k,j,i) )     &
383                                               * ddzu(k+1) * rmask(j,i,sr)
384                ENDIF
385
386             ENDDO
387
388!
389!--          Subgridscale fluxes in the Prandtl layer
390             IF ( use_surface_fluxes )  THEN
391                sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + &
392                                    usws(j,i) * rmask(j,i,sr)     ! w"u"
393                sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + &
394                                    vsws(j,i) * rmask(j,i,sr)     ! w"v"
395                sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + &
396                                    shf(j,i)  * rmask(j,i,sr)     ! w"pt"
397                sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + &
398                                    0.0 * rmask(j,i,sr)           ! u"pt"
399                sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + &
400                                    0.0 * rmask(j,i,sr)           ! v"pt"
401                IF ( moisture )  THEN
402                   sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + &
403                                       qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
404                   IF ( cloud_physics )  THEN
405                      sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (           &
406                                          ( 1.0 + 0.61 * q(nzb,j,i) ) *   &
407                                          shf(j,i) + 0.61 * pt(nzb,j,i) * &
408                                                     qsws(j,i)            &
409                                                              )
410!
411!--                   Formula does not work if ql(nzb) /= 0.0
412                      sums_l(nzb,51,tn) = sums_l(nzb,51,tn) + &   ! w"q" (w"qv")
413                                          qsws(j,i) * rmask(j,i,sr)
414                   ENDIF
415                ENDIF
416                IF ( passive_scalar )  THEN
417                   sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + &
418                                       qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
419                ENDIF
420             ENDIF
421
422!
423!--          Resolved fluxes (can be computed for all horizontal points)
424!--          NOTE: for simplicity, nzb_s_outer is used below, although strictly
425!--          ----  speaking the following k-loop would have to be split up and
426!--                rearranged according to the staggered grid.
427             DO  k = nzb_s_outer(j,i), nzt
428                ust = 0.5 * ( u(k,j,i)   - hom(k,1,1,sr) + &
429                              u(k+1,j,i) - hom(k+1,1,1,sr) )
430                vst = 0.5 * ( v(k,j,i)   - hom(k,1,2,sr) + &
431                              v(k+1,j,i) - hom(k+1,1,2,sr) )
432                pts = 0.5 * ( pt(k,j,i)   - hom(k,1,4,sr) + &
433                              pt(k+1,j,i) - hom(k+1,1,4,sr) )
434!
435!--             Momentum flux w*u*
436                sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5 *                     &
437                                                    ( w(k,j,i-1) + w(k,j,i) ) &
438                                                    * ust * rmask(j,i,sr)
439!
440!--             Momentum flux w*v*
441                sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5 *                     &
442                                                    ( w(k,j-1,i) + w(k,j,i) ) &
443                                                    * vst * rmask(j,i,sr)
444!
445!--             Heat flux w*pt*
446!--             The following formula (comment line, not executed) does not
447!--             work if applied to subregions
448!                sums_l(k,17,tn) = sums_l(k,17,tn) + 0.5 *                     &
449!                                                    ( pt(k,j,i)+pt(k+1,j,i) ) &
450!                                                    * w(k,j,i) * rmask(j,i,sr)
451                sums_l(k,17,tn) = sums_l(k,17,tn) + pts * w(k,j,i) * &
452                                                    rmask(j,i,sr)
453!
454!--             Higher moments
455                sums_l(k,35,tn) = sums_l(k,35,tn) + pts * w(k,j,i)**2 * &
456                                                    rmask(j,i,sr)
457                sums_l(k,36,tn) = sums_l(k,36,tn) + pts**2 * w(k,j,i) * &
458                                                    rmask(j,i,sr)
459
460!
461!--             Buoyancy flux, water flux, humidity flux and liquid water
462!--             content
463                IF ( moisture )  THEN
464                   pts = 0.5 * ( vpt(k,j,i)   - hom(k,1,44,sr) + &
465                                 vpt(k+1,j,i) - hom(k+1,1,44,sr) )
466                   sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * &
467                                                       rmask(j,i,sr)
468                   pts = 0.5 * ( q(k,j,i)   - hom(k,1,41,sr) + &
469                                 q(k+1,j,i) - hom(k+1,1,41,sr) )
470                   sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) * &
471                                                       rmask(j,i,sr)
472                   IF ( cloud_physics  .OR.  cloud_droplets )  THEN
473                      pts = 0.5 *                                            &
474                             ( ( q(k,j,i)   - ql(k,j,i)   ) - hom(k,1,42,sr) &
475                             + ( q(k+1,j,i) - ql(k+1,j,i) ) - hom(k+1,1,42,sr) )
476                      sums_l(k,52,tn) = sums_l(k,52,tn) + pts * w(k,j,i) * &
477                                                          rmask(j,i,sr)
478                      sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * &
479                                                          rmask(j,i,sr)
480                   ENDIF
481                ENDIF
482
483!
484!--             Passive scalar flux
485                IF ( passive_scalar )  THEN
486                   pts = 0.5 * ( q(k,j,i)   - hom(k,1,41,sr) + &
487                                 q(k+1,j,i) - hom(k+1,1,41,sr) )
488                   sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) * &
489                                                       rmask(j,i,sr)
490                ENDIF
491
492!
493!--             Energy flux w*e*
494                sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5 *           &
495                                              ( ust**2 + vst**2 + w(k,j,i)**2 )&
496                                              * rmask(j,i,sr)
497         
498             ENDDO
499          ENDDO
500       ENDDO
501
502!
503!--    Divergence of vertical flux of resolved scale energy and pressure
504!--    fluctuations. First calculate the products, then the divergence.
505!--    Calculation is time consuming. Do it only, if profiles shall be plotted.
506       IF ( hom(nzb+1,2,55,0) /= 0.0 )  THEN
507
508          sums_ll = 0.0  ! local array
509
510          !$OMP DO
511          DO  i = nxl, nxr
512             DO  j = nys, nyn
513                DO  k = nzb_s_outer(j,i)+1, nzt
514
515                   sums_ll(k,1) = sums_ll(k,1) + 0.5 * w(k,j,i) * (        &
516                  ( 0.25 * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1)   &
517                           - 2.0 * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) )     &
518                           ) )**2                                          &
519                + ( 0.25 * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i)   &
520                           - 2.0 * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) )     &
521                           ) )**2                                          &
522                   + w(k,j,i)**2                                  )
523
524                   sums_ll(k,2) = sums_ll(k,2) + 0.5 * w(k,j,i) &
525                                               * ( p(k,j,i) + p(k+1,j,i) )
526
527                ENDDO
528             ENDDO
529          ENDDO
530          sums_ll(0,1)     = 0.0    ! because w is zero at the bottom
531          sums_ll(nzt+1,1) = 0.0
532          sums_ll(0,2)     = 0.0
533          sums_ll(nzt+1,2) = 0.0
534
535          DO  k = nzb_s_outer(j,i)+1, nzt
536             sums_l(k,55,tn) = ( sums_ll(k,1) - sums_ll(k-1,1) ) * ddzw(k)
537             sums_l(k,56,tn) = ( sums_ll(k,2) - sums_ll(k-1,2) ) * ddzw(k)
538          ENDDO
539          sums_l(nzb,55,tn) = sums_l(nzb+1,55,tn)
540          sums_l(nzb,56,tn) = sums_l(nzb+1,56,tn)
541
542       ENDIF
543
544!
545!--    Divergence of vertical flux of SGS TKE
546       IF ( hom(nzb+1,2,57,0) /= 0.0 )  THEN
547
548          !$OMP DO
549          DO  i = nxl, nxr
550             DO  j = nys, nyn
551                DO  k = nzb_s_outer(j,i)+1, nzt
552
553                   sums_l(k,57,tn) = sums_l(k,57,tn) + (                       &
554                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
555                 - (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k)   &
556                                                 ) * ddzw(k)
557
558                ENDDO
559             ENDDO
560          ENDDO
561          sums_l(nzb,57,tn) = sums_l(nzb+1,57,tn)
562
563       ENDIF
564
565!
566!--    Horizontal heat fluxes (subgrid, resolved, total).
567!--    Do it only, if profiles shall be plotted.
568       IF ( hom(nzb+1,2,58,0) /= 0.0 ) THEN
569
570          !$OMP DO
571          DO  i = nxl, nxr
572             DO  j = nys, nyn
573                DO  k = nzb_s_outer(j,i)+1, nzt
574!
575!--                Subgrid horizontal heat fluxes u"pt", v"pt"
576                   sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5 *                   &
577                                                   ( kh(k,j,i) + kh(k,j,i-1) ) &
578                                                 * ( pt(k,j,i-1) - pt(k,j,i) ) &
579                                                 * ddx * rmask(j,i,sr)
580                   sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5 *                   &
581                                                   ( kh(k,j,i) + kh(k,j-1,i) ) &
582                                                 * ( pt(k,j-1,i) - pt(k,j,i) ) &
583                                                 * ddy * rmask(j,i,sr)
584!
585!--                Resolved horizontal heat fluxes u*pt*, v*pt*
586                   sums_l(k,59,tn) = sums_l(k,59,tn) +                         &
587                                                  ( u(k,j,i) - hom(k,1,1,sr) ) &
588                                       * 0.5 * ( pt(k,j,i-1) - hom(k,1,4,sr) + &
589                                                 pt(k,j,i)   - hom(k,1,4,sr) )
590                   pts = 0.5 * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
591                                 pt(k,j,i)   - hom(k,1,4,sr) )
592                   sums_l(k,62,tn) = sums_l(k,62,tn) +                         &
593                                                  ( v(k,j,i) - hom(k,1,2,sr) ) &
594                                       * 0.5 * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
595                                                 pt(k,j,i)   - hom(k,1,4,sr) )
596                ENDDO
597             ENDDO
598          ENDDO
599!
600!--       Fluxes at the surface must be zero (e.g. due to the Prandtl-layer)
601          sums(nzb,58) = 0.0
602          sums(nzb,59) = 0.0
603          sums(nzb,60) = 0.0
604          sums(nzb,61) = 0.0
605          sums(nzb,62) = 0.0
606          sums(nzb,63) = 0.0
607
608       ENDIF
609       !$OMP END PARALLEL
610
611!
612!--    Summation of thread sums
613       IF ( threads_per_task > 1 )  THEN
614          DO  i = 1, threads_per_task-1
615             sums_l(:,3,0)          = sums_l(:,3,0) + sums_l(:,3,i)
616             sums_l(:,4:40,0)       = sums_l(:,4:40,0) + sums_l(:,4:40,i)
617             sums_l(:,45:var_sum,0) = sums_l(:,45:var_sum,0) + &
618                                      sums_l(:,45:var_sum,i)
619          ENDDO
620       ENDIF
621
622#if defined( __parallel )
623!
624!--    Compute total sum from local sums
625       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), ngp_sums, MPI_REAL, &
626                           MPI_SUM, comm2d, ierr )
627#else
628       sums = sums_l(:,:,0)
629#endif
630
631!
632!--    Final values are obtained by division by the total number of grid points
633!--    used for summation. After that store profiles.
634!--    Profiles:
635       DO  k = nzb, nzt+1
636          sums(k,:var_sum-2)      = sums(k,:var_sum-2) / ngp_2dh_outer(k,sr)
637       ENDDO
638!--    Upstream-parts
639       sums(nzb:nzb+11,var_sum-1) = sums(nzb:nzb+11,var_sum-1) / ngp_3d(sr)
640!--    u* and so on
641!--    As sums(nzb:nzb+3,var_sum) are full 2D arrays (us, usws, vsws, ts) whose
642!--    size is always ( nx + 1 ) * ( ny + 1 ), defined at the first grid layer
643!--    above the topography, they are being divided by ngp_2dh(sr)
644       sums(nzb:nzb+3,var_sum)    = sums(nzb:nzb+3,var_sum)    / &
645                                    ngp_2dh(sr)
646!--    eges, e*
647       sums(nzb+4:nzb+5,var_sum)  = sums(nzb+4:nzb+5,var_sum)  / &
648                                    ngp_3d_inner(sr)
649!--    Old and new divergence
650       sums(nzb+9:nzb+10,var_sum) = sums(nzb+9:nzb+10,var_sum) / &
651                                    ngp_3d_inner(sr)
652
653!
654!--    Collect horizontal average in hom.
655!--    Compute deduced averages (e.g. total heat flux)
656       hom(:,1,3,sr)  = sums(:,3)      ! w
657       hom(:,1,8,sr)  = sums(:,8)      ! e     profiles 5-7 are initial profiles
658       hom(:,1,9,sr)  = sums(:,9)      ! km
659       hom(:,1,10,sr) = sums(:,10)     ! kh
660       hom(:,1,11,sr) = sums(:,11)     ! l
661       hom(:,1,12,sr) = sums(:,12)     ! w"u"
662       hom(:,1,13,sr) = sums(:,13)     ! w*u*
663       hom(:,1,14,sr) = sums(:,14)     ! w"v"
664       hom(:,1,15,sr) = sums(:,15)     ! w*v*
665       hom(:,1,16,sr) = sums(:,16)     ! w"pt"
666       hom(:,1,17,sr) = sums(:,17)     ! w*pt*
667       hom(:,1,18,sr) = sums(:,16) + sums(:,17)    ! wpt
668       hom(:,1,19,sr) = sums(:,12) + sums(:,13)    ! wu
669       hom(:,1,20,sr) = sums(:,14) + sums(:,15)    ! wv
670       hom(:,1,21,sr) = sums(:,21)     ! w*pt*BC
671       hom(:,1,22,sr) = sums(:,16) + sums(:,21)    ! wptBC
672                                       ! profiles 23-29 left empty for initial
673                                       ! profiles
674       hom(:,1,30,sr) = sums(:,30)     ! u*2
675       hom(:,1,31,sr) = sums(:,31)     ! v*2
676       hom(:,1,32,sr) = sums(:,32)     ! w*2
677       hom(:,1,33,sr) = sums(:,33)     ! pt*2
678       hom(:,1,34,sr) = sums(:,34)     ! e*
679       hom(:,1,35,sr) = sums(:,35)     ! w*2pt*
680       hom(:,1,36,sr) = sums(:,36)     ! w*pt*2
681       hom(:,1,37,sr) = sums(:,37)     ! w*e*
682       hom(:,1,38,sr) = sums(:,38)     ! w*3
683       hom(:,1,39,sr) = sums(:,38) / ( sums(:,32) + 1E-20 )**1.5    ! Sw
684       hom(:,1,40,sr) = sums(:,40)     ! p
685       hom(:,1,45,sr) = sums(:,45)     ! w"q"
686       hom(:,1,46,sr) = sums(:,46)     ! w*vpt*       
687       hom(:,1,47,sr) = sums(:,45) + sums(:,46)    ! wvpt
688       hom(:,1,48,sr) = sums(:,48)     ! w"q" (w"qv")
689       hom(:,1,49,sr) = sums(:,49)     ! w*q* (w*qv*)
690       hom(:,1,50,sr) = sums(:,48) + sums(:,49)    ! wq (wqv)
691       hom(:,1,51,sr) = sums(:,51)     ! w"qv"
692       hom(:,1,52,sr) = sums(:,52)     ! w*qv*       
693       hom(:,1,53,sr) = sums(:,52) + sums(:,51)    ! wq (wqv)
694       hom(:,1,54,sr) = sums(:,54)     ! ql
695       hom(:,1,55,sr) = sums(:,55)     ! w*u*u*/dz
696       hom(:,1,56,sr) = sums(:,56)     ! w*p*/dz
697       hom(:,1,57,sr) = sums(:,57)     ! w"e/dz
698       hom(:,1,58,sr) = sums(:,58)     ! u"pt"
699       hom(:,1,59,sr) = sums(:,59)     ! u*pt*
700       hom(:,1,60,sr) = sums(:,58) + sums(:,59)    ! upt_t
701       hom(:,1,61,sr) = sums(:,61)     ! v"pt"
702       hom(:,1,62,sr) = sums(:,62)     ! v*pt*
703       hom(:,1,63,sr) = sums(:,61) + sums(:,62)    ! vpt_t
704
705       hom(:,1,var_hom-1,sr) = sums(:,var_sum-1)
706                                       ! upstream-parts u_x, u_y, u_z, v_x,
707                                       ! v_y, usw. (in last but one profile)
708       hom(:,1,var_hom,sr) = sums(:,var_sum) 
709                                       ! u*, w'u', w'v', t* (in last profile)
710
711!
712!--    Determine the boundary layer height using two different schemes.
713!--    First scheme: Starting from the Earth's surface, look for the first
714!--    relative minimum of the total heat flux. The corresponding height is
715!--    accepted as the boundary layer height, if it is less than 1.5 times the
716!--    height where the heat flux becomes negative for the first time.
717!--    NOTE: This criterion is still capable of improving!
718       z_i(1) = 0.0
719       first = .TRUE.
720       DO  k = nzb, nzt-1
721          IF ( first .AND. hom(k,1,18,sr) < 0.0 )  THEN
722             first = .FALSE.
723             height = zw(k)
724          ENDIF
725          IF ( hom(k,1,18,sr) < 0.0  .AND. &
726               hom(k+1,1,18,sr) > hom(k,1,18,sr) )  THEN
727             IF ( zw(k) < 1.5 * height )  THEN
728                z_i(1) = zw(k)
729             ELSE
730                z_i(1) = height
731             ENDIF
732             EXIT
733          ENDIF
734       ENDDO
735
736!
737!--    Second scheme: Starting from the top model boundary, look for the first
738!--    characteristic kink in the temperature profile, where the originally
739!--    stable stratification notably weakens.
740       z_i(2) = 0.0
741       DO  k = nzt-1, nzb+1, -1
742          IF ( ( hom(k+1,1,4,sr) - hom(k,1,4,sr) ) > &
743               2.0 * ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) )  THEN
744             z_i(2) = zu(k)
745             EXIT
746          ENDIF
747       ENDDO
748
749       hom(nzb+6,1,var_hom,sr) = z_i(1)
750       hom(nzb+7,1,var_hom,sr) = z_i(2)
751
752!
753!--    Computation of both the characteristic vertical velocity and
754!--    the characteristic convective boundary layer temperature.
755!--    The horizontal average at nzb+1 is input for the average temperature.
756       IF ( hom(nzb,1,18,sr) > 0.0  .AND.  z_i(1) /= 0.0 )  THEN
757          hom(nzb+8,1,var_hom,sr)  = ( g / hom(nzb+1,1,4,sr) * &
758                                       hom(nzb,1,18,sr) * z_i(1) )**0.333333333
759!--       so far this only works if Prandtl layer is used
760          hom(nzb+11,1,var_hom,sr) = hom(nzb,1,16,sr) / hom(nzb+8,1,var_hom,sr)
761       ELSE
762          hom(nzb+8,1,var_hom,sr)  = 0.0
763          hom(nzb+11,1,var_hom,sr) = 0.0
764       ENDIF
765
766    ENDDO    ! loop of the subregions
767
768!
769!-- Calculate additional statistics provided by the user interface
770    CALL user_statistics
771
772!
773!-- If required, sum up horizontal averages for subsequent time averaging
774    IF ( do_sum )  THEN
775       IF ( average_count_pr == 0 )  hom_sum = 0.0
776       hom_sum = hom_sum + hom(:,1,:,:)
777       average_count_pr = average_count_pr + 1
778       do_sum = .FALSE.
779    ENDIF
780
781!
782!-- Set flag for other UPs (e.g. output routines, but also buoyancy).
783!-- This flag is reset after each time step in time_integration.
784    flow_statistics_called = .TRUE.
785
786    CALL cpu_log( log_point(10), 'flow_statistics', 'stop' )
787
788
789 END SUBROUTINE flow_statistics
790
791
792
Note: See TracBrowser for help on using the repository browser.