Changeset 2026 for palm

Timestamp:

Oct 18, 2016 10:27:02 AM (8 years ago)

Author:

suehring

Message:

Bugfix, enable output of s*2; calculation of domain_averaged perturbation energy; some formatting adjustments

Location:

palm/trunk/SOURCE

Files:

• check_parameters.f90 (modified) (2 diffs)
• data_output_profiles.f90 (modified) (2 diffs)
• flow_statistics.f90 (modified) (9 diffs)

palm/trunk/SOURCE/check_parameters.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-! 
+! Bugfix, enable output of s*2
 ! 
 ! Former revisions:
@@ -2799 +2799 @@
              ENDIF
 
+          CASE ( 's*2' )
+             IF (  .NOT.  passive_scalar )  THEN
+                message_string = 'data_output_pr = ' //                        &
+                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
+                                 'lemented for passive_scalar = .FALSE.'
+                CALL message( 'check_parameters', 'PA0185', 1, 2, 0, 6, 0 )
+             ELSE
+                dopr_index(i) = 118
+                dopr_unit(i)  = 'kg2/kg2'
+                hom(:,2,118,:) = SPREAD( zu, 2, statistic_regions+1 )
+             ENDIF
+
  
 

palm/trunk/SOURCE/data_output_profiles.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-! 
+! Formatting adjustment
 ! 
 ! Former revisions:
@@ -164 +164 @@
 !--          runs, has a negative value
              DO  i = 1, dopr_n
-             IF ( dopr_initial_index(i) /= 0 )  THEN
-                nc_stat = NF90_PUT_VAR( id_set_pr, id_var_time_pr,  &
-                                        (/ -coupling_start_time /), &
+                IF ( dopr_initial_index(i) /= 0 )  THEN
+                   nc_stat = NF90_PUT_VAR( id_set_pr, id_var_time_pr,          &
+                                        (/ -coupling_start_time /),            &
                                         start = (/ 1 /), count = (/ 1 /) )
                    CALL netcdf_handle_error( 'data_output_profiles', 329 )

palm/trunk/SOURCE/flow_statistics.f90

@@ -20 +20 @@
 ! Current revisions:
 ! -----------------
-! 
+! Bugfix, enable output of s*2.
+! Change, calculation of domain-averaged perturbation energy.
+! Some formatting adjustments.
 ! 
 ! Former revisions:
@@ -358 +360 @@
                               ( sums_us2_ws_l(:,i) + sums_vs2_ws_l(:,i) +      &
                                 sums_ws2_ws_l(:,i) )    ! e*
-             DO  k = nzb, nzt
-                sums_l(nzb+5,pr_palm,i) = sums_l(nzb+5,pr_palm,i) + 0.5_wp * ( &
-                                                      sums_us2_ws_l(k,i) +     &
-                                                      sums_vs2_ws_l(k,i) +     &
-                                                      sums_ws2_ws_l(k,i) )
-             ENDDO
           ENDDO
 
@@ -657 +653 @@
           DO  i = nxl, nxr
              DO  j =  nys, nyn
-                sums_l_eper = 0.0_wp
                 DO  k = nzb_s_inner(j,i), nzt+1
                    u2   = u(k,j,i)**2
@@ -669 +664 @@
                    sums_l(k,32,tn) = sums_l(k,32,tn) + w2   * rmask(j,i,sr)
 !
-!--             Perturbation energy
+!--                Perturbation energy
 
                    sums_l(k,34,tn) = sums_l(k,34,tn) + 0.5_wp *                &
                                   ( ust2 + vst2 + w2 ) * rmask(j,i,sr)
-                   sums_l_eper     = sums_l_eper +                             &
-                                     0.5_wp * ( ust2+vst2+w2 ) * rmask(j,i,sr)
-
-                ENDDO
+                ENDDO
+             ENDDO
+          ENDDO
+       ENDIF
+!
+!--    Computaion of domain-averaged perturbation energy. Please note, 
+!--    to prevent that perturbation energy is larger (even if only slightly) 
+!--    than the total kinetic energy, calculation is based on deviations from
+!--    the horizontal mean, instead of spatial descretization of the advection
+!--    term. 
+       !$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
+                ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2
+                vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2
+                w2   = w(k,j,i)**2
+
                 sums_l(nzb+5,pr_palm,tn) = sums_l(nzb+5,pr_palm,tn)            &
-                     + sums_l_eper
-             ENDDO
-          ENDDO
-       ENDIF
+                                 + 0.5_wp * ( ust2 + vst2 + w2 ) * rmask(j,i,sr)
+             ENDDO
+          ENDDO
+       ENDDO
 
 !
@@ -758 +768 @@
                    sums_l(k,119,tn) = sums_l(k,119,tn)                         &
                                          - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
-                                               * ( s(k+1,j,i) - s(k,j,i) )     &
-                                               * ddzu(k+1) * rmask(j,i,sr)
+                                                  * ( s(k+1,j,i) - s(k,j,i) )  &
+                                                  * ddzu(k+1) * rmask(j,i,sr)
                 ENDIF
 
@@ -803 +813 @@
                 IF ( passive_scalar )  THEN
                    sums_l(nzb,119,tn) = sums_l(nzb,119,tn) +                     &
-                                       ssws(j,i) * rmask(j,i,sr) ! w"s"
+                                        ssws(j,i) * rmask(j,i,sr) ! w"s"
                 ENDIF
              ENDIF
@@ -880 +890 @@
                 IF ( passive_scalar )  THEN
                    sums_l(nzt,119,tn) = sums_l(nzt,119,tn) + &
-                                       sswst(j,i) * rmask(j,i,sr) ! w"s"
+                                        sswst(j,i) * rmask(j,i,sr) ! w"s"
                 ENDIF
              ENDIF
@@ -1032 +1042 @@
                      pts = 0.5_wp * ( s(k,j,i)   - hom(k,1,117,sr) +            &
                                       s(k+1,j,i) - hom(k+1,1,117,sr) )
-                     sums_l(k,116,tn) = sums_l(k,116,tn) + pts * w(k,j,i) *      &
-                                       rmask(j,i,sr)
+                     sums_l(k,116,tn) = sums_l(k,116,tn) + pts * w(k,j,i) *     &
+                                        rmask(j,i,sr)
                   ENDIF
                ENDDO
@@ -1470 +1480 @@
           hom(:,1,116,sr) = sums(:,116)     ! w*s*
           hom(:,1,120,sr) = sums(:,119) + sums(:,116)    ! ws 
+          hom(:,1,118,sr) = sums(:,118)     ! s*2
        ENDIF