Changeset 1007 for palm/trunk/SOURCE/flow_statistics.f90

Timestamp:

Sep 19, 2012 2:30:36 PM (12 years ago)

Author:

franke

Message:

Bugfixes

---

Missing calculation of mean particle weighting factor for output added. (data_output_2d, data_output_3d, data_output_mask, sum_up_3d_data)
Calculation of mean particle radius for output now considers the weighting factor. (data_output_mask)
Calculation of sugrid-scale buoyancy flux for humidity and cloud droplets corrected. (flow_statistics)
Factor in calculation of enhancement factor for collision efficencies corrected. (lpm_collision_kernels)
Calculation of buoyancy production now considers the liquid water mixing ratio in case of cloud droplets. (production_e)

Changes

---

Calculation of buoyancy flux for humidity in case of WS-scheme is now using turbulent fluxes of WS-scheme. (flow_statistics)
Calculation of the collision kernels now in SI units. (lpm_collision_kernels)

File:

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

palm/trunk/SOURCE/flow_statistics.f90

@@ -4 +4 @@
 ! Current revisions:
 ! -----------------
+! Calculation of buoyancy flux for humidity in case of WS-scheme is now using
+! turbulent fluxes of WS-scheme
+! Bugfix: Calculation of subgridscale buoyancy flux for humidity and cloud
+! droplets at nzb and nzt added
 !
 ! Former revisions:
@@ -151 +155 @@
                         'timestep'
        CALL message( 'flow_statistics', 'PA0190', 1, 2, 0, 6, 0 )
-      
+
     ENDIF
 
@@ -176 +180 @@
        IF ( ws_scheme_mom .AND. sr == 0 )  THEN
 
-!       
+!
 !--       According to the Neumann bc for the horizontal velocity components,
 !--       the corresponding fluxes has to satisfiy the same bc.
           IF ( ocean )  THEN
              sums_us2_ws_l(nzt+1,:) = sums_us2_ws_l(nzt,:)
-             sums_vs2_ws_l(nzt+1,:) = sums_vs2_ws_l(nzt,:)    
+             sums_vs2_ws_l(nzt+1,:) = sums_vs2_ws_l(nzt,:)
           ENDIF
 
           DO  i = 0, threads_per_task-1
-!          
+!
 !--          Swap the turbulent quantities evaluated in advec_ws.
              sums_l(:,13,i) = sums_wsus_ws_l(:,i)       ! w*u* 
@@ -453 +457 @@
                                   ( q(k,j,i)-hom(k,1,41,sr) )**2 * rmask(j,i,sr)
                 ENDIF
-                
+
 !
 !--             Higher moments
@@ -578 +582 @@
                                                * ( q(k+1,j,i) - q(k,j,i) )     &
                                                * ddzu(k+1) * rmask(j,i,sr)
+
                    IF ( cloud_physics ) THEN
                       sums_l(k,51,tn) = sums_l(k,51,tn)                        &
@@ -618 +623 @@
                    sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + &
                                        qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
+                   sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (           &
+                                       ( 1.0 + 0.61 * q(nzb,j,i) ) *   &
+                                       shf(j,i) + 0.61 * pt(nzb,j,i) * &
+                                                  qsws(j,i) )
+                   IF ( cloud_droplets )  THEN
+                      sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (        &
+                                         ( 1.0 + 0.61 * q(nzb,j,i) -   &
+                                           ql(nzb,j,i) ) * shf(j,i) +  &
+                                           0.61 * pt(nzb,j,i) * qsws(j,i) )
+                   ENDIF
                    IF ( cloud_physics )  THEN
-                      sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (           &
-                                          ( 1.0 + 0.61 * q(nzb,j,i) ) *   &
-                                          shf(j,i) + 0.61 * pt(nzb,j,i) * &
-                                                     qsws(j,i)            &
-                                                              )
 !
 !--                   Formula does not work if ql(nzb) /= 0.0
@@ -657 +667 @@
                    sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
                                        qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
+                   sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (             &
+                                       ( 1.0 + 0.61 * q(nzt,j,i) ) *     &
+                                       tswst(j,i) + 0.61 * pt(nzt,j,i) * &
+                                                  qswst(j,i) )
+                   IF ( cloud_droplets )  THEN
+                      sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (           &
+                                          ( 1.0 + 0.61 * q(nzt,j,i) -     &
+                                            ql(nzt,j,i) ) * tswst(j,i) +  &
+                                            0.61 * pt(nzt,j,i) * qswst(j,i) )
+                   ENDIF
                    IF ( cloud_physics )  THEN
-                      sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (           &
-                                          ( 1.0 + 0.61 * q(nzt,j,i) ) *   &
-                                          tswst(j,i) + 0.61 * pt(nzt,j,i) * &
-                                                     qsws(j,i)            &
-                                                              )
 !
 !--                   Formula does not work if ql(nzb) /= 0.0
@@ -714 +729 @@
 !--             content
                 IF ( humidity )  THEN
-                   pts = 0.5 * ( vpt(k,j,i)   - hom(k,1,44,sr) + &
-                                 vpt(k+1,j,i) - hom(k+1,1,44,sr) )
-                   sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * &
-                                                       rmask(j,i,sr)
-
-                   IF ( cloud_physics  .OR.  cloud_droplets )  THEN
-                      pts = 0.5 *                                            &
-                             ( ( q(k,j,i)   - ql(k,j,i)   ) - hom(k,1,42,sr) &
-                             + ( q(k+1,j,i) - ql(k+1,j,i) ) - hom(k+1,1,42,sr) )
-                      sums_l(k,52,tn) = sums_l(k,52,tn) + pts * w(k,j,i) * &
+                   IF ( cloud_physics .OR. cloud_droplets )  THEN
+                      pts = 0.5 * ( vpt(k,j,i)   - hom(k,1,44,sr) + &
+                                    vpt(k+1,j,i) - hom(k+1,1,44,sr) )
+                      sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * &
                                                           rmask(j,i,sr)
                       sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * &
                                                           rmask(j,i,sr)
-                   ENDIF
-                ENDIF
-
+                   ELSE
+                      IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
+                         pts = 0.5 * ( vpt(k,j,i)   - hom(k,1,44,sr) + &
+                                       vpt(k+1,j,i) - hom(k+1,1,44,sr) )
+                         sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * &
+                                                             rmask(j,i,sr)
+                      ELSE IF ( ws_scheme_sca .AND. sr == 0 )  THEN
+                         sums_l(k,46,tn) = ( 1.0 + 0.61 * hom(k,1,41,sr) ) * &
+                                             sums_l(k,17,tn) +      &
+                                          0.61 * hom(k,1,4,sr) * sums_l(k,49,tn)
+                      END IF
+                   END IF
+                ENDIF
 !
 !--             Passive scalar flux
@@ -762 +781 @@
                   vst = 0.5 * ( v(k,j,i)   - hom(k,1,2,sr) + &
                               v(k+1,j,i) - hom(k+1,1,2,sr) )
-!                              
+!
 !--               Momentum flux w*u*
                   sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5 *                   &
@@ -1002 +1021 @@
           ENDDO
        ENDIF
+
 !
 !--    Collect horizontal average in hom.