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

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

Initial repository layout and content

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