Changeset 699 for palm/trunk/SOURCE

Timestamp:

Mar 22, 2011 5:52:22 PM (14 years ago)

Author:

suehring

Message:

Bugfix in calculation of vertical velocity skewness.

File:

• palm/trunk/SOURCE/flow_statistics.f90 (modified) (5 diffs)

palm/trunk/SOURCE/flow_statistics.f90

@@ -4 +4 @@
 ! Current revisions:
 ! -----------------
-! 
+! Bugfix in calculation of vertical velocity skewness. The added absolute value
+! avoid negative values in the root. Negative values of w'w' can occur at the
+! top or bottom of the model domain due to degrading the order of advection 
+! scheme. Furthermore the calculation will be the same for all advection
+! schemes.
 !
 ! Former revisions:
@@ -437 +441 @@
                                   ( q(k,j,i)-hom(k,1,41,sr) )**2 * rmask(j,i,sr)
                 ENDIF
+                
+!
+!--             Higher moments
+!--             (Computation of the skewness of w further below)
+                sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i)**3 * rmask(j,i,sr)
 
                 sums_l_etot  = sums_l_etot + &
@@ -484 +493 @@
                    sums_l(k,32,tn) = sums_l(k,32,tn) + w2   * rmask(j,i,sr)
 !
-!--   Higher moments
-!--  (Computation of the skewness of w further below)
-                   sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i) * w2 * &
-                                                    rmask(j,i,sr)
-!
 !--             Perturbation energy
 
@@ -501 +505 @@
              ENDDO
           ENDDO
-       ELSE
-          !$OMP DO
-          DO  i = nxl, nxr
-             DO  j =  nys, nyn
-                DO  k = nzb_s_inner(j,i), nzt + 1
-                   w2   = w(k,j,i)**2
-!
-!--                Higher moments
-!--                (Computation of the skewness of w further below)
-                   sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i) * w2 * &
-                                                    rmask(j,i,sr)
-                ENDDO
-             ENDDO
-          ENDDO
-       ENDIF
-
+       ENDIF
 !
 !--    Horizontally averaged profiles of the vertical fluxes
@@ -1031 +1020 @@
        hom(:,1,37,sr) = sums(:,37)     ! w*e*
        hom(:,1,38,sr) = sums(:,38)     ! w*3
-       hom(:,1,39,sr) = sums(:,38) / ( sums(:,32) + 1E-20 )**1.5    ! Sw
+       hom(:,1,39,sr) = sums(:,38) / ( abs( sums(:,32) ) + 1E-20 )**1.5   ! Sw
        hom(:,1,40,sr) = sums(:,40)     ! p
        hom(:,1,45,sr) = sums(:,45)     ! w"vpt"