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

Last change on this file since 72 was 51, checked in by raasch, 17 years ago

preliminary version, several changes to be explained later

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