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

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

Preliminary update for user defined profiles

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