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

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

vorlaeufige Standalone-Version fuer Linux-Cluster

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