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