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

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

preliminary version of modified boundary conditions at top

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