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

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

preliminary version, several changes to be explained later

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