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

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

preliminary uncomplete changes for ocean version

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