Index: palm/trunk/SOURCE/advec_ws.f90
===================================================================
--- palm/trunk/SOURCE/advec_ws.f90	(revision 4315)
+++ palm/trunk/SOURCE/advec_ws.f90	(revision 4317)
@@ -25,4 +25,9 @@
 ! -----------------
 ! $Id$
+! Comments revised/added, formatting improved, fluxes for u,v, and scalars are 
+! explicitly set to zero at nzt+1, fluxes of w-component are now calculated only
+! until nzt-1 (Prognostic equation for w-velocity component ends at nzt-1)
+! 
+! 4204 2019-08-30 12:30:17Z knoop
 ! Bugfix: Changed sk_num initialization default to avoid implicit SAVE-Attribut
 ! 
@@ -450,5 +455,5 @@
                 ENDIF                                                          
 !                                                                              
-!--             u component - z-direction                                      
+!--             u component - z-direction. Fluxes are calculated on w-grid level                                      
 !--             WS1 (6), WS3 (7), WS5 (8)                                      
                 IF ( k == nzb+1 )  THEN                                        
@@ -479,9 +484,9 @@
                 ELSEIF ( ( .NOT. BTEST(wall_flags_0(k_mm,j,i),1)    .OR.       &
                            .NOT. BTEST(wall_flags_0(k_ppp,j,i),1) ) .AND.      & 
-                               BTEST(wall_flags_0(k-1,j,i),1)  .AND.           &
-                               BTEST(wall_flags_0(k,j,i),1)    .AND.           &
-                               BTEST(wall_flags_0(k+1,j,i),1)  .AND.           &
-                               BTEST(wall_flags_0(k_pp,j,i),1) .OR.            &
-                               k == nzt - 1 )                                  &
+                                 BTEST(wall_flags_0(k-1,j,i),1)     .AND.      &
+                                 BTEST(wall_flags_0(k,j,i),1)       .AND.      &
+                                 BTEST(wall_flags_0(k+1,j,i),1)     .AND.      &
+                                 BTEST(wall_flags_0(k_pp,j,i),1)    .OR.       &
+                         k == nzt - 1 )                                        &
                 THEN                                                           
                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 7 )       
@@ -579,5 +584,5 @@
                 ENDIF
 !
-!--             v component - z-direction
+!--             v component - z-direction. Fluxes are calculated on w-grid level
 !--             WS1 (15), WS3 (16), WS5 (17)
                 IF ( k == nzb+1 )  THEN
@@ -608,9 +613,9 @@
                 ELSEIF ( ( .NOT. BTEST(wall_flags_0(k_mm,j,i),2)    .OR.       &
                            .NOT. BTEST(wall_flags_0(k_ppp,j,i),2) ) .AND.      & 
-                               BTEST(wall_flags_0(k-1,j,i),2)  .AND.           &
-                               BTEST(wall_flags_0(k,j,i),2)    .AND.           &
-                               BTEST(wall_flags_0(k+1,j,i),2)  .AND.           &
-                               BTEST(wall_flags_0(k_pp,j,i),2)  .OR.           &
-                               k == nzt - 1 )                                  &
+                                 BTEST(wall_flags_0(k-1,j,i),2)     .AND.      &
+                                 BTEST(wall_flags_0(k,j,i),2)       .AND.      &
+                                 BTEST(wall_flags_0(k+1,j,i),2)     .AND.      &
+                                 BTEST(wall_flags_0(k_pp,j,i),2)    .OR.       &
+                         k == nzt - 1 )                                        &
                 THEN                                                           
                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 16 )      
@@ -706,7 +711,6 @@
                 ENDIF
 !
-!--             w component - z-direction
-!--             WS1 (24), WS3 (25), WS5 (26)
-                flag_set = .FALSE.
+!--             w component - z-direction. Fluxes are calculated on scalar grid
+!--             level. WS1 (24), WS3 (25), WS5 (26)
                 IF ( k == nzb+1 )  THEN
                    k_mm = nzb
@@ -725,15 +729,16 @@
                 ENDIF  
                 
-                IF ( ( .NOT. BTEST(wall_flags_0(k-1,j,i),3)  .AND.             &
-                       .NOT. BTEST(wall_flags_0(k,j,i),3)    .AND.             &
-                             BTEST(wall_flags_0(k+1,j,i),3) )  .OR.            &
-                     ( .NOT. BTEST(wall_flags_0(k-1,j,i),3)  .AND.             &
-                             BTEST(wall_flags_0(k,j,i),3) )  .OR.              &
-                     ( .NOT. BTEST(wall_flags_0(k+1,j,i),3)  .AND.             &
-                             BTEST(wall_flags_0(k,j,i),3) )  .OR.              &        
+                flag_set = .FALSE.
+                IF ( ( .NOT. BTEST(wall_flags_0(k-1,j,i),3)   .AND.            &
+                       .NOT. BTEST(wall_flags_0(k,j,i),3)     .AND.            &
+                             BTEST(wall_flags_0(k+1,j,i),3) ) .OR.             &
+                     ( .NOT. BTEST(wall_flags_0(k-1,j,i),3)   .AND.            &
+                             BTEST(wall_flags_0(k,j,i),3) )   .OR.             &
+                     ( .NOT. BTEST(wall_flags_0(k+1,j,i),3)   .AND.            &
+                             BTEST(wall_flags_0(k,j,i),3) )   .OR.             &        
                      k == nzt )                                                &
                 THEN
 !
-!--                Please note, at k == nzb_w_inner(j,i) a flag is explictely
+!--                Please note, at k == nzb_w_inner(j,i) a flag is explicitly
 !--                set, although this is not a prognostic level. However, 
 !--                contrary to the advection of u,v and s this is necessary
@@ -742,9 +747,9 @@
                    advc_flags_m(k,j,i) = IBSET( advc_flags_m(k,j,i), 24 )
                    flag_set = .TRUE.
-                ELSEIF ( ( .NOT. BTEST(wall_flags_0(k_mm,j,i),3)     .OR.      &
+                ELSEIF ( ( .NOT. BTEST(wall_flags_0(k_mm,j,i),3)    .OR.       &
                            .NOT. BTEST(wall_flags_0(k_ppp,j,i),3) ) .AND.      &
-                                 BTEST(wall_flags_0(k-1,j,i),3)  .AND.         &
-                                 BTEST(wall_flags_0(k,j,i),3)    .AND.         &
-                                 BTEST(wall_flags_0(k+1,j,i),3)  .OR.          &
+                                 BTEST(wall_flags_0(k-1,j,i),3)     .AND.      &
+                                 BTEST(wall_flags_0(k,j,i),3)       .AND.      &
+                                 BTEST(wall_flags_0(k+1,j,i),3)     .OR.       &
                          k == nzt - 1 )                                        &
                 THEN                                                           
@@ -1027,5 +1032,5 @@
                 
 !                                                                              
-!--             scalar - z-direction                                           
+!--             scalar - z-direction. Fluxes are calculated on w-grid level                                           
 !--             WS1 (6), WS3 (7), WS5 (8)                                      
                 IF ( k == nzb+1 )  THEN                                        
@@ -1056,19 +1061,19 @@
                 ELSEIF ( ( .NOT. BTEST(wall_flags_0(k_mm,j,i),0)    .OR.       &
                            .NOT. BTEST(wall_flags_0(k_ppp,j,i),0) ) .AND.      & 
-                               BTEST(wall_flags_0(k-1,j,i),0)  .AND.           &
-                               BTEST(wall_flags_0(k,j,i),0)    .AND.           &
-                               BTEST(wall_flags_0(k+1,j,i),0)  .AND.           &
-                               BTEST(wall_flags_0(k_pp,j,i),0) .OR.            &    
+                                 BTEST(wall_flags_0(k-1,j,i),0)     .AND.      &
+                                 BTEST(wall_flags_0(k,j,i),0)       .AND.      &
+                                 BTEST(wall_flags_0(k+1,j,i),0)     .AND.      &
+                                 BTEST(wall_flags_0(k_pp,j,i),0)    .OR.       &    
                          k == nzt - 1 )                                        &
                 THEN                                                           
                    advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 7 )             
                    flag_set = .TRUE.                                           
-                ELSEIF ( BTEST(wall_flags_0(k_mm,j,i),0)                       &
-                  .AND.  BTEST(wall_flags_0(k-1,j,i),0)                        &
-                  .AND.  BTEST(wall_flags_0(k,j,i),0)                          &
-                  .AND.  BTEST(wall_flags_0(k+1,j,i),0)                        &
-                  .AND.  BTEST(wall_flags_0(k_pp,j,i),0)                       &
-                  .AND.  BTEST(wall_flags_0(k_ppp,j,i),0)                      &
-                  .AND. .NOT. flag_set )                                       &
+                ELSEIF ( BTEST(wall_flags_0(k_mm,j,i),0)  .AND.                &
+                         BTEST(wall_flags_0(k-1,j,i),0)   .AND.                &
+                         BTEST(wall_flags_0(k,j,i),0)     .AND.                &
+                         BTEST(wall_flags_0(k+1,j,i),0)   .AND.                &
+                         BTEST(wall_flags_0(k_pp,j,i),0)  .AND.                &
+                         BTEST(wall_flags_0(k_ppp,j,i),0) .AND.                &
+                        .NOT. flag_set )                                       &
                 THEN
                    advc_flag(k,j,i) = IBSET( advc_flag(k,j,i), 8 )
@@ -1480,4 +1485,5 @@
        flux_t(nzb) = 0.0_wp
        diss_t(nzb) = 0.0_wp
+       
        DO  k = nzb+1, nzb+2
           ibit8 = REAL( IBITS(advc_flag(k,j,i),8,1), KIND = wp )
@@ -1598,5 +1604,10 @@
                                                          )
        ENDDO
-
+       
+!
+!--    Set resolved/turbulent flux at model top to zero (w-level)
+       flux_t(nzt+1) = 0.0_wp
+       diss_t(nzt+1) = 0.0_wp
+       
        IF ( limiter )  THEN
 !
@@ -2191,5 +2202,5 @@
 !--    indirect indexing. This allows better vectorization for the main loop.
 !--    First, compute the flux at model surface, which need has to be 
-!--    calculated explicetely for the tendency at
+!--    calculated explicitly for the tendency at
 !--    the first w-level. For topography wall this is done implicitely by
 !--    advc_flags_m.
@@ -2197,4 +2208,5 @@
        diss_t(nzb) = 0.0_wp
        w_comp(nzb) = 0.0_wp
+       
        DO  k = nzb+1, nzb+2
 !
@@ -2322,4 +2334,10 @@
                                                            )
        ENDDO
+       
+!
+!--    Set resolved/turbulent flux at model top to zero (w-level)
+       flux_t(nzt+1) = 0.0_wp
+       diss_t(nzt+1) = 0.0_wp
+       w_comp(nzt+1) = 0.0_wp
        
        DO  k = nzb+1, nzb_max_l
@@ -2748,5 +2766,5 @@
 !--    indirect indexing. This allows better vectorization for the main loop.
 !--    First, compute the flux at model surface, which need has to be 
-!--    calculated explicetely for the tendency at
+!--    calculated explicitly for the tendency at
 !--    the first w-level. For topography wall this is done implicitely by
 !--    advc_flags_m.
@@ -2880,4 +2898,10 @@
                                                            )
        ENDDO
+       
+!
+!--    Set resolved/turbulent flux at model top to zero (w-level)
+       flux_t(nzt+1) = 0.0_wp
+       diss_t(nzt+1) = 0.0_wp
+       w_comp(nzt+1) = 0.0_wp
        
        DO  k = nzb+1, nzb_max_l
@@ -3125,5 +3149,5 @@
           ENDDO
 
-          DO  k = nzb_max_l+1, nzt
+          DO  k = nzb_max_l+1, nzt-1
 
              v_comp(k)      = v(k+1,j,i) + v(k,j,i) - gv
@@ -3183,5 +3207,5 @@
           ENDDO
 
-          DO  k = nzb_max_l+1, nzt
+          DO  k = nzb_max_l+1, nzt-1
 
              u_comp(k)        = u(k+1,j,i) + u(k,j,i) - gu
@@ -3274,5 +3298,5 @@
        ENDDO
 
-       DO  k = nzb_max_l+1, nzt
+       DO  k = nzb_max_l+1, nzt-1
 
           u_comp(k) = u(k+1,j,i+1) + u(k,j,i+1) - gu
@@ -3305,5 +3329,5 @@
 !--    indirect indexing. This allows better vectorization for the main loop.
 !--    First, compute the flux at model surface, which need has to be 
-!--    calculated explicetely for the tendency at
+!--    calculated explicitly for the tendency at
 !--    the first w-level. For topography wall this is done implicitely by
 !--    advc_flags_m.
@@ -3395,5 +3419,5 @@
        ENDDO
        
-       DO  k = nzt-1, nzt
+       DO  k = nzt-1, nzt-1
 !
 !--       k index has to be modified near bottom and top, else array
@@ -3513,5 +3537,5 @@
        ENDDO
        
-       DO  k = nzb_max_l+1, nzt
+       DO  k = nzb_max_l+1, nzt-1
 
           flux_d    = flux_t(k-1)
@@ -5870,5 +5894,5 @@
           ENDDO
 
-          DO  k = nzb_max_l+1, nzt
+          DO  k = nzb_max_l+1, nzt-1
 
              u_comp                   = u(k+1,j,i) + u(k,j,i) - gu
@@ -5947,5 +5971,5 @@
           ENDDO
 
-          DO  k = nzb_max_l+1, nzt
+          DO  k = nzb_max_l+1, nzt-1
 
              v_comp                 = v(k+1,j,i) + v(k,j,i) - gv
@@ -5965,5 +5989,5 @@
 
 !
-!--          The lower flux has to be calculated explicetely for the tendency
+!--          The lower flux has to be calculated explicitly for the tendency
 !--          at the first w-level. For topography wall this is done implicitely
 !--          by advc_flags_m.
@@ -6233,5 +6257,5 @@
              ENDDO
 
-             DO  k = nzb_max_l+1, nzt
+             DO  k = nzb_max_l+1, nzt-1
 
                 u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu
